module_bl_ysu.F
References to this file elsewhere.
1 !wrf:model_layer:physics
2 !
3 !
4 !
5 !
6 module module_bl_ysu
7 contains
8 !
9 !-------------------------------------------------------------------
10 !
11 subroutine ysu(u3d,v3d,th3d,t3d,qv3d,qc3d,qi3d,p3d,p3di,pi3d, &
12 rublten,rvblten,rthblten, &
13 rqvblten,rqcblten,rqiblten,flag_qi, &
14 cp,g,rovcp,rd,rovg, &
15 dz8w,z,xlv,rv,psfc, &
16 znt,ust,zol,hol,hpbl,psim,psih, &
17 xland,hfx,qfx,tsk,gz1oz0,wspd,br, &
18 dt,dtmin,kpbl2d, &
19 svp1,svp2,svp3,svpt0,ep1,ep2,karman,eomeg,stbolt, &
20 exch_h, &
21 u10,v10, &
22 ids,ide, jds,jde, kds,kde, &
23 ims,ime, jms,jme, kms,kme, &
24 its,ite, jts,jte, kts,kte, &
25 !optional
26 regime )
27 !-------------------------------------------------------------------
28 implicit none
29 !-------------------------------------------------------------------
30 !-- u3d 3d u-velocity interpolated to theta points (m/s)
31 !-- v3d 3d v-velocity interpolated to theta points (m/s)
32 !-- th3d 3d potential temperature (k)
33 !-- t3d temperature (k)
34 !-- qv3d 3d water vapor mixing ratio (kg/kg)
35 !-- qc3d 3d cloud mixing ratio (kg/kg)
36 !-- qi3d 3d ice mixing ratio (kg/kg)
37 ! (note: if P_QI<PARAM_FIRST_SCALAR this should be zero filled)
38 !-- p3d 3d pressure (pa)
39 !-- p3di 3d pressure (pa) at interface level
40 !-- pi3d 3d exner function (dimensionless)
41 !-- rr3d 3d dry air density (kg/m^3)
42 !-- rublten u tendency due to
43 ! pbl parameterization (m/s/s)
44 !-- rvblten v tendency due to
45 ! pbl parameterization (m/s/s)
46 !-- rthblten theta tendency due to
47 ! pbl parameterization (K/s)
48 !-- rqvblten qv tendency due to
49 ! pbl parameterization (kg/kg/s)
50 !-- rqcblten qc tendency due to
51 ! pbl parameterization (kg/kg/s)
52 !-- rqiblten qi tendency due to
53 ! pbl parameterization (kg/kg/s)
54 !-- cp heat capacity at constant pressure for dry air (j/kg/k)
55 !-- g acceleration due to gravity (m/s^2)
56 !-- rovcp r/cp
57 !-- rd gas constant for dry air (j/kg/k)
58 !-- rovg r/g
59 !-- dz8w dz between full levels (m)
60 !-- z height above sea level (m)
61 !-- xlv latent heat of vaporization (j/kg)
62 !-- rv gas constant for water vapor (j/kg/k)
63 !-- psfc pressure at the surface (pa)
64 !-- znt roughness length (m)
65 !-- ust u* in similarity theory (m/s)
66 !-- zol z/l height over monin-obukhov length
67 !-- hol pbl height over monin-obukhov length
68 !-- hpbl pbl height (m)
69 !-- regime flag indicating pbl regime (stable, unstable, etc.)
70 !-- psim similarity stability function for momentum
71 !-- psih similarity stability function for heat
72 !-- xland land mask (1 for land, 2 for water)
73 !-- hfx upward heat flux at the surface (w/m^2)
74 !-- qfx upward moisture flux at the surface (kg/m^2/s)
75 !-- tsk surface temperature (k)
76 !-- gz1oz0 log(z/z0) where z0 is roughness length
77 !-- wspd wind speed at lowest model level (m/s)
78 !-- u10 u-wind speed at 10 m (m/s)
79 !-- v10 v-wind speed at 10 m (m/s)
80 !-- br bulk richardson number in surface layer
81 !-- dt time step (s)
82 !-- dtmin time step (minute)
83 !-- rvovrd r_v divided by r_d (dimensionless)
84 !-- svp1 constant for saturation vapor pressure (kpa)
85 !-- svp2 constant for saturation vapor pressure (dimensionless)
86 !-- svp3 constant for saturation vapor pressure (k)
87 !-- svpt0 constant for saturation vapor pressure (k)
88 !-- ep1 constant for virtual temperature (r_v/r_d - 1) (dimensionless)
89 !-- ep2 constant for specific humidity calculation
90 !-- karman von karman constant
91 !-- eomeg angular velocity of earths rotation (rad/s)
92 !-- stbolt stefan-boltzmann constant (w/m^2/k^4)
93 !-- ids start index for i in domain
94 !-- ide end index for i in domain
95 !-- jds start index for j in domain
96 !-- jde end index for j in domain
97 !-- kds start index for k in domain
98 !-- kde end index for k in domain
99 !-- ims start index for i in memory
100 !-- ime end index for i in memory
101 !-- jms start index for j in memory
102 !-- jme end index for j in memory
103 !-- kms start index for k in memory
104 !-- kme end index for k in memory
105 !-- its start index for i in tile
106 !-- ite end index for i in tile
107 !-- jts start index for j in tile
108 !-- jte end index for j in tile
109 !-- kts start index for k in tile
110 !-- kte end index for k in tile
111 !-------------------------------------------------------------------
112 !
113 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
114 ims,ime, jms,jme, kms,kme, &
115 its,ite, jts,jte, kts,kte
116 !
117 real, intent(in ) :: dt,dtmin,cp,g,rovcp,rovg,rd,xlv,rv
118 !
119 real, intent(in ) :: svp1,svp2,svp3,svpt0
120 real, intent(in ) :: ep1,ep2,karman,eomeg,stbolt
121 !
122 real, dimension( ims:ime, kms:kme, jms:jme ) , &
123 intent(in ) :: qv3d, &
124 qc3d, &
125 qi3d, &
126 p3d, &
127 pi3d, &
128 th3d, &
129 t3d, &
130 dz8w, &
131 z
132 real, dimension( ims:ime, kms:kme, jms:jme ) , &
133 intent(in ) :: p3di
134 !
135 real, dimension( ims:ime, kms:kme, jms:jme ) , &
136 intent(inout) :: rublten, &
137 rvblten, &
138 rthblten, &
139 rqvblten, &
140 rqcblten
141 !
142 real, dimension( ims:ime, kms:kme, jms:jme ) , &
143 intent(inout) :: exch_h
144 real, dimension( ims:ime, jms:jme ) , &
145 intent(in ) :: u10, &
146 v10
147 !
148 real, dimension( ims:ime, jms:jme ) , &
149 intent(in ) :: xland, &
150 hfx, &
151 qfx, &
152 psim, &
153 psih, &
154 gz1oz0, &
155 br, &
156 psfc, &
157 tsk
158 !
159 real, dimension( ims:ime, jms:jme ) , &
160 intent(inout) :: hol, &
161 ust, &
162 hpbl, &
163 znt, &
164 wspd, &
165 zol
166 !
167 real, dimension( ims:ime, kms:kme, jms:jme ) , &
168 intent(in ) :: u3d, &
169 v3d
170 !
171 integer, dimension( ims:ime, jms:jme ) , &
172 intent(out ) :: kpbl2d
173 !
174 logical, intent(in) :: flag_qi
175 !optional
176 real, dimension( ims:ime, jms:jme ) , &
177 optional , &
178 intent(inout) :: regime
179 !
180 real, dimension( ims:ime, kms:kme, jms:jme ) , &
181 optional , &
182 intent(inout) :: rqiblten
183 !
184 !local
185 integer :: i,j,k
186 !
187 do j = jts,jte
188 call ysu2d(J=j,ux=u3d(ims,kms,j),vx=v3d(ims,kms,j) &
189 ,tx=t3d(ims,kms,j) &
190 ,qx=qv3d(ims,kms,j),qcx=qc3d(ims,kms,j) &
191 ,qix=qi3d(ims,kms,j) &
192 ,p2d=p3d(ims,kms,j),p2di=p3di(ims,kms,j) &
193 ,pi2d=pi3d(ims,kms,j) &
194 ,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j) &
195 ,ttnp=rthblten(ims,kms,j),qtnp=rqvblten(ims,kms,j) &
196 ,qctnp=rqcblten(ims,kms,j),qitnp=rqiblten(ims,kms,j) &
197 ,cp=cp,g=g,rovcp=rovcp,rd=rd,rovg=rovg &
198 ,xlv=xlv,rv=rv &
199 ,psfcpa=psfc(ims,j),znt=znt(ims,j),ust=ust(ims,j) &
200 ,zol=zol(ims,j),hol=hol(ims,j),hpbl=hpbl(ims,j) &
201 ,regime=regime(ims,j),psim=psim(ims,j) &
202 ,psih=psih(ims,j),xland=xland(ims,j) &
203 ,hfx=hfx(ims,j),qfx=qfx(ims,j) &
204 ,tsk=tsk(ims,j),gz1oz0=gz1oz0(ims,j) &
205 ,wspd=wspd(ims,j),br=br(ims,j) &
206 ,dt=dt,dtmin=dtmin,kpbl1d=kpbl2d(ims,j) &
207 ,svp1=svp1,svp2=svp2,svp3=svp3,svpt0=svpt0 &
208 ,ep1=ep1,ep2=ep2,karman=karman,eomeg=eomeg &
209 ,stbolt=stbolt &
210 ,exch_hx=exch_h(ims,kms,j) &
211 ,flag_qi=flag_qi &
212 ,u10=u10(ims,j),v10=v10(ims,j) &
213 ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
214 ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
215 ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
216 do k = kts,kte
217 do i = its,ite
218 rthblten(i,k,j) = rthblten(i,k,j)/pi3d(i,k,j)
219 enddo
220 enddo
221 enddo
222 !
223 end subroutine ysu
224 !
225 !-------------------------------------------------------------------
226 !
227 subroutine ysu2d(j,ux,vx,tx,qx,qcx,qix,p2d,p2di,pi2d, &
228 utnp,vtnp,ttnp, &
229 qtnp,qctnp,qitnp, &
230 cp,g,rovcp,rd,rovg, &
231 xlv,rv,psfcpa, &
232 znt,ust,zol,hol,hpbl,psim,psih, &
233 xland,hfx,qfx,tsk,gz1oz0,wspd,br, &
234 dt,dtmin,kpbl1d, &
235 svp1,svp2,svp3,svpt0,ep1,ep2,karman,eomeg,stbolt, &
236 exch_hx,flag_qi, &
237 u10,v10, &
238 ids,ide, jds,jde, kds,kde, &
239 ims,ime, jms,jme, kms,kme, &
240 its,ite, jts,jte, kts,kte, &
241 !optional
242 regime )
243 !-------------------------------------------------------------------
244 implicit none
245 !-------------------------------------------------------------------
246 !
247 ! this code is a revised vertical diffusion package ("ysupbl")
248 ! with a nonlocal turbulent mixing in the pbl after "mrfpbl".
249 ! the ysupbl (hong et al. 2006) is based on the study of noh
250 ! et al.(2003) and accumulated realism of the behavior of the
251 ! troen and mahrt (1986) concept implemented by hong and pan(1996).
252 ! the major ingredient of the ysupbl is the inclusion of an explicit
253 ! treatment of the entrainment processes at the entrainment layer.
254 ! this routine uses an implicit approach for vertical flux
255 ! divergence and does not require "miter" timesteps.
256 ! it includes vertical diffusion in the stable atmosphere
257 ! and moist vertical diffusion in clouds.
258 !
259 ! mrfpbl:
260 ! coded by song-you hong (ncep), implemented by jimy dudhia (ncar)
261 ! fall 1996
262 !
263 ! ysupbl:
264 ! coded by song-you hong (yonsei university) and implemented by
265 ! song-you hong (yonsei university) and jimy dudhia (ncar)
266 ! summer 2002
267 !
268 ! references:
269 !
270 ! hong, noh, and dudhia (2006), mon. wea. rev.
271 ! hong and pan (1996), mon. wea. rev.
272 ! noh, chun, hong, and raasch (2003), boundary layer met.
273 ! troen and mahrt (1986), boundary layer met.
274 !
275 !-------------------------------------------------------------------
276 !
277 integer,parameter :: ncloud = 3
278 real,parameter :: xkzmin = 0.01,xkzmax = 1000.,rimin = -100.
279 real,parameter :: rlam = 30.,prmin = 0.25,prmax = 4.
280 real,parameter :: brcr_ub = 0.0,brcr_sb = 0.25,cori = 1.e-4
281 real,parameter :: afac = 6.8,bfac = 6.8,pfac = 2.0
282 real,parameter :: phifac = 8.,sfcfrac = 0.1
283 real,parameter :: d1 = 0.02, d2 = 0.05, d3 = 0.001
284 real,parameter :: h1 = 0.33333335, h2 = 0.6666667
285 real,parameter :: ckz = 0.001,zfmin = 1.e-8,aphi5 = 5.,aphi16 = 16.
286 real,parameter :: tmin=1.e-2
287 real,parameter :: gamcrt = 3.,gamcrq = 2.e-3
288 real,parameter :: xka = 2.4e-5
289 !
290 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
291 ims,ime, jms,jme, kms,kme, &
292 its,ite, jts,jte, kts,kte, j
293 !
294 real, intent(in ) :: dt,dtmin,cp,g,rovcp,rovg,rd,xlv,rv
295 !
296 real, intent(in ) :: svp1,svp2,svp3,svpt0
297 real, intent(in ) :: ep1,ep2,karman,eomeg,stbolt
298 !
299 real, dimension( ims:ime, kms:kme ) , &
300 intent(in ) :: tx, &
301 qx, &
302 qcx, &
303 qix, &
304 p2d, &
305 pi2d
306 real, dimension( ims:ime, kms:kme ) , &
307 intent(in ) :: p2di
308 !
309 real, dimension( ims:ime, kms:kme ) , &
310 intent(inout) :: utnp, &
311 vtnp, &
312 ttnp, &
313 qtnp, &
314 qctnp, &
315 qitnp
316 !
317 real, dimension( ims:ime ) , &
318 intent(inout) :: hol, &
319 ust, &
320 hpbl, &
321 znt
322 real, dimension( ims:ime ) , &
323 intent(in ) :: xland, &
324 hfx, &
325 qfx
326 !
327 real, dimension( ims:ime ), intent(inout) :: wspd
328 real, dimension( ims:ime ), intent(in ) :: br
329 !
330 real, dimension( ims:ime ), intent(in ) :: psim, &
331 psih
332 real, dimension( ims:ime ), intent(in ) :: gz1oz0
333 !
334 real, dimension( ims:ime ), intent(in ) :: psfcpa
335 real, dimension( ims:ime ), intent(in ) :: tsk
336 real, dimension( ims:ime ), intent(inout) :: zol
337 integer, dimension( ims:ime ), intent(out ) :: kpbl1d
338 !
339 logical, intent(in) :: flag_qi
340 !
341 !
342 real, dimension( ims:ime, kms:kme ) , &
343 intent(in ) :: ux, &
344 vx
345 !optional
346 real, dimension( ims:ime ) , &
347 optional , &
348 intent(inout) :: regime
349 !
350 ! local vars
351 !
352 real, dimension( its:ite, kts:kte+1 ) :: zq
353 !
354 real, dimension( its:ite, kts:kte ) :: &
355 thx,thvx, &
356 del, &
357 dza, &
358 dzq, &
359 za, &
360 tvx, &
361 uxs,vxs, &
362 thxs,qxs, &
363 qcxs,qixs
364 !
365 real, dimension( its:ite ) :: qixsv,rhox, &
366 govrth, &
367 thxsv, &
368 uxsv,vxsv, &
369 qxsv,qcxsv, &
370 qgh,tgdsa,ps
371 !
372 real, dimension( its:ite ) :: &
373 zl1,thermal, &
374 wscale,hgamt, &
375 hgamq,brdn, &
376 brup,phim, &
377 phih,cpm, &
378 dusfc,dvsfc, &
379 dtsfc,dqsfc, &
380 thgb,prpbl, &
381 wspd1
382 !
383 real, dimension( its:ite, kts:kte ) :: xkzm,xkzh, &
384 f1,f2, &
385 r1,r2, &
386 ad,au, &
387 cu, &
388 al, &
389 zfac
390 !
391 !jdf added exch_hx
392 real, dimension( ims:ime, kms:kme ) , &
393 intent(inout) :: exch_hx
394 !
395 real, dimension( ims:ime ) , &
396 intent(in ) :: u10, &
397 v10
398 real, dimension( its:ite ) :: &
399 brcr_sbro
400 !
401 real, dimension( its:ite, kts:kte, ncloud) :: r3,f3
402 !
403 logical, dimension( its:ite ) :: pblflg, &
404 sfcflg, &
405 stable
406 !
407 integer :: n,i,k,l,nzol,imvdif,ic
408 integer :: klpbl
409 !
410 integer, dimension( its:ite ) :: kpbl
411 !
412 real :: zoln,x,y,tvcon,e1,dtstep,brcr
413 real :: zl,tskv,dthvdz,dthvm,vconv,rzol
414 real :: dtthx,psix,dtg,psiq,ustm
415 real :: dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum
416 real :: xkzo,ss,ri,qmean,tmean,alph,chi,zk,rl2,dk,sri
417 real :: brint,dtodsd,rdz,dsdzt,dsdzq,dsdz2,ttend,qtend
418 real :: utend,vtend,qctend,qitend,tgc,dtodsu,govrthv
419 !
420 real, dimension( its:ite, kts:kte ) :: wscalek, &
421 xkzml,xkzhl, &
422 zfacent,entfac
423 real, dimension( its:ite ) :: ust3, &
424 wstar3,wstar, &
425 hgamu,hgamv, &
426 wm2, we, &
427 bfxpbl, &
428 hfxpbl,qfxpbl, &
429 ufxpbl,vfxpbl, &
430 delta,dthvx
431 real :: prnumfac,bfx0,hfx0,qfx0,delb,dux,dvx, &
432 dsdzu,dsdzv,wm3,dthx,dqx,wspd10,ross,tem1,dsig
433 !
434 !----------------------------------------------------------------------
435 !
436 klpbl = kte
437 !
438 !-- imvdif imvdif = 1 for moist adiabat vertical diffusion
439 imvdif = 1
440 !
441 !----convert ground temperature to potential temperature:
442 !
443 do i = its,ite
444 tgdsa(i) = tsk(i)
445 ps(i) = psfcpa(i)/1000. ! ps psfc cb
446 thgb(i) = tsk(i)*(100./ps(i))**rovcp
447 enddo
448 !
449 do k = kts,kte
450 do i = its,ite
451 thx(i,k) = tx(i,k)/pi2d(i,k)
452 enddo
453 enddo
454 !
455 do i = its,ite
456 qgh(i) = 0.
457 cpm(i) = cp
458 enddo
459 !
460 do k = kts,kte
461 do i = its,ite
462 tvcon = (1.+ep1*qx(i,k))
463 thvx(i,k) = thx(i,k)*tvcon
464 enddo
465 enddo
466 !
467 do i = its,ite
468 e1 = svp1*exp(svp2*(tgdsa(i)-svpt0)/(tgdsa(i)-svp3))
469 qgh(i) = ep2*e1/(ps(i)-e1)
470 cpm(i) = cp*(1.+0.8*qx(i,1))
471 enddo
472 !
473 !-----compute the height of full- and half-sigma levels above ground
474 ! level, and the layer thicknesses.
475 !
476 do i = its,ite
477 zq(i,1) = 0.
478 rhox(i) = ps(i)*1000./(rd*tx(i,1))
479 enddo
480 !
481 do k = kts,kte
482 do i = its,ite
483 tvx(i,k) = tx(i,k)*(1.+ep1*qx(i,k))
484 del(i,k) = p2di(i,k)-p2di(i,k+1)
485 dzq(i,k) = del(i,k)/p2d(i,k)*tvx(i,k)*rovg
486 zq(i,k+1) = dzq(i,k)+zq(i,k)
487 enddo
488 enddo
489 !
490 do k = kts,kte
491 do i = its,ite
492 za(i,1) = 0.5*(zq(i,1)+zq(i,2))
493 dza(i,1) = za(i,1)
494 enddo
495 enddo
496 !
497 do k = kts,kte-1
498 do i = its,ite
499 dza(i,k+1) = (p2d(i,k)-p2d(i,k+1))/p2di(i,k+1)* &
500 0.5*(tvx(i,k)+tvx(i,k+1))*rovg
501 za(i,k+1) = za(i,k)+dza(i,k+1)
502 enddo
503 enddo
504 !
505 do i = its,ite
506 govrth(i) = g/thx(i,1)
507 enddo
508 !
509 !-----initialize vertical tendencies and
510 !
511 do i = its,ite
512 do k = kts,kte
513 utnp(i,k) = 0.
514 vtnp(i,k) = 0.
515 ttnp(i,k) = 0.
516 enddo
517 enddo
518 !
519 do k = kts,kte
520 do i = its,ite
521 qtnp(i,k) = 0.
522 enddo
523 enddo
524 !
525 do k = kts,kte
526 do i = its,ite
527 qctnp(i,k) = 0.
528 qitnp(i,k) = 0.
529 enddo
530 enddo
531 !
532 do i = its,ite
533 wspd1(i) = sqrt(ux(i,1)*ux(i,1)+vx(i,1)*vx(i,1))+1.e-9
534 enddo
535 !
536 !---- compute vertical diffusion
537 !
538 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
539 ! compute preliminary variables
540 !
541 dtstep = dt
542 dt2 = 2.*dtstep
543 rdt = 1./dt2
544 !
545 do i = its,ite
546 bfxpbl(i) = 0.0
547 hfxpbl(i) = 0.0
548 qfxpbl(i) = 0.0
549 ufxpbl(i) = 0.0
550 vfxpbl(i) = 0.0
551 hgamu(i) = 0.0
552 hgamv(i) = 0.0
553 delta(i) = 0.0
554 enddo
555 !
556 do k = kts,klpbl
557 do i = its,ite
558 wscalek(i,k) = 0.0
559 enddo
560 enddo
561 !
562 do k = kts,klpbl
563 do i = its,ite
564 zfac(i,k) = 0.0
565 enddo
566 enddo
567 !
568 do i = its,ite
569 hgamt(i) = 0.
570 hgamq(i) = 0.
571 wscale(i) = 0.
572 kpbl(i) = 1
573 hpbl(i) = zq(i,1)
574 zl1(i) = za(i,1)
575 thermal(i)= thvx(i,1)
576 pblflg(i) = .true.
577 sfcflg(i) = .true.
578 if(br(i).gt.0.0) sfcflg(i) = .false.
579 enddo
580 !
581 ! compute the first guess of pbl height
582 !
583 do i = its,ite
584 stable(i) = .false.
585 brup(i) = br(i)
586 enddo
587 !
588 brcr = brcr_ub
589 do k = 2,klpbl
590 do i = its,ite
591 if(.not.stable(i))then
592 brdn(i) = brup(i)
593 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
594 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
595 kpbl(i) = k
596 stable(i) = brup(i).gt.brcr
597 endif
598 enddo
599 enddo
600 !
601 do i = its,ite
602 k = kpbl(i)
603 if(brdn(i).ge.brcr)then
604 brint = 0.
605 elseif(brup(i).le.brcr)then
606 brint = 1.
607 else
608 brint = (brcr-brdn(i))/(brup(i)-brdn(i))
609 endif
610 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
611 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
612 if(kpbl(i).le.1) pblflg(i) = .false.
613 enddo
614 !
615 do i = its,ite
616 fm = gz1oz0(i)-psim(i)
617 fh = gz1oz0(i)-psih(i)
618 hol(i) = max(br(i)*fm*fm/fh,rimin)
619 if(sfcflg(i))then
620 hol(i) = min(hol(i),-zfmin)
621 else
622 hol(i) = max(hol(i),zfmin)
623 endif
624 hol1 = hol(i)*hpbl(i)/zl1(i)*sfcfrac
625 hol(i) = -hol(i)*hpbl(i)/zl1(i)
626 if(sfcflg(i))then
627 phim(i) = (1.-aphi16*hol1)**(-1./4.)
628 phih(i) = (1.-aphi16*hol1)**(-1./2.)
629 bfx0 = max(hfx(i)/rhox(i)/cpm(i)+ep1*thx(i,1)*qfx(i)/rhox(i),0.)
630 hfx0 = max(hfx(i)/rhox(i)/cpm(i),0.)
631 qfx0 = max(ep1*thx(i,1)*qfx(i)/rhox(i),0.)
632 wstar3(i) = (govrth(i)*bfx0*hpbl(i))
633 wstar(i) = (wstar3(i))**h1
634 else
635 phim(i) = (1.+aphi5*hol1)
636 phih(i) = phim(i)
637 wstar(i) = 0.
638 wstar3(i) = 0.
639 endif
640 ust3(i) = ust(i)**3.
641 wscale(i) = (ust3(i)+phifac*karman*wstar3(i)*0.5)**h1
642 wscale(i) = min(wscale(i),ust(i)*aphi16)
643 wscale(i) = max(wscale(i),ust(i)/aphi5)
644 enddo
645 !
646 ! compute the surface variables for pbl height estimation
647 ! under unstable conditions
648 !
649 do i = its,ite
650 if(sfcflg(i))then
651 gamfac = bfac/rhox(i)/wscale(i)
652 hgamt(i) = min(gamfac*hfx(i)/cpm(i),gamcrt)
653 hgamq(i) = min(gamfac*qfx(i),gamcrq)
654 vpert = (hgamt(i)+ep1*thx(i,1)*hgamq(i))/bfac*afac
655 thermal(i) = thermal(i)+max(vpert,0.)
656 hgamt(i) = max(hgamt(i),0.0)
657 hgamq(i) = max(hgamq(i),0.0)
658 brint = -15.9*ust(i)*ust(i)/wspd(i)*wstar3(i)/(wscale(i)**4.)
659 hgamu(i) = brint*ux(i,1)
660 hgamv(i) = brint*vx(i,1)
661 else
662 pblflg(i) = .false.
663 endif
664 enddo
665 !
666 ! enhance the pbl height by considering the thermal
667 !
668 do i = its,ite
669 if(pblflg(i))then
670 kpbl(i) = 1
671 hpbl(i) = zq(i,1)
672 endif
673 enddo
674 !
675 do i = its,ite
676 if(pblflg(i))then
677 stable(i) = .false.
678 brup(i) = br(i)
679 endif
680 enddo
681 !
682 brcr = brcr_ub
683 do k = 2,klpbl
684 do i = its,ite
685 if(.not.stable(i).and.pblflg(i))then
686 brdn(i) = brup(i)
687 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
688 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
689 kpbl(i) = k
690 stable(i) = brup(i).gt.brcr
691 endif
692 enddo
693 enddo
694 !
695 do i = its,ite
696 if(pblflg(i)) then
697 k = kpbl(i)
698 if(brdn(i).ge.brcr)then
699 brint = 0.
700 elseif(brup(i).le.brcr)then
701 brint = 1.
702 else
703 brint = (brcr-brdn(i))/(brup(i)-brdn(i))
704 endif
705 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
706 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
707 if(kpbl(i).le.1) pblflg(i) = .false.
708 endif
709 enddo
710 !
711 ! stable boundary layer
712 !
713 do i = its,ite
714 if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then
715 brup(i) = br(i)
716 stable(i) = .false.
717 else
718 stable(i) = .true.
719 endif
720 enddo
721 !
722 do i = its,ite
723 if((.not.stable(i)).and.((xland(i)-1.5).ge.0))then
724 wspd10 = u10(i)*u10(i) + v10(i)*v10(i)
725 wspd10 = sqrt(wspd10)
726 ross = wspd10 / (cori*znt(i))
727 brcr_sbro(i) = min(0.16*(1.e-7*ross)**(-0.18),.3)
728 endif
729 enddo
730 !
731 do k = 2,klpbl
732 do i = its,ite
733 if((xland(i)-1.5).ge.0)then
734 brcr = brcr_sbro(i)
735 else
736 brcr = brcr_sb
737 endif
738 if(.not.stable(i))then
739 brdn(i) = brup(i)
740 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
741 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
742 kpbl(i) = k
743 stable(i) = brup(i).gt.brcr
744 endif
745 enddo
746 enddo
747 !
748 do i = its,ite
749 if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then
750 if((xland(i)-1.5).ge.0)then
751 brcr = brcr_sbro(i)
752 else
753 brcr = brcr_sb
754 endif
755 k = kpbl(i)
756 if(brdn(i).ge.brcr)then
757 brint = 0.
758 elseif(brup(i).le.brcr)then
759 brint = 1.
760 else
761 brint = (brcr-brdn(i))/(brup(i)-brdn(i))
762 endif
763 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
764 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
765 if(kpbl(i).le.1) pblflg(i) = .false.
766 endif
767 enddo
768 !
769 ! estimate the entrainment parameters
770 !
771 do i = its,ite
772 if(pblflg(i)) then
773 k = kpbl(i) - 1
774 prpbl(i) = 1.0
775 wm3 = wstar3(i) + 5. * ust3(i)
776 wm2(i) = wm3**h2
777 bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i)
778 dthvx(i) = max(thvx(i,k+1)-thvx(i,k),tmin)
779 dthx = max(thx(i,k+1)-thx(i,k),tmin)
780 dqx = min(qx(i,k+1)-qx(i,k),0.0)
781 we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i)))
782 hfxpbl(i) = we(i)*dthx
783 qfxpbl(i) = we(i)*dqx
784 !
785 dux = ux(i,k+1)-ux(i,k)
786 dvx = vx(i,k+1)-vx(i,k)
787 if(dux.gt.tmin) then
788 ufxpbl(i) = max(prpbl(i)*we(i)*dux,-ust(i)*ust(i))
789 elseif(dux.lt.-tmin) then
790 ufxpbl(i) = min(prpbl(i)*we(i)*dux,ust(i)*ust(i))
791 else
792 ufxpbl(i) = 0.0
793 endif
794 if(dvx.gt.tmin) then
795 vfxpbl(i) = max(prpbl(i)*we(i)*dvx,-ust(i)*ust(i))
796 elseif(dvx.lt.-tmin) then
797 vfxpbl(i) = min(prpbl(i)*we(i)*dvx,ust(i)*ust(i))
798 else
799 vfxpbl(i) = 0.0
800 endif
801 delb = govrth(i)*d3*hpbl(i)
802 delta(i) = min(d1*hpbl(i) + d2*wm2(i)/delb,100.)
803 endif
804 enddo
805 !
806 do k = kts,klpbl
807 do i = its,ite
808 if(pblflg(i).and.k.ge.kpbl(i))then
809 entfac(i,k) = ((zq(i,k+1)-hpbl(i))/delta(i))**2.
810 else
811 entfac(i,k) = 1.e30
812 endif
813 enddo
814 enddo
815 !
816 ! compute diffusion coefficients below pbl
817 !
818 do k = kts,klpbl
819 do i = its,ite
820 if(k.lt.kpbl(i)) then
821 zfac(i,k) = min(max((1.-(zq(i,k+1)-zl1(i))/(hpbl(i)-zl1(i))),zfmin),1.)
822 xkzo = ckz*dza(i,k+1)
823 zfacent(i,k) = (1.-zfac(i,k))**3.
824 prnumfac = -3.*(max(zq(i,k+1)-sfcfrac*hpbl(i),0.))**2./hpbl(i)**2.
825 prnum = (phih(i)/phim(i)+bfac*karman*sfcfrac)
826 prnum = 1. + (prnum-1.)*exp(prnumfac)
827 prnum = min(prnum,prmax)
828 prnum = max(prnum,prmin)
829 wscalek(i,k) = (ust3(i)+phifac*karman*wstar3(i)*(1.-zfac(i,k)))**h1
830 xkzm(i,k) = xkzo+wscalek(i,k)*karman*zq(i,k+1)*zfac(i,k)**pfac
831 xkzh(i,k) = xkzm(i,k)/prnum
832 xkzm(i,k) = min(xkzm(i,k),xkzmax)
833 xkzm(i,k) = max(xkzm(i,k),xkzmin)
834 xkzh(i,k) = min(xkzh(i,k),xkzmax)
835 xkzh(i,k) = max(xkzh(i,k),xkzmin)
836 endif
837 enddo
838 enddo
839 !
840 ! compute diffusion coefficients over pbl (free atmosphere)
841 !
842 do k = kts,kte-1
843 do i = its,ite
844 xkzo = ckz*dza(i,k+1)
845 if(k.ge.kpbl(i)) then
846 ss = ((ux(i,k+1)-ux(i,k))*(ux(i,k+1)-ux(i,k)) &
847 +(vx(i,k+1)-vx(i,k))*(vx(i,k+1)-vx(i,k))) &
848 /(dza(i,k+1)*dza(i,k+1))+1.e-9
849 govrthv = g/(0.5*(thvx(i,k+1)+thvx(i,k)))
850 ri = govrthv*(thvx(i,k+1)-thvx(i,k))/(ss*dza(i,k+1))
851 if(imvdif.eq.1)then
852 if((qcx(i,k)+qix(i,k)).gt.0.01e-3.and.(qcx(i,k+1)+ &
853 qix(i,k+1)).gt.0.01e-3)then
854 ! in cloud
855 qmean = 0.5*(qx(i,k)+qx(i,k+1))
856 tmean = 0.5*(tx(i,k)+tx(i,k+1))
857 alph = xlv*qmean/rd/tmean
858 chi = xlv*xlv*qmean/cp/rv/tmean/tmean
859 ri = (1.+alph)*(ri-g*g/ss/tmean/cp*((chi-alph)/(1.+chi)))
860 endif
861 endif
862 zk = karman*zq(i,k+1)
863 rl2 = (zk*rlam/(rlam+zk))**2
864 dk = rl2*sqrt(ss)
865 if(ri.lt.0.)then
866 ! unstable regime
867 sri = sqrt(-ri)
868 xkzm(i,k) = xkzo+dk*(1+8.*(-ri)/(1+1.746*sri))
869 xkzh(i,k) = xkzo+dk*(1+8.*(-ri)/(1+1.286*sri))
870 else
871 ! stable regime
872 xkzh(i,k) = xkzo+dk/(1+5.*ri)**2
873 prnum = 1.0+2.1*ri
874 prnum = min(prnum,prmax)
875 xkzm(i,k) = (xkzh(i,k)-xkzo)*prnum+xkzo
876 endif
877 !
878 xkzm(i,k) = min(xkzm(i,k),xkzmax)
879 xkzm(i,k) = max(xkzm(i,k),xkzmin)
880 xkzh(i,k) = min(xkzh(i,k),xkzmax)
881 xkzh(i,k) = max(xkzh(i,k),xkzmin)
882 xkzml(i,k) = xkzm(i,k)
883 xkzhl(i,k) = xkzh(i,k)
884 endif
885 enddo
886 enddo
887 !
888 ! compute tridiagonal matrix elements for heat and moisture, and clouds
889 !
890 do k = kts,kte
891 do i = its,ite
892 au(i,k) = 0.
893 al(i,k) = 0.
894 ad(i,k) = 0.
895 f1(i,k) = 0.
896 enddo
897 enddo
898 !
899 do ic = 1,ncloud
900 do i = its,ite
901 do k = kts,kte
902 f3(i,k,ic) = 0.
903 enddo
904 enddo
905 enddo
906 !
907 do i = its,ite
908 ad(i,1) = 1.
909 f1(i,1) = thx(i,1)-300.+hfx(i)/(rhox(i)*cpm(i))/zq(i,2)*dt2
910 f3(i,1,1) = qx(i,1)+qfx(i)/(rhox(i))/zq(i,2)*dt2
911 enddo
912 !
913 if(ncloud.ge.2) then
914 do ic = 2,ncloud
915 do i = its,ite
916 if(ic.eq.2) then
917 f3(i,1,ic) = qcx(i,1)
918 elseif(ic.eq.3) then
919 f3(i,1,ic) = qix(i,1)
920 endif
921 enddo
922 enddo
923 endif
924 !
925 do k = kts,kte-1
926 do i = its,ite
927 dtodsd = dt2/del(i,k)
928 dtodsu = dt2/del(i,k+1)
929 dsig = p2d(i,k)-p2d(i,k+1)
930 rdz = 1./dza(i,k+1)
931 tem1 = dsig*xkzh(i,k)*rdz
932 if(pblflg(i).and.k.lt.kpbl(i)) then
933 dsdzt = tem1*(-hgamt(i)/hpbl(i)-hfxpbl(i)*zfacent(i,k)/xkzh(i,k))
934 dsdzq = tem1*(-qfxpbl(i)*zfacent(i,k)/xkzh(i,k))
935 f1(i,k) = f1(i,k)+dtodsd*dsdzt
936 f1(i,k+1) = thx(i,k+1)-300.-dtodsu*dsdzt
937 f3(i,k,1) = f3(i,k,1)+dtodsd*dsdzq
938 f3(i,k+1,1) = qx(i,k+1)-dtodsu*dsdzq
939 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
940 xkzh(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
941 xkzh(i,k) = sqrt(xkzh(i,k)*xkzhl(i,k))
942 xkzh(i,k) = min(xkzh(i,k),xkzmax)
943 xkzh(i,k) = max(xkzh(i,k),xkzmin)
944 f1(i,k+1) = thx(i,k+1)-300.
945 f3(i,k+1,1) = qx(i,k+1)
946 else
947 f1(i,k+1) = thx(i,k+1)-300.
948 f3(i,k+1,1) = qx(i,k+1)
949 endif
950 dsdz2 = tem1*rdz
951 au(i,k) = -dtodsd*dsdz2
952 al(i,k) = -dtodsu*dsdz2
953 ad(i,k) = ad(i,k)-au(i,k)
954 ad(i,k+1) = 1.-al(i,k)
955 exch_hx(i,k) = xkzh(i,k)
956 enddo
957 enddo
958 !
959 if(ncloud.ge.2) then
960 do ic = 2,ncloud
961 do k = kts,kte-1
962 do i = its,ite
963 if(ic.eq.2) then
964 f3(i,k+1,ic) = qcx(i,k+1)
965 elseif(ic.eq.3) then
966 f3(i,k+1,ic) = qix(i,k+1)
967 endif
968 enddo
969 enddo
970 enddo
971 endif
972 !
973 ! copies here to avoid duplicate input args for tridin
974 !
975 do k = kts,kte
976 do i = its,ite
977 cu(i,k) = au(i,k)
978 r1(i,k) = f1(i,k)
979 enddo
980 enddo
981 !
982 do ic = 1,ncloud
983 do k = kts,kte
984 do i = its,ite
985 r3(i,k,ic) = f3(i,k,ic)
986 enddo
987 enddo
988 enddo
989 !
990 ! solve tridiagonal problem for heat and moisture, and clouds
991 !
992 call tridin(al,ad,cu,r1,r3,au,f1,f3,its,ite,kts,kte,ncloud)
993 !
994 ! recover tendencies of heat and moisture
995 !
996 do k = kte,kts,-1
997 do i = its,ite
998 ttend = (f1(i,k)-thx(i,k)+300.)*rdt*pi2d(i,k)
999 qtend = (f3(i,k,1)-qx(i,k))*rdt
1000 ttnp(i,k) = ttnp(i,k)+ttend
1001 qtnp(i,k) = qtnp(i,k)+qtend
1002 enddo
1003 enddo
1004 !
1005 if(ncloud.ge.2) then
1006 do ic = 2,ncloud
1007 do k = kte,kts,-1
1008 do i = its,ite
1009 if(ic.eq.2) then
1010 qctend = (f3(i,k,ic)-qcx(i,k))*rdt
1011 qctnp(i,k) = qctnp(i,k)+qctend
1012 elseif(ic.eq.3) then
1013 qitend = (f3(i,k,ic)-qix(i,k))*rdt
1014 if(flag_qi)qitnp(i,k) = qitnp(i,k)+qitend
1015 endif
1016 enddo
1017 enddo
1018 enddo
1019 endif
1020 !
1021 ! compute tridiagonal matrix elements for momentum
1022 !
1023 do i = its,ite
1024 do k = kts,kte
1025 au(i,k) = 0.
1026 al(i,k) = 0.
1027 ad(i,k) = 0.
1028 f1(i,k) = 0.
1029 f2(i,k) = 0.
1030 enddo
1031 enddo
1032 !
1033 do i = its,ite
1034 ad(i,1) = 1.
1035 f1(i,1) = ux(i,1)-ux(i,1)/wspd1(i)*ust(i)*ust(i)/zq(i,2)*dt2 &
1036 *(wspd1(i)/wspd(i))**2
1037 f2(i,1) = vx(i,1)-vx(i,1)/wspd1(i)*ust(i)*ust(i)/zq(i,2)*dt2 &
1038 *(wspd1(i)/wspd(i))**2
1039 enddo
1040 !
1041 do k = kts,kte-1
1042 do i = its,ite
1043 dtodsd = dt2/del(i,k)
1044 dtodsu = dt2/del(i,k+1)
1045 dsig = p2d(i,k)-p2d(i,k+1)
1046 rdz = 1./dza(i,k+1)
1047 tem1 = dsig*xkzm(i,k)*rdz
1048 if(pblflg(i).and.k.lt.kpbl(i))then
1049 dsdzu = tem1*(-hgamu(i)/hpbl(i)-ufxpbl(i)*zfacent(i,k)/xkzm(i,k))
1050 dsdzv = tem1*(-hgamv(i)/hpbl(i)-vfxpbl(i)*zfacent(i,k)/xkzm(i,k))
1051 f1(i,k) = f1(i,k)+dtodsd*dsdzu
1052 f1(i,k+1) = ux(i,k+1)-dtodsu*dsdzu
1053 f2(i,k) = f2(i,k)+dtodsd*dsdzv
1054 f2(i,k+1) = vx(i,k+1)-dtodsu*dsdzv
1055 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
1056 xkzm(i,k) = prpbl(i)*xkzh(i,k)
1057 xkzm(i,k) = sqrt(xkzm(i,k)*xkzml(i,k))
1058 xkzm(i,k) = min(xkzm(i,k),xkzmax)
1059 xkzm(i,k) = max(xkzm(i,k),xkzmin)
1060 f1(i,k+1) = ux(i,k+1)
1061 f2(i,k+1) = vx(i,k+1)
1062 else
1063 f1(i,k+1) = ux(i,k+1)
1064 f2(i,k+1) = vx(i,k+1)
1065 endif
1066 dsdz2 = tem1*rdz
1067 au(i,k) = -dtodsd*dsdz2
1068 al(i,k) = -dtodsu*dsdz2
1069 ad(i,k) = ad(i,k)-au(i,k)
1070 ad(i,k+1) = 1.-al(i,k)
1071 enddo
1072 enddo
1073 !
1074 ! copies here to avoid duplicate input args for tridin
1075 !
1076 do k = kts,kte
1077 do i = its,ite
1078 cu(i,k) = au(i,k)
1079 r1(i,k) = f1(i,k)
1080 r2(i,k) = f2(i,k)
1081 enddo
1082 enddo
1083 !
1084 ! solve tridiagonal problem for momentum
1085 !
1086 call tridin(al,ad,cu,r1,r2,au,f1,f2,its,ite,kts,kte,1)
1087 !
1088 ! recover tendencies of momentum
1089 !
1090 do k = kte,kts,-1
1091 do i = its,ite
1092 utend = (f1(i,k)-ux(i,k))*rdt
1093 vtend = (f2(i,k)-vx(i,k))*rdt
1094 utnp(i,k) = utnp(i,k)+utend
1095 vtnp(i,k) = vtnp(i,k)+vtend
1096 enddo
1097 enddo
1098 !
1099 !---- end of vertical diffusion
1100 !
1101 do i = its,ite
1102 kpbl1d(i) = kpbl(i)
1103 enddo
1104 !
1105 end subroutine ysu2d
1106 !
1107 subroutine tridin(cl,cm,cu,r1,r2,au,f1,f2,its,ite,kts,kte,nt)
1108 !----------------------------------------------------------------
1109 implicit none
1110 !----------------------------------------------------------------
1111 !
1112 integer, intent(in ) :: its,ite, kts,kte, nt
1113 !
1114 real, dimension( its:ite, kts+1:kte+1 ) , &
1115 intent(in ) :: cl
1116 !
1117 real, dimension( its:ite, kts:kte ) , &
1118 intent(in ) :: cm, &
1119 r1
1120 real, dimension( its:ite, kts:kte,nt ) , &
1121 intent(in ) :: r2
1122 !
1123 real, dimension( its:ite, kts:kte ) , &
1124 intent(inout) :: au, &
1125 cu, &
1126 f1
1127 real, dimension( its:ite, kts:kte,nt ) , &
1128 intent(inout) :: f2
1129 !
1130 real :: fk
1131 integer :: i,k,l,n,it
1132 !
1133 !----------------------------------------------------------------
1134 !
1135 l = ite
1136 n = kte
1137 !
1138 do i = its,l
1139 fk = 1./cm(i,1)
1140 au(i,1) = fk*cu(i,1)
1141 f1(i,1) = fk*r1(i,1)
1142 enddo
1143 do it = 1,nt
1144 do i = its,l
1145 fk = 1./cm(i,1)
1146 f2(i,1,it) = fk*r2(i,1,it)
1147 enddo
1148 enddo
1149 do k = kts+1,n-1
1150 do i = its,l
1151 fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1152 au(i,k) = fk*cu(i,k)
1153 f1(i,k) = fk*(r1(i,k)-cl(i,k)*f1(i,k-1))
1154 enddo
1155 enddo
1156 do it = 1,nt
1157 do k = kts+1,n-1
1158 do i = its,l
1159 fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1160 f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it))
1161 enddo
1162 enddo
1163 enddo
1164 do i = its,l
1165 fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1166 f1(i,n) = fk*(r1(i,n)-cl(i,n)*f1(i,n-1))
1167 enddo
1168 do it = 1,nt
1169 do i = its,l
1170 fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1171 f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it))
1172 enddo
1173 enddo
1174 do k = n-1,kts,-1
1175 do i = its,l
1176 f1(i,k) = f1(i,k)-au(i,k)*f1(i,k+1)
1177 enddo
1178 enddo
1179 do it = 1,nt
1180 do k = n-1,kts,-1
1181 do i = its,l
1182 f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it)
1183 enddo
1184 enddo
1185 enddo
1186 !
1187 end subroutine tridin
1188 !
1189 subroutine ysuinit(rublten,rvblten,rthblten,rqvblten, &
1190 rqcblten,rqiblten,p_qi,p_first_scalar, &
1191 restart, allowed_to_read, &
1192 ids, ide, jds, jde, kds, kde, &
1193 ims, ime, jms, jme, kms, kme, &
1194 its, ite, jts, jte, kts, kte )
1195 !-------------------------------------------------------------------
1196 implicit none
1197 !-------------------------------------------------------------------
1198 !
1199 logical , intent(in) :: restart, allowed_to_read
1200 integer , intent(in) :: ids, ide, jds, jde, kds, kde, &
1201 ims, ime, jms, jme, kms, kme, &
1202 its, ite, jts, jte, kts, kte
1203 integer , intent(in) :: p_qi,p_first_scalar
1204 real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) :: &
1205 rublten, &
1206 rvblten, &
1207 rthblten, &
1208 rqvblten, &
1209 rqcblten, &
1210 rqiblten
1211 integer :: i, j, k, itf, jtf, ktf
1212 !
1213 jtf = min0(jte,jde-1)
1214 ktf = min0(kte,kde-1)
1215 itf = min0(ite,ide-1)
1216 !
1217 if(.not.restart)then
1218 do j = jts,jtf
1219 do k = kts,ktf
1220 do i = its,itf
1221 rublten(i,k,j) = 0.
1222 rvblten(i,k,j) = 0.
1223 rthblten(i,k,j) = 0.
1224 rqvblten(i,k,j) = 0.
1225 rqcblten(i,k,j) = 0.
1226 enddo
1227 enddo
1228 enddo
1229 endif
1230 !
1231 if (p_qi .ge. p_first_scalar .and. .not.restart) then
1232 do j = jts,jtf
1233 do k = kts,ktf
1234 do i = its,itf
1235 rqiblten(i,k,j) = 0.
1236 enddo
1237 enddo
1238 enddo
1239 endif
1240 !
1241 end subroutine ysuinit
1242 !-------------------------------------------------------------------
1243 end module module_bl_ysu