!+---+-----------------------------------------------------------------+
!.. 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 et al. (2004, 2007).
!.. 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: 26 Oct 2007
!+---+-----------------------------------------------------------------+
!wrft:model_layer:physics
!+---+-----------------------------------------------------------------+
!
MODULE module_mp_thompson07 
(docs) 2
USE module_wrf_error
IMPLICIT NONE
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 parameters for rain & graupel. However, these are not
!.. constant and vary depending on mixing ratio. Furthermore, when
!.. mu is non-zero, these become equiv y-intercept for an exponential
!.. distrib and proper values computed based on assumed mu value.
REAL, PARAMETER, PRIVATE:: gonv_min = 1.E4
REAL, PARAMETER, PRIVATE:: gonv_max = 5.E6
REAL, PARAMETER, PRIVATE:: ronv_min = 2.E6
REAL, PARAMETER, PRIVATE:: ronv_max = 9.E9
REAL, PARAMETER, PRIVATE:: ronv_sl = 1./4.
REAL, PARAMETER, PRIVATE:: ronv_r0 = 0.10E-3
REAL, PARAMETER, PRIVATE:: ronv_c0 = ronv_sl/ronv_r0
REAL, PARAMETER, PRIVATE:: ronv_c1 = (ronv_max-ronv_min)*0.5
REAL, PARAMETER, PRIVATE:: ronv_c2 = (ronv_max+ronv_min)*0.5
!..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 (2006). 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.1
REAL, PARAMETER, PRIVATE:: Ef_rs = 0.95
REAL, PARAMETER, PRIVATE:: Ef_rg = 0.95
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_r1 = 37
INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55
INTEGER, PARAMETER, PRIVATE:: ntb_t = 9
INTEGER, PRIVATE:: nic2, nii2, nii3, nir2, nir3, nis2, nig2
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 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 tp/tc/tm
!.. represent lookup tables.
DOUBLE PRECISION, DIMENSION(ntb_g,ntb_r1,ntb_r):: &
tcg_racg, tmr_racg, tcr_gacr, tmg_gacr
DOUBLE PRECISION, DIMENSION(ntb_s,ntb_t,ntb_r1,ntb_r):: &
tcs_racs1, tmr_racs1, tcs_racs2, tmr_racs2, &
tcr_sacr1, tms_sacr1, tcr_sacr2, tms_sacr2
DOUBLE PRECISION, DIMENSION(ntb_c,45):: &
tpi_qcfz, tni_qcfz
DOUBLE PRECISION, DIMENSION(ntb_r,ntb_r1,45):: &
tpi_qrfz, tpg_qrfz, tni_qrfz
DOUBLE PRECISION, DIMENSION(ntb_i,ntb_i1):: &
tps_iaus, tni_iaus, tpi_ide
REAL, DIMENSION(nbr,nbc):: t_Efrw
REAL, DIMENSION(nbs,nbc):: t_Efsw
!..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(12), 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
!+---+-----------------------------------------------------------------+
!+---+
!
CONTAINS
SUBROUTINE thompson07_init (docs) 1,26
IMPLICIT NONE
INTEGER:: i, j, k, m, n
!..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) = bm_r + mu_r + 3.
cre(6) = bm_r + mu_r + bv_r + 1.
cre(7) = bm_r + mu_r + bv_r + 2.
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 + mu_r + 4.
do n = 1, 12
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 + bv_s + 2.
cse(6) = bm_s + bv_s + 3.
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 + 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 + 2.
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 + mu_g + 3.
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)))
!..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 k = 1, ntb_r
do j = 1, ntb_r1
do i = 1, ntb_g
tcg_racg(i,j,k) = 0.0d0
tmr_racg(i,j,k) = 0.0d0
tcr_gacr(i,j,k) = 0.0d0
tmg_gacr(i,j,k) = 0.0d0
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
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
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
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
if (.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')
END SUBROUTINE thompson07_init
!+---+-----------------------------------------------------------------+
!
!+---+-----------------------------------------------------------------+
!..This is a wrapper routine designed to transfer values from 3D to 1D.
!+---+-----------------------------------------------------------------+
SUBROUTINE mp_gt_driver07 (docs) (qv, qc, qr, qi, qs, qg, ni, & 1,9
th, pii, p, dz, dt_in, itimestep, &
RAINNC, RAINNCV, SR, &
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, 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, INTENT(IN):: dt_in
INTEGER, INTENT(IN):: itimestep
!..Local variables
REAL, DIMENSION(kts:kte):: &
qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
t1d, p1d, dz1d
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
INTEGER:: i, j, k
INTEGER:: imax_qc, imax_qr, imax_qi, imax_qs, imax_qg, imax_ni
INTEGER:: jmax_qc, jmax_qr, jmax_qi, jmax_qs, jmax_qg, jmax_ni
INTEGER:: kmax_qc, kmax_qr, kmax_qi, kmax_qs, kmax_qg, kmax_ni
INTEGER:: i_start, j_start, i_end, j_end
!+---+
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 + 1
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 + 1
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.
imax_qc = 0
imax_qr = 0
imax_qi = 0
imax_qs = 0
imax_qg = 0
imax_ni = 0
jmax_qc = 0
jmax_qr = 0
jmax_qi = 0
jmax_qs = 0
jmax_qg = 0
jmax_ni = 0
kmax_qc = 0
kmax_qr = 0
kmax_qi = 0
kmax_qs = 0
kmax_qg = 0
kmax_ni = 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.
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)
enddo
call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
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
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)
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 (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
write(mp_debug,*) 'WARNING, negative qv ', qv1d(k), &
' at i,j,k=', i,j,k
CALL wrf_debug(150, mp_debug)
endif
enddo
enddo i_loop
enddo j_loop
! DEBUG - GT
write(mp_debug,'(a,6(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, ')'
CALL wrf_debug(150, mp_debug)
! END DEBUG - GT
do i = 1, 256
mp_debug(i:i) = char(0)
enddo
END SUBROUTINE mp_gt_driver07
!+---+-----------------------------------------------------------------+
!
!+---+-----------------------------------------------------------------+
!+---+-----------------------------------------------------------------+
!.. 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, 2006).
!+---+-----------------------------------------------------------------+
!
subroutine mp_thompson (docs) (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & 2,8
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, &
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
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
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
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, &
smoc, smod, smoe, smof
REAL, DIMENSION(kts:kte):: sed_r, sed_s, sed_g, sed_i, sed_n
REAL:: rgvm, delta_tp, orho, onstep, lfus2
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
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, vtnk, vtrk, 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
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
REAL:: dtsave, odts, odt, odzq
INTEGER:: i, k, k2, ksed1, ku, n, nn, nstep, k_0, kbot, IT, iexfrq
INTEGER:: nir, nis, nig, nii, nic
INTEGER:: idx_tc,idx_t,idx_s,idx_g,idx_r1,idx_r,idx_i1,idx_i,idx_c
INTEGER:: idx
LOGICAL:: melti, no_micro
LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg
!+---+
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.
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.
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. 30.E-6) then
lami = cie(2)/30.E-6
ni(k) = MIN(250.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)
L_qr(k) = .true.
else
qr1d(k) = 0.0
rr(k) = R1
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 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.
!+---+-----------------------------------------------------------------+
N0_min = gonv_max
do k = kte, kts, -1
if (.not. L_qg(k)) CYCLE
N0_exp = 200.0*rho(k)/rg(k)
N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max)))
N0_min = MIN(N0_exp, N0_min)
N0_exp = N0_min
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.
!.. New treatment for variable y-intercept of rain. When rain comes
!.. from melted snow/graupel, compute mass-weighted mean size, melt
!.. into water, compute its mvd and recompute slope/intercept.
!.. If rain not from melted snow, use old relation but hold N0_r
!.. constant at its lowest value. While doing all this, ensure rain
!.. mvd does not exceed reasonable size like 2.5 mm.
!+---+-----------------------------------------------------------------+
N0_min = ronv_max
do k = kte, kts, -1
! if (.not. L_qr(k)) CYCLE
N0_exp = ronv_c1*tanh(ronv_c0*(ronv_r0-rr(k))) + ronv_c2
N0_min = MIN(N0_exp, N0_min)
N0_exp = N0_min
lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
lamr = lam_exp * (crg(3)*org2*org1)**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) / 2.5e-3
lam_exp = lamr * (crg(3)*org2*org1)**bm_r
N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
endif
N0_r(k) = N0_exp/(crg(2)*lam_exp) * lamr**cre(2)
ilamr(k) = 1./lamr
enddo
if (.not. iiwarm) then
k_0 = kts
melti = .false.
do k = kte-1, kts, -1
if ( (temp(k).gt. T_0) .and. (rr(k).gt. 0.001e-3) &
.and. ((rs(k+1)+rg(k+1)).gt. 0.01e-3) ) then
k_0 = MAX(k+1, k_0)
melti=.true.
goto 135
endif
enddo
135 continue
if (melti) then
!.. Locate bottom of melting layer (if any).
kbot = kts
do k = k_0-1, kts, -1
if ( (rs(k)+rg(k)).lt. 0.01e-3) goto 136
enddo
136 continue
kbot = MAX(k, kts)
!.. Compute melted snow/graupel equiv water diameter one K-level above
!.. melting. Set starting rain mvd to either 50 microns or max from
!.. higher up in column.
if (L_qs(k_0)) then
xDs = smoc(k_0) / smob(k_0)
Ds_m = (am_s*xDs**bm_s / am_r)**obmr
else
Ds_m = 1.0e-6
endif
if (L_qg(k_0)) then
xDg = (bm_g + mu_g + 1.) * ilamg(k_0)
Dg_m = (am_g*xDg**bm_g / am_r)**obmr
else
Dg_m = 1.0e-6
endif
r_mvd1 = mvd_r(k_0)
r_mvd2 = MIN(MAX(Ds_m, Dg_m, r_mvd1+1.e-6, mvd_r(kbot)), &
2.5e-3)
!.. Within melting layer, apply linear increase of rain mvd from r_mvd1
!.. to equiv melted snow/graupel value (r_mvd2). So, by the bottom of
!.. the melting layer, the rain will have an mvd that matches that from
!.. melted snow and/or graupel.
if (kbot.gt. 2) then
do k = k_0-1, kbot, -1
if (.not. L_qr(k)) CYCLE
xkrat = REAL(k_0-k)/REAL(k_0-kbot)
mvd_r(k) = MAX(mvd_r(k), xkrat*(r_mvd2-r_mvd1)+r_mvd1)
lamr = (4.0+mu_r) / mvd_r(k)
lam_exp = lamr * (crg(3)*org2*org1)**bm_r
N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
N0_exp = MAX(DBLE(ronv_min), MIN(N0_exp, DBLE(ronv_max)))
lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
lamr = lam_exp * (crg(3)*org2*org1)**obmr
mvd_r(k) = (3.0+mu_r+0.672) / lamr
N0_r(k) = rr(k)*lamr**cre(3) / (am_r*crg(3))
ilamr(k) = 1./lamr
enddo
!.. Below melting layer, hold N0_r constant unless changes to mixing
!.. ratio increase mvd beyond 2.5 mm threshold, then adjust slope and
!.. intercept to cap mvd at 2.5 mm. In future, we could lower N0_r to
!.. account for self-collection or other sinks.
do k = kbot-1, kts, -1
if (.not. L_qr(k)) CYCLE
N0_r(k) = MIN(N0_r(k), N0_r(kbot))
lamr = (N0_r(k)*am_r*crg(3)/rr(k))**(1./cre(3))
lam_exp = lamr * (crg(3)*org2*org1)**bm_r
N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
N0_exp = MAX(DBLE(ronv_min), MIN(N0_exp, DBLE(ronv_max)))
lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
lamr = lam_exp * (crg(3)*org2*org1)**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)
endif
N0_r(k) = rr(k)*lamr**cre(3) / (am_r*crg(3))
ilamr(k) = 1./lamr
enddo
endif
endif
endif
!+---+-----------------------------------------------------------------+
!..Compute warm-rain process terms (except evap done later).
!+---+-----------------------------------------------------------------+
do k = kts, kte
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))
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))
else
idx_g = 1
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))
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)
endif
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_g,idx_r1,idx_r) &
+ tcr_gacr(idx_g,idx_r1,idx_r)
prg_rcg(k) = MIN(DBLE(rr(k)*odts), prg_rcg(k))
prr_rcg(k) = -prg_rcg(k)
else
prr_rcg(k) = tcg_racg(idx_g,idx_r1,idx_r) &
+ 0.5*tmg_gacr(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
elseif (rr(k).gt. R1 .and. temp(k).lt.HGFR) then
pri_rfz(k) = rr(k)*odts
pni_rfz(k) = N0_r(k)*crg(2) * ilamr(k)**cre(2)
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<-10C or 20%+ ice supersaturated.
if ( (ssati(k).ge. 0.20) .or. (ssatw(k).gt. eps &
.and. temp(k).lt.263.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)))
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))
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))
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
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
!.. 30 and 300 microns. Also no more than 250 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. 30.E-6) then
lami = cie(2)/30.E-6
xni = MIN(250.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.250.E3) &
niten(k) = (250.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
!..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)
L_qr(k) = .true.
else
rr(k) = R1
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.
!+---+-----------------------------------------------------------------+
N0_min = gonv_max
do k = kte, kts, -1
if (.not. L_qg(k)) CYCLE
N0_exp = 200.0*rho(k)/rg(k)
N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max)))
N0_min = MIN(N0_exp, N0_min)
N0_exp = N0_min
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.
!+---+-----------------------------------------------------------------+
N0_min = ronv_max
do k = kte, kts, -1
! if (.not. L_qr(k)) CYCLE
N0_exp = ronv_c1*tanh(ronv_c0*(ronv_r0-rr(k))) + ronv_c2
N0_min = MIN(N0_exp, N0_min)
N0_exp = N0_min
lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
lamr = lam_exp * (crg(3)*org2*org1)**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) / 2.5e-3
lam_exp = lamr * (crg(3)*org2*org1)**bm_r
N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
endif
N0_r(k) = N0_exp/(crg(2)*lam_exp) * lamr**cre(2)
ilamr(k) = 1./lamr
enddo
if (.not. iiwarm) then
k_0 = kts
melti = .false.
do k = kte-1, kts, -1
if ( (temp(k).gt. T_0) .and. (rr(k).gt. 0.001e-3) &
.and. ((rs(k+1)+rg(k+1)).gt. 0.01e-3) ) then
k_0 = MAX(k+1, k_0)
melti=.true.
goto 137
endif
enddo
137 continue
if (melti) then
!.. Locate bottom of melting layer (if any).
kbot = kts
do k = k_0-1, kts, -1
if ( (rs(k)+rg(k)).lt. 0.01e-3) goto 138
enddo
138 continue
kbot = MAX(k, kts)
!.. Compute melted snow/graupel equiv water diameter one K-level above
!.. melting. Set starting rain mvd to either 50 microns or max from
!.. higher up in column.
if (L_qs(k_0)) then
xDs = smoc(k_0) / smob(k_0)
Ds_m = (am_s*xDs**bm_s / am_r)**obmr
else
Ds_m = 1.0e-6
endif
if (L_qg(k_0)) then
xDg = (bm_g + mu_g + 1.) * ilamg(k_0)
Dg_m = (am_g*xDg**bm_g / am_r)**obmr
else
Dg_m = 1.0e-6
endif
r_mvd1 = mvd_r(k_0)
r_mvd2 = MIN(MAX(Ds_m, Dg_m, r_mvd1+1.e-6, mvd_r(kbot)), &
2.5e-3)
!.. Within melting layer, apply linear increase of rain mvd from r_mvd1
!.. to equiv melted snow/graupel value (r_mvd2). So, by the bottom of
!.. the melting layer, the rain will have an mvd that matches that from
!.. melted snow and/or graupel.
if (kbot.gt. 2) then
do k = k_0-1, kbot, -1
if (.not. L_qr(k)) CYCLE
xkrat = REAL(k_0-k)/REAL(k_0-kbot)
mvd_r(k) = MAX(mvd_r(k), xkrat*(r_mvd2-r_mvd1)+r_mvd1)
lamr = (4.0+mu_r) / mvd_r(k)
lam_exp = lamr * (crg(3)*org2*org1)**bm_r
N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
N0_exp = MAX(DBLE(ronv_min), MIN(N0_exp, DBLE(ronv_max)))
lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
lamr = lam_exp * (crg(3)*org2*org1)**obmr
mvd_r(k) = (3.0+mu_r+0.672) / lamr
N0_r(k) = rr(k)*lamr**cre(3) / (am_r*crg(3))
ilamr(k) = 1./lamr
enddo
!.. Below melting layer, hold N0_r constant unless changes to mixing
!.. ratio increase mvd beyond 2.5 mm threshold, then adjust slope and
!.. intercept to cap mvd at 2.5 mm. In future, we could lower N0_r to
!.. account for self-collection or other sinks.
do k = kbot-1, kts, -1
if (.not. L_qr(k)) CYCLE
N0_r(k) = MIN(N0_r(k), N0_r(kbot))
lamr = (N0_r(k)*am_r*crg(3)/rr(k))**(1./cre(3))
lam_exp = lamr * (crg(3)*org2*org1)**bm_r
N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
N0_exp = MAX(DBLE(ronv_min), MIN(N0_exp, DBLE(ronv_max)))
lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
lamr = lam_exp * (crg(3)*org2*org1)**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)
endif
N0_r(k) = rr(k)*lamr**cre(3) / (am_r*crg(3))
ilamr(k) = 1./lamr
enddo
endif
endif
endif
!+---+-----------------------------------------------------------------+
!..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))
qrten(k) = qrten(k) - prv_rev(k)
qvten(k) = qvten(k) + prv_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))
temp(k) = t1d(k) + DT*tten(k)
rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
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).
!+---+-----------------------------------------------------------------+
nstep = 0
ksed1 = 0
do k = kte+1, kts, -1
vtrk(k) = 0.
vtik(k) = 0.
vtnk(k) = 0.
vtsk(k) = 0.
vtgk(k) = 0.
enddo
do k = kte, kts, -1
vtr = 0.
vti = 0.
vts = 0.
vtg = 0.
rhof(k) = SQRT(RHO_NOT/rho(k))
if (rr(k).gt. R2) then
lamr = 1./ilamr(k)
vtr = rhof(k)*av_r*crg(6)*org3 * (lamr/(lamr+fv_r))**cre(3) &
*((lamr+fv_r)**(-bv_r))
vtrk(k) = vtr
endif
if (.not. iiwarm) then
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
vtnk(k) = vti
endif
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)
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 (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
endif
rgvm = MAX(vtik(k), vtrk(k), vtsk(k), vtgk(k))
if (rgvm .gt. 1.E-3) then
ksed1 = MAX(ksed1, k)
delta_tp = dzq(k)/rgvm
nstep = MAX(nstep, INT(DT/delta_tp + 1.))
endif
enddo
if (ksed1 .eq. kte) ksed1 = kte-1
if (nstep .gt. 0) onstep = 1./REAL(nstep)
!+---+-----------------------------------------------------------------+
!..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.
!+---+-----------------------------------------------------------------+
do n = 1, nstep
do k = kte, kts, -1
sed_r(k) = vtrk(k)*rr(k)
sed_i(k) = vtik(k)*ri(k)
sed_n(k) = vtnk(k)*ni(k)
sed_g(k) = vtgk(k)*rg(k)
sed_s(k) = vtsk(k)*rs(k)
enddo
k = kte
odzq = 1./dzq(k)
orho = 1./rho(k)
qrten(k) = qrten(k) - sed_r(k)*odzq*onstep*orho
qiten(k) = qiten(k) - sed_i(k)*odzq*onstep*orho
niten(k) = niten(k) - sed_n(k)*odzq*onstep*orho
qgten(k) = qgten(k) - sed_g(k)*odzq*onstep*orho
qsten(k) = qsten(k) - sed_s(k)*odzq*onstep*orho
rr(k) = MAX(R1, rr(k) - sed_r(k)*odzq*DT*onstep)
ri(k) = MAX(R1, ri(k) - sed_i(k)*odzq*DT*onstep)
ni(k) = MAX(1., ni(k) - sed_n(k)*odzq*DT*onstep)
rg(k) = MAX(R1, rg(k) - sed_g(k)*odzq*DT*onstep)
rs(k) = MAX(R1, rs(k) - sed_s(k)*odzq*DT*onstep)
do k = ksed1, kts, -1
odzq = 1./dzq(k)
orho = 1./rho(k)
qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) &
*odzq*onstep*orho
qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) &
*odzq*onstep*orho
niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) &
*odzq*onstep*orho
qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) &
*odzq*onstep*orho
qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) &
*odzq*onstep*orho
rr(k) = MAX(R1, rr(k) + (sed_r(k+1)-sed_r(k)) &
*odzq*DT*onstep)
ri(k) = MAX(R1, ri(k) + (sed_i(k+1)-sed_i(k)) &
*odzq*DT*onstep)
ni(k) = MAX(1., ni(k) + (sed_n(k+1)-sed_n(k)) &
*odzq*DT*onstep)
rg(k) = MAX(R1, rg(k) + (sed_g(k+1)-sed_g(k)) &
*odzq*DT*onstep)
rs(k) = MAX(R1, rs(k) + (sed_s(k+1)-sed_s(k)) &
*odzq*DT*onstep)
enddo
!+---+-----------------------------------------------------------------+
!..Precipitation reaching the ground.
!+---+-----------------------------------------------------------------+
pptrain = pptrain + sed_r(kts)*DT*onstep
pptsnow = pptsnow + sed_s(kts)*DT*onstep
pptgraul = pptgraul + sed_g(kts)*DT*onstep
pptice = pptice + sed_i(kts)*DT*onstep
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) = -qi1d(k)*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/(2.*xm0i)*odt
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. 30.E-6) then
lami = cie(2)/30.E-6
ni1d(k) = MIN(250.D3, cig(1)*oig2*qi1d(k)/am_i*lami**bm_i)
elseif (xDi.gt. 300.E-6) then
lami = cie(2)/300.E-6
ni1d(k) = cig(1)*oig2*qi1d(k)/am_i*lami**bm_i
endif
else
lami = cie(2)/30.E-6
ni1d(k) = MIN(250.D3, cig(1)*oig2*qi1d(k)/am_i*lami**bm_i)
endif
endif
qr1d(k) = qr1d(k) + qrten(k)*DT
if (qr1d(k) .le. R1) qr1d(k) = 0.0
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
!+---+-----------------------------------------------------------------+
!
!+---+-----------------------------------------------------------------+
!..Creation of the lookup tables and support functions found below here.
!+---+-----------------------------------------------------------------+
!..Rain collecting graupel (and inverse). Explicit CE integration.
!+---+-----------------------------------------------------------------+
subroutine qr_acr_qg (docs) 2
implicit none
!..Local variables
INTEGER:: i, j, k, n, n2
DOUBLE PRECISION, DIMENSION(nbg):: vg, N_g
DOUBLE PRECISION, DIMENSION(nbr):: vr, N_r
DOUBLE PRECISION:: N0_exp, N0_r, N0_g, lam_exp, lamg, lamr, N0_s
DOUBLE PRECISION:: massg, massr, dvg, dvr, t1, t2, z1, z2
!+---+
do n2 = 1, nbr
vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
enddo
do n = 1, nbg
vg(n) = av_g*Dg(n)**bv_g
enddo
do k = 1, ntb_r
do j = 1, ntb_r1
lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(k))**ore1
lamr = lam_exp * (crg(3)*org2*org1)**obmr
N0_r = N0r_exp(j)/(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 i = 1, ntb_g
N0_exp = 200.0d0/r_g(i)
N0_exp = DMAX1(gonv_min*1.d0,DMIN1(N0_exp,gonv_max*1.d0))
lam_exp = (N0_exp*am_g*cgg(1)/r_g(i))**oge1
lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
N0_g = N0_exp/(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
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)
t2 = t2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
*dvr*massr * 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) = t1
tmr_racg(i,j,k) = DMIN1(z1, r_r(k)*1.0d0)
tcr_gacr(i,j,k) = t2
tmg_gacr(i,j,k) = z2
enddo
enddo
enddo
end subroutine qr_acr_qg
!+---+-----------------------------------------------------------------+
!
!+---+-----------------------------------------------------------------+
!..Rain collecting snow (and inverse). Explicit CE integration.
!+---+-----------------------------------------------------------------+
subroutine qr_acr_qs (docs) 2
implicit none
!..Local variables
INTEGER:: i, j, k, m, n, n2
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
!+---+
do n2 = 1, nbr
vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*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
do m = 1, ntb_r
do k = 1, ntb_r1
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
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. 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)
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)
endif
if (massr .gt. masss) then
t2 = t2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
*dvr*massr * 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)
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
enddo
enddo
enddo
enddo
end subroutine qr_acr_qs
!+---+-----------------------------------------------------------------+
!
!+---+-----------------------------------------------------------------+
!..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 (docs) 2
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
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
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
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
!+---+-----------------------------------------------------------------+
!
!+---+-----------------------------------------------------------------+
!..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 (docs) 2,2
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
tpi_ide(i,j) = GAMMP(mu_i+2.0, SNGL(lami)*D0s) * 1.0D0
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
!
!+---+-----------------------------------------------------------------+
!..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 (docs) 2
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
!+---+-----------------------------------------------------------------+
!..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 (docs) 2
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
!+---+-----------------------------------------------------------------+
!+---+-----------------------------------------------------------------+
SUBROUTINE GCF (docs) (GAMMCF,A,X,GLN) 2,2
! --- 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 (docs) (GAMSER,A,X,GLN) 2,2
! --- 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 (docs) (XX) 4
! --- 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 (docs) (A,X) 3,4
! --- 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 (docs) (y) 24
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 (docs) (P,T) 6
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)
END FUNCTION RSLF
!
!+---+-----------------------------------------------------------------+
! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A
! FUNCTION OF TEMPERATURE AND PRESSURE
!
REAL FUNCTION RSIF (docs) (P,T) 2
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)
END FUNCTION RSIF
!+---+-----------------------------------------------------------------+
END MODULE module_mp_thompson07
!+---+-----------------------------------------------------------------+
!
! MODIFICATIONS TO MAKE IN OTHER MODULES (pre v2.2 code only)
!
! Use this new code by changing the "THOMPSON" section of code found
! in "module_microphysics_driver.F" with this section. [Of course
! remove the leading comment character that you see here.]
!
! CASE (THOMPSON)
! CALL wrf_debug ( 100 , 'microphysics_driver: calling thompson et al' )
! IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. &
! PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. &
! PRESENT( QS_CURR ) .AND. PRESENT ( QG_CURR ) .AND. &
! PRESENT ( QNI_CURR ).AND. &
! PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) ) THEN
! CALL mp_gt_driver( &
! QV=qv_curr, &
! QC=qc_curr, &
! QR=qr_curr, &
! QI=qi_curr, &
! QS=qs_curr, &
! QG=qg_curr, &
! NI=qni_curr, &
! TH=th, &
! PII=pi_phy, &
! P=p, &
! DZ=dz8w, &
! DT_IN=dt, &
! ITIMESTEP=itimestep, &
! RAINNC=RAINNC, &
! RAINNCV=RAINNCV, &
! SR=SR &
! ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
! ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
! ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte)
! ELSE
! CALL wrf_error_fatal ( 'arguments not present for calling thompson_et_al' )
! ENDIF
!
! Then rename the call from "thomp_init" to "thompson_init" in the file
! "module_physics_init.F" (seen below):
!
! CASE (THOMPSON)
! CALL thompson_init