subroutine sascnvn(im,ix,km,jcap,delt,del,prsl,ps,phil,ql,
& q1,t1,u1,v1,rcs,cldwrk,rn,kbot,ktop,kcnv,slimsk,
& dot,ncloud,ud_mf,dd_mf,dt_mf)
! & dot,ncloud,ud_mf,dd_mf,dt_mf,me)
!
use machine , only : kind_phys
use funcphys , only : fpvs
use physcons, grav => con_g, cp => con_cp, hvap => con_hvap
&, rv => con_rv, fv => con_fvirt, t0c => con_t0c
&, cvap => con_cvap, cliq => con_cliq
&, eps => con_eps, epsm1 => con_epsm1
implicit none
!
integer im, ix, km, jcap, ncloud,
& kbot(im), ktop(im), kcnv(im)
! &, me
real(kind=kind_phys) delt
real(kind=kind_phys) ps(im), del(ix,km), prsl(ix,km),
& ql(ix,km,2),q1(ix,km), t1(ix,km),
& u1(ix,km), v1(ix,km), rcs(im),
& cldwrk(im), rn(im), slimsk(im),
& dot(ix,km), phil(ix,km)
! hchuang code change mass flux output
&, ud_mf(ix,km),dd_mf(ix,km),dt_mf(ix,km)
!
integer i, indx, jmn, k, latd, lond, km1
!
real(kind=kind_phys) clam, cxlamu, xlamde, xlamdd
!
real(kind=kind_phys) adw, alpha, alphal, alphas,
& aup, beta, betal, betas,
& c0, cpoel, dellat, delta,
& desdt, deta, detad, dg,
& dh, dhh, dlnsig, dp,
& dq, dqsdp, dqsdt, dt,
& dt2, dtmax, dtmin, dv1h,
& dv1q, dv2h, dv2q, dv1u,
& dv1v, dv2u, dv2v, dv3q,
& dv3h, dv3u, dv3v,
& dz, dz1, e1, edtmax,
& edtmaxl, edtmaxs, el2orc, elocp,
& es, etah, cthk, dthk,
& evef, evfact, evfactl, fact1,
& fact2, factor, fjcap, fkm,
& g, gamma, pprime, c1,
& qc, qlk, qrch, qs,
& rain, rfact, shear, tem1,
& tem2, terr, val, val1,
& val2, w1, w1l, w1s,
& w2, w2l, w2s, w3,
& w3l, w3s, w4, w4l,
& w4s, xdby, xpw, xpwd,
& xqc, xqrch, mbdt, tem,
& ptem, ptem1, pgcon
!
integer kb(im), kbcon(im), kbcon1(im),
& ktcon(im), jmin(im),
& lmin(im), kbmax(im),kbm(im), kmax(im)
!
real(kind=kind_phys) aa1(im), acrt(im), acrtfct(im),
& delhbar(im), delq(im), delq2(im),
& delqbar(im), delqev(im), deltbar(im),
& deltv(im), dtconv(im), edt(im),
& edto(im), edtx(im), fld(im),
& hcdo(im,km), hmax(im), hmin(im),
& ucdo(im,km), vcdo(im,km),
& pbcdif(im), pdot(im), po(im,km),
& pwavo(im), pwevo(im), xlamud(im),
& qcdo(im,km), qcond(im), qevap(im),
& rntot(im), vshear(im), xaa0(im),
& xk(im), xlamd(im),
& xmb(im), xmbmax(im), xpwav(im),
& xpwev(im), delubar(im),delvbar(im)
cj
real(kind=kind_phys) cincr, cincrmax, cincrmin
cj
c physical parameters
parameter(g=grav)
parameter(cpoel=cp/hvap,elocp=hvap/cp,
& el2orc=hvap*hvap/(rv*cp))
parameter(terr=0.,c0=.002,c1=3.e-4,delta=fv)
parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
parameter(cthk=150.,cincrmax=180.,cincrmin=120.,dthk=25.)
c local variables and arrays
real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km),
& uo(im,km), vo(im,km), qeso(im,km)
c cloud water
real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km),
& dbyo(im,km), zo(im,km), xlamue(im,km),
& fent1(im,km), fent2(im,km), frh(im,km),
& heo(im,km), heso(im,km),
& qrcd(im,km), dellah(im,km), dellaq(im,km),
& dellau(im,km), dellav(im,km), hcko(im,km),
& ucko(im,km), vcko(im,km), qcko(im,km),
& eta(im,km), etad(im,km), zi(im,km),
& qrcdo(im,km), pwo(im,km), pwdo(im,km),
& tx1(im), sumx(im)
! &, rhbar(im)
!
logical totflg, cnvflg(im), flg(im)
!
real(kind=kind_phys) pcrit(15), acritt(15), acrit(15)
! save pcrit, acritt
data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,
& 350.,300.,250.,200.,150./
data acritt/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216,
& .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
c gdas derived acrit
c data acritt/.203,.515,.521,.566,.625,.665,.659,.688,
c & .743,.813,.886,.947,1.138,1.377,1.896/
real(kind=kind_phys) tf, tcr, tcrf
parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
!
c-----------------------------------------------------------------------
!
km1 = km - 1
c
c initialize arrays
c
do i=1,im
cnvflg(i) = .true.
rn(i)=0.
kbot(i)=km+1
ktop(i)=0
kbcon(i)=km
ktcon(i)=1
dtconv(i) = 3600.
cldwrk(i) = 0.
pdot(i) = 0.
pbcdif(i)= 0.
lmin(i) = 1
jmin(i) = 1
qlko_ktcon(i) = 0.
edt(i) = 0.
edto(i) = 0.
edtx(i) = 0.
acrt(i) = 0.
acrtfct(i) = 1.
aa1(i) = 0.
xaa0(i) = 0.
pwavo(i)= 0.
pwevo(i)= 0.
xpwav(i)= 0.
xpwev(i)= 0.
vshear(i) = 0.
enddo
! hchuang code change
do k = 1, km
do i = 1, im
ud_mf(i,k) = 0.
dd_mf(i,k) = 0.
dt_mf(i,k) = 0.
enddo
enddo
c
do k = 1, 15
acrit(k) = acritt(k) * (975. - pcrit(k))
enddo
dt2 = delt
val = 1200.
dtmin = max(dt2, val )
val = 3600.
dtmax = max(dt2, val )
c model tunable parameters are all here
mbdt = 10.
edtmaxl = .3
edtmaxs = .3
alphal = .5
alphas = .5
clam = .1
! betal = .15
! betas = .15
betal = .05
betas = .05
c evef = 0.07
evfact = 0.3
evfactl = 0.3
!
cxlamu = 1.0e-4
xlamde = 1.0e-4
xlamdd = 1.0e-4
!
! pgcon = 0.7 ! Gregory et al. (1997, QJRMS)
pgcon = 0.55 ! Zhang & Wu (2003,JAS)
fjcap = (float(jcap) / 126.) ** 2
val = 1.
fjcap = max(fjcap,val)
fkm = (float(km) / 28.) ** 2
fkm = max(fkm,val)
w1l = -8.e-3
w2l = -4.e-2
w3l = -5.e-3
w4l = -5.e-4
w1s = -2.e-4
w2s = -2.e-3
w3s = -1.e-3
w4s = -2.e-5
c
c define top layer for search of the downdraft originating layer
c and the maximum thetae for updraft
c
do i=1,im
kbmax(i) = km
kbm(i) = km
kmax(i) = km
tx1(i) = 1.0 / ps(i)
enddo
!
do k = 1, km
do i=1,im
if (prsl(i,k)*tx1(i) .gt. 0.04) kmax(i) = k + 1
if (prsl(i,k)*tx1(i) .gt. 0.45) kbmax(i) = k + 1
if (prsl(i,k)*tx1(i) .gt. 0.70) kbm(i) = k + 1
enddo
enddo
do i=1,im
kbmax(i) = min(kbmax(i),kmax(i))
kbm(i) = min(kbm(i),kmax(i))
enddo
c
c hydrostatic height assume zero terr and initially assume
c updraft entrainment rate as an inverse function of height
c
do k = 1, km
do i=1,im
zo(i,k) = phil(i,k) / g
enddo
enddo
do k = 1, km1
do i=1,im
zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
xlamue(i,k) = clam / zi(i,k)
enddo
enddo
c
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c convert surface pressure to mb from cb
c
do k = 1, km
do i = 1, im
if (k .le. kmax(i)) then
pfld(i,k) = prsl(i,k) * 10.0
eta(i,k) = 1.
fent1(i,k)= 1.
fent2(i,k)= 1.
frh(i,k) = 0.
hcko(i,k) = 0.
qcko(i,k) = 0.
ucko(i,k) = 0.
vcko(i,k) = 0.
etad(i,k) = 1.
hcdo(i,k) = 0.
qcdo(i,k) = 0.
ucdo(i,k) = 0.
vcdo(i,k) = 0.
qrcd(i,k) = 0.
qrcdo(i,k)= 0.
dbyo(i,k) = 0.
pwo(i,k) = 0.
pwdo(i,k) = 0.
dellal(i,k) = 0.
to(i,k) = t1(i,k)
qo(i,k) = q1(i,k)
uo(i,k) = u1(i,k) * rcs(i)
vo(i,k) = v1(i,k) * rcs(i)
endif
enddo
enddo
c
c column variables
c p is pressure of the layer (mb)
c t is temperature at t-dt (k)..tn
c q is mixing ratio at t-dt (kg/kg)..qn
c to is temperature at t+dt (k)... this is after advection and turbulan
c qo is mixing ratio at t+dt (kg/kg)..q1
c
do k = 1, km
do i=1,im
if (k .le. kmax(i)) then
qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
val1 = 1.e-8
qeso(i,k) = max(qeso(i,k), val1)
val2 = 1.e-10
qo(i,k) = max(qo(i,k), val2 )
! qo(i,k) = min(qo(i,k),qeso(i,k))
! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
endif
enddo
enddo
c
c compute moist static energy
c
do k = 1, km
do i=1,im
if (k .le. kmax(i)) then
! tem = g * zo(i,k) + cp * to(i,k)
tem = phil(i,k) + cp * to(i,k)
heo(i,k) = tem + hvap * qo(i,k)
heso(i,k) = tem + hvap * qeso(i,k)
c heo(i,k) = min(heo(i,k),heso(i,k))
endif
enddo
enddo
c
c determine level with largest moist static energy
c this is the level where updraft starts
c
do i=1,im
hmax(i) = heo(i,1)
kb(i) = 1
enddo
do k = 2, km
do i=1,im
if (k .le. kbm(i)) then
if(heo(i,k).gt.hmax(i)) then
kb(i) = k
hmax(i) = heo(i,k)
endif
endif
enddo
enddo
c
do k = 1, km1
do i=1,im
if (k .le. kmax(i)-1) then
dz = .5 * (zo(i,k+1) - zo(i,k))
dp = .5 * (pfld(i,k+1) - pfld(i,k))
es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
pprime = pfld(i,k+1) + epsm1 * es
qs = eps * es / pprime
dqsdp = - qs / pprime
desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
dq = dqsdt * dt + dqsdp * dp
to(i,k) = to(i,k+1) + dt
qo(i,k) = qo(i,k+1) + dq
po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
endif
enddo
enddo
!
do k = 1, km1
do i=1,im
if (k .le. kmax(i)-1) then
qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
val1 = 1.e-8
qeso(i,k) = max(qeso(i,k), val1)
val2 = 1.e-10
qo(i,k) = max(qo(i,k), val2 )
! qo(i,k) = min(qo(i,k),qeso(i,k))
frh(i,k) = 1. - min(qo(i,k)/qeso(i,k), 1.)
heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
& cp * to(i,k) + hvap * qo(i,k)
heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
& cp * to(i,k) + hvap * qeso(i,k)
uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
endif
enddo
enddo
c
c look for the level of free convection as cloud base
c
do i=1,im
flg(i) = .true.
kbcon(i) = kmax(i)
enddo
do k = 1, km1
do i=1,im
if (flg(i).and.k.le.kbmax(i)) then
if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
kbcon(i) = k
flg(i) = .false.
endif
endif
enddo
enddo
c
do i=1,im
if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
enddo
!!
totflg = .true.
do i=1,im
totflg = totflg .and. (.not. cnvflg(i))
enddo
if(totflg) return
!!
c
c determine critical convective inhibition
c as a function of vertical velocity at cloud base.
c
do i=1,im
if(cnvflg(i)) then
pdot(i) = 10.* dot(i,kbcon(i))
endif
enddo
do i=1,im
if(cnvflg(i)) then
if(slimsk(i).eq.1.) then
w1 = w1l
w2 = w2l
w3 = w3l
w4 = w4l
else
w1 = w1s
w2 = w2s
w3 = w3s
w4 = w4s
endif
if(pdot(i).le.w4) then
tem = (pdot(i) - w4) / (w3 - w4)
elseif(pdot(i).ge.-w4) then
tem = - (pdot(i) + w4) / (w4 - w3)
else
tem = 0.
endif
val1 = -1.
tem = max(tem,val1)
val2 = 1.
tem = min(tem,val2)
tem = 1. - tem
tem1= .5*(cincrmax-cincrmin)
cincr = cincrmax - tem * tem1
pbcdif(i) = pfld(i,kb(i)) - pfld(i,kbcon(i))
if(pbcdif(i).gt.cincr) then
cnvflg(i) = .false.
endif
endif
enddo
!!
totflg = .true.
do i=1,im
totflg = totflg .and. (.not. cnvflg(i))
enddo
if(totflg) return
!!
c
c assume that updraft entrainment rate above cloud base is
c same as that at cloud base
c
do k = 1, km1
do i=1,im
if(cnvflg(i).and.
& (k.gt.kbcon(i).and.k.le.kmax(i))) then
xlamue(i,k) = xlamue(i,kbcon(i))
endif
enddo
enddo
c
c assume the detrainment rate for the updrafts to be same as
c the entrainment rate at cloud base
c
do i = 1, im
if(cnvflg(i)) then
xlamud(i) = xlamue(i,kbcon(i))
endif
enddo
c
c functions rapidly decreasing with height, mimicking a cloud ensemble
c (Bechtold et al., 2008)
c
do k = 1, km1
do i=1,im
if(cnvflg(i).and.
& (k.gt.kbcon(i).and.k.le.kmax(i))) then
tem = qeso(i,k)/qeso(i,kbcon(i))
fent1(i,k) = tem**2
fent2(i,k) = tem**3
endif
enddo
enddo
c
c final entrainment rate as the sum of turbulent part and organized entrainment
c depending on the environmental relative humidity
c (Bechtold et al., 2008)
c
do k = 1, km1
do i=1,im
if(cnvflg(i).and.
& (k.ge.kbcon(i).and.k.le.kmax(i))) then
tem = cxlamu * frh(i,k) * fent2(i,k)
xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
endif
enddo
enddo
c
c determine updraft mass flux for the subcloud layers
c
do k = km1, 1, -1
do i = 1, im
if (cnvflg(i)) then
if(k.lt.kbcon(i).and.k.ge.kb(i)) then
dz = zi(i,k+1) - zi(i,k)
ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
endif
endif
enddo
enddo
c
c compute updraft cloud properties
c
do i = 1, im
if(cnvflg(i)) then
indx = kb(i)
hcko(i,indx) = heo(i,indx)
ucko(i,indx) = uo(i,indx)
vcko(i,indx) = vo(i,indx)
pwavo(i) = 0.
endif
enddo
c
c cloud property is modified by the entrainment process
c
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k.gt.kb(i).and.k.lt.kmax(i)) then
dz = zi(i,k) - zi(i,k-1)
tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
tem1 = 0.5 * xlamud(i) * dz
factor = 1. + tem - tem1
ptem = 0.5 * tem + pgcon
ptem1= 0.5 * tem - pgcon
hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
& (heo(i,k)+heo(i,k-1)))/factor
ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k)
& +ptem1*uo(i,k-1))/factor
vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k)
& +ptem1*vo(i,k-1))/factor
dbyo(i,k) = hcko(i,k) - heso(i,k)
endif
endif
enddo
enddo
c
c taking account into convection inhibition due to existence of
c dry layers below cloud base
c
do i=1,im
flg(i) = cnvflg(i)
kbcon1(i) = kmax(i)
enddo
do k = 2, km1
do i=1,im
if (flg(i).and.k.lt.kmax(i)) then
if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
kbcon1(i) = k
flg(i) = .false.
endif
endif
enddo
enddo
do i=1,im
if(cnvflg(i)) then
if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
endif
enddo
do i=1,im
if(cnvflg(i)) then
tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
if(tem.gt.dthk) then
cnvflg(i) = .false.
endif
endif
enddo
!!
totflg = .true.
do i = 1, im
totflg = totflg .and. (.not. cnvflg(i))
enddo
if(totflg) return
!!
c
c determine cloud top as the level of zero buoyancy
c
do i = 1, im
flg(i) = cnvflg(i)
ktcon(i) = 1
enddo
do k = 2, km
do i = 1, im
if (flg(i).and.k .lt. kmax(i)) then
if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
ktcon(i) = k
flg(i) = .false.
endif
endif
enddo
enddo
c
do i = 1, im
if(cnvflg(i).and.(pfld(i,kbcon(i))-pfld(i,ktcon(i))).lt.cthk)
& cnvflg(i) = .false.
enddo
!!
totflg = .true.
do i = 1, im
totflg = totflg .and. (.not. cnvflg(i))
enddo
if(totflg) return
!!
c
c search for downdraft originating level above theta-e minimum
c
do i = 1, im
if(cnvflg(i)) then
hmin(i) = heo(i,kbcon1(i))
lmin(i) = kbmax(i)
jmin(i) = kbmax(i)
endif
enddo
do k = 2, km1
do i = 1, im
if (cnvflg(i) .and. k .le. kbmax(i)) then
if(k.gt.kbcon1(i).and.heo(i,k).lt.hmin(i)) then
lmin(i) = k + 1
hmin(i) = heo(i,k)
endif
endif
enddo
enddo
c
c make sure that jmin(i) is within the cloud
c
do i = 1, im
if(cnvflg(i)) then
jmin(i) = min(lmin(i),ktcon(i)-1)
jmin(i) = max(jmin(i),kbcon1(i)+1)
if(jmin(i).ge.ktcon(i)) cnvflg(i) = .false.
endif
enddo
c
c specify upper limit of mass flux at cloud base
c
do i = 1, im
if(cnvflg(i)) then
! xmbmax(i) = .1
!
k = kbcon(i)
dp = 1000. * del(i,k)
xmbmax(i) = dp / (g * dt2)
!
! tem = dp / (g * dt2)
! xmbmax(i) = min(tem, xmbmax(i))
endif
enddo
c
c compute mass flux within cloud layer
c
do k = 2, km1
do i = 1, im
if(cnvflg(i))then
if(k.gt.kbcon(i).and.k.le.ktcon(i)) then
dz = zi(i,k) - zi(i,k-1)
ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
endif
endif
enddo
enddo
c
c compute cloud moisture property and precipitation
c
do i = 1, im
if (cnvflg(i)) then
aa1(i) = 0.
qcko(i,kb(i)) = qo(i,kb(i))
! rhbar(i) = 0.
endif
enddo
do k = 1, km
do i = 1, im
if (cnvflg(i)) then
if(k.gt.kb(i).and.k.lt.ktcon(i)) then
dz = zi(i,k) - zi(i,k-1)
dz1 = zo(i,k+1) - zo(i,k)
gamma = el2orc * qeso(i,k) / (to(i,k)**2)
qrch = qeso(i,k)
& + gamma * dbyo(i,k) / (hvap * (1. + gamma))
cj
tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
tem1 = 0.5 * xlamud(i) * dz
factor = 1. + tem - tem1
qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
& (qo(i,k)+qo(i,k-1)))/factor
cj
dq = eta(i,k) * (qcko(i,k) - qrch)
c
! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
c
c check if there is excess moisture to release latent heat
c
if(dq.gt.0.) then
dp = 1000. * del(i,k)
etah = .5 * (eta(i,k) + eta(i,k-1))
qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
aa1(i) = aa1(i) - dz1 * g * qlk
qc = qlk + qrch
pwo(i,k) = etah * c0 * dz * qlk
dellal(i,k) = etah * c1 * dz * qlk * g / dp
qcko(i,k) = qc
pwavo(i) = pwavo(i) + pwo(i,k)
endif
endif
endif
enddo
enddo
c
! do i = 1, im
! if(cnvflg(i)) then
! indx = ktcon(i) - kb(i) - 1
! rhbar(i) = rhbar(i) / float(indx)
! endif
! enddo
c
c this section is ready for cloud water
c
if(ncloud.gt.0) then
c
c compute liquid and vapor separation at cloud top
c
do i = 1, im
if(cnvflg(i)) then
k = ktcon(i)
gamma = el2orc * qeso(i,k) / (to(i,k)**2)
qrch = qeso(i,k)
& + gamma * dbyo(i,k) / (hvap * (1. + gamma))
dq = qcko(i,k-1) - qrch
c
c check if there is excess moisture to release latent heat
c
if(dq.gt.0.) then
qlko_ktcon(i) = dq
qcko(i,k-1) = qrch
endif
endif
enddo
endif
c
c calculate cloud work function at t+dt
c
do k = 1, km
do i = 1, im
if (cnvflg(i)) then
if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
dz1 = zo(i,k+1) - zo(i,k)
gamma = el2orc * qeso(i,k) / (to(i,k)**2)
rfact = 1. + delta * cp * gamma
& * to(i,k) / hvap
aa1(i) = aa1(i) +
& dz1 * (g / (cp * to(i,k)))
& * dbyo(i,k) / (1. + gamma)
& * rfact
val = 0.
aa1(i)=aa1(i)+
& dz1 * g * delta *
& max(val,(qeso(i,k) - qo(i,k)))
endif
endif
enddo
enddo
do i = 1, im
if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
enddo
!!
totflg = .true.
do i=1,im
totflg = totflg .and. (.not. cnvflg(i))
enddo
if(totflg) return
!!
c
ccccc if(lat.eq.latd.and.lon.eq.lond.and.cnvflg(i)) then
ccccc print *, ' aa1(i) before dwndrft =', aa1(i)
ccccc endif
c
c------- downdraft calculations
c
c--- compute precipitation efficiency in terms of windshear
c
do i = 1, im
if(cnvflg(i)) then
vshear(i) = 0.
endif
enddo
do k = 2, km
do i = 1, im
if (cnvflg(i)) then
if(k.gt.kb(i).and.k.le.ktcon(i)) then
shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2
& + (vo(i,k)-vo(i,k-1)) ** 2)
vshear(i) = vshear(i) + shear
endif
endif
enddo
enddo
do i = 1, im
if(cnvflg(i)) then
vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
e1=1.591-.639*vshear(i)
& +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
edt(i)=1.-e1
val = .9
edt(i) = min(edt(i),val)
val = .0
edt(i) = max(edt(i),val)
edto(i)=edt(i)
edtx(i)=edt(i)
endif
enddo
c
c determine detrainment rate between 1 and kbcon
c
do i = 1, im
if(cnvflg(i)) then
sumx(i) = 0.
endif
enddo
do k = 1, km1
do i = 1, im
if(cnvflg(i).and.k.ge.1.and.k.lt.kbcon(i)) then
dz = zi(i,k+1) - zi(i,k)
sumx(i) = sumx(i) + dz
endif
enddo
enddo
do i = 1, im
beta = betas
if(slimsk(i).eq.1.) beta = betal
if(cnvflg(i)) then
dz = (sumx(i)+zi(i,1))/float(kbcon(i))
tem = 1./float(kbcon(i))
xlamd(i) = (1.-beta**tem)/dz
endif
enddo
c
c determine downdraft mass flux
c
do k = km1, 1, -1
do i = 1, im
if (cnvflg(i) .and. k .le. kmax(i)-1) then
if(k.lt.jmin(i).and.k.ge.kbcon(i)) then
dz = zi(i,k+1) - zi(i,k)
ptem = xlamdd - xlamde
etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
else if(k.lt.kbcon(i)) then
dz = zi(i,k+1) - zi(i,k)
ptem = xlamd(i) + xlamdd - xlamde
etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
endif
endif
enddo
enddo
c
c--- downdraft moisture properties
c
do i = 1, im
if(cnvflg(i)) then
jmn = jmin(i)
hcdo(i,jmn) = heo(i,jmn)
qcdo(i,jmn) = qo(i,jmn)
qrcdo(i,jmn)= qeso(i,jmn)
ucdo(i,jmn) = uo(i,jmn)
vcdo(i,jmn) = vo(i,jmn)
pwevo(i) = 0.
endif
enddo
cj
do k = km1, 1, -1
do i = 1, im
if (cnvflg(i) .and. k.lt.jmin(i)) then
dz = zi(i,k+1) - zi(i,k)
if(k.ge.kbcon(i)) then
tem = xlamde * dz
tem1 = 0.5 * xlamdd * dz
else
tem = xlamde * dz
tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
endif
factor = 1. + tem - tem1
ptem = 0.5 * tem - pgcon
ptem1= 0.5 * tem + pgcon
hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5*
& (heo(i,k)+heo(i,k+1)))/factor
ucdo(i,k) = ((1.-tem1)*ucdo(i,k+1)+ptem*uo(i,k+1)
& +ptem1*uo(i,k))/factor
vcdo(i,k) = ((1.-tem1)*vcdo(i,k+1)+ptem*vo(i,k+1)
& +ptem1*vo(i,k))/factor
dbyo(i,k) = hcdo(i,k) - heso(i,k)
endif
enddo
enddo
c
do k = km1, 1, -1
do i = 1, im
if (cnvflg(i).and.k.lt.jmin(i)) then
gamma = el2orc * qeso(i,k) / (to(i,k)**2)
qrcdo(i,k) = qeso(i,k)+
& (1./hvap)*(gamma/(1.+gamma))*dbyo(i,k)
! detad = etad(i,k+1) - etad(i,k)
cj
dz = zi(i,k+1) - zi(i,k)
if(k.ge.kbcon(i)) then
tem = xlamde * dz
tem1 = 0.5 * xlamdd * dz
else
tem = xlamde * dz
tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
endif
factor = 1. + tem - tem1
qcdo(i,k) = ((1.-tem1)*qcdo(i,k+1)+tem*0.5*
& (qo(i,k)+qo(i,k+1)))/factor
cj
! pwdo(i,k) = etad(i,k+1) * qcdo(i,k+1) -
! & etad(i,k) * qrcdo(i,k)
! pwdo(i,k) = pwdo(i,k) - detad *
! & .5 * (qrcdo(i,k) + qrcdo(i,k+1))
cj
pwdo(i,k) = etad(i,k+1) * (qcdo(i,k) - qrcdo(i,k))
qcdo(i,k) = qrcdo(i,k)
pwevo(i) = pwevo(i) + pwdo(i,k)
endif
enddo
enddo
c
c--- final downdraft strength dependent on precip
c--- efficiency (edt), normalized condensate (pwav), and
c--- evaporate (pwev)
c
do i = 1, im
edtmax = edtmaxl
if(slimsk(i).eq.0.) edtmax = edtmaxs
if(cnvflg(i)) then
if(pwevo(i).lt.0.) then
edto(i) = -edto(i) * pwavo(i) / pwevo(i)
edto(i) = min(edto(i),edtmax)
else
edto(i) = 0.
endif
endif
enddo
c
c--- downdraft cloudwork functions
c
do k = km1, 1, -1
do i = 1, im
if (cnvflg(i) .and. k .lt. jmin(i)) then
gamma = el2orc * qeso(i,k) / to(i,k)**2
dhh=hcdo(i,k)
dt=to(i,k)
dg=gamma
dh=heso(i,k)
dz=-1.*(zo(i,k+1)-zo(i,k))
aa1(i)=aa1(i)+edto(i)*dz*(g/(cp*dt))*((dhh-dh)/(1.+dg))
& *(1.+delta*cp*dg*dt/hvap)
val=0.
aa1(i)=aa1(i)+edto(i)*
& dz*g*delta*max(val,(qeso(i,k)-qo(i,k)))
endif
enddo
enddo
do i = 1, im
if(cnvflg(i).and.aa1(i).le.0.) then
cnvflg(i) = .false.
endif
enddo
!!
totflg = .true.
do i=1,im
totflg = totflg .and. (.not. cnvflg(i))
enddo
if(totflg) return
!!
c
c--- what would the change be, that a cloud with unit mass
c--- will do to the environment?
c
do k = 1, km
do i = 1, im
if(cnvflg(i) .and. k .le. kmax(i)) then
dellah(i,k) = 0.
dellaq(i,k) = 0.
dellau(i,k) = 0.
dellav(i,k) = 0.
endif
enddo
enddo
do i = 1, im
if(cnvflg(i)) then
dp = 1000. * del(i,1)
dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1)
& - heo(i,1)) * g / dp
dellaq(i,1) = edto(i) * etad(i,1) * (qcdo(i,1)
& - qo(i,1)) * g / dp
dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1)
& - uo(i,1)) * g / dp
dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1)
& - vo(i,1)) * g / dp
endif
enddo
c
c--- changed due to subsidence and entrainment
c
do k = 2, km1
do i = 1, im
if (cnvflg(i).and.k.lt.ktcon(i)) then
aup = 1.
if(k.le.kb(i)) aup = 0.
adw = 1.
if(k.gt.jmin(i)) adw = 0.
dp = 1000. * del(i,k)
dz = zi(i,k) - zi(i,k-1)
c
dv1h = heo(i,k)
dv2h = .5 * (heo(i,k) + heo(i,k-1))
dv3h = heo(i,k-1)
dv1q = qo(i,k)
dv2q = .5 * (qo(i,k) + qo(i,k-1))
dv3q = qo(i,k-1)
dv1u = uo(i,k)
dv2u = .5 * (uo(i,k) + uo(i,k-1))
dv3u = uo(i,k-1)
dv1v = vo(i,k)
dv2v = .5 * (vo(i,k) + vo(i,k-1))
dv3v = vo(i,k-1)
c
tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
tem1 = xlamud(i)
c
if(k.le.kbcon(i)) then
ptem = xlamde
ptem1 = xlamd(i)+xlamdd
else
ptem = xlamde
ptem1 = xlamdd
endif
cj
dellah(i,k) = dellah(i,k) +
& ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1h
& - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3h
& - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2h*dz
& + aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz
& + adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1))*dz
& ) *g/dp
cj
dellaq(i,k) = dellaq(i,k) +
& ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1q
& - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3q
& - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz
& + aup*tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz
& + adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qrcdo(i,k-1))*dz
& ) *g/dp
cj
dellau(i,k) = dellau(i,k) +
& ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1u
& - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3u
& - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2u*dz
& + aup*tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz
& + adw*edto(i)*ptem1*etad(i,k)*.5*(ucdo(i,k)+ucdo(i,k-1))*dz
& - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1u-dv3u)
& ) *g/dp
cj
dellav(i,k) = dellav(i,k) +
& ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1v
& - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3v
& - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2v*dz
& + aup*tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz
& + adw*edto(i)*ptem1*etad(i,k)*.5*(vcdo(i,k)+vcdo(i,k-1))*dz
& - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1v-dv3v)
& ) *g/dp
cj
endif
enddo
enddo
c
c------- cloud top
c
do i = 1, im
if(cnvflg(i)) then
indx = ktcon(i)
dp = 1000. * del(i,indx)
dv1h = heo(i,indx-1)
dellah(i,indx) = eta(i,indx-1) *
& (hcko(i,indx-1) - dv1h) * g / dp
dv1q = qo(i,indx-1)
dellaq(i,indx) = eta(i,indx-1) *
& (qcko(i,indx-1) - dv1q) * g / dp
dv1u = uo(i,indx-1)
dellau(i,indx) = eta(i,indx-1) *
& (ucko(i,indx-1) - dv1u) * g / dp
dv1v = vo(i,indx-1)
dellav(i,indx) = eta(i,indx-1) *
& (vcko(i,indx-1) - dv1v) * g / dp
c
c cloud water
c
dellal(i,indx) = eta(i,indx-1) *
& qlko_ktcon(i) * g / dp
endif
enddo
c
c------- final changed variable per unit mass flux
c
do k = 1, km
do i = 1, im
if (cnvflg(i).and.k .le. kmax(i)) then
if(k.gt.ktcon(i)) then
qo(i,k) = q1(i,k)
to(i,k) = t1(i,k)
endif
if(k.le.ktcon(i)) then
qo(i,k) = dellaq(i,k) * mbdt + q1(i,k)
dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
to(i,k) = dellat * mbdt + t1(i,k)
val = 1.e-10
qo(i,k) = max(qo(i,k), val )
endif
endif
enddo
enddo
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
c--- the above changed environment is now used to calulate the
c--- effect the arbitrary cloud (with unit mass flux)
c--- would have on the stability,
c--- which then is used to calculate the real mass flux,
c--- necessary to keep this change in balance with the large-scale
c--- destabilization.
c
c--- environmental conditions again, first heights
c
do k = 1, km
do i = 1, im
if(cnvflg(i) .and. k .le. kmax(i)) then
qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
qeso(i,k) = eps * qeso(i,k) / (pfld(i,k)+epsm1*qeso(i,k))
val = 1.e-8
qeso(i,k) = max(qeso(i,k), val )
! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
endif
enddo
enddo
c
c--- moist static energy
c
do k = 1, km1
do i = 1, im
if(cnvflg(i) .and. k .le. kmax(i)-1) then
dz = .5 * (zo(i,k+1) - zo(i,k))
dp = .5 * (pfld(i,k+1) - pfld(i,k))
es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
pprime = pfld(i,k+1) + epsm1 * es
qs = eps * es / pprime
dqsdp = - qs / pprime
desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
dq = dqsdt * dt + dqsdp * dp
to(i,k) = to(i,k+1) + dt
qo(i,k) = qo(i,k+1) + dq
po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
endif
enddo
enddo
do k = 1, km1
do i = 1, im
if(cnvflg(i) .and. k .le. kmax(i)-1) then
qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1 * qeso(i,k))
val1 = 1.e-8
qeso(i,k) = max(qeso(i,k), val1)
val2 = 1.e-10
qo(i,k) = max(qo(i,k), val2 )
! qo(i,k) = min(qo(i,k),qeso(i,k))
heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
& cp * to(i,k) + hvap * qo(i,k)
heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) +
& cp * to(i,k) + hvap * qeso(i,k)
endif
enddo
enddo
do i = 1, im
if(cnvflg(i)) then
k = kmax(i)
heo(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qo(i,k)
heso(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qeso(i,k)
c heo(i,k) = min(heo(i,k),heso(i,k))
endif
enddo
c
c**************************** static control
c
c------- moisture and cloud work functions
c
do i = 1, im
if(cnvflg(i)) then
xaa0(i) = 0.
xpwav(i) = 0.
endif
enddo
c
do i = 1, im
if(cnvflg(i)) then
indx = kb(i)
hcko(i,indx) = heo(i,indx)
qcko(i,indx) = qo(i,indx)
endif
enddo
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k.gt.kb(i).and.k.le.ktcon(i)) then
dz = zi(i,k) - zi(i,k-1)
tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
tem1 = 0.5 * xlamud(i) * dz
factor = 1. + tem - tem1
hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*
& (heo(i,k)+heo(i,k-1)))/factor
endif
endif
enddo
enddo
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k.gt.kb(i).and.k.lt.ktcon(i)) then
dz = zi(i,k) - zi(i,k-1)
gamma = el2orc * qeso(i,k) / (to(i,k)**2)
xdby = hcko(i,k) - heso(i,k)
val = 0.
xdby = max(xdby,val)
xqrch = qeso(i,k)
& + gamma * xdby / (hvap * (1. + gamma))
cj
tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
tem1 = 0.5 * xlamud(i) * dz
factor = 1. + tem - tem1
qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*
& (qo(i,k)+qo(i,k-1)))/factor
cj
dq = eta(i,k) * (qcko(i,k) - xqrch)
c
if(dq.gt.0.) then
etah = .5 * (eta(i,k) + eta(i,k-1))
qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
xaa0(i) = xaa0(i) - (zo(i,k) - zo(i,k-1)) * g * qlk
xqc = qlk + xqrch
xpw = etah * c0 * dz * qlk
qcko(i,k) = xqc
xpwav(i) = xpwav(i) + xpw
endif
endif
if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
dz1 = zo(i,k+1) - zo(i,k)
gamma = el2orc * qeso(i,k) / (to(i,k)**2)
rfact = 1. + delta * cp * gamma
& * to(i,k) / hvap
xdby = hcko(i,k) - heso(i,k)
xaa0(i) = xaa0(i)
& + dz1 * (g / (cp * to(i,k)))
& * xdby / (1. + gamma)
& * rfact
val=0.
xaa0(i)=xaa0(i)+
& dz1 * g * delta *
& max(val,(qeso(i,k) - qo(i,k)))
endif
endif
enddo
enddo
c
c------- downdraft calculations
c
c--- downdraft moisture properties
c
do i = 1, im
if(cnvflg(i)) then
jmn = jmin(i)
hcdo(i,jmn) = heo(i,jmn)
qcdo(i,jmn) = qo(i,jmn)
qrcd(i,jmn) = qeso(i,jmn)
xpwev(i) = 0.
endif
enddo
cj
do k = km1, 1, -1
do i = 1, im
if (cnvflg(i) .and. k.lt.jmin(i)) then
dz = zi(i,k+1) - zi(i,k)
if(k.ge.kbcon(i)) then
tem = xlamde * dz
tem1 = 0.5 * xlamdd * dz
else
tem = xlamde * dz
tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
endif
factor = 1. + tem - tem1
hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5*
& (heo(i,k)+heo(i,k+1)))/factor
endif
enddo
enddo
cj
do k = km1, 1, -1
do i = 1, im
if (cnvflg(i) .and. k .lt. jmin(i)) then
dq = qeso(i,k)
dt = to(i,k)
gamma = el2orc * dq / dt**2
dh = hcdo(i,k) - heso(i,k)
qrcd(i,k)=dq+(1./hvap)*(gamma/(1.+gamma))*dh
! detad = etad(i,k+1) - etad(i,k)
cj
dz = zi(i,k+1) - zi(i,k)
if(k.ge.kbcon(i)) then
tem = xlamde * dz
tem1 = 0.5 * xlamdd * dz
else
tem = xlamde * dz
tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
endif
factor = 1. + tem - tem1
qcdo(i,k) = ((1.-tem1)*qcdo(i,k+1)+tem*0.5*
& (qo(i,k)+qo(i,k+1)))/factor
cj
! xpwd = etad(i,k+1) * qcdo(i,k+1) -
! & etad(i,k) * qrcd(i,k)
! xpwd = xpwd - detad *
! & .5 * (qrcd(i,k) + qrcd(i,k+1))
cj
xpwd = etad(i,k+1) * (qcdo(i,k) - qrcd(i,k))
qcdo(i,k)= qrcd(i,k)
xpwev(i) = xpwev(i) + xpwd
endif
enddo
enddo
c
do i = 1, im
edtmax = edtmaxl
if(slimsk(i).eq.0.) edtmax = edtmaxs
if(cnvflg(i)) then
if(xpwev(i).ge.0.) then
edtx(i) = 0.
else
edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
edtx(i) = min(edtx(i),edtmax)
endif
endif
enddo
c
c
c--- downdraft cloudwork functions
c
c
do k = km1, 1, -1
do i = 1, im
if (cnvflg(i) .and. k.lt.jmin(i)) then
gamma = el2orc * qeso(i,k) / to(i,k)**2
dhh=hcdo(i,k)
dt= to(i,k)
dg= gamma
dh= heso(i,k)
dz=-1.*(zo(i,k+1)-zo(i,k))
xaa0(i)=xaa0(i)+edtx(i)*dz*(g/(cp*dt))*((dhh-dh)/(1.+dg))
& *(1.+delta*cp*dg*dt/hvap)
val=0.
xaa0(i)=xaa0(i)+edtx(i)*
& dz*g*delta*max(val,(qeso(i,k)-qo(i,k)))
endif
enddo
enddo
c
c calculate critical cloud work function
c
do i = 1, im
if(cnvflg(i)) then
if(pfld(i,ktcon(i)).lt.pcrit(15))then
acrt(i)=acrit(15)*(975.-pfld(i,ktcon(i)))
& /(975.-pcrit(15))
else if(pfld(i,ktcon(i)).gt.pcrit(1))then
acrt(i)=acrit(1)
else
k = int((850. - pfld(i,ktcon(i)))/50.) + 2
k = min(k,15)
k = max(k,2)
acrt(i)=acrit(k)+(acrit(k-1)-acrit(k))*
& (pfld(i,ktcon(i))-pcrit(k))/(pcrit(k-1)-pcrit(k))
endif
endif
enddo
do i = 1, im
if(cnvflg(i)) then
if(slimsk(i).eq.1.) then
w1 = w1l
w2 = w2l
w3 = w3l
w4 = w4l
else
w1 = w1s
w2 = w2s
w3 = w3s
w4 = w4s
endif
c
c modify critical cloud workfunction by cloud base vertical velocity
c
if(pdot(i).le.w4) then
acrtfct(i) = (pdot(i) - w4) / (w3 - w4)
elseif(pdot(i).ge.-w4) then
acrtfct(i) = - (pdot(i) + w4) / (w4 - w3)
else
acrtfct(i) = 0.
endif
val1 = -1.
acrtfct(i) = max(acrtfct(i),val1)
val2 = 1.
acrtfct(i) = min(acrtfct(i),val2)
acrtfct(i) = 1. - acrtfct(i)
c
c modify acrtfct(i) by colume mean rh if rhbar(i) is greater than 80 percent
c
c if(rhbar(i).ge..8) then
c acrtfct(i) = acrtfct(i) * (.9 - min(rhbar(i),.9)) * 10.
c endif
c
c modify adjustment time scale by cloud base vertical velocity
c
dtconv(i) = dt2 + max((1800. - dt2),0.) *
& (pdot(i) - w2) / (w1 - w2)
c dtconv(i) = max(dtconv(i), dt2)
c dtconv(i) = 1800. * (pdot(i) - w2) / (w1 - w2)
dtconv(i) = max(dtconv(i),dtmin)
dtconv(i) = min(dtconv(i),dtmax)
c
endif
enddo
c
c--- large scale forcing
c
do i= 1, im
if(cnvflg(i)) then
fld(i)=(aa1(i)-acrt(i)* acrtfct(i))/dtconv(i)
if(fld(i).le.0.) cnvflg(i) = .false.
endif
if(cnvflg(i)) then
c xaa0(i) = max(xaa0(i),0.)
xk(i) = (xaa0(i) - aa1(i)) / mbdt
if(xk(i).ge.0.) cnvflg(i) = .false.
endif
c
c--- kernel, cloud base mass flux
c
if(cnvflg(i)) then
xmb(i) = -fld(i) / xk(i)
xmb(i) = min(xmb(i),xmbmax(i))
endif
enddo
!!
totflg = .true.
do i=1,im
totflg = totflg .and. (.not. cnvflg(i))
enddo
if(totflg) return
!!
c
c restore to,qo,uo,vo to t1,q1,u1,v1 in case convection stops
c
do k = 1, km
do i = 1, im
if (cnvflg(i) .and. k .le. kmax(i)) then
to(i,k) = t1(i,k)
qo(i,k) = q1(i,k)
uo(i,k) = u1(i,k)
vo(i,k) = v1(i,k)
qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
val = 1.e-8
qeso(i,k) = max(qeso(i,k), val )
endif
enddo
enddo
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
c--- feedback: simply the changes from the cloud with unit mass flux
c--- multiplied by the mass flux necessary to keep the
c--- equilibrium with the larger-scale.
c
do i = 1, im
delhbar(i) = 0.
delqbar(i) = 0.
deltbar(i) = 0.
delubar(i) = 0.
delvbar(i) = 0.
qcond(i) = 0.
enddo
do k = 1, km
do i = 1, im
if (cnvflg(i) .and. k .le. kmax(i)) then
if(k.le.ktcon(i)) then
dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
tem = 1./rcs(i)
u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
dp = 1000. * del(i,k)
delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
endif
endif
enddo
enddo
do k = 1, km
do i = 1, im
if (cnvflg(i) .and. k .le. kmax(i)) then
if(k.le.ktcon(i)) then
qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
val = 1.e-8
qeso(i,k) = max(qeso(i,k), val )
endif
endif
enddo
enddo
c
do i = 1, im
rntot(i) = 0.
delqev(i) = 0.
delq2(i) = 0.
flg(i) = cnvflg(i)
enddo
do k = km, 1, -1
do i = 1, im
if (cnvflg(i) .and. k .le. kmax(i)) then
if(k.lt.ktcon(i)) then
aup = 1.
if(k.le.kb(i)) aup = 0.
adw = 1.
if(k.ge.jmin(i)) adw = 0.
rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
rntot(i) = rntot(i) + rain * xmb(i) * .001 * dt2
endif
endif
enddo
enddo
do k = km, 1, -1
do i = 1, im
if (k .le. kmax(i)) then
deltv(i) = 0.
delq(i) = 0.
qevap(i) = 0.
if(cnvflg(i).and.k.lt.ktcon(i)) then
aup = 1.
if(k.le.kb(i)) aup = 0.
adw = 1.
if(k.ge.jmin(i)) adw = 0.
rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
rn(i) = rn(i) + rain * xmb(i) * .001 * dt2
endif
if(flg(i).and.k.lt.ktcon(i)) then
evef = edt(i) * evfact
if(slimsk(i).eq.1.) evef=edt(i) * evfactl
! if(slimsk(i).eq.1.) evef=.07
c if(slimsk(i).ne.1.) evef = 0.
qcond(i) = evef * (q1(i,k) - qeso(i,k))
& / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
dp = 1000. * del(i,k)
if(rn(i).gt.0..and.qcond(i).lt.0.) then
qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
endif
if(rn(i).gt.0..and.qcond(i).lt.0..and.
& delq2(i).gt.rntot(i)) then
qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
flg(i) = .false.
endif
if(rn(i).gt.0..and.qevap(i).gt.0.) then
q1(i,k) = q1(i,k) + qevap(i)
t1(i,k) = t1(i,k) - elocp * qevap(i)
rn(i) = rn(i) - .001 * qevap(i) * dp / g
deltv(i) = - elocp*qevap(i)/dt2
delq(i) = + qevap(i)/dt2
delqev(i) = delqev(i) + .001*dp*qevap(i)/g
endif
dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
delqbar(i) = delqbar(i) + delq(i)*dp/g
deltbar(i) = deltbar(i) + deltv(i)*dp/g
endif
endif
enddo
enddo
cj
! do i = 1, im
! if(me.eq.31.and.cnvflg(i)) then
! if(cnvflg(i)) then
! print *, ' deep delhbar, delqbar, deltbar = ',
! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
! print *, ' deep delubar, delvbar = ',delubar(i),delvbar(i)
! print *, ' precip =', hvap*rn(i)*1000./dt2
! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
! endif
! enddo
c
c precipitation rate converted to actual precip
c in unit of m instead of kg
c
do i = 1, im
if(cnvflg(i)) then
c
c in the event of upper level rain evaporation and lower level downdraft
c moistening, rn can become negative, in this case, we back out of the
c heating and the moistening
c
if(rn(i).lt.0..and..not.flg(i)) rn(i) = 0.
if(rn(i).le.0.) then
rn(i) = 0.
else
ktop(i) = ktcon(i)
kbot(i) = kbcon(i)
kcnv(i) = 1
cldwrk(i) = aa1(i)
endif
endif
enddo
c
c cloud water
c
if (ncloud.gt.0) then
!
do k = 1, km
do i = 1, im
if (cnvflg(i) .and. rn(i).gt.0.) then
if (k.gt.kb(i).and.k.le.ktcon(i)) then
tem = dellal(i,k) * xmb(i) * dt2
tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
if (ql(i,k,2) .gt. -999.0) then
ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice
ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water
else
ql(i,k,1) = ql(i,k,1) + tem
endif
endif
endif
enddo
enddo
!
endif
c
do k = 1, km
do i = 1, im
if(cnvflg(i).and.rn(i).le.0.) then
if (k .le. kmax(i)) then
t1(i,k) = to(i,k)
q1(i,k) = qo(i,k)
u1(i,k) = uo(i,k)
v1(i,k) = vo(i,k)
endif
endif
enddo
enddo
!
! hchuang code change
!
do k = 1, km
do i = 1, im
if(cnvflg(i).and.rn(i).gt.0.) then
if(k.ge.kb(i) .and. k.lt.ktop(i)) then
ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
endif
endif
enddo
enddo
do i = 1, im
if(cnvflg(i).and.rn(i).gt.0.) then
k = ktop(i)-1
dt_mf(i,k) = ud_mf(i,k)
endif
enddo
do k = 1, km
do i = 1, im
if(cnvflg(i).and.rn(i).gt.0.) then
if(k.ge.1 .and. k.le.jmin(i)) then
dd_mf(i,k) = edto(i) * etad(i,k) * xmb(i) * dt2
endif
endif
enddo
enddo
!!
return
end