510 & qn, qtend, cwtend, qc, qi, nc, ni,fprcp,qrn,qsnw,nrn,nsnw, &
511 & p, pdel, cldn, liqcldf, &
512 & icecldf, cldo, pint, rpdel, zm, rate1ord_cw2pr_st, naai, &
513 & npccnin, rndst,nacon, rhdfda, rhu00, fice, tlat, qvlat, qctend, &
514 & qitend, nctend, nitend, effc, effc_fn, effi, prect, preci, &
515 & nevapr, evapsnow, prain, prodsnow, cmeout, deffi, pgamrad, &
516 & lamcrad,qsout2,dsout2,qrout2, drout2, qcsevap,qisevap,qvres, &
517 & cmeiout, vtrmc,vtrmi,qcsedten,qisedten, prao, &
518 & prco,mnuccco,mnuccto, &
519 & msacwio,psacwso, bergso,bergo,melto,homoo,qcreso,prcio,praio, &
520 & qireso, mnuccro,pracso,meltsdt,frzrdt, ncal, ncai, mnuccdo, &
522 & nsout2, nrout2, ncnst, ninst, nimm, miu_disp, nsoot, rnsoot, &
523 & ui_scale, dcrit, nnuccdo, nnuccco, nsacwio, nsubio, nprcio, &
524 & npraio, npccno, npsacwso, nsubco, nprao, nprc1o, tlataux, &
525 & nbincontactdust, lprint, xlat, xlon, rhc)
535 use constituents,
only: pcnst
536 real(r8),
parameter :: ncnst = 100.0e6
537 real(r8),
parameter :: ninst = 0.10e6
542 integer,
intent (in) :: pcols, pver,fprcp
543 real(r8),
intent (in) :: ncnst
544 real(r8),
intent (in) :: ninst
545 real(r8),
intent (in) :: nimm (pcols,pver)
546 real(r8),
intent (in) :: miu_disp , ui_scale, dcrit
548 real(r8),
intent (in) :: nsoot (pcols,pver) , rnsoot (pcols,pver)
556 integer,
intent(in) :: lchnk
557 integer,
intent(in) :: ncol
558 real(r8),
intent(in) :: deltatin
559 real(r8),
intent(in) :: tn(pcols,pver)
560 real(r8),
intent(in) :: ttend(pcols,pver)
561 real(r8),
intent(in) :: qn(pcols,pver)
562 real(r8),
intent(in) :: qtend(pcols,pver)
563 real(r8),
intent(in) :: cwtend(pcols,pver)
565 real(r8),
intent(inout) :: qc(pcols,pver)
566 real(r8),
intent(inout) :: qi(pcols,pver)
567 real(r8),
intent(inout) :: nc(pcols,pver)
568 real(r8),
intent(inout) :: ni(pcols,pver)
569 real(r8),
intent(inout) :: qrn(pcols,pver)
570 real(r8),
intent(inout) :: qsnw(pcols,pver)
571 real(r8),
intent(inout) :: nrn(pcols,pver)
572 real(r8),
intent(inout) :: nsnw(pcols,pver)
573 real(r8),
intent(in) :: p(pcols,pver)
574 real(r8),
intent(in) :: pdel(pcols,pver)
575 real(r8),
intent(in) :: cldn(pcols,pver)
576 real(r8),
intent(in) :: icecldf(pcols,pver)
577 real(r8),
intent(in) :: liqcldf(pcols,pver)
578 real(r8),
intent(inout) :: cldo(pcols,pver)
579 real(r8),
intent(in) :: pint(pcols,pver+1)
581 real(r8),
intent(in) :: rpdel(pcols,pver)
582 real(r8),
intent(in) :: zm(pcols,pver)
584 real(r8),
intent(in) :: rhc(pcols,pver)
586 real(r8),
intent(out) :: rate1ord_cw2pr_st(pcols,pver)
589 real(r8),
intent(in) :: naai(pcols,pver)
590 real(r8),
intent(inout) :: npccnin(pcols,pver)
591 integer :: nbincontactdust
592 real(r8),
intent(in),
dimension(pcols,pver, 10) :: rndst, nacon
595 real(r8),
intent(in) :: rhdfda(pcols,pver)
596 real(r8),
intent(in) :: rhu00(pcols,pver)
597 real(r8),
intent(in) :: fice(pcols,pver)
599 real(r8),
intent(out) :: tlat(pcols,pver)
601 real(r8),
intent(out) :: tlataux(pcols,pver)
603 real(r8),
intent(out) :: qvlat(pcols,pver)
604 real(r8),
intent(out) :: qctend(pcols,pver)
605 real(r8),
intent(out) :: qitend(pcols,pver)
606 real(r8),
intent(out) :: nctend(pcols,pver)
607 real(r8),
intent(out) :: nitend(pcols,pver)
608 real(r8),
intent(out) :: effc(pcols,pver)
609 real(r8),
intent(out) :: effc_fn(pcols,pver)
610 real(r8),
intent(out) :: effi(pcols,pver)
611 real(r8),
intent(out) :: prect(pcols)
612 real(r8),
intent(out) :: preci(pcols)
613 real(r8),
intent(out) :: nevapr(pcols,pver)
614 real(r8),
intent(out) :: evapsnow(pcols,pver)
615 real(r8),
intent(out) :: prain(pcols,pver)
616 real(r8),
intent(out) :: prodsnow(pcols,pver)
617 real(r8),
intent(out) :: cmeout(pcols,pver)
618 real(r8),
intent(out) :: deffi(pcols,pver)
619 real(r8),
intent(out) :: pgamrad(pcols,pver)
620 real(r8),
intent(out) :: lamcrad(pcols,pver)
622 real(r8),
intent(out) :: qcsevap(pcols,pver)
623 real(r8),
intent(out) :: qisevap(pcols,pver)
624 real(r8),
intent(out) :: qvres(pcols,pver)
625 real(r8),
intent(out) :: cmeiout(pcols,pver)
626 real(r8),
intent(out) :: vtrmc(pcols,pver)
627 real(r8),
intent(out) :: vtrmi(pcols,pver)
628 real(r8),
intent(out) :: qcsedten(pcols,pver)
629 real(r8),
intent(out) :: qisedten(pcols,pver)
631 real(r8),
intent(out) :: prao(pcols,pver)
632 real(r8),
intent(out) :: prco(pcols,pver)
633 real(r8),
intent(out) :: mnuccco(pcols,pver)
634 real(r8),
intent(out) :: mnuccto(pcols,pver)
635 real(r8),
intent(out) :: msacwio(pcols,pver)
636 real(r8),
intent(out) :: psacwso(pcols,pver)
637 real(r8),
intent(out) :: bergso(pcols,pver)
638 real(r8),
intent(out) :: bergo(pcols,pver)
639 real(r8),
intent(out) :: melto(pcols,pver)
640 real(r8),
intent(out) :: homoo(pcols,pver)
641 real(r8),
intent(out) :: qcreso(pcols,pver)
642 real(r8),
intent(out) :: prcio(pcols,pver)
643 real(r8),
intent(out) :: praio(pcols,pver)
644 real(r8),
intent(out) :: qireso(pcols,pver)
645 real(r8),
intent(out) :: mnuccro(pcols,pver)
646 real(r8),
intent(out) :: pracso (pcols,pver)
647 real(r8),
intent(out) :: meltsdt(pcols,pver)
648 real(r8),
intent(out) :: frzrdt (pcols,pver)
649 real(r8),
intent(out) :: mnuccdo(pcols,pver)
650 real(r8),
intent(out) :: nnuccto(pcols,pver)
652 real(r8),
intent(out) :: nnuccdo(pcols,pver)
653 real(r8),
intent(out) :: nnuccco(pcols,pver)
654 real(r8),
intent(out) :: nsacwio(pcols,pver)
655 real(r8),
intent(out) :: nsubio(pcols,pver)
656 real(r8),
intent(out) :: nprcio(pcols,pver)
657 real(r8),
intent(out) :: npraio(pcols,pver)
659 real(r8),
intent(out) :: npccno(pcols,pver)
660 real(r8),
intent(out) :: npsacwso(pcols,pver)
661 real(r8),
intent(out) :: nsubco(pcols,pver)
662 real(r8),
intent(out) :: nprao(pcols,pver)
663 real(r8),
intent(out) :: nprc1o(pcols,pver)
670 real(r8) :: t1(pcols,pver)
671 real(r8) :: q1(pcols,pver)
672 real(r8) :: qc1(pcols,pver)
673 real(r8) :: qi1(pcols,pver)
674 real(r8) :: nc1(pcols,pver)
675 real(r8) :: ni1(pcols,pver)
676 real(r8) :: tlat1(pcols,pver)
677 real(r8) :: qvlat1(pcols,pver)
678 real(r8) :: qctend1(pcols,pver)
679 real(r8) :: qitend1(pcols,pver)
680 real(r8) :: nctend1(pcols,pver)
681 real(r8) :: nitend1(pcols,pver)
682 real(r8) :: prect1(pcols)
683 real(r8) :: preci1(pcols)
685 real(r8) :: tlat1_aux(pcols,pver)
690 real(r8) :: ncold(pcols,pver)
693 real(r8) :: deltat, dti, deltam
697 real(r8),
dimension(pcols,pver) :: q, t, rho, irho, rhof, dv
698 &, mu, sc, kap, cldmax, cldm
699 &, icldm, lcldm, cme, cmei
700 &, cwml, cwmi, lcldn, lcldo
701 &, nctend_mixnuc, npccn
703 real(r8),
dimension(pcols) :: icwc, calpha, cbeta, cbetah
704 &, cgamma, cgamah, rcgama
705 &, cmec1, cmec2, cmec3, cmec4
707 real(r8),
dimension(pver) :: nnuccd, mnuccd
708 &, qcsinksum_rate1ord
711 real(r8) :: qtmp, dum, dum1, dum2, qcld
715 real(r8) :: qcic(pcols,pver)
716 real(r8) :: qiic(pcols,pver)
717 real(r8) :: qniic(pcols,pver)
718 real(r8) :: qric(pcols,pver)
719 real(r8) :: ncic(pcols,pver)
720 real(r8) :: niic(pcols,pver)
721 real(r8) :: nsic(pcols,pver)
722 real(r8) :: nric(pcols,pver)
723 real(r8) :: lami(pver)
724 real(r8) :: n0i(pver)
725 real(r8) :: lamc(pver)
726 real(r8) :: n0c(pver)
727 real(r8) :: lams(pver)
728 real(r8) :: n0s(pver)
729 real(r8) :: lamr(pver)
730 real(r8) :: n0r(pver)
731 real(r8) :: cdist1(pver)
733 real(r8) :: rercld(pcols,pver)
734 real(r8) :: arcld(pcols,pver)
738 real(r8) :: pgam(pver)
742 real(r8) :: mnuccc(pver)
743 real(r8) :: nnuccc(pver)
745 real(r8) :: mnucct(pver)
746 real(r8) :: nnucct(pver)
747 real(r8) :: msacwi(pver)
748 real(r8) :: nsacwi(pver)
750 real(r8) :: prc(pver)
751 real(r8) :: nprc(pver)
752 real(r8) :: nprc1(pver)
753 real(r8) :: nsagg(pver)
757 real(r8) :: psacws(pver)
758 real(r8) :: npsacws(pver)
761 real(r8) :: uns(pver)
762 real(r8) :: ums(pver)
763 real(r8) :: unr(pver)
764 real(r8) :: umr(pver)
767 real(r8) :: pracs(pver)
768 real(r8) :: npracs(pver)
769 real(r8) :: mnuccr(pver)
770 real(r8) :: nnuccr(pver)
771 real(r8) :: pra(pver)
772 real(r8) :: npra(pver)
773 real(r8) :: nragg(pver)
774 real(r8) :: prci(pver)
775 real(r8) :: nprci(pver)
776 real(r8) :: prai(pver)
777 real(r8) :: nprai(pver)
784 real(r8) :: abi,oneoabi
787 real(r8) :: pre(pver)
788 real(r8) :: prds(pver)
794 real(r8) :: dumc(pcols,pver)
795 real(r8) :: dumnc(pcols,pver)
796 real(r8) :: dumi(pcols,pver)
797 real(r8) :: dumni(pcols,pver)
798 real(r8) :: dums(pcols,pver)
799 real(r8) :: dumns(pcols,pver)
800 real(r8) :: dumr(pcols,pver)
801 real(r8) :: dumnr(pcols,pver)
804 real(r8) :: fnr(pver)
806 real(r8) :: fnc(pver)
808 real(r8) :: fni(pver)
810 real(r8) :: fns(pver)
811 real(r8) :: faloutr(pver)
812 real(r8) :: faloutnr(pver)
813 real(r8) :: faloutc(pver)
814 real(r8) :: faloutnc(pver)
815 real(r8) :: falouti(pver)
816 real(r8) :: faloutni(pver)
817 real(r8) :: falouts(pver)
818 real(r8) :: faloutns(pver)
827 real(r8) :: faltndqie
828 real(r8) :: faltndqce
830 real(r8) :: relhum(pcols,pver)
831 real(r8) :: csigma(pcols)
833 real(r8) :: arn(pcols,pver)
834 real(r8) :: asn(pcols,pver)
835 real(r8) :: acn(pcols,pver)
836 real(r8) :: ain(pcols,pver)
837 real(r8) :: nsubi(pver)
838 real(r8) :: nsubc(pver)
839 real(r8) :: nsubs(pver)
840 real(r8) :: nsubr(pver)
842 real(r8) :: dz(pcols,pver)
845 real(r8) :: nfice(pcols,pver)
848 real(r8) :: rflx(pcols,pver+1)
849 real(r8) ::
sflx(pcols,pver+1)
851 real(r8) :: rflx1(pcols,pver+1)
852 real(r8) :: sflx1(pcols,pver+1)
855 real(r8) :: tsp(pcols,pver)
856 real(r8) :: qsp(pcols,pver)
857 real(r8) :: qsphy(pcols,pver)
858 real(r8) :: qs(pcols)
859 real(r8) :: es(pcols)
860 real(r8) :: esl(pcols,pver)
861 real(r8) :: esi(pcols,pver)
866 real(r8) :: qnitend(pcols,pver)
867 real(r8) :: nstend(pcols,pver)
868 real(r8) :: qrtend(pcols,pver)
869 real(r8) :: nrtend(pcols,pver)
883 real(r8) :: berg(pcols,pver)
884 real(r8) :: bergs(pver)
893 real(r8) :: qrout(pcols,pver)
894 real(r8) :: nrout(pcols,pver)
895 real(r8) :: nsout(pcols,pver)
896 real(r8) :: dsout(pcols,pver)
897 real(r8) :: qsout(pcols,pver)
901 real(r8) ,
intent(out) :: qrout2(pcols,pver)
902 real(r8) ,
intent(out) :: qsout2(pcols,pver)
903 real(r8) ,
intent(out) :: nrout2(pcols,pver)
904 real(r8) ,
intent(out) :: nsout2(pcols,pver)
905 real(r8) :: freqs(pcols,pver)
906 real(r8) :: freqr(pcols,pver)
908 real(r8),
intent(out),
dimension(pcols,pver) :: drout2, dsout2
911 real(r8) :: dum2i(pcols,pver)
912 real(r8) :: dum2l(pcols,pver)
913 real(r8) :: ncmax(pcols,pver)
917 real(r8) :: ncai(pcols,pver)
918 real(r8) :: ncal(pcols,pver)
921 integer i,k,nstep,n, l
922 integer ii,kk, m, ind_aux, km, kp
925 integer iter,it,ltrue(pcols)
928 real(r8) tcnt, viscosity, mfp, nslipsoot, ndfaersoot
929 real(r8),
dimension(nbincontactdust) :: ndfaer, nslip, slip
933 real(r8) bbi, cci, ak, iciwc, rvi, riter
936 real(r8) tk, deles, aprpr, bprpr, cice, qi0, crate, qidep
939 real(r8) cldmw(pcols,pver)
951 real(r8) :: refl(pcols,pver)
952 real(r8) :: rainrt(pcols,pver)
953 real(r8) :: rainrt1(pcols,pver)
954 real(r8) :: csrfl(pcols,pver)
955 real(r8) :: arefl(pcols,pver)
956 real(r8) :: acsrfl(pcols,pver)
957 real(r8) :: frefl(pcols,pver)
958 real(r8) :: fcsrfl(pcols,pver)
959 real(r8) :: areflz(pcols,pver)
960 real(r8) :: tmp, miu_ice(pver)
964 real(r8) dmc,ssmc,dstrn
966 real(r8),
parameter :: cdnl = 0.e6_r8
970 integer,
parameter :: auto_option=1
972 real(r8) :: beta6, xs, nssoot, nsdust, taux, psc, bh, vaux, aux,
973 & lw, nw, tx1, tx2, tx3, tx4, tx5, omeps, esloesi
986 logical :: nccons,nicons
1020 qcsedten(i,k) = zero
1021 qisedten(i,k) = zero
1054 npsacwso(i,k) = zero
1064 deltam =
one / max(deltat, 150.0_r8)
1069 omsm = 0.99999999_r8
1070 dto2 = 0.5_r8*deltat
1081 rho(i,k) = p(i,k) / (
r*t(i,k))
1082 irho(i,k) =
one / rho(i,k)
1083 dv(i,k) = 8.794e-5_r8*t(i,k)**1.81_r8/p(i,k)
1084 tx1 = t(i,k) * sqrt(t(i,k)) / (t(i,k)+120._r8)
1085 mu(i,k) = 1.496e-6_r8 * tx1
1086 sc(i,k) = mu(i,k)/(rho(i,k)*dv(i,k))
1087 kap(i,k) = 1.414e3_r8 * 1.496e-6_r8 * tx1
1093 rhof(i,k) = (
rhosu*irho(i,k))**0.54
1095 arn(i,k) =
ar * rhof(i,k)
1096 asn(i,k) =
as * rhof(i,k)
1097 acn(i,k) =
ac * rhof(i,k)
1098 ain(i,k) =
ai * rhof(i,k)
1103 dz(i,k) = pdel(i,k) * irho(i,k) *
ginv
1118 lcldm(i,k) = max(liqcldf(i,k),mincld)
1119 if (qc(i,k) >=
qsmall)
then
1120 npccn(i,k) = max( (lcldm(i,k)*npccnin(i,k)-nc(i,k))
1122 ncold(i,k) = nc(i,k)
1123 nc(i,k) = nc(i,k) + npccn(i,k)*deltat
1138 tlat1_aux(i,k) = zero
1155 evapsnow(i,k) = zero
1157 prodsnow(i,k) = zero
1164 cldmax(i,k) = mincld
1202 esl(i,k) = min(
fpvsl(t(i,k)), p(i,k))
1203 esi(i,k) = min(
fpvsi(t(i,k)), p(i,k))
1205 esl(i,k) = min(
polysvp(t(i,k),0), p(i,k))
1206 esi(i,k) = min(
polysvp(t(i,k),1), p(i,k))
1209 esloesi = esl(i,k) / esi(i,k)
1211 if (t(i,k) >
tmelt) esi(i,k) = esl(i,k)
1213 qs(i) = epsqs*esl(i,k)/(p(i,k)-omeps*esl(i,k))
1215 relhum(i,k) = min(q(i,k)/qs(i),
one)
1223 cldm(i,k) = max(cldn(i,k), mincld)
1224 cldmw(i,k) = max(cldn(i,k), mincld)
1226 icldm(i,k) = max(icecldf(i,k), mincld)
1227 lcldm(i,k) = max(liqcldf(i,k), mincld)
1230 if (qc(i,k) >=
qsmall)
then
1231 tx1 =
one / lcldm(i,k)
1232 dum2l(i,k) = (ncold(i,k)+npccn(i,k)*deltat) * tx1
1233 dum2l(i,k) = max(dum2l(i,k),cdnl*irho(i,k))
1234 ncmax(i,k) = max(dum2l(i,k)*lcldm(i,k), zero)
1235 dum2l(i,k) = npccn(i,k)*deltat*tx1
1245 dumfice = qc(i,k) + qi(i,k)
1247 nfice(i,k) = qi(i,k)/dumfice
1258 dum2 = 0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k)))*1000._r8
1261 dum2 = min(dum2, 208.9e3_r8)*irho(i,k)
1264 dumnnuc = max((dum2*icldm(i,k)-ni(i,k))*deltam, zero)
1269 ninew = ni(i,k) + dumnnuc*deltat
1270 qinew = qi(i,k) + dumnnuc*deltat*
mi0
1292 if (icldm(i,k) > zero)
then
1293 tx1 =
one / icldm(i,k)
1294 qiic(i,k) = qinew * tx1
1295 niic(i,k) = ninew * tx1
1303 niic(i,k) = ninst*irho(i,k)
1309 if (t(i,k) < 273.15)
then
1311 if (qi(i,k) >
qsmall)
then
1315 qvi = epsqs*esi(i,k) / (p(i,k)-omeps*esi(i,k))
1316 qvl = epsqs*esl(i,k) / (p(i,k)-omeps*esl(i,k))
1318 dqsidt =
xxls*qvi / (
rv*t(i,k)*t(i,k))
1324 if (qiic(i,k) >=
qsmall)
then
1329 miu_ice(k) = max(min(0.008_r8*(lami(k)*0.01)**0.87_r8,
1331 tx1 =
one + miu_ice(k)
1334 lami(k) = aux*lami(k)
1335 n0i(k) = niic(i,k) * lami(k)**tx1 * tx2
1340 if (lami(k) <
lammini*aux)
then
1344 niic(i,k) = n0i(k)/lami(k)
1345 else if (lami(k) >
lammaxi*aux)
then
1349 niic(i,k) = n0i(k)/lami(k)
1352 epsi = (miu_ice(k)+
two)*(miu_ice(k)+
one)*pi*
1353 & niic(i,k)*rho(i,k)*dv(i,k)/lami(k)
1356 if (qc(i,k) >
qsmall)
then
1364 prd = epsi*(qvl-qvi)*oneoabi
1372 berg(i,k) = max(zero, prd*min(icldm(i,k),lcldm(i,k)))
1377 if (berg(i,k) > zero)
then
1379 bergtsf = max(zero, tx1/berg(i,k))
1380 if(bergtsf <
one) berg(i,k) = max(zero, tx1)
1385 if (bergtsf <
one .or. icldm(i,k) > lcldm(i,k))
then
1387 if (qiic(i,k) >=
qsmall)
then
1391 if (qc(i,k) >=
qsmall)
then
1392 rhin = (
one + relhum(i,k)) *
half
1393 if (rhin*esloesi >
one)
then
1394 prd = epsi*(rhin*qvl-qvi)*oneoabi
1397 cmei(i,k) = cmei(i,k)
1398 & + prd * min(icldm(i,k),lcldm(i,k))
1409 if (qc(i,k) <
qsmall)
then
1416 if (qc(i,k) <
qsmall)
then
1430 if (rhin*esloesi >
one)
then
1432 prd = epsi*(rhin*qvl-qvi)*oneoabi
1436 cmei(i,k) = cmei(i,k)
1437 & + prd * max((icldm(i,k)-dum), zero)
1447 if(cmei(i,k) > zero .and. relhum(i,k)*esloesi >
one)
1448 & cmei(i,k) = min(cmei(i,k),(q(i,k)-qs(i)/esloesi)
1463 if (-berg(i,k) < -tx1) berg(i,k) = max(tx1, zero)
1469 if (relhum(i,k)*esloesi <
one .and.
1470 & qiic(i,k) >=
qsmall )
then
1472 qvi = epsqs*esi(i,k)/(p(i,k)-omeps*esi(i,k))
1473 qvl = epsqs*esl(i,k)/(p(i,k)-omeps*esl(i,k))
1474 dqsidt =
xxls*qvi/(
rv*t(i,k)*t(i,k))
1483 miu_ice(k) = max(min(0.008_r8*(lami(k)*0.01)**0.87_r8,
1485 tx1 =
one + miu_ice(k)
1488 lami(k) = aux*lami(k)
1490 n0i(k) = niic(i,k) * lami(k)**tx1 * tx2
1493 if (lami(k) <
lammini*aux)
then
1497 else if (lami(k) >
lammaxi*aux)
then
1503 epsi = (miu_ice(k)+
two)*(miu_ice(k)+
one)*pi*
1504 & niic(i,k)*rho(i,k)*dv(i,k) / lami(k)
1507 prd = epsi*(relhum(i,k)*qvl-qvi)*oneoabi * icldm(i,k)
1508 cmei(i,k) = min(prd, zero)
1516 cmei(i,k) = max(cmei(i,k), -qi(i,k)*dti)
1522 if(cmei(i,k) < zero .and. relhum(i,k)*esloesi <
one)
1523 & cmei(i,k) = min(zero, max(cmei(i,k),(q(i,k)-qs(i)
1524 & /esloesi) * oneoabi * dti ))
1531 cmei(i,k) = cmei(i,k)*omsm
1543 dum2i(i,k) = naai(i,k)
1547 dum2i(i,k) = 0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k)))
1551 dum2i(i,k) = min(dum2i(i,k),208.9e3_r8)*irho(i,k)
1571 cldo(i,k) = cldn(i,k)
1584 rate1ord_cw2pr_st(i,k) = zero
1593 if (fprcp == 1)
then
1599 riter =
one / float(iter)
1601 deltat = deltat * riter
1616 if (ltrue(i) == 0)
then
1640 qcsinksum_rate1ord(k) = zero
1641 qcsum_rate1ord(k) = zero
1707 cldmax(i,k) = cldm(i,k)
1711 if (qric(i,km) >=
qsmall .or.
1712 & qniic(i,km) >=
qsmall)
then
1713 cldmax(i,k) = max(cldmax(i,km),cldm(i,k))
1715 cldmax(i,k) = cldm(i,k)
1719 rdz = rho(i,k) * dz(i,k)
1727 if (dum2i(i,k) > zero .and. t(i,k) < (
tmelt-
five) .and.
1728 & relhum(i,k)*esl(i,k)/esi(i,k) >
rhmini+0.05_r8)
then
1733 nimax = dum2i(i,k)*icldm(i,k)
1735 nnuccd(k) = max(zero, (nimax-ni(i,k))*dti)
1740 mnuccd(k) = nnuccd(k) *
mi0
1743 cmei(i,k)= cmei(i,k) + mnuccd(k) * mtime
1749 qvi = epsqs*esi(i,k)/(p(i,k)-omeps*esi(i,k))
1750 dqsidt =
xxls*qvi/(
rv*t(i,k)*t(i,k))
1752 cmei(i,k) = min(cmei(i,k),(q(i,k)-qvi)/(abi*deltat))
1755 cmei(i,k) = cmei(i,k)*omsm
1773 tx1 =
one / lcldm(i,k)
1774 tx2 =
one / icldm(i,k)
1775 qcic(i,k) = min(cwml(i,k)*tx1, 5.e-3_r8)
1776 qiic(i,k) = min(cwmi(i,k)*tx2, 5.e-3_r8)
1777 ncic(i,k) = max(nc(i,k)*tx1, zero)
1778 niic(i,k) = max(ni(i,k)*tx2, zero)
1782 ncic(i,k) = ncnst*irho(i,k)
1786 niic(i,k) = ninst*irho(i,k)
1789 tx1 = qc(i,k) - berg(i,k)*deltat
1793 if (tx1 < zero)
then
1794 berg(i,k) = qc(i,k)*dti*omsm
1798 tx1 = qi(i,k) + (cmei(i,k)+berg(i,k))*deltat
1802 if (tx1 < zero)
then
1803 cmei(i,k) = (-qi(i,k)*dti-berg(i,k))*omsm
1810 cmeout(i,k) = cmeout(i,k) + cmei(i,k)
1817 if (cmei(i,k) < zero .and. qi(i,k) >
qsmall
1818 & .and. cldm(i,k) > mincld)
then
1819 nsubi(k) = cmei(i,k)*ni(i,k) / (qi(i,k)*cldm(i,k))
1837 if (qiic(i,k) >=
qsmall)
then
1840 niic(i,k) = min(niic(i,k),qiic(i,k)*1.e20_r8)
1846 miu_ice(k) = max(min(0.008_r8*(lami(k)*0.01)**0.87_r8,
1848 tx1 =
one + miu_ice(k)
1851 lami(k) = aux*lami(k)
1853 n0i(k) = niic(i,k)*lami(k)**tx1 * tx2
1861 niic(i,k) = n0i(k)/lami(k)
1862 else if (lami(k) >
lammaxi)
then
1865 niic(i,k) = n0i(k)/lami(k)
1873 if (qcic(i,k) >=
qsmall)
then
1876 ncic(i,k) = min(ncic(i,k),qcic(i,k)*1.e20_r8)
1878 ncic(i,k) = max(ncic(i,k),cdnl*irho(i,k))
1882 pgam(k) = 0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))
1887 if ((ncic(i,k) > 1.0e-3) .and.
1888 & (qcic(i,k) > 1.0e-11))
then
1889 xs = 0.07_r8*(1000._r8*qcic(i,k)/ncic(i,k))
1895 xs = max(min(xs, 1.7_r8), 1.1_r8)
1897 xs = (xs + sqrt(xs+8.0_r8)*sqrt(xs) -
four)/8.0_r8
1904 pgam(k) =
one / (pgam(k)*pgam(k)) -
one
1905 pgam(k) = max(
two, min(15._r8, pgam(k)))
1911 tx2 = 6._r8 * qcic(i,k) * tx3
1913 lamc(k) = (tx1*ncic(i,k)/tx2) **
oneb3
1917 lammin = (pgam(k)+
one) / 50.e-6_r8
1918 lammax = (pgam(k)+
one) / 2.e-6_r8
1920 if (lamc(k) < lammin)
then
1922 ncic(i,k) = tx2*lamc(k)*lamc(k)*lamc(k) / tx1
1923 else if (lamc(k) > lammax)
then
1925 ncic(i,k) = tx2*lamc(k)*lamc(k)*lamc(k) / tx1
1930 cdist1(k) = ncic(i,k) / tx3
1947 if (qcic(i,k) >= 1.e-8_r8)
then
1960 if (auto_option == 1)
then
1962 tx1 = qcic(i,k)*rho(i, k)/3.0e-4
1963 prc(k) = 1.0e-3 * qcic(i,k) * (
one-exp(-tx1*tx1))
1966 elseif (auto_option == 2)
then
1968 tx1 = qcic(i,k)*rho(i, k)/7.5e-4
1970 prc(k) = 1.0e-4 * qcic(i,k) * (
one-exp(-tx1*tx1))
1973 elseif (auto_option == 3)
then
1975 tx1 = qcic(i,k) / 3.0e-4
1977 prc(k) = 1.0e-3 * (
one-exp(-tx1*tx1))
1979 & / (1.e-6_r8/50.0*ncic(i,k)*rho(i,k))**1.79_r8
1981 elseif (auto_option == 4)
then
1985 beta6 = (
one+3.0*xs)*(
one+4.0*xs)*(
one+5.0*xs)
1986 & / ((
one+xs)*(
one+xs+xs))
1988 lw = 1.0e-3_r8 * qcic(i,k) * rho(i,k)
1989 nw = ncic(i,k) * rho(i,k) * 1.e-6_r8
1991 xs = min(20.0, 1.03e16*(lw*lw)/(nw*sqrt(nw)))
1992 prc(k) = 1.1e10*beta6*lw*lw*lw
1993 & * (
one-exp(-(xs**miu_disp))) / nw
1994 prc(k) = prc(k)*1.0e3*irho(i,k)
1998 prc(k) = prc(k)*dcrit
2003 & * qcic(i,k)**2.47_r8
2004 & * (1.e-6_r8*ncic(i,k)*rho(i,k))**(-1.79_r8)
2007 nprc(k) = prc(k)/cons22
2009 nprc1(k) = ncic(i,k)*prc(k) / (qcic(i,k)*(
one+xs))
2024 if (fprcp == 1)
then
2025 tx1 =
one / lcldm(i,k)
2026 qric(i,k) = min(qrn(i,k)*tx1, 10.e-3_r8)
2027 nric(i,k) = max(nrn(i,k)*tx1, zero)
2033 tx1 = lcldm(i,k)*dz(i,k)/(cldmax(i,k)*dum)
2034 qric(i,k) = prc(k) * tx1
2035 nric(i,k) = nprc(k) * tx1
2037 if (qric(i,km) >=
qsmall)
then
2041 dum = max(umr(km),dum)
2042 dum1 = max(unr(km),dum1)
2049 if (qric(i,km) >= 1.e-9_r8 .or.
2050 & qniic(i,km) >= 1.e-9_r8)
then
2054 tx1 = rho(i,km) * cldmax(i,km)
2055 tx3 = rho(i,k) * cldmax(i,k)
2056 qric(i,k) = (tx1*umr(km)*qric(i,km)
2057 & + (rdz*((pra(km)+prc(k))*lcldm(i,k)
2058 & + (pre(km)-pracs(km)-mnuccr(km))*cldmax(i,k))))
2068 nric(i,k) = (tx1*unr(km)*nric(i,km)
2069 & + (rdz*(nprc(k)*lcldm(i,k)
2070 & +(-npracs(km)-nnuccr(km)+nragg(km))*cldmax(i,k))))
2080 if (t(i,k) <= 273.15_r8 .and. qiic(i,k) >=
qsmall)
then
2088 nprci(k) = (niic(i, k)/vaux)*exp(-lami(k)*
dcs)
2091 prci(k) = pi*irho(i,k)*niic(i,k)*lami(k)
2093 & * (cons23*tx1+
three*cons24*tx2
2094 & + 6._r8*
dcs*tx1*tx2+6._r8*tx2*tx2)
2095 & * exp(-lami(k)*
dcs)
2100 miu_ice(k) =max(min(0.008_r8*(lami(k)*0.01)**0.87_r8
2101 &, 10.0_r8), 0.1_r8)
2119 if (fprcp == 1)
then
2120 tx1 =
one / icldm(i,k)
2121 qniic(i,k) = min(qsnw(i,k)*tx1, 10.e-3_r8)
2122 nsic(i,k) = max(nsnw(i,k)*tx1, 0._r8)
2125 dum = (asn(i,k)*cons25)
2126 dum1 = (asn(i,k)*cons25)
2129 tx1 = icldm(i,k)*dz(i,k)/(cldmax(i,k)*dum)
2130 qniic(i,k) = prci(k) * tx1
2131 nsic(i,k) = nprci(k) * tx1
2133 if (qniic(i,km) >=
qsmall)
then
2138 tx1 = rho(i,km) * cldmax(i,km)
2139 tx3 = rho(i,k) * cldmax(i,k)
2140 qniic(i,k) = (tx1*ums(km)*qniic(i,km) + (rdz*
2141 & ((prci(k)+prai(km)+psacws(km)+bergs(km))*icldm(i,k)
2142 & +(prds(km)+pracs(km)+mnuccr(km))*cldmax(i,k))))
2150 nsic(i,k) = (tx1*uns(km)*nsic(i,km)
2151 & + (rdz*(nprci(k)*icldm(i,k)+(nsagg(km)
2152 & + nnuccr(km))*cldmax(i,k)))) /(dum1*tx3)
2159 if (qniic(i,k) <
qsmall)
then
2164 if (qric(i,k) <
qsmall)
then
2172 nric(i,k) = max(nric(i,k),zero)
2173 nsic(i,k) = max(nsic(i,k),zero)
2180 if (qric(i,k) >=
qsmall)
then
2182 n0r(k) = nric(i,k)*lamr(k)
2190 n0r(k) = tx1 * tx1 * qric(i,k)/
pirhow
2191 nric(i,k) = n0r(k)/lamr(k)
2192 else if (lamr(k) >
lammaxr)
then
2195 n0r(k) = tx1 * tx1 * qric(i,k)/
pirhow
2196 nric(i,k) = n0r(k)/lamr(k)
2201 tx1 = arn(i,k) / lamr(k) **
br
2202 tx2 = 9.1_r8*rhof(i,k)
2203 unr(k) = min(tx1*
cons4, tx2)
2204 umr(k) = min(tx1*(
cons5/6._r8), tx2)
2216 if (qniic(i,k) >=
qsmall)
then
2217 lams(k) = (cons6*
cs*nsic(i,k) / qniic(i,k))**(
one/
ds)
2218 n0s(k) = nsic(i,k)*lams(k)
2223 if (lams(k) < lammins)
then
2225 n0s(k) = lams(k)**(
ds+
one)*qniic(i,k)/(
cs*cons6)
2226 nsic(i,k) = n0s(k)/lams(k)
2227 else if (lams(k) > lammaxs)
then
2229 n0s(k) = lams(k)**(
ds+
one)*qniic(i,k)/(
cs*cons6)
2230 nsic(i,k) = n0s(k)/lams(k)
2235 tx1 = asn(i,k) / lams(k)**
bs
2236 tx2 = 1.2_r8*rhof(i,k)
2237 ums(k) = min(tx1*(cons8/6._r8), tx2)
2238 uns(k) = min(tx1*cons7, tx2)
2248checked up to here moorthi
2252 if (qcic(i,k) >=
qsmall .and. t(i,k) < 269.15_r8)
then
2256 tx2 =
one / (lamc(k) * lamc(k) * lamc(k))
2257 tx3 = min(
aimm*(273.15_r8-t(i,k)), 25.0)
2258 tx1 = cdist1(k) * tx2 *
bimm * exp(tx3)
2259 mnuccc(k) = cons9/(
cons3*cons19)* pi*
pirhow/36._r8
2260 & *
gamma(7._r8+pgam(k))*tx1*tx2
2266 nnuccc(k) = nimm(i,k)
2267 mnuccc(k) = nimm(i,k)*qcic(i,k)/max(cdist1(k),
one)
2278 tcnt = (270.16_r8-t(i,k))**1.3_r8
2279 viscosity = 1.8e-5_r8*(t(i,k)/298.0_r8)**0.85_r8
2281 & / (p(i,k) *sqrt(8.0_r8*28.96e-3_r8/(pi*
2282 & 8.314409_r8*t(i,k))))
2283 taux = t(i,k) -
three
2284 taux = max(taux-273.15, -40.0)
2286 nsdust = max(exp(-0.517*taux + 8.934) -3.76e6, 0.0)
2287 nssoot = max(1.0e4*exp((-0.0101*taux-0.8525)*taux
2288 & +0.7667) -3.77e9, 0.0)
2290 do ind_aux = 1, nbincontactdust
2291 tx1 = rndst(i,k,ind_aux)
2292 nslip(ind_aux) =
one + (mfp/tx1) *
2293 & (1.257_r8+(0.4_r8*exp(-(1.1_r8*tx1/mfp))))
2294 ndfaer(ind_aux) = 1.381e-23_r8*t(i,k)* nslip(ind_aux)
2295 & / (6._r8*pi*viscosity*tx1)
2297 ndfaer(ind_aux) = ndfaer(ind_aux)
2298 & * (1.0-exp(-nsdust*12.5664*tx1*tx1))
2301 nslipsoot =
one + (mfp/rnsoot(i,k)) *
2302 & (1.257_r8+(0.4_r8*exp(-(1.1_r8*rnsoot(i,k)/mfp))))
2303 ndfaersoot = 1.381e-23_r8*t(i,k)*nslipsoot
2304 & / (6._r8*pi*viscosity*rnsoot(i,k))
2305 ndfaersoot = ndfaersoot *
2306 & (1.0-exp(-nssoot*12.5664*rnsoot(i,k)*rnsoot(i,k)))
2310 tx1 = sum(ndfaer(:)*nacon(i,k,:))+ndfaersoot*nsoot(i,k)
2311 tx1 = tx1 * pi * cdist1(k)
2318 nnucct(k) = (tx1+tx1) *
gamma(pgam(k)+
two) * tx2
2323 if (nnuccc(k)*lcldm(i,k) > nnuccd(k))
then
2324 dum = (nnuccd(k)/(nnuccc(k)*lcldm(i,k)))
2326 mnuccc(k) = mnuccc(k)*dum
2327 nnuccc(k) = nnuccd(k)/lcldm(i,k)
2332 if ((nnucct(k)+nnuccc(k))*deltat > ncic(i, k))
then
2335 nnuccc(k) = ncic(i, k)*dti
2336 mnuccc(k) = cons9/(
cons3*cons19)* pi*
pirhow/36._r8
2337 & * cdist1(k)*
gamma(7._r8+pgam(k))*ncic(i,k)
2357 if (qniic(i,k) >=
qsmall .and. t(i,k) <= 273.15_r8)
then
2360 & * (rho(i,k)*qniic(i,k)/
rhosn) ** tx1
2362 & / (
four*720._r8*rho(i,k))
2377 & .and. qcic(i,k) >=
qsmall)
then
2384 dc0 = (pgam(k)+
one)/lamc(k)
2386 dum = dc0*dc0*uns(k)*
rhow/(9._r8*mu(i,k)*ds0)
2387 eci = dum*dum/((dum+0.4_r8)*(dum+0.4_r8))
2389 eci = min(
one, max(eci, zero))
2395 tx1 = pi/
four*asn(i,k)*rho(i,k)*n0s(k)*eci*cons11
2397 psacws(k) = tx1 * qcic(i,k)
2398 npsacws(k) = tx1 * ncic(i,k)
2407 if((t(i,k) < 270.16_r8) .and. (t(i,k) >= 268.16_r8))
then
2408 ni_secp = 0.5*3.5e8_r8*(270.16_r8-t(i,k))*psacws(k)
2410 msacwi(k) = min(ni_secp*
mi0, psacws(k))
2411 else if((t(i,k) < 268.16_r8) .and.
2412 & (t(i,k) >= 265.16_r8))
then
2413 ni_secp =
oneb3*3.5e8_r8*(t(i,k)-265.16_r8)*psacws(k)
2415 msacwi(k) = min(ni_secp*
mi0, psacws(k))
2421 psacws(k) = max(zero, psacws(k)-ni_secp*
mi0)
2427 if (qric(i,k) >= 1.e-8_r8 .and. qniic(i,k) >= 1.e-8_r8
2428 & .and. t(i,k) <= 273.15_r8)
then
2430 tx1 = 1.2_r8*umr(k) - 0.95_r8*ums(k)
2431 tx2 = sqrt(tx1*tx1+0.08_r8*ums(k)*umr(k))
2435 tx5 = pi *
ecr * rho(i,k) *n0r(k) * n0s(k)
2437 pracs(k) =
pirhow*tx2*tx5 *
2441 tx2 = unr(k) - uns(k)
2442 tx2 = sqrt(1.7_r8*tx2*tx2 + 0.3_r8*unr(k)*uns(k))
2444 npracs(k) =
half*tx2*tx5*tx1*tx3*(tx4+tx3*(tx1+tx3))
2458 if (t(i,k) < 269.15_r8 .and. qric(i,k) >=
qsmall
2459 & .and. t(i,k) > 260.0_r8)
then
2462 tx1 = exp(min(
aimm*(273.15_r8-t(i,k))-
one, 25.0))
2463 tx2 = 1.0 / (lamr(k)*lamr(k)*lamr(k))
2465 nnuccr(k) = pi * nric(i,k) *
bimm * tx1 * tx2
2466 mnuccr(k) = 20._r8 *
pirhow * nnuccr(k) * tx2
2479 if (qric(i,k) >=
qsmall .and. qcic(i,k) >=
qsmall)
then
2483 pra(k) = cons12/(
cons3*cons20)
2484 & * 67._r8*(qcic(i,k)*qric(i,k))**1.15_r8
2485 npra(k) = pra(k) * (ncic(i,k)/qcic(i,k))
2497 if (qric(i,k) >=
qsmall)
then
2498 nragg(k) = -8._r8*nric(i,k)*qric(i,k)*rho(i,k)
2509 & .and.t(i,k) <= 273.15_r8)
then
2511 tx1 = (pi/
four)*asn(i,k)*rho(i,k)*n0s(k)*
eii*cons11
2513 prai(k) = tx1 * qiic(i,k)
2514 nprai(k) = tx1 * niic(i,k)
2516 nprai(k)= min(nprai(k), 1.0e10)
2536 if (qcic(i,k)+qiic(i,k) < 1.e-7_r8 .or.
2537 & cldmax(i,k) > lcldm(i,k))
then
2542 if (qcic(i,k)+qiic(i,k) < 1.e-7_r8)
then
2550 esn = min(
fpvsl(t(i,k)), p(i,k))
2551 qsn = min(epsqs*esn/(p(i,k)-omeps*esn),
one)
2556 esi(i,k) = min(
fpvsi(t(i,k)), p(i,k))
2558 if (t(i,k) >
tmelt) esi(i,k) = esl(i,k)
2562 qclr = (q(i,k)-dum*qsn) / (
one-dum)
2564 if (qric(i,k) >=
qsmall)
then
2566 qvs = epsqs*esl(i,k)/(p(i,k)-omeps*esl(i,k))
2567 dqsdt =
xxlv*qvs/(
rv*t(i,k)*t(i,k))
2569 epsr = (pi+pi)*n0r(k)*rho(i,k)*dv(i,k)
2570 & * (
f1r/(lamr(k)*lamr(k))
2571 & +
f2r*sqrt(arn(i,k)*rho(i,k)/mu(i,k))
2575 pre(k) = epsr*(qclr-qvs)/ab
2579 pre(k) = min(pre(k)*(cldmax(i,k)-dum), zero)
2580 pre(k) = pre(k) / cldmax(i,k)
2585 if (qniic(i,k) >=
qsmall)
then
2586 qvi = epsqs*esi(i,k)/(p(i,k)-omeps*esi(i,k))
2587 dqsidt =
xxls*qvi/(
rv*t(i,k)*t(i,k))
2589 epss = (pi+pi)*n0s(k)*rho(i,k)*dv(i,k)
2590 & * (
f1s/(lams(k)*lams(k))
2591 & +
f2s*sqrt((asn(i,k)*rho(i,k)/mu(i,k)))
2594 prds(k) = epss*(qclr-qvi)/abi
2597 prds(k) = min(prds(k)*(cldmax(i,k)-dum), zero)
2598 prds(k) = prds(k)/cldmax(i,k)
2605 tx1 = pre(k) * cldmax(i,k)
2606 tx2 = prds(k) * cldmax(i,k)
2607 qtmp = q(i,k) - (cmei(i,k)+tx1+tx2) * deltat
2608 ttmp = t(i,k) + (tx1*
xxlv + (cmei(i,k)+tx2)*
xxls)
2612 ttmp = max(180._r8,min(ttmp,323._r8))
2615 esn = min(
fpvsl(ttmp), p(i,k))
2616 qsn = min(epsqs*esn/(p(i,k)-omeps*esn),
one)
2619 if (qtmp > qsn)
then
2620 if (pre(k)+prds(k) < -1.e-20)
then
2621 dum1 = pre(k) / (pre(k)+prds(k))
2623 qtmp = q(i,k) - cmei(i,k)*deltat
2624 ttmp = t(i,k) + cmei(i,k)*
xxls*deltat/
cpp
2626 esn = min(
fpvsl(ttmp), p(i,k))
2627 qsn = min(epsqs*esn/(p(i,k)-omeps*esn),
one)
2629 dum = min(zero, (qtmp-qsn)/(
one+cons27*qsn*tx1))
2632 pre(k) = dum*dum1/(deltat*cldmax(i,k))
2636 esn = min(
fpvsi(ttmp), p(i,k))
2637 qsn = min(epsqs*esn/(p(i,k)-omeps*esn),
one)
2638 dum = min(zero, (qtmp-qsn)/(
one+cons28*qsn*tx1))
2641 prds(k) = dum*(
one-dum1)/(deltat*cldmax(i,k))
2649 & .and. t(i,k) <
tmelt)
then
2650 qvi = epsqs*esi(i,k)/(p(i,k)-omeps*esi(i,k))
2651 qvs = epsqs*esl(i,k)/(p(i,k)-omeps*esl(i,k))
2652 dqsidt =
xxls*qvi/(
rv*t(i,k)*t(i,k))
2654 epss = (pi+pi)*n0s(k)*rho(i,k)*dv(i,k)
2655 & * (
f1s/(lams(k)*lams(k))
2656 & +
f2s*sqrt(asn(i,k)*rho(i,k)/mu(i,k))
2659 bergs(k) = epss*(qvs-qvi)/abi
2678 qce = max(qc(i,k)-berg(i,k)*deltat, zero)
2679 nce = nc(i,k) + npccn(i,k)*deltat*mtime
2680 qie = qi(i,k) + (cmei(i,k)+berg(i,k))*deltat
2681 nie = ni(i,k) + nnuccd(k)*deltat*mtime
2685 tx1 = lcldm(i,k) * deltat
2686 dum = (prc(k) + pra(k) + mnuccc(k) + mnucct(k)
2687 & + msacwi(k) + psacws(k) + bergs(k)) * tx1
2691 ratio = qce/dum * omsm
2693 prc(k) = prc(k) * ratio
2694 pra(k) = pra(k) * ratio
2695 mnuccc(k) = mnuccc(k) * ratio
2696 mnucct(k) = mnucct(k) * ratio
2697 msacwi(k) = msacwi(k) * ratio
2698 psacws(k) = psacws(k) * ratio
2699 bergs(k) = bergs(k) * ratio
2705 dum = (nprc1(k) + npra(k) + nnuccc(k) + nnucct(k)
2706 & + npsacws(k) - nsubc(k)) * tx1
2709 ratio = nce/dum * omsm
2711 nprc1(k) = nprc1(k) * ratio
2712 npra(k) = npra(k) * ratio
2713 nnuccc(k) = nnuccc(k) * ratio
2714 nnucct(k) = nnucct(k) * ratio
2715 npsacws(k) = npsacws(k) * ratio
2716 nsubc(k) = nsubc(k) * ratio
2721 dum = ((-mnuccc(k)-mnucct(k)-msacwi(k))*lcldm(i,k)
2722 & + (prci(k)+prai(k))*icldm(i,k)) * deltat
2725 if (prci(k)+prai(k) > zero)
then
2727 & +(mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k))
2728 & / ((prci(k)+prai(k))*icldm(i,k))*omsm
2733 prci(k) = prci(k)*ratio
2734 prai(k) = prai(k)*ratio
2739 dum = ((-nnucct(k)-nsacwi(k))*lcldm(i,k)
2740 & + (nprci(k)+ nprai(k)-nsubi(k))*icldm(i,k))*deltat
2743 if ( abs(nprci(k)+nprai(k)-nsubi(k)) > zero)
then
2744 ratio = (nie*dti+(nnucct(k)+nsacwi(k))*lcldm(i,k))
2745 & / ((nprci(k)+nprai(k)-nsubi(k))*icldm(i,k))*omsm
2750 nprci(k) = nprci(k) * ratio
2751 nprai(k) = nprai(k) * ratio
2752 nsubi(k) = nsubi(k) * ratio
2760 if ( ((prc(k)+pra(k))*lcldm(i,k)
2761 & + (-mnuccr(k)+pre(k)-pracs(k))*cldmax(i,k))*rdz
2762 & + qrtot < zero)
then
2764 if (-pre(k)+pracs(k)+mnuccr(k) >=
qsmall)
then
2766 ratio = (qrtot*rdzi + (prc(k)+pra(k))*lcldm(i,k))
2767 & / ((-pre(k)+pracs(k)+mnuccr(k))*cldmax(i,k))*omsm
2769 pre(k) = pre(k) * ratio
2770 pracs(k) = pracs(k) * ratio
2771 mnuccr(k) = mnuccr(k) * ratio
2779 if ((nprc(k)*lcldm(i,k)+(-nnuccr(k)+nsubr(k)-npracs(k) +
2780 & nragg(k))*cldmax(i,k))*rdz + nrtot < zero)
then
2782 if (-nsubr(k)-nragg(k)+npracs(k)+nnuccr(k) >=
qsmall)
2784 ratio = (nrtot*rdzi+nprc(k)*lcldm(i,k))
2785 & / ((-nsubr(k)-nragg(k)+npracs(k)
2786 & +nnuccr(k))*cldmax(i,k))*omsm
2788 nsubr(k) = nsubr(k) * ratio
2789 npracs(k) = npracs(k) * ratio
2790 nnuccr(k) = nnuccr(k) * ratio
2791 nragg(k) = nragg(k) * ratio
2797 tx1 = (bergs(k)+psacws(k))*lcldm(i,k)
2798 & + (prai(k)+prci(k))*icldm(i,k)
2799 if ((tx1+(pracs(k)+ mnuccr(k)+prds(k))*cldmax(i,k))
2800 & *rdz+qstot < zero)
then
2802 if (-prds(k) >=
qsmall)
then
2804 ratio = (qstot*rdzi + tx1
2805 & + (pracs(k)+mnuccr(k))*cldmax(i,k))
2806 & / (-prds(k)*cldmax(i,k))*omsm
2808 prds(k) = prds(k) * ratio
2818 if ((nprci(k)*icldm(i,k)
2819 & +(nnuccr(k)+nsubs(k)+nsagg(k))*cldmax(i,k))
2820 & * rdz + nstot < zero)
then
2822 if (-nsubs(k)-nsagg(k) >=
qsmall)
then
2825 & + nprci(k)*icldm(i,k)+ nnuccr(k)*cldmax(i,k))
2826 & / ((-nsubs(k)-nsagg(k))*cldmax(i,k))*omsm
2828 nsubs(k) = nsubs(k) * ratio
2829 nsagg(k) = nsagg(k) * ratio
2838 qvlat(i,k) = qvlat(i,k)
2839 & - (pre(k)+prds(k))*cldmax(i,k) - cmei(i,k)
2844 tlat(i,k) = tlat(i,k)
2845 & + pre(k)*cldmax(i,k) *
xxlv
2846 & + (prds(k)*cldmax(i,k)+cmei(i,k)) *
xxls
2847 & + ((bergs(k)+psacws(k)+mnuccc(k)+
2848 & mnucct(k)+msacwi(k))*lcldm(i,k)
2849 & + (mnuccr(k)+ pracs(k))*cldmax(i,k)
2850 & + berg(i,k)) *
xlf
2859 qctend(i,k) = qctend(i,k)
2860 & + (-pra(k)-prc(k)-mnuccc(k)-mnucct(k)
2861 & -msacwi(k)-psacws(k)-bergs(k))*lcldm(i,k)
2864 qitend(i,k) = qitend(i,k)
2865 & + (mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k)
2866 & + (-prci(k)-prai(k))*icldm(i,k)
2867 & + cmei(i,k) + berg(i,k)
2869 qrtend(i,k) = qrtend(i,k)
2870 & + (pra(k)+prc(k))*lcldm(i,k)
2871 & + (pre(k)-pracs(k)-mnuccr(k))*cldmax(i,k)
2873 qnitend(i,k) = qnitend(i,k)
2874 & + (prai(k)+prci(k))*icldm(i,k)
2875 & + (psacws(k)+bergs(k))*lcldm(i,k)
2876 & + (prds(k)+pracs(k)+mnuccr(k))*cldmax(i,k)
2883 cmeiout(i,k) = cmeiout(i,k) + cmei(i,k)
2892 evapsnow(i,k) = evapsnow(i,k) + prds(k) * cldmax(i,k)
2893 nevapr(i,k) = nevapr(i,k) + pre(k) * cldmax(i,k)
2897 prain(i,k) = prain(i,k)
2898 & + (pra(k)+prc(k))*lcldm(i,k)
2899 & + (-pracs(k)-mnuccr(k))*cldmax(i,k)
2901 prodsnow(i,k) = prodsnow(i,k)
2902 & + (prai(k)+prci(k))*icldm(i,k)
2903 & + (psacws(k)+bergs(k))*lcldm(i,k)
2904 & + (pracs(k)+mnuccr(k))*cldmax(i,k)
2914 qcsinksum_rate1ord(k) = qcsinksum_rate1ord(k)
2915 & +(pra(k)+prc(k)+psacws(k))*lcldm(i,k)
2916 qcsum_rate1ord(k) = qcsum_rate1ord(k) + qc(i,k)
2919 prao(i,k) = prao(i,k) + pra(k) * lcldm(i,k)
2920 prco(i,k) = prco(i,k) + prc(k) * lcldm(i,k)
2921 mnuccco(i,k) = mnuccco(i,k) + mnuccc(k) * lcldm(i,k)
2922 mnuccto(i,k) = mnuccto(i,k) + mnucct(k) * lcldm(i,k)
2923 msacwio(i,k) = msacwio(i,k) + msacwi(k) * lcldm(i,k)
2924 psacwso(i,k) = psacwso(i,k) + psacws(k) * lcldm(i,k)
2925 bergso(i,k) = bergso(i,k) + bergs(k) * lcldm(i,k)
2926 bergo(i,k) = bergo(i,k) + berg(i,k)
2927 prcio(i,k) = prcio(i,k) + prci(k) * icldm(i,k)
2928 praio(i,k) = praio(i,k) + prai(k) * icldm(i,k)
2929 mnuccro(i,k) = mnuccro(i,k) + mnuccr(k) * cldmax(i,k)
2930 pracso(i,k) = pracso(i,k) + pracs(k) * cldmax(i,k)
2932 mnuccdo(i,k) = mnuccdo(i,k) + mnuccd(k)
2933 nnuccto(i,k) = nnuccto(i,k) + nnucct(k) * lcldm(i,k)
2935 nnuccdo(i,k) = nnuccdo(i,k) + nnuccd(k)
2936 nnuccco(i,k) = nnuccco(i,k) + nnuccc(k) * lcldm(i,k)
2937 nsacwio(i,k) = nsacwio(i,k) + nsacwi(k) * icldm(i,k)
2938 nsubio(i,k) = nsubio(i,k) + nsubi(k) * icldm(i,k)
2939 nprcio(i,k) = nprcio(i,k) + nprci(k) * icldm(i,k)
2940 npraio(i,k) = npraio(i,k) + nprai(k) * icldm(i,k)
2942 npccno(i, k) = npccno(i,k) + npccn(i,k)*lcldm(i,k)
2943 npsacwso(i,k) = npsacwso(i,k) + npsacws(k)*lcldm(i,k)
2944 nsubco(i,k) = nsubco(i,k) + nsubc(k) * lcldm(i,k)
2945 nprao(i,k) = nprao(i,k) + npra(k) * lcldm(i,k)
2946 nprc1o(i,k) = nprc1o(i,k) + nprc1(k) * lcldm(i,k)
2951 nctend(i,k) = nctend(i,k) + npccn(i,k)*mtime
2952 & + (-nnuccc(k)-nnucct(k)-npsacws(k)+nsubc(k)
2953 & -npra(k)-nprc1(k))*lcldm(i,k)
2955 nitend(i,k) = nitend(i,k) + nnuccd(k)*mtime
2956 & + (nnuccc(k)+nnucct(k)+nsacwi(k))*lcldm(i,k)
2957 & + (nsubi(k)-nprci(k)- nprai(k))*icldm(i,k)
2959 nstend(i,k) = nstend(i,k)
2960 & + (nsubs(k)+ nsagg(k)+nnuccr(k))*cldmax(i,k)
2961 & + nprci(k)*icldm(i,k)
2963 nrtend(i,k) = nrtend(i,k) + nprc(k)*lcldm(i,k)
2964 & + (nsubr(k)-npracs(k)-nnuccr(k)+nragg(k))
2971 if (nctend(i,k) > zero .and.
2972 & nc(i,k)+nctend(i,k)*deltat > ncmax(i,k))
then
2973 nctend(i,k) = max(zero, (ncmax(i,k)-nc(i,k))*dti)
2975 if (nitend(i,k) > zero .and.
2976 & ni(i,k)+nitend(i,k)*deltat > nimax)
then
2977 nitend(i,k) = max(zero, (nimax-ni(i,k))*dti)
2984 if (fprcp == 0)
then
2987 if (qric(i,k) >=
qsmall)
then
2989 qric(i,k) = qrtend(i,k)*dz(i,k)/(cldmax(i,k)*umr(k))
2990 nric(i,k) = nrtend(i,k)*dz(i,k)/(cldmax(i,k)*unr(k))
2992 tx1 = rho(i,km) * cldmax(i,km)
2993 tx3 = rho(i,k) * cldmax(i,k)
2995 qric(i,k) = (tx1*umr(km)*qric(i,km)
2996 & + rdz*qrtend(i,k)) / (umr(k)*tx3)
2998 nric(i,k) = (tx1*unr(km)*nric(i,km)
2999 & + rdz*nrtend(i,k)) / (unr(k)*tx3)
3009 if (qniic(i,k) >=
qsmall)
then
3011 tx1 = dz(i,k)/cldmax(i,k)
3012 qniic(i,k) = qnitend(i,k)*tx1/ums(k)
3013 nsic(i,k) = nstend(i,k)*tx1/uns(k)
3015 tx1 = rho(i,km) * cldmax(i,km)
3016 tx3 = rho(i,k) * cldmax(i,k)
3018 qniic(i,k) = (tx1*ums(km)*qniic(i,km)
3019 & + rdz*qnitend(i,k)) / (ums(k)*tx3)
3020 nsic(i,k) = (tx1*uns(km)*nsic(i,km)
3021 & + rdz*nstend(i,k)) / (uns(k)*tx3)
3032 prect(i) = prect(i) + (qrtend(i,k)+ qnitend(i,k))*tx1
3033 preci(i) = preci(i) + qnitend(i,k)*tx1
3041 rainrt(i,k) = qric(i,k)*rho(i,k)*umr(k)
3042 & /
rhow*3600._r8*1000._r8
3046 qrtot = max(qrtot + qrtend(i,k)*rdz, zero)
3047 qstot = max(qstot + qnitend(i,k)*rdz, zero)
3048 nrtot = max(nrtot + nrtend(i,k)*rdz, zero)
3049 nstot = max(nstot + nstend(i,k)*rdz, zero)
3058 tx1 = t(i,k) + tlat(i,k)*
onebcp*deltat
3059 if (tx1 > 275.15_r8)
then
3061 if (qstot > zero)
then
3065 if (tx1-dum*deltat < 273.15_r8)
then
3066 tx2 = (tx1-275.15_r8) * dti
3067 dum = min(
one,max(zero,tx2/dum))
3072 qric(i,k) = qric(i,k) + dum*qniic(i,k)
3073 nric(i,k) = nric(i,k) + dum*nsic(i,k)
3074 qniic(i,k) = (
one-dum)*qniic(i,k)
3075 nsic(i,k) = (
one-dum)*nsic(i,k)
3077 tmp = -
xlf*dum*qstot*rdzi
3078 meltsdt(i,k) = meltsdt(i,k) + tmp
3080 tlat(i,k) = tlat(i,k) + tmp
3086 qrtot = qrtot + dum*qstot
3087 nrtot = nrtot + dum*nstot
3088 qstot = (
one-dum)*qstot
3089 nstot = (
one-dum)*nstot
3090 preci(i) = (
one-dum)*preci(i)
3095 tlataux(i, k) = tlat(i, k)
3100 tx1 = t(i,k) + tlat(i,k)*
onebcp*deltat
3103 if (qrtot > zero)
then
3109 dum = min(
one,max(zero,tx2/dum))
3113 qniic(i,k) = qniic(i,k) + dum*qric(i,k)
3114 nsic(i,k) = nsic(i,k) + dum*nric(i,k)
3115 qric(i,k) = (
one-dum)*qric(i,k)
3116 nric(i,k) = (
one-dum)*nric(i,k)
3118 tmp =
xlf*dum*qrtot*rdzi
3119 frzrdt(i,k) = frzrdt(i,k) + tmp
3121 tlat(i,k) = tlat(i,k) + tmp
3127 qstot = qstot + dum*qrtot
3128 qrtot = (
one-dum)*qrtot
3129 nstot = nstot + dum*nrtot
3130 nrtot = (
one-dum)*nrtot
3131 preci(i) = preci(i) + dum*(prect(i)-preci(i))
3144 if (qniic(i,k) <
qsmall)
then
3149 if (qric(i,k) <
qsmall)
then
3157 nric(i,k) = max(nric(i,k), zero)
3158 nsic(i,k) = max(nsic(i,k), zero)
3165 if (qric(i,k) >=
qsmall)
then
3167 n0r(k) = nric(i,k)*lamr(k)
3175 tx1 = lamr(k) * lamr(k)
3176 n0r(k) = tx1*tx1*qric(i,k)/
pirhow
3177 nric(i,k) = n0r(k)/lamr(k)
3178 else if (lamr(k) >
lammaxr)
then
3180 tx1 = lamr(k) * lamr(k)
3181 n0r(k) = tx1*tx1*qric(i,k)/
pirhow
3182 nric(i,k) = n0r(k)/lamr(k)
3188 tx1 = arn(i,k) / lamr(k)**
br
3189 tx2 = 9.1_r8*rhof(i,k)
3190 unr(k) = min(tx1*
cons4, tx2)
3191 umr(k) = min(tx1*(
cons5/6._r8),tx2)
3202 if (lamr(k) > zero)
then
3203 artmp = n0r(k) * (0.5*pi) / (lamr(k)*lamr(k)*lamr(k))
3208 if (lamc(k) > zero)
then
3209 actmp = cdist1(k) * pi *
gamma(pgam(k)+
three)
3210 & / (
four * lamc(k)*lamc(k))
3215 if (actmp > zero .or.artmp > zero)
then
3216 rercld(i,k) = rercld(i,k)+
three*(qric(i,k)+qcic(i,k))
3218 arcld(i,k) = arcld(i,k) +
one
3224 if (qniic(i,k) >=
qsmall)
then
3225 lams(k) = (cons6*
cs*nsic(i,k) / qniic(i,k))**(
one/
ds)
3226 n0s(k) = nsic(i,k)*lams(k)
3231 if (lams(k) < lammins)
then
3233 n0s(k) = lams(k)**(
ds+
one)*qniic(i,k)/(
cs*cons6)
3234 nsic(i,k) = n0s(k)/lams(k)
3236 else if (lams(k) > lammaxs)
then
3238 n0s(k) = lams(k)**(
ds+
one)*qniic(i,k)/(
cs*cons6)
3239 nsic(i,k) = n0s(k)/lams(k)
3244 tx1 = asn(i,k) / lams(k)**
bs
3245 tx2 = 1.2_r8*rhof(i,k)
3246 ums(k) = min(tx1*(cons8/6._r8), tx2)
3247 uns(k) = min(tx1*cons7, tx2)
3262 qrout(i,k) = qrout(i,k) + qric(i,k)*cldmax(i,k)
3263 qsout(i,k) = qsout(i,k) + qniic(i,k)*cldmax(i,k)
3264 tx1 = rho(i,k)*cldmax(i,k)
3265 nrout(i,k) = nrout(i,k) + nric(i,k)*tx1
3266 nsout(i,k) = nsout(i,k) + nsic(i,k)*tx1
3269 tlat1(i,k) = tlat1(i,k) + tlat(i,k)
3270 tlat1_aux(i,k) = tlat1_aux(i,k) + tlataux(i,k)
3272 qvlat1(i,k) = qvlat1(i,k) + qvlat(i,k)
3273 qctend1(i,k) = qctend1(i,k) + qctend(i,k)
3274 qitend1(i,k) = qitend1(i,k) + qitend(i,k)
3275 nctend1(i,k) = nctend1(i,k) + nctend(i,k)
3276 nitend1(i,k) = nitend1(i,k) + nitend(i,k)
3278 t(i,k) = t(i,k) + tlat(i,k) * deltat/
cpp
3279 q(i,k) = q(i,k) + qvlat(i,k) * deltat
3280 qc(i,k) = qc(i,k) + qctend(i,k) * deltat
3281 qi(i,k) = qi(i,k) + qitend(i,k) * deltat
3282 nc(i,k) = nc(i,k) + nctend(i,k) * deltat
3283 ni(i,k) = ni(i,k) + nitend(i,k) * deltat
3285 rainrt1(i,k) = rainrt1(i,k) + rainrt(i,k)
3288 if (arcld(i,k) > zero)
then
3289 rercld(i,k) = rercld(i,k) / arcld(i,k)
3297 rflx(i,k+1) = qrout(i,k)*rho(i,k)*umr(k)
3298 sflx(i,k+1) = qsout(i,k)*rho(i,k)*ums(k)
3301 rflx1(i,k+1) = rflx1(i,k+1) + rflx(i,k+1)
3302 sflx1(i,k+1) = sflx1(i,k+1) +
sflx(i,k+1)
3308 prect1(i) = prect1(i) + prect(i)
3309 preci1(i) = preci1(i) + preci(i)
3317 rate1ord_cw2pr_st(i,k) = qcsinksum_rate1ord(k)
3319 & / max(qcsum_rate1ord(k),1.0e-30_r8)
3327 deltat = deltat*real(iter)
3335 if (ltrue(i) == 0)
then
3341 effc_fn(i,k) = 10._r8
3353 prect(i) = prect1(i) * riter
3354 preci(i) = preci1(i) * riter
3372 tlat(i,k) = tlat1(i,k) * riter
3373 qvlat(i,k) = qvlat1(i,k) * riter
3374 qctend(i,k) = qctend1(i,k) * riter
3375 qitend(i,k) = qitend1(i,k) * riter
3376 nctend(i,k) = nctend1(i,k) * riter
3377 nitend(i,k) = nitend1(i,k) * riter
3379 tlataux(i,k) = tlat1_aux(i,k) * riter
3382 rainrt(i,k) = rainrt1(i,k) * riter
3385 rflx(i,k+1) = rflx1(i,k+1) * riter
3386 sflx(i,k+1) = sflx1(i,k+1) * riter
3390 qrout(i,k) = qrout(i,k) * riter
3391 qsout(i,k) = qsout(i,k) * riter
3392 nrout(i,k) = nrout(i,k) * riter
3393 nsout(i,k) = nsout(i,k) * riter
3397 nevapr(i,k) = nevapr(i,k) * riter
3398 evapsnow(i,k) = evapsnow(i,k) * riter
3399 prain(i,k) = prain(i,k) * riter
3400 prodsnow(i,k) = prodsnow(i,k) * riter
3401 cmeout(i,k) = cmeout(i,k) * riter
3403 cmeiout(i,k) = cmeiout(i,k) * riter
3404 meltsdt(i,k) = meltsdt(i,k) * riter
3405 frzrdt(i,k) = frzrdt(i,k) * riter
3409 prao(i,k) = prao(i,k) * riter
3410 prco(i,k) = prco(i,k) * riter
3411 mnuccco(i,k) = mnuccco(i,k) * riter
3412 mnuccto(i,k) = mnuccto(i,k) * riter
3413 msacwio(i,k) = msacwio(i,k) * riter
3414 psacwso(i,k) = psacwso(i,k) * riter
3415 bergso(i,k) = bergso(i,k) * riter
3416 bergo(i,k) = bergo(i,k) * riter
3417 prcio(i,k) = prcio(i,k) * riter
3418 praio(i,k) = praio(i,k) * riter
3420 mnuccdo(i,k) = mnuccdo(i,k) * riter
3421 mnuccto(i,k) = mnuccto(i,k) * riter
3423 mnuccro(i,k) = mnuccro(i,k) * riter
3424 pracso(i,k) = pracso(i,k) * riter
3427 nnuccdo(i,k) = nnuccdo(i,k) * riter
3428 nnuccco(i,k) = nnuccco(i,k) * riter
3429 nsacwio(i,k) = nsacwio(i,k) * riter
3430 nsubio(i,k) = nsubio(i,k) * riter
3431 nprcio(i,k) = nprcio(i,k) * riter
3432 npraio(i,k) = npraio(i,k) * riter
3434 npccno(i,k) = npccno(i,k) * riter
3435 npsacwso(i,k) = npsacwso(i,k) * riter
3436 nsubco(i,k) = nsubco(i,k) * riter
3437 nprao(i,k) = nprao(i,k) * riter
3438 nprc1o(i,k) = nprc1o(i,k) * riter
3445 prain(i,k) = prain(i,k) + prodsnow(i,k)
3455 tx1 =
one / lcldm(i,k)
3456 tx2 =
one / icldm(i,k)
3457 dumc(i,k) = (qc(i,k) + qctend(i,k)*deltat) * tx1
3458 dumi(i,k) = (qi(i,k) + qitend(i,k)*deltat) * tx2
3459 dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat)*tx1, zero)
3460 dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat)*tx2, zero)
3463 dumnc(i,k) = ncnst*irho(i,k)
3467 dumni(i,k) = ninst*irho(i,k)
3472 if (dumi(i,k) >=
qsmall)
then
3474 dumni(i,k) = min(dumni(i,k),dumi(i,k)*1.e20_r8)
3478 miu_ice(k) = max(min(0.008_r8*(lami(k)*0.01)**0.87_r8,
3480 tx1 =
one + miu_ice(k)
3483 lami(k) = aux*lami(k)
3485 n0i(k) = niic(i,k)*lami(k)**tx1 * tx2
3488 if (lami(k) <
lammini*aux)
then
3491 niic(i,k) = n0i(k)/lami(k)
3493 if (lami(k) >
lammaxi*aux)
then
3496 niic(i,k) = n0i(k)/lami(k)
3503 if (dumc(i,k) >=
qsmall)
then
3505 dumnc(i,k) = min(dumnc(i,k),dumc(i,k)*1.e20_r8)
3507 dumnc(i,k) = max(dumnc(i,k),cdnl*irho(i,k))
3508 tx1 = 0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
3509 pgam(k) = max(
two, min(15._r8,
one/(tx1*tx1)-
one))
3513 lammin = (pgam(k)+
one) / 50.e-6_r8
3514 lammax = (pgam(k)+
one) / 2.e-6_r8
3515 lamc(k) = min(lammax, max(lamc(k),lammin))
3524 if (dumc(i,k) >=
qsmall)
then
3541 if (dumi(i,k) >=
qsmall)
then
3545 tx1 = ain(i,k) / lami(k)**
bi
3546 tx2 = 1.2_r8*rhof(i,k)
3547 uni = min(tx1*cons16, tx2)
3548 umi = min(tx1*cons17, tx2)
3570 rgvm = max(fi(k),fc(k),fni(k),fnc(k))
3571 nstep = max(int(rgvm*deltat/pdel(i,k)+
one),nstep)
3576 dumc(i,k) = qc(i,k) + qctend(i,k)*deltat
3577 dumi(i,k) = qi(i,k) + qitend(i,k)*deltat
3578 dumnc(i,k) = max(nc(i,k) + nctend(i,k)*deltat, zero)
3579 dumni(i,k) = max(ni(i,k) + nitend(i,k)*deltat, zero)
3581 if (dumc(i,k) <
qsmall) dumnc(i,k) = zero
3582 if (dumi(i,k) <
qsmall) dumni(i,k) = zero
3589 falouti(k) = max(fi(k) * dumi(i,k), zero)
3590 faloutni(k) = max(fni(k) * dumni(i,k), zero)
3591 faloutc(k) = max(fc(k) * dumc(i,k), zero)
3592 faloutnc(k) = max(fnc(k) * dumnc(i,k), zero)
3598 tx1 =
one / pdel(i,k)
3599 faltndi = falouti(k) * tx1
3600 faltndni = faloutni(k) * tx1
3601 faltndc = faloutc(k) * tx1
3602 faltndnc = faloutnc(k) * tx1
3606 tx2 =
one / float(nstep)
3609 qitend(i,k) = qitend(i,k) - faltndi * tx2
3610 nitend(i,k) = nitend(i,k) - faltndni * tx2
3611 qctend(i,k) = qctend(i,k) - faltndc * tx2
3612 nctend(i,k) = nctend(i,k) - faltndnc * tx2
3615 qcsedten(i,k) = qcsedten(i,k) - faltndc * tx2
3616 qisedten(i,k) = qisedten(i,k) - faltndi * tx2
3618 dumi(i,k) = dumi(i,k) - faltndi *tx3
3619 dumni(i,k) = dumni(i,k) - faltndni*tx3
3620 dumc(i,k) = dumc(i,k) - faltndc *tx3
3621 dumnc(i,k) = dumnc(i,k) - faltndnc*tx3
3630 tx1 =
one / pdel(i,k)
3632 if (lcldm(i,k-1) > zero)
then
3633 dum = min(
one, lcldm(i,k)/lcldm(i,k-1))
3635 dum = min(
one, lcldm(i,k))
3637 if (icldm(i,k-1) > zero)
then
3638 dum1 = min(
one, icldm(i,k)/icldm(i,k-1))
3640 dum1 = min(
one, icldm(i,k))
3643 faltndqie = (falouti(k) - falouti(k-1)) * tx1
3644 faltndi = (falouti(k) - dum1*falouti(k-1)) * tx1
3645 faltndni = (faloutni(k) - dum1*faloutni(k-1)) * tx1
3647 faltndqce = (faloutc(k) - faloutc(k-1)) * tx1
3648 faltndc = (faloutc(k) - dum*faloutc(k-1)) * tx1
3649 faltndnc = (faloutnc(k)- dum*faloutnc(k-1)) * tx1
3653 qitend(i,k) = qitend(i,k) - faltndi * tx2
3654 nitend(i,k) = nitend(i,k) - faltndni * tx2
3655 qctend(i,k) = qctend(i,k) - faltndc * tx2
3656 nctend(i,k) = nctend(i,k) - faltndnc * tx2
3660 qcsedten(i,k) = qcsedten(i,k) - faltndc * tx2
3661 qisedten(i,k) = qisedten(i,k) - faltndi * tx2
3665 qvlat(i,k) = qvlat(i,k) - (faltndqie-faltndi) * tx2
3670 qisevap(i,k) = qisevap(i,k) + (faltndqie-faltndi) * tx2
3673 qvlat(i,k) = qvlat(i,k) - (faltndqce-faltndc) * tx2
3677 qcsevap(i,k) = qcsevap(i,k) + (faltndqce-faltndc) * tx2
3680 tlat(i,k) = tlat(i,k) + ((faltndqie-faltndi)*
xxls
3681 & + (faltndqce-faltndc)*
xxlv) * tx2
3684 dumi(i,k) = dumi(i,k) - faltndi * tx3
3685 dumni(i,k) = dumni(i,k) - faltndni * tx3
3686 dumc(i,k) = dumc(i,k) - faltndc * tx3
3687 dumnc(i,k) = dumnc(i,k) - faltndnc * tx3
3689 fni(k) = max(fni(k)*tx1, fni(k-1)*tx4) * pdel(i,k)
3690 fi(k) = max(fi(k)*tx1, fi(k-1)*tx4) * pdel(i,k)
3691 fnc(k) = max(fnc(k)*tx1, fnc(k-1)*tx4) * pdel(i,k)
3692 fc(k) = max(fc(k)*tx1, fc(k-1)*tx4) * pdel(i,k)
3702 tx3 = tx2 / (
g*1000._r8)
3703 prect(i) = prect(i) + (faloutc(pver)+falouti(pver)) * tx3
3704 preci(i) = preci(i) + falouti(pver) * tx3
3720 if (fprcp == 1)
then
3732 dumc(i,k) = max(qrn(i,k)+qrtend(i,k)*deltat, zero)
3733 dumi(i,k) = max(qsnw(i,k)+qnitend(i,k)*deltat,zero)
3734 dumnc(i,k) = max(nrn(i,k)+nrtend(i,k)*deltat, zero)
3735 dumni(i,k) = max(nsnw(i,k)+nstend(i,k)*deltat, zero)
3739 if (dumi(i,k) <
qsmall)
then
3744 if (dumc(i,k) <
qsmall)
then
3752 dumnc(i,k) = max(dumnc(i,k),zero)
3753 dumni(i,k) = max(dumni(i,k),zero)
3760 if (dumc(i,k) >=
qsmall)
then
3761 lamr(k) = (pi*
rhow*dumnc(i,k)/dumc(i,k))**
oneb3
3762 n0r(k) = dumnc(i,k)*lamr(k)
3772 n0r(k) = lamr(k)**4*dumc(i,k)/(pi*
rhow)
3773 dumnc(i,k) = n0r(k)/lamr(k)
3774 else if (lamr(k) >
lammaxr)
then
3776 n0r(k) = lamr(k)**4*dumc(i,k)/(pi*
rhow)
3777 dumnc(i,k) = n0r(k)/lamr(k)
3783 tx1 = arn(i,k) / lamr(k)**
br
3784 tx2 = 9.1_r8*rhof(i,k)
3785 unr(k) = min(tx1*
cons4, tx2)
3786 umr(k) = min(tx1*(
cons5/6._r8),tx2)
3798 if (dumi(i,k) >=
qsmall)
then
3799 lams(k) = (cons6*
cs*dumni(i,k)/ dumi(i,k))**(
one/
ds)
3800 n0s(k) = dumni(i,k)*lams(k)
3805 if (lams(k) < lammins)
then
3807 n0s(k) = lams(k)**(
ds+
one)*dumi(i,k)/(
cs*cons6)
3808 dumni(i,k) = n0s(k)/lams(k)
3809 else if (lams(k) > lammaxs)
then
3811 n0s(k) = lams(k)**(
ds+
one)*dumi(i,k)/(
cs*cons6)
3812 dumni(i,k) = n0s(k)/lams(k)
3817 tx1 = asn(i,k) / lams(k)**
bs
3818 tx2 = 1.2_r8*rhof(i,k)
3819 ums(k) = min(tx1*(cons8/6._r8), tx2)
3820 uns(k) = min(tx1*cons7, tx2)
3838 rgvm = max(fi(k),fc(k),fni(k),fnc(k))
3839 nstep = max(int(rgvm*deltat/pdel(i,k)+
one),nstep)
3844 qrn(i,k) = (qrn(i,k) + qrtend(i,k)*deltat)
3845 qsnw(i,k) = (qsnw(i,k) + qnitend(i,k)*deltat)
3846 nrn(i,k) = max((nrn(i,k) + nrtend(i,k)*deltat),zero)
3847 nsnw(i,k) = max((nsnw(i,k) + nstend(i,k)*deltat),zero)
3849 if (qrn(i,k) <
qsmall) nrn(i,k) = zero
3850 if (qsnw(i,k) <
qsmall) nsnw(i,k) = zero
3854 tx2 =
one / float(nstep)
3859 falouti(k) = max(fi(k) * qsnw(i,k), zero)
3860 faloutni(k) = max(fni(k) * nsnw(i,k), zero)
3861 faloutc(k) = max(fc(k) * qrn(i,k), zero)
3862 faloutnc(k) = max(fnc(k) * nrn(i,k), zero)
3868 tx1 =
one / pdel(i,k)
3869 faltndi = falouti(k) * tx1
3870 faltndni = faloutni(k) * tx1
3871 faltndc = faloutc(k) * tx1
3872 faltndnc = faloutnc(k) * tx1
3883 qsnw(i,k) = qsnw(i,k) - faltndi * tx3
3884 nsnw(i,k) = nsnw(i,k) - faltndni * tx3
3885 qrn(i,k) = qrn(i,k) - faltndc * tx3
3886 nrn(i,k) = nrn(i,k) - faltndnc * tx3
3891 tx1 =
one / pdel(i,k)
3895 faltndc = (faloutc(k) - faloutc(k-1)) * tx1
3896 faltndnc = (faloutnc(k) - faloutnc(k-1)) * tx1
3898 faltndi = (falouti(k) - falouti(k-1)) * tx1
3899 faltndni = (faloutni(k) - faloutni(k-1)) * tx1
3901 qsnw(i,k) = qsnw(i,k) - faltndi * tx3
3902 nsnw(i,k) = nsnw(i,k) - faltndni * tx3
3903 qrn(i,k) = qrn(i,k) - faltndc * tx3
3904 nrn(i,k) = nrn(i,k) - faltndnc * tx3
3909 tx5 = tx2 / (
g*1000._r8)
3910 prect(i) = prect(i) + (faloutc(pver)+falouti(pver)) * tx5
3911 preci(i) = preci(i) + (falouti(pver)) * tx5
3922 if (fprcp == 1)
then
3927 tx1 = t(i,k) + tlat(i,k)*
onebcp*deltat
3928 if (tx1 > 275.15_r8)
then
3929 if (qsnw(i,k) > zero)
then
3934 if (tx1-dum < 273.15_r8)
then
3936 dum = min(
one,max(zero,tx2/dum))
3941 qrn(i,k) = qrn(i,k) + dum*qsnw(i,k)
3942 nrn(i,k) = nrn(i,k) + dum*nsnw(i,k)
3943 qsnw(i,k) = (
one-dum)*qsnw(i,k)
3944 nsnw(i,k) = (
one-dum)*nsnw(i,k)
3946 tmp = -
xlf*dum*qsnw(i,k)/deltat
3948 tlat(i,k) = tlat(i,k) + tmp
3950 preci(i) = (
one-dum)*preci(i)
3956 tx1 = t(i,k) + tlat(i,k)*
onebcp*deltat
3959 if (qrn(i,k) > zero)
then
3964 tx2 = -(tx1-(
tmelt-5._r8))
3965 dum = min(
one,max(zero,tx2/dum))
3970 qsnw(i,k) = qsnw(i,k) + dum*qrn(i,k)
3971 nsnw(i,k) = nsnw(i,k) + dum*nrn(i,k)
3972 qrn(i,k) = (
one-dum)*qrn(i,k)
3973 nrn(i,k) = (
one-dum)*nrn(i,k)
3975 tmp =
xlf*dum*qrn(i,k)/deltat
3977 tlat(i,k) = tlat(i,k) + tmp
3979 preci(i) = preci(i) + dum*(prect(i)-preci(i))
3985 if (qsnw(i,k) <
qsmall)
then
3989 if (qrn(i,k) <
qsmall)
then
3996 nrn(i,k) = max(nrn(i,k),zero)
3997 nsnw(i,k) = max(nsnw(i,k),zero)
3999 qrout(i,k) = qrout(i,k) + qrn(i,k)
4000 qsout(i,k) = qsout(i,k) + qsnw(i,k)
4001 nrout(i,k) = nrout(i,k) + nrn(i,k)*rho(i,k)
4002 nsout(i,k) = nsout(i,k) + nsnw(i,k)*rho(i,k)
4008 dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat, zero)
4009 dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat, zero)
4010 dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat, zero)
4011 dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat, zero)
4015 dumnc(i,k) = ncnst * irho(i,k) * lcldm(i,k)
4019 dumni(i,k) = ninst * irho(i,k) * icldm(i,k)
4022 if (dumc(i,k) <
qsmall) dumnc(i,k) = zero
4023 if (dumi(i,k) <
qsmall) dumni(i,k) = zero
4027 tx1 = t(i,k) + tlat(i,k)*
onebcp*deltat
4028 if (tx1 >
tmelt)
then
4029 if (dumi(i,k) > zero)
then
4033 if (tx1-dum <
tmelt)
then
4035 dum = max(zero, min(
one, dum))
4040 tx2 = dum*dumi(i,k)*dti
4041 qctend(i,k) = qctend(i,k) + tx2
4048 nctend(i,k) = nctend(i,k) +
three*tx2
4051 qitend(i,k) = ((
one-dum)*dumi(i,k)-qi(i,k)) * dti
4052 nitend(i,k) = ((
one-dum)*dumni(i,k)-ni(i,k)) * dti
4053 tlat(i,k) = tlat(i,k) -
xlf*tx2
4060 tx1 = t(i,k) + tlat(i,k)*
onebcp*deltat
4061 if (tx1 < 233.15_r8)
then
4063 if (dumc(i,k) > zero .and. qc(i,k) > zero)
then
4067 if (tx1+dum > 233.15_r8)
then
4068 dum = -(tx1-233.15_r8)*
cpoxlf / dumc(i,k)
4069 dum = max(zero, min(
one, dum))
4074 tx2 = dum*dumc(i,k)*dti
4075 qitend(i,k) = qitend(i,k) + tx2
4082 nitend(i,k) = nitend(i,k) + dum*
three*dumc(i,k)
4083 & / (
four*3.14_r8*1.563e-14_r8* 500._r8) * dti
4084 qctend(i,k) = ((
one-dum)*dumc(i,k)-qc(i,k)) * dti
4085 nctend(i,k) = ((
one-dum)*dumnc(i,k)-nc(i,k)) * dti
4086 tlat(i,k) = tlat(i,k) +
xlf*tx2
4130 tx1 =
one / lcldm(i,k)
4131 tx2 =
one / icldm(i,k)
4132 dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat, zero) * tx1
4133 dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat, zero) * tx2
4134 dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat, zero) * tx1
4135 dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat, zero) * tx2
4146 dumnc(i,k) = ncnst * irho(i,k)
4150 dumni(i,k) = ninst * irho(i,k)
4156 dumc(i,k) = min(dumc(i,k),5.e-3_r8)
4157 dumi(i,k) = min(dumi(i,k),5.e-3_r8)
4162 if (dumi(i,k) >=
qsmall)
then
4164 dumni(i,k) = min(dumni(i,k),dumi(i,k)*1.e20_r8)
4168 miu_ice(k) = max(min(0.008_r8*(lami(k)*0.01)**0.87_r8,
4170 tx1 =
one + miu_ice(k)
4173 lami(k) = aux*lami(k)
4174 n0i(k) = niic(i,k) * lami(k)**tx1 * tx2
4176 if (lami(k) <
lammini*aux)
then
4180 niic(i,k) = n0i(k)/lami(k)
4182 nitend(i,k) = (niic(i,k)*icldm(i,k)-ni(i,k))/deltat
4183 else if (lami(k) >
lammaxi*aux)
then
4187 niic(i,k) = n0i(k)/lami(k)
4188 nitend(i,k) = (niic(i,k)*icldm(i,k)-ni(i,k))/deltat
4191 effi(i,k) = 1.e6_r8*
gamma(
four+miu_ice(k))
4201 if (dumc(i,k) >=
qsmall)
then
4204 dumnc(i,k) = min(dumnc(i,k),dumc(i,k)*1.e20_r8)
4210 if (dumnc(i,k) < cdnl*irho(i,k))
then
4211 nctend(i,k) = (cdnl*irho(i,k)*cldm(i,k)-nc(i,k))*dti
4213 dumnc(i,k) = max(dumnc(i,k),cdnl*irho(i,k))
4221 nctend(i,k) = (ncnst*irho(i,k)*lcldm(i,k)-nc(i,k))
4227 tx1 = 0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
4228 pgam(k) = max(
two, min(15._r8,
one/(tx1*tx1)-
one))
4232 lamc(k) = (
cr*dumnc(i,k)*tx2
4233 & / (dumc(i,k)*tx1))**
oneb3
4234 lammin = (pgam(k)+
one) / 50.e-6_r8
4235 lammax = (pgam(k)+
one) / 2.e-6_r8
4236 if (lamc(k) < lammin)
then
4238 ncic(i,k) = 6._r8*lamc(k)*lamc(k)*lamc(k)*dumc(i,k)
4241 nctend(i,k) = (ncic(i,k)*lcldm(i,k)-nc(i,k))*dti
4243 else if (lamc(k) > lammax)
then
4245 ncic(i,k) = 6._r8*lamc(k)*lamc(k)*lamc(k)*dumc(i,k)
4248 nctend(i,k) = (ncic(i,k)*lcldm(i,k)-nc(i,k))*dti
4251 effc(i,k) = tx2 / (
gamma(pgam(k)+
three)*lamc(k)*2.e6_r8)
4253 lamcrad(i,k) = lamc(k)
4254 pgamrad(i,k) = pgam(k)
4264 deffi(i,k) = effi(i,k) * (
rhoi / 917._r8*
two)
4273 if (dumc(i,k) >=
qsmall)
then
4274 tx1 = 0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
4275 pgam(k) = max(
two, min(15._r8,
one/(tx1*tx1)-
one))
4277 lamc(k) = (
cr*dumnc(i,k)*tx2
4279 lammin = (pgam(k)+
one) / 50.e-6_r8
4280 lammax = (pgam(k)+
one) / 2.e-6_r8
4281 if (lamc(k) < lammin)
then
4283 else if (lamc(k) > lammax)
then
4286 effc_fn(i,k) = tx2/(
gamma(pgam(k)+
three)*lamc(k)*2.e6_r8)
4289 effc_fn(i,k) = 10._r8
4302 deltat = deltat*real(iter)
4309 if (qc(i,k)+qctend(i,k)*deltat <
qsmall)
4310 & nctend(i,k) = -nc(i,k) * dti
4311 if (qi(i,k)+qitend(i,k)*deltat <
qsmall)
4312 & nitend(i,k) = -ni(i,k) * dti
4324 call outfld(
'QRAIN',qrout, pcols, lchnk)
4325 call outfld(
'QSNOW',qsout, pcols, lchnk)
4326 call outfld(
'NRAIN',nrout, pcols, lchnk)
4327 call outfld(
'NSNOW',nsout, pcols, lchnk)
4333 if (qsout(i,k) > 1.e-7_r8 .and. nsout(i,k) > zero)
then
4334 dsout(i,k) = (
pirhosn * nsout(i,k)/qsout(i,k)) **(-
oneb3)
4340 call outfld(
'DSNOW',dsout, pcols, lchnk)
4343 call outfld(
'MGFLXPRC',rflx, pcols, lchnk)
4344 call outfld(
'MGFLXSNW',
sflx, pcols, lchnk)
4349 if ((qc(i,k)+qctend(i,k)*deltat >=
qsmall) .and.
4350 & (cldmax(i,k) > zero) .and. (lcldm(i, k) > zero) .and.
4351 & (nc(i,k)+nctend(i,k)*deltat > 0.0))
then
4353 tx1 = rho(i,k) / lcldm(i,k)
4354 tx2 = (qc(i,k)+qctend(i,k)*deltat)*tx1*1000._r8
4355 dum = tx2 * tx2 * lcldm(i,k)
4356 & / (0.109_r8*(nc(i,k)+nctend(i,k)*deltat)
4357 & *tx1*1.e-6_r8*cldmax(i,k))
4361 if (qi(i,k)+qitend(i,k)*deltat >=
qsmall)
then
4362 dum1 = ((qi(i,k)+qitend(i,k)*deltat)*rho(i,k)/icldm(i,k)
4363 & * 1000._r8/0.1_r8)**(
one/0.63_r8)
4364 & * icldm(i,k)/cldmax(i,k)
4369 if (qsout(i,k) >=
qsmall)
then
4370 dum1 = dum1 + (qsout(i,k)*rho(i,k)*1000._r8/0.1_r8)
4374 refl(i,k) = dum + dum1
4376 if (rainrt(i,k) >= 0.001)
then
4377 dum = log10(rainrt(i,k)**6._r8) + 16._r8
4378 dum = 10._r8**(dum/10._r8)
4383 refl(i,k) = refl(i,k) + dum
4385 areflz(i,k) = refl(i,k)
4388 refl(i,k) = 10._r8*log10(refl(i,k))
4390 refl(i,k) = -9999._r8
4393 if (refl(i,k) >
mindbz)
then
4394 arefl(i,k) = refl(i,k)
4402 csrfl(i,k) = min(
csmax,refl(i,k))
4404 if (csrfl(i,k) >
csmin)
then
4405 acsrfl(i,k) = refl(i,k)
4416 call outfld(
'REFL',refl, pcols, lchnk)
4417 call outfld(
'AREFL',arefl, pcols, lchnk)
4418 call outfld(
'AREFLZ',areflz, pcols, lchnk)
4419 call outfld(
'FREFL',frefl, pcols, lchnk)
4420 call outfld(
'CSRFL',csrfl, pcols, lchnk)
4421 call outfld(
'ACSRFL',acsrfl, pcols, lchnk)
4422 call outfld(
'FCSRFL',fcsrfl, pcols, lchnk)
4424 call outfld(
'RERCLD',rercld, pcols, lchnk)
4441 if (qrout(i,k) > 1.e-7_r8 .and. nrout(i,k) > zero)
then
4442 qrout2(i,k) = qrout(i,k)
4443 nrout2(i,k) = nrout(i,k)
4444 drout2(i,k) = (
pirhow*nrout(i,k) / qrout(i,k))**(-
oneb3)
4447 if (qsout(i,k) > 1.e-7_r8 .and. nsout(i,k) > zero)
then
4448 qsout2(i,k) = qsout(i,k)
4449 nsout2(i,k) = nsout(i,k)
4450 dsout2(i,k) = (
pirhosn*nsout(i,k) / qsout(i,k))**(-
oneb3)
4455 ncai(i,k) = dum2i(i,k)*rho(i,k)
4456 ncal(i,k) = dum2l(i,k)*rho(i,k)
4462 call outfld(
'NCAL',ncal, pcols,lchnk)
4463 call outfld(
'NCAI',ncai, pcols,lchnk)
4466 call outfld(
'AQRAIN',qrout2, pcols,lchnk)
4467 call outfld(
'AQSNOW',qsout2, pcols,lchnk)
4468 call outfld(
'ANRAIN',nrout2, pcols,lchnk)
4469 call outfld(
'ANSNOW',nsout2, pcols,lchnk)
4470 call outfld(
'ADRAIN',drout2, pcols,lchnk)
4471 call outfld(
'ADSNOW',dsout2, pcols,lchnk)
4472 call outfld(
'FREQR',freqr, pcols,lchnk)
4473 call outfld(
'FREQS',freqs, pcols,lchnk)
4480 tx1 = qc(i,k) + qctend(i,k)*deltat
4481 tx2 = qi(i,k) + qitend(i,k)*deltat
4482 dumfice = qsout(i,k) + qrout(i,k) + tx1 + tx2
4483 if (dumfice > zero)
then
4484 nfice(i,k) = (qsout(i,k) + tx2) / dumfice
4490 call outfld(
'FICE',nfice, pcols, lchnk)