Line data Source code
1 : module hb_diff
2 : !---------------------------------------------------------------------------------
3 : ! Module to compute mixing coefficients associated with turbulence in the
4 : ! planetary boundary layer and elsewhere. PBL coefficients are based on Holtslag
5 : ! and Boville, 1991.
6 : !
7 : ! Public interfaces:
8 : ! init_hb_diff initializes time independent coefficients
9 : ! compute_hb_diff computes eddy diffusivities and counter-gradient fluxes
10 : !
11 : ! Private methods:
12 : ! trbintd initializes time dependent variables
13 : ! pblintd initializes time dependent variables that depend pbl depth
14 : ! austausch_pbl computes pbl exchange coefficients
15 : !
16 : !---------------------------Code history--------------------------------
17 : ! Standardized: J. Rosinski, June 1992
18 : ! Reviewed: P. Rasch, B. Boville, August 1992
19 : ! Reviewed: P. Rasch, April 1996
20 : ! Reviewed: B. Boville, April 1996
21 : ! rewritten: B. Boville, May 2000
22 : ! rewritten: B. Stevens, August 2000
23 : ! modularized: J. McCaa, September 2004
24 : !---------------------------------------------------------------------------------
25 : use shr_kind_mod, only: r8 => shr_kind_r8
26 : use spmd_utils, only: masterproc ! output from hb_init should be eliminated
27 : use ppgrid, only: pver, pverp, pcols ! these should be passed in
28 : use cam_logfile, only: iulog
29 :
30 : implicit none
31 : private
32 : save
33 :
34 : ! Public interfaces
35 : public init_hb_diff
36 : public compute_hb_diff
37 : public compute_hb_free_atm_diff
38 : public pblintd
39 : !
40 : ! PBL limits
41 : !
42 : real(r8), parameter :: pblmaxp = 4.e4_r8 ! pbl max depth in pressure units
43 : real(r8), parameter :: zkmin = 0.01_r8 ! Minimum kneutral*f(ri)
44 : !
45 : ! PBL Parameters
46 : !
47 : real(r8), parameter :: onet = 1._r8/3._r8 ! 1/3 power in wind gradient expression
48 : real(r8), parameter :: betam = 15.0_r8 ! Constant in wind gradient expression
49 : real(r8), parameter :: betas = 5.0_r8 ! Constant in surface layer gradient expression
50 : real(r8), parameter :: betah = 15.0_r8 ! Constant in temperature gradient expression
51 : real(r8), parameter :: fakn = 7.2_r8 ! Constant in turbulent prandtl number
52 : real(r8), parameter :: fak = 8.5_r8 ! Constant in surface temperature excess
53 : real(r8), parameter :: ricr = 0.3_r8 ! Critical richardson number
54 : real(r8), parameter :: sffrac= 0.1_r8 ! Surface layer fraction of boundary layer
55 : real(r8), parameter :: binm = betam*sffrac ! betam * sffrac
56 : real(r8), parameter :: binh = betah*sffrac ! betah * sffrac
57 :
58 : ! Pbl constants set using values from other parts of code
59 :
60 : real(r8) :: cpair ! Specific heat of dry air
61 : real(r8) :: g ! Gravitational acceleration
62 : real(r8) :: ml2(pverp) ! Mixing lengths squared
63 : real(r8) :: vk ! Von Karman's constant
64 : real(r8) :: ccon ! fak * sffrac * vk
65 :
66 : integer :: npbl ! Maximum number of levels in pbl from surface
67 : integer :: ntop_turb ! Top level to which turbulent vertical diffusion is applied.
68 : integer :: nbot_turb ! Bottom level to which turbulent vertical diff is applied.
69 :
70 : !===============================================================================
71 : CONTAINS
72 : !===============================================================================
73 :
74 1536 : subroutine init_hb_diff(gravx, cpairx, ntop_eddy, nbot_eddy, pref_mid, &
75 1536 : vkx, eddy_scheme)
76 :
77 : !-----------------------------------------------------------------------
78 : !
79 : ! Initialize time independent variables of turbulence/pbl package.
80 : !
81 : !-----------------------------------------------------------------------
82 :
83 : !------------------------------Arguments--------------------------------
84 : real(r8), intent(in) :: gravx ! acceleration of gravity
85 : real(r8), intent(in) :: cpairx ! specific heat of dry air
86 : real(r8), intent(in) :: pref_mid(pver)! reference pressures at midpoints
87 : real(r8), intent(in) :: vkx ! Von Karman's constant
88 : integer, intent(in) :: ntop_eddy ! Top level to which eddy vert diff is applied.
89 : integer, intent(in) :: nbot_eddy ! Bottom level to which eddy vert diff is applied.
90 : character(len=16), intent(in) :: eddy_scheme
91 :
92 : !---------------------------Local workspace-----------------------------
93 : integer :: k ! vertical loop index
94 : !-----------------------------------------------------------------------
95 :
96 : ! Basic constants
97 1536 : cpair = cpairx
98 1536 : g = gravx
99 1536 : vk = vkx
100 1536 : ccon = fak*sffrac*vk
101 1536 : ntop_turb = ntop_eddy
102 1536 : nbot_turb = nbot_eddy
103 :
104 : ! Set the square of the mixing lengths.
105 1536 : ml2(ntop_turb) = 0._r8
106 142848 : do k = ntop_turb+1, nbot_turb
107 141312 : ml2(k) = 30.0_r8**2 ! HB scheme: length scale = 30m
108 142848 : if ( eddy_scheme .eq. 'HBR' ) then
109 0 : ml2(k) = 1.0_r8**2 ! HBR scheme: length scale = 1m
110 : end if
111 : end do
112 1536 : ml2(nbot_turb+1) = 0._r8
113 :
114 : ! Limit pbl height to regions below 400 mb
115 : ! npbl = max number of levels (from bottom) in pbl
116 :
117 1536 : npbl = 0
118 144384 : do k=nbot_turb,ntop_turb,-1
119 144384 : if (pref_mid(k) >= pblmaxp) then
120 43008 : npbl = npbl + 1
121 : end if
122 : end do
123 1536 : npbl = max(npbl,1)
124 :
125 1536 : if (masterproc) then
126 2 : write(iulog,*)'INIT_HB_DIFF: PBL height will be limited to bottom ',npbl, &
127 4 : ' model levels. Top is ',pref_mid(pverp-npbl),' pascals'
128 : end if
129 :
130 1536 : end subroutine init_hb_diff
131 :
132 : !===============================================================================
133 :
134 0 : subroutine compute_hb_diff(ncol , &
135 : th ,t ,q ,z ,zi , &
136 : pmid ,u ,v ,taux ,tauy , &
137 : shflx ,qflx ,obklen ,ustar ,pblh , &
138 : kvm ,kvh ,kvq ,cgh ,cgs , &
139 : tpert ,qpert ,cldn ,ocnfrac ,tke , &
140 : ri , &
141 0 : eddy_scheme)
142 : !-----------------------------------------------------------------------
143 : !
144 : ! Purpose:
145 : ! Interface routines for calcualtion and diatnostics of turbulence related
146 : ! coefficients
147 : !
148 : ! Author: B. Stevens (rewrite August 2000)
149 : !
150 : !-----------------------------------------------------------------------
151 :
152 : use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length, &
153 : calc_eddy_flux_coefficient, calc_ideal_gas_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, &
154 : calc_kinematic_buoyancy_flux
155 : use physconst, only: zvir, rair, gravit, karman
156 :
157 : !------------------------------Arguments--------------------------------
158 : !
159 : ! Input arguments
160 : !
161 : integer, intent(in) :: ncol ! number of atmospheric columns
162 :
163 : real(r8), intent(in) :: th(pcols,pver) ! potential temperature [K]
164 : real(r8), intent(in) :: t(pcols,pver) ! temperature (used for density)
165 : real(r8), intent(in) :: q(pcols,pver) ! specific humidity [kg/kg]
166 : real(r8), intent(in) :: z(pcols,pver) ! height above surface [m]
167 : real(r8), intent(in) :: zi(pcols,pverp) ! height above surface [m]
168 : real(r8), intent(in) :: u(pcols,pver) ! zonal velocity
169 : real(r8), intent(in) :: v(pcols,pver) ! meridional velocity
170 : real(r8), intent(in) :: taux(pcols) ! zonal stress [N/m2]
171 : real(r8), intent(in) :: tauy(pcols) ! meridional stress [N/m2]
172 : real(r8), intent(in) :: shflx(pcols) ! sensible heat flux
173 : real(r8), intent(in) :: qflx(pcols) ! water vapor flux
174 : real(r8), intent(in) :: pmid(pcols,pver) ! midpoint pressures
175 : real(r8), intent(in) :: cldn(pcols,pver) ! new cloud fraction
176 : real(r8), intent(in) :: ocnfrac(pcols) ! Land fraction
177 : character(len=16), intent(in) :: eddy_scheme
178 :
179 : !
180 : ! Output arguments
181 : !
182 : real(r8), intent(out) :: kvm(pcols,pverp) ! eddy diffusivity for momentum [m2/s]
183 : real(r8), intent(out) :: kvh(pcols,pverp) ! eddy diffusivity for heat [m2/s]
184 : real(r8), intent(out) :: kvq(pcols,pverp) ! eddy diffusivity for constituents [m2/s]
185 : real(r8), intent(out) :: cgh(pcols,pverp) ! counter-gradient term for heat [J/kg/m]
186 : real(r8), intent(out) :: cgs(pcols,pverp) ! counter-gradient star (cg/flux)
187 : real(r8), intent(out) :: tpert(pcols) ! convective temperature excess
188 : real(r8), intent(out) :: qpert(pcols) ! convective humidity excess
189 : real(r8), intent(out) :: ustar(pcols) ! surface friction velocity [m/s]
190 : real(r8), intent(out) :: obklen(pcols) ! Obukhov length
191 : real(r8), intent(out) :: pblh(pcols) ! boundary-layer height [m]
192 : real(r8), intent(out) :: tke(pcols,pverp) ! turbulent kinetic energy (estimated)
193 : real(r8), intent(out) :: ri(pcols,pver) ! richardson number: n2/s2
194 : !
195 : !---------------------------Local workspace-----------------------------
196 : !
197 : real(r8) :: thv(pcols,pver) ! virtual temperature
198 : real(r8) :: rrho(pcols) ! 1./bottom level density
199 : real(r8) :: wstar(pcols) ! convective velocity scale [m/s]
200 : real(r8) :: kqfs(pcols) ! kinematic surf constituent flux (kg/m2/s)
201 : real(r8) :: khfs(pcols) ! kinimatic surface heat flux
202 : real(r8) :: kbfs(pcols) ! surface buoyancy flux
203 : real(r8) :: kvf(pcols,pverp) ! free atmospheric eddy diffsvty [m2/s]
204 : real(r8) :: s2(pcols,pver) ! shear squared
205 : real(r8) :: n2(pcols,pver) ! brunt vaisaila frequency
206 : real(r8) :: bge(pcols) ! buoyancy gradient enhancment
207 : integer :: ktopbl(pcols) ! index of first midpoint inside pbl
208 : integer :: i,k
209 : !
210 : ! Initialize time dependent variables that do not depend on pbl height
211 : !
212 :
213 : ! virtual temperature
214 0 : thv(:ncol,ntop_turb:) = calc_virtual_temperature(th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), zvir)
215 :
216 : ! Compute ustar, Obukhov length, and kinematic surface fluxes.
217 0 : rrho(:ncol) = calc_ideal_gas_rrho(rair, t(:ncol,pver), pmid(:ncol,pver))
218 0 : ustar(:ncol) = calc_friction_velocity(taux(:ncol),tauy(:ncol), rrho(:ncol))
219 0 : khfs(:ncol) = calc_kinematic_heat_flux(shflx(:ncol), rrho(:ncol), cpair)
220 0 : kqfs(:ncol) = calc_kinematic_water_vapor_flux(qflx(:ncol), rrho(:ncol))
221 0 : kbfs(:ncol) = calc_kinematic_buoyancy_flux(khfs(:ncol), zvir, th(:ncol,pver), kqfs(:ncol))
222 0 : obklen(:ncol) = calc_obukhov_length(thv(:ncol,pver), ustar(:ncol), gravit, karman, kbfs(:ncol))
223 :
224 : ! Calculate s2, n2, and Richardson number.
225 : call trbintd(ncol , &
226 : thv ,z ,u ,v , &
227 0 : s2 ,n2 ,ri )
228 : !
229 : ! Initialize time dependent variables that do depend on pbl height
230 : !
231 : call pblintd(ncol , &
232 : thv ,z ,u ,v , &
233 : ustar ,obklen ,kbfs ,pblh ,wstar , &
234 0 : zi ,cldn ,ocnfrac ,bge )
235 : !
236 : ! Get atmosphere exchange coefficients
237 : !
238 0 : kvf(:ncol,:) = 0.0_r8
239 0 : do k = ntop_turb, nbot_turb-1
240 0 : do i = 1, ncol
241 0 : kvf(i,k+1) = calc_eddy_flux_coefficient(ml2(k), ri(i, k), s2(i, k))
242 : end do
243 : end do
244 :
245 : !
246 : ! Get pbl exchange coefficients
247 : !
248 : call austausch_pbl(ncol , &
249 : z ,kvf ,kqfs ,khfs ,kbfs , &
250 : obklen ,ustar ,wstar ,pblh ,kvm , &
251 : kvh ,cgh ,cgs ,tpert ,qpert , &
252 0 : ktopbl ,tke ,bge ,eddy_scheme)
253 : !
254 :
255 0 : kvq(:ncol,:) = kvh(:ncol,:)
256 0 : end subroutine compute_hb_diff
257 :
258 58824 : subroutine compute_hb_free_atm_diff(ncol, &
259 : th ,t ,q ,z , &
260 : pmid ,u ,v ,taux ,tauy , &
261 : shflx ,qflx ,obklen ,ustar , &
262 : kvm ,kvh ,kvq ,cgh ,cgs , &
263 : ri )
264 : !-----------------------------------------------------------------------
265 : !
266 : ! This is a version of compute_hb_diff that only computes free
267 : ! atmosphere exchange (no PBL computations)
268 : !
269 : ! Author: B. Stevens (rewrite August 2000)
270 : ! Modified by Thomas Toniazzo and Peter H. Lauritzen (June 2023)
271 : !
272 : !-----------------------------------------------------------------------
273 :
274 : use atmos_phys_pbl_utils, only: calc_virtual_temperature, calc_friction_velocity, calc_obukhov_length, &
275 : calc_free_atm_eddy_flux_coefficient, calc_ideal_gas_rrho, calc_kinematic_heat_flux, calc_kinematic_water_vapor_flux, &
276 : calc_kinematic_buoyancy_flux
277 : use physconst, only: zvir, rair, gravit, karman
278 :
279 : !------------------------------Arguments--------------------------------
280 : !
281 : ! Input arguments
282 : !
283 : integer, intent(in) :: ncol ! number of atmospheric columns
284 :
285 : real(r8), intent(in) :: th(pcols,pver) ! potential temperature [K]
286 : real(r8), intent(in) :: t(pcols,pver) ! temperature (used for density)
287 : real(r8), intent(in) :: q(pcols,pver) ! specific humidity [kg/kg]
288 : real(r8), intent(in) :: z(pcols,pver) ! height above surface [m]
289 : real(r8), intent(in) :: u(pcols,pver) ! zonal velocity
290 : real(r8), intent(in) :: v(pcols,pver) ! meridional velocity
291 : real(r8), intent(in) :: taux(pcols) ! zonal stress [N/m2]
292 : real(r8), intent(in) :: tauy(pcols) ! meridional stress [N/m2]
293 : real(r8), intent(in) :: shflx(pcols) ! sensible heat flux
294 : real(r8), intent(in) :: qflx(pcols) ! water vapor flux
295 : real(r8), intent(in) :: pmid(pcols,pver) ! midpoint pressures
296 : !
297 : ! Output arguments
298 : !
299 : real(r8), intent(out) :: kvm(pcols,pverp) ! eddy diffusivity for momentum [m2/s]
300 : real(r8), intent(out) :: kvh(pcols,pverp) ! eddy diffusivity for heat [m2/s]
301 : real(r8), intent(out) :: kvq(pcols,pverp) ! eddy diffusivity for constituents [m2/s]
302 : real(r8), intent(out) :: cgh(pcols,pverp) ! counter-gradient term for heat [J/kg/m]
303 : real(r8), intent(out) :: cgs(pcols,pverp) ! counter-gradient star (cg/flux)
304 : real(r8), intent(out) :: ustar(pcols) ! surface friction velocity [m/s]
305 : real(r8), intent(out) :: obklen(pcols) ! Obukhov length
306 : real(r8), intent(out) :: ri(pcols,pver) ! richardson number: n2/s2
307 : !
308 : !---------------------------Local workspace-----------------------------
309 : !
310 : real(r8) :: thv(pcols,pver) ! virtual potential temperature
311 : real(r8) :: rrho(pcols) ! 1./bottom level density
312 : real(r8) :: kqfs(pcols) ! kinematic surf constituent flux (kg/m2/s)
313 : real(r8) :: khfs(pcols) ! kinimatic surface heat flux
314 : real(r8) :: kbfs(pcols) ! surface buoyancy flux
315 : real(r8) :: kvf(pcols,pverp) ! free atmospheric eddy diffsvty [m2/s]
316 : real(r8) :: s2(pcols,pver) ! shear squared
317 : real(r8) :: n2(pcols,pver) ! brunt vaisaila frequency
318 : integer :: i, k
319 :
320 : ! virtual potential temperature
321 91405656 : thv(:ncol,ntop_turb:) = calc_virtual_temperature(th(:ncol,ntop_turb:),q(:ncol,ntop_turb:), zvir)
322 :
323 : ! Compute ustar, Obukhov length, and kinematic surface fluxes.
324 982224 : rrho(:ncol) = calc_ideal_gas_rrho(rair, t(:ncol,pver), pmid(:ncol,pver))
325 982224 : ustar(:ncol) = calc_friction_velocity(taux(:ncol),tauy(:ncol), rrho(:ncol))
326 982224 : khfs(:ncol) = calc_kinematic_heat_flux(shflx(:ncol), rrho(:ncol), cpair)
327 982224 : kqfs(:ncol) = calc_kinematic_water_vapor_flux(qflx(:ncol), rrho(:ncol))
328 982224 : kbfs(:ncol) = calc_kinematic_buoyancy_flux(khfs(:ncol), zvir, th(:ncol,pver), kqfs(:ncol))
329 982224 : obklen(:ncol) = calc_obukhov_length(thv(:ncol,pver), ustar(:ncol), gravit, karman, kbfs(:ncol))
330 : ! Calculate s2, n2, and Richardson number.
331 : call trbintd(ncol , &
332 : thv ,z ,u ,v , &
333 58824 : s2 ,n2 ,ri )
334 : !
335 : ! Get free atmosphere exchange coefficients
336 : !
337 92387880 : kvf(:ncol,:) = 0.0_r8
338 5470632 : do k = ntop_turb, nbot_turb - 1
339 90423432 : do i = 1, ncol
340 90364608 : kvf(i,k+1) = calc_free_atm_eddy_flux_coefficient(ml2(k), ri(i, k), s2(i, k))
341 : end do
342 : end do
343 :
344 92387880 : kvq(:ncol,:) = kvf(:ncol,:)
345 92387880 : kvm(:ncol,:) = kvf(:ncol,:)
346 92387880 : kvh(:ncol,:) = kvf(:ncol,:)
347 92387880 : cgh(:ncol,:) = 0._r8
348 92387880 : cgs(:ncol,:) = 0._r8
349 :
350 58824 : end subroutine compute_hb_free_atm_diff
351 :
352 :
353 : !
354 : !===============================================================================
355 58824 : subroutine trbintd(ncol , &
356 : thv ,z ,u ,v , &
357 : s2 ,n2 ,ri )
358 :
359 : !-----------------------------------------------------------------------
360 : !
361 : ! Purpose:
362 : ! Time dependent initialization
363 : !
364 : ! Method:
365 : ! Diagnosis of variables that do not depend on mixing assumptions or
366 : ! PBL depth
367 : !
368 : ! Author: B. Stevens (extracted from pbldiff, August, 2000)
369 : !
370 : !-----------------------------------------------------------------------
371 : !------------------------------Arguments--------------------------------
372 : !
373 : ! Input arguments
374 : !
375 : integer, intent(in) :: ncol ! number of atmospheric columns
376 :
377 : real(r8), intent(in) :: thv(pcols,pver) ! virtual temperature
378 : real(r8), intent(in) :: z(pcols,pver) ! height above surface [m]
379 : real(r8), intent(in) :: u(pcols,pver) ! windspeed x-direction [m/s]
380 : real(r8), intent(in) :: v(pcols,pver) ! windspeed y-direction [m/s]
381 :
382 : !
383 : ! Output arguments
384 : !
385 : real(r8), intent(out) :: s2(pcols,pver) ! shear squared
386 : real(r8), intent(out) :: n2(pcols,pver) ! brunt vaisaila frequency
387 : real(r8), intent(out) :: ri(pcols,pver) ! richardson number: n2/s2
388 : !
389 : !---------------------------Local workspace-----------------------------
390 : !
391 : integer :: i ! longitude index
392 : integer :: k ! level index
393 :
394 : real(r8) :: vvk ! velocity magnitude squared
395 : real(r8) :: dvdz2 ! velocity shear squared
396 : real(r8) :: dz ! delta z between midpoints
397 : !
398 : ! Compute shear squared (s2), brunt vaisaila frequency (n2) and related richardson
399 : ! number (ri). Use virtual temperature to compute n2.
400 : !
401 :
402 5470632 : do k=ntop_turb,nbot_turb-1
403 90423432 : do i=1,ncol
404 84952800 : dvdz2 = (u(i,k)-u(i,k+1))**2 + (v(i,k)-v(i,k+1))**2
405 84952800 : dvdz2 = max(dvdz2,1.e-36_r8)
406 84952800 : dz = z(i,k) - z(i,k+1)
407 84952800 : s2(i,k) = dvdz2/(dz**2)
408 84952800 : n2(i,k) = g*2.0_r8*( thv(i,k) - thv(i,k+1))/((thv(i,k) + thv(i,k+1))*dz)
409 90364608 : ri(i,k) = n2(i,k)/s2(i,k)
410 : end do
411 : end do
412 :
413 58824 : return
414 : end subroutine trbintd
415 : !
416 : !===============================================================================
417 176472 : subroutine pblintd(ncol , &
418 : thv ,z ,u ,v , &
419 : ustar ,obklen ,kbfs ,pblh ,wstar , &
420 : zi ,cldn ,ocnfrac ,bge )
421 : !-----------------------------------------------------------------------
422 : !
423 : ! Purpose:
424 : ! Diagnose standard PBL variables
425 : !
426 : ! Method:
427 : ! Diagnosis of PBL depth and related variables. In this case only wstar.
428 : ! The PBL depth follows:
429 : ! Holtslag, A.A.M., and B.A. Boville, 1993:
430 : ! Local versus Nonlocal Boundary-Layer Diffusion in a Global Climate
431 : ! Model. J. Clim., vol. 6., p. 1825--1842.
432 : !
433 : ! Updated by Holtslag and Hack to exclude the surface layer from the
434 : ! definition of the boundary layer Richardson number. Ri is now defined
435 : ! across the outer layer of the pbl (between the top of the surface
436 : ! layer and the pbl top) instead of the full pbl (between the surface and
437 : ! the pbl top). For simiplicity, the surface layer is assumed to be the
438 : ! region below the first model level (otherwise the boundary layer depth
439 : ! determination would require iteration).
440 : !
441 : ! Modified for boundary layer height diagnosis: Bert Holtslag, june 1994
442 : ! >>>>>>>>> (Use ricr = 0.3 in this formulation)
443 : !
444 : ! Author: B. Stevens (extracted from pbldiff, August 2000)
445 : !
446 : !-----------------------------------------------------------------------
447 : !------------------------------Arguments--------------------------------
448 : !
449 : ! Input arguments
450 : !
451 : integer, intent(in) :: ncol ! number of atmospheric columns
452 :
453 : real(r8), intent(in) :: thv(pcols,pver) ! virtual temperature
454 : real(r8), intent(in) :: z(pcols,pver) ! height above surface [m]
455 : real(r8), intent(in) :: u(pcols,pver) ! windspeed x-direction [m/s]
456 : real(r8), intent(in) :: v(pcols,pver) ! windspeed y-direction [m/s]
457 : real(r8), intent(in) :: ustar(pcols) ! surface friction velocity [m/s]
458 : real(r8), intent(in) :: obklen(pcols) ! Obukhov length
459 : real(r8), intent(in) :: kbfs(pcols) ! sfc kinematic buoyancy flux [m^2/s^3]
460 : real(r8), intent(in) :: zi(pcols,pverp) ! height above surface [m]
461 : real(r8), intent(in) :: cldn(pcols,pver) ! new cloud fraction
462 : real(r8), intent(in) :: ocnfrac(pcols) ! Land fraction
463 :
464 : !
465 : ! Output arguments
466 : !
467 : real(r8), intent(out) :: wstar(pcols) ! convective sclae velocity [m/s]
468 : real(r8), intent(out) :: pblh(pcols) ! boundary-layer height [m]
469 : real(r8), intent(out) :: bge(pcols) ! buoyancy gradient enhancment
470 : !
471 : !---------------------------Local parameters----------------------------
472 : !
473 : real(r8), parameter :: tiny = 1.e-36_r8 ! lower bound for wind magnitude
474 : real(r8), parameter :: fac = 100._r8 ! ustar parameter in height diagnosis
475 : !
476 : !---------------------------Local workspace-----------------------------
477 : !
478 : integer :: i ! longitude index
479 : integer :: k ! level index
480 :
481 : real(r8) :: phiminv(pcols) ! inverse phi function for momentum
482 : real(r8) :: phihinv(pcols) ! inverse phi function for heat
483 : real(r8) :: rino(pcols,pver) ! bulk Richardson no. from level to ref lev
484 : real(r8) :: tlv(pcols) ! ref. level pot tmp + tmp excess
485 : real(r8) :: vvk ! velocity magnitude squared
486 :
487 : logical :: unstbl(pcols) ! pts w/unstbl pbl (positive virtual ht flx)
488 : logical :: check(pcols) ! True=>chk if Richardson no.>critcal
489 : logical :: ocncldcheck(pcols) ! True=>if ocean surface and cloud in lowest layer
490 : !
491 : ! Compute Obukhov length virtual temperature flux and various arrays for use later:
492 : !
493 2946672 : do i=1,ncol
494 2770200 : check(i) = .true.
495 2770200 : rino(i,pver) = 0.0_r8
496 2946672 : pblh(i) = z(i,pver)
497 : end do
498 : !
499 : !
500 : ! PBL height calculation: Scan upward until the Richardson number between
501 : ! the first level and the current level exceeds the "critical" value.
502 : !
503 4941216 : do k=pver-1,pver-npbl+1,-1
504 79736616 : do i=1,ncol
505 79560144 : if (check(i)) then
506 18572607 : vvk = (u(i,k) - u(i,pver))**2 + (v(i,k) - v(i,pver))**2 + fac*ustar(i)**2
507 18572607 : vvk = max(vvk,tiny)
508 18572607 : rino(i,k) = g*(thv(i,k) - thv(i,pver))*(z(i,k)-z(i,pver))/(thv(i,pver)*vvk)
509 18572607 : if (rino(i,k) >= ricr) then
510 2770200 : pblh(i) = z(i,k+1) + (ricr - rino(i,k+1))/(rino(i,k) - rino(i,k+1)) * &
511 2770200 : (z(i,k) - z(i,k+1))
512 2770200 : check(i) = .false.
513 : end if
514 : end if
515 : end do
516 : end do
517 : !
518 : ! Estimate an effective surface temperature to account for surface fluctuations
519 : !
520 2946672 : do i=1,ncol
521 2770200 : if (check(i)) pblh(i) = z(i,pverp-npbl)
522 2770200 : unstbl(i) = (kbfs(i) > 0._r8)
523 2770200 : check(i) = (kbfs(i) > 0._r8)
524 2946672 : if (check(i)) then
525 2003613 : phiminv(i) = (1._r8 - binm*pblh(i)/obklen(i))**onet
526 2003613 : rino(i,pver) = 0.0_r8
527 2003613 : tlv(i) = thv(i,pver) + kbfs(i)*fak/( ustar(i)*phiminv(i) )
528 : end if
529 : end do
530 : !
531 : ! Improve pblh estimate for unstable conditions using the convective temperature excess:
532 : !
533 2946672 : do i = 1,ncol
534 2946672 : bge(i) = 1.e-8_r8
535 : end do
536 4941216 : do k=pver-1,pver-npbl+1,-1
537 79736616 : do i=1,ncol
538 79560144 : if (check(i)) then
539 17648377 : vvk = (u(i,k) - u(i,pver))**2 + (v(i,k) - v(i,pver))**2 + fac*ustar(i)**2
540 17648377 : vvk = max(vvk,tiny)
541 17648377 : rino(i,k) = g*(thv(i,k) - tlv(i))*(z(i,k)-z(i,pver))/(thv(i,pver)*vvk)
542 17648377 : if (rino(i,k) >= ricr) then
543 2003613 : pblh(i) = z(i,k+1) + (ricr - rino(i,k+1))/(rino(i,k) - rino(i,k+1))* &
544 2003613 : (z(i,k) - z(i,k+1))
545 2003613 : bge(i) = 2._r8*g/(thv(i,k)+thv(i,k+1))*(thv(i,k)-thv(i,k+1))/(z(i,k)-z(i,k+1))*pblh(i)
546 2003613 : if (bge(i).lt.0._r8) then
547 141 : bge(i) = 1.e-8_r8
548 : endif
549 2003613 : check(i) = .false.
550 : end if
551 : end if
552 : end do
553 : end do
554 : !
555 : ! PBL height must be greater than some minimum mechanical mixing depth
556 : ! Several investigators have proposed minimum mechanical mixing depth
557 : ! relationships as a function of the local friction velocity, u*. We
558 : ! make use of a linear relationship of the form h = c u* where c=700.
559 : ! The scaling arguments that give rise to this relationship most often
560 : ! represent the coefficient c as some constant over the local coriolis
561 : ! parameter. Here we make use of the experimental results of Koracin
562 : ! and Berkowicz (1988) [BLM, Vol 43] for wich they recommend 0.07/f
563 : ! where f was evaluated at 39.5 N and 52 N. Thus we use a typical mid
564 : ! latitude value for f so that c = 0.07/f = 700. Also, do not allow
565 : ! PBL to exceed some maximum (npbl) number of allowable points
566 : !
567 2946672 : do i=1,ncol
568 2770200 : if (check(i)) pblh(i) = z(i,pverp-npbl)
569 2770200 : pblh(i) = max(pblh(i),700.0_r8*ustar(i))
570 2946672 : wstar(i) = (max(0._r8,kbfs(i))*g*pblh(i)/thv(i,pver))**onet
571 : end do
572 : !
573 : ! Final requirement on PBL heightis that it must be greater than the depth
574 : ! of the lowest model level over ocean if there is any cloud diagnosed in
575 : ! the lowest model level. This is to deal with the inadequacies of the
576 : ! current "dry" formulation of the boundary layer, where this test is
577 : ! used to identify circumstances where there is marine stratus in the
578 : ! lowest level, and to provide a weak ventilation of the layer to avoid
579 : ! a pathology in the cloud scheme (locking in low-level stratiform cloud)
580 : ! If over an ocean surface, and any cloud is diagnosed in the
581 : ! lowest level, set pblh to 50 meters higher than top interface of lowest level
582 : !
583 : ! jrm This is being applied everywhere (not just ocean)!
584 2946672 : do i=1,ncol
585 2770200 : ocncldcheck(i) = .false.
586 2770200 : if (cldn(i,pver).ge.0.0_r8) ocncldcheck(i) = .true.
587 2946672 : if (ocncldcheck(i)) pblh(i) = max(pblh(i),zi(i,pver) + 50._r8)
588 : end do
589 : !
590 176472 : return
591 : end subroutine pblintd
592 : !
593 : !===============================================================================
594 0 : subroutine austausch_pbl(ncol , &
595 : z ,kvf ,kqfs ,khfs ,kbfs , &
596 : obklen ,ustar ,wstar ,pblh ,kvm , &
597 : kvh ,cgh ,cgs ,tpert ,qpert , &
598 0 : ktopbl ,tke ,bge ,eddy_scheme)
599 : !-----------------------------------------------------------------------
600 : !
601 : ! Purpose:
602 : ! Atmospheric Boundary Layer Computation
603 : !
604 : ! Method:
605 : ! Nonlocal scheme that determines eddy diffusivities based on a
606 : ! specified boundary layer height and a turbulent velocity scale;
607 : ! also, countergradient effects for heat and moisture, and constituents
608 : ! are included, along with temperature and humidity perturbations which
609 : ! measure the strength of convective thermals in the lower part of the
610 : ! atmospheric boundary layer.
611 : !
612 : ! For more information, see Holtslag, A.A.M., and B.A. Boville, 1993:
613 : ! Local versus Nonlocal Boundary-Layer Diffusion in a Global Climate
614 : ! Model. J. Clim., vol. 6., p. 1825--1842.
615 : !
616 : ! Updated by Holtslag and Hack to exclude the surface layer from the
617 : ! definition of the boundary layer Richardson number. Ri is now defined
618 : ! across the outer layer of the pbl (between the top of the surface
619 : ! layer and the pbl top) instead of the full pbl (between the surface and
620 : ! the pbl top). For simiplicity, the surface layer is assumed to be the
621 : ! region below the first model level (otherwise the boundary layer depth
622 : ! determination would require iteration).
623 : !
624 : ! Author: B. Boville, B. Stevens (rewrite August 2000)
625 : !
626 : !------------------------------Arguments--------------------------------
627 : !
628 : ! Input arguments
629 : !
630 : integer, intent(in) :: ncol ! number of atmospheric columns
631 :
632 : real(r8), intent(in) :: z(pcols,pver) ! height above surface [m]
633 : real(r8), intent(in) :: kvf(pcols,pverp) ! free atmospheric eddy diffsvty [m2/s]
634 : real(r8), intent(in) :: kqfs(pcols) ! kinematic surf cnstituent flux (kg/m2/s)
635 : real(r8), intent(in) :: khfs(pcols) ! kinimatic surface heat flux
636 : real(r8), intent(in) :: kbfs(pcols) ! surface buoyancy flux
637 : real(r8), intent(in) :: pblh(pcols) ! boundary-layer height [m]
638 : real(r8), intent(in) :: obklen(pcols) ! Obukhov length
639 : real(r8), intent(in) :: ustar(pcols) ! surface friction velocity [m/s]
640 : real(r8), intent(in) :: wstar(pcols) ! convective velocity scale [m/s]
641 : real(r8), intent(in) :: bge(pcols) ! buoyancy gradient enhancment
642 : character(len=16), intent(in) :: eddy_scheme
643 :
644 : !
645 : ! Output arguments
646 : !
647 : real(r8), intent(out) :: kvm(pcols,pverp) ! eddy diffusivity for momentum [m2/s]
648 : real(r8), intent(out) :: kvh(pcols,pverp) ! eddy diffusivity for heat [m2/s]
649 : real(r8), intent(out) :: cgh(pcols,pverp) ! counter-gradient term for heat [J/kg/m]
650 : real(r8), intent(out) :: cgs(pcols,pverp) ! counter-gradient star (cg/flux)
651 : real(r8), intent(out) :: tpert(pcols) ! convective temperature excess
652 : real(r8), intent(out) :: qpert(pcols) ! convective humidity excess
653 :
654 : integer, intent(out) :: ktopbl(pcols) ! index of first midpoint inside pbl
655 : real(r8), intent(out) :: tke(pcols,pverp) ! turbulent kinetic energy (estimated)
656 : !
657 : !---------------------------Local workspace-----------------------------
658 : !
659 : integer :: i ! longitude index
660 : integer :: k ! level index
661 :
662 : real(r8) :: phiminv(pcols) ! inverse phi function for momentum
663 : real(r8) :: phihinv(pcols) ! inverse phi function for heat
664 : real(r8) :: wm(pcols) ! turbulent velocity scale for momentum
665 : real(r8) :: zp(pcols) ! current level height + one level up
666 : real(r8) :: fak1(pcols) ! k*ustar*pblh
667 : real(r8) :: fak2(pcols) ! k*wm*pblh
668 : real(r8) :: fak3(pcols) ! fakn*wstar/wm
669 : real(r8) :: pblk(pcols) ! level eddy diffusivity for momentum
670 : real(r8) :: pr(pcols) ! Prandtl number for eddy diffusivities
671 : real(r8) :: zl(pcols) ! zmzp / Obukhov length
672 : real(r8) :: zh(pcols) ! zmzp / pblh
673 : real(r8) :: zzh(pcols) ! (1-(zmzp/pblh))**2
674 : real(r8) :: zmzp ! level height halfway between zm and zp
675 : real(r8) :: term ! intermediate calculation
676 : real(r8) :: kve ! diffusivity at entrainment layer in unstable cases
677 :
678 : logical :: unstbl(pcols) ! pts w/unstbl pbl (positive virtual ht flx)
679 : logical :: pblpt(pcols) ! pts within pbl
680 : !
681 : ! Initialize height independent arrays
682 : !
683 :
684 : !drb initialize variables for runtime error checking
685 0 : kvm = 0._r8
686 0 : kvh = 0._r8
687 0 : kve = 0._r8
688 0 : cgh = 0._r8
689 0 : cgs = 0._r8
690 0 : tpert = 0._r8
691 0 : qpert = 0._r8
692 0 : ktopbl = 0._r8
693 0 : tke = 0._r8
694 :
695 0 : do i=1,ncol
696 0 : unstbl(i) = (kbfs(i) > 0._r8)
697 0 : pblk(i) = 0.0_r8
698 0 : fak1(i) = ustar(i)*pblh(i)*vk
699 0 : if (unstbl(i)) then
700 0 : phiminv(i) = (1._r8 - binm*pblh(i)/obklen(i))**onet
701 0 : phihinv(i) = sqrt(1._r8 - binh*pblh(i)/obklen(i))
702 0 : wm(i) = ustar(i)*phiminv(i)
703 0 : fak2(i) = wm(i)*pblh(i)*vk
704 0 : fak3(i) = fakn*wstar(i)/wm(i)
705 0 : tpert(i) = max(khfs(i)*fak/wm(i),0._r8)
706 0 : qpert(i) = max(kqfs(i)*fak/wm(i),0._r8)
707 : else
708 0 : tpert(i) = max(khfs(i)*fak/ustar(i),0._r8)
709 0 : qpert(i) = max(kqfs(i)*fak/ustar(i),0._r8)
710 : end if
711 : end do
712 : !
713 : ! Initialize output arrays with free atmosphere values
714 : !
715 0 : do k=1,pverp
716 0 : do i=1,ncol
717 0 : kvm(i,k) = kvf(i,k)
718 0 : kvh(i,k) = kvf(i,k)
719 0 : cgh(i,k) = 0.0_r8
720 0 : cgs(i,k) = 0.0_r8
721 : end do
722 : end do
723 : !
724 : ! Main level loop to compute the diffusivities and counter-gradient terms. These terms are
725 : ! only calculated at points determined to be in the interior of the pbl (pblpt(i)==.true.),
726 : ! and then calculations are directed toward regime: stable vs unstable, surface vs outer
727 : ! layer.
728 : !
729 0 : do k=pver,pver-npbl+2,-1
730 0 : do i=1,ncol
731 0 : pblpt(i) = (z(i,k) < pblh(i))
732 0 : if (pblpt(i)) then
733 0 : ktopbl(i) = k
734 0 : zp(i) = z(i,k-1)
735 : if (zkmin == 0.0_r8 .and. zp(i) > pblh(i)) zp(i) = pblh(i)
736 0 : zmzp = 0.5_r8*(z(i,k) + zp(i)) ! we think this is an approximation to the interface height (where KVs are calculated)
737 0 : zh(i) = zmzp/pblh(i)
738 0 : zl(i) = zmzp/obklen(i)
739 0 : zzh(i) = zh(i)*max(0._r8,(1._r8 - zh(i)))**2
740 0 : if (unstbl(i)) then
741 0 : if (zh(i) < sffrac) then
742 0 : term = (1._r8 - betam*zl(i))**onet
743 0 : pblk(i) = fak1(i)*zzh(i)*term
744 0 : pr(i) = term/sqrt(1._r8 - betah*zl(i))
745 : else
746 0 : pblk(i) = fak2(i)*zzh(i)
747 0 : pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak
748 0 : cgs(i,k) = fak3(i)/(pblh(i)*wm(i))
749 0 : cgh(i,k) = khfs(i)*cgs(i,k)*cpair
750 : end if
751 : else
752 0 : if (zl(i) <= 1._r8) then
753 0 : pblk(i) = fak1(i)*zzh(i)/(1._r8 + betas*zl(i))
754 : else
755 0 : pblk(i) = fak1(i)*zzh(i)/(betas + zl(i))
756 : end if
757 0 : pr(i) = 1._r8
758 : end if
759 0 : kvm(i,k) = max(pblk(i),kvf(i,k))
760 0 : kvh(i,k) = max(pblk(i)/pr(i),kvf(i,k))
761 : end if
762 : end do
763 : end do
764 :
765 : !
766 : ! Check whether last allowed midpoint is within pbl
767 : !
768 :
769 0 : if ( eddy_scheme .eq. 'HBR' ) then
770 : ! apply new diffusivity at entrainment zone
771 0 : do i = 1,ncol
772 0 : if (bge(i) > 1.e-7_r8) then
773 0 : k = ktopbl(i)
774 0 : kve = 0.2_r8*(wstar(i)**3+5._r8*ustar(i)**3)/bge(i)
775 0 : kvm(i,k) = kve
776 0 : kvh(i,k) = kve
777 : end if
778 : end do
779 : end if
780 :
781 : ! Crude estimate of tke (tke=0 above boundary layer)
782 0 : do k = max(pverp-npbl,2),pverp
783 0 : do i = 1, ncol
784 0 : if (z(i,k-1) < pblh(i)) then
785 0 : tke(i,k) = ( kvm(i,k) / pblh(i) ) ** 2
786 : endif
787 : end do
788 : end do
789 0 : return
790 : end subroutine austausch_pbl
791 :
792 : end module hb_diff
|
