module_mp_thompson.F
References to this file elsewhere.
1 !+---+-----------------------------------------------------------------+
2 !.. This subroutine computes the moisture tendencies of water vapor,
3 !.. cloud droplets, rain, cloud ice (pristine), snow, and graupel.
4 !.. Previously this code was based on Reisner et al (1998), but few of
5 !.. those pieces remain. A complete description is now found in
6 !.. Thompson et al. (2004, 2006).
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.
21 !..
22 !..Author: Greg Thompson, NCAR-RAL, gthompsn@ucar.edu, 303-497-2805
23 !..Last modified: 14 Feb 2007
24 !+---+-----------------------------------------------------------------+
25 !wrft:model_layer:physics
26 !+---+-----------------------------------------------------------------+
27 !
28 MODULE module_mp_thompson
29 USE module_wrf_error
30
31 IMPLICIT NONE
32
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
37
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
43
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
48 !.. scheme.
49 REAL, PARAMETER, PRIVATE:: Nt_c = 100.E6
50
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
56 REAL, PRIVATE:: mu_c
57
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
68
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 = 1.E6
75 REAL, PARAMETER, PRIVATE:: ronv_min = 2.E6
76 REAL, PARAMETER, PRIVATE:: ronv_max = 2.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
82
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 = 0.069
92 REAL, PARAMETER, PRIVATE:: bm_i = 2.0
93
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 = 2247.0
106 REAL, PARAMETER, PRIVATE:: bv_i = 1.0
107
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
111
112 !..Collection efficiencies. Rain/snow/graupel collection of cloud
113 !.. droplets use variables (Ef_rw, Ef_ri, Ef_sw, Ef_gw respectively) and
114 !.. get computed elsewhere because they are dependent on stokes
115 !.. number.
116 REAL, PARAMETER, PRIVATE:: Ef_si = 0.1
117 REAL, PARAMETER, PRIVATE:: Ef_rs = 0.99
118 REAL, PARAMETER, PRIVATE:: Ef_rg = 0.99
119
120 !..Minimum microphys values
121 !.. R1 value, 1.E-12, cannot be set lower because of numerical
122 !.. problems with Paul Field's moments and should not be set larger
123 !.. because of truncation problems in snow/ice growth.
124 REAL, PARAMETER, PRIVATE:: R1 = 1.E-12
125 REAL, PARAMETER, PRIVATE:: R2 = 1.E-8
126 REAL, PARAMETER, PRIVATE:: eps = 1.E-29
127
128 !..Constants in Cooper curve relation for cloud ice number.
129 REAL, PARAMETER, PRIVATE:: TNO = 5.0
130 REAL, PARAMETER, PRIVATE:: ATO = 0.304
131
132 !..Rho_not used in fallspeed relations (rho_not/rho)**.5 adjustment.
133 REAL, PARAMETER, PRIVATE:: rho_not = 101325.0/(287.05*298.0)
134
135 !..Schmidt number
136 REAL, PARAMETER, PRIVATE:: Sc = 0.632
137 REAL, PRIVATE:: Sc3
138
139 !..Homogeneous freezing temperature
140 REAL, PARAMETER, PRIVATE:: HGFR = 235.16
141
142 !..Water vapor and air gas constants at constant pressure
143 REAL, PARAMETER, PRIVATE:: Rv = 461.5
144 REAL, PARAMETER, PRIVATE:: oRv = 1./Rv
145 REAL, PARAMETER, PRIVATE:: R = 287.04
146 REAL, PARAMETER, PRIVATE:: Cp = 1004.0
147
148 !..Enthalpy of sublimation, vaporization, and fusion at 0C.
149 REAL, PARAMETER, PRIVATE:: lsub = 2.834E6
150 REAL, PARAMETER, PRIVATE:: lvap0 = 2.5E6
151 REAL, PARAMETER, PRIVATE:: lfus = lsub - lvap0
152 REAL, PARAMETER, PRIVATE:: olfus = 1./lfus
153
154 !..Ice initiates with this mass (kg), corresponding diameter calc.
155 !..Min diameters and mass of cloud, rain, snow, and graupel (m, kg).
156 REAL, PARAMETER, PRIVATE:: xm0i = 1.E-12
157 REAL, PARAMETER, PRIVATE:: D0c = 1.E-6
158 REAL, PARAMETER, PRIVATE:: D0r = 50.E-6
159 REAL, PARAMETER, PRIVATE:: D0s = 125.E-6
160 REAL, PARAMETER, PRIVATE:: D0g = 150.E-6
161 REAL, PRIVATE:: D0i, xm0s, xm0g
162
163 !..Lookup table dimensions
164 INTEGER, PARAMETER, PRIVATE:: ntb_c = 37
165 INTEGER, PARAMETER, PRIVATE:: ntb_i = 64
166 INTEGER, PARAMETER, PRIVATE:: ntb_r = 37
167 INTEGER, PARAMETER, PRIVATE:: ntb_s = 28
168 INTEGER, PARAMETER, PRIVATE:: ntb_g = 28
169 INTEGER, PARAMETER, PRIVATE:: ntb_r1 = 37
170 INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55
171 INTEGER, PARAMETER, PRIVATE:: ntb_t = 9
172 INTEGER, PRIVATE:: nic2, nii2, nii3, nir2, nir3, nis2, nig2
173
174 !..Lookup tables for cloud water content (kg/m**3).
175 REAL, DIMENSION(ntb_c), PARAMETER, PRIVATE:: &
176 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, &
177 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, &
178 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, &
179 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, &
180 1.e-2/)
181
182 !..Lookup tables for cloud ice content (kg/m**3).
183 REAL, DIMENSION(ntb_i), PARAMETER, PRIVATE:: &
184 r_i = (/1.e-10,2.e-10,3.e-10,4.e-10, &
185 5.e-10,6.e-10,7.e-10,8.e-10,9.e-10, &
186 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, &
187 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, &
188 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, &
189 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, &
190 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, &
191 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, &
192 1.e-3/)
193
194 !..Lookup tables for rain content (kg/m**3).
195 REAL, DIMENSION(ntb_r), PARAMETER, PRIVATE:: &
196 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, &
197 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, &
198 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, &
199 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, &
200 1.e-2/)
201
202 !..Lookup tables for graupel content (kg/m**3).
203 REAL, DIMENSION(ntb_g), PARAMETER, PRIVATE:: &
204 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, &
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, &
206 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, &
207 1.e-2/)
208
209 !..Lookup tables for snow content (kg/m**3).
210 REAL, DIMENSION(ntb_s), PARAMETER, PRIVATE:: &
211 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, &
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, &
214 1.e-2/)
215
216 !..Lookup tables for rain y-intercept parameter (/m**4).
217 REAL, DIMENSION(ntb_r1), PARAMETER, PRIVATE:: &
218 N0r_exp = (/1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, &
219 1.e7,2.e7,3.e7,4.e7,5.e7,6.e7,7.e7,8.e7,9.e7, &
220 1.e8,2.e8,3.e8,4.e8,5.e8,6.e8,7.e8,8.e8,9.e8, &
221 1.e9,2.e9,3.e9,4.e9,5.e9,6.e9,7.e9,8.e9,9.e9, &
222 1.e10/)
223
224 !..Lookup tables for ice number concentration (/m**3).
225 REAL, DIMENSION(ntb_i1), PARAMETER, PRIVATE:: &
226 Nt_i = (/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0, &
227 1.e1,2.e1,3.e1,4.e1,5.e1,6.e1,7.e1,8.e1,9.e1, &
228 1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, &
229 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, &
230 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
231 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
232 1.e6/)
233
234 !..For snow moments conversions (from Field et al. 2005)
235 REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
236 sa = (/ 5.065339, -0.062659, -3.032362, 0.029469, -0.000285, &
237 0.31255, 0.000204, 0.003199, 0.0, -0.015952/)
238 REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
239 sb = (/ 0.476221, -0.015896, 0.165977, 0.007468, -0.000141, &
240 0.060366, 0.000079, 0.000594, 0.0, -0.003577/)
241
242 !..Lookup table for A1 and A2 in Bergeron process (Koenig 1971)
243 REAL, DIMENSION(31), PARAMETER, PRIVATE:: &
244 ABER1 = (/.7939E-07,.7841E-06,.3369E-05,.4336E-05, &
245 .5285E-05,.3728E-05,.1852E-05,.2991E-06,.4248E-06, &
246 .7434E-06,.1812E-05,.4394E-05,.9145E-05,.1725E-04, &
247 .3348E-04,.1725E-04,.9175E-05,.4412E-05,.2252E-05, &
248 .9115E-06,.4876E-06,.3473E-06,.4758E-06,.6306E-06, &
249 .8573E-06,.7868E-06,.7192E-06,.6513E-06,.5956E-06, &
250 .5333E-06,.4834E-06/)
251 REAL, DIMENSION(31), PARAMETER, PRIVATE:: &
252 ABER2 = (/.4006,.4831,.5320,.5307,.5319,.5249, &
253 .4888,.3894,.4047,.4318,.4771,.5183,.5463,.5651, &
254 .5813,.5655,.5478,.5203,.4906,.4447,.4126,.3960, &
255 .4149,.4320,.4506,.4483,.4460,.4433,.4413,.4382, &
256 .4361/)
257 REAL, DIMENSION(31), PRIVATE:: CBG
258
259 !..Temperatures (5 C interval 0 to -40) used in lookup tables.
260 REAL, DIMENSION(ntb_t), PARAMETER, PRIVATE:: &
261 Tc = (/-0.01, -5., -10., -15., -20., -25., -30., -35., -40./)
262
263 !..Lookup tables for various accretion/collection terms.
264 !.. ntb_x refers to the number of elements for rain, snow, graupel,
265 !.. and temperature array indices. Variables beginning with tp/tc/tm
266 !.. represent lookup tables.
267 DOUBLE PRECISION, DIMENSION(ntb_g,ntb_r1,ntb_r):: &
268 tcg_racg, tmr_racg, tcr_gacr, tmg_gacr
269 DOUBLE PRECISION, DIMENSION(ntb_s,ntb_t,ntb_r1,ntb_r):: &
270 tcs_racs1, tmr_racs1, tcs_racs2, tmr_racs2, &
271 tcr_sacr1, tms_sacr1, tcr_sacr2, tms_sacr2
272 DOUBLE PRECISION, DIMENSION(ntb_c,45):: &
273 tpi_qcfz, tni_qcfz
274 DOUBLE PRECISION, DIMENSION(ntb_r,ntb_r1,45):: &
275 tpi_qrfz, tpg_qrfz, tni_qrfz
276 DOUBLE PRECISION, DIMENSION(ntb_i,ntb_i1):: &
277 tps_iaus, tni_iaus, tpi_ide
278
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
289 REAL, DIMENSION(12), PRIVATE:: cge, cgg
290 REAL, PRIVATE:: oge1, ogg1, ogg2, ogg3, obmg
291
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
297
298 !+---+
299 !+---+-----------------------------------------------------------------+
300 !..END DECLARATIONS
301 !+---+-----------------------------------------------------------------+
302 !+---+
303 !
304
305 CONTAINS
306
307 SUBROUTINE thompson_init
308
309 IMPLICIT NONE
310
311 INTEGER:: i, j, k, m, n
312
313 !..From Martin et al. (1994), assign gamma shape parameter mu for cloud
314 !.. drops according to general dispersion characteristics (disp=~0.25
315 !.. for Maritime and 0.45 for Continental).
316 !.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime
317 !.. to 2 for really dirty air.
318 mu_c = MIN(15., (1000.E6/Nt_c + 2.))
319
320 !..Schmidt number to one-third used numerous times.
321 Sc3 = Sc**(1./3.)
322
323 !..Compute min ice diam from mass, min snow/graupel mass from diam.
324 D0i = (xm0i/am_i)**(1./bm_i)
325 xm0s = am_s * D0s**bm_s
326 xm0g = am_g * D0g**bm_g
327
328 !..These constants various exponents and gamma() assoc with cloud,
329 !.. rain, snow, and graupel.
330 cce(1) = mu_c + 1.
331 cce(2) = bm_r + mu_c + 1.
332 cce(3) = bm_r + mu_c + 4.
333 ccg(1) = WGAMMA(cce(1))
334 ccg(2) = WGAMMA(cce(2))
335 ccg(3) = WGAMMA(cce(3))
336 ocg1 = 1./ccg(1)
337 ocg2 = 1./ccg(2)
338
339 cie(1) = mu_i + 1.
340 cie(2) = bm_i + mu_i + 1.
341 cie(3) = bm_i + mu_i + bv_i + 1.
342 cie(4) = mu_i + bv_i + 1.
343 cie(5) = mu_i + 2.
344 cie(6) = bm_i + bv_i
345 cig(1) = WGAMMA(cie(1))
346 cig(2) = WGAMMA(cie(2))
347 cig(3) = WGAMMA(cie(3))
348 cig(4) = WGAMMA(cie(4))
349 cig(5) = WGAMMA(cie(5))
350 cig(6) = WGAMMA(cie(6))
351 oig1 = 1./cig(1)
352 oig2 = 1./cig(2)
353 obmi = 1./bm_i
354
355 do n = 1, 31
356 CBG(n) = WGAMMA(mu_i + 1 + ABER2(n)*bm_i)
357 enddo
358
359 cre(1) = bm_r + 1.
360 cre(2) = mu_r + 1.
361 cre(3) = bm_r + mu_r + 1.
362 cre(4) = bm_r + mu_r + 2.
363 cre(5) = bm_r + mu_r + 3.
364 cre(6) = bm_r + mu_r + bv_r + 1.
365 cre(7) = bm_r + mu_r + bv_r + 2.
366 cre(8) = bm_r + mu_r + bv_r + 3.
367 cre(9) = mu_r + bv_r + 3.
368 cre(10) = mu_r + 2.
369 cre(11) = 0.5*(bv_r + 5. + 2.*mu_r)
370 cre(12) = bm_r + mu_r + 4.
371 do n = 1, 12
372 crg(n) = WGAMMA(cre(n))
373 enddo
374 obmr = 1./bm_r
375 ore1 = 1./cre(1)
376 org1 = 1./crg(1)
377 org2 = 1./crg(2)
378 org3 = 1./crg(3)
379
380 cse(1) = bm_s + 1.
381 cse(2) = bm_s + 2.
382 cse(3) = bm_s + 3.
383 cse(4) = bm_s + bv_s + 1.
384 cse(5) = bm_s + bv_s + 2.
385 cse(6) = bm_s + bv_s + 3.
386 cse(7) = bm_s + mu_s + 1.
387 cse(8) = bm_s + mu_s + 2.
388 cse(9) = bm_s + mu_s + 3.
389 cse(10) = bm_s + mu_s + bv_s + 1.
390 cse(11) = bm_s + mu_s + bv_s + 2.
391 cse(12) = bm_s + mu_s + bv_s + 3.
392 cse(13) = bv_s + 2.
393 cse(14) = bm_s + bv_s
394 cse(15) = bm_s - 1.
395 cse(16) = 1.0 + (1.0 + bv_s)/2.
396 cse(17) = bv_s + 3.
397 cse(18) = bv_s + mu_s + 3.
398 do n = 1, 18
399 csg(n) = WGAMMA(cse(n))
400 enddo
401 oams = 1./am_s
402 obms = 1./bm_s
403 ! sg_lam = csg(2) / csg(1)
404 ! sg_nos = csg(2)**cse(1) / csg(1)**cse(2)
405
406 cge(1) = bm_g + 1.
407 cge(2) = mu_g + 1.
408 cge(3) = bm_g + mu_g + 1.
409 cge(4) = bm_g + mu_g + 2.
410 cge(5) = bm_g + mu_g + 3.
411 cge(6) = bm_g + mu_g + bv_g + 1.
412 cge(7) = bm_g + mu_g + bv_g + 2.
413 cge(8) = bm_g + mu_g + bv_g + 3.
414 cge(9) = mu_g + bv_g + 3.
415 cge(10) = mu_g + 2.
416 cge(11) = 0.5*(bv_g + 5. + 2.*mu_g)
417 cge(12) = 0.5*(bv_g + 5.) + mu_g
418 do n = 1, 12
419 cgg(n) = WGAMMA(cge(n))
420 enddo
421 obmg = 1./bm_g
422 oge1 = 1./cge(1)
423 ogg1 = 1./cgg(1)
424 ogg2 = 1./cgg(2)
425 ogg3 = 1./cgg(3)
426
427 !+---+-----------------------------------------------------------------+
428 !..Simplify various rate eqns the best we can now.
429 !+---+-----------------------------------------------------------------+
430
431 !..Rain collecting cloud water and cloud ice
432 t1_qr_qc = PI*.25*av_r * crg(9)
433 t1_qr_qi = PI*.25*av_r * crg(9)
434 t2_qr_qi = PI*.25*am_r*av_r * crg(8)
435
436 !..Graupel collecting cloud water
437 t1_qg_qc = PI*.25*av_g * cgg(9)
438
439 !..Snow collecting cloud water
440 t1_qs_qc = PI*.25*av_s
441
442 !..Snow collecting cloud ice
443 t1_qs_qi = PI*.25*av_s
444
445 !..Evaporation of rain; ignore depositional growth of rain.
446 t1_qr_ev = 0.78 * crg(10)
447 t2_qr_ev = 0.308*Sc3*SQRT(av_r) * crg(11)
448
449 !..Sublimation/depositional growth of snow
450 t1_qs_sd = 0.86
451 t2_qs_sd = 0.28*Sc3*SQRT(av_s)
452
453 !..Melting of snow
454 t1_qs_me = PI*4.*C_sqrd*olfus * 0.86
455 t2_qs_me = PI*4.*C_sqrd*olfus * 0.28*Sc3*SQRT(av_s)
456
457 !..Sublimation/depositional growth of graupel
458 t1_qg_sd = 0.86 * cgg(10)
459 t2_qg_sd = 0.28*Sc3*SQRT(av_g) * cgg(11)
460
461 !..Melting of graupel
462 t1_qg_me = PI*4.*C_cube*olfus * 0.86 * cgg(10)
463 t2_qg_me = PI*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g) * cgg(11)
464
465 !..Constants for helping find lookup table indexes.
466 nic2 = NINT(ALOG10(r_c(1)))
467 nii2 = NINT(ALOG10(r_i(1)))
468 nii3 = NINT(ALOG10(Nt_i(1)))
469 nir2 = NINT(ALOG10(r_r(1)))
470 nir3 = NINT(ALOG10(N0r_exp(1)))
471 nis2 = NINT(ALOG10(r_s(1)))
472 nig2 = NINT(ALOG10(r_g(1)))
473
474 !+---+-----------------------------------------------------------------+
475 !..Create lookup tables for most costly calculations.
476 !+---+-----------------------------------------------------------------+
477
478 do k = 1, ntb_r
479 do j = 1, ntb_r1
480 do i = 1, ntb_g
481 tcg_racg(i,j,k) = 0.0d0
482 tmr_racg(i,j,k) = 0.0d0
483 tcr_gacr(i,j,k) = 0.0d0
484 tmg_gacr(i,j,k) = 0.0d0
485 enddo
486 enddo
487 enddo
488
489 do m = 1, ntb_r
490 do k = 1, ntb_r1
491 do j = 1, ntb_t
492 do i = 1, ntb_s
493 tcs_racs1(i,j,k,m) = 0.0d0
494 tmr_racs1(i,j,k,m) = 0.0d0
495 tcs_racs2(i,j,k,m) = 0.0d0
496 tmr_racs2(i,j,k,m) = 0.0d0
497 tcr_sacr1(i,j,k,m) = 0.0d0
498 tms_sacr1(i,j,k,m) = 0.0d0
499 tcr_sacr2(i,j,k,m) = 0.0d0
500 tms_sacr2(i,j,k,m) = 0.0d0
501 enddo
502 enddo
503 enddo
504 enddo
505
506 do k = 1, 45
507 do j = 1, ntb_r1
508 do i = 1, ntb_r
509 tpi_qrfz(i,j,k) = 0.0d0
510 tni_qrfz(i,j,k) = 0.0d0
511 tpg_qrfz(i,j,k) = 0.0d0
512 enddo
513 enddo
514 do i = 1, ntb_c
515 tpi_qcfz(i,k) = 0.0d0
516 tni_qcfz(i,k) = 0.0d0
517 enddo
518 enddo
519
520 do j = 1, ntb_i1
521 do i = 1, ntb_i
522 tps_iaus(i,j) = 0.0d0
523 tni_iaus(i,j) = 0.0d0
524 tpi_ide(i,j) = 0.0d0
525 enddo
526 enddo
527
528 if (.not. iiwarm) then
529 CALL wrf_debug(150, 'CREATING MICROPHYSICS LOOKUP TABLES ... ')
530 WRITE (wrf_err_message, '(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') &
531 ' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g
532 CALL wrf_debug(150, wrf_err_message)
533
534 !..Rain collecting graupel & graupel collecting rain.
535 CALL wrf_debug(200, ' creating rain collecting graupel table')
536 call qr_acr_qg
537
538 !..Rain collecting snow & snow collecting rain.
539 CALL wrf_debug(200, ' creating rain collecting snow table')
540 call qr_acr_qs
541
542 !..Cloud water and rain freezing (Bigg, 1953).
543 CALL wrf_debug(200, ' creating freezing of water drops table')
544 call freezeH2O
545
546 !..Conversion of some ice mass into snow category.
547 CALL wrf_debug(200, ' creating ice converting to snow table')
548 call qi_aut_qs
549
550 CALL wrf_debug(150, ' ... DONE microphysical lookup tables')
551 endif
552
553 END SUBROUTINE thompson_init
554 !+---+-----------------------------------------------------------------+
555 !
556 !+---+-----------------------------------------------------------------+
557 !..This is a wrapper routine designed to transfer values from 3D to 1D.
558 !+---+-----------------------------------------------------------------+
559 SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, &
560 th, pii, p, dz, dt_in, itimestep, &
561 RAINNC, RAINNCV, SR, &
562 ids,ide, jds,jde, kds,kde, & ! domain dims
563 ims,ime, jms,jme, kms,kme, & ! memory dims
564 its,ite, jts,jte, kts,kte) ! tile dims
565
566 implicit none
567
568 !..Subroutine arguments
569 INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, &
570 ims,ime, jms,jme, kms,kme, &
571 its,ite, jts,jte, kts,kte
572 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
573 qv, qc, qr, qi, qs, qg, ni, th
574 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
575 pii, p, dz
576 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: &
577 RAINNC, RAINNCV, SR
578 REAL, INTENT(IN):: dt_in
579 INTEGER, INTENT(IN):: itimestep
580
581 !..Local variables
582 REAL, DIMENSION(kts:kte):: &
583 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
584 t1d, p1d, dz1d
585 REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
586 REAL:: dt, pptrain, pptsnow, pptgraul, pptice
587 REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max
588 INTEGER:: i, j, k
589 INTEGER:: imax_qc, imax_qr, imax_qi, imax_qs, imax_qg, imax_ni
590 INTEGER:: jmax_qc, jmax_qr, jmax_qi, jmax_qs, jmax_qg, jmax_ni
591 INTEGER:: kmax_qc, kmax_qr, kmax_qi, kmax_qs, kmax_qg, kmax_ni
592 CHARACTER*256:: mp_debug
593
594 !+---+
595
596 dt = dt_in
597
598 qc_max = 0.
599 qr_max = 0.
600 qs_max = 0.
601 qi_max = 0.
602 qg_max = 0
603 ni_max = 0.
604 imax_qc = 0
605 imax_qr = 0
606 imax_qi = 0
607 imax_qs = 0
608 imax_qg = 0
609 imax_ni = 0
610 jmax_qc = 0
611 jmax_qr = 0
612 jmax_qi = 0
613 jmax_qs = 0
614 jmax_qg = 0
615 jmax_ni = 0
616 kmax_qc = 0
617 kmax_qr = 0
618 kmax_qi = 0
619 kmax_qs = 0
620 kmax_qg = 0
621 kmax_ni = 0
622 do i = 1, 256
623 mp_debug(i:i) = char(0)
624 enddo
625
626 j_loop: do j = jts, jte
627 i_loop: do i = its, ite
628
629 pptrain = 0.
630 pptsnow = 0.
631 pptgraul = 0.
632 pptice = 0.
633 RAINNCV(i,j) = 0.
634 SR(i,j) = 0.
635
636 do k = kts, kte
637 t1d(k) = th(i,k,j)*pii(i,k,j)
638 p1d(k) = p(i,k,j)
639 dz1d(k) = dz(i,k,j)
640 qv1d(k) = qv(i,k,j)
641 qc1d(k) = qc(i,k,j)
642 qi1d(k) = qi(i,k,j)
643 qr1d(k) = qr(i,k,j)
644 qs1d(k) = qs(i,k,j)
645 qg1d(k) = qg(i,k,j)
646 ni1d(k) = ni(i,k,j)
647 enddo
648
649 call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
650 t1d, p1d, dz1d, &
651 pptrain, pptsnow, pptgraul, pptice, &
652 kts, kte, dt, i, j)
653
654 pcp_ra(i,j) = pptrain
655 pcp_sn(i,j) = pptsnow
656 pcp_gr(i,j) = pptgraul
657 pcp_ic(i,j) = pptice
658 RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice
659 RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
660 SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12)
661
662 do k = kts, kte
663 qv(i,k,j) = qv1d(k)
664 qc(i,k,j) = qc1d(k)
665 qi(i,k,j) = qi1d(k)
666 qr(i,k,j) = qr1d(k)
667 qs(i,k,j) = qs1d(k)
668 qg(i,k,j) = qg1d(k)
669 ni(i,k,j) = ni1d(k)
670 th(i,k,j) = t1d(k)/pii(i,k,j)
671 if (qc1d(k) .gt. qc_max) then
672 imax_qc = i
673 jmax_qc = j
674 kmax_qc = k
675 qc_max = qc1d(k)
676 elseif (qc1d(k) .lt. 0.0) then
677 write(mp_debug,*) 'WARNING, negative qc ', qc1d(k), &
678 ' at i,j,k=', i,j,k
679 CALL wrf_debug(150, mp_debug)
680 endif
681 if (qr1d(k) .gt. qr_max) then
682 imax_qr = i
683 jmax_qr = j
684 kmax_qr = k
685 qr_max = qr1d(k)
686 elseif (qr1d(k) .lt. 0.0) then
687 write(mp_debug,*) 'WARNING, negative qr ', qr1d(k), &
688 ' at i,j,k=', i,j,k
689 CALL wrf_debug(150, mp_debug)
690 endif
691 if (qs1d(k) .gt. qs_max) then
692 imax_qs = i
693 jmax_qs = j
694 kmax_qs = k
695 qs_max = qs1d(k)
696 elseif (qs1d(k) .lt. 0.0) then
697 write(mp_debug,*) 'WARNING, negative qs ', qs1d(k), &
698 ' at i,j,k=', i,j,k
699 CALL wrf_debug(150, mp_debug)
700 endif
701 if (qi1d(k) .gt. qi_max) then
702 imax_qi = i
703 jmax_qi = j
704 kmax_qi = k
705 qi_max = qi1d(k)
706 elseif (qi1d(k) .lt. 0.0) then
707 write(mp_debug,*) 'WARNING, negative qi ', qi1d(k), &
708 ' at i,j,k=', i,j,k
709 CALL wrf_debug(150, mp_debug)
710 endif
711 if (qg1d(k) .gt. qg_max) then
712 imax_qg = i
713 jmax_qg = j
714 kmax_qg = k
715 qg_max = qg1d(k)
716 elseif (qg1d(k) .lt. 0.0) then
717 write(mp_debug,*) 'WARNING, negative qg ', qg1d(k), &
718 ' at i,j,k=', i,j,k
719 CALL wrf_debug(150, mp_debug)
720 endif
721 if (ni1d(k) .gt. ni_max) then
722 imax_ni = i
723 jmax_ni = j
724 kmax_ni = k
725 ni_max = ni1d(k)
726 elseif (ni1d(k) .lt. 0.0) then
727 write(mp_debug,*) 'WARNING, negative ni ', ni1d(k), &
728 ' at i,j,k=', i,j,k
729 CALL wrf_debug(150, mp_debug)
730 endif
731 enddo
732
733 enddo i_loop
734 enddo j_loop
735
736 ! DEBUG - GT
737 write(mp_debug,'(a,6(a,e13.6,1x,a,i3,a,i3,a,i3,a,1x))') 'MP-GT:', &
738 'qc: ', qc_max, '(', imax_qc, ',', jmax_qc, ',', kmax_qc, ')', &
739 'qr: ', qr_max, '(', imax_qr, ',', jmax_qr, ',', kmax_qr, ')', &
740 'qi: ', qi_max, '(', imax_qi, ',', jmax_qi, ',', kmax_qi, ')', &
741 'qs: ', qs_max, '(', imax_qs, ',', jmax_qs, ',', kmax_qs, ')', &
742 'qg: ', qg_max, '(', imax_qg, ',', jmax_qg, ',', kmax_qg, ')', &
743 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')'
744 CALL wrf_debug(150, mp_debug)
745 ! END DEBUG - GT
746
747 END SUBROUTINE mp_gt_driver
748
749 !+---+-----------------------------------------------------------------+
750 !
751 !+---+-----------------------------------------------------------------+
752 !+---+-----------------------------------------------------------------+
753 !.. This subroutine computes the moisture tendencies of water vapor,
754 !.. cloud droplets, rain, cloud ice (pristine), snow, and graupel.
755 !.. Previously this code was based on Reisner et al (1998), but few of
756 !.. those pieces remain. A complete description is now found in
757 !.. Thompson et al. (2004, 2006).
758 !+---+-----------------------------------------------------------------+
759 !
760 subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
761 t1d, p1d, dzq, &
762 pptrain, pptsnow, pptgraul, pptice, &
763 kts, kte, dt, ii, jj)
764
765 implicit none
766
767 !..Sub arguments
768 INTEGER, INTENT(IN):: kts, kte, ii, jj
769 REAL, DIMENSION(kts:kte), INTENT(INOUT):: &
770 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
771 t1d, p1d
772 REAL, DIMENSION(kts:kte), INTENT(IN):: dzq
773 REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice
774 REAL, INTENT(IN):: dt
775
776 !..Local variables
777 REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, &
778 qrten, qsten, qgten, niten
779
780 DOUBLE PRECISION, DIMENSION(kts:kte):: prw_vcd
781
782 DOUBLE PRECISION, DIMENSION(kts:kte):: prr_wau, prr_rcw, prr_rcs, &
783 prr_rcg, prr_sml, prr_gml, &
784 prr_rci, prv_rev
785
786 DOUBLE PRECISION, DIMENSION(kts:kte):: pri_inu, pni_inu, pri_ihm, &
787 pni_ihm, pri_wfz, pni_wfz, &
788 pri_rfz, pni_rfz, pri_ide, &
789 pni_ide, pri_rci, pni_rci, &
790 pni_sci, pni_iau
791
792 DOUBLE PRECISION, DIMENSION(kts:kte):: prs_iau, prs_sci, prs_rcs, &
793 prs_scw, prs_sde, prs_ihm, &
794 prs_ide
795
796 DOUBLE PRECISION, DIMENSION(kts:kte):: prg_scw, prg_rfz, prg_gde, &
797 prg_gcw, prg_rci, prg_rcs, &
798 prg_rcg, prg_ihm
799
800 REAL, DIMENSION(kts:kte):: temp, pres, qv
801 REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni
802 REAL, DIMENSION(kts:kte):: rho, rhof, rhof2
803 REAL, DIMENSION(kts:kte):: qvs, qvsi
804 REAL, DIMENSION(kts:kte):: satw, sati, ssatw, ssati
805 REAL, DIMENSION(kts:kte):: diffu, visco, vsc2, &
806 tcond, lvap, ocp, lvt2
807
808 DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g
809 REAL, DIMENSION(kts:kte):: mvd_r, mvd_c
810 REAL, DIMENSION(kts:kte):: smob, smo2, smo1, &
811 smoc, smod, smoe, smof
812
813 REAL, DIMENSION(kts:kte):: sed_r, sed_s, sed_g, sed_i, sed_n
814
815 REAL:: rgvm, delta_tp, orho, onstep, lfus2
816 DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamc, lamr, lamg
817 DOUBLE PRECISION:: lami, ilami
818 REAL:: Dc, Dc_b, Dc_g, Di, Dr, Ds, Dg, Ds_m, Dg_m
819 REAL:: zeta1, zeta, taud, tau
820 REAL:: stoke_r, stoke_s, stoke_g, stoke_i
821 REAL:: vti, vtr, vts, vtg
822 REAL, DIMENSION(kts:kte+1):: vtik, vtnk, vtrk, vtsk, vtgk
823 REAL, DIMENSION(kts:kte):: vts_boost
824 REAL:: Mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, C_snow
825 REAL:: a_, b_, loga_, A1, A2, tf
826 REAL:: tempc, tc0, r_mvd1, r_mvd2, xkrat
827 REAL:: xnc, xri, xni, xmi, oxmi, xrc
828 REAL:: xsat, rate_max, sump, ratio
829 REAL:: clap, fcd, dfcd
830 REAL:: otemp, rvs, rvs_p, rvs_pp, gamsc, alphsc, t1_evap, t1_subl
831 REAL:: r_frac, g_frac
832 REAL:: Ef_rw, Ef_ri, Ef_sw, Ef_gw
833 REAL:: dts, odts, odt, odzq
834 INTEGER:: i, k, k2, ksed1, ku, n, nn, nstep, k_0, kbot, IT, iexfrq
835 INTEGER:: nir, nis, nig, nii, nic
836 INTEGER:: idx_tc,idx_t,idx_s,idx_g,idx_r1,idx_r,idx_i1,idx_i,idx_c
837 LOGICAL:: melti, no_micro
838 LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg
839
840 !+---+
841
842 no_micro = .true.
843 dts = dt
844 odt = 1./dt
845 odts = 1./dts
846 iexfrq = 1
847
848 !+---+-----------------------------------------------------------------+
849 !.. Source/sink terms. First 2 chars: "pr" represents source/sink of
850 !.. mass while "pn" represents source/sink of number. Next char is one
851 !.. of "v" for water vapor, "r" for rain, "i" for cloud ice, "w" for
852 !.. cloud water, "s" for snow, and "g" for graupel. Next chars
853 !.. represent processes: "de" for sublimation/deposition, "ev" for
854 !.. evaporation, "fz" for freezing, "ml" for melting, "au" for
855 !.. autoconversion, "nu" for ice nucleation, "hm" for Hallet/Mossop
856 !.. secondary ice production, and "c" for collection followed by the
857 !.. character for the species being collected. ALL of these terms are
858 !.. positive (except for deposition/sublimation terms which can switch
859 !.. signs based on super/subsaturation) and are treated as negatives
860 !.. where necessary in the tendency equations.
861 !+---+-----------------------------------------------------------------+
862
863 do k = kts, kte
864 tten(k) = 0.
865 qvten(k) = 0.
866 qcten(k) = 0.
867 qiten(k) = 0.
868 qrten(k) = 0.
869 qsten(k) = 0.
870 qgten(k) = 0.
871 niten(k) = 0.
872
873 prw_vcd(k) = 0.
874
875 prv_rev(k) = 0.
876 prr_wau(k) = 0.
877 prr_rcw(k) = 0.
878 prr_rcs(k) = 0.
879 prr_rcg(k) = 0.
880 prr_sml(k) = 0.
881 prr_gml(k) = 0.
882 prr_rci(k) = 0.
883
884 pri_inu(k) = 0.
885 pni_inu(k) = 0.
886 pri_ihm(k) = 0.
887 pni_ihm(k) = 0.
888 pri_wfz(k) = 0.
889 pni_wfz(k) = 0.
890 pri_rfz(k) = 0.
891 pni_rfz(k) = 0.
892 pri_ide(k) = 0.
893 pni_ide(k) = 0.
894 pri_rci(k) = 0.
895 pni_rci(k) = 0.
896 pni_sci(k) = 0.
897 pni_iau(k) = 0.
898
899 prs_iau(k) = 0.
900 prs_sci(k) = 0.
901 prs_rcs(k) = 0.
902 prs_scw(k) = 0.
903 prs_sde(k) = 0.
904 prs_ihm(k) = 0.
905 prs_ide(k) = 0.
906
907 prg_scw(k) = 0.
908 prg_rfz(k) = 0.
909 prg_gde(k) = 0.
910 prg_gcw(k) = 0.
911 prg_rci(k) = 0.
912 prg_rcs(k) = 0.
913 prg_rcg(k) = 0.
914 prg_ihm(k) = 0.
915 enddo
916
917 !+---+-----------------------------------------------------------------+
918 !..Put column of data into local arrays.
919 !+---+-----------------------------------------------------------------+
920 do k = kts, kte
921 temp(k) = t1d(k)
922 qv(k) = MAX(1.E-10, qv1d(k))
923 pres(k) = p1d(k)
924 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
925 if (qc1d(k) .gt. R1) then
926 no_micro = .false.
927 rc(k) = qc1d(k)*rho(k)
928 L_qc(k) = .true.
929 else
930 qc1d(k) = 0.0
931 rc(k) = R1
932 L_qc(k) = .false.
933 endif
934 if (qi1d(k) .gt. R1) then
935 no_micro = .false.
936 ri(k) = qi1d(k)*rho(k)
937 ni(k) = MAX(1., ni1d(k)*rho(k))
938 L_qi(k) = .true.
939 else
940 qi1d(k) = 0.0
941 ni1d(k) = 0.0
942 ri(k) = R1
943 ni(k) = 0.01
944 L_qi(k) = .false.
945 endif
946 if (qr1d(k) .gt. R1) then
947 no_micro = .false.
948 rr(k) = qr1d(k)*rho(k)
949 L_qr(k) = .true.
950 else
951 qr1d(k) = 0.0
952 rr(k) = R1
953 L_qr(k) = .false.
954 endif
955 if (qs1d(k) .gt. R1) then
956 no_micro = .false.
957 rs(k) = qs1d(k)*rho(k)
958 L_qs(k) = .true.
959 else
960 qs1d(k) = 0.0
961 rs(k) = R1
962 L_qs(k) = .false.
963 endif
964 if (qg1d(k) .gt. R1) then
965 no_micro = .false.
966 rg(k) = qg1d(k)*rho(k)
967 L_qg(k) = .true.
968 else
969 qg1d(k) = 0.0
970 rg(k) = R1
971 L_qg(k) = .false.
972 endif
973 enddo
974
975
976 !+---+-----------------------------------------------------------------+
977 !..Derive various thermodynamic variables frequently used.
978 !.. Saturation vapor pressure (mixing ratio) over liquid/ice comes from
979 !.. Flatau et al. 1992; enthalpy (latent heat) of vaporization from
980 !.. Bohren & Albrecht 1998; others from Pruppacher & Klett 1978.
981 !+---+-----------------------------------------------------------------+
982 do k = kts, kte
983 tempc = temp(k) - 273.15
984 rhof(k) = SQRT(RHO_NOT/rho(k))
985 rhof2(k) = SQRT(rhof(k))
986 qvs(k) = rslf(pres(k), temp(k))
987 if (tempc .le. 0.0) then
988 qvsi(k) = rsif(pres(k), temp(k))
989 else
990 qvsi(k) = qvs(k)
991 endif
992 satw(k) = qv(k)/qvs(k)
993 sati(k) = qv(k)/qvsi(k)
994 ssatw(k) = satw(k) - 1.
995 ssati(k) = sati(k) - 1.
996 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
997 if (abs(ssati(k)).lt. eps) ssati(k) = 0.0
998 if (no_micro .and. ssati(k).gt. 0.0) no_micro = .false.
999 diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
1000 if (tempc .ge. 0.0) then
1001 visco(k) = (1.718+0.0049*tempc)*1.0E-5
1002 else
1003 visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
1004 endif
1005 ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
1006 vsc2(k) = SQRT(rho(k)/visco(k))
1007 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
1008 tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
1009 enddo
1010
1011 !+---+-----------------------------------------------------------------+
1012 !..If no existing hydrometeor species and no chance to initiate ice or
1013 !.. condense cloud water, just exit quickly!
1014 !+---+-----------------------------------------------------------------+
1015
1016 if (no_micro) return
1017
1018 !+---+-----------------------------------------------------------------+
1019 !..Calculate y-intercept, slope, and useful moments for snow.
1020 !+---+-----------------------------------------------------------------+
1021 if (.not. iiwarm) then
1022 do k = kts, kte
1023 if (.not. L_qs(k)) CYCLE
1024 tc0 = MIN(-0.1, temp(k)-273.15)
1025 smob(k) = rs(k)*oams
1026
1027 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
1028 !.. then we must compute actual 2nd moment and use as reference.
1029 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
1030 smo2(k) = smob(k)
1031 else
1032 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
1033 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
1034 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
1035 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
1036 + sa(10)*bm_s*bm_s*bm_s
1037 a_ = 10.0**loga_
1038 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
1039 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
1040 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
1041 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
1042 + sb(10)*bm_s*bm_s*bm_s
1043 smo2(k) = (smob(k)/a_)**(1./b_)
1044 endif
1045
1046 !..Calculate 1st moment. Useful for depositional growth and melting.
1047 loga_ = sa(1) + sa(2)*tc0 + sa(3) &
1048 + sa(4)*tc0 + sa(5)*tc0*tc0 &
1049 + sa(6) + sa(7)*tc0*tc0 &
1050 + sa(8)*tc0 + sa(9)*tc0*tc0*tc0 &
1051 + sa(10)
1052 a_ = 10.0**loga_
1053 b_ = sb(1)+ sb(2)*tc0 + sb(3) + sb(4)*tc0 &
1054 + sb(5)*tc0*tc0 + sb(6) &
1055 + sb(7)*tc0*tc0 + sb(8)*tc0 &
1056 + sb(9)*tc0*tc0*tc0 + sb(10)
1057 smo1(k) = a_ * smo2(k)**b_
1058
1059 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
1060 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
1061 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
1062 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
1063 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
1064 + sa(10)*cse(1)*cse(1)*cse(1)
1065 a_ = 10.0**loga_
1066 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
1067 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
1068 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
1069 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
1070 smoc(k) = a_ * smo2(k)**b_
1071
1072 !..Calculate bv_s+2 (th) moment. Useful for riming.
1073 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(13) &
1074 + sa(4)*tc0*cse(13) + sa(5)*tc0*tc0 &
1075 + sa(6)*cse(13)*cse(13) + sa(7)*tc0*tc0*cse(13) &
1076 + sa(8)*tc0*cse(13)*cse(13) + sa(9)*tc0*tc0*tc0 &
1077 + sa(10)*cse(13)*cse(13)*cse(13)
1078 a_ = 10.0**loga_
1079 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(13) + sb(4)*tc0*cse(13) &
1080 + sb(5)*tc0*tc0 + sb(6)*cse(13)*cse(13) &
1081 + sb(7)*tc0*tc0*cse(13) + sb(8)*tc0*cse(13)*cse(13) &
1082 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(13)*cse(13)*cse(13)
1083 smoe(k) = a_ * smo2(k)**b_
1084
1085 !..Calculate 1+(bv_s+1)/2 (th) moment. Useful for depositional growth.
1086 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(16) &
1087 + sa(4)*tc0*cse(16) + sa(5)*tc0*tc0 &
1088 + sa(6)*cse(16)*cse(16) + sa(7)*tc0*tc0*cse(16) &
1089 + sa(8)*tc0*cse(16)*cse(16) + sa(9)*tc0*tc0*tc0 &
1090 + sa(10)*cse(16)*cse(16)*cse(16)
1091 a_ = 10.0**loga_
1092 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(16) + sb(4)*tc0*cse(16) &
1093 + sb(5)*tc0*tc0 + sb(6)*cse(16)*cse(16) &
1094 + sb(7)*tc0*tc0*cse(16) + sb(8)*tc0*cse(16)*cse(16) &
1095 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(16)*cse(16)*cse(16)
1096 smof(k) = a_ * smo2(k)**b_
1097
1098 enddo
1099
1100 !+---+-----------------------------------------------------------------+
1101 !..Calculate y-intercept, slope values for graupel.
1102 !+---+-----------------------------------------------------------------+
1103 N0_min = gonv_max
1104 do k = kte, kts, -1
1105 if (.not. L_qg(k)) CYCLE
1106 N0_exp = 100.0*rho(k)/rg(k)
1107 N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max)))
1108 N0_min = MIN(N0_exp, N0_min)
1109 N0_exp = N0_min
1110 lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
1111 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
1112 ilamg(k) = 1./lamg
1113 N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2)
1114 enddo
1115
1116 endif
1117
1118 !+---+-----------------------------------------------------------------+
1119 !..Calculate y-intercept & slope values for rain.
1120 !.. New treatment for variable y-intercept of rain. When rain comes
1121 !.. from melted snow/graupel, compute mass-weighted mean size, melt
1122 !.. into water, compute its mvd and recompute slope/intercept.
1123 !.. If rain not from melted snow, use old relation but hold N0_r
1124 !.. constant at its lowest value. While doing all this, ensure rain
1125 !.. mvd does not exceed reasonable size like 3 mm.
1126 !+---+-----------------------------------------------------------------+
1127 N0_min = ronv_max
1128 do k = kte, kts, -1
1129 ! if (.not. L_qr(k)) CYCLE
1130 N0_exp = ronv_c1*tanh(ronv_c0*(ronv_r0-rr(k))) + ronv_c2
1131 N0_min = MIN(N0_exp, N0_min)
1132 N0_exp = N0_min
1133 lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
1134 lamr = lam_exp * (crg(3)*org2*org1)**obmr
1135 mvd_r(k) = (3.0+mu_r+0.672) / lamr
1136 if (mvd_r(k) .gt. 3.e-3) then
1137 mvd_r(k) = 3.e-3
1138 lamr = (3.0+mu_r+0.672) / 3.e-3
1139 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
1140 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
1141 endif
1142 N0_r(k) = N0_exp/(crg(2)*lam_exp) * lamr**cre(2)
1143 ilamr(k) = 1./lamr
1144 enddo
1145
1146 if (.not. iiwarm) then
1147 k_0 = kts
1148 melti = .false.
1149 do k = kte-1, kts, -1
1150 if ( (temp(k).gt. T_0) .and. (rr(k).gt. 0.001e-3) &
1151 .and. ((rs(k+1)+rg(k+1)).gt. 0.01e-3) ) then
1152 k_0 = MAX(k+1, k_0)
1153 melti=.true.
1154 goto 135
1155 endif
1156 enddo
1157 135 continue
1158
1159 if (melti) then
1160 !.. Locate bottom of melting layer (if any).
1161 kbot = kts
1162 do k = k_0-1, kts, -1
1163 if ( (rs(k)+rg(k)).lt. 0.01e-3) goto 136
1164 enddo
1165 136 continue
1166 kbot = MAX(k, kts)
1167
1168 !.. Compute melted snow/graupel equiv water diameter one K-level above
1169 !.. melting. Set starting rain mvd to either 50 microns or max from
1170 !.. higher up in column.
1171 if (L_qs(k_0)) then
1172 Ds = smoc(k_0) / smob(k_0)
1173 Ds_m = (am_s*Ds**bm_s / am_r)**obmr
1174 else
1175 Ds_m = 1.0e-6
1176 endif
1177 if (L_qg(k_0)) then
1178 Dg = (bm_g + mu_g + 1.) * ilamg(k_0)
1179 Dg_m = (am_g*Dg**bm_g / am_r)**obmr
1180 else
1181 Dg_m = 1.0e-6
1182 endif
1183 r_mvd1 = mvd_r(k_0)
1184 r_mvd2 = MIN(MAX(Ds_m, Dg_m, r_mvd1+1.e-6, mvd_r(kbot)), &
1185 3.e-3)
1186
1187 !.. Within melting layer, apply linear increase of rain mvd from r_mvd1
1188 !.. to equiv melted snow/graupel value (r_mvd2). So, by the bottom of
1189 !.. the melting layer, the rain will have an mvd that matches that from
1190 !.. melted snow and/or graupel.
1191 if (kbot.gt. 2) then
1192 do k = k_0-1, kbot, -1
1193 if (.not. L_qr(k)) CYCLE
1194 xkrat = REAL(k_0-k)/REAL(k_0-kbot)
1195 mvd_r(k) = MAX(mvd_r(k), xkrat*(r_mvd2-r_mvd1)+r_mvd1)
1196 lamr = (4.0+mu_r) / mvd_r(k)
1197 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
1198 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
1199 N0_exp = MAX(DBLE(ronv_min), MIN(N0_exp, DBLE(ronv_max)))
1200 lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
1201 lamr = lam_exp * (crg(3)*org2*org1)**obmr
1202 mvd_r(k) = (3.0+mu_r+0.672) / lamr
1203 N0_r(k) = rr(k)*lamr**cre(3) / (am_r*crg(3))
1204 ilamr(k) = 1./lamr
1205 enddo
1206
1207 !.. Below melting layer, hold N0_r constant unless changes to mixing
1208 !.. ratio increase mvd beyond 3 mm threshold, then adjust slope and
1209 !.. intercept to cap mvd at 3 mm. In future, we could lower N0_r to
1210 !.. account for self-collection or other sinks.
1211 do k = kbot-1, kts, -1
1212 if (.not. L_qr(k)) CYCLE
1213 N0_r(k) = MIN(N0_r(k), N0_r(kbot))
1214 lamr = (N0_r(k)*am_r*crg(3)/rr(k))**(1./cre(3))
1215 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
1216 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
1217 N0_exp = MAX(DBLE(ronv_min), MIN(N0_exp, DBLE(ronv_max)))
1218 lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
1219 lamr = lam_exp * (crg(3)*org2*org1)**obmr
1220 mvd_r(k) = (3.0+mu_r+0.672) / lamr
1221 if (mvd_r(k) .gt. 3.e-3) then
1222 mvd_r(k) = 3.e-3
1223 lamr = (3.0+mu_r+0.672) / mvd_r(k)
1224 endif
1225 N0_r(k) = rr(k)*lamr**cre(3) / (am_r*crg(3))
1226 ilamr(k) = 1./lamr
1227 enddo
1228 endif
1229
1230 endif
1231 endif
1232
1233 !+---+-----------------------------------------------------------------+
1234 !..Compute warm-rain process terms (except evap done later).
1235 !+---+-----------------------------------------------------------------+
1236
1237 do k = kts, kte
1238 if (.not. L_qc(k)) CYCLE
1239 Dc = MAX(D0c*1.E6, ((rc(k)/(am_r*Nt_c))**obmr) * 1.E6)
1240 lamc = (Nt_c*am_r* ccg(2) * ocg1 / rc(k))**obmr
1241 mvd_c(k) = (3.0+mu_c+0.672) / lamc
1242
1243 !..Autoconversion follows Berry & Reinhardt (1974) with characteristic
1244 !.. diameters correctly computed from gamma distrib of cloud droplets.
1245 if (rc(k).gt. 0.01e-3) then
1246 Dc_g = ((ccg(3)*ocg2)**obmr / lamc) * 1.E6
1247 Dc_b = (Dc*Dc*Dc*Dc_g*Dc_g*Dc_g - Dc*Dc*Dc*Dc*Dc*Dc)**(1./6.)
1248 zeta1 = 0.5*((6.25E-6*Dc*Dc_b*Dc_b*Dc_b - 0.4) &
1249 + abs(6.25E-6*Dc*Dc_b*Dc_b*Dc_b - 0.4))
1250 zeta = 0.027*rc(k)*zeta1
1251 taud = 0.5*((0.5*Dc_b - 7.5) + abs(0.5*Dc_b - 7.5)) + R1
1252 tau = 3.72/(rc(k)*taud)
1253 prr_wau(k) = zeta/tau
1254 prr_wau(k) = MIN(DBLE(rc(k)*odts), prr_wau(k))
1255 endif
1256
1257 !..Rain collecting cloud water. In CE, assume Dc<<Dr and vtc=~0.
1258 if (L_qr(k) .and. mvd_r(k).gt. D0r .and. mvd_c(k).gt. D0c) then
1259 Dr = (bm_r + mu_r + 1.) * ilamr(k)
1260 lamr = 1./ilamr(k)
1261 vtr = rhof(k)*av_r*crg(6)*org3 * (lamr/(lamr+fv_r))**cre(3) &
1262 / (lamr+fv_r)**bv_r
1263 stoke_r = mvd_c(k)*mvd_c(k)*vtr*rho_w/(9.*visco(k)*Dr)
1264 Ef_rw = stoke_r*stoke_r/((stoke_r+0.5)*(stoke_r+0.5))
1265 Ef_rw = MAX(0.0, MIN(0.99, Ef_rw))
1266 prr_rcw(k) = rhof(k)*t1_qr_qc*Ef_rw*rc(k)*N0_r(k) &
1267 *((lamr+fv_r)**(-cre(9)))
1268 prr_rcw(k) = MIN(DBLE(rc(k)*odts), prr_rcw(k))
1269 endif
1270 enddo
1271
1272 !+---+-----------------------------------------------------------------+
1273 !..Compute all frozen hydrometeor species' process terms.
1274 !+---+-----------------------------------------------------------------+
1275 if (.not. iiwarm) then
1276 do k = kts, kte
1277 vts_boost(k) = 1.5
1278
1279 !..Temperature lookup table indexes.
1280 tempc = temp(k) - 273.15
1281 idx_tc = MAX(1, MIN(NINT(-tempc), 45) )
1282 idx_t = INT( (tempc-2.5)/5. ) - 1
1283 idx_t = MAX(1, -idx_t)
1284 idx_t = MIN(idx_t, ntb_t)
1285 IT = MAX(1, MIN(NINT(-tempc), 31) )
1286
1287 !..Cloud water lookup table index.
1288 if (rc(k).gt. r_c(1)) then
1289 nic = NINT(ALOG10(rc(k)))
1290 do nn = nic-1, nic+1
1291 n = nn
1292 if ( (rc(k)/10.**nn).ge.1.0 .and. &
1293 (rc(k)/10.**nn).lt.10.0) goto 141
1294 enddo
1295 141 continue
1296 idx_c = INT(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
1297 idx_c = MAX(1, MIN(idx_c, ntb_c))
1298 else
1299 idx_c = 1
1300 endif
1301
1302 !..Cloud ice lookup table indexes.
1303 if (ri(k).gt. r_i(1)) then
1304 nii = NINT(ALOG10(ri(k)))
1305 do nn = nii-1, nii+1
1306 n = nn
1307 if ( (ri(k)/10.**nn).ge.1.0 .and. &
1308 (ri(k)/10.**nn).lt.10.0) goto 142
1309 enddo
1310 142 continue
1311 idx_i = INT(ri(k)/10.**n) + 10*(n-nii2) - (n-nii2)
1312 idx_i = MAX(1, MIN(idx_i, ntb_i))
1313 else
1314 idx_i = 1
1315 endif
1316
1317 if (ni(k).gt. Nt_i(1)) then
1318 nii = NINT(ALOG10(ni(k)))
1319 do nn = nii-1, nii+1
1320 n = nn
1321 if ( (ni(k)/10.**nn).ge.1.0 .and. &
1322 (ni(k)/10.**nn).lt.10.0) goto 143
1323 enddo
1324 143 continue
1325 idx_i1 = INT(ni(k)/10.**n) + 10*(n-nii3) - (n-nii3)
1326 idx_i1 = MAX(1, MIN(idx_i1, ntb_i1))
1327 else
1328 idx_i1 = 1
1329 endif
1330
1331 !..Rain lookup table indexes.
1332 if (rr(k).gt. r_r(1)) then
1333 nir = NINT(ALOG10(rr(k)))
1334 do nn = nir-1, nir+1
1335 n = nn
1336 if ( (rr(k)/10.**nn).ge.1.0 .and. &
1337 (rr(k)/10.**nn).lt.10.0) goto 144
1338 enddo
1339 144 continue
1340 idx_r = INT(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2)
1341 idx_r = MAX(1, MIN(idx_r, ntb_r))
1342
1343 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
1344 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
1345 nir = NINT(DLOG10(N0_exp))
1346 do nn = nir-1, nir+1
1347 n = nn
1348 if ( (N0_exp/10.**nn).ge.1.0 .and. &
1349 (N0_exp/10.**nn).lt.10.0) goto 145
1350 enddo
1351 145 continue
1352 idx_r1 = INT(N0_exp/10.**n) + 10*(n-nir3) - (n-nir3)
1353 idx_r1 = MAX(1, MIN(idx_r1, ntb_r1))
1354 else
1355 idx_r = 1
1356 idx_r1 = ntb_r1
1357 endif
1358
1359 !..Snow lookup table index.
1360 if (rs(k).gt. r_s(1)) then
1361 nis = NINT(ALOG10(rs(k)))
1362 do nn = nis-1, nis+1
1363 n = nn
1364 if ( (rs(k)/10.**nn).ge.1.0 .and. &
1365 (rs(k)/10.**nn).lt.10.0) goto 146
1366 enddo
1367 146 continue
1368 idx_s = INT(rs(k)/10.**n) + 10*(n-nis2) - (n-nis2)
1369 idx_s = MAX(1, MIN(idx_s, ntb_s))
1370 else
1371 idx_s = 1
1372 endif
1373
1374 !..Graupel lookup table index.
1375 if (rg(k).gt. r_g(1)) then
1376 nig = NINT(ALOG10(rg(k)))
1377 do nn = nig-1, nig+1
1378 n = nn
1379 if ( (rg(k)/10.**nn).ge.1.0 .and. &
1380 (rg(k)/10.**nn).lt.10.0) goto 147
1381 enddo
1382 147 continue
1383 idx_g = INT(rg(k)/10.**n) + 10*(n-nig2) - (n-nig2)
1384 idx_g = MAX(1, MIN(idx_g, ntb_g))
1385 else
1386 idx_g = 1
1387 endif
1388
1389 !..Deposition/sublimation prefactor (from Srivastava & Coen 1992).
1390 otemp = 1./temp(k)
1391 rvs = rho(k)*qvsi(k)
1392 rvs_p = rvs*otemp*(lsub*otemp*oRv - 1.)
1393 rvs_pp = rvs * ( otemp*(lsub*otemp*oRv - 1.) &
1394 *otemp*(lsub*otemp*oRv - 1.) &
1395 + (-2.*lsub*otemp*otemp*otemp*oRv) &
1396 + otemp*otemp)
1397 gamsc = lsub*diffu(k)/tcond(k) * rvs_p
1398 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
1399 * rvs_pp/rvs_p * rvs/rvs_p
1400 alphsc = MAX(1.E-9, alphsc)
1401 xsat = ssati(k)
1402 if (abs(xsat).lt. 1.E-9) xsat=0.
1403 t1_subl = 4.*PI*( 1.0 - alphsc*xsat &
1404 + 2.*alphsc*alphsc*xsat*xsat &
1405 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
1406 / (1.+gamsc)
1407
1408 !..Snow collecting cloud water. In CE, assume Dc<<Ds and vtc=~0.
1409 if (L_qc(k) .and. mvd_c(k).gt. D0c) then
1410 Ds = 0.0
1411 if (L_qs(k)) Ds = smoc(k) / smob(k)
1412 if (Ds .gt. D0s) then
1413 Mrat = 1./Ds
1414 ils1 = 1./(Mrat*Lam0 + fv_s)
1415 ils2 = 1./(Mrat*Lam1 + fv_s)
1416 t1_vts = Kap0*csg(4)*ils1**cse(4)
1417 t2_vts = Kap1*Mrat**mu_s*csg(10)*ils2**cse(10)
1418 t3_vts = Kap0*csg(1)*ils1**cse(1)
1419 t4_vts = Kap1*Mrat**mu_s*csg(7)*ils2**cse(7)
1420 vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts)
1421 stoke_s = mvd_c(k)*mvd_c(k)*vts*rho_w/(9.*visco(k)*Ds)
1422 Ef_sw = stoke_s*stoke_s/((stoke_s+0.4)*(stoke_s+0.4))
1423 Ef_sw = MAX(0.0, MIN(0.99, Ef_sw))
1424 prs_scw(k) = rhof(k)*t1_qs_qc*Ef_sw*rc(k)*smoe(k)
1425 endif
1426
1427 !..Graupel collecting cloud water. In CE, assume Dc<<Dg and vtc=~0.
1428 if (rg(k).ge. r_g(1) .and. mvd_c(k).gt. D0c) then
1429 Dg = (bm_g + mu_g + 1.) * ilamg(k)
1430 vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
1431 stoke_g = mvd_c(k)*mvd_c(k)*vtg*rho_w/(9.*visco(k)*Dg)
1432 Ef_gw = stoke_g*stoke_g/((stoke_g+0.5)*(stoke_g+0.5))
1433 Ef_gw = MAX(0.0, MIN(0.99, Ef_gw))
1434 prg_gcw(k) = rhof(k)*t1_qg_qc*Ef_gw*rc(k)*N0_g(k) &
1435 *ilamg(k)**cge(9)
1436 endif
1437 endif
1438
1439 !..Rain collecting snow. Cannot assume Wisner (1972) approximation
1440 !.. or Mizuno (1990) approach so we solve the CE explicitly and store
1441 !.. results in lookup table.
1442 if (rr(k).ge. r_r(1)) then
1443 if (rs(k).ge. r_s(1)) then
1444 if (temp(k).lt.T_0) then
1445 prr_rcs(k) = -(tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
1446 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
1447 + tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
1448 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r))
1449 prs_rcs(k) = tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
1450 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
1451 - tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
1452 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
1453 prg_rcs(k) = tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
1454 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
1455 + tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
1456 + tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
1457 prr_rcs(k) = MAX(DBLE(-rr(k)*odts), prr_rcs(k))
1458 prs_rcs(k) = MAX(DBLE(-rs(k)*odts), prs_rcs(k))
1459 prg_rcs(k) = MIN(DBLE((rr(k)+rs(k))*odts), prg_rcs(k))
1460 else
1461 prs_rcs(k) = -(tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
1462 + tcs_racs2(idx_s,idx_t,idx_r1,idx_r))
1463 prs_rcs(k) = MAX(DBLE(-rs(k)*odts), prs_rcs(k))
1464 prr_rcs(k) = -prs_rcs(k)
1465 endif
1466 endif
1467
1468 !..Rain collecting graupel. Cannot assume Wisner (1972) approximation
1469 !.. or Mizuno (1990) approach so we solve the CE explicitly and store
1470 !.. results in lookup table.
1471 if (rg(k).ge. r_g(1)) then
1472 if (temp(k).lt.T_0) then
1473 prg_rcg(k) = tmr_racg(idx_g,idx_r1,idx_r) &
1474 + tcr_gacr(idx_g,idx_r1,idx_r)
1475 prg_rcg(k) = MIN(DBLE(rr(k)*odts), prg_rcg(k))
1476 prr_rcg(k) = -prg_rcg(k)
1477 else
1478 prr_rcg(k) = tcg_racg(idx_g,idx_r1,idx_r)
1479 prr_rcg(k) = MIN(DBLE(rg(k)*odts), prr_rcg(k))
1480 prg_rcg(k) = -prr_rcg(k)
1481 endif
1482 endif
1483 endif
1484
1485 !+---+-----------------------------------------------------------------+
1486 !..Next IF block handles only those processes below 0C.
1487 !+---+-----------------------------------------------------------------+
1488
1489 if (temp(k).lt.T_0) then
1490
1491 vts_boost(k) = 1.0
1492 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
1493
1494 !..Freezing of water drops into graupel/cloud ice (Bigg 1953).
1495 if (rr(k).gt. r_r(1)) then
1496 prg_rfz(k) = tpg_qrfz(idx_r,idx_r1,idx_tc)*odts
1497 pri_rfz(k) = tpi_qrfz(idx_r,idx_r1,idx_tc)*odts
1498 pni_rfz(k) = tni_qrfz(idx_r,idx_r1,idx_tc)*odts
1499 elseif (rr(k).gt. R1 .and. temp(k).lt.HGFR) then
1500 pri_rfz(k) = rr(k)*odts
1501 pni_rfz(k) = N0_r(k)*crg(2) * ilamr(k)**cre(2)
1502 endif
1503 if (rc(k).gt. r_c(1)) then
1504 pri_wfz(k) = tpi_qcfz(idx_c,idx_tc)*odts
1505 pri_wfz(k) = MIN(DBLE(rc(k)*odts), pri_wfz(k))
1506 pni_wfz(k) = tni_qcfz(idx_c,idx_tc)*odts
1507 pni_wfz(k) = MIN(DBLE(Nt_c*odts), pri_wfz(k)/(2.*xm0i), &
1508 pni_wfz(k))
1509 endif
1510
1511 !..Nucleate ice from deposition & condensation freezing (Cooper 1986)
1512 !.. but only if water or 8%+ ice supersaturated.
1513 if ( (ssati(k).ge. 0.08) .or. (ssatw(k).gt. eps &
1514 .and. temp(k).lt.265.15) ) then
1515 xnc = MIN(5.E5, TNO*EXP(ATO*(T_0-temp(k))))
1516 xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dts
1517 pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts
1518 pri_inu(k) = MIN(DBLE(rate_max), xm0i*pni_inu(k))
1519 pni_inu(k) = pri_inu(k)/xm0i
1520 endif
1521
1522 !..Deposition/sublimation of cloud ice (Srivastava & Coen 1992).
1523 if (L_qi(k)) then
1524 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
1525 ilami = 1./lami
1526 Di = MAX(DBLE(D0i), (bm_i + mu_i + 1.) * ilami)
1527 xmi = am_i*Di**bm_i
1528 oxmi = 1./xmi
1529 pri_ide(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
1530 *oig1*cig(5)*ni(k)*ilami
1531
1532 if (pri_ide(k) .lt. 0.0) then
1533 pri_ide(k) = MAX(DBLE(-ri(k)*odts), pri_ide(k), DBLE(rate_max))
1534 pni_ide(k) = pri_ide(k)*oxmi
1535 pni_ide(k) = MAX(DBLE(-ni(k)*odts), pni_ide(k))
1536 else
1537 pri_ide(k) = MIN(pri_ide(k), DBLE(rate_max))
1538 prs_ide(k) = (1.0D0-tpi_ide(idx_i,idx_i1))*pri_ide(k)
1539 pri_ide(k) = tpi_ide(idx_i,idx_i1)*pri_ide(k)
1540 endif
1541
1542 !..Some cloud ice needs to move into the snow category. Use lookup
1543 !.. table that resulted from explicit bin representation of distrib.
1544 if ( (idx_i.eq. ntb_i) .or. (Di.gt. 5.0*D0s) ) then
1545 prs_iau(k) = ri(k)*.99*odts
1546 pni_iau(k) = ni(k)*.95*odts
1547 elseif (Di.lt. 0.1*D0s) then
1548 prs_iau(k) = 0.
1549 pni_iau(k) = 0.
1550 else
1551 prs_iau(k) = tps_iaus(idx_i,idx_i1)*odts
1552 prs_iau(k) = MIN(DBLE(ri(k)*.99*odts), prs_iau(k))
1553 pni_iau(k) = tni_iaus(idx_i,idx_i1)*odts
1554 pni_iau(k) = MIN(DBLE(ni(k)*.95*odts), pni_iau(k))
1555 endif
1556 endif
1557
1558 !..Deposition/sublimation of snow/graupel follows Srivastava & Coen
1559 !.. (1992).
1560 if (L_qs(k)) then
1561 C_snow = C_sqrd + (tempc+15.)*(C_cube-C_sqrd)/(-30.+15.)
1562 C_snow = MAX(C_sqrd, MIN(C_snow, C_cube))
1563 prs_sde(k) = C_snow*t1_subl*diffu(k)*ssati(k)*rvs &
1564 * (t1_qs_sd*smo1(k) &
1565 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
1566 if (prs_sde(k).lt. 0.) then
1567 prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k), DBLE(rate_max))
1568 else
1569 prs_sde(k) = MIN(prs_sde(k), DBLE(rate_max))
1570 endif
1571 endif
1572
1573 if (L_qg(k) .and. ssati(k).lt. -eps) then
1574 prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
1575 * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
1576 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
1577 if (prg_gde(k).lt. 0.) then
1578 prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k), DBLE(rate_max))
1579 else
1580 prg_gde(k) = MIN(prg_gde(k), DBLE(rate_max))
1581 endif
1582 endif
1583
1584 !..Snow collecting cloud ice. In CE, assume Di<<Ds and vti=~0.
1585 if (L_qi(k)) then
1586 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
1587 ilami = 1./lami
1588 Di = MAX(DBLE(D0i), (bm_i + mu_i + 1.) * ilami)
1589 xmi = am_i*Di**bm_i
1590 oxmi = 1./xmi
1591 if (rs(k).ge. r_s(1)) then
1592 prs_sci(k) = t1_qs_qi*rhof(k)*Ef_si*ri(k)*smoe(k)
1593 pni_sci(k) = prs_sci(k) * oxmi
1594 endif
1595
1596 !..Rain collecting cloud ice. In CE, assume Di<<Dr and vti=~0.
1597 if (rr(k).ge. r_r(1) .and. mvd_r(k).gt. 4.*Di) then
1598 lamr = 1./ilamr(k)
1599 Dr = (bm_r + mu_r + 1.) * ilamr(k)
1600 lamr = 1./ilamr(k)
1601 vtr = rhof(k)*av_r*crg(6)*org3*(lamr/(lamr+fv_r))**cre(3) &
1602 /(lamr+fv_r)**bv_r
1603 stoke_i = Di*Di*vtr*rho_i/(9.*visco(k)*Dr)
1604 Ef_ri = stoke_i*stoke_i/((stoke_i+0.5)*(stoke_i+0.5))
1605 Ef_ri = MAX(0.0, MIN(0.99, Ef_ri))
1606 pri_rci(k) = rhof(k)*t1_qr_qi*Ef_ri*ri(k)*N0_r(k) &
1607 *((lamr+fv_r)**(-cre(9)))
1608 pni_rci(k) = pri_rci(k) * oxmi
1609 prr_rci(k) = rhof(k)*t2_qr_qi*Ef_ri*ni(k)*N0_r(k) &
1610 *((lamr+fv_r)**(-cre(8)))
1611 prr_rci(k) = MIN(DBLE(rr(k)*odts), prr_rci(k))
1612 prg_rci(k) = pri_rci(k) + prr_rci(k)
1613 endif
1614 endif
1615
1616 !..Ice multiplication from rime-splinters (Hallet & Mossop 1974).
1617 if (prg_gcw(k).gt. eps .and. tempc.gt.-8.0) then
1618 tf = 0.
1619 if (tempc.ge.-5.0 .and. tempc.lt.-3.0) then
1620 tf = 0.5*(-3.0 - tempc)
1621 elseif (tempc.gt.-8.0 .and. tempc.lt.-5.0) then
1622 tf = 0.33333333*(8.0 + tempc)
1623 endif
1624 pni_ihm(k) = 3.5E8*tf*prg_gcw(k)
1625 pri_ihm(k) = xm0i*pni_ihm(k)
1626 prs_ihm(k) = prs_scw(k)/(prs_scw(k)+prg_gcw(k)) &
1627 * pri_ihm(k)
1628 prg_ihm(k) = prg_gcw(k)/(prs_scw(k)+prg_gcw(k)) &
1629 * pri_ihm(k)
1630 endif
1631
1632 !..A portion of rimed snow converts to graupel but some remains snow.
1633 !.. Interp from 5 to 75% as riming factor increases from 5.0 to 30.0
1634 !.. 0.028 came from (.75-.05)/(30.-5.). This remains ad-hoc and should
1635 !.. be revisited.
1636 if (prs_scw(k).gt.5.0*prs_sde(k) .and. &
1637 prs_sde(k).gt.eps) then
1638 r_frac = MIN(30.0D0, prs_scw(k)/prs_sde(k))
1639 g_frac = MIN(0.75, 0.05 + (r_frac-5.)*.028)
1640 vts_boost(k) = MIN(1.5, 1.1 + (r_frac-5.)*.016)
1641 prg_scw(k) = g_frac*prs_scw(k)
1642 prs_scw(k) = (1. - g_frac)*prs_scw(k)
1643 endif
1644
1645 else
1646
1647 !..Melt snow and graupel and enhance from collisions with liquid.
1648 !.. We also need to sublimate snow and graupel if subsaturated.
1649 if (L_qs(k)) then
1650 prr_sml(k) = tempc*tcond(k)*(t1_qs_me*smo1(k) &
1651 + t2_qs_me*rhof2(k)*vsc2(k)*smof(k))
1652 prr_sml(k) = prr_sml(k) + 4218.*olfus*tempc &
1653 * (prr_rcs(k)+prs_scw(k))
1654 prr_sml(k) = MIN(DBLE(rs(k)*odts), prr_sml(k))
1655
1656 if (ssati(k).lt. 0.) then
1657 prs_sde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
1658 * (t1_qs_sd*smo1(k) &
1659 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
1660 prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k))
1661 endif
1662 endif
1663
1664 if (L_qg(k)) then
1665 prr_gml(k) = tempc*N0_g(k)*tcond(k) &
1666 *(t1_qg_me*ilamg(k)**cge(10) &
1667 + t2_qg_me*rhof2(k)*vsc2(k)*ilamg(k)**cge(11))
1668 prr_gml(k) = prr_gml(k) + 4218.*olfus*tempc &
1669 * (prr_rcg(k)+prg_gcw(k))
1670 prr_gml(k) = MIN(DBLE(rg(k)*odts), prr_gml(k))
1671
1672 if (ssati(k).lt. 0.) then
1673 prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
1674 * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
1675 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
1676 prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k))
1677 endif
1678 endif
1679
1680 endif
1681
1682 enddo
1683 endif
1684
1685 !+---+-----------------------------------------------------------------+
1686 !..Ensure we do not deplete more hydrometeor species than exists.
1687 !+---+-----------------------------------------------------------------+
1688 do k = kts, kte
1689
1690 !..If ice supersaturated, ensure sum of depos growth terms does not
1691 !.. deplete more vapor than possibly exists. If subsaturated, limit
1692 !.. sum of sublimation terms such that vapor does not reproduce ice
1693 !.. supersat again.
1694 sump = pri_inu(k) + pri_ide(k) + prs_ide(k) &
1695 + prs_sde(k) + prg_gde(k)
1696 rate_max = (qv(k)-qvsi(k))*odts*0.999
1697 if ( (sump.gt. eps .and. sump.gt. rate_max) .or. &
1698 (sump.lt. -eps .and. sump.lt. rate_max) ) then
1699 ratio = rate_max/sump
1700 pri_inu(k) = pri_inu(k) * ratio
1701 pri_ide(k) = pri_ide(k) * ratio
1702 pni_ide(k) = pni_ide(k) * ratio
1703 prs_ide(k) = prs_ide(k) * ratio
1704 prs_sde(k) = prs_sde(k) * ratio
1705 prg_gde(k) = prg_gde(k) * ratio
1706 endif
1707
1708 !..Cloud water conservation.
1709 sump = -prr_wau(k) - pri_wfz(k) - prr_rcw(k) &
1710 - prs_scw(k) - prg_scw(k) - prg_gcw(k)
1711 rate_max = -rc(k)*odts
1712 if (sump.lt. rate_max .and. L_qc(k)) then
1713 ratio = rate_max/sump
1714 prr_wau(k) = prr_wau(k) * ratio
1715 pri_wfz(k) = pri_wfz(k) * ratio
1716 prr_rcw(k) = prr_rcw(k) * ratio
1717 prs_scw(k) = prs_scw(k) * ratio
1718 prg_scw(k) = prg_scw(k) * ratio
1719 prg_gcw(k) = prg_gcw(k) * ratio
1720 endif
1721
1722 !..Cloud ice conservation.
1723 sump = pri_ide(k) - prs_iau(k) - prs_sci(k) &
1724 - pri_rci(k)
1725 rate_max = -ri(k)*odts
1726 if (sump.lt. rate_max .and. L_qi(k)) then
1727 ratio = rate_max/sump
1728 pri_ide(k) = pri_ide(k) * ratio
1729 prs_iau(k) = prs_iau(k) * ratio
1730 prs_sci(k) = prs_sci(k) * ratio
1731 pri_rci(k) = pri_rci(k) * ratio
1732 endif
1733
1734 !..Rain conservation.
1735 sump = -prg_rfz(k) - pri_rfz(k) - prr_rci(k) &
1736 + prr_rcs(k) + prr_rcg(k)
1737 rate_max = -rr(k)*odts
1738 if (sump.lt. rate_max .and. L_qr(k)) then
1739 ratio = rate_max/sump
1740 prg_rfz(k) = prg_rfz(k) * ratio
1741 pri_rfz(k) = pri_rfz(k) * ratio
1742 prr_rci(k) = prr_rci(k) * ratio
1743 prr_rcs(k) = prr_rcs(k) * ratio
1744 prr_rcg(k) = prr_rcg(k) * ratio
1745 endif
1746
1747 !..Snow conservation.
1748 sump = prs_sde(k) - prs_ihm(k) - prr_sml(k) &
1749 + prs_rcs(k)
1750 rate_max = -rs(k)*odts
1751 if (sump.lt. rate_max .and. L_qs(k)) then
1752 ratio = rate_max/sump
1753 prs_sde(k) = prs_sde(k) * ratio
1754 prs_ihm(k) = prs_ihm(k) * ratio
1755 prr_sml(k) = prr_sml(k) * ratio
1756 prs_rcs(k) = prs_rcs(k) * ratio
1757 endif
1758
1759 !..Graupel conservation.
1760 sump = prg_gde(k) - prg_ihm(k) - prr_gml(k) &
1761 + prg_rcg(k)
1762 rate_max = -rg(k)*odts
1763 if (sump.lt. rate_max .and. L_qg(k)) then
1764 ratio = rate_max/sump
1765 prg_gde(k) = prg_gde(k) * ratio
1766 prg_ihm(k) = prg_ihm(k) * ratio
1767 prr_gml(k) = prr_gml(k) * ratio
1768 prg_rcg(k) = prg_rcg(k) * ratio
1769 endif
1770
1771 enddo
1772
1773 !+---+-----------------------------------------------------------------+
1774 !..Calculate tendencies of all species but constrain the number of ice
1775 !.. to reasonable values.
1776 !+---+-----------------------------------------------------------------+
1777 do k = kts, kte
1778 orho = 1./rho(k)
1779 lfus2 = lsub - lvap(k)
1780
1781 !..Water vapor tendency
1782 qvten(k) = qvten(k) + (-pri_inu(k) - pri_ide(k) &
1783 - prs_ide(k) - prs_sde(k) - prg_gde(k)) &
1784 * orho
1785
1786 !..Cloud water tendency
1787 qcten(k) = qcten(k) + (-prr_wau(k) - pri_wfz(k) &
1788 - prr_rcw(k) - prs_scw(k) - prg_scw(k) &
1789 - prg_gcw(k)) &
1790 * orho
1791
1792 !..Cloud ice mixing ratio tendency
1793 qiten(k) = qiten(k) + (pri_inu(k) + pri_ihm(k) &
1794 + pri_wfz(k) + pri_rfz(k) + pri_ide(k) &
1795 - prs_iau(k) - prs_sci(k) - pri_rci(k)) &
1796 * orho
1797
1798 !..Cloud ice number tendency.
1799 niten(k) = niten(k) + (pni_inu(k) + pni_ihm(k) &
1800 + pni_wfz(k) + pni_rfz(k) + pni_ide(k) &
1801 - pni_iau(k) - pni_sci(k) - pni_rci(k)) &
1802 * orho
1803
1804 !..Cloud ice mass/number balance; keep mass-wt mean size between 0.30*D0s
1805 !.. and 5*D0s. Also no more than 500 xtals per liter.
1806 xri=MAX(R1,(qi1d(k) + qiten(k)*dts)*rho(k))
1807 xni=MAX(1.,(ni1d(k) + niten(k)*dts)*rho(k))
1808 if (xri.gt. R1) then
1809 lami = (am_i*cig(2)*oig1*xni/xri)**obmi
1810 ilami = 1./lami
1811 Di = (bm_i + mu_i + 1.) * ilami
1812 if (Di.lt. 0.30*D0s) then
1813 lami = cie(2)/(0.30*D0s)
1814 xni = MIN(5.D5, cig(1)*oig2*xri/am_i*lami**bm_i)
1815 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
1816 elseif (Di.gt. 5.0*D0s) then
1817 lami = cie(2)/(5.0*D0s)
1818 xni = cig(1)*oig2*xri/am_i*lami**bm_i
1819 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
1820 endif
1821 else
1822 niten(k) = -ni1d(k)*odts
1823 endif
1824 xni=MAX(0.,(ni1d(k) + niten(k)*dts)*rho(k))
1825 if (xni.gt.5.E5) &
1826 niten(k) = (5.E5-ni1d(k)*rho(k))*odts*orho
1827
1828 !..Rain tendency
1829 qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) &
1830 + prr_sml(k) + prr_gml(k) + prr_rcs(k) &
1831 + prr_rcg(k) - prg_rfz(k) &
1832 - pri_rfz(k) - prr_rci(k)) &
1833 * orho
1834
1835 !..Snow tendency
1836 qsten(k) = qsten(k) + (prs_iau(k) + prs_sde(k) &
1837 + prs_sci(k) + prs_scw(k) + prs_rcs(k) &
1838 + prs_ide(k) - prs_ihm(k) - prr_sml(k)) &
1839 * orho
1840
1841 !..Graupel tendency
1842 qgten(k) = qgten(k) + (prg_scw(k) + prg_rfz(k) &
1843 + prg_gde(k) + prg_rcg(k) + prg_gcw(k) &
1844 + prg_rci(k) + prg_rcs(k) - prg_ihm(k) &
1845 - prr_gml(k)) &
1846 * orho
1847
1848 !..Temperature tendency
1849 if (temp(k).lt.T_0) then
1850 tten(k) = tten(k) &
1851 + ( lsub*ocp(k)*(pri_inu(k) + pri_ide(k) &
1852 + prs_ide(k) + prs_sde(k) &
1853 + prg_gde(k)) &
1854 + lfus2*ocp(k)*(pri_wfz(k) + pri_rfz(k) &
1855 + prg_rfz(k) + prs_scw(k) &
1856 + prg_scw(k) + prg_gcw(k) &
1857 + prg_rcs(k) + prs_rcs(k) &
1858 + prr_rci(k) + prg_rcg(k)) &
1859 )*orho * (1-IFDRY)
1860 else
1861 tten(k) = tten(k) &
1862 + ( lfus*ocp(k)*(-prr_sml(k) - prr_gml(k) &
1863 - prr_rcg(k) - prr_rcs(k)) &
1864 + lsub*ocp(k)*(prs_sde(k) + prg_gde(k)) &
1865 )*orho * (1-IFDRY)
1866 endif
1867
1868 enddo
1869
1870 !+---+-----------------------------------------------------------------+
1871 !..Update variables for TAU+1 before condensation & sedimention.
1872 !+---+-----------------------------------------------------------------+
1873 do k = kts, kte
1874 temp(k) = t1d(k) + DT*tten(k)
1875 otemp = 1./temp(k)
1876 tempc = temp(k) - 273.15
1877 qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
1878 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
1879 rhof(k) = SQRT(RHO_NOT/rho(k))
1880 rhof2(k) = SQRT(rhof(k))
1881 qvs(k) = rslf(pres(k), temp(k))
1882 ssatw(k) = qv(k)/qvs(k) - 1.
1883 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
1884 diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
1885 if (tempc .ge. 0.0) then
1886 visco(k) = (1.718+0.0049*tempc)*1.0E-5
1887 else
1888 visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
1889 endif
1890 vsc2(k) = SQRT(rho(k)/visco(k))
1891 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
1892 tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
1893 ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
1894 lvt2(k)=lvap(k)*lvap(k)*ocp(k)*oRv*otemp*otemp
1895
1896 if ((qc1d(k) + qcten(k)*DT) .gt. R1) then
1897 rc(k) = (qc1d(k) + qcten(k)*DT)*rho(k)
1898 L_qc(k) = .true.
1899 else
1900 rc(k) = R1
1901 L_qc(k) = .false.
1902 endif
1903
1904 if ((qi1d(k) + qiten(k)*DT) .gt. R1) then
1905 ri(k) = (qi1d(k) + qiten(k)*DT)*rho(k)
1906 ni(k) = MAX(1., (ni1d(k) + niten(k)*DT)*rho(k))
1907 L_qi(k) = .true.
1908 else
1909 ri(k) = R1
1910 ni(k) = 1.
1911 L_qi(k) = .false.
1912 endif
1913
1914 if ((qr1d(k) + qrten(k)*DT) .gt. R1) then
1915 rr(k) = (qr1d(k) + qrten(k)*DT)*rho(k)
1916 L_qr(k) = .true.
1917 else
1918 rr(k) = R1
1919 L_qr(k) = .false.
1920 endif
1921
1922 if ((qs1d(k) + qsten(k)*DT) .gt. R1) then
1923 rs(k) = (qs1d(k) + qsten(k)*DT)*rho(k)
1924 L_qs(k) = .true.
1925 else
1926 rs(k) = R1
1927 L_qs(k) = .false.
1928 endif
1929
1930 if ((qg1d(k) + qgten(k)*DT) .gt. R1) then
1931 rg(k) = (qg1d(k) + qgten(k)*DT)*rho(k)
1932 L_qg(k) = .true.
1933 else
1934 rg(k) = R1
1935 L_qg(k) = .false.
1936 endif
1937 enddo
1938
1939 !+---+-----------------------------------------------------------------+
1940 !..With tendency-updated mixing ratios, recalculate snow moments and
1941 !.. intercepts/slopes of graupel and rain.
1942 !+---+-----------------------------------------------------------------+
1943 if (.not. iiwarm) then
1944 do k = kts, kte
1945 if (.not. L_qs(k)) CYCLE
1946 tc0 = MIN(-0.1, temp(k)-273.15)
1947 smob(k) = rs(k)*oams
1948
1949 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
1950 !.. then we must compute actual 2nd moment and use as reference.
1951 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
1952 smo2(k) = smob(k)
1953 else
1954 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
1955 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
1956 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
1957 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
1958 + sa(10)*bm_s*bm_s*bm_s
1959 a_ = 10.0**loga_
1960 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
1961 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
1962 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
1963 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
1964 + sb(10)*bm_s*bm_s*bm_s
1965 smo2(k) = (smob(k)/a_)**(1./b_)
1966 endif
1967
1968 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
1969 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
1970 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
1971 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
1972 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
1973 + sa(10)*cse(1)*cse(1)*cse(1)
1974 a_ = 10.0**loga_
1975 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
1976 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
1977 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
1978 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
1979 smoc(k) = a_ * smo2(k)**b_
1980
1981 !..Calculate bm_s+bv_s (th) moment. Useful for sedimentation.
1982 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(14) &
1983 + sa(4)*tc0*cse(14) + sa(5)*tc0*tc0 &
1984 + sa(6)*cse(14)*cse(14) + sa(7)*tc0*tc0*cse(14) &
1985 + sa(8)*tc0*cse(14)*cse(14) + sa(9)*tc0*tc0*tc0 &
1986 + sa(10)*cse(14)*cse(14)*cse(14)
1987 a_ = 10.0**loga_
1988 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(14) + sb(4)*tc0*cse(14) &
1989 + sb(5)*tc0*tc0 + sb(6)*cse(14)*cse(14) &
1990 + sb(7)*tc0*tc0*cse(14) + sb(8)*tc0*cse(14)*cse(14) &
1991 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(14)*cse(14)*cse(14)
1992 smod(k) = a_ * smo2(k)**b_
1993 enddo
1994
1995 !+---+-----------------------------------------------------------------+
1996 !..Calculate y-intercept, slope values for graupel.
1997 !+---+-----------------------------------------------------------------+
1998 N0_min = gonv_max
1999 do k = kte, kts, -1
2000 if (.not. L_qg(k)) CYCLE
2001 N0_exp = 100.0*rho(k)/rg(k)
2002 N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max)))
2003 N0_min = MIN(N0_exp, N0_min)
2004 N0_exp = N0_min
2005 lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
2006 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
2007 ilamg(k) = 1./lamg
2008 N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2)
2009 enddo
2010
2011 endif
2012
2013 !+---+-----------------------------------------------------------------+
2014 !..Calculate y-intercept, slope values for rain.
2015 !+---+-----------------------------------------------------------------+
2016 N0_min = ronv_max
2017 do k = kte, kts, -1
2018 ! if (.not. L_qr(k)) CYCLE
2019 N0_exp = ronv_c1*tanh(ronv_c0*(ronv_r0-rr(k))) + ronv_c2
2020 N0_min = MIN(N0_exp, N0_min)
2021 N0_exp = N0_min
2022 lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
2023 lamr = lam_exp * (crg(3)*org2*org1)**obmr
2024 mvd_r(k) = (3.0+mu_r+0.672) / lamr
2025 if (mvd_r(k) .gt. 3.e-3) then
2026 mvd_r(k) = 3.e-3
2027 lamr = (3.0+mu_r+0.672) / 3.e-3
2028 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
2029 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
2030 endif
2031 N0_r(k) = N0_exp/(crg(2)*lam_exp) * lamr**cre(2)
2032 ilamr(k) = 1./lamr
2033 enddo
2034
2035 if (.not. iiwarm) then
2036 k_0 = kts
2037 melti = .false.
2038 do k = kte-1, kts, -1
2039 if ( (temp(k).gt. T_0) .and. (rr(k).gt. 0.001e-3) &
2040 .and. ((rs(k+1)+rg(k+1)).gt. 0.01e-3) ) then
2041 k_0 = MAX(k+1, k_0)
2042 melti=.true.
2043 goto 137
2044 endif
2045 enddo
2046 137 continue
2047
2048 if (melti) then
2049 !.. Locate bottom of melting layer (if any).
2050 kbot = kts
2051 do k = k_0-1, kts, -1
2052 if ( (rs(k)+rg(k)).lt. 0.01e-3) goto 138
2053 enddo
2054 138 continue
2055 kbot = MAX(k, kts)
2056
2057 !.. Compute melted snow/graupel equiv water diameter one K-level above
2058 !.. melting. Set starting rain mvd to either 50 microns or max from
2059 !.. higher up in column.
2060 if (L_qs(k_0)) then
2061 Ds = smoc(k_0) / smob(k_0)
2062 Ds_m = (am_s*Ds**bm_s / am_r)**obmr
2063 else
2064 Ds_m = 1.0e-6
2065 endif
2066 if (L_qg(k_0)) then
2067 Dg = (bm_g + mu_g + 1.) * ilamg(k_0)
2068 Dg_m = (am_g*Dg**bm_g / am_r)**obmr
2069 else
2070 Dg_m = 1.0e-6
2071 endif
2072 r_mvd1 = mvd_r(k_0)
2073 r_mvd2 = MIN(MAX(Ds_m, Dg_m, r_mvd1+1.e-6, mvd_r(kbot)), &
2074 3.e-3)
2075
2076 !.. Within melting layer, apply linear increase of rain mvd from r_mvd1
2077 !.. to equiv melted snow/graupel value (r_mvd2). So, by the bottom of
2078 !.. the melting layer, the rain will have an mvd that matches that from
2079 !.. melted snow and/or graupel.
2080 if (kbot.gt. 2) then
2081 do k = k_0-1, kbot, -1
2082 if (.not. L_qr(k)) CYCLE
2083 xkrat = REAL(k_0-k)/REAL(k_0-kbot)
2084 mvd_r(k) = MAX(mvd_r(k), xkrat*(r_mvd2-r_mvd1)+r_mvd1)
2085 lamr = (4.0+mu_r) / mvd_r(k)
2086 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
2087 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
2088 N0_exp = MAX(DBLE(ronv_min), MIN(N0_exp, DBLE(ronv_max)))
2089 lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
2090 lamr = lam_exp * (crg(3)*org2*org1)**obmr
2091 mvd_r(k) = (3.0+mu_r+0.672) / lamr
2092 N0_r(k) = rr(k)*lamr**cre(3) / (am_r*crg(3))
2093 ilamr(k) = 1./lamr
2094 enddo
2095
2096 !.. Below melting layer, hold N0_r constant unless changes to mixing
2097 !.. ratio increase mvd beyond 3 mm threshold, then adjust slope and
2098 !.. intercept to cap mvd at 3 mm. In future, we could lower N0_r to
2099 !.. account for self-collection or other sinks.
2100 do k = kbot-1, kts, -1
2101 if (.not. L_qr(k)) CYCLE
2102 N0_r(k) = MIN(N0_r(k), N0_r(kbot))
2103 lamr = (N0_r(k)*am_r*crg(3)/rr(k))**(1./cre(3))
2104 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
2105 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
2106 N0_exp = MAX(DBLE(ronv_min), MIN(N0_exp, DBLE(ronv_max)))
2107 lam_exp = (N0_exp*am_r*crg(1)/rr(k))**ore1
2108 lamr = lam_exp * (crg(3)*org2*org1)**obmr
2109 mvd_r(k) = (3.0+mu_r+0.672) / lamr
2110 if (mvd_r(k) .gt. 3.e-3) then
2111 mvd_r(k) = 3.e-3
2112 lamr = (3.0+mu_r+0.672) / mvd_r(k)
2113 endif
2114 N0_r(k) = rr(k)*lamr**cre(3) / (am_r*crg(3))
2115 ilamr(k) = 1./lamr
2116 enddo
2117 endif
2118
2119
2120 endif
2121 endif
2122
2123 !+---+-----------------------------------------------------------------+
2124 !..Cloud water condensation and evaporation. Newly formulated using
2125 !.. Newton-Raphson iterations (3 should suffice) as provided by B. Hall.
2126 !+---+-----------------------------------------------------------------+
2127 do k = kts, kte
2128 if ( (ssatw(k).gt. eps) .or. (ssatw(k).lt. -eps .and. &
2129 L_qc(k)) ) then
2130 clap = (qv(k)-qvs(k))/(1. + lvt2(k)*qvs(k))
2131 do n = 1, 3
2132 fcd = qvs(k)* EXP(lvt2(k)*clap) - qv(k) + clap
2133 dfcd = qvs(k)*lvt2(k)* EXP(lvt2(k)*clap) + 1.
2134 clap = clap - fcd/dfcd
2135 enddo
2136 xrc = rc(k) + clap
2137 if (xrc.gt. 0.0) then
2138 prw_vcd(k) = clap*odt
2139 else
2140 prw_vcd(k) = -rc(k)/rho(k)*odt
2141 endif
2142
2143 qcten(k) = qcten(k) + prw_vcd(k)
2144 qvten(k) = qvten(k) - prw_vcd(k)
2145 tten(k) = tten(k) + lvap(k)*ocp(k)*prw_vcd(k)*(1-IFDRY)
2146 rc(k) = MAX(R1, (qc1d(k) + DT*qcten(k))*rho(k))
2147 qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
2148 temp(k) = t1d(k) + DT*tten(k)
2149 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
2150 qvs(k) = rslf(pres(k), temp(k))
2151 ssatw(k) = qv(k)/qvs(k) - 1.
2152 endif
2153 enddo
2154
2155 !+---+-----------------------------------------------------------------+
2156 !.. If still subsaturated, allow rain to evaporate, following
2157 !.. Srivastava & Coen (1992).
2158 !+---+-----------------------------------------------------------------+
2159 do k = kts, kte
2160 if ( (ssatw(k).lt. -eps) .and. L_qr(k) &
2161 .and. (.not.(prw_vcd(k).gt. 0.)) ) then
2162 tempc = temp(k) - 273.15
2163 otemp = 1./temp(k)
2164 rhof(k) = SQRT(RHO_NOT/rho(k))
2165 rhof2(k) = SQRT(rhof(k))
2166 diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
2167 if (tempc .ge. 0.0) then
2168 visco(k) = (1.718+0.0049*tempc)*1.0E-5
2169 else
2170 visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
2171 endif
2172 vsc2(k) = SQRT(rho(k)/visco(k))
2173 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
2174 tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
2175 ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
2176
2177 rvs = rho(k)*qvs(k)
2178 rvs_p = rvs*otemp*(lvap(k)*otemp*oRv - 1.)
2179 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*oRv - 1.) &
2180 *otemp*(lvap(k)*otemp*oRv - 1.) &
2181 + (-2.*lvap(k)*otemp*otemp*otemp*oRv) &
2182 + otemp*otemp)
2183 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
2184 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
2185 * rvs_pp/rvs_p * rvs/rvs_p
2186 alphsc = MAX(1.E-9, alphsc)
2187 xsat = ssatw(k)
2188 if (xsat.lt. -1.E-9) xsat = -1.E-9
2189 t1_evap = 2.*PI*( 1.0 - alphsc*xsat &
2190 + 2.*alphsc*alphsc*xsat*xsat &
2191 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
2192 / (1.+gamsc)
2193
2194 lamr = 1./ilamr(k)
2195 prv_rev(k) = t1_evap*diffu(k)*(-ssatw(k))*N0_r(k)*rvs &
2196 * (t1_qr_ev*ilamr(k)**cre(10) &
2197 + t2_qr_ev*vsc2(k)*rhof2(k)*((lamr+0.5*fv_r)**(-cre(11))))
2198 prv_rev(k) = MIN(DBLE(rr(k)/rho(k)*odts), &
2199 prv_rev(k)/rho(k))
2200
2201 qrten(k) = qrten(k) - prv_rev(k)
2202 qvten(k) = qvten(k) + prv_rev(k)
2203 tten(k) = tten(k) - lvap(k)*ocp(k)*prv_rev(k)*(1-IFDRY)
2204
2205 rr(k) = MAX(R1, (qr1d(k) + DT*qrten(k))*rho(k))
2206 qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
2207 temp(k) = t1d(k) + DT*tten(k)
2208 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
2209 endif
2210 enddo
2211
2212 !+---+-----------------------------------------------------------------+
2213 !..Find max terminal fallspeed (distribution mass-weighted mean
2214 !.. velocity) and use it to determine if we need to split the timestep
2215 !.. (var nstep>1). Either way, only bother to do sedimentation below
2216 !.. 1st level that contains any sedimenting particles (k=ksed1 on down).
2217 !+---+-----------------------------------------------------------------+
2218 nstep = 0
2219 ksed1 = 0
2220 do k = kte+1, kts, -1
2221 vtrk(k) = 0.
2222 vtik(k) = 0.
2223 vtnk(k) = 0.
2224 vtsk(k) = 0.
2225 vtgk(k) = 0.
2226 enddo
2227 do k = kte, kts, -1
2228 vtr = 0.
2229 vti = 0.
2230 vts = 0.
2231 vtg = 0.
2232 rhof(k) = SQRT(RHO_NOT/rho(k))
2233
2234 if (rr(k).gt. R2) then
2235 lamr = 1./ilamr(k)
2236 vtr = rhof(k)*av_r*crg(6)*org3 * (lamr/(lamr+fv_r))**cre(3) &
2237 *((lamr+fv_r)**(-bv_r))
2238 vtrk(k) = vtr
2239 endif
2240
2241 if (.not. iiwarm) then
2242 if (ri(k).gt. R2) then
2243 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
2244 ilami = 1./lami
2245 vti = rhof(k)*av_i*cig(3)*oig2 * ilami**bv_i
2246 vtik(k) = vti
2247 vti = rhof(k)*av_i*cig(4)*oig1 * ilami**bv_i
2248 vtnk(k) = vti
2249 endif
2250
2251 if (rs(k).gt. R2) then
2252 Ds = smoc(k) / smob(k)
2253 Mrat = 1./Ds
2254 ils1 = 1./(Mrat*Lam0 + fv_s)
2255 ils2 = 1./(Mrat*Lam1 + fv_s)
2256 t1_vts = Kap0*csg(4)*ils1**cse(4)
2257 t2_vts = Kap1*Mrat**mu_s*csg(10)*ils2**cse(10)
2258 t3_vts = Kap0*csg(1)*ils1**cse(1)
2259 t4_vts = Kap1*Mrat**mu_s*csg(7)*ils2**cse(7)
2260 vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts)
2261 if (temp(k).gt. T_0) then
2262 vtsk(k) = MAX(vts*vts_boost(k), vtrk(k))
2263 else
2264 vtsk(k) = vts*vts_boost(k)
2265 endif
2266 endif
2267
2268 if (rg(k).gt. R2) then
2269 vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
2270 if (temp(k).gt. T_0) then
2271 vtgk(k) = MAX(vtg, vtrk(k))
2272 else
2273 vtgk(k) = vtg
2274 endif
2275 endif
2276 endif
2277
2278 rgvm = MAX(vtik(k), vtrk(k), vtsk(k), vtgk(k))
2279 if (rgvm .gt. 1.E-3) then
2280 ksed1 = MAX(ksed1, k)
2281 delta_tp = dzq(k)/rgvm
2282 nstep = MAX(nstep, INT(DT/delta_tp + 1.))
2283 endif
2284 enddo
2285 if (ksed1 .eq. kte) ksed1 = kte-1
2286 if (nstep .gt. 0) onstep = 1./REAL(nstep)
2287
2288 !+---+-----------------------------------------------------------------+
2289 !..Sedimentation of mixing ratio is the integral of v(D)*m(D)*N(D)*dD,
2290 !.. whereas neglect m(D) term for number concentration. Therefore,
2291 !.. cloud ice has proper differential sedimentation.
2292 !+---+-----------------------------------------------------------------+
2293 do n = 1, nstep
2294 do k = kte, kts, -1
2295 sed_r(k) = vtrk(k)*rr(k)
2296 sed_i(k) = vtik(k)*ri(k)
2297 sed_n(k) = vtnk(k)*ni(k)
2298 sed_g(k) = vtgk(k)*rg(k)
2299 sed_s(k) = vtsk(k)*rs(k)
2300 enddo
2301 k = kte
2302 odzq = 1./dzq(k)
2303 orho = 1./rho(k)
2304 qrten(k) = qrten(k) - sed_r(k)*odzq*onstep*orho
2305 qiten(k) = qiten(k) - sed_i(k)*odzq*onstep*orho
2306 niten(k) = niten(k) - sed_n(k)*odzq*onstep*orho
2307 qgten(k) = qgten(k) - sed_g(k)*odzq*onstep*orho
2308 qsten(k) = qsten(k) - sed_s(k)*odzq*onstep*orho
2309 rr(k) = MAX(R1, rr(k) - sed_r(k)*odzq*DT*onstep)
2310 ri(k) = MAX(R1, ri(k) - sed_i(k)*odzq*DT*onstep)
2311 ni(k) = MAX(1., ni(k) - sed_n(k)*odzq*DT*onstep)
2312 rg(k) = MAX(R1, rg(k) - sed_g(k)*odzq*DT*onstep)
2313 rs(k) = MAX(R1, rs(k) - sed_s(k)*odzq*DT*onstep)
2314 do k = ksed1, kts, -1
2315 odzq = 1./dzq(k)
2316 orho = 1./rho(k)
2317 qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) &
2318 *odzq*onstep*orho
2319 qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) &
2320 *odzq*onstep*orho
2321 niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) &
2322 *odzq*onstep*orho
2323 qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) &
2324 *odzq*onstep*orho
2325 qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) &
2326 *odzq*onstep*orho
2327 rr(k) = MAX(R1, rr(k) + (sed_r(k+1)-sed_r(k)) &
2328 *odzq*DT*onstep)
2329 ri(k) = MAX(R1, ri(k) + (sed_i(k+1)-sed_i(k)) &
2330 *odzq*DT*onstep)
2331 ni(k) = MAX(1., ni(k) + (sed_n(k+1)-sed_n(k)) &
2332 *odzq*DT*onstep)
2333 rg(k) = MAX(R1, rg(k) + (sed_g(k+1)-sed_g(k)) &
2334 *odzq*DT*onstep)
2335 rs(k) = MAX(R1, rs(k) + (sed_s(k+1)-sed_s(k)) &
2336 *odzq*DT*onstep)
2337 enddo
2338
2339 !+---+-----------------------------------------------------------------+
2340 !..Precipitation reaching the ground.
2341 !+---+-----------------------------------------------------------------+
2342 pptrain = pptrain + sed_r(kts)*DT*onstep
2343 pptsnow = pptsnow + sed_s(kts)*DT*onstep
2344 pptgraul = pptgraul + sed_g(kts)*DT*onstep
2345 pptice = pptice + sed_i(kts)*DT*onstep
2346
2347 enddo
2348
2349 !+---+-----------------------------------------------------------------+
2350 !.. Instantly melt any cloud ice into cloud water if above 0C and
2351 !.. instantly freeze any cloud water found below HGFR.
2352 !+---+-----------------------------------------------------------------+
2353 if (.not. iiwarm) then
2354 do k = kts, kte
2355 xri = MAX(0.0, qi1d(k) + qiten(k)*DT)
2356 if ( (temp(k).gt. T_0) .and. (xri.gt. 0.0) ) then
2357 qcten(k) = qcten(k) + xri*odt
2358 qiten(k) = -qi1d(k)*odt
2359 niten(k) = -ni1d(k)*odt
2360 tten(k) = tten(k) - lfus*ocp(k)*xri*odt*(1-IFDRY)
2361 endif
2362
2363 xrc = MAX(0.0, qc1d(k) + qcten(k)*DT)
2364 if ( (temp(k).lt. HGFR) .and. (xrc.gt. 0.0) ) then
2365 lfus2 = lsub - lvap(k)
2366 qiten(k) = qiten(k) + xrc*odt
2367 niten(k) = niten(k) + xrc/(2.*xm0i)*odt
2368 qcten(k) = -xrc*odt
2369 tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-IFDRY)
2370 endif
2371 enddo
2372 endif
2373
2374 !+---+-----------------------------------------------------------------+
2375 !.. All tendencies computed, apply and pass back final values to parent.
2376 !+---+-----------------------------------------------------------------+
2377 do k = kts, kte
2378 t1d(k) = t1d(k) + tten(k)*DT
2379 qv1d(k) = MAX(1.E-10, qv1d(k) + qvten(k)*DT)
2380 qc1d(k) = qc1d(k) + qcten(k)*DT
2381 if (qc1d(k) .le. R1) qc1d(k) = 0.0
2382 qi1d(k) = qi1d(k) + qiten(k)*DT
2383 ni1d(k) = ni1d(k) + niten(k)*DT
2384 if (qi1d(k) .le. R1) then
2385 qi1d(k) = 0.0
2386 ni1d(k) = 0.0
2387 else
2388 if (ni1d(k) .gt. 1.0) then
2389 lami = (am_i*cig(2)*oig1*ni1d(k)/qi1d(k))**obmi
2390 ilami = 1./lami
2391 Di = (bm_i + mu_i + 1.) * ilami
2392 if (Di.lt. 0.30*D0s) then
2393 lami = cie(2)/(0.30*D0s)
2394 ni1d(k) = MIN(5.D5, cig(1)*oig2*qi1d(k)/am_i*lami**bm_i)
2395 elseif (Di.gt. 5.0*D0s) then
2396 lami = cie(2)/(5.0*D0s)
2397 ni1d(k) = cig(1)*oig2*qi1d(k)/am_i*lami**bm_i
2398 endif
2399 else
2400 lami = cie(2)/(0.30*D0s)
2401 ni1d(k) = MIN(5.D5, cig(1)*oig2*qi1d(k)/am_i*lami**bm_i)
2402 endif
2403 endif
2404 qr1d(k) = qr1d(k) + qrten(k)*DT
2405 if (qr1d(k) .le. R1) qr1d(k) = 0.0
2406 qs1d(k) = qs1d(k) + qsten(k)*DT
2407 if (qs1d(k) .le. R1) qs1d(k) = 0.0
2408 qg1d(k) = qg1d(k) + qgten(k)*DT
2409 if (qg1d(k) .le. R1) qg1d(k) = 0.0
2410 enddo
2411
2412 end subroutine mp_thompson
2413 !+---+-----------------------------------------------------------------+
2414 !
2415 !+---+-----------------------------------------------------------------+
2416 !..Creation of the lookup tables and support functions found below here.
2417 !+---+-----------------------------------------------------------------+
2418 !..Rain collecting graupel (and inverse). Explicit CE integration.
2419 !+---+-----------------------------------------------------------------+
2420
2421 subroutine qr_acr_qg
2422
2423 implicit none
2424
2425 !..Local variables
2426 INTEGER:: i, j, k, n, n2
2427 INTEGER, PARAMETER:: nbr = 100
2428 INTEGER, PARAMETER:: nbg = 100
2429 DOUBLE PRECISION, DIMENSION(nbg):: Dg, vg, N_g, dtg
2430 DOUBLE PRECISION, DIMENSION(nbr):: Dr, vr, N_r, dtr
2431 DOUBLE PRECISION, DIMENSION(nbg+1):: Dx
2432 DOUBLE PRECISION, DIMENSION(nbr+1):: Dy
2433 DOUBLE PRECISION:: N0_exp, N0_r, N0_g, lam_exp, lamg, lamr, N0_s
2434 DOUBLE PRECISION:: massg, massr, dvg, dvr, t1, t2, z1, z2
2435
2436 !+---+
2437
2438 !..Create bins of rain (from min diameter up to 5 mm).
2439 Dy(1) = D0r*1.0d0
2440 Dy(nbr+1) = 0.005d0
2441 do n2 = 2, nbr
2442 Dy(n2) = DEXP(DFLOAT(n2-1)/DFLOAT(nbr) &
2443 *DLOG(Dy(nbr+1)/Dy(1)) +DLOG(Dy(1)))
2444 enddo
2445 do n2 = 1, nbr
2446 Dr(n2) = DSQRT(Dy(n2)*Dy(n2+1))
2447 vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
2448 dtr(n2) = Dy(n2+1) - Dy(n2)
2449 enddo
2450
2451 !..Create bins of graupel (from min diameter up to 5 cm).
2452 Dx(1) = D0g*1.0d0
2453 Dx(nbg+1) = 0.05d0
2454 do n = 2, nbg
2455 Dx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbg) &
2456 *DLOG(Dx(nbg+1)/Dx(1)) +DLOG(Dx(1)))
2457 enddo
2458 do n = 1, nbg
2459 Dg(n) = DSQRT(Dx(n)*Dx(n+1))
2460 vg(n) = av_g*Dg(n)**bv_g
2461 dtg(n) = Dx(n+1) - Dx(n)
2462 enddo
2463
2464 do k = 1, ntb_r
2465 do j = 1, ntb_r1
2466
2467 lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(k))**ore1
2468 lamr = lam_exp * (crg(3)*org2*org1)**obmr
2469 N0_r = N0r_exp(j)/(crg(2)*lam_exp) * lamr**cre(2)
2470 do n2 = 1, nbr
2471 N_r(n2) = N0_r*Dr(n2)**mu_r *DEXP(-lamr*Dr(n2))*dtr(n2)
2472 enddo
2473
2474 do i = 1, ntb_g
2475 N0_exp = 100.0d0/r_g(i)
2476 N0_exp = DMAX1(gonv_min*1.d0,DMIN1(N0_exp,gonv_max*1.d0))
2477 lam_exp = (N0_exp*am_g*cgg(1)/r_g(i))**oge1
2478 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
2479 N0_g = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2)
2480 do n = 1, nbg
2481 N_g(n) = N0_g*Dg(n)**mu_g * DEXP(-lamg*Dg(n))*dtg(n)
2482 enddo
2483
2484 t1 = 0.0d0
2485 t2 = 0.0d0
2486 z1 = 0.0d0
2487 z2 = 0.0d0
2488 do n2 = 1, nbr
2489 massr = am_r * Dr(n2)**bm_r
2490 do n = 1, nbg
2491 massg = am_g * Dg(n)**bm_g
2492
2493 dvg = 0.5d0*((vr(n2) - vg(n)) + DABS(vr(n2)-vg(n)))
2494 dvr = 0.5d0*((vg(n) - vr(n2)) + DABS(vg(n)-vr(n2)))
2495
2496 t1 = t1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
2497 *dvg*massg * N_g(n)* N_r(n2)
2498 z1 = z1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
2499 *dvg*massr * N_g(n)* N_r(n2)
2500
2501 t2 = t2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
2502 *dvr*massr * N_g(n)* N_r(n2)
2503 z2 = z2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
2504 *dvr*massg * N_g(n)* N_r(n2)
2505 enddo
2506 97 continue
2507 enddo
2508 tcg_racg(i,j,k) = t1
2509 tmr_racg(i,j,k) = DMIN1(z1, r_r(k)*1.0d0)
2510 tcr_gacr(i,j,k) = t2
2511 tmg_gacr(i,j,k) = z2
2512 enddo
2513 enddo
2514 enddo
2515
2516 end subroutine qr_acr_qg
2517 !+---+-----------------------------------------------------------------+
2518 !
2519 !+---+-----------------------------------------------------------------+
2520 !..Rain collecting snow (and inverse). Explicit CE integration.
2521 !+---+-----------------------------------------------------------------+
2522
2523 subroutine qr_acr_qs
2524
2525 implicit none
2526
2527 !..Local variables
2528 INTEGER:: i, j, k, m, n, n2
2529 INTEGER, PARAMETER:: nbr = 100
2530 INTEGER, PARAMETER:: nbs = 100
2531 DOUBLE PRECISION, DIMENSION(nbr):: Dr, vr, D1, N_r, dtr
2532 DOUBLE PRECISION, DIMENSION(nbs):: Ds, vs, N_s, dts
2533 DOUBLE PRECISION, DIMENSION(nbr+1):: Dy
2534 DOUBLE PRECISION, DIMENSION(nbs+1):: Dx
2535 DOUBLE PRECISION:: loga_, a_, b_, second, M0, M2, M3, Mrat, oM3
2536 DOUBLE PRECISION:: N0_r, lam_exp, lamr, slam1, slam2
2537 DOUBLE PRECISION:: dvs, dvr, masss, massr
2538 DOUBLE PRECISION:: t1, t2, t3, t4, z1, z2, z3, z4
2539
2540 !+---+
2541
2542 !..Create bins of rain (from min diameter up to 5 mm).
2543 Dy(1) = D0r*1.0d0
2544 Dy(nbr+1) = 0.005d0
2545 do n2 = 2, nbr
2546 Dy(n2) = DEXP(DFLOAT(n2-1)/DFLOAT(nbr) &
2547 *DLOG(Dy(nbr+1)/Dy(1)) +DLOG(Dy(1)))
2548 enddo
2549 do n2 = 1, nbr
2550 Dr(n2) = DSQRT(Dy(n2)*Dy(n2+1))
2551 vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
2552 D1(n2) = (vr(n2)/av_s)**(1./bv_s)
2553 dtr(n2) = Dy(n2+1) - Dy(n2)
2554 enddo
2555
2556 !..Create bins of snow (from min diameter up to 2 cm).
2557 Dx(1) = D0s*1.0d0
2558 Dx(nbs+1) = 0.02d0
2559 do n = 2, nbs
2560 Dx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbs) &
2561 *DLOG(Dx(nbs+1)/Dx(1)) +DLOG(Dx(1)))
2562 enddo
2563 do n = 1, nbs
2564 Ds(n) = DSQRT(Dx(n)*Dx(n+1))
2565 vs(n) = 1.5*av_s*Ds(n)**bv_s * DEXP(-fv_s*Ds(n))
2566 dts(n) = Dx(n+1) - Dx(n)
2567 enddo
2568
2569 do m = 1, ntb_r
2570 do k = 1, ntb_r1
2571 lam_exp = (N0r_exp(k)*am_r*crg(1)/r_r(m))**ore1
2572 lamr = lam_exp * (crg(3)*org2*org1)**obmr
2573 N0_r = N0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2)
2574 do n2 = 1, nbr
2575 N_r(n2) = N0_r*Dr(n2)**mu_r * DEXP(-lamr*Dr(n2))*dtr(n2)
2576 enddo
2577
2578 do j = 1, ntb_t
2579 do i = 1, ntb_s
2580
2581 !..From the bm_s moment, compute plus one moment. If we are not
2582 !.. using bm_s=2, then we must transform to the pure 2nd moment
2583 !.. (variable called "second") and then to the bm_s+1 moment.
2584
2585 M2 = r_s(i)*oams *1.0d0
2586 if (bm_s.gt.2.0-1.E-3 .and. bm_s.lt.2.0+1.E-3) then
2587 loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*bm_s &
2588 + sa(4)*Tc(j)*bm_s + sa(5)*Tc(j)*Tc(j) &
2589 + sa(6)*bm_s*bm_s + sa(7)*Tc(j)*Tc(j)*bm_s &
2590 + sa(8)*Tc(j)*bm_s*bm_s + sa(9)*Tc(j)*Tc(j)*Tc(j) &
2591 + sa(10)*bm_s*bm_s*bm_s
2592 a_ = 10.0**loga_
2593 b_ = sb(1) + sb(2)*Tc(j) + sb(3)*bm_s &
2594 + sb(4)*Tc(j)*bm_s + sb(5)*Tc(j)*Tc(j) &
2595 + sb(6)*bm_s*bm_s + sb(7)*Tc(j)*Tc(j)*bm_s &
2596 + sb(8)*Tc(j)*bm_s*bm_s + sb(9)*Tc(j)*Tc(j)*Tc(j) &
2597 + sb(10)*bm_s*bm_s*bm_s
2598 second = (M2/a_)**(1./b_)
2599 else
2600 second = M2
2601 endif
2602
2603 loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*cse(1) &
2604 + sa(4)*Tc(j)*cse(1) + sa(5)*Tc(j)*Tc(j) &
2605 + sa(6)*cse(1)*cse(1) + sa(7)*Tc(j)*Tc(j)*cse(1) &
2606 + sa(8)*Tc(j)*cse(1)*cse(1) + sa(9)*Tc(j)*Tc(j)*Tc(j) &
2607 + sa(10)*cse(1)*cse(1)*cse(1)
2608 a_ = 10.0**loga_
2609 b_ = sb(1)+sb(2)*Tc(j)+sb(3)*cse(1) + sb(4)*Tc(j)*cse(1) &
2610 + sb(5)*Tc(j)*Tc(j) + sb(6)*cse(1)*cse(1) &
2611 + sb(7)*Tc(j)*Tc(j)*cse(1) + sb(8)*Tc(j)*cse(1)*cse(1) &
2612 + sb(9)*Tc(j)*Tc(j)*Tc(j)+sb(10)*cse(1)*cse(1)*cse(1)
2613 M3 = a_ * second**b_
2614
2615 oM3 = 1./M3
2616 Mrat = M2*(M2*oM3)*(M2*oM3)*(M2*oM3)
2617 M0 = (M2*oM3)**mu_s
2618 slam1 = M2 * oM3 * Lam0
2619 slam2 = M2 * oM3 * Lam1
2620
2621 do n = 1, nbs
2622 N_s(n) = Mrat*(Kap0*DEXP(-slam1*Ds(n)) &
2623 + Kap1*M0*Ds(n)**mu_s * DEXP(-slam2*Ds(n)))*dts(n)
2624 enddo
2625
2626 t1 = 0.0d0
2627 t2 = 0.0d0
2628 t3 = 0.0d0
2629 t4 = 0.0d0
2630 z1 = 0.0d0
2631 z2 = 0.0d0
2632 z3 = 0.0d0
2633 z4 = 0.0d0
2634 do n2 = 1, nbr
2635 massr = am_r * Dr(n2)**bm_r
2636 do n = 1, nbs
2637 masss = am_s * Ds(n)**bm_s
2638
2639 dvs = 0.5d0*((vr(n2) - vs(n)) + DABS(vr(n2)-vs(n)))
2640 dvr = 0.5d0*((vs(n) - vr(n2)) + DABS(vs(n)-vr(n2)))
2641
2642 if (massr .gt. masss) then
2643 t1 = t1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2644 *dvs*masss * N_s(n)* N_r(n2)
2645 z1 = z1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2646 *dvs*massr * N_s(n)* N_r(n2)
2647 else
2648 t3 = t3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2649 *dvs*masss * N_s(n)* N_r(n2)
2650 z3 = z3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2651 *dvs*massr * N_s(n)* N_r(n2)
2652 endif
2653
2654 if (massr .gt. masss) then
2655 t2 = t2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2656 *dvr*massr * N_s(n)* N_r(n2)
2657 z2 = z2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2658 *dvr*masss * N_s(n)* N_r(n2)
2659 else
2660 t4 = t4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2661 *dvr*massr * N_s(n)* N_r(n2)
2662 z4 = z4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
2663 *dvr*masss * N_s(n)* N_r(n2)
2664 endif
2665
2666 enddo
2667 enddo
2668 tcs_racs1(i,j,k,m) = t1
2669 tmr_racs1(i,j,k,m) = DMIN1(z1, r_r(m)*1.0d0)
2670 tcs_racs2(i,j,k,m) = t3
2671 tmr_racs2(i,j,k,m) = z3
2672 tcr_sacr1(i,j,k,m) = t2
2673 tms_sacr1(i,j,k,m) = z2
2674 tcr_sacr2(i,j,k,m) = t4
2675 tms_sacr2(i,j,k,m) = z4
2676 enddo
2677 enddo
2678 enddo
2679 enddo
2680
2681 end subroutine qr_acr_qs
2682 !+---+-----------------------------------------------------------------+
2683 !
2684 !+---+-----------------------------------------------------------------+
2685 !..This is a literal adaptation of Bigg (1954) probability of drops of
2686 !..a particular volume freezing. Given this probability, simply freeze
2687 !..the proportion of drops summing their masses.
2688 !+---+-----------------------------------------------------------------+
2689
2690 subroutine freezeH2O
2691
2692 implicit none
2693
2694 !..Local variables
2695 INTEGER:: i, j, k, n, n2
2696 INTEGER, PARAMETER:: nbr = 100
2697 INTEGER, PARAMETER:: nbc = 50
2698 DOUBLE PRECISION, DIMENSION(nbr):: Dr, N_r, dtr, massr
2699 DOUBLE PRECISION, DIMENSION(nbr+1):: Dy
2700 DOUBLE PRECISION, DIMENSION(nbc):: Dc, N_c, dtc, massc
2701 DOUBLE PRECISION:: sum1, sum2, sumn1, sumn2, &
2702 prob, vol, Texp, orho_w, &
2703 lam_exp, lamr, N0_r, lamc, N0_c, y
2704
2705 !+---+
2706
2707 orho_w = 1./rho_w
2708
2709 !..Create bins of rain (from min diameter up to 5 mm).
2710 Dy(1) = D0r*1.0d0
2711 Dy(nbr+1) = 0.005d0
2712 do n2 = 2, nbr
2713 Dy(n2) = DEXP(DFLOAT(n2-1)/DFLOAT(nbr) &
2714 *DLOG(Dy(nbr+1)/Dy(1)) +DLOG(Dy(1)))
2715 enddo
2716 do n2 = 1, nbr
2717 Dr(n2) = DSQRT(Dy(n2)*Dy(n2+1))
2718 massr(n2) = am_r*Dr(n2)**bm_r
2719 dtr(n2) = Dy(n2+1) - Dy(n2)
2720 enddo
2721
2722 !..Create bins of cloud water (from min diameter up to 50 microns).
2723 Dc(1) = D0c*1.0d0
2724 massc(1) = am_r*Dc(1)**bm_r
2725 dtc(1) = D0c*1.0D6
2726 do n = 2, nbc
2727 Dc(n) = Dc(n-1) + 1.0D-6
2728 massc(n) = am_r*Dc(n)**bm_r
2729 dtc(n) = (Dc(n) - Dc(n-1)) * 1.0D6
2730 enddo
2731
2732 !..Freeze water (smallest drops become cloud ice, otherwise graupel).
2733 do k = 1, 45
2734 ! print*, ' Freezing water for temp = ', -k
2735 Texp = DEXP( DFLOAT(k) ) - 1.0D0
2736 do j = 1, ntb_r1
2737 do i = 1, ntb_r
2738 lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(i))**ore1
2739 lamr = lam_exp * (crg(3)*org2*org1)**obmr
2740 N0_r = N0r_exp(j)/(crg(2)*lam_exp) * lamr**cre(2)
2741 sum1 = 0.0d0
2742 sum2 = 0.0d0
2743 sumn1 = 0.0d0
2744 do n2 = 1, nbr
2745 N_r(n2) = N0_r*Dr(n2)**mu_r*DEXP(-lamr*Dr(n2))*dtr(n2)
2746 vol = massr(n2)*orho_w
2747 prob = 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp)
2748 if (massr(n2) .lt. xm0g) then
2749 sumn1 = sumn1 + prob*N_r(n2)
2750 sum1 = sum1 + prob*N_r(n2)*massr(n2)
2751 else
2752 sum2 = sum2 + prob*N_r(n2)*massr(n2)
2753 endif
2754 enddo
2755 tpi_qrfz(i,j,k) = sum1
2756 tni_qrfz(i,j,k) = sumn1
2757 tpg_qrfz(i,j,k) = sum2
2758 enddo
2759 enddo
2760 do i = 1, ntb_c
2761 lamc = 1.0D-6 * (Nt_c*am_r* ccg(2) * ocg1 / r_c(i))**obmr
2762 N0_c = 1.0D-18 * Nt_c*ocg1 * lamc**cce(1)
2763 sum1 = 0.0d0
2764 sumn2 = 0.0d0
2765 do n = 1, nbc
2766 y = Dc(n)*1.0D6
2767 vol = massc(n)*orho_w
2768 prob = 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp)
2769 N_c(n) = N0_c* y**mu_c * EXP(-lamc*y)*dtc(n)
2770 N_c(n) = 1.0D18 * N_c(n)
2771 sumn2 = sumn2 + prob*N_c(n)
2772 sum1 = sum1 + prob*N_c(n)*massc(n)
2773 enddo
2774 tpi_qcfz(i,k) = sum1
2775 tni_qcfz(i,k) = sumn2
2776 enddo
2777 enddo
2778
2779 end subroutine freezeH2O
2780 !+---+-----------------------------------------------------------------+
2781 !
2782 !+---+-----------------------------------------------------------------+
2783 !..Cloud ice converting to snow since portion greater than min snow
2784 !.. size. Given cloud ice content (kg/m**3), number concentration
2785 !.. (#/m**3) and gamma shape parameter, mu_i, break the distrib into
2786 !.. bins and figure out the mass/number of ice with sizes larger than
2787 !.. D0s. Also, compute incomplete gamma function for the integration
2788 !.. of ice depositional growth from diameter=0 to D0s. Amount of
2789 !.. ice depositional growth is this portion of distrib while larger
2790 !.. diameters contribute to snow growth (as in Harrington et al. 1995).
2791 !+---+-----------------------------------------------------------------+
2792
2793 subroutine qi_aut_qs
2794
2795 implicit none
2796
2797 !..Local variables
2798 INTEGER:: i, j, n2
2799 INTEGER, PARAMETER:: nbi = 100
2800 DOUBLE PRECISION, DIMENSION(nbi):: Di, N_i, dti
2801 DOUBLE PRECISION, DIMENSION(nbi+1):: Dy
2802 DOUBLE PRECISION:: N0_i, lami, Di_mean, t1, t2
2803
2804 !+---+
2805
2806 !..Create bins of cloud ice (from min diameter up to 10x min snow size).
2807 Dy(1) = D0i*1.0d0
2808 Dy(nbi+1) = 10.0d0*D0s
2809 do n2 = 2, nbi
2810 Dy(n2) = DEXP(DFLOAT(n2-1)/DFLOAT(nbi) &
2811 *DLOG(Dy(nbi+1)/Dy(1)) +DLOG(Dy(1)))
2812 enddo
2813 do n2 = 1, nbi
2814 Di(n2) = DSQRT(Dy(n2)*Dy(n2+1))
2815 dti(n2) = Dy(n2+1) - Dy(n2)
2816 enddo
2817
2818 do j = 1, ntb_i1
2819 do i = 1, ntb_i
2820 lami = (am_i*cig(2)*oig1*Nt_i(j)/r_i(i))**obmi
2821 Di_mean = (bm_i + mu_i + 1.) / lami
2822 N0_i = Nt_i(j)*oig1 * lami**cie(1)
2823 t1 = 0.0d0
2824 t2 = 0.0d0
2825 if (SNGL(Di_mean) .gt. 5.*D0s) then
2826 t1 = r_i(i)
2827 t2 = Nt_i(j)
2828 tpi_ide(i,j) = 0.0D0
2829 elseif (SNGL(Di_mean) .lt. D0i) then
2830 t1 = 0.0D0
2831 t2 = 0.0D0
2832 tpi_ide(i,j) = 1.0D0
2833 else
2834 tpi_ide(i,j) = GAMMP(mu_i+2.0, SNGL(lami)*D0s) * 1.0D0
2835 do n2 = 1, nbi
2836 N_i(n2) = N0_i*Di(n2)**mu_i * DEXP(-lami*Di(n2))*dti(n2)
2837 if (Di(n2).ge.D0s) then
2838 t1 = t1 + N_i(n2) * am_i*Di(n2)**bm_i
2839 t2 = t2 + N_i(n2)
2840 endif
2841 enddo
2842 endif
2843 tps_iaus(i,j) = t1
2844 tni_iaus(i,j) = t2
2845 enddo
2846 enddo
2847
2848 end subroutine qi_aut_qs
2849 !+---+-----------------------------------------------------------------+
2850 !+---+-----------------------------------------------------------------+
2851 SUBROUTINE GCF(GAMMCF,A,X,GLN)
2852 ! --- RETURNS THE INCOMPLETE GAMMA FUNCTION Q(A,X) EVALUATED BY ITS
2853 ! --- CONTINUED FRACTION REPRESENTATION AS GAMMCF. ALSO RETURNS
2854 ! --- LN(GAMMA(A)) AS GLN. THE CONTINUED FRACTION IS EVALUATED BY
2855 ! --- A MODIFIED LENTZ METHOD.
2856 ! --- USES GAMMLN
2857 IMPLICIT NONE
2858 INTEGER, PARAMETER:: ITMAX=100
2859 REAL, PARAMETER:: gEPS=3.E-7
2860 REAL, PARAMETER:: FPMIN=1.E-30
2861 REAL, INTENT(IN):: A, X
2862 REAL:: GAMMCF,GLN
2863 INTEGER:: I
2864 REAL:: AN,B,C,D,DEL,H
2865 GLN=GAMMLN(A)
2866 B=X+1.-A
2867 C=1./FPMIN
2868 D=1./B
2869 H=D
2870 DO 11 I=1,ITMAX
2871 AN=-I*(I-A)
2872 B=B+2.
2873 D=AN*D+B
2874 IF(ABS(D).LT.FPMIN)D=FPMIN
2875 C=B+AN/C
2876 IF(ABS(C).LT.FPMIN)C=FPMIN
2877 D=1./D
2878 DEL=D*C
2879 H=H*DEL
2880 IF(ABS(DEL-1.).LT.gEPS)GOTO 1
2881 11 CONTINUE
2882 PRINT *, 'A TOO LARGE, ITMAX TOO SMALL IN GCF'
2883 1 GAMMCF=EXP(-X+A*LOG(X)-GLN)*H
2884 END SUBROUTINE GCF
2885 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
2886 !+---+-----------------------------------------------------------------+
2887 SUBROUTINE GSER(GAMSER,A,X,GLN)
2888 ! --- RETURNS THE INCOMPLETE GAMMA FUNCTION P(A,X) EVALUATED BY ITS
2889 ! --- ITS SERIES REPRESENTATION AS GAMSER. ALSO RETURNS LN(GAMMA(A))
2890 ! --- AS GLN.
2891 ! --- USES GAMMLN
2892 IMPLICIT NONE
2893 INTEGER, PARAMETER:: ITMAX=100
2894 REAL, PARAMETER:: gEPS=3.E-7
2895 REAL, INTENT(IN):: A, X
2896 REAL:: GAMSER,GLN
2897 INTEGER:: N
2898 REAL:: AP,DEL,SUM
2899 GLN=GAMMLN(A)
2900 IF(X.LE.0.)THEN
2901 IF(X.LT.0.) PRINT *, 'X < 0 IN GSER'
2902 GAMSER=0.
2903 RETURN
2904 ENDIF
2905 AP=A
2906 SUM=1./A
2907 DEL=SUM
2908 DO 11 N=1,ITMAX
2909 AP=AP+1.
2910 DEL=DEL*X/AP
2911 SUM=SUM+DEL
2912 IF(ABS(DEL).LT.ABS(SUM)*gEPS)GOTO 1
2913 11 CONTINUE
2914 PRINT *,'A TOO LARGE, ITMAX TOO SMALL IN GSER'
2915 1 GAMSER=SUM*EXP(-X+A*LOG(X)-GLN)
2916 END SUBROUTINE GSER
2917 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
2918 !+---+-----------------------------------------------------------------+
2919 REAL FUNCTION GAMMLN(XX)
2920 ! --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 0.
2921 IMPLICIT NONE
2922 REAL, INTENT(IN):: XX
2923 DOUBLE PRECISION, PARAMETER:: STP = 2.5066282746310005D0
2924 DOUBLE PRECISION, DIMENSION(6), PARAMETER:: &
2925 COF = (/76.18009172947146D0, -86.50532032941677D0, &
2926 24.01409824083091D0, -1.231739572450155D0, &
2927 .1208650973866179D-2, -.5395239384953D-5/)
2928 DOUBLE PRECISION:: SER,TMP,X,Y
2929 INTEGER:: J
2930
2931 X=XX
2932 Y=X
2933 TMP=X+5.5D0
2934 TMP=(X+0.5D0)*LOG(TMP)-TMP
2935 SER=1.000000000190015D0
2936 DO 11 J=1,6
2937 Y=Y+1.D0
2938 SER=SER+COF(J)/Y
2939 11 CONTINUE
2940 GAMMLN=TMP+LOG(STP*SER/X)
2941 END FUNCTION GAMMLN
2942 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
2943 !+---+-----------------------------------------------------------------+
2944 REAL FUNCTION GAMMP(A,X)
2945 ! --- COMPUTES THE INCOMPLETE GAMMA FUNCTION P(A,X)
2946 ! --- SEE ABRAMOWITZ AND STEGUN 6.5.1
2947 ! --- USES GCF,GSER
2948 IMPLICIT NONE
2949 REAL, INTENT(IN):: A,X
2950 REAL:: GAMMCF,GAMSER,GLN
2951 GAMMP = 0.
2952 IF((X.LT.0.) .OR. (A.LE.0.)) THEN
2953 PRINT *, 'BAD ARGUMENTS IN GAMMP'
2954 RETURN
2955 ELSEIF(X.LT.A+1.)THEN
2956 CALL GSER(GAMSER,A,X,GLN)
2957 GAMMP=GAMSER
2958 ELSE
2959 CALL GCF(GAMMCF,A,X,GLN)
2960 GAMMP=1.-GAMMCF
2961 ENDIF
2962 END FUNCTION GAMMP
2963 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
2964 !+---+-----------------------------------------------------------------+
2965 REAL FUNCTION WGAMMA(y)
2966
2967 IMPLICIT NONE
2968 REAL, INTENT(IN):: y
2969
2970 WGAMMA = EXP(GAMMLN(y))
2971
2972 END FUNCTION WGAMMA
2973 !+---+-----------------------------------------------------------------+
2974 ! THIS FUNCTION CALCULATES THE LIQUID SATURATION VAPOR MIXING RATIO AS
2975 ! A FUNCTION OF TEMPERATURE AND PRESSURE
2976 !
2977 REAL FUNCTION RSLF(P,T)
2978
2979 IMPLICIT NONE
2980 REAL, INTENT(IN):: P, T
2981 REAL:: ESL,X
2982 REAL, PARAMETER:: C0= .611583699E03
2983 REAL, PARAMETER:: C1= .444606896E02
2984 REAL, PARAMETER:: C2= .143177157E01
2985 REAL, PARAMETER:: C3= .264224321E-1
2986 REAL, PARAMETER:: C4= .299291081E-3
2987 REAL, PARAMETER:: C5= .203154182E-5
2988 REAL, PARAMETER:: C6= .702620698E-8
2989 REAL, PARAMETER:: C7= .379534310E-11
2990 REAL, PARAMETER:: C8=-.321582393E-13
2991
2992 X=MAX(-80.,T-273.16)
2993
2994 ! ESL=612.2*EXP(17.67*X/(T-29.65))
2995 ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
2996 RSLF=.622*ESL/(P-ESL)
2997
2998 END FUNCTION RSLF
2999 !
3000 !+---+-----------------------------------------------------------------+
3001 ! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A
3002 ! FUNCTION OF TEMPERATURE AND PRESSURE
3003 !
3004 REAL FUNCTION RSIF(P,T)
3005
3006 IMPLICIT NONE
3007 REAL, INTENT(IN):: P, T
3008 REAL:: ESI,X
3009 REAL, PARAMETER:: C0= .609868993E03
3010 REAL, PARAMETER:: C1= .499320233E02
3011 REAL, PARAMETER:: C2= .184672631E01
3012 REAL, PARAMETER:: C3= .402737184E-1
3013 REAL, PARAMETER:: C4= .565392987E-3
3014 REAL, PARAMETER:: C5= .521693933E-5
3015 REAL, PARAMETER:: C6= .307839583E-7
3016 REAL, PARAMETER:: C7= .105785160E-9
3017 REAL, PARAMETER:: C8= .161444444E-12
3018
3019 X=MAX(-80.,T-273.16)
3020 ESI=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
3021 RSIF=.622*ESI/(P-ESI)
3022
3023 END FUNCTION RSIF
3024 !+---+-----------------------------------------------------------------+
3025 END MODULE module_mp_thompson
3026 !+---+-----------------------------------------------------------------+
3027 !
3028 ! MODIFICATIONS TO MAKE IN OTHER MODULES
3029 !
3030 ! Use this new code by changing the "THOMPSON" section of code found
3031 ! in "module_microphysics_driver.F" with this section. [Of course
3032 ! remove the leading comment character that you see here.]
3033 !
3034 ! CASE (THOMPSON)
3035 ! CALL wrf_debug ( 100 , 'microphysics_driver: calling thompson et al' )
3036 ! IF ( PRESENT( QV_CURR ) .AND. PRESENT ( QC_CURR ) .AND. &
3037 ! PRESENT( QR_CURR ) .AND. PRESENT ( QI_CURR ) .AND. &
3038 ! PRESENT( QS_CURR ) .AND. PRESENT ( QG_CURR ) .AND. &
3039 ! PRESENT ( QNI_CURR ).AND. &
3040 ! PRESENT( RAINNC ) .AND. PRESENT ( RAINNCV ) ) THEN
3041 ! CALL mp_gt_driver( &
3042 ! QV=qv_curr, &
3043 ! QC=qc_curr, &
3044 ! QR=qr_curr, &
3045 ! QI=qi_curr, &
3046 ! QS=qs_curr, &
3047 ! QG=qg_curr, &
3048 ! NI=qni_curr, &
3049 ! TH=th, &
3050 ! PII=pi_phy, &
3051 ! P=p, &
3052 ! DZ=dz8w, &
3053 ! DT_IN=dt, &
3054 ! ITIMESTEP=itimestep, &
3055 ! RAINNC=RAINNC, &
3056 ! RAINNCV=RAINNCV, &
3057 ! SR=SR &
3058 ! ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
3059 ! ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
3060 ! ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte)
3061 ! ELSE
3062 ! CALL wrf_error_fatal ( 'arguments not present for calling thompson_et_al' )
3063 ! ENDIF
3064 !
3065 ! Then rename the call from "thomp_init" to "thompson_init" in the file
3066 ! "module_physics_init.F" (seen below):
3067 !
3068 ! CASE (THOMPSON)
3069 ! CALL thompson_init