MODULE module_mp_wsm6
USE module_utility, ONLY: WRFU_Clock, WRFU_Alarm
USE module_domain, ONLY : HISTORY_ALARM, Is_alarm_tstep
USE module_mp_radar
REAL, PARAMETER, PRIVATE :: dtcldcr = 120.
REAL, PARAMETER, PRIVATE :: n0r = 8.e6
REAL, PARAMETER, PRIVATE :: avtr = 841.9
REAL, PARAMETER, PRIVATE :: bvtr = 0.8
REAL, PARAMETER, PRIVATE :: r0 = .8e-5
REAL, PARAMETER, PRIVATE :: peaut = .55
REAL, PARAMETER, PRIVATE :: xncr = 3.e8
REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5
REAL, PARAMETER, PRIVATE :: avts = 11.72
REAL, PARAMETER, PRIVATE :: bvts = .41
REAL, PARAMETER, PRIVATE :: n0smax = 1.e11
REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4
REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5
REAL, PARAMETER, PRIVATE :: dicon = 11.9
REAL, PARAMETER, PRIVATE :: dimax = 500.e-6
REAL, PARAMETER, PRIVATE :: n0s = 2.e6
REAL, PARAMETER, PRIVATE :: alpha = .12
REAL, PARAMETER, PRIVATE :: pfrz1 = 100.
REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66
REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9
REAL, PARAMETER, PRIVATE :: eacrc = 1.0
REAL, PARAMETER, PRIVATE :: dens = 100.0
REAL, PARAMETER, PRIVATE :: qs0 = 6.e-4
REAL, SAVE :: &
qc0, qck1, pidnc, &
bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
bvtr6,g6pbr, &
precr1,precr2,roqimax,bvts1, &
bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
n0g,avtg,bvtg,deng,lamdagmax, &
g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
pidn0s,xlv1,pacrc,pi, &
bvtg1,bvtg2,bvtg3,bvtg4,g1pbg, &
g3pbg,g4pbg,g5pbgo2,pvtg,pacrg, &
precg1,precg2,pidn0g, &
rslopermax,rslopesmax,rslopegmax, &
rsloperbmax,rslopesbmax,rslopegbmax, &
rsloper2max,rslopes2max,rslopeg2max, &
rsloper3max,rslopes3max,rslopeg3max
CONTAINS
SUBROUTINE wsm6(th, q, qc, qr, qi, qs, qg &
,den, pii, p, delz &
,delt,g, cpd, cpv, rd, rv, t0c &
,ep1, ep2, qmin &
,XLS, XLV0, XLF0, den0, denr &
,cliq,cice,psat &
,rain, rainncv &
,snow, snowncv &
,sr &
,refl_10cm, diagflag, do_radar_ref &
,graupel, graupelncv &
,has_reqc, has_reqi, has_reqs &
,re_cloud, re_ice, re_snow &
,ids,ide, jds,jde, kds,kde &
,ims,ime, jms,jme, kms,kme &
,its,ite, jts,jte, kts,kte &
)
IMPLICIT NONE
INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
ims,ime, jms,jme, kms,kme , &
its,ite, jts,jte, kts,kte
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(INOUT) :: &
th, &
q, &
qc, &
qi, &
qr, &
qs, &
qg
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(IN ) :: &
den, &
pii, &
p, &
delz
REAL, INTENT(IN ) :: delt, &
g, &
rd, &
rv, &
t0c, &
den0, &
cpd, &
cpv, &
ep1, &
ep2, &
qmin, &
XLS, &
XLV0, &
XLF0, &
cliq, &
cice, &
psat, &
denr
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(INOUT) :: rain, &
rainncv, &
sr
INTEGER, INTENT(IN):: &
has_reqc, &
has_reqi, &
has_reqs
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), &
INTENT(INOUT):: &
re_cloud, &
re_ice, &
re_snow
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
refl_10cm
REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
INTENT(INOUT) :: snow, &
snowncv
REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
INTENT(INOUT) :: graupel, &
graupelncv
REAL, DIMENSION( its:ite , kts:kte ) :: t
REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci
REAL, DIMENSION( its:ite , kts:kte, 3 ) :: qrs
INTEGER :: i,j,k
REAL, DIMENSION(kts:kte):: qv1d, t1d, p1d, qr1d, qs1d, qg1d, dBZ
LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
REAL, DIMENSION( kts:kte ) :: den1d
REAL, DIMENSION( kts:kte ) :: qc1d
REAL, DIMENSION( kts:kte ) :: qi1d
REAL, DIMENSION( kts:kte ) :: re_qc, re_qi, re_qs
DO j=jts,jte
DO k=kts,kte
DO i=its,ite
t(i,k)=th(i,k,j)*pii(i,k,j)
qci(i,k,1) = qc(i,k,j)
qci(i,k,2) = qi(i,k,j)
qrs(i,k,1) = qr(i,k,j)
qrs(i,k,2) = qs(i,k,j)
qrs(i,k,3) = qg(i,k,j)
ENDDO
ENDDO
CALL wsm62D(t, q(ims,kms,j), qci, qrs &
,den(ims,kms,j) &
,p(ims,kms,j), delz(ims,kms,j) &
,delt,g, cpd, cpv, rd, rv, t0c &
,ep1, ep2, qmin &
,XLS, XLV0, XLF0, den0, denr &
,cliq,cice,psat &
,j &
,rain(ims,j),rainncv(ims,j) &
,sr(ims,j) &
,ids,ide, jds,jde, kds,kde &
,ims,ime, jms,jme, kms,kme &
,its,ite, jts,jte, kts,kte &
,snow,snowncv &
,graupel,graupelncv &
)
DO K=kts,kte
DO I=its,ite
th(i,k,j)=t(i,k)/pii(i,k,j)
qc(i,k,j) = qci(i,k,1)
qi(i,k,j) = qci(i,k,2)
qr(i,k,j) = qrs(i,k,1)
qs(i,k,j) = qrs(i,k,2)
qg(i,k,j) = qrs(i,k,3)
ENDDO
ENDDO
IF ( PRESENT (diagflag) ) THEN
if (diagflag .and. do_radar_ref == 1) then
DO I=its,ite
DO K=kts,kte
t1d(k)=th(i,k,j)*pii(i,k,j)
p1d(k)=p(i,k,j)
qv1d(k)=q(i,k,j)
qr1d(k)=qr(i,k,j)
qs1d(k)=qs(i,k,j)
qg1d(k)=qg(i,k,j)
ENDDO
call refl10cm_wsm6 (qv1d, qr1d, qs1d, qg1d, &
t1d, p1d, dBZ, kts, kte, i, j)
do k = kts, kte
refl_10cm(i,k,j) = MAX(-35., dBZ(k))
enddo
ENDDO
endif
ENDIF
if (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) then
do i=its,ite
do k=kts,kte
re_qc(k) = 2.51E-6
re_qi(k) = 10.01E-6
re_qs(k) = 25.E-6
t1d(k) = th(i,k,j)*pii(i,k,j)
den1d(k)= den(i,k,j)
qc1d(k) = qc(i,k,j)
qi1d(k) = qi(i,k,j)
qs1d(k) = qs(i,k,j)
enddo
call effectRad_wsm6(t1d, qc1d, qi1d, qs1d, den1d, &
qmin, t0c, re_qc, re_qi, re_qs, &
kts, kte, i, j)
do k=kts,kte
re_cloud(i,k,j) = MAX(2.51E-6, MIN(re_qc(k), 50.E-6))
re_ice(i,k,j) = MAX(10.01E-6, MIN(re_qi(k), 125.E-6))
re_snow(i,k,j) = MAX(25.E-6, MIN(re_qs(k), 999.E-6))
enddo
enddo
endif
ENDDO
END SUBROUTINE wsm6
SUBROUTINE wsm62D(t, q &
,qci, qrs, den, p, delz &
,delt,g, cpd, cpv, rd, rv, t0c &
,ep1, ep2, qmin &
,XLS, XLV0, XLF0, den0, denr &
,cliq,cice,psat &
,lat &
,rain,rainncv &
,sr &
,ids,ide, jds,jde, kds,kde &
,ims,ime, jms,jme, kms,kme &
,its,ite, jts,jte, kts,kte &
,snow,snowncv &
,graupel,graupelncv &
)
IMPLICIT NONE
INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
ims,ime, jms,jme, kms,kme , &
its,ite, jts,jte, kts,kte, &
lat
REAL, DIMENSION( its:ite , kts:kte ), &
INTENT(INOUT) :: &
t
REAL, DIMENSION( its:ite , kts:kte, 2 ), &
INTENT(INOUT) :: &
qci
REAL, DIMENSION( its:ite , kts:kte, 3 ), &
INTENT(INOUT) :: &
qrs
REAL, DIMENSION( ims:ime , kms:kme ), &
INTENT(INOUT) :: &
q
REAL, DIMENSION( ims:ime , kms:kme ), &
INTENT(IN ) :: &
den, &
p, &
delz
REAL, INTENT(IN ) :: delt, &
g, &
cpd, &
cpv, &
t0c, &
den0, &
rd, &
rv, &
ep1, &
ep2, &
qmin, &
XLS, &
XLV0, &
XLF0, &
cliq, &
cice, &
psat, &
denr
REAL, DIMENSION( ims:ime ), &
INTENT(INOUT) :: rain, &
rainncv, &
sr
REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
INTENT(INOUT) :: snow, &
snowncv
REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
INTENT(INOUT) :: graupel, &
graupelncv
REAL, DIMENSION( its:ite , kts:kte , 3) :: &
rh, &
qs, &
rslope, &
rslope2, &
rslope3, &
rslopeb, &
qrs_tmp, &
falk, &
fall, &
work1
REAL, DIMENSION( its:ite , kts:kte ) :: &
fallc, &
falkc, &
work1c, &
work2c, &
workr, &
worka
REAL, DIMENSION( its:ite , kts:kte ) :: &
den_tmp, &
delz_tmp
REAL, DIMENSION( its:ite , kts:kte ) :: &
pigen, &
pidep, &
pcond, &
prevp, &
psevp, &
pgevp, &
psdep, &
pgdep, &
praut, &
psaut, &
pgaut, &
piacr, &
pracw, &
praci, &
pracs, &
psacw, &
psaci, &
psacr, &
pgacw, &
pgaci, &
pgacr, &
pgacs, &
paacw, &
psmlt, &
pgmlt, &
pseml, &
pgeml
REAL, DIMENSION( its:ite , kts:kte ) :: &
qsum, &
xl, &
cpm, &
work2, &
denfac, &
xni, &
denqrs1, &
denqrs2, &
denqrs3, &
denqci, &
n0sfac
REAL, DIMENSION( its:ite ) :: delqrs1, &
delqrs2, &
delqrs3, &
delqi
REAL, DIMENSION( its:ite ) :: tstepsnow, &
tstepgraup
INTEGER, DIMENSION( its:ite ) :: mstep, &
numdt
LOGICAL, DIMENSION( its:ite ) :: flgcld
REAL :: &
cpmcal, xlcal, diffus, &
viscos, xka, venfac, conden, diffac, &
x, y, z, a, b, c, d, e, &
qdt, holdrr, holdrs, holdrg, supcol, supcolt, pvt, &
coeres, supsat, dtcld, xmi, eacrs, satdt, &
qimax, diameter, xni0, roqi0, &
fallsum, fallsum_qsi, fallsum_qg, &
vt2i,vt2r,vt2s,vt2g,acrfac,egs,egi, &
xlwork2, factor, source, value, &
xlf, pfrzdtc, pfrzdtr, supice, alpha2, delta2, delta3
REAL :: vt2ave
REAL :: holdc, holdci
INTEGER :: i, j, k, mstepmax, &
iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
REAL, DIMENSION( its:ite ) :: tvec1
REAL :: temp
cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
xlcal(x) = xlv0-xlv1*(x-t0c)
diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y
viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y
xka(x,y) = 1.414e3*viscos(x,y)*y
diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
/sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
idim = ite-its+1
kdim = kte-kts+1
do k = kts, kte
do i = its, ite
qci(i,k,1) = max(qci(i,k,1),0.0)
qrs(i,k,1) = max(qrs(i,k,1),0.0)
qci(i,k,2) = max(qci(i,k,2),0.0)
qrs(i,k,2) = max(qrs(i,k,2),0.0)
qrs(i,k,3) = max(qrs(i,k,3),0.0)
enddo
enddo
do k = kts, kte
do i = its, ite
cpm(i,k) = cpmcal(q(i,k))
xl(i,k) = xlcal(t(i,k))
enddo
enddo
do k = kts, kte
do i = its, ite
delz_tmp(i,k) = delz(i,k)
den_tmp(i,k) = den(i,k)
enddo
enddo
do i = its, ite
rainncv(i) = 0.
if(PRESENT (snowncv) .AND. PRESENT (snow)) snowncv(i,lat) = 0.
if(PRESENT (graupelncv) .AND. PRESENT (graupel)) graupelncv(i,lat) = 0.
sr(i) = 0.
tstepsnow(i) = 0.
tstepgraup(i) = 0.
enddo
loops = max(nint(delt/dtcldcr),1)
dtcld = delt/loops
if(delt.le.dtcldcr) dtcld = delt
do loop = 1,loops
do i = its, ite
mstep(i) = 1
flgcld(i) = .true.
enddo
do k = kts, kte
CALL vsrec( tvec1(its), den(its,k), ite-its+1)
do i = its, ite
tvec1(i) = tvec1(i)*den0
enddo
CALL vssqrt( denfac(its,k), tvec1(its), ite-its+1)
enddo
hsub = xls
hvap = xlv0
cvap = cpv
ttp=t0c+0.01
dldt=cvap-cliq
xa=-dldt/rv
xb=xa+hvap/(rv*ttp)
dldti=cvap-cice
xai=-dldti/rv
xbi=xai+hsub/(rv*ttp)
do k = kts, kte
do i = its, ite
tr=ttp/t(i,k)
qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
qs(i,k,1) = max(qs(i,k,1),qmin)
rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
tr=ttp/t(i,k)
if(t(i,k).lt.ttp) then
qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
else
qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
endif
qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
qs(i,k,2) = max(qs(i,k,2),qmin)
rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
enddo
enddo
do k = kts, kte
do i = its, ite
prevp(i,k) = 0.
psdep(i,k) = 0.
pgdep(i,k) = 0.
praut(i,k) = 0.
psaut(i,k) = 0.
pgaut(i,k) = 0.
pracw(i,k) = 0.
praci(i,k) = 0.
piacr(i,k) = 0.
psaci(i,k) = 0.
psacw(i,k) = 0.
pracs(i,k) = 0.
psacr(i,k) = 0.
pgacw(i,k) = 0.
paacw(i,k) = 0.
pgaci(i,k) = 0.
pgacr(i,k) = 0.
pgacs(i,k) = 0.
pigen(i,k) = 0.
pidep(i,k) = 0.
pcond(i,k) = 0.
psmlt(i,k) = 0.
pgmlt(i,k) = 0.
pseml(i,k) = 0.
pgeml(i,k) = 0.
psevp(i,k) = 0.
pgevp(i,k) = 0.
falk(i,k,1) = 0.
falk(i,k,2) = 0.
falk(i,k,3) = 0.
fall(i,k,1) = 0.
fall(i,k,2) = 0.
fall(i,k,3) = 0.
fallc(i,k) = 0.
falkc(i,k) = 0.
xni(i,k) = 1.e3
enddo
enddo
do k = kts, kte
do i = its, ite
temp = (den(i,k)*max(qci(i,k,2),qmin))
temp = sqrt(sqrt(temp*temp*temp))
xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
enddo
enddo
do k = kts, kte
do i = its, ite
qrs_tmp(i,k,1) = qrs(i,k,1)
qrs_tmp(i,k,2) = qrs(i,k,2)
qrs_tmp(i,k,3) = qrs(i,k,3)
enddo
enddo
call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
work1,its,ite,kts,kte)
do k = kte, kts, -1
do i = its, ite
workr(i,k) = work1(i,k,1)
qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
IF ( qsum(i,k) .gt. 1.e-15 ) THEN
worka(i,k) = (work1(i,k,2)*qrs(i,k,2) + work1(i,k,3)*qrs(i,k,3)) &
/qsum(i,k)
ELSE
worka(i,k) = 0.
ENDIF
denqrs1(i,k) = den(i,k)*qrs(i,k,1)
denqrs2(i,k) = den(i,k)*qrs(i,k,2)
denqrs3(i,k) = den(i,k)*qrs(i,k,3)
if(qrs(i,k,1).le.0.0) workr(i,k) = 0.0
enddo
enddo
call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,workr,denqrs1, &
delqrs1,dtcld,1,1)
call nislfv_rain_plm6(idim,kdim,den_tmp,denfac,t,delz_tmp,worka, &
denqrs2,denqrs3,delqrs2,delqrs3,dtcld,1,1)
do k = kts, kte
do i = its, ite
qrs(i,k,1) = max(denqrs1(i,k)/den(i,k),0.)
qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
qrs(i,k,3) = max(denqrs3(i,k)/den(i,k),0.)
fall(i,k,1) = denqrs1(i,k)*workr(i,k)/delz(i,k)
fall(i,k,2) = denqrs2(i,k)*worka(i,k)/delz(i,k)
fall(i,k,3) = denqrs3(i,k)*worka(i,k)/delz(i,k)
enddo
enddo
do i = its, ite
fall(i,1,1) = delqrs1(i)/delz(i,1)/dtcld
fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
fall(i,1,3) = delqrs3(i)/delz(i,1)/dtcld
enddo
do k = kts, kte
do i = its, ite
qrs_tmp(i,k,1) = qrs(i,k,1)
qrs_tmp(i,k,2) = qrs(i,k,2)
qrs_tmp(i,k,3) = qrs(i,k,3)
enddo
enddo
call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
work1,its,ite,kts,kte)
do k = kte, kts, -1
do i = its, ite
supcol = t0c-t(i,k)
n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
if(t(i,k).gt.t0c) then
xlf = xlf0
work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
if(qrs(i,k,2).gt.0.) then
coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. &
*n0sfac(i,k)*(precs1*rslope2(i,k,2) &
+precs2*work2(i,k)*coeres)/den(i,k)
psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
-qrs(i,k,2)/mstep(i)),0.)
qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
endif
if(qrs(i,k,3).gt.0.) then
coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf &
*(t0c-t(i,k))*(precg1*rslope2(i,k,3) &
+precg2*work2(i,k)*coeres)/den(i,k)
pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i), &
-qrs(i,k,3)/mstep(i)),0.)
qrs(i,k,3) = qrs(i,k,3) + pgmlt(i,k)
qrs(i,k,1) = qrs(i,k,1) - pgmlt(i,k)
t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
endif
endif
enddo
enddo
do k = kte, kts, -1
do i = its, ite
if(qci(i,k,2).le.0.) then
work1c(i,k) = 0.
else
xmi = den(i,k)*qci(i,k,2)/xni(i,k)
diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
endif
enddo
enddo
do k = kte, kts, -1
do i = its, ite
denqci(i,k) = den(i,k)*qci(i,k,2)
enddo
enddo
call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci, &
delqi,dtcld,1,0)
do k = kts, kte
do i = its, ite
qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
enddo
enddo
do i = its, ite
fallc(i,1) = delqi(i)/delz(i,1)/dtcld
enddo
do i = its, ite
fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
fallsum_qg = fall(i,kts,3)
if(fallsum.gt.0.) then
rainncv(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rainncv(i)
rain(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rain(i)
endif
if(fallsum_qsi.gt.0.) then
tstepsnow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
+tstepsnow(i)
IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
snowncv(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
+snowncv(i,lat)
snow(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i,lat)
ENDIF
endif
if(fallsum_qg.gt.0.) then
tstepgraup(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. &
+tstepgraup(i)
IF ( PRESENT (graupelncv) .AND. PRESENT (graupel)) THEN
graupelncv(i,lat) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. &
+ graupelncv(i,lat)
graupel(i,lat) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. + graupel(i,lat)
ENDIF
endif
IF ( PRESENT (snowncv)) THEN
if(fallsum.gt.0.)sr(i)=(snowncv(i,lat) + graupelncv(i,lat))/(rainncv(i)+1.e-12)
ELSE
if(fallsum.gt.0.)sr(i)=(tstepsnow(i) + tstepgraup(i))/(rainncv(i)+1.e-12)
ENDIF
enddo
do k = kts, kte
do i = its, ite
supcol = t0c-t(i,k)
xlf = xls-xl(i,k)
if(supcol.lt.0.) xlf = xlf0
if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
qci(i,k,2) = 0.
endif
if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
qci(i,k,1) = 0.
endif
if(supcol.gt.0..and.qci(i,k,1).gt.qmin) then
supcolt=min(supcol,50.)
pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.) &
*den(i,k)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
qci(i,k,2) = qci(i,k,2) + pfrzdtc
t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
qci(i,k,1) = qci(i,k,1)-pfrzdtc
endif
if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
temp = rslope3(i,k,1)
temp = temp*temp*rslope(i,k,1)
supcolt=min(supcol,50.)
pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k) &
*(exp(pfrz2*supcolt)-1.)*temp*dtcld, &
qrs(i,k,1))
qrs(i,k,3) = qrs(i,k,3) + pfrzdtr
t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
endif
enddo
enddo
do k = kts, kte
do i = its, ite
qrs_tmp(i,k,1) = qrs(i,k,1)
qrs_tmp(i,k,2) = qrs(i,k,2)
qrs_tmp(i,k,3) = qrs(i,k,3)
enddo
enddo
call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
work1,its,ite,kts,kte)
do k = kts, kte
do i = its, ite
work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
enddo
enddo
do k = kts, kte
do i = its, ite
supsat = max(q(i,k),qmin)-qs(i,k,1)
satdt = supsat/dtcld
if(qci(i,k,1).gt.qc0) then
praut(i,k) = qck1*qci(i,k,1)**(7./3.)
praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
endif
if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1) &
*qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
endif
if(qrs(i,k,1).gt.0.) then
coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1) &
+precr2*work2(i,k)*coeres)/work1(i,k,1)
if(prevp(i,k).lt.0.) then
prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
prevp(i,k) = max(prevp(i,k),satdt/2)
else
prevp(i,k) = min(prevp(i,k),satdt/2)
endif
endif
enddo
enddo
do k = kts, kte
do i = its, ite
supcol = t0c-t(i,k)
n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
supsat = max(q(i,k),qmin)-qs(i,k,2)
satdt = supsat/dtcld
ifsat = 0
temp = (den(i,k)*max(qci(i,k,2),qmin))
temp = sqrt(sqrt(temp*temp*temp))
xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
eacrs = exp(0.07*(-supcol))
xmi = den(i,k)*qci(i,k,2)/xni(i,k)
diameter = min(dicon * sqrt(xmi),dimax)
vt2i = 1.49e4*diameter**1.31
vt2r=pvtr*rslopeb(i,k,1)*denfac(i,k)
vt2s=pvts*rslopeb(i,k,2)*denfac(i,k)
vt2g=pvtg*rslopeb(i,k,3)*denfac(i,k)
qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
if(qsum(i,k) .gt. 1.e-15) then
vt2ave=(vt2s*qrs(i,k,2)+vt2g*qrs(i,k,3))/(qsum(i,k))
else
vt2ave=0.
endif
if(supcol.gt.0.and.qci(i,k,2).gt.qmin) then
if(qrs(i,k,1).gt.qcrmin) then
acrfac = 2.*rslope3(i,k,1)+2.*diameter*rslope2(i,k,1) &
+diameter**2*rslope(i,k,1)
praci(i,k) = pi*qci(i,k,2)*n0r*abs(vt2r-vt2i)*acrfac/4.
praci(i,k) = min(praci(i,k),qci(i,k,2)/dtcld)
piacr(i,k) = pi**2*avtr*n0r*denr*xni(i,k)*denfac(i,k) &
*g6pbr*rslope3(i,k,1)*rslope3(i,k,1) &
*rslopeb(i,k,1)/24./den(i,k)
piacr(i,k) = min(piacr(i,k),qrs(i,k,1)/dtcld)
endif
if(qrs(i,k,2).gt.qcrmin) then
acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
+diameter**2*rslope(i,k,2)
psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k) &
*abs(vt2ave-vt2i)*acrfac/4.
psaci(i,k) = min(psaci(i,k),qci(i,k,2)/dtcld)
endif
if(qrs(i,k,3).gt.qcrmin) then
egi = exp(0.07*(-supcol))
acrfac = 2.*rslope3(i,k,3)+2.*diameter*rslope2(i,k,3) &
+diameter**2*rslope(i,k,3)
pgaci(i,k) = pi*egi*qci(i,k,2)*n0g*abs(vt2ave-vt2i)*acrfac/4.
pgaci(i,k) = min(pgaci(i,k),qci(i,k,2)/dtcld)
endif
endif
if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) &
*qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
endif
if(qrs(i,k,3).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
pgacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3) &
*qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
endif
if(qsum(i,k) .gt. 1.e-15) then
paacw(i,k) = (qrs(i,k,2)*psacw(i,k)+qrs(i,k,3)*pgacw(i,k)) &
/(qsum(i,k))
endif
if(qrs(i,k,2).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
if(supcol.gt.0) then
acrfac = 5.*rslope3(i,k,2)*rslope3(i,k,2)*rslope(i,k,1) &
+2.*rslope3(i,k,2)*rslope2(i,k,2)*rslope2(i,k,1) &
+.5*rslope2(i,k,2)*rslope2(i,k,2)*rslope3(i,k,1)
pracs(i,k) = pi**2*n0r*n0s*n0sfac(i,k)*abs(vt2r-vt2ave) &
*(dens/den(i,k))*acrfac
pracs(i,k) = min(pracs(i,k),qrs(i,k,2)/dtcld)
endif
acrfac = 5.*rslope3(i,k,1)*rslope3(i,k,1)*rslope(i,k,2) &
+2.*rslope3(i,k,1)*rslope2(i,k,1)*rslope2(i,k,2) &
+.5*rslope2(i,k,1)*rslope2(i,k,1)*rslope3(i,k,2)
psacr(i,k) = pi**2*n0r*n0s*n0sfac(i,k)*abs(vt2ave-vt2r) &
*(denr/den(i,k))*acrfac
psacr(i,k) = min(psacr(i,k),qrs(i,k,1)/dtcld)
endif
if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
acrfac = 5.*rslope3(i,k,1)*rslope3(i,k,1)*rslope(i,k,3) &
+2.*rslope3(i,k,1)*rslope2(i,k,1)*rslope2(i,k,3) &
+.5*rslope2(i,k,1)*rslope2(i,k,1)*rslope3(i,k,3)
pgacr(i,k) = pi**2*n0r*n0g*abs(vt2ave-vt2r)*(denr/den(i,k)) &
*acrfac
pgacr(i,k) = min(pgacr(i,k),qrs(i,k,1)/dtcld)
endif
if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,2).gt.qcrmin) then
pgacs(i,k) = 0.
endif
if(supcol.le.0) then
xlf = xlf0
if(qrs(i,k,2).gt.0.) &
pseml(i,k) = min(max(cliq*supcol*(paacw(i,k)+psacr(i,k)) &
/xlf,-qrs(i,k,2)/dtcld),0.)
if(qrs(i,k,3).gt.0.) &
pgeml(i,k) = min(max(cliq*supcol*(paacw(i,k)+pgacr(i,k)) &
/xlf,-qrs(i,k,3)/dtcld),0.)
endif
if(supcol.gt.0) then
if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
supice = satdt-prevp(i,k)
if(pidep(i,k).lt.0.) then
pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
else
pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
endif
if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
endif
if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2) &
+ precs2*work2(i,k)*coeres)/work1(i,k,2)
supice = satdt-prevp(i,k)-pidep(i,k)
if(psdep(i,k).lt.0.) then
psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
else
psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
endif
if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
ifsat = 1
endif
if(qrs(i,k,3).gt.0..and.ifsat.ne.1) then
coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
pgdep(i,k) = (rh(i,k,2)-1.)*(precg1*rslope2(i,k,3) &
+precg2*work2(i,k)*coeres)/work1(i,k,2)
supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
if(pgdep(i,k).lt.0.) then
pgdep(i,k) = max(pgdep(i,k),-qrs(i,k,3)/dtcld)
pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
else
pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
endif
if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge. &
abs(satdt)) ifsat = 1
endif
if(supsat.gt.0.and.ifsat.ne.1) then
supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
xni0 = 1.e3*exp(0.1*supcol)
roqi0 = 4.92e-11*xni0**1.33
pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))/dtcld)
pigen(i,k) = min(min(pigen(i,k),satdt),supice)
endif
if(qci(i,k,2).gt.0.) then
qimax = roqimax/den(i,k)
psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
endif
if(qrs(i,k,2).gt.0.) then
alpha2 = 1.e-3*exp(0.09*(-supcol))
pgaut(i,k) = min(max(0.,alpha2*(qrs(i,k,2)-qs0)),qrs(i,k,2)/dtcld)
endif
endif
if(supcol.lt.0.) then
if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) then
coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
psevp(i,k) = (rh(i,k,1)-1.)*n0sfac(i,k)*(precs1 &
*rslope2(i,k,2)+precs2*work2(i,k) &
*coeres)/work1(i,k,1)
psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
endif
if(qrs(i,k,3).gt.0..and.rh(i,k,1).lt.1.) then
coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
pgevp(i,k) = (rh(i,k,1)-1.)*(precg1*rslope2(i,k,3) &
+precg2*work2(i,k)*coeres)/work1(i,k,1)
pgevp(i,k) = min(max(pgevp(i,k),-qrs(i,k,3)/dtcld),0.)
endif
endif
enddo
enddo
do k = kts, kte
do i = its, ite
delta2=0.
delta3=0.
if(qrs(i,k,1).lt.1.e-4.and.qrs(i,k,2).lt.1.e-4) delta2=1.
if(qrs(i,k,1).lt.1.e-4) delta3=1.
if(t(i,k).le.t0c) then
value = max(qmin,qci(i,k,1))
source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
if (source.gt.value) then
factor = value/source
praut(i,k) = praut(i,k)*factor
pracw(i,k) = pracw(i,k)*factor
paacw(i,k) = paacw(i,k)*factor
endif
value = max(qmin,qci(i,k,2))
source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k)+psaci(i,k) &
+pgaci(i,k))*dtcld
if (source.gt.value) then
factor = value/source
psaut(i,k) = psaut(i,k)*factor
pigen(i,k) = pigen(i,k)*factor
pidep(i,k) = pidep(i,k)*factor
praci(i,k) = praci(i,k)*factor
psaci(i,k) = psaci(i,k)*factor
pgaci(i,k) = pgaci(i,k)*factor
endif
value = max(qmin,qrs(i,k,1))
source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k)+psacr(i,k) &
+pgacr(i,k))*dtcld
if (source.gt.value) then
factor = value/source
praut(i,k) = praut(i,k)*factor
prevp(i,k) = prevp(i,k)*factor
pracw(i,k) = pracw(i,k)*factor
piacr(i,k) = piacr(i,k)*factor
psacr(i,k) = psacr(i,k)*factor
pgacr(i,k) = pgacr(i,k)*factor
endif
value = max(qmin,qrs(i,k,2))
source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k)+piacr(i,k) &
*delta3+praci(i,k)*delta3-pracs(i,k)*(1.-delta2) &
+psacr(i,k)*delta2+psaci(i,k)-pgacs(i,k) )*dtcld
if (source.gt.value) then
factor = value/source
psdep(i,k) = psdep(i,k)*factor
psaut(i,k) = psaut(i,k)*factor
pgaut(i,k) = pgaut(i,k)*factor
paacw(i,k) = paacw(i,k)*factor
piacr(i,k) = piacr(i,k)*factor
praci(i,k) = praci(i,k)*factor
psaci(i,k) = psaci(i,k)*factor
pracs(i,k) = pracs(i,k)*factor
psacr(i,k) = psacr(i,k)*factor
pgacs(i,k) = pgacs(i,k)*factor
endif
value = max(qmin,qrs(i,k,3))
source = -(pgdep(i,k)+pgaut(i,k) &
+piacr(i,k)*(1.-delta3)+praci(i,k)*(1.-delta3) &
+psacr(i,k)*(1.-delta2)+pracs(i,k)*(1.-delta2) &
+pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
if (source.gt.value) then
factor = value/source
pgdep(i,k) = pgdep(i,k)*factor
pgaut(i,k) = pgaut(i,k)*factor
piacr(i,k) = piacr(i,k)*factor
praci(i,k) = praci(i,k)*factor
psacr(i,k) = psacr(i,k)*factor
pracs(i,k) = pracs(i,k)*factor
paacw(i,k) = paacw(i,k)*factor
pgaci(i,k) = pgaci(i,k)*factor
pgacr(i,k) = pgacr(i,k)*factor
pgacs(i,k) = pgacs(i,k)*factor
endif
work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k)+pidep(i,k))
q(i,k) = q(i,k)+work2(i,k)*dtcld
qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
+paacw(i,k)+paacw(i,k))*dtcld,0.)
qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
+prevp(i,k)-piacr(i,k)-pgacr(i,k) &
-psacr(i,k))*dtcld,0.)
qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+praci(i,k) &
+psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k)) &
*dtcld,0.)
qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+paacw(i,k) &
-pgaut(i,k)+piacr(i,k)*delta3 &
+praci(i,k)*delta3+psaci(i,k)-pgacs(i,k) &
-pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2) &
*dtcld,0.)
qrs(i,k,3) = max(qrs(i,k,3)+(pgdep(i,k)+pgaut(i,k) &
+piacr(i,k)*(1.-delta3) &
+praci(i,k)*(1.-delta3)+psacr(i,k)*(1.-delta2) &
+pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k) &
+pgacr(i,k)+pgacs(i,k))*dtcld,0.)
xlf = xls-xl(i,k)
xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k)) &
-xl(i,k)*prevp(i,k)-xlf*(piacr(i,k)+paacw(i,k) &
+paacw(i,k)+pgacr(i,k)+psacr(i,k))
t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
else
value = max(qmin,qci(i,k,1))
source=(praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
if (source.gt.value) then
factor = value/source
praut(i,k) = praut(i,k)*factor
pracw(i,k) = pracw(i,k)*factor
paacw(i,k) = paacw(i,k)*factor
endif
value = max(qmin,qrs(i,k,1))
source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k)-pracw(i,k) &
-paacw(i,k)-prevp(i,k))*dtcld
if (source.gt.value) then
factor = value/source
praut(i,k) = praut(i,k)*factor
prevp(i,k) = prevp(i,k)*factor
pracw(i,k) = pracw(i,k)*factor
paacw(i,k) = paacw(i,k)*factor
pseml(i,k) = pseml(i,k)*factor
pgeml(i,k) = pgeml(i,k)*factor
endif
value = max(qcrmin,qrs(i,k,2))
source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
if (source.gt.value) then
factor = value/source
pgacs(i,k) = pgacs(i,k)*factor
psevp(i,k) = psevp(i,k)*factor
pseml(i,k) = pseml(i,k)*factor
endif
value = max(qcrmin,qrs(i,k,3))
source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
if (source.gt.value) then
factor = value/source
pgacs(i,k) = pgacs(i,k)*factor
pgevp(i,k) = pgevp(i,k)*factor
pgeml(i,k) = pgeml(i,k)*factor
endif
work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
q(i,k) = q(i,k)+work2(i,k)*dtcld
qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
+paacw(i,k)+paacw(i,k))*dtcld,0.)
qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
+prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k) &
-pgeml(i,k))*dtcld,0.)
qrs(i,k,2) = max(qrs(i,k,2)+(psevp(i,k)-pgacs(i,k) &
+pseml(i,k))*dtcld,0.)
qrs(i,k,3) = max(qrs(i,k,3)+(pgacs(i,k)+pgevp(i,k) &
+pgeml(i,k))*dtcld,0.)
xlf = xls-xl(i,k)
xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k)) &
-xlf*(pseml(i,k)+pgeml(i,k))
t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
endif
enddo
enddo
hsub = xls
hvap = xlv0
cvap = cpv
ttp=t0c+0.01
dldt=cvap-cliq
xa=-dldt/rv
xb=xa+hvap/(rv*ttp)
dldti=cvap-cice
xai=-dldti/rv
xbi=xai+hsub/(rv*ttp)
do k = kts, kte
do i = its, ite
tr=ttp/t(i,k)
qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
qs(i,k,1) = max(qs(i,k,1),qmin)
tr=ttp/t(i,k)
if(t(i,k).lt.ttp) then
qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
else
qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
endif
qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
qs(i,k,2) = max(qs(i,k,2),qmin)
enddo
enddo
do k = kts, kte
do i = its, ite
work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
work2(i,k) = qci(i,k,1)+work1(i,k,1)
pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
q(i,k) = q(i,k)-pcond(i,k)*dtcld
qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
enddo
enddo
do k = kts, kte
do i = its, ite
if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
enddo
enddo
enddo
END SUBROUTINE wsm62d
REAL FUNCTION rgmma(x)
IMPLICIT NONE
REAL :: euler
PARAMETER (euler=0.577215664901532)
REAL :: x, y
INTEGER :: i
if(x.eq.1.)then
rgmma=0.
else
rgmma=x*exp(euler*x)
do i=1,10000
y=float(i)
rgmma=rgmma*(1.000+x/y)*exp(-x/y)
enddo
rgmma=1./rgmma
endif
END FUNCTION rgmma
REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
IMPLICIT NONE
REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
xai,xbi,ttp,tr
INTEGER ice
ttp=t0c+0.01
dldt=cvap-cliq
xa=-dldt/rv
xb=xa+hvap/(rv*ttp)
dldti=cvap-cice
xai=-dldti/rv
xbi=xai+hsub/(rv*ttp)
tr=ttp/t
if(t.lt.ttp.and.ice.eq.1) then
fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
else
fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
endif
END FUNCTION fpvs
SUBROUTINE wsm6init(den0,denr,dens,cl,cpv,hail_opt,allowed_to_read)
IMPLICIT NONE
REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
INTEGER, INTENT(IN) :: hail_opt
LOGICAL, INTENT(IN) :: allowed_to_read
IF (hail_opt .eq. 1) THEN
n0g = 4.e4
deng = 700.
avtg = 285.0
bvtg = 0.8
lamdagmax = 2.e4
ELSE
n0g = 4.e6
deng = 500
avtg = 330.0
bvtg = 0.8
lamdagmax = 6.e4
ENDIF
pi = 4.*atan(1.)
xlv1 = cl-cpv
qc0 = 4./3.*pi*denr*r0**3*xncr/den0
qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.)
pidnc = pi*denr/6.
bvtr1 = 1.+bvtr
bvtr2 = 2.5+.5*bvtr
bvtr3 = 3.+bvtr
bvtr4 = 4.+bvtr
bvtr6 = 6.+bvtr
g1pbr = rgmma(bvtr1)
g3pbr = rgmma(bvtr3)
g4pbr = rgmma(bvtr4)
g6pbr = rgmma(bvtr6)
g5pbro2 = rgmma(bvtr2)
pvtr = avtr*g4pbr/6.
eacrr = 1.0
pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
precr1 = 2.*pi*n0r*.78
precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
roqimax = 2.08e22*dimax**8
bvts1 = 1.+bvts
bvts2 = 2.5+.5*bvts
bvts3 = 3.+bvts
bvts4 = 4.+bvts
g1pbs = rgmma(bvts1)
g3pbs = rgmma(bvts3)
g4pbs = rgmma(bvts4)
g5pbso2 = rgmma(bvts2)
pvts = avts*g4pbs/6.
pacrs = pi*n0s*avts*g3pbs*.25
precs1 = 4.*n0s*.65
precs2 = 4.*n0s*.44*avts**.5*g5pbso2
pidn0r = pi*denr*n0r
pidn0s = pi*dens*n0s
pacrc = pi*n0s*avts*g3pbs*.25*eacrc
bvtg1 = 1.+bvtg
bvtg2 = 2.5+.5*bvtg
bvtg3 = 3.+bvtg
bvtg4 = 4.+bvtg
g1pbg = rgmma(bvtg1)
g3pbg = rgmma(bvtg3)
g4pbg = rgmma(bvtg4)
pacrg = pi*n0g*avtg*g3pbg*.25
g5pbgo2 = rgmma(bvtg2)
pvtg = avtg*g4pbg/6.
precg1 = 2.*pi*n0g*.78
precg2 = 2.*pi*n0g*.31*avtg**.5*g5pbgo2
pidn0g = pi*deng*n0g
rslopermax = 1./lamdarmax
rslopesmax = 1./lamdasmax
rslopegmax = 1./lamdagmax
rsloperbmax = rslopermax ** bvtr
rslopesbmax = rslopesmax ** bvts
rslopegbmax = rslopegmax ** bvtg
rsloper2max = rslopermax * rslopermax
rslopes2max = rslopesmax * rslopesmax
rslopeg2max = rslopegmax * rslopegmax
rsloper3max = rsloper2max * rslopermax
rslopes3max = rslopes2max * rslopesmax
rslopeg3max = rslopeg2max * rslopegmax
xam_r = PI*denr/6.
xbm_r = 3.
xmu_r = 0.
xam_s = PI*dens/6.
xbm_s = 3.
xmu_s = 0.
xam_g = PI*deng/6.
xbm_g = 3.
xmu_g = 0.
call radar_init
END SUBROUTINE wsm6init
subroutine slope_wsm6(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
vt,its,ite,kts,kte)
IMPLICIT NONE
INTEGER :: its,ite, jts,jte, kts,kte
REAL, DIMENSION( its:ite , kts:kte,3) :: &
qrs, &
rslope, &
rslopeb, &
rslope2, &
rslope3, &
vt
REAL, DIMENSION( its:ite , kts:kte) :: &
den, &
denfac, &
t
REAL, PARAMETER :: t0c = 273.15
REAL, DIMENSION( its:ite , kts:kte ) :: &
n0sfac
REAL :: lamdar, lamdas, lamdag, x, y, z, supcol
integer :: i, j, k
lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y)))
lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))
lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y)))
do k = kts, kte
do i = its, ite
supcol = t0c-t(i,k)
n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
if(qrs(i,k,1).le.qcrmin)then
rslope(i,k,1) = rslopermax
rslopeb(i,k,1) = rsloperbmax
rslope2(i,k,1) = rsloper2max
rslope3(i,k,1) = rsloper3max
else
rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
rslopeb(i,k,1) = rslope(i,k,1)**bvtr
rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
endif
if(qrs(i,k,2).le.qcrmin)then
rslope(i,k,2) = rslopesmax
rslopeb(i,k,2) = rslopesbmax
rslope2(i,k,2) = rslopes2max
rslope3(i,k,2) = rslopes3max
else
rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
rslopeb(i,k,2) = rslope(i,k,2)**bvts
rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
endif
if(qrs(i,k,3).le.qcrmin)then
rslope(i,k,3) = rslopegmax
rslopeb(i,k,3) = rslopegbmax
rslope2(i,k,3) = rslopeg2max
rslope3(i,k,3) = rslopeg3max
else
rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
rslopeb(i,k,3) = rslope(i,k,3)**bvtg
rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
endif
vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
vt(i,k,3) = pvtg*rslopeb(i,k,3)*denfac(i,k)
if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
if(qrs(i,k,3).le.0.0) vt(i,k,3) = 0.0
enddo
enddo
END subroutine slope_wsm6
subroutine slope_rain(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
vt,its,ite,kts,kte)
IMPLICIT NONE
INTEGER :: its,ite, jts,jte, kts,kte
REAL, DIMENSION( its:ite , kts:kte) :: &
qrs, &
rslope, &
rslopeb, &
rslope2, &
rslope3, &
vt, &
den, &
denfac, &
t
REAL, PARAMETER :: t0c = 273.15
REAL, DIMENSION( its:ite , kts:kte ) :: &
n0sfac
REAL :: lamdar, x, y, z, supcol
integer :: i, j, k
lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y)))
do k = kts, kte
do i = its, ite
if(qrs(i,k).le.qcrmin)then
rslope(i,k) = rslopermax
rslopeb(i,k) = rsloperbmax
rslope2(i,k) = rsloper2max
rslope3(i,k) = rsloper3max
else
rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
rslopeb(i,k) = rslope(i,k)**bvtr
rslope2(i,k) = rslope(i,k)*rslope(i,k)
rslope3(i,k) = rslope2(i,k)*rslope(i,k)
endif
vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
if(qrs(i,k).le.0.0) vt(i,k) = 0.0
enddo
enddo
END subroutine slope_rain
subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
vt,its,ite,kts,kte)
IMPLICIT NONE
INTEGER :: its,ite, jts,jte, kts,kte
REAL, DIMENSION( its:ite , kts:kte) :: &
qrs, &
rslope, &
rslopeb, &
rslope2, &
rslope3, &
vt, &
den, &
denfac, &
t
REAL, PARAMETER :: t0c = 273.15
REAL, DIMENSION( its:ite , kts:kte ) :: &
n0sfac
REAL :: lamdas, x, y, z, supcol
integer :: i, j, k
lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y)))
do k = kts, kte
do i = its, ite
supcol = t0c-t(i,k)
n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
if(qrs(i,k).le.qcrmin)then
rslope(i,k) = rslopesmax
rslopeb(i,k) = rslopesbmax
rslope2(i,k) = rslopes2max
rslope3(i,k) = rslopes3max
else
rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
rslopeb(i,k) = rslope(i,k)**bvts
rslope2(i,k) = rslope(i,k)*rslope(i,k)
rslope3(i,k) = rslope2(i,k)*rslope(i,k)
endif
vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
if(qrs(i,k).le.0.0) vt(i,k) = 0.0
enddo
enddo
END subroutine slope_snow
subroutine slope_graup(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
vt,its,ite,kts,kte)
IMPLICIT NONE
INTEGER :: its,ite, jts,jte, kts,kte
REAL, DIMENSION( its:ite , kts:kte) :: &
qrs, &
rslope, &
rslopeb, &
rslope2, &
rslope3, &
vt, &
den, &
denfac, &
t
REAL, PARAMETER :: t0c = 273.15
REAL, DIMENSION( its:ite , kts:kte ) :: &
n0sfac
REAL :: lamdag, x, y, z, supcol
integer :: i, j, k
lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y)))
do k = kts, kte
do i = its, ite
if(qrs(i,k).le.qcrmin)then
rslope(i,k) = rslopegmax
rslopeb(i,k) = rslopegbmax
rslope2(i,k) = rslopeg2max
rslope3(i,k) = rslopeg3max
else
rslope(i,k) = 1./lamdag(qrs(i,k),den(i,k))
rslopeb(i,k) = rslope(i,k)**bvtg
rslope2(i,k) = rslope(i,k)*rslope(i,k)
rslope3(i,k) = rslope2(i,k)*rslope(i,k)
endif
vt(i,k) = pvtg*rslopeb(i,k)*denfac(i,k)
if(qrs(i,k).le.0.0) vt(i,k) = 0.0
enddo
enddo
END subroutine slope_graup
SUBROUTINE nislfv_rain_plm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
implicit none
integer im,km,id
real dt
real dzl(im,km),wwl(im,km),rql(im,km),precip(im)
real denl(im,km),denfacl(im,km),tkl(im,km)
integer i,k,n,m,kk,kb,kt,iter
real tl,tl2,qql,dql,qqd
real th,th2,qqh,dqh
real zsum,qsum,dim,dip,c1,con1,fa1,fa2
real allold, allnew, zz, dzamin, cflmax, decfl
real dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
real den(km), denfac(km), tk(km)
real wi(km+1), zi(km+1), za(km+1)
real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
precip(:) = 0.0
i_loop : do i=1,im
dz(:) = dzl(i,:)
qq(:) = rql(i,:)
ww(:) = wwl(i,:)
den(:) = denl(i,:)
denfac(:) = denfacl(i,:)
tk(:) = tkl(i,:)
allold = 0.0
do k=1,km
allold = allold + qq(k)
enddo
if(allold.le.0.0) then
cycle i_loop
endif
zi(1)=0.0
do k=1,km
zi(k+1) = zi(k)+dz(k)
enddo
wd(:) = ww(:)
n=1
100 continue
wi(1) = ww(1)
wi(km+1) = ww(km)
do k=2,km
wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
enddo
fa1 = 9./16.
fa2 = 1./16.
wi(1) = ww(1)
wi(2) = 0.5*(ww(2)+ww(1))
do k=3,km-1
wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
enddo
wi(km) = 0.5*(ww(km)+ww(km-1))
wi(km+1) = ww(km)
do k=2,km
if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
enddo
con1 = 0.05
do k=km,1,-1
decfl = (wi(k+1)-wi(k))*dt/dz(k)
if( decfl .gt. con1 ) then
wi(k) = wi(k+1) - con1*dz(k)/dt
endif
enddo
do k=1,km+1
za(k) = zi(k) - wi(k)*dt
enddo
do k=1,km
dza(k) = za(k+1)-za(k)
enddo
dza(km+1) = zi(km+1) - za(km+1)
do k=1,km
qa(k) = qq(k)*dz(k)/dza(k)
qr(k) = qa(k)/den(k)
enddo
qa(km+1) = 0.0
if( n.le.iter ) then
call slope_rain(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
do k=1,km
ww(k) = 0.5* ( wd(k)+wa(k) )
enddo
was(:) = wa(:)
n=n+1
go to 100
endif
do k=2,km
dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
if( dip*dim.le.0.0 ) then
qmi(k)=qa(k)
qpi(k)=qa(k)
else
qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
qmi(k)=2.0*qa(k)-qpi(k)
if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
qpi(k) = qa(k)
qmi(k) = qa(k)
endif
endif
enddo
qpi(1)=qa(1)
qmi(1)=qa(1)
qmi(km+1)=qa(km+1)
qpi(km+1)=qa(km+1)
qn = 0.0
kb=1
kt=1
intp : do k=1,km
kb=max(kb-1,1)
kt=max(kt-1,1)
if( zi(k).ge.za(km+1) ) then
exit intp
else
find_kb : do kk=kb,km
if( zi(k).le.za(kk+1) ) then
kb = kk
exit find_kb
else
cycle find_kb
endif
enddo find_kb
find_kt : do kk=kt,km
if( zi(k+1).le.za(kk) ) then
kt = kk
exit find_kt
else
cycle find_kt
endif
enddo find_kt
kt = kt - 1
if( kt.eq.kb ) then
tl=(zi(k)-za(kb))/dza(kb)
th=(zi(k+1)-za(kb))/dza(kb)
tl2=tl*tl
th2=th*th
qqd=0.5*(qpi(kb)-qmi(kb))
qqh=qqd*th2+qmi(kb)*th
qql=qqd*tl2+qmi(kb)*tl
qn(k) = (qqh-qql)/(th-tl)
else if( kt.gt.kb ) then
tl=(zi(k)-za(kb))/dza(kb)
tl2=tl*tl
qqd=0.5*(qpi(kb)-qmi(kb))
qql=qqd*tl2+qmi(kb)*tl
dql = qa(kb)-qql
zsum = (1.-tl)*dza(kb)
qsum = dql*dza(kb)
if( kt-kb.gt.1 ) then
do m=kb+1,kt-1
zsum = zsum + dza(m)
qsum = qsum + qa(m) * dza(m)
enddo
endif
th=(zi(k+1)-za(kt))/dza(kt)
th2=th*th
qqd=0.5*(qpi(kt)-qmi(kt))
dqh=qqd*th2+qmi(kt)*th
zsum = zsum + th*dza(kt)
qsum = qsum + dqh*dza(kt)
qn(k) = qsum/zsum
endif
cycle intp
endif
enddo intp
sum_precip: do k=1,km
if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
precip(i) = precip(i) + qa(k)*dza(k)
cycle sum_precip
else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
precip(i) = precip(i) + qa(k)*(0.0-za(k))
exit sum_precip
endif
exit sum_precip
enddo sum_precip
rql(i,:) = qn(:)
enddo i_loop
END SUBROUTINE nislfv_rain_plm
SUBROUTINE nislfv_rain_plm6(im,km,denl,denfacl,tkl,dzl,wwl,rql,rql2, precip1, precip2,dt,id,iter)
implicit none
integer im,km,id
real dt
real dzl(im,km),wwl(im,km),rql(im,km),rql2(im,km),precip(im),precip1(im),precip2(im)
real denl(im,km),denfacl(im,km),tkl(im,km)
integer i,k,n,m,kk,kb,kt,iter,ist
real tl,tl2,qql,dql,qqd
real th,th2,qqh,dqh
real zsum,qsum,dim,dip,c1,con1,fa1,fa2
real allold, allnew, zz, dzamin, cflmax, decfl
real dz(km), ww(km), qq(km), qq2(km), wd(km), wa(km), wa2(km), was(km)
real den(km), denfac(km), tk(km)
real wi(km+1), zi(km+1), za(km+1)
real qn(km), qr(km),qr2(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
real dza(km+1), qa(km+1), qa2(km+1),qmi(km+1), qpi(km+1)
precip(:) = 0.0
precip1(:) = 0.0
precip2(:) = 0.0
i_loop : do i=1,im
dz(:) = dzl(i,:)
qq(:) = rql(i,:)
qq2(:) = rql2(i,:)
ww(:) = wwl(i,:)
den(:) = denl(i,:)
denfac(:) = denfacl(i,:)
tk(:) = tkl(i,:)
allold = 0.0
do k=1,km
allold = allold + qq(k) + qq2(k)
enddo
if(allold.le.0.0) then
cycle i_loop
endif
zi(1)=0.0
do k=1,km
zi(k+1) = zi(k)+dz(k)
enddo
wd(:) = ww(:)
n=1
100 continue
wi(1) = ww(1)
wi(km+1) = ww(km)
do k=2,km
wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
enddo
fa1 = 9./16.
fa2 = 1./16.
wi(1) = ww(1)
wi(2) = 0.5*(ww(2)+ww(1))
do k=3,km-1
wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
enddo
wi(km) = 0.5*(ww(km)+ww(km-1))
wi(km+1) = ww(km)
do k=2,km
if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
enddo
con1 = 0.05
do k=km,1,-1
decfl = (wi(k+1)-wi(k))*dt/dz(k)
if( decfl .gt. con1 ) then
wi(k) = wi(k+1) - con1*dz(k)/dt
endif
enddo
do k=1,km+1
za(k) = zi(k) - wi(k)*dt
enddo
do k=1,km
dza(k) = za(k+1)-za(k)
enddo
dza(km+1) = zi(km+1) - za(km+1)
do k=1,km
qa(k) = qq(k)*dz(k)/dza(k)
qa2(k) = qq2(k)*dz(k)/dza(k)
qr(k) = qa(k)/den(k)
qr2(k) = qa2(k)/den(k)
enddo
qa(km+1) = 0.0
qa2(km+1) = 0.0
if( n.le.iter ) then
call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
call slope_graup(qr2,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa2,1,1,1,km)
do k = 1, km
tmp(k) = max((qr(k)+qr2(k)), 1.E-15)
IF ( tmp(k) .gt. 1.e-15 ) THEN
wa(k) = (wa(k)*qr(k) + wa2(k)*qr2(k))/tmp(k)
ELSE
wa(k) = 0.
ENDIF
enddo
if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
do k=1,km
ww(k) = 0.5* ( wd(k)+wa(k) )
enddo
was(:) = wa(:)
n=n+1
go to 100
endif
ist_loop : do ist = 1, 2
if (ist.eq.2) then
qa(:) = qa2(:)
endif
precip(i) = 0.
do k=2,km
dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
if( dip*dim.le.0.0 ) then
qmi(k)=qa(k)
qpi(k)=qa(k)
else
qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
qmi(k)=2.0*qa(k)-qpi(k)
if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
qpi(k) = qa(k)
qmi(k) = qa(k)
endif
endif
enddo
qpi(1)=qa(1)
qmi(1)=qa(1)
qmi(km+1)=qa(km+1)
qpi(km+1)=qa(km+1)
qn = 0.0
kb=1
kt=1
intp : do k=1,km
kb=max(kb-1,1)
kt=max(kt-1,1)
if( zi(k).ge.za(km+1) ) then
exit intp
else
find_kb : do kk=kb,km
if( zi(k).le.za(kk+1) ) then
kb = kk
exit find_kb
else
cycle find_kb
endif
enddo find_kb
find_kt : do kk=kt,km
if( zi(k+1).le.za(kk) ) then
kt = kk
exit find_kt
else
cycle find_kt
endif
enddo find_kt
kt = kt - 1
if( kt.eq.kb ) then
tl=(zi(k)-za(kb))/dza(kb)
th=(zi(k+1)-za(kb))/dza(kb)
tl2=tl*tl
th2=th*th
qqd=0.5*(qpi(kb)-qmi(kb))
qqh=qqd*th2+qmi(kb)*th
qql=qqd*tl2+qmi(kb)*tl
qn(k) = (qqh-qql)/(th-tl)
else if( kt.gt.kb ) then
tl=(zi(k)-za(kb))/dza(kb)
tl2=tl*tl
qqd=0.5*(qpi(kb)-qmi(kb))
qql=qqd*tl2+qmi(kb)*tl
dql = qa(kb)-qql
zsum = (1.-tl)*dza(kb)
qsum = dql*dza(kb)
if( kt-kb.gt.1 ) then
do m=kb+1,kt-1
zsum = zsum + dza(m)
qsum = qsum + qa(m) * dza(m)
enddo
endif
th=(zi(k+1)-za(kt))/dza(kt)
th2=th*th
qqd=0.5*(qpi(kt)-qmi(kt))
dqh=qqd*th2+qmi(kt)*th
zsum = zsum + th*dza(kt)
qsum = qsum + dqh*dza(kt)
qn(k) = qsum/zsum
endif
cycle intp
endif
enddo intp
sum_precip: do k=1,km
if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
precip(i) = precip(i) + qa(k)*dza(k)
cycle sum_precip
else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
precip(i) = precip(i) + qa(k)*(0.0-za(k))
exit sum_precip
endif
exit sum_precip
enddo sum_precip
if(ist.eq.1) then
rql(i,:) = qn(:)
precip1(i) = precip(i)
else
rql2(i,:) = qn(:)
precip2(i) = precip(i)
endif
enddo ist_loop
enddo i_loop
END SUBROUTINE nislfv_rain_plm6
subroutine refl10cm_wsm6 (qv1d, qr1d, qs1d, qg1d, &
t1d, p1d, dBZ, kts, kte, ii, jj)
IMPLICIT NONE
INTEGER, INTENT(IN):: kts, kte, ii, jj
REAL, DIMENSION(kts:kte), INTENT(IN):: &
qv1d, qr1d, qs1d, qg1d, t1d, p1d
REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
REAL, DIMENSION(kts:kte):: rr, rs, rg
REAL:: temp_C
DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams, ilamg
DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s, N0_g
DOUBLE PRECISION:: lamr, lams, lamg
LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
DOUBLE PRECISION:: fmelt_s, fmelt_g
INTEGER:: i, k, k_0, kbot, n
LOGICAL:: melti
DOUBLE PRECISION:: cback, x, eta, f_d
REAL, PARAMETER:: R=287.
do k = kts, kte
dBZ(k) = -35.0
enddo
do k = kts, kte
temp(k) = t1d(k)
temp_C = min(-0.001, temp(K)-273.15)
qv(k) = MAX(1.E-10, qv1d(k))
pres(k) = p1d(k)
rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
if (qr1d(k) .gt. 1.E-9) then
rr(k) = qr1d(k)*rho(k)
N0_r(k) = n0r
lamr = (xam_r*xcrg(3)*N0_r(k)/rr(k))**(1./xcre(1))
ilamr(k) = 1./lamr
L_qr(k) = .true.
else
rr(k) = 1.E-12
L_qr(k) = .false.
endif
if (qs1d(k) .gt. 1.E-9) then
rs(k) = qs1d(k)*rho(k)
N0_s(k) = min(n0smax, n0s*exp(-alpha*temp_C))
lams = (xam_s*xcsg(3)*N0_s(k)/rs(k))**(1./xcse(1))
ilams(k) = 1./lams
L_qs(k) = .true.
else
rs(k) = 1.E-12
L_qs(k) = .false.
endif
if (qg1d(k) .gt. 1.E-9) then
rg(k) = qg1d(k)*rho(k)
N0_g(k) = n0g
lamg = (xam_g*xcgg(3)*N0_g(k)/rg(k))**(1./xcge(1))
ilamg(k) = 1./lamg
L_qg(k) = .true.
else
rg(k) = 1.E-12
L_qg(k) = .false.
endif
enddo
melti = .false.
k_0 = kts
do k = kte-1, kts, -1
if ( (temp(k).gt.273.15) .and. L_qr(k) &
.and. (L_qs(k+1).or.L_qg(k+1)) ) then
k_0 = MAX(k+1, k_0)
melti=.true.
goto 195
endif
enddo
195 continue
do k = kts, kte
ze_rain(k) = 1.e-22
ze_snow(k) = 1.e-22
ze_graupel(k) = 1.e-22
if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) &
* (xam_s/900.0)*(xam_s/900.0) &
* N0_s(k)*xcsg(4)*ilams(k)**xcse(4)
if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) &
* (xam_g/900.0)*(xam_g/900.0) &
* N0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
enddo
if (melti .and. k_0.ge.kts+1) then
do k = k_0-1, kts, -1
if (L_qs(k) .and. L_qs(k_0) ) then
fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
eta = 0.d0
lams = 1./ilams(k)
do n = 1, nrbins
x = xam_s * xxDs(n)**xbm_s
call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), &
fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
CBACK, mixingrulestring_s, matrixstring_s, &
inclusionstring_s, hoststring_s, &
hostmatrixstring_s, hostinclusionstring_s)
f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n))
eta = eta + f_d * CBACK * simpson(n) * xdts(n)
enddo
ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
endif
if (L_qg(k) .and. L_qg(k_0) ) then
fmelt_g = MAX(0.005d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0))
eta = 0.d0
lamg = 1./ilamg(k)
do n = 1, nrbins
x = xam_g * xxDg(n)**xbm_g
call rayleigh_soak_wetgraupel (x,DBLE(xocmg),DBLE(xobmg), &
fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
CBACK, mixingrulestring_g, matrixstring_g, &
inclusionstring_g, hoststring_g, &
hostmatrixstring_g, hostinclusionstring_g)
f_d = N0_g(k)*xxDg(n)**xmu_g * DEXP(-lamg*xxDg(n))
eta = eta + f_d * CBACK * simpson(n) * xdtg(n)
enddo
ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
endif
enddo
endif
do k = kte, kts, -1
dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
enddo
end subroutine refl10cm_wsm6
subroutine effectRad_wsm6 (t, qc, qi, qs, rho, qmin, t0c, &
re_qc, re_qi, re_qs, kts, kte, ii, jj)
implicit none
integer, intent(in) :: kts, kte, ii, jj
real, intent(in) :: qmin
real, intent(in) :: t0c
real, dimension( kts:kte ), intent(in):: t
real, dimension( kts:kte ), intent(in):: qc
real, dimension( kts:kte ), intent(in):: qi
real, dimension( kts:kte ), intent(in):: qs
real, dimension( kts:kte ), intent(in):: rho
real, dimension( kts:kte ), intent(inout):: re_qc
real, dimension( kts:kte ), intent(inout):: re_qi
real, dimension( kts:kte ), intent(inout):: re_qs
integer:: i,k
integer :: inu_c
real, dimension( kts:kte ):: ni
real, dimension( kts:kte ):: rqc
real, dimension( kts:kte ):: rqi
real, dimension( kts:kte ):: rni
real, dimension( kts:kte ):: rqs
real :: temp
real :: lamdac
real :: supcol, n0sfac, lamdas
real :: diai
logical :: has_qc, has_qi, has_qs
real, parameter :: R1 = 1.E-12
real, parameter :: R2 = 1.E-6
real, parameter :: bm_r = 3.0
real, parameter :: obmr = 1.0/bm_r
real, parameter :: nc0 = 3.E8
has_qc = .false.
has_qi = .false.
has_qs = .false.
do k = kts, kte
rqc(k) = max(R1, qc(k)*rho(k))
if (rqc(k).gt.R1) has_qc = .true.
rqi(k) = max(R1, qi(k)*rho(k))
temp = (rho(k)*max(qi(k),qmin))
temp = sqrt(sqrt(temp*temp*temp))
ni(k) = min(max(5.38e7*temp,1.e3),1.e6)
rni(k)= max(R2, ni(k)*rho(k))
if (rqi(k).gt.R1 .and. rni(k).gt.R2) has_qi = .true.
rqs(k) = max(R1, qs(k)*rho(k))
if (rqs(k).gt.R1) has_qs = .true.
enddo
if (has_qc) then
do k=kts,kte
if (rqc(k).le.R1) CYCLE
lamdac = (pidnc*nc0/rqc(k))**obmr
re_qc(k) = max(2.51E-6,min(1.5*(1.0/lamdac),50.E-6))
enddo
endif
if (has_qi) then
do k=kts,kte
if (rqi(k).le.R1 .or. rni(k).le.R2) CYCLE
diai = 11.9*sqrt(rqi(k)/ni(k))
re_qi(k) = max(10.01E-6,min(0.75*0.163*diai,125.E-6))
enddo
endif
if (has_qs) then
do k=kts,kte
if (rqs(k).le.R1) CYCLE
supcol = t0c-t(k)
n0sfac = max(min(exp(alpha*supcol),n0smax/n0s),1.)
lamdas = sqrt(sqrt(pidn0s*n0sfac/rqs(k)))
re_qs(k) = max(25.E-6,min(0.5*(1./lamdas), 999.E-6))
enddo
endif
end subroutine effectRad_wsm6
END MODULE module_mp_wsm6