442 integer ,
intent(in ) :: &
443 ids,ide, jds,jde, kds,kde, &
444 ims,ime, jms,jme, kms,kme, &
445 its,ite, jts,jte, kts,kte
446 integer ,
intent(in ) :: lat
447 integer ,
optional,
intent(in ) :: microphysics_debug, diag_i_dbg, diag_j_dbg
448 real(kind=kind_phys) ,
intent(in ) :: delt
449 real(kind=kind_phys) ,
intent(in ) :: g, cpd, cpv, t0c, &
455 real(kind=kind_phys) ,
intent(in ) :: ccn0
456 real(kind=kind_phys),
dimension(ims:ime) ,
intent(in ) :: slmsk
457 real(kind=kind_phys),
dimension(ims:ime,kms:kme) ,
intent(in ) :: p
458 real(kind=kind_phys),
dimension(ims:ime,kms:kme) ,
intent(in ) :: delz
459 real(kind=kind_phys),
dimension(ims:ime,kms:kme) ,
intent(in ) :: den
460 real(kind=kind_phys),
dimension(its:ite,kts:kte) ,
intent(inout) :: t
461 real(kind=kind_phys),
dimension(its:ite,kts:kte,2) ,
intent(inout) :: qci
462 real(kind=kind_phys),
dimension(its:ite,kts:kte,3) ,
intent(inout) :: qrs
463 real(kind=kind_phys),
dimension(its:ite,kts:kte,3) ,
intent(inout) :: ncr
464 real(kind=kind_phys),
dimension(ims:ime,kms:kme) ,
intent(inout) :: q
465 real(kind=kind_phys),
dimension(ims:ime) ,
intent(inout) :: rain
466 real(kind=kind_phys),
dimension(ims:ime) ,
intent(inout) :: rainncv
467 real(kind=kind_phys),
dimension(ims:ime) ,
intent(inout) :: sr
468 real(kind=kind_phys),
dimension(ims:ime),
optional ,
intent(inout) :: snow
469 real(kind=kind_phys),
dimension(ims:ime),
optional ,
intent(inout) :: snowncv
470 real(kind=kind_phys),
dimension(ims:ime),
optional ,
intent(inout) :: graupel
471 real(kind=kind_phys),
dimension(ims:ime),
optional ,
intent(inout) :: graupelncv
474 integer :: debug_local, i_dbg_local, j_dbg_local
476 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: dend
477 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: qcr
478 real(kind=kind_phys),
dimension(its:ite,kts:kte,3) :: rh
479 real(kind=kind_phys),
dimension(its:ite,kts:kte,3) :: qs
480 real(kind=kind_phys),
dimension(its:ite,kts:kte,3) :: rslope, rslope2, rslope3, rslopeb
481 real(kind=kind_phys),
dimension(its:ite,kts:kte,3) :: falk, fall
482 real(kind=kind_phys),
dimension(its:ite,kts:kte,3) :: work1
483 real(kind=kind_phys),
dimension(its:ite,kts:kte,3) :: qrs_tmp
484 real(kind=kind_phys),
dimension(its:ite,kts:kte,2) :: avedia
485 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: rslopec, rslopec2,rslopec3
486 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: workn, falln, falkn
487 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: worka, workr
488 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: den_tmp, delz_tmp, ncr_tmp
489 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: lamdr_tmp
490 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: lamdc_tmp
491 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: falkc, work1c, work2c, fallc
492 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: dqr,dnr
493 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: pcact, prevp, psdep, pgdep, praut, &
494 psaut, pgaut, pracw, psacw, pgacw, &
495 pgacr, pgacs, psaci, pgmlt, praci, &
496 piacr, pracs, psacr, pgaci, pseml, &
497 pgeml, paacw, pigen, pidep, pcond, &
499 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: nraut, nracw, ncevp, nccol, nrcol, &
500 nsacw, ngacw, niacr, nsacr, ngacr, &
501 naacw, nseml, ngeml, ncact
502 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: xl, cpm, work2, denfac, xni, n0sfac, &
504 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: denqrs1, denqrs2, denqrs3
505 real(kind=kind_phys),
dimension(its:ite,kts:kte) :: denqr1, denncr3, denqci
506 real(kind=kind_phys),
dimension(its:ite) :: delqrs1, delqrs2, delqrs3, delncr3
507 real(kind=kind_phys),
dimension(its:ite) :: delqi
508 real(kind=kind_phys),
dimension(its:ite) :: tstepsnow, tstepgraup
509 real(kind=kind_phys) :: gfac, sfac
517 real(kind=kind_phys),
dimension(its:ite) :: tvec1
518 integer,
dimension(its:ite) :: mnstep, numndt
519 integer,
dimension(its:ite) :: mstep, numdt
520 logical,
dimension(its:ite) :: flgcld
521 real(kind=kind_phys) :: temp
522 real(kind=kind_phys) :: vt2ave
523 real(kind=kind_phys) :: holdc, holdci
524 real(kind=kind_phys) :: cpmcal, xlcal, lamdac, diffus, &
525 viscos, xka, venfac, conden, diffac, &
526 x, y, z, a, b, c, d, e, &
527 ndt, qdt, holdrr, holdrs, holdrg, &
528 supcol, supcolt, pvt, coeres, supsat, &
529 dtcld, xmi, eacrs, satdt, qimax, &
530 diameter, xni0, roqi0, &
531 fallsum, fallsum_qsi, fallsum_qg, &
532 vt2i, vt2r, vt2s, vt2g, acrfac, &
534 xlwork2, factor, source,
value, &
536 taucon, lencon, lenconcr, &
537 xlf, pfrzdtc, pfrzdtr, supice, &
538 alpha2, delta2, delta3
540 integer :: i, j, k, mstepmax, &
541 iprt, latd, lond, loop, loops, ifsat, &
546 real(kind=kind_phys) :: dldti, xb, xai, tr, xbi, xa, hvap, &
547 cvap, hsub, dldt, ttp
552 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
553 xlcal(x) = xlv0-xlv1*(x-t0c)
560 lamdac(x,y,z)= exp(log(((pidnc*z)/(x*y)))*((.33333333)))
565 diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y
566 viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y
567 xka(x,y) = 1.414e3*viscos(x,y)*y
568 diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
569 venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
570 /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
571 conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
579 if (
present(microphysics_debug))
then
580 debug_local = microphysics_debug
584 if (
present(diag_i_dbg))
then
585 i_dbg_local = diag_i_dbg
589 if (
present(diag_j_dbg))
then
590 j_dbg_local = diag_j_dbg
597 qci(i,k,1) = max(qci(i,k,1),0.0)
598 qrs(i,k,1) = max(qrs(i,k,1),0.0)
599 qci(i,k,2) = max(qci(i,k,2),0.0)
600 qrs(i,k,2) = max(qrs(i,k,2),0.0)
601 qrs(i,k,3) = max(qrs(i,k,3),0.0)
602 ncr(i,k,1) = min(max(ncr(i,k,1),1.e8),2.e10)
603 ncr(i,k,2) = max(ncr(i,k,2),0.0)
604 ncr(i,k,3) = max(ncr(i,k,3),0.0)
620 cpm(i,k) = cpmcal(q(i,k))
621 xl(i,k) = xlcal(t(i,k))
627 if(slmsk(i).eq.2)
then
636 delz_tmp(i,k) = delz(i,k)
637 den_tmp(i,k) = den(i,k)
645 if(
present(snowncv) .and.
present(snow)) snowncv(i) = 0.
646 if(
present (graupelncv) .and.
present(graupel)) graupelncv(i) = 0.
657 loops = max(nint(delt/dtcldcr),1)
659 if(delt.le.dtcldcr) dtcld = delt
672 call vrec_d( tvec1(its), den(its,k), ite-its+1)
674 tvec1(i) = tvec1(i)*den0
676 call vsqrt_d( denfac(its,k), tvec1(its), ite-its+1)
689 xb = xa+hvap/(rv*ttp)
692 xbi = xai+hsub/(rv*ttp)
700 qs(i,k,1) = psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
701 qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
703 qs(i,k,1) = ep2*qs(i,k,1)/(p(i,k)-qs(i,k,1))
705 qs(i,k,1) = max(qs(i,k,1),qmin)
707 rh(i,k,1) = max(q(i,k)/qs(i,k,1),qmin)
711 if(t(i,k).lt.ttp)
then
713 qs(i,k,2) = psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
716 qs(i,k,2) = psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
718 qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
720 qs(i,k,2) = ep2*qs(i,k,2)/(p(i,k)-qs(i,k,2))
722 qs(i,k,2) = max(qs(i,k,2),qmin)
724 rh(i,k,2) = max(q(i,k)/qs(i,k,2),qmin)
796 if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin )
then
797 rslopec(i,k) = rslopecmax
798 rslopec2(i,k) = rslopec2max
799 rslopec3(i,k) = rslopec3max
801 rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
802 rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
803 rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
808 temp = (den(i,k)*max(qci(i,k,2),qmin))
809 temp = sqrt(sqrt(temp*temp*temp))
810 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
823 qrs_tmp(i,k,1) = qrs(i,k,1)
824 qrs_tmp(i,k,2) = qrs(i,k,2)
825 qrs_tmp(i,k,3) = qrs(i,k,3)
826 ncr_tmp(i,k) = ncr(i,k,3)
829 call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
830 rslope3,work1,workn,its,ite,kts,kte)
842 work1(i,k,1) = work1(i,k,1)/delz(i,k)
843 workn(i,k) = workn(i,k)/delz(i,k)
844 numdt(i) = max(nint(max(work1(i,k,1),workn(i,k))*dtcld+.5),1)
845 if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
849 if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
859 if(n.le.mstep(i))
then
860 falk(i,k,1) = dend(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
861 falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
862 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
863 falln(i,k) = falln(i,k)+falkn(i,k)
864 qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/dend(i,k),0.)
865 ncr(i,k,3) = max(ncr(i,k,3)-falkn(i,k)*dtcld,0.)
869 do k = kte_in-1,kts,-1
871 if(n.le.mstep(i))
then
872 falk(i,k,1) = dend(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
873 falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
874 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
875 falln(i,k) = falln(i,k)+falkn(i,k)
876 dqr(i,k) = min(falk(i,k,1)*dtcld/dend(i,k),qrs(i,k,1))
877 dqr(i,k+1) = min(falk(i,k+1,1)*delz(i,k+1)/delz(i,k) &
878 *dtcld/dend(i,k),qrs(i,k+1,1))
879 dnr(i,k) = min(falkn(i,k)*dtcld,ncr(i,k,3))
880 dnr(i,k+1) = min(falkn(i,k+1)*delz(i,k+1)/delz(i,k)*dtcld, &
882 qrs(i,k,1) = max(qrs(i,k,1)-dqr(i,k)+dqr(i,k+1),0.)
883 ncr(i,k,3) = max(ncr(i,k,3)-dnr(i,k)+dnr(i,k+1),0.)
889 qrs_tmp(i,k,1) = qrs(i,k,1)
890 ncr_tmp(i,k) = ncr(i,k,3)
893 call slope_rain(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
894 rslope3,work1,workn,its,ite,kts,kte)
897 work1(i,k,1) = work1(i,k,1)/delz(i,k)
898 workn(i,k) = workn(i,k)/delz(i,k)
910 qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.e-15)
911 if(qsum(i,k) .gt. 1.e-15 )
then
912 worka(i,k) = (work1(i,k,2)*qrs(i,k,2) + work1(i,k,3)*qrs(i,k,3)) &
917 denqrs2(i,k) = den(i,k)*qrs(i,k,2)
918 denqrs3(i,k) = den(i,k)*qrs(i,k,3)
921 call nislfv_rain_plm6(idim,kdim,den_tmp,denfac,t,delz_tmp,worka, &
922 denqrs2,denqrs3,delqrs2,delqrs3,dtcld,1,1)
925 qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
926 qrs(i,k,3) = max(denqrs3(i,k)/den(i,k),0.)
927 fall(i,k,2) = denqrs2(i,k)*worka(i,k)/delz(i,k)
928 fall(i,k,3) = denqrs3(i,k)*worka(i,k)/delz(i,k)
932 fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
933 fall(i,1,3) = delqrs3(i)/delz(i,1)/dtcld
942 qrs_tmp(i,k,1) = qrs(i,k,1)
943 qrs_tmp(i,k,2) = qrs(i,k,2)
944 qrs_tmp(i,k,3) = qrs(i,k,3)
945 ncr_tmp(i,k) = ncr(i,k,3)
948 call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
949 rslope3,work1,workn,its,ite,kts,kte)
958 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
959 if(t(i,k).gt.t0c)
then
965 work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
966 if(qrs(i,k,2).gt.0.)
then
967 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
968 psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. &
969 *n0sfac(i,k)*(precs1*rslope2(i,k,2) &
970 +precs2*work2(i,k)*coeres)/den(i,k)
971 psmlt(i,k) = min(max(psmlt(i,k)*dtcld,-qrs(i,k,2)),0.)
976 if(qrs(i,k,2).gt.qcrmin)
then
977 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2)
978 ncr(i,k,3) = ncr(i,k,3) - sfac*psmlt(i,k)
981 qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
982 qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
983 t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
989 if(qrs(i,k,3).gt.0.)
then
990 coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
991 pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*(precg1 &
992 *rslope2(i,k,3) + precg2*work2(i,k)*coeres) &
994 pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld,-qrs(i,k,3)),0.)
999 if(qrs(i,k,3).gt.qcrmin)
then
1000 gfac = rslope(i,k,3)*n0g/qrs(i,k,3)
1001 ncr(i,k,3) = ncr(i,k,3) - gfac*pgmlt(i,k)
1004 qrs(i,k,3) = qrs(i,k,3) + pgmlt(i,k)
1005 qrs(i,k,1) = qrs(i,k,1) - pgmlt(i,k)
1006 t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
1018 if(qci(i,k,2).le.0.)
then
1021 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
1022 diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
1023 work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
1034 denqci(i,k) = den(i,k)*qci(i,k,2)
1037 call nislfv_rain_plmr(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci,denqci, &
1041 qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
1045 fallc(i,1) = delqi(i)/delz(i,1)/dtcld
1056 fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
1057 fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
1058 fallsum_qg = fall(i,kts,3)
1059 if(fallsum.gt.0.)
then
1060 rainncv(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rainncv(i)
1061 rain(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rain(i)
1063 If(fallsum_qsi.gt.0.)
then
1064 tstepsnow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + tstepsnow(i)
1065 IF(
PRESENT (snowncv) .AND.
PRESENT (snow))
THEN
1066 snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snowncv(i)
1067 snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i)
1070 IF(fallsum_qg.gt.0.)
then
1071 tstepgraup(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. &
1073 IF(
PRESENT (graupelncv) .AND.
PRESENT (graupel))
THEN
1074 graupelncv(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. &
1076 graupel(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. + graupel(i)
1079 IF (
PRESENT (snowncv))
THEN
1080 if(fallsum.gt.0.)sr(i)=(snowncv(i) + graupelncv(i))/(rainncv(i)+1.e-12)
1082 if(fallsum.gt.0.)sr(i)=(tstepsnow(i) + tstepgraup(i))/(rainncv(i)+1.e-12)
1098 if(supcol.lt.0.) xlf = xlf0
1099 if(supcol.lt.0 .and. qci(i,k,2).gt.0.)
then
1100 qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
1105 if(qci(i,k,2).gt.qmin)
then
1106 ncr(i,k,2) = ncr(i,k,2) + xni(i,k)
1108 t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
1116 if(supcol.gt.40. .and. qci(i,k,1).gt.0.)
then
1117 qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
1122 if(ncr(i,k,2).gt.0.) ncr(i,k,2) = 0.
1123 t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
1131 if(supcol.gt.0. .and. qci(i,k,1).gt.qmin)
then
1132 supcolt=min(supcol,70.)
1133 pfrzdtc = min(pi*pi*pfrz1*(exp(pfrz2*supcolt)-1.)*denr/den(i,k) &
1134 *ncr(i,k,2)*rslopec3(i,k)*rslopec3(i,k)/18.*dtcld &
1140 if(ncr(i,k,2).gt.ncmin)
then
1141 nfrzdtc = min(pi*pfrz1*(exp(pfrz2*supcolt)-1.)*ncr(i,k,2) &
1142 *rslopec3(i,k)/6.*dtcld,ncr(i,k,2))
1143 ncr(i,k,2) = ncr(i,k,2) - nfrzdtc
1145 qci(i,k,2) = qci(i,k,2) + pfrzdtc
1146 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
1147 qci(i,k,1) = qci(i,k,1)-pfrzdtc
1154 if(supcol.gt.0. .and. qrs(i,k,1).gt.0.)
then
1155 supcolt=min(supcol,70.)
1156 pfrzdtr = min(140.*(pi*pi)*pfrz1*ncr(i,k,3)*denr/den(i,k) &
1157 *(exp(pfrz2*supcolt)-1.)*rslope3(i,k,1)*rslope3(i,k,1) &
1165 if(ncr(i,k,3).gt.nrmin)
then
1166 nfrzdtr = min(4.*pi*pfrz1*ncr(i,k,3)*(exp(pfrz2*supcolt)-1.) &
1167 *rslope3(i,k,1)*dtcld, ncr(i,k,3))
1168 ncr(i,k,3) = ncr(i,k,3) - nfrzdtr
1170 qrs(i,k,3) = qrs(i,k,3) + pfrzdtr
1171 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
1172 qrs(i,k,1) = qrs(i,k,1) - pfrzdtr
1183 ncr(i,k,2) = max(ncr(i,k,2),0.0)
1184 ncr(i,k,3) = max(ncr(i,k,3),0.0)
1194 qrs_tmp(i,k,1) = qrs(i,k,1)
1195 qrs_tmp(i,k,2) = qrs(i,k,2)
1196 qrs_tmp(i,k,3) = qrs(i,k,3)
1197 ncr_tmp(i,k) = ncr(i,k,3)
1200 call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
1201 rslope3,work1,workn,its,ite,kts,kte)
1208 avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333))
1210 if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin)
then
1211 rslopec(i,k) = rslopecmax
1212 rslopec2(i,k) = rslopec2max
1213 rslopec3(i,k) = rslopec3max
1215 rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
1216 rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
1217 rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
1223 avedia(i,k,1) = rslopec(i,k)
1232 work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
1235 work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
1238 work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
1256 supsat = max(q(i,k),qmin)-qs(i,k,1)
1257 satdt = supsat/dtcld
1262 lencon = 2.7e-2*den(i,k)*qci(i,k,1)*(1.e20/16.*rslopec2(i,k) &
1264 lenconcr = max(1.2*lencon, qcrmin)
1265 if(qci(i,k,1).gt.qcr(i,k).and.ncr(i,k,2).gt.ncmin)
then
1266 praut(i,k) = qck1*qci(i,k,1)**(7./3.)*ncr(i,k,2)**(-1./3.)
1267 praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
1272 nraut(i,k) = 3.5e9*den(i,k)*praut(i,k)
1273 if(qrs(i,k,1).gt.lenconcr) &
1274 nraut(i,k) = ncr(i,k,3)/qrs(i,k,1)*praut(i,k)
1275 nraut(i,k) = min(nraut(i,k),ncr(i,k,2)/dtcld)
1283 if(qrs(i,k,1).ge.lenconcr)
then
1284 if(avedia(i,k,2).ge.di100)
then
1285 nracw(i,k) = min(ncrk1*ncr(i,k,2)*ncr(i,k,3)*(rslopec3(i,k) &
1286 + 24.*rslope3(i,k,1)),ncr(i,k,2)/dtcld)
1287 pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk1*ncr(i,k,2) &
1288 *ncr(i,k,3)*rslopec3(i,k)*(2.*rslopec3(i,k) &
1289 + 24.*rslope3(i,k,1)),qci(i,k,1)/dtcld)
1291 nracw(i,k) = min(ncrk2*ncr(i,k,2)*ncr(i,k,3)*(2.*rslopec3(i,k) &
1292 *rslopec3(i,k)+5040.*rslope3(i,k,1) &
1293 *rslope3(i,k,1)),ncr(i,k,2)/dtcld)
1294 pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk2*ncr(i,k,2) &
1295 *ncr(i,k,3)*rslopec3(i,k)*(6.*rslopec3(i,k) &
1296 *rslopec3(i,k)+5040.*rslope3(i,k,1)*rslope3(i,k,1)) &
1304 if(avedia(i,k,1).ge.di100)
then
1305 nccol(i,k) = ncrk1*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k)
1307 nccol(i,k) = 2.*ncrk2*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k) &
1314 if(qrs(i,k,1).ge.lenconcr)
then
1315 if(avedia(i,k,2).lt.di100)
then
1316 nrcol(i,k) = 5040.*ncrk2*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1) &
1318 elseif(avedia(i,k,2).ge.di100 .and. avedia(i,k,2).lt.di600)
then
1319 nrcol(i,k) = 24.*ncrk1*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1)
1320 elseif(avedia(i,k,2).ge.di600 .and. avedia(i,k,2).lt.di2000)
then
1321 coecol = -2.5e3*(avedia(i,k,2)-di600)
1322 nrcol(i,k) = 24.*exp(coecol)*ncrk1*ncr(i,k,3)*ncr(i,k,3) &
1332 if(qrs(i,k,1).gt.0.)
then
1333 coeres = rslope(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
1334 prevp(i,k) = (rh(i,k,1)-1.)*ncr(i,k,3)*(precr1*rslope(i,k,1) &
1335 + precr2*work2(i,k)*coeres)/work1(i,k,1)
1336 if(prevp(i,k).lt.0.)
then
1337 prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
1338 prevp(i,k) = max(prevp(i,k),satdt/2)
1343 if(prevp(i,k).eq.-qrs(i,k,1)/dtcld)
then
1344 ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,3)
1347 else if(prevp(i,k).eq.0.)
then
1354 prevp(i,k) = min(prevp(i,k),satdt/2)
1373 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1374 supsat = max(q(i,k),qmin)-qs(i,k,2)
1375 satdt = supsat/dtcld
1382 temp = (den(i,k)*max(qci(i,k,2),qmin))
1383 temp = sqrt(sqrt(temp*temp*temp))
1384 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
1385 eacrs = exp(0.07*(-supcol))
1387 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
1388 diameter = min(dicon * sqrt(xmi),dimax)
1389 vt2i = 1.49e4*diameter**1.31
1390 vt2r=pvtr*rslopeb(i,k,1)*denfac(i,k)
1391 vt2s=pvts*rslopeb(i,k,2)*denfac(i,k)
1392 vt2g=pvtg*rslopeb(i,k,3)*denfac(i,k)
1393 qsum(i,k) = max((qrs(i,k,2)+qrs(i,k,3)),1.e-15)
1394 if(qsum(i,k) .gt. 1.e-15)
then
1395 vt2ave=(vt2s*qrs(i,k,2)+vt2g*qrs(i,k,3))/(qsum(i,k))
1399 if(supcol.gt.0. .and. qci(i,k,2).gt.qmin)
then
1400 if(qrs(i,k,1).gt.qcrmin)
then
1405 acrfac = 6.*rslope2(i,k,1)+4.*diameter*rslope(i,k,1) + diameter**2
1406 praci(i,k) = pi*qci(i,k,2)*ncr(i,k,3)*abs(vt2r-vt2i)*acrfac/4.
1408 praci(i,k) = praci(i,k)*min(max(0.0,qrs(i,k,1)/qci(i,k,2)),1.)**2
1409 praci(i,k) = min(praci(i,k),qci(i,k,2)/dtcld)
1414 piacr(i,k) = pi*pi*avtr*ncr(i,k,3)*denr*xni(i,k)*denfac(i,k) &
1415 *g7pbr*rslope3(i,k,1)*rslope2(i,k,1)*rslopeb(i,k,1) &
1418 piacr(i,k) = piacr(i,k)*min(max(0.0,qci(i,k,2)/qrs(i,k,1)),1.)**2
1419 piacr(i,k) = min(piacr(i,k),qrs(i,k,1)/dtcld)
1425 if(ncr(i,k,3).gt.nrmin)
then
1426 niacr(i,k) = pi*avtr*ncr(i,k,3)*xni(i,k)*denfac(i,k)*g4pbr &
1427 *rslope2(i,k,1)*rslopeb(i,k,1)/4.
1429 niacr(i,k) = niacr(i,k)*min(max(0.0,qci(i,k,2)/qrs(i,k,1)),1.)**2
1430 niacr(i,k) = min(niacr(i,k),ncr(i,k,3)/dtcld)
1436 if(qrs(i,k,2).gt.qcrmin)
then
1437 acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
1438 + diameter**2*rslope(i,k,2)
1439 psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k) &
1440 *abs(vt2ave-vt2i)*acrfac/4.
1441 psaci(i,k) = min(psaci(i,k),qci(i,k,2)/dtcld)
1449 if(qrs(i,k,3).gt.qcrmin)
then
1450 egi = exp(0.07*(-supcol))
1451 acrfac = 2.*rslope3(i,k,3)+2.*diameter*rslope2(i,k,3) &
1452 + diameter**2*rslope(i,k,3)
1453 pgaci(i,k) = pi*egi*qci(i,k,2)*n0g*abs(vt2ave-vt2i)*acrfac/4.
1454 pgaci(i,k) = min(pgaci(i,k),qci(i,k,2)/dtcld)
1465 if(qrs(i,k,2).gt.qcrmin .and. qci(i,k,1).gt.qmin)
then
1466 psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) &
1468 *min(max(0.0,qrs(i,k,2)/qci(i,k,1)),1.)**2 &
1469 *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
1475 if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,2).gt.ncmin)
then
1476 nsacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) &
1478 *min(max(0.0,qrs(i,k,2)/qci(i,k,1)),1.)**2 &
1479 *ncr(i,k,2)*denfac(i,k),ncr(i,k,2)/dtcld)
1485 if(qrs(i,k,3).gt.qcrmin .and. qci(i,k,1).gt.qmin)
then
1486 pgacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3)*qci(i,k,1) &
1488 *min(max(0.0,qrs(i,k,3)/qci(i,k,1)),1.)**2 &
1489 *denfac(i,k),qci(i,k,1)/dtcld)
1495 if(qrs(i,k,3).gt.qcrmin .and. ncr(i,k,2).gt.ncmin)
then
1496 ngacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3)*ncr(i,k,2) &
1498 *min(max(0.0,qrs(i,k,3)/qci(i,k,1)),1.)**2 &
1499 *denfac(i,k),ncr(i,k,2)/dtcld)
1505 if(qsum(i,k) .gt. 1.e-15 )
then
1506 paacw(i,k) = (qrs(i,k,2)*psacw(i,k)+qrs(i,k,3)*pgacw(i,k))/(qsum(i,k))
1511 naacw(i,k) = (qrs(i,k,2)*nsacw(i,k)+qrs(i,k,3)*ngacw(i,k))/(qsum(i,k))
1518 if(qrs(i,k,2).gt.qcrmin .and. qrs(i,k,1).gt.qcrmin)
then
1519 if(supcol.gt.0)
then
1520 acrfac = 5.*rslope3(i,k,2)*rslope3(i,k,2) &
1521 + 4.*rslope3(i,k,2)*rslope2(i,k,2)*rslope(i,k,1) &
1522 + 1.5*rslope2(i,k,2)*rslope2(i,k,2)*rslope2(i,k,1)
1523 pracs(i,k) = pi*pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2r-vt2ave) &
1524 *(dens/den(i,k))*acrfac
1526 pracs(i,k) = pracs(i,k)*min(max(0.0,qrs(i,k,1)/qrs(i,k,2)),1.)**2
1527 pracs(i,k) = min(pracs(i,k),qrs(i,k,2)/dtcld)
1533 acrfac = 30.*rslope3(i,k,1)*rslope2(i,k,1)*rslope(i,k,2) &
1534 +10.*rslope2(i,k,1)*rslope2(i,k,1)*rslope2(i,k,2) &
1535 + 2.*rslope3(i,k,1)*rslope3(i,k,2)
1536 psacr(i,k) = pi*pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2ave-vt2r) &
1537 *(denr/den(i,k))*acrfac
1539 psacr(i,k) = psacr(i,k)*min(max(0.0,qrs(i,k,2)/qrs(i,k,1)),1.)**2
1540 psacr(i,k) = min(psacr(i,k),qrs(i,k,1)/dtcld)
1542 if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,3).gt.nrmin)
then
1547 acrfac = 1.5*rslope2(i,k,1)*rslope(i,k,2) &
1548 + 1.0*rslope(i,k,1)*rslope2(i,k,2)+.5*rslope3(i,k,2)
1549 nsacr(i,k) = pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2ave-vt2r) &
1552 nsacr(i,k) = nsacr(i,k)*min(max(0.0,qrs(i,k,2)/qrs(i,k,1)),1.)**2
1553 nsacr(i,k) = min(nsacr(i,k),ncr(i,k,3)/dtcld)
1559 if(qrs(i,k,3).gt.qcrmin .and. qrs(i,k,1).gt.qcrmin)
then
1560 acrfac = 30.*rslope3(i,k,1)*rslope2(i,k,1)*rslope(i,k,3) &
1561 +10.*rslope2(i,k,1)*rslope2(i,k,1)*rslope2(i,k,3) &
1562 + 2.*rslope3(i,k,1)*rslope3(i,k,3)
1563 pgacr(i,k) = pi*pi*ncr(i,k,3)*n0g*abs(vt2ave-vt2r)*(denr/den(i,k)) &
1566 pgacr(i,k) = pgacr(i,k)*min(max(0.0,qrs(i,k,3)/qrs(i,k,1)),1.)**2
1567 pgacr(i,k) = min(pgacr(i,k),qrs(i,k,1)/dtcld)
1573 if(qrs(i,k,3).gt.qcrmin .and. ncr(i,k,3).gt.nrmin)
then
1574 acrfac = 1.5*rslope2(i,k,1)*rslope(i,k,3) &
1575 + 1.0*rslope(i,k,1)*rslope2(i,k,3) + .5*rslope3(i,k,3)
1576 ngacr(i,k) = pi*ncr(i,k,3)*n0g*abs(vt2ave-vt2r)*acrfac
1578 ngacr(i,k) = ngacr(i,k)*min(max(0.0,qrs(i,k,3)/qrs(i,k,1)),1.)**2
1579 ngacr(i,k) = min(ngacr(i,k),ncr(i,k,3)/dtcld)
1587 if(qrs(i,k,3).gt.qcrmin .and. qrs(i,k,2).gt.qcrmin)
then
1590 if(supcol.le.0)
then
1596 if(qrs(i,k,2).gt.0.) &
1597 pseml(i,k) = min(max(cliq*supcol*(paacw(i,k)+psacr(i,k)) &
1598 /xlf,-qrs(i,k,2)/dtcld),0.)
1603 if (qrs(i,k,2).gt.qcrmin)
then
1604 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2)
1605 nseml(i,k) = -sfac*pseml(i,k)
1611 if(qrs(i,k,3).gt.0.) &
1612 pgeml(i,k) = min(max(cliq*supcol*(paacw(i,k)+pgacr(i,k))/xlf &
1613 ,-qrs(i,k,3)/dtcld),0.)
1618 if (qrs(i,k,3).gt.qcrmin)
then
1619 gfac = rslope(i,k,3)*n0g/qrs(i,k,3)
1620 ngeml(i,k) = -gfac*pgeml(i,k)
1625 if(supcol.gt.0)
then
1630 if(qci(i,k,2).gt.0. .and. ifsat.ne.1)
then
1631 pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1632 supice = satdt-prevp(i,k)
1633 if(pidep(i,k).lt.0.)
then
1634 pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1635 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
1637 pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1639 if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1646 if(qrs(i,k,2).gt.0. .and. ifsat.ne.1)
then
1647 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1648 psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2) &
1649 + precs2*work2(i,k)*coeres)/work1(i,k,2)
1650 supice = satdt-prevp(i,k)-pidep(i,k)
1651 if(psdep(i,k).lt.0.)
then
1652 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
1653 psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1655 psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1657 if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) ifsat = 1
1664 if(qrs(i,k,3).gt.0. .and. ifsat.ne.1)
then
1665 coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
1666 pgdep(i,k) = (rh(i,k,2)-1.)*(precg1*rslope2(i,k,3) &
1667 + precg2*work2(i,k)*coeres)/work1(i,k,2)
1668 supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1669 if(pgdep(i,k).lt.0.)
then
1670 pgdep(i,k) = max(pgdep(i,k),-qrs(i,k,3)/dtcld)
1671 pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
1673 pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
1675 if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge. &
1676 abs(satdt)) ifsat = 1
1682 if(supsat.gt.0. .and. ifsat.ne.1)
then
1683 supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
1684 xni0 = 1.e3*exp(0.1*supcol)
1685 roqi0 = 4.92e-11*xni0**1.33
1687 pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))/dtcld)
1689 pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1698 if(qci(i,k,2).gt.0.)
then
1699 qimax = roqimax/den(i,k)
1700 psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1707 if(qrs(i,k,2).gt.0.)
then
1708 alpha2 = 1.e-3*exp(0.09*(-supcol))
1709 pgaut(i,k) = min(max(0.,alpha2*(qrs(i,k,2)-qs0)),qrs(i,k,2)/dtcld)
1719 if(supcol.lt.0.)
then
1720 if(qrs(i,k,2).gt.0. .and. rh(i,k,1).lt.1.)
then
1721 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1722 psevp(i,k) = (rh(i,k,1)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2) &
1723 +precs2*work2(i,k)*coeres)/work1(i,k,1)
1724 psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1730 if(qrs(i,k,3).gt.0. .and. rh(i,k,1).lt.1.)
then
1731 coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
1732 pgevp(i,k) = (rh(i,k,1)-1.)*(precg1*rslope2(i,k,3) &
1733 + precg2*work2(i,k)*coeres)/work1(i,k,1)
1734 pgevp(i,k) = min(max(pgevp(i,k),-qrs(i,k,3)/dtcld),0.)
1752 if(qrs(i,k,1).lt.1.e-4 .and. qrs(i,k,2).lt.1.e-4) delta2=1.
1753 if(qrs(i,k,1).lt.1.e-4) delta3=1.
1754 if(t(i,k).le.t0c)
then
1758 value = max(qmin,qci(i,k,1))
1759 source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k)) &
1761 if (source.gt.
value)
then
1762 factor =
value/source
1763 praut(i,k) = praut(i,k)*factor
1764 pracw(i,k) = pracw(i,k)*factor
1765 paacw(i,k) = paacw(i,k)*factor
1770 value = max(qmin,qci(i,k,2))
1771 source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k)+psaci(i,k) &
1773 if (source.gt.
value)
then
1774 factor =
value/source
1775 psaut(i,k) = psaut(i,k)*factor
1776 pigen(i,k) = pigen(i,k)*factor
1777 pidep(i,k) = pidep(i,k)*factor
1778 praci(i,k) = praci(i,k)*factor
1779 psaci(i,k) = psaci(i,k)*factor
1780 pgaci(i,k) = pgaci(i,k)*factor
1785 value = max(qmin,qrs(i,k,1))
1786 source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k) &
1787 +psacr(i,k)+pgacr(i,k))*dtcld
1788 if (source.gt.
value)
then
1789 factor =
value/source
1790 praut(i,k) = praut(i,k)*factor
1791 prevp(i,k) = prevp(i,k)*factor
1792 pracw(i,k) = pracw(i,k)*factor
1793 piacr(i,k) = piacr(i,k)*factor
1794 psacr(i,k) = psacr(i,k)*factor
1795 pgacr(i,k) = pgacr(i,k)*factor
1800 value = max(qmin,qrs(i,k,2))
1801 source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k) &
1802 +piacr(i,k)*delta3+praci(i,k)*delta3 &
1803 -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2 &
1804 +psaci(i,k)-pgacs(i,k) )*dtcld
1805 if (source.gt.
value)
then
1806 factor =
value/source
1807 psdep(i,k) = psdep(i,k)*factor
1808 psaut(i,k) = psaut(i,k)*factor
1809 pgaut(i,k) = pgaut(i,k)*factor
1810 paacw(i,k) = paacw(i,k)*factor
1811 piacr(i,k) = piacr(i,k)*factor
1812 praci(i,k) = praci(i,k)*factor
1813 psaci(i,k) = psaci(i,k)*factor
1814 pracs(i,k) = pracs(i,k)*factor
1815 psacr(i,k) = psacr(i,k)*factor
1816 pgacs(i,k) = pgacs(i,k)*factor
1821 value = max(qmin,qrs(i,k,3))
1822 source = -(pgdep(i,k)+pgaut(i,k) &
1823 +piacr(i,k)*(1.-delta3)+praci(i,k)*(1.-delta3) &
1824 +psacr(i,k)*(1.-delta2)+pracs(i,k)*(1.-delta2) &
1825 +pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
1826 if (source.gt.
value)
then
1827 factor =
value/source
1828 pgdep(i,k) = pgdep(i,k)*factor
1829 pgaut(i,k) = pgaut(i,k)*factor
1830 piacr(i,k) = piacr(i,k)*factor
1831 praci(i,k) = praci(i,k)*factor
1832 psacr(i,k) = psacr(i,k)*factor
1833 pracs(i,k) = pracs(i,k)*factor
1834 paacw(i,k) = paacw(i,k)*factor
1835 pgaci(i,k) = pgaci(i,k)*factor
1836 pgacr(i,k) = pgacr(i,k)*factor
1837 pgacs(i,k) = pgacs(i,k)*factor
1842 value = max(ncmin,ncr(i,k,2))
1843 source = (nraut(i,k)+nccol(i,k)+nracw(i,k) &
1844 +naacw(i,k)+naacw(i,k))*dtcld
1845 if (source.gt.
value)
then
1846 factor =
value/source
1847 nraut(i,k) = nraut(i,k)*factor
1848 nccol(i,k) = nccol(i,k)*factor
1849 nracw(i,k) = nracw(i,k)*factor
1850 naacw(i,k) = naacw(i,k)*factor
1855 value = max(nrmin,ncr(i,k,3))
1856 source = (-nraut(i,k)+nrcol(i,k)+niacr(i,k)+nsacr(i,k)+ngacr(i,k) &
1858 if (source.gt.
value)
then
1859 factor =
value/source
1860 nraut(i,k) = nraut(i,k)*factor
1861 nrcol(i,k) = nrcol(i,k)*factor
1862 niacr(i,k) = niacr(i,k)*factor
1863 nsacr(i,k) = nsacr(i,k)*factor
1864 ngacr(i,k) = ngacr(i,k)*factor
1867 work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k)+pidep(i,k))
1870 q(i,k) = q(i,k)+work2(i,k)*dtcld
1871 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1872 +paacw(i,k)+paacw(i,k))*dtcld,0.)
1873 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1874 +prevp(i,k)-piacr(i,k)-pgacr(i,k) &
1875 -psacr(i,k))*dtcld,0.)
1876 qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+praci(i,k) &
1877 +psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k)) &
1879 qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+paacw(i,k) &
1880 -pgaut(i,k)+piacr(i,k)*delta3 &
1881 +praci(i,k)*delta3+psaci(i,k)-pgacs(i,k) &
1882 -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2) &
1884 qrs(i,k,3) = max(qrs(i,k,3)+(pgdep(i,k)+pgaut(i,k) &
1885 +piacr(i,k)*(1.-delta3) &
1886 +praci(i,k)*(1.-delta3)+psacr(i,k)*(1.-delta2) &
1887 +pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k) &
1888 +pgacr(i,k)+pgacs(i,k))*dtcld,0.)
1889 ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k) &
1890 -naacw(i,k)-naacw(i,k))*dtcld,0.)
1891 ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)-niacr(i,k) &
1892 -nsacr(i,k)-ngacr(i,k))*dtcld,0.)
1894 xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k)) &
1895 -xl(i,k)*prevp(i,k)-xlf*(piacr(i,k)+paacw(i,k) &
1896 +paacw(i,k)+pgacr(i,k)+psacr(i,k))
1898 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1903 value = max(qmin,qci(i,k,1))
1904 source= (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k)) &
1906 if (source.gt.
value)
then
1907 factor =
value/source
1908 praut(i,k) = praut(i,k)*factor
1909 pracw(i,k) = pracw(i,k)*factor
1910 paacw(i,k) = paacw(i,k)*factor
1915 value = max(qmin,qrs(i,k,1))
1916 source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k) &
1917 -pracw(i,k)-paacw(i,k)-prevp(i,k))*dtcld
1918 if (source.gt.
value)
then
1919 factor =
value/source
1920 praut(i,k) = praut(i,k)*factor
1921 prevp(i,k) = prevp(i,k)*factor
1922 pracw(i,k) = pracw(i,k)*factor
1923 paacw(i,k) = paacw(i,k)*factor
1924 pseml(i,k) = pseml(i,k)*factor
1925 pgeml(i,k) = pgeml(i,k)*factor
1930 value = max(qcrmin,qrs(i,k,2))
1931 source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
1932 if (source.gt.
value)
then
1933 factor =
value/source
1934 pgacs(i,k) = pgacs(i,k)*factor
1935 psevp(i,k) = psevp(i,k)*factor
1936 pseml(i,k) = pseml(i,k)*factor
1941 value = max(qcrmin,qrs(i,k,3))
1942 source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
1943 if (source.gt.
value)
then
1944 factor =
value/source
1945 pgacs(i,k) = pgacs(i,k)*factor
1946 pgevp(i,k) = pgevp(i,k)*factor
1947 pgeml(i,k) = pgeml(i,k)*factor
1952 value = max(ncmin,ncr(i,k,2))
1953 source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+naacw(i,k) &
1955 if (source.gt.
value)
then
1956 factor =
value/source
1957 nraut(i,k) = nraut(i,k)*factor
1958 nccol(i,k) = nccol(i,k)*factor
1959 nracw(i,k) = nracw(i,k)*factor
1960 naacw(i,k) = naacw(i,k)*factor
1965 value = max(nrmin,ncr(i,k,3))
1966 source = (-nraut(i,k)+nrcol(i,k)-nseml(i,k)-ngeml(i,k) &
1968 if (source.gt.
value)
then
1969 factor =
value/source
1970 nraut(i,k) = nraut(i,k)*factor
1971 nrcol(i,k) = nrcol(i,k)*factor
1972 nseml(i,k) = nseml(i,k)*factor
1973 ngeml(i,k) = ngeml(i,k)*factor
1976 work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
1978 q(i,k) = q(i,k)+work2(i,k)*dtcld
1979 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1980 +paacw(i,k)+paacw(i,k))*dtcld,0.)
1981 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1982 +prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k) &
1983 -pgeml(i,k))*dtcld,0.)
1984 qrs(i,k,2) = max(qrs(i,k,2)+(psevp(i,k)-pgacs(i,k) &
1985 +pseml(i,k))*dtcld,0.)
1986 qrs(i,k,3) = max(qrs(i,k,3)+(pgacs(i,k)+pgevp(i,k) &
1987 +pgeml(i,k))*dtcld,0.)
1988 ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k) &
1989 -naacw(i,k)-naacw(i,k))*dtcld,0.)
1990 ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)+nseml(i,k) &
1991 +ngeml(i,k))*dtcld,0.)
1993 xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k)) &
1994 -xlf*(pseml(i,k)+pgeml(i,k))
1995 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
2015 xbi=xai+hsub/(rv*ttp)
2019 qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
2020 qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
2021 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
2022 qs(i,k,1) = max(qs(i,k,1),qmin)
2024 if(t(i,k).lt.ttp)
then
2025 qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
2027 qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
2029 qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
2030 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
2031 qs(i,k,2) = max(qs(i,k,2),qmin)
2032 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
2042 qrs_tmp(i,k,1) = qrs(i,k,1)
2043 qrs_tmp(i,k,2) = qrs(i,k,2)
2044 qrs_tmp(i,k,3) = qrs(i,k,3)
2045 ncr_tmp(i,k) = ncr(i,k,3)
2049 call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
2050 rslope3,work1,workn,its,ite,kts,kte)
2057 avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333))
2062 if(avedia(i,k,2).le.di82)
then
2063 ncr(i,k,2) = ncr(i,k,2)+ncr(i,k,3)
2069 qci(i,k,1) = qci(i,k,1)+qrs(i,k,1)
2084 if(rh(i,k,1).gt.1.)
then
2086 ncact(i,k) = max(0.,((ncr(i,k,1)+ncr(i,k,2)) &
2087 *min(1.,(rh(i,k,1)/satmax)**actk) - ncr(i,k,2)))/dtcld
2088 ncact(i,k) =min(ncact(i,k),max(ncr(i,k,1),0.)/dtcld)
2089 pcact(i,k) = min(4.*pi*denr*(actr*1.e-6)**3*ncact(i,k)/ &
2090 (3.*den(i,k)),max(q(i,k),0.)/dtcld)
2091 q(i,k) = max(q(i,k)-pcact(i,k)*dtcld,0.)
2092 qci(i,k,1) = max(qci(i,k,1)+pcact(i,k)*dtcld,0.)
2093 ncr(i,k,1) = max(ncr(i,k,1)-ncact(i,k)*dtcld,0.)
2094 ncr(i,k,2) = max(ncr(i,k,2)+ncact(i,k)*dtcld,0.)
2095 t(i,k) = t(i,k)+pcact(i,k)*xl(i,k)/cpm(i,k)*dtcld
2106 qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
2107 qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
2108 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
2109 qs(i,k,1) = max(qs(i,k,1),qmin)
2110 work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
2111 work2(i,k) = qci(i,k,1)+work1(i,k,1)
2112 pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
2113 if(qci(i,k,1).gt.0. .and. work1(i,k,1).lt.0.) &
2114 pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
2119 if(pcond(i,k).eq.-qci(i,k,1)/dtcld)
then
2121 ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,2)
2126 q(i,k) = q(i,k)-pcond(i,k)*dtcld
2127 qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
2128 t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
2141 if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
2142 if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
2143 if(qrs(i,k,1).ge.qcrmin .and. ncr(i,k,3) .ge. nrmin)
then
2144 lamdr_tmp(i,k) = exp(log(((pidnr*ncr(i,k,3)) &
2145 /(den(i,k)*qrs(i,k,1))))*((.33333333)))
2146 if(lamdr_tmp(i,k) .le. lamdarmin)
then
2147 lamdr_tmp(i,k) = lamdarmin
2148 ncr(i,k,3) = den(i,k)*qrs(i,k,1)*lamdr_tmp(i,k)**3/pidnr
2149 elseif(lamdr_tmp(i,k) .ge. lamdarmax)
then
2150 lamdr_tmp(i,k) = lamdarmax
2151 ncr(i,k,3) = den(i,k)*qrs(i,k,1)*lamdr_tmp(i,k)**3/pidnr
2154 if(qci(i,k,1).ge.qmin .and. ncr(i,k,2) .ge. ncmin )
then
2155 lamdc_tmp(i,k) = exp(log(((pidnc*ncr(i,k,2)) &
2156 /(den(i,k)*qci(i,k,1))))*((.33333333)))
2157 if(lamdc_tmp(i,k) .le. lamdacmin)
then
2158 lamdc_tmp(i,k) = lamdacmin
2159 ncr(i,k,2) = den(i,k)*qci(i,k,1)*lamdc_tmp(i,k)**3/pidnc
2160 elseif(lamdc_tmp(i,k) .ge. lamdacmax)
then
2161 lamdc_tmp(i,k) = lamdacmax
2162 ncr(i,k,2) = den(i,k)*qci(i,k,1)*lamdc_tmp(i,k)**3/pidnc