GFS Operational Physics Documentation  gsm/branches/DTC/phys-doc-all phys-doc-all R82971
moninedmf.f
Go to the documentation of this file.
1 
4 
15 
89  subroutine moninedmf(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, &
90  & u1,v1,t1,q1,swh,hlw,xmu, &
91  & psk,rbsoil,zorl,u10m,v10m,fm,fh, &
92  & tsea,qss,heat,evap,stress,spd1,kpbl, &
93  & prsi,del,prsl,prslk,phii,phil,delt,dspheat, &
94  & dusfc,dvsfc,dtsfc,dqsfc,hpbl,hgamt,hgamq,dkt, &
95  & kinver,xkzm_m,xkzm_h,xkzm_s,lprnt,ipr)
96 !
97  use machine , only : kind_phys
98  use funcphys , only : fpvs
99  use physcons, grav => con_g, rd => con_rd, cp => con_cp
100  &, hvap => con_hvap, fv => con_fvirt
101  implicit none
102 !
103 ! arguments
104 !
105  logical lprnt
106  integer ipr
107  integer ix, im, km, ntrac, ntcw, kpbl(im), kinver(im)
108 !
109  real(kind=kind_phys) delt, xkzm_m, xkzm_h, xkzm_s
110  real(kind=kind_phys) dv(im,km), du(im,km), &
111  & tau(im,km), rtg(im,km,ntrac), &
112  & u1(ix,km), v1(ix,km), &
113  & t1(ix,km), q1(ix,km,ntrac), &
114  & swh(ix,km), hlw(ix,km), &
115  & xmu(im), psk(im), &
116  & rbsoil(im), zorl(im), &
117  & u10m(im), v10m(im), &
118  & fm(im), fh(im), &
119  & tsea(im), qss(im), &
120  & spd1(im), &
121  & prsi(ix,km+1), del(ix,km), &
122  & prsl(ix,km), prslk(ix,km), &
123  & phii(ix,km+1), phil(ix,km), &
124  & dusfc(im), dvsfc(im), &
125  & dtsfc(im), dqsfc(im), &
126  & hpbl(im), hpblx(im), &
127  & hgamt(im), hgamq(im)
128 !
129  logical dspheat
130 ! flag for tke dissipative heating
131 !
132 ! locals
133 !
134  integer i,iprt,is,iun,k,kk,km1,kmpbl,latd,lond
135  integer lcld(im),icld(im),kcld(im),krad(im)
136  integer kx1(im), kpblx(im)
137 !
138 ! real(kind=kind_phys) betaq(im), betat(im), betaw(im),
139  real(kind=kind_phys) evap(im), heat(im), phih(im), &
140  & phim(im), rbdn(im), rbup(im), &
141  & stress(im),beta(im), sflux(im), &
142  & z0(im), crb(im), wstar(im), &
143  & zol(im), ustmin(im), ustar(im), &
144  & thermal(im),wscale(im), wscaleu(im)
145 !
146  real(kind=kind_phys) theta(im,km),thvx(im,km), thlvx(im,km), &
147  & qlx(im,km), thetae(im,km), &
148  & qtx(im,km), bf(im,km-1), diss(im,km), &
149  & radx(im,km-1), &
150  & govrth(im), hrad(im), &
151 ! & hradm(im), radmin(im), vrad(im), &
152  & radmin(im), vrad(im), &
153  & zd(im), zdd(im), thlvx1(im)
154 !
155  real(kind=kind_phys) rdzt(im,km-1),dktx(im,km-1), &
156  & zi(im,km+1), zl(im,km), xkzo(im,km-1), &
157  & dku(im,km-1), dkt(im,km-1), xkzmo(im,km-1), &
158  & cku(im,km-1), ckt(im,km-1), &
159  & ti(im,km-1), shr2(im,km-1), &
160  & al(im,km-1), ad(im,km), &
161  & au(im,km-1), a1(im,km), &
162  & a2(im,km*ntrac)
163 !
164  real(kind=kind_phys) tcko(im,km), qcko(im,km,ntrac), &
165  & ucko(im,km), vcko(im,km), xmf(im,km)
166 !
167  real(kind=kind_phys) prinv(im), rent(im)
168 !
169  logical pblflg(im), sfcflg(im), scuflg(im), flg(im)
170  logical ublflg(im), pcnvflg(im)
171 !
172 ! pcnvflg: true for convective(strongly unstable) pbl
173 ! ublflg: true for unstable but not convective(strongly unstable) pbl
174 !
175  real(kind=kind_phys) aphi16, aphi5, bvf2, wfac,
176  & cfac, conq, cont, conw,
177  & dk, dkmax, dkmin,
178  & dq1, dsdz2, dsdzq, dsdzt,
179  & dsdzu, dsdzv,
180  & dsig, dt2, dthe1, dtodsd,
181  & dtodsu, dw2, dw2min, g,
182  & gamcrq, gamcrt, gocp,
183  & gravi, f0,
184  & prnum, prmax, prmin, pfac, crbcon,
185  & qmin, tdzmin, qtend, crbmin,crbmax,
186  & rbint, rdt, rdz, qlmin,
187  & ri, rimin, rl2, rlam, rlamun,
188  & rone, rzero, sfcfrac,
189  & spdk2, sri, zol1, zolcr, zolcru,
190  & robn, ttend,
191  & utend, vk, vk2,
192  & ust3, wst3,
193  & vtend, zfac, vpert, cteit,
194  & rentf1, rentf2, radfac,
195  & zfmin, zk, tem, tem1, tem2,
196  & xkzm, xkzmu, xkzminv,
197  & ptem, ptem1, ptem2, tx1(im), tx2(im)
198 !
199  real(kind=kind_phys) zstblmax,h1, h2, qlcr, actei,
200  & cldtime
201 cc
202  parameter(gravi=1.0/grav)
203  parameter(g=grav)
204  parameter(gocp=g/cp)
205  parameter(cont=cp/g,conq=hvap/g,conw=1.0/g) ! for del in pa
206 ! parameter(cont=1000.*cp/g,conq=1000.*hvap/g,conw=1000./g) ! for del in kpa
207  parameter(rlam=30.0,vk=0.4,vk2=vk*vk)
208  parameter(prmin=0.25,prmax=4.,zolcr=0.2,zolcru=-0.5)
209  parameter(dw2min=0.0001,dkmin=0.0,dkmax=1000.,rimin=-100.)
210  parameter(crbcon=0.25,crbmin=0.15,crbmax=0.35)
211  parameter(wfac=7.0,cfac=6.5,pfac=2.0,sfcfrac=0.1)
212 ! parameter(qmin=1.e-8,xkzm=1.0,zfmin=1.e-8,aphi5=5.,aphi16=16.)
213  parameter(qmin=1.e-8, zfmin=1.e-8,aphi5=5.,aphi16=16.)
214  parameter(tdzmin=1.e-3,qlmin=1.e-12,f0=1.e-4)
215  parameter(h1=0.33333333,h2=0.66666667)
216  parameter(cldtime=500.,xkzminv=0.3)
217 ! parameter(cldtime=500.,xkzmu=3.0,xkzminv=0.3)
218 ! parameter(gamcrt=3.,gamcrq=2.e-3,rlamun=150.0)
219  parameter(gamcrt=3.,gamcrq=0.,rlamun=150.0)
220  parameter(rentf1=0.2,rentf2=1.0,radfac=0.85)
221  parameter(iun=84)
222 !
223 ! parameter (zstblmax = 2500., qlcr=1.0e-5)
224 ! parameter (zstblmax = 2500., qlcr=3.0e-5)
225 ! parameter (zstblmax = 2500., qlcr=3.5e-5)
226 ! parameter (zstblmax = 2500., qlcr=1.0e-4)
227  parameter(zstblmax = 2500., qlcr=3.5e-5)
228 ! parameter (actei = 0.23)
229  parameter(actei = 0.7)
230 c
231 c-----------------------------------------------------------------------
232 c
233  601 format(1x,' moninp lat lon step hour ',3i6,f6.1)
234  602 format(1x,' k',' z',' t',' th',
235  1 ' tvh',' q',' u',' v',
236  2 ' sp')
237  603 format(1x,i5,8f9.1)
238  604 format(1x,' sfc',9x,f9.1,18x,f9.1)
239  605 format(1x,' k zl spd2 thekv the1v'
240  1 ,' thermal rbup')
241  606 format(1x,i5,6f8.2)
242  607 format(1x,' kpbl hpbl fm fh hgamt',
243  1 ' hgamq ws ustar cd ch')
244  608 format(1x,i5,9f8.2)
245  609 format(1x,' k pr dkt dku ',i5,3f8.2)
246  610 format(1x,' k pr dkt dku ',i5,3f8.2,' l2 ri t2',
247  1 ' sr2 ',2f8.2,2e10.2)
248 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
250 
251 ! compute preliminary variables
252 !
253  if (ix .lt. im) stop
254 !
255 ! iprt = 0
256 ! if(iprt.eq.1) then
257 !cc latd = 0
258 ! lond = 0
259 ! else
260 !cc latd = 0
261 ! lond = 0
262 ! endif
263 !
264  dt2 = delt
265  rdt = 1. / dt2
266  km1 = km - 1
267  kmpbl = km / 2
269  do k=1,km
270  do i=1,im
271  zi(i,k) = phii(i,k) * gravi
272  zl(i,k) = phil(i,k) * gravi
273  enddo
274  enddo
275  do i=1,im
276  zi(i,km+1) = phii(i,km+1) * gravi
277  enddo
279  do k = 1,km1
280  do i=1,im
281  rdzt(i,k) = 1.0 / (zl(i,k+1) - zl(i,k))
282  enddo
283  enddo
285  do i=1,im
286  kx1(i) = 1
287  tx1(i) = 1.0 / prsi(i,1)
288  tx2(i) = tx1(i)
289  enddo
291  do k = 1,km1
292  do i=1,im
293  xkzo(i,k) = 0.0
294  xkzmo(i,k) = 0.0
295  if (k < kinver(i)) then
296 ! vertical background diffusivity
297  ptem = prsi(i,k+1) * tx1(i)
298  tem1 = 1.0 - ptem
299  tem1 = tem1 * tem1 * 10.0
300  xkzo(i,k) = xkzm_h * min(1.0, exp(-tem1))
301 
302 ! vertical background diffusivity for momentum
303  if (ptem >= xkzm_s) then
304  xkzmo(i,k) = xkzm_m
305  kx1(i) = k + 1
306  else
307  if (k == kx1(i) .and. k > 1) tx2(i) = 1.0 / prsi(i,k)
308  tem1 = 1.0 - prsi(i,k+1) * tx2(i)
309  tem1 = tem1 * tem1 * 5.0
310  xkzmo(i,k) = xkzm_m * min(1.0, exp(-tem1))
311  endif
312  endif
313  enddo
314  enddo
315 ! if (lprnt) then
316 ! print *,' xkzo=',(xkzo(ipr,k),k=1,km1)
317 ! print *,' xkzmo=',(xkzmo(ipr,k),k=1,km1)
318 ! endif
319 !
320 ! diffusivity in the inversion layer is set to be xkzminv (m^2/s)
322  do k = 1,kmpbl
323  do i=1,im
324 ! if(zi(i,k+1) > 200..and.zi(i,k+1) < zstblmax) then
325  if(zi(i,k+1) > 250.) then
326  tem1 = (t1(i,k+1)-t1(i,k)) * rdzt(i,k)
327  if(tem1 > 1.e-5) then
328  xkzo(i,k) = min(xkzo(i,k),xkzminv)
329  endif
330  endif
331  enddo
332  enddo
334  do i = 1,im
335  z0(i) = 0.01 * zorl(i)
336  dusfc(i) = 0.
337  dvsfc(i) = 0.
338  dtsfc(i) = 0.
339  dqsfc(i) = 0.
340  wscale(i)= 0.
341  wscaleu(i)= 0.
342  kpbl(i) = 1
343  hpbl(i) = zi(i,1)
344  hpblx(i) = zi(i,1)
345  pblflg(i)= .true.
346  sfcflg(i)= .true.
347  if(rbsoil(i) > 0.) sfcflg(i) = .false.
348  ublflg(i)= .false.
349  pcnvflg(i)= .false.
350  scuflg(i)= .true.
351  if(scuflg(i)) then
352  radmin(i)= 0.
353  rent(i) = rentf1
354  hrad(i) = zi(i,1)
355 ! hradm(i) = zi(i,1)
356  krad(i) = 1
357  icld(i) = 0
358  lcld(i) = km1
359  kcld(i) = km1
360  zd(i) = 0.
361  endif
362  enddo
364  do k = 1,km
365  do i = 1,im
366  theta(i,k) = t1(i,k) * psk(i) / prslk(i,k)
367  qlx(i,k) = max(q1(i,k,ntcw),qlmin)
368  qtx(i,k) = max(q1(i,k,1),qmin)+qlx(i,k)
369  ptem = qlx(i,k)
370  ptem1 = hvap*max(q1(i,k,1),qmin)/(cp*t1(i,k))
371  thetae(i,k)= theta(i,k)*(1.+ptem1)
372  thvx(i,k) = theta(i,k)*(1.+fv*max(q1(i,k,1),qmin)-ptem)
373  ptem2 = theta(i,k)-(hvap/cp)*ptem
374  thlvx(i,k) = ptem2*(1.+fv*qtx(i,k))
375  enddo
376  enddo
378  do k = 1,km1
379  do i = 1,im
380  dku(i,k) = 0.
381  dkt(i,k) = 0.
382  dktx(i,k) = 0.
383  cku(i,k) = 0.
384  ckt(i,k) = 0.
385  tem = zi(i,k+1)-zi(i,k)
386  radx(i,k) = tem*(swh(i,k)*xmu(i)+hlw(i,k))
387  enddo
388  enddo
390  do i=1,im
391  flg(i) = scuflg(i)
392  enddo
393  do k = 1, km1
394  do i=1,im
395  if(flg(i).and.zl(i,k) >= zstblmax) then
396  lcld(i)=k
397  flg(i)=.false.
398  endif
399  enddo
400  enddo
401 !
402 ! compute virtual potential temp gradient (bf) and winshear square
404  do k = 1, km1
405  do i = 1, im
406  rdz = rdzt(i,k)
407  bf(i,k) = (thvx(i,k+1)-thvx(i,k))*rdz
408  ti(i,k) = 2./(t1(i,k)+t1(i,k+1))
409  dw2 = (u1(i,k)-u1(i,k+1))**2
410  & + (v1(i,k)-v1(i,k+1))**2
411  shr2(i,k) = max(dw2,dw2min)*rdz*rdz
412  enddo
413  enddo
415  do i = 1,im
416  govrth(i) = g/theta(i,1)
417  enddo
418 !
419  do i=1,im
420  beta(i) = dt2 / (zi(i,2)-zi(i,1))
421  enddo
422 !
423  do i=1,im
424  ustar(i) = sqrt(stress(i))
425  enddo
426 !
427  do i = 1,im
428  sflux(i) = heat(i) + evap(i)*fv*theta(i,1)
429  if(.not.sfcflg(i) .or. sflux(i) <= 0.) pblflg(i)=.false.
430  enddo
435 ! compute the pbl height
436 !
437  do i=1,im
438  flg(i) = .false.
439  rbup(i) = rbsoil(i)
440 !
441  if(pblflg(i)) then
442  thermal(i) = thvx(i,1)
443  crb(i) = crbcon
444  else
445  thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin))
446  tem = sqrt(u10m(i)**2+v10m(i)**2)
447  tem = max(tem, 1.)
448  robn = tem / (f0 * z0(i))
449  tem1 = 1.e-7 * robn
450  crb(i) = 0.16 * (tem1 ** (-0.18))
451  crb(i) = max(min(crb(i), crbmax), crbmin)
452  endif
453  enddo
462  do k = 1, kmpbl
463  do i = 1, im
464  if(.not.flg(i)) then
465  rbdn(i) = rbup(i)
466  spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
467  rbup(i) = (thvx(i,k)-thermal(i))*
468  & (g*zl(i,k)/thvx(i,1))/spdk2
469  kpbl(i) = k
470  flg(i) = rbup(i) > crb(i)
471  endif
472  enddo
473  enddo
475  do i = 1,im
476  if(kpbl(i) > 1) then
477  k = kpbl(i)
478  if(rbdn(i) >= crb(i)) then
479  rbint = 0.
480  elseif(rbup(i) <= crb(i)) then
481  rbint = 1.
482  else
483  rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
484  endif
485  hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
486  if(hpbl(i) < zi(i,kpbl(i))) kpbl(i) = kpbl(i) - 1
487  else
488  hpbl(i) = zl(i,1)
489  kpbl(i) = 1
490  endif
491  kpblx(i) = kpbl(i)
492  hpblx(i) = hpbl(i)
493  enddo
494 !
495 ! compute similarity parameters
511  do i=1,im
512  zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin)
513  if(sfcflg(i)) then
514  zol(i) = min(zol(i),-zfmin)
515  else
516  zol(i) = max(zol(i),zfmin)
517  endif
518  zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1)
519  if(sfcflg(i)) then
520 ! phim(i) = (1.-aphi16*zol1)**(-1./4.)
521 ! phih(i) = (1.-aphi16*zol1)**(-1./2.)
522  tem = 1.0 / (1. - aphi16*zol1)
523  phih(i) = sqrt(tem)
524  phim(i) = sqrt(phih(i))
525  else
526  phim(i) = 1. + aphi5*zol1
527  phih(i) = phim(i)
528  endif
529  wscale(i) = ustar(i)/phim(i)
530  ustmin(i) = ustar(i)/aphi5
531  wscale(i) = max(wscale(i),ustmin(i))
532  enddo
533  do i=1,im
534  if(pblflg(i)) then
535  if(zol(i) < zolcru .and. kpbl(i) > 1) then
536  pcnvflg(i) = .true.
537  else
538  ublflg(i) = .true.
539  endif
540  wst3 = govrth(i)*sflux(i)*hpbl(i)
541  wstar(i)= wst3**h1
542  ust3 = ustar(i)**3.
543  wscaleu(i) = (ust3+wfac*vk*wst3*sfcfrac)**h1
544  wscaleu(i) = max(wscaleu(i),ustmin(i))
545  endif
546  enddo
547 !
548 ! compute counter-gradient mixing term for heat and moisture
551  do i = 1,im
552  if(ublflg(i)) then
553  hgamt(i) = min(cfac*heat(i)/wscaleu(i),gamcrt)
554  hgamq(i) = min(cfac*evap(i)/wscaleu(i),gamcrq)
555  vpert = hgamt(i) + hgamq(i)*fv*theta(i,1)
556  vpert = min(vpert,gamcrt)
557  thermal(i)= thermal(i)+max(vpert,0.)
558  hgamt(i) = max(hgamt(i),0.0)
559  hgamq(i) = max(hgamq(i),0.0)
560  endif
561  enddo
562 !
563 ! enhance the pbl height by considering the thermal excess
565  do i=1,im
566  flg(i) = .true.
567  if(ublflg(i)) then
568  flg(i) = .false.
569  rbup(i) = rbsoil(i)
570  endif
571  enddo
572  do k = 2, kmpbl
573  do i = 1, im
574  if(.not.flg(i)) then
575  rbdn(i) = rbup(i)
576  spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
577  rbup(i) = (thvx(i,k)-thermal(i))*
578  & (g*zl(i,k)/thvx(i,1))/spdk2
579  kpbl(i) = k
580  flg(i) = rbup(i) > crb(i)
581  endif
582  enddo
583  enddo
584  do i = 1,im
585  if(ublflg(i)) then
586  k = kpbl(i)
587  if(rbdn(i) >= crb(i)) then
588  rbint = 0.
589  elseif(rbup(i) <= crb(i)) then
590  rbint = 1.
591  else
592  rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
593  endif
594  hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
595  if(hpbl(i) < zi(i,kpbl(i))) kpbl(i) = kpbl(i) - 1
596  if(kpbl(i) <= 1) then
597  ublflg(i) = .false.
598  pblflg(i) = .false.
599  endif
600  endif
601  enddo
602 !
603 ! look for stratocumulus
606  do i = 1, im
607  flg(i)=scuflg(i)
608  enddo
609  do k = kmpbl,1,-1
610  do i = 1, im
611  if(flg(i) .and. k <= lcld(i)) then
612  if(qlx(i,k).ge.qlcr) then
613  kcld(i)=k
614  flg(i)=.false.
615  endif
616  endif
617  enddo
618  enddo
619  do i = 1, im
620  if(scuflg(i) .and. kcld(i)==km1) scuflg(i)=.false.
621  enddo
623  do i = 1, im
624  flg(i)=scuflg(i)
625  enddo
626  do k = kmpbl,1,-1
627  do i = 1, im
628  if(flg(i) .and. k <= kcld(i)) then
629  if(qlx(i,k) >= qlcr) then
630  if(radx(i,k) < radmin(i)) then
631  radmin(i)=radx(i,k)
632  krad(i)=k
633  endif
634  else
635  flg(i)=.false.
636  endif
637  endif
638  enddo
639  enddo
640  do i = 1, im
641  if(scuflg(i) .and. krad(i) <= 1) scuflg(i)=.false.
642  if(scuflg(i) .and. radmin(i)>=0.) scuflg(i)=.false.
643  enddo
645  do i = 1, im
646  flg(i)=scuflg(i)
647  enddo
648  do k = kmpbl,2,-1
649  do i = 1, im
650  if(flg(i) .and. k <= krad(i)) then
651  if(qlx(i,k) >= qlcr) then
652  icld(i)=icld(i)+1
653  else
654  flg(i)=.false.
655  endif
656  endif
657  enddo
658  enddo
659  do i = 1, im
660  if(scuflg(i) .and. icld(i) < 1) scuflg(i)=.false.
661  enddo
663  do i = 1, im
664  if(scuflg(i)) then
665  hrad(i) = zi(i,krad(i)+1)
666 ! hradm(i)= zl(i,krad(i))
667  endif
668  enddo
669 !
670  do i = 1, im
671  if(scuflg(i) .and. hrad(i)<zi(i,2)) scuflg(i)=.false.
672  enddo
674  do i = 1, im
675  if(scuflg(i)) then
676  k = krad(i)
677  tem = zi(i,k+1)-zi(i,k)
678  tem1 = cldtime*radmin(i)/tem
679  thlvx1(i) = thlvx(i,k)+tem1
680 ! if(thlvx1(i) > thlvx(i,k-1)) scuflg(i)=.false.
681  endif
682  enddo
684  do i = 1, im
685  flg(i)=scuflg(i)
686  enddo
687  do k = kmpbl,1,-1
688  do i = 1, im
689  if(flg(i) .and. k <= krad(i))then
690  if(thlvx1(i) <= thlvx(i,k))then
691  tem=zi(i,k+1)-zi(i,k)
692  zd(i)=zd(i)+tem
693  else
694  flg(i)=.false.
695  endif
696  endif
697  enddo
698  enddo
700  do i = 1, im
701  if(scuflg(i))then
702  kk = max(1, krad(i)+1-icld(i))
703  zdd(i) = hrad(i)-zi(i,kk)
704  endif
705  enddo
707  do i = 1, im
708  if(scuflg(i))then
709  zd(i) = max(zd(i),zdd(i))
710  zd(i) = min(zd(i),hrad(i))
711  tem = govrth(i)*zd(i)*(-radmin(i))
712  vrad(i)= tem**h1
713  endif
714  enddo
715 !
716 ! compute inverse prandtl number
719  do i = 1, im
720  if(ublflg(i)) then
721  tem = phih(i)/phim(i)+cfac*vk*sfcfrac
722  else
723  tem = phih(i)/phim(i)
724  endif
725  prinv(i) = 1.0 / tem
726  prinv(i) = min(prinv(i),prmax)
727  prinv(i) = max(prinv(i),prmin)
728  enddo
729  do i = 1, im
730  if(zol(i) > zolcr) then
731  kpbl(i) = 1
732  endif
733  enddo
734 !
735 ! compute diffusion coefficients below pbl
738  do k = 1, kmpbl
739  do i=1,im
740  if(k < kpbl(i)) then
741 ! zfac = max((1.-(zi(i,k+1)-zl(i,1))/
742 ! 1 (hpbl(i)-zl(i,1))), zfmin)
743  zfac = max((1.-zi(i,k+1)/hpbl(i)), zfmin)
744  tem = zi(i,k+1) * (zfac**pfac)
745  if(pblflg(i)) then
746  tem1 = vk * wscaleu(i) * tem
747 ! dku(i,k) = xkzmo(i,k) + tem1
748 ! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i)
749  dku(i,k) = tem1
750  dkt(i,k) = tem1 * prinv(i)
751  else
752  tem1 = vk * wscale(i) * tem
753 ! dku(i,k) = xkzmo(i,k) + tem1
754 ! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i)
755  dku(i,k) = tem1
756  dkt(i,k) = tem1 * prinv(i)
757  endif
758  dku(i,k) = min(dku(i,k),dkmax)
759  dku(i,k) = max(dku(i,k),xkzmo(i,k))
760  dkt(i,k) = min(dkt(i,k),dkmax)
761  dkt(i,k) = max(dkt(i,k),xkzo(i,k))
762  dktx(i,k)= dkt(i,k)
763  endif
764  enddo
765  enddo
766 !
767 ! compute diffusion coefficients based on local scheme above pbl
802  do k = 1, km1
803  do i=1,im
804  if(k >= kpbl(i)) then
805  bvf2 = g*bf(i,k)*ti(i,k)
806  ri = max(bvf2/shr2(i,k),rimin)
807  zk = vk*zi(i,k+1)
808  if(ri < 0.) then ! unstable regime
809  rl2 = zk*rlamun/(rlamun+zk)
810  dk = rl2*rl2*sqrt(shr2(i,k))
811  sri = sqrt(-ri)
812 ! dku(i,k) = xkzmo(i,k) + dk*(1+8.*(-ri)/(1+1.746*sri))
813 ! dkt(i,k) = xkzo(i,k) + dk*(1+8.*(-ri)/(1+1.286*sri))
814  dku(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri))
815  dkt(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri))
816  else ! stable regime
817  rl2 = zk*rlam/(rlam+zk)
818 !! tem = rlam * sqrt(0.01*prsi(i,k))
819 !! rl2 = zk*tem/(tem+zk)
820  dk = rl2*rl2*sqrt(shr2(i,k))
821  tem1 = dk/(1+5.*ri)**2
822 !
823  if(k >= kpblx(i)) then
824  prnum = 1.0 + 2.1*ri
825  prnum = min(prnum,prmax)
826  else
827  prnum = 1.0
828  endif
829 ! dku(i,k) = xkzmo(i,k) + tem1 * prnum
830 ! dkt(i,k) = xkzo(i,k) + tem1
831  dku(i,k) = tem1 * prnum
832  dkt(i,k) = tem1
833  endif
834 !
835  dku(i,k) = min(dku(i,k),dkmax)
836  dku(i,k) = max(dku(i,k),xkzmo(i,k))
837  dkt(i,k) = min(dkt(i,k),dkmax)
838  dkt(i,k) = max(dkt(i,k),xkzo(i,k))
839 !
840  endif
841 !
842  enddo
843  enddo
844 !
845 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
846 ! compute components for mass flux mixing by large thermals
849  do k = 1, km
850  do i = 1, im
851  if(pcnvflg(i)) then
852  tcko(i,k) = t1(i,k)
853  ucko(i,k) = u1(i,k)
854  vcko(i,k) = v1(i,k)
855  xmf(i,k) = 0.
856  endif
857  enddo
858  enddo
859  do kk = 1, ntrac
860  do k = 1, km
861  do i = 1, im
862  if(pcnvflg(i)) then
863  qcko(i,k,kk) = q1(i,k,kk)
864  endif
865  enddo
866  enddo
867  enddo
869  call mfpbl(im,ix,km,ntrac,dt2,pcnvflg,
870  & zl,zi,thvx,q1,t1,u1,v1,hpbl,kpbl,
871  & sflux,ustar,wstar,xmf,tcko,qcko,ucko,vcko)
872 !
873 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
874 ! compute diffusion coefficients for cloud-top driven diffusion
875 ! if the condition for cloud-top instability is met,
876 ! increase entrainment flux at cloud top
877 !
888  do i = 1, im
889  if(scuflg(i)) then
890  k = krad(i)
891  tem = thetae(i,k) - thetae(i,k+1)
892  tem1 = qtx(i,k) - qtx(i,k+1)
893  if (tem > 0. .and. tem1 > 0.) then
894  cteit= cp*tem/(hvap*tem1)
895  if(cteit > actei) rent(i) = rentf2
896  endif
897  endif
898  enddo
899  do i = 1, im
900  if(scuflg(i)) then
901  k = krad(i)
902  tem1 = max(bf(i,k),tdzmin)
903  ckt(i,k) = -rent(i)*radmin(i)/tem1
904  cku(i,k) = ckt(i,k)
905  endif
906  enddo
907 !
908  do k = 1, kmpbl
909  do i=1,im
910  if(scuflg(i) .and. k < krad(i)) then
911  tem1=hrad(i)-zd(i)
912  tem2=zi(i,k+1)-tem1
913  if(tem2 > 0.) then
914  ptem= tem2/zd(i)
915  if(ptem.ge.1.) ptem= 1.
916  ptem= tem2*ptem*sqrt(1.-ptem)
917  ckt(i,k) = radfac*vk*vrad(i)*ptem
918  cku(i,k) = 0.75*ckt(i,k)
919  ckt(i,k) = max(ckt(i,k),dkmin)
920  ckt(i,k) = min(ckt(i,k),dkmax)
921  cku(i,k) = max(cku(i,k),dkmin)
922  cku(i,k) = min(cku(i,k),dkmax)
923  endif
924  endif
925  enddo
926  enddo
927 !
928 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
929 !
931  do k = 1, kmpbl
932  do i=1,im
933  if(scuflg(i)) then
934  dkt(i,k) = dkt(i,k)+ckt(i,k)
935  dku(i,k) = dku(i,k)+cku(i,k)
936  dkt(i,k) = min(dkt(i,k),dkmax)
937  dku(i,k) = min(dku(i,k),dkmax)
938  endif
939  enddo
940  enddo
941 !
942 ! compute tridiagonal matrix elements for heat and moisture
943 !
946  do i=1,im
947  ad(i,1) = 1.
948  a1(i,1) = t1(i,1) + beta(i) * heat(i)
949  a2(i,1) = q1(i,1,1) + beta(i) * evap(i)
950  enddo
951 
952  if(ntrac >= 2) then
953  do k = 2, ntrac
954  is = (k-1) * km
955  do i = 1, im
956  a2(i,1+is) = q1(i,1,k)
957  enddo
958  enddo
959  endif
960 !
961  do k = 1,km1
962  do i = 1,im
963  dtodsd = dt2/del(i,k)
964  dtodsu = dt2/del(i,k+1)
965  dsig = prsl(i,k)-prsl(i,k+1)
966  rdz = rdzt(i,k)
967  tem1 = dsig * dkt(i,k) * rdz
968  dsdz2 = tem1 * rdz
969  au(i,k) = -dtodsd*dsdz2
970  al(i,k) = -dtodsu*dsdz2
971 !
972  if(pcnvflg(i) .and. k < kpbl(i)) then
973  tem2 = dsig * rdz
974  ptem = 0.5 * tem2 * xmf(i,k)
975  ptem1 = dtodsd * ptem
976  ptem2 = dtodsu * ptem
977  ad(i,k) = ad(i,k)-au(i,k)-ptem1
978  ad(i,k+1) = 1.-al(i,k)+ptem2
979  au(i,k) = au(i,k)-ptem1
980  al(i,k) = al(i,k)+ptem2
981  ptem = tcko(i,k) + tcko(i,k+1)
982  dsdzt = tem1 * gocp
983  a1(i,k) = a1(i,k)+dtodsd*dsdzt-ptem1*ptem
984  a1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt+ptem2*ptem
985  ptem = qcko(i,k,1) + qcko(i,k+1,1)
986  a2(i,k) = a2(i,k) - ptem1 * ptem
987  a2(i,k+1) = q1(i,k+1,1) + ptem2 * ptem
988  elseif(ublflg(i) .and. k < kpbl(i)) then
989  ptem1 = dsig * dktx(i,k) * rdz
990  tem = 1.0 / hpbl(i)
991  dsdzt = tem1 * gocp - ptem1 * hgamt(i) * tem
992  dsdzq = - ptem1 * hgamq(i) * tem
993  ad(i,k) = ad(i,k)-au(i,k)
994  ad(i,k+1) = 1.-al(i,k)
995  a1(i,k) = a1(i,k)+dtodsd*dsdzt
996  a1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt
997  a2(i,k) = a2(i,k)+dtodsd*dsdzq
998  a2(i,k+1) = q1(i,k+1,1)-dtodsu*dsdzq
999  else
1000  ad(i,k) = ad(i,k)-au(i,k)
1001  ad(i,k+1) = 1.-al(i,k)
1002  dsdzt = tem1 * gocp
1003  a1(i,k) = a1(i,k)+dtodsd*dsdzt
1004  a1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt
1005  a2(i,k+1) = q1(i,k+1,1)
1006  endif
1007 !
1008  enddo
1009  enddo
1010 !
1011  if(ntrac >= 2) then
1012  do kk = 2, ntrac
1013  is = (kk-1) * km
1014  do k = 1, km1
1015  do i = 1, im
1016  if(pcnvflg(i) .and. k < kpbl(i)) then
1017  dtodsd = dt2/del(i,k)
1018  dtodsu = dt2/del(i,k+1)
1019  dsig = prsl(i,k)-prsl(i,k+1)
1020  tem = dsig * rdzt(i,k)
1021  ptem = 0.5 * tem * xmf(i,k)
1022  ptem1 = dtodsd * ptem
1023  ptem2 = dtodsu * ptem
1024  tem1 = qcko(i,k,kk) + qcko(i,k+1,kk)
1025  a2(i,k+is) = a2(i,k+is) - ptem1*tem1
1026  a2(i,k+1+is)= q1(i,k+1,kk) + ptem2*tem1
1027  else
1028  a2(i,k+1+is) = q1(i,k+1,kk)
1029  endif
1030  enddo
1031  enddo
1032  enddo
1033  endif
1034 !
1035 ! solve tridiagonal problem for heat and moisture
1036 !
1038  call tridin(im,km,ntrac,al,ad,au,a1,a2,au,a1,a2)
1039 
1040 !
1041 ! recover tendencies of heat and moisture
1042 !
1044  do k = 1,km
1045  do i = 1,im
1046  ttend = (a1(i,k)-t1(i,k)) * rdt
1047  qtend = (a2(i,k)-q1(i,k,1))*rdt
1048  tau(i,k) = tau(i,k)+ttend
1049  rtg(i,k,1) = rtg(i,k,1)+qtend
1050  dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend
1051  dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend
1052  enddo
1053  enddo
1054  if(ntrac >= 2) then
1055  do kk = 2, ntrac
1056  is = (kk-1) * km
1057  do k = 1, km
1058  do i = 1, im
1059  qtend = (a2(i,k+is)-q1(i,k,kk))*rdt
1060  rtg(i,k,kk) = rtg(i,k,kk)+qtend
1061  enddo
1062  enddo
1063  enddo
1064  endif
1065 !
1066 ! compute tke dissipation rate
1067 !
1070  if(dspheat) then
1071 !
1072  do k = 1,km1
1073  do i = 1,im
1074  diss(i,k) = dku(i,k)*shr2(i,k)-g*ti(i,k)*dkt(i,k)*bf(i,k)
1075 ! diss(i,k) = dku(i,k)*shr2(i,k)
1076  enddo
1077  enddo
1078 !
1079 ! add dissipative heating at the first model layer
1080 !
1082  do i = 1,im
1083  tem = govrth(i)*sflux(i)
1084  tem1 = tem + stress(i)*spd1(i)/zl(i,1)
1085  tem2 = 0.5 * (tem1+diss(i,1))
1086  tem2 = max(tem2, 0.)
1087  ttend = tem2 / cp
1088  tau(i,1) = tau(i,1)+0.5*ttend
1089  enddo
1090 !
1091 ! add dissipative heating above the first model layer
1092 !
1093  do k = 2,km1
1094  do i = 1,im
1095  tem = 0.5 * (diss(i,k-1)+diss(i,k))
1096  tem = max(tem, 0.)
1097  ttend = tem / cp
1098  tau(i,k) = tau(i,k) + 0.5*ttend
1099  enddo
1100  enddo
1101 !
1102  endif
1103 !
1104 ! compute tridiagonal matrix elements for momentum
1105 !
1108  do i=1,im
1109  ad(i,1) = 1.0 + beta(i) * stress(i) / spd1(i)
1110  a1(i,1) = u1(i,1)
1111  a2(i,1) = v1(i,1)
1112  enddo
1113 !
1114  do k = 1,km1
1115  do i=1,im
1116  dtodsd = dt2/del(i,k)
1117  dtodsu = dt2/del(i,k+1)
1118  dsig = prsl(i,k)-prsl(i,k+1)
1119  rdz = rdzt(i,k)
1120  tem1 = dsig*dku(i,k)*rdz
1121  dsdz2 = tem1 * rdz
1122  au(i,k) = -dtodsd*dsdz2
1123  al(i,k) = -dtodsu*dsdz2
1124 !
1125  if(pcnvflg(i) .and. k < kpbl(i)) then
1126  tem2 = dsig * rdz
1127  ptem = 0.5 * tem2 * xmf(i,k)
1128  ptem1 = dtodsd * ptem
1129  ptem2 = dtodsu * ptem
1130  ad(i,k) = ad(i,k)-au(i,k)-ptem1
1131  ad(i,k+1) = 1.-al(i,k)+ptem2
1132  au(i,k) = au(i,k)-ptem1
1133  al(i,k) = al(i,k)+ptem2
1134  ptem = ucko(i,k) + ucko(i,k+1)
1135  a1(i,k) = a1(i,k) - ptem1 * ptem
1136  a1(i,k+1) = u1(i,k+1) + ptem2 * ptem
1137  ptem = vcko(i,k) + vcko(i,k+1)
1138  a2(i,k) = a2(i,k) - ptem1 * ptem
1139  a2(i,k+1) = v1(i,k+1) + ptem2 * ptem
1140  else
1141  ad(i,k) = ad(i,k)-au(i,k)
1142  ad(i,k+1) = 1.-al(i,k)
1143  a1(i,k+1) = u1(i,k+1)
1144  a2(i,k+1) = v1(i,k+1)
1145  endif
1146 !
1147  enddo
1148  enddo
1149 !
1150 ! solve tridiagonal problem for momentum
1151 !
1152  call tridi2(im,km,al,ad,au,a1,a2,au,a1,a2)
1153 !
1154 ! recover tendencies of momentum
1155 !
1157  do k = 1,km
1158  do i = 1,im
1159  utend = (a1(i,k)-u1(i,k))*rdt
1160  vtend = (a2(i,k)-v1(i,k))*rdt
1161  du(i,k) = du(i,k) + utend
1162  dv(i,k) = dv(i,k) + vtend
1163  dusfc(i) = dusfc(i) + conw*del(i,k)*utend
1164  dvsfc(i) = dvsfc(i) + conw*del(i,k)*vtend
1165 !
1166 ! for dissipative heating for ecmwf model
1167 !
1168 ! tem1 = 0.5*(a1(i,k)+u1(i,k))
1169 ! tem2 = 0.5*(a2(i,k)+v1(i,k))
1170 ! diss(i,k) = -(tem1*utend+tem2*vtend)
1171 ! diss(i,k) = max(diss(i,k),0.)
1172 ! ttend = diss(i,k) / cp
1173 ! tau(i,k) = tau(i,k) + ttend
1174 !
1175  enddo
1176  enddo
1177 !
1178 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1179 !
1180  do i = 1, im
1181  hpbl(i) = hpblx(i)
1182  kpbl(i) = kpblx(i)
1183  enddo
1184 !
1185 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1186  return
1187  end
1188 
1189 
1190 c-----------------------------------------------------------------------
1195  subroutine tridi2(l,n,cl,cm,cu,r1,r2,au,a1,a2)
1196 cc
1197  use machine , only : kind_phys
1198  implicit none
1199  integer k,n,l,i
1200  real(kind=kind_phys) fk
1201 cc
1202  real(kind=kind_phys) cl(l,2:n),cm(l,n),cu(l,n-1),r1(l,n),r2(l,n), &
1203  & au(l,n-1),a1(l,n),a2(l,n)
1204 c-----------------------------------------------------------------------
1205  do i=1,l
1206  fk = 1./cm(i,1)
1207  au(i,1) = fk*cu(i,1)
1208  a1(i,1) = fk*r1(i,1)
1209  a2(i,1) = fk*r2(i,1)
1210  enddo
1211  do k=2,n-1
1212  do i=1,l
1213  fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1214  au(i,k) = fk*cu(i,k)
1215  a1(i,k) = fk*(r1(i,k)-cl(i,k)*a1(i,k-1))
1216  a2(i,k) = fk*(r2(i,k)-cl(i,k)*a2(i,k-1))
1217  enddo
1218  enddo
1219  do i=1,l
1220  fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1221  a1(i,n) = fk*(r1(i,n)-cl(i,n)*a1(i,n-1))
1222  a2(i,n) = fk*(r2(i,n)-cl(i,n)*a2(i,n-1))
1223  enddo
1224  do k=n-1,1,-1
1225  do i=1,l
1226  a1(i,k) = a1(i,k)-au(i,k)*a1(i,k+1)
1227  a2(i,k) = a2(i,k)-au(i,k)*a2(i,k+1)
1228  enddo
1229  enddo
1230 c-----------------------------------------------------------------------
1231  return
1232  end
1233 c-----------------------------------------------------------------------
1238  subroutine tridin(l,n,nt,cl,cm,cu,r1,r2,au,a1,a2)
1239 cc
1240  use machine , only : kind_phys
1241  implicit none
1242  integer is,k,kk,n,nt,l,i
1243  real(kind=kind_phys) fk(l)
1244 cc
1245  real(kind=kind_phys) cl(l,2:n), cm(l,n), cu(l,n-1), &
1246  & r1(l,n), r2(l,n*nt), &
1247  & au(l,n-1), a1(l,n), a2(l,n*nt), &
1248  & fkk(l,2:n-1)
1249 c-----------------------------------------------------------------------
1250  do i=1,l
1251  fk(i) = 1./cm(i,1)
1252  au(i,1) = fk(i)*cu(i,1)
1253  a1(i,1) = fk(i)*r1(i,1)
1254  enddo
1255  do k = 1, nt
1256  is = (k-1) * n
1257  do i = 1, l
1258  a2(i,1+is) = fk(i) * r2(i,1+is)
1259  enddo
1260  enddo
1261  do k=2,n-1
1262  do i=1,l
1263  fkk(i,k) = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1264  au(i,k) = fkk(i,k)*cu(i,k)
1265  a1(i,k) = fkk(i,k)*(r1(i,k)-cl(i,k)*a1(i,k-1))
1266  enddo
1267  enddo
1268  do kk = 1, nt
1269  is = (kk-1) * n
1270  do k=2,n-1
1271  do i=1,l
1272  a2(i,k+is) = fkk(i,k)*(r2(i,k+is)-cl(i,k)*a2(i,k+is-1))
1273  enddo
1274  enddo
1275  enddo
1276  do i=1,l
1277  fk(i) = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1278  a1(i,n) = fk(i)*(r1(i,n)-cl(i,n)*a1(i,n-1))
1279  enddo
1280  do k = 1, nt
1281  is = (k-1) * n
1282  do i = 1, l
1283  a2(i,n+is) = fk(i)*(r2(i,n+is)-cl(i,n)*a2(i,n+is-1))
1284  enddo
1285  enddo
1286  do k=n-1,1,-1
1287  do i=1,l
1288  a1(i,k) = a1(i,k) - au(i,k)*a1(i,k+1)
1289  enddo
1290  enddo
1291  do kk = 1, nt
1292  is = (kk-1) * n
1293  do k=n-1,1,-1
1294  do i=1,l
1295  a2(i,k+is) = a2(i,k+is) - au(i,k)*a2(i,k+is+1)
1296  enddo
1297  enddo
1298  enddo
1299 c-----------------------------------------------------------------------
1300  return
1301  end
1302 
physcons::con_g
real(kind=kind_phys), parameter con_g
gravity ( )
Definition: physcons.f:59
tridi2
subroutine tridi2(l, n, cl, cm, cu, r1, r2, au, a1, a2)
Routine to solve the tridiagonal system to calculate temperature and moisture at ; part of two-part p...
Definition: moninedmf.f:1196
moninedmf
subroutine moninedmf(ix, im, km, ntrac, ntcw, dv, du, tau, rtg, u1, v1, t1, q1, swh, hlw, xmu, psk, rbsoil, zorl, u10m, v10m, fm, fh, tsea, qss, heat, evap, stress, spd1, kpbl, prsi, del, prsl, prslk, phii, phil, delt, dspheat, dusfc, dvsfc, dtsfc, dqsfc, hpbl, hgamt, hgamq, dkt, kinver, xkzm_m, xkzm_h, xkzm_s, lprnt, ipr)
This subroutine contains all of logic for the Hybrid EDMF PBL scheme except for the calculation of th...
Definition: moninedmf.f:96
tridin
subroutine tridin(l, n, nt, cl, cm, cu, r1, r2, au, a1, a2)
Routine to solve the tridiagonal system to calculate u- and v-momentum at ; part of two-part process ...
Definition: moninedmf.f:1239
mfpbl
subroutine mfpbl(im, ix, km, ntrac, delt, cnvflg, zl, zm, thvx, q1, t1, u1, v1, hpbl, kpbl, sflx, ustar, wstar, xmf, tcko, qcko, ucko, vcko)
This subroutine is used for calculating the mass flux and updraft properties.
Definition: mfpbl.f:41
physcons::con_cp
real(kind=kind_phys), parameter con_cp
spec heat air at p ( )
Definition: physcons.f:80
physcons::con_hvap
real(kind=kind_phys), parameter con_hvap
lat heat H2O cond ( )
Definition: physcons.f:90
physcons::con_rd
real(kind=kind_phys), parameter con_rd
gas constant air ( )
Definition: physcons.f:76