Line data Source code
1 : module pbl_utils
2 : !-----------------------------------------------------------------------!
3 : ! Module to hold PBL-related subprograms that may be used with multiple !
4 : ! different vertical diffusion schemes. !
5 : ! !
6 : ! Public subroutines: !
7 : !
8 : ! calc_obklen !
9 : ! !
10 : !------------------ History --------------------------------------------!
11 : ! Created: Apr. 2012, by S. Santos !
12 : !-----------------------------------------------------------------------!
13 :
14 : use shr_kind_mod, only: r8 => shr_kind_r8
15 :
16 : implicit none
17 : private
18 :
19 : ! Public Procedures
20 : !----------------------------------------------------------------------!
21 : ! Excepting the initialization procedure, these are elemental
22 : ! procedures, so they can accept scalars or any dimension of array as
23 : ! arguments, as long as all arguments have the same number of
24 : ! elements.
25 : public pbl_utils_init
26 : public calc_ustar
27 : public calc_obklen
28 : public virtem
29 : public compute_radf
30 : public austausch_atm, austausch_atm_free
31 :
32 : real(r8), parameter :: ustar_min = 0.01_r8
33 :
34 : real(r8) :: g ! acceleration of gravity
35 : real(r8) :: vk ! Von Karman's constant
36 : real(r8) :: cpair ! specific heat of dry air
37 : real(r8) :: rair ! gas constant for dry air
38 : real(r8) :: zvir ! rh2o/rair - 1
39 :
40 :
41 : !------------------------------------------------------------------------!
42 : ! Purpose: Compilers aren't creating optimized vector versions of !
43 : ! elemental routines, so we'll explicitly create them and bind !
44 : ! them via an interface for transparent use !
45 : !------------------------------------------------------------------------!
46 : interface calc_ustar
47 : module procedure calc_ustar_scalar
48 : module procedure calc_ustar_vector
49 : end interface
50 :
51 : interface calc_obklen
52 : module procedure calc_obklen_scalar
53 : module procedure calc_obklen_vector
54 : end interface
55 :
56 : interface virtem
57 : module procedure virtem_vector1D
58 : module procedure virtem_vector2D ! Used in hb_diff.F90
59 : end interface
60 :
61 :
62 :
63 : contains
64 :
65 1536 : subroutine pbl_utils_init(g_in,vk_in,cpair_in,rair_in,zvir_in)
66 :
67 : !-----------------------------------------------------------------------!
68 : ! Purpose: Set constants to be used in calls to later functions !
69 : !-----------------------------------------------------------------------!
70 :
71 : real(r8), intent(in) :: g_in ! acceleration of gravity
72 : real(r8), intent(in) :: vk_in ! Von Karman's constant
73 : real(r8), intent(in) :: cpair_in ! specific heat of dry air
74 : real(r8), intent(in) :: rair_in ! gas constant for dry air
75 : real(r8), intent(in) :: zvir_in ! rh2o/rair - 1
76 :
77 1536 : g = g_in
78 1536 : vk = vk_in
79 1536 : cpair = cpair_in
80 1536 : rair = rair_in
81 1536 : zvir = zvir_in
82 :
83 1536 : end subroutine pbl_utils_init
84 :
85 0 : subroutine calc_ustar_scalar( t, pmid, taux, tauy, &
86 : rrho, ustar)
87 :
88 : !-----------------------------------------------------------------------!
89 : ! Purpose: Calculate ustar and bottom level density (necessary for !
90 : ! Obukhov length calculation). !
91 : !-----------------------------------------------------------------------!
92 :
93 : real(r8), intent(in) :: t ! surface temperature
94 : real(r8), intent(in) :: pmid ! midpoint pressure (bottom level)
95 : real(r8), intent(in) :: taux ! surface u stress [N/m2]
96 : real(r8), intent(in) :: tauy ! surface v stress [N/m2]
97 :
98 : real(r8), intent(out) :: rrho ! 1./bottom level density
99 : real(r8), intent(out) :: ustar ! surface friction velocity [m/s]
100 :
101 0 : rrho = rair * t / pmid
102 0 : ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min )
103 :
104 0 : end subroutine calc_ustar_scalar
105 :
106 1489176 : subroutine calc_ustar_vector(n, t, pmid, taux, tauy, &
107 1489176 : rrho, ustar)
108 :
109 : !-----------------------------------------------------------------------!
110 : ! Purpose: Calculate ustar and bottom level density (necessary for !
111 : ! Obukhov length calculation). !
112 : !-----------------------------------------------------------------------!
113 : integer, intent(in) :: n ! Length of vectors
114 :
115 : real(r8), intent(in) :: t(n) ! surface temperature
116 : real(r8), intent(in) :: pmid(n) ! midpoint pressure (bottom level)
117 : real(r8), intent(in) :: taux(n) ! surface u stress [N/m2]
118 : real(r8), intent(in) :: tauy(n) ! surface v stress [N/m2]
119 :
120 :
121 : real(r8), intent(out) :: rrho(n) ! 1./bottom level density
122 : real(r8), intent(out) :: ustar(n) ! surface friction velocity [m/s]
123 :
124 :
125 24865776 : rrho = rair * t / pmid
126 24865776 : ustar = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), ustar_min )
127 :
128 1489176 : end subroutine calc_ustar_vector
129 :
130 0 : subroutine calc_obklen_scalar( ths, thvs, qflx, shflx, rrho, ustar, &
131 : khfs, kqfs, kbfs, obklen)
132 :
133 : !-----------------------------------------------------------------------!
134 : ! Purpose: Calculate Obukhov length and kinematic fluxes. !
135 : !-----------------------------------------------------------------------!
136 :
137 : real(r8), intent(in) :: ths ! potential temperature at surface [K]
138 : real(r8), intent(in) :: thvs ! virtual potential temperature at surface
139 : real(r8), intent(in) :: qflx ! water vapor flux (kg/m2/s)
140 : real(r8), intent(in) :: shflx ! surface heat flux (W/m2)
141 :
142 : real(r8), intent(in) :: rrho ! 1./bottom level density [ m3/kg ]
143 : real(r8), intent(in) :: ustar ! Surface friction velocity [ m/s ]
144 :
145 : real(r8), intent(out) :: khfs ! sfc kinematic heat flux [mK/s]
146 : real(r8), intent(out) :: kqfs ! sfc kinematic water vapor flux [m/s]
147 : real(r8), intent(out) :: kbfs ! sfc kinematic buoyancy flux [m^2/s^3]
148 : real(r8), intent(out) :: obklen ! Obukhov length
149 :
150 : ! Need kinematic fluxes for Obukhov:
151 0 : khfs = shflx*rrho/cpair
152 0 : kqfs = qflx*rrho
153 0 : kbfs = khfs + zvir*ths*kqfs
154 :
155 : ! Compute Obukhov length:
156 0 : obklen = -thvs * ustar**3 / (g*vk*(kbfs + sign(1.e-10_r8,kbfs)))
157 :
158 0 : end subroutine calc_obklen_scalar
159 :
160 1489176 : subroutine calc_obklen_vector(n, ths, thvs, qflx, shflx, rrho, ustar, &
161 1489176 : khfs, kqfs, kbfs, obklen)
162 :
163 : !-----------------------------------------------------------------------!
164 : ! Purpose: Calculate Obukhov length and kinematic fluxes. !
165 : !-----------------------------------------------------------------------!
166 : integer, intent(in) :: n ! Length of vectors
167 :
168 : real(r8), intent(in) :: ths(n) ! potential temperature at surface [K]
169 : real(r8), intent(in) :: thvs(n) ! virtual potential temperature at surface
170 : real(r8), intent(in) :: qflx(n) ! water vapor flux (kg/m2/s)
171 : real(r8), intent(in) :: shflx(n) ! surface heat flux (W/m2)
172 :
173 : real(r8), intent(in) :: rrho(n) ! 1./bottom level density [ m3/kg ]
174 : real(r8), intent(in) :: ustar(n) ! Surface friction velocity [ m/s ]
175 :
176 : real(r8), intent(out) :: khfs(n) ! sfc kinematic heat flux [mK/s]
177 : real(r8), intent(out) :: kqfs(n) ! sfc kinematic water vapor flux [m/s]
178 : real(r8), intent(out) :: kbfs(n) ! sfc kinematic buoyancy flux [m^2/s^3]
179 : real(r8), intent(out) :: obklen(n) ! Obukhov length
180 :
181 :
182 : ! Need kinematic fluxes for Obukhov:
183 24865776 : khfs = shflx*rrho/cpair
184 24865776 : kqfs = qflx*rrho
185 24865776 : kbfs = khfs + zvir*ths*kqfs
186 :
187 : ! Compute Obukhov length:
188 24865776 : obklen = -thvs * ustar**3 / (g*vk*(kbfs + sign(1.e-10_r8,kbfs)))
189 :
190 1489176 : end subroutine calc_obklen_vector
191 :
192 0 : subroutine virtem_vector1D(n, t,q, virtem)
193 :
194 : !-----------------------------------------------------------------------!
195 : ! Purpose: Calculate virtual temperature from temperature and specific !
196 : ! humidity. !
197 : !-----------------------------------------------------------------------!
198 :
199 : integer, intent(in) :: n ! vector length
200 :
201 : real(r8), intent(in) :: t(n), q(n)
202 : real(r8), intent(out):: virtem(n)
203 :
204 0 : virtem = t * (1.0_r8 + zvir*q)
205 :
206 0 : end subroutine virtem_vector1D
207 :
208 1489176 : subroutine virtem_vector2D(n, m, t, q, virtem)
209 :
210 : !-----------------------------------------------------------------------!
211 : ! Purpose: Calculate virtual temperature from temperature and specific !
212 : ! humidity. !
213 : !-----------------------------------------------------------------------!
214 :
215 : integer, intent(in) :: n, m ! vector lengths
216 :
217 : real(r8), intent(in) :: t(n,m), q(n,m)
218 : real(r8), intent(out):: virtem(n,m)
219 :
220 647999352 : virtem = t * (1.0_r8 + zvir*q)
221 :
222 1489176 : end subroutine virtem_vector2D
223 :
224 :
225 0 : subroutine compute_radf( choice_radf, i, pcols, pver, ncvmax, ncvfin, ktop, qmin, &
226 0 : ql, pi, qrlw, g, cldeff, zi, chs, lwp_CL, opt_depth_CL, &
227 0 : radinvfrac_CL, radf_CL )
228 : ! -------------------------------------------------------------------------- !
229 : ! Purpose: !
230 : ! Calculate cloud-top radiative cooling contribution to buoyancy production. !
231 : ! Here, 'radf' [m2/s3] is additional buoyancy flux at the CL top interface !
232 : ! associated with cloud-top LW cooling being mainly concentrated near the CL !
233 : ! top interface ( just below CL top interface ). Contribution of SW heating !
234 : ! within the cloud is not included in this radiative buoyancy production !
235 : ! since SW heating is more broadly distributed throughout the CL top layer. !
236 : ! -------------------------------------------------------------------------- !
237 :
238 : !-----------------!
239 : ! Input variables !
240 : !-----------------!
241 : character(len=6), intent(in) :: choice_radf ! Method for calculating radf
242 : integer, intent(in) :: i ! Index of current column
243 : integer, intent(in) :: pcols ! Number of atmospheric columns
244 : integer, intent(in) :: pver ! Number of atmospheric layers
245 : integer, intent(in) :: ncvmax ! Max numbers of CLs (perhaps equal to pver)
246 : integer, intent(in) :: ncvfin(pcols) ! Total number of CL in column
247 : integer, intent(in) :: ktop(pcols, ncvmax) ! ktop for current column
248 : real(r8), intent(in) :: qmin ! Minimum grid-mean LWC counted as clouds [kg/kg]
249 : real(r8), intent(in) :: ql(pcols, pver) ! Liquid water specific humidity [ kg/kg ]
250 : real(r8), intent(in) :: pi(pcols, pver+1) ! Interface pressures [ Pa ]
251 : real(r8), intent(in) :: qrlw(pcols, pver) ! Input grid-mean LW heating rate : [ K/s ] * cpair * dp = [ W/kg*Pa ]
252 : real(r8), intent(in) :: g ! Gravitational acceleration
253 : real(r8), intent(in) :: cldeff(pcols,pver) ! Effective Cloud Fraction [fraction]
254 : real(r8), intent(in) :: zi(pcols, pver+1) ! Interface heights [ m ]
255 : real(r8), intent(in) :: chs(pcols, pver+1) ! Buoyancy coeffi. saturated sl (heat) coef. at all interfaces.
256 :
257 : !------------------!
258 : ! Output variables !
259 : !------------------!
260 : real(r8), intent(out) :: lwp_CL(ncvmax) ! LWP in the CL top layer [ kg/m2 ]
261 : real(r8), intent(out) :: opt_depth_CL(ncvmax) ! Optical depth of the CL top layer
262 : real(r8), intent(out) :: radinvfrac_CL(ncvmax) ! Fraction of LW radiative cooling confined in the top portion of CL
263 : real(r8), intent(out) :: radf_CL(ncvmax) ! Buoyancy production at the CL top due to radiative cooling [ m2/s3 ]
264 :
265 : !-----------------!
266 : ! Local variables !
267 : !-----------------!
268 : integer :: kt, ncv
269 : real(r8) :: lwp, opt_depth, radinvfrac, radf
270 :
271 :
272 : !-----------------!
273 : ! Begin main code !
274 : !-----------------!
275 0 : lwp_CL = 0._r8
276 0 : opt_depth_CL = 0._r8
277 0 : radinvfrac_CL = 0._r8
278 0 : radf_CL = 0._r8
279 :
280 : ! ---------------------------------------- !
281 : ! Perform do loop for individual CL regime !
282 : ! ---------------------------------------- !
283 0 : do ncv = 1, ncvfin(i)
284 0 : kt = ktop(i,ncv)
285 : !-----------------------------------------------------!
286 : ! Compute radf for each CL regime and for each column !
287 : !-----------------------------------------------------!
288 0 : if( choice_radf .eq. 'orig' ) then
289 0 : if( ql(i,kt) .gt. qmin .and. ql(i,kt-1) .lt. qmin ) then
290 0 : lwp = ql(i,kt) * ( pi(i,kt+1) - pi(i,kt) ) / g
291 0 : opt_depth = 156._r8 * lwp ! Estimated LW optical depth in the CL top layer
292 : ! Approximate LW cooling fraction concentrated at the inversion by using
293 : ! polynomial approx to exact formula 1-2/opt_depth+2/(exp(opt_depth)-1))
294 :
295 0 : radinvfrac = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth ) + opt_depth**2 )
296 0 : radf = qrlw(i,kt) / ( pi(i,kt) - pi(i,kt+1) ) ! Cp*radiative cooling = [ W/kg ]
297 0 : radf = max( radinvfrac * radf * ( zi(i,kt) - zi(i,kt+1) ), 0._r8 ) * chs(i,kt)
298 : ! We can disable cloud LW cooling contribution to turbulence by uncommenting:
299 : ! radf = 0._r8
300 : end if
301 :
302 0 : elseif( choice_radf .eq. 'ramp' ) then
303 :
304 0 : lwp = ql(i,kt) * ( pi(i,kt+1) - pi(i,kt) ) / g
305 0 : opt_depth = 156._r8 * lwp ! Estimated LW optical depth in the CL top layer
306 0 : radinvfrac = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth ) + opt_depth**2 )
307 0 : radinvfrac = max(cldeff(i,kt)-cldeff(i,kt-1),0._r8) * radinvfrac
308 0 : radf = qrlw(i,kt) / ( pi(i,kt) - pi(i,kt+1) ) ! Cp*radiative cooling [W/kg]
309 0 : radf = max( radinvfrac * radf * ( zi(i,kt) - zi(i,kt+1) ), 0._r8 ) * chs(i,kt)
310 :
311 0 : elseif( choice_radf .eq. 'maxi' ) then
312 :
313 : ! Radiative flux divergence both in 'kt' and 'kt-1' layers are included
314 : ! 1. From 'kt' layer
315 0 : lwp = ql(i,kt) * ( pi(i,kt+1) - pi(i,kt) ) / g
316 0 : opt_depth = 156._r8 * lwp ! Estimated LW optical depth in the CL top layer
317 0 : radinvfrac = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth ) + opt_depth**2 )
318 0 : radf = max( radinvfrac * qrlw(i,kt) / ( pi(i,kt) - pi(i,kt+1) ) * ( zi(i,kt) - zi(i,kt+1) ), 0._r8 )
319 : ! 2. From 'kt-1' layer and add the contribution from 'kt' layer
320 0 : lwp = ql(i,kt-1) * ( pi(i,kt) - pi(i,kt-1) ) / g
321 0 : opt_depth = 156._r8 * lwp ! Estimated LW optical depth in the CL top layer
322 0 : radinvfrac = opt_depth * ( 4._r8 + opt_depth ) / ( 6._r8 * ( 4._r8 + opt_depth) + opt_depth**2 )
323 0 : radf = radf + max( radinvfrac * qrlw(i,kt-1) / ( pi(i,kt-1) - pi(i,kt) ) * ( zi(i,kt-1) - zi(i,kt) ), 0._r8 )
324 0 : radf = max( radf, 0._r8 ) * chs(i,kt)
325 :
326 : endif
327 :
328 0 : lwp_CL(ncv) = lwp
329 0 : opt_depth_CL(ncv) = opt_depth
330 0 : radinvfrac_CL(ncv) = radinvfrac
331 0 : radf_CL(ncv) = radf
332 : end do ! ncv = 1, ncvfin(i)
333 0 : end subroutine compute_radf
334 :
335 1489176 : subroutine austausch_atm(pcols, ncol, pver, ntop, nbot, ml2, ri, s2, kvf)
336 :
337 : !---------------------------------------------------------------------- !
338 : ! !
339 : ! Purpose: Computes exchange coefficients for free turbulent flows. !
340 : ! !
341 : ! Method: !
342 : ! !
343 : ! The free atmosphere diffusivities are based on standard mixing length !
344 : ! forms for the neutral diffusivity multiplied by functns of Richardson !
345 : ! number. K = l^2 * |dV/dz| * f(Ri). The same functions are used for !
346 : ! momentum, potential temperature, and constitutents. !
347 : ! !
348 : ! The stable Richardson num function (Ri>0) is taken from Holtslag and !
349 : ! Beljaars (1989), ECMWF proceedings. f = 1 / (1 + 10*Ri*(1 + 8*Ri)) !
350 : ! The unstable Richardson number function (Ri<0) is taken from CCM1. !
351 : ! f = sqrt(1 - 18*Ri) !
352 : ! !
353 : ! Author: B. Stevens (rewrite, August 2000) !
354 : ! !
355 : !---------------------------------------------------------------------- !
356 :
357 : ! --------------- !
358 : ! Input arguments !
359 : ! --------------- !
360 :
361 : integer, intent(in) :: pcols ! Atmospheric columns dimension size
362 : integer, intent(in) :: ncol ! Number of atmospheric columns
363 : integer, intent(in) :: pver ! Number of atmospheric layers
364 : integer, intent(in) :: ntop ! Top layer for calculation
365 : integer, intent(in) :: nbot ! Bottom layer for calculation
366 :
367 : real(r8), intent(in) :: ml2(pver+1) ! Mixing lengths squared
368 : real(r8), intent(in) :: s2(pcols,pver) ! Shear squared
369 : real(r8), intent(in) :: ri(pcols,pver) ! Richardson no
370 :
371 : ! ---------------- !
372 : ! Output arguments !
373 : ! ---------------- !
374 :
375 : real(r8), intent(out) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers
376 :
377 : ! --------------- !
378 : ! Local Variables !
379 : ! --------------- !
380 :
381 : real(r8) :: fofri ! f(ri)
382 : real(r8) :: kvn ! Neutral Kv
383 :
384 : integer :: i ! Longitude index
385 : integer :: k ! Vertical index
386 :
387 : real(r8), parameter :: zkmin = 0.01_r8 ! Minimum kneutral*f(ri).
388 :
389 : ! ----------------------- !
390 : ! Main Computation Begins !
391 : ! ----------------------- !
392 :
393 672865128 : kvf(:ncol,:) = 0.0_r8
394 :
395 : ! Compute the free atmosphere vertical diffusion coefficients: kvh = kvq = kvm.
396 :
397 38718576 : do k = ntop, nbot - 1
398 623133576 : do i = 1, ncol
399 584415000 : if( ri(i,k) < 0.0_r8 ) then
400 12022687 : fofri = sqrt( max( 1._r8 - 18._r8 * ri(i,k), 0._r8 ) )
401 : else
402 572392313 : fofri = 1.0_r8 / ( 1.0_r8 + 10.0_r8 * ri(i,k) * ( 1.0_r8 + 8.0_r8 * ri(i,k) ) )
403 : end if
404 584415000 : kvn = ml2(k) * sqrt(s2(i,k))
405 621644400 : kvf(i,k+1) = max( zkmin, kvn * fofri )
406 : end do
407 : end do
408 :
409 1489176 : end subroutine austausch_atm
410 :
411 0 : subroutine austausch_atm_free(pcols, ncol, pver, ntop, nbot, ml2, ri, s2, kvf)
412 :
413 : !---------------------------------------------------------------------- !
414 : ! !
415 : ! same as austausch_atm but only mixing for Ri<0 !
416 : ! i.e. no background mixing and mixing for Ri>0 !
417 : ! !
418 : !---------------------------------------------------------------------- !
419 :
420 : ! --------------- !
421 : ! Input arguments !
422 : ! --------------- !
423 :
424 : integer, intent(in) :: pcols ! Atmospheric columns dimension size
425 : integer, intent(in) :: ncol ! Number of atmospheric columns
426 : integer, intent(in) :: pver ! Number of atmospheric layers
427 : integer, intent(in) :: ntop ! Top layer for calculation
428 : integer, intent(in) :: nbot ! Bottom layer for calculation
429 :
430 : real(r8), intent(in) :: ml2(pver+1) ! Mixing lengths squared
431 : real(r8), intent(in) :: s2(pcols,pver) ! Shear squared
432 : real(r8), intent(in) :: ri(pcols,pver) ! Richardson no
433 :
434 : ! ---------------- !
435 : ! Output arguments !
436 : ! ---------------- !
437 :
438 : real(r8), intent(out) :: kvf(pcols,pver+1) ! Eddy diffusivity for heat and tracers
439 :
440 : ! --------------- !
441 : ! Local Variables !
442 : ! --------------- !
443 :
444 : real(r8) :: fofri ! f(ri)
445 : real(r8) :: kvn ! Neutral Kv
446 :
447 : integer :: i ! Longitude index
448 : integer :: k ! Vertical index
449 :
450 : ! ----------------------- !
451 : ! Main Computation Begins !
452 : ! ----------------------- !
453 :
454 0 : kvf(:ncol,:) = 0.0_r8
455 : ! Compute the free atmosphere vertical diffusion coefficients: kvh = kvq = kvm.
456 0 : do k = ntop, nbot - 1
457 0 : do i = 1, ncol
458 0 : if( ri(i,k) < 0.0_r8 ) then
459 0 : fofri = sqrt( max( 1._r8 - 18._r8 * ri(i,k), 0._r8 ) )
460 : else
461 : fofri = 0.0_r8
462 : end if
463 0 : kvn = ml2(k) * sqrt(s2(i,k))
464 0 : kvf(i,k+1) = kvn * fofri
465 : end do
466 : end do
467 0 : end subroutine austausch_atm_free
468 :
469 : end module pbl_utils
|
