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