54 real :: missing_value = - 1.e10
56 logical :: module_is_initialized = .false.
57 logical :: qsmith_tables_initialized = .false.
59 character (len = 17) :: mod_name =
'gfdl_cloud_microphys'
61 real,
parameter :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6
62 real,
parameter :: rhos = 0.1e3, rhog = 0.4e3
63 real,
parameter :: grav = 9.80665
64 real,
parameter :: rdgas = 287.05
65 real,
parameter :: rvgas = 461.50
66 real,
parameter :: cp_air = 1004.6
67 real,
parameter :: hlv = 2.5e6
68 real,
parameter :: hlf = 3.3358e5
69 real,
parameter :: pi = 3.1415926535897931
74 real,
parameter :: cp_vap = 4.0 * rvgas
76 real,
parameter :: cv_air = cp_air - rdgas
78 real,
parameter :: cv_vap = 3.0 * rvgas
84 real,
parameter :: c_ice = 1972.0
85 real,
parameter :: c_liq = 4185.5
88 real,
parameter :: eps = rdgas / rvgas
89 real,
parameter :: zvir = rvgas / rdgas - 1.
91 real,
parameter :: t_ice = 273.16
92 real,
parameter :: table_ice = 273.16
95 real,
parameter :: e00 = 611.21
97 real,
parameter :: dc_vap = cp_vap - c_liq
98 real,
parameter :: dc_ice = c_liq - c_ice
100 real,
parameter :: hlv0 = hlv
102 real,
parameter :: hlf0 = hlf
105 real,
parameter :: lv0 = hlv0 - dc_vap * t_ice
106 real,
parameter :: li00 = hlf0 - dc_ice * t_ice
108 real,
parameter :: d2ice = dc_vap + dc_ice
109 real,
parameter :: li2 = lv0 + li00
111 real,
parameter :: qrmin = 1.e-8
112 real,
parameter :: qvmin = 1.e-20
113 real,
parameter :: qcmin = 1.e-12
115 real,
parameter :: vr_min = 1.e-3
116 real,
parameter :: vf_min = 1.e-5
118 real,
parameter :: dz_min = 1.e-2
120 real,
parameter :: sfcrho = 1.2
121 real,
parameter :: rhor = 1.e3
124 real,
parameter :: rnzr = 8.0e6
125 real,
parameter :: rnzs = 3.0e6
126 real,
parameter :: rnzg = 4.0e6
127 real,
parameter :: rnzh = 4.0e4
131 real,
parameter :: rhoh = 0.917e3
133 public rhor, rhos, rhog, rhoh, rnzr, rnzs, rnzg, rnzh
134 real :: cracs, csacr, cgacr, cgacs, csacw, craci, csaci, cgacw, cgaci, cracw
136 real :: cssub (5), cgsub (5), crevp (5), cgfr (2), csmlt (5), cgmlt (5)
139 real :: pie, rgrav, fac_rc
142 real :: lati, latv, lats, lat2, lcp, icp, tcp
149 integer :: icloud_f = 0
150 integer :: irain_f = 0
152 logical :: de_ice = .false.
153 logical :: sedi_transport = .true.
154 logical :: do_sedi_w = .false.
155 logical :: do_sedi_heat = .true.
156 logical :: prog_ccn = .false.
157 logical :: do_qa = .true.
158 logical :: rad_snow = .true.
159 logical :: rad_graupel = .true.
160 logical :: rad_rain = .true.
161 logical :: fix_negative = .false.
162 logical :: do_setup = .true.
163 logical :: p_nonhydro = .false.
165 real,
allocatable :: table (:), table2 (:), table3 (:), tablew (:)
166 real,
allocatable :: des (:), des2 (:), des3 (:), desw (:)
168 logical :: tables_are_initialized = .false.
173 real,
parameter :: dt_fr = 8.
190 real :: cld_min = 0.05
191 real :: tice = 273.16
195 real :: mp_time = 150.
199 real :: rh_inc = 0.25
200 real :: rh_inr = 0.25
201 real :: rh_ins = 0.25
205 real :: tau_r2g = 900.
206 real :: tau_smlt = 900.
207 real :: tau_g2r = 600.
208 real :: tau_imlt = 600.
209 real :: tau_i2s = 1000.
210 real :: tau_l2r = 900.
211 real :: tau_v2l = 150.
212 real :: tau_l2v = 300.
213 real :: tau_g2v = 900.
214 real :: tau_v2g = 21600.
218 real :: dw_land = 0.20
219 real :: dw_ocean = 0.10
226 real :: rthresh = 10.0e-6
236 real :: sat_adj0 = 0.90
238 real :: qc_crt = 5.0e-8
242 real :: ql_mlt = 2.0e-3
243 real :: qs_mlt = 1.0e-6
245 real :: ql_gen = 1.0e-3
246 real :: qi_gen = 1.82e-6
250 real :: ql0_max = 2.0e-3
251 real :: qi0_max = 1.0e-4
253 real :: qi0_crt = 1.0e-4
255 real :: qr0_crt = 1.0e-4
257 real :: qs0_crt = 1.0e-3
259 real :: c_paut = 0.55
260 real :: c_psaci = 0.02
261 real :: c_piacr = 5.0
262 real :: c_cracw = 0.9
263 real :: c_pgacs = 2.0e-3
272 logical :: const_vi = .false.
273 logical :: const_vs = .false.
274 logical :: const_vg = .false.
275 logical :: const_vr = .false.
293 logical :: fast_sat_adj = .false.
294 logical :: z_slope_liq = .true.
295 logical :: z_slope_ice = .false.
296 logical :: use_ccn = .false.
297 logical :: use_ppm = .false.
298 logical :: mono_prof = .true.
299 logical :: mp_print = .false.
300 logical :: do_hail = .false.
304 real :: log_10, tice0, t_wfr
306 integer :: reiflag = 1
310 logical :: tintqs = .false.
312 real :: rewmin = 5.0, rewmax = 10.0
313 real :: reimin = 10.0, reimax = 150.0
314 real :: rermin = 10.0, rermax = 10000.0
315 real :: resmin = 150.0, resmax = 10000.0
316 real :: regmin = 300.0, regmax = 10000.0
322 namelist / gfdl_cloud_microphysics_nml / &
323 mp_time, t_min, t_sub, tau_r2g, tau_smlt, tau_g2r, dw_land, dw_ocean, &
324 vi_fac, vr_fac, vs_fac, vg_fac, ql_mlt, do_qa, fix_negative, vi_max, &
325 vs_max, vg_max, vr_max, qs_mlt, qs0_crt, qi_gen, ql0_max, qi0_max, &
326 qi0_crt, qr0_crt, fast_sat_adj, rh_inc, rh_ins, rh_inr, const_vi, &
327 const_vs, const_vg, const_vr, use_ccn, rthresh, ccn_l, ccn_o, qc_crt, &
328 tau_g2v, tau_v2g, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v, &
329 tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs, &
330 z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice, &
331 rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof, &
332 do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, &
333 mp_print, reiflag, rewmin, rewmax, reimin, reimax, rermin, rermax, &
334 resmin, resmax, regmin, regmax, tintqs, do_hail
337 mp_time, t_min, t_sub, tau_r2g, tau_smlt, tau_g2r, dw_land, dw_ocean, &
338 vi_fac, vr_fac, vs_fac, vg_fac, ql_mlt, do_qa, fix_negative, vi_max, &
339 vs_max, vg_max, vr_max, qs_mlt, qs0_crt, qi_gen, ql0_max, qi0_max, &
340 qi0_crt, qr0_crt, fast_sat_adj, rh_inc, rh_ins, rh_inr, const_vi, &
341 const_vs, const_vg, const_vr, use_ccn, rthresh, ccn_l, ccn_o, qc_crt, &
342 tau_g2v, tau_v2g, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v, &
343 tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs, &
344 z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice, &
345 rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof, &
346 do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, &
347 mp_print, reiflag, rewmin, rewmax, reimin, reimax, rermin, rermax, &
348 resmin, resmax, regmin, regmax, tintqs, do_hail
359 iis, iie, jjs, jje, kks, kke, ktop, kbot, &
360 qv, ql, qr, qi, qs, qg, qa, qn, &
361 qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, pt_dt, pt, w, &
362 uin, vin, udt, vdt, dz, delp, area, dt_in, land, &
363 rain, snow, ice, graupel, hydrostatic, phys_hydrostatic, &
364 p, lradar, refl_10cm, reset, pfils, pflls)
369 integer,
intent (in) :: iis, iie, jjs, jje
370 integer,
intent (in) :: kks, kke
371 integer,
intent (in) :: ktop, kbot
373 real,
intent (in) :: dt_in
375 real,
intent (in),
dimension (iis:iie, jjs:jje) :: area
376 real,
intent (in),
dimension (iis:iie, jjs:jje) :: land
378 real,
intent (in),
dimension (iis:iie, jjs:jje, kks:kke) :: delp, dz, uin, vin
379 real,
intent (in),
dimension (iis:iie, jjs:jje, kks:kke) :: pt, qv, ql, qr, qg, qa, qn
381 real,
intent (inout),
dimension (iis:iie, jjs:jje, kks:kke) :: qi, qs
382 real,
intent (inout),
dimension (iis:iie, jjs:jje, kks:kke) :: pt_dt, qa_dt, udt, vdt, w
383 real,
intent (inout),
dimension (iis:iie, jjs:jje, kks:kke) :: qv_dt, ql_dt, qr_dt
384 real,
intent (inout),
dimension (iis:iie, jjs:jje, kks:kke) :: qi_dt, qs_dt, qg_dt
386 real,
intent (out),
dimension (iis:iie, jjs:jje) :: rain, snow, ice, graupel
388 logical,
intent (in) :: hydrostatic, phys_hydrostatic
391 real,
intent (in),
dimension (iis:iie, jjs:jje, kks:kke) :: p
392 logical,
intent (in) :: lradar
393 real,
intent (out),
dimension (iis:iie, jjs:jje, kks:kke) :: refl_10cm
394 logical,
intent (in) :: reset
395 real,
intent (out),
dimension (iis:iie, jjs:jje, kks:kke) :: pfils, pflls
398 logical :: melti = .false.
400 real :: mpdt, rdt, dts, convt, tot_prec
403 integer :: is, ie, js, je
405 integer :: days, ntimes, kflip
407 real,
dimension (iie-iis+1, jje-jjs+1) :: prec_mp, prec1, cond, w_var, rh0
409 real,
dimension (iie-iis+1, jje-jjs+1,kke-kks+1) :: vt_r, vt_s, vt_g, vt_i, qn2
411 real,
dimension (size(pt,1), size(pt,3)) :: m2_rain, m2_sol
416 real,
dimension(ktop:kbot):: qv1d, t1d, p1d, qr1d, qs1d, qg1d, dbz
431 if (phys_hydrostatic .or. hydrostatic)
then
440 d0_vap = c_vap - c_liq
441 lv00 = hlv0 - d0_vap * t_ice
443 if (hydrostatic) do_sedi_w = .false.
456 tcp = (latv + lati) / cp_air
464 mpdt = min(dt_in, mp_time)
466 ntimes = nint(dt_in / mpdt)
469 dts = dt_in / real(ntimes)
495 call mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, qg,&
496 qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes, &
497 rain(:, j), snow(:, j), graupel(:, j), ice(:, j), m2_rain, &
498 m2_sol, cond(:, j), area(:, j), land(:, j), udt, vdt, pt_dt, &
499 qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, w_var, vt_r, &
500 vt_s, vt_g, vt_i, qn2)
503 pfils(i, j, k) = m2_sol(i, k)
504 pflls(i, j, k) = m2_rain(i, k)
562 convt = 86400. * rdt * rgrav
565 rain(i, j) = rain(i, j) * convt
566 snow(i, j) = snow(i, j) * convt
567 ice(i, j) = ice(i, j) * convt
568 graupel(i, j) = graupel(i, j) * convt
569 prec_mp(i, j) = rain(i, j) + snow(i, j) + ice(i, j) + graupel(i, j)
650 kflip = kbot-ktop+1-k+1
651 t1d(k) = pt(i,j,kflip)
652 p1d(k) = p(i,j,kflip)
653 qv1d(k) = qv(i,j,kflip)/(1-qv(i,j,kflip))
654 qr1d(k) = qr(i,j,kflip)
655 qs1d(k) = qs(i,j,kflip)
656 qg1d(k) = qg(i,j,kflip)
659 t1d, p1d, dbz, ktop, kbot, i, j, melti)
661 kflip = kbot-ktop+1-k+1
662 refl_10cm(i,j,kflip) = max(-35., dbz(k))
677subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, &
678 qg, qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes, &
679 rain, snow, graupel, ice, m2_rain, m2_sol, cond, area1, land, &
680 u_dt, v_dt, pt_dt, qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, &
681 w_var, vt_r, vt_s, vt_g, vt_i, qn2)
685 logical,
intent (in) :: hydrostatic
687 integer,
intent (in) :: j, is, ie, js, je, ks, ke
688 integer,
intent (in) :: ntimes, ktop, kbot
690 real,
intent (in) :: dt_in
692 real,
intent (in),
dimension (is:) :: area1, land
694 real,
intent (in),
dimension (is:, js:, ks:) :: uin, vin, delp, pt, dz
695 real,
intent (in),
dimension (is:, js:, ks:) :: qv, ql, qr, qg, qa, qn
697 real,
intent (inout),
dimension (is:, js:, ks:) :: qi, qs
698 real,
intent (inout),
dimension (is:, js:, ks:) :: u_dt, v_dt, w, pt_dt, qa_dt
699 real,
intent (inout),
dimension (is:, js:, ks:) :: qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt
701 real,
intent (inout),
dimension (is:) :: rain, snow, ice, graupel, cond
703 real,
intent (out),
dimension (is:, js:) :: w_var
705 real,
intent (out),
dimension (is:, js:, ks:) :: vt_r, vt_s, vt_g, vt_i, qn2
707 real,
intent (out),
dimension (is:, ks:) :: m2_rain, m2_sol
709 real,
dimension (ktop:kbot) :: qvz, qlz, qrz, qiz, qsz, qgz, qaz
710 real,
dimension (ktop:kbot) :: vtiz, vtsz, vtgz, vtrz
711 real,
dimension (ktop:kbot) :: dp0, dp1, dz0, dz1
712 real,
dimension (ktop:kbot) :: qv0, ql0, qr0, qi0, qs0, qg0, qa0
713 real,
dimension (ktop:kbot) :: t0, den, den0, tz, p1, denfac
714 real,
dimension (ktop:kbot) :: ccn, c_praut, m1_rain, m1_sol, m1
715 real,
dimension (ktop:kbot) :: u0, v0, u1, v1, w1
717 real :: cpaut, rh_adj, rh_rain
718 real :: r1, s1, i1, g1, rdt, ccn0
720 real :: s_leng, t_land, t_ocean, h_var
721 real :: cvm, tmp, omq
722 real :: dqi, qio, qin
726 dts = dt_in / real(ntimes)
756 qio = qiz(k) - dt_in * qi_dt(i, j, k)
757 qin = max(qio, qi0_max)
758 if (qiz(k) > qin)
then
759 qsz(k) = qsz(k) + qiz(k) - qin
761 dqi = (qin - qio) * rdt
762 qs_dt(i, j, k) = qs_dt(i, j, k) + qi_dt(i, j, k) - dqi
774 dp1(k) = delp(i, j, k)
788 dp1(k) = dp1(k) * (1. - qvz(k))
789 omq = dp0(k) / dp1(k)
791 qvz(k) = qvz(k) * omq
792 qlz(k) = qlz(k) * omq
793 qrz(k) = qrz(k) * omq
794 qiz(k) = qiz(k) * omq
795 qsz(k) = qsz(k) * omq
796 qgz(k) = qgz(k) * omq
802 den0(k) = - dp1(k) / (grav * dz0(k))
803 p1(k) = den0(k) * rdgas * t0(k)
842 cpaut = c_paut * 0.104 * grav / 1.717e-5
847 ccn(k) = qn(i, j, k) * 1.e6
848 c_praut(k) = cpaut * (ccn(k) * rhor) ** (- 1. / 3.)
852 ccn0 = (ccn_l * land(i) + ccn_o * (1. - land(i))) * 1.e6
857 ccn0 = ccn0 * rdgas * tz(kbot) / p1(kbot)
859 tmp = cpaut * (ccn0 * rhor) ** (- 1. / 3.)
892 s_leng = sqrt(sqrt(area1(i) / 1.e10))
893 t_land = dw_land * s_leng
894 t_ocean = dw_ocean * s_leng
895 h_var = t_land * land(i) + t_ocean * (1. - land(i))
896 h_var = min(0.20, max(0.01, h_var))
903 rh_adj = 1. - h_var - rh_inc
904 rh_rain = max(0.35, rh_adj - rh_inr)
911 call neg_adj (ktop, kbot, tz, dp1, qvz, qlz, qrz, qiz, qsz, qgz)
927 denfac(k) = sqrt(sfcrho / den(k))
931 dz1(k) = dz0(k) * tz(k) / t0(k)
932 den(k) = den0(k) * dz0(k) / dz1(k)
933 denfac(k) = sqrt(sfcrho / den(k))
941 call warm_rain (dt_rain, ktop, kbot, dp1, dz1, tz, qvz, qlz, qrz, qiz, qsz, &
942 qgz, den, denfac, ccn, c_praut, rh_rain, vtrz, r1, m1_rain, w1, h_var)
944 rain(i) = rain(i) + r1
947 m2_rain(i, k) = m2_rain(i, k) + m1_rain(k)
948 m1(k) = m1(k) + m1_rain(k)
957 call fall_speed (ktop, kbot, den, qsz, qiz, qgz, qlz, tz, vtsz, vtiz, vtgz)
960 call terminal_fall (dts, ktop, kbot, tz, qvz, qlz, qrz, qgz, qsz, qiz, &
961 dz1, dp1, den, vtgz, vtsz, vtiz, r1, g1, s1, i1, m1_sol, w1)
963 rain(i) = rain(i) + r1
964 snow(i) = snow(i) + s1
965 graupel(i) = graupel(i) + g1
973 call sedi_heat (ktop, kbot, dp1, m1_sol, dz1, tz, qvz, qlz, qrz, qiz, &
980 call warm_rain (dt_rain, ktop, kbot, dp1, dz1, tz, qvz, qlz, qrz, qiz, qsz, &
981 qgz, den, denfac, ccn, c_praut, rh_rain, vtrz, r1, m1_rain, w1, h_var)
983 rain(i) = rain(i) + r1
986 m2_rain(i, k) = m2_rain(i, k) + m1_rain(k)
987 m2_sol(i, k) = m2_sol(i, k) + m1_sol(k)
988 m1(k) = m1(k) + m1_rain(k) + m1_sol(k)
995 call icloud (ktop, kbot, tz, p1, qvz, qlz, qrz, qiz, qsz, qgz, dp1, den, &
996 denfac, vtsz, vtgz, vtrz, qaz, rh_adj, rh_rain, dts, h_var)
1001 m2_rain(i, :) = m2_rain(i, :) * rdt * rgrav
1002 m2_sol(i, :) = m2_sol(i, :) * rdt * rgrav
1009 if (sedi_transport)
then
1010 do k = ktop + 1, kbot
1011 u1(k) = (dp0(k) * u1(k) + m1(k - 1) * u1(k - 1)) / (dp0(k) + m1(k - 1))
1012 v1(k) = (dp0(k) * v1(k) + m1(k - 1) * v1(k - 1)) / (dp0(k) + m1(k - 1))
1013 u_dt(i, j, k) = u_dt(i, j, k) + (u1(k) - u0(k)) * rdt
1014 v_dt(i, j, k) = v_dt(i, j, k) + (v1(k) - v0(k)) * rdt
1030 omq = dp1(k) / dp0(k)
1031 qv_dt(i, j, k) = qv_dt(i, j, k) + rdt * (qvz(k) - qv0(k)) * omq
1032 ql_dt(i, j, k) = ql_dt(i, j, k) + rdt * (qlz(k) - ql0(k)) * omq
1033 qr_dt(i, j, k) = qr_dt(i, j, k) + rdt * (qrz(k) - qr0(k)) * omq
1034 qi_dt(i, j, k) = qi_dt(i, j, k) + rdt * (qiz(k) - qi0(k)) * omq
1035 qs_dt(i, j, k) = qs_dt(i, j, k) + rdt * (qsz(k) - qs0(k)) * omq
1036 qg_dt(i, j, k) = qg_dt(i, j, k) + rdt * (qgz(k) - qg0(k)) * omq
1037 cvm = c_air + qvz(k) * c_vap + (qrz(k) + qlz(k)) * c_liq + (qiz(k) + qsz(k) + qgz(k)) * c_ice
1038 pt_dt(i, j, k) = pt_dt(i, j, k) + rdt * (tz(k) - t0(k)) * cvm / cp_air
1049 qa_dt(i, j, k) = qa_dt(i, j, k) + rdt * (qaz(k) / real(ntimes) - qa0(k))
1100subroutine sedi_heat (ktop, kbot, dm, m1, dz, tz, qv, ql, qr, qi, qs, qg, cw)
1106 integer,
intent (in) :: ktop, kbot
1108 real,
intent (in),
dimension (ktop:kbot) :: dm, m1, dz, qv, ql, qr, qi, qs, qg
1110 real,
intent (inout),
dimension (ktop:kbot) :: tz
1112 real,
intent (in) :: cw
1114 real,
dimension (ktop:kbot) :: dgz, cvn
1121 dgz(k) = - 0.5 * grav * dz(k)
1122 cvn(k) = dm(k) * (cv_air + qv(k) * cv_vap + (qr(k) + ql(k)) * &
1123 c_liq + (qi(k) + qs(k) + qg(k)) * c_ice)
1136 tmp = cvn(k) + m1(k) * cw
1137 tz(k) = (tmp * tz(k) + m1(k) * dgz(k)) / tmp
1144 do k = ktop + 1, kbot
1145 tz(k) = ((cvn(k) + cw * (m1(k) - m1(k - 1))) * tz(k) + m1(k - 1) * &
1146 cw * tz(k - 1) + dgz(k) * (m1(k - 1) + m1(k))) / (cvn(k) + cw * m1(k))
1155subroutine warm_rain (dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, &
1156 den, denfac, ccn, c_praut, rh_rain, vtr, r1, m1_rain, w1, h_var)
1160 integer,
intent (in) :: ktop, kbot
1162 real,
intent (in) :: dt
1163 real,
intent (in) :: rh_rain, h_var
1165 real,
intent (in),
dimension (ktop:kbot) :: dp, dz, den
1166 real,
intent (in),
dimension (ktop:kbot) :: denfac, ccn, c_praut
1168 real,
intent (inout),
dimension (ktop:kbot) :: tz, vtr
1169 real,
intent (inout),
dimension (ktop:kbot) :: qv, ql, qr, qi, qs, qg
1170 real,
intent (inout),
dimension (ktop:kbot) :: m1_rain, w1
1172 real,
intent (out) :: r1
1174 real,
parameter :: so3 = 7. / 3.
1176 real,
dimension (ktop:kbot) :: dl, dm
1177 real,
dimension (ktop:kbot + 1) :: ze, zt
1179 real :: sink, dq, qc0, qc
1188 real,
parameter :: vconr = 2503.23638966667
1189 real,
parameter :: normr = 25132741228.7183
1190 real,
parameter :: thr = 1.e-8
1217 qden = qr(k) * den(k)
1218 if (qr(k) < thr)
then
1221 vtr(k) = vr_fac * vconr * sqrt(min(10., sfcrho / den(k))) * &
1222 exp(0.2 * log(qden / normr))
1223 vtr(k) = min(vr_max, max(vr_min, vtr(k)))
1229 do k = kbot, ktop, - 1
1230 ze(k) = ze(k + 1) - dz(k)
1239 call revap_racc (ktop, kbot, dt5, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
1243 dm(k) = dp(k) * (1. + qv(k) + ql(k) + qr(k) + qi(k) + qs(k) + qg(k))
1254 do k = ktop + 1, kbot
1255 zt(k) = ze(k) - dt5 * (vtr(k - 1) + vtr(k))
1257 zt(kbot + 1) = zs - dt * vtr(kbot)
1260 if (zt(k + 1) >= zt(k)) zt(k + 1) = zt(k) - dz_min
1264 call implicit_fall (dt, ktop, kbot, ze, vtr, dp, qr, r1, m1_rain)
1274 w1(ktop) = (dm(ktop) * w1(ktop) + m1_rain(ktop) * vtr(ktop)) / (dm(ktop) - m1_rain(ktop))
1275 do k = ktop + 1, kbot
1276 w1(k) = (dm(k) * w1(k) - m1_rain(k - 1) * vtr(k - 1) + m1_rain(k) * vtr(k)) &
1277 / (dm(k) + m1_rain(k - 1) - m1_rain(k))
1286 call sedi_heat (ktop, kbot, dp, m1_rain, dz, tz, qv, ql, qr, qi, qs, qg, c_liq)
1293 call revap_racc (ktop, kbot, dt5, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
1303 if (irain_f /= 0)
then
1310 qc0 = fac_rc * ccn(k)
1311 if (tz(k) > t_wfr)
then
1322 sink = min(dq, dt * c_praut(k) * den(k) * exp(so3 * log(ql(k))))
1323 ql(k) = ql(k) - sink
1324 qr(k) = qr(k) + sink
1335 call linear_prof (kbot - ktop + 1, ql(ktop), dl(ktop), z_slope_liq, h_var)
1338 qc0 = fac_rc * ccn(k)
1339 if (tz(k) > t_wfr + dt_fr)
then
1340 dl(k) = min(max(1.e-6, dl(k)), 0.5 * ql(k))
1352 dq = 0.5 * (ql(k) + dl(k) - qc)
1361 sink = min(1., dq / dl(k)) * dt * c_praut(k) * den(k) * exp(so3 * log(ql(k)))
1362 ql(k) = ql(k) - sink
1363 qr(k) = qr(k) + sink
1375subroutine revap_racc (ktop, kbot, dt, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
1379 integer,
intent (in) :: ktop, kbot
1381 real,
intent (in) :: dt
1382 real,
intent (in) :: rh_rain, h_var
1384 real,
intent (in),
dimension (ktop:kbot) :: den, denfac
1386 real,
intent (inout),
dimension (ktop:kbot) :: tz, qv, qr, ql, qi, qs, qg
1388 real,
dimension (ktop:kbot) :: lhl, cvm, q_liq, q_sol, lcpk
1390 real :: dqv, qsat, dqsdt, evap, t2, qden, q_plus, q_minus, sink
1391 real :: qpz, dq, dqh, tin
1397 if (tz(k) > t_wfr .and. qr(k) > qrmin)
then
1403 lhl(k) = lv00 + d0_vap * tz(k)
1404 q_liq(k) = ql(k) + qr(k)
1405 q_sol(k) = qi(k) + qs(k) + qg(k)
1406 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1407 lcpk(k) = lhl(k) / cvm(k)
1409 tin = tz(k) - lcpk(k) * ql(k)
1411 qsat = wqs2(tin, den(k), dqsdt)
1412 dqh = max(ql(k), h_var * max(qpz, qcmin))
1413 dqh = min(dqh, 0.2 * qpz)
1427 if (dqv > qvmin .and. qsat > q_minus)
then
1428 if (qsat > q_plus)
then
1435 dq = 0.25 * (q_minus - qsat) ** 2 / dqh
1437 qden = qr(k) * den(k)
1439 evap = crevp(1) * t2 * dq * (crevp(2) * sqrt(qden) + crevp(3) * &
1440 exp(0.725 * log(qden))) / (crevp(4) * t2 + crevp(5) * qsat * den(k))
1441 evap = min(qr(k), dt * evap, dqv / (1. + lcpk(k) * dqsdt))
1447 qr(k) = qr(k) - evap
1448 qv(k) = qv(k) + evap
1449 q_liq(k) = q_liq(k) - evap
1450 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1451 tz(k) = tz(k) - evap * lhl(k) / cvm(k)
1459 if (qr(k) > qrmin .and. ql(k) > 1.e-6 .and. qsat < q_minus)
then
1460 sink = dt * denfac(k) * cracw * exp(0.95 * log(qr(k) * den(k)))
1461 sink = sink / (1. + sink) * ql(k)
1462 ql(k) = ql(k) - sink
1463 qr(k) = qr(k) + sink
1482 integer,
intent (in) :: km
1484 real,
intent (in) :: q (km), h_var
1486 real,
intent (out) :: dm (km)
1488 logical,
intent (in) :: z_var
1496 dq(k) = 0.5 * (q(k) - q(k - 1))
1506 dm(k) = 0.5 * min(abs(dq(k) + dq(k + 1)), 0.5 * q(k))
1507 if (dq(k) * dq(k + 1) <= 0.)
then
1508 if (dq(k) > 0.)
then
1509 dm(k) = min(dm(k), dq(k), - dq(k + 1))
1523 dm(k) = max(dm(k), qvmin, h_var * q(k))
1527 dm(k) = max(qvmin, h_var * q(k))
1543subroutine icloud (ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, &
1544 den, denfac, vts, vtg, vtr, qak, rh_adj, rh_rain, dts, h_var)
1548 integer,
intent (in) :: ktop, kbot
1550 real,
intent (in),
dimension (ktop:kbot) :: p1, dp1, den, denfac, vts, vtg, vtr
1552 real,
intent (inout),
dimension (ktop:kbot) :: tzk, qvk, qlk, qrk, qik, qsk, qgk, qak
1554 real,
intent (in) :: rh_adj, rh_rain, dts, h_var
1556 real,
dimension (ktop:kbot) :: lcpk, icpk, tcpk, di, lhl, lhi
1557 real,
dimension (ktop:kbot) :: cvm, q_liq, q_sol
1559 real :: rdts, fac_g2v, fac_v2g, fac_i2s, fac_imlt
1560 real :: tz, qv, ql, qr, qi, qs, qg, melt
1561 real :: pracs, psacw, pgacw, psacr, pgacr, pgaci, praci, psaci
1562 real :: pgmlt, psmlt, pgfr, pgaut, psaut, pgsub
1563 real :: tc, tsq, dqs0, qden, qim, qsm
1564 real :: dt5, factor, sink, qi_crt
1565 real :: tmp, qsw, qsi, dqsdt, dq
1566 real :: dtmp, qc, q_plus, q_minus
1578 fac_i2s = 1. - exp(- dts / tau_i2s)
1579 fac_g2v = 1. - exp(- dts / tau_g2v)
1580 fac_v2g = 1. - exp(- dts / tau_v2g)
1582 fac_imlt = 1. - exp(- dt5 / tau_imlt)
1589 lhi(k) = li00 + dc_ice * tzk(k)
1590 q_liq(k) = qlk(k) + qrk(k)
1591 q_sol(k) = qik(k) + qsk(k) + qgk(k)
1592 cvm(k) = c_air + qvk(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1593 icpk(k) = lhi(k) / cvm(k)
1604 if (tzk(k) > tice .and. qik(k) > qcmin)
then
1610 melt = min(qik(k), fac_imlt * (tzk(k) - tice) / icpk(k))
1611 tmp = min(melt, dim(ql_mlt, qlk(k)))
1612 qlk(k) = qlk(k) + tmp
1613 qrk(k) = qrk(k) + melt - tmp
1614 qik(k) = qik(k) - melt
1615 q_liq(k) = q_liq(k) + melt
1616 q_sol(k) = q_sol(k) - melt
1617 cvm(k) = c_air + qvk(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1618 tzk(k) = tzk(k) - melt * lhi(k) / cvm(k)
1620 elseif (tzk(k) < t_wfr .and. qlk(k) > qcmin)
then
1627 dtmp = t_wfr - tzk(k)
1628 factor = min(1., dtmp / dt_fr)
1629 sink = min(qlk(k) * factor, dtmp / icpk(k))
1630 qi_crt = qi_gen * min(qi_lim, 0.1 * (tice - tzk(k))) / den(k)
1631 tmp = min(sink, dim(qi_crt, qik(k)))
1632 qlk(k) = qlk(k) - sink
1633 qsk(k) = qsk(k) + sink - tmp
1634 qik(k) = qik(k) + tmp
1635 q_liq(k) = q_liq(k) - sink
1636 q_sol(k) = q_sol(k) + sink
1637 cvm(k) = c_air + qvk(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1638 tzk(k) = tzk(k) + sink * lhi(k) / cvm(k)
1647 call linear_prof (kbot - ktop + 1, qik(ktop), di(ktop), z_slope_ice, h_var)
1654 lhl(k) = lv00 + d0_vap * tzk(k)
1655 lhi(k) = li00 + dc_ice * tzk(k)
1656 lcpk(k) = lhl(k) / cvm(k)
1657 icpk(k) = lhi(k) / cvm(k)
1658 tcpk(k) = lcpk(k) + icpk(k)
1667 if (p1(k) < p_min) cycle
1681 if (tc .ge. 0.)
then
1687 dqs0 = ces0 / p1(k) - qv
1689 if (qs > qcmin)
then
1696 if (ql > qrmin)
then
1697 factor = denfac(k) * csacw * exp(0.8125 * log(qs * den(k)))
1698 psacw = factor / (1. + dts * factor) * ql
1708 if (qr > qrmin)
then
1709 psacr = min(
acr3d(vts(k), vtr(k), qr, qs, csacr, acco(1, 2), &
1711 pracs =
acr3d(vtr(k), vts(k), qs, qr, cracs, acco(1, 1), den(k))
1722 psmlt = max(0.,
smlt(tc, dqs0, qs * den(k), psacw, psacr, csmlt, &
1724 sink = min(qs, dts * (psmlt + pracs), tc / icpk(k))
1727 tmp = min(sink, dim(qs_mlt, ql))
1729 qr = qr + sink - tmp
1732 q_liq(k) = q_liq(k) + sink
1733 q_sol(k) = q_sol(k) - sink
1734 cvm(k) = c_air + qv * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1735 tz = tz - sink * lhi(k) / cvm(k)
1744 lhi(k) = li00 + dc_ice * tz
1745 icpk(k) = lhi(k) / cvm(k)
1751 if (qg > qcmin .and. tc > 0.)
then
1758 pgacr = min(
acr3d(vtg(k), vtr(k), qr, qg, cgacr, acco(1, 3), &
1766 if (ql > qrmin)
then
1767 factor = cgacw * qden / sqrt(den(k) * sqrt(sqrt(qden)))
1768 pgacw = factor / (1. + dts * factor) * ql
1775 pgmlt = dts * gmlt(tc, dqs0, qden, pgacw, pgacr, cgmlt, den(k))
1776 pgmlt = min(max(0., pgmlt), qg, tc / icpk(k))
1779 q_liq(k) = q_liq(k) + pgmlt
1780 q_sol(k) = q_sol(k) - pgmlt
1781 cvm(k) = c_air + qv * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1782 tz = tz - pgmlt * lhi(k) / cvm(k)
1796 if (qi > 3.e-7)
then
1798 if (qs > 1.e-7)
then
1803 factor = dts * denfac(k) * csaci * exp(0.05 * tc + 0.8125 * log(qs * den(k)))
1804 psaci = factor / (1. + factor) * qi
1815 if (qi0_crt < 0.)
then
1818 qim = qi0_crt / den(k)
1828 tmp = fac_i2s * exp(0.025 * tc)
1831 di(k) = max(di(k), qrmin)
1833 if (q_plus > (qim + qrmin))
then
1834 if (qim > (qi - di(k)))
then
1835 dq = (0.25 * (q_plus - qim) ** 2) / di(k)
1846 sink = min(0.75 * qi, psaci + psaut)
1854 if (qg > 1.e-6)
then
1859 factor = dts * cgaci * sqrt(den(k)) * qg
1860 pgaci = factor / (1. + factor) * qi
1877 if (qr > 1.e-7 .and. tc < 0.)
then
1889 if (qs > 1.e-7)
then
1890 psacr = dts *
acr3d(vts(k), vtr(k), qr, qs, csacr, acco(1, 2), den(k))
1899 pgfr = dts * cgfr(1) / den(k) * (exp(- cgfr(2) * tc) - 1.) * &
1900 exp(1.75 * log(qr * den(k)))
1910 factor = min(sink, qr, - tc / icpk(k)) / max(sink, qrmin)
1912 psacr = factor * psacr
1913 pgfr = factor * pgfr
1919 q_liq(k) = q_liq(k) - sink
1920 q_sol(k) = q_sol(k) + sink
1921 cvm(k) = c_air + qv * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1922 tz = tz + sink * lhi(k) / cvm(k)
1930 lhi(k) = li00 + dc_ice * tz
1931 icpk(k) = lhi(k) / cvm(k)
1937 if (qs > 1.e-7)
then
1943 if (qg > qrmin)
then
1944 sink = dts *
acr3d(vtg(k), vts(k), qs, qg, cgacs, acco(1, 4), den(k))
1953 qsm = qs0_crt / den(k)
1955 factor = dts * 1.e-3 * exp(0.09 * (tz - tice))
1956 sink = sink + factor / (1. + factor) * (qs - qsm)
1958 sink = min(qs, sink)
1964 if (qg > 1.e-7 .and. tz < tice0)
then
1970 if (ql > 1.e-6)
then
1972 factor = dts * cgacw * qden / sqrt(den(k) * sqrt(sqrt(qden)))
1973 pgacw = factor / (1. + factor) * ql
1982 if (qr > 1.e-6)
then
1983 pgacr = min(dts *
acr3d(vtg(k), vtr(k), qr, qg, cgacr, acco(1, 3), &
1989 sink = pgacr + pgacw
1990 factor = min(sink, dim(tice, tz) / icpk(k)) / max(sink, qrmin)
1991 pgacr = factor * pgacr
1992 pgacw = factor * pgacw
1994 sink = pgacr + pgacw
1998 q_liq(k) = q_liq(k) - sink
1999 q_sol(k) = q_sol(k) + sink
2000 cvm(k) = c_air + qv * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2001 tz = tz + sink * lhi(k) / cvm(k)
2021 call subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, rh_adj, tzk, qvk, &
2022 qlk, qrk, qik, qsk, qgk, qak, h_var, rh_rain)
2031 ql, qr, qi, qs, qg, qa, h_var, rh_rain)
2035 integer,
intent (in) :: ktop, kbot
2037 real,
intent (in),
dimension (ktop:kbot) :: p1, den, denfac
2039 real,
intent (in) :: dts, rh_adj, h_var, rh_rain
2041 real,
intent (inout),
dimension (ktop:kbot) :: tz, qv, ql, qr, qi, qs, qg, qa
2043 real,
dimension (ktop:kbot) :: lcpk, icpk, tcpk, tcp3, lhl, lhi
2044 real,
dimension (ktop:kbot) :: cvm, q_liq, q_sol, q_cond
2046 real :: fac_v2l, fac_l2v
2048 real :: pidep, qi_crt
2055 real :: rh, rqi, tin, qsw, qsi, qpz, qstar
2056 real :: dqsdt, dwsdt, dq, dq0, factor, tmp
2057 real :: q_plus, q_minus, dt_evap, dt_pisub
2058 real :: evap, sink, tc, pisub, q_adj, dtmp
2059 real :: pssub, pgsub, tsq, qden, fac_g2v, fac_v2g
2063 if (fast_sat_adj)
then
2073 fac_v2l = 1. - exp(- dt_evap / tau_v2l)
2074 fac_l2v = 1. - exp(- dt_evap / tau_l2v)
2076 fac_g2v = 1. - exp(- dts / tau_g2v)
2077 fac_v2g = 1. - exp(- dts / tau_v2g)
2084 lhl(k) = lv00 + d0_vap * tz(k)
2085 lhi(k) = li00 + dc_ice * tz(k)
2086 q_liq(k) = ql(k) + qr(k)
2087 q_sol(k) = qi(k) + qs(k) + qg(k)
2088 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2089 lcpk(k) = lhl(k) / cvm(k)
2090 icpk(k) = lhi(k) / cvm(k)
2091 tcpk(k) = lcpk(k) + icpk(k)
2092 tcp3(k) = lcpk(k) + icpk(k) * min(1., dim(tice, tz(k)) / (tice - t_wfr))
2097 if (p1(k) < p_min) cycle
2103 if (tz(k) < t_min)
then
2104 sink = dim(qv(k), 1.e-7)
2105 qv(k) = qv(k) - sink
2106 qi(k) = qi(k) + sink
2107 q_sol(k) = q_sol(k) + sink
2108 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2109 tz(k) = tz(k) + sink * (lhl(k) + lhi(k)) / cvm(k)
2110 if (.not. do_qa) qa(k) = qa(k) + 1.
2118 lhl(k) = lv00 + d0_vap * tz(k)
2119 lhi(k) = li00 + dc_ice * tz(k)
2120 lcpk(k) = lhl(k) / cvm(k)
2121 icpk(k) = lhi(k) / cvm(k)
2122 tcpk(k) = lcpk(k) + icpk(k)
2123 tcp3(k) = lcpk(k) + icpk(k) * min(1., dim(tice, tz(k)) / (tice - t_wfr))
2129 qpz = qv(k) + ql(k) + qi(k)
2130 tin = tz(k) - (lhl(k) * (ql(k) + qi(k)) + lhi(k) * qi(k)) / (c_air + &
2131 qpz * c_vap + qr(k) * c_liq + (qs(k) + qg(k)) * c_ice)
2132 if (tin > t_sub + 6.)
then
2133 rh = qpz / iqs1(tin, den(k))
2134 if (rh < rh_adj)
then
2147 qsw = wqs2(tz(k), den(k), dwsdt)
2154 factor = min(1., fac_l2v * (10. * dq0 / qsw))
2155 evap = min(ql(k), factor * dq0 / (1. + tcp3(k) * dwsdt))
2161 evap = dq0 / (1. + tcp3(k) * dwsdt)
2163 qv(k) = qv(k) + evap
2164 ql(k) = ql(k) - evap
2165 q_liq(k) = q_liq(k) - evap
2166 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2167 tz(k) = tz(k) - evap * lhl(k) / cvm(k)
2173 lhi(k) = li00 + dc_ice * tz(k)
2174 icpk(k) = lhi(k) / cvm(k)
2180 dtmp = t_wfr - tz(k)
2181 if (dtmp > 0. .and. ql(k) > qcmin)
then
2182 sink = min(ql(k), ql(k) * dtmp * 0.125, dtmp / icpk(k))
2183 ql(k) = ql(k) - sink
2184 qi(k) = qi(k) + sink
2185 q_liq(k) = q_liq(k) - sink
2186 q_sol(k) = q_sol(k) + sink
2187 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2188 tz(k) = tz(k) + sink * lhi(k) / cvm(k)
2195 lhi(k) = li00 + dc_ice * tz(k)
2196 icpk(k) = lhi(k) / cvm(k)
2202 if (fast_sat_adj)
then
2203 dt_pisub = 0.5 * dts
2207 if (ql(k) > qrmin .and. tc > 0.)
then
2208 sink = 3.3333e-10 * dts * (exp(0.66 * tc) - 1.) * den(k) * ql(k) * ql(k)
2209 sink = min(ql(k), tc / icpk(k), sink)
2210 ql(k) = ql(k) - sink
2211 qi(k) = qi(k) + sink
2212 q_liq(k) = q_liq(k) - sink
2213 q_sol(k) = q_sol(k) + sink
2214 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2215 tz(k) = tz(k) + sink * lhi(k) / cvm(k)
2223 lhl(k) = lv00 + d0_vap * tz(k)
2224 lhi(k) = li00 + dc_ice * tz(k)
2225 lcpk(k) = lhl(k) / cvm(k)
2226 icpk(k) = lhi(k) / cvm(k)
2227 tcpk(k) = lcpk(k) + icpk(k)
2233 if (tz(k) < tice)
then
2234 qsi = iqs2(tz(k), den(k), dqsdt)
2236 sink = dq / (1. + tcpk(k) * dqsdt)
2237 if (qi(k) > qrmin)
then
2240 pidep = dt_pisub * dq * 349138.78 * exp(0.875 * log(qi(k) * den(k))) &
2241 / (qsi * den(k) * lat2 / (0.0243 * rvgas * tz(k) ** 2) + 4.42478e4)
2249 qi_crt = qi_gen * min(qi_lim, 0.1 * tmp) / den(k)
2250 sink = min(sink, max(qi_crt - qi(k), pidep), tmp / tcpk(k))
2252 pidep = pidep * min(1., dim(tz(k), t_sub) * 0.2)
2253 sink = max(pidep, sink, - qi(k))
2255 qv(k) = qv(k) - sink
2256 qi(k) = qi(k) + sink
2257 q_sol(k) = q_sol(k) + sink
2258 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2259 tz(k) = tz(k) + sink * (lhl(k) + lhi(k)) / cvm(k)
2266 lhl(k) = lv00 + d0_vap * tz(k)
2267 lhi(k) = li00 + dc_ice * tz(k)
2268 lcpk(k) = lhl(k) / cvm(k)
2269 icpk(k) = lhi(k) / cvm(k)
2270 tcpk(k) = lcpk(k) + icpk(k)
2277 if (qs(k) > qrmin)
then
2278 qsi = iqs2(tz(k), den(k), dqsdt)
2279 qden = qs(k) * den(k)
2280 tmp = exp(0.65625 * log(qden))
2282 dq = (qsi - qv(k)) / (1. + tcpk(k) * dqsdt)
2283 pssub = cssub(1) * tsq * (cssub(2) * sqrt(qden) + cssub(3) * tmp * &
2284 sqrt(denfac(k))) / (cssub(4) * tsq + cssub(5) * qsi * den(k))
2285 pssub = (qsi - qv(k)) * dts * pssub
2286 if (pssub > 0.)
then
2287 pssub = min(pssub * min(1., dim(tz(k), t_sub) * 0.2), qs(k))
2289 if (tz(k) > tice)
then
2292 pssub = max(pssub, dq, (tz(k) - tice) / tcpk(k))
2295 qs(k) = qs(k) - pssub
2296 qv(k) = qv(k) + pssub
2297 q_sol(k) = q_sol(k) - pssub
2298 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2299 tz(k) = tz(k) - pssub * (lhl(k) + lhi(k)) / cvm(k)
2306 lhl(k) = lv00 + d0_vap * tz(k)
2307 lhi(k) = li00 + dc_ice * tz(k)
2308 lcpk(k) = lhl(k) / cvm(k)
2309 icpk(k) = lhi(k) / cvm(k)
2310 tcpk(k) = lcpk(k) + icpk(k)
2316 if (qg(k) > qrmin)
then
2317 qsi = iqs2(tz(k), den(k), dqsdt)
2318 dq = (qv(k) - qsi) / (1. + tcpk(k) * dqsdt)
2319 pgsub = (qv(k) / qsi - 1.) * qg(k)
2320 if (pgsub > 0.)
then
2321 if (tz(k) > tice)
then
2324 pgsub = min(fac_v2g * pgsub, 0.2 * dq, ql(k) + qr(k), &
2325 (tice - tz(k)) / tcpk(k))
2328 pgsub = max(fac_g2v * pgsub, dq) * min(1., dim(tz(k), t_sub) * 0.1)
2330 qg(k) = qg(k) + pgsub
2331 qv(k) = qv(k) - pgsub
2332 q_sol(k) = q_sol(k) + pgsub
2333 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2334 tz(k) = tz(k) + pgsub * (lhl(k) + lhi(k)) / cvm(k)
2342 lhl(k) = lv00 + d0_vap * tz(k)
2343 lcpk(k) = lhl(k) / cvm(k)
2349 if (qr(k) > qcmin)
then
2350 qsw = wqs2(tz(k), den(k), dqsdt)
2351 sink = min(qr(k), dim(rh_rain * qsw, qv(k)) / (1. + lcpk(k) * dqsdt))
2352 qv(k) = qv(k) + sink
2353 qr(k) = qr(k) - sink
2354 q_liq(k) = q_liq(k) - sink
2355 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2356 tz(k) = tz(k) - sink * lhl(k) / cvm(k)
2364 lhl(k) = lv00 + d0_vap * tz(k)
2365 cvm(k) = c_air + (qv(k) + q_liq(k) + q_sol(k)) * c_vap
2366 lcpk(k) = lhl(k) / cvm(k)
2379 q_sol(k) = qi(k) + qs(k)
2384 q_liq(k) = ql(k) + qr(k)
2388 q_cond(k) = q_liq(k) + q_sol(k)
2390 qpz = qv(k) + q_cond(k)
2396 tin = tz(k) - (lcpk(k) * q_cond(k) + icpk(k) * q_sol(k))
2404 if (tin <= t_wfr)
then
2406 qstar = iqs1(tin, den(k))
2407 elseif (tin >= tice)
then
2409 qstar =
wqs1(tin, den(k))
2412 qsi = iqs1(tin, den(k))
2413 qsw =
wqs1(tin, den(k))
2414 if (q_cond(k) > 3.e-6)
then
2415 rqi = q_sol(k) / q_cond(k)
2420 rqi = (tice - tin) / (tice - t_wfr)
2422 qstar = rqi * qsi + (1. - rqi) * qsw
2430 if (qpz > qrmin)
then
2432 dq = max(qcmin, h_var * qpz)
2435 if (qstar < q_minus)
then
2437 elseif (qstar < q_plus .and. q_cond(k) > qc_crt)
then
2438 qa(k) = qa(k) + (q_plus - qstar) / (dq + dq)
2450subroutine revap_rac1 (hydrostatic, is, ie, dt, tz, qv, ql, qr, qi, qs, qg, den, hvar)
2454 logical,
intent (in) :: hydrostatic
2456 integer,
intent (in) :: is, ie
2458 real,
intent (in) :: dt
2460 real,
intent (in),
dimension (is:ie) :: den, hvar, qi, qs, qg
2462 real,
intent (inout),
dimension (is:ie) :: tz, qv, qr, ql
2464 real,
dimension (is:ie) :: lcp2, denfac, q_liq, q_sol, cvm, lhl
2466 real :: dqv, qsat, dqsdt, evap, qden, q_plus, q_minus, sink
2467 real :: tin, t2, qpz, dq, dqh
2476 lhl(i) = lv00 + d0_vap * tz(i)
2477 q_liq(i) = ql(i) + qr(i)
2478 q_sol(i) = qi(i) + qs(i) + qg(i)
2479 cvm(i) = c_air + qv(i) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
2480 lcp2(i) = lhl(i) / cvm(i)
2485 if (qr(i) > qrmin .and. tz(i) > t_wfr)
then
2487 tin = tz(i) - lcp2(i) * ql(i)
2488 qsat = wqs2(tin, den(i), dqsdt)
2489 dqh = max(ql(i), hvar(i) * max(qpz, qcmin))
2503 if (dqv > qvmin .and. qsat > q_minus)
then
2504 if (qsat > q_plus)
then
2509 dq = 0.25 * (q_minus - qsat) ** 2 / dqh
2511 qden = qr(i) * den(i)
2513 evap = crevp(1) * t2 * dq * (crevp(2) * sqrt(qden) + crevp(3) * exp(0.725 * log(qden))) &
2514 / (crevp(4) * t2 + crevp(5) * qsat * den(i))
2515 evap = min(qr(i), dt * evap, dqv / (1. + lcp2(i) * dqsdt))
2516 qr(i) = qr(i) - evap
2517 qv(i) = qv(i) + evap
2518 q_liq(i) = q_liq(i) - evap
2519 cvm(i) = c_air + qv(i) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
2520 tz(i) = tz(i) - evap * lhl(i) / cvm(i)
2527 if (qr(i) > qrmin .and. ql(i) > 1.e-8 .and. qsat < q_plus)
then
2528 denfac(i) = sqrt(sfcrho / den(i))
2529 sink = dt * denfac(i) * cracw * exp(0.95 * log(qr(i) * den(i)))
2530 sink = sink / (1. + sink) * ql(i)
2531 ql(i) = ql(i) - sink
2532 qr(i) = qr(i) + sink
2543subroutine terminal_fall (dtm, ktop, kbot, tz, qv, ql, qr, qg, qs, qi, dz, dp, &
2544 den, vtg, vts, vti, r1, g1, s1, i1, m1_sol, w1)
2548 integer,
intent (in) :: ktop, kbot
2550 real,
intent (in) :: dtm
2552 real,
intent (in),
dimension (ktop:kbot) :: vtg, vts, vti, den, dp, dz
2554 real,
intent (inout),
dimension (ktop:kbot) :: qv, ql, qr, qg, qs, qi, tz, m1_sol, w1
2556 real,
intent (out) :: r1, g1, s1, i1
2558 real,
dimension (ktop:kbot + 1) :: ze, zt
2560 real :: qsat, dqsdt, dt5, evap, dtime
2561 real :: factor, frac
2562 real :: tmp, precip, tc, sink
2564 real,
dimension (ktop:kbot) :: lcpk, icpk, cvm, q_liq, q_sol, lhl, lhi
2565 real,
dimension (ktop:kbot) :: m1, dm
2575 fac_imlt = 1. - exp(- dt5 / tau_imlt)
2583 lhl(k) = lv00 + d0_vap * tz(k)
2584 lhi(k) = li00 + dc_ice * tz(k)
2585 q_liq(k) = ql(k) + qr(k)
2586 q_sol(k) = qi(k) + qs(k) + qg(k)
2587 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2588 lcpk(k) = lhl(k) / cvm(k)
2589 icpk(k) = lhi(k) / cvm(k)
2597 do k = ktop, kbot - 1
2598 if (tz(k) > tice)
then
2610 if (qi(k) > qcmin .and. tc > 0.)
then
2611 sink = min(qi(k), fac_imlt * tc / icpk(k))
2612 tmp = min(sink, dim(ql_mlt, ql(k)))
2614 qr(k) = qr(k) + sink - tmp
2615 qi(k) = qi(k) - sink
2616 q_liq(k) = q_liq(k) + sink
2617 q_sol(k) = q_sol(k) - sink
2618 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2619 tz(k) = tz(k) - sink * lhi(k) / cvm(k)
2628 if (dtm < 60.) k0 = kbot
2635 do k = kbot, ktop, - 1
2636 ze(k) = ze(k + 1) - dz(k)
2646 lhi(k) = li00 + dc_ice * tz(k)
2647 icpk(k) = lhi(k) / cvm(k)
2656 if (vi_fac < 1.e-5 .or. no_fall)
then
2660 do k = ktop + 1, kbot
2661 zt(k) = ze(k) - dt5 * (vti(k - 1) + vti(k))
2663 zt(kbot + 1) = zs - dtm * vti(kbot)
2666 if (zt(k + 1) >= zt(k)) zt(k + 1) = zt(k) - dz_min
2670 do k = kbot - 1, k0, - 1
2671 if (qi(k) > qrmin)
then
2673 if (zt(k + 1) >= ze(m))
exit
2674 if (zt(k) < ze(m + 1) .and. tz(m) > tice)
then
2675 dtime = min(1.0, (ze(m) - ze(m + 1)) / (max(vr_min, vti(k)) * tau_imlt))
2676 sink = min(qi(k) * dp(k) / dp(m), dtime * (tz(m) - tice) / icpk(m))
2677 tmp = min(sink, dim(ql_mlt, ql(m)))
2679 qr(m) = qr(m) - tmp + sink
2680 tz(m) = tz(m) - sink * icpk(m)
2681 qi(k) = qi(k) - sink * dp(m) / dp(k)
2690 dm(k) = dp(k) * (1. + qv(k) + ql(k) + qr(k) + qi(k) + qs(k) + qg(k))
2697 call implicit_fall (dtm, ktop, kbot, ze, vti, dp, qi, i1, m1_sol)
2701 w1(ktop) = (dm(ktop) * w1(ktop) + m1_sol(ktop) * vti(ktop)) / (dm(ktop) - m1_sol(ktop))
2702 do k = ktop + 1, kbot
2703 w1(k) = (dm(k) * w1(k) - m1_sol(k - 1) * vti(k - 1) + m1_sol(k) * vti(k)) &
2704 / (dm(k) + m1_sol(k - 1) - m1_sol(k))
2722 do k = ktop + 1, kbot
2723 zt(k) = ze(k) - dt5 * (vts(k - 1) + vts(k))
2725 zt(kbot + 1) = zs - dtm * vts(kbot)
2728 if (zt(k + 1) >= zt(k)) zt(k + 1) = zt(k) - dz_min
2732 do k = kbot - 1, k0, - 1
2733 if (qs(k) > qrmin)
then
2735 if (zt(k + 1) >= ze(m))
exit
2736 dtime = min(dtm, (ze(m) - ze(m + 1)) / (vr_min + vts(k)))
2737 if (zt(k) < ze(m + 1) .and. tz(m) > tice)
then
2738 dtime = min(1.0, dtime / tau_smlt)
2739 sink = min(qs(k) * dp(k) / dp(m), dtime * (tz(m) - tice) / icpk(m))
2740 tz(m) = tz(m) - sink * icpk(m)
2741 qs(k) = qs(k) - sink * dp(m) / dp(k)
2742 if (zt(k) < zs)
then
2743 r1 = r1 + sink * dp(m)
2746 qr(m) = qr(m) + sink
2749 if (qs(k) < qrmin)
exit
2757 dm(k) = dp(k) * (1. + qv(k) + ql(k) + qr(k) + qi(k) + qs(k) + qg(k))
2764 call implicit_fall (dtm, ktop, kbot, ze, vts, dp, qs, s1, m1)
2768 m1_sol(k) = m1_sol(k) + m1(k)
2772 w1(ktop) = (dm(ktop) * w1(ktop) + m1(ktop) * vts(ktop)) / (dm(ktop) - m1(ktop))
2773 do k = ktop + 1, kbot
2774 w1(k) = (dm(k) * w1(k) - m1(k - 1) * vts(k - 1) + m1(k) * vts(k)) &
2775 / (dm(k) + m1(k - 1) - m1(k))
2791 do k = ktop + 1, kbot
2792 zt(k) = ze(k) - dt5 * (vtg(k - 1) + vtg(k))
2794 zt(kbot + 1) = zs - dtm * vtg(kbot)
2797 if (zt(k + 1) >= zt(k)) zt(k + 1) = zt(k) - dz_min
2801 do k = kbot - 1, k0, - 1
2802 if (qg(k) > qrmin)
then
2804 if (zt(k + 1) >= ze(m))
exit
2805 dtime = min(dtm, (ze(m) - ze(m + 1)) / vtg(k))
2806 if (zt(k) < ze(m + 1) .and. tz(m) > tice)
then
2807 dtime = min(1., dtime / tau_g2r)
2808 sink = min(qg(k) * dp(k) / dp(m), dtime * (tz(m) - tice) / icpk(m))
2809 tz(m) = tz(m) - sink * icpk(m)
2810 qg(k) = qg(k) - sink * dp(m) / dp(k)
2811 if (zt(k) < zs)
then
2812 r1 = r1 + sink * dp(m)
2814 qr(m) = qr(m) + sink
2817 if (qg(k) < qrmin)
exit
2825 dm(k) = dp(k) * (1. + qv(k) + ql(k) + qr(k) + qi(k) + qs(k) + qg(k))
2832 call implicit_fall (dtm, ktop, kbot, ze, vtg, dp, qg, g1, m1)
2836 m1_sol(k) = m1_sol(k) + m1(k)
2840 w1(ktop) = (dm(ktop) * w1(ktop) + m1(ktop) * vtg(ktop)) / (dm(ktop) - m1(ktop))
2841 do k = ktop + 1, kbot
2842 w1(k) = (dm(k) * w1(k) - m1(k - 1) * vtg(k - 1) + m1(k) * vtg(k)) &
2843 / (dm(k) + m1(k - 1) - m1(k))
2859 integer,
intent (in) :: ktop, kbot
2861 real,
intent (in) :: q (ktop:kbot)
2863 logical,
intent (out) :: no_fall
2870 if (q(k) > qrmin)
then
2887 integer,
intent (in) :: ktop, kbot
2889 real,
intent (in) :: dt
2891 real,
intent (in),
dimension (ktop:kbot + 1) :: ze
2893 real,
intent (in),
dimension (ktop:kbot) :: vt, dp
2895 real,
intent (inout),
dimension (ktop:kbot) :: q
2897 real,
intent (out),
dimension (ktop:kbot) :: m1
2899 real,
intent (out) :: precip
2901 real,
dimension (ktop:kbot) :: dz, qm, dd
2906 dz(k) = ze(k) - ze(k + 1)
2915 qm(ktop) = q(ktop) / (dz(ktop) + dd(ktop))
2916 do k = ktop + 1, kbot
2917 qm(k) = (q(k) + dd(k - 1) * qm(k - 1)) / (dz(k) + dd(k))
2925 qm(k) = qm(k) * dz(k)
2932 m1(ktop) = q(ktop) - qm(ktop)
2933 do k = ktop + 1, kbot
2934 m1(k) = m1(k - 1) + q(k) - qm(k)
2943 q(k) = qm(k) / dp(k)
2956 integer,
intent (in) :: ktop, kbot
2958 real,
intent (in) :: zs
2960 logical,
intent (in) :: mono
2962 real,
intent (in),
dimension (ktop:kbot + 1) :: ze, zt
2964 real,
intent (in),
dimension (ktop:kbot) :: dp
2967 real,
intent (inout),
dimension (ktop:kbot) :: q, m1
2969 real,
intent (out) :: precip
2971 real,
dimension (ktop:kbot) :: qm, dz
2973 real :: a4 (4, ktop:kbot)
2975 real :: pl, pr, delz, esl
2977 integer :: k, k0, n, m
2979 real,
parameter :: r3 = 1. / 3., r23 = 2. / 3.
2986 dz(k) = zt(k) - zt(k + 1)
2988 a4(1, k) = q(k) / dz(k)
2996 call cs_profile (a4(1, ktop), dz(ktop), kbot - ktop + 1, mono)
3001 if (ze(k) <= zt(n) .and. ze(k) >= zt(n + 1))
then
3002 pl = (zt(n) - ze(k)) / dz(n)
3003 if (zt(n + 1) <= ze(k + 1))
then
3005 pr = (zt(n) - ze(k + 1)) / dz(n)
3006 qm(k) = a4(2, n) + 0.5 * (a4(4, n) + a4(3, n) - a4(2, n)) * (pr + pl) - &
3007 a4(4, n) * r3 * (pr * (pr + pl) + pl ** 2)
3008 qm(k) = qm(k) * (ze(k) - ze(k + 1))
3012 qm(k) = (ze(k) - zt(n + 1)) * (a4(2, n) + 0.5 * (a4(4, n) + &
3013 a4(3, n) - a4(2, n)) * (1. + pl) - a4(4, n) * (r3 * (1. + pl * (1. + pl))))
3017 if (ze(k + 1) < zt(m + 1))
then
3018 qm(k) = qm(k) + q(m)
3020 delz = zt(m) - ze(k + 1)
3022 qm(k) = qm(k) + delz * (a4(2, m) + 0.5 * esl * &
3023 (a4(3, m) - a4(2, m) + a4(4, m) * (1. - r23 * esl)))
3036 m1(ktop) = q(ktop) - qm(ktop)
3037 do k = ktop + 1, kbot
3038 m1(k) = m1(k - 1) + q(k) - qm(k)
3046 q(k) = qm(k) / dp(k)
3056 integer,
intent (in) :: km
3058 real,
intent (in) :: del (km)
3060 logical,
intent (in) :: do_mono
3062 real,
intent (inout) :: a4 (4, km)
3064 real,
parameter :: qp_min = 1.e-6
3068 real :: d4, bet, a_bot, grat, pmp, lac
3069 real :: pmp_1, lac_1, pmp_2, lac_2
3070 real :: da1, da2, a6da
3076 grat = del(2) / del(1)
3077 bet = grat * (grat + 0.5)
3078 q(1) = (2. * grat * (grat + 1.) * a4(1, 1) + a4(1, 2)) / bet
3079 gam(1) = (1. + grat * (grat + 1.5)) / bet
3082 d4 = del(k - 1) / del(k)
3083 bet = 2. + 2. * d4 - gam(k - 1)
3084 q(k) = (3. * (a4(1, k - 1) + d4 * a4(1, k)) - q(k - 1)) / bet
3088 a_bot = 1. + d4 * (d4 + 1.5)
3089 q(km + 1) = (2. * d4 * (d4 + 1.) * a4(1, km) + a4(1, km - 1) - a_bot * q(km)) &
3090 / (d4 * (d4 + 0.5) - a_bot * gam(km))
3093 q(k) = q(k) - gam(k) * q(k + 1)
3101 gam(k) = a4(1, k) - a4(1, k - 1)
3112 q(1) = max(q(1), 0.)
3113 q(2) = min(q(2), max(a4(1, 1), a4(1, 2)))
3114 q(2) = max(q(2), min(a4(1, 1), a4(1, 2)), 0.)
3121 if (gam(k - 1) * gam(k + 1) > 0.)
then
3122 q(k) = min(q(k), max(a4(1, k - 1), a4(1, k)))
3123 q(k) = max(q(k), min(a4(1, k - 1), a4(1, k)))
3125 if (gam(k - 1) > 0.)
then
3127 q(k) = max(q(k), min(a4(1, k - 1), a4(1, k)))
3130 q(k) = min(q(k), max(a4(1, k - 1), a4(1, k)))
3131 q(k) = max(q(k), 0.0)
3140 q(km) = min(q(km), max(a4(1, km - 1), a4(1, km)))
3141 q(km) = max(q(km), min(a4(1, km - 1), a4(1, km)), 0.)
3154 if (gam(k) * gam(k + 1) > 0.0)
then
3165 if (a4(1, k) < qp_min .or. extm(k - 1) .or. extm(k + 1))
then
3170 a4(4, k) = 6. * a4(1, k) - 3. * (a4(2, k) + a4(3, k))
3171 if (abs(a4(4, k)) > abs(a4(2, k) - a4(3, k)))
then
3173 pmp_1 = a4(1, k) - 2.0 * gam(k + 1)
3174 lac_1 = pmp_1 + 1.5 * gam(k + 2)
3175 a4(2, k) = min(max(a4(2, k), min(a4(1, k), pmp_1, lac_1)), &
3176 max(a4(1, k), pmp_1, lac_1))
3177 pmp_2 = a4(1, k) + 2.0 * gam(k)
3178 lac_2 = pmp_2 - 1.5 * gam(k - 1)
3179 a4(3, k) = min(max(a4(3, k), min(a4(1, k), pmp_2, lac_2)), &
3180 max(a4(1, k), pmp_2, lac_2))
3187 if (a4(1, k) < qp_min .or. extm(k - 1) .or. extm(k + 1))
then
3196 a4(4, k) = 6. * a4(1, k) - 3. * (a4(2, k) + a4(3, k))
3205 da1 = a4(3, k) - a4(2, k)
3207 a6da = a4(4, k) * da1
3208 if (a6da < - da2)
then
3209 a4(4, k) = 3. * (a4(2, k) - a4(1, k))
3210 a4(3, k) = a4(2, k) - a4(4, k)
3211 elseif (a6da > da2)
then
3212 a4(4, k) = 3. * (a4(3, k) - a4(1, k))
3213 a4(2, k) = a4(3, k) - a4(4, k)
3223 a4(2, km) = a4(1, km)
3224 a4(3, km) = a4(1, km)
3235 integer,
intent (in) :: km
3237 real,
intent (inout) :: a4 (4, km)
3239 real,
parameter :: r12 = 1. / 12.
3248 if (abs(a4(3, k) - a4(2, k)) < - a4(4, k))
then
3249 if ((a4(1, k) + 0.25 * (a4(3, k) - a4(2, k)) ** 2 / a4(4, k) + a4(4, k) * r12) < 0.)
then
3250 if (a4(1, k) < a4(3, k) .and. a4(1, k) < a4(2, k))
then
3254 elseif (a4(3, k) > a4(2, k))
then
3255 a4(4, k) = 3. * (a4(2, k) - a4(1, k))
3256 a4(3, k) = a4(2, k) - a4(4, k)
3258 a4(4, k) = 3. * (a4(3, k) - a4(1, k))
3259 a4(2, k) = a4(3, k) - a4(4, k)
3270subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg)
3274 integer,
intent (in) :: ktop, kbot
3276 real,
intent (in),
dimension (ktop:kbot) :: den, qs, qi, qg, ql, tk
3277 real,
intent (out),
dimension (ktop:kbot) :: vts, vti, vtg
3281 real,
parameter :: thi = 1.0e-8
3282 real,
parameter :: thg = 1.0e-8
3283 real,
parameter :: ths = 1.0e-8
3288 real,
parameter :: aa = - 4.14122e-5
3289 real,
parameter :: bb = - 0.00538922
3290 real,
parameter :: cc = - 0.0516344
3291 real,
parameter :: dd = 0.00216078
3292 real,
parameter :: ee = 1.9714
3296 real,
parameter :: vcons = 6.6280504
3297 real,
parameter :: vcong = 87.2382675
3298 real,
parameter :: vconh = vcong * sqrt(rhoh / rhog)
3299 real,
parameter :: norms = 942477796.076938
3300 real,
parameter :: normg = 5026548245.74367
3301 real,
parameter :: normh = pi * rhoh * rnzh
3303 real,
dimension (ktop:kbot) :: qden, tc, rhof
3319 rhof(k) = sqrt(min(10., sfcrho / den(k)))
3335 if (qi(k) < thi)
then
3338 tc(k) = tk(k) - tice
3339 vti(k) = (3. + log10(qi(k) * den(k))) * (tc(k) * (aa * tc(k) + bb) + cc) + dd * tc(k) + ee
3340 vti(k) = vi0 * exp(log_10 * vti(k)) * 0.9
3341 vti(k) = min(vi_max, max(vf_min, vti(k)))
3354 if (qs(k) < ths)
then
3357 vts(k) = vs_fac * vcons * rhof(k) * exp(0.0625 * log(qs(k) * den(k) / norms))
3358 vts(k) = min(vs_max, max(vf_min, vts(k)))
3371 if (qg(k) < thg)
then
3374 vtg(k) = vg_fac * vconh * rhof(k) * sqrt(sqrt(sqrt(qg(k) * den(k) / normh)))
3375 vtg(k) = min(vg_max, max(vf_min, vtg(k)))
3380 if (qg(k) < thg)
then
3383 vtg(k) = vg_fac * vcong * rhof(k) * sqrt(sqrt(sqrt(qg(k) * den(k) / normg)))
3384 vtg(k) = min(vg_max, max(vf_min, vtg(k)))
3399 real :: gcon, cd, scm3, pisq, act (8)
3400 real :: vdifu, tcond
3403 real :: hlts, hltc, ri50
3405 real,
parameter :: gam263 = 1.456943, gam275 = 1.608355, gam290 = 1.827363, &
3406 gam325 = 2.54925, gam350 = 3.323363, gam380 = 4.694155, &
3407 gam425 = 8.285063, gam450 = 11.631769, gam480 = 17.837789, &
3408 gam625 = 184.860962, gam680 = 496.604067
3420 real,
parameter :: acc (3) = (/ 5.0, 2.0, 0.5 /)
3426 pie = 4. * atan(1.0)
3430 fac_rc = (4. / 3.) * pie * rhor * rthresh ** 3
3435 den_rc = fac_rc * ccn_o * 1.e6
3437 den_rc = fac_rc * ccn_l * 1.e6
3453 scm3 = (visk / vdifu) ** (1. / 3.)
3455 cracs = pisq * rnzr * rnzs * rhos
3456 csacr = pisq * rnzr * rnzs * rhor
3458 cgacr = pisq * rnzr * rnzh * rhor
3459 cgacs = pisq * rnzh * rnzs * rhos
3461 cgacr = pisq * rnzr * rnzg * rhor
3462 cgacs = pisq * rnzg * rnzs * rhos
3464 cgacs = cgacs * c_pgacs
3469 act(1) = pie * rnzs * rhos
3470 act(2) = pie * rnzr * rhor
3472 act(6) = pie * rnzh * rhoh
3474 act(6) = pie * rnzg * rhog
3484 acco(i, k) = acc(i) / (act(2 * k - 1) ** ((7 - i) * 0.25) * act(2 * k) ** (i * 0.25))
3488 gcon = 40.74 * sqrt(sfcrho)
3490 csacw = pie * rnzs * clin * gam325 / (4. * act(1) ** 0.8125)
3493 craci = pie * rnzr * alin * gam380 / (4. * act(2) ** 0.95)
3494 csaci = csacw * c_psaci
3497 cgacw = pie * rnzh * gam350 * gcon / (4. * act(6) ** 0.875)
3499 cgacw = pie * rnzg * gam350 * gcon / (4. * act(6) ** 0.875)
3504 cgaci = cgacw * 0.05
3508 cracw = c_cracw * cracw
3512 cssub(1) = 2. * pie * vdifu * tcond * rvgas * rnzs
3514 cgsub(1) = 2. * pie * vdifu * tcond * rvgas * rnzh
3516 cgsub(1) = 2. * pie * vdifu * tcond * rvgas * rnzg
3518 crevp(1) = 2. * pie * vdifu * tcond * rvgas * rnzr
3519 cssub(2) = 0.78 / sqrt(act(1))
3520 cgsub(2) = 0.78 / sqrt(act(6))
3521 crevp(2) = 0.78 / sqrt(act(2))
3522 cssub(3) = 0.31 * scm3 * gam263 * sqrt(clin / visk) / act(1) ** 0.65625
3523 cgsub(3) = 0.31 * scm3 * gam275 * sqrt(gcon / visk) / act(6) ** 0.6875
3524 crevp(3) = 0.31 * scm3 * gam290 * sqrt(alin / visk) / act(2) ** 0.725
3525 cssub(4) = tcond * rvgas
3526 cssub(5) = hlts ** 2 * vdifu
3530 crevp(5) = hltc ** 2 * vdifu
3532 cgfr(1) = 20.e2 * pisq * rnzr * rhor / act(2) ** 1.75
3537 csmlt(1) = 2. * pie * tcond * rnzs / hltf
3538 csmlt(2) = 2. * pie * vdifu * rnzs * hltc / hltf
3541 csmlt(5) = ch2o / hltf
3546 cgmlt(1) = 2. * pie * tcond * rnzh / hltf
3547 cgmlt(2) = 2. * pie * vdifu * rnzh * hltc / hltf
3549 cgmlt(1) = 2. * pie * tcond * rnzg / hltf
3550 cgmlt(2) = 2. * pie * vdifu * rnzg * hltc / hltf
3554 cgmlt(5) = ch2o / hltf
3567 fn_nml, errmsg, errflg)
3571 integer,
intent (in) :: me
3572 integer,
intent (in) :: master
3573 integer,
intent (in) :: nlunit
3574 integer,
intent (in) :: logunit
3576 character (len = 64),
intent (in) :: fn_nml
3577 character (len = *),
intent (in) :: input_nml_file(:)
3578 character(len=*),
intent(out) :: errmsg
3579 integer,
intent(out) :: errflg
3598#ifdef INTERNAL_FILE_NML
3599 read (input_nml_file, nml = gfdl_cloud_microphysics_nml)
3601 inquire (file = trim(fn_nml), exist = exists)
3602 if (.not. exists)
then
3603 write (6, *)
'gfdl - mp :: namelist file: ', trim(fn_nml),
' does not exist'
3605 errmsg =
'ERROR(gfdl_cloud_microphys_mod_init): namelist file '//trim(fn_nml)//
' does not exist'
3608 open (unit = nlunit, file = fn_nml, action =
'read' , status =
'old', iostat = ios)
3611 read (nlunit, nml = gfdl_cloud_microphysics_nml)
3616 if (me == master)
then
3617 write (logunit, *)
" ================================================================== "
3618 write (logunit, *)
"gfdl_cloud_microphys_mod"
3619 write (logunit, nml = gfdl_cloud_microphysics_nml)
3684 module_is_initialized = .true.
3723 tables_are_initialized = .false.
3739 if (.not. qsmith_tables_initialized)
call qsmith_init
3741 qsmith_tables_initialized = .true.
3751real function acr3d (v1, v2, q1, q2, c, cac, rho)
3755 real,
intent (in) :: v1, v2, c, rho
3756 real,
intent (in) :: q1, q2
3757 real,
intent (in) :: cac (3)
3776 acr3d = c * abs(v1 - v2) * q1 * s2 * (cac(1) * t1 + cac(2) * sqrt(t1) * s2 + cac(3) * s1)
3786real function smlt (tc, dqs, qsrho, psacw, psacr, c, rho, rhofac)
3790 real,
intent (in) :: tc, dqs, qsrho, psacw, psacr, c (5), rho, rhofac
3792 smlt = (c(1) * tc / rho - c(2) * dqs) * (c(3) * sqrt(qsrho) + &
3793 c(4) * qsrho ** 0.65625 * sqrt(rhofac)) + c(5) * tc * (psacw + psacr)
3802real function gmlt (tc, dqs, qgrho, pgacw, pgacr, c, rho)
3806 real,
intent (in) :: tc, dqs, qgrho, pgacw, pgacr, c (5), rho
3808 gmlt = (c(1) * tc / rho - c(2) * dqs) * (c(3) * sqrt(qgrho) + &
3809 c(4) * qgrho ** 0.6875 / rho ** 0.25) + c(5) * tc * (pgacw + pgacr)
3828 integer,
parameter :: length = 2621
3832 if (.not. tables_are_initialized)
then
3845 allocate (table(length))
3846 allocate (table2(length))
3847 allocate (table3(length))
3848 allocate (tablew(length))
3849 allocate (des(length))
3850 allocate (des2(length))
3851 allocate (des3(length))
3852 allocate (desw(length))
3859 do i = 1, length - 1
3860 des(i) = max(0., table(i + 1) - table(i))
3861 des2(i) = max(0., table2(i + 1) - table2(i))
3862 des3(i) = max(0., table3(i + 1) - table3(i))
3863 desw(i) = max(0., tablew(i + 1) - tablew(i))
3865 des(length) = des(length - 1)
3866 des2(length) = des2(length - 1)
3867 des3(length) = des3(length - 1)
3868 desw(length) = desw(length - 1)
3870 tables_are_initialized = .true.
3888 real,
intent (in) :: ta, den
3890 real :: es, ap1, tmin
3894 tmin = table_ice - 160.
3895 ap1 = 10. * dim(ta, tmin) + 1.
3896 ap1 = min(2621., ap1)
3898 es = tablew(it) + (ap1 - it) * desw(it)
3899 wqs1 = es / (rvgas * ta * den)
3910real function wqs2 (ta, den, dqdt)
3917 real,
intent (in) :: ta, den
3919 real,
intent (out) :: dqdt
3921 real :: es, ap1, tmin
3925 tmin = table_ice - 160.
3927 if (.not. tables_are_initialized)
call qsmith_init
3929 ap1 = 10. * dim(ta, tmin) + 1.
3930 ap1 = min(2621., ap1)
3932 es = tablew(it) + (ap1 - it) * desw(it)
3933 wqs2 = es / (rvgas * ta * den)
3936 dqdt = 10. * (desw(it) + (ap1 - it) * (desw(it + 1) - desw(it))) / (rvgas * ta * den)
3946real function wet_bulb (q, t, den)
3950 real,
intent (in) :: t, q, den
3952 real :: qs, tp, dqdt
3955 qs = wqs2(wet_bulb, den, dqdt)
3956 tp = 0.5 * (qs - q) / (1. + lcp * dqdt) * lcp
3957 wet_bulb = wet_bulb - tp
3961 qs = wqs2(wet_bulb, den, dqdt)
3962 tp = (qs - q) / (1. + lcp * dqdt) * lcp
3963 wet_bulb = wet_bulb - tp
3966end function wet_bulb
3973real function iqs1 (ta, den)
3980 real,
intent (in) :: ta, den
3982 real :: es, ap1, tmin
3986 tmin = table_ice - 160.
3987 ap1 = 10. * dim(ta, tmin) + 1.
3988 ap1 = min(2621., ap1)
3990 es = table2(it) + (ap1 - it) * des2(it)
3991 iqs1 = es / (rvgas * ta * den)
4000real function iqs2 (ta, den, dqdt)
4007 real,
intent (in) :: ta, den
4009 real,
intent (out) :: dqdt
4011 real :: es, ap1, tmin
4015 tmin = table_ice - 160.
4016 ap1 = 10. * dim(ta, tmin) + 1.
4017 ap1 = min(2621., ap1)
4019 es = table2(it) + (ap1 - it) * des2(it)
4020 iqs2 = es / (rvgas * ta * den)
4022 dqdt = 10. * (des2(it) + (ap1 - it) * (des2(it + 1) - des2(it))) / (rvgas * ta * den)
4031real function qs1d_moist (ta, qv, pa, dqdt)
4035 real,
intent (in) :: ta, pa, qv
4037 real,
intent (out) :: dqdt
4039 real :: es, ap1, tmin, eps10
4043 tmin = table_ice - 160.
4045 ap1 = 10. * dim(ta, tmin) + 1.
4046 ap1 = min(2621., ap1)
4048 es = table2(it) + (ap1 - it) * des2(it)
4049 qs1d_moist = eps * es * (1. + zvir * qv) / pa
4051 dqdt = eps10 * (des2(it) + (ap1 - it) * (des2(it + 1) - des2(it))) * (1. + zvir * qv) / pa
4053end function qs1d_moist
4061real function wqsat2_moist (ta, qv, pa, dqdt)
4065 real,
intent (in) :: ta, pa, qv
4067 real,
intent (out) :: dqdt
4069 real :: es, ap1, tmin, eps10
4073 tmin = table_ice - 160.
4075 ap1 = 10. * dim(ta, tmin) + 1.
4076 ap1 = min(2621., ap1)
4078 es = tablew(it) + (ap1 - it) * desw(it)
4079 wqsat2_moist = eps * es * (1. + zvir * qv) / pa
4081 dqdt = eps10 * (desw(it) + (ap1 - it) * (desw(it + 1) - desw(it))) * (1. + zvir * qv) / pa
4083end function wqsat2_moist
4091real function wqsat_moist (ta, qv, pa)
4095 real,
intent (in) :: ta, pa, qv
4097 real :: es, ap1, tmin
4101 tmin = table_ice - 160.
4102 ap1 = 10. * dim(ta, tmin) + 1.
4103 ap1 = min(2621., ap1)
4105 es = tablew(it) + (ap1 - it) * desw(it)
4106 wqsat_moist = eps * es * (1. + zvir * qv) / pa
4108end function wqsat_moist
4115real function qs1d_m (ta, qv, pa)
4119 real,
intent (in) :: ta, pa, qv
4121 real :: es, ap1, tmin
4125 tmin = table_ice - 160.
4126 ap1 = 10. * dim(ta, tmin) + 1.
4127 ap1 = min(2621., ap1)
4129 es = table2(it) + (ap1 - it) * des2(it)
4130 qs1d_m = eps * es * (1. + zvir * qv) / pa
4139real function d_sat (ta, den)
4143 real,
intent (in) :: ta, den
4145 real :: es_w, es_i, ap1, tmin
4149 tmin = table_ice - 160.
4150 ap1 = 10. * dim(ta, tmin) + 1.
4151 ap1 = min(2621., ap1)
4153 es_w = tablew(it) + (ap1 - it) * desw(it)
4154 es_i = table2(it) + (ap1 - it) * des2(it)
4155 d_sat = dim(es_w, es_i) / (rvgas * ta * den)
4164real function esw_table (ta)
4168 real,
intent (in) :: ta
4174 tmin = table_ice - 160.
4175 ap1 = 10. * dim(ta, tmin) + 1.
4176 ap1 = min(2621., ap1)
4178 esw_table = tablew(it) + (ap1 - it) * desw(it)
4180end function esw_table
4187real function es2_table (ta)
4191 real,
intent (in) :: ta
4197 tmin = table_ice - 160.
4198 ap1 = 10. * dim(ta, tmin) + 1.
4199 ap1 = min(2621., ap1)
4201 es2_table = table2(it) + (ap1 - it) * des2(it)
4203end function es2_table
4213 integer,
intent (in) :: n
4215 real,
intent (in) :: ta (n)
4217 real,
intent (out) :: es (n)
4223 tmin = table_ice - 160.
4226 ap1 = 10. * dim(ta(i), tmin) + 1.
4227 ap1 = min(2621., ap1)
4229 es(i) = tablew(it) + (ap1 - it) * desw(it)
4242 integer,
intent (in) :: n
4244 real,
intent (in) :: ta (n)
4246 real,
intent (out) :: es (n)
4252 tmin = table_ice - 160.
4255 ap1 = 10. * dim(ta(i), tmin) + 1.
4256 ap1 = min(2621., ap1)
4258 es(i) = table2(it) + (ap1 - it) * des2(it)
4271 integer,
intent (in) :: n
4273 real,
intent (in) :: ta (n)
4275 real,
intent (out) :: es (n)
4281 tmin = table_ice - 160.
4284 ap1 = 10. * dim(ta(i), tmin) + 1.
4285 ap1 = min(2621., ap1)
4287 es(i) = table3(it) + (ap1 - it) * des3(it)
4300 integer,
intent (in) :: n
4303 real :: tmin, tem, fac0, fac1, fac2
4307 tmin = table_ice - 160.
4314 tem = tmin + delt * real(i - 1)
4315 fac0 = (tem - t_ice) / (tem * t_ice)
4317 fac2 = (dc_vap * log(tem / t_ice) + fac1) / rvgas
4318 tablew(i) = e00 * exp(fac2)
4331 integer,
intent (in) :: n
4334 real :: tmin, tem0, tem1, fac0, fac1, fac2
4336 integer :: i, i0, i1
4338 tmin = table_ice - 160.
4341 tem0 = tmin + delt * real(i - 1)
4342 fac0 = (tem0 - t_ice) / (tem0 * t_ice)
4348 fac2 = (d2ice * log(tem0 / t_ice) + fac1) / rvgas
4354 fac2 = (dc_vap * log(tem0 / t_ice) + fac1) / rvgas
4356 table2(i) = e00 * exp(fac2)
4365 tem0 = 0.25 * (table2(i0 - 1) + 2. * table(i0) + table2(i0 + 1))
4366 tem1 = 0.25 * (table2(i1 - 1) + 2. * table(i1) + table2(i1 + 1))
4380 integer,
intent (in) :: n
4383 real :: esbasw, tbasw, esbasi, tmin, tem, aa, b, c, d, e
4386 integer :: i, i0, i1
4389 tbasw = table_ice + 100.
4391 tmin = table_ice - 160.
4394 tem = tmin + delt * real(i - 1)
4401 aa = - 9.09718 * (table_ice / tem - 1.)
4402 b = - 3.56654 * alog10(table_ice / tem)
4403 c = 0.876793 * (1. - tem / table_ice)
4405 table3(i) = 0.1 * 10 ** (aa + b + c + e)
4411 aa = - 7.90298 * (tbasw / tem - 1.)
4412 b = 5.02808 * alog10(tbasw / tem)
4413 c = - 1.3816e-7 * (10 ** ((1. - tem / tbasw) * 11.344) - 1.)
4414 d = 8.1328e-3 * (10 ** ((tbasw / tem - 1.) * (- 3.49149)) - 1.)
4416 table3(i) = 0.1 * 10 ** (aa + b + c + d + e)
4426 tem0 = 0.25 * (table3(i0 - 1) + 2. * table(i0) + table3(i0 + 1))
4427 tem1 = 0.25 * (table3(i1 - 1) + 2. * table(i1) + table3(i1 + 1))
4443 real,
intent (in) :: t, p, q
4445 real :: es, ap1, tmin
4449 tmin = table_ice - 160.
4450 ap1 = 10. * dim(t, tmin) + 1.
4451 ap1 = min(2621., ap1)
4453 es = table(it) + (ap1 - it) * des(it)
4454 qs_blend = eps * es * (1. + zvir * q) / p
4466 integer,
intent (in) :: n
4469 real :: tmin, tem, esh20
4470 real :: wice, wh2o, fac0, fac1, fac2
4475 tmin = table_ice - 160.
4482 tem = tmin + delt * real(i - 1)
4483 fac0 = (tem - t_ice) / (tem * t_ice)
4485 fac2 = (d2ice * log(tem / t_ice) + fac1) / rvgas
4486 table(i) = e00 * exp(fac2)
4494 tem = 253.16 + delt * real(i - 1)
4495 fac0 = (tem - t_ice) / (tem * t_ice)
4497 fac2 = (dc_vap * log(tem / t_ice) + fac1) / rvgas
4498 esh20 = e00 * exp(fac2)
4502 table(i + 1400) = esh20
4511 tem = 253.16 + delt * real(i - 1)
4512 wice = 0.05 * (table_ice - tem)
4513 wh2o = 0.05 * (tem - 253.16)
4514 table(i + 1400) = wice * table(i + 1400) + wh2o * esupc(i)
4526subroutine qsmith (im, km, ks, t, p, q, qs, dqdt)
4530 integer,
intent (in) :: im, km, ks
4532 real,
intent (in),
dimension (im, km) :: t, p, q
4534 real,
intent (out),
dimension (im, km) :: qs
4536 real,
intent (out),
dimension (im, km),
optional :: dqdt
4538 real :: eps10, ap1, tmin
4540 real,
dimension (im, km) :: es
4544 tmin = table_ice - 160.
4547 if (.not. tables_are_initialized)
then
4553 ap1 = 10. * dim(t(i, k), tmin) + 1.
4554 ap1 = min(2621., ap1)
4556 es(i, k) = table(it) + (ap1 - it) * des(it)
4557 qs(i, k) = eps * es(i, k) * (1. + zvir * q(i, k)) / p(i, k)
4561 if (
present (dqdt))
then
4564 ap1 = 10. * dim(t(i, k), tmin) + 1.
4565 ap1 = min(2621., ap1) - 0.5
4567 dqdt(i, k) = eps10 * (des(it) + (ap1 - it) * (des(it + 1) - des(it))) * (1. + zvir * q(i, k)) / p(i, k)
4578subroutine neg_adj (ktop, kbot, pt, dp, qv, ql, qr, qi, qs, qg)
4582 integer,
intent (in) :: ktop, kbot
4584 real,
intent (in),
dimension (ktop:kbot) :: dp
4586 real,
intent (inout),
dimension (ktop:kbot) :: pt, qv, ql, qr, qi, qs, qg
4588 real,
dimension (ktop:kbot) :: lcpk, icpk
4599 cvm = c_air + qv(k) * c_vap + (qr(k) + ql(k)) * c_liq + (qi(k) + qs(k) + qg(k)) * c_ice
4600 lcpk(k) = (lv00 + d0_vap * pt(k)) / cvm
4601 icpk(k) = (li00 + dc_ice * pt(k)) / cvm
4611 if (qi(k) < 0.)
then
4612 qs(k) = qs(k) + qi(k)
4616 if (qs(k) < 0.)
then
4617 qg(k) = qg(k) + qs(k)
4621 if (qg(k) < 0.)
then
4622 qr(k) = qr(k) + qg(k)
4623 pt(k) = pt(k) - qg(k) * icpk(k)
4632 if (qr(k) < 0.)
then
4633 ql(k) = ql(k) + qr(k)
4637 if (ql(k) < 0.)
then
4638 qv(k) = qv(k) + ql(k)
4639 pt(k) = pt(k) - ql(k) * lcpk(k)
4649 do k = ktop, kbot - 1
4650 if (qv(k) < 0.)
then
4651 qv(k + 1) = qv(k + 1) + qv(k) * dp(k) / dp(k + 1)
4660 if (qv(kbot) < 0. .and. qv(kbot - 1) > 0.)
then
4661 dq = min(- qv(kbot) * dp(kbot), qv(kbot - 1) * dp(kbot - 1))
4662 qv(kbot - 1) = qv(kbot - 1) - dq / dp(kbot - 1)
4663 qv(kbot) = qv(kbot) + dq / dp(kbot)
4721 integer,
intent (in) :: is, ie, js, je, km
4723 real,
intent (in),
dimension (is:ie, js:je, km) :: a3
4725 real,
intent (in),
dimension (is:ie, js:je, km + 1) :: hgt
4727 real,
intent (in) :: zl
4729 real,
intent (out),
dimension (is:ie, js:je) :: a2
4731 real,
dimension (km) :: zm
4740 zm(k) = 0.5 * (hgt(i, j, k) + hgt(i, j, k + 1))
4742 if (zl >= zm(1))
then
4743 a2(i, j) = a3(i, j, 1)
4744 elseif (zl <= zm(km))
then
4745 a2(i, j) = a3(i, j, km)
4748 if (zl <= zm(k) .and. zl >= zm(k + 1))
then
4749 a2(i, j) = a3(i, j, k) + (a3(i, j, k + 1) - a3(i, j, k)) * (zm(k) - zl) / (zm(k) - zm(k + 1))
4765subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms, qmg, t, &
4766 rew, rei, rer, res, reg)
4770 integer,
intent (in) :: is, ie, ks, ke
4771 integer,
intent (in),
dimension (is:ie) :: lsm
4773 real,
intent (in),
dimension (is:ie, ks:ke) :: den, delp, t
4774 real,
intent (in),
dimension (is:ie, ks:ke) :: qmw, qmi, qmr, qms, qmg
4776 real,
intent (out),
dimension (is:ie, ks:ke) :: rew, rei, rer, res, reg
4778 real,
dimension (is:ie, ks:ke) :: qcw, qci, qcr, qcs, qcg
4782 real :: lambdar, lambdas, lambdag
4783 real :: dpg, rei_fac, mask, ccn, bw
4784 real,
parameter :: rho_0 = 50.e-3
4786 real :: rhow = 1.0e3, rhor = 1.0e3, rhos = 1.0e2, rhog = 4.0e2
4787 real :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6
4788 real :: alphar = 0.8, alphas = 0.25, alphag = 0.5
4789 real :: gammar = 17.837789, gammas = 8.2850630, gammag = 11.631769
4790 real :: qmin = 1.0e-12, beta = 1.22, qmin1 = 9.e-6
4795 dpg = abs(delp(i, k)) / grav
4796 mask = min(max(real(lsm(i)), 0.0), 2.0)
4802 ccn = 0.80 * (- 1.15e-3 * (ccn_o ** 2) + 0.963 * ccn_o + 5.30) * abs(mask - 1.0) + &
4803 0.67 * (- 2.10e-4 * (ccn_l ** 2) + 0.568 * ccn_l - 27.9) * (1.0 - abs(mask - 1.0))
4805 if (qmw(i, k) .gt. qmin)
then
4806 qcw(i, k) = dpg * qmw(i, k) * 1.0e3
4807 rew(i, k) = exp(1.0 / 3.0 * log((3.0 * den(i, k) * qmw(i, k)) / (4.0 * pi * rhow * ccn))) * 1.0e4
4808 rew(i, k) = max(rewmin, min(rewmax, rew(i, k)))
4814 if (reiflag .eq. 1)
then
4820 if (qmi(i, k) .gt. qmin1)
then
4821 qci(i, k) = dpg * qmi(i, k) * 1.0e3
4822 rei_fac = log(1.0e3 * qmi(i, k) * den(i, k))
4823 if (t(i, k) - tice .lt. - 50)
then
4824 rei(i, k) = beta / 9.917 * exp(0.109 * rei_fac) * 1.0e3
4825 elseif (t(i, k) - tice .lt. - 40)
then
4826 rei(i, k) = beta / 9.337 * exp(0.080 * rei_fac) * 1.0e3
4827 elseif (t(i, k) - tice .lt. - 30)
then
4828 rei(i, k) = beta / 9.208 * exp(0.055 * rei_fac) * 1.0e3
4830 rei(i, k) = beta / 9.387 * exp(0.031 * rei_fac) * 1.0e3
4832 rei(i, k) = max(reimin, min(reimax, rei(i, k)))
4840 if (reiflag .eq. 2)
then
4846 if (qmi(i, k) .gt. qmin1)
then
4847 qci(i, k) = dpg * qmi(i, k) * 1.0e3
4848 bw = - 2. + 1.e-3 * log10(den(i, k) * qmi(i, k) / rho_0) * max(0.0, tice - t(i, k)) ** 1.5
4849 rei(i, k) = 377.4 + bw * (203.3 + bw * (37.91 + 2.3696 * bw))
4850 rei(i, k) = max(reimin, min(reimax, rei(i, k)))
4862 if (qmr(i, k) .gt. qmin)
then
4863 qcr(i, k) = dpg * qmr(i, k) * 1.0e3
4864 lambdar = exp(0.25 * log(pi * rhor * n0r / qmr(i, k) / den(i, k)))
4865 rer(i, k) = 0.5 * exp(log(gammar / 6) / alphar) / lambdar * 1.0e6
4866 rer(i, k) = max(rermin, min(rermax, rer(i, k)))
4876 if (qms(i, k) .gt. qmin1)
then
4877 qcs(i, k) = dpg * qms(i, k) * 1.0e3
4878 lambdas = exp(0.25 * log(pi * rhos * n0s / qms(i, k) / den(i, k)))
4879 res(i, k) = 0.5 * exp(log(gammas / 6) / alphas) / lambdas * 1.0e6
4880 res(i, k) = max(resmin, min(resmax, res(i, k)))
4890 if (qmg(i, k) .gt. qmin)
then
4891 qcg(i, k) = dpg * qmg(i, k) * 1.0e3
4892 lambdag = exp(0.25 * log(pi * rhog * n0g / qmg(i, k) / den(i, k)))
4893 reg(i, k) = 0.5 * exp(log(gammag / 6) / alphag) / lambdag * 1.0e6
4894 reg(i, k) = max(regmin, min(regmax, reg(i, k)))
4909 t1d, p1d, dBZ, kts, kte, ii, jj, melti)
4914 INTEGER,
INTENT(IN):: kts, kte, ii,jj
4915 REAL,
DIMENSION(kts:kte),
INTENT(IN):: &
4916 qv1d, qr1d, qs1d, qg1d, t1d, p1d
4917 REAL,
DIMENSION(kts:kte),
INTENT(INOUT):: dBZ
4920 REAL,
DIMENSION(kts:kte):: temp, pres, qv, rho
4921 REAL,
DIMENSION(kts:kte):: rr, rs, rg
4924 DOUBLE PRECISION,
DIMENSION(kts:kte):: ilamr, ilams, ilamg
4925 DOUBLE PRECISION,
DIMENSION(kts:kte):: N0_r, N0_s, N0_g
4926 DOUBLE PRECISION:: lamr, lams, lamg
4927 LOGICAL,
DIMENSION(kts:kte):: L_qr, L_qs, L_qg
4929 REAL,
DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
4930 DOUBLE PRECISION:: fmelt_s, fmelt_g
4932 INTEGER:: i, k, k_0, kbot, n
4933 LOGICAL,
INTENT(IN):: melti
4934 DOUBLE PRECISION:: cback, x, eta, f_d
4947 qv(k) = max(1.e-10, qv1d(k))
4949 rho(k) = 0.622*pres(k)/(rdgas*temp(k)*(qv(k)+0.622))
4951 if (qr1d(k) .gt. 1.e-9)
then
4952 rr(k) = qr1d(k)*rho(k)
4954 lamr = (xam_r*xcrg(3)*n0_r(k)/rr(k))**(1./xcre(1))
4962 if (qs1d(k) .gt. 1.e-9)
then
4963 rs(k) = qs1d(k)*rho(k)
4965 lams = (xam_s*xcsg(3)*n0_s(k)/rs(k))**(1./xcse(1))
4973 if (qg1d(k) .gt. 1.e-9)
then
4974 rg(k) = qg1d(k)*rho(k)
4976 lamg = (xam_g*xcgg(3)*n0_g(k)/rg(k))**(1./xcge(1))
4989 k_loop:
do k = kte-1, kts, -1
4990 if ( melti .and. (temp(k).gt.273.15) .and. l_qr(k) &
4991 .and. (l_qs(k+1).or.l_qg(k+1)) )
then
5004 ze_graupel(k) = 1.e-22
5005 if (l_qr(k)) ze_rain(k) = n0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
5006 if (l_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/pi)*(6.0/pi) &
5007 * (xam_s/900.0)*(xam_s/900.0) &
5008 * n0_s(k)*xcsg(4)*ilams(k)**xcse(4)
5009 if (l_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/pi)*(6.0/pi) &
5010 * (xam_g/900.0)*(xam_g/900.0) &
5011 * n0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
5023 if (melti .and. k_0.ge.kts+1)
then
5024 do k = k_0-1, kts, -1
5027 if (l_qs(k) .and. l_qs(k_0) )
then
5028 fmelt_s = max(0.005d0, min(1.0d0-rs(k)/rs(k_0), 0.99d0))
5032 x = xam_s * xxds(n)**xbm_s
5033 call rayleigh_soak_wetgraupel (x,dble(xocms),dble(xobms), &
5034 fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
5035 cback, mixingrulestring_s, matrixstring_s, &
5036 inclusionstring_s, hoststring_s, &
5037 hostmatrixstring_s, hostinclusionstring_s)
5038 f_d = n0_s(k)*xxds(n)**xmu_s * dexp(-lams*xxds(n))
5039 eta = eta + f_d * cback * simpson(n) * xdts(n)
5041 ze_snow(k) = sngl(lamda4 / (pi5 * k_w) * eta)
5047 if (l_qg(k) .and. l_qg(k_0) )
then
5048 fmelt_g = max(0.005d0, min(1.0d0-rg(k)/rg(k_0), 0.99d0))
5052 x = xam_g * xxdg(n)**xbm_g
5053 call rayleigh_soak_wetgraupel (x,dble(xocmg),dble(xobmg), &
5054 fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
5055 cback, mixingrulestring_g, matrixstring_g, &
5056 inclusionstring_g, hoststring_g, &
5057 hostmatrixstring_g, hostinclusionstring_g)
5058 f_d = n0_g(k)*xxdg(n)**xmu_g * dexp(-lamg*xxdg(n))
5059 eta = eta + f_d * cback * simpson(n) * xdtg(n)
5061 ze_graupel(k) = sngl(lamda4 / (pi5 * k_w) * eta)
5068 dbz(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
real function acr3d(v1, v2, q1, q2, c, cac, rho)
The function is an accretion function (Lin et al.(1983) lin_et_al_1983 )
real function smlt(tc, dqs, qsrho, psacw, psacr, c, rho, rhofac)
Melting of snow function (Lin et al.(1983) lin_et_al_1983) note: psacw and psacr must be calc before ...
subroutine esw_table1d(ta, es, n)
The subroutine 'esw_table1d' computes the saturated water vapor pressure for table ii.
subroutine interpolate_z(is, ie, js, je, km, zl, hgt, a3, a2)
quick local sum algorithm
subroutine qs_table(n)
saturation water vapor pressure table i
subroutine qs_table3(n)
saturation water vapor pressure table iv
subroutine refl10cm_gfdl(qv1d, qr1d, qs1d, qg1d, t1d, p1d, dbz, kts, kte, ii, jj, melti)
This subroutine calculates radar reflectivity.
subroutine implicit_fall(dt, ktop, kbot, ze, vt, dp, q, precip, m1)
The subroutine computes the time-implicit monotonic fall scheme.
subroutine qsmith(im, km, ks, t, p, q, qs, dqdt)
The function 'qsmith' computes the saturated specific humidity with a blend of water and ice dependin...
subroutine, public cloud_diagnosis(is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms, qmg, t, rew, rei, rer, res, reg)
The subroutine 'cloud_diagnosis' diagnoses the radius of cloud species.
subroutine cs_limiters(km, a4)
This subroutine perform positive definite constraint.
subroutine linear_prof(km, q, dm, z_var, h_var)
Definition of vertical subgrid variability used for cloud ice and cloud water autoconversion.
subroutine qs_table2(n)
saturation water vapor pressure table iii
subroutine fall_speed(ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg)
The subroutine calculates vertical fall speed of snow/ice/graupel.
subroutine check_column(ktop, kbot, q, no_fall)
The subroutine 'check_column' checks if the water species is large enough to fall.
subroutine revap_racc(ktop, kbot, dt, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
This subroutine calculates evaporation of rain and accretion of rain.
subroutine icloud(ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, den, denfac, vts, vtg, vtr, qak, rh_adj, rh_rain, dts, h_var)
This subroutine includes cloud ice microphysics processes.
real function wqs1(ta, den)
The function 'wqs1' returns the saturation vapor pressure over pure liquid water for a given temperat...
subroutine revap_rac1(hydrostatic, is, ie, dt, tz, qv, ql, qr, qi, qs, qg, den, hvar)
This subroutine calculates rain evaporation.
subroutine terminal_fall(dtm, ktop, kbot, tz, qv, ql, qr, qg, qs, qi, dz, dp, den, vtg, vts, vti, r1, g1, s1, i1, m1_sol, w1)
The subroutine 'terminal_fall' computes terminal fall speed.
subroutine neg_adj(ktop, kbot, pt, dp, qv, ql, qr, qi, qs, qg)
The subroutine 'neg_adj' fixes negative water species.
subroutine mpdrv(hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, qg, qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes, rain, snow, graupel, ice, m2_rain, m2_sol, cond, area1, land, u_dt, v_dt, pt_dt, qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, w_var, vt_r, vt_s, vt_g, vt_i, qn2)
GFDL cloud microphysics, major program, and is based on Lin et al.(1983) lin_et_al_1983 and Rutledge ...
subroutine setupm
The subroutine sets up gfdl cloud microphysics parameters.
subroutine es3_table1d(ta, es, n)
The subroutine 'es3_table1d' computes the saturated water vapor pressure for table iv.
subroutine, public gfdl_cloud_microphys_mod_end()
The subroutine 'gfdl_cloud_microphys_init' terminates the GFDL cloud microphysics.
subroutine, public gfdl_cloud_microphys_mod_driver(iis, iie, jjs, jje, kks, kke, ktop, kbot, qv, ql, qr, qi, qs, qg, qa, qn, qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, pt_dt, pt, w, uin, vin, udt, vdt, dz, delp, area, dt_in, land, rain, snow, ice, graupel, hydrostatic, phys_hydrostatic, p, lradar, refl_10cm, reset, pfils, pflls)
This subroutine is the driver of the GFDL cloud microphysics.
subroutine lagrangian_fall_ppm(ktop, kbot, zs, ze, zt, dp, q, precip, m1, mono)
Lagrangian scheme.
subroutine qsmith_init
The subroutine 'qsmith_init' initializes lookup tables for saturation water vapor pressure for the fo...
subroutine, public gfdl_cloud_microphys_mod_init(me, master, nlunit, input_nml_file, logunit, fn_nml, errmsg, errflg)
The subroutine 'gfdl_cloud_microphys_init' initializes the GFDL cloud microphysics.
subroutine warm_rain(dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, den, denfac, ccn, c_praut, rh_rain, vtr, r1, m1_rain, w1, h_var)
This subroutine includes warm rain cloud microphysics.
subroutine qs_tablew(n)
saturation water vapor pressure table ii
subroutine sedi_heat(ktop, kbot, dm, m1, dz, tz, qv, ql, qr, qi, qs, qg, cw)
This subroutine calculates sedimentation of heat.
real function qs_blend(t, p, q)
The function 'qs_blend' computes the saturated specific humidity with a blend of water and ice depend...
subroutine setup_con
The subroutine 'setup_con' sets up constants and calls 'qsmith_init'.
subroutine cs_profile(a4, del, km, do_mono)
subroutine subgrid_z_proc(ktop, kbot, p1, den, denfac, dts, rh_adj, tz, qv, ql, qr, qi, qs, qg, qa, h_var, rh_rain)
This subroutine calculates temperature sentive high vertical resolution processes.
subroutine es2_table1d(ta, es, n)
The subroutine 'es3_table1d' computes the saturated water vapor pressure for table iii.
This module is more library code whereas the individual microphysics schemes contains specific detail...