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