Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_thompson.F @ 3094

Last change on this file since 3094 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 118.0 KB

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

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

Download in other formats:

Original Format