30 (nd, nl, delt, iflag, &
31 precip, wd, tprime, qprime, cbmf )
182 integer :: iflag, nd, nl
184 real :: cbmf, delt, precip, qprime, tprime, wd
188 integer :: i, icb, ihmin, inb, inb1, j, jtt, k
191 real :: ad, afac, ahmax, ahmin, alt, altem
192 real :: am, amp1, anum, asij, awat, b6, bf2, bsum, by
193 real :: byp, c6, cape, capem, cbmfold, chi, coeff
194 real :: cpinv, cwat, damps, dbo, dbosum
195 real :: defrac, dei, delm, delp, delt0, delti, denom, dhdp
196 real :: dpinv, dtma, dtmin, dtpbl, elacrit, ents
197 real :: epmax, fac, fqold, frac, ftold
198 real :: plcl, qp1, qsm, qstm, qti, rat
199 real :: rdcp, revap, rh, scrit, sigt, sjmax
200 real :: sjmin, smid, smin, stemp, tca
201 real :: tvaplcl, tvpplcl, tvx, tvy, wdtrain
206 real :: fup(na),fdown(na)
211 REAL :: m(na),mp(na),ment(na,na),qent(na,na),elij(na,na)
212 REAL :: sij(na,na),tvp(na),tv(na),water(na)
213 REAL :: qp(na),ep(na),th(na),wt(na),evap(na),clw(na)
214 REAL :: sigp(na),tp(na),cpn(na)
215 REAL :: lv(na),lvcp(na),h(na),hp(na),gz(na),hm(na)
233 INTEGER,
PARAMETER :: ipbl=0
234 INTEGER,
PARAMETER :: minorig=1
265 REAL,
PARAMETER :: elcrit=.0011
266 REAL,
PARAMETER :: tlcrit=-55.0
267 REAL,
PARAMETER :: entp=1.5
268 REAL,
PARAMETER :: sigd=0.05
269 REAL,
PARAMETER :: sigs=0.12
270 REAL,
PARAMETER :: omtrain=50.0
271 REAL,
PARAMETER :: omtsnow=5.5
272 REAL,
PARAMETER :: coeffr=1.0
273 REAL,
PARAMETER :: coeffs=0.8
274 REAL,
PARAMETER :: cu=0.7
275 REAL,
PARAMETER :: beta=10.0
276 REAL,
PARAMETER :: dtmax=0.9
277 REAL,
PARAMETER :: alpha=0.025
278 REAL,
PARAMETER :: damp=0.1
286 REAL,
PARAMETER :: cpd=1005.7
287 REAL,
PARAMETER :: cpv=1870.0
288 REAL,
PARAMETER :: cl=2500.0
289 REAL,
PARAMETER :: rv=461.5
290 REAL,
PARAMETER :: rd=287.04
291 REAL,
PARAMETER :: lv0=2.501e6
292 REAL,
PARAMETER :: g=9.81
293 REAL,
PARAMETER :: rowl=1000.0
295 REAL,
PARAMETER :: cpvmcl=cl-cpv
296 REAL,
PARAMETER :: eps0=rd/rv
297 REAL,
PARAMETER :: epsi=1./eps0
298 REAL,
PARAMETER :: ginv=1.0/g
299 REAL,
PARAMETER :: epsilon=1.e-20
322 rdcp=(rd*(1.-qconv(i))+qconv(i)*rv)/ &
323 (cpd*(1.-qconv(i))+qconv(i)*cpv)
324 th(i)=tconv(i)*(1000.0/pconv_hpa(i))**rdcp
412 cpn(1)=cpd*(1.-qconv(1))+qconv(1)*cpv
414 lv(1)=lv0-cpvmcl*(tconv(1)-273.15)
416 tv(1)=tconv(1)*(1.+qconv(1)*epsi-qconv(1))
420 tvx=tconv(i)*(1.+qconv(i)*epsi-qconv(i))
421 tvy=tconv(i-1)*(1.+qconv(i-1)*epsi-qconv(i-1))
422 gz(i)=gz(i-1)+0.5*rd*(tvx+tvy)*(pconv_hpa(i-1)-pconv_hpa(i))/ &
424 cpn(i)=cpd*(1.-qconv(i))+cpv*qconv(i)
425 h(i)=tconv(i)*cpn(i)+gz(i)
426 lv(i)=lv0-cpvmcl*(tconv(i)-273.15)
427 hm(i)=(cpd*(1.-qconv(i))+cl*qconv(i))*(tconv(i)-tconv(1))+ &
429 tv(i)=tconv(i)*(1.+qconv(i)*epsi-qconv(i))
433 IF(i.GE.minorig.AND.hm(i).LT.ahmin.AND.hm(i).LT.hm(i-1))
THEN
438 ihmin=min(ihmin, nl-1)
448 IF(hm(i).GT.ahmax)
THEN
458 IF(tconv(nk).LT.250.0.OR.qconv(nk).LE.0.0.OR.ihmin.EQ.(nl-1)) &
468 rh=qconv(nk)/qsconv(nk)
469 chi=tconv(nk)/(1669.0-122.0*rh-tconv(nk))
470 plcl=pconv_hpa(nk)*(rh**chi)
471 IF(plcl.LT.200.0.OR.plcl.GE.2000.0)
THEN
481 IF(pconv_hpa(i).LT.plcl)
THEN
485 IF(icb.GE.(nl-1))
THEN
497 CALL
tlift(gz,icb,nk,tvp,tp,clw,nd,nl,1)
499 tvp(i)=tvp(i)-tp(i)*qconv(nk)
505 IF(cbmf.EQ.0.0.AND.tvp(icb).LE.(tv(icb)-dtmax))
THEN
512 IF(iflag.NE.4)iflag=1
516 CALL
tlift(gz,icb,nk,tvp,tp,clw,nd,nl,2)
531 elacrit=elcrit*(1.0-tca/tlcrit)
533 elacrit=max(elacrit,0.0)
535 ep(i)=epmax*(1.0-elacrit/max(clw(i),1.0e-8))
537 ep(i)=min(ep(i),epmax)
545 tvp(i)=tvp(i)-tp(i)*qconv(nk)
547 tvp(nl+1)=tvp(nl)-(gz(nl+1)-gz(nl))/cpd
579 by=(tvp(i)-tv(i))*(phconv_hpa(i)-phconv_hpa(i+1))/pconv_hpa(i)
581 IF(by.GE.0.0)inb1=i+1
584 byp=(tvp(i+1)-tv(i+1))*(phconv_hpa(i+1)-phconv_hpa(i+2))/ &
592 defrac=max(defrac,0.001)
600 hp(i)=h(nk)+(lv(i)+(cpd-cpv)*tconv(i))*ep(i)*clw(i)
611 tvpplcl=tvp(icb-1)-rd*tvp(icb-1)*(pconv_hpa(icb-1)-plcl)/ &
612 (cpn(icb-1)*pconv_hpa(icb-1))
613 tvaplcl=tv(icb)+(tvp(icb)-tvp(icb+1))*(plcl-pconv_hpa(icb))/ &
614 (pconv_hpa(icb)-pconv_hpa(icb+1))
617 dtpbl=dtpbl+(tvp(i)-tv(i))*(phconv_hpa(i)-phconv_hpa(i+1))
619 dtpbl=dtpbl/(phconv_hpa(nk)-phconv_hpa(icb))
620 dtmin=tvpplcl-tvaplcl+dtmax+dtpbl
628 damps=damp*delt/delt0
629 cbmf=(1.-damps)*cbmf+0.1*alpha*dtma
634 IF(cbmf.EQ.0.0.AND.cbmfold.EQ.0.0)
THEN
644 dbo=abs(tv(k)-tvp(k))+ &
645 entp*0.02*(phconv_hpa(k)-phconv_hpa(k+1))
658 qti=qconv(nk)-ep(i)*clw(i)
660 bf2=1.+lv(j)*lv(j)*qsconv(j)/(rv*tconv(j)*tconv(j)*cpd)
661 anum=h(j)-hp(i)+(cpv-cpd)*tconv(j)*(qti-qconv(j))
662 denom=h(i)-hp(i)+(cpd-cpv)*(qconv(i)-qti)*tconv(j)
664 IF(abs(dei).LT.0.01)dei=0.01
667 altem=sij(i,j)*qconv(i)+(1.-sij(i,j))*qti-qsconv(j)
669 cwat=clw(j)*(1.-ep(j))
671 IF((stemp.LT.0.0.OR.stemp.GT.1.0.OR. &
672 altem.GT.cwat).AND.j.GT.i)
THEN
673 anum=anum-lv(j)*(qti-qsconv(j)-cwat*bf2)
674 denom=denom+lv(j)*(qconv(i)-qti)
675 IF(abs(denom).LT.0.01)denom=0.01
677 altem=sij(i,j)*qconv(i)+(1.-sij(i,j))*qti-qsconv(j)
678 altem=altem-(bf2-1.)*cwat
680 IF(sij(i,j).GT.0.0.AND.sij(i,j).LT.0.9)
THEN
681 qent(i,j)=sij(i,j)*qconv(i)+(1.-sij(i,j))*qti
683 elij(i,j)=max(0.0,elij(i,j))
684 ment(i,j)=m(i)/(1.-sij(i,j))
687 sij(i,j)=max(0.0,sij(i,j))
688 sij(i,j)=min(1.0,sij(i,j))
696 qent(i,i)=qconv(nk)-ep(i)*clw(i)
708 qp1=qconv(nk)-ep(i)*clw(i)
709 anum=h(i)-hp(i)-lv(i)*(qp1-qsconv(i))
710 denom=h(i)-hp(i)+lv(i)*(qconv(i)-qp1)
711 IF(abs(denom).LT.0.01)denom=0.01
713 alt=qp1-qsconv(i)+scrit*(qconv(i)-qp1)
714 IF(alt.LT.0.0)scrit=1.0
719 IF(sij(i,j).GT.0.0.AND.sij(i,j).LT.0.9)
THEN
721 smid=min(sij(i,j),scrit)
724 IF(smid.LT.smin.AND.sij(i,j+1).LT.smid)
THEN
726 sjmax=min(sij(i,j+1),sij(i,j),scrit)
727 sjmin=max(sij(i,j-1),sij(i,j))
728 sjmin=min(sjmin,scrit)
731 sjmax=max(sij(i,j+1),scrit)
732 smid=max(sij(i,j),scrit)
734 IF(j.GT.1)sjmin=sij(i,j-1)
735 sjmin=max(sjmin,scrit)
739 asij=asij+(delp+delm)*(phconv_hpa(j)-phconv_hpa(j+1))
740 ment(i,j)=ment(i,j)*(delp+delm)* &
741 (phconv_hpa(j)-phconv_hpa(j+1))
744 asij=max(1.0e-21,asij)
747 ment(i,j)=ment(i,j)*asij
753 IF(bsum.LT.1.0e-18)
THEN
756 qent(i,i)=qconv(nk)-ep(i)*clw(i)
766 IF(ep(inb).LT.0.0001)goto 405
779 wdtrain=g*ep(i)*m(i)*clw(i)
782 awat=elij(j,i)-(1.-ep(i))*clw(i)
784 wdtrain=wdtrain+g*awat*ment(j,i)
799 IF(tconv(i).GT.273.0)
THEN
803 qsm=0.5*(qconv(i)+qp(i+1))
804 afac=coeff*phconv_hpa(i)*(qsconv(i)-qsm)/ &
805 (1.0e4+2.0e3*phconv_hpa(i)*qsconv(i))
810 b6=100.*(phconv_hpa(i)-phconv_hpa(i+1))*sigt*afac/wt(i)
811 c6=(water(i+1)*wt(i+1)+wdtrain/sigd)/wt(i)
812 revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
813 evap(i)=sigt*afac*revap
820 dhdp=(h(i)-h(i-1))/(pconv_hpa(i-1)-pconv_hpa(i))
822 mp(i)=100.*ginv*lv(i)*sigd*evap(i)/dhdp
827 fac=20.0/(phconv_hpa(i-1)-phconv_hpa(i))
828 mp(i)=(fac*mp(i+1)+mp(i))/(1.+fac)
833 IF(pconv_hpa(i).GT.(0.949*pconv_hpa(1)))
THEN
835 mp(i)=mp(jtt)*(pconv_hpa(1)-pconv_hpa(i))/(pconv_hpa(1)- &
848 IF(mp(i).GT.mp(i+1))
THEN
850 qp(i)=qp(i+1)*rat+qconv(i)*(1.0-rat)+100.*ginv* &
851 sigd*(phconv_hpa(i)-phconv_hpa(i+1))*(evap(i)/mp(i))
853 IF(mp(i+1).GT.0.0)
THEN
854 qp(i)=(gz(i+1)-gz(i)+qp(i+1)*(lv(i+1)+tconv(i+1)*(cl-cpd))+ &
855 cpd*(tconv(i+1)-tconv(i)))/(lv(i)+tconv(i)*(cl-cpd))
858 qp(i)=min(qp(i),qstm)
865 precip=precip+wt(1)*sigd*water(1)*3600.*24000./(rowl*g)
872 wd=beta*abs(mp(icb))*0.01*rd*tconv(icb)/(sigd*pconv_hpa(icb))
873 qprime=0.5*(qp(1)-qconv(1))
874 tprime=lv0*qprime/cpd
879 dpinv=0.01/(phconv_hpa(1)-phconv_hpa(2))
888 IF((2.*g*dpinv*am).GE.delti)iflag=4
889 ft(1)=ft(1)+g*dpinv*am*(tconv(2)-tconv(1)+(gz(2)-gz(1))/cpn(1))
890 ft(1)=ft(1)-lvcp(1)*sigd*evap(1)
891 ft(1)=ft(1)+sigd*wt(2)*(cl-cpd)*water(2)*(tconv(2)- &
892 tconv(1))*dpinv/cpn(1)
893 fq(1)=fq(1)+g*mp(2)*(qp(2)-qconv(1))* &
895 fq(1)=fq(1)+g*am*(qconv(2)-qconv(1))*dpinv
897 fq(1)=fq(1)+g*dpinv*ment(j,1)*(qent(j,1)-qconv(1))
907 dpinv=0.01/(phconv_hpa(i)-phconv_hpa(i+1))
923 IF((2.*g*dpinv*amp1).GE.delti)iflag=4
931 ft(i)=ft(i)+g*dpinv*(amp1*(tconv(i+1)-tconv(i)+(gz(i+1)-gz(i))* &
932 cpinv)-ad*(tconv(i)-tconv(i-1)+(gz(i)-gz(i-1))*cpinv)) &
933 -sigd*lvcp(i)*evap(i)
934 ft(i)=ft(i)+g*dpinv*ment(i,i)*(hp(i)-h(i)+ &
935 tconv(i)*(cpv-cpd)*(qconv(i)-qent(i,i)))*cpinv
936 ft(i)=ft(i)+sigd*wt(i+1)*(cl-cpd)*water(i+1)* &
937 (tconv(i+1)-tconv(i))*dpinv*cpinv
938 fq(i)=fq(i)+g*dpinv*(amp1*(qconv(i+1)-qconv(i))- &
939 ad*(qconv(i)-qconv(i-1)))
941 awat=elij(k,i)-(1.-ep(i))*clw(i)
943 fq(i)=fq(i)+g*dpinv*ment(k,i)*(qent(k,i)-awat-qconv(i))
946 fq(i)=fq(i)+g*dpinv*ment(k,i)*(qent(k,i)-qconv(i))
948 fq(i)=fq(i)+sigd*evap(i)+g*(mp(i+1)* &
949 (qp(i+1)-qconv(i))-mp(i)*(qp(i)-qconv(i-1)))*dpinv
956 fq(inb)=fq(inb)*(1.-frac)
957 fq(inb-1)=fq(inb-1)+frac*fqold*((phconv_hpa(inb)- &
958 phconv_hpa(inb+1))/ &
959 (phconv_hpa(inb-1)-phconv_hpa(inb)))*lv(inb)/lv(inb-1)
961 ft(inb)=ft(inb)*(1.-frac)
962 ft(inb-1)=ft(inb-1)+frac*ftold*((phconv_hpa(inb)- &
963 phconv_hpa(inb+1))/ &
964 (phconv_hpa(inb-1)-phconv_hpa(inb)))*cpn(inb)/cpn(inb-1)
971 ents=ents+(cpn(i)*ft(i)+lv(i)*fq(i))* &
972 (phconv_hpa(i)-phconv_hpa(i+1))
974 ents=ents/(phconv_hpa(1)-phconv_hpa(inb+1))
976 ft(i)=ft(i)-ents/cpn(i)
998 fmass(j,i)=fmass(j,i)+m(i)
1000 fmass(j,i)=fmass(j,i)+ment(j,i)
1001 IF (fmass(j,i).GT.epsilon) nconvtop=max(nconvtop,i,j)
1004 sub(i)=fup(i-1)-fdown(i)
1015 SUBROUTINE tlift(GZ,ICB,NK,TVP,TPK,CLW,ND,NL,KK)
1027 integer :: icb, kk, nd, nk, nl
1031 integer :: i, j, nsb, nst
1033 real :: ah0, ahg, alv, cpinv, cpp, denom
1034 real :: es, qg, rg, s, tc, tg
1038 REAL :: gz(nd),tpk(nd),clw(nd)
1043 REAL,
PARAMETER :: cpd=1005.7
1044 REAL,
PARAMETER :: cpv=1870.0
1045 REAL,
PARAMETER :: cl=2500.0
1046 REAL,
PARAMETER :: rv=461.5
1047 REAL,
PARAMETER :: rd=287.04
1048 REAL,
PARAMETER :: lv0=2.501e6
1050 REAL,
PARAMETER :: cpvmcl=cl-cpv
1051 REAL,
PARAMETER :: eps0=rd/rv
1052 REAL,
PARAMETER :: epsi=1./eps0
1056 ah0=(cpd*(1.-qconv(nk))+cl*qconv(nk))*tconv(nk)+qconv(nk)* &
1058 tconv(nk)-273.15))+gz(nk)
1059 cpp=cpd*(1.-qconv(nk))+qconv(nk)*cpv
1070 tpk(i)=tconv(nk)-(gz(i)-gz(nk))*cpinv
1071 tvp(i)=tpk(i)*(1.+qconv(nk)*epsi)
1086 alv=lv0-cpvmcl*(tconv(i)-273.15)
1088 s=cpd+alv*alv*qg/(rv*tconv(i)*tconv(i))
1090 ahg=cpd*tg+(cl-cpd)*qconv(nk)*tconv(i)+alv*qg+gz(i)
1096 es=6.112*exp(17.67*tc/denom)
1098 es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
1100 qg=eps0*es/(pconv_hpa(i)-es*(1.-eps0))
1102 alv=lv0-cpvmcl*(tconv(i)-273.15)
1103 tpk(i)=(ah0-(cl-cpd)*qconv(nk)*tconv(i)-gz(i)-alv*qg)/cpd
1105 clw(i)=max(0.0,clw(i))
1106 rg=qg/(1.-qconv(nk))
1107 tvp(i)=tpk(i)*(1.+rg*epsi)
1110 END SUBROUTINE tlift
subroutine convect(ND, NL, DELT, IFLAG, PRECIP, WD, TPRIME, QPRIME, CBMF)
subroutine tlift(GZ, ICB, NK, TVP, TPK, CLW, ND, NL, KK)