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