From 5d10dfd8373fc59b23594f140c30bcbf62c60ec3 Mon Sep 17 00:00:00 2001 From: John Truesdale Date: Tue, 30 Apr 2024 00:24:20 -0600 Subject: [PATCH] dadadj ccpp updates --- bld/configure | 1 + src/control/cam_snapshot_common.F90 | 2 +- src/physics/cam/dadadj.F90 | 174 ---------------------------- src/physics/cam/dadadj_cam.F90 | 98 +++++++++------- src/physics/cam/physpkg.F90 | 4 +- src/physics/cam_dev/physpkg.F90 | 4 +- 6 files changed, 61 insertions(+), 222 deletions(-) delete mode 100644 src/physics/cam/dadadj.F90 diff --git a/bld/configure b/bld/configure index e2fb784495..05f1d7155c 100755 --- a/bld/configure +++ b/bld/configure @@ -2304,6 +2304,7 @@ sub write_filepath #Add the CCPP'ized subdirectories print $fh "$camsrcdir/src/atmos_phys/zhang_mcfarlane\n"; + print $fh "$camsrcdir/src/atmos_phys/dadadj\n"; # Dynamics package and test utilities print $fh "$camsrcdir/src/dynamics/$dyn\n"; diff --git a/src/control/cam_snapshot_common.F90 b/src/control/cam_snapshot_common.F90 index f2a4780619..81e8694006 100644 --- a/src/control/cam_snapshot_common.F90 +++ b/src/control/cam_snapshot_common.F90 @@ -86,7 +86,7 @@ module cam_snapshot_common type (snapshot_type) :: tend_snapshot(6) type (snapshot_type) :: cam_in_snapshot(30) type (snapshot_type) :: cam_out_snapshot(30) -type (snapshot_type_nd) :: pbuf_snapshot(250) +type (snapshot_type_nd) :: pbuf_snapshot(300) contains diff --git a/src/physics/cam/dadadj.F90 b/src/physics/cam/dadadj.F90 deleted file mode 100644 index b9762f8f5f..0000000000 --- a/src/physics/cam/dadadj.F90 +++ /dev/null @@ -1,174 +0,0 @@ -module dadadj -!----------------------------------------------------------------------- -! -! Purpose: -! GFDL style dry adiabatic adjustment -! -! Method: -! if stratification is unstable, adjustment to the dry adiabatic lapse -! rate is forced subject to the condition that enthalpy is conserved. -! -! Author: J.Hack -! -!----------------------------------------------------------------------- - -use shr_kind_mod, only: r8 => shr_kind_r8 - -implicit none -private -save - -public :: & - dadadj_initial, & - dadadj_calc - -integer :: nlvdry ! number of layers from top of model to apply the adjustment -integer :: niter ! number of iterations for convergence - -!=============================================================================== -contains -!=============================================================================== - -subroutine dadadj_initial(nlvdry_in, niter_in) - - integer, intent(in) :: nlvdry_in - integer, intent(in) :: niter_in - - nlvdry = nlvdry_in - niter = niter_in - -end subroutine dadadj_initial - -!=============================================================================== - -subroutine dadadj_calc( & - ncol, pmid, pint, pdel, cappav, t, & - q, dadpdf, icol_err) - - ! Arguments - - integer, intent(in) :: ncol ! number of atmospheric columns - - real(r8), intent(in) :: pmid(:,:) ! pressure at model levels - real(r8), intent(in) :: pint(:,:) ! pressure at model interfaces - real(r8), intent(in) :: pdel(:,:) ! vertical delta-p - real(r8), intent(in) :: cappav(:,:) ! variable Kappa - - real(r8), intent(inout) :: t(:,:) ! temperature (K) - real(r8), intent(inout) :: q(:,:) ! specific humidity - - real(r8), intent(out) :: dadpdf(:,:) ! PDF of where adjustments happened - - integer, intent(out) :: icol_err ! index of column in which error occurred - - !---------------------------Local workspace----------------------------- - - integer :: i,k ! longitude, level indices - integer :: jiter ! iteration index - - real(r8), allocatable :: c1dad(:) ! intermediate constant - real(r8), allocatable :: c2dad(:) ! intermediate constant - real(r8), allocatable :: c3dad(:) ! intermediate constant - real(r8), allocatable :: c4dad(:) ! intermediate constant - real(r8) :: gammad ! dry adiabatic lapse rate (deg/Pa) - real(r8) :: zeps ! convergence criterion (deg/Pa) - real(r8) :: rdenom ! reciprocal of denominator of expression - real(r8) :: dtdp ! delta-t/delta-p - real(r8) :: zepsdp ! zeps*delta-p - real(r8) :: zgamma ! intermediate constant - real(r8) :: qave ! mean q between levels - real(r8) :: cappa ! Kappa at level intefaces - - logical :: ilconv ! .TRUE. ==> convergence was attained - logical :: dodad(ncol) ! .TRUE. ==> do dry adjustment - - !----------------------------------------------------------------------- - - icol_err = 0 - zeps = 2.0e-5_r8 ! set convergence criteria - - allocate(c1dad(nlvdry), c2dad(nlvdry), c3dad(nlvdry), c4dad(nlvdry)) - - ! Find gridpoints with unstable stratification - - do i = 1, ncol - cappa = 0.5_r8*(cappav(i,2) + cappav(i,1)) - gammad = cappa*0.5_r8*(t(i,2) + t(i,1))/pint(i,2) - dtdp = (t(i,2) - t(i,1))/(pmid(i,2) - pmid(i,1)) - dodad(i) = (dtdp + zeps) .gt. gammad - end do - - dadpdf(:ncol,:) = 0._r8 - do k= 2, nlvdry - do i = 1, ncol - cappa = 0.5_r8*(cappav(i,k+1) + cappav(i,k)) - gammad = cappa*0.5_r8*(t(i,k+1) + t(i,k))/pint(i,k+1) - dtdp = (t(i,k+1) - t(i,k))/(pmid(i,k+1) - pmid(i,k)) - dodad(i) = dodad(i) .or. (dtdp + zeps).gt.gammad - if ((dtdp + zeps).gt.gammad) then - dadpdf(i,k) = 1._r8 - end if - end do - end do - - ! Make a dry adiabatic adjustment - ! Note: nlvdry ****MUST**** be < pver - - COL: do i = 1, ncol - - if (dodad(i)) then - - zeps = 2.0e-5_r8 - - do k = 1, nlvdry - c1dad(k) = cappa*0.5_r8*(pmid(i,k+1)-pmid(i,k))/pint(i,k+1) - c2dad(k) = (1._r8 - c1dad(k))/(1._r8 + c1dad(k)) - rdenom = 1._r8/(pdel(i,k)*c2dad(k) + pdel(i,k+1)) - c3dad(k) = rdenom*pdel(i,k) - c4dad(k) = rdenom*pdel(i,k+1) - end do - -50 continue - - do jiter = 1, niter - ilconv = .true. - - do k = 1, nlvdry - zepsdp = zeps*(pmid(i,k+1) - pmid(i,k)) - zgamma = c1dad(k)*(t(i,k) + t(i,k+1)) - - if ((t(i,k+1)-t(i,k)) >= (zgamma+zepsdp)) then - ilconv = .false. - t(i,k+1) = t(i,k)*c3dad(k) + t(i,k+1)*c4dad(k) - t(i,k) = c2dad(k)*t(i,k+1) - qave = (pdel(i,k+1)*q(i,k+1) + pdel(i,k)*q(i,k))/(pdel(i,k+1)+ pdel(i,k)) - q(i,k+1) = qave - q(i,k) = qave - end if - - end do - - if (ilconv) cycle COL ! convergence => next longitude - end do - - ! Double convergence criterion if no convergence in niter iterations - - zeps = zeps + zeps - if (zeps > 1.e-4_r8) then - icol_err = i - return ! error return - else - go to 50 - end if - - end if - - end do COL - - deallocate(c1dad, c2dad, c3dad, c4dad) - -end subroutine dadadj_calc - -!=============================================================================== - -end module dadadj diff --git a/src/physics/cam/dadadj_cam.F90 b/src/physics/cam/dadadj_cam.F90 index 0717865ca8..d12c307ac9 100644 --- a/src/physics/cam/dadadj_cam.F90 +++ b/src/physics/cam/dadadj_cam.F90 @@ -2,7 +2,7 @@ module dadadj_cam ! CAM interfaces for the dry adiabatic adjustment parameterization -use shr_kind_mod, only: r8=>shr_kind_r8, cs=>shr_kind_cs +use shr_kind_mod, only: r8=>shr_kind_r8, cs=>shr_kind_cs, cm=>shr_kind_cm use ppgrid, only: pcols, pver, pverp use constituents, only: pcnst use air_composition, only: cappav, cpairv @@ -17,7 +17,7 @@ module dadadj_cam use namelist_utils, only: find_group_name use units, only: getunit, freeunit -use dadadj, only: dadadj_initial, dadadj_calc +use dadadj, only: dadadj_init, dadadj_run implicit none private @@ -25,7 +25,7 @@ module dadadj_cam public :: & dadadj_readnl, & - dadadj_init, & + dadadj_cam_init, & dadadj_tend ! Namelist variables @@ -42,8 +42,10 @@ subroutine dadadj_readnl(filein) namelist /dadadj_nl/ dadadj_nlvdry, dadadj_niter - integer :: unitn, ierr - character(len=*), parameter :: sub='dadadj_readnl' + integer :: unitn, ierr + integer :: errflg ! CCPP physics scheme error flag + character(len=512) :: errmsg ! CCPP physics scheme error message + character(len=*), parameter :: sub='dadadj_readnl' !------------------------------------------------------------------ ! Read namelist @@ -67,13 +69,16 @@ subroutine dadadj_readnl(filein) call mpibcast(dadadj_niter, 1, mpi_integer, masterprocid, mpicom) #endif - call dadadj_initial(dadadj_nlvdry, dadadj_niter) + call dadadj_init(dadadj_nlvdry, dadadj_niter, pver, errmsg, errflg) + if (errflg /=0) then + call endrun('dadadj_readnl: Error returned from dadadj_init: '//trim(errmsg)) + end if if (masterproc .and. .not. use_simple_phys) then write(iulog,*)'Dry adiabatic adjustment applied to top N layers; N=', & - dadadj_nlvdry + dadadj_nlvdry write(iulog,*)'Dry adiabatic adjustment number of iterations for convergence =', & - dadadj_niter + dadadj_niter end if end subroutine dadadj_readnl @@ -81,12 +86,12 @@ end subroutine dadadj_readnl !=============================================================================== -subroutine dadadj_init() +subroutine dadadj_cam_init() use cam_history, only: addfld call addfld('DADADJ_PD', (/ 'lev' /), 'A', 'probability', 'dry adiabatic adjustment probability') -end subroutine dadadj_init +end subroutine dadadj_cam_init !=============================================================================== @@ -98,39 +103,46 @@ subroutine dadadj_tend(dt, state, ptend) type(physics_state), intent(in) :: state ! Physics state variables type(physics_ptend), intent(out) :: ptend ! parameterization tendencies - logical :: lq(pcnst) - real(r8) :: dadpdf(pcols, pver) - integer :: ncol, lchnk, icol_err - character(len=128) :: errstring ! Error string - - ncol = state%ncol - lchnk = state%lchnk - lq(:) = .FALSE. - lq(1) = .TRUE. - call physics_ptend_init(ptend, state%psetcols, 'dadadj', ls=.true., lq=lq) - - ! use the ptend components for temporary storate and copy state info for input to - ! dadadj_calc which directly updates the temperature and moisture input arrays. - - ptend%s(:ncol,:pver) = state%t(:ncol,:pver) - ptend%q(:ncol,:pver,1) = state%q(:ncol,:pver,1) - - call dadadj_calc( & - ncol, state%pmid, state%pint, state%pdel, cappav(:,:,lchnk), ptend%s, & - ptend%q(:,:,1), dadpdf, icol_err) - - call outfld('DADADJ_PD', dadpdf(:ncol,:), ncol, lchnk) - - if (icol_err > 0) then - ! error exit - write(errstring, *) & - 'dadadj_calc: No convergence in column at lat,lon:', & - state%lat(icol_err)*180._r8/pi, state%lon(icol_err)*180._r8/pi - call handle_errmsg(errstring, subname="dadadj_tend") - end if - - ptend%s(:ncol,:) = (ptend%s(:ncol,:) - state%t(:ncol,:) )/dt * cpairv(:ncol,:,lchnk) - ptend%q(:ncol,:,1) = (ptend%q(:ncol,:,1) - state%q(:ncol,:,1))/dt + character(len=512) :: errstring ! Error string + character(len=512) :: errmsg ! CCPP physics scheme error message + character(len=64) :: scheme_name! CCPP physics scheme name (not used in CAM) + integer :: icol_err + integer :: lchnk + integer :: ncol + integer :: errflg ! CCPP physics scheme error flag + logical :: lq(pcnst) + real(r8) :: dadpdf(pcols, pver) + + !------------------------------------------------------------------ + ncol = state%ncol + lchnk = state%lchnk + lq(:) = .FALSE. + lq(1) = .TRUE. + call physics_ptend_init(ptend, state%psetcols, 'dadadj', ls=.true., lq=lq) + + !REMOVECAM - no longer need these when CAM is retired and pcols no longer exists + dadpdf = 0._r8 + ptend%s = 0._r8 + ptend%q = 0._r8 + !REMOVECAM_END + + ! dadadj_run returns t tend, we are passing the ptend%s array to receive the t tendency and will convert it to s + ! before it is returned to CAM.. + call dadadj_run( & + ncol, dt, state%pmid(:ncol,:), state%pint(:ncol,:), state%pdel(:ncol,:), state%t(:ncol,:), state%q(:ncol,:,1), cappav(:ncol,:,lchnk), & + ptend%s(:ncol,:), ptend%q(:ncol,:,1), dadpdf(:ncol,:), scheme_name, errmsg, errflg) + + if (errflg /= 0) then + ! error exit + write(errstring, *) errmsg,' at lat,lon:', & + state%lat(errflg)*180._r8/pi, state%lon(errflg)*180._r8/pi + call endrun('dadadj_tend: Error returned from dadadj_run: '//trim(errstring)) + end if + + call outfld('DADADJ_PD', dadpdf(:ncol,:), ncol, lchnk) + + ! convert the t tendency to an s tendency for cam + ptend%s(:ncol,:) = ptend%s(:ncol,:) * cpairv(:ncol,:,lchnk) end subroutine dadadj_tend diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index cb7322254f..883b774fd1 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -777,7 +777,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use sslt_rebin, only: sslt_rebin_init use tropopause, only: tropopause_init use solar_data, only: solar_data_init - use dadadj_cam, only: dadadj_init + use dadadj_cam, only: dadadj_cam_init use cam_abortutils, only: endrun use nudging, only: Nudge_Model, nudging_init use cam_snapshot, only: cam_snapshot_init @@ -952,7 +952,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) #endif call sslt_rebin_init() call tropopause_init() - call dadadj_init() + call dadadj_cam_init() prec_dp_idx = pbuf_get_index('PREC_DP') snow_dp_idx = pbuf_get_index('SNOW_DP') diff --git a/src/physics/cam_dev/physpkg.F90 b/src/physics/cam_dev/physpkg.F90 index aef997716f..45d195de2e 100644 --- a/src/physics/cam_dev/physpkg.F90 +++ b/src/physics/cam_dev/physpkg.F90 @@ -761,7 +761,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use clubb_intr, only: clubb_ini_cam use tropopause, only: tropopause_init use solar_data, only: solar_data_init - use dadadj_cam, only: dadadj_init + use dadadj_cam, only: dadadj_cam_init use cam_abortutils, only: endrun use nudging, only: Nudge_Model, nudging_init use cam_snapshot, only: cam_snapshot_init @@ -920,7 +920,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) call metdata_phys_init() #endif call tropopause_init() - call dadadj_init() + call dadadj_cam_init() prec_dp_idx = pbuf_get_index('PREC_DP') snow_dp_idx = pbuf_get_index('SNOW_DP')