53 & eps,epsm1,fv,grav,hvap,rd,rv, &
54 & t0c,delt,ntk,ntr,delp,first_time_step,restart, &
55 & tmf,qmicro,progsigma, &
56 & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, &
57 & rn,kbot,ktop,kcnv,islimsk,garea, &
58 & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, &
59 & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, &
60 & sigmain,sigmaout,betadcu,betamcu,betascu,errmsg,errflg)
62 use machine ,
only : kind_phys
63 use funcphys ,
only : fpvs
67 integer,
intent(in) :: im, km, itc, ntc, ntk, ntr, ncloud
68 integer,
intent(in) :: islimsk(:)
69 real(kind=kind_phys),
intent(in) :: cliq, cp, cvap, &
70 & eps, epsm1, fv, grav, hvap, rd, rv, t0c, betascu, betadcu, &
72 real(kind=kind_phys),
intent(in) :: delt
73 real(kind=kind_phys),
intent(in) :: psp(:), delp(:,:), &
74 & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), &
76 real(kind=kind_phys),
intent(in),
optional :: qmicro(:,:), &
78 real(kind=kind_phys),
intent(in),
optional :: sigmain(:,:)
80 real(kind=kind_phys),
dimension(:),
intent(in) :: fscav
81 integer,
intent(inout) :: kcnv(:)
83 real(kind=kind_phys),
intent(inout) :: qtr(:,:,:), &
84 & q1(:,:), t1(:,:), u1(:,:), v1(:,:)
86 integer,
intent(out) :: kbot(:), ktop(:)
87 real(kind=kind_phys),
intent(out) :: rn(:), &
88 & cnvw(:,:), cnvc(:,:), dt_mf(:,:)
90 real(kind=kind_phys),
intent(out),
optional :: ud_mf(:,:), &
92 real(kind=kind_phys),
intent(in) :: clam, c0s, c1, &
93 & asolfac, evef, pgcon
94 logical,
intent(in) :: hwrf_samfshal,first_time_step, &
96 character(len=*),
intent(out) :: errmsg
97 integer,
intent(out) :: errflg
100 integer i,j,indx, k, kk, km1, n
103 real(kind=kind_phys) clamd, tkemx, tkemn, dtke
105 real(kind=kind_phys) dellat,
108 & dq, dqsdp, dqsdt, dt,
113 & el2orc, elocp, aafac,
115 & es, etah, h1, shevf,
117 & fact1, fact2, factor,
118 & cthk, cthkmn, dthk,
119 & gamma, pprime, betaw, tauadv,
121 & rfact, shear, tfac,
126 & rho, tem, tem1, tem2,
129 integer kb(im), kb1(im), kbcon(im), kbcon1(im),
130 & ktcon(im), ktcon1(im),
133 real(kind=kind_phys) aa1(im), cina(im),
134 & tkemean(im), clamt(im),
135 & ps(im), del(im,km), prsl(im,km),
136 & umean(im), advfac(im), gdx(im),
137 & delhbar(im), delq(im), delq2(im),
138 & delqbar(im), delqev(im), deltbar(im),
140 & deltv(im), dtconv(im),
141 & pdot(im), po(im,km),
142 & qcond(im), qevap(im), hmax(im),
145 & xlamud(im), xmb(im), xmbmax(im),
147 & delubar(im), delvbar(im)
149 real(kind=kind_phys) c0(im), sfcpbl(im)
151 real(kind=kind_phys) crtlame, crtlamd
153 real(kind=kind_phys) cinpcr, cinpcrmx, cinpcrmn,
154 & cinacr, cinacrmx, cinacrmn,
159 real(kind=kind_phys) bb1, bb2, csmf, wucb
163 real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km),
164 & omegac(im),zeta(im,km),dbyo1(im,km),
165 & sigmab(im),qadv(im,km)
166 real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins,
168 logical flag_shallow,flag_mid
186 parameter(cm=1.0,cq=1.0)
188 parameter(clamd=0.1,tkemx=0.65,tkemn=0.05)
189 parameter(dtke=tkemx-tkemn)
190 parameter(cthk=200.,cthkmn=0.,dthk=25.)
191 parameter(sfclfac=0.2,rhcrt=0.75)
192 parameter(cinpcrmx=180.,cinpcrmn=120.)
194 parameter(cinacrmx=-120.,shevf=2.0)
195 parameter(dtmax=10800.,dtmin=600.)
196 parameter(bb1=4.0,bb2=0.8,csmf=0.2)
197 parameter(tkcrt=2.,cmxfac=10.)
199 parameter(betaw=.03,dxcrtc0=9.e3)
200 parameter(h1=0.33333333)
202 parameter(dxcrtas=30.e3,sigmind=0.01,sigmins=0.03,sigminm=0.01)
203c local variables and arrays
204 real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km),
205 & uo(im,km), vo(im,km), qeso(im,km),
206 & ctr(im,km,ntr), ctro(im,km,ntr)
209c variables for tracer wet deposition,
210 real(kind=kind_phys),
dimension(im,km,ntc) :: chem_c, chem_pw,
212 real(kind=kind_phys),
parameter :: escav = 0.8
215 real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km),
216 & wush(im,km), wc(im)
219 real(kind=kind_phys) scaldfunc(im), sigmagfm(im)
223 real(kind=kind_phys) qlko_ktcon(im), dellal(im,km),
224 & dbyo(im,km), zo(im,km), xlamue(im,km),
226 & heo(im,km), heso(im,km),
227 & dellah(im,km), dellaq(im,km),
229 & dellau(im,km), dellav(im,km), hcko(im,km),
230 & ucko(im,km), vcko(im,km), qcko(im,km),
231 & qrcko(im,km), ecko(im,km,ntr),
232 & ercko(im,km,ntr), eta(im,km),
233 & zi(im,km), pwo(im,km), c0t(im,km),
234 & sumx(im), tx1(im), cnvwt(im,km)
240 real(kind=kind_phys) q_diff(im,0:km-1), e_diff(im,0:km-1,ntr),
242 real(kind=kind_phys) rrkp, phkp
243 real(kind=kind_phys) tsumn(im), tsump(im), rtnp(im)
245 logical do_aerosols, totflg, cnvflg(im), flg(im)
247 real(kind=kind_phys) tf, tcr, tcrf
248 parameter(tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
251c-----------------------------------------------------------------------
261 el2orc = hvap*hvap/(rv*cp)
263 fact1 = (cvap-cliq)/rv
264 fact2 = hvap/rv-fact1*t0c
266 if (.not.hwrf_samfshal)
then
276c-----------------------------------------------------------------------
277 if (.not.hwrf_samfshal)
then
279 do_aerosols = (itc > 0) .and. (ntc > 0) .and. (ntr > 0)
280 if (do_aerosols) do_aerosols = (ntr >= itc + ntc - 3)
302 if(hwrf_samfshal)
then
305 if(kcnv(i) == 1) cnvflg(i) = .false.
310 sfcpbl(i) = sfclfac * hpbl(i)
322 gdx(i) = sqrt(garea(i))
331 if(kcnv(i) == 1) cnvflg(i) = .false.
336 sfcpbl(i) = sfclfac * hpbl(i)
347 gdx(i) = sqrt(garea(i))
356 totflg = totflg .and. (.not. cnvflg(i))
362 if(islimsk(i) == 1)
then
371 if(gdx(i) < dxcrtc0)
then
372 tem = gdx(i) / dxcrtc0
381 if(t1(i,k) > 273.16)
then
384 tem = d0 * (t1(i,k) - 273.16)
386 c0t(i,k) = c0(i) * tem1
409c model tunable parameters are all here
426c define top layer for search of the downdraft originating layer
427c and the maximum thetae for updraft
438 if (prsl(i,k)*tx1(i) > 0.70) kbm(i) = k + 1
439 if (prsl(i,k)*tx1(i) > 0.60) kmax(i) = k + 1
443 kbm(i) = min(kbm(i),kmax(i))
446c hydrostatic height assume zero terr and compute
447c updraft entrainment rate as an inverse
function of height
452 zo(i,k) = phil(i,k) / grav
456 if(hwrf_samfshal)
then
459 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
460 xlamue(i,k) = clam / zi(i,k)
464 xlamue(i,km) = xlamue(i,km1)
469 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
483 if (flg(i) .and. zo(i,k) <= hpbl(i))
then
491 kpbl(i)= min(kpbl(i),kbm(i))
495c convert surface pressure to mb from cb
500 if (cnvflg(i) .and. k <= kmax(i))
then
501 pfld(i,k) = prsl(i,k) * 10.0
542 if (.not.hwrf_samfshal)
then
547 if (cnvflg(i) .and. k <= kmax(i))
then
548 ctr(i,k,kk) = qtr(i,k,n)
549 ctro(i,k,kk) = qtr(i,k,n)
560 if (cnvflg(i) .and. k <= kmax(i))
then
561 qeso(i,k) = 0.01 * fpvs(to(i,k))
562 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
564 qeso(i,k) = max(qeso(i,k), val1)
566 qo(i,k) = max(qo(i,k), val2 )
573c compute moist static
energy
578 if (cnvflg(i) .and. k <= kmax(i))
then
580 tem = phil(i,k) + cp * to(i,k)
581 heo(i,k) = tem + hvap * qo(i,k)
582 heso(i,k) = tem + hvap * qeso(i,k)
583c heo(i,k) = min(heo(i,k),heso(i,k))
588c determine level with largest moist static
energy within pbl
589c this is the level
where updraft starts
599 if (flg(i) .and. zo(i,k) <= sfcpbl(i))
then
607 kb1(i) = min(kb1(i),kpbl(i))
613 hmax(i) = heo(i,kb1(i))
619 if(cnvflg(i) .and. (k > kb1(i) .and. k <= kpbl(i)))
then
620 if(heo(i,k) > hmax(i))
then
631 if (cnvflg(i) .and. k <= kmax(i)-1)
then
632 dz = .5 * (zo(i,k+1) - zo(i,k))
633 dp = .5 * (pfld(i,k+1) - pfld(i,k))
634 es = 0.01 * fpvs(to(i,k+1))
635 pprime = pfld(i,k+1) + epsm1 * es
636 qs = eps * es / pprime
637 dqsdp = - qs / pprime
638 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
639 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
640 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
641 dt = (grav*dz + hvap*dqsdp*dp) / (cp*(1. + gamma))
642 dq = dqsdt * dt + dqsdp * dp
643 to(i,k) = to(i,k+1) + dt
644 qo(i,k) = qo(i,k+1) + dq
645 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
653 if (cnvflg(i) .and. k <= kmax(i)-1)
then
654 qeso(i,k) = 0.01 * fpvs(to(i,k))
655 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
657 qeso(i,k) = max(qeso(i,k), val1)
659 qo(i,k) = max(qo(i,k), val2 )
661 rh(i,k) = min(qo(i,k)/qeso(i,k), 1.)
662 heo(i,k) = .5 * grav * (zo(i,k) + zo(i,k+1)) +
663 & cp * to(i,k) + hvap * qo(i,k)
664 heso(i,k) = .5 * grav * (zo(i,k) + zo(i,k+1)) +
665 & cp * to(i,k) + hvap * qeso(i,k)
666 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
667 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
672 if (.not.hwrf_samfshal)
then
676 if (cnvflg(i) .and. k <= kmax(i)-1)
then
677 ctro(i,k,n) = .5 * (ctro(i,k,n) + ctro(i,k+1,n))
684c look for the level of free convection as cloud base
689 if(flg(i)) kbcon(i) = kmax(i)
693 if (flg(i) .and. k < kbm(i))
then
694 if(k > kb(i) .and. heo(i,kb(i)) > heso(i,k))
then
704 if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
711 totflg = totflg .and. (.not. cnvflg(i))
719 pdot(i) = 0.01 * dot(i,kbcon(i))
723c turn off convection
if pressure depth between parcel source level
724c and cloud base is larger than a critical
value, cinpcr
728 if(islimsk(i) == 1)
then
739 if(pdot(i) <= w4)
then
740 tem = (pdot(i) - w4) / (w3 - w4)
741 elseif(pdot(i) >= -w4)
then
742 tem = - (pdot(i) + w4) / (w4 - w3)
751 ptem1= .5*(cinpcrmx-cinpcrmn)
752 cinpcr = cinpcrmx - ptem * ptem1
753 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
755 if(tem1 > cinpcr .and.
756 & zi(i,kbcon(i)) > hpbl(i))
then
764 totflg = totflg .and. (.not. cnvflg(i))
778 if (cnvflg(i) .and. k <= kpbl(i))
then
779 if(heo(i,k) > hmax(i))
then
789 if(flg(i)) kbcon(i) = kmax(i)
793 if (flg(i) .and. k < kbm(i))
then
794 if(k > kb(i) .and. heo(i,kb(i)) > heso(i,k))
then
804 if(kbcon(i) == kmax(i)) cnvflg(i) = .false.
810 totflg = totflg .and. (.not. cnvflg(i))
817 pdot(i) = 0.01 * dot(i,kbcon(i))
830 if(k >= kb(i) .and. k < kbcon(i))
then
831 dz = zo(i,k+1) - zo(i,k)
832 rhbar(i) = rhbar(i) + rh(i,k) * dz
833 sumx(i) = sumx(i) + dz
840 rhbar(i) = rhbar(i) / sumx(i)
841 if(rhbar(i) < rhcrt)
then
849 totflg = totflg .and. (.not. cnvflg(i))
860 if (hwrf_samfshal)
then
863 xlamud(i) = xlamue(i,kbcon(i))
879 if(k >= kb(i) .and. k < kbcon(i))
then
880 dz = zo(i,k+1) - zo(i,k)
881 tem = 0.5 * (qtr(i,k,ntk)+qtr(i,k+1,ntk))
882 tkemean(i) = tkemean(i) + tem * dz
883 sumx(i) = sumx(i) + dz
891 tkemean(i) = tkemean(i) / sumx(i)
892 if(tkemean(i) > tkemx)
then
893 clamt(i) = clam + clamd
894 else if(tkemean(i) < tkemn)
then
895 clamt(i) = clam - clamd
897 tem = tkemx - tkemean(i)
898 tem1 = 1. - 2. * tem / dtke
899 clamt(i) = clam + clamd * tem1
906 if(tkemean(i) > tkcrt)
then
907 tem = 1. + tkemean(i)/tkcrt
908 tem1 = min(tem, cmxfac)
909 clamt(i) = tem1 * clam
931 dz = zo(i,k+1) - zo(i,k)
932 xlamue(i,k) = clamt(i) / (zi(i,k) + dz)
933 xlamue(i,k) = max(xlamue(i,k), crtlame)
939 xlamue(i,km) = xlamue(i,km1)
943c specify the detrainment rate for the updrafts
952 xlamud(i) = 0.001 * clamt(i)
957c determine updraft mass flux for the subcloud layers
967 if(k < kbcon(i) .and. k >= kb(i))
then
968 dz = zi(i,k+1) - zi(i,k)
969 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
970 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
976c compute mass flux above cloud base
984 if(k > kbcon(i) .and. k < kmax(i))
then
985 dz = zi(i,k) - zi(i,k-1)
986 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
987 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
988 if(eta(i,k) <= 0.)
then
990 kbm(i) = min(kbm(i),kmax(i))
998c compute updraft cloud property
1004 hcko(i,indx) = heo(i,indx)
1005 ucko(i,indx) = uo(i,indx)
1006 vcko(i,indx) = vo(i,indx)
1010 if (.not. hwrf_samfshal)
then
1015 ecko(i,indx,n) = ctro(i,indx,n)
1016 ercko(i,indx,n) = ctro(i,indx,n)
1028 if(k > kb(i) .and. k < kmax(i))
then
1029 dz = zi(i,k) - zi(i,k-1)
1030 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1031 tem1 = 0.5 * xlamud(i) * dz
1032 factor = 1. + tem - tem1
1033 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
1034 & (heo(i,k)+heo(i,k-1)))/factor
1035 dbyo(i,k) = hcko(i,k) - heso(i,k)
1037 tem = 0.5 * cm * tem
1041 ucko(i,k) = ((1.-tem)*ucko(i,k-1)+ptem*uo(i,k)
1042 & +ptem1*uo(i,k-1))/factor
1043 vcko(i,k) = ((1.-tem)*vcko(i,k-1)+ptem*vo(i,k)
1044 & +ptem1*vo(i,k-1))/factor
1050 if (.not.hwrf_samfshal)
then
1051 if (do_aerosols)
then
1060 if(k > kb(i) .and. k < kmax(i))
then
1061 dz = zi(i,k) - zi(i,k-1)
1062 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1065 ecko(i,k,n) = ((1.-tem)*ecko(i,k-1,n)+tem*
1066 & (ctro(i,k,n)+ctro(i,k-1,n)))/factor
1067 ercko(i,k,n) = ecko(i,k,n)
1073 if (do_aerosols)
then
1079 if(k > kb(i) .and. k < kmax(i))
then
1080 dz = zi(i,k) - zi(i,k-1)
1081 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1084 ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem *
1085 & (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor
1086 ercko(i,k,kk) = ecko(i,k,kk)
1087 chem_c(i,k,n) = escav * fscav(n) * ecko(i,k,kk)
1088 tem = chem_c(i,k,n) / (1. + c0t(i,k) * dz)
1089 chem_pw(i,k,n) = c0t(i,k) * dz * tem * eta(i,k-1)
1090 ecko(i,k,kk) = tem + ecko(i,k,kk) - chem_c(i,k,n)
1099c taking account into convection inhibition due to existence of
1100c dry layers below cloud base
1109 if (flg(i) .and. k < kbm(i))
then
1110 if(k >= kbcon(i) .and. dbyo(i,k) > 0.)
then
1119 if(kbcon1(i) == kmax(i)) cnvflg(i) = .false.
1124 tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
1133 totflg = totflg .and. (.not. cnvflg(i))
1138c calculate convective inhibition
1144 if(k > kb(i) .and. k < kbcon1(i))
then
1145 dz1 = zo(i,k+1) - zo(i,k)
1146 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1147 rfact = 1. + fv * cp * gamma
1151 & dz1 * (grav / (cp * to(i,k)))
1152 & * dbyo(i,k) / (1. + gamma)
1158 & max(val,(qeso(i,k) - qo(i,k)))
1165 if (hwrf_samfshal)
then
1169 if(cina(i) < cinacr) cnvflg(i) = .false.
1175 if(islimsk(i) == 1)
then
1186 if(pdot(i) <= w4)
then
1187 tem = (pdot(i) - w4) / (w3 - w4)
1188 elseif(pdot(i) >= -w4)
then
1189 tem = - (pdot(i) + w4) / (w4 - w3)
1199 tem1= .5*(cinacrmx-cinacrmn)
1200 cinacr = cinacrmx - tem * tem1
1201 if(cina(i) < cinacr) cnvflg(i) = .false.
1208 totflg = totflg .and. (.not. cnvflg(i))
1213c determine first guess cloud top as the level of zero buoyancy
1214c limited to the level of p/ps=0.7
1219 if(flg(i)) ktcon(i) = 1
1223 if (flg(i) .and. k < kbm(i))
then
1224 if(k > kbcon1(i) .and. dbyo(i,k) < 0.)
then
1232c turn off shallow convection
if cloud depth is larger than cthk or less than cthkmn
1236 tem = pfld(i,kbcon(i))-pfld(i,ktcon(i))
1237 if(tem > cthk .or. tem < cthkmn)
then
1244c specify upper limit of mass flux at cloud base
1250 dp = 1000. * del(i,k)
1251 xmbmax(i) = dp / (grav * dt2)
1255c compute cloud moisture property and precipitation
1261 qcko(i,kb(i)) = qo(i,kb(i))
1262 qrcko(i,kb(i)) = qo(i,kb(i))
1269 if(k > kb(i) .and. k < ktcon(i))
then
1270 dz = zi(i,k) - zi(i,k-1)
1271 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1273 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1275 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1278 qcko(i,k) = ((1.-tem)*qcko(i,k-1)+tem*
1279 & (qo(i,k)+qo(i,k-1)))/factor
1280 qrcko(i,k) = qcko(i,k)
1282 dq = eta(i,k) * (qcko(i,k) - qrch)
1286c below lfc check
if there is excess moisture to release latent heat
1288 if(k >= kbcon(i) .and. dq > 0.)
then
1289 etah = .5 * (eta(i,k) + eta(i,k-1))
1290 dp = 1000. * del(i,k)
1292 ptem = c0t(i,k) + c1
1293 qlk = dq / (eta(i,k) + etah * ptem * dz)
1294 dellal(i,k) = etah * c1 * dz * qlk * grav / dp
1296 qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1298 buo(i,k) = buo(i,k) - grav * qlk
1299 qcko(i,k)= qlk + qrch
1300 pwo(i,k) = etah * c0t(i,k) * dz * qlk
1301 cnvwt(i,k) = etah * qlk * grav / dp
1302 zdqca(i,k)=dq/eta(i,k)
1307 if(k >= kbcon(i))
then
1308 rfact = 1. + fv * cp * gamma
1310 buo(i,k) = buo(i,k) + (grav / (cp * to(i,k)))
1311 & * dbyo(i,k) / (1. + gamma)
1314 buo(i,k) = buo(i,k) + grav * fv *
1315 & max(val,(qeso(i,k) - qo(i,k)))
1316 drag(i,k) = max(xlamue(i,k),xlamud(i))
1318 tem = ((uo(i,k)-uo(i,k-1))/dz)**2
1319 tem = tem+((vo(i,k)-vo(i,k-1))/dz)**2
1320 wush(i,k) = csmf * sqrt(tem)
1329c calculate cloud work function
1372 if(k >= kbcon(i) .and. k < ktcon(i))
then
1373 dz1 = zo(i,k+1) - zo(i,k)
1374 aa1(i) = aa1(i) + buo(i,k) * dz1
1375 dbyo1(i,k) = hcko(i,k) - heso(i,k)
1381 if(cnvflg(i) .and. aa1(i) <= 0.) cnvflg(i) = .false.
1387 totflg = totflg .and. (.not. cnvflg(i))
1392c estimate the onvective overshooting as the level
1393c
where the [aafac * cloud work function] becomes zero,
1394c which is the final cloud top
1395c limited to the level of p/ps=0.7
1400 aa1(i) = aafac * aa1(i)
1411 if(k >= ktcon(i) .and. k < kbm(i))
then
1412 dz1 = zo(i,k+1) - zo(i,k)
1413 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1414 rfact = 1. + fv * cp * gamma
1418 & dz1 * (grav / (cp * to(i,k)))
1419 & * dbyo(i,k) / (1. + gamma)
1426 if(aa1(i) < 0.)
then
1435c compute cloud moisture property, detraining cloud
water
1436c and precipitation in overshooting layers
1442 if(k >= ktcon(i) .and. k < ktcon1(i))
then
1443 dz = zi(i,k) - zi(i,k-1)
1444 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1446 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1448 tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1451 qcko(i,k) = ((1.-tem)*qcko(i,k-1)+tem*
1452 & (qo(i,k)+qo(i,k-1)))/factor
1453 qrcko(i,k) = qcko(i,k)
1455 dq = eta(i,k) * (qcko(i,k) - qrch)
1457c check
if there is excess moisture to release latent heat
1460 etah = .5 * (eta(i,k) + eta(i,k-1))
1461 dp = 1000. * del(i,k)
1463 ptem = c0t(i,k) + c1
1464 qlk = dq / (eta(i,k) + etah * ptem * dz)
1465 dellal(i,k) = etah * c1 * dz * qlk * grav / dp
1467 qlk = dq / (eta(i,k) + etah * c0t(i,k) * dz)
1469 qcko(i,k) = qlk + qrch
1470 pwo(i,k) = etah * c0t(i,k) * dz * qlk
1471 cnvwt(i,k) = etah * qlk * grav / dp
1472 zdqca(i,k)=dq/eta(i,k)
1482 if (hwrf_samfshal)
then
1486 tem = po(i,k) / (rd * to(i,k))
1487 wucb = -0.01 * dot(i,k) / (tem * grav)
1489 wu2(i,k) = wucb * wucb
1499 if(k > kbcon1(i) .and. k < ktcon(i))
then
1500 dz = zi(i,k) - zi(i,k-1)
1501 tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz
1502 tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k))
1503 tem2 = wush(i,k) * sqrt(wu2(i,k-1))
1504 tem2 = (tem1 - tem2) * dz
1505 ptem = (1. - tem) * wu2(i,k-1)
1507 wu2(i,k) = (ptem + tem2) / ptem1
1508 wu2(i,k) = max(wu2(i,k), 0.)
1518 if(k > kbcon1(i) .and. k < ktcon(i))
then
1519 rho = po(i,k)*100. / (rd * to(i,k))
1520 omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav
1521 omega_u(i,k)=max(omega_u(i,k),-80.)
1538 if(k > kbcon1(i) .and. k < ktcon(i))
then
1539 dz = zi(i,k) - zi(i,k-1)
1540 tem = 0.5 * (sqrt(wu2(i,k)) + sqrt(wu2(i,k-1)))
1541 wc(i) = wc(i) + tem * dz
1542 sumx(i) = sumx(i) + dz
1549 if(sumx(i) == 0.)
then
1552 wc(i) = wc(i) / sumx(i)
1555 if (wc(i) < val) cnvflg(i)=.false.
1568 if(k > kbcon1(i) .and. k < ktcon(i))
then
1569 dp = 1000. * del(i,k)
1570 tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1))
1571 omegac(i) = omegac(i) + tem * dp
1572 sumx(i) = sumx(i) + dp
1579 if(sumx(i) == 0.)
then
1582 omegac(i) = omegac(i) / sumx(i)
1585 if (omegac(i) > val) cnvflg(i)=.false.
1593 if(k > kbcon1(i) .and. k < ktcon(i))
then
1594 if(omega_u(i,k) .ne. 0.)
then
1595 zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k))
1599 zeta(i,k)=max(0.,zeta(i,k))
1600 zeta(i,k)=min(1.,zeta(i,k))
1607c exchange ktcon with ktcon1
1612 ktcon(i) = ktcon1(i)
1617c this section is ready for cloud
water
1621c compute liquid and vapor separation at cloud top
1627 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1629 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
1630 dq = qcko(i,k) - qrch
1632c check
if there is excess moisture to release latent heat
1644c--- compute precipitation efficiency in terms of windshear
1680c--- what would the change be, that a cloud with unit mass
1681c--- will
do to the environment?
1687 if(cnvflg(i) .and. k <= kmax(i))
then
1695 if (.not.hwrf_samfshal)
then
1699 if(cnvflg(i) .and. k <= kmax(i))
then
1707c--- changed due to subsidence and entrainment
1712 if(k > kb(i) .and. k < ktcon(i))
then
1713 dp = 1000. * del(i,k)
1714 dz = zi(i,k) - zi(i,k-1)
1717 dv2h = .5 * (heo(i,k) + heo(i,k-1))
1720 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
1725 dellah(i,k) = dellah(i,k) +
1726 & ( eta(i,k)*dv1h - eta(i,k-1)*dv3h
1727 & - tem*eta(i,k-1)*dv2h*dz
1728 & + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz
1731 tem1 = -eta(i,k) * qrcko(i,k)
1732 tem2 = -eta(i,k-1) * qcko(i,k-1)
1733 dellaq(i,k) = dellaq(i,k) + (tem1-tem2) * factor
1735 tem1=eta(i,k)*(uo(i,k)-ucko(i,k))
1736 tem2=eta(i,k-1)*(uo(i,k-1)-ucko(i,k-1))
1737 dellau(i,k) = dellau(i,k) + (tem1-tem2) * factor
1739 tem1=eta(i,k)*(vo(i,k)-vcko(i,k))
1740 tem2=eta(i,k-1)*(vo(i,k-1)-vcko(i,k-1))
1741 dellav(i,k) = dellav(i,k) + (tem1-tem2) * factor
1747 if(.not.hwrf_samfshal)
then
1752 if(k > kb(i) .and. k < ktcon(i))
then
1753 dp = 1000. * del(i,k)
1755 tem1 = -eta(i,k) * ercko(i,k,n)
1756 tem2 = -eta(i,k-1) * ecko(i,k-1,n)
1757 dellae(i,k,n) = dellae(i,k,n) + (tem1-tem2) * grav/dp
1771 dp = 1000. * del(i,indx)
1772 tem = eta(i,indx-1) * grav / dp
1773 dellah(i,indx) = tem * (hcko(i,indx-1) - heo(i,indx-1))
1774 dellaq(i,indx) = tem * qcko(i,indx-1)
1775 dellau(i,indx) = tem * (ucko(i,indx-1) - uo(i,indx-1))
1776 dellav(i,indx) = tem * (vcko(i,indx-1) - vo(i,indx-1))
1780 dellal(i,indx) = tem * qlko_ktcon(i)
1783 if (.not.hwrf_samfshal)
then
1788 dp = 1000. * del(i,indx)
1789 dellae(i,indx,n) = eta(i,indx-1) *
1790 & ecko(i,indx-1,n) * grav / dp
1803 if(cnvflg(i) .and. k <= ktcon(i))
then
1804 q_diff(i,k) = q1(i,k) - q1(i,k+1)
1810 if(q1(i,1) >= 0.)
then
1811 q_diff(i,0) = max(0.,2.*q1(i,1)-q1(i,2))-
1814 q_diff(i,0) = min(0.,2.*q1(i,1)-q1(i,2))-
1824 & (k >= kb(i) .and. k < ktcon(i)))
then
1825 if(eta(i,k) > 0.)
then
1827 if(abs(q_diff(i,k)) > 1.e-22)
1828 & rrkp = q_diff(i,k+1) / q_diff(i,k)
1829 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
1831 & phkp*(qo(i,k)-q1(i,k+1))
1832 flxtvd(i,k) = eta(i,k) * tem1
1841 & (k > kb(i) .and. k <= ktcon(i)))
then
1842 dp = 1000. * del(i,k)
1843 dellaq(i,k) = dellaq(i,k) +
1844 & (flxtvd(i,k) - flxtvd(i,k-1)) * grav/dp
1851 if (.not.hwrf_samfshal)
then
1856 if(cnvflg(i) .and. k <= ktcon(i))
then
1857 e_diff(i,k,n) = ctr(i,k,n) - ctr(i,k+1,n)
1863 if(ctr(i,1,n) >= 0.)
then
1864 e_diff(i,0,n) = max(0.,2.*ctr(i,1,n)-ctr(i,2,n))-
1867 e_diff(i,0,n) = min(0.,2.*ctr(i,1,n)-ctr(i,2,n))-
1880 & (k >= kb(i) .and. k < ktcon(i)))
then
1881 if(eta(i,k) > 0.)
then
1883 if(abs(e_diff(i,k,n)) > 1.e-22)
1884 & rrkp = e_diff(i,k+1,n) / e_diff(i,k,n)
1885 phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp))
1886 tem1 = ctr(i,k+1,n) +
1887 & phkp*(ctro(i,k,n)-ctr(i,k+1,n))
1888 flxtvd(i,k) = eta(i,k) * tem1
1897 & (k > kb(i) .and. k <= ktcon(i)))
then
1898 dp = 1000. * del(i,k)
1899 dellae(i,k,n) = dellae(i,k,n) +
1900 & (flxtvd(i,k) - flxtvd(i,k-1)) * grav/dp
1914 tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
1915 dtconv(i) = tem / wc(i)
1916 if (.not.hwrf_samfshal)
then
1917 tfac = 1. + gdx(i) / 75000.
1918 dtconv(i) = tfac * dtconv(i)
1920 dtconv(i) = max(dtconv(i),dtmin)
1921 dtconv(i) = max(dtconv(i),dt2)
1922 dtconv(i) = min(dtconv(i),dtmax)
1936 if(k >= kbcon1(i) .and. k < ktcon1(i))
then
1937 dz = zi(i,k) - zi(i,k-1)
1938 tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k))
1939 umean(i) = umean(i) + tem * dz
1940 sumx(i) = sumx(i) + dz
1947 umean(i) = umean(i) / sumx(i)
1948 umean(i) = max(umean(i), 1.)
1949 tauadv = gdx(i) / umean(i)
1950 advfac(i) = tauadv / dtconv(i)
1951 advfac(i) = min(advfac(i), 1.)
1955c compute cloud base mass flux as a
function of the mean
1961 if(first_time_step .and. .not.restart)
then
1970 qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt
1977 tmfq(i,k)=tmf(i,k,1)
1981 flag_shallow = .true.
1983 call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
1984 & flag_mid,del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,
1985 & delt,qadv,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu,
1986 & sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
1995 rho = po(i,k)*100. / (rd*to(i,k))
1996 if(
progsigma .and. gdx(i) < dxcrtas)
then
1997 xmb(i) = advfac(i)*sigmab(i)*((-1.0*omegac(i))*gravinv)
1999 xmb(i) = advfac(i)*betaw*rho*wc(i)
2007 tem = min(max(xlamue(i,kbcon(i)), 2.e-4), 6.e-4)
2009 tem1 = 3.14 * tem * tem
2010 sigmagfm(i) = tem1 / garea(i)
2011 sigmagfm(i) = max(sigmagfm(i), 0.001)
2012 sigmagfm(i) = min(sigmagfm(i), 0.999)
2019 if (gdx(i) < dxcrt)
then
2021 scaldfunc(i)=(1.-sigmab(i))*(1.-sigmab(i))
2023 scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i))
2025 scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.)
2029 xmb(i) = xmb(i) * scaldfunc(i)
2030 xmb(i) = min(xmb(i),xmbmax(i))
2053 if (cnvflg(i) .and. k <= kmax(i))
then
2054 qeso(i,k) = 0.01 * fpvs(t1(i,k))
2055 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
2057 qeso(i,k) = max(qeso(i,k), val )
2074 if (.not. hwrf_samfshal)
then
2084 if(k > kb(i) .and. k <= ktcon(i))
then
2085 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
2086 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
2087 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
2091 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2
2092 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2
2093 dp = 1000. * del(i,k)
2094 tem = xmb(i) * dp / grav
2095 delhbar(i) = delhbar(i) + tem * dellah(i,k)
2096 delqbar(i) = delqbar(i) + tem * dellaq(i,k)
2097 deltbar(i) = deltbar(i) + tem * dellat
2098 delubar(i) = delubar(i) + tem * dellau(i,k)
2099 delvbar(i) = delvbar(i) + tem * dellav(i,k)
2116 if(k > kb(i) .and. k <= ktcon(i))
then
2117 tem = q1(i,k) * delp(i,k) / grav
2118 if(q1(i,k) < 0.) tsumn(i) = tsumn(i) + tem
2119 if(q1(i,k) > 0.) tsump(i) = tsump(i) + tem
2126 if(tsump(i) > 0. .and. tsumn(i) < 0.)
then
2127 if(tsump(i) > abs(tsumn(i)))
then
2128 rtnp(i) = tsumn(i) / tsump(i)
2130 rtnp(i) = tsump(i) / tsumn(i)
2138 if(k > kb(i) .and. k <= ktcon(i))
then
2139 if(rtnp(i) < 0.)
then
2140 if(tsump(i) > abs(tsumn(i)))
then
2141 if(q1(i,k) < 0.) q1(i,k)= 0.
2142 if(q1(i,k) > 0.) q1(i,k)=(1.+rtnp(i))*q1(i,k)
2144 if(q1(i,k) < 0.) q1(i,k)=(1.+rtnp(i))*q1(i,k)
2145 if(q1(i,k) > 0.) q1(i,k)=0.
2153 if (.not.hwrf_samfshal)
then
2161 if(k > kb(i) .and. k <= ktcon(i))
then
2162 ctr(i,k,n) = ctr(i,k,n)+dellae(i,k,n)*xmb(i)*dt2
2163 dp = 1000. * del(i,k)
2164 delebar(i,n)=delebar(i,n)+dellae(i,k,n)*xmb(i)*dp/grav
2181 if(k > kb(i) .and. k <= ktcon(i))
then
2184 dz = zi(i,k) - zi(i,k-1)
2188 tem = ctr(i,k,n) * dz
2190 tem = ctr(i,k,n) * delp(i,k) / grav
2192 if(ctr(i,k,n) < 0.) tsumn(i) = tsumn(i) + tem
2193 if(ctr(i,k,n) > 0.) tsump(i) = tsump(i) + tem
2200 if(tsump(i) > 0. .and. tsumn(i) < 0.)
then
2201 if(tsump(i) > abs(tsumn(i)))
then
2202 rtnp(i) = tsumn(i) / tsump(i)
2204 rtnp(i) = tsump(i) / tsumn(i)
2212 if(k > kb(i) .and. k <= ktcon(i))
then
2213 if(rtnp(i) < 0.)
then
2214 if(tsump(i) > abs(tsumn(i)))
then
2215 if(ctr(i,k,n)<0.) ctr(i,k,n)=0.
2216 if(ctr(i,k,n)>0.) ctr(i,k,n)=(1.+rtnp(i))*ctr(i,k,n)
2218 if(ctr(i,k,n)<0.) ctr(i,k,n)=(1.+rtnp(i))*ctr(i,k,n)
2219 if(ctr(i,k,n)>0.) ctr(i,k,n)=0.
2231 if(k > kb(i) .and. k <= ktcon(i))
then
2232 qtr(i,k,kk) = ctr(i,k,n)
2240 if (do_aerosols)
then
2248 if(k > kb(i) .and. k < ktcon(i))
then
2249 dp = 1000. * del(i,k)
2250 wet_dep(i,k,n) = chem_pw(i,k,n)*grav/dp
2251 wet_dep(i,k,n) = wet_dep(i,k,n)*xmb(i)*dt2*dp
2261 if(k > kb(i) .and. k < ktcon(i))
then
2262 dp = 1000. * del(i,k)
2263 if (qtr(i,k,kk) < 0.)
then
2265 tem = -qtr(i,k,kk)*dp
2266 if(wet_dep(i,k,n) >= tem)
then
2267 wet_dep(i,k,n) = wet_dep(i,k,n) - tem
2271 qtr(i,k,kk) = qtr(i,k,kk)+wet_dep(i,k,n)/dp
2289 if(k > kb(i) .and. k <= ktcon(i))
then
2290 qeso(i,k) = 0.01 * fpvs(t1(i,k))
2291 qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
2293 qeso(i,k) = max(qeso(i,k), val )
2309 if(k < ktcon(i) .and. k > kb(i))
then
2310 rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
2323 if (k <= kmax(i))
then
2328 if(k < ktcon(i) .and. k > kb(i))
then
2329 rn(i) = rn(i) + pwo(i,k) * xmb(i) * .001 * dt2
2332 if(flg(i) .and. k < ktcon(i))
then
2336 qcond(i) = shevf * evef * (q1(i,k) - qeso(i,k))
2337 & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
2338 dp = 1000. * del(i,k)
2340 if(rn(i) > 0. .and. qcond(i) < 0.)
then
2341 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
2342 qevap(i) = min(qevap(i), rn(i)*1000.*grav/dp)
2343 delq2(i) = delqev(i) + .001 * qevap(i) * factor
2345 if(rn(i) > 0. .and. qcond(i) < 0. .and.
2346 & delq2(i) > rntot(i))
then
2347 qevap(i) = 1000.* grav * (rntot(i) - delqev(i)) / dp
2350 if(rn(i) > 0. .and. qevap(i) > 0.)
then
2352 tem1 = qevap(i) * tem
2353 if(tem1 > rn(i))
then
2354 qevap(i) = rn(i) / tem
2357 rn(i) = rn(i) - tem1
2359 q1(i,k) = q1(i,k) + qevap(i)
2360 t1(i,k) = t1(i,k) - elocp * qevap(i)
2361 deltv(i) = - elocp*qevap(i)/dt2
2362 delq(i) = + qevap(i)/dt2
2363 delqev(i) = delqev(i) + tem * qevap(i)
2365 delqbar(i) = delqbar(i) + delq(i) * factor
2366 deltbar(i) = deltbar(i) + deltv(i) * factor
2393 if(rn(i) < 0. .or. .not.flg(i)) rn(i) = 0.
2400c convective cloud
water
2406 if (k >= kbcon(i) .and. k < ktcon(i))
then
2407 cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
2414c convective cloud cover
2420 if (k >= kbcon(i) .and. k < ktcon(i))
then
2421 cnvc(i,k) = 0.04 * log(1. + 675. * eta(i,k) * xmb(i))
2422 cnvc(i,k) = min(cnvc(i,k), 0.2)
2423 cnvc(i,k) = max(cnvc(i,k), 0.0)
2432 if (ncloud > 0)
then
2438 if (k >= kbcon(i) .and. k <= ktcon(i))
then
2439 tem = dellal(i,k) * xmb(i) * dt2
2440 tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
2441 if (qtr(i,k,2) > -999.0)
then
2442 qtr(i,k,1) = qtr(i,k,1) + tem * tem1
2443 qtr(i,k,2) = qtr(i,k,2) + tem *(1.0-tem1)
2445 qtr(i,k,1) = qtr(i,k,1) + tem
2477 if(k >= kb(i) .and. k < ktop(i))
then
2478 ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
2487 dt_mf(i,k) = ud_mf(i,k)
2493 if (.not.hwrf_samfshal)
then
2499 if(k > kb(i) .and. k < ktop(i))
then
2500 tem = 0.5 * (eta(i,k-1) + eta(i,k)) * xmb(i)
2501 tem1 = pfld(i,k) * 100. / (rd * t1(i,k))
2505 tem2 = max(sigmagfm(i), betaw)
2507 ptem = tem / (tem2 * tem1)
2508 qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*tem2*ptem*ptem