40 subroutine ysuvdif_run(im,km,ux,vx,tx,qx,p2d,p2di,pi2d,karman, &
41 utnp,vtnp,ttnp,qtnp, &
42 swh,hlw,xmu,ntrac,ndiff,ntcw,ntiw, &
44 zorl,stress,hpbl,psim,psih, &
45 landmask,heat,evap,wspd,br, &
46 g,rd,cp,rv,ep1,ep2,xlv, &
47 dusfc,dvsfc,dtsfc,dqsfc, &
48 dt,kpbl1d,u10,v10,lssav,ldiag3d,qdiag3d, &
49 flag_for_pbl_generic_tend,ntoz,ntqv,dtend,dtidx, &
50 index_of_temperature,index_of_x_wind,index_of_y_wind, &
51 index_of_process_pbl,errmsg,errflg )
58 real(kind=kind_phys),
parameter :: xkzminm = 0.1,xkzminh = 0.01
59 real(kind=kind_phys),
parameter :: xkzmin = 0.01,xkzmax = 1000.,rimin = -100.
60 real(kind=kind_phys),
parameter :: rlam = 30.,prmin = 0.25,prmax = 4.
61 real(kind=kind_phys),
parameter :: brcr_ub = 0.0,brcr_sb = 0.25,cori = 1.e-4
62 real(kind=kind_phys),
parameter :: afac = 6.8,bfac = 6.8,pfac = 2.0,pfac_q = 2.0
63 real(kind=kind_phys),
parameter :: phifac = 8.,sfcfrac = 0.1
64 real(kind=kind_phys),
parameter :: d1 = 0.02, d2 = 0.05, d3 = 0.001
65 real(kind=kind_phys),
parameter :: h1 = 0.33333335, h2 = 0.6666667
66 real(kind=kind_phys),
parameter :: zfmin = 1.e-8,aphi5 = 5.,aphi16 = 16.
67 real(kind=kind_phys),
parameter :: tmin=1.e-2
68 real(kind=kind_phys),
parameter :: gamcrt = 3.,gamcrq = 2.e-3
69 real(kind=kind_phys),
parameter :: xka = 2.4e-5
70 real(kind=kind_phys),
parameter :: rcl = 1.0
71 real(kind=kind_phys),
intent(in) :: karman
72 integer,
parameter :: imvdif = 1
73 integer,
parameter :: ysu_topdown_pblmix = 1
77 integer,
intent(in ) :: im,km,ntrac,ndiff,ntcw,ntiw,ntoz
78 real(kind=kind_phys),
intent(in ) :: g,cp,rd,rv,ep1,ep2,xlv,dt
80 real(kind=kind_phys),
dimension( :,: ), &
81 intent(in) :: pi2d,p2d,phil,ux,vx,swh,hlw,tx
83 real(kind=kind_phys),
dimension( :,:,: ) , &
86 real(kind=kind_phys),
dimension( :,: ) , &
87 intent(in ) :: p2di,phii
89 real(kind=kind_phys),
dimension( : ) , &
90 intent(in) :: stress,zorl,heat,evap,wspd,br,psim,psih,psfcpa, &
92 integer,
dimension(:) ,&
93 intent(in ) :: landmask
94 logical,
intent(in ) :: lssav, ldiag3d, qdiag3d, &
95 flag_for_pbl_generic_tend
100 real(kind=kind_phys),
dimension( :,: ) , &
101 intent(inout) :: utnp,vtnp,ttnp
102 real(kind=kind_phys),
dimension( :,:,: ) , &
103 intent(inout) :: qtnp
104 real(kind=kind_phys),
optional,
intent(inout),
optional :: dtend(:,:,:)
105 integer,
intent(in) :: dtidx(:,:), ntqv, index_of_temperature, &
106 index_of_x_wind, index_of_y_wind, index_of_process_pbl
110 integer,
dimension( : ),
intent(out ) :: kpbl1d
111 real(kind=kind_phys),
dimension( : ), &
113 real(kind=kind_phys),
dimension( : ), &
114 intent(out) :: dusfc,dvsfc, dtsfc,dqsfc
117 character(len=*),
intent(out) :: errmsg
118 integer,
intent(out) :: errflg
124 real(kind=kind_phys),
dimension( im ) :: hol
125 real(kind=kind_phys),
dimension( im, km+1 ) :: zq
127 real(kind=kind_phys),
dimension( im, km ) :: &
136 real(kind=kind_phys),
dimension( im ) :: &
147 real(kind=kind_phys),
dimension( im, km ) :: xkzm,xkzh, &
158 real(kind=kind_phys),
dimension( im ) :: &
164 real(kind=kind_phys),
dimension( im ) :: xland
165 real(kind=kind_phys),
dimension( im ) :: ust
166 real(kind=kind_phys),
dimension( im ) :: hfx
167 real(kind=kind_phys),
dimension( im ) :: qfx
168 real(kind=kind_phys),
dimension( im ) :: znt
169 real(kind=kind_phys),
dimension( im ) :: uox
170 real(kind=kind_phys),
dimension( im ) :: vox
172 real(kind=kind_phys),
dimension( im, km, ndiff) :: r3,f3
173 integer,
dimension( im ) :: kpbl,kpblold
175 logical,
dimension( im ) :: pblflg, &
180 logical :: definebrup
182 integer :: n,i,k,l,ic,is,kk
183 integer :: klpbl, ktrace1, ktrace2, ktrace3
186 real(kind=kind_phys) :: dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum,prnum0
187 real(kind=kind_phys) :: ss,ri,qmean,tmean,alph,chi,zk,rl2,dk,sri
188 real(kind=kind_phys) :: brint,dtodsd,dtodsu,rdz,dsdzt,dsdzq,dsdz2,rlamdz
189 real(kind=kind_phys) :: utend,vtend,ttend,qtend
190 real(kind=kind_phys) :: dtstep,govrthv
191 real(kind=kind_phys) :: cont, conq, conw, conwrc, rovcp
194 real(kind=kind_phys),
dimension( im, km ) :: wscalek,wscalek2
195 real(kind=kind_phys),
dimension( im ) :: wstar
196 real(kind=kind_phys),
dimension( im ) :: delta
197 real(kind=kind_phys),
dimension( im, km ) :: xkzml,xkzhl, &
199 real(kind=kind_phys),
dimension( im ) :: ust3, &
208 real(kind=kind_phys) :: prnumfac,bfx0,hfx0,qfx0,delb,dux,dvx, &
209 dsdzu,dsdzv,wm3,dthx,dqx,wspd10,ross,tem1,dsig,tvcon,conpr, &
210 prfac,prfac2,phim8z,radsum,tmp1,templ,rvls,temps,ent_eff, &
211 rcldb,bruptmp,radflux
226 conwrc = conw*sqrt(rcl)
227 conpr = bfac*karman*sfcfrac
231 if(landmask(i).eq.0)
then
240 thx(i,k) = tx(i,k)/pi2d(i,k)
241 thlix(i,k) = (tx(i,k)-xlv*qx(i,k,ntcw)/cp-2.834e6*qx(i,k,ntiw)/cp)/pi2d(i,k)
247 tvcon = (1.+ep1*qx(i,k,1))
248 thvx(i,k) = thx(i,k)*tvcon
253 tvcon = (1.+ep1*qx(i,1,1))
254 rhox(i) = psfcpa(i)/(rd*tx(i,1)*tvcon)
255 govrth(i) = g/thx(i,1)
256 hfx(i) = heat(i)*rhox(i)*cp
257 qfx(i) = evap(i)*rhox(i)
258 ust(i) = sqrt(stress(i))
259 znt(i) = 0.01*zorl(i)
273 zq(i,k+1) = phii(i,k+1)*conw
274 tvcon = (1.+ep1*qx(i,k,1))
275 rhox2(i,k) = p2d(i,k)/(rd*tx(i,k)*tvcon)
281 za(i,k) = phil(i,k)*conw
282 dzq(i,k) = zq(i,k+1)-zq(i,k)
283 del(i,k) = p2di(i,k)-p2di(i,k+1)
293 dza(i,k) = za(i,k)-za(i,k-1)
326 wspd1(i) = sqrt( (ux(i,1)-uox(i))*(ux(i,1)-uox(i)) + (vx(i,1)-vox(i))*(vx(i,1)-vox(i)) )+1.e-9
383 thermal(i)= thvx(i,1)
384 thermalli(i) = thlix(i,1)
387 sflux(i) = hfx(i)/rhox(i)/cp + qfx(i)/rhox(i)*ep1*thx(i,1)
388 if(br(i).gt.0.0) sfcflg(i) = .false.
401 if(.not.stable(i))
then
403 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
404 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
406 stable(i) = brup(i).gt.brcr(i)
413 if(brdn(i).ge.brcr(i))
then
415 elseif(brup(i).le.brcr(i))
then
418 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
420 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
421 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
422 if(kpbl(i).le.1) pblflg(i) = .false.
428 zol1(i) = max(br(i)*fm*fm/fh,rimin)
430 zol1(i) = min(zol1(i),-zfmin)
432 zol1(i) = max(zol1(i),zfmin)
434 hol1 = zol1(i)*hpbl(i)/zl1(i)*sfcfrac
436 phim(i) = (1.-aphi16*hol1)**(-1./4.)
437 phih(i) = (1.-aphi16*hol1)**(-1./2.)
438 bfx0 = max(sflux(i),0.)
439 hfx0 = max(hfx(i)/rhox(i)/cp,0.)
440 qfx0 = max(ep1*thx(i,1)*qfx(i)/rhox(i),0.)
441 wstar3(i) = (govrth(i)*bfx0*hpbl(i))
442 wstar(i) = (wstar3(i))**h1
444 phim(i) = (1.+aphi5*hol1)
450 wscale(i) = (ust3(i)+phifac*karman*wstar3(i)*0.5)**h1
451 wscale(i) = min(wscale(i),ust(i)*aphi16)
452 wscale(i) = max(wscale(i),ust(i)/aphi5)
459 if(sfcflg(i).and.sflux(i).gt.0.0)
then
460 gamfac = bfac/rhox(i)/wscale(i)
461 hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt)
462 hgamq(i) = min(gamfac*qfx(i),gamcrq)
463 vpert = (hgamt(i)+ep1*thx(i,1)*hgamq(i))/bfac*afac
464 thermal(i) = thermal(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0)
465 thermalli(i)= thermalli(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0)
466 hgamt(i) = max(hgamt(i),0.0)
467 hgamq(i) = max(hgamq(i),0.0)
468 brint = -15.9*ust(i)*ust(i)/wspd(i)*wstar3(i)/(wscale(i)**4.)
469 hgamu(i) = brint*ux(i,1)
470 hgamv(i) = brint*vx(i,1)
495 if(.not.stable(i).and.pblflg(i))
then
497 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
498 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
500 stable(i) = brup(i).gt.brcr(i)
507 if (ysu_topdown_pblmix.eq.1)
then
511 do k = kpblold(i), km-1
512 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
513 bruptmp = (thlix(i,k)-thermalli(i))*(g*za(i,k)/thlix(i,1))/spdk2
514 stable(i) = bruptmp.ge.brcr(i)
520 if (.not.stable(i))
then
532 if(brdn(i).ge.brcr(i))
then
534 elseif(brup(i).le.brcr(i))
then
537 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
539 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
540 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
541 if(kpbl(i).le.1) pblflg(i) = .false.
548 if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2))
then
557 if((.not.stable(i)).and.((xland(i)-1.5).ge.0))
then
558 wspd10 = u10(i)*u10(i) + v10(i)*v10(i)
559 wspd10 = sqrt(wspd10)
560 ross = wspd10 / (cori*znt(i))
561 brcr_sbro(i) = min(0.16*(1.e-7*ross)**(-0.18),.3)
566 if(.not.stable(i))
then
567 if((xland(i)-1.5).ge.0)
then
568 brcr(i) = brcr_sbro(i)
577 if(.not.stable(i))
then
579 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
580 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
582 stable(i) = brup(i).gt.brcr(i)
588 if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2))
then
590 if(brdn(i).ge.brcr(i))
then
592 elseif(brup(i).le.brcr(i))
then
595 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
597 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
598 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
599 if(kpbl(i).le.1) pblflg(i) = .false.
609 wm3 = wstar3(i) + 5. * ust3(i)
611 bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i)
612 dthvx(i) = max(thvx(i,k+1)-thvx(i,k),tmin)
613 we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i)))
614 if((qx(i,k,ntcw)+qx(i,k,ntiw)).gt.0.01e-3.and.ysu_topdown_pblmix.eq.1)
then
615 if ( kpbl(i) .ge. 2)
then
617 templ=thlix(i,k)*(p2di(i,k+1)/100000)**rovcp
619 rvls=100.*6.112*exp(17.67*(templ-273.16)/(templ-29.65))*(ep2/p2di(i,k+1))
620 temps=templ + ((qx(i,k,1)+qx(i,k,ntcw))-rvls)/(cp/xlv + &
621 ep2*xlv*rvls/(rd*templ**2))
622 rvls=100.*6.112*exp(17.67*(temps-273.15)/(temps-29.65))*(ep2/p2di(i,k+1))
623 rcldb=max((qx(i,k,1)+qx(i,k,ntcw))-rvls,0.)
625 dthvx(i) = (thlix(i,k+2)+thx(i,k+2)*ep1*(qx(i,k+2,1)+qx(i,k+2,ntcw))) &
626 - (thlix(i,k) + thx(i,k) *ep1*(qx(i,k,1) +qx(i,k,ntcw)))
627 dthvx(i) = max(dthvx(i),0.1)
628 tmp1 = xlv/cp * rcldb/(pi2d(i,k)*dthvx(i))
629 ent_eff = 0.2 * 8. * tmp1 +0.2
633 radflux=swh(i,kk)*xmu(i)+hlw(i,kk)
634 radflux=radflux*cp/g*(p2di(i,kk)-p2di(i,kk+1))
635 if (radflux < 0.0 ) radsum=abs(radflux)+radsum
637 radsum=max(radsum,0.0)
640 bfx0 = max(max(sflux(i),0.0)-radsum/rhox2(i,k)/cp,0.)
641 bfx0 = max(sflux(i),0.0)
642 wm3 = (govrth(i)*bfx0*hpbl(i))+5. * ust3(i)
644 bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i)
645 dthvx(i) = max(thvx(i,k+1)-thvx(i,k),tmin)
646 we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i)))
649 bfx0 = max(radsum/rhox2(i,k)/cp-max(sflux(i),0.0),0.)
650 bfx0 = max(radsum/rhox2(i,k)/cp,0.)
651 wm3 = (g/thvx(i,k)*bfx0*hpbl(i))
652 wm2(i) = wm2(i)+wm3**h2
653 bfxpbl(i) = - ent_eff * bfx0
654 dthvx(i) = max(thvx(i,k+1)-thvx(i,k),0.1)
655 we(i) = we(i) + max(bfxpbl(i)/dthvx(i),-sqrt(wm3**h2))
658 bfx0 = max(radsum/rhox2(i,k)/cp,0.)
659 wstar3_2(i) = (g/thvx(i,k)*bfx0*hpbl(i))
661 wscale(i) = (ust3(i)+phifac*karman*(wstar3(i)+wstar3_2(i))*0.5)**h1
662 wscale(i) = min(wscale(i),ust(i)*aphi16)
663 wscale(i) = max(wscale(i),ust(i)/aphi5)
664 gamfac = bfac/rhox(i)/wscale(i)
665 hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt)
666 hgamq(i) = min(gamfac*qfx(i),gamcrq)
667 gamfac = bfac/rhox2(i,k)/wscale(i)
668 hgamt2(i,k) = min(gamfac*radsum/cp,gamcrt)
669 hgamt(i) = max(hgamt(i),0.0) + max(hgamt2(i,k),0.0)
670 brint = -15.9*ust(i)*ust(i)/wspd(i)*(wstar3(i)+wstar3_2(i))/(wscale(i)**4.)
671 hgamu(i) = brint*ux(i,1)
672 hgamv(i) = brint*vx(i,1)
676 dthx = max(thx(i,k+1)-thx(i,k),tmin)
677 dqx = min(qx(i,k+1,1)-qx(i,k,1),0.0)
678 hfxpbl(i) = we(i)*dthx
679 qfxpbl(i) = we(i)*dqx
681 dux = ux(i,k+1)-ux(i,k)
682 dvx = vx(i,k+1)-vx(i,k)
684 ufxpbl(i) = max(prpbl(i)*we(i)*dux,-ust(i)*ust(i))
685 elseif(dux.lt.-tmin)
then
686 ufxpbl(i) = min(prpbl(i)*we(i)*dux,ust(i)*ust(i))
691 vfxpbl(i) = max(prpbl(i)*we(i)*dvx,-ust(i)*ust(i))
692 elseif(dvx.lt.-tmin)
then
693 vfxpbl(i) = min(prpbl(i)*we(i)*dvx,ust(i)*ust(i))
697 delb = govrth(i)*d3*hpbl(i)
698 delta(i) = min(d1*hpbl(i) + d2*wm2(i)/delb,100.)
704 if(pblflg(i).and.k.ge.kpbl(i))
then
705 entfac(i,k) = ((zq(i,k+1)-hpbl(i))/delta(i))**2.
716 if(k.lt.kpbl(i))
then
717 zfac(i,k) = min(max((1.-(zq(i,k+1)-zl1(i))/(hpbl(i)-zl1(i))),zfmin),1.)
718 zfacent(i,k) = (1.-zfac(i,k))**3.
719 wscalek(i,k) = (ust3(i)+phifac*karman*wstar3(i)*(1.-zfac(i,k)))**h1
720 wscalek2(i,k) = (phifac*karman*wstar3_2(i)*(zfac(i,k)))**h1
723 prfac2 = 15.9*(wstar3(i)+wstar3_2(i))/ust3(i)/(1.+4.*karman*(wstar3(i)+wstar3_2(i))/ust3(i))
724 prnumfac = -3.*(max(zq(i,k+1)-sfcfrac*hpbl(i),0.))**2./hpbl(i)**2.
729 phim8z = 1.+aphi5*zol1(i)*zq(i,k+1)/zl1(i)
730 wscalek(i,k) = ust(i)/phim8z
731 wscalek(i,k) = max(wscalek(i,k),0.001)
733 prnum0 = (phih(i)/phim(i)+prfac)
734 prnum0 = max(min(prnum0,prmax),prmin)
735 xkzm(i,k) = wscalek(i,k) *karman* zq(i,k+1) * zfac(i,k)**pfac+ &
736 wscalek2(i,k)*karman*(hpbl(i)-zq(i,k+1))*(1-zfac(i,k))**pfac
738 if (k.eq.kpbl(i)-1.and.cloudflg(i).and.we(i).lt.0.0)
then
741 prnum = 1. + (prnum0-1.)*exp(prnumfac)
742 xkzq(i,k) = xkzm(i,k)/prnum*zfac(i,k)**(pfac_q-pfac)
743 prnum0 = prnum0/(1.+prfac2*karman*sfcfrac)
744 prnum = 1. + (prnum0-1.)*exp(prnumfac)
745 xkzh(i,k) = xkzm(i,k)/prnum
746 xkzm(i,k) = xkzm(i,k)+xkzom(i,k)
747 xkzh(i,k) = xkzh(i,k)+xkzoh(i,k)
748 xkzq(i,k) = xkzq(i,k)+xkzoh(i,k)
749 xkzm(i,k) = min(xkzm(i,k),xkzmax)
750 xkzh(i,k) = min(xkzh(i,k),xkzmax)
751 xkzq(i,k) = min(xkzq(i,k),xkzmax)
760 if(k.ge.kpbl(i))
then
761 ss = ((ux(i,k+1)-ux(i,k))*(ux(i,k+1)-ux(i,k)) &
762 +(vx(i,k+1)-vx(i,k))*(vx(i,k+1)-vx(i,k))) &
763 /(dza(i,k+1)*dza(i,k+1))+1.e-9
764 govrthv = g/(0.5*(thvx(i,k+1)+thvx(i,k)))
765 ri = govrthv*(thvx(i,k+1)-thvx(i,k))/(ss*dza(i,k+1))
766 if(imvdif.eq.1.and.ntcw.ge.2.and.ntiw.ge.2)
then
767 if((qx(i,k,ntcw)+qx(i,k,ntiw)).gt.0.01e-3.and.(qx(i &
768 ,k+1,ntcw)+qx(i,k+1,ntiw)).gt.0.01e-3)
then
770 qmean = 0.5*(qx(i,k,1)+qx(i,k+1,1))
771 tmean = 0.5*(tx(i,k)+tx(i,k+1))
772 alph = xlv*qmean/rd/tmean
773 chi = xlv*xlv*qmean/cp/rv/tmean/tmean
774 ri = (1.+alph)*(ri-g*g/ss/tmean/cp*((chi-alph)/(1.+chi)))
777 zk = karman*zq(i,k+1)
778 rlamdz = min(max(0.1*dza(i,k+1),rlam),300.)
779 rlamdz = min(dza(i,k+1),rlamdz)
780 rl2 = (zk*rlamdz/(rlamdz+zk))**2
786 xkzm(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri))
787 xkzh(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri))
790 xkzh(i,k) = dk/(1+5.*ri)**2
792 prnum = min(prnum,prmax)
793 xkzm(i,k) = xkzh(i,k)*prnum
796 xkzm(i,k) = xkzm(i,k)+xkzom(i,k)
797 xkzh(i,k) = xkzh(i,k)+xkzoh(i,k)
798 xkzm(i,k) = min(xkzm(i,k),xkzmax)
799 xkzh(i,k) = min(xkzh(i,k),xkzmax)
800 xkzml(i,k) = xkzm(i,k)
801 xkzhl(i,k) = xkzh(i,k)
819 f1(i,1) = thx(i,1)-300.+hfx(i)/cont/del(i,1)*dt2
824 dtodsd = dt2/del(i,k)
825 dtodsu = dt2/del(i,k+1)
826 dsig = p2d(i,k)-p2d(i,k+1)
828 tem1 = dsig*xkzh(i,k)*rdz
829 if(pblflg(i).and.k.lt.kpbl(i))
then
830 dsdzt = tem1*(-hgamt(i)/hpbl(i)-hfxpbl(i)*zfacent(i,k)/xkzh(i,k))
831 f1(i,k) = f1(i,k)+dtodsd*dsdzt
832 f1(i,k+1) = thx(i,k+1)-300.-dtodsu*dsdzt
833 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6)
then
834 xkzh(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
835 xkzh(i,k) = sqrt(xkzh(i,k)*xkzhl(i,k))
836 xkzh(i,k) = max(xkzh(i,k),xkzoh(i,k))
837 xkzh(i,k) = min(xkzh(i,k),xkzmax)
838 f1(i,k+1) = thx(i,k+1)-300.
840 f1(i,k+1) = thx(i,k+1)-300.
842 tem1 = dsig*xkzh(i,k)*rdz
844 au(i,k) = -dtodsd*dsdz2
845 al(i,k) = -dtodsu*dsdz2
846 ad(i,k) = ad(i,k)-au(i,k)
847 ad(i,k+1) = 1.-al(i,k)
860 call tridin_ysu(al,ad,cu,r1,au,f1,im,km,1)
866 ttend = (f1(i,k)-thx(i,k)+300.)*rdt*pi2d(i,k)
867 ttnp(i,k) = ttnp(i,k)+ttend
868 dtsfc(i) = dtsfc(i)+ttend*cont*del(i,k)
871 if(lssav .and. ldiag3d .and. .not. flag_for_pbl_generic_tend)
then
872 idtend = dtidx(index_of_temperature,index_of_process_pbl)
874 dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f1-thx+300.)*rdt*pi2d
898 f3(i,1,1) = qx(i,1,1)+qfx(i)*g/del(i,1)*dt2
904 f3(i,1,ic) = qx(i,1,ic)
911 if(k.ge.kpbl(i))
then
912 xkzq(i,k) = xkzh(i,k)
919 dtodsd = dt2/del(i,k)
920 dtodsu = dt2/del(i,k+1)
921 dsig = p2d(i,k)-p2d(i,k+1)
923 tem1 = dsig*xkzq(i,k)*rdz
924 if(pblflg(i).and.k.lt.kpbl(i))
then
925 dsdzq = tem1*(-qfxpbl(i)*zfacent(i,k)/xkzq(i,k))
926 f3(i,k,1) = f3(i,k,1)+dtodsd*dsdzq
927 f3(i,k+1,1) = qx(i,k+1,1)-dtodsu*dsdzq
928 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6)
then
929 xkzq(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
930 xkzq(i,k) = sqrt(xkzq(i,k)*xkzhl(i,k))
931 xkzq(i,k) = max(xkzq(i,k),xkzoh(i,k))
932 xkzq(i,k) = min(xkzq(i,k),xkzmax)
933 f3(i,k+1,1) = qx(i,k+1,1)
935 f3(i,k+1,1) = qx(i,k+1,1)
937 tem1 = dsig*xkzq(i,k)*rdz
939 au(i,k) = -dtodsd*dsdz2
940 al(i,k) = -dtodsu*dsdz2
941 ad(i,k) = ad(i,k)-au(i,k)
942 ad(i,k+1) = 1.-al(i,k)
950 f3(i,k+1,ic) = qx(i,k+1,ic)
967 r3(i,k,ic) = f3(i,k,ic)
974 call tridin_ysu(al,ad,cu,r3,au,f3,im,km,ndiff)
980 qtend = (f3(i,k,1)-qx(i,k,1))*rdt
981 qtnp(i,k,1) = qtnp(i,k,1)+qtend
982 dqsfc(i) = dqsfc(i)+qtend*conq*del(i,k)
985 if(lssav .and. ldiag3d .and. qdiag3d .and. .not. flag_for_pbl_generic_tend)
then
986 idtend = dtidx(ntqv+100,index_of_process_pbl)
988 dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f3(:,:,1)-qx(:,:,1))*rdt
996 qtend = (f3(i,k,ic)-qx(i,k,ic))*rdt
997 qtnp(i,k,ic) = qtnp(i,k,ic)+qtend
1001 if(lssav .and. ldiag3d .and. ntoz>0 .and. qdiag3d .and. &
1002 & .not. flag_for_pbl_generic_tend)
then
1003 idtend = dtidx(100+ntoz,index_of_process_pbl)
1005 dtend(:,:,idtend) = dtend(:,:,idtend) + f3(:,:,ntoz)-qx(:,:,ntoz)
1023 ad(i,1) = 1.+ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2 &
1024 *(wspd1(i)/wspd(i))**2
1025 f1(i,1) = ux(i,1)+uox(i)*ust(i)**2*g/del(i,1)*dt2/wspd1(i)
1026 f2(i,1) = vx(i,1)+vox(i)*ust(i)**2*g/del(i,1)*dt2/wspd1(i)
1031 dtodsd = dt2/del(i,k)
1032 dtodsu = dt2/del(i,k+1)
1033 dsig = p2d(i,k)-p2d(i,k+1)
1035 tem1 = dsig*xkzm(i,k)*rdz
1036 if(pblflg(i).and.k.lt.kpbl(i))
then
1037 dsdzu = tem1*(-hgamu(i)/hpbl(i)-ufxpbl(i)*zfacent(i,k)/xkzm(i,k))
1038 dsdzv = tem1*(-hgamv(i)/hpbl(i)-vfxpbl(i)*zfacent(i,k)/xkzm(i,k))
1039 f1(i,k) = f1(i,k)+dtodsd*dsdzu
1040 f1(i,k+1) = ux(i,k+1)-dtodsu*dsdzu
1041 f2(i,k) = f2(i,k)+dtodsd*dsdzv
1042 f2(i,k+1) = vx(i,k+1)-dtodsu*dsdzv
1043 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6)
then
1044 xkzm(i,k) = prpbl(i)*xkzh(i,k)
1045 xkzm(i,k) = sqrt(xkzm(i,k)*xkzml(i,k))
1046 xkzm(i,k) = max(xkzm(i,k),xkzom(i,k))
1047 xkzm(i,k) = min(xkzm(i,k),xkzmax)
1048 f1(i,k+1) = ux(i,k+1)
1049 f2(i,k+1) = vx(i,k+1)
1051 f1(i,k+1) = ux(i,k+1)
1052 f2(i,k+1) = vx(i,k+1)
1054 tem1 = dsig*xkzm(i,k)*rdz
1056 au(i,k) = -dtodsd*dsdz2
1057 al(i,k) = -dtodsu*dsdz2
1058 ad(i,k) = ad(i,k)-au(i,k)
1059 ad(i,k+1) = 1.-al(i,k)
1075 call tridi1n(al,ad,cu,r1,r2,au,f1,f2,im,km,1)
1081 utend = (f1(i,k)-ux(i,k))*rdt
1082 vtend = (f2(i,k)-vx(i,k))*rdt
1083 utnp(i,k) = utnp(i,k)+utend
1084 vtnp(i,k) = vtnp(i,k)+vtend
1085 dusfc(i) = dusfc(i) + utend*conwrc*del(i,k)
1086 dvsfc(i) = dvsfc(i) + vtend*conwrc*del(i,k)
1089 if(lssav .and. ldiag3d .and. .not. flag_for_pbl_generic_tend)
then
1090 idtend = dtidx(index_of_x_wind,index_of_process_pbl)
1092 dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f1-ux)*rdt
1095 idtend = dtidx(index_of_y_wind,index_of_process_pbl)
1097 dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f2-vx)*rdt