From 2e063bcd3eec4b96b735be25dcb14bfb3dbd47e7 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 14 Aug 2024 19:06:21 -0400 Subject: [PATCH 01/25] Copy interpolatae_data.F90 utility routine from CAM; initial tropopause_read_file.F90 utility implementation --- src/physics/utils/tropopause_read_file.F90 | 160 +++ src/utils/interpolate_data.F90 | 1230 ++++++++++++++++++++ 2 files changed, 1390 insertions(+) create mode 100644 src/physics/utils/tropopause_read_file.F90 create mode 100644 src/utils/interpolate_data.F90 diff --git a/src/physics/utils/tropopause_read_file.F90 b/src/physics/utils/tropopause_read_file.F90 new file mode 100644 index 00000000..e240fb64 --- /dev/null +++ b/src/physics/utils/tropopause_read_file.F90 @@ -0,0 +1,160 @@ +module tropopause_read_file + !------------------------------------------------------------------- + ! Support module for CAM-SIMA tropopause_find to read in + ! climatological tropopause data. + ! + ! Remarks: This module is not CCPP-ized, but has been cleaned up + ! for use within CAM-SIMA, particularly removal of chunk support. + !------------------------------------------------------------------- + + implicit none + private + + public :: tropopause_read_file + + ! These variables are used to store the climatology data. + real(kind_phys) :: days(12) ! days in the climatology + real(kind_phys), pointer :: tropp_p_loc(:,:,:) ! climatological tropopause pressures (pcols,1,ntimes) + +contains + subroutine tropopause_read_file(tropopause_climo_file) + !------------------------------------------------------------------ + ! ... initialize upper boundary values + !------------------------------------------------------------------ + use interpolate_data, only : lininterp_init, lininterp, interp_type, lininterp_finish + use physics_grid, only : get_rlat_all_p, get_rlon_all_p + use physics_grid, only : pcols => columns_on_task + use physconst, only : pi + use time_manager, only : get_calday + use ioFileMod, only : getfil + use cam_pio_utils, only: cam_pio_openfile + use pio, only : file_desc_t, var_desc_t, pio_inq_dimid, pio_inq_dimlen, & + pio_inq_varid, pio_get_var, pio_closefile, pio_nowrite + + character(len=256), intent(in) :: tropopause_climo_file ! absolute path of climatology file + + !------------------------------------------------------------------ + ! ... local variables + !------------------------------------------------------------------ + integer :: i, j, n + integer :: ierr + type(file_desc_t) :: pio_id + integer :: dimid + type(var_desc_t) :: vid + integer :: nlon, nlat, ntimes + integer :: start(3) + integer :: count(3) + integer, parameter :: dates(12) = (/ 116, 214, 316, 415, 516, 615, & + 716, 816, 915, 1016, 1115, 1216 /) + type(interp_type) :: lon_wgts, lat_wgts + real(kind_phys), allocatable :: tropp_p_in(:,:,:) + real(kind_phys), allocatable :: lat(:) + real(kind_phys), allocatable :: lon(:) + real(kind_phys) :: to_lats(pcols), to_lons(pcols) + real(kind_phys), parameter :: d2r=pi/180._kind_phys, zero=0._kind_phys, twopi=pi*2._kind_phys + character(len=256) :: locfn + + !----------------------------------------------------------------------- + ! ... open netcdf file + !----------------------------------------------------------------------- + call getfil (tropopause_climo_file, locfn, 0) + call cam_pio_openfile(pio_id, trim(locfn), PIO_NOWRITE) + + !----------------------------------------------------------------------- + ! ... get time dimension + !----------------------------------------------------------------------- + ierr = pio_inq_dimid( pio_id, 'time', dimid ) + ierr = pio_inq_dimlen( pio_id, dimid, ntimes ) + if( ntimes /= 12 )then + write(iulog,*) 'tropopause_find_init: number of months = ',ntimes,'; expecting 12' + call endrun + end if + !----------------------------------------------------------------------- + ! ... get latitudes + !----------------------------------------------------------------------- + ierr = pio_inq_dimid( pio_id, 'lat', dimid ) + ierr = pio_inq_dimlen( pio_id, dimid, nlat ) + allocate( lat(nlat), stat=ierr ) + if( ierr /= 0 ) then + write(iulog,*) 'tropopause_find_init: lat allocation error = ',ierr + call endrun + end if + ierr = pio_inq_varid( pio_id, 'lat', vid ) + ierr = pio_get_var( pio_id, vid, lat ) + lat(:nlat) = lat(:nlat) * d2r + !----------------------------------------------------------------------- + ! ... get longitudes + !----------------------------------------------------------------------- + ierr = pio_inq_dimid( pio_id, 'lon', dimid ) + ierr = pio_inq_dimlen( pio_id, dimid, nlon ) + allocate( lon(nlon), stat=ierr ) + if( ierr /= 0 ) then + write(iulog,*) 'tropopause_find_init: lon allocation error = ',ierr + call endrun + end if + ierr = pio_inq_varid( pio_id, 'lon', vid ) + ierr = pio_get_var( pio_id, vid, lon ) + lon(:nlon) = lon(:nlon) * d2r + + !------------------------------------------------------------------ + ! ... allocate arrays + !------------------------------------------------------------------ + allocate( tropp_p_in(nlon,nlat,ntimes), stat=ierr ) + if( ierr /= 0 ) then + write(iulog,*) 'tropopause_find_init: tropp_p_in allocation error = ',ierr + call endrun + end if + !------------------------------------------------------------------ + ! ... read in the tropopause pressure + !------------------------------------------------------------------ + ierr = pio_inq_varid( pio_id, 'trop_p', vid ) + start = (/ 1, 1, 1 /) + count = (/ nlon, nlat, ntimes /) + ierr = pio_get_var( pio_id, vid, start, count, tropp_p_in ) + + !------------------------------------------------------------------ + ! ... close the netcdf file + !------------------------------------------------------------------ + call pio_closefile( pio_id ) + + !-------------------------------------------------------------------- + ! ... regrid + !-------------------------------------------------------------------- + + ! tropp_p_loc is allocated with dimensions (pcols, begchunk:endchunk, ntimes) in CAM. + ! to maintain backwards compatibility with existing CAM, the second dimension here + ! must exist but is set to 1 in CAM-SIMA. (hplin, 8/14/24) + allocate( tropp_p_loc(pcols,1,ntimes), stat=ierr ) + + if( ierr /= 0 ) then + write(iulog,*) 'tropopause_find_init: tropp_p_loc allocation error = ',ierr + call endrun + end if + + call get_rlat_all_p(pcols, to_lats) + call get_rlon_all_p(pcols, to_lons) + call lininterp_init(lon, nlon, to_lons, pcols, 2, lon_wgts, zero, twopi) + call lininterp_init(lat, nlat, to_lats, pcols, 1, lat_wgts) + do n=1,ntimes + call lininterp(tropp_p_in(:,:,n), nlon, nlat, tropp_p_loc(1:pcols,c,n), pcols, lon_wgts, lat_wgts) + end do + call lininterp_finish(lon_wgts) + call lininterp_finish(lat_wgts) + deallocate(lon) + deallocate(lat) + deallocate(tropp_p_in) + + !-------------------------------------------------------- + ! ... initialize the monthly day of year times + !-------------------------------------------------------- + + do n = 1,12 + days(n) = get_calday( dates(n), 0 ) + end do + if (masterproc) then + write(iulog,*) 'tropopause_find_init : days' + write(iulog,'(1p,5g15.8)') days(:) + endif + + end subroutine tropopause_read_file +end module tropopause_read_file \ No newline at end of file diff --git a/src/utils/interpolate_data.F90 b/src/utils/interpolate_data.F90 new file mode 100644 index 00000000..142e3600 --- /dev/null +++ b/src/utils/interpolate_data.F90 @@ -0,0 +1,1230 @@ +module interpolate_data +! Description: +! Routines for interpolation of data in latitude, longitude and time. +! +! Author: Gathered from a number of places and put into the current format by Jim Edwards +! +! Modules Used: +! + use shr_kind_mod, only: r8 => shr_kind_r8 + use cam_abortutils, only: endrun + use cam_logfile, only: iulog + implicit none + private +! +! Public Methods: +! + + public :: interp_type, lininterp, vertinterp, bilin, get_timeinterp_factors + public :: lininterp_init, lininterp_finish + type interp_type + real(r8), pointer :: wgts(:) + real(r8), pointer :: wgtn(:) + integer, pointer :: jjm(:) + integer, pointer :: jjp(:) + end type interp_type + interface lininterp + module procedure lininterp_original + module procedure lininterp_full1d + module procedure lininterp1d + module procedure lininterp2d2d + module procedure lininterp2d1d + module procedure lininterp3d2d + end interface + +integer, parameter, public :: extrap_method_zero = 0 +integer, parameter, public :: extrap_method_bndry = 1 +integer, parameter, public :: extrap_method_cycle = 2 + +contains + subroutine lininterp_full1d (arrin, yin, nin, arrout, yout, nout) + integer, intent(in) :: nin, nout + real(r8), intent(in) :: arrin(nin), yin(nin), yout(nout) + real(r8), intent(out) :: arrout(nout) + type (interp_type) :: interp_wgts + + call lininterp_init(yin, nin, yout, nout, extrap_method_bndry, interp_wgts) + call lininterp1d(arrin, nin, arrout, nout, interp_wgts) + call lininterp_finish(interp_wgts) + + end subroutine lininterp_full1d + + subroutine lininterp_init(yin, nin, yout, nout, extrap_method, interp_wgts, & + cyclicmin, cyclicmax) +! +! Description: +! Initialize a variable of type(interp_type) with weights for linear interpolation. +! this variable can then be used in calls to lininterp1d and lininterp2d. +! yin is a 1d array of length nin of locations to interpolate from - this array must +! be monotonic but can be increasing or decreasing +! yout is a 1d array of length nout of locations to interpolate to, this array need +! not be ordered +! extrap_method determines how to handle yout points beyond the bounds of yin +! if 0 set values outside output grid to 0 +! if 1 set to boundary value +! if 2 set to cyclic boundaries +! optional values cyclicmin and cyclicmax can be used to set the bounds of the +! cyclic mapping - these default to 0 and 360. +! + + integer, intent(in) :: nin + integer, intent(in) :: nout + real(r8), intent(in) :: yin(:) ! input mesh + real(r8), intent(in) :: yout(:) ! output mesh + integer, intent(in) :: extrap_method ! if 0 set values outside output grid to 0 + ! if 1 set to boundary value + ! if 2 set to cyclic boundaries + real(r8), intent(in), optional :: cyclicmin, cyclicmax + + type (interp_type), intent(out) :: interp_wgts + real(r8) :: cmin, cmax + real(r8) :: extrap + real(r8) :: dyinwrap + real(r8) :: ratio + real(r8) :: avgdyin + integer :: i, j, icount + integer :: jj + real(r8), pointer :: wgts(:) + real(r8), pointer :: wgtn(:) + integer, pointer :: jjm(:) + integer, pointer :: jjp(:) + logical :: increasing + ! + ! Check validity of input coordinate arrays: must be monotonically increasing, + ! and have a total of at least 2 elements + ! + if (nin.lt.2) then + call endrun('LININTERP: Must have at least 2 input points for interpolation') + end if + if(present(cyclicmin)) then + cmin=cyclicmin + else + cmin=0_r8 + end if + if(present(cyclicmax)) then + cmax=cyclicmax + else + cmax=360_r8 + end if + if(cmax<=cmin) then + call endrun('LININTERP: cyclic min value must be < max value') + end if + increasing=.true. + icount = 0 + do j=1,nin-1 + if (yin(j).gt.yin(j+1)) icount = icount + 1 + end do + if(icount.eq.nin-1) then + increasing = .false. + icount=0 + endif + if (icount.gt.0) then + call endrun('LININTERP: Non-monotonic input coordinate array found') + end if + allocate(interp_wgts%jjm(nout), & + interp_wgts%jjp(nout), & + interp_wgts%wgts(nout), & + interp_wgts%wgtn(nout)) + + jjm => interp_wgts%jjm + jjp => interp_wgts%jjp + wgts => interp_wgts%wgts + wgtn => interp_wgts%wgtn + + ! + ! Initialize index arrays for later checking + ! + jjm = 0 + jjp = 0 + + extrap = 0._r8 + select case (extrap_method) + case (extrap_method_zero) + ! + ! For values which extend beyond boundaries, set weights + ! such that values will be 0. + ! + do j=1,nout + if(increasing) then + if (yout(j).lt.yin(1)) then + jjm(j) = 1 + jjp(j) = 1 + wgts(j) = 0._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + else if (yout(j).gt.yin(nin)) then + jjm(j) = nin + jjp(j) = nin + wgts(j) = 0._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + end if + else + if (yout(j).gt.yin(1)) then + jjm(j) = 1 + jjp(j) = 1 + wgts(j) = 0._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + else if (yout(j).lt.yin(nin)) then + jjm(j) = nin + jjp(j) = nin + wgts(j) = 0._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + end if + end if + end do + case (extrap_method_bndry) + ! + ! For values which extend beyond boundaries, set weights + ! such that values will just be copied. + ! + do j=1,nout + if(increasing) then + if (yout(j).le.yin(1)) then + jjm(j) = 1 + jjp(j) = 1 + wgts(j) = 1._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + else if (yout(j).gt.yin(nin)) then + jjm(j) = nin + jjp(j) = nin + wgts(j) = 1._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + end if + else + if (yout(j).gt.yin(1)) then + jjm(j) = 1 + jjp(j) = 1 + wgts(j) = 1._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + else if (yout(j).le.yin(nin)) then + jjm(j) = nin + jjp(j) = nin + wgts(j) = 1._r8 + wgtn(j) = 0._r8 + extrap = extrap + 1._r8 + end if + end if + end do + case (extrap_method_cycle) + ! + ! For values which extend beyond boundaries, set weights + ! for circular boundaries + ! + dyinwrap = yin(1) + (cmax-cmin) - yin(nin) + avgdyin = abs(yin(nin)-yin(1))/(nin-1._r8) + ratio = dyinwrap/avgdyin + if (ratio < 0.9_r8 .or. ratio > 1.1_r8) then + write(iulog,*) 'Lininterp: Bad dyinwrap value =',dyinwrap,& + ' avg=', avgdyin, yin(1),yin(nin) + call endrun('interpolate_data') + end if + + do j=1,nout + if(increasing) then + if (yout(j) <= yin(1)) then + jjm(j) = nin + jjp(j) = 1 + wgts(j) = (yin(1)-yout(j))/dyinwrap + wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap + else if (yout(j) > yin(nin)) then + jjm(j) = nin + jjp(j) = 1 + wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap + wgtn(j) = (yout(j)-yin(nin))/dyinwrap + end if + else + if (yout(j) > yin(1)) then + jjm(j) = nin + jjp(j) = 1 + wgts(j) = (yin(1)-yout(j))/dyinwrap + wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap + else if (yout(j) <= yin(nin)) then + jjm(j) = nin + jjp(j) = 1 + wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap + wgtn(j) = (yout(j)+(cmax-cmin)-yin(nin))/dyinwrap + end if + + endif + end do + end select + + ! + ! Loop though output indices finding input indices and weights + ! + if(increasing) then + do j=1,nout + do jj=1,nin-1 + if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then + jjm(j) = jj + jjp(j) = jj + 1 + wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj)) + wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj)) + exit + end if + end do + end do + else + do j=1,nout + do jj=1,nin-1 + if (yout(j).le.yin(jj) .and. yout(j).gt.yin(jj+1)) then + jjm(j) = jj + jjp(j) = jj + 1 + wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj)) + wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj)) + exit + end if + end do + end do + end if + + ! + ! Check that interp/extrap points have been found for all outputs + ! + icount = 0 + do j=1,nout + if (jjm(j).eq.0 .or. jjp(j).eq.0) icount = icount + 1 + ratio=wgts(j)+wgtn(j) + if((ratio<0.9_r8.or.ratio>1.1_r8).and.extrap_method.ne.0) then + write(iulog,*) j, wgts(j),wgtn(j),jjm(j),jjp(j), increasing,extrap_method + call endrun('Bad weight computed in LININTERP_init') + end if + end do + if (icount.gt.0) then + call endrun('LININTERP: Point found without interp indices') + end if + + end subroutine lininterp_init + + subroutine lininterp1d (arrin, n1, arrout, m1, interp_wgts) + !----------------------------------------------------------------------- + ! + ! Purpose: Do a linear interpolation from input mesh to output + ! mesh with weights as set in lininterp_init. + ! + ! + ! Author: Jim Edwards + ! + !----------------------------------------------------------------------- + !----------------------------------------------------------------------- + implicit none + !----------------------------------------------------------------------- + ! + ! Arguments + ! + integer, intent(in) :: n1 ! number of input latitudes + integer, intent(in) :: m1 ! number of output latitudes + + real(r8), intent(in) :: arrin(n1) ! input array of values to interpolate + type(interp_type), intent(in) :: interp_wgts + real(r8), intent(out) :: arrout(m1) ! interpolated array + + ! + ! Local workspace + ! + integer j ! latitude indices + integer, pointer :: jjm(:) + integer, pointer :: jjp(:) + + real(r8), pointer :: wgts(:) + real(r8), pointer :: wgtn(:) + + + jjm => interp_wgts%jjm + jjp => interp_wgts%jjp + wgts => interp_wgts%wgts + wgtn => interp_wgts%wgtn + + ! + ! Do the interpolation + ! + do j=1,m1 + arrout(j) = arrin(jjm(j))*wgts(j) + arrin(jjp(j))*wgtn(j) + end do + + return + end subroutine lininterp1d + + subroutine lininterp2d2d(arrin, n1, n2, arrout, m1, m2, wgt1, wgt2) + implicit none + !----------------------------------------------------------------------- + ! + ! Arguments + ! + integer, intent(in) :: n1, n2, m1, m2 + real(r8), intent(in) :: arrin(n1,n2) ! input array of values to interpolate + type(interp_type), intent(in) :: wgt1, wgt2 + real(r8), intent(out) :: arrout(m1,m2) ! interpolated array + ! + ! locals + ! + integer i,j ! indices + integer, pointer :: iim(:), jjm(:) + integer, pointer :: iip(:), jjp(:) + + real(r8), pointer :: wgts1(:), wgts2(:) + real(r8), pointer :: wgtn1(:), wgtn2(:) + + real(r8) :: arrtmp(n1,m2) + + + jjm => wgt2%jjm + jjp => wgt2%jjp + wgts2 => wgt2%wgts + wgtn2 => wgt2%wgtn + + iim => wgt1%jjm + iip => wgt1%jjp + wgts1 => wgt1%wgts + wgtn1 => wgt1%wgtn + + do j=1,m2 + do i=1,n1 + arrtmp(i,j) = arrin(i,jjm(j))*wgts2(j) + arrin(i,jjp(j))*wgtn2(j) + end do + end do + + do j=1,m2 + do i=1,m1 + arrout(i,j) = arrtmp(iim(i),j)*wgts1(i) + arrtmp(iip(i),j)*wgtn1(i) + end do + end do + + end subroutine lininterp2d2d + subroutine lininterp2d1d(arrin, n1, n2, arrout, m1, wgt1, wgt2, fldname) + implicit none + !----------------------------------------------------------------------- + ! + ! Arguments + ! + integer, intent(in) :: n1, n2, m1 + real(r8), intent(in) :: arrin(n1,n2) ! input array of values to interpolate + type(interp_type), intent(in) :: wgt1, wgt2 + real(r8), intent(out) :: arrout(m1) ! interpolated array + character(len=*), intent(in), optional :: fldname(:) + ! + ! locals + ! + integer i ! indices + integer, pointer :: iim(:), jjm(:) + integer, pointer :: iip(:), jjp(:) + + real(r8), pointer :: wgts(:), wgte(:) + real(r8), pointer :: wgtn(:), wgtw(:) + + jjm => wgt2%jjm + jjp => wgt2%jjp + wgts => wgt2%wgts + wgtn => wgt2%wgtn + + iim => wgt1%jjm + iip => wgt1%jjp + wgtw => wgt1%wgts + wgte => wgt1%wgtn + + do i=1,m1 + arrout(i) = arrin(iim(i),jjm(i))*wgtw(i)*wgts(i)+arrin(iip(i),jjm(i))*wgte(i)*wgts(i) + & + arrin(iim(i),jjp(i))*wgtw(i)*wgtn(i)+arrin(iip(i),jjp(i))*wgte(i)*wgtn(i) + end do + + + end subroutine lininterp2d1d + subroutine lininterp3d2d(arrin, n1, n2, n3, arrout, m1, len1, wgt1, wgt2) + implicit none + !----------------------------------------------------------------------- + ! + ! Arguments + ! + integer, intent(in) :: n1, n2, n3, m1, len1 ! m1 is to len1 as ncols is to pcols + real(r8), intent(in) :: arrin(n1,n2,n3) ! input array of values to interpolate + type(interp_type), intent(in) :: wgt1, wgt2 + real(r8), intent(out) :: arrout(len1, n3) ! interpolated array + + ! + ! locals + ! + integer i, k ! indices + integer, pointer :: iim(:), jjm(:) + integer, pointer :: iip(:), jjp(:) + + real(r8), pointer :: wgts(:), wgte(:) + real(r8), pointer :: wgtn(:), wgtw(:) + + jjm => wgt2%jjm + jjp => wgt2%jjp + wgts => wgt2%wgts + wgtn => wgt2%wgtn + + iim => wgt1%jjm + iip => wgt1%jjp + wgtw => wgt1%wgts + wgte => wgt1%wgtn + + do k=1,n3 + do i=1,m1 + arrout(i,k) = arrin(iim(i),jjm(i),k)*wgtw(i)*wgts(i)+arrin(iip(i),jjm(i),k)*wgte(i)*wgts(i) + & + arrin(iim(i),jjp(i),k)*wgtw(i)*wgtn(i)+arrin(iip(i),jjp(i),k)*wgte(i)*wgtn(i) + end do + end do + + end subroutine lininterp3d2d + + + + + subroutine lininterp_finish(interp_wgts) + type(interp_type) :: interp_wgts + + deallocate(interp_wgts%jjm, & + interp_wgts%jjp, & + interp_wgts%wgts, & + interp_wgts%wgtn) + + nullify(interp_wgts%jjm, & + interp_wgts%jjp, & + interp_wgts%wgts, & + interp_wgts%wgtn) + end subroutine lininterp_finish + + subroutine lininterp_original (arrin, yin, nlev, nlatin, arrout, & + yout, nlatout) + !----------------------------------------------------------------------- + ! + ! Purpose: Do a linear interpolation from input mesh defined by yin to output + ! mesh defined by yout. Where extrapolation is necessary, values will + ! be copied from the extreme edge of the input grid. Vectorization is over + ! the vertical (nlev) dimension. + ! + ! Method: Check validity of input, then determine weights, then do the N-S interpolation. + ! + ! Author: Jim Rosinski + ! Modified: Jim Edwards so that there is no requirement of monotonacity on the yout array + ! + !----------------------------------------------------------------------- + implicit none + !----------------------------------------------------------------------- + ! + ! Arguments + ! + integer, intent(in) :: nlev ! number of vertical levels + integer, intent(in) :: nlatin ! number of input latitudes + integer, intent(in) :: nlatout ! number of output latitudes + + real(r8), intent(in) :: arrin(nlev,nlatin) ! input array of values to interpolate + real(r8), intent(in) :: yin(nlatin) ! input mesh + real(r8), intent(in) :: yout(nlatout) ! output mesh + + real(r8), intent(out) :: arrout(nlev,nlatout) ! interpolated array + ! + ! Local workspace + ! + integer j, jj ! latitude indices + integer jjprev ! latitude indices + integer k ! level index + integer icount ! number of values + + real(r8) extrap ! percent grid non-overlap + ! + ! Dynamic + ! + integer :: jjm(nlatout) + integer :: jjp(nlatout) + + real(r8) :: wgts(nlatout) + real(r8) :: wgtn(nlatout) + ! + ! Check validity of input coordinate arrays: must be monotonically increasing, + ! and have a total of at least 2 elements + ! + if (nlatin.lt.2) then + call endrun('LININTERP: Must have at least 2 input points for interpolation') + end if + + icount = 0 + do j=1,nlatin-1 + if (yin(j).gt.yin(j+1)) icount = icount + 1 + end do + + + if (icount.gt.0) then + call endrun('LININTERP: Non-monotonic coordinate array(s) found') + end if + ! + ! Initialize index arrays for later checking + ! + do j=1,nlatout + jjm(j) = 0 + jjp(j) = 0 + end do + ! + ! For values which extend beyond N and S boundaries, set weights + ! such that values will just be copied. + ! + extrap = 0._r8 + + do j=1,nlatout + if (yout(j).le.yin(1)) then + jjm(j) = 1 + jjp(j) = 1 + wgts(j) = 1._r8 + wgtn(j) = 0._r8 + extrap=extrap+1._r8 + else if (yout(j).gt.yin(nlatin)) then + jjm(j) = nlatin + jjp(j) = nlatin + wgts(j) = 1._r8 + wgtn(j) = 0._r8 + extrap=extrap+1._r8 + endif + end do + + ! + ! Loop though output indices finding input indices and weights + ! + do j=1,nlatout + do jj=1,nlatin-1 + if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then + jjm(j) = jj + jjp(j) = jj + 1 + wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj)) + wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj)) + exit + end if + end do + end do + ! + ! Check that interp/extrap points have been found for all outputs + ! + icount = 0 + do j=1,nlatout + if (jjm(j).eq.0 .or. jjp(j).eq.0) then + icount = icount + 1 + end if + end do + if (icount.gt.0) then + call endrun('LININTERP: Point found without interp indices') + end if + ! + ! Do the interpolation + ! + do j=1,nlatout + do k=1,nlev + arrout(k,j) = arrin(k,jjm(j))*wgts(j) + arrin(k,jjp(j))*wgtn(j) + end do + end do + + return + end subroutine lininterp_original + + + subroutine bilin (arrin, xin, yin, nlondin, nlonin, & + nlevdin, nlev, nlatin, arrout, xout, & + yout, nlondout, nlonout, nlevdout, nlatout) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! + ! Do a bilinear interpolation from input mesh defined by xin, yin to output + ! mesh defined by xout, yout. Circularity is assumed in the x-direction so + ! input x-grid must be in degrees east and must start from Greenwich. When + ! extrapolation is necessary in the N-S direction, values will be copied + ! from the extreme edge of the input grid. Vectorization is over the + ! longitude dimension. Input grid is assumed rectangular. Output grid + ! is assumed ragged, i.e. xout is a 2-d mesh. + ! + ! Method: Interpolate first in longitude, then in latitude. + ! + ! Author: Jim Rosinski + ! + !----------------------------------------------------------------------- + use shr_kind_mod, only: r8 => shr_kind_r8 + use cam_abortutils, only: endrun + !----------------------------------------------------------------------- + implicit none + !----------------------------------------------------------------------- + ! + ! Input arguments + ! + integer, intent(in) :: nlondin ! longitude dimension of input grid + integer, intent(in) :: nlonin ! number of real longitudes (input) + integer, intent(in) :: nlevdin ! vertical dimension of input grid + integer, intent(in) :: nlev ! number of vertical levels + integer, intent(in) :: nlatin ! number of input latitudes + integer, intent(in) :: nlatout ! number of output latitudes + integer, intent(in) :: nlondout ! longitude dimension of output grid + integer, intent(in) :: nlonout(nlatout) ! number of output longitudes per lat + integer, intent(in) :: nlevdout ! vertical dimension of output grid + + real(r8), intent(in) :: arrin(nlondin,nlevdin,nlatin) ! input array of values to interpolate + real(r8), intent(in) :: xin(nlondin) ! input x mesh + real(r8), intent(in) :: yin(nlatin) ! input y mesh + real(r8), intent(in) :: xout(nlondout,nlatout) ! output x mesh + real(r8), intent(in) :: yout(nlatout) ! output y mesh + ! + ! Output arguments + ! + real(r8), intent(out) :: arrout(nlondout,nlevdout,nlatout) ! interpolated array + ! + ! Local workspace + ! + integer :: i, ii, iw, ie, iiprev ! longitude indices + integer :: j, jj, js, jn, jjprev ! latitude indices + integer :: k ! level index + integer :: icount ! number of bad values + + real(r8) :: extrap ! percent grid non-overlap + real(r8) :: dxinwrap ! delta-x on input grid for 2-pi + real(r8) :: avgdxin ! avg input delta-x + real(r8) :: ratio ! compare dxinwrap to avgdxin + real(r8) :: sum ! sum of weights (used for testing) + ! + ! Dynamic + ! + integer :: iim(nlondout) ! interpolation index to the left + integer :: iip(nlondout) ! interpolation index to the right + integer :: jjm(nlatout) ! interpolation index to the south + integer :: jjp(nlatout) ! interpolation index to the north + + real(r8) :: wgts(nlatout) ! interpolation weight to the north + real(r8) :: wgtn(nlatout) ! interpolation weight to the north + real(r8) :: wgte(nlondout) ! interpolation weight to the north + real(r8) :: wgtw(nlondout) ! interpolation weight to the north + real(r8) :: igrid(nlonin) ! interpolation weight to the north + ! + ! Check validity of input coordinate arrays: must be monotonically increasing, + ! and have a total of at least 2 elements + ! + if (nlonin < 2 .or. nlatin < 2) then + call endrun ('BILIN: Must have at least 2 input points for interpolation') + end if + + if (xin(1) < 0._r8 .or. xin(nlonin) > 360._r8) then + call endrun ('BILIN: Input x-grid must be between 0 and 360') + end if + + icount = 0 + do i=1,nlonin-1 + if (xin(i) >= xin(i+1)) icount = icount + 1 + end do + + do j=1,nlatin-1 + if (yin(j) >= yin(j+1)) icount = icount + 1 + end do + + do j=1,nlatout-1 + if (yout(j) >= yout(j+1)) icount = icount + 1 + end do + + do j=1,nlatout + do i=1,nlonout(j)-1 + if (xout(i,j) >= xout(i+1,j)) icount = icount + 1 + end do + end do + + if (icount > 0) then + call endrun ('BILIN: Non-monotonic coordinate array(s) found') + end if + + if (yout(nlatout) <= yin(1) .or. yout(1) >= yin(nlatin)) then + call endrun ('BILIN: No overlap in y-coordinates') + end if + + do j=1,nlatout + if (xout(1,j) < 0._r8 .or. xout(nlonout(j),j) > 360._r8) then + call endrun ('BILIN: Output x-grid must be between 0 and 360') + end if + + if (xout(nlonout(j),j) <= xin(1) .or. & + xout(1,j) >= xin(nlonin)) then + call endrun ('BILIN: No overlap in x-coordinates') + end if + end do + ! + ! Initialize index arrays for later checking + ! + do j=1,nlatout + jjm(j) = 0 + jjp(j) = 0 + end do + ! + ! For values which extend beyond N and S boundaries, set weights + ! such that values will just be copied. + ! + do js=1,nlatout + if (yout(js) > yin(1)) exit + jjm(js) = 1 + jjp(js) = 1 + wgts(js) = 1._r8 + wgtn(js) = 0._r8 + end do + + do jn=nlatout,1,-1 + if (yout(jn) <= yin(nlatin)) exit + jjm(jn) = nlatin + jjp(jn) = nlatin + wgts(jn) = 1._r8 + wgtn(jn) = 0._r8 + end do + ! + ! Loop though output indices finding input indices and weights + ! + jjprev = 1 + do j=js,jn + do jj=jjprev,nlatin-1 + if (yout(j) > yin(jj) .and. yout(j) <= yin(jj+1)) then + jjm(j) = jj + jjp(j) = jj + 1 + wgts(j) = (yin(jj+1) - yout(j)) / (yin(jj+1) - yin(jj)) + wgtn(j) = (yout(j) - yin(jj)) / (yin(jj+1) - yin(jj)) + goto 30 + end if + end do + call endrun ('BILIN: Failed to find interp values') +30 jjprev = jj + end do + + dxinwrap = xin(1) + 360._r8 - xin(nlonin) + ! + ! Check for sane dxinwrap values. Allow to differ no more than 10% from avg + ! + avgdxin = (xin(nlonin)-xin(1))/(nlonin-1._r8) + ratio = dxinwrap/avgdxin + if (ratio < 0.9_r8 .or. ratio > 1.1_r8) then + write(iulog,*)'BILIN: Insane dxinwrap value =',dxinwrap,' avg=', avgdxin + call endrun + end if + ! + ! Check that interp/extrap points have been found for all outputs, and that + ! they are reasonable (i.e. within 32-bit roundoff) + ! + icount = 0 + do j=1,nlatout + if (jjm(j) == 0 .or. jjp(j) == 0) icount = icount + 1 + sum = wgts(j) + wgtn(j) + if (sum < 0.99999_r8 .or. sum > 1.00001_r8) icount = icount + 1 + if (wgts(j) < 0._r8 .or. wgts(j) > 1._r8) icount = icount + 1 + if (wgtn(j) < 0._r8 .or. wgtn(j) > 1._r8) icount = icount + 1 + end do + + if (icount > 0) then + call endrun ('BILIN: something bad in latitude indices or weights') + end if + ! + ! Do the bilinear interpolation + ! + do j=1,nlatout + ! + ! Initialize index arrays for later checking + ! + do i=1,nlondout + iim(i) = 0 + iip(i) = 0 + end do + ! + ! For values which extend beyond E and W boundaries, set weights such that + ! values will be interpolated between E and W edges of input grid. + ! + do iw=1,nlonout(j) + if (xout(iw,j) > xin(1)) exit + iim(iw) = nlonin + iip(iw) = 1 + wgtw(iw) = (xin(1) - xout(iw,j)) /dxinwrap + wgte(iw) = (xout(iw,j)+360._r8 - xin(nlonin))/dxinwrap + end do + + do ie=nlonout(j),1,-1 + if (xout(ie,j) <= xin(nlonin)) exit + iim(ie) = nlonin + iip(ie) = 1 + wgtw(ie) = (xin(1)+360._r8 - xout(ie,j)) /dxinwrap + wgte(ie) = (xout(ie,j) - xin(nlonin))/dxinwrap + end do + ! + ! Loop though output indices finding input indices and weights + ! + iiprev = 1 + do i=iw,ie + do ii=iiprev,nlonin-1 + if (xout(i,j) > xin(ii) .and. xout(i,j) <= xin(ii+1)) then + iim(i) = ii + iip(i) = ii + 1 + wgtw(i) = (xin(ii+1) - xout(i,j)) / (xin(ii+1) - xin(ii)) + wgte(i) = (xout(i,j) - xin(ii)) / (xin(ii+1) - xin(ii)) + goto 60 + end if + end do + call endrun ('BILIN: Failed to find interp values') +60 iiprev = ii + end do + + icount = 0 + do i=1,nlonout(j) + if (iim(i) == 0 .or. iip(i) == 0) icount = icount + 1 + sum = wgtw(i) + wgte(i) + if (sum < 0.99999_r8 .or. sum > 1.00001_r8) icount = icount + 1 + if (wgtw(i) < 0._r8 .or. wgtw(i) > 1._r8) icount = icount + 1 + if (wgte(i) < 0._r8 .or. wgte(i) > 1._r8) icount = icount + 1 + end do + + if (icount > 0) then + write(iulog,*)'BILIN: j=',j,' Something bad in longitude indices or weights' + call endrun + end if + ! + ! Do the interpolation, 1st in longitude then latitude + ! + do k=1,nlev + do i=1,nlonin + igrid(i) = arrin(i,k,jjm(j))*wgts(j) + arrin(i,k,jjp(j))*wgtn(j) + end do + + do i=1,nlonout(j) + arrout(i,k,j) = igrid(iim(i))*wgtw(i) + igrid(iip(i))*wgte(i) + end do + end do + end do + + + return + end subroutine bilin + +!========================================================================================= + +subroutine vertinterp(ncol, ncold, nlev, pin, pout, arrin, arrout, & + extrapolate, ln_interp, ps, phis, tbot) + + ! Vertically interpolate input array to output pressure level. The + ! interpolation is linear in either pressure (the default) or ln pressure. + ! + ! If above the top boundary then return top boundary value. + ! + ! If below bottom boundary then use bottom boundary value, or optionally + ! for T or Z use the extrapolation algorithm from ECMWF (which was taken + ! from the NCL code base). + + use physconst, only: rair, rga + + !------------------------------Arguments-------------------------------- + integer, intent(in) :: ncol ! number active columns + integer, intent(in) :: ncold ! column dimension + integer, intent(in) :: nlev ! vertical dimension + real(r8), intent(in) :: pin(ncold,nlev) ! input pressure levels + real(r8), intent(in) :: pout ! output pressure level + real(r8), intent(in) :: arrin(ncold,nlev) ! input array + real(r8), intent(out) :: arrout(ncold) ! output array (interpolated) + + character(len=1), optional, intent(in) :: extrapolate ! either 'T' or 'Z' + logical, optional, intent(in) :: ln_interp ! set true to interolate + ! in ln(pressure) + real(r8), optional, intent(in) :: ps(ncold) ! surface pressure + real(r8), optional, intent(in) :: phis(ncold) ! surface geopotential + real(r8), optional, intent(in) :: tbot(ncold) ! temperature at bottom level + + !---------------------------Local variables----------------------------- + real(r8) :: alpha + logical :: linear_interp + logical :: do_extrapolate_T, do_extrapolate_Z + + integer :: i,k ! indices + integer :: kupper(ncold) ! Level indices for interpolation + real(r8) :: dpu ! upper level pressure difference + real(r8) :: dpl ! lower level pressure difference + logical :: found(ncold) ! true if input levels found + logical :: error ! error flag + !---------------------------------------------------------------------------- + + alpha = 0.0065_r8*rair*rga + + do_extrapolate_T = .false. + do_extrapolate_Z = .false. + if (present(extrapolate)) then + + if (extrapolate == 'T') do_extrapolate_T = .true. + if (extrapolate == 'Z') do_extrapolate_Z = .true. + + if (.not. do_extrapolate_T .and. .not. do_extrapolate_Z) then + call endrun ('VERTINTERP: extrapolate must be T or Z') + end if + if (.not. present(ps)) then + call endrun ('VERTINTERP: ps required for extrapolation') + end if + if (.not. present(phis)) then + call endrun ('VERTINTERP: phis required for extrapolation') + end if + if (do_extrapolate_Z) then + if (.not. present(tbot)) then + call endrun ('VERTINTERP: extrapolate must be T or Z') + end if + end if + end if + + linear_interp = .true. + if (present(ln_interp)) then + if (ln_interp) linear_interp = .false. + end if + + do i = 1, ncol + found(i) = .false. + kupper(i) = 1 + end do + error = .false. + + ! Find indices of upper pressure bound + do k = 1, nlev - 1 + do i = 1, ncol + if ((.not. found(i)) .and. pin(i,k)= pin(i,nlev)) then + + if (do_extrapolate_T) then + ! use ECMWF algorithm to extrapolate T + arrout(i) = extrapolate_T() + + else if (do_extrapolate_Z) then + ! use ECMWF algorithm to extrapolate Z + arrout(i) = extrapolate_Z() + + else + ! use bottom value + arrout(i) = arrin(i,nlev) + end if + + else if (found(i)) then + if (linear_interp) then + ! linear interpolation in p + dpu = pout - pin(i,kupper(i)) + dpl = pin(i,kupper(i)+1) - pout + arrout(i) = (arrin(i,kupper(i))*dpl + arrin(i,kupper(i)+1)*dpu)/(dpl + dpu) + else + ! linear interpolation in ln(p) + arrout(i) = arrin(i,kupper(i)) + (arrin(i,kupper(i)+1) - arrin(i,kupper(i)))* & + log(pout/pin(i,kupper(i))) / & + log(pin(i,kupper(i)+1)/pin(i,kupper(i))) + end if + else + error = .true. + end if + end do + + ! Error check + if (error) then + call endrun ('VERTINTERP: ERROR FLAG') + end if + +contains + + !====================================================================================== + + real(r8) function extrapolate_T() + + ! local variables + real(r8) :: sixth + real(r8) :: tstar + real(r8) :: hgt + real(r8) :: alnp + real(r8) :: t0 + real(r8) :: tplat + real(r8) :: tprime0 + !------------------------------------------------------------------------- + + sixth = 1._r8 / 6._r8 + tstar = arrin(i,nlev)*(1._r8 + alpha*(ps(i)/pin(i,nlev) - 1._r8)) + hgt = phis(i)*rga + + if (hgt < 2000._r8) then + alnp = alpha*log(pout/ps(i)) + else + t0 = tstar + 0.0065_r8*hgt + tplat = min(t0, 298._r8) + + if (hgt <= 2500._r8) then + tprime0 = 0.002_r8*((2500._r8 - hgt)*t0 + & + (hgt - 2000._r8)*tplat) + else + tprime0 = tplat + end if + + if (tprime0 < tstar) then + alnp = 0._r8 + else + alnp = rair*(tprime0 - tstar)/phis(i) * log(pout/ps(i)) + end if + + end if + + extrapolate_T = tstar*(1._r8 + alnp + 0.5_r8*alnp**2 + sixth*alnp**3) + + end function extrapolate_T + + !====================================================================================== + + real(r8) function extrapolate_Z() + + ! local variables + real(r8) :: sixth + real(r8) :: hgt + real(r8) :: tstar + real(r8) :: t0 + real(r8) :: alph + real(r8) :: alnp + !------------------------------------------------------------------------- + + sixth = 1._r8 / 6._r8 + hgt = phis(i)*rga + tstar = tbot(i)*(1._r8 + alpha*(ps(i)/pin(i,nlev) - 1._r8)) + t0 = tstar + 0.0065_r8*hgt + + if (tstar <= 290.5_r8 .and. t0 > 290.5_r8) then + alph = rair/phis(i) * (290.5_r8 - tstar) + + else if (tstar > 290.5_r8 .and. t0 > 290.5_r8) then + alph = 0._r8 + tstar = 0.5_r8*(290.5_r8 + tstar) + + else + alph = alpha + + end if + + if (tstar < 255._r8) then + tstar = 0.5_r8*(tstar + 255._r8) + end if + alnp = alph*log(pout/ps(i)) + + extrapolate_Z = hgt - rair*tstar*rga*log(pout/ps(i))* & + (1._r8 + 0.5_r8*alnp + sixth*alnp**2) + + end function extrapolate_Z + +end subroutine vertinterp + +!========================================================================================= + + subroutine get_timeinterp_factors (cycflag, np1, cdayminus, cdayplus, cday, & + fact1, fact2, str) + !--------------------------------------------------------------------------- + ! + ! Purpose: Determine time interpolation factors (normally for a boundary dataset) + ! for linear interpolation. + ! + ! Method: Assume 365 days per year. Output variable fact1 will be the weight to + ! apply to data at calendar time "cdayminus", and fact2 the weight to apply + ! to data at time "cdayplus". Combining these values will produce a result + ! valid at time "cday". Output arguments fact1 and fact2 will be between + ! 0 and 1, and fact1 + fact2 = 1 to roundoff. + ! + ! Author: Jim Rosinski + ! + !--------------------------------------------------------------------------- + implicit none + ! + ! Arguments + ! + logical, intent(in) :: cycflag ! flag indicates whether dataset is being cycled yearly + + integer, intent(in) :: np1 ! index points to forward time slice matching cdayplus + + real(r8), intent(in) :: cdayminus ! calendar day of rearward time slice + real(r8), intent(in) :: cdayplus ! calendar day of forward time slice + real(r8), intent(in) :: cday ! calenar day to be interpolated to + real(r8), intent(out) :: fact1 ! time interpolation factor to apply to rearward time slice + real(r8), intent(out) :: fact2 ! time interpolation factor to apply to forward time slice + + character(len=*), intent(in) :: str ! string to be added to print in case of error (normally the callers name) + ! + ! Local workspace + ! + real(r8) :: deltat ! time difference (days) between cdayminus and cdayplus + real(r8), parameter :: daysperyear = 365._r8 ! number of days in a year + ! + ! Initial sanity checks + ! + if (np1 == 1 .and. .not. cycflag) then + call endrun ('GETFACTORS:'//str//' cycflag false and forward month index = Jan. not allowed') + end if + + if (np1 < 1) then + call endrun ('GETFACTORS:'//str//' input arg np1 must be > 0') + end if + + if (cycflag) then + if ((cday < 1._r8) .or. (cday > (daysperyear+1._r8))) then + write(iulog,*) 'GETFACTORS:', str, ' bad cday=',cday + call endrun () + end if + else + if (cday < 1._r8) then + write(iulog,*) 'GETFACTORS:', str, ' bad cday=',cday + call endrun () + end if + end if + ! + ! Determine time interpolation factors. Account for December-January + ! interpolation if dataset is being cycled yearly. + ! + if (cycflag .and. np1 == 1) then ! Dec-Jan interpolation + deltat = cdayplus + daysperyear - cdayminus + if (cday > cdayplus) then ! We are in December + fact1 = (cdayplus + daysperyear - cday)/deltat + fact2 = (cday - cdayminus)/deltat + else ! We are in January + fact1 = (cdayplus - cday)/deltat + fact2 = (cday + daysperyear - cdayminus)/deltat + end if + else + deltat = cdayplus - cdayminus + fact1 = (cdayplus - cday)/deltat + fact2 = (cday - cdayminus)/deltat + end if + + if (.not. valid_timeinterp_factors (fact1, fact2)) then + write(iulog,*) 'GETFACTORS: ', str, ' bad fact1 and/or fact2=', fact1, fact2 + call endrun () + end if + + return + end subroutine get_timeinterp_factors + + logical function valid_timeinterp_factors (fact1, fact2) + !--------------------------------------------------------------------------- + ! + ! Purpose: check sanity of time interpolation factors to within 32-bit roundoff + ! + !--------------------------------------------------------------------------- + implicit none + + real(r8), intent(in) :: fact1, fact2 ! time interpolation factors + + valid_timeinterp_factors = .true. + + ! The fact1 .ne. fact1 and fact2 .ne. fact2 comparisons are to detect NaNs. + if (abs(fact1+fact2-1._r8) > 1.e-6_r8 .or. & + fact1 > 1.000001_r8 .or. fact1 < -1.e-6_r8 .or. & + fact2 > 1.000001_r8 .or. fact2 < -1.e-6_r8 .or. & + fact1 .ne. fact1 .or. fact2 .ne. fact2) then + + valid_timeinterp_factors = .false. + end if + + return + end function valid_timeinterp_factors + +end module interpolate_data From a8986d88d139f1e06fd1b5f74bc9a2dc750f42ea Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 15 Aug 2024 14:23:41 -0400 Subject: [PATCH 02/25] Add cam_comp bindings for test build --- src/control/cam_comp.F90 | 16 ++++++ src/data/registry.xml | 63 ++++++++++++++++++++++ src/physics/utils/tropopause_read_file.F90 | 20 ++++--- 3 files changed, 88 insertions(+), 11 deletions(-) diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index db0e943e..a5829e00 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -23,10 +23,12 @@ module cam_comp use time_manager, only: timemgr_init, get_step_size use time_manager, only: get_nstep use time_manager, only: is_first_step, is_first_restart_step + use time_manager, only: get_curr_calday use camsrfexch, only: cam_out_t, cam_in_t use physics_types, only: phys_state, phys_tend use physics_types, only: dtime_phys + use physics_types, only: calday use dyn_comp, only: dyn_import_t, dyn_export_t use perf_mod, only: t_barrierf, t_startf, t_stopf @@ -98,6 +100,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & use physics_grid, only: columns_on_task use vert_coord, only: pver use phys_vars_init_check, only: mark_as_initialized + use tropopause_read_file, only: tropopause_read_file ! Arguments character(len=cl), intent(in) :: caseid ! case ID @@ -168,6 +171,10 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & ref_tod, stop_ymd, stop_tod, curr_ymd, curr_tod, & perpetual_run, perpetual_ymd, initial_run_in) + ! Get current fractional calendar day. Needs to be updated at every timestep. + calday = get_curr_calday() + mark_as_initialized('calday') + ! Read CAM namelists. filein = "atm_in" // trim(inst_suffix) call read_namelist(filein, single_column, scmlat, scmlon) @@ -224,6 +231,12 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & !!XXgoldyXX: ^ need to import this end if + ! Read tropopause climatology + ! FIXME hplin 8/15/24: needs to get data from tropopause_nl, how to pass from CCPP? + call tropopause_read_file('/fs/cgd/csm/inputdata/atm/cam/chem/trop_mozart/ub/clim_p_trop.nc') + mark_as_initialized('tropp_p_loc') + mark_as_initialized('tropp_days') + call phys_init() !!XXgoldyXX: v need to import this @@ -270,6 +283,9 @@ subroutine cam_timestep_init() ! call phys_timestep_init() + ! Update current fractional calendar day. Needs to be updated at every timestep. + calday = get_curr_calday() + end subroutine cam_timestep_init ! !----------------------------------------------------------------------- diff --git a/src/data/registry.xml b/src/data/registry.xml index 57a9ae3b..3c232fd5 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -326,6 +326,69 @@ current timestep number 0 + + + fractional calendar day at end of current timestep + + + + Climatological tropopause pressures from file + horizontal_dimension 12 + + + Climatological tropopause calendar day indices from file + 12 + + + Tropopause primary find method + 5 + + + + Tropopause secondary find method + 3 + + + + Tropopause model level + horizontal_dimension + + + Tropopause level pressure + horizontal_dimension + + + Tropopause level temperature + horizontal_dimension + + + Tropopause level altitude + horizontal_dimension + Date: Thu, 15 Aug 2024 15:37:45 -0400 Subject: [PATCH 03/25] Add metadata for tropopause_read_file for inclusion in registry --- src/data/registry.xml | 15 +------------- src/physics/utils/tropopause_read_file.F90 | 14 ++++++++++--- src/physics/utils/tropopause_read_file.meta | 23 +++++++++++++++++++++ 3 files changed, 35 insertions(+), 17 deletions(-) create mode 100644 src/physics/utils/tropopause_read_file.meta diff --git a/src/data/registry.xml b/src/data/registry.xml index 3c232fd5..e88a0959 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -10,6 +10,7 @@ $SRCROOT/src/data/physconst.meta $SRCROOT/src/physics/utils/physics_grid.meta $SRCROOT/src/physics/utils/cam_constituents.meta + $SRCROOT/src/physics/utils/tropopause_read_file.meta $SRCROOT/src/data/air_composition.meta $SRCROOT/src/data/cam_thermo.meta $SRCROOT/src/data/ref_pres.meta @@ -333,20 +334,6 @@ fractional calendar day at end of current timestep - - Climatological tropopause pressures from file - horizontal_dimension 12 - - - Climatological tropopause calendar day indices from file - 12 - diff --git a/src/physics/utils/tropopause_read_file.F90 b/src/physics/utils/tropopause_read_file.F90 index d0fa6ebb..b73238ee 100644 --- a/src/physics/utils/tropopause_read_file.F90 +++ b/src/physics/utils/tropopause_read_file.F90 @@ -7,14 +7,22 @@ module tropopause_read_file ! for use within CAM-SIMA, particularly removal of chunk support. !------------------------------------------------------------------- - ! climatological tropopause pressures (pcols,1,ntimes) and monthly day-of-year times (12) - use physics_types, only: tropp_p_loc, tropp_days - implicit none private public :: tropopause_read_file +!> \section arg_table_tropopause_read_file Argument Table +!! \htmlinclude tropopause_read_file.html + ! days in year for climatological tropopause pressure data + integer, public, parameter :: tropp_slices = 12 + + ! climatological tropopause pressures (pcols,ntimes) + real(kind_phys), public, pointer :: tropp_p_loc(:,:) + + ! monthly day-of-year times corresponding to climatological data (12) + integer, public, pointer :: tropp_days(:) + contains subroutine tropopause_read_file(tropopause_climo_file) !------------------------------------------------------------------ diff --git a/src/physics/utils/tropopause_read_file.meta b/src/physics/utils/tropopause_read_file.meta new file mode 100644 index 00000000..9545f1ca --- /dev/null +++ b/src/physics/utils/tropopause_read_file.meta @@ -0,0 +1,23 @@ +[ccpp-table-properties] + name = tropopause_read_file + type = module + +[ccpp-arg-table] + name = tropopause_read_file + type = module +[ tropp_slices ] + standard_name = number_of_months_in_year + units = 1 + type = integer + dimensions = () +[ tropp_p_loc ] + standard_name = tropopause_air_pressure_from_climatology + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_dimension, number_of_months_in_year) +[ tropp_days ] + standard_name = tropopause_calendar_days_from_climatology + long_name = Climatological tropopause calendar day indices from file + units = 1 + type = integer + dimensions = (number_of_months_in_year) From c801df4118660623482a94c0104f41a29dd44df1 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 15 Aug 2024 16:05:27 -0400 Subject: [PATCH 04/25] Fix registry units for model_level_number_at_tropopause --- src/data/registry.xml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/data/registry.xml b/src/data/registry.xml index e88a0959..2e532be2 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -336,21 +336,21 @@ + units="1" type="integer"> Tropopause primary find method 5 + units="1" type="integer"> Tropopause secondary find method 3 Tropopause model level horizontal_dimension From 20fc03002b583e899db564b4d95fbfdcfde3e415 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 15 Aug 2024 16:51:13 -0400 Subject: [PATCH 05/25] Add messages for endrun calls in interpolate_data.F90 --- src/utils/interpolate_data.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/utils/interpolate_data.F90 b/src/utils/interpolate_data.F90 index 142e3600..c02f5582 100644 --- a/src/utils/interpolate_data.F90 +++ b/src/utils/interpolate_data.F90 @@ -797,7 +797,7 @@ subroutine bilin (arrin, xin, yin, nlondin, nlonin, & ratio = dxinwrap/avgdxin if (ratio < 0.9_r8 .or. ratio > 1.1_r8) then write(iulog,*)'BILIN: Insane dxinwrap value =',dxinwrap,' avg=', avgdxin - call endrun + call endrun('BILIN: Insane dxinwrap value') end if ! ! Check that interp/extrap points have been found for all outputs, and that @@ -874,7 +874,7 @@ subroutine bilin (arrin, xin, yin, nlondin, nlonin, & if (icount > 0) then write(iulog,*)'BILIN: j=',j,' Something bad in longitude indices or weights' - call endrun + call endrun('BILIN: Something bad in longitude indices or weights') end if ! ! Do the interpolation, 1st in longitude then latitude @@ -1168,12 +1168,12 @@ subroutine get_timeinterp_factors (cycflag, np1, cdayminus, cdayplus, cday, & if (cycflag) then if ((cday < 1._r8) .or. (cday > (daysperyear+1._r8))) then write(iulog,*) 'GETFACTORS:', str, ' bad cday=',cday - call endrun () + call endrun ('get_timeinterp_factors GETFACTORS bad cday (1)') end if else if (cday < 1._r8) then write(iulog,*) 'GETFACTORS:', str, ' bad cday=',cday - call endrun () + call endrun ('get_timeinterp_factors GETFACTORS bad cday (2)') end if end if ! @@ -1197,7 +1197,7 @@ subroutine get_timeinterp_factors (cycflag, np1, cdayminus, cdayplus, cday, & if (.not. valid_timeinterp_factors (fact1, fact2)) then write(iulog,*) 'GETFACTORS: ', str, ' bad fact1 and/or fact2=', fact1, fact2 - call endrun () + call endrun ('get_timeinterp_factors GETFACTORS bad fact1 and/or fact2') end if return From 9c5fc0a188c19646aef9b610202ae1ac91cd59a2 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 15 Aug 2024 17:07:52 -0400 Subject: [PATCH 06/25] Update tropopause_read_file subroutine name to avoid conflict with module naming --- src/control/cam_comp.F90 | 4 +--- src/physics/utils/tropopause_read_file.F90 | 8 +++++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index a5829e00..8c465413 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -233,9 +233,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & ! Read tropopause climatology ! FIXME hplin 8/15/24: needs to get data from tropopause_nl, how to pass from CCPP? - call tropopause_read_file('/fs/cgd/csm/inputdata/atm/cam/chem/trop_mozart/ub/clim_p_trop.nc') - mark_as_initialized('tropp_p_loc') - mark_as_initialized('tropp_days') + call tropopause_read_climo_file('/fs/cgd/csm/inputdata/atm/cam/chem/trop_mozart/ub/clim_p_trop.nc') call phys_init() diff --git a/src/physics/utils/tropopause_read_file.F90 b/src/physics/utils/tropopause_read_file.F90 index b73238ee..23ff6d62 100644 --- a/src/physics/utils/tropopause_read_file.F90 +++ b/src/physics/utils/tropopause_read_file.F90 @@ -7,10 +7,12 @@ module tropopause_read_file ! for use within CAM-SIMA, particularly removal of chunk support. !------------------------------------------------------------------- + use ccpp_kinds, only: kind_phys + implicit none private - public :: tropopause_read_file + public :: tropopause_read_climo_file !> \section arg_table_tropopause_read_file Argument Table !! \htmlinclude tropopause_read_file.html @@ -24,7 +26,7 @@ module tropopause_read_file integer, public, pointer :: tropp_days(:) contains - subroutine tropopause_read_file(tropopause_climo_file) + subroutine tropopause_read_climo_file(tropopause_climo_file) !------------------------------------------------------------------ ! ... initialize upper boundary values !------------------------------------------------------------------ @@ -162,5 +164,5 @@ subroutine tropopause_read_file(tropopause_climo_file) write(iulog,'(1p,5g15.8)') tropp_days(:) endif - end subroutine tropopause_read_file + end subroutine tropopause_read_climo_file end module tropopause_read_file \ No newline at end of file From 6aa43a649fdd2a127b288fd365eca000da235618 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 15 Aug 2024 17:14:37 -0400 Subject: [PATCH 07/25] Update tropopause_read_file to CAM-SIMA subroutines --- src/physics/utils/tropopause_read_file.F90 | 37 ++++++++++++++-------- 1 file changed, 23 insertions(+), 14 deletions(-) diff --git a/src/physics/utils/tropopause_read_file.F90 b/src/physics/utils/tropopause_read_file.F90 index 23ff6d62..4ec967f7 100644 --- a/src/physics/utils/tropopause_read_file.F90 +++ b/src/physics/utils/tropopause_read_file.F90 @@ -7,7 +7,10 @@ module tropopause_read_file ! for use within CAM-SIMA, particularly removal of chunk support. !------------------------------------------------------------------- - use ccpp_kinds, only: kind_phys + use ccpp_kinds, only: kind_phys + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + use spmd_utils, only: masterproc implicit none private @@ -35,7 +38,7 @@ subroutine tropopause_read_climo_file(tropopause_climo_file) use physics_grid, only : pcols => columns_on_task use physconst, only : pi use time_manager, only : get_calday - use ioFileMod, only : getfil + use ioFileMod, only : cam_get_file use cam_pio_utils, only: cam_pio_openfile use pio, only : file_desc_t, var_desc_t, pio_inq_dimid, pio_inq_dimlen, & pio_inq_varid, pio_get_var, pio_closefile, pio_nowrite @@ -66,7 +69,7 @@ subroutine tropopause_read_climo_file(tropopause_climo_file) !----------------------------------------------------------------------- ! ... open netcdf file !----------------------------------------------------------------------- - call getfil (tropopause_climo_file, locfn, 0) + call cam_get_file (tropopause_climo_file, locfn, allow_fail=.false.) call cam_pio_openfile(pio_id, trim(locfn), PIO_NOWRITE) !----------------------------------------------------------------------- @@ -75,8 +78,8 @@ subroutine tropopause_read_climo_file(tropopause_climo_file) ierr = pio_inq_dimid( pio_id, 'time', dimid ) ierr = pio_inq_dimlen( pio_id, dimid, ntimes ) if( ntimes /= 12 )then - write(iulog,*) 'tropopause_find_init: number of months = ',ntimes,'; expecting 12' - call endrun + write(iulog,*) 'tropopause_read_climo_file: number of months = ',ntimes,'; expecting 12' + call endrun('tropopause_read_climo_file: number of months in file incorrect') end if !----------------------------------------------------------------------- ! ... get latitudes @@ -85,8 +88,8 @@ subroutine tropopause_read_climo_file(tropopause_climo_file) ierr = pio_inq_dimlen( pio_id, dimid, nlat ) allocate( lat(nlat), stat=ierr ) if( ierr /= 0 ) then - write(iulog,*) 'tropopause_find_init: lat allocation error = ',ierr - call endrun + write(iulog,*) 'tropopause_read_climo_file: lat allocation error = ',ierr + call endrun('tropopause_read_climo_file: failed to allocate lat') end if ierr = pio_inq_varid( pio_id, 'lat', vid ) ierr = pio_get_var( pio_id, vid, lat ) @@ -98,8 +101,8 @@ subroutine tropopause_read_climo_file(tropopause_climo_file) ierr = pio_inq_dimlen( pio_id, dimid, nlon ) allocate( lon(nlon), stat=ierr ) if( ierr /= 0 ) then - write(iulog,*) 'tropopause_find_init: lon allocation error = ',ierr - call endrun + write(iulog,*) 'tropopause_read_climo_file: lon allocation error = ',ierr + call endrun('tropopause_read_climo_file: failed to allocate lon') end if ierr = pio_inq_varid( pio_id, 'lon', vid ) ierr = pio_get_var( pio_id, vid, lon ) @@ -110,8 +113,8 @@ subroutine tropopause_read_climo_file(tropopause_climo_file) !------------------------------------------------------------------ allocate( tropp_p_in(nlon,nlat,ntimes), stat=ierr ) if( ierr /= 0 ) then - write(iulog,*) 'tropopause_find_init: tropp_p_in allocation error = ',ierr - call endrun + write(iulog,*) 'tropopause_read_climo_file: tropp_p_in allocation error = ',ierr + call endrun('tropopause_read_climo_file: failed to allocate tropp_p_in') end if !------------------------------------------------------------------ ! ... read in the tropopause pressure @@ -135,8 +138,8 @@ subroutine tropopause_read_climo_file(tropopause_climo_file) allocate( tropp_p_loc(pcols,ntimes), stat=ierr ) if( ierr /= 0 ) then - write(iulog,*) 'tropopause_find_init: tropp_p_loc allocation error = ',ierr - call endrun + write(iulog,*) 'tropopause_read_climo_file: tropp_p_loc allocation error = ',ierr + call endrun('tropopause_read_climo_file: failed to allocate tropp_p_loc') end if call get_rlat_all_p(pcols, to_lats) @@ -156,11 +159,17 @@ subroutine tropopause_read_climo_file(tropopause_climo_file) ! ... initialize the monthly day of year times !-------------------------------------------------------- + allocate( tropp_days(tropp_slices), stat=ierr ) + if( ierr /= 0 ) then + write(iulog,*) 'tropopause_read_climo_file: tropp_days allocation error = ',ierr + call endrun('tropopause_read_climo_file: failed to allocate tropp_days') + end if + do n = 1,12 tropp_days(n) = get_calday( dates(n), 0 ) end do if (masterproc) then - write(iulog,*) 'tropopause_find_init : tropp_days' + write(iulog,*) 'tropopause_read_climo_file : tropp_days' write(iulog,'(1p,5g15.8)') tropp_days(:) endif From e442cf52815c984e888726a4819654bfb24bb176 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 15 Aug 2024 17:26:50 -0400 Subject: [PATCH 08/25] Removal of POINTERs per CCPP-convention --- src/physics/utils/tropopause_read_file.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/physics/utils/tropopause_read_file.F90 b/src/physics/utils/tropopause_read_file.F90 index 4ec967f7..a05677f6 100644 --- a/src/physics/utils/tropopause_read_file.F90 +++ b/src/physics/utils/tropopause_read_file.F90 @@ -20,13 +20,13 @@ module tropopause_read_file !> \section arg_table_tropopause_read_file Argument Table !! \htmlinclude tropopause_read_file.html ! days in year for climatological tropopause pressure data - integer, public, parameter :: tropp_slices = 12 + integer, public, parameter :: tropp_slices = 12 ! climatological tropopause pressures (pcols,ntimes) - real(kind_phys), public, pointer :: tropp_p_loc(:,:) + real(kind_phys), public, allocatable :: tropp_p_loc(:,:) ! monthly day-of-year times corresponding to climatological data (12) - integer, public, pointer :: tropp_days(:) + integer, public, allocatable :: tropp_days(:) contains subroutine tropopause_read_climo_file(tropopause_climo_file) From c22676bfbae4642a30b67946cf9b39320a85305c Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 15 Aug 2024 17:32:19 -0400 Subject: [PATCH 09/25] Fix call to mark_as_initialized for calday (not read from initial file) --- src/control/cam_comp.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index 8c465413..49dfaaf4 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -100,7 +100,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & use physics_grid, only: columns_on_task use vert_coord, only: pver use phys_vars_init_check, only: mark_as_initialized - use tropopause_read_file, only: tropopause_read_file + use tropopause_read_file, only: tropopause_read_climo_file ! Arguments character(len=cl), intent(in) :: caseid ! case ID @@ -173,7 +173,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & ! Get current fractional calendar day. Needs to be updated at every timestep. calday = get_curr_calday() - mark_as_initialized('calday') + call mark_as_initialized('calday') ! Read CAM namelists. filein = "atm_in" // trim(inst_suffix) From cda55758af1d9769e5f7eebcd6f4965377296f26 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 16 Aug 2024 13:10:33 -0400 Subject: [PATCH 10/25] Read tropopause_nl within cam_comp --- cime_config/namelist_definition_cam.xml | 18 +++++++- src/control/cam_comp.F90 | 9 +++- src/control/runtime_opts.F90 | 3 ++ src/physics/utils/tropopause_read_file.F90 | 49 ++++++++++++++++++++- src/physics/utils/tropopause_read_file.meta | 6 +++ 5 files changed, 80 insertions(+), 5 deletions(-) diff --git a/cime_config/namelist_definition_cam.xml b/cime_config/namelist_definition_cam.xml index c12d23a1..d49ef015 100644 --- a/cime_config/namelist_definition_cam.xml +++ b/cime_config/namelist_definition_cam.xml @@ -333,6 +333,22 @@ + + + char*256 + + tropo + tropopause_nl + + Full pathname of boundary dataset for tropopause climatology. + + Default: set by build-namelist. + + + ${DIN_LOC_ROOT}/atm/cam/chem/trop_mozart/ub/clim_p_trop.nc + + + @@ -364,7 +380,7 @@ cam_logfile_nl 0,1,2,3 - Include exra checks and logging output. + Include extra checks and logging output. 0: No extra checks or output 1: Extra informational output 2: Level 1 plus extra checks and informational output diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index 49dfaaf4..ea299053 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -166,6 +166,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & call cam_ctrl_set_orbit(eccen, obliqr, lambm0, mvelpp) + call timemgr_init( & dtime, calendar, start_ymd, start_tod, ref_ymd, & ref_tod, stop_ymd, stop_tod, curr_ymd, curr_tod, & @@ -232,8 +233,12 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & end if ! Read tropopause climatology - ! FIXME hplin 8/15/24: needs to get data from tropopause_nl, how to pass from CCPP? - call tropopause_read_climo_file('/fs/cgd/csm/inputdata/atm/cam/chem/trop_mozart/ub/clim_p_trop.nc') + call tropopause_read_climo_file() + + ! Set tropopause find method as initialized + ! FIXME hplin 8/16/24: probably has to move somewhere else or split into other functions + call mark_as_initialized('tropopause_method') + call mark_as_initialized('tropopause_method_secondary') call phys_init() diff --git a/src/control/runtime_opts.F90 b/src/control/runtime_opts.F90 index 12081a81..11743c37 100644 --- a/src/control/runtime_opts.F90 +++ b/src/control/runtime_opts.F90 @@ -44,6 +44,8 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) ! use cam_diagnostics, only: diag_readnl use inic_analytic_utils, only: analytic_ic_readnl + use tropopause_read_file, only: tropopause_readnl + ! use tracers, only: tracers_readnl ! use nudging, only: nudging_readnl @@ -99,6 +101,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) ! call diag_readnl(nlfilename) ! call check_energy_readnl(nlfilename) call analytic_ic_readnl(nlfilename) + call tropopause_readnl(nlfilename) ! call scam_readnl(nlfilename, single_column, scmlat, scmlon) ! call nudging_readnl(nlfilename) diff --git a/src/physics/utils/tropopause_read_file.F90 b/src/physics/utils/tropopause_read_file.F90 index a05677f6..38bd9fc5 100644 --- a/src/physics/utils/tropopause_read_file.F90 +++ b/src/physics/utils/tropopause_read_file.F90 @@ -8,6 +8,7 @@ module tropopause_read_file !------------------------------------------------------------------- use ccpp_kinds, only: kind_phys + use runtime_obj, only: unset_str use cam_logfile, only: iulog use cam_abortutils, only: endrun use spmd_utils, only: masterproc @@ -15,6 +16,7 @@ module tropopause_read_file implicit none private + public :: tropopause_readnl public :: tropopause_read_climo_file !> \section arg_table_tropopause_read_file Argument Table @@ -28,8 +30,51 @@ module tropopause_read_file ! monthly day-of-year times corresponding to climatological data (12) integer, public, allocatable :: tropp_days(:) + ! Private module data + character(len=256) :: tropopause_climo_file = unset_str + contains - subroutine tropopause_read_climo_file(tropopause_climo_file) + ! Read namelist variable tropopause_climo_file. + ! Containing this within CAM as otherwise the climo file has to be passed from physics -> here + ! then the data from here -> physics. Keeping it at CAM level for now. (hplin, 8/16/24) + subroutine tropopause_readnl(nlfile) + use shr_nl_mod, only: find_group_name => shr_nl_find_group_name + use mpi, only: mpi_character + use spmd_utils, only: mpicom + + ! filepath for file containing namelist input + character(len=*), intent(in) :: nlfile + + ! Local variables + integer :: unitn, ierr + character(len=*), parameter :: subname = 'tropopause_readnl' + + namelist /tropopause_nl/ tropopause_climo_file + + if (masterproc) then + open(newunit=unitn, file=trim(nlfile), status='old') + call find_group_name(unitn, 'tropopause_nl', status=ierr) + if (ierr == 0) then + read(unitn, tropopause_nl, iostat=ierr) + if (ierr /= 0) then + call endrun(subname // ':: ERROR reading namelist') + end if + end if + close(unitn) + end if + + ! Broadcast namelist variables + call mpi_bcast(tropopause_climo_file, len(tropopause_climo_file), mpi_character, 0, mpicom, ierr) + + ! Print out namelist variables + if (masterproc) then + write(iulog,*) subname, ' options:' + write(iulog,*) ' Tropopause climatology file will be read: ', trim(tropopause_climo_file) + endif + + end subroutine tropopause_readnl + + subroutine tropopause_read_climo_file() !------------------------------------------------------------------ ! ... initialize upper boundary values !------------------------------------------------------------------ @@ -43,7 +88,7 @@ subroutine tropopause_read_climo_file(tropopause_climo_file) use pio, only : file_desc_t, var_desc_t, pio_inq_dimid, pio_inq_dimlen, & pio_inq_varid, pio_get_var, pio_closefile, pio_nowrite - character(len=256), intent(in) :: tropopause_climo_file ! absolute path of climatology file + ! character(len=256), intent(in) :: tropopause_climo_file ! absolute path of climatology file !------------------------------------------------------------------ ! ... local variables diff --git a/src/physics/utils/tropopause_read_file.meta b/src/physics/utils/tropopause_read_file.meta index 9545f1ca..1f799672 100644 --- a/src/physics/utils/tropopause_read_file.meta +++ b/src/physics/utils/tropopause_read_file.meta @@ -21,3 +21,9 @@ units = 1 type = integer dimensions = (number_of_months_in_year) +[ tropopause_climo_file ] + standard_name = filename_of_tropopause_climatology + long_name = File path to tropopause climatology file + units = none + type = character | kind = len=256 + dimensions = () \ No newline at end of file From f1d54097a86a0da02ea6efaedf08f22c76e74937 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 16 Aug 2024 13:42:24 -0400 Subject: [PATCH 11/25] Registry updates for tropopause_find --- src/control/cam_comp.F90 | 13 ++++++++++--- src/data/registry.xml | 12 ++++++------ 2 files changed, 16 insertions(+), 9 deletions(-) diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index ea299053..5bf261ba 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -174,7 +174,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & ! Get current fractional calendar day. Needs to be updated at every timestep. calday = get_curr_calday() - call mark_as_initialized('calday') + call mark_as_initialized('fractional_calendar_days_on_end_of_current_timestep') ! Read CAM namelists. filein = "atm_in" // trim(inst_suffix) @@ -237,8 +237,15 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & ! Set tropopause find method as initialized ! FIXME hplin 8/16/24: probably has to move somewhere else or split into other functions - call mark_as_initialized('tropopause_method') - call mark_as_initialized('tropopause_method_secondary') + call mark_as_initialized('control_for_tropopause_find_method') + call mark_as_initialized('control_for_tropopause_find_method_secondary') + + ! FIXME hplin 8/16/24: These are new state variables in CAM-SIMA and updated at every timestep + ! should not be read from ncdata for now. + ! call mark_as_initialized('model_level_number_at_tropopause') + ! call mark_as_initialized('tropopause_air_pressure') + ! call mark_as_initialized('tropopause_air_temperature') + ! call mark_as_initialized('tropopause_altitude') call phys_init() diff --git a/src/data/registry.xml b/src/data/registry.xml index 2e532be2..d9567279 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -348,34 +348,34 @@ 3 - + allocatable="allocatable"> Tropopause model level horizontal_dimension + allocatable="allocatable"> Tropopause level pressure horizontal_dimension + allocatable="allocatable"> Tropopause level temperature horizontal_dimension + allocatable="allocatable"> Tropopause level altitude horizontal_dimension - + --> Date: Mon, 19 Aug 2024 12:12:54 -0400 Subject: [PATCH 12/25] Update for calculating methods needed by CAM and History output (registry only) --- src/data/registry.xml | 43 ------------------------------------------- 1 file changed, 43 deletions(-) diff --git a/src/data/registry.xml b/src/data/registry.xml index d9567279..ef37313e 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -333,49 +333,6 @@ units="1" type="real" kind="kind_phys"> fractional calendar day at end of current timestep - - - Tropopause primary find method - 5 - - - - Tropopause secondary find method - 3 - - - Date: Mon, 19 Aug 2024 13:48:25 -0400 Subject: [PATCH 13/25] Update variable initialization --- src/control/cam_comp.F90 | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index 5bf261ba..8017af76 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -235,17 +235,10 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & ! Read tropopause climatology call tropopause_read_climo_file() - ! Set tropopause find method as initialized - ! FIXME hplin 8/16/24: probably has to move somewhere else or split into other functions - call mark_as_initialized('control_for_tropopause_find_method') - call mark_as_initialized('control_for_tropopause_find_method_secondary') - - ! FIXME hplin 8/16/24: These are new state variables in CAM-SIMA and updated at every timestep - ! should not be read from ncdata for now. - ! call mark_as_initialized('model_level_number_at_tropopause') - ! call mark_as_initialized('tropopause_air_pressure') - ! call mark_as_initialized('tropopause_air_temperature') - ! call mark_as_initialized('tropopause_altitude') + ! Mark variables as initialized so they are not read from initial conditions + call mark_as_initialized('tropopause_air_pressure_from_climatology') + call mark_as_initialized('tropopause_calendar_days_from_climatology') + call mark_as_initialized('filename_of_tropopause_climatology') call phys_init() From 8c1523494ba16d604f22b00aba03441ad94051e5 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 20 Aug 2024 18:21:47 -0400 Subject: [PATCH 14/25] Fix tropp_days is fractional day-of-year --- src/physics/utils/tropopause_read_file.F90 | 2 +- src/physics/utils/tropopause_read_file.meta | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/physics/utils/tropopause_read_file.F90 b/src/physics/utils/tropopause_read_file.F90 index 38bd9fc5..efcad3e8 100644 --- a/src/physics/utils/tropopause_read_file.F90 +++ b/src/physics/utils/tropopause_read_file.F90 @@ -28,7 +28,7 @@ module tropopause_read_file real(kind_phys), public, allocatable :: tropp_p_loc(:,:) ! monthly day-of-year times corresponding to climatological data (12) - integer, public, allocatable :: tropp_days(:) + real(kind_phys), public, allocatable :: tropp_days(:) ! Private module data character(len=256) :: tropopause_climo_file = unset_str diff --git a/src/physics/utils/tropopause_read_file.meta b/src/physics/utils/tropopause_read_file.meta index 1f799672..836532f9 100644 --- a/src/physics/utils/tropopause_read_file.meta +++ b/src/physics/utils/tropopause_read_file.meta @@ -19,7 +19,7 @@ standard_name = tropopause_calendar_days_from_climatology long_name = Climatological tropopause calendar day indices from file units = 1 - type = integer + type = real | kind = kind_phys dimensions = (number_of_months_in_year) [ tropopause_climo_file ] standard_name = filename_of_tropopause_climatology From d727997c743750d0eb0340b123171e27651202ce Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 23 Aug 2024 11:30:24 -0400 Subject: [PATCH 15/25] Update git submodule to jimmielin/hplin/tropopause_find for inclusion of CCPP-ized tropopause_find --- .gitmodules | 4 ++-- src/physics/ncar_ccpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.gitmodules b/.gitmodules index 787b14c9..9d09b7b8 100644 --- a/.gitmodules +++ b/.gitmodules @@ -13,8 +13,8 @@ fxDONOTUSEurl = https://github.com/MPAS-Dev/MPAS-Model.git [submodule "ncar-physics"] path = src/physics/ncar_ccpp - url = https://github.com/ESCOMP/atmospheric_physics - fxtag = atmos_phys0_03_000 + url = https://github.com/jimmielin/atmospheric_physics + fxtag = 2b1d98ef fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/src/physics/ncar_ccpp b/src/physics/ncar_ccpp index f4c09618..2b1d98ef 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit f4c09618eaaa19eaf3382f0473a531e20aa9f808 +Subproject commit 2b1d98ef54403e41310f9f49d82f4ef04d8bbe24 From 2777faee1545e452c7fcb0af6640d139be5477b9 Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Fri, 6 Sep 2024 15:27:52 -0600 Subject: [PATCH 16/25] Add diagnostic/history fillvalue to SIMA registry. --- src/data/registry.xml | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/data/registry.xml b/src/data/registry.xml index ef37313e..aabc6bee 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -346,6 +346,14 @@ ccpp error message '' + + + Fill value for diagnostic/history outputs + 1.e36 + Date: Mon, 9 Sep 2024 11:37:37 -0400 Subject: [PATCH 17/25] Add check for valid numeric string in initial_value in generate_registry_data.py --- src/data/generate_registry_data.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/data/generate_registry_data.py b/src/data/generate_registry_data.py index 33adc954..3ae56567 100755 --- a/src/data/generate_registry_data.py +++ b/src/data/generate_registry_data.py @@ -26,6 +26,7 @@ sys.path.append(__SPINSCRIPTS) # end if _ALL_STRINGS_REGEX = re.compile(r'[A-Za-z][A-Za-z_0-9]+') +_FORTRAN_NUMERIC_REGEX = r'^[+-]?(\d+\.?\d*|\.\d+)([eE][+-]?\d+)?$' # CCPP framework imports # pylint: disable=wrong-import-position @@ -940,9 +941,9 @@ def add_variable(self, newvar): all_strings = _ALL_STRINGS_REGEX.findall(newvar.initial_value) init_val_vars = set() excluded_initializations = {'null', 'nan', 'false', 'true'} - # Exclude NULL and nan variables + # Exclude NULL and nan variables and valid Fortran numeric values that pass the string regex (e.g. 1.e36, -3.2e5) for var in all_strings: - if var.lower() not in excluded_initializations: + if var.lower() not in excluded_initializations and not bool(re.match(_FORTRAN_NUMERIC_REGEX, newvar.initial_value)): init_val_vars.add(var) # end if # end if From 4944a98ad7b0a8c6fc38b4365342d930217ee839 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 9 Sep 2024 13:57:30 -0400 Subject: [PATCH 18/25] Update precision of fillvalue to kind_phys; allow this in the numeric regex match This resolves bit-for-bit issues with fillvalue being written as 9.99e+35 without this fix. --- src/data/generate_registry_data.py | 2 +- src/data/registry.xml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/data/generate_registry_data.py b/src/data/generate_registry_data.py index 3ae56567..c2ece2e2 100755 --- a/src/data/generate_registry_data.py +++ b/src/data/generate_registry_data.py @@ -26,7 +26,7 @@ sys.path.append(__SPINSCRIPTS) # end if _ALL_STRINGS_REGEX = re.compile(r'[A-Za-z][A-Za-z_0-9]+') -_FORTRAN_NUMERIC_REGEX = r'^[+-]?(\d+\.?\d*|\.\d+)([eE][+-]?\d+)?$' +_FORTRAN_NUMERIC_REGEX = r'^[+-]?(\d+\.?\d*|\.\d+)([eE][+-]?\d+)?(_kind_phys)?$' # CCPP framework imports # pylint: disable=wrong-import-position diff --git a/src/data/registry.xml b/src/data/registry.xml index aabc6bee..4ce6fa88 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -352,7 +352,7 @@ units="1" type="real" kind="kind_phys" allocatable="parameter"> Fill value for diagnostic/history outputs - 1.e36 + 1.e36_kind_phys Date: Mon, 9 Sep 2024 16:14:20 -0400 Subject: [PATCH 19/25] Update atmospheric_physics submodule to latest version from PR --- .gitmodules | 2 +- src/physics/ncar_ccpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index 9d09b7b8..a68a275f 100644 --- a/.gitmodules +++ b/.gitmodules @@ -14,7 +14,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/jimmielin/atmospheric_physics - fxtag = 2b1d98ef + fxtag = 5d4a18d6 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/src/physics/ncar_ccpp b/src/physics/ncar_ccpp index 2b1d98ef..5d4a18d6 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit 2b1d98ef54403e41310f9f49d82f4ef04d8bbe24 +Subproject commit 5d4a18d658b494f8d473fd3bc1780c16b4b3942a From 26319f1699b00d1d374be524b48d77cf574d58d2 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 18 Sep 2024 19:08:57 -0400 Subject: [PATCH 20/25] Address review comments; update standard name to tropopause_air_pressure_from_climatology_dataset for tropp_p_loc --- src/control/cam_comp.F90 | 2 +- src/physics/utils/tropopause_read_file.F90 | 4 +--- src/physics/utils/tropopause_read_file.meta | 2 +- 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index 8017af76..2b4b9ae2 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -236,7 +236,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & call tropopause_read_climo_file() ! Mark variables as initialized so they are not read from initial conditions - call mark_as_initialized('tropopause_air_pressure_from_climatology') + call mark_as_initialized('tropopause_air_pressure_from_climatology_dataset') call mark_as_initialized('tropopause_calendar_days_from_climatology') call mark_as_initialized('filename_of_tropopause_climatology') diff --git a/src/physics/utils/tropopause_read_file.F90 b/src/physics/utils/tropopause_read_file.F90 index efcad3e8..be9e3ca6 100644 --- a/src/physics/utils/tropopause_read_file.F90 +++ b/src/physics/utils/tropopause_read_file.F90 @@ -24,7 +24,7 @@ module tropopause_read_file ! days in year for climatological tropopause pressure data integer, public, parameter :: tropp_slices = 12 - ! climatological tropopause pressures (pcols,ntimes) + ! climatological tropopause pressures (ncol,ntimes) real(kind_phys), public, allocatable :: tropp_p_loc(:,:) ! monthly day-of-year times corresponding to climatological data (12) @@ -178,8 +178,6 @@ subroutine tropopause_read_climo_file() ! ... regrid !-------------------------------------------------------------------- - ! tropp_p_loc is allocated with dimensions (pcols, begchunk:endchunk, ntimes) in CAM. - ! in CAM-SIMA, the chunk dimension is collapsed as it is unused. allocate( tropp_p_loc(pcols,ntimes), stat=ierr ) if( ierr /= 0 ) then diff --git a/src/physics/utils/tropopause_read_file.meta b/src/physics/utils/tropopause_read_file.meta index 836532f9..ecb117f8 100644 --- a/src/physics/utils/tropopause_read_file.meta +++ b/src/physics/utils/tropopause_read_file.meta @@ -11,7 +11,7 @@ type = integer dimensions = () [ tropp_p_loc ] - standard_name = tropopause_air_pressure_from_climatology + standard_name = tropopause_air_pressure_from_climatology_dataset units = Pa type = real | kind = kind_phys dimensions = (horizontal_dimension, number_of_months_in_year) From c74a58a2d886cc4309a5671dbad3ff740f430c48 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 26 Sep 2024 11:45:49 -0400 Subject: [PATCH 21/25] Update atmospheric_physics external --- .gitmodules | 2 +- src/physics/ncar_ccpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index a68a275f..41bcfa7c 100644 --- a/.gitmodules +++ b/.gitmodules @@ -14,7 +14,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/jimmielin/atmospheric_physics - fxtag = 5d4a18d6 + fxtag = 886c8952 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/src/physics/ncar_ccpp b/src/physics/ncar_ccpp index 5d4a18d6..886c8952 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit 5d4a18d658b494f8d473fd3bc1780c16b4b3942a +Subproject commit 886c8952eb296c4eea3107f2498a242d1321c490 From 214c9188bd149c42c125d2f0baceee3d66117e20 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Tue, 1 Oct 2024 14:45:11 -0400 Subject: [PATCH 22/25] Update atmospheric_physics external to atmos_phys0_05_000 --- .gitmodules | 4 ++-- src/physics/ncar_ccpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.gitmodules b/.gitmodules index 41bcfa7c..1fd0df45 100644 --- a/.gitmodules +++ b/.gitmodules @@ -13,8 +13,8 @@ fxDONOTUSEurl = https://github.com/MPAS-Dev/MPAS-Model.git [submodule "ncar-physics"] path = src/physics/ncar_ccpp - url = https://github.com/jimmielin/atmospheric_physics - fxtag = 886c8952 + url = https://github.com/ESCOMP/atmospheric_physics + fxtag = atmos_phys0_05_000 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/src/physics/ncar_ccpp b/src/physics/ncar_ccpp index 886c8952..93a1dbf9 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit 886c8952eb296c4eea3107f2498a242d1321c490 +Subproject commit 93a1dbf9c47ccedb8d8a48eba640e48ab2048774 From f85cdf870032a8117179256bd39dd5f7e396e956 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 2 Oct 2024 22:45:33 -0400 Subject: [PATCH 23/25] Update to address review comments: - rename tropopause_read_file -> tropopause_climo_read and corrresponding subroutine names - move ccpp variable initialization calls into read subroutine - improved error handling with errmsg - compile Fortran numeric regex in generate_registry_data.py --- cime_config/namelist_definition_cam.xml | 4 +- src/control/cam_comp.F90 | 9 +- src/control/runtime_opts.F90 | 4 +- src/data/generate_registry_data.py | 4 +- src/data/registry.xml | 2 +- ...ead_file.F90 => tropopause_climo_read.F90} | 126 ++++++++++-------- ...d_file.meta => tropopause_climo_read.meta} | 4 +- 7 files changed, 83 insertions(+), 70 deletions(-) rename src/physics/utils/{tropopause_read_file.F90 => tropopause_climo_read.F90} (59%) rename src/physics/utils/{tropopause_read_file.meta => tropopause_climo_read.meta} (92%) diff --git a/cime_config/namelist_definition_cam.xml b/cime_config/namelist_definition_cam.xml index d49ef015..4181c5d2 100644 --- a/cime_config/namelist_definition_cam.xml +++ b/cime_config/namelist_definition_cam.xml @@ -337,12 +337,10 @@ char*256 - tropo + tropopause tropopause_nl Full pathname of boundary dataset for tropopause climatology. - - Default: set by build-namelist. ${DIN_LOC_ROOT}/atm/cam/chem/trop_mozart/ub/clim_p_trop.nc diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index 2b4b9ae2..e8304af8 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -100,7 +100,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & use physics_grid, only: columns_on_task use vert_coord, only: pver use phys_vars_init_check, only: mark_as_initialized - use tropopause_read_file, only: tropopause_read_climo_file + use tropopause_climo_read, only: tropopause_climo_read_file ! Arguments character(len=cl), intent(in) :: caseid ! case ID @@ -233,12 +233,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url, & end if ! Read tropopause climatology - call tropopause_read_climo_file() - - ! Mark variables as initialized so they are not read from initial conditions - call mark_as_initialized('tropopause_air_pressure_from_climatology_dataset') - call mark_as_initialized('tropopause_calendar_days_from_climatology') - call mark_as_initialized('filename_of_tropopause_climatology') + call tropopause_climo_read_file() call phys_init() diff --git a/src/control/runtime_opts.F90 b/src/control/runtime_opts.F90 index 11743c37..7031be1f 100644 --- a/src/control/runtime_opts.F90 +++ b/src/control/runtime_opts.F90 @@ -44,7 +44,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) ! use cam_diagnostics, only: diag_readnl use inic_analytic_utils, only: analytic_ic_readnl - use tropopause_read_file, only: tropopause_readnl + use tropopause_climo_read, only: tropopause_climo_readnl ! use tracers, only: tracers_readnl ! use nudging, only: nudging_readnl @@ -101,7 +101,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) ! call diag_readnl(nlfilename) ! call check_energy_readnl(nlfilename) call analytic_ic_readnl(nlfilename) - call tropopause_readnl(nlfilename) + call tropopause_climo_readnl(nlfilename) ! call scam_readnl(nlfilename, single_column, scmlat, scmlon) ! call nudging_readnl(nlfilename) diff --git a/src/data/generate_registry_data.py b/src/data/generate_registry_data.py index c2ece2e2..187cf7be 100755 --- a/src/data/generate_registry_data.py +++ b/src/data/generate_registry_data.py @@ -26,7 +26,7 @@ sys.path.append(__SPINSCRIPTS) # end if _ALL_STRINGS_REGEX = re.compile(r'[A-Za-z][A-Za-z_0-9]+') -_FORTRAN_NUMERIC_REGEX = r'^[+-]?(\d+\.?\d*|\.\d+)([eE][+-]?\d+)?(_kind_phys)?$' +_FORTRAN_NUMERIC_REGEX = re.compile(r'^[+-]?(\d+\.?\d*|\.\d+)([eE][+-]?\d+)?(_kind_phys)?$') # CCPP framework imports # pylint: disable=wrong-import-position @@ -943,7 +943,7 @@ def add_variable(self, newvar): excluded_initializations = {'null', 'nan', 'false', 'true'} # Exclude NULL and nan variables and valid Fortran numeric values that pass the string regex (e.g. 1.e36, -3.2e5) for var in all_strings: - if var.lower() not in excluded_initializations and not bool(re.match(_FORTRAN_NUMERIC_REGEX, newvar.initial_value)): + if var.lower() not in excluded_initializations and not _FORTRAN_NUMERIC_REGEX.match(newvar.initial_value): init_val_vars.add(var) # end if # end if diff --git a/src/data/registry.xml b/src/data/registry.xml index 4ce6fa88..0a971e30 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -10,7 +10,7 @@ $SRCROOT/src/data/physconst.meta $SRCROOT/src/physics/utils/physics_grid.meta $SRCROOT/src/physics/utils/cam_constituents.meta - $SRCROOT/src/physics/utils/tropopause_read_file.meta + $SRCROOT/src/physics/utils/tropopause_climo_read.meta $SRCROOT/src/data/air_composition.meta $SRCROOT/src/data/cam_thermo.meta $SRCROOT/src/data/ref_pres.meta diff --git a/src/physics/utils/tropopause_read_file.F90 b/src/physics/utils/tropopause_climo_read.F90 similarity index 59% rename from src/physics/utils/tropopause_read_file.F90 rename to src/physics/utils/tropopause_climo_read.F90 index be9e3ca6..65f90ea6 100644 --- a/src/physics/utils/tropopause_read_file.F90 +++ b/src/physics/utils/tropopause_climo_read.F90 @@ -1,4 +1,4 @@ -module tropopause_read_file +module tropopause_climo_read !------------------------------------------------------------------- ! Support module for CAM-SIMA tropopause_find to read in ! climatological tropopause data. @@ -9,19 +9,16 @@ module tropopause_read_file use ccpp_kinds, only: kind_phys use runtime_obj, only: unset_str - use cam_logfile, only: iulog - use cam_abortutils, only: endrun - use spmd_utils, only: masterproc implicit none private - public :: tropopause_readnl - public :: tropopause_read_climo_file + public :: tropopause_climo_readnl + public :: tropopause_climo_read_file -!> \section arg_table_tropopause_read_file Argument Table -!! \htmlinclude tropopause_read_file.html - ! days in year for climatological tropopause pressure data +!> \section arg_table_tropopause_climo_read Argument Table +!! \htmlinclude tropopause_climo_read.html + ! months in year for climatological tropopause pressure data integer, public, parameter :: tropp_slices = 12 ! climatological tropopause pressures (ncol,ntimes) @@ -35,29 +32,36 @@ module tropopause_read_file contains ! Read namelist variable tropopause_climo_file. - ! Containing this within CAM as otherwise the climo file has to be passed from physics -> here - ! then the data from here -> physics. Keeping it at CAM level for now. (hplin, 8/16/24) - subroutine tropopause_readnl(nlfile) + ! Containing this within CAM-SIMA instead of within scheme as otherwise the climo filepath + ! has to be passed from physics -> here then the data from here -> physics. + subroutine tropopause_climo_readnl(nlfile) use shr_nl_mod, only: find_group_name => shr_nl_find_group_name + use shr_kind_mod, only: shr_kind_cm use mpi, only: mpi_character use spmd_utils, only: mpicom + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + use spmd_utils, only: masterproc ! filepath for file containing namelist input character(len=*), intent(in) :: nlfile ! Local variables - integer :: unitn, ierr - character(len=*), parameter :: subname = 'tropopause_readnl' + integer :: unitn, ierr + character(len=*), parameter :: subname = 'tropopause_climo_readnl' + character(len=shr_kind_cm) :: errmsg namelist /tropopause_nl/ tropopause_climo_file + errmsg = '' + if (masterproc) then open(newunit=unitn, file=trim(nlfile), status='old') call find_group_name(unitn, 'tropopause_nl', status=ierr) if (ierr == 0) then - read(unitn, tropopause_nl, iostat=ierr) + read(unitn, tropopause_nl, iostat=ierr, iomsg=errmsg) if (ierr /= 0) then - call endrun(subname // ':: ERROR reading namelist') + call endrun(subname // ':: ERROR reading namelist:' // errmsg) end if end if close(unitn) @@ -72,23 +76,26 @@ subroutine tropopause_readnl(nlfile) write(iulog,*) ' Tropopause climatology file will be read: ', trim(tropopause_climo_file) endif - end subroutine tropopause_readnl + end subroutine tropopause_climo_readnl - subroutine tropopause_read_climo_file() + subroutine tropopause_climo_read_file() !------------------------------------------------------------------ - ! ... initialize upper boundary values + ! ... read tropopause climatology dataset file !------------------------------------------------------------------ - use interpolate_data, only : lininterp_init, lininterp, interp_type, lininterp_finish - use physics_grid, only : get_rlat_all_p, get_rlon_all_p - use physics_grid, only : pcols => columns_on_task - use physconst, only : pi - use time_manager, only : get_calday - use ioFileMod, only : cam_get_file - use cam_pio_utils, only: cam_pio_openfile - use pio, only : file_desc_t, var_desc_t, pio_inq_dimid, pio_inq_dimlen, & - pio_inq_varid, pio_get_var, pio_closefile, pio_nowrite - - ! character(len=256), intent(in) :: tropopause_climo_file ! absolute path of climatology file + use shr_kind_mod, only: shr_kind_cm + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + use spmd_utils, only: masterproc + use interpolate_data, only: lininterp_init, lininterp, interp_type, lininterp_finish + use physics_grid, only: get_rlat_all_p, get_rlon_all_p + use physics_grid, only: pcols => columns_on_task + use physconst, only: pi + use time_manager, only: get_calday + use ioFileMod, only: cam_get_file + use cam_pio_utils, only: cam_pio_openfile + use pio, only: file_desc_t, var_desc_t, pio_inq_dimid, pio_inq_dimlen + use pio, only: pio_inq_varid, pio_get_var, pio_closefile, pio_nowrite + use phys_vars_init_check, only: mark_as_initialized !------------------------------------------------------------------ ! ... local variables @@ -101,7 +108,10 @@ subroutine tropopause_read_climo_file() integer :: nlon, nlat, ntimes integer :: start(3) integer :: count(3) - integer, parameter :: dates(12) = (/ 116, 214, 316, 415, 516, 615, & + + ! YMD (01-16, 02-14, ...) corresponding to the time slices of the climatological + ! tropopause dataset. Will be converted to day-of-year and stored in tropp_days. + integer, parameter :: dates(tropp_slices) = (/ 116, 214, 316, 415, 516, 615, & 716, 816, 915, 1016, 1115, 1216 /) type(interp_type) :: lon_wgts, lat_wgts real(kind_phys), allocatable :: tropp_p_in(:,:,:) @@ -110,6 +120,9 @@ subroutine tropopause_read_climo_file() real(kind_phys) :: to_lats(pcols), to_lons(pcols) real(kind_phys), parameter :: d2r=pi/180._kind_phys, zero=0._kind_phys, twopi=pi*2._kind_phys character(len=256) :: locfn + character(len=shr_kind_cm) :: errmsg + + errmsg = '' !----------------------------------------------------------------------- ! ... open netcdf file @@ -122,19 +135,19 @@ subroutine tropopause_read_climo_file() !----------------------------------------------------------------------- ierr = pio_inq_dimid( pio_id, 'time', dimid ) ierr = pio_inq_dimlen( pio_id, dimid, ntimes ) - if( ntimes /= 12 )then - write(iulog,*) 'tropopause_read_climo_file: number of months = ',ntimes,'; expecting 12' - call endrun('tropopause_read_climo_file: number of months in file incorrect') + if( ntimes /= tropp_slices ) then + write(iulog,*) 'tropopause_climo_read_file: number of months = ',ntimes,'; expecting ',tropp_slices + call endrun('tropopause_climo_read_file: number of months in file incorrect') end if !----------------------------------------------------------------------- ! ... get latitudes !----------------------------------------------------------------------- ierr = pio_inq_dimid( pio_id, 'lat', dimid ) ierr = pio_inq_dimlen( pio_id, dimid, nlat ) - allocate( lat(nlat), stat=ierr ) + allocate( lat(nlat), stat=ierr, errmsg=errmsg ) if( ierr /= 0 ) then - write(iulog,*) 'tropopause_read_climo_file: lat allocation error = ',ierr - call endrun('tropopause_read_climo_file: failed to allocate lat') + write(iulog,*) 'tropopause_climo_read_file: lat allocation error = ',ierr + call endrun('tropopause_climo_read_file: failed to allocate lat, error = ' // errmsg) end if ierr = pio_inq_varid( pio_id, 'lat', vid ) ierr = pio_get_var( pio_id, vid, lat ) @@ -144,10 +157,10 @@ subroutine tropopause_read_climo_file() !----------------------------------------------------------------------- ierr = pio_inq_dimid( pio_id, 'lon', dimid ) ierr = pio_inq_dimlen( pio_id, dimid, nlon ) - allocate( lon(nlon), stat=ierr ) + allocate( lon(nlon), stat=ierr, errmsg=errmsg ) if( ierr /= 0 ) then - write(iulog,*) 'tropopause_read_climo_file: lon allocation error = ',ierr - call endrun('tropopause_read_climo_file: failed to allocate lon') + write(iulog,*) 'tropopause_climo_read_file: lon allocation error = ',ierr + call endrun('tropopause_climo_read_file: failed to allocate lon, error = ' // errmsg) end if ierr = pio_inq_varid( pio_id, 'lon', vid ) ierr = pio_get_var( pio_id, vid, lon ) @@ -156,10 +169,10 @@ subroutine tropopause_read_climo_file() !------------------------------------------------------------------ ! ... allocate arrays !------------------------------------------------------------------ - allocate( tropp_p_in(nlon,nlat,ntimes), stat=ierr ) + allocate( tropp_p_in(nlon,nlat,ntimes), stat=ierr, errmsg=errmsg ) if( ierr /= 0 ) then - write(iulog,*) 'tropopause_read_climo_file: tropp_p_in allocation error = ',ierr - call endrun('tropopause_read_climo_file: failed to allocate tropp_p_in') + write(iulog,*) 'tropopause_climo_read_file: tropp_p_in allocation error = ',ierr + call endrun('tropopause_climo_read_file: failed to allocate tropp_p_in, error = ' // errmsg) end if !------------------------------------------------------------------ ! ... read in the tropopause pressure @@ -178,11 +191,11 @@ subroutine tropopause_read_climo_file() ! ... regrid !-------------------------------------------------------------------- - allocate( tropp_p_loc(pcols,ntimes), stat=ierr ) + allocate( tropp_p_loc(pcols,ntimes), stat=ierr, errmsg=errmsg ) if( ierr /= 0 ) then - write(iulog,*) 'tropopause_read_climo_file: tropp_p_loc allocation error = ',ierr - call endrun('tropopause_read_climo_file: failed to allocate tropp_p_loc') + write(iulog,*) 'tropopause_climo_read_file: tropp_p_loc allocation error = ',ierr + call endrun('tropopause_climo_read_file: failed to allocate tropp_p_loc, error = ' // errmsg) end if call get_rlat_all_p(pcols, to_lats) @@ -202,19 +215,26 @@ subroutine tropopause_read_climo_file() ! ... initialize the monthly day of year times !-------------------------------------------------------- - allocate( tropp_days(tropp_slices), stat=ierr ) + allocate( tropp_days(tropp_slices), stat=ierr, errmsg=errmsg ) if( ierr /= 0 ) then - write(iulog,*) 'tropopause_read_climo_file: tropp_days allocation error = ',ierr - call endrun('tropopause_read_climo_file: failed to allocate tropp_days') + write(iulog,*) 'tropopause_climo_read_file: tropp_days allocation error = ',ierr + call endrun('tropopause_climo_read_file: failed to allocate tropp_days, error = ' // errmsg) end if - do n = 1,12 + do n = 1,tropp_slices tropp_days(n) = get_calday( dates(n), 0 ) end do if (masterproc) then - write(iulog,*) 'tropopause_read_climo_file : tropp_days' + write(iulog,*) 'tropopause_climo_read_file: tropp_days (fractional day-of-year in climatology dataset) =' write(iulog,'(1p,5g15.8)') tropp_days(:) endif - end subroutine tropopause_read_climo_file -end module tropopause_read_file \ No newline at end of file + !-------------------------------------------------------- + ! Mark variables as initialized so they are not read from initial conditions + !-------------------------------------------------------- + call mark_as_initialized('tropopause_air_pressure_from_climatology_dataset') + call mark_as_initialized('tropopause_calendar_days_from_climatology') + call mark_as_initialized('filename_of_tropopause_climatology') + + end subroutine tropopause_climo_read_file +end module tropopause_climo_read \ No newline at end of file diff --git a/src/physics/utils/tropopause_read_file.meta b/src/physics/utils/tropopause_climo_read.meta similarity index 92% rename from src/physics/utils/tropopause_read_file.meta rename to src/physics/utils/tropopause_climo_read.meta index ecb117f8..440198b4 100644 --- a/src/physics/utils/tropopause_read_file.meta +++ b/src/physics/utils/tropopause_climo_read.meta @@ -1,9 +1,9 @@ [ccpp-table-properties] - name = tropopause_read_file + name = tropopause_climo_read type = module [ccpp-arg-table] - name = tropopause_read_file + name = tropopause_climo_read type = module [ tropp_slices ] standard_name = number_of_months_in_year From d226e34a949373c9d14e7cfbbcb004edda519324 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 3 Oct 2024 16:58:13 -0400 Subject: [PATCH 24/25] Remove commented out code in namelist_definition_cam.xml --- cime_config/namelist_definition_cam.xml | 1 - 1 file changed, 1 deletion(-) diff --git a/cime_config/namelist_definition_cam.xml b/cime_config/namelist_definition_cam.xml index 4181c5d2..ba8fbe40 100644 --- a/cime_config/namelist_definition_cam.xml +++ b/cime_config/namelist_definition_cam.xml @@ -336,7 +336,6 @@ char*256 - tropopause tropopause_nl From d8f13a5641022ebb9bc94ad33d8c790ffd729aea Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 7 Oct 2024 11:20:48 -0600 Subject: [PATCH 25/25] Change string length --- src/physics/utils/tropopause_climo_read.F90 | 7 ++++--- src/physics/utils/tropopause_climo_read.meta | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/physics/utils/tropopause_climo_read.F90 b/src/physics/utils/tropopause_climo_read.F90 index 65f90ea6..e7dea239 100644 --- a/src/physics/utils/tropopause_climo_read.F90 +++ b/src/physics/utils/tropopause_climo_read.F90 @@ -9,6 +9,7 @@ module tropopause_climo_read use ccpp_kinds, only: kind_phys use runtime_obj, only: unset_str + use shr_kind_mod, only: shr_kind_cl implicit none private @@ -28,7 +29,7 @@ module tropopause_climo_read real(kind_phys), public, allocatable :: tropp_days(:) ! Private module data - character(len=256) :: tropopause_climo_file = unset_str + character(len=shr_kind_cl) :: tropopause_climo_file = unset_str contains ! Read namelist variable tropopause_climo_file. @@ -119,7 +120,7 @@ subroutine tropopause_climo_read_file() real(kind_phys), allocatable :: lon(:) real(kind_phys) :: to_lats(pcols), to_lons(pcols) real(kind_phys), parameter :: d2r=pi/180._kind_phys, zero=0._kind_phys, twopi=pi*2._kind_phys - character(len=256) :: locfn + character(len=shr_kind_cl) :: locfn character(len=shr_kind_cm) :: errmsg errmsg = '' @@ -237,4 +238,4 @@ subroutine tropopause_climo_read_file() call mark_as_initialized('filename_of_tropopause_climatology') end subroutine tropopause_climo_read_file -end module tropopause_climo_read \ No newline at end of file +end module tropopause_climo_read diff --git a/src/physics/utils/tropopause_climo_read.meta b/src/physics/utils/tropopause_climo_read.meta index 440198b4..65dfbcf6 100644 --- a/src/physics/utils/tropopause_climo_read.meta +++ b/src/physics/utils/tropopause_climo_read.meta @@ -26,4 +26,4 @@ long_name = File path to tropopause climatology file units = none type = character | kind = len=256 - dimensions = () \ No newline at end of file + dimensions = ()