standard WRF version 3.0.1.1
[wrffire.git] / wrfv2_fire / phys / module_bl_ysu.F
blob324397c13d315b62fc5f575887bb5bd65931b40d
1 !wrf:model_layer:physics
6 module module_bl_ysu
7 contains
9 !-------------------------------------------------------------------
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 !-------------------------------------------------------------------
114    integer,  intent(in   )   ::      ids,ide, jds,jde, kds,kde,                &
115                                      ims,ime, jms,jme, kms,kme,                &
116                                      its,ite, jts,jte, kts,kte
118    real,     intent(in   )   ::      dt,dtmin,cp,g,rovcp,rovg,rd,xlv,rv
120    real,     intent(in )     ::      svp1,svp2,svp3,svpt0
121    real,     intent(in )     ::      ep1,ep2,karman,eomeg,stbolt
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
136    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &
137              intent(inout)   ::                                       rublten, &
138                                                                       rvblten, &
139                                                                      rthblten, &
140                                                                      rqvblten, &
141                                                                      rqcblten
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
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
160    real,     dimension( ims:ime, jms:jme )                                   , &
161              intent(inout)   ::                                           hol, &
162                                                                           ust, &
163                                                                          hpbl, &
164                                                                           znt, &
165                                                                          wspd, &
166                                                                           zol
168   real,      dimension( ims:ime, kms:kme, jms:jme )                          , &
169              intent(in   )   ::                                           u3d, &
170                                                                           v3d
172   integer,   dimension( ims:ime, jms:jme )                                   , &
173              intent(out  )   ::                                        kpbl2d
175      logical, intent(in)        :: flag_qi
176 !optional
177    real,     dimension( ims:ime, jms:jme )                                   , &
178              optional                                                        , &
179              intent(inout)   ::                                        regime
181    real,     dimension( ims:ime, kms:kme, jms:jme )                          , &
182              optional                                                        , &
183              intent(inout)   ::                                       rqiblten
185    real,     dimension( kms:kme )                                            , &
186              optional                                                        , &
187              intent(in   )   ::                                           znu, &
188                                                                           znw
190    real,     dimension( ims:ime, jms:jme )                                   , &
191              optional                                                        , &
192              intent(in   )   ::                                           mut
194    real,     optional, intent(in   )   ::                               p_top
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
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
255    end subroutine ysu
257 !-------------------------------------------------------------------
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 !-------------------------------------------------------------------
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.
291 !     mrfpbl:
292 !     coded by song-you hong (ncep), implemented by jimy dudhia (ncar)
293 !     fall 1996
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
300 !     references:
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.
307 !-------------------------------------------------------------------
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
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
326    real,     intent(in   )   ::      dt,dtmin,cp,g,rovcp,rovg,rd,xlv,rv
328    real,     intent(in )     ::     svp1,svp2,svp3,svpt0
329    real,     intent(in )     ::     ep1,ep2,karman,eomeg,stbolt
331    real,     dimension( ims:ime, kms:kme ),                                    &
332              intent(in)      ::                                        dz8w2d, &
333                                                                           z2d
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
344    real,     dimension( its:ite, kts:kte )                                   , &
345              intent(in   )   ::                                           p2d
347    real,     dimension( its:ite, kts:kte )                                   , &
348              intent(inout)   ::                                         qitnp
350    real,     dimension( ims:ime, kms:kme )                                   , &
351              intent(inout)   ::                                          utnp, &
352                                                                          vtnp, &
353                                                                          ttnp, &
354                                                                          qtnp, &
355                                                                         qctnp
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
367    real,     dimension( ims:ime ), intent(inout)   ::                    wspd
368    real,     dimension( ims:ime ), intent(in  )    ::                      br
370    real,     dimension( ims:ime ), intent(in   )   ::                    psim, &
371                                                                          psih
372    real,     dimension( ims:ime ), intent(in   )   ::                  gz1oz0
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
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
388 ! local vars
390    real,     dimension( its:ite, kts:kte+1 ) ::                            zq
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
403    real,    dimension( its:ite )             ::                    qixsv,rhox, &
404                                                                        govrth, &
405                                                                         thxsv, &
406                                                                     uxsv,vxsv, &
407                                                                    qxsv,qcxsv, &
408                                                                  qgh,tgdsa,ps
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
421    real,    dimension( its:ite, kts:kte )    ::                     xkzm,xkzh, &
422                                                                         f1,f2, &
423                                                                         r1,r2, &
424                                                                         ad,au, &
425                                                                            cu, &
426                                                                            al, &
427                                                                          zfac
429 !jdf added exch_hx
430    real,    dimension( ims:ime, kms:kme )                                    , &
431             intent(inout)   ::                                        exch_hx
433    real,    dimension( ims:ime )                                             , & 
434             intent(in  )    ::                                            u10, &
435                                                                           v10
436    real,    dimension( its:ite )    ::                                         &
437                                                                     brcr_sbro
439    real,    dimension( its:ite, kts:kte, ncloud)  ::                    r3,f3
441    logical, dimension( its:ite )             ::                        pblflg, &
442                                                                        sfcflg, &
443                                                                        stable
445    integer ::  n,i,k,l,nzol,imvdif,ic
446    integer ::  klpbl
448    integer, dimension( its:ite )             ::      kpbl
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
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
472 !----------------------------------------------------------------------
474    klpbl = kte
476 !-- imvdif      imvdif = 1 for moist adiabat vertical diffusion
477    imvdif = 1
479 !----convert ground temperature to potential temperature:
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
487    do k = kts,kte
488      do i = its,ite
489        thx(i,k) = tx(i,k)/pi2d(i,k)
490      enddo
491    enddo
493    do i = its,ite
494      qgh(i) = 0.
495      cpm(i) = cp
496    enddo
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
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
511 !-----compute the height of full- and half-sigma levels above ground
512 !     level, and the layer thicknesses.
514    do i = its,ite
515      zq(i,1) = 0.
516      rhox(i) = ps(i)*1000./(rd*tx(i,1))
517    enddo
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
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
533       do i = its,ite
534         dza(i,1) = za(i,1)
535       enddo
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
543    do i = its,ite
544      govrth(i) = g/thx(i,1)
545    enddo
547 !-----initialize vertical tendencies and
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
557    do k = kts,kte
558      do i = its,ite
559        qtnp(i,k) = 0.
560      enddo
561    enddo
563    do k = kts,kte
564      do i = its,ite
565        qctnp(i,k) = 0.
566        qitnp(i,k) = 0.
567      enddo
568    enddo
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
574 !---- compute vertical diffusion
576 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
577 !     compute preliminary variables
579    dtstep = dt
580    dt2 = 2.*dtstep
581    rdt = 1./dt2
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
594    do k = kts,klpbl
595      do i = its,ite
596        wscalek(i,k) = 0.0
597      enddo
598    enddo
600    do k = kts,klpbl
601      do i = its,ite
602        zfac(i,k) = 0.0
603      enddo
604    enddo
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
619 !     compute the first guess of pbl height
621    do i = its,ite
622      stable(i) = .false.
623      brup(i) = br(i)
624    enddo
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
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
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
684 !     compute the surface variables for pbl height estimation
685 !     under unstable conditions
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
704 !     enhance the pbl height by considering the thermal
706    do i = its,ite
707      if(pblflg(i))then
708        kpbl(i) = 1
709        hpbl(i) = zq(i,1)
710      endif
711    enddo
713    do i = its,ite
714      if(pblflg(i))then
715        stable(i) = .false.
716        brup(i) = br(i)
717      endif
718    enddo
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
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
749 !     stable boundary layer
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
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
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
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
807 !     estimate the entrainment parameters
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
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
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
854 !     compute diffusion coefficients below pbl
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
878 !     compute diffusion coefficients over pbl (free atmosphere)
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
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
926 !     compute tridiagonal matrix elements for heat and moisture, and clouds
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
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
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
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
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
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
1011 ! copies here to avoid duplicate input args for tridin
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
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
1028 !     solve tridiagonal problem for heat and moisture, and clouds
1030    call tridin(al,ad,cu,r1,r3,au,f1,f3,its,ite,kts,kte,ncloud)
1032 !     recover tendencies of heat and moisture
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
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
1059 !     compute tridiagonal matrix elements for momentum
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
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
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
1112 ! copies here to avoid duplicate input args for tridin
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
1122 !     solve tridiagonal problem for momentum
1124    call tridin(al,ad,cu,r1,r2,au,f1,f2,its,ite,kts,kte,1)
1126 !     recover tendencies of momentum
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
1137 !---- end of vertical diffusion
1139    do i = its,ite
1140      kpbl1d(i) = kpbl(i)
1141    enddo
1143    end subroutine ysu2d
1145    subroutine tridin(cl,cm,cu,r1,r2,au,f1,f2,its,ite,kts,kte,nt)
1146 !----------------------------------------------------------------
1147    implicit none
1148 !----------------------------------------------------------------
1150    integer, intent(in )      ::     its,ite, kts,kte, nt
1152    real, dimension( its:ite, kts+1:kte+1 )                                   , &
1153          intent(in   )  ::                                                 cl
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
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
1168    real    :: fk
1169    integer :: i,k,l,n,it
1171 !----------------------------------------------------------------
1173    l = ite
1174    n = kte
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
1225    end subroutine tridin
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 !-------------------------------------------------------------------
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
1251    jtf = min0(jte,jde-1)
1252    ktf = min0(kte,kde-1)
1253    itf = min0(ite,ide-1)
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
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
1279    end subroutine ysuinit
1280 !-------------------------------------------------------------------
1281 end module module_bl_ysu