From 2e063bcd3eec4b96b735be25dcb14bfb3dbd47e7 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 14 Aug 2024 19:06:21 -0400 Subject: [PATCH 01/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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 d6f830aa1e90921e00448a82605450b8dfaf52ab Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 6 Sep 2024 11:36:47 -0400 Subject: [PATCH 16/32] Read total energy using dycore formula from fixed CAM snapshots - read ncdata variables te_ini_phys, te_ini_dyn, te_cur_phys, te_cur_dyn - update to standard names per discussion in https://github.com/ESCOMP/CAM/issues/1141 --- src/data/registry.xml | 28 ++++++++++++++++----- tools/stdnames_to_inputnames_dictionary.xml | 24 +++++++++++++----- 2 files changed, 40 insertions(+), 12 deletions(-) diff --git a/src/data/registry.xml b/src/data/registry.xml index 57a9ae3b..0b607f09 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -202,21 +202,37 @@ horizontal_dimension vertical_interface_dimension zi state_zi - horizontal_dimension - te_ini state_te_ini + te_ini_phys state_te_ini_phys - horizontal_dimension - te_cur state_te_cur + te_cur_phys state_te_cur_phys + + + horizontal_dimension + te_ini_dyn state_te_ini_dyn + + + horizontal_dimension + te_cur_dyn state_te_cur_dyn state_zi - + - te_ini - state_te_ini + te_ini_phys + state_te_ini_phys - + - te_cur - state_te_cur + te_cur_phys + state_te_cur_phys + + + + + te_ini_dyn + state_te_ini_dyn + + + + + te_cur_dyn + state_te_cur_dyn From 5e4dafff6812c1ee7687e070a4932a7c9a3c0931 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 6 Sep 2024 11:39:03 -0400 Subject: [PATCH 17/32] Update tw_ini, tw_cur standard names --- src/data/registry.xml | 4 ++-- tools/stdnames_to_inputnames_dictionary.xml | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/data/registry.xml b/src/data/registry.xml index 0b607f09..590ea2a2 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -235,7 +235,7 @@ te_cur_dyn state_te_cur_dyn @@ -243,7 +243,7 @@ tw_ini state_tw_ini diff --git a/tools/stdnames_to_inputnames_dictionary.xml b/tools/stdnames_to_inputnames_dictionary.xml index 0c71e351..77e81cd2 100644 --- a/tools/stdnames_to_inputnames_dictionary.xml +++ b/tools/stdnames_to_inputnames_dictionary.xml @@ -184,13 +184,13 @@ state_te_cur_dyn - + tw_ini state_tw_ini - + tw_cur state_tw_cur From 2777faee1545e452c7fcb0af6640d139be5477b9 Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Fri, 6 Sep 2024 15:27:52 -0600 Subject: [PATCH 18/32] 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: Fri, 6 Sep 2024 15:48:36 -0600 Subject: [PATCH 19/32] add file to track test failures --- test/existing-test-failures.txt | 1 + 1 file changed, 1 insertion(+) create mode 100644 test/existing-test-failures.txt diff --git a/test/existing-test-failures.txt b/test/existing-test-failures.txt new file mode 100644 index 00000000..92c15920 --- /dev/null +++ b/test/existing-test-failures.txt @@ -0,0 +1 @@ +No test failures From a6cede25add298a43a6db5973a493aa70e7bc963 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 9 Sep 2024 11:37:37 -0400 Subject: [PATCH 20/32] 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 21/32] 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 22/32] 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 3112ab620ab36c0755bd1e3dc8d682d1b81ed730 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 16 Sep 2024 18:37:37 -0400 Subject: [PATCH 23/32] Update standard names for tw_ini, tw_cur; add to physics_state ddt --- src/data/registry.xml | 10 ++++++++-- tools/stdnames_to_inputnames_dictionary.xml | 4 ++-- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/data/registry.xml b/src/data/registry.xml index 590ea2a2..9bca0f62 100644 --- a/src/data/registry.xml +++ b/src/data/registry.xml @@ -235,7 +235,7 @@ te_cur_dyn state_te_cur_dyn @@ -243,7 +243,7 @@ tw_ini state_tw_ini @@ -313,6 +313,12 @@ inverse_exner_function_wrt_surface_pressure frontogenesis_function frontogenesis_angle + vertically_integrated_total_energy_of_initial_state_using_physics_energy_formula + vertically_integrated_total_energy_of_current_state_using_physics_energy_formula + vertically_integrated_total_energy_of_initial_state_using_dycore_energy_formula + vertically_integrated_total_energy_of_current_state_using_dycore_energy_formula + vertically_integrated_water_vapor_and_condensed_water_of_initial_state + vertically_integrated_water_vapor_and_condensed_water_of_current_state tendency_of_air_temperature_due_to_model_physics diff --git a/tools/stdnames_to_inputnames_dictionary.xml b/tools/stdnames_to_inputnames_dictionary.xml index 77e81cd2..0ac72a21 100644 --- a/tools/stdnames_to_inputnames_dictionary.xml +++ b/tools/stdnames_to_inputnames_dictionary.xml @@ -184,13 +184,13 @@ state_te_cur_dyn - + tw_ini state_tw_ini - + tw_cur state_tw_cur From 26319f1699b00d1d374be524b48d77cf574d58d2 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Wed, 18 Sep 2024 19:08:57 -0400 Subject: [PATCH 24/32] 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 f8698e08c7a28acb08b5cadb7f378d99bead43f4 Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Tue, 24 Sep 2024 13:44:27 -0600 Subject: [PATCH 25/32] Create LICENSE Add Apache 2.0 license to CAM-SIMA. --- LICENSE | 201 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 201 insertions(+) create mode 100644 LICENSE diff --git a/LICENSE b/LICENSE new file mode 100644 index 00000000..261eeb9e --- /dev/null +++ b/LICENSE @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. From c74a58a2d886cc4309a5671dbad3ff740f430c48 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 26 Sep 2024 11:45:49 -0400 Subject: [PATCH 26/32] 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 27/32] 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 28/32] 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 29/32] 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 30/32] 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 = () From 9f1f6691370456420c2b243653cef431c54cc267 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 7 Oct 2024 21:06:20 -0600 Subject: [PATCH 31/32] Fix compile error with tropopause_find and CCPP fortran parser The CCPP Fortran parser is not able to read module variables of character type defined with the shr_kind_cl parameter size. The tropopause_climo_file variable does not need to be a CCPP variable so it has been removed. This will fix the build process again. --- src/physics/utils/tropopause_climo_read.F90 | 7 +++---- src/physics/utils/tropopause_climo_read.meta | 6 ------ 2 files changed, 3 insertions(+), 10 deletions(-) diff --git a/src/physics/utils/tropopause_climo_read.F90 b/src/physics/utils/tropopause_climo_read.F90 index e7dea239..7033d91b 100644 --- a/src/physics/utils/tropopause_climo_read.F90 +++ b/src/physics/utils/tropopause_climo_read.F90 @@ -17,6 +17,9 @@ module tropopause_climo_read public :: tropopause_climo_readnl public :: tropopause_climo_read_file + ! Private module data + character(len=shr_kind_cl) :: tropopause_climo_file = unset_str + !> \section arg_table_tropopause_climo_read Argument Table !! \htmlinclude tropopause_climo_read.html ! months in year for climatological tropopause pressure data @@ -28,9 +31,6 @@ module tropopause_climo_read ! monthly day-of-year times corresponding to climatological data (12) real(kind_phys), public, allocatable :: tropp_days(:) - ! Private module data - character(len=shr_kind_cl) :: tropopause_climo_file = unset_str - contains ! Read namelist variable tropopause_climo_file. ! Containing this within CAM-SIMA instead of within scheme as otherwise the climo filepath @@ -235,7 +235,6 @@ subroutine tropopause_climo_read_file() !-------------------------------------------------------- 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 diff --git a/src/physics/utils/tropopause_climo_read.meta b/src/physics/utils/tropopause_climo_read.meta index 65dfbcf6..6d4d8538 100644 --- a/src/physics/utils/tropopause_climo_read.meta +++ b/src/physics/utils/tropopause_climo_read.meta @@ -21,9 +21,3 @@ units = 1 type = real | kind = kind_phys 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 = () From 112cc2b8fd3ef067c8f5293b2b2c24e962b2e312 Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Tue, 15 Oct 2024 09:09:20 -0600 Subject: [PATCH 32/32] Implement re-organized CCPP physics external (#306) Originator(s): nusbaume Summary (include the keyword ['closes', 'fixes', 'resolves'] and issue number): This PR bring in a new atmospheric_physics version with a re-organized directory structure. The way this impacts CAM-SIMA is that now when a physics suite is listed it will first look in the case's SourceMods first, `ncar_ccpp/suites` second, and `ncar_ccpp/test/test_suites` third, after which it will error if it hasn't found the Suite Definition File (SDF). The CAM-SIMA build system also now looks in `ncar_ccpp/schemes` for any CCPP physics source code and metadata files. Fixes #305 The associated atmospehric_physics PR can be found here: ESCOMP/atmospheric_physics#126 Describe any changes made to build system: The CAM-SIMA configuration routines will now look under `suites` or `test_suites` for SDFs, with the `suites` entries taking precedence (unless the case has SourceMods, in which those always take precedence). All of the relevant source code and metadata files for atmospheric physics must also now be present under the `schemes` directory in the atmospheric_physics repo. Describe any changes made to the namelist: N/A List any changes to the defaults for the input datasets (e.g. boundary datasets): N/A List all files eliminated and why: Remove unused "test/include" directory: D test/include/Makefile D test/include/cam_abortutils.F90 D test/include/cam_logfile.F90 D test/include/ccpp_kinds.F90 D test/include/shr_infnan_mod.F90 D test/include/shr_kind_mod.F90 D test/include/spmd_utils.F90 List all files added and what they do: N/A List all existing files that have been modified, and describe the changes: (Helpful git command: git diff --name-status development...) Update atmospheric_physics external: M .gitmodules M src/physics/ncar_ccpp Update CCPP SDF, source, and metadata file search locations and precedence: M cime_config/cam_autogen.py Add FTJ16 compset and cleanup simple physics configuration options: M cime_config/config_component.xml If there are new failures (compare to the existing-test-failures.txt file), have them OK'd by the gatekeeper, note them here, and add them to the file. If there are baseline differences, include the test and the reason for the diff. What is the nature of the change? Roundoff? derecho/intel/aux_sima: ALL PASS derecho/gnu/aux_sima: ALL PASS CAM-SIMA date used for the baseline comparison tests if different than latest: --- .gitmodules | 2 +- cime_config/cam_autogen.py | 21 +- cime_config/config_component.xml | 10 +- src/physics/ncar_ccpp | 2 +- test/include/Makefile | 14 - test/include/cam_abortutils.F90 | 17 - test/include/cam_logfile.F90 | 96 -- test/include/ccpp_kinds.F90 | 10 - test/include/shr_infnan_mod.F90 | 1907 ------------------------------ test/include/shr_kind_mod.F90 | 19 - test/include/spmd_utils.F90 | 11 - 11 files changed, 21 insertions(+), 2088 deletions(-) delete mode 100644 test/include/Makefile delete mode 100644 test/include/cam_abortutils.F90 delete mode 100644 test/include/cam_logfile.F90 delete mode 100644 test/include/ccpp_kinds.F90 delete mode 100644 test/include/shr_infnan_mod.F90 delete mode 100644 test/include/shr_kind_mod.F90 delete mode 100644 test/include/spmd_utils.F90 diff --git a/.gitmodules b/.gitmodules index 949b88cd..1fc97042 100644 --- a/.gitmodules +++ b/.gitmodules @@ -14,7 +14,7 @@ [submodule "ncar-physics"] path = src/physics/ncar_ccpp url = https://github.com/ESCOMP/atmospheric_physics - fxtag = atmos_phys0_05_000 + fxtag = atmos_phys0_05_001 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics [submodule "ccs_config"] diff --git a/cime_config/cam_autogen.py b/cime_config/cam_autogen.py index eeb31229..863806ec 100644 --- a/cime_config/cam_autogen.py +++ b/cime_config/cam_autogen.py @@ -435,10 +435,15 @@ def generate_physics_suites(build_cache, preproc_defs, host_name, if not os.path.exists(physics_blddir): os.makedirs(physics_blddir) # End if - # Collect all source directories - atm_phys_src_dir = os.path.join(atm_root, "src", "physics", "ncar_ccpp") - source_search = [source_mods_dir, atm_phys_src_dir] - # Find all metadata files, organize by scheme name + # Set top-level CCPP physics directory + atm_phys_top_dir = os.path.join(atm_root, "src", "physics", "ncar_ccpp") + # Collect all possible Suite Definition File (SDF) locations + atm_suites_path = os.path.join(atm_phys_top_dir, "suites") + atm_test_suites_path = os.path.join(atm_phys_top_dir, "test", "test_suites") + suite_search = [source_mods_dir, atm_suites_path, atm_test_suites_path] + # Find all scheme metadata files, organized by scheme name + atm_schemes_path = os.path.join(atm_phys_top_dir, "schemes") + source_search = [source_mods_dir, atm_schemes_path] all_scheme_files = _find_metadata_files(source_search, find_scheme_names) # Find the SDFs specified for this model build @@ -446,11 +451,15 @@ def generate_physics_suites(build_cache, preproc_defs, host_name, scheme_files = [] xml_files = {} # key is scheme, value is xml file path for sdf in phys_suites_str.split(';'): - sdf_path = _find_file(f"suite_{sdf}.xml", source_search) + sdf_path = _find_file(f"suite_{sdf}.xml", suite_search) if not sdf_path: emsg = f"ERROR: Unable to find SDF for suite '{sdf}'" raise CamAutoGenError(emsg) # End if + if os.path.dirname(os.path.abspath(sdf_path)) == atm_test_suites_path: + #Notify user that a test suite is being used + _LOGGER.info("Using non-standard test suite: %s", sdf) + # End if sdfs.append(sdf_path) # Given an SDF, find all the schemes it calls _, suite = read_xml_file(sdf_path) @@ -587,7 +596,7 @@ def generate_physics_suites(build_cache, preproc_defs, host_name, # there to the bld directory: if do_gen_ccpp: # Set CCPP physics "utilities" path - atm_phys_util_dir = os.path.join(atm_phys_src_dir, "utilities") + atm_phys_util_dir = os.path.join(atm_schemes_path, "utilities") # Check that directory exists if not os.path.isdir(atm_phys_util_dir): diff --git a/cime_config/config_component.xml b/cime_config/config_component.xml index 9b715c97..c93a9de8 100644 --- a/cime_config/config_component.xml +++ b/cime_config/config_component.xml @@ -159,12 +159,10 @@ -nlev 145 --> - + + --physics-suites tj2016 --analytic_ic + --physics-suites kessler --analytic_ic --physics-suites held_suarez_1994 --analytic_ic --dyn none --physics-suites adiabatic diff --git a/src/physics/ncar_ccpp b/src/physics/ncar_ccpp index 93a1dbf9..f8ce60bf 160000 --- a/src/physics/ncar_ccpp +++ b/src/physics/ncar_ccpp @@ -1 +1 @@ -Subproject commit 93a1dbf9c47ccedb8d8a48eba640e48ab2048774 +Subproject commit f8ce60bf40f800623f8eb3065021ec5dfa9e6b45 diff --git a/test/include/Makefile b/test/include/Makefile deleted file mode 100644 index 699930a6..00000000 --- a/test/include/Makefile +++ /dev/null @@ -1,14 +0,0 @@ -FC = gfortran -FFLAGS = -c -DCPRGNU - -SOURCES = shr_kind_mod.F90 shr_infnan_mod.F90 ccpp_kinds.F90 cam_abortutils.F90 -SOURCES += spmd_utils.F90 cam_logfile.F90 -OBJS = $(SOURCES:.F90=.o) - -all: objs - -objs: $(SOURCES) - $(FC) $(FFLAGS) $(SOURCES) - -clean: - ${RM} *.o *.mod diff --git a/test/include/cam_abortutils.F90 b/test/include/cam_abortutils.F90 deleted file mode 100644 index 8db9729e..00000000 --- a/test/include/cam_abortutils.F90 +++ /dev/null @@ -1,17 +0,0 @@ -module cam_abortutils - - implicit none - private - - public endrun - -CONTAINS - - subroutine endrun(msg) - character(len=*), intent(in) :: msg - - write(6, *) msg - STOP - end subroutine endrun - -end module cam_abortutils diff --git a/test/include/cam_logfile.F90 b/test/include/cam_logfile.F90 deleted file mode 100644 index 8e1a8998..00000000 --- a/test/include/cam_logfile.F90 +++ /dev/null @@ -1,96 +0,0 @@ -module cam_logfile - -!----------------------------------------------------------------------- -! -! Purpose: This module is responsible for managing the logical unit -! of CAM's output log -! -! Author: mvr, Sep 2007 -! -!----------------------------------------------------------------------- - -!----------------------------------------------------------------------- -!- use statements ------------------------------------------------------ -!----------------------------------------------------------------------- -!----------------------------------------------------------------------- -!- module boilerplate -------------------------------------------------- -!----------------------------------------------------------------------- - implicit none - private - save - -!----------------------------------------------------------------------- -! Public interfaces ---------------------------------------------------- -!----------------------------------------------------------------------- - public :: cam_set_log_unit - public :: cam_logfile_readnl - public :: cam_log_multiwrite -!----------------------------------------------------------------------- -! Public data ---------------------------------------------------------- -!----------------------------------------------------------------------- - integer, public, protected :: iulog = 6 - integer, public, parameter :: DEBUGOUT_NONE = 0 - integer, public, parameter :: DEBUGOUT_INFO = 1 - integer, public, parameter :: DEBUGOUT_VERBOSE = 2 - integer, public, parameter :: DEBUGOUT_DEBUG = 3 - integer, public, protected :: debug_output = DEBUGOUT_NONE - -!----------------------------------------------------------------------- -! Private data --------------------------------------------------------- -!----------------------------------------------------------------------- - logical :: iulog_set = .true. - - interface cam_log_multiwrite - module procedure cam_log_multiwrite_ni ! Multiple integers - end interface cam_log_multiwrite - -CONTAINS - -!----------------------------------------------------------------------- -! Subroutines and functions -------------------------------------------- -!----------------------------------------------------------------------- - - subroutine cam_set_log_unit(unit_num) - - integer, intent(in) :: unit_num - - ! Change iulog to unit_num on this PE or log a waring - ! The log unit number can be set at most once per run - if (iulog_set) then - write(iulog, *) 'cam_set_log_unit: Cannot change log unit during run' - else - iulog = unit_num - iulog_set = .true. - end if - end subroutine cam_set_log_unit - - subroutine cam_logfile_readnl(nlfile) - - ! nlfile: filepath for file containing namelist input - character(len=*), intent(in) :: nlfile - - end subroutine cam_logfile_readnl - - subroutine cam_log_multiwrite_ni(subname, headers, fmt_string, values) - ! Print out values from every task - use spmd_utils, only: masterproc - - ! Dummy arguments - character(len=*), intent(in) :: subname - character(len=*), intent(in) :: headers - character(len=*), intent(in) :: fmt_string - integer, intent(in) :: values(:) - ! Local variables - integer :: num_fields - integer :: fnum - - num_fields = size(values, 1) - - if (masterproc) then - write(iulog, '(2a)') trim(subname), trim(headers) - write(iulog, fmt_string) subname, 0, & - (values(fnum), fnum = 1, num_fields) - end if - end subroutine cam_log_multiwrite_ni - -end module cam_logfile diff --git a/test/include/ccpp_kinds.F90 b/test/include/ccpp_kinds.F90 deleted file mode 100644 index c90c9cae..00000000 --- a/test/include/ccpp_kinds.F90 +++ /dev/null @@ -1,10 +0,0 @@ -module ccpp_kinds - - use ISO_FORTRAN_ENV, only: kind_phys => REAL64 - - implicit none - private - - public kind_phys - -end module ccpp_kinds diff --git a/test/include/shr_infnan_mod.F90 b/test/include/shr_infnan_mod.F90 deleted file mode 100644 index 8863882d..00000000 --- a/test/include/shr_infnan_mod.F90 +++ /dev/null @@ -1,1907 +0,0 @@ -! This file is a stand-in for CIME's shr_infnan_mod.F90.in -!=================================================== - -! Flag representing compiler support of Fortran 2003's -! ieee_arithmetic intrinsic module. -#if defined CPRIBM || defined CPRPGI || defined CPRINTEL || defined CPRCRAY || defined CPRNAG -#define HAVE_IEEE_ARITHMETIC -#endif - -module shr_infnan_mod -!--------------------------------------------------------------------- -! Module to test for IEEE Inf and NaN values, which also provides a -! method of setting +/-Inf and signaling or quiet NaN. -! -! All functions are elemental, and thus work on arrays. -!--------------------------------------------------------------------- -! To test for these values, just call the corresponding function, e.g: -! -! var_is_nan = shr_infnan_isnan(x) -! -! You can also use it on arrays: -! -! array_contains_nan = any(shr_infnan_isnan(my_array)) -! -!--------------------------------------------------------------------- -! To generate these values, assign one of the provided derived-type -! variables to a real: -! -! use shr_infnan_mod, only: nan => shr_infnan_nan, & -! inf => shr_infnan_inf, & -! assignment(=) -! real(r4) :: my_nan -! real(r8) :: my_inf_array(2,2) -! my_nan = nan -! my_inf_array = inf -! -! Keep in mind that "shr_infnan_nan" and "shr_infnan_inf" cannot be -! passed to functions that expect real arguments. To pass a real -! NaN, you will have to use shr_infnan_nan to set a local real of -! the correct kind. -!--------------------------------------------------------------------- - -use shr_kind_mod, only: & - r4 => SHR_KIND_R4, & - r8 => SHR_KIND_R8 - -#ifdef HAVE_IEEE_ARITHMETIC - -! If we have IEEE_ARITHMETIC, the NaN test is provided for us. -use, intrinsic :: ieee_arithmetic, only: & - shr_infnan_isnan => ieee_is_nan - -#else - -! Integers of correct size for bit patterns below. -use shr_kind_mod, only: i4 => shr_kind_i4, i8 => shr_kind_i8 - -#endif - -implicit none -private -save - -! Test functions for NaN/Inf values. -public :: shr_infnan_isnan -public :: shr_infnan_isinf -public :: shr_infnan_isposinf -public :: shr_infnan_isneginf - -! Locally defined isnan. -#ifndef HAVE_IEEE_ARITHMETIC - -interface shr_infnan_isnan - ! TYPE double,real - module procedure shr_infnan_isnan_double - ! TYPE double,real - module procedure shr_infnan_isnan_real -end interface -#endif - - -interface shr_infnan_isinf - ! TYPE double,real - module procedure shr_infnan_isinf_double - ! TYPE double,real - module procedure shr_infnan_isinf_real -end interface - - -interface shr_infnan_isposinf - ! TYPE double,real - module procedure shr_infnan_isposinf_double - ! TYPE double,real - module procedure shr_infnan_isposinf_real -end interface - - -interface shr_infnan_isneginf - ! TYPE double,real - module procedure shr_infnan_isneginf_double - ! TYPE double,real - module procedure shr_infnan_isneginf_real -end interface - -! Derived types for generation of NaN/Inf -! Even though there's no reason to "use" the types directly, some compilers -! might have trouble with an object being used without its type. -public :: shr_infnan_nan_type -public :: shr_infnan_inf_type -public :: assignment(=) -public :: shr_infnan_to_r4 -public :: shr_infnan_to_r8 - -! Type representing Not A Number. -type :: shr_infnan_nan_type - logical :: quiet = .false. -end type shr_infnan_nan_type - -! Type representing +/-Infinity. -type :: shr_infnan_inf_type - logical :: positive = .true. -end type shr_infnan_inf_type - -! Allow assigning reals to NaN or Inf. - -interface assignment(=) - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_0d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_1d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_2d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_3d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_4d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_5d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_6d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_7d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_0d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_1d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_2d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_3d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_4d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_5d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_6d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_nan_7d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_0d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_1d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_2d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_3d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_4d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_5d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_6d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_7d_double - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_0d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_1d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_2d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_3d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_4d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_5d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_6d_real - ! TYPE double,real - ! DIMS 0,1,2,3,4,5,6,7 - module procedure set_inf_7d_real -end interface - -! Conversion functions. - -interface shr_infnan_to_r8 - module procedure nan_r8 - module procedure inf_r8 -end interface - - -interface shr_infnan_to_r4 - module procedure nan_r4 - module procedure inf_r4 -end interface - -! Initialize objects of NaN/Inf type for other modules to use. - -! Default NaN is signaling, but also provide snan and qnan to choose -! explicitly. -type(shr_infnan_nan_type), public, parameter :: shr_infnan_nan = & - shr_infnan_nan_type(.false.) -type(shr_infnan_nan_type), public, parameter :: shr_infnan_snan = & - shr_infnan_nan_type(.false.) -type(shr_infnan_nan_type), public, parameter :: shr_infnan_qnan = & - shr_infnan_nan_type(.true.) - -! Default Inf is positive, but provide posinf to go with neginf. -type(shr_infnan_inf_type), public, parameter :: shr_infnan_inf = & - shr_infnan_inf_type(.true.) -type(shr_infnan_inf_type), public, parameter :: shr_infnan_posinf = & - shr_infnan_inf_type(.true.) -type(shr_infnan_inf_type), public, parameter :: shr_infnan_neginf = & - shr_infnan_inf_type(.false.) - -! Bit patterns for implementation without ieee_arithmetic. -! Note that in order to satisfy gfortran's range check, we have to use -! ibset to set the sign bit from a BOZ pattern. -#ifndef HAVE_IEEE_ARITHMETIC -! Single precision. -integer(i4), parameter :: ssnan_pat = int(Z'7FA00000',i4) -integer(i4), parameter :: sqnan_pat = int(Z'7FC00000',i4) -integer(i4), parameter :: sposinf_pat = int(Z'7F800000',i4) -integer(i4), parameter :: sneginf_pat = ibset(sposinf_pat,bit_size(1_i4)-1) -! Double precision. -integer(i8), parameter :: dsnan_pat = int(Z'7FF4000000000000',i8) -integer(i8), parameter :: dqnan_pat = int(Z'7FF8000000000000',i8) -integer(i8), parameter :: dposinf_pat = int(Z'7FF0000000000000',i8) -integer(i8), parameter :: dneginf_pat = ibset(dposinf_pat,bit_size(1_i8)-1) -#endif - - -contains - -!--------------------------------------------------------------------- -! TEST FUNCTIONS -!--------------------------------------------------------------------- -! The "isinf" function simply calls "isposinf" and "isneginf". -!--------------------------------------------------------------------- - -! TYPE double,real - -elemental function shr_infnan_isinf_double(x) result(isinf) - real(r8), intent(in) :: x - logical :: isinf - - isinf = shr_infnan_isposinf(x) .or. shr_infnan_isneginf(x) - - -end function shr_infnan_isinf_double -! TYPE double,real - -elemental function shr_infnan_isinf_real(x) result(isinf) - real(r4), intent(in) :: x - logical :: isinf - - isinf = shr_infnan_isposinf(x) .or. shr_infnan_isneginf(x) - - -end function shr_infnan_isinf_real - -#ifdef HAVE_IEEE_ARITHMETIC - -!--------------------------------------------------------------------- -! The "isposinf" and "isneginf" functions get the IEEE class of a -! real, and test to see if the class is equal to ieee_positive_inf -! or ieee_negative_inf. -!--------------------------------------------------------------------- - -! TYPE double,real - -elemental function shr_infnan_isposinf_double(x) result(isposinf) - use, intrinsic :: ieee_arithmetic, only: & - ieee_class, & - ieee_positive_inf, & - operator(==) - real(r8), intent(in) :: x - logical :: isposinf - - isposinf = (ieee_positive_inf == ieee_class(x)) - - -end function shr_infnan_isposinf_double -! TYPE double,real - -elemental function shr_infnan_isposinf_real(x) result(isposinf) - use, intrinsic :: ieee_arithmetic, only: & - ieee_class, & - ieee_positive_inf, & - operator(==) - real(r4), intent(in) :: x - logical :: isposinf - - isposinf = (ieee_positive_inf == ieee_class(x)) - - -end function shr_infnan_isposinf_real - -! TYPE double,real - -elemental function shr_infnan_isneginf_double(x) result(isneginf) - use, intrinsic :: ieee_arithmetic, only: & - ieee_class, & - ieee_negative_inf, & - operator(==) - real(r8), intent(in) :: x - logical :: isneginf - - isneginf = (ieee_negative_inf == ieee_class(x)) - - -end function shr_infnan_isneginf_double -! TYPE double,real - -elemental function shr_infnan_isneginf_real(x) result(isneginf) - use, intrinsic :: ieee_arithmetic, only: & - ieee_class, & - ieee_negative_inf, & - operator(==) - real(r4), intent(in) :: x - logical :: isneginf - - isneginf = (ieee_negative_inf == ieee_class(x)) - - -end function shr_infnan_isneginf_real - -#else -! Don't have ieee_arithmetic. - -#ifdef CPRGNU -! NaN testing on gfortran. -! TYPE double,real - -elemental function shr_infnan_isnan_double(x) result(is_nan) - real(r8), intent(in) :: x - logical :: is_nan - - is_nan = isnan(x) - - -end function shr_infnan_isnan_double -! TYPE double,real - -elemental function shr_infnan_isnan_real(x) result(is_nan) - real(r4), intent(in) :: x - logical :: is_nan - - is_nan = isnan(x) - - -end function shr_infnan_isnan_real -! End GNU section. -#endif - -!--------------------------------------------------------------------- -! The "isposinf" and "isneginf" functions just test against a known -! bit pattern if we don't have ieee_arithmetic. -!--------------------------------------------------------------------- - -! TYPE double,real - -elemental function shr_infnan_isposinf_double(x) result(isposinf) - real(r8), intent(in) :: x - logical :: isposinf -#if (102 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat -#endif - - isposinf = (x == transfer(posinf_pat,x)) - - -end function shr_infnan_isposinf_double -! TYPE double,real - -elemental function shr_infnan_isposinf_real(x) result(isposinf) - real(r4), intent(in) :: x - logical :: isposinf -#if (101 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat -#endif - - isposinf = (x == transfer(posinf_pat,x)) - - -end function shr_infnan_isposinf_real - -! TYPE double,real - -elemental function shr_infnan_isneginf_double(x) result(isneginf) - real(r8), intent(in) :: x - logical :: isneginf -#if (102 == TYPEREAL) - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif - - isneginf = (x == transfer(neginf_pat,x)) - - -end function shr_infnan_isneginf_double -! TYPE double,real - -elemental function shr_infnan_isneginf_real(x) result(isneginf) - real(r4), intent(in) :: x - logical :: isneginf -#if (101 == TYPEREAL) - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif - - isneginf = (x == transfer(neginf_pat,x)) - - -end function shr_infnan_isneginf_real - -! End ieee_arithmetic conditional. -#endif - -!--------------------------------------------------------------------- -! GENERATION FUNCTIONS -!--------------------------------------------------------------------- -! Two approaches for generation of NaN and Inf values: -! 1. With Fortran 2003, use the ieee_value intrinsic to get a value -! from the corresponding class. These are: -! - ieee_signaling_nan -! - ieee_quiet_nan -! - ieee_positive_inf -! - ieee_negative_inf -! 2. Without Fortran 2003, set the IEEE bit patterns directly. -! Use BOZ literals to get an integer with the correct bit -! pattern, then use "transfer" to transfer those bits into a -! real. -!--------------------------------------------------------------------- - -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_0d_double(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r8), intent(out) :: output - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_0d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_1d_double(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r8), intent(out) :: output(:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_1d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_2d_double(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r8), intent(out) :: output(:,:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_2d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_3d_double(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r8), intent(out) :: output(:,:,:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_3d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_4d_double(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r8), intent(out) :: output(:,:,:,:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_4d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_5d_double(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r8), intent(out) :: output(:,:,:,:,:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_5d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_6d_double(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r8), intent(out) :: output(:,:,:,:,:,:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_6d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_7d_double(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r8), intent(out) :: output(:,:,:,:,:,:,:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_7d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_0d_real(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r4), intent(out) :: output - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_0d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_1d_real(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r4), intent(out) :: output(:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_1d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_2d_real(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r4), intent(out) :: output(:,:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_2d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_3d_real(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r4), intent(out) :: output(:,:,:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_3d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_4d_real(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r4), intent(out) :: output(:,:,:,:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_4d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_5d_real(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r4), intent(out) :: output(:,:,:,:,:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_5d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_6d_real(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r4), intent(out) :: output(:,:,:,:,:,:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_6d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_nan_7d_real(output, nan) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_signaling_nan, & - ieee_quiet_nan, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: snan_pat = ssnan_pat - integer(i4), parameter :: qnan_pat = sqnan_pat -#else - integer(i8), parameter :: snan_pat = dsnan_pat - integer(i8), parameter :: qnan_pat = dqnan_pat -#endif -#endif - real(r4), intent(out) :: output(:,:,:,:,:,:,:) - type(shr_infnan_nan_type), intent(in) :: nan - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (nan%quiet) then - tmp = ieee_value(tmp, ieee_quiet_nan) - else - tmp = ieee_value(tmp, ieee_signaling_nan) - end if -#else - if (nan%quiet) then - tmp = transfer(qnan_pat, tmp) - else - tmp = transfer(snan_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_nan_7d_real - -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_0d_double(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r8), intent(out) :: output - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_0d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_1d_double(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r8), intent(out) :: output(:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_1d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_2d_double(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r8), intent(out) :: output(:,:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_2d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_3d_double(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r8), intent(out) :: output(:,:,:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_3d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_4d_double(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r8), intent(out) :: output(:,:,:,:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_4d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_5d_double(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r8), intent(out) :: output(:,:,:,:,:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_5d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_6d_double(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r8), intent(out) :: output(:,:,:,:,:,:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_6d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_7d_double(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (102 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r8), intent(out) :: output(:,:,:,:,:,:,:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r8) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_7d_double -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_0d_real(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r4), intent(out) :: output - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_0d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_1d_real(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r4), intent(out) :: output(:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_1d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_2d_real(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r4), intent(out) :: output(:,:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_2d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_3d_real(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r4), intent(out) :: output(:,:,:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_3d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_4d_real(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r4), intent(out) :: output(:,:,:,:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_4d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_5d_real(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r4), intent(out) :: output(:,:,:,:,:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_5d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_6d_real(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r4), intent(out) :: output(:,:,:,:,:,:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_6d_real -! TYPE double,real -! DIMS 0,1,2,3,4,5,6,7 - -pure subroutine set_inf_7d_real(output, inf) -#ifdef HAVE_IEEE_ARITHMETIC - use, intrinsic :: ieee_arithmetic, only: & - ieee_positive_inf, & - ieee_negative_inf, & - ieee_value -#else -#if (101 == TYPEREAL) - integer(i4), parameter :: posinf_pat = sposinf_pat - integer(i4), parameter :: neginf_pat = sneginf_pat -#else - integer(i8), parameter :: posinf_pat = dposinf_pat - integer(i8), parameter :: neginf_pat = dneginf_pat -#endif -#endif - real(r4), intent(out) :: output(:,:,:,:,:,:,:) - type(shr_infnan_inf_type), intent(in) :: inf - - ! Use scalar temporary for performance reasons, to reduce the cost of - ! the ieee_value call. - real(r4) :: tmp - -#ifdef HAVE_IEEE_ARITHMETIC - if (inf%positive) then - tmp = ieee_value(tmp,ieee_positive_inf) - else - tmp = ieee_value(tmp,ieee_negative_inf) - end if -#else - if (inf%positive) then - tmp = transfer(posinf_pat, tmp) - else - tmp = transfer(neginf_pat, tmp) - end if -#endif - - output = tmp - - -end subroutine set_inf_7d_real - -!--------------------------------------------------------------------- -! CONVERSION INTERFACES. -!--------------------------------------------------------------------- -! Function methods to get reals from nan/inf types. -!--------------------------------------------------------------------- - - -pure function nan_r8(nan) result(output) - class(shr_infnan_nan_type), intent(in) :: nan - real(r8) :: output - - output = nan - - -end function nan_r8 - - -pure function nan_r4(nan) result(output) - class(shr_infnan_nan_type), intent(in) :: nan - real(r4) :: output - - output = nan - - -end function nan_r4 - - -pure function inf_r8(inf) result(output) - class(shr_infnan_inf_type), intent(in) :: inf - real(r8) :: output - - output = inf - - -end function inf_r8 - - -pure function inf_r4(inf) result(output) - class(shr_infnan_inf_type), intent(in) :: inf - real(r4) :: output - - output = inf - - -end function inf_r4 - -end module shr_infnan_mod diff --git a/test/include/shr_kind_mod.F90 b/test/include/shr_kind_mod.F90 deleted file mode 100644 index e9e7d170..00000000 --- a/test/include/shr_kind_mod.F90 +++ /dev/null @@ -1,19 +0,0 @@ -MODULE shr_kind_mod - - !---------------------------------------------------------------------------- - ! precision/kind constants add data public - !---------------------------------------------------------------------------- - public - integer,parameter :: SHR_KIND_R8 = selected_real_kind(12) ! 8 byte real - integer,parameter :: SHR_KIND_R4 = selected_real_kind( 6) ! 4 byte real - integer,parameter :: SHR_KIND_RN = kind(1.0) ! native real - integer,parameter :: SHR_KIND_I8 = selected_int_kind (13) ! 8 byte integer - integer,parameter :: SHR_KIND_I4 = selected_int_kind ( 6) ! 4 byte integer - integer,parameter :: SHR_KIND_IN = kind(1) ! native integer - integer,parameter :: SHR_KIND_CS = 80 ! short char - integer,parameter :: SHR_KIND_CM = 160 ! mid-sized char - integer,parameter :: SHR_KIND_CL = 256 ! long char - integer,parameter :: SHR_KIND_CX = 512 ! extra-long char - integer,parameter :: SHR_KIND_CXX= 4096 ! extra-extra-long char - -END MODULE shr_kind_mod diff --git a/test/include/spmd_utils.F90 b/test/include/spmd_utils.F90 deleted file mode 100644 index c827ac56..00000000 --- a/test/include/spmd_utils.F90 +++ /dev/null @@ -1,11 +0,0 @@ -module spmd_utils - - implicit none - private - - integer, parameter, public :: masterprocid = 0 - integer, parameter, public :: iam = 0 - integer, parameter, public :: npes = 1 - logical, parameter, public :: masterproc = .true. - -end module spmd_utils