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).
14 !.. Beginning with WRFv3.6, this is also the "aerosol-aware" scheme as
15 !.. described in Thompson, G. and T. Eidhammer, 2014: A study of
16 !.. aerosol impacts on clouds and precipitation development in a large
17 !.. winter cyclone. J. Atmos. Sci., 71, 3636-3658. Setting WRF
18 !.. namelist option mp_physics=8 utilizes the older one-moment cloud
19 !.. water with constant droplet concentration set as Nt_c (found below)
20 !.. while mp_physics=28 uses double-moment cloud droplet number
21 !.. concentration, which is not permitted to exceed Nt_c_max below.
23 !.. Most importantly, users may wish to modify the prescribed number of
24 !.. cloud droplets (Nt_c; see guidelines mentioned below). Otherwise,
25 !.. users may alter the rain and graupel size distribution parameters
26 !.. to use exponential (Marshal-Palmer) or generalized gamma shape.
27 !.. The snow field assumes a combination of two gamma functions (from
28 !.. Field et al. 2005) and would require significant modifications
29 !.. throughout the entire code to alter its shape as well as accretion
30 !.. rates. Users may also alter the constants used for density of rain,
31 !.. graupel, ice, and snow, but the latter is not constant when using
32 !.. Paul Field's snow distribution and moments methods. Other values
33 !.. users can modify include the constants for mass and/or velocity
34 !.. power law relations and assumed capacitances used in deposition/
35 !.. sublimation/evaporation/melting.
36 !.. Remaining values should probably be left alone.
38 !..Author: Greg Thompson, NCAR-RAL, gthompsn@ucar.edu, 303-497-2805
39 !..Last modified: 24 Jan 2018 Aerosol additions to v3.5.1 code 9/2013
40 !.. Cloud fraction additions 11/2014 part of pre-v3.7
41 !.. 3-moment graupel/hail additions Aug. 2018 (released in v4.4)
42 !+---+-----------------------------------------------------------------+
43 !wrft:model_layer:physics
44 !+---+-----------------------------------------------------------------+
48 MODULE module_mp_thompson
52 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
53 USE module_dm, ONLY : wrf_dm_max_real
55 USE module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG
59 LOGICAL, PARAMETER, PRIVATE:: iiwarm = .false.
60 LOGICAL, PRIVATE:: is_aerosol_aware = .false.
61 LOGICAL, PARAMETER, PRIVATE:: dustyIce = .true.
62 LOGICAL, PARAMETER, PRIVATE:: homogIce = .true.
63 LOGICAL, PRIVATE:: is_hail_aware = .false.
65 INTEGER, PARAMETER, PRIVATE:: IFDRY = 0
66 REAL, PARAMETER, PRIVATE:: T_0 = 273.15
67 REAL, PARAMETER, PRIVATE:: PI = 3.1415926536
69 !..Densities of rain, snow, graupel, and cloud ice.
70 REAL, PARAMETER, PRIVATE:: rho_w = 1000.0
71 REAL, PRIVATE:: rho_s = 100.0
72 REAL, PARAMETER, PRIVATE:: rho_i = 890.0
73 INTEGER, PARAMETER, PRIVATE:: NRHG = 9
74 INTEGER, PARAMETER, PRIVATE:: NRHG1 = 1
77 REAL, DIMENSION(NRHG), PARAMETER, PRIVATE:: &
78 rho_g = (/50., 100., 200., 300., 400., 500., 600., 700., 800./)
79 INTEGER, PARAMETER :: idx_bg1 = 5 ! index for rhog when mp=8 or 28
81 !..Prescribed number of cloud droplets. Set according to known data or
82 !.. roughly 100 per cc (100.E6 m^-3) for Maritime cases and
83 !.. 300 per cc (300.E6 m^-3) for Continental. Gamma shape parameter,
84 !.. mu_c, calculated based on Nt_c is important in autoconversion
85 !.. scheme. In 2-moment cloud water, Nt_c represents a maximum of
86 !.. droplet concentration and nu_c is also variable depending on local
87 !.. droplet number concentration.
88 REAL, PARAMETER, PRIVATE:: Nt_c = 100.E6
89 REAL, PARAMETER, PRIVATE:: Nt_c_max = 1999.E6
91 !..Declaration of constants for assumed CCN/IN aerosols when none in
92 !.. the input data. Look inside the init routine for modifications
93 !.. due to surface land-sea points or vegetation characteristics.
94 REAL, PARAMETER, PRIVATE:: naIN0 = 1.5E6
95 REAL, PARAMETER, PRIVATE:: naIN1 = 0.5E6
96 REAL, PARAMETER, PRIVATE:: naCCN0 = 300.0E6
97 REAL, PARAMETER, PRIVATE:: naCCN1 = 50.0E6
98 REAL, PARAMETER, PRIVATE:: naBC0 = 150.0E6
99 REAL, PARAMETER, PRIVATE:: naBC1 = 25.0E6
101 !..Generalized gamma distributions for rain, graupel and cloud ice.
102 !.. N(D) = N_0 * D**mu * exp(-lamda*D); mu=0 is exponential.
103 REAL, PARAMETER, PRIVATE:: mu_r = 0.0
104 REAL, PARAMETER, PRIVATE:: mu_g = 0.0
105 REAL, PARAMETER, PRIVATE:: mu_i = 0.0
108 !..Sum of two gamma distrib for snow (Field et al. 2005).
109 !.. N(D) = M2**4/M3**3 * [Kap0*exp(-M2*Lam0*D/M3)
110 !.. + Kap1*(M2/M3)**mu_s * D**mu_s * exp(-M2*Lam1*D/M3)]
111 !.. M2 and M3 are the (bm_s)th and (bm_s+1)th moments respectively
112 !.. calculated as function of ice water content and temperature.
113 REAL, PARAMETER, PRIVATE:: mu_s = 0.6357
114 REAL, PARAMETER, PRIVATE:: Kap0 = 490.6
115 REAL, PARAMETER, PRIVATE:: Kap1 = 17.46
116 REAL, PARAMETER, PRIVATE:: Lam0 = 20.78
117 REAL, PARAMETER, PRIVATE:: Lam1 = 3.29
119 !..Y-intercept parameter for graupel is not constant and depends on
120 !.. mixing ratio. Also, when mu_g is non-zero, these become equiv
121 !.. y-intercept for an exponential distrib and proper values are
122 !.. computed based on same mixing ratio and total number concentration.
123 REAL, PARAMETER, PRIVATE:: gonv_min = 1.E2
124 REAL, PARAMETER, PRIVATE:: gonv_max = 1.E6
126 !..Mass power law relations: mass = am*D**bm
127 !.. Snow from Field et al. (2005), others assume spherical form.
128 REAL, PARAMETER, PRIVATE:: am_r = PI*rho_w/6.0
129 REAL, PARAMETER, PRIVATE:: bm_r = 3.0
130 REAL, PARAMETER, PRIVATE:: am_s = 0.069
131 REAL, PARAMETER, PRIVATE:: bm_s = 2.0
132 REAL, DIMENSION(NRHG), PARAMETER,PRIVATE:: am_g = (/ &
133 & PI*rho_g(1)/6.0, PI*rho_g(2)/6.0, PI*rho_g(3)/6.0, &
134 & PI*rho_g(4)/6.0, PI*rho_g(5)/6.0, PI*rho_g(6)/6.0, &
135 & PI*rho_g(7)/6.0, PI*rho_g(8)/6.0, PI*rho_g(9)/6.0 /)
136 REAL, PARAMETER, PRIVATE:: bm_g = 3.0
137 REAL, PARAMETER, PRIVATE:: am_i = PI*rho_i/6.0
138 REAL, PARAMETER, PRIVATE:: bm_i = 3.0
140 !..Fallspeed power laws relations: v = (av*D**bv)*exp(-fv*D)
141 !.. Rain from Ferrier (1994), ice, snow, and graupel from
142 !.. Thompson et al (2008). Coefficient fv is zero for graupel/ice.
143 REAL, PARAMETER, PRIVATE:: av_r = 4854.0
144 REAL, PARAMETER, PRIVATE:: bv_r = 1.0
145 REAL, PARAMETER, PRIVATE:: fv_r = 195.0
146 REAL, PARAMETER, PRIVATE:: av_s = 40.0
147 REAL, PARAMETER, PRIVATE:: bv_s = 0.55
148 REAL, PARAMETER, PRIVATE:: fv_s = 100.0
149 REAL, PARAMETER, PRIVATE:: av_g_old = 442.0
150 REAL, PARAMETER, PRIVATE:: bv_g_old = 0.89
151 REAL, DIMENSION(NRHG), PRIVATE:: & ! Computed from A. Heymsfield: Best - Reynolds relation
152 & av_g = (/ 45.9173813, 67.0867386, 98.0158463, 122.353378, &
153 & 143.204224, 161.794724, 178.762115, 194.488785, &
155 REAL, DIMENSION(NRHG), PRIVATE:: & ! Computed from A. Heymsfield: Best - Reynolds relation
156 & bv_g = (/0.640961647, 0.640961647, 0.640961647, 0.640961647, &
157 & 0.640961647, 0.640961647, 0.640961647, 0.640961647, &
159 REAL, PARAMETER, PRIVATE:: a_coeff = 0.47244157
160 REAL, PARAMETER, PRIVATE:: b_coeff = 0.54698726
161 REAL, PARAMETER, PRIVATE:: av_i = 1493.9
162 REAL, PARAMETER, PRIVATE:: bv_i = 1.0
163 REAL, PARAMETER, PRIVATE:: av_c = 0.316946E8
164 REAL, PARAMETER, PRIVATE:: bv_c = 2.0
166 !..Capacitance of sphere and plates/aggregates: D**3, D**2
167 REAL, PARAMETER, PRIVATE:: C_cube = 0.5
168 REAL, PARAMETER, PRIVATE:: C_sqrd = 0.15
170 !..Collection efficiencies. Rain/snow/graupel collection of cloud
171 !.. droplets use variables (Ef_rw, Ef_sw, Ef_gw respectively) and
172 !.. get computed elsewhere because they are dependent on stokes
174 REAL, PARAMETER, PRIVATE:: Ef_si = 0.05
175 REAL, PARAMETER, PRIVATE:: Ef_rs = 0.95
176 REAL, PARAMETER, PRIVATE:: Ef_rg = 0.75
177 REAL, PARAMETER, PRIVATE:: Ef_ri = 0.95
179 !..Minimum microphys values
180 !.. R1 value, 1.E-12, cannot be set lower because of numerical
181 !.. problems with Paul Field's moments and should not be set larger
182 !.. because of truncation problems in snow/ice growth.
183 REAL, PARAMETER, PRIVATE:: R1 = 1.E-12
184 REAL, PARAMETER, PRIVATE:: R2 = 1.E-6
185 REAL, PARAMETER, PRIVATE:: eps = 1.E-15
187 !..Constants in Cooper curve relation for cloud ice number.
188 REAL, PARAMETER, PRIVATE:: TNO = 5.0
189 REAL, PARAMETER, PRIVATE:: ATO = 0.304
191 !..Rho_not used in fallspeed relations (rho_not/rho)**.5 adjustment.
192 REAL, PARAMETER, PRIVATE:: rho_not = 101325.0/(287.05*298.0)
195 REAL, PARAMETER, PRIVATE:: Sc = 0.632
198 !..Homogeneous freezing temperature
199 REAL, PARAMETER, PRIVATE:: HGFR = 235.16
201 !..Water vapor and air gas constants at constant pressure
202 REAL, PARAMETER, PRIVATE:: Rv = 461.5
203 REAL, PARAMETER, PRIVATE:: oRv = 1./Rv
204 REAL, PARAMETER, PRIVATE:: R = 287.04
205 REAL, PARAMETER, PRIVATE:: Cp = 1004.0
206 REAL, PARAMETER, PRIVATE:: R_uni = 8.314 ! J (mol K)-1
208 DOUBLE PRECISION, PARAMETER, PRIVATE:: k_b = 1.38065E-23 ! Boltzmann constant [J/K]
209 DOUBLE PRECISION, PARAMETER, PRIVATE:: M_w = 18.01528E-3 ! molecular mass of water [kg/mol]
210 DOUBLE PRECISION, PARAMETER, PRIVATE:: M_a = 28.96E-3 ! molecular mass of air [kg/mol]
211 DOUBLE PRECISION, PARAMETER, PRIVATE:: N_avo = 6.022E23 ! Avogadro number [1/mol]
212 DOUBLE PRECISION, PARAMETER, PRIVATE:: ma_w = M_w / N_avo ! mass of water molecule [kg]
213 REAL, PARAMETER, PRIVATE:: ar_volume = 4./3.*PI*(2.5e-6)**3 ! assume radius of 0.025 micrometer, 2.5e-6 cm
215 !..Enthalpy of sublimation, vaporization, and fusion at 0C.
216 REAL, PARAMETER, PRIVATE:: lsub = 2.834E6
217 REAL, PARAMETER, PRIVATE:: lvap0 = 2.5E6
218 REAL, PARAMETER, PRIVATE:: lfus = lsub - lvap0
219 REAL, PARAMETER, PRIVATE:: olfus = 1./lfus
221 !..Ice initiates with this mass (kg), corresponding diameter calc.
222 !..Min diameters and mass of cloud, rain, snow, and graupel (m, kg).
223 REAL, PARAMETER, PRIVATE:: xm0i = 1.E-12
224 REAL, PARAMETER, PRIVATE:: D0c = 1.E-6
225 REAL, PARAMETER, PRIVATE:: D0r = 50.E-6
226 REAL, PARAMETER, PRIVATE:: D0s = 300.E-6
227 REAL, PARAMETER, PRIVATE:: D0g = 350.E-6
228 REAL, PRIVATE:: D0i, xm0s, xm0g
230 !..Lookup table dimensions
231 INTEGER, PARAMETER, PRIVATE:: nbins = 100
232 INTEGER, PARAMETER, PRIVATE:: nbc = nbins
233 INTEGER, PARAMETER, PRIVATE:: nbi = nbins
234 INTEGER, PARAMETER, PRIVATE:: nbr = nbins
235 INTEGER, PARAMETER, PRIVATE:: nbs = nbins
236 INTEGER, PARAMETER, PRIVATE:: nbg = nbins
237 INTEGER, PARAMETER, PRIVATE:: ntb_c = 37
238 INTEGER, PARAMETER, PRIVATE:: ntb_i = 64
239 INTEGER, PARAMETER, PRIVATE:: ntb_r = 37
240 INTEGER, PARAMETER, PRIVATE:: ntb_s = 37
241 INTEGER, PARAMETER, PRIVATE:: ntb_g = 37
242 INTEGER, PARAMETER, PRIVATE:: ntb_g1 = 37
243 INTEGER, PARAMETER, PRIVATE:: ntb_r1 = 37
244 INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55
245 INTEGER, PARAMETER, PRIVATE:: ntb_t = 9
246 INTEGER, PRIVATE:: nic1, nic2, nii2, nii3, nir2, nir3, nis2, nig2, nig3
247 INTEGER, PARAMETER, PRIVATE:: ntb_arc = 7
248 INTEGER, PARAMETER, PRIVATE:: ntb_arw = 9
249 INTEGER, PARAMETER, PRIVATE:: ntb_art = 7
250 INTEGER, PARAMETER, PRIVATE:: ntb_arr = 5
251 INTEGER, PARAMETER, PRIVATE:: ntb_ark = 4
252 INTEGER, PARAMETER, PRIVATE:: ntb_IN = 55
253 INTEGER, PRIVATE:: niIN2
255 DOUBLE PRECISION, DIMENSION(nbins+1):: xDx
256 DOUBLE PRECISION, DIMENSION(nbc):: Dc, dtc
257 DOUBLE PRECISION, DIMENSION(nbi):: Di, dti
258 DOUBLE PRECISION, DIMENSION(nbr):: Dr, dtr
259 DOUBLE PRECISION, DIMENSION(nbs):: Ds, dts
260 DOUBLE PRECISION, DIMENSION(nbg):: Dg, dtg
261 DOUBLE PRECISION, DIMENSION(nbc):: t_Nc
263 !..Lookup tables for cloud water content (kg/m**3).
264 REAL, DIMENSION(ntb_c), PARAMETER, PRIVATE:: &
265 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, &
266 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, &
267 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, &
268 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, &
271 !..Lookup tables for cloud ice content (kg/m**3).
272 REAL, DIMENSION(ntb_i), PARAMETER, PRIVATE:: &
273 r_i = (/1.e-10,2.e-10,3.e-10,4.e-10, &
274 5.e-10,6.e-10,7.e-10,8.e-10,9.e-10, &
275 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, &
276 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, &
277 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, &
278 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, &
279 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, &
280 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, &
283 !..Lookup tables for rain content (kg/m**3).
284 REAL, DIMENSION(ntb_r), PARAMETER, PRIVATE:: &
285 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, &
286 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, &
287 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, &
288 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, &
291 !..Lookup tables for graupel content (kg/m**3).
292 REAL, DIMENSION(ntb_g), PARAMETER, PRIVATE:: &
293 r_g = (/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, &
294 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, &
295 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, &
296 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, &
299 !..Lookup tables for snow content (kg/m**3).
300 REAL, DIMENSION(ntb_s), PARAMETER, PRIVATE:: &
301 r_s = (/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, &
302 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, &
303 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, &
304 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, &
307 !..Lookup tables for rain y-intercept parameter (/m**4).
308 REAL, DIMENSION(ntb_r1), PARAMETER, PRIVATE:: &
309 N0r_exp = (/1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, &
310 1.e7,2.e7,3.e7,4.e7,5.e7,6.e7,7.e7,8.e7,9.e7, &
311 1.e8,2.e8,3.e8,4.e8,5.e8,6.e8,7.e8,8.e8,9.e8, &
312 1.e9,2.e9,3.e9,4.e9,5.e9,6.e9,7.e9,8.e9,9.e9, &
315 !..Lookup tables for graupel y-intercept parameter (/m**4).
316 REAL, DIMENSION(ntb_g1), PARAMETER, PRIVATE:: &
317 N0g_exp = (/1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, &
318 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, &
319 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
320 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
323 !..Lookup tables for ice number concentration (/m**3).
324 REAL, DIMENSION(ntb_i1), PARAMETER, PRIVATE:: &
325 Nt_i = (/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0, &
326 1.e1,2.e1,3.e1,4.e1,5.e1,6.e1,7.e1,8.e1,9.e1, &
327 1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, &
328 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, &
329 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
330 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
333 !..Aerosol table parameter: Number of available aerosols, vertical
334 !.. velocity, temperature, aerosol mean radius, and hygroscopicity.
335 REAL, DIMENSION(ntb_arc), PARAMETER, PRIVATE:: &
336 ta_Na = (/10.0, 31.6, 100.0, 316.0, 1000.0, 3160.0, 10000.0/)
337 REAL, DIMENSION(ntb_arw), PARAMETER, PRIVATE:: &
338 ta_Ww = (/0.01, 0.0316, 0.1, 0.316, 1.0, 3.16, 10.0, 31.6, 100.0/)
339 REAL, DIMENSION(ntb_art), PARAMETER, PRIVATE:: &
340 ta_Tk = (/243.15, 253.15, 263.15, 273.15, 283.15, 293.15, 303.15/)
341 REAL, DIMENSION(ntb_arr), PARAMETER, PRIVATE:: &
342 ta_Ra = (/0.01, 0.02, 0.04, 0.08, 0.16/)
343 REAL, DIMENSION(ntb_ark), PARAMETER, PRIVATE:: &
344 ta_Ka = (/0.2, 0.4, 0.6, 0.8/)
346 !..Lookup tables for IN concentration (/m**3) from 0.001 to 1000/Liter.
347 REAL, DIMENSION(ntb_IN), PARAMETER, PRIVATE:: &
348 Nt_IN = (/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0, &
349 1.e1,2.e1,3.e1,4.e1,5.e1,6.e1,7.e1,8.e1,9.e1, &
350 1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, &
351 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, &
352 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
353 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
356 !..For snow moments conversions (from Field et al. 2005)
357 REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
358 sa = (/ 5.065339, -0.062659, -3.032362, 0.029469, -0.000285, &
359 0.31255, 0.000204, 0.003199, 0.0, -0.015952/)
360 REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
361 sb = (/ 0.476221, -0.015896, 0.165977, 0.007468, -0.000141, &
362 0.060366, 0.000079, 0.000594, 0.0, -0.003577/)
364 !..Temperatures (5 C interval 0 to -40) used in lookup tables.
365 REAL, DIMENSION(ntb_t), PARAMETER, PRIVATE:: &
366 Tc = (/-0.01, -5., -10., -15., -20., -25., -30., -35., -40./)
368 !..Lookup tables for various accretion/collection terms.
369 !.. ntb_x refers to the number of elements for rain, snow, graupel,
370 !.. and temperature array indices. Variables beginning with t-p/c/m/n
371 !.. represent lookup tables. Save compile-time memory by making
372 !.. allocatable (2009Jun12, J. Michalakes).
373 INTEGER, PARAMETER, PRIVATE:: R8SIZE = 8
374 INTEGER, PARAMETER, PRIVATE:: R4SIZE = 4
375 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:,:):: &
376 tcg_racg, tmr_racg, tcr_gacr, & ! tmg_gacr
378 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: &
379 tcs_racs1, tmr_racs1, tcs_racs2, tmr_racs2, &
380 tcr_sacr1, tms_sacr1, tcr_sacr2, tms_sacr2, &
381 tnr_racs1, tnr_racs2, tnr_sacr1, tnr_sacr2
382 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: &
384 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: &
385 tpi_qrfz, tpg_qrfz, tni_qrfz, tnr_qrfz
386 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: &
387 tps_iaus, tni_iaus, tpi_ide
388 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efrw
389 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efsw
390 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: tnr_rev
391 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: &
393 REAL (KIND=R4SIZE), ALLOCATABLE, DIMENSION(:,:,:,:,:):: tnccn_act
395 !..Variables holding a bunch of exponents and gamma values (cloud water,
396 !.. cloud ice, rain, snow, then graupel).
397 REAL, DIMENSION(5,15), PRIVATE:: cce, ccg
398 REAL, DIMENSION(15), PRIVATE:: ocg1, ocg2
399 REAL, DIMENSION(7), PRIVATE:: cie, cig
400 REAL, PRIVATE:: oig1, oig2, obmi
401 REAL, DIMENSION(13), PRIVATE:: cre, crg
402 REAL, PRIVATE:: ore1, org1, org2, org3, obmr
403 REAL, DIMENSION(17), PRIVATE:: cse, csg
404 REAL, PRIVATE:: oams, obms, ocms
405 REAL, DIMENSION(12,NRHG), PRIVATE:: cge, cgg
406 REAL, PRIVATE:: oge1, ogg1, ogg2, ogg3, obmg
407 REAL, DIMENSION(NRHG), PRIVATE:: oamg, ocmg
409 !..Declaration of precomputed constants in various rate eqns.
410 REAL:: t1_qr_qc, t1_qr_qi, t2_qr_qi, t1_qg_qc, t1_qs_qc, t1_qs_qi
411 REAL:: t1_qr_ev, t2_qr_ev
412 REAL:: t1_qs_sd, t2_qs_sd, t1_qg_sd, t2_qg_sd
413 REAL:: t1_qs_me, t2_qs_me, t1_qg_me, t2_qg_me
416 !+---+-----------------------------------------------------------------+
418 !+---+-----------------------------------------------------------------+
424 SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, ng, &
425 nwfa, nifa, nbca, wif_input_opt, frc_urb2d, &
427 ids, ide, jds, jde, kds, kde, &
428 ims, ime, jms, jme, kms, kme, &
429 its, ite, jts, jte, kts, kte)
433 INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, &
434 ims,ime, jms,jme, kms,kme, &
435 its,ite, jts,jte, kts,kte
436 REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN):: hgt
438 !..OPTIONAL variables that control application of hail-aware scheme
439 REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(INOUT) :: ng
441 !..OPTIONAL variables that control application of aerosol-aware scheme
443 REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(INOUT) :: nwfa, nifa, nbca
444 REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(INOUT) :: nwfa2d, nbca2d
445 REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(IN) :: orho
446 REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(IN) :: frc_urb2d
447 REAL, OPTIONAL, INTENT(IN) :: DX, DY
448 LOGICAL, OPTIONAL, INTENT(IN) :: is_start
449 INTEGER, OPTIONAL, INTENT(IN) :: wif_input_opt
450 CHARACTER*256:: mp_debug
453 INTEGER:: i, j, k, l, m, n
454 REAL:: h_01, niIN3, niCCN3, niBC3, max_test, z1
455 LOGICAL:: micro_init, has_CCN, has_IN
459 if (PRESENT(ng)) then
460 is_hail_aware = .TRUE.
463 av_g(idx_bg1) = av_g_old
464 bv_g(idx_bg1) = bv_g_old
468 is_aerosol_aware = .FALSE.
473 write(mp_debug,*) ' DEBUG checking column of hgt ', its+1,jts+1
474 CALL wrf_debug(250, mp_debug)
476 write(mp_debug,*) ' DEBUGT k, hgt = ', k, hgt(its+1,k,jts+1)
477 CALL wrf_debug(250, mp_debug)
480 if (PRESENT(nwfa2d) .AND. PRESENT(nwfa) .AND. PRESENT(nifa)) is_aerosol_aware = .TRUE.
482 if (is_aerosol_aware) then
484 !..Check for existing aerosol data, both CCN and IN aerosols. If missing
485 !.. fill in just a basic vertical profile, somewhat boundary-layer following.
487 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
488 max_test = wrf_dm_max_real ( MAXVAL(nwfa(its:ite-1,:,jts:jte-1)) )
490 max_test = MAXVAL ( nwfa(its:ite-1,:,jts:jte-1) )
493 if (max_test .lt. eps) then
494 write(mp_debug,*) ' Apparently there are no initial CCN aerosols.'
495 CALL wrf_debug(100, mp_debug)
496 write(mp_debug,*) ' checked column at point (i,j) = ', its,jts
497 CALL wrf_debug(100, mp_debug)
498 do j = jts, min(jde-1,jte)
499 do i = its, min(ide-1,ite)
500 if (hgt(i,1,j).le.1000.0) then
502 elseif (hgt(i,1,j).ge.2500.0) then
505 h_01 = 0.8*cos(hgt(i,1,j)*0.001 - 1.0)
507 niCCN3 = -1.0*ALOG(naCCN1/naCCN0)/h_01
508 nwfa(i,1,j) = naCCN1+naCCN0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niCCN3)
509 z1=hgt(i,2,j)-hgt(i,1,j)
510 nwfa2d(i,j) = nwfa(i,1,j) * 0.000196 * (50./z1)
512 nwfa(i,k,j) = naCCN1+naCCN0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niCCN3)
518 write(mp_debug,*) ' Apparently initial CCN aerosols are present.'
519 CALL wrf_debug(100, mp_debug)
520 write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nwfa(its,:,jts))
521 CALL wrf_debug(100, mp_debug)
525 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
526 max_test = wrf_dm_max_real ( MAXVAL(nifa(its:ite-1,:,jts:jte-1)) )
528 max_test = MAXVAL ( nifa(its:ite-1,:,jts:jte-1) )
531 if (max_test .lt. eps) then
532 write(mp_debug,*) ' Apparently there are no initial IN aerosols.'
533 CALL wrf_debug(100, mp_debug)
534 write(mp_debug,*) ' checked column at point (i,j) = ', its,jts
535 CALL wrf_debug(100, mp_debug)
536 do j = jts, min(jde-1,jte)
537 do i = its, min(ide-1,ite)
538 if (hgt(i,1,j).le.1000.0) then
540 elseif (hgt(i,1,j).ge.2500.0) then
543 h_01 = 0.8*cos(hgt(i,1,j)*0.001 - 1.0)
545 niIN3 = -1.0*ALOG(naIN1/naIN0)/h_01
546 nifa(i,1,j) = naIN1+naIN0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niIN3)
548 nifa(i,k,j) = naIN1+naIN0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niIN3)
554 write(mp_debug,*) ' Apparently initial IN aerosols are present.'
555 CALL wrf_debug(100, mp_debug)
556 write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nifa(its,:,jts))
557 CALL wrf_debug(100, mp_debug)
561 if ( wif_input_opt .eq. 2 ) then
563 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
564 max_test = wrf_dm_max_real ( MAXVAL(nbca(its:ite-1,:,jts:jte-1)) )
566 max_test = MAXVAL ( nbca(its:ite-1,:,jts:jte-1) )
569 if (max_test .lt. eps) then
570 write(mp_debug,*) ' Apparently there are no initial BC aerosols.'
571 CALL wrf_debug(100, mp_debug)
572 write(mp_debug,*) ' checked column at point (i,j) = ', its,jts
573 CALL wrf_debug(100, mp_debug)
574 do j = jts, min(jde-1,jte)
575 do i = its, min(ide-1,ite)
576 if (hgt(i,1,j).le.1000.0) then
578 elseif (hgt(i,1,j).ge.2500.0) then
581 h_01 = 0.8*cos(hgt(i,1,j)*0.001 - 1.0)
583 niBC3 = -1.0*ALOG(naBC1/naBC0)/h_01
584 nbca(i,1,j) = naBC1+naBC0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niBC3)
585 z1=hgt(i,2,j)-hgt(i,1,j)
586 nbca2d(i,j) = nbca(i,1,j) * 0.000098 * (50./z1) * (1. + frc_urb2d(i,j))
588 nbca(i,k,j) = naBC1+naBC0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niBC3)
593 write(mp_debug,*) ' Apparently initial BC aerosols are present.'
594 CALL wrf_debug(100, mp_debug)
595 write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nbca(its,:,jts))
596 CALL wrf_debug(100, mp_debug)
604 !..Allocate space for lookup tables (J. Michalakes 2009Jun08).
606 if (.NOT. ALLOCATED(tcg_racg) ) then
607 ALLOCATE(tcg_racg(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r))
611 if (.NOT. ALLOCATED(tmr_racg)) ALLOCATE(tmr_racg(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r))
612 if (.NOT. ALLOCATED(tcr_gacr)) ALLOCATE(tcr_gacr(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r))
613 ! if (.NOT. ALLOCATED(tmg_gacr)) ALLOCATE(tmg_gacr(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r))
614 if (.NOT. ALLOCATED(tnr_racg)) ALLOCATE(tnr_racg(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r))
615 if (.NOT. ALLOCATED(tnr_gacr)) ALLOCATE(tnr_gacr(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r))
617 if (.NOT. ALLOCATED(tcs_racs1)) ALLOCATE(tcs_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
618 if (.NOT. ALLOCATED(tmr_racs1)) ALLOCATE(tmr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
619 if (.NOT. ALLOCATED(tcs_racs2)) ALLOCATE(tcs_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
620 if (.NOT. ALLOCATED(tmr_racs2)) ALLOCATE(tmr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
621 if (.NOT. ALLOCATED(tcr_sacr1)) ALLOCATE(tcr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
622 if (.NOT. ALLOCATED(tms_sacr1)) ALLOCATE(tms_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
623 if (.NOT. ALLOCATED(tcr_sacr2)) ALLOCATE(tcr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
624 if (.NOT. ALLOCATED(tms_sacr2)) ALLOCATE(tms_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
625 if (.NOT. ALLOCATED(tnr_racs1)) ALLOCATE(tnr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
626 if (.NOT. ALLOCATED(tnr_racs2)) ALLOCATE(tnr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
627 if (.NOT. ALLOCATED(tnr_sacr1)) ALLOCATE(tnr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
628 if (.NOT. ALLOCATED(tnr_sacr2)) ALLOCATE(tnr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
630 if (.NOT. ALLOCATED(tpi_qcfz)) ALLOCATE(tpi_qcfz(ntb_c,nbc,45,ntb_IN))
631 if (.NOT. ALLOCATED(tni_qcfz)) ALLOCATE(tni_qcfz(ntb_c,nbc,45,ntb_IN))
633 if (.NOT. ALLOCATED(tpi_qrfz)) ALLOCATE(tpi_qrfz(ntb_r,ntb_r1,45,ntb_IN))
634 if (.NOT. ALLOCATED(tpg_qrfz)) ALLOCATE(tpg_qrfz(ntb_r,ntb_r1,45,ntb_IN))
635 if (.NOT. ALLOCATED(tni_qrfz)) ALLOCATE(tni_qrfz(ntb_r,ntb_r1,45,ntb_IN))
636 if (.NOT. ALLOCATED(tnr_qrfz)) ALLOCATE(tnr_qrfz(ntb_r,ntb_r1,45,ntb_IN))
638 if (.NOT. ALLOCATED(tps_iaus)) ALLOCATE(tps_iaus(ntb_i,ntb_i1))
639 if (.NOT. ALLOCATED(tni_iaus)) ALLOCATE(tni_iaus(ntb_i,ntb_i1))
640 if (.NOT. ALLOCATED(tpi_ide)) ALLOCATE(tpi_ide(ntb_i,ntb_i1))
642 if (.NOT. ALLOCATED(t_Efrw)) ALLOCATE(t_Efrw(nbr,nbc))
643 if (.NOT. ALLOCATED(t_Efsw)) ALLOCATE(t_Efsw(nbs,nbc))
645 if (.NOT. ALLOCATED(tnr_rev)) ALLOCATE(tnr_rev(nbr, ntb_r1, ntb_r))
646 if (.NOT. ALLOCATED(tpc_wev)) ALLOCATE(tpc_wev(nbc,ntb_c,nbc))
647 if (.NOT. ALLOCATED(tnc_wev)) ALLOCATE(tnc_wev(nbc,ntb_c,nbc))
649 if (.NOT. ALLOCATED(tnccn_act)) &
650 ALLOCATE(tnccn_act(ntb_arc,ntb_arw,ntb_art,ntb_arr,ntb_ark))
654 !..From Martin et al. (1994), assign gamma shape parameter mu for cloud
655 !.. drops according to general dispersion characteristics (disp=~0.25
656 !.. for Maritime and 0.45 for Continental).
657 !.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime
658 !.. to 2 for really dirty air. This not used in 2-moment cloud water
659 !.. scheme and nu_c used instead and varies from 2 to 15 (integer-only).
660 mu_c = MIN(15., (1000.E6/Nt_c + 2.))
662 !..Schmidt number to one-third used numerous times.
665 !..Compute min ice diam from mass, min snow/graupel mass from diam.
666 D0i = (xm0i/am_i)**(1./bm_i)
667 xm0s = am_s * D0s**bm_s
668 xm0g = am_g(NRHG) * D0g**bm_g
670 !..These constants various exponents and gamma() assoc with cloud,
671 !.. rain, snow, and graupel.
674 cce(2,n) = bm_r + n + 1.
675 cce(3,n) = bm_r + n + 4.
676 cce(4,n) = n + bv_c + 1.
677 cce(5,n) = bm_r + n + bv_c + 1.
678 ccg(1,n) = WGAMMA(cce(1,n))
679 ccg(2,n) = WGAMMA(cce(2,n))
680 ccg(3,n) = WGAMMA(cce(3,n))
681 ccg(4,n) = WGAMMA(cce(4,n))
682 ccg(5,n) = WGAMMA(cce(5,n))
683 ocg1(n) = 1./ccg(1,n)
684 ocg2(n) = 1./ccg(2,n)
688 cie(2) = bm_i + mu_i + 1.
689 cie(3) = bm_i + mu_i + bv_i + 1.
690 cie(4) = mu_i + bv_i + 1.
692 cie(6) = bm_i*0.5 + mu_i + bv_i + 1.
693 cie(7) = bm_i*0.5 + mu_i + 1.
694 cig(1) = WGAMMA(cie(1))
695 cig(2) = WGAMMA(cie(2))
696 cig(3) = WGAMMA(cie(3))
697 cig(4) = WGAMMA(cie(4))
698 cig(5) = WGAMMA(cie(5))
699 cig(6) = WGAMMA(cie(6))
700 cig(7) = WGAMMA(cie(7))
707 cre(3) = bm_r + mu_r + 1.
708 cre(4) = bm_r*2. + mu_r + 1.
709 cre(5) = mu_r + bv_r + 1.
710 cre(6) = bm_r + mu_r + bv_r + 1.
711 cre(7) = bm_r*0.5 + mu_r + bv_r + 1.
712 cre(8) = bm_r + mu_r + bv_r + 3.
713 cre(9) = mu_r + bv_r + 3.
715 cre(11) = 0.5*(bv_r + 5. + 2.*mu_r)
716 cre(12) = bm_r*0.5 + mu_r + 1.
717 cre(13) = bm_r*2. + mu_r + bv_r + 1.
719 crg(n) = WGAMMA(cre(n))
730 cse(4) = bm_s + bv_s + 1.
731 cse(5) = bm_s*2. + bv_s + 1.
732 cse(6) = bm_s*2. + 1.
733 cse(7) = bm_s + mu_s + 1.
734 cse(8) = bm_s + mu_s + 2.
735 cse(9) = bm_s + mu_s + 3.
736 cse(10) = bm_s + mu_s + bv_s + 1.
737 cse(11) = bm_s*2. + mu_s + bv_s + 1.
738 cse(12) = bm_s*2. + mu_s + 1.
740 cse(14) = bm_s + bv_s
742 cse(16) = 1.0 + (1.0 + bv_s)/2.
743 cse(17) = bm_s + bv_s + 2.
745 csg(n) = WGAMMA(cse(n))
753 cge(3,:) = bm_g + mu_g + 1.
754 cge(4,:) = bm_g*2. + mu_g + 1.
755 cge(10,:) = mu_g + 2.
756 cge(12,:) = bm_g*0.5 + mu_g + 1.
758 cge(5,m) = bm_g*2. + mu_g + bv_g(m) + 1.
759 cge(6,m) = bm_g + mu_g + bv_g(m) + 1.
760 cge(7,m) = bm_g*0.5 + mu_g + bv_g(m) + 1.
761 cge(8,m) = mu_g + bv_g(m) + 1. ! not used
762 cge(9,m) = mu_g + bv_g(m) + 3.
763 cge(11,m) = 0.5*(bv_g(m) + 5. + 2.*mu_g)
767 cgg(n,m) = WGAMMA(cge(n,m))
774 ocmg(m) = oamg(m)**obmg
781 !+---+-----------------------------------------------------------------+
782 !..Simplify various rate eqns the best we can now.
783 !+---+-----------------------------------------------------------------+
785 !..Rain collecting cloud water and cloud ice
786 t1_qr_qc = PI*.25*av_r * crg(9)
787 t1_qr_qi = PI*.25*av_r * crg(9)
788 t2_qr_qi = PI*.25*am_r*av_r * crg(8)
790 !..Graupel collecting cloud water
791 ! t1_qg_qc = PI*.25*av_g * cgg(9)
793 !..Snow collecting cloud water
794 t1_qs_qc = PI*.25*av_s
796 !..Snow collecting cloud ice
797 t1_qs_qi = PI*.25*av_s
799 !..Evaporation of rain; ignore depositional growth of rain.
800 t1_qr_ev = 0.78 * crg(10)
801 t2_qr_ev = 0.308*Sc3*SQRT(av_r) * crg(11)
803 !..Sublimation/depositional growth of snow
805 t2_qs_sd = 0.28*Sc3*SQRT(av_s)
808 t1_qs_me = PI*4.*C_sqrd*olfus * 0.86
809 t2_qs_me = PI*4.*C_sqrd*olfus * 0.28*Sc3*SQRT(av_s)
811 !..Sublimation/depositional growth of graupel
812 t1_qg_sd = 0.86 * cgg(10,1)
813 ! t2_qg_sd = 0.28*Sc3*SQRT(av_g) * cgg(11)
815 !..Melting of graupel
816 t1_qg_me = PI*4.*C_cube*olfus * 0.86 * cgg(10,1)
817 ! t2_qg_me = PI*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g) * cgg(11)
819 !..Constants for helping find lookup table indexes.
820 nic2 = NINT(ALOG10(r_c(1)))
821 nii2 = NINT(ALOG10(r_i(1)))
822 nii3 = NINT(ALOG10(Nt_i(1)))
823 nir2 = NINT(ALOG10(r_r(1)))
824 nir3 = NINT(ALOG10(N0r_exp(1)))
825 nis2 = NINT(ALOG10(r_s(1)))
826 nig2 = NINT(ALOG10(r_g(1)))
827 nig3 = NINT(ALOG10(N0g_exp(1)))
828 niIN2 = NINT(ALOG10(Nt_IN(1)))
830 !..Create bins of cloud water (from min diameter up to 100 microns).
834 Dc(n) = Dc(n-1) + 1.0D-6
835 dtc(n) = (Dc(n) - Dc(n-1))
838 !..Create bins of cloud ice (from min diameter up to 5x min snow size).
840 xDx(nbi+1) = 5.0d0*D0s
842 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbi) &
843 *DLOG(xDx(nbi+1)/xDx(1)) +DLOG(xDx(1)))
846 Di(n) = DSQRT(xDx(n)*xDx(n+1))
847 dti(n) = xDx(n+1) - xDx(n)
850 !..Create bins of rain (from min diameter up to 5 mm).
854 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbr) &
855 *DLOG(xDx(nbr+1)/xDx(1)) +DLOG(xDx(1)))
858 Dr(n) = DSQRT(xDx(n)*xDx(n+1))
859 dtr(n) = xDx(n+1) - xDx(n)
862 !..Create bins of snow (from min diameter up to 2 cm).
866 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbs) &
867 *DLOG(xDx(nbs+1)/xDx(1)) +DLOG(xDx(1)))
870 Ds(n) = DSQRT(xDx(n)*xDx(n+1))
871 dts(n) = xDx(n+1) - xDx(n)
874 !..Create bins of graupel (from min diameter up to 5 cm).
878 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbg) &
879 *DLOG(xDx(nbg+1)/xDx(1)) +DLOG(xDx(1)))
882 Dg(n) = DSQRT(xDx(n)*xDx(n+1))
883 dtg(n) = xDx(n+1) - xDx(n)
886 !..Create bins of cloud droplet number concentration (1 to 3000 per cc).
888 xDx(nbc+1) = 3000.0d0
890 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbc) &
891 *DLOG(xDx(nbc+1)/xDx(1)) +DLOG(xDx(1)))
894 t_Nc(n) = DSQRT(xDx(n)*xDx(n+1)) * 1.D6
896 nic1 = DLOG(t_Nc(nbc)/t_Nc(1))
898 !+---+-----------------------------------------------------------------+
899 !..Create lookup tables for most costly calculations.
900 !+---+-----------------------------------------------------------------+
907 tcg_racg(i,j,n,k,m) = 0.0d0
908 tmr_racg(i,j,n,k,m) = 0.0d0
909 tcr_gacr(i,j,n,k,m) = 0.0d0
910 !tmg_gacr(i,j,n,k,m) = 0.0d0
911 tnr_racg(i,j,n,k,m) = 0.0d0
912 tnr_gacr(i,j,n,k,m) = 0.0d0
923 tcs_racs1(i,j,k,m) = 0.0d0
924 tmr_racs1(i,j,k,m) = 0.0d0
925 tcs_racs2(i,j,k,m) = 0.0d0
926 tmr_racs2(i,j,k,m) = 0.0d0
927 tcr_sacr1(i,j,k,m) = 0.0d0
928 tms_sacr1(i,j,k,m) = 0.0d0
929 tcr_sacr2(i,j,k,m) = 0.0d0
930 tms_sacr2(i,j,k,m) = 0.0d0
931 tnr_racs1(i,j,k,m) = 0.0d0
932 tnr_racs2(i,j,k,m) = 0.0d0
933 tnr_sacr1(i,j,k,m) = 0.0d0
934 tnr_sacr2(i,j,k,m) = 0.0d0
944 tpi_qrfz(i,j,k,m) = 0.0d0
945 tni_qrfz(i,j,k,m) = 0.0d0
946 tpg_qrfz(i,j,k,m) = 0.0d0
947 tnr_qrfz(i,j,k,m) = 0.0d0
952 tpi_qcfz(i,j,k,m) = 0.0d0
953 tni_qcfz(i,j,k,m) = 0.0d0
961 tps_iaus(i,j) = 0.0d0
962 tni_iaus(i,j) = 0.0d0
979 tnr_rev(i,j,k) = 0.0d0
987 tpc_wev(i,j,k) = 0.0d0
988 tnc_wev(i,j,k) = 0.0d0
998 tnccn_act(i,j,k,l,m) = 1.0
1005 CALL wrf_debug(150, 'CREATING MICROPHYSICS LOOKUP TABLES ... ')
1006 WRITE (wrf_err_message, '(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') &
1007 ' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g
1008 CALL wrf_debug(150, wrf_err_message)
1010 !..Read a static file containing CCN activation of aerosols. The
1011 !.. data were created from a parcel model by Feingold & Heymsfield with
1012 !.. further changes by Eidhammer and Kriedenweis.
1013 if (is_aerosol_aware) then
1014 CALL wrf_debug(200, ' calling table_ccnAct routine')
1018 !..Collision efficiency between rain/snow and cloud water.
1019 CALL wrf_debug(200, ' creating qc collision eff tables')
1023 !..Drop evaporation.
1024 CALL wrf_debug(200, ' creating rain evap table')
1027 !..Initialize various constants for computing radar reflectivity.
1034 xam_g = am_g(idx_bg1)
1039 if (.not. iiwarm) then
1041 !..Rain collecting graupel & graupel collecting rain.
1042 CALL wrf_debug(200, ' creating rain collecting graupel table')
1044 call qr_acr_qg(dimNRHG)
1046 !..Rain collecting snow & snow collecting rain.
1047 CALL wrf_debug(200, ' creating rain collecting snow table')
1050 !..Cloud water and rain freezing (Bigg, 1953).
1051 CALL wrf_debug(200, ' creating freezing of water drops table')
1054 !..Conversion of some ice mass into snow category.
1055 CALL wrf_debug(200, ' creating ice converting to snow table')
1060 CALL wrf_debug(150, ' ... DONE microphysical lookup tables')
1064 END SUBROUTINE thompson_init
1065 !+---+-----------------------------------------------------------------+
1067 !+---+-----------------------------------------------------------------+
1068 !..This is a wrapper routine designed to transfer values from 3D to 1D.
1069 !+---+-----------------------------------------------------------------+
1070 SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, qb, ni, nr, nc, ng,&
1071 nwfa, nifa, nbca, nwfa2d, nifa2d, nbca2d, &
1072 aer_init_opt, wif_input_opt, &
1073 th, pii, p, w, dz, &
1077 GRAUPELNC, GRAUPELNCV, SR, &
1078 #if ( WRF_CHEM == 1 )
1079 wetscav_on, rainprod, evapprod, &
1081 refl_10cm, diagflag, ke_diag, do_radar_ref, &
1082 re_cloud, re_ice, re_snow, &
1083 has_reqc, has_reqi, has_reqs, &
1084 ids,ide, jds,jde, kds,kde, & ! domain dims
1085 ims,ime, jms,jme, kms,kme, & ! memory dims
1086 its,ite, jts,jte, kts,kte) ! tile dims
1090 !..Subroutine arguments
1091 INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, &
1092 ims,ime, jms,jme, kms,kme, &
1093 its,ite, jts,jte, kts,kte
1094 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
1095 qv, qc, qr, qi, qs, qg, ni, nr, th
1096 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: &
1097 nc, nwfa, nifa, nbca, qb, ng
1098 REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d, &
1100 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
1101 re_cloud, re_ice, re_snow
1102 INTEGER, INTENT(IN):: has_reqc, has_reqi, has_reqs
1103 #if ( WRF_CHEM == 1 )
1104 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
1107 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
1109 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: &
1111 REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(INOUT):: &
1112 SNOWNC, SNOWNCV, GRAUPELNC, GRAUPELNCV
1113 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
1115 REAL, INTENT(IN):: dt_in
1116 INTEGER, INTENT(IN):: itimestep
1117 INTEGER, OPTIONAL, INTENT(IN):: aer_init_opt, wif_input_opt
1118 #if ( WRF_CHEM == 1 )
1119 LOGICAL, INTENT(in) :: wetscav_on
1123 REAL, DIMENSION(kts:kte):: &
1124 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, &
1125 ni1d, nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d,&
1126 t1d, p1d, w1d, dz1d, rho, dBZ
1127 REAL, DIMENSION(kts:kte):: re_qc1d, re_qi1d, re_qs1d
1128 #if ( WRF_CHEM == 1 )
1129 REAL, DIMENSION(kts:kte):: &
1130 rainprod1d, evapprod1d
1132 REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
1133 REAL:: dt, pptrain, pptsnow, pptgraul, pptice
1134 REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max
1137 DOUBLE PRECISION:: lamg, lam_exp, lamr, N0_min, N0_exp
1138 INTEGER:: i, j, k, k_0, k_inj
1139 INTEGER:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr
1140 INTEGER:: jmax_qc,jmax_qr,jmax_qi,jmax_qs,jmax_qg,jmax_ni,jmax_nr
1141 INTEGER:: kmax_qc,kmax_qr,kmax_qi,kmax_qs,kmax_qg,kmax_ni,kmax_nr
1142 INTEGER:: i_start, j_start, i_end, j_end
1143 LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
1144 INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref, ke_diag
1145 CHARACTER*256:: mp_debug
1147 integer :: kediagloc
1153 i_end = MIN(ite, ide-1)
1154 j_end = MIN(jte, jde-1)
1156 !..For idealized testing by developer.
1157 ! if ( (ide-ids+1).gt.4 .and. (jde-jds+1).lt.4 .and. &
1158 ! ids.eq.its.and.ide.eq.ite.and.jds.eq.jts.and.jde.eq.jte) then
1196 mp_debug(i:i) = char(0)
1199 if (.NOT. is_aerosol_aware .AND. PRESENT(nc) .AND. PRESENT(nwfa) &
1200 .AND. PRESENT(nifa) .AND. PRESENT(nwfa2d)) then
1201 write(mp_debug,*) 'WARNING, nc-nwfa-nifa-nwfa2d present but is_aerosol_aware is FALSE'
1202 CALL wrf_debug(0, mp_debug)
1205 j_loop: do j = j_start, j_end
1206 i_loop: do i = i_start, i_end
1213 IF ( PRESENT (snowncv) ) THEN
1216 IF ( PRESENT (graupelncv) ) THEN
1217 GRAUPELNCV(i,j) = 0.
1222 t1d(k) = th(i,k,j)*pii(i,k,j)
1234 rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622))
1236 if (is_aerosol_aware) then
1239 nwfa1d(k) = nwfa(i,k,j)
1240 nifa1d(k) = nifa(i,k,j)
1241 if (wif_input_opt .eq. 2) then
1242 nbca1d(k) = nbca(i,k,j)
1250 nc1d(k) = Nt_c/rho(k)
1251 nwfa1d(k) = 11.1E6/rho(k)
1252 nifa1d(k) = naIN1*0.01/rho(k)
1253 nbca1d(k) = 5.55E6/rho(k)
1257 !..If not the variable-density graupel-hail hybrid, then set the vol mixing
1258 !.. ratio to mass mixing ratio divided by constant density (500kg/m3) value.
1260 if (is_hail_aware) then
1268 if (qg1d(k).gt.R1) then
1269 ygra1 = alog10(max(1.E-9, qg1d(k)*rho(k)))
1270 zans1 = 3.0 + 2./7.*(ygra1+8.)
1271 zans1 = MAX(2., MIN(zans1, 6.))
1272 N0_exp = 10.**(zans1)
1273 lam_exp = (N0_exp*am_g(idx_bg1)*cgg(1,1)/(rho(k)*qg1d(k)))**oge1
1274 lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg
1275 ng1d(k) = cgg(2,1)*ogg3*rho(k)*qg1d(k)*lamg**bm_g / am_g(idx_bg1)
1276 ng1d(k) = MAX(R2, ng1d(k)/rho(k))
1277 qb1d(k) = qg1d(k)/rho_g(idx_bg1)
1285 call mp_thompson(aer_init_opt, wif_input_opt, &
1286 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, ni1d, &
1287 nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d, t1d, p1d, w1d, dz1d, &
1288 pptrain, pptsnow, pptgraul, pptice, &
1289 #if ( WRF_CHEM == 1 )
1290 wetscav_on, rainprod1d, evapprod1d, &
1294 pcp_ra(i,j) = pptrain
1295 pcp_sn(i,j) = pptsnow
1296 pcp_gr(i,j) = pptgraul
1297 pcp_ic(i,j) = pptice
1298 RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice
1299 RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
1300 IF ( PRESENT(snowncv) .AND. PRESENT(snownc) ) THEN
1301 SNOWNCV(i,j) = pptsnow + pptice
1302 SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + pptice
1304 IF ( PRESENT(graupelncv) .AND. PRESENT(graupelnc) ) THEN
1305 GRAUPELNCV(i,j) = pptgraul
1306 GRAUPELNC(i,j) = GRAUPELNC(i,j) + pptgraul
1308 SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12)
1312 !..BEGIN AEROSOL EMISSIONS
1314 !..Reset lowest model level to initial state aerosols (fake sfc source).
1315 !.. Changed 13 May 2013 to fake emissions in which nwfa2d is aerosol
1316 !.. number tendency (number per kg per second).
1317 if (is_aerosol_aware) then
1318 !..Add anthropogenic emissions
1319 !-GT nwfa1d(kts) = nwfa1
1320 nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt_in
1321 nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt_in
1322 if (wif_input_opt .eq. 2) then
1323 nbca1d(kts) = nbca1d(kts) + nbca2d(i,j)*dt_in
1328 !..END AEROSOL EMISSIONS
1333 nwfa(i,k,j) = nwfa1d(k)
1334 nifa(i,k,j) = nifa1d(k)
1335 if (wif_input_opt .eq. 2) then
1336 nbca(i,k,j) = nbca1d(k)
1342 if (is_hail_aware) then
1358 th(i,k,j) = t1d(k)/pii(i,k,j)
1359 #if ( WRF_CHEM == 1 )
1360 IF ( wetscav_on ) THEN
1361 rainprod(i,k,j) = rainprod1d(k)
1362 evapprod(i,k,j) = evapprod1d(k)
1365 if (qc1d(k) .gt. qc_max) then
1370 elseif (qc1d(k) .lt. 0.0) then
1371 write(mp_debug,*) 'WARNING, negative qc ', qc1d(k), &
1373 CALL wrf_debug(150, mp_debug)
1375 if (qr1d(k) .gt. qr_max) then
1380 elseif (qr1d(k) .lt. 0.0) then
1381 write(mp_debug,*) 'WARNING, negative qr ', qr1d(k), &
1383 CALL wrf_debug(150, mp_debug)
1385 if (nr1d(k) .gt. nr_max) then
1390 elseif (nr1d(k) .lt. 0.0) then
1391 write(mp_debug,*) 'WARNING, negative nr ', nr1d(k), &
1393 CALL wrf_debug(150, mp_debug)
1395 if (qs1d(k) .gt. qs_max) then
1400 elseif (qs1d(k) .lt. 0.0) then
1401 write(mp_debug,*) 'WARNING, negative qs ', qs1d(k), &
1403 CALL wrf_debug(150, mp_debug)
1405 if (qi1d(k) .gt. qi_max) then
1410 elseif (qi1d(k) .lt. 0.0) then
1411 write(mp_debug,*) 'WARNING, negative qi ', qi1d(k), &
1413 CALL wrf_debug(150, mp_debug)
1415 if (qg1d(k) .gt. qg_max) then
1420 elseif (qg1d(k) .lt. 0.0) then
1421 write(mp_debug,*) 'WARNING, negative qg ', qg1d(k), &
1423 CALL wrf_debug(150, mp_debug)
1425 if (ni1d(k) .gt. ni_max) then
1430 elseif (ni1d(k) .lt. 0.0) then
1431 write(mp_debug,*) 'WARNING, negative ni ', ni1d(k), &
1433 CALL wrf_debug(150, mp_debug)
1435 if (qv1d(k) .lt. 0.0) then
1436 write(mp_debug,*) 'WARNING, negative qv ', qv1d(k), &
1438 CALL wrf_debug(150, mp_debug)
1439 if (k.lt.kte-2 .and. k.gt.kts+1) then
1440 write(mp_debug,*) ' below and above are: ', qv(i,k-1,j), qv(i,k+1,j)
1441 CALL wrf_debug(150, mp_debug)
1442 qv(i,k,j) = MAX(1.E-7, 0.5*(qv(i,k-1,j) + qv(i,k+1,j)))
1449 IF ( PRESENT (diagflag) ) THEN
1450 if (diagflag .and. do_radar_ref == 1) then
1452 IF ( present(ke_diag) ) THEN
1458 call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, &
1459 ng1d, qb1d, t1d, p1d, dBZ, kts, kte, i, j, kediagloc)
1461 refl_10cm(i,k,j) = MAX(-35., dBZ(k))
1466 IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN
1468 re_qc1d(k) = RE_QC_BG
1469 re_qi1d(k) = RE_QI_BG
1470 re_qs1d(k) = RE_QS_BG
1472 call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
1473 re_qc1d, re_qi1d, re_qs1d, kts, kte)
1475 re_cloud(i,k,j) = MAX(RE_QC_BG, MIN(re_qc1d(k), 50.E-6))
1476 re_ice(i,k,j) = MAX(RE_QI_BG, MIN(re_qi1d(k), 125.E-6))
1477 re_snow(i,k,j) = MAX(RE_QS_BG, MIN(re_qs1d(k), 999.E-6))
1485 write(mp_debug,'(a,7(a,e13.6,1x,a,i3,a,i3,a,i3,a,1x))') 'MP-GT:', &
1486 'qc: ', qc_max, '(', imax_qc, ',', jmax_qc, ',', kmax_qc, ')', &
1487 'qr: ', qr_max, '(', imax_qr, ',', jmax_qr, ',', kmax_qr, ')', &
1488 'qi: ', qi_max, '(', imax_qi, ',', jmax_qi, ',', kmax_qi, ')', &
1489 'qs: ', qs_max, '(', imax_qs, ',', jmax_qs, ',', kmax_qs, ')', &
1490 'qg: ', qg_max, '(', imax_qg, ',', jmax_qg, ',', kmax_qg, ')', &
1491 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')', &
1492 'nr: ', nr_max, '(', imax_nr, ',', jmax_nr, ',', kmax_nr, ')'
1493 CALL wrf_debug(150, mp_debug)
1497 mp_debug(i:i) = char(0)
1500 END SUBROUTINE mp_gt_driver
1502 !+---+-----------------------------------------------------------------+
1504 !+---+-----------------------------------------------------------------+
1505 !+---+-----------------------------------------------------------------+
1506 !.. This subroutine computes the moisture tendencies of water vapor,
1507 !.. cloud droplets, rain, cloud ice (pristine), snow, and graupel.
1508 !.. Previously this code was based on Reisner et al (1998), but few of
1509 !.. those pieces remain. A complete description is now found in
1510 !.. Thompson et al. (2004, 2008).
1511 !+---+-----------------------------------------------------------------+
1513 subroutine mp_thompson (aer_init_opt, wif_input_opt, &
1514 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, &
1515 ni1d, nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d, &
1516 t1d, p1d, w1d, dzq, &
1517 pptrain, pptsnow, pptgraul, pptice, &
1518 #if ( WRF_CHEM == 1 )
1519 wetscav_on, rainprod, evapprod, &
1521 kts, kte, dt, ii, jj)
1526 INTEGER, OPTIONAL, INTENT(IN):: aer_init_opt, wif_input_opt
1527 INTEGER, INTENT(IN):: kts, kte, ii, jj
1528 REAL, DIMENSION(kts:kte), INTENT(INOUT):: &
1529 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, &
1530 ni1d, nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d, t1d
1531 REAL, DIMENSION(kts:kte), INTENT(IN):: p1d, w1d, dzq
1532 REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice
1533 REAL, INTENT(IN):: dt
1534 #if ( WRF_CHEM == 1 )
1535 REAL, DIMENSION(kts:kte), INTENT(INOUT):: &
1537 LOGICAL, INTENT(IN) :: wetscav_on
1541 REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, qrten, &
1542 qsten, qgten, qbten, niten, nrten, ncten, ngten, nwfaten, nifaten, &
1545 DOUBLE PRECISION, DIMENSION(kts:kte):: prw_vcd
1547 DOUBLE PRECISION, DIMENSION(kts:kte):: pnc_wcd, pnc_wau, pnc_rcw, &
1550 DOUBLE PRECISION, DIMENSION(kts:kte):: pna_rca, pna_sca, pna_gca, &
1551 pnd_rcd, pnd_scd, pnd_gcd, pnb_rcb, pnb_scb, pnb_gcb
1553 DOUBLE PRECISION, DIMENSION(kts:kte):: prr_wau, prr_rcw, prr_rcs, &
1554 prr_rcg, prr_sml, prr_gml, &
1556 pnr_wau, pnr_rcs, pnr_rcg, &
1557 pnr_rci, pnr_sml, pnr_gml, &
1558 pnr_rev, pnr_rcr, pnr_rfz
1560 DOUBLE PRECISION, DIMENSION(kts:kte):: pri_inu, pni_inu, pri_ihm, &
1561 pni_ihm, pri_wfz, pni_wfz, &
1562 pri_rfz, pni_rfz, pri_ide, &
1563 pni_ide, pri_rci, pni_rci, &
1564 pni_sci, pni_iau, pri_iha, pni_iha
1566 DOUBLE PRECISION, DIMENSION(kts:kte):: prs_iau, prs_sci, prs_rcs, &
1567 prs_scw, prs_sde, prs_ihm, &
1570 DOUBLE PRECISION, DIMENSION(kts:kte):: prg_scw, prg_rfz, prg_gde, &
1571 prg_gcw, prg_rci, prg_rcs, prg_rcg, prg_ihm, &
1572 png_rcs, png_rcg, png_scw, png_gde, &
1573 pbg_scw, pbg_rfz, pbg_gcw, pbg_rci, pbg_rcs, pbg_rcg, &
1576 DOUBLE PRECISION, PARAMETER:: zeroD0 = 0.0d0
1578 REAL, DIMENSION(kts:kte):: temp, twet, pres, qv
1579 REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, rb
1580 REAL, DIMENSION(kts:kte):: ni, nr, nc, ns, ng, nwfa, nifa, nbca
1581 REAL, DIMENSION(kts:kte):: rho, rhof, rhof2
1582 REAL, DIMENSION(kts:kte):: qvs, qvsi, delQvs
1583 REAL, DIMENSION(kts:kte):: satw, sati, ssatw, ssati
1584 REAL, DIMENSION(kts:kte):: diffu, visco, vsc2, &
1585 tcond, lvap, ocp, lvt2
1587 DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g
1588 DOUBLE PRECISION:: N0_melt
1589 REAL, DIMENSION(kts:kte):: mvd_r, mvd_c, mvd_g
1590 REAL, DIMENSION(kts:kte):: smob, smo2, smo1, smo0, &
1591 smoc, smod, smoe, smof, smog
1593 REAL, DIMENSION(kts:kte):: sed_r,sed_s,sed_g,sed_i,sed_n,sed_c,sed_b
1595 REAL:: rgvm, delta_tp, orho, lfus2
1596 REAL, DIMENSION(5):: onstep
1597 DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamc, lamr, lamg
1598 DOUBLE PRECISION:: lami, ilami, ilamc
1599 REAL:: xDc, Dc_b, Dc_g, xDi, xDr, xDs, xDg, Ds_m, Dg_m
1600 DOUBLE PRECISION:: Dr_star, Dc_star
1601 REAL:: zeta1, zeta, taud, tau
1602 REAL:: stoke_r, stoke_s, stoke_g, stoke_i
1603 REAL:: vti, vtr, vts, vtg, vtc
1604 REAL:: xrho_g, afall, vtg1, vtg2
1605 REAL:: bfall = 3*b_coeff - 1
1606 REAL, DIMENSION(kts:kte+1):: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk, &
1608 REAL, DIMENSION(kts:kte):: vts_boost
1609 REAL:: M0, slam1, slam2
1610 REAL:: Mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, C_snow
1611 REAL:: a_, b_, loga_, A1, A2, tf
1612 REAL:: tempc, tc0, r_mvd1, r_mvd2, xkrat
1613 REAL:: dew_t, Tlcl, The
1614 REAL:: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr, xrg, xng, xrb
1615 REAL:: xsat, rate_max, sump, ratio
1616 REAL:: clap, fcd, dfcd
1617 REAL:: otemp, rvs, rvs_p, rvs_pp, gamsc, alphsc, t1_evap, t1_subl
1618 REAL:: r_frac, g_frac, const_Ri, rime_dens
1619 REAL:: Ef_rw, Ef_sw, Ef_gw, Ef_rr
1620 REAL:: Ef_ra, Ef_sa, Ef_ga
1621 REAL:: dtsave, odts, odt, odzq, hgt_agl
1622 REAL:: xslw1, ygra1, zans1, eva_factor
1624 INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq, k_melting
1625 INTEGER, DIMENSION(5):: ksed1
1626 INTEGER:: nir, nis, nig, nii, nic, niin
1627 INTEGER:: idx_tc, idx_t, idx_s, idx_g1, idx_g, idx_r1, idx_r, &
1628 idx_i1, idx_i, idx_c, idx, idx_d, idx_n, idx_in
1629 INTEGER, DIMENSION(kts:kte):: idx_bg
1631 LOGICAL:: melti, no_micro
1632 LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg
1633 LOGICAL:: debug_flag
1634 CHARACTER*256:: mp_debug
1639 debug_flag = .false.
1640 ! if (ii.eq.901 .and. jj.eq.379) debug_flag = .true.
1642 write(mp_debug, *) 'DEBUG INFO, mp_thompson at (i,j) ', ii, ', ', jj
1643 CALL wrf_debug(550, mp_debug)
1652 !+---+-----------------------------------------------------------------+
1653 !.. Source/sink terms. First 2 chars: "pr" represents source/sink of
1654 !.. mass while "pn" represents source/sink of number. Next char is one
1655 !.. of "v" for water vapor, "r" for rain, "i" for cloud ice, "w" for
1656 !.. cloud water, "s" for snow, and "g" for graupel. Next chars
1657 !.. represent processes: "de" for sublimation/deposition, "ev" for
1658 !.. evaporation, "fz" for freezing, "ml" for melting, "au" for
1659 !.. autoconversion, "nu" for ice nucleation, "hm" for Hallet/Mossop
1660 !.. secondary ice production, and "c" for collection followed by the
1661 !.. character for the species being collected. ALL of these terms are
1662 !.. positive (except for deposition/sublimation terms which can switch
1663 !.. signs based on super/subsaturation) and are treated as negatives
1664 !.. where necessary in the tendency equations.
1665 !+---+-----------------------------------------------------------------+
1743 ! new source/sink terms for 3-moment graupel
1770 #if ( WRF_CHEM == 1 )
1771 if ( wetscav_on ) then
1779 !..Bug fix (2016Jun15), prevent use of uninitialized value(s) of snow moments.
1795 !+---+-----------------------------------------------------------------+
1796 !..Put column of data into local arrays.
1797 !+---+-----------------------------------------------------------------+
1800 qv(k) = MAX(1.E-10, qv1d(k))
1802 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
1803 if (is_aerosol_aware) then
1804 if (aer_init_opt .lt. 2) then ! Constant or climatology (e.g., GOCART)
1805 nwfa(k) = MAX(11.1E6, MIN(9999.E6, nwfa1d(k)*rho(k)))
1806 nifa(k) = MAX(naIN1*0.01, MIN(9999.E6, nifa1d(k)*rho(k)))
1807 if (wif_input_opt .eq. 2) then ! Considering BC aerosol
1808 nbca(k) = MAX(5.55E6, MIN(9999.E6, nbca1d(k)*rho(k)))
1812 else ! First guess (e.g., GEOS-5)
1813 nwfa(k) = MAX(0.0, nwfa1d(k)*rho(k))
1814 nifa(k) = MAX(0.0, nifa1d(k)*rho(k))
1815 if (wif_input_opt .eq. 2) then ! Considering BC aerosol
1816 nbca(k) = MAX(0.0, nbca1d(k)*rho(k))
1822 nwfa(k) = MAX(11.1E6, MIN(9999.E6, nwfa1d(k)*rho(k)))
1823 nifa(k) = MAX(naIN1*0.01, MIN(9999.E6, nifa1d(k)*rho(k)))
1824 nbca(k) = MAX(5.55E6, MIN(9999.E6, nbca1d(k)*rho(k)))
1827 if (qc1d(k) .gt. R1) then
1829 rc(k) = qc1d(k)*rho(k)
1830 nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max))
1832 nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2)
1833 lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
1834 xDc = (bm_r + nu_c + 1.) / lamc
1835 if (xDc.lt. D0c) then
1836 lamc = cce(2,nu_c)/D0c
1837 elseif (xDc.gt. D0r*2.) then
1838 lamc = cce(2,nu_c)/(D0r*2.)
1840 nc(k) = MIN( DBLE(Nt_c_max), ccg(1,nu_c)*ocg2(nu_c)*rc(k) &
1842 if (.NOT. is_aerosol_aware) nc(k) = Nt_c
1851 if (qi1d(k) .gt. R1) then
1853 ri(k) = qi1d(k)*rho(k)
1854 ni(k) = MAX(R2, ni1d(k)*rho(k))
1855 if (ni(k).le. R2) then
1857 ni(k) = MIN(999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
1860 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
1862 xDi = (bm_i + mu_i + 1.) * ilami
1863 if (xDi.lt. 5.E-6) then
1865 ni(k) = MIN(999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
1866 elseif (xDi.gt. 300.E-6) then
1867 lami = cie(2)/300.E-6
1868 ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i
1878 if (qr1d(k) .gt. R1) then
1880 rr(k) = qr1d(k)*rho(k)
1881 nr(k) = MAX(R2, nr1d(k)*rho(k))
1882 if (nr(k).le. R2) then
1884 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
1885 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
1888 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
1889 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
1890 if (mvd_r(k) .gt. 2.5E-3) then
1892 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
1893 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
1894 elseif (mvd_r(k) .lt. D0r*0.75) then
1896 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
1897 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
1906 if (qs1d(k) .gt. R1) then
1908 rs(k) = qs1d(k)*rho(k)
1915 if (qg1d(k) .gt. R1) then
1918 rg(k) = qg1d(k)*rho(k)
1919 ng(k) = MAX(R2, ng1d(k)*rho(k))
1920 rb(k) = MAX(qg1d(k)/rho_g(NRHG), qb1d(k))
1921 rb(k) = MIN(qg1d(k)/rho_g(1), rb(k))
1923 idx_bg(k) = MAX(1,MIN(NINT(qg1d(k)/rb(k) *0.01)+1,NRHG))
1924 if (ng(k).le. R2) then
1926 lamg = (3.0 + mu_g + 0.672) / mvd_g(k)
1927 ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k))
1929 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg
1930 mvd_g(k) = (3.0 + mu_g + 0.672) / lamg
1931 if (mvd_g(k) .gt. 25.4E-3) then
1933 lamg = (3.0 + mu_g + 0.672) / mvd_g(k)
1934 ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k))
1935 elseif (mvd_g(k) .lt. D0r) then
1937 lamg = (3.0 + mu_g + 0.672) / mvd_g(k)
1938 ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k))
1947 rb(k) = R1/rho(k)/rho_g(NRHG)
1950 if (.not. is_hail_aware) idx_bg(k) = idx_bg1
1953 !+---+-----------------------------------------------------------------+
1954 ! if (debug_flag) then
1955 ! write(mp_debug,*) 'DEBUG-VERBOSE at (i,j) ', ii, ', ', jj
1956 ! CALL wrf_debug(550, mp_debug)
1958 ! write(mp_debug, '(a,i3,f8.2,1x,f7.2,1x, 13(1x,e13.6))') &
1959 ! & 'VERBOSE: ', k, pres(k)*0.01, temp(k)-273.15, qv(k), rc(k), rr(k), ri(k), rs(k), rg(k), nc(k), nr(k), ni(k), ng(k), rb(k), nwfa(k), nifa(k)
1960 ! CALL wrf_debug(550, mp_debug)
1963 !+---+-----------------------------------------------------------------+
1965 !+---+-----------------------------------------------------------------+
1966 !..Derive various thermodynamic variables frequently used.
1967 !.. Saturation vapor pressure (mixing ratio) over liquid/ice comes from
1968 !.. Flatau et al. 1992; enthalpy (latent heat) of vaporization from
1969 !.. Bohren & Albrecht 1998; others from Pruppacher & Klett 1978.
1970 !+---+-----------------------------------------------------------------+
1973 tempc = temp(k) - 273.15
1974 rhof(k) = SQRT(RHO_NOT/rho(k))
1975 rhof2(k) = SQRT(rhof(k))
1976 qvs(k) = rslf(pres(k), temp(k))
1977 delQvs(k) = MAX(0.0, rslf(pres(k), 273.15)-qv(k))
1978 if (tempc .le. 0.0) then
1979 qvsi(k) = rsif(pres(k), temp(k))
1982 k_melting = MAX(k, k_melting)
1984 satw(k) = qv(k)/qvs(k)
1985 sati(k) = qv(k)/qvsi(k)
1986 ssatw(k) = satw(k) - 1.
1987 ssati(k) = sati(k) - 1.
1988 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
1989 if (abs(ssati(k)).lt. eps) ssati(k) = 0.0
1990 if (no_micro .and. ssati(k).gt. 0.0) no_micro = .false.
1991 diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
1992 if (tempc .ge. 0.0) then
1993 visco(k) = (1.718+0.0049*tempc)*1.0E-5
1995 visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
1997 ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
1998 vsc2(k) = SQRT(rho(k)/visco(k))
1999 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
2000 tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
2004 if (k_melting .gt. 0) then
2005 do k = kts, k_melting
2006 if (satw(k) .lt. 0.999) then
2007 dew_t = MIN(temp(k)-0.001, t_dew(pres(k), qv(k)))
2008 Tlcl = t_lcl(temp(k), dew_t)
2009 The = theta_e(pres(k), temp(k), qv(k), Tlcl)
2010 twet(k) = MIN(temp(k), compT_fr_The(The, pres(k)))
2015 !+---+-----------------------------------------------------------------+
2016 !..If no existing hydrometeor species and no chance to initiate ice or
2017 !.. condense cloud water, just exit quickly!
2018 !+---+-----------------------------------------------------------------+
2020 if (no_micro) return
2022 !+---+-----------------------------------------------------------------+
2023 !..Calculate y-intercept, slope, and useful moments for snow.
2024 !+---+-----------------------------------------------------------------+
2025 if (.not. iiwarm) then
2027 if (.not. L_qs(k)) CYCLE
2028 tc0 = MIN(-0.1, temp(k)-273.15)
2029 smob(k) = rs(k)*oams
2031 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
2032 !.. then we must compute actual 2nd moment and use as reference.
2033 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
2036 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
2037 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
2038 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
2039 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
2040 + sa(10)*bm_s*bm_s*bm_s
2042 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
2043 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
2044 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
2045 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
2046 + sb(10)*bm_s*bm_s*bm_s
2047 smo2(k) = (smob(k)/a_)**(1./b_)
2050 !..Calculate 0th moment. Represents snow number concentration.
2051 loga_ = sa(1) + sa(2)*tc0 + sa(5)*tc0*tc0 + sa(9)*tc0*tc0*tc0
2053 b_ = sb(1) + sb(2)*tc0 + sb(5)*tc0*tc0 + sb(9)*tc0*tc0*tc0
2054 smo0(k) = a_ * smo2(k)**b_
2056 !..Calculate 1st moment. Useful for depositional growth and melting.
2057 loga_ = sa(1) + sa(2)*tc0 + sa(3) &
2058 + sa(4)*tc0 + sa(5)*tc0*tc0 &
2059 + sa(6) + sa(7)*tc0*tc0 &
2060 + sa(8)*tc0 + sa(9)*tc0*tc0*tc0 &
2063 b_ = sb(1)+ sb(2)*tc0 + sb(3) + sb(4)*tc0 &
2064 + sb(5)*tc0*tc0 + sb(6) &
2065 + sb(7)*tc0*tc0 + sb(8)*tc0 &
2066 + sb(9)*tc0*tc0*tc0 + sb(10)
2067 smo1(k) = a_ * smo2(k)**b_
2069 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
2070 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
2071 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
2072 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
2073 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
2074 + sa(10)*cse(1)*cse(1)*cse(1)
2076 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
2077 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
2078 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
2079 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
2080 smoc(k) = a_ * smo2(k)**b_
2082 !..Calculate snow number concentration (explicit integral, not smo0)
2083 M0 = smob(k)/smoc(k)
2084 Mrat = smob(k)*M0*M0*M0
2087 ns(k) = Mrat*Kap0/slam1 &
2088 + Mrat*Kap1*M0**mu_s*csg(15)/slam2**cse(15)
2090 !..Calculate bv_s+2 (th) moment. Useful for riming.
2091 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(13) &
2092 + sa(4)*tc0*cse(13) + sa(5)*tc0*tc0 &
2093 + sa(6)*cse(13)*cse(13) + sa(7)*tc0*tc0*cse(13) &
2094 + sa(8)*tc0*cse(13)*cse(13) + sa(9)*tc0*tc0*tc0 &
2095 + sa(10)*cse(13)*cse(13)*cse(13)
2097 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(13) + sb(4)*tc0*cse(13) &
2098 + sb(5)*tc0*tc0 + sb(6)*cse(13)*cse(13) &
2099 + sb(7)*tc0*tc0*cse(13) + sb(8)*tc0*cse(13)*cse(13) &
2100 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(13)*cse(13)*cse(13)
2101 smoe(k) = a_ * smo2(k)**b_
2103 !..Calculate 1+(bv_s+1)/2 (th) moment. Useful for depositional growth.
2104 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(16) &
2105 + sa(4)*tc0*cse(16) + sa(5)*tc0*tc0 &
2106 + sa(6)*cse(16)*cse(16) + sa(7)*tc0*tc0*cse(16) &
2107 + sa(8)*tc0*cse(16)*cse(16) + sa(9)*tc0*tc0*tc0 &
2108 + sa(10)*cse(16)*cse(16)*cse(16)
2110 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(16) + sb(4)*tc0*cse(16) &
2111 + sb(5)*tc0*tc0 + sb(6)*cse(16)*cse(16) &
2112 + sb(7)*tc0*tc0*cse(16) + sb(8)*tc0*cse(16)*cse(16) &
2113 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(16)*cse(16)*cse(16)
2114 smof(k) = a_ * smo2(k)**b_
2116 !..Calculate bm_s + bv_s+2 (th) moment. Useful for riming into graupel.
2117 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(17) &
2118 + sa(4)*tc0*cse(17) + sa(5)*tc0*tc0 &
2119 + sa(6)*cse(17)*cse(17) + sa(7)*tc0*tc0*cse(17) &
2120 + sa(8)*tc0*cse(17)*cse(17) + sa(9)*tc0*tc0*tc0 &
2121 + sa(10)*cse(17)*cse(17)*cse(17)
2123 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(17) + sb(4)*tc0*cse(17) &
2124 + sb(5)*tc0*tc0 + sb(6)*cse(17)*cse(17) &
2125 + sb(7)*tc0*tc0*cse(17) + sb(8)*tc0*cse(17)*cse(17) &
2126 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(17)*cse(17)*cse(17)
2127 smog(k) = a_ * smo2(k)**b_
2131 !+---+-----------------------------------------------------------------+
2132 !..Calculate y-intercept, slope values for graupel.
2133 !+---+-----------------------------------------------------------------+
2136 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg
2138 N0_g(k) = ng(k)*ogg2*lamg**cge(2,1)
2143 !+---+-----------------------------------------------------------------+
2144 !..Calculate y-intercept, slope values for rain.
2145 !+---+-----------------------------------------------------------------+
2147 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
2149 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
2150 N0_r(k) = nr(k)*org2*lamr**cre(2)
2153 !+---+-----------------------------------------------------------------+
2154 !..Compute warm-rain process terms (except evap done later).
2155 !+---+-----------------------------------------------------------------+
2159 !..Rain self-collection follows Seifert, 1994 and drop break-up
2160 !.. follows Verlinde and Cotton, 1993. RAIN2M
2161 if (L_qr(k) .and. mvd_r(k).gt. D0r) then
2163 !-GT if (mvd_r(k) .gt. 1500.0E-6) then
2164 Ef_rr = 1.0 - EXP(2300.0*(mvd_r(k)-1950.0E-6))
2166 pnr_rcr(k) = Ef_rr * 2.0*nr(k)*rr(k)
2171 nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2)
2172 xDc = MAX(D0c*1.E6, ((rc(k)/(am_r*nc(k)))**obmr) * 1.E6)
2173 lamc = (nc(k)*am_r* ccg(2,nu_c) * ocg1(nu_c) / rc(k))**obmr
2174 mvd_c(k) = (3.0+nu_c+0.672) / lamc
2175 mvd_c(k) = MAX(D0c, MIN(mvd_c(k), D0r))
2178 !..Autoconversion follows Berry & Reinhardt (1974) with characteristic
2179 !.. diameters correctly computed from gamma distrib of cloud droplets.
2180 if (rc(k).gt. 0.01e-3) then
2181 Dc_g = ((ccg(3,nu_c)*ocg2(nu_c))**obmr / lamc) * 1.E6
2182 Dc_b = (xDc*xDc*xDc*Dc_g*Dc_g*Dc_g - xDc*xDc*xDc*xDc*xDc*xDc) &
2184 zeta1 = 0.5*((6.25E-6*xDc*Dc_b*Dc_b*Dc_b - 0.4) &
2185 + abs(6.25E-6*xDc*Dc_b*Dc_b*Dc_b - 0.4))
2186 zeta = 0.027*rc(k)*zeta1
2187 taud = 0.5*((0.5*Dc_b - 7.5) + abs(0.5*Dc_b - 7.5)) + R1
2188 tau = 3.72/(rc(k)*taud)
2189 prr_wau(k) = zeta/tau
2190 prr_wau(k) = MIN(DBLE(rc(k)*odts), prr_wau(k))
2191 pnr_wau(k) = prr_wau(k) / (am_r*nu_c*10.*D0r*D0r*D0r) ! RAIN2M
2192 pnc_wau(k) = MIN(DBLE(nc(k)*odts), prr_wau(k) &
2193 / (am_r*mvd_c(k)*mvd_c(k)*mvd_c(k))) ! Qc2M
2196 !..Rain collecting cloud water. In CE, assume Dc<<Dr and vtc=~0.
2197 if (L_qr(k) .and. mvd_r(k).gt. D0r .and. mvd_c(k).gt. D0c) then
2199 idx = 1 + INT(nbr*DLOG(mvd_r(k)/Dr(1))/DLOG(Dr(nbr)/Dr(1)))
2201 Ef_rw = t_Efrw(idx, INT(mvd_c(k)*1.E6))
2202 prr_rcw(k) = rhof(k)*t1_qr_qc*Ef_rw*rc(k)*N0_r(k) &
2203 *((lamr+fv_r)**(-cre(9)))
2204 prr_rcw(k) = MIN(DBLE(rc(k)*odts), prr_rcw(k))
2205 pnc_rcw(k) = rhof(k)*t1_qr_qc*Ef_rw*nc(k)*N0_r(k) &
2206 *((lamr+fv_r)**(-cre(9))) ! Qc2M
2207 pnc_rcw(k) = MIN(DBLE(nc(k)*odts), pnc_rcw(k))
2210 !..Rain collecting aerosols, wet scavenging.
2211 if (L_qr(k) .and. mvd_r(k).gt. D0r) then
2212 Ef_ra = Eff_aero(mvd_r(k),0.04E-6,visco(k),rho(k),temp(k),'r')
2214 pna_rca(k) = rhof(k)*t1_qr_qc*Ef_ra*nwfa(k)*N0_r(k) &
2215 *((lamr+fv_r)**(-cre(9)))
2216 pna_rca(k) = MIN(DBLE(nwfa(k)*odts), pna_rca(k))
2218 Ef_ra = Eff_aero(mvd_r(k),0.8E-6,visco(k),rho(k),temp(k),'r')
2219 pnd_rcd(k) = rhof(k)*t1_qr_qc*Ef_ra*nifa(k)*N0_r(k) &
2220 *((lamr+fv_r)**(-cre(9)))
2221 pnd_rcd(k) = MIN(DBLE(nifa(k)*odts), pnd_rcd(k))
2223 if (present(wif_input_opt)) then
2224 if (wif_input_opt .eq. 2) then
2225 Ef_ra = Eff_aero(mvd_r(k),0.0236E-6, &
2226 visco(k),rho(k),temp(k),'r')
2227 pnb_rcb(k) = rhof(k)*t1_qr_qc*Ef_ra*nbca(k)*N0_r(k) &
2228 *((lamr+fv_r)**(-cre(9)))
2229 pnb_rcb(k) = MIN(DBLE(nbca(k)*odts), pnb_rcb(k))
2236 !+---+-----------------------------------------------------------------+
2237 !..Compute all frozen hydrometeor species' process terms.
2238 !+---+-----------------------------------------------------------------+
2239 if (.not. iiwarm) then
2240 !..vts_boost is the factor applied to snow terminal
2241 !..fallspeed due to riming of snow
2247 xDs = smoc(k) / smob(k)
2248 rho_s = MAX(rho_g(1), MIN(0.13/xDs, rho_i-100.))
2254 !..Temperature lookup table indexes.
2255 tempc = temp(k) - 273.15
2256 idx_tc = MAX(1, MIN(NINT(-tempc), 45) )
2257 idx_t = INT( (tempc-2.5)/5. ) - 1
2258 idx_t = MAX(1, -idx_t)
2259 idx_t = MIN(idx_t, ntb_t)
2260 IT = MAX(1, MIN(NINT(-tempc), 31) )
2262 !..Cloud water lookup table index.
2263 if (rc(k).gt. r_c(1)) then
2264 nic = NINT(ALOG10(rc(k)))
2265 do nn = nic-1, nic+1
2267 if ( (rc(k)/10.**nn).ge.1.0 .and. &
2268 (rc(k)/10.**nn).lt.10.0) goto 141
2271 idx_c = INT(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
2272 idx_c = MAX(1, MIN(idx_c, ntb_c))
2277 !..Cloud droplet number lookup table index.
2278 idx_n = NINT(1.0 + FLOAT(nbc) * DLOG(nc(k)/t_Nc(1)) / nic1)
2279 idx_n = MAX(1, MIN(idx_n, nbc))
2281 !..Cloud ice lookup table indexes.
2282 if (ri(k).gt. r_i(1)) then
2283 nii = NINT(ALOG10(ri(k)))
2284 do nn = nii-1, nii+1
2286 if ( (ri(k)/10.**nn).ge.1.0 .and. &
2287 (ri(k)/10.**nn).lt.10.0) goto 142
2290 idx_i = INT(ri(k)/10.**n) + 10*(n-nii2) - (n-nii2)
2291 idx_i = MAX(1, MIN(idx_i, ntb_i))
2296 if (ni(k).gt. Nt_i(1)) then
2297 nii = NINT(ALOG10(ni(k)))
2298 do nn = nii-1, nii+1
2300 if ( (ni(k)/10.**nn).ge.1.0 .and. &
2301 (ni(k)/10.**nn).lt.10.0) goto 143
2304 idx_i1 = INT(ni(k)/10.**n) + 10*(n-nii3) - (n-nii3)
2305 idx_i1 = MAX(1, MIN(idx_i1, ntb_i1))
2310 !..Rain lookup table indexes.
2311 if (rr(k).gt. r_r(1)) then
2312 nir = NINT(ALOG10(rr(k)))
2313 do nn = nir-1, nir+1
2315 if ( (rr(k)/10.**nn).ge.1.0 .and. &
2316 (rr(k)/10.**nn).lt.10.0) goto 144
2319 idx_r = INT(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2)
2320 idx_r = MAX(1, MIN(idx_r, ntb_r))
2323 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
2324 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
2325 nir = NINT(DLOG10(N0_exp))
2326 do nn = nir-1, nir+1
2328 if ( (N0_exp/10.**nn).ge.1.0 .and. &
2329 (N0_exp/10.**nn).lt.10.0) goto 145
2332 idx_r1 = INT(N0_exp/10.**n) + 10*(n-nir3) - (n-nir3)
2333 idx_r1 = MAX(1, MIN(idx_r1, ntb_r1))
2339 !..Snow lookup table index.
2340 if (rs(k).gt. r_s(1)) then
2341 nis = NINT(ALOG10(rs(k)))
2342 do nn = nis-1, nis+1
2344 if ( (rs(k)/10.**nn).ge.1.0 .and. &
2345 (rs(k)/10.**nn).lt.10.0) goto 146
2348 idx_s = INT(rs(k)/10.**n) + 10*(n-nis2) - (n-nis2)
2349 idx_s = MAX(1, MIN(idx_s, ntb_s))
2354 !..Graupel lookup table index.
2355 if (rg(k).gt. r_g(1)) then
2356 nig = NINT(ALOG10(rg(k)))
2357 do nn = nig-1, nig+1
2359 if ( (rg(k)/10.**nn).ge.1.0 .and. &
2360 (rg(k)/10.**nn).lt.10.0) goto 147
2363 idx_g = INT(rg(k)/10.**n) + 10*(n-nig2) - (n-nig2)
2364 idx_g = MAX(1, MIN(idx_g, ntb_g))
2367 lam_exp = lamg * (cgg(3,1)*ogg2*ogg1)**bm_g
2368 N0_exp = ogg1*rg(k)/am_g(idx_bg(k)) * lam_exp**cge(1,1)
2369 nig = NINT(DLOG10(N0_exp))
2370 do nn = nig-1, nig+1
2372 if ( (N0_exp/10.**nn).ge.1.0 .and. &
2373 (N0_exp/10.**nn).lt.10.0) goto 148
2376 idx_g1 = INT(N0_exp/10.**n) + 10*(n-nig3) - (n-nig3)
2377 idx_g1 = MAX(1, MIN(idx_g1, ntb_g1))
2383 !..Deposition/sublimation prefactor (from Srivastava & Coen 1992).
2385 rvs = rho(k)*qvsi(k)
2386 rvs_p = rvs*otemp*(lsub*otemp*oRv - 1.)
2387 rvs_pp = rvs * ( otemp*(lsub*otemp*oRv - 1.) &
2388 *otemp*(lsub*otemp*oRv - 1.) &
2389 + (-2.*lsub*otemp*otemp*otemp*oRv) &
2391 gamsc = lsub*diffu(k)/tcond(k) * rvs_p
2392 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
2393 * rvs_pp/rvs_p * rvs/rvs_p
2394 alphsc = MAX(1.E-9, alphsc)
2396 if (abs(xsat).lt. 1.E-9) xsat=0.
2397 t1_subl = 4.*PI*( 1.0 - alphsc*xsat &
2398 + 2.*alphsc*alphsc*xsat*xsat &
2399 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
2402 !..Snow collecting cloud water. In CE, assume Dc<<Ds and vtc=~0.
2403 if (L_qc(k) .and. mvd_c(k).gt. D0c) then
2404 if (xDs .gt. D0s) then
2405 idx = 1 + INT(nbs*DLOG(xDs/Ds(1))/DLOG(Ds(nbs)/Ds(1)))
2407 Ef_sw = t_Efsw(idx, INT(mvd_c(k)*1.E6))
2408 prs_scw(k) = rhof(k)*t1_qs_qc*Ef_sw*rc(k)*smoe(k)
2409 prs_scw(k) = MIN(DBLE(rc(k)*odts), prs_scw(k))
2410 pnc_scw(k) = rhof(k)*t1_qs_qc*Ef_sw*nc(k)*smoe(k) ! Qc2M
2411 pnc_scw(k) = MIN(DBLE(nc(k)*odts), pnc_scw(k))
2414 !..Graupel collecting cloud water. In CE, assume Dc<<Dg and vtc=~0.
2415 if (rg(k).ge. r_g(1) .and. mvd_c(k).gt. D0c) then
2416 xDg = (bm_g + mu_g + 1.) * ilamg(k)
2417 vtg = rhof(k)*av_g(idx_bg(k))*cgg(6,idx_bg(k))*ogg3 * ilamg(k)**bv_g(idx_bg(k))
2418 stoke_g = mvd_c(k)*mvd_c(k)*vtg*rho_w/(9.*visco(k)*xDg)
2419 !..Rime density formula of Cober and List (1993) also used by Milbrandt and Morrison (2014).
2420 const_Ri = -1.*(mvd_c(k)*0.5E6)*vtg/MIN(-0.1,tempc)
2421 const_Ri = MAX(0.1, MIN(const_Ri, 10.))
2422 rime_dens = (0.051 + 0.114*const_Ri - 0.0055*const_Ri*const_Ri)*1000.
2423 if (stoke_g.ge.0.4 .and. stoke_g.le.10.) then
2424 Ef_gw = 0.55*ALOG10(2.51*stoke_g)
2425 elseif (stoke_g.lt.0.4) then
2427 elseif (stoke_g.gt.10) then
2430 !!! Not sure what to do here - hail increases size rapidly here below melting level.
2431 if (twet(k).gt.T_0) Ef_gw = Ef_gw*0.1
2432 t1_qg_qc = PI*.25*av_g(idx_bg(k)) * cgg(9,idx_bg(k))
2433 prg_gcw(k) = rhof(k)*t1_qg_qc*Ef_gw*rc(k)*N0_g(k) &
2434 *ilamg(k)**cge(9,idx_bg(k))
2435 pnc_gcw(k) = rhof(k)*t1_qg_qc*Ef_gw*nc(k)*N0_g(k) &
2436 *ilamg(k)**cge(9,idx_bg(k)) ! Qc2M
2437 pnc_gcw(k) = MIN(DBLE(nc(k)*odts), pnc_gcw(k))
2438 if (twet(k).lt.T_0) pbg_gcw(k) = prg_gcw(k)/rime_dens
2442 !..Snow and graupel collecting aerosols, wet scavenging.
2443 if (rs(k) .gt. r_s(1)) then
2444 Ef_sa = Eff_aero(xDs,0.04E-6,visco(k),rho(k),temp(k),'s')
2445 pna_sca(k) = rhof(k)*t1_qs_qc*Ef_sa*nwfa(k)*smoe(k)
2446 pna_sca(k) = MIN(DBLE(nwfa(k)*odts), pna_sca(k))
2448 Ef_sa = Eff_aero(xDs,0.8E-6,visco(k),rho(k),temp(k),'s')
2449 pnd_scd(k) = rhof(k)*t1_qs_qc*Ef_sa*nifa(k)*smoe(k)
2450 pnd_scd(k) = MIN(DBLE(nifa(k)*odts), pnd_scd(k))
2452 if (present(wif_input_opt)) then
2453 if (wif_input_opt .eq. 2) then
2454 Ef_sa = Eff_aero(xDs,0.0236E-6,visco(k),rho(k),temp(k),'s')
2455 pnb_scb(k) = rhof(k)*t1_qs_qc*Ef_sa*nbca(k)*smoe(k)
2456 pnb_scb(k) = MIN(DBLE(nbca(k)*odts), pnb_scb(k))
2460 if (rg(k) .gt. r_g(1)) then
2461 xDg = (bm_g + mu_g + 1.) * ilamg(k)
2462 Ef_ga = Eff_aero(xDg,0.04E-6,visco(k),rho(k),temp(k),'g')
2463 t1_qg_qc = PI*.25*av_g(idx_bg(k)) * cgg(9,idx_bg(k))
2464 pna_gca(k) = rhof(k)*t1_qg_qc*Ef_ga*nwfa(k)*N0_g(k) &
2465 *ilamg(k)**cge(9,idx_bg(k))
2466 pna_gca(k) = MIN(DBLE(nwfa(k)*odts), pna_gca(k))
2468 Ef_ga = Eff_aero(xDg,0.8E-6,visco(k),rho(k),temp(k),'g')
2469 pnd_gcd(k) = rhof(k)*t1_qg_qc*Ef_ga*nifa(k)*N0_g(k) &
2470 *ilamg(k)**cge(9,idx_bg(k))
2471 pnd_gcd(k) = MIN(DBLE(nifa(k)*odts), pnd_gcd(k))
2473 if (present(wif_input_opt)) then
2474 if (wif_input_opt .eq. 2) then
2475 Ef_ga = Eff_aero(xDg,0.0236E-6,visco(k),rho(k),temp(k),'g')
2476 pnb_gcb(k) = rhof(k)*t1_qg_qc*Ef_ga*nbca(k)*N0_g(k) &
2477 *ilamg(k)**cge(9,idx_bg(k))
2478 pnb_gcb(k) = MIN(DBLE(nbca(k)*odts), pnb_gcb(k))
2483 !..Rain collecting snow. Cannot assume Wisner (1972) approximation
2484 !.. or Mizuno (1990) approach so we solve the CE explicitly and store
2485 !.. results in lookup table.
2486 if (rr(k).ge. r_r(1)) then
2487 if (rs(k).ge. r_s(1)) then
2488 if (twet(k).lt.T_0) then
2489 prr_rcs(k) = -(tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2490 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
2491 + tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
2492 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r))
2493 prs_rcs(k) = tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2494 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
2495 - tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2496 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
2497 prg_rcs(k) = tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
2498 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2499 + tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2500 + tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
2501 prr_rcs(k) = MAX(DBLE(-rr(k)*odts), prr_rcs(k))
2502 prs_rcs(k) = MAX(DBLE(-rs(k)*odts), prs_rcs(k))
2503 prg_rcs(k) = MIN(DBLE((rr(k)+rs(k))*odts), prg_rcs(k))
2504 pnr_rcs(k) = tnr_racs1(idx_s,idx_t,idx_r1,idx_r) & ! RAIN2M
2505 + tnr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2506 + tnr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2507 + tnr_sacr2(idx_s,idx_t,idx_r1,idx_r)
2508 pnr_rcs(k) = MIN(DBLE(nr(k)*odts), pnr_rcs(k))
2509 png_rcs(k) = pnr_rcs(k)
2510 !-GT pbg_rcs(k) = prg_rcs(k)/(0.5*(rho_i+rho_s))
2511 pbg_rcs(k) = prg_rcs(k)/rho_i
2513 prs_rcs(k) = -tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2514 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2515 + tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2516 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r)
2517 prs_rcs(k) = MAX(DBLE(-rs(k)*odts), prs_rcs(k))
2518 prr_rcs(k) = -prs_rcs(k)
2522 !..Rain collecting graupel. Cannot assume Wisner (1972) approximation
2523 !.. or Mizuno (1990) approach so we solve the CE explicitly and store
2524 !.. results in lookup table.
2525 if (rg(k).ge. r_g(1)) then
2526 if (twet(k).lt.T_0) then
2527 prg_rcg(k) = tmr_racg(idx_g1,idx_g,idx_bg(k),idx_r1,idx_r) &
2528 + tcr_gacr(idx_g1,idx_g,idx_bg(k),idx_r1,idx_r)
2529 prg_rcg(k) = MIN(DBLE(rr(k)*odts), prg_rcg(k))
2530 prr_rcg(k) = -prg_rcg(k)
2531 pnr_rcg(k) = tnr_racg(idx_g1,idx_g,idx_bg(k),idx_r1,idx_r) & ! RAIN2M
2532 + tnr_gacr(idx_g1,idx_g,idx_bg(k),idx_r1,idx_r)
2533 pnr_rcg(k) = MIN(DBLE(nr(k)*odts), pnr_rcg(k))
2534 !-GT pbg_rcg(k) = prg_rcg(k)/(0.5*(rho_i+rho_g(idx_bg(k))))
2535 pbg_rcg(k) = prg_rcg(k)/rho_i
2537 prr_rcg(k) = tcg_racg(idx_g1,idx_g,idx_bg(k),idx_r1,idx_r)
2538 prr_rcg(k) = MIN(DBLE(rg(k)*odts), prr_rcg(k))
2539 prg_rcg(k) = -prr_rcg(k)
2540 png_rcg(k) = tnr_racg(idx_g1,idx_g,idx_bg(k),idx_r1,idx_r)
2541 !!! + tnr_gacr(idx_g1,idx_g,idx_bg(k),idx_r1,idx_r)
2542 png_rcg(k) = MIN(DBLE(ng(k)*odts), png_rcg(k))
2543 pbg_rcg(k) = prg_rcg(k)/rho_g(idx_bg(k))
2544 !..Put in explicit drop break-up due to collisions.
2545 pnr_rcg(k) = -5.*tnr_gacr(idx_g1,idx_g,idx_bg(k),idx_r1,idx_r) ! RAIN2M
2550 !+---+-----------------------------------------------------------------+
2551 !..Next IF block handles only those processes below 0C.
2552 !+---+-----------------------------------------------------------------+
2554 if (temp(k).lt.T_0) then
2557 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
2559 !+---+---------------- BEGIN NEW ICE NUCLEATION -----------------------+
2560 !..Freezing of supercooled water (rain or cloud) is influenced by dust
2561 !.. but still using Bigg 1953 with a temperature adjustment of a few
2562 !.. degrees depending on dust concentration. A default value by way
2563 !.. of idx_IN is 1.0 per Liter of air is used when dustyIce flag is
2564 !.. false. Next, a combination of deposition/condensation freezing
2565 !.. using DeMott et al (2010) dust nucleation when water saturated or
2566 !.. Phillips et al (2008) when below water saturation; else, without
2567 !.. dustyIce flag, use the previous Cooper (1986) temperature-dependent
2568 !.. value. Lastly, allow homogeneous freezing of deliquesced aerosols
2569 !.. following Koop et al. (2001, Nature).
2570 !.. Implemented by T. Eidhammer and G. Thompson 2012Dec18
2571 !+---+-----------------------------------------------------------------+
2573 if (dustyIce .AND. is_aerosol_aware) then
2574 xni = iceDeMott(tempc,qvs(k),qvs(k),qvsi(k),rho(k),nifa(k))
2576 xni = 1.0 *1000. ! Default is 1.0 per Liter
2579 !..Ice nuclei lookup table index.
2580 if (xni.gt. Nt_IN(1)) then
2581 niin = NINT(ALOG10(xni))
2582 do nn = niin-1, niin+1
2584 if ( (xni/10.**nn).ge.1.0 .and. &
2585 (xni/10.**nn).lt.10.0) goto 149
2588 idx_IN = INT(xni/10.**n) + 10*(n-niin2) - (n-niin2)
2589 idx_IN = MAX(1, MIN(idx_IN, ntb_IN))
2594 !..Freezing of water drops into graupel/cloud ice (Bigg 1953).
2595 if (rr(k).gt. r_r(1)) then
2596 prg_rfz(k) = tpg_qrfz(idx_r,idx_r1,idx_tc,idx_IN)*odts
2597 pri_rfz(k) = tpi_qrfz(idx_r,idx_r1,idx_tc,idx_IN)*odts
2598 pni_rfz(k) = tni_qrfz(idx_r,idx_r1,idx_tc,idx_IN)*odts
2599 pnr_rfz(k) = tnr_qrfz(idx_r,idx_r1,idx_tc,idx_IN)*odts ! RAIN2M
2600 pnr_rfz(k) = MIN(DBLE(nr(k)*odts), pnr_rfz(k))
2601 elseif (rr(k).gt. R1 .and. temp(k).lt.HGFR) then
2602 pri_rfz(k) = rr(k)*odts
2603 pni_rfz(k) = nr(k)*odts ! RAIN2M
2605 pbg_rfz(k) = prg_rfz(k)/rho_i
2607 if (rc(k).gt. r_c(1)) then
2608 pri_wfz(k) = tpi_qcfz(idx_c,idx_n,idx_tc,idx_IN)*odts
2609 pri_wfz(k) = MIN(DBLE(rc(k)*odts), pri_wfz(k))
2610 pni_wfz(k) = tni_qcfz(idx_c,idx_n,idx_tc,idx_IN)*odts
2611 pni_wfz(k) = MIN(DBLE(nc(k)*odts), pri_wfz(k)/(2.*xm0i), &
2613 elseif (rc(k).gt. R1 .and. temp(k).lt.HGFR) then
2614 pri_wfz(k) = rc(k)*odts
2615 pni_wfz(k) = nc(k)*odts
2618 !..Deposition nucleation of dust/mineral from DeMott et al (2010)
2619 !.. we may need to relax the temperature and ssati constraints.
2620 if ( (ssati(k).ge. 0.25) .or. (ssatw(k).gt. eps &
2621 .and. temp(k).lt.253.15) ) then
2622 if (dustyIce .AND. is_aerosol_aware) then
2623 xnc = iceDeMott(tempc,qv(k),qvs(k),qvsi(k),rho(k),nifa(k))
2625 xnc = MIN(250.E3, TNO*EXP(ATO*(T_0-temp(k))))
2627 xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dtsave
2628 pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts
2629 pri_inu(k) = MIN(DBLE(rate_max), xm0i*pni_inu(k))
2630 pni_inu(k) = pri_inu(k)/xm0i
2633 !..Freezing of aqueous aerosols based on Koop et al (2001, Nature)
2634 xni = ns(k)+ni(k) + (pni_rfz(k)+pni_wfz(k)+pni_inu(k))*dtsave
2635 if (is_aerosol_aware .AND. homogIce .AND. (xni.le.999.E3) &
2636 & .AND.(temp(k).lt.238).AND.(ssati(k).ge.0.4) ) then
2637 xnc = iceKoop(temp(k),qv(k),qvs(k),nwfa(k), dtsave)
2638 pni_iha(k) = xnc*odts
2639 pri_iha(k) = MIN(DBLE(rate_max), xm0i*0.1*pni_iha(k))
2640 pni_iha(k) = pri_iha(k)/(xm0i*0.1)
2642 !+---+------------------ END NEW ICE NUCLEATION -----------------------+
2645 !..Deposition/sublimation of cloud ice (Srivastava & Coen 1992).
2647 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
2649 xDi = MAX(DBLE(D0i), (bm_i + mu_i + 1.) * ilami)
2650 xmi = am_i*xDi**bm_i
2652 pri_ide(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
2653 *oig1*cig(5)*ni(k)*ilami
2655 if (pri_ide(k) .lt. 0.0) then
2656 pri_ide(k) = MAX(DBLE(-ri(k)*odts), pri_ide(k), DBLE(rate_max))
2657 pni_ide(k) = pri_ide(k)*oxmi
2658 pni_ide(k) = MAX(DBLE(-ni(k)*odts), pni_ide(k))
2660 pri_ide(k) = MIN(pri_ide(k), DBLE(rate_max))
2661 prs_ide(k) = (1.0D0-tpi_ide(idx_i,idx_i1))*pri_ide(k)
2662 pri_ide(k) = tpi_ide(idx_i,idx_i1)*pri_ide(k)
2665 !..Some cloud ice needs to move into the snow category. Use lookup
2666 !.. table that resulted from explicit bin representation of distrib.
2667 if ( (idx_i.eq. ntb_i) .or. (xDi.gt. 5.0*D0s) ) then
2668 prs_iau(k) = ri(k)*.99*odts
2669 pni_iau(k) = ni(k)*.95*odts
2670 elseif (xDi.lt. 0.1*D0s) then
2674 prs_iau(k) = tps_iaus(idx_i,idx_i1)*odts
2675 prs_iau(k) = MIN(DBLE(ri(k)*.99*odts), prs_iau(k))
2676 pni_iau(k) = tni_iaus(idx_i,idx_i1)*odts
2677 pni_iau(k) = MIN(DBLE(ni(k)*.95*odts), pni_iau(k))
2681 !..Deposition/sublimation of snow/graupel follows Srivastava & Coen
2684 C_snow = C_sqrd + (tempc+1.5)*(C_cube-C_sqrd)/(-30.+1.5)
2685 C_snow = MAX(C_sqrd, MIN(C_snow, C_cube))
2686 prs_sde(k) = C_snow*t1_subl*diffu(k)*ssati(k)*rvs &
2687 * (t1_qs_sd*smo1(k) &
2688 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
2689 if (prs_sde(k).lt. 0.) then
2690 prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k), DBLE(rate_max))
2692 prs_sde(k) = MIN(prs_sde(k), DBLE(rate_max))
2696 if (L_qg(k) .and. ssati(k).lt. -eps) then
2697 t2_qg_sd = 0.28*Sc3*SQRT(av_g(idx_bg(k))) * cgg(11,idx_bg(k))
2698 prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
2699 * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10,1) &
2700 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11,idx_bg(k)))
2701 if (prg_gde(k).lt. 0.) then
2702 prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k), DBLE(rate_max))
2703 png_gde(k) = prg_gde(k) * ng(k)/rg(k)
2705 prg_gde(k) = MIN(prg_gde(k), DBLE(rate_max))
2709 !..Snow collecting cloud ice. In CE, assume Di<<Ds and vti=~0.
2711 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
2713 xDi = MAX(DBLE(D0i), (bm_i + mu_i + 1.) * ilami)
2714 xmi = am_i*xDi**bm_i
2716 if (rs(k).ge. r_s(1)) then
2717 prs_sci(k) = t1_qs_qi*rhof(k)*Ef_si*ri(k)*smoe(k)
2718 pni_sci(k) = prs_sci(k) * oxmi
2721 !..Rain collecting cloud ice. In CE, assume Di<<Dr and vti=~0.
2722 if (rr(k).ge. r_r(1) .and. mvd_r(k).gt. 4.*xDi) then
2724 pri_rci(k) = rhof(k)*t1_qr_qi*Ef_ri*ri(k)*N0_r(k) &
2725 *((lamr+fv_r)**(-cre(9)))
2726 pnr_rci(k) = rhof(k)*t1_qr_qi*Ef_ri*ni(k)*N0_r(k) & ! RAIN2M
2727 *((lamr+fv_r)**(-cre(9)))
2728 pnr_rci(k) = MIN(DBLE(nr(k)*odts), pnr_rci(k))
2729 pni_rci(k) = pri_rci(k) * oxmi
2730 prr_rci(k) = rhof(k)*t2_qr_qi*Ef_ri*ni(k)*N0_r(k) &
2731 *((lamr+fv_r)**(-cre(8)))
2732 prr_rci(k) = MIN(DBLE(rr(k)*odts), prr_rci(k))
2733 prg_rci(k) = pri_rci(k) + prr_rci(k)
2734 pbg_rci(k) = prg_rci(k)/rho_i
2738 !..Ice multiplication from rime-splinters (Hallet & Mossop 1974).
2739 if (prg_gcw(k).gt. eps .and. tempc.gt.-8.0) then
2741 if (tempc.ge.-5.0 .and. tempc.lt.-3.0) then
2742 tf = 0.5*(-3.0 - tempc)
2743 elseif (tempc.gt.-8.0 .and. tempc.lt.-5.0) then
2744 tf = 0.33333333*(8.0 + tempc)
2746 pni_ihm(k) = 3.5E8*tf*prg_gcw(k)
2747 pri_ihm(k) = xm0i*pni_ihm(k)
2748 prs_ihm(k) = prs_scw(k)/(prs_scw(k)+prg_gcw(k)) &
2750 prg_ihm(k) = prg_gcw(k)/(prs_scw(k)+prg_gcw(k)) &
2754 !..A portion of rimed snow converts to graupel but some remains snow.
2755 !.. Interp from 15 to 95% as riming factor increases from 2.0 to 30.0
2756 !.. 0.028 came from (.95-.15)/(30.-2.). This remains ad-hoc and should
2758 if (prs_scw(k).gt.2.0*prs_sde(k) .and. &
2759 prs_sde(k).gt.eps) then
2760 r_frac = MIN(30.0D0, prs_scw(k)/prs_sde(k))
2761 g_frac = MIN(0.95, 0.15 + (r_frac-2.)*.028)
2762 vts_boost(k) = MIN(1.5, 1.1 + (r_frac-2.)*.014)
2763 prg_scw(k) = g_frac*prs_scw(k)
2764 png_scw(k) = prg_scw(k)*smo0(k)/rs(k)
2765 !..GT png_scw(k) = prg_scw(k)*ns(k)/rs(k)
2766 vts = av_s*xDs**bv_s * EXP(-fv_s*xDs)
2767 const_Ri = -1.*(mvd_c(k)*0.5E6)*vts/MIN(-0.1,tempc)
2768 const_Ri = MAX(0.1, MIN(const_Ri, 10.))
2769 rime_dens = (0.051 + 0.114*const_Ri - 0.0055*const_Ri*const_Ri)*1000.
2770 if(rime_dens .lt. 150.) then ! Idea of A. Jensen
2775 pbg_scw(k) = prg_scw(k)/(0.5*(rime_dens+rho_s))
2776 prs_scw(k) = (1. - g_frac)*prs_scw(k)
2781 !..Melt snow and graupel and enhance from collisions with liquid.
2782 !.. We also need to sublimate snow and graupel if subsaturated.
2784 prr_sml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delQvs(k)) &
2785 * (t1_qs_me*smo1(k) + t2_qs_me*rhof2(k)*vsc2(k)*smof(k))
2786 if (prr_sml(k) .gt. 0.) then
2787 prr_sml(k) = prr_sml(k) + 4218.*olfus*(twet(k)-T_0) &
2788 * (prr_rcs(k)+prs_scw(k))
2790 prr_sml(k) = MIN(DBLE(rs(k)*odts), MAX(0.D0, prr_sml(k)))
2792 if (prr_sml(k).gt.0.0) then
2793 pnr_sml(k) = smo0(k)/rs(k)*prr_sml(k) * 10.0**(-0.25*(twet(k)-T_0))
2794 elseif (ssati(k).lt. 0.) then
2795 prs_sde(k) = C_sqrd*t1_subl*diffu(k)*ssati(k)*rvs &
2796 * (t1_qs_sd*smo1(k) &
2797 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
2798 prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k))
2804 if ((rg(k)*ng(k)) .lt. 1.E-4) then
2806 N0_melt = (1.E-4/rg(k))*ogg2*lamg**cge(2,1)
2808 t2_qg_me = PI*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g(idx_bg(k))) * cgg(11,idx_bg(k))
2809 prr_gml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delQvs(k)) &
2810 * N0_melt*(t1_qg_me*ilamg(k)**cge(10,1) &
2811 + t2_qg_me*rhof2(k)*vsc2(k)*ilamg(k)**cge(11,idx_bg(k)))
2812 ! if (prr_gml(k) .gt. 0.) then
2813 ! prr_gml(k) = prr_gml(k) + 4218.*olfus*(twet(k)-T_0) &
2814 ! * (prr_rcg(k)+prg_gcw(k))
2816 prr_gml(k) = MIN(DBLE(rg(k)*odts), MAX(0.D0, prr_gml(k)))
2817 ! pnr_gml(k) = N0_g(k)*cgg(2)*ilamg(k)**cge(2) / rg(k) & ! RAIN2M
2818 ! * prr_gml(k) * 10.0**(-0.5*tempc)
2820 if (prr_gml(k) .gt. 0.0) then
2821 melt_f = MAX(0.05, MIN(prr_gml(k)*dt/rg(k),1.0))
2822 !..1000 is density water, 50 is lower limit (max ice density is 800)
2823 pbg_gml(k) = prr_gml(k) / MAX(MIN(melt_f*rho_g(idx_bg(k)),1000.),50.)
2824 !-GT pnr_gml(k) = prr_gml(k)*ng(k)/rg(k)
2825 pnr_gml(k) = prr_gml(k)*ng(k)/rg(k) * 10.0**(-0.33*(twet(k)-T_0))
2826 elseif (ssati(k).lt. 0.) then
2827 t2_qg_sd = 0.28*Sc3*SQRT(av_g(idx_bg(k))) * cgg(11,idx_bg(k))
2828 prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
2829 * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10,1) &
2830 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11,idx_bg(k)))
2831 prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k))
2832 png_gde(k) = prg_gde(k) * ng(k)/rg(k)
2836 !.. This change will be required if users run adaptive time step that
2837 !.. results in delta-t that is generally too long to allow cloud water
2838 !.. collection by snow/graupel above melting temperature.
2839 !.. Credit to Bjorn-Egil Nygaard for discovering.
2840 if (dt .gt. 120.) then
2841 prr_rcw(k)=prr_rcw(k)+prs_scw(k)+prg_gcw(k)
2848 if (.not. is_hail_aware) idx_bg(k) = idx_bg1
2853 !+---+-----------------------------------------------------------------+
2854 !..Ensure we do not deplete more hydrometeor species than exists.
2855 !+---+-----------------------------------------------------------------+
2858 !..If ice supersaturated, ensure sum of depos growth terms does not
2859 !.. deplete more vapor than possibly exists. If subsaturated, limit
2860 !.. sum of sublimation terms such that vapor does not reproduce ice
2862 sump = pri_inu(k) + pri_ide(k) + prs_ide(k) &
2863 + prs_sde(k) + prg_gde(k) + pri_iha(k)
2864 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
2865 if ( (sump.gt. eps .and. sump.gt. rate_max) .or. &
2866 (sump.lt. -eps .and. sump.lt. rate_max) ) then
2867 ratio = rate_max/sump
2868 pri_inu(k) = pri_inu(k) * ratio
2869 pri_ide(k) = pri_ide(k) * ratio
2870 pni_ide(k) = pni_ide(k) * ratio
2871 prs_ide(k) = prs_ide(k) * ratio
2872 prs_sde(k) = prs_sde(k) * ratio
2873 prg_gde(k) = prg_gde(k) * ratio
2874 pri_iha(k) = pri_iha(k) * ratio
2877 !..Cloud water conservation.
2878 sump = -prr_wau(k) - pri_wfz(k) - prr_rcw(k) &
2879 - prs_scw(k) - prg_scw(k) - prg_gcw(k)
2880 rate_max = -rc(k)*odts
2881 if (sump.lt. rate_max .and. L_qc(k)) then
2882 ratio = rate_max/sump
2883 prr_wau(k) = prr_wau(k) * ratio
2884 pri_wfz(k) = pri_wfz(k) * ratio
2885 prr_rcw(k) = prr_rcw(k) * ratio
2886 prs_scw(k) = prs_scw(k) * ratio
2887 prg_scw(k) = prg_scw(k) * ratio
2888 prg_gcw(k) = prg_gcw(k) * ratio
2891 !..Cloud ice conservation.
2892 sump = pri_ide(k) - prs_iau(k) - prs_sci(k) &
2894 rate_max = -ri(k)*odts
2895 if (sump.lt. rate_max .and. L_qi(k)) then
2896 ratio = rate_max/sump
2897 pri_ide(k) = pri_ide(k) * ratio
2898 prs_iau(k) = prs_iau(k) * ratio
2899 prs_sci(k) = prs_sci(k) * ratio
2900 pri_rci(k) = pri_rci(k) * ratio
2903 !..Rain conservation.
2904 sump = -prg_rfz(k) - pri_rfz(k) - prr_rci(k) &
2905 + prr_rcs(k) + prr_rcg(k)
2906 rate_max = -rr(k)*odts
2907 if (sump.lt. rate_max .and. L_qr(k)) then
2908 ratio = rate_max/sump
2909 prg_rfz(k) = prg_rfz(k) * ratio
2910 pri_rfz(k) = pri_rfz(k) * ratio
2911 prr_rci(k) = prr_rci(k) * ratio
2912 prr_rcs(k) = prr_rcs(k) * ratio
2913 prr_rcg(k) = prr_rcg(k) * ratio
2916 !..Snow conservation.
2917 sump = prs_sde(k) - prs_ihm(k) - prr_sml(k) &
2919 rate_max = -rs(k)*odts
2920 if (sump.lt. rate_max .and. L_qs(k)) then
2921 ratio = rate_max/sump
2922 prs_sde(k) = prs_sde(k) * ratio
2923 prs_ihm(k) = prs_ihm(k) * ratio
2924 prr_sml(k) = prr_sml(k) * ratio
2925 prs_rcs(k) = prs_rcs(k) * ratio
2928 !..Graupel conservation.
2929 sump = prg_gde(k) - prg_ihm(k) - prr_gml(k) &
2931 rate_max = -rg(k)*odts
2932 if (sump.lt. rate_max .and. L_qg(k)) then
2933 ratio = rate_max/sump
2934 prg_gde(k) = prg_gde(k) * ratio
2935 prg_ihm(k) = prg_ihm(k) * ratio
2936 prr_gml(k) = prr_gml(k) * ratio
2937 prg_rcg(k) = prg_rcg(k) * ratio
2940 !..Re-enforce proper mass conservation for subsequent elements in case
2941 !.. any of the above terms were altered. Thanks P. Blossey. 2009Sep28
2942 pri_ihm(k) = prs_ihm(k) + prg_ihm(k)
2943 ratio = MIN( ABS(prr_rcg(k)), ABS(prg_rcg(k)) )
2944 prr_rcg(k) = ratio * SIGN(1.0, SNGL(prr_rcg(k)))
2945 prg_rcg(k) = -prr_rcg(k)
2946 if (twet(k).gt.T_0) then
2947 ratio = MIN( ABS(prr_rcs(k)), ABS(prs_rcs(k)) )
2948 prr_rcs(k) = ratio * SIGN(1.0, SNGL(prr_rcs(k)))
2949 prs_rcs(k) = -prr_rcs(k)
2954 !+---+-----------------------------------------------------------------+
2955 !..Calculate tendencies of all species but constrain the number of ice
2956 !.. to reasonable values.
2957 !+---+-----------------------------------------------------------------+
2960 lfus2 = lsub - lvap(k)
2962 !..Aerosol number tendency
2963 if (is_aerosol_aware) then
2964 nwfaten(k) = nwfaten(k) - (pna_rca(k) + pna_sca(k) &
2965 + pna_gca(k) + pni_iha(k)) * orho
2966 nifaten(k) = nifaten(k) - (pnd_rcd(k) + pnd_scd(k) &
2967 + pnd_gcd(k)) * orho
2968 if (wif_input_opt .eq. 2) then
2969 nbcaten(k) = nbcaten(k) - (pnb_rcb(k) + pnb_scb(k) &
2970 + pnb_gcb(k)) * orho
2975 nifaten(k) = nifaten(k) - pni_inu(k)*orho
2981 !..Water vapor tendency
2982 qvten(k) = qvten(k) + (-pri_inu(k) - pri_iha(k) - pri_ide(k) &
2983 - prs_ide(k) - prs_sde(k) - prg_gde(k)) &
2986 !..Cloud water tendency
2987 qcten(k) = qcten(k) + (-prr_wau(k) - pri_wfz(k) &
2988 - prr_rcw(k) - prs_scw(k) - prg_scw(k) &
2992 !..Cloud water number tendency
2993 ncten(k) = ncten(k) + (-pnc_wau(k) - pnc_rcw(k) &
2994 - pni_wfz(k) - pnc_scw(k) - pnc_gcw(k)) &
2997 !..Cloud water mass/number balance; keep mass-wt mean size between
2998 !.. 1 and 50 microns. Also no more than Nt_c_max drops total.
2999 xrc=MAX(R1, (qc1d(k) + qcten(k)*dtsave)*rho(k))
3000 xnc=MAX(2., (nc1d(k) + ncten(k)*dtsave)*rho(k))
3001 if (xrc .gt. R1) then
3002 nu_c = MIN(15, NINT(1000.E6/xnc) + 2)
3003 lamc = (xnc*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
3004 xDc = (bm_r + nu_c + 1.) / lamc
3005 if (xDc.lt. D0c) then
3006 lamc = cce(2,nu_c)/D0c
3007 xnc = ccg(1,nu_c)*ocg2(nu_c)*xrc/am_r*lamc**bm_r
3008 ncten(k) = (xnc-nc1d(k)*rho(k))*odts*orho
3009 elseif (xDc.gt. D0r*2.) then
3010 lamc = cce(2,nu_c)/(D0r*2.)
3011 xnc = ccg(1,nu_c)*ocg2(nu_c)*xrc/am_r*lamc**bm_r
3012 ncten(k) = (xnc-nc1d(k)*rho(k))*odts*orho
3015 ncten(k) = -nc1d(k)*odts
3017 xnc=MAX(0.,(nc1d(k) + ncten(k)*dtsave)*rho(k))
3018 if (xnc.gt.Nt_c_max) &
3019 ncten(k) = (Nt_c_max-nc1d(k)*rho(k))*odts*orho
3021 !..Cloud ice mixing ratio tendency
3022 qiten(k) = qiten(k) + (pri_inu(k) + pri_iha(k) + pri_ihm(k) &
3023 + pri_wfz(k) + pri_rfz(k) + pri_ide(k) &
3024 - prs_iau(k) - prs_sci(k) - pri_rci(k)) &
3027 !..Cloud ice number tendency.
3028 niten(k) = niten(k) + (pni_inu(k) + pni_iha(k) + pni_ihm(k) &
3029 + pni_wfz(k) + pni_rfz(k) + pni_ide(k) &
3030 - pni_iau(k) - pni_sci(k) - pni_rci(k)) &
3033 !..Cloud ice mass/number balance; keep mass-wt mean size between
3034 !.. 5 and 300 microns. Also no more than 500 xtals per liter.
3035 xri=MAX(R1,(qi1d(k) + qiten(k)*dtsave)*rho(k))
3036 xni=MAX(R2,(ni1d(k) + niten(k)*dtsave)*rho(k))
3037 if (xri.gt. R1) then
3038 lami = (am_i*cig(2)*oig1*xni/xri)**obmi
3040 xDi = (bm_i + mu_i + 1.) * ilami
3041 if (xDi.lt. 5.E-6) then
3043 xni = MIN(999.D3, cig(1)*oig2*xri/am_i*lami**bm_i)
3044 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
3045 elseif (xDi.gt. 300.E-6) then
3046 lami = cie(2)/300.E-6
3047 xni = cig(1)*oig2*xri/am_i*lami**bm_i
3048 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
3051 niten(k) = -ni1d(k)*odts
3053 xni=MAX(0.,(ni1d(k) + niten(k)*dtsave)*rho(k))
3054 if (xni.gt.999.E3) &
3055 niten(k) = (999.E3-ni1d(k)*rho(k))*odts*orho
3058 qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) &
3059 + prr_sml(k) + prr_gml(k) + prr_rcs(k) &
3060 + prr_rcg(k) - prg_rfz(k) &
3061 - pri_rfz(k) - prr_rci(k)) &
3064 !..Rain number tendency
3065 nrten(k) = nrten(k) + (pnr_wau(k) + pnr_sml(k) + pnr_gml(k) &
3066 - (pnr_rfz(k) + pnr_rcr(k) + pnr_rcg(k) &
3067 + pnr_rcs(k) + pnr_rci(k) + pni_rfz(k)) ) &
3070 !..Rain mass/number balance; keep median volume diameter between
3071 !.. 37 microns (D0r*0.75) and 2.5 mm.
3072 xrr=MAX(R1,(qr1d(k) + qrten(k)*dtsave)*rho(k))
3073 xnr=MAX(R2,(nr1d(k) + nrten(k)*dtsave)*rho(k))
3074 if (xrr.gt. R1) then
3075 lamr = (am_r*crg(3)*org2*xnr/xrr)**obmr
3076 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3077 if (mvd_r(k) .gt. 2.5E-3) then
3079 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3080 xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
3081 nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
3082 elseif (mvd_r(k) .lt. D0r*0.75) then
3084 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3085 xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
3086 nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
3089 qrten(k) = -qr1d(k)*odts
3090 nrten(k) = -nr1d(k)*odts
3094 qsten(k) = qsten(k) + (prs_iau(k) + prs_sde(k) + prs_sci(k) &
3095 + prs_scw(k) + prs_rcs(k) + prs_ide(k) &
3096 - prs_ihm(k) - prr_sml(k)) &
3100 qgten(k) = qgten(k) + (prg_scw(k) + prg_rfz(k) &
3101 + prg_gde(k) + prg_rcg(k) + prg_gcw(k) &
3102 + prg_rci(k) + prg_rcs(k) - prg_ihm(k) &
3106 !..Graupel number tendency
3107 ngten(k) = ngten(k) + (png_scw(k) + pnr_rfz(k) - png_rcg(k) &
3108 + pnr_rci(k) + png_rcs(k) + png_gde(k) &
3109 - pnr_gml(k)) * orho
3111 !..Graupel volume mixing ratio tendency
3112 qbten(k) = qbten(k) + (pbg_scw(k) + pbg_rfz(k) &
3113 + pbg_gcw(k) + pbg_rci(k) + pbg_rcs(k) &
3114 + pbg_rcg(k) + pbg_sml(k) - pbg_gml(k) &
3115 + (prg_gde(k) - prg_ihm(k)) /rho_g(idx_bg(k)) ) &
3118 !..Graupel mass/number balance; keep its median volume diameter between
3119 !.. 3.0 times minimum size (D0g) and 25 mm.
3120 xrg=MAX(R1,(qg1d(k) + qgten(k)*dtsave)*rho(k))
3121 xng=MAX(R2,(ng1d(k) + ngten(k)*dtsave)*rho(k))
3122 xrb=MAX(xrg/rho(k)/rho_g(NRHG),(qb1d(k) + qbten(k)*dtsave))
3123 xrb=MIN(xrg/rho(k)/rho_g(1), xrb)
3125 if (xrg .gt. R1) then
3126 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*xng/xrg)**obmg
3127 mvd_g(k) = (3.0 + mu_g + 0.672) / lamg
3130 xrho_g = MAX(rho_g(1),MIN(rg(k)/rho(k)/rb(k),rho_g(NRHG)))
3131 afall = a_coeff*((4.0*xrho_g*9.8)/(3.0*rho(k)))**b_coeff
3132 afall = afall * visco(k)**(1.0-2.0*b_coeff)
3133 vtg1 = rhof(k)*afall*cgg(6,idx_bg(k))*ogg3 / lamg**bfall
3134 vtg2 = rhof(k)*afall*cgg(7,idx_bg(k))/cgg(12,idx_bg(k)) / lamg**bfall
3135 write(mp_debug, *) 'DEBUG: xrg, xng, xrb, ns, idx_bg, rho_g, mvd, vtg: ', k, xrg, xng, xrb, ns(k), idx_bg(k), xrg/rho(k)/xrb, mvd_g(k)*1000., vtg1, vtg2
3136 CALL wrf_debug(0, mp_debug)
3137 write(mp_debug, *) 'DEBUG: qgten: ', qgten(k), prg_scw(k), prg_rfz(k), prg_gde(k), prg_rcg(k), prg_gcw(k), prg_rci(k), prg_rcs(k), prg_ihm(k), prr_gml(k)
3138 CALL wrf_debug(0, mp_debug)
3139 write(mp_debug, *) 'DEBUG: ngten: ', ngten(k), png_scw(k), pnr_rfz(k), png_rcg(k), pnr_rci(k), png_rcs(k), pnr_gml(k)
3140 CALL wrf_debug(0, mp_debug)
3141 write(mp_debug, *) 'DEBUG: qbten: ', qbten(k), pbg_scw(k), pbg_rfz(k), pbg_gcw(k), pbg_rci(k), pbg_rcs(k), pbg_rcg(k), pbg_sml(k), pbg_gml(k)
3142 CALL wrf_debug(0, mp_debug)
3145 if (mvd_g(k) .gt. 25.4E-3) then
3147 lamg = (3.0 + mu_g + 0.672) / mvd_g(k)
3148 xng = cgg(2,1)*ogg3*xrg*lamg**bm_g / am_g(idx_bg(k))
3149 ngten(k) = (xng-ng1d(k)*rho(k))*odts*orho
3150 elseif (mvd_g(k) .lt. D0r) then
3152 lamg = (3.0 + mu_g + 0.672) / mvd_g(k)
3153 xng = cgg(2,1)*ogg3*xrg*lamg**bm_g / am_g(idx_bg(k))
3154 ngten(k) = (xng-ng1d(k)*rho(k))*odts*orho
3158 qgten(k) = -qg1d(k)*odts
3159 ngten(k) = -ng1d(k)*odts
3160 qbten(k) = -qb1d(k)*odts
3163 !..Temperature tendency
3164 if (temp(k).lt.T_0) then
3166 + ( lsub*ocp(k)*(pri_inu(k) + pri_ide(k) &
3167 + prs_ide(k) + prs_sde(k) &
3168 + prg_gde(k) + pri_iha(k)) &
3169 + lfus2*ocp(k)*(pri_wfz(k) + pri_rfz(k) &
3170 + prg_rfz(k) + prs_scw(k) &
3171 + prg_scw(k) + prg_gcw(k) &
3172 + prg_rcs(k) + prs_rcs(k) &
3173 + prr_rci(k) + prg_rcg(k)) &
3177 + ( lfus*ocp(k)*(-prr_sml(k) - prr_gml(k) &
3178 - prr_rcg(k) - prr_rcs(k)) &
3179 + lsub*ocp(k)*(prs_sde(k) + prg_gde(k)) &
3185 !+---+-----------------------------------------------------------------+
3186 !..Update variables for TAU+1 before condensation & sedimention.
3187 !+---+-----------------------------------------------------------------+
3189 temp(k) = t1d(k) + DT*tten(k)
3191 tempc = temp(k) - 273.15
3192 qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
3193 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
3194 rhof(k) = SQRT(RHO_NOT/rho(k))
3195 rhof2(k) = SQRT(rhof(k))
3196 qvs(k) = rslf(pres(k), temp(k))
3197 ssatw(k) = qv(k)/qvs(k) - 1.
3198 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
3199 diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
3200 if (tempc .ge. 0.0) then
3201 visco(k) = (1.718+0.0049*tempc)*1.0E-5
3203 visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
3205 vsc2(k) = SQRT(rho(k)/visco(k))
3206 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
3207 tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
3208 ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
3209 lvt2(k)=lvap(k)*lvap(k)*ocp(k)*oRv*otemp*otemp
3211 nwfa(k) = MAX(11.1E6, (nwfa1d(k) + nwfaten(k)*DT)*rho(k))
3215 if ((qc1d(k) + qcten(k)*DT) .gt. R1) then
3216 rc(k) = (qc1d(k) + qcten(k)*DT)*rho(k)
3217 nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max))
3218 if (.NOT. is_aerosol_aware) nc(k) = Nt_c
3226 if ((qi1d(k) + qiten(k)*DT) .gt. R1) then
3227 ri(k) = (qi1d(k) + qiten(k)*DT)*rho(k)
3228 ni(k) = MAX(R2, (ni1d(k) + niten(k)*DT)*rho(k))
3236 if ((qr1d(k) + qrten(k)*DT) .gt. R1) then
3237 rr(k) = (qr1d(k) + qrten(k)*DT)*rho(k)
3238 nr(k) = MAX(R2, (nr1d(k) + nrten(k)*DT)*rho(k))
3240 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3241 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3242 if (mvd_r(k) .gt. 2.5E-3) then
3244 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3245 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
3246 elseif (mvd_r(k) .lt. D0r*0.75) then
3248 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3249 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
3257 if ((qs1d(k) + qsten(k)*DT) .gt. R1) then
3258 rs(k) = (qs1d(k) + qsten(k)*DT)*rho(k)
3266 if (is_hail_aware) then
3268 if ((qg1d(k) + qgten(k)*DT) .gt. R1) then
3270 rg(k) = (qg1d(k) + qgten(k)*DT)*rho(k)
3271 ng(k) = MAX(R2, (ng1d(k) + ngten(k)*DT)*rho(k))
3272 rb(k) = MAX(rg(k)/rho(k)/rho_g(NRHG), qb1d(k) + qbten(k)*DT)
3273 rb(k) = MIN(rg(k)/rho(k)/rho_g(1), rb(k))
3274 idx_bg(k) = MAX(1,MIN(NINT(rg(k)/rho(k)/rb(k) *0.01)+1,NRHG))
3278 rb(k) = R1/rho(k)/rho_g(NRHG)
3288 if ((qg1d(k) + qgten(k)*DT) .gt. R1) then
3289 rg(k) = (qg1d(k) + qgten(k)*DT)*rho(k)
3290 ygra1 = alog10(max(1.E-9, rg(k)))
3291 zans1 = 3.0 + 2./7.*(ygra1+8.)
3292 zans1 = MAX(2., MIN(zans1, 6.))
3293 N0_exp = 10.**(zans1)
3294 lam_exp = (N0_exp*am_g(idx_bg(k))*cgg(1,1)/rg(k))**oge1
3295 lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg
3296 ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k))
3297 rb(k) = rg(k)/rho(k)/rho_g(idx_bg(k))
3301 rb(k) = R1/rho(k)/rho_g(NRHG)
3307 !+---+-----------------------------------------------------------------+
3308 !..With tendency-updated mixing ratios, recalculate snow moments and
3309 !.. intercepts/slopes of graupel and rain.
3310 !+---+-----------------------------------------------------------------+
3311 if (.not. iiwarm) then
3319 if (.not. L_qs(k)) CYCLE
3320 tc0 = MIN(-0.1, temp(k)-273.15)
3321 smob(k) = rs(k)*oams
3323 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
3324 !.. then we must compute actual 2nd moment and use as reference.
3325 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
3328 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
3329 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
3330 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
3331 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
3332 + sa(10)*bm_s*bm_s*bm_s
3334 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
3335 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
3336 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
3337 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
3338 + sb(10)*bm_s*bm_s*bm_s
3339 smo2(k) = (smob(k)/a_)**(1./b_)
3342 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
3343 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
3344 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
3345 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
3346 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
3347 + sa(10)*cse(1)*cse(1)*cse(1)
3349 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
3350 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
3351 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
3352 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
3353 smoc(k) = a_ * smo2(k)**b_
3355 !..Calculate bm_s+bv_s (th) moment. Useful for sedimentation.
3356 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(14) &
3357 + sa(4)*tc0*cse(14) + sa(5)*tc0*tc0 &
3358 + sa(6)*cse(14)*cse(14) + sa(7)*tc0*tc0*cse(14) &
3359 + sa(8)*tc0*cse(14)*cse(14) + sa(9)*tc0*tc0*tc0 &
3360 + sa(10)*cse(14)*cse(14)*cse(14)
3362 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(14) + sb(4)*tc0*cse(14) &
3363 + sb(5)*tc0*tc0 + sb(6)*cse(14)*cse(14) &
3364 + sb(7)*tc0*tc0*cse(14) + sb(8)*tc0*cse(14)*cse(14) &
3365 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(14)*cse(14)*cse(14)
3366 smod(k) = a_ * smo2(k)**b_
3369 !+---+-----------------------------------------------------------------+
3370 !..Calculate y-intercept, slope values for graupel.
3371 !+---+-----------------------------------------------------------------+
3374 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg
3376 N0_g(k) = ng(k)*ogg2*lamg**cge(2,1)
3381 !+---+-----------------------------------------------------------------+
3382 !..Calculate y-intercept, slope values for rain.
3383 !+---+-----------------------------------------------------------------+
3385 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3387 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3388 N0_r(k) = nr(k)*org2*lamr**cre(2)
3391 !+---+-----------------------------------------------------------------+
3392 !..Cloud water condensation and evaporation. Nucleate cloud droplets
3393 !.. using explicit CCN aerosols with hygroscopicity like sulfates using
3394 !.. parcel model lookup table results provided by T. Eidhammer. Evap
3395 !.. drops using calculation of max drop size capable of evaporating in
3396 !.. single timestep and explicit number of drops smaller than Dc_star
3397 !.. from lookup table.
3398 !+---+-----------------------------------------------------------------+
3401 if ( (ssatw(k).gt. eps) .or. (ssatw(k).lt. -eps .and. &
3403 clap = (qv(k)-qvs(k))/(1. + lvt2(k)*qvs(k))
3405 fcd = qvs(k)* EXP(lvt2(k)*clap) - qv(k) + clap
3406 dfcd = qvs(k)*lvt2(k)* EXP(lvt2(k)*clap) + 1.
3407 clap = clap - fcd/dfcd
3409 xrc = rc(k) + clap*rho(k)
3411 if (xrc.gt. R1) then
3412 prw_vcd(k) = clap*odt
3413 !+---+-----------------------------------------------------------------+ ! DROPLET NUCLEATION
3414 if (clap .gt. eps) then
3415 if (is_aerosol_aware) then
3416 xnc = MAX(2., activ_ncloud(temp(k), w1d(k), nwfa(k)))
3420 pnc_wcd(k) = 0.5*(xnc-nc(k) + abs(xnc-nc(k)))*odts*orho
3422 !+---+-----------------------------------------------------------------+ ! EVAPORATION
3423 elseif (clap .lt. -eps .AND. ssatw(k).lt.-1.E-6 .AND. is_aerosol_aware) then
3424 tempc = temp(k) - 273.15
3427 rvs_p = rvs*otemp*(lvap(k)*otemp*oRv - 1.)
3428 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*oRv - 1.) &
3429 *otemp*(lvap(k)*otemp*oRv - 1.) &
3430 + (-2.*lvap(k)*otemp*otemp*otemp*oRv) &
3432 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
3433 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
3434 * rvs_pp/rvs_p * rvs/rvs_p
3435 alphsc = MAX(1.E-9, alphsc)
3437 if (abs(xsat).lt. 1.E-9) xsat=0.
3438 t1_evap = 2.*PI*( 1.0 - alphsc*xsat &
3439 + 2.*alphsc*alphsc*xsat*xsat &
3440 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
3443 Dc_star = DSQRT(-2.D0*DT * t1_evap/(2.*PI) &
3444 * 4.*diffu(k)*ssatw(k)*rvs/rho_w)
3445 idx_d = MAX(1, MIN(INT(1.E6*Dc_star), nbc))
3447 idx_n = NINT(1.0 + FLOAT(nbc) * DLOG(nc(k)/t_Nc(1)) / nic1)
3448 idx_n = MAX(1, MIN(idx_n, nbc))
3450 !..Cloud water lookup table index.
3451 if (rc(k).gt. r_c(1)) then
3452 nic = NINT(ALOG10(rc(k)))
3453 do nn = nic-1, nic+1
3455 if ( (rc(k)/10.**nn).ge.1.0 .and. &
3456 (rc(k)/10.**nn).lt.10.0) goto 159
3459 idx_c = INT(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
3460 idx_c = MAX(1, MIN(idx_c, ntb_c))
3465 !prw_vcd(k) = MAX(DBLE(-rc(k)*orho*odt), &
3466 ! -tpc_wev(idx_d, idx_c, idx_n)*orho*odt)
3467 prw_vcd(k) = MAX(DBLE(-rc(k)*0.99*orho*odt), prw_vcd(k))
3468 pnc_wcd(k) = MAX(DBLE(-nc(k)*0.99*orho*odt), &
3469 DBLE(-tnc_wev(idx_d, idx_c, idx_n)*orho*odt))
3473 prw_vcd(k) = -rc(k)*orho*odt
3474 pnc_wcd(k) = -nc(k)*orho*odt
3477 !+---+-----------------------------------------------------------------+
3479 qvten(k) = qvten(k) - prw_vcd(k)
3480 qcten(k) = qcten(k) + prw_vcd(k)
3481 ncten(k) = ncten(k) + pnc_wcd(k)
3482 nwfaten(k) = nwfaten(k) - pnc_wcd(k)
3483 tten(k) = tten(k) + lvap(k)*ocp(k)*prw_vcd(k)*(1-IFDRY)
3484 rc(k) = MAX(R1, (qc1d(k) + DT*qcten(k))*rho(k))
3485 if (rc(k).eq.R1) L_qc(k) = .false.
3486 nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max))
3487 if (.NOT. is_aerosol_aware) nc(k) = Nt_c
3488 qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
3489 temp(k) = t1d(k) + DT*tten(k)
3490 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
3491 qvs(k) = rslf(pres(k), temp(k))
3492 ssatw(k) = qv(k)/qvs(k) - 1.
3496 !+---+-----------------------------------------------------------------+
3497 !.. If still subsaturated, allow rain to evaporate, following
3498 !.. Srivastava & Coen (1992).
3499 !+---+-----------------------------------------------------------------+
3501 if ( (ssatw(k).lt. -eps) .and. L_qr(k) &
3502 .and. (.not.(prw_vcd(k).gt. 0.)) ) then
3503 tempc = temp(k) - 273.15
3506 rhof(k) = SQRT(RHO_NOT*orho)
3507 rhof2(k) = SQRT(rhof(k))
3508 diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
3509 if (tempc .ge. 0.0) then
3510 visco(k) = (1.718+0.0049*tempc)*1.0E-5
3512 visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
3514 vsc2(k) = SQRT(rho(k)/visco(k))
3515 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
3516 tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
3517 ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
3520 rvs_p = rvs*otemp*(lvap(k)*otemp*oRv - 1.)
3521 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*oRv - 1.) &
3522 *otemp*(lvap(k)*otemp*oRv - 1.) &
3523 + (-2.*lvap(k)*otemp*otemp*otemp*oRv) &
3525 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
3526 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
3527 * rvs_pp/rvs_p * rvs/rvs_p
3528 alphsc = MAX(1.E-9, alphsc)
3529 xsat = MIN(-1.E-9, ssatw(k))
3530 t1_evap = 2.*PI*( 1.0 - alphsc*xsat &
3531 + 2.*alphsc*alphsc*xsat*xsat &
3532 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
3536 !..Rapidly eliminate near zero values when low humidity (<95%)
3537 if (qv(k)/qvs(k) .lt. 0.95 .AND. rr(k)*orho.le.1.E-8) then
3538 prv_rev(k) = rr(k)*orho*odts
3540 prv_rev(k) = t1_evap*diffu(k)*(-ssatw(k))*N0_r(k)*rvs &
3541 * (t1_qr_ev*ilamr(k)**cre(10) &
3542 + t2_qr_ev*vsc2(k)*rhof2(k)*((lamr+0.5*fv_r)**(-cre(11))))
3543 rate_max = MIN((rr(k)*orho*odts), (qvs(k)-qv(k))*odts)
3544 prv_rev(k) = MIN(DBLE(rate_max), prv_rev(k)*orho)
3546 !..TEST: G. Thompson 10 May 2013
3547 !..Reduce the rain evaporation in same places as melting graupel occurs.
3548 !..Rationale: falling and simultaneous melting graupel in subsaturated
3549 !..regions will not melt as fast because particle temperature stays
3550 !..at 0C. Also not much shedding of the water from the graupel so
3551 !..likely that the water-coated graupel evaporating much slower than
3552 !..if the water was immediately shed off.
3553 IF (prr_gml(k).gt.0.0) THEN
3554 eva_factor = MIN(1.0, 0.01+(0.99-0.01)*(tempc/20.0))
3555 prv_rev(k) = prv_rev(k)*eva_factor
3559 pnr_rev(k) = MIN(DBLE(nr(k)*0.99*orho*odts), & ! RAIN2M
3560 prv_rev(k) * nr(k)/rr(k))
3562 qrten(k) = qrten(k) - prv_rev(k)
3563 qvten(k) = qvten(k) + prv_rev(k)
3564 nrten(k) = nrten(k) - pnr_rev(k)
3565 nwfaten(k) = nwfaten(k) + pnr_rev(k)
3566 tten(k) = tten(k) - lvap(k)*ocp(k)*prv_rev(k)*(1-IFDRY)
3568 rr(k) = MAX(R1, (qr1d(k) + DT*qrten(k))*rho(k))
3569 qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
3570 nr(k) = MAX(R2, (nr1d(k) + DT*nrten(k))*rho(k))
3571 temp(k) = t1d(k) + DT*tten(k)
3572 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
3575 #if ( WRF_CHEM == 1 )
3576 if( wetscav_on ) then
3578 evapprod(k) = prv_rev(k) - (min(zeroD0,prs_sde(k)) + &
3579 min(zeroD0,prg_gde(k)))
3580 rainprod(k) = prr_wau(k) + prr_rcw(k) + prs_scw(k) + &
3581 prg_scw(k) + prs_iau(k) + &
3582 prg_gcw(k) + prs_sci(k) + &
3588 !+---+-----------------------------------------------------------------+
3589 !..Find max terminal fallspeed (distribution mass-weighted mean
3590 !.. velocity) and use it to determine if we need to split the timestep
3591 !.. (var nstep>1). Either way, only bother to do sedimentation below
3592 !.. 1st level that contains any sedimenting particles (k=ksed1 on down).
3593 !.. New in v3.0+ is computing separate for rain, ice, snow, and
3594 !.. graupel species thus making code faster with credit to J. Schmidt.
3595 !+---+-----------------------------------------------------------------+
3599 do k = kte+1, kts, -1
3611 if (ANY(L_qr .eqv. .true.)) then
3614 rhof(k) = SQRT(RHO_NOT/rho(k))
3616 if (rr(k).gt. R1) then
3617 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3618 vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) &
3619 *((lamr+fv_r)**(-cre(6)))
3621 ! First below is technically correct:
3622 ! vtr = rhof(k)*av_r*crg(5)*org2 * lamr**cre(2) &
3623 ! *((lamr+fv_r)**(-cre(5)))
3624 ! Test: make number fall faster (but still slower than mass)
3625 ! Goal: less prominent size sorting
3626 vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) &
3627 *((lamr+fv_r)**(-cre(7)))
3631 vtnrk(k) = vtnrk(k+1)
3634 if (MAX(vtrk(k),vtnrk(k)) .gt. 1.E-3) then
3635 ksed1(1) = MAX(ksed1(1), k)
3636 delta_tp = dzq(k)/(MAX(vtrk(k),vtnrk(k)))
3637 nstep = MAX(nstep, INT(DT/delta_tp + 1.))
3640 if (ksed1(1) .eq. kte) ksed1(1) = kte-1
3641 if (nstep .gt. 0) onstep(1) = 1./REAL(nstep)
3644 !+---+-----------------------------------------------------------------+
3646 if (ANY(L_qc .eqv. .true.)) then
3649 if (rc(k) .gt. R2) ksed1(5) = k
3650 hgt_agl = hgt_agl + dzq(k)
3651 if (hgt_agl .gt. 500.0) goto 151
3655 do k = ksed1(5), kts, -1
3657 if (rc(k) .gt. R1 .and. w1d(k) .lt. 1.E-1) then
3658 nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2)
3659 lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
3661 vtc = rhof(k)*av_c*ccg(5,nu_c)*ocg2(nu_c) * ilamc**bv_c
3663 vtc = rhof(k)*av_c*ccg(4,nu_c)*ocg1(nu_c) * ilamc**bv_c
3669 !+---+-----------------------------------------------------------------+
3671 if (.not. iiwarm) then
3673 if (ANY(L_qi .eqv. .true.)) then
3678 if (ri(k).gt. R1) then
3679 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
3681 vti = rhof(k)*av_i*cig(3)*oig2 * ilami**bv_i
3683 ! First below is technically correct:
3684 ! vti = rhof(k)*av_i*cig(4)*oig1 * ilami**bv_i
3685 ! Goal: less prominent size sorting
3686 vti = rhof(k)*av_i*cig(6)/cig(7) * ilami**bv_i
3690 vtnik(k) = vtnik(k+1)
3693 if (vtik(k) .gt. 1.E-3) then
3694 ksed1(2) = MAX(ksed1(2), k)
3695 delta_tp = dzq(k)/vtik(k)
3696 nstep = MAX(nstep, INT(DT/delta_tp + 1.))
3699 if (ksed1(2) .eq. kte) ksed1(2) = kte-1
3700 if (nstep .gt. 0) onstep(2) = 1./REAL(nstep)
3703 !+---+-----------------------------------------------------------------+
3705 if (ANY(L_qs .eqv. .true.)) then
3710 if (rs(k).gt. R1) then
3711 xDs = smoc(k) / smob(k)
3713 ils1 = 1./(Mrat*Lam0 + fv_s)
3714 ils2 = 1./(Mrat*Lam1 + fv_s)
3715 t1_vts = Kap0*csg(4)*ils1**cse(4)
3716 t2_vts = Kap1*Mrat**mu_s*csg(10)*ils2**cse(10)
3717 ils1 = 1./(Mrat*Lam0)
3718 ils2 = 1./(Mrat*Lam1)
3719 t3_vts = Kap0*csg(1)*ils1**cse(1)
3720 t4_vts = Kap1*Mrat**mu_s*csg(7)*ils2**cse(7)
3721 vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts)
3722 if (prr_sml(k) .gt. 0.0) then
3723 SR = rs(k)/(rs(k)+rr(k))
3724 vtsk(k) = vts*SR + (1.-SR)*vtrk(k)
3726 vtsk(k) = vts*vts_boost(k)
3732 if (vtsk(k) .gt. 1.E-3) then
3733 ksed1(3) = MAX(ksed1(3), k)
3734 delta_tp = dzq(k)/vtsk(k)
3735 nstep = MAX(nstep, INT(DT/delta_tp + 1.))
3738 if (ksed1(3) .eq. kte) ksed1(3) = kte-1
3739 if (nstep .gt. 0) onstep(3) = 1./REAL(nstep)
3742 !+---+-----------------------------------------------------------------+
3744 if (ANY(L_qg .eqv. .true.)) then
3749 if (rg(k).gt. R1) then
3750 if (is_hail_aware) then
3751 xrho_g = MAX(rho_g(1),MIN(rg(k)/rho(k)/rb(k),rho_g(NRHG)))
3752 afall = a_coeff*((4.0*xrho_g*9.8)/(3.0*rho(k)))**b_coeff
3753 afall = afall * visco(k)**(1.0-2.0*b_coeff)
3758 vtg = rhof(k)*afall*cgg(6,idx_bg(k))*ogg3 * ilamg(k)**bfall
3760 ! Goal: less prominent size sorting
3761 ! the ELSE section below is technically (mathematically) correct:
3762 if (mu_g .eq. 0) then
3763 vtg = rhof(k)*afall*cgg(7,idx_bg(k))/cgg(12,idx_bg(k)) * ilamg(k)**bfall
3765 vtg = rhof(k)*afall*cgg(8,idx_bg(k))*ogg2 * ilamg(k)**bfall
3770 vtngk(k) = vtngk(k+1)
3773 if (vtgk(k) .gt. 1.E-3) then
3774 ksed1(4) = MAX(ksed1(4), k)
3775 delta_tp = dzq(k)/vtgk(k)
3776 nstep = MAX(nstep, INT(DT/delta_tp + 1.))
3779 if (ksed1(4) .eq. kte) ksed1(4) = kte-1
3780 if (nstep .gt. 0) onstep(4) = 1./REAL(nstep)
3784 !+---+-----------------------------------------------------------------+
3785 !..Sedimentation of mixing ratio is the integral of v(D)*m(D)*N(D)*dD,
3786 !.. whereas neglect m(D) term for number concentration. Therefore,
3787 !.. cloud ice has proper differential sedimentation.
3788 !+---+-----------------------------------------------------------------+
3790 if (ANY(L_qr .eqv. .true.)) then
3791 nstep = NINT(1./onstep(1))
3794 sed_r(k) = vtrk(k)*rr(k)
3795 sed_n(k) = vtnrk(k)*nr(k)
3800 qrten(k) = qrten(k) - sed_r(k)*odzq*onstep(1)*orho
3801 nrten(k) = nrten(k) - sed_n(k)*odzq*onstep(1)*orho
3802 rr(k) = MAX(R1, rr(k) - sed_r(k)*odzq*DT*onstep(1))
3803 nr(k) = MAX(R2, nr(k) - sed_n(k)*odzq*DT*onstep(1))
3804 do k = ksed1(1), kts, -1
3807 qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) &
3808 *odzq*onstep(1)*orho
3809 nrten(k) = nrten(k) + (sed_n(k+1)-sed_n(k)) &
3810 *odzq*onstep(1)*orho
3811 rr(k) = MAX(R1, rr(k) + (sed_r(k+1)-sed_r(k)) &
3813 nr(k) = MAX(R2, nr(k) + (sed_n(k+1)-sed_n(k)) &
3817 if (rr(kts).gt.R1*1000.) &
3818 pptrain = pptrain + sed_r(kts)*DT*onstep(1)
3822 !+---+-----------------------------------------------------------------+
3824 if (ANY(L_qc .eqv. .true.)) then
3826 sed_c(k) = vtck(k)*rc(k)
3827 sed_n(k) = vtnck(k)*nc(k)
3829 do k = ksed1(5), kts, -1
3832 qcten(k) = qcten(k) + (sed_c(k+1)-sed_c(k)) *odzq*orho
3833 ncten(k) = ncten(k) + (sed_n(k+1)-sed_n(k)) *odzq*orho
3834 rc(k) = MAX(R1, rc(k) + (sed_c(k+1)-sed_c(k)) *odzq*DT)
3835 nc(k) = MAX(10., nc(k) + (sed_n(k+1)-sed_n(k)) *odzq*DT)
3839 !+---+-----------------------------------------------------------------+
3841 if (ANY(L_qi .eqv. .true.)) then
3842 nstep = NINT(1./onstep(2))
3845 sed_i(k) = vtik(k)*ri(k)
3846 sed_n(k) = vtnik(k)*ni(k)
3851 qiten(k) = qiten(k) - sed_i(k)*odzq*onstep(2)*orho
3852 niten(k) = niten(k) - sed_n(k)*odzq*onstep(2)*orho
3853 ri(k) = MAX(R1, ri(k) - sed_i(k)*odzq*DT*onstep(2))
3854 ni(k) = MAX(R2, ni(k) - sed_n(k)*odzq*DT*onstep(2))
3855 do k = ksed1(2), kts, -1
3858 qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) &
3859 *odzq*onstep(2)*orho
3860 niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) &
3861 *odzq*onstep(2)*orho
3862 ri(k) = MAX(R1, ri(k) + (sed_i(k+1)-sed_i(k)) &
3864 ni(k) = MAX(R2, ni(k) + (sed_n(k+1)-sed_n(k)) &
3868 if (ri(kts).gt.R1*1000.) &
3869 pptice = pptice + sed_i(kts)*DT*onstep(2)
3873 !+---+-----------------------------------------------------------------+
3875 if (ANY(L_qs .eqv. .true.)) then
3876 nstep = NINT(1./onstep(3))
3879 sed_s(k) = vtsk(k)*rs(k)
3884 qsten(k) = qsten(k) - sed_s(k)*odzq*onstep(3)*orho
3885 rs(k) = MAX(R1, rs(k) - sed_s(k)*odzq*DT*onstep(3))
3886 do k = ksed1(3), kts, -1
3889 qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) &
3890 *odzq*onstep(3)*orho
3891 rs(k) = MAX(R1, rs(k) + (sed_s(k+1)-sed_s(k)) &
3895 if (rs(kts).gt.R1*1000.) &
3896 pptsnow = pptsnow + sed_s(kts)*DT*onstep(3)
3900 !+---+-----------------------------------------------------------------+
3902 if (ANY(L_qg .eqv. .true.)) then
3903 nstep = NINT(1./onstep(4))
3906 sed_g(k) = vtgk(k)*rg(k)
3907 sed_n(k) = vtngk(k)*ng(k)
3908 sed_b(k) = vtgk(k)*rb(k)
3913 qgten(k) = qgten(k) - sed_g(k)*odzq*onstep(4)*orho
3914 ngten(k) = ngten(k) - sed_n(k)*odzq*onstep(4)*orho
3915 qbten(k) = qbten(k) - sed_b(k)*odzq*onstep(4)
3916 rg(k) = MAX(R1, rg(k) - sed_g(k)*odzq*DT*onstep(4))
3917 ng(k) = MAX(R2, ng(k) - sed_n(k)*odzq*DT*onstep(4))
3918 rb(k) = MAX(R1/rho(k)/rho_g(NRHG), rb(k) - sed_b(k)*odzq*DT*onstep(4))
3919 do k = ksed1(4), kts, -1
3922 qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) &
3923 *odzq*onstep(4)*orho
3924 ngten(k) = ngten(k) + (sed_n(k+1)-sed_n(k)) &
3925 *odzq*onstep(4)*orho
3926 qbten(k) = qbten(k) + (sed_b(k+1)-sed_b(k)) &
3928 rg(k) = MAX(R1, rg(k) + (sed_g(k+1)-sed_g(k)) &
3930 ng(k) = MAX(R2, ng(k) + (sed_n(k+1)-sed_n(k)) &
3932 rb(k) = MAX(rg(k)/rho(k)/rho_g(NRHG), rb(k) + (sed_b(k+1)-sed_b(k)) &
3936 if (rg(kts).gt.R1*1000.) &
3937 pptgraul = pptgraul + sed_g(kts)*DT*onstep(4)
3941 !+---+-----------------------------------------------------------------+
3942 !.. Instantly melt any cloud ice into cloud water if above 0C and
3943 !.. instantly freeze any cloud water found below HGFR.
3944 !+---+-----------------------------------------------------------------+
3945 if (.not. iiwarm) then
3947 xri = MAX(0.0, qi1d(k) + qiten(k)*DT)
3948 if ( (temp(k).gt. T_0) .and. (xri.gt. 0.0) ) then
3949 qcten(k) = qcten(k) + xri*odt
3950 ncten(k) = ncten(k) + ni1d(k)*odt
3951 qiten(k) = qiten(k) - xri*odt
3952 niten(k) = -ni1d(k)*odt
3953 tten(k) = tten(k) - lfus*ocp(k)*xri*odt*(1-IFDRY)
3956 xrc = MAX(0.0, qc1d(k) + qcten(k)*DT)
3957 if ( (temp(k).lt. HGFR) .and. (xrc.gt. 0.0) ) then
3958 lfus2 = lsub - lvap(k)
3959 xnc = nc1d(k) + ncten(k)*DT
3960 qiten(k) = qiten(k) + xrc*odt
3961 niten(k) = niten(k) + xnc*odt
3962 qcten(k) = qcten(k) - xrc*odt
3963 ncten(k) = ncten(k) - xnc*odt
3964 tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-IFDRY)
3969 !+---+-----------------------------------------------------------------+
3970 !.. All tendencies computed, apply and pass back final values to parent.
3971 !+---+-----------------------------------------------------------------+
3973 t1d(k) = t1d(k) + tten(k)*DT
3974 qv1d(k) = MAX(1.E-10, qv1d(k) + qvten(k)*DT)
3975 qc1d(k) = qc1d(k) + qcten(k)*DT
3976 nc1d(k) = MAX(2./rho(k), MIN(nc1d(k) + ncten(k)*DT, Nt_c_max))
3977 if (is_aerosol_aware) then
3978 if (aer_init_opt .lt. 2) then
3979 nwfa1d(k) = MAX(11.1E6, MIN(9999.E6, &
3980 (nwfa1d(k)+nwfaten(k)*DT)))
3981 nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6, &
3982 (nifa1d(k)+nifaten(k)*DT)))
3983 if (wif_input_opt .eq. 2) then
3984 nbca1d(k) = MAX(5.55E6, MIN(9999.E6, &
3985 (nbca1d(k)+nbcaten(k)*DT)))
3990 nwfa1d(k) = MAX(0.0, (nwfa1d(k)+nwfaten(k)*DT))
3991 nifa1d(k) = MAX(0.0, (nifa1d(k)+nifaten(k)*DT))
3992 if (wif_input_opt .eq. 2) then
3993 nbca1d(k) = MAX(0.0, (nbca1d(k)+nbcaten(k)*DT))
3999 nwfa1d(k) = MAX(11.1E6, MIN(9999.E6, &
4000 (nwfa1d(k)+nwfaten(k)*DT)))
4001 nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6, &
4002 (nifa1d(k)+nifaten(k)*DT)))
4003 nbca1d(k) = MAX(5.55E6, MIN(9999.E6, &
4004 (nbca1d(k)+nbcaten(k)*DT)))
4007 if (qc1d(k) .le. R1) then
4011 nu_c = MIN(15, NINT(1000.E6/(nc1d(k)*rho(k))) + 2)
4012 lamc = (am_r*ccg(2,nu_c)*ocg1(nu_c)*nc1d(k)/qc1d(k))**obmr
4013 xDc = (bm_r + nu_c + 1.) / lamc
4014 if (xDc.lt. D0c) then
4015 lamc = cce(2,nu_c)/D0c
4016 elseif (xDc.gt. D0r*2.) then
4017 lamc = cce(2,nu_c)/(D0r*2.)
4019 nc1d(k) = MIN(ccg(1,nu_c)*ocg2(nu_c)*qc1d(k)/am_r*lamc**bm_r,&
4020 DBLE(Nt_c_max)/rho(k))
4023 qi1d(k) = qi1d(k) + qiten(k)*DT
4024 ni1d(k) = MAX(R2/rho(k), ni1d(k) + niten(k)*DT)
4025 if (qi1d(k) .le. R1) then
4029 lami = (am_i*cig(2)*oig1*ni1d(k)/qi1d(k))**obmi
4031 xDi = (bm_i + mu_i + 1.) * ilami
4032 if (xDi.lt. 5.E-6) then
4034 elseif (xDi.gt. 300.E-6) then
4035 lami = cie(2)/300.E-6
4037 ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, &
4040 qr1d(k) = qr1d(k) + qrten(k)*DT
4041 nr1d(k) = MAX(R2/rho(k), nr1d(k) + nrten(k)*DT)
4042 if (qr1d(k) .le. R1) then
4046 lamr = (am_r*crg(3)*org2*nr1d(k)/qr1d(k))**obmr
4047 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
4048 if (mvd_r(k) .gt. 2.5E-3) then
4050 elseif (mvd_r(k) .lt. D0r*0.75) then
4053 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
4054 nr1d(k) = crg(2)*org3*qr1d(k)*lamr**bm_r / am_r
4056 qs1d(k) = qs1d(k) + qsten(k)*DT
4057 if (qs1d(k) .le. R1) qs1d(k) = 0.0
4058 qg1d(k) = qg1d(k) + qgten(k)*DT
4059 ng1d(k) = MAX(R2/rho(k), ng1d(k) + ngten(k)*DT)
4060 if (qg1d(k) .le. R1) then
4065 qb1d(k) = MAX(qg1d(k)/rho_g(NRHG), qb1d(k) + qbten(k)*DT)
4066 qb1d(k) = MIN(qg1d(k)/rho_g(1), qb1d(k))
4067 idx_bg(k) = MAX(1,MIN(NINT(qg1d(k)/qb1d(k) *0.01)+1,NRHG))
4068 if (.not. is_hail_aware) idx_bg(k) = idx_bg1
4069 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng1d(k)/qg1d(k))**obmg
4070 mvd_g(k) = (3.0 + mu_g + 0.672) / lamg
4071 if (mvd_g(k) .gt. 25.4E-3) then
4073 elseif (mvd_g(k) .lt. D0r) then
4076 lamg = (3.0 + mu_g + 0.672) / mvd_g(k)
4077 ng1d(k) = cgg(2,1)*ogg3*qg1d(k)*lamg**bm_g / am_g(idx_bg(k))
4082 end subroutine mp_thompson
4083 !+---+-----------------------------------------------------------------+
4085 !+---+-----------------------------------------------------------------+
4086 !..Creation of the lookup tables and support functions found below here.
4087 !+---+-----------------------------------------------------------------+
4088 !..Rain collecting graupel (and inverse). Explicit CE integration.
4089 !+---+-----------------------------------------------------------------+
4091 subroutine qr_acr_qg(NRHGtable)
4096 INTEGER, INTENT(IN) ::NRHGtable
4099 INTEGER:: i, j, k, m, n, n2, n3, idx_bg
4100 INTEGER:: km, km_s, km_e
4101 DOUBLE PRECISION, DIMENSION(nbg):: N_g
4102 DOUBLE PRECISION, DIMENSION(nbg,NRHGtable):: vg
4103 DOUBLE PRECISION, DIMENSION(nbr):: vr, N_r
4104 DOUBLE PRECISION:: N0_r, N0_g, lam_exp, lamg, lamr
4105 DOUBLE PRECISION:: massg, massr, dvg, dvr, t1, t2, z1, z2, y1, y2
4106 LOGICAL force_read_thompson, write_thompson_tables, write_thompson_mp38table
4107 LOGICAL lexist,lopen
4109 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
4110 CHARACTER (LEN=20):: qr_acr_qg_name
4113 CALL nl_get_force_read_thompson(1,force_read_thompson)
4114 CALL nl_get_write_thompson_tables(1,write_thompson_tables)
4115 CALL nl_get_write_thompson_mp38table(1,write_thompson_mp38table)
4117 if (is_hail_aware) then
4119 ! V1: new table dimension (NRHG), computes all, rho_g(1:NRHG)
4120 qr_acr_qg_name="qr_acr_qg_mp38V1.dat"
4122 ! Table for mp=8 or mp=28
4123 ! V4: new table dimension (NRHG), computes only rho_g(idx_bg1)
4124 qr_acr_qg_name="qr_acr_qg_V4.dat"
4128 IF ( wrf_dm_on_monitor() ) THEN
4129 INQUIRE(FILE=qr_acr_qg_name,EXIST=lexist)
4131 CALL wrf_message("ThompMP: read "//qr_acr_qg_name//" instead of computing")
4132 OPEN(63,file=qr_acr_qg_name,form="unformatted",err=1234)
4133 READ(63,err=1234) tcg_racg
4134 READ(63,err=1234) tmr_racg
4135 READ(63,err=1234) tcr_gacr
4136 !READ(63,err=1234) tmg_gacr
4137 READ(63,err=1234) tnr_racg
4138 READ(63,err=1234) tnr_gacr
4141 IF ( good .NE. 1 ) THEN
4142 INQUIRE(63,opened=lopen)
4144 IF( force_read_thompson ) THEN
4145 CALL wrf_error_fatal("Error reading "//qr_acr_qg_name//". Aborting because force_read_thompson is .true.")
4149 IF( force_read_thompson ) THEN
4150 CALL wrf_error_fatal("Error opening "//qr_acr_qg_name//". Aborting because force_read_thompson is .true.")
4154 INQUIRE(63,OPENED=lopen)
4160 IF( force_read_thompson ) THEN
4161 CALL wrf_error_fatal("Non-existent "//qr_acr_qg_name//". Aborting because force_read_thompson is .true.")
4163 IF( (is_hail_aware) .AND. (.NOT. write_thompson_mp38table)) THEN
4164 CALL wrf_error_fatal("Non-existent "//qr_acr_qg_name//". Aborting because write_thompson_mp38table is .false.")
4168 #if defined(DM_PARALLEL) && !defined(STUBMPI)
4169 CALL wrf_dm_bcast_integer(good,1)
4172 IF ( good .EQ. 1 ) THEN
4173 #if defined(DM_PARALLEL) && !defined(STUBMPI)
4174 CALL wrf_dm_bcast_double(tcg_racg,SIZE(tcg_racg))
4175 CALL wrf_dm_bcast_double(tmr_racg,SIZE(tmr_racg))
4176 CALL wrf_dm_bcast_double(tcr_gacr,SIZE(tcr_gacr))
4177 !CALL wrf_dm_bcast_double(tmg_gacr,SIZE(tmg_gacr))
4178 CALL wrf_dm_bcast_double(tnr_racg,SIZE(tnr_racg))
4179 CALL wrf_dm_bcast_double(tnr_gacr,SIZE(tnr_gacr))
4182 CALL wrf_message("ThompMP: computing qr_acr_qg table: "//qr_acr_qg_name)
4183 !+---+-----------------------------------------------------------------+
4187 ! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
4188 vr(n2) = -0.1021 + 4.932E3*Dr(n2) - 0.9551E6*Dr(n2)*Dr(n2) &
4189 + 0.07934E9*Dr(n2)*Dr(n2)*Dr(n2) &
4190 - 0.002362E12*Dr(n2)*Dr(n2)*Dr(n2)*Dr(n2)
4192 do n3 = 1, NRHGtable
4194 ! idx_bg indexes module coefficients, not vg
4195 if (.not. is_hail_aware) idx_bg = idx_bg1
4196 if (is_hail_aware) idx_bg = n3
4197 vg(n,n3) = av_g(idx_bg)*Dg(n)**bv_g(idx_bg)
4201 !..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for
4202 !.. fortran indices. J. Michalakes, 2009Oct30.
4204 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
4205 CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e )
4208 km_e = ntb_r*ntb_r1 - 1
4213 k = mod( km , ntb_r1 ) + 1
4215 lam_exp = (N0r_exp(k)*am_r*crg(1)/r_r(m))**ore1
4216 lamr = lam_exp * (crg(3)*org2*org1)**obmr
4217 N0_r = N0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2)
4219 N_r(n2) = N0_r*Dr(n2)**mu_r *DEXP(-lamr*Dr(n2))*dtr(n2)
4222 do n3 = 1, NRHGtable
4223 if (.not. is_hail_aware) idx_bg = idx_bg1
4224 if (is_hail_aware) idx_bg = n3
4228 lam_exp = (N0g_exp(i)*am_g(idx_bg)*cgg(1,1)/r_g(j))**oge1
4229 lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg
4230 N0_g = N0g_exp(i)/(cgg(2,1)*lam_exp) * lamg**cge(2,1)
4232 N_g(n) = N0_g*Dg(n)**mu_g * DEXP(-lamg*Dg(n))*dtg(n)
4242 massr = am_r * Dr(n2)**bm_r
4244 massg = am_g(idx_bg) * Dg(n)**bm_g
4246 dvg = 0.5d0*((vr(n2) - vg(n,n3)) + DABS(vr(n2)-vg(n,n3)))
4247 dvr = 0.5d0*((vg(n,n3) - vr(n2)) + DABS(vg(n,n3)-vr(n2)))
4249 t1 = t1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
4250 *dvg*massg * N_g(n)* N_r(n2)
4251 z1 = z1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
4252 *dvg*massr * N_g(n)* N_r(n2)
4253 y1 = y1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
4254 *dvg * N_g(n)* N_r(n2)
4256 t2 = t2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
4257 *dvr*massr * N_g(n)* N_r(n2)
4258 y2 = y2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
4259 *dvr * N_g(n)* N_r(n2)
4260 z2 = z2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
4261 *dvr*massg * N_g(n)* N_r(n2)
4265 tcg_racg(i,j,n3,k,m) = t1
4266 tmr_racg(i,j,n3,k,m) = DMIN1(z1, r_r(m)*1.0d0)
4267 tcr_gacr(i,j,n3,k,m) = t2
4268 !tmg_gacr(i,j,n3,k,m) = DMIN1(z2, r_g(j)*1.0d0)
4269 tnr_racg(i,j,n3,k,m) = y1
4270 tnr_gacr(i,j,n3,k,m) = y2
4276 !..Note wrf_dm_gatherv expects zero-based km_s, km_e (J. Michalakes, 2009Oct30).
4278 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
4279 CALL wrf_dm_gatherv(tcg_racg, ntb_g*ntb_g1*NRHGtable, km_s, km_e, R8SIZE)
4280 CALL wrf_dm_gatherv(tmr_racg, ntb_g*ntb_g1*NRHGtable, km_s, km_e, R8SIZE)
4281 CALL wrf_dm_gatherv(tcr_gacr, ntb_g*ntb_g1*NRHGtable, km_s, km_e, R8SIZE)
4282 !CALL wrf_dm_gatherv(tmg_gacr, ntb_g*ntb_g1*NRHGtable, km_s, km_e, R8SIZE)
4283 CALL wrf_dm_gatherv(tnr_racg, ntb_g*ntb_g1*NRHGtable, km_s, km_e, R8SIZE)
4284 CALL wrf_dm_gatherv(tnr_gacr, ntb_g*ntb_g1*NRHGtable, km_s, km_e, R8SIZE)
4286 CALL end_timing ( TRIM("table computation "//qr_acr_qg_name//"") )
4287 !+---+-----------------------------------------------------------------+
4290 IF ( write_thompson_tables .AND. wrf_dm_on_monitor() ) THEN
4291 CALL wrf_message("Writing "//qr_acr_qg_name//" in Thompson MP init")
4292 OPEN(63,file=qr_acr_qg_name,form="unformatted",err=9234)
4293 WRITE(63,err=9234) tcg_racg
4294 WRITE(63,err=9234) tmr_racg
4295 WRITE(63,err=9234) tcr_gacr
4296 !WRITE(63,err=9234) tmg_gacr
4297 WRITE(63,err=9234) tnr_racg
4298 WRITE(63,err=9234) tnr_gacr
4300 RETURN ! ----- RETURN
4302 CALL wrf_error_fatal("Error writing "//qr_acr_qg_name//"")
4306 end subroutine qr_acr_qg
4307 !+---+-----------------------------------------------------------------+
4309 !+---+-----------------------------------------------------------------+
4310 !..Rain collecting snow (and inverse). Explicit CE integration.
4311 !+---+-----------------------------------------------------------------+
4313 subroutine qr_acr_qs
4318 INTEGER:: i, j, k, m, n, n2
4319 INTEGER:: km, km_s, km_e
4320 DOUBLE PRECISION, DIMENSION(nbr):: vr, D1, N_r
4321 DOUBLE PRECISION, DIMENSION(nbs):: vs, N_s
4322 DOUBLE PRECISION:: loga_, a_, b_, second, M0, M2, M3, Mrat, oM3
4323 DOUBLE PRECISION:: N0_r, lam_exp, lamr, slam1, slam2
4324 DOUBLE PRECISION:: dvs, dvr, masss, massr
4325 DOUBLE PRECISION:: t1, t2, t3, t4, z1, z2, z3, z4
4326 DOUBLE PRECISION:: y1, y2, y3, y4
4327 LOGICAL force_read_thompson, write_thompson_tables
4328 LOGICAL lexist,lopen
4330 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
4334 CALL nl_get_force_read_thompson(1,force_read_thompson)
4335 CALL nl_get_write_thompson_tables(1,write_thompson_tables)
4338 IF ( wrf_dm_on_monitor() ) THEN
4339 INQUIRE(FILE="qr_acr_qsV2.dat",EXIST=lexist)
4341 CALL wrf_message("ThompMP: read qr_acr_qsV2.dat instead of computing")
4342 OPEN(63,file="qr_acr_qsV2.dat",form="unformatted",err=1234)
4343 READ(63,err=1234)tcs_racs1
4344 READ(63,err=1234)tmr_racs1
4345 READ(63,err=1234)tcs_racs2
4346 READ(63,err=1234)tmr_racs2
4347 READ(63,err=1234)tcr_sacr1
4348 READ(63,err=1234)tms_sacr1
4349 READ(63,err=1234)tcr_sacr2
4350 READ(63,err=1234)tms_sacr2
4351 READ(63,err=1234)tnr_racs1
4352 READ(63,err=1234)tnr_racs2
4353 READ(63,err=1234)tnr_sacr1
4354 READ(63,err=1234)tnr_sacr2
4357 IF ( good .NE. 1 ) THEN
4358 INQUIRE(63,opened=lopen)
4360 IF( force_read_thompson ) THEN
4361 CALL wrf_error_fatal("Error reading qr_acr_qsV2.dat. Aborting because force_read_thompson is .true.")
4365 IF( force_read_thompson ) THEN
4366 CALL wrf_error_fatal("Error opening qr_acr_qsV2.dat. Aborting because force_read_thompson is .true.")
4370 INQUIRE(63,opened=lopen)
4376 IF( force_read_thompson ) THEN
4377 CALL wrf_error_fatal("Non-existent qr_acr_qsV2.dat. Aborting because force_read_thompson is .true.")
4381 #if defined(DM_PARALLEL) && !defined(STUBMPI)
4382 CALL wrf_dm_bcast_integer(good,1)
4385 IF ( good .EQ. 1 ) THEN
4386 #if defined(DM_PARALLEL) && !defined(STUBMPI)
4387 CALL wrf_dm_bcast_double(tcs_racs1,SIZE(tcs_racs1))
4388 CALL wrf_dm_bcast_double(tmr_racs1,SIZE(tmr_racs1))
4389 CALL wrf_dm_bcast_double(tcs_racs2,SIZE(tcs_racs2))
4390 CALL wrf_dm_bcast_double(tmr_racs2,SIZE(tmr_racs2))
4391 CALL wrf_dm_bcast_double(tcr_sacr1,SIZE(tcr_sacr1))
4392 CALL wrf_dm_bcast_double(tms_sacr1,SIZE(tms_sacr1))
4393 CALL wrf_dm_bcast_double(tcr_sacr2,SIZE(tcr_sacr2))
4394 CALL wrf_dm_bcast_double(tms_sacr2,SIZE(tms_sacr2))
4395 CALL wrf_dm_bcast_double(tnr_racs1,SIZE(tnr_racs1))
4396 CALL wrf_dm_bcast_double(tnr_racs2,SIZE(tnr_racs2))
4397 CALL wrf_dm_bcast_double(tnr_sacr1,SIZE(tnr_sacr1))
4398 CALL wrf_dm_bcast_double(tnr_sacr2,SIZE(tnr_sacr2))
4401 CALL wrf_message("ThompMP: computing qr_acr_qs")
4403 ! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
4404 vr(n2) = -0.1021 + 4.932E3*Dr(n2) - 0.9551E6*Dr(n2)*Dr(n2) &
4405 + 0.07934E9*Dr(n2)*Dr(n2)*Dr(n2) &
4406 - 0.002362E12*Dr(n2)*Dr(n2)*Dr(n2)*Dr(n2)
4407 D1(n2) = (vr(n2)/av_s)**(1./bv_s)
4410 vs(n) = 1.5*av_s*Ds(n)**bv_s * DEXP(-fv_s*Ds(n))
4413 !..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for
4414 !.. fortran indices. J. Michalakes, 2009Oct30.
4416 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
4417 CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e )
4420 km_e = ntb_r*ntb_r1 - 1
4425 k = mod( km , ntb_r1 ) + 1
4427 lam_exp = (N0r_exp(k)*am_r*crg(1)/r_r(m))**ore1
4428 lamr = lam_exp * (crg(3)*org2*org1)**obmr
4429 N0_r = N0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2)
4431 N_r(n2) = N0_r*Dr(n2)**mu_r * DEXP(-lamr*Dr(n2))*dtr(n2)
4437 !..From the bm_s moment, compute plus one moment. If we are not
4438 !.. using bm_s=2, then we must transform to the pure 2nd moment
4439 !.. (variable called "second") and then to the bm_s+1 moment.
4441 M2 = r_s(i)*oams *1.0d0
4442 if (bm_s.gt.2.0-1.E-3 .and. bm_s.lt.2.0+1.E-3) then
4443 loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*bm_s &
4444 + sa(4)*Tc(j)*bm_s + sa(5)*Tc(j)*Tc(j) &
4445 + sa(6)*bm_s*bm_s + sa(7)*Tc(j)*Tc(j)*bm_s &
4446 + sa(8)*Tc(j)*bm_s*bm_s + sa(9)*Tc(j)*Tc(j)*Tc(j) &
4447 + sa(10)*bm_s*bm_s*bm_s
4449 b_ = sb(1) + sb(2)*Tc(j) + sb(3)*bm_s &
4450 + sb(4)*Tc(j)*bm_s + sb(5)*Tc(j)*Tc(j) &
4451 + sb(6)*bm_s*bm_s + sb(7)*Tc(j)*Tc(j)*bm_s &
4452 + sb(8)*Tc(j)*bm_s*bm_s + sb(9)*Tc(j)*Tc(j)*Tc(j) &
4453 + sb(10)*bm_s*bm_s*bm_s
4454 second = (M2/a_)**(1./b_)
4459 loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*cse(1) &
4460 + sa(4)*Tc(j)*cse(1) + sa(5)*Tc(j)*Tc(j) &
4461 + sa(6)*cse(1)*cse(1) + sa(7)*Tc(j)*Tc(j)*cse(1) &
4462 + sa(8)*Tc(j)*cse(1)*cse(1) + sa(9)*Tc(j)*Tc(j)*Tc(j) &
4463 + sa(10)*cse(1)*cse(1)*cse(1)
4465 b_ = sb(1)+sb(2)*Tc(j)+sb(3)*cse(1) + sb(4)*Tc(j)*cse(1) &
4466 + sb(5)*Tc(j)*Tc(j) + sb(6)*cse(1)*cse(1) &
4467 + sb(7)*Tc(j)*Tc(j)*cse(1) + sb(8)*Tc(j)*cse(1)*cse(1) &
4468 + sb(9)*Tc(j)*Tc(j)*Tc(j)+sb(10)*cse(1)*cse(1)*cse(1)
4469 M3 = a_ * second**b_
4472 Mrat = M2*(M2*oM3)*(M2*oM3)*(M2*oM3)
4474 slam1 = M2 * oM3 * Lam0
4475 slam2 = M2 * oM3 * Lam1
4478 N_s(n) = Mrat*(Kap0*DEXP(-slam1*Ds(n)) &
4479 + Kap1*M0*Ds(n)**mu_s * DEXP(-slam2*Ds(n)))*dts(n)
4495 massr = am_r * Dr(n2)**bm_r
4497 masss = am_s * Ds(n)**bm_s
4499 dvs = 0.5d0*((vr(n2) - vs(n)) + DABS(vr(n2)-vs(n)))
4500 dvr = 0.5d0*((vs(n) - vr(n2)) + DABS(vs(n)-vr(n2)))
4502 if (massr .gt. 1.5*masss) then
4503 t1 = t1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4504 *dvs*masss * N_s(n)* N_r(n2)
4505 z1 = z1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4506 *dvs*massr * N_s(n)* N_r(n2)
4507 y1 = y1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4508 *dvs * N_s(n)* N_r(n2)
4510 t3 = t3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4511 *dvs*masss * N_s(n)* N_r(n2)
4512 z3 = z3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4513 *dvs*massr * N_s(n)* N_r(n2)
4514 y3 = y3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4515 *dvs * N_s(n)* N_r(n2)
4518 if (massr .gt. 1.5*masss) then
4519 t2 = t2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4520 *dvr*massr * N_s(n)* N_r(n2)
4521 y2 = y2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4522 *dvr * N_s(n)* N_r(n2)
4523 z2 = z2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4524 *dvr*masss * N_s(n)* N_r(n2)
4526 t4 = t4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4527 *dvr*massr * N_s(n)* N_r(n2)
4528 y4 = y4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4529 *dvr * N_s(n)* N_r(n2)
4530 z4 = z4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4531 *dvr*masss * N_s(n)* N_r(n2)
4536 tcs_racs1(i,j,k,m) = t1
4537 tmr_racs1(i,j,k,m) = DMIN1(z1, r_r(m)*1.0d0)
4538 tcs_racs2(i,j,k,m) = t3
4539 tmr_racs2(i,j,k,m) = z3
4540 tcr_sacr1(i,j,k,m) = t2
4541 tms_sacr1(i,j,k,m) = z2
4542 tcr_sacr2(i,j,k,m) = t4
4543 tms_sacr2(i,j,k,m) = z4
4544 tnr_racs1(i,j,k,m) = y1
4545 tnr_racs2(i,j,k,m) = y3
4546 tnr_sacr1(i,j,k,m) = y2
4547 tnr_sacr2(i,j,k,m) = y4
4552 !..Note wrf_dm_gatherv expects zero-based km_s, km_e (J. Michalakes, 2009Oct30).
4554 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
4555 CALL wrf_dm_gatherv(tcs_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4556 CALL wrf_dm_gatherv(tmr_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4557 CALL wrf_dm_gatherv(tcs_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4558 CALL wrf_dm_gatherv(tmr_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4559 CALL wrf_dm_gatherv(tcr_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4560 CALL wrf_dm_gatherv(tms_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4561 CALL wrf_dm_gatherv(tcr_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4562 CALL wrf_dm_gatherv(tms_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4563 CALL wrf_dm_gatherv(tnr_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4564 CALL wrf_dm_gatherv(tnr_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4565 CALL wrf_dm_gatherv(tnr_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4566 CALL wrf_dm_gatherv(tnr_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4569 IF ( write_thompson_tables .AND. wrf_dm_on_monitor() ) THEN
4570 CALL wrf_message("Writing qr_acr_qsV2.dat in Thompson MP init")
4571 OPEN(63,file="qr_acr_qsV2.dat",form="unformatted",err=9234)
4572 WRITE(63,err=9234)tcs_racs1
4573 WRITE(63,err=9234)tmr_racs1
4574 WRITE(63,err=9234)tcs_racs2
4575 WRITE(63,err=9234)tmr_racs2
4576 WRITE(63,err=9234)tcr_sacr1
4577 WRITE(63,err=9234)tms_sacr1
4578 WRITE(63,err=9234)tcr_sacr2
4579 WRITE(63,err=9234)tms_sacr2
4580 WRITE(63,err=9234)tnr_racs1
4581 WRITE(63,err=9234)tnr_racs2
4582 WRITE(63,err=9234)tnr_sacr1
4583 WRITE(63,err=9234)tnr_sacr2
4585 RETURN ! ----- RETURN
4587 CALL wrf_error_fatal("Error writing qr_acr_qsV2.dat")
4591 end subroutine qr_acr_qs
4592 !+---+-----------------------------------------------------------------+
4594 !+---+-----------------------------------------------------------------+
4595 !..This is a literal adaptation of Bigg (1954) probability of drops of
4596 !..a particular volume freezing. Given this probability, simply freeze
4597 !..the proportion of drops summing their masses.
4598 !+---+-----------------------------------------------------------------+
4600 subroutine freezeH2O
4605 INTEGER:: i, j, k, m, n, n2
4606 INTEGER:: km, km_s, km_e
4607 DOUBLE PRECISION:: N_r, N_c
4608 DOUBLE PRECISION, DIMENSION(nbr):: massr
4609 DOUBLE PRECISION, DIMENSION(nbc):: massc
4610 DOUBLE PRECISION:: sum1, sum2, sumn1, sumn2, &
4611 prob, vol, Texp, orho_w, &
4612 lam_exp, lamr, N0_r, lamc, N0_c, y
4615 LOGICAL force_read_thompson, write_thompson_tables
4616 LOGICAL lexist,lopen
4618 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
4621 CALL nl_get_force_read_thompson(1,force_read_thompson)
4622 CALL nl_get_write_thompson_tables(1,write_thompson_tables)
4625 IF ( wrf_dm_on_monitor() ) THEN
4626 INQUIRE(FILE="freezeH2O.dat",EXIST=lexist)
4628 CALL wrf_message("ThompMP: read freezeH2O.dat instead of computing")
4629 OPEN(63,file="freezeH2O.dat",form="unformatted",err=1234)
4630 READ(63,err=1234)tpi_qrfz
4631 READ(63,err=1234)tni_qrfz
4632 READ(63,err=1234)tpg_qrfz
4633 READ(63,err=1234)tnr_qrfz
4634 READ(63,err=1234)tpi_qcfz
4635 READ(63,err=1234)tni_qcfz
4638 IF ( good .NE. 1 ) THEN
4639 INQUIRE(63,opened=lopen)
4641 IF( force_read_thompson ) THEN
4642 CALL wrf_error_fatal("Error reading freezeH2O.dat. Aborting because force_read_thompson is .true.")
4646 IF( force_read_thompson ) THEN
4647 CALL wrf_error_fatal("Error opening freezeH2O.dat. Aborting because force_read_thompson is .true.")
4651 INQUIRE(63,opened=lopen)
4657 IF( force_read_thompson ) THEN
4658 CALL wrf_error_fatal("Non-existent freezeH2O.dat. Aborting because force_read_thompson is .true.")
4663 #if defined(DM_PARALLEL) && !defined(STUBMPI)
4664 CALL wrf_dm_bcast_integer(good,1)
4667 IF ( good .EQ. 1 ) THEN
4668 #if defined(DM_PARALLEL) && !defined(STUBMPI)
4669 CALL wrf_dm_bcast_double(tpi_qrfz,SIZE(tpi_qrfz))
4670 CALL wrf_dm_bcast_double(tni_qrfz,SIZE(tni_qrfz))
4671 CALL wrf_dm_bcast_double(tpg_qrfz,SIZE(tpg_qrfz))
4672 CALL wrf_dm_bcast_double(tnr_qrfz,SIZE(tnr_qrfz))
4673 CALL wrf_dm_bcast_double(tpi_qcfz,SIZE(tpi_qcfz))
4674 CALL wrf_dm_bcast_double(tni_qcfz,SIZE(tni_qcfz))
4677 CALL wrf_message("ThompMP: computing freezeH2O")
4682 massr(n2) = am_r*Dr(n2)**bm_r
4685 massc(n) = am_r*Dc(n)**bm_r
4688 ! Need to split loops between MPI processes to speedup
4689 ! (2017Jul26, Jason Do)
4690 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
4691 CALL wrf_dm_decomp1d ( ntb_IN*45, km_s, km_e )
4694 km_e = ntb_IN*45 - 1
4697 !..Freeze water (smallest drops become cloud ice, otherwise graupel).
4700 k = mod( km , 45 ) + 1
4701 T_adjust = MAX(-3.0, MIN(3.0 - ALOG10(Nt_IN(m)), 3.0))
4702 ! print*, ' Freezing water for temp = ', -k
4703 Texp = DEXP( DFLOAT(k) - T_adjust*1.0D0 ) - 1.0D0
4705 !$OMP PARALLEL DO SCHEDULE(dynamic) &
4706 !$OMP PRIVATE(j,i,lam_exp,lamr,N0_r,sum1,sum2,sumn1,sumn2,n2,N_r,vol,prob)
4710 lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(i))**ore1
4711 lamr = lam_exp * (crg(3)*org2*org1)**obmr
4712 N0_r = N0r_exp(j)/(crg(2)*lam_exp) * lamr**cre(2)
4718 N_r = N0_r*Dr(n2)**mu_r*DEXP(-lamr*Dr(n2))*dtr(n2)
4719 vol = massr(n2)*orho_w
4720 prob = MAX(0.0D0, 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp))
4721 if (massr(n2) .lt. xm0g) then
4722 sumn1 = sumn1 + prob*N_r
4723 sum1 = sum1 + prob*N_r*massr(n2)
4725 sumn2 = sumn2 + prob*N_r
4726 sum2 = sum2 + prob*N_r*massr(n2)
4728 if ((sum1+sum2).ge.r_r(i)) EXIT
4730 tpi_qrfz(i,j,k,m) = sum1
4731 tni_qrfz(i,j,k,m) = sumn1
4732 tpg_qrfz(i,j,k,m) = sum2
4733 tnr_qrfz(i,j,k,m) = sumn2
4737 !$OMP END PARALLEL DO
4739 !$OMP PARALLEL DO SCHEDULE(dynamic) &
4740 !$OMP PRIVATE(j,i,nu_c,lamc,N0_c,sum1,sumn2,vol,prob,N_c)
4743 nu_c = MIN(15, NINT(1000.E6/t_Nc(j)) + 2)
4745 lamc = (t_Nc(j)*am_r* ccg(2,nu_c) * ocg1(nu_c) / r_c(i))**obmr
4746 N0_c = t_Nc(j)*ocg1(nu_c) * lamc**cce(1,nu_c)
4750 vol = massc(n)*orho_w
4751 prob = MAX(0.0D0, 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp))
4752 N_c = N0_c*Dc(n)**nu_c*EXP(-lamc*Dc(n))*dtc(n)
4753 sumn2 = MIN(t_Nc(j), sumn2 + prob*N_c)
4754 sum1 = sum1 + prob*N_c*massc(n)
4755 if (sum1 .ge. r_c(i)) EXIT
4757 tpi_qcfz(i,j,k,m) = sum1
4758 tni_qcfz(i,j,k,m) = sumn2
4762 !$OMP END PARALLEL DO
4766 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
4767 CALL wrf_dm_gatherv(tpi_qrfz, ntb_r*ntb_r1, km_s, km_e, R8SIZE)
4768 CALL wrf_dm_gatherv(tni_qrfz, ntb_r*ntb_r1, km_s, km_e, R8SIZE)
4769 CALL wrf_dm_gatherv(tpg_qrfz, ntb_r*ntb_r1, km_s, km_e, R8SIZE)
4770 CALL wrf_dm_gatherv(tnr_qrfz, ntb_r*ntb_r1, km_s, km_e, R8SIZE)
4771 CALL wrf_dm_gatherv(tpi_qcfz, ntb_c*nbc, km_s, km_e, R8SIZE)
4772 CALL wrf_dm_gatherv(tni_qcfz, ntb_c*nbc, km_s, km_e, R8SIZE)
4775 IF ( write_thompson_tables .AND. wrf_dm_on_monitor() ) THEN
4776 CALL wrf_message("Writing freezeH2O.dat in Thompson MP init")
4777 OPEN(63,file="freezeH2O.dat",form="unformatted",err=9234)
4778 WRITE(63,err=9234)tpi_qrfz
4779 WRITE(63,err=9234)tni_qrfz
4780 WRITE(63,err=9234)tpg_qrfz
4781 WRITE(63,err=9234)tnr_qrfz
4782 WRITE(63,err=9234)tpi_qcfz
4783 WRITE(63,err=9234)tni_qcfz
4785 RETURN ! ----- RETURN
4787 CALL wrf_error_fatal("Error writing freezeH2O.dat")
4791 end subroutine freezeH2O
4793 !+---+-----------------------------------------------------------------+
4795 !+---+-----------------------------------------------------------------+
4796 !..Cloud ice converting to snow since portion greater than min snow
4797 !.. size. Given cloud ice content (kg/m**3), number concentration
4798 !.. (#/m**3) and gamma shape parameter, mu_i, break the distrib into
4799 !.. bins and figure out the mass/number of ice with sizes larger than
4800 !.. D0s. Also, compute incomplete gamma function for the integration
4801 !.. of ice depositional growth from diameter=0 to D0s. Amount of
4802 !.. ice depositional growth is this portion of distrib while larger
4803 !.. diameters contribute to snow growth (as in Harrington et al. 1995).
4804 !+---+-----------------------------------------------------------------+
4806 subroutine qi_aut_qs
4812 DOUBLE PRECISION, DIMENSION(nbi):: N_i
4813 DOUBLE PRECISION:: N0_i, lami, Di_mean, t1, t2
4820 lami = (am_i*cig(2)*oig1*Nt_i(j)/r_i(i))**obmi
4821 Di_mean = (bm_i + mu_i + 1.) / lami
4822 N0_i = Nt_i(j)*oig1 * lami**cie(1)
4825 if (SNGL(Di_mean) .gt. 5.*D0s) then
4828 tpi_ide(i,j) = 0.0D0
4829 elseif (SNGL(Di_mean) .lt. D0i) then
4832 tpi_ide(i,j) = 1.0D0
4834 xlimit_intg = lami*D0s
4835 tpi_ide(i,j) = GAMMP(mu_i+2.0, xlimit_intg) * 1.0D0
4837 N_i(n2) = N0_i*Di(n2)**mu_i * DEXP(-lami*Di(n2))*dti(n2)
4838 if (Di(n2).ge.D0s) then
4839 t1 = t1 + N_i(n2) * am_i*Di(n2)**bm_i
4849 end subroutine qi_aut_qs
4851 !+---+-----------------------------------------------------------------+
4852 !..Variable collision efficiency for rain collecting cloud water using
4853 !.. method of Beard and Grover, 1974 if a/A less than 0.25; otherwise
4854 !.. uses polynomials to get close match of Pruppacher & Klett Fig 14-9.
4855 !+---+-----------------------------------------------------------------+
4857 subroutine table_Efrw
4862 DOUBLE PRECISION:: vtr, stokes, reynolds, Ef_rw
4863 DOUBLE PRECISION:: p, yc0, F, G, H, z, K0, X
4870 if (Dr(i).lt.50.E-6 .or. Dc(j).lt.3.E-6) then
4872 elseif (p.gt.0.25) then
4874 if (Dr(i) .lt. 75.e-6) then
4875 Ef_rw = 0.026794*X - 0.20604
4876 elseif (Dr(i) .lt. 125.e-6) then
4877 Ef_rw = -0.00066842*X*X + 0.061542*X - 0.37089
4878 elseif (Dr(i) .lt. 175.e-6) then
4879 Ef_rw = 4.091e-06*X*X*X*X - 0.00030908*X*X*X &
4880 + 0.0066237*X*X - 0.0013687*X - 0.073022
4881 elseif (Dr(i) .lt. 250.e-6) then
4882 Ef_rw = 9.6719e-5*X*X*X - 0.0068901*X*X + 0.17305*X &
4884 elseif (Dr(i) .lt. 350.e-6) then
4885 Ef_rw = 9.0488e-5*X*X*X - 0.006585*X*X + 0.16606*X &
4888 Ef_rw = 0.00010721*X*X*X - 0.0072962*X*X + 0.1704*X &
4892 vtr = -0.1021 + 4.932E3*Dr(i) - 0.9551E6*Dr(i)*Dr(i) &
4893 + 0.07934E9*Dr(i)*Dr(i)*Dr(i) &
4894 - 0.002362E12*Dr(i)*Dr(i)*Dr(i)*Dr(i)
4895 stokes = Dc(j)*Dc(j)*vtr*rho_w/(9.*1.718E-5*Dr(i))
4896 reynolds = 9.*stokes/(p*p*rho_w)
4899 G = -0.1007D0 - 0.358D0*F + 0.0261D0*F*F
4901 z = DLOG(stokes/(K0+1.D-15))
4902 H = 0.1465D0 + 1.302D0*z - 0.607D0*z*z + 0.293D0*z*z*z
4903 yc0 = 2.0D0/PI * ATAN(H)
4904 Ef_rw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p))
4908 t_Efrw(i,j) = MAX(0.0, MIN(SNGL(Ef_rw), 0.95))
4913 end subroutine table_Efrw
4915 !+---+-----------------------------------------------------------------+
4916 !..Variable collision efficiency for snow collecting cloud water using
4917 !.. method of Wang and Ji, 2000 except equate melted snow diameter to
4918 !.. their "effective collision cross-section."
4919 !+---+-----------------------------------------------------------------+
4921 subroutine table_Efsw
4926 DOUBLE PRECISION:: Ds_m, vts, vtc, stokes, reynolds, Ef_sw
4927 DOUBLE PRECISION:: p, yc0, F, G, H, z, K0
4931 vtc = 1.19D4 * (1.0D4*Dc(j)*Dc(j)*0.25D0)
4933 vts = av_s*Ds(i)**bv_s * DEXP(-fv_s*Ds(i)) - vtc
4934 Ds_m = (am_s*Ds(i)**bm_s / am_r)**obmr
4936 if (p.gt.0.25 .or. Ds(i).lt.D0s .or. Dc(j).lt.6.E-6 &
4937 .or. vts.lt.1.E-3) then
4940 stokes = Dc(j)*Dc(j)*vts*rho_w/(9.*1.718E-5*Ds_m)
4941 reynolds = 9.*stokes/(p*p*rho_w)
4944 G = -0.1007D0 - 0.358D0*F + 0.0261D0*F*F
4946 z = DLOG(stokes/(K0+1.D-15))
4947 H = 0.1465D0 + 1.302D0*z - 0.607D0*z*z + 0.293D0*z*z*z
4948 yc0 = 2.0D0/PI * ATAN(H)
4949 Ef_sw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p))
4951 t_Efsw(i,j) = MAX(0.0, MIN(SNGL(Ef_sw), 0.95))
4957 end subroutine table_Efsw
4959 !+---+-----------------------------------------------------------------+
4960 !..Function to compute collision efficiency of collector species (rain,
4961 !.. snow, graupel) of aerosols. Follows Wang et al, 2010, ACP, which
4962 !.. follows Slinn (1983).
4963 !+---+-----------------------------------------------------------------+
4965 real function Eff_aero(D, Da, visc,rhoa,Temp,species)
4968 real:: D, Da, visc, rhoa, Temp
4969 character(LEN=1):: species
4970 real:: aval, Cc, diff, Re, Sc, St, St2, vt, Eff
4971 real, parameter:: boltzman = 1.3806503E-23
4972 real, parameter:: meanPath = 0.0256E-6
4975 if (species .eq. 'r') then
4976 vt = -0.1021 + 4.932E3*D - 0.9551E6*D*D &
4977 + 0.07934E9*D*D*D - 0.002362E12*D*D*D*D
4978 elseif (species .eq. 's') then
4980 elseif (species .eq. 'g') then
4981 vt = av_g(idx_bg1)*D**bv_g(idx_bg1)
4984 Cc = 1. + 2.*meanPath/Da *(1.257+0.4*exp(-0.55*Da/meanPath))
4985 diff = boltzman*Temp*Cc/(3.*PI*visc*Da)
4987 Re = 0.5*rhoa*D*vt/visc
4988 Sc = visc/(rhoa*diff)
4990 St = Da*Da*vt*1000./(9.*visc*D)
4991 aval = 1.+LOG(1.+Re)
4992 St2 = (1.2 + 1./12.*aval)/(1.+aval)
4994 Eff = 4./(Re*Sc) * (1. + 0.4*SQRT(Re)*Sc**0.3333 &
4995 + 0.16*SQRT(Re)*SQRT(Sc)) &
4996 + 4.*Da/D * (0.02 + Da/D*(1.+2.*SQRT(Re)))
4998 if (St.gt.St2) Eff = Eff + ( (St-St2)/(St-St2+0.666667))**1.5
4999 Eff_aero = MAX(1.E-5, MIN(Eff, 1.0))
5001 end function Eff_aero
5004 !+---+-----------------------------------------------------------------+
5005 !..Integrate rain size distribution from zero to D-star to compute the
5006 !.. number of drops smaller than D-star that evaporate in a single
5007 !.. timestep. Drops larger than D-star dont evaporate entirely so do
5008 !.. not affect number concentration.
5009 !+---+-----------------------------------------------------------------+
5011 subroutine table_dropEvap
5016 INTEGER:: i, j, k, n
5017 DOUBLE PRECISION, DIMENSION(nbc):: N_c, massc
5018 DOUBLE PRECISION:: summ, summ2, lamc, N0_c
5020 ! DOUBLE PRECISION:: Nt_r, N0, lam_exp, lam
5021 ! REAL:: xlimit_intg
5024 massc(n) = am_r*Dc(n)**bm_r
5028 nu_c = MIN(15, NINT(1000.E6/t_Nc(k)) + 2)
5030 lamc = (t_Nc(k)*am_r* ccg(2,nu_c)*ocg1(nu_c) / r_c(j))**obmr
5031 N0_c = t_Nc(k)*ocg1(nu_c) * lamc**cce(1,nu_c)
5033 !-GT tnc_wev(i,j,k) = GAMMP(nu_c+1., SNGL(Dc(i)*lamc))*t_Nc(k)
5034 N_c(i) = N0_c* Dc(i)**nu_c*EXP(-lamc*Dc(i))*dtc(i)
5035 ! if(j.eq.18 .and. k.eq.50) print*, ' N_c = ', N_c(i)
5039 summ = summ + massc(n)*N_c(n)
5040 summ2 = summ2 + N_c(n)
5042 ! if(j.eq.18 .and. k.eq.50) print*, ' DEBUG-TABLE: ', r_c(j), t_Nc(k), summ2, summ
5043 tpc_wev(i,j,k) = summ
5044 tnc_wev(i,j,k) = summ2
5050 !..To do the same thing for rain.
5054 ! lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(k))**ore1
5055 ! lam = lam_exp * (crg(3)*org2*org1)**obmr
5056 ! N0 = N0r_exp(j)/(crg(2)*lam_exp) * lam**cre(2)
5057 ! Nt_r = N0 * crg(2) / lam**cre(2)
5059 ! xlimit_intg = lam*Dr(i)
5060 ! tnr_rev(i,j,k) = GAMMP(mu_r+1.0, xlimit_intg) * Nt_r
5065 ! TO APPLY TABLE ABOVE
5066 !..Rain lookup table indexes.
5067 ! Dr_star = DSQRT(-2.D0*DT * t1_evap/(2.*PI) &
5068 ! * 0.78*4.*diffu(k)*xsat*rvs/rho_w)
5069 ! idx_d = NINT(1.0 + FLOAT(nbr) * DLOG(Dr_star/D0r) &
5070 ! / DLOG(Dr(nbr)/D0r))
5071 ! idx_d = MAX(1, MIN(idx_d, nbr))
5073 ! nir = NINT(ALOG10(rr(k)))
5074 ! do nn = nir-1, nir+1
5076 ! if ( (rr(k)/10.**nn).ge.1.0 .and. &
5077 ! (rr(k)/10.**nn).lt.10.0) goto 154
5080 ! idx_r = INT(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2)
5081 ! idx_r = MAX(1, MIN(idx_r, ntb_r))
5083 ! lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
5084 ! lam_exp = lamr * (crg(3)*org2*org1)**bm_r
5085 ! N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
5086 ! nir = NINT(DLOG10(N0_exp))
5087 ! do nn = nir-1, nir+1
5089 ! if ( (N0_exp/10.**nn).ge.1.0 .and. &
5090 ! (N0_exp/10.**nn).lt.10.0) goto 155
5093 ! idx_r1 = INT(N0_exp/10.**n) + 10*(n-nir3) - (n-nir3)
5094 ! idx_r1 = MAX(1, MIN(idx_r1, ntb_r1))
5096 ! pnr_rev(k) = MIN(nr(k)*odts, SNGL(tnr_rev(idx_d,idx_r1,idx_r) & ! RAIN2M
5099 end subroutine table_dropEvap
5102 !+---+-----------------------------------------------------------------+
5103 !..Fill the table of CCN activation data created from parcel model run
5104 !.. by Trude Eidhammer with inputs of aerosol number concentration,
5105 !.. vertical velocity, temperature, lognormal mean aerosol radius, and
5106 !.. hygroscopicity, kappa. The data are read from external file and
5107 !.. contain activated fraction of CCN for given conditions.
5108 !+---+-----------------------------------------------------------------+
5110 subroutine table_ccnAct
5116 LOGICAL, EXTERNAL:: wrf_dm_on_monitor
5119 INTEGER:: iunit_mp_th1, i
5121 CHARACTER*64 errmess
5124 IF ( wrf_dm_on_monitor() ) THEN
5126 INQUIRE ( i , OPENED = opened )
5127 IF ( .NOT. opened ) THEN
5134 #if defined(DM_PARALLEL) && !defined(STUBMPI)
5135 CALL wrf_dm_bcast_bytes ( iunit_mp_th1 , IWORDSIZE )
5137 IF ( iunit_mp_th1 < 0 ) THEN
5138 CALL wrf_error_fatal ( 'module_mp_thompson: table_ccnAct: '// &
5139 'Can not find unused fortran unit to read in lookup table.')
5142 IF ( wrf_dm_on_monitor() ) THEN
5143 WRITE(errmess, '(A,I2)') 'module_mp_thompson: opening CCN_ACTIVATE.BIN on unit ',iunit_mp_th1
5144 CALL wrf_debug(150, errmess)
5145 OPEN(iunit_mp_th1,FILE='CCN_ACTIVATE.BIN', &
5146 FORM='UNFORMATTED',STATUS='OLD',ERR=9009)
5149 #define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes(A, size(A)*R4SIZE)
5151 IF ( wrf_dm_on_monitor() ) READ(iunit_mp_th1,ERR=9010) tnccn_act
5152 #if defined(DM_PARALLEL) && !defined(STUBMPI)
5153 DM_BCAST_MACRO(tnccn_act)
5159 WRITE( errmess , '(A,I2)' ) 'module_mp_thompson: error opening CCN_ACTIVATE.BIN on unit ',iunit_mp_th1
5160 CALL wrf_error_fatal(errmess)
5163 WRITE( errmess , '(A,I2)' ) 'module_mp_thompson: error reading CCN_ACTIVATE.BIN on unit ',iunit_mp_th1
5164 CALL wrf_error_fatal(errmess)
5166 end subroutine table_ccnAct
5168 !+---+-----------------------------------------------------------------+
5169 !..Retrieve fraction of CCN that gets activated given the model temp,
5170 !.. vertical velocity, and available CCN concentration. The lookup
5171 !.. table (read from external file) has CCN concentration varying the
5172 !.. quickest, then updraft, then temperature, then mean aerosol radius,
5173 !.. and finally hygroscopicity, kappa.
5174 !.. TO_DO ITEM: For radiation cooling producing fog, in which case the
5175 !.. updraft velocity could easily be negative, we could use the temp
5176 !.. and its tendency to diagnose a pretend postive updraft velocity.
5177 !+---+-----------------------------------------------------------------+
5178 real function activ_ncloud(Tt, Ww, NCCN)
5181 REAL, INTENT(IN):: Tt, Ww, NCCN
5182 REAL:: n_local, w_local
5183 INTEGER:: i, j, k, l, m, n
5184 REAL:: A, B, C, D, t, u, x1, x2, y1, y2, nx, wy, fraction
5187 ! ta_Na = (/10.0, 31.6, 100.0, 316.0, 1000.0, 3160.0, 10000.0/) ntb_arc
5188 ! ta_Ww = (/0.01, 0.0316, 0.1, 0.316, 1.0, 3.16, 10.0, 31.6, 100.0/) ntb_arw
5189 ! ta_Tk = (/243.15, 253.15, 263.15, 273.15, 283.15, 293.15, 303.15/) ntb_art
5190 ! ta_Ra = (/0.01, 0.02, 0.04, 0.08, 0.16/) ntb_arr
5191 ! ta_Ka = (/0.2, 0.4, 0.6, 0.8/) ntb_ark
5193 n_local = NCCN * 1.E-6
5196 if (n_local .ge. ta_Na(ntb_arc)) then
5197 n_local = ta_Na(ntb_arc) - 1.0
5198 elseif (n_local .le. ta_Na(1)) then
5199 n_local = ta_Na(1) + 1.0
5202 if (n_local.ge.ta_Na(n-1) .and. n_local.lt.ta_Na(n)) goto 8003
5206 x1 = LOG(ta_Na(i-1))
5209 if (w_local .ge. ta_Ww(ntb_arw)) then
5210 w_local = ta_Ww(ntb_arw) - 1.0
5211 elseif (w_local .le. ta_Ww(1)) then
5212 w_local = ta_Ww(1) + 0.001
5215 if (w_local.ge.ta_Ww(n-1) .and. w_local.lt.ta_Ww(n)) goto 8005
5219 y1 = LOG(ta_Ww(j-1))
5222 k = MAX(1, MIN( NINT( (Tt - ta_Tk(1))*0.1) + 1, ntb_art))
5224 !..The next two values are indexes of mean aerosol radius and
5225 !.. hygroscopicity. Currently these are constant but a future version
5226 !.. should implement other variables to allow more freedom such as
5227 !.. at least simple separation of tiny size sulfates from larger
5232 A = tnccn_act(i-1,j-1,k,l,m)
5233 B = tnccn_act(i,j-1,k,l,m)
5234 C = tnccn_act(i,j,k,l,m)
5235 D = tnccn_act(i-1,j,k,l,m)
5242 ! t = (n_local-ta(Na(i-1))/(ta_Na(i)-ta_Na(i-1))
5243 ! u = (w_local-ta_Ww(j-1))/(ta_Ww(j)-ta_Ww(j-1))
5245 fraction = (1.0-t)*(1.0-u)*A + t*(1.0-u)*B + t*u*C + (1.0-t)*u*D
5247 ! if (NCCN*fraction .gt. 0.75*Nt_c_max) then
5248 ! write(*,*) ' DEBUG-GT ', n_local, w_local, Tt, i, j, k
5251 activ_ncloud = NCCN*fraction
5253 end function activ_ncloud
5255 !+---+-----------------------------------------------------------------+
5256 !+---+-----------------------------------------------------------------+
5257 SUBROUTINE GCF(GAMMCF,A,X,GLN)
5258 ! --- RETURNS THE INCOMPLETE GAMMA FUNCTION Q(A,X) EVALUATED BY ITS
5259 ! --- CONTINUED FRACTION REPRESENTATION AS GAMMCF. ALSO RETURNS
5260 ! --- LN(GAMMA(A)) AS GLN. THE CONTINUED FRACTION IS EVALUATED BY
5261 ! --- A MODIFIED LENTZ METHOD.
5264 INTEGER, PARAMETER:: ITMAX=100
5265 REAL, PARAMETER:: gEPS=3.E-7
5266 REAL, PARAMETER:: FPMIN=1.E-30
5267 REAL, INTENT(IN):: A, X
5270 REAL:: AN,B,C,D,DEL,H
5280 IF(ABS(D).LT.FPMIN)D=FPMIN
5282 IF(ABS(C).LT.FPMIN)C=FPMIN
5286 IF(ABS(DEL-1.).LT.gEPS)GOTO 1
5288 PRINT *, 'A TOO LARGE, ITMAX TOO SMALL IN GCF'
5289 1 GAMMCF=EXP(-X+A*LOG(X)-GLN)*H
5291 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
5292 !+---+-----------------------------------------------------------------+
5293 SUBROUTINE GSER(GAMSER,A,X,GLN)
5294 ! --- RETURNS THE INCOMPLETE GAMMA FUNCTION P(A,X) EVALUATED BY ITS
5295 ! --- ITS SERIES REPRESENTATION AS GAMSER. ALSO RETURNS LN(GAMMA(A))
5299 INTEGER, PARAMETER:: ITMAX=100
5300 REAL, PARAMETER:: gEPS=3.E-7
5301 REAL, INTENT(IN):: A, X
5307 IF(X.LT.0.) PRINT *, 'X < 0 IN GSER'
5318 IF(ABS(DEL).LT.ABS(SUM)*gEPS)GOTO 1
5320 PRINT *,'A TOO LARGE, ITMAX TOO SMALL IN GSER'
5321 1 GAMSER=SUM*EXP(-X+A*LOG(X)-GLN)
5323 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
5324 !+---+-----------------------------------------------------------------+
5325 REAL FUNCTION GAMMLN(XX)
5326 ! --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 0.
5328 REAL, INTENT(IN):: XX
5329 DOUBLE PRECISION, PARAMETER:: STP = 2.5066282746310005D0
5330 DOUBLE PRECISION, DIMENSION(6), PARAMETER:: &
5331 COF = (/76.18009172947146D0, -86.50532032941677D0, &
5332 24.01409824083091D0, -1.231739572450155D0, &
5333 .1208650973866179D-2, -.5395239384953D-5/)
5334 DOUBLE PRECISION:: SER,TMP,X,Y
5340 TMP=(X+0.5D0)*LOG(TMP)-TMP
5341 SER=1.000000000190015D0
5346 GAMMLN=TMP+LOG(STP*SER/X)
5348 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
5349 !+---+-----------------------------------------------------------------+
5350 REAL FUNCTION GAMMP(A,X)
5351 ! --- COMPUTES THE INCOMPLETE GAMMA FUNCTION P(A,X)
5352 ! --- SEE ABRAMOWITZ AND STEGUN 6.5.1
5355 REAL, INTENT(IN):: A,X
5356 REAL:: GAMMCF,GAMSER,GLN
5358 IF((X.LT.0.) .OR. (A.LE.0.)) THEN
5359 PRINT *, 'BAD ARGUMENTS IN GAMMP'
5361 ELSEIF(X.LT.A+1.)THEN
5362 CALL GSER(GAMSER,A,X,GLN)
5365 CALL GCF(GAMMCF,A,X,GLN)
5369 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
5370 !+---+-----------------------------------------------------------------+
5371 REAL FUNCTION WGAMMA(y)
5374 REAL, INTENT(IN):: y
5376 WGAMMA = EXP(GAMMLN(y))
5379 !+---+-----------------------------------------------------------------+
5380 ! THIS FUNCTION CALCULATES THE LIQUID SATURATION VAPOR MIXING RATIO AS
5381 ! A FUNCTION OF TEMPERATURE AND PRESSURE
5383 REAL FUNCTION RSLF(P,T)
5386 REAL, INTENT(IN):: P, T
5388 REAL, PARAMETER:: C0= .611583699E03
5389 REAL, PARAMETER:: C1= .444606896E02
5390 REAL, PARAMETER:: C2= .143177157E01
5391 REAL, PARAMETER:: C3= .264224321E-1
5392 REAL, PARAMETER:: C4= .299291081E-3
5393 REAL, PARAMETER:: C5= .203154182E-5
5394 REAL, PARAMETER:: C6= .702620698E-8
5395 REAL, PARAMETER:: C7= .379534310E-11
5396 REAL, PARAMETER:: C8=-.321582393E-13
5398 X=MAX(-80.,T-273.16)
5400 ! ESL=612.2*EXP(17.67*X/(T-29.65))
5401 ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
5402 ESL=MIN(ESL, P*0.15) ! Even with P=1050mb and T=55C, the sat. vap. pres only contributes to ~15% of total pres.
5403 RSLF=.622*ESL/(P-ESL)
5406 ! ; Source: Murphy and Koop, Review of the vapour pressure of ice and
5407 ! supercooled water for atmospheric applications, Q. J. R.
5408 ! Meteorol. Soc (2005), 131, pp. 1539-1565.
5409 ! ESL = EXP(54.842763 - 6763.22 / T - 4.210 * ALOG(T) + 0.000367 * T
5410 ! + TANH(0.0415 * (T - 218.8)) * (53.878 - 1331.22
5411 ! / T - 9.44523 * ALOG(T) + 0.014025 * T))
5414 !+---+-----------------------------------------------------------------+
5415 ! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A
5416 ! FUNCTION OF TEMPERATURE AND PRESSURE
5418 REAL FUNCTION RSIF(P,T)
5421 REAL, INTENT(IN):: P, T
5423 REAL, PARAMETER:: C0= .609868993E03
5424 REAL, PARAMETER:: C1= .499320233E02
5425 REAL, PARAMETER:: C2= .184672631E01
5426 REAL, PARAMETER:: C3= .402737184E-1
5427 REAL, PARAMETER:: C4= .565392987E-3
5428 REAL, PARAMETER:: C5= .521693933E-5
5429 REAL, PARAMETER:: C6= .307839583E-7
5430 REAL, PARAMETER:: C7= .105785160E-9
5431 REAL, PARAMETER:: C8= .161444444E-12
5433 X=MAX(-80.,T-273.16)
5434 ESI=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
5435 ESI=MIN(ESI, P*0.15)
5436 RSIF=.622*ESI/max(1.e-4,(P-ESI))
5439 ! ; Source: Murphy and Koop, Review of the vapour pressure of ice and
5440 ! supercooled water for atmospheric applications, Q. J. R.
5441 ! Meteorol. Soc (2005), 131, pp. 1539-1565.
5442 ! ESI = EXP(9.550426 - 5723.265/T + 3.53068*ALOG(T) - 0.00728332*T)
5446 !+---+-----------------------------------------------------------------+
5447 real function iceDeMott(tempc, qv, qvs, qvsi, rho, nifa)
5450 REAL, INTENT(IN):: tempc, qv, qvs, qvsi, rho, nifa
5453 REAL:: satw, sati, siw, p_x, si0x, dtt, dsi, dsw, dab, fc, hx
5454 REAL:: ntilde, n_in, nmax, nhat, mux, xni, nifa_cc
5455 REAL, PARAMETER:: p_c1 = 1000.
5456 REAL, PARAMETER:: p_rho_c = 0.76
5457 REAL, PARAMETER:: p_alpha = 1.0
5458 REAL, PARAMETER:: p_gam = 2.
5459 REAL, PARAMETER:: delT = 5.
5460 REAL, PARAMETER:: T0x = -40.
5461 REAL, PARAMETER:: Sw0x = 0.97
5462 REAL, PARAMETER:: delSi = 0.1
5463 REAL, PARAMETER:: hdm = 0.15
5464 REAL, PARAMETER:: p_psi = 0.058707*p_gam/p_rho_c
5465 REAL, PARAMETER:: aap = 1.
5466 REAL, PARAMETER:: bbp = 0.
5467 REAL, PARAMETER:: y1p = -35.
5468 REAL, PARAMETER:: y2p = -25.
5469 REAL, PARAMETER:: rho_not0 = 101325./(287.05*273.15)
5477 ! p_x = -1.0261+(3.1656e-3*tempc)+(5.3938e-4*(tempc*tempc)) &
5478 ! + (8.2584e-6*(tempc*tempc*tempc))
5479 ! si0x = 1.+(10.**p_x)
5480 ! if (sati.ge.si0x .and. satw.lt.0.985) then
5481 ! dtt = delta_p (tempc, T0x, T0x+delT, 1., hdm)
5482 ! dsi = delta_p (sati, Si0x, Si0x+delSi, 0., 1.)
5483 ! dsw = delta_p (satw, Sw0x, 1., 0., 1.)
5485 ! hx = min(fc+((1.-fc)*dsw), 1.)
5486 ! ntilde = p_c1*p_gam*((exp(12.96*(sati-1.1)))**0.3) / p_rho_c
5487 ! if (tempc .le. y1p) then
5489 ! elseif (tempc .ge. y2p) then
5490 ! n_in = p_psi*p_c1*exp(12.96*(sati-1.)-0.639)
5492 ! if (tempc .le. -30.) then
5493 ! nmax = p_c1*p_gam*(exp(12.96*(siw-1.1)))**0.3/p_rho_c
5495 ! nmax = p_psi*p_c1*exp(12.96*(siw-1.)-0.639)
5497 ! ntilde = MIN(ntilde, nmax)
5498 ! nhat = MIN(p_psi*p_c1*exp(12.96*(sati-1.)-0.639), nmax)
5499 ! dab = delta_p (tempc, y1p, y2p, aap, bbp)
5500 ! n_in = MIN(nhat*(ntilde/nhat)**dab, nmax)
5502 ! mux = hx*p_alpha*n_in*rho
5503 ! xni = mux*((6700.*nifa)-200.)/((6700.*5.E5)-200.)
5504 ! elseif (satw.ge.0.985 .and. tempc.gt.HGFR-273.15) then
5505 nifa_cc = MAX(0.5, nifa*RHO_NOT0*1.E-6/rho)
5506 ! xni = 3.*nifa_cc**(1.25)*exp((0.46*(-tempc))-11.6) ! [DeMott, 2015]
5507 xni = (5.94e-5*(-tempc)**3.33) & ! [DeMott, 2010]
5508 * (nifa_cc**((-0.0264*(tempc))+0.0033))
5509 xni = xni*rho/RHO_NOT0 * 1000.
5512 iceDeMott = MAX(0., xni)
5514 end FUNCTION iceDeMott
5516 !+---+-----------------------------------------------------------------+
5517 !..Newer research since Koop et al (2001) suggests that the freezing
5518 !.. rate should be lower than original paper, so J_rate is reduced
5519 !.. by two orders of magnitude.
5521 real function iceKoop(temp, qv, qvs, naero, dt)
5524 REAL, INTENT(IN):: temp, qv, qvs, naero, DT
5525 REAL:: mu_diff, a_w_i, delta_aw, log_J_rate, J_rate, prob_h, satw
5530 mu_diff = 210368.0 + (131.438*temp) - (3.32373E6/temp) &
5531 & - (41729.1*alog(temp))
5532 a_w_i = exp(mu_diff/(R_uni*temp))
5533 delta_aw = satw - a_w_i
5534 log_J_rate = -906.7 + (8502.0*delta_aw) &
5535 & - (26924.0*delta_aw*delta_aw) &
5536 & + (29180.0*delta_aw*delta_aw*delta_aw)
5537 log_J_rate = MIN(20.0, log_J_rate)
5538 J_rate = 10.**log_J_rate ! cm-3 s-1
5539 prob_h = MIN(1.-exp(-J_rate*ar_volume*DT), 1.)
5540 if (prob_h .gt. 0.) then
5541 xni = MIN(prob_h*naero, 1000.E3)
5544 iceKoop = MAX(0.0, xni)
5546 end FUNCTION iceKoop
5548 !+---+-----------------------------------------------------------------+
5549 !.. Helper routine for Phillips et al (2008) ice nucleation. Trude
5551 REAL FUNCTION delta_p (yy, y1, y2, aa, bb)
5554 REAL, INTENT(IN):: yy, y1, y2, aa, bb
5555 REAL:: dab, A, B, a0, a1, a2, a3
5557 A = 6.*(aa-bb)/((y2-y1)*(y2-y1)*(y2-y1))
5558 B = aa+(A*y1*y1*y1/6.)-(A*y1*y1*y2*0.5)
5566 else if (yy.ge.y2) then
5569 dab = a0+(a1*yy)+(a2*yy*yy)+(a3*yy*yy*yy)
5580 END FUNCTION delta_p
5582 !+---+-----------------------------------------------------------------+
5585 !+---+-----------------------------------------------------------------+
5586 !..Compute _radiation_ effective radii of cloud water, ice, and snow.
5587 !.. These are entirely consistent with microphysics assumptions, not
5588 !.. constant or otherwise ad hoc as is internal to most radiation
5589 !.. schemes. Since only the smallest snowflakes should impact
5590 !.. radiation, compute from first portion of complicated Field number
5591 !.. distribution, not the second part, which is the larger sizes.
5592 !+---+-----------------------------------------------------------------+
5594 subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
5595 & re_qc1d, re_qi1d, re_qs1d, kts, kte)
5600 INTEGER, INTENT(IN):: kts, kte
5601 REAL, DIMENSION(kts:kte), INTENT(IN):: &
5602 & t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d
5603 REAL, DIMENSION(kts:kte), INTENT(INOUT):: re_qc1d, re_qi1d, re_qs1d
5606 REAL, DIMENSION(kts:kte):: rho, rc, nc, ri, ni, rs
5607 REAL:: smo2, smob, smoc
5608 REAL:: tc0, loga_, a_, b_
5609 DOUBLE PRECISION:: lamc, lami
5610 LOGICAL:: has_qc, has_qi, has_qs
5612 real, dimension(15), parameter:: g_ratio = (/24,60,120,210,336, &
5613 & 504,720,990,1320,1716,2184,2730,3360,4080,4896/)
5619 re_qc1d(:) = RE_QC_BG
5620 re_qi1d(:) = RE_QI_BG
5621 re_qs1d(:) = RE_QS_BG
5624 rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622))
5625 rc(k) = MAX(R1, qc1d(k)*rho(k))
5626 nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max))
5627 if (.NOT. is_aerosol_aware) nc(k) = Nt_c
5628 if (rc(k).gt.R1 .and. nc(k).gt.R2) has_qc = .true.
5629 ri(k) = MAX(R1, qi1d(k)*rho(k))
5630 ni(k) = MAX(R2, ni1d(k)*rho(k))
5631 if (ri(k).gt.R1 .and. ni(k).gt.R2) has_qi = .true.
5632 rs(k) = MAX(R1, qs1d(k)*rho(k))
5633 if (rs(k).gt.R1) has_qs = .true.
5638 if (rc(k).le.R1 .or. nc(k).le.R2) CYCLE
5639 if (nc(k).lt.100) then
5641 elseif (nc(k).gt.1.E10) then
5644 inu_c = MIN(15, NINT(1000.E6/nc(k)) + 2)
5646 lamc = (nc(k)*am_r*g_ratio(inu_c)/rc(k))**obmr
5647 re_qc1d(k) = MAX(2.51E-6, MIN(SNGL(0.5D0 * DBLE(3.+inu_c)/lamc), 50.E-6))
5653 if (ri(k).le.R1 .or. ni(k).le.R2) CYCLE
5654 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
5655 re_qi1d(k) = MAX(2.51E-6, MIN(SNGL(0.5D0 * DBLE(3.+mu_i)/lami), 125.E-6))
5661 if (rs(k).le.R1) CYCLE
5662 tc0 = MIN(-0.1, t1d(k)-273.15)
5665 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
5666 !.. then we must compute actual 2nd moment and use as reference.
5667 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
5670 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
5671 & + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
5672 & + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
5673 & + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
5674 & + sa(10)*bm_s*bm_s*bm_s
5676 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
5677 & + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
5678 & + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
5679 & + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
5680 & + sb(10)*bm_s*bm_s*bm_s
5681 smo2 = (smob/a_)**(1./b_)
5683 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
5684 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
5685 & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
5686 & + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
5687 & + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
5688 & + sa(10)*cse(1)*cse(1)*cse(1)
5690 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
5691 & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
5692 & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
5693 & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
5694 smoc = a_ * smo2**b_
5695 re_qs1d(k) = MAX(5.01E-6, MIN(0.5*(smoc/smob), 999.E-6))
5699 end subroutine calc_effectRad
5701 !+---+-----------------------------------------------------------------+
5702 !..Compute radar reflectivity assuming 10 cm wavelength radar and using
5703 !.. Rayleigh approximation. Only complication is melted snow/graupel
5704 !.. which we treat as water-coated ice spheres and use Uli Blahak's
5705 !.. library of routines. The meltwater fraction is simply the amount
5706 !.. of frozen species remaining from what initially existed at the
5707 !.. melting level interface.
5708 !+---+-----------------------------------------------------------------+
5710 subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, &
5711 ng1d, qb1d, t1d, p1d, dBZ, kts, kte, ii, jj, ke_diag)
5716 INTEGER, INTENT(IN):: kts, kte, ii, jj, ke_diag
5717 REAL, DIMENSION(kts:kte), INTENT(IN):: &
5718 qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, ng1d, qb1d, t1d, p1d
5719 REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
5720 ! REAL, DIMENSION(kts:kte), INTENT(INOUT):: vt_dBZ
5723 REAL, DIMENSION(kts:kte):: temp, pres, qv, rho, rhof
5724 REAL, DIMENSION(kts:kte):: rc, rr, nr, rs, rg, ng, rb
5726 DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g
5727 REAL, DIMENSION(kts:kte):: mvd_r
5728 REAL, DIMENSION(kts:kte):: smob, smo2, smoc, smoz
5729 REAL:: oM3, M0, Mrat, slam1, slam2, xDs
5730 REAL:: ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts
5731 REAL:: vtr_dbz_wt, vts_dbz_wt, vtg_dbz_wt
5733 REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
5734 INTEGER, DIMENSION(kts:kte):: idx_bg
5736 DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamr, lamg
5737 REAL:: a_, b_, loga_, tc0
5738 DOUBLE PRECISION:: fmelt_s, fmelt_g
5740 INTEGER:: i, k, k_0, kbot, n, ktop
5742 LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
5744 DOUBLE PRECISION:: cback, x, eta, f_d
5745 REAL:: xslw1, ygra1, zans1
5754 !+---+-----------------------------------------------------------------+
5755 !..Put column of data into local arrays.
5756 !+---+-----------------------------------------------------------------+
5760 qv(k) = MAX(1.E-10, qv1d(k))
5762 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
5763 rhof(k) = SQRT(RHO_NOT/rho(k))
5764 rc(k) = MAX(R1, qc1d(k)*rho(k))
5765 if (qr1d(k) .gt. R1) then
5766 rr(k) = qr1d(k)*rho(k)
5767 nr(k) = MAX(R2, nr1d(k)*rho(k))
5768 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
5770 N0_r(k) = nr(k)*org2*lamr**cre(2)
5771 mvd_r(k) = (3.0 + mu_r + 0.672) * ilamr(k)
5779 if (qs1d(k) .gt. R2) then
5780 rs(k) = qs1d(k)*rho(k)
5786 if (qg1d(k) .gt. R2) then
5787 rg(k) = qg1d(k)*rho(k)
5788 ng(k) = MAX(R2, ng1d(k)*rho(k))
5789 rb(k) = MAX(qg1d(k)/rho_g(NRHG), qb1d(k))
5790 rb(k) = MIN(qg1d(k)/rho_g(1), rb(k))
5791 idx_bg(k) = MAX(1,MIN(NINT(qg1d(k)/rb(k) *0.01)+1,NRHG))
5792 if (.not. is_hail_aware) idx_bg(k) = idx_bg1
5802 !+---+-----------------------------------------------------------------+
5803 !..Calculate y-intercept, slope, and useful moments for snow.
5804 !+---+-----------------------------------------------------------------+
5811 if ( ( ke_diag > kts .and. ANY(L_qs .eqv. .true.) ) .or. &
5812 (ke_diag == kts .and. L_qs(kts) .eqv. .true. ) ) then
5813 do k = kts, ke_diag ! kte
5814 if (.not. L_qs(k)) CYCLE
5815 tc0 = MIN(-0.1, temp(k)-273.15)
5816 smob(k) = rs(k)*oams
5818 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
5819 !.. then we must compute actual 2nd moment and use as reference.
5820 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
5823 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
5824 & + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
5825 & + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
5826 & + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
5827 & + sa(10)*bm_s*bm_s*bm_s
5829 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
5830 & + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
5831 & + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
5832 & + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
5833 & + sb(10)*bm_s*bm_s*bm_s
5834 smo2(k) = (smob(k)/a_)**(1./b_)
5837 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
5838 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
5839 & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
5840 & + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
5841 & + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
5842 & + sa(10)*cse(1)*cse(1)*cse(1)
5844 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
5845 & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
5846 & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
5847 & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
5848 smoc(k) = a_ * smo2(k)**b_
5850 !..Calculate bm_s*2 (th) moment. Useful for reflectivity.
5851 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(3) &
5852 & + sa(4)*tc0*cse(3) + sa(5)*tc0*tc0 &
5853 & + sa(6)*cse(3)*cse(3) + sa(7)*tc0*tc0*cse(3) &
5854 & + sa(8)*tc0*cse(3)*cse(3) + sa(9)*tc0*tc0*tc0 &
5855 & + sa(10)*cse(3)*cse(3)*cse(3)
5857 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(3) + sb(4)*tc0*cse(3) &
5858 & + sb(5)*tc0*tc0 + sb(6)*cse(3)*cse(3) &
5859 & + sb(7)*tc0*tc0*cse(3) + sb(8)*tc0*cse(3)*cse(3) &
5860 & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(3)*cse(3)*cse(3)
5861 smoz(k) = a_ * smo2(k)**b_
5865 !+---+-----------------------------------------------------------------+
5866 !..Calculate y-intercept, slope values for graupel.
5867 !+---+-----------------------------------------------------------------+
5869 if (ANY(L_qg .eqv. .true.)) then
5871 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg
5873 N0_g(k) = ng(k)*ogg2*lamg**cge(2,1)
5880 !+---+-----------------------------------------------------------------+
5881 !..Locate K-level of start of melting (k_0 is level above).
5882 !+---+-----------------------------------------------------------------+
5885 do k = kte-1, kts, -1
5886 if ( (temp(k).gt.273.15) .and. L_qr(k) &
5887 & .and. (L_qs(k+1).or.L_qg(k+1)) ) then
5895 ! Set loop limit for wet ice according to whether the full 3D field is needed or just k=1
5896 k_0loop = Min(k_0, ke_diag+1)
5898 !+---+-----------------------------------------------------------------+
5899 !..Do not do the so-called bright band if using the variable density
5900 !.. graupel-hail category since the density increases during melting and
5901 !.. will account for bright band behavior explicitly.
5902 !+---+-----------------------------------------------------------------+
5903 ! if (is_hail_aware) melti = .false.
5905 !+---+-----------------------------------------------------------------+
5906 !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
5907 !.. and non-water-coated snow and graupel when below freezing are
5908 !.. simple. Integrations of m(D)*m(D)*N(D)*dD.
5909 !+---+-----------------------------------------------------------------+
5911 do k = kts, ke_diag ! kte
5914 ze_graupel(k) = 1.e-22
5915 if (L_qr(k)) ze_rain(k) = N0_r(k)*crg(4)*ilamr(k)**cre(4)
5916 if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) &
5917 & * (am_s/900.0)*(am_s/900.0)*smoz(k)
5918 if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) &
5919 & * (am_g(idx_bg(k))/900.0)*(am_g(idx_bg(k))/900.0) &
5920 & * N0_g(k)*cgg(4,1)*ilamg(k)**cge(4,1)
5923 !+---+-----------------------------------------------------------------+
5924 !..Special case of melting ice (snow/graupel) particles. Assume the
5925 !.. ice is surrounded by the liquid water. Fraction of meltwater is
5926 !.. extremely simple based on amount found above the melting level.
5927 !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
5929 !+---+-----------------------------------------------------------------+
5931 if (.not. iiwarm .and. melti .and. k_0.ge.2) then
5932 do k = k_0loop-1, kts, -1
5934 !..Reflectivity contributed by melting snow
5935 if (L_qs(k) .and. L_qs(k_0) ) then
5936 fmelt_s = MAX(0.05d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
5940 Mrat = smob(k)*M0*M0*M0
5944 x = am_s * xxDs(n)**bm_s
5945 call rayleigh_soak_wetgraupel (x, DBLE(ocms), DBLE(obms), &
5946 & fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
5947 & CBACK, mixingrulestring_s, matrixstring_s, &
5948 & inclusionstring_s, hoststring_s, &
5949 & hostmatrixstring_s, hostinclusionstring_s)
5950 f_d = Mrat*(Kap0*DEXP(-slam1*xxDs(n)) &
5951 & + Kap1*(M0*xxDs(n))**mu_s * DEXP(-slam2*xxDs(n)))
5952 eta = eta + f_d * CBACK * simpson(n) * xdts(n)
5954 ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
5957 !..Reflectivity contributed by melting graupel
5959 ! if (L_qg(k) .and. L_qg(k_0) ) then
5960 ! fmelt_g = MAX(0.05d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0))
5962 ! lamg = 1./ilamg(k)
5964 ! x = am_g(idx_bg(k)) * xxDg(n)**bm_g
5965 ! call rayleigh_soak_wetgraupel (x, DBLE(ocmg(idx_bg(k))), DBLE(obmg), &
5966 ! & fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
5967 ! & CBACK, mixingrulestring_g, matrixstring_g, &
5968 ! & inclusionstring_g, hoststring_g, &
5969 ! & hostmatrixstring_g, hostinclusionstring_g)
5970 ! f_d = N0_g(k)*xxDg(n)**mu_g * DEXP(-lamg*xxDg(n))
5971 ! eta = eta + f_d * CBACK * simpson(n) * xdtg(n)
5973 ! ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
5979 do k = ke_diag, kts, -1
5980 dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
5984 !..Reflectivity-weighted terminal velocity (snow, rain, graupel, mix).
5985 ! do k = kte, kts, -1
5987 ! if (rs(k).gt.R2) then
5988 ! Mrat = smob(k) / smoc(k)
5989 ! ils1 = 1./(Mrat*Lam0 + fv_s)
5990 ! ils2 = 1./(Mrat*Lam1 + fv_s)
5991 ! t1_vts = Kap0*csg(5)*ils1**cse(5)
5992 ! t2_vts = Kap1*Mrat**mu_s*csg(11)*ils2**cse(11)
5993 ! ils1 = 1./(Mrat*Lam0)
5994 ! ils2 = 1./(Mrat*Lam1)
5995 ! t3_vts = Kap0*csg(6)*ils1**cse(6)
5996 ! t4_vts = Kap1*Mrat**mu_s*csg(12)*ils2**cse(12)
5997 ! vts_dbz_wt = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts)
5998 ! if (temp(k).ge.273.15 .and. temp(k).lt.275.15) then
5999 ! vts_dbz_wt = vts_dbz_wt*1.5
6000 ! elseif (temp(k).ge.275.15) then
6001 ! vts_dbz_wt = vts_dbz_wt*2.0
6004 ! vts_dbz_wt = 1.E-3
6007 ! if (rr(k).gt.R1) then
6008 ! lamr = 1./ilamr(k)
6009 ! vtr_dbz_wt = rhof(k)*av_r*crg(13)*(lamr+fv_r)**(-cre(13)) &
6010 ! & / (crg(4)*lamr**(-cre(4)))
6012 ! vtr_dbz_wt = 1.E-3
6015 ! if (rg(k).gt.R2) then
6016 ! lamg = 1./ilamg(k)
6017 ! vtg_dbz_wt = rhof(k)*av_g(idx_bg(k))*cgg(5,idx_bg(k))*lamg**(-cge(5,idx_bg(k))) &
6018 ! & / (cgg(4,1)*lamg**(-cge(4,1)))
6020 ! vtg_dbz_wt = 1.E-3
6023 ! vt_dBZ(k) = (vts_dbz_wt*ze_snow(k) + vtr_dbz_wt*ze_rain(k) &
6024 ! & + vtg_dbz_wt*ze_graupel(k)) &
6025 ! & / (ze_rain(k)+ze_snow(k)+ze_graupel(k))
6028 end subroutine calc_refl10cm
6030 !+---+-----------------------------------------------------------------+
6032 real function theta_e(pres_Pa,temp_K,w_non,tlcl_K)
6034 !.. The following code was based on Bolton (1980) eqn #43
6035 !.. and claims to have 0.3 K maximum error within -35 < T < 35 C
6036 !.. pres_Pa = Pressure in Pascals
6037 !.. temp_K = Temperature in Kelvin
6038 !.. w_non = mixing ratio (non-dimensional = kg/kg)
6039 !.. tlcl_K = Temperature at Lifting Condensation Level (K)
6043 real:: pres_Pa, temp_K, w_non, tlcl_K
6044 real:: pp, tt, rr, tlc, power, xx, p1, p2
6053 power=(0.2854*(1.0 - (0.28*rr) ))
6054 xx = tt * (100000.0/pp)**power
6056 p1 = (3.376/tlc) - 0.00254
6057 p2 = (rr*1000.0) * (1.0 + 0.81*rr)
6059 theta_e = xx * exp(p1*p2)
6064 !+---+-----------------------------------------------------------------+
6066 real function t_lcl(temp_K,tdew_K)
6068 !.. The following code was based on Bolton (1980) eqn #15
6069 !.. and claims to have 0.1 K maximum error within -35 < T < 35 C
6070 !.. temp_K = Temperature in Kelvin
6071 !.. tdew_K = Dewpoint T at Lifting Condensation Level (K)
6075 real:: temp_K, tdew_K
6076 real:: tt, tttd, denom
6082 denom= ( 1.0/(tttd-56.0) ) + (log(tt/tttd)/800.)
6083 t_lcl = ( 1.0 / denom ) + 56.0
6087 !+---+-----------------------------------------------------------------+
6089 real function t_dew(pres_Pa,w_non)
6091 !.. pres_Pa = Pressure in Pascals
6092 !.. w_non = mixing ratio (non-dimensional = kg/kg)
6096 real:: pres_Pa, w_non
6097 real:: p, RR, ES, ESLN
6105 T_Dew=(35.86*ESLN-4947.2325)/(ESLN-23.6837)
6109 !+---+-----------------------------------------------------------------+
6111 real function theta_wetb(thetae_K)
6113 !.. Eqn below was gotten from polynomial fit to data in
6114 !.. Smithsonian Meteorological Tables showing Theta-e
6124 real*8 c(0:6), d(0:6)
6125 data c/-1.00922292e-10, -1.47945344e-8, -1.7303757e-6 &
6126 & ,-0.00012709, 1.15849867e-6, -3.518296861e-9 &
6128 data d/0.00000000, -3.5223513e-10, -5.7250807e-8 &
6129 & ,-5.83975422e-6, 4.72445163e-8, -1.13402845e-10 &
6132 x=min(475.0,thetae_K)
6134 if( x .le. 335.5 ) then
6135 answer = c(0)+x*(c(1)+x*(c(2)+x*(c(3)+x*(c(4)+x*(c(5)+ &
6138 answer = d(0)+x*(d(1)+x*(d(2)+x*(d(3)+x*(d(4)+x*(d(5)+ &
6142 theta_wetb = answer + 273.15
6147 !+---+-----------------------------------------------------------------+
6149 real function compT_fr_The(thelcl_K,pres_Pa)
6151 !.. pres_Pa = Pressure in Pascals
6152 !.. thelcl = Theta-e at LCL (units in Kelvin)
6154 !.. Temperature (K) is returned given Theta-e at LCL
6155 !.. and a pressure. This describes a moist-adiabat.
6156 !.. This temperature is the parcel temp at level Pres
6157 !.. along moist adiabat described by theta-e.
6161 real:: thelcl_K, pres_Pa
6162 real:: guess, epsilon, w1, w2, tenu, tenup, cor, thwlcl_K
6167 guess= (thelcl_K - 0.5 * ( max(thelcl_K-270., 0.))**1.05) &
6168 & * (pres_Pa/100000.0)**.2
6171 w1 = rslf(pres_Pa,guess)
6172 w2 = rslf(pres_Pa,guess+1.)
6173 tenu = theta_e(pres_Pa,guess,w1,guess)
6174 tenup = theta_e(pres_Pa,guess+1,w2,guess+1.)
6175 cor = (thelcl_K - tenu) / (tenup - tenu)
6177 if( (cor.lt.epsilon) .and. (-cor.lt.epsilon) ) then
6182 ! print*, ' convergence not reached '
6183 thwlcl_K=theta_wetb(thelcl_K)
6184 compT_fr_The = thwlcl_K*((pres_Pa/100000.0)**0.286)
6189 !+---+-----------------------------------------------------------------+
6190 !+---+-----------------------------------------------------------------+
6191 END MODULE module_mp_thompson
6192 !+---+-----------------------------------------------------------------+