PBL utils transfer (#193)

Originator(s): @mwaxmonsky Summary (include the keyword ['closes', 'fixes', 'resolves'] and issue number): Adds pbl_utils capability from CAM. Describe any changes made to the namelist: N/A List all files eliminated and why: N/A List all files added and what they do: M phys_utils/atmos_phys_pbl_utils.F90 - Transferred all functionality (except for `compute_radf(...)`) from ESCOMP/cam/src/physics/cam/pbl_utils.F90 M .github/workflows/unit-tests.yaml - Added phys_utils to code coverage. M test/unit-test/CMakeLists.txt M test/unit-test/tests/CMakeLists.txt M test/unit-test/tests/phys_utils/CMakeLists.txt M test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf - Creating CMake target library for `phys_utils` and adding sample unit tests. List all existing files that have been modified, and describe the changes: N/A (Helpful git command: `git diff --name-status development...<your_branch_name>`) List any test failures: N/A Is this a science-changing update? New physics package, algorithm change, tuning changes, etc? No --------- Co-authored-by: Jesse Nusbaumer <[email protected]>
ESCOMP · Feb 19, 2025 · 603de70 · 603de70
1 parent 252b500
commit 603de70
Show file tree

Hide file tree

Showing 6 changed files with 256 additions and 2 deletions.
diff --git a/.github/workflows/unit-tests.yaml b/.github/workflows/unit-tests.yaml
@@ -62,7 +62,7 @@ jobs:
         run: |
           source venv/bin/activate
           cd build
-          gcovr -r .. --filter '\.\./schemes' --html atmospheric_physics_code_coverage.html --txt
+          gcovr -r .. --filter '\.\./schemes' --filter '\.\./phys_utils' --html atmospheric_physics_code_coverage.html --txt
 
       - name: Upload code coverage results
         uses: actions/upload-artifact@v4

diff --git a/phys_utils/atmos_phys_pbl_utils.F90 b/phys_utils/atmos_phys_pbl_utils.F90
@@ -0,0 +1,216 @@
+module atmos_phys_pbl_utils
+    ! Planetary boundary layer related functions used for vertical diffusion schemes.
+
+    use ccpp_kinds, only: kind_phys
+
+    implicit none
+    private
+
+    public :: calc_ideal_gas_rrho
+    public :: calc_friction_velocity
+    public :: calc_kinematic_heat_flux
+    public :: calc_kinematic_water_vapor_flux
+    public :: calc_kinematic_buoyancy_flux
+    public :: calc_obukhov_length
+    public :: calc_virtual_temperature
+    public :: calc_eddy_flux_coefficient
+    public :: calc_free_atm_eddy_flux_coefficient
+
+    real(kind_phys), parameter :: minimum_friction_velocity     = 0.01_kind_phys ! Assuming minimum for coarse grids
+    real(kind_phys), parameter :: minimum_eddy_flux_coefficient = 0.01_kind_phys ! CCM1 2.f.14
+
+contains
+
+    pure elemental function calc_ideal_gas_rrho(rair, surface_temperature, pmid) result(rrho)
+        ! air density reciprocal
+        ! Taken from https://glossary.ametsoc.org/wiki/Equation_of_state
+        ! where \alpha = rrho
+        real(kind_phys), intent(in) :: rair  ! gas constant for dry air         [ J kg-1 K-1 ]
+        real(kind_phys), intent(in) :: surface_temperature                    ! [ K          ]
+        real(kind_phys), intent(in) :: pmid  ! midpoint pressure (bottom level) [ Pa         ]
+
+        real(kind_phys)             :: rrho  ! 1./bottom level density          [ m3 kg-1    ]
+
+        rrho = rair * surface_temperature / pmid
+    end function calc_ideal_gas_rrho
+
+    pure elemental function calc_friction_velocity(taux, tauy, rrho) result(friction_velocity)
+        ! https://glossary.ametsoc.org/wiki/Friction_velocity
+        ! NOTE: taux,tauy come form the expansion of the Reynolds stress
+        !
+        ! Also found in:
+        ! Stull, Roland B. An Introduction to Boundary Layer Meteorology. Springer Kluwer Academic Publishers, 1988. Print.
+        ! DOI: https://doi.org/10.1007/978-94-009-3027-8
+        ! Equation 2.10b, page 67
+
+        real(kind_phys), intent(in) :: taux              ! surface u stress          [ N m-2   ]
+        real(kind_phys), intent(in) :: tauy              ! surface v stress          [ N m-2   ]
+        real(kind_phys), intent(in) :: rrho              ! 1./bottom level density   [ m3 kg-1 ]
+
+        real(kind_phys)             :: friction_velocity ! surface friction velocity [ m s-1   ]
+
+        friction_velocity = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), minimum_friction_velocity )
+    end function calc_friction_velocity
+
+    pure elemental function calc_kinematic_heat_flux(shflx, rrho, cpair) result(khfs)
+        real(kind_phys), intent(in)  :: shflx  ! surface heat flux            [ W m-2       ]
+        real(kind_phys), intent(in)  :: rrho   ! 1./bottom level density      [ m3 kg-1     ]
+        real(kind_phys), intent(in)  :: cpair  ! specific heat of dry air     [ J kg-1 K-1  ]
+
+        real(kind_phys)              :: khfs   ! surface kinematic heat flux  [ m K s-1     ]
+
+        khfs = shflx*rrho/cpair
+    end function calc_kinematic_heat_flux
+
+    pure elemental function calc_kinematic_water_vapor_flux(qflx, rrho) result(kqfs)
+        real(kind_phys), intent(in)  :: qflx ! water vapor flux                   [ kg m-2 s-1    ]
+        real(kind_phys), intent(in)  :: rrho ! 1./bottom level density            [ m3 kg-1       ]
+
+        real(kind_phys)              :: kqfs ! surface kinematic water vapor flux [ kg kg-1 m s-1 ]
+
+        kqfs = qflx*rrho
+    end function calc_kinematic_water_vapor_flux
+
+    pure elemental function calc_kinematic_buoyancy_flux(khfs, zvir, ths, kqfs) result(kbfs)
+        real(kind_phys), intent(in) :: khfs  ! surface kinematic heat flux          [ m K s-1       ]
+        real(kind_phys), intent(in) :: zvir  ! rh2o/rair - 1
+        real(kind_phys), intent(in) :: ths   ! potential temperature at surface     [ K             ]
+        real(kind_phys), intent(in) :: kqfs  ! surface kinematic water vapor flux   [ kg kg-1 m s-1 ]
+
+        ! (`kbfs = \overline{(w' \theta'_v)}_s`)
+        real(kind_phys)             :: kbfs  ! surface kinematic buoyancy flux      [ m K s-1 ]
+
+        kbfs = khfs + zvir*ths*kqfs
+    end function calc_kinematic_buoyancy_flux
+
+    pure elemental function calc_obukhov_length(thvs, ustar, g, karman, kbfs) result(obukhov_length)
+        ! Stull, Roland B. An Introduction to Boundary Layer Meteorology. Springer Kluwer Academic Publishers, 1988. Print.
+        ! DOI: https://doi.org/10.1007/978-94-009-3027-8
+        ! Equation 5.7c, page 181
+        ! \frac{-\theta*u_*^3}{g*k*\overline{(w' \theta_v')}_s} = frac{-\theta*u_*^3}{g*k*kbfs}
+
+        real(kind_phys), intent(in)  :: thvs              ! virtual potential temperature at surface [ K       ]
+        real(kind_phys), intent(in)  :: ustar             ! Surface friction velocity                [ m s-1   ]
+        real(kind_phys), intent(in)  :: g                 ! acceleration of gravity                  [ m s-2   ]
+        real(kind_phys), intent(in)  :: karman            ! Von Karman's constant (unitless)
+        real(kind_phys), intent(in)  :: kbfs              ! surface kinematic buoyancy flux          [ m K s-1 ]
+
+        real(kind_phys)              :: obukhov_length    ! Obukhov length                           [ m       ]
+
+        ! Added sign(...) term to prevent division by 0 and using the fact that
+        ! `kbfs = \overline{(w' \theta_v')}_s`
+        obukhov_length = -thvs * ustar**3 /                             &
+                         (g*karman*(kbfs + sign(1.e-10_kind_phys,kbfs)))
+    end function calc_obukhov_length
+
+    pure elemental function calc_virtual_temperature(temperature, specific_humidity, zvir) result(virtual_temperature)
+        ! Williamson, D., Kiehl, J., Ramanathan, V., Dickinson, R., & Hack, J. (1987).
+        ! Description of the NCAR Community Climate Model (CCM1).
+        ! University Corporation for Atmospheric Research. https://doi.org/10.5065/D6TB14WH (Original work published 1987)
+        ! Equation 2.a.7
+
+        real(kind_phys), intent(in) :: temperature
+        real(kind_phys), intent(in) :: specific_humidity
+        real(kind_phys), intent(in) :: zvir    ! rh2o/rair - 1
+
+        real(kind_phys)             :: virtual_temperature
+
+        virtual_temperature = temperature * (1.0_kind_phys + zvir*specific_humidity)
+    end function calc_virtual_temperature
+
+    pure elemental function calc_eddy_flux_coefficient(mixing_length_squared, &
+                                                       richardson_number,     &
+                                                       shear_squared)         &
+                                                       result(kvf)
+        ! Computes exchange coefficients for turbulent flows.
+        !
+        ! The stable case (Richardson number, Ri>0) is taken from Holtslag and
+        ! Beljaars (1989), ECMWF proceedings.
+        ! The unstable case (Richardson number, Ri<0) is taken from  CCM1.
+
+        real(kind_phys), intent(in)  :: mixing_length_squared     ! [ m2     ]
+        real(kind_phys), intent(in)  :: richardson_number         ! [ 1      ]
+        real(kind_phys), intent(in)  :: shear_squared             ! [ s-2    ]
+
+        real(kind_phys)              :: kvf    ! Eddy diffusivity ! [ m2 s-1 ]
+
+        real(kind_phys)              :: fofri  ! f(ri)
+        real(kind_phys)              :: kvn    ! Neutral Kv
+
+        if( richardson_number < 0.0_kind_phys ) then
+            fofri = unstable_gradient_richardson_stability_parameter(richardson_number)
+        else
+            fofri = stable_gradient_richardson_stability_parameter(richardson_number)
+        end if
+        kvn = neutral_exchange_coefficient(mixing_length_squared, shear_squared)
+        kvf = max( minimum_eddy_flux_coefficient, kvn * fofri )
+    end function calc_eddy_flux_coefficient
+
+    pure elemental function calc_free_atm_eddy_flux_coefficient(mixing_length_squared, &
+                                                                richardson_number,     &
+                                                                shear_squared)         &
+                                                                result(kvf)
+        ! same as austausch_atm but only mixing for Ri<0
+        ! i.e. no background mixing and mixing for Ri>0
+
+        real(kind_phys), intent(in)  :: mixing_length_squared ! [ m2     ]
+        real(kind_phys), intent(in)  :: richardson_number     ! [ 1      ]
+        real(kind_phys), intent(in)  :: shear_squared         ! [ s-2    ]
+
+        real(kind_phys)              :: kvf                   ! [ m2 s-1 ]
+
+        real(kind_phys)              :: fofri ! f(ri)
+        real(kind_phys)              :: kvn   ! Neutral Kv
+
+        kvf = 0.0_kind_phys
+        if( richardson_number < 0.0_kind_phys ) then
+            fofri = unstable_gradient_richardson_stability_parameter(richardson_number)
+            kvn = neutral_exchange_coefficient(mixing_length_squared, shear_squared)
+            kvf = kvn * fofri
+        end if
+    end function calc_free_atm_eddy_flux_coefficient
+
+    pure elemental function unstable_gradient_richardson_stability_parameter(richardson_number) result(modifier)
+        ! Williamson, D., Kiehl, J., Ramanathan, V., Dickinson, R., & Hack, J. (1987).
+        ! Description of the NCAR Community Climate Model (CCM1).
+        ! University Corporation for Atmospheric Research. https://doi.org/10.5065/D6TB14WH (Original work published 1987)
+        ! Equation 2.f.13
+        ! \sqrt{ 1-18*Ri }
+
+        real(kind_phys), intent(in)  :: richardson_number
+
+        real(kind_phys)              :: modifier
+
+        modifier = sqrt( 1._kind_phys - 18._kind_phys * richardson_number )
+    end function unstable_gradient_richardson_stability_parameter
+
+    pure elemental function stable_gradient_richardson_stability_parameter(richardson_number) result(modifier)
+        ! Holtslag, A. A. M., and Beljaars A. C. M. , 1989: Surface flux parameterization schemes: Developments and experiences at KNMI.
+        ! ECMWF Workshop on Parameterization of Fluxes and Land Surface, Reading, United Kingdom, ECMWF, 121–147.
+        ! equation 20, page 140
+        ! Originally used published equation from CCM1, 2.f.12, page 11
+        ! \frac{1}{1+10*Ri*(1+8*Ri)}
+
+        real(kind_phys), intent(in)  :: richardson_number
+
+        real(kind_phys)              :: modifier
+
+        modifier = 1.0_kind_phys /                                                                                              &
+                 ( 1.0_kind_phys + 10.0_kind_phys * richardson_number * ( 1.0_kind_phys + 8.0_kind_phys * richardson_number ) )
+    end function stable_gradient_richardson_stability_parameter
+
+    pure elemental function neutral_exchange_coefficient(mixing_length_squared, shear_squared) result(neutral_k)
+        ! Williamson, D., Kiehl, J., Ramanathan, V., Dickinson, R., & Hack, J. (1987).
+        ! Description of the NCAR Community Climate Model (CCM1).
+        ! University Corporation for Atmospheric Research. https://doi.org/10.5065/D6TB14WH (Original work published 1987)
+        ! Equation 2.f.15, page 12
+        ! NOTE: shear_squared variable currently (01/2025) computed in hb_diff.F90 (s2 in trbintd(...)) matches referenced equation.
+
+        real(kind_phys), intent(in) :: mixing_length_squared ! [ m2     ]
+        real(kind_phys), intent(in) :: shear_squared         ! [ s-2    ]
+
+        real(kind_phys)             :: neutral_k
+
+        neutral_k = mixing_length_squared * sqrt(shear_squared)
+    end function neutral_exchange_coefficient
+end module atmos_phys_pbl_utils
diff --git a/test/unit-test/CMakeLists.txt b/test/unit-test/CMakeLists.txt
@@ -32,6 +32,15 @@ add_library(utilities ${UTILITIES_SRC})
 target_compile_options(utilities PRIVATE -ffree-line-length-none)
 target_include_directories(utilities PUBLIC ${CMAKE_CURRENT_BINARY_DIR})
 
