Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_thompson.F @ 3089

Last change on this file since 3089 was 11, checked in by aslmd, 15 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 114.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format