Merge pull request #1235 from mwaxmonsky/feature/pbl-utils-atmos-phys-integration

cam6_4_071: PBL_utils atmospheric_physics integration

Author: mwaxmonsky

.gitmodules

@@ -36,7 +36,7 @@
 [submodule "atmos_phys"]
 	path = src/atmos_phys
 	url = https://github.com/ESCOMP/atmospheric_physics
-        fxtag = atmos_phys0_08_000
+        fxtag = atmos_phys0_09_000
         fxrequired = AlwaysRequired
 	fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics
 

bld/configure

@@ -2149,6 +2149,7 @@ sub write_filepath
     # in the list of filepaths.
     print $fh "$camsrcdir/src/physics/cam\n";
     print $fh "$camsrcdir/src/atmos_phys/to_be_ccppized\n";
+    print $fh "$camsrcdir/src/atmos_phys/phys_utils\n";
 
     #Add the CCPP'ized subdirectories
     print $fh "$camsrcdir/src/atmos_phys/schemes/tropopause_find\n";

doc/ChangeLog

@@ -1,5 +1,80 @@
 ===============================================================
 
+Tag name: cam6_4_071
+Originator(s): mwaxmonsky
+Date: 26 February 2025
+One-line Summary: PBL_utils atmospheric_physics integration
+Github PR URL: https://github.com/ESCOMP/CAM/pull/1235
+
+Purpose of changes (include the issue number and title text for each relevant GitHub issue):
+
+Closes #1235 - Integrates PBL now from atmospheric_physics to enable further
+vertical diffusion CCPP-ization
+
+Describe any changes made to build system: N/A
+
+Describe any changes made to the namelist: N/A
+
+List any changes to the defaults for the boundary datasets: N/A
+
+Describe any substantial timing or memory changes: N/A
+
+Code reviewed by: nusbaume
+
+List all files eliminated: N/A
+
+List all files added and what they do: N/A
+
+List all existing files that have been modified, and describe the changes:
+
+M       .gitmodules
+M       src/atmos_phys
+- Updates atmos_phys submodule
+
+M       bld/configure
+- Add phys utils which now contains refactored PBL utils code
+
+M       src/physics/cam/clubb_intr.F90
+M       src/physics/cam/eddy_diff_cam.F90
+M       src/physics/cam/hb_diff.F90
+M       src/physics/cam/pbl_utils.F90
+M       src/physics/cam/vertical_diffusion.F90
+- Updates old PBL references to new PBL interface API
+
+M       src/physics/cam/physpkg.F90
+M       src/physics/cam7/physpkg.F90
+ - Removing `pbl_utils_init()` call as no longer needed
+
+
+If there were any failures reported from running test_driver.sh on any test
+platform, and checkin with these failures has been OK'd by the gatekeeper,
+then copy the lines from the td.*.status files for the failed tests to the
+appropriate machine below.  All failed tests must be justified.
+
+derecho/intel/aux_cam:
+
+ERP_Ln9.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq9s (Overall: FAIL)
+SMS_Ld1.f09_f09_mg17.FCHIST_GC.derecho_intel.cam-outfrq1d (Overall: DIFF)
+ - pre-existing failures due to HEMCO not having reproducible results (issues #1018 and #856)
+
+SMS_D_Ln9.f19_f19_mg17.FXHIST.derecho_intel.cam-outfrq9s_amie (Overall: FAIL)
+SMS_D_Ln9_P1280x1.ne0CONUSne30x8_ne0CONUSne30x8_mt12.FCHIST.derecho_intel.cam-outfrq9s (Overall: FAIL)
+ - pre-existing failures due to build-namelist error requiring CLM/CTSM external update
+
+derecho/nvhpc/aux_cam: ALL PASS
+
+izumi/nag/aux_cam: ALL PASS
+
+izumi/gnu/aux_cam: ALL PASS
+
+CAM tag used for the baseline comparison tests if different than previous
+tag:
+
+Summarize any changes to answers: b4b
+
+===============================================================
+===============================================================
+
 Tag name: cam6_4_070
 Originator(s): patcal, nusbaume
 Date: 22 February 2025

src/atmos_phys

@@ -20,13 +20,14 @@ module clubb_intr
   use shr_kind_mod,        only: r8=>shr_kind_r8
   use ppgrid,              only: pver, pverp, pcols, begchunk, endchunk
   use phys_control,        only: phys_getopts
-  use physconst,           only: cpair, gravit, rga, latvap, latice, zvir, rh2o, karman, pi
+  use physconst,           only: cpair, gravit, rga, latvap, latice, zvir, rh2o, karman, pi, rair
   use air_composition,     only: rairv, cpairv
   use cam_history_support, only: max_fieldname_len
 
   use spmd_utils,          only: masterproc
   use constituents,        only: pcnst, cnst_add
-  use pbl_utils,           only: calc_ustar, calc_obklen
+  use atmos_phys_pbl_utils,only: calc_friction_velocity, calc_kinematic_heat_flux, calc_ideal_gas_rrho, &
+                                 calc_kinematic_water_vapor_flux, calc_kinematic_buoyancy_flux, calc_obukhov_length
   use ref_pres,            only: top_lev => trop_cloud_top_lev
 
 #ifdef CLUBB_SGS
@@ -3491,8 +3492,8 @@ subroutine clubb_tend_cam( state,   ptend_all,   pbuf,     hdtime, &
         ubar = sqrt(state1%u(i,pver)**2+state1%v(i,pver)**2)
         if (ubar <  0.25_r8) ubar = 0.25_r8
 
-        call calc_ustar( state1%t(i,pver), state1%pmid(i,pver), cam_in%wsx(i), cam_in%wsy(i), &
-              rrho(i), ustar )
+        rrho(i) = calc_ideal_gas_rrho(rair, state1%t(i,pver), state1%pmid(i,pver))
+        ustar   = calc_friction_velocity(cam_in%wsx(i), cam_in%wsy(i), rrho(i))
 
         upwp_sfc(i) = -state1%u(i,pver)*ustar**2/ubar
         vpwp_sfc(i) = -state1%v(i,pver)*ustar**2/ubar
@@ -4706,13 +4707,13 @@ subroutine clubb_tend_cam( state,   ptend_all,   pbuf,     hdtime, &
     enddo
 
     ! diagnose surface friction and obukhov length (inputs to diagnose PBL depth)
-    rrho(1:ncol) = (rga)*(state1%pdel(1:ncol,pver)/dz_g(1:ncol,pver))
-    call calc_ustar( ncol, state1%t(1:ncol,pver), state1%pmid(1:ncol,pver), cam_in%wsx(1:ncol), cam_in%wsy(1:ncol), &
-                    rrho(1:ncol), ustar2(1:ncol))
+    rrho   (1:ncol) = calc_ideal_gas_rrho(rair, state1%t(1:ncol,pver), state1%pmid(1:ncol,pver))
+    ustar2 (1:ncol) = calc_friction_velocity(cam_in%wsx(1:ncol), cam_in%wsy(1:ncol), rrho(1:ncol))
     ! use correct qflux from coupler
-    call calc_obklen( ncol, th(1:ncol,pver), thv(1:ncol,pver), cam_in%cflx(1:ncol,1), cam_in%shf(1:ncol), &
-                      rrho(1:ncol), ustar2(1:ncol), kinheat(1:ncol), kinwat(1:ncol), kbfs(1:ncol), &
-                      obklen(1:ncol))
+    kinheat(1:ncol) = calc_kinematic_heat_flux(cam_in%shf(1:ncol), rrho(1:ncol), cpair)
+    kinwat (1:ncol) = calc_kinematic_water_vapor_flux(cam_in%cflx(1:ncol,1), rrho(1:ncol))
+    kbfs   (1:ncol) = calc_kinematic_buoyancy_flux(kinheat(1:ncol), zvir, th(1:ncol,pver), kinwat(1:ncol))
+    obklen (1:ncol) = calc_obukhov_length(thv(1:ncol,pver), ustar2(1:ncol), gravit, karman, kbfs(1:ncol))
 
     dummy2(:) = 0._r8
     dummy3(:) = 0._r8

src/physics/cam/clubb_intr.F90

@@ -424,16 +424,16 @@ subroutine compute_eddy_diff( pbuf, lchnk  ,
   !                       May.    2008.                                 !
   !-------------------------------------------------------------------- !
 
-  use diffusion_solver, only: compute_vdiff
-  use cam_history,      only: outfld
-  use phys_debug_util,  only: phys_debug_col
-  use air_composition,  only: cpairv
-  use pbl_utils,        only: calc_ustar, austausch_atm
-  use error_messages,   only: handle_errmsg
-  use coords_1d,        only: Coords1D
-  use wv_saturation,    only: qsat
-  use eddy_diff,        only: trbintd, caleddy
-  use physics_buffer,   only: pbuf_get_field
+  use diffusion_solver,     only: compute_vdiff
+  use cam_history,          only: outfld
+  use phys_debug_util,      only: phys_debug_col
+  use air_composition,      only: cpairv
+  use atmos_phys_pbl_utils, only: calc_eddy_flux_coefficient, calc_ideal_gas_rrho, calc_friction_velocity
+  use error_messages,       only: handle_errmsg
+  use coords_1d,            only: Coords1D
+  use wv_saturation,        only: qsat
+  use eddy_diff,            only: trbintd, caleddy
+  use physics_buffer,       only: pbuf_get_field
 
   ! --------------- !
   ! Input Variables !
@@ -670,10 +670,11 @@ subroutine compute_eddy_diff( pbuf, lchnk  ,
      ! I am using updated wind, here.
 
      ! Compute ustar
-     call calc_ustar( ncol, tfd(:ncol,pver), pmid(:ncol,pver), &
-                      taux(:ncol) - ksrftms(:ncol) * ufd(:ncol,pver), & ! Zonal wind stress
-                      tauy(:ncol) - ksrftms(:ncol) * vfd(:ncol,pver), & ! Meridional wind stress
-                      rrho(:ncol), ustar(:ncol))
+     rrho(:ncol)   = calc_ideal_gas_rrho(rair, tfd(:ncol,pver), pmid(:ncol,pver))
+     ustar(:ncol)  = calc_friction_velocity(taux(:ncol) - ksrftms(:ncol) * ufd(:ncol,pver), & ! Zonal wind stress
+                                            tauy(:ncol) - ksrftms(:ncol) * vfd(:ncol,pver), & ! Meridional wind stress
+                                            rrho(:ncol))
+
      minpblh(:ncol) = 100.0_r8 * ustar(:ncol)   ! By construction, 'minpblh' is larger than 1 [m] when 'ustar_min = 0.01'.
 
      ! Calculate (qt,sl,n2,s2,ri) from a given set of (t,qv,ql,qi,u,v)
@@ -694,8 +695,12 @@ subroutine compute_eddy_diff( pbuf, lchnk  ,
 
      ! Get free atmosphere exchange coefficients. This 'kvf' is not used in UW moist PBL scheme
      if (use_kvf) then
-        call austausch_atm(pcols, ncol, pver, ntop_eddy, nbot_eddy, &
-             ml2, ri, s2, kvf )
+        kvf(:ncol,:) = 0.0_r8
+        do k = ntop_eddy, nbot_eddy-1
+           do i = 1, ncol
+              kvf(i,k+1) = calc_eddy_flux_coefficient(ml2(k), ri(i, k), s2(i, k))
+           end do
+        end do
      else
         kvf = 0._r8
      end if

src/physics/cam/eddy_diff_cam.F90

@@ -11,7 +11,6 @@ module hb_diff
   ! Private methods:
   !       trbintd         initializes time dependent variables
   !       pblintd         initializes time dependent variables that depend pbl depth
-  !       austausch_atm   computes free atmosphere exchange coefficients
   !       austausch_pbl   computes pbl exchange coefficients
   !
   !---------------------------Code history--------------------------------
@@ -150,13 +149,16 @@ subroutine compute_hb_diff(ncol                  , &
     ! 
     !-----------------------------------------------------------------------
 
-    use pbl_utils, only: virtem, calc_ustar, calc_obklen, austausch_atm
+    use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length,               &
+                                    calc_eddy_flux_coefficient, calc_ideal_gas_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, &
+                                    calc_kinematic_buoyancy_flux
+    use physconst,            only: zvir, rair, gravit, karman
 
     !------------------------------Arguments--------------------------------
     !
     ! Input arguments
     !
-    integer, intent(in) :: ncol                       ! number of atmospheric columns
+    integer, intent(in)   :: ncol                     ! number of atmospheric columns
 
     real(r8), intent(in)  :: th(pcols,pver)           ! potential temperature [K]
     real(r8), intent(in)  :: t(pcols,pver)            ! temperature (used for density)
@@ -203,20 +205,22 @@ subroutine compute_hb_diff(ncol                  , &
     real(r8) :: n2(pcols,pver)          ! brunt vaisaila frequency
     real(r8) :: bge(pcols)              ! buoyancy gradient enhancment
     integer  :: ktopbl(pcols)           ! index of first midpoint inside pbl
+    integer  :: i,k
     !
     ! Initialize time dependent variables that do not depend on pbl height
     !
 
     ! virtual temperature
-    call virtem(ncol, (pver-ntop_turb+1), th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), thv(:ncol,ntop_turb:))
+    thv(:ncol,ntop_turb:) = calc_virtual_temperature(th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), zvir)
 
     ! Compute ustar, Obukhov length, and kinematic surface fluxes.
-    call calc_ustar(ncol, t(:ncol,pver),pmid(:ncol,pver),taux(:ncol),tauy(:ncol), &
-         rrho(:ncol),ustar(:ncol))
-    call calc_obklen(ncol, th(:ncol,pver), thv(:ncol,pver), qflx(:ncol),  &
-                     shflx(:ncol),   rrho(:ncol),     ustar(:ncol), &
-                     khfs(:ncol),    kqfs(:ncol),     kbfs(:ncol),  &
-                     obklen(:ncol))
+    rrho(:ncol)   = calc_ideal_gas_rrho(rair, t(:ncol,pver), pmid(:ncol,pver))
+    ustar(:ncol)  = calc_friction_velocity(taux(:ncol),tauy(:ncol), rrho(:ncol))
+    khfs(:ncol)   = calc_kinematic_heat_flux(shflx(:ncol), rrho(:ncol), cpair)
+    kqfs(:ncol)   = calc_kinematic_water_vapor_flux(qflx(:ncol), rrho(:ncol))
+    kbfs(:ncol)   = calc_kinematic_buoyancy_flux(khfs(:ncol), zvir, th(:ncol,pver), kqfs(:ncol))
+    obklen(:ncol) = calc_obukhov_length(thv(:ncol,pver), ustar(:ncol), gravit, karman, kbfs(:ncol))
+
     ! Calculate s2, n2, and Richardson number.
     call trbintd(ncol    ,                            &
          thv     ,z       ,u       ,v       , &
@@ -229,10 +233,15 @@ subroutine compute_hb_diff(ncol                  , &
          ustar   ,obklen  ,kbfs    ,pblh    ,wstar   , &
          zi      ,cldn    ,ocnfrac ,bge     )
     !
-    ! Get free atmosphere exchange coefficients
+    ! Get atmosphere exchange coefficients
     !
-    call austausch_atm(pcols, ncol, pver, ntop_turb, nbot_turb, &
-         ml2, ri, s2, kvf)
+    kvf(:ncol,:) = 0.0_r8
+    do k = ntop_turb, nbot_turb-1
+       do i = 1, ncol
+          kvf(i,k+1) = calc_eddy_flux_coefficient(ml2(k), ri(i, k), s2(i, k))
+       end do
+    end do
+
     ! 
     ! Get pbl exchange coefficients
     !
@@ -262,7 +271,10 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     !
     !-----------------------------------------------------------------------
 
-    use pbl_utils, only: virtem, calc_ustar, calc_obklen, austausch_atm_free
+    use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length,                    &
+                                    calc_free_atm_eddy_flux_coefficient, calc_ideal_gas_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, &
+                                    calc_kinematic_buoyancy_flux
+    use physconst,            only: zvir, rair, gravit, karman
 
     !------------------------------Arguments--------------------------------
     !
@@ -303,26 +315,31 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     real(r8) :: kvf(pcols,pverp)        ! free atmospheric eddy diffsvty [m2/s]
     real(r8) :: s2(pcols,pver)          ! shear squared
     real(r8) :: n2(pcols,pver)          ! brunt vaisaila frequency
+    integer  :: i, k
 
     ! virtual potential temperature
-    call virtem(ncol, (pver-ntop_turb+1), th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), thv(:ncol,ntop_turb:))
+    thv(:ncol,ntop_turb:) = calc_virtual_temperature(th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), zvir)
 
     ! Compute ustar, Obukhov length, and kinematic surface fluxes.
-    call calc_ustar(ncol, t(:ncol,pver),pmid(:ncol,pver),taux(:ncol),tauy(:ncol), &
-         rrho(:ncol),ustar(:ncol))
-    call calc_obklen(ncol, th(:ncol,pver), thv(:ncol,pver), qflx(:ncol),  &
-                     shflx(:ncol),   rrho(:ncol),     ustar(:ncol), &
-                     khfs(:ncol),    kqfs(:ncol),     kbfs(:ncol),  &
-                     obklen(:ncol))
+    rrho(:ncol)   = calc_ideal_gas_rrho(rair, t(:ncol,pver), pmid(:ncol,pver))
+    ustar(:ncol)  = calc_friction_velocity(taux(:ncol),tauy(:ncol), rrho(:ncol))
+    khfs(:ncol)   = calc_kinematic_heat_flux(shflx(:ncol), rrho(:ncol), cpair)
+    kqfs(:ncol)   = calc_kinematic_water_vapor_flux(qflx(:ncol), rrho(:ncol))
+    kbfs(:ncol)   = calc_kinematic_buoyancy_flux(khfs(:ncol), zvir, th(:ncol,pver), kqfs(:ncol))
+    obklen(:ncol) = calc_obukhov_length(thv(:ncol,pver), ustar(:ncol), gravit, karman, kbfs(:ncol))
     ! Calculate s2, n2, and Richardson number.
     call trbintd(ncol    ,                            &
          thv     ,z       ,u       ,v       , &
          s2      ,n2      ,ri      )
     !
     ! Get free atmosphere exchange coefficients
     !
-    call austausch_atm_free(pcols, ncol, pver, ntop_turb, nbot_turb, &
-         ml2, ri, s2, kvf)
+    kvf(:ncol,:) = 0.0_r8
+    do k = ntop_turb, nbot_turb - 1
+       do i = 1, ncol
+          kvf(i,k+1) = calc_free_atm_eddy_flux_coefficient(ml2(k), ri(i, k), s2(i, k))
+       end do
+    end do
 
     kvq(:ncol,:) = kvf(:ncol,:)
     kvm(:ncol,:) = kvf(:ncol,:)