Line data Source code
1 : module atmos_phys_pbl_utils
2 : ! Planetary boundary layer related functions used for vertical diffusion schemes.
3 :
4 : use ccpp_kinds, only: kind_phys
5 :
6 : implicit none
7 : private
8 :
9 : public :: calc_ideal_gas_rrho
10 : public :: calc_friction_velocity
11 : public :: calc_kinematic_heat_flux
12 : public :: calc_kinematic_water_vapor_flux
13 : public :: calc_kinematic_buoyancy_flux
14 : public :: calc_obukhov_length
15 : public :: calc_virtual_temperature
16 : public :: calc_eddy_flux_coefficient
17 : public :: calc_free_atm_eddy_flux_coefficient
18 :
19 : real(kind_phys), parameter :: minimum_friction_velocity = 0.01_kind_phys ! Assuming minimum for coarse grids
20 : real(kind_phys), parameter :: minimum_eddy_flux_coefficient = 0.01_kind_phys ! CCM1 2.f.14
21 :
22 : contains
23 :
24 6463800 : pure elemental function calc_ideal_gas_rrho(rair, surface_temperature, pmid) result(rrho)
25 : ! air density reciprocal
26 : ! Taken from https://glossary.ametsoc.org/wiki/Equation_of_state
27 : ! where \alpha = rrho
28 : real(kind_phys), intent(in) :: rair ! gas constant for dry air [ J kg-1 K-1 ]
29 : real(kind_phys), intent(in) :: surface_temperature ! [ K ]
30 : real(kind_phys), intent(in) :: pmid ! midpoint pressure (bottom level) [ Pa ]
31 :
32 : real(kind_phys) :: rrho ! 1./bottom level density [ m3 kg-1 ]
33 :
34 6463800 : rrho = rair * surface_temperature / pmid
35 6463800 : end function calc_ideal_gas_rrho
36 :
37 6463800 : pure elemental function calc_friction_velocity(taux, tauy, rrho) result(friction_velocity)
38 : ! https://glossary.ametsoc.org/wiki/Friction_velocity
39 : ! NOTE: taux,tauy come form the expansion of the Reynolds stress
40 : !
41 : ! Also found in:
42 : ! Stull, Roland B. An Introduction to Boundary Layer Meteorology. Springer Kluwer Academic Publishers, 1988. Print.
43 : ! DOI: https://doi.org/10.1007/978-94-009-3027-8
44 : ! Equation 2.10b, page 67
45 :
46 : real(kind_phys), intent(in) :: taux ! surface u stress [ N m-2 ]
47 : real(kind_phys), intent(in) :: tauy ! surface v stress [ N m-2 ]
48 : real(kind_phys), intent(in) :: rrho ! 1./bottom level density [ m3 kg-1 ]
49 :
50 : real(kind_phys) :: friction_velocity ! surface friction velocity [ m s-1 ]
51 :
52 6463800 : friction_velocity = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), minimum_friction_velocity )
53 6463800 : end function calc_friction_velocity
54 :
55 3693600 : pure elemental function calc_kinematic_heat_flux(shflx, rrho, cpair) result(khfs)
56 : real(kind_phys), intent(in) :: shflx ! surface heat flux [ W m-2 ]
57 : real(kind_phys), intent(in) :: rrho ! 1./bottom level density [ m3 kg-1 ]
58 : real(kind_phys), intent(in) :: cpair ! specific heat of dry air [ J kg-1 K-1 ]
59 :
60 : real(kind_phys) :: khfs ! surface kinematic heat flux [ m K s-1 ]
61 :
62 3693600 : khfs = shflx*rrho/cpair
63 3693600 : end function calc_kinematic_heat_flux
64 :
65 3693600 : pure elemental function calc_kinematic_water_vapor_flux(qflx, rrho) result(kqfs)
66 : real(kind_phys), intent(in) :: qflx ! water vapor flux [ kg m-2 s-1 ]
67 : real(kind_phys), intent(in) :: rrho ! 1./bottom level density [ m3 kg-1 ]
68 :
69 : real(kind_phys) :: kqfs ! surface kinematic water vapor flux [ kg kg-1 m s-1 ]
70 :
71 3693600 : kqfs = qflx*rrho
72 3693600 : end function calc_kinematic_water_vapor_flux
73 :
74 3693600 : pure elemental function calc_kinematic_buoyancy_flux(khfs, zvir, ths, kqfs) result(kbfs)
75 : real(kind_phys), intent(in) :: khfs ! surface kinematic heat flux [ m K s-1 ]
76 : real(kind_phys), intent(in) :: zvir ! rh2o/rair - 1
77 : real(kind_phys), intent(in) :: ths ! potential temperature at surface [ K ]
78 : real(kind_phys), intent(in) :: kqfs ! surface kinematic water vapor flux [ kg kg-1 m s-1 ]
79 :
80 : ! (`kbfs = \overline{(w' \theta'_v)}_s`)
81 : real(kind_phys) :: kbfs ! surface kinematic buoyancy flux [ m K s-1 ]
82 :
83 3693600 : kbfs = khfs + zvir*ths*kqfs
84 3693600 : end function calc_kinematic_buoyancy_flux
85 :
86 3693600 : pure elemental function calc_obukhov_length(thvs, ustar, g, karman, kbfs) result(obukhov_length)
87 : ! Stull, Roland B. An Introduction to Boundary Layer Meteorology. Springer Kluwer Academic Publishers, 1988. Print.
88 : ! DOI: https://doi.org/10.1007/978-94-009-3027-8
89 : ! Equation 5.7c, page 181
90 : ! \frac{-\theta*u_*^3}{g*k*\overline{(w' \theta_v')}_s} = frac{-\theta*u_*^3}{g*k*kbfs}
91 :
92 : real(kind_phys), intent(in) :: thvs ! virtual potential temperature at surface [ K ]
93 : real(kind_phys), intent(in) :: ustar ! Surface friction velocity [ m s-1 ]
94 : real(kind_phys), intent(in) :: g ! acceleration of gravity [ m s-2 ]
95 : real(kind_phys), intent(in) :: karman ! Von Karman's constant (unitless)
96 : real(kind_phys), intent(in) :: kbfs ! surface kinematic buoyancy flux [ m K s-1 ]
97 :
98 : real(kind_phys) :: obukhov_length ! Obukhov length [ m ]
99 :
100 : ! Added sign(...) term to prevent division by 0 and using the fact that
101 : ! `kbfs = \overline{(w' \theta_v')}_s`
102 : obukhov_length = -thvs * ustar**3 / &
103 3693600 : (g*karman*(kbfs + sign(1.e-10_kind_phys,kbfs)))
104 3693600 : end function calc_obukhov_length
105 :
106 85876200 : pure elemental function calc_virtual_temperature(temperature, specific_humidity, zvir) result(virtual_temperature)
107 : ! Williamson, D., Kiehl, J., Ramanathan, V., Dickinson, R., & Hack, J. (1987).
108 : ! Description of the NCAR Community Climate Model (CCM1).
109 : ! University Corporation for Atmospheric Research. https://doi.org/10.5065/D6TB14WH (Original work published 1987)
110 : ! Equation 2.a.7
111 :
112 : real(kind_phys), intent(in) :: temperature
113 : real(kind_phys), intent(in) :: specific_humidity
114 : real(kind_phys), intent(in) :: zvir ! rh2o/rair - 1
115 :
116 : real(kind_phys) :: virtual_temperature
117 :
118 85876200 : virtual_temperature = temperature * (1.0_kind_phys + zvir*specific_humidity)
119 85876200 : end function calc_virtual_temperature
120 :
121 0 : pure elemental function calc_eddy_flux_coefficient(mixing_length_squared, &
122 : richardson_number, &
123 : shear_squared) &
124 : result(kvf)
125 : ! Computes exchange coefficients for turbulent flows.
126 : !
127 : ! The stable case (Richardson number, Ri>0) is taken from Holtslag and
128 : ! Beljaars (1989), ECMWF proceedings.
129 : ! The unstable case (Richardson number, Ri<0) is taken from CCM1.
130 :
131 : real(kind_phys), intent(in) :: mixing_length_squared ! [ m2 ]
132 : real(kind_phys), intent(in) :: richardson_number ! [ 1 ]
133 : real(kind_phys), intent(in) :: shear_squared ! [ s-2 ]
134 :
135 : real(kind_phys) :: kvf ! Eddy diffusivity ! [ m2 s-1 ]
136 :
137 : real(kind_phys) :: fofri ! f(ri)
138 : real(kind_phys) :: kvn ! Neutral Kv
139 :
140 0 : if( richardson_number < 0.0_kind_phys ) then
141 0 : fofri = unstable_gradient_richardson_stability_parameter(richardson_number)
142 : else
143 0 : fofri = stable_gradient_richardson_stability_parameter(richardson_number)
144 : end if
145 0 : kvn = neutral_exchange_coefficient(mixing_length_squared, shear_squared)
146 0 : kvf = max( minimum_eddy_flux_coefficient, kvn * fofri )
147 0 : end function calc_eddy_flux_coefficient
148 :
149 84952800 : pure elemental function calc_free_atm_eddy_flux_coefficient(mixing_length_squared, &
150 : richardson_number, &
151 : shear_squared) &
152 : result(kvf)
153 : ! same as austausch_atm but only mixing for Ri<0
154 : ! i.e. no background mixing and mixing for Ri>0
155 :
156 : real(kind_phys), intent(in) :: mixing_length_squared ! [ m2 ]
157 : real(kind_phys), intent(in) :: richardson_number ! [ 1 ]
158 : real(kind_phys), intent(in) :: shear_squared ! [ s-2 ]
159 :
160 : real(kind_phys) :: kvf ! [ m2 s-1 ]
161 :
162 : real(kind_phys) :: fofri ! f(ri)
163 : real(kind_phys) :: kvn ! Neutral Kv
164 :
165 84952800 : kvf = 0.0_kind_phys
166 84952800 : if( richardson_number < 0.0_kind_phys ) then
167 1947345 : fofri = unstable_gradient_richardson_stability_parameter(richardson_number)
168 1947345 : kvn = neutral_exchange_coefficient(mixing_length_squared, shear_squared)
169 1947345 : kvf = kvn * fofri
170 : end if
171 84952800 : end function calc_free_atm_eddy_flux_coefficient
172 :
173 1947345 : pure elemental function unstable_gradient_richardson_stability_parameter(richardson_number) result(modifier)
174 : ! Williamson, D., Kiehl, J., Ramanathan, V., Dickinson, R., & Hack, J. (1987).
175 : ! Description of the NCAR Community Climate Model (CCM1).
176 : ! University Corporation for Atmospheric Research. https://doi.org/10.5065/D6TB14WH (Original work published 1987)
177 : ! Equation 2.f.13
178 : ! \sqrt{ 1-18*Ri }
179 :
180 : real(kind_phys), intent(in) :: richardson_number
181 :
182 : real(kind_phys) :: modifier
183 :
184 1947345 : modifier = sqrt( 1._kind_phys - 18._kind_phys * richardson_number )
185 1947345 : end function unstable_gradient_richardson_stability_parameter
186 :
187 0 : pure elemental function stable_gradient_richardson_stability_parameter(richardson_number) result(modifier)
188 : ! Holtslag, A. A. M., and Beljaars A. C. M. , 1989: Surface flux parameterization schemes: Developments and experiences at KNMI.
189 : ! ECMWF Workshop on Parameterization of Fluxes and Land Surface, Reading, United Kingdom, ECMWF, 121–147.
190 : ! equation 20, page 140
191 : ! Originally used published equation from CCM1, 2.f.12, page 11
192 : ! \frac{1}{1+10*Ri*(1+8*Ri)}
193 :
194 : real(kind_phys), intent(in) :: richardson_number
195 :
196 : real(kind_phys) :: modifier
197 :
198 : modifier = 1.0_kind_phys / &
199 0 : ( 1.0_kind_phys + 10.0_kind_phys * richardson_number * ( 1.0_kind_phys + 8.0_kind_phys * richardson_number ) )
200 0 : end function stable_gradient_richardson_stability_parameter
201 :
202 1947345 : pure elemental function neutral_exchange_coefficient(mixing_length_squared, shear_squared) result(neutral_k)
203 : ! Williamson, D., Kiehl, J., Ramanathan, V., Dickinson, R., & Hack, J. (1987).
204 : ! Description of the NCAR Community Climate Model (CCM1).
205 : ! University Corporation for Atmospheric Research. https://doi.org/10.5065/D6TB14WH (Original work published 1987)
206 : ! Equation 2.f.15, page 12
207 : ! NOTE: shear_squared variable currently (01/2025) computed in hb_diff.F90 (s2 in trbintd(...)) matches referenced equation.
208 :
209 : real(kind_phys), intent(in) :: mixing_length_squared ! [ m2 ]
210 : real(kind_phys), intent(in) :: shear_squared ! [ s-2 ]
211 :
212 : real(kind_phys) :: neutral_k
213 :
214 1947345 : neutral_k = mixing_length_squared * sqrt(shear_squared)
215 1947345 : end function neutral_exchange_coefficient
216 : end module atmos_phys_pbl_utils
|
