985 nwfa, nifa, nwfa2d, nifa2d, &
987 p, w, dz, dt_in, dt_inner, &
988 sedi_semi, decfl, lsm, &
992 GRAUPELNC, GRAUPELNCV, SR, &
994 rainprod, evapprod, &
996 refl_10cm, diagflag, do_radar_ref, &
998 vt_dbz_wt, first_time_step, &
999 re_cloud, re_ice, re_snow, &
1000 has_reqc, has_reqi, has_reqs, &
1001 aero_ind_fdb, rand_perturb_on, &
1003 rand_pert, spp_prt_list, spp_var_list, &
1004 spp_stddev_cutoff, n_var_spp, &
1005 ids,ide, jds,jde, kds,kde, &
1006 ims,ime, jms,jme, kms,kme, &
1007 its,ite, jts,jte, kts,kte, &
1008 fullradar_diag, istep, nsteps, &
1015 prw_vcde, tpri_inu, tpri_ide_d, &
1016 tpri_ide_s, tprs_ide, tprs_sde_d, &
1017 tprs_sde_s, tprg_gde_d, &
1018 tprg_gde_s, tpri_iha, tpri_wfz, &
1019 tpri_rfz, tprg_rfz, tprs_scw, tprg_scw, &
1020 tprg_rcs, tprs_rcs, &
1021 tprr_rci, tprg_rcg, &
1022 tprw_vcd_c, tprw_vcd_e, tprr_sml, &
1023 tprr_gml, tprr_rcg, &
1024 tprr_rcs, tprv_rev, tten3, qvten3, &
1025 qrten3, qsten3, qgten3, qiten3, niten3, &
1026 nrten3, ncten3, qcten3, &
1032 INTEGER,
INTENT(IN):: ids,ide, jds,jde, kds,kde, &
1033 ims,ime, jms,jme, kms,kme, &
1034 its,ite, jts,jte, kts,kte
1035 REAL,
DIMENSION(ims:ime, kms:kme, jms:jme),
INTENT(INOUT):: &
1036 qv, qc, qr, qi, qs, qg, ni, nr
1037 REAL,
DIMENSION(ims:ime, kms:kme, jms:jme),
OPTIONAL,
INTENT(INOUT):: &
1039 REAL,
DIMENSION(ims:ime, kms:kme, jms:jme),
OPTIONAL,
INTENT(IN):: &
1041 REAL,
DIMENSION(ims:ime, kms:kme, jms:jme),
OPTIONAL,
INTENT(INOUT):: &
1043 REAL,
DIMENSION(ims:ime, jms:jme),
OPTIONAL,
INTENT(IN):: nwfa2d, nifa2d
1044 INTEGER,
DIMENSION(ims:ime, jms:jme),
INTENT(IN):: lsm
1045 REAL,
DIMENSION(ims:ime, kms:kme, jms:jme),
OPTIONAL,
INTENT(INOUT):: &
1046 re_cloud, re_ice, re_snow
1047 REAL,
DIMENSION(ims:ime, kms:kme, jms:jme),
INTENT(INOUT):: pfils, pflls
1048 INTEGER,
INTENT(IN) :: rand_perturb_on, kme_stoch, n_var_spp
1049 REAL,
DIMENSION(:,:),
INTENT(IN),
OPTIONAL :: rand_pert
1050 REAL,
DIMENSION(:),
INTENT(IN),
OPTIONAL :: spp_prt_list
1051 REAL,
DIMENSION(:),
INTENT(IN) :: spp_stddev_cutoff
1052 CHARACTER(len=10),
DIMENSION(:),
INTENT(IN),
OPTIONAL :: spp_var_list
1053 INTEGER,
INTENT(IN):: has_reqc, has_reqi, has_reqs
1054#if ( WRF_CHEM == 1 )
1055 REAL,
DIMENSION(ims:ime, kms:kme, jms:jme),
INTENT(INOUT):: &
1058 REAL,
DIMENSION(ims:ime, kms:kme, jms:jme),
INTENT(IN):: &
1060 REAL,
DIMENSION(ims:ime, jms:jme),
INTENT(INOUT):: &
1062 REAL,
DIMENSION(ims:ime, jms:jme),
OPTIONAL,
INTENT(INOUT):: &
1065 GRAUPELNC, GRAUPELNCV
1066 REAL,
DIMENSION(ims:ime, kms:kme, jms:jme),
INTENT(INOUT):: &
1068 REAL,
DIMENSION(ims:ime, jms:jme),
INTENT(INOUT):: &
1070 REAL,
DIMENSION(ims:ime, kms:kme, jms:jme),
OPTIONAL,
INTENT(INOUT):: &
1072 LOGICAL,
INTENT(IN) :: first_time_step
1073 REAL,
INTENT(IN):: dt_in, dt_inner
1074 LOGICAL,
INTENT(IN) :: sedi_semi
1075 INTEGER,
INTENT(IN) :: decfl
1077 INTEGER,
INTENT (IN) :: istep, nsteps
1078 LOGICAL,
INTENT (IN) :: fullradar_diag
1080 LOGICAL,
INTENT (IN) :: ext_diag
1081 LOGICAL,
OPTIONAL,
INTENT(IN):: aero_ind_fdb
1082 REAL,
DIMENSION(:,:,:),
INTENT(INOUT),
OPTIONAL :: &
1085 prw_vcde, tpri_inu, tpri_ide_d, &
1086 tpri_ide_s, tprs_ide, &
1087 tprs_sde_d, tprs_sde_s, tprg_gde_d, &
1088 tprg_gde_s, tpri_iha, tpri_wfz, &
1089 tpri_rfz, tprg_rfz, tprs_scw, tprg_scw, &
1090 tprg_rcs, tprs_rcs, &
1091 tprr_rci, tprg_rcg, &
1092 tprw_vcd_c, tprw_vcd_e, tprr_sml, &
1093 tprr_gml, tprr_rcg, &
1094 tprr_rcs, tprv_rev, tten3, qvten3, &
1095 qrten3, qsten3, qgten3, qiten3, niten3, &
1096 nrten3, ncten3, qcten3
1099 REAL,
DIMENSION(kts:kte):: &
1100 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
1101 nr1d, nc1d, nwfa1d, nifa1d, &
1102 t1d, p1d, w1d, dz1d, rho, dBZ, pfil1, pfll1
1104 REAL,
DIMENSION(:),
ALLOCATABLE:: &
1107 prw_vcde1, tpri_inu1, tpri_ide1_d, &
1108 tpri_ide1_s, tprs_ide1, &
1109 tprs_sde1_d, tprs_sde1_s, tprg_gde1_d, &
1110 tprg_gde1_s, tpri_iha1, tpri_wfz1, &
1111 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1,&
1112 tprg_rcs1, tprs_rcs1, &
1113 tprr_rci1, tprg_rcg1, &
1114 tprw_vcd1_c, tprw_vcd1_e, tprr_sml1, &
1115 tprr_gml1, tprr_rcg1, &
1116 tprr_rcs1, tprv_rev1, tten1, qvten1, &
1117 qrten1, qsten1, qgten1, qiten1, niten1, &
1118 nrten1, ncten1, qcten1
1120 REAL,
DIMENSION(kts:kte):: re_qc1d, re_qi1d, re_qs1d
1121#if ( WRF_CHEM == 1 )
1122 REAL,
DIMENSION(kts:kte):: &
1123 rainprod1d, evapprod1d
1125 REAL,
DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
1126 REAL:: dt, pptrain, pptsnow, pptgraul, pptice
1127 REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max
1129 REAL:: rand1, rand2, rand3, rand_pert_max
1130 INTEGER:: i, j, k, m
1131 INTEGER:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr
1132 INTEGER:: jmax_qc,jmax_qr,jmax_qi,jmax_qs,jmax_qg,jmax_ni,jmax_nr
1133 INTEGER:: kmax_qc,kmax_qr,kmax_qi,kmax_qs,kmax_qg,kmax_ni,kmax_nr
1134 INTEGER:: i_start, j_start, i_end, j_end
1135 LOGICAL,
OPTIONAL,
INTENT(IN) :: diagflag
1136 INTEGER,
OPTIONAL,
INTENT(IN) :: do_radar_ref
1137 logical :: melti = .false.
1141 character(len=*),
optional,
intent( out) :: errmsg
1142 integer,
optional,
intent( out) :: errflg
1145 if (
present(errmsg)) errmsg =
''
1146 if (
present(errflg)) errflg = 0
1149 test_only_once:
if (first_time_step .and. istep==1)
then
1152 if ( (
present(tt) .and. (
present(th) .or.
present(pii))) .or. &
1153 (.not.
present(tt) .and. .not.(
present(th) .and.
present(pii))) )
then
1154 if (
present(errmsg) .and.
present(errflg))
then
1155 write(errmsg,
'(a)')
'Logic error in mp_gt_driver: provide either tt or th+pii'
1159 write(*,
'(a)')
'Logic error in mp_gt_driver: provide either tt or th+pii'
1164 if (is_aerosol_aware .and. (.not.
present(nc) .or. &
1165 .not.
present(nwfa) .or. &
1166 .not.
present(nifa) .or. &
1167 .not.
present(nwfa2d) .or. &
1168 .not.
present(nifa2d) ))
then
1169 if (
present(errmsg) .and.
present(errflg))
then
1170 write(errmsg,
'(*(a))')
'Logic error in mp_gt_driver: provide nc, nwfa, nifa, nwfa2d', &
1171 ' and nifa2d for aerosol-aware version of Thompson microphysics'
1175 write(*,
'(*(a))')
'Logic error in mp_gt_driver: provide nc, nwfa, nifa, nwfa2d', &
1176 ' and nifa2d for aerosol-aware version of Thompson microphysics'
1179 else if (merra2_aerosol_aware .and. (.not.
present(nc) .or. &
1180 .not.
present(nwfa) .or. &
1181 .not.
present(nifa) ))
then
1182 if (
present(errmsg) .and.
present(errflg))
then
1183 write(errmsg,
'(*(a))')
'Logic error in mp_gt_driver: provide nc, nwfa, and nifa', &
1184 ' for merra2 aerosol-aware version of Thompson microphysics'
1188 write(*,
'(*(a))')
'Logic error in mp_gt_driver: provide nc, nwfa, and nifa', &
1189 ' for merra2 aerosol-aware version of Thompson microphysics'
1192 else if (.not.is_aerosol_aware .and. .not.merra2_aerosol_aware .and. &
1193 (
present(nwfa) .or.
present(nifa) .or.
present(nwfa2d) .or.
present(nifa2d)))
then
1194 write(*,*)
'WARNING, nc/nwfa/nifa/nwfa2d/nifa2d present but is_aerosol_aware/merra2_aerosol_aware are FALSE'
1196 end if test_only_once
1202 allocate_extended_diagnostics:
if (ext_diag)
then
1203 allocate (prw_vcdc1(kts:kte))
1204 allocate (prw_vcde1(kts:kte))
1205 allocate (tpri_inu1(kts:kte))
1206 allocate (tpri_ide1_d(kts:kte))
1207 allocate (tpri_ide1_s(kts:kte))
1208 allocate (tprs_ide1(kts:kte))
1209 allocate (tprs_sde1_d(kts:kte))
1210 allocate (tprs_sde1_s(kts:kte))
1211 allocate (tprg_gde1_d(kts:kte))
1212 allocate (tprg_gde1_s(kts:kte))
1213 allocate (tpri_iha1(kts:kte))
1214 allocate (tpri_wfz1(kts:kte))
1215 allocate (tpri_rfz1(kts:kte))
1216 allocate (tprg_rfz1(kts:kte))
1217 allocate (tprs_scw1(kts:kte))
1218 allocate (tprg_scw1(kts:kte))
1219 allocate (tprg_rcs1(kts:kte))
1220 allocate (tprs_rcs1(kts:kte))
1221 allocate (tprr_rci1(kts:kte))
1222 allocate (tprg_rcg1(kts:kte))
1223 allocate (tprw_vcd1_c(kts:kte))
1224 allocate (tprw_vcd1_e(kts:kte))
1225 allocate (tprr_sml1(kts:kte))
1226 allocate (tprr_gml1(kts:kte))
1227 allocate (tprr_rcg1(kts:kte))
1228 allocate (tprr_rcs1(kts:kte))
1229 allocate (tprv_rev1(kts:kte))
1230 allocate (tten1(kts:kte))
1231 allocate (qvten1(kts:kte))
1232 allocate (qrten1(kts:kte))
1233 allocate (qsten1(kts:kte))
1234 allocate (qgten1(kts:kte))
1235 allocate (qiten1(kts:kte))
1236 allocate (niten1(kts:kte))
1237 allocate (nrten1(kts:kte))
1238 allocate (ncten1(kts:kte))
1239 allocate (qcten1(kts:kte))
1240 end if allocate_extended_diagnostics
1261 graupelnc(:,:) = 0.0
1269 ndt = max(nint(dt_in/dt_inner),1)
1271 if(dt_in .le. dt_inner) dt= dt_in
1276 if (rand_perturb_on .ne. 0)
then
1278 select case (spp_var_list(k))
1280 rand_pert_max = spp_prt_list(k)*spp_stddev_cutoff(k)
1316 j_loop:
do j = j_start, j_end
1317 i_loop:
do i = i_start, i_end
1335 if (rand_perturb_on .ne. 0)
then
1336 if (mod(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,1)
1337 m = rshift(abs(rand_perturb_on),1)
1338 if (mod(m,2) .ne. 0) rand2 = rand_pert(i,1)*2.
1339 m = rshift(abs(rand_perturb_on),2)
1340 if (mod(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,1)+rand_pert_max)
1341 m = rshift(abs(rand_perturb_on),3)
1350 IF (
PRESENT (snowncv) )
THEN
1353 IF (
PRESENT (icencv) )
THEN
1356 IF (
PRESENT (graupelncv) )
THEN
1357 graupelncv(i,j) = 0.
1362 if (
present(tt))
then
1365 t1d(k) = th(i,k,j)*pii(i,k,j)
1378 rho(k) = 0.622*p1d(k)/(r*t1d(k)*(qv1d(k)+0.622))
1384 initialize_extended_diagnostics:
if (ext_diag)
then
1422 endif initialize_extended_diagnostics
1425 if (is_aerosol_aware .or. merra2_aerosol_aware)
then
1428 nwfa1d(k) = nwfa(i,k,j)
1429 nifa1d(k) = nifa(i,k,j)
1434 nc1d(k) = nt_c_l/rho(k)
1436 nc1d(k) = nt_c_o/rho(k)
1439 nifa1d(k) = nain1*0.01
1444 call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
1445 nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dz1d, &
1446 lsml, pptrain, pptsnow, pptgraul, pptice, &
1447#if ( WRF_CHEM == 1 )
1448 rainprod1d, evapprod1d, &
1450 rand1, rand2, rand3, &
1451 kts, kte, dt, i, j, ext_diag, &
1454 prw_vcdc1, prw_vcde1, &
1455 tpri_inu1, tpri_ide1_d, tpri_ide1_s, tprs_ide1, &
1456 tprs_sde1_d, tprs_sde1_s, &
1457 tprg_gde1_d, tprg_gde1_s, tpri_iha1, tpri_wfz1, &
1458 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1, &
1459 tprg_rcs1, tprs_rcs1, tprr_rci1, &
1460 tprg_rcg1, tprw_vcd1_c, &
1461 tprw_vcd1_e, tprr_sml1, tprr_gml1, tprr_rcg1, &
1462 tprr_rcs1, tprv_rev1, &
1463 tten1, qvten1, qrten1, qsten1, &
1464 qgten1, qiten1, niten1, nrten1, ncten1, qcten1, &
1467 pcp_ra(i,j) = pcp_ra(i,j) + pptrain
1468 pcp_sn(i,j) = pcp_sn(i,j) + pptsnow
1469 pcp_gr(i,j) = pcp_gr(i,j) + pptgraul
1470 pcp_ic(i,j) = pcp_ic(i,j) + pptice
1471 rainncv(i,j) = pptrain + pptsnow + pptgraul + pptice
1472 rainnc(i,j) = rainnc(i,j) + pptrain + pptsnow + pptgraul + pptice
1473 IF (
PRESENT(snowncv) .AND.
PRESENT(snownc) )
THEN
1475 IF ( .NOT.
PRESENT(icencv) .OR. .NOT.
PRESENT(icenc) )
THEN
1476 snowncv(i,j) = pptsnow + pptice
1477 snownc(i,j) = snownc(i,j) + pptsnow + pptice
1479 snowncv(i,j) = pptsnow
1480 snownc(i,j) = snownc(i,j) + pptsnow
1484 IF (
PRESENT(icencv) .AND.
PRESENT(icenc) )
THEN
1485 icencv(i,j) = pptice
1486 icenc(i,j) = icenc(i,j) + pptice
1488 IF (
PRESENT(graupelncv) .AND.
PRESENT(graupelnc) )
THEN
1489 graupelncv(i,j) = pptgraul
1490 graupelnc(i,j) = graupelnc(i,j) + pptgraul
1492 sr(i,j) = (pptsnow + pptgraul + pptice)/(rainncv(i,j)+1.e-12)
1499 if (is_aerosol_aware)
then
1500 if (
PRESENT (aero_ind_fdb) )
then
1501 if ( .not. aero_ind_fdb)
then
1502 nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt
1503 nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt
1506 nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt
1507 nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt
1512 nwfa(i,k,j) = nwfa1d(k)
1513 nifa(i,k,j) = nifa1d(k)
1517 if (merra2_aerosol_aware)
then
1520 nwfa(i,k,j) = nwfa1d(k)
1521 nifa(i,k,j) = nifa1d(k)
1534 pfils(i,k,j) = pfils(i,k,j) + pfil1(k)
1535 pflls(i,k,j) = pflls(i,k,j) + pfll1(k)
1536 if (
present(tt))
then
1539 th(i,k,j) = t1d(k)/pii(i,k,j)
1541#if ( WRF_CHEM == 1 )
1542 rainprod(i,k,j) = rainprod1d(k)
1543 evapprod(i,k,j) = evapprod1d(k)
1545 if (qc1d(k) .gt. qc_max)
then
1550 elseif (qc1d(k) .lt. 0.0)
then
1551 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qc ', qc1d(k), &
1554 if (qr1d(k) .gt. qr_max)
then
1559 elseif (qr1d(k) .lt. 0.0)
then
1560 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qr ', qr1d(k), &
1563 if (nr1d(k) .gt. nr_max)
then
1568 elseif (nr1d(k) .lt. 0.0)
then
1569 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative nr ', nr1d(k), &
1572 if (qs1d(k) .gt. qs_max)
then
1577 elseif (qs1d(k) .lt. 0.0)
then
1578 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qs ', qs1d(k), &
1581 if (qi1d(k) .gt. qi_max)
then
1586 elseif (qi1d(k) .lt. 0.0)
then
1587 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qi ', qi1d(k), &
1590 if (qg1d(k) .gt. qg_max)
then
1595 elseif (qg1d(k) .lt. 0.0)
then
1596 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qg ', qg1d(k), &
1599 if (ni1d(k) .gt. ni_max)
then
1604 elseif (ni1d(k) .lt. 0.0)
then
1605 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative ni ', ni1d(k), &
1608 if (qv1d(k) .lt. 0.0)
then
1609 write(*,
'(a,e16.7,a,3i8)')
'WARNING, negative qv ', qv1d(k), &
1611 if (k.lt.kte-2 .and. k.gt.kts+1)
then
1612 write(*,*)
' below and above are: ', qv(i,k-1,j), qv(i,k+1,j)
1613 qv(i,k,j) = max(1.e-7, 0.5*(qv(i,k-1,j) + qv(i,k+1,j)))
1620 assign_extended_diagnostics:
if (ext_diag)
then
1625 prw_vcdc(i,k,j) = prw_vcdc(i,k,j) + prw_vcdc1(k)
1626 prw_vcde(i,k,j) = prw_vcde(i,k,j) + prw_vcde1(k)
1627 tpri_inu(i,k,j) = tpri_inu(i,k,j) + tpri_inu1(k)
1628 tpri_ide_d(i,k,j) = tpri_ide_d(i,k,j) + tpri_ide1_d(k)
1629 tpri_ide_s(i,k,j) = tpri_ide_s(i,k,j) + tpri_ide1_s(k)
1630 tprs_ide(i,k,j) = tprs_ide(i,k,j) + tprs_ide1(k)
1631 tprs_sde_s(i,k,j) = tprs_sde_s(i,k,j) + tprs_sde1_s(k)
1632 tprs_sde_d(i,k,j) = tprs_sde_d(i,k,j) + tprs_sde1_d(k)
1633 tprg_gde_d(i,k,j) = tprg_gde_d(i,k,j) + tprg_gde1_d(k)
1634 tprg_gde_s(i,k,j) = tprg_gde_s(i,k,j) + tprg_gde1_s(k)
1635 tpri_iha(i,k,j) = tpri_iha(i,k,j) + tpri_iha1(k)
1636 tpri_wfz(i,k,j) = tpri_wfz(i,k,j) + tpri_wfz1(k)
1637 tpri_rfz(i,k,j) = tpri_rfz(i,k,j) + tpri_rfz1(k)
1638 tprg_rfz(i,k,j) = tprg_rfz(i,k,j) + tprg_rfz1(k)
1639 tprs_scw(i,k,j) = tprs_scw(i,k,j) + tprs_scw1(k)
1640 tprg_scw(i,k,j) = tprg_scw(i,k,j) + tprg_scw1(k)
1641 tprg_rcs(i,k,j) = tprg_rcs(i,k,j) + tprg_rcs1(k)
1642 tprs_rcs(i,k,j) = tprs_rcs(i,k,j) + tprs_rcs1(k)
1643 tprr_rci(i,k,j) = tprr_rci(i,k,j) + tprr_rci1(k)
1644 tprg_rcg(i,k,j) = tprg_rcg(i,k,j) + tprg_rcg1(k)
1645 tprw_vcd_c(i,k,j) = tprw_vcd_c(i,k,j) + tprw_vcd1_c(k)
1646 tprw_vcd_e(i,k,j) = tprw_vcd_e(i,k,j) + tprw_vcd1_e(k)
1647 tprr_sml(i,k,j) = tprr_sml(i,k,j) + tprr_sml1(k)
1648 tprr_gml(i,k,j) = tprr_gml(i,k,j) + tprr_gml1(k)
1649 tprr_rcg(i,k,j) = tprr_rcg(i,k,j) + tprr_rcg1(k)
1650 tprr_rcs(i,k,j) = tprr_rcs(i,k,j) + tprr_rcs1(k)
1651 tprv_rev(i,k,j) = tprv_rev(i,k,j) + tprv_rev1(k)
1652 tten3(i,k,j) = tten3(i,k,j) + tten1(k)
1653 qvten3(i,k,j) = qvten3(i,k,j) + qvten1(k)
1654 qrten3(i,k,j) = qrten3(i,k,j) + qrten1(k)
1655 qsten3(i,k,j) = qsten3(i,k,j) + qsten1(k)
1656 qgten3(i,k,j) = qgten3(i,k,j) + qgten1(k)
1657 qiten3(i,k,j) = qiten3(i,k,j) + qiten1(k)
1658 niten3(i,k,j) = niten3(i,k,j) + niten1(k)
1659 nrten3(i,k,j) = nrten3(i,k,j) + nrten1(k)
1660 ncten3(i,k,j) = ncten3(i,k,j) + ncten1(k)
1661 qcten3(i,k,j) = qcten3(i,k,j) + qcten1(k)
1664 endif assign_extended_diagnostics
1666 if (ndt>1 .and. it==ndt)
then
1668 sr(i,j) = (pcp_sn(i,j) + pcp_gr(i,j) + pcp_ic(i,j))/(rainnc(i,j)+1.e-12)
1669 rainncv(i,j) = rainnc(i,j)
1670 IF (
PRESENT (snowncv) )
THEN
1671 snowncv(i,j) = snownc(i,j)
1673 IF (
PRESENT (icencv) )
THEN
1674 icencv(i,j) = icenc(i,j)
1676 IF (
PRESENT (graupelncv) )
THEN
1677 graupelncv(i,j) = graupelnc(i,j)
1683 last_step_only:
IF ((ndt>1 .and. it==ndt) .or. &
1684 (nsteps>1 .and. istep==nsteps) .or. &
1685 (nsteps==1 .and. ndt==1))
THEN
1691 diagflag_present:
IF (
PRESENT (diagflag) )
THEN
1692 if (diagflag .and. do_radar_ref == 1)
then
1695 if (fullradar_diag)
then
1701 if (
present(vt_dbz_wt))
then
1703 t1d, p1d, dbz, rand1, kts, kte, i, j, &
1704 melti, vt_dbz_wt(i,:,j), &
1708 t1d, p1d, dbz, rand1, kts, kte, i, j, &
1712 refl_10cm(i,k,j) = max(-35., dbz(k))
1715 ENDIF diagflag_present
1717 IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0)
THEN
1719 re_qc1d(k) = re_qc_min
1720 re_qi1d(k) = re_qi_min
1721 re_qs1d(k) = re_qs_min
1724 call calc_effectrad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
1725 re_qc1d, re_qi1d, re_qs1d, lsml, kts, kte)
1727 re_cloud(i,k,j) = max(re_qc_min, min(re_qc1d(k), re_qc_max))
1728 re_ice(i,k,j) = max(re_qi_min, min(re_qi1d(k), re_qi_max))
1729 re_snow(i,k,j) = max(re_qs_min, min(re_qs1d(k), re_qs_max))
1732 ENDIF last_step_only
1749 do j = j_start, j_end
1751 do i = i_start, i_end
1752 pfils(i,k,j) = pfils(i,k,j)/dt_in
1753 pflls(i,k,j) = pflls(i,k,j)/dt_in
1762 deallocate_extended_diagnostics:
if (ext_diag)
then
1763 deallocate (prw_vcdc1)
1764 deallocate (prw_vcde1)
1765 deallocate (tpri_inu1)
1766 deallocate (tpri_ide1_d)
1767 deallocate (tpri_ide1_s)
1768 deallocate (tprs_ide1)
1769 deallocate (tprs_sde1_d)
1770 deallocate (tprs_sde1_s)
1771 deallocate (tprg_gde1_d)
1772 deallocate (tprg_gde1_s)
1773 deallocate (tpri_iha1)
1774 deallocate (tpri_wfz1)
1775 deallocate (tpri_rfz1)
1776 deallocate (tprg_rfz1)
1777 deallocate (tprs_scw1)
1778 deallocate (tprg_scw1)
1779 deallocate (tprg_rcs1)
1780 deallocate (tprs_rcs1)
1781 deallocate (tprr_rci1)
1782 deallocate (tprg_rcg1)
1783 deallocate (tprw_vcd1_c)
1784 deallocate (tprw_vcd1_e)
1785 deallocate (tprr_sml1)
1786 deallocate (tprr_gml1)
1787 deallocate (tprr_rcg1)
1788 deallocate (tprr_rcs1)
1789 deallocate (tprv_rev1)
1800 end if deallocate_extended_diagnostics
1867 nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dzq, &
1868 lsml, pptrain, pptsnow, pptgraul, pptice, &
1869#if ( WRF_CHEM == 1 )
1870 rainprod, evapprod, &
1872 rand1, rand2, rand3, &
1873 kts, kte, dt, ii, jj, &
1879 prw_vcdc1, prw_vcde1, &
1880 tpri_inu1, tpri_ide1_d, tpri_ide1_s, tprs_ide1, &
1881 tprs_sde1_d, tprs_sde1_s, &
1882 tprg_gde1_d, tprg_gde1_s, tpri_iha1, tpri_wfz1, &
1883 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1, &
1884 tprg_rcs1, tprs_rcs1, tprr_rci1, &
1885 tprg_rcg1, tprw_vcd1_c, &
1886 tprw_vcd1_e, tprr_sml1, tprr_gml1, tprr_rcg1, &
1887 tprr_rcs1, tprv_rev1, &
1888 tten1, qvten1, qrten1, qsten1, &
1889 qgten1, qiten1, niten1, nrten1, ncten1, qcten1, &
1898 INTEGER,
INTENT(IN):: kts, kte, ii, jj
1899 REAL,
DIMENSION(kts:kte),
INTENT(INOUT):: &
1900 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
1901 nr1d, nc1d, nwfa1d, nifa1d, t1d
1902 REAL,
DIMENSION(kts:kte),
INTENT(OUT):: pfil1, pfll1
1903 REAL,
DIMENSION(kts:kte),
INTENT(IN):: p1d, w1d, dzq
1904 REAL,
INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice
1905 REAL,
INTENT(IN):: dt
1906 INTEGER,
INTENT(IN):: lsml
1907 REAL,
INTENT(IN):: rand1, rand2, rand3
1909 LOGICAL,
INTENT(IN) :: ext_diag
1910 LOGICAL,
INTENT(IN) :: sedi_semi
1911 INTEGER,
INTENT(IN) :: decfl
1912 REAL,
DIMENSION(:),
INTENT(OUT),
OPTIONAL :: &
1915 prw_vcde1, tpri_inu1, tpri_ide1_d, &
1916 tpri_ide1_s, tprs_ide1, &
1917 tprs_sde1_d, tprs_sde1_s, tprg_gde1_d, &
1918 tprg_gde1_s, tpri_iha1, tpri_wfz1, &
1919 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1,&
1920 tprg_rcs1, tprs_rcs1, &
1921 tprr_rci1, tprg_rcg1, &
1922 tprw_vcd1_c, tprw_vcd1_e, tprr_sml1, &
1923 tprr_gml1, tprr_rcg1, &
1924 tprr_rcs1, tprv_rev1, tten1, qvten1, &
1925 qrten1, qsten1, qgten1, qiten1, niten1, &
1926 nrten1, ncten1, qcten1
1928#if ( WRF_CHEM == 1 )
1929 REAL,
DIMENSION(kts:kte),
INTENT(INOUT):: &
1934 REAL,
DIMENSION(kts:kte):: tten, qvten, qcten, qiten, &
1935 qrten, qsten, qgten, niten, nrten, ncten, nwfaten, nifaten
1937 DOUBLE PRECISION,
DIMENSION(kts:kte):: prw_vcd
1939 DOUBLE PRECISION,
DIMENSION(kts:kte):: pnc_wcd, pnc_wau, pnc_rcw, &
1942 DOUBLE PRECISION,
DIMENSION(kts:kte):: pna_rca, pna_sca, pna_gca, &
1943 pnd_rcd, pnd_scd, pnd_gcd
1945 DOUBLE PRECISION,
DIMENSION(kts:kte):: prr_wau, prr_rcw, prr_rcs, &
1946 prr_rcg, prr_sml, prr_gml, &
1948 pnr_wau, pnr_rcs, pnr_rcg, &
1949 pnr_rci, pnr_sml, pnr_gml, &
1950 pnr_rev, pnr_rcr, pnr_rfz
1952 DOUBLE PRECISION,
DIMENSION(kts:kte):: pri_inu, pni_inu, pri_ihm, &
1953 pni_ihm, pri_wfz, pni_wfz, &
1954 pri_rfz, pni_rfz, pri_ide, &
1955 pni_ide, pri_rci, pni_rci, &
1956 pni_sci, pni_iau, pri_iha, pni_iha
1958 DOUBLE PRECISION,
DIMENSION(kts:kte):: prs_iau, prs_sci, prs_rcs, &
1959 prs_scw, prs_sde, prs_ihm, &
1962 DOUBLE PRECISION,
DIMENSION(kts:kte):: prg_scw, prg_rfz, prg_gde, &
1963 prg_gcw, prg_rci, prg_rcs, &
1966 DOUBLE PRECISION,
PARAMETER:: zeroD0 = 0.0d0
1967 REAL :: dtcfl,rainsfc,graulsfc
1970 REAL,
DIMENSION(kts:kte):: temp, pres, qv, pfll, pfil, pdummy
1971 REAL,
DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni, nr, nc, nwfa, nifa
1972 REAL,
DIMENSION(kts:kte):: rr_tmp, nr_tmp, rg_tmp
1973 REAL,
DIMENSION(kts:kte):: rho, rhof, rhof2
1974 REAL,
DIMENSION(kts:kte):: qvs, qvsi, delQvs
1975 REAL,
DIMENSION(kts:kte):: satw, sati, ssatw, ssati
1976 REAL,
DIMENSION(kts:kte):: diffu, visco, vsc2, &
1977 tcond, lvap, ocp, lvt2
1979 DOUBLE PRECISION,
DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g
1980 REAL,
DIMENSION(kts:kte):: mvd_r, mvd_c
1981 REAL,
DIMENSION(kts:kte):: smob, smo2, smo1, smo0, &
1982 smoc, smod, smoe, smof
1984 REAL,
DIMENSION(kts:kte):: sed_r, sed_s, sed_g, sed_i, sed_n,sed_c
1986 REAL:: rgvm, delta_tp, orho, lfus2, orhodt
1987 REAL,
DIMENSION(5):: onstep
1988 DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamc, lamr, lamg
1989 DOUBLE PRECISION:: lami, ilami, ilamc
1990 REAL:: xDc, Dc_b, Dc_g, xDi, xDr, xDs, xDg, Ds_m, Dg_m
1991 DOUBLE PRECISION:: Dr_star, Dc_star
1992 REAL:: zeta1, zeta, taud, tau
1993 REAL:: stoke_r, stoke_s, stoke_g, stoke_i
1994 REAL:: vti, vtr, vts, vtg, vtc
1995 REAL,
DIMENSION(kts:kte+1):: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk, &
1997 REAL,
DIMENSION(kts:kte):: vts_boost
1998 REAL:: Mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, C_snow
1999 REAL:: a_, b_, loga_, A1, A2, tf
2000 REAL:: tempc, tc0, r_mvd1, r_mvd2, xkrat
2001 REAL:: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr
2002 REAL:: xsat, rate_max, sump, ratio
2003 REAL:: clap, fcd, dfcd
2004 REAL:: otemp, rvs, rvs_p, rvs_pp, gamsc, alphsc, t1_evap, t1_subl
2005 REAL:: r_frac, g_frac
2006 REAL:: Ef_rw, Ef_sw, Ef_gw, Ef_rr
2007 REAL:: Ef_ra, Ef_sa, Ef_ga
2008 REAL:: dtsave, odts, odt, odzq, hgt_agl, SR
2009 REAL:: xslw1, ygra1, zans1, eva_factor
2011 INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq
2012 INTEGER,
DIMENSION(5):: ksed1
2013 INTEGER:: nir, nis, nig, nii, nic, niin
2014 INTEGER:: idx_tc, idx_t, idx_s, idx_g1, idx_g, idx_r1, idx_r, &
2015 idx_i1, idx_i, idx_c, idx, idx_d, idx_n, idx_in
2018 LOGICAL,
DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg
2019 LOGICAL:: debug_flag
2024 debug_flag = .false.
2027 write(*, *)
'DEBUG INFO, mp_thompson at (i,j) ', ii,
', ', jj
2036 av_i = av_s * d0s ** (bv_s - bv_i)
2141#if ( WRF_CHEM == 1 )
2211 qv(k) = max(1.e-10, qv1d(k))
2213 rho(k) = 0.622*pres(k)/(r*temp(k)*(qv(k)+0.622))
2214 nwfa(k) = max(11.1e6*rho(k), min(9999.e6*rho(k), nwfa1d(k)*rho(k)))
2215 nifa(k) = max(nain1*0.01*rho(k), min(9999.e6*rho(k), nifa1d(k)*rho(k)))
2219 if (qc1d(k) .gt. r1)
then
2221 rc(k) = qc1d(k)*rho(k)
2222 nc(k) = max(2., min(nc1d(k)*rho(k), nt_c_max))
2224 if (nc(k).gt.10000.e6)
then
2226 elseif (nc(k).lt.100.)
then
2229 nu_c = nint(1000.e6/nc(k)) + 2
2230 nu_c = max(2, min(nu_c+nint(rand2), 15))
2232 lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
2233 xdc = (bm_r + nu_c + 1.) / lamc
2234 if (xdc.lt. d0c)
then
2235 lamc = cce(2,nu_c)/d0c
2236 elseif (xdc.gt. d0r*2.)
then
2237 lamc = cce(2,nu_c)/(d0r*2.)
2239 nc(k) = min( dble(nt_c_max), ccg(1,nu_c)*ocg2(nu_c)*rc(k) &
2241 if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware))
then
2256 if (qi1d(k) .gt. r1)
then
2258 ri(k) = qi1d(k)*rho(k)
2259 ni(k) = max(r2, ni1d(k)*rho(k))
2260 if (ni(k).le. r2)
then
2262 ni(k) = min(4999.d3, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
2265 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
2267 xdi = (bm_i + mu_i + 1.) * ilami
2268 if (xdi.lt. 5.e-6)
then
2270 ni(k) = min(4999.d3, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
2271 elseif (xdi.gt. 300.e-6)
then
2272 lami = cie(2)/300.e-6
2273 ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i
2283 if (qr1d(k) .gt. r1)
then
2285 rr(k) = qr1d(k)*rho(k)
2286 nr(k) = max(r2, nr1d(k)*rho(k))
2287 if (nr(k).le. r2)
then
2289 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
2290 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
2293 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
2294 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
2295 if (mvd_r(k) .gt. 2.5e-3)
then
2297 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
2298 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
2299 elseif (mvd_r(k) .lt. d0r*0.75)
then
2301 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
2302 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
2311 if (qs1d(k) .gt. r1)
then
2313 rs(k) = qs1d(k)*rho(k)
2320 if (qg1d(k) .gt. r1)
then
2322 rg(k) = qg1d(k)*rho(k)
2347 tempc = temp(k) - 273.15
2348 rhof(k) = sqrt(rho_not/rho(k))
2349 rhof2(k) = sqrt(rhof(k))
2350 qvs(k) =
rslf(pres(k), temp(k))
2351 delqvs(k) = max(0.0,
rslf(pres(k), 273.15)-qv(k))
2352 if (tempc .le. 0.0)
then
2353 qvsi(k) =
rsif(pres(k), temp(k))
2357 satw(k) = qv(k)/qvs(k)
2358 sati(k) = qv(k)/qvsi(k)
2359 ssatw(k) = satw(k) - 1.
2360 ssati(k) = sati(k) - 1.
2361 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
2362 if (abs(ssati(k)).lt. eps) ssati(k) = 0.0
2363 if (no_micro .and. ssati(k).gt. 0.0) no_micro = .false.
2364 diffu(k) = 2.11e-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
2365 if (tempc .ge. 0.0)
then
2366 visco(k) = (1.718+0.0049*tempc)*1.0e-5
2368 visco(k) = (1.718+0.0049*tempc-1.2e-5*tempc*tempc)*1.0e-5
2370 ocp(k) = 1./(cp*(1.+0.887*qv(k)))
2371 vsc2(k) = sqrt(rho(k)/visco(k))
2372 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
2373 tcond(k) = (5.69 + 0.0168*tempc)*1.0e-5 * 418.936
2381 if (no_micro)
return
2386 if (.not. iiwarm)
then
2388 if (.not. l_qs(k)) cycle
2389 tc0 = min(-0.1, temp(k)-273.15)
2390 smob(k) = rs(k)*oams
2394 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3))
then
2397 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
2398 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
2399 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
2400 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
2401 + sa(10)*bm_s*bm_s*bm_s
2403 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
2404 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
2405 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
2406 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
2407 + sb(10)*bm_s*bm_s*bm_s
2408 smo2(k) = (smob(k)/a_)**(1./b_)
2412 loga_ = sa(1) + sa(2)*tc0 + sa(5)*tc0*tc0 + sa(9)*tc0*tc0*tc0
2414 b_ = sb(1) + sb(2)*tc0 + sb(5)*tc0*tc0 + sb(9)*tc0*tc0*tc0
2415 smo0(k) = a_ * smo2(k)**b_
2418 loga_ = sa(1) + sa(2)*tc0 + sa(3) &
2419 + sa(4)*tc0 + sa(5)*tc0*tc0 &
2420 + sa(6) + sa(7)*tc0*tc0 &
2421 + sa(8)*tc0 + sa(9)*tc0*tc0*tc0 &
2424 b_ = sb(1)+ sb(2)*tc0 + sb(3) + sb(4)*tc0 &
2425 + sb(5)*tc0*tc0 + sb(6) &
2426 + sb(7)*tc0*tc0 + sb(8)*tc0 &
2427 + sb(9)*tc0*tc0*tc0 + sb(10)
2428 smo1(k) = a_ * smo2(k)**b_
2431 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
2432 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
2433 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
2434 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
2435 + sa(10)*cse(1)*cse(1)*cse(1)
2437 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
2438 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
2439 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
2440 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
2441 smoc(k) = a_ * smo2(k)**b_
2444 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(13) &
2445 + sa(4)*tc0*cse(13) + sa(5)*tc0*tc0 &
2446 + sa(6)*cse(13)*cse(13) + sa(7)*tc0*tc0*cse(13) &
2447 + sa(8)*tc0*cse(13)*cse(13) + sa(9)*tc0*tc0*tc0 &
2448 + sa(10)*cse(13)*cse(13)*cse(13)
2450 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(13) + sb(4)*tc0*cse(13) &
2451 + sb(5)*tc0*tc0 + sb(6)*cse(13)*cse(13) &
2452 + sb(7)*tc0*tc0*cse(13) + sb(8)*tc0*cse(13)*cse(13) &
2453 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(13)*cse(13)*cse(13)
2454 smoe(k) = a_ * smo2(k)**b_
2457 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(16) &
2458 + sa(4)*tc0*cse(16) + sa(5)*tc0*tc0 &
2459 + sa(6)*cse(16)*cse(16) + sa(7)*tc0*tc0*cse(16) &
2460 + sa(8)*tc0*cse(16)*cse(16) + sa(9)*tc0*tc0*tc0 &
2461 + sa(10)*cse(16)*cse(16)*cse(16)
2463 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(16) + sb(4)*tc0*cse(16) &
2464 + sb(5)*tc0*tc0 + sb(6)*cse(16)*cse(16) &
2465 + sb(7)*tc0*tc0*cse(16) + sb(8)*tc0*cse(16)*cse(16) &
2466 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(16)*cse(16)*cse(16)
2467 smof(k) = a_ * smo2(k)**b_
2481 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
2483 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
2484 n0_r(k) = nr(k)*org2*lamr**cre(2)
2495 if (l_qr(k) .and. mvd_r(k).gt. d0r)
then
2496 ef_rr = max(-0.1, 1.0 - exp(2300.0*(mvd_r(k)-1950.0e-6)))
2497 pnr_rcr(k) = ef_rr * 2.0*nr(k)*rr(k)
2501 if (nc(k).gt.10000.e6)
then
2503 elseif (nc(k).lt.100.)
then
2506 nu_c = nint(1000.e6/nc(k)) + 2
2507 nu_c = max(2, min(nu_c+nint(rand2), 15))
2509 xdc = max(d0c*1.e6, ((rc(k)/(am_r*nc(k)))**obmr) * 1.e6)
2510 lamc = (nc(k)*am_r* ccg(2,nu_c) * ocg1(nu_c) / rc(k))**obmr
2511 mvd_c(k) = (3.0+nu_c+0.672) / lamc
2512 mvd_c(k) = max(d0c, min(mvd_c(k), d0r))
2517 if (rc(k).gt. 0.01e-3)
then
2518 dc_g = ((ccg(3,nu_c)*ocg2(nu_c))**obmr / lamc) * 1.e6
2519 dc_b = (xdc*xdc*xdc*dc_g*dc_g*dc_g - xdc*xdc*xdc*xdc*xdc*xdc) &
2521 zeta1 = 0.5*((6.25e-6*xdc*dc_b*dc_b*dc_b - 0.4) &
2522 + abs(6.25e-6*xdc*dc_b*dc_b*dc_b - 0.4))
2523 zeta = 0.027*rc(k)*zeta1
2524 taud = 0.5*((0.5*dc_b - 7.5) + abs(0.5*dc_b - 7.5)) + r1
2525 tau = 3.72/(rc(k)*taud)
2526 prr_wau(k) = zeta/tau
2527 prr_wau(k) = min(dble(rc(k)*odts), prr_wau(k))
2528 pnr_wau(k) = prr_wau(k) / (am_r*nu_c*10.*d0r*d0r*d0r)
2529 pnc_wau(k) = min(dble(nc(k)*odts), prr_wau(k) &
2530 / (am_r*mvd_c(k)*mvd_c(k)*mvd_c(k)))
2534 if (l_qr(k) .and. mvd_r(k).gt. d0r .and. mvd_c(k).gt. d0c)
then
2536 idx = 1 + int(nbr*dlog(mvd_r(k)/dr(1))/dlog(dr(nbr)/dr(1)))
2538 ef_rw = t_efrw(idx, int(mvd_c(k)*1.e6))
2539 prr_rcw(k) = rhof(k)*t1_qr_qc*ef_rw*rc(k)*n0_r(k) &
2540 *((lamr+fv_r)**(-cre(9)))
2541 prr_rcw(k) = min(dble(rc(k)*odts), prr_rcw(k))
2542 pnc_rcw(k) = rhof(k)*t1_qr_qc*ef_rw*nc(k)*n0_r(k) &
2543 *((lamr+fv_r)**(-cre(9)))
2544 pnc_rcw(k) = min(dble(nc(k)*odts), pnc_rcw(k))
2548 if (l_qr(k) .and. mvd_r(k).gt. d0r)
then
2549 ef_ra =
eff_aero(mvd_r(k),0.04e-6,visco(k),rho(k),temp(k),
'r')
2551 pna_rca(k) = rhof(k)*t1_qr_qc*ef_ra*nwfa(k)*n0_r(k) &
2552 *((lamr+fv_r)**(-cre(9)))
2553 pna_rca(k) = min(dble(nwfa(k)*odts), pna_rca(k))
2555 ef_ra =
eff_aero(mvd_r(k),0.8e-6,visco(k),rho(k),temp(k),
'r')
2556 pnd_rcd(k) = rhof(k)*t1_qr_qc*ef_ra*nifa(k)*n0_r(k) &
2557 *((lamr+fv_r)**(-cre(9)))
2558 pnd_rcd(k) = min(dble(nifa(k)*odts), pnd_rcd(k))
2566 if (.not. iiwarm)
then
2570 if (l_qs(k)) xds = smoc(k) / smob(k)
2573 tempc = temp(k) - 273.15
2574 idx_tc = max(1, min(nint(-tempc), 45) )
2575 idx_t = int( (tempc-2.5)/5. ) - 1
2576 idx_t = max(1, -idx_t)
2577 idx_t = min(idx_t, ntb_t)
2578 it = max(1, min(nint(-tempc), 31) )
2581 if (rc(k).gt. r_c(1))
then
2582 nic = nint(alog10(rc(k)))
2583 do nn = nic-1, nic+1
2585 if ( (rc(k)/10.**nn).ge.1.0 .and. &
2586 (rc(k)/10.**nn).lt.10.0)
goto 141
2589 idx_c = int(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
2590 idx_c = max(1, min(idx_c, ntb_c))
2596 idx_n = nint(1.0 + float(nbc) * dlog(nc(k)/t_nc(1)) / nic1)
2597 idx_n = max(1, min(idx_n, nbc))
2600 if (ri(k).gt. r_i(1))
then
2601 nii = nint(alog10(ri(k)))
2602 do nn = nii-1, nii+1
2604 if ( (ri(k)/10.**nn).ge.1.0 .and. &
2605 (ri(k)/10.**nn).lt.10.0)
goto 142
2608 idx_i = int(ri(k)/10.**n) + 10*(n-nii2) - (n-nii2)
2609 idx_i = max(1, min(idx_i, ntb_i))
2614 if (ni(k).gt. nt_i(1))
then
2615 nii = nint(alog10(ni(k)))
2616 do nn = nii-1, nii+1
2618 if ( (ni(k)/10.**nn).ge.1.0 .and. &
2619 (ni(k)/10.**nn).lt.10.0)
goto 143
2622 idx_i1 = int(ni(k)/10.**n) + 10*(n-nii3) - (n-nii3)
2623 idx_i1 = max(1, min(idx_i1, ntb_i1))
2629 if (rr(k).gt. r_r(1))
then
2630 nir = nint(alog10(rr(k)))
2631 do nn = nir-1, nir+1
2633 if ( (rr(k)/10.**nn).ge.1.0 .and. &
2634 (rr(k)/10.**nn).lt.10.0)
goto 144
2637 idx_r = int(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2)
2638 idx_r = max(1, min(idx_r, ntb_r))
2641 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
2642 n0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
2643 nir = nint(dlog10(n0_exp))
2644 do nn = nir-1, nir+1
2646 if ( (n0_exp/10.**nn).ge.1.0 .and. &
2647 (n0_exp/10.**nn).lt.10.0)
goto 145
2650 idx_r1 = int(n0_exp/10.**n) + 10*(n-nir3) - (n-nir3)
2651 idx_r1 = max(1, min(idx_r1, ntb_r1))
2658 if (rs(k).gt. r_s(1))
then
2659 nis = nint(alog10(rs(k)))
2660 do nn = nis-1, nis+1
2662 if ( (rs(k)/10.**nn).ge.1.0 .and. &
2663 (rs(k)/10.**nn).lt.10.0)
goto 146
2666 idx_s = int(rs(k)/10.**n) + 10*(n-nis2) - (n-nis2)
2667 idx_s = max(1, min(idx_s, ntb_s))
2673 if (rg(k).gt. r_g(1))
then
2674 nig = nint(alog10(rg(k)))
2675 do nn = nig-1, nig+1
2677 if ( (rg(k)/10.**nn).ge.1.0 .and. &
2678 (rg(k)/10.**nn).lt.10.0)
goto 147
2681 idx_g = int(rg(k)/10.**n) + 10*(n-nig2) - (n-nig2)
2682 idx_g = max(1, min(idx_g, ntb_g))
2685 lam_exp = lamg * (cgg(3)*ogg2*ogg1)**bm_g
2686 n0_exp = ogg1*rg(k)/am_g * lam_exp**cge(1)
2687 nig = nint(dlog10(n0_exp))
2688 do nn = nig-1, nig+1
2690 if ( (n0_exp/10.**nn).ge.1.0 .and. &
2691 (n0_exp/10.**nn).lt.10.0)
goto 148
2694 idx_g1 = int(n0_exp/10.**n) + 10*(n-nig3) - (n-nig3)
2695 idx_g1 = max(1, min(idx_g1, ntb_g1))
2703 rvs = rho(k)*qvsi(k)
2704 rvs_p = rvs*otemp*(lsub*otemp*orv - 1.)
2705 rvs_pp = rvs * ( otemp*(lsub*otemp*orv - 1.) &
2706 *otemp*(lsub*otemp*orv - 1.) &
2707 + (-2.*lsub*otemp*otemp*otemp*orv) &
2709 gamsc = lsub*diffu(k)/tcond(k) * rvs_p
2710 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
2711 * rvs_pp/rvs_p * rvs/rvs_p
2712 alphsc = max(1.e-9, alphsc)
2714 if (abs(xsat).lt. 1.e-9) xsat=0.
2715 t1_subl = 4.*pi*( 1.0 - alphsc*xsat &
2716 + 2.*alphsc*alphsc*xsat*xsat &
2717 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
2721 if (l_qc(k) .and. mvd_c(k).gt. d0c)
then
2722 if (xds .gt. d0s)
then
2723 idx = 1 + int(nbs*dlog(xds/ds(1))/dlog(ds(nbs)/ds(1)))
2725 ef_sw = t_efsw(idx, int(mvd_c(k)*1.e6))
2726 prs_scw(k) = rhof(k)*t1_qs_qc*ef_sw*rc(k)*smoe(k)
2727 prs_scw(k) = min(dble(rc(k)*odts), prs_scw(k))
2728 pnc_scw(k) = rhof(k)*t1_qs_qc*ef_sw*nc(k)*smoe(k)
2729 pnc_scw(k) = min(dble(nc(k)*odts), pnc_scw(k))
2733 if (rg(k).ge. r_g(1) .and. mvd_c(k).gt. d0c)
then
2734 xdg = (bm_g + mu_g + 1.) * ilamg(k)
2735 vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
2736 stoke_g = mvd_c(k)*mvd_c(k)*vtg*rho_w/(9.*visco(k)*xdg)
2737 if (xdg.gt. d0g)
then
2738 if (stoke_g.ge.0.4 .and. stoke_g.le.10.)
then
2739 ef_gw = 0.55*alog10(2.51*stoke_g)
2740 elseif (stoke_g.lt.0.4)
then
2742 elseif (stoke_g.gt.10)
then
2745 prg_gcw(k) = rhof(k)*t1_qg_qc*ef_gw*rc(k)*n0_g(k) &
2747 pnc_gcw(k) = rhof(k)*t1_qg_qc*ef_gw*nc(k)*n0_g(k) &
2749 pnc_gcw(k) = min(dble(nc(k)*odts), pnc_gcw(k))
2755 if (rs(k) .gt. r_s(1))
then
2756 ef_sa =
eff_aero(xds,0.04e-6,visco(k),rho(k),temp(k),
's')
2757 pna_sca(k) = rhof(k)*t1_qs_qc*ef_sa*nwfa(k)*smoe(k)
2758 pna_sca(k) = min(dble(nwfa(k)*odts), pna_sca(k))
2760 ef_sa =
eff_aero(xds,0.8e-6,visco(k),rho(k),temp(k),
's')
2761 pnd_scd(k) = rhof(k)*t1_qs_qc*ef_sa*nifa(k)*smoe(k)
2762 pnd_scd(k) = min(dble(nifa(k)*odts), pnd_scd(k))
2764 if (rg(k) .gt. r_g(1))
then
2765 xdg = (bm_g + mu_g + 1.) * ilamg(k)
2766 ef_ga =
eff_aero(xdg,0.04e-6,visco(k),rho(k),temp(k),
'g')
2767 pna_gca(k) = rhof(k)*t1_qg_qc*ef_ga*nwfa(k)*n0_g(k) &
2769 pna_gca(k) = min(dble(nwfa(k)*odts), pna_gca(k))
2771 ef_ga =
eff_aero(xdg,0.8e-6,visco(k),rho(k),temp(k),
'g')
2772 pnd_gcd(k) = rhof(k)*t1_qg_qc*ef_ga*nifa(k)*n0_g(k) &
2774 pnd_gcd(k) = min(dble(nifa(k)*odts), pnd_gcd(k))
2780 if (rr(k).ge. r_r(1))
then
2781 if (rs(k).ge. r_s(1))
then
2782 if (temp(k).lt.t_0)
then
2783 prr_rcs(k) = -(tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2784 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
2785 + tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
2786 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r))
2787 prs_rcs(k) = tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2788 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
2789 - tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2790 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
2791 prg_rcs(k) = tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
2792 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2793 + tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2794 + tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
2795 prr_rcs(k) = max(dble(-rr(k)*odts), prr_rcs(k))
2796 prs_rcs(k) = max(dble(-rs(k)*odts), prs_rcs(k))
2797 prg_rcs(k) = min(dble((rr(k)+rs(k))*odts), prg_rcs(k))
2798 pnr_rcs(k) = tnr_racs1(idx_s,idx_t,idx_r1,idx_r) &
2799 + tnr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2800 + tnr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2801 + tnr_sacr2(idx_s,idx_t,idx_r1,idx_r)
2802 pnr_rcs(k) = min(dble(nr(k)*odts), pnr_rcs(k))
2804 prs_rcs(k) = -tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2805 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2806 + tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2807 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r)
2808 prs_rcs(k) = max(dble(-rs(k)*odts), prs_rcs(k))
2809 prr_rcs(k) = -prs_rcs(k)
2816 if (rg(k).ge. r_g(1))
then
2817 if (temp(k).lt.t_0)
then
2818 prg_rcg(k) = tmr_racg(idx_g1,idx_g,idx_r1,idx_r) &
2819 + tcr_gacr(idx_g1,idx_g,idx_r1,idx_r)
2820 prg_rcg(k) = min(dble(rr(k)*odts), prg_rcg(k))
2821 prr_rcg(k) = -prg_rcg(k)
2822 pnr_rcg(k) = tnr_racg(idx_g1,idx_g,idx_r1,idx_r) &
2823 + tnr_gacr(idx_g1,idx_g,idx_r1,idx_r)
2824 pnr_rcg(k) = min(dble(nr(k)*odts), pnr_rcg(k))
2826 prr_rcg(k) = tcg_racg(idx_g1,idx_g,idx_r1,idx_r)
2827 prr_rcg(k) = min(dble(rg(k)*odts), prr_rcg(k))
2828 prg_rcg(k) = -prr_rcg(k)
2830 pnr_rcg(k) = -1.5*tnr_gacr(idx_g1,idx_g,idx_r1,idx_r)
2835 if (temp(k).lt.t_0)
then
2836 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
2840 c_snow = c_sqrd + (tempc+1.5)*(c_cube-c_sqrd)/(-30.+1.5)
2841 c_snow = max(c_sqrd, min(c_snow, c_cube))
2842 prs_sde(k) = c_snow*t1_subl*diffu(k)*ssati(k)*rvs &
2843 * (t1_qs_sd*smo1(k) &
2844 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
2845 if (prs_sde(k).lt. 0.)
then
2846 prs_sde(k) = max(dble(-rs(k)*odts), prs_sde(k), dble(rate_max))
2848 prs_sde(k) = min(prs_sde(k), dble(rate_max))
2852 if (l_qg(k) .and. ssati(k).lt. -eps)
then
2853 prg_gde(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
2854 * n0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
2855 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
2856 if (prg_gde(k).lt. 0.)
then
2857 prg_gde(k) = max(dble(-rg(k)*odts), prg_gde(k), dble(rate_max))
2859 prg_gde(k) = min(prg_gde(k), dble(rate_max))
2867 if (prs_scw(k).gt.5.0*prs_sde(k) .and. &
2868 prs_sde(k).gt.eps)
then
2869 r_frac = min(30.0d0, prs_scw(k)/prs_sde(k))
2870 g_frac = min(0.75, 0.15 + (r_frac-5.)*.028)
2871 vts_boost(k) = min(1.5, 1.1 + (r_frac-5.)*.016)
2872 prg_scw(k) = g_frac*prs_scw(k)
2873 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(alog10(xni))
2909 do nn = niin-1, niin+1
2911 if ( (xni/10.**nn).ge.1.0 .and. &
2912 (xni/10.**nn).lt.10.0)
goto 149
2915 idx_in = int(xni/10.**n) + 10*(n-niin2) - (n-niin2)
2916 idx_in = max(1, min(idx_in, ntb_in))
2922 if (rr(k).gt. r_r(1))
then
2923 prg_rfz(k) = tpg_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
2924 pri_rfz(k) = tpi_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
2925 pni_rfz(k) = tni_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
2926 pnr_rfz(k) = tnr_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
2927 pnr_rfz(k) = min(dble(nr(k)*odts), pnr_rfz(k))
2928 elseif (rr(k).gt. r1 .and. temp(k).lt.hgfr)
then
2929 pri_rfz(k) = rr(k)*odts
2930 pni_rfz(k) = pnr_rfz(k)
2933 if (rc(k).gt. r_c(1))
then
2934 pri_wfz(k) = tpi_qcfz(idx_c,idx_n,idx_tc,idx_in)*odts
2935 pri_wfz(k) = min(dble(rc(k)*odts), pri_wfz(k))
2936 pni_wfz(k) = tni_qcfz(idx_c,idx_n,idx_tc,idx_in)*odts
2937 pni_wfz(k) = min(dble(nc(k)*odts), pri_wfz(k)/(2.*xm0i), &
2939 elseif (rc(k).gt. r1 .and. temp(k).lt.hgfr)
then
2940 pri_wfz(k) = rc(k)*odts
2941 pni_wfz(k) = nc(k)*odts
2946 if ( (ssati(k).ge. 0.15) .or. (ssatw(k).gt. eps &
2947 .and. temp(k).lt.253.15) )
then
2948 if (dustyice .AND. (is_aerosol_aware .or. merra2_aerosol_aware))
then
2949 xnc =
icedemott(tempc,qv(k),qvs(k),qvsi(k),rho(k),nifa(k))
2950 xnc = xnc*(1.0 + 50.*rand3)
2952 xnc = min(1000.e3, tno*exp(ato*(t_0-temp(k))))
2954 xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dtsave
2955 pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts
2956 pri_inu(k) = min(dble(rate_max), xm0i*pni_inu(k))
2957 pni_inu(k) = pri_inu(k)/xm0i
2961 xni = smo0(k)+ni(k) + (pni_rfz(k)+pni_wfz(k)+pni_inu(k))*dtsave
2962 if ((is_aerosol_aware .or. merra2_aerosol_aware) .AND. homogice .AND. (xni.le.4999.e3) &
2963 & .AND.(temp(k).lt.238).AND.(ssati(k).ge.0.4) )
then
2964 xnc =
icekoop(temp(k),qv(k),qvs(k),nwfa(k), dtsave)
2965 pni_iha(k) = xnc*odts
2966 pri_iha(k) = min(dble(rate_max), xm0i*0.1*pni_iha(k))
2967 pni_iha(k) = pri_iha(k)/(xm0i*0.1)
2974 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
2976 xdi = max(dble(d0i), (bm_i + mu_i + 1.) * ilami)
2977 xmi = am_i*xdi**bm_i
2979 pri_ide(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
2980 *oig1*cig(5)*ni(k)*ilami
2982 if (pri_ide(k) .lt. 0.0)
then
2983 pri_ide(k) = max(dble(-ri(k)*odts), pri_ide(k), dble(rate_max))
2984 pni_ide(k) = pri_ide(k)*oxmi
2985 pni_ide(k) = max(dble(-ni(k)*odts), pni_ide(k))
2987 pri_ide(k) = min(pri_ide(k), dble(rate_max))
2988 prs_ide(k) = (1.0d0-tpi_ide(idx_i,idx_i1))*pri_ide(k)
2989 pri_ide(k) = tpi_ide(idx_i,idx_i1)*pri_ide(k)
2994 if ( (idx_i.eq. ntb_i) .or. (xdi.gt. 5.0*d0s) )
then
2995 prs_iau(k) = ri(k)*.99*odts
2996 pni_iau(k) = ni(k)*.95*odts
2997 elseif (xdi.lt. 0.1*d0s)
then
3001 prs_iau(k) = tps_iaus(idx_i,idx_i1)*odts
3002 prs_iau(k) = min(dble(ri(k)*.99*odts), prs_iau(k))
3003 pni_iau(k) = tni_iaus(idx_i,idx_i1)*odts
3004 pni_iau(k) = min(dble(ni(k)*.95*odts), pni_iau(k))
3010 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
3012 xdi = max(dble(d0i), (bm_i + mu_i + 1.) * ilami)
3013 xmi = am_i*xdi**bm_i
3015 if (rs(k).ge. r_s(1))
then
3016 prs_sci(k) = t1_qs_qi*rhof(k)*ef_si*ri(k)*smoe(k)
3017 pni_sci(k) = prs_sci(k) * oxmi
3021 if (rr(k).ge. r_r(1) .and. mvd_r(k).gt. 4.*xdi)
then
3023 pri_rci(k) = rhof(k)*t1_qr_qi*ef_ri*ri(k)*n0_r(k) &
3024 *((lamr+fv_r)**(-cre(9)))
3025 pnr_rci(k) = rhof(k)*t1_qr_qi*ef_ri*ni(k)*n0_r(k) &
3026 *((lamr+fv_r)**(-cre(9)))
3027 pni_rci(k) = pri_rci(k) * oxmi
3028 prr_rci(k) = rhof(k)*t2_qr_qi*ef_ri*ni(k)*n0_r(k) &
3029 *((lamr+fv_r)**(-cre(8)))
3030 prr_rci(k) = min(dble(rr(k)*odts), prr_rci(k))
3031 prg_rci(k) = pri_rci(k) + prr_rci(k)
3036 if (prg_gcw(k).gt. eps .and. tempc.gt.-8.0)
then
3038 if (tempc.ge.-5.0 .and. tempc.lt.-3.0)
then
3039 tf = 0.5*(-3.0 - tempc)
3040 elseif (tempc.gt.-8.0 .and. tempc.lt.-5.0)
then
3041 tf = 0.33333333*(8.0 + tempc)
3043 pni_ihm(k) = 3.5e8*tf*prg_gcw(k)
3044 pri_ihm(k) = xm0i*pni_ihm(k)
3045 prs_ihm(k) = prs_scw(k)/(prs_scw(k)+prg_gcw(k)) &
3047 prg_ihm(k) = prg_gcw(k)/(prs_scw(k)+prg_gcw(k)) &
3056 prr_sml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delqvs(k)) &
3057 * (t1_qs_me*smo1(k) + t2_qs_me*rhof2(k)*vsc2(k)*smof(k))
3058 if (prr_sml(k) .gt. 0.)
then
3059 prr_sml(k) = prr_sml(k) + 4218.*olfus*tempc &
3060 * (prr_rcs(k)+prs_scw(k))
3061 prr_sml(k) = min(dble(rs(k)*odts), prr_sml(k))
3062 pnr_sml(k) = smo0(k)/rs(k)*prr_sml(k) * 10.0**(-0.25*tempc)
3063 pnr_sml(k) = min(dble(smo0(k)*odts), pnr_sml(k))
3064 elseif (ssati(k).lt. 0.)
then
3066 prs_sde(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
3067 * (t1_qs_sd*smo1(k) &
3068 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
3069 prs_sde(k) = max(dble(-rs(k)*odts), prs_sde(k))
3074 prr_gml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delqvs(k)) &
3075 * n0_g(k)*(t1_qg_me*ilamg(k)**cge(10) &
3076 + t2_qg_me*rhof2(k)*vsc2(k)*ilamg(k)**cge(11))
3077 if (prr_gml(k) .gt. 0.)
then
3078 prr_gml(k) = min(dble(rg(k)*odts), prr_gml(k))
3079 pnr_gml(k) = n0_g(k)*cgg(2)*ilamg(k)**cge(2) / rg(k) &
3080 * prr_gml(k) * 10.0**(-0.5*tempc)
3081 elseif (ssati(k).lt. 0.)
then
3083 prg_gde(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
3084 * n0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
3085 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
3086 prg_gde(k) = max(dble(-rg(k)*odts), prg_gde(k))
3094 if (dt .gt. 120.)
then
3095 prr_rcw(k)=prr_rcw(k)+prs_scw(k)+prg_gcw(k)
3114 sump = pri_inu(k) + pri_ide(k) + prs_ide(k) &
3115 + prs_sde(k) + prg_gde(k) + pri_iha(k)
3116 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
3117 if ( (sump.gt. eps .and. sump.gt. rate_max) .or. &
3118 (sump.lt. -eps .and. sump.lt. rate_max) )
then
3119 ratio = rate_max/sump
3120 pri_inu(k) = pri_inu(k) * ratio
3121 pri_ide(k) = pri_ide(k) * ratio
3122 pni_ide(k) = pni_ide(k) * ratio
3123 prs_ide(k) = prs_ide(k) * ratio
3124 prs_sde(k) = prs_sde(k) * ratio
3125 prg_gde(k) = prg_gde(k) * ratio
3126 pri_iha(k) = pri_iha(k) * ratio
3130 sump = -prr_wau(k) - pri_wfz(k) - prr_rcw(k) &
3131 - prs_scw(k) - prg_scw(k) - prg_gcw(k)
3132 rate_max = -rc(k)*odts
3133 if (sump.lt. rate_max .and. l_qc(k))
then
3134 ratio = rate_max/sump
3135 prr_wau(k) = prr_wau(k) * ratio
3136 pri_wfz(k) = pri_wfz(k) * ratio
3137 prr_rcw(k) = prr_rcw(k) * ratio
3138 prs_scw(k) = prs_scw(k) * ratio
3139 prg_scw(k) = prg_scw(k) * ratio
3140 prg_gcw(k) = prg_gcw(k) * ratio
3144 sump = pri_ide(k) - prs_iau(k) - prs_sci(k) &
3146 rate_max = -ri(k)*odts
3147 if (sump.lt. rate_max .and. l_qi(k))
then
3148 ratio = rate_max/sump
3149 pri_ide(k) = pri_ide(k) * ratio
3150 prs_iau(k) = prs_iau(k) * ratio
3151 prs_sci(k) = prs_sci(k) * ratio
3152 pri_rci(k) = pri_rci(k) * ratio
3156 sump = -prg_rfz(k) - pri_rfz(k) - prr_rci(k) &
3157 + prr_rcs(k) + prr_rcg(k)
3158 rate_max = -rr(k)*odts
3159 if (sump.lt. rate_max .and. l_qr(k))
then
3160 ratio = rate_max/sump
3161 prg_rfz(k) = prg_rfz(k) * ratio
3162 pri_rfz(k) = pri_rfz(k) * ratio
3163 prr_rci(k) = prr_rci(k) * ratio
3164 prr_rcs(k) = prr_rcs(k) * ratio
3165 prr_rcg(k) = prr_rcg(k) * ratio
3169 sump = prs_sde(k) - prs_ihm(k) - prr_sml(k) &
3171 rate_max = -rs(k)*odts
3172 if (sump.lt. rate_max .and. l_qs(k))
then
3173 ratio = rate_max/sump
3174 prs_sde(k) = prs_sde(k) * ratio
3175 prs_ihm(k) = prs_ihm(k) * ratio
3176 prr_sml(k) = prr_sml(k) * ratio
3177 prs_rcs(k) = prs_rcs(k) * ratio
3181 sump = prg_gde(k) - prg_ihm(k) - prr_gml(k) &
3183 rate_max = -rg(k)*odts
3184 if (sump.lt. rate_max .and. l_qg(k))
then
3185 ratio = rate_max/sump
3186 prg_gde(k) = prg_gde(k) * ratio
3187 prg_ihm(k) = prg_ihm(k) * ratio
3188 prr_gml(k) = prr_gml(k) * ratio
3189 prg_rcg(k) = prg_rcg(k) * ratio
3194 pri_ihm(k) = prs_ihm(k) + prg_ihm(k)
3195 ratio = min( abs(prr_rcg(k)), abs(prg_rcg(k)) )
3196 prr_rcg(k) = ratio * sign(1.0, sngl(prr_rcg(k)))
3197 prg_rcg(k) = -prr_rcg(k)
3198 if (temp(k).gt.t_0)
then
3199 ratio = min( abs(prr_rcs(k)), abs(prs_rcs(k)) )
3200 prr_rcs(k) = ratio * sign(1.0, sngl(prr_rcs(k)))
3201 prs_rcs(k) = -prr_rcs(k)
3212 lfus2 = lsub - lvap(k)
3215 if (is_aerosol_aware)
then
3216 nwfaten(k) = nwfaten(k) - (pna_rca(k) + pna_sca(k) &
3217 + pna_gca(k) + pni_iha(k)) * orho
3218 nifaten(k) = nifaten(k) - (pnd_rcd(k) + pnd_scd(k) &
3219 + pnd_gcd(k)) * orho
3221 nifaten(k) = nifaten(k) - pni_inu(k)*orho
3228 qvten(k) = qvten(k) + (-pri_inu(k) - pri_iha(k) - pri_ide(k) &
3229 - prs_ide(k) - prs_sde(k) - prg_gde(k)) &
3233 qcten(k) = qcten(k) + (-prr_wau(k) - pri_wfz(k) &
3234 - prr_rcw(k) - prs_scw(k) - prg_scw(k) &
3239 ncten(k) = ncten(k) + (-pnc_wau(k) - pnc_rcw(k) &
3240 - pni_wfz(k) - pnc_scw(k) - pnc_gcw(k)) &
3245 xrc=max(r1, (qc1d(k) + qcten(k)*dtsave)*rho(k))
3246 xnc=max(2., (nc1d(k) + ncten(k)*dtsave)*rho(k))
3247 if (xrc .gt. r1)
then
3248 if (xnc.gt.10000.e6)
then
3250 elseif (xnc.lt.100.)
then
3253 nu_c = nint(1000.e6/xnc) + 2
3254 nu_c = max(2, min(nu_c+nint(rand2), 15))
3256 lamc = (xnc*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
3257 xdc = (bm_r + nu_c + 1.) / lamc
3258 if (xdc.lt. d0c)
then
3259 lamc = cce(2,nu_c)/d0c
3260 xnc = ccg(1,nu_c)*ocg2(nu_c)*xrc/am_r*lamc**bm_r
3261 ncten(k) = (xnc-nc1d(k)*rho(k))*odts*orho
3262 elseif (xdc.gt. d0r*2.)
then
3263 lamc = cce(2,nu_c)/(d0r*2.)
3264 xnc = ccg(1,nu_c)*ocg2(nu_c)*xrc/am_r*lamc**bm_r
3265 ncten(k) = (xnc-nc1d(k)*rho(k))*odts*orho
3268 ncten(k) = -nc1d(k)*odts
3270 xnc=max(0.,(nc1d(k) + ncten(k)*dtsave)*rho(k))
3271 if (xnc.gt.nt_c_max) &
3272 ncten(k) = (nt_c_max-nc1d(k)*rho(k))*odts*orho
3275 qiten(k) = qiten(k) + (pri_inu(k) + pri_iha(k) + pri_ihm(k) &
3276 + pri_wfz(k) + pri_rfz(k) + pri_ide(k) &
3277 - prs_iau(k) - prs_sci(k) - pri_rci(k)) &
3281 niten(k) = niten(k) + (pni_inu(k) + pni_iha(k) + pni_ihm(k) &
3282 + pni_wfz(k) + pni_rfz(k) + pni_ide(k) &
3283 - pni_iau(k) - pni_sci(k) - pni_rci(k)) &
3288 xri=max(r1,(qi1d(k) + qiten(k)*dtsave)*rho(k))
3289 xni=max(r2,(ni1d(k) + niten(k)*dtsave)*rho(k))
3290 if (xri.gt. r1)
then
3291 lami = (am_i*cig(2)*oig1*xni/xri)**obmi
3293 xdi = (bm_i + mu_i + 1.) * ilami
3294 if (xdi.lt. 5.e-6)
then
3296 xni = min(4999.d3, cig(1)*oig2*xri/am_i*lami**bm_i)
3297 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
3298 elseif (xdi.gt. 300.e-6)
then
3299 lami = cie(2)/300.e-6
3300 xni = cig(1)*oig2*xri/am_i*lami**bm_i
3301 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
3304 niten(k) = -ni1d(k)*odts
3306 xni=max(0.,(ni1d(k) + niten(k)*dtsave)*rho(k))
3307 if (xni.gt.4999.e3) &
3308 niten(k) = (4999.e3-ni1d(k)*rho(k))*odts*orho
3311 qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) &
3312 + prr_sml(k) + prr_gml(k) + prr_rcs(k) &
3313 + prr_rcg(k) - prg_rfz(k) &
3314 - pri_rfz(k) - prr_rci(k)) &
3318 nrten(k) = nrten(k) + (pnr_wau(k) + pnr_sml(k) + pnr_gml(k) &
3319 - (pnr_rfz(k) + pnr_rcr(k) + pnr_rcg(k) &
3320 + pnr_rcs(k) + pnr_rci(k) + pni_rfz(k)) ) &
3325 xrr=max(r1,(qr1d(k) + qrten(k)*dtsave)*rho(k))
3326 xnr=max(r2,(nr1d(k) + nrten(k)*dtsave)*rho(k))
3327 if (xrr.gt. r1)
then
3328 lamr = (am_r*crg(3)*org2*xnr/xrr)**obmr
3329 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3330 if (mvd_r(k) .gt. 2.5e-3)
then
3332 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3333 xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
3334 nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
3335 elseif (mvd_r(k) .lt. d0r*0.75)
then
3337 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3338 xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
3339 nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
3342 qrten(k) = -qr1d(k)*odts
3343 nrten(k) = -nr1d(k)*odts
3347 qsten(k) = qsten(k) + (prs_iau(k) + prs_sde(k) &
3348 + prs_sci(k) + prs_scw(k) + prs_rcs(k) &
3349 + prs_ide(k) - prs_ihm(k) - prr_sml(k)) &
3353 qgten(k) = qgten(k) + (prg_scw(k) + prg_rfz(k) &
3354 + prg_gde(k) + prg_rcg(k) + prg_gcw(k) &
3355 + prg_rci(k) + prg_rcs(k) - prg_ihm(k) &
3360 if (temp(k).lt.t_0)
then
3362 + ( lsub*ocp(k)*(pri_inu(k) + pri_ide(k) &
3363 + prs_ide(k) + prs_sde(k) &
3364 + prg_gde(k) + pri_iha(k)) &
3365 + lfus2*ocp(k)*(pri_wfz(k) + pri_rfz(k) &
3366 + prg_rfz(k) + prs_scw(k) &
3367 + prg_scw(k) + prg_gcw(k) &
3368 + prg_rcs(k) + prs_rcs(k) &
3369 + prr_rci(k) + prg_rcg(k)) &
3373 + ( lfus*ocp(k)*(-prr_sml(k) - prr_gml(k) &
3374 - prr_rcg(k) - prr_rcs(k)) &
3375 + lsub*ocp(k)*(prs_sde(k) + prg_gde(k)) &
3385 temp(k) = t1d(k) + dt*tten(k)
3387 tempc = temp(k) - 273.15
3388 qv(k) = max(1.e-10, qv1d(k) + dt*qvten(k))
3389 rho(k) = 0.622*pres(k)/(r*temp(k)*(qv(k)+0.622))
3390 rhof(k) = sqrt(rho_not/rho(k))
3391 rhof2(k) = sqrt(rhof(k))
3392 qvs(k) =
rslf(pres(k), temp(k))
3393 ssatw(k) = qv(k)/qvs(k) - 1.
3394 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
3395 diffu(k) = 2.11e-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
3396 if (tempc .ge. 0.0)
then
3397 visco(k) = (1.718+0.0049*tempc)*1.0e-5
3399 visco(k) = (1.718+0.0049*tempc-1.2e-5*tempc*tempc)*1.0e-5
3401 vsc2(k) = sqrt(rho(k)/visco(k))
3402 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
3403 tcond(k) = (5.69 + 0.0168*tempc)*1.0e-5 * 418.936
3404 ocp(k) = 1./(cp*(1.+0.887*qv(k)))
3405 lvt2(k)=lvap(k)*lvap(k)*ocp(k)*orv*otemp*otemp
3406 if (is_aerosol_aware) &
3407 nwfa(k) = max(11.1e6*rho(k), (nwfa1d(k) + nwfaten(k)*dt)*rho(k))
3411 if ((qc1d(k) + qcten(k)*dt) .gt. r1)
then
3412 rc(k) = (qc1d(k) + qcten(k)*dt)*rho(k)
3413 nc(k) = max(2., min((nc1d(k)+ncten(k)*dt)*rho(k), nt_c_max))
3414 if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware))
then
3428 if ((qi1d(k) + qiten(k)*dt) .gt. r1)
then
3429 ri(k) = (qi1d(k) + qiten(k)*dt)*rho(k)
3430 ni(k) = max(r2, (ni1d(k) + niten(k)*dt)*rho(k))
3438 if ((qr1d(k) + qrten(k)*dt) .gt. r1)
then
3439 rr(k) = (qr1d(k) + qrten(k)*dt)*rho(k)
3440 nr(k) = max(r2, (nr1d(k) + nrten(k)*dt)*rho(k))
3442 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3443 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3444 if (mvd_r(k) .gt. 2.5e-3)
then
3446 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3447 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
3448 elseif (mvd_r(k) .lt. d0r*0.75)
then
3450 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3451 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
3459 if ((qs1d(k) + qsten(k)*dt) .gt. r1)
then
3460 rs(k) = (qs1d(k) + qsten(k)*dt)*rho(k)
3467 if ((qg1d(k) + qgten(k)*dt) .gt. r1)
then
3468 rg(k) = (qg1d(k) + qgten(k)*dt)*rho(k)
3480 if (.not. iiwarm)
then
3488 if (.not. l_qs(k)) cycle
3489 tc0 = min(-0.1, temp(k)-273.15)
3490 smob(k) = rs(k)*oams
3494 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3))
then
3497 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
3498 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
3499 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
3500 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
3501 + sa(10)*bm_s*bm_s*bm_s
3503 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
3504 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
3505 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
3506 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
3507 + sb(10)*bm_s*bm_s*bm_s
3508 smo2(k) = (smob(k)/a_)**(1./b_)
3512 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
3513 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
3514 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
3515 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
3516 + sa(10)*cse(1)*cse(1)*cse(1)
3518 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
3519 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
3520 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
3521 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
3522 smoc(k) = a_ * smo2(k)**b_
3525 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(14) &
3526 + sa(4)*tc0*cse(14) + sa(5)*tc0*tc0 &
3527 + sa(6)*cse(14)*cse(14) + sa(7)*tc0*tc0*cse(14) &
3528 + sa(8)*tc0*cse(14)*cse(14) + sa(9)*tc0*tc0*tc0 &
3529 + sa(10)*cse(14)*cse(14)*cse(14)
3531 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(14) + sb(4)*tc0*cse(14) &
3532 + sb(5)*tc0*tc0 + sb(6)*cse(14)*cse(14) &
3533 + sb(7)*tc0*tc0*cse(14) + sb(8)*tc0*cse(14)*cse(14) &
3534 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(14)*cse(14)*cse(14)
3535 smod(k) = a_ * smo2(k)**b_
3548 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3550 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3551 n0_r(k) = nr(k)*org2*lamr**cre(2)
3564 if ( (ssatw(k).gt. eps) .or. (ssatw(k).lt. -eps .and. &
3566 clap = (qv(k)-qvs(k))/(1. + lvt2(k)*qvs(k))
3568 fcd = qvs(k)* exp(lvt2(k)*clap) - qv(k) + clap
3569 dfcd = qvs(k)*lvt2(k)* exp(lvt2(k)*clap) + 1.
3570 clap = clap - fcd/dfcd
3572 xrc = rc(k) + clap*rho(k)
3574 if (xrc.gt. r1)
then
3575 prw_vcd(k) = clap*odt
3577 if (clap .gt. eps)
then
3578 if (is_aerosol_aware .or. merra2_aerosol_aware)
then
3579 xnc = max(2.,
activ_ncloud(temp(k), w1d(k)+rand3, nwfa(k), lsml))
3587 pnc_wcd(k) = 0.5*(xnc-nc(k) + abs(xnc-nc(k)))*odts*orho
3590 elseif (clap .lt. -eps .AND. ssatw(k).lt.-1.e-6 .AND. &
3591 is_aerosol_aware)
then
3592 tempc = temp(k) - 273.15
3595 rvs_p = rvs*otemp*(lvap(k)*otemp*orv - 1.)
3596 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*orv - 1.) &
3597 *otemp*(lvap(k)*otemp*orv - 1.) &
3598 + (-2.*lvap(k)*otemp*otemp*otemp*orv) &
3600 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
3601 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
3602 * rvs_pp/rvs_p * rvs/rvs_p
3603 alphsc = max(1.e-9, alphsc)
3605 if (abs(xsat).lt. 1.e-9) xsat=0.
3606 t1_evap = 2.*pi*( 1.0 - alphsc*xsat &
3607 + 2.*alphsc*alphsc*xsat*xsat &
3608 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
3611 dc_star = dsqrt(-2.d0*dt * t1_evap/(2.*pi) &
3612 * 4.*diffu(k)*ssatw(k)*rvs/rho_w)
3613 idx_d = max(1, min(int(1.e6*dc_star), nbc))
3615 idx_n = nint(1.0 + float(nbc) * dlog(nc(k)/t_nc(1)) / nic1)
3616 idx_n = max(1, min(idx_n, nbc))
3619 if (rc(k).gt. r_c(1))
then
3620 nic = nint(alog10(rc(k)))
3621 do nn = nic-1, nic+1
3623 if ( (rc(k)/10.**nn).ge.1.0 .and. &
3624 (rc(k)/10.**nn).lt.10.0)
goto 159
3627 idx_c = int(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
3628 idx_c = max(1, min(idx_c, ntb_c))
3635 prw_vcd(k) = max(dble(-rc(k)*0.99*orho*odt), prw_vcd(k))
3636 pnc_wcd(k) = max(dble(-nc(k)*0.99*orho*odt), &
3637 -tnc_wev(idx_d, idx_c, idx_n)*orho*odt)
3641 prw_vcd(k) = -rc(k)*orho*odt
3642 pnc_wcd(k) = -nc(k)*orho*odt
3647 qvten(k) = qvten(k) - prw_vcd(k)
3648 qcten(k) = qcten(k) + prw_vcd(k)
3649 ncten(k) = ncten(k) + pnc_wcd(k)
3650 if (is_aerosol_aware) &
3651 nwfaten(k) = nwfaten(k) - pnc_wcd(k)
3652 tten(k) = tten(k) + lvap(k)*ocp(k)*prw_vcd(k)*(1-ifdry)
3653 rc(k) = max(r1, (qc1d(k) + dt*qcten(k))*rho(k))
3654 if (rc(k).eq.r1) l_qc(k) = .false.
3655 nc(k) = max(2., min((nc1d(k)+ncten(k)*dt)*rho(k), nt_c_max))
3656 if (.NOT. (is_aerosol_aware .or. merra2_aerosol_aware))
then
3663 qv(k) = max(1.e-10, qv1d(k) + dt*qvten(k))
3664 temp(k) = t1d(k) + dt*tten(k)
3665 rho(k) = 0.622*pres(k)/(r*temp(k)*(qv(k)+0.622))
3666 qvs(k) =
rslf(pres(k), temp(k))
3667 ssatw(k) = qv(k)/qvs(k) - 1.
3676 if ( (ssatw(k).lt. -eps) .and. l_qr(k) &
3677 .and. (.not.(prw_vcd(k).gt. 0.)) )
then
3678 tempc = temp(k) - 273.15
3681 rhof(k) = sqrt(rho_not*orho)
3682 rhof2(k) = sqrt(rhof(k))
3683 diffu(k) = 2.11e-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
3684 if (tempc .ge. 0.0)
then
3685 visco(k) = (1.718+0.0049*tempc)*1.0e-5
3687 visco(k) = (1.718+0.0049*tempc-1.2e-5*tempc*tempc)*1.0e-5
3689 vsc2(k) = sqrt(rho(k)/visco(k))
3690 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
3691 tcond(k) = (5.69 + 0.0168*tempc)*1.0e-5 * 418.936
3692 ocp(k) = 1./(cp*(1.+0.887*qv(k)))
3695 rvs_p = rvs*otemp*(lvap(k)*otemp*orv - 1.)
3696 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*orv - 1.) &
3697 *otemp*(lvap(k)*otemp*orv - 1.) &
3698 + (-2.*lvap(k)*otemp*otemp*otemp*orv) &
3700 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
3701 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
3702 * rvs_pp/rvs_p * rvs/rvs_p
3703 alphsc = max(1.e-9, alphsc)
3704 xsat = min(-1.e-9, ssatw(k))
3705 t1_evap = 2.*pi*( 1.0 - alphsc*xsat &
3706 + 2.*alphsc*alphsc*xsat*xsat &
3707 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
3712 if (qv(k)/qvs(k) .lt. 0.95 .AND. rr(k)*orho.le.1.e-8)
then
3713 prv_rev(k) = rr(k)*orho*odts
3715 prv_rev(k) = t1_evap*diffu(k)*(-ssatw(k))*n0_r(k)*rvs &
3716 * (t1_qr_ev*ilamr(k)**cre(10) &
3717 + t2_qr_ev*vsc2(k)*rhof2(k)*((lamr+0.5*fv_r)**(-cre(11))))
3718 rate_max = min((rr(k)*orho*odts), (qvs(k)-qv(k))*odts)
3719 prv_rev(k) = min(dble(rate_max), prv_rev(k)*orho)
3728 IF (prr_gml(k).gt.0.0)
THEN
3729 eva_factor = min(1.0, 0.01+(0.99-0.01)*(tempc/20.0))
3730 prv_rev(k) = prv_rev(k)*eva_factor
3734 pnr_rev(k) = min(dble(nr(k)*0.99*orho*odts), &
3735 prv_rev(k) * nr(k)/rr(k))
3737 qrten(k) = qrten(k) - prv_rev(k)
3738 qvten(k) = qvten(k) + prv_rev(k)
3739 nrten(k) = nrten(k) - pnr_rev(k)
3740 if (is_aerosol_aware) &
3741 nwfaten(k) = nwfaten(k) + pnr_rev(k)
3742 tten(k) = tten(k) - lvap(k)*ocp(k)*prv_rev(k)*(1-ifdry)
3744 rr(k) = max(r1, (qr1d(k) + dt*qrten(k))*rho(k))
3745 qv(k) = max(1.e-10, qv1d(k) + dt*qvten(k))
3746 nr(k) = max(r2, (nr1d(k) + dt*nrten(k))*rho(k))
3747 temp(k) = t1d(k) + dt*tten(k)
3748 rho(k) = 0.622*pres(k)/(r*temp(k)*(qv(k)+0.622))
3751#if ( WRF_CHEM == 1 )
3753 evapprod(k) = prv_rev(k) - (min(zerod0,prs_sde(k)) + &
3754 min(zerod0,prg_gde(k)))
3755 rainprod(k) = prr_wau(k) + prr_rcw(k) + prs_scw(k) + &
3756 prg_scw(k) + prs_iau(k) + &
3757 prg_gcw(k) + prs_sci(k) + &
3773 do k = kte+1, kts, -1
3784 if (any(l_qr .eqv. .true.))
then
3787 rhof(k) = sqrt(rho_not/rho(k))
3789 if (rr(k).gt. r1)
then
3790 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3791 vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) &
3792 *((lamr+fv_r)**(-cre(6)))
3799 vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) &
3800 *((lamr+fv_r)**(-cre(7)))
3804 vtnrk(k) = vtnrk(k+1)
3807 if (max(vtrk(k),vtnrk(k)) .gt. 1.e-3)
then
3808 ksed1(1) = max(ksed1(1), k)
3809 delta_tp = dzq(k)/(max(vtrk(k),vtnrk(k)))
3810 nstep = max(nstep, int(dt/delta_tp + 1.))
3813 if (ksed1(1) .eq. kte) ksed1(1) = kte-1
3814 if (nstep .gt. 0) onstep(1) = 1./real(nstep)
3819 if (any(l_qc .eqv. .true.))
then
3822 if (rc(k) .gt. r2) ksed1(5) = k
3823 hgt_agl = hgt_agl + dzq(k)
3824 if (hgt_agl .gt. 500.0)
goto 151
3828 do k = ksed1(5), kts, -1
3830 if (rc(k) .gt. r1 .and. w1d(k) .lt. 1.e-1)
then
3831 if (nc(k).gt.10000.e6)
then
3833 elseif (nc(k).lt.100.)
then
3836 nu_c = nint(1000.e6/nc(k)) + 2
3837 nu_c = max(2, min(nu_c+nint(rand2), 15))
3839 lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
3841 vtc = rhof(k)*av_c*ccg(5,nu_c)*ocg2(nu_c) * ilamc**bv_c
3843 vtc = rhof(k)*av_c*ccg(4,nu_c)*ocg1(nu_c) * ilamc**bv_c
3851 if (.not. iiwarm)
then
3853 if (any(l_qi .eqv. .true.))
then
3858 if (ri(k).gt. r1)
then
3859 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
3861 vti = rhof(k)*av_i*cig(3)*oig2 * ilami**bv_i
3866 vti = rhof(k)*av_i*cig(6)/cig(7) * ilami**bv_i
3870 vtnik(k) = vtnik(k+1)
3873 if (vtik(k) .gt. 1.e-3)
then
3874 ksed1(2) = max(ksed1(2), k)
3875 delta_tp = dzq(k)/vtik(k)
3876 nstep = max(nstep, int(dt/delta_tp + 1.))
3879 if (ksed1(2) .eq. kte) ksed1(2) = kte-1
3880 if (nstep .gt. 0) onstep(2) = 1./real(nstep)
3885 if (any(l_qs .eqv. .true.))
then
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 if (prr_sml(k) .gt. 0.0)
then
3906 sr = rs(k)/(rs(k)+rr(k))
3907 vtsk(k) = vts*sr + (1.-sr)*vtrk(k)
3910 vtsk(k) = vts*vts_boost(k)
3918 if (vtsk(k) .gt. 1.e-3)
then
3919 ksed1(3) = max(ksed1(3), k)
3920 delta_tp = dzq(k)/vtsk(k)
3921 nstep = max(nstep, int(dt/delta_tp + 1.))
3924 if (ksed1(3) .eq. kte) ksed1(3) = kte-1
3925 if (nstep .gt. 0) onstep(3) = 1./real(nstep)
3930 if (any(l_qg .eqv. .true.))
then
3935 if (rg(k).gt. r1)
then
3936 vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
3937 if (temp(k).gt. t_0)
then
3938 vtgk(k) = max(vtg, vtrk(k))
3946 if (vtgk(k) .gt. 1.e-3)
then
3947 ksed1(4) = max(ksed1(4), k)
3948 delta_tp = dzq(k)/vtgk(k)
3949 nstep = max(nstep, int(dt/delta_tp + 1.))
3952 if (ksed1(4) .eq. kte) ksed1(4) = kte-1
3953 if (nstep .gt. 0) onstep(4) = 1./real(nstep)
3963 if (any(l_qr .eqv. .true.))
then
3964 nstep = nint(1./onstep(1))
3966 if(.not. sedi_semi)
then
3969 sed_r(k) = vtrk(k)*rr(k)
3970 sed_n(k) = vtnrk(k)*nr(k)
3975 qrten(k) = qrten(k) - sed_r(k)*odzq*onstep(1)*orho
3976 nrten(k) = nrten(k) - sed_n(k)*odzq*onstep(1)*orho
3977 rr(k) = max(r1, rr(k) - sed_r(k)*odzq*dt*onstep(1))
3978 nr(k) = max(r2, nr(k) - sed_n(k)*odzq*dt*onstep(1))
3979 pfll1(k) = pfll1(k) + sed_r(k)*dt*onstep(1)
3980 do k = ksed1(1), kts, -1
3983 qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) &
3984 *odzq*onstep(1)*orho
3985 nrten(k) = nrten(k) + (sed_n(k+1)-sed_n(k)) &
3986 *odzq*onstep(1)*orho
3987 rr(k) = max(r1, rr(k) + (sed_r(k+1)-sed_r(k)) &
3989 nr(k) = max(r2, nr(k) + (sed_n(k+1)-sed_n(k)) &
3991 pfll1(k) = pfll1(k) + sed_r(k)*dt*onstep(1)
3994 if (rr(kts).gt.r1*1000.) &
3995 pptrain = pptrain + sed_r(kts)*dt*onstep(1)
4000 niter = int(nstep/max(decfl,1)) + 1
4005 call semi_lagrange_sedim(kte,dzq,vtrk,rr,rainsfc,pfll,dtcfl,r1)
4006 call semi_lagrange_sedim(kte,dzq,vtnrk,nr,vtr,pdummy,dtcfl,r2)
4008 orhodt = 1./(rho(k)*dt)
4009 qrten(k) = qrten(k) + (rr(k) - rr_tmp(k)) * orhodt
4010 nrten(k) = nrten(k) + (nr(k) - nr_tmp(k)) * orhodt
4011 pfll1(k) = pfll1(k) + pfll(k)
4013 pptrain = pptrain + rainsfc
4015 do k = kte+1, kts, -1
4021 if (rr(k).gt. r1)
then
4022 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
4023 vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) &
4024 *((lamr+fv_r)**(-cre(6)))
4031 vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) &
4032 *((lamr+fv_r)**(-cre(7)))
4042 if (any(l_qc .eqv. .true.))
then
4044 sed_c(k) = vtck(k)*rc(k)
4045 sed_n(k) = vtnck(k)*nc(k)
4047 do k = ksed1(5), kts, -1
4050 qcten(k) = qcten(k) + (sed_c(k+1)-sed_c(k)) *odzq*orho
4051 ncten(k) = ncten(k) + (sed_n(k+1)-sed_n(k)) *odzq*orho
4052 rc(k) = max(r1, rc(k) + (sed_c(k+1)-sed_c(k)) *odzq*dt)
4053 nc(k) = max(10., nc(k) + (sed_n(k+1)-sed_n(k)) *odzq*dt)
4059 if (any(l_qi .eqv. .true.))
then
4060 nstep = nint(1./onstep(2))
4063 sed_i(k) = vtik(k)*ri(k)
4064 sed_n(k) = vtnik(k)*ni(k)
4069 qiten(k) = qiten(k) - sed_i(k)*odzq*onstep(2)*orho
4070 niten(k) = niten(k) - sed_n(k)*odzq*onstep(2)*orho
4071 ri(k) = max(r1, ri(k) - sed_i(k)*odzq*dt*onstep(2))
4072 ni(k) = max(r2, ni(k) - sed_n(k)*odzq*dt*onstep(2))
4073 pfil1(k) = pfil1(k) + sed_i(k)*dt*onstep(2)
4074 do k = ksed1(2), kts, -1
4077 qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) &
4078 *odzq*onstep(2)*orho
4079 niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) &
4080 *odzq*onstep(2)*orho
4081 ri(k) = max(r1, ri(k) + (sed_i(k+1)-sed_i(k)) &
4083 ni(k) = max(r2, ni(k) + (sed_n(k+1)-sed_n(k)) &
4085 pfil1(k) = pfil1(k) + sed_i(k)*dt*onstep(2)
4088 if (ri(kts).gt.r1*1000.) &
4089 pptice = pptice + sed_i(kts)*dt*onstep(2)
4095 if (any(l_qs .eqv. .true.))
then
4096 nstep = nint(1./onstep(3))
4099 sed_s(k) = vtsk(k)*rs(k)
4104 qsten(k) = qsten(k) - sed_s(k)*odzq*onstep(3)*orho
4105 rs(k) = max(r1, rs(k) - sed_s(k)*odzq*dt*onstep(3))
4106 pfil1(k) = pfil1(k) + sed_s(k)*dt*onstep(3)
4107 do k = ksed1(3), kts, -1
4110 qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) &
4111 *odzq*onstep(3)*orho
4112 rs(k) = max(r1, rs(k) + (sed_s(k+1)-sed_s(k)) &
4114 pfil1(k) = pfil1(k) + sed_s(k)*dt*onstep(3)
4117 if (rs(kts).gt.r1*1000.) &
4118 pptsnow = pptsnow + sed_s(kts)*dt*onstep(3)
4124 if (any(l_qg .eqv. .true.))
then
4125 nstep = nint(1./onstep(4))
4126 if(.not. sedi_semi)
then
4129 sed_g(k) = vtgk(k)*rg(k)
4134 qgten(k) = qgten(k) - sed_g(k)*odzq*onstep(4)*orho
4135 rg(k) = max(r1, rg(k) - sed_g(k)*odzq*dt*onstep(4))
4136 pfil1(k) = pfil1(k) + sed_g(k)*dt*onstep(4)
4137 do k = ksed1(4), kts, -1
4140 qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) &
4141 *odzq*onstep(4)*orho
4142 rg(k) = max(r1, rg(k) + (sed_g(k+1)-sed_g(k)) &
4144 pfil1(k) = pfil1(k) + sed_g(k)*dt*onstep(4)
4147 if (rg(kts).gt.r1*1000.) &
4148 pptgraul = pptgraul + sed_g(kts)*dt*onstep(4)
4153 niter = int(nstep/max(decfl,1)) + 1
4158 call semi_lagrange_sedim(kte,dzq,vtgk,rg,graulsfc,pfil,dtcfl,r1)
4160 orhodt = 1./(rho(k)*dt)
4161 qgten(k) = qgten(k) + (rg(k) - rg_tmp(k))*orhodt
4162 pfil1(k) = pfil1(k) + pfil(k)
4164 pptgraul = pptgraul + graulsfc
4165 do k = kte+1, kts, -1
4170 if (rg(k).gt. r1)
then
4171 ygra1 = alog10(max(1.e-9, rg(k)))
4172 zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1
4173 n0_exp = 10.**(zans1)
4174 n0_exp = max(dble(gonv_min), min(n0_exp, dble(gonv_max)))
4175 lam_exp = (n0_exp*am_g*cgg(1)/rg(k))**oge1
4176 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
4178 vtg = rhof(k)*av_g*cgg(6)*ogg3 * (1./lamg)**bv_g
4179 if (temp(k).gt. t_0)
then
4180 vtgk(k) = max(vtg, vtrk(k))
4194 if (.not. iiwarm)
then
4196 xri = max(0.0, qi1d(k) + qiten(k)*dt)
4197 if ( (temp(k).gt. t_0) .and. (xri.gt. 0.0) )
then
4198 qcten(k) = qcten(k) + xri*odt
4199 ncten(k) = ncten(k) + ni1d(k)*odt
4200 qiten(k) = qiten(k) - xri*odt
4201 niten(k) = -ni1d(k)*odt
4202 tten(k) = tten(k) - lfus*ocp(k)*xri*odt*(1-ifdry)
4207 xrc = max(0.0, qc1d(k) + qcten(k)*dt)
4208 if ( (temp(k).lt. hgfr) .and. (xrc.gt. 0.0) )
then
4209 lfus2 = lsub - lvap(k)
4210 xnc = nc1d(k) + ncten(k)*dt
4211 qiten(k) = qiten(k) + xrc*odt
4212 niten(k) = niten(k) + xnc*odt
4213 qcten(k) = qcten(k) - xrc*odt
4214 ncten(k) = ncten(k) - xnc*odt
4215 tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-ifdry)
4226 t1d(k) = t1d(k) + tten(k)*dt
4227 qv1d(k) = max(1.e-10, qv1d(k) + qvten(k)*dt)
4228 qc1d(k) = qc1d(k) + qcten(k)*dt
4229 nc1d(k) = max(2./rho(k), min(nc1d(k) + ncten(k)*dt, nt_c_max))
4230 if (is_aerosol_aware)
then
4231 nwfa1d(k) = max(11.1e6, min(9999.e6, &
4232 (nwfa1d(k)+nwfaten(k)*dt)))
4233 nifa1d(k) = max(nain1*0.01, min(9999.e6, &
4234 (nifa1d(k)+nifaten(k)*dt)))
4236 if (qc1d(k) .le. r1)
then
4240 if (nc1d(k)*rho(k).gt.10000.e6)
then
4242 elseif (nc1d(k)*rho(k).lt.100.)
then
4245 nu_c = nint(1000.e6/(nc1d(k)*rho(k))) + 2
4246 nu_c = max(2, min(nu_c+nint(rand2), 15))
4248 lamc = (am_r*ccg(2,nu_c)*ocg1(nu_c)*nc1d(k)/qc1d(k))**obmr
4249 xdc = (bm_r + nu_c + 1.) / lamc
4250 if (xdc.lt. d0c)
then
4251 lamc = cce(2,nu_c)/d0c
4252 elseif (xdc.gt. d0r*2.)
then
4253 lamc = cce(2,nu_c)/(d0r*2.)
4255 nc1d(k) = min(ccg(1,nu_c)*ocg2(nu_c)*qc1d(k)/am_r*lamc**bm_r,&
4256 dble(nt_c_max)/rho(k))
4259 qi1d(k) = qi1d(k) + qiten(k)*dt
4260 ni1d(k) = max(r2/rho(k), ni1d(k) + niten(k)*dt)
4261 if (qi1d(k) .le. r1)
then
4265 lami = (am_i*cig(2)*oig1*ni1d(k)/qi1d(k))**obmi
4267 xdi = (bm_i + mu_i + 1.) * ilami
4268 if (xdi.lt. 5.e-6)
then
4270 elseif (xdi.gt. 300.e-6)
then
4271 lami = cie(2)/300.e-6
4273 ni1d(k) = min(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, &
4276 qr1d(k) = qr1d(k) + qrten(k)*dt
4277 nr1d(k) = max(r2/rho(k), nr1d(k) + nrten(k)*dt)
4278 if (qr1d(k) .le. r1)
then
4282 lamr = (am_r*crg(3)*org2*nr1d(k)/qr1d(k))**obmr
4283 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
4284 if (mvd_r(k) .gt. 2.5e-3)
then
4286 elseif (mvd_r(k) .lt. d0r*0.75)
then
4289 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
4290 nr1d(k) = crg(2)*org3*qr1d(k)*lamr**bm_r / am_r
4292 qs1d(k) = qs1d(k) + qsten(k)*dt
4293 if (qs1d(k) .le. r1) qs1d(k) = 0.0
4294 qg1d(k) = qg1d(k) + qgten(k)*dt
4295 if (qg1d(k) .le. r1) qg1d(k) = 0.0
4299 calculate_extended_diagnostics:
if (ext_diag)
then
4301 if(prw_vcd(k).gt.0)
then
4302 prw_vcdc1(k) = prw_vcd(k)*dt
4303 elseif(prw_vcd(k).lt.0)
then
4304 prw_vcde1(k) = -1*prw_vcd(k)*dt
4307 tpri_inu1(k) = pri_inu(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4309 if(pri_ide(k).gt.0)
then
4310 tpri_ide1_d(k) = pri_ide(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4312 tpri_ide1_s(k) = -pri_ide(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4315 if(temp(k).lt.t_0)
then
4316 tprs_ide1(k) = prs_ide(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4319 if(prs_sde(k).gt.0)
then
4320 tprs_sde1_d(k) = prs_sde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4322 tprs_sde1_s(k) = -prs_sde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4325 if(prg_gde(k).gt.0)
then
4326 tprg_gde1_d(k) = prg_gde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4328 tprg_gde1_s(k) = -prg_gde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4331 tpri_iha1(k) = pri_iha(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
4332 tpri_wfz1(k) = pri_wfz(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4333 tpri_rfz1(k) = pri_rfz(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4334 tprg_rfz1(k) = prg_rfz(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4335 tprs_scw1(k) = prs_scw(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4336 tprg_scw1(k) = prg_scw(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4337 tprg_rcs1(k) = prg_rcs(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4339 if(temp(k).lt.t_0)
then
4340 tprs_rcs1(k) = prs_rcs(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4343 tprr_rci1(k) = prr_rci(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4345 if(temp(k).lt.t_0)
then
4346 tprg_rcg1(k) = prg_rcg(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
4349 if(prw_vcd(k).gt.0)
then
4350 tprw_vcd1_c(k) = lvap(k)*ocp(k)*prw_vcd(k)*(1-ifdry)*dt
4352 tprw_vcd1_e(k) = -lvap(k)*ocp(k)*prw_vcd(k)*(1-ifdry)*dt
4356 tprr_sml1(k) = prr_sml(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
4357 tprr_gml1(k) = prr_gml(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
4359 if(temp(k).ge.t_0)
then
4360 tprr_rcg1(k) = -prr_rcg(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
4363 if(temp(k).ge.t_0)
then
4364 tprr_rcs1(k) = -prr_rcs(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
4367 tprv_rev1(k) = lvap(k)*ocp(k)*prv_rev(k)*(1-ifdry)*dt
4368 tten1(k) = tten(k)*dt
4369 qvten1(k) = qvten(k)*dt
4370 qiten1(k) = qiten(k)*dt
4371 qrten1(k) = qrten(k)*dt
4372 qsten1(k) = qsten(k)*dt
4373 qgten1(k) = qgten(k)*dt
4374 niten1(k) = niten(k)*dt
4375 nrten1(k) = nrten(k)*dt
4376 ncten1(k) = ncten(k)*dt
4377 qcten1(k) = qcten(k)*dt
4379 endif calculate_extended_diagnostics