PBL utils transfer

.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

phys_utils/atmos_pbl_utils.F90

@@ -0,0 +1,213 @@
+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_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_rrho(rair, surface_temperature, pmid) result(rrho)
+        ! air density reciprocal
+        real(kind_phys), intent(in) :: rair  ! gas constant for dry air
+        real(kind_phys), intent(in) :: surface_temperature
+        real(kind_phys), intent(in) :: pmid  ! midpoint pressure (bottom level)
+
+        real(kind_phys)             :: rrho  ! 1./bottom level density
+
+        rrho = rair * surface_temperature / pmid
+    end function calc_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/m2]
+        real(kind_phys), intent(in) :: tauy              ! surface v stress [N/m2]
+        real(kind_phys), intent(in) :: rrho              ! 1./bottom level density
+
+        real(kind_phys)             :: friction_velocity ! surface friction velocity [m/s]
+
+        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/m2)
+        real(kind_phys), intent(in)  :: rrho   ! 1./bottom level density [ m3/kg ]
+        real(kind_phys), intent(in)  :: cpair  ! specific heat of dry air
+
+        real(kind_phys)              :: khfs   ! sfc kinematic heat flux [mK/s]
+
+        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/m2/s)
+        real(kind_phys), intent(in)  :: rrho ! 1./bottom level density [ m3/kg ]
+
+        real(kind_phys)              :: kqfs ! sfc kinematic water vapor flux [m/s]
+
+        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  ! sfc kinematic heat flux [mK/s]
+        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  ! sfc kinematic water vapor flux [m/s]
+
+        real(kind_phys)             :: kbfs  ! sfc kinematic buoyancy flux [mK/s] (`kbfs = \overline{(w' \theta'_v)}_s`)
+
+        kbfs = khfs + zvir*ths*kqfs
+    end function calc_kinematic_buoyancy_flux
+
+    pure elemental function calc_obukhov_length(thvs, ustar, g, vk, 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
+        real(kind_phys), intent(in)  :: ustar             ! Surface friction velocity [ m/s ]
+        real(kind_phys), intent(in)  :: g                 ! acceleration of gravity
+        real(kind_phys), intent(in)  :: vk                ! Von Karman's constant
+        real(kind_phys), intent(in)  :: kbfs              ! sfc kinematic buoyancy flux [m*K/s]
+
+        real(kind_phys)              :: obukhov_length    ! Obukhov length
+
+        ! 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*vk*(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
+        real(kind_phys), intent(in)  :: richardson_number
+        real(kind_phys), intent(in)  :: shear_squared
+
+        real(kind_phys)              :: kvf    ! Eddy diffusivity for heat and tracers
+
+        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
+        real(kind_phys), intent(in)  :: richardson_number
+        real(kind_phys), intent(in)  :: shear_squared
+
+        real(kind_phys)              :: kvf
+
+        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 vriable currently (01/2025) computed in hb_diff.F90 (trbintd) matches referenced equation.
+
+        real(kind_phys), intent(in) :: mixing_length_squared
+        real(kind_phys), intent(in) :: shear_squared
+
+        real(kind_phys)             :: neutral_k
+
+        neutral_k = mixing_length_squared * sqrt(shear_squared)
+    end function neutral_exchange_coefficient
+end module atmos_phys_pbl_utils

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_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)

test/unit-test/tests/CMakeLists.txt

@@ -1,2 +1,2 @@
 add_subdirectory(utilities)
-
+add_subdirectory(phys_utils)

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
+)

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