| blob | 3741e5cc93f1e452469fe6ccbeb0e595cbd056f3 |
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)
12 MODULE module_bl_mynn
14 USE module_model_constants, only: &
15 &karman, g, p1000mb, &
16 &cp, r_d, rcp, xlv, &
17 &svp1, svp2, svp3, svpt0, ep_1, ep_2
19 IMPLICIT NONE
21 ! The parameters below depend on stability functions of module_sf_mynn.
22 REAL, PARAMETER :: cphm_st=5.0, cphm_unst=16.0, &
23 cphh_st=5.0, cphh_unst=16.0
25 REAL, PARAMETER :: xlvcp=xlv/cp, ev=xlv, rd=r_d, rk=cp/rd, &
26 &svp11=svp1*1.e3, p608=ep_1, ep_3=1.-ep_2
28 REAL, PARAMETER :: tref=300.0 ! reference temperature (K)
29 REAL, PARAMETER :: tv0=p608*tref, tv1=(1.+p608)*tref, gtr=g/tref
31 ! Closure constants
32 REAL, PARAMETER :: &
33 &vk = karman, &
34 &pr = 0.74, &
35 &g1 = 0.235, &
36 &b1 = 24.0, &
37 &b2 = 15.0, &
38 &c2 = 0.75, &
39 &c3 = 0.352, &
40 &c4 = 0.0, &
41 &c5 = 0.2, &
42 &a1 = b1*( 1.0-3.0*g1 )/6.0, &
43 ! &c1 = g1 -1.0/( 3.0*a1*b1**(1.0/3.0) ), &
44 &c1 = g1 -1.0/( 3.0*a1*2.88449914061481660), &
45 &a2 = a1*( g1-c1 )/( g1*pr ), &
46 &g2 = b2/b1*( 1.0-c3 ) +2.0*a1/b1*( 3.0-2.0*c2 )
48 REAL, PARAMETER :: &
49 &cc2 = 1.0-c2, &
50 &cc3 = 1.0-c3, &
51 &e1c = 3.0*a2*b2*cc3, &
52 &e2c = 9.0*a1*a2*cc2, &
53 &e3c = 9.0*a2*a2*cc2*( 1.0-c5 ), &
54 &e4c = 12.0*a1*a2*cc2, &
55 &e5c = 6.0*a1*a1
57 ! Constants for length scale
58 REAL, PARAMETER :: qmin=0.0, zmax=1.0, cns=2.7, &
59 &alp1=0.23, alp2=1.0, alp3=5.0, alp4=100.0
61 ! Constants for gravitational settling
62 ! REAL, PARAMETER :: gno=1.e6/(1.e8)**(2./3.)
63 REAL, PARAMETER :: gno=4.64158883361278196
64 REAL, PARAMETER :: gpw=5./3., qcgmin=1.e-8
66 REAL, PARAMETER :: rr2=0.7071068, rrp=0.3989423
68 INTEGER :: mynn_level
70 CONTAINS
72 ! **********************************************************************
73 ! * An improved Mellor-Yamada turbulence closure model *
74 ! * *
75 ! * Aug/2005 M. Nakanishi (N.D.A) *
76 ! * Modified: Dec/2005 M. Nakanishi (N.D.A) *
77 ! * naka@nda.ac.jp *
78 ! * *
79 ! * Contents: *
80 ! * 1. mym_initialize (to be called once initially) *
81 ! * gives the closure constants and initializes the turbulent *
82 ! * quantities. *
83 ! * (2) mym_level2 (called in the other subroutines) *
84 ! * calculates the stability functions at Level 2. *
85 ! * (3) mym_length (called in the other subroutines) *
86 ! * calculates the master length scale. *
87 ! * 4. mym_turbulence *
88 ! * calculates the vertical diffusivity coefficients and the *
89 ! * production terms for the turbulent quantities. *
90 ! * 5. mym_predict *
91 ! * predicts the turbulent quantities at the next step. *
92 ! * 6. mym_condensation *
93 ! * determines the liquid water content and the cloud fraction *
94 ! * diagnostically. *
95 ! * *
96 ! * call mym_initialize *
97 ! * | *
98 ! * |<----------------+ *
99 ! * | | *
100 ! * call mym_condensation | *
101 ! * call mym_turbulence | *
102 ! * call mym_predict | *
103 ! * | | *
104 ! * |-----------------+ *
105 ! * | *
106 ! * end *
107 ! * *
108 ! * Variables worthy of special mention: *
109 ! * tref : Reference temperature *
110 ! * thl : Liquid water potential temperature *
111 ! * qw : Total water (water vapor+liquid water) content *
112 ! * ql : Liquid water content *
113 ! * vt, vq : Functions for computing the buoyancy flux *
114 ! * *
115 ! * If the water contents are unnecessary, e.g., in the case of *
116 ! * ocean models, thl is the potential temperature and qw, ql, vt *
117 ! * and vq are all zero. *
118 ! * *
119 ! * Grid arrangement: *
120 ! * k+1 +---------+ *
121 ! * | | i = 1 - nx *
122 ! * (k) | * | j = 1 - ny *
123 ! * | | k = 1 - nz *
124 ! * k +---------+ *
125 ! * i (i) i+1 *
126 ! * *
127 ! * All the predicted variables are defined at the center (*) of *
128 ! * the grid boxes. The diffusivity coefficients are, however, *
129 ! * defined on the walls of the grid boxes. *
130 ! * # Upper boundary values are given at k=nz. *
131 ! * *
132 ! * References: *
133 ! * 1. Nakanishi, M., 2001: *
134 ! * Boundary-Layer Meteor., 99, 349-378. *
135 ! * 2. Nakanishi, M. and H. Niino, 2004: *
136 ! * Boundary-Layer Meteor., 112, 1-31. *
137 ! * 3. Nakanishi, M. and H. Niino, 2006: *
138 ! * Boundary-Layer Meteor., (in press). *
139 ! **********************************************************************
140 !
141 !
142 ! SUBROUTINE mym_initialize:
143 !
144 ! Input variables:
145 ! iniflag : <>0; turbulent quantities will be initialized
146 ! = 0; turbulent quantities have been already
147 ! given, i.e., they will not be initialized
148 ! mx, my : Maximum numbers of grid boxes
149 ! in the x and y directions, respectively
150 ! nx, ny, nz : Numbers of the actual grid boxes
151 ! in the x, y and z directions, respectively
152 ! tref : Reference temperature (K)
153 ! dz(nz) : Vertical grid spacings (m)
154 ! # dz(nz)=dz(nz-1)
155 ! zw(nz+1) : Heights of the walls of the grid boxes (m)
156 ! # zw(1)=0.0 and zw(k)=zw(k-1)+dz(k-1)
157 ! h(mx,ny) : G^(1/2) in the terrain-following coordinate
158 ! # h=1-zg/zt, where zg is the height of the
159 ! terrain and zt the top of the model domain
160 ! pi0(mx,my,nz) : Exner function at zw*h+zg (J/kg K)
161 ! defined by c_p*( p_basic/1000hPa )^kappa
162 ! This is usually computed by integrating
163 ! d(pi0)/dz = -h*g/tref.
164 ! rmo(mx,ny) : Inverse of the Obukhov length (m^(-1))
165 ! flt, flq(mx,ny) : Turbulent fluxes of sensible and latent heat,
166 ! respectively, e.g., flt=-u_*Theta_* (K m/s)
167 !! flt - liquid water potential temperature surface flux
168 !! flq - total water flux surface flux
169 ! ust(mx,ny) : Friction velocity (m/s)
170 ! pmz(mx,ny) : phi_m-zeta at z1*h+z0, where z1 (=0.5*dz(1))
171 ! is the first grid point above the surafce, z0
172 ! the roughness length and zeta=(z1*h+z0)*rmo
173 ! phh(mx,ny) : phi_h at z1*h+z0
174 ! u, v(mx,my,nz): Components of the horizontal wind (m/s)
175 ! thl(mx,my,nz) : Liquid water potential temperature
176 ! (K)
177 ! qw(mx,my,nz) : Total water content Q_w (kg/kg)
178 !
179 ! Output variables:
180 ! ql(mx,my,nz) : Liquid water content (kg/kg)
181 ! v?(mx,my,nz) : Functions for computing the buoyancy flux
182 ! qke(mx,my,nz) : Twice the turbulent kineti! energy q^2
183 ! (m^2/s^2)
184 ! tsq(mx,my,nz) : Variance of Theta_l (K^2)
185 ! qsq(mx,my,nz) : Variance of Q_w
186 ! cov(mx,my,nz) : Covariance of Theta_l and Q_w (K)
187 ! el(mx,my,nz) : Master length scale L (m)
188 ! defined on the walls of the grid boxes
189 ! bsh : no longer used
190 ! via common : Closure constants
191 !
192 ! Work arrays: see subroutine mym_level2
193 ! pd?(mx,my,nz) : Half of the production terms at Level 2
194 ! defined on the walls of the grid boxes
195 ! qkw(mx,my,nz) : q on the walls of the grid boxes (m/s)
196 !
197 ! # As to dtl, ...gh, see subroutine mym_turbulence.
198 !
200 SUBROUTINE mym_initialize ( kts,kte,&
201 & dz, zw, &
202 & u, v, thl, qw, &
203 ! & ust, rmo, pmz, phh, flt, flq,&
204 & ust, rmo, &
205 & Qke, Tsq, Qsq, Cov)
207 !
209 INTEGER, INTENT(IN) :: kts,kte
210 ! REAL, INTENT(IN) :: ust, rmo, pmz, phh, flt, flq
211 REAL, INTENT(IN) :: ust, rmo
212 REAL, DIMENSION(kts:kte), INTENT(in) :: dz
213 REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw
214 REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw
216 REAL, DIMENSION(kts:kte), INTENT(out) :: qke,tsq,qsq,cov
219 REAL, DIMENSION(kts:kte) :: &
220 &ql,el,pdk,pdt,pdq,pdc,dtl,dqw,dtv,&
221 &gm,gh,sm,sh,qkw,vt,vq
222 INTEGER :: k,l,lmax
223 REAL :: phm,vkz,elq,elv,b1l,b2l,pmz=1.,phh=1.,flt=0.,flq=0.
226 ! ** At first ql, vt and vq are set to zero. **
227 DO k = kts,kte
228 ql(k) = 0.0
229 vt(k) = 0.0
230 vq(k) = 0.0
231 END DO
232 !
233 CALL mym_level2 ( kts,kte,&
234 & dz, &
235 & u, v, thl, qw, &
236 & ql, vt, vq, &
237 & dtl, dqw, dtv, gm, gh, sm, sh )
238 !
239 ! ** Preliminary setting **
241 el (kts) = 0.0
242 qke(kts) = ust**2 * ( b1*pmz )**(2.0/3.0)
243 !
244 phm = phh*b2 / ( b1*pmz )**(1.0/3.0)
245 tsq(kts) = phm*( flt/ust )**2
246 qsq(kts) = phm*( flq/ust )**2
247 cov(kts) = phm*( flt/ust )*( flq/ust )
248 !
249 DO k = kts+1,kte
250 vkz = vk*zw(k)
251 el (k) = vkz/( 1.0 + vkz/100.0 )
252 qke(k) = 0.0
253 !
254 tsq(k) = 0.0
255 qsq(k) = 0.0
256 cov(k) = 0.0
257 END DO
258 !
259 ! ** Initialization with an iterative manner **
260 ! ** lmax is the iteration time. This is arbitrary. **
261 lmax = kte +3
262 !
263 DO l = 1,lmax
264 !
265 CALL mym_length ( kts,kte,&
266 & dz, zw, &
267 & rmo, flt, flq, &
268 & vt, vq, &
269 & qke, &
270 & dtv, &
271 & el, &
272 & qkw)
273 !
274 DO k = kts+1,kte
275 elq = el(k)*qkw(k)
276 pdk(k) = elq*( sm(k)*gm (k)+&
277 &sh(k)*gh (k) )
278 pdt(k) = elq* sh(k)*dtl(k)**2
279 pdq(k) = elq* sh(k)*dqw(k)**2
280 pdc(k) = elq* sh(k)*dtl(k)*dqw(k)
281 END DO
282 !
283 ! ** Strictly, vkz*h(i,j) -> vk*( 0.5*dz(1)*h(i,j)+z0 ) **
284 vkz = vk*0.5*dz(kts)
285 !
286 elv = 0.5*( el(kts+1)+el(kts) ) / vkz
287 qke(kts) = ust**2 * ( b1*pmz*elv )**(2.0/3.0)
288 !
289 phm = phh*b2 / ( b1*pmz/elv**2 )**(1.0/3.0)
290 tsq(kts) = phm*( flt/ust )**2
291 qsq(kts) = phm*( flq/ust )**2
292 cov(kts) = phm*( flt/ust )*( flq/ust )
293 !
294 DO k = kts+1,kte-1
295 b1l = b1*0.25*( el(k+1)+el(k) )
296 qke(k) = &
297 &( b1l*( pdk(k+1)+pdk(k) ) )**(2.0/3.0)
298 !
299 IF ( qke(k) .LE. 0.0 ) THEN
300 b2l = 0.0
301 ELSE
302 b2l = b2*( b1l/b1 ) / SQRT( qke(k) )
303 END IF
304 !
305 tsq(k) = b2l*( pdt(k+1)+pdt(k) )
306 qsq(k) = b2l*( pdq(k+1)+pdq(k) )
307 cov(k) = b2l*( pdc(k+1)+pdc(k) )
308 END DO
310 !
311 END DO
313 !! qke(kts)=qke(kts+1)
314 !! tsq(kts)=tsq(kts+1)
315 !! qsq(kts)=qsq(kts+1)
316 !! cov(kts)=cov(kts+1)
318 qke(kte)=qke(kte-1)
319 tsq(kte)=tsq(kte-1)
320 qsq(kte)=qsq(kte-1)
321 cov(kte)=cov(kte-1)
323 !
324 ! RETURN
326 END SUBROUTINE mym_initialize
328 !
329 ! SUBROUTINE mym_level2:
330 !
331 ! Input variables: see subroutine mym_initialize
332 !
333 ! Output variables:
334 ! dtl(mx,my,nz) : Vertical gradient of Theta_l (K/m)
335 ! dqw(mx,my,nz) : Vertical gradient of Q_w
336 ! dtv(mx,my,nz) : Vertical gradient of Theta_V (K/m)
337 ! gm (mx,my,nz) : G_M divided by L^2/q^2 (s^(-2))
338 ! gh (mx,my,nz) : G_H divided by L^2/q^2 (s^(-2))
339 ! sm (mx,my,nz) : Stability function for momentum, at Level 2
340 ! sh (mx,my,nz) : Stability function for heat, at Level 2
341 !
342 ! These are defined on the walls of the grid boxes.
343 !
344 SUBROUTINE mym_level2 (kts,kte,&
345 & dz, &
346 & u, v, thl, qw, &
347 & ql, vt, vq, &
348 & dtl, dqw, dtv, gm, gh, sm, sh )
349 !
351 INTEGER, INTENT(IN) :: kts,kte
352 REAL, DIMENSION(kts:kte), INTENT(in) :: dz
353 REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw,ql,vt,vq
356 REAL, DIMENSION(kts:kte), INTENT(out) :: &
357 &dtl,dqw,dtv,gm,gh,sm,sh
359 INTEGER :: k
361 REAL :: rfc,f1,f2,rf1,rf2,smc,shc,&
362 &ri1,ri2,ri3,ri4,duz,dtz,dqz,vtt,vqq,dtq,dzk,afk,abk,ri,rf
364 ! ev = 2.5e6
365 ! tv0 = 0.61*tref
366 ! tv1 = 1.61*tref
367 ! gtr = 9.81/tref
368 !
369 rfc = g1/( g1+g2 )
370 f1 = b1*( g1-c1 ) +3.0*a2*( 1.0 -c2 )*( 1.0-c5 ) &
371 & +2.0*a1*( 3.0-2.0*c2 )
372 f2 = b1*( g1+g2 ) -3.0*a1*( 1.0 -c2 )
373 rf1 = b1*( g1-c1 )/f1
374 rf2 = b1* g1 /f2
375 smc = a1 /a2* f1/f2
376 shc = 3.0*a2*( g1+g2 )
377 !
378 ri1 = 0.5/smc
379 ri2 = rf1*smc
380 ri3 = 4.0*rf2*smc -2.0*ri2
381 ri4 = ri2**2
382 !
383 DO k = kts+1,kte
384 dzk = 0.5 *( dz(k)+dz(k-1) )
385 afk = dz(k)/( dz(k)+dz(k-1) )
386 abk = 1.0 -afk
387 duz = ( u(k)-u(k-1) )**2 +( v(k)-v(k-1) )**2
388 duz = duz /dzk**2
389 dtz = ( thl(k)-thl(k-1) )/( dzk )
390 dqz = ( qw(k)-qw(k-1) )/( dzk )
391 !
392 vtt = 1.0 +vt(k)*abk +vt(k-1)*afk
393 vqq = tv0 +vq(k)*abk +vq(k-1)*afk
394 dtq = vtt*dtz +vqq*dqz
395 !
396 dtl(k) = dtz
397 dqw(k) = dqz
398 dtv(k) = dtq
399 !? dtv(i,j,k) = dtz +tv0*dqz
400 !? : +( ev/pi0(i,j,k)-tv1 )
401 !? : *( ql(i,j,k)-ql(i,j,k-1) )/( dzk*h(i,j) )
402 !
403 gm (k) = duz
404 gh (k) = -dtq*gtr
405 !
406 ! ** Gradient Richardson number **
407 ri = -gh(k)/MAX( duz, 1.0e-10 )
408 ! ** Flux Richardson number **
409 rf = MIN( ri1*( ri+ri2-SQRT(ri**2-ri3*ri+ri4) ), rfc )
410 !
411 sh (k) = shc*( rfc-rf )/( 1.0-rf )
412 sm (k) = smc*( rf1-rf )/( rf2-rf ) * sh(k)
413 END DO
414 !
415 RETURN
417 END SUBROUTINE mym_level2
419 !
420 ! SUBROUTINE mym_length:
421 !
422 ! Input variables: see subroutine mym_initialize
423 !
424 ! Output variables: see subroutine mym_initialize
425 !
426 ! Work arrays:
427 ! elt(mx,ny) : Length scale depending on the PBL depth (m)
428 ! vsc(mx,ny) : Velocity scale q_c (m/s)
429 ! at first, used for computing elt
430 !
431 SUBROUTINE mym_length ( kts,kte,&
432 & dz, zw, &
433 & rmo, flt, flq, &
434 & vt, vq, &
435 & qke, &
436 & dtv, &
437 & el, &
438 & qkw)
440 INTEGER, INTENT(IN) :: kts,kte
441 REAL, DIMENSION(kts:kte), INTENT(in) :: dz
442 REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw
443 REAL, INTENT(in) :: rmo,flt,flq
444 REAL, DIMENSION(kts:kte), INTENT(IN) :: qke,vt,vq
447 REAL, DIMENSION(kts:kte), INTENT(out) :: &
448 &qkw, el
449 REAL, DIMENSION(kts:kte), INTENT(in) :: dtv
451 REAL :: elt,vsc
453 INTEGER :: i,j,k
454 REAL :: afk,abk,zwk,dzk,qdz,vflx,bv,elb,els,elf
456 ! tv0 = 0.61*tref
457 ! gtr = 9.81/tref
458 !
459 DO k = kts+1,kte
460 afk = dz(k)/( dz(k)+dz(k-1) )
461 abk = 1.0 -afk
462 qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*&
463 &afk,1.0e-10))
464 END DO
465 !
466 elt = 1.0e-5
467 vsc = 1.0e-5
468 !
469 ! ** Strictly, zwk*h(i,j) -> ( zwk*h(i,j)+z0 ) **
470 DO k = kts+1,kte-1
471 zwk = zw(k)
472 dzk = 0.5*( dz(k)+dz(k-1) )
473 qdz = MAX( qkw(k)-qmin, 0.0 )*dzk
474 elt = elt +qdz*zwk
475 vsc = vsc +qdz
476 END DO
477 !
478 vflx = ( vt(kts)+1.0 )*flt +( vq(kts)+tv0 )*flq
479 elt = alp1*elt/vsc
480 vsc = ( gtr*elt*MAX( vflx, 0.0 ) )**(1.0/3.0)
481 !
482 ! ** Strictly, el(i,j,1) is not zero. **
483 el(kts) = 0.0
484 !
485 DO k = kts+1,kte
486 zwk = zw(k)
487 ! ** Length scale limited by the buoyancy effect **
488 IF ( dtv(k) .GT. 0.0 ) THEN
489 bv = SQRT( gtr*dtv(k) )
490 elb = alp2*qkw(k) / bv &
491 & *( 1.0 + alp3/alp2*&
492 &SQRT( vsc/( bv*elt ) ) )
493 elf=elb
494 elf = alp2 * qkw(k)/bv
495 ELSE
496 elb = 1.0e10
497 elf =elb
498 END IF
499 !
500 ! ** Length scale in the surface layer **
501 IF ( rmo .GT. 0.0 ) THEN
502 els = vk*zwk &
503 & /( 1.0 + cns*MIN( zwk*rmo, zmax ) )
504 ELSE
505 els = vk*zwk &
506 & *( 1.0 - alp4* zwk*rmo )**0.2
507 END IF
508 !
509 ! el(k) = MIN(elb/( elb/elt+elb/els+1.0 ),elf)
510 el(k) = elb/( elb/elt+elb/els+1.0 )
512 END DO
513 !
514 RETURN
516 END SUBROUTINE mym_length
518 !
519 ! SUBROUTINE mym_turbulence:
520 !
521 ! Input variables: see subroutine mym_initialize
522 ! levflag : <>3; Level 2.5
523 ! = 3; Level 3
524 !
525 ! # ql, vt, vq, qke, tsq, qsq and cov are changed to input variables.
526 !
527 ! Output variables: see subroutine mym_initialize
528 ! dfm(mx,my,nz) : Diffusivity coefficient for momentum,
529 ! divided by dz (not dz*h(i,j)) (m/s)
530 ! dfh(mx,my,nz) : Diffusivity coefficient for heat,
531 ! divided by dz (not dz*h(i,j)) (m/s)
532 ! dfq(mx,my,nz) : Diffusivity coefficient for q^2,
533 ! divided by dz (not dz*h(i,j)) (m/s)
534 ! tcd(mx,my,nz) : Countergradient diffusion term for Theta_l
535 ! (K/s)
536 ! qcd(mx,my,nz) : Countergradient diffusion term for Q_w
537 ! (kg/kg s)
538 ! pd?(mx,my,nz) : Half of the production terms
539 !
540 ! Only tcd and qcd are defined at the center of the grid boxes
541 !
542 ! # DO NOT forget that tcd and qcd are added on the right-hand side
543 ! of the equations for Theta_l and Q_w, respectively.
544 !
545 ! Work arrays: see subroutine mym_initialize and level2
546 !
547 ! # dtl, dqw, dtv, gm and gh are allowed to share storage units with
548 ! dfm, dfh, dfq, tcd and qcd, respectively, for saving memory.
549 !
550 SUBROUTINE mym_turbulence ( kts,kte,&
551 & levflag, &
552 & dz, zw, &
553 & u, v, thl, ql, qw, &
554 & qke, tsq, qsq, cov, &
555 & vt, vq,&
556 & rmo, flt, flq, &
557 & El, Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc)
559 !
560 INTEGER, INTENT(IN) :: kts,kte
561 INTEGER, INTENT(IN) :: levflag
562 REAL, DIMENSION(kts:kte), INTENT(in) :: dz
563 REAL, DIMENSION(kts:kte+1), INTENT(in) :: zw
564 REAL, INTENT(in) :: rmo,flt,flq
565 REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,thl,qw,&
566 &ql,vt,vq,qke,tsq,qsq,cov
568 REAL, DIMENSION(kts:kte), INTENT(out) :: dfm,dfh,dfq,&
569 &pdk,pdt,pdq,pdc,tcd,qcd,el
571 REAL, DIMENSION(kts:kte) :: qkw,dtl,dqw,dtv,gm,gh,sm,sh
573 INTEGER :: k
574 ! REAL :: cc2,cc3,e1c,e2c,e3c,e4c,e5c
575 REAL :: e6c,dzk,afk,abk,vtt,vqq,&
576 &cw25,clow,cupp,gamt,gamq,smd,gamv,elq,elh
578 DOUBLE PRECISION q2sq, t2sq, r2sq, c2sq, elsq, gmel, ghel
579 DOUBLE PRECISION q3sq, t3sq, r3sq, c3sq, dlsq, qdiv
580 DOUBLE PRECISION e1, e2, e3, e4, enum, eden, wden
581 !
582 ! tv0 = 0.61*tref
583 ! gtr = 9.81/tref
584 !
585 ! cc2 = 1.0-c2
586 ! cc3 = 1.0-c3
587 ! e1c = 3.0*a2*b2*cc3
588 ! e2c = 9.0*a1*a2*cc2
589 ! e3c = 9.0*a2*a2*cc2*( 1.0-c5 )
590 ! e4c = 12.0*a1*a2*cc2
591 ! e5c = 6.0*a1*a1
592 !
594 CALL mym_level2 (kts,kte,&
595 & dz, &
596 & u, v, thl, qw, &
597 & ql, vt, vq, &
598 & dtl, dqw, dtv, gm, gh, sm, sh )
599 !
600 CALL mym_length (kts,kte, &
601 & dz, zw, &
602 & rmo, flt, flq, &
603 & vt, vq, &
604 & qke, &
605 & dtv, &
606 & el, &
607 & qkw)
608 !
609 DO k = kts+1,kte
610 dzk = 0.5 *( dz(k)+dz(k-1) )
611 afk = dz(k)/( dz(k)+dz(k-1) )
612 abk = 1.0 -afk
613 elsq = el (k)**2
614 q2sq = b1*elsq*( sm(k)*gm(k)+sh(k)*gh(k) )
615 q3sq = qkw(k)**2
616 !
617 ! Modified: Dec/22/2005, from here, (dlsq -> elsq)
618 gmel = gm (k)*elsq
619 ghel = gh (k)*elsq
620 ! Modified: Dec/22/2005, up to here
621 !
622 ! ** Since qkw is set to more than 0.0, q3sq > 0.0. **
623 IF ( q3sq .LT. q2sq ) THEN
624 qdiv = SQRT( q3sq/q2sq )
625 sm(k) = sm(k) * qdiv
626 sh(k) = sh(k) * qdiv
627 !
628 e1 = q3sq - e1c*ghel * qdiv**2
629 e2 = q3sq - e2c*ghel * qdiv**2
630 e3 = e1 + e3c*ghel * qdiv**2
631 e4 = e1 - e4c*ghel * qdiv**2
632 eden = e2*e4 + e3*e5c*gmel * qdiv**2
633 eden = MAX( eden, 1.0d-20 )
634 ELSE
635 e1 = q3sq - e1c*ghel
636 e2 = q3sq - e2c*ghel
637 e3 = e1 + e3c*ghel
638 e4 = e1 - e4c*ghel
639 eden = e2*e4 + e3*e5c*gmel
640 eden = MAX( eden, 1.0d-20 )
641 !
642 qdiv = 1.0
643 sm(k) = q3sq*a1*( e3-3.0*c1*e4 )/eden
644 sh(k) = q3sq*a2*( e2+3.0*c1*e5c*gmel )/eden
645 END IF
646 !
647 ! ** Level 3 : start **
648 IF ( levflag .EQ. 3 ) THEN
649 t2sq = qdiv*b2*elsq*sh(k)*dtl(k)**2
650 r2sq = qdiv*b2*elsq*sh(k)*dqw(k)**2
651 c2sq = qdiv*b2*elsq*sh(k)*dtl(k)*dqw(k)
652 t3sq = MAX( tsq(k)*abk+tsq(k-1)*afk, 0.0 )
653 r3sq = MAX( qsq(k)*abk+qsq(k-1)*afk, 0.0 )
654 c3sq = cov(k)*abk+cov(k-1)*afk
655 !
656 ! Modified: Dec/22/2005, from here
657 c3sq = SIGN( MIN( ABS(c3sq), SQRT(t3sq*r3sq) ), c3sq )
658 !
659 vtt = 1.0 +vt(k)*abk +vt(k-1)*afk
660 vqq = tv0 +vq(k)*abk +vq(k-1)*afk
661 t2sq = vtt*t2sq +vqq*c2sq
662 r2sq = vtt*c2sq +vqq*r2sq
663 c2sq = MAX( vtt*t2sq+vqq*r2sq, 0.0d0 )
664 t3sq = vtt*t3sq +vqq*c3sq
665 r3sq = vtt*c3sq +vqq*r3sq
666 c3sq = MAX( vtt*t3sq+vqq*r3sq, 0.0d0 )
667 !
668 cw25 = e1*( e2 + 3.0*c1*e5c*gmel*qdiv**2 )/( 3.0*eden )
669 !
670 ! ** Limitation on q, instead of L/q **
671 dlsq = elsq
672 IF ( q3sq/dlsq .LT. -gh(k) ) q3sq = -dlsq*gh(k)
673 !
674 ! ** Limitation on c3sq (0.12 =< cw =< 0.76) **
675 e2 = q3sq - e2c*ghel * qdiv**2
676 e3 = q3sq + e3c*ghel * qdiv**2
677 e4 = q3sq - e4c*ghel * qdiv**2
678 eden = e2*e4 + e3 *e5c*gmel * qdiv**2
679 !
680 wden = cc3*gtr**2 * dlsq**2/elsq * qdiv**2 &
681 & *( e2*e4c - e3c*e5c*gmel * qdiv**2 )
682 !
683 IF ( wden .NE. 0.0 ) THEN
684 clow = q3sq*( 0.12-cw25 )*eden/wden
685 cupp = q3sq*( 0.76-cw25 )*eden/wden
686 !
687 IF ( wden .GT. 0.0 ) THEN
688 c3sq = MIN( MAX( c3sq, c2sq+clow ), c2sq+cupp )
689 ELSE
690 c3sq = MAX( MIN( c3sq, c2sq+clow ), c2sq+cupp )
691 END IF
692 END IF
693 !
694 e1 = e2 + e5c*gmel * qdiv**2
695 eden = MAX( eden, 1.0d-20 )
696 ! Modified: Dec/22/2005, up to here
697 !
698 e6c = 3.0*a2*cc3*gtr * dlsq/elsq
700 IF ( c3sq .GT. c2sq ) THEN
701 enum = qdiv*e6c*( t3sq-t2sq )
702 gamt =-e1 *enum /eden
704 enum = qdiv*e6c*( r3sq-r2sq )
705 gamq =-e1 *enum /eden
707 enum = qdiv*e6c*( c3sq-c2sq )
708 smd = dlsq*enum*gtr/eden * qdiv**2 * (e3c+e4c)*a1/a2
709 gamv = e1 *enum*gtr/eden
711 ELSE
712 gamt = 0.0
713 gamq = 0.0
714 smd = 0.0
715 gamv = 0.0
716 END IF
718 !
719 ! ** for Gamma_theta **
720 ! enum = qdiv*e6c*( t3sq-t2sq )
721 ! enum = MAX( qdiv*e6c*( t3sq-t2sq ), 0.0 )
724 ! gamt =-e1 *enum /eden
725 !
726 ! ** for Gamma_q **
727 ! enum = qdiv*e6c*( r3sq-r2sq )
728 ! enum = MAX( qdiv*e6c*( r3sq-r2sq ), 0.0 )
730 ! gamq =-e1 *enum /eden
731 !
732 ! ** for Sm' and Sh'd(Theta_V)/dz **
733 ! enum = qdiv*e6c*( c3sq-c2sq )
734 ! enum = MAX( qdiv*e6c*( c3sq-c2sq ), 0.0 )
736 ! smd = dlsq*enum*gtr/eden * qdiv**2 * (e3c+e4c)*a1/a2
737 ! gamv = e1 *enum*gtr/eden
738 !
740 sm(k) = sm(k) +smd
741 !
742 ! ** For elh (see below), qdiv at Level 3 is reset to 1.0. **
743 qdiv = 1.0
744 ! ** Level 3 : end **
745 !
746 ELSE
747 ! ** At Level 2.5, qdiv is not reset. **
748 gamt = 0.0
749 gamq = 0.0
750 gamv = 0.0
751 END IF
752 !
753 elq = el(k)*qkw(k)
754 elh = elq*qdiv
755 !
756 pdk(k) = elq*( sm(k)*gm (k) &
757 & +sh(k)*gh (k)+gamv )
758 pdt(k) = elh*( sh(k)*dtl(k)+gamt )*dtl(k)
759 pdq(k) = elh*( sh(k)*dqw(k)+gamq )*dqw(k)
760 pdc(k) = elh*( sh(k)*dtl(k)+gamt )&
761 &*dqw(k)*0.5 &
762 &+elh*( sh(k)*dqw(k)+gamq )*dtl(k)*0.5
763 !
764 tcd(k) = elq*gamt
765 qcd(k) = elq*gamq
766 !
767 dfm(k) = elq*sm (k) / dzk
768 dfh(k) = elq*sh (k) / dzk
769 ! Modified: Dec/22/2005, from here
770 ! ** In sub.mym_predict, dfq for the TKE and scalar variance **
771 ! ** are set to 3.0*dfm and 1.0*dfm, respectively. **
772 dfq(k) = dfm(k)
773 ! Modified: Dec/22/2005, up to here
774 END DO
775 !
776 dfm(kts) = 0.0
777 dfh(kts) = 0.0
778 dfq(kts) = 0.0
779 tcd(kts) = 0.0
780 qcd(kts) = 0.0
782 !
783 DO k = kts,kte-1
784 dzk = dz(k)
785 tcd(k) = ( tcd(k+1)-tcd(k) )/( dzk )
786 qcd(k) = ( qcd(k+1)-qcd(k) )/( dzk )
787 END DO
788 !
789 RETURN
791 END SUBROUTINE mym_turbulence
793 !
794 ! SUBROUTINE mym_predict:
795 !
796 ! Input variables: see subroutine mym_initialize and turbulence
797 ! qke(mx,my,nz) : qke at (n)th time level
798 ! tsq, ...cov : ditto
799 !
800 ! Output variables:
801 ! qke(mx,my,nz) : qke at (n+1)th time level
802 ! tsq, ...cov : ditto
803 !
804 ! Work arrays:
805 ! qkw(mx,my,nz) : q at the center of the grid boxes (m/s)
806 ! bp (mx,my,nz) : = 1/2*F, see below
807 ! rp (mx,my,nz) : = P-1/2*F*Q, see below
808 !
809 ! # The equation for a turbulent quantity Q can be expressed as
810 ! dQ/dt + Ah + Av = Dh + Dv + P - F*Q, (1)
811 ! where A is the advection, D the diffusion, P the production,
812 ! F*Q the dissipation and h and v denote horizontal and vertical,
813 ! respectively. If Q is q^2, F is 2q/B_1L.
814 ! Using the Crank-Nicholson scheme for Av, Dv and F*Q, a finite
815 ! difference equation is written as
816 ! Q{n+1} - Q{n} = dt *( Dh{n} - Ah{n} + P{n} )
817 ! + dt/2*( Dv{n} - Av{n} - F*Q{n} )
818 ! + dt/2*( Dv{n+1} - Av{n+1} - F*Q{n+1} ), (2)
819 ! where n denotes the time level.
820 ! When the advection and diffusion terms are discretized as
821 ! dt/2*( Dv - Av ) = a(k)Q(k+1) - b(k)Q(k) + c(k)Q(k-1), (3)
822 ! Eq.(2) can be rewritten as
823 ! - a(k)Q(k+1) + [ 1 + b(k) + dt/2*F ]Q(k) - c(k)Q(k-1)
824 ! = Q{n} + dt *( Dh{n} - Ah{n} + P{n} )
825 ! + dt/2*( Dv{n} - Av{n} - F*Q{n} ), (4)
826 ! where Q on the left-hand side is at (n+1)th time level.
827 !
828 ! In this subroutine, a(k), b(k) and c(k) are obtained from
829 ! subprogram coefvu and are passed to subprogram tinteg via
830 ! common. 1/2*F and P-1/2*F*Q are stored in bp and rp,
831 ! respectively. Subprogram tinteg solves Eq.(4).
832 !
833 ! Modify this subroutine according to your numerical integration
834 ! scheme (program).
835 !
838 SUBROUTINE mym_predict (kts,kte,&
839 & levflag, &
840 & delt,&
841 & dz, &
842 & ust, flt, flq, pmz, phh, &
843 & el, dfq, &
844 & pdk, pdt, pdq, pdc,&
845 & qke, tsq, qsq, cov)
847 INTEGER, INTENT(IN) :: kts,kte
848 INTEGER, INTENT(IN) :: levflag
849 REAL, INTENT(IN) :: delt
850 REAL, DIMENSION(kts:kte), INTENT(IN) :: dz, dfq,el
851 REAL, DIMENSION(kts:kte), INTENT(INOUT) :: pdk, pdt, pdq, pdc
852 REAL, INTENT(IN) :: flt, flq, ust, pmz, phh
853 REAL, DIMENSION(kts:kte), INTENT(INOUT) :: qke,tsq, qsq, cov
855 INTEGER :: k,nz
856 REAL, DIMENSION(kts:kte) :: qkw, bp, rp, df3q
857 REAL :: vkz,pdk1,phm,pdt1,pdq1,pdc1,b1l,b2l
858 REAL, DIMENSION(kts:kte) :: dtz
859 REAL, DIMENSION(1:kte-kts+1) :: a,b,c,d
861 nz=kte-kts+1
863 ! ** Strictly, vkz*h(i,j) -> vk*( 0.5*dz(1)*h(i,j)+z0 ) **
864 vkz = vk*0.5*dz(kts)
865 !
866 ! Modified: Dec/22/2005, from here
867 ! ** dfq for the TKE is 3.0*dfm. **
868 ! CALL coefvu ( dfq, 3.0 ) ! make change here
869 ! Modified: Dec/22/2005, up to here
870 !
871 DO k = kts,kte
872 !! qke(k) = MAX(qke(k), 0.0)
873 qkw(k) = SQRT( MAX( qke(k), 0.0 ) )
874 df3q(k)=3.*dfq(k)
875 dtz(k)=delt/dz(k)
876 END DO
877 !
878 pdk1 = 2.0*ust**3*pmz/( vkz )
879 phm = 2.0/ust *phh/( vkz )
880 pdt1 = phm*flt**2
881 pdq1 = phm*flq**2
882 pdc1 = phm*flt*flq
883 !
884 ! ** pdk(i,j,1)+pdk(i,j,2) corresponds to pdk1. **
885 pdk(kts) = pdk1 -pdk(kts+1)
886 pdt(kts) = pdt1 -pdt(kts+1)
887 pdq(kts) = pdq1 -pdq(kts+1)
888 pdc(kts) = pdc1 -pdc(kts+1)
889 !
890 ! ** Prediction of twice the turbulent kinetic energy **
891 !! DO k = kts+1,kte-1
892 DO k = kts,kte-1
893 b1l = b1*0.5*( el(k+1)+el(k) )
894 bp(k) = 2.*qkw(k) / b1l
895 rp(k) = pdk(k+1) + pdk(k)
896 END DO
898 !! a(1)=0.
899 !! b(1)=1.
900 !! c(1)=-1.
901 !! d(1)=0.
903 ! Since df3q(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)*df3q(k+1)+bp(k)*delt.
904 DO k=kts,kte-1
905 a(k-kts+1)=-dtz(k)*df3q(k)
906 b(k-kts+1)=1.+dtz(k)*(df3q(k)+df3q(k+1))+bp(k)*delt
907 c(k-kts+1)=-dtz(k)*df3q(k+1)
908 d(k-kts+1)=rp(k)*delt + qke(k)
909 ENDDO
911 !! DO k=kts+1,kte-1
912 !! a(k-kts+1)=-dtz(k)*df3q(k)
913 !! b(k-kts+1)=1.+dtz(k)*(df3q(k)+df3q(k+1))
914 !! c(k-kts+1)=-dtz(k)*df3q(k+1)
915 !! d(k-kts+1)=rp(k)*delt + qke(k) - qke(k)*bp(k)*delt
916 !! ENDDO
918 a(nz)=-1. !0.
919 b(nz)=1.
920 c(nz)=0.
921 d(nz)=0.
923 CALL tridiag(nz,a,b,c,d)
925 DO k=kts,kte
926 qke(k)=d(k-kts+1)
927 ENDDO
930 IF ( levflag .EQ. 3 ) THEN
931 !
932 ! Modified: Dec/22/2005, from here
933 ! ** dfq for the scalar variance is 1.0*dfm. **
934 ! CALL coefvu ( dfq, 1.0 ) make change here
935 ! Modified: Dec/22/2005, up to here
936 !
937 ! ** Prediction of the temperature variance **
938 !! DO k = kts+1,kte-1
939 DO k = kts,kte-1
940 b2l = b2*0.5*( el(k+1)+el(k) )
941 bp(k) = 2.*qkw(k) / b2l
942 rp(k) = pdt(k+1) + pdt(k)
943 END DO
945 !zero gradient for tsq at bottom and top
947 !! a(1)=0.
948 !! b(1)=1.
949 !! c(1)=-1.
950 !! d(1)=0.
952 ! Since dfq(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)*dfq(k+1)+bp(k)*delt.
953 DO k=kts,kte-1
954 a(k-kts+1)=-dtz(k)*dfq(k)
955 b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))+bp(k)*delt
956 c(k-kts+1)=-dtz(k)*dfq(k+1)
957 d(k-kts+1)=rp(k)*delt + tsq(k)
958 ENDDO
960 !! DO k=kts+1,kte-1
961 !! a(k-kts+1)=-dtz(k)*dfq(k)
962 !! b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))
963 !! c(k-kts+1)=-dtz(k)*dfq(k+1)
964 !! d(k-kts+1)=rp(k)*delt + tsq(k) - tsq(k)*bp(k)*delt
965 !! ENDDO
967 a(nz)=-1. !0.
968 b(nz)=1.
969 c(nz)=0.
970 d(nz)=0.
972 CALL tridiag(nz,a,b,c,d)
974 DO k=kts,kte
975 tsq(k)=d(k-kts+1)
976 ENDDO
978 ! ** Prediction of the moisture variance **
979 !! DO k = kts+1,kte-1
980 DO k = kts,kte-1
981 b2l = b2*0.5*( el(k+1)+el(k) )
982 bp(k) = 2.*qkw(k) / b2l
983 rp(k) = pdq(k+1) +pdq(k)
984 END DO
986 !zero gradient for qsq at bottom and top
988 !! a(1)=0.
989 !! b(1)=1.
990 !! c(1)=-1.
991 !! d(1)=0.
993 ! Since dfq(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)*dfq(k+1)+bp(k)*delt.
994 DO k=kts,kte-1
995 a(k-kts+1)=-dtz(k)*dfq(k)
996 b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))+bp(k)*delt
997 c(k-kts+1)=-dtz(k)*dfq(k+1)
998 d(k-kts+1)=rp(k)*delt + qsq(k)
999 ENDDO
1001 !! DO k=kts+1,kte-1
1002 !! a(k-kts+1)=-dtz(k)*dfq(k)
1003 !! b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))
1004 !! c(k-kts+1)=-dtz(k)*dfq(k+1)
1005 !! d(k-kts+1)=rp(k)*delt + qsq(k) -qsq(k)*bp(k)*delt
1006 !! ENDDO
1008 a(nz)=-1. !0.
1009 b(nz)=1.
1010 c(nz)=0.
1011 d(nz)=0.
1013 CALL tridiag(nz,a,b,c,d)
1015 DO k=kts,kte
1016 qsq(k)=d(k-kts+1)
1017 ENDDO
1019 ! ** Prediction of the temperature-moisture covariance **
1020 !! DO k = kts+1,kte-1
1021 DO k = kts,kte-1
1022 b2l = b2*0.5*( el(k+1)+el(k) )
1023 bp(k) = 2.*qkw(k) / b2l
1024 rp(k) = pdc(k+1) + pdc(k)
1025 END DO
1027 !zero gradient for tqcov at bottom and top
1029 !! a(1)=0.
1030 !! b(1)=1.
1031 !! c(1)=-1.
1032 !! d(1)=0.
1034 ! Since dfq(kts)=0.0, a(1)=0.0 and b(1)=1.+dtz(k)*dfq(k+1)+bp(k)*delt.
1035 DO k=kts,kte-1
1036 a(k-kts+1)=-dtz(k)*dfq(k)
1037 b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))+bp(k)*delt
1038 c(k-kts+1)=-dtz(k)*dfq(k+1)
1039 d(k-kts+1)=rp(k)*delt + cov(k)
1040 ENDDO
1042 !! DO k=kts+1,kte-1
1043 !! a(k-kts+1)=-dtz(k)*dfq(k)
1044 !! b(k-kts+1)=1.+dtz(k)*(dfq(k)+dfq(k+1))
1045 !! c(k-kts+1)=-dtz(k)*dfq(k+1)
1046 !! d(k-kts+1)=rp(k)*delt + cov(k) - cov(k)*bp(k)*delt
1047 !! ENDDO
1049 a(nz)=-1. !0.
1050 b(nz)=1.
1051 c(nz)=0.
1052 d(nz)=0.
1054 CALL tridiag(nz,a,b,c,d)
1056 DO k=kts,kte
1057 cov(k)=d(k-kts+1)
1058 ENDDO
1060 ELSE
1061 !! DO k = kts+1,kte-1
1062 DO k = kts,kte-1
1063 IF ( qkw(k) .LE. 0.0 ) THEN
1064 b2l = 0.0
1065 ELSE
1066 b2l = b2*0.25*( el(k+1)+el(k) )/qkw(k)
1067 END IF
1068 !
1069 tsq(k) = b2l*( pdt(k+1)+pdt(k) )
1070 qsq(k) = b2l*( pdq(k+1)+pdq(k) )
1071 cov(k) = b2l*( pdc(k+1)+pdc(k) )
1072 END DO
1074 !! tsq(kts)=tsq(kts+1)
1075 !! qsq(kts)=qsq(kts+1)
1076 !! cov(kts)=cov(kts+1)
1078 tsq(kte)=tsq(kte-1)
1079 qsq(kte)=qsq(kte-1)
1080 cov(kte)=cov(kte-1)
1082 END IF
1084 END SUBROUTINE mym_predict
1086 ! SUBROUTINE mym_condensation:
1087 !
1088 ! Input variables: see subroutine mym_initialize and turbulence
1089 ! pi (mx,my,nz) : Perturbation of the Exner function (J/kg K)
1090 ! defined on the walls of the grid boxes
1091 ! This is usually computed by integrating
1092 ! d(pi)/dz = h*g*tv/tref**2
1093 ! from the upper boundary, where tv is the
1094 ! virtual potential temperature minus tref.
1095 !
1096 ! Output variables: see subroutine mym_initialize
1097 ! cld(mx,my,nz) : Cloud fraction
1098 !
1099 ! Work arrays:
1100 ! qmq(mx,my,nz) : Q_w-Q_{sl}, where Q_{sl} is the saturation
1101 ! specific humidity at T=Tl
1102 ! alp(mx,my,nz) : Functions in the condensation process
1103 ! bet(mx,my,nz) : ditto
1104 ! sgm(mx,my,nz) : Combined standard deviation sigma_s
1105 ! multiplied by 2/alp
1106 !
1107 ! # qmq, alp, bet and sgm are allowed to share storage units with
1108 ! any four of other work arrays for saving memory.
1109 !
1110 ! # Results are sensitive particularly to values of cp and rd.
1111 ! Set these values to those adopted by you.
1112 !
1116 SUBROUTINE mym_condensation (kts,kte, &
1117 & levflag, &
1118 & dz, &
1119 & thl, qw, &
1120 & p,exner, &
1121 & tsq, qsq, cov, &
1122 & Vt, Vq)
1124 INTEGER, INTENT(IN) :: kts,kte
1125 INTEGER, INTENT(in) :: levflag
1127 REAL, DIMENSION(kts:kte), INTENT(IN) :: dz
1128 REAL, DIMENSION(kts:kte), INTENT(IN) :: p,exner, thl, qw, &
1129 &tsq, qsq, cov
1131 REAL, DIMENSION(kts:kte), INTENT(OUT) :: vt,vq
1134 REAL, DIMENSION(kts:kte) :: qmq,alp,bet,sgm,ql,cld
1136 DOUBLE PRECISION :: t3sq, r3sq, c3sq
1137 !
1139 REAL :: p2a,t,esl,qsl,dqsl,q1,cld0,eq1,qll,&
1140 &q2p,pt,rac,qt
1141 INTEGER :: i,j,k
1143 REAL :: erf
1145 ! Note: kte needs to be larger than kts, i.e., kte >= kts+1.
1147 !! DO k = kts,kte-1
1148 DO k = kts,kte
1149 p2a = exner(k)
1150 t = thl(k)*p2a
1152 !x if ( ct .gt. 0.0 ) then
1153 ! a = 17.27
1154 ! b = 237.3
1155 !x else
1156 !x a = 21.87
1157 !x b = 265.5
1158 !x end if
1159 !
1160 ! ** 3.8 = 0.622*6.11 (hPa) **
1161 esl=svp11*EXP(svp2*(t-svpt0)/(t-svp3))
1162 qsl=ep_2*esl/(p(k)-ep_3*esl)
1163 ! qsl = 3.8*EXP( a*ct/( b+ct ) ) / ( 1000.0*p2a**rk )
1164 dqsl = qsl*SVP11*ev/( rd*t**2 )
1165 !
1166 qmq(k) = qw(k) -qsl
1168 alp(k) = 1.0/( 1.0+dqsl*xlvcp )
1169 bet(k) = dqsl*p2a
1170 !
1171 t3sq = MAX( tsq(k), 0.0 )
1172 r3sq = MAX( qsq(k), 0.0 )
1173 c3sq = cov(k)
1174 c3sq = SIGN( MIN( ABS(c3sq), SQRT(t3sq*r3sq) ), c3sq )
1175 !
1176 r3sq = r3sq +bet(k)**2*t3sq -2.0*bet(k)*c3sq
1177 sgm(k) = SQRT( MAX( r3sq, 1.0d-10 ) )
1178 END DO
1179 !
1181 IF ( levflag .NE. 3 ) THEN
1182 ! ** For some reason, sgm(i,j,1)=sgm(i,j,2) seems to be better. **
1183 sgm(kts) = sgm(kts+1)
1184 END IF
1185 !
1186 DO k = kts,kte-1
1187 q1 = qmq(k) / sgm(k)
1188 cld0 = 0.5*( 1.0+erf( q1*rr2 ) )
1189 ! q1=0.
1190 ! cld0=0.
1192 eq1 = rrp*EXP( -0.5*q1*q1 )
1193 qll = MAX( cld0*q1 + eq1, 0.0 )
1195 cld(k) = cld0
1196 ql (k) = alp(k)*sgm(k)*qll
1197 !
1198 q2p = xlvcp/exner( k )
1199 pt = thl(k) +q2p*ql(k)
1200 qt = 1.0 +p608*qw(k) -(1.+p608)*ql(k)
1201 rac = alp(k)*( cld0-qll*eq1 )*( q2p*qt-(1.+p608)*pt )
1202 !
1203 vt (k) = qt-1.0 -rac*bet(k)
1204 vq (k) = p608*pt-tv0 +rac
1205 END DO
1206 !
1208 cld(kte) = cld(kte-1)
1209 ql(kte) = ql(kte-1)
1210 vt(kte) = vt(kte-1)
1211 vq(kte) = vq(kte-1)
1213 RETURN
1215 END SUBROUTINE mym_condensation
1217 SUBROUTINE mynn_tendencies(kts,kte,&
1218 &levflag,grav_settling,&
1219 &delt,&
1220 &dz,&
1221 &u,v,th,qv,qc,p,exner,&
1222 &thl,sqv,sqc,sqw,&
1223 &ust,flt,flq,wspd,qcg,&
1224 &tsq,qsq,cov,&
1225 &tcd,qcd,&
1226 &dfm,dfh,dfq,&
1227 &Du,Dv,Dth,Dqv,Dqc)
1229 INTEGER, INTENT(in) :: kts,kte
1230 INTEGER, INTENT(in) :: grav_settling,levflag
1232 !! grav_settling = 1 for gravitational settling of droplets
1233 !! grav_settling = 0 otherwise
1234 ! thl - liquid water potential temperature
1235 ! qw - total water
1236 ! dfm,dfh,dfq - as above
1237 ! flt - surface flux of thl
1238 ! flq - surface flux of qw
1241 REAL, DIMENSION(kts:kte), INTENT(in) :: u,v,th,qv,qc,p,exner,&
1242 &dfm,dfh,dfq,dz,tsq,qsq,cov,tcd,qcd
1243 REAL, DIMENSION(kts:kte), INTENT(inout) :: thl,sqw,sqv,sqc
1244 REAL, DIMENSION(kts:kte), INTENT(out) :: du,dv,dth,dqv,dqc
1245 REAL, INTENT(IN) :: delt,ust,flt,flq,wspd,qcg
1247 ! REAL, INTENT(IN) :: delt,ust,flt,flq,qcg,&
1248 ! &gradu_top,gradv_top,gradth_top,gradqv_top
1251 !local vars
1253 REAL, DIMENSION(kts:kte) :: dtz,vt,vq
1255 REAL, DIMENSION(1:kte-kts+1) :: a,b,c,d
1257 REAL :: rhs,gfluxm,gfluxp,dztop
1258 INTEGER :: k,kk,nz
1260 nz=kte-kts+1
1262 dztop=.5*(dz(kte)+dz(kte-1))
1264 DO k=kts,kte
1265 dtz(k)=delt/dz(k)
1266 ENDDO
1268 !! u
1270 k=kts
1272 a(1)=0.
1273 b(1)=1.+dtz(k)*(dfm(k+1)+ust**2/wspd)
1274 c(1)=-dtz(k)*dfm(k+1)
1275 d(1)=u(k)
1277 !! a(1)=0.
1278 !! b(1)=1.+dtz(k)*dfm(k+1)
1279 !! c(1)=-dtz(k)*dfm(k+1)
1280 !! d(1)=u(k)*(1.-ust**2/wspd*dtz(k))
1282 DO k=kts+1,kte-1
1283 kk=k-kts+1
1284 a(kk)=-dtz(k)*dfm(k)
1285 b(kk)=1.+dtz(k)*(dfm(k)+dfm(k+1))
1286 c(kk)=-dtz(k)*dfm(k+1)
1287 d(kk)=u(k)
1288 ENDDO
1290 !! no flux at the top
1292 ! a(nz)=-1.
1293 ! b(nz)=1.
1294 ! c(nz)=0.
1295 ! d(nz)=0.
1297 !! specified gradient at the top
1299 ! a(nz)=-1.
1300 ! b(nz)=1.
1301 ! c(nz)=0.
1302 ! d(nz)=gradu_top*dztop
1304 !! prescribed value
1306 a(nz)=0
1307 b(nz)=1.
1308 c(nz)=0.
1309 d(nz)=u(kte)
1311 CALL tridiag(nz,a,b,c,d)
1313 DO k=kts,kte
1314 du(k)=(d(k-kts+1)-u(k))/delt
1315 ENDDO
1317 !! v
1319 k=kts
1321 a(1)=0.
1322 b(1)=1.+dtz(k)*(dfm(k+1)+ust**2/wspd)
1323 c(1)=-dtz(k)*dfm(k+1)
1324 d(1)=v(k)
1326 !! a(1)=0.
1327 !! b(1)=1.+dtz(k)*dfm(k+1)
1328 !! c(1)=-dtz(k)*dfm(k+1)
1329 !! d(1)=v(k)*(1.-ust**2/wspd*dtz(k))
1331 DO k=kts+1,kte-1
1332 kk=k-kts+1
1333 a(kk)=-dtz(k)*dfm(k)
1334 b(kk)=1.+dtz(k)*(dfm(k)+dfm(k+1))
1335 c(kk)=-dtz(k)*dfm(k+1)
1336 d(kk)=v(k)
1337 ENDDO
1339 !! no flux at the top
1341 ! a(nz)=-1.
1342 ! b(nz)=1.
1343 ! c(nz)=0.
1344 ! d(nz)=0.
1347 !! specified gradient at the top
1349 ! a(nz)=-1.
1350 ! b(nz)=1.
1351 ! c(nz)=0.
1352 ! d(nz)=gradv_top*dztop
1354 !! prescribed value
1356 a(nz)=0
1357 b(nz)=1.
1358 c(nz)=0.
1359 d(nz)=v(kte)
1361 CALL tridiag(nz,a,b,c,d)
1363 DO k=kts,kte
1364 dv(k)=(d(k-kts+1)-v(k))/delt
1365 ENDDO
1367 !! thl
1369 k=kts
1371 a(1)=0.
1372 b(1)=1.+dtz(k)*dfh(k+1)
1373 c(1)=-dtz(k)*dfh(k+1)
1375 ! if qcg not used then assume constant flux in the surface layer
1377 IF (qcg < qcgmin) THEN
1378 gfluxm=grav_settling*gno*sqc(k)**gpw
1379 ELSE
1380 gfluxm=grav_settling*gno*(qcg/(1.+qcg))**gpw
1381 ENDIF
1383 gfluxp=grav_settling*gno*(.5*(sqc(k+1)+sqc(k)))**gpw
1385 rhs=-xlvcp/exner(k)&
1386 &*( &
1387 &(gfluxp - gfluxm)/dz(k)&
1388 & ) + tcd(k)
1390 d(1)=thl(k)+dtz(k)*flt+rhs*delt
1392 DO k=kts+1,kte-1
1393 kk=k-kts+1
1394 a(kk)=-dtz(k)*dfh(k)
1395 b(kk)=1.+dtz(k)*(dfh(k)+dfh(k+1))
1396 c(kk)=-dtz(k)*dfh(k+1)
1397 gfluxp=grav_settling*gno*(.5*(sqc(k+1)+sqc(k)))**gpw
1398 gfluxm=grav_settling*gno*(.5*(sqc(k-1)+sqc(k)))**gpw
1399 rhs=-xlvcp/exner(k)&
1400 &*( &
1401 &(gfluxp - gfluxm)/dz(k)&
1402 & ) + tcd(k)
1403 d(kk)=thl(k)+rhs*delt
1404 ENDDO
1406 !! no flux at the top
1408 ! a(nz)=-1.
1409 ! b(nz)=1.
1410 ! c(nz)=0.
1411 ! d(nz)=0.
1413 !! specified gradient at the top
1415 !assume gradthl_top=gradth_top
1417 ! a(nz)=-1.
1418 ! b(nz)=1.
1419 ! c(nz)=0.
1420 ! d(nz)=gradth_top*dztop
1422 !! prescribed value
1424 a(nz)=0.
1425 b(nz)=1.
1426 c(nz)=0.
1427 d(nz)=thl(kte)
1429 CALL tridiag(nz,a,b,c,d)
1431 DO k=kts,kte
1432 thl(k)=d(k-kts+1)
1433 ENDDO
1435 !! total water
1437 k=kts
1439 a(1)=0.
1440 b(1)=1.+dtz(k)*dfh(k+1)
1441 c(1)=-dtz(k)*dfh(k+1)
1443 IF (qcg < qcgmin) THEN
1444 gfluxm=grav_settling*gno*sqc(k)**gpw
1445 ELSE
1446 gfluxm=grav_settling*gno*(qcg/(1.+qcg))**gpw
1447 ENDIF
1449 gfluxp=grav_settling*gno*(.5*(sqc(k+1)+sqc(k)))**gpw
1451 rhs=&
1452 &( &
1453 &(gfluxp - gfluxm)/dz(k)&
1454 & ) + qcd(k)
1456 d(1)=sqw(k)+dtz(k)*flq+rhs*delt
1458 DO k=kts+1,kte-1
1459 kk=k-kts+1
1460 a(kk)=-dtz(k)*dfh(k)
1461 b(kk)=1.+dtz(k)*(dfh(k)+dfh(k+1))
1462 c(kk)=-dtz(k)*dfh(k+1)
1463 gfluxp=grav_settling*gno*(.5*(sqc(k+1)+sqc(k)))**gpw
1464 gfluxm=grav_settling*gno*(.5*(sqc(k-1)+sqc(k)))**gpw
1466 rhs=&
1467 &( &
1468 &(gfluxp - gfluxm)/dz(k)&
1469 & ) + qcd(k)
1470 d(kk)=sqw(k) + rhs*delt
1471 ENDDO
1474 !! no flux at the top
1476 ! a(nz)=-1.
1477 ! b(nz)=1.
1478 ! c(nz)=0.
1479 ! d(nz)=0.
1481 !! specified gradient at the top
1482 !assume gradqw_top=gradqv_top
1484 ! a(nz)=-1.
1485 ! b(nz)=1.
1486 ! c(nz)=0.
1487 ! d(nz)=gradqv_top*dztop
1489 !! prescribed value
1491 a(nz)=0.
1492 b(nz)=1.
1493 c(nz)=0.
1494 d(nz)=sqw(kte)
1496 CALL tridiag(nz,a,b,c,d)
1498 !convert to mixing ratios for wrf
1500 DO k=kts,kte
1501 sqw(k)=d(k-kts+1)
1502 sqv(k)=sqw(k)-sqc(k)
1503 Dqv(k)=(sqv(k)/(1.-sqv(k))-qv(k))/delt
1504 ! Dqc(k)=(sqc(k)/(1.-sqc(k))-qc(k))/delt
1505 Dqc(k)=0.
1506 Dth(k)=(thl(k)+xlvcp/exner(k)*sqc(k)-th(k))/delt
1507 ENDDO
1509 END SUBROUTINE mynn_tendencies
1511 SUBROUTINE retrieve_exchange_coeffs(kts,kte,&
1512 &dfm,dfh,dfq,dz,&
1513 &K_m,K_h,K_q)
1515 INTEGER , INTENT(in) :: kts,kte
1517 REAL, DIMENSION(KtS:KtE), INTENT(in) :: dz,dfm,dfh,dfq
1519 REAL, DIMENSION(KtS:KtE), INTENT(out) :: &
1520 &K_m, K_h, K_q
1523 INTEGER :: k
1524 REAL :: dzk
1526 K_m(kts)=0.
1527 K_h(kts)=0.
1528 K_q(kts)=0.
1530 DO k=kts+1,kte
1531 dzk = 0.5 *( dz(k)+dz(k-1) )
1532 K_m(k)=dfm(k)*dzk
1533 K_h(k)=dfh(k)*dzk
1534 K_q(k)=dfq(k)*dzk
1535 ENDDO
1537 END SUBROUTINE retrieve_exchange_coeffs
1540 SUBROUTINE tridiag(n,a,b,c,d)
1542 !! to solve system of linear eqs on tridiagonal matrix n times n
1543 !! after Peaceman and Rachford, 1955
1544 !! a,b,c,d - are vectors of order n
1545 !! a,b,c - are coefficients on the LHS
1546 !! d - is initially RHS on the output becomes a solution vector
1548 INTEGER, INTENT(in):: n
1549 REAL, DIMENSION(n), INTENT(in) :: a,b
1550 REAL, DIMENSION(n), INTENT(inout) :: c,d
1552 INTEGER :: i
1553 REAL :: p
1554 REAL, DIMENSION(n) :: q
1556 c(n)=0.
1557 q(1)=-c(1)/b(1)
1558 d(1)=d(1)/b(1)
1560 DO i=2,n
1561 p=1./(b(i)+a(i)*q(i-1))
1562 q(i)=-c(i)*p
1563 d(i)=(d(i)-a(i)*d(i-1))*p
1564 ENDDO
1566 DO i=n-1,1,-1
1567 d(i)=d(i)+q(i)*d(i+1)
1568 ENDDO
1570 END SUBROUTINE tridiag
1572 SUBROUTINE mynn_bl_driver(&
1573 &initflag,&
1574 &grav_settling,&
1575 &delt,&
1576 &dz,&
1577 &u,v,th,qv,qc,&
1578 &p,exner,rho,&
1579 &xland,ts,qsfc,qcg,ps,&
1580 &ust,ch,hfx,qfx,rmol,wspd,&
1581 &Qke,Tsq,Qsq,Cov,&
1582 &Du,Dv,Dth,&
1583 &Dqv,Dqc,&
1584 ! &K_m,K_h,K_q&
1585 &K_h,k_m&
1586 &,IDS,IDE,JDS,JDE,KDS,KDE &
1587 &,IMS,IME,JMS,JME,KMS,KME &
1588 &,ITS,ITE,JTS,JTE,KTS,KTE)
1590 INTEGER, INTENT(in) :: initflag
1591 INTEGER, INTENT(in) :: grav_settling
1593 INTEGER,INTENT(IN) :: &
1594 & IDS,IDE,JDS,JDE,KDS,KDE &
1595 &,IMS,IME,JMS,JME,KMS,KME &
1596 &,ITS,ITE,JTS,JTE,KTS,KTE
1599 ! initflag > 0 for TRUE
1600 ! else for FALSE
1601 ! levflag : <>3; Level 2.5
1602 ! = 3; Level 3
1603 ! grav_settling = 1 when gravitational settling accounted for
1604 ! grav_settling = 0 when gravitational settling NOT accounted for
1606 REAL, INTENT(in) :: delt
1607 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(in) :: dz,&
1608 &u,v,th,qv,qc,p,exner,rho
1609 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(in) :: xland,ust,&
1610 &ch,rmol,ts,qsfc,qcg,ps,hfx,qfx, wspd
1611 !
1612 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(inout) :: &
1613 &Qke,Tsq,Qsq,Cov
1614 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(out) :: &
1615 &Du,Dv,Dth,Dqv,Dqc
1617 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(out) :: &
1618 &K_h,K_m
1622 !local vars
1623 INTEGER :: ITF,JTF,KTF
1624 INTEGER :: i,j,k
1625 REAL, DIMENSION(KMS:KME) :: thl,sqv,sqc,sqw,&
1626 &El, Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc, Vt, Vq
1628 REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME) :: K_q
1630 REAL, DIMENSION(KMS:KME+1) :: zw
1632 REAL :: cpm,sqcg,flt,flq,pmz,phh,exnerg,zet
1634 INTEGER, SAVE :: levflag
1636 JTF=MIN0(JTE,JDE-1)
1637 ITF=MIN0(ITE,IDE-1)
1638 KTF=MIN0(KTE,KDE-1)
1640 IF (initflag > 0) THEN
1642 levflag=mynn_level
1644 DO j=JTS,JTF
1645 DO i=ITS,ITF
1646 DO k=KTS,KTF
1647 sqv(k)=qv(i,k,j)/(1.+qv(i,k,j))
1648 thl(k)=th(i,k,j)
1650 IF (k==kts) THEN
1651 zw(k)=0.
1652 ELSE
1653 zw(k)=zw(k-1)+dz(i,k-1,j)
1654 ENDIF
1656 k_m(i,k,j)=0.
1657 k_h(i,k,j)=0.
1658 k_q(i,k,j)=0.
1660 ENDDO
1662 zw(ktf+1)=zw(ktf)+dz(i,ktf,j)
1664 CALL mym_initialize ( kts,kte,&
1665 &dz(i,kts:kte,j), zw, &
1666 &u(i,kts:kte,j), v(i,kts:kte,j), &
1667 &thl, sqv,&
1668 &ust(i,j), rmol(i,j),&
1669 &Qke(i,kts:kte,j), Tsq(i,kts:kte,j), &
1670 &Qsq(i,kts:kte,j), Cov(i,kts:kte,j))
1672 ENDDO
1673 ENDDO
1675 ENDIF
1677 DO j=JTS,JTF
1678 DO i=ITS,ITF
1679 DO k=KTS,KTF
1680 sqv(k)=qv(i,k,j)/(1.+qv(i,k,j))
1681 sqc(k)=qc(i,k,j)/(1.+qc(i,k,j))
1682 sqw(k)=sqv(k)+sqc(k)
1683 thl(k)=th(i,k,j)-xlvcp/exner(i,k,j)*sqc(k)
1684 IF (k==kts) THEN
1685 zw(k)=0.
1686 ELSE
1687 zw(k)=zw(k-1)+dz(i,k,j)
1688 ENDIF
1689 ENDDO
1691 zw(ktf+1)=zw(ktf)+dz(i,ktf,j)
1693 sqcg=qcg(i,j)/(1.+qcg(i,j))
1694 cpm=cp*(1.+0.8*qv(i,kts,j))
1696 ! The exchange coefficient for cloud water is assumed to be the same as
1697 ! that for heat. CH is multiplied by WSPD. See module_sf_mynn.F
1698 exnerg=(ps(i,j)/p1000mb)**rcp
1699 flt = hfx(i,j)/( rho(i,kts,j)*cpm ) &
1700 +xlvcp*ch(i,j)*(sqc(kts)/exner(i,kts,j)-sqcg/exnerg)
1701 flq = qfx(i,j)/ rho(i,kts,j) &
1702 -ch(i,j)*(sqc(kts) -sqcg )
1704 !!!!!
1705 zet = 0.5*dz(i,kts,j)*rmol(i,j)
1706 if ( zet >= 0.0 ) then
1707 pmz = 1.0 + (cphm_st-1.0) * zet
1708 phh = 1.0 + cphh_st * zet
1709 else
1710 pmz = 1.0/ (1.0-cphm_unst*zet)**0.25 - zet
1711 phh = 1.0/SQRT(1.0-cphh_unst*zet)
1712 end if
1713 !!!!!
1715 CALL mym_condensation ( kts,kte,&
1716 &levflag, &
1717 &dz(i,kts:kte,j), &
1718 &thl, sqw, &
1719 &p(i,kts:kte,j),exner(i,kts:kte,j), &
1720 &tsq(i,kts:kte,j), qsq(i,kts:kte,j), cov(i,kts:kte,j), &
1721 &Vt, Vq)
1723 CALL mym_turbulence ( kts,kte,&
1724 &levflag, &
1725 &dz(i,kts:kte,j), zw, &
1726 &u(i,kts:kte,j), v(i,kts:kte,j), thl, sqc, sqw, &
1727 &qke(i,kts:kte,j), tsq(i,kts:kte,j), &
1728 &qsq(i,kts:kte,j), cov(i,kts:kte,j), &
1729 &vt, vq,&
1730 &rmol(i,j), flt, flq, &
1731 &El, Dfm, Dfh, Dfq, Tcd, Qcd, Pdk, Pdt, Pdq, Pdc)
1734 CALL mym_predict (kts,kte,&
1735 &levflag, &
1736 &delt,&
1737 &dz(i,kts:kte,j), &
1738 &ust(i,j), flt, flq, pmz, phh, &
1739 &el, dfq, &
1740 &pdk, pdt, pdq, pdc,&
1741 &Qke(i,kts:kte,j), Tsq(i,kts:kte,j), &
1742 &Qsq(i,kts:kte,j), Cov(i,kts:kte,j))
1744 CALL mynn_tendencies(kts,kte,&
1745 &levflag,grav_settling,&
1746 &delt,&
1747 &dz(i,kts:kte,j),&
1748 &u(i,kts:kte,j),v(i,kts:kte,j),&
1749 &th(i,kts:kte,j),qv(i,kts:kte,j),qc(i,kts:kte,j),&
1750 &p(i,kts:kte,j),exner(i,kts:kte,j),&
1751 &thl,sqv,sqc,sqw,&
1752 &ust(i,j),flt,flq,wspd(i,j),qcg(i,j),&
1753 &tsq(i,kts:kte,j),qsq(i,kts:kte,j),cov(i,kts:kte,j),&
1754 &tcd,qcd,&
1755 &dfm,dfh,dfq,&
1756 &Du(i,kts:kte,j),Dv(i,kts:kte,j),Dth(i,kts:kte,j),&
1757 &Dqv(i,kts:kte,j),Dqc(i,kts:kte,j))
1759 CALL retrieve_exchange_coeffs(kts,kte,&
1760 &dfm,dfh,dfq,dz(i,kts:kte,j),&
1761 &K_m(i,kts:kte,j),K_h(i,kts:kte,j),K_q(i,kts:kte,j))
1764 ENDDO
1765 ENDDO
1767 END SUBROUTINE mynn_bl_driver
1769 SUBROUTINE mynn_bl_init_driver(&
1770 &Du,Dv,Dth,&
1771 &Dqv,Dqc&
1772 &,RESTART,ALLOWED_TO_READ,LEVEL&
1773 &,IDS,IDE,JDS,JDE,KDS,KDE &
1774 &,IMS,IME,JMS,JME,KMS,KME &
1775 &,ITS,ITE,JTS,JTE,KTS,KTE)
1777 LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
1778 INTEGER,INTENT(IN) :: LEVEL
1780 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
1781 & IMS,IME,JMS,JME,KMS,KME, &
1782 & ITS,ITE,JTS,JTE,KTS,KTE
1785 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: &
1786 &Du,Dv,Dth,Dqv,Dqc
1788 INTEGER :: I,J,K,ITF,JTF,KTF
1790 JTF=MIN0(JTE,JDE-1)
1791 KTF=MIN0(KTE,KDE-1)
1792 ITF=MIN0(ITE,IDE-1)
1794 IF(.NOT.RESTART)THEN
1795 DO J=JTS,JTF
1796 DO K=KTS,KTF
1797 DO I=ITS,ITF
1798 Du(i,k,j)=0.
1799 Dv(i,k,j)=0.
1800 Dth(i,k,j)=0.
1801 Dqv(i,k,j)=0.
1802 Dqc(i,k,j)=0.
1803 ENDDO
1804 ENDDO
1805 ENDDO
1806 ENDIF
1808 mynn_level=level
1810 END SUBROUTINE mynn_bl_init_driver
1812 END MODULE module_bl_mynn