1004 nwfa, nifa, nwfa2d, nifa2d, &
1006 p, w, dz, dt_in, dt_inner, &
1007 sedi_semi, decfl, lsm, &
1011 GRAUPELNC, GRAUPELNCV, SR, &
1012#if ( WRF_CHEM == 1 )
1013 rainprod, evapprod, &
1015 refl_10cm, diagflag, do_radar_ref, &
1016 max_hail_diam_sfc, &
1017 vt_dbz_wt, first_time_step, &
1018 re_cloud, re_ice, re_snow, &
1019 has_reqc, has_reqi, has_reqs, &
1020 aero_ind_fdb, rand_perturb_on, &
1022 rand_pert, spp_prt_list, spp_var_list, &
1023 spp_stddev_cutoff, n_var_spp, &
1024 ids,ide, jds,jde, kds,kde, &
1025 ims,ime, jms,jme, kms,kme, &
1026 its,ite, jts,jte, kts,kte, &
1027 fullradar_diag, istep, nsteps, &
1034 prw_vcde, tpri_inu, tpri_ide_d, &
1035 tpri_ide_s, tprs_ide, tprs_sde_d, &
1036 tprs_sde_s, tprg_gde_d, &
1037 tprg_gde_s, tpri_iha, tpri_wfz, &
1038 tpri_rfz, tprg_rfz, tprs_scw, tprg_scw, &
1039 tprg_rcs, tprs_rcs, &
1040 tprr_rci, tprg_rcg, &
1041 tprw_vcd_c, tprw_vcd_e, tprr_sml, &
1042 tprr_gml, tprr_rcg, &
1043 tprr_rcs, tprv_rev, tten3, qvten3, &
1044 qrten3, qsten3, qgten3, qiten3, niten3, &
1045 nrten3, ncten3, qcten3, &
1051 integer,
intent(in):: ids,ide, jds,jde, kds,kde, &
1052 ims,ime, jms,jme, kms,kme, &
1053 its,ite, jts,jte, kts,kte
1054 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
intent(inout):: &
1055 qv, qc, qr, qi, qs, qg, ni, nr
1056 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
optional,
intent(inout):: &
1058 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
optional,
intent(in):: &
1060 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
optional,
intent(inout):: &
1062 real(wp),
dimension(ims:ime, jms:jme),
optional,
intent(in):: nwfa2d, nifa2d
1063 integer,
dimension(ims:ime, jms:jme),
intent(in):: lsm
1064 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
optional,
intent(inout):: &
1065 re_cloud, re_ice, re_snow
1066 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
intent(inout):: pfils, pflls
1067 integer,
intent(in) :: rand_perturb_on, kme_stoch, n_var_spp
1068 real(wp),
dimension(:,:),
intent(in),
optional :: rand_pert
1069 real(wp),
dimension(:),
intent(in),
optional :: spp_prt_list, spp_stddev_cutoff
1070 character(len=10),
dimension(:),
intent(in),
optional :: spp_var_list
1071 integer,
intent(in):: has_reqc, has_reqi, has_reqs
1072#if ( WRF_CHEM == 1 )
1073 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
intent(inout):: &
1076 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
intent(in):: &
1078 real(wp),
dimension(ims:ime, jms:jme),
intent(inout):: &
1080 real(wp),
dimension(ims:ime, jms:jme),
optional,
intent(inout):: &
1083 GRAUPELNC, GRAUPELNCV
1084 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
intent(inout):: &
1086 real(wp),
dimension(ims:ime, jms:jme),
intent(inout):: &
1088 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
optional,
intent(inout):: &
1090 logical,
intent(in) :: first_time_step
1091 real(wp),
intent(in):: dt_in, dt_inner
1092 logical,
intent(in) :: sedi_semi
1093 integer,
intent(in) :: decfl
1095 integer,
intent (in) :: istep, nsteps
1096 logical,
intent (in) :: fullradar_diag
1098 logical,
intent (in) :: ext_diag
1099 logical,
optional,
intent(in):: aero_ind_fdb
1100 real(wp),
dimension(:,:,:),
optional,
intent(inout):: &
1103 prw_vcde, tpri_inu, tpri_ide_d, &
1104 tpri_ide_s, tprs_ide, &
1105 tprs_sde_d, tprs_sde_s, tprg_gde_d, &
1106 tprg_gde_s, tpri_iha, tpri_wfz, &
1107 tpri_rfz, tprg_rfz, tprs_scw, tprg_scw, &
1108 tprg_rcs, tprs_rcs, &
1109 tprr_rci, tprg_rcg, &
1110 tprw_vcd_c, tprw_vcd_e, tprr_sml, &
1111 tprr_gml, tprr_rcg, &
1112 tprr_rcs, tprv_rev, tten3, qvten3, &
1113 qrten3, qsten3, qgten3, qiten3, niten3, &
1114 nrten3, ncten3, qcten3
1117 real(wp),
dimension(kts:kte):: &
1118 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
1119 nr1d, nc1d, nwfa1d, nifa1d, &
1120 t1d, p1d, w1d, dz1d, rho, dBZ, pfil1, pfll1
1122 real(wp),
dimension(:),
allocatable:: &
1125 prw_vcde1, tpri_inu1, tpri_ide1_d, &
1126 tpri_ide1_s, tprs_ide1, &
1127 tprs_sde1_d, tprs_sde1_s, tprg_gde1_d, &
1128 tprg_gde1_s, tpri_iha1, tpri_wfz1, &
1129 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1,&
1130 tprg_rcs1, tprs_rcs1, &
1131 tprr_rci1, tprg_rcg1, &
1132 tprw_vcd1_c, tprw_vcd1_e, tprr_sml1, &
1133 tprr_gml1, tprr_rcg1, &
1134 tprr_rcs1, tprv_rev1, tten1, qvten1, &
1135 qrten1, qsten1, qgten1, qiten1, niten1, &
1136 nrten1, ncten1, qcten1
1138 real(wp),
dimension(kts:kte):: re_qc1d, re_qi1d, re_qs1d
1139#if ( WRF_CHEM == 1 )
1140 real(wp),
dimension(kts:kte):: &
1141 rainprod1d, evapprod1d
1143 real(wp),
dimension(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
1144 real(wp) :: dt, pptrain, pptsnow, pptgraul, pptice
1145 real(wp) :: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max
1147 real(wp) :: rand1, rand2, rand3, rand_pert_max
1148 integer:: i, j, k, m
1149 integer:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr
1150 integer:: jmax_qc,jmax_qr,jmax_qi,jmax_qs,jmax_qg,jmax_ni,jmax_nr
1151 integer:: kmax_qc,kmax_qr,kmax_qi,kmax_qs,kmax_qg,kmax_ni,kmax_nr
1152 integer:: i_start, j_start, i_end, j_end
1153 logical,
optional,
intent(in) :: diagflag
1154 integer,
optional,
intent(in) :: do_radar_ref
1155 logical :: melti = .false.
1159 character(len=*),
optional,
intent( out) :: errmsg
1160 integer,
optional,
intent( out) :: errflg
1163 if (
present(errmsg)) errmsg =
''
1164 if (
present(errflg)) errflg = 0
1167 test_only_once:
if (first_time_step .and. istep==1)
then
1170 if ( (
present(tt) .and. (
present(th) .or.
present(pii))) .or. &
1171 (.not.
present(tt) .and. .not.(
present(th) .and.
present(pii))) )
then
1172 if (
present(errmsg) .and.
present(errflg))
then
1173 write(errmsg,
'(a)')
'Logic error in mp_gt_driver: provide either tt or th+pii'
1177 write(*,
'(a)')
'Logic error in mp_gt_driver: provide either tt or th+pii'
1182 if (is_aerosol_aware .and. (.not.
present(nc) .or. &
1183 .not.
present(nwfa) .or. &
1184 .not.
present(nifa) .or. &
1185 .not.
present(nwfa2d) .or. &
1186 .not.
present(nifa2d) ))
then
1187 if (
present(errmsg) .and.
present(errflg))
then
1188 write(errmsg,
'(*(a))')
'Logic error in mp_gt_driver: provide nc, nwfa, nifa, nwfa2d', &
1189 ' and nifa2d for aerosol-aware version of Thompson microphysics'
1193 write(*,
'(*(a))')
'Logic error in mp_gt_driver: provide nc, nwfa, nifa, nwfa2d', &
1194 ' and nifa2d for aerosol-aware version of Thompson microphysics'
1197 else if (merra2_aerosol_aware .and. (.not.
present(nc) .or. &
1198 .not.
present(nwfa) .or. &
1199 .not.
present(nifa) ))
then
1200 if (
present(errmsg) .and.
present(errflg))
then
1201 write(errmsg,
'(*(a))')
'Logic error in mp_gt_driver: provide nc, nwfa, and nifa', &
1202 ' for merra2 aerosol-aware version of Thompson microphysics'
1206 write(*,
'(*(a))')
'Logic error in mp_gt_driver: provide nc, nwfa, and nifa', &
1207 ' for merra2 aerosol-aware version of Thompson microphysics'
1210 else if (.not.is_aerosol_aware .and. .not.merra2_aerosol_aware .and. &
1211 (
present(nwfa) .or.
present(nifa) .or.
present(nwfa2d) .or.
present(nifa2d)))
then
1212 write(*,*)
'WARNING, nc/nwfa/nifa/nwfa2d/nifa2d present but is_aerosol_aware/merra2_aerosol_aware are FALSE'
1214 end if test_only_once
1220 allocate_extended_diagnostics:
if (ext_diag)
then
1221 allocate (prw_vcdc1(kts:kte))
1222 allocate (prw_vcde1(kts:kte))
1223 allocate (tpri_inu1(kts:kte))
1224 allocate (tpri_ide1_d(kts:kte))
1225 allocate (tpri_ide1_s(kts:kte))
1226 allocate (tprs_ide1(kts:kte))
1227 allocate (tprs_sde1_d(kts:kte))
1228 allocate (tprs_sde1_s(kts:kte))
1229 allocate (tprg_gde1_d(kts:kte))
1230 allocate (tprg_gde1_s(kts:kte))
1231 allocate (tpri_iha1(kts:kte))
1232 allocate (tpri_wfz1(kts:kte))
1233 allocate (tpri_rfz1(kts:kte))
1234 allocate (tprg_rfz1(kts:kte))
1235 allocate (tprs_scw1(kts:kte))
1236 allocate (tprg_scw1(kts:kte))
1237 allocate (tprg_rcs1(kts:kte))
1238 allocate (tprs_rcs1(kts:kte))
1239 allocate (tprr_rci1(kts:kte))
1240 allocate (tprg_rcg1(kts:kte))
1241 allocate (tprw_vcd1_c(kts:kte))
1242 allocate (tprw_vcd1_e(kts:kte))
1243 allocate (tprr_sml1(kts:kte))
1244 allocate (tprr_gml1(kts:kte))
1245 allocate (tprr_rcg1(kts:kte))
1246 allocate (tprr_rcs1(kts:kte))
1247 allocate (tprv_rev1(kts:kte))
1248 allocate (tten1(kts:kte))
1249 allocate (qvten1(kts:kte))
1250 allocate (qrten1(kts:kte))
1251 allocate (qsten1(kts:kte))
1252 allocate (qgten1(kts:kte))
1253 allocate (qiten1(kts:kte))
1254 allocate (niten1(kts:kte))
1255 allocate (nrten1(kts:kte))
1256 allocate (ncten1(kts:kte))
1257 allocate (qcten1(kts:kte))
1258 end if allocate_extended_diagnostics
1279 graupelnc(:,:) = 0.0
1287 ndt = max(nint(dt_in/dt_inner),1)
1289 if(dt_in .le. dt_inner) dt= dt_in
1294 if (rand_perturb_on .ne. 0)
then
1296 select case (spp_var_list(k))
1298 rand_pert_max = spp_prt_list(k)*spp_stddev_cutoff(k)
1334 j_loop:
do j = j_start, j_end
1335 i_loop:
do i = i_start, i_end
1353 if (rand_perturb_on .ne. 0)
then
1354 if (mod(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,1)
1355 m = rshift(abs(rand_perturb_on),1)
1356 if (mod(m,2) .ne. 0) rand2 = rand_pert(i,1)*2.
1357 m = rshift(abs(rand_perturb_on),2)
1358 if (mod(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,1)+rand_pert_max)
1359 m = rshift(abs(rand_perturb_on),3)
1368 IF (
PRESENT (snowncv) )
THEN
1371 IF (
PRESENT (icencv) )
THEN
1374 IF (
PRESENT (graupelncv) )
THEN
1375 graupelncv(i,j) = 0.
1380 if (
present(tt))
then
1383 t1d(k) = th(i,k,j)*pii(i,k,j)
1396 rho(k) = roverrv*p1d(k) / (r*t1d(k)*(qv1d(k)+roverrv))
1402 initialize_extended_diagnostics:
if (ext_diag)
then
1440 endif initialize_extended_diagnostics
1444 if (is_aerosol_aware .or. merra2_aerosol_aware)
then
1447 nwfa1d(k) = nwfa(i,k,j)
1448 nifa1d(k) = nifa(i,k,j)
1453 nc1d(k) = nt_c_l/rho(k)
1455 nc1d(k) = nt_c_o/rho(k)
1458 nifa1d(k) = nain1*0.01
1463 call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
1464 nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dz1d, &
1465 lsml, pptrain, pptsnow, pptgraul, pptice, &
1466#if ( WRF_CHEM == 1 )
1467 rainprod1d, evapprod1d, &
1469 rand1, rand2, rand3, &
1470 kts, kte, dt, i, j, ext_diag, &
1473 prw_vcdc1, prw_vcde1, &
1474 tpri_inu1, tpri_ide1_d, tpri_ide1_s, tprs_ide1, &
1475 tprs_sde1_d, tprs_sde1_s, &
1476 tprg_gde1_d, tprg_gde1_s, tpri_iha1, tpri_wfz1, &
1477 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1, &
1478 tprg_rcs1, tprs_rcs1, tprr_rci1, &
1479 tprg_rcg1, tprw_vcd1_c, &
1480 tprw_vcd1_e, tprr_sml1, tprr_gml1, tprr_rcg1, &
1481 tprr_rcs1, tprv_rev1, &
1482 tten1, qvten1, qrten1, qsten1, &
1483 qgten1, qiten1, niten1, nrten1, ncten1, qcten1, &
1486 pcp_ra(i,j) = pcp_ra(i,j) + pptrain
1487 pcp_sn(i,j) = pcp_sn(i,j) + pptsnow
1488 pcp_gr(i,j) = pcp_gr(i,j) + pptgraul
1489 pcp_ic(i,j) = pcp_ic(i,j) + pptice
1490 rainncv(i,j) = pptrain + pptsnow + pptgraul + pptice
1491 rainnc(i,j) = rainnc(i,j) + pptrain + pptsnow + pptgraul + pptice
1492 IF (
PRESENT(snowncv) .AND.
PRESENT(snownc) )
THEN
1494 IF ( .NOT.
PRESENT(icencv) .OR. .NOT.
PRESENT(icenc) )
THEN
1495 snowncv(i,j) = pptsnow + pptice
1496 snownc(i,j) = snownc(i,j) + pptsnow + pptice
1498 snowncv(i,j) = pptsnow
1499 snownc(i,j) = snownc(i,j) + pptsnow
1503 IF (
PRESENT(icencv) .AND.
PRESENT(icenc) )
THEN
1504 icencv(i,j) = pptice
1505 icenc(i,j) = icenc(i,j) + pptice
1507 IF (
PRESENT(graupelncv) .AND.
PRESENT(graupelnc) )
THEN
1508 graupelncv(i,j) = pptgraul
1509 graupelnc(i,j) = graupelnc(i,j) + pptgraul
1511 sr(i,j) = (pptsnow + pptgraul + pptice) / (rainncv(i,j)+r1)
1516 if (is_aerosol_aware)
then
1517 if (
PRESENT (aero_ind_fdb) )
then
1518 if ( .not. aero_ind_fdb)
then
1519 nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt
1520 nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt
1523 nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt
1524 nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt
1529 nwfa(i,k,j) = nwfa1d(k)
1530 nifa(i,k,j) = nifa1d(k)
1534 if (merra2_aerosol_aware)
then
1537 nwfa(i,k,j) = nwfa1d(k)
1538 nifa(i,k,j) = nifa1d(k)
1551 pfils(i,k,j) = pfils(i,k,j) + pfil1(k)
1552 pflls(i,k,j) = pflls(i,k,j) + pfll1(k)
1553 if (
present(tt))
then
1556 th(i,k,j) = t1d(k)/pii(i,k,j)
1558#if ( WRF_CHEM == 1 )
1559 rainprod(i,k,j) = rainprod1d(k)
1560 evapprod(i,k,j) = evapprod1d(k)
1562 if (qc1d(k) .gt. qc_max)
then
1567 elseif (qc1d(k) .lt. 0.0)
then
1568 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qc ', qc1d(k), &
1571 if (qr1d(k) .gt. qr_max)
then
1576 elseif (qr1d(k) .lt. 0.0)
then
1577 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qr ', qr1d(k), &
1580 if (nr1d(k) .gt. nr_max)
then
1585 elseif (nr1d(k) .lt. 0.0)
then
1586 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative nr ', nr1d(k), &
1589 if (qs1d(k) .gt. qs_max)
then
1594 elseif (qs1d(k) .lt. 0.0)
then
1595 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qs ', qs1d(k), &
1598 if (qi1d(k) .gt. qi_max)
then
1603 elseif (qi1d(k) .lt. 0.0)
then
1604 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qi ', qi1d(k), &
1607 if (qg1d(k) .gt. qg_max)
then
1612 elseif (qg1d(k) .lt. 0.0)
then
1613 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qg ', qg1d(k), &
1616 if (ni1d(k) .gt. ni_max)
then
1621 elseif (ni1d(k) .lt. 0.0)
then
1622 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative ni ', ni1d(k), &
1625 if (qv1d(k) .lt. 0.0)
then
1626 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qv ', qv1d(k), &
1628 if (k.lt.kte-2 .and. k.gt.kts+1)
then
1629 write(*,*)
' below and above are: ', qv(i,k-1,j), qv(i,k+1,j)
1630 qv(i,k,j) = max(1.e-7, 0.5*(qv(i,k-1,j) + qv(i,k+1,j)))
1637 assign_extended_diagnostics:
if (ext_diag)
then
1642 prw_vcdc(i,k,j) = prw_vcdc(i,k,j) + prw_vcdc1(k)
1643 prw_vcde(i,k,j) = prw_vcde(i,k,j) + prw_vcde1(k)
1644 tpri_inu(i,k,j) = tpri_inu(i,k,j) + tpri_inu1(k)
1645 tpri_ide_d(i,k,j) = tpri_ide_d(i,k,j) + tpri_ide1_d(k)
1646 tpri_ide_s(i,k,j) = tpri_ide_s(i,k,j) + tpri_ide1_s(k)
1647 tprs_ide(i,k,j) = tprs_ide(i,k,j) + tprs_ide1(k)
1648 tprs_sde_s(i,k,j) = tprs_sde_s(i,k,j) + tprs_sde1_s(k)
1649 tprs_sde_d(i,k,j) = tprs_sde_d(i,k,j) + tprs_sde1_d(k)
1650 tprg_gde_d(i,k,j) = tprg_gde_d(i,k,j) + tprg_gde1_d(k)
1651 tprg_gde_s(i,k,j) = tprg_gde_s(i,k,j) + tprg_gde1_s(k)
1652 tpri_iha(i,k,j) = tpri_iha(i,k,j) + tpri_iha1(k)
1653 tpri_wfz(i,k,j) = tpri_wfz(i,k,j) + tpri_wfz1(k)
1654 tpri_rfz(i,k,j) = tpri_rfz(i,k,j) + tpri_rfz1(k)
1655 tprg_rfz(i,k,j) = tprg_rfz(i,k,j) + tprg_rfz1(k)
1656 tprs_scw(i,k,j) = tprs_scw(i,k,j) + tprs_scw1(k)
1657 tprg_scw(i,k,j) = tprg_scw(i,k,j) + tprg_scw1(k)
1658 tprg_rcs(i,k,j) = tprg_rcs(i,k,j) + tprg_rcs1(k)
1659 tprs_rcs(i,k,j) = tprs_rcs(i,k,j) + tprs_rcs1(k)
1660 tprr_rci(i,k,j) = tprr_rci(i,k,j) + tprr_rci1(k)
1661 tprg_rcg(i,k,j) = tprg_rcg(i,k,j) + tprg_rcg1(k)
1662 tprw_vcd_c(i,k,j) = tprw_vcd_c(i,k,j) + tprw_vcd1_c(k)
1663 tprw_vcd_e(i,k,j) = tprw_vcd_e(i,k,j) + tprw_vcd1_e(k)
1664 tprr_sml(i,k,j) = tprr_sml(i,k,j) + tprr_sml1(k)
1665 tprr_gml(i,k,j) = tprr_gml(i,k,j) + tprr_gml1(k)
1666 tprr_rcg(i,k,j) = tprr_rcg(i,k,j) + tprr_rcg1(k)
1667 tprr_rcs(i,k,j) = tprr_rcs(i,k,j) + tprr_rcs1(k)
1668 tprv_rev(i,k,j) = tprv_rev(i,k,j) + tprv_rev1(k)
1669 tten3(i,k,j) = tten3(i,k,j) + tten1(k)
1670 qvten3(i,k,j) = qvten3(i,k,j) + qvten1(k)
1671 qrten3(i,k,j) = qrten3(i,k,j) + qrten1(k)
1672 qsten3(i,k,j) = qsten3(i,k,j) + qsten1(k)
1673 qgten3(i,k,j) = qgten3(i,k,j) + qgten1(k)
1674 qiten3(i,k,j) = qiten3(i,k,j) + qiten1(k)
1675 niten3(i,k,j) = niten3(i,k,j) + niten1(k)
1676 nrten3(i,k,j) = nrten3(i,k,j) + nrten1(k)
1677 ncten3(i,k,j) = ncten3(i,k,j) + ncten1(k)
1678 qcten3(i,k,j) = qcten3(i,k,j) + qcten1(k)
1680 endif assign_extended_diagnostics
1682 if (ndt>1 .and. it==ndt)
then
1683 sr(i,j) = (pcp_sn(i,j) + pcp_gr(i,j) + pcp_ic(i,j)) / (rainnc(i,j)+r1)
1684 rainncv(i,j) = rainnc(i,j)
1685 IF (
PRESENT (snowncv) )
THEN
1686 snowncv(i,j) = snownc(i,j)
1688 IF (
PRESENT (icencv) )
THEN
1689 icencv(i,j) = icenc(i,j)
1691 IF (
PRESENT (graupelncv) )
THEN
1692 graupelncv(i,j) = graupelnc(i,j)
1698 last_step_only:
IF ((ndt>1 .and. it==ndt) .or. &
1699 (nsteps>1 .and. istep==nsteps) .or. &
1700 (nsteps==1 .and. ndt==1))
THEN
1706 diagflag_present:
IF (
PRESENT (diagflag) )
THEN
1707 if (diagflag .and. do_radar_ref == 1)
then
1710 if (fullradar_diag)
then
1716 if (
present(vt_dbz_wt))
then
1718 t1d, p1d, dbz, rand1, kts, kte, i, j, &
1719 melti, vt_dbz_wt(i,:,j), &
1723 t1d, p1d, dbz, rand1, kts, kte, i, j, &
1727 refl_10cm(i,k,j) = max(-35., dbz(k))
1730 ENDIF diagflag_present
1732 IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0)
THEN
1734 re_qc1d(k) = re_qc_min
1735 re_qi1d(k) = re_qi_min
1736 re_qs1d(k) = re_qs_min
1739 call calc_effectrad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
1740 re_qc1d, re_qi1d, re_qs1d, lsml, kts, kte)
1742 re_cloud(i,k,j) = max(re_qc_min, min(re_qc1d(k), re_qc_max))
1743 re_ice(i,k,j) = max(re_qi_min, min(re_qi1d(k), re_qi_max))
1744 re_snow(i,k,j) = max(re_qs_min, min(re_qs1d(k), re_qs_max))
1747 ENDIF last_step_only
1763 do j = j_start, j_end
1765 do i = i_start, i_end
1766 pfils(i,k,j) = pfils(i,k,j)/dt_in
1767 pflls(i,k,j) = pflls(i,k,j)/dt_in
1776 deallocate_extended_diagnostics:
if (ext_diag)
then
1777 deallocate (prw_vcdc1)
1778 deallocate (prw_vcde1)
1779 deallocate (tpri_inu1)
1780 deallocate (tpri_ide1_d)
1781 deallocate (tpri_ide1_s)
1782 deallocate (tprs_ide1)
1783 deallocate (tprs_sde1_d)
1784 deallocate (tprs_sde1_s)
1785 deallocate (tprg_gde1_d)
1786 deallocate (tprg_gde1_s)
1787 deallocate (tpri_iha1)
1788 deallocate (tpri_wfz1)
1789 deallocate (tpri_rfz1)
1790 deallocate (tprg_rfz1)
1791 deallocate (tprs_scw1)
1792 deallocate (tprg_scw1)
1793 deallocate (tprg_rcs1)
1794 deallocate (tprs_rcs1)
1795 deallocate (tprr_rci1)
1796 deallocate (tprg_rcg1)
1797 deallocate (tprw_vcd1_c)
1798 deallocate (tprw_vcd1_e)
1799 deallocate (tprr_sml1)
1800 deallocate (tprr_gml1)
1801 deallocate (tprr_rcg1)
1802 deallocate (tprr_rcs1)
1803 deallocate (tprv_rev1)
1814 end if deallocate_extended_diagnostics
1881 nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dzq, &
1882 lsml, pptrain, pptsnow, pptgraul, pptice, &
1883#if ( WRF_CHEM == 1 )
1884 rainprod, evapprod, &
1886 rand1, rand2, rand3, &
1887 kts, kte, dt, ii, jj, &
1893 prw_vcdc1, prw_vcde1, &
1894 tpri_inu1, tpri_ide1_d, tpri_ide1_s, tprs_ide1, &
1895 tprs_sde1_d, tprs_sde1_s, &
1896 tprg_gde1_d, tprg_gde1_s, tpri_iha1, tpri_wfz1, &
1897 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1, &
1898 tprg_rcs1, tprs_rcs1, tprr_rci1, &
1899 tprg_rcg1, tprw_vcd1_c, &
1900 tprw_vcd1_e, tprr_sml1, tprr_gml1, tprr_rcg1, &
1901 tprr_rcs1, tprv_rev1, &
1902 tten1, qvten1, qrten1, qsten1, &
1903 qgten1, qiten1, niten1, nrten1, ncten1, qcten1, &
1911 integer,
intent(in):: kts, kte, ii, jj
1912 real(wp),
dimension(kts:kte),
intent(inout) :: &
1913 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
1914 nr1d, nc1d, nwfa1d, nifa1d, t1d
1915 real(wp),
dimension(kts:kte),
intent(out) :: pfil1, pfll1
1916 real(wp),
dimension(kts:kte),
intent(in) :: p1d, w1d, dzq
1917 real(wp),
intent(inout) :: pptrain, pptsnow, pptgraul, pptice
1918 real(wp),
intent(in) :: dt
1919 integer,
intent(in) :: lsml
1920 real(wp),
intent(in) :: rand1, rand2, rand3
1922 logical,
intent(in) :: ext_diag
1923 logical,
intent(in) :: sedi_semi
1924 integer,
intent(in) :: decfl
1925 real(wp),
dimension(:),
intent(out),
optional :: &
1928 prw_vcde1, tpri_inu1, tpri_ide1_d, &
1929 tpri_ide1_s, tprs_ide1, &
1930 tprs_sde1_d, tprs_sde1_s, tprg_gde1_d, &
1931 tprg_gde1_s, tpri_iha1, tpri_wfz1, &
1932 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1,&
1933 tprg_rcs1, tprs_rcs1, &
1934 tprr_rci1, tprg_rcg1, &
1935 tprw_vcd1_c, tprw_vcd1_e, tprr_sml1, &
1936 tprr_gml1, tprr_rcg1, &
1937 tprr_rcs1, tprv_rev1, tten1, qvten1, &
1938 qrten1, qsten1, qgten1, qiten1, niten1, &
1939 nrten1, ncten1, qcten1
1941#if ( WRF_CHEM == 1 )
1942 real(wp),
dimension(kts:kte),
intent(inout) :: &
1947 real(wp),
dimension(kts:kte) :: tten, qvten, qcten, qiten, &
1948 qrten, qsten, qgten, niten, nrten, ncten, nwfaten, nifaten
1950 real(dp),
dimension(kts:kte) :: prw_vcd
1952 real(dp),
dimension(kts:kte) :: pnc_wcd, pnc_wau, pnc_rcw, &
1955 real(dp),
dimension(kts:kte) :: pna_rca, pna_sca, pna_gca, &
1956 pnd_rcd, pnd_scd, pnd_gcd
1958 real(dp),
dimension(kts:kte) :: prr_wau, prr_rcw, prr_rcs, &
1959 prr_rcg, prr_sml, prr_gml, &
1961 pnr_wau, pnr_rcs, pnr_rcg, &
1962 pnr_rci, pnr_sml, pnr_gml, &
1963 pnr_rev, pnr_rcr, pnr_rfz
1965 real(dp),
dimension(kts:kte) :: pri_inu, pni_inu, pri_ihm, &
1966 pni_ihm, pri_wfz, pni_wfz, &
1967 pri_rfz, pni_rfz, pri_ide, &
1968 pni_ide, pri_rci, pni_rci, &
1969 pni_sci, pni_iau, pri_iha, pni_iha
1971 real(dp),
dimension(kts:kte) :: prs_iau, prs_sci, prs_rcs, &
1972 prs_scw, prs_sde, prs_ihm, &
1975 real(dp),
dimension(kts:kte) :: prg_scw, prg_rfz, prg_gde, &
1976 prg_gcw, prg_rci, prg_rcs, &
1979 real(dp),
parameter:: zeroD0 = 0.0
1980 real(wp) :: dtcfl, rainsfc, graulsfc
1983 real(wp),
dimension(kts:kte) :: temp, pres, qv, pfll, pfil, pdummy
1984 real(wp),
dimension(kts:kte) :: rc, ri, rr, rs, rg, ni, nr, nc, nwfa, nifa
1985 real(wp),
dimension(kts:kte) :: rr_tmp, nr_tmp, rg_tmp
1986 real(wp),
dimension(kts:kte) :: rho, rhof, rhof2
1987 real(wp),
dimension(kts:kte) :: qvs, qvsi, delQvs
1988 real(wp),
dimension(kts:kte) :: satw, sati, ssatw, ssati
1989 real(wp),
dimension(kts:kte) :: diffu, visco, vsc2, &
1990 tcond, lvap, ocp, lvt2
1992 real(dp),
dimension(kts:kte) :: ilamr, ilamg, N0_r, N0_g
1993 real(wp),
dimension(kts:kte) :: mvd_r, mvd_c
1994 real(wp),
dimension(kts:kte) :: smob, smo2, smo1, smo0, &
1995 smoc, smod, smoe, smof
1997 real(wp),
dimension(kts:kte) :: sed_r, sed_s, sed_g, sed_i, sed_n,sed_c
1999 real(wp) :: rgvm, delta_tp, orho, lfus2, orhodt
2000 real(wp),
dimension(5):: onstep
2001 real(dp) :: N0_exp, N0_min, lam_exp, lamc, lamr, lamg
2002 real(dp) :: lami, ilami, ilamc
2003 real(wp) :: xDc, Dc_b, Dc_g, xDi, xDr, xDs, xDg, Ds_m, Dg_m
2004 real(dp) :: Dr_star, Dc_star
2005 real(wp) :: zeta1, zeta, taud, tau
2006 real(wp) :: stoke_r, stoke_s, stoke_g, stoke_i
2007 real(wp) :: vti, vtr, vts, vtg, vtc
2008 real(wp),
dimension(kts:kte+1):: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk, &
2010 real(wp),
dimension(kts:kte):: vts_boost
2011 real(wp) :: Mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, C_snow
2012 real(wp) :: a_, b_, loga_, A1, A2, tf
2013 real(wp) :: tempc, tc0, r_mvd1, r_mvd2, xkrat
2014 real(wp) :: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr
2015 real(wp) :: xsat, rate_max, sump, ratio
2016 real(wp) :: clap, fcd, dfcd
2017 real(wp) :: otemp, rvs, rvs_p, rvs_pp, gamsc, alphsc, t1_evap, t1_subl
2018 real(wp) :: r_frac, g_frac
2019 real(wp) :: Ef_rw, Ef_sw, Ef_gw, Ef_rr
2020 real(wp) :: Ef_ra, Ef_sa, Ef_ga
2021 real(wp) :: dtsave, odts, odt, odzq, hgt_agl, SR
2022 real(wp) :: xslw1, ygra1, zans1, eva_factor
2023 integer :: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq
2024 integer,
dimension(5) :: ksed1
2025 integer :: nir, nis, nig, nii, nic, niin
2026 integer :: idx_tc, idx_t, idx_s, idx_g1, idx_g, idx_r1, idx_r, &
2027 idx_i1, idx_i, idx_c, idx, idx_d, idx_n, idx_in
2030 logical,
dimension(kts:kte) :: L_qc, L_qi, L_qr, L_qs, L_qg
2031 logical :: debug_flag
2036 debug_flag = .false.
2039 write(*, *)
'DEBUG INFO, mp_thompson at (i,j) ', ii,
', ', jj
2151#if ( WRF_CHEM == 1 )
2221 qv(k) = max(1.e-10, qv1d(k))
2223 rho(k) = roverrv*pres(k) / (r*temp(k)*(qv(k)+roverrv))
2224 nwfa(k) = max(11.1e6*rho(k), min(9999.e6*rho(k), nwfa1d(k)*rho(k)))
2225 nifa(k) = max(nain1*0.01*rho(k), min(9999.e6*rho(k), nifa1d(k)*rho(k)))
2229 if (qc1d(k) .gt. r1)
then
2231 rc(k) = qc1d(k)*rho(k)
2232 nc(k) = max(2., min(nc1d(k)*rho(k), nt_c_max))
2234 if (nc(k).gt.10000.e6)
then
2236 elseif (nc(k).lt.100.)
then
2239 nu_c = nint(1000.e6/nc(k)) + 2
2240 nu_c = max(2, min(nu_c+nint(rand2), 15))
2242 lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
2243 xdc = (bm_r + nu_c + 1.) / lamc
2244 if (xdc.lt. d0c)
then
2245 lamc = cce(2,nu_c)/d0c
2246 elseif (xdc.gt. d0r*2.)
then
2247 lamc = cce(2,nu_c)/(d0r*2.)
2249 nc(k) = min(real(nt_c_max, kind=dp), ccg(1,nu_c)*ocg2(nu_c)*rc(k) &
2251 if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware))
then
2266 if (qi1d(k) .gt. r1)
then
2268 ri(k) = qi1d(k)*rho(k)
2269 ni(k) = max(r2, ni1d(k)*rho(k))
2270 if (ni(k).le. r2)
then
2272 ni(k) = min(nt_i_max, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
2275 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
2277 xdi = (bm_i + mu_i + 1.) * ilami
2278 if (xdi.lt. 5.e-6)
then
2280 ni(k) = min(nt_i_max, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
2281 elseif (xdi.gt. 300.e-6)
then
2282 lami = cie(2)/300.e-6
2283 ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i
2293 if (qr1d(k) .gt. r1)
then
2295 rr(k) = qr1d(k)*rho(k)
2296 nr(k) = max(r2, nr1d(k)*rho(k))
2297 if (nr(k).le. r2)
then
2299 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
2300 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
2303 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
2304 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
2305 if (mvd_r(k) .gt. 2.5e-3)
then
2307 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
2308 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
2309 elseif (mvd_r(k) .lt. d0r*0.75)
then
2311 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
2312 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
2321 if (qs1d(k) .gt. r1)
then
2323 rs(k) = qs1d(k)*rho(k)
2330 if (qg1d(k) .gt. r1)
then
2332 rg(k) = qg1d(k)*rho(k)
2357 tempc = temp(k) - 273.15
2358 rhof(k) = sqrt(rho_not/rho(k))
2359 rhof2(k) = sqrt(rhof(k))
2360 qvs(k) =
rslf(pres(k), temp(k))
2361 delqvs(k) = max(0.0,
rslf(pres(k), 273.15)-qv(k))
2362 if (tempc .le. 0.0)
then
2363 qvsi(k) =
rsif(pres(k), temp(k))
2367 satw(k) = qv(k)/qvs(k)
2368 sati(k) = qv(k)/qvsi(k)
2369 ssatw(k) = satw(k) - 1.
2370 ssati(k) = sati(k) - 1.
2371 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
2372 if (abs(ssati(k)).lt. eps) ssati(k) = 0.0
2373 if (no_micro .and. ssati(k).gt. 0.0) no_micro = .false.
2374 diffu(k) = 2.11e-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
2375 if (tempc .ge. 0.0)
then
2376 visco(k) = (1.718+0.0049*tempc)*1.0e-5
2378 visco(k) = (1.718+0.0049*tempc-1.2e-5*tempc*tempc)*1.0e-5
2380 ocp(k) = 1./(cp*(1.+0.887*qv(k)))
2381 vsc2(k) = sqrt(rho(k)/visco(k))
2382 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
2383 tcond(k) = (5.69 + 0.0168*tempc)*1.0e-5 * 418.936
2391 if (no_micro)
return
2396 if (.not. iiwarm)
then
2398 if (.not. l_qs(k)) cycle
2399 tc0 = min(-0.1, temp(k)-273.15)
2400 smob(k) = rs(k)*oams
2404 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3))
then
2407 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
2408 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
2409 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
2410 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
2411 + sa(10)*bm_s*bm_s*bm_s
2413 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
2414 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
2415 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
2416 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
2417 + sb(10)*bm_s*bm_s*bm_s
2418 smo2(k) = (smob(k)/a_)**(1./b_)
2422 loga_ = sa(1) + sa(2)*tc0 + sa(5)*tc0*tc0 + sa(9)*tc0*tc0*tc0
2424 b_ = sb(1) + sb(2)*tc0 + sb(5)*tc0*tc0 + sb(9)*tc0*tc0*tc0
2425 smo0(k) = a_ * smo2(k)**b_
2428 loga_ = sa(1) + sa(2)*tc0 + sa(3) &
2429 + sa(4)*tc0 + sa(5)*tc0*tc0 &
2430 + sa(6) + sa(7)*tc0*tc0 &
2431 + sa(8)*tc0 + sa(9)*tc0*tc0*tc0 &
2434 b_ = sb(1)+ sb(2)*tc0 + sb(3) + sb(4)*tc0 &
2435 + sb(5)*tc0*tc0 + sb(6) &
2436 + sb(7)*tc0*tc0 + sb(8)*tc0 &
2437 + sb(9)*tc0*tc0*tc0 + sb(10)
2438 smo1(k) = a_ * smo2(k)**b_
2441 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
2442 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
2443 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
2444 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
2445 + sa(10)*cse(1)*cse(1)*cse(1)
2447 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
2448 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
2449 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
2450 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
2451 smoc(k) = a_ * smo2(k)**b_
2454 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(13) &
2455 + sa(4)*tc0*cse(13) + sa(5)*tc0*tc0 &
2456 + sa(6)*cse(13)*cse(13) + sa(7)*tc0*tc0*cse(13) &
2457 + sa(8)*tc0*cse(13)*cse(13) + sa(9)*tc0*tc0*tc0 &
2458 + sa(10)*cse(13)*cse(13)*cse(13)
2460 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(13) + sb(4)*tc0*cse(13) &
2461 + sb(5)*tc0*tc0 + sb(6)*cse(13)*cse(13) &
2462 + sb(7)*tc0*tc0*cse(13) + sb(8)*tc0*cse(13)*cse(13) &
2463 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(13)*cse(13)*cse(13)
2464 smoe(k) = a_ * smo2(k)**b_
2467 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(16) &
2468 + sa(4)*tc0*cse(16) + sa(5)*tc0*tc0 &
2469 + sa(6)*cse(16)*cse(16) + sa(7)*tc0*tc0*cse(16) &
2470 + sa(8)*tc0*cse(16)*cse(16) + sa(9)*tc0*tc0*tc0 &
2471 + sa(10)*cse(16)*cse(16)*cse(16)
2473 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(16) + sb(4)*tc0*cse(16) &
2474 + sb(5)*tc0*tc0 + sb(6)*cse(16)*cse(16) &
2475 + sb(7)*tc0*tc0*cse(16) + sb(8)*tc0*cse(16)*cse(16) &
2476 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(16)*cse(16)*cse(16)
2477 smof(k) = a_ * smo2(k)**b_
2490 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
2492 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
2493 n0_r(k) = nr(k)*org2*lamr**cre(2)
2504 if (l_qr(k) .and. mvd_r(k).gt. d0r)
then
2505 ef_rr = max(-0.1, 1.0 - exp(2300.0*(mvd_r(k)-1950.0e-6)))
2506 pnr_rcr(k) = ef_rr * 2.0*nr(k)*rr(k)
2510 if (nc(k).gt.10000.e6)
then
2512 elseif (nc(k).lt.100.)
then
2515 nu_c = nint(1000.e6/nc(k)) + 2
2516 nu_c = max(2, min(nu_c+nint(rand2), 15))
2518 xdc = max(d0c*1.e6, ((rc(k)/(am_r*nc(k)))**obmr) * 1.e6)
2519 lamc = (nc(k)*am_r* ccg(2,nu_c) * ocg1(nu_c) / rc(k))**obmr
2520 mvd_c(k) = (3.0+nu_c+0.672) / lamc
2521 mvd_c(k) = max(d0c, min(mvd_c(k), d0r))
2526 if (rc(k).gt. 0.01e-3)
then
2527 dc_g = ((ccg(3,nu_c)*ocg2(nu_c))**obmr / lamc) * 1.e6
2528 dc_b = (xdc*xdc*xdc*dc_g*dc_g*dc_g - xdc*xdc*xdc*xdc*xdc*xdc) &
2530 zeta1 = 0.5*((6.25e-6*xdc*dc_b*dc_b*dc_b - 0.4) &
2531 + abs(6.25e-6*xdc*dc_b*dc_b*dc_b - 0.4))
2532 zeta = 0.027*rc(k)*zeta1
2533 taud = 0.5*((0.5*dc_b - 7.5) + abs(0.5*dc_b - 7.5)) + r1
2534 tau = 3.72/(rc(k)*taud)
2535 prr_wau(k) = zeta/tau
2536 prr_wau(k) = min(real(rc(k)*odts, kind=dp), prr_wau(k))
2537 pnr_wau(k) = prr_wau(k) / (am_r*nu_c*10.*d0r*d0r*d0r)
2538 pnc_wau(k) = min(real(nc(k)*odts, kind=dp), prr_wau(k) &
2539 / (am_r*mvd_c(k)*mvd_c(k)*mvd_c(k)))
2543 if (l_qr(k) .and. mvd_r(k).gt. d0r .and. mvd_c(k).gt. d0c)
then
2545 idx = 1 + int(nbr*log(real(mvd_r(k)/dr(1), kind=dp)) / log(dr(nbr)/dr(1)))
2547 ef_rw = t_efrw(idx, int(mvd_c(k)*1.e6))
2548 prr_rcw(k) = rhof(k)*t1_qr_qc*ef_rw*rc(k)*n0_r(k) &
2549 *((lamr+fv_r)**(-cre(9)))
2550 prr_rcw(k) = min(real(rc(k)*odts, kind=dp), prr_rcw(k))
2551 pnc_rcw(k) = rhof(k)*t1_qr_qc*ef_rw*nc(k)*n0_r(k) &
2552 *((lamr+fv_r)**(-cre(9)))
2553 pnc_rcw(k) = min(real(nc(k)*odts, kind=dp), pnc_rcw(k))
2557 if (l_qr(k) .and. mvd_r(k).gt. d0r)
then
2558 ef_ra =
eff_aero(mvd_r(k),0.04e-6,visco(k),rho(k),temp(k),
'r')
2560 pna_rca(k) = rhof(k)*t1_qr_qc*ef_ra*nwfa(k)*n0_r(k) &
2561 *((lamr+fv_r)**(-cre(9)))
2562 pna_rca(k) = min(real(nwfa(k)*odts, kind=dp), pna_rca(k))
2564 ef_ra =
eff_aero(mvd_r(k),0.8e-6,visco(k),rho(k),temp(k),
'r')
2565 pnd_rcd(k) = rhof(k)*t1_qr_qc*ef_ra*nifa(k)*n0_r(k) &
2566 *((lamr+fv_r)**(-cre(9)))
2567 pnd_rcd(k) = min(real(nifa(k)*odts, kind=dp), pnd_rcd(k))
2575 if (.not. iiwarm)
then
2579 if (l_qs(k)) xds = smoc(k) / smob(k)
2582 tempc = temp(k) - 273.15
2583 idx_tc = max(1, min(nint(-tempc), 45) )
2584 idx_t = int( (tempc-2.5)/5. ) - 1
2585 idx_t = max(1, -idx_t)
2586 idx_t = min(idx_t, ntb_t)
2587 it = max(1, min(nint(-tempc), 31) )
2590 if (rc(k).gt. r_c(1))
then
2591 nic = nint(log10(rc(k)))
2592 do_loop_rc:
do nn = nic-1, nic+1
2594 if ( (rc(k)/10.**nn).ge.1.0 .and. (rc(k)/10.**nn).lt.10.0 )
exit do_loop_rc
2596 idx_c = int(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
2597 idx_c = max(1, min(idx_c, ntb_c))
2603 idx_n = nint(1.0 + real(nbc, kind=wp) * log(real(nc(k)/t_nc(1), kind=dp)) / nic1)
2604 idx_n = max(1, min(idx_n, nbc))
2607 if (ri(k).gt. r_i(1))
then
2608 nii = nint(log10(ri(k)))
2609 do_loop_ri:
do nn = nii-1, nii+1
2611 if ( (ri(k)/10.**nn).ge.1.0 .and. (ri(k)/10.**nn).lt.10.0 )
exit do_loop_ri
2613 idx_i = int(ri(k)/10.**n) + 10*(n-nii2) - (n-nii2)
2614 idx_i = max(1, min(idx_i, ntb_i))
2619 if (ni(k).gt. nt_i(1))
then
2620 nii = nint(log10(ni(k)))
2621 do_loop_ni:
do nn = nii-1, nii+1
2623 if ( (ni(k)/10.**nn).ge.1.0 .and. (ni(k)/10.**nn).lt.10.0 )
exit do_loop_ni
2625 idx_i1 = int(ni(k)/10.**n) + 10*(n-nii3) - (n-nii3)
2626 idx_i1 = max(1, min(idx_i1, ntb_i1))
2632 if (rr(k).gt. r_r(1))
then
2633 nir = nint(log10(rr(k)))
2634 do_loop_rr:
do nn = nir-1, nir+1
2636 if ( (rr(k)/10.**nn).ge.1.0 .and. (rr(k)/10.**nn).lt.10.0 )
exit do_loop_rr
2638 idx_r = int(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2)
2639 idx_r = max(1, min(idx_r, ntb_r))
2642 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
2643 n0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
2644 nir = nint(log10(n0_exp))
2645 do_loop_nr:
do nn = nir-1, nir+1
2647 if ( (n0_exp/10.**nn).ge.1.0 .and. (n0_exp/10.**nn).lt.10.0 )
exit do_loop_nr
2649 idx_r1 = int(n0_exp/10.**n) + 10*(n-nir3) - (n-nir3)
2650 idx_r1 = max(1, min(idx_r1, ntb_r1))
2657 if (rs(k).gt. r_s(1))
then
2658 nis = nint(log10(rs(k)))
2659 do_loop_rs:
do nn = nis-1, nis+1
2661 if ( (rs(k)/10.**nn).ge.1.0 .and. (rs(k)/10.**nn).lt.10.0 )
exit do_loop_rs
2663 idx_s = int(rs(k)/10.**n) + 10*(n-nis2) - (n-nis2)
2664 idx_s = max(1, min(idx_s, ntb_s))
2670 if (rg(k).gt. r_g(1))
then
2671 nig = nint(log10(rg(k)))
2672 do_loop_rg:
do nn = nig-1, nig+1
2674 if ( (rg(k)/10.**nn).ge.1.0 .and. (rg(k)/10.**nn).lt.10.0 )
exit do_loop_rg
2676 idx_g = int(rg(k)/10.**n) + 10*(n-nig2) - (n-nig2)
2677 idx_g = max(1, min(idx_g, ntb_g))
2680 lam_exp = lamg * (cgg(3)*ogg2*ogg1)**bm_g
2681 n0_exp = ogg1*rg(k)/am_g * lam_exp**cge(1)
2682 nig = nint(log10(real(n0_exp, kind=dp)))
2683 do_loop_ng:
do nn = nig-1, nig+1
2685 if ( (n0_exp/10.**nn).ge.1.0 .and. (n0_exp/10.**nn).lt.10.0 )
exit do_loop_ng
2687 idx_g1 = int(n0_exp/10.**n) + 10*(n-nig3) - (n-nig3)
2688 idx_g1 = max(1, min(idx_g1, ntb_g1))
2696 rvs = rho(k)*qvsi(k)
2697 rvs_p = rvs*otemp*(lsub*otemp*orv - 1.)
2698 rvs_pp = rvs * ( otemp*(lsub*otemp*orv - 1.) &
2699 *otemp*(lsub*otemp*orv - 1.) &
2700 + (-2.*lsub*otemp*otemp*otemp*orv) &
2702 gamsc = lsub*diffu(k)/tcond(k) * rvs_p
2703 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
2704 * rvs_pp/rvs_p * rvs/rvs_p
2705 alphsc = max(1.e-9, alphsc)
2707 if (abs(xsat).lt. 1.e-9) xsat=0.
2708 t1_subl = 4.*pi*( 1.0 - alphsc*xsat &
2709 + 2.*alphsc*alphsc*xsat*xsat &
2710 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
2714 if (l_qc(k) .and. mvd_c(k).gt. d0c)
then
2715 if (xds .gt. d0s)
then
2716 idx = 1 + int(nbs*log(real(xds/ds(1), kind=dp)) / log(ds(nbs)/ds(1)))
2718 ef_sw = t_efsw(idx, int(mvd_c(k)*1.e6))
2719 prs_scw(k) = rhof(k)*t1_qs_qc*ef_sw*rc(k)*smoe(k)
2720 prs_scw(k) = min(real(rc(k)*odts, kind=dp), prs_scw(k))
2721 pnc_scw(k) = rhof(k)*t1_qs_qc*ef_sw*nc(k)*smoe(k)
2722 pnc_scw(k) = min(real(nc(k)*odts, kind=dp), pnc_scw(k))
2726 if (rg(k).ge. r_g(1) .and. mvd_c(k).gt. d0c)
then
2727 xdg = (bm_g + mu_g + 1.) * ilamg(k)
2728 vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
2729 stoke_g = mvd_c(k)*mvd_c(k)*vtg*rho_w/(9.*visco(k)*xdg)
2730 if (xdg.gt. d0g)
then
2731 if (stoke_g.ge.0.4 .and. stoke_g.le.10.)
then
2732 ef_gw = 0.55*log10(2.51*stoke_g)
2733 elseif (stoke_g.lt.0.4)
then
2735 elseif (stoke_g.gt.10)
then
2738 prg_gcw(k) = rhof(k)*t1_qg_qc*ef_gw*rc(k)*n0_g(k) &
2740 pnc_gcw(k) = rhof(k)*t1_qg_qc*ef_gw*nc(k)*n0_g(k) &
2742 pnc_gcw(k) = min(real(nc(k)*odts, kind=dp), pnc_gcw(k))
2748 if (rs(k) .gt. r_s(1))
then
2749 ef_sa =
eff_aero(xds,0.04e-6,visco(k),rho(k),temp(k),
's')
2750 pna_sca(k) = rhof(k)*t1_qs_qc*ef_sa*nwfa(k)*smoe(k)
2751 pna_sca(k) = min(real(nwfa(k)*odts, kind=dp), pna_sca(k))
2753 ef_sa =
eff_aero(xds,0.8e-6,visco(k),rho(k),temp(k),
's')
2754 pnd_scd(k) = rhof(k)*t1_qs_qc*ef_sa*nifa(k)*smoe(k)
2755 pnd_scd(k) = min(real(nifa(k)*odts, kind=dp), pnd_scd(k))
2757 if (rg(k) .gt. r_g(1))
then
2758 xdg = (bm_g + mu_g + 1.) * ilamg(k)
2759 ef_ga =
eff_aero(xdg,0.04e-6,visco(k),rho(k),temp(k),
'g')
2760 pna_gca(k) = rhof(k)*t1_qg_qc*ef_ga*nwfa(k)*n0_g(k) &
2762 pna_gca(k) = min(real(nwfa(k)*odts, kind=dp), pna_gca(k))
2764 ef_ga =
eff_aero(xdg,0.8e-6,visco(k),rho(k),temp(k),
'g')
2765 pnd_gcd(k) = rhof(k)*t1_qg_qc*ef_ga*nifa(k)*n0_g(k) &
2767 pnd_gcd(k) = min(real(nifa(k)*odts, kind=dp), pnd_gcd(k))
2773 if (rr(k).ge. r_r(1))
then
2774 if (rs(k).ge. r_s(1))
then
2775 if (temp(k).lt.t_0)
then
2776 prr_rcs(k) = -(tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2777 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
2778 + tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
2779 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r))
2780 prs_rcs(k) = tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2781 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
2782 - tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2783 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
2784 prg_rcs(k) = tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
2785 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2786 + tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2787 + tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
2788 prr_rcs(k) = max(real(-rr(k)*odts, kind=dp), prr_rcs(k))
2789 prs_rcs(k) = max(real(-rs(k)*odts, kind=dp), prs_rcs(k))
2790 prg_rcs(k) = min(real((rr(k)+rs(k))*odts, kind=dp), prg_rcs(k))
2791 pnr_rcs(k) = tnr_racs1(idx_s,idx_t,idx_r1,idx_r) &
2792 + tnr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2793 + tnr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2794 + tnr_sacr2(idx_s,idx_t,idx_r1,idx_r)
2795 pnr_rcs(k) = min(real(nr(k)*odts, kind=dp), pnr_rcs(k))
2797 prs_rcs(k) = -tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2798 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2799 + tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2800 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r)
2801 prs_rcs(k) = max(real(-rs(k)*odts, kind=dp), prs_rcs(k))
2802 prr_rcs(k) = -prs_rcs(k)
2809 if (rg(k).ge. r_g(1))
then
2810 if (temp(k).lt.t_0)
then
2811 prg_rcg(k) = tmr_racg(idx_g1,idx_g,idx_r1,idx_r) &
2812 + tcr_gacr(idx_g1,idx_g,idx_r1,idx_r)
2813 prg_rcg(k) = min(real(rr(k)*odts, kind=dp), prg_rcg(k))
2814 prr_rcg(k) = -prg_rcg(k)
2815 pnr_rcg(k) = tnr_racg(idx_g1,idx_g,idx_r1,idx_r) &
2816 + tnr_gacr(idx_g1,idx_g,idx_r1,idx_r)
2817 pnr_rcg(k) = min(real(nr(k)*odts, kind=dp), pnr_rcg(k))
2819 prr_rcg(k) = tcg_racg(idx_g1,idx_g,idx_r1,idx_r)
2820 prr_rcg(k) = min(real(rg(k)*odts, kind=dp), prr_rcg(k))
2821 prg_rcg(k) = -prr_rcg(k)
2823 pnr_rcg(k) = -1.5*tnr_gacr(idx_g1,idx_g,idx_r1,idx_r)
2828 if (temp(k).lt.t_0)
then
2829 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
2833 c_snow = c_sqrd + (tempc+1.5)*(c_cube-c_sqrd)/(-30.+1.5)
2834 c_snow = max(c_sqrd, min(c_snow, c_cube))
2835 prs_sde(k) = c_snow*t1_subl*diffu(k)*ssati(k)*rvs &
2836 * (t1_qs_sd*smo1(k) &
2837 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
2838 if (prs_sde(k).lt. 0.)
then
2839 prs_sde(k) = max(real(-rs(k)*odts, kind=dp), prs_sde(k), real(rate_max, kind=dp))
2841 prs_sde(k) = min(prs_sde(k), real(rate_max, kind=dp))
2845 if (l_qg(k) .and. ssati(k).lt. -eps)
then
2846 prg_gde(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
2847 * n0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
2848 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
2849 if (prg_gde(k).lt. 0.)
then
2850 prg_gde(k) = max(real(-rg(k)*odts, kind=dp), prg_gde(k), real(rate_max, kind=dp))
2852 prg_gde(k) = min(prg_gde(k), real(rate_max, kind=dp))
2860 if (prs_scw(k).gt.5.0*prs_sde(k) .and. &
2861 prs_sde(k).gt.eps)
then
2862 r_frac = min(30.0_dp, prs_scw(k)/prs_sde(k))
2863 g_frac = min(0.75, 0.15 + (r_frac-5.)*.028)
2864 vts_boost(k) = min(1.5, 1.1 + (r_frac-5.)*.016)
2865 prg_scw(k) = g_frac*prs_scw(k)
2866 prs_scw(k) = (1. - g_frac)*prs_scw(k)
2874 if (temp(k).lt.t_0)
then
2876 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
2892 if (dustyice .AND. (is_aerosol_aware .or. merra2_aerosol_aware))
then
2893 xni =
icedemott(tempc,qvs(k),qvs(k),qvsi(k),rho(k),nifa(k))
2899 if (xni.gt. nt_in(1))
then
2900 niin = nint(log10(xni))
2901 do_loop_xni:
do nn = niin-1, niin+1
2903 if ( (xni/10.**nn).ge.1.0 .and. (xni/10.**nn).lt.10.0 )
exit do_loop_xni
2905 idx_in = int(xni/10.**n) + 10*(n-niin2) - (n-niin2)
2906 idx_in = max(1, min(idx_in, ntb_in))
2912 if (rr(k).gt. r_r(1))
then
2913 prg_rfz(k) = tpg_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
2914 pri_rfz(k) = tpi_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
2915 pni_rfz(k) = tni_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
2916 pnr_rfz(k) = tnr_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
2917 pnr_rfz(k) = min(real(nr(k)*odts, kind=dp), pnr_rfz(k))
2918 elseif (rr(k).gt. r1 .and. temp(k).lt.hgfr)
then
2919 pri_rfz(k) = rr(k)*odts
2920 pni_rfz(k) = pnr_rfz(k)
2923 if (rc(k).gt. r_c(1))
then
2924 pri_wfz(k) = tpi_qcfz(idx_c,idx_n,idx_tc,idx_in)*odts
2925 pri_wfz(k) = min(real(rc(k)*odts, kind=dp), pri_wfz(k))
2926 pni_wfz(k) = tni_qcfz(idx_c,idx_n,idx_tc,idx_in)*odts
2927 pni_wfz(k) = min(real(nc(k)*odts, kind=dp), pri_wfz(k)/(2.0_dp*xm0i), &
2929 elseif (rc(k).gt. r1 .and. temp(k).lt.hgfr)
then
2930 pri_wfz(k) = rc(k)*odts
2931 pni_wfz(k) = nc(k)*odts
2936 if ( (ssati(k).ge. ssati_min) .or. (ssatw(k).gt. eps &
2937 .and. temp(k).lt.253.15) )
then
2938 if (dustyice .AND. (is_aerosol_aware .or. merra2_aerosol_aware))
then
2939 xnc =
icedemott(tempc,qv(k),qvs(k),qvsi(k),rho(k),nifa(k))
2940 xnc = xnc*(1.0 + 50.*rand3)
2942 xnc = min(xnc_max, tno*exp(ato*(t_0-temp(k))))
2944 xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dtsave
2945 pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts
2946 pri_inu(k) = min(real(rate_max, kind=dp), xm0i*pni_inu(k))
2947 pni_inu(k) = pri_inu(k)/xm0i
2951 xni = smo0(k)+ni(k) + (pni_rfz(k)+pni_wfz(k)+pni_inu(k))*dtsave
2952 if ((is_aerosol_aware .or. merra2_aerosol_aware) .AND. homogice .AND. (xni.le.nt_i_max) &
2953 .AND.(temp(k).lt.238).AND.(ssati(k).ge.0.4) )
then
2954 xnc =
icekoop(temp(k),qv(k),qvs(k),nwfa(k), dtsave)
2955 pni_iha(k) = xnc*odts
2956 pri_iha(k) = min(real(rate_max, kind=dp), xm0i*0.1*pni_iha(k))
2957 pni_iha(k) = pri_iha(k)/(xm0i*0.1)
2964 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
2966 xdi = max(real(d0i, kind=dp), (bm_i + mu_i + 1.) * ilami)
2967 xmi = am_i*xdi**bm_i
2969 pri_ide(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
2970 *oig1*cig(5)*ni(k)*ilami
2972 if (pri_ide(k) .lt. 0.0)
then
2973 pri_ide(k) = max(real(-ri(k)*odts, kind=dp), pri_ide(k), real(rate_max, kind=dp))
2974 pni_ide(k) = pri_ide(k)*oxmi
2975 pni_ide(k) = max(real(-ni(k)*odts, kind=dp), pni_ide(k))
2977 pri_ide(k) = min(pri_ide(k), real(rate_max, kind=dp))
2978 prs_ide(k) = (1.0_dp-tpi_ide(idx_i,idx_i1))*pri_ide(k)
2979 pri_ide(k) = tpi_ide(idx_i,idx_i1)*pri_ide(k)
2984 if ( (idx_i.eq. ntb_i) .or. (xdi.gt. 5.0*d0s) )
then
2985 prs_iau(k) = ri(k)*.99*odts
2986 pni_iau(k) = ni(k)*.95*odts
2987 elseif (xdi.lt. 0.1*d0s)
then
2991 prs_iau(k) = tps_iaus(idx_i,idx_i1)*odts
2992 prs_iau(k) = min(real(ri(k)*.99*odts, kind=dp), prs_iau(k))
2993 pni_iau(k) = tni_iaus(idx_i,idx_i1)*odts
2994 pni_iau(k) = min(real(ni(k)*.95*odts, kind=dp), pni_iau(k))
3000 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
3002 xdi = max(real(d0i, kind=dp), (bm_i + mu_i + 1.) * ilami)
3003 xmi = am_i*xdi**bm_i
3005 if (rs(k).ge. r_s(1))
then
3006 prs_sci(k) = t1_qs_qi*rhof(k)*ef_si*ri(k)*smoe(k)
3007 pni_sci(k) = prs_sci(k) * oxmi
3011 if (rr(k).ge. r_r(1) .and. mvd_r(k).gt. 4.*xdi)
then
3013 pri_rci(k) = rhof(k)*t1_qr_qi*ef_ri*ri(k)*n0_r(k) &
3014 *((lamr+fv_r)**(-cre(9)))
3015 pnr_rci(k) = rhof(k)*t1_qr_qi*ef_ri*ni(k)*n0_r(k) &
3016 *((lamr+fv_r)**(-cre(9)))
3017 pni_rci(k) = pri_rci(k) * oxmi
3018 prr_rci(k) = rhof(k)*t2_qr_qi*ef_ri*ni(k)*n0_r(k) &
3019 *((lamr+fv_r)**(-cre(8)))
3020 prr_rci(k) = min(real(rr(k)*odts, kind=dp), prr_rci(k))
3021 prg_rci(k) = pri_rci(k) + prr_rci(k)
3026 if (prg_gcw(k).gt. eps .and. tempc.gt.-8.0)
then
3028 if (tempc.ge.-5.0 .and. tempc.lt.-3.0)
then
3029 tf = 0.5*(-3.0 - tempc)
3030 elseif (tempc.gt.-8.0 .and. tempc.lt.-5.0)
then
3031 tf = 0.33333333*(8.0 + tempc)
3033 pni_ihm(k) = 3.5e8*tf*prg_gcw(k)
3034 pri_ihm(k) = xm0i*pni_ihm(k)
3035 prs_ihm(k) = prs_scw(k)/(prs_scw(k)+prg_gcw(k)) &
3037 prg_ihm(k) = prg_gcw(k)/(prs_scw(k)+prg_gcw(k)) &
3046 prr_sml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delqvs(k)) &
3047 * (t1_qs_me*smo1(k) + t2_qs_me*rhof2(k)*vsc2(k)*smof(k))
3048 if (prr_sml(k) .gt. 0.)
then
3049 prr_sml(k) = prr_sml(k) + 4218.*olfus*tempc &
3050 * (prr_rcs(k)+prs_scw(k))
3051 prr_sml(k) = min(real(rs(k)*odts, kind=dp), prr_sml(k))
3052 pnr_sml(k) = smo0(k)/rs(k)*prr_sml(k) * 10.0**(-0.25*tempc)
3053 pnr_sml(k) = min(real(smo0(k)*odts, kind=dp), pnr_sml(k))
3054 elseif (ssati(k).lt. 0.)
then
3056 prs_sde(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
3057 * (t1_qs_sd*smo1(k) &
3058 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
3059 prs_sde(k) = max(real(-rs(k)*odts, kind=dp), prs_sde(k))
3064 prr_gml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delqvs(k)) &
3065 * n0_g(k)*(t1_qg_me*ilamg(k)**cge(10) &
3066 + t2_qg_me*rhof2(k)*vsc2(k)*ilamg(k)**cge(11))
3067 if (prr_gml(k) .gt. 0.)
then
3068 prr_gml(k) = min(real(rg(k)*odts, kind=dp), prr_gml(k))
3069 pnr_gml(k) = n0_g(k)*cgg(2)*ilamg(k)**cge(2) / rg(k) &
3070 * prr_gml(k) * 10.0**(-0.5*tempc)
3071 elseif (ssati(k).lt. 0.)
then
3073 prg_gde(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
3074 * n0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
3075 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
3076 prg_gde(k) = max(real(-rg(k)*odts, kind=dp), prg_gde(k))
3084 if (dt .gt. 120.)
then
3085 prr_rcw(k)=prr_rcw(k)+prs_scw(k)+prg_gcw(k)
3103 sump = pri_inu(k) + pri_ide(k) + prs_ide(k) &
3104 + prs_sde(k) + prg_gde(k) + pri_iha(k)
3105 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
3106 if ( (sump.gt. eps .and. sump.gt. rate_max) .or. &
3107 (sump.lt. -eps .and. sump.lt. rate_max) )
then
3108 ratio = rate_max/sump
3109 pri_inu(k) = pri_inu(k) * ratio
3110 pri_ide(k) = pri_ide(k) * ratio
3111 pni_ide(k) = pni_ide(k) * ratio
3112 prs_ide(k) = prs_ide(k) * ratio
3113 prs_sde(k) = prs_sde(k) * ratio
3114 prg_gde(k) = prg_gde(k) * ratio
3115 pri_iha(k) = pri_iha(k) * ratio
3119 sump = -prr_wau(k) - pri_wfz(k) - prr_rcw(k) &
3120 - prs_scw(k) - prg_scw(k) - prg_gcw(k)
3121 rate_max = -rc(k)*odts
3122 if (sump.lt. rate_max .and. l_qc(k))
then
3123 ratio = rate_max/sump
3124 prr_wau(k) = prr_wau(k) * ratio
3125 pri_wfz(k) = pri_wfz(k) * ratio
3126 prr_rcw(k) = prr_rcw(k) * ratio
3127 prs_scw(k) = prs_scw(k) * ratio
3128 prg_scw(k) = prg_scw(k) * ratio
3129 prg_gcw(k) = prg_gcw(k) * ratio
3133 sump = pri_ide(k) - prs_iau(k) - prs_sci(k) &
3135 rate_max = -ri(k)*odts
3136 if (sump.lt. rate_max .and. l_qi(k))
then
3137 ratio = rate_max/sump
3138 pri_ide(k) = pri_ide(k) * ratio
3139 prs_iau(k) = prs_iau(k) * ratio
3140 prs_sci(k) = prs_sci(k) * ratio
3141 pri_rci(k) = pri_rci(k) * ratio
3145 sump = -prg_rfz(k) - pri_rfz(k) - prr_rci(k) &
3146 + prr_rcs(k) + prr_rcg(k)
3147 rate_max = -rr(k)*odts
3148 if (sump.lt. rate_max .and. l_qr(k))
then
3149 ratio = rate_max/sump
3150 prg_rfz(k) = prg_rfz(k) * ratio
3151 pri_rfz(k) = pri_rfz(k) * ratio
3152 prr_rci(k) = prr_rci(k) * ratio
3153 prr_rcs(k) = prr_rcs(k) * ratio
3154 prr_rcg(k) = prr_rcg(k) * ratio
3158 sump = prs_sde(k) - prs_ihm(k) - prr_sml(k) &
3160 rate_max = -rs(k)*odts
3161 if (sump.lt. rate_max .and. l_qs(k))
then
3162 ratio = rate_max/sump
3163 prs_sde(k) = prs_sde(k) * ratio
3164 prs_ihm(k) = prs_ihm(k) * ratio
3165 prr_sml(k) = prr_sml(k) * ratio
3166 prs_rcs(k) = prs_rcs(k) * ratio
3170 sump = prg_gde(k) - prg_ihm(k) - prr_gml(k) &
3172 rate_max = -rg(k)*odts
3173 if (sump.lt. rate_max .and. l_qg(k))
then
3174 ratio = rate_max/sump
3175 prg_gde(k) = prg_gde(k) * ratio
3176 prg_ihm(k) = prg_ihm(k) * ratio
3177 prr_gml(k) = prr_gml(k) * ratio
3178 prg_rcg(k) = prg_rcg(k) * ratio
3183 pri_ihm(k) = prs_ihm(k) + prg_ihm(k)
3184 ratio = min( abs(prr_rcg(k)), abs(prg_rcg(k)) )
3185 prr_rcg(k) = ratio * sign(1.0, sngl(prr_rcg(k)))
3186 prg_rcg(k) = -prr_rcg(k)
3187 if (temp(k).gt.t_0)
then
3188 ratio = min( abs(prr_rcs(k)), abs(prs_rcs(k)) )
3189 prr_rcs(k) = ratio * sign(1.0, sngl(prr_rcs(k)))
3190 prs_rcs(k) = -prr_rcs(k)
3201 lfus2 = lsub - lvap(k)
3204 if (is_aerosol_aware)
then
3205 nwfaten(k) = nwfaten(k) - (pna_rca(k) + pna_sca(k) &
3206 + pna_gca(k) + pni_iha(k)) * orho
3207 nifaten(k) = nifaten(k) - (pnd_rcd(k) + pnd_scd(k) &
3208 + pnd_gcd(k)) * orho
3210 nifaten(k) = nifaten(k) - pni_inu(k)*orho
3217 qvten(k) = qvten(k) + (-pri_inu(k) - pri_iha(k) - pri_ide(k) &
3218 - prs_ide(k) - prs_sde(k) - prg_gde(k)) &
3222 qcten(k) = qcten(k) + (-prr_wau(k) - pri_wfz(k) &
3223 - prr_rcw(k) - prs_scw(k) - prg_scw(k) &
3228 ncten(k) = ncten(k) + (-pnc_wau(k) - pnc_rcw(k) &
3229 - pni_wfz(k) - pnc_scw(k) - pnc_gcw(k)) &
3234 xrc=max(r1, (qc1d(k) + qcten(k)*dtsave)*rho(k))
3235 xnc=max(2., (nc1d(k) + ncten(k)*dtsave)*rho(k))
3236 if (xrc .gt. r1)
then
3237 if (xnc.gt.10000.e6)
then
3239 elseif (xnc.lt.100.)
then
3242 nu_c = nint(1000.e6/xnc) + 2
3243 nu_c = max(2, min(nu_c+nint(rand2), 15))
3245 lamc = (xnc*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
3246 xdc = (bm_r + nu_c + 1.) / lamc
3247 if (xdc.lt. d0c)
then
3248 lamc = cce(2,nu_c)/d0c
3249 xnc = ccg(1,nu_c)*ocg2(nu_c)*xrc/am_r*lamc**bm_r
3250 ncten(k) = (xnc-nc1d(k)*rho(k))*odts*orho
3251 elseif (xdc.gt. d0r*2.)
then
3252 lamc = cce(2,nu_c)/(d0r*2.)
3253 xnc = ccg(1,nu_c)*ocg2(nu_c)*xrc/am_r*lamc**bm_r
3254 ncten(k) = (xnc-nc1d(k)*rho(k))*odts*orho
3257 ncten(k) = -nc1d(k)*odts
3259 xnc=max(0.,(nc1d(k) + ncten(k)*dtsave)*rho(k))
3260 if (xnc.gt.nt_c_max) &
3261 ncten(k) = (nt_c_max-nc1d(k)*rho(k))*odts*orho
3264 qiten(k) = qiten(k) + (pri_inu(k) + pri_iha(k) + pri_ihm(k) &
3265 + pri_wfz(k) + pri_rfz(k) + pri_ide(k) &
3266 - prs_iau(k) - prs_sci(k) - pri_rci(k)) &
3270 niten(k) = niten(k) + (pni_inu(k) + pni_iha(k) + pni_ihm(k) &
3271 + pni_wfz(k) + pni_rfz(k) + pni_ide(k) &
3272 - pni_iau(k) - pni_sci(k) - pni_rci(k)) &
3277 xri=max(r1,(qi1d(k) + qiten(k)*dtsave)*rho(k))
3278 xni=max(r2,(ni1d(k) + niten(k)*dtsave)*rho(k))
3279 if (xri.gt. r1)
then
3280 lami = (am_i*cig(2)*oig1*xni/xri)**obmi
3282 xdi = (bm_i + mu_i + 1.) * ilami
3283 if (xdi.lt. 5.e-6)
then
3285 xni = min(nt_i_max, cig(1)*oig2*xri/am_i*lami**bm_i)
3286 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
3287 elseif (xdi.gt. 300.e-6)
then
3288 lami = cie(2)/300.e-6
3289 xni = cig(1)*oig2*xri/am_i*lami**bm_i
3290 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
3293 niten(k) = -ni1d(k)*odts
3295 xni=max(0.,(ni1d(k) + niten(k)*dtsave)*rho(k))
3296 if (xni.gt.nt_i_max) &
3297 niten(k) = (nt_i_max-ni1d(k)*rho(k))*odts*orho
3300 qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) &
3301 + prr_sml(k) + prr_gml(k) + prr_rcs(k) &
3302 + prr_rcg(k) - prg_rfz(k) &
3303 - pri_rfz(k) - prr_rci(k)) &
3307 nrten(k) = nrten(k) + (pnr_wau(k) + pnr_sml(k) + pnr_gml(k) &
3308 - (pnr_rfz(k) + pnr_rcr(k) + pnr_rcg(k) &
3309 + pnr_rcs(k) + pnr_rci(k) + pni_rfz(k)) ) &
3314 xrr=max(r1,(qr1d(k) + qrten(k)*dtsave)*rho(k))
3315 xnr=max(r2,(nr1d(k) + nrten(k)*dtsave)*rho(k))
3316 if (xrr.gt. r1)
then
3317 lamr = (am_r*crg(3)*org2*xnr/xrr)**obmr
3318 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3319 if (mvd_r(k) .gt. 2.5e-3)
then
3321 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3322 xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
3323 nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
3324 elseif (mvd_r(k) .lt. d0r*0.75)
then
3326 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3327 xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
3328 nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
3331 qrten(k) = -qr1d(k)*odts
3332 nrten(k) = -nr1d(k)*odts
3336 qsten(k) = qsten(k) + (prs_iau(k) + prs_sde(k) &
3337 + prs_sci(k) + prs_scw(k) + prs_rcs(k) &
3338 + prs_ide(k) - prs_ihm(k) - prr_sml(k)) &
3342 qgten(k) = qgten(k) + (prg_scw(k) + prg_rfz(k) &
3343 + prg_gde(k) + prg_rcg(k) + prg_gcw(k) &
3344 + prg_rci(k) + prg_rcs(k) - prg_ihm(k) &
3349 if (temp(k).lt.t_0)
then
3351 + ( lsub*ocp(k)*(pri_inu(k) + pri_ide(k) &
3352 + prs_ide(k) + prs_sde(k) &
3353 + prg_gde(k) + pri_iha(k)) &
3354 + lfus2*ocp(k)*(pri_wfz(k) + pri_rfz(k) &
3355 + prg_rfz(k) + prs_scw(k) &
3356 + prg_scw(k) + prg_gcw(k) &
3357 + prg_rcs(k) + prs_rcs(k) &
3358 + prr_rci(k) + prg_rcg(k)) &
3362 + ( lfus*ocp(k)*(-prr_sml(k) - prr_gml(k) &
3363 - prr_rcg(k) - prr_rcs(k)) &
3364 + lsub*ocp(k)*(prs_sde(k) + prg_gde(k)) &
3374 temp(k) = t1d(k) + dt*tten(k)
3376 tempc = temp(k) - 273.15
3377 qv(k) = max(1.e-10, qv1d(k) + dt*qvten(k))
3378 rho(k) = roverrv*pres(k) / (r*temp(k)*(qv(k)+roverrv))
3379 rhof(k) = sqrt(rho_not/rho(k))
3380 rhof2(k) = sqrt(rhof(k))
3381 qvs(k) =
rslf(pres(k), temp(k))
3382 ssatw(k) = qv(k)/qvs(k) - 1.
3383 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
3384 diffu(k) = 2.11e-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
3385 if (tempc .ge. 0.0)
then
3386 visco(k) = (1.718+0.0049*tempc)*1.0e-5
3388 visco(k) = (1.718+0.0049*tempc-1.2e-5*tempc*tempc)*1.0e-5
3390 vsc2(k) = sqrt(rho(k)/visco(k))
3391 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
3392 tcond(k) = (5.69 + 0.0168*tempc)*1.0e-5 * 418.936
3393 ocp(k) = 1./(cp*(1.+0.887*qv(k)))
3394 lvt2(k)=lvap(k)*lvap(k)*ocp(k)*orv*otemp*otemp
3395 if (is_aerosol_aware) &
3396 nwfa(k) = max(11.1e6*rho(k), (nwfa1d(k) + nwfaten(k)*dt)*rho(k))
3400 if ((qc1d(k) + qcten(k)*dt) .gt. r1)
then
3401 rc(k) = (qc1d(k) + qcten(k)*dt)*rho(k)
3402 nc(k) = max(2., min((nc1d(k)+ncten(k)*dt)*rho(k), nt_c_max))
3403 if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware))
then
3417 if ((qi1d(k) + qiten(k)*dt) .gt. r1)
then
3418 ri(k) = (qi1d(k) + qiten(k)*dt)*rho(k)
3419 ni(k) = max(r2, (ni1d(k) + niten(k)*dt)*rho(k))
3427 if ((qr1d(k) + qrten(k)*dt) .gt. r1)
then
3428 rr(k) = (qr1d(k) + qrten(k)*dt)*rho(k)
3429 nr(k) = max(r2, (nr1d(k) + nrten(k)*dt)*rho(k))
3431 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3432 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3433 if (mvd_r(k) .gt. 2.5e-3)
then
3435 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3436 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
3437 elseif (mvd_r(k) .lt. d0r*0.75)
then
3439 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3440 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
3448 if ((qs1d(k) + qsten(k)*dt) .gt. r1)
then
3449 rs(k) = (qs1d(k) + qsten(k)*dt)*rho(k)
3456 if ((qg1d(k) + qgten(k)*dt) .gt. r1)
then
3457 rg(k) = (qg1d(k) + qgten(k)*dt)*rho(k)
3469 if (.not. iiwarm)
then
3477 if (.not. l_qs(k)) cycle
3478 tc0 = min(-0.1, temp(k)-273.15)
3479 smob(k) = rs(k)*oams
3483 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3))
then
3486 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
3487 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
3488 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
3489 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
3490 + sa(10)*bm_s*bm_s*bm_s
3492 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
3493 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
3494 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
3495 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
3496 + sb(10)*bm_s*bm_s*bm_s
3497 smo2(k) = (smob(k)/a_)**(1./b_)
3501 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
3502 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
3503 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
3504 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
3505 + sa(10)*cse(1)*cse(1)*cse(1)
3507 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
3508 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
3509 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
3510 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
3511 smoc(k) = a_ * smo2(k)**b_
3514 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(14) &
3515 + sa(4)*tc0*cse(14) + sa(5)*tc0*tc0 &
3516 + sa(6)*cse(14)*cse(14) + sa(7)*tc0*tc0*cse(14) &
3517 + sa(8)*tc0*cse(14)*cse(14) + sa(9)*tc0*tc0*tc0 &
3518 + sa(10)*cse(14)*cse(14)*cse(14)
3520 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(14) + sb(4)*tc0*cse(14) &
3521 + sb(5)*tc0*tc0 + sb(6)*cse(14)*cse(14) &
3522 + sb(7)*tc0*tc0*cse(14) + sb(8)*tc0*cse(14)*cse(14) &
3523 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(14)*cse(14)*cse(14)
3524 smod(k) = a_ * smo2(k)**b_
3537 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3539 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3540 n0_r(k) = nr(k)*org2*lamr**cre(2)
3553 if ( (ssatw(k).gt. eps) .or. (ssatw(k).lt. -eps .and. &
3555 clap = (qv(k)-qvs(k))/(1. + lvt2(k)*qvs(k))
3557 fcd = qvs(k)* exp(lvt2(k)*clap) - qv(k) + clap
3558 dfcd = qvs(k)*lvt2(k)* exp(lvt2(k)*clap) + 1.
3559 clap = clap - fcd/dfcd
3561 xrc = rc(k) + clap*rho(k)
3563 if (xrc.gt. r1)
then
3564 prw_vcd(k) = clap*odt
3566 if (clap .gt. eps)
then
3567 if (is_aerosol_aware .or. merra2_aerosol_aware)
then
3568 xnc = max(2.,
activ_ncloud(temp(k), w1d(k)+rand3, nwfa(k), lsml))
3576 pnc_wcd(k) = 0.5*(xnc-nc(k) + abs(xnc-nc(k)))*odts*orho
3579 elseif (clap .lt. -eps .AND. ssatw(k).lt.-1.e-6 .AND. &
3580 is_aerosol_aware)
then
3581 tempc = temp(k) - 273.15
3584 rvs_p = rvs*otemp*(lvap(k)*otemp*orv - 1.)
3585 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*orv - 1.) &
3586 *otemp*(lvap(k)*otemp*orv - 1.) &
3587 + (-2.*lvap(k)*otemp*otemp*otemp*orv) &
3589 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
3590 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
3591 * rvs_pp/rvs_p * rvs/rvs_p
3592 alphsc = max(1.e-9, alphsc)
3594 if (abs(xsat).lt. 1.e-9) xsat=0.
3595 t1_evap = 2.*pi*( 1.0 - alphsc*xsat &
3596 + 2.*alphsc*alphsc*xsat*xsat &
3597 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
3600 dc_star = sqrt(-2.0_dp*dt * t1_evap/(2.*pi) &
3601 * 4.*diffu(k)*ssatw(k)*rvs/rho_w)
3602 idx_d = max(1, min(int(1.e6*dc_star), nbc))
3604 idx_n = nint(1.0 + real(nbc, kind=wp) * log(real(nc(k)/t_nc(1), kind=dp)) / nic1)
3605 idx_n = max(1, min(idx_n, nbc))
3608 if (rc(k).gt. r_c(1))
then
3609 nic = nint(log10(rc(k)))
3610 do_loop_rc_cond:
do nn = nic-1, nic+1
3612 if ( (rc(k)/10.**nn).ge.1.0 .and. (rc(k)/10.**nn).lt.10.0 )
exit do_loop_rc_cond
3613 enddo do_loop_rc_cond
3614 idx_c = int(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
3615 idx_c = max(1, min(idx_c, ntb_c))
3622 prw_vcd(k) = max(real(-rc(k)*0.99*orho*odt, kind=dp), prw_vcd(k))
3623 pnc_wcd(k) = max(real(-nc(k)*0.99*orho*odt, kind=dp), &
3624 -tnc_wev(idx_d, idx_c, idx_n)*orho*odt)
3628 prw_vcd(k) = -rc(k)*orho*odt
3629 pnc_wcd(k) = -nc(k)*orho*odt
3634 qvten(k) = qvten(k) - prw_vcd(k)
3635 qcten(k) = qcten(k) + prw_vcd(k)
3636 ncten(k) = ncten(k) + pnc_wcd(k)
3637 if (is_aerosol_aware) &
3638 nwfaten(k) = nwfaten(k) - pnc_wcd(k)
3639 tten(k) = tten(k) + lvap(k)*ocp(k)*prw_vcd(k)*(1-ifdry)
3640 rc(k) = max(r1, (qc1d(k) + dt*qcten(k))*rho(k))
3641 if (rc(k).eq.r1) l_qc(k) = .false.
3642 nc(k) = max(2., min((nc1d(k)+ncten(k)*dt)*rho(k), nt_c_max))
3643 if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware))
then
3650 qv(k) = max(1.e-10, qv1d(k) + dt*qvten(k))
3651 temp(k) = t1d(k) + dt*tten(k)
3652 rho(k) = roverrv*pres(k) / (r*temp(k)*(qv(k)+roverrv))
3653 qvs(k) =
rslf(pres(k), temp(k))
3654 ssatw(k) = qv(k)/qvs(k) - 1.
3663 if ( (ssatw(k).lt. -eps) .and. l_qr(k) &
3664 .and. (.not.(prw_vcd(k).gt. 0.)) )
then
3665 tempc = temp(k) - 273.15
3668 rhof(k) = sqrt(rho_not*orho)
3669 rhof2(k) = sqrt(rhof(k))
3670 diffu(k) = 2.11e-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
3671 if (tempc .ge. 0.0)
then
3672 visco(k) = (1.718+0.0049*tempc)*1.0e-5
3674 visco(k) = (1.718+0.0049*tempc-1.2e-5*tempc*tempc)*1.0e-5
3676 vsc2(k) = sqrt(rho(k)/visco(k))
3677 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
3678 tcond(k) = (5.69 + 0.0168*tempc)*1.0e-5 * 418.936
3679 ocp(k) = 1./(cp*(1.+0.887*qv(k)))
3682 rvs_p = rvs*otemp*(lvap(k)*otemp*orv - 1.)
3683 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*orv - 1.) &
3684 *otemp*(lvap(k)*otemp*orv - 1.) &
3685 + (-2.*lvap(k)*otemp*otemp*otemp*orv) &
3687 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
3688 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
3689 * rvs_pp/rvs_p * rvs/rvs_p
3690 alphsc = max(1.e-9, alphsc)
3691 xsat = min(-1.e-9, ssatw(k))
3692 t1_evap = 2.*pi*( 1.0 - alphsc*xsat &
3693 + 2.*alphsc*alphsc*xsat*xsat &
3694 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
3699 if (qv(k)/qvs(k) .lt. 0.95 .AND. rr(k)*orho.le.1.e-8)
then
3700 prv_rev(k) = rr(k)*orho*odts
3702 prv_rev(k) = t1_evap*diffu(k)*(-ssatw(k))*n0_r(k)*rvs &
3703 * (t1_qr_ev*ilamr(k)**cre(10) &
3704 + t2_qr_ev*vsc2(k)*rhof2(k)*((lamr+0.5*fv_r)**(-cre(11))))
3705 rate_max = min((rr(k)*orho*odts), (qvs(k)-qv(k))*odts)
3706 prv_rev(k) = min(real(rate_max, kind=dp), prv_rev(k)*orho)
3715 if (prr_gml(k).gt.0.0)
then
3716 eva_factor = min(1.0, 0.01+(0.99-0.01)*(tempc/20.0))
3717 prv_rev(k) = prv_rev(k)*eva_factor
3721 pnr_rev(k) = min(real(nr(k)*0.99*orho*odts, kind=dp), &
3722 prv_rev(k) * nr(k)/rr(k))
3724 qrten(k) = qrten(k) - prv_rev(k)
3725 qvten(k) = qvten(k) + prv_rev(k)
3726 nrten(k) = nrten(k) - pnr_rev(k)
3727 if (is_aerosol_aware) &
3728 nwfaten(k) = nwfaten(k) + pnr_rev(k)
3729 tten(k) = tten(k) - lvap(k)*ocp(k)*prv_rev(k)*(1-ifdry)
3731 rr(k) = max(r1, (qr1d(k) + dt*qrten(k))*rho(k))
3732 qv(k) = max(1.e-10, qv1d(k) + dt*qvten(k))
3733 nr(k) = max(r2, (nr1d(k) + dt*nrten(k))*rho(k))
3734 temp(k) = t1d(k) + dt*tten(k)
3735 rho(k) = roverrv*pres(k) / (r*temp(k)*(qv(k)+roverrv))
3738#if ( WRF_CHEM == 1 )
3740 evapprod(k) = prv_rev(k) - (min(zerod0,prs_sde(k)) + &
3741 min(zerod0,prg_gde(k)))
3742 rainprod(k) = prr_wau(k) + prr_rcw(k) + prs_scw(k) + &
3743 prg_scw(k) + prs_iau(k) + &
3744 prg_gcw(k) + prs_sci(k) + &
3760 do k = kte+1, kts, -1
3771 if (any(l_qr .eqv. .true.))
then
3774 rhof(k) = sqrt(rho_not/rho(k))
3776 if (rr(k).gt. r1)
then
3777 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3778 vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) &
3779 *((lamr+fv_r)**(-cre(6)))
3786 vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) &
3787 *((lamr+fv_r)**(-cre(7)))
3791 vtnrk(k) = vtnrk(k+1)
3794 if (max(vtrk(k),vtnrk(k)) .gt. 1.e-3)
then
3795 ksed1(1) = max(ksed1(1), k)
3796 delta_tp = dzq(k)/(max(vtrk(k),vtnrk(k)))
3797 nstep = max(nstep, int(dt/delta_tp + 1.))
3800 if (ksed1(1) .eq. kte) ksed1(1) = kte-1
3801 if (nstep .gt. 0) onstep(1) = 1./real(nstep, kind=wp)
3806 if (any(l_qc .eqv. .true.))
then
3808 do_loop_hgt_agl :
do k = kts, kte-1
3809 if (rc(k) .gt. r2) ksed1(5) = k
3810 hgt_agl = hgt_agl + dzq(k)
3811 if (hgt_agl .gt. 500.0)
exit do_loop_hgt_agl
3812 enddo do_loop_hgt_agl
3814 do k = ksed1(5), kts, -1
3816 if (rc(k) .gt. r1 .and. w1d(k) .lt. 1.e-1)
then
3817 if (nc(k).gt.10000.e6)
then
3819 elseif (nc(k).lt.100.)
then
3822 nu_c = nint(1000.e6/nc(k)) + 2
3823 nu_c = max(2, min(nu_c+nint(rand2), 15))
3825 lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
3827 vtc = rhof(k)*av_c*ccg(5,nu_c)*ocg2(nu_c) * ilamc**bv_c
3829 vtc = rhof(k)*av_c*ccg(4,nu_c)*ocg1(nu_c) * ilamc**bv_c
3837 if (.not. iiwarm)
then
3839 if (any(l_qi .eqv. .true.))
then
3844 if (ri(k).gt. r1)
then
3845 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
3847 vti = rhof(k)*av_i*cig(3)*oig2 * ilami**bv_i
3852 vti = rhof(k)*av_i*cig(6)/cig(7) * ilami**bv_i
3856 vtnik(k) = vtnik(k+1)
3859 if (vtik(k) .gt. 1.e-3)
then
3860 ksed1(2) = max(ksed1(2), k)
3861 delta_tp = dzq(k)/vtik(k)
3862 nstep = max(nstep, int(dt/delta_tp + 1.))
3865 if (ksed1(2) .eq. kte) ksed1(2) = kte-1
3866 if (nstep .gt. 0) onstep(2) = 1./real(nstep, kind=wp)
3871 if (any(l_qs .eqv. .true.))
then
3877 if (rs(k).gt. r1)
then
3878 xds = smoc(k) / smob(k)
3880 ils1 = 1./(mrat*lam0 + fv_s)
3881 ils2 = 1./(mrat*lam1 + fv_s)
3882 t1_vts = kap0*csg(4)*ils1**cse(4)
3883 t2_vts = kap1*mrat**mu_s*csg(10)*ils2**cse(10)
3884 ils1 = 1./(mrat*lam0)
3885 ils2 = 1./(mrat*lam1)
3886 t3_vts = kap0*csg(1)*ils1**cse(1)
3887 t4_vts = kap1*mrat**mu_s*csg(7)*ils2**cse(7)
3888 vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts)
3889 if (prr_sml(k) .gt. 0.0)
then
3892 sr = rs(k)/(rs(k)+rr(k))
3893 vtsk(k) = vts*sr + (1.-sr)*vtrk(k)
3896 vtsk(k) = vts*vts_boost(k)
3904 if (vtsk(k) .gt. 1.e-3)
then
3905 ksed1(3) = max(ksed1(3), k)
3906 delta_tp = dzq(k)/vtsk(k)
3907 nstep = max(nstep, int(dt/delta_tp + 1.))
3910 if (ksed1(3) .eq. kte) ksed1(3) = kte-1
3911 if (nstep .gt. 0) onstep(3) = 1./real(nstep, kind=wp)
3916 if (any(l_qg .eqv. .true.))
then
3921 if (rg(k).gt. r1)
then
3922 vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
3923 if (temp(k).gt. t_0)
then
3924 vtgk(k) = max(vtg, vtrk(k))
3932 if (vtgk(k) .gt. 1.e-3)
then
3933 ksed1(4) = max(ksed1(4), k)
3934 delta_tp = dzq(k)/vtgk(k)
3935 nstep = max(nstep, int(dt/delta_tp + 1.))
3938 if (ksed1(4) .eq. kte) ksed1(4) = kte-1
3939 if (nstep .gt. 0) onstep(4) = 1./real(nstep, kind=wp)
3949 if (any(l_qr .eqv. .true.))
then
3950 nstep = nint(1./onstep(1))
3952 if(.not. sedi_semi)
then
3955 sed_r(k) = vtrk(k)*rr(k)
3956 sed_n(k) = vtnrk(k)*nr(k)
3961 qrten(k) = qrten(k) - sed_r(k)*odzq*onstep(1)*orho
3962 nrten(k) = nrten(k) - sed_n(k)*odzq*onstep(1)*orho
3963 rr(k) = max(r1, rr(k) - sed_r(k)*odzq*dt*onstep(1))
3964 nr(k) = max(r2, nr(k) - sed_n(k)*odzq*dt*onstep(1))
3965 pfll1(k) = pfll1(k) + sed_r(k)*dt*onstep(1)
3966 do k = ksed1(1), kts, -1
3969 qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) &
3970 *odzq*onstep(1)*orho
3971 nrten(k) = nrten(k) + (sed_n(k+1)-sed_n(k)) &
3972 *odzq*onstep(1)*orho
3973 rr(k) = max(r1, rr(k) + (sed_r(k+1)-sed_r(k)) &
3975 nr(k) = max(r2, nr(k) + (sed_n(k+1)-sed_n(k)) &
3977 pfll1(k) = pfll1(k) + sed_r(k)*dt*onstep(1)
3980 if (rr(kts).gt.r1*rr_min)
then
3981 pptrain = pptrain + sed_r(kts)*dt*onstep(1)
3987 niter = int(nstep/max(decfl,1)) + 1
3992 call semi_lagrange_sedim(kte,dzq,vtrk,rr,rainsfc,pfll,dtcfl,r1)
3993 call semi_lagrange_sedim(kte,dzq,vtnrk,nr,vtr,pdummy,dtcfl,r2)
3995 orhodt = 1./(rho(k)*dt)
3996 qrten(k) = qrten(k) + (rr(k) - rr_tmp(k)) * orhodt
3997 nrten(k) = nrten(k) + (nr(k) - nr_tmp(k)) * orhodt
3998 pfll1(k) = pfll1(k) + pfll(k)
4000 pptrain = pptrain + rainsfc
4002 do k = kte+1, kts, -1
4008 if (rr(k).gt. r1)
then
4009 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
4010 vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) &
4011 *((lamr+fv_r)**(-cre(6)))
4018 vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) &
4019 *((lamr+fv_r)**(-cre(7)))
4029 if (any(l_qc .eqv. .true.))
then
4031 sed_c(k) = vtck(k)*rc(k)
4032 sed_n(k) = vtnck(k)*nc(k)
4034 do k = ksed1(5), kts, -1
4037 qcten(k) = qcten(k) + (sed_c(k+1)-sed_c(k)) *odzq*orho
4038 ncten(k) = ncten(k) + (sed_n(k+1)-sed_n(k)) *odzq*orho
4039 rc(k) = max(r1, rc(k) + (sed_c(k+1)-sed_c(k)) *odzq*dt)
4040 nc(k) = max(10., nc(k) + (sed_n(k+1)-sed_n(k)) *odzq*dt)
4046 if (any(l_qi .eqv. .true.))
then
4047 nstep = nint(1./onstep(2))
4050 sed_i(k) = vtik(k)*ri(k)
4051 sed_n(k) = vtnik(k)*ni(k)
4056 qiten(k) = qiten(k) - sed_i(k)*odzq*onstep(2)*orho
4057 niten(k) = niten(k) - sed_n(k)*odzq*onstep(2)*orho
4058 ri(k) = max(r1, ri(k) - sed_i(k)*odzq*dt*onstep(2))
4059 ni(k) = max(r2, ni(k) - sed_n(k)*odzq*dt*onstep(2))
4060 pfil1(k) = pfil1(k) + sed_i(k)*dt*onstep(2)
4061 do k = ksed1(2), kts, -1
4064 qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) &
4065 *odzq*onstep(2)*orho
4066 niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) &
4067 *odzq*onstep(2)*orho
4068 ri(k) = max(r1, ri(k) + (sed_i(k+1)-sed_i(k)) &
4070 ni(k) = max(r2, ni(k) + (sed_n(k+1)-sed_n(k)) &
4072 pfil1(k) = pfil1(k) + sed_i(k)*dt*onstep(2)
4075 if (ri(kts).gt.r1*rr_min)
then
4076 pptice = pptice + sed_i(kts)*dt*onstep(2)
4083 if (any(l_qs .eqv. .true.))
then
4084 nstep = nint(1./onstep(3))
4087 sed_s(k) = vtsk(k)*rs(k)
4092 qsten(k) = qsten(k) - sed_s(k)*odzq*onstep(3)*orho
4093 rs(k) = max(r1, rs(k) - sed_s(k)*odzq*dt*onstep(3))
4094 pfil1(k) = pfil1(k) + sed_s(k)*dt*onstep(3)
4095 do k = ksed1(3), kts, -1
4098 qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) &
4099 *odzq*onstep(3)*orho
4100 rs(k) = max(r1, rs(k) + (sed_s(k+1)-sed_s(k)) &
4102 pfil1(k) = pfil1(k) + sed_s(k)*dt*onstep(3)
4105 if (rs(kts).gt.r1*rr_min)
then
4106 pptsnow = pptsnow + sed_s(kts)*dt*onstep(3)
4113 if (any(l_qg .eqv. .true.))
then
4114 nstep = nint(1./onstep(4))
4115 if(.not. sedi_semi)
then
4118 sed_g(k) = vtgk(k)*rg(k)
4123 qgten(k) = qgten(k) - sed_g(k)*odzq*onstep(4)*orho
4124 rg(k) = max(r1, rg(k) - sed_g(k)*odzq*dt*onstep(4))
4125 pfil1(k) = pfil1(k) + sed_g(k)*dt*onstep(4)
4126 do k = ksed1(4), kts, -1
4129 qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) &
4130 *odzq*onstep(4)*orho
4131 rg(k) = max(r1, rg(k) + (sed_g(k+1)-sed_g(k)) &
4133 pfil1(k) = pfil1(k) + sed_g(k)*dt*onstep(4)
4136 if (rg(kts).gt.r1*rr_min)
then
4137 pptgraul = pptgraul + sed_g(kts)*dt*onstep(4)
4143 niter = int(nstep/max(decfl,1)) + 1
4148 call semi_lagrange_sedim(kte,dzq,vtgk,rg,graulsfc,pfil,dtcfl,r1)
4150 orhodt = 1./(rho(k)*dt)
4151 qgten(k) = qgten(k) + (rg(k) - rg_tmp(k))*orhodt
4152 pfil1(k) = pfil1(k) + pfil(k)
4154 pptgraul = pptgraul + graulsfc
4155 do k = kte+1, kts, -1
4160 if (rg(k).gt. r1)
then
4161 ygra1 = log10(max(1.e-9_wp, rg(k)))
4162 zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1
4163 n0_exp = 10.**(zans1)
4164 n0_exp = max(real(gonv_min, kind=dp), min(n0_exp, real(gonv_max, kind=dp)))
4165 lam_exp = (n0_exp*am_g*cgg(1)/rg(k))**oge1
4166 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
4168 vtg = rhof(k)*av_g*cgg(6)*ogg3 * (1./lamg)**bv_g
4169 if (temp(k).gt. t_0)
then
4170 vtgk(k) = max(vtg, vtrk(k))
4184 if (.not. iiwarm)
then
4186 xri = max(0.0, qi1d(k) + qiten(k)*dt)
4187 if ( (temp(k).gt. t_0) .and. (xri.gt. 0.0) )
then
4188 qcten(k) = qcten(k) + xri*odt
4189 ncten(k) = ncten(k) + ni1d(k)*odt
4190 qiten(k) = qiten(k) - xri*odt
4191 niten(k) = -ni1d(k)*odt
4192 tten(k) = tten(k) - lfus*ocp(k)*xri*odt*(1-ifdry)
4197 xrc = max(0.0, qc1d(k) + qcten(k)*dt)
4198 if ( (temp(k).lt. hgfr) .and. (xrc.gt. 0.0) )
then
4199 lfus2 = lsub - lvap(k)
4200 xnc = nc1d(k) + ncten(k)*dt
4201 qiten(k) = qiten(k) + xrc*odt
4202 niten(k) = niten(k) + xnc*odt
4203 qcten(k) = qcten(k) - xrc*odt
4204 ncten(k) = ncten(k) - xnc*odt
4205 tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-ifdry)
4216 t1d(k) = t1d(k) + tten(k)*dt
4217 qv1d(k) = max(1.e-10, qv1d(k) + qvten(k)*dt)
4218 qc1d(k) = qc1d(k) + qcten(k)*dt
4219 nc1d(k) = max(2./rho(k), min(nc1d(k) + ncten(k)*dt, nt_c_max))
4220 if (is_aerosol_aware)
then
4221 nwfa1d(k) = max(11.1e6, min(9999.e6, &
4222 (nwfa1d(k)+nwfaten(k)*dt)))
4223 nifa1d(k) = max(nain1*0.01, min(9999.e6, &
4224 (nifa1d(k)+nifaten(k)*dt)))
4226 if (qc1d(k) .le. r1)
then
4230 if (nc1d(k)*rho(k).gt.10000.e6)
then
4232 elseif (nc1d(k)*rho(k).lt.100.)
then
4235 nu_c = nint(1000.e6/(nc1d(k)*rho(k))) + 2
4236 nu_c = max(2, min(nu_c+nint(rand2), 15))
4238 lamc = (am_r*ccg(2,nu_c)*ocg1(nu_c)*nc1d(k)/qc1d(k))**obmr
4239 xdc = (bm_r + nu_c + 1.) / lamc
4240 if (xdc.lt. d0c)
then
4241 lamc = cce(2,nu_c)/d0c
4242 elseif (xdc.gt. d0r*2.)
then
4243 lamc = cce(2,nu_c)/(d0r*2.)
4245 nc1d(k) = min(ccg(1,nu_c)*ocg2(nu_c)*qc1d(k)/am_r*lamc**bm_r,&
4246 real(Nt_c_max, kind=dp)/rho(k))
4249 qi1d(k) = qi1d(k) + qiten(k)*dt
4250 ni1d(k) = max(r2/rho(k), ni1d(k) + niten(k)*dt)
4251 if (qi1d(k) .le. r1)
then
4255 lami = (am_i*cig(2)*oig1*ni1d(k)/qi1d(k))**obmi
4257 xdi = (bm_i + mu_i + 1.) * ilami
4258 if (xdi.lt. 5.e-6)
then
4260 elseif (xdi.gt. 300.e-6)
then
4261 lami = cie(2)/300.e-6
4263 ni1d(k) = min(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, &
4266 qr1d(k) = qr1d(k) + qrten(k)*dt
4267 nr1d(k) = max(r2/rho(k), nr1d(k) + nrten(k)*dt)
4268 if (qr1d(k) .le. r1)
then
4272 lamr = (am_r*crg(3)*org2*nr1d(k)/qr1d(k))**obmr
4273 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
4274 if (mvd_r(k) .gt. 2.5e-3)
then
4276 elseif (mvd_r(k) .lt. d0r*0.75)
then
4279 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
4280 nr1d(k) = crg(2)*org3*qr1d(k)*lamr**bm_r / am_r
4282 qs1d(k) = qs1d(k) + qsten(k)*dt
4283 if (qs1d(k) .le. r1) qs1d(k) = 0.0
4284 qg1d(k) = qg1d(k) + qgten(k)*dt
4285 if (qg1d(k) .le. r1) qg1d(k) = 0.0
4289 calculate_extended_diagnostics:
if (ext_diag)
then
4291 if(prw_vcd(k).gt.0)
then
4292 prw_vcdc1(k) = prw_vcd(k)*dt
4293 elseif(prw_vcd(k).lt.0)
then
4294 prw_vcde1(k) = -1*prw_vcd(k)*dt
4297 tpri_inu1(k) = pri_inu(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4299 if(pri_ide(k).gt.0)
then
4300 tpri_ide1_d(k) = pri_ide(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4302 tpri_ide1_s(k) = -pri_ide(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4305 if(temp(k).lt.t_0)
then
4306 tprs_ide1(k) = prs_ide(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4309 if(prs_sde(k).gt.0)
then
4310 tprs_sde1_d(k) = prs_sde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4312 tprs_sde1_s(k) = -prs_sde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4315 if(prg_gde(k).gt.0)
then
4316 tprg_gde1_d(k) = prg_gde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4318 tprg_gde1_s(k) = -prg_gde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4321 tpri_iha1(k) = pri_iha(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4322 tpri_wfz1(k) = pri_wfz(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4323 tpri_rfz1(k) = pri_rfz(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4324 tprg_rfz1(k) = prg_rfz(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4325 tprs_scw1(k) = prs_scw(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4326 tprg_scw1(k) = prg_scw(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4327 tprg_rcs1(k) = prg_rcs(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4329 if(temp(k).lt.t_0)
then
4330 tprs_rcs1(k) = prs_rcs(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4333 tprr_rci1(k) = prr_rci(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4335 if(temp(k).lt.t_0)
then
4336 tprg_rcg1(k) = prg_rcg(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4339 if(prw_vcd(k).gt.0)
then
4340 tprw_vcd1_c(k) = lvap(k)*ocp(k)*prw_vcd(k)*(1-ifdry)*dt
4342 tprw_vcd1_e(k) = -lvap(k)*ocp(k)*prw_vcd(k)*(1-ifdry)*dt
4346 tprr_sml1(k) = prr_sml(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
4347 tprr_gml1(k) = prr_gml(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
4349 if(temp(k).ge.t_0)
then
4350 tprr_rcg1(k) = -prr_rcg(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
4353 if(temp(k).ge.t_0)
then
4354 tprr_rcs1(k) = -prr_rcs(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
4357 tprv_rev1(k) = lvap(k)*ocp(k)*prv_rev(k)*(1-ifdry)*dt
4358 tten1(k) = tten(k)*dt
4359 qvten1(k) = qvten(k)*dt
4360 qiten1(k) = qiten(k)*dt
4361 qrten1(k) = qrten(k)*dt
4362 qsten1(k) = qsten(k)*dt
4363 qgten1(k) = qgten(k)*dt
4364 niten1(k) = niten(k)*dt
4365 nrten1(k) = nrten(k)*dt
4366 ncten1(k) = ncten(k)*dt
4367 qcten1(k) = qcten(k)*dt
4369 endif calculate_extended_diagnostics