Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_thompson.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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

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

Download in other formats:

Original Format