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