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, &
1047 fs_fac_rain, fs_fac_snow)
1052 integer,
intent(in):: ids,ide, jds,jde, kds,kde, &
1053 ims,ime, jms,jme, kms,kme, &
1054 its,ite, jts,jte, kts,kte
1055 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
intent(inout):: &
1056 qv, qc, qr, qi, qs, qg, ni, nr
1057 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
optional,
intent(inout):: &
1059 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
optional,
intent(in):: &
1061 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
optional,
intent(inout):: &
1063 real(wp),
dimension(ims:ime, jms:jme),
optional,
intent(in):: nwfa2d, nifa2d
1064 integer,
dimension(ims:ime, jms:jme),
intent(in):: lsm
1065 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
optional,
intent(inout):: &
1066 re_cloud, re_ice, re_snow
1067 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
intent(inout):: pfils, pflls
1068 integer,
intent(in) :: rand_perturb_on, kme_stoch, n_var_spp
1069 real(wp),
dimension(:,:),
intent(in),
optional :: rand_pert
1070 real(wp),
dimension(:),
intent(in),
optional :: spp_prt_list, spp_stddev_cutoff
1071 character(len=10),
dimension(:),
intent(in),
optional :: spp_var_list
1072 integer,
intent(in):: has_reqc, has_reqi, has_reqs
1073#if ( WRF_CHEM == 1 )
1074 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
intent(inout):: &
1077 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
intent(in):: &
1079 real(wp),
dimension(ims:ime, jms:jme),
intent(inout):: &
1081 real(wp),
dimension(ims:ime, jms:jme),
optional,
intent(inout):: &
1084 GRAUPELNC, GRAUPELNCV
1085 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
intent(inout):: &
1087 real(wp),
dimension(ims:ime, jms:jme),
intent(inout):: &
1089 real(wp),
dimension(ims:ime, kms:kme, jms:jme),
optional,
intent(inout):: &
1091 logical,
intent(in) :: first_time_step
1092 real(wp),
intent(in):: dt_in, dt_inner
1093 logical,
intent(in) :: sedi_semi
1094 integer,
intent(in) :: decfl
1096 integer,
intent (in) :: istep, nsteps
1097 logical,
intent (in) :: fullradar_diag
1099 logical,
intent (in) :: ext_diag
1100 logical,
optional,
intent(in):: aero_ind_fdb
1101 real(wp),
dimension(:,:,:),
optional,
intent(inout):: &
1104 prw_vcde, tpri_inu, tpri_ide_d, &
1105 tpri_ide_s, tprs_ide, &
1106 tprs_sde_d, tprs_sde_s, tprg_gde_d, &
1107 tprg_gde_s, tpri_iha, tpri_wfz, &
1108 tpri_rfz, tprg_rfz, tprs_scw, tprg_scw, &
1109 tprg_rcs, tprs_rcs, &
1110 tprr_rci, tprg_rcg, &
1111 tprw_vcd_c, tprw_vcd_e, tprr_sml, &
1112 tprr_gml, tprr_rcg, &
1113 tprr_rcs, tprv_rev, tten3, qvten3, &
1114 qrten3, qsten3, qgten3, qiten3, niten3, &
1115 nrten3, ncten3, qcten3
1117 real(wp),
INTENT (IN),
optional :: fs_fac_rain, fs_fac_snow
1120 real(wp),
dimension(kts:kte):: &
1121 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
1122 nr1d, nc1d, nwfa1d, nifa1d, &
1123 t1d, p1d, w1d, dz1d, rho, dBZ, pfil1, pfll1
1125 real(wp),
dimension(:),
allocatable:: &
1128 prw_vcde1, tpri_inu1, tpri_ide1_d, &
1129 tpri_ide1_s, tprs_ide1, &
1130 tprs_sde1_d, tprs_sde1_s, tprg_gde1_d, &
1131 tprg_gde1_s, tpri_iha1, tpri_wfz1, &
1132 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1,&
1133 tprg_rcs1, tprs_rcs1, &
1134 tprr_rci1, tprg_rcg1, &
1135 tprw_vcd1_c, tprw_vcd1_e, tprr_sml1, &
1136 tprr_gml1, tprr_rcg1, &
1137 tprr_rcs1, tprv_rev1, tten1, qvten1, &
1138 qrten1, qsten1, qgten1, qiten1, niten1, &
1139 nrten1, ncten1, qcten1
1141 real(wp),
dimension(kts:kte):: re_qc1d, re_qi1d, re_qs1d
1142#if ( WRF_CHEM == 1 )
1143 real(wp),
dimension(kts:kte):: &
1144 rainprod1d, evapprod1d
1146 real(wp),
dimension(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
1147 real(wp) :: dt, pptrain, pptsnow, pptgraul, pptice
1148 real(wp) :: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max
1150 real(wp) :: rand1, rand2, rand3, rand_pert_max
1151 integer:: i, j, k, m
1152 integer:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr
1153 integer:: jmax_qc,jmax_qr,jmax_qi,jmax_qs,jmax_qg,jmax_ni,jmax_nr
1154 integer:: kmax_qc,kmax_qr,kmax_qi,kmax_qs,kmax_qg,kmax_ni,kmax_nr
1155 integer:: i_start, j_start, i_end, j_end
1156 logical,
optional,
intent(in) :: diagflag
1157 integer,
optional,
intent(in) :: do_radar_ref
1158 logical :: melti = .false.
1162 character(len=*),
optional,
intent( out) :: errmsg
1163 integer,
optional,
intent( out) :: errflg
1166 if (
present(errmsg)) errmsg =
''
1167 if (
present(errflg)) errflg = 0
1170 test_only_once:
if (first_time_step .and. istep==1)
then
1173 if ( (
present(tt) .and. (
present(th) .or.
present(pii))) .or. &
1174 (.not.
present(tt) .and. .not.(
present(th) .and.
present(pii))) )
then
1175 if (
present(errmsg) .and.
present(errflg))
then
1176 write(errmsg,
'(a)')
'Logic error in mp_gt_driver: provide either tt or th+pii'
1180 write(*,
'(a)')
'Logic error in mp_gt_driver: provide either tt or th+pii'
1185 if (is_aerosol_aware .and. (.not.
present(nc) .or. &
1186 .not.
present(nwfa) .or. &
1187 .not.
present(nifa) .or. &
1188 .not.
present(nwfa2d) .or. &
1189 .not.
present(nifa2d) ))
then
1190 if (
present(errmsg) .and.
present(errflg))
then
1191 write(errmsg,
'(*(a))')
'Logic error in mp_gt_driver: provide nc, nwfa, nifa, nwfa2d', &
1192 ' and nifa2d for aerosol-aware version of Thompson microphysics'
1196 write(*,
'(*(a))')
'Logic error in mp_gt_driver: provide nc, nwfa, nifa, nwfa2d', &
1197 ' and nifa2d for aerosol-aware version of Thompson microphysics'
1200 else if (merra2_aerosol_aware .and. (.not.
present(nc) .or. &
1201 .not.
present(nwfa) .or. &
1202 .not.
present(nifa) ))
then
1203 if (
present(errmsg) .and.
present(errflg))
then
1204 write(errmsg,
'(*(a))')
'Logic error in mp_gt_driver: provide nc, nwfa, and nifa', &
1205 ' for merra2 aerosol-aware version of Thompson microphysics'
1209 write(*,
'(*(a))')
'Logic error in mp_gt_driver: provide nc, nwfa, and nifa', &
1210 ' for merra2 aerosol-aware version of Thompson microphysics'
1213 else if (.not.is_aerosol_aware .and. .not.merra2_aerosol_aware .and. &
1214 (
present(nwfa) .or.
present(nifa) .or.
present(nwfa2d) .or.
present(nifa2d)))
then
1215 write(*,*)
'WARNING, nc/nwfa/nifa/nwfa2d/nifa2d present but is_aerosol_aware/merra2_aerosol_aware are FALSE'
1217 end if test_only_once
1223 allocate_extended_diagnostics:
if (ext_diag)
then
1224 allocate (prw_vcdc1(kts:kte))
1225 allocate (prw_vcde1(kts:kte))
1226 allocate (tpri_inu1(kts:kte))
1227 allocate (tpri_ide1_d(kts:kte))
1228 allocate (tpri_ide1_s(kts:kte))
1229 allocate (tprs_ide1(kts:kte))
1230 allocate (tprs_sde1_d(kts:kte))
1231 allocate (tprs_sde1_s(kts:kte))
1232 allocate (tprg_gde1_d(kts:kte))
1233 allocate (tprg_gde1_s(kts:kte))
1234 allocate (tpri_iha1(kts:kte))
1235 allocate (tpri_wfz1(kts:kte))
1236 allocate (tpri_rfz1(kts:kte))
1237 allocate (tprg_rfz1(kts:kte))
1238 allocate (tprs_scw1(kts:kte))
1239 allocate (tprg_scw1(kts:kte))
1240 allocate (tprg_rcs1(kts:kte))
1241 allocate (tprs_rcs1(kts:kte))
1242 allocate (tprr_rci1(kts:kte))
1243 allocate (tprg_rcg1(kts:kte))
1244 allocate (tprw_vcd1_c(kts:kte))
1245 allocate (tprw_vcd1_e(kts:kte))
1246 allocate (tprr_sml1(kts:kte))
1247 allocate (tprr_gml1(kts:kte))
1248 allocate (tprr_rcg1(kts:kte))
1249 allocate (tprr_rcs1(kts:kte))
1250 allocate (tprv_rev1(kts:kte))
1251 allocate (tten1(kts:kte))
1252 allocate (qvten1(kts:kte))
1253 allocate (qrten1(kts:kte))
1254 allocate (qsten1(kts:kte))
1255 allocate (qgten1(kts:kte))
1256 allocate (qiten1(kts:kte))
1257 allocate (niten1(kts:kte))
1258 allocate (nrten1(kts:kte))
1259 allocate (ncten1(kts:kte))
1260 allocate (qcten1(kts:kte))
1261 end if allocate_extended_diagnostics
1282 graupelnc(:,:) = 0.0
1290 ndt = max(nint(dt_in/dt_inner),1)
1292 if(dt_in .le. dt_inner) dt= dt_in
1297 if (rand_perturb_on .ne. 0)
then
1299 select case (spp_var_list(k))
1301 rand_pert_max = spp_prt_list(k)*spp_stddev_cutoff(k)
1337 j_loop:
do j = j_start, j_end
1338 i_loop:
do i = i_start, i_end
1356 if (rand_perturb_on .ne. 0)
then
1357 if (mod(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,1)
1358 m = rshift(abs(rand_perturb_on),1)
1359 if (mod(m,2) .ne. 0) rand2 = rand_pert(i,1)*2.
1360 m = rshift(abs(rand_perturb_on),2)
1361 if (mod(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,1)+rand_pert_max)
1362 m = rshift(abs(rand_perturb_on),3)
1371 IF (
PRESENT (snowncv) )
THEN
1374 IF (
PRESENT (icencv) )
THEN
1377 IF (
PRESENT (graupelncv) )
THEN
1378 graupelncv(i,j) = 0.
1383 if (
present(tt))
then
1386 t1d(k) = th(i,k,j)*pii(i,k,j)
1399 rho(k) = roverrv*p1d(k) / (r*t1d(k)*(qv1d(k)+roverrv))
1405 initialize_extended_diagnostics:
if (ext_diag)
then
1443 endif initialize_extended_diagnostics
1447 if (is_aerosol_aware .or. merra2_aerosol_aware)
then
1450 nwfa1d(k) = nwfa(i,k,j)
1451 nifa1d(k) = nifa(i,k,j)
1456 nc1d(k) = nt_c_l/rho(k)
1458 nc1d(k) = nt_c_o/rho(k)
1461 nifa1d(k) = nain1*0.01
1466 call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
1467 nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dz1d, &
1468 lsml, pptrain, pptsnow, pptgraul, pptice, &
1469#if ( WRF_CHEM == 1 )
1470 rainprod1d, evapprod1d, &
1472 rand1, rand2, rand3, &
1473 kts, kte, dt, i, j, ext_diag, &
1476 prw_vcdc1, prw_vcde1, &
1477 tpri_inu1, tpri_ide1_d, tpri_ide1_s, tprs_ide1, &
1478 tprs_sde1_d, tprs_sde1_s, &
1479 tprg_gde1_d, tprg_gde1_s, tpri_iha1, tpri_wfz1, &
1480 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1, &
1481 tprg_rcs1, tprs_rcs1, tprr_rci1, &
1482 tprg_rcg1, tprw_vcd1_c, &
1483 tprw_vcd1_e, tprr_sml1, tprr_gml1, tprr_rcg1, &
1484 tprr_rcs1, tprv_rev1, &
1485 tten1, qvten1, qrten1, qsten1, &
1486 qgten1, qiten1, niten1, nrten1, ncten1, qcten1, &
1488 fs_fac_rain, fs_fac_snow)
1490 pcp_ra(i,j) = pcp_ra(i,j) + pptrain
1491 pcp_sn(i,j) = pcp_sn(i,j) + pptsnow
1492 pcp_gr(i,j) = pcp_gr(i,j) + pptgraul
1493 pcp_ic(i,j) = pcp_ic(i,j) + pptice
1494 rainncv(i,j) = pptrain + pptsnow + pptgraul + pptice
1495 rainnc(i,j) = rainnc(i,j) + pptrain + pptsnow + pptgraul + pptice
1496 IF (
PRESENT(snowncv) .AND.
PRESENT(snownc) )
THEN
1498 IF ( .NOT.
PRESENT(icencv) .OR. .NOT.
PRESENT(icenc) )
THEN
1499 snowncv(i,j) = pptsnow + pptice
1500 snownc(i,j) = snownc(i,j) + pptsnow + pptice
1502 snowncv(i,j) = pptsnow
1503 snownc(i,j) = snownc(i,j) + pptsnow
1507 IF (
PRESENT(icencv) .AND.
PRESENT(icenc) )
THEN
1508 icencv(i,j) = pptice
1509 icenc(i,j) = icenc(i,j) + pptice
1511 IF (
PRESENT(graupelncv) .AND.
PRESENT(graupelnc) )
THEN
1512 graupelncv(i,j) = pptgraul
1513 graupelnc(i,j) = graupelnc(i,j) + pptgraul
1515 sr(i,j) = (pptsnow + pptgraul + pptice) / (rainncv(i,j)+r1)
1520 if (is_aerosol_aware)
then
1521 if (
PRESENT (aero_ind_fdb) )
then
1522 if ( .not. aero_ind_fdb)
then
1523 nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt
1524 nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt
1527 nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt
1528 nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt
1533 nwfa(i,k,j) = nwfa1d(k)
1534 nifa(i,k,j) = nifa1d(k)
1538 if (merra2_aerosol_aware)
then
1541 nwfa(i,k,j) = nwfa1d(k)
1542 nifa(i,k,j) = nifa1d(k)
1555 pfils(i,k,j) = pfils(i,k,j) + pfil1(k)
1556 pflls(i,k,j) = pflls(i,k,j) + pfll1(k)
1557 if (
present(tt))
then
1560 th(i,k,j) = t1d(k)/pii(i,k,j)
1562#if ( WRF_CHEM == 1 )
1563 rainprod(i,k,j) = rainprod1d(k)
1564 evapprod(i,k,j) = evapprod1d(k)
1566 if (qc1d(k) .gt. qc_max)
then
1571 elseif (qc1d(k) .lt. 0.0)
then
1572 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qc ', qc1d(k), &
1575 if (qr1d(k) .gt. qr_max)
then
1580 elseif (qr1d(k) .lt. 0.0)
then
1581 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qr ', qr1d(k), &
1584 if (nr1d(k) .gt. nr_max)
then
1589 elseif (nr1d(k) .lt. 0.0)
then
1590 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative nr ', nr1d(k), &
1593 if (qs1d(k) .gt. qs_max)
then
1598 elseif (qs1d(k) .lt. 0.0)
then
1599 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qs ', qs1d(k), &
1602 if (qi1d(k) .gt. qi_max)
then
1607 elseif (qi1d(k) .lt. 0.0)
then
1608 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qi ', qi1d(k), &
1611 if (qg1d(k) .gt. qg_max)
then
1616 elseif (qg1d(k) .lt. 0.0)
then
1617 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qg ', qg1d(k), &
1620 if (ni1d(k) .gt. ni_max)
then
1625 elseif (ni1d(k) .lt. 0.0)
then
1626 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative ni ', ni1d(k), &
1629 if (qv1d(k) .lt. 0.0)
then
1630 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qv ', qv1d(k), &
1632 if (k.lt.kte-2 .and. k.gt.kts+1)
then
1633 write(*,*)
' below and above are: ', qv(i,k-1,j), qv(i,k+1,j)
1634 qv(i,k,j) = max(1.e-7, 0.5*(qv(i,k-1,j) + qv(i,k+1,j)))
1641 assign_extended_diagnostics:
if (ext_diag)
then
1646 prw_vcdc(i,k,j) = prw_vcdc(i,k,j) + prw_vcdc1(k)
1647 prw_vcde(i,k,j) = prw_vcde(i,k,j) + prw_vcde1(k)
1648 tpri_inu(i,k,j) = tpri_inu(i,k,j) + tpri_inu1(k)
1649 tpri_ide_d(i,k,j) = tpri_ide_d(i,k,j) + tpri_ide1_d(k)
1650 tpri_ide_s(i,k,j) = tpri_ide_s(i,k,j) + tpri_ide1_s(k)
1651 tprs_ide(i,k,j) = tprs_ide(i,k,j) + tprs_ide1(k)
1652 tprs_sde_s(i,k,j) = tprs_sde_s(i,k,j) + tprs_sde1_s(k)
1653 tprs_sde_d(i,k,j) = tprs_sde_d(i,k,j) + tprs_sde1_d(k)
1654 tprg_gde_d(i,k,j) = tprg_gde_d(i,k,j) + tprg_gde1_d(k)
1655 tprg_gde_s(i,k,j) = tprg_gde_s(i,k,j) + tprg_gde1_s(k)
1656 tpri_iha(i,k,j) = tpri_iha(i,k,j) + tpri_iha1(k)
1657 tpri_wfz(i,k,j) = tpri_wfz(i,k,j) + tpri_wfz1(k)
1658 tpri_rfz(i,k,j) = tpri_rfz(i,k,j) + tpri_rfz1(k)
1659 tprg_rfz(i,k,j) = tprg_rfz(i,k,j) + tprg_rfz1(k)
1660 tprs_scw(i,k,j) = tprs_scw(i,k,j) + tprs_scw1(k)
1661 tprg_scw(i,k,j) = tprg_scw(i,k,j) + tprg_scw1(k)
1662 tprg_rcs(i,k,j) = tprg_rcs(i,k,j) + tprg_rcs1(k)
1663 tprs_rcs(i,k,j) = tprs_rcs(i,k,j) + tprs_rcs1(k)
1664 tprr_rci(i,k,j) = tprr_rci(i,k,j) + tprr_rci1(k)
1665 tprg_rcg(i,k,j) = tprg_rcg(i,k,j) + tprg_rcg1(k)
1666 tprw_vcd_c(i,k,j) = tprw_vcd_c(i,k,j) + tprw_vcd1_c(k)
1667 tprw_vcd_e(i,k,j) = tprw_vcd_e(i,k,j) + tprw_vcd1_e(k)
1668 tprr_sml(i,k,j) = tprr_sml(i,k,j) + tprr_sml1(k)
1669 tprr_gml(i,k,j) = tprr_gml(i,k,j) + tprr_gml1(k)
1670 tprr_rcg(i,k,j) = tprr_rcg(i,k,j) + tprr_rcg1(k)
1671 tprr_rcs(i,k,j) = tprr_rcs(i,k,j) + tprr_rcs1(k)
1672 tprv_rev(i,k,j) = tprv_rev(i,k,j) + tprv_rev1(k)
1673 tten3(i,k,j) = tten3(i,k,j) + tten1(k)
1674 qvten3(i,k,j) = qvten3(i,k,j) + qvten1(k)
1675 qrten3(i,k,j) = qrten3(i,k,j) + qrten1(k)
1676 qsten3(i,k,j) = qsten3(i,k,j) + qsten1(k)
1677 qgten3(i,k,j) = qgten3(i,k,j) + qgten1(k)
1678 qiten3(i,k,j) = qiten3(i,k,j) + qiten1(k)
1679 niten3(i,k,j) = niten3(i,k,j) + niten1(k)
1680 nrten3(i,k,j) = nrten3(i,k,j) + nrten1(k)
1681 ncten3(i,k,j) = ncten3(i,k,j) + ncten1(k)
1682 qcten3(i,k,j) = qcten3(i,k,j) + qcten1(k)
1684 endif assign_extended_diagnostics
1686 if (ndt>1 .and. it==ndt)
then
1687 sr(i,j) = (pcp_sn(i,j) + pcp_gr(i,j) + pcp_ic(i,j)) / (rainnc(i,j)+r1)
1688 rainncv(i,j) = rainnc(i,j)
1689 IF (
PRESENT (snowncv) )
THEN
1690 snowncv(i,j) = snownc(i,j)
1692 IF (
PRESENT (icencv) )
THEN
1693 icencv(i,j) = icenc(i,j)
1695 IF (
PRESENT (graupelncv) )
THEN
1696 graupelncv(i,j) = graupelnc(i,j)
1702 last_step_only:
IF ((ndt>1 .and. it==ndt) .or. &
1703 (nsteps>1 .and. istep==nsteps) .or. &
1704 (nsteps==1 .and. ndt==1))
THEN
1710 diagflag_present:
IF (
PRESENT (diagflag) )
THEN
1711 if (diagflag .and. do_radar_ref == 1)
then
1714 if (fullradar_diag)
then
1720 if (
present(vt_dbz_wt))
then
1722 t1d, p1d, dbz, rand1, kts, kte, i, j, &
1723 melti, vt_dbz_wt(i,:,j), &
1727 t1d, p1d, dbz, rand1, kts, kte, i, j, &
1731 refl_10cm(i,k,j) = max(-35., dbz(k))
1734 ENDIF diagflag_present
1736 IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0)
THEN
1738 re_qc1d(k) = re_qc_min
1739 re_qi1d(k) = re_qi_min
1740 re_qs1d(k) = re_qs_min
1743 call calc_effectrad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
1744 re_qc1d, re_qi1d, re_qs1d, lsml, kts, kte)
1746 re_cloud(i,k,j) = max(re_qc_min, min(re_qc1d(k), re_qc_max))
1747 re_ice(i,k,j) = max(re_qi_min, min(re_qi1d(k), re_qi_max))
1748 re_snow(i,k,j) = max(re_qs_min, min(re_qs1d(k), re_qs_max))
1751 ENDIF last_step_only
1767 do j = j_start, j_end
1769 do i = i_start, i_end
1770 pfils(i,k,j) = pfils(i,k,j)/dt_in
1771 pflls(i,k,j) = pflls(i,k,j)/dt_in
1780 deallocate_extended_diagnostics:
if (ext_diag)
then
1781 deallocate (prw_vcdc1)
1782 deallocate (prw_vcde1)
1783 deallocate (tpri_inu1)
1784 deallocate (tpri_ide1_d)
1785 deallocate (tpri_ide1_s)
1786 deallocate (tprs_ide1)
1787 deallocate (tprs_sde1_d)
1788 deallocate (tprs_sde1_s)
1789 deallocate (tprg_gde1_d)
1790 deallocate (tprg_gde1_s)
1791 deallocate (tpri_iha1)
1792 deallocate (tpri_wfz1)
1793 deallocate (tpri_rfz1)
1794 deallocate (tprg_rfz1)
1795 deallocate (tprs_scw1)
1796 deallocate (tprg_scw1)
1797 deallocate (tprg_rcs1)
1798 deallocate (tprs_rcs1)
1799 deallocate (tprr_rci1)
1800 deallocate (tprg_rcg1)
1801 deallocate (tprw_vcd1_c)
1802 deallocate (tprw_vcd1_e)
1803 deallocate (tprr_sml1)
1804 deallocate (tprr_gml1)
1805 deallocate (tprr_rcg1)
1806 deallocate (tprr_rcs1)
1807 deallocate (tprv_rev1)
1818 end if deallocate_extended_diagnostics
1885 nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dzq, &
1886 lsml, pptrain, pptsnow, pptgraul, pptice, &
1887#if ( WRF_CHEM == 1 )
1888 rainprod, evapprod, &
1890 rand1, rand2, rand3, &
1891 kts, kte, dt, ii, jj, &
1897 prw_vcdc1, prw_vcde1, &
1898 tpri_inu1, tpri_ide1_d, tpri_ide1_s, tprs_ide1, &
1899 tprs_sde1_d, tprs_sde1_s, &
1900 tprg_gde1_d, tprg_gde1_s, tpri_iha1, tpri_wfz1, &
1901 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1, &
1902 tprg_rcs1, tprs_rcs1, tprr_rci1, &
1903 tprg_rcg1, tprw_vcd1_c, &
1904 tprw_vcd1_e, tprr_sml1, tprr_gml1, tprr_rcg1, &
1905 tprr_rcs1, tprv_rev1, &
1906 tten1, qvten1, qrten1, qsten1, &
1907 qgten1, qiten1, niten1, nrten1, ncten1, qcten1, &
1909 fs_fac_rain, fs_fac_snow)
1916 integer,
intent(in):: kts, kte, ii, jj
1917 real(wp),
dimension(kts:kte),
intent(inout) :: &
1918 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
1919 nr1d, nc1d, nwfa1d, nifa1d, t1d
1920 real(wp),
dimension(kts:kte),
intent(out) :: pfil1, pfll1
1921 real(wp),
dimension(kts:kte),
intent(in) :: p1d, w1d, dzq
1922 real(wp),
intent(inout) :: pptrain, pptsnow, pptgraul, pptice
1923 real(wp),
intent(in) :: dt
1924 integer,
intent(in) :: lsml
1925 real(wp),
intent(in) :: rand1, rand2, rand3
1927 logical,
intent(in) :: ext_diag
1928 logical,
intent(in) :: sedi_semi
1929 integer,
intent(in) :: decfl
1930 real(wp),
dimension(:),
intent(out),
optional :: &
1933 prw_vcde1, tpri_inu1, tpri_ide1_d, &
1934 tpri_ide1_s, tprs_ide1, &
1935 tprs_sde1_d, tprs_sde1_s, tprg_gde1_d, &
1936 tprg_gde1_s, tpri_iha1, tpri_wfz1, &
1937 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1,&
1938 tprg_rcs1, tprs_rcs1, &
1939 tprr_rci1, tprg_rcg1, &
1940 tprw_vcd1_c, tprw_vcd1_e, tprr_sml1, &
1941 tprr_gml1, tprr_rcg1, &
1942 tprr_rcs1, tprv_rev1, tten1, qvten1, &
1943 qrten1, qsten1, qgten1, qiten1, niten1, &
1944 nrten1, ncten1, qcten1
1946 real(wp),
intent(in),
optional :: fs_fac_rain, fs_fac_snow
1948#if ( WRF_CHEM == 1 )
1949 real(wp),
dimension(kts:kte),
intent(inout) :: &
1954 real(wp),
dimension(kts:kte) :: tten, qvten, qcten, qiten, &
1955 qrten, qsten, qgten, niten, nrten, ncten, nwfaten, nifaten
1957 real(dp),
dimension(kts:kte) :: prw_vcd
1959 real(dp),
dimension(kts:kte) :: pnc_wcd, pnc_wau, pnc_rcw, &
1962 real(dp),
dimension(kts:kte) :: pna_rca, pna_sca, pna_gca, &
1963 pnd_rcd, pnd_scd, pnd_gcd
1965 real(dp),
dimension(kts:kte) :: prr_wau, prr_rcw, prr_rcs, &
1966 prr_rcg, prr_sml, prr_gml, &
1968 pnr_wau, pnr_rcs, pnr_rcg, &
1969 pnr_rci, pnr_sml, pnr_gml, &
1970 pnr_rev, pnr_rcr, pnr_rfz
1972 real(dp),
dimension(kts:kte) :: pri_inu, pni_inu, pri_ihm, &
1973 pni_ihm, pri_wfz, pni_wfz, &
1974 pri_rfz, pni_rfz, pri_ide, &
1975 pni_ide, pri_rci, pni_rci, &
1976 pni_sci, pni_iau, pri_iha, pni_iha
1978 real(dp),
dimension(kts:kte) :: prs_iau, prs_sci, prs_rcs, &
1979 prs_scw, prs_sde, prs_ihm, &
1982 real(dp),
dimension(kts:kte) :: prg_scw, prg_rfz, prg_gde, &
1983 prg_gcw, prg_rci, prg_rcs, &
1986 real(dp),
parameter:: zeroD0 = 0.0
1987 real(wp) :: dtcfl, rainsfc, graulsfc
1990 real(wp),
dimension(kts:kte) :: temp, pres, qv, pfll, pfil, pdummy
1991 real(wp),
dimension(kts:kte) :: rc, ri, rr, rs, rg, ni, nr, nc, nwfa, nifa
1992 real(wp),
dimension(kts:kte) :: rr_tmp, nr_tmp, rg_tmp
1993 real(wp),
dimension(kts:kte) :: rho, rhof, rhof2
1994 real(wp),
dimension(kts:kte) :: qvs, qvsi, delQvs
1995 real(wp),
dimension(kts:kte) :: satw, sati, ssatw, ssati
1996 real(wp),
dimension(kts:kte) :: diffu, visco, vsc2, &
1997 tcond, lvap, ocp, lvt2
1999 real(dp),
dimension(kts:kte) :: ilamr, ilamg, N0_r, N0_g
2000 real(wp),
dimension(kts:kte) :: mvd_r, mvd_c
2001 real(wp),
dimension(kts:kte) :: smob, smo2, smo1, smo0, &
2002 smoc, smod, smoe, smof
2004 real(wp),
dimension(kts:kte) :: sed_r, sed_s, sed_g, sed_i, sed_n,sed_c
2006 real(wp) :: rgvm, delta_tp, orho, lfus2, orhodt
2007 real(wp),
dimension(5):: onstep
2008 real(dp) :: N0_exp, N0_min, lam_exp, lamc, lamr, lamg
2009 real(dp) :: lami, ilami, ilamc
2010 real(wp) :: xDc, Dc_b, Dc_g, xDi, xDr, xDs, xDg, Ds_m, Dg_m
2011 real(dp) :: Dr_star, Dc_star
2012 real(wp) :: zeta1, zeta, taud, tau
2013 real(wp) :: stoke_r, stoke_s, stoke_g, stoke_i
2014 real(wp) :: vti, vtr, vts, vtg, vtc
2015 real(wp),
dimension(kts:kte+1):: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk, &
2017 real(wp),
dimension(kts:kte):: vts_boost
2018 real(wp) :: Mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, C_snow
2019 real(wp) :: a_, b_, loga_, A1, A2, tf
2020 real(wp) :: tempc, tc0, r_mvd1, r_mvd2, xkrat
2021 real(wp) :: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr
2022 real(wp) :: xsat, rate_max, sump, ratio
2023 real(wp) :: clap, fcd, dfcd
2024 real(wp) :: otemp, rvs, rvs_p, rvs_pp, gamsc, alphsc, t1_evap, t1_subl
2025 real(wp) :: r_frac, g_frac
2026 real(wp) :: Ef_rw, Ef_sw, Ef_gw, Ef_rr
2027 real(wp) :: Ef_ra, Ef_sa, Ef_ga
2028 real(wp) :: dtsave, odts, odt, odzq, hgt_agl, SR
2029 real(wp) :: xslw1, ygra1, zans1, eva_factor
2030 integer :: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq
2031 integer,
dimension(5) :: ksed1
2032 integer :: nir, nis, nig, nii, nic, niin
2033 integer :: idx_tc, idx_t, idx_s, idx_g1, idx_g, idx_r1, idx_r, &
2034 idx_i1, idx_i, idx_c, idx, idx_d, idx_n, idx_in
2037 logical,
dimension(kts:kte) :: L_qc, L_qi, L_qr, L_qs, L_qg
2038 logical :: debug_flag
2041 real(wp) :: fallspeed_adjustment_factor
2044 debug_flag = .false.
2047 write(*, *)
'DEBUG INFO, mp_thompson at (i,j) ', ii,
', ', jj
2159#if ( WRF_CHEM == 1 )
2229 qv(k) = max(1.e-10, qv1d(k))
2231 rho(k) = roverrv*pres(k) / (r*temp(k)*(qv(k)+roverrv))
2232 nwfa(k) = max(11.1e6*rho(k), min(9999.e6*rho(k), nwfa1d(k)*rho(k)))
2233 nifa(k) = max(nain1*0.01*rho(k), min(9999.e6*rho(k), nifa1d(k)*rho(k)))
2237 if (qc1d(k) .gt. r1)
then
2239 rc(k) = qc1d(k)*rho(k)
2240 nc(k) = max(2., min(nc1d(k)*rho(k), nt_c_max))
2242 if (nc(k).gt.10000.e6)
then
2244 elseif (nc(k).lt.100.)
then
2247 nu_c = nint(1000.e6/nc(k)) + 2
2248 nu_c = max(2, min(nu_c+nint(rand2), 15))
2250 lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
2251 xdc = (bm_r + nu_c + 1.) / lamc
2252 if (xdc.lt. d0c)
then
2253 lamc = cce(2,nu_c)/d0c
2254 elseif (xdc.gt. d0r*2.)
then
2255 lamc = cce(2,nu_c)/(d0r*2.)
2257 nc(k) = min(real(nt_c_max, kind=dp), ccg(1,nu_c)*ocg2(nu_c)*rc(k) &
2259 if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware))
then
2274 if (qi1d(k) .gt. r1)
then
2276 ri(k) = qi1d(k)*rho(k)
2277 ni(k) = max(r2, ni1d(k)*rho(k))
2278 if (ni(k).le. r2)
then
2280 ni(k) = min(nt_i_max, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
2283 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
2285 xdi = (bm_i + mu_i + 1.) * ilami
2286 if (xdi.lt. 5.e-6)
then
2288 ni(k) = min(nt_i_max, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
2289 elseif (xdi.gt. 300.e-6)
then
2290 lami = cie(2)/300.e-6
2291 ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i
2301 if (qr1d(k) .gt. r1)
then
2303 rr(k) = qr1d(k)*rho(k)
2304 nr(k) = max(r2, nr1d(k)*rho(k))
2305 if (nr(k).le. r2)
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
2311 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
2312 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
2313 if (mvd_r(k) .gt. 2.5e-3)
then
2315 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
2316 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
2317 elseif (mvd_r(k) .lt. d0r*0.75)
then
2319 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
2320 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
2329 if (qs1d(k) .gt. r1)
then
2331 rs(k) = qs1d(k)*rho(k)
2338 if (qg1d(k) .gt. r1)
then
2340 rg(k) = qg1d(k)*rho(k)
2365 tempc = temp(k) - 273.15
2366 rhof(k) = sqrt(rho_not/rho(k))
2367 rhof2(k) = sqrt(rhof(k))
2368 qvs(k) =
rslf(pres(k), temp(k))
2369 delqvs(k) = max(0.0,
rslf(pres(k), 273.15)-qv(k))
2370 if (tempc .le. 0.0)
then
2371 qvsi(k) =
rsif(pres(k), temp(k))
2375 satw(k) = qv(k)/qvs(k)
2376 sati(k) = qv(k)/qvsi(k)
2377 ssatw(k) = satw(k) - 1.
2378 ssati(k) = sati(k) - 1.
2379 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
2380 if (abs(ssati(k)).lt. eps) ssati(k) = 0.0
2381 if (no_micro .and. ssati(k).gt. 0.0) no_micro = .false.
2382 diffu(k) = 2.11e-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
2383 if (tempc .ge. 0.0)
then
2384 visco(k) = (1.718+0.0049*tempc)*1.0e-5
2386 visco(k) = (1.718+0.0049*tempc-1.2e-5*tempc*tempc)*1.0e-5
2388 ocp(k) = 1./(cp*(1.+0.887*qv(k)))
2389 vsc2(k) = sqrt(rho(k)/visco(k))
2390 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
2391 tcond(k) = (5.69 + 0.0168*tempc)*1.0e-5 * 418.936
2399 if (no_micro)
return
2404 if (.not. iiwarm)
then
2406 if (.not. l_qs(k)) cycle
2407 tc0 = min(-0.1, temp(k)-273.15)
2408 smob(k) = rs(k)*oams
2412 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3))
then
2415 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
2416 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
2417 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
2418 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
2419 + sa(10)*bm_s*bm_s*bm_s
2421 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
2422 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
2423 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
2424 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
2425 + sb(10)*bm_s*bm_s*bm_s
2426 smo2(k) = (smob(k)/a_)**(1./b_)
2430 loga_ = sa(1) + sa(2)*tc0 + sa(5)*tc0*tc0 + sa(9)*tc0*tc0*tc0
2432 b_ = sb(1) + sb(2)*tc0 + sb(5)*tc0*tc0 + sb(9)*tc0*tc0*tc0
2433 smo0(k) = a_ * smo2(k)**b_
2436 loga_ = sa(1) + sa(2)*tc0 + sa(3) &
2437 + sa(4)*tc0 + sa(5)*tc0*tc0 &
2438 + sa(6) + sa(7)*tc0*tc0 &
2439 + sa(8)*tc0 + sa(9)*tc0*tc0*tc0 &
2442 b_ = sb(1)+ sb(2)*tc0 + sb(3) + sb(4)*tc0 &
2443 + sb(5)*tc0*tc0 + sb(6) &
2444 + sb(7)*tc0*tc0 + sb(8)*tc0 &
2445 + sb(9)*tc0*tc0*tc0 + sb(10)
2446 smo1(k) = a_ * smo2(k)**b_
2449 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
2450 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
2451 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
2452 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
2453 + sa(10)*cse(1)*cse(1)*cse(1)
2455 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
2456 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
2457 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
2458 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
2459 smoc(k) = a_ * smo2(k)**b_
2462 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(13) &
2463 + sa(4)*tc0*cse(13) + sa(5)*tc0*tc0 &
2464 + sa(6)*cse(13)*cse(13) + sa(7)*tc0*tc0*cse(13) &
2465 + sa(8)*tc0*cse(13)*cse(13) + sa(9)*tc0*tc0*tc0 &
2466 + sa(10)*cse(13)*cse(13)*cse(13)
2468 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(13) + sb(4)*tc0*cse(13) &
2469 + sb(5)*tc0*tc0 + sb(6)*cse(13)*cse(13) &
2470 + sb(7)*tc0*tc0*cse(13) + sb(8)*tc0*cse(13)*cse(13) &
2471 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(13)*cse(13)*cse(13)
2472 smoe(k) = a_ * smo2(k)**b_
2475 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(16) &
2476 + sa(4)*tc0*cse(16) + sa(5)*tc0*tc0 &
2477 + sa(6)*cse(16)*cse(16) + sa(7)*tc0*tc0*cse(16) &
2478 + sa(8)*tc0*cse(16)*cse(16) + sa(9)*tc0*tc0*tc0 &
2479 + sa(10)*cse(16)*cse(16)*cse(16)
2481 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(16) + sb(4)*tc0*cse(16) &
2482 + sb(5)*tc0*tc0 + sb(6)*cse(16)*cse(16) &
2483 + sb(7)*tc0*tc0*cse(16) + sb(8)*tc0*cse(16)*cse(16) &
2484 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(16)*cse(16)*cse(16)
2485 smof(k) = a_ * smo2(k)**b_
2498 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
2500 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
2501 n0_r(k) = nr(k)*org2*lamr**cre(2)
2512 if (l_qr(k) .and. mvd_r(k).gt. d0r)
then
2513 ef_rr = max(-0.1, 1.0 - exp(2300.0*(mvd_r(k)-1950.0e-6)))
2514 pnr_rcr(k) = ef_rr * 2.0*nr(k)*rr(k)
2518 if (nc(k).gt.10000.e6)
then
2520 elseif (nc(k).lt.100.)
then
2523 nu_c = nint(1000.e6/nc(k)) + 2
2524 nu_c = max(2, min(nu_c+nint(rand2), 15))
2526 xdc = max(d0c*1.e6, ((rc(k)/(am_r*nc(k)))**obmr) * 1.e6)
2527 lamc = (nc(k)*am_r* ccg(2,nu_c) * ocg1(nu_c) / rc(k))**obmr
2528 mvd_c(k) = (3.0+nu_c+0.672) / lamc
2529 mvd_c(k) = max(d0c, min(mvd_c(k), d0r))
2534 if (rc(k).gt. 0.01e-3)
then
2535 dc_g = ((ccg(3,nu_c)*ocg2(nu_c))**obmr / lamc) * 1.e6
2536 dc_b = (xdc*xdc*xdc*dc_g*dc_g*dc_g - xdc*xdc*xdc*xdc*xdc*xdc) &
2538 zeta1 = 0.5*((6.25e-6*xdc*dc_b*dc_b*dc_b - 0.4) &
2539 + abs(6.25e-6*xdc*dc_b*dc_b*dc_b - 0.4))
2540 zeta = 0.027*rc(k)*zeta1
2541 taud = 0.5*((0.5*dc_b - 7.5) + abs(0.5*dc_b - 7.5)) + r1
2542 tau = 3.72/(rc(k)*taud)
2543 prr_wau(k) = zeta/tau
2544 prr_wau(k) = min(real(rc(k)*odts, kind=dp), prr_wau(k))
2545 pnr_wau(k) = prr_wau(k) / (am_r*nu_c*10.*d0r*d0r*d0r)
2546 pnc_wau(k) = min(real(nc(k)*odts, kind=dp), prr_wau(k) &
2547 / (am_r*mvd_c(k)*mvd_c(k)*mvd_c(k)))
2551 if (l_qr(k) .and. mvd_r(k).gt. d0r .and. mvd_c(k).gt. d0c)
then
2553 idx = 1 + int(nbr*log(real(mvd_r(k)/dr(1), kind=dp)) / log(dr(nbr)/dr(1)))
2555 ef_rw = t_efrw(idx, int(mvd_c(k)*1.e6))
2556 prr_rcw(k) = rhof(k)*t1_qr_qc*ef_rw*rc(k)*n0_r(k) &
2557 *((lamr+fv_r)**(-cre(9)))
2558 prr_rcw(k) = min(real(rc(k)*odts, kind=dp), prr_rcw(k))
2559 pnc_rcw(k) = rhof(k)*t1_qr_qc*ef_rw*nc(k)*n0_r(k) &
2560 *((lamr+fv_r)**(-cre(9)))
2561 pnc_rcw(k) = min(real(nc(k)*odts, kind=dp), pnc_rcw(k))
2565 if (l_qr(k) .and. mvd_r(k).gt. d0r)
then
2566 ef_ra =
eff_aero(mvd_r(k),0.04e-6,visco(k),rho(k),temp(k),
'r')
2568 pna_rca(k) = rhof(k)*t1_qr_qc*ef_ra*nwfa(k)*n0_r(k) &
2569 *((lamr+fv_r)**(-cre(9)))
2570 pna_rca(k) = min(real(nwfa(k)*odts, kind=dp), pna_rca(k))
2572 ef_ra =
eff_aero(mvd_r(k),0.8e-6,visco(k),rho(k),temp(k),
'r')
2573 pnd_rcd(k) = rhof(k)*t1_qr_qc*ef_ra*nifa(k)*n0_r(k) &
2574 *((lamr+fv_r)**(-cre(9)))
2575 pnd_rcd(k) = min(real(nifa(k)*odts, kind=dp), pnd_rcd(k))
2583 if (.not. iiwarm)
then
2587 if (l_qs(k)) xds = smoc(k) / smob(k)
2590 tempc = temp(k) - 273.15
2591 idx_tc = max(1, min(nint(-tempc), 45) )
2592 idx_t = int( (tempc-2.5)/5. ) - 1
2593 idx_t = max(1, -idx_t)
2594 idx_t = min(idx_t, ntb_t)
2595 it = max(1, min(nint(-tempc), 31) )
2598 if (rc(k).gt. r_c(1))
then
2599 nic = nint(log10(rc(k)))
2600 do_loop_rc:
do nn = nic-1, nic+1
2602 if ( (rc(k)/10.**nn).ge.1.0 .and. (rc(k)/10.**nn).lt.10.0 )
exit do_loop_rc
2604 idx_c = int(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
2605 idx_c = max(1, min(idx_c, ntb_c))
2611 idx_n = nint(1.0 + real(nbc, kind=wp) * log(real(nc(k)/t_nc(1), kind=dp)) / nic1)
2612 idx_n = max(1, min(idx_n, nbc))
2615 if (ri(k).gt. r_i(1))
then
2616 nii = nint(log10(ri(k)))
2617 do_loop_ri:
do nn = nii-1, nii+1
2619 if ( (ri(k)/10.**nn).ge.1.0 .and. (ri(k)/10.**nn).lt.10.0 )
exit do_loop_ri
2621 idx_i = int(ri(k)/10.**n) + 10*(n-nii2) - (n-nii2)
2622 idx_i = max(1, min(idx_i, ntb_i))
2627 if (ni(k).gt. nt_i(1))
then
2628 nii = nint(log10(ni(k)))
2629 do_loop_ni:
do nn = nii-1, nii+1
2631 if ( (ni(k)/10.**nn).ge.1.0 .and. (ni(k)/10.**nn).lt.10.0 )
exit do_loop_ni
2633 idx_i1 = int(ni(k)/10.**n) + 10*(n-nii3) - (n-nii3)
2634 idx_i1 = max(1, min(idx_i1, ntb_i1))
2640 if (rr(k).gt. r_r(1))
then
2641 nir = nint(log10(rr(k)))
2642 do_loop_rr:
do nn = nir-1, nir+1
2644 if ( (rr(k)/10.**nn).ge.1.0 .and. (rr(k)/10.**nn).lt.10.0 )
exit do_loop_rr
2646 idx_r = int(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2)
2647 idx_r = max(1, min(idx_r, ntb_r))
2650 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
2651 n0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
2652 nir = nint(log10(n0_exp))
2653 do_loop_nr:
do nn = nir-1, nir+1
2655 if ( (n0_exp/10.**nn).ge.1.0 .and. (n0_exp/10.**nn).lt.10.0 )
exit do_loop_nr
2657 idx_r1 = int(n0_exp/10.**n) + 10*(n-nir3) - (n-nir3)
2658 idx_r1 = max(1, min(idx_r1, ntb_r1))
2665 if (rs(k).gt. r_s(1))
then
2666 nis = nint(log10(rs(k)))
2667 do_loop_rs:
do nn = nis-1, nis+1
2669 if ( (rs(k)/10.**nn).ge.1.0 .and. (rs(k)/10.**nn).lt.10.0 )
exit do_loop_rs
2671 idx_s = int(rs(k)/10.**n) + 10*(n-nis2) - (n-nis2)
2672 idx_s = max(1, min(idx_s, ntb_s))
2678 if (rg(k).gt. r_g(1))
then
2679 nig = nint(log10(rg(k)))
2680 do_loop_rg:
do nn = nig-1, nig+1
2682 if ( (rg(k)/10.**nn).ge.1.0 .and. (rg(k)/10.**nn).lt.10.0 )
exit do_loop_rg
2684 idx_g = int(rg(k)/10.**n) + 10*(n-nig2) - (n-nig2)
2685 idx_g = max(1, min(idx_g, ntb_g))
2688 lam_exp = lamg * (cgg(3)*ogg2*ogg1)**bm_g
2689 n0_exp = ogg1*rg(k)/am_g * lam_exp**cge(1)
2690 nig = nint(log10(real(n0_exp, kind=dp)))
2691 do_loop_ng:
do nn = nig-1, nig+1
2693 if ( (n0_exp/10.**nn).ge.1.0 .and. (n0_exp/10.**nn).lt.10.0 )
exit do_loop_ng
2695 idx_g1 = int(n0_exp/10.**n) + 10*(n-nig3) - (n-nig3)
2696 idx_g1 = max(1, min(idx_g1, ntb_g1))
2704 rvs = rho(k)*qvsi(k)
2705 rvs_p = rvs*otemp*(lsub*otemp*orv - 1.)
2706 rvs_pp = rvs * ( otemp*(lsub*otemp*orv - 1.) &
2707 *otemp*(lsub*otemp*orv - 1.) &
2708 + (-2.*lsub*otemp*otemp*otemp*orv) &
2710 gamsc = lsub*diffu(k)/tcond(k) * rvs_p
2711 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
2712 * rvs_pp/rvs_p * rvs/rvs_p
2713 alphsc = max(1.e-9, alphsc)
2715 if (abs(xsat).lt. 1.e-9) xsat=0.
2716 t1_subl = 4.*pi*( 1.0 - alphsc*xsat &
2717 + 2.*alphsc*alphsc*xsat*xsat &
2718 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
2722 if (l_qc(k) .and. mvd_c(k).gt. d0c)
then
2723 if (xds .gt. d0s)
then
2724 idx = 1 + int(nbs*log(real(xds/ds(1), kind=dp)) / log(ds(nbs)/ds(1)))
2726 ef_sw = t_efsw(idx, int(mvd_c(k)*1.e6))
2727 prs_scw(k) = rhof(k)*t1_qs_qc*ef_sw*rc(k)*smoe(k)
2728 prs_scw(k) = min(real(rc(k)*odts, kind=dp), prs_scw(k))
2729 pnc_scw(k) = rhof(k)*t1_qs_qc*ef_sw*nc(k)*smoe(k)
2730 pnc_scw(k) = min(real(nc(k)*odts, kind=dp), pnc_scw(k))
2734 if (rg(k).ge. r_g(1) .and. mvd_c(k).gt. d0c)
then
2735 xdg = (bm_g + mu_g + 1.) * ilamg(k)
2736 vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
2737 stoke_g = mvd_c(k)*mvd_c(k)*vtg*rho_w/(9.*visco(k)*xdg)
2738 if (xdg.gt. d0g)
then
2739 if (stoke_g.ge.0.4 .and. stoke_g.le.10.)
then
2740 ef_gw = 0.55*log10(2.51*stoke_g)
2741 elseif (stoke_g.lt.0.4)
then
2743 elseif (stoke_g.gt.10)
then
2746 prg_gcw(k) = rhof(k)*t1_qg_qc*ef_gw*rc(k)*n0_g(k) &
2748 pnc_gcw(k) = rhof(k)*t1_qg_qc*ef_gw*nc(k)*n0_g(k) &
2750 pnc_gcw(k) = min(real(nc(k)*odts, kind=dp), pnc_gcw(k))
2756 if (rs(k) .gt. r_s(1))
then
2757 ef_sa =
eff_aero(xds,0.04e-6,visco(k),rho(k),temp(k),
's')
2758 pna_sca(k) = rhof(k)*t1_qs_qc*ef_sa*nwfa(k)*smoe(k)
2759 pna_sca(k) = min(real(nwfa(k)*odts, kind=dp), pna_sca(k))
2761 ef_sa =
eff_aero(xds,0.8e-6,visco(k),rho(k),temp(k),
's')
2762 pnd_scd(k) = rhof(k)*t1_qs_qc*ef_sa*nifa(k)*smoe(k)
2763 pnd_scd(k) = min(real(nifa(k)*odts, kind=dp), pnd_scd(k))
2765 if (rg(k) .gt. r_g(1))
then
2766 xdg = (bm_g + mu_g + 1.) * ilamg(k)
2767 ef_ga =
eff_aero(xdg,0.04e-6,visco(k),rho(k),temp(k),
'g')
2768 pna_gca(k) = rhof(k)*t1_qg_qc*ef_ga*nwfa(k)*n0_g(k) &
2770 pna_gca(k) = min(real(nwfa(k)*odts, kind=dp), pna_gca(k))
2772 ef_ga =
eff_aero(xdg,0.8e-6,visco(k),rho(k),temp(k),
'g')
2773 pnd_gcd(k) = rhof(k)*t1_qg_qc*ef_ga*nifa(k)*n0_g(k) &
2775 pnd_gcd(k) = min(real(nifa(k)*odts, kind=dp), pnd_gcd(k))
2781 if (rr(k).ge. r_r(1))
then
2782 if (rs(k).ge. r_s(1))
then
2783 if (temp(k).lt.t_0)
then
2784 prr_rcs(k) = -(tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2785 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
2786 + tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
2787 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r))
2788 prs_rcs(k) = tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2789 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
2790 - tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2791 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
2792 prg_rcs(k) = tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
2793 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2794 + tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2795 + tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
2796 prr_rcs(k) = max(real(-rr(k)*odts, kind=dp), prr_rcs(k))
2797 prs_rcs(k) = max(real(-rs(k)*odts, kind=dp), prs_rcs(k))
2798 prg_rcs(k) = min(real((rr(k)+rs(k))*odts, kind=dp), prg_rcs(k))
2799 pnr_rcs(k) = tnr_racs1(idx_s,idx_t,idx_r1,idx_r) &
2800 + tnr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2801 + tnr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2802 + tnr_sacr2(idx_s,idx_t,idx_r1,idx_r)
2803 pnr_rcs(k) = min(real(nr(k)*odts, kind=dp), pnr_rcs(k))
2805 prs_rcs(k) = -tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2806 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2807 + tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2808 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r)
2809 prs_rcs(k) = max(real(-rs(k)*odts, kind=dp), prs_rcs(k))
2810 prr_rcs(k) = -prs_rcs(k)
2817 if (rg(k).ge. r_g(1))
then
2818 if (temp(k).lt.t_0)
then
2819 prg_rcg(k) = tmr_racg(idx_g1,idx_g,idx_r1,idx_r) &
2820 + tcr_gacr(idx_g1,idx_g,idx_r1,idx_r)
2821 prg_rcg(k) = min(real(rr(k)*odts, kind=dp), prg_rcg(k))
2822 prr_rcg(k) = -prg_rcg(k)
2823 pnr_rcg(k) = tnr_racg(idx_g1,idx_g,idx_r1,idx_r) &
2824 + tnr_gacr(idx_g1,idx_g,idx_r1,idx_r)
2825 pnr_rcg(k) = min(real(nr(k)*odts, kind=dp), pnr_rcg(k))
2827 prr_rcg(k) = tcg_racg(idx_g1,idx_g,idx_r1,idx_r)
2828 prr_rcg(k) = min(real(rg(k)*odts, kind=dp), prr_rcg(k))
2829 prg_rcg(k) = -prr_rcg(k)
2831 pnr_rcg(k) = -1.5*tnr_gacr(idx_g1,idx_g,idx_r1,idx_r)
2836 if (temp(k).lt.t_0)
then
2837 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
2841 c_snow = c_sqrd + (tempc+1.5)*(c_cube-c_sqrd)/(-30.+1.5)
2842 c_snow = max(c_sqrd, min(c_snow, c_cube))
2843 prs_sde(k) = c_snow*t1_subl*diffu(k)*ssati(k)*rvs &
2844 * (t1_qs_sd*smo1(k) &
2845 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
2846 if (prs_sde(k).lt. 0.)
then
2847 prs_sde(k) = max(real(-rs(k)*odts, kind=dp), prs_sde(k), real(rate_max, kind=dp))
2849 prs_sde(k) = min(prs_sde(k), real(rate_max, kind=dp))
2853 if (l_qg(k) .and. ssati(k).lt. -eps)
then
2854 prg_gde(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
2855 * n0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
2856 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
2857 if (prg_gde(k).lt. 0.)
then
2858 prg_gde(k) = max(real(-rg(k)*odts, kind=dp), prg_gde(k), real(rate_max, kind=dp))
2860 prg_gde(k) = min(prg_gde(k), real(rate_max, kind=dp))
2868 if (prs_scw(k).gt.5.0*prs_sde(k) .and. &
2869 prs_sde(k).gt.eps)
then
2870 r_frac = min(30.0_dp, prs_scw(k)/prs_sde(k))
2871 g_frac = min(0.75, 0.15 + (r_frac-5.)*.028)
2872 vts_boost(k) = min(1.5, 1.1 + (r_frac-5.)*.016)
2873 prg_scw(k) = g_frac*prs_scw(k)
2874 prs_scw(k) = (1. - g_frac)*prs_scw(k)
2882 if (temp(k).lt.t_0)
then
2884 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
2900 if (dustyice .AND. (is_aerosol_aware .or. merra2_aerosol_aware))
then
2901 xni = icedemott(tempc,qvs(k),qvs(k),qvsi(k),rho(k),nifa(k))
2907 if (xni.gt. nt_in(1))
then
2908 niin = nint(log10(xni))
2909 do_loop_xni:
do nn = niin-1, niin+1
2911 if ( (xni/10.**nn).ge.1.0 .and. (xni/10.**nn).lt.10.0 )
exit do_loop_xni
2913 idx_in = int(xni/10.**n) + 10*(n-niin2) - (n-niin2)
2914 idx_in = max(1, min(idx_in, ntb_in))
2920 if (rr(k).gt. r_r(1))
then
2921 prg_rfz(k) = tpg_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
2922 pri_rfz(k) = tpi_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
2923 pni_rfz(k) = tni_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
2924 pnr_rfz(k) = tnr_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
2925 pnr_rfz(k) = min(real(nr(k)*odts, kind=dp), pnr_rfz(k))
2926 elseif (rr(k).gt. r1 .and. temp(k).lt.hgfr)
then
2927 pri_rfz(k) = rr(k)*odts
2928 pni_rfz(k) = pnr_rfz(k)
2931 if (rc(k).gt. r_c(1))
then
2932 pri_wfz(k) = tpi_qcfz(idx_c,idx_n,idx_tc,idx_in)*odts
2933 pri_wfz(k) = min(real(rc(k)*odts, kind=dp), pri_wfz(k))
2934 pni_wfz(k) = tni_qcfz(idx_c,idx_n,idx_tc,idx_in)*odts
2935 pni_wfz(k) = min(real(nc(k)*odts, kind=dp), pri_wfz(k)/(2.0_dp*xm0i), &
2937 elseif (rc(k).gt. r1 .and. temp(k).lt.hgfr)
then
2938 pri_wfz(k) = rc(k)*odts
2939 pni_wfz(k) = nc(k)*odts
2944 if ( (ssati(k).ge. ssati_min) .or. (ssatw(k).gt. eps &
2945 .and. temp(k).lt.253.15) )
then
2946 if (dustyice .AND. (is_aerosol_aware .or. merra2_aerosol_aware))
then
2947 xnc = icedemott(tempc,qv(k),qvs(k),qvsi(k),rho(k),nifa(k))
2948 xnc = xnc*(1.0 + 50.*rand3)
2950 xnc = min(xnc_max, tno*exp(ato*(t_0-temp(k))))
2952 xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dtsave
2953 pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts
2954 pri_inu(k) = min(real(rate_max, kind=dp), xm0i*pni_inu(k))
2955 pni_inu(k) = pri_inu(k)/xm0i
2959 xni = smo0(k)+ni(k) + (pni_rfz(k)+pni_wfz(k)+pni_inu(k))*dtsave
2960 if ((is_aerosol_aware .or. merra2_aerosol_aware) .AND. homogice .AND. (xni.le.nt_i_max) &
2961 .AND.(temp(k).lt.238).AND.(ssati(k).ge.0.4) )
then
2962 xnc =
icekoop(temp(k),qv(k),qvs(k),nwfa(k), dtsave)
2963 pni_iha(k) = xnc*odts
2964 pri_iha(k) = min(real(rate_max, kind=dp), xm0i*0.1*pni_iha(k))
2965 pni_iha(k) = pri_iha(k)/(xm0i*0.1)
2972 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
2974 xdi = max(real(d0i, kind=dp), (bm_i + mu_i + 1.) * ilami)
2975 xmi = am_i*xdi**bm_i
2977 pri_ide(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
2978 *oig1*cig(5)*ni(k)*ilami
2980 if (pri_ide(k) .lt. 0.0)
then
2981 pri_ide(k) = max(real(-ri(k)*odts, kind=dp), pri_ide(k), real(rate_max, kind=dp))
2982 pni_ide(k) = pri_ide(k)*oxmi
2983 pni_ide(k) = max(real(-ni(k)*odts, kind=dp), pni_ide(k))
2985 pri_ide(k) = min(pri_ide(k), real(rate_max, kind=dp))
2986 prs_ide(k) = (1.0_dp-tpi_ide(idx_i,idx_i1))*pri_ide(k)
2987 pri_ide(k) = tpi_ide(idx_i,idx_i1)*pri_ide(k)
2992 if ( (idx_i.eq. ntb_i) .or. (xdi.gt. 5.0*d0s) )
then
2993 prs_iau(k) = ri(k)*.99*odts
2994 pni_iau(k) = ni(k)*.95*odts
2995 elseif (xdi.lt. 0.1*d0s)
then
2999 prs_iau(k) = tps_iaus(idx_i,idx_i1)*odts
3000 prs_iau(k) = min(real(ri(k)*.99*odts, kind=dp), prs_iau(k))
3001 pni_iau(k) = tni_iaus(idx_i,idx_i1)*odts
3002 pni_iau(k) = min(real(ni(k)*.95*odts, kind=dp), pni_iau(k))
3008 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
3010 xdi = max(real(d0i, kind=dp), (bm_i + mu_i + 1.) * ilami)
3011 xmi = am_i*xdi**bm_i
3013 if (rs(k).ge. r_s(1))
then
3014 prs_sci(k) = t1_qs_qi*rhof(k)*ef_si*ri(k)*smoe(k)
3015 pni_sci(k) = prs_sci(k) * oxmi
3019 if (rr(k).ge. r_r(1) .and. mvd_r(k).gt. 4.*xdi)
then
3021 pri_rci(k) = rhof(k)*t1_qr_qi*ef_ri*ri(k)*n0_r(k) &
3022 *((lamr+fv_r)**(-cre(9)))
3023 pnr_rci(k) = rhof(k)*t1_qr_qi*ef_ri*ni(k)*n0_r(k) &
3024 *((lamr+fv_r)**(-cre(9)))
3025 pni_rci(k) = pri_rci(k) * oxmi
3026 prr_rci(k) = rhof(k)*t2_qr_qi*ef_ri*ni(k)*n0_r(k) &
3027 *((lamr+fv_r)**(-cre(8)))
3028 prr_rci(k) = min(real(rr(k)*odts, kind=dp), prr_rci(k))
3029 prg_rci(k) = pri_rci(k) + prr_rci(k)
3034 if (prg_gcw(k).gt. eps .and. tempc.gt.-8.0)
then
3036 if (tempc.ge.-5.0 .and. tempc.lt.-3.0)
then
3037 tf = 0.5*(-3.0 - tempc)
3038 elseif (tempc.gt.-8.0 .and. tempc.lt.-5.0)
then
3039 tf = 0.33333333*(8.0 + tempc)
3041 pni_ihm(k) = 3.5e8*tf*prg_gcw(k)
3042 pri_ihm(k) = xm0i*pni_ihm(k)
3043 prs_ihm(k) = prs_scw(k)/(prs_scw(k)+prg_gcw(k)) &
3045 prg_ihm(k) = prg_gcw(k)/(prs_scw(k)+prg_gcw(k)) &
3054 prr_sml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delqvs(k)) &
3055 * (t1_qs_me*smo1(k) + t2_qs_me*rhof2(k)*vsc2(k)*smof(k))
3056 if (prr_sml(k) .gt. 0.)
then
3057 prr_sml(k) = prr_sml(k) + 4218.*olfus*tempc &
3058 * (prr_rcs(k)+prs_scw(k))
3059 prr_sml(k) = min(real(rs(k)*odts, kind=dp), prr_sml(k))
3060 pnr_sml(k) = smo0(k)/rs(k)*prr_sml(k) * 10.0**(-0.25*tempc)
3061 pnr_sml(k) = min(real(smo0(k)*odts, kind=dp), pnr_sml(k))
3062 elseif (ssati(k).lt. 0.)
then
3064 prs_sde(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
3065 * (t1_qs_sd*smo1(k) &
3066 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
3067 prs_sde(k) = max(real(-rs(k)*odts, kind=dp), prs_sde(k))
3072 prr_gml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delqvs(k)) &
3073 * n0_g(k)*(t1_qg_me*ilamg(k)**cge(10) &
3074 + t2_qg_me*rhof2(k)*vsc2(k)*ilamg(k)**cge(11))
3075 if (prr_gml(k) .gt. 0.)
then
3076 prr_gml(k) = min(real(rg(k)*odts, kind=dp), prr_gml(k))
3077 pnr_gml(k) = n0_g(k)*cgg(2)*ilamg(k)**cge(2) / rg(k) &
3078 * prr_gml(k) * 10.0**(-0.5*tempc)
3079 elseif (ssati(k).lt. 0.)
then
3081 prg_gde(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
3082 * n0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
3083 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
3084 prg_gde(k) = max(real(-rg(k)*odts, kind=dp), prg_gde(k))
3092 if (dt .gt. 120.)
then
3093 prr_rcw(k)=prr_rcw(k)+prs_scw(k)+prg_gcw(k)
3111 sump = pri_inu(k) + pri_ide(k) + prs_ide(k) &
3112 + prs_sde(k) + prg_gde(k) + pri_iha(k)
3113 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
3114 if ( (sump.gt. eps .and. sump.gt. rate_max) .or. &
3115 (sump.lt. -eps .and. sump.lt. rate_max) )
then
3116 ratio = rate_max/sump
3117 pri_inu(k) = pri_inu(k) * ratio
3118 pri_ide(k) = pri_ide(k) * ratio
3119 pni_ide(k) = pni_ide(k) * ratio
3120 prs_ide(k) = prs_ide(k) * ratio
3121 prs_sde(k) = prs_sde(k) * ratio
3122 prg_gde(k) = prg_gde(k) * ratio
3123 pri_iha(k) = pri_iha(k) * ratio
3127 sump = -prr_wau(k) - pri_wfz(k) - prr_rcw(k) &
3128 - prs_scw(k) - prg_scw(k) - prg_gcw(k)
3129 rate_max = -rc(k)*odts
3130 if (sump.lt. rate_max .and. l_qc(k))
then
3131 ratio = rate_max/sump
3132 prr_wau(k) = prr_wau(k) * ratio
3133 pri_wfz(k) = pri_wfz(k) * ratio
3134 prr_rcw(k) = prr_rcw(k) * ratio
3135 prs_scw(k) = prs_scw(k) * ratio
3136 prg_scw(k) = prg_scw(k) * ratio
3137 prg_gcw(k) = prg_gcw(k) * ratio
3141 sump = pri_ide(k) - prs_iau(k) - prs_sci(k) &
3143 rate_max = -ri(k)*odts
3144 if (sump.lt. rate_max .and. l_qi(k))
then
3145 ratio = rate_max/sump
3146 pri_ide(k) = pri_ide(k) * ratio
3147 prs_iau(k) = prs_iau(k) * ratio
3148 prs_sci(k) = prs_sci(k) * ratio
3149 pri_rci(k) = pri_rci(k) * ratio
3153 sump = -prg_rfz(k) - pri_rfz(k) - prr_rci(k) &
3154 + prr_rcs(k) + prr_rcg(k)
3155 rate_max = -rr(k)*odts
3156 if (sump.lt. rate_max .and. l_qr(k))
then
3157 ratio = rate_max/sump
3158 prg_rfz(k) = prg_rfz(k) * ratio
3159 pri_rfz(k) = pri_rfz(k) * ratio
3160 prr_rci(k) = prr_rci(k) * ratio
3161 prr_rcs(k) = prr_rcs(k) * ratio
3162 prr_rcg(k) = prr_rcg(k) * ratio
3166 sump = prs_sde(k) - prs_ihm(k) - prr_sml(k) &
3168 rate_max = -rs(k)*odts
3169 if (sump.lt. rate_max .and. l_qs(k))
then
3170 ratio = rate_max/sump
3171 prs_sde(k) = prs_sde(k) * ratio
3172 prs_ihm(k) = prs_ihm(k) * ratio
3173 prr_sml(k) = prr_sml(k) * ratio
3174 prs_rcs(k) = prs_rcs(k) * ratio
3178 sump = prg_gde(k) - prg_ihm(k) - prr_gml(k) &
3180 rate_max = -rg(k)*odts
3181 if (sump.lt. rate_max .and. l_qg(k))
then
3182 ratio = rate_max/sump
3183 prg_gde(k) = prg_gde(k) * ratio
3184 prg_ihm(k) = prg_ihm(k) * ratio
3185 prr_gml(k) = prr_gml(k) * ratio
3186 prg_rcg(k) = prg_rcg(k) * ratio
3191 pri_ihm(k) = prs_ihm(k) + prg_ihm(k)
3192 ratio = min( abs(prr_rcg(k)), abs(prg_rcg(k)) )
3193 prr_rcg(k) = ratio * sign(1.0, sngl(prr_rcg(k)))
3194 prg_rcg(k) = -prr_rcg(k)
3195 if (temp(k).gt.t_0)
then
3196 ratio = min( abs(prr_rcs(k)), abs(prs_rcs(k)) )
3197 prr_rcs(k) = ratio * sign(1.0, sngl(prr_rcs(k)))
3198 prs_rcs(k) = -prr_rcs(k)
3209 lfus2 = lsub - lvap(k)
3212 if (is_aerosol_aware)
then
3213 nwfaten(k) = nwfaten(k) - (pna_rca(k) + pna_sca(k) &
3214 + pna_gca(k) + pni_iha(k)) * orho
3215 nifaten(k) = nifaten(k) - (pnd_rcd(k) + pnd_scd(k) &
3216 + pnd_gcd(k)) * orho
3218 nifaten(k) = nifaten(k) - pni_inu(k)*orho
3225 qvten(k) = qvten(k) + (-pri_inu(k) - pri_iha(k) - pri_ide(k) &
3226 - prs_ide(k) - prs_sde(k) - prg_gde(k)) &
3230 qcten(k) = qcten(k) + (-prr_wau(k) - pri_wfz(k) &
3231 - prr_rcw(k) - prs_scw(k) - prg_scw(k) &
3236 ncten(k) = ncten(k) + (-pnc_wau(k) - pnc_rcw(k) &
3237 - pni_wfz(k) - pnc_scw(k) - pnc_gcw(k)) &
3242 xrc=max(r1, (qc1d(k) + qcten(k)*dtsave)*rho(k))
3243 xnc=max(2., (nc1d(k) + ncten(k)*dtsave)*rho(k))
3244 if (xrc .gt. r1)
then
3245 if (xnc.gt.10000.e6)
then
3247 elseif (xnc.lt.100.)
then
3250 nu_c = nint(1000.e6/xnc) + 2
3251 nu_c = max(2, min(nu_c+nint(rand2), 15))
3253 lamc = (xnc*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
3254 xdc = (bm_r + nu_c + 1.) / lamc
3255 if (xdc.lt. d0c)
then
3256 lamc = cce(2,nu_c)/d0c
3257 xnc = ccg(1,nu_c)*ocg2(nu_c)*xrc/am_r*lamc**bm_r
3258 ncten(k) = (xnc-nc1d(k)*rho(k))*odts*orho
3259 elseif (xdc.gt. d0r*2.)
then
3260 lamc = cce(2,nu_c)/(d0r*2.)
3261 xnc = ccg(1,nu_c)*ocg2(nu_c)*xrc/am_r*lamc**bm_r
3262 ncten(k) = (xnc-nc1d(k)*rho(k))*odts*orho
3265 ncten(k) = -nc1d(k)*odts
3267 xnc=max(0.,(nc1d(k) + ncten(k)*dtsave)*rho(k))
3268 if (xnc.gt.nt_c_max) &
3269 ncten(k) = (nt_c_max-nc1d(k)*rho(k))*odts*orho
3272 qiten(k) = qiten(k) + (pri_inu(k) + pri_iha(k) + pri_ihm(k) &
3273 + pri_wfz(k) + pri_rfz(k) + pri_ide(k) &
3274 - prs_iau(k) - prs_sci(k) - pri_rci(k)) &
3278 niten(k) = niten(k) + (pni_inu(k) + pni_iha(k) + pni_ihm(k) &
3279 + pni_wfz(k) + pni_rfz(k) + pni_ide(k) &
3280 - pni_iau(k) - pni_sci(k) - pni_rci(k)) &
3285 xri=max(r1,(qi1d(k) + qiten(k)*dtsave)*rho(k))
3286 xni=max(r2,(ni1d(k) + niten(k)*dtsave)*rho(k))
3287 if (xri.gt. r1)
then
3288 lami = (am_i*cig(2)*oig1*xni/xri)**obmi
3290 xdi = (bm_i + mu_i + 1.) * ilami
3291 if (xdi.lt. 5.e-6)
then
3293 xni = min(nt_i_max, cig(1)*oig2*xri/am_i*lami**bm_i)
3294 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
3295 elseif (xdi.gt. 300.e-6)
then
3296 lami = cie(2)/300.e-6
3297 xni = cig(1)*oig2*xri/am_i*lami**bm_i
3298 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
3301 niten(k) = -ni1d(k)*odts
3303 xni=max(0.,(ni1d(k) + niten(k)*dtsave)*rho(k))
3304 if (xni.gt.nt_i_max) &
3305 niten(k) = (nt_i_max-ni1d(k)*rho(k))*odts*orho
3308 qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) &
3309 + prr_sml(k) + prr_gml(k) + prr_rcs(k) &
3310 + prr_rcg(k) - prg_rfz(k) &
3311 - pri_rfz(k) - prr_rci(k)) &
3315 nrten(k) = nrten(k) + (pnr_wau(k) + pnr_sml(k) + pnr_gml(k) &
3316 - (pnr_rfz(k) + pnr_rcr(k) + pnr_rcg(k) &
3317 + pnr_rcs(k) + pnr_rci(k) + pni_rfz(k)) ) &
3322 xrr=max(r1,(qr1d(k) + qrten(k)*dtsave)*rho(k))
3323 xnr=max(r2,(nr1d(k) + nrten(k)*dtsave)*rho(k))
3324 if (xrr.gt. r1)
then
3325 lamr = (am_r*crg(3)*org2*xnr/xrr)**obmr
3326 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3327 if (mvd_r(k) .gt. 2.5e-3)
then
3329 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3330 xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
3331 nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
3332 elseif (mvd_r(k) .lt. d0r*0.75)
then
3334 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3335 xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
3336 nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
3339 qrten(k) = -qr1d(k)*odts
3340 nrten(k) = -nr1d(k)*odts
3344 qsten(k) = qsten(k) + (prs_iau(k) + prs_sde(k) &
3345 + prs_sci(k) + prs_scw(k) + prs_rcs(k) &
3346 + prs_ide(k) - prs_ihm(k) - prr_sml(k)) &
3350 qgten(k) = qgten(k) + (prg_scw(k) + prg_rfz(k) &
3351 + prg_gde(k) + prg_rcg(k) + prg_gcw(k) &
3352 + prg_rci(k) + prg_rcs(k) - prg_ihm(k) &
3357 if (temp(k).lt.t_0)
then
3359 + ( lsub*ocp(k)*(pri_inu(k) + pri_ide(k) &
3360 + prs_ide(k) + prs_sde(k) &
3361 + prg_gde(k) + pri_iha(k)) &
3362 + lfus2*ocp(k)*(pri_wfz(k) + pri_rfz(k) &
3363 + prg_rfz(k) + prs_scw(k) &
3364 + prg_scw(k) + prg_gcw(k) &
3365 + prg_rcs(k) + prs_rcs(k) &
3366 + prr_rci(k) + prg_rcg(k)) &
3370 + ( lfus*ocp(k)*(-prr_sml(k) - prr_gml(k) &
3371 - prr_rcg(k) - prr_rcs(k)) &
3372 + lsub*ocp(k)*(prs_sde(k) + prg_gde(k)) &
3382 temp(k) = t1d(k) + dt*tten(k)
3384 tempc = temp(k) - 273.15
3385 qv(k) = max(1.e-10, qv1d(k) + dt*qvten(k))
3386 rho(k) = roverrv*pres(k) / (r*temp(k)*(qv(k)+roverrv))
3387 rhof(k) = sqrt(rho_not/rho(k))
3388 rhof2(k) = sqrt(rhof(k))
3389 qvs(k) =
rslf(pres(k), temp(k))
3390 ssatw(k) = qv(k)/qvs(k) - 1.
3391 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
3392 diffu(k) = 2.11e-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
3393 if (tempc .ge. 0.0)
then
3394 visco(k) = (1.718+0.0049*tempc)*1.0e-5
3396 visco(k) = (1.718+0.0049*tempc-1.2e-5*tempc*tempc)*1.0e-5
3398 vsc2(k) = sqrt(rho(k)/visco(k))
3399 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
3400 tcond(k) = (5.69 + 0.0168*tempc)*1.0e-5 * 418.936
3401 ocp(k) = 1./(cp*(1.+0.887*qv(k)))
3402 lvt2(k)=lvap(k)*lvap(k)*ocp(k)*orv*otemp*otemp
3403 if (is_aerosol_aware) &
3404 nwfa(k) = max(11.1e6*rho(k), (nwfa1d(k) + nwfaten(k)*dt)*rho(k))
3408 if ((qc1d(k) + qcten(k)*dt) .gt. r1)
then
3409 rc(k) = (qc1d(k) + qcten(k)*dt)*rho(k)
3410 nc(k) = max(2., min((nc1d(k)+ncten(k)*dt)*rho(k), nt_c_max))
3411 if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware))
then
3425 if ((qi1d(k) + qiten(k)*dt) .gt. r1)
then
3426 ri(k) = (qi1d(k) + qiten(k)*dt)*rho(k)
3427 ni(k) = max(r2, (ni1d(k) + niten(k)*dt)*rho(k))
3435 if ((qr1d(k) + qrten(k)*dt) .gt. r1)
then
3436 rr(k) = (qr1d(k) + qrten(k)*dt)*rho(k)
3437 nr(k) = max(r2, (nr1d(k) + nrten(k)*dt)*rho(k))
3439 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3440 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3441 if (mvd_r(k) .gt. 2.5e-3)
then
3443 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3444 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
3445 elseif (mvd_r(k) .lt. d0r*0.75)
then
3447 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3448 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
3456 if ((qs1d(k) + qsten(k)*dt) .gt. r1)
then
3457 rs(k) = (qs1d(k) + qsten(k)*dt)*rho(k)
3464 if ((qg1d(k) + qgten(k)*dt) .gt. r1)
then
3465 rg(k) = (qg1d(k) + qgten(k)*dt)*rho(k)
3477 if (.not. iiwarm)
then
3485 if (.not. l_qs(k)) cycle
3486 tc0 = min(-0.1, temp(k)-273.15)
3487 smob(k) = rs(k)*oams
3491 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3))
then
3494 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
3495 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
3496 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
3497 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
3498 + sa(10)*bm_s*bm_s*bm_s
3500 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
3501 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
3502 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
3503 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
3504 + sb(10)*bm_s*bm_s*bm_s
3505 smo2(k) = (smob(k)/a_)**(1./b_)
3509 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
3510 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
3511 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
3512 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
3513 + sa(10)*cse(1)*cse(1)*cse(1)
3515 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
3516 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
3517 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
3518 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
3519 smoc(k) = a_ * smo2(k)**b_
3522 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(14) &
3523 + sa(4)*tc0*cse(14) + sa(5)*tc0*tc0 &
3524 + sa(6)*cse(14)*cse(14) + sa(7)*tc0*tc0*cse(14) &
3525 + sa(8)*tc0*cse(14)*cse(14) + sa(9)*tc0*tc0*tc0 &
3526 + sa(10)*cse(14)*cse(14)*cse(14)
3528 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(14) + sb(4)*tc0*cse(14) &
3529 + sb(5)*tc0*tc0 + sb(6)*cse(14)*cse(14) &
3530 + sb(7)*tc0*tc0*cse(14) + sb(8)*tc0*cse(14)*cse(14) &
3531 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(14)*cse(14)*cse(14)
3532 smod(k) = a_ * smo2(k)**b_
3545 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3547 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3548 n0_r(k) = nr(k)*org2*lamr**cre(2)
3561 if ( (ssatw(k).gt. eps) .or. (ssatw(k).lt. -eps .and. &
3563 clap = (qv(k)-qvs(k))/(1. + lvt2(k)*qvs(k))
3565 fcd = qvs(k)* exp(lvt2(k)*clap) - qv(k) + clap
3566 dfcd = qvs(k)*lvt2(k)* exp(lvt2(k)*clap) + 1.
3567 clap = clap - fcd/dfcd
3569 xrc = rc(k) + clap*rho(k)
3571 if (xrc.gt. r1)
then
3572 prw_vcd(k) = clap*odt
3574 if (clap .gt. eps)
then
3575 if (is_aerosol_aware .or. merra2_aerosol_aware)
then
3576 xnc = max(2.,
activ_ncloud(temp(k), w1d(k)+rand3, nwfa(k), lsml))
3584 pnc_wcd(k) = 0.5*(xnc-nc(k) + abs(xnc-nc(k)))*odts*orho
3587 elseif (clap .lt. -eps .AND. ssatw(k).lt.-1.e-6 .AND. &
3588 is_aerosol_aware)
then
3589 tempc = temp(k) - 273.15
3592 rvs_p = rvs*otemp*(lvap(k)*otemp*orv - 1.)
3593 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*orv - 1.) &
3594 *otemp*(lvap(k)*otemp*orv - 1.) &
3595 + (-2.*lvap(k)*otemp*otemp*otemp*orv) &
3597 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
3598 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
3599 * rvs_pp/rvs_p * rvs/rvs_p
3600 alphsc = max(1.e-9, alphsc)
3602 if (abs(xsat).lt. 1.e-9) xsat=0.
3603 t1_evap = 2.*pi*( 1.0 - alphsc*xsat &
3604 + 2.*alphsc*alphsc*xsat*xsat &
3605 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
3608 dc_star = sqrt(-2.0_dp*dt * t1_evap/(2.*pi) &
3609 * 4.*diffu(k)*ssatw(k)*rvs/rho_w)
3610 idx_d = max(1, min(int(1.e6*dc_star), nbc))
3612 idx_n = nint(1.0 + real(nbc, kind=wp) * log(real(nc(k)/t_nc(1), kind=dp)) / nic1)
3613 idx_n = max(1, min(idx_n, nbc))
3616 if (rc(k).gt. r_c(1))
then
3617 nic = nint(log10(rc(k)))
3618 do_loop_rc_cond:
do nn = nic-1, nic+1
3620 if ( (rc(k)/10.**nn).ge.1.0 .and. (rc(k)/10.**nn).lt.10.0 )
exit do_loop_rc_cond
3621 enddo do_loop_rc_cond
3622 idx_c = int(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
3623 idx_c = max(1, min(idx_c, ntb_c))
3630 prw_vcd(k) = max(real(-rc(k)*0.99*orho*odt, kind=dp), prw_vcd(k))
3631 pnc_wcd(k) = max(real(-nc(k)*0.99*orho*odt, kind=dp), &
3632 -tnc_wev(idx_d, idx_c, idx_n)*orho*odt)
3636 prw_vcd(k) = -rc(k)*orho*odt
3637 pnc_wcd(k) = -nc(k)*orho*odt
3642 qvten(k) = qvten(k) - prw_vcd(k)
3643 qcten(k) = qcten(k) + prw_vcd(k)
3644 ncten(k) = ncten(k) + pnc_wcd(k)
3645 if (is_aerosol_aware) &
3646 nwfaten(k) = nwfaten(k) - pnc_wcd(k)
3647 tten(k) = tten(k) + lvap(k)*ocp(k)*prw_vcd(k)*(1-ifdry)
3648 rc(k) = max(r1, (qc1d(k) + dt*qcten(k))*rho(k))
3649 if (rc(k).eq.r1) l_qc(k) = .false.
3650 nc(k) = max(2., min((nc1d(k)+ncten(k)*dt)*rho(k), nt_c_max))
3651 if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware))
then
3658 qv(k) = max(1.e-10, qv1d(k) + dt*qvten(k))
3659 temp(k) = t1d(k) + dt*tten(k)
3660 rho(k) = roverrv*pres(k) / (r*temp(k)*(qv(k)+roverrv))
3661 qvs(k) =
rslf(pres(k), temp(k))
3662 ssatw(k) = qv(k)/qvs(k) - 1.
3671 if ( (ssatw(k).lt. -eps) .and. l_qr(k) &
3672 .and. (.not.(prw_vcd(k).gt. 0.)) )
then
3673 tempc = temp(k) - 273.15
3676 rhof(k) = sqrt(rho_not*orho)
3677 rhof2(k) = sqrt(rhof(k))
3678 diffu(k) = 2.11e-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
3679 if (tempc .ge. 0.0)
then
3680 visco(k) = (1.718+0.0049*tempc)*1.0e-5
3682 visco(k) = (1.718+0.0049*tempc-1.2e-5*tempc*tempc)*1.0e-5
3684 vsc2(k) = sqrt(rho(k)/visco(k))
3685 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
3686 tcond(k) = (5.69 + 0.0168*tempc)*1.0e-5 * 418.936
3687 ocp(k) = 1./(cp*(1.+0.887*qv(k)))
3690 rvs_p = rvs*otemp*(lvap(k)*otemp*orv - 1.)
3691 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*orv - 1.) &
3692 *otemp*(lvap(k)*otemp*orv - 1.) &
3693 + (-2.*lvap(k)*otemp*otemp*otemp*orv) &
3695 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
3696 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
3697 * rvs_pp/rvs_p * rvs/rvs_p
3698 alphsc = max(1.e-9, alphsc)
3699 xsat = min(-1.e-9, ssatw(k))
3700 t1_evap = 2.*pi*( 1.0 - alphsc*xsat &
3701 + 2.*alphsc*alphsc*xsat*xsat &
3702 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
3707 if (qv(k)/qvs(k) .lt. 0.95 .AND. rr(k)*orho.le.1.e-8)
then
3708 prv_rev(k) = rr(k)*orho*odts
3710 prv_rev(k) = t1_evap*diffu(k)*(-ssatw(k))*n0_r(k)*rvs &
3711 * (t1_qr_ev*ilamr(k)**cre(10) &
3712 + t2_qr_ev*vsc2(k)*rhof2(k)*((lamr+0.5*fv_r)**(-cre(11))))
3713 rate_max = min((rr(k)*orho*odts), (qvs(k)-qv(k))*odts)
3714 prv_rev(k) = min(real(rate_max, kind=dp), prv_rev(k)*orho)
3723 if (prr_gml(k).gt.0.0)
then
3724 eva_factor = min(1.0, 0.01+(0.99-0.01)*(tempc/20.0))
3725 prv_rev(k) = prv_rev(k)*eva_factor
3729 pnr_rev(k) = min(real(nr(k)*0.99*orho*odts, kind=dp), &
3730 prv_rev(k) * nr(k)/rr(k))
3732 qrten(k) = qrten(k) - prv_rev(k)
3733 qvten(k) = qvten(k) + prv_rev(k)
3734 nrten(k) = nrten(k) - pnr_rev(k)
3735 if (is_aerosol_aware) &
3736 nwfaten(k) = nwfaten(k) + pnr_rev(k)
3737 tten(k) = tten(k) - lvap(k)*ocp(k)*prv_rev(k)*(1-ifdry)
3739 rr(k) = max(r1, (qr1d(k) + dt*qrten(k))*rho(k))
3740 qv(k) = max(1.e-10, qv1d(k) + dt*qvten(k))
3741 nr(k) = max(r2, (nr1d(k) + dt*nrten(k))*rho(k))
3742 temp(k) = t1d(k) + dt*tten(k)
3743 rho(k) = roverrv*pres(k) / (r*temp(k)*(qv(k)+roverrv))
3746#if ( WRF_CHEM == 1 )
3748 evapprod(k) = prv_rev(k) - (min(zerod0,prs_sde(k)) + &
3749 min(zerod0,prg_gde(k)))
3750 rainprod(k) = prr_wau(k) + prr_rcw(k) + prs_scw(k) + &
3751 prg_scw(k) + prs_iau(k) + &
3752 prg_gcw(k) + prs_sci(k) + &
3768 do k = kte+1, kts, -1
3779 if (any(l_qr .eqv. .true.))
then
3780 fallspeed_adjustment_factor=1.0
3781 if (
present(fs_fac_rain) ) fallspeed_adjustment_factor=fs_fac_rain
3785 rhof(k) = sqrt(rho_not/rho(k))
3787 if (rr(k).gt. r1)
then
3788 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3789 vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) &
3790 *((lamr+fv_r)**(-cre(6)))
3791 vtrk(k) = vtr*fallspeed_adjustment_factor
3797 vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) &
3798 *((lamr+fv_r)**(-cre(7)))
3799 vtnrk(k) = vtr*fallspeed_adjustment_factor
3802 vtnrk(k) = vtnrk(k+1)
3805 if (max(vtrk(k),vtnrk(k)) .gt. 1.e-3)
then
3806 ksed1(1) = max(ksed1(1), k)
3807 delta_tp = dzq(k)/(max(vtrk(k),vtnrk(k)))
3808 nstep = max(nstep, int(dt/delta_tp + 1.))
3811 if (ksed1(1) .eq. kte) ksed1(1) = kte-1
3812 if (nstep .gt. 0) onstep(1) = 1./real(nstep, kind=wp)
3817 if (any(l_qc .eqv. .true.))
then
3819 do_loop_hgt_agl :
do k = kts, kte-1
3820 if (rc(k) .gt. r2) ksed1(5) = k
3821 hgt_agl = hgt_agl + dzq(k)
3822 if (hgt_agl .gt. 500.0)
exit do_loop_hgt_agl
3823 enddo do_loop_hgt_agl
3825 do k = ksed1(5), kts, -1
3827 if (rc(k) .gt. r1 .and. w1d(k) .lt. 1.e-1)
then
3828 if (nc(k).gt.10000.e6)
then
3830 elseif (nc(k).lt.100.)
then
3833 nu_c = nint(1000.e6/nc(k)) + 2
3834 nu_c = max(2, min(nu_c+nint(rand2), 15))
3836 lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
3838 vtc = rhof(k)*av_c*ccg(5,nu_c)*ocg2(nu_c) * ilamc**bv_c
3840 vtc = rhof(k)*av_c*ccg(4,nu_c)*ocg1(nu_c) * ilamc**bv_c
3848 if (.not. iiwarm)
then
3850 if (any(l_qi .eqv. .true.))
then
3855 if (ri(k).gt. r1)
then
3856 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
3858 vti = rhof(k)*av_i*cig(3)*oig2 * ilami**bv_i
3863 vti = rhof(k)*av_i*cig(6)/cig(7) * ilami**bv_i
3867 vtnik(k) = vtnik(k+1)
3870 if (vtik(k) .gt. 1.e-3)
then
3871 ksed1(2) = max(ksed1(2), k)
3872 delta_tp = dzq(k)/vtik(k)
3873 nstep = max(nstep, int(dt/delta_tp + 1.))
3876 if (ksed1(2) .eq. kte) ksed1(2) = kte-1
3877 if (nstep .gt. 0) onstep(2) = 1./real(nstep, kind=wp)
3882 if (any(l_qs .eqv. .true.))
then
3883 fallspeed_adjustment_factor=1.0
3884 if (
present(fs_fac_snow) ) fallspeed_adjustment_factor=fs_fac_snow
3891 if (rs(k).gt. r1)
then
3892 xds = smoc(k) / smob(k)
3894 ils1 = 1./(mrat*lam0 + fv_s)
3895 ils2 = 1./(mrat*lam1 + fv_s)
3896 t1_vts = kap0*csg(4)*ils1**cse(4)
3897 t2_vts = kap1*mrat**mu_s*csg(10)*ils2**cse(10)
3898 ils1 = 1./(mrat*lam0)
3899 ils2 = 1./(mrat*lam1)
3900 t3_vts = kap0*csg(1)*ils1**cse(1)
3901 t4_vts = kap1*mrat**mu_s*csg(7)*ils2**cse(7)
3902 vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts)
3903 vts=vts*fallspeed_adjustment_factor
3904 if (prr_sml(k) .gt. 0.0)
then
3907 sr = rs(k)/(rs(k)+rr(k))
3908 vtsk(k) = vts*sr + (1.-sr)*vtrk(k)
3911 vtsk(k) = vts*vts_boost(k)
3919 if (vtsk(k) .gt. 1.e-3)
then
3920 ksed1(3) = max(ksed1(3), k)
3921 delta_tp = dzq(k)/vtsk(k)
3922 nstep = max(nstep, int(dt/delta_tp + 1.))
3925 if (ksed1(3) .eq. kte) ksed1(3) = kte-1
3926 if (nstep .gt. 0) onstep(3) = 1./real(nstep, kind=wp)
3931 if (any(l_qg .eqv. .true.))
then
3936 if (rg(k).gt. r1)
then
3937 vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
3938 if (temp(k).gt. t_0)
then
3939 vtgk(k) = max(vtg, vtrk(k))
3947 if (vtgk(k) .gt. 1.e-3)
then
3948 ksed1(4) = max(ksed1(4), k)
3949 delta_tp = dzq(k)/vtgk(k)
3950 nstep = max(nstep, int(dt/delta_tp + 1.))
3953 if (ksed1(4) .eq. kte) ksed1(4) = kte-1
3954 if (nstep .gt. 0) onstep(4) = 1./real(nstep, kind=wp)
3964 if (any(l_qr .eqv. .true.))
then
3965 nstep = nint(1./onstep(1))
3967 if(.not. sedi_semi)
then
3970 sed_r(k) = vtrk(k)*rr(k)
3971 sed_n(k) = vtnrk(k)*nr(k)
3976 qrten(k) = qrten(k) - sed_r(k)*odzq*onstep(1)*orho
3977 nrten(k) = nrten(k) - sed_n(k)*odzq*onstep(1)*orho
3978 rr(k) = max(r1, rr(k) - sed_r(k)*odzq*dt*onstep(1))
3979 nr(k) = max(r2, nr(k) - sed_n(k)*odzq*dt*onstep(1))
3980 pfll1(k) = pfll1(k) + sed_r(k)*dt*onstep(1)
3981 do k = ksed1(1), kts, -1
3984 qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) &
3985 *odzq*onstep(1)*orho
3986 nrten(k) = nrten(k) + (sed_n(k+1)-sed_n(k)) &
3987 *odzq*onstep(1)*orho
3988 rr(k) = max(r1, rr(k) + (sed_r(k+1)-sed_r(k)) &
3990 nr(k) = max(r2, nr(k) + (sed_n(k+1)-sed_n(k)) &
3992 pfll1(k) = pfll1(k) + sed_r(k)*dt*onstep(1)
3995 if (rr(kts).gt.r1*rr_min)
then
3996 pptrain = pptrain + sed_r(kts)*dt*onstep(1)
4002 niter = int(nstep/max(decfl,1)) + 1
4007 call semi_lagrange_sedim(kte,dzq,vtrk,rr,rainsfc,pfll,dtcfl,r1)
4008 call semi_lagrange_sedim(kte,dzq,vtnrk,nr,vtr,pdummy,dtcfl,r2)
4010 orhodt = 1./(rho(k)*dt)
4011 qrten(k) = qrten(k) + (rr(k) - rr_tmp(k)) * orhodt
4012 nrten(k) = nrten(k) + (nr(k) - nr_tmp(k)) * orhodt
4013 pfll1(k) = pfll1(k) + pfll(k)
4015 pptrain = pptrain + rainsfc
4017 do k = kte+1, kts, -1
4023 if (rr(k).gt. r1)
then
4024 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
4025 vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) &
4026 *((lamr+fv_r)**(-cre(6)))
4033 vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) &
4034 *((lamr+fv_r)**(-cre(7)))
4044 if (any(l_qc .eqv. .true.))
then
4046 sed_c(k) = vtck(k)*rc(k)
4047 sed_n(k) = vtnck(k)*nc(k)
4049 do k = ksed1(5), kts, -1
4052 qcten(k) = qcten(k) + (sed_c(k+1)-sed_c(k)) *odzq*orho
4053 ncten(k) = ncten(k) + (sed_n(k+1)-sed_n(k)) *odzq*orho
4054 rc(k) = max(r1, rc(k) + (sed_c(k+1)-sed_c(k)) *odzq*dt)
4055 nc(k) = max(10., nc(k) + (sed_n(k+1)-sed_n(k)) *odzq*dt)
4061 if (any(l_qi .eqv. .true.))
then
4062 nstep = nint(1./onstep(2))
4065 sed_i(k) = vtik(k)*ri(k)
4066 sed_n(k) = vtnik(k)*ni(k)
4071 qiten(k) = qiten(k) - sed_i(k)*odzq*onstep(2)*orho
4072 niten(k) = niten(k) - sed_n(k)*odzq*onstep(2)*orho
4073 ri(k) = max(r1, ri(k) - sed_i(k)*odzq*dt*onstep(2))
4074 ni(k) = max(r2, ni(k) - sed_n(k)*odzq*dt*onstep(2))
4075 pfil1(k) = pfil1(k) + sed_i(k)*dt*onstep(2)
4076 do k = ksed1(2), kts, -1
4079 qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) &
4080 *odzq*onstep(2)*orho
4081 niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) &
4082 *odzq*onstep(2)*orho
4083 ri(k) = max(r1, ri(k) + (sed_i(k+1)-sed_i(k)) &
4085 ni(k) = max(r2, ni(k) + (sed_n(k+1)-sed_n(k)) &
4087 pfil1(k) = pfil1(k) + sed_i(k)*dt*onstep(2)
4090 if (ri(kts).gt.r1*rr_min)
then
4091 pptice = pptice + sed_i(kts)*dt*onstep(2)
4098 if (any(l_qs .eqv. .true.))
then
4099 nstep = nint(1./onstep(3))
4102 sed_s(k) = vtsk(k)*rs(k)
4107 qsten(k) = qsten(k) - sed_s(k)*odzq*onstep(3)*orho
4108 rs(k) = max(r1, rs(k) - sed_s(k)*odzq*dt*onstep(3))
4109 pfil1(k) = pfil1(k) + sed_s(k)*dt*onstep(3)
4110 do k = ksed1(3), kts, -1
4113 qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) &
4114 *odzq*onstep(3)*orho
4115 rs(k) = max(r1, rs(k) + (sed_s(k+1)-sed_s(k)) &
4117 pfil1(k) = pfil1(k) + sed_s(k)*dt*onstep(3)
4120 if (rs(kts).gt.r1*rr_min)
then
4121 pptsnow = pptsnow + sed_s(kts)*dt*onstep(3)
4128 if (any(l_qg .eqv. .true.))
then
4129 nstep = nint(1./onstep(4))
4130 if(.not. sedi_semi)
then
4133 sed_g(k) = vtgk(k)*rg(k)
4138 qgten(k) = qgten(k) - sed_g(k)*odzq*onstep(4)*orho
4139 rg(k) = max(r1, rg(k) - sed_g(k)*odzq*dt*onstep(4))
4140 pfil1(k) = pfil1(k) + sed_g(k)*dt*onstep(4)
4141 do k = ksed1(4), kts, -1
4144 qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) &
4145 *odzq*onstep(4)*orho
4146 rg(k) = max(r1, rg(k) + (sed_g(k+1)-sed_g(k)) &
4148 pfil1(k) = pfil1(k) + sed_g(k)*dt*onstep(4)
4151 if (rg(kts).gt.r1*rr_min)
then
4152 pptgraul = pptgraul + sed_g(kts)*dt*onstep(4)
4158 niter = int(nstep/max(decfl,1)) + 1
4163 call semi_lagrange_sedim(kte,dzq,vtgk,rg,graulsfc,pfil,dtcfl,r1)
4165 orhodt = 1./(rho(k)*dt)
4166 qgten(k) = qgten(k) + (rg(k) - rg_tmp(k))*orhodt
4167 pfil1(k) = pfil1(k) + pfil(k)
4169 pptgraul = pptgraul + graulsfc
4170 do k = kte+1, kts, -1
4175 if (rg(k).gt. r1)
then
4176 ygra1 = log10(max(1.e-9_wp, rg(k)))
4177 zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1
4178 n0_exp = 10.**(zans1)
4179 n0_exp = max(real(gonv_min, kind=dp), min(n0_exp, real(gonv_max, kind=dp)))
4180 lam_exp = (n0_exp*am_g*cgg(1)/rg(k))**oge1
4181 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
4183 vtg = rhof(k)*av_g*cgg(6)*ogg3 * (1./lamg)**bv_g
4184 if (temp(k).gt. t_0)
then
4185 vtgk(k) = max(vtg, vtrk(k))
4199 if (.not. iiwarm)
then
4201 xri = max(0.0, qi1d(k) + qiten(k)*dt)
4202 if ( (temp(k).gt. t_0) .and. (xri.gt. 0.0) )
then
4203 qcten(k) = qcten(k) + xri*odt
4204 ncten(k) = ncten(k) + ni1d(k)*odt
4205 qiten(k) = qiten(k) - xri*odt
4206 niten(k) = -ni1d(k)*odt
4207 tten(k) = tten(k) - lfus*ocp(k)*xri*odt*(1-ifdry)
4212 xrc = max(0.0, qc1d(k) + qcten(k)*dt)
4213 if ( (temp(k).lt. hgfr) .and. (xrc.gt. 0.0) )
then
4214 lfus2 = lsub - lvap(k)
4215 xnc = nc1d(k) + ncten(k)*dt
4216 qiten(k) = qiten(k) + xrc*odt
4217 niten(k) = niten(k) + xnc*odt
4218 qcten(k) = qcten(k) - xrc*odt
4219 ncten(k) = ncten(k) - xnc*odt
4220 tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-ifdry)
4231 t1d(k) = t1d(k) + tten(k)*dt
4232 qv1d(k) = max(1.e-10, qv1d(k) + qvten(k)*dt)
4233 qc1d(k) = qc1d(k) + qcten(k)*dt
4234 nc1d(k) = max(2./rho(k), min(nc1d(k) + ncten(k)*dt, nt_c_max))
4235 if (is_aerosol_aware)
then
4236 nwfa1d(k) = max(11.1e6, min(9999.e6, &
4237 (nwfa1d(k)+nwfaten(k)*dt)))
4238 nifa1d(k) = max(nain1*0.01, min(9999.e6, &
4239 (nifa1d(k)+nifaten(k)*dt)))
4241 if (qc1d(k) .le. r1)
then
4245 if (nc1d(k)*rho(k).gt.10000.e6)
then
4247 elseif (nc1d(k)*rho(k).lt.100.)
then
4250 nu_c = nint(1000.e6/(nc1d(k)*rho(k))) + 2
4251 nu_c = max(2, min(nu_c+nint(rand2), 15))
4253 lamc = (am_r*ccg(2,nu_c)*ocg1(nu_c)*nc1d(k)/qc1d(k))**obmr
4254 xdc = (bm_r + nu_c + 1.) / lamc
4255 if (xdc.lt. d0c)
then
4256 lamc = cce(2,nu_c)/d0c
4257 elseif (xdc.gt. d0r*2.)
then
4258 lamc = cce(2,nu_c)/(d0r*2.)
4260 nc1d(k) = min(ccg(1,nu_c)*ocg2(nu_c)*qc1d(k)/am_r*lamc**bm_r,&
4261 real(Nt_c_max, kind=dp)/rho(k))
4264 qi1d(k) = qi1d(k) + qiten(k)*dt
4265 ni1d(k) = max(r2/rho(k), ni1d(k) + niten(k)*dt)
4266 if (qi1d(k) .le. r1)
then
4270 lami = (am_i*cig(2)*oig1*ni1d(k)/qi1d(k))**obmi
4272 xdi = (bm_i + mu_i + 1.) * ilami
4273 if (xdi.lt. 5.e-6)
then
4275 elseif (xdi.gt. 300.e-6)
then
4276 lami = cie(2)/300.e-6
4278 ni1d(k) = min(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, &
4281 qr1d(k) = qr1d(k) + qrten(k)*dt
4282 nr1d(k) = max(r2/rho(k), nr1d(k) + nrten(k)*dt)
4283 if (qr1d(k) .le. r1)
then
4287 lamr = (am_r*crg(3)*org2*nr1d(k)/qr1d(k))**obmr
4288 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
4289 if (mvd_r(k) .gt. 2.5e-3)
then
4291 elseif (mvd_r(k) .lt. d0r*0.75)
then
4294 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
4295 nr1d(k) = crg(2)*org3*qr1d(k)*lamr**bm_r / am_r
4297 qs1d(k) = qs1d(k) + qsten(k)*dt
4298 if (qs1d(k) .le. r1) qs1d(k) = 0.0
4299 qg1d(k) = qg1d(k) + qgten(k)*dt
4300 if (qg1d(k) .le. r1) qg1d(k) = 0.0
4304 calculate_extended_diagnostics:
if (ext_diag)
then
4306 if(prw_vcd(k).gt.0)
then
4307 prw_vcdc1(k) = prw_vcd(k)*dt
4308 elseif(prw_vcd(k).lt.0)
then
4309 prw_vcde1(k) = -1*prw_vcd(k)*dt
4312 tpri_inu1(k) = pri_inu(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4314 if(pri_ide(k).gt.0)
then
4315 tpri_ide1_d(k) = pri_ide(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4317 tpri_ide1_s(k) = -pri_ide(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4320 if(temp(k).lt.t_0)
then
4321 tprs_ide1(k) = prs_ide(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4324 if(prs_sde(k).gt.0)
then
4325 tprs_sde1_d(k) = prs_sde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4327 tprs_sde1_s(k) = -prs_sde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4330 if(prg_gde(k).gt.0)
then
4331 tprg_gde1_d(k) = prg_gde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4333 tprg_gde1_s(k) = -prg_gde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4336 tpri_iha1(k) = pri_iha(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4337 tpri_wfz1(k) = pri_wfz(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4338 tpri_rfz1(k) = pri_rfz(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4339 tprg_rfz1(k) = prg_rfz(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4340 tprs_scw1(k) = prs_scw(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4341 tprg_scw1(k) = prg_scw(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4342 tprg_rcs1(k) = prg_rcs(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4344 if(temp(k).lt.t_0)
then
4345 tprs_rcs1(k) = prs_rcs(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4348 tprr_rci1(k) = prr_rci(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4350 if(temp(k).lt.t_0)
then
4351 tprg_rcg1(k) = prg_rcg(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4354 if(prw_vcd(k).gt.0)
then
4355 tprw_vcd1_c(k) = lvap(k)*ocp(k)*prw_vcd(k)*(1-ifdry)*dt
4357 tprw_vcd1_e(k) = -lvap(k)*ocp(k)*prw_vcd(k)*(1-ifdry)*dt
4361 tprr_sml1(k) = prr_sml(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
4362 tprr_gml1(k) = prr_gml(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
4364 if(temp(k).ge.t_0)
then
4365 tprr_rcg1(k) = -prr_rcg(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
4368 if(temp(k).ge.t_0)
then
4369 tprr_rcs1(k) = -prr_rcs(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
4372 tprv_rev1(k) = lvap(k)*ocp(k)*prv_rev(k)*(1-ifdry)*dt
4373 tten1(k) = tten(k)*dt
4374 qvten1(k) = qvten(k)*dt
4375 qiten1(k) = qiten(k)*dt
4376 qrten1(k) = qrten(k)*dt
4377 qsten1(k) = qsten(k)*dt
4378 qgten1(k) = qgten(k)*dt
4379 niten1(k) = niten(k)*dt
4380 nrten1(k) = nrten(k)*dt
4381 ncten1(k) = ncten(k)*dt
4382 qcten1(k) = qcten(k)*dt
4384 endif calculate_extended_diagnostics