From adbad1c6aac264278fdff07d9e2b5661f430a378 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 13 Jan 2025 14:11:09 -0700 Subject: [PATCH 01/29] Initial PBL utils refactor. --- phys_utils/pbl_utils.F90 | 244 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 244 insertions(+) create mode 100644 phys_utils/pbl_utils.F90 diff --git a/phys_utils/pbl_utils.F90 b/phys_utils/pbl_utils.F90 new file mode 100644 index 00000000..c26ed14d --- /dev/null +++ b/phys_utils/pbl_utils.F90 @@ -0,0 +1,244 @@ +module pbl_utils + + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: calc_rrho + public :: calc_ustar + 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 :: austausch_atm + public :: austausch_atm_free + + real(r8), parameter :: ustar_min = 0.01_r8 + real(r8), parameter :: zkmin = 0.01_r8 ! Minimum kneutral*f(ri). CCM1 2.f.14 + +contains + + pure elemental function calc_rrho(rair, t, pmid) result(rrho) + ! air density reciprocal + + real(r8), intent(in) :: t ! surface temperature + real(r8), intent(in) :: pmid ! midpoint pressure (bottom level) + real(r8), intent(in) :: rair ! gas constant for dry air + + real(r8) :: rrho ! 1./bottom level density + + rrho = rair * t / pmid + end function calc_rrho + + pure elemental function calc_ustar(taux, tauy, rrho) result(ustar) + ! 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 181 + + real(r8), intent(in) :: taux ! surface u stress [N/m2] + real(r8), intent(in) :: tauy ! surface v stress [N/m2] + real(r8), intent(in) :: rrho ! 1./bottom level density + + real(r8) :: ustar ! surface friction velocity [m/s] + + ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min ) + end function calc_ustar + + pure elemental function calc_kinematic_heat_flux(shflx, rrho, cpair) result(khfs) + 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) :: cpair ! specific heat of dry air + + real(r8) :: 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(r8), intent(in) :: qflx ! water vapor flux (kg/m2/s) + real(r8), intent(in) :: rrho ! 1./bottom level density [ m3/kg ] + + real(r8) :: 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(r8), intent(in) :: khfs ! sfc kinematic heat flux [mK/s] + real(r8), intent(in) :: zvir ! rh2o/rair - 1 + real(r8), intent(in) :: ths ! potential temperature at surface [K] + real(r8), intent(in) :: kqfs ! sfc kinematic water vapor flux [m/s] + + real(r8) :: kbfs ! sfc kinematic buoyancy flux [m^2/s^3] + + 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 + + real(r8), intent(in) :: thvs ! virtual potential temperature at surface + real(r8), intent(in) :: ustar ! Surface friction velocity [ m/s ] + real(r8), intent(in) :: g ! acceleration of gravity + real(r8), intent(in) :: vk ! Von Karman's constant + real(r8), intent(in) :: kbfs ! sfc kinematic buoyancy flux [m*K/s] + + real(r8) :: obukhov_length(n) ! 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_r8,kbfs))) + end function calc_obukhov_length + + pure elemental function calc_virtual_temperature(t, q, zvir) result(virtem) + ! 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(r8), intent(in) :: t, q ! temperature and specific humidity + real(r8), intent(in) :: zvir ! rh2o/rair - 1 + + real(r8) :: virtual_temperature + + virtual_temperature = t * (1.0_r8 + zvir*q) + end function calc_virtual_temperature + + pure function austausch_atm(pcols, ncol, pver, ntop, nbot, ml2, ri, s2) result(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) ! + ! ! + !---------------------------------------------------------------------- ! + + 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 number + + real(r8) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers + + real(r8) :: fofri ! f(ri) + real(r8) :: kvn ! Neutral Kv + integer :: i ! Longitude index + integer :: k ! Vertical index + + kvf(:ncol,:) = 0.0_r8 + + ! Compute the atmosphere free vertical diffusion coefficients: kvh = kvq = kvm. + do k = ntop, nbot - 1 + do i = 1, ncol + if( ri(i,k) < 0.0_r8 ) then + fofri = unstable_gradient_richardson_stability_parameter(ri(i,k)) + else + fofri = stable_gradient_richardson_stability_parameter(ri(i,k)) + end if + kvn = neutral_exchange_coefficient(ml2(k), s2(i,k)) + kvf(i,k+1) = max( zkmin, kvn * fofri ) + end do + end do + + end function austausch_atm + + pure function austausch_atm_free(pcols, ncol, pver, ntop, nbot, ml2, ri, s2) result(kvf) + !---------------------------------------------------------------------- ! + ! ! + ! same as austausch_atm but only mixing for Ri<0 ! + ! i.e. no background mixing and mixing for Ri>0 ! + ! ! + !---------------------------------------------------------------------- ! + + 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 + + real(r8) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers + + real(r8) :: fofri ! f(ri) + real(r8) :: kvn ! Neutral Kv + integer :: i ! Longitude index + integer :: k ! Vertical index + + 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 = unstable_gradient_richardson_stability_parameter(ri(i,k)) + kvn = neutral_exchange_coefficient(ml2(k), s2(i,k)) + kvf(i,k+1) = kvn * fofri + end if + end do + end do + end function austausch_atm_free + + 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 + + real(r8), intent(in) :: richardson_number + + real(r8) :: modifier + + modifier = sqrt( 1._r8 - 18._r8 * richardson_number ) + end function unstable_gradient_richardson_stability_parameter + + pure elemental function stable_gradient_richardson_stability_parameter(richardson_number) result(modifier) + ! ECMWF 74888 Surface flux parameterization schemes: developments and experiences at KNMI (eq. 20) + ! Originally used published equation from in CCM1 2.f.12 + + real(r8), intent(in) :: richardson_number + + real(r8) :: modifier + + modifier = 1.0_r8 / ( 1.0_r8 + 10.0_r8 * richardson_number * ( 1.0_r8 + 8.0_r8 * richardson_number ) ) + end function stable_gradient_richardson_stability_parameter + + pure elemental function neutral_exchange_coefficient(mixing_length, 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 + + real(r8), intent(in) :: mixing_length + real(r8), intent(in) :: shear_squared + + real(r8) :: neutral_k + + neutral_k = mixing_length * sqrt(shear_squared) + end function neutral_exchange_coefficient +end module pbl_utils \ No newline at end of file From 06e03285b3ae9772b746d46226f189d7f100e635 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Tue, 14 Jan 2025 08:12:23 -0700 Subject: [PATCH 02/29] Updating documentation and clarifying variable names. --- phys_utils/pbl_utils.F90 | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/phys_utils/pbl_utils.F90 b/phys_utils/pbl_utils.F90 index c26ed14d..4ed6a6cd 100644 --- a/phys_utils/pbl_utils.F90 +++ b/phys_utils/pbl_utils.F90 @@ -218,8 +218,10 @@ pure elemental function unstable_gradient_richardson_stability_parameter(richard end function unstable_gradient_richardson_stability_parameter pure elemental function stable_gradient_richardson_stability_parameter(richardson_number) result(modifier) - ! ECMWF 74888 Surface flux parameterization schemes: developments and experiences at KNMI (eq. 20) - ! Originally used published equation from in CCM1 2.f.12 + ! 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 real(r8), intent(in) :: richardson_number @@ -228,17 +230,18 @@ pure elemental function stable_gradient_richardson_stability_parameter(richardso modifier = 1.0_r8 / ( 1.0_r8 + 10.0_r8 * richardson_number * ( 1.0_r8 + 8.0_r8 * richardson_number ) ) end function stable_gradient_richardson_stability_parameter - pure elemental function neutral_exchange_coefficient(mixing_length, shear_squared) result(neutral_k) + 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 + ! Equation 2.f.15, page 12 + ! NOTE: shear_squared vriable currently (01/2025) computed in hbdiff matches references equation. - real(r8), intent(in) :: mixing_length + real(r8), intent(in) :: mixing_length_squared real(r8), intent(in) :: shear_squared real(r8) :: neutral_k - neutral_k = mixing_length * sqrt(shear_squared) + neutral_k = mixing_length_squared * sqrt(shear_squared) end function neutral_exchange_coefficient end module pbl_utils \ No newline at end of file From aba5ec1fef094574ba2b8fb06157586e28790fb6 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Tue, 14 Jan 2025 08:16:12 -0700 Subject: [PATCH 03/29] Minor clarification. --- phys_utils/pbl_utils.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phys_utils/pbl_utils.F90 b/phys_utils/pbl_utils.F90 index 4ed6a6cd..6d3cca93 100644 --- a/phys_utils/pbl_utils.F90 +++ b/phys_utils/pbl_utils.F90 @@ -235,7 +235,7 @@ pure elemental function neutral_exchange_coefficient(mixing_length_squared, shea ! 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 hbdiff matches references equation. + ! NOTE: shear_squared vriable currently (01/2025) computed in hb_diff.F90 (trbintd) matches references equation. real(r8), intent(in) :: mixing_length_squared real(r8), intent(in) :: shear_squared From 9fbafbbbe48c6be38ffab7485c2161cd09baad76 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Wed, 15 Jan 2025 07:27:00 -0700 Subject: [PATCH 04/29] Replacing r8 with kind_phys. Additional whitespace and spelling corrections. --- phys_utils/pbl_utils.F90 | 152 +++++++++++++++++++-------------------- 1 file changed, 76 insertions(+), 76 deletions(-) diff --git a/phys_utils/pbl_utils.F90 b/phys_utils/pbl_utils.F90 index 6d3cca93..8cc639b3 100644 --- a/phys_utils/pbl_utils.F90 +++ b/phys_utils/pbl_utils.F90 @@ -1,4 +1,4 @@ -module pbl_utils +module atmos_phys_pbl_utils use ccpp_kinds, only: kind_phys @@ -15,19 +15,19 @@ module pbl_utils public :: austausch_atm public :: austausch_atm_free - real(r8), parameter :: ustar_min = 0.01_r8 - real(r8), parameter :: zkmin = 0.01_r8 ! Minimum kneutral*f(ri). CCM1 2.f.14 + real(kind_phys), parameter :: ustar_min = 0.01_kind_phys + real(kind_phys), parameter :: zkmin = 0.01_kind_phys ! Minimum kneutral*f(ri). CCM1 2.f.14 contains pure elemental function calc_rrho(rair, t, pmid) result(rrho) ! air density reciprocal - real(r8), intent(in) :: t ! surface temperature - real(r8), intent(in) :: pmid ! midpoint pressure (bottom level) - real(r8), intent(in) :: rair ! gas constant for dry air + real(kind_phys), intent(in) :: t ! surface temperature + real(kind_phys), intent(in) :: pmid ! midpoint pressure (bottom level) + real(kind_phys), intent(in) :: rair ! gas constant for dry air - real(r8) :: rrho ! 1./bottom level density + real(kind_phys) :: rrho ! 1./bottom level density rrho = rair * t / pmid end function calc_rrho @@ -41,41 +41,41 @@ pure elemental function calc_ustar(taux, tauy, rrho) result(ustar) ! DOI: https://doi.org/10.1007/978-94-009-3027-8 ! Equation 2.10b, page 181 - real(r8), intent(in) :: taux ! surface u stress [N/m2] - real(r8), intent(in) :: tauy ! surface v stress [N/m2] - real(r8), intent(in) :: rrho ! 1./bottom level density + 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(r8) :: ustar ! surface friction velocity [m/s] + real(kind_phys) :: ustar ! surface friction velocity [m/s] ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min ) end function calc_ustar pure elemental function calc_kinematic_heat_flux(shflx, rrho, cpair) result(khfs) - 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) :: cpair ! specific heat of dry air + 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(r8) :: khfs ! sfc kinematic heat flux [mK/s] + 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(r8), intent(in) :: qflx ! water vapor flux (kg/m2/s) - real(r8), intent(in) :: rrho ! 1./bottom level density [ m3/kg ] + 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(r8) :: kqfs ! sfc kinematic water vapor flux [m/s] + 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(r8), intent(in) :: khfs ! sfc kinematic heat flux [mK/s] - real(r8), intent(in) :: zvir ! rh2o/rair - 1 - real(r8), intent(in) :: ths ! potential temperature at surface [K] - real(r8), intent(in) :: kqfs ! sfc kinematic water vapor flux [m/s] + 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(r8) :: kbfs ! sfc kinematic buoyancy flux [m^2/s^3] + real(kind_phys) :: kbfs ! sfc kinematic buoyancy flux [m^2/s^3] kbfs = khfs + zvir*ths*kqfs end function calc_kinematic_buoyancy_flux @@ -85,31 +85,31 @@ pure elemental function calc_obukhov_length(thvs, ustar, g, vk, kbfs) result(obu ! DOI: https://doi.org/10.1007/978-94-009-3027-8 ! Equation 5.7c, page 181 - real(r8), intent(in) :: thvs ! virtual potential temperature at surface - real(r8), intent(in) :: ustar ! Surface friction velocity [ m/s ] - real(r8), intent(in) :: g ! acceleration of gravity - real(r8), intent(in) :: vk ! Von Karman's constant - real(r8), intent(in) :: kbfs ! sfc kinematic buoyancy flux [m*K/s] + 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(r8) :: obukhov_length(n) ! Obukhov length + 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_r8,kbfs))) + 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(t, q, zvir) result(virtem) + pure elemental function calc_virtual_temperature(t, q, 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(r8), intent(in) :: t, q ! temperature and specific humidity - real(r8), intent(in) :: zvir ! rh2o/rair - 1 + real(kind_phys), intent(in) :: t, q ! temperature and specific humidity + real(kind_phys), intent(in) :: zvir ! rh2o/rair - 1 - real(r8) :: virtual_temperature + real(kind_phys) :: virtual_temperature - virtual_temperature = t * (1.0_r8 + zvir*q) + virtual_temperature = t * (1.0_kind_phys + zvir*q) end function calc_virtual_temperature pure function austausch_atm(pcols, ncol, pver, ntop, nbot, ml2, ri, s2) result(kvf) @@ -133,28 +133,28 @@ pure function austausch_atm(pcols, ncol, pver, ntop, nbot, ml2, ri, s2) result(k ! ! !---------------------------------------------------------------------- ! - 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 number + 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(kind_phys), intent(in) :: ml2(pver+1) ! Mixing lengths squared + real(kind_phys), intent(in) :: s2(pcols,pver) ! Shear squared + real(kind_phys), intent(in) :: ri(pcols,pver) ! Richardson number - real(r8) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers + real(kind_phys) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers - real(r8) :: fofri ! f(ri) - real(r8) :: kvn ! Neutral Kv - integer :: i ! Longitude index - integer :: k ! Vertical index + real(kind_phys) :: fofri ! f(ri) + real(kind_phys) :: kvn ! Neutral Kv + integer :: i ! Longitude index + integer :: k ! Vertical index - kvf(:ncol,:) = 0.0_r8 + kvf(:ncol,:) = 0.0_kind_phys ! Compute the atmosphere free vertical diffusion coefficients: kvh = kvq = kvm. do k = ntop, nbot - 1 do i = 1, ncol - if( ri(i,k) < 0.0_r8 ) then + if( ri(i,k) < 0.0_kind_phys ) then fofri = unstable_gradient_richardson_stability_parameter(ri(i,k)) else fofri = stable_gradient_richardson_stability_parameter(ri(i,k)) @@ -174,28 +174,28 @@ pure function austausch_atm_free(pcols, ncol, pver, ntop, nbot, ml2, ri, s2) res ! ! !---------------------------------------------------------------------- ! - 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 + 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(kind_phys), intent(in) :: ml2(pver+1) ! Mixing lengths squared + real(kind_phys), intent(in) :: s2(pcols,pver) ! Shear squared + real(kind_phys), intent(in) :: ri(pcols,pver) ! Richardson no - real(r8) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers + real(kind_phys) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers - real(r8) :: fofri ! f(ri) - real(r8) :: kvn ! Neutral Kv - integer :: i ! Longitude index - integer :: k ! Vertical index + real(kind_phys) :: fofri ! f(ri) + real(kind_phys) :: kvn ! Neutral Kv + integer :: i ! Longitude index + integer :: k ! Vertical index - kvf(:ncol,:) = 0.0_r8 + kvf(:ncol,:) = 0.0_kind_phys ! 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 + if( ri(i,k) < 0.0_kind_phys ) then fofri = unstable_gradient_richardson_stability_parameter(ri(i,k)) kvn = neutral_exchange_coefficient(ml2(k), s2(i,k)) kvf(i,k+1) = kvn * fofri @@ -210,11 +210,11 @@ pure elemental function unstable_gradient_richardson_stability_parameter(richard ! University Corporation for Atmospheric Research. https://doi.org/10.5065/D6TB14WH (Original work published 1987) ! Equation 2.f.13 - real(r8), intent(in) :: richardson_number + real(kind_phys), intent(in) :: richardson_number - real(r8) :: modifier + real(kind_phys) :: modifier - modifier = sqrt( 1._r8 - 18._r8 * richardson_number ) + 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) @@ -223,11 +223,11 @@ pure elemental function stable_gradient_richardson_stability_parameter(richardso ! equation 20, page 140 ! Originally used published equation from CCM1, 2.f.12, page 11 - real(r8), intent(in) :: richardson_number + real(kind_phys), intent(in) :: richardson_number - real(r8) :: modifier + real(kind_phys) :: modifier - modifier = 1.0_r8 / ( 1.0_r8 + 10.0_r8 * richardson_number * ( 1.0_r8 + 8.0_r8 * richardson_number ) ) + 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) @@ -235,13 +235,13 @@ pure elemental function neutral_exchange_coefficient(mixing_length_squared, shea ! 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 references equation. + ! NOTE: shear_squared vriable currently (01/2025) computed in hb_diff.F90 (trbintd) matches referenced equation. - real(r8), intent(in) :: mixing_length_squared - real(r8), intent(in) :: shear_squared + real(kind_phys), intent(in) :: mixing_length_squared + real(kind_phys), intent(in) :: shear_squared - real(r8) :: neutral_k + real(kind_phys) :: neutral_k neutral_k = mixing_length_squared * sqrt(shear_squared) end function neutral_exchange_coefficient -end module pbl_utils \ No newline at end of file +end module atmos_phys_pbl_utils \ No newline at end of file From 1fcc007ada10267d8b483a50ba59bbea463eeff0 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Wed, 15 Jan 2025 11:22:15 -0700 Subject: [PATCH 05/29] Update function name and fix equation page reference. --- phys_utils/pbl_utils.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/phys_utils/pbl_utils.F90 b/phys_utils/pbl_utils.F90 index 8cc639b3..c1c9aa17 100644 --- a/phys_utils/pbl_utils.F90 +++ b/phys_utils/pbl_utils.F90 @@ -6,7 +6,7 @@ module atmos_phys_pbl_utils private public :: calc_rrho - public :: calc_ustar + public :: calc_friction_velocity public :: calc_kinematic_heat_flux public :: calc_kinematic_water_vapor_flux public :: calc_kinematic_buoyancy_flux @@ -32,14 +32,14 @@ pure elemental function calc_rrho(rair, t, pmid) result(rrho) rrho = rair * t / pmid end function calc_rrho - pure elemental function calc_ustar(taux, tauy, rrho) result(ustar) + pure elemental function calc_friction_velocity(taux, tauy, rrho) result(ustar) ! 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 181 + ! 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] @@ -48,7 +48,7 @@ pure elemental function calc_ustar(taux, tauy, rrho) result(ustar) real(kind_phys) :: ustar ! surface friction velocity [m/s] ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min ) - end function calc_ustar + 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) From 8e106fbf10e5075fab8b557f67ad34f1f333d004 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Fri, 17 Jan 2025 12:51:08 -0700 Subject: [PATCH 06/29] Elemental refactoring of austausch functions. Refactored variable names. --- phys_utils/pbl_utils.F90 | 125 ++++++++++++++++----------------------- 1 file changed, 51 insertions(+), 74 deletions(-) diff --git a/phys_utils/pbl_utils.F90 b/phys_utils/pbl_utils.F90 index c1c9aa17..daf5f7b6 100644 --- a/phys_utils/pbl_utils.F90 +++ b/phys_utils/pbl_utils.F90 @@ -20,19 +20,19 @@ module atmos_phys_pbl_utils contains - pure elemental function calc_rrho(rair, t, pmid) result(rrho) + pure elemental function calc_rrho(rair, surface_temperature, pmid) result(rrho) ! air density reciprocal - real(kind_phys), intent(in) :: t ! surface temperature + real(kind_phys), intent(in) :: surface_temperature real(kind_phys), intent(in) :: pmid ! midpoint pressure (bottom level) real(kind_phys), intent(in) :: rair ! gas constant for dry air real(kind_phys) :: rrho ! 1./bottom level density - rrho = rair * t / pmid + rrho = rair * surface_temperature / pmid end function calc_rrho - pure elemental function calc_friction_velocity(taux, tauy, rrho) result(ustar) + 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 ! @@ -41,13 +41,13 @@ pure elemental function calc_friction_velocity(taux, tauy, rrho) result(ustar) ! 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), 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) :: ustar ! surface friction velocity [m/s] + real(kind_phys) :: friction_velocity ! surface friction velocity [m/s] - ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min ) + friction_velocity = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min ) end function calc_friction_velocity pure elemental function calc_kinematic_heat_flux(shflx, rrho, cpair) result(khfs) @@ -98,21 +98,25 @@ pure elemental function calc_obukhov_length(thvs, ustar, g, vk, kbfs) result(obu 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(t, q, zvir) result(virtual_temperature) + 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) :: t, q ! temperature and specific humidity + 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 = t * (1.0_kind_phys + zvir*q) + virtual_temperature = temperature * (1.0_kind_phys + zvir*specific_humidity) end function calc_virtual_temperature - pure function austausch_atm(pcols, ncol, pver, ntop, nbot, ml2, ri, s2) result(kvf) + pure elemental function austausch_atm(mixing_length_squared, & + richardson_number, & + shear_squared) & + result(kvf) !---------------------------------------------------------------------- ! ! ! ! Purpose: Computes exchange coefficients for free turbulent flows. ! @@ -133,75 +137,48 @@ pure function austausch_atm(pcols, ncol, pver, ntop, nbot, ml2, ri, s2) result(k ! ! !---------------------------------------------------------------------- ! - 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(kind_phys), intent(in) :: ml2(pver+1) ! Mixing lengths squared - real(kind_phys), intent(in) :: s2(pcols,pver) ! Shear squared - real(kind_phys), intent(in) :: ri(pcols,pver) ! Richardson number - - real(kind_phys) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers - - real(kind_phys) :: fofri ! f(ri) - real(kind_phys) :: kvn ! Neutral Kv - integer :: i ! Longitude index - integer :: k ! Vertical index - - kvf(:ncol,:) = 0.0_kind_phys - - ! Compute the atmosphere free vertical diffusion coefficients: kvh = kvq = kvm. - do k = ntop, nbot - 1 - do i = 1, ncol - if( ri(i,k) < 0.0_kind_phys ) then - fofri = unstable_gradient_richardson_stability_parameter(ri(i,k)) - else - fofri = stable_gradient_richardson_stability_parameter(ri(i,k)) - end if - kvn = neutral_exchange_coefficient(ml2(k), s2(i,k)) - kvf(i,k+1) = max( zkmin, kvn * fofri ) - end do - end do - + 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( zkmin, kvn * fofri ) end function austausch_atm - - pure function austausch_atm_free(pcols, ncol, pver, ntop, nbot, ml2, ri, s2) result(kvf) + + pure elemental function austausch_atm_free(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 - 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(kind_phys), intent(in) :: ml2(pver+1) ! Mixing lengths squared - real(kind_phys), intent(in) :: s2(pcols,pver) ! Shear squared - real(kind_phys), intent(in) :: ri(pcols,pver) ! Richardson no - - real(kind_phys) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers - - real(kind_phys) :: fofri ! f(ri) - real(kind_phys) :: kvn ! Neutral Kv - integer :: i ! Longitude index - integer :: k ! Vertical index - - kvf(:ncol,:) = 0.0_kind_phys - - ! 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_kind_phys ) then - fofri = unstable_gradient_richardson_stability_parameter(ri(i,k)) - kvn = neutral_exchange_coefficient(ml2(k), s2(i,k)) - kvf(i,k+1) = kvn * fofri - end if - end do - end do + 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 austausch_atm_free pure elemental function unstable_gradient_richardson_stability_parameter(richardson_number) result(modifier) From 71f5e041ee84a8a981a764064961b660789442ef Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Fri, 17 Jan 2025 12:51:53 -0700 Subject: [PATCH 07/29] Renaming file to not conflict with CAM. --- phys_utils/{pbl_utils.F90 => atmos_pbl_utils.F90} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename phys_utils/{pbl_utils.F90 => atmos_pbl_utils.F90} (100%) diff --git a/phys_utils/pbl_utils.F90 b/phys_utils/atmos_pbl_utils.F90 similarity index 100% rename from phys_utils/pbl_utils.F90 rename to phys_utils/atmos_pbl_utils.F90 From 483f46bf82295e4bae8e433d3136e54419093a0c Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Fri, 17 Jan 2025 12:59:43 -0700 Subject: [PATCH 08/29] Whitespace fixes. --- phys_utils/pbl_utils.F90 | 224 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 224 insertions(+) create mode 100644 phys_utils/pbl_utils.F90 diff --git a/phys_utils/pbl_utils.F90 b/phys_utils/pbl_utils.F90 new file mode 100644 index 00000000..dc00aeac --- /dev/null +++ b/phys_utils/pbl_utils.F90 @@ -0,0 +1,224 @@ +module atmos_phys_pbl_utils + + 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 :: austausch_atm + public :: austausch_atm_free + + real(kind_phys), parameter :: ustar_min = 0.01_kind_phys + real(kind_phys), parameter :: zkmin = 0.01_kind_phys ! Minimum kneutral*f(ri). CCM1 2.f.14 + +contains + + pure elemental function calc_rrho(rair, surface_temperature, pmid) result(rrho) + ! air density reciprocal + + real(kind_phys), intent(in) :: surface_temperature + real(kind_phys), intent(in) :: pmid ! midpoint pressure (bottom level) + real(kind_phys), intent(in) :: rair ! gas constant for dry air + + 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 ), ustar_min ) + 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 [m^2/s^3] + + 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 + + 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 austausch_atm(mixing_length_squared, & + richardson_number, & + shear_squared) & + result(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) ! + ! ! + !---------------------------------------------------------------------- ! + + 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( zkmin, kvn * fofri ) + end function austausch_atm + + pure elemental function austausch_atm_free(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 + + 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 austausch_atm_free + + 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 + + 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 + + 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 \ No newline at end of file From d8c00aa75ad6bb01b2f5b719ba4b8c414e475f89 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Wed, 22 Jan 2025 14:07:57 -0700 Subject: [PATCH 09/29] Removing file added back accidentally --- phys_utils/pbl_utils.F90 | 224 --------------------------------------- 1 file changed, 224 deletions(-) delete mode 100644 phys_utils/pbl_utils.F90 diff --git a/phys_utils/pbl_utils.F90 b/phys_utils/pbl_utils.F90 deleted file mode 100644 index dc00aeac..00000000 --- a/phys_utils/pbl_utils.F90 +++ /dev/null @@ -1,224 +0,0 @@ -module atmos_phys_pbl_utils - - 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 :: austausch_atm - public :: austausch_atm_free - - real(kind_phys), parameter :: ustar_min = 0.01_kind_phys - real(kind_phys), parameter :: zkmin = 0.01_kind_phys ! Minimum kneutral*f(ri). CCM1 2.f.14 - -contains - - pure elemental function calc_rrho(rair, surface_temperature, pmid) result(rrho) - ! air density reciprocal - - real(kind_phys), intent(in) :: surface_temperature - real(kind_phys), intent(in) :: pmid ! midpoint pressure (bottom level) - real(kind_phys), intent(in) :: rair ! gas constant for dry air - - 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 ), ustar_min ) - 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 [m^2/s^3] - - 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 - - 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 austausch_atm(mixing_length_squared, & - richardson_number, & - shear_squared) & - result(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) ! - ! ! - !---------------------------------------------------------------------- ! - - 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( zkmin, kvn * fofri ) - end function austausch_atm - - pure elemental function austausch_atm_free(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 - - 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 austausch_atm_free - - 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 - - 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 - - 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 \ No newline at end of file From 2e6a946f1674e98b22c574c5c1dbde63fa862b84 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Wed, 22 Jan 2025 14:08:44 -0700 Subject: [PATCH 10/29] Whitespace fix. --- phys_utils/atmos_pbl_utils.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/phys_utils/atmos_pbl_utils.F90 b/phys_utils/atmos_pbl_utils.F90 index daf5f7b6..b3976684 100644 --- a/phys_utils/atmos_pbl_utils.F90 +++ b/phys_utils/atmos_pbl_utils.F90 @@ -114,9 +114,9 @@ pure elemental function calc_virtual_temperature(temperature, specific_humidity, end function calc_virtual_temperature pure elemental function austausch_atm(mixing_length_squared, & - richardson_number, & - shear_squared) & - result(kvf) + richardson_number, & + shear_squared) & + result(kvf) !---------------------------------------------------------------------- ! ! ! ! Purpose: Computes exchange coefficients for free turbulent flows. ! @@ -156,9 +156,9 @@ pure elemental function austausch_atm(mixing_length_squared, & end function austausch_atm pure elemental function austausch_atm_free(mixing_length_squared, & - richardson_number, & - shear_squared) & - result(kvf) + richardson_number, & + shear_squared) & + result(kvf) !---------------------------------------------------------------------- ! ! ! ! same as austausch_atm but only mixing for Ri<0 ! @@ -221,4 +221,4 @@ pure elemental function neutral_exchange_coefficient(mixing_length_squared, shea neutral_k = mixing_length_squared * sqrt(shear_squared) end function neutral_exchange_coefficient -end module atmos_phys_pbl_utils \ No newline at end of file +end module atmos_phys_pbl_utils From 9193552576698433d612536985fc534594902f4c Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Wed, 22 Jan 2025 14:09:46 -0700 Subject: [PATCH 11/29] Fixing runtime error. --- phys_utils/atmos_pbl_utils.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/phys_utils/atmos_pbl_utils.F90 b/phys_utils/atmos_pbl_utils.F90 index b3976684..198f1c51 100644 --- a/phys_utils/atmos_pbl_utils.F90 +++ b/phys_utils/atmos_pbl_utils.F90 @@ -173,7 +173,8 @@ pure elemental function austausch_atm_free(mixing_length_squared, & 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) From 860abf9cc7019df8cdbf9a6e11d16c942a109f2b Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Thu, 23 Jan 2025 15:06:51 -0700 Subject: [PATCH 12/29] Reorder parameter declarations to match function parameter list and add note about variable equivalency. --- phys_utils/atmos_pbl_utils.F90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/phys_utils/atmos_pbl_utils.F90 b/phys_utils/atmos_pbl_utils.F90 index 198f1c51..8399fb8e 100644 --- a/phys_utils/atmos_pbl_utils.F90 +++ b/phys_utils/atmos_pbl_utils.F90 @@ -22,10 +22,9 @@ module atmos_phys_pbl_utils 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), intent(in) :: rair ! gas constant for dry air real(kind_phys) :: rrho ! 1./bottom level density @@ -75,7 +74,7 @@ pure elemental function calc_kinematic_buoyancy_flux(khfs, zvir, ths, kqfs) resu 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 [m^2/s^3] + 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 From 8dc628411fc437039391a2cde935b37c1848c0fb Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Thu, 23 Jan 2025 15:13:56 -0700 Subject: [PATCH 13/29] Fixing line spacing. --- phys_utils/atmos_pbl_utils.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/phys_utils/atmos_pbl_utils.F90 b/phys_utils/atmos_pbl_utils.F90 index 8399fb8e..84ed0a33 100644 --- a/phys_utils/atmos_pbl_utils.F90 +++ b/phys_utils/atmos_pbl_utils.F90 @@ -204,7 +204,8 @@ pure elemental function stable_gradient_richardson_stability_parameter(richardso 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 ) ) + 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) From e71de9a01b901f3b30ae83e53568996d7b79ae29 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Thu, 23 Jan 2025 15:17:59 -0700 Subject: [PATCH 14/29] Whitespace fixing --- phys_utils/atmos_pbl_utils.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/phys_utils/atmos_pbl_utils.F90 b/phys_utils/atmos_pbl_utils.F90 index 84ed0a33..1baa7e67 100644 --- a/phys_utils/atmos_pbl_utils.F90 +++ b/phys_utils/atmos_pbl_utils.F90 @@ -94,7 +94,8 @@ pure elemental function calc_obukhov_length(thvs, ustar, g, vk, kbfs) result(obu ! 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))) + 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) @@ -204,7 +205,7 @@ pure elemental function stable_gradient_richardson_stability_parameter(richardso real(kind_phys) :: modifier - modifier = 1.0_kind_phys / & + 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 From bcc145ad9cc06185e58054ec01b008c2d9b9fa67 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Thu, 23 Jan 2025 15:20:00 -0700 Subject: [PATCH 15/29] Fix minimum parameter name. --- phys_utils/atmos_pbl_utils.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/phys_utils/atmos_pbl_utils.F90 b/phys_utils/atmos_pbl_utils.F90 index 1baa7e67..c0804063 100644 --- a/phys_utils/atmos_pbl_utils.F90 +++ b/phys_utils/atmos_pbl_utils.F90 @@ -15,8 +15,8 @@ module atmos_phys_pbl_utils public :: austausch_atm public :: austausch_atm_free - real(kind_phys), parameter :: ustar_min = 0.01_kind_phys - real(kind_phys), parameter :: zkmin = 0.01_kind_phys ! Minimum kneutral*f(ri). CCM1 2.f.14 + real(kind_phys), parameter :: minimum_friction_velocity = 0.01_kind_phys + real(kind_phys), parameter :: zkmin = 0.01_kind_phys ! Minimum kneutral*f(ri). CCM1 2.f.14 contains @@ -46,7 +46,7 @@ pure elemental function calc_friction_velocity(taux, tauy, rrho) result(friction real(kind_phys) :: friction_velocity ! surface friction velocity [m/s] - friction_velocity = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min ) + 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) From 1bf16333a752f5930a4e937bf473778356c88bf1 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Fri, 24 Jan 2025 08:21:21 -0700 Subject: [PATCH 16/29] Adding example unit test. --- test/unit-test/CMakeLists.txt | 9 ++++++ test/unit-test/tests/CMakeLists.txt | 2 +- .../unit-test/tests/phys_utils/CMakeLists.txt | 5 ++++ .../tests/phys_utils/test_atmos_pbl_utils.pf | 29 +++++++++++++++++++ 4 files changed, 44 insertions(+), 1 deletion(-) create mode 100644 test/unit-test/tests/phys_utils/CMakeLists.txt create mode 100644 test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf diff --git a/test/unit-test/CMakeLists.txt b/test/unit-test/CMakeLists.txt index d307934e..590d80e1 100644 --- 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_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 index b4c4714c..705189f1 100644 --- 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 new file mode 100644 index 00000000..1ed0887c --- /dev/null +++ b/test/unit-test/tests/phys_utils/CMakeLists.txt @@ -0,0 +1,5 @@ +add_pfunit_ctest(phys_utils_tests + TEST_SOURCES test_atmos_pbl_utils.pf + LINK_LIBRARIES phys_utils +) +# target_compile_options(phys_utils_tests PRIVATE -std=f2008 -fmax-identifier-length=64) \ No newline at end of file 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 new file mode 100644 index 00000000..bc5bb8e4 --- /dev/null +++ b/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf @@ -0,0 +1,29 @@ +#define GFC_MAX_MANGLED_SYMBOL_LEN 200 + +@test +subroutine test_austausch_free_is_zero_when_richardson_stable_at_zero() + use funit + use atmos_phys_pbl_utils, only : austausch_atm_free + use ccpp_kinds, only: kind_phys + + real(kind_phys) :: kvf + + kvf = austausch_atm_free(30.0_kind_phys, 0.0_kind_phys, 0.01_kind_phys) + + @assertEqual(0.0_kind_phys, kvf) +end subroutine test_austausch_free_is_zero_when_richardson_stable_at_zero + +#define GFC_MAX_MANGLED_SYMBOL_LEN 200 + +@test +subroutine test_austausch_free_is_zero_when_richardson_stable_near_zero() + use funit + use atmos_phys_pbl_utils, only : austausch_atm_free + use ccpp_kinds, only: kind_phys + + real(kind_phys) :: kvf + + kvf = austausch_atm_free(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_austausch_free_is_zero_when_richardson_stable_near_zero \ No newline at end of file From 6ad515cc70d873e1611126bb3020e9be0b491738 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Fri, 24 Jan 2025 08:52:34 -0700 Subject: [PATCH 17/29] Removing unused ifdef and fix test name. --- test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) 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 index bc5bb8e4..cb93b1d2 100644 --- a/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf +++ b/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf @@ -1,7 +1,5 @@ -#define GFC_MAX_MANGLED_SYMBOL_LEN 200 - @test -subroutine test_austausch_free_is_zero_when_richardson_stable_at_zero() +subroutine test_austausch_free_is_zero_when_richardson_equals_zero() use funit use atmos_phys_pbl_utils, only : austausch_atm_free use ccpp_kinds, only: kind_phys @@ -11,9 +9,7 @@ subroutine test_austausch_free_is_zero_when_richardson_stable_at_zero() kvf = austausch_atm_free(30.0_kind_phys, 0.0_kind_phys, 0.01_kind_phys) @assertEqual(0.0_kind_phys, kvf) -end subroutine test_austausch_free_is_zero_when_richardson_stable_at_zero - -#define GFC_MAX_MANGLED_SYMBOL_LEN 200 +end subroutine test_austausch_free_is_zero_when_richardson_equals_zero @test subroutine test_austausch_free_is_zero_when_richardson_stable_near_zero() From 31632c6306bead0c925cdae7afa27c0113c40af0 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Fri, 24 Jan 2025 08:57:01 -0700 Subject: [PATCH 18/29] Removing unused compiler flags. --- test/unit-test/tests/phys_utils/CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/test/unit-test/tests/phys_utils/CMakeLists.txt b/test/unit-test/tests/phys_utils/CMakeLists.txt index 1ed0887c..5e3020cc 100644 --- a/test/unit-test/tests/phys_utils/CMakeLists.txt +++ b/test/unit-test/tests/phys_utils/CMakeLists.txt @@ -2,4 +2,3 @@ add_pfunit_ctest(phys_utils_tests TEST_SOURCES test_atmos_pbl_utils.pf LINK_LIBRARIES phys_utils ) -# target_compile_options(phys_utils_tests PRIVATE -std=f2008 -fmax-identifier-length=64) \ No newline at end of file From 933affd0d1484ebf56323eeb96a71c9d7399cacc Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Fri, 24 Jan 2025 08:57:46 -0700 Subject: [PATCH 19/29] Adding missing newline. --- test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 index cb93b1d2..8d15a733 100644 --- a/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf +++ b/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf @@ -22,4 +22,4 @@ subroutine test_austausch_free_is_zero_when_richardson_stable_near_zero() kvf = austausch_atm_free(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_austausch_free_is_zero_when_richardson_stable_near_zero \ No newline at end of file +end subroutine test_austausch_free_is_zero_when_richardson_stable_near_zero From 937a0873456611a452f8bfbe25b27c189689aca6 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Fri, 24 Jan 2025 09:08:54 -0700 Subject: [PATCH 20/29] Add phys_utils to code coverage. --- .github/workflows/unit-tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/unit-tests.yaml b/.github/workflows/unit-tests.yaml index d9003cee..45ef3a57 100644 --- 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 From 8ad300aa631440dec56525ab581716923a8b541b Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 27 Jan 2025 08:42:29 -0700 Subject: [PATCH 21/29] Updating variable names to attempt modernization. --- phys_utils/atmos_pbl_utils.F90 | 66 +++++++++++++--------------------- 1 file changed, 25 insertions(+), 41 deletions(-) diff --git a/phys_utils/atmos_pbl_utils.F90 b/phys_utils/atmos_pbl_utils.F90 index c0804063..9bfafdfe 100644 --- a/phys_utils/atmos_pbl_utils.F90 +++ b/phys_utils/atmos_pbl_utils.F90 @@ -1,4 +1,5 @@ module atmos_phys_pbl_utils + ! Planetary boundary layer related functions used for vertical diffusion schemes. use ccpp_kinds, only: kind_phys @@ -12,11 +13,11 @@ module atmos_phys_pbl_utils public :: calc_kinematic_buoyancy_flux public :: calc_obukhov_length public :: calc_virtual_temperature - public :: austausch_atm - public :: austausch_atm_free + public :: calc_eddy_flux_coefficient + public :: calc_free_eddy_flux_coefficient - real(kind_phys), parameter :: minimum_friction_velocity = 0.01_kind_phys - real(kind_phys), parameter :: zkmin = 0.01_kind_phys ! Minimum kneutral*f(ri). CCM1 2.f.14 + real(kind_phys), parameter :: minimum_friction_velocity = 0.01_kind_phys + real(kind_phys), parameter :: eddy_flux_coefficient_minimum_value = 0.01_kind_phys ! CCM1 2.f.14 contains @@ -113,29 +114,15 @@ pure elemental function calc_virtual_temperature(temperature, specific_humidity, virtual_temperature = temperature * (1.0_kind_phys + zvir*specific_humidity) end function calc_virtual_temperature - pure elemental function austausch_atm(mixing_length_squared, & - richardson_number, & - shear_squared) & - result(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) ! - ! ! - !---------------------------------------------------------------------- ! + 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 @@ -152,19 +139,16 @@ pure elemental function austausch_atm(mixing_length_squared, & fofri = stable_gradient_richardson_stability_parameter(richardson_number) end if kvn = neutral_exchange_coefficient(mixing_length_squared, shear_squared) - kvf = max( zkmin, kvn * fofri ) - end function austausch_atm - - pure elemental function austausch_atm_free(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 ! - ! ! - !---------------------------------------------------------------------- ! + kvf = max( eddy_flux_coefficient_minimum_value, kvn * fofri ) + end function calc_eddy_flux_coefficient + + pure elemental function calc_free_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 @@ -180,7 +164,7 @@ pure elemental function austausch_atm_free(mixing_length_squared, & kvn = neutral_exchange_coefficient(mixing_length_squared, shear_squared) kvf = kvn * fofri end if - end function austausch_atm_free + end function calc_free_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). From 0d9e0be5e3de06a2f02c0c02ae3a83e1427db3d0 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 27 Jan 2025 08:43:19 -0700 Subject: [PATCH 22/29] Updating unit tests missed in previous commit. --- test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) 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 index 8d15a733..3e408bab 100644 --- a/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf +++ b/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf @@ -1,12 +1,12 @@ @test subroutine test_austausch_free_is_zero_when_richardson_equals_zero() use funit - use atmos_phys_pbl_utils, only : austausch_atm_free + use atmos_phys_pbl_utils, only : calc_free_eddy_flux_coefficient use ccpp_kinds, only: kind_phys real(kind_phys) :: kvf - kvf = austausch_atm_free(30.0_kind_phys, 0.0_kind_phys, 0.01_kind_phys) + kvf = calc_free_eddy_flux_coefficient(30.0_kind_phys, 0.0_kind_phys, 0.01_kind_phys) @assertEqual(0.0_kind_phys, kvf) end subroutine test_austausch_free_is_zero_when_richardson_equals_zero @@ -14,12 +14,12 @@ end subroutine test_austausch_free_is_zero_when_richardson_equals_zero @test subroutine test_austausch_free_is_zero_when_richardson_stable_near_zero() use funit - use atmos_phys_pbl_utils, only : austausch_atm_free + use atmos_phys_pbl_utils, only : calc_free_eddy_flux_coefficient use ccpp_kinds, only: kind_phys real(kind_phys) :: kvf - kvf = austausch_atm_free(30.0_kind_phys, nearest(0.0_kind_phys, 1.0_kind_phys), 0.01_kind_phys) + kvf = calc_free_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_austausch_free_is_zero_when_richardson_stable_near_zero From 7538c8687b101ba90f38f3fd7f18c9791e6e0dac Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 27 Jan 2025 08:45:47 -0700 Subject: [PATCH 23/29] Update test names to reflect function name changes. --- test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) 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 index 3e408bab..b6530b4b 100644 --- a/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf +++ b/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf @@ -1,5 +1,5 @@ @test -subroutine test_austausch_free_is_zero_when_richardson_equals_zero() +subroutine test_free_eddy_coef_is_zero_when_richardson_equals_zero() use funit use atmos_phys_pbl_utils, only : calc_free_eddy_flux_coefficient use ccpp_kinds, only: kind_phys @@ -9,10 +9,10 @@ subroutine test_austausch_free_is_zero_when_richardson_equals_zero() kvf = calc_free_eddy_flux_coefficient(30.0_kind_phys, 0.0_kind_phys, 0.01_kind_phys) @assertEqual(0.0_kind_phys, kvf) -end subroutine test_austausch_free_is_zero_when_richardson_equals_zero +end subroutine test_free_eddy_coef_is_zero_when_richardson_equals_zero @test -subroutine test_austausch_free_is_zero_when_richardson_stable_near_zero() +subroutine test_free_eddy_coef_is_zero_when_richardson_stable_near_zero() use funit use atmos_phys_pbl_utils, only : calc_free_eddy_flux_coefficient use ccpp_kinds, only: kind_phys @@ -22,4 +22,4 @@ subroutine test_austausch_free_is_zero_when_richardson_stable_near_zero() kvf = calc_free_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_austausch_free_is_zero_when_richardson_stable_near_zero +end subroutine test_free_eddy_coef_is_zero_when_richardson_stable_near_zero From a3483929e8862f7c02e81890a4cfc0c2921bc2bf Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 27 Jan 2025 08:48:56 -0700 Subject: [PATCH 24/29] Updating minimum value member name. --- phys_utils/atmos_pbl_utils.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/phys_utils/atmos_pbl_utils.F90 b/phys_utils/atmos_pbl_utils.F90 index 9bfafdfe..fd6867e6 100644 --- a/phys_utils/atmos_pbl_utils.F90 +++ b/phys_utils/atmos_pbl_utils.F90 @@ -16,8 +16,8 @@ module atmos_phys_pbl_utils public :: calc_eddy_flux_coefficient public :: calc_free_eddy_flux_coefficient - real(kind_phys), parameter :: minimum_friction_velocity = 0.01_kind_phys - real(kind_phys), parameter :: eddy_flux_coefficient_minimum_value = 0.01_kind_phys ! CCM1 2.f.14 + real(kind_phys), parameter :: minimum_friction_velocity = 0.01_kind_phys + real(kind_phys), parameter :: minimum_eddy_flux_coefficient = 0.01_kind_phys ! CCM1 2.f.14 contains @@ -139,7 +139,7 @@ pure elemental function calc_eddy_flux_coefficient(mixing_length_squared, & fofri = stable_gradient_richardson_stability_parameter(richardson_number) end if kvn = neutral_exchange_coefficient(mixing_length_squared, shear_squared) - kvf = max( eddy_flux_coefficient_minimum_value, kvn * fofri ) + kvf = max( minimum_eddy_flux_coefficient, kvn * fofri ) end function calc_eddy_flux_coefficient pure elemental function calc_free_eddy_flux_coefficient(mixing_length_squared, & From a74defee441529f1fad352dc22cba6bffca3ffaf Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 27 Jan 2025 09:55:36 -0700 Subject: [PATCH 25/29] Updating free function to free_atm. --- phys_utils/atmos_pbl_utils.F90 | 9 ++++++--- .../tests/phys_utils/test_atmos_pbl_utils.pf | 12 ++++++------ 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/phys_utils/atmos_pbl_utils.F90 b/phys_utils/atmos_pbl_utils.F90 index fd6867e6..428c8c90 100644 --- a/phys_utils/atmos_pbl_utils.F90 +++ b/phys_utils/atmos_pbl_utils.F90 @@ -14,7 +14,7 @@ module atmos_phys_pbl_utils public :: calc_obukhov_length public :: calc_virtual_temperature public :: calc_eddy_flux_coefficient - public :: calc_free_eddy_flux_coefficient + public :: calc_free_atm_eddy_flux_coefficient real(kind_phys), parameter :: minimum_friction_velocity = 0.01_kind_phys real(kind_phys), parameter :: minimum_eddy_flux_coefficient = 0.01_kind_phys ! CCM1 2.f.14 @@ -84,6 +84,7 @@ pure elemental function calc_obukhov_length(thvs, ustar, g, vk, kbfs) result(obu ! 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{ real(kind_phys), intent(in) :: thvs ! virtual potential temperature at surface real(kind_phys), intent(in) :: ustar ! Surface friction velocity [ m/s ] @@ -142,7 +143,7 @@ pure elemental function calc_eddy_flux_coefficient(mixing_length_squared, & kvf = max( minimum_eddy_flux_coefficient, kvn * fofri ) end function calc_eddy_flux_coefficient - pure elemental function calc_free_eddy_flux_coefficient(mixing_length_squared, & + pure elemental function calc_free_atm_eddy_flux_coefficient(mixing_length_squared, & richardson_number, & shear_squared) & result(kvf) @@ -164,13 +165,14 @@ pure elemental function calc_free_eddy_flux_coefficient(mixing_length_squared, & kvn = neutral_exchange_coefficient(mixing_length_squared, shear_squared) kvf = kvn * fofri end if - end function calc_free_eddy_flux_coefficient + 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 @@ -184,6 +186,7 @@ pure elemental function stable_gradient_richardson_stability_parameter(richardso ! 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 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 index b6530b4b..0f40ac59 100644 --- a/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf +++ b/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf @@ -1,5 +1,5 @@ @test -subroutine test_free_eddy_coef_is_zero_when_richardson_equals_zero() +subroutine test_free_eddy_coef_is_zero_when_ri_equals_zero() use funit use atmos_phys_pbl_utils, only : calc_free_eddy_flux_coefficient use ccpp_kinds, only: kind_phys @@ -9,17 +9,17 @@ subroutine test_free_eddy_coef_is_zero_when_richardson_equals_zero() kvf = calc_free_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_richardson_equals_zero +end subroutine test_free_eddy_coef_is_zero_when_ri_equals_zero @test -subroutine test_free_eddy_coef_is_zero_when_richardson_stable_near_zero() +subroutine test_free_eddy_atm_coef_is_zero_when_ri_stable_near_zero() use funit - use atmos_phys_pbl_utils, only : calc_free_eddy_flux_coefficient + 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_eddy_flux_coefficient(30.0_kind_phys, nearest(0.0_kind_phys, 1.0_kind_phys), 0.01_kind_phys) + 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_coef_is_zero_when_richardson_stable_near_zero +end subroutine test_free_eddy_atm_coef_is_zero_when_ri_stable_near_zero From 7a70379f881deab0479c8c0630274589864cfc73 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 27 Jan 2025 10:10:17 -0700 Subject: [PATCH 26/29] Fixing function name under test. --- test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 index 0f40ac59..d7bb2683 100644 --- a/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf +++ b/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf @@ -1,12 +1,12 @@ @test subroutine test_free_eddy_coef_is_zero_when_ri_equals_zero() use funit - use atmos_phys_pbl_utils, only : calc_free_eddy_flux_coefficient + 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_eddy_flux_coefficient(30.0_kind_phys, 0.0_kind_phys, 0.01_kind_phys) + 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 From 1d0801a8da920d567636d5653b2e7d99b5c86bee Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 27 Jan 2025 10:20:19 -0700 Subject: [PATCH 27/29] adding latex form of obukhov length --- phys_utils/atmos_pbl_utils.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phys_utils/atmos_pbl_utils.F90 b/phys_utils/atmos_pbl_utils.F90 index 428c8c90..743bc43d 100644 --- a/phys_utils/atmos_pbl_utils.F90 +++ b/phys_utils/atmos_pbl_utils.F90 @@ -84,7 +84,7 @@ pure elemental function calc_obukhov_length(thvs, ustar, g, vk, kbfs) result(obu ! 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{ + ! \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 ] From 330262591a964d09e4d22bad462e767f0950b2dc Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Mon, 27 Jan 2025 14:33:26 -0700 Subject: [PATCH 28/29] Adding comment on unknown variable. --- phys_utils/atmos_pbl_utils.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phys_utils/atmos_pbl_utils.F90 b/phys_utils/atmos_pbl_utils.F90 index 743bc43d..f3626d95 100644 --- a/phys_utils/atmos_pbl_utils.F90 +++ b/phys_utils/atmos_pbl_utils.F90 @@ -16,7 +16,7 @@ module atmos_phys_pbl_utils public :: calc_eddy_flux_coefficient public :: calc_free_atm_eddy_flux_coefficient - real(kind_phys), parameter :: minimum_friction_velocity = 0.01_kind_phys + 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 From e76c482d392a9bfd163925fc0c006c23c372b4e9 Mon Sep 17 00:00:00 2001 From: Michael Waxmonsky Date: Tue, 28 Jan 2025 16:15:52 -0700 Subject: [PATCH 29/29] Fixing whitespace --- phys_utils/atmos_pbl_utils.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/phys_utils/atmos_pbl_utils.F90 b/phys_utils/atmos_pbl_utils.F90 index f3626d95..24104e76 100644 --- a/phys_utils/atmos_pbl_utils.F90 +++ b/phys_utils/atmos_pbl_utils.F90 @@ -144,9 +144,9 @@ pure elemental function calc_eddy_flux_coefficient(mixing_length_squared, & end function calc_eddy_flux_coefficient pure elemental function calc_free_atm_eddy_flux_coefficient(mixing_length_squared, & - richardson_number, & - shear_squared) & - result(kvf) + 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