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 !+---+-----------------------------------------------------------------+
42 !wrft:model_layer:physics
43 !+---+-----------------------------------------------------------------+
45 MODULE module_mp_thompson
49 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
50 USE module_dm, ONLY : wrf_dm_max_real
52 USE module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG
56 LOGICAL, PARAMETER, PRIVATE:: iiwarm = .false.
57 LOGICAL, PRIVATE:: is_aerosol_aware = .false.
58 LOGICAL, PARAMETER, PRIVATE:: dustyIce = .true.
59 LOGICAL, PARAMETER, PRIVATE:: homogIce = .true.
61 INTEGER, PARAMETER, PRIVATE:: IFDRY = 0
62 REAL, PARAMETER, PRIVATE:: T_0 = 273.15
63 REAL, PARAMETER, PRIVATE:: PI = 3.1415926536
65 !..Densities of rain, snow, graupel, and cloud ice.
66 REAL, PARAMETER, PRIVATE:: rho_w = 1000.0
67 REAL, PARAMETER, PRIVATE:: rho_s = 100.0
68 REAL, PARAMETER, PRIVATE:: rho_g = 500.0
69 REAL, PARAMETER, PRIVATE:: rho_i = 890.0
71 !..Prescribed number of cloud droplets. Set according to known data or
72 !.. roughly 100 per cc (100.E6 m^-3) for Maritime cases and
73 !.. 300 per cc (300.E6 m^-3) for Continental. Gamma shape parameter,
74 !.. mu_c, calculated based on Nt_c is important in autoconversion
75 !.. scheme. In 2-moment cloud water, Nt_c represents a maximum of
76 !.. droplet concentration and nu_c is also variable depending on local
77 !.. droplet number concentration.
78 REAL, PARAMETER, PRIVATE:: Nt_c = 100.E6
79 REAL, PARAMETER, PRIVATE:: Nt_c_max = 1999.E6
81 !..Declaration of constants for assumed CCN/IN aerosols when none in
82 !.. the input data. Look inside the init routine for modifications
83 !.. due to surface land-sea points or vegetation characteristics.
84 REAL, PARAMETER, PRIVATE:: naIN0 = 1.5E6
85 REAL, PARAMETER, PRIVATE:: naIN1 = 0.5E6
86 REAL, PARAMETER, PRIVATE:: naCCN0 = 300.0E6
87 REAL, PARAMETER, PRIVATE:: naCCN1 = 50.0E6
88 REAL, PARAMETER, PRIVATE:: naBC0 = 150.0E6
89 REAL, PARAMETER, PRIVATE:: naBC1 = 25.0E6
91 !..Generalized gamma distributions for rain, graupel and cloud ice.
92 !.. N(D) = N_0 * D**mu * exp(-lamda*D); mu=0 is exponential.
93 REAL, PARAMETER, PRIVATE:: mu_r = 0.0
94 REAL, PARAMETER, PRIVATE:: mu_g = 0.0
95 REAL, PARAMETER, PRIVATE:: mu_i = 0.0
98 !..Sum of two gamma distrib for snow (Field et al. 2005).
99 !.. N(D) = M2**4/M3**3 * [Kap0*exp(-M2*Lam0*D/M3)
100 !.. + Kap1*(M2/M3)**mu_s * D**mu_s * exp(-M2*Lam1*D/M3)]
101 !.. M2 and M3 are the (bm_s)th and (bm_s+1)th moments respectively
102 !.. calculated as function of ice water content and temperature.
103 REAL, PARAMETER, PRIVATE:: mu_s = 0.6357
104 REAL, PARAMETER, PRIVATE:: Kap0 = 490.6
105 REAL, PARAMETER, PRIVATE:: Kap1 = 17.46
106 REAL, PARAMETER, PRIVATE:: Lam0 = 20.78
107 REAL, PARAMETER, PRIVATE:: Lam1 = 3.29
109 !..Y-intercept parameter for graupel is not constant and depends on
110 !.. mixing ratio. Also, when mu_g is non-zero, these become equiv
111 !.. y-intercept for an exponential distrib and proper values are
112 !.. computed based on same mixing ratio and total number concentration.
113 REAL, PARAMETER, PRIVATE:: gonv_min = 1.E2
114 REAL, PARAMETER, PRIVATE:: gonv_max = 1.E6
116 !..Mass power law relations: mass = am*D**bm
117 !.. Snow from Field et al. (2005), others assume spherical form.
118 REAL, PARAMETER, PRIVATE:: am_r = PI*rho_w/6.0
119 REAL, PARAMETER, PRIVATE:: bm_r = 3.0
120 REAL, PARAMETER, PRIVATE:: am_s = 0.069
121 REAL, PARAMETER, PRIVATE:: bm_s = 2.0
122 REAL, PARAMETER, PRIVATE:: am_g = PI*rho_g/6.0
123 REAL, PARAMETER, PRIVATE:: bm_g = 3.0
124 REAL, PARAMETER, PRIVATE:: am_i = PI*rho_i/6.0
125 REAL, PARAMETER, PRIVATE:: bm_i = 3.0
127 !..Fallspeed power laws relations: v = (av*D**bv)*exp(-fv*D)
128 !.. Rain from Ferrier (1994), ice, snow, and graupel from
129 !.. Thompson et al (2008). Coefficient fv is zero for graupel/ice.
130 REAL, PARAMETER, PRIVATE:: av_r = 4854.0
131 REAL, PARAMETER, PRIVATE:: bv_r = 1.0
132 REAL, PARAMETER, PRIVATE:: fv_r = 195.0
133 REAL, PARAMETER, PRIVATE:: av_s = 40.0
134 REAL, PARAMETER, PRIVATE:: bv_s = 0.55
135 REAL, PARAMETER, PRIVATE:: fv_s = 100.0
136 REAL, PARAMETER, PRIVATE:: av_g = 442.0
137 REAL, PARAMETER, PRIVATE:: bv_g = 0.89
138 REAL, PARAMETER, PRIVATE:: av_i = 1847.5
139 REAL, PARAMETER, PRIVATE:: bv_i = 1.0
140 REAL, PARAMETER, PRIVATE:: av_c = 0.316946E8
141 REAL, PARAMETER, PRIVATE:: bv_c = 2.0
143 !..Capacitance of sphere and plates/aggregates: D**3, D**2
144 REAL, PARAMETER, PRIVATE:: C_cube = 0.5
145 REAL, PARAMETER, PRIVATE:: C_sqrd = 0.15
147 !..Collection efficiencies. Rain/snow/graupel collection of cloud
148 !.. droplets use variables (Ef_rw, Ef_sw, Ef_gw respectively) and
149 !.. get computed elsewhere because they are dependent on stokes
151 REAL, PARAMETER, PRIVATE:: Ef_si = 0.05
152 REAL, PARAMETER, PRIVATE:: Ef_rs = 0.95
153 REAL, PARAMETER, PRIVATE:: Ef_rg = 0.75
154 REAL, PARAMETER, PRIVATE:: Ef_ri = 0.95
156 !..Minimum microphys values
157 !.. R1 value, 1.E-12, cannot be set lower because of numerical
158 !.. problems with Paul Field's moments and should not be set larger
159 !.. because of truncation problems in snow/ice growth.
160 REAL, PARAMETER, PRIVATE:: R1 = 1.E-12
161 REAL, PARAMETER, PRIVATE:: R2 = 1.E-6
162 REAL, PARAMETER, PRIVATE:: eps = 1.E-15
164 !..Constants in Cooper curve relation for cloud ice number.
165 REAL, PARAMETER, PRIVATE:: TNO = 5.0
166 REAL, PARAMETER, PRIVATE:: ATO = 0.304
168 !..Rho_not used in fallspeed relations (rho_not/rho)**.5 adjustment.
169 REAL, PARAMETER, PRIVATE:: rho_not = 101325.0/(287.05*298.0)
172 REAL, PARAMETER, PRIVATE:: Sc = 0.632
175 !..Homogeneous freezing temperature
176 REAL, PARAMETER, PRIVATE:: HGFR = 235.16
178 !..Water vapor and air gas constants at constant pressure
179 REAL, PARAMETER, PRIVATE:: Rv = 461.5
180 REAL, PARAMETER, PRIVATE:: oRv = 1./Rv
181 REAL, PARAMETER, PRIVATE:: R = 287.04
182 REAL, PARAMETER, PRIVATE:: Cp = 1004.0
183 REAL, PARAMETER, PRIVATE:: R_uni = 8.314 ! J (mol K)-1
185 DOUBLE PRECISION, PARAMETER, PRIVATE:: k_b = 1.38065E-23 ! Boltzmann constant [J/K]
186 DOUBLE PRECISION, PARAMETER, PRIVATE:: M_w = 18.01528E-3 ! molecular mass of water [kg/mol]
187 DOUBLE PRECISION, PARAMETER, PRIVATE:: M_a = 28.96E-3 ! molecular mass of air [kg/mol]
188 DOUBLE PRECISION, PARAMETER, PRIVATE:: N_avo = 6.022E23 ! Avogadro number [1/mol]
189 DOUBLE PRECISION, PARAMETER, PRIVATE:: ma_w = M_w / N_avo ! mass of water molecule [kg]
190 REAL, PARAMETER, PRIVATE:: ar_volume = 4./3.*PI*(2.5e-6)**3 ! assume radius of 0.025 micrometer, 2.5e-6 cm
192 !..Enthalpy of sublimation, vaporization, and fusion at 0C.
193 REAL, PARAMETER, PRIVATE:: lsub = 2.834E6
194 REAL, PARAMETER, PRIVATE:: lvap0 = 2.5E6
195 REAL, PARAMETER, PRIVATE:: lfus = lsub - lvap0
196 REAL, PARAMETER, PRIVATE:: olfus = 1./lfus
198 !..Ice initiates with this mass (kg), corresponding diameter calc.
199 !..Min diameters and mass of cloud, rain, snow, and graupel (m, kg).
200 REAL, PARAMETER, PRIVATE:: xm0i = 1.E-12
201 REAL, PARAMETER, PRIVATE:: D0c = 1.E-6
202 REAL, PARAMETER, PRIVATE:: D0r = 50.E-6
203 REAL, PARAMETER, PRIVATE:: D0s = 300.E-6
204 REAL, PARAMETER, PRIVATE:: D0g = 350.E-6
205 REAL, PRIVATE:: D0i, xm0s, xm0g
207 !..Lookup table dimensions
208 INTEGER, PARAMETER, PRIVATE:: nbins = 100
209 INTEGER, PARAMETER, PRIVATE:: nbc = nbins
210 INTEGER, PARAMETER, PRIVATE:: nbi = nbins
211 INTEGER, PARAMETER, PRIVATE:: nbr = nbins
212 INTEGER, PARAMETER, PRIVATE:: nbs = nbins
213 INTEGER, PARAMETER, PRIVATE:: nbg = nbins
214 INTEGER, PARAMETER, PRIVATE:: ntb_c = 37
215 INTEGER, PARAMETER, PRIVATE:: ntb_i = 64
216 INTEGER, PARAMETER, PRIVATE:: ntb_r = 37
217 INTEGER, PARAMETER, PRIVATE:: ntb_s = 37
218 INTEGER, PARAMETER, PRIVATE:: ntb_g = 37
219 INTEGER, PARAMETER, PRIVATE:: ntb_g1 = 37
220 INTEGER, PARAMETER, PRIVATE:: ntb_r1 = 37
221 INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55
222 INTEGER, PARAMETER, PRIVATE:: ntb_t = 9
223 INTEGER, PRIVATE:: nic1, nic2, nii2, nii3, nir2, nir3, nis2, nig2, nig3
224 INTEGER, PARAMETER, PRIVATE:: ntb_arc = 7
225 INTEGER, PARAMETER, PRIVATE:: ntb_arw = 9
226 INTEGER, PARAMETER, PRIVATE:: ntb_art = 7
227 INTEGER, PARAMETER, PRIVATE:: ntb_arr = 5
228 INTEGER, PARAMETER, PRIVATE:: ntb_ark = 4
229 INTEGER, PARAMETER, PRIVATE:: ntb_IN = 55
230 INTEGER, PRIVATE:: niIN2
232 DOUBLE PRECISION, DIMENSION(nbins+1):: xDx
233 DOUBLE PRECISION, DIMENSION(nbc):: Dc, dtc
234 DOUBLE PRECISION, DIMENSION(nbi):: Di, dti
235 DOUBLE PRECISION, DIMENSION(nbr):: Dr, dtr
236 DOUBLE PRECISION, DIMENSION(nbs):: Ds, dts
237 DOUBLE PRECISION, DIMENSION(nbg):: Dg, dtg
238 DOUBLE PRECISION, DIMENSION(nbc):: t_Nc
240 !..Lookup tables for cloud water content (kg/m**3).
241 REAL, DIMENSION(ntb_c), PARAMETER, PRIVATE:: &
242 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, &
243 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, &
244 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, &
245 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, &
248 !..Lookup tables for cloud ice content (kg/m**3).
249 REAL, DIMENSION(ntb_i), PARAMETER, PRIVATE:: &
250 r_i = (/1.e-10,2.e-10,3.e-10,4.e-10, &
251 5.e-10,6.e-10,7.e-10,8.e-10,9.e-10, &
252 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, &
253 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, &
254 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, &
255 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, &
256 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, &
257 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, &
260 !..Lookup tables for rain content (kg/m**3).
261 REAL, DIMENSION(ntb_r), PARAMETER, PRIVATE:: &
262 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, &
263 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, &
264 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, &
265 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, &
268 !..Lookup tables for graupel content (kg/m**3).
269 REAL, DIMENSION(ntb_g), PARAMETER, PRIVATE:: &
270 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, &
271 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, &
272 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, &
273 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, &
276 !..Lookup tables for snow content (kg/m**3).
277 REAL, DIMENSION(ntb_s), PARAMETER, PRIVATE:: &
278 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, &
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, &
281 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, &
284 !..Lookup tables for rain y-intercept parameter (/m**4).
285 REAL, DIMENSION(ntb_r1), PARAMETER, PRIVATE:: &
286 N0r_exp = (/1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, &
287 1.e7,2.e7,3.e7,4.e7,5.e7,6.e7,7.e7,8.e7,9.e7, &
288 1.e8,2.e8,3.e8,4.e8,5.e8,6.e8,7.e8,8.e8,9.e8, &
289 1.e9,2.e9,3.e9,4.e9,5.e9,6.e9,7.e9,8.e9,9.e9, &
292 !..Lookup tables for graupel y-intercept parameter (/m**4).
293 REAL, DIMENSION(ntb_g1), PARAMETER, PRIVATE:: &
294 N0g_exp = (/1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, &
295 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, &
296 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
297 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
300 !..Lookup tables for ice number concentration (/m**3).
301 REAL, DIMENSION(ntb_i1), PARAMETER, PRIVATE:: &
302 Nt_i = (/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0, &
303 1.e1,2.e1,3.e1,4.e1,5.e1,6.e1,7.e1,8.e1,9.e1, &
304 1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, &
305 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, &
306 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
307 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
310 !..Aerosol table parameter: Number of available aerosols, vertical
311 !.. velocity, temperature, aerosol mean radius, and hygroscopicity.
312 REAL, DIMENSION(ntb_arc), PARAMETER, PRIVATE:: &
313 ta_Na = (/10.0, 31.6, 100.0, 316.0, 1000.0, 3160.0, 10000.0/)
314 REAL, DIMENSION(ntb_arw), PARAMETER, PRIVATE:: &
315 ta_Ww = (/0.01, 0.0316, 0.1, 0.316, 1.0, 3.16, 10.0, 31.6, 100.0/)
316 REAL, DIMENSION(ntb_art), PARAMETER, PRIVATE:: &
317 ta_Tk = (/243.15, 253.15, 263.15, 273.15, 283.15, 293.15, 303.15/)
318 REAL, DIMENSION(ntb_arr), PARAMETER, PRIVATE:: &
319 ta_Ra = (/0.01, 0.02, 0.04, 0.08, 0.16/)
320 REAL, DIMENSION(ntb_ark), PARAMETER, PRIVATE:: &
321 ta_Ka = (/0.2, 0.4, 0.6, 0.8/)
323 !..Lookup tables for IN concentration (/m**3) from 0.001 to 1000/Liter.
324 REAL, DIMENSION(ntb_IN), PARAMETER, PRIVATE:: &
325 Nt_IN = (/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 !..For snow moments conversions (from Field et al. 2005)
334 REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
335 sa = (/ 5.065339, -0.062659, -3.032362, 0.029469, -0.000285, &
336 0.31255, 0.000204, 0.003199, 0.0, -0.015952/)
337 REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
338 sb = (/ 0.476221, -0.015896, 0.165977, 0.007468, -0.000141, &
339 0.060366, 0.000079, 0.000594, 0.0, -0.003577/)
341 !..Temperatures (5 C interval 0 to -40) used in lookup tables.
342 REAL, DIMENSION(ntb_t), PARAMETER, PRIVATE:: &
343 Tc = (/-0.01, -5., -10., -15., -20., -25., -30., -35., -40./)
345 !..Lookup tables for various accretion/collection terms.
346 !.. ntb_x refers to the number of elements for rain, snow, graupel,
347 !.. and temperature array indices. Variables beginning with t-p/c/m/n
348 !.. represent lookup tables. Save compile-time memory by making
349 !.. allocatable (2009Jun12, J. Michalakes).
350 INTEGER, PARAMETER, PRIVATE:: R8SIZE = 8
351 INTEGER, PARAMETER, PRIVATE:: R4SIZE = 4
352 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: &
353 tcg_racg, tmr_racg, tcr_gacr, tmg_gacr, &
355 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: &
356 tcs_racs1, tmr_racs1, tcs_racs2, tmr_racs2, &
357 tcr_sacr1, tms_sacr1, tcr_sacr2, tms_sacr2, &
358 tnr_racs1, tnr_racs2, tnr_sacr1, tnr_sacr2
359 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: &
361 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: &
362 tpi_qrfz, tpg_qrfz, tni_qrfz, tnr_qrfz
363 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: &
364 tps_iaus, tni_iaus, tpi_ide
365 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efrw
366 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efsw
367 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: tnr_rev
368 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: &
370 REAL (KIND=R4SIZE), ALLOCATABLE, DIMENSION(:,:,:,:,:):: tnccn_act
372 !..Variables holding a bunch of exponents and gamma values (cloud water,
373 !.. cloud ice, rain, snow, then graupel).
374 REAL, DIMENSION(5,15), PRIVATE:: cce, ccg
375 REAL, DIMENSION(15), PRIVATE:: ocg1, ocg2
376 REAL, DIMENSION(7), PRIVATE:: cie, cig
377 REAL, PRIVATE:: oig1, oig2, obmi
378 REAL, DIMENSION(13), PRIVATE:: cre, crg
379 REAL, PRIVATE:: ore1, org1, org2, org3, obmr
380 REAL, DIMENSION(18), PRIVATE:: cse, csg
381 REAL, PRIVATE:: oams, obms, ocms
382 REAL, DIMENSION(12), PRIVATE:: cge, cgg
383 REAL, PRIVATE:: oge1, ogg1, ogg2, ogg3, oamg, obmg, ocmg
385 !..Declaration of precomputed constants in various rate eqns.
386 REAL:: t1_qr_qc, t1_qr_qi, t2_qr_qi, t1_qg_qc, t1_qs_qc, t1_qs_qi
387 REAL:: t1_qr_ev, t2_qr_ev
388 REAL:: t1_qs_sd, t2_qs_sd, t1_qg_sd, t2_qg_sd
389 REAL:: t1_qs_me, t2_qs_me, t1_qg_me, t2_qg_me
392 !+---+-----------------------------------------------------------------+
394 !+---+-----------------------------------------------------------------+
400 SUBROUTINE thompson_init(hgt, orho, nwfa2d, nbca2d, &
401 nwfa, nifa, nbca, wif_input_opt, frc_urb2d, &
403 ids, ide, jds, jde, kds, kde, &
404 ims, ime, jms, jme, kms, kme, &
405 its, ite, jts, jte, kts, kte)
409 INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, &
410 ims,ime, jms,jme, kms,kme, &
411 its,ite, jts,jte, kts,kte
412 REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN):: hgt
414 !..OPTIONAL variables that control application of aerosol-aware scheme
416 REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(INOUT) :: nwfa, nifa, nbca
417 REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(INOUT) :: nwfa2d, nbca2d
418 REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(IN) :: orho
419 REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(IN) :: frc_urb2d
420 REAL, OPTIONAL, INTENT(IN) :: DX, DY
421 LOGICAL, OPTIONAL, INTENT(IN) :: is_start
422 INTEGER, OPTIONAL, INTENT(IN) :: wif_input_opt
423 CHARACTER*256:: mp_debug
426 INTEGER:: i, j, k, l, m, n
427 REAL:: h_01, niIN3, niCCN3, niBC3, max_test, z1
428 LOGICAL:: micro_init, has_CCN, has_IN
430 is_aerosol_aware = .FALSE.
435 write(mp_debug,*) ' DEBUG checking column of hgt ', its+1,jts+1
436 CALL wrf_debug(250, mp_debug)
438 write(mp_debug,*) ' DEBUGT k, hgt = ', k, hgt(its+1,k,jts+1)
439 CALL wrf_debug(250, mp_debug)
442 if (PRESENT(nwfa2d) .AND. PRESENT(nwfa) .AND. PRESENT(nifa)) is_aerosol_aware = .TRUE.
444 if (is_aerosol_aware) then
446 !..Check for existing aerosol data, both CCN and IN aerosols. If missing
447 !.. fill in just a basic vertical profile, somewhat boundary-layer following.
449 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
450 max_test = wrf_dm_max_real ( MAXVAL(nwfa(its:ite-1,:,jts:jte-1)) )
452 max_test = MAXVAL ( nwfa(its:ite-1,:,jts:jte-1) )
455 if (max_test .lt. eps) then
456 write(mp_debug,*) ' Apparently there are no initial CCN aerosols.'
457 CALL wrf_debug(100, mp_debug)
458 write(mp_debug,*) ' checked column at point (i,j) = ', its,jts
459 CALL wrf_debug(100, mp_debug)
460 do j = jts, min(jde-1,jte)
461 do i = its, min(ide-1,ite)
462 if (hgt(i,1,j).le.1000.0) then
464 elseif (hgt(i,1,j).ge.2500.0) then
467 h_01 = 0.8*cos(hgt(i,1,j)*0.001 - 1.0)
469 niCCN3 = -1.0*ALOG(naCCN1/naCCN0)/h_01
470 nwfa(i,1,j) = naCCN1+naCCN0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niCCN3)
471 z1=hgt(i,2,j)-hgt(i,1,j)
472 nwfa2d(i,j) = nwfa(i,1,j) * 0.000196 * (50./z1)
474 nwfa(i,k,j) = naCCN1+naCCN0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niCCN3)
480 write(mp_debug,*) ' Apparently initial CCN aerosols are present.'
481 CALL wrf_debug(100, mp_debug)
482 write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nwfa(its,:,jts))
483 CALL wrf_debug(100, mp_debug)
487 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
488 max_test = wrf_dm_max_real ( MAXVAL(nifa(its:ite-1,:,jts:jte-1)) )
490 max_test = MAXVAL ( nifa(its:ite-1,:,jts:jte-1) )
493 if (max_test .lt. eps) then
494 write(mp_debug,*) ' Apparently there are no initial IN 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 niIN3 = -1.0*ALOG(naIN1/naIN0)/h_01
508 nifa(i,1,j) = naIN1+naIN0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niIN3)
510 nifa(i,k,j) = naIN1+naIN0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niIN3)
516 write(mp_debug,*) ' Apparently initial IN aerosols are present.'
517 CALL wrf_debug(100, mp_debug)
518 write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nifa(its,:,jts))
519 CALL wrf_debug(100, mp_debug)
523 if ( wif_input_opt .eq. 2 ) then
525 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
526 max_test = wrf_dm_max_real ( MAXVAL(nbca(its:ite-1,:,jts:jte-1)) )
528 max_test = MAXVAL ( nbca(its:ite-1,:,jts:jte-1) )
531 if (max_test .lt. eps) then
532 write(mp_debug,*) ' Apparently there are no initial BC 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 niBC3 = -1.0*ALOG(naBC1/naBC0)/h_01
546 nbca(i,1,j) = naBC1+naBC0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niBC3)
547 z1=hgt(i,2,j)-hgt(i,1,j)
548 nbca2d(i,j) = nbca(i,1,j) * 0.000098 * (50./z1) * (1. + frc_urb2d(i,j))
550 nbca(i,k,j) = naBC1+naBC0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niBC3)
555 write(mp_debug,*) ' Apparently initial BC aerosols are present.'
556 CALL wrf_debug(100, mp_debug)
557 write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nbca(its,:,jts))
558 CALL wrf_debug(100, mp_debug)
566 !..Allocate space for lookup tables (J. Michalakes 2009Jun08).
568 if (.NOT. ALLOCATED(tcg_racg) ) then
569 ALLOCATE(tcg_racg(ntb_g1,ntb_g,ntb_r1,ntb_r))
573 if (.NOT. ALLOCATED(tmr_racg)) ALLOCATE(tmr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r))
574 if (.NOT. ALLOCATED(tcr_gacr)) ALLOCATE(tcr_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r))
575 if (.NOT. ALLOCATED(tmg_gacr)) ALLOCATE(tmg_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r))
576 if (.NOT. ALLOCATED(tnr_racg)) ALLOCATE(tnr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r))
577 if (.NOT. ALLOCATED(tnr_gacr)) ALLOCATE(tnr_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r))
579 if (.NOT. ALLOCATED(tcs_racs1)) ALLOCATE(tcs_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
580 if (.NOT. ALLOCATED(tmr_racs1)) ALLOCATE(tmr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
581 if (.NOT. ALLOCATED(tcs_racs2)) ALLOCATE(tcs_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
582 if (.NOT. ALLOCATED(tmr_racs2)) ALLOCATE(tmr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
583 if (.NOT. ALLOCATED(tcr_sacr1)) ALLOCATE(tcr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
584 if (.NOT. ALLOCATED(tms_sacr1)) ALLOCATE(tms_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
585 if (.NOT. ALLOCATED(tcr_sacr2)) ALLOCATE(tcr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
586 if (.NOT. ALLOCATED(tms_sacr2)) ALLOCATE(tms_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
587 if (.NOT. ALLOCATED(tnr_racs1)) ALLOCATE(tnr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
588 if (.NOT. ALLOCATED(tnr_racs2)) ALLOCATE(tnr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
589 if (.NOT. ALLOCATED(tnr_sacr1)) ALLOCATE(tnr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
590 if (.NOT. ALLOCATED(tnr_sacr2)) ALLOCATE(tnr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
592 if (.NOT. ALLOCATED(tpi_qcfz)) ALLOCATE(tpi_qcfz(ntb_c,nbc,45,ntb_IN))
593 if (.NOT. ALLOCATED(tni_qcfz)) ALLOCATE(tni_qcfz(ntb_c,nbc,45,ntb_IN))
595 if (.NOT. ALLOCATED(tpi_qrfz)) ALLOCATE(tpi_qrfz(ntb_r,ntb_r1,45,ntb_IN))
596 if (.NOT. ALLOCATED(tpg_qrfz)) ALLOCATE(tpg_qrfz(ntb_r,ntb_r1,45,ntb_IN))
597 if (.NOT. ALLOCATED(tni_qrfz)) ALLOCATE(tni_qrfz(ntb_r,ntb_r1,45,ntb_IN))
598 if (.NOT. ALLOCATED(tnr_qrfz)) ALLOCATE(tnr_qrfz(ntb_r,ntb_r1,45,ntb_IN))
600 if (.NOT. ALLOCATED(tps_iaus)) ALLOCATE(tps_iaus(ntb_i,ntb_i1))
601 if (.NOT. ALLOCATED(tni_iaus)) ALLOCATE(tni_iaus(ntb_i,ntb_i1))
602 if (.NOT. ALLOCATED(tpi_ide)) ALLOCATE(tpi_ide(ntb_i,ntb_i1))
604 if (.NOT. ALLOCATED(t_Efrw)) ALLOCATE(t_Efrw(nbr,nbc))
605 if (.NOT. ALLOCATED(t_Efsw)) ALLOCATE(t_Efsw(nbs,nbc))
607 if (.NOT. ALLOCATED(tnr_rev)) ALLOCATE(tnr_rev(nbr, ntb_r1, ntb_r))
608 if (.NOT. ALLOCATED(tpc_wev)) ALLOCATE(tpc_wev(nbc,ntb_c,nbc))
609 if (.NOT. ALLOCATED(tnc_wev)) ALLOCATE(tnc_wev(nbc,ntb_c,nbc))
611 if (.NOT. ALLOCATED(tnccn_act)) &
612 ALLOCATE(tnccn_act(ntb_arc,ntb_arw,ntb_art,ntb_arr,ntb_ark))
616 !..From Martin et al. (1994), assign gamma shape parameter mu for cloud
617 !.. drops according to general dispersion characteristics (disp=~0.25
618 !.. for Maritime and 0.45 for Continental).
619 !.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime
620 !.. to 2 for really dirty air. This not used in 2-moment cloud water
621 !.. scheme and nu_c used instead and varies from 2 to 15 (integer-only).
622 mu_c = MIN(15., (1000.E6/Nt_c + 2.))
624 !..Schmidt number to one-third used numerous times.
627 !..Compute min ice diam from mass, min snow/graupel mass from diam.
628 D0i = (xm0i/am_i)**(1./bm_i)
629 xm0s = am_s * D0s**bm_s
630 xm0g = am_g * D0g**bm_g
632 !..These constants various exponents and gamma() assoc with cloud,
633 !.. rain, snow, and graupel.
636 cce(2,n) = bm_r + n + 1.
637 cce(3,n) = bm_r + n + 4.
638 cce(4,n) = n + bv_c + 1.
639 cce(5,n) = bm_r + n + bv_c + 1.
640 ccg(1,n) = WGAMMA(cce(1,n))
641 ccg(2,n) = WGAMMA(cce(2,n))
642 ccg(3,n) = WGAMMA(cce(3,n))
643 ccg(4,n) = WGAMMA(cce(4,n))
644 ccg(5,n) = WGAMMA(cce(5,n))
645 ocg1(n) = 1./ccg(1,n)
646 ocg2(n) = 1./ccg(2,n)
650 cie(2) = bm_i + mu_i + 1.
651 cie(3) = bm_i + mu_i + bv_i + 1.
652 cie(4) = mu_i + bv_i + 1.
654 cie(6) = bm_i*0.5 + mu_i + bv_i + 1.
655 cie(7) = bm_i*0.5 + mu_i + 1.
656 cig(1) = WGAMMA(cie(1))
657 cig(2) = WGAMMA(cie(2))
658 cig(3) = WGAMMA(cie(3))
659 cig(4) = WGAMMA(cie(4))
660 cig(5) = WGAMMA(cie(5))
661 cig(6) = WGAMMA(cie(6))
662 cig(7) = WGAMMA(cie(7))
669 cre(3) = bm_r + mu_r + 1.
670 cre(4) = bm_r*2. + mu_r + 1.
671 cre(5) = mu_r + bv_r + 1.
672 cre(6) = bm_r + mu_r + bv_r + 1.
673 cre(7) = bm_r*0.5 + mu_r + bv_r + 1.
674 cre(8) = bm_r + mu_r + bv_r + 3.
675 cre(9) = mu_r + bv_r + 3.
677 cre(11) = 0.5*(bv_r + 5. + 2.*mu_r)
678 cre(12) = bm_r*0.5 + mu_r + 1.
679 cre(13) = bm_r*2. + mu_r + bv_r + 1.
681 crg(n) = WGAMMA(cre(n))
692 cse(4) = bm_s + bv_s + 1.
693 cse(5) = bm_s*2. + bv_s + 1.
694 cse(6) = bm_s*2. + 1.
695 cse(7) = bm_s + mu_s + 1.
696 cse(8) = bm_s + mu_s + 2.
697 cse(9) = bm_s + mu_s + 3.
698 cse(10) = bm_s + mu_s + bv_s + 1.
699 cse(11) = bm_s*2. + mu_s + bv_s + 1.
700 cse(12) = bm_s*2. + mu_s + 1.
702 cse(14) = bm_s + bv_s
704 cse(16) = 1.0 + (1.0 + bv_s)/2.
705 cse(17) = cse(16) + mu_s + 1.
706 cse(18) = bv_s + mu_s + 3.
708 csg(n) = WGAMMA(cse(n))
716 cge(3) = bm_g + mu_g + 1.
717 cge(4) = bm_g*2. + mu_g + 1.
718 cge(5) = bm_g*2. + mu_g + bv_g + 1.
719 cge(6) = bm_g + mu_g + bv_g + 1.
720 cge(7) = bm_g + mu_g + bv_g + 2.
721 cge(8) = bm_g + mu_g + bv_g + 3.
722 cge(9) = mu_g + bv_g + 3.
724 cge(11) = 0.5*(bv_g + 5. + 2.*mu_g)
725 cge(12) = 0.5*(bv_g + 5.) + mu_g
727 cgg(n) = WGAMMA(cge(n))
737 !+---+-----------------------------------------------------------------+
738 !..Simplify various rate eqns the best we can now.
739 !+---+-----------------------------------------------------------------+
741 !..Rain collecting cloud water and cloud ice
742 t1_qr_qc = PI*.25*av_r * crg(9)
743 t1_qr_qi = PI*.25*av_r * crg(9)
744 t2_qr_qi = PI*.25*am_r*av_r * crg(8)
746 !..Graupel collecting cloud water
747 t1_qg_qc = PI*.25*av_g * cgg(9)
749 !..Snow collecting cloud water
750 t1_qs_qc = PI*.25*av_s
752 !..Snow collecting cloud ice
753 t1_qs_qi = PI*.25*av_s
755 !..Evaporation of rain; ignore depositional growth of rain.
756 t1_qr_ev = 0.78 * crg(10)
757 t2_qr_ev = 0.308*Sc3*SQRT(av_r) * crg(11)
759 !..Sublimation/depositional growth of snow
761 t2_qs_sd = 0.28*Sc3*SQRT(av_s)
764 t1_qs_me = PI*4.*C_sqrd*olfus * 0.86
765 t2_qs_me = PI*4.*C_sqrd*olfus * 0.28*Sc3*SQRT(av_s)
767 !..Sublimation/depositional growth of graupel
768 t1_qg_sd = 0.86 * cgg(10)
769 t2_qg_sd = 0.28*Sc3*SQRT(av_g) * cgg(11)
771 !..Melting of graupel
772 t1_qg_me = PI*4.*C_cube*olfus * 0.86 * cgg(10)
773 t2_qg_me = PI*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g) * cgg(11)
775 !..Constants for helping find lookup table indexes.
776 nic2 = NINT(ALOG10(r_c(1)))
777 nii2 = NINT(ALOG10(r_i(1)))
778 nii3 = NINT(ALOG10(Nt_i(1)))
779 nir2 = NINT(ALOG10(r_r(1)))
780 nir3 = NINT(ALOG10(N0r_exp(1)))
781 nis2 = NINT(ALOG10(r_s(1)))
782 nig2 = NINT(ALOG10(r_g(1)))
783 nig3 = NINT(ALOG10(N0g_exp(1)))
784 niIN2 = NINT(ALOG10(Nt_IN(1)))
786 !..Create bins of cloud water (from min diameter up to 100 microns).
790 Dc(n) = Dc(n-1) + 1.0D-6
791 dtc(n) = (Dc(n) - Dc(n-1))
794 !..Create bins of cloud ice (from min diameter up to 5x min snow size).
796 xDx(nbi+1) = 5.0d0*D0s
798 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbi) &
799 *DLOG(xDx(nbi+1)/xDx(1)) +DLOG(xDx(1)))
802 Di(n) = DSQRT(xDx(n)*xDx(n+1))
803 dti(n) = xDx(n+1) - xDx(n)
806 !..Create bins of rain (from min diameter up to 5 mm).
810 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbr) &
811 *DLOG(xDx(nbr+1)/xDx(1)) +DLOG(xDx(1)))
814 Dr(n) = DSQRT(xDx(n)*xDx(n+1))
815 dtr(n) = xDx(n+1) - xDx(n)
818 !..Create bins of snow (from min diameter up to 2 cm).
822 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbs) &
823 *DLOG(xDx(nbs+1)/xDx(1)) +DLOG(xDx(1)))
826 Ds(n) = DSQRT(xDx(n)*xDx(n+1))
827 dts(n) = xDx(n+1) - xDx(n)
830 !..Create bins of graupel (from min diameter up to 5 cm).
834 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbg) &
835 *DLOG(xDx(nbg+1)/xDx(1)) +DLOG(xDx(1)))
838 Dg(n) = DSQRT(xDx(n)*xDx(n+1))
839 dtg(n) = xDx(n+1) - xDx(n)
842 !..Create bins of cloud droplet number concentration (1 to 3000 per cc).
844 xDx(nbc+1) = 3000.0d0
846 xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbc) &
847 *DLOG(xDx(nbc+1)/xDx(1)) +DLOG(xDx(1)))
850 t_Nc(n) = DSQRT(xDx(n)*xDx(n+1)) * 1.D6
852 nic1 = DLOG(t_Nc(nbc)/t_Nc(1))
854 !+---+-----------------------------------------------------------------+
855 !..Create lookup tables for most costly calculations.
856 !+---+-----------------------------------------------------------------+
862 tcg_racg(i,j,k,m) = 0.0d0
863 tmr_racg(i,j,k,m) = 0.0d0
864 tcr_gacr(i,j,k,m) = 0.0d0
865 tmg_gacr(i,j,k,m) = 0.0d0
866 tnr_racg(i,j,k,m) = 0.0d0
867 tnr_gacr(i,j,k,m) = 0.0d0
877 tcs_racs1(i,j,k,m) = 0.0d0
878 tmr_racs1(i,j,k,m) = 0.0d0
879 tcs_racs2(i,j,k,m) = 0.0d0
880 tmr_racs2(i,j,k,m) = 0.0d0
881 tcr_sacr1(i,j,k,m) = 0.0d0
882 tms_sacr1(i,j,k,m) = 0.0d0
883 tcr_sacr2(i,j,k,m) = 0.0d0
884 tms_sacr2(i,j,k,m) = 0.0d0
885 tnr_racs1(i,j,k,m) = 0.0d0
886 tnr_racs2(i,j,k,m) = 0.0d0
887 tnr_sacr1(i,j,k,m) = 0.0d0
888 tnr_sacr2(i,j,k,m) = 0.0d0
898 tpi_qrfz(i,j,k,m) = 0.0d0
899 tni_qrfz(i,j,k,m) = 0.0d0
900 tpg_qrfz(i,j,k,m) = 0.0d0
901 tnr_qrfz(i,j,k,m) = 0.0d0
906 tpi_qcfz(i,j,k,m) = 0.0d0
907 tni_qcfz(i,j,k,m) = 0.0d0
915 tps_iaus(i,j) = 0.0d0
916 tni_iaus(i,j) = 0.0d0
933 tnr_rev(i,j,k) = 0.0d0
941 tpc_wev(i,j,k) = 0.0d0
942 tnc_wev(i,j,k) = 0.0d0
952 tnccn_act(i,j,k,l,m) = 1.0
959 CALL wrf_debug(150, 'CREATING MICROPHYSICS LOOKUP TABLES ... ')
960 WRITE (wrf_err_message, '(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') &
961 ' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g
962 CALL wrf_debug(150, wrf_err_message)
964 !..Read a static file containing CCN activation of aerosols. The
965 !.. data were created from a parcel model by Feingold & Heymsfield with
966 !.. further changes by Eidhammer and Kriedenweis.
967 if (is_aerosol_aware) then
968 CALL wrf_debug(200, ' calling table_ccnAct routine')
972 !..Collision efficiency between rain/snow and cloud water.
973 CALL wrf_debug(200, ' creating qc collision eff tables')
978 CALL wrf_debug(200, ' creating rain evap table')
981 !..Initialize various constants for computing radar reflectivity.
993 if (.not. iiwarm) then
995 !..Rain collecting graupel & graupel collecting rain.
996 CALL wrf_debug(200, ' creating rain collecting graupel table')
999 !..Rain collecting snow & snow collecting rain.
1000 CALL wrf_debug(200, ' creating rain collecting snow table')
1003 !..Cloud water and rain freezing (Bigg, 1953).
1004 CALL wrf_debug(200, ' creating freezing of water drops table')
1007 !..Conversion of some ice mass into snow category.
1008 CALL wrf_debug(200, ' creating ice converting to snow table')
1013 CALL wrf_debug(150, ' ... DONE microphysical lookup tables')
1017 END SUBROUTINE thompson_init
1018 !+---+-----------------------------------------------------------------+
1020 !+---+-----------------------------------------------------------------+
1021 !..This is a wrapper routine designed to transfer values from 3D to 1D.
1022 !+---+-----------------------------------------------------------------+
1023 SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
1024 nwfa, nifa, nbca, nwfa2d, nifa2d, nbca2d, &
1025 aer_init_opt, wif_input_opt, &
1026 th, pii, p, w, dz, dt_in, itimestep, &
1029 GRAUPELNC, GRAUPELNCV, SR, &
1030 #if ( WRF_CHEM == 1 )
1031 wetscav_on, rainprod, evapprod, &
1033 refl_10cm, diagflag, ke_diag, do_radar_ref, &
1034 re_cloud, re_ice, re_snow, &
1035 has_reqc, has_reqi, has_reqs, &
1036 ids,ide, jds,jde, kds,kde, & ! domain dims
1037 ims,ime, jms,jme, kms,kme, & ! memory dims
1038 its,ite, jts,jte, kts,kte) ! tile dims
1042 !..Subroutine arguments
1043 INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, &
1044 ims,ime, jms,jme, kms,kme, &
1045 its,ite, jts,jte, kts,kte
1046 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
1047 qv, qc, qr, qi, qs, qg, ni, nr, th
1048 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: &
1049 nc, nwfa, nifa, nbca
1050 REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d, &
1052 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
1053 re_cloud, re_ice, re_snow
1054 INTEGER, INTENT(IN):: has_reqc, has_reqi, has_reqs
1055 #if ( WRF_CHEM == 1 )
1056 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
1059 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
1061 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: &
1063 REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(INOUT):: &
1064 SNOWNC, SNOWNCV, GRAUPELNC, GRAUPELNCV
1065 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
1067 REAL, INTENT(IN):: dt_in
1068 INTEGER, INTENT(IN):: itimestep
1069 INTEGER, OPTIONAL, INTENT(IN):: aer_init_opt, wif_input_opt
1070 #if ( WRF_CHEM == 1 )
1071 LOGICAL, INTENT(in) :: wetscav_on
1075 REAL, DIMENSION(kts:kte):: &
1076 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
1077 nr1d, nc1d, nwfa1d, nifa1d, nbca1d, &
1078 t1d, p1d, w1d, dz1d, rho, dBZ
1079 REAL, DIMENSION(kts:kte):: re_qc1d, re_qi1d, re_qs1d
1080 #if ( WRF_CHEM == 1 )
1081 REAL, DIMENSION(kts:kte):: &
1082 rainprod1d, evapprod1d
1084 REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
1085 REAL:: dt, pptrain, pptsnow, pptgraul, pptice
1086 REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max
1089 INTEGER:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr
1090 INTEGER:: jmax_qc,jmax_qr,jmax_qi,jmax_qs,jmax_qg,jmax_ni,jmax_nr
1091 INTEGER:: kmax_qc,kmax_qr,kmax_qi,kmax_qs,kmax_qg,kmax_ni,kmax_nr
1092 INTEGER:: i_start, j_start, i_end, j_end
1093 LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
1094 INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref, ke_diag
1095 CHARACTER*256:: mp_debug
1097 integer :: kediagloc
1103 i_end = MIN(ite, ide-1)
1104 j_end = MIN(jte, jde-1)
1106 !..For idealized testing by developer.
1107 ! if ( (ide-ids+1).gt.4 .and. (jde-jds+1).lt.4 .and. &
1108 ! ids.eq.its.and.ide.eq.ite.and.jds.eq.jts.and.jde.eq.jte) then
1146 mp_debug(i:i) = char(0)
1149 if (.NOT. is_aerosol_aware .AND. PRESENT(nc) .AND. PRESENT(nwfa) &
1150 .AND. PRESENT(nifa) .AND. PRESENT(nwfa2d)) then
1151 write(mp_debug,*) 'WARNING, nc-nwfa-nifa-nwfa2d present but is_aerosol_aware is FALSE'
1152 CALL wrf_debug(0, mp_debug)
1155 j_loop: do j = j_start, j_end
1156 i_loop: do i = i_start, i_end
1163 IF ( PRESENT (snowncv) ) THEN
1166 IF ( PRESENT (graupelncv) ) THEN
1167 GRAUPELNCV(i,j) = 0.
1172 t1d(k) = th(i,k,j)*pii(i,k,j)
1184 rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622))
1186 if (is_aerosol_aware) then
1189 nwfa1d(k) = nwfa(i,k,j)
1190 nifa1d(k) = nifa(i,k,j)
1191 if (wif_input_opt .eq. 2) then
1192 nbca1d(k) = nbca(i,k,j)
1200 nc1d(k) = Nt_c/rho(k)
1201 nwfa1d(k) = 11.1E6/rho(k)
1202 nifa1d(k) = naIN1*0.01/rho(k)
1203 nbca1d(k) = 5.55E6/rho(k)
1208 call mp_thompson(aer_init_opt, wif_input_opt, &
1209 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
1210 nr1d, nc1d, nwfa1d, nifa1d, nbca1d, t1d, p1d, w1d, dz1d, &
1211 pptrain, pptsnow, pptgraul, pptice, &
1212 #if ( WRF_CHEM == 1 )
1213 wetscav_on, rainprod1d, evapprod1d, &
1217 pcp_ra(i,j) = pptrain
1218 pcp_sn(i,j) = pptsnow
1219 pcp_gr(i,j) = pptgraul
1220 pcp_ic(i,j) = pptice
1221 RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice
1222 RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
1223 IF ( PRESENT(snowncv) .AND. PRESENT(snownc) ) THEN
1224 SNOWNCV(i,j) = pptsnow + pptice
1225 SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + pptice
1227 IF ( PRESENT(graupelncv) .AND. PRESENT(graupelnc) ) THEN
1228 GRAUPELNCV(i,j) = pptgraul
1229 GRAUPELNC(i,j) = GRAUPELNC(i,j) + pptgraul
1231 SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12)
1235 !..BEGIN AEROSOL EMISSIONS
1237 !..Reset lowest model level to initial state aerosols (fake sfc source).
1238 !.. Changed 13 May 2013 to fake emissions in which nwfa2d is aerosol
1239 !.. number tendency (number per kg per second).
1240 if (is_aerosol_aware) then
1241 !..Add anthropogenic emissions
1242 !-GT nwfa1d(kts) = nwfa1
1243 nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt_in
1244 nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt_in
1245 if (wif_input_opt .eq. 2) then
1246 nbca1d(kts) = nbca1d(kts) + nbca2d(i,j)*dt_in
1251 !..END AEROSOL EMISSIONS
1256 nwfa(i,k,j) = nwfa1d(k)
1257 nifa(i,k,j) = nifa1d(k)
1258 if (wif_input_opt .eq. 2) then
1259 nbca(i,k,j) = nbca1d(k)
1275 th(i,k,j) = t1d(k)/pii(i,k,j)
1276 #if ( WRF_CHEM == 1 )
1277 IF ( wetscav_on ) THEN
1278 rainprod(i,k,j) = rainprod1d(k)
1279 evapprod(i,k,j) = evapprod1d(k)
1282 if (qc1d(k) .gt. qc_max) then
1287 elseif (qc1d(k) .lt. 0.0) then
1288 write(mp_debug,*) 'WARNING, negative qc ', qc1d(k), &
1290 CALL wrf_debug(150, mp_debug)
1292 if (qr1d(k) .gt. qr_max) then
1297 elseif (qr1d(k) .lt. 0.0) then
1298 write(mp_debug,*) 'WARNING, negative qr ', qr1d(k), &
1300 CALL wrf_debug(150, mp_debug)
1302 if (nr1d(k) .gt. nr_max) then
1307 elseif (nr1d(k) .lt. 0.0) then
1308 write(mp_debug,*) 'WARNING, negative nr ', nr1d(k), &
1310 CALL wrf_debug(150, mp_debug)
1312 if (qs1d(k) .gt. qs_max) then
1317 elseif (qs1d(k) .lt. 0.0) then
1318 write(mp_debug,*) 'WARNING, negative qs ', qs1d(k), &
1320 CALL wrf_debug(150, mp_debug)
1322 if (qi1d(k) .gt. qi_max) then
1327 elseif (qi1d(k) .lt. 0.0) then
1328 write(mp_debug,*) 'WARNING, negative qi ', qi1d(k), &
1330 CALL wrf_debug(150, mp_debug)
1332 if (qg1d(k) .gt. qg_max) then
1337 elseif (qg1d(k) .lt. 0.0) then
1338 write(mp_debug,*) 'WARNING, negative qg ', qg1d(k), &
1340 CALL wrf_debug(150, mp_debug)
1342 if (ni1d(k) .gt. ni_max) then
1347 elseif (ni1d(k) .lt. 0.0) then
1348 write(mp_debug,*) 'WARNING, negative ni ', ni1d(k), &
1350 CALL wrf_debug(150, mp_debug)
1352 if (qv1d(k) .lt. 0.0) then
1353 write(mp_debug,*) 'WARNING, negative qv ', qv1d(k), &
1355 CALL wrf_debug(150, mp_debug)
1356 if (k.lt.kte-2 .and. k.gt.kts+1) then
1357 write(mp_debug,*) ' below and above are: ', qv(i,k-1,j), qv(i,k+1,j)
1358 CALL wrf_debug(150, mp_debug)
1359 qv(i,k,j) = MAX(1.E-7, 0.5*(qv(i,k-1,j) + qv(i,k+1,j)))
1366 IF ( PRESENT (diagflag) ) THEN
1367 if (diagflag .and. do_radar_ref == 1) then
1369 IF ( present(ke_diag) ) THEN
1375 call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, &
1376 t1d, p1d, dBZ, kts, kte, i, j, kediagloc)
1378 refl_10cm(i,k,j) = MAX(-35., dBZ(k))
1383 IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN
1385 re_qc1d(k) = RE_QC_BG
1386 re_qi1d(k) = RE_QI_BG
1387 re_qs1d(k) = RE_QS_BG
1389 call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
1390 re_qc1d, re_qi1d, re_qs1d, kts, kte)
1392 re_cloud(i,k,j) = MAX(RE_QC_BG, MIN(re_qc1d(k), 50.E-6))
1393 re_ice(i,k,j) = MAX(RE_QI_BG, MIN(re_qi1d(k), 125.E-6))
1394 re_snow(i,k,j) = MAX(RE_QS_BG, MIN(re_qs1d(k), 999.E-6))
1402 write(mp_debug,'(a,7(a,e13.6,1x,a,i3,a,i3,a,i3,a,1x))') 'MP-GT:', &
1403 'qc: ', qc_max, '(', imax_qc, ',', jmax_qc, ',', kmax_qc, ')', &
1404 'qr: ', qr_max, '(', imax_qr, ',', jmax_qr, ',', kmax_qr, ')', &
1405 'qi: ', qi_max, '(', imax_qi, ',', jmax_qi, ',', kmax_qi, ')', &
1406 'qs: ', qs_max, '(', imax_qs, ',', jmax_qs, ',', kmax_qs, ')', &
1407 'qg: ', qg_max, '(', imax_qg, ',', jmax_qg, ',', kmax_qg, ')', &
1408 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')', &
1409 'nr: ', nr_max, '(', imax_nr, ',', jmax_nr, ',', kmax_nr, ')'
1410 CALL wrf_debug(150, mp_debug)
1414 mp_debug(i:i) = char(0)
1417 END SUBROUTINE mp_gt_driver
1419 !+---+-----------------------------------------------------------------+
1421 !+---+-----------------------------------------------------------------+
1422 !+---+-----------------------------------------------------------------+
1423 !.. This subroutine computes the moisture tendencies of water vapor,
1424 !.. cloud droplets, rain, cloud ice (pristine), snow, and graupel.
1425 !.. Previously this code was based on Reisner et al (1998), but few of
1426 !.. those pieces remain. A complete description is now found in
1427 !.. Thompson et al. (2004, 2008).
1428 !+---+-----------------------------------------------------------------+
1430 subroutine mp_thompson (aer_init_opt, wif_input_opt, &
1431 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
1432 nr1d, nc1d, nwfa1d, nifa1d, nbca1d, t1d, p1d, w1d, dzq, &
1433 pptrain, pptsnow, pptgraul, pptice, &
1434 #if ( WRF_CHEM == 1 )
1435 wetscav_on, rainprod, evapprod, &
1437 kts, kte, dt, ii, jj)
1442 INTEGER, OPTIONAL, INTENT(IN):: aer_init_opt, wif_input_opt
1443 INTEGER, INTENT(IN):: kts, kte, ii, jj
1444 REAL, DIMENSION(kts:kte), INTENT(INOUT):: &
1445 qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
1446 nr1d, nc1d, nwfa1d, nifa1d, nbca1d, t1d
1447 REAL, DIMENSION(kts:kte), INTENT(IN):: p1d, w1d, dzq
1448 REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice
1449 REAL, INTENT(IN):: dt
1450 #if ( WRF_CHEM == 1 )
1451 REAL, DIMENSION(kts:kte), INTENT(INOUT):: &
1453 LOGICAL, INTENT(IN) :: wetscav_on
1457 REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, &
1458 qrten, qsten, qgten, niten, nrten, ncten, nwfaten, nifaten, &
1461 DOUBLE PRECISION, DIMENSION(kts:kte):: prw_vcd
1463 DOUBLE PRECISION, DIMENSION(kts:kte):: pnc_wcd, pnc_wau, pnc_rcw, &
1466 DOUBLE PRECISION, DIMENSION(kts:kte):: pna_rca, pna_sca, pna_gca, &
1467 pnd_rcd, pnd_scd, pnd_gcd, pnb_rcb, pnb_scb, pnb_gcb
1469 DOUBLE PRECISION, DIMENSION(kts:kte):: prr_wau, prr_rcw, prr_rcs, &
1470 prr_rcg, prr_sml, prr_gml, &
1472 pnr_wau, pnr_rcs, pnr_rcg, &
1473 pnr_rci, pnr_sml, pnr_gml, &
1474 pnr_rev, pnr_rcr, pnr_rfz
1476 DOUBLE PRECISION, DIMENSION(kts:kte):: pri_inu, pni_inu, pri_ihm, &
1477 pni_ihm, pri_wfz, pni_wfz, &
1478 pri_rfz, pni_rfz, pri_ide, &
1479 pni_ide, pri_rci, pni_rci, &
1480 pni_sci, pni_iau, pri_iha, pni_iha
1482 DOUBLE PRECISION, DIMENSION(kts:kte):: prs_iau, prs_sci, prs_rcs, &
1483 prs_scw, prs_sde, prs_ihm, &
1486 DOUBLE PRECISION, DIMENSION(kts:kte):: prg_scw, prg_rfz, prg_gde, &
1487 prg_gcw, prg_rci, prg_rcs, &
1490 DOUBLE PRECISION, PARAMETER:: zeroD0 = 0.0d0
1492 REAL, DIMENSION(kts:kte):: temp, pres, qv
1493 REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni, nr, nc, nwfa, nifa, nbca
1494 REAL, DIMENSION(kts:kte):: rho, rhof, rhof2
1495 REAL, DIMENSION(kts:kte):: qvs, qvsi, delQvs
1496 REAL, DIMENSION(kts:kte):: satw, sati, ssatw, ssati
1497 REAL, DIMENSION(kts:kte):: diffu, visco, vsc2, &
1498 tcond, lvap, ocp, lvt2
1500 DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g
1501 REAL, DIMENSION(kts:kte):: mvd_r, mvd_c
1502 REAL, DIMENSION(kts:kte):: smob, smo2, smo1, smo0, &
1503 smoc, smod, smoe, smof
1505 REAL, DIMENSION(kts:kte):: sed_r, sed_s, sed_g, sed_i, sed_n,sed_c
1507 REAL:: rgvm, delta_tp, orho, lfus2
1508 REAL, DIMENSION(5):: onstep
1509 DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamc, lamr, lamg
1510 DOUBLE PRECISION:: lami, ilami, ilamc
1511 REAL:: xDc, Dc_b, Dc_g, xDi, xDr, xDs, xDg, Ds_m, Dg_m
1512 DOUBLE PRECISION:: Dr_star, Dc_star
1513 REAL:: zeta1, zeta, taud, tau
1514 REAL:: stoke_r, stoke_s, stoke_g, stoke_i
1515 REAL:: vti, vtr, vts, vtg, vtc
1516 REAL, DIMENSION(kts:kte+1):: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk, &
1518 REAL, DIMENSION(kts:kte):: vts_boost
1519 REAL:: Mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, C_snow
1520 REAL:: a_, b_, loga_, A1, A2, tf
1521 REAL:: tempc, tc0, r_mvd1, r_mvd2, xkrat
1522 REAL:: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr
1523 REAL:: xsat, rate_max, sump, ratio
1524 REAL:: clap, fcd, dfcd
1525 REAL:: otemp, rvs, rvs_p, rvs_pp, gamsc, alphsc, t1_evap, t1_subl
1526 REAL:: r_frac, g_frac
1527 REAL:: Ef_rw, Ef_sw, Ef_gw, Ef_rr
1528 REAL:: Ef_ra, Ef_sa, Ef_ga
1529 REAL:: dtsave, odts, odt, odzq, hgt_agl, SR
1530 REAL:: xslw1, ygra1, zans1, eva_factor
1531 INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq
1532 INTEGER, DIMENSION(5):: ksed1
1533 INTEGER:: nir, nis, nig, nii, nic, niin
1534 INTEGER:: idx_tc, idx_t, idx_s, idx_g1, idx_g, idx_r1, idx_r, &
1535 idx_i1, idx_i, idx_c, idx, idx_d, idx_n, idx_in
1537 LOGICAL:: melti, no_micro
1538 LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg
1539 LOGICAL:: debug_flag
1540 CHARACTER*256:: mp_debug
1545 debug_flag = .false.
1546 ! if (ii.eq.901 .and. jj.eq.379) debug_flag = .true.
1548 write(mp_debug, *) 'DEBUG INFO, mp_thompson at (i,j) ', ii, ', ', jj
1549 CALL wrf_debug(550, mp_debug)
1558 !+---+-----------------------------------------------------------------+
1559 !.. Source/sink terms. First 2 chars: "pr" represents source/sink of
1560 !.. mass while "pn" represents source/sink of number. Next char is one
1561 !.. of "v" for water vapor, "r" for rain, "i" for cloud ice, "w" for
1562 !.. cloud water, "s" for snow, and "g" for graupel. Next chars
1563 !.. represent processes: "de" for sublimation/deposition, "ev" for
1564 !.. evaporation, "fz" for freezing, "ml" for melting, "au" for
1565 !.. autoconversion, "nu" for ice nucleation, "hm" for Hallet/Mossop
1566 !.. secondary ice production, and "c" for collection followed by the
1567 !.. character for the species being collected. ALL of these terms are
1568 !.. positive (except for deposition/sublimation terms which can switch
1569 !.. signs based on super/subsaturation) and are treated as negatives
1570 !.. where necessary in the tendency equations.
1571 !+---+-----------------------------------------------------------------+
1660 #if ( WRF_CHEM == 1 )
1661 if ( wetscav_on ) then
1669 !..Bug fix (2016Jun15), prevent use of uninitialized value(s).
1683 !+---+-----------------------------------------------------------------+
1684 !..Put column of data into local arrays.
1685 !+---+-----------------------------------------------------------------+
1688 qv(k) = MAX(1.E-10, qv1d(k))
1690 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
1691 if (is_aerosol_aware) then
1692 if (aer_init_opt .lt. 2) then ! Constant or climatology (e.g., GOCART)
1693 nwfa(k) = MAX(11.1E6, MIN(9999.E6, nwfa1d(k)*rho(k)))
1694 nifa(k) = MAX(naIN1*0.01, MIN(9999.E6, nifa1d(k)*rho(k)))
1695 if (wif_input_opt .eq. 2) then ! Considering BC aerosol
1696 nbca(k) = MAX(5.55E6, MIN(9999.E6, nbca1d(k)*rho(k)))
1700 else ! First guess (e.g., GEOS-5)
1701 nwfa(k) = MAX(0.0, nwfa1d(k)*rho(k))
1702 nifa(k) = MAX(0.0, nifa1d(k)*rho(k))
1703 if (wif_input_opt .eq. 2) then ! Considering BC aerosol
1704 nbca(k) = MAX(0.0, nbca1d(k)*rho(k))
1710 nwfa(k) = MAX(11.1E6, MIN(9999.E6, nwfa1d(k)*rho(k)))
1711 nifa(k) = MAX(naIN1*0.01, MIN(9999.E6, nifa1d(k)*rho(k)))
1712 nbca(k) = MAX(5.55E6, MIN(9999.E6, nbca1d(k)*rho(k)))
1715 if (qc1d(k) .gt. R1) then
1717 rc(k) = qc1d(k)*rho(k)
1718 nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max))
1720 nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2)
1721 lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
1722 xDc = (bm_r + nu_c + 1.) / lamc
1723 if (xDc.lt. D0c) then
1724 lamc = cce(2,nu_c)/D0c
1725 elseif (xDc.gt. D0r*2.) then
1726 lamc = cce(2,nu_c)/(D0r*2.)
1728 nc(k) = MIN( DBLE(Nt_c_max), ccg(1,nu_c)*ocg2(nu_c)*rc(k) &
1730 if (.NOT. is_aerosol_aware) nc(k) = Nt_c
1739 if (qi1d(k) .gt. R1) then
1741 ri(k) = qi1d(k)*rho(k)
1742 ni(k) = MAX(R2, ni1d(k)*rho(k))
1743 if (ni(k).le. R2) then
1745 ni(k) = MIN(999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
1748 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
1750 xDi = (bm_i + mu_i + 1.) * ilami
1751 if (xDi.lt. 5.E-6) then
1753 ni(k) = MIN(999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
1754 elseif (xDi.gt. 300.E-6) then
1755 lami = cie(2)/300.E-6
1756 ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i
1766 if (qr1d(k) .gt. R1) then
1768 rr(k) = qr1d(k)*rho(k)
1769 nr(k) = MAX(R2, nr1d(k)*rho(k))
1770 if (nr(k).le. R2) then
1772 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
1773 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
1776 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
1777 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
1778 if (mvd_r(k) .gt. 2.5E-3) then
1780 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
1781 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
1782 elseif (mvd_r(k) .lt. D0r*0.75) then
1784 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
1785 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
1794 if (qs1d(k) .gt. R1) then
1796 rs(k) = qs1d(k)*rho(k)
1803 if (qg1d(k) .gt. R1) then
1805 rg(k) = qg1d(k)*rho(k)
1814 !+---+-----------------------------------------------------------------+
1815 ! if (debug_flag) then
1816 ! write(mp_debug,*) 'DEBUG-VERBOSE at (i,j) ', ii, ', ', jj
1817 ! CALL wrf_debug(550, mp_debug)
1819 ! write(mp_debug, '(a,i3,f8.2,1x,f7.2,1x, 11(1x,e13.6))') &
1820 ! & '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), nwfa(k), nifa(k)
1821 ! CALL wrf_debug(550, mp_debug)
1824 !+---+-----------------------------------------------------------------+
1826 !+---+-----------------------------------------------------------------+
1827 !..Derive various thermodynamic variables frequently used.
1828 !.. Saturation vapor pressure (mixing ratio) over liquid/ice comes from
1829 !.. Flatau et al. 1992; enthalpy (latent heat) of vaporization from
1830 !.. Bohren & Albrecht 1998; others from Pruppacher & Klett 1978.
1831 !+---+-----------------------------------------------------------------+
1833 tempc = temp(k) - 273.15
1834 rhof(k) = SQRT(RHO_NOT/rho(k))
1835 rhof2(k) = SQRT(rhof(k))
1836 qvs(k) = rslf(pres(k), temp(k))
1837 delQvs(k) = MAX(0.0, rslf(pres(k), 273.15)-qv(k))
1838 if (tempc .le. 0.0) then
1839 qvsi(k) = rsif(pres(k), temp(k))
1843 satw(k) = qv(k)/qvs(k)
1844 sati(k) = qv(k)/qvsi(k)
1845 ssatw(k) = satw(k) - 1.
1846 ssati(k) = sati(k) - 1.
1847 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
1848 if (abs(ssati(k)).lt. eps) ssati(k) = 0.0
1849 if (no_micro .and. ssati(k).gt. 0.0) no_micro = .false.
1850 diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
1851 if (tempc .ge. 0.0) then
1852 visco(k) = (1.718+0.0049*tempc)*1.0E-5
1854 visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
1856 ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
1857 vsc2(k) = SQRT(rho(k)/visco(k))
1858 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
1859 tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
1862 !+---+-----------------------------------------------------------------+
1863 !..If no existing hydrometeor species and no chance to initiate ice or
1864 !.. condense cloud water, just exit quickly!
1865 !+---+-----------------------------------------------------------------+
1867 if (no_micro) return
1869 !+---+-----------------------------------------------------------------+
1870 !..Calculate y-intercept, slope, and useful moments for snow.
1871 !+---+-----------------------------------------------------------------+
1872 if (.not. iiwarm) then
1874 if (.not. L_qs(k)) CYCLE
1875 tc0 = MIN(-0.1, temp(k)-273.15)
1876 smob(k) = rs(k)*oams
1878 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
1879 !.. then we must compute actual 2nd moment and use as reference.
1880 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
1883 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
1884 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
1885 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
1886 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
1887 + sa(10)*bm_s*bm_s*bm_s
1889 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
1890 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
1891 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
1892 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
1893 + sb(10)*bm_s*bm_s*bm_s
1894 smo2(k) = (smob(k)/a_)**(1./b_)
1897 !..Calculate 0th moment. Represents snow number concentration.
1898 loga_ = sa(1) + sa(2)*tc0 + sa(5)*tc0*tc0 + sa(9)*tc0*tc0*tc0
1900 b_ = sb(1) + sb(2)*tc0 + sb(5)*tc0*tc0 + sb(9)*tc0*tc0*tc0
1901 smo0(k) = a_ * smo2(k)**b_
1903 !..Calculate 1st moment. Useful for depositional growth and melting.
1904 loga_ = sa(1) + sa(2)*tc0 + sa(3) &
1905 + sa(4)*tc0 + sa(5)*tc0*tc0 &
1906 + sa(6) + sa(7)*tc0*tc0 &
1907 + sa(8)*tc0 + sa(9)*tc0*tc0*tc0 &
1910 b_ = sb(1)+ sb(2)*tc0 + sb(3) + sb(4)*tc0 &
1911 + sb(5)*tc0*tc0 + sb(6) &
1912 + sb(7)*tc0*tc0 + sb(8)*tc0 &
1913 + sb(9)*tc0*tc0*tc0 + sb(10)
1914 smo1(k) = a_ * smo2(k)**b_
1916 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
1917 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
1918 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
1919 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
1920 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
1921 + sa(10)*cse(1)*cse(1)*cse(1)
1923 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
1924 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
1925 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
1926 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
1927 smoc(k) = a_ * smo2(k)**b_
1929 !..Calculate bv_s+2 (th) moment. Useful for riming.
1930 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(13) &
1931 + sa(4)*tc0*cse(13) + sa(5)*tc0*tc0 &
1932 + sa(6)*cse(13)*cse(13) + sa(7)*tc0*tc0*cse(13) &
1933 + sa(8)*tc0*cse(13)*cse(13) + sa(9)*tc0*tc0*tc0 &
1934 + sa(10)*cse(13)*cse(13)*cse(13)
1936 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(13) + sb(4)*tc0*cse(13) &
1937 + sb(5)*tc0*tc0 + sb(6)*cse(13)*cse(13) &
1938 + sb(7)*tc0*tc0*cse(13) + sb(8)*tc0*cse(13)*cse(13) &
1939 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(13)*cse(13)*cse(13)
1940 smoe(k) = a_ * smo2(k)**b_
1942 !..Calculate 1+(bv_s+1)/2 (th) moment. Useful for depositional growth.
1943 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(16) &
1944 + sa(4)*tc0*cse(16) + sa(5)*tc0*tc0 &
1945 + sa(6)*cse(16)*cse(16) + sa(7)*tc0*tc0*cse(16) &
1946 + sa(8)*tc0*cse(16)*cse(16) + sa(9)*tc0*tc0*tc0 &
1947 + sa(10)*cse(16)*cse(16)*cse(16)
1949 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(16) + sb(4)*tc0*cse(16) &
1950 + sb(5)*tc0*tc0 + sb(6)*cse(16)*cse(16) &
1951 + sb(7)*tc0*tc0*cse(16) + sb(8)*tc0*cse(16)*cse(16) &
1952 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(16)*cse(16)*cse(16)
1953 smof(k) = a_ * smo2(k)**b_
1957 !+---+-----------------------------------------------------------------+
1958 !..Calculate y-intercept, slope values for graupel.
1959 !+---+-----------------------------------------------------------------+
1962 ygra1 = alog10(max(1.E-9, rg(k)))
1963 zans1 = 3.0 + 2./7.*(ygra1+8.)
1964 zans1 = MAX(2., MIN(zans1, 6.))
1965 N0_exp = 10.**(zans1)
1966 lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
1967 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
1969 N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2)
1974 !+---+-----------------------------------------------------------------+
1975 !..Calculate y-intercept, slope values for rain.
1976 !+---+-----------------------------------------------------------------+
1978 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
1980 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
1981 N0_r(k) = nr(k)*org2*lamr**cre(2)
1984 !+---+-----------------------------------------------------------------+
1985 !..Compute warm-rain process terms (except evap done later).
1986 !+---+-----------------------------------------------------------------+
1990 !..Rain self-collection follows Seifert, 1994 and drop break-up
1991 !.. follows Verlinde and Cotton, 1993. RAIN2M
1992 if (L_qr(k) .and. mvd_r(k).gt. D0r) then
1994 !-GT if (mvd_r(k) .gt. 1500.0E-6) then
1995 Ef_rr = 1.0 - EXP(2300.0*(mvd_r(k)-1950.0E-6))
1997 pnr_rcr(k) = Ef_rr * 2.0*nr(k)*rr(k)
2002 nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2)
2003 xDc = MAX(D0c*1.E6, ((rc(k)/(am_r*nc(k)))**obmr) * 1.E6)
2004 lamc = (nc(k)*am_r* ccg(2,nu_c) * ocg1(nu_c) / rc(k))**obmr
2005 mvd_c(k) = (3.0+nu_c+0.672) / lamc
2008 !..Autoconversion follows Berry & Reinhardt (1974) with characteristic
2009 !.. diameters correctly computed from gamma distrib of cloud droplets.
2010 if (rc(k).gt. 0.01e-3) then
2011 Dc_g = ((ccg(3,nu_c)*ocg2(nu_c))**obmr / lamc) * 1.E6
2012 Dc_b = (xDc*xDc*xDc*Dc_g*Dc_g*Dc_g - xDc*xDc*xDc*xDc*xDc*xDc) &
2014 zeta1 = 0.5*((6.25E-6*xDc*Dc_b*Dc_b*Dc_b - 0.4) &
2015 + abs(6.25E-6*xDc*Dc_b*Dc_b*Dc_b - 0.4))
2016 zeta = 0.027*rc(k)*zeta1
2017 taud = 0.5*((0.5*Dc_b - 7.5) + abs(0.5*Dc_b - 7.5)) + R1
2018 tau = 3.72/(rc(k)*taud)
2019 prr_wau(k) = zeta/tau
2020 prr_wau(k) = MIN(DBLE(rc(k)*odts), prr_wau(k))
2021 pnr_wau(k) = prr_wau(k) / (am_r*nu_c*200.*D0r*D0r*D0r) ! RAIN2M
2022 pnc_wau(k) = MIN(DBLE(nc(k)*odts), prr_wau(k) &
2023 / (am_r*mvd_c(k)*mvd_c(k)*mvd_c(k))) ! Qc2M
2026 !..Rain collecting cloud water. In CE, assume Dc<<Dr and vtc=~0.
2027 if (L_qr(k) .and. mvd_r(k).gt. D0r .and. mvd_c(k).gt. D0c) then
2029 idx = 1 + INT(nbr*DLOG(mvd_r(k)/Dr(1))/DLOG(Dr(nbr)/Dr(1)))
2031 Ef_rw = t_Efrw(idx, INT(mvd_c(k)*1.E6))
2032 prr_rcw(k) = rhof(k)*t1_qr_qc*Ef_rw*rc(k)*N0_r(k) &
2033 *((lamr+fv_r)**(-cre(9)))
2034 prr_rcw(k) = MIN(DBLE(rc(k)*odts), prr_rcw(k))
2035 pnc_rcw(k) = rhof(k)*t1_qr_qc*Ef_rw*nc(k)*N0_r(k) &
2036 *((lamr+fv_r)**(-cre(9))) ! Qc2M
2037 pnc_rcw(k) = MIN(DBLE(nc(k)*odts), pnc_rcw(k))
2040 !..Rain collecting aerosols, wet scavenging.
2041 if (L_qr(k) .and. mvd_r(k).gt. D0r) then
2042 Ef_ra = Eff_aero(mvd_r(k),0.04E-6,visco(k),rho(k),temp(k),'r')
2044 pna_rca(k) = rhof(k)*t1_qr_qc*Ef_ra*nwfa(k)*N0_r(k) &
2045 *((lamr+fv_r)**(-cre(9)))
2046 pna_rca(k) = MIN(DBLE(nwfa(k)*odts), pna_rca(k))
2048 Ef_ra = Eff_aero(mvd_r(k),0.8E-6,visco(k),rho(k),temp(k),'r')
2049 pnd_rcd(k) = rhof(k)*t1_qr_qc*Ef_ra*nifa(k)*N0_r(k) &
2050 *((lamr+fv_r)**(-cre(9)))
2051 pnd_rcd(k) = MIN(DBLE(nifa(k)*odts), pnd_rcd(k))
2053 if (present(wif_input_opt)) then
2054 if (wif_input_opt .eq. 2) then
2055 Ef_ra = Eff_aero(mvd_r(k),0.0236E-6, &
2056 visco(k),rho(k),temp(k),'r')
2057 pnb_rcb(k) = rhof(k)*t1_qr_qc*Ef_ra*nbca(k)*N0_r(k) &
2058 *((lamr+fv_r)**(-cre(9)))
2059 pnb_rcb(k) = MIN(DBLE(nbca(k)*odts), pnb_rcb(k))
2066 !+---+-----------------------------------------------------------------+
2067 !..Compute all frozen hydrometeor species' process terms.
2068 !+---+-----------------------------------------------------------------+
2069 if (.not. iiwarm) then
2070 !..vts_boost is the factor applied to snow terminal
2071 !..fallspeed due to riming of snow
2075 if (L_qs(k)) xDs = smoc(k) / smob(k)
2077 !..Temperature lookup table indexes.
2078 tempc = temp(k) - 273.15
2079 idx_tc = MAX(1, MIN(NINT(-tempc), 45) )
2080 idx_t = INT( (tempc-2.5)/5. ) - 1
2081 idx_t = MAX(1, -idx_t)
2082 idx_t = MIN(idx_t, ntb_t)
2083 IT = MAX(1, MIN(NINT(-tempc), 31) )
2085 !..Cloud water lookup table index.
2086 if (rc(k).gt. r_c(1)) then
2087 nic = NINT(ALOG10(rc(k)))
2088 do nn = nic-1, nic+1
2090 if ( (rc(k)/10.**nn).ge.1.0 .and. &
2091 (rc(k)/10.**nn).lt.10.0) goto 141
2094 idx_c = INT(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
2095 idx_c = MAX(1, MIN(idx_c, ntb_c))
2100 !..Cloud droplet number lookup table index.
2101 idx_n = NINT(1.0 + FLOAT(nbc) * DLOG(nc(k)/t_Nc(1)) / nic1)
2102 idx_n = MAX(1, MIN(idx_n, nbc))
2104 !..Cloud ice lookup table indexes.
2105 if (ri(k).gt. r_i(1)) then
2106 nii = NINT(ALOG10(ri(k)))
2107 do nn = nii-1, nii+1
2109 if ( (ri(k)/10.**nn).ge.1.0 .and. &
2110 (ri(k)/10.**nn).lt.10.0) goto 142
2113 idx_i = INT(ri(k)/10.**n) + 10*(n-nii2) - (n-nii2)
2114 idx_i = MAX(1, MIN(idx_i, ntb_i))
2119 if (ni(k).gt. Nt_i(1)) then
2120 nii = NINT(ALOG10(ni(k)))
2121 do nn = nii-1, nii+1
2123 if ( (ni(k)/10.**nn).ge.1.0 .and. &
2124 (ni(k)/10.**nn).lt.10.0) goto 143
2127 idx_i1 = INT(ni(k)/10.**n) + 10*(n-nii3) - (n-nii3)
2128 idx_i1 = MAX(1, MIN(idx_i1, ntb_i1))
2133 !..Rain lookup table indexes.
2134 if (rr(k).gt. r_r(1)) then
2135 nir = NINT(ALOG10(rr(k)))
2136 do nn = nir-1, nir+1
2138 if ( (rr(k)/10.**nn).ge.1.0 .and. &
2139 (rr(k)/10.**nn).lt.10.0) goto 144
2142 idx_r = INT(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2)
2143 idx_r = MAX(1, MIN(idx_r, ntb_r))
2146 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
2147 N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
2148 nir = NINT(DLOG10(N0_exp))
2149 do nn = nir-1, nir+1
2151 if ( (N0_exp/10.**nn).ge.1.0 .and. &
2152 (N0_exp/10.**nn).lt.10.0) goto 145
2155 idx_r1 = INT(N0_exp/10.**n) + 10*(n-nir3) - (n-nir3)
2156 idx_r1 = MAX(1, MIN(idx_r1, ntb_r1))
2162 !..Snow lookup table index.
2163 if (rs(k).gt. r_s(1)) then
2164 nis = NINT(ALOG10(rs(k)))
2165 do nn = nis-1, nis+1
2167 if ( (rs(k)/10.**nn).ge.1.0 .and. &
2168 (rs(k)/10.**nn).lt.10.0) goto 146
2171 idx_s = INT(rs(k)/10.**n) + 10*(n-nis2) - (n-nis2)
2172 idx_s = MAX(1, MIN(idx_s, ntb_s))
2177 !..Graupel lookup table index.
2178 if (rg(k).gt. r_g(1)) then
2179 nig = NINT(ALOG10(rg(k)))
2180 do nn = nig-1, nig+1
2182 if ( (rg(k)/10.**nn).ge.1.0 .and. &
2183 (rg(k)/10.**nn).lt.10.0) goto 147
2186 idx_g = INT(rg(k)/10.**n) + 10*(n-nig2) - (n-nig2)
2187 idx_g = MAX(1, MIN(idx_g, ntb_g))
2190 lam_exp = lamg * (cgg(3)*ogg2*ogg1)**bm_g
2191 N0_exp = ogg1*rg(k)/am_g * lam_exp**cge(1)
2192 nig = NINT(DLOG10(N0_exp))
2193 do nn = nig-1, nig+1
2195 if ( (N0_exp/10.**nn).ge.1.0 .and. &
2196 (N0_exp/10.**nn).lt.10.0) goto 148
2199 idx_g1 = INT(N0_exp/10.**n) + 10*(n-nig3) - (n-nig3)
2200 idx_g1 = MAX(1, MIN(idx_g1, ntb_g1))
2206 !..Deposition/sublimation prefactor (from Srivastava & Coen 1992).
2208 rvs = rho(k)*qvsi(k)
2209 rvs_p = rvs*otemp*(lsub*otemp*oRv - 1.)
2210 rvs_pp = rvs * ( otemp*(lsub*otemp*oRv - 1.) &
2211 *otemp*(lsub*otemp*oRv - 1.) &
2212 + (-2.*lsub*otemp*otemp*otemp*oRv) &
2214 gamsc = lsub*diffu(k)/tcond(k) * rvs_p
2215 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
2216 * rvs_pp/rvs_p * rvs/rvs_p
2217 alphsc = MAX(1.E-9, alphsc)
2219 if (abs(xsat).lt. 1.E-9) xsat=0.
2220 t1_subl = 4.*PI*( 1.0 - alphsc*xsat &
2221 + 2.*alphsc*alphsc*xsat*xsat &
2222 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
2225 !..Snow collecting cloud water. In CE, assume Dc<<Ds and vtc=~0.
2226 if (L_qc(k) .and. mvd_c(k).gt. D0c) then
2227 if (xDs .gt. D0s) then
2228 idx = 1 + INT(nbs*DLOG(xDs/Ds(1))/DLOG(Ds(nbs)/Ds(1)))
2230 Ef_sw = t_Efsw(idx, INT(mvd_c(k)*1.E6))
2231 prs_scw(k) = rhof(k)*t1_qs_qc*Ef_sw*rc(k)*smoe(k)
2232 prs_scw(k) = MIN(DBLE(rc(k)*odts), prs_scw(k))
2233 pnc_scw(k) = rhof(k)*t1_qs_qc*Ef_sw*nc(k)*smoe(k) ! Qc2M
2234 pnc_scw(k) = MIN(DBLE(nc(k)*odts), pnc_scw(k))
2237 !..Graupel collecting cloud water. In CE, assume Dc<<Dg and vtc=~0.
2238 if (rg(k).ge. r_g(1) .and. mvd_c(k).gt. D0c) then
2239 xDg = (bm_g + mu_g + 1.) * ilamg(k)
2240 vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
2241 stoke_g = mvd_c(k)*mvd_c(k)*vtg*rho_w/(9.*visco(k)*xDg)
2242 if (xDg.gt. D0g) then
2243 if (stoke_g.ge.0.4 .and. stoke_g.le.10.) then
2244 Ef_gw = 0.55*ALOG10(2.51*stoke_g)
2245 elseif (stoke_g.lt.0.4) then
2247 elseif (stoke_g.gt.10) then
2250 prg_gcw(k) = rhof(k)*t1_qg_qc*Ef_gw*rc(k)*N0_g(k) &
2252 pnc_gcw(k) = rhof(k)*t1_qg_qc*Ef_gw*nc(k)*N0_g(k) &
2253 *ilamg(k)**cge(9) ! Qc2M
2254 pnc_gcw(k) = MIN(DBLE(nc(k)*odts), pnc_gcw(k))
2259 !..Snow and graupel collecting aerosols, wet scavenging.
2260 if (rs(k) .gt. r_s(1)) then
2261 Ef_sa = Eff_aero(xDs,0.04E-6,visco(k),rho(k),temp(k),'s')
2262 pna_sca(k) = rhof(k)*t1_qs_qc*Ef_sa*nwfa(k)*smoe(k)
2263 pna_sca(k) = MIN(DBLE(nwfa(k)*odts), pna_sca(k))
2265 Ef_sa = Eff_aero(xDs,0.8E-6,visco(k),rho(k),temp(k),'s')
2266 pnd_scd(k) = rhof(k)*t1_qs_qc*Ef_sa*nifa(k)*smoe(k)
2267 pnd_scd(k) = MIN(DBLE(nifa(k)*odts), pnd_scd(k))
2269 if (present(wif_input_opt)) then
2270 if (wif_input_opt .eq. 2) then
2271 Ef_sa = Eff_aero(xDs,0.0236E-6,visco(k),rho(k),temp(k),'s')
2272 pnb_scb(k) = rhof(k)*t1_qs_qc*Ef_sa*nbca(k)*smoe(k)
2273 pnb_scb(k) = MIN(DBLE(nbca(k)*odts), pnb_scb(k))
2277 if (rg(k) .gt. r_g(1)) then
2278 xDg = (bm_g + mu_g + 1.) * ilamg(k)
2279 Ef_ga = Eff_aero(xDg,0.04E-6,visco(k),rho(k),temp(k),'g')
2280 pna_gca(k) = rhof(k)*t1_qg_qc*Ef_ga*nwfa(k)*N0_g(k) &
2282 pna_gca(k) = MIN(DBLE(nwfa(k)*odts), pna_gca(k))
2284 Ef_ga = Eff_aero(xDg,0.8E-6,visco(k),rho(k),temp(k),'g')
2285 pnd_gcd(k) = rhof(k)*t1_qg_qc*Ef_ga*nifa(k)*N0_g(k) &
2287 pnd_gcd(k) = MIN(DBLE(nifa(k)*odts), pnd_gcd(k))
2289 if (present(wif_input_opt)) then
2290 if (wif_input_opt .eq. 2) then
2291 Ef_ga = Eff_aero(xDg,0.0236E-6,visco(k),rho(k),temp(k),'g')
2292 pnb_gcb(k) = rhof(k)*t1_qg_qc*Ef_ga*nbca(k)*N0_g(k) &
2294 pnb_gcb(k) = MIN(DBLE(nbca(k)*odts), pnb_gcb(k))
2299 !..Rain collecting snow. Cannot assume Wisner (1972) approximation
2300 !.. or Mizuno (1990) approach so we solve the CE explicitly and store
2301 !.. results in lookup table.
2302 if (rr(k).ge. r_r(1)) then
2303 if (rs(k).ge. r_s(1)) then
2304 if (temp(k).lt.T_0) then
2305 prr_rcs(k) = -(tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2306 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
2307 + tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
2308 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r))
2309 prs_rcs(k) = tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2310 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
2311 - tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2312 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
2313 prg_rcs(k) = tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
2314 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2315 + tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2316 + tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
2317 prr_rcs(k) = MAX(DBLE(-rr(k)*odts), prr_rcs(k))
2318 prs_rcs(k) = MAX(DBLE(-rs(k)*odts), prs_rcs(k))
2319 prg_rcs(k) = MIN(DBLE((rr(k)+rs(k))*odts), prg_rcs(k))
2320 pnr_rcs(k) = tnr_racs1(idx_s,idx_t,idx_r1,idx_r) & ! RAIN2M
2321 + tnr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2322 + tnr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2323 + tnr_sacr2(idx_s,idx_t,idx_r1,idx_r)
2324 pnr_rcs(k) = MIN(DBLE(nr(k)*odts), pnr_rcs(k))
2326 prs_rcs(k) = -tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
2327 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r) &
2328 + tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
2329 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r)
2330 prs_rcs(k) = MAX(DBLE(-rs(k)*odts), prs_rcs(k))
2331 prr_rcs(k) = -prs_rcs(k)
2335 !..Rain collecting graupel. Cannot assume Wisner (1972) approximation
2336 !.. or Mizuno (1990) approach so we solve the CE explicitly and store
2337 !.. results in lookup table.
2338 if (rg(k).ge. r_g(1)) then
2339 if (temp(k).lt.T_0) then
2340 prg_rcg(k) = tmr_racg(idx_g1,idx_g,idx_r1,idx_r) &
2341 + tcr_gacr(idx_g1,idx_g,idx_r1,idx_r)
2342 prg_rcg(k) = MIN(DBLE(rr(k)*odts), prg_rcg(k))
2343 prr_rcg(k) = -prg_rcg(k)
2344 pnr_rcg(k) = tnr_racg(idx_g1,idx_g,idx_r1,idx_r) & ! RAIN2M
2345 + tnr_gacr(idx_g1,idx_g,idx_r1,idx_r)
2346 pnr_rcg(k) = MIN(DBLE(nr(k)*odts), pnr_rcg(k))
2348 prr_rcg(k) = tcg_racg(idx_g1,idx_g,idx_r1,idx_r)
2349 prr_rcg(k) = MIN(DBLE(rg(k)*odts), prr_rcg(k))
2350 prg_rcg(k) = -prr_rcg(k)
2351 !..Put in explicit drop break-up due to collisions.
2352 pnr_rcg(k) = -5.*tnr_gacr(idx_g1,idx_g,idx_r1,idx_r) ! RAIN2M
2357 !+---+-----------------------------------------------------------------+
2358 !..Next IF block handles only those processes below 0C.
2359 !+---+-----------------------------------------------------------------+
2361 if (temp(k).lt.T_0) then
2364 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
2366 !+---+---------------- BEGIN NEW ICE NUCLEATION -----------------------+
2367 !..Freezing of supercooled water (rain or cloud) is influenced by dust
2368 !.. but still using Bigg 1953 with a temperature adjustment of a few
2369 !.. degrees depending on dust concentration. A default value by way
2370 !.. of idx_IN is 1.0 per Liter of air is used when dustyIce flag is
2371 !.. false. Next, a combination of deposition/condensation freezing
2372 !.. using DeMott et al (2010) dust nucleation when water saturated or
2373 !.. Phillips et al (2008) when below water saturation; else, without
2374 !.. dustyIce flag, use the previous Cooper (1986) temperature-dependent
2375 !.. value. Lastly, allow homogeneous freezing of deliquesced aerosols
2376 !.. following Koop et al. (2001, Nature).
2377 !.. Implemented by T. Eidhammer and G. Thompson 2012Dec18
2378 !+---+-----------------------------------------------------------------+
2380 if (dustyIce .AND. is_aerosol_aware) then
2381 xni = iceDeMott(tempc,qvs(k),qvs(k),qvsi(k),rho(k),nifa(k))
2383 xni = 1.0 *1000. ! Default is 1.0 per Liter
2386 !..Ice nuclei lookup table index.
2387 if (xni.gt. Nt_IN(1)) then
2388 niin = NINT(ALOG10(xni))
2389 do nn = niin-1, niin+1
2391 if ( (xni/10.**nn).ge.1.0 .and. &
2392 (xni/10.**nn).lt.10.0) goto 149
2395 idx_IN = INT(xni/10.**n) + 10*(n-niin2) - (n-niin2)
2396 idx_IN = MAX(1, MIN(idx_IN, ntb_IN))
2401 !..Freezing of water drops into graupel/cloud ice (Bigg 1953).
2402 if (rr(k).gt. r_r(1)) then
2403 prg_rfz(k) = tpg_qrfz(idx_r,idx_r1,idx_tc,idx_IN)*odts
2404 pri_rfz(k) = tpi_qrfz(idx_r,idx_r1,idx_tc,idx_IN)*odts
2405 pni_rfz(k) = tni_qrfz(idx_r,idx_r1,idx_tc,idx_IN)*odts
2406 pnr_rfz(k) = tnr_qrfz(idx_r,idx_r1,idx_tc,idx_IN)*odts ! RAIN2M
2407 pnr_rfz(k) = MIN(DBLE(nr(k)*odts), pnr_rfz(k))
2408 elseif (rr(k).gt. R1 .and. temp(k).lt.HGFR) then
2409 pri_rfz(k) = rr(k)*odts
2410 pni_rfz(k) = nr(k)*odts ! RAIN2M
2413 if (rc(k).gt. r_c(1)) then
2414 pri_wfz(k) = tpi_qcfz(idx_c,idx_n,idx_tc,idx_IN)*odts
2415 pri_wfz(k) = MIN(DBLE(rc(k)*odts), pri_wfz(k))
2416 pni_wfz(k) = tni_qcfz(idx_c,idx_n,idx_tc,idx_IN)*odts
2417 pni_wfz(k) = MIN(DBLE(nc(k)*odts), pri_wfz(k)/(2.*xm0i), &
2419 elseif (rc(k).gt. R1 .and. temp(k).lt.HGFR) then
2420 pri_wfz(k) = rc(k)*odts
2421 pni_wfz(k) = nc(k)*odts
2424 !..Deposition nucleation of dust/mineral from DeMott et al (2010)
2425 !.. we may need to relax the temperature and ssati constraints.
2426 if ( (ssati(k).ge. 0.25) .or. (ssatw(k).gt. eps &
2427 .and. temp(k).lt.253.15) ) then
2428 if (dustyIce .AND. is_aerosol_aware) then
2429 xnc = iceDeMott(tempc,qv(k),qvs(k),qvsi(k),rho(k),nifa(k))
2431 xnc = MIN(250.E3, TNO*EXP(ATO*(T_0-temp(k))))
2433 xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dtsave
2434 pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts
2435 pri_inu(k) = MIN(DBLE(rate_max), xm0i*pni_inu(k))
2436 pni_inu(k) = pri_inu(k)/xm0i
2439 !..Freezing of aqueous aerosols based on Koop et al (2001, Nature)
2440 xni = smo0(k)+ni(k) + (pni_rfz(k)+pni_wfz(k)+pni_inu(k))*dtsave
2441 if (is_aerosol_aware .AND. homogIce .AND. (xni.le.999.E3) &
2442 & .AND.(temp(k).lt.238).AND.(ssati(k).ge.0.4) ) then
2443 xnc = iceKoop(temp(k),qv(k),qvs(k),nwfa(k), dtsave)
2444 pni_iha(k) = xnc*odts
2445 pri_iha(k) = MIN(DBLE(rate_max), xm0i*0.1*pni_iha(k))
2446 pni_iha(k) = pri_iha(k)/(xm0i*0.1)
2448 !+---+------------------ END NEW ICE NUCLEATION -----------------------+
2451 !..Deposition/sublimation of cloud ice (Srivastava & Coen 1992).
2453 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
2455 xDi = MAX(DBLE(D0i), (bm_i + mu_i + 1.) * ilami)
2456 xmi = am_i*xDi**bm_i
2458 pri_ide(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
2459 *oig1*cig(5)*ni(k)*ilami
2461 if (pri_ide(k) .lt. 0.0) then
2462 pri_ide(k) = MAX(DBLE(-ri(k)*odts), pri_ide(k), DBLE(rate_max))
2463 pni_ide(k) = pri_ide(k)*oxmi
2464 pni_ide(k) = MAX(DBLE(-ni(k)*odts), pni_ide(k))
2466 pri_ide(k) = MIN(pri_ide(k), DBLE(rate_max))
2467 prs_ide(k) = (1.0D0-tpi_ide(idx_i,idx_i1))*pri_ide(k)
2468 pri_ide(k) = tpi_ide(idx_i,idx_i1)*pri_ide(k)
2471 !..Some cloud ice needs to move into the snow category. Use lookup
2472 !.. table that resulted from explicit bin representation of distrib.
2473 if ( (idx_i.eq. ntb_i) .or. (xDi.gt. 5.0*D0s) ) then
2474 prs_iau(k) = ri(k)*.99*odts
2475 pni_iau(k) = ni(k)*.95*odts
2476 elseif (xDi.lt. 0.1*D0s) then
2480 prs_iau(k) = tps_iaus(idx_i,idx_i1)*odts
2481 prs_iau(k) = MIN(DBLE(ri(k)*.99*odts), prs_iau(k))
2482 pni_iau(k) = tni_iaus(idx_i,idx_i1)*odts
2483 pni_iau(k) = MIN(DBLE(ni(k)*.95*odts), pni_iau(k))
2487 !..Deposition/sublimation of snow/graupel follows Srivastava & Coen
2490 C_snow = C_sqrd + (tempc+1.5)*(C_cube-C_sqrd)/(-30.+1.5)
2491 C_snow = MAX(C_sqrd, MIN(C_snow, C_cube))
2492 prs_sde(k) = C_snow*t1_subl*diffu(k)*ssati(k)*rvs &
2493 * (t1_qs_sd*smo1(k) &
2494 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
2495 if (prs_sde(k).lt. 0.) then
2496 prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k), DBLE(rate_max))
2498 prs_sde(k) = MIN(prs_sde(k), DBLE(rate_max))
2502 if (L_qg(k) .and. ssati(k).lt. -eps) then
2503 prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
2504 * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
2505 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
2506 if (prg_gde(k).lt. 0.) then
2507 prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k), DBLE(rate_max))
2509 prg_gde(k) = MIN(prg_gde(k), DBLE(rate_max))
2513 !..Snow collecting cloud ice. In CE, assume Di<<Ds and vti=~0.
2515 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
2517 xDi = MAX(DBLE(D0i), (bm_i + mu_i + 1.) * ilami)
2518 xmi = am_i*xDi**bm_i
2520 if (rs(k).ge. r_s(1)) then
2521 prs_sci(k) = t1_qs_qi*rhof(k)*Ef_si*ri(k)*smoe(k)
2522 pni_sci(k) = prs_sci(k) * oxmi
2525 !..Rain collecting cloud ice. In CE, assume Di<<Dr and vti=~0.
2526 if (rr(k).ge. r_r(1) .and. mvd_r(k).gt. 4.*xDi) then
2528 pri_rci(k) = rhof(k)*t1_qr_qi*Ef_ri*ri(k)*N0_r(k) &
2529 *((lamr+fv_r)**(-cre(9)))
2530 pnr_rci(k) = rhof(k)*t1_qr_qi*Ef_ri*ni(k)*N0_r(k) & ! RAIN2M
2531 *((lamr+fv_r)**(-cre(9)))
2532 pni_rci(k) = pri_rci(k) * oxmi
2533 prr_rci(k) = rhof(k)*t2_qr_qi*Ef_ri*ni(k)*N0_r(k) &
2534 *((lamr+fv_r)**(-cre(8)))
2535 prr_rci(k) = MIN(DBLE(rr(k)*odts), prr_rci(k))
2536 prg_rci(k) = pri_rci(k) + prr_rci(k)
2540 !..Ice multiplication from rime-splinters (Hallet & Mossop 1974).
2541 if (prg_gcw(k).gt. eps .and. tempc.gt.-8.0) then
2543 if (tempc.ge.-5.0 .and. tempc.lt.-3.0) then
2544 tf = 0.5*(-3.0 - tempc)
2545 elseif (tempc.gt.-8.0 .and. tempc.lt.-5.0) then
2546 tf = 0.33333333*(8.0 + tempc)
2548 pni_ihm(k) = 3.5E8*tf*prg_gcw(k)
2549 pri_ihm(k) = xm0i*pni_ihm(k)
2550 prs_ihm(k) = prs_scw(k)/(prs_scw(k)+prg_gcw(k)) &
2552 prg_ihm(k) = prg_gcw(k)/(prs_scw(k)+prg_gcw(k)) &
2556 !..A portion of rimed snow converts to graupel but some remains snow.
2557 !.. Interp from 15 to 95% as riming factor increases from 2.0 to 30.0
2558 !.. 0.028 came from (.95-.15)/(30.-2.). This remains ad-hoc and should
2560 if (prs_scw(k).gt.2.0*prs_sde(k) .and. &
2561 prs_sde(k).gt.eps) then
2562 r_frac = MIN(30.0D0, prs_scw(k)/prs_sde(k))
2563 g_frac = MIN(0.95, 0.15 + (r_frac-2.)*.028)
2564 vts_boost(k) = MIN(1.5, 1.1 + (r_frac-2.)*.014)
2565 prg_scw(k) = g_frac*prs_scw(k)
2566 prs_scw(k) = (1. - g_frac)*prs_scw(k)
2571 !..Melt snow and graupel and enhance from collisions with liquid.
2572 !.. We also need to sublimate snow and graupel if subsaturated.
2574 prr_sml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delQvs(k)) &
2575 * (t1_qs_me*smo1(k) + t2_qs_me*rhof2(k)*vsc2(k)*smof(k))
2576 if (prr_sml(k) .gt. 0.) then
2577 prr_sml(k) = prr_sml(k) + 4218.*olfus*tempc &
2578 * (prr_rcs(k)+prs_scw(k))
2580 prr_sml(k) = MIN(DBLE(rs(k)*odts), MAX(0.D0, prr_sml(k)))
2581 pnr_sml(k) = smo0(k)/rs(k)*prr_sml(k) * 10.0**(-0.25*tempc) ! RAIN2M
2582 pnr_sml(k) = MIN(DBLE(smo0(k)*odts), pnr_sml(k))
2584 if (ssati(k).lt. 0.) then
2585 prs_sde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
2586 * (t1_qs_sd*smo1(k) &
2587 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
2588 prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k))
2593 prr_gml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delQvs(k)) &
2594 * N0_g(k)*(t1_qg_me*ilamg(k)**cge(10) &
2595 + t2_qg_me*rhof2(k)*vsc2(k)*ilamg(k)**cge(11))
2596 !-GT prr_gml(k) = prr_gml(k) + 4218.*olfus*tempc &
2597 !-GT * (prr_rcg(k)+prg_gcw(k))
2598 prr_gml(k) = MIN(DBLE(rg(k)*odts), MAX(0.D0, prr_gml(k)))
2599 pnr_gml(k) = N0_g(k)*cgg(2)*ilamg(k)**cge(2) / rg(k) & ! RAIN2M
2600 * prr_gml(k) * 10.0**(-0.5*tempc)
2602 if (ssati(k).lt. 0.) then
2603 prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
2604 * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
2605 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
2606 prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k))
2610 !.. This change will be required if users run adaptive time step that
2611 !.. results in delta-t that is generally too long to allow cloud water
2612 !.. collection by snow/graupel above melting temperature.
2613 !.. Credit to Bjorn-Egil Nygaard for discovering.
2614 if (dt .gt. 120.) then
2615 prr_rcw(k)=prr_rcw(k)+prs_scw(k)+prg_gcw(k)
2625 !+---+-----------------------------------------------------------------+
2626 !..Ensure we do not deplete more hydrometeor species than exists.
2627 !+---+-----------------------------------------------------------------+
2630 !..If ice supersaturated, ensure sum of depos growth terms does not
2631 !.. deplete more vapor than possibly exists. If subsaturated, limit
2632 !.. sum of sublimation terms such that vapor does not reproduce ice
2634 sump = pri_inu(k) + pri_ide(k) + prs_ide(k) &
2635 + prs_sde(k) + prg_gde(k) + pri_iha(k)
2636 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
2637 if ( (sump.gt. eps .and. sump.gt. rate_max) .or. &
2638 (sump.lt. -eps .and. sump.lt. rate_max) ) then
2639 ratio = rate_max/sump
2640 pri_inu(k) = pri_inu(k) * ratio
2641 pri_ide(k) = pri_ide(k) * ratio
2642 pni_ide(k) = pni_ide(k) * ratio
2643 prs_ide(k) = prs_ide(k) * ratio
2644 prs_sde(k) = prs_sde(k) * ratio
2645 prg_gde(k) = prg_gde(k) * ratio
2646 pri_iha(k) = pri_iha(k) * ratio
2649 !..Cloud water conservation.
2650 sump = -prr_wau(k) - pri_wfz(k) - prr_rcw(k) &
2651 - prs_scw(k) - prg_scw(k) - prg_gcw(k)
2652 rate_max = -rc(k)*odts
2653 if (sump.lt. rate_max .and. L_qc(k)) then
2654 ratio = rate_max/sump
2655 prr_wau(k) = prr_wau(k) * ratio
2656 pri_wfz(k) = pri_wfz(k) * ratio
2657 prr_rcw(k) = prr_rcw(k) * ratio
2658 prs_scw(k) = prs_scw(k) * ratio
2659 prg_scw(k) = prg_scw(k) * ratio
2660 prg_gcw(k) = prg_gcw(k) * ratio
2663 !..Cloud ice conservation.
2664 sump = pri_ide(k) - prs_iau(k) - prs_sci(k) &
2666 rate_max = -ri(k)*odts
2667 if (sump.lt. rate_max .and. L_qi(k)) then
2668 ratio = rate_max/sump
2669 pri_ide(k) = pri_ide(k) * ratio
2670 prs_iau(k) = prs_iau(k) * ratio
2671 prs_sci(k) = prs_sci(k) * ratio
2672 pri_rci(k) = pri_rci(k) * ratio
2675 !..Rain conservation.
2676 sump = -prg_rfz(k) - pri_rfz(k) - prr_rci(k) &
2677 + prr_rcs(k) + prr_rcg(k)
2678 rate_max = -rr(k)*odts
2679 if (sump.lt. rate_max .and. L_qr(k)) then
2680 ratio = rate_max/sump
2681 prg_rfz(k) = prg_rfz(k) * ratio
2682 pri_rfz(k) = pri_rfz(k) * ratio
2683 prr_rci(k) = prr_rci(k) * ratio
2684 prr_rcs(k) = prr_rcs(k) * ratio
2685 prr_rcg(k) = prr_rcg(k) * ratio
2688 !..Snow conservation.
2689 sump = prs_sde(k) - prs_ihm(k) - prr_sml(k) &
2691 rate_max = -rs(k)*odts
2692 if (sump.lt. rate_max .and. L_qs(k)) then
2693 ratio = rate_max/sump
2694 prs_sde(k) = prs_sde(k) * ratio
2695 prs_ihm(k) = prs_ihm(k) * ratio
2696 prr_sml(k) = prr_sml(k) * ratio
2697 prs_rcs(k) = prs_rcs(k) * ratio
2700 !..Graupel conservation.
2701 sump = prg_gde(k) - prg_ihm(k) - prr_gml(k) &
2703 rate_max = -rg(k)*odts
2704 if (sump.lt. rate_max .and. L_qg(k)) then
2705 ratio = rate_max/sump
2706 prg_gde(k) = prg_gde(k) * ratio
2707 prg_ihm(k) = prg_ihm(k) * ratio
2708 prr_gml(k) = prr_gml(k) * ratio
2709 prg_rcg(k) = prg_rcg(k) * ratio
2712 !..Re-enforce proper mass conservation for subsequent elements in case
2713 !.. any of the above terms were altered. Thanks P. Blossey. 2009Sep28
2714 pri_ihm(k) = prs_ihm(k) + prg_ihm(k)
2715 ratio = MIN( ABS(prr_rcg(k)), ABS(prg_rcg(k)) )
2716 prr_rcg(k) = ratio * SIGN(1.0, SNGL(prr_rcg(k)))
2717 prg_rcg(k) = -prr_rcg(k)
2718 if (temp(k).gt.T_0) then
2719 ratio = MIN( ABS(prr_rcs(k)), ABS(prs_rcs(k)) )
2720 prr_rcs(k) = ratio * SIGN(1.0, SNGL(prr_rcs(k)))
2721 prs_rcs(k) = -prr_rcs(k)
2726 !+---+-----------------------------------------------------------------+
2727 !..Calculate tendencies of all species but constrain the number of ice
2728 !.. to reasonable values.
2729 !+---+-----------------------------------------------------------------+
2732 lfus2 = lsub - lvap(k)
2734 !..Aerosol number tendency
2735 if (is_aerosol_aware) then
2736 nwfaten(k) = nwfaten(k) - (pna_rca(k) + pna_sca(k) &
2737 + pna_gca(k) + pni_iha(k)) * orho
2738 nifaten(k) = nifaten(k) - (pnd_rcd(k) + pnd_scd(k) &
2739 + pnd_gcd(k)) * orho
2740 if (wif_input_opt .eq. 2) then
2741 nbcaten(k) = nbcaten(k) - (pnb_rcb(k) + pnb_scb(k) &
2742 + pnb_gcb(k)) * orho
2747 nifaten(k) = nifaten(k) - pni_inu(k)*orho
2753 !..Water vapor tendency
2754 qvten(k) = qvten(k) + (-pri_inu(k) - pri_iha(k) - pri_ide(k) &
2755 - prs_ide(k) - prs_sde(k) - prg_gde(k)) &
2758 !..Cloud water tendency
2759 qcten(k) = qcten(k) + (-prr_wau(k) - pri_wfz(k) &
2760 - prr_rcw(k) - prs_scw(k) - prg_scw(k) &
2764 !..Cloud water number tendency
2765 ncten(k) = ncten(k) + (-pnc_wau(k) - pnc_rcw(k) &
2766 - pni_wfz(k) - pnc_scw(k) - pnc_gcw(k)) &
2769 !..Cloud water mass/number balance; keep mass-wt mean size between
2770 !.. 1 and 50 microns. Also no more than Nt_c_max drops total.
2771 xrc=MAX(R1, (qc1d(k) + qcten(k)*dtsave)*rho(k))
2772 xnc=MAX(2., (nc1d(k) + ncten(k)*dtsave)*rho(k))
2773 if (xrc .gt. R1) then
2774 nu_c = MIN(15, NINT(1000.E6/xnc) + 2)
2775 lamc = (xnc*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
2776 xDc = (bm_r + nu_c + 1.) / lamc
2777 if (xDc.lt. D0c) then
2778 lamc = cce(2,nu_c)/D0c
2779 xnc = ccg(1,nu_c)*ocg2(nu_c)*xrc/am_r*lamc**bm_r
2780 ncten(k) = (xnc-nc1d(k)*rho(k))*odts*orho
2781 elseif (xDc.gt. D0r*2.) then
2782 lamc = cce(2,nu_c)/(D0r*2.)
2783 xnc = ccg(1,nu_c)*ocg2(nu_c)*xrc/am_r*lamc**bm_r
2784 ncten(k) = (xnc-nc1d(k)*rho(k))*odts*orho
2787 ncten(k) = -nc1d(k)*odts
2789 xnc=MAX(0.,(nc1d(k) + ncten(k)*dtsave)*rho(k))
2790 if (xnc.gt.Nt_c_max) &
2791 ncten(k) = (Nt_c_max-nc1d(k)*rho(k))*odts*orho
2793 !..Cloud ice mixing ratio tendency
2794 qiten(k) = qiten(k) + (pri_inu(k) + pri_iha(k) + pri_ihm(k) &
2795 + pri_wfz(k) + pri_rfz(k) + pri_ide(k) &
2796 - prs_iau(k) - prs_sci(k) - pri_rci(k)) &
2799 !..Cloud ice number tendency.
2800 niten(k) = niten(k) + (pni_inu(k) + pni_iha(k) + pni_ihm(k) &
2801 + pni_wfz(k) + pni_rfz(k) + pni_ide(k) &
2802 - pni_iau(k) - pni_sci(k) - pni_rci(k)) &
2805 !..Cloud ice mass/number balance; keep mass-wt mean size between
2806 !.. 5 and 300 microns. Also no more than 500 xtals per liter.
2807 xri=MAX(R1,(qi1d(k) + qiten(k)*dtsave)*rho(k))
2808 xni=MAX(R2,(ni1d(k) + niten(k)*dtsave)*rho(k))
2809 if (xri.gt. R1) then
2810 lami = (am_i*cig(2)*oig1*xni/xri)**obmi
2812 xDi = (bm_i + mu_i + 1.) * ilami
2813 if (xDi.lt. 5.E-6) then
2815 xni = MIN(999.D3, cig(1)*oig2*xri/am_i*lami**bm_i)
2816 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
2817 elseif (xDi.gt. 300.E-6) then
2818 lami = cie(2)/300.E-6
2819 xni = cig(1)*oig2*xri/am_i*lami**bm_i
2820 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
2823 niten(k) = -ni1d(k)*odts
2825 xni=MAX(0.,(ni1d(k) + niten(k)*dtsave)*rho(k))
2826 if (xni.gt.999.E3) &
2827 niten(k) = (999.E3-ni1d(k)*rho(k))*odts*orho
2830 qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) &
2831 + prr_sml(k) + prr_gml(k) + prr_rcs(k) &
2832 + prr_rcg(k) - prg_rfz(k) &
2833 - pri_rfz(k) - prr_rci(k)) &
2836 !..Rain number tendency
2837 nrten(k) = nrten(k) + (pnr_wau(k) + pnr_sml(k) + pnr_gml(k) &
2838 - (pnr_rfz(k) + pnr_rcr(k) + pnr_rcg(k) &
2839 + pnr_rcs(k) + pnr_rci(k) + pni_rfz(k)) ) &
2842 !..Rain mass/number balance; keep median volume diameter between
2843 !.. 37 microns (D0r*0.75) and 2.5 mm.
2844 xrr=MAX(R1,(qr1d(k) + qrten(k)*dtsave)*rho(k))
2845 xnr=MAX(R2,(nr1d(k) + nrten(k)*dtsave)*rho(k))
2846 if (xrr.gt. R1) then
2847 lamr = (am_r*crg(3)*org2*xnr/xrr)**obmr
2848 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
2849 if (mvd_r(k) .gt. 2.5E-3) then
2851 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
2852 xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
2853 nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
2854 elseif (mvd_r(k) .lt. D0r*0.75) then
2856 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
2857 xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
2858 nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
2861 qrten(k) = -qr1d(k)*odts
2862 nrten(k) = -nr1d(k)*odts
2866 qsten(k) = qsten(k) + (prs_iau(k) + prs_sde(k) &
2867 + prs_sci(k) + prs_scw(k) + prs_rcs(k) &
2868 + prs_ide(k) - prs_ihm(k) - prr_sml(k)) &
2872 qgten(k) = qgten(k) + (prg_scw(k) + prg_rfz(k) &
2873 + prg_gde(k) + prg_rcg(k) + prg_gcw(k) &
2874 + prg_rci(k) + prg_rcs(k) - prg_ihm(k) &
2878 !..Temperature tendency
2879 if (temp(k).lt.T_0) then
2881 + ( lsub*ocp(k)*(pri_inu(k) + pri_ide(k) &
2882 + prs_ide(k) + prs_sde(k) &
2883 + prg_gde(k) + pri_iha(k)) &
2884 + lfus2*ocp(k)*(pri_wfz(k) + pri_rfz(k) &
2885 + prg_rfz(k) + prs_scw(k) &
2886 + prg_scw(k) + prg_gcw(k) &
2887 + prg_rcs(k) + prs_rcs(k) &
2888 + prr_rci(k) + prg_rcg(k)) &
2892 + ( lfus*ocp(k)*(-prr_sml(k) - prr_gml(k) &
2893 - prr_rcg(k) - prr_rcs(k)) &
2894 + lsub*ocp(k)*(prs_sde(k) + prg_gde(k)) &
2900 !+---+-----------------------------------------------------------------+
2901 !..Update variables for TAU+1 before condensation & sedimention.
2902 !+---+-----------------------------------------------------------------+
2904 temp(k) = t1d(k) + DT*tten(k)
2906 tempc = temp(k) - 273.15
2907 qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
2908 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
2909 rhof(k) = SQRT(RHO_NOT/rho(k))
2910 rhof2(k) = SQRT(rhof(k))
2911 qvs(k) = rslf(pres(k), temp(k))
2912 ssatw(k) = qv(k)/qvs(k) - 1.
2913 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
2914 diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
2915 if (tempc .ge. 0.0) then
2916 visco(k) = (1.718+0.0049*tempc)*1.0E-5
2918 visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
2920 vsc2(k) = SQRT(rho(k)/visco(k))
2921 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
2922 tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
2923 ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
2924 lvt2(k)=lvap(k)*lvap(k)*ocp(k)*oRv*otemp*otemp
2926 nwfa(k) = MAX(11.1E6, (nwfa1d(k) + nwfaten(k)*DT)*rho(k))
2930 if ((qc1d(k) + qcten(k)*DT) .gt. R1) then
2931 rc(k) = (qc1d(k) + qcten(k)*DT)*rho(k)
2932 nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max))
2933 if (.NOT. is_aerosol_aware) nc(k) = Nt_c
2941 if ((qi1d(k) + qiten(k)*DT) .gt. R1) then
2942 ri(k) = (qi1d(k) + qiten(k)*DT)*rho(k)
2943 ni(k) = MAX(R2, (ni1d(k) + niten(k)*DT)*rho(k))
2951 if ((qr1d(k) + qrten(k)*DT) .gt. R1) then
2952 rr(k) = (qr1d(k) + qrten(k)*DT)*rho(k)
2953 nr(k) = MAX(R2, (nr1d(k) + nrten(k)*DT)*rho(k))
2955 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
2956 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
2957 if (mvd_r(k) .gt. 2.5E-3) then
2959 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
2960 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
2961 elseif (mvd_r(k) .lt. D0r*0.75) then
2963 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
2964 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
2972 if ((qs1d(k) + qsten(k)*DT) .gt. R1) then
2973 rs(k) = (qs1d(k) + qsten(k)*DT)*rho(k)
2980 if ((qg1d(k) + qgten(k)*DT) .gt. R1) then
2981 rg(k) = (qg1d(k) + qgten(k)*DT)*rho(k)
2989 !+---+-----------------------------------------------------------------+
2990 !..With tendency-updated mixing ratios, recalculate snow moments and
2991 !.. intercepts/slopes of graupel and rain.
2992 !+---+-----------------------------------------------------------------+
2993 if (.not. iiwarm) then
3001 if (.not. L_qs(k)) CYCLE
3002 tc0 = MIN(-0.1, temp(k)-273.15)
3003 smob(k) = rs(k)*oams
3005 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
3006 !.. then we must compute actual 2nd moment and use as reference.
3007 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
3010 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
3011 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
3012 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
3013 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
3014 + sa(10)*bm_s*bm_s*bm_s
3016 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
3017 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
3018 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
3019 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
3020 + sb(10)*bm_s*bm_s*bm_s
3021 smo2(k) = (smob(k)/a_)**(1./b_)
3024 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
3025 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
3026 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
3027 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
3028 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
3029 + sa(10)*cse(1)*cse(1)*cse(1)
3031 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
3032 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
3033 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
3034 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
3035 smoc(k) = a_ * smo2(k)**b_
3037 !..Calculate bm_s+bv_s (th) moment. Useful for sedimentation.
3038 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(14) &
3039 + sa(4)*tc0*cse(14) + sa(5)*tc0*tc0 &
3040 + sa(6)*cse(14)*cse(14) + sa(7)*tc0*tc0*cse(14) &
3041 + sa(8)*tc0*cse(14)*cse(14) + sa(9)*tc0*tc0*tc0 &
3042 + sa(10)*cse(14)*cse(14)*cse(14)
3044 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(14) + sb(4)*tc0*cse(14) &
3045 + sb(5)*tc0*tc0 + sb(6)*cse(14)*cse(14) &
3046 + sb(7)*tc0*tc0*cse(14) + sb(8)*tc0*cse(14)*cse(14) &
3047 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(14)*cse(14)*cse(14)
3048 smod(k) = a_ * smo2(k)**b_
3051 !+---+-----------------------------------------------------------------+
3052 !..Calculate y-intercept, slope values for graupel.
3053 !+---+-----------------------------------------------------------------+
3056 ygra1 = alog10(max(1.E-9, rg(k)))
3057 zans1 = 3.0 + 2./7.*(ygra1+8.)
3058 zans1 = MAX(2., MIN(zans1, 6.))
3059 N0_exp = 10.**(zans1)
3060 lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
3061 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
3063 N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2)
3068 !+---+-----------------------------------------------------------------+
3069 !..Calculate y-intercept, slope values for rain.
3070 !+---+-----------------------------------------------------------------+
3072 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3074 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3075 N0_r(k) = nr(k)*org2*lamr**cre(2)
3078 !+---+-----------------------------------------------------------------+
3079 !..Cloud water condensation and evaporation. Nucleate cloud droplets
3080 !.. using explicit CCN aerosols with hygroscopicity like sulfates using
3081 !.. parcel model lookup table results provided by T. Eidhammer. Evap
3082 !.. drops using calculation of max drop size capable of evaporating in
3083 !.. single timestep and explicit number of drops smaller than Dc_star
3084 !.. from lookup table.
3085 !+---+-----------------------------------------------------------------+
3088 if ( (ssatw(k).gt. eps) .or. (ssatw(k).lt. -eps .and. &
3090 clap = (qv(k)-qvs(k))/(1. + lvt2(k)*qvs(k))
3092 fcd = qvs(k)* EXP(lvt2(k)*clap) - qv(k) + clap
3093 dfcd = qvs(k)*lvt2(k)* EXP(lvt2(k)*clap) + 1.
3094 clap = clap - fcd/dfcd
3096 xrc = rc(k) + clap*rho(k)
3098 if (xrc.gt. R1) then
3099 prw_vcd(k) = clap*odt
3100 !+---+-----------------------------------------------------------------+ ! DROPLET NUCLEATION
3101 if (clap .gt. eps) then
3102 if (is_aerosol_aware) then
3103 xnc = MAX(2., activ_ncloud(temp(k), w1d(k), nwfa(k)))
3107 pnc_wcd(k) = 0.5*(xnc-nc(k) + abs(xnc-nc(k)))*odts*orho
3109 !+---+-----------------------------------------------------------------+ ! EVAPORATION
3110 elseif (clap .lt. -eps .AND. ssatw(k).lt.-1.E-6 .AND. is_aerosol_aware) then
3111 tempc = temp(k) - 273.15
3114 rvs_p = rvs*otemp*(lvap(k)*otemp*oRv - 1.)
3115 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*oRv - 1.) &
3116 *otemp*(lvap(k)*otemp*oRv - 1.) &
3117 + (-2.*lvap(k)*otemp*otemp*otemp*oRv) &
3119 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
3120 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
3121 * rvs_pp/rvs_p * rvs/rvs_p
3122 alphsc = MAX(1.E-9, alphsc)
3124 if (abs(xsat).lt. 1.E-9) xsat=0.
3125 t1_evap = 2.*PI*( 1.0 - alphsc*xsat &
3126 + 2.*alphsc*alphsc*xsat*xsat &
3127 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
3130 Dc_star = DSQRT(-2.D0*DT * t1_evap/(2.*PI) &
3131 * 4.*diffu(k)*ssatw(k)*rvs/rho_w)
3132 idx_d = MAX(1, MIN(INT(1.E6*Dc_star), nbc))
3134 idx_n = NINT(1.0 + FLOAT(nbc) * DLOG(nc(k)/t_Nc(1)) / nic1)
3135 idx_n = MAX(1, MIN(idx_n, nbc))
3137 !..Cloud water lookup table index.
3138 if (rc(k).gt. r_c(1)) then
3139 nic = NINT(ALOG10(rc(k)))
3140 do nn = nic-1, nic+1
3142 if ( (rc(k)/10.**nn).ge.1.0 .and. &
3143 (rc(k)/10.**nn).lt.10.0) goto 159
3146 idx_c = INT(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
3147 idx_c = MAX(1, MIN(idx_c, ntb_c))
3152 !prw_vcd(k) = MAX(DBLE(-rc(k)*orho*odt), &
3153 ! -tpc_wev(idx_d, idx_c, idx_n)*orho*odt)
3154 prw_vcd(k) = MAX(DBLE(-rc(k)*0.99*orho*odt), prw_vcd(k))
3155 pnc_wcd(k) = MAX(DBLE(-nc(k)*0.99*orho*odt), &
3156 DBLE(-tnc_wev(idx_d, idx_c, idx_n)*orho*odt))
3160 prw_vcd(k) = -rc(k)*orho*odt
3161 pnc_wcd(k) = -nc(k)*orho*odt
3164 !+---+-----------------------------------------------------------------+
3166 qvten(k) = qvten(k) - prw_vcd(k)
3167 qcten(k) = qcten(k) + prw_vcd(k)
3168 ncten(k) = ncten(k) + pnc_wcd(k)
3169 nwfaten(k) = nwfaten(k) - pnc_wcd(k)
3170 tten(k) = tten(k) + lvap(k)*ocp(k)*prw_vcd(k)*(1-IFDRY)
3171 rc(k) = MAX(R1, (qc1d(k) + DT*qcten(k))*rho(k))
3172 if (rc(k).eq.R1) L_qc(k) = .false.
3173 nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max))
3174 if (.NOT. is_aerosol_aware) nc(k) = Nt_c
3175 qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
3176 temp(k) = t1d(k) + DT*tten(k)
3177 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
3178 qvs(k) = rslf(pres(k), temp(k))
3179 ssatw(k) = qv(k)/qvs(k) - 1.
3183 !+---+-----------------------------------------------------------------+
3184 !.. If still subsaturated, allow rain to evaporate, following
3185 !.. Srivastava & Coen (1992).
3186 !+---+-----------------------------------------------------------------+
3188 if ( (ssatw(k).lt. -eps) .and. L_qr(k) &
3189 .and. (.not.(prw_vcd(k).gt. 0.)) ) then
3190 tempc = temp(k) - 273.15
3193 rhof(k) = SQRT(RHO_NOT*orho)
3194 rhof2(k) = SQRT(rhof(k))
3195 diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
3196 if (tempc .ge. 0.0) then
3197 visco(k) = (1.718+0.0049*tempc)*1.0E-5
3199 visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5
3201 vsc2(k) = SQRT(rho(k)/visco(k))
3202 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
3203 tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936
3204 ocp(k) = 1./(Cp*(1.+0.887*qv(k)))
3207 rvs_p = rvs*otemp*(lvap(k)*otemp*oRv - 1.)
3208 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*oRv - 1.) &
3209 *otemp*(lvap(k)*otemp*oRv - 1.) &
3210 + (-2.*lvap(k)*otemp*otemp*otemp*oRv) &
3212 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
3213 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
3214 * rvs_pp/rvs_p * rvs/rvs_p
3215 alphsc = MAX(1.E-9, alphsc)
3216 xsat = MIN(-1.E-9, ssatw(k))
3217 t1_evap = 2.*PI*( 1.0 - alphsc*xsat &
3218 + 2.*alphsc*alphsc*xsat*xsat &
3219 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
3223 !..Rapidly eliminate near zero values when low humidity (<95%)
3224 if (qv(k)/qvs(k) .lt. 0.95 .AND. rr(k)*orho.le.1.E-8) then
3225 prv_rev(k) = rr(k)*orho*odts
3227 prv_rev(k) = t1_evap*diffu(k)*(-ssatw(k))*N0_r(k)*rvs &
3228 * (t1_qr_ev*ilamr(k)**cre(10) &
3229 + t2_qr_ev*vsc2(k)*rhof2(k)*((lamr+0.5*fv_r)**(-cre(11))))
3230 rate_max = MIN((rr(k)*orho*odts), (qvs(k)-qv(k))*odts)
3231 prv_rev(k) = MIN(DBLE(rate_max), prv_rev(k)*orho)
3233 !..TEST: G. Thompson 10 May 2013
3234 !..Reduce the rain evaporation in same places as melting graupel occurs.
3235 !..Rationale: falling and simultaneous melting graupel in subsaturated
3236 !..regions will not melt as fast because particle temperature stays
3237 !..at 0C. Also not much shedding of the water from the graupel so
3238 !..likely that the water-coated graupel evaporating much slower than
3239 !..if the water was immediately shed off.
3240 IF (prr_gml(k).gt.0.0) THEN
3241 eva_factor = MIN(1.0, 0.01+(0.99-0.01)*(tempc/20.0))
3242 prv_rev(k) = prv_rev(k)*eva_factor
3246 pnr_rev(k) = MIN(DBLE(nr(k)*0.99*orho*odts), & ! RAIN2M
3247 prv_rev(k) * nr(k)/rr(k))
3249 qrten(k) = qrten(k) - prv_rev(k)
3250 qvten(k) = qvten(k) + prv_rev(k)
3251 nrten(k) = nrten(k) - pnr_rev(k)
3252 nwfaten(k) = nwfaten(k) + pnr_rev(k)
3253 tten(k) = tten(k) - lvap(k)*ocp(k)*prv_rev(k)*(1-IFDRY)
3255 rr(k) = MAX(R1, (qr1d(k) + DT*qrten(k))*rho(k))
3256 qv(k) = MAX(1.E-10, qv1d(k) + DT*qvten(k))
3257 nr(k) = MAX(R2, (nr1d(k) + DT*nrten(k))*rho(k))
3258 temp(k) = t1d(k) + DT*tten(k)
3259 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
3262 #if ( WRF_CHEM == 1 )
3263 if( wetscav_on ) then
3265 evapprod(k) = prv_rev(k) - (min(zeroD0,prs_sde(k)) + &
3266 min(zeroD0,prg_gde(k)))
3267 rainprod(k) = prr_wau(k) + prr_rcw(k) + prs_scw(k) + &
3268 prg_scw(k) + prs_iau(k) + &
3269 prg_gcw(k) + prs_sci(k) + &
3275 !+---+-----------------------------------------------------------------+
3276 !..Find max terminal fallspeed (distribution mass-weighted mean
3277 !.. velocity) and use it to determine if we need to split the timestep
3278 !.. (var nstep>1). Either way, only bother to do sedimentation below
3279 !.. 1st level that contains any sedimenting particles (k=ksed1 on down).
3280 !.. New in v3.0+ is computing separate for rain, ice, snow, and
3281 !.. graupel species thus making code faster with credit to J. Schmidt.
3282 !+---+-----------------------------------------------------------------+
3286 do k = kte+1, kts, -1
3297 if (ANY(L_qr .eqv. .true.)) then
3300 rhof(k) = SQRT(RHO_NOT/rho(k))
3302 if (rr(k).gt. R1) then
3303 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
3304 vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) &
3305 *((lamr+fv_r)**(-cre(6)))
3307 ! First below is technically correct:
3308 ! vtr = rhof(k)*av_r*crg(5)*org2 * lamr**cre(2) &
3309 ! *((lamr+fv_r)**(-cre(5)))
3310 ! Test: make number fall faster (but still slower than mass)
3311 ! Goal: less prominent size sorting
3312 vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) &
3313 *((lamr+fv_r)**(-cre(7)))
3317 vtnrk(k) = vtnrk(k+1)
3320 if (MAX(vtrk(k),vtnrk(k)) .gt. 1.E-3) then
3321 ksed1(1) = MAX(ksed1(1), k)
3322 delta_tp = dzq(k)/(MAX(vtrk(k),vtnrk(k)))
3323 nstep = MAX(nstep, INT(DT/delta_tp + 1.))
3326 if (ksed1(1) .eq. kte) ksed1(1) = kte-1
3327 if (nstep .gt. 0) onstep(1) = 1./REAL(nstep)
3330 !+---+-----------------------------------------------------------------+
3332 if (ANY(L_qc .eqv. .true.)) then
3335 if (rc(k) .gt. R2) ksed1(5) = k
3336 hgt_agl = hgt_agl + dzq(k)
3337 if (hgt_agl .gt. 500.0) goto 151
3341 do k = ksed1(5), kts, -1
3343 if (rc(k) .gt. R1 .and. w1d(k) .lt. 1.E-1) then
3344 nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2)
3345 lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
3347 vtc = rhof(k)*av_c*ccg(5,nu_c)*ocg2(nu_c) * ilamc**bv_c
3349 vtc = rhof(k)*av_c*ccg(4,nu_c)*ocg1(nu_c) * ilamc**bv_c
3355 !+---+-----------------------------------------------------------------+
3357 if (.not. iiwarm) then
3359 if (ANY(L_qi .eqv. .true.)) then
3364 if (ri(k).gt. R1) then
3365 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
3367 vti = rhof(k)*av_i*cig(3)*oig2 * ilami**bv_i
3369 ! First below is technically correct:
3370 ! vti = rhof(k)*av_i*cig(4)*oig1 * ilami**bv_i
3371 ! Goal: less prominent size sorting
3372 vti = rhof(k)*av_i*cig(6)/cig(7) * ilami**bv_i
3376 vtnik(k) = vtnik(k+1)
3379 if (vtik(k) .gt. 1.E-3) then
3380 ksed1(2) = MAX(ksed1(2), k)
3381 delta_tp = dzq(k)/vtik(k)
3382 nstep = MAX(nstep, INT(DT/delta_tp + 1.))
3385 if (ksed1(2) .eq. kte) ksed1(2) = kte-1
3386 if (nstep .gt. 0) onstep(2) = 1./REAL(nstep)
3389 !+---+-----------------------------------------------------------------+
3391 if (ANY(L_qs .eqv. .true.)) then
3396 if (rs(k).gt. R1) then
3397 xDs = smoc(k) / smob(k)
3399 ils1 = 1./(Mrat*Lam0 + fv_s)
3400 ils2 = 1./(Mrat*Lam1 + fv_s)
3401 t1_vts = Kap0*csg(4)*ils1**cse(4)
3402 t2_vts = Kap1*Mrat**mu_s*csg(10)*ils2**cse(10)
3403 ils1 = 1./(Mrat*Lam0)
3404 ils2 = 1./(Mrat*Lam1)
3405 t3_vts = Kap0*csg(1)*ils1**cse(1)
3406 t4_vts = Kap1*Mrat**mu_s*csg(7)*ils2**cse(7)
3407 vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts)
3408 if (temp(k).gt. (T_0+0.1)) then
3409 SR = rs(k)/(rs(k)+rr(k))
3410 vtsk(k) = vts*SR + (1.-SR)*vtrk(k)
3412 vtsk(k) = vts*vts_boost(k)
3418 if (vtsk(k) .gt. 1.E-3) then
3419 ksed1(3) = MAX(ksed1(3), k)
3420 delta_tp = dzq(k)/vtsk(k)
3421 nstep = MAX(nstep, INT(DT/delta_tp + 1.))
3424 if (ksed1(3) .eq. kte) ksed1(3) = kte-1
3425 if (nstep .gt. 0) onstep(3) = 1./REAL(nstep)
3428 !+---+-----------------------------------------------------------------+
3430 if (ANY(L_qg .eqv. .true.)) then
3435 if (rg(k).gt. R1) then
3436 vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g
3437 if (temp(k).gt. T_0) then
3438 vtgk(k) = MAX(vtg, vtrk(k))
3446 if (vtgk(k) .gt. 1.E-3) then
3447 ksed1(4) = MAX(ksed1(4), k)
3448 delta_tp = dzq(k)/vtgk(k)
3449 nstep = MAX(nstep, INT(DT/delta_tp + 1.))
3452 if (ksed1(4) .eq. kte) ksed1(4) = kte-1
3453 if (nstep .gt. 0) onstep(4) = 1./REAL(nstep)
3457 !+---+-----------------------------------------------------------------+
3458 !..Sedimentation of mixing ratio is the integral of v(D)*m(D)*N(D)*dD,
3459 !.. whereas neglect m(D) term for number concentration. Therefore,
3460 !.. cloud ice has proper differential sedimentation.
3461 !+---+-----------------------------------------------------------------+
3463 if (ANY(L_qr .eqv. .true.)) then
3464 nstep = NINT(1./onstep(1))
3467 sed_r(k) = vtrk(k)*rr(k)
3468 sed_n(k) = vtnrk(k)*nr(k)
3473 qrten(k) = qrten(k) - sed_r(k)*odzq*onstep(1)*orho
3474 nrten(k) = nrten(k) - sed_n(k)*odzq*onstep(1)*orho
3475 rr(k) = MAX(R1, rr(k) - sed_r(k)*odzq*DT*onstep(1))
3476 nr(k) = MAX(R2, nr(k) - sed_n(k)*odzq*DT*onstep(1))
3477 do k = ksed1(1), kts, -1
3480 qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) &
3481 *odzq*onstep(1)*orho
3482 nrten(k) = nrten(k) + (sed_n(k+1)-sed_n(k)) &
3483 *odzq*onstep(1)*orho
3484 rr(k) = MAX(R1, rr(k) + (sed_r(k+1)-sed_r(k)) &
3486 nr(k) = MAX(R2, nr(k) + (sed_n(k+1)-sed_n(k)) &
3490 if (rr(kts).gt.R1*1000.) &
3491 pptrain = pptrain + sed_r(kts)*DT*onstep(1)
3495 !+---+-----------------------------------------------------------------+
3497 if (ANY(L_qc .eqv. .true.)) then
3499 sed_c(k) = vtck(k)*rc(k)
3500 sed_n(k) = vtnck(k)*nc(k)
3502 do k = ksed1(5), kts, -1
3505 qcten(k) = qcten(k) + (sed_c(k+1)-sed_c(k)) *odzq*orho
3506 ncten(k) = ncten(k) + (sed_n(k+1)-sed_n(k)) *odzq*orho
3507 rc(k) = MAX(R1, rc(k) + (sed_c(k+1)-sed_c(k)) *odzq*DT)
3508 nc(k) = MAX(10., nc(k) + (sed_n(k+1)-sed_n(k)) *odzq*DT)
3512 !+---+-----------------------------------------------------------------+
3514 if (ANY(L_qi .eqv. .true.)) then
3515 nstep = NINT(1./onstep(2))
3518 sed_i(k) = vtik(k)*ri(k)
3519 sed_n(k) = vtnik(k)*ni(k)
3524 qiten(k) = qiten(k) - sed_i(k)*odzq*onstep(2)*orho
3525 niten(k) = niten(k) - sed_n(k)*odzq*onstep(2)*orho
3526 ri(k) = MAX(R1, ri(k) - sed_i(k)*odzq*DT*onstep(2))
3527 ni(k) = MAX(R2, ni(k) - sed_n(k)*odzq*DT*onstep(2))
3528 do k = ksed1(2), kts, -1
3531 qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) &
3532 *odzq*onstep(2)*orho
3533 niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) &
3534 *odzq*onstep(2)*orho
3535 ri(k) = MAX(R1, ri(k) + (sed_i(k+1)-sed_i(k)) &
3537 ni(k) = MAX(R2, ni(k) + (sed_n(k+1)-sed_n(k)) &
3541 if (ri(kts).gt.R1*1000.) &
3542 pptice = pptice + sed_i(kts)*DT*onstep(2)
3546 !+---+-----------------------------------------------------------------+
3548 if (ANY(L_qs .eqv. .true.)) then
3549 nstep = NINT(1./onstep(3))
3552 sed_s(k) = vtsk(k)*rs(k)
3557 qsten(k) = qsten(k) - sed_s(k)*odzq*onstep(3)*orho
3558 rs(k) = MAX(R1, rs(k) - sed_s(k)*odzq*DT*onstep(3))
3559 do k = ksed1(3), kts, -1
3562 qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) &
3563 *odzq*onstep(3)*orho
3564 rs(k) = MAX(R1, rs(k) + (sed_s(k+1)-sed_s(k)) &
3568 if (rs(kts).gt.R1*1000.) &
3569 pptsnow = pptsnow + sed_s(kts)*DT*onstep(3)
3573 !+---+-----------------------------------------------------------------+
3575 if (ANY(L_qg .eqv. .true.)) then
3576 nstep = NINT(1./onstep(4))
3579 sed_g(k) = vtgk(k)*rg(k)
3584 qgten(k) = qgten(k) - sed_g(k)*odzq*onstep(4)*orho
3585 rg(k) = MAX(R1, rg(k) - sed_g(k)*odzq*DT*onstep(4))
3586 do k = ksed1(4), kts, -1
3589 qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) &
3590 *odzq*onstep(4)*orho
3591 rg(k) = MAX(R1, rg(k) + (sed_g(k+1)-sed_g(k)) &
3595 if (rg(kts).gt.R1*1000.) &
3596 pptgraul = pptgraul + sed_g(kts)*DT*onstep(4)
3600 !+---+-----------------------------------------------------------------+
3601 !.. Instantly melt any cloud ice into cloud water if above 0C and
3602 !.. instantly freeze any cloud water found below HGFR.
3603 !+---+-----------------------------------------------------------------+
3604 if (.not. iiwarm) then
3606 xri = MAX(0.0, qi1d(k) + qiten(k)*DT)
3607 if ( (temp(k).gt. T_0) .and. (xri.gt. 0.0) ) then
3608 qcten(k) = qcten(k) + xri*odt
3609 ncten(k) = ncten(k) + ni1d(k)*odt
3610 qiten(k) = qiten(k) - xri*odt
3611 niten(k) = -ni1d(k)*odt
3612 tten(k) = tten(k) - lfus*ocp(k)*xri*odt*(1-IFDRY)
3615 xrc = MAX(0.0, qc1d(k) + qcten(k)*DT)
3616 if ( (temp(k).lt. HGFR) .and. (xrc.gt. 0.0) ) then
3617 lfus2 = lsub - lvap(k)
3618 xnc = nc1d(k) + ncten(k)*DT
3619 qiten(k) = qiten(k) + xrc*odt
3620 niten(k) = niten(k) + xnc*odt
3621 qcten(k) = qcten(k) - xrc*odt
3622 ncten(k) = ncten(k) - xnc*odt
3623 tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-IFDRY)
3628 !+---+-----------------------------------------------------------------+
3629 !.. All tendencies computed, apply and pass back final values to parent.
3630 !+---+-----------------------------------------------------------------+
3632 t1d(k) = t1d(k) + tten(k)*DT
3633 qv1d(k) = MAX(1.E-10, qv1d(k) + qvten(k)*DT)
3634 qc1d(k) = qc1d(k) + qcten(k)*DT
3635 nc1d(k) = MAX(2./rho(k), MIN(nc1d(k) + ncten(k)*DT, Nt_c_max))
3636 if (is_aerosol_aware) then
3637 if (aer_init_opt .lt. 2) then
3638 nwfa1d(k) = MAX(11.1E6, MIN(9999.E6, &
3639 (nwfa1d(k)+nwfaten(k)*DT)))
3640 nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6, &
3641 (nifa1d(k)+nifaten(k)*DT)))
3642 if (wif_input_opt .eq. 2) then
3643 nbca1d(k) = MAX(5.55E6, MIN(9999.E6, &
3644 (nbca1d(k)+nbcaten(k)*DT)))
3649 nwfa1d(k) = MAX(0.0, (nwfa1d(k)+nwfaten(k)*DT))
3650 nifa1d(k) = MAX(0.0, (nifa1d(k)+nifaten(k)*DT))
3651 if (wif_input_opt .eq. 2) then
3652 nbca1d(k) = MAX(0.0, (nbca1d(k)+nbcaten(k)*DT))
3658 nwfa1d(k) = MAX(11.1E6, MIN(9999.E6, &
3659 (nwfa1d(k)+nwfaten(k)*DT)))
3660 nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6, &
3661 (nifa1d(k)+nifaten(k)*DT)))
3662 nbca1d(k) = MAX(5.55E6, MIN(9999.E6, &
3663 (nbca1d(k)+nbcaten(k)*DT)))
3666 if (qc1d(k) .le. R1) then
3670 nu_c = MIN(15, NINT(1000.E6/(nc1d(k)*rho(k))) + 2)
3671 lamc = (am_r*ccg(2,nu_c)*ocg1(nu_c)*nc1d(k)/qc1d(k))**obmr
3672 xDc = (bm_r + nu_c + 1.) / lamc
3673 if (xDc.lt. D0c) then
3674 lamc = cce(2,nu_c)/D0c
3675 elseif (xDc.gt. D0r*2.) then
3676 lamc = cce(2,nu_c)/(D0r*2.)
3678 nc1d(k) = MIN(ccg(1,nu_c)*ocg2(nu_c)*qc1d(k)/am_r*lamc**bm_r,&
3679 DBLE(Nt_c_max)/rho(k))
3682 qi1d(k) = qi1d(k) + qiten(k)*DT
3683 ni1d(k) = MAX(R2/rho(k), ni1d(k) + niten(k)*DT)
3684 if (qi1d(k) .le. R1) then
3688 lami = (am_i*cig(2)*oig1*ni1d(k)/qi1d(k))**obmi
3690 xDi = (bm_i + mu_i + 1.) * ilami
3691 if (xDi.lt. 5.E-6) then
3693 elseif (xDi.gt. 300.E-6) then
3694 lami = cie(2)/300.E-6
3696 ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, &
3699 qr1d(k) = qr1d(k) + qrten(k)*DT
3700 nr1d(k) = MAX(R2/rho(k), nr1d(k) + nrten(k)*DT)
3701 if (qr1d(k) .le. R1) then
3705 lamr = (am_r*crg(3)*org2*nr1d(k)/qr1d(k))**obmr
3706 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
3707 if (mvd_r(k) .gt. 2.5E-3) then
3709 elseif (mvd_r(k) .lt. D0r*0.75) then
3712 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
3713 nr1d(k) = crg(2)*org3*qr1d(k)*lamr**bm_r / am_r
3715 qs1d(k) = qs1d(k) + qsten(k)*DT
3716 if (qs1d(k) .le. R1) qs1d(k) = 0.0
3717 qg1d(k) = qg1d(k) + qgten(k)*DT
3718 if (qg1d(k) .le. R1) qg1d(k) = 0.0
3721 end subroutine mp_thompson
3722 !+---+-----------------------------------------------------------------+
3724 !+---+-----------------------------------------------------------------+
3725 !..Creation of the lookup tables and support functions found below here.
3726 !+---+-----------------------------------------------------------------+
3727 !..Rain collecting graupel (and inverse). Explicit CE integration.
3728 !+---+-----------------------------------------------------------------+
3730 subroutine qr_acr_qg
3735 INTEGER:: i, j, k, m, n, n2
3736 INTEGER:: km, km_s, km_e
3737 DOUBLE PRECISION, DIMENSION(nbg):: vg, N_g
3738 DOUBLE PRECISION, DIMENSION(nbr):: vr, N_r
3739 DOUBLE PRECISION:: N0_r, N0_g, lam_exp, lamg, lamr
3740 DOUBLE PRECISION:: massg, massr, dvg, dvr, t1, t2, z1, z2, y1, y2
3741 LOGICAL force_read_thompson, write_thompson_tables
3742 LOGICAL lexist,lopen
3744 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
3748 CALL nl_get_force_read_thompson(1,force_read_thompson)
3749 CALL nl_get_write_thompson_tables(1,write_thompson_tables)
3752 IF ( wrf_dm_on_monitor() ) THEN
3753 INQUIRE(FILE="qr_acr_qgV3.dat",EXIST=lexist)
3755 CALL wrf_message("ThompMP: read qr_acr_qgV3.dat instead of computing")
3756 OPEN(63,file="qr_acr_qgV3.dat",form="unformatted",err=1234)
3757 READ(63,err=1234) tcg_racg
3758 READ(63,err=1234) tmr_racg
3759 READ(63,err=1234) tcr_gacr
3760 READ(63,err=1234) tmg_gacr
3761 READ(63,err=1234) tnr_racg
3762 READ(63,err=1234) tnr_gacr
3765 IF ( good .NE. 1 ) THEN
3766 INQUIRE(63,opened=lopen)
3768 IF( force_read_thompson ) THEN
3769 CALL wrf_error_fatal("Error reading qr_acr_qgV3.dat. Aborting because force_read_thompson is .true.")
3773 IF( force_read_thompson ) THEN
3774 CALL wrf_error_fatal("Error opening qr_acr_qgV3.dat. Aborting because force_read_thompson is .true.")
3778 INQUIRE(63,OPENED=lopen)
3784 IF( force_read_thompson ) THEN
3785 CALL wrf_error_fatal("Non-existent qr_acr_qgV3.dat. Aborting because force_read_thompson is .true.")
3789 #if defined(DM_PARALLEL) && !defined(STUBMPI)
3790 CALL wrf_dm_bcast_integer(good,1)
3793 IF ( good .EQ. 1 ) THEN
3794 #if defined(DM_PARALLEL) && !defined(STUBMPI)
3795 CALL wrf_dm_bcast_double(tcg_racg,SIZE(tcg_racg))
3796 CALL wrf_dm_bcast_double(tmr_racg,SIZE(tmr_racg))
3797 CALL wrf_dm_bcast_double(tcr_gacr,SIZE(tcr_gacr))
3798 CALL wrf_dm_bcast_double(tmg_gacr,SIZE(tmg_gacr))
3799 CALL wrf_dm_bcast_double(tnr_racg,SIZE(tnr_racg))
3800 CALL wrf_dm_bcast_double(tnr_gacr,SIZE(tnr_gacr))
3803 CALL wrf_message("ThompMP: computing qr_acr_qg")
3805 ! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
3806 vr(n2) = -0.1021 + 4.932E3*Dr(n2) - 0.9551E6*Dr(n2)*Dr(n2) &
3807 + 0.07934E9*Dr(n2)*Dr(n2)*Dr(n2) &
3808 - 0.002362E12*Dr(n2)*Dr(n2)*Dr(n2)*Dr(n2)
3811 vg(n) = av_g*Dg(n)**bv_g
3814 !..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for
3815 !.. fortran indices. J. Michalakes, 2009Oct30.
3817 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
3818 CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e )
3821 km_e = ntb_r*ntb_r1 - 1
3826 k = mod( km , ntb_r1 ) + 1
3828 lam_exp = (N0r_exp(k)*am_r*crg(1)/r_r(m))**ore1
3829 lamr = lam_exp * (crg(3)*org2*org1)**obmr
3830 N0_r = N0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2)
3832 N_r(n2) = N0_r*Dr(n2)**mu_r *DEXP(-lamr*Dr(n2))*dtr(n2)
3837 lam_exp = (N0g_exp(i)*am_g*cgg(1)/r_g(j))**oge1
3838 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
3839 N0_g = N0g_exp(i)/(cgg(2)*lam_exp) * lamg**cge(2)
3841 N_g(n) = N0_g*Dg(n)**mu_g * DEXP(-lamg*Dg(n))*dtg(n)
3851 massr = am_r * Dr(n2)**bm_r
3853 massg = am_g * Dg(n)**bm_g
3855 dvg = 0.5d0*((vr(n2) - vg(n)) + DABS(vr(n2)-vg(n)))
3856 dvr = 0.5d0*((vg(n) - vr(n2)) + DABS(vg(n)-vr(n2)))
3858 t1 = t1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
3859 *dvg*massg * N_g(n)* N_r(n2)
3860 z1 = z1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
3861 *dvg*massr * N_g(n)* N_r(n2)
3862 y1 = y1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
3863 *dvg * N_g(n)* N_r(n2)
3865 t2 = t2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
3866 *dvr*massr * N_g(n)* N_r(n2)
3867 y2 = y2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
3868 *dvr * N_g(n)* N_r(n2)
3869 z2 = z2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) &
3870 *dvr*massg * N_g(n)* N_r(n2)
3874 tcg_racg(i,j,k,m) = t1
3875 tmr_racg(i,j,k,m) = DMIN1(z1, r_r(m)*1.0d0)
3876 tcr_gacr(i,j,k,m) = t2
3877 tmg_gacr(i,j,k,m) = DMIN1(z2, r_g(j)*1.0d0)
3878 !DAVE tmg_gacr(i,j,k,m) = DMIN1(z2, DBLE(r_g(j)))
3879 tnr_racg(i,j,k,m) = y1
3880 tnr_gacr(i,j,k,m) = y2
3885 !..Note wrf_dm_gatherv expects zero-based km_s, km_e (J. Michalakes, 2009Oct30).
3887 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
3888 CALL wrf_dm_gatherv(tcg_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE)
3889 CALL wrf_dm_gatherv(tmr_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE)
3890 CALL wrf_dm_gatherv(tcr_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE)
3891 CALL wrf_dm_gatherv(tmg_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE)
3892 CALL wrf_dm_gatherv(tnr_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE)
3893 CALL wrf_dm_gatherv(tnr_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE)
3896 IF ( write_thompson_tables .AND. wrf_dm_on_monitor() ) THEN
3897 CALL wrf_message("Writing qr_acr_qgV3.dat in Thompson MP init")
3898 OPEN(63,file="qr_acr_qgV3.dat",form="unformatted",err=9234)
3899 WRITE(63,err=9234) tcg_racg
3900 WRITE(63,err=9234) tmr_racg
3901 WRITE(63,err=9234) tcr_gacr
3902 WRITE(63,err=9234) tmg_gacr
3903 WRITE(63,err=9234) tnr_racg
3904 WRITE(63,err=9234) tnr_gacr
3906 RETURN ! ----- RETURN
3908 CALL wrf_error_fatal("Error writing qr_acr_qgV2.dat")
3912 end subroutine qr_acr_qg
3913 !+---+-----------------------------------------------------------------+
3915 !+---+-----------------------------------------------------------------+
3916 !..Rain collecting snow (and inverse). Explicit CE integration.
3917 !+---+-----------------------------------------------------------------+
3919 subroutine qr_acr_qs
3924 INTEGER:: i, j, k, m, n, n2
3925 INTEGER:: km, km_s, km_e
3926 DOUBLE PRECISION, DIMENSION(nbr):: vr, D1, N_r
3927 DOUBLE PRECISION, DIMENSION(nbs):: vs, N_s
3928 DOUBLE PRECISION:: loga_, a_, b_, second, M0, M2, M3, Mrat, oM3
3929 DOUBLE PRECISION:: N0_r, lam_exp, lamr, slam1, slam2
3930 DOUBLE PRECISION:: dvs, dvr, masss, massr
3931 DOUBLE PRECISION:: t1, t2, t3, t4, z1, z2, z3, z4
3932 DOUBLE PRECISION:: y1, y2, y3, y4
3933 LOGICAL force_read_thompson, write_thompson_tables
3934 LOGICAL lexist,lopen
3936 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
3940 CALL nl_get_force_read_thompson(1,force_read_thompson)
3941 CALL nl_get_write_thompson_tables(1,write_thompson_tables)
3944 IF ( wrf_dm_on_monitor() ) THEN
3945 INQUIRE(FILE="qr_acr_qsV2.dat",EXIST=lexist)
3947 CALL wrf_message("ThompMP: read qr_acr_qsV2.dat instead of computing")
3948 OPEN(63,file="qr_acr_qsV2.dat",form="unformatted",err=1234)
3949 READ(63,err=1234)tcs_racs1
3950 READ(63,err=1234)tmr_racs1
3951 READ(63,err=1234)tcs_racs2
3952 READ(63,err=1234)tmr_racs2
3953 READ(63,err=1234)tcr_sacr1
3954 READ(63,err=1234)tms_sacr1
3955 READ(63,err=1234)tcr_sacr2
3956 READ(63,err=1234)tms_sacr2
3957 READ(63,err=1234)tnr_racs1
3958 READ(63,err=1234)tnr_racs2
3959 READ(63,err=1234)tnr_sacr1
3960 READ(63,err=1234)tnr_sacr2
3963 IF ( good .NE. 1 ) THEN
3964 INQUIRE(63,opened=lopen)
3966 IF( force_read_thompson ) THEN
3967 CALL wrf_error_fatal("Error reading qr_acr_qsV2.dat. Aborting because force_read_thompson is .true.")
3971 IF( force_read_thompson ) THEN
3972 CALL wrf_error_fatal("Error opening qr_acr_qsV2.dat. Aborting because force_read_thompson is .true.")
3976 INQUIRE(63,opened=lopen)
3982 IF( force_read_thompson ) THEN
3983 CALL wrf_error_fatal("Non-existent qr_acr_qsV2.dat. Aborting because force_read_thompson is .true.")
3987 #if defined(DM_PARALLEL) && !defined(STUBMPI)
3988 CALL wrf_dm_bcast_integer(good,1)
3991 IF ( good .EQ. 1 ) THEN
3992 #if defined(DM_PARALLEL) && !defined(STUBMPI)
3993 CALL wrf_dm_bcast_double(tcs_racs1,SIZE(tcs_racs1))
3994 CALL wrf_dm_bcast_double(tmr_racs1,SIZE(tmr_racs1))
3995 CALL wrf_dm_bcast_double(tcs_racs2,SIZE(tcs_racs2))
3996 CALL wrf_dm_bcast_double(tmr_racs2,SIZE(tmr_racs2))
3997 CALL wrf_dm_bcast_double(tcr_sacr1,SIZE(tcr_sacr1))
3998 CALL wrf_dm_bcast_double(tms_sacr1,SIZE(tms_sacr1))
3999 CALL wrf_dm_bcast_double(tcr_sacr2,SIZE(tcr_sacr2))
4000 CALL wrf_dm_bcast_double(tms_sacr2,SIZE(tms_sacr2))
4001 CALL wrf_dm_bcast_double(tnr_racs1,SIZE(tnr_racs1))
4002 CALL wrf_dm_bcast_double(tnr_racs2,SIZE(tnr_racs2))
4003 CALL wrf_dm_bcast_double(tnr_sacr1,SIZE(tnr_sacr1))
4004 CALL wrf_dm_bcast_double(tnr_sacr2,SIZE(tnr_sacr2))
4007 CALL wrf_message("ThompMP: computing qr_acr_qs")
4009 ! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
4010 vr(n2) = -0.1021 + 4.932E3*Dr(n2) - 0.9551E6*Dr(n2)*Dr(n2) &
4011 + 0.07934E9*Dr(n2)*Dr(n2)*Dr(n2) &
4012 - 0.002362E12*Dr(n2)*Dr(n2)*Dr(n2)*Dr(n2)
4013 D1(n2) = (vr(n2)/av_s)**(1./bv_s)
4016 vs(n) = 1.5*av_s*Ds(n)**bv_s * DEXP(-fv_s*Ds(n))
4019 !..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for
4020 !.. fortran indices. J. Michalakes, 2009Oct30.
4022 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
4023 CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e )
4026 km_e = ntb_r*ntb_r1 - 1
4031 k = mod( km , ntb_r1 ) + 1
4033 lam_exp = (N0r_exp(k)*am_r*crg(1)/r_r(m))**ore1
4034 lamr = lam_exp * (crg(3)*org2*org1)**obmr
4035 N0_r = N0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2)
4037 N_r(n2) = N0_r*Dr(n2)**mu_r * DEXP(-lamr*Dr(n2))*dtr(n2)
4043 !..From the bm_s moment, compute plus one moment. If we are not
4044 !.. using bm_s=2, then we must transform to the pure 2nd moment
4045 !.. (variable called "second") and then to the bm_s+1 moment.
4047 M2 = r_s(i)*oams *1.0d0
4048 if (bm_s.gt.2.0-1.E-3 .and. bm_s.lt.2.0+1.E-3) then
4049 loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*bm_s &
4050 + sa(4)*Tc(j)*bm_s + sa(5)*Tc(j)*Tc(j) &
4051 + sa(6)*bm_s*bm_s + sa(7)*Tc(j)*Tc(j)*bm_s &
4052 + sa(8)*Tc(j)*bm_s*bm_s + sa(9)*Tc(j)*Tc(j)*Tc(j) &
4053 + sa(10)*bm_s*bm_s*bm_s
4055 b_ = sb(1) + sb(2)*Tc(j) + sb(3)*bm_s &
4056 + sb(4)*Tc(j)*bm_s + sb(5)*Tc(j)*Tc(j) &
4057 + sb(6)*bm_s*bm_s + sb(7)*Tc(j)*Tc(j)*bm_s &
4058 + sb(8)*Tc(j)*bm_s*bm_s + sb(9)*Tc(j)*Tc(j)*Tc(j) &
4059 + sb(10)*bm_s*bm_s*bm_s
4060 second = (M2/a_)**(1./b_)
4065 loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*cse(1) &
4066 + sa(4)*Tc(j)*cse(1) + sa(5)*Tc(j)*Tc(j) &
4067 + sa(6)*cse(1)*cse(1) + sa(7)*Tc(j)*Tc(j)*cse(1) &
4068 + sa(8)*Tc(j)*cse(1)*cse(1) + sa(9)*Tc(j)*Tc(j)*Tc(j) &
4069 + sa(10)*cse(1)*cse(1)*cse(1)
4071 b_ = sb(1)+sb(2)*Tc(j)+sb(3)*cse(1) + sb(4)*Tc(j)*cse(1) &
4072 + sb(5)*Tc(j)*Tc(j) + sb(6)*cse(1)*cse(1) &
4073 + sb(7)*Tc(j)*Tc(j)*cse(1) + sb(8)*Tc(j)*cse(1)*cse(1) &
4074 + sb(9)*Tc(j)*Tc(j)*Tc(j)+sb(10)*cse(1)*cse(1)*cse(1)
4075 M3 = a_ * second**b_
4078 Mrat = M2*(M2*oM3)*(M2*oM3)*(M2*oM3)
4080 slam1 = M2 * oM3 * Lam0
4081 slam2 = M2 * oM3 * Lam1
4084 N_s(n) = Mrat*(Kap0*DEXP(-slam1*Ds(n)) &
4085 + Kap1*M0*Ds(n)**mu_s * DEXP(-slam2*Ds(n)))*dts(n)
4101 massr = am_r * Dr(n2)**bm_r
4103 masss = am_s * Ds(n)**bm_s
4105 dvs = 0.5d0*((vr(n2) - vs(n)) + DABS(vr(n2)-vs(n)))
4106 dvr = 0.5d0*((vs(n) - vr(n2)) + DABS(vs(n)-vr(n2)))
4108 if (massr .gt. 1.5*masss) then
4109 t1 = t1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4110 *dvs*masss * N_s(n)* N_r(n2)
4111 z1 = z1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4112 *dvs*massr * N_s(n)* N_r(n2)
4113 y1 = y1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4114 *dvs * N_s(n)* N_r(n2)
4116 t3 = t3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4117 *dvs*masss * N_s(n)* N_r(n2)
4118 z3 = z3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4119 *dvs*massr * N_s(n)* N_r(n2)
4120 y3 = y3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4121 *dvs * N_s(n)* N_r(n2)
4124 if (massr .gt. 1.5*masss) then
4125 t2 = t2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4126 *dvr*massr * N_s(n)* N_r(n2)
4127 y2 = y2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4128 *dvr * N_s(n)* N_r(n2)
4129 z2 = z2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4130 *dvr*masss * N_s(n)* N_r(n2)
4132 t4 = t4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4133 *dvr*massr * N_s(n)* N_r(n2)
4134 y4 = y4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4135 *dvr * N_s(n)* N_r(n2)
4136 z4 = z4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) &
4137 *dvr*masss * N_s(n)* N_r(n2)
4142 tcs_racs1(i,j,k,m) = t1
4143 tmr_racs1(i,j,k,m) = DMIN1(z1, r_r(m)*1.0d0)
4144 tcs_racs2(i,j,k,m) = t3
4145 tmr_racs2(i,j,k,m) = z3
4146 tcr_sacr1(i,j,k,m) = t2
4147 tms_sacr1(i,j,k,m) = z2
4148 tcr_sacr2(i,j,k,m) = t4
4149 tms_sacr2(i,j,k,m) = z4
4150 tnr_racs1(i,j,k,m) = y1
4151 tnr_racs2(i,j,k,m) = y3
4152 tnr_sacr1(i,j,k,m) = y2
4153 tnr_sacr2(i,j,k,m) = y4
4158 !..Note wrf_dm_gatherv expects zero-based km_s, km_e (J. Michalakes, 2009Oct30).
4160 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
4161 CALL wrf_dm_gatherv(tcs_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4162 CALL wrf_dm_gatherv(tmr_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4163 CALL wrf_dm_gatherv(tcs_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4164 CALL wrf_dm_gatherv(tmr_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4165 CALL wrf_dm_gatherv(tcr_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4166 CALL wrf_dm_gatherv(tms_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4167 CALL wrf_dm_gatherv(tcr_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4168 CALL wrf_dm_gatherv(tms_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4169 CALL wrf_dm_gatherv(tnr_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4170 CALL wrf_dm_gatherv(tnr_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4171 CALL wrf_dm_gatherv(tnr_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4172 CALL wrf_dm_gatherv(tnr_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE)
4175 IF ( write_thompson_tables .AND. wrf_dm_on_monitor() ) THEN
4176 CALL wrf_message("Writing qr_acr_qsV2.dat in Thompson MP init")
4177 OPEN(63,file="qr_acr_qsV2.dat",form="unformatted",err=9234)
4178 WRITE(63,err=9234)tcs_racs1
4179 WRITE(63,err=9234)tmr_racs1
4180 WRITE(63,err=9234)tcs_racs2
4181 WRITE(63,err=9234)tmr_racs2
4182 WRITE(63,err=9234)tcr_sacr1
4183 WRITE(63,err=9234)tms_sacr1
4184 WRITE(63,err=9234)tcr_sacr2
4185 WRITE(63,err=9234)tms_sacr2
4186 WRITE(63,err=9234)tnr_racs1
4187 WRITE(63,err=9234)tnr_racs2
4188 WRITE(63,err=9234)tnr_sacr1
4189 WRITE(63,err=9234)tnr_sacr2
4191 RETURN ! ----- RETURN
4193 CALL wrf_error_fatal("Error writing qr_acr_qsV2.dat")
4197 end subroutine qr_acr_qs
4198 !+---+-----------------------------------------------------------------+
4200 !+---+-----------------------------------------------------------------+
4201 !..This is a literal adaptation of Bigg (1954) probability of drops of
4202 !..a particular volume freezing. Given this probability, simply freeze
4203 !..the proportion of drops summing their masses.
4204 !+---+-----------------------------------------------------------------+
4206 subroutine freezeH2O
4211 INTEGER:: i, j, k, m, n, n2
4212 INTEGER:: km, km_s, km_e
4213 DOUBLE PRECISION:: N_r, N_c
4214 DOUBLE PRECISION, DIMENSION(nbr):: massr
4215 DOUBLE PRECISION, DIMENSION(nbc):: massc
4216 DOUBLE PRECISION:: sum1, sum2, sumn1, sumn2, &
4217 prob, vol, Texp, orho_w, &
4218 lam_exp, lamr, N0_r, lamc, N0_c, y
4221 LOGICAL force_read_thompson, write_thompson_tables
4222 LOGICAL lexist,lopen
4224 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
4227 CALL nl_get_force_read_thompson(1,force_read_thompson)
4228 CALL nl_get_write_thompson_tables(1,write_thompson_tables)
4231 IF ( wrf_dm_on_monitor() ) THEN
4232 INQUIRE(FILE="freezeH2O.dat",EXIST=lexist)
4234 CALL wrf_message("ThompMP: read freezeH2O.dat instead of computing")
4235 OPEN(63,file="freezeH2O.dat",form="unformatted",err=1234)
4236 READ(63,err=1234)tpi_qrfz
4237 READ(63,err=1234)tni_qrfz
4238 READ(63,err=1234)tpg_qrfz
4239 READ(63,err=1234)tnr_qrfz
4240 READ(63,err=1234)tpi_qcfz
4241 READ(63,err=1234)tni_qcfz
4244 IF ( good .NE. 1 ) THEN
4245 INQUIRE(63,opened=lopen)
4247 IF( force_read_thompson ) THEN
4248 CALL wrf_error_fatal("Error reading freezeH2O.dat. Aborting because force_read_thompson is .true.")
4252 IF( force_read_thompson ) THEN
4253 CALL wrf_error_fatal("Error opening freezeH2O.dat. Aborting because force_read_thompson is .true.")
4257 INQUIRE(63,opened=lopen)
4263 IF( force_read_thompson ) THEN
4264 CALL wrf_error_fatal("Non-existent freezeH2O.dat. Aborting because force_read_thompson is .true.")
4269 #if defined(DM_PARALLEL) && !defined(STUBMPI)
4270 CALL wrf_dm_bcast_integer(good,1)
4273 IF ( good .EQ. 1 ) THEN
4274 #if defined(DM_PARALLEL) && !defined(STUBMPI)
4275 CALL wrf_dm_bcast_double(tpi_qrfz,SIZE(tpi_qrfz))
4276 CALL wrf_dm_bcast_double(tni_qrfz,SIZE(tni_qrfz))
4277 CALL wrf_dm_bcast_double(tpg_qrfz,SIZE(tpg_qrfz))
4278 CALL wrf_dm_bcast_double(tnr_qrfz,SIZE(tnr_qrfz))
4279 CALL wrf_dm_bcast_double(tpi_qcfz,SIZE(tpi_qcfz))
4280 CALL wrf_dm_bcast_double(tni_qcfz,SIZE(tni_qcfz))
4283 CALL wrf_message("ThompMP: computing freezeH2O")
4288 massr(n2) = am_r*Dr(n2)**bm_r
4291 massc(n) = am_r*Dc(n)**bm_r
4294 ! Need to split loops between MPI processes to speedup
4295 ! (2017Jul26, Jason Do)
4296 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
4297 CALL wrf_dm_decomp1d ( ntb_IN*45, km_s, km_e )
4300 km_e = ntb_IN*45 - 1
4303 !..Freeze water (smallest drops become cloud ice, otherwise graupel).
4306 k = mod( km , 45 ) + 1
4307 T_adjust = MAX(-3.0, MIN(3.0 - ALOG10(Nt_IN(m)), 3.0))
4308 ! print*, ' Freezing water for temp = ', -k
4309 Texp = DEXP( DFLOAT(k) - T_adjust*1.0D0 ) - 1.0D0
4311 !$OMP PARALLEL DO SCHEDULE(dynamic) &
4312 !$OMP PRIVATE(j,i,lam_exp,lamr,N0_r,sum1,sum2,sumn1,sumn2,n2,N_r,vol,prob)
4316 lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(i))**ore1
4317 lamr = lam_exp * (crg(3)*org2*org1)**obmr
4318 N0_r = N0r_exp(j)/(crg(2)*lam_exp) * lamr**cre(2)
4324 N_r = N0_r*Dr(n2)**mu_r*DEXP(-lamr*Dr(n2))*dtr(n2)
4325 vol = massr(n2)*orho_w
4326 prob = MAX(0.0D0, 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp))
4327 if (massr(n2) .lt. xm0g) then
4328 sumn1 = sumn1 + prob*N_r
4329 sum1 = sum1 + prob*N_r*massr(n2)
4331 sumn2 = sumn2 + prob*N_r
4332 sum2 = sum2 + prob*N_r*massr(n2)
4334 if ((sum1+sum2).ge.r_r(i)) EXIT
4336 tpi_qrfz(i,j,k,m) = sum1
4337 tni_qrfz(i,j,k,m) = sumn1
4338 tpg_qrfz(i,j,k,m) = sum2
4339 tnr_qrfz(i,j,k,m) = sumn2
4343 !$OMP END PARALLEL DO
4345 !$OMP PARALLEL DO SCHEDULE(dynamic) &
4346 !$OMP PRIVATE(j,i,nu_c,lamc,N0_c,sum1,sumn2,vol,prob,N_c)
4349 nu_c = MIN(15, NINT(1000.E6/t_Nc(j)) + 2)
4351 lamc = (t_Nc(j)*am_r* ccg(2,nu_c) * ocg1(nu_c) / r_c(i))**obmr
4352 N0_c = t_Nc(j)*ocg1(nu_c) * lamc**cce(1,nu_c)
4356 vol = massc(n)*orho_w
4357 prob = MAX(0.0D0, 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp))
4358 N_c = N0_c*Dc(n)**nu_c*EXP(-lamc*Dc(n))*dtc(n)
4359 sumn2 = MIN(t_Nc(j), sumn2 + prob*N_c)
4360 sum1 = sum1 + prob*N_c*massc(n)
4361 if (sum1 .ge. r_c(i)) EXIT
4363 tpi_qcfz(i,j,k,m) = sum1
4364 tni_qcfz(i,j,k,m) = sumn2
4368 !$OMP END PARALLEL DO
4372 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
4373 CALL wrf_dm_gatherv(tpi_qrfz, ntb_r*ntb_r1, km_s, km_e, R8SIZE)
4374 CALL wrf_dm_gatherv(tni_qrfz, ntb_r*ntb_r1, km_s, km_e, R8SIZE)
4375 CALL wrf_dm_gatherv(tpg_qrfz, ntb_r*ntb_r1, km_s, km_e, R8SIZE)
4376 CALL wrf_dm_gatherv(tnr_qrfz, ntb_r*ntb_r1, km_s, km_e, R8SIZE)
4377 CALL wrf_dm_gatherv(tpi_qcfz, ntb_c*nbc, km_s, km_e, R8SIZE)
4378 CALL wrf_dm_gatherv(tni_qcfz, ntb_c*nbc, km_s, km_e, R8SIZE)
4381 IF ( write_thompson_tables .AND. wrf_dm_on_monitor() ) THEN
4382 CALL wrf_message("Writing freezeH2O.dat in Thompson MP init")
4383 OPEN(63,file="freezeH2O.dat",form="unformatted",err=9234)
4384 WRITE(63,err=9234)tpi_qrfz
4385 WRITE(63,err=9234)tni_qrfz
4386 WRITE(63,err=9234)tpg_qrfz
4387 WRITE(63,err=9234)tnr_qrfz
4388 WRITE(63,err=9234)tpi_qcfz
4389 WRITE(63,err=9234)tni_qcfz
4391 RETURN ! ----- RETURN
4393 CALL wrf_error_fatal("Error writing freezeH2O.dat")
4397 end subroutine freezeH2O
4399 !+---+-----------------------------------------------------------------+
4401 !+---+-----------------------------------------------------------------+
4402 !..Cloud ice converting to snow since portion greater than min snow
4403 !.. size. Given cloud ice content (kg/m**3), number concentration
4404 !.. (#/m**3) and gamma shape parameter, mu_i, break the distrib into
4405 !.. bins and figure out the mass/number of ice with sizes larger than
4406 !.. D0s. Also, compute incomplete gamma function for the integration
4407 !.. of ice depositional growth from diameter=0 to D0s. Amount of
4408 !.. ice depositional growth is this portion of distrib while larger
4409 !.. diameters contribute to snow growth (as in Harrington et al. 1995).
4410 !+---+-----------------------------------------------------------------+
4412 subroutine qi_aut_qs
4418 DOUBLE PRECISION, DIMENSION(nbi):: N_i
4419 DOUBLE PRECISION:: N0_i, lami, Di_mean, t1, t2
4426 lami = (am_i*cig(2)*oig1*Nt_i(j)/r_i(i))**obmi
4427 Di_mean = (bm_i + mu_i + 1.) / lami
4428 N0_i = Nt_i(j)*oig1 * lami**cie(1)
4431 if (SNGL(Di_mean) .gt. 5.*D0s) then
4434 tpi_ide(i,j) = 0.0D0
4435 elseif (SNGL(Di_mean) .lt. D0i) then
4438 tpi_ide(i,j) = 1.0D0
4440 xlimit_intg = lami*D0s
4441 tpi_ide(i,j) = GAMMP(mu_i+2.0, xlimit_intg) * 1.0D0
4443 N_i(n2) = N0_i*Di(n2)**mu_i * DEXP(-lami*Di(n2))*dti(n2)
4444 if (Di(n2).ge.D0s) then
4445 t1 = t1 + N_i(n2) * am_i*Di(n2)**bm_i
4455 end subroutine qi_aut_qs
4457 !+---+-----------------------------------------------------------------+
4458 !..Variable collision efficiency for rain collecting cloud water using
4459 !.. method of Beard and Grover, 1974 if a/A less than 0.25; otherwise
4460 !.. uses polynomials to get close match of Pruppacher & Klett Fig 14-9.
4461 !+---+-----------------------------------------------------------------+
4463 subroutine table_Efrw
4468 DOUBLE PRECISION:: vtr, stokes, reynolds, Ef_rw
4469 DOUBLE PRECISION:: p, yc0, F, G, H, z, K0, X
4476 if (Dr(i).lt.50.E-6 .or. Dc(j).lt.3.E-6) then
4478 elseif (p.gt.0.25) then
4480 if (Dr(i) .lt. 75.e-6) then
4481 Ef_rw = 0.026794*X - 0.20604
4482 elseif (Dr(i) .lt. 125.e-6) then
4483 Ef_rw = -0.00066842*X*X + 0.061542*X - 0.37089
4484 elseif (Dr(i) .lt. 175.e-6) then
4485 Ef_rw = 4.091e-06*X*X*X*X - 0.00030908*X*X*X &
4486 + 0.0066237*X*X - 0.0013687*X - 0.073022
4487 elseif (Dr(i) .lt. 250.e-6) then
4488 Ef_rw = 9.6719e-5*X*X*X - 0.0068901*X*X + 0.17305*X &
4490 elseif (Dr(i) .lt. 350.e-6) then
4491 Ef_rw = 9.0488e-5*X*X*X - 0.006585*X*X + 0.16606*X &
4494 Ef_rw = 0.00010721*X*X*X - 0.0072962*X*X + 0.1704*X &
4498 vtr = -0.1021 + 4.932E3*Dr(i) - 0.9551E6*Dr(i)*Dr(i) &
4499 + 0.07934E9*Dr(i)*Dr(i)*Dr(i) &
4500 - 0.002362E12*Dr(i)*Dr(i)*Dr(i)*Dr(i)
4501 stokes = Dc(j)*Dc(j)*vtr*rho_w/(9.*1.718E-5*Dr(i))
4502 reynolds = 9.*stokes/(p*p*rho_w)
4505 G = -0.1007D0 - 0.358D0*F + 0.0261D0*F*F
4507 z = DLOG(stokes/(K0+1.D-15))
4508 H = 0.1465D0 + 1.302D0*z - 0.607D0*z*z + 0.293D0*z*z*z
4509 yc0 = 2.0D0/PI * ATAN(H)
4510 Ef_rw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p))
4514 t_Efrw(i,j) = MAX(0.0, MIN(SNGL(Ef_rw), 0.95))
4519 end subroutine table_Efrw
4521 !+---+-----------------------------------------------------------------+
4522 !..Variable collision efficiency for snow collecting cloud water using
4523 !.. method of Wang and Ji, 2000 except equate melted snow diameter to
4524 !.. their "effective collision cross-section."
4525 !+---+-----------------------------------------------------------------+
4527 subroutine table_Efsw
4532 DOUBLE PRECISION:: Ds_m, vts, vtc, stokes, reynolds, Ef_sw
4533 DOUBLE PRECISION:: p, yc0, F, G, H, z, K0
4537 vtc = 1.19D4 * (1.0D4*Dc(j)*Dc(j)*0.25D0)
4539 vts = av_s*Ds(i)**bv_s * DEXP(-fv_s*Ds(i)) - vtc
4540 Ds_m = (am_s*Ds(i)**bm_s / am_r)**obmr
4542 if (p.gt.0.25 .or. Ds(i).lt.D0s .or. Dc(j).lt.6.E-6 &
4543 .or. vts.lt.1.E-3) then
4546 stokes = Dc(j)*Dc(j)*vts*rho_w/(9.*1.718E-5*Ds_m)
4547 reynolds = 9.*stokes/(p*p*rho_w)
4550 G = -0.1007D0 - 0.358D0*F + 0.0261D0*F*F
4552 z = DLOG(stokes/(K0+1.D-15))
4553 H = 0.1465D0 + 1.302D0*z - 0.607D0*z*z + 0.293D0*z*z*z
4554 yc0 = 2.0D0/PI * ATAN(H)
4555 Ef_sw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p))
4557 t_Efsw(i,j) = MAX(0.0, MIN(SNGL(Ef_sw), 0.95))
4563 end subroutine table_Efsw
4565 !+---+-----------------------------------------------------------------+
4566 !..Function to compute collision efficiency of collector species (rain,
4567 !.. snow, graupel) of aerosols. Follows Wang et al, 2010, ACP, which
4568 !.. follows Slinn (1983).
4569 !+---+-----------------------------------------------------------------+
4571 real function Eff_aero(D, Da, visc,rhoa,Temp,species)
4574 real:: D, Da, visc, rhoa, Temp
4575 character(LEN=1):: species
4576 real:: aval, Cc, diff, Re, Sc, St, St2, vt, Eff
4577 real, parameter:: boltzman = 1.3806503E-23
4578 real, parameter:: meanPath = 0.0256E-6
4581 if (species .eq. 'r') then
4582 vt = -0.1021 + 4.932E3*D - 0.9551E6*D*D &
4583 + 0.07934E9*D*D*D - 0.002362E12*D*D*D*D
4584 elseif (species .eq. 's') then
4586 elseif (species .eq. 'g') then
4590 Cc = 1. + 2.*meanPath/Da *(1.257+0.4*exp(-0.55*Da/meanPath))
4591 diff = boltzman*Temp*Cc/(3.*PI*visc*Da)
4593 Re = 0.5*rhoa*D*vt/visc
4594 Sc = visc/(rhoa*diff)
4596 St = Da*Da*vt*1000./(9.*visc*D)
4597 aval = 1.+LOG(1.+Re)
4598 St2 = (1.2 + 1./12.*aval)/(1.+aval)
4600 Eff = 4./(Re*Sc) * (1. + 0.4*SQRT(Re)*Sc**0.3333 &
4601 + 0.16*SQRT(Re)*SQRT(Sc)) &
4602 + 4.*Da/D * (0.02 + Da/D*(1.+2.*SQRT(Re)))
4604 if (St.gt.St2) Eff = Eff + ( (St-St2)/(St-St2+0.666667))**1.5
4605 Eff_aero = MAX(1.E-5, MIN(Eff, 1.0))
4607 end function Eff_aero
4610 !+---+-----------------------------------------------------------------+
4611 !..Integrate rain size distribution from zero to D-star to compute the
4612 !.. number of drops smaller than D-star that evaporate in a single
4613 !.. timestep. Drops larger than D-star dont evaporate entirely so do
4614 !.. not affect number concentration.
4615 !+---+-----------------------------------------------------------------+
4617 subroutine table_dropEvap
4622 INTEGER:: i, j, k, n
4623 DOUBLE PRECISION, DIMENSION(nbc):: N_c, massc
4624 DOUBLE PRECISION:: summ, summ2, lamc, N0_c
4626 ! DOUBLE PRECISION:: Nt_r, N0, lam_exp, lam
4627 ! REAL:: xlimit_intg
4630 massc(n) = am_r*Dc(n)**bm_r
4634 nu_c = MIN(15, NINT(1000.E6/t_Nc(k)) + 2)
4636 lamc = (t_Nc(k)*am_r* ccg(2,nu_c)*ocg1(nu_c) / r_c(j))**obmr
4637 N0_c = t_Nc(k)*ocg1(nu_c) * lamc**cce(1,nu_c)
4639 !-GT tnc_wev(i,j,k) = GAMMP(nu_c+1., SNGL(Dc(i)*lamc))*t_Nc(k)
4640 N_c(i) = N0_c* Dc(i)**nu_c*EXP(-lamc*Dc(i))*dtc(i)
4641 ! if(j.eq.18 .and. k.eq.50) print*, ' N_c = ', N_c(i)
4645 summ = summ + massc(n)*N_c(n)
4646 summ2 = summ2 + N_c(n)
4648 ! if(j.eq.18 .and. k.eq.50) print*, ' DEBUG-TABLE: ', r_c(j), t_Nc(k), summ2, summ
4649 tpc_wev(i,j,k) = summ
4650 tnc_wev(i,j,k) = summ2
4656 !..To do the same thing for rain.
4660 ! lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(k))**ore1
4661 ! lam = lam_exp * (crg(3)*org2*org1)**obmr
4662 ! N0 = N0r_exp(j)/(crg(2)*lam_exp) * lam**cre(2)
4663 ! Nt_r = N0 * crg(2) / lam**cre(2)
4665 ! xlimit_intg = lam*Dr(i)
4666 ! tnr_rev(i,j,k) = GAMMP(mu_r+1.0, xlimit_intg) * Nt_r
4671 ! TO APPLY TABLE ABOVE
4672 !..Rain lookup table indexes.
4673 ! Dr_star = DSQRT(-2.D0*DT * t1_evap/(2.*PI) &
4674 ! * 0.78*4.*diffu(k)*xsat*rvs/rho_w)
4675 ! idx_d = NINT(1.0 + FLOAT(nbr) * DLOG(Dr_star/D0r) &
4676 ! / DLOG(Dr(nbr)/D0r))
4677 ! idx_d = MAX(1, MIN(idx_d, nbr))
4679 ! nir = NINT(ALOG10(rr(k)))
4680 ! do nn = nir-1, nir+1
4682 ! if ( (rr(k)/10.**nn).ge.1.0 .and. &
4683 ! (rr(k)/10.**nn).lt.10.0) goto 154
4686 ! idx_r = INT(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2)
4687 ! idx_r = MAX(1, MIN(idx_r, ntb_r))
4689 ! lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
4690 ! lam_exp = lamr * (crg(3)*org2*org1)**bm_r
4691 ! N0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
4692 ! nir = NINT(DLOG10(N0_exp))
4693 ! do nn = nir-1, nir+1
4695 ! if ( (N0_exp/10.**nn).ge.1.0 .and. &
4696 ! (N0_exp/10.**nn).lt.10.0) goto 155
4699 ! idx_r1 = INT(N0_exp/10.**n) + 10*(n-nir3) - (n-nir3)
4700 ! idx_r1 = MAX(1, MIN(idx_r1, ntb_r1))
4702 ! pnr_rev(k) = MIN(nr(k)*odts, SNGL(tnr_rev(idx_d,idx_r1,idx_r) & ! RAIN2M
4705 end subroutine table_dropEvap
4708 !+---+-----------------------------------------------------------------+
4709 !..Fill the table of CCN activation data created from parcel model run
4710 !.. by Trude Eidhammer with inputs of aerosol number concentration,
4711 !.. vertical velocity, temperature, lognormal mean aerosol radius, and
4712 !.. hygroscopicity, kappa. The data are read from external file and
4713 !.. contain activated fraction of CCN for given conditions.
4714 !+---+-----------------------------------------------------------------+
4716 subroutine table_ccnAct
4722 LOGICAL, EXTERNAL:: wrf_dm_on_monitor
4725 INTEGER:: iunit_mp_th1, i
4727 CHARACTER*64 errmess
4730 IF ( wrf_dm_on_monitor() ) THEN
4732 INQUIRE ( i , OPENED = opened )
4733 IF ( .NOT. opened ) THEN
4740 #if defined(DM_PARALLEL) && !defined(STUBMPI)
4741 CALL wrf_dm_bcast_bytes ( iunit_mp_th1 , IWORDSIZE )
4743 IF ( iunit_mp_th1 < 0 ) THEN
4744 CALL wrf_error_fatal ( 'module_mp_thompson: table_ccnAct: '// &
4745 'Can not find unused fortran unit to read in lookup table.')
4748 IF ( wrf_dm_on_monitor() ) THEN
4749 WRITE(errmess, '(A,I2)') 'module_mp_thompson: opening CCN_ACTIVATE.BIN on unit ',iunit_mp_th1
4750 CALL wrf_debug(150, errmess)
4751 OPEN(iunit_mp_th1,FILE='CCN_ACTIVATE.BIN', &
4752 FORM='UNFORMATTED',STATUS='OLD',ERR=9009)
4755 #define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes(A, size(A)*R4SIZE)
4757 IF ( wrf_dm_on_monitor() ) READ(iunit_mp_th1,ERR=9010) tnccn_act
4758 #if defined(DM_PARALLEL) && !defined(STUBMPI)
4759 DM_BCAST_MACRO(tnccn_act)
4765 WRITE( errmess , '(A,I2)' ) 'module_mp_thompson: error opening CCN_ACTIVATE.BIN on unit ',iunit_mp_th1
4766 CALL wrf_error_fatal(errmess)
4769 WRITE( errmess , '(A,I2)' ) 'module_mp_thompson: error reading CCN_ACTIVATE.BIN on unit ',iunit_mp_th1
4770 CALL wrf_error_fatal(errmess)
4772 end subroutine table_ccnAct
4774 !+---+-----------------------------------------------------------------+
4775 !..Retrieve fraction of CCN that gets activated given the model temp,
4776 !.. vertical velocity, and available CCN concentration. The lookup
4777 !.. table (read from external file) has CCN concentration varying the
4778 !.. quickest, then updraft, then temperature, then mean aerosol radius,
4779 !.. and finally hygroscopicity, kappa.
4780 !.. TO_DO ITEM: For radiation cooling producing fog, in which case the
4781 !.. updraft velocity could easily be negative, we could use the temp
4782 !.. and its tendency to diagnose a pretend postive updraft velocity.
4783 !+---+-----------------------------------------------------------------+
4784 real function activ_ncloud(Tt, Ww, NCCN)
4787 REAL, INTENT(IN):: Tt, Ww, NCCN
4788 REAL:: n_local, w_local
4789 INTEGER:: i, j, k, l, m, n
4790 REAL:: A, B, C, D, t, u, x1, x2, y1, y2, nx, wy, fraction
4793 ! ta_Na = (/10.0, 31.6, 100.0, 316.0, 1000.0, 3160.0, 10000.0/) ntb_arc
4794 ! ta_Ww = (/0.01, 0.0316, 0.1, 0.316, 1.0, 3.16, 10.0, 31.6, 100.0/) ntb_arw
4795 ! ta_Tk = (/243.15, 253.15, 263.15, 273.15, 283.15, 293.15, 303.15/) ntb_art
4796 ! ta_Ra = (/0.01, 0.02, 0.04, 0.08, 0.16/) ntb_arr
4797 ! ta_Ka = (/0.2, 0.4, 0.6, 0.8/) ntb_ark
4799 n_local = NCCN * 1.E-6
4802 if (n_local .ge. ta_Na(ntb_arc)) then
4803 n_local = ta_Na(ntb_arc) - 1.0
4804 elseif (n_local .le. ta_Na(1)) then
4805 n_local = ta_Na(1) + 1.0
4808 if (n_local.ge.ta_Na(n-1) .and. n_local.lt.ta_Na(n)) goto 8003
4812 x1 = LOG(ta_Na(i-1))
4815 if (w_local .ge. ta_Ww(ntb_arw)) then
4816 w_local = ta_Ww(ntb_arw) - 1.0
4817 elseif (w_local .le. ta_Ww(1)) then
4818 w_local = ta_Ww(1) + 0.001
4821 if (w_local.ge.ta_Ww(n-1) .and. w_local.lt.ta_Ww(n)) goto 8005
4825 y1 = LOG(ta_Ww(j-1))
4828 k = MAX(1, MIN( NINT( (Tt - ta_Tk(1))*0.1) + 1, ntb_art))
4830 !..The next two values are indexes of mean aerosol radius and
4831 !.. hygroscopicity. Currently these are constant but a future version
4832 !.. should implement other variables to allow more freedom such as
4833 !.. at least simple separation of tiny size sulfates from larger
4838 A = tnccn_act(i-1,j-1,k,l,m)
4839 B = tnccn_act(i,j-1,k,l,m)
4840 C = tnccn_act(i,j,k,l,m)
4841 D = tnccn_act(i-1,j,k,l,m)
4848 ! t = (n_local-ta(Na(i-1))/(ta_Na(i)-ta_Na(i-1))
4849 ! u = (w_local-ta_Ww(j-1))/(ta_Ww(j)-ta_Ww(j-1))
4851 fraction = (1.0-t)*(1.0-u)*A + t*(1.0-u)*B + t*u*C + (1.0-t)*u*D
4853 ! if (NCCN*fraction .gt. 0.75*Nt_c_max) then
4854 ! write(*,*) ' DEBUG-GT ', n_local, w_local, Tt, i, j, k
4857 activ_ncloud = NCCN*fraction
4859 end function activ_ncloud
4861 !+---+-----------------------------------------------------------------+
4862 !+---+-----------------------------------------------------------------+
4863 SUBROUTINE GCF(GAMMCF,A,X,GLN)
4864 ! --- RETURNS THE INCOMPLETE GAMMA FUNCTION Q(A,X) EVALUATED BY ITS
4865 ! --- CONTINUED FRACTION REPRESENTATION AS GAMMCF. ALSO RETURNS
4866 ! --- LN(GAMMA(A)) AS GLN. THE CONTINUED FRACTION IS EVALUATED BY
4867 ! --- A MODIFIED LENTZ METHOD.
4870 INTEGER, PARAMETER:: ITMAX=100
4871 REAL, PARAMETER:: gEPS=3.E-7
4872 REAL, PARAMETER:: FPMIN=1.E-30
4873 REAL, INTENT(IN):: A, X
4876 REAL:: AN,B,C,D,DEL,H
4886 IF(ABS(D).LT.FPMIN)D=FPMIN
4888 IF(ABS(C).LT.FPMIN)C=FPMIN
4892 IF(ABS(DEL-1.).LT.gEPS)GOTO 1
4894 PRINT *, 'A TOO LARGE, ITMAX TOO SMALL IN GCF'
4895 1 GAMMCF=EXP(-X+A*LOG(X)-GLN)*H
4897 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
4898 !+---+-----------------------------------------------------------------+
4899 SUBROUTINE GSER(GAMSER,A,X,GLN)
4900 ! --- RETURNS THE INCOMPLETE GAMMA FUNCTION P(A,X) EVALUATED BY ITS
4901 ! --- ITS SERIES REPRESENTATION AS GAMSER. ALSO RETURNS LN(GAMMA(A))
4905 INTEGER, PARAMETER:: ITMAX=100
4906 REAL, PARAMETER:: gEPS=3.E-7
4907 REAL, INTENT(IN):: A, X
4913 IF(X.LT.0.) PRINT *, 'X < 0 IN GSER'
4924 IF(ABS(DEL).LT.ABS(SUM)*gEPS)GOTO 1
4926 PRINT *,'A TOO LARGE, ITMAX TOO SMALL IN GSER'
4927 1 GAMSER=SUM*EXP(-X+A*LOG(X)-GLN)
4929 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
4930 !+---+-----------------------------------------------------------------+
4931 REAL FUNCTION GAMMLN(XX)
4932 ! --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 0.
4934 REAL, INTENT(IN):: XX
4935 DOUBLE PRECISION, PARAMETER:: STP = 2.5066282746310005D0
4936 DOUBLE PRECISION, DIMENSION(6), PARAMETER:: &
4937 COF = (/76.18009172947146D0, -86.50532032941677D0, &
4938 24.01409824083091D0, -1.231739572450155D0, &
4939 .1208650973866179D-2, -.5395239384953D-5/)
4940 DOUBLE PRECISION:: SER,TMP,X,Y
4946 TMP=(X+0.5D0)*LOG(TMP)-TMP
4947 SER=1.000000000190015D0
4952 GAMMLN=TMP+LOG(STP*SER/X)
4954 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
4955 !+---+-----------------------------------------------------------------+
4956 REAL FUNCTION GAMMP(A,X)
4957 ! --- COMPUTES THE INCOMPLETE GAMMA FUNCTION P(A,X)
4958 ! --- SEE ABRAMOWITZ AND STEGUN 6.5.1
4961 REAL, INTENT(IN):: A,X
4962 REAL:: GAMMCF,GAMSER,GLN
4964 IF((X.LT.0.) .OR. (A.LE.0.)) THEN
4965 PRINT *, 'BAD ARGUMENTS IN GAMMP'
4967 ELSEIF(X.LT.A+1.)THEN
4968 CALL GSER(GAMSER,A,X,GLN)
4971 CALL GCF(GAMMCF,A,X,GLN)
4975 ! (C) Copr. 1986-92 Numerical Recipes Software 2.02
4976 !+---+-----------------------------------------------------------------+
4977 REAL FUNCTION WGAMMA(y)
4980 REAL, INTENT(IN):: y
4982 WGAMMA = EXP(GAMMLN(y))
4985 !+---+-----------------------------------------------------------------+
4986 ! THIS FUNCTION CALCULATES THE LIQUID SATURATION VAPOR MIXING RATIO AS
4987 ! A FUNCTION OF TEMPERATURE AND PRESSURE
4989 REAL FUNCTION RSLF(P,T)
4992 REAL, INTENT(IN):: P, T
4994 REAL, PARAMETER:: C0= .611583699E03
4995 REAL, PARAMETER:: C1= .444606896E02
4996 REAL, PARAMETER:: C2= .143177157E01
4997 REAL, PARAMETER:: C3= .264224321E-1
4998 REAL, PARAMETER:: C4= .299291081E-3
4999 REAL, PARAMETER:: C5= .203154182E-5
5000 REAL, PARAMETER:: C6= .702620698E-8
5001 REAL, PARAMETER:: C7= .379534310E-11
5002 REAL, PARAMETER:: C8=-.321582393E-13
5004 X=MAX(-80.,T-273.16)
5006 ! ESL=612.2*EXP(17.67*X/(T-29.65))
5007 ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
5008 ESL=MIN(ESL, P*0.15) ! Even with P=1050mb and T=55C, the sat. vap. pres only contributes to ~15% of total pres.
5009 RSLF=.622*ESL/(P-ESL)
5012 ! ; Source: Murphy and Koop, Review of the vapour pressure of ice and
5013 ! supercooled water for atmospheric applications, Q. J. R.
5014 ! Meteorol. Soc (2005), 131, pp. 1539-1565.
5015 ! ESL = EXP(54.842763 - 6763.22 / T - 4.210 * ALOG(T) + 0.000367 * T
5016 ! + TANH(0.0415 * (T - 218.8)) * (53.878 - 1331.22
5017 ! / T - 9.44523 * ALOG(T) + 0.014025 * T))
5020 !+---+-----------------------------------------------------------------+
5021 ! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A
5022 ! FUNCTION OF TEMPERATURE AND PRESSURE
5024 REAL FUNCTION RSIF(P,T)
5027 REAL, INTENT(IN):: P, T
5029 REAL, PARAMETER:: C0= .609868993E03
5030 REAL, PARAMETER:: C1= .499320233E02
5031 REAL, PARAMETER:: C2= .184672631E01
5032 REAL, PARAMETER:: C3= .402737184E-1
5033 REAL, PARAMETER:: C4= .565392987E-3
5034 REAL, PARAMETER:: C5= .521693933E-5
5035 REAL, PARAMETER:: C6= .307839583E-7
5036 REAL, PARAMETER:: C7= .105785160E-9
5037 REAL, PARAMETER:: C8= .161444444E-12
5039 X=MAX(-80.,T-273.16)
5040 ESI=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
5041 ESI=MIN(ESI, P*0.15)
5042 RSIF=.622*ESI/(P-ESI)
5045 ! ; Source: Murphy and Koop, Review of the vapour pressure of ice and
5046 ! supercooled water for atmospheric applications, Q. J. R.
5047 ! Meteorol. Soc (2005), 131, pp. 1539-1565.
5048 ! ESI = EXP(9.550426 - 5723.265/T + 3.53068*ALOG(T) - 0.00728332*T)
5052 !+---+-----------------------------------------------------------------+
5053 real function iceDeMott(tempc, qv, qvs, qvsi, rho, nifa)
5056 REAL, INTENT(IN):: tempc, qv, qvs, qvsi, rho, nifa
5059 REAL:: satw, sati, siw, p_x, si0x, dtt, dsi, dsw, dab, fc, hx
5060 REAL:: ntilde, n_in, nmax, nhat, mux, xni, nifa_cc
5061 REAL, PARAMETER:: p_c1 = 1000.
5062 REAL, PARAMETER:: p_rho_c = 0.76
5063 REAL, PARAMETER:: p_alpha = 1.0
5064 REAL, PARAMETER:: p_gam = 2.
5065 REAL, PARAMETER:: delT = 5.
5066 REAL, PARAMETER:: T0x = -40.
5067 REAL, PARAMETER:: Sw0x = 0.97
5068 REAL, PARAMETER:: delSi = 0.1
5069 REAL, PARAMETER:: hdm = 0.15
5070 REAL, PARAMETER:: p_psi = 0.058707*p_gam/p_rho_c
5071 REAL, PARAMETER:: aap = 1.
5072 REAL, PARAMETER:: bbp = 0.
5073 REAL, PARAMETER:: y1p = -35.
5074 REAL, PARAMETER:: y2p = -25.
5075 REAL, PARAMETER:: rho_not0 = 101325./(287.05*273.15)
5083 ! p_x = -1.0261+(3.1656e-3*tempc)+(5.3938e-4*(tempc*tempc)) &
5084 ! + (8.2584e-6*(tempc*tempc*tempc))
5085 ! si0x = 1.+(10.**p_x)
5086 ! if (sati.ge.si0x .and. satw.lt.0.985) then
5087 ! dtt = delta_p (tempc, T0x, T0x+delT, 1., hdm)
5088 ! dsi = delta_p (sati, Si0x, Si0x+delSi, 0., 1.)
5089 ! dsw = delta_p (satw, Sw0x, 1., 0., 1.)
5091 ! hx = min(fc+((1.-fc)*dsw), 1.)
5092 ! ntilde = p_c1*p_gam*((exp(12.96*(sati-1.1)))**0.3) / p_rho_c
5093 ! if (tempc .le. y1p) then
5095 ! elseif (tempc .ge. y2p) then
5096 ! n_in = p_psi*p_c1*exp(12.96*(sati-1.)-0.639)
5098 ! if (tempc .le. -30.) then
5099 ! nmax = p_c1*p_gam*(exp(12.96*(siw-1.1)))**0.3/p_rho_c
5101 ! nmax = p_psi*p_c1*exp(12.96*(siw-1.)-0.639)
5103 ! ntilde = MIN(ntilde, nmax)
5104 ! nhat = MIN(p_psi*p_c1*exp(12.96*(sati-1.)-0.639), nmax)
5105 ! dab = delta_p (tempc, y1p, y2p, aap, bbp)
5106 ! n_in = MIN(nhat*(ntilde/nhat)**dab, nmax)
5108 ! mux = hx*p_alpha*n_in*rho
5109 ! xni = mux*((6700.*nifa)-200.)/((6700.*5.E5)-200.)
5110 ! elseif (satw.ge.0.985 .and. tempc.gt.HGFR-273.15) then
5111 nifa_cc = MAX(0.5, nifa*RHO_NOT0*1.E-6/rho)
5112 ! xni = 3.*nifa_cc**(1.25)*exp((0.46*(-tempc))-11.6) ! [DeMott, 2015]
5113 xni = (5.94e-5*(-tempc)**3.33) & ! [DeMott, 2010]
5114 * (nifa_cc**((-0.0264*(tempc))+0.0033))
5115 xni = xni*rho/RHO_NOT0 * 1000.
5118 iceDeMott = MAX(0., xni)
5120 end FUNCTION iceDeMott
5122 !+---+-----------------------------------------------------------------+
5123 !..Newer research since Koop et al (2001) suggests that the freezing
5124 !.. rate should be lower than original paper, so J_rate is reduced
5125 !.. by two orders of magnitude.
5127 real function iceKoop(temp, qv, qvs, naero, dt)
5130 REAL, INTENT(IN):: temp, qv, qvs, naero, DT
5131 REAL:: mu_diff, a_w_i, delta_aw, log_J_rate, J_rate, prob_h, satw
5136 mu_diff = 210368.0 + (131.438*temp) - (3.32373E6/temp) &
5137 & - (41729.1*alog(temp))
5138 a_w_i = exp(mu_diff/(R_uni*temp))
5139 delta_aw = satw - a_w_i
5140 log_J_rate = -906.7 + (8502.0*delta_aw) &
5141 & - (26924.0*delta_aw*delta_aw) &
5142 & + (29180.0*delta_aw*delta_aw*delta_aw)
5143 log_J_rate = MIN(20.0, log_J_rate)
5144 J_rate = 10.**log_J_rate ! cm-3 s-1
5145 prob_h = MIN(1.-exp(-J_rate*ar_volume*DT), 1.)
5146 if (prob_h .gt. 0.) then
5147 xni = MIN(prob_h*naero, 1000.E3)
5150 iceKoop = MAX(0.0, xni)
5152 end FUNCTION iceKoop
5154 !+---+-----------------------------------------------------------------+
5155 !.. Helper routine for Phillips et al (2008) ice nucleation. Trude
5157 REAL FUNCTION delta_p (yy, y1, y2, aa, bb)
5160 REAL, INTENT(IN):: yy, y1, y2, aa, bb
5161 REAL:: dab, A, B, a0, a1, a2, a3
5163 A = 6.*(aa-bb)/((y2-y1)*(y2-y1)*(y2-y1))
5164 B = aa+(A*y1*y1*y1/6.)-(A*y1*y1*y2*0.5)
5172 else if (yy.ge.y2) then
5175 dab = a0+(a1*yy)+(a2*yy*yy)+(a3*yy*yy*yy)
5186 END FUNCTION delta_p
5188 !+---+-----------------------------------------------------------------+
5191 !+---+-----------------------------------------------------------------+
5192 !..Compute _radiation_ effective radii of cloud water, ice, and snow.
5193 !.. These are entirely consistent with microphysics assumptions, not
5194 !.. constant or otherwise ad hoc as is internal to most radiation
5195 !.. schemes. Since only the smallest snowflakes should impact
5196 !.. radiation, compute from first portion of complicated Field number
5197 !.. distribution, not the second part, which is the larger sizes.
5198 !+---+-----------------------------------------------------------------+
5200 subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
5201 & re_qc1d, re_qi1d, re_qs1d, kts, kte)
5206 INTEGER, INTENT(IN):: kts, kte
5207 REAL, DIMENSION(kts:kte), INTENT(IN):: &
5208 & t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d
5209 REAL, DIMENSION(kts:kte), INTENT(INOUT):: re_qc1d, re_qi1d, re_qs1d
5212 REAL, DIMENSION(kts:kte):: rho, rc, nc, ri, ni, rs
5213 REAL:: smo2, smob, smoc
5214 REAL:: tc0, loga_, a_, b_
5215 DOUBLE PRECISION:: lamc, lami
5216 LOGICAL:: has_qc, has_qi, has_qs
5218 real, dimension(15), parameter:: g_ratio = (/24,60,120,210,336, &
5219 & 504,720,990,1320,1716,2184,2730,3360,4080,4896/)
5226 rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622))
5227 rc(k) = MAX(R1, qc1d(k)*rho(k))
5228 nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max))
5229 if (.NOT. is_aerosol_aware) nc(k) = Nt_c
5230 if (rc(k).gt.R1 .and. nc(k).gt.R2) has_qc = .true.
5231 ri(k) = MAX(R1, qi1d(k)*rho(k))
5232 ni(k) = MAX(R2, ni1d(k)*rho(k))
5233 if (ri(k).gt.R1 .and. ni(k).gt.R2) has_qi = .true.
5234 rs(k) = MAX(R1, qs1d(k)*rho(k))
5235 if (rs(k).gt.R1) has_qs = .true.
5240 re_qc1d(k) = RE_QC_BG
5241 if (rc(k).le.R1 .or. nc(k).le.R2) CYCLE
5242 if (nc(k).lt.100) then
5244 elseif (nc(k).gt.1.E10) then
5247 inu_c = MIN(15, NINT(1000.E6/nc(k)) + 2)
5249 lamc = (nc(k)*am_r*g_ratio(inu_c)/rc(k))**obmr
5250 re_qc1d(k) = MAX(2.51E-6, MIN(SNGL(0.5D0 * DBLE(3.+inu_c)/lamc), 50.E-6))
5256 re_qi1d(k) = RE_QI_BG
5257 if (ri(k).le.R1 .or. ni(k).le.R2) CYCLE
5258 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
5259 re_qi1d(k) = MAX(5.01E-6, MIN(SNGL(0.5D0 * DBLE(3.+mu_i)/lami), 125.E-6))
5265 re_qs1d(k) = RE_QS_BG
5266 if (rs(k).le.R1) CYCLE
5267 tc0 = MIN(-0.1, t1d(k)-273.15)
5270 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
5271 !.. then we must compute actual 2nd moment and use as reference.
5272 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
5275 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
5276 & + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
5277 & + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
5278 & + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
5279 & + sa(10)*bm_s*bm_s*bm_s
5281 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
5282 & + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
5283 & + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
5284 & + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
5285 & + sb(10)*bm_s*bm_s*bm_s
5286 smo2 = (smob/a_)**(1./b_)
5288 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
5289 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
5290 & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
5291 & + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
5292 & + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
5293 & + sa(10)*cse(1)*cse(1)*cse(1)
5295 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
5296 & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
5297 & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
5298 & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
5299 smoc = a_ * smo2**b_
5300 re_qs1d(k) = MAX(10.01E-6, MIN(0.5*(smoc/smob), 999.E-6))
5304 end subroutine calc_effectRad
5306 !+---+-----------------------------------------------------------------+
5307 !..Compute radar reflectivity assuming 10 cm wavelength radar and using
5308 !.. Rayleigh approximation. Only complication is melted snow/graupel
5309 !.. which we treat as water-coated ice spheres and use Uli Blahak's
5310 !.. library of routines. The meltwater fraction is simply the amount
5311 !.. of frozen species remaining from what initially existed at the
5312 !.. melting level interface.
5313 !+---+-----------------------------------------------------------------+
5315 subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, &
5316 t1d, p1d, dBZ, kts, kte, ii, jj, ke_diag)
5321 INTEGER, INTENT(IN):: kts, kte, ii, jj, ke_diag
5322 REAL, DIMENSION(kts:kte), INTENT(IN):: &
5323 qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, t1d, p1d
5324 REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
5325 ! REAL, DIMENSION(kts:kte), INTENT(INOUT):: vt_dBZ
5328 REAL, DIMENSION(kts:kte):: temp, pres, qv, rho, rhof
5329 REAL, DIMENSION(kts:kte):: rc, rr, nr, rs, rg
5331 DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g
5332 REAL, DIMENSION(kts:kte):: mvd_r
5333 REAL, DIMENSION(kts:kte):: smob, smo2, smoc, smoz
5334 REAL:: oM3, M0, Mrat, slam1, slam2, xDs
5335 REAL:: ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts
5336 REAL:: vtr_dbz_wt, vts_dbz_wt, vtg_dbz_wt
5338 REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
5340 DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamr, lamg
5341 REAL:: a_, b_, loga_, tc0
5342 DOUBLE PRECISION:: fmelt_s, fmelt_g
5344 INTEGER:: i, k, k_0, kbot, n, ktop
5346 LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
5348 DOUBLE PRECISION:: cback, x, eta, f_d
5349 REAL:: xslw1, ygra1, zans1
5358 !+---+-----------------------------------------------------------------+
5359 !..Put column of data into local arrays.
5360 !+---+-----------------------------------------------------------------+
5364 qv(k) = MAX(1.E-10, qv1d(k))
5366 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
5367 rhof(k) = SQRT(RHO_NOT/rho(k))
5368 rc(k) = MAX(R1, qc1d(k)*rho(k))
5369 if (qr1d(k) .gt. R1) then
5370 rr(k) = qr1d(k)*rho(k)
5371 nr(k) = MAX(R2, nr1d(k)*rho(k))
5372 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
5374 N0_r(k) = nr(k)*org2*lamr**cre(2)
5375 mvd_r(k) = (3.0 + mu_r + 0.672) * ilamr(k)
5383 if (qs1d(k) .gt. R2) then
5384 rs(k) = qs1d(k)*rho(k)
5390 if (qg1d(k) .gt. R2) then
5391 rg(k) = qg1d(k)*rho(k)
5399 !+---+-----------------------------------------------------------------+
5400 !..Calculate y-intercept, slope, and useful moments for snow.
5401 !+---+-----------------------------------------------------------------+
5408 if ( ( ke_diag > kts .and. ANY(L_qs .eqv. .true.) ) .or. &
5409 (ke_diag == kts .and. L_qs(kts) .eqv. .true. ) ) then
5410 do k = kts, ke_diag ! kte
5411 if (.not. L_qs(k)) CYCLE
5412 tc0 = MIN(-0.1, temp(k)-273.15)
5413 smob(k) = rs(k)*oams
5415 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
5416 !.. then we must compute actual 2nd moment and use as reference.
5417 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
5420 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
5421 & + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
5422 & + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
5423 & + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
5424 & + sa(10)*bm_s*bm_s*bm_s
5426 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
5427 & + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
5428 & + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
5429 & + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
5430 & + sb(10)*bm_s*bm_s*bm_s
5431 smo2(k) = (smob(k)/a_)**(1./b_)
5434 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
5435 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
5436 & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
5437 & + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
5438 & + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
5439 & + sa(10)*cse(1)*cse(1)*cse(1)
5441 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
5442 & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
5443 & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
5444 & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
5445 smoc(k) = a_ * smo2(k)**b_
5447 !..Calculate bm_s*2 (th) moment. Useful for reflectivity.
5448 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(3) &
5449 & + sa(4)*tc0*cse(3) + sa(5)*tc0*tc0 &
5450 & + sa(6)*cse(3)*cse(3) + sa(7)*tc0*tc0*cse(3) &
5451 & + sa(8)*tc0*cse(3)*cse(3) + sa(9)*tc0*tc0*tc0 &
5452 & + sa(10)*cse(3)*cse(3)*cse(3)
5454 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(3) + sb(4)*tc0*cse(3) &
5455 & + sb(5)*tc0*tc0 + sb(6)*cse(3)*cse(3) &
5456 & + sb(7)*tc0*tc0*cse(3) + sb(8)*tc0*cse(3)*cse(3) &
5457 & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(3)*cse(3)*cse(3)
5458 smoz(k) = a_ * smo2(k)**b_
5462 !+---+-----------------------------------------------------------------+
5463 !..Calculate y-intercept, slope values for graupel.
5464 !+---+-----------------------------------------------------------------+
5466 if (ANY(L_qg .eqv. .true.)) then
5468 ygra1 = alog10(max(1.E-9, rg(k)))
5469 zans1 = 3.0 + 2./7.*(ygra1+8.)
5470 zans1 = MAX(2., MIN(zans1, 6.))
5471 N0_exp = 10.**(zans1)
5472 lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
5473 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
5475 N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2)
5479 !+---+-----------------------------------------------------------------+
5480 !..Locate K-level of start of melting (k_0 is level above).
5481 !+---+-----------------------------------------------------------------+
5484 do k = kte-1, kts, -1
5485 if ( (temp(k).gt.273.15) .and. L_qr(k) &
5486 & .and. (L_qs(k+1).or.L_qg(k+1)) ) then
5494 ! Set loop limit for wet ice according to whether the full 3D field is needed or just k=1
5495 k_0loop = Min(k_0, ke_diag+1)
5497 !+---+-----------------------------------------------------------------+
5498 !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
5499 !.. and non-water-coated snow and graupel when below freezing are
5500 !.. simple. Integrations of m(D)*m(D)*N(D)*dD.
5501 !+---+-----------------------------------------------------------------+
5503 do k = kts, ke_diag ! kte
5506 ze_graupel(k) = 1.e-22
5507 if (L_qr(k)) ze_rain(k) = N0_r(k)*crg(4)*ilamr(k)**cre(4)
5508 if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) &
5509 & * (am_s/900.0)*(am_s/900.0)*smoz(k)
5510 if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) &
5511 & * (am_g/900.0)*(am_g/900.0) &
5512 & * N0_g(k)*cgg(4)*ilamg(k)**cge(4)
5515 !+---+-----------------------------------------------------------------+
5516 !..Special case of melting ice (snow/graupel) particles. Assume the
5517 !.. ice is surrounded by the liquid water. Fraction of meltwater is
5518 !.. extremely simple based on amount found above the melting level.
5519 !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
5521 !+---+-----------------------------------------------------------------+
5523 if (.not. iiwarm .and. melti .and. k_0.ge.2) then
5524 do k = k_0loop-1, kts, -1
5526 !..Reflectivity contributed by melting snow
5527 if (L_qs(k) .and. L_qs(k_0) ) then
5528 fmelt_s = MAX(0.05d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
5532 Mrat = smob(k)*M0*M0*M0
5536 x = am_s * xxDs(n)**bm_s
5537 call rayleigh_soak_wetgraupel (x, DBLE(ocms), DBLE(obms), &
5538 & fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
5539 & CBACK, mixingrulestring_s, matrixstring_s, &
5540 & inclusionstring_s, hoststring_s, &
5541 & hostmatrixstring_s, hostinclusionstring_s)
5542 f_d = Mrat*(Kap0*DEXP(-slam1*xxDs(n)) &
5543 & + Kap1*(M0*xxDs(n))**mu_s * DEXP(-slam2*xxDs(n)))
5544 eta = eta + f_d * CBACK * simpson(n) * xdts(n)
5546 ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
5549 !..Reflectivity contributed by melting graupel
5551 if (L_qg(k) .and. L_qg(k_0) ) then
5552 fmelt_g = MAX(0.05d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0))
5556 x = am_g * xxDg(n)**bm_g
5557 call rayleigh_soak_wetgraupel (x, DBLE(ocmg), DBLE(obmg), &
5558 & fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
5559 & CBACK, mixingrulestring_g, matrixstring_g, &
5560 & inclusionstring_g, hoststring_g, &
5561 & hostmatrixstring_g, hostinclusionstring_g)
5562 f_d = N0_g(k)*xxDg(n)**mu_g * DEXP(-lamg*xxDg(n))
5563 eta = eta + f_d * CBACK * simpson(n) * xdtg(n)
5565 ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
5571 do k = ke_diag, kts, -1
5572 dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
5576 !..Reflectivity-weighted terminal velocity (snow, rain, graupel, mix).
5577 ! do k = kte, kts, -1
5579 ! if (rs(k).gt.R2) then
5580 ! Mrat = smob(k) / smoc(k)
5581 ! ils1 = 1./(Mrat*Lam0 + fv_s)
5582 ! ils2 = 1./(Mrat*Lam1 + fv_s)
5583 ! t1_vts = Kap0*csg(5)*ils1**cse(5)
5584 ! t2_vts = Kap1*Mrat**mu_s*csg(11)*ils2**cse(11)
5585 ! ils1 = 1./(Mrat*Lam0)
5586 ! ils2 = 1./(Mrat*Lam1)
5587 ! t3_vts = Kap0*csg(6)*ils1**cse(6)
5588 ! t4_vts = Kap1*Mrat**mu_s*csg(12)*ils2**cse(12)
5589 ! vts_dbz_wt = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts)
5590 ! if (temp(k).ge.273.15 .and. temp(k).lt.275.15) then
5591 ! vts_dbz_wt = vts_dbz_wt*1.5
5592 ! elseif (temp(k).ge.275.15) then
5593 ! vts_dbz_wt = vts_dbz_wt*2.0
5596 ! vts_dbz_wt = 1.E-3
5599 ! if (rr(k).gt.R1) then
5600 ! lamr = 1./ilamr(k)
5601 ! vtr_dbz_wt = rhof(k)*av_r*crg(13)*(lamr+fv_r)**(-cre(13)) &
5602 ! & / (crg(4)*lamr**(-cre(4)))
5604 ! vtr_dbz_wt = 1.E-3
5607 ! if (rg(k).gt.R2) then
5608 ! lamg = 1./ilamg(k)
5609 ! vtg_dbz_wt = rhof(k)*av_g*cgg(5)*lamg**(-cge(5)) &
5610 ! & / (cgg(4)*lamg**(-cge(4)))
5612 ! vtg_dbz_wt = 1.E-3
5615 ! vt_dBZ(k) = (vts_dbz_wt*ze_snow(k) + vtr_dbz_wt*ze_rain(k) &
5616 ! & + vtg_dbz_wt*ze_graupel(k)) &
5617 ! & / (ze_rain(k)+ze_snow(k)+ze_graupel(k))
5620 end subroutine calc_refl10cm
5622 !+---+-----------------------------------------------------------------+
5624 !+---+-----------------------------------------------------------------+
5625 END MODULE module_mp_thompson
5626 !+---+-----------------------------------------------------------------+