Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_mp_thompson.F @ 3094

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

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

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

Download in other formats:

Original Format