1 !+---+-----------------------------------------------------------------+
2 !.. This subroutine computes the moisture tendencies of water vapor,
3 !.. cloud droplets, rain, cloud ice (pristine), snow, and graupel.
4 !.. Prior to WRFv2.2 this code was based on Reisner et al (1998), but
5 !.. few of those pieces remain. A complete description is now found
6 !.. in Thompson et al. (2004, 2007).
7 !.. Most importantly, users may wish to modify the prescribed number of
8 !.. cloud droplets (Nt_c; see guidelines mentioned below). Otherwise,
9 !.. users may alter the rain and graupel size distribution parameters
10 !.. to use exponential (Marshal-Palmer) or generalized gamma shape.
11 !.. The snow field assumes a combination of two gamma functions (from
12 !.. Field et al. 2005) and would require significant modifications
13 !.. throughout the entire code to alter its shape as well as accretion
14 !.. rates. Users may also alter the constants used for density of rain,
15 !.. graupel, ice, and snow, but the latter is not constant when using
16 !.. Paul Field's snow distribution and moments methods. Other values
17 !.. users can modify include the constants for mass and/or velocity
18 !.. power law relations and assumed capacitances used in deposition/
19 !.. sublimation/evaporation/melting.
20 !.. Remaining values should probably be left alone.
22 !..Author: Greg Thompson, NCAR-RAL, gthompsn@ucar.edu, 303-497-2805
23 !..Last modified: 26 Oct 2007
24 !+---+-----------------------------------------------------------------+
25 !wrft:model_layer:physics
26 !+---+-----------------------------------------------------------------+
28 MODULE module_mp_thompson07
33 LOGICAL, PARAMETER, PRIVATE:: iiwarm = .false.
34 INTEGER, PARAMETER, PRIVATE:: IFDRY = 0
35 REAL, PARAMETER, PRIVATE:: T_0 = 273.15
36 REAL, PARAMETER, PRIVATE:: PI = 3.1415926536
38 !..Densities of rain, snow, graupel, and cloud ice.
39 REAL, PARAMETER, PRIVATE:: rho_w = 1000.0
40 REAL, PARAMETER, PRIVATE:: rho_s = 100.0
41 REAL, PARAMETER, PRIVATE:: rho_g = 400.0
42 REAL, PARAMETER, PRIVATE:: rho_i = 890.0
44 !..Prescribed number of cloud droplets. Set according to known data or
45 !.. roughly 100 per cc (100.E6 m^-3) for Maritime cases and
46 !.. 300 per cc (300.E6 m^-3) for Continental. Gamma shape parameter,
47 !.. mu_c, calculated based on Nt_c is important in autoconversion
49 REAL, PARAMETER, PRIVATE:: Nt_c = 100.E6
51 !..Generalized gamma distributions for rain, graupel and cloud ice.
52 !.. N(D) = N_0 * D**mu * exp(-lamda*D); mu=0 is exponential.
53 REAL, PARAMETER, PRIVATE:: mu_r = 0.0
54 REAL, PARAMETER, PRIVATE:: mu_g = 0.0
55 REAL, PARAMETER, PRIVATE:: mu_i = 0.0
58 !..Sum of two gamma distrib for snow (Field et al. 2005).
59 !.. N(D) = M2**4/M3**3 * [Kap0*exp(-M2*Lam0*D/M3)
60 !.. + Kap1*(M2/M3)**mu_s * D**mu_s * exp(-M2*Lam1*D/M3)]
61 !.. M2 and M3 are the (bm_s)th and (bm_s+1)th moments respectively
62 !.. calculated as function of ice water content and temperature.
63 REAL, PARAMETER, PRIVATE:: mu_s = 0.6357
64 REAL, PARAMETER, PRIVATE:: Kap0 = 490.6
65 REAL, PARAMETER, PRIVATE:: Kap1 = 17.46
66 REAL, PARAMETER, PRIVATE:: Lam0 = 20.78
67 REAL, PARAMETER, PRIVATE:: Lam1 = 3.29
69 !..Y-intercept parameters for rain & graupel. However, these are not
70 !.. constant and vary depending on mixing ratio. Furthermore, when
71 !.. mu is non-zero, these become equiv y-intercept for an exponential
72 !.. distrib and proper values computed based on assumed mu value.
73 REAL, PARAMETER, PRIVATE:: gonv_min = 1.E4
74 REAL, PARAMETER, PRIVATE:: gonv_max = 5.E6
75 REAL, PARAMETER, PRIVATE:: ronv_min = 2.E6
76 REAL, PARAMETER, PRIVATE:: ronv_max = 9.E9
77 REAL, PARAMETER, PRIVATE:: ronv_sl = 1./4.
78 REAL, PARAMETER, PRIVATE:: ronv_r0 = 0.10E-3
79 REAL, PARAMETER, PRIVATE:: ronv_c0 = ronv_sl/ronv_r0
80 REAL, PARAMETER, PRIVATE:: ronv_c1 = (ronv_max-ronv_min)*0.5
81 REAL, PARAMETER, PRIVATE:: ronv_c2 = (ronv_max+ronv_min)*0.5
83 !..Mass power law relations: mass = am*D**bm
84 !.. Snow from Field et al. (2005), others assume spherical form.
85 REAL, PARAMETER, PRIVATE:: am_r = PI*rho_w/6.0
86 REAL, PARAMETER, PRIVATE:: bm_r = 3.0
87 REAL, PARAMETER, PRIVATE:: am_s = 0.069
88 REAL, PARAMETER, PRIVATE:: bm_s = 2.0
89 REAL, PARAMETER, PRIVATE:: am_g = PI*rho_g/6.0
90 REAL, PARAMETER, PRIVATE:: bm_g = 3.0
91 REAL, PARAMETER, PRIVATE:: am_i = PI*rho_i/6.0
92 REAL, PARAMETER, PRIVATE:: bm_i = 3.0
94 !..Fallspeed power laws relations: v = (av*D**bv)*exp(-fv*D)
95 !.. Rain from Ferrier (1994), ice, snow, and graupel from
96 !.. Thompson et al (2006). Coefficient fv is zero for graupel/ice.
97 REAL, PARAMETER, PRIVATE:: av_r = 4854.0
98 REAL, PARAMETER, PRIVATE:: bv_r = 1.0
99 REAL, PARAMETER, PRIVATE:: fv_r = 195.0
100 REAL, PARAMETER, PRIVATE:: av_s = 40.0
101 REAL, PARAMETER, PRIVATE:: bv_s = 0.55
102 REAL, PARAMETER, PRIVATE:: fv_s = 125.0
103 REAL, PARAMETER, PRIVATE:: av_g = 442.0
104 REAL, PARAMETER, PRIVATE:: bv_g = 0.89
105 REAL, PARAMETER, PRIVATE:: av_i = 1847.5
106 REAL, PARAMETER, PRIVATE:: bv_i = 1.0
108 !..Capacitance of sphere and plates/aggregates: D**3, D**2
109 REAL, PARAMETER, PRIVATE:: C_cube = 0.5
110 REAL, PARAMETER, PRIVATE:: C_sqrd = 0.3
112 !..Collection efficiencies. Rain/snow/graupel collection of cloud
113 !.. droplets use variables (Ef_rw, Ef_sw, Ef_gw respectively) and
114 !.. get computed elsewhere because they are dependent on stokes
116 REAL, PARAMETER, PRIVATE:: Ef_si = 0.1
117 REAL, PARAMETER, PRIVATE:: Ef_rs = 0.95
118 REAL, PARAMETER, PRIVATE:: Ef_rg = 0.95
119 REAL, PARAMETER, PRIVATE:: Ef_ri = 0.95
121 !..Minimum microphys values
122 !.. R1 value, 1.E-12, cannot be set lower because of numerical
123 !.. problems with Paul Field's moments and should not be set larger
124 !.. because of truncation problems in snow/ice growth.
125 REAL, PARAMETER, PRIVATE:: R1 = 1.E-12
126 REAL, PARAMETER, PRIVATE:: R2 = 1.E-8
127 REAL, PARAMETER, PRIVATE:: eps = 1.E-29
129 !..Constants in Cooper curve relation for cloud ice number.
130 REAL, PARAMETER, PRIVATE:: TNO = 5.0
131 REAL, PARAMETER, PRIVATE:: ATO = 0.304
133 !..Rho_not used in fallspeed relations (rho_not/rho)**.5 adjustment.
134 REAL, PARAMETER, PRIVATE:: rho_not = 101325.0/(287.05*298.0)
137 REAL, PARAMETER, PRIVATE:: Sc = 0.632
140 !..Homogeneous freezing temperature
141 REAL, PARAMETER, PRIVATE:: HGFR = 235.16
143 !..Water vapor and air gas constants at constant pressure
144 REAL, PARAMETER, PRIVATE:: Rv = 461.5
145 REAL, PARAMETER, PRIVATE:: oRv = 1./Rv
146 REAL, PARAMETER, PRIVATE:: R = 287.04
147 REAL, PARAMETER, PRIVATE:: Cp = 1004.0
149 !..Enthalpy of sublimation, vaporization, and fusion at 0C.
150 REAL, PARAMETER, PRIVATE:: lsub = 2.834E6
151 REAL, PARAMETER, PRIVATE:: lvap0 = 2.5E6
152 REAL, PARAMETER, PRIVATE:: lfus = lsub - lvap0
153 REAL, PARAMETER, PRIVATE:: olfus = 1./lfus
155 !..Ice initiates with this mass (kg), corresponding diameter calc.
156 !..Min diameters and mass of cloud, rain, snow, and graupel (m, kg).
157 REAL, PARAMETER, PRIVATE:: xm0i = 1.E-12
158 REAL, PARAMETER, PRIVATE:: D0c = 1.E-6
159 REAL, PARAMETER, PRIVATE:: D0r = 50.E-6
160 REAL, PARAMETER, PRIVATE:: D0s = 200.E-6
161 REAL, PARAMETER, PRIVATE:: D0g = 250.E-6
162 REAL, PRIVATE:: D0i, xm0s, xm0g
164 !..Lookup table dimensions
165 INTEGER, PARAMETER, PRIVATE:: nbins = 100
166 INTEGER, PARAMETER, PRIVATE:: nbc = nbins
167 INTEGER, PARAMETER, PRIVATE:: nbi = nbins
168 INTEGER, PARAMETER, PRIVATE:: nbr = nbins
169 INTEGER, PARAMETER, PRIVATE:: nbs = nbins
170 INTEGER, PARAMETER, PRIVATE:: nbg = nbins
171 INTEGER, PARAMETER, PRIVATE:: ntb_c = 37
172 INTEGER, PARAMETER, PRIVATE:: ntb_i = 64
173 INTEGER, PARAMETER, PRIVATE:: ntb_r = 37
174 INTEGER, PARAMETER, PRIVATE:: ntb_s = 28
175 INTEGER, PARAMETER, PRIVATE:: ntb_g = 28
176 INTEGER, PARAMETER, PRIVATE:: ntb_r1 = 37
177 INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55
178 INTEGER, PARAMETER, PRIVATE:: ntb_t = 9
179 INTEGER, PRIVATE:: nic2, nii2, nii3, nir2, nir3, nis2, nig2
181 DOUBLE PRECISION, DIMENSION(nbins+1):: xDx
182 DOUBLE PRECISION, DIMENSION(nbc):: Dc, dtc
183 DOUBLE PRECISION, DIMENSION(nbi):: Di, dti
184 DOUBLE PRECISION, DIMENSION(nbr):: Dr, dtr
185 DOUBLE PRECISION, DIMENSION(nbs):: Ds, dts
186 DOUBLE PRECISION, DIMENSION(nbg):: Dg, dtg
188 !..Lookup tables for cloud water content (kg/m**3).
189 REAL, DIMENSION(ntb_c), PARAMETER, PRIVATE:: &
190 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, &
191 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, &
192 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, &
193 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, &
196 !..Lookup tables for cloud ice content (kg/m**3).
197 REAL, DIMENSION(ntb_i), PARAMETER, PRIVATE:: &
198 r_i = (/1.e-10,2.e-10,3.e-10,4.e-10, &
199 5.e-10,6.e-10,7.e-10,8.e-10,9.e-10, &
200 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, &
201 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, &
202 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, &
203 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, &
204 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, &
205 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, &
208 !..Lookup tables for rain content (kg/m**3).
209 REAL, DIMENSION(ntb_r), PARAMETER, PRIVATE:: &
210 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, &
211 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, &
212 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, &
213 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, &
216 !..Lookup tables for graupel content (kg/m**3).
217 REAL, DIMENSION(ntb_g), PARAMETER, PRIVATE:: &
218 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, &
219 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, &
220 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, &
223 !..Lookup tables for snow content (kg/m**3).
224 REAL, DIMENSION(ntb_s), PARAMETER, PRIVATE:: &
225 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, &
226 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, &
227 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, &
230 !..Lookup tables for rain y-intercept parameter (/m**4).
231 REAL, DIMENSION(ntb_r1), PARAMETER, PRIVATE:: &
232 N0r_exp = (/1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, &
233 1.e7,2.e7,3.e7,4.e7,5.e7,6.e7,7.e7,8.e7,9.e7, &
234 1.e8,2.e8,3.e8,4.e8,5.e8,6.e8,7.e8,8.e8,9.e8, &
235 1.e9,2.e9,3.e9,4.e9,5.e9,6.e9,7.e9,8.e9,9.e9, &
238 !..Lookup tables for ice number concentration (/m**3).
239 REAL, DIMENSION(ntb_i1), PARAMETER, PRIVATE:: &
240 Nt_i = (/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0, &
241 1.e1,2.e1,3.e1,4.e1,5.e1,6.e1,7.e1,8.e1,9.e1, &
242 1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, &
243 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, &
244 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
245 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
248 !..For snow moments conversions (from Field et al. 2005)
249 REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
250 sa = (/ 5.065339, -0.062659, -3.032362, 0.029469, -0.000285, &
251 0.31255, 0.000204, 0.003199, 0.0, -0.015952/)
252 REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
253 sb = (/ 0.476221, -0.015896, 0.165977, 0.007468, -0.000141, &
254 0.060366, 0.000079, 0.000594, 0.0, -0.003577/)
256 !..Temperatures (5 C interval 0 to -40) used in lookup tables.
257 REAL, DIMENSION(ntb_t), PARAMETER, PRIVATE:: &
258 Tc = (/-0.01, -5., -10., -15., -20., -25., -30., -35., -40./)
260 !..Lookup tables for various accretion/collection terms.
261 !.. ntb_x refers to the number of elements for rain, snow, graupel,
262 !.. and temperature array indices. Variables beginning with tp/tc/tm
263 !.. represent lookup tables.
264 ! These are allocatable, 20090612, JM
265 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:):: &
266 tcg_racg, tmr_racg, tcr_gacr, tmg_gacr
267 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:,:):: &
268 tcs_racs1, tmr_racs1, tcs_racs2, tmr_racs2, &
269 tcr_sacr1, tms_sacr1, tcr_sacr2, tms_sacr2
270 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:):: &
272 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:):: &
273 tpi_qrfz, tpg_qrfz, tni_qrfz
274 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:):: &
275 tps_iaus, tni_iaus, tpi_ide
276 REAL, ALLOCATABLE, DIMENSION(:,:):: t_Efrw
277 REAL, ALLOCATABLE, DIMENSION(:,:):: t_Efsw
279 !..Variables holding a bunch of exponents and gamma values (cloud water,
280 !.. cloud ice, rain, snow, then graupel).
281 REAL, DIMENSION(3), PRIVATE:: cce, ccg
282 REAL, PRIVATE:: ocg1, ocg2
283 REAL, DIMENSION(6), PRIVATE:: cie, cig
284 REAL, PRIVATE:: oig1, oig2, obmi
285 REAL, DIMENSION(12), PRIVATE:: cre, crg
286 REAL, PRIVATE:: ore1, org1, org2, org3, obmr
287 REAL, DIMENSION(18), PRIVATE:: cse, csg
288 REAL, PRIVATE:: oams, obms, ocms
289 REAL, DIMENSION(12), PRIVATE:: cge, cgg
290 REAL, PRIVATE:: oge1, ogg1, ogg2, ogg3, oamg, obmg, ocmg
292 !..Declaration of precomputed constants in various rate eqns.
293 REAL:: t1_qr_qc, t1_qr_qi, t2_qr_qi, t1_qg_qc, t1_qs_qc, t1_qs_qi
294 REAL:: t1_qr_ev, t2_qr_ev
295 REAL:: t1_qs_sd, t2_qs_sd, t1_qg_sd, t2_qg_sd
296 REAL:: t1_qs_me, t2_qs_me, t1_qg_me, t2_qg_me
298 CHARACTER*256:: mp_debug
301 !+---+-----------------------------------------------------------------+
303 !+---+-----------------------------------------------------------------+
309 SUBROUTINE thompson07_init
313 INTEGER:: i, j, k, m, n
316 ! jm allocate lookup tables
322 !jm allocate the lookup tables
324 ! use whether the first lookup table has been allocated to also determine whether
325 ! to initialize them all.
328 IF ( .NOT. ALLOCATED( tcg_racg ) ) THEN
329 ALLOCATE( tcg_racg(ntb_g,ntb_r1,ntb_r) )
332 IF ( .NOT. ALLOCATED( tmr_racg) ) ALLOCATE( tmr_racg(ntb_g,ntb_r1,ntb_r) )
333 IF ( .NOT. ALLOCATED( tcr_gacr) ) ALLOCATE( tcr_gacr(ntb_g,ntb_r1,ntb_r) )
334 IF ( .NOT. ALLOCATED( tmg_gacr) ) ALLOCATE( tmg_gacr(ntb_g,ntb_r1,ntb_r) )
336 IF ( .NOT. ALLOCATED( tcs_racs1) ) ALLOCATE( tcs_racs1( ntb_s,ntb_t,ntb_r1,ntb_r ) )
337 IF ( .NOT. ALLOCATED( tmr_racs1) ) ALLOCATE( tmr_racs1( ntb_s,ntb_t,ntb_r1,ntb_r ) )
338 IF ( .NOT. ALLOCATED( tcs_racs2) ) ALLOCATE( tcs_racs2( ntb_s,ntb_t,ntb_r1,ntb_r ) )
339 IF ( .NOT. ALLOCATED( tmr_racs2) ) ALLOCATE( tmr_racs2( ntb_s,ntb_t,ntb_r1,ntb_r ) )
341 IF ( .NOT. ALLOCATED( tcr_sacr1) ) ALLOCATE( tcr_sacr1( ntb_s,ntb_t,ntb_r1,ntb_r ) )
342 IF ( .NOT. ALLOCATED( tms_sacr1) ) ALLOCATE( tms_sacr1( ntb_s,ntb_t,ntb_r1,ntb_r ) )
343 IF ( .NOT. ALLOCATED( tcr_sacr2) ) ALLOCATE( tcr_sacr2( ntb_s,ntb_t,ntb_r1,ntb_r ) )
344 IF ( .NOT. ALLOCATED( tms_sacr2) ) ALLOCATE( tms_sacr2( ntb_s,ntb_t,ntb_r1,ntb_r ) )
346 IF ( .NOT. ALLOCATED( tpi_qcfz) ) ALLOCATE( tpi_qcfz( ntb_c,45 ) )
347 IF ( .NOT. ALLOCATED( tni_qcfz) ) ALLOCATE( tni_qcfz( ntb_c,45 ) )
349 IF ( .NOT. ALLOCATED( tpi_qrfz) ) ALLOCATE( tpi_qrfz( ntb_r,ntb_r1,45 ) )
350 IF ( .NOT. ALLOCATED( tpg_qrfz) ) ALLOCATE( tpg_qrfz( ntb_r,ntb_r1,45 ) )
351 IF ( .NOT. ALLOCATED( tni_qrfz) ) ALLOCATE( tni_qrfz( ntb_r,ntb_r1,45 ) )
353 IF ( .NOT. ALLOCATED( tps_iaus) ) ALLOCATE( tps_iaus(ntb_i,ntb_i1) )
354 IF ( .NOT. ALLOCATED( tni_iaus) ) ALLOCATE( tni_iaus(ntb_i,ntb_i1) )
355 IF ( .NOT. ALLOCATED( tpi_ide) ) ALLOCATE( tpi_ide(ntb_i,ntb_i1) )
357 IF ( .NOT. ALLOCATED( t_Efrw) ) ALLOCATE( t_Efrw(nbr,nbc) )
358 IF ( .NOT. ALLOCATED( t_Efsw) ) ALLOCATE( t_Efsw(nbs,nbc) )
366 !..From Martin et al. (1994), assign gamma shape parameter mu for cloud
367 !.. drops according to general dispersion characteristics (disp=~0.25
368 !.. for Maritime and 0.45 for Continental).
369 !.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime
370 !.. to 2 for really dirty air.
371 mu_c = MIN(15., (1000.E6/Nt_c + 2.))
373 !..Schmidt number to one-third used numerous times.
376 !..Compute min ice diam from mass, min snow/graupel mass from diam.
377 D0i = (xm0i/am_i)**(1./bm_i)
378 xm0s = am_s * D0s**bm_s
379 xm0g = am_g * D0g**bm_g
381 !..These constants various exponents and gamma() assoc with cloud,
382 !.. rain, snow, and graupel.
384 cce(2) = bm_r + mu_c + 1.
385 cce(3) = bm_r + mu_c + 4.
386 ccg(1) = WGAMMA(cce(1))
387 ccg(2) = WGAMMA(cce(2))
388 ccg(3) = WGAMMA(cce(3))
393 cie(2) = bm_i + mu_i + 1.
394 cie(3) = bm_i + mu_i + bv_i + 1.
395 cie(4) = mu_i + bv_i + 1.
398 cig(1) = WGAMMA(cie(1))
399 cig(2) = WGAMMA(cie(2))
400 cig(3) = WGAMMA(cie(3))
401 cig(4) = WGAMMA(cie(4))
402 cig(5) = WGAMMA(cie(5))
403 cig(6) = WGAMMA(cie(6))
410 cre(3) = bm_r + mu_r + 1.
411 cre(4) = bm_r*2. + mu_r + 1.
412 cre(5) = bm_r + mu_r + 3.
413 cre(6) = bm_r + mu_r + bv_r + 1.
414 cre(7) = bm_r + mu_r + bv_r + 2.
415 cre(8) = bm_r + mu_r + bv_r + 3.
416 cre(9) = mu_r + bv_r + 3.
418 cre(11) = 0.5*(bv_r + 5. + 2.*mu_r)
419 cre(12) = bm_r + mu_r + 4.
421 crg(n) = WGAMMA(cre(n))
432 cse(4) = bm_s + bv_s + 1.
433 cse(5) = bm_s + bv_s + 2.
434 cse(6) = bm_s + bv_s + 3.
435 cse(7) = bm_s + mu_s + 1.
436 cse(8) = bm_s + mu_s + 2.
437 cse(9) = bm_s + mu_s + 3.
438 cse(10) = bm_s + mu_s + bv_s + 1.
439 cse(11) = bm_s + mu_s + bv_s + 2.
440 cse(12) = bm_s*2. + mu_s + 1.
442 cse(14) = bm_s + bv_s
444 cse(16) = 1.0 + (1.0 + bv_s)/2.
445 cse(17) = cse(16) + mu_s + 1.
446 cse(18) = bv_s + mu_s + 3.
448 csg(n) = WGAMMA(cse(n))
456 cge(3) = bm_g + mu_g + 1.
457 cge(4) = bm_g*2. + mu_g + 1.
458 cge(5) = bm_g + mu_g + 3.
459 cge(6) = bm_g + mu_g + bv_g + 1.
460 cge(7) = bm_g + mu_g + bv_g + 2.
461 cge(8) = bm_g + mu_g + bv_g + 3.
462 cge(9) = mu_g + bv_g + 3.
464 cge(11) = 0.5*(bv_g + 5. + 2.*mu_g)
465 cge(12) = 0.5*(bv_g + 5.) + mu_g
467 cgg(n) = WGAMMA(cge(n))
477 !+---+-----------------------------------------------------------------+
478 !..Simplify various rate eqns the best we can now.
479 !+---+-----------------------------------------------------------------+
481 !..Rain collecting cloud water and cloud ice
482 t1_qr_qc = PI*.25*av_r * crg(9)
483 t1_qr_qi = PI*.25*av_r * crg(9)
484 t2_qr_qi = PI*.25*am_r*av_r * crg(8)
486 !..Graupel collecting cloud water
487 t1_qg_qc = PI*.25*av_g * cgg(9)
489 !..Snow collecting cloud water
490 t1_qs_qc = PI*.25*av_s
492 !..Snow collecting cloud ice
493 t1_qs_qi = PI*.25*av_s
495 !..Evaporation of rain; ignore depositional growth of rain.
496 t1_qr_ev = 0.78 * crg(10)
497 t2_qr_ev = 0.308*Sc3*SQRT(av_r) * crg(11)
499 !..Sublimation/depositional growth of snow
501 t2_qs_sd = 0.28*Sc3*SQRT(av_s)
504 t1_qs_me = PI*4.*C_sqrd*olfus * 0.86
505 t2_qs_me = PI*4.*C_sqrd*olfus * 0.28*Sc3*SQRT(av_s)
507 !..Sublimation/depositional growth of graupel
508 t1_qg_sd = 0.86 * cgg(10)
509 t2_qg_sd = 0.28*Sc3*SQRT(av_g) * cgg(11)
511 !..Melting of graupel
512 t1_qg_me = PI*4.*C_cube*olfus * 0.86 * cgg(10)
513 t2_qg_me = PI*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g) * cgg(11)
515 !..Constants for helping find lookup table indexes.
516 nic2 = NINT(ALOG10(r_c(1)))
517 nii2 = NINT(ALOG10(r_i(1)))
518 nii3 = NINT(ALOG10(Nt_i(1)))
519 nir2 = NINT(ALOG10(r_r(1)))
520 nir3 = NINT(ALOG10(N0r_exp(1)))
521 nis2 = NINT(ALOG10(r_s(1)))
522 nig2 = NINT(ALOG10(r_g(1)))
524 !..Create bins of cloud water (from min diameter up to 100 microns).
528 Dc(n) = Dc(n-1) + 1.0D-6
529 dtc(n) = (Dc(n) - Dc(n-1))
532 !..Create bins of cloud ice (from min diameter up to 5x min snow size).
534 xDx(nbi+1) = 5.0d0*D0s
536 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbi) &
537 *DLOG(xDx(nbi+1)/xDx(1)) +DLOG(xDx(1)))
540 Di(n) = DSQRT(xDx(n)*xDx(n+1))
541 dti(n) = xDx(n+1) - xDx(n)
544 !..Create bins of rain (from min diameter up to 5 mm).
548 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbr) &
549 *DLOG(xDx(nbr+1)/xDx(1)) +DLOG(xDx(1)))
552 Dr(n) = DSQRT(xDx(n)*xDx(n+1))
553 dtr(n) = xDx(n+1) - xDx(n)
556 !..Create bins of snow (from min diameter up to 2 cm).
560 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbs) &
561 *DLOG(xDx(nbs+1)/xDx(1)) +DLOG(xDx(1)))
564 Ds(n) = DSQRT(xDx(n)*xDx(n+1))
565 dts(n) = xDx(n+1) - xDx(n)
568 !..Create bins of graupel (from min diameter up to 5 cm).
572 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbg) &
573 *DLOG(xDx(nbg+1)/xDx(1)) +DLOG(xDx(1)))
576 Dg(n) = DSQRT(xDx(n)*xDx(n+1))
577 dtg(n) = xDx(n+1) - xDx(n)
580 !+---+-----------------------------------------------------------------+
581 !..Create lookup tables for most costly calculations.
582 !+---+-----------------------------------------------------------------+
587 tcg_racg(i,j,k) = 0.0d0
588 tmr_racg(i,j,k) = 0.0d0
589 tcr_gacr(i,j,k) = 0.0d0
590 tmg_gacr(i,j,k) = 0.0d0
599 tcs_racs1(i,j,k,m) = 0.0d0
600 tmr_racs1(i,j,k,m) = 0.0d0
601 tcs_racs2(i,j,k,m) = 0.0d0
602 tmr_racs2(i,j,k,m) = 0.0d0
603 tcr_sacr1(i,j,k,m) = 0.0d0
604 tms_sacr1(i,j,k,m) = 0.0d0
605 tcr_sacr2(i,j,k,m) = 0.0d0
606 tms_sacr2(i,j,k,m) = 0.0d0
615 tpi_qrfz(i,j,k) = 0.0d0
616 tni_qrfz(i,j,k) = 0.0d0
617 tpg_qrfz(i,j,k) = 0.0d0
621 tpi_qcfz(i,k) = 0.0d0
622 tni_qcfz(i,k) = 0.0d0
628 tps_iaus(i,j) = 0.0d0
629 tni_iaus(i,j) = 0.0d0
643 CALL wrf_debug(150, 'CREATING MICROPHYSICS LOOKUP TABLES ... ')
644 WRITE (wrf_err_message, '(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') &
645 ' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g
646 CALL wrf_debug(150, wrf_err_message)
648 !..Collision efficiency between rain/snow and cloud water.
649 CALL wrf_debug(200, ' creating qc collision eff tables')
653 IF (.not. iiwarm) THEN
654 !..Rain collecting graupel & graupel collecting rain.
655 CALL wrf_debug(200, ' creating rain collecting graupel table')
658 !..Rain collecting snow & snow collecting rain.
659 CALL wrf_debug(200, ' creating rain collecting snow table')
662 !..Cloud water and rain freezing (Bigg, 1953).
663 CALL wrf_debug(200, ' creating freezing of water drops table')
666 !..Conversion of some ice mass into snow category.
667 CALL wrf_debug(200, ' creating ice converting to snow table')
672 CALL wrf_debug(150, ' ... DONE microphysical lookup tables')
676 END SUBROUTINE thompson07_init
677 !+---+-----------------------------------------------------------------+
679 !+---+-----------------------------------------------------------------+
680 !..This is a wrapper routine designed to transfer values from 3D to 1D.
681 !+---+-----------------------------------------------------------------+
682 SUBROUTINE mp_gt_driver07(qv, qc, qr, qi, qs, qg, ni, &
683 th, pii, p, dz, dt_in, itimestep, &
684 RAINNC, RAINNCV, SR, &
685 ids,ide, jds,jde, kds,kde, & ! domain dims
686 ims,ime, jms,jme, kms,kme, & ! memory dims
687 its,ite, jts,jte, kts,kte) ! tile dims
691 !..Subroutine arguments
692 INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, &
693 ims,ime, jms,jme, kms,kme, &
694 its,ite, jts,jte, kts,kte
695 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
696 qv, qc, qr, qi, qs, qg, ni, th
697 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
699 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: &
701 REAL, INTENT(IN):: dt_in
702 INTEGER, INTENT(IN):: itimestep
705 REAL, DIMENSION(kts:kte):: &
706 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
708 REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
709 REAL:: dt, pptrain, pptsnow, pptgraul, pptice
710 REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max
712 INTEGER:: imax_qc, imax_qr, imax_qi, imax_qs, imax_qg, imax_ni
713 INTEGER:: jmax_qc, jmax_qr, jmax_qi, jmax_qs, jmax_qg, jmax_ni
714 INTEGER:: kmax_qc, kmax_qr, kmax_qi, kmax_qs, kmax_qg, kmax_ni
715 INTEGER:: i_start, j_start, i_end, j_end
723 if ( (ite-its+1).gt.4 .and. (jte-jts+1).lt.4) then
728 elseif ( (ite-its+1).lt.4 .and. (jte-jts+1).gt.4) then
762 mp_debug(i:i) = char(0)
765 j_loop: do j = j_start, j_end
766 i_loop: do i = i_start, i_end
776 t1d(k) = th(i,k,j)*pii(i,k,j)
788 call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
790 pptrain, pptsnow, pptgraul, pptice, &
793 pcp_ra(i,j) = pptrain
794 pcp_sn(i,j) = pptsnow
795 pcp_gr(i,j) = pptgraul
797 RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice
798 RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
799 SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12)
809 th(i,k,j) = t1d(k)/pii(i,k,j)
810 if (qc1d(k) .gt. qc_max) then
815 elseif (qc1d(k) .lt. 0.0) then
816 write(mp_debug,*) 'WARNING, negative qc ', qc1d(k), &
818 CALL wrf_debug(150, mp_debug)
820 if (qr1d(k) .gt. qr_max) then
825 elseif (qr1d(k) .lt. 0.0) then
826 write(mp_debug,*) 'WARNING, negative qr ', qr1d(k), &
828 CALL wrf_debug(150, mp_debug)
830 if (qs1d(k) .gt. qs_max) then
835 elseif (qs1d(k) .lt. 0.0) then
836 write(mp_debug,*) 'WARNING, negative qs ', qs1d(k), &
838 CALL wrf_debug(150, mp_debug)
840 if (qi1d(k) .gt. qi_max) then
845 elseif (qi1d(k) .lt. 0.0) then
846 write(mp_debug,*) 'WARNING, negative qi ', qi1d(k), &
848 CALL wrf_debug(150, mp_debug)
850 if (qg1d(k) .gt. qg_max) then
855 elseif (qg1d(k) .lt. 0.0) then
856 write(mp_debug,*) 'WARNING, negative qg ', qg1d(k), &
858 CALL wrf_debug(150, mp_debug)
860 if (ni1d(k) .gt. ni_max) then
865 elseif (ni1d(k) .lt. 0.0) then
866 write(mp_debug,*) 'WARNING, negative ni ', ni1d(k), &
868 CALL wrf_debug(150, mp_debug)
870 if (qv1d(k) .lt. 0.0) then
871 write(mp_debug,*) 'WARNING, negative qv ', qv1d(k), &
873 CALL wrf_debug(150, mp_debug)
881 write(mp_debug,'(a,6(a,e13.6,1x,a,i3,a,i3,a,i3,a,1x))') 'MP-GT:', &
882 'qc: ', qc_max, '(', imax_qc, ',', jmax_qc, ',', kmax_qc, ')', &
883 'qr: ', qr_max, '(', imax_qr, ',', jmax_qr, ',', kmax_qr, ')', &
884 'qi: ', qi_max, '(', imax_qi, ',', jmax_qi, ',', kmax_qi, ')', &
885 'qs: ', qs_max, '(', imax_qs, ',', jmax_qs, ',', kmax_qs, ')', &
886 'qg: ', qg_max, '(', imax_qg, ',', jmax_qg, ',', kmax_qg, ')', &
887 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')'
888 CALL wrf_debug(150, mp_debug)
892 mp_debug(i:i) = char(0)
895 END SUBROUTINE mp_gt_driver07
897 !+---+-----------------------------------------------------------------+
899 !+---+-----------------------------------------------------------------+
900 !+---+-----------------------------------------------------------------+
901 !.. This subroutine computes the moisture tendencies of water vapor,
902 !.. cloud droplets, rain, cloud ice (pristine), snow, and graupel.
903 !.. Previously this code was based on Reisner et al (1998), but few of
904 !.. those pieces remain. A complete description is now found in
905 !.. Thompson et al. (2004, 2006).
906 !+---+-----------------------------------------------------------------+
908 subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
910 pptrain, pptsnow, pptgraul, pptice, &
911 kts, kte, dt, ii, jj)
916 INTEGER, INTENT(IN):: kts, kte, ii, jj
917 REAL, DIMENSION(kts:kte), INTENT(INOUT):: &
918 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
920 REAL, DIMENSION(kts:kte), INTENT(IN):: dzq
921 REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice
922 REAL, INTENT(IN):: dt
925 REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, &
926 qrten, qsten, qgten, niten
928 DOUBLE PRECISION, DIMENSION(kts:kte):: prw_vcd
930 DOUBLE PRECISION, DIMENSION(kts:kte):: prr_wau, prr_rcw, prr_rcs, &
931 prr_rcg, prr_sml, prr_gml, &
934 DOUBLE PRECISION, DIMENSION(kts:kte):: pri_inu, pni_inu, pri_ihm, &
935 pni_ihm, pri_wfz, pni_wfz, &
936 pri_rfz, pni_rfz, pri_ide, &
937 pni_ide, pri_rci, pni_rci, &
940 DOUBLE PRECISION, DIMENSION(kts:kte):: prs_iau, prs_sci, prs_rcs, &
941 prs_scw, prs_sde, prs_ihm, &
944 DOUBLE PRECISION, DIMENSION(kts:kte):: prg_scw, prg_rfz, prg_gde, &
945 prg_gcw, prg_rci, prg_rcs, &
948 REAL, DIMENSION(kts:kte):: temp, pres, qv
949 REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni
950 REAL, DIMENSION(kts:kte):: rho, rhof, rhof2
951 REAL, DIMENSION(kts:kte):: qvs, qvsi
952 REAL, DIMENSION(kts:kte):: satw, sati, ssatw, ssati
953 REAL, DIMENSION(kts:kte):: diffu, visco, vsc2, &
954 tcond, lvap, ocp, lvt2
956 DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g
957 REAL, DIMENSION(kts:kte):: mvd_r, mvd_c
958 REAL, DIMENSION(kts:kte):: smob, smo2, smo1, &
959 smoc, smod, smoe, smof
961 REAL, DIMENSION(kts:kte):: sed_r, sed_s, sed_g, sed_i, sed_n
963 REAL:: rgvm, delta_tp, orho, onstep, lfus2
964 DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamc, lamr, lamg
965 DOUBLE PRECISION:: lami, ilami
966 REAL:: xDc, Dc_b, Dc_g, xDi, xDr, xDs, xDg, Ds_m, Dg_m
967 REAL:: zeta1, zeta, taud, tau
968 REAL:: stoke_r, stoke_s, stoke_g, stoke_i
969 REAL:: vti, vtr, vts, vtg
970 REAL, DIMENSION(kts:kte+1):: vtik, vtnk, vtrk, vtsk, vtgk
971 REAL, DIMENSION(kts:kte):: vts_boost
972 REAL:: Mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, C_snow
973 REAL:: a_, b_, loga_, A1, A2, tf
974 REAL:: tempc, tc0, r_mvd1, r_mvd2, xkrat
975 REAL:: xnc, xri, xni, xmi, oxmi, xrc
976 REAL:: xsat, rate_max, sump, ratio
977 REAL:: clap, fcd, dfcd
978 REAL:: otemp, rvs, rvs_p, rvs_pp, gamsc, alphsc, t1_evap, t1_subl
979 REAL:: r_frac, g_frac
980 REAL:: Ef_rw, Ef_sw, Ef_gw
981 REAL:: dtsave, odts, odt, odzq
982 INTEGER:: i, k, k2, ksed1, ku, n, nn, nstep, k_0, kbot, IT, iexfrq
983 INTEGER:: nir, nis, nig, nii, nic
984 INTEGER:: idx_tc,idx_t,idx_s,idx_g,idx_r1,idx_r,idx_i1,idx_i,idx_c
986 LOGICAL:: melti, no_micro
987 LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg
997 !+---+-----------------------------------------------------------------+
998 !.. Source/sink terms. First 2 chars: "pr" represents source/sink of
999 !.. mass while "pn" represents source/sink of number. Next char is one
1000 !.. of "v" for water vapor, "r" for rain, "i" for cloud ice, "w" for
1001 !.. cloud water, "s" for snow, and "g" for graupel. Next chars
1002 !.. represent processes: "de" for sublimation/deposition, "ev" for
1003 !.. evaporation, "fz" for freezing, "ml" for melting, "au" for
1004 !.. autoconversion, "nu" for ice nucleation, "hm" for Hallet/Mossop
1005 !.. secondary ice production, and "c" for collection followed by the
1006 !.. character for the species being collected. ALL of these terms are
1007 !.. positive (except for deposition/sublimation terms which can switch
1008 !.. signs based on super/subsaturation) and are treated as negatives
1009 !.. where necessary in the tendency equations.
1010 !+---+-----------------------------------------------------------------+
1066 !+---+-----------------------------------------------------------------+
1067 !..Put column of data into local arrays.
1068 !+---+-----------------------------------------------------------------+
1071 qv(k) = MAX(1.E-10, qv1d(k))
1073 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
1074 if (qc1d(k) .gt. R1) then
1076 rc(k) = qc1d(k)*rho(k)
1083 if (qi1d(k) .gt. R1) then
1085 ri(k) = qi1d(k)*rho(k)
1086 ni(k) = MAX(1., ni1d(k)*rho(k))
1088 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
1090 xDi = (bm_i + mu_i + 1.) * ilami
1091 if (xDi.lt. 30.E-6) then
1092 lami = cie(2)/30.E-6
1093 ni(k) = MIN(250.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
1094 elseif (xDi.gt. 300.E-6) then
1095 lami = cie(2)/300.E-6
1096 ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i
1105 if (qr1d(k) .gt. R1) then
1107 rr(k) = qr1d(k)*rho(k)
1114 if (qs1d(k) .gt. R1) then
1116 rs(k) = qs1d(k)*rho(k)
1123 if (qg1d(k) .gt. R1) then
1125 rg(k) = qg1d(k)*rho(k)
1135 !+---+-----------------------------------------------------------------+
1136 !..Derive various thermodynamic variables frequently used.
1137 !.. Saturation vapor pressure (mixing ratio) over liquid/ice comes from
1138 !.. Flatau et al. 1992; enthalpy (latent heat) of vaporization from
1139 !.. Bohren & Albrecht 1998; others from Pruppacher & Klett 1978.
1140 !+---+-----------------------------------------------------------------+
1142 tempc = temp(k) - 273.15
1143 rhof(k) = SQRT(RHO_NOT/rho(k))
1144 rhof2(k) = SQRT(rhof(k))
1145 qvs(k) = rslf(pres(k), temp(k))
1146 if (tempc .le. 0.0) then
1147 qvsi(k) = rsif(pres(k), temp(k))
1151 satw(k) = qv(k)/qvs(k)
1152 sati(k) = qv(k)/qvsi(k)
1153 ssatw(k) = satw(k) - 1.
1154 ssati(k) = sati(k) - 1.
1155 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
1156 if (abs(ssati(k)).lt. eps) ssati(k) = 0.0
1157 if (no_micro .and. ssati(k).gt. 0.0) no_micro = .false.
1158 diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
1159 if (tempc .ge. 0.0) then
1160 visco(k) = (1.718+0.0049*tempc)*1.0E-5
1162 visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
1164 ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
1165 vsc2(k) = SQRT(rho(k)/visco(k))
1166 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
1167 tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
1170 !+---+-----------------------------------------------------------------+
1171 !..If no existing hydrometeor species and no chance to initiate ice or
1172 !.. condense cloud water, just exit quickly!
1173 !+---+-----------------------------------------------------------------+
1175 if (no_micro) return
1177 !+---+-----------------------------------------------------------------+
1178 !..Calculate y-intercept, slope, and useful moments for snow.
1179 !+---+-----------------------------------------------------------------+
1180 if (.not. iiwarm) then
1182 if (.not. L_qs(k)) CYCLE
1183 tc0 = MIN(-0.1, temp(k)-273.15)
1184 smob(k) = rs(k)*oams
1186 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
1187 !.. then we must compute actual 2nd moment and use as reference.
1188 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
1191 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
1192 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
1193 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
1194 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
1195 + sa(10)*bm_s*bm_s*bm_s
1197 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
1198 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
1199 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
1200 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
1201 + sb(10)*bm_s*bm_s*bm_s
1202 smo2(k) = (smob(k)/a_)**(1./b_)
1205 !..Calculate 1st moment. Useful for depositional growth and melting.
1206 loga_ = sa(1) + sa(2)*tc0 + sa(3) &
1207 + sa(4)*tc0 + sa(5)*tc0*tc0 &
1208 + sa(6) + sa(7)*tc0*tc0 &
1209 + sa(8)*tc0 + sa(9)*tc0*tc0*tc0 &
1212 b_ = sb(1)+ sb(2)*tc0 + sb(3) + sb(4)*tc0 &
1213 + sb(5)*tc0*tc0 + sb(6) &
1214 + sb(7)*tc0*tc0 + sb(8)*tc0 &
1215 + sb(9)*tc0*tc0*tc0 + sb(10)
1216 smo1(k) = a_ * smo2(k)**b_
1218 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
1219 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
1220 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
1221 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
1222 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
1223 + sa(10)*cse(1)*cse(1)*cse(1)
1225 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
1226 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
1227 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
1228 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
1229 smoc(k) = a_ * smo2(k)**b_
1231 !..Calculate bv_s+2 (th) moment. Useful for riming.
1232 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(13) &
1233 + sa(4)*tc0*cse(13) + sa(5)*tc0*tc0 &
1234 + sa(6)*cse(13)*cse(13) + sa(7)*tc0*tc0*cse(13) &
1235 + sa(8)*tc0*cse(13)*cse(13) + sa(9)*tc0*tc0*tc0 &
1236 + sa(10)*cse(13)*cse(13)*cse(13)
1238 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(13) + sb(4)*tc0*cse(13) &
1239 + sb(5)*tc0*tc0 + sb(6)*cse(13)*cse(13) &
1240 + sb(7)*tc0*tc0*cse(13) + sb(8)*tc0*cse(13)*cse(13) &
1241 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(13)*cse(13)*cse(13)
1242 smoe(k) = a_ * smo2(k)**b_
1244 !..Calculate 1+(bv_s+1)/2 (th) moment. Useful for depositional growth.
1245 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(16) &
1246 + sa(4)*tc0*cse(16) + sa(5)*tc0*tc0 &
1247 + sa(6)*cse(16)*cse(16) + sa(7)*tc0*tc0*cse(16) &
1248 + sa(8)*tc0*cse(16)*cse(16) + sa(9)*tc0*tc0*tc0 &
1249 + sa(10)*cse(16)*cse(16)*cse(16)
1251 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(16) + sb(4)*tc0*cse(16) &
1252 + sb(5)*tc0*tc0 + sb(6)*cse(16)*cse(16) &
1253 + sb(7)*tc0*tc0*cse(16) + sb(8)*tc0*cse(16)*cse(16) &
1254 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(16)*cse(16)*cse(16)
1255 smof(k) = a_ * smo2(k)**b_
1259 !+---+-----------------------------------------------------------------+
1260 !..Calculate y-intercept, slope values for graupel.
1261 !+---+-----------------------------------------------------------------+
1264 if (.not. L_qg(k)) CYCLE
1265 N0_exp = 200.0*rho(k)/rg(k)
1266 N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max)))
1267 N0_min = MIN(N0_exp, N0_min)
1269 lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
1270 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
1272 N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2)
1277 !+---+-----------------------------------------------------------------+
1278 !..Calculate y-intercept & slope values for rain.
1279 !.. New treatment for variable y-intercept of rain. When rain comes
1280 !.. from melted snow/graupel, compute mass-weighted mean size, melt
1281 !.. into water, compute its mvd and recompute slope/intercept.
1282 !.. If rain not from melted snow, use old relation but hold N0_r
1283 !.. constant at its lowest value. While doing all this, ensure rain
1284 !.. mvd does not exceed reasonable size like 2.5 mm.
1285 !+---+-----------------------------------------------------------------+
1288 ! if (.not. L_qr(k)) CYCLE
1289 N0_exp = ronv_c1*tanh(ronv_c0*(ronv_r0-rr(k))) + ronv_c2
1290 N0_min = MIN(N0_exp, N0_min)
1292 lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
1293 lamr = lam_exp * (crg(3)*org2*org1)**obmr
1294 mvd_r(k) = (3.0+mu_r+0.672) / lamr
1295 if (mvd_r(k) .gt. 2.5e-3) then
1297 lamr = (3.0+mu_r+0.672) / 2.5e-3
1298 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
1299 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
1301 N0_r(k) = N0_exp/(crg(2)*lam_exp) * lamr**cre(2)
1305 if (.not. iiwarm) then
1308 do k = kte-1, kts, -1
1309 if ( (temp(k).gt. T_0) .and. (rr(k).gt. 0.001e-3) &
1310 .and. ((rs(k+1)+rg(k+1)).gt. 0.01e-3) ) then
1319 !.. Locate bottom of melting layer (if any).
1321 do k = k_0-1, kts, -1
1322 if ( (rs(k)+rg(k)).lt. 0.01e-3) goto 136
1327 !.. Compute melted snow/graupel equiv water diameter one K-level above
1328 !.. melting. Set starting rain mvd to either 50 microns or max from
1329 !.. higher up in column.
1331 xDs = smoc(k_0) / smob(k_0)
1332 Ds_m = (am_s*xDs**bm_s / am_r)**obmr
1337 xDg = (bm_g + mu_g + 1.) * ilamg(k_0)
1338 Dg_m = (am_g*xDg**bm_g / am_r)**obmr
1343 r_mvd2 = MIN(MAX(Ds_m, Dg_m, r_mvd1+1.e-6, mvd_r(kbot)), &
1346 !.. Within melting layer, apply linear increase of rain mvd from r_mvd1
1347 !.. to equiv melted snow/graupel value (r_mvd2). So, by the bottom of
1348 !.. the melting layer, the rain will have an mvd that matches that from
1349 !.. melted snow and/or graupel.
1350 if (kbot.gt. 2) then
1351 do k = k_0-1, kbot, -1
1352 if (.not. L_qr(k)) CYCLE
1353 xkrat = REAL(k_0-k)/REAL(k_0-kbot)
1354 mvd_r(k) = MAX(mvd_r(k), xkrat*(r_mvd2-r_mvd1)+r_mvd1)
1355 lamr = (4.0+mu_r) / mvd_r(k)
1356 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
1357 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
1358 N0_exp = MAX(DBLE(ronv_min), MIN(N0_exp, DBLE(ronv_max)))
1359 lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
1360 lamr = lam_exp * (crg(3)*org2*org1)**obmr
1361 mvd_r(k) = (3.0+mu_r+0.672) / lamr
1362 N0_r(k) = rr(k)*lamr**cre(3) / (am_r*crg(3))
1366 !.. Below melting layer, hold N0_r constant unless changes to mixing
1367 !.. ratio increase mvd beyond 2.5 mm threshold, then adjust slope and
1368 !.. intercept to cap mvd at 2.5 mm. In future, we could lower N0_r to
1369 !.. account for self-collection or other sinks.
1370 do k = kbot-1, kts, -1
1371 if (.not. L_qr(k)) CYCLE
1372 N0_r(k) = MIN(N0_r(k), N0_r(kbot))
1373 lamr = (N0_r(k)*am_r*crg(3)/rr(k))**(1./cre(3))
1374 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
1375 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
1376 N0_exp = MAX(DBLE(ronv_min), MIN(N0_exp, DBLE(ronv_max)))
1377 lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
1378 lamr = lam_exp * (crg(3)*org2*org1)**obmr
1379 mvd_r(k) = (3.0+mu_r+0.672) / lamr
1380 if (mvd_r(k) .gt. 2.5e-3) then
1382 lamr = (3.0+mu_r+0.672) / mvd_r(k)
1384 N0_r(k) = rr(k)*lamr**cre(3) / (am_r*crg(3))
1392 !+---+-----------------------------------------------------------------+
1393 !..Compute warm-rain process terms (except evap done later).
1394 !+---+-----------------------------------------------------------------+
1397 if (.not. L_qc(k)) CYCLE
1398 xDc = MAX(D0c*1.E6, ((rc(k)/(am_r*Nt_c))**obmr) * 1.E6)
1399 lamc = (Nt_c*am_r* ccg(2) * ocg1 / rc(k))**obmr
1400 mvd_c(k) = (3.0+mu_c+0.672) / lamc
1402 !..Autoconversion follows Berry & Reinhardt (1974) with characteristic
1403 !.. diameters correctly computed from gamma distrib of cloud droplets.
1404 if (rc(k).gt. 0.01e-3) then
1405 Dc_g = ((ccg(3)*ocg2)**obmr / lamc) * 1.E6
1406 Dc_b = (xDc*xDc*xDc*Dc_g*Dc_g*Dc_g - xDc*xDc*xDc*xDc*xDc*xDc) &
1408 zeta1 = 0.5*((6.25E-6*xDc*Dc_b*Dc_b*Dc_b - 0.4) &
1409 + abs(6.25E-6*xDc*Dc_b*Dc_b*Dc_b - 0.4))
1410 zeta = 0.027*rc(k)*zeta1
1411 taud = 0.5*((0.5*Dc_b - 7.5) + abs(0.5*Dc_b - 7.5)) + R1
1412 tau = 3.72/(rc(k)*taud)
1413 prr_wau(k) = zeta/tau
1414 prr_wau(k) = MIN(DBLE(rc(k)*odts), prr_wau(k))
1417 !..Rain collecting cloud water. In CE, assume Dc<<Dr and vtc=~0.
1418 if (L_qr(k) .and. mvd_r(k).gt. D0r .and. mvd_c(k).gt. D0c) then
1420 idx = 1 + INT(nbr*DLOG(mvd_r(k)/Dr(1))/DLOG(Dr(nbr)/Dr(1)))
1422 Ef_rw = t_Efrw(idx, INT(mvd_c(k)*1.E6))
1423 prr_rcw(k) = rhof(k)*t1_qr_qc*Ef_rw*rc(k)*N0_r(k) &
1424 *((lamr+fv_r)**(-cre(9)))
1425 prr_rcw(k) = MIN(DBLE(rc(k)*odts), prr_rcw(k))
1429 !+---+-----------------------------------------------------------------+
1430 !..Compute all frozen hydrometeor species' process terms.
1431 !+---+-----------------------------------------------------------------+
1432 if (.not. iiwarm) then
1436 !..Temperature lookup table indexes.
1437 tempc = temp(k) - 273.15
1438 idx_tc = MAX(1, MIN(NINT(-tempc), 45) )
1439 idx_t = INT( (tempc-2.5)/5. ) - 1
1440 idx_t = MAX(1, -idx_t)
1441 idx_t = MIN(idx_t, ntb_t)
1442 IT = MAX(1, MIN(NINT(-tempc), 31) )
1444 !..Cloud water lookup table index.
1445 if (rc(k).gt. r_c(1)) then
1446 nic = NINT(ALOG10(rc(k)))
1447 do nn = nic-1, nic+1
1449 if ( (rc(k)/10.**nn).ge.1.0 .and. &
1450 (rc(k)/10.**nn).lt.10.0) goto 141
1453 idx_c = INT(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
1454 idx_c = MAX(1, MIN(idx_c, ntb_c))
1459 !..Cloud ice lookup table indexes.
1460 if (ri(k).gt. r_i(1)) then
1461 nii = NINT(ALOG10(ri(k)))
1462 do nn = nii-1, nii+1
1464 if ( (ri(k)/10.**nn).ge.1.0 .and. &
1465 (ri(k)/10.**nn).lt.10.0) goto 142
1468 idx_i = INT(ri(k)/10.**n) + 10*(n-nii2) - (n-nii2)
1469 idx_i = MAX(1, MIN(idx_i, ntb_i))
1474 if (ni(k).gt. Nt_i(1)) then
1475 nii = NINT(ALOG10(ni(k)))
1476 do nn = nii-1, nii+1
1478 if ( (ni(k)/10.**nn).ge.1.0 .and. &
1479 (ni(k)/10.**nn).lt.10.0) goto 143
1482 idx_i1 = INT(ni(k)/10.**n) + 10*(n-nii3) - (n-nii3)
1483 idx_i1 = MAX(1, MIN(idx_i1, ntb_i1))
1488 !..Rain lookup table indexes.
1489 if (rr(k).gt. r_r(1)) then
1490 nir = NINT(ALOG10(rr(k)))
1491 do nn = nir-1, nir+1
1493 if ( (rr(k)/10.**nn).ge.1.0 .and. &
1494 (rr(k)/10.**nn).lt.10.0) goto 144
1497 idx_r = INT(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2)
1498 idx_r = MAX(1, MIN(idx_r, ntb_r))
1501 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
1502 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
1503 nir = NINT(DLOG10(N0_exp))
1504 do nn = nir-1, nir+1
1506 if ( (N0_exp/10.**nn).ge.1.0 .and. &
1507 (N0_exp/10.**nn).lt.10.0) goto 145
1510 idx_r1 = INT(N0_exp/10.**n) + 10*(n-nir3) - (n-nir3)
1511 idx_r1 = MAX(1, MIN(idx_r1, ntb_r1))
1517 !..Snow lookup table index.
1518 if (rs(k).gt. r_s(1)) then
1519 nis = NINT(ALOG10(rs(k)))
1520 do nn = nis-1, nis+1
1522 if ( (rs(k)/10.**nn).ge.1.0 .and. &
1523 (rs(k)/10.**nn).lt.10.0) goto 146
1526 idx_s = INT(rs(k)/10.**n) + 10*(n-nis2) - (n-nis2)
1527 idx_s = MAX(1, MIN(idx_s, ntb_s))
1532 !..Graupel lookup table index.
1533 if (rg(k).gt. r_g(1)) then
1534 nig = NINT(ALOG10(rg(k)))
1535 do nn = nig-1, nig+1
1537 if ( (rg(k)/10.**nn).ge.1.0 .and. &
1538 (rg(k)/10.**nn).lt.10.0) goto 147
1541 idx_g = INT(rg(k)/10.**n) + 10*(n-nig2) - (n-nig2)
1542 idx_g = MAX(1, MIN(idx_g, ntb_g))
1547 !..Deposition/sublimation prefactor (from Srivastava & Coen 1992).
1549 rvs = rho(k)*qvsi(k)
1550 rvs_p = rvs*otemp*(lsub*otemp*oRv - 1.)
1551 rvs_pp = rvs * ( otemp*(lsub*otemp*oRv - 1.) &
1552 *otemp*(lsub*otemp*oRv - 1.) &
1553 + (-2.*lsub*otemp*otemp*otemp*oRv) &
1555 gamsc = lsub*diffu(k)/tcond(k) * rvs_p
1556 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
1557 * rvs_pp/rvs_p * rvs/rvs_p
1558 alphsc = MAX(1.E-9, alphsc)
1560 if (abs(xsat).lt. 1.E-9) xsat=0.
1561 t1_subl = 4.*PI*( 1.0 - alphsc*xsat &
1562 + 2.*alphsc*alphsc*xsat*xsat &
1563 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
1566 !..Snow collecting cloud water. In CE, assume Dc<<Ds and vtc=~0.
1567 if (L_qc(k) .and. mvd_c(k).gt. D0c) then
1569 if (L_qs(k)) xDs = smoc(k) / smob(k)
1570 if (xDs .gt. D0s) then
1571 idx = 1 + INT(nbs*DLOG(xDs/Ds(1))/DLOG(Ds(nbs)/Ds(1)))
1573 Ef_sw = t_Efsw(idx, INT(mvd_c(k)*1.E6))
1574 prs_scw(k) = rhof(k)*t1_qs_qc*Ef_sw*rc(k)*smoe(k)
1577 !..Graupel collecting cloud water. In CE, assume Dc<<Dg and vtc=~0.
1578 if (rg(k).ge. r_g(1) .and. mvd_c(k).gt. D0c) then
1579 xDg = (bm_g + mu_g + 1.) * ilamg(k)
1580 vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
1581 stoke_g = mvd_c(k)*mvd_c(k)*vtg*rho_w/(9.*visco(k)*xDg)
1582 if (xDg.gt. D0g) then
1583 if (stoke_g.ge.0.4 .and. stoke_g.le.10.) then
1584 Ef_gw = 0.55*ALOG10(2.51*stoke_g)
1585 elseif (stoke_g.lt.0.4) then
1587 elseif (stoke_g.gt.10) then
1590 prg_gcw(k) = rhof(k)*t1_qg_qc*Ef_gw*rc(k)*N0_g(k) &
1596 !..Rain collecting snow. Cannot assume Wisner (1972) approximation
1597 !.. or Mizuno (1990) approach so we solve the CE explicitly and store
1598 !.. results in lookup table.
1599 if (rr(k).ge. r_r(1)) then
1600 if (rs(k).ge. r_s(1)) then
1601 if (temp(k).lt.T_0) then
1602 prr_rcs(k) = -(tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
1603 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
1604 + tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
1605 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r))
1606 prs_rcs(k) = tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
1607 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
1608 - tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
1609 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
1610 prg_rcs(k) = tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
1611 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
1612 + tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
1613 + tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
1614 prr_rcs(k) = MAX(DBLE(-rr(k)*odts), prr_rcs(k))
1615 prs_rcs(k) = MAX(DBLE(-rs(k)*odts), prs_rcs(k))
1616 prg_rcs(k) = MIN(DBLE((rr(k)+rs(k))*odts), prg_rcs(k))
1618 prs_rcs(k) = -(tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
1619 + tcs_racs2(idx_s,idx_t,idx_r1,idx_r))
1620 prs_rcs(k) = MAX(DBLE(-rs(k)*odts), prs_rcs(k))
1621 prr_rcs(k) = -prs_rcs(k)
1625 !..Rain collecting graupel. Cannot assume Wisner (1972) approximation
1626 !.. or Mizuno (1990) approach so we solve the CE explicitly and store
1627 !.. results in lookup table.
1628 if (rg(k).ge. r_g(1)) then
1629 if (temp(k).lt.T_0) then
1630 prg_rcg(k) = tmr_racg(idx_g,idx_r1,idx_r) &
1631 + tcr_gacr(idx_g,idx_r1,idx_r)
1632 prg_rcg(k) = MIN(DBLE(rr(k)*odts), prg_rcg(k))
1633 prr_rcg(k) = -prg_rcg(k)
1635 prr_rcg(k) = tcg_racg(idx_g,idx_r1,idx_r) &
1636 + 0.5*tmg_gacr(idx_g,idx_r1,idx_r)
1637 prr_rcg(k) = MIN(DBLE(rg(k)*odts), prr_rcg(k))
1638 prg_rcg(k) = -prr_rcg(k)
1643 !+---+-----------------------------------------------------------------+
1644 !..Next IF block handles only those processes below 0C.
1645 !+---+-----------------------------------------------------------------+
1647 if (temp(k).lt.T_0) then
1650 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
1652 !..Freezing of water drops into graupel/cloud ice (Bigg 1953).
1653 if (rr(k).gt. r_r(1)) then
1654 prg_rfz(k) = tpg_qrfz(idx_r,idx_r1,idx_tc)*odts
1655 pri_rfz(k) = tpi_qrfz(idx_r,idx_r1,idx_tc)*odts
1656 pni_rfz(k) = tni_qrfz(idx_r,idx_r1,idx_tc)*odts
1657 elseif (rr(k).gt. R1 .and. temp(k).lt.HGFR) then
1658 pri_rfz(k) = rr(k)*odts
1659 pni_rfz(k) = N0_r(k)*crg(2) * ilamr(k)**cre(2)
1661 if (rc(k).gt. r_c(1)) then
1662 pri_wfz(k) = tpi_qcfz(idx_c,idx_tc)*odts
1663 pri_wfz(k) = MIN(DBLE(rc(k)*odts), pri_wfz(k))
1664 pni_wfz(k) = tni_qcfz(idx_c,idx_tc)*odts
1665 pni_wfz(k) = MIN(DBLE(Nt_c*odts), pri_wfz(k)/(2.*xm0i), &
1669 !..Nucleate ice from deposition & condensation freezing (Cooper 1986)
1670 !.. but only if water sat and T<-10C or 20%+ ice supersaturated.
1671 if ( (ssati(k).ge. 0.20) .or. (ssatw(k).gt. eps &
1672 .and. temp(k).lt.263.15) ) then
1673 xnc = MIN(250.E3, TNO*EXP(ATO*(T_0-temp(k))))
1674 xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dtsave
1675 pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts
1676 pri_inu(k) = MIN(DBLE(rate_max), xm0i*pni_inu(k))
1677 pni_inu(k) = pri_inu(k)/xm0i
1680 !..Deposition/sublimation of cloud ice (Srivastava & Coen 1992).
1682 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
1684 xDi = MAX(DBLE(D0i), (bm_i + mu_i + 1.) * ilami)
1685 xmi = am_i*xDi**bm_i
1687 pri_ide(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
1688 *oig1*cig(5)*ni(k)*ilami
1690 if (pri_ide(k) .lt. 0.0) then
1691 pri_ide(k) = MAX(DBLE(-ri(k)*odts), pri_ide(k), DBLE(rate_max))
1692 pni_ide(k) = pri_ide(k)*oxmi
1693 pni_ide(k) = MAX(DBLE(-ni(k)*odts), pni_ide(k))
1695 pri_ide(k) = MIN(pri_ide(k), DBLE(rate_max))
1696 prs_ide(k) = (1.0D0-tpi_ide(idx_i,idx_i1))*pri_ide(k)
1697 pri_ide(k) = tpi_ide(idx_i,idx_i1)*pri_ide(k)
1700 !..Some cloud ice needs to move into the snow category. Use lookup
1701 !.. table that resulted from explicit bin representation of distrib.
1702 if ( (idx_i.eq. ntb_i) .or. (xDi.gt. 5.0*D0s) ) then
1703 prs_iau(k) = ri(k)*.99*odts
1704 pni_iau(k) = ni(k)*.95*odts
1705 elseif (xDi.lt. 0.1*D0s) then
1709 prs_iau(k) = tps_iaus(idx_i,idx_i1)*odts
1710 prs_iau(k) = MIN(DBLE(ri(k)*.99*odts), prs_iau(k))
1711 pni_iau(k) = tni_iaus(idx_i,idx_i1)*odts
1712 pni_iau(k) = MIN(DBLE(ni(k)*.95*odts), pni_iau(k))
1716 !..Deposition/sublimation of snow/graupel follows Srivastava & Coen
1719 C_snow = C_sqrd + (tempc+15.)*(C_cube-C_sqrd)/(-30.+15.)
1720 C_snow = MAX(C_sqrd, MIN(C_snow, C_cube))
1721 prs_sde(k) = C_snow*t1_subl*diffu(k)*ssati(k)*rvs &
1722 * (t1_qs_sd*smo1(k) &
1723 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
1724 if (prs_sde(k).lt. 0.) then
1725 prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k), DBLE(rate_max))
1727 prs_sde(k) = MIN(prs_sde(k), DBLE(rate_max))
1731 if (L_qg(k) .and. ssati(k).lt. -eps) then
1732 prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
1733 * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
1734 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
1735 if (prg_gde(k).lt. 0.) then
1736 prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k), DBLE(rate_max))
1738 prg_gde(k) = MIN(prg_gde(k), DBLE(rate_max))
1742 !..Snow collecting cloud ice. In CE, assume Di<<Ds and vti=~0.
1744 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
1746 xDi = MAX(DBLE(D0i), (bm_i + mu_i + 1.) * ilami)
1747 xmi = am_i*xDi**bm_i
1749 if (rs(k).ge. r_s(1)) then
1750 prs_sci(k) = t1_qs_qi*rhof(k)*Ef_si*ri(k)*smoe(k)
1751 pni_sci(k) = prs_sci(k) * oxmi
1754 !..Rain collecting cloud ice. In CE, assume Di<<Dr and vti=~0.
1755 if (rr(k).ge. r_r(1) .and. mvd_r(k).gt. 4.*xDi) then
1757 pri_rci(k) = rhof(k)*t1_qr_qi*Ef_ri*ri(k)*N0_r(k) &
1758 *((lamr+fv_r)**(-cre(9)))
1759 pni_rci(k) = pri_rci(k) * oxmi
1760 prr_rci(k) = rhof(k)*t2_qr_qi*Ef_ri*ni(k)*N0_r(k) &
1761 *((lamr+fv_r)**(-cre(8)))
1762 prr_rci(k) = MIN(DBLE(rr(k)*odts), prr_rci(k))
1763 prg_rci(k) = pri_rci(k) + prr_rci(k)
1767 !..Ice multiplication from rime-splinters (Hallet & Mossop 1974).
1768 if (prg_gcw(k).gt. eps .and. tempc.gt.-8.0) then
1770 if (tempc.ge.-5.0 .and. tempc.lt.-3.0) then
1771 tf = 0.5*(-3.0 - tempc)
1772 elseif (tempc.gt.-8.0 .and. tempc.lt.-5.0) then
1773 tf = 0.33333333*(8.0 + tempc)
1775 pni_ihm(k) = 3.5E8*tf*prg_gcw(k)
1776 pri_ihm(k) = xm0i*pni_ihm(k)
1777 prs_ihm(k) = prs_scw(k)/(prs_scw(k)+prg_gcw(k)) &
1779 prg_ihm(k) = prg_gcw(k)/(prs_scw(k)+prg_gcw(k)) &
1783 !..A portion of rimed snow converts to graupel but some remains snow.
1784 !.. Interp from 5 to 75% as riming factor increases from 5.0 to 30.0
1785 !.. 0.028 came from (.75-.05)/(30.-5.). This remains ad-hoc and should
1787 if (prs_scw(k).gt.5.0*prs_sde(k) .and. &
1788 prs_sde(k).gt.eps) then
1789 r_frac = MIN(30.0D0, prs_scw(k)/prs_sde(k))
1790 g_frac = MIN(0.75, 0.05 + (r_frac-5.)*.028)
1791 vts_boost(k) = MIN(1.5, 1.1 + (r_frac-5.)*.016)
1792 prg_scw(k) = g_frac*prs_scw(k)
1793 prs_scw(k) = (1. - g_frac)*prs_scw(k)
1798 !..Melt snow and graupel and enhance from collisions with liquid.
1799 !.. We also need to sublimate snow and graupel if subsaturated.
1801 prr_sml(k) = tempc*tcond(k)*(t1_qs_me*smo1(k) &
1802 + t2_qs_me*rhof2(k)*vsc2(k)*smof(k))
1803 prr_sml(k) = prr_sml(k) + 4218.*olfus*tempc &
1804 * (prr_rcs(k)+prs_scw(k))
1805 prr_sml(k) = MIN(DBLE(rs(k)*odts), prr_sml(k))
1807 if (ssati(k).lt. 0.) then
1808 prs_sde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
1809 * (t1_qs_sd*smo1(k) &
1810 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
1811 prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k))
1816 prr_gml(k) = tempc*N0_g(k)*tcond(k) &
1817 *(t1_qg_me*ilamg(k)**cge(10) &
1818 + t2_qg_me*rhof2(k)*vsc2(k)*ilamg(k)**cge(11))
1819 prr_gml(k) = prr_gml(k) + 4218.*olfus*tempc &
1820 * (prr_rcg(k)+prg_gcw(k))
1821 prr_gml(k) = MIN(DBLE(rg(k)*odts), prr_gml(k))
1823 if (ssati(k).lt. 0.) then
1824 prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
1825 * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
1826 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
1827 prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k))
1836 !+---+-----------------------------------------------------------------+
1837 !..Ensure we do not deplete more hydrometeor species than exists.
1838 !+---+-----------------------------------------------------------------+
1841 !..If ice supersaturated, ensure sum of depos growth terms does not
1842 !.. deplete more vapor than possibly exists. If subsaturated, limit
1843 !.. sum of sublimation terms such that vapor does not reproduce ice
1845 sump = pri_inu(k) + pri_ide(k) + prs_ide(k) &
1846 + prs_sde(k) + prg_gde(k)
1847 rate_max = (qv(k)-qvsi(k))*odts*0.999
1848 if ( (sump.gt. eps .and. sump.gt. rate_max) .or. &
1849 (sump.lt. -eps .and. sump.lt. rate_max) ) then
1850 ratio = rate_max/sump
1851 pri_inu(k) = pri_inu(k) * ratio
1852 pri_ide(k) = pri_ide(k) * ratio
1853 pni_ide(k) = pni_ide(k) * ratio
1854 prs_ide(k) = prs_ide(k) * ratio
1855 prs_sde(k) = prs_sde(k) * ratio
1856 prg_gde(k) = prg_gde(k) * ratio
1859 !..Cloud water conservation.
1860 sump = -prr_wau(k) - pri_wfz(k) - prr_rcw(k) &
1861 - prs_scw(k) - prg_scw(k) - prg_gcw(k)
1862 rate_max = -rc(k)*odts
1863 if (sump.lt. rate_max .and. L_qc(k)) then
1864 ratio = rate_max/sump
1865 prr_wau(k) = prr_wau(k) * ratio
1866 pri_wfz(k) = pri_wfz(k) * ratio
1867 prr_rcw(k) = prr_rcw(k) * ratio
1868 prs_scw(k) = prs_scw(k) * ratio
1869 prg_scw(k) = prg_scw(k) * ratio
1870 prg_gcw(k) = prg_gcw(k) * ratio
1873 !..Cloud ice conservation.
1874 sump = pri_ide(k) - prs_iau(k) - prs_sci(k) &
1876 rate_max = -ri(k)*odts
1877 if (sump.lt. rate_max .and. L_qi(k)) then
1878 ratio = rate_max/sump
1879 pri_ide(k) = pri_ide(k) * ratio
1880 prs_iau(k) = prs_iau(k) * ratio
1881 prs_sci(k) = prs_sci(k) * ratio
1882 pri_rci(k) = pri_rci(k) * ratio
1885 !..Rain conservation.
1886 sump = -prg_rfz(k) - pri_rfz(k) - prr_rci(k) &
1887 + prr_rcs(k) + prr_rcg(k)
1888 rate_max = -rr(k)*odts
1889 if (sump.lt. rate_max .and. L_qr(k)) then
1890 ratio = rate_max/sump
1891 prg_rfz(k) = prg_rfz(k) * ratio
1892 pri_rfz(k) = pri_rfz(k) * ratio
1893 prr_rci(k) = prr_rci(k) * ratio
1894 prr_rcs(k) = prr_rcs(k) * ratio
1895 prr_rcg(k) = prr_rcg(k) * ratio
1898 !..Snow conservation.
1899 sump = prs_sde(k) - prs_ihm(k) - prr_sml(k) &
1901 rate_max = -rs(k)*odts
1902 if (sump.lt. rate_max .and. L_qs(k)) then
1903 ratio = rate_max/sump
1904 prs_sde(k) = prs_sde(k) * ratio
1905 prs_ihm(k) = prs_ihm(k) * ratio
1906 prr_sml(k) = prr_sml(k) * ratio
1907 prs_rcs(k) = prs_rcs(k) * ratio
1910 !..Graupel conservation.
1911 sump = prg_gde(k) - prg_ihm(k) - prr_gml(k) &
1913 rate_max = -rg(k)*odts
1914 if (sump.lt. rate_max .and. L_qg(k)) then
1915 ratio = rate_max/sump
1916 prg_gde(k) = prg_gde(k) * ratio
1917 prg_ihm(k) = prg_ihm(k) * ratio
1918 prr_gml(k) = prr_gml(k) * ratio
1919 prg_rcg(k) = prg_rcg(k) * ratio
1924 !+---+-----------------------------------------------------------------+
1925 !..Calculate tendencies of all species but constrain the number of ice
1926 !.. to reasonable values.
1927 !+---+-----------------------------------------------------------------+
1930 lfus2 = lsub - lvap(k)
1932 !..Water vapor tendency
1933 qvten(k) = qvten(k) + (-pri_inu(k) - pri_ide(k) &
1934 - prs_ide(k) - prs_sde(k) - prg_gde(k)) &
1937 !..Cloud water tendency
1938 qcten(k) = qcten(k) + (-prr_wau(k) - pri_wfz(k) &
1939 - prr_rcw(k) - prs_scw(k) - prg_scw(k) &
1943 !..Cloud ice mixing ratio tendency
1944 qiten(k) = qiten(k) + (pri_inu(k) + pri_ihm(k) &
1945 + pri_wfz(k) + pri_rfz(k) + pri_ide(k) &
1946 - prs_iau(k) - prs_sci(k) - pri_rci(k)) &
1949 !..Cloud ice number tendency.
1950 niten(k) = niten(k) + (pni_inu(k) + pni_ihm(k) &
1951 + pni_wfz(k) + pni_rfz(k) + pni_ide(k) &
1952 - pni_iau(k) - pni_sci(k) - pni_rci(k)) &
1955 !..Cloud ice mass/number balance; keep mass-wt mean size between
1956 !.. 30 and 300 microns. Also no more than 250 xtals per liter.
1957 xri=MAX(R1,(qi1d(k) + qiten(k)*dtsave)*rho(k))
1958 xni=MAX(1.,(ni1d(k) + niten(k)*dtsave)*rho(k))
1959 if (xri.gt. R1) then
1960 lami = (am_i*cig(2)*oig1*xni/xri)**obmi
1962 xDi = (bm_i + mu_i + 1.) * ilami
1963 if (xDi.lt. 30.E-6) then
1964 lami = cie(2)/30.E-6
1965 xni = MIN(250.D3, cig(1)*oig2*xri/am_i*lami**bm_i)
1966 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
1967 elseif (xDi.gt. 300.E-6) then
1968 lami = cie(2)/300.E-6
1969 xni = cig(1)*oig2*xri/am_i*lami**bm_i
1970 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
1973 niten(k) = -ni1d(k)*odts
1975 xni=MAX(0.,(ni1d(k) + niten(k)*dtsave)*rho(k))
1976 if (xni.gt.250.E3) &
1977 niten(k) = (250.E3-ni1d(k)*rho(k))*odts*orho
1980 qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) &
1981 + prr_sml(k) + prr_gml(k) + prr_rcs(k) &
1982 + prr_rcg(k) - prg_rfz(k) &
1983 - pri_rfz(k) - prr_rci(k)) &
1987 qsten(k) = qsten(k) + (prs_iau(k) + prs_sde(k) &
1988 + prs_sci(k) + prs_scw(k) + prs_rcs(k) &
1989 + prs_ide(k) - prs_ihm(k) - prr_sml(k)) &
1993 qgten(k) = qgten(k) + (prg_scw(k) + prg_rfz(k) &
1994 + prg_gde(k) + prg_rcg(k) + prg_gcw(k) &
1995 + prg_rci(k) + prg_rcs(k) - prg_ihm(k) &
1999 !..Temperature tendency
2000 if (temp(k).lt.T_0) then
2002 + ( lsub*ocp(k)*(pri_inu(k) + pri_ide(k) &
2003 + prs_ide(k) + prs_sde(k) &
2005 + lfus2*ocp(k)*(pri_wfz(k) + pri_rfz(k) &
2006 + prg_rfz(k) + prs_scw(k) &
2007 + prg_scw(k) + prg_gcw(k) &
2008 + prg_rcs(k) + prs_rcs(k) &
2009 + prr_rci(k) + prg_rcg(k)) &
2013 + ( lfus*ocp(k)*(-prr_sml(k) - prr_gml(k) &
2014 - prr_rcg(k) - prr_rcs(k)) &
2015 + lsub*ocp(k)*(prs_sde(k) + prg_gde(k)) &
2021 !+---+-----------------------------------------------------------------+
2022 !..Update variables for TAU+1 before condensation & sedimention.
2023 !+---+-----------------------------------------------------------------+
2025 temp(k) = t1d(k) + DT*tten(k)
2027 tempc = temp(k) - 273.15
2028 qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
2029 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
2030 rhof(k) = SQRT(RHO_NOT/rho(k))
2031 rhof2(k) = SQRT(rhof(k))
2032 qvs(k) = rslf(pres(k), temp(k))
2033 ssatw(k) = qv(k)/qvs(k) - 1.
2034 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
2035 diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
2036 if (tempc .ge. 0.0) then
2037 visco(k) = (1.718+0.0049*tempc)*1.0E-5
2039 visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
2041 vsc2(k) = SQRT(rho(k)/visco(k))
2042 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
2043 tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
2044 ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
2045 lvt2(k)=lvap(k)*lvap(k)*ocp(k)*oRv*otemp*otemp
2047 if ((qc1d(k) + qcten(k)*DT) .gt. R1) then
2048 rc(k) = (qc1d(k) + qcten(k)*DT)*rho(k)
2055 if ((qi1d(k) + qiten(k)*DT) .gt. R1) then
2056 ri(k) = (qi1d(k) + qiten(k)*DT)*rho(k)
2057 ni(k) = MAX(1., (ni1d(k) + niten(k)*DT)*rho(k))
2065 if ((qr1d(k) + qrten(k)*DT) .gt. R1) then
2066 rr(k) = (qr1d(k) + qrten(k)*DT)*rho(k)
2073 if ((qs1d(k) + qsten(k)*DT) .gt. R1) then
2074 rs(k) = (qs1d(k) + qsten(k)*DT)*rho(k)
2081 if ((qg1d(k) + qgten(k)*DT) .gt. R1) then
2082 rg(k) = (qg1d(k) + qgten(k)*DT)*rho(k)
2090 !+---+-----------------------------------------------------------------+
2091 !..With tendency-updated mixing ratios, recalculate snow moments and
2092 !.. intercepts/slopes of graupel and rain.
2093 !+---+-----------------------------------------------------------------+
2094 if (.not. iiwarm) then
2096 if (.not. L_qs(k)) CYCLE
2097 tc0 = MIN(-0.1, temp(k)-273.15)
2098 smob(k) = rs(k)*oams
2100 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
2101 !.. then we must compute actual 2nd moment and use as reference.
2102 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
2105 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
2106 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
2107 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
2108 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
2109 + sa(10)*bm_s*bm_s*bm_s
2111 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
2112 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
2113 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
2114 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
2115 + sb(10)*bm_s*bm_s*bm_s
2116 smo2(k) = (smob(k)/a_)**(1./b_)
2119 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
2120 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
2121 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
2122 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
2123 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
2124 + sa(10)*cse(1)*cse(1)*cse(1)
2126 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
2127 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
2128 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
2129 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
2130 smoc(k) = a_ * smo2(k)**b_
2132 !..Calculate bm_s+bv_s (th) moment. Useful for sedimentation.
2133 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(14) &
2134 + sa(4)*tc0*cse(14) + sa(5)*tc0*tc0 &
2135 + sa(6)*cse(14)*cse(14) + sa(7)*tc0*tc0*cse(14) &
2136 + sa(8)*tc0*cse(14)*cse(14) + sa(9)*tc0*tc0*tc0 &
2137 + sa(10)*cse(14)*cse(14)*cse(14)
2139 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(14) + sb(4)*tc0*cse(14) &
2140 + sb(5)*tc0*tc0 + sb(6)*cse(14)*cse(14) &
2141 + sb(7)*tc0*tc0*cse(14) + sb(8)*tc0*cse(14)*cse(14) &
2142 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(14)*cse(14)*cse(14)
2143 smod(k) = a_ * smo2(k)**b_
2146 !+---+-----------------------------------------------------------------+
2147 !..Calculate y-intercept, slope values for graupel.
2148 !+---+-----------------------------------------------------------------+
2151 if (.not. L_qg(k)) CYCLE
2152 N0_exp = 200.0*rho(k)/rg(k)
2153 N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max)))
2154 N0_min = MIN(N0_exp, N0_min)
2156 lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
2157 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
2159 N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2)
2164 !+---+-----------------------------------------------------------------+
2165 !..Calculate y-intercept, slope values for rain.
2166 !+---+-----------------------------------------------------------------+
2169 ! if (.not. L_qr(k)) CYCLE
2170 N0_exp = ronv_c1*tanh(ronv_c0*(ronv_r0-rr(k))) + ronv_c2
2171 N0_min = MIN(N0_exp, N0_min)
2173 lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
2174 lamr = lam_exp * (crg(3)*org2*org1)**obmr
2175 mvd_r(k) = (3.0+mu_r+0.672) / lamr
2176 if (mvd_r(k) .gt. 2.5e-3) then
2178 lamr = (3.0+mu_r+0.672) / 2.5e-3
2179 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
2180 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
2182 N0_r(k) = N0_exp/(crg(2)*lam_exp) * lamr**cre(2)
2186 if (.not. iiwarm) then
2189 do k = kte-1, kts, -1
2190 if ( (temp(k).gt. T_0) .and. (rr(k).gt. 0.001e-3) &
2191 .and. ((rs(k+1)+rg(k+1)).gt. 0.01e-3) ) then
2200 !.. Locate bottom of melting layer (if any).
2202 do k = k_0-1, kts, -1
2203 if ( (rs(k)+rg(k)).lt. 0.01e-3) goto 138
2208 !.. Compute melted snow/graupel equiv water diameter one K-level above
2209 !.. melting. Set starting rain mvd to either 50 microns or max from
2210 !.. higher up in column.
2212 xDs = smoc(k_0) / smob(k_0)
2213 Ds_m = (am_s*xDs**bm_s / am_r)**obmr
2218 xDg = (bm_g + mu_g + 1.) * ilamg(k_0)
2219 Dg_m = (am_g*xDg**bm_g / am_r)**obmr
2224 r_mvd2 = MIN(MAX(Ds_m, Dg_m, r_mvd1+1.e-6, mvd_r(kbot)), &
2227 !.. Within melting layer, apply linear increase of rain mvd from r_mvd1
2228 !.. to equiv melted snow/graupel value (r_mvd2). So, by the bottom of
2229 !.. the melting layer, the rain will have an mvd that matches that from
2230 !.. melted snow and/or graupel.
2231 if (kbot.gt. 2) then
2232 do k = k_0-1, kbot, -1
2233 if (.not. L_qr(k)) CYCLE
2234 xkrat = REAL(k_0-k)/REAL(k_0-kbot)
2235 mvd_r(k) = MAX(mvd_r(k), xkrat*(r_mvd2-r_mvd1)+r_mvd1)
2236 lamr = (4.0+mu_r) / mvd_r(k)
2237 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
2238 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
2239 N0_exp = MAX(DBLE(ronv_min), MIN(N0_exp, DBLE(ronv_max)))
2240 lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
2241 lamr = lam_exp * (crg(3)*org2*org1)**obmr
2242 mvd_r(k) = (3.0+mu_r+0.672) / lamr
2243 N0_r(k) = rr(k)*lamr**cre(3) / (am_r*crg(3))
2247 !.. Below melting layer, hold N0_r constant unless changes to mixing
2248 !.. ratio increase mvd beyond 2.5 mm threshold, then adjust slope and
2249 !.. intercept to cap mvd at 2.5 mm. In future, we could lower N0_r to
2250 !.. account for self-collection or other sinks.
2251 do k = kbot-1, kts, -1
2252 if (.not. L_qr(k)) CYCLE
2253 N0_r(k) = MIN(N0_r(k), N0_r(kbot))
2254 lamr = (N0_r(k)*am_r*crg(3)/rr(k))**(1./cre(3))
2255 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
2256 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
2257 N0_exp = MAX(DBLE(ronv_min), MIN(N0_exp, DBLE(ronv_max)))
2258 lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
2259 lamr = lam_exp * (crg(3)*org2*org1)**obmr
2260 mvd_r(k) = (3.0+mu_r+0.672) / lamr
2261 if (mvd_r(k) .gt. 2.5e-3) then
2263 lamr = (3.0+mu_r+0.672) / mvd_r(k)
2265 N0_r(k) = rr(k)*lamr**cre(3) / (am_r*crg(3))
2274 !+---+-----------------------------------------------------------------+
2275 !..Cloud water condensation and evaporation. Newly formulated using
2276 !.. Newton-Raphson iterations (3 should suffice) as provided by B. Hall.
2277 !+---+-----------------------------------------------------------------+
2279 if ( (ssatw(k).gt. eps) .or. (ssatw(k).lt. -eps .and. &
2281 clap = (qv(k)-qvs(k))/(1. + lvt2(k)*qvs(k))
2283 fcd = qvs(k)* EXP(lvt2(k)*clap) - qv(k) + clap
2284 dfcd = qvs(k)*lvt2(k)* EXP(lvt2(k)*clap) + 1.
2285 clap = clap - fcd/dfcd
2288 if (xrc.gt. 0.0) then
2289 prw_vcd(k) = clap*odt
2291 prw_vcd(k) = -rc(k)/rho(k)*odt
2294 qcten(k) = qcten(k) + prw_vcd(k)
2295 qvten(k) = qvten(k) - prw_vcd(k)
2296 tten(k) = tten(k) + lvap(k)*ocp(k)*prw_vcd(k)*(1-IFDRY)
2297 rc(k) = MAX(R1, (qc1d(k) + DT*qcten(k))*rho(k))
2298 qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
2299 temp(k) = t1d(k) + DT*tten(k)
2300 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
2301 qvs(k) = rslf(pres(k), temp(k))
2302 ssatw(k) = qv(k)/qvs(k) - 1.
2306 !+---+-----------------------------------------------------------------+
2307 !.. If still subsaturated, allow rain to evaporate, following
2308 !.. Srivastava & Coen (1992).
2309 !+---+-----------------------------------------------------------------+
2311 if ( (ssatw(k).lt. -eps) .and. L_qr(k) &
2312 .and. (.not.(prw_vcd(k).gt. 0.)) ) then
2313 tempc = temp(k) - 273.15
2315 rhof(k) = SQRT(RHO_NOT/rho(k))
2316 rhof2(k) = SQRT(rhof(k))
2317 diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
2318 if (tempc .ge. 0.0) then
2319 visco(k) = (1.718+0.0049*tempc)*1.0E-5
2321 visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
2323 vsc2(k) = SQRT(rho(k)/visco(k))
2324 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
2325 tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
2326 ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
2329 rvs_p = rvs*otemp*(lvap(k)*otemp*oRv - 1.)
2330 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*oRv - 1.) &
2331 *otemp*(lvap(k)*otemp*oRv - 1.) &
2332 + (-2.*lvap(k)*otemp*otemp*otemp*oRv) &
2334 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
2335 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
2336 * rvs_pp/rvs_p * rvs/rvs_p
2337 alphsc = MAX(1.E-9, alphsc)
2339 if (xsat.lt. -1.E-9) xsat = -1.E-9
2340 t1_evap = 2.*PI*( 1.0 - alphsc*xsat &
2341 + 2.*alphsc*alphsc*xsat*xsat &
2342 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
2346 prv_rev(k) = t1_evap*diffu(k)*(-ssatw(k))*N0_r(k)*rvs &
2347 * (t1_qr_ev*ilamr(k)**cre(10) &
2348 + t2_qr_ev*vsc2(k)*rhof2(k)*((lamr+0.5*fv_r)**(-cre(11))))
2349 prv_rev(k) = MIN(DBLE(rr(k)/rho(k)*odts), &
2352 qrten(k) = qrten(k) - prv_rev(k)
2353 qvten(k) = qvten(k) + prv_rev(k)
2354 tten(k) = tten(k) - lvap(k)*ocp(k)*prv_rev(k)*(1-IFDRY)
2356 rr(k) = MAX(R1, (qr1d(k) + DT*qrten(k))*rho(k))
2357 qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
2358 temp(k) = t1d(k) + DT*tten(k)
2359 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
2363 !+---+-----------------------------------------------------------------+
2364 !..Find max terminal fallspeed (distribution mass-weighted mean
2365 !.. velocity) and use it to determine if we need to split the timestep
2366 !.. (var nstep>1). Either way, only bother to do sedimentation below
2367 !.. 1st level that contains any sedimenting particles (k=ksed1 on down).
2368 !+---+-----------------------------------------------------------------+
2371 do k = kte+1, kts, -1
2383 rhof(k) = SQRT(RHO_NOT/rho(k))
2385 if (rr(k).gt. R2) then
2387 vtr = rhof(k)*av_r*crg(6)*org3 * (lamr/(lamr+fv_r))**cre(3) &
2388 *((lamr+fv_r)**(-bv_r))
2392 if (.not. iiwarm) then
2393 if (ri(k).gt. R2) then
2394 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
2396 vti = rhof(k)*av_i*cig(3)*oig2 * ilami**bv_i
2398 vti = rhof(k)*av_i*cig(4)*oig1 * ilami**bv_i
2402 if (rs(k).gt. R2) then
2403 xDs = smoc(k) / smob(k)
2405 ils1 = 1./(Mrat*Lam0 + fv_s)
2406 ils2 = 1./(Mrat*Lam1 + fv_s)
2407 t1_vts = Kap0*csg(4)*ils1**cse(4)
2408 t2_vts = Kap1*Mrat**mu_s*csg(10)*ils2**cse(10)
2409 t3_vts = Kap0*csg(1)*ils1**cse(1)
2410 t4_vts = Kap1*Mrat**mu_s*csg(7)*ils2**cse(7)
2411 vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts)
2412 if (temp(k).gt. T_0) then
2413 vtsk(k) = MAX(vts*vts_boost(k), vtrk(k))
2415 vtsk(k) = vts*vts_boost(k)
2419 if (rg(k).gt. R2) then
2420 vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
2421 if (temp(k).gt. T_0) then
2422 vtgk(k) = MAX(vtg, vtrk(k))
2429 rgvm = MAX(vtik(k), vtrk(k), vtsk(k), vtgk(k))
2430 if (rgvm .gt. 1.E-3) then
2431 ksed1 = MAX(ksed1, k)
2432 delta_tp = dzq(k)/rgvm
2433 nstep = MAX(nstep, INT(DT/delta_tp + 1.))
2436 if (ksed1 .eq. kte) ksed1 = kte-1
2437 if (nstep .gt. 0) onstep = 1./REAL(nstep)
2439 !+---+-----------------------------------------------------------------+
2440 !..Sedimentation of mixing ratio is the integral of v(D)*m(D)*N(D)*dD,
2441 !.. whereas neglect m(D) term for number concentration. Therefore,
2442 !.. cloud ice has proper differential sedimentation.
2443 !+---+-----------------------------------------------------------------+
2446 sed_r(k) = vtrk(k)*rr(k)
2447 sed_i(k) = vtik(k)*ri(k)
2448 sed_n(k) = vtnk(k)*ni(k)
2449 sed_g(k) = vtgk(k)*rg(k)
2450 sed_s(k) = vtsk(k)*rs(k)
2455 qrten(k) = qrten(k) - sed_r(k)*odzq*onstep*orho
2456 qiten(k) = qiten(k) - sed_i(k)*odzq*onstep*orho
2457 niten(k) = niten(k) - sed_n(k)*odzq*onstep*orho
2458 qgten(k) = qgten(k) - sed_g(k)*odzq*onstep*orho
2459 qsten(k) = qsten(k) - sed_s(k)*odzq*onstep*orho
2460 rr(k) = MAX(R1, rr(k) - sed_r(k)*odzq*DT*onstep)
2461 ri(k) = MAX(R1, ri(k) - sed_i(k)*odzq*DT*onstep)
2462 ni(k) = MAX(1., ni(k) - sed_n(k)*odzq*DT*onstep)
2463 rg(k) = MAX(R1, rg(k) - sed_g(k)*odzq*DT*onstep)
2464 rs(k) = MAX(R1, rs(k) - sed_s(k)*odzq*DT*onstep)
2465 do k = ksed1, kts, -1
2468 qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) &
2470 qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) &
2472 niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) &
2474 qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) &
2476 qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) &
2478 rr(k) = MAX(R1, rr(k) + (sed_r(k+1)-sed_r(k)) &
2480 ri(k) = MAX(R1, ri(k) + (sed_i(k+1)-sed_i(k)) &
2482 ni(k) = MAX(1., ni(k) + (sed_n(k+1)-sed_n(k)) &
2484 rg(k) = MAX(R1, rg(k) + (sed_g(k+1)-sed_g(k)) &
2486 rs(k) = MAX(R1, rs(k) + (sed_s(k+1)-sed_s(k)) &
2490 !+---+-----------------------------------------------------------------+
2491 !..Precipitation reaching the ground.
2492 !+---+-----------------------------------------------------------------+
2493 pptrain = pptrain + sed_r(kts)*DT*onstep
2494 pptsnow = pptsnow + sed_s(kts)*DT*onstep
2495 pptgraul = pptgraul + sed_g(kts)*DT*onstep
2496 pptice = pptice + sed_i(kts)*DT*onstep
2500 !+---+-----------------------------------------------------------------+
2501 !.. Instantly melt any cloud ice into cloud water if above 0C and
2502 !.. instantly freeze any cloud water found below HGFR.
2503 !+---+-----------------------------------------------------------------+
2504 if (.not. iiwarm) then
2506 xri = MAX(0.0, qi1d(k) + qiten(k)*DT)
2507 if ( (temp(k).gt. T_0) .and. (xri.gt. 0.0) ) then
2508 qcten(k) = qcten(k) + xri*odt
2509 qiten(k) = -qi1d(k)*odt
2510 niten(k) = -ni1d(k)*odt
2511 tten(k) = tten(k) - lfus*ocp(k)*xri*odt*(1-IFDRY)
2514 xrc = MAX(0.0, qc1d(k) + qcten(k)*DT)
2515 if ( (temp(k).lt. HGFR) .and. (xrc.gt. 0.0) ) then
2516 lfus2 = lsub - lvap(k)
2517 qiten(k) = qiten(k) + xrc*odt
2518 niten(k) = niten(k) + xrc/(2.*xm0i)*odt
2520 tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-IFDRY)
2525 !+---+-----------------------------------------------------------------+
2526 !.. All tendencies computed, apply and pass back final values to parent.
2527 !+---+-----------------------------------------------------------------+
2529 t1d(k) = t1d(k) + tten(k)*DT
2530 qv1d(k) = MAX(1.E-10, qv1d(k) + qvten(k)*DT)
2531 qc1d(k) = qc1d(k) + qcten(k)*DT
2532 if (qc1d(k) .le. R1) qc1d(k) = 0.0
2533 qi1d(k) = qi1d(k) + qiten(k)*DT
2534 ni1d(k) = ni1d(k) + niten(k)*DT
2535 if (qi1d(k) .le. R1) then
2539 if (ni1d(k) .gt. 1.0) then
2540 lami = (am_i*cig(2)*oig1*ni1d(k)/qi1d(k))**obmi
2542 xDi = (bm_i + mu_i + 1.) * ilami
2543 if (xDi.lt. 30.E-6) then
2544 lami = cie(2)/30.E-6
2545 ni1d(k) = MIN(250.D3, cig(1)*oig2*qi1d(k)/am_i*lami**bm_i)
2546 elseif (xDi.gt. 300.E-6) then
2547 lami = cie(2)/300.E-6
2548 ni1d(k) = cig(1)*oig2*qi1d(k)/am_i*lami**bm_i
2551 lami = cie(2)/30.E-6
2552 ni1d(k) = MIN(250.D3, cig(1)*oig2*qi1d(k)/am_i*lami**bm_i)
2555 qr1d(k) = qr1d(k) + qrten(k)*DT
2556 if (qr1d(k) .le. R1) qr1d(k) = 0.0
2557 qs1d(k) = qs1d(k) + qsten(k)*DT
2558 if (qs1d(k) .le. R1) qs1d(k) = 0.0
2559 qg1d(k) = qg1d(k) + qgten(k)*DT
2560 if (qg1d(k) .le. R1) qg1d(k) = 0.0
2563 end subroutine mp_thompson
2564 !+---+-----------------------------------------------------------------+
2566 !+---+-----------------------------------------------------------------+
2567 !..Creation of the lookup tables and support functions found below here.
2568 !+---+-----------------------------------------------------------------+
2569 !..Rain collecting graupel (and inverse). Explicit CE integration.
2570 !+---+-----------------------------------------------------------------+
2572 subroutine qr_acr_qg
2577 INTEGER:: i, j, k, n, n2
2578 DOUBLE PRECISION, DIMENSION(nbg):: vg, N_g
2579 DOUBLE PRECISION, DIMENSION(nbr):: vr, N_r
2580 DOUBLE PRECISION:: N0_exp, N0_r, N0_g, lam_exp, lamg, lamr, N0_s
2581 DOUBLE PRECISION:: massg, massr, dvg, dvr, t1, t2, z1, z2
2586 vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
2589 vg(n) = av_g*Dg(n)**bv_g
2595 lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(k))**ore1
2596 lamr = lam_exp * (crg(3)*org2*org1)**obmr
2597 N0_r = N0r_exp(j)/(crg(2)*lam_exp) * lamr**cre(2)
2599 N_r(n2) = N0_r*Dr(n2)**mu_r *DEXP(-lamr*Dr(n2))*dtr(n2)
2603 N0_exp = 200.0d0/r_g(i)
2604 N0_exp = DMAX1(gonv_min*1.d0,DMIN1(N0_exp,gonv_max*1.d0))
2605 lam_exp = (N0_exp*am_g*cgg(1)/r_g(i))**oge1
2606 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
2607 N0_g = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2)
2609 N_g(n) = N0_g*Dg(n)**mu_g * DEXP(-lamg*Dg(n))*dtg(n)
2617 massr = am_r * Dr(n2)**bm_r
2619 massg = am_g * Dg(n)**bm_g
2621 dvg = 0.5d0*((vr(n2) - vg(n)) + DABS(vr(n2)-vg(n)))
2622 dvr = 0.5d0*((vg(n) - vr(n2)) + DABS(vg(n)-vr(n2)))
2624 t1 = t1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
2625 *dvg*massg * N_g(n)* N_r(n2)
2626 z1 = z1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
2627 *dvg*massr * N_g(n)* N_r(n2)
2629 t2 = t2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
2630 *dvr*massr * N_g(n)* N_r(n2)
2631 z2 = z2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
2632 *dvr*massg * N_g(n)* N_r(n2)
2636 tcg_racg(i,j,k) = t1
2637 tmr_racg(i,j,k) = DMIN1(z1, r_r(k)*1.0d0)
2638 tcr_gacr(i,j,k) = t2
2639 tmg_gacr(i,j,k) = z2
2644 end subroutine qr_acr_qg
2645 !+---+-----------------------------------------------------------------+
2647 !+---+-----------------------------------------------------------------+
2648 !..Rain collecting snow (and inverse). Explicit CE integration.
2649 !+---+-----------------------------------------------------------------+
2651 subroutine qr_acr_qs
2656 INTEGER:: i, j, k, m, n, n2
2657 DOUBLE PRECISION, DIMENSION(nbr):: vr, D1, N_r
2658 DOUBLE PRECISION, DIMENSION(nbs):: vs, N_s
2659 DOUBLE PRECISION:: loga_, a_, b_, second, M0, M2, M3, Mrat, oM3
2660 DOUBLE PRECISION:: N0_r, lam_exp, lamr, slam1, slam2
2661 DOUBLE PRECISION:: dvs, dvr, masss, massr
2662 DOUBLE PRECISION:: t1, t2, t3, t4, z1, z2, z3, z4
2667 vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
2668 D1(n2) = (vr(n2)/av_s)**(1./bv_s)
2671 vs(n) = 1.5*av_s*Ds(n)**bv_s * DEXP(-fv_s*Ds(n))
2676 lam_exp = (N0r_exp(k)*am_r*crg(1)/r_r(m))**ore1
2677 lamr = lam_exp * (crg(3)*org2*org1)**obmr
2678 N0_r = N0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2)
2680 N_r(n2) = N0_r*Dr(n2)**mu_r * DEXP(-lamr*Dr(n2))*dtr(n2)
2686 !..From the bm_s moment, compute plus one moment. If we are not
2687 !.. using bm_s=2, then we must transform to the pure 2nd moment
2688 !.. (variable called "second") and then to the bm_s+1 moment.
2690 M2 = r_s(i)*oams *1.0d0
2691 if (bm_s.gt.2.0-1.E-3 .and. bm_s.lt.2.0+1.E-3) then
2692 loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*bm_s &
2693 + sa(4)*Tc(j)*bm_s + sa(5)*Tc(j)*Tc(j) &
2694 + sa(6)*bm_s*bm_s + sa(7)*Tc(j)*Tc(j)*bm_s &
2695 + sa(8)*Tc(j)*bm_s*bm_s + sa(9)*Tc(j)*Tc(j)*Tc(j) &
2696 + sa(10)*bm_s*bm_s*bm_s
2698 b_ = sb(1) + sb(2)*Tc(j) + sb(3)*bm_s &
2699 + sb(4)*Tc(j)*bm_s + sb(5)*Tc(j)*Tc(j) &
2700 + sb(6)*bm_s*bm_s + sb(7)*Tc(j)*Tc(j)*bm_s &
2701 + sb(8)*Tc(j)*bm_s*bm_s + sb(9)*Tc(j)*Tc(j)*Tc(j) &
2702 + sb(10)*bm_s*bm_s*bm_s
2703 second = (M2/a_)**(1./b_)
2708 loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*cse(1) &
2709 + sa(4)*Tc(j)*cse(1) + sa(5)*Tc(j)*Tc(j) &
2710 + sa(6)*cse(1)*cse(1) + sa(7)*Tc(j)*Tc(j)*cse(1) &
2711 + sa(8)*Tc(j)*cse(1)*cse(1) + sa(9)*Tc(j)*Tc(j)*Tc(j) &
2712 + sa(10)*cse(1)*cse(1)*cse(1)
2714 b_ = sb(1)+sb(2)*Tc(j)+sb(3)*cse(1) + sb(4)*Tc(j)*cse(1) &
2715 + sb(5)*Tc(j)*Tc(j) + sb(6)*cse(1)*cse(1) &
2716 + sb(7)*Tc(j)*Tc(j)*cse(1) + sb(8)*Tc(j)*cse(1)*cse(1) &
2717 + sb(9)*Tc(j)*Tc(j)*Tc(j)+sb(10)*cse(1)*cse(1)*cse(1)
2718 M3 = a_ * second**b_
2721 Mrat = M2*(M2*oM3)*(M2*oM3)*(M2*oM3)
2723 slam1 = M2 * oM3 * Lam0
2724 slam2 = M2 * oM3 * Lam1
2727 N_s(n) = Mrat*(Kap0*DEXP(-slam1*Ds(n)) &
2728 + Kap1*M0*Ds(n)**mu_s * DEXP(-slam2*Ds(n)))*dts(n)
2740 massr = am_r * Dr(n2)**bm_r
2742 masss = am_s * Ds(n)**bm_s
2744 dvs = 0.5d0*((vr(n2) - vs(n)) + DABS(vr(n2)-vs(n)))
2745 dvr = 0.5d0*((vs(n) - vr(n2)) + DABS(vs(n)-vr(n2)))
2747 if (massr .gt. masss) then
2748 t1 = t1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2749 *dvs*masss * N_s(n)* N_r(n2)
2750 z1 = z1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2751 *dvs*massr * N_s(n)* N_r(n2)
2753 t3 = t3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2754 *dvs*masss * N_s(n)* N_r(n2)
2755 z3 = z3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2756 *dvs*massr * N_s(n)* N_r(n2)
2759 if (massr .gt. masss) then
2760 t2 = t2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2761 *dvr*massr * N_s(n)* N_r(n2)
2762 z2 = z2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2763 *dvr*masss * N_s(n)* N_r(n2)
2765 t4 = t4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2766 *dvr*massr * N_s(n)* N_r(n2)
2767 z4 = z4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2768 *dvr*masss * N_s(n)* N_r(n2)
2773 tcs_racs1(i,j,k,m) = t1
2774 tmr_racs1(i,j,k,m) = DMIN1(z1, r_r(m)*1.0d0)
2775 tcs_racs2(i,j,k,m) = t3
2776 tmr_racs2(i,j,k,m) = z3
2777 tcr_sacr1(i,j,k,m) = t2
2778 tms_sacr1(i,j,k,m) = z2
2779 tcr_sacr2(i,j,k,m) = t4
2780 tms_sacr2(i,j,k,m) = z4
2786 end subroutine qr_acr_qs
2787 !+---+-----------------------------------------------------------------+
2789 !+---+-----------------------------------------------------------------+
2790 !..This is a literal adaptation of Bigg (1954) probability of drops of
2791 !..a particular volume freezing. Given this probability, simply freeze
2792 !..the proportion of drops summing their masses.
2793 !+---+-----------------------------------------------------------------+
2795 subroutine freezeH2O
2800 INTEGER:: i, j, k, n, n2
2801 DOUBLE PRECISION, DIMENSION(nbr):: N_r, massr
2802 DOUBLE PRECISION, DIMENSION(nbc):: N_c, massc
2803 DOUBLE PRECISION:: sum1, sum2, sumn1, sumn2, &
2804 prob, vol, Texp, orho_w, &
2805 lam_exp, lamr, N0_r, lamc, N0_c, y
2812 massr(n2) = am_r*Dr(n2)**bm_r
2815 massc(n) = am_r*Dc(n)**bm_r
2818 !..Freeze water (smallest drops become cloud ice, otherwise graupel).
2820 ! print*, ' Freezing water for temp = ', -k
2821 Texp = DEXP( DFLOAT(k) ) - 1.0D0
2824 lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(i))**ore1
2825 lamr = lam_exp * (crg(3)*org2*org1)**obmr
2826 N0_r = N0r_exp(j)/(crg(2)*lam_exp) * lamr**cre(2)
2831 N_r(n2) = N0_r*Dr(n2)**mu_r*DEXP(-lamr*Dr(n2))*dtr(n2)
2832 vol = massr(n2)*orho_w
2833 prob = 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp)
2834 if (massr(n2) .lt. xm0g) then
2835 sumn1 = sumn1 + prob*N_r(n2)
2836 sum1 = sum1 + prob*N_r(n2)*massr(n2)
2838 sum2 = sum2 + prob*N_r(n2)*massr(n2)
2841 tpi_qrfz(i,j,k) = sum1
2842 tni_qrfz(i,j,k) = sumn1
2843 tpg_qrfz(i,j,k) = sum2
2847 lamc = 1.0D-6 * (Nt_c*am_r* ccg(2) * ocg1 / r_c(i))**obmr
2848 N0_c = 1.0D-18 * Nt_c*ocg1 * lamc**cce(1)
2853 vol = massc(n)*orho_w
2854 prob = 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp)
2855 N_c(n) = N0_c* y**mu_c * EXP(-lamc*y)*dtc(n)
2856 N_c(n) = 1.0D24 * N_c(n)
2857 sumn2 = sumn2 + prob*N_c(n)
2858 sum1 = sum1 + prob*N_c(n)*massc(n)
2860 tpi_qcfz(i,k) = sum1
2861 tni_qcfz(i,k) = sumn2
2865 end subroutine freezeH2O
2866 !+---+-----------------------------------------------------------------+
2868 !+---+-----------------------------------------------------------------+
2869 !..Cloud ice converting to snow since portion greater than min snow
2870 !.. size. Given cloud ice content (kg/m**3), number concentration
2871 !.. (#/m**3) and gamma shape parameter, mu_i, break the distrib into
2872 !.. bins and figure out the mass/number of ice with sizes larger than
2873 !.. D0s. Also, compute incomplete gamma function for the integration
2874 !.. of ice depositional growth from diameter=0 to D0s. Amount of
2875 !.. ice depositional growth is this portion of distrib while larger
2876 !.. diameters contribute to snow growth (as in Harrington et al. 1995).
2877 !+---+-----------------------------------------------------------------+
2879 subroutine qi_aut_qs
2885 DOUBLE PRECISION, DIMENSION(nbi):: N_i
2886 DOUBLE PRECISION:: N0_i, lami, Di_mean, t1, t2
2892 lami = (am_i*cig(2)*oig1*Nt_i(j)/r_i(i))**obmi
2893 Di_mean = (bm_i + mu_i + 1.) / lami
2894 N0_i = Nt_i(j)*oig1 * lami**cie(1)
2897 if (SNGL(Di_mean) .gt. 5.*D0s) then
2900 tpi_ide(i,j) = 0.0D0
2901 elseif (SNGL(Di_mean) .lt. D0i) then
2904 tpi_ide(i,j) = 1.0D0
2906 #if (DWORDSIZE == 8 && RWORDSIZE == 8)
2907 tpi_ide(i,j) = GAMMP(mu_i+2.0, REAL(lami,KIND=8)*D0s) * 1.0D0
2908 #elif (DWORDSIZE == 8 && RWORDSIZE == 4)
2909 tpi_ide(i,j) = GAMMP(mu_i+2.0, REAL(lami,KIND=4)*D0s) * 1.0D0
2911 This is a temporary hack assuming double precision is 8 bytes.
2914 N_i(n2) = N0_i*Di(n2)**mu_i * DEXP(-lami*Di(n2))*dti(n2)
2915 if (Di(n2).ge.D0s) then
2916 t1 = t1 + N_i(n2) * am_i*Di(n2)**bm_i
2926 end subroutine qi_aut_qs
2928 !+---+-----------------------------------------------------------------+
2929 !..Variable collision efficiency for rain collecting cloud water using
2930 !.. method of Beard and Grover, 1974 if a/A less than 0.25; otherwise
2931 !.. uses polynomials to get close match of Pruppacher & Klett Fig 14-9.
2932 !+---+-----------------------------------------------------------------+
2934 subroutine table_Efrw
2939 DOUBLE PRECISION:: vtr, stokes, reynolds, Ef_rw
2940 DOUBLE PRECISION:: p, yc0, F, G, H, z, K0, X
2947 if (Dr(i).lt.50.E-6 .or. Dc(j).lt.3.E-6) then
2949 elseif (p.gt.0.25) then
2951 if (Dr(i) .lt. 75.e-6) then
2952 Ef_rw = 0.026794*X - 0.20604
2953 elseif (Dr(i) .lt. 125.e-6) then
2954 Ef_rw = -0.00066842*X*X + 0.061542*X - 0.37089
2955 elseif (Dr(i) .lt. 175.e-6) then
2956 Ef_rw = 4.091e-06*X*X*X*X - 0.00030908*X*X*X &
2957 + 0.0066237*X*X - 0.0013687*X - 0.073022
2958 elseif (Dr(i) .lt. 250.e-6) then
2959 Ef_rw = 9.6719e-5*X*X*X - 0.0068901*X*X + 0.17305*X &
2961 elseif (Dr(i) .lt. 350.e-6) then
2962 Ef_rw = 9.0488e-5*X*X*X - 0.006585*X*X + 0.16606*X &
2965 Ef_rw = 0.00010721*X*X*X - 0.0072962*X*X + 0.1704*X &
2969 vtr = -0.1021 + 4.932E3*Dr(i) - 0.9551E6*Dr(i)*Dr(i) &
2970 + 0.07934E9*Dr(i)*Dr(i)*Dr(i) &
2971 - 0.002362E12*Dr(i)*Dr(i)*Dr(i)*Dr(i)
2972 stokes = Dc(j)*Dc(j)*vtr*rho_w/(9.*1.718E-5*Dr(i))
2973 reynolds = 9.*stokes/(p*p*rho_w)
2976 G = -0.1007D0 - 0.358D0*F + 0.0261D0*F*F
2978 z = DLOG(stokes/(K0+1.D-15))
2979 H = 0.1465D0 + 1.302D0*z - 0.607D0*z*z + 0.293D0*z*z*z
2980 yc0 = 2.0D0/PI * ATAN(H)
2981 Ef_rw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p))
2985 t_Efrw(i,j) = MAX(0.0, MIN(SNGL(Ef_rw), 0.95))
2990 end subroutine table_Efrw
2992 !+---+-----------------------------------------------------------------+
2993 !..Variable collision efficiency for snow collecting cloud water using
2994 !.. method of Wang and Ji, 2000 except equate melted snow diameter to
2995 !.. their "effective collision cross-section."
2996 !+---+-----------------------------------------------------------------+
2998 subroutine table_Efsw
3003 DOUBLE PRECISION:: Ds_m, vts, vtc, stokes, reynolds, Ef_sw
3004 DOUBLE PRECISION:: p, yc0, F, G, H, z, K0
3008 vtc = 1.19D4 * (1.0D4*Dc(j)*Dc(j)*0.25D0)
3010 vts = av_s*Ds(i)**bv_s * DEXP(-fv_s*Ds(i)) - vtc
3011 Ds_m = (am_s*Ds(i)**bm_s / am_r)**obmr
3013 if (p.gt.0.25 .or. Ds(i).lt.D0s .or. Dc(j).lt.6.E-6 &
3014 .or. vts.lt.1.E-3) then
3017 stokes = Dc(j)*Dc(j)*vts*rho_w/(9.*1.718E-5*Ds_m)
3018 reynolds = 9.*stokes/(p*p*rho_w)
3021 G = -0.1007D0 - 0.358D0*F + 0.0261D0*F*F
3023 z = DLOG(stokes/(K0+1.D-15))
3024 H = 0.1465D0 + 1.302D0*z - 0.607D0*z*z + 0.293D0*z*z*z
3025 yc0 = 2.0D0/PI * ATAN(H)
3026 Ef_sw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p))
3028 t_Efsw(i,j) = MAX(0.0, MIN(SNGL(Ef_sw), 0.95))
3034 end subroutine table_Efsw
3036 !+---+-----------------------------------------------------------------+
3037 !+---+-----------------------------------------------------------------+
3038 SUBROUTINE GCF(GAMMCF,A,X,GLN)
3039 ! --- RETURNS THE INCOMPLETE GAMMA FUNCTION Q(A,X) EVALUATED BY ITS
3040 ! --- CONTINUED FRACTION REPRESENTATION AS GAMMCF. ALSO RETURNS
3041 ! --- LN(GAMMA(A)) AS GLN. THE CONTINUED FRACTION IS EVALUATED BY
3042 ! --- A MODIFIED LENTZ METHOD.
3045 INTEGER, PARAMETER:: ITMAX=100
3046 REAL, PARAMETER:: gEPS=3.E-7
3047 REAL, PARAMETER:: FPMIN=1.E-30
3048 REAL, INTENT(IN):: A, X
3051 REAL:: AN,B,C,D,DEL,H
3061 IF(ABS(D).LT.FPMIN)D=FPMIN
3063 IF(ABS(C).LT.FPMIN)C=FPMIN
3067 IF(ABS(DEL-1.).LT.gEPS)GOTO 1
3069 PRINT *, 'A TOO LARGE, ITMAX TOO SMALL IN GCF'
3070 1 GAMMCF=EXP(-X+A*LOG(X)-GLN)*H
3072 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
3073 !+---+-----------------------------------------------------------------+
3074 SUBROUTINE GSER(GAMSER,A,X,GLN)
3075 ! --- RETURNS THE INCOMPLETE GAMMA FUNCTION P(A,X) EVALUATED BY ITS
3076 ! --- ITS SERIES REPRESENTATION AS GAMSER. ALSO RETURNS LN(GAMMA(A))
3080 INTEGER, PARAMETER:: ITMAX=100
3081 REAL, PARAMETER:: gEPS=3.E-7
3082 REAL, INTENT(IN):: A, X
3088 IF(X.LT.0.) PRINT *, 'X < 0 IN GSER'
3099 IF(ABS(DEL).LT.ABS(SUM)*gEPS)GOTO 1
3101 PRINT *,'A TOO LARGE, ITMAX TOO SMALL IN GSER'
3102 1 GAMSER=SUM*EXP(-X+A*LOG(X)-GLN)
3104 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
3105 !+---+-----------------------------------------------------------------+
3106 REAL FUNCTION GAMMLN(XX)
3107 ! --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 0.
3109 REAL, INTENT(IN):: XX
3110 DOUBLE PRECISION, PARAMETER:: STP = 2.5066282746310005D0
3111 DOUBLE PRECISION, DIMENSION(6), PARAMETER:: &
3112 COF = (/76.18009172947146D0, -86.50532032941677D0, &
3113 24.01409824083091D0, -1.231739572450155D0, &
3114 .1208650973866179D-2, -.5395239384953D-5/)
3115 DOUBLE PRECISION:: SER,TMP,X,Y
3121 TMP=(X+0.5D0)*LOG(TMP)-TMP
3122 SER=1.000000000190015D0
3127 GAMMLN=TMP+LOG(STP*SER/X)
3129 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
3130 !+---+-----------------------------------------------------------------+
3131 REAL FUNCTION GAMMP(A,X)
3132 ! --- COMPUTES THE INCOMPLETE GAMMA FUNCTION P(A,X)
3133 ! --- SEE ABRAMOWITZ AND STEGUN 6.5.1
3136 REAL, INTENT(IN):: A,X
3137 REAL:: GAMMCF,GAMSER,GLN
3139 IF((X.LT.0.) .OR. (A.LE.0.)) THEN
3140 PRINT *, 'BAD ARGUMENTS IN GAMMP'
3142 ELSEIF(X.LT.A+1.)THEN
3143 CALL GSER(GAMSER,A,X,GLN)
3146 CALL GCF(GAMMCF,A,X,GLN)
3150 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
3151 !+---+-----------------------------------------------------------------+
3152 REAL FUNCTION WGAMMA(y)
3155 REAL, INTENT(IN):: y
3157 WGAMMA = EXP(GAMMLN(y))
3160 !+---+-----------------------------------------------------------------+
3161 ! THIS FUNCTION CALCULATES THE LIQUID SATURATION VAPOR MIXING RATIO AS
3162 ! A FUNCTION OF TEMPERATURE AND PRESSURE
3164 REAL FUNCTION RSLF(P,T)
3167 REAL, INTENT(IN):: P, T
3169 REAL, PARAMETER:: C0= .611583699E03
3170 REAL, PARAMETER:: C1= .444606896E02
3171 REAL, PARAMETER:: C2= .143177157E01
3172 REAL, PARAMETER:: C3= .264224321E-1
3173 REAL, PARAMETER:: C4= .299291081E-3
3174 REAL, PARAMETER:: C5= .203154182E-5
3175 REAL, PARAMETER:: C6= .702620698E-8
3176 REAL, PARAMETER:: C7= .379534310E-11
3177 REAL, PARAMETER:: C8=-.321582393E-13
3179 X=MAX(-80.,T-273.16)
3181 ! ESL=612.2*EXP(17.67*X/(T-29.65))
3182 ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
3183 RSLF=.622*ESL/(P-ESL)
3187 !+---+-----------------------------------------------------------------+
3188 ! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A
3189 ! FUNCTION OF TEMPERATURE AND PRESSURE
3191 REAL FUNCTION RSIF(P,T)
3194 REAL, INTENT(IN):: P, T
3196 REAL, PARAMETER:: C0= .609868993E03
3197 REAL, PARAMETER:: C1= .499320233E02
3198 REAL, PARAMETER:: C2= .184672631E01
3199 REAL, PARAMETER:: C3= .402737184E-1
3200 REAL, PARAMETER:: C4= .565392987E-3
3201 REAL, PARAMETER:: C5= .521693933E-5
3202 REAL, PARAMETER:: C6= .307839583E-7
3203 REAL, PARAMETER:: C7= .105785160E-9
3204 REAL, PARAMETER:: C8= .161444444E-12
3206 X=MAX(-80.,T-273.16)
3207 ESI=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
3208 RSIF=.622*ESI/(P-ESI)
3211 !+---+-----------------------------------------------------------------+
3212 END MODULE module_mp_thompson07
3213 !+---+-----------------------------------------------------------------+
3215 ! MODIFICATIONS TO MAKE IN OTHER MODULES (pre v2.2 code only)
3217 ! Use this new code by changing the "THOMPSON" section of code found
3218 ! in "module_microphysics_driver.F" with this section. [Of course
3219 ! remove the leading comment character that you see here.]
3222 ! CALL wrf_debug ( 100 , 'microphysics_driver: calling thompson et al' )
3223 ! IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. &
3224 ! PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. &
3225 ! PRESENT( QS_CURR ) .AND. PRESENT ( QG_CURR ) .AND. &
3226 ! PRESENT ( QNI_CURR ).AND. &
3227 ! PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) ) THEN
3228 ! CALL mp_gt_driver( &
3241 ! ITIMESTEP=itimestep, &
3243 ! RAINNCV=RAINNCV, &
3245 ! ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
3246 ! ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
3247 ! ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte)
3249 ! CALL wrf_error_fatal ( 'arguments not present for calling thompson_et_al' )
3252 ! Then rename the call from "thomp_init" to "thompson_init" in the file
3253 ! "module_physics_init.F" (seen below):
3256 ! CALL thompson_init