39 use machine,
only: kind_phys, r8 => kind_dbl_prec
41 t_min, t_sub, tau_r2g, tau_smlt, tau_gmlt, dw_land, dw_ocean, vw_fac, vi_fac, &
42 vr_fac, vs_fac, vg_fac, ql_mlt, do_qa, fix_negative, vw_max, vi_max, vs_max, &
43 vg_max, vr_max, qs_mlt, qs0_crt, ql0_max, qi0_max, qi0_crt, ifflag, rh_inc, rh_ins,&
44 rh_inr, const_vw, const_vi, const_vs, const_vg, const_vr, rthresh, ccn_l, ccn_o, &
45 igflag, c_paut, tau_imlt, tau_v2l, tau_l2v, tau_i2s, tau_l2r, qi_lim, ql_gen, &
46 do_hail, inflag, c_psacw, c_psaci, c_pracs, c_psacr, c_pgacr, c_pgacs, c_pgacw, &
47 c_pgaci, z_slope_liq, z_slope_ice, prog_ccn, c_pracw, c_praci, rad_snow, &
48 rad_graupel, rad_rain, cld_min, sedflag, sed_fac, do_sedi_uv, do_sedi_w, &
49 do_sedi_heat, icloud_f, irain_f, xr_a, xr_b, xr_c, ntimes, tau_revp, tice_mlt, &
50 do_cond_timescale, mp_time, consv_checker, te_err, tw_err, use_rhc_cevap, &
51 use_rhc_revap, tau_wbf, do_warm_rain_mp, rh_thres, f_dq_p, f_dq_m, do_cld_adj, &
52 rhc_cevap, rhc_revap, beta, liq_ice_combine, rewflag, reiflag, rerflag, resflag, &
53 regflag, rewmin, rewmax, reimin, reimax, rermin, rermax, resmin, resmax, regmin, &
54 regmax, fs2g_fac, fi2s_fac, fi2g_fac, do_sedi_melt, radr_flag, rads_flag, &
55 radg_flag, do_wbf, do_psd_water_fall, do_psd_ice_fall, n0w_sig, n0i_sig, n0r_sig, &
56 n0s_sig, n0g_sig, n0h_sig, n0w_exp, n0i_exp, n0r_exp, n0s_exp, n0g_exp, n0h_exp, &
57 muw, mui, mur, mus, mug, muh, alinw, alini, alinr, alins, aling, alinh, blinw, &
58 blini, blinr, blins, bling, blinh, do_new_acc_water, do_new_acc_ice, is_fac, &
59 ss_fac, gs_fac, rh_fac_evap, rh_fac_cond, snow_grauple_combine, do_psd_water_num, &
60 do_psd_ice_num, vdiffflag, rewfac, reifac, cp_heating, nconds, do_evap_timescale, &
61 delay_cond_evap, do_subgrid_proc, fast_fr_mlt, fast_dep_sub, qi_gen, tice
62 use physcons,
only: grav => con_g, &
65 boltzmann => con_boltz, &
66 avogadro => con_sbc, &
70 runiver => con_runiver, &
80 con_amd, con_amw, visd, &
81 visk, vdifu, tcond, cdg, &
86 rhoh => rhoch, qcmin, qfmin
118 procedure wet_bulb_dry
119 procedure wet_bulb_moist
362 real(kind_phys) :: gcon, hcon, scm3, pisq, act (20), ace (20), occ (3), aone
368 if (do_warm_rain_mp)
then
378 fac_rc = (4. / 3.) * pi * rhow * rthresh ** 3
380 aone = 2. / 9. * (3. / 4.) ** (4. / 3.) / pi ** (1. / 3.)
381 cpaut = c_paut * aone * grav / visd
387 gcon = (4. * grav * rhog / (3. * cdg * rho0)) ** 0.5
388 hcon = (4. * grav * rhoh / (3. * cdh * rho0)) ** 0.5
394 normw = pi * rhow * n0w_sig * gamma(muw + 3)
395 normi = pi * rhoi * n0i_sig * gamma(mui + 3)
396 normr = pi * rhor * n0r_sig * gamma(mur + 3)
397 norms = pi * rhos * n0s_sig * gamma(mus + 3)
398 normg = pi * rhog * n0g_sig * gamma(mug + 3)
399 normh = pi * rhoh * n0h_sig * gamma(muh + 3)
401 expow = exp(n0w_exp / (muw + 3) * log(10.))
402 expoi = exp(n0i_exp / (mui + 3) * log(10.))
403 expor = exp(n0r_exp / (mur + 3) * log(10.))
404 expos = exp(n0s_exp / (mus + 3) * log(10.))
405 expog = exp(n0g_exp / (mug + 3) * log(10.))
406 expoh = exp(n0h_exp / (muh + 3) * log(10.))
414 pcaw = exp(3 / (muw + 3) * log(n0w_sig)) * gamma(muw) * exp(3 * n0w_exp / (muw + 3) * log(10.))
415 pcai = exp(3 / (mui + 3) * log(n0i_sig)) * gamma(mui) * exp(3 * n0i_exp / (mui + 3) * log(10.))
416 pcar = exp(3 / (mur + 3) * log(n0r_sig)) * gamma(mur) * exp(3 * n0r_exp / (mur + 3) * log(10.))
417 pcas = exp(3 / (mus + 3) * log(n0s_sig)) * gamma(mus) * exp(3 * n0s_exp / (mus + 3) * log(10.))
418 pcag = exp(3 / (mug + 3) * log(n0g_sig)) * gamma(mug) * exp(3 * n0g_exp / (mug + 3) * log(10.))
419 pcah = exp(3 / (muh + 3) * log(n0h_sig)) * gamma(muh) * exp(3 * n0h_exp / (muh + 3) * log(10.))
421 pcbw = exp(muw / (muw + 3) * log(pi * rhow * gamma(muw + 3)))
422 pcbi = exp(mui / (mui + 3) * log(pi * rhoi * gamma(mui + 3)))
423 pcbr = exp(mur / (mur + 3) * log(pi * rhor * gamma(mur + 3)))
424 pcbs = exp(mus / (mus + 3) * log(pi * rhos * gamma(mus + 3)))
425 pcbg = exp(mug / (mug + 3) * log(pi * rhog * gamma(mug + 3)))
426 pcbh = exp(muh / (muh + 3) * log(pi * rhoh * gamma(muh + 3)))
428 edaw = exp(- 1. / (muw + 3) * log(n0w_sig)) * (muw + 2) * exp(- n0w_exp / (muw + 3) * log(10.))
429 edai = exp(- 1. / (mui + 3) * log(n0i_sig)) * (mui + 2) * exp(- n0i_exp / (mui + 3) * log(10.))
430 edar = exp(- 1. / (mur + 3) * log(n0r_sig)) * (mur + 2) * exp(- n0r_exp / (mur + 3) * log(10.))
431 edas = exp(- 1. / (mus + 3) * log(n0s_sig)) * (mus + 2) * exp(- n0s_exp / (mus + 3) * log(10.))
432 edag = exp(- 1. / (mug + 3) * log(n0g_sig)) * (mug + 2) * exp(- n0g_exp / (mug + 3) * log(10.))
433 edah = exp(- 1. / (muh + 3) * log(n0h_sig)) * (muh + 2) * exp(- n0h_exp / (muh + 3) * log(10.))
435 edbw = exp(1. / (muw + 3) * log(pi * rhow * gamma(muw + 3)))
436 edbi = exp(1. / (mui + 3) * log(pi * rhoi * gamma(mui + 3)))
437 edbr = exp(1. / (mur + 3) * log(pi * rhor * gamma(mur + 3)))
438 edbs = exp(1. / (mus + 3) * log(pi * rhos * gamma(mus + 3)))
439 edbg = exp(1. / (mug + 3) * log(pi * rhog * gamma(mug + 3)))
440 edbh = exp(1. / (muh + 3) * log(pi * rhoh * gamma(muh + 3)))
442 oeaw = exp(1. / (muw + 3) * log(n0w_sig)) * pi * gamma(muw + 2) * &
443 exp(n0w_exp / (muw + 3) * log(10.))
444 oeai = exp(1. / (mui + 3) * log(n0i_sig)) * pi * gamma(mui + 2) * &
445 exp(n0i_exp / (mui + 3) * log(10.))
446 oear = exp(1. / (mur + 3) * log(n0r_sig)) * pi * gamma(mur + 2) * &
447 exp(n0r_exp / (mur + 3) * log(10.))
448 oeas = exp(1. / (mus + 3) * log(n0s_sig)) * pi * gamma(mus + 2) * &
449 exp(n0s_exp / (mus + 3) * log(10.))
450 oeag = exp(1. / (mug + 3) * log(n0g_sig)) * pi * gamma(mug + 2) * &
451 exp(n0g_exp / (mug + 3) * log(10.))
452 oeah = exp(1. / (muh + 3) * log(n0h_sig)) * pi * gamma(muh + 2) * &
453 exp(n0h_exp / (muh + 3) * log(10.))
455 oebw = 2 * exp((muw + 2) / (muw + 3) * log(pi * rhow * gamma(muw + 3)))
456 oebi = 2 * exp((mui + 2) / (mui + 3) * log(pi * rhoi * gamma(mui + 3)))
457 oebr = 2 * exp((mur + 2) / (mur + 3) * log(pi * rhor * gamma(mur + 3)))
458 oebs = 2 * exp((mus + 2) / (mus + 3) * log(pi * rhos * gamma(mus + 3)))
459 oebg = 2 * exp((mug + 2) / (mug + 3) * log(pi * rhog * gamma(mug + 3)))
460 oebh = 2 * exp((muh + 2) / (muh + 3) * log(pi * rhoh * gamma(muh + 3)))
462 rraw = exp(- 3 / (muw + 3) * log(n0w_sig)) * gamma(muw + 6) * &
463 exp(- 3 * n0w_exp / (muw + 3) * log(10.))
464 rrai = exp(- 3 / (mui + 3) * log(n0i_sig)) * gamma(mui + 6) * &
465 exp(- 3 * n0i_exp / (mui + 3) * log(10.))
466 rrar = exp(- 3 / (mur + 3) * log(n0r_sig)) * gamma(mur + 6) * &
467 exp(- 3 * n0r_exp / (mur + 3) * log(10.))
468 rras = exp(- 3 / (mus + 3) * log(n0s_sig)) * gamma(mus + 6) * &
469 exp(- 3 * n0s_exp / (mus + 3) * log(10.))
470 rrag = exp(- 3 / (mug + 3) * log(n0g_sig)) * gamma(mug + 6) * &
471 exp(- 3 * n0g_exp / (mug + 3) * log(10.))
472 rrah = exp(- 3 / (muh + 3) * log(n0h_sig)) * gamma(muh + 6) * &
473 exp(- 3 * n0h_exp / (muh + 3) * log(10.))
475 rrbw = exp((muw + 6) / (muw + 3) * log(pi * rhow * gamma(muw + 3)))
476 rrbi = exp((mui + 6) / (mui + 3) * log(pi * rhoi * gamma(mui + 3)))
477 rrbr = exp((mur + 6) / (mur + 3) * log(pi * rhor * gamma(mur + 3)))
478 rrbs = exp((mus + 6) / (mus + 3) * log(pi * rhos * gamma(mus + 3)))
479 rrbg = exp((mug + 6) / (mug + 3) * log(pi * rhog * gamma(mug + 3)))
480 rrbh = exp((muh + 6) / (muh + 3) * log(pi * rhoh * gamma(muh + 3)))
482 tvaw = exp(- blinw / (muw + 3) * log(n0w_sig)) * alinw * gamma(muw + blinw + 3) * &
483 exp(- blinw * n0w_exp / (muw + 3) * log(10.))
484 tvai = exp(- blini / (mui + 3) * log(n0i_sig)) * alini * gamma(mui + blini + 3) * &
485 exp(- blini * n0i_exp / (mui + 3) * log(10.))
486 tvar = exp(- blinr / (mur + 3) * log(n0r_sig)) * alinr * gamma(mur + blinr + 3) * &
487 exp(- blinr * n0r_exp / (mur + 3) * log(10.))
488 tvas = exp(- blins / (mus + 3) * log(n0s_sig)) * alins * gamma(mus + blins + 3) * &
489 exp(- blins * n0s_exp / (mus + 3) * log(10.))
490 tvag = exp(- bling / (mug + 3) * log(n0g_sig)) * aling * gamma(mug + bling + 3) * &
491 exp(- bling * n0g_exp / (mug + 3) * log(10.)) * gcon
492 tvah = exp(- blinh / (muh + 3) * log(n0h_sig)) * alinh * gamma(muh + blinh + 3) * &
493 exp(- blinh * n0h_exp / (muh + 3) * log(10.)) * hcon
495 tvbw = exp(blinw / (muw + 3) * log(pi * rhow * gamma(muw + 3))) * gamma(muw + 3)
496 tvbi = exp(blini / (mui + 3) * log(pi * rhoi * gamma(mui + 3))) * gamma(mui + 3)
497 tvbr = exp(blinr / (mur + 3) * log(pi * rhor * gamma(mur + 3))) * gamma(mur + 3)
498 tvbs = exp(blins / (mus + 3) * log(pi * rhos * gamma(mus + 3))) * gamma(mus + 3)
499 tvbg = exp(bling / (mug + 3) * log(pi * rhog * gamma(mug + 3))) * gamma(mug + 3)
500 tvbh = exp(blinh / (muh + 3) * log(pi * rhoh * gamma(muh + 3))) * gamma(muh + 3)
506 scm3 = exp(1. / 3. * log(visk / vdifu))
514 cracw = pi * n0r_sig * alinr * gamma(2 + mur + blinr) / &
515 (4. * exp((2 + mur + blinr) / (mur + 3) * log(normr))) * &
516 exp((1 - blinr) * log(expor))
517 craci = pi * n0r_sig * alinr * gamma(2 + mur + blinr) / &
518 (4. * exp((2 + mur + blinr) / (mur + 3) * log(normr))) * &
519 exp((1 - blinr) * log(expor))
520 csacw = pi * n0s_sig * alins * gamma(2 + mus + blins) / &
521 (4. * exp((2 + mus + blins) / (mus + 3) * log(norms))) * &
522 exp((1 - blins) * log(expos))
523 csaci = pi * n0s_sig * alins * gamma(2 + mus + blins) / &
524 (4. * exp((2 + mus + blins) / (mus + 3) * log(norms))) * &
525 exp((1 - blins) * log(expos))
527 cgacw = pi * n0h_sig * alinh * gamma(2 + muh + blinh) * hcon / &
528 (4. * exp((2 + muh + blinh) / (muh + 3) * log(normh))) * &
529 exp((1 - blinh) * log(expoh))
530 cgaci = pi * n0h_sig * alinh * gamma(2 + muh + blinh) * hcon / &
531 (4. * exp((2 + muh + blinh) / (muh + 3) * log(normh))) * &
532 exp((1 - blinh) * log(expoh))
534 cgacw = pi * n0g_sig * aling * gamma(2 + mug + bling) * gcon / &
535 (4. * exp((2 + mug + bling) / (mug + 3) * log(normg))) * &
536 exp((1 - bling) * log(expog))
537 cgaci = pi * n0g_sig * aling * gamma(2 + mug + bling) * gcon / &
538 (4. * exp((2 + mug + bling) / (mug + 3) * log(normg))) * &
539 exp((1 - bling) * log(expog))
542 if (do_new_acc_water)
then
544 cracw = pisq * n0r_sig * n0w_sig * rhow / 24.
545 csacw = pisq * n0s_sig * n0w_sig * rhow / 24.
547 cgacw = pisq * n0h_sig * n0w_sig * rhow / 24.
549 cgacw = pisq * n0g_sig * n0w_sig * rhow / 24.
554 if (do_new_acc_ice)
then
556 craci = pisq * n0r_sig * n0i_sig * rhoi / 24.
557 csaci = pisq * n0s_sig * n0i_sig * rhoi / 24.
559 cgaci = pisq * n0h_sig * n0i_sig * rhoi / 24.
561 cgaci = pisq * n0g_sig * n0i_sig * rhoi / 24.
566 cracw = cracw * c_pracw
567 craci = craci * c_praci
568 csacw = csacw * c_psacw
569 csaci = csaci * c_psaci
570 cgacw = cgacw * c_pgacw
571 cgaci = cgaci * c_pgaci
577 cracs = pisq * n0r_sig * n0s_sig * rhos / 24.
578 csacr = pisq * n0s_sig * n0r_sig * rhor / 24.
580 cgacr = pisq * n0h_sig * n0r_sig * rhor / 24.
581 cgacs = pisq * n0h_sig * n0s_sig * rhos / 24.
583 cgacr = pisq * n0g_sig * n0r_sig * rhor / 24.
584 cgacs = pisq * n0g_sig * n0s_sig * rhos / 24.
587 cracs = cracs * c_pracs
588 csacr = csacr * c_psacr
589 cgacr = cgacr * c_pgacr
590 cgacs = cgacs * c_pgacs
685 acco(i, k) = occ(i) * gamma(6 + acc(2 * k - 1) - i) * gamma(acc(2 * k) + i - 1) / &
686 (exp((6 + acc(2 * k - 1) - i) / (acc(2 * k - 1) + 3) * log(act(2 * k - 1))) * &
687 exp((acc(2 * k) + i - 1) / (acc(2 * k) + 3) * log(act(2 * k)))) * &
688 exp((i - 3) * log(ace(2 * k - 1))) * exp((4 - i) * log(ace(2 * k)))
696 crevp(1) = 2. * pi * vdifu * tcond * rvgas * n0r_sig * gamma(1 + mur) / &
697 exp((1 + mur) / (mur + 3) * log(normr)) * exp(2.0 * log(expor))
699 crevp(3) = 0.31 * scm3 * sqrt(alinr / visk) * gamma((3 + 2 * mur + blinr) / 2) / &
700 exp((3 + 2 * mur + blinr) / (mur + 3) / 2 * log(normr)) * &
701 exp((1 + mur) / (mur + 3) * log(normr)) / gamma(1 + mur) * &
702 exp((- 1 - blinr) / 2. * log(expor))
703 crevp(4) = tcond * rvgas
706 cssub(1) = 2. * pi * vdifu * tcond * rvgas * n0s_sig * gamma(1 + mus) / &
707 exp((1 + mus) / (mus + 3) * log(norms)) * exp(2.0 * log(expos))
709 cssub(3) = 0.31 * scm3 * sqrt(alins / visk) * gamma((3 + 2 * mus + blins) / 2) / &
710 exp((3 + 2 * mus + blins) / (mus + 3) / 2 * log(norms)) * &
711 exp((1 + mus) / (mus + 3) * log(norms)) / gamma(1 + mus) * &
712 exp((- 1 - blins) / 2. * log(expos))
713 cssub(4) = tcond * rvgas
717 cgsub(1) = 2. * pi * vdifu * tcond * rvgas * n0h_sig * gamma(1 + muh) / &
718 exp((1 + muh) / (muh + 3) * log(normh)) * exp(2.0 * log(expoh))
720 cgsub(3) = 0.31 * scm3 * sqrt(alinh * hcon / visk) * gamma((3 + 2 * muh + blinh) / 2) / &
721 exp(1. / (muh + 3) * (3 + 2 * muh + blinh) / 2 * log(normh)) * &
722 exp(1. / (muh + 3) * (1 + muh) * log(normh)) / gamma(1 + muh) * &
723 exp((- 1 - blinh) / 2. * log(expoh))
725 cgsub(1) = 2. * pi * vdifu * tcond * rvgas * n0g_sig * gamma(1 + mug) / &
726 exp((1 + mug) / (mug + 3) * log(normg)) * exp(2.0 * log(expog))
728 cgsub(3) = 0.31 * scm3 * sqrt(aling * gcon / visk) * gamma((3 + 2 * mug + bling) / 2) / &
729 exp((3 + 2 * mug + bling) / (mug + 3) / 2 * log(normg)) * &
730 exp((1 + mug) / (mug + 3) * log(normg)) / gamma(1 + mug) * &
731 exp((- 1 - bling) / 2. * log(expog))
733 cgsub(4) = tcond * rvgas
740 csmlt(1) = 2. * pi * tcond * n0s_sig * gamma(1 + mus) / &
741 exp((1 + mus) / (mus + 3) * log(norms)) * exp(2.0 * log(expos))
742 csmlt(2) = 2. * pi * vdifu * n0s_sig * gamma(1 + mus) / &
743 exp((1 + mus) / (mus + 3) * log(norms)) * exp(2.0 * log(expos))
752 cgmlt(1) = 2. * pi * tcond * n0h_sig * gamma(1 + muh) / &
753 exp((1 + muh) / (muh + 3) * log(normh)) * exp(2.0 * log(expoh))
754 cgmlt(2) = 2. * pi * vdifu * n0h_sig * gamma(1 + muh) / &
755 exp((1 + muh) / (muh + 3) * log(normh)) * exp(2.0 * log(expoh))
757 cgmlt(1) = 2. * pi * tcond * n0g_sig * gamma(1 + mug) / &
758 exp((1 + mug) / (mug + 3) * log(normg)) * exp(2.0 * log(expog))
759 cgmlt(2) = 2. * pi * vdifu * n0g_sig * gamma(1 + mug) / &
760 exp((1 + mug) / (mug + 3) * log(normg)) * exp(2.0 * log(expog))
769 cgfr(1) = 1.e2 / 36 * pisq * n0r_sig * rhor * gamma(6 + mur) / &
770 exp((6 + mur) / (mur + 3) * log(normr)) * exp(- 3.0 * log(expor))
818subroutine mpdrv (hydrostatic, ua, va, wa, delp, pt, qv, ql, qr, qi, qs, qg, &
819 qa, qnl, qni, delz, is, ie, ks, ke, dtm, water, rain, ice, snow, graupel, &
820 gsize, hs, q_con, cappa, consv_te, adj_vmr, te, dte, prefluxw, prefluxr, &
821 prefluxi, prefluxs, prefluxg, last_step, do_inline_mp, do_mp_fast, do_mp_full)
829 integer,
intent (in) :: is, ie, ks, ke
831 logical,
intent (in) :: hydrostatic, last_step, consv_te, do_inline_mp
832 logical,
intent (in) :: do_mp_fast, do_mp_full
834 real(kind_phys),
intent (in) :: dtm
836 real(kind_phys),
intent (in),
dimension (is:ie) :: gsize, hs
838 real(kind_phys),
intent (in),
dimension (is:ie, ks:ke) :: qnl, qni
840 real(kind_phys),
intent (inout),
dimension (is:ie, ks:ke) :: delp, delz, pt, ua, va, wa
841 real(kind_phys),
intent (inout),
dimension (is:ie, ks:ke) :: qv, ql, qr, qi, qs, qg, qa
842 real(kind_phys),
intent (inout),
dimension (is:ie, ks:ke) :: prefluxw, prefluxr, prefluxi, prefluxs, prefluxg
844 real(kind_phys),
intent (inout),
dimension (is:, ks:) :: q_con, cappa
846 real(kind_phys),
intent (inout),
dimension (is:ie) :: water, rain, ice, snow, graupel
848 real(kind_phys),
intent (out),
dimension (is:ie, ks:ke) :: te, adj_vmr
850 real (kind = r8),
intent (out),
dimension (is:ie) :: dte
858 real(kind_phys) :: rh_adj, rh_rain, ccn0, cin0, cond, q1, q2
859 real(kind_phys) :: convt, dts, q_cond, t_lnd, t_ocn, h_var, tmp, nl, ni
861 real(kind_phys),
dimension (ks:ke) :: q_liq, q_sol, dp, dz, dp0
862 real(kind_phys),
dimension (ks:ke) :: qvz, qlz, qrz, qiz, qsz, qgz, qaz
863 real(kind_phys),
dimension (ks:ke) :: den, pz, denfac, ccn, cin
864 real(kind_phys),
dimension (ks:ke) :: u, v, w
866 real(kind_phys),
dimension (is:ie, ks:ke) :: pcw, edw, oew, rrw, tvw
867 real(kind_phys),
dimension (is:ie, ks:ke) :: pci, edi, oei, rri, tvi
868 real(kind_phys),
dimension (is:ie, ks:ke) :: pcr, edr, oer, rrr, tvr
869 real(kind_phys),
dimension (is:ie, ks:ke) :: pcs, eds, oes, rrs, tvs
870 real(kind_phys),
dimension (is:ie, ks:ke) :: pcg, edg, oeg, rrg, tvg
872 real(kind_phys),
dimension (is:ie) :: condensation, deposition
873 real(kind_phys),
dimension (is:ie) :: evaporation, sublimation
875 real (kind = r8) :: con_r8, c8, cp8
877 real (kind = r8),
dimension (is:ie, ks:ke) :: te_beg_d, te_end_d, tw_beg_d, tw_end_d
878 real (kind = r8),
dimension (is:ie, ks:ke) :: te_beg_m, te_end_m, tw_beg_m, tw_end_m
880 real (kind = r8),
dimension (is:ie) :: te_b_beg_d, te_b_end_d, tw_b_beg_d, tw_b_end_d, te_loss
881 real (kind = r8),
dimension (is:ie) :: te_b_beg_m, te_b_end_m, tw_b_beg_m, tw_b_end_m
883 real (kind = r8),
dimension (ks:ke) :: tz, tzuv, tzw
889 ntimes = max(ntimes, int(dtm / min(dtm, mp_time)))
890 dts = dtm / real(ntimes, kind=kind_phys)
909 convt = 86400. * rgrav / dtm
917 if (do_inline_mp)
then
919 q_cond = ql(i, k) + qr(i, k) + qi(i, k) + qs(i, k) + qg(i, k)
920 tz(k) = pt(i, k) / ((1. + zvir * qv(i, k)) * (1. - q_cond))
933 if (hydrostatic)
then
935 te(i, k) = - c_air * tz(k) * delp(i, k)
939 te(i, k) = - mte(qv(i, k), ql(i, k), qr(i, k), qi(i, k), &
940 qs(i, k), qg(i, k), tz(k), delp(i, k), .true.) * grav
949 if (consv_checker)
then
950 call mtetw (ks, ke, qv(i, :), ql(i, :), qr(i, :), qi(i, :), &
951 qs(i, :), qg(i, :), tz, ua(i, :), va(i, :), wa(i, :), &
952 delp(i, :), dte(i), 0.0, water(i), rain(i), ice(i), &
953 snow(i), graupel(i), 0.0, 0.0, dtm, te_beg_m(i, :), &
954 tw_beg_m(i, :), te_b_beg_m(i), tw_b_beg_m(i), .true., hydrostatic)
971 if (do_inline_mp)
then
972 q_cond = qlz(k) + qrz(k) + qiz(k) + qsz(k) + qgz(k)
973 con_r8 = one_r8 - (qvz(k) + q_cond)
975 con_r8 = one_r8 - qvz(k)
979 dp(k) = delp(i, k) * con_r8
980 con_r8 = one_r8 / con_r8
981 qvz(k) = qvz(k) * con_r8
982 qlz(k) = qlz(k) * con_r8
983 qrz(k) = qrz(k) * con_r8
984 qiz(k) = qiz(k) * con_r8
985 qsz(k) = qsz(k) * con_r8
986 qgz(k) = qgz(k) * con_r8
993 den(k) = - dp(k) / (grav * dz(k))
994 pz(k) = den(k) * rdgas * tz(k)
1002 if (.not. hydrostatic)
then
1009 denfac(k) = sqrt(den(ke) / den(k))
1016 if (consv_checker)
then
1017 call mtetw (ks, ke, qvz, qlz, qrz, qiz, qsz, qgz, tz, u, v, w, &
1018 dp, dte(i), 0.0, water(i), rain(i), ice(i), snow(i), &
1019 graupel(i), 0.0, 0.0, dtm, te_beg_d(i, :), tw_beg_d(i, :), &
1020 te_b_beg_d(i), tw_b_beg_d(i), .false., hydrostatic)
1030 nl = min(1., abs(hs(i)) / (10. * grav)) * &
1031 (10. ** 2.24 * (qnl(i, k) * den(k) * 1.e9) ** 0.257) + &
1032 (1. - min(1., abs(hs(i)) / (10. * grav))) * &
1033 (10. ** 2.06 * (qnl(i, k) * den(k) * 1.e9) ** 0.48)
1035 ccn(k) = max(10.0, nl) * 1.e6
1036 cin(k) = max(10.0, ni) * 1.e6
1037 ccn(k) = ccn(k) / den(k)
1038 cin(k) = cin(k) / den(k)
1041 ccn0 = (ccn_l * min(1., abs(hs(i)) / (10. * grav)) + &
1042 ccn_o * (1. - min(1., abs(hs(i)) / (10. * grav)))) * 1.e6
1045 ccn(k) = ccn0 / den(k)
1046 cin(k) = cin0 / den(k)
1055 t_lnd = dw_land * sqrt(gsize(i) / 1.e5)
1056 t_ocn = dw_ocean * sqrt(gsize(i) / 1.e5)
1057 tmp = min(1., abs(hs(i)) / (10. * grav))
1058 h_var = t_lnd * tmp + t_ocn * (1. - tmp)
1059 h_var = min(0.20, max(0.01, h_var))
1065 rh_adj = 1. - h_var - rh_inc
1066 rh_rain = max(0.35, rh_adj - rh_inr)
1073 call neg_adj (ks, ke, tz, dp, qvz, qlz, qrz, qiz, qsz, qgz, cond)
1075 condensation(i) = condensation(i) + cond * convt
1081 if (do_mp_fast)
then
1083 call mp_fast (ks, ke, tz, qvz, qlz, qrz, qiz, qsz, qgz, dtm, dp, den, &
1084 ccn, cin, condensation(i), deposition(i), evaporation(i), &
1085 sublimation(i), denfac, convt, last_step)
1093 if (do_mp_full)
then
1095 call mp_full (ks, ke, ntimes, tz, qvz, qlz, qrz, qiz, qsz, qgz, dp, dz, &
1096 u, v, w, den, denfac, ccn, cin, dts, rh_adj, rh_rain, h_var, dte(i), &
1097 water(i), rain(i), ice(i), snow(i), graupel(i), prefluxw(i, :), &
1098 prefluxr(i, :), prefluxi(i, :), prefluxs(i, :), prefluxg(i, :), &
1099 condensation(i), deposition(i), evaporation(i), sublimation(i), &
1108 if (do_qa .and. last_step)
then
1109 call cloud_fraction (ks, ke, pz, den, qvz, qlz, qrz, qiz, qsz, qgz, qaz, &
1110 tz, h_var, gsize(i))
1146 if (qlz(k) .gt. qcmin)
then
1147 call cal_pc_ed_oe_rr_tv (qlz(k), den(k), blinw, muw, pcaw, pcbw, pcw(i, k), &
1148 edaw, edbw, edw(i, k), oeaw, oebw, oew(i, k), rraw, rrbw, rrw(i, k), &
1149 tvaw, tvbw, tvw(i, k))
1151 if (qiz(k) .gt. qcmin)
then
1152 call cal_pc_ed_oe_rr_tv (qiz(k), den(k), blini, mui, pcai, pcbi, pci(i, k), &
1153 edai, edbi, edi(i, k), oeai, oebi, oei(i, k), rrai, rrbi, rri(i, k), &
1154 tvai, tvbi, tvi(i, k))
1156 if (qrz(k) .gt. qcmin)
then
1157 call cal_pc_ed_oe_rr_tv (qrz(k), den(k), blinr, mur, pcar, pcbr, pcr(i, k), &
1158 edar, edbr, edr(i, k), oear, oebr, oer(i, k), rrar, rrbr, rrr(i, k), &
1159 tvar, tvbr, tvr(i, k))
1161 if (qsz(k) .gt. qcmin)
then
1162 call cal_pc_ed_oe_rr_tv (qsz(k), den(k), blins, mus, pcas, pcbs, pcs(i, k), &
1163 edas, edbs, eds(i, k), oeas, oebs, oes(i, k), rras, rrbs, rrs(i, k), &
1164 tvas, tvbs, tvs(i, k))
1167 if (qgz(k) .gt. qcmin)
then
1168 call cal_pc_ed_oe_rr_tv (qgz(k), den(k), blinh, muh, pcah, pcbh, pcg(i, k), &
1169 edah, edbh, edg(i, k), oeah, oebh, oeg(i, k), rrah, rrbh, rrg(i, k), &
1170 tvah, tvbh, tvg(i, k))
1173 if (qgz(k) .gt. qcmin)
then
1174 call cal_pc_ed_oe_rr_tv (qgz(k), den(k), bling, mug, pcag, pcbg, pcg(i, k), &
1175 edag, edbg, edg(i, k), oeag, oebg, oeg(i, k), rrag, rrbg, rrg(i, k), &
1176 tvag, tvbg, tvg(i, k))
1186 if (do_sedi_uv)
then
1188 c8 =
mhc(qvz(k), qlz(k), qrz(k), qiz(k), qsz(k), qgz(k)) * c_air
1189 tzuv(k) = 0.5 * (ua(i, k) ** 2 + va(i, k) ** 2 - (u(k) ** 2 + v(k) ** 2)) / c8
1190 tz(k) = tz(k) + tzuv(k)
1196 c8 =
mhc(qvz(k), qlz(k), qrz(k), qiz(k), qsz(k), qgz(k)) * c_air
1197 tzw(k) = 0.5 * (wa(i, k) ** 2 - w(k) ** 2) / c8
1198 tz(k) = tz(k) + tzw(k)
1206 if (consv_checker)
then
1207 call mtetw (ks, ke, qvz, qlz, qrz, qiz, qsz, qgz, tz, u, v, w, &
1208 dp, dte(i), 0.0, water(i), rain(i), ice(i), snow(i), &
1209 graupel(i), 0.0, 0.0, dtm, te_end_d(i, :), tw_end_d(i, :), &
1210 te_b_end_d(i), tw_b_end_d(i), .false., hydrostatic, te_loss(i))
1219 if (do_inline_mp)
then
1220 q_cond = qlz(k) + qrz(k) + qiz(k) + qsz(k) + qgz(k)
1221 con_r8 = one_r8 + qvz(k) + q_cond
1223 con_r8 = one_r8 + qvz(k)
1226 delp(i, k) = dp(k) * con_r8
1227 con_r8 = one_r8 / con_r8
1228 qvz(k) = qvz(k) * con_r8
1229 qlz(k) = qlz(k) * con_r8
1230 qrz(k) = qrz(k) * con_r8
1231 qiz(k) = qiz(k) * con_r8
1232 qsz(k) = qsz(k) * con_r8
1233 qgz(k) = qgz(k) * con_r8
1235 q1 = qv(i, k) + ql(i, k) + qr(i, k) + qi(i, k) + qs(i, k) + qg(i, k)
1236 q2 = qvz(k) + qlz(k) + qrz(k) + qiz(k) + qsz(k) + qgz(k)
1237 adj_vmr(i, k) = ((one_r8 - q1) / (one_r8 - q2)) / (one_r8 + q2 - q1)
1251 q_liq(k) = qlz(k) + qrz(k)
1252 q_sol(k) = qiz(k) + qsz(k) + qgz(k)
1253 q_cond = q_liq(k) + q_sol(k)
1254 con_r8 = one_r8 - (qvz(k) + q_cond)
1255 c8 =
mhc(con_r8, qvz(k), q_liq(k), q_sol(k)) * c_air
1258 q_con(i, k) = q_cond
1261 tmp = rdgas * (1. + zvir * qvz(k))
1262 cappa(i, k) = tmp / (tmp + c8)
1272 if (do_sedi_uv)
then
1274 tz(k) = tz(k) - tzuv(k)
1275 q_liq(k) = qlz(k) + qrz(k)
1276 q_sol(k) = qiz(k) + qsz(k) + qgz(k)
1277 q_cond = q_liq(k) + q_sol(k)
1278 con_r8 = one_r8 - (qvz(k) + q_cond)
1279 c8 =
mhc(con_r8, qvz(k), q_liq(k), q_sol(k)) * c_air
1280 tzuv(k) = (0.5 * (ua(i, k) ** 2 + va(i, k) ** 2) * dp0(k) - &
1281 0.5 * (u(k) ** 2 + v(k) ** 2) * delp(i, k)) / c8 / delp(i, k)
1282 tz(k) = tz(k) + tzuv(k)
1292 tz(k) = tz(k) - tzw(k)
1293 q_liq(k) = qlz(k) + qrz(k)
1294 q_sol(k) = qiz(k) + qsz(k) + qgz(k)
1295 q_cond = q_liq(k) + q_sol(k)
1296 con_r8 = one_r8 - (qvz(k) + q_cond)
1297 c8 =
mhc(con_r8, qvz(k), q_liq(k), q_sol(k)) * c_air
1298 tzw(k) = (0.5 * (wa(i, k) ** 2) * dp0(k) - &
1299 0.5 * (w(k) ** 2) * delp(i, k)) / c8 / delp(i, k)
1300 tz(k) = tz(k) + tzw(k)
1311 if (consv_checker)
then
1312 call mtetw (ks, ke, qv(i, :), ql(i, :), qr(i, :), qi(i, :), &
1313 qs(i, :), qg(i, :), tz, ua(i, :), va(i, :), wa(i, :), &
1314 delp(i, :), dte(i), 0.0, water(i), rain(i), ice(i), &
1315 snow(i), graupel(i), 0.0, 0.0, dtm, te_end_m(i, :), &
1316 tw_end_m(i, :), te_b_end_m(i), tw_b_end_m(i), .true., hydrostatic)
1324 if (hydrostatic)
then
1326 te(i, k) = te(i, k) + c_air * tz(k) * delp(i, k)
1330 te(i, k) = te(i, k) + mte(qv(i, k), ql(i, k), qr(i, k), qi(i, k), &
1331 qs(i, k), qg(i, k), tz(k), delp(i, k), .true.) * grav
1340 if (do_inline_mp)
then
1342 q_cond = qlz(k) + qrz(k) + qiz(k) + qsz(k) + qgz(k)
1343 if (cp_heating)
then
1344 con_r8 = one_r8 - (qvz(k) + q_cond)
1345 c8 =
mhc(con_r8, qvz(k), q_liq(k), q_sol(k)) * c_air
1346 cp8 = con_r8 * cp_air + qvz(k) * cp_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1347 delz(i, k) = delz(i, k) / pt(i, k)
1348 pt(i, k) = pt(i, k) + (tz(k) * ((1. + zvir * qvz(k)) * (1. - q_cond)) - pt(i, k)) * c8 / cp8
1349 delz(i, k) = delz(i, k) * pt(i, k)
1351 pt(i, k) = tz(k) * ((1. + zvir * qvz(k)) * (1. - q_cond))
1356 q_liq(k) = qlz(k) + qrz(k)
1357 q_sol(k) = qiz(k) + qsz(k) + qgz(k)
1358 q_cond = q_liq(k) + q_sol(k)
1359 con_r8 = one_r8 - (qvz(k) + q_cond)
1360 c8 =
mhc(con_r8, qvz(k), q_liq(k), q_sol(k)) * c_air
1361 pt(i, k) = pt(i, k) + (tz(k) - pt(i, k)) * c8 / cp_air
1369 if (consv_checker)
then
1370 if (abs(sum(te_end_d(i, :)) + te_b_end_d(i) - sum(te_beg_d(i, :)) - te_b_beg_d(i)) / &
1371 (sum(te_beg_d(i, :)) + te_b_beg_d(i)) .gt. te_err)
then
1372 print*,
"GFDL-MP-DRY TE: ", &
1375 (sum(te_end_d(i, :)) + te_b_end_d(i) - sum(te_beg_d(i, :)) - te_b_beg_d(i)) / &
1376 (sum(te_beg_d(i, :)) + te_b_beg_d(i))
1378 if (abs(sum(tw_end_d(i, :)) + tw_b_end_d(i) - sum(tw_beg_d(i, :)) - tw_b_beg_d(i)) / &
1379 (sum(tw_beg_d(i, :)) + tw_b_beg_d(i)) .gt. tw_err)
then
1380 print*,
"GFDL-MP-DRY TW: ", &
1383 (sum(tw_end_d(i, :)) + tw_b_end_d(i) - sum(tw_beg_d(i, :)) - tw_b_beg_d(i)) / &
1384 (sum(tw_beg_d(i, :)) + tw_b_beg_d(i))
1387 if (abs(sum(te_end_m(i, :)) + te_b_end_m(i) - sum(te_beg_m(i, :)) - te_b_beg_m(i)) / &
1388 (sum(te_beg_m(i, :)) + te_b_beg_m(i)) .gt. te_err)
then
1389 print*,
"GFDL-MP-WET TE: ", &
1392 (sum(te_end_m(i, :)) + te_b_end_m(i) - sum(te_beg_m(i, :)) - te_b_beg_m(i)) / &
1393 (sum(te_beg_m(i, :)) + te_b_beg_m(i))
1395 if (abs(sum(tw_end_m(i, :)) + tw_b_end_m(i) - sum(tw_beg_m(i, :)) - tw_b_beg_m(i)) / &
1396 (sum(tw_beg_m(i, :)) + tw_b_beg_m(i)) .gt. tw_err)
then
1397 print*,
"GFDL-MP-WET TW: ", &
1400 (sum(tw_end_m(i, :)) + tw_b_end_m(i) - sum(tw_beg_m(i, :)) - tw_b_beg_m(i)) / &
1401 (sum(tw_beg_m(i, :)) + tw_b_beg_m(i))
1532subroutine mp_full (ks, ke, ntimes, tz, qv, ql, qr, qi, qs, qg, dp, dz, u, v, w, &
1533 den, denfac, ccn, cin, dts, rh_adj, rh_rain, h_var, dte, water, rain, ice, &
1534 snow, graupel, prefluxw, prefluxr, prefluxi, prefluxs, prefluxg, &
1535 condensation, deposition, evaporation, sublimation, convt, last_step)
1543 logical,
intent (in) :: last_step
1545 integer,
intent (in) :: ks, ke, ntimes
1547 real(kind_phys),
intent (in) :: dts, rh_adj, rh_rain, h_var, convt
1549 real(kind_phys),
intent (in),
dimension (ks:ke) :: dp, dz, den, denfac
1551 real(kind_phys),
intent (inout),
dimension (ks:ke) :: qv, ql, qr, qi, qs, qg, u, v, w, ccn, cin
1552 real(kind_phys),
intent (inout),
dimension (ks:ke) :: prefluxw, prefluxr, prefluxi, prefluxs, prefluxg
1554 real (kind = r8),
intent (inout),
dimension (ks:ke) :: tz
1556 real(kind_phys),
intent (inout) :: water, rain, ice, snow, graupel
1557 real(kind_phys),
intent (inout) :: condensation, deposition
1558 real(kind_phys),
intent (inout) :: evaporation, sublimation
1560 real (kind = r8),
intent (inout) :: dte
1568 real(kind_phys) :: w1, r1, i1, s1, g1, cond, dep, reevap, sub
1570 real(kind_phys),
dimension (ks:ke) :: vtw, vtr, vti, vts, vtg, pfw, pfr, pfi, pfs, pfg
1578 call sedimentation (dts, ks, ke, tz, qv, ql, qr, qi, qs, qg, &
1579 dz, dp, vtw, vtr, vti, vts, vtg, w1, r1, i1, s1, g1, pfw, pfr, pfi, pfs, pfg, &
1580 u, v, w, den, denfac, dte)
1582 water = water + w1 * convt
1583 rain = rain + r1 * convt
1584 ice = ice + i1 * convt
1585 snow = snow + s1 * convt
1586 graupel = graupel + g1 * convt
1588 prefluxw = prefluxw + pfw * convt
1589 prefluxr = prefluxr + pfr * convt
1590 prefluxi = prefluxi + pfi * convt
1591 prefluxs = prefluxs + pfs * convt
1592 prefluxg = prefluxg + pfg * convt
1598 call warm_rain (dts, ks, ke, dp, dz, tz, qv, ql, qr, qi, qs, qg, &
1599 den, denfac, vtw, vtr, ccn, rh_rain, h_var, reevap)
1601 evaporation = evaporation + reevap * convt
1607 call ice_cloud (ks, ke, tz, qv, ql, qr, qi, qs, qg, den, &
1608 denfac, vtw, vtr, vti, vts, vtg, dts, h_var)
1610 if (do_subgrid_proc)
then
1616 call subgrid_z_proc (ks, ke, den, denfac, dts, rh_adj, tz, qv, ql, &
1617 qr, qi, qs, qg, dp, ccn, cin, cond, dep, reevap, sub, last_step)
1619 condensation = condensation + cond * convt
1620 deposition = deposition + dep * convt
1621 evaporation = evaporation + reevap * convt
1622 sublimation = sublimation + sub * convt
1634subroutine mp_fast (ks, ke, tz, qv, ql, qr, qi, qs, qg, dtm, dp, den, &
1635 ccn, cin, condensation, deposition, evaporation, sublimation, &
1636 denfac, convt, last_step)
1644 logical,
intent (in) :: last_step
1646 integer,
intent (in) :: ks, ke
1648 real(kind_phys),
intent (in) :: dtm, convt
1650 real(kind_phys),
intent (in),
dimension (ks:ke) :: dp, den, denfac
1652 real(kind_phys),
intent (inout),
dimension (ks:ke) :: qv, ql, qr, qi, qs, qg, ccn, cin
1654 real (kind = r8),
intent (inout),
dimension (ks:ke) :: tz
1656 real(kind_phys),
intent (inout) :: condensation, deposition
1657 real(kind_phys),
intent (inout) :: evaporation, sublimation
1663 logical :: cond_evap
1667 real(kind_phys) :: cond, dep, reevap, sub
1669 real(kind_phys),
dimension (ks:ke) :: q_liq, q_sol, lcpk, icpk, tcpk, tcp3
1671 real (kind = r8),
dimension (ks:ke) :: cvm, te8
1686 call cal_mhc_lhc (ks, ke, qv, ql, qr, qi, qs, qg, q_liq, q_sol, cvm, te8, tz, &
1687 lcpk, icpk, tcpk, tcp3)
1689 if (.not. do_warm_rain_mp .and. fast_fr_mlt)
then
1695 call pimlt (ks, ke, dtm, qv, ql, qr, qi, qs, qg, tz, cvm, te8, &
1696 lcpk, icpk, tcpk, tcp3)
1702 call pcomp (ks, ke, qv, ql, qr, qi, qs, qg, tz, cvm, te8, &
1703 lcpk, icpk, tcpk, tcp3)
1711 if (delay_cond_evap)
then
1712 cond_evap = last_step
1719 call pcond_pevap (ks, ke, dtm, qv, ql, qr, qi, qs, qg, tz, dp, cvm, te8, den, &
1720 lcpk, icpk, tcpk, tcp3, cond, reevap)
1724 condensation = condensation + cond * convt
1725 evaporation = evaporation + reevap * convt
1727 if (.not. do_warm_rain_mp .and. fast_fr_mlt)
then
1733 call pifr (ks, ke, qv, ql, qr, qi, qs, qg, tz, cvm, te8, den, &
1734 lcpk, icpk, tcpk, tcp3)
1740 call pwbf (ks, ke, dtm, qv, ql, qr, qi, qs, qg, tz, cvm, te8, den, &
1741 lcpk, icpk, tcpk, tcp3)
1747 call pbigg (ks, ke, dtm, qv, ql, qr, qi, qs, qg, tz, cvm, te8, den, ccn, &
1748 lcpk, icpk, tcpk, tcp3)
1754 call pgfr_simp (ks, ke, dtm, qv, ql, qr, qi, qs, qg, tz, cvm, te8, &
1755 lcpk, icpk, tcpk, tcp3)
1761 call psmlt_simp (ks, ke, dtm, qv, ql, qr, qi, qs, qg, tz, cvm, te8, &
1762 lcpk, icpk, tcpk, tcp3)
1770 call praut_simp (ks, ke, dtm, tz, qv, ql, qr, qi, qs, qg)
1772 if (.not. do_warm_rain_mp .and. fast_dep_sub)
then
1778 call pidep_pisub (ks, ke, dtm, qv, ql, qr, qi, qs, qg, tz, dp, cvm, te8, den, &
1779 lcpk, icpk, tcpk, tcp3, cin, dep, sub)
1781 deposition = deposition + dep * convt
1782 sublimation = sublimation + sub * convt
1788 call psaut_simp (ks, ke, dtm, qv, ql, qr, qi, qs, qg, tz, den)
1794 call psdep_pssub (ks, ke, dtm, qv, ql, qr, qi, qs, qg, tz, dp, cvm, te8, den, &
1795 denfac, lcpk, icpk, tcpk, tcp3, dep, sub)
1801 call pgdep_pgsub (ks, ke, dtm, qv, ql, qr, qi, qs, qg, tz, dp, cvm, te8, den, &
1802 denfac, lcpk, icpk, tcpk, tcp3, dep, sub)
1812subroutine sedimentation (dts, ks, ke, tz, qv, ql, qr, qi, qs, qg, dz, dp, &
1813 vtw, vtr, vti, vts, vtg, w1, r1, i1, s1, g1, pfw, pfr, pfi, pfs, pfg, &
1814 u, v, w, den, denfac, dte)
1822 integer,
intent (in) :: ks, ke
1824 real(kind_phys),
intent (in) :: dts
1826 real(kind_phys),
intent (in),
dimension (ks:ke) :: dp, dz, den, denfac
1828 real(kind_phys),
intent (inout),
dimension (ks:ke) :: qv, ql, qr, qi, qs, qg, u, v, w
1830 real(kind_phys),
intent (out) :: w1, r1, i1, s1, g1
1832 real(kind_phys),
intent (out),
dimension (ks:ke) :: vtw, vtr, vti, vts, vtg, pfw, pfr, pfi, pfs, pfg
1834 real (kind = r8),
intent (inout) :: dte
1836 real (kind = r8),
intent (inout),
dimension (ks:ke) :: tz
1844 real(kind_phys),
dimension (ks:ke) :: q_liq, q_sol, lcpk, icpk, tcpk, tcp3
1846 real (kind = r8),
dimension (ks:ke) :: te8, cvm
1870 call cal_mhc_lhc (ks, ke, qv, ql, qr, qi, qs, qg, q_liq, q_sol, cvm, te8, tz, &
1871 lcpk, icpk, tcpk, tcp3)
1877 if (do_psd_ice_fall)
then
1878 call term_rsg (ks, ke, qi, den, denfac, vi_fac, blini, mui, tvai, tvbi, vi_max, const_vi, vti)
1880 call term_ice (ks, ke, tz, qi, den, vi_fac, vi_max, const_vi, vti)
1883 if (do_sedi_melt)
then
1884 call sedi_melt (dts, ks, ke, tz, qv, ql, qr, qi, qs, qg, dz, dp, &
1885 vti, r1, tau_imlt, icpk,
"qi")
1888 call terminal_fall (dts, ks, ke, tz, qv, ql, qr, qi, qs, qg, dz, dp, &
1889 vti, i1, pfi, u, v, w, dte,
"qi")
1891 pfi(ks) = max(0.0, pfi(ks))
1892 do k = ke, ks + 1, -1
1893 pfi(k) = max(0.0, pfi(k) - pfi(k - 1))
1900 call term_rsg (ks, ke, qs, den, denfac, vs_fac, blins, mus, tvas, tvbs, vs_max, const_vs, vts)
1902 if (do_sedi_melt)
then
1903 call sedi_melt (dts, ks, ke, tz, qv, ql, qr, qi, qs, qg, dz, dp, &
1904 vts, r1, tau_smlt, icpk,
"qs")
1907 call terminal_fall (dts, ks, ke, tz, qv, ql, qr, qi, qs, qg, dz, dp, &
1908 vts, s1, pfs, u, v, w, dte,
"qs")
1910 pfs(ks) = max(0.0, pfs(ks))
1911 do k = ke, ks + 1, -1
1912 pfs(k) = max(0.0, pfs(k) - pfs(k - 1))
1920 call term_rsg (ks, ke, qg, den, denfac, vg_fac, blinh, muh, tvah, tvbh, vg_max, const_vg, vtg)
1922 call term_rsg (ks, ke, qg, den, denfac, vg_fac, bling, mug, tvag, tvbg, vg_max, const_vg, vtg)
1925 if (do_sedi_melt)
then
1926 call sedi_melt (dts, ks, ke, tz, qv, ql, qr, qi, qs, qg, dz, dp, &
1927 vtg, r1, tau_gmlt, icpk,
"qg")
1930 call terminal_fall (dts, ks, ke, tz, qv, ql, qr, qi, qs, qg, dz, dp, &
1931 vtg, g1, pfg, u, v, w, dte,
"qg")
1933 pfg(ks) = max(0.0, pfg(ks))
1934 do k = ke, ks + 1, -1
1935 pfg(k) = max(0.0, pfg(k) - pfg(k - 1))
1942 if (do_psd_water_fall)
then
1944 call term_rsg (ks, ke, ql, den, denfac, vw_fac, blinw, muw, tvaw, tvbw, vw_max, const_vw, vtw)
1946 call terminal_fall (dts, ks, ke, tz, qv, ql, qr, qi, qs, qg, dz, dp, &
1947 vtw, w1, pfw, u, v, w, dte,
"ql")
1949 pfw(ks) = max(0.0, pfw(ks))
1950 do k = ke, ks + 1, -1
1951 pfw(k) = max(0.0, pfw(k) - pfw(k - 1))
1960 call term_rsg (ks, ke, qr, den, denfac, vr_fac, blinr, mur, tvar, tvbr, vr_max, const_vr, vtr)
1962 call terminal_fall (dts, ks, ke, tz, qv, ql, qr, qi, qs, qg, dz, dp, &
1963 vtr, r1, pfr, u, v, w, dte,
"qr")
1965 pfr(ks) = max(0.0, pfr(ks))
1966 do k = ke, ks + 1, -1
1967 pfr(k) = max(0.0, pfr(k) - pfr(k - 1))
2768subroutine ice_cloud (ks, ke, tz, qv, ql, qr, qi, qs, qg, den, &
2769 denfac, vtw, vtr, vti, vts, vtg, dts, h_var)
2777 integer,
intent (in) :: ks, ke
2779 real(kind_phys),
intent (in) :: dts, h_var
2781 real(kind_phys),
intent (in),
dimension (ks:ke) :: den, denfac, vtw, vtr, vti, vts, vtg
2783 real(kind_phys),
intent (inout),
dimension (ks:ke) :: qv, ql, qr, qi, qs, qg
2785 real (kind = r8),
intent (inout),
dimension (ks:ke) :: tz
2791 real(kind_phys),
dimension (ks:ke) :: di, q_liq, q_sol, lcpk, icpk, tcpk, tcp3
2793 real (kind = r8),
dimension (ks:ke) :: cvm, te8
2799 call cal_mhc_lhc (ks, ke, qv, ql, qr, qi, qs, qg, q_liq, q_sol, cvm, te8, tz, &
2800 lcpk, icpk, tcpk, tcp3)
2802 if (.not. do_warm_rain_mp)
then
2808 call pimlt (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, cvm, te8, lcpk, icpk, tcpk, tcp3)
2814 call pifr (ks, ke, qv, ql, qr, qi, qs, qg, tz, cvm, te8, den, lcpk, icpk, tcpk, tcp3)
2820 call linear_prof (ke - ks + 1, qi, di, z_slope_ice, h_var)
2826 call psmlt (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, cvm, te8, den, denfac, &
2827 vtw, vtr, vts, lcpk, icpk, tcpk, tcp3)
2833 call pgmlt (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, cvm, te8, den, denfac, &
2834 vtw, vtr, vtg, lcpk, icpk, tcpk, tcp3)
2840 call psaci (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, den, denfac, vti, vts)
2846 call psaut (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, den, di)
2852 call pgaci (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, den, denfac, vti, vtg)
2858 call psacr_pgfr (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, cvm, te8, den, denfac, &
2859 vtr, vts, lcpk, icpk, tcpk, tcp3)
2865 call pgacs (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, den, vts, vtg)
2871 call pgaut (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, den)
2877 call pgacw_pgacr (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, cvm, te8, den, denfac, &
2878 vtr, vtg, lcpk, icpk, tcpk, tcp3)
3590subroutine subgrid_z_proc (ks, ke, den, denfac, dts, rh_adj, tz, qv, ql, qr, &
3591 qi, qs, qg, dp, ccn, cin, cond, dep, reevap, sub, last_step)
3599 logical,
intent (in) :: last_step
3601 integer,
intent (in) :: ks, ke
3603 real(kind_phys),
intent (in) :: dts, rh_adj
3605 real(kind_phys),
intent (in),
dimension (ks:ke) :: den, denfac, dp
3607 real(kind_phys),
intent (inout),
dimension (ks:ke) :: qv, ql, qr, qi, qs, qg, ccn, cin
3609 real(kind_phys),
intent (out) :: cond, dep, reevap, sub
3611 real (kind = r8),
intent (inout),
dimension (ks:ke) :: tz
3617 logical :: cond_evap
3621 real(kind_phys),
dimension (ks:ke) :: q_liq, q_sol, q_cond, lcpk, icpk, tcpk, tcp3
3623 real (kind = r8),
dimension (ks:ke) :: cvm, te8
3638 call cal_mhc_lhc (ks, ke, qv, ql, qr, qi, qs, qg, q_liq, q_sol, cvm, te8, tz, &
3639 lcpk, icpk, tcpk, tcp3)
3645 if (.not. do_warm_rain_mp)
then
3647 call pinst (ks, ke, qv, ql, qr, qi, qs, qg, tz, dp, cvm, te8, den, &
3648 lcpk, icpk, tcpk, tcp3, rh_adj, dep, sub, reevap)
3656 if (delay_cond_evap)
then
3657 cond_evap = last_step
3664 call pcond_pevap (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, dp, cvm, te8, den, &
3665 lcpk, icpk, tcpk, tcp3, cond, reevap)
3669 if (.not. do_warm_rain_mp)
then
3675 call pcomp (ks, ke, qv, ql, qr, qi, qs, qg, tz, cvm, te8, lcpk, icpk, tcpk, tcp3)
3681 call pwbf (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, cvm, te8, den, lcpk, icpk, tcpk, tcp3)
3687 call pbigg (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, cvm, te8, den, ccn, lcpk, icpk, tcpk, tcp3)
3693 call pidep_pisub (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, dp, cvm, te8, den, &
3694 lcpk, icpk, tcpk, tcp3, cin, dep, sub)
3700 call psdep_pssub (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, dp, cvm, te8, den, &
3701 denfac, lcpk, icpk, tcpk, tcp3, dep, sub)
3707 call pgdep_pgsub (ks, ke, dts, qv, ql, qr, qi, qs, qg, tz, dp, cvm, te8, den, &
3708 denfac, lcpk, icpk, tcpk, tcp3, dep, sub)
4603subroutine cs_profile (a4, del, km)
4611 integer,
intent (in) :: km
4613 real(kind_phys),
intent (in) :: del (km)
4615 real(kind_phys),
intent (inout) :: a4 (4, km)
4623 logical :: extm (km)
4625 real(kind_phys) :: gam (km), q (km + 1), d4, bet, a_bot, grat, pmp, lac
4626 real(kind_phys) :: pmp_1, lac_1, pmp_2, lac_2, da1, da2, a6da
4628 grat = del(2) / del(1)
4629 bet = grat * (grat + 0.5)
4630 q(1) = (2. * grat * (grat + 1.) * a4(1, 1) + a4(1, 2)) / bet
4631 gam(1) = (1. + grat * (grat + 1.5)) / bet
4634 d4 = del(k - 1) / del(k)
4635 bet = 2. + 2. * d4 - gam(k - 1)
4636 q(k) = (3. * (a4(1, k - 1) + d4 * a4(1, k)) - q(k - 1)) / bet
4640 a_bot = 1. + d4 * (d4 + 1.5)
4641 q(km + 1) = (2. * d4 * (d4 + 1.) * a4(1, km) + a4(1, km - 1) - a_bot * q(km)) &
4642 / (d4 * (d4 + 0.5) - a_bot * gam(km))
4645 q(k) = q(k) - gam(k) * q(k + 1)
4653 gam(k) = a4(1, k) - a4(1, k - 1)
4660 q(1) = max(q(1), 0.)
4661 q(2) = min(q(2), max(a4(1, 1), a4(1, 2)))
4662 q(2) = max(q(2), min(a4(1, 1), a4(1, 2)), 0.)
4669 if (gam(k - 1) * gam(k + 1) .gt. 0.)
then
4671 q(k) = min(q(k), max(a4(1, k - 1), a4(1, k)))
4672 q(k) = max(q(k), min(a4(1, k - 1), a4(1, k)))
4674 if (gam(k - 1) .gt. 0.)
then
4676 q(k) = max(q(k), min(a4(1, k - 1), a4(1, k)))
4679 q(k) = min(q(k), max(a4(1, k - 1), a4(1, k)))
4681 q(k) = max(q(k), 0.0)
4690 q(km) = min(q(km), max(a4(1, km - 1), a4(1, km)))
4691 q(km) = max(q(km), min(a4(1, km - 1), a4(1, km)), 0.)
4692 q(km + 1) = max(q(km + 1), 0.)
4700 if (k .eq. 1 .or. k .eq. km)
then
4701 extm(k) = (a4(2, k) - a4(1, k)) * (a4(3, k) - a4(1, k)) .gt. 0.
4703 extm(k) = gam(k) * gam(k + 1) .lt. 0.
4717 a4(2, 1) = max(0., a4(2, 1))
4726 if (a4(1, k) .lt. qcmin .or. extm(k - 1) .or. extm(k + 1))
then
4731 a4(4, k) = 6. * a4(1, k) - 3. * (a4(2, k) + a4(3, k))
4732 if (abs(a4(4, k)) .gt. abs(a4(2, k) - a4(3, k)))
then
4734 pmp_1 = a4(1, k) - 2.0 * gam(k + 1)
4735 lac_1 = pmp_1 + 1.5 * gam(k + 2)
4736 a4(2, k) = min(max(a4(2, k), min(a4(1, k), pmp_1, lac_1)), &
4737 max(a4(1, k), pmp_1, lac_1))
4738 pmp_2 = a4(1, k) + 2.0 * gam(k)
4739 lac_2 = pmp_2 - 1.5 * gam(k - 1)
4740 a4(3, k) = min(max(a4(3, k), min(a4(1, k), pmp_2, lac_2)), &
4741 max(a4(1, k), pmp_2, lac_2))
4747 a4(4, k) = 6. * a4(1, k) - 3. * (a4(2, k) + a4(3, k))
4756 da1 = a4(3, k) - a4(2, k)
4758 a6da = a4(4, k) * da1
4759 if (a6da .lt. - da2)
then
4760 a4(4, k) = 3. * (a4(2, k) - a4(1, k))
4761 a4(3, k) = a4(2, k) - a4(4, k)
4762 elseif (a6da .gt. da2)
then
4763 a4(4, k) = 3. * (a4(3, k) - a4(1, k))
4764 a4(2, k) = a4(3, k) - a4(4, k)
4768 call cs_limiters (km - 1, a4)
4774 a4(2, km) = a4(1, km)
4775 a4(3, km) = a4(1, km)
5592subroutine cld_eff_rad (is, ie, ks, ke, lsm, p, delp, t, qv, qw, qi, qr, qs, qg, qa, &
5593 rew, rei, rer, res, reg, snowd, cnvw, cnvi)
5601 integer,
intent (in) :: is, ie, ks, ke
5603 real(kind_phys),
intent (in),
dimension (is:ie) :: lsm, snowd
5605 real(kind_phys),
intent (in),
dimension (is:ie, ks:ke) :: delp, t, p
5606 real(kind_phys),
intent (in),
dimension (is:ie, ks:ke) :: qv, qw, qi, qr, qs, qg, qa
5608 real(kind_phys),
intent (in),
dimension (is:ie, ks:ke),
optional :: cnvw, cnvi
5610 real(kind_phys),
intent (inout),
dimension (is:ie, ks:ke) :: rew, rei, rer, res, reg
5616 integer :: i, k, ind
5618 real(kind_phys),
dimension (is:ie, ks:ke) :: qcw, qci, qcr, qcs, qcg
5619 real(kind_phys),
dimension (is:ie, ks:ke) :: qmw, qmr, qmi, qms, qmg
5621 real(kind_phys) :: dpg, rho, ccnw, mask, cor, tc, bw
5622 real(kind_phys) :: lambdaw, lambdar, lambdai, lambdas, lambdag, rei_fac
5624 real(kind_phys) :: retab (138) = (/ &
5625 0.05000, 0.05000, 0.05000, 0.05000, 0.05000, 0.05000, &
5626 0.05500, 0.06000, 0.07000, 0.08000, 0.09000, 0.10000, &
5627 0.20000, 0.30000, 0.40000, 0.50000, 0.60000, 0.70000, &
5628 0.80000, 0.90000, 1.00000, 1.10000, 1.20000, 1.30000, &
5629 1.40000, 1.50000, 1.60000, 1.80000, 2.00000, 2.20000, &
5630 2.40000, 2.60000, 2.80000, 3.00000, 3.20000, 3.50000, &
5631 3.80000, 4.10000, 4.40000, 4.70000, 5.00000, 5.30000, &
5632 5.60000, 5.92779, 6.26422, 6.61973, 6.99539, 7.39234, &
5633 7.81177, 8.25496, 8.72323, 9.21800, 9.74075, 10.2930, &
5634 10.8765, 11.4929, 12.1440, 12.8317, 13.5581, 14.2319, &
5635 15.0351, 15.8799, 16.7674, 17.6986, 18.6744, 19.6955, &
5636 20.7623, 21.8757, 23.0364, 24.2452, 25.5034, 26.8125, &
5637 27.7895, 28.6450, 29.4167, 30.1088, 30.7306, 31.2943, &
5638 31.8151, 32.3077, 32.7870, 33.2657, 33.7540, 34.2601, &
5639 34.7892, 35.3442, 35.9255, 36.5316, 37.1602, 37.8078, &
5640 38.4720, 39.1508, 39.8442, 40.5552, 41.2912, 42.0635, &
5641 42.8876, 43.7863, 44.7853, 45.9170, 47.2165, 48.7221, &
5642 50.4710, 52.4980, 54.8315, 57.4898, 60.4785, 63.7898, &
5643 65.5604, 71.2885, 75.4113, 79.7368, 84.2351, 88.8833, &
5644 93.6658, 98.5739, 103.603, 108.752, 114.025, 119.424, &
5645 124.954, 130.630, 136.457, 142.446, 148.608, 154.956, &
5646 161.503, 168.262, 175.248, 182.473, 189.952, 197.699, &
5647 205.728, 214.055, 222.694, 231.661, 240.971, 250.639 /)
5659 if (
present (cnvw))
then
5662 if (
present (cnvi))
then
5670 if (liq_ice_combine)
then
5673 qmw(i, k) = qmw(i, k) + qmr(i, k)
5675 qmi(i, k) = qmi(i, k) + qms(i, k) + qmg(i, k)
5686 if (snow_grauple_combine)
then
5689 qms(i, k) = qms(i, k) + qmg(i, k)
5699 qmw(i, k) = max(qmw(i, k), qcmin)
5700 qmi(i, k) = max(qmi(i, k), qcmin)
5701 qmr(i, k) = max(qmr(i, k), qcmin)
5702 qms(i, k) = max(qms(i, k), qcmin)
5703 qmg(i, k) = max(qmg(i, k), qcmin)
5706 mask = min(max(lsm(i), 0.0), 2.0)
5708 dpg = abs(delp(i, k)) / grav
5709 rho = p(i, k) / (rdgas * t(i, k) * (1. + zvir * qv(i, k)))
5713 if (rewflag .eq. 1)
then
5721 ccnw = (1.0 - abs(mask - 1.0)) * &
5722 (10. ** 2.24 * (qa(i, k) * rho * 1.e9) ** 0.257) + &
5724 (10. ** 2.06 * (qa(i, k) * rho * 1.e9) ** 0.48)
5726 ccnw = ccn_o * abs(mask - 1.0) + ccn_l * (1.0 - abs(mask - 1.0))
5729 if (qmw(i, k) .gt. qcmin)
then
5730 qcw(i, k) = dpg * qmw(i, k) * 1.0e3
5731 rew(i, k) = exp(1.0 / 3.0 * log((3.0 * qmw(i, k) * rho) / &
5732 (4.0 * pi * rhow * ccnw))) * 1.0e4
5733 rew(i, k) = max(rewmin, min(rewmax, rew(i, k)))
5741 if (rewflag .eq. 2)
then
5749 ccnw = (1.0 - abs(mask - 1.0)) * &
5750 (10. ** 2.24 * (qa(i, k) * rho * 1.e9) ** 0.257) + &
5752 (10. ** 2.06 * (qa(i, k) * rho * 1.e9) ** 0.48)
5754 ccnw = 1.077 * ccn_o * abs(mask - 1.0) + 1.143 * ccn_l * (1.0 - abs(mask - 1.0))
5757 if (qmw(i, k) .gt. qcmin)
then
5758 qcw(i, k) = dpg * qmw(i, k) * 1.0e3
5759 rew(i, k) = exp(1.0 / 3.0 * log((3.0 * qmw(i, k) * rho) / &
5760 (4.0 * pi * rhow * ccnw))) * 1.0e4
5761 rew(i, k) = max(rewmin, min(rewmax, rew(i, k)))
5769 if (rewflag .eq. 3)
then
5775 if (qmw(i, k) .gt. qcmin)
then
5776 qcw(i, k) = dpg * qmw(i, k) * 1.0e3
5777 rew(i, k) = 14.0 * abs(mask - 1.0) + &
5778 (8.0 + (14.0 - 8.0) * min(1.0, max(0.0, - tc / 30.0))) * &
5779 (1.0 - abs(mask - 1.0))
5780 rew(i, k) = rew(i, k) + (14.0 - rew(i, k)) * &
5781 min(1.0, max(0.0, snowd(i) / 1000.0))
5782 rew(i, k) = max(rewmin, min(rewmax, rew(i, k)))
5790 if (rewflag .eq. 4)
then
5796 if (qmw(i, k) .gt. qcmin)
then
5797 qcw(i, k) = dpg * qmw(i, k) * 1.0e3
5798 call cal_pc_ed_oe_rr_tv (qmw(i, k), rho, blinw, muw, &
5799 eda = edaw, edb = edbw, ed = rew(i, k))
5800 rew(i, k) = rewfac * 0.5 * rew(i, k) * 1.0e6
5801 rew(i, k) = max(rewmin, min(rewmax, rew(i, k)))
5809 if (reiflag .eq. 1)
then
5815 if (qmi(i, k) .gt. qcmin)
then
5816 qci(i, k) = dpg * qmi(i, k) * 1.0e3
5817 rei_fac = log(1.0e3 * qmi(i, k) * rho)
5818 if (tc .lt. - 50)
then
5819 rei(i, k) = beta / 9.917 * exp(0.109 * rei_fac) * 1.0e3
5820 elseif (tc .lt. - 40)
then
5821 rei(i, k) = beta / 9.337 * exp(0.080 * rei_fac) * 1.0e3
5822 elseif (tc .lt. - 30)
then
5823 rei(i, k) = beta / 9.208 * exp(0.055 * rei_fac) * 1.0e3
5825 rei(i, k) = beta / 9.387 * exp(0.031 * rei_fac) * 1.0e3
5827 rei(i, k) = max(reimin, min(reimax, rei(i, k)))
5835 if (reiflag .eq. 2)
then
5841 if (qmi(i, k) .gt. qcmin)
then
5842 qci(i, k) = dpg * qmi(i, k) * 1.0e3
5843 if (tc .le. - 55)
then
5844 rei(i, k) = 15.41627
5845 elseif (tc .le. - 50)
then
5846 rei(i, k) = 16.60895
5847 elseif (tc .le. - 45)
then
5848 rei(i, k) = 32.89967
5849 elseif (tc .le. - 40)
then
5850 rei(i, k) = 35.29989
5851 elseif (tc .le. - 35)
then
5852 rei(i, k) = 55.65818
5853 elseif (tc .le. - 30)
then
5854 rei(i, k) = 85.19071
5855 elseif (tc .le. - 25)
then
5856 rei(i, k) = 72.35392
5858 rei(i, k) = 92.46298
5860 rei(i, k) = max(reimin, min(reimax, rei(i, k)))
5868 if (reiflag .eq. 3)
then
5874 if (qmi(i, k) .gt. qcmin)
then
5875 qci(i, k) = dpg * qmi(i, k) * 1.0e3
5876 rei(i, k) = 47.05 + tc * (0.6624 + 0.001741 * tc)
5877 rei(i, k) = max(reimin, min(reimax, rei(i, k)))
5885 if (reiflag .eq. 4)
then
5891 if (qmi(i, k) .gt. qcmin)
then
5892 qci(i, k) = dpg * qmi(i, k) * 1.0e3
5893 ind = min(max(int(t(i, k) - 136.0), 44), 138 - 1)
5894 cor = t(i, k) - int(t(i, k))
5895 rei(i, k) = retab(ind) * (1. - cor) + retab(ind + 1) * cor
5896 rei(i, k) = max(reimin, min(reimax, rei(i, k)))
5904 if (reiflag .eq. 5)
then
5910 if (qmi(i, k) .gt. qcmin)
then
5911 qci(i, k) = dpg * qmi(i, k) * 1.0e3
5912 bw = - 2. + 1.e-3 * log10(rho * qmi(i, k) / 50.e-3) * &
5913 exp(1.5 * log(max(1.e-10, - tc)))
5914 rei(i, k) = 377.4 + bw * (203.3 + bw * (37.91 + 2.3696 * bw))
5915 rei(i, k) = max(reimin, min(reimax, rei(i, k)))
5923 if (reiflag .eq. 6)
then
5929 if (qmi(i, k) .gt. qcmin)
then
5930 qci(i, k) = dpg * qmi(i, k) * 1.0e3
5931 rei_fac = log(1.0e3 * qmi(i, k) * rho)
5932 rei(i, k) = 45.8966 * exp(0.2214 * rei_fac) + &
5933 0.7957 * exp(0.2535 * rei_fac) * (tc + 190.0)
5934 rei(i, k) = (1.2351 + 0.0105 * tc) * rei(i, k)
5935 rei(i, k) = max(reimin, min(reimax, rei(i, k)))
5943 if (reiflag .eq. 7)
then
5949 if (qmi(i, k) .gt. qcmin)
then
5950 qci(i, k) = dpg * qmi(i, k) * 1.0e3
5951 call cal_pc_ed_oe_rr_tv (qmi(i, k), rho, blini, mui, &
5952 eda = edai, edb = edbi, ed = rei(i, k))
5953 rei(i, k) = reifac * 0.5 * rei(i, k) * 1.0e6
5954 rei(i, k) = max(reimin, min(reimax, rei(i, k)))
5962 if (rerflag .eq. 1)
then
5968 if (qmr(i, k) .gt. qcmin)
then
5969 qcr(i, k) = dpg * qmr(i, k) * 1.0e3
5970 call cal_pc_ed_oe_rr_tv (qmr(i, k), rho, blinr, mur, &
5971 eda = edar, edb = edbr, ed = rer(i, k))
5972 rer(i, k) = 0.5 * rer(i, k) * 1.0e6
5973 rer(i, k) = max(rermin, min(rermax, rer(i, k)))
5981 if (resflag .eq. 1)
then
5987 if (qms(i, k) .gt. qcmin)
then
5988 qcs(i, k) = dpg * qms(i, k) * 1.0e3
5989 call cal_pc_ed_oe_rr_tv (qms(i, k), rho, blins, mus, &
5990 eda = edas, edb = edbs, ed = res(i, k))
5991 res(i, k) = 0.5 * res(i, k) * 1.0e6
5992 res(i, k) = max(resmin, min(resmax, res(i, k)))
6000 if (regflag .eq. 1)
then
6006 if (qmg(i, k) .gt. qcmin)
then
6007 qcg(i, k) = dpg * qmg(i, k) * 1.0e3
6009 call cal_pc_ed_oe_rr_tv (qmg(i, k), rho, blinh, muh, &
6010 eda = edah, edb = edbh, ed = reg(i, k))
6012 call cal_pc_ed_oe_rr_tv (qmg(i, k), rho, bling, mug, &
6013 eda = edag, edb = edbg, ed = reg(i, k))
6015 reg(i, k) = 0.5 * reg(i, k) * 1.0e6
6016 reg(i, k) = max(regmin, min(regmax, reg(i, k)))
6034subroutine rad_ref (is, ie, js, je, qv, qr, qs, qg, pt, delp, &
6035 delz, dbz, npz, hydrostatic, do_inline_mp, mp_top)
6043 logical,
intent (in) :: hydrostatic, do_inline_mp
6045 integer,
intent (in) :: is, ie, js, je
6046 integer,
intent (in) :: npz, mp_top
6051 real(kind_phys),
intent (in),
dimension (is:ie, js:je, npz) :: delz
6053 real(kind_phys),
intent (in),
dimension (is:ie, js:je, npz) :: pt, delp
6055 real(kind_phys),
intent (in),
dimension (is:ie, js:je, npz) :: qv, qr, qs, qg
6063 real(kind_phys),
intent (out),
dimension (is:ie, js:je, npz) :: dbz
6071 real(kind_phys),
parameter :: alpha = 0.224, mp_const = 200 * exp(1.6 * log(3.6e6))
6073 real (kind = r8) :: qden, z_e
6074 real(kind_phys) :: fac_r, fac_s, fac_g
6075 real(kind_phys) :: allmax
6076 real(kind_phys),
dimension (is:ie, js:je) :: maxdbz
6078 real(kind_phys),
dimension (npz) :: den, denfac, qmr, qms, qmg, vtr, vts, vtg
6113 den(k) = - delp(i, j, k) / (grav * delz(i, j, k))
6114 qmr(k) = max(qcmin, qr(i, j, k))
6115 qms(k) = max(qcmin, qs(i, j, k))
6116 qmg(k) = max(qcmin, qg(i, j, k))
6120 denfac(k) = sqrt(den(npz) / den(k))
6127 if (radr_flag .eq. 3)
then
6128 call term_rsg (1, npz, qmr, den, denfac, vr_fac, blinr, &
6129 mur, tvar, tvbr, vr_max, const_vr, vtr)
6133 if (rads_flag .eq. 3)
then
6134 call term_rsg (1, npz, qms, den, denfac, vs_fac, blins, &
6135 mus, tvas, tvbs, vs_max, const_vs, vts)
6139 if (radg_flag .eq. 3)
then
6140 if (do_hail .and. .not. do_inline_mp)
then
6141 call term_rsg (1, npz, qmg, den, denfac, vg_fac, blinh, &
6142 muh, tvah, tvbh, vg_max, const_vg, vtg)
6145 call term_rsg (1, npz, qmg, den, denfac, vg_fac, bling, &
6146 mug, tvag, tvbg, vg_max, const_vg, vtg)
6155 do k = mp_top + 1, npz
6159 qden = den(k) * qmr(k)
6160 if (qmr(k) .gt. qcmin)
then
6161 call cal_pc_ed_oe_rr_tv (qmr(k), den(k), blinr, mur, &
6162 rra = rrar, rrb = rrbr, rr = fac_r)
6166 if (radr_flag .eq. 1 .or. radr_flag .eq. 2)
then
6167 z_e = z_e + fac_r * 1.e18
6169 if (radr_flag .eq. 3)
then
6170 z_e = z_e + mp_const * exp(1.6 * log(qden * vtr(k)))
6175 qden = den(k) * qms(k)
6176 if (qms(k) .gt. qcmin)
then
6177 call cal_pc_ed_oe_rr_tv (qms(k), den(k), blins, mus, &
6178 rra = rras, rrb = rrbs, rr = fac_s)
6182 if (rads_flag .eq. 1)
then
6183 if (pt(i, j, k) .lt. tice)
then
6184 z_e = z_e + fac_s * 1.e18 * alpha * (rhos / rhor) ** 2
6186 z_e = z_e + fac_s * 1.e18 * alpha * (rhos / rhor) ** 2 / alpha
6189 if (rads_flag .eq. 2)
then
6190 if (pt(i, j, k) .lt. tice)
then
6191 z_e = z_e + fac_s * 1.e18 * alpha * (rhos / rhoi) ** 2
6193 z_e = z_e + fac_s * 1.e18
6196 if (rads_flag .eq. 3)
then
6197 z_e = z_e + mp_const * exp(1.6 * log(qden * vts(k)))
6202 qden = den(k) * qmg(k)
6203 if (do_hail .and. .not. do_inline_mp)
then
6204 if (qmg(k) .gt. qcmin)
then
6205 call cal_pc_ed_oe_rr_tv (qmg(k), den(k), blinh, muh, &
6206 rra = rrah, rrb = rrbh, rr = fac_g)
6210 if (radg_flag .eq. 1)
then
6211 if (pt(i, j, k) .lt. tice)
then
6212 z_e = z_e + fac_g * 1.e18 * alpha * (rhoh / rhor) ** 2
6214 z_e = z_e + fac_g * 1.e18 * alpha * (rhoh / rhor) ** 2 / alpha
6217 if (radg_flag .eq. 2)
then
6218 z_e = z_e + fac_g * 1.e18
6221 if (qmg(k) .gt. qcmin)
then
6222 call cal_pc_ed_oe_rr_tv (qmg(k), den(k), bling, mug, &
6223 rra = rrag, rrb = rrbg, rr = fac_g)
6227 if (radg_flag .eq. 1)
then
6228 if (pt(i, j, k) .lt. tice)
then
6229 z_e = z_e + fac_g * 1.e18 * alpha * (rhog / rhor) ** 2
6231 z_e = z_e + fac_g * 1.e18 * alpha * (rhog / rhor) ** 2 / alpha
6234 if (radg_flag .eq. 2)
then
6235 z_e = z_e + fac_g * 1.e18
6238 if (radg_flag .eq. 3)
then
6239 z_e = z_e + mp_const * exp(1.6 * log(qden * vtg(k)))
6243 dbz(i, j, k) = 10. * log10(max(0.01, z_e))
6246 do k = mp_top + 1, npz
6247 maxdbz(i, j) = max(dbz(i, j, k), maxdbz(i, j))
6250 allmax = max(maxdbz(i, j), allmax)