<p><b>laura@ucar.edu</b> 2010-05-20 13:46:32 -0600 (Thu, 20 May 2010)</p><p>Thompson cloud microphysics scheme. As in WRF version 3.2, except for commenting out calls to wrf_debug. Added capabilities to read pre-made look-up tables to speed up initialization.<br>
</p><hr noshade><pre><font color="gray">Added: branches/atmos_physics/src/core_hyd_phys/module_mp_thompson.F
===================================================================
--- branches/atmos_physics/src/core_hyd_phys/module_mp_thompson.F (rev 0)
+++ branches/atmos_physics/src/core_hyd_phys/module_mp_thompson.F 2010-05-20 19:46:32 UTC (rev 284)
@@ -0,0 +1,3653 @@
+!+---+-----------------------------------------------------------------+
+!.. This subroutine computes the moisture tendencies of water vapor,
+!.. cloud droplets, rain, cloud ice (pristine), snow, and graupel.
+!.. Prior to WRFv2.2 this code was based on Reisner et al (1998), but
+!.. few of those pieces remain. A complete description is now found in
+!.. Thompson, G., P. R. Field, R. M. Rasmussen, and W. D. Hall, 2008:
+!.. Explicit Forecasts of winter precipitation using an improved bulk
+!.. microphysics scheme. Part II: Implementation of a new snow
+!.. parameterization. Mon. Wea. Rev., 136, 5095-5115.
+!.. Prior to WRFv3.1, this code was single-moment rain prediction as
+!.. described in the reference above, but in v3.1 and higher, the
+!.. scheme is two-moment rain (predicted rain number concentration).
+!..
+!.. Most importantly, users may wish to modify the prescribed number of
+!.. cloud droplets (Nt_c; see guidelines mentioned below). Otherwise,
+!.. users may alter the rain and graupel size distribution parameters
+!.. to use exponential (Marshal-Palmer) or generalized gamma shape.
+!.. The snow field assumes a combination of two gamma functions (from
+!.. Field et al. 2005) and would require significant modifications
+!.. throughout the entire code to alter its shape as well as accretion
+!.. rates. Users may also alter the constants used for density of rain,
+!.. graupel, ice, and snow, but the latter is not constant when using
+!.. Paul Field's snow distribution and moments methods. Other values
+!.. users can modify include the constants for mass and/or velocity
+!.. power law relations and assumed capacitances used in deposition/
+!.. sublimation/evaporation/melting.
+!.. Remaining values should probably be left alone.
+!..
+!..Author: Greg Thompson, NCAR-RAL, gthompsn@ucar.edu, 303-497-2805
+!..Last modified: 09 Nov 2009
+!+---+-----------------------------------------------------------------+
+!wrft:model_layer:physics
+!+---+-----------------------------------------------------------------+
+!
+ MODULE module_mp_thompson
+! USE module_wrf_error
+! USE module_utility, ONLY: WRFU_Clock, WRFU_Alarm
+! USE module_domain, ONLY : HISTORY_ALARM, Is_alarm_tstep
+
+ IMPLICIT NONE
+
+!LDF begin (05-13-2010): Added the capabilities to read pre-calculated
+!look-up tables to speed up initialization.
+ LOGICAL, PRIVATE:: iiwarm
+ LOGICAL, PRIVATE:: l_qr_acr_qg
+ LOGICAL, PRIVATE:: l_qr_acr_qs
+ LOGICAL, PRIVATE:: l_qi_aut_qs
+ LOGICAL, PRIVATE:: l_freezeH2O
+!LDF end.
+! LOGICAL, PARAMETER, PRIVATE:: iiwarm = .false.
+ INTEGER, PARAMETER, PRIVATE:: IFDRY = 0
+ REAL, PARAMETER, PRIVATE:: T_0 = 273.15
+ REAL, PARAMETER, PRIVATE:: PI = 3.1415926536
+
+!..Densities of rain, snow, graupel, and cloud ice.
+ REAL, PARAMETER, PRIVATE:: rho_w = 1000.0
+ REAL, PARAMETER, PRIVATE:: rho_s = 100.0
+ REAL, PARAMETER, PRIVATE:: rho_g = 400.0
+ REAL, PARAMETER, PRIVATE:: rho_i = 890.0
+
+!..Prescribed number of cloud droplets. Set according to known data or
+!.. roughly 100 per cc (100.E6 m^-3) for Maritime cases and
+!.. 300 per cc (300.E6 m^-3) for Continental. Gamma shape parameter,
+!.. mu_c, calculated based on Nt_c is important in autoconversion
+!.. scheme.
+ REAL, PARAMETER, PRIVATE:: Nt_c = 100.E6
+
+!..Generalized gamma distributions for rain, graupel and cloud ice.
+!.. N(D) = N_0 * D**mu * exp(-lamda*D); mu=0 is exponential.
+ REAL, PARAMETER, PRIVATE:: mu_r = 0.0
+ REAL, PARAMETER, PRIVATE:: mu_g = 0.0
+ REAL, PARAMETER, PRIVATE:: mu_i = 0.0
+ REAL, PRIVATE:: mu_c
+
+!..Sum of two gamma distrib for snow (Field et al. 2005).
+!.. N(D) = M2**4/M3**3 * [Kap0*exp(-M2*Lam0*D/M3)
+!.. + Kap1*(M2/M3)**mu_s * D**mu_s * exp(-M2*Lam1*D/M3)]
+!.. M2 and M3 are the (bm_s)th and (bm_s+1)th moments respectively
+!.. calculated as function of ice water content and temperature.
+ REAL, PARAMETER, PRIVATE:: mu_s = 0.6357
+ REAL, PARAMETER, PRIVATE:: Kap0 = 490.6
+ REAL, PARAMETER, PRIVATE:: Kap1 = 17.46
+ REAL, PARAMETER, PRIVATE:: Lam0 = 20.78
+ REAL, PARAMETER, PRIVATE:: Lam1 = 3.29
+
+!..Y-intercept parameter for graupel is not constant and depends on
+!.. mixing ratio. Also, when mu_g is non-zero, these become equiv
+!.. y-intercept for an exponential distrib and proper values are
+!.. computed based on same mixing ratio and total number concentration.
+ REAL, PARAMETER, PRIVATE:: gonv_min = 1.E4
+ REAL, PARAMETER, PRIVATE:: gonv_max = 3.E6
+
+!..Mass power law relations: mass = am*D**bm
+!.. Snow from Field et al. (2005), others assume spherical form.
+ REAL, PARAMETER, PRIVATE:: am_r = PI*rho_w/6.0
+ REAL, PARAMETER, PRIVATE:: bm_r = 3.0
+ REAL, PARAMETER, PRIVATE:: am_s = 0.069
+ REAL, PARAMETER, PRIVATE:: bm_s = 2.0
+ REAL, PARAMETER, PRIVATE:: am_g = PI*rho_g/6.0
+ REAL, PARAMETER, PRIVATE:: bm_g = 3.0
+ REAL, PARAMETER, PRIVATE:: am_i = PI*rho_i/6.0
+ REAL, PARAMETER, PRIVATE:: bm_i = 3.0
+
+!..Fallspeed power laws relations: v = (av*D**bv)*exp(-fv*D)
+!.. Rain from Ferrier (1994), ice, snow, and graupel from
+!.. Thompson et al (2008). Coefficient fv is zero for graupel/ice.
+ REAL, PARAMETER, PRIVATE:: av_r = 4854.0
+ REAL, PARAMETER, PRIVATE:: bv_r = 1.0
+ REAL, PARAMETER, PRIVATE:: fv_r = 195.0
+ REAL, PARAMETER, PRIVATE:: av_s = 40.0
+ REAL, PARAMETER, PRIVATE:: bv_s = 0.55
+ REAL, PARAMETER, PRIVATE:: fv_s = 125.0
+ REAL, PARAMETER, PRIVATE:: av_g = 442.0
+ REAL, PARAMETER, PRIVATE:: bv_g = 0.89
+ REAL, PARAMETER, PRIVATE:: av_i = 1847.5
+ REAL, PARAMETER, PRIVATE:: bv_i = 1.0
+
+!..Capacitance of sphere and plates/aggregates: D**3, D**2
+ REAL, PARAMETER, PRIVATE:: C_cube = 0.5
+ REAL, PARAMETER, PRIVATE:: C_sqrd = 0.3
+
+!..Collection efficiencies. Rain/snow/graupel collection of cloud
+!.. droplets use variables (Ef_rw, Ef_sw, Ef_gw respectively) and
+!.. get computed elsewhere because they are dependent on stokes
+!.. number.
+ REAL, PARAMETER, PRIVATE:: Ef_si = 0.05
+ REAL, PARAMETER, PRIVATE:: Ef_rs = 0.95
+ REAL, PARAMETER, PRIVATE:: Ef_rg = 0.75
+ REAL, PARAMETER, PRIVATE:: Ef_ri = 0.95
+
+!..Minimum microphys values
+!.. R1 value, 1.E-12, cannot be set lower because of numerical
+!.. problems with Paul Field's moments and should not be set larger
+!.. because of truncation problems in snow/ice growth.
+ REAL, PARAMETER, PRIVATE:: R1 = 1.E-12
+ REAL, PARAMETER, PRIVATE:: R2 = 1.E-8
+ REAL, PARAMETER, PRIVATE:: eps = 1.E-29
+
+!..Constants in Cooper curve relation for cloud ice number.
+ REAL, PARAMETER, PRIVATE:: TNO = 5.0
+ REAL, PARAMETER, PRIVATE:: ATO = 0.304
+
+!..Rho_not used in fallspeed relations (rho_not/rho)**.5 adjustment.
+ REAL, PARAMETER, PRIVATE:: rho_not = 101325.0/(287.05*298.0)
+
+!..Schmidt number
+ REAL, PARAMETER, PRIVATE:: Sc = 0.632
+ REAL, PRIVATE:: Sc3
+
+!..Homogeneous freezing temperature
+ REAL, PARAMETER, PRIVATE:: HGFR = 235.16
+
+!..Water vapor and air gas constants at constant pressure
+ REAL, PARAMETER, PRIVATE:: Rv = 461.5
+ REAL, PARAMETER, PRIVATE:: oRv = 1./Rv
+ REAL, PARAMETER, PRIVATE:: R = 287.04
+ REAL, PARAMETER, PRIVATE:: Cp = 1004.0
+
+!..Enthalpy of sublimation, vaporization, and fusion at 0C.
+ REAL, PARAMETER, PRIVATE:: lsub = 2.834E6
+ REAL, PARAMETER, PRIVATE:: lvap0 = 2.5E6
+ REAL, PARAMETER, PRIVATE:: lfus = lsub - lvap0
+ REAL, PARAMETER, PRIVATE:: olfus = 1./lfus
+
+!..Ice initiates with this mass (kg), corresponding diameter calc.
+!..Min diameters and mass of cloud, rain, snow, and graupel (m, kg).
+ REAL, PARAMETER, PRIVATE:: xm0i = 1.E-12
+ REAL, PARAMETER, PRIVATE:: D0c = 1.E-6
+ REAL, PARAMETER, PRIVATE:: D0r = 50.E-6
+ REAL, PARAMETER, PRIVATE:: D0s = 200.E-6
+ REAL, PARAMETER, PRIVATE:: D0g = 250.E-6
+ REAL, PRIVATE:: D0i, xm0s, xm0g
+
+!..Lookup table dimensions
+ INTEGER, PARAMETER, PRIVATE:: nbins = 100
+ INTEGER, PARAMETER, PRIVATE:: nbc = nbins
+ INTEGER, PARAMETER, PRIVATE:: nbi = nbins
+ INTEGER, PARAMETER, PRIVATE:: nbr = nbins
+ INTEGER, PARAMETER, PRIVATE:: nbs = nbins
+ INTEGER, PARAMETER, PRIVATE:: nbg = nbins
+ INTEGER, PARAMETER, PRIVATE:: ntb_c = 37
+ INTEGER, PARAMETER, PRIVATE:: ntb_i = 64
+ INTEGER, PARAMETER, PRIVATE:: ntb_r = 37
+ INTEGER, PARAMETER, PRIVATE:: ntb_s = 28
+ INTEGER, PARAMETER, PRIVATE:: ntb_g = 28
+ INTEGER, PARAMETER, PRIVATE:: ntb_g1 = 28
+ INTEGER, PARAMETER, PRIVATE:: ntb_r1 = 37
+ INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55
+ INTEGER, PARAMETER, PRIVATE:: ntb_t = 9
+ INTEGER, PRIVATE:: nic2, nii2, nii3, nir2, nir3, nis2, nig2, nig3
+
+ DOUBLE PRECISION, DIMENSION(nbins+1):: xDx
+ DOUBLE PRECISION, DIMENSION(nbc):: Dc, dtc
+ DOUBLE PRECISION, DIMENSION(nbi):: Di, dti
+ DOUBLE PRECISION, DIMENSION(nbr):: Dr, dtr
+ DOUBLE PRECISION, DIMENSION(nbs):: Ds, dts
+ DOUBLE PRECISION, DIMENSION(nbg):: Dg, dtg
+
+!..Lookup tables for cloud water content (kg/m**3).
+ REAL, DIMENSION(ntb_c), PARAMETER, PRIVATE:: &
+ r_c = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, &
+ 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
+ 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, &
+ 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, &
+ 1.e-2/)
+
+!..Lookup tables for cloud ice content (kg/m**3).
+ REAL, DIMENSION(ntb_i), PARAMETER, PRIVATE:: &
+ r_i = (/1.e-10,2.e-10,3.e-10,4.e-10, &
+ 5.e-10,6.e-10,7.e-10,8.e-10,9.e-10, &
+ 1.e-9,2.e-9,3.e-9,4.e-9,5.e-9,6.e-9,7.e-9,8.e-9,9.e-9, &
+ 1.e-8,2.e-8,3.e-8,4.e-8,5.e-8,6.e-8,7.e-8,8.e-8,9.e-8, &
+ 1.e-7,2.e-7,3.e-7,4.e-7,5.e-7,6.e-7,7.e-7,8.e-7,9.e-7, &
+ 1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, &
+ 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
+ 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, &
+ 1.e-3/)
+
+!..Lookup tables for rain content (kg/m**3).
+ REAL, DIMENSION(ntb_r), PARAMETER, PRIVATE:: &
+ r_r = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, &
+ 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
+ 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, &
+ 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, &
+ 1.e-2/)
+
+!..Lookup tables for graupel content (kg/m**3).
+ REAL, DIMENSION(ntb_g), PARAMETER, PRIVATE:: &
+ r_g = (/1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
+ 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, &
+ 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, &
+ 1.e-2/)
+
+!..Lookup tables for snow content (kg/m**3).
+ REAL, DIMENSION(ntb_s), PARAMETER, PRIVATE:: &
+ r_s = (/1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
+ 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, &
+ 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, &
+ 1.e-2/)
+
+!..Lookup tables for rain y-intercept parameter (/m**4).
+ REAL, DIMENSION(ntb_r1), PARAMETER, PRIVATE:: &
+ N0r_exp = (/1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, &
+ 1.e7,2.e7,3.e7,4.e7,5.e7,6.e7,7.e7,8.e7,9.e7, &
+ 1.e8,2.e8,3.e8,4.e8,5.e8,6.e8,7.e8,8.e8,9.e8, &
+ 1.e9,2.e9,3.e9,4.e9,5.e9,6.e9,7.e9,8.e9,9.e9, &
+ 1.e10/)
+
+!..Lookup tables for graupel y-intercept parameter (/m**4).
+ REAL, DIMENSION(ntb_g1), PARAMETER, PRIVATE:: &
+ N0g_exp = (/1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
+ 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
+ 1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, &
+ 1.e7/)
+
+!..Lookup tables for ice number concentration (/m**3).
+ REAL, DIMENSION(ntb_i1), PARAMETER, PRIVATE:: &
+ Nt_i = (/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0, &
+ 1.e1,2.e1,3.e1,4.e1,5.e1,6.e1,7.e1,8.e1,9.e1, &
+ 1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, &
+ 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, &
+ 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
+ 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
+ 1.e6/)
+
+!..For snow moments conversions (from Field et al. 2005)
+ REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
+ sa = (/ 5.065339, -0.062659, -3.032362, 0.029469, -0.000285, &
+ 0.31255, 0.000204, 0.003199, 0.0, -0.015952/)
+ REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
+ sb = (/ 0.476221, -0.015896, 0.165977, 0.007468, -0.000141, &
+ 0.060366, 0.000079, 0.000594, 0.0, -0.003577/)
+
+!..Temperatures (5 C interval 0 to -40) used in lookup tables.
+ REAL, DIMENSION(ntb_t), PARAMETER, PRIVATE:: &
+ Tc = (/-0.01, -5., -10., -15., -20., -25., -30., -35., -40./)
+
+!..Lookup tables for various accretion/collection terms.
+!.. ntb_x refers to the number of elements for rain, snow, graupel,
+!.. and temperature array indices. Variables beginning with t-p/c/m/n
+!.. represent lookup tables. Save compile-time memory by making
+!.. allocatable (2009Jun12, J. Michalakes).
+ INTEGER, PARAMETER, PRIVATE:: R8SIZE = 8
+ REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: &
+ tcg_racg, tmr_racg, tcr_gacr, tmg_gacr, &
+ tnr_racg, tnr_gacr
+ REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: &
+ tcs_racs1, tmr_racs1, tcs_racs2, tmr_racs2, &
+ tcr_sacr1, tms_sacr1, tcr_sacr2, tms_sacr2, &
+ tnr_racs1, tnr_racs2, tnr_sacr1, tnr_sacr2
+ REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: &
+ tpi_qcfz, tni_qcfz
+ REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: &
+ tpi_qrfz, tpg_qrfz, tni_qrfz, tnr_qrfz
+ REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: &
+ tps_iaus, tni_iaus, tpi_ide
+ REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efrw
+ REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efsw
+ REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: tnr_rev
+
+!..Variables holding a bunch of exponents and gamma values (cloud water,
+!.. cloud ice, rain, snow, then graupel).
+ REAL, DIMENSION(3), PRIVATE:: cce, ccg
+ REAL, PRIVATE:: ocg1, ocg2
+ REAL, DIMENSION(6), PRIVATE:: cie, cig
+ REAL, PRIVATE:: oig1, oig2, obmi
+ REAL, DIMENSION(13), PRIVATE:: cre, crg
+ REAL, PRIVATE:: ore1, org1, org2, org3, obmr
+ REAL, DIMENSION(18), PRIVATE:: cse, csg
+ REAL, PRIVATE:: oams, obms, ocms
+ REAL, DIMENSION(12), PRIVATE:: cge, cgg
+ REAL, PRIVATE:: oge1, ogg1, ogg2, ogg3, oamg, obmg, ocmg
+
+!..Declaration of precomputed constants in various rate eqns.
+ REAL:: t1_qr_qc, t1_qr_qi, t2_qr_qi, t1_qg_qc, t1_qs_qc, t1_qs_qi
+ REAL:: t1_qr_ev, t2_qr_ev
+ REAL:: t1_qs_sd, t2_qs_sd, t1_qg_sd, t2_qg_sd
+ REAL:: t1_qs_me, t2_qs_me, t1_qg_me, t2_qg_me
+
+ CHARACTER*256:: mp_debug
+
+!+---+
+!+---+-----------------------------------------------------------------+
+!..END DECLARATIONS
+!+---+-----------------------------------------------------------------+
+!+---+
+!ctrlL
+
+ CONTAINS
+
+ SUBROUTINE thompson_init
+
+ IMPLICIT NONE
+
+ INTEGER:: i, j, k, m, n
+ LOGICAL:: micro_init
+
+!..Allocate space for lookup tables (J. Michalakes 2009Jun08).
+ micro_init = .FALSE.
+
+ if (.NOT. ALLOCATED(tcg_racg) ) then
+ ALLOCATE(tcg_racg(ntb_g1,ntb_g,ntb_r1,ntb_r))
+ micro_init = .TRUE.
+ endif
+
+ if (.NOT. ALLOCATED(tmr_racg)) ALLOCATE(tmr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tcr_gacr)) ALLOCATE(tcr_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tmg_gacr)) ALLOCATE(tmg_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tnr_racg)) ALLOCATE(tnr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tnr_gacr)) ALLOCATE(tnr_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r))
+
+ if (.NOT. ALLOCATED(tcs_racs1)) ALLOCATE(tcs_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tmr_racs1)) ALLOCATE(tmr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tcs_racs2)) ALLOCATE(tcs_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tmr_racs2)) ALLOCATE(tmr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tcr_sacr1)) ALLOCATE(tcr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tms_sacr1)) ALLOCATE(tms_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tcr_sacr2)) ALLOCATE(tcr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tms_sacr2)) ALLOCATE(tms_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tnr_racs1)) ALLOCATE(tnr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tnr_racs2)) ALLOCATE(tnr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tnr_sacr1)) ALLOCATE(tnr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
+ if (.NOT. ALLOCATED(tnr_sacr2)) ALLOCATE(tnr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
+
+ if (.NOT. ALLOCATED(tpi_qcfz)) ALLOCATE(tpi_qcfz(ntb_c,45))
+ if (.NOT. ALLOCATED(tni_qcfz)) ALLOCATE(tni_qcfz(ntb_c,45))
+
+ if (.NOT. ALLOCATED(tpi_qrfz)) ALLOCATE(tpi_qrfz(ntb_r,ntb_r1,45))
+ if (.NOT. ALLOCATED(tpg_qrfz)) ALLOCATE(tpg_qrfz(ntb_r,ntb_r1,45))
+ if (.NOT. ALLOCATED(tni_qrfz)) ALLOCATE(tni_qrfz(ntb_r,ntb_r1,45))
+ if (.NOT. ALLOCATED(tnr_qrfz)) ALLOCATE(tnr_qrfz(ntb_r,ntb_r1,45))
+
+ if (.NOT. ALLOCATED(tps_iaus)) ALLOCATE(tps_iaus(ntb_i,ntb_i1))
+ if (.NOT. ALLOCATED(tni_iaus)) ALLOCATE(tni_iaus(ntb_i,ntb_i1))
+ if (.NOT. ALLOCATED(tpi_ide)) ALLOCATE(tpi_ide(ntb_i,ntb_i1))
+
+ if (.NOT. ALLOCATED(t_Efrw)) ALLOCATE(t_Efrw(nbr,nbc))
+ if (.NOT. ALLOCATED(t_Efsw)) ALLOCATE(t_Efsw(nbs,nbc))
+
+ if (.NOT. ALLOCATED(tnr_rev)) ALLOCATE(tnr_rev(nbr, ntb_r1, ntb_r))
+
+ if (micro_init) then
+
+!..From Martin et al. (1994), assign gamma shape parameter mu for cloud
+!.. drops according to general dispersion characteristics (disp=~0.25
+!.. for Maritime and 0.45 for Continental).
+!.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime
+!.. to 2 for really dirty air.
+ mu_c = MIN(15., (1000.E6/Nt_c + 2.))
+
+!..Schmidt number to one-third used numerous times.
+ Sc3 = Sc**(1./3.)
+
+!..Compute min ice diam from mass, min snow/graupel mass from diam.
+ D0i = (xm0i/am_i)**(1./bm_i)
+ xm0s = am_s * D0s**bm_s
+ xm0g = am_g * D0g**bm_g
+
+!..These constants various exponents and gamma() assoc with cloud,
+!.. rain, snow, and graupel.
+ cce(1) = mu_c + 1.
+ cce(2) = bm_r + mu_c + 1.
+ cce(3) = bm_r + mu_c + 4.
+ ccg(1) = WGAMMA(cce(1))
+ ccg(2) = WGAMMA(cce(2))
+ ccg(3) = WGAMMA(cce(3))
+ ocg1 = 1./ccg(1)
+ ocg2 = 1./ccg(2)
+
+ cie(1) = mu_i + 1.
+ cie(2) = bm_i + mu_i + 1.
+ cie(3) = bm_i + mu_i + bv_i + 1.
+ cie(4) = mu_i + bv_i + 1.
+ cie(5) = mu_i + 2.
+ cie(6) = bm_i + bv_i
+ cig(1) = WGAMMA(cie(1))
+ cig(2) = WGAMMA(cie(2))
+ cig(3) = WGAMMA(cie(3))
+ cig(4) = WGAMMA(cie(4))
+ cig(5) = WGAMMA(cie(5))
+ cig(6) = WGAMMA(cie(6))
+ oig1 = 1./cig(1)
+ oig2 = 1./cig(2)
+ obmi = 1./bm_i
+
+ cre(1) = bm_r + 1.
+ cre(2) = mu_r + 1.
+ cre(3) = bm_r + mu_r + 1.
+ cre(4) = bm_r*2. + mu_r + 1.
+ cre(5) = mu_r + bv_r + 1.
+ cre(6) = bm_r + mu_r + bv_r + 1.
+ cre(7) = bm_r*0.5 + mu_r + bv_r + 1.
+ cre(8) = bm_r + mu_r + bv_r + 3.
+ cre(9) = mu_r + bv_r + 3.
+ cre(10) = mu_r + 2.
+ cre(11) = 0.5*(bv_r + 5. + 2.*mu_r)
+ cre(12) = bm_r*0.5 + mu_r + 1.
+ cre(13) = bm_r*2. + mu_r + bv_r + 1.
+ do n = 1, 13
+ crg(n) = WGAMMA(cre(n))
+ enddo
+ obmr = 1./bm_r
+ ore1 = 1./cre(1)
+ org1 = 1./crg(1)
+ org2 = 1./crg(2)
+ org3 = 1./crg(3)
+
+ cse(1) = bm_s + 1.
+ cse(2) = bm_s + 2.
+ cse(3) = bm_s*2.
+ cse(4) = bm_s + bv_s + 1.
+ cse(5) = bm_s*2. + bv_s + 1.
+ cse(6) = bm_s*2. + 1.
+ cse(7) = bm_s + mu_s + 1.
+ cse(8) = bm_s + mu_s + 2.
+ cse(9) = bm_s + mu_s + 3.
+ cse(10) = bm_s + mu_s + bv_s + 1.
+ cse(11) = bm_s*2. + mu_s + bv_s + 2.
+ cse(12) = bm_s*2. + mu_s + 1.
+ cse(13) = bv_s + 2.
+ cse(14) = bm_s + bv_s
+ cse(15) = mu_s + 1.
+ cse(16) = 1.0 + (1.0 + bv_s)/2.
+ cse(17) = cse(16) + mu_s + 1.
+ cse(18) = bv_s + mu_s + 3.
+ do n = 1, 18
+ csg(n) = WGAMMA(cse(n))
+ enddo
+ oams = 1./am_s
+ obms = 1./bm_s
+ ocms = oams**obms
+
+ cge(1) = bm_g + 1.
+ cge(2) = mu_g + 1.
+ cge(3) = bm_g + mu_g + 1.
+ cge(4) = bm_g*2. + mu_g + 1.
+ cge(5) = bm_g*2. + mu_g + bv_g + 1.
+ cge(6) = bm_g + mu_g + bv_g + 1.
+ cge(7) = bm_g + mu_g + bv_g + 2.
+ cge(8) = bm_g + mu_g + bv_g + 3.
+ cge(9) = mu_g + bv_g + 3.
+ cge(10) = mu_g + 2.
+ cge(11) = 0.5*(bv_g + 5. + 2.*mu_g)
+ cge(12) = 0.5*(bv_g + 5.) + mu_g
+ do n = 1, 12
+ cgg(n) = WGAMMA(cge(n))
+ enddo
+ oamg = 1./am_g
+ obmg = 1./bm_g
+ ocmg = oamg**obmg
+ oge1 = 1./cge(1)
+ ogg1 = 1./cgg(1)
+ ogg2 = 1./cgg(2)
+ ogg3 = 1./cgg(3)
+
+!+---+-----------------------------------------------------------------+
+!..Simplify various rate eqns the best we can now.
+!+---+-----------------------------------------------------------------+
+
+!..Rain collecting cloud water and cloud ice
+ t1_qr_qc = PI*.25*av_r * crg(9)
+ t1_qr_qi = PI*.25*av_r * crg(9)
+ t2_qr_qi = PI*.25*am_r*av_r * crg(8)
+
+!..Graupel collecting cloud water
+ t1_qg_qc = PI*.25*av_g * cgg(9)
+
+!..Snow collecting cloud water
+ t1_qs_qc = PI*.25*av_s
+
+!..Snow collecting cloud ice
+ t1_qs_qi = PI*.25*av_s
+
+!..Evaporation of rain; ignore depositional growth of rain.
+ t1_qr_ev = 0.78 * crg(10)
+ t2_qr_ev = 0.308*Sc3*SQRT(av_r) * crg(11)
+
+!..Sublimation/depositional growth of snow
+ t1_qs_sd = 0.86
+ t2_qs_sd = 0.28*Sc3*SQRT(av_s)
+
+!..Melting of snow
+ t1_qs_me = PI*4.*C_sqrd*olfus * 0.86
+ t2_qs_me = PI*4.*C_sqrd*olfus * 0.28*Sc3*SQRT(av_s)
+
+!..Sublimation/depositional growth of graupel
+ t1_qg_sd = 0.86 * cgg(10)
+ t2_qg_sd = 0.28*Sc3*SQRT(av_g) * cgg(11)
+
+!..Melting of graupel
+ t1_qg_me = PI*4.*C_cube*olfus * 0.86 * cgg(10)
+ t2_qg_me = PI*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g) * cgg(11)
+
+!..Constants for helping find lookup table indexes.
+ nic2 = NINT(ALOG10(r_c(1)))
+ nii2 = NINT(ALOG10(r_i(1)))
+ nii3 = NINT(ALOG10(Nt_i(1)))
+ nir2 = NINT(ALOG10(r_r(1)))
+ nir3 = NINT(ALOG10(N0r_exp(1)))
+ nis2 = NINT(ALOG10(r_s(1)))
+ nig2 = NINT(ALOG10(r_g(1)))
+ nig3 = NINT(ALOG10(N0g_exp(1)))
+
+!..Create bins of cloud water (from min diameter up to 100 microns).
+ Dc(1) = D0c*1.0d0
+ dtc(1) = D0c*1.0d0
+ do n = 2, nbc
+ Dc(n) = Dc(n-1) + 1.0D-6
+ dtc(n) = (Dc(n) - Dc(n-1))
+ enddo
+
+!..Create bins of cloud ice (from min diameter up to 5x min snow size).
+ xDx(1) = D0i*1.0d0
+ xDx(nbi+1) = 5.0d0*D0s
+ do n = 2, nbi
+ xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbi) &
+ *DLOG(xDx(nbi+1)/xDx(1)) +DLOG(xDx(1)))
+ enddo
+ do n = 1, nbi
+ Di(n) = DSQRT(xDx(n)*xDx(n+1))
+ dti(n) = xDx(n+1) - xDx(n)
+ enddo
+
+!..Create bins of rain (from min diameter up to 5 mm).
+ xDx(1) = D0r*1.0d0
+ xDx(nbr+1) = 0.005d0
+ do n = 2, nbr
+ xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbr) &
+ *DLOG(xDx(nbr+1)/xDx(1)) +DLOG(xDx(1)))
+ enddo
+ do n = 1, nbr
+ Dr(n) = DSQRT(xDx(n)*xDx(n+1))
+ dtr(n) = xDx(n+1) - xDx(n)
+ enddo
+
+!..Create bins of snow (from min diameter up to 2 cm).
+ xDx(1) = D0s*1.0d0
+ xDx(nbs+1) = 0.02d0
+ do n = 2, nbs
+ xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbs) &
+ *DLOG(xDx(nbs+1)/xDx(1)) +DLOG(xDx(1)))
+ enddo
+ do n = 1, nbs
+ Ds(n) = DSQRT(xDx(n)*xDx(n+1))
+ dts(n) = xDx(n+1) - xDx(n)
+ enddo
+
+!..Create bins of graupel (from min diameter up to 5 cm).
+ xDx(1) = D0g*1.0d0
+ xDx(nbg+1) = 0.05d0
+ do n = 2, nbg
+ xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbg) &
+ *DLOG(xDx(nbg+1)/xDx(1)) +DLOG(xDx(1)))
+ enddo
+ do n = 1, nbg
+ Dg(n) = DSQRT(xDx(n)*xDx(n+1))
+ dtg(n) = xDx(n+1) - xDx(n)
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..Create lookup tables for most costly calculations.
+!+---+-----------------------------------------------------------------+
+
+ do m = 1, ntb_r
+ do k = 1, ntb_r1
+ do j = 1, ntb_g
+ do i = 1, ntb_g1
+ tcg_racg(i,j,k,m) = 0.0d0
+ tmr_racg(i,j,k,m) = 0.0d0
+ tcr_gacr(i,j,k,m) = 0.0d0
+ tmg_gacr(i,j,k,m) = 0.0d0
+ tnr_racg(i,j,k,m) = 0.0d0
+ tnr_gacr(i,j,k,m) = 0.0d0
+ enddo
+ enddo
+ enddo
+ enddo
+
+ do m = 1, ntb_r
+ do k = 1, ntb_r1
+ do j = 1, ntb_t
+ do i = 1, ntb_s
+ tcs_racs1(i,j,k,m) = 0.0d0
+ tmr_racs1(i,j,k,m) = 0.0d0
+ tcs_racs2(i,j,k,m) = 0.0d0
+ tmr_racs2(i,j,k,m) = 0.0d0
+ tcr_sacr1(i,j,k,m) = 0.0d0
+ tms_sacr1(i,j,k,m) = 0.0d0
+ tcr_sacr2(i,j,k,m) = 0.0d0
+ tms_sacr2(i,j,k,m) = 0.0d0
+ tnr_racs1(i,j,k,m) = 0.0d0
+ tnr_racs2(i,j,k,m) = 0.0d0
+ tnr_sacr1(i,j,k,m) = 0.0d0
+ tnr_sacr2(i,j,k,m) = 0.0d0
+ enddo
+ enddo
+ enddo
+ enddo
+
+ do k = 1, 45
+ do j = 1, ntb_r1
+ do i = 1, ntb_r
+ tpi_qrfz(i,j,k) = 0.0d0
+ tni_qrfz(i,j,k) = 0.0d0
+ tpg_qrfz(i,j,k) = 0.0d0
+ tnr_qrfz(i,j,k) = 0.0d0
+ enddo
+ enddo
+ do i = 1, ntb_c
+ tpi_qcfz(i,k) = 0.0d0
+ tni_qcfz(i,k) = 0.0d0
+ enddo
+ enddo
+
+ do j = 1, ntb_i1
+ do i = 1, ntb_i
+ tps_iaus(i,j) = 0.0d0
+ tni_iaus(i,j) = 0.0d0
+ tpi_ide(i,j) = 0.0d0
+ enddo
+ enddo
+
+ do j = 1, nbc
+ do i = 1, nbr
+ t_Efrw(i,j) = 0.0
+ enddo
+ do i = 1, nbs
+ t_Efsw(i,j) = 0.0
+ enddo
+ enddo
+
+ do k = 1, ntb_r
+ do j = 1, ntb_r1
+ do i = 1, nbr
+ tnr_rev(i,j,k) = 0.0d0
+ enddo
+ enddo
+ enddo
+
+! CALL wrf_debug(150, 'CREATING MICROPHYSICS LOOKUP TABLES ... ')
+! WRITE (wrf_err_message, '(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') &
+! ' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g
+! CALL wrf_debug(150, wrf_err_message)
+
+!..Collision efficiency between rain/snow and cloud water.
+! CALL wrf_debug(200, ' creating qc collision eff tables')
+ call table_Efrw
+ call table_Efsw
+
+!..Drop evaporation.
+! CALL wrf_debug(200, ' creating rain evap table')
+! call table_dropEvap
+
+!..Initialize various constants for computing radar reflectivity.
+! call radar_init
+
+!LDF begin (05-13-2010): read pre-calculated look-up tables.
+ iiwarm = .false.
+ l_qr_acr_qg = .false.
+ l_qr_acr_qs = .false.
+ l_qi_aut_qs = .false.
+ l_freezeH2O = .false.
+
+ inquire(file='./LOOKUP_TABLES/table_qr_acr_qg.dat',exist=l_qr_acr_qg)
+ inquire(file='./LOOKUP_TABLES/table_qr_acr_qs.dat',exist=l_qr_acr_qs)
+ inquire(file='./LOOKUP_TABLES/table_qi_aut_qs.dat',exist=l_qi_aut_qs)
+ inquire(file='./LOOKUP_TABLES/table_freezeH2O.dat',exist=l_freezeH2O)
+
+ IF(l_qr_acr_qg .AND. l_qr_acr_qs .AND. l_qi_aut_qs .AND. &
+ l_freezeH2O) iiwarm = .true.
+
+ if(iiwarm) then
+ write(6,*) '--- BEGIN READ PRE-CALCULATED LOOK-UP TABLES'
+!..Rain collecting graupel & graupel collecting rain.
+ open(unit=11,file='./LOOKUP_TABLES/table_qr_acr_qg.dat', &
+ form='unformatted',status='old',readonly)
+ read(11) tcg_racg
+ read(11) tmr_racg
+ read(11) tcr_gacr
+ read(11) tmg_gacr
+ read(11) tnr_racg
+ read(11) tnr_gacr
+ close(unit=11)
+
+!..Rain collecting snow & snow collecting rain.
+ open(unit=11,file='./LOOKUP_TABLES/table_qr_acr_qs.dat', &
+ form='unformatted',status='old',readonly)
+ read(11) tcs_racs1
+ read(11) tmr_racs1
+ read(11) tcs_racs2
+ read(11) tmr_racs2
+ read(11) tcr_sacr1
+ read(11) tms_sacr1
+ read(11) tcr_sacr2
+ read(11) tms_sacr2
+ read(11) tnr_racs1
+ read(11) tnr_racs2
+ read(11) tnr_sacr1
+ read(11) tnr_sacr2
+ close(unit=11)
+
+!..Cloud water and rain freezing (Bigg, 1953).
+ open(unit=11,file='./LOOKUP_TABLES/table_freezeH2O.dat', &
+ form='unformatted',status='old',readonly)
+ read(11) tpi_qrfz
+ read(11) tni_qrfz
+ read(11) tpg_qrfz
+ read(11) tnr_qrfz
+ read(11) tpi_qcfz
+ read(11) tni_qcfz
+ close(unit=11)
+
+!..Conversion of some ice mass into snow category.
+ open(unit=11,file='./LOOKUP_TABLES/table_qi_aut_qs.dat', &
+ form='unformatted',status='old',readonly)
+ read(11) tpi_ide
+ read(11) tps_iaus
+ read(11) tni_iaus
+ close(unit=11)
+
+ write(6,*) '--- END PRE-CALCULATED LOOK-UP TABLES'
+ iiwarm = .false.
+
+ elseif (.not. iiwarm) then
+
+!..Rain collecting graupel & graupel collecting rain.
+! CALL wrf_debug(200, ' creating rain collecting graupel table')
+ call qr_acr_qg
+
+!..Rain collecting snow & snow collecting rain.
+! CALL wrf_debug(200, ' creating rain collecting snow table')
+ call qr_acr_qs
+
+!..Cloud water and rain freezing (Bigg, 1953).
+! CALL wrf_debug(200, ' creating freezing of water drops table')
+ call freezeH2O
+
+!..Conversion of some ice mass into snow category.
+! CALL wrf_debug(200, ' creating ice converting to snow table')
+ call qi_aut_qs
+
+ endif
+
+! CALL wrf_debug(150, ' ... DONE microphysical lookup tables')
+ endif
+
+ END SUBROUTINE thompson_init
+!+---+-----------------------------------------------------------------+
+!ctrlL
+!+---+-----------------------------------------------------------------+
+!..This is a wrapper routine designed to transfer values from 3D to 1D.
+!+---+-----------------------------------------------------------------+
+ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, &
+ th, pii, p, dz, dt_in, itimestep, &
+ RAINNC, RAINNCV, &
+ SNOWNC, SNOWNCV, &
+ GRAUPELNC, GRAUPELNCV, &
+ SR, &
+! refl_10cm, grid_clock, grid_alarms, &
+ ids,ide, jds,jde, kds,kde, & ! domain dims
+ ims,ime, jms,jme, kms,kme, & ! memory dims
+ its,ite, jts,jte, kts,kte) ! tile dims
+
+ implicit none
+
+!..Subroutine arguments
+ INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, &
+ ims,ime, jms,jme, kms,kme, &
+ its,ite, jts,jte, kts,kte
+ REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
+ qv, qc, qr, qi, qs, qg, ni, nr, th
+ REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
+ pii, p, dz
+ REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: &
+ RAINNC, RAINNCV, SR
+ REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(INOUT):: &
+ SNOWNC, SNOWNCV, GRAUPELNC, GRAUPELNCV
+! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
+! refl_10cm
+ REAL, INTENT(IN):: dt_in
+ INTEGER, INTENT(IN):: itimestep
+
+! TYPE (WRFU_Clock):: grid_clock
+! TYPE (WRFU_Alarm), POINTER:: grid_alarms(:)
+
+!..Local variables
+ REAL, DIMENSION(kts:kte):: &
+ qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
+ nr1d, t1d, p1d, dz1d, dBZ
+ REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
+ REAL:: dt, pptrain, pptsnow, pptgraul, pptice
+ REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max
+ INTEGER:: i, j, k
+ INTEGER:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr
+ INTEGER:: jmax_qc,jmax_qr,jmax_qi,jmax_qs,jmax_qg,jmax_ni,jmax_nr
+ INTEGER:: kmax_qc,kmax_qr,kmax_qi,kmax_qs,kmax_qg,kmax_ni,kmax_nr
+ INTEGER:: i_start, j_start, i_end, j_end
+ LOGICAL:: dBZ_tstep
+
+!+---+
+
+ dBZ_tstep = .false.
+! if ( Is_alarm_tstep(grid_clock, grid_alarms(HISTORY_ALARM)) ) then
+! dBZ_tstep = .true.
+! endif
+
+ i_start = its
+ j_start = jts
+ i_end = ite
+ j_end = jte
+! if ( (ite-its+1).gt.4 .and. (jte-jts+1).lt.4) then
+! i_start = its + 2
+! i_end = ite - 1
+! j_start = jts
+! j_end = jte
+! elseif ( (ite-its+1).lt.4 .and. (jte-jts+1).gt.4) then
+! i_start = its
+! i_end = ite
+! j_start = jts + 2
+! j_end = jte - 1
+! endif
+
+ dt = dt_in
+
+ qc_max = 0.
+ qr_max = 0.
+ qs_max = 0.
+ qi_max = 0.
+ qg_max = 0
+ ni_max = 0.
+ nr_max = 0.
+ imax_qc = 0
+ imax_qr = 0
+ imax_qi = 0
+ imax_qs = 0
+ imax_qg = 0
+ imax_ni = 0
+ imax_nr = 0
+ jmax_qc = 0
+ jmax_qr = 0
+ jmax_qi = 0
+ jmax_qs = 0
+ jmax_qg = 0
+ jmax_ni = 0
+ jmax_nr = 0
+ kmax_qc = 0
+ kmax_qr = 0
+ kmax_qi = 0
+ kmax_qs = 0
+ kmax_qg = 0
+ kmax_ni = 0
+ kmax_nr = 0
+ do i = 1, 256
+ mp_debug(i:i) = char(0)
+ enddo
+
+ j_loop: do j = j_start, j_end
+ i_loop: do i = i_start, i_end
+
+ pptrain = 0.
+ pptsnow = 0.
+ pptgraul = 0.
+ pptice = 0.
+ RAINNCV(i,j) = 0.
+ IF ( PRESENT (snowncv) ) THEN
+ SNOWNCV(i,j) = 0.
+ ENDIF
+ IF ( PRESENT (graupelncv) ) THEN
+ GRAUPELNCV(i,j) = 0.
+ ENDIF
+ SR(i,j) = 0.
+
+ do k = kts, kte
+ t1d(k) = th(i,k,j)*pii(i,k,j)
+ p1d(k) = p(i,k,j)
+ dz1d(k) = dz(i,k,j)
+ qv1d(k) = qv(i,k,j)
+ qc1d(k) = qc(i,k,j)
+ qi1d(k) = qi(i,k,j)
+ qr1d(k) = qr(i,k,j)
+ qs1d(k) = qs(i,k,j)
+ qg1d(k) = qg(i,k,j)
+ ni1d(k) = ni(i,k,j)
+ nr1d(k) = nr(i,k,j)
+ enddo
+
+ call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
+ nr1d, t1d, p1d, dz1d, &
+ pptrain, pptsnow, pptgraul, pptice, &
+ kts, kte, dt, i, j)
+
+ pcp_ra(i,j) = pptrain
+ pcp_sn(i,j) = pptsnow
+ pcp_gr(i,j) = pptgraul
+ pcp_ic(i,j) = pptice
+ RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice
+ RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
+ IF ( PRESENT(snowncv) .AND. PRESENT(snownc) ) THEN
+ SNOWNCV(i,j) = pptsnow + pptice
+ SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + pptice
+ ENDIF
+ IF ( PRESENT(graupelncv) .AND. PRESENT(graupelnc) ) THEN
+ GRAUPELNCV(i,j) = pptgraul
+ GRAUPELNC(i,j) = GRAUPELNC(i,j) + pptgraul
+ ENDIF
+ SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12)
+
+ do k = kts, kte
+ qv(i,k,j) = qv1d(k)
+ qc(i,k,j) = qc1d(k)
+ qi(i,k,j) = qi1d(k)
+ qr(i,k,j) = qr1d(k)
+ qs(i,k,j) = qs1d(k)
+ qg(i,k,j) = qg1d(k)
+ ni(i,k,j) = ni1d(k)
+ nr(i,k,j) = nr1d(k)
+ th(i,k,j) = t1d(k)/pii(i,k,j)
+ if (qc1d(k) .gt. qc_max) then
+ imax_qc = i
+ jmax_qc = j
+ kmax_qc = k
+ qc_max = qc1d(k)
+ elseif (qc1d(k) .lt. 0.0) then
+ write(mp_debug,*) 'WARNING, negative qc ', qc1d(k), &
+ ' at i,j,k=', i,j,k
+! CALL wrf_debug(150, mp_debug)
+ endif
+ if (qr1d(k) .gt. qr_max) then
+ imax_qr = i
+ jmax_qr = j
+ kmax_qr = k
+ qr_max = qr1d(k)
+ elseif (qr1d(k) .lt. 0.0) then
+ write(mp_debug,*) 'WARNING, negative qr ', qr1d(k), &
+ ' at i,j,k=', i,j,k
+! CALL wrf_debug(150, mp_debug)
+ endif
+ if (nr1d(k) .gt. nr_max) then
+ imax_nr = i
+ jmax_nr = j
+ kmax_nr = k
+ nr_max = nr1d(k)
+ elseif (nr1d(k) .lt. 0.0) then
+ write(mp_debug,*) 'WARNING, negative nr ', nr1d(k), &
+ ' at i,j,k=', i,j,k
+! CALL wrf_debug(150, mp_debug)
+ endif
+ if (qs1d(k) .gt. qs_max) then
+ imax_qs = i
+ jmax_qs = j
+ kmax_qs = k
+ qs_max = qs1d(k)
+ elseif (qs1d(k) .lt. 0.0) then
+ write(mp_debug,*) 'WARNING, negative qs ', qs1d(k), &
+ ' at i,j,k=', i,j,k
+! CALL wrf_debug(150, mp_debug)
+ endif
+ if (qi1d(k) .gt. qi_max) then
+ imax_qi = i
+ jmax_qi = j
+ kmax_qi = k
+ qi_max = qi1d(k)
+ elseif (qi1d(k) .lt. 0.0) then
+ write(mp_debug,*) 'WARNING, negative qi ', qi1d(k), &
+ ' at i,j,k=', i,j,k
+! CALL wrf_debug(150, mp_debug)
+ endif
+ if (qg1d(k) .gt. qg_max) then
+ imax_qg = i
+ jmax_qg = j
+ kmax_qg = k
+ qg_max = qg1d(k)
+ elseif (qg1d(k) .lt. 0.0) then
+ write(mp_debug,*) 'WARNING, negative qg ', qg1d(k), &
+ ' at i,j,k=', i,j,k
+! CALL wrf_debug(150, mp_debug)
+ endif
+ if (ni1d(k) .gt. ni_max) then
+ imax_ni = i
+ jmax_ni = j
+ kmax_ni = k
+ ni_max = ni1d(k)
+ elseif (ni1d(k) .lt. 0.0) then
+ write(mp_debug,*) 'WARNING, negative ni ', ni1d(k), &
+ ' at i,j,k=', i,j,k
+! CALL wrf_debug(150, mp_debug)
+ endif
+ if (qv1d(k) .lt. 0.0) then
+ if (k.lt.kte-2 .and. k.gt.kts+1) then
+ qv(i,k,j) = 0.5*(qv(i,k-1,j) + qv(i,k+1,j))
+ else
+ qv(i,k,j) = 1.E-7
+ endif
+ write(mp_debug,*) 'WARNING, negative qv ', qv1d(k), &
+ ' at i,j,k=', i,j,k
+! CALL wrf_debug(150, mp_debug)
+ endif
+ enddo
+
+! if (dBZ_tstep) then
+! call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, &
+! t1d, p1d, dBZ, kts, kte, i, j)
+! do k = kts, kte
+! refl_10cm(i,k,j) = MAX(-35., dBZ(k))
+! enddo
+! endif
+
+ enddo i_loop
+ enddo j_loop
+
+! DEBUG - GT
+ write(mp_debug,'(a,7(a,e13.6,1x,a,i3,a,i3,a,i3,a,1x))') 'MP-GT:', &
+ 'qc: ', qc_max, '(', imax_qc, ',', jmax_qc, ',', kmax_qc, ')', &
+ 'qr: ', qr_max, '(', imax_qr, ',', jmax_qr, ',', kmax_qr, ')', &
+ 'qi: ', qi_max, '(', imax_qi, ',', jmax_qi, ',', kmax_qi, ')', &
+ 'qs: ', qs_max, '(', imax_qs, ',', jmax_qs, ',', kmax_qs, ')', &
+ 'qg: ', qg_max, '(', imax_qg, ',', jmax_qg, ',', kmax_qg, ')', &
+ 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')', &
+ 'nr: ', nr_max, '(', imax_nr, ',', jmax_nr, ',', kmax_nr, ')'
+! CALL wrf_debug(150, mp_debug)
+! END DEBUG - GT
+
+ do i = 1, 256
+ mp_debug(i:i) = char(0)
+ enddo
+
+ END SUBROUTINE mp_gt_driver
+
+!+---+-----------------------------------------------------------------+
+!ctrlL
+!+---+-----------------------------------------------------------------+
+!+---+-----------------------------------------------------------------+
+!.. This subroutine computes the moisture tendencies of water vapor,
+!.. cloud droplets, rain, cloud ice (pristine), snow, and graupel.
+!.. Previously this code was based on Reisner et al (1998), but few of
+!.. those pieces remain. A complete description is now found in
+!.. Thompson et al. (2004, 2008).
+!+---+-----------------------------------------------------------------+
+!
+ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
+ nr1d, t1d, p1d, dzq, &
+ pptrain, pptsnow, pptgraul, pptice, &
+ kts, kte, dt, ii, jj)
+
+ implicit none
+
+!..Sub arguments
+ INTEGER, INTENT(IN):: kts, kte, ii, jj
+ REAL, DIMENSION(kts:kte), INTENT(INOUT):: &
+ qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
+ nr1d, t1d, p1d
+ REAL, DIMENSION(kts:kte), INTENT(IN):: dzq
+ REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice
+ REAL, INTENT(IN):: dt
+
+!..Local variables
+ REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, &
+ qrten, qsten, qgten, niten, nrten
+
+ DOUBLE PRECISION, DIMENSION(kts:kte):: prw_vcd
+
+ DOUBLE PRECISION, DIMENSION(kts:kte):: prr_wau, prr_rcw, prr_rcs, &
+ prr_rcg, prr_sml, prr_gml, &
+ prr_rci, prv_rev, &
+ pnr_wau, pnr_rcs, pnr_rcg, &
+ pnr_rci, pnr_sml, pnr_gml, &
+ pnr_rev, pnr_rcr, pnr_rfz
+
+ DOUBLE PRECISION, DIMENSION(kts:kte):: pri_inu, pni_inu, pri_ihm, &
+ pni_ihm, pri_wfz, pni_wfz, &
+ pri_rfz, pni_rfz, pri_ide, &
+ pni_ide, pri_rci, pni_rci, &
+ pni_sci, pni_iau
+
+ DOUBLE PRECISION, DIMENSION(kts:kte):: prs_iau, prs_sci, prs_rcs, &
+ prs_scw, prs_sde, prs_ihm, &
+ prs_ide
+
+ DOUBLE PRECISION, DIMENSION(kts:kte):: prg_scw, prg_rfz, prg_gde, &
+ prg_gcw, prg_rci, prg_rcs, &
+ prg_rcg, prg_ihm
+
+ REAL, DIMENSION(kts:kte):: temp, pres, qv
+ REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni, nr
+ REAL, DIMENSION(kts:kte):: rho, rhof, rhof2
+ REAL, DIMENSION(kts:kte):: qvs, qvsi
+ REAL, DIMENSION(kts:kte):: satw, sati, ssatw, ssati
+ REAL, DIMENSION(kts:kte):: diffu, visco, vsc2, &
+ tcond, lvap, ocp, lvt2
+
+ DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g
+ REAL, DIMENSION(kts:kte):: mvd_r, mvd_c
+ REAL, DIMENSION(kts:kte):: smob, smo2, smo1, smo0, &
+ smoc, smod, smoe, smof
+
+ REAL, DIMENSION(kts:kte):: sed_r, sed_s, sed_g, sed_i, sed_n
+
+ REAL:: rgvm, delta_tp, orho, lfus2
+ REAL, DIMENSION(4):: onstep
+ DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamc, lamr, lamg
+ DOUBLE PRECISION:: lami, ilami
+ REAL:: xDc, Dc_b, Dc_g, xDi, xDr, xDs, xDg, Ds_m, Dg_m
+ DOUBLE PRECISION:: Dr_star
+ REAL:: zeta1, zeta, taud, tau
+ REAL:: stoke_r, stoke_s, stoke_g, stoke_i
+ REAL:: vti, vtr, vts, vtg
+ REAL, DIMENSION(kts:kte+1):: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk
+ REAL, DIMENSION(kts:kte):: vts_boost
+ REAL:: Mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, C_snow
+ REAL:: a_, b_, loga_, A1, A2, tf
+ REAL:: tempc, tc0, r_mvd1, r_mvd2, xkrat
+ REAL:: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr
+ REAL:: xsat, rate_max, sump, ratio
+ REAL:: clap, fcd, dfcd
+ REAL:: otemp, rvs, rvs_p, rvs_pp, gamsc, alphsc, t1_evap, t1_subl
+ REAL:: r_frac, g_frac
+ REAL:: Ef_rw, Ef_sw, Ef_gw, Ef_rr
+ REAL:: dtsave, odts, odt, odzq
+ INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq
+ INTEGER, DIMENSION(4):: ksed1
+ INTEGER:: nir, nis, nig, nii, nic
+ INTEGER:: idx_tc, idx_t, idx_s, idx_g1, idx_g, idx_r1, idx_r, &
+ idx_i1, idx_i, idx_c, idx, idx_d
+ LOGICAL:: melti, no_micro
+ LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg
+ LOGICAL:: debug_flag
+
+!+---+
+
+ debug_flag = .false.
+! if (ii.eq.319 .and. jj.eq.39) debug_flag = .true.
+
+ no_micro = .true.
+ dtsave = dt
+ odt = 1./dt
+ odts = 1./dtsave
+ iexfrq = 1
+
+!+---+-----------------------------------------------------------------+
+!.. Source/sink terms. First 2 chars: "pr" represents source/sink of
+!.. mass while "pn" represents source/sink of number. Next char is one
+!.. of "v" for water vapor, "r" for rain, "i" for cloud ice, "w" for
+!.. cloud water, "s" for snow, and "g" for graupel. Next chars
+!.. represent processes: "de" for sublimation/deposition, "ev" for
+!.. evaporation, "fz" for freezing, "ml" for melting, "au" for
+!.. autoconversion, "nu" for ice nucleation, "hm" for Hallet/Mossop
+!.. secondary ice production, and "c" for collection followed by the
+!.. character for the species being collected. ALL of these terms are
+!.. positive (except for deposition/sublimation terms which can switch
+!.. signs based on super/subsaturation) and are treated as negatives
+!.. where necessary in the tendency equations.
+!+---+-----------------------------------------------------------------+
+
+ do k = kts, kte
+ tten(k) = 0.
+ qvten(k) = 0.
+ qcten(k) = 0.
+ qiten(k) = 0.
+ qrten(k) = 0.
+ qsten(k) = 0.
+ qgten(k) = 0.
+ niten(k) = 0.
+ nrten(k) = 0.
+
+ prw_vcd(k) = 0.
+
+ prv_rev(k) = 0.
+ prr_wau(k) = 0.
+ prr_rcw(k) = 0.
+ prr_rcs(k) = 0.
+ prr_rcg(k) = 0.
+ prr_sml(k) = 0.
+ prr_gml(k) = 0.
+ prr_rci(k) = 0.
+ pnr_wau(k) = 0.
+ pnr_rcs(k) = 0.
+ pnr_rcg(k) = 0.
+ pnr_rci(k) = 0.
+ pnr_sml(k) = 0.
+ pnr_gml(k) = 0.
+ pnr_rev(k) = 0.
+ pnr_rcr(k) = 0.
+ pnr_rfz(k) = 0.
+
+ pri_inu(k) = 0.
+ pni_inu(k) = 0.
+ pri_ihm(k) = 0.
+ pni_ihm(k) = 0.
+ pri_wfz(k) = 0.
+ pni_wfz(k) = 0.
+ pri_rfz(k) = 0.
+ pni_rfz(k) = 0.
+ pri_ide(k) = 0.
+ pni_ide(k) = 0.
+ pri_rci(k) = 0.
+ pni_rci(k) = 0.
+ pni_sci(k) = 0.
+ pni_iau(k) = 0.
+
+ prs_iau(k) = 0.
+ prs_sci(k) = 0.
+ prs_rcs(k) = 0.
+ prs_scw(k) = 0.
+ prs_sde(k) = 0.
+ prs_ihm(k) = 0.
+ prs_ide(k) = 0.
+
+ prg_scw(k) = 0.
+ prg_rfz(k) = 0.
+ prg_gde(k) = 0.
+ prg_gcw(k) = 0.
+ prg_rci(k) = 0.
+ prg_rcs(k) = 0.
+ prg_rcg(k) = 0.
+ prg_ihm(k) = 0.
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..Put column of data into local arrays.
+!+---+-----------------------------------------------------------------+
+ do k = kts, kte
+ temp(k) = t1d(k)
+ qv(k) = MAX(1.E-10, qv1d(k))
+ pres(k) = p1d(k)
+ rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
+ if (qc1d(k) .gt. R1) then
+ no_micro = .false.
+ rc(k) = qc1d(k)*rho(k)
+ L_qc(k) = .true.
+ else
+ qc1d(k) = 0.0
+ rc(k) = R1
+ L_qc(k) = .false.
+ endif
+ if (qi1d(k) .gt. R1) then
+ no_micro = .false.
+ ri(k) = qi1d(k)*rho(k)
+ ni(k) = MAX(1., ni1d(k)*rho(k))
+ L_qi(k) = .true.
+ lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
+ ilami = 1./lami
+ xDi = (bm_i + mu_i + 1.) * ilami
+ if (xDi.lt. 20.E-6) then
+ lami = cie(2)/20.E-6
+ ni(k) = MIN(500.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
+ elseif (xDi.gt. 300.E-6) then
+ lami = cie(2)/300.E-6
+ ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i
+ endif
+ else
+ qi1d(k) = 0.0
+ ni1d(k) = 0.0
+ ri(k) = R1
+ ni(k) = 0.01
+ L_qi(k) = .false.
+ endif
+
+ if (qr1d(k) .gt. R1) then
+ no_micro = .false.
+ rr(k) = qr1d(k)*rho(k)
+ nr(k) = MAX(1., nr1d(k)*rho(k))
+ L_qr(k) = .true.
+ if (nr(k) .gt. 1.0) then
+ lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
+ mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
+ if (mvd_r(k) .gt. 2.5E-3) then
+ mvd_r(k) = 2.5E-3
+ lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
+ nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
+ elseif (mvd_r(k) .lt. D0r*0.75) then
+ mvd_r(k) = D0r*0.75
+ lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
+ nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
+ endif
+ else
+ if (qr1d(k) .gt. R2) then
+ mvd_r(k) = 2.5E-3
+ else
+ mvd_r(k) = 2.5E-3 / 3.0**(ALOG10(R2)-ALOG10(qr1d(k)))
+ endif
+ lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
+ nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
+ endif
+ else
+ qr1d(k) = 0.0
+ nr1d(k) = 0.0
+ rr(k) = R1
+ nr(k) = 1.0
+ L_qr(k) = .false.
+ endif
+ if (qs1d(k) .gt. R1) then
+ no_micro = .false.
+ rs(k) = qs1d(k)*rho(k)
+ L_qs(k) = .true.
+ else
+ qs1d(k) = 0.0
+ rs(k) = R1
+ L_qs(k) = .false.
+ endif
+ if (qg1d(k) .gt. R1) then
+ no_micro = .false.
+ rg(k) = qg1d(k)*rho(k)
+ L_qg(k) = .true.
+ else
+ qg1d(k) = 0.0
+ rg(k) = R1
+ L_qg(k) = .false.
+ endif
+ enddo
+
+
+!+---+-----------------------------------------------------------------+
+!..Derive various thermodynamic variables frequently used.
+!.. Saturation vapor pressure (mixing ratio) over liquid/ice comes from
+!.. Flatau et al. 1992; enthalpy (latent heat) of vaporization from
+!.. Bohren & Albrecht 1998; others from Pruppacher & Klett 1978.
+!+---+-----------------------------------------------------------------+
+ do k = kts, kte
+ tempc = temp(k) - 273.15
+ rhof(k) = SQRT(RHO_NOT/rho(k))
+ rhof2(k) = SQRT(rhof(k))
+ qvs(k) = rslf(pres(k), temp(k))
+ if (tempc .le. 0.0) then
+ qvsi(k) = rsif(pres(k), temp(k))
+ else
+ qvsi(k) = qvs(k)
+ endif
+ satw(k) = qv(k)/qvs(k)
+ sati(k) = qv(k)/qvsi(k)
+ ssatw(k) = satw(k) - 1.
+ ssati(k) = sati(k) - 1.
+ if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
+ if (abs(ssati(k)).lt. eps) ssati(k) = 0.0
+ if (no_micro .and. ssati(k).gt. 0.0) no_micro = .false.
+ diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
+ if (tempc .ge. 0.0) then
+ visco(k) = (1.718+0.0049*tempc)*1.0E-5
+ else
+ visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
+ endif
+ ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
+ vsc2(k) = SQRT(rho(k)/visco(k))
+ lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
+ tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..If no existing hydrometeor species and no chance to initiate ice or
+!.. condense cloud water, just exit quickly!
+!+---+-----------------------------------------------------------------+
+
+ if (no_micro) return
+
+!+---+-----------------------------------------------------------------+
+!..Calculate y-intercept, slope, and useful moments for snow.
+!+---+-----------------------------------------------------------------+
+ if (.not. iiwarm) then
+ do k = kts, kte
+ if (.not. L_qs(k)) CYCLE
+ tc0 = MIN(-0.1, temp(k)-273.15)
+ smob(k) = rs(k)*oams
+
+!..All other moments based on reference, 2nd moment. If bm_s.ne.2,
+!.. then we must compute actual 2nd moment and use as reference.
+ if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
+ smo2(k) = smob(k)
+ else
+ loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
+ + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
+ + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
+ + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
+ + sa(10)*bm_s*bm_s*bm_s
+ a_ = 10.0**loga_
+ b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
+ + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
+ + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
+ + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
+ + sb(10)*bm_s*bm_s*bm_s
+ smo2(k) = (smob(k)/a_)**(1./b_)
+ endif
+
+!..Calculate 0th moment. Represents snow number concentration.
+ loga_ = sa(1) + sa(2)*tc0 + sa(5)*tc0*tc0 + sa(9)*tc0*tc0*tc0
+ a_ = 10.0**loga_
+ b_ = sb(1) + sb(2)*tc0 + sb(5)*tc0*tc0 + sb(9)*tc0*tc0*tc0
+ smo0(k) = a_ * smo2(k)**b_
+
+!..Calculate 1st moment. Useful for depositional growth and melting.
+ loga_ = sa(1) + sa(2)*tc0 + sa(3) &
+ + sa(4)*tc0 + sa(5)*tc0*tc0 &
+ + sa(6) + sa(7)*tc0*tc0 &
+ + sa(8)*tc0 + sa(9)*tc0*tc0*tc0 &
+ + sa(10)
+ a_ = 10.0**loga_
+ b_ = sb(1)+ sb(2)*tc0 + sb(3) + sb(4)*tc0 &
+ + sb(5)*tc0*tc0 + sb(6) &
+ + sb(7)*tc0*tc0 + sb(8)*tc0 &
+ + sb(9)*tc0*tc0*tc0 + sb(10)
+ smo1(k) = a_ * smo2(k)**b_
+
+!..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
+ loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
+ + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
+ + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
+ + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
+ + sa(10)*cse(1)*cse(1)*cse(1)
+ a_ = 10.0**loga_
+ b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
+ + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
+ + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
+ + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
+ smoc(k) = a_ * smo2(k)**b_
+
+!..Calculate bv_s+2 (th) moment. Useful for riming.
+ loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(13) &
+ + sa(4)*tc0*cse(13) + sa(5)*tc0*tc0 &
+ + sa(6)*cse(13)*cse(13) + sa(7)*tc0*tc0*cse(13) &
+ + sa(8)*tc0*cse(13)*cse(13) + sa(9)*tc0*tc0*tc0 &
+ + sa(10)*cse(13)*cse(13)*cse(13)
+ a_ = 10.0**loga_
+ b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(13) + sb(4)*tc0*cse(13) &
+ + sb(5)*tc0*tc0 + sb(6)*cse(13)*cse(13) &
+ + sb(7)*tc0*tc0*cse(13) + sb(8)*tc0*cse(13)*cse(13) &
+ + sb(9)*tc0*tc0*tc0 + sb(10)*cse(13)*cse(13)*cse(13)
+ smoe(k) = a_ * smo2(k)**b_
+
+!..Calculate 1+(bv_s+1)/2 (th) moment. Useful for depositional growth.
+ loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(16) &
+ + sa(4)*tc0*cse(16) + sa(5)*tc0*tc0 &
+ + sa(6)*cse(16)*cse(16) + sa(7)*tc0*tc0*cse(16) &
+ + sa(8)*tc0*cse(16)*cse(16) + sa(9)*tc0*tc0*tc0 &
+ + sa(10)*cse(16)*cse(16)*cse(16)
+ a_ = 10.0**loga_
+ b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(16) + sb(4)*tc0*cse(16) &
+ + sb(5)*tc0*tc0 + sb(6)*cse(16)*cse(16) &
+ + sb(7)*tc0*tc0*cse(16) + sb(8)*tc0*cse(16)*cse(16) &
+ + sb(9)*tc0*tc0*tc0 + sb(10)*cse(16)*cse(16)*cse(16)
+ smof(k) = a_ * smo2(k)**b_
+
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..Calculate y-intercept, slope values for graupel.
+!+---+-----------------------------------------------------------------+
+ do k = kte, kts, -1
+ N0_exp = (gonv_max-gonv_min)*0.5D0 &
+ * tanh((0.01E-3-(rc(k)+rr(k)))/0.75E-3) &
+ + (gonv_max+gonv_min)*0.5D0
+! N0_exp = (gonv_max-gonv_min)*0.5D0 &
+! * tanh((-15.-(temp(k)-273.15))/7.5) &
+! + (gonv_max+gonv_min)*0.5D0
+ lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
+ lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
+ ilamg(k) = 1./lamg
+ N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2)
+ enddo
+
+ endif
+
+!+---+-----------------------------------------------------------------+
+!..Calculate y-intercept, slope values for rain.
+!+---+-----------------------------------------------------------------+
+ do k = kte, kts, -1
+ lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
+ ilamr(k) = 1./lamr
+ mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
+ N0_r(k) = nr(k)*org2*lamr**cre(2)
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..Compute warm-rain process terms (except evap done later).
+!+---+-----------------------------------------------------------------+
+
+ do k = kts, kte
+
+!..Rain self-collection follows Seifert, 1994 and drop break-up
+!.. follows Verlinde and Cotton, 1993. RAIN2M
+ if (L_qr(k) .and. mvd_r(k).gt. D0r) then
+ if (mvd_r(k) .le. 1750.0E-6) then
+ Ef_rr = 1.0
+ else
+ Ef_rr = 2.0 - EXP(2300.0*(mvd_r(k)-1750.0E-6))
+ endif
+ pnr_rcr(k) = Ef_rr * 8.*nr(k)*rr(k)
+ endif
+
+ if (.not. L_qc(k)) CYCLE
+ xDc = MAX(D0c*1.E6, ((rc(k)/(am_r*Nt_c))**obmr) * 1.E6)
+ lamc = (Nt_c*am_r* ccg(2) * ocg1 / rc(k))**obmr
+ mvd_c(k) = (3.0+mu_c+0.672) / lamc
+
+!..Autoconversion follows Berry & Reinhardt (1974) with characteristic
+!.. diameters correctly computed from gamma distrib of cloud droplets.
+ if (rc(k).gt. 0.01e-3) then
+ Dc_g = ((ccg(3)*ocg2)**obmr / lamc) * 1.E6
+ Dc_b = (xDc*xDc*xDc*Dc_g*Dc_g*Dc_g - xDc*xDc*xDc*xDc*xDc*xDc) &
+ **(1./6.)
+ zeta1 = 0.5*((6.25E-6*xDc*Dc_b*Dc_b*Dc_b - 0.4) &
+ + abs(6.25E-6*xDc*Dc_b*Dc_b*Dc_b - 0.4))
+ zeta = 0.027*rc(k)*zeta1
+ taud = 0.5*((0.5*Dc_b - 7.5) + abs(0.5*Dc_b - 7.5)) + R1
+ tau = 3.72/(rc(k)*taud)
+ prr_wau(k) = zeta/tau
+ prr_wau(k) = MIN(DBLE(rc(k)*odts), prr_wau(k))
+ pnr_wau(k) = prr_wau(k) / (am_r*mu_c/3.*D0r*D0r*D0r/rho(k)) ! RAIN2M
+ endif
+
+!..Rain collecting cloud water. In CE, assume Dc<<Dr and vtc=~0.
+ if (L_qr(k) .and. mvd_r(k).gt. D0r .and. mvd_c(k).gt. D0c) then
+ lamr = 1./ilamr(k)
+ idx = 1 + INT(nbr*DLOG(mvd_r(k)/Dr(1))/DLOG(Dr(nbr)/Dr(1)))
+ idx = MIN(idx, nbr)
+ Ef_rw = t_Efrw(idx, INT(mvd_c(k)*1.E6))
+ prr_rcw(k) = rhof(k)*t1_qr_qc*Ef_rw*rc(k)*N0_r(k) &
+ *((lamr+fv_r)**(-cre(9)))
+ prr_rcw(k) = MIN(DBLE(rc(k)*odts), prr_rcw(k))
+ endif
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..Compute all frozen hydrometeor species' process terms.
+!+---+-----------------------------------------------------------------+
+ if (.not. iiwarm) then
+ do k = kts, kte
+ vts_boost(k) = 1.5
+
+!..Temperature lookup table indexes.
+ tempc = temp(k) - 273.15
+ idx_tc = MAX(1, MIN(NINT(-tempc), 45) )
+ idx_t = INT( (tempc-2.5)/5. ) - 1
+ idx_t = MAX(1, -idx_t)
+ idx_t = MIN(idx_t, ntb_t)
+ IT = MAX(1, MIN(NINT(-tempc), 31) )
+
+!..Cloud water lookup table index.
+ if (rc(k).gt. r_c(1)) then
+ nic = NINT(ALOG10(rc(k)))
+ do nn = nic-1, nic+1
+ n = nn
+ if ( (rc(k)/10.**nn).ge.1.0 .and. &
+ (rc(k)/10.**nn).lt.10.0) goto 141
+ enddo
+ 141 continue
+ idx_c = INT(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
+ idx_c = MAX(1, MIN(idx_c, ntb_c))
+ else
+ idx_c = 1
+ endif
+
+!..Cloud ice lookup table indexes.
+ if (ri(k).gt. r_i(1)) then
+ nii = NINT(ALOG10(ri(k)))
+ do nn = nii-1, nii+1
+ n = nn
+ if ( (ri(k)/10.**nn).ge.1.0 .and. &
+ (ri(k)/10.**nn).lt.10.0) goto 142
+ enddo
+ 142 continue
+ idx_i = INT(ri(k)/10.**n) + 10*(n-nii2) - (n-nii2)
+ idx_i = MAX(1, MIN(idx_i, ntb_i))
+ else
+ idx_i = 1
+ endif
+
+ if (ni(k).gt. Nt_i(1)) then
+ nii = NINT(ALOG10(ni(k)))
+ do nn = nii-1, nii+1
+ n = nn
+ if ( (ni(k)/10.**nn).ge.1.0 .and. &
+ (ni(k)/10.**nn).lt.10.0) goto 143
+ enddo
+ 143 continue
+ idx_i1 = INT(ni(k)/10.**n) + 10*(n-nii3) - (n-nii3)
+ idx_i1 = MAX(1, MIN(idx_i1, ntb_i1))
+ else
+ idx_i1 = 1
+ endif
+
+!..Rain lookup table indexes.
+ if (rr(k).gt. r_r(1)) then
+ nir = NINT(ALOG10(rr(k)))
+ do nn = nir-1, nir+1
+ n = nn
+ if ( (rr(k)/10.**nn).ge.1.0 .and. &
+ (rr(k)/10.**nn).lt.10.0) goto 144
+ enddo
+ 144 continue
+ idx_r = INT(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2)
+ idx_r = MAX(1, MIN(idx_r, ntb_r))
+
+ lamr = 1./ilamr(k)
+ lam_exp = lamr * (crg(3)*org2*org1)**bm_r
+ N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
+ nir = NINT(DLOG10(N0_exp))
+ do nn = nir-1, nir+1
+ n = nn
+ if ( (N0_exp/10.**nn).ge.1.0 .and. &
+ (N0_exp/10.**nn).lt.10.0) goto 145
+ enddo
+ 145 continue
+ idx_r1 = INT(N0_exp/10.**n) + 10*(n-nir3) - (n-nir3)
+ idx_r1 = MAX(1, MIN(idx_r1, ntb_r1))
+ else
+ idx_r = 1
+ idx_r1 = ntb_r1
+ endif
+
+!..Snow lookup table index.
+ if (rs(k).gt. r_s(1)) then
+ nis = NINT(ALOG10(rs(k)))
+ do nn = nis-1, nis+1
+ n = nn
+ if ( (rs(k)/10.**nn).ge.1.0 .and. &
+ (rs(k)/10.**nn).lt.10.0) goto 146
+ enddo
+ 146 continue
+ idx_s = INT(rs(k)/10.**n) + 10*(n-nis2) - (n-nis2)
+ idx_s = MAX(1, MIN(idx_s, ntb_s))
+ else
+ idx_s = 1
+ endif
+
+!..Graupel lookup table index.
+ if (rg(k).gt. r_g(1)) then
+ nig = NINT(ALOG10(rg(k)))
+ do nn = nig-1, nig+1
+ n = nn
+ if ( (rg(k)/10.**nn).ge.1.0 .and. &
+ (rg(k)/10.**nn).lt.10.0) goto 147
+ enddo
+ 147 continue
+ idx_g = INT(rg(k)/10.**n) + 10*(n-nig2) - (n-nig2)
+ idx_g = MAX(1, MIN(idx_g, ntb_g))
+
+ lamg = 1./ilamg(k)
+ lam_exp = lamg * (cgg(3)*ogg2*ogg1)**bm_g
+ N0_exp = ogg1*rg(k)/am_g * lam_exp**cge(1)
+ nig = NINT(DLOG10(N0_exp))
+ do nn = nig-1, nig+1
+ n = nn
+ if ( (N0_exp/10.**nn).ge.1.0 .and. &
+ (N0_exp/10.**nn).lt.10.0) goto 148
+ enddo
+ 148 continue
+ idx_g1 = INT(N0_exp/10.**n) + 10*(n-nig3) - (n-nig3)
+ idx_g1 = MAX(1, MIN(idx_g1, ntb_g1))
+ else
+ idx_g = 1
+ idx_g1 = ntb_g1
+ endif
+
+!..Deposition/sublimation prefactor (from Srivastava & Coen 1992).
+ otemp = 1./temp(k)
+ rvs = rho(k)*qvsi(k)
+ rvs_p = rvs*otemp*(lsub*otemp*oRv - 1.)
+ rvs_pp = rvs * ( otemp*(lsub*otemp*oRv - 1.) &
+ *otemp*(lsub*otemp*oRv - 1.) &
+ + (-2.*lsub*otemp*otemp*otemp*oRv) &
+ + otemp*otemp)
+ gamsc = lsub*diffu(k)/tcond(k) * rvs_p
+ alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
+ * rvs_pp/rvs_p * rvs/rvs_p
+ alphsc = MAX(1.E-9, alphsc)
+ xsat = ssati(k)
+ if (abs(xsat).lt. 1.E-9) xsat=0.
+ t1_subl = 4.*PI*( 1.0 - alphsc*xsat &
+ + 2.*alphsc*alphsc*xsat*xsat &
+ - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
+ / (1.+gamsc)
+
+!..Snow collecting cloud water. In CE, assume Dc<<Ds and vtc=~0.
+ if (L_qc(k) .and. mvd_c(k).gt. D0c) then
+ xDs = 0.0
+ if (L_qs(k)) xDs = smoc(k) / smob(k)
+ if (xDs .gt. D0s) then
+ idx = 1 + INT(nbs*DLOG(xDs/Ds(1))/DLOG(Ds(nbs)/Ds(1)))
+ idx = MIN(idx, nbs)
+ Ef_sw = t_Efsw(idx, INT(mvd_c(k)*1.E6))
+ prs_scw(k) = rhof(k)*t1_qs_qc*Ef_sw*rc(k)*smoe(k)
+ endif
+
+!..Graupel collecting cloud water. In CE, assume Dc<<Dg and vtc=~0.
+ if (rg(k).ge. r_g(1) .and. mvd_c(k).gt. D0c) then
+ xDg = (bm_g + mu_g + 1.) * ilamg(k)
+ vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
+ stoke_g = mvd_c(k)*mvd_c(k)*vtg*rho_w/(9.*visco(k)*xDg)
+ if (xDg.gt. D0g) then
+ if (stoke_g.ge.0.4 .and. stoke_g.le.10.) then
+ Ef_gw = 0.55*ALOG10(2.51*stoke_g)
+ elseif (stoke_g.lt.0.4) then
+ Ef_gw = 0.0
+ elseif (stoke_g.gt.10) then
+ Ef_gw = 0.77
+ endif
+ prg_gcw(k) = rhof(k)*t1_qg_qc*Ef_gw*rc(k)*N0_g(k) &
+ *ilamg(k)**cge(9)
+ endif
+ endif
+ endif
+
+!..Rain collecting snow. Cannot assume Wisner (1972) approximation
+!.. or Mizuno (1990) approach so we solve the CE explicitly and store
+!.. results in lookup table.
+ if (rr(k).ge. r_r(1)) then
+ if (rs(k).ge. r_s(1)) then
+ if (temp(k).lt.T_0) then
+ prr_rcs(k) = -(tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
+ + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
+ + tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
+ + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r))
+ prs_rcs(k) = tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
+ + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
+ - tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
+ - tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
+ prg_rcs(k) = tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
+ + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
+ + tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
+ + tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
+ prr_rcs(k) = MAX(DBLE(-rr(k)*odts), prr_rcs(k))
+ prs_rcs(k) = MAX(DBLE(-rs(k)*odts), prs_rcs(k))
+ prg_rcs(k) = MIN(DBLE((rr(k)+rs(k))*odts), prg_rcs(k))
+ pnr_rcs(k) = tnr_racs1(idx_s,idx_t,idx_r1,idx_r) & ! RAIN2M
+ + tnr_racs2(idx_s,idx_t,idx_r1,idx_r) &
+ + tnr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
+ + tnr_sacr2(idx_s,idx_t,idx_r1,idx_r)
+ else
+ prs_rcs(k) = -(tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
+ + tcs_racs2(idx_s,idx_t,idx_r1,idx_r))
+ prs_rcs(k) = MAX(DBLE(-rs(k)*odts), prs_rcs(k))
+ prr_rcs(k) = -prs_rcs(k)
+ pnr_rcs(k) = tnr_racs2(idx_s,idx_t,idx_r1,idx_r) & ! RAIN2M
+ + tnr_sacr2(idx_s,idx_t,idx_r1,idx_r)
+ endif
+ pnr_rcs(k) = MIN(DBLE(nr(k)*odts), pnr_rcs(k))
+ endif
+
+!..Rain collecting graupel. Cannot assume Wisner (1972) approximation
+!.. or Mizuno (1990) approach so we solve the CE explicitly and store
+!.. results in lookup table.
+ if (rg(k).ge. r_g(1)) then
+ if (temp(k).lt.T_0) then
+ prg_rcg(k) = tmr_racg(idx_g1,idx_g,idx_r1,idx_r) &
+ + tcr_gacr(idx_g1,idx_g,idx_r1,idx_r)
+ prg_rcg(k) = MIN(DBLE(rr(k)*odts), prg_rcg(k))
+ prr_rcg(k) = -prg_rcg(k)
+ pnr_rcg(k) = tnr_racg(idx_g1,idx_g,idx_r1,idx_r) & ! RAIN2M
+ + tnr_gacr(idx_g1,idx_g,idx_r1,idx_r)
+ pnr_rcg(k) = MIN(DBLE(nr(k)*odts), pnr_rcg(k))
+ else
+ prr_rcg(k) = tcg_racg(idx_g1,idx_g,idx_r1,idx_r)
+ prr_rcg(k) = MIN(DBLE(rg(k)*odts), prr_rcg(k))
+ prg_rcg(k) = -prr_rcg(k)
+ endif
+ endif
+ endif
+
+!+---+-----------------------------------------------------------------+
+!..Next IF block handles only those processes below 0C.
+!+---+-----------------------------------------------------------------+
+
+ if (temp(k).lt.T_0) then
+
+ vts_boost(k) = 1.0
+ rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
+
+!..Freezing of water drops into graupel/cloud ice (Bigg 1953).
+ if (rr(k).gt. r_r(1)) then
+ prg_rfz(k) = tpg_qrfz(idx_r,idx_r1,idx_tc)*odts
+ pri_rfz(k) = tpi_qrfz(idx_r,idx_r1,idx_tc)*odts
+ pni_rfz(k) = tni_qrfz(idx_r,idx_r1,idx_tc)*odts
+ pnr_rfz(k) = tnr_qrfz(idx_r,idx_r1,idx_tc)*odts ! RAIN2M
+ pnr_rfz(k) = MIN(DBLE(nr(k)*odts), pnr_rfz(k))
+ elseif (rr(k).gt. R1 .and. temp(k).lt.HGFR) then
+ pri_rfz(k) = rr(k)*odts
+ pnr_rfz(k) = nr(k)*odts ! RAIN2M
+ pni_rfz(k) = pnr_rfz(k)
+ endif
+ if (rc(k).gt. r_c(1)) then
+ pri_wfz(k) = tpi_qcfz(idx_c,idx_tc)*odts
+ pri_wfz(k) = MIN(DBLE(rc(k)*odts), pri_wfz(k))
+ pni_wfz(k) = tni_qcfz(idx_c,idx_tc)*odts
+ pni_wfz(k) = MIN(DBLE(Nt_c*odts), pri_wfz(k)/(2.*xm0i), &
+ pni_wfz(k))
+ endif
+
+!..Nucleate ice from deposition & condensation freezing (Cooper 1986)
+!.. but only if water sat and T<-12C or 25%+ ice supersaturated.
+ if ( (ssati(k).ge. 0.25) .or. (ssatw(k).gt. eps &
+ .and. temp(k).lt.261.15) ) then
+ xnc = MIN(250.E3, TNO*EXP(ATO*(T_0-temp(k))))
+ xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dtsave
+ pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts
+ pri_inu(k) = MIN(DBLE(rate_max), xm0i*pni_inu(k))
+ pni_inu(k) = pri_inu(k)/xm0i
+ endif
+
+!..Deposition/sublimation of cloud ice (Srivastava & Coen 1992).
+ if (L_qi(k)) then
+ lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
+ ilami = 1./lami
+ xDi = MAX(DBLE(D0i), (bm_i + mu_i + 1.) * ilami)
+ xmi = am_i*xDi**bm_i
+ oxmi = 1./xmi
+ pri_ide(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
+ *oig1*cig(5)*ni(k)*ilami
+
+ if (pri_ide(k) .lt. 0.0) then
+ pri_ide(k) = MAX(DBLE(-ri(k)*odts), pri_ide(k), DBLE(rate_max))
+ pni_ide(k) = pri_ide(k)*oxmi
+ pni_ide(k) = MAX(DBLE(-ni(k)*odts), pni_ide(k))
+ else
+ pri_ide(k) = MIN(pri_ide(k), DBLE(rate_max))
+ prs_ide(k) = (1.0D0-tpi_ide(idx_i,idx_i1))*pri_ide(k)
+ pri_ide(k) = tpi_ide(idx_i,idx_i1)*pri_ide(k)
+ endif
+
+!..Some cloud ice needs to move into the snow category. Use lookup
+!.. table that resulted from explicit bin representation of distrib.
+ if ( (idx_i.eq. ntb_i) .or. (xDi.gt. 5.0*D0s) ) then
+ prs_iau(k) = ri(k)*.99*odts
+ pni_iau(k) = ni(k)*.95*odts
+ elseif (xDi.lt. 0.1*D0s) then
+ prs_iau(k) = 0.
+ pni_iau(k) = 0.
+ else
+ prs_iau(k) = tps_iaus(idx_i,idx_i1)*odts
+ prs_iau(k) = MIN(DBLE(ri(k)*.99*odts), prs_iau(k))
+ pni_iau(k) = tni_iaus(idx_i,idx_i1)*odts
+ pni_iau(k) = MIN(DBLE(ni(k)*.95*odts), pni_iau(k))
+ endif
+ endif
+
+!..Deposition/sublimation of snow/graupel follows Srivastava & Coen
+!.. (1992).
+ if (L_qs(k)) then
+ C_snow = C_sqrd + (tempc+15.)*(C_cube-C_sqrd)/(-30.+15.)
+ C_snow = MAX(C_sqrd, MIN(C_snow, C_cube))
+ prs_sde(k) = C_snow*t1_subl*diffu(k)*ssati(k)*rvs &
+ * (t1_qs_sd*smo1(k) &
+ + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
+ if (prs_sde(k).lt. 0.) then
+ prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k), DBLE(rate_max))
+ else
+ prs_sde(k) = MIN(prs_sde(k), DBLE(rate_max))
+ endif
+ endif
+
+ if (L_qg(k) .and. ssati(k).lt. -eps) then
+ prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
+ * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
+ + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
+ if (prg_gde(k).lt. 0.) then
+ prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k), DBLE(rate_max))
+ else
+ prg_gde(k) = MIN(prg_gde(k), DBLE(rate_max))
+ endif
+ endif
+
+!..Snow collecting cloud ice. In CE, assume Di<<Ds and vti=~0.
+ if (L_qi(k)) then
+ lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
+ ilami = 1./lami
+ xDi = MAX(DBLE(D0i), (bm_i + mu_i + 1.) * ilami)
+ xmi = am_i*xDi**bm_i
+ oxmi = 1./xmi
+ if (rs(k).ge. r_s(1)) then
+ prs_sci(k) = t1_qs_qi*rhof(k)*Ef_si*ri(k)*smoe(k)
+ pni_sci(k) = prs_sci(k) * oxmi
+ endif
+
+!..Rain collecting cloud ice. In CE, assume Di<<Dr and vti=~0.
+ if (rr(k).ge. r_r(1) .and. mvd_r(k).gt. 4.*xDi) then
+ lamr = 1./ilamr(k)
+ pri_rci(k) = rhof(k)*t1_qr_qi*Ef_ri*ri(k)*N0_r(k) &
+ *((lamr+fv_r)**(-cre(9)))
+ pnr_rci(k) = rhof(k)*t1_qr_qi*Ef_ri*ni(k)*N0_r(k) & ! RAIN2M
+ *((lamr+fv_r)**(-cre(9)))
+ pni_rci(k) = pri_rci(k) * oxmi
+ prr_rci(k) = rhof(k)*t2_qr_qi*Ef_ri*ni(k)*N0_r(k) &
+ *((lamr+fv_r)**(-cre(8)))
+ prr_rci(k) = MIN(DBLE(rr(k)*odts), prr_rci(k))
+ prg_rci(k) = pri_rci(k) + prr_rci(k)
+ endif
+ endif
+
+!..Ice multiplication from rime-splinters (Hallet & Mossop 1974).
+ if (prg_gcw(k).gt. eps .and. tempc.gt.-8.0) then
+ tf = 0.
+ if (tempc.ge.-5.0 .and. tempc.lt.-3.0) then
+ tf = 0.5*(-3.0 - tempc)
+ elseif (tempc.gt.-8.0 .and. tempc.lt.-5.0) then
+ tf = 0.33333333*(8.0 + tempc)
+ endif
+ pni_ihm(k) = 3.5E8*tf*prg_gcw(k)
+ pri_ihm(k) = xm0i*pni_ihm(k)
+ prs_ihm(k) = prs_scw(k)/(prs_scw(k)+prg_gcw(k)) &
+ * pri_ihm(k)
+ prg_ihm(k) = prg_gcw(k)/(prs_scw(k)+prg_gcw(k)) &
+ * pri_ihm(k)
+ endif
+
+!..A portion of rimed snow converts to graupel but some remains snow.
+!.. Interp from 5 to 75% as riming factor increases from 5.0 to 30.0
+!.. 0.028 came from (.75-.05)/(30.-5.). This remains ad-hoc and should
+!.. be revisited.
+ if (prs_scw(k).gt.5.0*prs_sde(k) .and. &
+ prs_sde(k).gt.eps) then
+ r_frac = MIN(30.0D0, prs_scw(k)/prs_sde(k))
+ g_frac = MIN(0.75, 0.05 + (r_frac-5.)*.028)
+ vts_boost(k) = MIN(1.5, 1.1 + (r_frac-5.)*.016)
+ prg_scw(k) = g_frac*prs_scw(k)
+ prs_scw(k) = (1. - g_frac)*prs_scw(k)
+ endif
+
+ else
+
+!..Melt snow and graupel and enhance from collisions with liquid.
+!.. We also need to sublimate snow and graupel if subsaturated.
+ if (L_qs(k)) then
+ prr_sml(k) = tempc*tcond(k)*(t1_qs_me*smo1(k) &
+ + t2_qs_me*rhof2(k)*vsc2(k)*smof(k))
+ prr_sml(k) = prr_sml(k) + 4218.*olfus*tempc &
+ * (prr_rcs(k)+prs_scw(k))
+ prr_sml(k) = MIN(DBLE(rs(k)*odts), prr_sml(k))
+ pnr_sml(k) = smo0(k)/rs(k)*prr_sml(k) * 10.0**(-0.50*tempc) ! RAIN2M
+ pnr_sml(k) = MIN(DBLE(smo0(k)*odts), pnr_sml(k))
+ if (tempc.gt.3.5 .or. rs(k).lt.0.005E-3) pnr_sml(k)=0.0
+
+ if (ssati(k).lt. 0.) then
+ prs_sde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
+ * (t1_qs_sd*smo1(k) &
+ + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
+ prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k))
+ endif
+ endif
+
+ if (L_qg(k)) then
+ prr_gml(k) = tempc*N0_g(k)*tcond(k) &
+ *(t1_qg_me*ilamg(k)**cge(10) &
+ + t2_qg_me*rhof2(k)*vsc2(k)*ilamg(k)**cge(11))
+ prr_gml(k) = prr_gml(k) + 4218.*olfus*tempc &
+ * (prr_rcg(k)+prg_gcw(k))
+ prr_gml(k) = MIN(DBLE(rg(k)*odts), prr_gml(k))
+ pnr_gml(k) = (N0_g(k) / (cgg(1)*am_g*N0_g(k)/rg(k))**oge1) & ! RAIN2M
+ / rg(k) * prr_gml(k) * 10.0**(-0.35*tempc)
+ if (rg(k).lt.0.005E-3) pnr_gml(k)=0.0
+
+ if (ssati(k).lt. 0.) then
+ prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
+ * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
+ + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
+ prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k))
+ endif
+ endif
+
+ endif
+
+ enddo
+ endif
+
+!+---+-----------------------------------------------------------------+
+!..Ensure we do not deplete more hydrometeor species than exists.
+!+---+-----------------------------------------------------------------+
+ do k = kts, kte
+
+!..If ice supersaturated, ensure sum of depos growth terms does not
+!.. deplete more vapor than possibly exists. If subsaturated, limit
+!.. sum of sublimation terms such that vapor does not reproduce ice
+!.. supersat again.
+ sump = pri_inu(k) + pri_ide(k) + prs_ide(k) &
+ + prs_sde(k) + prg_gde(k)
+ rate_max = (qv(k)-qvsi(k))*odts*0.999
+ if ( (sump.gt. eps .and. sump.gt. rate_max) .or. &
+ (sump.lt. -eps .and. sump.lt. rate_max) ) then
+ ratio = rate_max/sump
+ pri_inu(k) = pri_inu(k) * ratio
+ pri_ide(k) = pri_ide(k) * ratio
+ pni_ide(k) = pni_ide(k) * ratio
+ prs_ide(k) = prs_ide(k) * ratio
+ prs_sde(k) = prs_sde(k) * ratio
+ prg_gde(k) = prg_gde(k) * ratio
+ endif
+
+!..Cloud water conservation.
+ sump = -prr_wau(k) - pri_wfz(k) - prr_rcw(k) &
+ - prs_scw(k) - prg_scw(k) - prg_gcw(k)
+ rate_max = -rc(k)*odts
+ if (sump.lt. rate_max .and. L_qc(k)) then
+ ratio = rate_max/sump
+ prr_wau(k) = prr_wau(k) * ratio
+ pri_wfz(k) = pri_wfz(k) * ratio
+ prr_rcw(k) = prr_rcw(k) * ratio
+ prs_scw(k) = prs_scw(k) * ratio
+ prg_scw(k) = prg_scw(k) * ratio
+ prg_gcw(k) = prg_gcw(k) * ratio
+ endif
+
+!..Cloud ice conservation.
+ sump = pri_ide(k) - prs_iau(k) - prs_sci(k) &
+ - pri_rci(k)
+ rate_max = -ri(k)*odts
+ if (sump.lt. rate_max .and. L_qi(k)) then
+ ratio = rate_max/sump
+ pri_ide(k) = pri_ide(k) * ratio
+ prs_iau(k) = prs_iau(k) * ratio
+ prs_sci(k) = prs_sci(k) * ratio
+ pri_rci(k) = pri_rci(k) * ratio
+ endif
+
+!..Rain conservation.
+ sump = -prg_rfz(k) - pri_rfz(k) - prr_rci(k) &
+ + prr_rcs(k) + prr_rcg(k)
+ rate_max = -rr(k)*odts
+ if (sump.lt. rate_max .and. L_qr(k)) then
+ ratio = rate_max/sump
+ prg_rfz(k) = prg_rfz(k) * ratio
+ pri_rfz(k) = pri_rfz(k) * ratio
+ prr_rci(k) = prr_rci(k) * ratio
+ prr_rcs(k) = prr_rcs(k) * ratio
+ prr_rcg(k) = prr_rcg(k) * ratio
+ endif
+
+!..Snow conservation.
+ sump = prs_sde(k) - prs_ihm(k) - prr_sml(k) &
+ + prs_rcs(k)
+ rate_max = -rs(k)*odts
+ if (sump.lt. rate_max .and. L_qs(k)) then
+ ratio = rate_max/sump
+ prs_sde(k) = prs_sde(k) * ratio
+ prs_ihm(k) = prs_ihm(k) * ratio
+ prr_sml(k) = prr_sml(k) * ratio
+ prs_rcs(k) = prs_rcs(k) * ratio
+ endif
+
+!..Graupel conservation.
+ sump = prg_gde(k) - prg_ihm(k) - prr_gml(k) &
+ + prg_rcg(k)
+ rate_max = -rg(k)*odts
+ if (sump.lt. rate_max .and. L_qg(k)) then
+ ratio = rate_max/sump
+ prg_gde(k) = prg_gde(k) * ratio
+ prg_ihm(k) = prg_ihm(k) * ratio
+ prr_gml(k) = prr_gml(k) * ratio
+ prg_rcg(k) = prg_rcg(k) * ratio
+ endif
+
+!..Re-enforce proper mass conservation for subsequent elements in case
+!.. any of the above terms were altered. Thanks P. Blossey. 2009Sep28
+ pri_ihm(k) = prs_ihm(k) + prg_ihm(k)
+ ratio = MIN( ABS(prr_rcg(k)), ABS(prg_rcg(k)) )
+ prr_rcg(k) = ratio * SIGN(1.0, SNGL(prr_rcg(k)))
+ prg_rcg(k) = -prr_rcg(k)
+ if (temp(k).lt.T_0) then
+ prg_rcs(k) = prs_rcs(k) + prr_rcs(k)
+ else
+ ratio = MIN( ABS(prr_rcs(k)), ABS(prs_rcs(k)) )
+ prr_rcs(k) = ratio * SIGN(1.0, SNGL(prr_rcs(k)))
+ prs_rcs(k) = -prr_rcs(k)
+ endif
+
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..Calculate tendencies of all species but constrain the number of ice
+!.. to reasonable values.
+!+---+-----------------------------------------------------------------+
+ do k = kts, kte
+ orho = 1./rho(k)
+ lfus2 = lsub - lvap(k)
+
+!..Water vapor tendency
+ qvten(k) = qvten(k) + (-pri_inu(k) - pri_ide(k) &
+ - prs_ide(k) - prs_sde(k) - prg_gde(k)) &
+ * orho
+
+!..Cloud water tendency
+ qcten(k) = qcten(k) + (-prr_wau(k) - pri_wfz(k) &
+ - prr_rcw(k) - prs_scw(k) - prg_scw(k) &
+ - prg_gcw(k)) &
+ * orho
+
+!..Cloud ice mixing ratio tendency
+ qiten(k) = qiten(k) + (pri_inu(k) + pri_ihm(k) &
+ + pri_wfz(k) + pri_rfz(k) + pri_ide(k) &
+ - prs_iau(k) - prs_sci(k) - pri_rci(k)) &
+ * orho
+
+!..Cloud ice number tendency.
+ niten(k) = niten(k) + (pni_inu(k) + pni_ihm(k) &
+ + pni_wfz(k) + pni_rfz(k) + pni_ide(k) &
+ - pni_iau(k) - pni_sci(k) - pni_rci(k)) &
+ * orho
+
+!..Cloud ice mass/number balance; keep mass-wt mean size between
+!.. 20 and 300 microns. Also no more than 500 xtals per liter.
+ xri=MAX(R1,(qi1d(k) + qiten(k)*dtsave)*rho(k))
+ xni=MAX(1.,(ni1d(k) + niten(k)*dtsave)*rho(k))
+ if (xri.gt. R1) then
+ lami = (am_i*cig(2)*oig1*xni/xri)**obmi
+ ilami = 1./lami
+ xDi = (bm_i + mu_i + 1.) * ilami
+ if (xDi.lt. 20.E-6) then
+ lami = cie(2)/20.E-6
+ xni = MIN(500.D3, cig(1)*oig2*xri/am_i*lami**bm_i)
+ niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
+ elseif (xDi.gt. 300.E-6) then
+ lami = cie(2)/300.E-6
+ xni = cig(1)*oig2*xri/am_i*lami**bm_i
+ niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
+ endif
+ else
+ niten(k) = -ni1d(k)*odts
+ endif
+ xni=MAX(0.,(ni1d(k) + niten(k)*dtsave)*rho(k))
+ if (xni.gt.500.E3) &
+ niten(k) = (500.E3-ni1d(k)*rho(k))*odts*orho
+
+!..Rain tendency
+ qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) &
+ + prr_sml(k) + prr_gml(k) + prr_rcs(k) &
+ + prr_rcg(k) - prg_rfz(k) &
+ - pri_rfz(k) - prr_rci(k)) &
+ * orho
+
+!..Rain number tendency
+ nrten(k) = nrten(k) + (pnr_wau(k) + pnr_sml(k) + pnr_gml(k) &
+ - (pnr_rfz(k) + pnr_rcr(k) + pnr_rcg(k) &
+ + pnr_rcs(k) + pnr_rci(k)) ) &
+ * orho
+
+!..Rain mass/number balance; keep median volume diameter between
+!.. 37 microns (D0r*0.75) and 2.5 mm.
+ xrr=MAX(R1,(qr1d(k) + qrten(k)*dtsave)*rho(k))
+ xnr=MAX(1.,(nr1d(k) + nrten(k)*dtsave)*rho(k))
+ if (xrr.gt. R1) then
+ lamr = (am_r*crg(3)*org2*xnr/xrr)**obmr
+ mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
+ if (mvd_r(k) .gt. 2.5E-3) then
+ mvd_r(k) = 2.5E-3
+ lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
+ xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
+ nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
+ elseif (mvd_r(k) .lt. D0r*0.75) then
+ mvd_r(k) = D0r*0.75
+ lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
+ xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
+ nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
+ endif
+ else
+ qrten(k) = -qr1d(k)*odts
+ nrten(k) = -nr1d(k)*odts
+ endif
+
+!..Snow tendency
+ qsten(k) = qsten(k) + (prs_iau(k) + prs_sde(k) &
+ + prs_sci(k) + prs_scw(k) + prs_rcs(k) &
+ + prs_ide(k) - prs_ihm(k) - prr_sml(k)) &
+ * orho
+
+!..Graupel tendency
+ qgten(k) = qgten(k) + (prg_scw(k) + prg_rfz(k) &
+ + prg_gde(k) + prg_rcg(k) + prg_gcw(k) &
+ + prg_rci(k) + prg_rcs(k) - prg_ihm(k) &
+ - prr_gml(k)) &
+ * orho
+
+!..Temperature tendency
+ if (temp(k).lt.T_0) then
+ tten(k) = tten(k) &
+ + ( lsub*ocp(k)*(pri_inu(k) + pri_ide(k) &
+ + prs_ide(k) + prs_sde(k) &
+ + prg_gde(k)) &
+ + lfus2*ocp(k)*(pri_wfz(k) + pri_rfz(k) &
+ + prg_rfz(k) + prs_scw(k) &
+ + prg_scw(k) + prg_gcw(k) &
+ + prg_rcs(k) + prs_rcs(k) &
+ + prr_rci(k) + prg_rcg(k)) &
+ )*orho * (1-IFDRY)
+ else
+ tten(k) = tten(k) &
+ + ( lfus*ocp(k)*(-prr_sml(k) - prr_gml(k) &
+ - prr_rcg(k) - prr_rcs(k)) &
+ + lsub*ocp(k)*(prs_sde(k) + prg_gde(k)) &
+ )*orho * (1-IFDRY)
+ endif
+
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..Update variables for TAU+1 before condensation & sedimention.
+!+---+-----------------------------------------------------------------+
+ do k = kts, kte
+ temp(k) = t1d(k) + DT*tten(k)
+ otemp = 1./temp(k)
+ tempc = temp(k) - 273.15
+ qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
+ rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
+ rhof(k) = SQRT(RHO_NOT/rho(k))
+ rhof2(k) = SQRT(rhof(k))
+ qvs(k) = rslf(pres(k), temp(k))
+ ssatw(k) = qv(k)/qvs(k) - 1.
+ if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
+ diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
+ if (tempc .ge. 0.0) then
+ visco(k) = (1.718+0.0049*tempc)*1.0E-5
+ else
+ visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
+ endif
+ vsc2(k) = SQRT(rho(k)/visco(k))
+ lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
+ tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
+ ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
+ lvt2(k)=lvap(k)*lvap(k)*ocp(k)*oRv*otemp*otemp
+
+ if ((qc1d(k) + qcten(k)*DT) .gt. R1) then
+ rc(k) = (qc1d(k) + qcten(k)*DT)*rho(k)
+ L_qc(k) = .true.
+ else
+ rc(k) = R1
+ L_qc(k) = .false.
+ endif
+
+ if ((qi1d(k) + qiten(k)*DT) .gt. R1) then
+ ri(k) = (qi1d(k) + qiten(k)*DT)*rho(k)
+ ni(k) = MAX(1., (ni1d(k) + niten(k)*DT)*rho(k))
+ L_qi(k) = .true.
+ else
+ ri(k) = R1
+ ni(k) = 1.
+ L_qi(k) = .false.
+ endif
+
+ if ((qr1d(k) + qrten(k)*DT) .gt. R1) then
+ rr(k) = (qr1d(k) + qrten(k)*DT)*rho(k)
+ nr(k) = MAX(1., (nr1d(k) + nrten(k)*DT)*rho(k))
+ L_qr(k) = .true.
+ lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
+ mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
+ if (mvd_r(k) .gt. 2.5E-3) then
+ mvd_r(k) = 2.5E-3
+ lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
+ nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
+ elseif (mvd_r(k) .lt. D0r*0.75) then
+ mvd_r(k) = D0r*0.75
+ lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
+ nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
+ endif
+ else
+ rr(k) = R1
+ nr(k) = 1.
+ L_qr(k) = .false.
+ endif
+
+ if ((qs1d(k) + qsten(k)*DT) .gt. R1) then
+ rs(k) = (qs1d(k) + qsten(k)*DT)*rho(k)
+ L_qs(k) = .true.
+ else
+ rs(k) = R1
+ L_qs(k) = .false.
+ endif
+
+ if ((qg1d(k) + qgten(k)*DT) .gt. R1) then
+ rg(k) = (qg1d(k) + qgten(k)*DT)*rho(k)
+ L_qg(k) = .true.
+ else
+ rg(k) = R1
+ L_qg(k) = .false.
+ endif
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..With tendency-updated mixing ratios, recalculate snow moments and
+!.. intercepts/slopes of graupel and rain.
+!+---+-----------------------------------------------------------------+
+ if (.not. iiwarm) then
+ do k = kts, kte
+ if (.not. L_qs(k)) CYCLE
+ tc0 = MIN(-0.1, temp(k)-273.15)
+ smob(k) = rs(k)*oams
+
+!..All other moments based on reference, 2nd moment. If bm_s.ne.2,
+!.. then we must compute actual 2nd moment and use as reference.
+ if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
+ smo2(k) = smob(k)
+ else
+ loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
+ + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
+ + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
+ + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
+ + sa(10)*bm_s*bm_s*bm_s
+ a_ = 10.0**loga_
+ b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
+ + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
+ + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
+ + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
+ + sb(10)*bm_s*bm_s*bm_s
+ smo2(k) = (smob(k)/a_)**(1./b_)
+ endif
+
+!..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
+ loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
+ + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
+ + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
+ + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
+ + sa(10)*cse(1)*cse(1)*cse(1)
+ a_ = 10.0**loga_
+ b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
+ + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
+ + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
+ + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
+ smoc(k) = a_ * smo2(k)**b_
+
+!..Calculate bm_s+bv_s (th) moment. Useful for sedimentation.
+ loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(14) &
+ + sa(4)*tc0*cse(14) + sa(5)*tc0*tc0 &
+ + sa(6)*cse(14)*cse(14) + sa(7)*tc0*tc0*cse(14) &
+ + sa(8)*tc0*cse(14)*cse(14) + sa(9)*tc0*tc0*tc0 &
+ + sa(10)*cse(14)*cse(14)*cse(14)
+ a_ = 10.0**loga_
+ b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(14) + sb(4)*tc0*cse(14) &
+ + sb(5)*tc0*tc0 + sb(6)*cse(14)*cse(14) &
+ + sb(7)*tc0*tc0*cse(14) + sb(8)*tc0*cse(14)*cse(14) &
+ + sb(9)*tc0*tc0*tc0 + sb(10)*cse(14)*cse(14)*cse(14)
+ smod(k) = a_ * smo2(k)**b_
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..Calculate y-intercept, slope values for graupel.
+!+---+-----------------------------------------------------------------+
+ do k = kte, kts, -1
+ N0_exp = (gonv_max-gonv_min)*0.5D0 &
+ * tanh((0.01E-3-(rc(k)+rr(k)))/0.75E-3) &
+ + (gonv_max+gonv_min)*0.5D0
+ lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
+ lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
+ ilamg(k) = 1./lamg
+ N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2)
+ enddo
+
+ endif
+
+!+---+-----------------------------------------------------------------+
+!..Calculate y-intercept, slope values for rain.
+!+---+-----------------------------------------------------------------+
+ do k = kte, kts, -1
+ lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
+ ilamr(k) = 1./lamr
+ mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
+ N0_r(k) = nr(k)*org2*lamr**cre(2)
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..Cloud water condensation and evaporation. Newly formulated using
+!.. Newton-Raphson iterations (3 should suffice) as provided by B. Hall.
+!+---+-----------------------------------------------------------------+
+ do k = kts, kte
+ if ( (ssatw(k).gt. eps) .or. (ssatw(k).lt. -eps .and. &
+ L_qc(k)) ) then
+ clap = (qv(k)-qvs(k))/(1. + lvt2(k)*qvs(k))
+ do n = 1, 3
+ fcd = qvs(k)* EXP(lvt2(k)*clap) - qv(k) + clap
+ dfcd = qvs(k)*lvt2(k)* EXP(lvt2(k)*clap) + 1.
+ clap = clap - fcd/dfcd
+ enddo
+ xrc = rc(k) + clap
+ if (xrc.gt. 0.0) then
+ prw_vcd(k) = clap*odt
+ else
+ prw_vcd(k) = -rc(k)/rho(k)*odt
+ endif
+
+ qcten(k) = qcten(k) + prw_vcd(k)
+ qvten(k) = qvten(k) - prw_vcd(k)
+ tten(k) = tten(k) + lvap(k)*ocp(k)*prw_vcd(k)*(1-IFDRY)
+ rc(k) = MAX(R1, (qc1d(k) + DT*qcten(k))*rho(k))
+ qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
+ temp(k) = t1d(k) + DT*tten(k)
+ rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
+ qvs(k) = rslf(pres(k), temp(k))
+ ssatw(k) = qv(k)/qvs(k) - 1.
+ endif
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!.. If still subsaturated, allow rain to evaporate, following
+!.. Srivastava & Coen (1992).
+!+---+-----------------------------------------------------------------+
+ do k = kts, kte
+ if ( (ssatw(k).lt. -eps) .and. L_qr(k) &
+ .and. (.not.(prw_vcd(k).gt. 0.)) ) then
+ tempc = temp(k) - 273.15
+ otemp = 1./temp(k)
+ rhof(k) = SQRT(RHO_NOT/rho(k))
+ rhof2(k) = SQRT(rhof(k))
+ diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
+ if (tempc .ge. 0.0) then
+ visco(k) = (1.718+0.0049*tempc)*1.0E-5
+ else
+ visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
+ endif
+ vsc2(k) = SQRT(rho(k)/visco(k))
+ lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
+ tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
+ ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
+
+ rvs = rho(k)*qvs(k)
+ rvs_p = rvs*otemp*(lvap(k)*otemp*oRv - 1.)
+ rvs_pp = rvs * ( otemp*(lvap(k)*otemp*oRv - 1.) &
+ *otemp*(lvap(k)*otemp*oRv - 1.) &
+ + (-2.*lvap(k)*otemp*otemp*otemp*oRv) &
+ + otemp*otemp)
+ gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
+ alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
+ * rvs_pp/rvs_p * rvs/rvs_p
+ alphsc = MAX(1.E-9, alphsc)
+ xsat = ssatw(k)
+ if (xsat.lt. -1.E-9) xsat = -1.E-9
+ t1_evap = 2.*PI*( 1.0 - alphsc*xsat &
+ + 2.*alphsc*alphsc*xsat*xsat &
+ - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
+ / (1.+gamsc)
+
+ lamr = 1./ilamr(k)
+ prv_rev(k) = t1_evap*diffu(k)*(-ssatw(k))*N0_r(k)*rvs &
+ * (t1_qr_ev*ilamr(k)**cre(10) &
+ + t2_qr_ev*vsc2(k)*rhof2(k)*((lamr+0.5*fv_r)**(-cre(11))))
+ prv_rev(k) = MIN(DBLE(rr(k)/rho(k)*odts), &
+ prv_rev(k)/rho(k))
+ pnr_rev(k) = MIN(DBLE(nr(k)*0.99/rho(k)*odts), & ! RAIN2M
+ prv_rev(k) * nr(k)/rr(k))
+
+ qrten(k) = qrten(k) - prv_rev(k)
+ qvten(k) = qvten(k) + prv_rev(k)
+ nrten(k) = nrten(k) - pnr_rev(k)
+ tten(k) = tten(k) - lvap(k)*ocp(k)*prv_rev(k)*(1-IFDRY)
+
+ rr(k) = MAX(R1, (qr1d(k) + DT*qrten(k))*rho(k))
+ qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
+ nr(k) = MAX(1.,(nr1d(k) + DT*nrten(k))*rho(k))
+ temp(k) = t1d(k) + DT*tten(k)
+ rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
+ endif
+
+!+---+-----------------------------------------------------------------+
+! if( debug_flag .and. k.lt.42) then
+! if (k.eq.1) write(mp_debug,*) 'DEBUG-GT: prg_scw, prg_rfz, prg_gde, prg_rcg, prg_gcw, prg_rci, prg_rcs, prg_ihm, prr_gml, rg, N0_g, ilamg'
+! if (k.eq.1) CALL wrf_debug(0, mp_debug)
+! write(mp_debug, 'a, i2, 1x, 12(1x,e13.6,1x)') ' GT,k= ', k, &
+! prg_scw(k), prg_rfz(k), prg_gde(k), prg_rcg(k), prg_gcw(k), &
+! prg_rci(k), prg_rcs(k), prg_ihm(k), prr_gml(k), &
+! rg(k), N0_g(k), ilamg(k)
+! CALL wrf_debug(0, mp_debug)
+! endif
+!+---+-----------------------------------------------------------------+
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!..Find max terminal fallspeed (distribution mass-weighted mean
+!.. velocity) and use it to determine if we need to split the timestep
+!.. (var nstep>1). Either way, only bother to do sedimentation below
+!.. 1st level that contains any sedimenting particles (k=ksed1 on down).
+!.. New in v3.0+ is computing separate for rain, ice, snow, and
+!.. graupel species thus making code faster with credit to J. Schmidt.
+!+---+-----------------------------------------------------------------+
+ nstep = 0
+ onstep(:) = 1.0
+ ksed1(:) = 1
+ do k = kte+1, kts, -1
+ vtrk(k) = 0.
+ vtnrk(k) = 0.
+ vtik(k) = 0.
+ vtnik(k) = 0.
+ vtsk(k) = 0.
+ vtgk(k) = 0.
+ enddo
+ do k = kte, kts, -1
+ vtr = 0.
+ rhof(k) = SQRT(RHO_NOT/rho(k))
+
+ if (rr(k).gt. R2) then
+ lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
+ vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) &
+ *((lamr+fv_r)**(-cre(6)))
+ vtrk(k) = vtr
+! First below is technically correct:
+! vtr = rhof(k)*av_r*crg(5)*org2 * lamr**cre(2) &
+! *((lamr+fv_r)**(-cre(5)))
+! Test: make number fall faster (but still slower than mass)
+! Goal: less prominent size sorting
+ vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) &
+ *((lamr+fv_r)**(-cre(7)))
+ vtnrk(k) = vtr
+ endif
+
+ if (MAX(vtrk(k),vtnrk(k)) .gt. 1.E-3) then
+ ksed1(1) = MAX(ksed1(1), k)
+ delta_tp = dzq(k)/(MAX(vtrk(k),vtnrk(k)))
+ nstep = MAX(nstep, INT(DT/delta_tp + 1.))
+ endif
+ enddo
+ if (ksed1(1) .eq. kte) ksed1(1) = kte-1
+ if (nstep .gt. 0) onstep(1) = 1./REAL(nstep)
+
+!+---+-----------------------------------------------------------------+
+
+ if (.not. iiwarm) then
+
+ nstep = 0
+ do k = kte, kts, -1
+ vti = 0.
+
+ if (ri(k).gt. R2) then
+ lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
+ ilami = 1./lami
+ vti = rhof(k)*av_i*cig(3)*oig2 * ilami**bv_i
+ vtik(k) = vti
+ vti = rhof(k)*av_i*cig(4)*oig1 * ilami**bv_i
+ vtnik(k) = vti
+ endif
+
+ if (vtik(k) .gt. 1.E-3) then
+ ksed1(2) = MAX(ksed1(2), k)
+ delta_tp = dzq(k)/vtik(k)
+ nstep = MAX(nstep, INT(DT/delta_tp + 1.))
+ endif
+ enddo
+ if (ksed1(2) .eq. kte) ksed1(2) = kte-1
+ if (nstep .gt. 0) onstep(2) = 1./REAL(nstep)
+
+!+---+-----------------------------------------------------------------+
+
+ nstep = 0
+ do k = kte, kts, -1
+ vts = 0.
+
+ if (rs(k).gt. R2) then
+ xDs = smoc(k) / smob(k)
+ Mrat = 1./xDs
+ ils1 = 1./(Mrat*Lam0 + fv_s)
+ ils2 = 1./(Mrat*Lam1 + fv_s)
+ t1_vts = Kap0*csg(4)*ils1**cse(4)
+ t2_vts = Kap1*Mrat**mu_s*csg(10)*ils2**cse(10)
+ ils1 = 1./(Mrat*Lam0)
+ ils2 = 1./(Mrat*Lam1)
+ t3_vts = Kap0*csg(1)*ils1**cse(1)
+ t4_vts = Kap1*Mrat**mu_s*csg(7)*ils2**cse(7)
+ vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts)
+ if (temp(k).gt. T_0) then
+ vtsk(k) = MAX(vts*vts_boost(k), vtrk(k))
+ else
+ vtsk(k) = vts*vts_boost(k)
+ endif
+ endif
+
+ if (vtsk(k) .gt. 1.E-3) then
+ ksed1(3) = MAX(ksed1(3), k)
+ delta_tp = dzq(k)/vtsk(k)
+ nstep = MAX(nstep, INT(DT/delta_tp + 1.))
+ endif
+ enddo
+ if (ksed1(3) .eq. kte) ksed1(3) = kte-1
+ if (nstep .gt. 0) onstep(3) = 1./REAL(nstep)
+
+!+---+-----------------------------------------------------------------+
+
+ nstep = 0
+ do k = kte, kts, -1
+ vtg = 0.
+
+ if (rg(k).gt. R2) then
+ vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
+ if (temp(k).gt. T_0) then
+ vtgk(k) = MAX(vtg, vtrk(k))
+ else
+ vtgk(k) = vtg
+ endif
+ endif
+
+ if (vtgk(k) .gt. 1.E-3) then
+ ksed1(4) = MAX(ksed1(4), k)
+ delta_tp = dzq(k)/vtgk(k)
+ nstep = MAX(nstep, INT(DT/delta_tp + 1.))
+ endif
+ enddo
+ if (ksed1(4) .eq. kte) ksed1(4) = kte-1
+ if (nstep .gt. 0) onstep(4) = 1./REAL(nstep)
+ endif
+
+!+---+-----------------------------------------------------------------+
+!..Sedimentation of mixing ratio is the integral of v(D)*m(D)*N(D)*dD,
+!.. whereas neglect m(D) term for number concentration. Therefore,
+!.. cloud ice has proper differential sedimentation.
+!.. New in v3.0+ is computing separate for rain, ice, snow, and
+!.. graupel species thus making code faster with credit to J. Schmidt.
+!+---+-----------------------------------------------------------------+
+
+ nstep = NINT(1./onstep(1))
+ do n = 1, nstep
+ do k = kte, kts, -1
+ sed_r(k) = vtrk(k)*rr(k)
+ sed_n(k) = vtnrk(k)*nr(k)
+ enddo
+ k = kte
+ odzq = 1./dzq(k)
+ orho = 1./rho(k)
+ qrten(k) = qrten(k) - sed_r(k)*odzq*onstep(1)*orho
+ nrten(k) = nrten(k) - sed_n(k)*odzq*onstep(1)*orho
+ rr(k) = MAX(R1, rr(k) - sed_r(k)*odzq*DT*onstep(1))
+ nr(k) = MAX(1., nr(k) - sed_n(k)*odzq*DT*onstep(1))
+ do k = ksed1(1), kts, -1
+ odzq = 1./dzq(k)
+ orho = 1./rho(k)
+ qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) &
+ *odzq*onstep(1)*orho
+ nrten(k) = nrten(k) + (sed_n(k+1)-sed_n(k)) &
+ *odzq*onstep(1)*orho
+ rr(k) = MAX(R1, rr(k) + (sed_r(k+1)-sed_r(k)) &
+ *odzq*DT*onstep(1))
+ nr(k) = MAX(1., nr(k) + (sed_n(k+1)-sed_n(k)) &
+ *odzq*DT*onstep(1))
+ enddo
+
+ pptrain = pptrain + sed_r(kts)*DT*onstep(1)
+ enddo
+
+!+---+-----------------------------------------------------------------+
+
+ nstep = NINT(1./onstep(2))
+ do n = 1, nstep
+ do k = kte, kts, -1
+ sed_i(k) = vtik(k)*ri(k)
+ sed_n(k) = vtnik(k)*ni(k)
+ enddo
+ k = kte
+ odzq = 1./dzq(k)
+ orho = 1./rho(k)
+ qiten(k) = qiten(k) - sed_i(k)*odzq*onstep(2)*orho
+ niten(k) = niten(k) - sed_n(k)*odzq*onstep(2)*orho
+ ri(k) = MAX(R1, ri(k) - sed_i(k)*odzq*DT*onstep(2))
+ ni(k) = MAX(1., ni(k) - sed_n(k)*odzq*DT*onstep(2))
+ do k = ksed1(2), kts, -1
+ odzq = 1./dzq(k)
+ orho = 1./rho(k)
+ qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) &
+ *odzq*onstep(2)*orho
+ niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) &
+ *odzq*onstep(2)*orho
+ ri(k) = MAX(R1, ri(k) + (sed_i(k+1)-sed_i(k)) &
+ *odzq*DT*onstep(2))
+ ni(k) = MAX(1., ni(k) + (sed_n(k+1)-sed_n(k)) &
+ *odzq*DT*onstep(2))
+ enddo
+
+ pptice = pptice + sed_i(kts)*DT*onstep(2)
+ enddo
+
+!+---+-----------------------------------------------------------------+
+
+ nstep = NINT(1./onstep(3))
+ do n = 1, nstep
+ do k = kte, kts, -1
+ sed_s(k) = vtsk(k)*rs(k)
+ enddo
+ k = kte
+ odzq = 1./dzq(k)
+ orho = 1./rho(k)
+ qsten(k) = qsten(k) - sed_s(k)*odzq*onstep(3)*orho
+ rs(k) = MAX(R1, rs(k) - sed_s(k)*odzq*DT*onstep(3))
+ do k = ksed1(3), kts, -1
+ odzq = 1./dzq(k)
+ orho = 1./rho(k)
+ qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) &
+ *odzq*onstep(3)*orho
+ rs(k) = MAX(R1, rs(k) + (sed_s(k+1)-sed_s(k)) &
+ *odzq*DT*onstep(3))
+ enddo
+
+ pptsnow = pptsnow + sed_s(kts)*DT*onstep(3)
+ enddo
+
+!+---+-----------------------------------------------------------------+
+
+ nstep = NINT(1./onstep(4))
+ do n = 1, nstep
+ do k = kte, kts, -1
+ sed_g(k) = vtgk(k)*rg(k)
+ enddo
+ k = kte
+ odzq = 1./dzq(k)
+ orho = 1./rho(k)
+ qgten(k) = qgten(k) - sed_g(k)*odzq*onstep(4)*orho
+ rg(k) = MAX(R1, rg(k) - sed_g(k)*odzq*DT*onstep(4))
+ do k = ksed1(4), kts, -1
+ odzq = 1./dzq(k)
+ orho = 1./rho(k)
+ qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) &
+ *odzq*onstep(4)*orho
+ rg(k) = MAX(R1, rg(k) + (sed_g(k+1)-sed_g(k)) &
+ *odzq*DT*onstep(4))
+ enddo
+
+ pptgraul = pptgraul + sed_g(kts)*DT*onstep(4)
+ enddo
+
+!+---+-----------------------------------------------------------------+
+!.. Instantly melt any cloud ice into cloud water if above 0C and
+!.. instantly freeze any cloud water found below HGFR.
+!+---+-----------------------------------------------------------------+
+ if (.not. iiwarm) then
+ do k = kts, kte
+ xri = MAX(0.0, qi1d(k) + qiten(k)*DT)
+ if ( (temp(k).gt. T_0) .and. (xri.gt. 0.0) ) then
+ qcten(k) = qcten(k) + xri*odt
+ qiten(k) = qiten(k) - xri*odt
+ niten(k) = -ni1d(k)*odt
+ tten(k) = tten(k) - lfus*ocp(k)*xri*odt*(1-IFDRY)
+ endif
+
+ xrc = MAX(0.0, qc1d(k) + qcten(k)*DT)
+ if ( (temp(k).lt. HGFR) .and. (xrc.gt. 0.0) ) then
+ lfus2 = lsub - lvap(k)
+ qiten(k) = qiten(k) + xrc*odt
+ niten(k) = niten(k) + xrc/xm0i * odt
+ qcten(k) = qcten(k) - xrc*odt
+ tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-IFDRY)
+ endif
+ enddo
+ endif
+
+!+---+-----------------------------------------------------------------+
+!.. All tendencies computed, apply and pass back final values to parent.
+!+---+-----------------------------------------------------------------+
+ do k = kts, kte
+ t1d(k) = t1d(k) + tten(k)*DT
+ qv1d(k) = MAX(1.E-10, qv1d(k) + qvten(k)*DT)
+ qc1d(k) = qc1d(k) + qcten(k)*DT
+ if (qc1d(k) .le. R1) qc1d(k) = 0.0
+ qi1d(k) = qi1d(k) + qiten(k)*DT
+ ni1d(k) = ni1d(k) + niten(k)*DT
+ if (qi1d(k) .le. R1) then
+ qi1d(k) = 0.0
+ ni1d(k) = 0.0
+ else
+ if (ni1d(k) .gt. 1.0) then
+ lami = (am_i*cig(2)*oig1*ni1d(k)/qi1d(k))**obmi
+ ilami = 1./lami
+ xDi = (bm_i + mu_i + 1.) * ilami
+ if (xDi.lt. 20.E-6) then
+ lami = cie(2)/20.E-6
+ elseif (xDi.gt. 300.E-6) then
+ lami = cie(2)/300.E-6
+ endif
+ else
+ lami = cie(2)/D0s
+ endif
+ ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, &
+ 500.D3/rho(k))
+ endif
+ qr1d(k) = qr1d(k) + qrten(k)*DT
+ nr1d(k) = nr1d(k) + nrten(k)*DT
+ if (qr1d(k) .le. R1) then
+ qr1d(k) = 0.0
+ nr1d(k) = 0.0
+ else
+ if (nr1d(k) .gt. 1.0) then
+ lamr = (am_r*crg(3)*org2*nr1d(k)/qr1d(k))**obmr
+ mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
+ if (mvd_r(k) .gt. 2.5E-3) then
+ mvd_r(k) = 2.5E-3
+ elseif (mvd_r(k) .lt. D0r*0.75) then
+ mvd_r(k) = D0r*0.75
+ endif
+ else
+ if (qr1d(k) .gt. R2) then
+ mvd_r(k) = 2.5E-3
+ else
+ mvd_r(k) = 2.5E-3 / 3.0**(ALOG10(R2)-ALOG10(qr1d(k)))
+ endif
+ endif
+ lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
+ nr1d(k) = crg(2)*org3*qr1d(k)*lamr**bm_r / am_r
+ endif
+ qs1d(k) = qs1d(k) + qsten(k)*DT
+ if (qs1d(k) .le. R1) qs1d(k) = 0.0
+ qg1d(k) = qg1d(k) + qgten(k)*DT
+ if (qg1d(k) .le. R1) qg1d(k) = 0.0
+ enddo
+
+ end subroutine mp_thompson
+!+---+-----------------------------------------------------------------+
+!ctrlL
+!+---+-----------------------------------------------------------------+
+!..Creation of the lookup tables and support functions found below here.
+!+---+-----------------------------------------------------------------+
+!..Rain collecting graupel (and inverse). Explicit CE integration.
+!+---+-----------------------------------------------------------------+
+
+ subroutine qr_acr_qg
+
+ implicit none
+
+!..Local variables
+ INTEGER:: i, j, k, m, n, n2
+ INTEGER:: km, km_s, km_e
+ DOUBLE PRECISION, DIMENSION(nbg):: vg, N_g
+ DOUBLE PRECISION, DIMENSION(nbr):: vr, N_r
+ DOUBLE PRECISION:: N0_r, N0_g, lam_exp, lamg, lamr
+ DOUBLE PRECISION:: massg, massr, dvg, dvr, t1, t2, z1, z2, y1, y2
+
+!+---+
+
+ do n2 = 1, nbr
+! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
+ vr(n2) = -0.1021 + 4.932E3*Dr(n2) - 0.9551E6*Dr(n2)*Dr(n2) &
+ + 0.07934E9*Dr(n2)*Dr(n2)*Dr(n2) &
+ - 0.002362E12*Dr(n2)*Dr(n2)*Dr(n2)*Dr(n2)
+ enddo
+ do n = 1, nbg
+ vg(n) = av_g*Dg(n)**bv_g
+ enddo
+
+!..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for
+!.. fortran indices. J. Michalakes, 2009Oct30.
+
+#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
+ CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e )
+#else
+ km_s = 0
+ km_e = ntb_r*ntb_r1 - 1
+#endif
+
+ do km = km_s, km_e
+ m = km / ntb_r1 + 1
+ k = mod( km , ntb_r1 ) + 1
+
+ lam_exp = (N0r_exp(k)*am_r*crg(1)/r_r(m))**ore1
+ lamr = lam_exp * (crg(3)*org2*org1)**obmr
+ N0_r = N0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2)
+ do n2 = 1, nbr
+ N_r(n2) = N0_r*Dr(n2)**mu_r *DEXP(-lamr*Dr(n2))*dtr(n2)
+ enddo
+
+ do j = 1, ntb_g
+ do i = 1, ntb_g1
+ lam_exp = (N0g_exp(i)*am_g*cgg(1)/r_g(j))**oge1
+ lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
+ N0_g = N0g_exp(i)/(cgg(2)*lam_exp) * lamg**cge(2)
+ do n = 1, nbg
+ N_g(n) = N0_g*Dg(n)**mu_g * DEXP(-lamg*Dg(n))*dtg(n)
+ enddo
+
+ t1 = 0.0d0
+ t2 = 0.0d0
+ z1 = 0.0d0
+ z2 = 0.0d0
+ y1 = 0.0d0
+ y2 = 0.0d0
+ do n2 = 1, nbr
+ massr = am_r * Dr(n2)**bm_r
+ do n = 1, nbg
+ massg = am_g * Dg(n)**bm_g
+
+ dvg = 0.5d0*((vr(n2) - vg(n)) + DABS(vr(n2)-vg(n)))
+ dvr = 0.5d0*((vg(n) - vr(n2)) + DABS(vg(n)-vr(n2)))
+
+ t1 = t1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
+ *dvg*massg * N_g(n)* N_r(n2)
+ z1 = z1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
+ *dvg*massr * N_g(n)* N_r(n2)
+ y1 = y1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
+ *dvg * N_g(n)* N_r(n2)
+
+ t2 = t2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
+ *dvr*massr * N_g(n)* N_r(n2)
+ y2 = y2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
+ *dvr * N_g(n)* N_r(n2)
+ z2 = z2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
+ *dvr*massg * N_g(n)* N_r(n2)
+ enddo
+ 97 continue
+ enddo
+ tcg_racg(i,j,k,m) = t1
+ tmr_racg(i,j,k,m) = DMIN1(z1, r_r(m)*1.0d0)
+ tcr_gacr(i,j,k,m) = t2
+ tmg_gacr(i,j,k,m) = z2
+ tnr_racg(i,j,k,m) = y1
+ tnr_gacr(i,j,k,m) = y2
+ enddo
+ enddo
+ enddo
+
+!..Note wrf_dm_gatherv expects zero-based km_s, km_e (J. Michalakes, 2009Oct30).
+
+#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
+ CALL wrf_dm_gatherv(tcg_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tmr_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tcr_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tmg_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tnr_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tnr_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE)
+#endif
+
+
+ end subroutine qr_acr_qg
+!+---+-----------------------------------------------------------------+
+!ctrlL
+!+---+-----------------------------------------------------------------+
+!..Rain collecting snow (and inverse). Explicit CE integration.
+!+---+-----------------------------------------------------------------+
+
+ subroutine qr_acr_qs
+
+ implicit none
+
+!..Local variables
+ INTEGER:: i, j, k, m, n, n2
+ INTEGER:: km, km_s, km_e
+ DOUBLE PRECISION, DIMENSION(nbr):: vr, D1, N_r
+ DOUBLE PRECISION, DIMENSION(nbs):: vs, N_s
+ DOUBLE PRECISION:: loga_, a_, b_, second, M0, M2, M3, Mrat, oM3
+ DOUBLE PRECISION:: N0_r, lam_exp, lamr, slam1, slam2
+ DOUBLE PRECISION:: dvs, dvr, masss, massr
+ DOUBLE PRECISION:: t1, t2, t3, t4, z1, z2, z3, z4
+ DOUBLE PRECISION:: y1, y2, y3, y4
+
+!+---+
+
+ do n2 = 1, nbr
+! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
+ vr(n2) = -0.1021 + 4.932E3*Dr(n2) - 0.9551E6*Dr(n2)*Dr(n2) &
+ + 0.07934E9*Dr(n2)*Dr(n2)*Dr(n2) &
+ - 0.002362E12*Dr(n2)*Dr(n2)*Dr(n2)*Dr(n2)
+ D1(n2) = (vr(n2)/av_s)**(1./bv_s)
+ enddo
+ do n = 1, nbs
+ vs(n) = 1.5*av_s*Ds(n)**bv_s * DEXP(-fv_s*Ds(n))
+ enddo
+
+!..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for
+!.. fortran indices. J. Michalakes, 2009Oct30.
+
+#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
+ CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e )
+#else
+ km_s = 0
+ km_e = ntb_r*ntb_r1 - 1
+#endif
+
+ do km = km_s, km_e
+ m = km / ntb_r1 + 1
+ k = mod( km , ntb_r1 ) + 1
+
+ lam_exp = (N0r_exp(k)*am_r*crg(1)/r_r(m))**ore1
+ lamr = lam_exp * (crg(3)*org2*org1)**obmr
+ N0_r = N0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2)
+ do n2 = 1, nbr
+ N_r(n2) = N0_r*Dr(n2)**mu_r * DEXP(-lamr*Dr(n2))*dtr(n2)
+ enddo
+
+ do j = 1, ntb_t
+ do i = 1, ntb_s
+
+!..From the bm_s moment, compute plus one moment. If we are not
+!.. using bm_s=2, then we must transform to the pure 2nd moment
+!.. (variable called "second") and then to the bm_s+1 moment.
+
+ M2 = r_s(i)*oams *1.0d0
+ if (bm_s.gt.2.0-1.E-3 .and. bm_s.lt.2.0+1.E-3) then
+ loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*bm_s &
+ + sa(4)*Tc(j)*bm_s + sa(5)*Tc(j)*Tc(j) &
+ + sa(6)*bm_s*bm_s + sa(7)*Tc(j)*Tc(j)*bm_s &
+ + sa(8)*Tc(j)*bm_s*bm_s + sa(9)*Tc(j)*Tc(j)*Tc(j) &
+ + sa(10)*bm_s*bm_s*bm_s
+ a_ = 10.0**loga_
+ b_ = sb(1) + sb(2)*Tc(j) + sb(3)*bm_s &
+ + sb(4)*Tc(j)*bm_s + sb(5)*Tc(j)*Tc(j) &
+ + sb(6)*bm_s*bm_s + sb(7)*Tc(j)*Tc(j)*bm_s &
+ + sb(8)*Tc(j)*bm_s*bm_s + sb(9)*Tc(j)*Tc(j)*Tc(j) &
+ + sb(10)*bm_s*bm_s*bm_s
+ second = (M2/a_)**(1./b_)
+ else
+ second = M2
+ endif
+
+ loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*cse(1) &
+ + sa(4)*Tc(j)*cse(1) + sa(5)*Tc(j)*Tc(j) &
+ + sa(6)*cse(1)*cse(1) + sa(7)*Tc(j)*Tc(j)*cse(1) &
+ + sa(8)*Tc(j)*cse(1)*cse(1) + sa(9)*Tc(j)*Tc(j)*Tc(j) &
+ + sa(10)*cse(1)*cse(1)*cse(1)
+ a_ = 10.0**loga_
+ b_ = sb(1)+sb(2)*Tc(j)+sb(3)*cse(1) + sb(4)*Tc(j)*cse(1) &
+ + sb(5)*Tc(j)*Tc(j) + sb(6)*cse(1)*cse(1) &
+ + sb(7)*Tc(j)*Tc(j)*cse(1) + sb(8)*Tc(j)*cse(1)*cse(1) &
+ + sb(9)*Tc(j)*Tc(j)*Tc(j)+sb(10)*cse(1)*cse(1)*cse(1)
+ M3 = a_ * second**b_
+
+ oM3 = 1./M3
+ Mrat = M2*(M2*oM3)*(M2*oM3)*(M2*oM3)
+ M0 = (M2*oM3)**mu_s
+ slam1 = M2 * oM3 * Lam0
+ slam2 = M2 * oM3 * Lam1
+
+ do n = 1, nbs
+ N_s(n) = Mrat*(Kap0*DEXP(-slam1*Ds(n)) &
+ + Kap1*M0*Ds(n)**mu_s * DEXP(-slam2*Ds(n)))*dts(n)
+ enddo
+
+ t1 = 0.0d0
+ t2 = 0.0d0
+ t3 = 0.0d0
+ t4 = 0.0d0
+ z1 = 0.0d0
+ z2 = 0.0d0
+ z3 = 0.0d0
+ z4 = 0.0d0
+ y1 = 0.0d0
+ y2 = 0.0d0
+ y3 = 0.0d0
+ y4 = 0.0d0
+ do n2 = 1, nbr
+ massr = am_r * Dr(n2)**bm_r
+ do n = 1, nbs
+ masss = am_s * Ds(n)**bm_s
+
+ dvs = 0.5d0*((vr(n2) - vs(n)) + DABS(vr(n2)-vs(n)))
+ dvr = 0.5d0*((vs(n) - vr(n2)) + DABS(vs(n)-vr(n2)))
+
+ if (massr .gt. 2.5*masss) then
+ t1 = t1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
+ *dvs*masss * N_s(n)* N_r(n2)
+ z1 = z1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
+ *dvs*massr * N_s(n)* N_r(n2)
+ y1 = y1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
+ *dvs * N_s(n)* N_r(n2)
+ else
+ t3 = t3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
+ *dvs*masss * N_s(n)* N_r(n2)
+ z3 = z3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
+ *dvs*massr * N_s(n)* N_r(n2)
+ y3 = y3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
+ *dvs * N_s(n)* N_r(n2)
+ endif
+
+ if (massr .gt. 2.5*masss) then
+ t2 = t2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
+ *dvr*massr * N_s(n)* N_r(n2)
+ y2 = y2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
+ *dvr * N_s(n)* N_r(n2)
+ z2 = z2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
+ *dvr*masss * N_s(n)* N_r(n2)
+ else
+ t4 = t4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
+ *dvr*massr * N_s(n)* N_r(n2)
+ y4 = y4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
+ *dvr * N_s(n)* N_r(n2)
+ z4 = z4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
+ *dvr*masss * N_s(n)* N_r(n2)
+ endif
+
+ enddo
+ enddo
+ tcs_racs1(i,j,k,m) = t1
+ tmr_racs1(i,j,k,m) = DMIN1(z1, r_r(m)*1.0d0)
+ tcs_racs2(i,j,k,m) = t3
+ tmr_racs2(i,j,k,m) = z3
+ tcr_sacr1(i,j,k,m) = t2
+ tms_sacr1(i,j,k,m) = z2
+ tcr_sacr2(i,j,k,m) = t4
+ tms_sacr2(i,j,k,m) = z4
+ tnr_racs1(i,j,k,m) = y1
+ tnr_racs2(i,j,k,m) = y3
+ tnr_sacr1(i,j,k,m) = y2
+ tnr_sacr2(i,j,k,m) = y4
+ enddo
+ enddo
+ enddo
+
+!..Note wrf_dm_gatherv expects zero-based km_s, km_e (J. Michalakes, 2009Oct30).
+
+#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
+ CALL wrf_dm_gatherv(tcs_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tmr_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tcs_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tmr_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tcr_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tms_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tcr_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tms_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tnr_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tnr_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tnr_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
+ CALL wrf_dm_gatherv(tnr_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
+#endif
+
+
+ end subroutine qr_acr_qs
+!+---+-----------------------------------------------------------------+
+!ctrlL
+!+---+-----------------------------------------------------------------+
+!..This is a literal adaptation of Bigg (1954) probability of drops of
+!..a particular volume freezing. Given this probability, simply freeze
+!..the proportion of drops summing their masses.
+!+---+-----------------------------------------------------------------+
+
+ subroutine freezeH2O
+
+ implicit none
+
+!..Local variables
+ INTEGER:: i, j, k, n, n2
+ DOUBLE PRECISION, DIMENSION(nbr):: N_r, massr
+ DOUBLE PRECISION, DIMENSION(nbc):: N_c, massc
+ DOUBLE PRECISION:: sum1, sum2, sumn1, sumn2, &
+ prob, vol, Texp, orho_w, &
+ lam_exp, lamr, N0_r, lamc, N0_c, y
+
+!+---+
+
+ orho_w = 1./rho_w
+
+ do n2 = 1, nbr
+ massr(n2) = am_r*Dr(n2)**bm_r
+ enddo
+ do n = 1, nbc
+ massc(n) = am_r*Dc(n)**bm_r
+ enddo
+
+!..Freeze water (smallest drops become cloud ice, otherwise graupel).
+ do k = 1, 45
+! print*, ' Freezing water for temp = ', -k
+ Texp = DEXP( DFLOAT(k) ) - 1.0D0
+ do j = 1, ntb_r1
+ do i = 1, ntb_r
+ lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(i))**ore1
+ lamr = lam_exp * (crg(3)*org2*org1)**obmr
+ N0_r = N0r_exp(j)/(crg(2)*lam_exp) * lamr**cre(2)
+ sum1 = 0.0d0
+ sum2 = 0.0d0
+ sumn1 = 0.0d0
+ sumn2 = 0.0d0
+ do n2 = 1, nbr
+ N_r(n2) = N0_r*Dr(n2)**mu_r*DEXP(-lamr*Dr(n2))*dtr(n2)
+ vol = massr(n2)*orho_w
+ prob = 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp)
+ if (massr(n2) .lt. xm0g) then
+ sumn1 = sumn1 + prob*N_r(n2)
+ sum1 = sum1 + prob*N_r(n2)*massr(n2)
+ else
+ sumn2 = sumn2 + prob*N_r(n2)
+ sum2 = sum2 + prob*N_r(n2)*massr(n2)
+ endif
+ enddo
+ tpi_qrfz(i,j,k) = sum1
+ tni_qrfz(i,j,k) = sumn1
+ tpg_qrfz(i,j,k) = sum2
+ tnr_qrfz(i,j,k) = sumn2
+ enddo
+ enddo
+ do i = 1, ntb_c
+ lamc = 1.0D-6 * (Nt_c*am_r* ccg(2) * ocg1 / r_c(i))**obmr
+ N0_c = 1.0D-18 * Nt_c*ocg1 * lamc**cce(1)
+ sum1 = 0.0d0
+ sumn2 = 0.0d0
+ do n = 1, nbc
+ y = Dc(n)*1.0D6
+ vol = massc(n)*orho_w
+ prob = 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp)
+ N_c(n) = N0_c* y**mu_c * EXP(-lamc*y)*dtc(n)
+ N_c(n) = 1.0D24 * N_c(n)
+ sumn2 = sumn2 + prob*N_c(n)
+ sum1 = sum1 + prob*N_c(n)*massc(n)
+ enddo
+ tpi_qcfz(i,k) = sum1
+ tni_qcfz(i,k) = sumn2
+ enddo
+ enddo
+
+ end subroutine freezeH2O
+!+---+-----------------------------------------------------------------+
+!ctrlL
+!+---+-----------------------------------------------------------------+
+!..Cloud ice converting to snow since portion greater than min snow
+!.. size. Given cloud ice content (kg/m**3), number concentration
+!.. (#/m**3) and gamma shape parameter, mu_i, break the distrib into
+!.. bins and figure out the mass/number of ice with sizes larger than
+!.. D0s. Also, compute incomplete gamma function for the integration
+!.. of ice depositional growth from diameter=0 to D0s. Amount of
+!.. ice depositional growth is this portion of distrib while larger
+!.. diameters contribute to snow growth (as in Harrington et al. 1995).
+!+---+-----------------------------------------------------------------+
+
+ subroutine qi_aut_qs
+
+ implicit none
+
+!..Local variables
+ INTEGER:: i, j, n2
+ DOUBLE PRECISION, DIMENSION(nbi):: N_i
+ DOUBLE PRECISION:: N0_i, lami, Di_mean, t1, t2
+
+!+---+
+
+ do j = 1, ntb_i1
+ do i = 1, ntb_i
+ lami = (am_i*cig(2)*oig1*Nt_i(j)/r_i(i))**obmi
+ Di_mean = (bm_i + mu_i + 1.) / lami
+ N0_i = Nt_i(j)*oig1 * lami**cie(1)
+ t1 = 0.0d0
+ t2 = 0.0d0
+ if (SNGL(Di_mean) .gt. 5.*D0s) then
+ t1 = r_i(i)
+ t2 = Nt_i(j)
+ tpi_ide(i,j) = 0.0D0
+ elseif (SNGL(Di_mean) .lt. D0i) then
+ t1 = 0.0D0
+ t2 = 0.0D0
+ tpi_ide(i,j) = 1.0D0
+ else
+#if (DWORDSIZE == 8 && RWORDSIZE == 8)
+ tpi_ide(i,j) = GAMMP(mu_i+2.0, REAL(lami,KIND=8)*D0s) * 1.0D0
+#elif (DWORDSIZE == 8 && RWORDSIZE == 4)
+ tpi_ide(i,j) = GAMMP(mu_i+2.0, REAL(lami,KIND=4)*D0s) * 1.0D0
+#else
+! This is a temporary hack assuming double precision is 8 bytes.
+#endif
+ do n2 = 1, nbi
+ N_i(n2) = N0_i*Di(n2)**mu_i * DEXP(-lami*Di(n2))*dti(n2)
+ if (Di(n2).ge.D0s) then
+ t1 = t1 + N_i(n2) * am_i*Di(n2)**bm_i
+ t2 = t2 + N_i(n2)
+ endif
+ enddo
+ endif
+ tps_iaus(i,j) = t1
+ tni_iaus(i,j) = t2
+ enddo
+ enddo
+
+ end subroutine qi_aut_qs
+!ctrlL
+!+---+-----------------------------------------------------------------+
+!..Variable collision efficiency for rain collecting cloud water using
+!.. method of Beard and Grover, 1974 if a/A less than 0.25; otherwise
+!.. uses polynomials to get close match of Pruppacher & Klett Fig 14-9.
+!+---+-----------------------------------------------------------------+
+
+ subroutine table_Efrw
+
+ implicit none
+
+!..Local variables
+ DOUBLE PRECISION:: vtr, stokes, reynolds, Ef_rw
+ DOUBLE PRECISION:: p, yc0, F, G, H, z, K0, X
+ INTEGER:: i, j
+
+ do j = 1, nbc
+ do i = 1, nbr
+ Ef_rw = 0.0
+ p = Dc(j)/Dr(i)
+ if (Dr(i).lt.50.E-6 .or. Dc(j).lt.3.E-6) then
+ t_Efrw(i,j) = 0.0
+ elseif (p.gt.0.25) then
+ X = Dc(j)*1.D6
+ if (Dr(i) .lt. 75.e-6) then
+ Ef_rw = 0.026794*X - 0.20604
+ elseif (Dr(i) .lt. 125.e-6) then
+ Ef_rw = -0.00066842*X*X + 0.061542*X - 0.37089
+ elseif (Dr(i) .lt. 175.e-6) then
+ Ef_rw = 4.091e-06*X*X*X*X - 0.00030908*X*X*X &
+ + 0.0066237*X*X - 0.0013687*X - 0.073022
+ elseif (Dr(i) .lt. 250.e-6) then
+ Ef_rw = 9.6719e-5*X*X*X - 0.0068901*X*X + 0.17305*X &
+ - 0.65988
+ elseif (Dr(i) .lt. 350.e-6) then
+ Ef_rw = 9.0488e-5*X*X*X - 0.006585*X*X + 0.16606*X &
+ - 0.56125
+ else
+ Ef_rw = 0.00010721*X*X*X - 0.0072962*X*X + 0.1704*X &
+ - 0.46929
+ endif
+ else
+ vtr = -0.1021 + 4.932E3*Dr(i) - 0.9551E6*Dr(i)*Dr(i) &
+ + 0.07934E9*Dr(i)*Dr(i)*Dr(i) &
+ - 0.002362E12*Dr(i)*Dr(i)*Dr(i)*Dr(i)
+ stokes = Dc(j)*Dc(j)*vtr*rho_w/(9.*1.718E-5*Dr(i))
+ reynolds = 9.*stokes/(p*p*rho_w)
+
+ F = DLOG(reynolds)
+ G = -0.1007D0 - 0.358D0*F + 0.0261D0*F*F
+ K0 = DEXP(G)
+ z = DLOG(stokes/(K0+1.D-15))
+ H = 0.1465D0 + 1.302D0*z - 0.607D0*z*z + 0.293D0*z*z*z
+ yc0 = 2.0D0/PI * ATAN(H)
+ Ef_rw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p))
+
+ endif
+
+ t_Efrw(i,j) = MAX(0.0, MIN(SNGL(Ef_rw), 0.95))
+
+ enddo
+ enddo
+
+ end subroutine table_Efrw
+!ctrlL
+!+---+-----------------------------------------------------------------+
+!..Variable collision efficiency for snow collecting cloud water using
+!.. method of Wang and Ji, 2000 except equate melted snow diameter to
+!.. their "effective collision cross-section."
+!+---+-----------------------------------------------------------------+
+
+ subroutine table_Efsw
+
+ implicit none
+
+!..Local variables
+ DOUBLE PRECISION:: Ds_m, vts, vtc, stokes, reynolds, Ef_sw
+ DOUBLE PRECISION:: p, yc0, F, G, H, z, K0
+ INTEGER:: i, j
+
+ do j = 1, nbc
+ vtc = 1.19D4 * (1.0D4*Dc(j)*Dc(j)*0.25D0)
+ do i = 1, nbs
+ vts = av_s*Ds(i)**bv_s * DEXP(-fv_s*Ds(i)) - vtc
+ Ds_m = (am_s*Ds(i)**bm_s / am_r)**obmr
+ p = Dc(j)/Ds_m
+ if (p.gt.0.25 .or. Ds(i).lt.D0s .or. Dc(j).lt.6.E-6 &
+ .or. vts.lt.1.E-3) then
+ t_Efsw(i,j) = 0.0
+ else
+ stokes = Dc(j)*Dc(j)*vts*rho_w/(9.*1.718E-5*Ds_m)
+ reynolds = 9.*stokes/(p*p*rho_w)
+
+ F = DLOG(reynolds)
+ G = -0.1007D0 - 0.358D0*F + 0.0261D0*F*F
+ K0 = DEXP(G)
+ z = DLOG(stokes/(K0+1.D-15))
+ H = 0.1465D0 + 1.302D0*z - 0.607D0*z*z + 0.293D0*z*z*z
+ yc0 = 2.0D0/PI * ATAN(H)
+ Ef_sw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p))
+
+ t_Efsw(i,j) = MAX(0.0, MIN(SNGL(Ef_sw), 0.95))
+ endif
+
+ enddo
+ enddo
+
+ end subroutine table_Efsw
+!ctrlL
+!+---+-----------------------------------------------------------------+
+!..Integrate rain size distribution from zero to D-star to compute the
+!.. number of drops smaller than D-star that evaporate in a single
+!.. timestep. Drops larger than D-star dont evaporate entirely so do
+!.. not affect number concentration.
+!+---+-----------------------------------------------------------------+
+
+ subroutine table_dropEvap
+
+ implicit none
+
+!..Local variables
+ DOUBLE PRECISION:: Nt_r, N0, lam_exp, lam
+ INTEGER:: i, j, k
+
+ do k = 1, ntb_r
+ do j = 1, ntb_r1
+ lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(k))**ore1
+ lam = lam_exp * (crg(3)*org2*org1)**obmr
+ N0 = N0r_exp(j)/(crg(2)*lam_exp) * lam**cre(2)
+ Nt_r = N0 * crg(2) / lam**cre(2)
+
+ do i = 1, nbr
+#if (DWORDSIZE == 8 && RWORDSIZE == 8)
+ tnr_rev(i,j,k) = GAMMP(mu_r+1.0, REAL(Dr(i)*lam,KIND=8)) * Nt_r
+#elif (DWORDSIZE == 8 && RWORDSIZE == 4)
+ tnr_rev(i,j,k) = GAMMP(mu_r+1.0, REAL(Dr(i)*lam,KIND=4)) * Nt_r
+#else
+! This is a temporary hack assuming double precision is 8 bytes.
+#endif
+ enddo
+
+ enddo
+ enddo
+
+ end subroutine table_dropEvap
+
+! TO APPLY TABLE ABOVE
+!..Rain lookup table indexes.
+! Dr_star = DSQRT(-2.D0*DT * t1_evap/(2.*PI) &
+! * 0.78*4.*diffu(k)*xsat*rvs/rho_w)
+! idx_d = NINT(1.0 + FLOAT(nbr) * DLOG(Dr_star/D0r) &
+! / DLOG(Dr(nbr)/D0r))
+! idx_d = MAX(1, MIN(idx_d, nbr))
+!
+! nir = NINT(ALOG10(rr(k)))
+! do nn = nir-1, nir+1
+! n = nn
+! if ( (rr(k)/10.**nn).ge.1.0 .and. &
+! (rr(k)/10.**nn).lt.10.0) goto 154
+! enddo
+!154 continue
+! idx_r = INT(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2)
+! idx_r = MAX(1, MIN(idx_r, ntb_r))
+!
+! lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
+! lam_exp = lamr * (crg(3)*org2*org1)**bm_r
+! N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
+! nir = NINT(DLOG10(N0_exp))
+! do nn = nir-1, nir+1
+! n = nn
+! if ( (N0_exp/10.**nn).ge.1.0 .and. &
+! (N0_exp/10.**nn).lt.10.0) goto 155
+! enddo
+!155 continue
+! idx_r1 = INT(N0_exp/10.**n) + 10*(n-nir3) - (n-nir3)
+! idx_r1 = MAX(1, MIN(idx_r1, ntb_r1))
+!
+! pnr_rev(k) = MIN(nr(k)*odts, SNGL(tnr_rev(idx_d,idx_r1,idx_r) & ! RAIN2M
+! * odts))
+!
+!ctrlL
+!+---+-----------------------------------------------------------------+
+!+---+-----------------------------------------------------------------+
+ SUBROUTINE GCF(GAMMCF,A,X,GLN)
+! --- RETURNS THE INCOMPLETE GAMMA FUNCTION Q(A,X) EVALUATED BY ITS
+! --- CONTINUED FRACTION REPRESENTATION AS GAMMCF. ALSO RETURNS
+! --- LN(GAMMA(A)) AS GLN. THE CONTINUED FRACTION IS EVALUATED BY
+! --- A MODIFIED LENTZ METHOD.
+! --- USES GAMMLN
+ IMPLICIT NONE
+ INTEGER, PARAMETER:: ITMAX=100
+ REAL, PARAMETER:: gEPS=3.E-7
+ REAL, PARAMETER:: FPMIN=1.E-30
+ REAL, INTENT(IN):: A, X
+ REAL:: GAMMCF,GLN
+ INTEGER:: I
+ REAL:: AN,B,C,D,DEL,H
+ GLN=GAMMLN(A)
+ B=X+1.-A
+ C=1./FPMIN
+ D=1./B
+ H=D
+ DO 11 I=1,ITMAX
+ AN=-I*(I-A)
+ B=B+2.
+ D=AN*D+B
+ IF(ABS(D).LT.FPMIN)D=FPMIN
+ C=B+AN/C
+ IF(ABS(C).LT.FPMIN)C=FPMIN
+ D=1./D
+ DEL=D*C
+ H=H*DEL
+ IF(ABS(DEL-1.).LT.gEPS)GOTO 1
+ 11 CONTINUE
+ PRINT *, 'A TOO LARGE, ITMAX TOO SMALL IN GCF'
+ 1 GAMMCF=EXP(-X+A*LOG(X)-GLN)*H
+ END SUBROUTINE GCF
+! (C) Copr. 1986-92 Numerical Recipes Software 2.02
+!+---+-----------------------------------------------------------------+
+ SUBROUTINE GSER(GAMSER,A,X,GLN)
+! --- RETURNS THE INCOMPLETE GAMMA FUNCTION P(A,X) EVALUATED BY ITS
+! --- ITS SERIES REPRESENTATION AS GAMSER. ALSO RETURNS LN(GAMMA(A))
+! --- AS GLN.
+! --- USES GAMMLN
+ IMPLICIT NONE
+ INTEGER, PARAMETER:: ITMAX=100
+ REAL, PARAMETER:: gEPS=3.E-7
+ REAL, INTENT(IN):: A, X
+ REAL:: GAMSER,GLN
+ INTEGER:: N
+ REAL:: AP,DEL,SUM
+ GLN=GAMMLN(A)
+ IF(X.LE.0.)THEN
+ IF(X.LT.0.) PRINT *, 'X < 0 IN GSER'
+ GAMSER=0.
+ RETURN
+ ENDIF
+ AP=A
+ SUM=1./A
+ DEL=SUM
+ DO 11 N=1,ITMAX
+ AP=AP+1.
+ DEL=DEL*X/AP
+ SUM=SUM+DEL
+ IF(ABS(DEL).LT.ABS(SUM)*gEPS)GOTO 1
+ 11 CONTINUE
+ PRINT *,'A TOO LARGE, ITMAX TOO SMALL IN GSER'
+ 1 GAMSER=SUM*EXP(-X+A*LOG(X)-GLN)
+ END SUBROUTINE GSER
+! (C) Copr. 1986-92 Numerical Recipes Software 2.02
+!+---+-----------------------------------------------------------------+
+ REAL FUNCTION GAMMLN(XX)
+! --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 0.
+ IMPLICIT NONE
+ REAL, INTENT(IN):: XX
+ DOUBLE PRECISION, PARAMETER:: STP = 2.5066282746310005D0
+ DOUBLE PRECISION, DIMENSION(6), PARAMETER:: &
+ COF = (/76.18009172947146D0, -86.50532032941677D0, &
+ 24.01409824083091D0, -1.231739572450155D0, &
+ .1208650973866179D-2, -.5395239384953D-5/)
+ DOUBLE PRECISION:: SER,TMP,X,Y
+ INTEGER:: J
+
+ X=XX
+ Y=X
+ TMP=X+5.5D0
+ TMP=(X+0.5D0)*LOG(TMP)-TMP
+ SER=1.000000000190015D0
+ DO 11 J=1,6
+ Y=Y+1.D0
+ SER=SER+COF(J)/Y
+11 CONTINUE
+ GAMMLN=TMP+LOG(STP*SER/X)
+ END FUNCTION GAMMLN
+! (C) Copr. 1986-92 Numerical Recipes Software 2.02
+!+---+-----------------------------------------------------------------+
+ REAL FUNCTION GAMMP(A,X)
+! --- COMPUTES THE INCOMPLETE GAMMA FUNCTION P(A,X)
+! --- SEE ABRAMOWITZ AND STEGUN 6.5.1
+! --- USES GCF,GSER
+ IMPLICIT NONE
+ REAL, INTENT(IN):: A,X
+ REAL:: GAMMCF,GAMSER,GLN
+ GAMMP = 0.
+ IF((X.LT.0.) .OR. (A.LE.0.)) THEN
+ PRINT *, 'BAD ARGUMENTS IN GAMMP'
+ RETURN
+ ELSEIF(X.LT.A+1.)THEN
+ CALL GSER(GAMSER,A,X,GLN)
+ GAMMP=GAMSER
+ ELSE
+ CALL GCF(GAMMCF,A,X,GLN)
+ GAMMP=1.-GAMMCF
+ ENDIF
+ END FUNCTION GAMMP
+! (C) Copr. 1986-92 Numerical Recipes Software 2.02
+!+---+-----------------------------------------------------------------+
+ REAL FUNCTION WGAMMA(y)
+
+ IMPLICIT NONE
+ REAL, INTENT(IN):: y
+
+ WGAMMA = EXP(GAMMLN(y))
+
+ END FUNCTION WGAMMA
+!+---+-----------------------------------------------------------------+
+! THIS FUNCTION CALCULATES THE LIQUID SATURATION VAPOR MIXING RATIO AS
+! A FUNCTION OF TEMPERATURE AND PRESSURE
+!
+ REAL FUNCTION RSLF(P,T)
+
+ IMPLICIT NONE
+ REAL, INTENT(IN):: P, T
+ REAL:: ESL,X
+ REAL, PARAMETER:: C0= .611583699E03
+ REAL, PARAMETER:: C1= .444606896E02
+ REAL, PARAMETER:: C2= .143177157E01
+ REAL, PARAMETER:: C3= .264224321E-1
+ REAL, PARAMETER:: C4= .299291081E-3
+ REAL, PARAMETER:: C5= .203154182E-5
+ REAL, PARAMETER:: C6= .702620698E-8
+ REAL, PARAMETER:: C7= .379534310E-11
+ REAL, PARAMETER:: C8=-.321582393E-13
+
+ X=MAX(-80.,T-273.16)
+
+! ESL=612.2*EXP(17.67*X/(T-29.65))
+ ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
+ RSLF=.622*ESL/(P-ESL)
+
+! ALTERNATIVE
+! ; Source: Murphy and Koop, Review of the vapour pressure of ice and
+! supercooled water for atmospheric applications, Q. J. R.
+! Meteorol. Soc (2005), 131, pp. 1539-1565.
+! ESL = EXP(54.842763 - 6763.22 / T - 4.210 * ALOG(T) + 0.000367 * T
+! + TANH(0.0415 * (T - 218.8)) * (53.878 - 1331.22
+! / T - 9.44523 * ALOG(T) + 0.014025 * T))
+
+ END FUNCTION RSLF
+!+---+-----------------------------------------------------------------+
+! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A
+! FUNCTION OF TEMPERATURE AND PRESSURE
+!
+ REAL FUNCTION RSIF(P,T)
+
+ IMPLICIT NONE
+ REAL, INTENT(IN):: P, T
+ REAL:: ESI,X
+ REAL, PARAMETER:: C0= .609868993E03
+ REAL, PARAMETER:: C1= .499320233E02
+ REAL, PARAMETER:: C2= .184672631E01
+ REAL, PARAMETER:: C3= .402737184E-1
+ REAL, PARAMETER:: C4= .565392987E-3
+ REAL, PARAMETER:: C5= .521693933E-5
+ REAL, PARAMETER:: C6= .307839583E-7
+ REAL, PARAMETER:: C7= .105785160E-9
+ REAL, PARAMETER:: C8= .161444444E-12
+
+ X=MAX(-80.,T-273.16)
+ ESI=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
+ RSIF=.622*ESI/(P-ESI)
+
+! ALTERNATIVE
+! ; Source: Murphy and Koop, Review of the vapour pressure of ice and
+! supercooled water for atmospheric applications, Q. J. R.
+! Meteorol. Soc (2005), 131, pp. 1539-1565.
+! ESI = EXP(9.550426 - 5723.265/T + 3.53068*ALOG(T) - 0.00728332*T)
+
+ END FUNCTION RSIF
+!+---+-----------------------------------------------------------------+
+
+!+---+-----------------------------------------------------------------+
+END MODULE module_mp_thompson
+!+---+-----------------------------------------------------------------+
</font>
</pre>