From 514718ef9b0da5cd2d819e46e2bc6ca3167dd899 Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Wed, 15 Jan 2025 16:04:15 -0700
Subject: [PATCH 01/19] Initial attempt at calling atmos_phys physics utils
 layer.

---
 src/physics/cam/hb_diff.F90 | 47 ++++++++++++++++++++-----------------
 1 file changed, 25 insertions(+), 22 deletions(-)

diff --git a/src/physics/cam/hb_diff.F90 b/src/physics/cam/hb_diff.F90
index a3bb11a17d..306caa95e8 100644
--- a/src/physics/cam/hb_diff.F90
+++ b/src/physics/cam/hb_diff.F90
@@ -150,13 +150,14 @@ 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_ustar, calc_obukhov_length, austausch_atm, calc_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)
@@ -208,15 +209,16 @@ subroutine compute_hb_diff(ncol                  , &
     !
 
     ! 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_rrho(rair, t(:ncol,pver), pmid(:ncol,pver))
+    ustar(:ncol)  = calc_ustar(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 +231,10 @@ 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 = austausch_atm(pcols, ncol, pver, ntop_turb, nbot_turb, &
+                        ml2, ri, s2)
     ! 
     ! Get pbl exchange coefficients
     !
@@ -262,7 +264,8 @@ 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_ustar, calc_obukhov_length, austausch_atm_free, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, calc_kinematic_buoyancy_flux
+    use physconst,            only: zvir, rair, gravit, karman
 
     !------------------------------Arguments--------------------------------
     !
@@ -305,15 +308,15 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     real(r8) :: n2(pcols,pver)          ! brunt vaisaila frequency
 
     ! 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_rrho(rair, t(:ncol,pver), pmid(:ncol,pver))
+    ustar(:ncol)  = calc_ustar(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       , &
@@ -321,8 +324,8 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     !
     ! Get free atmosphere exchange coefficients
     !
-    call austausch_atm_free(pcols, ncol, pver, ntop_turb, nbot_turb, &
-         ml2, ri, s2, kvf)
+    kvf = austausch_atm_free(pcols, ncol, pver, ntop_turb, nbot_turb, &
+                             ml2, ri, s2)
 
     kvq(:ncol,:) = kvf(:ncol,:)
     kvm(:ncol,:) = kvf(:ncol,:)

From beb93e1a2357b87bba2924ea314e138a4d6375c4 Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Tue, 21 Jan 2025 06:46:35 -0700
Subject: [PATCH 02/19] Initial integration of elemental austausch.

---
 src/physics/cam/hb_diff.F90 | 46 ++++++++++++++++++++++++++++++-------
 1 file changed, 38 insertions(+), 8 deletions(-)

diff --git a/src/physics/cam/hb_diff.F90 b/src/physics/cam/hb_diff.F90
index 306caa95e8..87796ecfd0 100644
--- a/src/physics/cam/hb_diff.F90
+++ b/src/physics/cam/hb_diff.F90
@@ -150,7 +150,7 @@ subroutine compute_hb_diff(ncol                  , &
     ! 
     !-----------------------------------------------------------------------
 
-    use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_ustar, calc_obukhov_length, austausch_atm, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, calc_kinematic_buoyancy_flux
+    use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length, austausch_atm, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, calc_kinematic_buoyancy_flux
     use physconst,            only: zvir, rair, gravit, karman
 
     !------------------------------Arguments--------------------------------
@@ -204,6 +204,7 @@ 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
     !
@@ -213,7 +214,7 @@ subroutine compute_hb_diff(ncol                  , &
 
     ! Compute ustar, Obukhov length, and kinematic surface fluxes.
     rrho(:ncol)   = calc_rrho(rair, t(:ncol,pver), pmid(:ncol,pver))
-    ustar(:ncol)  = calc_ustar(taux(:ncol),tauy(:ncol), rrho(:ncol))
+    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))
@@ -233,8 +234,22 @@ subroutine compute_hb_diff(ncol                  , &
     !
     ! Get atmosphere exchange coefficients
     !
-    kvf = austausch_atm(pcols, ncol, pver, ntop_turb, nbot_turb, &
-                        ml2, ri, s2)
+    kvf(:ncol,:) = 0.0_r8
+    do k = ntop_turb, nbot_turb-1
+      do i=1,ncol
+        kvf(i,k+1) = austausch_atm(ml2(k), ri(i, k), s2(i, k))
+      end do
+    end do
+   !  do i=1,ncol
+   !     kvf(i,ntop_turb+1:nbot_turb) = austausch_atm(ml2(   ntop_turb:nbot_turb-1),  &
+   !                                                   ri(i, ntop_turb:nbot_turb-1),  &
+   !                                                   s2(i, ntop_turb:nbot_turb-1))
+   !  end do
+   !  kvf(1:ncol,ntop_turb+1:nbot_turb) = austausch_atm(spread(ml2(        ntop_turb:nbot_turb-1), dim=1, ncopies=ncol), &
+   !                                                            ri(1:ncol, ntop_turb:nbot_turb-1),                       &
+   !                                                            s2(1:ncol, ntop_turb:nbot_turb-1))
+   !  kvf = austausch_atm(pcols, ncol, pver, ntop_turb, nbot_turb, &
+                        ! ml2, ri, s2)
     ! 
     ! Get pbl exchange coefficients
     !
@@ -264,7 +279,7 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     !
     !-----------------------------------------------------------------------
 
-    use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_ustar, calc_obukhov_length, austausch_atm_free, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, calc_kinematic_buoyancy_flux
+    use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length, austausch_atm_free, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, calc_kinematic_buoyancy_flux
     use physconst,            only: zvir, rair, gravit, karman
 
     !------------------------------Arguments--------------------------------
@@ -306,13 +321,14 @@ 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
     thv(:ncol,ntop_turb:) = calc_virtual_temperature(th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), zvir)
 
     ! Compute ustar, Obukhov length, and kinematic surface fluxes.
     rrho(:ncol)   = calc_rrho(rair, t(:ncol,pver), pmid(:ncol,pver))
-    ustar(:ncol)  = calc_ustar(taux(:ncol),tauy(:ncol), rrho(:ncol))
+    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))
@@ -324,8 +340,22 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     !
     ! Get free atmosphere exchange coefficients
     !
-    kvf = austausch_atm_free(pcols, ncol, pver, ntop_turb, nbot_turb, &
-                             ml2, ri, s2)
+    kvf(:ncol,:) = 0.0_r8
+    do k = ntop_turb, nbot_turb - 1
+      do i=1,ncol
+        kvf(i,k+1) = austausch_atm_free(ml2(k), ri(i, k), s2(i, k))
+      end do
+    end do
+   !  do i=1,ncol
+   !     kvf(i,ntop_turb+1:nbot_turb) = austausch_atm_free(ml2(ntop_turb:nbot_turb-1),     &
+   !                                                        ri(i, ntop_turb:nbot_turb-1),  &
+   !                                                        s2(i, ntop_turb:nbot_turb-1))
+   !  end do
+   !  kvf(1:ncol,ntop_turb+1:nbot_turb) = austausch_atm_free(spread(ml2(        ntop_turb:nbot_turb-1), dim=1, ncopies=ncol), &
+   !                                                                 ri(1:ncol, ntop_turb:nbot_turb-1),                       &
+   !                                                                 s2(1:ncol, ntop_turb:nbot_turb-1))
+   !  kvf = austausch_atm_free(pcols, ncol, pver, ntop_turb, nbot_turb, &
+                           !   ml2, ri, s2)
 
     kvq(:ncol,:) = kvf(:ncol,:)
     kvm(:ncol,:) = kvf(:ncol,:)

From f48e48f88ede852c9b69b062ec3ab868d603425a Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Tue, 21 Jan 2025 06:47:38 -0700
Subject: [PATCH 03/19] Updating atmos_phys to pbl utils branch and updating
 CIME ref to fix test bug.

---
 .gitmodules | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/.gitmodules b/.gitmodules
index 752c3bf4f7..b6e0371076 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -35,10 +35,10 @@
 
 [submodule "atmos_phys"]
 	path = src/atmos_phys
-	url = https://github.com/ESCOMP/atmospheric_physics
-        fxtag = atmos_phys0_07_001
+	url = https://github.com/mwaxmonsky/atmospheric_physics
+        fxtag = 71f5e041ee84a8a981a764064961b660789442ef
         fxrequired = AlwaysRequired
-	fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics
+	fxDONOTUSEurl = https://github.com/mwaxmonsky/atmospheric_physics
 
 [submodule "fv3"]
 	path = src/dynamics/fv3
@@ -151,7 +151,7 @@ fxDONOTUSEurl = https://github.com/ESMCI/ccs_config_cesm.git
 [submodule "cime"]
 path = cime
 url = https://github.com/ESMCI/cime
-fxtag = cime6.1.56
+fxtag = ef81824baa69947354f807e83924daf756e12937
 fxrequired = ToplevelRequired
 fxDONOTUSEurl = https://github.com/ESMCI/cime
 

From e9b851609f319ac89ace6b53f0c7d25ecb8aca39 Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Wed, 22 Jan 2025 14:14:38 -0700
Subject: [PATCH 04/19] Updating to latest atmos_phys fixes.  Adding atmos_phys
 phys_utils directory to build and integrating new api calls into hb_diff.F90

---
 .gitmodules                 |  2 +-
 bld/configure               |  1 +
 src/physics/cam/hb_diff.F90 | 33 +++++++--------------------------
 3 files changed, 9 insertions(+), 27 deletions(-)

diff --git a/.gitmodules b/.gitmodules
index b6e0371076..d404745936 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -36,7 +36,7 @@
 [submodule "atmos_phys"]
 	path = src/atmos_phys
 	url = https://github.com/mwaxmonsky/atmospheric_physics
-        fxtag = 71f5e041ee84a8a981a764064961b660789442ef
+        fxtag = 9193552576698433d612536985fc534594902f4c
         fxrequired = AlwaysRequired
 	fxDONOTUSEurl = https://github.com/mwaxmonsky/atmospheric_physics
 
diff --git a/bld/configure b/bld/configure
index 8f9f435c8c..9a537b837b 100755
--- a/bld/configure
+++ b/bld/configure
@@ -2336,6 +2336,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";
diff --git a/src/physics/cam/hb_diff.F90 b/src/physics/cam/hb_diff.F90
index 87796ecfd0..b5b509876d 100644
--- a/src/physics/cam/hb_diff.F90
+++ b/src/physics/cam/hb_diff.F90
@@ -235,21 +235,12 @@ subroutine compute_hb_diff(ncol                  , &
     ! Get atmosphere exchange coefficients
     !
     kvf(:ncol,:) = 0.0_r8
-    do k = ntop_turb, nbot_turb-1
-      do i=1,ncol
-        kvf(i,k+1) = austausch_atm(ml2(k), ri(i, k), s2(i, k))
-      end do
-    end do
-   !  do i=1,ncol
-   !     kvf(i,ntop_turb+1:nbot_turb) = austausch_atm(ml2(   ntop_turb:nbot_turb-1),  &
-   !                                                   ri(i, ntop_turb:nbot_turb-1),  &
-   !                                                   s2(i, ntop_turb:nbot_turb-1))
-   !  end do
-   !  kvf(1:ncol,ntop_turb+1:nbot_turb) = austausch_atm(spread(ml2(        ntop_turb:nbot_turb-1), dim=1, ncopies=ncol), &
-   !                                                            ri(1:ncol, ntop_turb:nbot_turb-1),                       &
-   !                                                            s2(1:ncol, ntop_turb:nbot_turb-1))
-   !  kvf = austausch_atm(pcols, ncol, pver, ntop_turb, nbot_turb, &
-                        ! ml2, ri, s2)
+   do k = ntop_turb, nbot_turb-1
+     do i = 1, ncol
+       kvf(i,k+1) = austausch_atm(ml2(k), ri(i, k), s2(i, k))
+     end do
+   end do
+
     ! 
     ! Get pbl exchange coefficients
     !
@@ -342,20 +333,10 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     !
     kvf(:ncol,:) = 0.0_r8
     do k = ntop_turb, nbot_turb - 1
-      do i=1,ncol
+      do i = 1, ncol
         kvf(i,k+1) = austausch_atm_free(ml2(k), ri(i, k), s2(i, k))
       end do
     end do
-   !  do i=1,ncol
-   !     kvf(i,ntop_turb+1:nbot_turb) = austausch_atm_free(ml2(ntop_turb:nbot_turb-1),     &
-   !                                                        ri(i, ntop_turb:nbot_turb-1),  &
-   !                                                        s2(i, ntop_turb:nbot_turb-1))
-   !  end do
-   !  kvf(1:ncol,ntop_turb+1:nbot_turb) = austausch_atm_free(spread(ml2(        ntop_turb:nbot_turb-1), dim=1, ncopies=ncol), &
-   !                                                                 ri(1:ncol, ntop_turb:nbot_turb-1),                       &
-   !                                                                 s2(1:ncol, ntop_turb:nbot_turb-1))
-   !  kvf = austausch_atm_free(pcols, ncol, pver, ntop_turb, nbot_turb, &
-                           !   ml2, ri, s2)
 
     kvq(:ncol,:) = kvf(:ncol,:)
     kvm(:ncol,:) = kvf(:ncol,:)

From e33da715c054acf5881fe71d8d4ece9f68608541 Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Thu, 23 Jan 2025 08:41:06 -0700
Subject: [PATCH 05/19] Initial PBL utils refactor for ccpp-ization.

---
 src/physics/cam/clubb_intr.F90         |  21 +-
 src/physics/cam/eddy_diff_cam.F90      |  19 +-
 src/physics/cam/pbl_utils.F90          | 330 -------------------------
 src/physics/cam/physpkg.F90            |   2 -
 src/physics/cam/vertical_diffusion.F90 |  26 +-
 src/physics/cam7/physpkg.F90           |   2 -
 6 files changed, 37 insertions(+), 363 deletions(-)

diff --git a/src/physics/cam/clubb_intr.F90 b/src/physics/cam/clubb_intr.F90
index a9f25f0256..2db987ab92 100644
--- a/src/physics/cam/clubb_intr.F90
+++ b/src/physics/cam/clubb_intr.F90
@@ -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_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
@@ -3493,8 +3494,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_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
@@ -4708,13 +4709,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) = (rga)*(state1%pdel(1:ncol,pver)/dz_g(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
diff --git a/src/physics/cam/eddy_diff_cam.F90 b/src/physics/cam/eddy_diff_cam.F90
index 1742bf5038..1316e72e3b 100644
--- a/src/physics/cam/eddy_diff_cam.F90
+++ b/src/physics/cam/eddy_diff_cam.F90
@@ -428,7 +428,7 @@ subroutine compute_eddy_diff( pbuf, lchnk  ,
   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 atmos_phys_pbl_utils, only: austausch_atm, calc_rrho, calc_friction_velocity
   use error_messages,   only: handle_errmsg
   use coords_1d,        only: Coords1D
   use wv_saturation,    only: qsat
@@ -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_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) = austausch_atm(ml2(k), ri(i, k), s2(i, k))
+           end do
+        end do
      else
         kvf = 0._r8
      end if
diff --git a/src/physics/cam/pbl_utils.F90 b/src/physics/cam/pbl_utils.F90
index 66759e295d..cd3af2ec60 100644
--- a/src/physics/cam/pbl_utils.F90
+++ b/src/physics/cam/pbl_utils.F90
@@ -22,206 +22,10 @@ module pbl_utils
 ! procedures, so they can accept scalars or any dimension of array as
 ! arguments, as long as all arguments have the same number of
 ! elements.
-public pbl_utils_init
-public calc_ustar
-public calc_obklen
-public virtem
 public compute_radf
-public austausch_atm, austausch_atm_free
-
-real(r8), parameter :: ustar_min = 0.01_r8
-
-real(r8) :: g         ! acceleration of gravity
-real(r8) :: vk        ! Von Karman's constant
-real(r8) :: cpair     ! specific heat of dry air
-real(r8) :: rair      ! gas constant for dry air
-real(r8) :: zvir      ! rh2o/rair - 1
-
-
-!------------------------------------------------------------------------!
-! Purpose: Compilers aren't creating optimized vector versions of        !
-!          elemental routines, so we'll explicitly create them and bind  !
-!          them via an interface for transparent use                     !
-!------------------------------------------------------------------------!
-interface calc_ustar
-  module procedure calc_ustar_scalar
-  module procedure calc_ustar_vector
-end interface 
-
-interface calc_obklen
-  module procedure calc_obklen_scalar
-  module procedure calc_obklen_vector
-end interface
-
-interface virtem
-  module procedure virtem_vector1D
-  module procedure virtem_vector2D  ! Used in hb_diff.F90
-end interface
-
-
 
 contains
 
-subroutine pbl_utils_init(g_in,vk_in,cpair_in,rair_in,zvir_in)
-
-  !-----------------------------------------------------------------------!
-  ! Purpose: Set constants to be used in calls to later functions         !
-  !-----------------------------------------------------------------------!
-
-  real(r8), intent(in) :: g_in       ! acceleration of gravity
-  real(r8), intent(in) :: vk_in      ! Von Karman's constant
-  real(r8), intent(in) :: cpair_in   ! specific heat of dry air
-  real(r8), intent(in) :: rair_in    ! gas constant for dry air
-  real(r8), intent(in) :: zvir_in    ! rh2o/rair - 1
-
-  g = g_in
-  vk = vk_in
-  cpair = cpair_in
-  rair = rair_in
-  zvir = zvir_in
-
-end subroutine pbl_utils_init
-
-subroutine calc_ustar_scalar( t,    pmid, taux, tauy, &
-                                 rrho, ustar)
-
-  !-----------------------------------------------------------------------!
-  ! Purpose: Calculate ustar and bottom level density (necessary for      !
-  !  Obukhov length calculation).                                         !
-  !-----------------------------------------------------------------------!
-
-  real(r8), intent(in) :: t         ! surface temperature
-  real(r8), intent(in) :: pmid      ! midpoint pressure (bottom level)
-  real(r8), intent(in) :: taux      ! surface u stress [N/m2]
-  real(r8), intent(in) :: tauy      ! surface v stress [N/m2]
-
-  real(r8), intent(out) :: rrho     ! 1./bottom level density
-  real(r8), intent(out) :: ustar    ! surface friction velocity [m/s]
-
-  rrho = rair * t / pmid
-  ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min )
-
-end subroutine calc_ustar_scalar
-
-subroutine calc_ustar_vector(n, t, pmid, taux, tauy, &
-                                 rrho, ustar)
-
-  !-----------------------------------------------------------------------!
-  ! Purpose: Calculate ustar and bottom level density (necessary for      !
-  !  Obukhov length calculation).                                         !
-  !-----------------------------------------------------------------------!
-  integer, intent(in) :: n             ! Length of vectors
-
-  real(r8), intent(in) :: t(n)         ! surface temperature
-  real(r8), intent(in) :: pmid(n)      ! midpoint pressure (bottom level)
-  real(r8), intent(in) :: taux(n)      ! surface u stress [N/m2]
-  real(r8), intent(in) :: tauy(n)      ! surface v stress [N/m2]
-
-
-  real(r8), intent(out) :: rrho(n)     ! 1./bottom level density
-  real(r8), intent(out) :: ustar(n)    ! surface friction velocity [m/s]
-
-
-  rrho = rair * t / pmid
-  ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min )
-
-end subroutine calc_ustar_vector
-
-subroutine calc_obklen_scalar( ths,  thvs, qflx, shflx, rrho, ustar, &
-                                  khfs, kqfs, kbfs, obklen)
-
-  !-----------------------------------------------------------------------!
-  ! Purpose: Calculate Obukhov length and kinematic fluxes.               !
-  !-----------------------------------------------------------------------!
-
-  real(r8), intent(in)  :: ths           ! potential temperature at surface [K]
-  real(r8), intent(in)  :: thvs          ! virtual potential temperature at surface
-  real(r8), intent(in)  :: qflx          ! water vapor flux (kg/m2/s)
-  real(r8), intent(in)  :: shflx         ! surface heat flux (W/m2)
-
-  real(r8), intent(in)  :: rrho          ! 1./bottom level density [ m3/kg ]
-  real(r8), intent(in)  :: ustar         ! Surface friction velocity [ m/s ]
-
-  real(r8), intent(out) :: khfs          ! sfc kinematic heat flux [mK/s]
-  real(r8), intent(out) :: kqfs          ! sfc kinematic water vapor flux [m/s]
-  real(r8), intent(out) :: kbfs          ! sfc kinematic buoyancy flux [m^2/s^3]
-  real(r8), intent(out) :: obklen        ! Obukhov length
-
-  ! Need kinematic fluxes for Obukhov:
-  khfs = shflx*rrho/cpair
-  kqfs = qflx*rrho
-  kbfs = khfs + zvir*ths*kqfs
-
-  ! Compute Obukhov length:
-  obklen = -thvs * ustar**3 / (g*vk*(kbfs + sign(1.e-10_r8,kbfs)))
-
-end subroutine calc_obklen_scalar
-
-subroutine calc_obklen_vector(n, ths,  thvs, qflx, shflx, rrho, ustar, &
-                                  khfs, kqfs, kbfs, obklen)
-
-  !-----------------------------------------------------------------------!
-  ! Purpose: Calculate Obukhov length and kinematic fluxes.               !
-  !-----------------------------------------------------------------------!
-  integer, intent(in) :: n                  ! Length of vectors
-
-  real(r8), intent(in)  :: ths(n)           ! potential temperature at surface [K]
-  real(r8), intent(in)  :: thvs(n)          ! virtual potential temperature at surface
-  real(r8), intent(in)  :: qflx(n)          ! water vapor flux (kg/m2/s)
-  real(r8), intent(in)  :: shflx(n)         ! surface heat flux (W/m2)
-
-  real(r8), intent(in)  :: rrho(n)          ! 1./bottom level density [ m3/kg ]
-  real(r8), intent(in)  :: ustar(n)         ! Surface friction velocity [ m/s ]
-
-  real(r8), intent(out) :: khfs(n)          ! sfc kinematic heat flux [mK/s]
-  real(r8), intent(out) :: kqfs(n)          ! sfc kinematic water vapor flux [m/s]
-  real(r8), intent(out) :: kbfs(n)          ! sfc kinematic buoyancy flux [m^2/s^3]
-  real(r8), intent(out) :: obklen(n)        ! Obukhov length
-
-
-  ! Need kinematic fluxes for Obukhov:
-  khfs = shflx*rrho/cpair
-  kqfs = qflx*rrho
-  kbfs = khfs + zvir*ths*kqfs
-
-  ! Compute Obukhov length:
-  obklen = -thvs * ustar**3 / (g*vk*(kbfs + sign(1.e-10_r8,kbfs)))
-
-end subroutine calc_obklen_vector
-
-subroutine virtem_vector1D(n, t,q, virtem)
-
-  !-----------------------------------------------------------------------!
-  ! Purpose: Calculate virtual temperature from temperature and specific  !
-  !  humidity.                                                            !
-  !-----------------------------------------------------------------------!
-
-  integer,  intent(in) :: n              ! vector length
-
-  real(r8), intent(in) :: t(n), q(n)
-  real(r8), intent(out):: virtem(n)
-
-  virtem = t * (1.0_r8 + zvir*q)
-
-end subroutine virtem_vector1D
-
-subroutine virtem_vector2D(n, m, t, q, virtem)
-
-  !-----------------------------------------------------------------------!
-  ! Purpose: Calculate virtual temperature from temperature and specific  !
-  !  humidity.                                                            !
-  !-----------------------------------------------------------------------!
-
-  integer,  intent(in) :: n, m            ! vector lengths
-
-  real(r8), intent(in) :: t(n,m), q(n,m)
-  real(r8), intent(out):: virtem(n,m)
-
-  virtem = t * (1.0_r8 + zvir*q)
-
-end subroutine virtem_vector2D
-
-
 subroutine compute_radf( choice_radf, i, pcols, pver, ncvmax, ncvfin, ktop, qmin, &
                          ql, pi, qrlw, g, cldeff, zi, chs, lwp_CL, opt_depth_CL,  &
                          radinvfrac_CL, radf_CL )
@@ -332,138 +136,4 @@ subroutine compute_radf( choice_radf, i, pcols, pver, ncvmax, ncvfin, ktop, qmin
   end do ! ncv = 1, ncvfin(i)
 end subroutine compute_radf
 
-subroutine austausch_atm(pcols, ncol, pver, ntop, nbot, ml2, ri, s2, kvf)
-
-  !---------------------------------------------------------------------- !
-  !                                                                       !
-  ! Purpose: Computes exchange coefficients for free turbulent flows.     !
-  !                                                                       !
-  ! Method:                                                               !
-  !                                                                       !
-  ! The free atmosphere diffusivities are based on standard mixing length !
-  ! forms for the neutral diffusivity multiplied by functns of Richardson !
-  ! number. K = l^2 * |dV/dz| * f(Ri). The same functions are used for    !
-  ! momentum, potential temperature, and constitutents.                   !
-  !                                                                       !
-  ! The stable Richardson num function (Ri>0) is taken from Holtslag and  !
-  ! Beljaars (1989), ECMWF proceedings. f = 1 / (1 + 10*Ri*(1 + 8*Ri))    !
-  ! The unstable Richardson number function (Ri<0) is taken from  CCM1.   !
-  ! f = sqrt(1 - 18*Ri)                                                   !
-  !                                                                       !
-  ! Author: B. Stevens (rewrite, August 2000)                             !
-  !                                                                       !
-  !---------------------------------------------------------------------- !
-
-  ! --------------- !
-  ! Input arguments !
-  ! --------------- !
-
-  integer,  intent(in)  :: pcols                ! Atmospheric columns dimension size
-  integer,  intent(in)  :: ncol                 ! Number of atmospheric columns
-  integer,  intent(in)  :: pver                 ! Number of atmospheric layers
-  integer,  intent(in)  :: ntop                 ! Top layer for calculation
-  integer,  intent(in)  :: nbot                 ! Bottom layer for calculation
-
-  real(r8), intent(in)  :: ml2(pver+1)          ! Mixing lengths squared
-  real(r8), intent(in)  :: s2(pcols,pver)       ! Shear squared
-  real(r8), intent(in)  :: ri(pcols,pver)       ! Richardson no
-
-  ! ---------------- !
-  ! Output arguments !
-  ! ---------------- !
-
-  real(r8), intent(out) :: kvf(pcols,pver+1)    ! Eddy diffusivity for heat and tracers
-
-  ! --------------- !
-  ! Local Variables !
-  ! --------------- !
-
-  real(r8)              :: fofri                ! f(ri)
-  real(r8)              :: kvn                  ! Neutral Kv
-
-  integer               :: i                    ! Longitude index
-  integer               :: k                    ! Vertical index
-
-  real(r8), parameter :: zkmin =  0.01_r8       ! Minimum kneutral*f(ri).
-
-  ! ----------------------- !
-  ! Main Computation Begins !
-  ! ----------------------- !
-
-  kvf(:ncol,:)           = 0.0_r8
-
-  ! Compute the free atmosphere vertical diffusion coefficients: kvh = kvq = kvm.
-
-  do k = ntop, nbot - 1
-     do i = 1, ncol
-        if( ri(i,k) < 0.0_r8 ) then
-           fofri = sqrt( max( 1._r8 - 18._r8 * ri(i,k), 0._r8 ) )
-        else
-           fofri = 1.0_r8 / ( 1.0_r8 + 10.0_r8 * ri(i,k) * ( 1.0_r8 + 8.0_r8 * ri(i,k) ) )
-        end if
-        kvn = ml2(k) * sqrt(s2(i,k))
-        kvf(i,k+1) = max( zkmin, kvn * fofri )
-     end do
-  end do
-
-end subroutine austausch_atm
-
-subroutine austausch_atm_free(pcols, ncol, pver, ntop, nbot, ml2, ri, s2, kvf)
-
-  !---------------------------------------------------------------------- !
-  !                                                                       !
-  ! same as austausch_atm but only mixing for Ri<0                        !
-  ! i.e. no background mixing and mixing for Ri>0                         !
-  !                                                                       !
-  !---------------------------------------------------------------------- !
-
-  ! --------------- !
-  ! Input arguments !
-  ! --------------- !
-
-  integer,  intent(in)  :: pcols                ! Atmospheric columns dimension size
-  integer,  intent(in)  :: ncol                 ! Number of atmospheric columns
-  integer,  intent(in)  :: pver                 ! Number of atmospheric layers
-  integer,  intent(in)  :: ntop                 ! Top layer for calculation
-  integer,  intent(in)  :: nbot                 ! Bottom layer for calculation
-
-  real(r8), intent(in)  :: ml2(pver+1)          ! Mixing lengths squared
-  real(r8), intent(in)  :: s2(pcols,pver)       ! Shear squared
-  real(r8), intent(in)  :: ri(pcols,pver)       ! Richardson no
-
-  ! ---------------- !
-  ! Output arguments !
-  ! ---------------- !
-
-  real(r8), intent(out) :: kvf(pcols,pver+1)    ! Eddy diffusivity for heat and tracers
-
-  ! --------------- !
-  ! Local Variables !
-  ! --------------- !
-
-  real(r8)              :: fofri                ! f(ri)
-  real(r8)              :: kvn                  ! Neutral Kv
-
-  integer               :: i                    ! Longitude index
-  integer               :: k                    ! Vertical index
-
-  ! ----------------------- !
-  ! Main Computation Begins !
-  ! ----------------------- !
-
-  kvf(:ncol,:)           = 0.0_r8
-  ! Compute the free atmosphere vertical diffusion coefficients: kvh = kvq = kvm.
-  do k = ntop, nbot - 1
-     do i = 1, ncol
-        if( ri(i,k) < 0.0_r8 ) then
-           fofri = sqrt( max( 1._r8 - 18._r8 * ri(i,k), 0._r8 ) )
-        else
-           fofri = 0.0_r8
-        end if
-        kvn = ml2(k) * sqrt(s2(i,k))
-        kvf(i,k+1) = kvn * fofri
-     end do
-  end do
-end subroutine austausch_atm_free
-
 end module pbl_utils
diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90
index 609006e75a..a204872a2a 100644
--- a/src/physics/cam/physpkg.F90
+++ b/src/physics/cam/physpkg.F90
@@ -746,7 +746,6 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )
     use tracers,            only: tracers_init
     use aoa_tracers,        only: aoa_tracers_init
     use rayleigh_friction,  only: rayleigh_friction_init
-    use pbl_utils,          only: pbl_utils_init
     use vertical_diffusion, only: vertical_diffusion_init
     use phys_debug_util,    only: phys_debug_init
     use rad_constituents,   only: rad_cnst_init
@@ -881,7 +880,6 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )
 
     call rayleigh_friction_init()
 
