MODULE module_mp_wsm3 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 :: lamdagmax = 6.e4 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 :: qcrmin = 1.e-9 REAL, SAVE :: & qc0, qck1, pidnc, & bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, & g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, & precr1,precr2,xmmax,roqimax,bvts1, & bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, & g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, & pidn0s,xlv1,pi, & rslopermax,rslopesmax,rslopegmax, & rsloperbmax,rslopesbmax,rslopegbmax, & rsloper2max,rslopes2max,rslopeg2max, & rsloper3max,rslopes3max,rslopeg3max CONTAINS SUBROUTINE wsm3(th, q, qci, qrs & , w, 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 & , 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, & qci, & qrs REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & INTENT(IN ) :: w, & 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 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, & INTENT(INOUT) :: snow, & snowncv, & 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( its:ite , kts:kte ) :: t INTEGER :: i,j,k REAL, DIMENSION( kts:kte ) :: t1d REAL, DIMENSION( kts:kte ) :: den1d REAL, DIMENSION( kts:kte ) :: qc1d REAL, DIMENSION( kts:kte ) :: qi1d REAL, DIMENSION( kts:kte ) :: qs1d 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) ENDDO ENDDO CALL wsm32D(t, q(ims,kms,j), qci(ims,kms,j) & ,qrs(ims,kms,j),w(ims,kms,j), 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) & ,snow(ims,j),snowncv(ims,j) & ,sr(ims,j) & ,ids,ide, jds,jde, kds,kde & ,ims,ime, jms,jme, kms,kme & ,its,ite, jts,jte, kts,kte & ) DO K=kts,kte DO I=its,ite th(i,k,j)=t(i,k)/pii(i,k,j) ENDDO ENDDO 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) if(t(i,k).ge.t0c) then qc1d(k) = qci(i,k,j) qi1d(k) = 0.0 qs1d(k) = 0.0 else qc1d(k) = 0.0 qi1d(k) = qci(i,k,j) qs1d(k) = qrs(i,k,j) endif enddo call effectRad_wsm3(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 wsm3 SUBROUTINE wsm32D(t, q & ,qci, qrs,w, den, p, delz & ,delt,g, cpd, cpv, rd, rv, t0c & ,ep1, ep2, qmin & ,XLS, XLV0, XLF0, den0, denr & ,cliq,cice,psat & ,lat & ,rain, rainncv & ,snow,snowncv & ,sr & ,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, & lat REAL, DIMENSION( its:ite , kts:kte ), & INTENT(INOUT) :: & t REAL, DIMENSION( ims:ime , kms:kme ), & INTENT(INOUT) :: & q, & qci, & qrs REAL, DIMENSION( ims:ime , kms:kme ), & INTENT(IN ) :: w, & 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 REAL, DIMENSION( ims:ime ), OPTIONAL, & INTENT(INOUT) :: snow, & snowncv, & sr REAL, DIMENSION( its:ite , kts:kte ) :: & rh, & qs, & denfac, & rslope, & rslope2, & rslope3, & qrs_tmp, & den_tmp, & delz_tmp, & rslopeb REAL, DIMENSION( its:ite , kts:kte ) :: & pgen, & pisd, & paut, & pacr, & pres, & pcon REAL, DIMENSION( its:ite , kts:kte ) :: & fall, & falk, & xl, & cpm, & work1, & work2, & xni, & qs0, & denqci, & denqrs, & n0sfac, & falkc, & work1c, & work2c, & fallc REAL, DIMENSION( its:ite ) :: delqrs,& delqi REAL, DIMENSION(its:ite) :: tstepsnow INTEGER, DIMENSION( its:ite ) :: kwork1,& kwork2 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, & fallsum, fallsum_qsi, vt2i,vt2s,acrfac, & qdt, pvt, qik, delq, facq, qrsci, frzmlt, & snomlt, hold, holdrs, facqci, supcol, coeres, & supsat, dtcld, xmi, qciik, delqci, eacrs, satdt, & qimax, diameter, xni0, roqi0, supice,holdc, holdci INTEGER :: i, j, k, mstepmax, & iprt, latd, lond, loop, loops, ifsat, kk, n, idim, kdim REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp REAL, DIMENSION( its:ite ) :: tvec1 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) = max(qci(i,k),0.0) qrs(i,k) = max(qrs(i,k),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) = 0. sr(i) = 0. tstepsnow(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 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 cvap = cpv hvap=xlv0 hsub=xls 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) if(t(i,k).lt.ttp) then qs(i,k) =psat*(exp(log(tr)*(xai)))*exp(xbi*(1.-tr)) else qs(i,k) =psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr)) endif qs0(i,k) =psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr)) qs0(i,k) = (qs0(i,k)-qs(i,k))/qs(i,k) qs(i,k) = min(qs(i,k),0.99*p(i,k)) qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k)) qs(i,k) = max(qs(i,k),qmin) rh(i,k) = max(q(i,k) / qs(i,k),qmin) enddo enddo do k = kts, kte do i = its, ite pres(i,k) = 0. paut(i,k) = 0. pacr(i,k) = 0. pgen(i,k) = 0. pisd(i,k) = 0. pcon(i,k) = 0. fall(i,k) = 0. falk(i,k) = 0. fallc(i,k) = 0. falkc(i,k) = 0. xni(i,k) = 1.e3 enddo enddo do k = kts, kte do i = its, ite xni(i,k) = min(max(5.38e7 & *exp(log((den(i,k)*max(qci(i,k),qmin)))*(0.75)),1.e3),1.e6) enddo enddo do k = kts, kte do i = its, ite qrs_tmp(i,k) = qrs(i,k) enddo enddo call slope_wsm3(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, & work1,its,ite,kts,kte) do k = kte, kts, -1 do i = its, ite denqrs(i,k) = den(i,k)*qrs(i,k) enddo enddo call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1,denqrs, & delqrs,dtcld,1,1) do k = kts, kte do i = its, ite qrs(i,k) = max(denqrs(i,k)/den(i,k),0.) fall(i,k) = denqrs(i,k)*work1(i,k)/delz(i,k) enddo enddo do i = its, ite fall(i,1) = delqrs(i)/delz(i,1)/dtcld enddo do k = kte, kts, -1 do i = its, ite if(t(i,k).lt.t0c.and.qci(i,k).gt.0.) then xmi = den(i,k)*qci(i,k)/xni(i,k) diameter = max(dicon * sqrt(xmi), 1.e-25) work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31)) else work1c(i,k) = 0. endif enddo enddo do k = kte, kts, -1 do i = its, ite denqci(i,k) = den(i,k)*qci(i,k) 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) = 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 mstep(i) = 0 enddo do k = kts, kte do i = its, ite if(t(i,k).ge.t0c) then mstep(i) = k endif enddo enddo do i = its, ite kwork2(i) = mstep(i) kwork1(i) = mstep(i) if(mstep(i).ne.0) then if (w(i,mstep(i)).gt.0.) then kwork1(i) = mstep(i) + 1 endif endif enddo do i = its, ite k = kwork1(i) kk = kwork2(i) if(k*kk.ge.1) then qrsci = qrs(i,k) + qci(i,k) if(qrsci.gt.0..or.fall(i,kk).gt.0.) then frzmlt = min(max(-w(i,k)*qrsci/delz(i,k),-qrsci/dtcld), & qrsci/dtcld) snomlt = min(max(fall(i,kk)/den(i,kk),-qrs(i,k)/dtcld), & qrs(i,k)/dtcld) if(k.eq.kk) then t(i,k) = t(i,k) - xlf0/cpm(i,k)*(frzmlt+snomlt)*dtcld else t(i,k) = t(i,k) - xlf0/cpm(i,k)*frzmlt*dtcld t(i,kk) = t(i,kk) - xlf0/cpm(i,kk)*snomlt*dtcld endif endif endif enddo do i = its, ite fallsum = fall(i,1) fallsum_qsi = 0. if((t0c-t(i,1)).gt.0) then fallsum = fallsum+fallc(i,1) fallsum_qsi = fall(i,1)+fallc(i,1) endif if(fallsum.gt.0.) then rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000. + rainncv(i) rain(i) = fallsum*delz(i,1)/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) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snowncv(i) snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i) ENDIF endif IF ( PRESENT (snowncv) ) THEN if(fallsum.gt.0.) sr(i) = snowncv(i)/(rainncv(i)+1.e-12) ELSE if(fallsum.gt.0.) sr(i) = tstepsnow(i)/(rainncv(i)+1.e-12) ENDIF enddo do k = kts, kte do i = its, ite qrs_tmp(i,k) = qrs(i,k) enddo enddo call slope_wsm3(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, & work1,its,ite,kts,kte) do k = kts, kte do i = its, ite if(t(i,k).ge.t0c) then work1(i,k) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k)) else work1(i,k) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k)) endif 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) satdt = supsat/dtcld if(t(i,k).ge.t0c) then if(qci(i,k).gt.qc0) then paut(i,k) = qck1*exp(log(qci(i,k))*((7./3.))) paut(i,k) = min(paut(i,k),qci(i,k)/dtcld) endif if(qrs(i,k).gt.qcrmin.and.qci(i,k).gt.qmin) then pacr(i,k) = min(pacrr*rslope3(i,k)*rslopeb(i,k) & *qci(i,k)*denfac(i,k),qci(i,k)/dtcld) endif if(qrs(i,k).gt.0.) then coeres = rslope2(i,k)*sqrt(rslope(i,k)*rslopeb(i,k)) pres(i,k) = (rh(i,k)-1.)*(precr1*rslope2(i,k) & +precr2*work2(i,k)*coeres)/work1(i,k) if(pres(i,k).lt.0.) then pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld) pres(i,k) = max(pres(i,k),satdt/2) else pres(i,k) = min(pres(i,k),satdt/2) endif endif else supcol = t0c-t(i,k) n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.) ifsat = 0 xni(i,k) = min(max(5.38e7 & *exp(log((den(i,k)*max(qci(i,k),qmin)))*(0.75)),1.e3),1.e6) eacrs = exp(0.07*(-supcol)) if(qrs(i,k).gt.qcrmin.and.qci(i,k).gt.qmin) then xmi = den(i,k)*qci(i,k)/xni(i,k) diameter = min(dicon * sqrt(xmi),dimax) vt2i = 1.49e4*diameter**1.31 vt2s = pvts*rslopeb(i,k)*denfac(i,k) acrfac = 2.*rslope3(i,k)+2.*diameter*rslope2(i,k) & +diameter**2*rslope(i,k) pacr(i,k) = min(pi*qci(i,k)*eacrs*n0s*n0sfac(i,k) & *abs(vt2s-vt2i)*acrfac/4.,qci(i,k)/dtcld) endif if(qci(i,k).gt.0.) then xmi = den(i,k)*qci(i,k)/xni(i,k) diameter = dicon * sqrt(xmi) pisd(i,k) = 4.*diameter*xni(i,k)*(rh(i,k)-1.)/work1(i,k) if(pisd(i,k).lt.0.) then pisd(i,k) = max(pisd(i,k),satdt/2) pisd(i,k) = max(pisd(i,k),-qci(i,k)/dtcld) else pisd(i,k) = min(pisd(i,k),satdt/2) endif if(abs(pisd(i,k)).ge.abs(satdt)) ifsat = 1 endif if(qrs(i,k).gt.0..and.ifsat.ne.1) then coeres = rslope2(i,k)*sqrt(rslope(i,k)*rslopeb(i,k)) pres(i,k) = (rh(i,k)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k) & +precs2*work2(i,k)*coeres)/work1(i,k) supice = satdt-pisd(i,k) if(pres(i,k).lt.0.) then pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld) pres(i,k) = max(max(pres(i,k),satdt/2),supice) else pres(i,k) = min(min(pres(i,k),satdt/2),supice) endif if(abs(pisd(i,k)+pres(i,k)).ge.abs(satdt)) ifsat = 1 endif if(supsat.gt.0.and.ifsat.ne.1) then supice = satdt-pisd(i,k)-pres(i,k) xni0 = 1.e3*exp(0.1*supcol) roqi0 = 4.92e-11*exp(log(xni0)*(1.33)) pgen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k),0.))/dtcld) pgen(i,k) = min(min(pgen(i,k),satdt),supice) endif if(qci(i,k).gt.0.) then qimax = roqimax/den(i,k) paut(i,k) = max(0.,(qci(i,k)-qimax)/dtcld) endif endif enddo enddo do k = kts, kte do i = its, ite qciik = max(qmin,qci(i,k)) delqci = (paut(i,k)+pacr(i,k)-pgen(i,k)-pisd(i,k))*dtcld if(delqci.ge.qciik) then facqci = qciik/delqci paut(i,k) = paut(i,k)*facqci pacr(i,k) = pacr(i,k)*facqci pgen(i,k) = pgen(i,k)*facqci pisd(i,k) = pisd(i,k)*facqci endif qik = max(qmin,q(i,k)) delq = (pres(i,k)+pgen(i,k)+pisd(i,k))*dtcld if(delq.ge.qik) then facq = qik/delq pres(i,k) = pres(i,k)*facq pgen(i,k) = pgen(i,k)*facq pisd(i,k) = pisd(i,k)*facq endif work2(i,k) = -pres(i,k)-pgen(i,k)-pisd(i,k) q(i,k) = q(i,k)+work2(i,k)*dtcld qci(i,k) = max(qci(i,k)-(paut(i,k)+pacr(i,k)-pgen(i,k)-pisd(i,k)) & *dtcld,0.) qrs(i,k) = max(qrs(i,k)+(paut(i,k)+pacr(i,k)+pres(i,k))*dtcld,0.) if(t(i,k).lt.t0c) then t(i,k) = t(i,k)-xls*work2(i,k)/cpm(i,k)*dtcld else t(i,k) = t(i,k)-xl(i,k)*work2(i,k)/cpm(i,k)*dtcld endif enddo enddo cvap = cpv hvap = xlv0 hsub = xls 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)=psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr)) qs(i,k) = min(qs(i,k),0.99*p(i,k)) qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k)) qs(i,k) = max(qs(i,k),qmin) denfac(i,k) = sqrt(den0/den(i,k)) enddo enddo do k = kts, kte do i = its, ite work1(i,k) = conden(t(i,k),q(i,k),qs(i,k),xl(i,k),cpm(i,k)) work2(i,k) = qci(i,k)+work1(i,k) pcon(i,k) = min(max(work1(i,k),0.),max(q(i,k),0.))/dtcld if(qci(i,k).gt.0..and.work1(i,k).lt.0.and.t(i,k).gt.t0c) & pcon(i,k) = max(work1(i,k),-qci(i,k))/dtcld q(i,k) = q(i,k)-pcon(i,k)*dtcld qci(i,k) = max(qci(i,k)+pcon(i,k)*dtcld,0.) t(i,k) = t(i,k)+pcon(i,k)*xl(i,k)/cpm(i,k)*dtcld enddo enddo do k = kts, kte do i = its, ite if(qci(i,k).le.qmin) qci(i,k) = 0.0 if(qrs(i,k).le.qcrmin) qrs(i,k) = 0.0 enddo enddo enddo END SUBROUTINE wsm32D 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 wsm3init(den0,denr,dens,cl,cpv,allowed_to_read) IMPLICIT NONE REAL, INTENT(IN) :: den0,denr,dens,cl,cpv LOGICAL, INTENT(IN) :: allowed_to_read 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 g1pbr = rgmma(bvtr1) g3pbr = rgmma(bvtr3) g4pbr = rgmma(bvtr4) 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 xmmax = (dimax/dicon)**2 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 rslopermax = 1./lamdarmax rslopesmax = 1./lamdasmax rsloperbmax = rslopermax ** bvtr rslopesbmax = rslopesmax ** bvts rsloper2max = rslopermax * rslopermax rslopes2max = rslopesmax * rslopesmax rsloper3max = rsloper2max * rslopermax rslopes3max = rslopes2max * rslopesmax END SUBROUTINE wsm3init subroutine slope_wsm3(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, & den, & denfac, & t, & rslope, & rslopeb, & rslope2, & rslope3, & vt REAL, PARAMETER :: t0c = 273.15 REAL, DIMENSION( its:ite , kts:kte ) :: & n0sfac REAL :: lamdar,lamdas,x, y, z, supcol, pvt integer :: i, j, k lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) do k = kts, kte do i = its, ite if(t(i,k).ge.t0c) then pvt = pvtr 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) = exp(log(rslope(i,k))*(bvtr)) rslope2(i,k) = rslope(i,k)*rslope(i,k) rslope3(i,k) = rslope2(i,k)*rslope(i,k) endif else supcol = t0c-t(i,k) n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.) pvt = pvts 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) = exp(log(rslope(i,k))*(bvts)) rslope2(i,k) = rslope(i,k)*rslope(i,k) rslope3(i,k) = rslope2(i,k)*rslope(i,k) endif endif vt(i,k) = pvt*rslopeb(i,k)*denfac(i,k) if(qrs(i,k).le.0.0) vt(i,k) = 0.0 enddo enddo END subroutine slope_wsm3 SUBROUTINE nislfv_rain_pcm(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 zsumt,qsumt,zsumb,qsumb 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) do k=2,km wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k)) enddo 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_wsm3(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km) if( n.eq.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 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 if( kt-kb.eq.1 ) then qn(k) = qa(kb) else if( kt-kb.ge.2 ) then zsumb = za(kb+1)-zi(k) qsumb = qa(kb) * zsumb zsumt = zi(k+1)-za(kt-1) qsumt = qa(kt-1) * zsumt qsum = 0.0 zsum = 0.0 if( kt-kb.ge.3 ) then do m=kb+1,kt-2 qsum = qsum + qa(m) * dza(m) zsum = zsum + dza(m) enddo endif qn(k) = (qsumb+qsum+qsumt)/(zsumb+zsum+zsumt) 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_pcm 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_wsm3(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 effectRad_wsm3 (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_wsm3 END MODULE module_mp_wsm3