| blob | 1e029c02bb2229a47756fa2d4a7e6a51fd34a8fe |
1 ! translated from NN f77 to F90 and put into WRF by Mariusz Pagowski
2 ! NOAA/GSD & CIRA/CSU, Feb 2008
3 ! changes to original code:
4 ! 1. code is 1d (in z)
5 ! 2. no advection of TKE, covariances and variances
6 ! 3. Cranck-Nicholson replaced with the implicit scheme
7 ! 4. removed terrain dependent grid since input in WRF in actual
8 ! distances in z[m]
9 ! 5. cosmetic changes to adhere to WRF standard (remove common blocks,
10 ! intent etc)
11 !-------------------------------------------------------------------
12 !Modifications implemented by Joseph Olson NOAA/GSD/AMB - CU/CIRES
13 !(approved by Mikio Nakanishi):
14 ! 1. Addition of BouLac mixing length in the free atmosphere.
15 ! 2. Changed the turbulent mixing length to be integrated from the
16 ! surface to the top of the BL + a transition layer depth.
17 ! 3. Option to use Kitamura/Canuto modification which removes the
18 ! critical Richardson number and negative TKE (default).
19 ! 4. Hybrid PBL height diagnostic, which blends a theta-v-based
20 ! definition in neutral/convective BL and a TKE-based definition
21 ! in stable conditions.
22 ! For changes 1 and 3, see "JOE's mods" below:
23 !-------------------------------------------------------------------
25 MODULE module_bl_mynn
27 USE module_model_constants, only: &
28 &karman, g, p1000mb, &
29 &cp, r_d, rcp, xlv, &
30 &svp1, svp2, svp3, svpt0, ep_1, ep_2
32 !-------------------------------------------------------------------
33 IMPLICIT NONE
34 !-------------------------------------------------------------------
36 ! The parameters below depend on stability functions of module_sf_mynn.
37 REAL, PARAMETER :: cphm_st=5.0, cphm_unst=16.0, &
38 cphh_st=5.0, cphh_unst=16.0
40 REAL, PARAMETER :: xlvcp=xlv/cp, ev=xlv, rd=r_d, rk=cp/rd, &
41 &svp11=svp1*1.e3, p608=ep_1, ep_3=1.-ep_2
43 REAL, PARAMETER :: tref=300.0 ! reference temperature (K)
44 REAL, PARAMETER :: tv0=p608*tref, tv1=(1.+p608)*tref, gtr=g/tref
46 ! Closure constants
47 REAL, PARAMETER :: &
48 &vk = karman, &
49 &pr = 0.74, &
50 &g1 = 0.235, &
51 &b1 = 24.0, &
52 &b2 = 15.0, & ! CKmod NN2009
53 &c2 = 0.729, & ! 0.729, & !0.75, &
54 &c3 = 0.340, & ! 0.340, & !0.352, &
55 &c4 = 0.0, &
56 &c5 = 0.2, &
57 &a1 = b1*( 1.0-3.0*g1 )/6.0, &
58 ! &c1 = g1 -1.0/( 3.0*a1*b1**(1.0/3.0) ), &
59 &c1 = g1 -1.0/( 3.0*a1*2.88449914061481660), &
60 &a2 = a1*( g1-c1 )/( g1*pr ), &
61 &g2 = b2/b1*( 1.0-c3 ) +2.0*a1/b1*( 3.0-2.0*c2 )
63 REAL, PARAMETER :: &
64 &cc2 = 1.0-c2, &
65 &cc3 = 1.0-c3, &
66 &e1c = 3.0*a2*b2*cc3, &
67 &e2c = 9.0*a1*a2*cc2, &
68 &e3c = 9.0*a2*a2*cc2*( 1.0-c5 ), &
69 &e4c = 12.0*a1*a2*cc2, &
70 &e5c = 6.0*a1*a1
72 ! Constants for length scale
73 REAL, PARAMETER :: qmin=0.0, zmax=1.0, cns=2.7, &
74 !NN2009: &alp1=0.23, alp2=1.0, alp3=5.0, alp4=100.0, Sqfac=3.0
75 &alp1=0.23, alp2=0.53, alp3=5.0, alp4=100.0, Sqfac=2.0
77 ! Constants for gravitational settling
78 ! REAL, PARAMETER :: gno=1.e6/(1.e8)**(2./3.), gpw=5./3., qcgmin=1.e-8
79 REAL, PARAMETER :: gno=4.64158883361278196
80 REAL, PARAMETER :: gpw=5./3., qcgmin=1.e-8,qkemin=1.e-12
81 ! REAL, PARAMETER :: pblh_ref=1500.
83 REAL, PARAMETER :: rr2=0.7071068, rrp=0.3989423
85 !JOE's mods
86 !Use Canuto/Kitamura mod (remove Ric and negative TKE) (1:yes, 0:no)
87 !For more info, see Canuto et al. (2008 JAS) and Kitamura (Journal of the
88 !Meteorological Society of Japan, Vol. 88, No. 5, pp. 857-864, 2010).
89 !Note that this change required further modification of other parameters
90 !above (c2, c3, alp2, and Sqfac). If you want to remove this option, set these
91 !parameters back to NN2009 values (see commented out lines next to the
92 !parameters above). This only removes the negative TKE problem
93 !but does not necessarily improve performance - neutral impact.
94 REAL, PARAMETER :: CKmod=1.
96 !Use BouLac mixing length in free atmosphere (1:yes, 0:no)
97 !This helps remove excessively large mixing in unstable layers aloft.
98 REAL, PARAMETER :: BLmod=1.
99 !JOE-end
101 INTEGER :: mynn_level
103 CONTAINS
105 ! **********************************************************************
106 ! * An improved Mellor-Yamada turbulence closure model *
107 ! * *
108 ! * Aug/2005 M. Nakanishi (N.D.A) *
109 ! * Modified: Dec/2005 M. Nakanishi (N.D.A) *
110 ! * naka@nda.ac.jp *
111 ! * *
112 ! * Contents: *
113 ! * 1. mym_initialize (to be called once initially) *
114 ! * gives the closure constants and initializes the turbulent *
115 ! * quantities. *
116 ! * (2) mym_level2 (called in the other subroutines) *
117 ! * calculates the stability functions at Level 2. *
118 ! * (3) mym_length (called in the other subroutines) *
119 ! * calculates the master length scale. *
120 ! * 4. mym_turbulence *
121 ! * calculates the vertical diffusivity coefficients and the *
122 ! * production terms for the turbulent quantities. *
123 ! * 5. mym_predict *
124 ! * predicts the turbulent quantities at the next step. *
125 ! * 6. mym_condensation *
126 ! * determines the liquid water content and the cloud fraction *
127 ! * diagnostically. *
128 ! * *
129 ! * call mym_initialize *
130 ! * | *
131 ! * |<----------------+ *
132 ! * | | *
133 ! * call mym_condensation | *
134 ! * call mym_turbulence | *
135 ! * call mym_predict | *
136 ! * | | *
137 ! * |-----------------+ *
138 ! * | *
139 ! * end *
140 ! * *
141 ! * Variables worthy of special mention: *
142 ! * tref : Reference temperature *
143 ! * thl : Liquid water potential temperature *
144 ! * qw : Total water (water vapor+liquid water) content *
145 ! * ql : Liquid water content *
146 ! * vt, vq : Functions for computing the buoyancy flux *
147 ! * *
148 ! * If the water contents are unnecessary, e.g., in the case of *
149 ! * ocean models, thl is the potential temperature and qw, ql, vt *
150 ! * and vq are all zero. *
151 ! * *
152 ! * Grid arrangement: *
153 ! * k+1 +---------+ *
154 ! * | | i = 1 - nx *
155 ! * (k) | * | j = 1 - ny *
156 ! * | | k = 1 - nz *
157 ! * k +---------+ *
158 ! * i (i) i+1 *
159 ! * *
160 ! * All the predicted variables are defined at the center (*) of *
161 ! * the grid boxes. The diffusivity coefficients are, however, *
162 ! * defined on the walls of the grid boxes. *
163 ! * # Upper boundary values are given at k=nz. *
164 ! * *
165 ! * References: *
166 ! * 1. Nakanishi, M., 2001: *
167 ! * Boundary-Layer Meteor., 99, 349-378. *
168 ! * 2. Nakanishi, M. and H. Niino, 2004: *
169 ! * Boundary-Layer Meteor., 112, 1-31. *
170 ! * 3. Nakanishi, M. and H. Niino, 2006: *
171 ! * Boundary-Layer Meteor., (in press). *
172 ! * 4. Nakanishi, M. and H. Niino, 2009: *
173 ! * Jour. Meteor. Soc. Japan, 87, 895-912. *
174 ! **********************************************************************
175 !
176 ! SUBROUTINE mym_initialize:
177 !
178 ! Input variables:
179 ! iniflag : <>0; turbulent quantities will be initialized
180 ! = 0; turbulent quantities have been already
181 ! given, i.e., they will not be initialized
182 ! mx, my : Maximum numbers of grid boxes
183 ! in the x and y directions, respectively
184 ! nx, ny, nz : Numbers of the actual grid boxes
185 ! in the x, y and z directions, respectively
186 ! tref : Reference temperature (K)
187 ! dz(nz) : Vertical grid spacings (m)
188 ! # dz(nz)=dz(nz-1)
189 ! zw(nz+1) : Heights of the walls of the grid boxes (m)
190 ! # zw(1)=0.0 and zw(k)=zw(k-1)+dz(k-1)
191 ! h(mx,ny) : G^(1/2) in the terrain-following coordinate
192 ! # h=1-zg/zt, where zg is the height of the
193 ! terrain and zt the top of the model domain
194 ! pi0(mx,my,nz) : Exner function at zw*h+zg (J/kg K)
195 ! defined by c_p*( p_basic/1000hPa )^kappa
196 ! This is usually computed by integrating
197 ! d(pi0)/dz = -h*g/tref.
198 ! rmo(mx,ny) : Inverse of the Obukhov length (m^(-1))
199 ! flt, flq(mx,ny) : Turbulent fluxes of sensible and latent heat,
200 ! respectively, e.g., flt=-u_*Theta_* (K m/s)
201 !! flt - liquid water potential temperature surface flux
202 !! flq - total water flux surface flux
203 ! ust(mx,ny) : Friction velocity (m/s)
204 ! pmz(mx,ny) : phi_m-zeta at z1*h+z0, where z1 (=0.5*dz(1))
205 ! is the first grid point above the surafce, z0
206 ! the roughness length and zeta=(z1*h+z0)*rmo
207 ! phh(mx,ny) : phi_h at z1*h+z0
208 ! u, v(mx,my,nz): Components of the horizontal wind (m/s)
209 ! thl(mx,my,nz) : Liquid water potential temperature
210 ! (K)
211 ! qw(mx,my,nz) : Total water content Q_w (kg/kg)
212 !
213 ! Output variables:
214 ! ql(mx,my,nz) : Liquid water content (kg/kg)
215 ! v?(mx,my,nz) : Functions for computing the buoyancy flux
216 ! qke(mx,my,nz) : Twice the turbulent kinetic energy q^2
217 ! (m^2/s^2)
218 ! tsq(mx,my,nz) : Variance of Theta_l (K^2)
219 ! qsq(mx,my,nz) : Variance of Q_w
220 ! cov(mx,my,nz) : Covariance of Theta_l and Q_w (K)
221 ! el(mx,my,nz) : Master length scale L (m)
222 ! defined on the walls of the grid boxes
223 ! bsh : no longer used
224 ! via common : Closure constants
225 !
226 ! Work arrays: see subroutine mym_level2
227 ! pd?(mx,my,nz) : Half of the production terms at Level 2
228 ! defined on the walls of the grid boxes
229 ! qkw(mx,my,nz) : q on the walls of the grid boxes (m/s)
230 !
231 ! # As to dtl, ...gh, see subroutine mym_turbulence.
232 !
233 !-------------------------------------------------------------------
234 SUBROUTINE mym_initialize ( kts,kte,&
235 & dz, zw, &
236 & u, v, thl, qw, &
237 ! & ust, rmo, pmz, phh, flt, flq,&
238 !JOE-BouLac/PBLH mod
239 & zi,theta,&
240 !JOE-end
241 & ust, rmo, &
242 & Qke, Tsq, Qsq, Cov)
243 !
244 !-------------------------------------------------------------------
246 INTEGER, INTENT(IN) :: kts,kte
247 ! REAL, INTENT(IN) :: ust, rmo, pmz, phh, flt, flq
248 REAL, INTENT(IN) :: ust, rmo
249 REAL, DIMENSION(kts:kte), INTENT(in) :: dz
250 REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw
251 REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw
253 REAL, DIMENSION(kts:kte), INTENT(out) :: qke,tsq,qsq,cov
256 REAL, DIMENSION(kts:kte) :: &
257 &ql,el,pdk,pdt,pdq,pdc,dtl,dqw,dtv,&
258 &gm,gh,sm,sh,qkw,vt,vq
259 INTEGER :: k,l,lmax
260 REAL :: phm,vkz,elq,elv,b1l,b2l,pmz=1.,phh=1.,flt=0.,flq=0.,tmpq
261 !JOE-BouLac and PBLH mod
262 REAL :: zi
263 REAL, DIMENSION(kts:kte) :: theta
264 !JOE-end
267 ! ** At first ql, vt and vq are set to zero. **
268 DO k = kts,kte
269 ql(k) = 0.0
270 vt(k) = 0.0
271 vq(k) = 0.0
272 END DO
273 !
274 CALL mym_level2 ( kts,kte,&
275 & dz, &
276 & u, v, thl, qw, &
277 & ql, vt, vq, &
278 & dtl, dqw, dtv, gm, gh, sm, sh )
279 !
280 ! ** Preliminary setting **
282 el (kts) = 0.0
283 qke(kts) = ust**2 * ( b1*pmz )**(2.0/3.0)
284 !
285 phm = phh*b2 / ( b1*pmz )**(1.0/3.0)
286 tsq(kts) = phm*( flt/ust )**2
287 qsq(kts) = phm*( flq/ust )**2
288 cov(kts) = phm*( flt/ust )*( flq/ust )
289 !
290 DO k = kts+1,kte
291 vkz = vk*zw(k)
292 el (k) = vkz/( 1.0 + vkz/100.0 )
293 qke(k) = 0.0
294 !
295 tsq(k) = 0.0
296 qsq(k) = 0.0
297 cov(k) = 0.0
298 END DO
299 !
300 ! ** Initialization with an iterative manner **
301 ! ** lmax is the iteration count. This is arbitrary. **
302 lmax = 5 !!kte +3
303 !
304 DO l = 1,lmax
305 !
306 CALL mym_length ( kts,kte,&
307 & dz, zw, &
308 & rmo, flt, flq, &
309 & vt, vq, &
310 & qke, &
311 & dtv, &
312 & el, &
313 !JOE-added for BouLac/PBHL Test
314 & zi,theta,&
315 !JOE-end
316 & qkw)
317 !
318 DO k = kts+1,kte
319 elq = el(k)*qkw(k)
320 pdk(k) = elq*( sm(k)*gm (k)+&
321 &sh(k)*gh (k) )
322 pdt(k) = elq* sh(k)*dtl(k)**2
323 pdq(k) = elq* sh(k)*dqw(k)**2
324 pdc(k) = elq* sh(k)*dtl(k)*dqw(k)
325 END DO
326 !
327 ! ** Strictly, vkz*h(i,j) -> vk*( 0.5*dz(1)*h(i,j)+z0 ) **
328 vkz = vk*0.5*dz(kts)
329 !
330 elv = 0.5*( el(kts+1)+el(kts) ) / vkz
331 qke(kts) = ust**2 * ( b1*pmz*elv )**(2.0/3.0)
332 !
333 phm = phh*b2 / ( b1*pmz/elv**2 )**(1.0/3.0)
334 tsq(kts) = phm*( flt/ust )**2
335 qsq(kts) = phm*( flq/ust )**2
336 cov(kts) = phm*( flt/ust )*( flq/ust )
337 !
338 DO k = kts+1,kte-1
339 b1l = b1*0.25*( el(k+1)+el(k) )
340 tmpq=MAX(b1l*( pdk(k+1)+pdk(k) ),qkemin)
341 ! PRINT *,'tmpqqqqq',tmpq,pdk(k+1),pdk(k)
342 qke(k) = tmpq**(2.0/3.0)
344 !
345 IF ( qke(k) .LE. 0.0 ) THEN
346 b2l = 0.0
347 ELSE
348 b2l = b2*( b1l/b1 ) / SQRT( qke(k) )
349 END IF
350 !
351 tsq(k) = b2l*( pdt(k+1)+pdt(k) )
352 qsq(k) = b2l*( pdq(k+1)+pdq(k) )
353 cov(k) = b2l*( pdc(k+1)+pdc(k) )
354 END DO
356 !
357 END DO
359 !! qke(kts)=qke(kts+1)
360 !! tsq(kts)=tsq(kts+1)
361 !! qsq(kts)=qsq(kts+1)
362 !! cov(kts)=cov(kts+1)
364 qke(kte)=qke(kte-1)
365 tsq(kte)=tsq(kte-1)
366 qsq(kte)=qsq(kte-1)
367 cov(kte)=cov(kte-1)
369 !
370 ! RETURN
372 END SUBROUTINE mym_initialize
374 !
375 ! ==================================================================
376 ! SUBROUTINE mym_level2:
377 !
378 ! Input variables: see subroutine mym_initialize
379 !
380 ! Output variables:
381 ! dtl(mx,my,nz) : Vertical gradient of Theta_l (K/m)
382 ! dqw(mx,my,nz) : Vertical gradient of Q_w
383 ! dtv(mx,my,nz) : Vertical gradient of Theta_V (K/m)
384 ! gm (mx,my,nz) : G_M divided by L^2/q^2 (s^(-2))
385 ! gh (mx,my,nz) : G_H divided by L^2/q^2 (s^(-2))
386 ! sm (mx,my,nz) : Stability function for momentum, at Level 2
387 ! sh (mx,my,nz) : Stability function for heat, at Level 2
388 !
389 ! These are defined on the walls of the grid boxes.
390 !
391 SUBROUTINE mym_level2 (kts,kte,&
392 & dz, &
393 & u, v, thl, qw, &
394 & ql, vt, vq, &
395 & dtl, dqw, dtv, gm, gh, sm, sh )
396 !
397 !-------------------------------------------------------------------
399 INTEGER, INTENT(IN) :: kts,kte
400 REAL, DIMENSION(kts:kte), INTENT(in) :: dz
401 REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw,ql,vt,vq
403 REAL, DIMENSION(kts:kte), INTENT(out) :: &
404 &dtl,dqw,dtv,gm,gh,sm,sh
406 INTEGER :: k
408 REAL :: rfc,f1,f2,rf1,rf2,smc,shc,&
409 &ri1,ri2,ri3,ri4,duz,dtz,dqz,vtt,vqq,dtq,dzk,afk,abk,ri,rf
411 !JOE-Canuto/Kitamura mod
412 REAL :: a2den
413 !JOE-end
415 ! ev = 2.5e6
416 ! tv0 = 0.61*tref
417 ! tv1 = 1.61*tref
418 ! gtr = 9.81/tref
419 !
420 rfc = g1/( g1+g2 )
421 f1 = b1*( g1-c1 ) +3.0*a2*( 1.0 -c2 )*( 1.0-c5 ) &
422 & +2.0*a1*( 3.0-2.0*c2 )
423 f2 = b1*( g1+g2 ) -3.0*a1*( 1.0 -c2 )
424 rf1 = b1*( g1-c1 )/f1
425 rf2 = b1* g1 /f2
426 smc = a1 /a2* f1/f2
427 shc = 3.0*a2*( g1+g2 )
428 !
429 ri1 = 0.5/smc
430 ri2 = rf1*smc
431 ri3 = 4.0*rf2*smc -2.0*ri2
432 ri4 = ri2**2
433 !
434 DO k = kts+1,kte
435 dzk = 0.5 *( dz(k)+dz(k-1) )
436 afk = dz(k)/( dz(k)+dz(k-1) )
437 abk = 1.0 -afk
438 duz = ( u(k)-u(k-1) )**2 +( v(k)-v(k-1) )**2
439 duz = duz /dzk**2
440 dtz = ( thl(k)-thl(k-1) )/( dzk )
441 dqz = ( qw(k)-qw(k-1) )/( dzk )
442 !
443 vtt = 1.0 +vt(k)*abk +vt(k-1)*afk
444 vqq = tv0 +vq(k)*abk +vq(k-1)*afk
445 dtq = vtt*dtz +vqq*dqz
446 !
447 dtl(k) = dtz
448 dqw(k) = dqz
449 dtv(k) = dtq
450 !? dtv(i,j,k) = dtz +tv0*dqz
451 !? : +( ev/pi0(i,j,k)-tv1 )
452 !? : *( ql(i,j,k)-ql(i,j,k-1) )/( dzk*h(i,j) )
453 !
454 gm (k) = duz
455 gh (k) = -dtq*gtr
456 !
457 ! ** Gradient Richardson number **
458 ri = -gh(k)/MAX( duz, 1.0e-10 )
460 !JOE-Canuto/Kitamura mod
461 IF (CKmod .eq. 1) THEN
462 a2den = 1. + MAX(ri,0.0)
463 ELSE
464 a2den = 1. + 0.0
465 ENDIF
467 rfc = g1/( g1+g2 )
468 f1 = b1*( g1-c1 ) +3.0*(a2/a2den)*( 1.0 -c2 )*( 1.0-c5 ) &
469 & +2.0*a1*( 3.0-2.0*c2 )
470 f2 = b1*( g1+g2 ) -3.0*a1*( 1.0 -c2 )
471 rf1 = b1*( g1-c1 )/f1
472 rf2 = b1* g1 /f2
473 smc = a1 /(a2/a2den)* f1/f2
474 shc = 3.0*(a2/a2den)*( g1+g2 )
476 ri1 = 0.5/smc
477 ri2 = rf1*smc
478 ri3 = 4.0*rf2*smc -2.0*ri2
479 ri4 = ri2**2
480 !JOE-end
482 ! ** Flux Richardson number **
483 rf = MIN( ri1*( ri+ri2-SQRT(ri**2-ri3*ri+ri4) ), rfc )
484 !
485 sh (k) = shc*( rfc-rf )/( 1.0-rf )
486 sm (k) = smc*( rf1-rf )/( rf2-rf ) * sh(k)
487 END DO
488 !
489 RETURN
491 END SUBROUTINE mym_level2
493 ! ==================================================================
494 ! SUBROUTINE mym_length:
495 !
496 ! Input variables: see subroutine mym_initialize
497 !
498 ! Output variables: see subroutine mym_initialize
499 !
500 ! Work arrays:
501 ! elt(mx,ny) : Length scale depending on the PBL depth (m)
502 ! vsc(mx,ny) : Velocity scale q_c (m/s)
503 ! at first, used for computing elt
504 !
505 SUBROUTINE mym_length ( kts,kte,&
506 & dz, zw, &
507 & rmo, flt, flq, &
508 & vt, vq, &
509 & qke, &
510 & dtv, &
511 & el, &
512 !JOE-added for BouLac ML (PBLH)
513 & zi,theta,&
514 !JOE-end
515 & qkw)
517 !-------------------------------------------------------------------
519 INTEGER, INTENT(IN) :: kts,kte
520 REAL, DIMENSION(kts:kte), INTENT(in) :: dz
521 REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw
522 REAL, INTENT(in) :: rmo,flt,flq
523 REAL, DIMENSION(kts:kte), INTENT(IN) :: qke,vt,vq
525 REAL, DIMENSION(kts:kte), INTENT(out) :: qkw, el
526 REAL, DIMENSION(kts:kte), INTENT(in) :: dtv
528 REAL :: elt,vsc
529 !JOE-added for BouLac ML
530 REAL, DIMENSION(kts:kte), INTENT(IN) :: theta
531 REAL, DIMENSION(kts:kte) :: qtke,elBLmin,elBLavg
532 REAL :: wt,zi,zi2,elt0,h1,h2,alp1mod
534 !THE FOLLOWING LIMITS DO NOT DIRECTLY AFFECT THE ACTUAL PBLH.
535 !THEY ONLY IMPOSE LIMITS ON THE CALCULATION OF THE MIXING LENGTH
536 !SCALES SO THAT THE BOULAC MIXING LENGTH (IN FREE ATMOS) DOES
537 !NOT ENCROACH UPON THE BOUNDARY LAYER MIXING LENGTH (els, elb & elt).
538 REAL, PARAMETER :: minzi = 300. !min mixed-layer height
539 REAL, PARAMETER :: maxdz = 750. !max (half) transition layer depth
540 !=0.3*2500 m PBLH, so the transition
541 !layer stops growing for PBLHs > 2.5 km.
542 REAL, PARAMETER :: mindz = 300. !min (half) transition layer depth
543 !Joe-end
545 INTEGER :: i,j,k
546 REAL :: afk,abk,zwk,dzk,qdz,vflx,bv,elb,els,elf
548 ! tv0 = 0.61*tref
549 ! gtr = 9.81/tref
550 !
551 !JOE-added to impose limits on the height integration for lt as well
552 ! as the transition layer depth - limit artificial gradients.
553 zi2=MAX(zi,minzi)
554 h1=MAX(0.3*zi2,mindz)
555 h1=MIN(h1,maxdz) !limit (1/2) transition layer depth
556 h2=h1/2.0 !~1/4 of the transition layer depth
557 !Joe-end
558 !JOE-added for BouLac ML
559 qtke(kts)=MAX(qke(kts)/2.,0.01)
560 !JOE-end
562 DO k = kts+1,kte
563 afk = dz(k)/( dz(k)+dz(k-1) )
564 abk = 1.0 -afk
565 qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*&
566 &afk,1.0e-10))
568 !JOE- BouLac Start
569 qtke(k) = MAX(qke(k)/2.,0.001)
570 !JOE- BouLac End
572 END DO
573 !
574 elt = 1.0e-5
575 vsc = 1.0e-5
576 !
577 ! ** Strictly, zwk*h(i,j) -> ( zwk*h(i,j)+z0 ) **
578 !JOE-Lt mod: only integrate to top of PBL (+ transition/entrainment
579 ! layer), since TKE aloft is not relevant. Make WHILE loop, so it
580 ! exits after looping through the boundary layer.
581 !
582 k = kts+1
583 zwk = zw(k)
584 DO WHILE (zwk .LE. (zi2+h1))
585 dzk = 0.5*( dz(k)+dz(k-1) )
586 qdz = MAX( qkw(k)-qmin, 0.0 )*dzk
587 elt = elt +qdz*zwk
588 vsc = vsc +qdz
589 k = k+1
590 zwk = zw(k)
591 END DO
592 !
593 elt = alp1*elt/vsc
594 vflx = ( vt(kts)+1.0 )*flt +( vq(kts)+tv0 )*flq
595 vsc = ( gtr*elt*MAX( vflx, 0.0 ) )**(1.0/3.0)
596 !
597 ! ** Strictly, el(i,j,1) is not zero. **
598 el(kts) = 0.0
599 !
600 !JOE- BouLac Start
601 IF ( BLmod .GT. 0. ) THEN
602 ! COMPUTE BouLac mixing length
603 CALL boulac_length(kts,kte,zw,dz,qtke,theta,elBLmin,elBLavg)
604 ENDIF
605 !JOE- BouLac END
607 DO k = kts+1,kte
608 zwk = zw(k)
610 !JOE- BouLac Start - add blending to only use elt in the boundary
611 ! layer and use the BouLac mixing length in free atmos
612 ! [defined relative to the PBLH (zi) + transition layer (h1)].
613 IF ( BLmod .GT. 0. ) THEN
614 wt=.5*TANH((zwk - (zi2+h1))/h2) + .5
615 elt0=elt*(1.-wt) + elBLmin(k)*wt
616 ELSE
617 elt0=elt
618 ENDIF
619 !JOE- BouLac END
621 ! ** Length scale limited by the buoyancy effect **
622 IF ( dtv(k) .GT. 0.0 ) THEN
623 bv = SQRT( gtr*dtv(k) )
624 elb = alp2*qkw(k) / bv &
625 & *( 1.0 + alp3/alp2*&
626 &SQRT( vsc/( bv*elt ) ) )
628 elf = alp2 * qkw(k)/bv
629 ELSE
630 elb = 1.0e10
631 elf = elb
632 END IF
633 !
634 !JOE- BouLac Start - only use BL ML in free atmos.
635 IF ( BLmod .GT. 0. ) THEN
636 wt=.5*TANH((zwk - (zi2+h1))/h2) + .5
637 elb=elb*(1.-wt) + elBLmin(k)*wt
638 !TEST: turn off mixing aloft
639 !elb=elb*(1.-wt) + 0.01*wt
640 ENDIF
641 !!JOE- BouLac END
643 ! ** Length scale in the surface layer **
644 IF ( rmo .GT. 0.0 ) THEN
645 els = vk*zwk &
646 & /( 1.0 + cns*MIN( zwk*rmo, zmax ) )
647 ELSE
648 els = vk*zwk &
649 & *( 1.0 - alp4* zwk*rmo )**0.2
650 END IF
651 !
653 !JOE- BouLac Start
654 ! el(k) = MIN(elb/( elb/elt+elb/els+1.0 ),elf)
655 el(k) = MIN(elb/( elb/elt0+elb/els+1.0 ),elf)
656 ! el(k) = elb/( elb/elt+elb/els+1.0 )
657 !JOE- BouLac End
659 END DO
660 !
661 RETURN
663 END SUBROUTINE mym_length
665 !JOE- BouLac Code Start -
666 ! ==================================================================
667 SUBROUTINE boulac_length(kts,kte,zw,dz,qtke,theta,lb1,lb2)
668 !
669 ! NOTE: This subroutine was taken from the BouLac scheme in WRF-ARW
670 ! and modified for integration into the MYNN PBL scheme.
671 ! WHILE loops were added to reduce the computational expense.
672 ! This subroutine computes the length scales up and down
673 ! and then computes the min, average of the up/down
674 ! length scales, and also considers the distance to the
675 ! surface.
676 !
677 ! dlu = the distance a parcel can be lifted upwards give a finite
678 ! amount of TKE.
679 ! dld = the distance a parcel can be displaced downwards given a
680 ! finite amount of TKE.
681 ! lb1 = the minimum of the length up and length down
682 ! lb2 = the average of the length up and length down
683 !-------------------------------------------------------------------
685 INTEGER, INTENT(IN) :: kts,kte
686 REAL, DIMENSION(kts:kte), INTENT(IN) :: qtke,dz,theta
687 REAL, DIMENSION(kts:kte), INTENT(OUT) :: lb1,lb2
688 REAL, DIMENSION(kts:kte+1), INTENT(IN) :: zw
690 !LOCAL VARS
691 INTEGER :: iz, izz, found
692 REAL, DIMENSION(kts:kte) :: dlu,dld
693 !REAL, PARAMETER :: g=9.81
694 REAL :: dzt, zup, beta, zup_inf, bbb, tl, zdo, zdo_sup, zzz
696 !print*,"IN MYNN-BouLac",kts, kte
698 do iz=kts,kte
700 !----------------------------------
701 ! FIND DISTANCE UPWARD
702 !----------------------------------
703 zup=0.
704 dlu(iz)=zw(kte+1)-zw(iz)-dz(iz)/2.
705 zzz=0.
706 zup_inf=0.
707 beta=g/theta(iz) !Buoyancy coefficient
709 if (iz .lt. kte) then !cant integrate upwards from highest level
711 !do izz=iz,kte-1 !integrate upwards to find dlu
712 found = 0
713 izz=iz
714 DO WHILE (found .EQ. 0)
716 if (izz .lt. kte) then
717 dzt=(dz(izz+1)+dz(izz))/2. ! avg layer depth of above and below layer
718 zup=zup-beta*theta(iz)*dzt ! initial PE the parcel has at iz
719 !print*," ",iz,izz,theta(izz)
720 zup=zup+beta*(theta(izz+1)+theta(izz))*dzt/2. ! PE gained by lifting a parcel to izz+1
721 zzz=zzz+dzt ! depth of layer iz to izz+1
722 if (qtke(iz).lt.zup .and. qtke(iz).ge.zup_inf) then
723 bbb=(theta(izz+1)-theta(izz))/dzt
724 if (bbb .ne. 0.) then
725 tl=(-beta*(theta(izz)-theta(iz)) + &
726 & sqrt( max(0.,(beta*(theta(izz)-theta(iz)))**2. + &
727 & 2.*bbb*beta*(qtke(iz)-zup_inf))))/bbb/beta
728 else
729 if (theta(izz) .ne. theta(iz))then
730 tl=(qtke(iz)-zup_inf)/(beta*(theta(izz)-theta(iz)))
731 else
732 tl=0.
733 endif
734 endif
735 dlu(iz)=zzz-dzt+tl
736 found = 1
737 endif
738 zup_inf=zup
739 izz=izz+1
740 ELSE
741 found = 1
742 ENDIF
744 ENDDO
746 endif
748 !----------------------------------
749 ! FIND DISTANCE DOWN
750 !----------------------------------
751 zdo=0.
752 zdo_sup=0.
753 dld(iz)=zw(iz)+dz(iz)/2.
754 zzz=0.
756 if (iz .gt. kts) then !cant integrate downwards from lowest level
758 !do izz=iz,kts+1,-1 !integrate downwards to find dld
759 found = 0
760 izz=iz
761 DO WHILE (found .EQ. 0)
763 if (izz .gt. kts) then
764 dzt=(dz(izz-1)+dz(izz))/2.
765 zdo=zdo+beta*theta(iz)*dzt
766 zdo=zdo-beta*(theta(izz-1)+theta(izz))*dzt/2.
767 zzz=zzz+dzt
768 if (qtke(iz).lt.zdo .and. qtke(iz).ge.zdo_sup) then
769 bbb=(theta(izz)-theta(izz-1))/dzt
770 if (bbb .ne. 0.) then
771 tl=(beta*(theta(izz)-theta(iz))+ &
772 & sqrt( max(0.,(beta*(theta(izz)-theta(iz)))**2. + &
773 & 2.*bbb*beta*(qtke(iz)-zdo_sup))))/bbb/beta
774 else
775 if (theta(izz) .ne. theta(iz)) then
776 tl=(qtke(iz)-zdo_sup)/(beta*(theta(izz)-theta(iz)))
777 else
778 tl=0.
779 endif
780 endif
781 dld(iz)=zzz-dzt+tl
782 found = 1
783 endif
784 zdo_sup=zdo
785 izz=izz-1
786 ELSE
787 found = 1
788 ENDIF
789 ENDDO
791 endif
793 !----------------------------------
794 ! GET MINIMUM (OR AVERAGE)
795 !----------------------------------
796 !The surface layer length scale can exceed z for large z/L,
797 !so keep maximum distance down > z.
798 dld(iz) = min(dld(iz),zw(iz+1))
799 lb1(iz) = min(dlu(iz),dld(iz)) !minimum
800 lb2(iz) = sqrt(dlu(iz)*dld(iz)) !average - biased towards smallest
801 !lb2(iz) = 0.5*(dlu(iz)+dld(iz)) !average
803 if (iz .eq. kte) then
804 lb1(kte) = lb1(kte-1)
805 lb2(kte) = lb2(kte-1)
806 endif
807 !print*,"IN MYNN-BouLac",kts, kte,lb1(iz)
808 !print*,"IN MYNN-BouLac",iz,dld(iz),dlu(iz)
810 ENDDO
812 END SUBROUTINE boulac_length
813 !
814 !JOE-END BOULAC CODE
816 ! ==================================================================
817 ! SUBROUTINE mym_turbulence:
818 !
819 ! Input variables: see subroutine mym_initialize
820 ! levflag : <>3; Level 2.5
821 ! = 3; Level 3
822 !
823 ! # ql, vt, vq, qke, tsq, qsq and cov are changed to input variables.
824 !
825 ! Output variables: see subroutine mym_initialize
826 ! dfm(mx,my,nz) : Diffusivity coefficient for momentum,
827 ! divided by dz (not dz*h(i,j)) (m/s)
828 ! dfh(mx,my,nz) : Diffusivity coefficient for heat,
829 ! divided by dz (not dz*h(i,j)) (m/s)
830 ! dfq(mx,my,nz) : Diffusivity coefficient for q^2,
831 ! divided by dz (not dz*h(i,j)) (m/s)
832 ! tcd(mx,my,nz) : Countergradient diffusion term for Theta_l
833 ! (K/s)
834 ! qcd(mx,my,nz) : Countergradient diffusion term for Q_w
835 ! (kg/kg s)
836 ! pd?(mx,my,nz) : Half of the production terms
837 !
838 ! Only tcd and qcd are defined at the center of the grid boxes
839 !
840 ! # DO NOT forget that tcd and qcd are added on the right-hand side
841 ! of the equations for Theta_l and Q_w, respectively.
842 !
843 ! Work arrays: see subroutine mym_initialize and level2
844 !
845 ! # dtl, dqw, dtv, gm and gh are allowed to share storage units with
846 ! dfm, dfh, dfq, tcd and qcd, respectively, for saving memory.
847 !
848 SUBROUTINE mym_turbulence ( kts,kte,&
849 & levflag, &
850 & dz, zw, &
851 & u, v, thl, ql, qw, &
852 & qke, tsq, qsq, cov, &
853 & vt, vq,&
854 & rmo, flt, flq, &
855 !JOE-BouLac/PBLH test
856 & zi,theta,&
857 !JOE-end
858 & El,&
859 & Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc, &
860 !JOE-TKE_BUDGET
861 & qWT1D,qSHEAR1D,qBUOY1D,qDISS1D)
862 !JOE-end
864 !-------------------------------------------------------------------
865 !
866 INTEGER, INTENT(IN) :: kts,kte
867 INTEGER, INTENT(IN) :: levflag
868 REAL, DIMENSION(kts:kte), INTENT(in) :: dz
869 REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw
870 REAL, INTENT(in) :: rmo,flt,flq
871 REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw,&
872 &ql,vt,vq,qke,tsq,qsq,cov
874 REAL, DIMENSION(kts:kte), INTENT(out) :: dfm,dfh,dfq,&
875 &pdk,pdt,pdq,pdc,tcd,qcd,el
876 !JOE-TKE-BUDGET
877 REAL, DIMENSION(kts:kte), INTENT(out) :: qWT1D,&
878 &qSHEAR1D,qBUOY1D,qDISS1D
879 REAL :: q3sq_old,dlsq1,qWTP_old,qWTP_new
880 REAL :: dudz,dvdz,dTdz,&
881 upwp,vpwp,Tpwp
882 !JOE-end
884 REAL, DIMENSION(kts:kte) :: qkw,dtl,dqw,dtv,gm,gh,sm,sh
886 INTEGER :: k
887 ! REAL :: cc2,cc3,e1c,e2c,e3c,e4c,e5c
888 REAL :: e6c,dzk,afk,abk,vtt,vqq,&
889 &cw25,clow,cupp,gamt,gamq,smd,gamv,elq,elh
891 !JOE-added for BouLac/PBLH test
892 REAL :: zi
893 REAL, DIMENSION(kts:kte), INTENT(in) :: theta
894 !JOE-end
895 !JOE-Canuto/Kitamura mod
896 REAL :: a2den, duz, ri
897 !JOE-end
899 DOUBLE PRECISION q2sq, t2sq, r2sq, c2sq, elsq, gmel, ghel
900 DOUBLE PRECISION q3sq, t3sq, r3sq, c3sq, dlsq, qdiv
901 DOUBLE PRECISION e1, e2, e3, e4, enum, eden, wden
902 !
903 ! tv0 = 0.61*tref
904 ! gtr = 9.81/tref
905 !
906 ! cc2 = 1.0-c2
907 ! cc3 = 1.0-c3
908 ! e1c = 3.0*a2*b2*cc3
909 ! e2c = 9.0*a1*a2*cc2
910 ! e3c = 9.0*a2*a2*cc2*( 1.0-c5 )
911 ! e4c = 12.0*a1*a2*cc2
912 ! e5c = 6.0*a1*a1
913 !
915 CALL mym_level2 (kts,kte,&
916 & dz, &
917 & u, v, thl, qw, &
918 & ql, vt, vq, &
919 & dtl, dqw, dtv, gm, gh, sm, sh )
920 !
921 CALL mym_length (kts,kte, &
922 & dz, zw, &
923 & rmo, flt, flq, &
924 & vt, vq, &
925 & qke, &
926 & dtv, &
927 & el, &
928 !JOE BouLac/PBLH test
929 & zi,theta,&
930 !JOE-end
931 & qkw)
932 !
934 DO k = kts+1,kte
935 dzk = 0.5 *( dz(k)+dz(k-1) )
936 afk = dz(k)/( dz(k)+dz(k-1) )
937 abk = 1.0 -afk
938 elsq = el (k)**2
939 q2sq = b1*elsq*( sm(k)*gm(k)+sh(k)*gh(k) )
940 q3sq = qkw(k)**2
942 !JOE-Canuto/Kitamura mod
943 duz = ( u(k)-u(k-1) )**2 +( v(k)-v(k-1) )**2
944 duz = duz /dzk**2
945 ! ** Gradient Richardson number **
946 ri = -gh(k)/MAX( duz, 1.0e-10 )
947 IF (CKmod .eq. 1) THEN
948 a2den = 1. + MAX(ri,0.0)
949 ELSE
950 a2den = 1. + 0.0
951 ENDIF
952 !JOE-end
953 !
954 ! Modified: Dec/22/2005, from here, (dlsq -> elsq)
955 gmel = gm (k)*elsq
956 ghel = gh (k)*elsq
957 ! Modified: Dec/22/2005, up to here
958 !
959 ! ** Since qkw is set to more than 0.0, q3sq > 0.0. **
960 IF ( q3sq .LT. q2sq ) THEN
961 qdiv = SQRT( q3sq/q2sq )
962 sm(k) = sm(k) * qdiv
963 sh(k) = sh(k) * qdiv
964 !
965 !JOE-Canuto/Kitamura mod
966 ! e1 = q3sq - e1c*ghel * qdiv**2
967 ! e2 = q3sq - e2c*ghel * qdiv**2
968 ! e3 = e1 + e3c*ghel * qdiv**2
969 ! e4 = e1 - e4c*ghel * qdiv**2
970 e1 = q3sq - e1c*ghel/a2den * qdiv**2
971 e2 = q3sq - e2c*ghel/a2den * qdiv**2
972 e3 = e1 + e3c*ghel/(a2den**2) * qdiv**2
973 e4 = e1 - e4c*ghel/a2den * qdiv**2
974 !JOE-end
975 eden = e2*e4 + e3*e5c*gmel * qdiv**2
976 eden = MAX( eden, 1.0d-20 )
977 ELSE
978 !JOE-Canuto/Kitamura mod
979 ! e1 = q3sq - e1c*ghel
980 ! e2 = q3sq - e2c*ghel
981 ! e3 = e1 + e3c*ghel
982 ! e4 = e1 - e4c*ghel
983 e1 = q3sq - e1c*ghel/a2den
984 e2 = q3sq - e2c*ghel/a2den
985 e3 = e1 + e3c*ghel/(a2den**2)
986 e4 = e1 - e4c*ghel/a2den
987 !JOE-end
988 eden = e2*e4 + e3*e5c*gmel
989 eden = MAX( eden, 1.0d-20 )
990 !
991 qdiv = 1.0
992 sm(k) = q3sq*a1*( e3-3.0*c1*e4 )/eden
993 !JOE-Canuto/Kitamura mod
994 ! sh(k) = q3sq*a2*( e2+3.0*c1*e5c*gmel )/eden
995 sh(k) = q3sq*(a2/a2den)*( e2+3.0*c1*e5c*gmel )/eden
996 !JOE-end
997 END IF
998 !
999 ! ** Level 3 : start **
1000 IF ( levflag .EQ. 3 ) THEN
1001 t2sq = qdiv*b2*elsq*sh(k)*dtl(k)**2
1002 r2sq = qdiv*b2*elsq*sh(k)*dqw(k)**2
1003 c2sq = qdiv*b2*elsq*sh(k)*dtl(k)*dqw(k)
1004 t3sq = MAX( tsq(k)*abk+tsq(k-1)*afk, 0.0 )
1005 r3sq = MAX( qsq(k)*abk+qsq(k-1)*afk, 0.0 )
1006 c3sq = cov(k)*abk+cov(k-1)*afk
1007 !
1008 ! Modified: Dec/22/2005, from here
1009 c3sq = SIGN( MIN( ABS(c3sq), SQRT(t3sq*r3sq) ), c3sq )
1010 !
1011 vtt = 1.0 +vt(k)*abk +vt(k-1)*afk
1012 vqq = tv0 +vq(k)*abk +vq(k-1)*afk
1013 t2sq = vtt*t2sq +vqq*c2sq
1014 r2sq = vtt*c2sq +vqq*r2sq
1015 c2sq = MAX( vtt*t2sq+vqq*r2sq, 0.0d0 )
1016 t3sq = vtt*t3sq +vqq*c3sq
1017 r3sq = vtt*c3sq +vqq*r3sq
1018 c3sq = MAX( vtt*t3sq+vqq*r3sq, 0.0d0 )
1019 !
1020 cw25 = e1*( e2 + 3.0*c1*e5c*gmel*qdiv**2 )/( 3.0*eden )
1021 !
1022 ! ** Limitation on q, instead of L/q **
1023 dlsq = elsq
1024 IF ( q3sq/dlsq .LT. -gh(k) ) q3sq = -dlsq*gh(k)
1025 !
1026 ! ** Limitation on c3sq (0.12 =< cw =< 0.76) **
1027 !JOE-Canuto/Kitamura mod
1028 ! e2 = q3sq - e2c*ghel * qdiv**2
1029 ! e3 = q3sq + e3c*ghel * qdiv**2
1030 ! e4 = q3sq - e4c*ghel * qdiv**2
1031 e2 = q3sq - e2c*ghel/a2den * qdiv**2
1032 e3 = q3sq + e3c*ghel/(a2den**2) * qdiv**2
1033 e4 = q3sq - e4c*ghel/a2den * qdiv**2
1034 !JOE-end
1035 eden = e2*e4 + e3 *e5c*gmel * qdiv**2
1036 !
1037 !JOE-Canuto/Kitamura mod
1038 ! wden = cc3*gtr**2 * dlsq**2/elsq * qdiv**2 &
1039 ! & *( e2*e4c - e3c*e5c*gmel * qdiv**2 )
1040 wden = cc3*gtr**2 * dlsq**2/elsq * qdiv**2 &
1041 & *( e2*e4c/a2den - e3c*e5c*gmel/(a2den**2) * qdiv**2 )
1042 !JOE-end
1043 !
1044 IF ( wden .NE. 0.0 ) THEN
1045 clow = q3sq*( 0.12-cw25 )*eden/wden
1046 cupp = q3sq*( 0.76-cw25 )*eden/wden
1047 !
1048 IF ( wden .GT. 0.0 ) THEN
1049 c3sq = MIN( MAX( c3sq, c2sq+clow ), c2sq+cupp )
1050 ELSE
1051 c3sq = MAX( MIN( c3sq, c2sq+clow ), c2sq+cupp )
1052 END IF
1053 END IF
1054 !
1055 e1 = e2 + e5c*gmel * qdiv**2
1056 eden = MAX( eden, 1.0d-20 )
1057 ! Modified: Dec/22/2005, up to here
1058 !
1059 !JOE-Canuto/Kitamura mod
1060 ! e6c = 3.0*a2*cc3*gtr * dlsq/elsq
1061 e6c = 3.0*(a2/a2den)*cc3*gtr * dlsq/elsq
1062 !JOE-end
1063 !
1064 ! ** for Gamma_theta **
1065 !! enum = qdiv*e6c*( t3sq-t2sq )
1066 IF ( t2sq .GE. 0.0 ) THEN
1067 enum = MAX( qdiv*e6c*( t3sq-t2sq ), 0.0d0 )
1068 ELSE
1069 enum = MIN( qdiv*e6c*( t3sq-t2sq ), 0.0d0 )
1070 ENDIF
1072 gamt =-e1 *enum /eden
1073 !
1074 ! ** for Gamma_q **
1075 !! enum = qdiv*e6c*( r3sq-r2sq )
1076 IF ( r2sq .GE. 0.0 ) THEN
1077 enum = MAX( qdiv*e6c*( r3sq-r2sq ), 0.0d0 )
1078 ELSE
1079 enum = MIN( qdiv*e6c*( r3sq-r2sq ), 0.0d0 )
1080 ENDIF
1082 gamq =-e1 *enum /eden
1083 !
1084 ! ** for Sm' and Sh'd(Theta_V)/dz **
1085 !! enum = qdiv*e6c*( c3sq-c2sq )
1086 enum = MAX( qdiv*e6c*( c3sq-c2sq ), 0.0d0)
1088 !JOE-Canuto/Kitamura mod
1089 ! smd = dlsq*enum*gtr/eden * qdiv**2 * (e3c+e4c)*a1/a2
1090 smd = dlsq*enum*gtr/eden * qdiv**2 * (e3c/(a2den**2) + e4c/a2den)*a1/(a2/a2den)
1091 !JOE-end
1092 gamv = e1 *enum*gtr/eden
1093 !
1095 sm(k) = sm(k) +smd
1096 !
1097 ! ** For elh (see below), qdiv at Level 3 is reset to 1.0. **
1098 qdiv = 1.0
1099 ! ** Level 3 : end **
1100 !
1101 ELSE
1102 ! ** At Level 2.5, qdiv is not reset. **
1103 gamt = 0.0
1104 gamq = 0.0
1105 gamv = 0.0
1106 END IF
1107 !
1108 elq = el(k)*qkw(k)
1109 elh = elq*qdiv
1110 !
1111 pdk(k) = elq*( sm(k)*gm (k) &
1112 & +sh(k)*gh (k)+gamv )
1113 pdt(k) = elh*( sh(k)*dtl(k)+gamt )*dtl(k)
1114 pdq(k) = elh*( sh(k)*dqw(k)+gamq )*dqw(k)
1115 pdc(k) = elh*( sh(k)*dtl(k)+gamt )&
1116 &*dqw(k)*0.5 &
1117 &+elh*( sh(k)*dqw(k)+gamq )*dtl(k)*0.5
1118 !
1119 tcd(k) = elq*gamt
1120 qcd(k) = elq*gamq
1121 !
1122 dfm(k) = elq*sm (k) / dzk
1123 dfh(k) = elq*sh (k) / dzk
1124 ! Modified: Dec/22/2005, from here
1125 ! ** In sub.mym_predict, dfq for the TKE and scalar variance **
1126 ! ** are set to 3.0*dfm and 1.0*dfm, respectively. (Sqfac) **
1127 dfq(k) = dfm(k)
1128 ! Modified: Dec/22/2005, up to here
1130 !TKE BUDGET-start
1131 dudz = ( u(k)-u(k-1) )/dzk
1132 dvdz = ( v(k)-v(k-1) )/dzk
1133 dTdz = ( thl(k)-thl(k-1) )/dzk
1135 upwp = -elq*sm(k)*dudz
1136 vpwp = -elq*sm(k)*dvdz
1137 Tpwp = -elq*sh(k)*dTdz
1138 Tpwp = SIGN(MAX(ABS(Tpwp),1.E-6),Tpwp)
1140 IF ( k .EQ. kts+1 ) THEN
1141 qWT1D(kts)=0.
1142 q3sq_old =0.
1143 qWTP_old =0.
1144 !** Limitation on q, instead of L/q **
1145 dlsq1 = MAX(el(kts)**2,1.0)
1146 IF ( q3sq_old/dlsq1 .LT. -gh(k) ) q3sq_old = -dlsq1*gh(k)
1147 ENDIF
1149 !!!Vertical Transport Term
1150 qWTP_new = elq*Sqfac*sm(k)*(q3sq - q3sq_old)/dzk
1151 qWT1D(k) = 0.5*(qWTP_new - qWTP_old)/dzk
1152 qWTP_old = elq*Sqfac*sm(k)*(q3sq - q3sq_old)/dzk
1153 q3sq_old = q3sq
1155 !!!Shear Term
1156 !!!qSHEAR1D(k)=-(upwp*dudz + vpwp*dvdz)
1157 qSHEAR1D(k) = elq*sm(k)*gm(k)
1159 !!!Buoyancy Term
1160 !!!qBUOY1D(k)=g*Tpwp/thl(k)
1161 !qBUOY1D(k)= elq*(sh(k)*gh(k) + gamv)
1162 qBUOY1D(k) = elq*(sh(k)*(-dTdz*g/thl(k)) + gamv)
1164 !!!Dissipation Term
1165 qDISS1D(k) = (q3sq**(3./2.))/(b1*MAX(el(k),1.))
1166 !TKE BUDGET-end
1168 END DO
1169 !
1171 dfm(kts) = 0.0
1172 dfh(kts) = 0.0
1173 dfq(kts) = 0.0
1174 tcd(kts) = 0.0
1175 qcd(kts) = 0.0
1177 tcd(kte) = 0.0
1178 qcd(kte) = 0.0
1180 !
1181 DO k = kts,kte-1
1182 dzk = dz(k)
1183 tcd(k) = ( tcd(k+1)-tcd(k) )/( dzk )
1184 qcd(k) = ( qcd(k+1)-qcd(k) )/( dzk )
1185 END DO
1186 !
1187 !JOE-TKE-BUDGET
1188 qWT1D(kts)=0.
1189 qSHEAR1D(kts)=qSHEAR1D(kts+1)
1190 qBUOY1D(kts)=qBUOY1D(kts+1)
1191 qDISS1D(kts)=qDISS1D(kts+1)
1192 !JOE-end
1194 RETURN
1196 END SUBROUTINE mym_turbulence
1198 ! ==================================================================
1199 ! SUBROUTINE mym_predict:
1200 !
1201 ! Input variables: see subroutine mym_initialize and turbulence
1202 ! qke(mx,my,nz) : qke at (n)th time level
1203 ! tsq, ...cov : ditto
1204 !
1205 ! Output variables:
1206 ! qke(mx,my,nz) : qke at (n+1)th time level
1207 ! tsq, ...cov : ditto
1208 !
1209 ! Work arrays:
1210 ! qkw(mx,my,nz) : q at the center of the grid boxes (m/s)
1211 ! bp (mx,my,nz) : = 1/2*F, see below
1212 ! rp (mx,my,nz) : = P-1/2*F*Q, see below
1213 !
1214 ! # The equation for a turbulent quantity Q can be expressed as
1215 ! dQ/dt + Ah + Av = Dh + Dv + P - F*Q, (1)
1216 ! where A is the advection, D the diffusion, P the production,
1217 ! F*Q the dissipation and h and v denote horizontal and vertical,
1218 ! respectively. If Q is q^2, F is 2q/B_1L.
1219 ! Using the Crank-Nicholson scheme for Av, Dv and F*Q, a finite
1220 ! difference equation is written as
1221 ! Q{n+1} - Q{n} = dt *( Dh{n} - Ah{n} + P{n} )
1222 ! + dt/2*( Dv{n} - Av{n} - F*Q{n} )
1223 ! + dt/2*( Dv{n+1} - Av{n+1} - F*Q{n+1} ), (2)
1224 ! where n denotes the time level.
1225 ! When the advection and diffusion terms are discretized as
1226 ! dt/2*( Dv - Av ) = a(k)Q(k+1) - b(k)Q(k) + c(k)Q(k-1), (3)
1227 ! Eq.(2) can be rewritten as
1228 ! - a(k)Q(k+1) + [ 1 + b(k) + dt/2*F ]Q(k) - c(k)Q(k-1)
1229 ! = Q{n} + dt *( Dh{n} - Ah{n} + P{n} )
1230 ! + dt/2*( Dv{n} - Av{n} - F*Q{n} ), (4)
1231 ! where Q on the left-hand side is at (n+1)th time level.
1232 !
1233 ! In this subroutine, a(k), b(k) and c(k) are obtained from
1234 ! subprogram coefvu and are passed to subprogram tinteg via
1235 ! common. 1/2*F and P-1/2*F*Q are stored in bp and rp,
1236 ! respectively. Subprogram tinteg solves Eq.(4).
1237 !
1238 ! Modify this subroutine according to your numerical integration
1239 ! scheme (program).
1240 !
1241 !-------------------------------------------------------------------
1242 SUBROUTINE mym_predict (kts,kte,&
1243 & levflag, &
1244 & delt,&
1245 & dz, &
1246 & ust, flt, flq, pmz, phh, &
1247 & el, dfq, &
1248 & pdk, pdt, pdq, pdc,&
1249 & qke, tsq, qsq, cov,&
1250 !JOE-TKE-BUDGET
1251 & dqke1D)
1252 !JOE-end
1254 !-------------------------------------------------------------------
1255 INTEGER, INTENT(IN) :: kts,kte
1256 INTEGER, INTENT(IN) :: levflag
1257 REAL, INTENT(IN) :: delt
1258 REAL, DIMENSION(kts:kte), INTENT(IN) :: dz, dfq,el
1259 REAL, DIMENSION(kts:kte), INTENT(INOUT) :: pdk, pdt, pdq, pdc
1260 REAL, INTENT(IN) :: flt, flq, ust, pmz, phh
1261 REAL, DIMENSION(kts:kte), INTENT(INOUT) :: qke,tsq, qsq, cov
1262 !JOE-TKE-BUDGET
1263 REAL, DIMENSION(kts:kte), INTENT(OUT) :: dqke1D
1264 !JOE-end
1266 INTEGER :: k,nz
1267 REAL, DIMENSION(kts:kte) :: qkw, bp, rp, df3q
1268 REAL :: vkz,pdk1,phm,pdt1,pdq1,pdc1,b1l,b2l
1269 REAL, DIMENSION(kts:kte) :: dtz
1270 REAL, DIMENSION(1:kte-kts+1) :: a,b,c,d
1272 nz=kte-kts+1
1274 ! ** Strictly, vkz*h(i,j) -> vk*( 0.5*dz(1)*h(i,j)+z0 ) **
1275 vkz = vk*0.5*dz(kts)
1276 !
1277 ! Modified: Dec/22/2005, from here
1278 ! ** dfq for the TKE is 3.0*dfm. **
1279 ! CALL coefvu ( dfq, 3.0 ) ! make change here
1280 ! Modified: Dec/22/2005, up to here
1281 !
1282 DO k = kts,kte
1283 !! qke(k) = MAX(qke(k), 0.0)
1284 qkw(k) = SQRT( MAX( qke(k), 0.0 ) )
1285 !df3q(k)=3.*dfq(k)
1286 df3q(k)=Sqfac*dfq(k)
1287 dtz(k)=delt/dz(k)
1288 END DO
1289 !
1290 pdk1 = 2.0*ust**3*pmz/( vkz )
1291 phm = 2.0/ust *phh/( vkz )
1292 pdt1 = phm*flt**2
1293 pdq1 = phm*flq**2
1294 pdc1 = phm*flt*flq
1295 !
1296 ! ** pdk(i,j,1)+pdk(i,j,2) corresponds to pdk1. **
1297 pdk(kts) = pdk1 -pdk(kts+1)
1299 !! pdt(kts) = pdt1 -pdt(kts+1)
1300 !! pdq(kts) = pdq1 -pdq(kts+1)
1301 !! pdc(kts) = pdc1 -pdc(kts+1)
1302 pdt(kts) = pdt(kts+1)
1303 pdq(kts) = pdq(kts+1)
1304 pdc(kts) = pdc(kts+1)
1305 !
1306 ! ** Prediction of twice the turbulent kinetic energy **
1307 !! DO k = kts+1,kte-1
1308 DO k = kts,kte-1
1309 b1l = b1*0.5*( el(k+1)+el(k) )
1310 bp(k) = 2.*qkw(k) / b1l
1311 rp(k) = pdk(k+1) + pdk(k)
1312 END DO
1314 !! a(1)=0.
1315 !! b(1)=1.
1316 !! c(1)=-1.
1317 !! d(1)=0.
1319 ! Since df3q(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)*df3q(k+1)+bp(k)*delt.
1320 DO k=kts,kte-1
1321 a(k-kts+1)=-dtz(k)*df3q(k)
1322 b(k-kts+1)=1.+dtz(k)*(df3q(k)+df3q(k+1))+bp(k)*delt
1323 c(k-kts+1)=-dtz(k)*df3q(k+1)
1324 d(k-kts+1)=rp(k)*delt + qke(k)
1325 ENDDO
1327 !! DO k=kts+1,kte-1
1328 !! a(k-kts+1)=-dtz(k)*df3q(k)
1329 !! b(k-kts+1)=1.+dtz(k)*(df3q(k)+df3q(k+1))
1330 !! c(k-kts+1)=-dtz(k)*df3q(k+1)
1331 !! d(k-kts+1)=rp(k)*delt + qke(k) - qke(k)*bp(k)*delt
1332 !! ENDDO
1334 a(nz)=-1. !0.
1335 b(nz)=1.
1336 c(nz)=0.
1337 d(nz)=0.
1339 CALL tridiag(nz,a,b,c,d)
1341 DO k=kts,kte
1342 !JOE-TKE_BUDGET
1343 dqke1D(k)=qke(k)
1344 !JOE-end
1345 qke(k)=d(k-kts+1)
1346 !JOE-TKE_BUDGET
1347 dqke1D(k)=(qke(k)-dqke1D(k))*0.5
1348 !JOE-end
1349 ENDDO
1352 IF ( levflag .EQ. 3 ) THEN
1353 !
1354 ! Modified: Dec/22/2005, from here
1355 ! ** dfq for the scalar variance is 1.0*dfm. **
1356 ! CALL coefvu ( dfq, 1.0 ) make change here
1357 ! Modified: Dec/22/2005, up to here
1358 !
1359 ! ** Prediction of the temperature variance **
1360 !! DO k = kts+1,kte-1
1361 DO k = kts,kte-1
1362 b2l = b2*0.5*( el(k+1)+el(k) )
1363 bp(k) = 2.*qkw(k) / b2l
1364 rp(k) = pdt(k+1) + pdt(k)
1365 END DO
1367 !zero gradient for tsq at bottom and top
1369 !! a(1)=0.
1370 !! b(1)=1.
1371 !! c(1)=-1.
1372 !! d(1)=0.
1374 ! Since dfq(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)*dfq(k+1)+bp(k)*delt.
1375 DO k=kts,kte-1
1376 a(k-kts+1)=-dtz(k)*dfq(k)
1377 b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))+bp(k)*delt
1378 c(k-kts+1)=-dtz(k)*dfq(k+1)
1379 d(k-kts+1)=rp(k)*delt + tsq(k)
1380 ENDDO
1382 !! DO k=kts+1,kte-1
1383 !! a(k-kts+1)=-dtz(k)*dfq(k)
1384 !! b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))
1385 !! c(k-kts+1)=-dtz(k)*dfq(k+1)
1386 !! d(k-kts+1)=rp(k)*delt + tsq(k) - tsq(k)*bp(k)*delt
1387 !! ENDDO
1389 a(nz)=-1. !0.
1390 b(nz)=1.
1391 c(nz)=0.
1392 d(nz)=0.
1394 CALL tridiag(nz,a,b,c,d)
1396 DO k=kts,kte
1397 tsq(k)=d(k-kts+1)
1398 ENDDO
1400 ! ** Prediction of the moisture variance **
1401 !! DO k = kts+1,kte-1
1402 DO k = kts,kte-1
1403 b2l = b2*0.5*( el(k+1)+el(k) )
1404 bp(k) = 2.*qkw(k) / b2l
1405 rp(k) = pdq(k+1) +pdq(k)
1406 END DO
1408 !zero gradient for qsq at bottom and top
1410 !! a(1)=0.
1411 !! b(1)=1.
1412 !! c(1)=-1.
1413 !! d(1)=0.
1415 ! Since dfq(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)*dfq(k+1)+bp(k)*delt.
1416 DO k=kts,kte-1
1417 a(k-kts+1)=-dtz(k)*dfq(k)
1418 b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))+bp(k)*delt
1419 c(k-kts+1)=-dtz(k)*dfq(k+1)
1420 d(k-kts+1)=rp(k)*delt + qsq(k)
1421 ENDDO
1423 !! DO k=kts+1,kte-1
1424 !! a(k-kts+1)=-dtz(k)*dfq(k)
1425 !! b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))
1426 !! c(k-kts+1)=-dtz(k)*dfq(k+1)
1427 !! d(k-kts+1)=rp(k)*delt + qsq(k) -qsq(k)*bp(k)*delt
1428 !! ENDDO
1430 a(nz)=-1. !0.
1431 b(nz)=1.
1432 c(nz)=0.
1433 d(nz)=0.
1435 CALL tridiag(nz,a,b,c,d)
1437 DO k=kts,kte
1438 qsq(k)=d(k-kts+1)
1439 ENDDO
1441 ! ** Prediction of the temperature-moisture covariance **
1442 !! DO k = kts+1,kte-1
1443 DO k = kts,kte-1
1444 b2l = b2*0.5*( el(k+1)+el(k) )
1445 bp(k) = 2.*qkw(k) / b2l
1446 rp(k) = pdc(k+1) + pdc(k)
1447 END DO
1449 !zero gradient for tqcov at bottom and top
1451 !! a(1)=0.
1452 !! b(1)=1.
1453 !! c(1)=-1.
1454 !! d(1)=0.
1456 ! Since dfq(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)*dfq(k+1)+bp(k)*delt.
1457 DO k=kts,kte-1
1458 a(k-kts+1)=-dtz(k)*dfq(k)
1459 b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))+bp(k)*delt
1460 c(k-kts+1)=-dtz(k)*dfq(k+1)
1461 d(k-kts+1)=rp(k)*delt + cov(k)
1462 ENDDO
1464 !! DO k=kts+1,kte-1
1465 !! a(k-kts+1)=-dtz(k)*dfq(k)
1466 !! b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))
1467 !! c(k-kts+1)=-dtz(k)*dfq(k+1)
1468 !! d(k-kts+1)=rp(k)*delt + cov(k) - cov(k)*bp(k)*delt
1469 !! ENDDO
1471 a(nz)=-1. !0.
1472 b(nz)=1.
1473 c(nz)=0.
1474 d(nz)=0.
1476 CALL tridiag(nz,a,b,c,d)
1478 DO k=kts,kte
1479 cov(k)=d(k-kts+1)
1480 ENDDO
1482 ELSE
1483 !! DO k = kts+1,kte-1
1484 DO k = kts,kte-1
1485 IF ( qkw(k) .LE. 0.0 ) THEN
1486 b2l = 0.0
1487 ELSE
1488 b2l = b2*0.25*( el(k+1)+el(k) )/qkw(k)
1489 END IF
1490 !
1491 tsq(k) = b2l*( pdt(k+1)+pdt(k) )
1492 qsq(k) = b2l*( pdq(k+1)+pdq(k) )
1493 cov(k) = b2l*( pdc(k+1)+pdc(k) )
1494 END DO
1496 !! tsq(kts)=tsq(kts+1)
1497 !! qsq(kts)=qsq(kts+1)
1498 !! cov(kts)=cov(kts+1)
1500 tsq(kte)=tsq(kte-1)
1501 qsq(kte)=qsq(kte-1)
1502 cov(kte)=cov(kte-1)
1504 END IF
1506 END SUBROUTINE mym_predict
1508 ! ==================================================================
1509 ! SUBROUTINE mym_condensation:
1510 !
1511 ! Input variables: see subroutine mym_initialize and turbulence
1512 ! pi (mx,my,nz) : Perturbation of the Exner function (J/kg K)
1513 ! defined on the walls of the grid boxes
1514 ! This is usually computed by integrating
1515 ! d(pi)/dz = h*g*tv/tref**2
1516 ! from the upper boundary, where tv is the
1517 ! virtual potential temperature minus tref.
1518 !
1519 ! Output variables: see subroutine mym_initialize
1520 ! cld(mx,my,nz) : Cloud fraction
1521 !
1522 ! Work arrays:
1523 ! qmq(mx,my,nz) : Q_w-Q_{sl}, where Q_{sl} is the saturation
1524 ! specific humidity at T=Tl
1525 ! alp(mx,my,nz) : Functions in the condensation process
1526 ! bet(mx,my,nz) : ditto
1527 ! sgm(mx,my,nz) : Combined standard deviation sigma_s
1528 ! multiplied by 2/alp
1529 !
1530 ! # qmq, alp, bet and sgm are allowed to share storage units with
1531 ! any four of other work arrays for saving memory.
1532 !
1533 ! # Results are sensitive particularly to values of cp and rd.
1534 ! Set these values to those adopted by you.
1535 !
1536 !-------------------------------------------------------------------
1537 SUBROUTINE mym_condensation (kts,kte, &
1538 & dz, &
1539 & thl, qw, &
1540 & p,exner, &
1541 & tsq, qsq, cov, &
1542 & Vt, Vq)
1544 !-------------------------------------------------------------------
1545 INTEGER, INTENT(IN) :: kts,kte
1547 REAL, DIMENSION(kts:kte), INTENT(IN) :: dz
1548 REAL, DIMENSION(kts:kte), INTENT(IN) :: p,exner, thl, qw, &
1549 &tsq, qsq, cov
1551 REAL, DIMENSION(kts:kte), INTENT(OUT) :: vt,vq
1553 REAL, DIMENSION(kts:kte) :: qmq,alp,bet,sgm,ql,cld
1555 DOUBLE PRECISION :: t3sq, r3sq, c3sq
1556 !
1558 REAL :: p2a,t,esl,qsl,dqsl,q1,cld0,eq1,qll,&
1559 &q2p,pt,rac,qt
1560 INTEGER :: i,j,k
1562 REAL :: erf
1564 ! Note: kte needs to be larger than kts, i.e., kte >= kts+1.
1566 DO k = kts,kte-1
1567 p2a = exner(k)
1568 t = thl(k)*p2a
1570 !x if ( ct .gt. 0.0 ) then
1571 ! a = 17.27
1572 ! b = 237.3
1573 !x else
1574 !x a = 21.87
1575 !x b = 265.5
1576 !x end if
1577 !
1578 ! ** 3.8 = 0.622*6.11 (hPa) **
1579 esl=svp11*EXP(svp2*(t-svpt0)/(t-svp3))
1580 qsl=ep_2*esl/(p(k)-ep_3*esl)
1581 ! qsl = 3.8*EXP( a*ct/( b+ct ) ) / ( 1000.0*p2a**rk )
1582 dqsl = qsl*ep_2*ev/( rd*t**2 )
1583 !
1584 qmq(k) = qw(k) -qsl
1586 alp(k) = 1.0/( 1.0+dqsl*xlvcp )
1587 bet(k) = dqsl*p2a
1588 !
1589 t3sq = MAX( tsq(k), 0.0 )
1590 r3sq = MAX( qsq(k), 0.0 )
1591 c3sq = cov(k)
1592 c3sq = SIGN( MIN( ABS(c3sq), SQRT(t3sq*r3sq) ), c3sq )
1593 !
1594 r3sq = r3sq +bet(k)**2*t3sq -2.0*bet(k)*c3sq
1595 sgm(k) = SQRT( MAX( r3sq, 1.0d-10 ) )
1596 END DO
1597 !
1598 DO k = kts,kte-1
1599 q1 = qmq(k) / sgm(k)
1600 cld0 = 0.5*( 1.0+erf( q1*rr2 ) )
1601 ! q1=0.
1602 ! cld0=0.
1604 eq1 = rrp*EXP( -0.5*q1*q1 )
1605 qll = MAX( cld0*q1 + eq1, 0.0 )
1607 cld(k) = cld0
1608 ql (k) = alp(k)*sgm(k)*qll
1609 !
1610 q2p = xlvcp/exner( k )
1611 pt = thl(k) +q2p*ql(k)
1612 qt = 1.0 +p608*qw(k) -(1.+p608)*ql(k)
1613 rac = alp(k)*( cld0-qll*eq1 )*( q2p*qt-(1.+p608)*pt )
1614 !
1615 vt (k) = qt-1.0 -rac*bet(k)
1616 vq (k) = p608*pt-tv0 +rac
1617 END DO
1618 !
1620 cld(kte) = cld(kte-1)
1621 ql(kte) = ql(kte-1)
1622 vt(kte) = vt(kte-1)
1623 vq(kte) = vq(kte-1)
1625 RETURN
1627 END SUBROUTINE mym_condensation
1629 ! ==================================================================
1630 SUBROUTINE mynn_tendencies(kts,kte,&
1631 &levflag,grav_settling,&
1632 &delt,&
1633 &dz,&
1634 &u,v,th,qv,qc,p,exner,&
1635 &thl,sqv,sqc,sqw,&
1636 &ust,flt,flq,wspd,qcg,&
1637 &tsq,qsq,cov,&
1638 &tcd,qcd,&
1639 &dfm,dfh,dfq,&
1640 &Du,Dv,Dth,Dqv,Dqc)
1642 !-------------------------------------------------------------------
1643 INTEGER, INTENT(in) :: kts,kte
1644 INTEGER, INTENT(in) :: grav_settling,levflag
1646 !! grav_settling = 1 for gravitational settling of droplets
1647 !! grav_settling = 0 otherwise
1648 ! thl - liquid water potential temperature
1649 ! qw - total water
1650 ! dfm,dfh,dfq - as above
1651 ! flt - surface flux of thl
1652 ! flq - surface flux of qw
1654 REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,th,qv,qc,p,exner,&
1655 &dfm,dfh,dfq,dz,tsq,qsq,cov,tcd,qcd
1656 REAL, DIMENSION(kts:kte), INTENT(inout) :: thl,sqw,sqv,sqc
1657 REAL, DIMENSION(kts:kte), INTENT(out) :: du,dv,dth,dqv,dqc
1658 REAL, INTENT(IN) :: delt,ust,flt,flq,wspd,qcg
1660 ! REAL, INTENT(IN) :: delt,ust,flt,flq,qcg,&
1661 ! &gradu_top,gradv_top,gradth_top,gradqv_top
1663 !local vars
1665 REAL, DIMENSION(kts:kte) :: dtz,vt,vq
1667 REAL, DIMENSION(1:kte-kts+1) :: a,b,c,d
1669 REAL :: rhs,gfluxm,gfluxp,dztop
1670 INTEGER :: k,kk,nz
1672 nz=kte-kts+1
1674 dztop=.5*(dz(kte)+dz(kte-1))
1676 DO k=kts,kte
1677 dtz(k)=delt/dz(k)
1678 ENDDO
1680 !! u
1682 k=kts
1684 a(1)=0.
1685 b(1)=1.+dtz(k)*(dfm(k+1)+ust**2/wspd)
1686 c(1)=-dtz(k)*dfm(k+1)
1687 d(1)=u(k)
1689 !! a(1)=0.
1690 !! b(1)=1.+dtz(k)*dfm(k+1)
1691 !! c(1)=-dtz(k)*dfm(k+1)
1692 !! d(1)=u(k)*(1.-ust**2/wspd*dtz(k))
1694 DO k=kts+1,kte-1
1695 kk=k-kts+1
1696 a(kk)=-dtz(k)*dfm(k)
1697 b(kk)=1.+dtz(k)*(dfm(k)+dfm(k+1))
1698 c(kk)=-dtz(k)*dfm(k+1)
1699 d(kk)=u(k)
1700 ENDDO
1702 !! no flux at the top
1704 ! a(nz)=-1.
1705 ! b(nz)=1.
1706 ! c(nz)=0.
1707 ! d(nz)=0.
1709 !! specified gradient at the top
1711 ! a(nz)=-1.
1712 ! b(nz)=1.
1713 ! c(nz)=0.
1714 ! d(nz)=gradu_top*dztop
1716 !! prescribed value
1718 a(nz)=0
1719 b(nz)=1.
1720 c(nz)=0.
1721 d(nz)=u(kte)
1723 CALL tridiag(nz,a,b,c,d)
1725 DO k=kts,kte
1726 du(k)=(d(k-kts+1)-u(k))/delt
1727 ENDDO
1729 !! v
1731 k=kts
1733 a(1)=0.
1734 b(1)=1.+dtz(k)*(dfm(k+1)+ust**2/wspd)
1735 c(1)=-dtz(k)*dfm(k+1)
1736 d(1)=v(k)
1738 !! a(1)=0.
1739 !! b(1)=1.+dtz(k)*dfm(k+1)
1740 !! c(1)=-dtz(k)*dfm(k+1)
1741 !! d(1)=v(k)*(1.-ust**2/wspd*dtz(k))
1743 DO k=kts+1,kte-1
1744 kk=k-kts+1
1745 a(kk)=-dtz(k)*dfm(k)
1746 b(kk)=1.+dtz(k)*(dfm(k)+dfm(k+1))
1747 c(kk)=-dtz(k)*dfm(k+1)
1748 d(kk)=v(k)
1749 ENDDO
1751 !! no flux at the top
1753 ! a(nz)=-1.
1754 ! b(nz)=1.
1755 ! c(nz)=0.
1756 ! d(nz)=0.
1759 !! specified gradient at the top
1761 ! a(nz)=-1.
1762 ! b(nz)=1.
1763 ! c(nz)=0.
1764 ! d(nz)=gradv_top*dztop
1766 !! prescribed value
1768 a(nz)=0
1769 b(nz)=1.
1770 c(nz)=0.
1771 d(nz)=v(kte)
1773 CALL tridiag(nz,a,b,c,d)
1775 DO k=kts,kte
1776 dv(k)=(d(k-kts+1)-v(k))/delt
1777 ENDDO
1779 !! thl
1781 k=kts
1783 a(1)=0.
1784 b(1)=1.+dtz(k)*dfh(k+1)
1785 c(1)=-dtz(k)*dfh(k+1)
1787 ! if qcg not used then assume constant flux in the surface layer
1789 IF (qcg < qcgmin) THEN
1790 IF (sqc(k) > qcgmin) THEN
1791 gfluxm=grav_settling*gno*sqc(k)**gpw
1792 ELSE
1793 gfluxm=0.
1794 ENDIF
1795 ELSE
1796 gfluxm=grav_settling*gno*(qcg/(1.+qcg))**gpw
1797 ENDIF
1799 IF (.5*(sqc(k+1)+sqc(k)) > qcgmin) THEN
1800 gfluxp=grav_settling*gno*(.5*(sqc(k+1)+sqc(k)))**gpw
1801 ELSE
1802 gfluxp=0.
1803 ENDIF
1805 rhs=-xlvcp/exner(k)&
1806 &*( &
1807 &(gfluxp - gfluxm)/dz(k)&
1808 & ) + tcd(k)
1810 d(1)=thl(k)+dtz(k)*flt+rhs*delt
1812 DO k=kts+1,kte-1
1813 kk=k-kts+1
1814 a(kk)=-dtz(k)*dfh(k)
1815 b(kk)=1.+dtz(k)*(dfh(k)+dfh(k+1))
1816 c(kk)=-dtz(k)*dfh(k+1)
1818 IF (.5*(sqc(k+1)+sqc(k)) > qcgmin) THEN
1819 gfluxp=grav_settling*gno*(.5*(sqc(k+1)+sqc(k)))**gpw
1820 ELSE
1821 gfluxp=0.
1822 ENDIF
1824 IF (.5*(sqc(k-1)+sqc(k)) > qcgmin) THEN
1825 gfluxm=grav_settling*gno*(.5*(sqc(k-1)+sqc(k)))**gpw
1826 ELSE
1827 gfluxm=0.
1828 ENDIF
1830 rhs=-xlvcp/exner(k)&
1831 &*( &
1832 &(gfluxp - gfluxm)/dz(k)&
1833 & ) + tcd(k)
1834 d(kk)=thl(k)+rhs*delt
1835 ENDDO
1837 !! no flux at the top
1839 ! a(nz)=-1.
1840 ! b(nz)=1.
1841 ! c(nz)=0.
1842 ! d(nz)=0.
1844 !! specified gradient at the top
1846 !assume gradthl_top=gradth_top
1848 ! a(nz)=-1.
1849 ! b(nz)=1.
1850 ! c(nz)=0.
1851 ! d(nz)=gradth_top*dztop
1853 !! prescribed value
1855 a(nz)=0.
1856 b(nz)=1.
1857 c(nz)=0.
1858 d(nz)=thl(kte)
1860 CALL tridiag(nz,a,b,c,d)
1862 DO k=kts,kte
1863 thl(k)=d(k-kts+1)
1864 ENDDO
1866 !! total water
1868 k=kts
1870 a(1)=0.
1871 b(1)=1.+dtz(k)*dfh(k+1)
1872 c(1)=-dtz(k)*dfh(k+1)
1874 IF (qcg < qcgmin) THEN
1875 IF (sqc(k) > qcgmin) THEN
1876 gfluxm=grav_settling*gno*sqc(k)**gpw
1877 ELSE
1878 gfluxm=0.
1879 ENDIF
1880 ELSE
1881 gfluxm=grav_settling*gno*(qcg/(1.+qcg))**gpw
1882 ENDIF
1884 IF (.5*(sqc(k+1)+sqc(k)) > qcgmin) THEN
1885 gfluxp=grav_settling*gno*(.5*(sqc(k+1)+sqc(k)))**gpw
1886 ELSE
1887 gfluxp=0.
1888 ENDIF
1890 rhs=&
1891 &( &
1892 &(gfluxp - gfluxm)/dz(k)&
1893 & ) + qcd(k)
1895 d(1)=sqw(k)+dtz(k)*flq+rhs*delt
1897 DO k=kts+1,kte-1
1898 kk=k-kts+1
1899 a(kk)=-dtz(k)*dfh(k)
1900 b(kk)=1.+dtz(k)*(dfh(k)+dfh(k+1))
1901 c(kk)=-dtz(k)*dfh(k+1)
1903 IF (.5*(sqc(k+1)+sqc(k)) > qcgmin) THEN
1904 gfluxp=grav_settling*gno*(.5*(sqc(k+1)+sqc(k)))**gpw
1905 ELSE
1906 gfluxp=0.
1907 ENDIF
1909 IF (.5*(sqc(k-1)+sqc(k)) > qcgmin) THEN
1910 gfluxm=grav_settling*gno*(.5*(sqc(k-1)+sqc(k)))**gpw
1911 ELSE
1912 gfluxm=0.
1913 ENDIF
1915 rhs=&
1916 &( &
1917 &(gfluxp - gfluxm)/dz(k)&
1918 & ) + qcd(k)
1919 d(kk)=sqw(k) + rhs*delt
1920 ENDDO
1923 !! no flux at the top
1925 ! a(nz)=-1.
1926 ! b(nz)=1.
1927 ! c(nz)=0.
1928 ! d(nz)=0.
1930 !! specified gradient at the top
1931 !assume gradqw_top=gradqv_top
1933 ! a(nz)=-1.
1934 ! b(nz)=1.
1935 ! c(nz)=0.
1936 ! d(nz)=gradqv_top*dztop
1938 !! prescribed value
1940 a(nz)=0.
1941 b(nz)=1.
1942 c(nz)=0.
1943 d(nz)=sqw(kte)
1945 CALL tridiag(nz,a,b,c,d)
1947 !convert to mixing ratios for wrf
1949 DO k=kts,kte
1950 sqw(k)=d(k-kts+1)
1951 sqv(k)=sqw(k)-sqc(k)
1952 Dqv(k)=(sqv(k)/(1.-sqv(k))-qv(k))/delt
1953 ! Dqc(k)=(sqc(k)/(1.-sqc(k))-qc(k))/delt
1954 Dqc(k)=0.
1955 Dth(k)=(thl(k)+xlvcp/exner(k)*sqc(k)-th(k))/delt
1956 ENDDO
1958 END SUBROUTINE mynn_tendencies
1960 ! ==================================================================
1961 SUBROUTINE retrieve_exchange_coeffs(kts,kte,&
1962 &dfm,dfh,dfq,dz,&
1963 &K_m,K_h,K_q)
1965 !-------------------------------------------------------------------
1967 INTEGER , INTENT(in) :: kts,kte
1969 REAL, DIMENSION(KtS:KtE), INTENT(in) :: dz,dfm,dfh,dfq
1971 REAL, DIMENSION(KtS:KtE), INTENT(out) :: &
1972 &K_m, K_h, K_q
1975 INTEGER :: k
1976 REAL :: dzk
1978 K_m(kts)=0.
1979 K_h(kts)=0.
1980 K_q(kts)=0.
1982 DO k=kts+1,kte
1983 dzk = 0.5 *( dz(k)+dz(k-1) )
1984 K_m(k)=dfm(k)*dzk
1985 K_h(k)=dfh(k)*dzk
1986 K_q(k)=dfq(k)*dzk
1987 ENDDO
1989 END SUBROUTINE retrieve_exchange_coeffs
1991 ! ==================================================================
1992 SUBROUTINE tridiag(n,a,b,c,d)
1994 !! to solve system of linear eqs on tridiagonal matrix n times n
1995 !! after Peaceman and Rachford, 1955
1996 !! a,b,c,d - are vectors of order n
1997 !! a,b,c - are coefficients on the LHS
1998 !! d - is initially RHS on the output becomes a solution vector
2000 !-------------------------------------------------------------------
2002 INTEGER, INTENT(in):: n
2003 REAL, DIMENSION(n), INTENT(in) :: a,b
2004 REAL, DIMENSION(n), INTENT(inout) :: c,d
2006 INTEGER :: i
2007 REAL :: p
2008 REAL, DIMENSION(n) :: q
2010 c(n)=0.
2011 q(1)=-c(1)/b(1)
2012 d(1)=d(1)/b(1)
2014 DO i=2,n
2015 p=1./(b(i)+a(i)*q(i-1))
2016 q(i)=-c(i)*p
2017 d(i)=(d(i)-a(i)*d(i-1))*p
2018 ENDDO
2020 DO i=n-1,1,-1
2021 d(i)=d(i)+q(i)*d(i+1)
2022 ENDDO
2024 END SUBROUTINE tridiag
2026 ! ==================================================================
2027 SUBROUTINE mynn_bl_driver(&
2028 &initflag,&
2029 &grav_settling,&
2030 &delt,&
2031 &dz,&
2032 &u,v,th,qv,qc,&
2033 &p,exner,rho,&
2034 &xland,ts,qsfc,qcg,ps,&
2035 &ust,ch,hfx,qfx,rmol,wspd,&
2036 &Qke,Tsq,Qsq,Cov,&
2037 &Du,Dv,Dth,&
2038 &Dqv,Dqc,&
2039 ! &K_m,K_h,K_q&
2040 &K_h,k_m,&
2041 &Pblh&
2042 !JOE-added for extra ouput
2043 &,el_mynn&
2044 !JOE-end
2045 !JOE-TKE-BUDGET
2046 &,dqke,qWT,qSHEAR,qBUOY,qDISS&
2047 !JOE-end
2048 &,IDS,IDE,JDS,JDE,KDS,KDE &
2049 &,IMS,IME,JMS,JME,KMS,KME &
2050 &,ITS,ITE,JTS,JTE,KTS,KTE)
2052 !-------------------------------------------------------------------
2054 INTEGER, INTENT(in) :: initflag
2055 INTEGER, INTENT(in) :: grav_settling
2057 INTEGER,INTENT(IN) :: &
2058 & IDS,IDE,JDS,JDE,KDS,KDE &
2059 &,IMS,IME,JMS,JME,KMS,KME &
2060 &,ITS,ITE,JTS,JTE,KTS,KTE
2063 ! initflag > 0 for TRUE
2064 ! else for FALSE
2065 ! levflag : <>3; Level 2.5
2066 ! = 3; Level 3
2067 ! grav_settling = 1 when gravitational settling accounted for
2068 ! grav_settling = 0 when gravitational settling NOT accounted for
2070 REAL, INTENT(in) :: delt
2071 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(in) :: dz,&
2072 &u,v,th,qv,qc,p,exner,rho
2073 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(in) :: xland,ust,&
2074 &ch,rmol,ts,qsfc,qcg,ps,hfx,qfx, wspd
2075 !
2076 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(inout) :: &
2077 &Qke,Tsq,Qsq,Cov
2078 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(out) :: &
2079 &Du,Dv,Dth,Dqv,Dqc
2081 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(out) :: &
2082 &K_h,K_m
2084 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(inout) :: &
2085 &Pblh
2087 !JOE-added for extra output
2088 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(inout) :: &
2089 &el_mynn
2090 !JOE-end
2091 !JOE-TKE-BUDGET
2092 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(inout) :: &
2093 &qWT,qSHEAR,qBUOY,qDISS,dqke
2094 !JOE-end
2096 !local vars
2097 INTEGER :: ITF,JTF,KTF
2098 INTEGER :: i,j,k
2099 REAL, DIMENSION(KMS:KME) :: thl,sqv,sqc,sqw,&
2100 &El, Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc, Vt, Vq
2102 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME) :: K_q
2104 REAL, DIMENSION(KMS:KME+1) :: zw
2106 REAL :: cpm,sqcg,flt,flq,pmz,phh,exnerg,zet
2108 REAL, DIMENSION(KMS:KME) :: thetav
2110 INTEGER, SAVE :: levflag
2112 JTF=MIN0(JTE,JDE-1)
2113 ITF=MIN0(ITE,IDE-1)
2114 KTF=MIN0(KTE,KDE-1)
2116 levflag=mynn_level
2118 IF (initflag > 0) THEN
2120 DO j=JTS,JTF
2121 DO i=ITS,ITF
2122 DO k=KTS,KTF
2123 sqv(k)=qv(i,k,j)/(1.+qv(i,k,j))
2124 thl(k)=th(i,k,j)
2125 !JOE-for PBLH
2126 thetav(k)=th(i,k,j)*(1.+0.61*qv(i,k,j))
2127 !JOE-end
2128 IF (k==kts) THEN
2129 zw(k)=0.
2130 ELSE
2131 zw(k)=zw(k-1)+dz(i,k-1,j)
2132 ENDIF
2134 k_m(i,k,j)=0.
2135 k_h(i,k,j)=0.
2136 k_q(i,k,j)=0.
2137 !JOE-added for extra output
2138 el_mynn(i,k,j)=0.
2139 !JOE-end
2141 !JOE-TKE-BUDGET
2142 qWT(i,k,j)=0.
2143 qSHEAR(i,k,j)=0.
2144 qBUOY(i,k,j)=0.
2145 qDISS(i,k,j)=0.
2146 dqke(i,k,j)=0.
2147 !JOE-end
2148 ENDDO
2150 zw(ktf+1)=zw(ktf)+dz(i,ktf,j)
2152 !JOE-PBLH begin
2153 CALL GET_PBLH(KTS,KTE,PBLH(i,j),thetav(kts:kte),&
2154 & Qke(i,kts:kte,j),zw(kts:kte+1),dz(i,kts:kte,j),xland(i,j))
2155 !JOE-PBLH end
2157 CALL mym_initialize ( kts,kte,&
2158 &dz(i,kts:kte,j), zw(kts:kte+1), &
2159 &u(i,kts:kte,j), v(i,kts:kte,j), &
2160 &thl(kts:kte), sqv(kts:kte),&
2161 !JOE-BouLac mod
2162 &PBLH(i,j),th(i,kts:kte,j),&
2163 !JOE-end
2164 &ust(i,j), rmol(i,j),&
2165 &Qke(i,kts:kte,j), Tsq(i,kts:kte,j), &
2166 &Qsq(i,kts:kte,j), Cov(i,kts:kte,j))
2168 ENDDO
2169 ENDDO
2171 ENDIF ! end initflag
2173 DO j=JTS,JTF
2174 DO i=ITS,ITF
2175 DO k=KTS,KTF
2176 sqv(k)=qv(i,k,j)/(1.+qv(i,k,j))
2177 sqc(k)=qc(i,k,j)/(1.+qc(i,k,j))
2178 sqw(k)=sqv(k)+sqc(k)
2179 thl(k)=th(i,k,j)-xlvcp/exner(i,k,j)*sqc(k)
2180 !JOE-for PBLH
2181 thetav(k)=th(i,k,j)*(1.+0.61*qv(i,k,j))
2182 !JOE-end
2184 IF (k==kts) THEN
2185 zw(k)=0.
2186 ELSE
2187 zw(k)=zw(k-1)+dz(i,k-1,j)
2188 ENDIF
2189 ENDDO
2191 zw(ktf+1)=zw(ktf)+dz(i,ktf,j)
2193 !JOE-PBLH begin
2194 CALL GET_PBLH(KTS,KTE,PBLH(i,j),thetav(kts:kte),&
2195 & Qke(i,kts:kte,j),zw(kts:kte+1),dz(i,kts:kte,j),xland(i,j))
2196 !JOE-PBLH end
2198 sqcg=qcg(i,j)/(1.+qcg(i,j))
2199 cpm=cp*(1.+0.8*qv(i,kts,j))
2201 ! The exchange coefficient for cloud water is assumed to be the same as
2202 ! that for heat. CH is multiplied by WSPD. See module_sf_mynn.F
2203 exnerg=(ps(i,j)/p1000mb)**rcp
2204 flt = hfx(i,j)/( rho(i,kts,j)*cpm ) &
2205 +xlvcp*ch(i,j)*(sqc(kts)/exner(i,kts,j)-sqcg/exnerg)
2206 flq = qfx(i,j)/ rho(i,kts,j) &
2207 -ch(i,j)*(sqc(kts) -sqcg )
2209 !!!!!
2210 zet = 0.5*dz(i,kts,j)*rmol(i,j)
2211 if ( zet >= 0.0 ) then
2212 pmz = 1.0 + (cphm_st-1.0) * zet
2213 phh = 1.0 + cphh_st * zet
2214 else
2215 pmz = 1.0/ (1.0-cphm_unst*zet)**0.25 - zet
2216 phh = 1.0/SQRT(1.0-cphh_unst*zet)
2217 end if
2218 !!!!!
2220 CALL mym_condensation ( kts,kte,&
2221 &dz(i,kts:kte,j), &
2222 &thl(kts:kte), sqw(kts:kte), &
2223 &p(i,kts:kte,j),exner(i,kts:kte,j), &
2224 &tsq(i,kts:kte,j), qsq(i,kts:kte,j), cov(i,kts:kte,j), &
2225 &Vt(kts:kte), Vq(kts:kte))
2227 CALL mym_turbulence ( kts,kte,&
2228 &levflag, &
2229 &dz(i,kts:kte,j), zw(kts:kte+1), &
2230 &u(i,kts:kte,j), v(i,kts:kte,j), thl(kts:kte),&
2231 &sqc(kts:kte), sqw(kts:kte), &
2232 &qke(i,kts:kte,j), tsq(i,kts:kte,j), &
2233 &qsq(i,kts:kte,j), cov(i,kts:kte,j), &
2234 &vt(kts:kte), vq(kts:kte),&
2235 &rmol(i,j), flt, flq, &
2236 !JOE-BouLac mod
2237 &PBLH(i,j),th(i,kts:kte,j),&
2238 !JOE-end
2239 &el_mynn(i,kts:kte,j), &
2240 &Dfm(kts:kte),Dfh(kts:kte),Dfq(kts:kte), &
2241 &Tcd(kts:kte),Qcd(kts:kte),Pdk(kts:kte), &
2242 &Pdt(kts:kte),Pdq(kts:kte),Pdc(kts:kte),&
2243 !JOE-TKE-BUDGET
2244 &qWT(i,kts:kte,j),qSHEAR(i,kts:kte,j),&
2245 &qBUOY(i,kts:kte,j),qDISS(i,kts:kte,j))
2246 !JOE-end
2248 CALL mym_predict (kts,kte,&
2249 &levflag, &
2250 &delt,&
2251 &dz(i,kts:kte,j), &
2252 &ust(i,j), flt, flq, pmz, phh, &
2253 &el_mynn(i,kts:kte,j), dfq(kts:kte), pdk(kts:kte), &
2254 &pdt(kts:kte), pdq(kts:kte), pdc(kts:kte),&
2255 &Qke(i,kts:kte,j), Tsq(i,kts:kte,j), &
2256 &Qsq(i,kts:kte,j), Cov(i,kts:kte,j), &
2257 !JOE-TKE-BUDGET
2258 &dqke(i,kts:kte,j))
2259 !JOE-end
2261 CALL mynn_tendencies(kts,kte,&
2262 &levflag,grav_settling,&
2263 &delt,&
2264 &dz(i,kts:kte,j),&
2265 &u(i,kts:kte,j),v(i,kts:kte,j),&
2266 &th(i,kts:kte,j),qv(i,kts:kte,j),qc(i,kts:kte,j),&
2267 &p(i,kts:kte,j),exner(i,kts:kte,j),&
2268 &thl(kts:kte),sqv(kts:kte),sqc(kts:kte),sqw(kts:kte),&
2269 &ust(i,j),flt,flq,wspd(i,j),qcg(i,j),&
2270 &tsq(i,kts:kte,j),qsq(i,kts:kte,j),cov(i,kts:kte,j),&
2271 &tcd(kts:kte),qcd(kts:kte),&
2272 &dfm(kts:kte),dfh(kts:kte),dfq(kts:kte),&
2273 &Du(i,kts:kte,j),Dv(i,kts:kte,j),Dth(i,kts:kte,j),&
2274 &Dqv(i,kts:kte,j),Dqc(i,kts:kte,j))
2276 CALL retrieve_exchange_coeffs(kts,kte,&
2277 &dfm(kts:kte),dfh(kts:kte),dfq(kts:kte),dz(i,kts:kte,j),&
2278 &K_m(i,kts:kte,j),K_h(i,kts:kte,j),K_q(i,kts:kte,j))
2280 DO k=KTS,KTF
2281 qWT(i,k,j)= qWT(i,k,j)*delt
2282 qSHEAR(i,k,j)= qSHEAR(i,k,j)*delt
2283 qBUOY(i,k,j)= qBUOY(i,k,j)*delt
2284 qDISS(i,k,j)= qDISS(i,k,j)*delt
2285 ENDDO
2286 ENDDO
2287 ENDDO
2289 END SUBROUTINE mynn_bl_driver
2291 ! ==================================================================
2292 SUBROUTINE mynn_bl_init_driver(&
2293 &Du,Dv,Dth,&
2294 &Dqv,Dqc&
2295 &,RESTART,ALLOWED_TO_READ,LEVEL&
2296 &,IDS,IDE,JDS,JDE,KDS,KDE &
2297 &,IMS,IME,JMS,JME,KMS,KME &
2298 &,ITS,ITE,JTS,JTE,KTS,KTE)
2300 !---------------------------------------------------------------
2301 LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
2302 INTEGER,INTENT(IN) :: LEVEL
2304 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
2305 & IMS,IME,JMS,JME,KMS,KME, &
2306 & ITS,ITE,JTS,JTE,KTS,KTE
2309 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: &
2310 &Du,Dv,Dth,Dqv,Dqc
2312 INTEGER :: I,J,K,ITF,JTF,KTF
2314 JTF=MIN0(JTE,JDE-1)
2315 KTF=MIN0(KTE,KDE-1)
2316 ITF=MIN0(ITE,IDE-1)
2318 IF(.NOT.RESTART)THEN
2319 DO J=JTS,JTF
2320 DO K=KTS,KTF
2321 DO I=ITS,ITF
2322 Du(i,k,j)=0.
2323 Dv(i,k,j)=0.
2324 Dth(i,k,j)=0.
2325 Dqv(i,k,j)=0.
2326 Dqc(i,k,j)=0.
2327 ENDDO
2328 ENDDO
2329 ENDDO
2330 ENDIF
2332 mynn_level=level
2334 END SUBROUTINE mynn_bl_init_driver
2336 ! ==================================================================
2338 SUBROUTINE GET_PBLH(KTS,KTE,zi,thetav1D,qke1D,zw1D,dz1D,landsea)
2340 !---------------------------------------------------------------
2341 ! NOTES ON THE PBLH FORMULATION
2342 !
2343 !The 1.5-theta-increase method defines PBL heights as the level at
2344 !which the potential temperature first exceeds the minimum potential
2345 !temperature within the boundary layer by 1.5 K. When applied to
2346 !observed temperatures, this method has been shown to produce PBL-
2347 !height estimates that are unbiased relative to profiler-based
2348 !estimates (Nielsen-Gammon et al. 2008). However, their study did not
2349 !include LLJs. Banta and Pichugina (2008) show that a TKE-based
2350 !threshold is a good estimate of the PBL height in LLJs. Therefore,
2351 !a hybrid definition is implemented that uses both methods, weighting
2352 !the TKE-method more during stable conditions (PBLH < 400 m).
2353 !A variable tke threshold (TKEeps) is used since no hard-wired
2354 !value could be found to work best in all conditions.
2355 !---------------------------------------------------------------
2357 INTEGER,INTENT(IN) :: KTS,KTE
2358 REAL, INTENT(OUT) :: zi
2359 REAL, INTENT(IN) :: landsea
2360 REAL, DIMENSION(KTS:KTE), INTENT(IN) :: thetav1D, qke1D, dz1D
2361 REAL, DIMENSION(KTS:KTE+1), INTENT(IN) :: zw1D
2362 !LOCAL VARS
2363 REAL :: PBLH_TKE,qtke,qtkem1,wt,maxqke,TKEeps,minthv
2364 REAL :: delt_thv !delta theta-v; dependent on land/sea point
2365 REAL, PARAMETER :: sbl_lim = 200. !Theta-v PBL lower limit of trust (m).
2366 REAL, PARAMETER :: sbl_damp = 400. !Damping range for averaging with TKE-based PBLH (m).
2367 INTEGER :: I,J,K,kthv,ktke
2369 !FIND MAX TKE AND MIN THETAV IN THE LOWEST 500 M
2370 k = kts+1
2371 kthv = 1
2372 ktke = 1
2373 maxqke = 0.
2374 minthv = 9.E9
2375 DO WHILE (zw1D(k) .LE. 500.)
2376 qtke =MAX(Qke1D(k),0.) ! maximum QKE
2377 IF (maxqke < qtke) then
2378 maxqke = qtke
2379 ktke = k
2380 ENDIF
2381 IF (minthv > thetav1D(k)) then
2382 minthv = thetav1D(k)
2383 kthv = k
2384 ENDIF
2385 k = k+1
2386 ENDDO
2387 !TKEeps = maxtke/20. = maxqke/40.
2388 TKEeps = maxqke/40.
2389 TKEeps = MAX(TKEeps,0.025)
2390 TKEeps = MIN(TKEeps,0.25)
2392 !FIND THETAV-BASED PBLH (BEST FOR DAYTIME).
2393 zi=0.
2394 k = kthv+1
2395 IF((landsea-1.5).GE.0)THEN
2396 ! WATER
2397 delt_thv = 0.75
2398 ELSE
2399 ! LAND
2400 delt_thv = 1.5
2401 ENDIF
2403 zi=0.
2404 k = kthv+1
2405 DO WHILE (zi .EQ. 0.)
2406 IF (thetav1D(k) .GE. (minthv + delt_thv))THEN
2407 zi = zw1D(k) - dz1D(k-1)* &
2408 & MIN((thetav1D(k)-(minthv + delt_thv))/MAX(thetav1D(k)-thetav1D(k-1),1E-6),1.0)
2409 ENDIF
2410 k = k+1
2411 IF (k .EQ. kte-1) zi = zw1D(kts+1) !EXIT SAFEGUARD
2412 ENDDO
2413 !print*,"IN GET_PBLH:",thsfc,zi
2415 !FOR STABLE BOUNDARY LAYERS, USE TKE METHOD TO COMPLEMENT THE
2416 !THETAV-BASED DEFINITION (WHEN THE THETA-V BASED PBLH IS BELOW ~0.5 KM).
2417 !THE TANH WEIGHTING FUNCTION WILL MAKE THE TKE-BASED DEFINITION NEGLIGIBLE
2418 !WHEN THE THETA-V-BASED DEFINITION IS ABOVE ~1 KM.
2419 !FIND TKE-BASED PBLH (BEST FOR NOCTURNAL/STABLE CONDITIONS).
2421 PBLH_TKE=0.
2422 k = ktke+1
2423 DO WHILE (PBLH_TKE .EQ. 0.)
2424 !QKE CAN BE NEGATIVE (IF CKmod == 0)... MAKE TKE NON-NEGATIVE.
2425 qtke =MAX(Qke1D(k)/2.,0.) ! maximum TKE
2426 qtkem1=MAX(Qke1D(k-1)/2.,0.)
2427 IF (qtke .LE. TKEeps) THEN
2428 PBLH_TKE = zw1D(k) - dz1D(k-1)* &
2429 & MIN((TKEeps-qtke)/MAX(qtkem1-qtke, 1E-6), 1.0)
2430 !IN CASE OF NEAR ZERO TKE, SET PBLH = LOWEST LEVEL.
2431 PBLH_TKE = MAX(PBLH_TKE,zw1D(kts+1))
2432 !print *,"PBLH_TKE:",i,j,PBLH_TKE, Qke1D(k)/2., zw1D(kts+1)
2433 ENDIF
2434 k = k+1
2435 IF (k .EQ. kte-1) PBLH_TKE = zw1D(kts+1) !EXIT SAFEGUARD
2436 ENDDO
2438 !BLEND THE TWO PBLH TYPES HERE:
2440 wt=.5*TANH((zi - sbl_lim)/sbl_damp) + .5
2441 zi=PBLH_TKE*(1.-wt) + zi*wt
2443 END SUBROUTINE GET_PBLH
2445 ! ==================================================================
2447 END MODULE module_bl_mynn