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