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