69 subroutine hedmf_run (im,km,ntrac,ntcw,dv,du,tau,rtg, &
70 & u1,v1,t1,q1,swh,hlw,xmu, &
71 & psk,rbsoil,zorl,u10m,v10m,fm,fh, &
72 & tsea,heat,evap,stress,spd1,kpbl, &
73 & prsi,del,prsl,prslk,phii,phil,delt,dspheat, &
74 & dusfc,dvsfc,dtsfc,dqsfc,hpbl,hgamt,hgamq,dkt,dku, &
75 & kinver,xkzm_m,xkzm_h,xkzm_s,lprnt,ipr, &
76 & xkzminv,moninq_fac,hurr_pbl,islimsk,var_ric, &
77 & coef_ric_l,coef_ric_s,ldiag3d,ntqv,rtg_ozone_index,ntoz, &
78 & dtend,dtidx,index_of_process_pbl,index_of_x_wind, &
79 & index_of_y_wind,index_of_temperature, &
80 & flag_for_pbl_generic_tend,errmsg,errflg)
87 use physcons, grav => con_g, cp => con_cp,
88 & hvap => con_hvap, fv => con_fvirt
94 logical,
intent(in) :: lprnt, hurr_pbl, ldiag3d
95 logical,
intent(in) :: flag_for_pbl_generic_tend
96 integer,
intent(in) :: ipr, islimsk(:)
97 integer,
intent(in) :: im, km, ntrac, ntcw, kinver(:), ntoz
98 integer,
intent(out) :: kpbl(:)
101 real(kind=kind_phys),
intent(in) :: delt, xkzm_m, xkzm_h, xkzm_s
102 real(kind=kind_phys),
intent(in) :: xkzminv, moninq_fac, var_ric, &
103 & coef_ric_l, coef_ric_s
104 real(kind=kind_phys),
intent(inout) :: dv(:,:), du(:,:), &
105 & tau(:,:), rtg(:,:,:)
107 real(kind=kind_phys),
intent(inout),
optional :: dtend(:,:,:)
108 integer,
intent(in) :: dtidx(:,:)
109 integer,
intent(in) :: index_of_x_wind, index_of_y_wind, &
110 & index_of_process_pbl, index_of_temperature, ntqv, rtg_ozone_index
111 real(kind=kind_phys),
intent(in) :: &
112 & u1(:,:), v1(:,:), &
113 & t1(:,:), q1(:,:,:), &
114 & swh(:,:), hlw(:,:), &
116 & rbsoil(:), zorl(:), &
117 & u10m(:), v10m(:), &
120 & heat(:), evap(:), &
122 real(kind=kind_phys),
intent(in) :: &
123 & prsi(:,:), del(:,:), &
124 & prsl(:,:), prslk(:,:), &
125 & phii(:,:), phil(:,:)
126 real(kind=kind_phys),
intent(out) :: &
127 & dusfc(:), dvsfc(:), &
128 & dtsfc(:), dqsfc(:), &
130 real(kind=kind_phys),
intent(out) :: &
132 real(kind=kind_phys),
intent(inout) :: &
135 logical,
intent(in) :: dspheat
137 character(len=*),
intent(out) :: errmsg
138 integer,
intent(out) :: errflg
143 integer i,iprt,is,iun,k,kk,km1,kmpbl,latd,lond
144 integer lcld(im),icld(im),kcld(im),krad(im)
145 integer kx1(im), kpblx(im)
148 real(kind=kind_phys) phih(im), phim(im), hpblx(im), &
149 & rbdn(im), rbup(im), &
150 & beta(im), sflux(im), &
151 & z0(im), crb(im), wstar(im), &
152 & zol(im), ustmin(im), ustar(im), &
153 & thermal(im),wscale(im), wscaleu(im)
155 real(kind=kind_phys) theta(im,km),thvx(im,km), thlvx(im,km), &
156 & qlx(im,km), thetae(im,km), &
157 & qtx(im,km), bf(im,km-1), diss(im,km), &
159 & govrth(im), hrad(im), &
161 & radmin(im), vrad(im), &
162 & zd(im), zdd(im), thlvx1(im)
164 real(kind=kind_phys) rdzt(im,km-1),dktx(im,km-1), &
165 & zi(im,km+1), zl(im,km), &
166 & xkzo(im,km-1), xkzmo(im,km-1), &
167 & cku(im,km-1), ckt(im,km-1), &
168 & ti(im,km-1), shr2(im,km-1), &
169 & al(im,km-1), ad(im,km), &
170 & au(im,km-1), a1(im,km), &
173 real(kind=kind_phys) tcko(im,km), qcko(im,km,ntrac), &
174 & ucko(im,km), vcko(im,km), xmf(im,km)
176 real(kind=kind_phys) prinv(im), rent(im)
178 logical pblflg(im), sfcflg(im), scuflg(im), flg(im)
179 logical ublflg(im), pcnvflg(im)
184 real(kind=kind_phys) aphi16, aphi5, bvf2, wfac,
185 & cfac, conq, cont, conw,
187 & dq1, dsdz2, dsdzq, dsdzt,
189 & dsig, dt2, dthe1, dtodsd,
190 & dtodsu, dw2, dw2min,
191 & gamcrq, gamcrt, gocp,
193 & prnum, prmax, prmin, pfac, crbcon,
194 & qmin, tdzmin, qtend, crbmin,crbmax,
195 & rbint, rdt, rdz, qlmin,
196 & ri, rimin, rl2, rlam, rlamun,
197 & rone, rzero, sfcfrac,
198 & spdk2, sri, zol1, zolcr, zolcru,
202 & vtend, zfac, vpert, cteit,
203 & rentf1, rentf2, radfac,
204 & zfmin, zk, tem, tem1, tem2,
206 & ptem, ptem1, ptem2, tx1(im), tx2(im)
208 real(kind=kind_phys) zstblmax,h1, h2, qlcr, actei,
212 integer :: idtend1, idtend2
215 real(kind=kind_phys) wspm(im,km-1)
219 integer,
parameter :: useshape=2
220 real :: smax,ashape,sz2h, sksfc,skmax,ashape1,skminusk0, hmax
222 parameter(gravi=1.0/grav)
223 parameter(gocp=grav/cp)
224 parameter(cont=cp/grav,conq=hvap/grav,conw=1.0/grav)
226 parameter(rlam=30.0,vk=0.4,vk2=vk*vk)
227 parameter(prmin=0.25,prmax=4.,zolcr=0.2,zolcru=-0.5)
228 parameter(dw2min=0.0001,dkmin=0.0,dkmax=1000.,rimin=-100.)
229 parameter(crbcon=0.25,crbmin=0.15,crbmax=0.35)
230 parameter(wfac=7.0,cfac=6.5,pfac=2.0,sfcfrac=0.1)
232 parameter(qmin=1.e-8, zfmin=1.e-8,aphi5=5.,aphi16=16.)
233 parameter(tdzmin=1.e-3,qlmin=1.e-12,f0=1.e-4)
234 parameter(h1=0.33333333,h2=0.66666667)
236 parameter(cldtime=500.)
239 parameter(gamcrt=3.,gamcrq=0.,rlamun=150.0)
240 parameter(rentf1=0.2,rentf2=1.0,radfac=0.85)
247 parameter(zstblmax = 2500., qlcr=3.5e-5)
249 parameter(actei = 0.7)
251c-----------------------------------------------------------------------
253 601
format(1x,
' moninp lat lon step hour ',3i6,f6.1)
254 602
format(1x,
' k',
' z',
' t',
' th',
255 1
' tvh',
' q',
' u',
' v',
257 603
format(1x,i5,8f9.1)
258 604
format(1x,
' sfc',9x,f9.1,18x,f9.1)
259 605
format(1x,
' k zl spd2 thekv the1v'
261 606
format(1x,i5,6f8.2)
262 607
format(1x,
' kpbl hpbl fm fh hgamt',
263 1
' hgamq ws ustar cd ch')
264 608
format(1x,i5,9f8.2)
265 609
format(1x,
' k pr dkt dku ',i5,3f8.2)
266 610
format(1x,
' k pr dkt dku ',i5,3f8.2,
' l2 ri t2',
267 1
' sr2 ',2f8.2,2e10.2)
295 zi(i,k) = phii(i,k) * gravi
296 zl(i,k) = phil(i,k) * gravi
300 zi(i,km+1) = phii(i,km+1) * gravi
305 rdzt(i,k) = 1.0 / (zl(i,k+1) - zl(i,k))
311 tx1(i) = 1.0 / prsi(i,1)
319 if (k < kinver(i))
then
321 ptem = prsi(i,k+1) * tx1(i)
323 tem1 = tem1 * tem1 * 10.0
324 xkzo(i,k) = xkzm_h * min(1.0, exp(-tem1))
327 if (ptem >= xkzm_s)
then
331 if (k == kx1(i) .and. k > 1) tx2(i) = 1.0 / prsi(i,k)
332 tem1 = 1.0 - prsi(i,k+1) * tx2(i)
333 tem1 = tem1 * tem1 * 5.0
334 xkzmo(i,k) = xkzm_m * min(1.0, exp(-tem1))
349 if(zi(i,k+1) > 250.)
then
350 tem1 = (t1(i,k+1)-t1(i,k)) * rdzt(i,k)
351 if(tem1 > 1.e-5)
then
352 xkzo(i,k) = min(xkzo(i,k),xkzminv)
359 z0(i) = 0.01 * zorl(i)
371 if(rbsoil(i) > 0.) sfcflg(i) = .false.
390 theta(i,k) = t1(i,k) * psk(i) / prslk(i,k)
391 qlx(i,k) = max(q1(i,k,ntcw),qlmin)
392 qtx(i,k) = max(q1(i,k,1),qmin)+qlx(i,k)
394 ptem1 = hvap*max(q1(i,k,1),qmin)/(cp*t1(i,k))
395 thetae(i,k)= theta(i,k)*(1.+ptem1)
396 thvx(i,k) = theta(i,k)*(1.+fv*max(q1(i,k,1),qmin)-ptem)
397 ptem2 = theta(i,k)-(hvap/cp)*ptem
398 thlvx(i,k) = ptem2*(1.+fv*qtx(i,k))
413 tem = zi(i,k+1)-zi(i,k)
414 radx(i,k) = tem*(swh(i,k)*xmu(i)+hlw(i,k))
423 if(flg(i).and.zl(i,k) >= zstblmax)
then
435 bf(i,k) = (thvx(i,k+1)-thvx(i,k))*rdz
436 ti(i,k) = 2./(t1(i,k)+t1(i,k+1))
437 dw2 = (u1(i,k)-u1(i,k+1))**2
438 & + (v1(i,k)-v1(i,k+1))**2
439 shr2(i,k) = max(dw2,dw2min)*rdz*rdz
444 govrth(i) = grav/theta(i,1)
448 beta(i) = dt2 / (zi(i,2)-zi(i,1))
452 ustar(i) = sqrt(stress(i))
456 sflux(i) = heat(i) + evap(i)*fv*theta(i,1)
457 if(.not.sfcflg(i) .or. sflux(i) <= 0.) pblflg(i)=.false.
465 if (.not. (hurr_pbl .and. moninq_fac < 0.0))
then
471 thermal(i) = thvx(i,1)
474 thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin))
475 tem = sqrt(u10m(i)**2+v10m(i)**2)
477 robn = tem / (f0 * z0(i))
479 crb(i) = 0.16 * (tem1 ** (-0.18))
480 crb(i) = max(min(crb(i), crbmax), crbmin)
490 thermal(i) = thvx(i,1)
492 thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin))
494 tem = sqrt(u10m(i)**2+v10m(i)**2)
496 robn = tem / (f0 * z0(i))
499 if (var_ric .eq. 1.)
then
500 if (islimsk(i) .eq. 1) crb(i) = coef_ric_l*(tem1)**(-0.18)
501 if (islimsk(i) .eq. 0) crb(i) = coef_ric_s*(tem1)**(-0.18)
503 crb(i) = max(min(crb(i), crbmax), crbmin)
519 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
520 rbup(i) = (thvx(i,k)-thermal(i))*
521 & (grav*zl(i,k)/thvx(i,1))/spdk2
523 flg(i) = rbup(i) > crb(i)
531 if(rbdn(i) >= crb(i))
then
533 elseif(rbup(i) <= crb(i))
then
536 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
538 hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
539 if(hpbl(i) < zi(i,kpbl(i))) kpbl(i) = kpbl(i) - 1
565 zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin)
567 zol(i) = min(zol(i),-zfmin)
569 zol(i) = max(zol(i),zfmin)
571 zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1)
575 tem = 1.0 / (1. - aphi16*zol1)
577 phim(i) = sqrt(phih(i))
579 phim(i) = 1. + aphi5*zol1
582 wscale(i) = ustar(i)/phim(i)
583 ustmin(i) = ustar(i)/aphi5
584 wscale(i) = max(wscale(i),ustmin(i))
588 if(zol(i) < zolcru .and. kpbl(i) > 1)
then
593 wst3 = govrth(i)*sflux(i)*hpbl(i)
596 wscaleu(i) = (ust3+wfac*vk*wst3*sfcfrac)**h1
597 wscaleu(i) = max(wscaleu(i),ustmin(i))
606 hgamt(i) = min(cfac*heat(i)/wscaleu(i),gamcrt)
607 hgamq(i) = min(cfac*evap(i)/wscaleu(i),gamcrq)
608 vpert = hgamt(i) + hgamq(i)*fv*theta(i,1)
609 vpert = min(vpert,gamcrt)
610 thermal(i)= thermal(i)+max(vpert,0.)
611 hgamt(i) = max(hgamt(i),0.0)
612 hgamq(i) = max(hgamq(i),0.0)
629 spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
630 rbup(i) = (thvx(i,k)-thermal(i))*
631 & (grav*zl(i,k)/thvx(i,1))/spdk2
633 flg(i) = rbup(i) > crb(i)
640 if(rbdn(i) >= crb(i))
then
642 elseif(rbup(i) <= crb(i))
then
645 rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i))
647 hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1))
648 if(hpbl(i) < zi(i,kpbl(i))) kpbl(i) = kpbl(i) - 1
649 if(kpbl(i) <= 1)
then
664 if(flg(i) .and. k <= lcld(i))
then
665 if(qlx(i,k).ge.qlcr)
then
673 if(scuflg(i) .and. kcld(i)==km1) scuflg(i)=.false.
681 if(flg(i) .and. k <= kcld(i))
then
682 if(qlx(i,k) >= qlcr)
then
683 if(radx(i,k) < radmin(i))
then
694 if(scuflg(i) .and. krad(i) <= 1) scuflg(i)=.false.
695 if(scuflg(i) .and. radmin(i)>=0.) scuflg(i)=.false.
703 if(flg(i) .and. k <= krad(i))
then
704 if(qlx(i,k) >= qlcr)
then
713 if(scuflg(i) .and. icld(i) < 1) scuflg(i)=.false.
718 hrad(i) = zi(i,krad(i)+1)
724 if(scuflg(i) .and. hrad(i)<zi(i,2)) scuflg(i)=.false.
730 tem = zi(i,k+1)-zi(i,k)
731 tem1 = cldtime*radmin(i)/tem
732 thlvx1(i) = thlvx(i,k)+tem1
742 if(flg(i) .and. k <= krad(i))
then
743 if(thlvx1(i) <= thlvx(i,k))
then
744 tem=zi(i,k+1)-zi(i,k)
755 kk = max(1, krad(i)+1-icld(i))
756 zdd(i) = hrad(i)-zi(i,kk)
762 zd(i) = max(zd(i),zdd(i))
763 zd(i) = min(zd(i),hrad(i))
764 tem = govrth(i)*zd(i)*(-radmin(i))
774 tem = phih(i)/phim(i)+cfac*vk*sfcfrac
776 tem = phih(i)/phim(i)
779 prinv(i) = min(prinv(i),prmax)
780 prinv(i) = max(prinv(i),prmin)
783 if(zol(i) > zolcr)
then
797 if(hurr_pbl .and. moninq_fac < 0.0)
then
802 if (zi(i,k) .le. 500. .and. zi(i,k+1) .gt. 500.)
then
803 spdk2 = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k))
804 wspm(i,1) = spdk2/0.6
815 if (.not. (hurr_pbl .and. moninq_fac < 0.0))
then
821 zfac = max((1.-zi(i,k+1)/hpbl(i)), zfmin)
822 tem = zi(i,k+1) * (zfac**pfac) * moninq_fac
824 tem1 = vk * wscaleu(i) * tem
828 dkt(i,k) = tem1 * prinv(i)
830 tem1 = vk * wscale(i) * tem
834 dkt(i,k) = tem1 * prinv(i)
836 dku(i,k) = min(dku(i,k),dkmax)
837 dku(i,k) = max(dku(i,k),xkzmo(i,k))
838 dkt(i,k) = min(dkt(i,k),dkmax)
839 dkt(i,k) = max(dkt(i,k),xkzo(i,k))
848 if (k < kpbl(i))
then
851 zfac = max((1.-zi(i,k+1)/hpbl(i)), zfmin)
852 tem = zi(i,k+1) * (zfac**pfac) * abs(moninq_fac)
855 if (useshape .ge. 1)
then
856 sz2h=(zi(i,k+1)-zl(i,1))/(hpbl(i)-zl(i,1))
859 zfac=(1.0-sz2h)**pfac
863 skmax=hmax*(1.0-hmax)**pfac
864 sksfc=min(zi(i,2)/hpbl(i),0.05)
865 sksfc=sksfc*(1-sksfc)**pfac
868 ashape=max(abs(moninq_fac),0.2)
869 if (useshape == 1)
then
870 ashape=(1.0 - ((sz2h*zfac/smax)**0.25) *(1.0 - ashape))
871 tem = zi(i,k+1) * (zfac) * ashape
872 elseif (useshape == 2)
then
874 if (skmax > sksfc)
then
875 ashape1=(skmax*ashape-sksfc)/(skmax-sksfc)
877 skminusk0 = zi(i,k+1)*zfac - hpbl(i)*sksfc
878 tem = zi(i,k+1) * (zfac)
879 if (skminusk0 > 0)
then
880 tem = skminusk0*ashape1 + hpbl(i)*sksfc
894 tem1 = vk * wscaleu(i) * tem
898 dkt(i,k) = tem1 * prinv(i)
900 tem1 = vk * wscale(i) * tem
904 dkt(i,k) = tem1 * prinv(i)
906 dku(i,k) = min(dku(i,k),dkmax)
907 dku(i,k) = max(dku(i,k),xkzmo(i,k))
908 dkt(i,k) = min(dkt(i,k),dkmax)
909 dkt(i,k) = max(dkt(i,k),xkzo(i,k))
916 if (islimsk(i) .eq. 0)
then
922 if (moninq_fac .lt. 0.)
then
924 kloc = int(wspm(i,2))
929 if(kpbl(i) .gt. kloc)
then
937 if(xdku .ge. wspm(i,1))
then
938 wspm(i,3) = wspm(i,1)/xdku
940 wspm(i,4) = min(wspm(i,3),1.0)
946 zfac = max((1.-zi(i,k+1)/hpbl(i)), zfmin)
947 tem = zi(i,k+1) * (zfac**pfac) * wspm(i,4)
950 if(useshape .ge. 1)
then
951 sz2h=(zi(i,k+1)-zl(i,1))/(hpbl(i)-zl(i,1))
954 zfac=(1.0-sz2h)**pfac
957 skmax=hmax*(1.0-hmax)**pfac
958 sksfc=min(zi(i,2)/hpbl(i),0.05)
959 sksfc=sksfc*(1-sksfc)**pfac
962 ashape=max(wspm(i,4),0.2)
963 if(useshape ==1)
then
964 ashape=(1.0 - ((sz2h*zfac/smax)**0.25)*
966 tem = zi(i,k+1) * (zfac) * ashape
967 elseif (useshape == 2)
then
969 if (skmax > sksfc)
then
970 ashape1=(skmax*ashape-sksfc)/(skmax-sksfc)
972 skminusk0=zi(i,k+1)*zfac - hpbl(i)*sksfc
973 tem = zi(i,k+1) * (zfac)
974 if (skminusk0 > 0)
then
975 tem = skminusk0*ashape1 + hpbl(i)*sksfc
981 tem1 = vk * wscaleu(i) * tem
985 dkt(i,k) = tem1 * prinv(i)
987 tem1 = vk * wscale(i) * tem
991 dkt(i,k) = tem1 * prinv(i)
993 dku(i,k) = min(dku(i,k),dkmax)
994 dku(i,k) = max(dku(i,k),xkzmo(i,k))
995 dkt(i,k) = min(dkt(i,k),dkmax)
996 dkt(i,k) = max(dkt(i,k),xkzo(i,k))
1044 if(k >= kpbl(i))
then
1045 bvf2 = grav*bf(i,k)*ti(i,k)
1046 ri = max(bvf2/shr2(i,k),rimin)
1049 rl2 = zk*rlamun/(rlamun+zk)
1050 dk = rl2*rl2*sqrt(shr2(i,k))
1054 dku(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri))
1055 dkt(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri))
1057 rl2 = zk*rlam/(rlam+zk)
1060 dk = rl2*rl2*sqrt(shr2(i,k))
1061 tem1 = dk/(1+5.*ri)**2
1063 if(k >= kpblx(i))
then
1064 prnum = 1.0 + 2.1*ri
1065 prnum = min(prnum,prmax)
1071 dku(i,k) = tem1 * prnum
1075 dku(i,k) = min(dku(i,k),dkmax)
1076 dku(i,k) = max(dku(i,k),xkzmo(i,k))
1077 dkt(i,k) = min(dkt(i,k),dkmax)
1078 dkt(i,k) = max(dkt(i,k),xkzo(i,k))
1103 qcko(i,k,kk) = q1(i,k,kk)
1109 call mfpbl(im,im,km,ntrac,dt2,pcnvflg,
1110 & zl,zi,thvx,q1,t1,u1,v1,hpbl,kpbl,
1111 & sflux,ustar,wstar,xmf,tcko,qcko,ucko,vcko)
1131 tem = thetae(i,k) - thetae(i,k+1)
1132 tem1 = qtx(i,k) - qtx(i,k+1)
1133 if (tem > 0. .and. tem1 > 0.)
then
1134 cteit= cp*tem/(hvap*tem1)
1135 if(cteit > actei) rent(i) = rentf2
1142 tem1 = max(bf(i,k),tdzmin)
1143 ckt(i,k) = -rent(i)*radmin(i)/tem1
1150 if(scuflg(i) .and. k < krad(i))
then
1155 if(ptem.ge.1.) ptem= 1.
1156 ptem= tem2*ptem*sqrt(1.-ptem)
1157 ckt(i,k) = radfac*vk*vrad(i)*ptem
1158 cku(i,k) = 0.75*ckt(i,k)
1159 ckt(i,k) = max(ckt(i,k),dkmin)
1160 ckt(i,k) = min(ckt(i,k),dkmax)
1161 cku(i,k) = max(cku(i,k),dkmin)
1162 cku(i,k) = min(cku(i,k),dkmax)
1171 if (.not. hurr_pbl)
then
1175 dkt(i,k) = dkt(i,k)+ckt(i,k)
1176 dku(i,k) = dku(i,k)+cku(i,k)
1177 dkt(i,k) = min(dkt(i,k),dkmax)
1178 dku(i,k) = min(dku(i,k),dkmax)
1187 if (.not. (hurr_pbl .and. moninq_fac < 0.0))
then
1188 dkt(i,k) = dkt(i,k)+ckt(i,k)
1189 dku(i,k) = dku(i,k)+cku(i,k)
1191 dkt(i,k) = min(dkt(i,k),dkmax)
1192 dku(i,k) = min(dku(i,k),dkmax)
1204 a1(i,1) = t1(i,1) + beta(i) * heat(i)
1205 a2(i,1) = q1(i,1,1) + beta(i) * evap(i)
1212 a2(i,1+is) = q1(i,1,k)
1219 dtodsd = dt2/del(i,k)
1220 dtodsu = dt2/del(i,k+1)
1221 dsig = prsl(i,k)-prsl(i,k+1)
1223 tem1 = dsig * dkt(i,k) * rdz
1225 au(i,k) = -dtodsd*dsdz2
1226 al(i,k) = -dtodsu*dsdz2
1228 if(pcnvflg(i) .and. k < kpbl(i))
then
1230 ptem = 0.5 * tem2 * xmf(i,k)
1231 ptem1 = dtodsd * ptem
1232 ptem2 = dtodsu * ptem
1233 ad(i,k) = ad(i,k)-au(i,k)-ptem1
1234 ad(i,k+1) = 1.-al(i,k)+ptem2
1235 au(i,k) = au(i,k)-ptem1
1236 al(i,k) = al(i,k)+ptem2
1237 ptem = tcko(i,k) + tcko(i,k+1)
1239 a1(i,k) = a1(i,k)+dtodsd*dsdzt-ptem1*ptem
1240 a1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt+ptem2*ptem
1241 ptem = qcko(i,k,1) + qcko(i,k+1,1)
1242 a2(i,k) = a2(i,k) - ptem1 * ptem
1243 a2(i,k+1) = q1(i,k+1,1) + ptem2 * ptem
1244 elseif(ublflg(i) .and. k < kpbl(i))
then
1245 ptem1 = dsig * dktx(i,k) * rdz
1247 dsdzt = tem1 * gocp - ptem1 * hgamt(i) * tem
1248 dsdzq = - ptem1 * hgamq(i) * tem
1249 ad(i,k) = ad(i,k)-au(i,k)
1250 ad(i,k+1) = 1.-al(i,k)
1251 a1(i,k) = a1(i,k)+dtodsd*dsdzt
1252 a1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt
1253 a2(i,k) = a2(i,k)+dtodsd*dsdzq
1254 a2(i,k+1) = q1(i,k+1,1)-dtodsu*dsdzq
1256 ad(i,k) = ad(i,k)-au(i,k)
1257 ad(i,k+1) = 1.-al(i,k)
1259 a1(i,k) = a1(i,k)+dtodsd*dsdzt
1260 a1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt
1261 a2(i,k+1) = q1(i,k+1,1)
1272 if(pcnvflg(i) .and. k < kpbl(i))
then
1273 dtodsd = dt2/del(i,k)
1274 dtodsu = dt2/del(i,k+1)
1275 dsig = prsl(i,k)-prsl(i,k+1)
1276 tem = dsig * rdzt(i,k)
1277 ptem = 0.5 * tem * xmf(i,k)
1278 ptem1 = dtodsd * ptem
1279 ptem2 = dtodsu * ptem
1280 tem1 = qcko(i,k,kk) + qcko(i,k+1,kk)
1281 a2(i,k+is) = a2(i,k+is) - ptem1*tem1
1282 a2(i,k+1+is)= q1(i,k+1,kk) + ptem2*tem1
1284 a2(i,k+1+is) = q1(i,k+1,kk)
1294 call tridin(im,km,ntrac,al,ad,au,a1,a2,au,a1,a2)
1302 ttend = (a1(i,k)-t1(i,k)) * rdt
1303 qtend = (a2(i,k)-q1(i,k,1))*rdt
1304 tau(i,k) = tau(i,k)+ttend
1305 rtg(i,k,1) = rtg(i,k,1)+qtend
1306 dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend
1307 dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend
1310 if(.not.flag_for_pbl_generic_tend)
then
1311 idtend1=dtidx(index_of_temperature,index_of_process_pbl)
1312 idtend2=dtidx(ntqv+100,index_of_process_pbl)
1316 ttend = (a1(i,k)-t1(i,k)) * rdt
1317 dtend(i,k,idtend1) = dtend(i,k,idtend1) + ttend*delt
1324 qtend = (a2(i,k)-q1(i,k,1))*rdt
1325 dtend(i,k,idtend2) = dtend(i,k,idtend2) + qtend*delt
1335 qtend = (a2(i,k+is)-q1(i,k,kk))*rdt
1336 rtg(i,k,kk) = rtg(i,k,kk)+qtend
1340 if(.not.flag_for_pbl_generic_tend .and. ldiag3d .and. &
1341 & rtg_ozone_index>0)
then
1342 idtend1 = dtidx(100+ntoz,index_of_process_pbl)
1344 kk = rtg_ozone_index
1348 qtend = (a2(i,k+is)-q1(i,k,kk))
1349 dtend(i,k,idtend1) = dtend(i,k,idtend1)+qtend
1364 diss(i,k) = dku(i,k)*shr2(i,k)-grav*ti(i,k)*dkt(i,k)*bf(i,k)
1372 if (hurr_pbl .and. moninq_fac < 0.0)
then
1379 tem = govrth(i)*sflux(i)
1380 tem1 = tem + stress(i)*spd1(i)/zl(i,1)
1381 tem2 = 0.5 * (tem1+diss(i,1))
1382 tem2 = max(tem2, 0.)
1384 tau(i,1) = tau(i,1)+ttend_fac*ttend
1391 tem = 0.5 * (diss(i,k-1)+diss(i,k))
1394 tau(i,k) = tau(i,k) + ttend_fac*ttend
1405 ad(i,1) = 1.0 + beta(i) * stress(i) / spd1(i)
1412 dtodsd = dt2/del(i,k)
1413 dtodsu = dt2/del(i,k+1)
1414 dsig = prsl(i,k)-prsl(i,k+1)
1416 tem1 = dsig*dku(i,k)*rdz
1418 au(i,k) = -dtodsd*dsdz2
1419 al(i,k) = -dtodsu*dsdz2
1421 if(pcnvflg(i) .and. k < kpbl(i))
then
1423 ptem = 0.5 * tem2 * xmf(i,k)
1424 ptem1 = dtodsd * ptem
1425 ptem2 = dtodsu * ptem
1426 ad(i,k) = ad(i,k)-au(i,k)-ptem1
1427 ad(i,k+1) = 1.-al(i,k)+ptem2
1428 au(i,k) = au(i,k)-ptem1
1429 al(i,k) = al(i,k)+ptem2
1430 ptem = ucko(i,k) + ucko(i,k+1)
1431 a1(i,k) = a1(i,k) - ptem1 * ptem
1432 a1(i,k+1) = u1(i,k+1) + ptem2 * ptem
1433 ptem = vcko(i,k) + vcko(i,k+1)
1434 a2(i,k) = a2(i,k) - ptem1 * ptem
1435 a2(i,k+1) = v1(i,k+1) + ptem2 * ptem
1437 ad(i,k) = ad(i,k)-au(i,k)
1438 ad(i,k+1) = 1.-al(i,k)
1439 a1(i,k+1) = u1(i,k+1)
1440 a2(i,k+1) = v1(i,k+1)
1449 call tridi2(im,km,al,ad,au,a1,a2,au,a1,a2)
1456 utend = (a1(i,k)-u1(i,k))*rdt
1457 vtend = (a2(i,k)-v1(i,k))*rdt
1458 du(i,k) = du(i,k) + utend
1459 dv(i,k) = dv(i,k) + vtend
1460 dusfc(i) = dusfc(i) + conw*del(i,k)*utend
1461 dvsfc(i) = dvsfc(i) + conw*del(i,k)*vtend
1474 if(.not.flag_for_pbl_generic_tend)
then
1475 idtend1 = dtidx(index_of_x_wind,index_of_process_pbl)
1479 utend = (a1(i,k)-u1(i,k))*rdt
1480 dtend(i,k,idtend1) = dtend(i,k,idtend1) + utend*delt
1485 idtend2 = dtidx(index_of_y_wind,index_of_process_pbl)
1489 vtend = (a2(i,k)-v1(i,k))*rdt
1490 dtend(i,k,idtend2) = dtend(i,k,idtend2) + vtend*delt