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