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