62 subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, &
63 & grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, &
64 & dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea, &
65 & psk,rbsoil,zorl,u10m,v10m,fm,fh, &
66 & tsea,heat,evap,stress,spd1,kpbl, &
67 & prsi,del,prsl,prslk,phii,phil,delt, &
68 & dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl, &
69 & kinver,xkzm_m,xkzm_h,xkzm_s, &
70 & index_of_temperature,index_of_x_wind,index_of_y_wind, &
71 & index_of_process_pbl,ntqv,ntoz,dtend,dtidx, &
72 & gen_tend,ldiag3d,errmsg,errflg)
74 use machine ,
only : kind_phys
75 use funcphys ,
only : fpvs
80 integer,
intent(in) :: im, km, ntrac, ntcw, ntiw, ntke
81 integer,
intent(in) :: kinver(:)
82 integer,
intent(out) :: kpbl(:)
84 logical,
intent(in) :: gen_tend, ldiag3d
85 real(kind=kind_phys),
intent(inout),
dimension(:,:,:),
optional ::&
87 integer,
intent(in) :: index_of_temperature,index_of_x_wind, &
88 & index_of_y_wind, ntqv, ntoz, dtidx(:,:), index_of_process_pbl
90 real(kind=kind_phys),
intent(in) :: grav,rd,cp,rv,hvap,hfus,fv, &
92 real(kind=kind_phys),
intent(in) :: delt, xkzm_m, xkzm_h, xkzm_s
93 real(kind=kind_phys),
intent(inout) :: dv(:,:), du(:,:), &
94 & tdt(:,:), rtg(:,:,:)
95 real(kind=kind_phys),
intent(in) :: &
97 & t1(:,:), q1(:,:,:), &
98 & swh(:,:), hlw(:,:), &
100 & psk(:), rbsoil(:), &
101 & zorl(:), tsea(:), &
102 & u10m(:), v10m(:), &
104 & evap(:), heat(:), &
105 & stress(:), spd1(:), &
106 & prsi(:,:), del(:,:), &
107 & prsl(:,:), prslk(:,:), &
108 & phii(:,:), phil(:,:)
109 real(kind=kind_phys),
intent(out) :: &
110 & dusfc(:), dvsfc(:), &
111 & dtsfc(:), dqsfc(:), &
114 logical,
intent(in) :: dspheat
115 character(len=*),
intent(out) :: errmsg
116 integer,
intent(out) :: errflg
124 integer i,is,k,kk,n,km1,kmpbl,kmscu,ntrac1
125 integer lcld(im),kcld(im),krad(im),mrad(im)
126 integer kx1(im), kpblx(im)
128 real(kind=kind_phys) tke(im,km), tkeh(im,km-1)
130 real(kind=kind_phys) theta(im,km),thvx(im,km), thlvx(im,km),
131 & qlx(im,km), thetae(im,km),thlx(im,km),
132 & slx(im,km), svx(im,km), qtx(im,km),
133 & tvx(im,km), pix(im,km), radx(im,km-1),
134 & dku(im,km-1),dkt(im,km-1), dkq(im,km-1),
135 & cku(im,km-1),ckt(im,km-1)
137 real(kind=kind_phys) plyr(im,km), rhly(im,km), cfly(im,km),
140 real(kind=kind_phys) dtdz1(im), gdx(im),
141 & phih(im), phim(im), prn(im,km-1),
142 & rbdn(im), rbup(im), thermal(im),
143 & ustar(im), wstar(im), hpblx(im),
144 & ust3(im), wst3(im),
146 & hgamt(im), hgamq(im),
147 & wscale(im),vpert(im),
148 & zol(im), sflux(im), radj(im),
151 real(kind=kind_phys) radmin(im)
153 real(kind=kind_phys) zi(im,km+1), zl(im,km), zm(im,km),
154 & xkzo(im,km-1),xkzmo(im,km-1),
155 & xkzm_hx(im), xkzm_mx(im),
157 & al(im,km-1), ad(im,km), au(im,km-1),
158 & f1(im,km), f2(im,km*(ntrac-1))
160 real(kind=kind_phys) elm(im,km), ele(im,km), rle(im,km-1),
161 & ckz(im,km), chz(im,km),
162 & diss(im,km-1),prod(im,km-1),
163 & bf(im,km-1), shr2(im,km-1),
164 & xlamue(im,km-1), xlamde(im,km-1),
165 & gotvx(im,km), rlam(im,km-1)
169 real(kind=kind_phys) tcko(im,km), qcko(im,km,ntrac),
170 & ucko(im,km), vcko(im,km),
171 & buou(im,km), xmf(im,km)
175 real(kind=kind_phys) tcdo(im,km), qcdo(im,km,ntrac),
176 & ucdo(im,km), vcdo(im,km),
177 & buod(im,km), xmfd(im,km)
179 logical pblflg(im), sfcflg(im), flg(im)
180 logical scuflg(im), pcnvflg(im)
185 real(kind=kind_phys) aphi16, aphi5,
187 & gamcrt, gamcrq, sfcfrac,
189 & dsdz2, dsdzt, dkmax,
191 & dtodsu, g, factor, dz,
192 & gocp, gravi, zol1, zolcru,
193 & buop, shrp, dtn, cdtn,
194 & prnum, prmax, prmin, prtke,
195 & prscu, dw2, dw2min, zk,
196 & elmfac, elefac, dspmax,
198 & f0, robn, crbmin, crbmax,
199 & es, qs,
value, onemrh,
200 & cfh, gamma, elocp, el2orc,
201 & epsi, beta, chx, cqx,
202 & rdt, rdz, qmin, qlmin,
204 & rbcr, rbint, tdzmin,
206 & ttend, utend, vtend, qtend,
207 & zfac, zfmin, vk, spdk2,
208 & tkmin, xkzinv, dspfac, xkgdx,
211 & ptem, ptem0, ptem1, ptem2
213 real(kind=kind_phys) ck0, ck1, ch0, ch1, ce0, rchck
215 real(kind=kind_phys) qlcr, zstblmax
217 real(kind=kind_phys) h1
220 parameter(wfac=7.0,cfac=4.5)
221 parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1)
222 parameter(vk=0.4,rimin=-100.)
223 parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3)
224 parameter(rlmn=30.,rlmx=500.,elmx=500.)
225 parameter(prmin=0.25,prmax=4.0,prtke=1.0,prscu=0.67)
226 parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35)
227 parameter(tkmin=1.e-9,dspfac=0.5,dspmax=10.0)
228 parameter(qmin=1.e-8,qlmin=1.e-12,zfmin=1.e-8)
229 parameter(aphi5=5.,aphi16=16.)
230 parameter(elmfac=1.0,elefac=1.0,cql=100.)
231 parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=25000.)
232 parameter(qlcr=3.5e-5,zstblmax=2500.,xkzinv=0.15)
233 parameter(h1=0.33333333)
234 parameter(ck0=0.4,ck1=0.15,ch0=0.4,ch1=0.15,ce0=0.4)
235 parameter(rchck=1.5,cdtn=25.)
245 el2orc=hvap*hvap/(rv*cp)
268 zi(i,k) = phii(i,k) * gravi
269 zl(i,k) = phil(i,k) * gravi
279 zi(i,km+1) = phii(i,km+1) * gravi
288 gdx(i) = sqrt(garea(i))
294 tke(i,k) = max(q1(i,k,ntke), tkmin)
299 tkeh(i,k) = 0.5 * (tke(i,k) + tke(i,k+1))
305 rdzt(i,k) = 1.0 / (zl(i,k+1) - zl(i,k))
320 tx1(i) = 1.0 / prsi(i,1)
322 if(gdx(i) >= xkgdx)
then
326 tem = 1. / (xkgdx - 5.)
327 tem1 = (xkzm_h - 0.01) * tem
328 tem2 = (xkzm_m - 0.01) * tem
330 xkzm_hx(i) = 0.01 + tem1 * ptem
331 xkzm_mx(i) = 0.01 + tem2 * ptem
338 if (k < kinver(i))
then
340 ptem = prsi(i,k+1) * tx1(i)
342 tem1 = tem1 * tem1 * 10.0
343 xkzo(i,k) = xkzm_hx(i) * min(1.0, exp(-tem1))
345 if (ptem >= xkzm_s)
then
346 xkzmo(i,k) = xkzm_mx(i)
349 if (k == kx1(i) .and. k > 1) tx2(i) = 1.0 / prsi(i,k)
350 tem1 = 1.0 - prsi(i,k+1) * tx2(i)
351 tem1 = tem1 * tem1 * 5.0
352 xkzmo(i,k) = xkzm_mx(i) * min(1.0, exp(-tem1))
360 z0(i) = 0.01 * zorl(i)
371 if(rbsoil(i) > 0.) sfcflg(i) = .false.
387 pix(i,k) = psk(i) / prslk(i,k)
388 theta(i,k) = t1(i,k) * pix(i,k)
390 tem = max(q1(i,k,ntcw),qlmin)
391 tem1 = max(q1(i,k,ntiw),qlmin)
392 qlx(i,k) = tem + tem1
393 ptem = hvap*tem + (hvap+hfus)*tem1
394 slx(i,k) = cp * t1(i,k) + phil(i,k) - ptem
396 qlx(i,k) = max(q1(i,k,ntcw),qlmin)
397 slx(i,k) = cp * t1(i,k) + phil(i,k) - hvap*qlx(i,k)
399 tem2 = 1.+fv*max(q1(i,k,1),qmin)-qlx(i,k)
400 thvx(i,k) = theta(i,k) * tem2
401 tvx(i,k) = t1(i,k) * tem2
402 qtx(i,k) = max(q1(i,k,1),qmin)+qlx(i,k)
403 thlx(i,k) = theta(i,k) - pix(i,k)*elocp*qlx(i,k)
404 thlvx(i,k) = thlx(i,k) * (1. + fv * qtx(i,k))
405 svx(i,k) = cp * tvx(i,k)
406 ptem1 = elocp * pix(i,k) * max(q1(i,k,1),qmin)
407 thetae(i,k)= theta(i,k) + ptem1
408 gotvx(i,k) = g / tvx(i,k)
417 tem1 = (tvx(i,k+1)-tvx(i,k)) * rdzt(i,k)
418 if(tem1 > 1.e-5)
then
419 xkzo(i,k) = min(xkzo(i,k),xkzinv)
420 xkzmo(i,k) = min(xkzmo(i,k),xkzinv)
429 plyr(i,k) = 0.01 * prsl(i,k)
431 es = 0.01 * fpvs(t1(i,k))
432 qs = max(qmin, eps * es / (plyr(i,k) + epsm1*es))
433 rhly(i,k) = max(0.0, min(1.0, max(qmin, q1(i,k,1))/qs))
441 clwt = 1.0e-6 * (plyr(i,k)*0.001)
442 if (qlx(i,k) > clwt)
then
443 onemrh= max(1.e-10, 1.0-rhly(i,k))
444 tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0)
446 value = max(min( tem1*qlx(i,k), 50.0), 0.0)
447 tem2 = sqrt(sqrt(rhly(i,k)))
448 cfly(i,k) = min(max(tem2*(1.0-exp(-
value)), 0.0), 1.0)
457 tem = 0.5 * (svx(i,k) + svx(i,k+1))
458 tem1 = 0.5 * (t1(i,k) + t1(i,k+1))
459 tem2 = 0.5 * (qstl(i,k) + qstl(i,k+1))
460 cfh = min(cfly(i,k+1),0.5*(cfly(i,k)+cfly(i,k+1)))
462 gamma = el2orc * tem2 / (tem1**2)
464 beta = (1. + gamma*epsi*(1.+fv)) / (1. + gamma)
465 chx = cfh * alp * beta + (1. - cfh) * alp
466 cqx = cfh * alp * hvap * (beta - epsi)
467 cqx = cqx + (1. - cfh) * fv * g
468 ptem1 = (slx(i,k+1)-slx(i,k))*rdzt(i,k)
469 ptem2 = (qtx(i,k+1)-qtx(i,k))*rdzt(i,k)
470 bf(i,k) = chx * ptem1 + cqx * ptem2
485 tem = zi(i,k+1)-zi(i,k)
486 radx(i,k) = tem*(swh(i,k)*xmu(i)+hlw(i,k))
493 sflux(i) = heat(i) + evap(i)*fv*theta(i,1)
494 if(.not.sfcflg(i) .or. sflux(i) <= 0.) pblflg(i)=.false.
515 thermal(i) = thlvx(i,1)
518 thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin))
519 tem = sqrt(u10m(i)**2+v10m(i)**2)
521 robn = tem / (f0 * z0(i))
523 crb(i) = 0.16 * (tem1 ** (-0.18))
524 crb(i) = max(min(crb(i), crbmax), crbmin)
530 dtdz1(i) = dt2 / (zi(i,2)-zi(i,1))
534 ustar(i) = sqrt(stress(i))
544 dw2 = (u1(i,k)-u1(i,k+1))**2
545 & + (v1(i,k)-v1(i,k+1))**2
546 shr2(i,k) = max(dw2,dw2min)*rdz*rdz
567 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
570 rbup(i) = (thlvx(i,k)-thermal(i))*
571 & (g*zl(i,k)/thlvx(i,1))/spdk2
573 flg(i) = rbup(i) > crb(i)
581 if(kpblx(i) > 1)
then
583 if(rbdn(i) >= crb(i))
then
585 elseif(rbup(i) <= crb(i))
then
588 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
590 hpblx(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
591 if(hpblx(i) < zi(i,kpblx(i))) kpblx(i)=kpblx(i)-1
598 if(kpbl(i) <= 1) pblflg(i)=.false.
611 zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin)
613 zol(i) = min(zol(i),-zfmin)
615 zol(i) = max(zol(i),zfmin)
629 zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1)
631 tem = 1.0 / (1. - aphi16*zol1)
633 phim(i) = sqrt(phih(i))
635 phim(i) = 1. + aphi5*zol1
657 if(zol(i) < zolcru)
then
660 wst3(i) = gotvx(i,1)*sflux(i)*hpbl(i)
661 wstar(i)= wst3(i)**h1
662 ust3(i) = ustar(i)**3.
663 wscale(i)=(ust3(i)+wfac*vk*wst3(i)*sfcfrac)**h1
664 ptem = ustar(i)/aphi5
665 wscale(i) = max(wscale(i),ptem)
673 hgamt(i) = heat(i)/wscale(i)
674 hgamq(i) = evap(i)/wscale(i)
675 vpert(i) = hgamt(i) + hgamq(i)*fv*theta(i,1)
676 vpert(i) = max(vpert(i),0.)
677 tem = min(cfac*vpert(i),gamcrt)
678 thermal(i)= thermal(i) + tem
698 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
699 rbup(i) = (thlvx(i,k)-thermal(i))*
700 & (g*zl(i,k)/thlvx(i,1))/spdk2
702 flg(i) = rbup(i) > crb(i)
709 if(rbdn(i) >= crb(i))
then
711 elseif(rbup(i) <= crb(i))
then
714 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
716 hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
717 if(hpbl(i) < zi(i,kpbl(i)))
then
718 kpbl(i) = kpbl(i) - 1
720 if(kpbl(i) <= 1)
then
738 if(flg(i).and.zl(i,k) >= zstblmax)
then
749 if(flg(i) .and. k <= lcld(i))
then
750 if(qlx(i,k) >= qlcr)
then
758 if(scuflg(i) .and. kcld(i)==km1) scuflg(i)=.false.
770 if(flg(i) .and. k <= kcld(i))
then
771 if(qlx(i,k) >= qlcr)
then
772 if(radx(i,k) < radmin(i))
then
783 if(scuflg(i) .and. krad(i) <= 1) scuflg(i)=.false.
784 if(scuflg(i) .and. radmin(i)>=0.) scuflg(i)=.false.
810 qcko(i,k,kk) = q1(i,k,kk)
813 qcdo(i,k,kk) = q1(i,k,kk)
820 call mfpblt(im,im,km,kmpbl,ntcw,ntrac1,dt2,
821 & pcnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx,
822 & gdx,hpbl,kpbl,vpert,buou,xmf,
823 & tcko,qcko,ucko,vcko,xlamue)
826 call mfscu(im,im,km,kmscu,ntcw,ntrac1,dt2,
827 & scuflg,zl,zm,q1,t1,u1,v1,plyr,pix,
828 & thlx,thvx,thlvx,gdx,thetae,radj,
829 & krad,mrad,radmin,buod,xmfd,
830 & tcdo,qcdo,ucdo,vcdo,xlamde)
845 tem = phih(i)/phim(i)
846 ptem = -3.*(max(zi(i,k+1)-sfcfrac*hpbl(i),0.))**2.
849 prn(i,k) = 1. + (tem-1.)*exp(ptem)
853 prn(i,k) = min(prn(i,k),prmax)
854 prn(i,k) = max(prn(i,k),prmin)
856 ckz(i,k) = ck1 + (ck0-ck1)*exp(ptem)
857 ckz(i,k) = min(ckz(i,k),ck0)
858 ckz(i,k) = max(ckz(i,k),ck1)
859 chz(i,k) = ch1 + (ch0-ch1)*exp(ptem)
860 chz(i,k) = min(chz(i,k),ch0)
861 chz(i,k) = max(chz(i,k),ch1)
877 dz = zl(i,n+1) - zl(i,n)
878 ptem = gotvx(i,n)*(thvx(i,n+1)-thvx(i,k))*dz
882 if(bsum >= tke(i,k))
then
884 tem2 = max(ptem, zfmin)
886 tem2 = min(ptem, -zfmin)
888 ptem1 = (bsum - tke(i,k)) / tem2
889 zlup = zlup - ptem1 * dz
902 tem1 = tsea(i)*(1.+fv*max(q1(i,1,1),qmin))
904 dz = zl(i,n) - zl(i,n-1)
908 ptem = gotvx(i,n)*(thvx(i,k)-tem1)*dz
912 if(bsum >= tke(i,k))
then
914 tem2 = max(ptem, zfmin)
916 tem2 = min(ptem, -zfmin)
918 ptem1 = (bsum - tke(i,k)) / tem2
919 zldn = zldn - ptem1 * dz
926 tem = 0.5 * (zi(i,k+1)-zi(i,k))
927 tem1 = min(tem, rlmn)
940 ptem2 = min(zlup,zldn)
941 rlam(i,k) = elmfac * ptem2
942 rlam(i,k) = max(rlam(i,k), tem1)
943 rlam(i,k) = min(rlam(i,k), rlmx)
945 ptem2 = sqrt(zlup*zldn)
946 ele(i,k) = elefac * ptem2
947 ele(i,k) = max(ele(i,k), tem1)
948 ele(i,k) = min(ele(i,k), elmx)
957 if (zol(i) < 0.)
then
958 ptem = 1. - 100. * zol(i)
961 elseif (zol(i) >= 1.)
then
964 ptem = 1. + 2.7 * zol(i)
967 elm(i,k) = zk*rlam(i,k)/(rlam(i,k)+zk)
969 dz = zi(i,k+1) - zi(i,k)
971 elm(i,k) = min(elm(i,k), tem)
972 ele(i,k) = min(ele(i,k), tem)
977 elm(i,km) = elm(i,km1)
978 ele(i,km) = ele(i,km1)
985 tem = 0.5 * (elm(i,k) + elm(i,k+1))
986 tem = tem * sqrt(tkeh(i,k))
989 dku(i,k) = ckz(i,k) * tem
990 dkt(i,k) = dku(i,k) / prn(i,k)
992 dkt(i,k) = chz(i,k) * tem
993 dku(i,k) = dkt(i,k) * prn(i,k)
996 ri = max(bf(i,k)/shr2(i,k),rimin)
999 dkt(i,k) = rchck * dku(i,k)
1001 dkt(i,k) = ch1 * tem
1002 prnum = 1.0 + 2.1*ri
1003 prnum = min(prnum,prmax)
1004 dku(i,k) = dkt(i,k) * prnum
1009 if(k >= mrad(i) .and. k < krad(i))
then
1010 tem1 = ckz(i,k) * tem
1011 ptem1 = tem1 / prscu
1012 dku(i,k) = max(dku(i,k), tem1)
1013 dkt(i,k) = max(dkt(i,k), ptem1)
1017 dkq(i,k) = prtke * dkt(i,k)
1019 dkt(i,k) = min(dkt(i,k),dkmax)
1020 dkt(i,k) = max(dkt(i,k),xkzo(i,k))
1021 dkq(i,k) = min(dkq(i,k),dkmax)
1022 dkq(i,k) = max(dkq(i,k),xkzo(i,k))
1023 dku(i,k) = min(dku(i,k),dkmax)
1024 dku(i,k) = max(dku(i,k),xkzmo(i,k))
1032 tem = bf(i,k) / gotvx(i,k)
1033 tem1 = max(tem, tdzmin)
1034 ptem = radj(i) / tem1
1035 dkt(i,k) = dkt(i,k) + ptem
1036 dku(i,k) = dku(i,k) + ptem
1037 dkq(i,k) = dkq(i,k) + ptem
1047 tem = -dkt(i,1) * bf(i,1)
1053 if(scuflg(i) .and. mrad(i) == 1)
then
1054 ptem2 = xmfd(i,1) * buod(i,1)
1058 tem = tem + ptem1 + ptem2
1059 buop = 0.5 * (gotvx(i,1) * sflux(i) + tem)
1061 tem1 = dku(i,1) * shr2(i,1)
1063 tem = (u1(i,2)-u1(i,1))*rdzt(i,1)
1070 if(scuflg(i) .and. mrad(i) == 1)
then
1071 ptem = ucdo(i,1)+ucdo(i,2)-u1(i,1)-u1(i,2)
1072 ptem = 0.5 * tem * xmfd(i,1) * ptem
1076 ptem1 = ptem1 + ptem
1078 tem = (v1(i,2)-v1(i,1))*rdzt(i,1)
1085 if(scuflg(i) .and. mrad(i) == 1)
then
1086 ptem = vcdo(i,1)+vcdo(i,2)-v1(i,1)-v1(i,2)
1087 ptem = 0.5 * tem * xmfd(i,1) * ptem
1091 ptem2 = ptem2 + ptem
1094 tem2 = stress(i)*ustar(i)*phim(i)/(vk*zl(i,1))
1095 shrp = 0.5 * (tem1 + ptem1 + ptem2 + tem2)
1097 tem1 = -dkt(i,k-1) * bf(i,k-1)
1098 tem2 = -dkt(i,k) * bf(i,k)
1099 tem = 0.5 * (tem1 + tem2)
1100 if(pcnvflg(i) .and. k <= kpbl(i))
then
1101 ptem = 0.5 * (xmf(i,k-1) + xmf(i,k))
1102 ptem1 = ptem * buou(i,k)
1107 if(k >= mrad(i) .and. k < krad(i))
then
1108 ptem0 = 0.5 * (xmfd(i,k-1) + xmfd(i,k))
1109 ptem2 = ptem0 * buod(i,k)
1116 buop = tem + ptem1 + ptem2
1118 tem1 = dku(i,k-1) * shr2(i,k-1)
1119 tem2 = dku(i,k) * shr2(i,k)
1120 tem = 0.5 * (tem1 + tem2)
1121 tem1 = (u1(i,k+1)-u1(i,k))*rdzt(i,k)
1122 tem2 = (u1(i,k)-u1(i,k-1))*rdzt(i,k-1)
1123 if(pcnvflg(i) .and. k <= kpbl(i))
then
1124 ptem = xmf(i,k) * tem1 + xmf(i,k-1) * tem2
1125 ptem1 = 0.5 * ptem * (u1(i,k)-ucko(i,k))
1130 if(k >= mrad(i) .and. k < krad(i))
then
1131 ptem0 = xmfd(i,k) * tem1 + xmfd(i,k-1) * tem2
1132 ptem2 = 0.5 * ptem0 * (ucdo(i,k)-u1(i,k))
1139 shrp = tem + ptem1 + ptem2
1140 tem1 = (v1(i,k+1)-v1(i,k))*rdzt(i,k)
1141 tem2 = (v1(i,k)-v1(i,k-1))*rdzt(i,k-1)
1142 if(pcnvflg(i) .and. k <= kpbl(i))
then
1143 ptem = xmf(i,k) * tem1 + xmf(i,k-1) * tem2
1144 ptem1 = 0.5 * ptem * (v1(i,k)-vcko(i,k))
1149 if(k >= mrad(i) .and. k < krad(i))
then
1150 ptem0 = xmfd(i,k) * tem1 + xmfd(i,k-1) * tem2
1151 ptem2 = 0.5 * ptem0 * (vcdo(i,k)-v1(i,k))
1158 shrp = shrp + ptem1 + ptem2
1160 prod(i,k) = buop + shrp
1169 rle(i,k) = ce0 / ele(i,k)
1172 kk = max(nint(dt2/cdtn), 1)
1173 dtn = dt2 / float(kk)
1177 tem = sqrt(tke(i,k))
1178 diss(i,k) = rle(i,k) * tke(i,k) * tem
1179 tem1 = prod(i,k) + tke(i,k) / dtn
1180 diss(i,k)=max(min(diss(i,k), tem1), 0.)
1181 tke(i,k) = tke(i,k) + dtn * (prod(i,k)-diss(i,k))
1182 tke(i,k) = max(tke(i,k), tkmin)
1192 qcko(i,k,ntke) = tke(i,k)
1195 qcdo(i,k,ntke) = tke(i,k)
1201 if (pcnvflg(i) .and. k <= kpbl(i))
then
1202 dz = zl(i,k) - zl(i,k-1)
1203 tem = 0.5 * xlamue(i,k-1) * dz
1205 qcko(i,k,ntke)=((1.-tem)*qcko(i,k-1,ntke)+tem*
1206 & (tke(i,k)+tke(i,k-1)))/factor
1212 if (scuflg(i) .and. k < krad(i))
then
1213 if(k >= mrad(i))
then
1214 dz = zl(i,k+1) - zl(i,k)
1215 tem = 0.5 * xlamde(i,k) * dz
1217 qcdo(i,k,ntke)=((1.-tem)*qcdo(i,k+1,ntke)+tem*
1218 & (tke(i,k)+tke(i,k+1)))/factor
1234 dtodsd = dt2/del(i,k)
1235 dtodsu = dt2/del(i,k+1)
1236 dsig = prsl(i,k)-prsl(i,k+1)
1238 tem1 = dsig * dkq(i,k) * rdz
1240 au(i,k) = -dtodsd*dsdz2
1241 al(i,k) = -dtodsu*dsdz2
1242 ad(i,k) = ad(i,k)-au(i,k)
1243 ad(i,k+1)= 1.-al(i,k)
1246 if(pcnvflg(i) .and. k < kpbl(i))
then
1247 ptem = 0.5 * tem2 * xmf(i,k)
1248 ptem1 = dtodsd * ptem
1249 ptem2 = dtodsu * ptem
1250 tem = tke(i,k) + tke(i,k+1)
1251 ptem = qcko(i,k,ntke) + qcko(i,k+1,ntke)
1252 f1(i,k) = f1(i,k)-(ptem-tem)*ptem1
1253 f1(i,k+1) = tke(i,k+1)+(ptem-tem)*ptem2
1255 f1(i,k+1) = tke(i,k+1)
1259 if(k >= mrad(i) .and. k < krad(i))
then
1260 ptem = 0.5 * tem2 * xmfd(i,k)
1261 ptem1 = dtodsd * ptem
1262 ptem2 = dtodsu * ptem
1263 tem = tke(i,k) + tke(i,k+1)
1264 ptem = qcdo(i,k,ntke) + qcdo(i,k+1,ntke)
1265 f1(i,k) = f1(i,k) + (ptem - tem) * ptem1
1266 f1(i,k+1) = f1(i,k+1) - (ptem - tem) * ptem2
1275 call tridit(im,km,1,al,ad,au,f1,au,f1)
1282 qtend = (f1(i,k)-q1(i,k,ntke))*rdt
1283 rtg(i,k,ntke) = rtg(i,k,ntke)+qtend
1291 f1(i,1) = t1(i,1) + dtdz1(i) * heat(i)
1292 f2(i,1) = q1(i,1,1) + dtdz1(i) * evap(i)
1294 if(ntrac1 >= 2)
then
1298 f2(i,1+is) = q1(i,1,kk)
1305 dtodsd = dt2/del(i,k)
1306 dtodsu = dt2/del(i,k+1)
1307 dsig = prsl(i,k)-prsl(i,k+1)
1309 tem1 = dsig * dkt(i,k) * rdz
1312 au(i,k) = -dtodsd*dsdz2
1313 al(i,k) = -dtodsu*dsdz2
1314 ad(i,k) = ad(i,k)-au(i,k)
1315 ad(i,k+1)= 1.-al(i,k)
1318 if(pcnvflg(i) .and. k < kpbl(i))
then
1319 ptem = 0.5 * tem2 * xmf(i,k)
1320 ptem1 = dtodsd * ptem
1321 ptem2 = dtodsu * ptem
1322 tem = t1(i,k) + t1(i,k+1)
1323 ptem = tcko(i,k) + tcko(i,k+1)
1324 f1(i,k) = f1(i,k)+dtodsd*dsdzt-(ptem-tem)*ptem1
1325 f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt+(ptem-tem)*ptem2
1326 tem = q1(i,k,1) + q1(i,k+1,1)
1327 ptem = qcko(i,k,1) + qcko(i,k+1,1)
1328 f2(i,k) = f2(i,k) - (ptem - tem) * ptem1
1329 f2(i,k+1) = q1(i,k+1,1) + (ptem - tem) * ptem2
1331 f1(i,k) = f1(i,k)+dtodsd*dsdzt
1332 f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt
1333 f2(i,k+1) = q1(i,k+1,1)
1337 if(k >= mrad(i) .and. k < krad(i))
then
1338 ptem = 0.5 * tem2 * xmfd(i,k)
1339 ptem1 = dtodsd * ptem
1340 ptem2 = dtodsu * ptem
1341 ptem = tcdo(i,k) + tcdo(i,k+1)
1342 tem = t1(i,k) + t1(i,k+1)
1343 f1(i,k) = f1(i,k) + (ptem - tem) * ptem1
1344 f1(i,k+1) = f1(i,k+1) - (ptem - tem) * ptem2
1345 tem = q1(i,k,1) + q1(i,k+1,1)
1346 ptem = qcdo(i,k,1) + qcdo(i,k+1,1)
1347 f2(i,k) = f2(i,k) + (ptem - tem) * ptem1
1348 f2(i,k+1) = f2(i,k+1) - (ptem - tem) * ptem2
1354 if(ntrac1 >= 2)
then
1359 if(pcnvflg(i) .and. k < kpbl(i))
then
1360 dtodsd = dt2/del(i,k)
1361 dtodsu = dt2/del(i,k+1)
1362 dsig = prsl(i,k)-prsl(i,k+1)
1363 tem = dsig * rdzt(i,k)
1364 ptem = 0.5 * tem * xmf(i,k)
1365 ptem1 = dtodsd * ptem
1366 ptem2 = dtodsu * ptem
1367 tem1 = qcko(i,k,kk) + qcko(i,k+1,kk)
1368 tem2 = q1(i,k,kk) + q1(i,k+1,kk)
1369 f2(i,k+is) = f2(i,k+is) - (tem1 - tem2) * ptem1
1370 f2(i,k+1+is)= q1(i,k+1,kk) + (tem1 - tem2) * ptem2
1372 f2(i,k+1+is) = q1(i,k+1,kk)
1376 if(k >= mrad(i) .and. k < krad(i))
then
1377 dtodsd = dt2/del(i,k)
1378 dtodsu = dt2/del(i,k+1)
1379 dsig = prsl(i,k)-prsl(i,k+1)
1380 tem = dsig * rdzt(i,k)
1381 ptem = 0.5 * tem * xmfd(i,k)
1382 ptem1 = dtodsd * ptem
1383 ptem2 = dtodsu * ptem
1384 tem1 = qcdo(i,k,kk) + qcdo(i,k+1,kk)
1385 tem2 = q1(i,k,kk) + q1(i,k+1,kk)
1386 f2(i,k+is) = f2(i,k+is) + (tem1 - tem2) * ptem1
1387 f2(i,k+1+is)= f2(i,k+1+is) - (tem1 - tem2) * ptem2
1398 call tridin(im,km,ntrac1,al,ad,au,f1,f2,au,f1,f2)
1404 ttend = (f1(i,k)-t1(i,k))*rdt
1405 qtend = (f2(i,k)-q1(i,k,1))*rdt
1406 tdt(i,k) = tdt(i,k)+ttend
1407 rtg(i,k,1) = rtg(i,k,1)+qtend
1408 dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend
1409 dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend
1412 if (ldiag3d .and. .not. gen_tend)
then
1413 idtend = dtidx(index_of_temperature,index_of_process_pbl)
1415 dtend(:,:,idtend) = dtend(:,:,idtend) + delt*rdt*(f1-t1)
1417 idtend = dtidx(100+ntqv,index_of_process_pbl)
1419 dtend(:,:,idtend) = dtend(:,:,idtend) +delt*rdt*(f2-q1(:,:,1))
1423 if(ntrac1 >= 2)
then
1428 qtend = (f2(i,k+is)-q1(i,k,kk))*rdt
1429 rtg(i,k,kk) = rtg(i,k,kk)+qtend
1442 ttend = diss(i,k) / cp
1443 tdt(i,k) = tdt(i,k) + dspfac * ttend
1451 ad(i,1) = 1.0 + dtdz1(i) * stress(i) / spd1(i)
1458 dtodsd = dt2/del(i,k)
1459 dtodsu = dt2/del(i,k+1)
1460 dsig = prsl(i,k)-prsl(i,k+1)
1462 tem1 = dsig * dku(i,k) * rdz
1464 au(i,k) = -dtodsd*dsdz2
1465 al(i,k) = -dtodsu*dsdz2
1466 ad(i,k) = ad(i,k)-au(i,k)
1467 ad(i,k+1)= 1.-al(i,k)
1470 if(pcnvflg(i) .and. k < kpbl(i))
then
1471 ptem = 0.5 * tem2 * xmf(i,k)
1472 ptem1 = dtodsd * ptem
1473 ptem2 = dtodsu * ptem
1474 tem = u1(i,k) + u1(i,k+1)
1475 ptem = ucko(i,k) + ucko(i,k+1)
1476 f1(i,k) = f1(i,k) - (ptem - tem) * ptem1
1477 f1(i,k+1) = u1(i,k+1) + (ptem - tem) * ptem2
1478 tem = v1(i,k) + v1(i,k+1)
1479 ptem = vcko(i,k) + vcko(i,k+1)
1480 f2(i,k) = f2(i,k) - (ptem - tem) * ptem1
1481 f2(i,k+1) = v1(i,k+1) + (ptem - tem) * ptem2
1483 f1(i,k+1) = u1(i,k+1)
1484 f2(i,k+1) = v1(i,k+1)
1488 if(k >= mrad(i) .and. k < krad(i))
then
1489 ptem = 0.5 * tem2 * xmfd(i,k)
1490 ptem1 = dtodsd * ptem
1491 ptem2 = dtodsu * ptem
1492 tem = u1(i,k) + u1(i,k+1)
1493 ptem = ucdo(i,k) + ucdo(i,k+1)
1494 f1(i,k) = f1(i,k) + (ptem - tem) *ptem1
1495 f1(i,k+1) = f1(i,k+1) - (ptem - tem) *ptem2
1496 tem = v1(i,k) + v1(i,k+1)
1497 ptem = vcdo(i,k) + vcdo(i,k+1)
1498 f2(i,k) = f2(i,k) + (ptem - tem) * ptem1
1499 f2(i,k+1) = f2(i,k+1) - (ptem - tem) * ptem2
1508 call tridi2(im,km,al,ad,au,f1,f2,au,f1,f2)
1514 utend = (f1(i,k)-u1(i,k))*rdt
1515 vtend = (f2(i,k)-v1(i,k))*rdt
1516 du(i,k) = du(i,k)+utend
1517 dv(i,k) = dv(i,k)+vtend
1518 dusfc(i) = dusfc(i)+conw*del(i,k)*utend
1519 dvsfc(i) = dvsfc(i)+conw*del(i,k)*vtend
1522 if (ldiag3d .and. .not. gen_tend)
then
1523 idtend=dtidx(index_of_x_wind,index_of_process_pbl)
1525 dtend(:,:,idtend) = dtend(:,:,idtend) + delt*rdt*(f1-u1)
1527 idtend=dtidx(index_of_y_wind,index_of_process_pbl)
1529 dtend(:,:,idtend) = dtend(:,:,idtend) + delt*rdt*(f2-v1)