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