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