+set(PHYS_UTILS_SRC
+    ../../phys_utils/atmos_phys_pbl_utils.F90
+    include/ccpp_kinds.F90
+)
+
+add_library(phys_utils ${PHYS_UTILS_SRC})
+target_compile_options(phys_utils PRIVATE -ffree-line-length-none)
+target_include_directories(phys_utils PUBLIC ${CMAKE_CURRENT_BINARY_DIR})
+
 if(ATMOSPHERIC_PHYSICS_ENABLE_TESTS OR ATMOSPHERIC_PHYSICS_ENABLE_CODE_COVERAGE)
   enable_testing()
   add_subdirectory(tests)

diff --git a/test/unit-test/tests/CMakeLists.txt b/test/unit-test/tests/CMakeLists.txt
@@ -1,2 +1,2 @@
 add_subdirectory(utilities)
-
+add_subdirectory(phys_utils)
diff --git a/test/unit-test/tests/phys_utils/CMakeLists.txt b/test/unit-test/tests/phys_utils/CMakeLists.txt
@@ -0,0 +1,4 @@
+add_pfunit_ctest(phys_utils_tests
+  TEST_SOURCES test_atmos_pbl_utils.pf
+  LINK_LIBRARIES phys_utils
+)
diff --git a/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf b/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf
@@ -0,0 +1,25 @@
+@test
+subroutine test_free_eddy_coef_is_zero_when_ri_equals_zero()
+    use funit
+    use atmos_phys_pbl_utils, only : calc_free_atm_eddy_flux_coefficient
+    use ccpp_kinds, only: kind_phys
+
+    real(kind_phys) :: kvf
+
+    kvf = calc_free_atm_eddy_flux_coefficient(30.0_kind_phys, 0.0_kind_phys, 0.01_kind_phys)
+
+    @assertEqual(0.0_kind_phys, kvf)
+end subroutine test_free_eddy_coef_is_zero_when_ri_equals_zero
+
+@test
+subroutine test_free_eddy_atm_coef_is_zero_when_ri_stable_near_zero()
+    use funit
+    use atmos_phys_pbl_utils, only : calc_free_atm_eddy_flux_coefficient
+    use ccpp_kinds, only: kind_phys
+
+    real(kind_phys) :: kvf
+
+    kvf = calc_free_atm_eddy_flux_coefficient(30.0_kind_phys, nearest(0.0_kind_phys, 1.0_kind_phys), 0.01_kind_phys)
+
+    @assertEqual(0.0_kind_phys, kvf)
+end subroutine test_free_eddy_atm_coef_is_zero_when_ri_stable_near_zero
Original file line number	Diff line number	Diff line change
		@@ -1,2 +1,2 @@
		add_subdirectory(utilities)

		add_subdirectory(phys_utils)