40 subroutine shinhongvdif_run(im,km,ux,vx,tx,qx,p2d,p2di,pi2d,karman, &
41 utnp,vtnp,ttnp,qtnp,ntrac,ndiff,ntcw,ntiw, &
43 zorl,stress,hpbl,psim,psih, &
44 landmask,heat,evap,wspd,br, &
45 g,rd,cp,rv,ep1,ep2,xlv, &
46 dusfc,dvsfc,dtsfc,dqsfc, &
50 flag_for_pbl_generic_tend,ntoz,ntqv,dtend,dtidx, &
51 index_of_process_pbl,index_of_temperature,index_of_x_wind, &
52 index_of_y_wind,errmsg,errflg )
91 real(kind=kind_phys),
parameter :: xkzminm = 0.1,xkzminh = 0.01
92 real(kind=kind_phys),
parameter :: xkzmax = 1000.,rimin = -100.
93 real(kind=kind_phys),
parameter :: rlam = 30.,prmin = 0.25,prmax = 4.
94 real(kind=kind_phys),
parameter :: brcr_ub = 0.0,brcr_sb = 0.25,cori = 1.e-4
95 real(kind=kind_phys),
parameter :: afac = 6.8,bfac = 6.8,pfac = 2.0,pfac_q = 2.0
96 real(kind=kind_phys),
parameter :: phifac = 8.,sfcfrac = 0.1
97 real(kind=kind_phys),
parameter :: d1 = 0.02, d2 = 0.05, d3 = 0.001
98 real(kind=kind_phys),
parameter :: h1 = 0.33333335, h2 = 0.6666667
99 real(kind=kind_phys),
parameter :: ckz = 0.001,zfmin = 1.e-8,aphi5 = 5.,aphi16 = 16.
100 real(kind=kind_phys),
parameter :: tmin=1.e-2
101 real(kind=kind_phys),
parameter :: gamcrt = 3.,gamcrq = 2.e-3
102 real(kind=kind_phys),
parameter :: xka = 2.4e-5
103 real(kind=kind_phys),
intent(in) :: karman
104 real(kind=kind_phys),
parameter :: corf=0.000073
105 real(kind=kind_phys),
parameter :: rcl = 1.0
106 integer,
parameter :: imvdif = 1
107 integer,
parameter :: shinhong_tke_diag = 0
111 real(kind=kind_phys),
parameter :: epsq2l = 0.01,c_1 = 1.0,gamcre = 0.224
115 real(kind=kind_phys),
parameter :: mltop = 1.0,sfcfracn1 = 0.075
116 real(kind=kind_phys),
parameter :: nlfrac = 0.7,enlfrac = -0.4
117 real(kind=kind_phys),
parameter :: a11 = 1.0,a12 = -1.15
118 real(kind=kind_phys),
parameter :: ezfac = 1.5
119 real(kind=kind_phys),
parameter :: cpent = -0.4,rigsmax = 100.
120 real(kind=kind_phys),
parameter :: entfmin = 1.0, entfmax = 5.0
122 integer,
intent(in ) :: im,km,ntrac,ndiff,ntcw,ntiw,ntoz
123 real(kind=kind_phys),
intent(in ) :: g,cp,rd,rv,ep1,ep2,xlv,dt
124 logical,
intent(in ) :: lssav, ldiag3d, flag_for_pbl_generic_tend
126 real(kind=kind_phys),
dimension(:,:) , &
127 intent(in ) :: phil, &
133 real(kind=kind_phys),
dimension(:,:,:) , &
136 real(kind=kind_phys),
dimension(:,:) , &
137 intent(in ) :: p2di, &
140 real(kind=kind_phys),
dimension(:,:) , &
141 intent(inout) :: utnp, &
144 real(kind=kind_phys),
dimension(:,:,:) , &
145 intent(inout) :: qtnp
147 integer,
dimension(:) , &
148 intent(in ) :: landmask
150 real(kind=kind_phys),
dimension(:) , &
151 intent(in ) :: heat, &
164 integer,
dimension(:) , &
165 intent(out ) :: kpbl1d
167 real(kind=kind_phys),
dimension(:) , &
168 intent(out ) :: hpbl, &
174 real(kind=kind_phys),
intent(inout),
optional :: dtend(:,:,:)
175 integer,
intent(in) :: dtidx(:,:), index_of_process_pbl, ntqv, &
176 index_of_x_wind, index_of_y_wind, index_of_temperature
182 character(len=*),
intent(out) :: errmsg
183 integer,
intent(out) :: errflg
187 integer :: n,i,k,l,ic
189 integer :: lmh,lmxl,kts,kte,its,ite
191 real(kind=kind_phys) :: dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum,prnum0
192 real(kind=kind_phys) :: ss,ri,qmean,tmean,alpha,chi,zk,rl2,dk,sri
193 real(kind=kind_phys) :: brint,dtodsd,dtodsu,rdz,dsdzt,dsdzq,dsdz2,rlamdz
194 real(kind=kind_phys) :: utend,vtend,ttend,qtend
195 real(kind=kind_phys) :: dtstep,govrthv
196 real(kind=kind_phys) :: cont, conq, conw, conwrc
197 real(kind=kind_phys) :: delxy,pu1,pth1,pq1
198 real(kind=kind_phys) :: dex,hgame_c
199 real(kind=kind_phys) :: zfacdx
200 real(kind=kind_phys) :: amf1,amf2,bmf2,amf3,bmf3,amf4,bmf4,sflux0,snlflux0
201 real(kind=kind_phys) :: mlfrac,ezfrac,sfcfracn
202 real(kind=kind_phys) :: uwst,uwstx,csfac
203 real(kind=kind_phys) :: prnumfac,bfx0,hfx0,qfx0,delb,dux,dvx, &
204 dsdzu,dsdzv,wm3,dthx,dqx,wspd10,ross,tem1,dsig,tvcon,conpr, &
207 integer,
dimension(im) :: kpbl
208 real(kind=kind_phys),
dimension(im) :: hol
209 real(kind=kind_phys),
dimension(im) :: deltaoh
210 real(kind=kind_phys),
dimension(im) :: rigs, &
213 real(kind=kind_phys),
dimension(im) :: &
225 real(kind=kind_phys),
dimension(im) :: &
235 real(kind=kind_phys),
dimension(im) :: &
240 real(kind=kind_phys),
dimension(im) :: &
246 real(kind=kind_phys),
dimension(im,km) :: &
255 real(kind=kind_phys),
dimension(im,km) :: &
263 real(kind=kind_phys),
dimension(im,km) :: &
265 real(kind=kind_phys),
dimension(im,km) :: &
268 real(kind=kind_phys),
dimension(im,km) :: &
272 real(kind=kind_phys),
dimension(im,km) :: &
278 real(kind=kind_phys),
dimension( im, km+1 ) :: zq
279 real(kind=kind_phys),
dimension( im, km, ndiff ) :: r3,f3
281 real(kind=kind_phys),
dimension( km ) :: &
286 real(kind=kind_phys),
dimension( km ) :: &
287 ps1d,pb1d,eps1d,pt1d, &
288 xkze1d,eflx_l1d,eflx_nl1d, &
290 real(kind=kind_phys),
dimension( 2:km ) :: &
293 mfk,ufxpblk,vfxpblk,qfxpblk
294 real(kind=kind_phys),
dimension( km+1 ) :: zqk
296 real(kind=kind_phys),
dimension(im,km) :: dz8w2d
298 logical,
dimension(im) :: pblflg, &
301 logical,
dimension( ndiff ) :: ifvmix
320 conwrc = conw*sqrt(rcl)
321 conpr = bfac*karman*sfcfrac
324 if(landmask(i).eq.0)
then
337 thx(i,k) = tx(i,k)/pi2d(i,k)
343 tvcon = (1.+ep1*qx(i,k,1))
344 thvx(i,k) = thx(i,k)*tvcon
349 tvcon = (1.+ep1*qx(i,1,1))
350 rhox(i) = psfcpa(i)/(rd*tx(i,1)*tvcon)
351 govrth(i) = g/thx(i,1)
352 hfx(i) = heat(i)*rhox(i)*cp
353 qfx(i) = evap(i)*rhox(i)
354 ust(i) = sqrt(stress(i))
355 znt(i) = 0.01*zorl(i)
367 zq(i,k+1) = phii(i,k+1)*conw
368 za(i,k) = phil(i,k)*conw
374 dzq(i,k) = zq(i,k+1)-zq(i,k)
375 del(i,k) = p2di(i,k)-p2di(i,k+1)
386 dza(i,k) = za(i,k)-za(i,k-1)
391 wspd1(i) = sqrt(ux(i,1)*ux(i,1)+vx(i,1)*vx(i,1))+1.e-9
503 hpbl_cbl(i) = zq(i,1)
505 thermal(i)= thvx(i,1)
508 sflux(i) = hfx(i)/rhox(i)/cp + qfx(i)/rhox(i)*ep1*thx(i,1)
509 if(br(i).gt.0.0) sfcflg(i) = .false.
522 if(.not.stable(i))
then
524 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
525 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
527 stable(i) = brup(i).gt.brcr(i)
534 if(brdn(i).ge.brcr(i))
then
536 elseif(brup(i).le.brcr(i))
then
539 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
541 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
542 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
543 if(kpbl(i).le.1) pblflg(i) = .false.
549 zol1(i) = max(br(i)*fm*fm/fh,rimin)
551 zol1(i) = min(zol1(i),-zfmin)
553 zol1(i) = max(zol1(i),zfmin)
555 hol1 = zol1(i)*hpbl(i)/zl1(i)*sfcfrac
558 phim(i) = (1.-aphi16*hol1)**(-1./4.)
559 phih(i) = (1.-aphi16*hol1)**(-1./2.)
560 bfx0 = max(sflux(i),0.)
561 hfx0 = max(hfx(i)/rhox(i)/cp,0.)
562 qfx0 = max(ep1*thx(i,1)*qfx(i)/rhox(i),0.)
563 wstar3(i) = (govrth(i)*bfx0*hpbl(i))
564 wstar(i) = (wstar3(i))**h1
566 phim(i) = (1.+aphi5*hol1)
572 wscale(i) = (ust3(i)+phifac*karman*wstar3(i)*0.5)**h1
573 wscale(i) = min(wscale(i),ust(i)*aphi16)
574 wscale(i) = max(wscale(i),ust(i)/aphi5)
581 if(sfcflg(i).and.sflux(i).gt.0.0)
then
582 gamfac = bfac/rhox(i)/wscale(i)
583 hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt)
584 hgamq(i) = min(gamfac*qfx(i),gamcrq)
585 vpert = (hgamt(i)+ep1*thx(i,1)*hgamq(i))/bfac*afac
586 thermal(i) = thermal(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0)
587 hgamt(i) = max(hgamt(i),0.0)
588 hgamq(i) = max(hgamq(i),0.0)
589 brint = -15.9*ust(i)*ust(i)/wspd(i)*wstar3(i)/(wscale(i)**4.)
590 hgamu(i) = brint*ux(i,1)
591 hgamv(i) = brint*vx(i,1)
616 if(.not.stable(i).and.pblflg(i))
then
618 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
619 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
621 stable(i) = brup(i).gt.brcr(i)
629 if(brdn(i).ge.brcr(i))
then
631 elseif(brup(i).le.brcr(i))
then
634 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
636 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
637 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
638 if(kpbl(i).le.1) pblflg(i) = .false.
639 uwst = abs(ust(i)/wstar(i)-0.5)
640 uwstx = -80.*uwst+14.
641 csfac = 0.5*(tanh(uwstx)+3.)
642 cslen(i) = csfac*hpbl(i)
649 hpbl_cbl(i) = hpbl(i)
650 if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2))
then
659 if((.not.stable(i)).and.((xland(i)-1.5).ge.0))
then
660 wspd10 = u10(i)*u10(i) + v10(i)*v10(i)
661 wspd10 = sqrt(wspd10)
662 ross = wspd10 / (cori*znt(i))
663 brcr_sbro(i) = min(0.16*(1.e-7*ross)**(-0.18),.3)
668 if(.not.stable(i))
then
669 if((xland(i)-1.5).ge.0)
then
670 brcr(i) = brcr_sbro(i)
679 if(.not.stable(i))
then
681 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
682 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
684 stable(i) = brup(i).gt.brcr(i)
690 if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2))
then
692 if(brdn(i).ge.brcr(i))
then
694 elseif(brup(i).le.brcr(i))
then
697 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
699 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
700 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
701 if(kpbl(i).le.1) pblflg(i) = .false.
708 pu1=pu(dx(i),cslen(i))
709 pq1=pq(dx(i),cslen(i))
711 hgamu(i) = hgamu(i)*pu1
712 hgamv(i) = hgamv(i)*pu1
713 hgamq(i) = hgamq(i)*pq1
723 wm3 = wstar3(i) + 5. * ust3(i)
725 bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i)
726 dthvx(i) = max(thvx(i,k+1)-thvx(i,k),tmin)
727 dthx = max(thx(i,k+1)-thx(i,k),tmin)
728 dqx = min(qx(i,k+1,1)-qx(i,k,1),0.0)
729 we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i)))
730 hfxpbl(i) = we(i)*dthx
731 pq1=pq(dx(i),cslen(i))
732 qfxpbl(i) = we(i)*dqx*pq1
734 pu1=pu(dx(i),cslen(i))
735 dux = ux(i,k+1)-ux(i,k)
736 dvx = vx(i,k+1)-vx(i,k)
738 ufxpbl(i) = max(prpbl(i)*we(i)*dux*pu1,-ust(i)*ust(i))
739 elseif(dux.lt.-tmin)
then
740 ufxpbl(i) = min(prpbl(i)*we(i)*dux*pu1,ust(i)*ust(i))
745 vfxpbl(i) = max(prpbl(i)*we(i)*dvx*pu1,-ust(i)*ust(i))
746 elseif(dvx.lt.-tmin)
then
747 vfxpbl(i) = min(prpbl(i)*we(i)*dvx*pu1,ust(i)*ust(i))
751 delb = govrth(i)*d3*hpbl(i)
752 delta(i) = min(d1*hpbl(i) + d2*wm2(i)/delb,100.)
753 delb = govrth(i)*dthvx(i)
754 deltaoh(i) = d1*hpbl(i) + d2*wm2(i)/delb
755 deltaoh(i) = max(ezfac*deltaoh(i),hpbl(i)-za(i,kpbl(i)-1)-1.)
756 deltaoh(i) = min(deltaoh(i) ,hpbl(i))
757 rigs(i) = govrth(i)*dthvx(i)*deltaoh(i)/(dux**2.+dvx**2.)
758 rigs(i) = max(min(rigs(i), rigsmax),rimin)
759 enlfrac2(i) = max(min(wm3/wstar3(i)/(1.+cpent/rigs(i)),entfmax), entfmin)
760 enlfrac2(i) = enlfrac2(i)*enlfrac
767 entfacmf(i,k) = sqrt(((zq(i,k+1)-hpbl(i))/deltaoh(i))**2.)
769 if(pblflg(i).and.k.ge.kpbl(i))
then
770 entfac(i,k) = ((zq(i,k+1)-hpbl(i))/deltaoh(i))**2.
781 if(k.lt.kpbl(i))
then
782 zfac(i,k) = min(max((1.-(zq(i,k+1)-zl1(i))/(hpbl(i)-zl1(i))),zfmin),1.)
783 zfacent(i,k) = (1.-zfac(i,k))**3.
784 wscalek(i,k) = (ust3(i)+phifac*karman*wstar3(i)*(1.-zfac(i,k)))**h1
787 prfac2 = 15.9*wstar3(i)/ust3(i)/(1.+4.*karman*wstar3(i)/ust3(i))
788 prnumfac = -3.*(max(zq(i,k+1)-sfcfrac*hpbl(i),0.))**2./hpbl(i)**2.
793 phim8z = 1.+aphi5*zol1(i)*zq(i,k+1)/zl1(i)
794 wscalek(i,k) = ust(i)/phim8z
795 wscalek(i,k) = max(wscalek(i,k),0.001)
797 prnum0 = (phih(i)/phim(i)+prfac)
798 prnum0 = max(min(prnum0,prmax),prmin)
799 xkzm(i,k) = wscalek(i,k)*karman*zq(i,k+1)*zfac(i,k)**pfac
800 prnum = 1. + (prnum0-1.)*exp(prnumfac)
801 xkzq(i,k) = xkzm(i,k)/prnum*zfac(i,k)**(pfac_q-pfac)
802 prnum0 = prnum0/(1.+prfac2*karman*sfcfrac)
803 prnum = 1. + (prnum0-1.)*exp(prnumfac)
804 xkzh(i,k) = xkzm(i,k)/prnum
805 xkzm(i,k) = xkzm(i,k)+xkzom(i,k)
806 xkzh(i,k) = xkzh(i,k)+xkzoh(i,k)
807 xkzq(i,k) = xkzq(i,k)+xkzoh(i,k)
808 xkzm(i,k) = min(xkzm(i,k),xkzmax)
809 xkzh(i,k) = min(xkzh(i,k),xkzmax)
810 xkzq(i,k) = min(xkzq(i,k),xkzmax)
819 if(k.ge.kpbl(i))
then
820 ss = ((ux(i,k+1)-ux(i,k))*(ux(i,k+1)-ux(i,k)) &
821 +(vx(i,k+1)-vx(i,k))*(vx(i,k+1)-vx(i,k))) &
822 /(dza(i,k+1)*dza(i,k+1))+1.e-9
823 govrthv = g/(0.5*(thvx(i,k+1)+thvx(i,k)))
824 ri = govrthv*(thvx(i,k+1)-thvx(i,k))/(ss*dza(i,k+1))
826 if(imvdif.eq.1.and.ntcw.ge.2.and.ntiw.ge.2)
then
827 if((qx(i,k,ntcw)+qx(i,k,ntiw)).gt.0.01e-3 &
828 .and.(qx(i,k+1,ntcw)+qx(i,k+1,ntiw)).gt.0.01e-3)
then
829 qmean = 0.5*(qx(i,k,1)+qx(i,k+1,1))
830 tmean = 0.5*(tx(i,k)+tx(i,k+1))
831 alpha = xlv*qmean/rd/tmean
832 chi = xlv*xlv*qmean/cp/rv/tmean/tmean
833 ri = (1.+alpha)*(ri-g*g/ss/tmean/cp*((chi-alpha)/(1.+chi)))
836 zk = karman*zq(i,k+1)
837 rlamdz = min(max(0.1*dza(i,k+1),rlam),300.)
838 rlamdz = min(dza(i,k+1),rlamdz)
839 rl2 = (zk*rlamdz/(rlamdz+zk))**2
845 xkzm(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri))
846 xkzh(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri))
849 xkzh(i,k) = dk/(1+5.*ri)**2
851 prnum = min(prnum,prmax)
852 xkzm(i,k) = xkzh(i,k)*prnum
855 xkzm(i,k) = xkzm(i,k)+xkzom(i,k)
856 xkzh(i,k) = xkzh(i,k)+xkzoh(i,k)
857 xkzm(i,k) = min(xkzm(i,k),xkzmax)
858 xkzh(i,k) = min(xkzh(i,k),xkzmax)
859 xkzml(i,k) = xkzm(i,k)
860 xkzhl(i,k) = xkzh(i,k)
868 deltaoh(i) = deltaoh(i)/hpbl(i)
872 mlfrac = mltop-deltaoh(i)
873 ezfrac = mltop+deltaoh(i)
874 zfacmf(i,1) = min(max((zq(i,2)/hpbl(i)),zfmin),1.)
875 sfcfracn = max(sfcfracn1,zfacmf(i,1))
877 sflux0 = (a11+a12*sfcfracn)*sflux(i)
878 snlflux0 = nlfrac*sflux0
879 amf1 = snlflux0/sfcfracn
880 amf2 = -snlflux0/(mlfrac-sfcfracn)
882 amf3 = snlflux0*enlfrac2(i)/deltaoh(i)
884 hfxpbl(i) = amf3+bmf3
885 pth1=pthnl(dx(i),cslen(i))
886 hfxpbl(i) = hfxpbl(i)*pth1
889 zfacmf(i,k) = max((zq(i,k+1)/hpbl(i)),zfmin)
890 if(pblflg(i).and.k.lt.kpbl(i))
then
891 if(zfacmf(i,k).le.sfcfracn)
then
892 mf(i,k) = amf1*zfacmf(i,k)
893 else if (zfacmf(i,k).le.mlfrac)
then
894 mf(i,k) = amf2*zfacmf(i,k)+bmf2
896 mf(i,k) = mf(i,k)+hfxpbl(i)*exp(-entfacmf(i,k))
897 mf(i,k) = mf(i,k)*pth1
915 f1(i,1) = thx(i,1)-300.+hfx(i)/cont/del(i,1)*dt2
920 dtodsd = dt2/del(i,k)
921 dtodsu = dt2/del(i,k+1)
922 dsig = p2d(i,k)-p2d(i,k+1)
924 tem1 = dsig*xkzh(i,k)*rdz
925 if(pblflg(i).and.k.lt.kpbl(i))
then
926 dsdzt = tem1*(-mf(i,k)/xkzh(i,k))
927 f1(i,k) = f1(i,k)+dtodsd*dsdzt
928 f1(i,k+1) = thx(i,k+1)-300.-dtodsu*dsdzt
929 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6)
then
930 xkzh(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
931 xkzh(i,k) = sqrt(xkzh(i,k)*xkzhl(i,k))
932 xkzh(i,k) = max(xkzh(i,k),xkzoh(i,k))
933 xkzh(i,k) = min(xkzh(i,k),xkzmax)
934 f1(i,k+1) = thx(i,k+1)-300.
936 f1(i,k+1) = thx(i,k+1)-300.
938 tem1 = dsig*xkzh(i,k)*rdz
940 au(i,k) = -dtodsd*dsdz2
941 al(i,k) = -dtodsu*dsdz2
945 zfacdx=0.2*hpbl(i)/zq(i,k+1)
946 delxy=dx(i)*max(zfacdx,1.0)
947 pth1=pthl(delxy,hpbl(i))
948 if(pblflg(i).and.k.lt.kpbl(i))
then
949 au(i,k) = au(i,k)*pth1
950 al(i,k) = al(i,k)*pth1
952 ad(i,k) = ad(i,k)-au(i,k)
953 ad(i,k+1) = 1.-al(i,k)
966 call tridin_ysu(al,ad,cu,r1,au,f1,its,ite,kts,kte,1)
972 ttend = (f1(i,k)-thx(i,k)+300.)*rdt*pi2d(i,k)
973 ttnp(i,k) = ttnp(i,k)+ttend
974 dtsfc(i) = dtsfc(i)+ttend*cont*del(i,k)
976 tflux_e(i,k) = ttend*dz8w2d(i,k)
978 tflux_e(i,k) = tflux_e(i,k+1) + ttend*dz8w2d(i,k)
982 if(lssav .and. ldiag3d .and. .not. flag_for_pbl_generic_tend)
then
983 idtend = dtidx(index_of_temperature,index_of_process_pbl)
985 dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f1-thx+300.)*rdt*pi2d
1009 f3(i,1,1) = qx(i,1,1)+qfx(i)*g/del(i,1)*dt2
1015 f3(i,1,ic) = qx(i,1,ic)
1022 if(k.ge.kpbl(i))
then
1023 xkzq(i,k) = xkzh(i,k)
1030 dtodsd = dt2/del(i,k)
1031 dtodsu = dt2/del(i,k+1)
1032 dsig = p2d(i,k)-p2d(i,k+1)
1034 tem1 = dsig*xkzq(i,k)*rdz
1035 if(pblflg(i).and.k.lt.kpbl(i))
then
1036 dsdzq = tem1*(-qfxpbl(i)*zfacent(i,k)/xkzq(i,k))
1037 f3(i,k,1) = f3(i,k,1)+dtodsd*dsdzq
1038 f3(i,k+1,1) = qx(i,k+1,1)-dtodsu*dsdzq
1039 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6)
then
1040 xkzq(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
1041 xkzq(i,k) = sqrt(xkzq(i,k)*xkzhl(i,k))
1042 xkzq(i,k) = max(xkzq(i,k),xkzoh(i,k))
1043 xkzq(i,k) = min(xkzq(i,k),xkzmax)
1044 f3(i,k+1,1) = qx(i,k+1,1)
1046 f3(i,k+1,1) = qx(i,k+1,1)
1048 tem1 = dsig*xkzq(i,k)*rdz
1050 au(i,k) = -dtodsd*dsdz2
1051 al(i,k) = -dtodsu*dsdz2
1055 zfacdx=0.2*hpbl(i)/zq(i,k+1)
1056 delxy=dx(i)*max(zfacdx,1.0)
1057 pq1=pq(delxy,hpbl(i))
1058 if(pblflg(i).and.k.lt.kpbl(i))
then
1059 au(i,k) = au(i,k)*pq1
1060 al(i,k) = al(i,k)*pq1
1062 ad(i,k) = ad(i,k)-au(i,k)
1063 ad(i,k+1) = 1.-al(i,k)
1071 f3(i,k+1,ic) = qx(i,k+1,ic)
1088 r3(i,k,ic) = f3(i,k,ic)
1095 call tridin_ysu(al,ad,cu,r3,au,f3,its,ite,kts,kte,ndiff)
1101 qtend = (f3(i,k,1)-qx(i,k,1))*rdt
1102 qtnp(i,k,1) = qtnp(i,k,1)+qtend
1103 dqsfc(i) = dqsfc(i)+qtend*conq*del(i,k)
1105 qflux_e(i,k) = qtend*dz8w2d(i,k)
1107 qflux_e(i,k) = qflux_e(i,k+1) + qtend*dz8w2d(i,k)
1109 tvflux_e(i,k) = tflux_e(i,k) + qflux_e(i,k)*ep1*thx(i,k)
1112 if(lssav .and. ldiag3d .and. .not. flag_for_pbl_generic_tend)
then
1113 idtend = dtidx(ntqv+100,index_of_process_pbl)
1115 dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*rdt*(f3(:,:,1)-qx(:,:,1))
1122 if(pblflg(i).and.k.lt.kpbl(i))
then
1123 hgame_c=c_1*0.2*2.5*(g/thvx(i,k))*wstar(i)/(0.25*(q2x(i,k+1)+q2x(i,k)))
1124 hgame_c=min(hgame_c,gamcre)
1126 hgame2d(i,k)=hgame_c*0.5*tvflux_e(i,k)*hpbl(i)
1127 hgame2d(i,k)=max(hgame2d(i,k),0.0)
1129 hgame2d(i,k)=hgame_c*0.5*(tvflux_e(i,k)+tvflux_e(i,k+1))*hpbl(i)
1130 hgame2d(i,k)=max(hgame2d(i,k),0.0)
1141 qtend = (f3(i,k,ic)-qx(i,k,ic))*rdt
1142 qtnp(i,k,ic) = qtnp(i,k,ic)+qtend
1147 if(lssav .and. ldiag3d .and. ntoz>0 .and. &
1148 & .not. flag_for_pbl_generic_tend)
then
1149 idtend=dtidx(ntoz+100,index_of_process_pbl)
1151 dtend(:,:,idtend) = dtend(:,:,idtend) + qtend*(f3(:,:,ntoz)-qx(:,:,ntoz))
1169 ad(i,1) = 1.+ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2 &
1170 *(wspd1(i)/wspd(i))**2
1177 dtodsd = dt2/del(i,k)
1178 dtodsu = dt2/del(i,k+1)
1179 dsig = p2d(i,k)-p2d(i,k+1)
1181 tem1 = dsig*xkzm(i,k)*rdz
1182 if(pblflg(i).and.k.lt.kpbl(i))
then
1183 dsdzu = tem1*(-hgamu(i)/hpbl(i)-ufxpbl(i)*zfacent(i,k)/xkzm(i,k))
1184 dsdzv = tem1*(-hgamv(i)/hpbl(i)-vfxpbl(i)*zfacent(i,k)/xkzm(i,k))
1185 f1(i,k) = f1(i,k)+dtodsd*dsdzu
1186 f1(i,k+1) = ux(i,k+1)-dtodsu*dsdzu
1187 f2(i,k) = f2(i,k)+dtodsd*dsdzv
1188 f2(i,k+1) = vx(i,k+1)-dtodsu*dsdzv
1189 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6)
then
1190 xkzm(i,k) = prpbl(i)*xkzh(i,k)
1191 xkzm(i,k) = sqrt(xkzm(i,k)*xkzml(i,k))
1192 xkzm(i,k) = max(xkzm(i,k),xkzom(i,k))
1193 xkzm(i,k) = min(xkzm(i,k),xkzmax)
1194 f1(i,k+1) = ux(i,k+1)
1195 f2(i,k+1) = vx(i,k+1)
1197 f1(i,k+1) = ux(i,k+1)
1198 f2(i,k+1) = vx(i,k+1)
1200 tem1 = dsig*xkzm(i,k)*rdz
1202 au(i,k) = -dtodsd*dsdz2
1203 al(i,k) = -dtodsu*dsdz2
1207 zfacdx=0.2*hpbl(i)/zq(i,k+1)
1208 delxy=dx(i)*max(zfacdx,1.0)
1209 pu1=pu(delxy,hpbl(i))
1210 if(pblflg(i).and.k.lt.kpbl(i))
then
1211 au(i,k) = au(i,k)*pu1
1212 al(i,k) = al(i,k)*pu1
1214 ad(i,k) = ad(i,k)-au(i,k)
1215 ad(i,k+1) = 1.-al(i,k)
1231 call tridi1n(al,ad,cu,r1,r2,au,f1,f2,its,ite,kts,kte,1)
1237 utend = (f1(i,k)-ux(i,k))*rdt
1238 vtend = (f2(i,k)-vx(i,k))*rdt
1239 utnp(i,k) = utnp(i,k)+utend
1240 vtnp(i,k) = vtnp(i,k)+vtend
1241 dusfc(i) = dusfc(i) + utend*conwrc*del(i,k)
1242 dvsfc(i) = dvsfc(i) + vtend*conwrc*del(i,k)
1245 if(lssav .and. ldiag3d .and. .not. flag_for_pbl_generic_tend)
then
1246 idtend=dtidx(index_of_x_wind,index_of_process_pbl)
1248 dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*rdt*(f1-ux)
1250 idtend=dtidx(index_of_y_wind,index_of_process_pbl)
1252 dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*rdt*(f2-vx)
1262 if (shinhong_tke_diag.eq.1)
then
1264 tke_calculation:
do i = its,ite
1305 thvxk(k) = thvx(i,k)
1307 hgame(k) = hgame2d(i,k)
1311 if(pblflg(i).and.k.le.kpbl(i))
then
1312 zfacdx = 0.2*hpbl(i)/za(i,k)
1313 delxy = dx(i)*max(zfacdx,1.0)
1314 ptke1(k+1) = ptke(delxy,hpbl(i))
1323 akmk(k) = xkzm(i,k-1)
1324 akhk(k) = xkzh(i,k-1)
1325 mfk(k) = mf(i,k-1)/xkzh(i,k-1)
1326 ufxpblk(k) = ufxpbl(i)*zfacent(i,k-1)/xkzm(i,k-1)
1327 vfxpblk(k) = vfxpbl(i)*zfacent(i,k-1)/xkzm(i,k-1)
1328 qfxpblk(k) = qfxpbl(i)*zfacent(i,k-1)/xkzq(i,k-1)
1333 dex = 0.25*(q2xk(k+2)-q2xk(k))
1334 efxpbl(i) = we(i)*dex
1339 call mixlen(lmh,uxk,vxk,txk,thxk,qx(i,kts:kte,1),qx(i,kts:kte,ntcw) &
1340 ,q2xk,zqk,ust(i),corf,epshol(i) &
1342 ,hpbl(i),kpbl(i),lmxl,ct(i) &
1343 ,hgamu(i),hgamv(i),hgamq(i),pblflg(i) &
1344 ,mfk,ufxpblk,vfxpblk,qfxpblk &
1350 call prodq2(lmh,dt,ust(i),s2,rig,q2xk,el,zqk,akmk,akhk &
1351 ,uxk,vxk,thxk,thvxk &
1352 ,hgamu(i),hgamv(i),hgamq(i),dx(i) &
1353 ,hpbl(i),pblflg(i),kpbl(i) &
1354 ,mfk,ufxpblk,vfxpblk,qfxpblk &
1361 call vdifq(lmh,dt,q2xk,el,zqk &
1363 ,hgame,hpbl(i),pblflg(i),kpbl(i) &
1370 q2x(i,k) = amax1(q2xk(k),epsq2l)
1373 enddo tke_calculation