90 subroutine shalcnv_run( &
91 & grav,cp,hvap,rv,fv,t0c,rd,cvap,cliq,eps,epsm1, &
92 & im,km,jcap,delt,delp,prslp,psp,phil,qlc,qli, &
93 & q1,t1,u1,v1,rn,kbot,ktop,kcnv,islimsk, &
94 & dot,ncloud,hpbl,heat,evap,ud_mf,dt_mf,cnvw,cnvc, &
95 & clam,c0,c1,pgcon,errmsg,errflg)
107 real(kind=kind_phys),
intent(in) :: grav, cp, hvap, rv, fv, t0c, &
108 & rd, cvap, cliq, eps, epsm1
109 integer,
intent(in) :: im, km, jcap, ncloud
110 integer,
intent(inout) :: kbot(:), ktop(:), kcnv(:)
111 integer,
intent(in) :: islimsk(:)
112 real(kind=kind_phys),
intent(in) :: delt, clam, c0, c1, pgcon
113 real(kind=kind_phys),
intent(in) :: psp(:), delp(:,:), &
114 & prslp(:,:), dot(:,:), &
115 & phil(:,:), hpbl(:), &
117 real(kind=kind_phys),
intent(inout) :: &
118 & qlc(:,:), qli(:,:), &
119 & q1(:,:), t1(:,:), &
120 & u1(:,:), v1(:,:), &
121 & cnvw(:,:), cnvc(:,:)
122 real(kind=kind_phys),
intent(out) :: rn(:), dt_mf(:,:)
123 real(kind=kind_phys),
intent(out),
optional :: ud_mf(:,:)
124 character(len=*),
intent(out) :: errmsg
125 integer,
intent(out) :: errflg
129 integer i,j,indx, k, kk, km1
132 real(kind=kind_phys) dellat, delta,
135 & dq, dqsdp, dqsdt, dt,
136 & dt2, dtmax, dtmin, dv1h,
137 & dv1q, dv2h, dv2q, dv1u,
138 & dv1v, dv2u, dv2v, dv3q,
141 & el2orc, elocp, aafac,
142 & es, etah, h1, dthk,
143 & evef, evfact, evfactl, fact1,
144 & fact2, factor, fjcap,
145 & g, gamma, pprime, betaw,
147 & rfact, shear, tem1,
149 & val2, w1, w1l, w1s,
152 & w4s, tem, ptem, ptem1
154 integer kb(im), kbcon(im), kbcon1(im),
155 & ktcon(im), ktcon1(im),
158 real(kind=kind_phys) aa1(im),
159 & delhbar(im), delq(im), delq2(im),
160 & delqbar(im), delqev(im), deltbar(im),
161 & deltv(im), edt(im),
162 & wstar(im),
sflx(im),
163 & pdot(im), po(im,km),
164 & qcond(im), qevap(im), hmax(im),
165 & rntot(im), vshear(im),
166 & xlamud(im), xmb(im), xmbmax(im),
167 & delubar(im), delvbar(im),
168 & ps(im), del(im,km), prsl(im,km)
170 real(kind=kind_phys) cincr, cincrmax, cincrmin
179 parameter(cincrmax=180.,cincrmin=120.,dthk=25.)
180 parameter(h1=0.33333333)
182 real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km),
183 & uo(im,km), vo(im,km), qeso(im,km)
186 real(kind=kind_phys) qlko_ktcon(im), dellal(im,km),
187 & dbyo(im,km), zo(im,km), xlamue(im,km),
188 & heo(im,km), heso(im,km),
189 & dellah(im,km), dellaq(im,km),
190 & dellau(im,km), dellav(im,km), hcko(im,km),
191 & ucko(im,km), vcko(im,km), qcko(im,km),
192 & qrcko(im,km), eta(im,km),
193 & zi(im,km), pwo(im,km),
194 & tx1(im), cnvwt(im,km)
196 logical totflg, cnvflg(im), flg(im)
198 real(kind=kind_phys) tf, tcr, tcrf
199 parameter(tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
207 el2orc = hvap*hvap/(rv*cp)
209 fact1 = (cvap-cliq)/rv
210 fact2 = hvap/rv-fact1*t0c
234 sflx(i) = heat(i)+fv*t1(i,1)*evap(i)
242 if(kcnv(i).eq.1) cnvflg(i) = .false.
243 if(
sflx(i).le.0.) cnvflg(i) = .false.
270 totflg = totflg .and. (.not. cnvflg(i))
278 dtmin = max(dt2, val )
280 dtmax = max(dt2, val )
291 fjcap = (float(jcap) / 126.) ** 2
293 fjcap = max(fjcap,val)
315 if (prsl(i,k)*tx1(i) .gt. 0.70) kbm(i) = k + 1
316 if (prsl(i,k)*tx1(i) .gt. 0.60) kmax(i) = k + 1
320 kbm(i) = min(kbm(i),kmax(i))
329 zo(i,k) = phil(i,k) / g
335 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
336 xlamue(i,k) = clam / zi(i,k)
340 xlamue(i,km) = xlamue(i,km1)
352 if (flg(i).and.zo(i,k).le.hpbl(i))
then
360 kpbl(i)= min(kpbl(i),kbm(i))
369 if (cnvflg(i) .and. k .le. kmax(i))
then
370 pfld(i,k) = prsl(i,k) * 10.0
401 if (cnvflg(i) .and. k .le. kmax(i))
then
402 qeso(i,k) = 0.01 * fpvs(to(i,k))
403 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
405 qeso(i,k) = max(qeso(i,k), val1)
407 qo(i,k) = max(qo(i,k), val2 )
419 if (cnvflg(i) .and. k .le. kmax(i))
then
421 tem = phil(i,k) + cp * to(i,k)
422 heo(i,k) = tem + hvap * qo(i,k)
423 heso(i,k) = tem + hvap * qeso(i,k)
442 if (cnvflg(i).and.k.le.kpbl(i))
then
443 if(heo(i,k).gt.hmax(i))
then
454 if (cnvflg(i) .and. k .le. kmax(i)-1)
then
455 dz = .5 * (zo(i,k+1) - zo(i,k))
456 dp = .5 * (pfld(i,k+1) - pfld(i,k))
457 es = 0.01 * fpvs(to(i,k+1))
458 pprime = pfld(i,k+1) + epsm1 * es
459 qs = eps * es / pprime
460 dqsdp = - qs / pprime
461 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
462 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
463 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
464 dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
465 dq = dqsdt * dt + dqsdp * dp
466 to(i,k) = to(i,k+1) + dt
467 qo(i,k) = qo(i,k+1) + dq
468 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
476 if (cnvflg(i) .and. k .le. kmax(i)-1)
then
477 qeso(i,k) = 0.01 * fpvs(to(i,k))
478 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
480 qeso(i,k) = max(qeso(i,k), val1)
482 qo(i,k) = max(qo(i,k), val2 )
484 heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
485 & cp * to(i,k) + hvap * qo(i,k)
486 heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
487 & cp * to(i,k) + hvap * qeso(i,k)
488 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
489 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
498 if(flg(i)) kbcon(i) = kmax(i)
502 if (flg(i).and.k.lt.kbm(i))
then
503 if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k))
then
514 if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
520 totflg = totflg .and. (.not. cnvflg(i))
532 pdot(i) = 0.01 * dot(i,kbcon(i))
537 if(islimsk(i) == 1)
then
548 if(pdot(i).le.w4)
then
549 ptem = (pdot(i) - w4) / (w3 - w4)
550 elseif(pdot(i).ge.-w4)
then
551 ptem = - (pdot(i) + w4) / (w4 - w3)
556 ptem = max(ptem,val1)
558 ptem = min(ptem,val2)
560 ptem1= .5*(cincrmax-cincrmin)
561 cincr = cincrmax - ptem * ptem1
562 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
563 if(tem1.gt.cincr)
then
571 totflg = totflg .and. (.not. cnvflg(i))
582 xlamud(i) = xlamue(i,kbcon(i))
596 if(k.lt.kbcon(i).and.k.ge.kb(i))
then
597 dz = zi(i,k+1) - zi(i,k)
598 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
599 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
610 if(k.gt.kbcon(i).and.k.lt.kmax(i))
then
611 dz = zi(i,k) - zi(i,k-1)
612 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
613 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
625 hcko(i,indx) = heo(i,indx)
626 ucko(i,indx) = uo(i,indx)
627 vcko(i,indx) = vo(i,indx)
635 if(k.gt.kb(i).and.k.lt.kmax(i))
then
636 dz = zi(i,k) - zi(i,k-1)
637 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
638 tem1 = 0.5 * xlamud(i) * dz
639 factor = 1. + tem - tem1
640 ptem = 0.5 * tem + pgcon
641 ptem1= 0.5 * tem - pgcon
642 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
643 & (heo(i,k)+heo(i,k-1)))/factor
644 ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k)
645 & +ptem1*uo(i,k-1))/factor
646 vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k)
647 & +ptem1*vo(i,k-1))/factor
648 dbyo(i,k) = hcko(i,k) - heso(i,k)
664 if (flg(i).and.k.lt.kbm(i))
then
665 if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.)
then
674 if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
679 tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
688 totflg = totflg .and. (.not. cnvflg(i))
699 if(flg(i)) ktcon(i) = kbm(i)
703 if (flg(i).and.k .lt. kbm(i))
then
704 if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.)
then
736 dp = 1000. * del(i,k)
737 xmbmax(i) = dp / (g * dt2)
750 qcko(i,kb(i)) = qo(i,kb(i))
751 qrcko(i,kb(i)) = qo(i,kb(i))
758 if(k.gt.kb(i).and.k.lt.ktcon(i))
then
759 dz = zi(i,k) - zi(i,k-1)
760 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
762 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
764 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
765 tem1 = 0.5 * xlamud(i) * dz
766 factor = 1. + tem - tem1
767 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
768 & (qo(i,k)+qo(i,k-1)))/factor
769 qrcko(i,k) = qcko(i,k)
771 dq = eta(i,k) * (qcko(i,k) - qrch)
777 if(k.ge.kbcon(i).and.dq.gt.0.)
then
778 etah = .5 * (eta(i,k) + eta(i,k-1))
779 if(ncloud.gt.0.)
then
780 dp = 1000. * del(i,k)
781 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
782 dellal(i,k) = etah * c1 * dz * qlk * g / dp
784 qlk = dq / (eta(i,k) + etah * c0 * dz)
786 aa1(i) = aa1(i) - dz * g * qlk
787 qcko(i,k)= qlk + qrch
788 pwo(i,k) = etah * c0 * dz * qlk
789 cnvwt(i,k) = etah * qlk * g / dp
806 if(k.ge.kbcon(i).and.k.lt.ktcon(i))
then
807 dz1 = zo(i,k+1) - zo(i,k)
808 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
809 rfact = 1. + delta * cp * gamma
812 & dz1 * (g / (cp * to(i,k)))
813 & * dbyo(i,k) / (1. + gamma)
818 & max(val,(qeso(i,k) - qo(i,k)))
825 if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
830 totflg = totflg .and. (.not. cnvflg(i))
843 aa1(i) = aafac * aa1(i)
854 if(k.ge.ktcon(i).and.k.lt.kbm(i))
then
855 dz1 = zo(i,k+1) - zo(i,k)
856 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
857 rfact = 1. + delta * cp * gamma
860 & dz1 * (g / (cp * to(i,k)))
861 & * dbyo(i,k) / (1. + gamma)
863 if(aa1(i).lt.0.)
then
879 if(k.ge.ktcon(i).and.k.lt.ktcon1(i))
then
880 dz = zi(i,k) - zi(i,k-1)
881 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
883 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
885 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
886 tem1 = 0.5 * xlamud(i) * dz
887 factor = 1. + tem - tem1
888 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
889 & (qo(i,k)+qo(i,k-1)))/factor
890 qrcko(i,k) = qcko(i,k)
892 dq = eta(i,k) * (qcko(i,k) - qrch)
897 etah = .5 * (eta(i,k) + eta(i,k-1))
898 if(ncloud.gt.0.)
then
899 dp = 1000. * del(i,k)
900 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
901 dellal(i,k) = etah * c1 * dz * qlk * g / dp
903 qlk = dq / (eta(i,k) + etah * c0 * dz)
905 qcko(i,k) = qlk + qrch
906 pwo(i,k) = etah * c0 * dz * qlk
907 cnvwt(i,k) = etah * qlk * g / dp
934 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
936 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
937 dq = qcko(i,k) - qrch
964 if(k.gt.kb(i).and.k.le.ktcon(i))
then
965 shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2
966 & + (vo(i,k)-vo(i,k-1)) ** 2)
967 vshear(i) = vshear(i) + shear
974 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
975 e1=1.591-.639*vshear(i)
976 & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
979 edt(i) = min(edt(i),val)
981 edt(i) = max(edt(i),val)
992 if(cnvflg(i) .and. k .le. kmax(i))
then
1006 if(k.gt.kb(i).and.k.lt.ktcon(i))
then
1007 dp = 1000. * del(i,k)
1008 dz = zi(i,k) - zi(i,k-1)
1011 dv2h = .5 * (heo(i,k) + heo(i,k-1))
1014 dv2q = .5 * (qo(i,k) + qo(i,k-1))
1017 dv2u = .5 * (uo(i,k) + uo(i,k-1))
1020 dv2v = .5 * (vo(i,k) + vo(i,k-1))
1023 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
1026 dellah(i,k) = dellah(i,k) +
1027 & ( eta(i,k)*dv1h - eta(i,k-1)*dv3h
1028 & - tem*eta(i,k-1)*dv2h*dz
1029 & + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz
1032 dellaq(i,k) = dellaq(i,k) +
1033 & ( eta(i,k)*dv1q - eta(i,k-1)*dv3q
1034 & - tem*eta(i,k-1)*dv2q*dz
1035 & + tem1*eta(i,k-1)*.5*(qrcko(i,k)+qcko(i,k-1))*dz
1038 dellau(i,k) = dellau(i,k) +
1039 & ( eta(i,k)*dv1u - eta(i,k-1)*dv3u
1040 & - tem*eta(i,k-1)*dv2u*dz
1041 & + tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz
1042 & - pgcon*eta(i,k-1)*(dv1u-dv3u)
1045 dellav(i,k) = dellav(i,k) +
1046 & ( eta(i,k)*dv1v - eta(i,k-1)*dv3v
1047 & - tem*eta(i,k-1)*dv2v*dz
1048 & + tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz
1049 & - pgcon*eta(i,k-1)*(dv1v-dv3v)
1062 dp = 1000. * del(i,indx)
1063 dv1h = heo(i,indx-1)
1064 dellah(i,indx) = eta(i,indx-1) *
1065 & (hcko(i,indx-1) - dv1h) * g / dp
1067 dellaq(i,indx) = eta(i,indx-1) *
1068 & (qcko(i,indx-1) - dv1q) * g / dp
1070 dellau(i,indx) = eta(i,indx-1) *
1071 & (ucko(i,indx-1) - dv1u) * g / dp
1073 dellav(i,indx) = eta(i,indx-1) *
1074 & (vcko(i,indx-1) - dv1v) * g / dp
1078 dellal(i,indx) = eta(i,indx-1) *
1079 & qlko_ktcon(i) * g / dp
1095 ptem = g*
sflx(i)*hpbl(i)/t1(i,1)
1097 tem = po(i,k)*100. / (rd*t1(i,k))
1098 xmb(i) = betaw*tem*wstar(i)
1099 xmb(i) = min(xmb(i),xmbmax(i))
1109 if (cnvflg(i) .and. k .le. kmax(i))
then
1110 qeso(i,k) = 0.01 * fpvs(t1(i,k))
1111 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
1113 qeso(i,k) = max(qeso(i,k), val )
1133 if(k.gt.kb(i).and.k.le.ktcon(i))
then
1134 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
1135 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
1136 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
1140 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2
1141 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2
1142 dp = 1000. * del(i,k)
1143 delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
1144 delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
1145 deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
1146 delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
1147 delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
1156 if(k.gt.kb(i).and.k.le.ktcon(i))
then
1157 qeso(i,k) = 0.01 * fpvs(t1(i,k))
1158 qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
1160 qeso(i,k) = max(qeso(i,k), val )
1176 if(k.lt.ktcon(i).and.k.gt.kb(i))
then
1177 rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
1190 if (k .le. kmax(i))
then
1195 if(k.lt.ktcon(i).and.k.gt.kb(i))
then
1196 rn(i) = rn(i) + pwo(i,k) * xmb(i) * .001 * dt2
1199 if(flg(i).and.k.lt.ktcon(i))
then
1200 evef = edt(i) * evfact
1201 if(islimsk(i) == 1) evef=edt(i) * evfactl
1204 qcond(i) = evef * (q1(i,k) - qeso(i,k))
1205 & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
1206 dp = 1000. * del(i,k)
1207 if(rn(i).gt.0..and.qcond(i).lt.0.)
then
1208 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
1209 qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
1210 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
1212 if(rn(i).gt.0..and.qcond(i).lt.0..and.
1213 & delq2(i).gt.rntot(i))
then
1214 qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
1217 if(rn(i).gt.0..and.qevap(i).gt.0.)
then
1219 tem1 = qevap(i) * tem
1220 if(tem1.gt.rn(i))
then
1221 qevap(i) = rn(i) / tem
1224 rn(i) = rn(i) - tem1
1226 q1(i,k) = q1(i,k) + qevap(i)
1227 t1(i,k) = t1(i,k) - elocp * qevap(i)
1228 deltv(i) = - elocp*qevap(i)/dt2
1229 delq(i) = + qevap(i)/dt2
1230 delqev(i) = delqev(i) + .001*dp*qevap(i)/g
1232 dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
1233 delqbar(i) = delqbar(i) + delq(i)*dp/g
1234 deltbar(i) = deltbar(i) + deltv(i)*dp/g
1253 if(rn(i).lt.0..or..not.flg(i)) rn(i) = 0.
1265 if (cnvflg(i) .and. rn(i).gt.0.)
then
1266 if (k.ge.kbcon(i).and.k.lt.ktcon(i))
then
1267 cnvw(i,k) = cnvwt(i,k) * xmb(i) * dt2
1279 if (cnvflg(i) .and. rn(i).gt.0.)
then
1280 if (k.ge.kbcon(i).and.k.lt.ktcon(i))
then
1281 cnvc(i,k) = 0.04 * log(1. + 675. * eta(i,k) * xmb(i))
1282 cnvc(i,k) = min(cnvc(i,k), 0.2)
1283 cnvc(i,k) = max(cnvc(i,k), 0.0)
1293 if (ncloud.gt.0)
then
1298 if (k.gt.kb(i).and.k.le.ktcon(i))
then
1299 tem = dellal(i,k) * xmb(i) * dt2
1300 tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
1301 if (qlc(i,k) .gt. -999.0)
then
1302 qli(i,k) = qli(i,k) + tem * tem1
1303 qlc(i,k) = qlc(i,k) + tem *(1.0-tem1)
1305 qli(i,k) = qli(i,k) + tem
1320 if(k.ge.kb(i) .and. k.lt.ktop(i))
then
1321 ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
1330 dt_mf(i,k) = ud_mf(i,k)