File: C:\NOAA\NEMS_11731\src\atmos\phys\moninq1.f

1     cfpp$ noconcur r
2           subroutine moninq1(ix,im,km,ntrac,dv,du,tau,rtg,
3          &     u1,v1,t1,q1,swh,hlw,xmu,
4          &     psk,rbsoil,fm,fh,tsea,qss,heat,evap,stress,spd1,kpbl,
5          &     prsi,del,prsl,prslk,phii,phil,deltim,
6          &     dusfc,dvsfc,dtsfc,dqsfc,hpbl,hgamt,hgamq,dkt,kinver)
7     !
8           use machine     , only : kind_phys
9           use funcphys , only : fpvs
10           use physcons, grav => con_g, rd => con_rd, cp => con_cp
11          &,             hvap => con_hvap, rog => con_rog, fv => con_fvirt
12          &,             eps => con_eps, epsm1 => con_epsm1
13           implicit none
14     !
15     !     include 'constant.h'
16     !
17     !
18     !     arguments
19     !
20           integer ix, im, km, ntrac, kpbl(im), kpblx(im)
21           integer kinver(im)
22     !
23           real(kind=kind_phys) deltim
24           real(kind=kind_phys) dv(im,km),     du(im,km),
25          &                     tau(im,km),    rtg(im,km,ntrac),
26          &                     u1(ix,km),     v1(ix,km),
27          &                     t1(ix,km),     q1(ix,km,ntrac),
28          &                     swh(ix,km),    hlw(ix,km),
29          &                     xmu(im), 
30          &                     psk(im),       rbsoil(im),
31     !    &                     cd(im),        ch(im),
32          &                     fm(im),        fh(im),
33          &                     tsea(im),      qss(im),
34          &                                    spd1(im),
35     !    &                     dphi(im),      spd1(im),
36          &                     prsi(ix,km+1), del(ix,km),
37          &                     prsl(ix,km),   prslk(ix,km),
38          &                     phii(ix,km+1), phil(ix,km),
39          &                     dusfc(im),
40          &                     dvsfc(im),     dtsfc(im),
41          &                     dqsfc(im),     hpbl(im),      hpblx(im),
42          &                     hgamt(im),     hgamq(im),
43          &                     hgamu(im),     hgamv(im),     hgams(im)
44     !
45     !    locals
46     !
47           integer i,iprt,is,iun,k,kk,km1,kmpbl,latd,lond
48           integer lcld(im),icld(im),kcld(im),krad(im)
49           integer kemx(im)
50     !
51     !     real(kind=kind_phys) betaq(im), betat(im),   betaw(im),
52           real(kind=kind_phys) evap(im),  heat(im),    phih(im),
53          &                     phim(im),  rbdn(im),    rbup(im),
54          &                     stress(im),beta(im),
55          &                     ustar(im), wscale(im),  thermal(im),
56          &                     ust3(im),  wstar3(im),  
57          &                     sumz(im),  sumt(im)
58          &,                    entflx(im),sflux(im)
59     !
60           real(kind=kind_phys) thlx(im,km), thlvx(im,km), tlx(im,km),
61          &                     thvx(im,km), qlx(im,km),
62          &                     qtx(im,km),  bf(im,km-1),
63          &                     govrth(im),  hrad(im),     radx(im,km-1),
64          &                     hradm(im),   radmin(im),   vrad(im),
65          &                     zd(im),      thlvx1(im)
66     !
67           real(kind=kind_phys) rdzt(im,km-1),dktx(im,km-1),dkux(im,km-1),
68          &                     zi(im,km+1),  zl(im,km),  xkzo(im,km),
69          &                     dku(im,km-1), dkt(im,km-1),
70          &                     dkuy(im,km-1),dkty(im,km-1),
71          &                     cku(im,km-1), ckt(im,km-1),
72          &                     al(im,km-1),  ad(im,km),
73          &                     au(im,km-1),  a1(im,km), 
74          &                     a2(im,km*ntrac), theta(im,km)
75     !
76           real(kind=kind_phys) hol(im), prinv(im), hpbl01(im)
77     !
78           logical  pblflg(im), sfcflg(im), scuflg(im), flg(im)
79     !
80           real(kind=kind_phys) aphi16,  aphi5,  bvf2,   wfac,
81          &                     cfac,    conq,   cont,   conw,
82          &                     conwrc,  dk,     dkmax,  dkmin,
83          &                     dq1,     dsdz2,  dsdzq,  dsdzt,
84          &                     dsdzu,   dsdzv,  sfac,
85          &                     dsig,    dt,     dthe1,  dtodsd,
86          &                     dtodsu,  dw2,    dw2min, g,
87          &                     gamcrq,  gamcrt, gocp,   gor, gravi,
88          &                     hol1,    pfac,   prmax,  prmin,
89          &                     prnum,   qmin,   tdzmin, qtend, rbcr,
90          &                     rbint,   rdt,    rdz,    qlmin, 
91     !    &                     rbint,   rdt,    rdz,    rdzt1,
92          &                     ri,      rimin,  rl2,    rlam,  rlamun,
93          &                     rone,    rzero,  sfcfrac,
94          &                     shr2,    spdk2,  sri,
95          &                     tem,     ti,     ttend,  tvd,
96          &                     tvu,     utend,  vk,     vk2,
97          &                     vtend,   zfac,   vpert,  cpert,
98          &                     entfac,  rentfac,radfac,
99          &                     zfmin,   zk,     tem1,   tem2,  xkzm,
100          &                     ptem,    ptem1,  ptem2
101     !
102           real(kind=kind_phys) zstblmax,h1,     h2,
103          &                     qlcr,    cldtime,alpri,  chiri,
104          &                     u01,     v01,    delu,   delv
105     cc
106           parameter(gravi=1.0/grav)
107           parameter(g=grav)
108           parameter(gor=g/rd,gocp=g/cp)
109           parameter(cont=cp/g,conq=hvap/g,conw=1.0/g)
110           parameter(alpri=hvap/rd,chiri=eps*hvap*hvap/cp/rd)
111           parameter(rlam=30.0,vk=0.4,vk2=vk*vk)
112           parameter(prmin=0.25,prmax=4.)
113           parameter(dw2min=0.0001,dkmin=0.0,dkmax=1000.,rimin=-100.)
114           parameter(rbcr=0.25,wfac=7.0,cfac=6.5,pfac=2.0,sfcfrac=0.1)
115     !     parameter(qmin=1.e-8,xkzm=1.0,zfmin=1.e-8,aphi5=5.,aphi16=16.)
116           parameter(qmin=1.e-8,xkzm=0.25,zfmin=1.e-8,aphi5=5.,aphi16=16.)
117           parameter(tdzmin=1.e-3,qlmin=1.e-12,cpert=0.25,sfac=5.4)
118           parameter(h1=0.33333333,h2=0.66666667)
119           parameter(cldtime=500.)
120     !     parameter(gamcrt=3.,gamcrq=2.e-3,rlamun=150.0)
121           parameter(gamcrt=3.,gamcrq=0.,rlamun=150.0)
122           parameter(entfac=0.2,rentfac=0.2,radfac=0.85)
123           parameter(iun=84)
124     !
125     !     parameter (zstblmax = 2500., qlcr=3.0e-5)
126     !     parameter (zstblmax = 2500., qlcr=3.5e-5)
127     !     parameter (zstblmax = 2500., qlcr=4.5e-5)
128           parameter (zstblmax = 2500., qlcr=5.0e-5)
129     c
130     c-----------------------------------------------------------------------
131     c
132      601  format(1x,' moninp lat lon step hour ',3i6,f6.1)
133      602      format(1x,'    k','        z','        t','       th',
134          1     '      tvh','        q','        u','        v',
135          2     '       sp')
136      603      format(1x,i5,8f9.1)
137      604      format(1x,'  sfc',9x,f9.1,18x,f9.1)
138      605      format(1x,'    k      zl    spd2   thekv   the1v'
139          1         ,' thermal    rbup')
140      606      format(1x,i5,6f8.2)
141      607      format(1x,' kpbl    hpbl      fm      fh   hgamt',
142          1         '   hgamq      ws   ustar      cd      ch')
143      608      format(1x,i5,9f8.2)
144      609      format(1x,' k pr dkt dku ',i5,3f8.2)
145      610      format(1x,' k pr dkt dku ',i5,3f8.2,' l2 ri t2',
146          1         ' sr2  ',2f8.2,2e10.2)
147     c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148     c     compute preliminary variables
149     c
150           if (ix .lt. im) stop
151     !
152     !     iprt = 0
153     !     if(iprt.eq.1) then
154     ccc   latd = 0
155     !     lond = 0
156     !     else
157     ccc   latd = 0
158     !     lond = 0
159     !     endif
160     c
161           dt    = 2. * deltim
162           rdt   = 1. / dt
163           km1   = km - 1
164           kmpbl = km / 2
165     !
166           do k=1,km
167             do i=1,im
168               zi(i,k) = phii(i,k) * gravi
169               zl(i,k) = phil(i,k) * gravi
170             enddo
171           enddo
172           do i=1,im
173              zi(i,km+1) = phii(i,km+1) * gravi
174           enddo
175     c
176           do k = 1,km1
177             do i=1,im
178               rdzt(i,k) = 1.0 / (zl(i,k+1) - zl(i,k))
179             enddo
180           enddo
181     c
182           do k = 1,km1
183             do i=1,im
184               if (k .lt. kinver(i)) then
185                 tem1      = 1.0 - prsi(i,k+1) / prsi(i,1)
186                 tem1      = tem1 * tem1 * 10.0
187                 xkzo(i,k) = xkzm * min(1.0, exp(-tem1))
188               else
189                 xkzo(i,k) = 0.0
190               endif
191             enddo
192           enddo
193     c
194           do i = 1,im
195              dusfc(i) = 0.
196              dvsfc(i) = 0.
197              dtsfc(i) = 0.
198              dqsfc(i) = 0.
199              hgamt(i) = 0.
200              hgamq(i) = 0.
201              hgamu(i) = 0.
202              hgamv(i) = 0.
203              hgams(i) = 0.
204              wscale(i)= 0.
205              entflx(i)= 0.
206              kpbl(i)  = 1
207              hpbl(i)  = zi(i,1)
208              pblflg(i)= .true.
209              sfcflg(i)= .true.
210              if(rbsoil(i).gt.0.0) sfcflg(i) = .false.
211              sumz(i)  = 0.
212              sumt(i)  = 0.
213              scuflg(i)= .true.
214              if(scuflg(i)) then
215                radmin(i)= 0.
216                hrad(i)  = zi(i,1)
217                icld(i)  = 0
218                lcld(i)  = km1
219                kcld(i)  = km1
220                krad(i)  = km1
221                zd(i)    = 0.
222             endif
223           enddo
224     !
225           do k = 1,km
226             do i = 1,im
227               theta(i,k) = t1(i,k) * psk(i) / prslk(i,k)
228               qlx(i,k)   = max(q1(i,k,ntrac),qlmin)
229               qtx(i,k)   = max(q1(i,k,1),qmin)+qlx(i,k)
230               ptem       = qlx(i,k)
231               ptem1      = hvap*max(q1(i,k,1),qmin)/(cp*t1(i,k))
232               thvx(i,k)  = theta(i,k)*(1.+fv*max(q1(i,k,1),qmin)-ptem)
233               tlx(i,k)   = t1(i,k)-(hvap/cp)*ptem
234               thlx(i,k)  = theta(i,k)-(hvap/cp)*ptem
235               thlvx(i,k) = thlx(i,k)*(1.+fv*qtx(i,k))
236             enddo
237           enddo
238           do k = 1,km1
239             do i = 1,im
240               dku(i,k)  = 0.
241               dkt(i,k)  = 0.
242               dktx(i,k) = 0.
243               dkux(i,k) = 0.
244               dkty(i,k) = 0.
245               dkuy(i,k) = 0.
246               cku(i,k)  = 0.
247               ckt(i,k)  = 0.
248               tem       = zi(i,k+1)-zi(i,k)
249               radx(i,k) = tem*(swh(i,k)*xmu(i)+hlw(i,k))
250             enddo
251           enddo
252     c
253           do i=1,im
254              flg(i)  = scuflg(i)
255           enddo
256           do k = 1, km1
257             do i=1,im
258               if(flg(i).and.zl(i,k).ge.zstblmax) then
259                  lcld(i)=k
260                  flg(i)=.false.
261               endif
262           enddo
263           enddo
264     c
265     c  compute buoyancy flux
266     c
267           do k = 1, km1
268           do i = 1, im
269              bf(i,k) = (thvx(i,k+1)-thvx(i,k))*rdzt(i,k)
270           enddo
271           enddo
272     c
273           do i = 1,im
274             govrth(i) = g/theta(i,1)
275           enddo
276     c
277           do i=1,im
278              beta(i)  = dt / (zi(i,2)-zi(i,1))
279           enddo
280     c
281           do i=1,im
282              ustar(i) = sqrt(stress(i))
283           enddo
284     c
285     c  compute the first guess pbl height
286     c
287           do i = 1,im
288              sflux(i)  = heat(i) + evap(i)*fv*theta(i,1)
289              if(.not.sfcflg(i).or.sflux(i).le.0.0) pblflg(i)=.false.
290           enddo
291     c
292           do i=1,im
293              flg(i) = .true.
294              if(pblflg(i)) then
295                 flg(i)  = .false.
296                 sumz(i) = zl(i,1)
297                 if(scuflg(i)) then
298                   rbup(i) = thlvx(i,1)
299                 else
300                   rbup(i) = thvx(i,1)
301                 endif
302                 sumt(i) = rbup(i)*zl(i,1)
303              endif
304           enddo
305           do k = 2, kmpbl
306           do i = 1, im
307             if(.not.flg(i)) then
308               rbdn(i) = rbup(i)
309               tem     = zl(i,k)-zl(i,k-1)
310               sumz(i) = sumz(i)+tem
311               if(scuflg(i)) then
312                 tem1    = 0.5*(thlvx(i,k)+thlvx(i,k-1))
313                 rbup(i) = thlvx(i,k)
314               else
315                 tem1    = 0.5*(thvx(i,k)+thvx(i,k-1))
316                 rbup(i) = thvx(i,k)
317               endif
318               sumt(i) = sumt(i)+tem1*tem
319               thermal(i)= sumt(i)/sumz(i)+cpert
320               kpbl(i) = k
321               flg(i)  = rbup(i).gt.thermal(i)
322             endif
323           enddo
324           enddo
325           do i = 1,im
326             if(pblflg(i)) then
327               k = kpbl(i)
328               if(rbdn(i).ge.thermal(i)) then
329                  rbint = 0.
330               elseif(rbup(i).le.thermal(i)) then
331                  rbint = 1.
332               else
333                  rbint = (thermal(i)-rbdn(i))/(rbup(i)-rbdn(i))
334               endif
335               hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
336               kpbl(i) = kpbl(i) - 1
337             endif
338           enddo
339     c
340           do i = 1, im
341             if(pblflg(i)) then
342                hpbl01(i) = sfcfrac*hpbl(i)
343                if(scuflg(i)) then
344                  thermal(i) = thlvx(i,1)
345                else
346                  thermal(i) = thvx(i,1)
347                endif
348             endif
349           enddo
350     c
351     c     compute similarity parameters 
352     c
353           do i=1,im
354              if(pblflg(i)) then
355                hol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin)
356                hol(i) = min(hol(i),-zfmin)
357     c
358                hol1 = hol(i)*hpbl(i)/zl(i,1)*sfcfrac
359     !          phim(i) = (1.-aphi16*hol1)**(-1./4.)
360     !          phih(i) = (1.-aphi16*hol1)**(-1./2.)
361                tem     = 1.0 / (1. - aphi16*hol1)
362                phih(i) = sqrt(tem)
363                phim(i) = sqrt(phih(i))
364                ptem = max(heat(i)+fv*theta(i,1)*evap(i),0.)
365                wstar3(i) = govrth(i)*ptem*hpbl(i)
366                ust3(i)   = ustar(i)**3.
367                wscale(i) = (ust3(i)+wfac*vk*wstar3(i)*sfcfrac)**h1
368     !          wscale(i) = ustar(i)/phim(i)
369                wscale(i) = min(wscale(i),ustar(i)*aphi16)
370                wscale(i) = max(wscale(i),ustar(i)/aphi5)
371              endif
372           enddo
373     c
374     c compute counter-gradient mixing term for heat and moisture
375     c
376           do i = 1,im
377              if(pblflg(i)) then
378                hgamt(i)  = min(cfac*heat(i)/wscale(i),gamcrt)
379                hgamq(i)  = min(cfac*evap(i)/wscale(i),gamcrq)
380                vpert     = hgamt(i) + hgamq(i)*fv*theta(i,1)
381                vpert     = min(vpert,gamcrt)
382                thermal(i)= thermal(i)+max(vpert,0.)
383                hgamt(i)  = max(hgamt(i),0.0)
384                hgamq(i)  = max(hgamq(i),0.0)
385              endif
386           enddo
387     c
388     c compute large-scale mixing term for momentum
389     c
390           do i = 1,im
391             flg(i) = pblflg(i)
392             kemx(i)= 1
393           enddo
394           do k = 1, kmpbl
395           do i = 1, im
396             if(flg(i).and.zl(i,k).gt.hpbl01(i)) then
397               kemx(i) = k
398               flg(i)  = .false.
399             endif
400           enddo
401           enddo
402           do i = 1, im
403             if(pblflg(i)) then
404               kk = kpbl(i)
405               if(kemx(i).le.1) then
406                 ptem  = u1(i,1)/zl(i,1)
407                 ptem1 = v1(i,1)/zl(i,1)
408                 u01   = ptem*hpbl01(i)
409                 v01   = ptem1*hpbl01(i)
410               else
411                 tem   = zl(i,kemx(i))-zl(i,kemx(i)-1)
412                 ptem  = (u1(i,kemx(i))-u1(i,kemx(i)-1))/tem
413                 ptem1 = (v1(i,kemx(i))-v1(i,kemx(i)-1))/tem
414                 tem1  = hpbl01(i)-zl(i,kemx(i)-1)
415                 u01   = u1(i,kemx(i)-1)+ptem*tem1
416                 v01   = v1(i,kemx(i)-1)+ptem1*tem1
417               endif
418               delu  = u1(i,kk)-u01
419               delv  = v1(i,kk)-v01
420               tem2  = sqrt(delu**2+delv**2)
421               tem2  = max(tem2,0.1)
422               ptem2 = -sfac*ustar(i)*ustar(i)*wstar3(i)
423          1                /(wscale(i)**4.)
424               hgamu(i) = ptem2*delu/tem2 
425               hgamv(i) = ptem2*delv/tem2
426               tem  = sqrt(u1(i,kk)**2+v1(i,kk)**2)
427               tem1 = sqrt(u01**2+v01**2)
428               ptem = tem - tem1
429               if(ptem.gt.0.) then
430                 hgams(i) = -sfac*vk*sfcfrac*wstar3(i)/(wscale(i)**3.)
431               else
432                 hgams(i) = sfac*vk*sfcfrac*wstar3(i)/(wscale(i)**3.)
433               endif
434             endif
435           enddo
436     c
437     c  enhance the pbl height by considering the thermal excess
438     c
439           do i=1,im
440              flg(i)  = .true.
441              if(pblflg(i)) then
442                flg(i)  = .false.
443                if(scuflg(i)) then
444                  rbup(i) = thlvx(i,1)
445                else
446                  rbup(i) = thvx(i,1)
447                endif
448              endif
449           enddo
450           do k = 2, kmpbl
451           do i = 1, im
452             if(.not.flg(i)) then
453               rbdn(i) = rbup(i)
454               if(scuflg(i)) then
455                 rbup(i) = thlvx(i,k)
456               else
457                 rbup(i) = thvx(i,k)
458               endif
459               kpblx(i) = k
460               flg(i)  = rbup(i).gt.thermal(i)
461             endif
462           enddo
463           enddo
464           do i = 1,im
465             if(pblflg(i)) then
466                k = kpblx(i)
467                if(rbdn(i).ge.thermal(i)) then
468                   rbint = 0.
469                elseif(rbup(i).le.thermal(i)) then
470                   rbint = 1.
471                else
472                   rbint = (thermal(i)-rbdn(i))/(rbup(i)-rbdn(i))
473                endif
474                hpblx(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
475                kpblx(i) = kpblx(i) - 1
476                if(hpblx(i).gt.hpbl(i)) then
477                   hpbl(i) = hpblx(i)
478                   kpbl(i) = kpblx(i)
479                   if(kpbl(i).le.1) pblflg(i)=.false.
480                endif
481             endif
482           enddo
483     c
484     c  look for stratocumulus-topped pbl
485     c
486           do i = 1, im
487             flg(i)=scuflg(i)
488           enddo
489           do k = kmpbl,1,-1
490           do i = 1, im
491             if(flg(i).and.k.le.lcld(i)) then
492               if(qlx(i,k).ge.qlcr) then
493                  kcld(i)=k
494                  flg(i)=.false.
495               endif
496             endif
497           enddo
498           enddo
499           do i = 1, im
500             if(scuflg(i).and.kcld(i).eq.km1) scuflg(i)=.false.
501           enddo
502     c
503           do i = 1, im
504             flg(i)=scuflg(i)
505           enddo
506           do k = kmpbl,1,-1
507           do i = 1, im
508             if(flg(i).and.k.le.kcld(i)) then
509               if(qlx(i,k).ge.qlcr) then
510                 if(radx(i,k).lt.radmin(i)) then
511                   radmin(i)=radx(i,k)
512                   krad(i)=k
513                 endif
514               else
515                 flg(i)=.false.
516               endif
517             endif
518           enddo
519           enddo
520           do i = 1, im
521             if(scuflg(i).and.krad(i).eq.km1) scuflg(i)=.false.
522             if(scuflg(i).and.krad(i).le.1) scuflg(i)=.false.
523             if(scuflg(i).and.radmin(i).ge.0.) scuflg(i)=.false.
524           enddo
525     c
526           do i = 1, im
527             flg(i)=scuflg(i)
528           enddo
529           do k = kmpbl,2,-1
530           do i = 1, im
531             if(flg(i).and.k.le.krad(i)) then
532               if(qlx(i,k).ge.qlcr) then
533                 icld(i)=icld(i)+1
534               else
535                 flg(i)=.false.
536               endif
537             endif
538           enddo
539           enddo
540           do i = 1, im
541             if(scuflg(i).and.icld(i).lt.2) scuflg(i)=.false.
542           enddo
543     c
544           do i = 1, im
545             if(scuflg(i)) then
546                hrad(i) = zi(i,krad(i)+1)
547                hradm(i)= zl(i,krad(i))
548             endif
549           enddo
550     c
551           do i = 1, im
552             if(scuflg(i).and.hrad(i).lt.200.) scuflg(i)=.false.
553           enddo
554     c
555           do i = 1, im
556             if(scuflg(i)) then
557               k    = krad(i)
558               tem  = zi(i,k+1)-zi(i,k)
559               tem1 = cldtime*radmin(i)/tem
560               thlvx1(i) = thlvx(i,k)+tem1
561               if(thlvx1(i).gt.thlvx(i,k-1)) scuflg(i)=.false.
562             endif
563           enddo
564     c 
565           do i = 1, im
566              flg(i)=scuflg(i)
567           enddo
568           do k = kmpbl,1,-1
569           do i = 1, im
570             if(flg(i).and.k.le.krad(i))then
571               if(thlvx1(i).le.thlvx(i,k))then
572                  tem=zi(i,k+1)-zi(i,k)
573                  zd(i)=zd(i)+tem
574               else
575                  flg(i)=.false.
576               endif
577             endif
578           enddo
579           enddo
580           do i = 1, im
581             if(scuflg(i))then
582               zd(i) = min(zd(i),hrad(i))
583               tem   = govrth(i)*zd(i)*(-radmin(i))
584               vrad(i)= tem**h1
585             endif
586           enddo
587     c
588           do i = 1, im
589             if(scuflg(i).and.pblflg(i)) then
590               if(hpbl(i).ge.hradm(i)) then
591                  hpbl(i)=hrad(i)
592                  kpbl(i)=krad(i)
593               endif
594             endif
595           enddo
596     c
597     c     compute inverse Prandtl number
598     c
599           do i = 1, im
600             if(pblflg(i)) then
601               tem = phih(i)/phim(i)+cfac*vk*sfcfrac
602               prinv(i) = (1.0-hgams(i))/tem
603               prinv(i) = min(prinv(i),prmax)
604               prinv(i) = max(prinv(i),prmin)
605             endif
606           enddo
607     c
608     c     compute entrainment flux at pbl top
609     c
610           do i = 1, im
611             if(pblflg(i)) then
612                ptem=(theta(i,1)/g)*ust3(i)
613                entflx(i)=entfac*(sflux(i)+5.*ptem/hpbl(i))
614             endif
615           enddo
616     c
617     c     compute diffusion coefficients below pbl
618     c
619           do k = 1, kmpbl
620           do i=1,im
621              if(pblflg(i).and.k.lt.kpbl(i)) then
622     !           zfac = max((1.-(zi(i,k+1)-zl(i,1))/
623     !    1             (hpbl(i)-zl(i,1))), zfmin)
624                 zfac = max((1.-zi(i,k+1)/hpbl(i)), zfmin)
625                 dku(i,k) = wscale(i)*vk*zi(i,k+1)
626          1                 *zfac**pfac
627     !           dku(i,k) = xkzo(i,k)+wscale(i)*vk*zi(i,k+1)
628     !    1                 *zfac**pfac
629                 dkt(i,k) = dku(i,k)*prinv(i)
630                 dku(i,k) = min(dku(i,k),dkmax)
631                 dku(i,k) = max(dku(i,k),xkzo(i,k))
632                 dkt(i,k) = min(dkt(i,k),dkmax)
633                 dkt(i,k) = max(dkt(i,k),xkzo(i,k))
634                 dktx(i,k)= dkt(i,k)
635                 dkux(i,k)= dku(i,k)
636              endif
637           enddo
638           enddo
639     c
640           do i = 1, im
641             if(pblflg(i)) then
642                k = kpbl(i)
643                if(bf(i,k).gt.0.) then
644                  ptem = max(bf(i,k),tdzmin)
645                  dkt(i,k) = entflx(i)/ptem
646                  dkt(i,k) = min(dkt(i,k),dkmax)
647                  dku(i,k) = dkt(i,k)
648                endif
649             endif
650           enddo
651     c
652     c compute diffusion coefficients based on local scheme,
653     c which are applied to nighttime stable boundary layer and
654     c free atmosphere over pbl in the daytime unstable condition 
655     c
656           do k = 1, km1
657              do i=1,im
658     !!          if(k.ge.kpbl(i)) then
659                    rdz  = rdzt(i,k)
660                    dw2  = ((u1(i,k)-u1(i,k+1))**2
661          &              + (v1(i,k)-v1(i,k+1))**2)
662                    shr2 = max(dw2,dw2min)*rdz*rdz
663     !              if(qlx(i,k).ge.qlcr.and.qlx(i,k+1).ge.qlcr) then
664                    if(qlx(i,k).ge.qlcr.or.qlx(i,k+1).ge.qlcr) then
665                      ti   = 2./(t1(i,k)+t1(i,k+1))
666                      tem  = .5*(max(q1(i,k,1),qmin)+max(q1(i,k+1,1),qmin))
667                      tem1 = alpri*tem*ti
668                      tem2 = chiri*tem*ti*ti
669                      ptem = (tem2-tem1)/(1.+tem2)
670                      ptem1= bf(i,k)/thvx(i,k+1)-g*ptem*ti/cp
671                      bvf2 = (1.+tem1)*g*ptem1
672                    else
673                      bvf2 = g*bf(i,k)/thvx(i,k+1)
674                    endif 
675                    ri   = max(bvf2/shr2,rimin)
676                    zk   = vk*zi(i,k+1)
677                    if(ri.lt.0.) then ! unstable regime
678                       rl2      = zk*rlamun/(rlamun+zk)
679                       dk       = rl2*rl2*sqrt(shr2)
680                       sri      = sqrt(-ri)
681     !                 dkuy(i,k) = xkzo(i,k) + dk*(1+8.*(-ri)/(1+1.746*sri))
682     !                 dkty(i,k) = xkzo(i,k) + dk*(1+8.*(-ri)/(1+1.286*sri))
683                       dkuy(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri))
684                       dkty(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri))
685                    else             ! stable regime
686                       rl2      = zk*rlam/(rlam+zk)
687     !!                tem      = rlam * sqrt(0.01*prsi(i,k))
688     !!                rl2      = zk*tem/(tem+zk)
689                       dk       = rl2*rl2*sqrt(shr2)
690     !                 dkty(i,k)= xkzo(i,k) + dk/(1+5.*ri)**2
691                       dkty(i,k)= dk/(1.+5.*ri)**2
692                       prnum = 1.0 + 2.1*ri
693                       prnum = min(prnum,prmax)
694     !                 dkuy(i,k)= (dkty(i,k)-xkzo(i,k))*prnum + xkzo(i,k)
695                       dkuy(i,k)= dkty(i,k)*prnum
696                    endif
697     c
698                    dkuy(i,k) = min(dkuy(i,k),dkmax)
699                    dkuy(i,k) = max(dkuy(i,k),xkzo(i,k))
700                    dkty(i,k) = min(dkty(i,k),dkmax)
701                    dkty(i,k) = max(dkty(i,k),xkzo(i,k))
702     c
703     !!          endif
704     c
705              enddo
706           enddo
707     c
708     c  compute diffusion coefficients for cloud-top driven diffusion
709     c
710           do i = 1, im
711             if(scuflg(i)) then
712                k = krad(i)
713                if(bf(i,k).gt.0.) then
714                  tem1  = max(bf(i,k),tdzmin)
715                  ckt(i,k) = -rentfac*radmin(i)/tem1
716                  cku(i,k) = ckt(i,k)
717                endif
718             endif
719           enddo
720     c
721           do k = 1, kmpbl
722              do i=1,im
723                 if(scuflg(i).and.k.lt.krad(i)) then
724                    tem1=hrad(i)-zd(i)
725                    tem2=zi(i,k+1)-tem1
726                    if(tem2.gt.0.) then
727                       ptem= tem2/zd(i)
728                       if(ptem.ge.1.) ptem= 1.
729                       ptem= tem2*ptem*sqrt(1.-ptem)
730                       ckt(i,k) = radfac*vk*vrad(i)*ptem
731                       cku(i,k) = 0.75*ckt(i,k)
732                       ckt(i,k) = max(ckt(i,k),dkmin)
733                       ckt(i,k) = min(ckt(i,k),dkmax)
734                       cku(i,k) = max(cku(i,k),dkmin)
735                       cku(i,k) = min(cku(i,k),dkmax)
736                    endif
737                 endif
738              enddo
739           enddo
740     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
741     !
742           do k = 1, kmpbl
743             do i=1,im
744               if(scuflg(i)) then
745                  dkt(i,k) = dkt(i,k)+ckt(i,k)
746                  dku(i,k) = dku(i,k)+cku(i,k)
747                  dkt(i,k) = min(dkt(i,k),dkmax)
748                  dku(i,k) = min(dku(i,k),dkmax)
749               endif
750             enddo
751           enddo
752     c
753           do k = 1, km1
754             do i=1,im
755                dkt(i,k) = max(dkt(i,k),dkty(i,k))
756                dku(i,k) = max(dku(i,k),dkuy(i,k))
757             enddo
758           enddo
759     c
760     c     compute tridiagonal matrix elements for heat and moisture
761     c
762           do i=1,im
763              ad(i,1) = 1.
764              a1(i,1) = t1(i,1)   + beta(i) * heat(i)
765              a2(i,1) = q1(i,1,1) + beta(i) * evap(i)
766           enddo
767           if(ntrac.ge.2) then
768             do k = 2, ntrac
769               is = (k-1) * km
770               do i = 1, im
771                 a2(i,1+is) = q1(i,1,k)
772               enddo
773             enddo
774           endif
775     c
776           do k = 1,km1
777             do i = 1,im
778               dtodsd = dt/del(i,k)
779               dtodsu = dt/del(i,k+1)
780               dsig   = prsl(i,k)-prsl(i,k+1)
781     !         rdz    = rdzt(i,k)*2./(t1(i,k)+t1(i,k+1))
782               rdz    = rdzt(i,k)
783               tem1   = dsig * dkt(i,k) * rdz
784               if(pblflg(i).and.k.lt.kpbl(i)) then
785     !            dsdzt = dsig*dkt(i,k)*rdz*(gocp-hgamt(i)/hpbl(i))
786     !            dsdzq = dsig*dkt(i,k)*rdz*(-hgamq(i)/hpbl(i))
787                  ptem1 = dsig * dktx(i,k) * rdz
788                  tem   = 1.0 / hpbl(i)
789                  dsdzt = tem1 * gocp - ptem1*hgamt(i)*tem
790                  dsdzq = ptem1 * (-hgamq(i)*tem)
791                  a2(i,k)   = a2(i,k)+dtodsd*dsdzq
792                  a2(i,k+1) = q1(i,k+1,1)-dtodsu*dsdzq
793               else
794     !            dsdzt = dsig*dkt(i,k)*rdz*(gocp)
795                  dsdzt = tem1 * gocp
796                  a2(i,k+1) = q1(i,k+1,1)
797               endif
798     !         dsdz2 = dsig*dkt(i,k)*rdz*rdz
799               dsdz2     = tem1 * rdz
800               au(i,k)   = -dtodsd*dsdz2
801               al(i,k)   = -dtodsu*dsdz2
802               ad(i,k)   = ad(i,k)-au(i,k)
803               ad(i,k+1) = 1.-al(i,k)
804               a1(i,k)   = a1(i,k)+dtodsd*dsdzt
805               a1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt
806             enddo
807           enddo
808           if(ntrac.ge.2) then
809             do kk = 2, ntrac
810               is = (kk-1) * km
811               do k = 1, km1
812                 do i = 1, im
813                   a2(i,k+1+is) = q1(i,k+1,kk)
814                 enddo
815               enddo
816             enddo
817           endif
818     c
819     c     solve tridiagonal problem for heat and moisture
820     c
821           call tridin(im,km,ntrac,al,ad,au,a1,a2,au,a1,a2)
822     c
823     c     recover tendencies of heat and moisture
824     c
825           do  k = 1,km
826              do i = 1,im
827                 ttend      = (a1(i,k)-t1(i,k))*rdt
828                 qtend      = (a2(i,k)-q1(i,k,1))*rdt
829                 tau(i,k)   = tau(i,k)+ttend
830                 rtg(i,k,1) = rtg(i,k,1)+qtend
831                 dtsfc(i)   = dtsfc(i)+cont*del(i,k)*ttend
832                 dqsfc(i)   = dqsfc(i)+conq*del(i,k)*qtend
833              enddo
834           enddo
835           if(ntrac.ge.2) then
836             do kk = 2, ntrac
837               is = (kk-1) * km
838               do k = 1, km 
839                 do i = 1, im
840                   qtend = (a2(i,k+is)-q1(i,k,kk))*rdt
841                   rtg(i,k,kk) = rtg(i,k,kk)+qtend
842                 enddo
843               enddo
844             enddo
845           endif
846     c
847     c     compute tridiagonal matrix elements for momentum
848     c
849           do i=1,im
850              ad(i,1) = 1.0 + beta(i) * stress(i) / spd1(i)
851              a1(i,1) = u1(i,1)
852              a2(i,1) = v1(i,1)
853           enddo
854     c
855           do k = 1,km1
856             do i=1,im
857               dtodsd = dt/del(i,k)
858               dtodsu = dt/del(i,k+1)
859               dsig   = prsl(i,k)-prsl(i,k+1)
860               rdz    = rdzt(i,k)
861               tem1   = dsig*dku(i,k)*rdz
862               if(pblflg(i).and.k.lt.kpbl(i))then
863                 ptem1 = dsig*dkux(i,k)*rdz
864                 dsdzu = ptem1*(-hgamu(i)/hpbl(i))
865                 dsdzv = ptem1*(-hgamv(i)/hpbl(i))
866                 a1(i,k)   = a1(i,k)+dtodsd*dsdzu
867                 a1(i,k+1) = u1(i,k+1)-dtodsu*dsdzu
868                 a2(i,k)   = a2(i,k)+dtodsd*dsdzv
869                 a2(i,k+1) = v1(i,k+1)-dtodsu*dsdzv
870               else
871                 a1(i,k+1) = u1(i,k+1)
872                 a2(i,k+1) = v1(i,k+1)
873               endif
874     !         dsdz2     = dsig*dku(i,k)*rdz*rdz
875               dsdz2     = tem1*rdz
876               au(i,k)   = -dtodsd*dsdz2
877               al(i,k)   = -dtodsu*dsdz2
878               ad(i,k)   = ad(i,k)-au(i,k)
879               ad(i,k+1) = 1.-al(i,k)
880             enddo
881           enddo
882     c
883     c     solve tridiagonal problem for momentum
884     c
885           call tridi2(im,km,al,ad,au,a1,a2,au,a1,a2)
886     c
887     c     recover tendencies of momentum
888     c
889           do k = 1,km
890              do i = 1,im
891                 utend    = (a1(i,k)-u1(i,k))*rdt
892                 vtend    = (a2(i,k)-v1(i,k))*rdt
893                 du(i,k)  = du(i,k)  + utend
894                 dv(i,k)  = dv(i,k)  + vtend
895                 dusfc(i) = dusfc(i) + conw*del(i,k)*utend
896                 dvsfc(i) = dvsfc(i) + conw*del(i,k)*vtend
897              enddo
898           enddo
899     c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
900     c  estimate the pbl height for diagnostic purpose
901     c
902           do i = 1, im
903              if(scuflg(i)) then
904                 thermal(i) = thlvx(i,1)
905              else
906                 thermal(i) = thvx(i,1)
907              endif
908              flg(i) = .false.
909              rbup(i)= rbsoil(i)
910           enddo
911           do k = 2, kmpbl
912           do i = 1, im
913             if(.not.flg(i)) then
914               rbdn(i)   = rbup(i)
915               spdk2     = max((u1(i,k)**2+v1(i,k)**2),1.)
916               if(scuflg(i)) then
917                  rbup(i)   = (thlvx(i,k)-thermal(i))*
918          &                   (g*zl(i,k)/thlvx(i,1))/spdk2
919               else
920                  rbup(i)   = (thvx(i,k)-thermal(i))*
921          &                   (g*zl(i,k)/thvx(i,1))/spdk2
922               endif
923               kpbl(i)= k
924               flg(i) = rbup(i).gt.rbcr
925             endif
926           enddo
927           enddo
928     c
929           do i = 1,im
930              k = kpbl(i)
931              if(rbdn(i).ge.rbcr) then
932                 rbint = 0.
933              elseif(rbup(i).le.rbcr) then
934                 rbint = 1.
935              else
936                 rbint = (rbcr-rbdn(i))/(rbup(i)-rbdn(i))
937              endif
938              hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
939              if(hpbl(i).lt.zi(i,kpbl(i))) kpbl(i) = kpbl(i) - 1
940           enddo
941     c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
942           return
943           end
944