| blob | e81ae17df76b484f2d8cb76b5b4201a230e460af |
1 !+---+-----------------------------------------------------------------+
2 !.. This subroutine computes the moisture tendencies of water vapor,
3 !.. cloud droplets, rain, cloud ice (pristine), snow, and graupel.
4 !.. Prior to WRFv2.2 this code was based on Reisner et al (1998), but
5 !.. few of those pieces remain. A complete description is now found
6 !.. in Thompson et al. (2004, 2007).
7 !.. Most importantly, users may wish to modify the prescribed number of
8 !.. cloud droplets (Nt_c; see guidelines mentioned below). Otherwise,
9 !.. users may alter the rain and graupel size distribution parameters
10 !.. to use exponential (Marshal-Palmer) or generalized gamma shape.
11 !.. The snow field assumes a combination of two gamma functions (from
12 !.. Field et al. 2005) and would require significant modifications
13 !.. throughout the entire code to alter its shape as well as accretion
14 !.. rates. Users may also alter the constants used for density of rain,
15 !.. graupel, ice, and snow, but the latter is not constant when using
16 !.. Paul Field's snow distribution and moments methods. Other values
17 !.. users can modify include the constants for mass and/or velocity
18 !.. power law relations and assumed capacitances used in deposition/
19 !.. sublimation/evaporation/melting.
20 !.. Remaining values should probably be left alone.
21 !..
22 !..Author: Greg Thompson, NCAR-RAL, gthompsn@ucar.edu, 303-497-2805
23 !..Last modified: 26 Oct 2007
24 !+---+-----------------------------------------------------------------+
25 !wrft:model_layer:physics
26 !+---+-----------------------------------------------------------------+
27 !
28 MODULE module_mp_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 !+---+
305 !
307 CONTAINS
309 SUBROUTINE thompson07_init
311 IMPLICIT NONE
313 INTEGER:: i, j, k, m, n
315 !
316 ! jm allocate lookup tables
317 !
319 LOGICAL :: do_init
321 !
322 !jm allocate the lookup tables
323 !
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
362 !
363 ! jm
364 !
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 !+---+-----------------------------------------------------------------+
678 !
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
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
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 !+---+-----------------------------------------------------------------+
898 !
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 !+---+-----------------------------------------------------------------+
907 !
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
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 !+---+-----------------------------------------------------------------+
2565 !
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 !+---+-----------------------------------------------------------------+
2646 !
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
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 !+---+-----------------------------------------------------------------+
2788 !
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 !+---+-----------------------------------------------------------------+
2867 !
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
2927 !
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
3163 !
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
3186 !
3187 !+---+-----------------------------------------------------------------+
3188 ! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A
3189 ! FUNCTION OF TEMPERATURE AND PRESSURE
3190 !
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 !+---+-----------------------------------------------------------------+
3214 !
3215 ! MODIFICATIONS TO MAKE IN OTHER MODULES (pre v2.2 code only)
3216 !
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.]
3220 !
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
3251 !
3252 ! Then rename the call from "thomp_init" to "thompson_init" in the file
3253 ! "module_physics_init.F" (seen below):
3254 !
3255 ! CASE (THOMPSON)
3256 ! CALL thompson_init