-    call pbl_utils_init(gravit, karman, cpair, rair, zvir)
     call vertical_diffusion_init(pbuf2d)
 
     if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then
diff --git a/src/physics/cam/vertical_diffusion.F90 b/src/physics/cam/vertical_diffusion.F90
index e2d571a4de..c9798534b2 100644
--- a/src/physics/cam/vertical_diffusion.F90
+++ b/src/physics/cam/vertical_diffusion.F90
@@ -727,7 +727,9 @@ subroutine vertical_diffusion_tend( &
   use constituents,       only : cnst_get_type_byind, cnst_name, &
        cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx
   use physconst,          only : pi
-  use pbl_utils,          only : virtem, calc_obklen, calc_ustar
+  use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_rrho, calc_friction_velocity, &
+                                  calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, calc_kinematic_buoyancy_flux, &
+                                  calc_obukhov_length
   use upper_bc,           only : ubc_get_vals, ubc_fixed_temp
   use upper_bc,           only : ubc_get_flxs
   use coords_1d,          only : Coords1D
@@ -1024,10 +1026,10 @@ subroutine vertical_diffusion_tend( &
 
      ! The diag_TKE scheme does not calculate the Monin-Obukhov length, which is used in dry deposition calculations.
      ! Use the routines from pbl_utils to accomplish this. Assumes ustar and rrho have been set.
-     call virtem(ncol, th(:ncol,pver),state%q(:ncol,pver,1), thvs(:ncol))
-     call calc_obklen(ncol, th(:ncol,pver), thvs(:ncol), cam_in%cflx(:ncol,1), &
-          cam_in%shf(:ncol), rrho(:ncol), ustar(:ncol), &
-          khfs(:ncol),    kqfs(:ncol), kbfs(:ncol),   obklen(:ncol))
+   !   call virtem(ncol, th(:ncol,pver),state%q(:ncol,pver,1), thvs(:ncol))
+   !   call calc_obklen(ncol, th(:ncol,pver), thvs(:ncol), cam_in%cflx(:ncol,1), &
+   !        cam_in%shf(:ncol), rrho(:ncol), ustar(:ncol), &
+   !        khfs(:ncol),    kqfs(:ncol), kbfs(:ncol),   obklen(:ncol))
 
 
   case ( 'HB', 'HBR', 'SPCAM_sam1mom' )
@@ -1080,14 +1082,14 @@ subroutine vertical_diffusion_tend( &
       ! is only handling other things, e.g. some boundary conditions, tms,
       ! and molecular diffusion.
 
-      call virtem(ncol, th(:ncol,pver),state%q(:ncol,pver,1), thvs(:ncol))
+      thvs  (:ncol) = calc_virtual_temperature(th(:ncol,pver), state%q(:ncol,pver,1), zvir)
+      rrho  (:ncol) = calc_rrho(rair, state%t(:ncol,pver), state%pmid(:ncol,pver))
+      ustar (:ncol) = calc_friction_velocity(cam_in%wsx(:ncol), cam_in%wsy(:ncol), rrho(:ncol))
+      khfs  (:ncol) = calc_kinematic_heat_flux(cam_in%shf(:ncol), rrho(:ncol), cpair)
+      kqfs  (:ncol) = calc_kinematic_water_vapor_flux(cam_in%cflx(:ncol,1), rrho(:ncol))
+      kbfs  (:ncol) = calc_kinematic_buoyancy_flux(khfs(:ncol), zvir, th(:ncol,pver), kqfs(:ncol))
+      obklen(:ncol) = calc_obukhov_length(thvs(:ncol), ustar(:ncol), gravit, karman, kbfs(:ncol))
 
-      call calc_ustar( ncol, state%t(:ncol,pver), state%pmid(:ncol,pver), &
-           cam_in%wsx(:ncol), cam_in%wsy(:ncol), rrho(:ncol), ustar(:ncol))
-      ! Use actual qflux, not lhf/latvap as was done previously
-      call calc_obklen( ncol, th(:ncol,pver), thvs(:ncol), cam_in%cflx(:ncol,1), &
-           cam_in%shf(:ncol), rrho(:ncol), ustar(:ncol),  &
-           khfs(:ncol), kqfs(:ncol), kbfs(:ncol), obklen(:ncol))
       ! These tendencies all applied elsewhere.
       kvm = 0._r8
       kvh = 0._r8
diff --git a/src/physics/cam7/physpkg.F90 b/src/physics/cam7/physpkg.F90
index e58f6399b8..54de7c3cdc 100644
--- a/src/physics/cam7/physpkg.F90
+++ b/src/physics/cam7/physpkg.F90
@@ -744,7 +744,6 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )
     use tracers,            only: tracers_init
     use aoa_tracers,        only: aoa_tracers_init
     use rayleigh_friction,  only: rayleigh_friction_init
-    use pbl_utils,          only: pbl_utils_init
     use vertical_diffusion, only: vertical_diffusion_init
     use phys_debug_util,    only: phys_debug_init
     use phys_debug,         only: phys_debug_state_init
@@ -879,7 +878,6 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out )
 
     call rayleigh_friction_init()
 
-    call pbl_utils_init(gravit, karman, cpair, rair, zvir)
     call vertical_diffusion_init(pbuf2d)
 
     if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then

From e295c91b43f213715fc3153e9ab7bfe37afb83c2 Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Thu, 23 Jan 2025 11:01:54 -0700
Subject: [PATCH 06/19] Remove dead code and add missing rrho function that was
 in calc_ustar.

---
 src/physics/cam/clubb_intr.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/physics/cam/clubb_intr.F90 b/src/physics/cam/clubb_intr.F90
index 2db987ab92..5878af1cbb 100644
--- a/src/physics/cam/clubb_intr.F90
+++ b/src/physics/cam/clubb_intr.F90
@@ -4709,7 +4709,7 @@ 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))
+    rrho   (1:ncol) = calc_rrho(rair, state1%t(i,pver), state1%pmid(i,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
     kinheat(1:ncol) = calc_kinematic_heat_flux(cam_in%shf(1:ncol), rrho(1:ncol), cpair)

From 6987d6a43945238a0384a0428c4b8a648d07ac17 Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Thu, 23 Jan 2025 11:02:32 -0700
Subject: [PATCH 07/19] Re-implement virtual temperature and obukhov length to
 fix diag_TKE tests.

---
 src/physics/cam/vertical_diffusion.F90 | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/physics/cam/vertical_diffusion.F90 b/src/physics/cam/vertical_diffusion.F90
index c9798534b2..4787b585ac 100644
--- a/src/physics/cam/vertical_diffusion.F90
+++ b/src/physics/cam/vertical_diffusion.F90
@@ -1026,11 +1026,12 @@ subroutine vertical_diffusion_tend( &
 
      ! The diag_TKE scheme does not calculate the Monin-Obukhov length, which is used in dry deposition calculations.
      ! Use the routines from pbl_utils to accomplish this. Assumes ustar and rrho have been set.
-   !   call virtem(ncol, th(:ncol,pver),state%q(:ncol,pver,1), thvs(:ncol))
-   !   call calc_obklen(ncol, th(:ncol,pver), thvs(:ncol), cam_in%cflx(:ncol,1), &
-   !        cam_in%shf(:ncol), rrho(:ncol), ustar(:ncol), &
-   !        khfs(:ncol),    kqfs(:ncol), kbfs(:ncol),   obklen(:ncol))
+      thvs  (:ncol) = calc_virtual_temperature(th(:ncol,pver), state%q(:ncol,pver,1), zvir)
 
+      khfs  (:ncol) = calc_kinematic_heat_flux(cam_in%shf(:ncol), rrho(:ncol), cpair)
+      kqfs  (:ncol) = calc_kinematic_water_vapor_flux(cam_in%cflx(:ncol,1), rrho(:ncol))
+      kbfs  (:ncol) = calc_kinematic_buoyancy_flux(khfs(:ncol), zvir, th(:ncol,pver), kqfs(:ncol))
+      obklen(:ncol) = calc_obukhov_length(thvs(:ncol), ustar(:ncol), gravit, karman, kbfs(:ncol))
 
   case ( 'HB', 'HBR', 'SPCAM_sam1mom' )
 

From b25f801430be40bf4f7a1986e75923ef1b9770be Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Thu, 23 Jan 2025 14:03:02 -0700
Subject: [PATCH 08/19] Fix incorrect parameters.

---
 src/physics/cam/clubb_intr.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/physics/cam/clubb_intr.F90 b/src/physics/cam/clubb_intr.F90
index 5878af1cbb..ef49915830 100644
--- a/src/physics/cam/clubb_intr.F90
+++ b/src/physics/cam/clubb_intr.F90
@@ -4709,7 +4709,7 @@ subroutine clubb_tend_cam( state,   ptend_all,   pbuf,     hdtime, &
     enddo
 
     ! diagnose surface friction and obukhov length (inputs to diagnose PBL depth)
-    rrho   (1:ncol) = calc_rrho(rair, state1%t(i,pver), state1%pmid(i,pver))
+    rrho   (1:ncol) = calc_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
     kinheat(1:ncol) = calc_kinematic_heat_flux(cam_in%shf(1:ncol), rrho(1:ncol), cpair)

From 73846518ae7470f70e1a6615054476fac54bbc5c Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Thu, 23 Jan 2025 16:38:15 -0700
Subject: [PATCH 09/19] Updating to latest atmos_phys branch.

---
 .gitmodules | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.gitmodules b/.gitmodules
index d404745936..a3e22f352a 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -36,7 +36,7 @@
 [submodule "atmos_phys"]
 	path = src/atmos_phys
 	url = https://github.com/mwaxmonsky/atmospheric_physics
-        fxtag = 9193552576698433d612536985fc534594902f4c
+        fxtag = bcc145ad9cc06185e58054ec01b008c2d9b9fa67
         fxrequired = AlwaysRequired
 	fxDONOTUSEurl = https://github.com/mwaxmonsky/atmospheric_physics
 

From ca5df2fcbad6e0d42a3031d456fc59db27b6a10e Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Thu, 23 Jan 2025 20:44:10 -0700
Subject: [PATCH 10/19] Updating loop whitespace.

---
 src/physics/cam/hb_diff.F90 | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/physics/cam/hb_diff.F90 b/src/physics/cam/hb_diff.F90
index b5b509876d..e97e30d0f5 100644
--- a/src/physics/cam/hb_diff.F90
+++ b/src/physics/cam/hb_diff.F90
@@ -235,11 +235,11 @@ subroutine compute_hb_diff(ncol                  , &
     ! Get atmosphere exchange coefficients
     !
     kvf(:ncol,:) = 0.0_r8
-   do k = ntop_turb, nbot_turb-1
-     do i = 1, ncol
-       kvf(i,k+1) = austausch_atm(ml2(k), ri(i, k), s2(i, k))
-     end do
-   end do
+    do k = ntop_turb, nbot_turb-1
+       do i = 1, ncol
+          kvf(i,k+1) = austausch_atm(ml2(k), ri(i, k), s2(i, k))
+       end do
+    end do
 
     ! 
     ! Get pbl exchange coefficients

From eedf005aec179b463c3c0fd2900819f711c09458 Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Thu, 23 Jan 2025 20:44:46 -0700
Subject: [PATCH 11/19] Updating loop whitespace.

---
 src/physics/cam/hb_diff.F90 | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/physics/cam/hb_diff.F90 b/src/physics/cam/hb_diff.F90
index e97e30d0f5..9d7917f904 100644
--- a/src/physics/cam/hb_diff.F90
+++ b/src/physics/cam/hb_diff.F90
@@ -333,9 +333,9 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     !
     kvf(:ncol,:) = 0.0_r8
     do k = ntop_turb, nbot_turb - 1
-      do i = 1, ncol
-        kvf(i,k+1) = austausch_atm_free(ml2(k), ri(i, k), s2(i, k))
-      end do
+       do i = 1, ncol
+          kvf(i,k+1) = austausch_atm_free(ml2(k), ri(i, k), s2(i, k))
+       end do
     end do
 
     kvq(:ncol,:) = kvf(:ncol,:)

From 049f6575bebe77ace1c97e15447b0d1c4d5fb767 Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Thu, 23 Jan 2025 21:39:22 -0700
Subject: [PATCH 12/19] Spacing updates.

---
 src/physics/cam/hb_diff.F90            |  8 +++-
 src/physics/cam/vertical_diffusion.F90 | 52 +++++++++++++-------------
 2 files changed, 32 insertions(+), 28 deletions(-)

diff --git a/src/physics/cam/hb_diff.F90 b/src/physics/cam/hb_diff.F90
index 9d7917f904..20658de54b 100644
--- a/src/physics/cam/hb_diff.F90
+++ b/src/physics/cam/hb_diff.F90
@@ -150,7 +150,9 @@ subroutine compute_hb_diff(ncol                  , &
     ! 
     !-----------------------------------------------------------------------
 
-    use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length, austausch_atm, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, calc_kinematic_buoyancy_flux
+    use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length,               &
+                                    austausch_atm, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, &
+                                    calc_kinematic_buoyancy_flux
     use physconst,            only: zvir, rair, gravit, karman
 
     !------------------------------Arguments--------------------------------
@@ -270,7 +272,9 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     !
     !-----------------------------------------------------------------------
 
-    use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length, austausch_atm_free, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, calc_kinematic_buoyancy_flux
+    use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length,                    &
+                                    austausch_atm_free, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, &
+                                    calc_kinematic_buoyancy_flux
     use physconst,            only: zvir, rair, gravit, karman
 
     !------------------------------Arguments--------------------------------
diff --git a/src/physics/cam/vertical_diffusion.F90 b/src/physics/cam/vertical_diffusion.F90
index 4787b585ac..be1b5be586 100644
--- a/src/physics/cam/vertical_diffusion.F90
+++ b/src/physics/cam/vertical_diffusion.F90
@@ -707,34 +707,34 @@ subroutine vertical_diffusion_tend( &
   !---------------------------------------------------- !
   ! This is an interface routine for vertical diffusion !
   !---------------------------------------------------- !
-  use physics_buffer,     only : physics_buffer_desc, pbuf_get_field, pbuf_set_field
-  use physics_types,      only : physics_state, physics_ptend, physics_ptend_init
-  use physics_types,      only : set_dry_to_wet, set_wet_to_dry
-
-  use camsrfexch,         only : cam_in_t
-  use cam_history,        only : outfld
-
-  use trb_mtn_stress_cam, only : trb_mtn_stress_tend
-  use beljaars_drag_cam,  only : beljaars_drag_tend
-  use eddy_diff_cam,      only : eddy_diff_tend
-  use hb_diff,            only : compute_hb_diff, compute_hb_free_atm_diff
-  use wv_saturation,      only : qsat
-  use molec_diff,         only : compute_molec_diff, vd_lu_qdecomp
-  use constituents,       only : qmincg, qmin, cnst_type
-  use diffusion_solver,   only : compute_vdiff, any, operator(.not.)
-  use air_composition,    only : cpairv, rairv !Needed for calculation of upward H flux
-  use time_manager,       only : get_nstep
-  use constituents,       only : cnst_get_type_byind, cnst_name, &
-       cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx
-  use physconst,          only : pi
-  use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_rrho, calc_friction_velocity, &
+  use physics_buffer,       only : physics_buffer_desc, pbuf_get_field, pbuf_set_field
+  use physics_types,        only : physics_state, physics_ptend, physics_ptend_init
+  use physics_types,        only : set_dry_to_wet, set_wet_to_dry
+
+  use camsrfexch,           only : cam_in_t
+  use cam_history,          only : outfld
+
+  use trb_mtn_stress_cam,   only : trb_mtn_stress_tend
+  use beljaars_drag_cam,    only : beljaars_drag_tend
+  use eddy_diff_cam,        only : eddy_diff_tend
+  use hb_diff,              only : compute_hb_diff, compute_hb_free_atm_diff
+  use wv_saturation,        only : qsat
+  use molec_diff,           only : compute_molec_diff, vd_lu_qdecomp
+  use constituents,         only : qmincg, qmin, cnst_type
+  use diffusion_solver,     only : compute_vdiff, any, operator(.not.)
+  use air_composition,      only : cpairv, rairv !Needed for calculation of upward H flux
+  use time_manager,         only : get_nstep
+  use constituents,         only : cnst_get_type_byind, cnst_name, &
+                                   cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx
+  use physconst,            only : pi
+  use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_rrho, calc_friction_velocity,                             &
                                   calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, calc_kinematic_buoyancy_flux, &
                                   calc_obukhov_length
-  use upper_bc,           only : ubc_get_vals, ubc_fixed_temp
-  use upper_bc,           only : ubc_get_flxs
-  use coords_1d,          only : Coords1D
-  use phys_control,       only : cam_physpkg_is
-  use ref_pres,           only : ptop_ref
+  use upper_bc,             only : ubc_get_vals, ubc_fixed_temp
+  use upper_bc,             only : ubc_get_flxs
+  use coords_1d,            only : Coords1D
+  use phys_control,         only : cam_physpkg_is
+  use ref_pres,             only : ptop_ref
 
   ! --------------- !
   ! Input Arguments !

From c551de7b06c8efc6355bb820e3eefb4e64597a1d Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Thu, 23 Jan 2025 21:41:45 -0700
Subject: [PATCH 13/19] Spacing update.

---
 src/physics/cam/eddy_diff_cam.F90 | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/src/physics/cam/eddy_diff_cam.F90 b/src/physics/cam/eddy_diff_cam.F90
index 1316e72e3b..cd6b9fb959 100644
--- a/src/physics/cam/eddy_diff_cam.F90
+++ b/src/physics/cam/eddy_diff_cam.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 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: austausch_atm, calc_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
+  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 !

From 445bcbe5edf54a0e28b2b1417c6dbbef6c950d0f Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Tue, 28 Jan 2025 14:39:57 -0700
Subject: [PATCH 14/19] Updating austausch_atm references to new function
 names.

---
 .gitmodules                       | 2 +-
 src/atmos_phys                    | 2 +-
 src/physics/cam/eddy_diff_cam.F90 | 4 ++--
 src/physics/cam/hb_diff.F90       | 9 ++++-----
 4 files changed, 8 insertions(+), 9 deletions(-)

diff --git a/.gitmodules b/.gitmodules
index a3e22f352a..a7135f27b3 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -36,7 +36,7 @@
 [submodule "atmos_phys"]
 	path = src/atmos_phys
 	url = https://github.com/mwaxmonsky/atmospheric_physics
-        fxtag = bcc145ad9cc06185e58054ec01b008c2d9b9fa67
+        fxtag = 330262591a964d09e4d22bad462e767f0950b2dc
         fxrequired = AlwaysRequired
 	fxDONOTUSEurl = https://github.com/mwaxmonsky/atmospheric_physics
 
diff --git a/src/atmos_phys b/src/atmos_phys
index c3de8468f7..330262591a 160000
--- a/src/atmos_phys
+++ b/src/atmos_phys
@@ -1 +1 @@
-Subproject commit c3de8468f7b245a939448f4ca6d3ef386584e92d
+Subproject commit 330262591a964d09e4d22bad462e767f0950b2dc
diff --git a/src/physics/cam/eddy_diff_cam.F90 b/src/physics/cam/eddy_diff_cam.F90
index cd6b9fb959..c3ab0653ba 100644
--- a/src/physics/cam/eddy_diff_cam.F90
+++ b/src/physics/cam/eddy_diff_cam.F90
@@ -428,7 +428,7 @@ subroutine compute_eddy_diff( pbuf, lchnk  ,
   use cam_history,          only: outfld
   use phys_debug_util,      only: phys_debug_col
   use air_composition,      only: cpairv
-  use atmos_phys_pbl_utils, only: austausch_atm, calc_rrho, calc_friction_velocity
+  use atmos_phys_pbl_utils, only: calc_eddy_flux_coefficient, calc_rrho, calc_friction_velocity
   use error_messages,       only: handle_errmsg
   use coords_1d,            only: Coords1D
   use wv_saturation,        only: qsat
@@ -698,7 +698,7 @@ subroutine compute_eddy_diff( pbuf, lchnk  ,
         kvf(:ncol,:) = 0.0_r8
         do k = ntop_eddy, nbot_eddy-1
            do i = 1, ncol
-              kvf(i,k+1) = austausch_atm(ml2(k), ri(i, k), s2(i, k))
+              kvf(i,k+1) = calc_eddy_flux_coefficient(ml2(k), ri(i, k), s2(i, k))
            end do
         end do
      else
diff --git a/src/physics/cam/hb_diff.F90 b/src/physics/cam/hb_diff.F90
index 20658de54b..300cc873a1 100644
--- a/src/physics/cam/hb_diff.F90
+++ b/src/physics/cam/hb_diff.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--------------------------------
@@ -151,7 +150,7 @@ subroutine compute_hb_diff(ncol                  , &
     !-----------------------------------------------------------------------
 
     use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length,               &
-                                    austausch_atm, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, &
+                                    calc_eddy_flux_coefficient, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, &
                                     calc_kinematic_buoyancy_flux
     use physconst,            only: zvir, rair, gravit, karman
 
@@ -239,7 +238,7 @@ subroutine compute_hb_diff(ncol                  , &
     kvf(:ncol,:) = 0.0_r8
     do k = ntop_turb, nbot_turb-1
        do i = 1, ncol
-          kvf(i,k+1) = austausch_atm(ml2(k), ri(i, k), s2(i, k))
+          kvf(i,k+1) = calc_eddy_flux_coefficient(ml2(k), ri(i, k), s2(i, k))
        end do
     end do
 
@@ -273,7 +272,7 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     !-----------------------------------------------------------------------
 
     use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length,                    &
-                                    austausch_atm_free, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, &
+                                    calc_free_atm_eddy_flux_coefficient, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, &
                                     calc_kinematic_buoyancy_flux
     use physconst,            only: zvir, rair, gravit, karman
 
@@ -338,7 +337,7 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     kvf(:ncol,:) = 0.0_r8
     do k = ntop_turb, nbot_turb - 1
        do i = 1, ncol
-          kvf(i,k+1) = austausch_atm_free(ml2(k), ri(i, k), s2(i, k))
+          kvf(i,k+1) = calc_free_atm_eddy_flux_coefficient(ml2(k), ri(i, k), s2(i, k))
        end do
     end do
 

From e8e5d55d197a81dd5a023f3cdbbb437c89f9111a Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Tue, 28 Jan 2025 14:56:13 -0700
Subject: [PATCH 15/19] Merging in latest changes from atmos_phys.

---
 .gitmodules    | 2 +-
 src/atmos_phys | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/.gitmodules b/.gitmodules
index a7135f27b3..09c2c67ba0 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -36,7 +36,7 @@
 [submodule "atmos_phys"]
 	path = src/atmos_phys
 	url = https://github.com/mwaxmonsky/atmospheric_physics
-        fxtag = 330262591a964d09e4d22bad462e767f0950b2dc
+        fxtag = 426f759e1a49d1e0d10e7476e4bf69b02da812a6
         fxrequired = AlwaysRequired
 	fxDONOTUSEurl = https://github.com/mwaxmonsky/atmospheric_physics
 
diff --git a/src/atmos_phys b/src/atmos_phys
index 330262591a..426f759e1a 160000
--- a/src/atmos_phys
+++ b/src/atmos_phys
@@ -1 +1 @@
-Subproject commit 330262591a964d09e4d22bad462e767f0950b2dc
+Subproject commit 426f759e1a49d1e0d10e7476e4bf69b02da812a6

From e776da57f8a8fe5c43bf15050702fa66af5ecb4e Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Mon, 24 Feb 2025 09:42:48 -0700
Subject: [PATCH 16/19] Updating atmos_phys to latest tag.

---
 .gitmodules    | 2 +-
 src/atmos_phys | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/.gitmodules b/.gitmodules
index 8c3dd6f72c..b0765c5b8b 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -36,7 +36,7 @@
 [submodule "atmos_phys"]
 	path = src/atmos_phys
 	url = https://github.com/mwaxmonsky/atmospheric_physics
-        fxtag = 426f759e1a49d1e0d10e7476e4bf69b02da812a6
+        fxtag = 7031edf10a3b2b63b92573cf03277ce8d65c073e
         fxrequired = AlwaysRequired
 	fxDONOTUSEurl = https://github.com/mwaxmonsky/atmospheric_physics
 
diff --git a/src/atmos_phys b/src/atmos_phys
index 426f759e1a..7031edf10a 160000
--- a/src/atmos_phys
+++ b/src/atmos_phys
@@ -1 +1 @@
-Subproject commit 426f759e1a49d1e0d10e7476e4bf69b02da812a6
+Subproject commit 7031edf10a3b2b63b92573cf03277ce8d65c073e

From 29f776818398383c68456adcedaecb54d5fb03ce Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Mon, 24 Feb 2025 09:59:17 -0700
Subject: [PATCH 17/19] Converting reference from 'calc_rrho' to
 'calc_ideal_gas_rrho' from latest atmos_phys.

---
 src/physics/cam/clubb_intr.F90         | 6 +++---
 src/physics/cam/eddy_diff_cam.F90      | 4 ++--
 src/physics/cam/hb_diff.F90            | 8 ++++----
 src/physics/cam/vertical_diffusion.F90 | 4 ++--
 4 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/src/physics/cam/clubb_intr.F90 b/src/physics/cam/clubb_intr.F90
index 7e5a11d7b8..ccd655d2ab 100644
--- a/src/physics/cam/clubb_intr.F90
+++ b/src/physics/cam/clubb_intr.F90
@@ -26,7 +26,7 @@ module clubb_intr
 
   use spmd_utils,          only: masterproc
   use constituents,        only: pcnst, cnst_add
-  use atmos_phys_pbl_utils,only: calc_friction_velocity, calc_kinematic_heat_flux, calc_rrho, &
+  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
 
@@ -3492,7 +3492,7 @@ 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
 
-        rrho(i) = calc_rrho(rair, state1%t(i,pver), state1%pmid(i,pver))
+        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
@@ -4707,7 +4707,7 @@ subroutine clubb_tend_cam( state,   ptend_all,   pbuf,     hdtime, &
     enddo
 
     ! diagnose surface friction and obukhov length (inputs to diagnose PBL depth)
-    rrho   (1:ncol) = calc_rrho(rair, state1%t(1:ncol,pver), state1%pmid(1:ncol,pver))
+    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
     kinheat(1:ncol) = calc_kinematic_heat_flux(cam_in%shf(1:ncol), rrho(1:ncol), cpair)
diff --git a/src/physics/cam/eddy_diff_cam.F90 b/src/physics/cam/eddy_diff_cam.F90
index c3ab0653ba..05ff18d9cd 100644
--- a/src/physics/cam/eddy_diff_cam.F90
+++ b/src/physics/cam/eddy_diff_cam.F90
@@ -428,7 +428,7 @@ subroutine compute_eddy_diff( pbuf, lchnk  ,
   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_rrho, calc_friction_velocity
+  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
@@ -670,7 +670,7 @@ subroutine compute_eddy_diff( pbuf, lchnk  ,
      ! I am using updated wind, here.
 
      ! Compute ustar
-     rrho(:ncol)   = calc_rrho(rair, tfd(:ncol,pver), pmid(:ncol,pver))
+     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))
diff --git a/src/physics/cam/hb_diff.F90 b/src/physics/cam/hb_diff.F90
index 300cc873a1..81ad8ff7bf 100644
--- a/src/physics/cam/hb_diff.F90
+++ b/src/physics/cam/hb_diff.F90
@@ -150,7 +150,7 @@ subroutine compute_hb_diff(ncol                  , &
     !-----------------------------------------------------------------------
 
     use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length,               &
-                                    calc_eddy_flux_coefficient, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, &
+                                    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
 
@@ -214,7 +214,7 @@ subroutine compute_hb_diff(ncol                  , &
     thv(:ncol,ntop_turb:) = calc_virtual_temperature(th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), zvir)
 
     ! Compute ustar, Obukhov length, and kinematic surface fluxes.
-    rrho(:ncol)   = calc_rrho(rair, t(:ncol,pver), pmid(:ncol,pver))
+    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))
@@ -272,7 +272,7 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     !-----------------------------------------------------------------------
 
     use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length,                    &
-                                    calc_free_atm_eddy_flux_coefficient, calc_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, &
+                                    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
 
@@ -321,7 +321,7 @@ subroutine compute_hb_free_atm_diff(ncol,          &
     thv(:ncol,ntop_turb:) = calc_virtual_temperature(th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), zvir)
 
     ! Compute ustar, Obukhov length, and kinematic surface fluxes.
-    rrho(:ncol)   = calc_rrho(rair, t(:ncol,pver), pmid(:ncol,pver))
+    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))
diff --git a/src/physics/cam/vertical_diffusion.F90 b/src/physics/cam/vertical_diffusion.F90
index 1f03592172..d44f952559 100644
--- a/src/physics/cam/vertical_diffusion.F90
+++ b/src/physics/cam/vertical_diffusion.F90
@@ -726,7 +726,7 @@ subroutine vertical_diffusion_tend( &
   use constituents,         only : cnst_get_type_byind, cnst_name, &
                                    cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx
   use physconst,            only : pi
-  use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_rrho, calc_friction_velocity,                             &
+  use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_ideal_gas_rrho, calc_friction_velocity,                             &
                                   calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, calc_kinematic_buoyancy_flux, &
                                   calc_obukhov_length
   use upper_bc,             only : ubc_get_vals, ubc_fixed_temp
@@ -1080,7 +1080,7 @@ subroutine vertical_diffusion_tend( &
       ! and molecular diffusion.
 
       thvs  (:ncol) = calc_virtual_temperature(th(:ncol,pver), state%q(:ncol,pver,1), zvir)
-      rrho  (:ncol) = calc_rrho(rair, state%t(:ncol,pver), state%pmid(:ncol,pver))
+      rrho  (:ncol) = calc_ideal_gas_rrho(rair, state%t(:ncol,pver), state%pmid(:ncol,pver))
       ustar (:ncol) = calc_friction_velocity(cam_in%wsx(:ncol), cam_in%wsy(:ncol), rrho(:ncol))
       khfs  (:ncol) = calc_kinematic_heat_flux(cam_in%shf(:ncol), rrho(:ncol), cpair)
       kqfs  (:ncol) = calc_kinematic_water_vapor_flux(cam_in%cflx(:ncol,1), rrho(:ncol))

From bb7e7d5b1685f90517e0fa0bdff59e9d9c8a04be Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Mon, 24 Feb 2025 13:18:22 -0700
Subject: [PATCH 18/19] Updating atmos_phys to latest tag.

---
 .gitmodules | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/.gitmodules b/.gitmodules
index b0765c5b8b..6c2dbe2e68 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -35,10 +35,10 @@
 
 [submodule "atmos_phys"]
 	path = src/atmos_phys
-	url = https://github.com/mwaxmonsky/atmospheric_physics
-        fxtag = 7031edf10a3b2b63b92573cf03277ce8d65c073e
+	url = https://github.com/ESCOMP/atmospheric_physics
+        fxtag = atmos_phys0_09_000
         fxrequired = AlwaysRequired
-	fxDONOTUSEurl = https://github.com/mwaxmonsky/atmospheric_physics
+	fxDONOTUSEurl = https://github.com/ESCOMP/atmospheric_physics
 
 [submodule "fv3"]
 	path = src/dynamics/fv3

From a1774df8b2c83c5a6d3549a91b09f483c726d5b0 Mon Sep 17 00:00:00 2001
From: Michael Waxmonsky <mwaxmonsky@ucar.edu>
Date: Wed, 26 Feb 2025 15:41:37 -0700
Subject: [PATCH 19/19] Updating ChangeLog.

---
 doc/ChangeLog | 75 +++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 75 insertions(+)

diff --git a/doc/ChangeLog b/doc/ChangeLog
index e57931b4cd..84e989f6de 100644
--- a/doc/ChangeLog
+++ b/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