CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
module_mp_tempo_main.F90
1! 1D TEMPO microphysics scheme
2!=================================================================================================================
3module module_mp_tempo_main
4
5 use module_mp_tempo_params
6 use module_mp_tempo_utils, only: rslf, rsif
7
8#if defined(mpas)
9 use mpas_kind_types, only: wp => rkind, sp => r4kind, dp => r8kind
10#elif defined(standalone)
11 use machine, only: wp => kind_phys, sp => kind_sngl_prec, dp => kind_dbl_prec
12#else
13 use machine, only: wp => kind_phys, sp => kind_sngl_prec, dp => kind_dbl_prec
14#define ccpp_default 1
15#endif
16
17 implicit none
18
19contains
20 !=================================================================================================================
21 ! This subroutine computes the moisture tendencies of water vapor, cloud droplets, rain, cloud ice (pristine),
22 ! snow, and graupel. Previously this code was based on Reisner et al (1998), but few of those pieces remain.
23 ! A complete description is now found in Thompson et al. (2004, 2008), Thompson and Eidhammer (2014),
24 ! and Jensen et al. (2023).
25
26 subroutine mp_tempo_main(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, ni1d, nr1d, nc1d, ng1d, &
27 nwfa1d, nifa1d, t1d, p1d, w1d, dzq, pptrain, pptsnow, pptgraul, pptice, &
28#if defined(mpas)
29 rainprod, evapprod, &
30#endif
31
32#if defined(ccpp_default)
33 ! Extended diagnostics, most arrays only
34 ! allocated if ext_diag flag is .true.
35 ext_diag, sedi_semi, &
36 prw_vcdc1, prw_vcde1, &
37 tpri_inu1, tpri_ide1_d, tpri_ide1_s, tprs_ide1, &
38 tprs_sde1_d, tprs_sde1_s, &
39 tprg_gde1_d, tprg_gde1_s, tpri_iha1, tpri_wfz1, &
40 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1,&
41 tprg_rcs1, tprs_rcs1, tprr_rci1, &
42 tprg_rcg1, tprw_vcd1_c, &
43 tprw_vcd1_e, tprr_sml1, tprr_gml1, tprr_rcg1, &
44 tprr_rcs1, tprv_rev1, &
45 tten1, qvten1, qrten1, qsten1, &
46 qgten1, qiten1, niten1, nrten1, ncten1, qcten1, &
47#endif
48 decfl, pfil1, pfll1, &
49 lsml, rand1, rand2, rand3, &
50 kts, kte, dt, ii, jj, &
51 configs)
52
53#if defined(ccpp_default) && defined (MPI)
54 use mpi_f08
55#endif
56
57 ! Subroutine arguments
58 integer, intent(in) :: kts, kte, ii, jj
59 real(wp), dimension(kts:kte), intent(inout) :: qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, &
60 ni1d, nr1d, nc1d, ng1d, nwfa1d, nifa1d, t1d
61 real(wp), dimension(kts:kte), intent(in) :: p1d, w1d, dzq
62 real(wp), intent(inout) :: pptrain, pptsnow, pptgraul, pptice
63 real(wp), intent(in) :: dt
64
65 type(ty_tempo_cfg), intent(in) :: configs
66
67 integer, intent(in), optional :: lsml
68 real(wp), intent(in), optional :: rand1, rand2, rand3
69 real(wp), dimension(kts:kte), intent(out), optional :: pfil1, pfll1
70 integer, intent(in), optional :: decfl
71#if defined(ccpp_default)
72 ! Extended diagnostics, most arrays only allocated if ext_diag is true
73 logical, intent(in) :: ext_diag
74 logical, intent(in) :: sedi_semi
75 real(wp), optional, dimension(:), intent(out) :: &
76 prw_vcdc1, &
77 prw_vcde1, tpri_inu1, tpri_ide1_d, &
78 tpri_ide1_s, tprs_ide1, &
79 tprs_sde1_d, tprs_sde1_s, tprg_gde1_d, &
80 tprg_gde1_s, tpri_iha1, tpri_wfz1, &
81 tpri_rfz1, tprg_rfz1, tprs_scw1, tprg_scw1, &
82 tprg_rcs1, tprs_rcs1, &
83 tprr_rci1, tprg_rcg1, &
84 tprw_vcd1_c, tprw_vcd1_e, tprr_sml1, &
85 tprr_gml1, tprr_rcg1, &
86 tprr_rcs1, tprv_rev1, tten1, qvten1, &
87 qrten1, qsten1, qgten1, qiten1, niten1, &
88 nrten1, ncten1, qcten1
89#endif
90
91#if defined(mpas)
92 real(wp), dimension(kts:kte), intent(inout) :: rainprod, evapprod
93#endif
94
95 !=================================================================================================================
96 ! Local variables
97 real(wp), dimension(kts:kte) :: tten, qvten, qcten, qiten, qrten, &
98 qsten, qgten, qbten, niten, nrten, ncten, ngten, nwfaten, nifaten
99
100 real(dp), dimension(kts:kte) :: prw_vcd
101
102 real(dp), dimension(kts:kte) :: pnc_wcd, pnc_wau, pnc_rcw, pnc_scw, pnc_gcw
103
104 real(dp), dimension(kts:kte) :: pna_rca, pna_sca, pna_gca, pnd_rcd, pnd_scd, pnd_gcd
105
106 real(dp), dimension(kts:kte) :: prr_wau, prr_rcw, prr_rcs, &
107 prr_rcg, prr_sml, prr_gml, &
108 prr_rci, prv_rev, &
109 pnr_wau, pnr_rcs, pnr_rcg, &
110 pnr_rci, pnr_sml, pnr_gml, &
111 pnr_rev, pnr_rcr, pnr_rfz
112
113 real(dp), dimension(kts:kte) :: pri_inu, pni_inu, pri_ihm, &
114 pni_ihm, pri_wfz, pni_wfz, &
115 pri_rfz, pni_rfz, pri_ide, &
116 pni_ide, pri_rci, pni_rci, &
117 pni_sci, pni_iau, pri_iha, pni_iha
118
119 real(dp), dimension(kts:kte) :: prs_iau, prs_sci, prs_rcs, prs_scw, prs_sde, prs_ihm, prs_ide
120
121 real(dp), dimension(kts:kte) :: prg_scw, prg_rfz, prg_gde, &
122 prg_gcw, prg_rci, prg_rcs, prg_rcg, prg_ihm, &
123 png_rcs, png_rcg, png_scw, png_gde, &
124 pbg_scw, pbg_rfz, pbg_gcw, pbg_rci, pbg_rcs, pbg_rcg, &
125 pbg_sml, pbg_gml
126
127 real(dp), parameter :: zerod0 = 0.0d0
128
129 real(wp), dimension(kts:kte) :: pfll, pfil, pdummy
130 real(wp) :: dtcfl, rainsfc, graulsfc, orhodt
131 integer :: niter
132 real(wp), dimension(kts:kte) :: rr_tmp, nr_tmp, rg_tmp
133
134 real(wp), dimension(kts:kte) :: temp, twet, pres, qv
135 real(wp), dimension(kts:kte) :: rc, ri, rr, rs, rg, rb
136 real(wp), dimension(kts:kte) :: ni, nr, nc, ns, ng, nwfa, nifa
137 real(wp), dimension(kts:kte) :: rho, rhof, rhof2
138 real(wp), dimension(kts:kte) :: qvs, qvsi, delqvs
139 real(wp), dimension(kts:kte) :: satw, sati, ssatw, ssati
140 real(wp), dimension(kts:kte) :: diffu, visco, vsc2, &
141 tcond, lvap, ocp, lvt2
142
143 real(dp), dimension(kts:kte) :: ilamr, ilamg, n0_r, n0_g
144 real(dp) :: n0_melt
145 real(wp), dimension(kts:kte) :: mvd_r, mvd_c, mvd_g
146 real(wp), dimension(kts:kte) :: smob, smo2, smo1, smo0, &
147 smoc, smod, smoe, smof, smog
148
149 real(wp), dimension(kts:kte) :: sed_r, sed_s, sed_g, sed_i, sed_n, sed_c, sed_b
150
151 real(wp) :: rgvm, delta_tp, orho, lfus2
152 real(wp), dimension(5):: onstep
153 real(dp) :: n0_exp, n0_min, lam_exp, lamc, lamr, lamg
154 real(dp) :: lami, ilami, ilamc
155 real(wp) :: xdc, dc_b, dc_g, xdi, xdr, xds, xdg, ds_m, dg_m
156 real(dp) :: dr_star, dc_star
157 real(wp) :: zeta1, zeta, taud, tau
158 real(wp) :: stoke_r, stoke_s, stoke_g, stoke_i
159 real(wp) :: vti, vtr, vts, vtg, vtc
160 real(wp) :: xrho_g, afall, vtg1, vtg2
161 real(wp) :: bfall = 3*b_coeff - 1
162 real(wp), dimension(kts:kte+1) :: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk, vtngk, vtck, vtnck
163 real(wp), dimension(kts:kte) :: vts_boost
164 real(wp) :: m0, slam1, slam2
165 real(wp) :: mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, c_snow
166 real(wp) :: a_, b_, loga_, a1, a2, tf
167 real(wp) :: tempc, tc0, r_mvd1, r_mvd2, xkrat
168 real(wp) :: dew_t, tlcl, the
169 real(wp) :: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr, xrg, xng, xrb
170 real(wp) :: xsat, rate_max, sump, ratio
171 real(wp) :: clap, fcd, dfcd
172 real(wp) :: otemp, rvs, rvs_p, rvs_pp, gamsc, alphsc, t1_evap, t1_subl
173 real(wp) :: r_frac, g_frac, const_Ri, rime_dens
174 real(wp) :: Ef_rw, Ef_sw, Ef_gw, Ef_rr
175 real(wp) :: Ef_ra, Ef_sa, Ef_ga
176 real(wp) :: dtsave, odts, odt, odzq, hgt_agl, SR
177 real(wp) :: xslw1, ygra1, zans1, eva_factor
178 real(wp) :: melt_f, rand
179 integer :: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq, k_melting
180 integer, dimension(5) :: ksed1
181 integer :: nir, nis, nig, nii, nic, niin
182 integer :: idx_tc, idx_t, idx_s, idx_g1, idx_g, idx_r1, idx_r, &
183 idx_i1, idx_i, idx_c, idx, idx_d, idx_n, idx_in
184 integer, dimension(kts:kte) :: idx_bg, idx_table
185
186 logical :: melti, no_micro
187 logical, dimension(kts:kte) :: l_qc, l_qi, l_qr, l_qs, l_qg
188 logical :: debug_flag
189 character*256 :: mp_debug
190 integer :: nu_c, decfl_
191
192 !=================================================================================================================
193
194 debug_flag = .false.
195
196 no_micro = .true.
197 dtsave = dt
198 odt = 1./dt
199 odts = 1./dtsave
200 iexfrq = 1
201 rand = 0.0
202 decfl_ = 10
203 if (present(decfl)) decfl_ = decfl
204
205#if defined(ccpp_default)
206 ! Transition value of coefficient matching at crossover from cloud ice to snow
207 av_i = av_s * d0s ** (bv_s - bv_i)
208#endif
209
210 !=================================================================================================================
211 ! Source/sink terms. First 2 chars: "pr" represents source/sink of
212 ! mass while "pn" represents source/sink of number. Next char is one
213 ! of "v" for water vapor, "r" for rain, "i" for cloud ice, "w" for
214 ! cloud water, "s" for snow, and "g" for graupel. Next chars
215 ! represent processes: "de" for sublimation/deposition, "ev" for
216 ! evaporation, "fz" for freezing, "ml" for melting, "au" for
217 ! autoconversion, "nu" for ice nucleation, "hm" for Hallet/Mossop
218 ! secondary ice production, and "c" for collection followed by the
219 ! character for the species being collected. ALL of these terms are
220 ! positive (except for deposition/sublimation terms which can switch
221 ! signs based on super/subsaturation) and are treated as negatives
222 ! where necessary in the tendency equations.
223 !=================================================================================================================
224 ! TODO: Put these in derived data type and add initialization subroutine
225 do k = kts, kte
226 tten(k) = 0.
227 qvten(k) = 0.
228 qcten(k) = 0.
229 qiten(k) = 0.
230 qrten(k) = 0.
231 qsten(k) = 0.
232 qgten(k) = 0.
233 ngten(k) = 0.
234 qbten(k) = 0.
235 niten(k) = 0.
236 nrten(k) = 0.
237 ncten(k) = 0.
238 nwfaten(k) = 0.
239 nifaten(k) = 0.
240
241 prw_vcd(k) = 0.
242
243 pnc_wcd(k) = 0.
244 pnc_wau(k) = 0.
245 pnc_rcw(k) = 0.
246 pnc_scw(k) = 0.
247 pnc_gcw(k) = 0.
248
249 prv_rev(k) = 0.
250 prr_wau(k) = 0.
251 prr_rcw(k) = 0.
252 prr_rcs(k) = 0.
253 prr_rcg(k) = 0.
254 prr_sml(k) = 0.
255 prr_gml(k) = 0.
256 prr_rci(k) = 0.
257 pnr_wau(k) = 0.
258 pnr_rcs(k) = 0.
259 pnr_rcg(k) = 0.
260 pnr_rci(k) = 0.
261 pnr_sml(k) = 0.
262 pnr_gml(k) = 0.
263 pnr_rev(k) = 0.
264 pnr_rcr(k) = 0.
265 pnr_rfz(k) = 0.
266
267 pri_inu(k) = 0.
268 pni_inu(k) = 0.
269 pri_ihm(k) = 0.
270 pni_ihm(k) = 0.
271 pri_wfz(k) = 0.
272 pni_wfz(k) = 0.
273 pri_rfz(k) = 0.
274 pni_rfz(k) = 0.
275 pri_ide(k) = 0.
276 pni_ide(k) = 0.
277 pri_rci(k) = 0.
278 pni_rci(k) = 0.
279 pni_sci(k) = 0.
280 pni_iau(k) = 0.
281 pri_iha(k) = 0.
282 pni_iha(k) = 0.
283
284 prs_iau(k) = 0.
285 prs_sci(k) = 0.
286 prs_rcs(k) = 0.
287 prs_scw(k) = 0.
288 prs_sde(k) = 0.
289 prs_ihm(k) = 0.
290 prs_ide(k) = 0.
291
292 prg_scw(k) = 0.
293 prg_rfz(k) = 0.
294 prg_gde(k) = 0.
295 prg_gcw(k) = 0.
296 prg_rci(k) = 0.
297 prg_rcs(k) = 0.
298 prg_rcg(k) = 0.
299 prg_ihm(k) = 0.
300 ! new source/sink terms for 3-moment graupel
301 png_scw(k) = 0.
302 png_rcs(k) = 0.
303 png_rcg(k) = 0.
304 png_gde(k) = 0.
305
306 pbg_scw(k) = 0.
307 pbg_rfz(k) = 0.
308 pbg_gcw(k) = 0.
309 pbg_rci(k) = 0.
310 pbg_rcs(k) = 0.
311 pbg_rcg(k) = 0.
312 pbg_sml(k) = 0.
313 pbg_gml(k) = 0.
314
315 pna_rca(k) = 0.
316 pna_sca(k) = 0.
317 pna_gca(k) = 0.
318
319 pnd_rcd(k) = 0.
320 pnd_scd(k) = 0.
321 pnd_gcd(k) = 0.
322
323 if (present(pfil1)) pfil1(k) = 0.
324 if (present(pfll1)) pfll1(k) = 0.
325 pfil(k) = 0.
326 pfll(k) = 0.
327 pdummy(k) = 0.
328 enddo
329#if defined(mpas)
330 do k = kts, kte
331 rainprod(k) = 0.
332 evapprod(k) = 0.
333 enddo
334#endif
335
336#if defined(ccpp_default)
337 !Diagnostics
338 if (ext_diag) then
339 do k = kts, kte
340 !vtsk1(k) = 0.
341 !txrc1(k) = 0.
342 !txri1(k) = 0.
343 prw_vcdc1(k) = 0.
344 prw_vcde1(k) = 0.
345 tpri_inu1(k) = 0.
346 tpri_ide1_d(k) = 0.
347 tpri_ide1_s(k) = 0.
348 tprs_ide1(k) = 0.
349 tprs_sde1_d(k) = 0.
350 tprs_sde1_s(k) = 0.
351 tprg_gde1_d(k) = 0.
352 tprg_gde1_s(k) = 0.
353 tpri_iha1(k) = 0.
354 tpri_wfz1(k) = 0.
355 tpri_rfz1(k) = 0.
356 tprg_rfz1(k) = 0.
357 tprg_scw1(k) = 0.
358 tprs_scw1(k) = 0.
359 tprg_rcs1(k) = 0.
360 tprs_rcs1(k) = 0.
361 tprr_rci1(k) = 0.
362 tprg_rcg1(k) = 0.
363 tprw_vcd1_c(k) = 0.
364 tprw_vcd1_e(k) = 0.
365 tprr_sml1(k) = 0.
366 tprr_gml1(k) = 0.
367 tprr_rcg1(k) = 0.
368 tprr_rcs1(k) = 0.
369 tprv_rev1(k) = 0.
370 tten1(k) = 0.
371 qvten1(k) = 0.
372 qrten1(k) = 0.
373 qsten1(k) = 0.
374 qgten1(k) = 0.
375 qiten1(k) = 0.
376 niten1(k) = 0.
377 nrten1(k) = 0.
378 ncten1(k) = 0.
379 qcten1(k) = 0.
380 enddo
381 endif
382#endif
383 !..Bug fix (2016Jun15), prevent use of uninitialized value(s) of snow moments.
384 do k = kts, kte
385 smo0(k) = 0.
386 smo1(k) = 0.
387 smo2(k) = 0.
388 smob(k) = 0.
389 smoc(k) = 0.
390 smod(k) = 0.
391 smoe(k) = 0.
392 smof(k) = 0.
393 smog(k) = 0.
394 ns(k) = 0.
395 mvd_r(k) = 0.
396 mvd_c(k) = 0.
397 enddo
398
399 !=================================================================================================================
400 ! Convert microphysics variables to concentrations (kg / m^3 and # / m^3)
401 do k = kts, kte
402 temp(k) = t1d(k)
403 qv(k) = max(min_qv, qv1d(k))
404 pres(k) = p1d(k)
405 rho(k) = roverrv*pres(k)/(r*temp(k)*(qv(k)+roverrv))
406
407 ! CCPP version has rho(k) multiplier for min and max
408 ! nwfa(k) = max(11.1e6, min(9999.e6, nwfa1d(k)*rho(k)))
409 ! nifa(k) = max(nain1*0.01, min(9999.e6, nifa1d(k)*rho(k)))
410 nwfa(k) = max(nwfa_default*rho(k), min(aero_max*rho(k), nwfa1d(k)*rho(k)))
411 nifa(k) = max(nifa_default*rho(k), min(aero_max*rho(k), nifa1d(k)*rho(k)))
412
413 ! From CCPP version
414 mvd_r(k) = d0r
415 mvd_c(k) = d0c
416
417 if (qc1d(k) .gt. r1) then
418 no_micro = .false.
419 rc(k) = qc1d(k)*rho(k)
420 nc(k) = max(2., min(nc1d(k)*rho(k), nt_c_max))
421 l_qc(k) = .true.
422 if (nc(k).gt.10000.e6) then
423 nu_c = 2
424 elseif (nc(k).lt.100.) then
425 nu_c = 15
426 else
427 nu_c = nint(nu_c_scale/nc(k)) + 2
428 rand = 0.0
429 if (present(rand2)) then
430 rand = rand2
431 endif
432 nu_c = max(2, min(nu_c+nint(rand), 15))
433 endif
434 lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
435 xdc = (bm_r + nu_c + 1.) / lamc
436 if (xdc.lt. d0c) then
437 lamc = cce(2,nu_c)/d0c
438 elseif (xdc.gt. d0r*2.) then
439 lamc = cce(2,nu_c)/(d0r*2.)
440 endif
441 nc(k) = min(real(nt_c_max, kind=dp), ccg(1,nu_c)*ocg2(nu_c)*rc(k) / am_r*lamc**bm_r)
442 ! CCPP version has different values of Nt_c for land/ocean
443 if (.not.(configs%aerosol_aware .or. merra2_aerosol_aware)) then
444 nc(k) = nt_c
445 if (present(lsml)) then
446 if (lsml == 1) then
447 nc(k) = nt_c_l
448 else
449 nc(k) = nt_c_o
450 endif
451 endif
452 endif
453 else
454 qc1d(k) = 0.0
455 nc1d(k) = 0.0
456 rc(k) = r1
457 nc(k) = 2.
458 l_qc(k) = .false.
459 endif
460
461 if (qi1d(k) .gt. r1) then
462 no_micro = .false.
463 ri(k) = qi1d(k)*rho(k)
464 ni(k) = max(r2, ni1d(k)*rho(k))
465 if (ni(k).le. r2) then
466 lami = cie(2)/5.e-6
467 ni(k) = min(max_ni, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
468 endif
469 l_qi(k) = .true.
470 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
471 ilami = 1./lami
472 xdi = (bm_i + mu_i + 1.) * ilami
473 if (xdi.lt. 5.e-6) then
474 lami = cie(2)/5.e-6
475 ni(k) = min(max_ni, cig(1)*oig2*ri(k)/am_i*lami**bm_i)
476 elseif (xdi.gt. 300.e-6) then
477 lami = cie(2)/300.e-6
478 ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i
479 endif
480 else
481 qi1d(k) = 0.0
482 ni1d(k) = 0.0
483 ri(k) = r1
484 ni(k) = r2
485 l_qi(k) = .false.
486 endif
487
488 if (qr1d(k) .gt. r1) then
489 no_micro = .false.
490 rr(k) = qr1d(k)*rho(k)
491 nr(k) = max(r2, nr1d(k)*rho(k))
492 if (nr(k).le. r2) then
493 mvd_r(k) = 1.0e-3
494 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
495 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
496 endif
497 l_qr(k) = .true.
498 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
499 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
500 if (mvd_r(k) .gt. 2.5e-3) then
501 mvd_r(k) = 2.5e-3
502 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
503 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
504 elseif (mvd_r(k) .lt. d0r*0.75) then
505 mvd_r(k) = d0r*0.75
506 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
507 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
508 endif
509 else
510 qr1d(k) = 0.0
511 nr1d(k) = 0.0
512 rr(k) = r1
513 nr(k) = r2
514 l_qr(k) = .false.
515 endif
516 if (qs1d(k) .gt. r1) then
517 no_micro = .false.
518 rs(k) = qs1d(k)*rho(k)
519 l_qs(k) = .true.
520 else
521 qs1d(k) = 0.0
522 rs(k) = r1
523 l_qs(k) = .false.
524 endif
525 if (qg1d(k) .gt. r1) then
526 no_micro = .false.
527 l_qg(k) = .true.
528 rg(k) = qg1d(k)*rho(k)
529 ng(k) = max(r2, ng1d(k)*rho(k))
530 rb(k) = max(qg1d(k)/rho_g(nrhg), qb1d(k))
531 rb(k) = min(qg1d(k)/rho_g(1), rb(k))
532 qb1d(k) = rb(k)
533 idx_bg(k) = max(1,min(nint(qg1d(k)/rb(k) *0.01)+1,nrhg))
534 idx_table(k) = idx_bg(k)
535 if (ng(k).le. r2) then
536 mvd_g(k) = 1.5e-3
537 lamg = (3.0 + mu_g + 0.672) / mvd_g(k)
538 ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k))
539 endif
540 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg
541 mvd_g(k) = (3.0 + mu_g + 0.672) / lamg
542 if (mvd_g(k) .gt. 25.4e-3) then
543 mvd_g(k) = 25.4e-3
544 lamg = (3.0 + mu_g + 0.672) / mvd_g(k)
545 ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k))
546 elseif (mvd_g(k) .lt. d0r) then
547 mvd_g(k) = d0r
548 lamg = (3.0 + mu_g + 0.672) / mvd_g(k)
549 ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k))
550 endif
551 else
552 qg1d(k) = 0.0
553 ng1d(k) = 0.0
554 qb1d(k) = 0.0
555 idx_bg(k) = idx_bg1
556 idx_table(k) = idx_bg(k)
557 rg(k) = r1
558 ng(k) = r2
559 rb(k) = r1/rho(k)/rho_g(nrhg)
560 l_qg(k) = .false.
561 endif
562 if (.not. configs%hail_aware) then
563 idx_bg(k) = idx_bg1
564 idx_table(k) = idx_bg(k)
565 ! If dimNRHG = 1, set idx_table(k) = 1,
566 ! otherwise idx_bg1
567 if(.not. using_hail_aware_table) then
568 idx_table(k) = 1
569 endif
570 endif
571 enddo
572
573 ! if (debug_flag) then
574 ! write(mp_debug,*) 'DEBUG-VERBOSE at (i,j) ', ii, ', ', jj
575 ! CALL wrf_debug(550, mp_debug)
576 ! do k = kts, kte
577 ! write(mp_debug, '(a,i3,f8.2,1x,f7.2,1x, 11(1x,e13.6))') &
578 ! & 'VERBOSE: ', k, pres(k)*0.01, temp(k)-273.15, qv(k), rc(k), rr(k), ri(k), rs(k), rg(k), nc(k), nr(k), ni(k), nwfa(k), nifa(k)
579 ! CALL wrf_debug(550, mp_debug)
580 ! enddo
581 ! endif
582
583 !=================================================================================================================
584 ! Derive various thermodynamic variables frequently used.
585 ! Saturation vapor pressure (mixing ratio) over liquid/ice comes from
586 ! Flatau et al. 1992; enthalpy (latent heat) of vaporization from
587 ! Bohren & Albrecht 1998; others from Pruppacher & Klett 1978.
588 do k = kts, kte
589 tempc = temp(k) - 273.15
590 rhof(k) = sqrt(rho_not/rho(k))
591 rhof2(k) = sqrt(rhof(k))
592 qvs(k) = rslf(pres(k), temp(k))
593 delqvs(k) = max(0.0, rslf(pres(k), 273.15)-qv(k))
594 if (tempc .le. 0.0) then
595 qvsi(k) = rsif(pres(k), temp(k))
596 else
597 qvsi(k) = qvs(k)
598 endif
599 satw(k) = qv(k)/qvs(k)
600 sati(k) = qv(k)/qvsi(k)
601 ssatw(k) = satw(k) - 1.
602 ssati(k) = sati(k) - 1.
603 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
604 if (abs(ssati(k)).lt. eps) ssati(k) = 0.0
605 if (no_micro .and. ssati(k).gt. 0.0) no_micro = .false.
606 diffu(k) = 2.11e-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
607 if (tempc .ge. 0.0) then
608 visco(k) = (1.718+0.0049*tempc)*1.0e-5
609 else
610 visco(k) = (1.718+0.0049*tempc-1.2e-5*tempc*tempc)*1.0e-5
611 endif
612 ocp(k) = 1./(cp2*(1.+0.887*qv(k)))
613 vsc2(k) = sqrt(rho(k)/visco(k))
614 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
615 tcond(k) = (5.69 + 0.0168*tempc)*1.0e-5 * 418.936
616 enddo
617
618 !=================================================================================================================
619 !..If no existing hydrometeor species and no chance to initiate ice or
620 !.. condense cloud water, just exit quickly!
621 if (no_micro) return
622
623 !..Calculate y-intercept, slope, and useful moments for snow.
624 if (.not. iiwarm) then
625 do k = kts, kte
626 if (.not. l_qs(k)) cycle
627 tc0 = min(-0.1, temp(k)-273.15)
628 smob(k) = rs(k)*oams
629
630 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
631 !.. then we must compute actual 2nd moment and use as reference.
632 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
633 smo2(k) = smob(k)
634 else
635 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
636 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
637 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
638 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
639 + sa(10)*bm_s*bm_s*bm_s
640 a_ = 10.0**loga_
641 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
642 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
643 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
644 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
645 + sb(10)*bm_s*bm_s*bm_s
646 smo2(k) = (smob(k)/a_)**(1./b_)
647 endif
648
649 !..Calculate 0th moment. Represents snow number concentration.
650 loga_ = sa(1) + sa(2)*tc0 + sa(5)*tc0*tc0 + sa(9)*tc0*tc0*tc0
651 a_ = 10.0**loga_
652 b_ = sb(1) + sb(2)*tc0 + sb(5)*tc0*tc0 + sb(9)*tc0*tc0*tc0
653 smo0(k) = a_ * smo2(k)**b_
654
655 !..Calculate 1st moment. Useful for depositional growth and melting.
656 loga_ = sa(1) + sa(2)*tc0 + sa(3) &
657 + sa(4)*tc0 + sa(5)*tc0*tc0 &
658 + sa(6) + sa(7)*tc0*tc0 &
659 + sa(8)*tc0 + sa(9)*tc0*tc0*tc0 &
660 + sa(10)
661 a_ = 10.0**loga_
662 b_ = sb(1)+ sb(2)*tc0 + sb(3) + sb(4)*tc0 &
663 + sb(5)*tc0*tc0 + sb(6) &
664 + sb(7)*tc0*tc0 + sb(8)*tc0 &
665 + sb(9)*tc0*tc0*tc0 + sb(10)
666 smo1(k) = a_ * smo2(k)**b_
667
668 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
669 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
670 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
671 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
672 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
673 + sa(10)*cse(1)*cse(1)*cse(1)
674 a_ = 10.0**loga_
675 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
676 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
677 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
678 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
679 smoc(k) = a_ * smo2(k)**b_
680 !..Calculate snow number concentration (explicit integral, not smo0)
681 m0 = smob(k)/smoc(k)
682 mrat = smob(k)*m0*m0*m0
683 slam1 = m0 * lam0
684 slam2 = m0 * lam1
685 ns(k) = mrat*kap0/slam1 &
686 + mrat*kap1*m0**mu_s*csg(15)/slam2**cse(15)
687
688 !..Calculate bv_s+2 (th) moment. Useful for riming.
689 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(13) &
690 + sa(4)*tc0*cse(13) + sa(5)*tc0*tc0 &
691 + sa(6)*cse(13)*cse(13) + sa(7)*tc0*tc0*cse(13) &
692 + sa(8)*tc0*cse(13)*cse(13) + sa(9)*tc0*tc0*tc0 &
693 + sa(10)*cse(13)*cse(13)*cse(13)
694 a_ = 10.0**loga_
695 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(13) + sb(4)*tc0*cse(13) &
696 + sb(5)*tc0*tc0 + sb(6)*cse(13)*cse(13) &
697 + sb(7)*tc0*tc0*cse(13) + sb(8)*tc0*cse(13)*cse(13) &
698 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(13)*cse(13)*cse(13)
699 smoe(k) = a_ * smo2(k)**b_
700
701 !..Calculate 1+(bv_s+1)/2 (th) moment. Useful for depositional growth.
702 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(16) &
703 + sa(4)*tc0*cse(16) + sa(5)*tc0*tc0 &
704 + sa(6)*cse(16)*cse(16) + sa(7)*tc0*tc0*cse(16) &
705 + sa(8)*tc0*cse(16)*cse(16) + sa(9)*tc0*tc0*tc0 &
706 + sa(10)*cse(16)*cse(16)*cse(16)
707 a_ = 10.0**loga_
708 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(16) + sb(4)*tc0*cse(16) &
709 + sb(5)*tc0*tc0 + sb(6)*cse(16)*cse(16) &
710 + sb(7)*tc0*tc0*cse(16) + sb(8)*tc0*cse(16)*cse(16) &
711 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(16)*cse(16)*cse(16)
712 smof(k) = a_ * smo2(k)**b_
713 !..Calculate bm_s + bv_s+2 (th) moment. Useful for riming into graupel.
714 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(17) &
715 + sa(4)*tc0*cse(17) + sa(5)*tc0*tc0 &
716 + sa(6)*cse(17)*cse(17) + sa(7)*tc0*tc0*cse(17) &
717 + sa(8)*tc0*cse(17)*cse(17) + sa(9)*tc0*tc0*tc0 &
718 + sa(10)*cse(17)*cse(17)*cse(17)
719 a_ = 10.0**loga_
720 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(17) + sb(4)*tc0*cse(17) &
721 + sb(5)*tc0*tc0 + sb(6)*cse(17)*cse(17) &
722 + sb(7)*tc0*tc0*cse(17) + sb(8)*tc0*cse(17)*cse(17) &
723 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(17)*cse(17)*cse(17)
724 smog(k) = a_ * smo2(k)**b_
725
726 enddo
727
728 !+---+-----------------------------------------------------------------+
729 !..Calculate y-intercept, slope values for graupel.
730 !+---+-----------------------------------------------------------------+
731 do k = kte, kts, -1
732 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg
733 ilamg(k) = 1./lamg
734 n0_g(k) = ng(k)*ogg2*lamg**cge(2,1)
735 enddo
736 ! do k = kte, kts, -1
737 ! ygra1 = alog10(max(1.e-9, rg(k)))
738 ! zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1
739 ! N0_exp = 10.**(zans1)
740 ! N0_exp = max(dble(gonv_min), min(N0_exp, dble(gonv_max)))
741 ! lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
742 ! lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
743 ! ilamg(k) = 1./lamg
744 ! N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2)
745 ! enddo
746 endif
747
748 !+---+-----------------------------------------------------------------+
749 !..Calculate y-intercept, slope values for rain.
750 !+---+-----------------------------------------------------------------+
751 do k = kte, kts, -1
752 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
753 ilamr(k) = 1./lamr
754 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
755 n0_r(k) = nr(k)*org2*lamr**cre(2)
756 enddo
757
758 !=================================================================================================================
759 !..Compute warm-rain process terms (except evap done later).
760 !+---+-----------------------------------------------------------------+
761
762 do k = kts, kte
763
764 !..Rain self-collection follows Seifert, 1994 and drop break-up
765 !.. follows Verlinde and Cotton, 1993. RAIN2M
766 if (l_qr(k) .and. mvd_r(k).gt. d0r) then
767
768 ef_rr = max(-0.1, 1.0 - exp(2300.0*(mvd_r(k)-1950.0e-6)))
769!!! Ef_rr = 1.0 - exp(2300.0*(mvd_r(k)-1950.0E-6))
770 pnr_rcr(k) = ef_rr * 2.0*nr(k)*rr(k)
771 endif
772
773 mvd_c(k) = d0c
774 if (l_qc(k)) then
775 if (nc(k).gt.10000.e6) then
776 nu_c = 2
777 elseif (nc(k).lt.100.) then
778 nu_c = 15
779 else
780 nu_c = nint(nu_c_scale/nc(k)) + 2
781 rand = 0.0
782 if (present(rand2)) then
783 rand = rand2
784 endif
785 nu_c = max(2, min(nu_c+nint(rand), 15))
786 endif
787 xdc = max(d0c*1.e6, ((rc(k)/(am_r*nc(k)))**obmr) * 1.e6)
788 lamc = (nc(k)*am_r* ccg(2,nu_c) * ocg1(nu_c) / rc(k))**obmr
789 mvd_c(k) = (3.0+nu_c+0.672) / lamc
790 mvd_c(k) = max(d0c, min(mvd_c(k), d0r))
791 endif
792
793 !..Autoconversion follows Berry & Reinhardt (1974) with characteristic
794 !.. diameters correctly computed from gamma distrib of cloud droplets.
795 if (rc(k).gt. 0.01e-3) then
796 dc_g = ((ccg(3,nu_c)*ocg2(nu_c))**obmr / lamc) * 1.e6
797 dc_b = (xdc*xdc*xdc*dc_g*dc_g*dc_g - xdc*xdc*xdc*xdc*xdc*xdc) &
798 **(1./6.)
799 zeta1 = 0.5*((6.25e-6*xdc*dc_b*dc_b*dc_b - 0.4) &
800 + abs(6.25e-6*xdc*dc_b*dc_b*dc_b - 0.4))
801 zeta = 0.027*rc(k)*zeta1
802 taud = 0.5*((0.5*dc_b - 7.5) + abs(0.5*dc_b - 7.5)) + r1
803 tau = 3.72/(rc(k)*taud)
804 prr_wau(k) = zeta/tau
805 prr_wau(k) = min(real(rc(k)*odts, kind=dp), prr_wau(k))
806 pnr_wau(k) = prr_wau(k) / (am_r*nu_c*10.*d0r*d0r*d0r) ! RAIN2M
807 pnc_wau(k) = min(real(nc(k)*odts, kind=dp), prr_wau(k) / (am_r*mvd_c(k)*mvd_c(k)*mvd_c(k))) ! Qc2M
808 endif
809
811 if (l_qr(k) .and. mvd_r(k).gt. d0r .and. mvd_c(k).gt. d0c) then
812 lamr = 1./ilamr(k)
813 idx = 1 + int(nbr*log(real(mvd_r(k)/dr(1), kind=dp)) / log(real(dr(nbr)/dr(1), kind=dp)))
814 idx = min(idx, nbr)
815 ef_rw = t_efrw(idx, int(mvd_c(k)*1.e6))
816 prr_rcw(k) = rhof(k)*t1_qr_qc*ef_rw*rc(k)*n0_r(k) &
817 *((lamr+fv_r)**(-cre(9)))
818 prr_rcw(k) = min(real(rc(k)*odts, kind=dp), prr_rcw(k))
819 pnc_rcw(k) = rhof(k)*t1_qr_qc*ef_rw*nc(k)*n0_r(k) &
820 *((lamr+fv_r)**(-cre(9))) ! Qc2M
821 pnc_rcw(k) = min(real(nc(k)*odts, kind=dp), pnc_rcw(k))
822 endif
823
825 if (l_qr(k) .and. mvd_r(k).gt. d0r) then
826 ef_ra = eff_aero(mvd_r(k),0.04e-6,visco(k),rho(k),temp(k),'r')
827 lamr = 1./ilamr(k)
828 pna_rca(k) = rhof(k)*t1_qr_qc*ef_ra*nwfa(k)*n0_r(k) &
829 *((lamr+fv_r)**(-cre(9)))
830 pna_rca(k) = min(real(nwfa(k)*odts, kind=dp), pna_rca(k))
831
832 ef_ra = eff_aero(mvd_r(k),0.8e-6,visco(k),rho(k),temp(k),'r')
833 pnd_rcd(k) = rhof(k)*t1_qr_qc*ef_ra*nifa(k)*n0_r(k) &
834 *((lamr+fv_r)**(-cre(9)))
835 pnd_rcd(k) = min(real(nifa(k)*odts, kind=dp), pnd_rcd(k))
836 endif
837
838 enddo
839
840 !=================================================================================================================
841 !..Compute all frozen hydrometeor species' process terms.
842 !+---+-----------------------------------------------------------------+
843 if (.not. iiwarm) then
844 do k = kts, kte
845 vts_boost(k) = 1.0
846 xds = 0.0
847 if (l_qs(k)) xds = smoc(k) / smob(k)
848 ! orho = 1./rho(k)
849
850 ! if (L_qs(k)) then
851 ! xDs = smoc(k) / smob(k)
852 ! rho_s2 = max(rho_g(1), min(0.13/xDs, rho_i-100.))
853 ! else
854 ! xDs = 0.
855 ! rho_s2 = 100.
856 ! endif
857
858 !..Temperature lookup table indexes.
859 tempc = temp(k) - 273.15
860 idx_tc = max(1, min(nint(-tempc), 45) )
861 idx_t = int( (tempc-2.5)/5. ) - 1
862 idx_t = max(1, -idx_t)
863 idx_t = min(idx_t, ntb_t)
864 it = max(1, min(nint(-tempc), 31) )
865
867 if (rc(k).gt. r_c(1)) then
868 nic = nint(log10(rc(k)))
869 do_loop_rc: do nn = nic-1, nic+1
870 n = nn
871 if ( (rc(k)/10.**nn).ge.1.0 .and. (rc(k)/10.**nn).lt.10.0 ) exit do_loop_rc
872 enddo do_loop_rc
873 idx_c = int(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
874 idx_c = max(1, min(idx_c, ntb_c))
875 else
876 idx_c = 1
877 endif
878
880 idx_n = nint(1.0 + real(nbc, kind=wp) * log(real(nc(k)/t_nc(1), kind=dp)) / nic1)
881 idx_n = max(1, min(idx_n, nbc))
882
884 if (ri(k).gt. r_i(1)) then
885 nii = nint(log10(ri(k)))
886 do_loop_ri: do nn = nii-1, nii+1
887 n = nn
888 if ( (ri(k)/10.**nn).ge.1.0 .and. (ri(k)/10.**nn).lt.10.0 ) exit do_loop_ri
889 enddo do_loop_ri
890 idx_i = int(ri(k)/10.**n) + 10*(n-nii2) - (n-nii2)
891 idx_i = max(1, min(idx_i, ntb_i))
892 else
893 idx_i = 1
894 endif
895
896 if (ni(k).gt. nt_i(1)) then
897 nii = nint(log10(ni(k)))
898 do_loop_ni: do nn = nii-1, nii+1
899 n = nn
900 if ( (ni(k)/10.**nn).ge.1.0 .and. (ni(k)/10.**nn).lt.10.0 ) exit do_loop_ni
901 enddo do_loop_ni
902 idx_i1 = int(ni(k)/10.**n) + 10*(n-nii3) - (n-nii3)
903 idx_i1 = max(1, min(idx_i1, ntb_i1))
904 else
905 idx_i1 = 1
906 endif
907
909 if (rr(k).gt. r_r(1)) then
910 nir = nint(log10(rr(k)))
911 do_loop_rr: do nn = nir-1, nir+1
912 n = nn
913 if ( (rr(k)/10.**nn).ge.1.0 .and. (rr(k)/10.**nn).lt.10.0 ) exit do_loop_rr
914 enddo do_loop_rr
915 idx_r = int(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2)
916 idx_r = max(1, min(idx_r, ntb_r))
917
918 lamr = 1./ilamr(k)
919 lam_exp = lamr * (crg(3)*org2*org1)**bm_r
920 n0_exp = org1*rr(k)/am_r * lam_exp**cre(1)
921 nir = nint(log10(real(n0_exp, kind=dp)))
922 do_loop_nr: do nn = nir-1, nir+1
923 n = nn
924 if ( (n0_exp/10.**nn).ge.1.0 .and. (n0_exp/10.**nn).lt.10.0 ) exit do_loop_nr
925 enddo do_loop_nr
926 idx_r1 = int(n0_exp/10.**n) + 10*(n-nir3) - (n-nir3)
927 idx_r1 = max(1, min(idx_r1, ntb_r1))
928 else
929 idx_r = 1
930 idx_r1 = ntb_r1
931 endif
932
934 if (rs(k).gt. r_s(1)) then
935 nis = nint(log10(rs(k)))
936 do_loop_rs: do nn = nis-1, nis+1
937 n = nn
938 if ( (rs(k)/10.**nn).ge.1.0 .and. (rs(k)/10.**nn).lt.10.0 ) exit do_loop_rs
939 enddo do_loop_rs
940 idx_s = int(rs(k)/10.**n) + 10*(n-nis2) - (n-nis2)
941 idx_s = max(1, min(idx_s, ntb_s))
942 else
943 idx_s = 1
944 endif
945
947 if (rg(k).gt. r_g(1)) then
948 nig = nint(log10(rg(k)))
949 do_loop_rg: do nn = nig-1, nig+1
950 n = nn
951 if ( (rg(k)/10.**nn).ge.1.0 .and. (rg(k)/10.**nn).lt.10.0 ) exit do_loop_rg
952 enddo do_loop_rg
953 idx_g = int(rg(k)/10.**n) + 10*(n-nig2) - (n-nig2)
954 idx_g = max(1, min(idx_g, ntb_g))
955
956 lamg = 1./ilamg(k)
957 lam_exp = lamg * (cgg(3,1)*ogg2*ogg1)**bm_g
958 n0_exp = ogg1*rg(k)/am_g(idx_bg(k)) * lam_exp**cge(1,1)
959 nig = nint(log10(real(n0_exp, kind=dp)))
960 do_loop_ng: do nn = nig-1, nig+1
961 n = nn
962 if ( (n0_exp/10.**nn).ge.1.0 .and. (n0_exp/10.**nn).lt.10.0 ) exit do_loop_ng
963 enddo do_loop_ng
964 idx_g1 = int(n0_exp/10.**n) + 10*(n-nig3) - (n-nig3)
965 idx_g1 = max(1, min(idx_g1, ntb_g1))
966 else
967 idx_g = 1
968 idx_g1 = ntb_g1
969 endif
970
971 !..Deposition/sublimation prefactor (from Srivastava & Coen 1992).
972 otemp = 1./temp(k)
973 rvs = rho(k)*qvsi(k)
974 rvs_p = rvs*otemp*(lsub*otemp*orv - 1.)
975 rvs_pp = rvs * ( otemp*(lsub*otemp*orv - 1.) &
976 *otemp*(lsub*otemp*orv - 1.) &
977 + (-2.*lsub*otemp*otemp*otemp*orv) &
978 + otemp*otemp)
979 gamsc = lsub*diffu(k)/tcond(k) * rvs_p
980 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
981 * rvs_pp/rvs_p * rvs/rvs_p
982 alphsc = max(1.e-9, alphsc)
983 xsat = ssati(k)
984 if (abs(xsat).lt. 1.e-9) xsat=0.
985 t1_subl = 4.*pi*( 1.0 - alphsc*xsat &
986 + 2.*alphsc*alphsc*xsat*xsat &
987 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
988 / (1.+gamsc)
989
990 !..Snow collecting cloud water. In CE, assume Dc<<Ds and vtc=~0.
991 if (l_qc(k) .and. mvd_c(k).gt. d0c) then
992 ! xDs = 0.0
993 ! if (L_qs(k)) xDs = smoc(k) / smob(k)
994 if (xds > d0s) then
995 idx = 1 + int(nbs*log(real(xds/ds(1), kind=dp)) / log(real(ds(nbs)/ds(1), kind=dp)))
996 idx = min(idx, nbs)
997 ef_sw = t_efsw(idx, int(mvd_c(k)*1.e6))
998 prs_scw(k) = rhof(k)*t1_qs_qc*ef_sw*rc(k)*smoe(k)
999 prs_scw(k) = min(real(rc(k)*odts, kind=dp), prs_scw(k))
1000 pnc_scw(k) = rhof(k)*t1_qs_qc*ef_sw*nc(k)*smoe(k) ! Qc2M
1001 pnc_scw(k) = min(real(nc(k)*odts, kind=dp), pnc_scw(k))
1002 endif
1003
1004 !..Graupel collecting cloud water. In CE, assume Dc<<Dg and vtc=~0.
1005 if (rg(k).ge. r_g(1) .and. mvd_c(k).gt. d0c) then
1006 xdg = (bm_g + mu_g + 1.) * ilamg(k)
1007 vtg = rhof(k)*av_g(idx_bg(k))*cgg(6,idx_bg(k))*ogg3 * ilamg(k)**bv_g(idx_bg(k))
1008 stoke_g = mvd_c(k)*mvd_c(k)*vtg*rho_w2/(9.*visco(k)*xdg)
1009 !..Rime density formula of Cober and List (1993) also used by Milbrandt and Morrison (2014).
1010 const_ri = -1.*(mvd_c(k)*0.5e6)*vtg/min(-0.1,tempc)
1011 const_ri = max(0.1, min(const_ri, 10.))
1012 rime_dens = (0.051 + 0.114*const_ri - 0.0055*const_ri*const_ri)*1000.
1013 ! CCPP version has check on xDg > D0g
1014 if (xdg > d0g) then
1015 if (stoke_g.ge.0.4 .and. stoke_g.le.10.) then
1016 ef_gw = 0.55*log10(2.51*stoke_g)
1017 elseif (stoke_g.lt.0.4) then
1018 ef_gw = 0.0
1019 elseif (stoke_g.gt.10) then
1020 ef_gw = 0.77
1021 endif
1022 ! Not sure what to do here - hail increases size rapidly here below melting level.
1023 if (temp(k).gt.t_0) ef_gw = ef_gw*0.1
1024 t1_qg_qc = pi*.25*av_g(idx_bg(k)) * cgg(9,idx_bg(k))
1025 prg_gcw(k) = rhof(k)*t1_qg_qc*ef_gw*rc(k)*n0_g(k) &
1026 *ilamg(k)**cge(9,idx_bg(k))
1027 pnc_gcw(k) = rhof(k)*t1_qg_qc*ef_gw*nc(k)*n0_g(k) &
1028 *ilamg(k)**cge(9,idx_bg(k)) ! Qc2M
1029 pnc_gcw(k) = min(real(nc(k)*odts, kind=dp), pnc_gcw(k))
1030 if (temp(k).lt.t_0) pbg_gcw(k) = prg_gcw(k)/rime_dens
1031 ! CCPP version has end check on xDg > D0g
1032 endif
1033 endif
1034 endif
1035
1037 if (rs(k) .gt. r_s(1)) then
1038 ef_sa = eff_aero(xds,0.04e-6,visco(k),rho(k),temp(k),'s')
1039 pna_sca(k) = rhof(k)*t1_qs_qc*ef_sa*nwfa(k)*smoe(k)
1040 pna_sca(k) = min(real(nwfa(k)*odts, kind=dp), pna_sca(k))
1041
1042 ef_sa = eff_aero(xds,0.8e-6,visco(k),rho(k),temp(k),'s')
1043 pnd_scd(k) = rhof(k)*t1_qs_qc*ef_sa*nifa(k)*smoe(k)
1044 pnd_scd(k) = min(real(nifa(k)*odts, kind=dp), pnd_scd(k))
1045 endif
1046 if (rg(k) .gt. r_g(1)) then
1047 xdg = (bm_g + mu_g + 1.) * ilamg(k)
1048 ef_ga = eff_aero(xdg,0.04e-6,visco(k),rho(k),temp(k),'g')
1049 t1_qg_qc = pi*.25*av_g(idx_bg(k)) * cgg(9,idx_bg(k))
1050 pna_gca(k) = rhof(k)*t1_qg_qc*ef_ga*nwfa(k)*n0_g(k) &
1051 *ilamg(k)**cge(9,idx_bg(k))
1052 pna_gca(k) = min(real(nwfa(k)*odts, kind=dp), pna_gca(k))
1053
1054 ef_ga = eff_aero(xdg,0.8e-6,visco(k),rho(k),temp(k),'g')
1055 pnd_gcd(k) = rhof(k)*t1_qg_qc*ef_ga*nifa(k)*n0_g(k) &
1056 *ilamg(k)**cge(9,idx_bg(k))
1057 pnd_gcd(k) = min(real(nifa(k)*odts, kind=dp), pnd_gcd(k))
1058 endif
1059
1060 !..Rain collecting snow. Cannot assume Wisner (1972) approximation
1061 !.. or Mizuno (1990) approach so we solve the CE explicitly and store
1062 !.. results in lookup table.
1063 if (rr(k).ge. r_r(1)) then
1064 if (rs(k).ge. r_s(1)) then
1065 if (temp(k).lt.t_0) then
1066 prr_rcs(k) = -(tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
1067 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
1068 + tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
1069 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r))
1070 prs_rcs(k) = tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
1071 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) &
1072 - tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
1073 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
1074 prg_rcs(k) = tmr_racs1(idx_s,idx_t,idx_r1,idx_r) &
1075 + tcr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
1076 + tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
1077 + tms_sacr1(idx_s,idx_t,idx_r1,idx_r)
1078 prr_rcs(k) = max(real(-rr(k)*odts, kind=dp), prr_rcs(k))
1079 prs_rcs(k) = max(real(-rs(k)*odts, kind=dp), prs_rcs(k))
1080 prg_rcs(k) = min(real((rr(k)+rs(k))*odts, kind=dp), prg_rcs(k))
1081 pnr_rcs(k) = tnr_racs1(idx_s,idx_t,idx_r1,idx_r) & ! RAIN2M
1082 + tnr_racs2(idx_s,idx_t,idx_r1,idx_r) &
1083 + tnr_sacr1(idx_s,idx_t,idx_r1,idx_r) &
1084 + tnr_sacr2(idx_s,idx_t,idx_r1,idx_r)
1085 pnr_rcs(k) = min(real(nr(k)*odts, kind=dp), pnr_rcs(k))
1086 png_rcs(k) = pnr_rcs(k)
1087 !-GT pbg_rcs(k) = prg_rcs(k)/(0.5*(rho_i+rho_s))
1088 pbg_rcs(k) = prg_rcs(k)/rho_i
1089 else
1090 prs_rcs(k) = -tcs_racs1(idx_s,idx_t,idx_r1,idx_r) &
1091 - tms_sacr1(idx_s,idx_t,idx_r1,idx_r) &
1092 + tmr_racs2(idx_s,idx_t,idx_r1,idx_r) &
1093 + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r)
1094 prs_rcs(k) = max(real(-rs(k)*odts, kind=dp), prs_rcs(k))
1095 prr_rcs(k) = -prs_rcs(k)
1096 !pnr_rcs(k) = tnr_racs2(idx_s,idx_t,idx_r1,idx_r) & ! RAIN2M
1097 ! + tnr_sacr2(idx_s,idx_t,idx_r1,idx_r)
1098 endif
1099 !pnr_rcs(k) = MIN(DBLE(nr(k)*odts), pnr_rcs(k))
1100 endif
1101
1102 !..Rain collecting graupel. Cannot assume Wisner (1972) approximation
1103 !.. or Mizuno (1990) approach so we solve the CE explicitly and store
1104 !.. results in lookup table.
1105 if (rg(k).ge. r_g(1)) then
1106 if (temp(k).lt.t_0) then
1107 prg_rcg(k) = tmr_racg(idx_g1,idx_g,idx_table(k),idx_r1,idx_r) &
1108 + tcr_gacr(idx_g1,idx_g,idx_table(k),idx_r1,idx_r)
1109 prg_rcg(k) = min(real(rr(k)*odts, kind=dp), prg_rcg(k))
1110 prr_rcg(k) = -prg_rcg(k)
1111 pnr_rcg(k) = tnr_racg(idx_g1,idx_g,idx_table(k),idx_r1,idx_r) & ! RAIN2M
1112 + tnr_gacr(idx_g1,idx_g,idx_table(k),idx_r1,idx_r)
1113 pnr_rcg(k) = min(real(nr(k)*odts, kind=dp), pnr_rcg(k))
1114 !-GT pbg_rcg(k) = prg_rcg(k)/(0.5*(rho_i+rho_g(idx_bg(k))))
1115 pbg_rcg(k) = prg_rcg(k)/rho_i
1116 else
1117 prr_rcg(k) = tcg_racg(idx_g1,idx_g,idx_table(k),idx_r1,idx_r)
1118 prr_rcg(k) = min(real(rg(k)*odts, kind=dp), prr_rcg(k))
1119 prg_rcg(k) = -prr_rcg(k)
1120 png_rcg(k) = tnr_racg(idx_g1,idx_g,idx_table(k),idx_r1,idx_r)
1121!!! + tnr_gacr(idx_g1,idx_g,idx_table(k),idx_r1,idx_r)
1122 png_rcg(k) = min(real(ng(k)*odts, kind=dp), png_rcg(k))
1123 pbg_rcg(k) = prg_rcg(k)/rho_g(idx_bg(k))
1124 !..Put in explicit drop break-up due to collisions.
1125 pnr_rcg(k) = -1.5*tnr_gacr(idx_g1,idx_g,idx_table(k),idx_r1,idx_r) ! RAIN2M
1126 endif
1127 endif
1128 endif
1129
1130 !+---+-----------------------------------------------------------------+
1131 !..Next IF block handles only those processes below 0C.
1132 !+---+-----------------------------------------------------------------+
1133
1134 if (temp(k).lt.t_0) then
1135
1136 vts_boost(k) = 1.0
1137 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
1138
1139 !+---+---------------- BEGIN NEW ICE NUCLEATION -----------------------+
1140 !..Freezing of supercooled water (rain or cloud) is influenced by dust
1141 !.. but still using Bigg 1953 with a temperature adjustment of a few
1142 !.. degrees depending on dust concentration. A default value by way
1143 !.. of idx_IN is 1.0 per Liter of air is used when dustyIce flag is
1144 !.. false. Next, a combination of deposition/condensation freezing
1145 !.. using DeMott et al (2010) dust nucleation when water saturated or
1146 !.. Phillips et al (2008) when below water saturation; else, without
1147 !.. dustyIce flag, use the previous Cooper (1986) temperature-dependent
1148 !.. value. Lastly, allow homogeneous freezing of deliquesced aerosols
1149 !.. following Koop et al. (2001, Nature).
1150 !.. Implemented by T. Eidhammer and G. Thompson 2012Dec18
1151 !+---+-----------------------------------------------------------------+
1152
1153 ! if (dustyIce) then
1154 if (dustyice .and. (configs%aerosol_aware .or. merra2_aerosol_aware)) then
1155 xni = icedemott(tempc,qvs(k),qvs(k),qvsi(k),rho(k),nifa(k))
1156 else
1157 xni = 1.0 *1000. ! Default is 1.0 per Liter
1158 endif
1159
1161 if (xni.gt. nt_in(1)) then
1162 niin = nint(log10(xni))
1163 do_loop_xni: do nn = niin-1, niin+1
1164 n = nn
1165 if ( (xni/10.**nn).ge.1.0 .and. (xni/10.**nn).lt.10.0 ) exit do_loop_xni
1166 enddo do_loop_xni
1167 idx_in = int(xni/10.**n) + 10*(n-niin2) - (n-niin2)
1168 idx_in = max(1, min(idx_in, ntb_in))
1169 else
1170 idx_in = 1
1171 endif
1172
1173 !..Freezing of water drops into graupel/cloud ice (Bigg 1953).
1174 if (rr(k).gt. r_r(1)) then
1175 prg_rfz(k) = tpg_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
1176 pri_rfz(k) = tpi_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
1177 pni_rfz(k) = tni_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts
1178 pnr_rfz(k) = tnr_qrfz(idx_r,idx_r1,idx_tc,idx_in)*odts ! RAIN2M
1179 pnr_rfz(k) = min(real(nr(k)*odts, kind=dp), pnr_rfz(k))
1180 elseif (rr(k).gt. r1 .and. temp(k).lt.hgfr) then
1181 pri_rfz(k) = rr(k)*odts
1182 pni_rfz(k) = nr(k)*odts
1183 endif
1184 pbg_rfz(k) = prg_rfz(k)/rho_i
1185
1186 if (rc(k).gt. r_c(1)) then
1187 pri_wfz(k) = tpi_qcfz(idx_c,idx_n,idx_tc,idx_in)*odts
1188 pri_wfz(k) = min(real(rc(k)*odts, kind=dp), pri_wfz(k))
1189 pni_wfz(k) = tni_qcfz(idx_c,idx_n,idx_tc,idx_in)*odts
1190 pni_wfz(k) = min(real(nc(k)*odts, kind=dp), pri_wfz(k)/(2.0_dp*xm0i), pni_wfz(k))
1191 elseif (rc(k).gt. r1 .and. temp(k).lt.hgfr) then
1192 pri_wfz(k) = rc(k)*odts
1193 pni_wfz(k) = nc(k)*odts
1194 endif
1195
1196 !..Deposition nucleation of dust/mineral from DeMott et al (2010)
1197 !.. we may need to relax the temperature and ssati constraints.
1198 if ( (ssati(k).ge. demott_nuc_ssati) .or. (ssatw(k).gt. eps &
1199 .and. temp(k).lt.253.15) ) then
1200
1201 if (dustyice .and. (configs%aerosol_aware .or. merra2_aerosol_aware)) then
1202 xnc = icedemott(tempc,qv(k),qvs(k),qvsi(k),rho(k),nifa(k))
1203 rand = 0.0
1204 if (present(rand3)) then
1205 rand = rand3
1206 endif
1207 xnc = xnc*(1.0 + 50.*rand)
1208 else
1209 xnc = min(icenuc_max, tno*exp(ato*(t_0-temp(k))))
1210 endif
1211 xni = ni(k) + (pni_rfz(k)+pni_wfz(k))*dtsave
1212 pni_inu(k) = 0.5*(xnc-xni + abs(xnc-xni))*odts
1213 pri_inu(k) = min(real(rate_max, kind=dp), xm0i*pni_inu(k))
1214 pni_inu(k) = pri_inu(k)/xm0i
1215 endif
1216
1217 !..Freezing of aqueous aerosols based on Koop et al (2001, Nature)
1218 xni = smo0(k)+ni(k) + (pni_rfz(k)+pni_wfz(k)+pni_inu(k))*dtsave
1219 if ((configs%aerosol_aware .or. merra2_aerosol_aware) .and. homogice .and. &
1220 (xni.le.max_ni) .and.(temp(k).lt.238.).and.(ssati(k).ge.0.4)) then
1221
1222 xnc = icekoop(temp(k),qv(k),qvs(k),nwfa(k), dtsave)
1223 pni_iha(k) = xnc*odts
1224 pri_iha(k) = min(real(rate_max, kind=dp), xm0i*0.1*pni_iha(k))
1225 pni_iha(k) = pri_iha(k)/(xm0i*0.1)
1226 endif
1227 !+---+------------------ END NEW ICE NUCLEATION -----------------------+
1228
1229
1230 !..Deposition/sublimation of cloud ice (Srivastava & Coen 1992).
1231 if (l_qi(k)) then
1232 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
1233 ilami = 1./lami
1234 xdi = max(real(d0i, kind=dp), (bm_i + mu_i + 1.) * ilami)
1235 xmi = am_i*xdi**bm_i
1236 oxmi = 1./xmi
1237 pri_ide(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
1238 *oig1*cig(5)*ni(k)*ilami
1239
1240 if (pri_ide(k) .lt. 0.0) then
1241 pri_ide(k) = max(real(-ri(k)*odts, kind=dp), pri_ide(k), real(rate_max, kind=dp))
1242 pni_ide(k) = pri_ide(k)*oxmi
1243 pni_ide(k) = max(real(-ni(k)*odts, kind=dp), pni_ide(k))
1244 else
1245 pri_ide(k) = min(pri_ide(k), real(rate_max, kind=dp))
1246 prs_ide(k) = (1.0_dp-tpi_ide(idx_i,idx_i1))*pri_ide(k)
1247 pri_ide(k) = tpi_ide(idx_i,idx_i1)*pri_ide(k)
1248 endif
1249
1250 !..Some cloud ice needs to move into the snow category. Use lookup
1251 !.. table that resulted from explicit bin representation of distrib.
1252 if ( (idx_i.eq. ntb_i) .or. (xdi.gt. 5.0*d0s) ) then
1253 prs_iau(k) = ri(k)*.99*odts
1254 pni_iau(k) = ni(k)*.95*odts
1255 elseif (xdi.lt. 0.1*d0s) then
1256 prs_iau(k) = 0.
1257 pni_iau(k) = 0.
1258 else
1259 prs_iau(k) = tps_iaus(idx_i,idx_i1)*odts
1260 prs_iau(k) = min(real(ri(k)*.99*odts, kind=dp), prs_iau(k))
1261 pni_iau(k) = tni_iaus(idx_i,idx_i1)*odts
1262 pni_iau(k) = min(real(ni(k)*.95*odts, kind=dp), pni_iau(k))
1263 endif
1264 endif
1265
1266 !..Deposition/sublimation of snow/graupel follows Srivastava & Coen
1267 !.. (1992).
1268 if (l_qs(k)) then
1269 c_snow = c_sqrd + (tempc+1.5)*(c_cube-c_sqrd)/(-30.+1.5)
1270 c_snow = max(c_sqrd, min(c_snow, c_cube))
1271 prs_sde(k) = c_snow*t1_subl*diffu(k)*ssati(k)*rvs &
1272 * (t1_qs_sd*smo1(k) &
1273 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
1274 if (prs_sde(k).lt. 0.) then
1275 prs_sde(k) = max(real(-rs(k)*odts, kind=dp), prs_sde(k), real(rate_max, kind=dp))
1276 else
1277 prs_sde(k) = min(prs_sde(k), real(rate_max, kind=dp))
1278 endif
1279 endif
1280
1281 if (l_qg(k) .and. ssati(k).lt. -eps) then
1282 t2_qg_sd = 0.28*sc3*sqrt(av_g(idx_bg(k))) * cgg(11,idx_bg(k))
1283 prg_gde(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
1284 * n0_g(k) * (t1_qg_sd*ilamg(k)**cge(10,1) &
1285 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11,idx_bg(k)))
1286 if (prg_gde(k).lt. 0.) then
1287 prg_gde(k) = max(real(-rg(k)*odts, kind=dp), prg_gde(k), real(rate_max, kind=dp))
1288 png_gde(k) = prg_gde(k) * ng(k)/rg(k)
1289 else
1290 prg_gde(k) = min(prg_gde(k), real(rate_max, kind=dp))
1291 endif
1292 endif
1293
1294 !..Snow collecting cloud ice. In CE, assume Di<<Ds and vti=~0.
1295 if (l_qi(k)) then
1296 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
1297 ilami = 1./lami
1298 xdi = max(real(d0i, kind=dp), (bm_i + mu_i + 1.) * ilami)
1299 xmi = am_i*xdi**bm_i
1300 oxmi = 1./xmi
1301 if (rs(k).ge. r_s(1)) then
1302 prs_sci(k) = t1_qs_qi*rhof(k)*ef_si*ri(k)*smoe(k)
1303 pni_sci(k) = prs_sci(k) * oxmi
1304 endif
1305
1306 !..Rain collecting cloud ice. In CE, assume Di<<Dr and vti=~0.
1307 if (rr(k).ge. r_r(1) .and. mvd_r(k).gt. 4.*xdi) then
1308 lamr = 1./ilamr(k)
1309 pri_rci(k) = rhof(k)*t1_qr_qi*ef_ri*ri(k)*n0_r(k) &
1310 *((lamr+fv_r)**(-cre(9)))
1311 pnr_rci(k) = rhof(k)*t1_qr_qi*ef_ri*ni(k)*n0_r(k) & ! RAIN2M
1312 *((lamr+fv_r)**(-cre(9)))
1313 pnr_rci(k) = min(real(nr(k)*odts, kind=dp), pnr_rci(k))
1314 pni_rci(k) = pri_rci(k) * oxmi
1315 prr_rci(k) = rhof(k)*t2_qr_qi*ef_ri*ni(k)*n0_r(k) &
1316 *((lamr+fv_r)**(-cre(8)))
1317 prr_rci(k) = min(real(rr(k)*odts, kind=dp), prr_rci(k))
1318 prg_rci(k) = pri_rci(k) + prr_rci(k)
1319 pbg_rci(k) = prg_rci(k)/rho_i
1320 endif
1321 endif
1322
1323 !..Ice multiplication from rime-splinters (Hallet & Mossop 1974).
1324 if (prg_gcw(k).gt. eps .and. tempc.gt.-8.0) then
1325 tf = 0.
1326 if (tempc.ge.-5.0 .and. tempc.lt.-3.0) then
1327 tf = 0.5*(-3.0 - tempc)
1328 elseif (tempc.gt.-8.0 .and. tempc.lt.-5.0) then
1329 tf = 0.33333333*(8.0 + tempc)
1330 endif
1331 pni_ihm(k) = 3.5e8*tf*prg_gcw(k)
1332 pri_ihm(k) = xm0i*pni_ihm(k)
1333 prs_ihm(k) = prs_scw(k)/(prs_scw(k)+prg_gcw(k)) &
1334 * pri_ihm(k)
1335 prg_ihm(k) = prg_gcw(k)/(prs_scw(k)+prg_gcw(k)) &
1336 * pri_ihm(k)
1337 endif
1338
1339 !..A portion of rimed snow converts to graupel but some remains snow.
1340 !.. Interp from 15 to 95% as riming factor increases from 2.0 to 30.0
1341 !.. 0.028 came from (.95-.15)/(30.-2.). This remains ad-hoc and should
1342 !.. be revisited.
1343 if (prs_scw(k).gt.rime_threshold*prs_sde(k) .and. prs_sde(k).gt.eps) then
1344 r_frac = min(30.0_dp, prs_scw(k)/prs_sde(k))
1345 g_frac = min(rime_conversion, 0.15 + (r_frac-2.)*.028)
1346 vts_boost(k) = min(1.5, 1.1 + (r_frac-2.)*.016)
1347 prg_scw(k) = g_frac*prs_scw(k)
1348 png_scw(k) = prg_scw(k)*smo0(k)/rs(k)
1349 !..gt png_scw(k) = prg_scw(k)*ns(k)/rs(k)
1350 vts = av_s*xds**bv_s * exp(-fv_s*xds)
1351 const_ri = -1.*(mvd_c(k)*0.5e6)*vts/min(-0.1,tempc)
1352 const_ri = max(0.1, min(const_ri, 10.))
1353 rime_dens = (0.051 + 0.114*const_ri - 0.0055*const_ri*const_ri)*1000.
1354 if(rime_dens .lt. 150.) then ! Idea of A. Jensen
1355 g_frac = 0.
1356 prg_scw(k)=0.
1357 png_scw(k)=0.
1358 endif
1359 pbg_scw(k) = prg_scw(k)/(0.5*(rime_dens+rho_s2))
1360 prs_scw(k) = (1. - g_frac)*prs_scw(k)
1361 endif
1362 else
1363
1364 !..Melt snow and graupel and enhance from collisions with liquid.
1365 !.. We also need to sublimate snow and graupel if subsaturated.
1366 if (l_qs(k)) then
1367 prr_sml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delqvs(k)) &
1368 * (t1_qs_me*smo1(k) + t2_qs_me*rhof2(k)*vsc2(k)*smof(k))
1369 if (prr_sml(k) .gt. 0.) then
1370 prr_sml(k) = prr_sml(k) + 4218.*olfus*tempc &
1371 * (prr_rcs(k)+prs_scw(k))
1372 prr_sml(k) = min(real(rs(k)*odts, kind=dp), max(zerod0, prr_sml(k)))
1373
1374 pnr_sml(k) = smo0(k)/rs(k)*prr_sml(k) * 10.0**(-0.25*tempc) ! RAIN2M
1375 pnr_sml(k) = min(real(smo0(k)*odts, kind=dp), pnr_sml(k))
1376 else
1377 prr_sml(k) = 0.0
1378 pnr_sml(k) = 0.0
1379 if (ssati(k).lt. 0.) then
1380 prs_sde(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
1381 * (t1_qs_sd*smo1(k) &
1382 + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
1383 prs_sde(k) = max(real(-rs(k)*odts, kind=dp), prs_sde(k))
1384
1385 endif
1386 endif
1387 endif
1388
1389 if (l_qg(k)) then
1390 n0_melt = n0_g(k)
1391 if ((rg(k)*ng(k)) .lt. 1.e-4) then
1392 lamg = 1./ilamg(k)
1393 n0_melt = (1.e-4/rg(k))*ogg2*lamg**cge(2,1)
1394 endif
1395 t2_qg_me = pi*4.*c_cube*olfus * 0.28*sc3*sqrt(av_g(idx_bg(k))) * cgg(11,idx_bg(k))
1396 prr_gml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delqvs(k)) &
1397 * n0_melt*(t1_qg_me*ilamg(k)**cge(10,1) &
1398 + t2_qg_me*rhof2(k)*vsc2(k)*ilamg(k)**cge(11,idx_bg(k)))
1399 ! if (prr_gml(k) .gt. 0.) then
1400 ! prr_gml(k) = prr_gml(k) + 4218.*olfus*(twet(k)-T_0) &
1401 ! * (prr_rcg(k)+prg_gcw(k))
1402 ! endif
1403 prr_gml(k) = min(real(rg(k)*odts, kind=dp), max(0.d0, prr_gml(k)))
1404 ! pnr_gml(k) = N0_g(k)*cgg(2)*ilamg(k)**cge(2) / rg(k) & ! RAIN2M
1405 ! * prr_gml(k) * 10.0**(-0.5*tempc)
1406
1407 if (prr_gml(k) .gt. 0.0) then
1408 melt_f = max(0.05, min(prr_gml(k)*dt/rg(k),1.0))
1409 !..1000 is density water, 50 is lower limit (max ice density is 800)
1410 pbg_gml(k) = prr_gml(k) / max(min(melt_f*rho_g(idx_bg(k)),1000.),50.)
1411 !-GT pnr_gml(k) = prr_gml(k)*ng(k)/rg(k)
1412 pnr_gml(k) = prr_gml(k)*ng(k)/rg(k) * 10.0**(-0.33*(temp(k)-t_0))
1413 else
1414 prr_gml(k) = 0.0
1415 pnr_gml(k) = 0.0
1416 pbg_gml(k) = 0.0
1417 if (ssati(k).lt. 0.) then
1418 t2_qg_sd = 0.28*sc3*sqrt(av_g(idx_bg(k))) * cgg(11,idx_bg(k))
1419 prg_gde(k) = c_cube*t1_subl*diffu(k)*ssati(k)*rvs &
1420 * n0_g(k) * (t1_qg_sd*ilamg(k)**cge(10,1) &
1421 + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11,idx_bg(k)))
1422 prg_gde(k) = max(real(-rg(k)*odts, kind=dp), prg_gde(k))
1423 png_gde(k) = prg_gde(k) * ng(k)/rg(k)
1424 endif
1425 endif
1426 endif
1427
1428 !.. This change will be required if users run adaptive time step that
1429 !.. results in delta-t that is generally too long to allow cloud water
1430 !.. collection by snow/graupel above melting temperature.
1431 !.. Credit to Bjorn-Egil Nygaard for discovering.
1432 if (dt .gt. 120.) then
1433 prr_rcw(k)=prr_rcw(k)+prs_scw(k)+prg_gcw(k)
1434 prs_scw(k)=0.
1435 prg_gcw(k)=0.
1436 endif
1437
1438 endif
1439 if (.not. configs%hail_aware) idx_bg(k) = idx_bg1
1440 enddo
1441 endif
1442
1443 !=================================================================================================================
1444 !..Ensure we do not deplete more hydrometeor species than exists.
1445 !+---+-----------------------------------------------------------------+
1446 do k = kts, kte
1447
1448 !..If ice supersaturated, ensure sum of depos growth terms does not
1449 !.. deplete more vapor than possibly exists. If subsaturated, limit
1450 !.. sum of sublimation terms such that vapor does not reproduce ice
1451 !.. supersat again.
1452 sump = pri_inu(k) + pri_ide(k) + prs_ide(k) &
1453 + prs_sde(k) + prg_gde(k) + pri_iha(k)
1454 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999
1455 if ( (sump.gt. eps .and. sump.gt. rate_max) .or. &
1456 (sump.lt. -eps .and. sump.lt. rate_max) ) then
1457 ratio = rate_max/sump
1458 pri_inu(k) = pri_inu(k) * ratio
1459 pri_ide(k) = pri_ide(k) * ratio
1460 pni_ide(k) = pni_ide(k) * ratio
1461 prs_ide(k) = prs_ide(k) * ratio
1462 prs_sde(k) = prs_sde(k) * ratio
1463 prg_gde(k) = prg_gde(k) * ratio
1464 pri_iha(k) = pri_iha(k) * ratio
1465 endif
1466
1467 !..Cloud water conservation.
1468 sump = -prr_wau(k) - pri_wfz(k) - prr_rcw(k) &
1469 - prs_scw(k) - prg_scw(k) - prg_gcw(k)
1470 rate_max = -rc(k)*odts
1471 if (sump.lt. rate_max .and. l_qc(k)) then
1472 ratio = rate_max/sump
1473 prr_wau(k) = prr_wau(k) * ratio
1474 pri_wfz(k) = pri_wfz(k) * ratio
1475 prr_rcw(k) = prr_rcw(k) * ratio
1476 prs_scw(k) = prs_scw(k) * ratio
1477 prg_scw(k) = prg_scw(k) * ratio
1478 prg_gcw(k) = prg_gcw(k) * ratio
1479 endif
1480
1481 !..Cloud ice conservation.
1482 sump = pri_ide(k) - prs_iau(k) - prs_sci(k) &
1483 - pri_rci(k)
1484 rate_max = -ri(k)*odts
1485 if (sump.lt. rate_max .and. l_qi(k)) then
1486 ratio = rate_max/sump
1487 pri_ide(k) = pri_ide(k) * ratio
1488 prs_iau(k) = prs_iau(k) * ratio
1489 prs_sci(k) = prs_sci(k) * ratio
1490 pri_rci(k) = pri_rci(k) * ratio
1491 endif
1492
1493 !..Rain conservation.
1494 sump = -prg_rfz(k) - pri_rfz(k) - prr_rci(k) &
1495 + prr_rcs(k) + prr_rcg(k)
1496 rate_max = -rr(k)*odts
1497 if (sump.lt. rate_max .and. l_qr(k)) then
1498 ratio = rate_max/sump
1499 prg_rfz(k) = prg_rfz(k) * ratio
1500 pri_rfz(k) = pri_rfz(k) * ratio
1501 prr_rci(k) = prr_rci(k) * ratio
1502 prr_rcs(k) = prr_rcs(k) * ratio
1503 prr_rcg(k) = prr_rcg(k) * ratio
1504 endif
1505
1506 !..Snow conservation.
1507 sump = prs_sde(k) - prs_ihm(k) - prr_sml(k) &
1508 + prs_rcs(k)
1509 rate_max = -rs(k)*odts
1510 if (sump.lt. rate_max .and. l_qs(k)) then
1511 ratio = rate_max/sump
1512 prs_sde(k) = prs_sde(k) * ratio
1513 prs_ihm(k) = prs_ihm(k) * ratio
1514 prr_sml(k) = prr_sml(k) * ratio
1515 prs_rcs(k) = prs_rcs(k) * ratio
1516 endif
1517
1518 !..Graupel conservation.
1519 sump = prg_gde(k) - prg_ihm(k) - prr_gml(k) &
1520 + prg_rcg(k)
1521 rate_max = -rg(k)*odts
1522 if (sump.lt. rate_max .and. l_qg(k)) then
1523 ratio = rate_max/sump
1524 prg_gde(k) = prg_gde(k) * ratio
1525 prg_ihm(k) = prg_ihm(k) * ratio
1526 prr_gml(k) = prr_gml(k) * ratio
1527 prg_rcg(k) = prg_rcg(k) * ratio
1528 endif
1529
1530 !..Re-enforce proper mass conservation for subsequent elements in case
1531 !.. any of the above terms were altered. Thanks P. Blossey. 2009Sep28
1532 pri_ihm(k) = prs_ihm(k) + prg_ihm(k)
1533 ratio = min(abs(prr_rcg(k)), abs(prg_rcg(k)) )
1534 prr_rcg(k) = ratio * sign(1.0, sngl(prr_rcg(k)))
1535 prg_rcg(k) = -prr_rcg(k)
1536 if (temp(k).gt.t_0) then
1537 ratio = min(abs(prr_rcs(k)), abs(prs_rcs(k)) )
1538 prr_rcs(k) = ratio * sign(1.0, sngl(prr_rcs(k)))
1539 prs_rcs(k) = -prr_rcs(k)
1540 endif
1541
1542 enddo
1543
1544 !=================================================================================================================
1545 !..Calculate tendencies of all species but constrain the number of ice
1546 !.. to reasonable values.
1547 !+---+-----------------------------------------------------------------+
1548 do k = kts, kte
1549 orho = 1./rho(k)
1550 lfus2 = lsub - lvap(k)
1551
1552 !..Aerosol number tendency
1553 if (configs%aerosol_aware) then
1554 nwfaten(k) = nwfaten(k) - (pna_rca(k) + pna_sca(k) &
1555 + pna_gca(k) + pni_iha(k)) * orho
1556 nifaten(k) = nifaten(k) - (pnd_rcd(k) + pnd_scd(k) &
1557 + pnd_gcd(k)) * orho
1558 if (dustyice) then
1559 nifaten(k) = nifaten(k) - pni_inu(k)*orho
1560 else
1561 nifaten(k) = 0.
1562 endif
1563 endif
1564
1565 !..Water vapor tendency
1566 qvten(k) = qvten(k) + (-pri_inu(k) - pri_iha(k) - pri_ide(k) &
1567 - prs_ide(k) - prs_sde(k) - prg_gde(k)) &
1568 * orho
1569
1570 !..Cloud water tendency
1571 qcten(k) = qcten(k) + (-prr_wau(k) - pri_wfz(k) &
1572 - prr_rcw(k) - prs_scw(k) - prg_scw(k) &
1573 - prg_gcw(k)) &
1574 * orho
1575
1576 !..Cloud water number tendency
1577 ncten(k) = ncten(k) + (-pnc_wau(k) - pnc_rcw(k) &
1578 - pni_wfz(k) - pnc_scw(k) - pnc_gcw(k)) &
1579 * orho
1580
1581 !..Cloud water mass/number balance; keep mass-wt mean size between
1582 !.. 1 and 50 microns. Also no more than Nt_c_max drops total.
1583 xrc=max(r1, (qc1d(k) + qcten(k)*dtsave)*rho(k))
1584 xnc=max(2., (nc1d(k) + ncten(k)*dtsave)*rho(k))
1585 if (xrc .gt. r1) then
1586 if (xnc.gt.10000.e6) then
1587 nu_c = 2
1588 elseif (xnc.lt.100.) then
1589 nu_c = 15
1590 else
1591 nu_c = nint(nu_c_scale/xnc) + 2
1592 rand = 0.0
1593 if (present(rand2)) then
1594 rand = rand2
1595 endif
1596 nu_c = max(2, min(nu_c+nint(rand), 15))
1597 endif
1598 lamc = (xnc*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
1599 xdc = (bm_r + nu_c + 1.) / lamc
1600 if (xdc.lt. d0c) then
1601 lamc = cce(2,nu_c)/d0c
1602 xnc = ccg(1,nu_c)*ocg2(nu_c)*xrc/am_r*lamc**bm_r
1603 ncten(k) = (xnc-nc1d(k)*rho(k))*odts*orho
1604 elseif (xdc.gt. d0r*2.) then
1605 lamc = cce(2,nu_c)/(d0r*2.)
1606 xnc = ccg(1,nu_c)*ocg2(nu_c)*xrc/am_r*lamc**bm_r
1607 ncten(k) = (xnc-nc1d(k)*rho(k))*odts*orho
1608 endif
1609 else
1610 ncten(k) = -nc1d(k)*odts
1611 endif
1612 xnc=max(0., (nc1d(k) + ncten(k)*dtsave)*rho(k))
1613 if (xnc.gt.nt_c_max) &
1614 ncten(k) = (nt_c_max-nc1d(k)*rho(k))*odts*orho
1615
1616 !..Cloud ice mixing ratio tendency
1617 qiten(k) = qiten(k) + (pri_inu(k) + pri_iha(k) + pri_ihm(k) &
1618 + pri_wfz(k) + pri_rfz(k) + pri_ide(k) &
1619 - prs_iau(k) - prs_sci(k) - pri_rci(k)) &
1620 * orho
1621
1622 !..Cloud ice number tendency.
1623 niten(k) = niten(k) + (pni_inu(k) + pni_iha(k) + pni_ihm(k) &
1624 + pni_wfz(k) + pni_rfz(k) + pni_ide(k) &
1625 - pni_iau(k) - pni_sci(k) - pni_rci(k)) &
1626 * orho
1627
1628 !..Cloud ice mass/number balance; keep mass-wt mean size between
1629 !.. 5 and 300 microns. Also no more than 500 xtals per liter.
1630 xri=max(r1,(qi1d(k) + qiten(k)*dtsave)*rho(k))
1631 xni=max(r2,(ni1d(k) + niten(k)*dtsave)*rho(k))
1632 if (xri.gt. r1) then
1633 lami = (am_i*cig(2)*oig1*xni/xri)**obmi
1634 ilami = 1./lami
1635 xdi = (bm_i + mu_i + 1.) * ilami
1636 if (xdi.lt. 5.e-6) then
1637 lami = cie(2)/5.e-6
1638 xni = min(max_ni, cig(1)*oig2*xri/am_i*lami**bm_i)
1639 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
1640 elseif (xdi.gt. 300.e-6) then
1641 lami = cie(2)/300.e-6
1642 xni = cig(1)*oig2*xri/am_i*lami**bm_i
1643 niten(k) = (xni-ni1d(k)*rho(k))*odts*orho
1644 endif
1645 else
1646 niten(k) = -ni1d(k)*odts
1647 endif
1648 xni=max(0.,(ni1d(k) + niten(k)*dtsave)*rho(k))
1649 if (xni.gt.max_ni) &
1650 niten(k) = (max_ni-ni1d(k)*rho(k))*odts*orho
1651
1652 !..Rain tendency
1653 qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) &
1654 + prr_sml(k) + prr_gml(k) + prr_rcs(k) &
1655 + prr_rcg(k) - prg_rfz(k) &
1656 - pri_rfz(k) - prr_rci(k)) &
1657 * orho
1658
1659 !..Rain number tendency
1660 nrten(k) = nrten(k) + (pnr_wau(k) + pnr_sml(k) + pnr_gml(k) &
1661 - (pnr_rfz(k) + pnr_rcr(k) + pnr_rcg(k) &
1662 + pnr_rcs(k) + pnr_rci(k) + pni_rfz(k)) ) &
1663 * orho
1664
1665 !..Rain mass/number balance; keep median volume diameter between
1666 !.. 37 microns (D0r*0.75) and 2.5 mm.
1667 xrr=max(r1,(qr1d(k) + qrten(k)*dtsave)*rho(k))
1668 xnr=max(r2,(nr1d(k) + nrten(k)*dtsave)*rho(k))
1669 if (xrr.gt. r1) then
1670 lamr = (am_r*crg(3)*org2*xnr/xrr)**obmr
1671 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
1672 if (mvd_r(k) .gt. 2.5e-3) then
1673 mvd_r(k) = 2.5e-3
1674 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
1675 xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
1676 nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
1677 elseif (mvd_r(k) .lt. d0r*0.75) then
1678 mvd_r(k) = d0r*0.75
1679 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
1680 xnr = crg(2)*org3*xrr*lamr**bm_r / am_r
1681 nrten(k) = (xnr-nr1d(k)*rho(k))*odts*orho
1682 endif
1683 else
1684 qrten(k) = -qr1d(k)*odts
1685 nrten(k) = -nr1d(k)*odts
1686 endif
1687
1688 !..Snow tendency
1689 qsten(k) = qsten(k) + (prs_iau(k) + prs_sde(k) + prs_sci(k) &
1690 + prs_scw(k) + prs_rcs(k) + prs_ide(k) &
1691 - prs_ihm(k) - prr_sml(k)) &
1692 * orho
1693
1694 !..Graupel tendency
1695 qgten(k) = qgten(k) + (prg_scw(k) + prg_rfz(k) &
1696 + prg_gde(k) + prg_rcg(k) + prg_gcw(k) &
1697 + prg_rci(k) + prg_rcs(k) - prg_ihm(k) &
1698 - prr_gml(k)) &
1699 * orho
1700
1701 !..Graupel number tendency
1702 ngten(k) = ngten(k) + (png_scw(k) + pnr_rfz(k) - png_rcg(k) &
1703 + pnr_rci(k) + png_rcs(k) + png_gde(k) &
1704 - pnr_gml(k)) * orho
1705
1706 !..Graupel volume mixing ratio tendency
1707 qbten(k) = qbten(k) + (pbg_scw(k) + pbg_rfz(k) &
1708 + pbg_gcw(k) + pbg_rci(k) + pbg_rcs(k) &
1709 + pbg_rcg(k) + pbg_sml(k) - pbg_gml(k) &
1710 + (prg_gde(k) - prg_ihm(k)) /rho_g(idx_bg(k)) ) &
1711 * orho
1712
1713 !..Graupel mass/number balance; keep its median volume diameter between
1714 !.. 3.0 times minimum size (D0g) and 25 mm.
1715 xrg=max(r1,(qg1d(k) + qgten(k)*dtsave)*rho(k))
1716 xng=max(r2,(ng1d(k) + ngten(k)*dtsave)*rho(k))
1717 xrb=max(xrg/rho(k)/rho_g(nrhg),(qb1d(k) + qbten(k)*dtsave))
1718 xrb=min(xrg/rho(k)/rho_g(1), xrb)
1719
1720 if (xrg .gt. r1) then
1721 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*xng/xrg)**obmg
1722 mvd_g(k) = (3.0 + mu_g + 0.672) / lamg
1723
1724 if (mvd_g(k) .gt. 25.4e-3) then
1725 mvd_g(k) = 25.4e-3
1726 lamg = (3.0 + mu_g + 0.672) / mvd_g(k)
1727 xng = cgg(2,1)*ogg3*xrg*lamg**bm_g / am_g(idx_bg(k))
1728 ngten(k) = (xng-ng1d(k)*rho(k))*odts*orho
1729 elseif (mvd_g(k) .lt. d0r) then
1730 mvd_g(k) = d0r
1731 lamg = (3.0 + mu_g + 0.672) / mvd_g(k)
1732 xng = cgg(2,1)*ogg3*xrg*lamg**bm_g / am_g(idx_bg(k))
1733 ngten(k) = (xng-ng1d(k)*rho(k))*odts*orho
1734 endif
1735
1736 else
1737 qgten(k) = -qg1d(k)*odts
1738 ngten(k) = -ng1d(k)*odts
1739 qbten(k) = -qb1d(k)*odts
1740 endif
1741
1742 !..Temperature tendency
1743 if (temp(k).lt.t_0) then
1744 tten(k) = tten(k) &
1745 + ( lsub*ocp(k)*(pri_inu(k) + pri_ide(k) &
1746 + prs_ide(k) + prs_sde(k) &
1747 + prg_gde(k) + pri_iha(k)) &
1748 + lfus2*ocp(k)*(pri_wfz(k) + pri_rfz(k) &
1749 + prg_rfz(k) + prs_scw(k) &
1750 + prg_scw(k) + prg_gcw(k) &
1751 + prg_rcs(k) + prs_rcs(k) &
1752 + prr_rci(k) + prg_rcg(k)) &
1753 )*orho * (1-ifdry)
1754 else
1755 tten(k) = tten(k) &
1756 + ( lfus*ocp(k)*(-prr_sml(k) - prr_gml(k) &
1757 - prr_rcg(k) - prr_rcs(k)) &
1758 + lsub*ocp(k)*(prs_sde(k) + prg_gde(k)) &
1759 )*orho * (1-ifdry)
1760 endif
1761
1762 enddo
1763
1764 !=================================================================================================================
1765 !..Update variables for TAU+1 before condensation & sedimention.
1766 !+---+-----------------------------------------------------------------+
1767 do k = kts, kte
1768 temp(k) = t1d(k) + dt*tten(k)
1769 otemp = 1./temp(k)
1770 tempc = temp(k) - 273.15
1771 qv(k) = max(min_qv, qv1d(k) + dt*qvten(k))
1772 rho(k) = roverrv*pres(k)/(r*temp(k)*(qv(k)+roverrv))
1773 rhof(k) = sqrt(rho_not/rho(k))
1774 rhof2(k) = sqrt(rhof(k))
1775 qvs(k) = rslf(pres(k), temp(k))
1776 ssatw(k) = qv(k)/qvs(k) - 1.
1777 if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0
1778 diffu(k) = 2.11e-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
1779 if (tempc .ge. 0.0) then
1780 visco(k) = (1.718+0.0049*tempc)*1.0e-5
1781 else
1782 visco(k) = (1.718+0.0049*tempc-1.2e-5*tempc*tempc)*1.0e-5
1783 endif
1784 vsc2(k) = sqrt(rho(k)/visco(k))
1785 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
1786 tcond(k) = (5.69 + 0.0168*tempc)*1.0e-5 * 418.936
1787 ocp(k) = 1./(cp2*(1.+0.887*qv(k)))
1788 lvt2(k)=lvap(k)*lvap(k)*ocp(k)*orv*otemp*otemp
1789
1790 if (configs%aerosol_aware) then
1791 nwfa(k) = max(nwfa_default*rho(k), (nwfa1d(k) + nwfaten(k)*dt)*rho(k))
1792 endif
1793 enddo
1794
1795 do k = kts, kte
1796 if ((qc1d(k) + qcten(k)*dt) .gt. r1) then
1797 rc(k) = (qc1d(k) + qcten(k)*dt)*rho(k)
1798 nc(k) = max(2., min((nc1d(k)+ncten(k)*dt)*rho(k), nt_c_max))
1799 if (.not.(configs%aerosol_aware .or. merra2_aerosol_aware)) then
1800 nc(k) = nt_c
1801 if (present(lsml)) then
1802 if(lsml == 1) then
1803 nc(k) = nt_c_l
1804 else
1805 nc(k) = nt_c_o
1806 endif
1807 endif
1808 endif
1809 l_qc(k) = .true.
1810 else
1811 rc(k) = r1
1812 nc(k) = 2.
1813 l_qc(k) = .false.
1814 endif
1815
1816 if ((qi1d(k) + qiten(k)*dt) .gt. r1) then
1817 ri(k) = (qi1d(k) + qiten(k)*dt)*rho(k)
1818 ni(k) = max(r2, (ni1d(k) + niten(k)*dt)*rho(k))
1819 l_qi(k) = .true.
1820 else
1821 ri(k) = r1
1822 ni(k) = r2
1823 l_qi(k) = .false.
1824 endif
1825
1826 if ((qr1d(k) + qrten(k)*dt) .gt. r1) then
1827 rr(k) = (qr1d(k) + qrten(k)*dt)*rho(k)
1828 nr(k) = max(r2, (nr1d(k) + nrten(k)*dt)*rho(k))
1829 l_qr(k) = .true.
1830 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
1831 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
1832 if (mvd_r(k) .gt. 2.5e-3) then
1833 mvd_r(k) = 2.5e-3
1834 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
1835 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
1836 elseif (mvd_r(k) .lt. d0r*0.75) then
1837 mvd_r(k) = d0r*0.75
1838 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
1839 nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
1840 endif
1841 else
1842 rr(k) = r1
1843 nr(k) = r2
1844 l_qr(k) = .false.
1845 endif
1846
1847 if ((qs1d(k) + qsten(k)*dt) .gt. r1) then
1848 rs(k) = (qs1d(k) + qsten(k)*dt)*rho(k)
1849 l_qs(k) = .true.
1850 else
1851 rs(k) = r1
1852 l_qs(k) = .false.
1853 endif
1854 enddo
1855
1856 if (configs%hail_aware) then
1857 do k = kts, kte
1858 if ((qg1d(k) + qgten(k)*dt) .gt. r1) then
1859 l_qg(k) = .true.
1860 rg(k) = (qg1d(k) + qgten(k)*dt)*rho(k)
1861 ng(k) = max(r2, (ng1d(k) + ngten(k)*dt)*rho(k))
1862 rb(k) = max(rg(k)/rho(k)/rho_g(nrhg), qb1d(k) + qbten(k)*dt)
1863 rb(k) = min(rg(k)/rho(k)/rho_g(1), rb(k))
1864 idx_bg(k) = max(1,min(nint(rg(k)/rho(k)/rb(k) *0.01)+1,nrhg))
1865 else
1866 rg(k) = r1
1867 ng(k) = r2
1868 rb(k) = r1/rho(k)/rho_g(nrhg)
1869 idx_bg(k) = idx_bg1
1870 l_qg(k) = .false.
1871 endif
1872 enddo
1873 else
1874 do k = kte, kts, -1
1875 idx_bg(k) = idx_bg1
1876 enddo
1877 do k = kte, kts, -1
1878 if ((qg1d(k) + qgten(k)*dt) .gt. r1) then
1879 rg(k) = (qg1d(k) + qgten(k)*dt)*rho(k)
1880 ygra1 = log10(max(1.e-9, rg(k)))
1881 zans1 = 3.4 + 2./7.*(ygra1+8.)
1882 ! zans1 = max(2., min(zans1, 6.))
1883 n0_exp = max(gonv_min, min(10.0**(zans1), gonv_max))
1884 lam_exp = (n0_exp*am_g(idx_bg(k))*cgg(1,1)/rg(k))**oge1
1885 lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg
1886 ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k))
1887 rb(k) = rg(k)/rho(k)/rho_g(idx_bg(k))
1888 else
1889 rg(k) = r1
1890 ng(k) = r2
1891 rb(k) = r1/rho(k)/rho_g(nrhg)
1892 l_qg(k) = .false.
1893 endif
1894 enddo
1895 endif
1896
1897 !=================================================================================================================
1898 !..With tendency-updated mixing ratios, recalculate snow moments and
1899 !.. intercepts/slopes of graupel and rain.
1900 !+---+-----------------------------------------------------------------+
1901 if (.not. iiwarm) then
1902 do k = kts, kte
1903 smo2(k) = 0.
1904 smob(k) = 0.
1905 smoc(k) = 0.
1906 smod(k) = 0.
1907 enddo
1908 do k = kts, kte
1909 if (.not. l_qs(k)) cycle
1910 tc0 = min(-0.1, temp(k)-273.15)
1911 smob(k) = rs(k)*oams
1912
1913 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
1914 !.. then we must compute actual 2nd moment and use as reference.
1915 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
1916 smo2(k) = smob(k)
1917 else
1918 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
1919 + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
1920 + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
1921 + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
1922 + sa(10)*bm_s*bm_s*bm_s
1923 a_ = 10.0**loga_
1924 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
1925 + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
1926 + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
1927 + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
1928 + sb(10)*bm_s*bm_s*bm_s
1929 smo2(k) = (smob(k)/a_)**(1./b_)
1930 endif
1931
1932 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
1933 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
1934 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
1935 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
1936 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
1937 + sa(10)*cse(1)*cse(1)*cse(1)
1938 a_ = 10.0**loga_
1939 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
1940 + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
1941 + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
1942 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
1943 smoc(k) = a_ * smo2(k)**b_
1944
1945 !..Calculate bm_s+bv_s (th) moment. Useful for sedimentation.
1946 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(14) &
1947 + sa(4)*tc0*cse(14) + sa(5)*tc0*tc0 &
1948 + sa(6)*cse(14)*cse(14) + sa(7)*tc0*tc0*cse(14) &
1949 + sa(8)*tc0*cse(14)*cse(14) + sa(9)*tc0*tc0*tc0 &
1950 + sa(10)*cse(14)*cse(14)*cse(14)
1951 a_ = 10.0**loga_
1952 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(14) + sb(4)*tc0*cse(14) &
1953 + sb(5)*tc0*tc0 + sb(6)*cse(14)*cse(14) &
1954 + sb(7)*tc0*tc0*cse(14) + sb(8)*tc0*cse(14)*cse(14) &
1955 + sb(9)*tc0*tc0*tc0 + sb(10)*cse(14)*cse(14)*cse(14)
1956 smod(k) = a_ * smo2(k)**b_
1957 enddo
1958
1959 !+---+-----------------------------------------------------------------+
1960 !..Calculate y-intercept, slope values for graupel.
1961 !+---+-----------------------------------------------------------------+
1962 do k = kte, kts, -1
1963 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg
1964 ilamg(k) = 1./lamg
1965 n0_g(k) = ng(k)*ogg2*lamg**cge(2,1)
1966 enddo
1967
1968 endif
1969
1970 !+---+-----------------------------------------------------------------+
1971 !..Calculate y-intercept, slope values for rain.
1972 !+---+-----------------------------------------------------------------+
1973 do k = kte, kts, -1
1974 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
1975 ilamr(k) = 1./lamr
1976 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
1977 n0_r(k) = nr(k)*org2*lamr**cre(2)
1978 enddo
1979
1980 !=================================================================================================================
1981 !..Cloud water condensation and evaporation. Nucleate cloud droplets
1982 !.. using explicit CCN aerosols with hygroscopicity like sulfates using
1983 !.. parcel model lookup table results provided by T. Eidhammer. Evap
1984 !.. drops using calculation of max drop size capable of evaporating in
1985 !.. single timestep and explicit number of drops smaller than Dc_star
1986 !.. from lookup table.
1987 !+---+-----------------------------------------------------------------+
1988 do k = kts, kte
1989 orho = 1./rho(k)
1990 if ((ssatw(k).gt. eps) .or. (ssatw(k).lt. -eps .and. l_qc(k))) then
1991 clap = (qv(k)-qvs(k))/(1. + lvt2(k)*qvs(k))
1992 do n = 1, 3
1993 fcd = qvs(k)* exp(lvt2(k)*clap) - qv(k) + clap
1994 dfcd = qvs(k)*lvt2(k)* exp(lvt2(k)*clap) + 1.
1995 clap = clap - fcd/dfcd
1996 enddo
1997 xrc = rc(k) + clap*rho(k)
1998 xnc = 0.
1999 if (xrc > r1) then
2000 prw_vcd(k) = clap*odt
2001 !+---+-----------------------------------------------------------------+ ! DROPLET NUCLEATION
2002 if (clap .gt. eps) then
2003 if (configs%aerosol_aware .or. merra2_aerosol_aware) then
2004 rand = 0.0
2005 if (present(rand3)) then
2006 rand = rand3
2007 endif
2008 xnc = max(2., activ_ncloud(temp(k), w1d(k)+rand, nwfa(k), lsml))
2009 else
2010 xnc = nt_c
2011 if (present(lsml)) then
2012 if(lsml == 1) then
2013 xnc = nt_c_l
2014 else
2015 xnc = nt_c_o
2016 endif
2017 endif
2018 endif
2019 pnc_wcd(k) = 0.5*(xnc-nc(k) + abs(xnc-nc(k)))*odts*orho
2020
2021 !+---+-----------------------------------------------------------------+ ! EVAPORATION
2022 elseif (clap .lt. -eps .AND. ssatw(k).lt.-1.e-6 .and. configs%aerosol_aware) then
2023 tempc = temp(k) - 273.15
2024 otemp = 1./temp(k)
2025 rvs = rho(k)*qvs(k)
2026 rvs_p = rvs*otemp*(lvap(k)*otemp*orv - 1.)
2027 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*orv - 1.) &
2028 *otemp*(lvap(k)*otemp*orv - 1.) &
2029 + (-2.*lvap(k)*otemp*otemp*otemp*orv) &
2030 + otemp*otemp)
2031 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
2032 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
2033 * rvs_pp/rvs_p * rvs/rvs_p
2034 alphsc = max(1.e-9, alphsc)
2035 xsat = ssatw(k)
2036 if (abs(xsat).lt. 1.e-9) xsat=0.
2037 t1_evap = 2.*pi*( 1.0 - alphsc*xsat &
2038 + 2.*alphsc*alphsc*xsat*xsat &
2039 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
2040 / (1.+gamsc)
2041
2042 dc_star = dsqrt(-2.d0*dt * t1_evap/(2.*pi) &
2043 * 4.*diffu(k)*ssatw(k)*rvs/rho_w2)
2044 idx_d = max(1, min(int(1.e6*dc_star), nbc))
2045
2046 idx_n = nint(1.0 + float(nbc) * dlog(nc(k)/t_nc(1)) / nic1)
2047 idx_n = max(1, min(idx_n, nbc))
2048
2050 if (rc(k).gt. r_c(1)) then
2051 nic = nint(log10(rc(k)))
2052 do_loop_rc_cond: do nn = nic-1, nic+1
2053 n = nn
2054 if ( (rc(k)/10.**nn).ge.1.0 .and. (rc(k)/10.**nn).lt.10.0 ) exit do_loop_rc_cond
2055 enddo do_loop_rc_cond
2056 idx_c = int(rc(k)/10.**n) + 10*(n-nic2) - (n-nic2)
2057 idx_c = max(1, min(idx_c, ntb_c))
2058 else
2059 idx_c = 1
2060 endif
2061
2062 !prw_vcd(k) = MAX(DBLE(-rc(k)*orho*odt), &
2063 ! -tpc_wev(idx_d, idx_c, idx_n)*orho*odt)
2064 prw_vcd(k) = max(real(-rc(k)*0.99*orho*odt, kind=dp), prw_vcd(k))
2065 pnc_wcd(k) = max(real(-nc(k)*0.99*orho*odt, &
2066 kind=dp), real(-tnc_wev(idx_d, idx_c, idx_n)*orho*odt, kind=dp))
2067 endif
2068 else
2069 prw_vcd(k) = -rc(k)*orho*odt
2070 pnc_wcd(k) = -nc(k)*orho*odt
2071 endif
2072
2073 !+---+-----------------------------------------------------------------+
2074
2075 qvten(k) = qvten(k) - prw_vcd(k)
2076 qcten(k) = qcten(k) + prw_vcd(k)
2077 ncten(k) = ncten(k) + pnc_wcd(k)
2078 ! Be careful here: depending on initial conditions,
2079 ! cloud evaporation can increase aerosols
2080 if (configs%aerosol_aware) nwfaten(k) = nwfaten(k) - pnc_wcd(k)
2081
2082 tten(k) = tten(k) + lvap(k)*ocp(k)*prw_vcd(k)*(1-ifdry)
2083 rc(k) = max(r1, (qc1d(k) + dt*qcten(k))*rho(k))
2084 if (rc(k).eq.r1) l_qc(k) = .false.
2085 nc(k) = max(2., min((nc1d(k)+ncten(k)*dt)*rho(k), nt_c_max))
2086 if (.not.(configs%aerosol_aware .or. merra2_aerosol_aware)) then
2087 nc(k) = nt_c
2088 if (present(lsml)) then
2089 if(lsml == 1) then
2090 nc(k) = nt_c_l
2091 else
2092 nc(k) = nt_c_o
2093 endif
2094 endif
2095 endif
2096 qv(k) = max(min_qv, qv1d(k) + dt*qvten(k))
2097 temp(k) = t1d(k) + dt*tten(k)
2098 rho(k) = roverrv*pres(k)/(r*temp(k)*(qv(k)+roverrv))
2099 qvs(k) = rslf(pres(k), temp(k))
2100 ssatw(k) = qv(k)/qvs(k) - 1.
2101 endif
2102 enddo
2103
2104 !=================================================================================================================
2105 !.. If still subsaturated, allow rain to evaporate, following
2106 !.. Srivastava & Coen (1992).
2107 !+---+-----------------------------------------------------------------+
2108 do k = kts, kte
2109 if ( (ssatw(k).lt. -eps) .and. l_qr(k) &
2110 .and. (.not.(prw_vcd(k).gt. 0.)) ) then
2111 tempc = temp(k) - 273.15
2112 otemp = 1./temp(k)
2113 orho = 1./rho(k)
2114 rhof(k) = sqrt(rho_not*orho)
2115 rhof2(k) = sqrt(rhof(k))
2116 diffu(k) = 2.11e-5*(temp(k)/273.15)**1.94 * (101325./pres(k))
2117 if (tempc .ge. 0.0) then
2118 visco(k) = (1.718+0.0049*tempc)*1.0e-5
2119 else
2120 visco(k) = (1.718+0.0049*tempc-1.2e-5*tempc*tempc)*1.0e-5
2121 endif
2122 vsc2(k) = sqrt(rho(k)/visco(k))
2123 lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc
2124 tcond(k) = (5.69 + 0.0168*tempc)*1.0e-5 * 418.936
2125 ocp(k) = 1./(cp2*(1.+0.887*qv(k)))
2126
2127 rvs = rho(k)*qvs(k)
2128 rvs_p = rvs*otemp*(lvap(k)*otemp*orv - 1.)
2129 rvs_pp = rvs * ( otemp*(lvap(k)*otemp*orv - 1.) &
2130 *otemp*(lvap(k)*otemp*orv - 1.) &
2131 + (-2.*lvap(k)*otemp*otemp*otemp*orv) &
2132 + otemp*otemp)
2133 gamsc = lvap(k)*diffu(k)/tcond(k) * rvs_p
2134 alphsc = 0.5*(gamsc/(1.+gamsc))*(gamsc/(1.+gamsc)) &
2135 * rvs_pp/rvs_p * rvs/rvs_p
2136 alphsc = max(1.e-9, alphsc)
2137 xsat = min(-1.e-9, ssatw(k))
2138 t1_evap = 2.*pi*( 1.0 - alphsc*xsat &
2139 + 2.*alphsc*alphsc*xsat*xsat &
2140 - 5.*alphsc*alphsc*alphsc*xsat*xsat*xsat ) &
2141 / (1.+gamsc)
2142
2143 lamr = 1./ilamr(k)
2144 !..Rapidly eliminate near zero values when low humidity (<95%)
2145 if (qv(k)/qvs(k) .lt. 0.95 .AND. rr(k)*orho.le.1.e-8) then
2146 prv_rev(k) = rr(k)*orho*odts
2147 else
2148 prv_rev(k) = t1_evap*diffu(k)*(-ssatw(k))*n0_r(k)*rvs &
2149 * (t1_qr_ev*ilamr(k)**cre(10) &
2150 + t2_qr_ev*vsc2(k)*rhof2(k)*((lamr+0.5*fv_r)**(-cre(11))))
2151 rate_max = min((rr(k)*orho*odts), (qvs(k)-qv(k))*odts)
2152 prv_rev(k) = min(real(rate_max, kind=dp), prv_rev(k)*orho)
2153
2154 !..TEST: G. Thompson 10 May 2013
2155 !..Reduce the rain evaporation in same places as melting graupel occurs.
2156 !..Rationale: falling and simultaneous melting graupel in subsaturated
2157 !..regions will not melt as fast because particle temperature stays
2158 !..at 0C. Also not much shedding of the water from the graupel so
2159 !..likely that the water-coated graupel evaporating much slower than
2160 !..if the water was immediately shed off.
2161 IF (prr_gml(k).gt.0.0) THEN
2162 eva_factor = min(1.0, 0.01+(0.99-0.01)*(tempc/20.0))
2163 prv_rev(k) = prv_rev(k)*eva_factor
2164 ENDIF
2165 endif
2166
2167 pnr_rev(k) = min(real(nr(k)*0.99*orho*odts, kind=dp), & ! RAIN2M
2168 prv_rev(k) * nr(k)/rr(k))
2169
2170 qrten(k) = qrten(k) - prv_rev(k)
2171 qvten(k) = qvten(k) + prv_rev(k)
2172 nrten(k) = nrten(k) - pnr_rev(k)
2173 nwfaten(k) = nwfaten(k) + pnr_rev(k)
2174 tten(k) = tten(k) - lvap(k)*ocp(k)*prv_rev(k)*(1-ifdry)
2175
2176 rr(k) = max(r1, (qr1d(k) + dt*qrten(k))*rho(k))
2177 qv(k) = max(1.e-10, qv1d(k) + dt*qvten(k))
2178 nr(k) = max(r2, (nr1d(k) + dt*nrten(k))*rho(k))
2179 temp(k) = t1d(k) + dt*tten(k)
2180 rho(k) = roverrv*pres(k)/(r*temp(k)*(qv(k)+roverrv))
2181 endif
2182 enddo
2183#if defined(mpas)
2184 do k = kts, kte
2185 evapprod(k) = prv_rev(k) - (min(zerod0,prs_sde(k)) + &
2186 min(zerod0,prg_gde(k)))
2187 rainprod(k) = prr_wau(k) + prr_rcw(k) + prs_scw(k) + &
2188 prg_scw(k) + prs_iau(k) + &
2189 prg_gcw(k) + prs_sci(k) + &
2190 pri_rci(k)
2191 enddo
2192#endif
2193
2194 !=================================================================================================================
2195 !..Find max terminal fallspeed (distribution mass-weighted mean
2196 !.. velocity) and use it to determine if we need to split the timestep
2197 !.. (var nstep>1). Either way, only bother to do sedimentation below
2198 !.. 1st level that contains any sedimenting particles (k=ksed1 on down).
2199 !.. New in v3.0+ is computing separate for rain, ice, snow, and
2200 !.. graupel species thus making code faster with credit to J. Schmidt.
2201 !+---+-----------------------------------------------------------------+
2202 nstep = 0
2203 onstep(:) = 1.0
2204 ksed1(:) = 1
2205 do k = kte+1, kts, -1
2206 vtrk(k) = 0.
2207 vtnrk(k) = 0.
2208 vtik(k) = 0.
2209 vtnik(k) = 0.
2210 vtsk(k) = 0.
2211 vtgk(k) = 0.
2212 vtngk(k) = 0.
2213 vtck(k) = 0.
2214 vtnck(k) = 0.
2215 enddo
2216 if (any(l_qr .eqv. .true.)) then
2217 do k = kte, kts, -1
2218 vtr = 0.
2219 rhof(k) = sqrt(rho_not/rho(k))
2220
2221 if (rr(k).gt. r1) then
2222 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
2223 vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) &
2224 *((lamr+fv_r)**(-cre(6)))
2225 vtrk(k) = vtr
2226 ! First below is technically correct:
2227 ! vtr = rhof(k)*av_r*crg(5)*org2 * lamr**cre(2) &
2228 ! *((lamr+fv_r)**(-cre(5)))
2229 ! Test: make number fall faster (but still slower than mass)
2230 ! Goal: less prominent size sorting
2231 vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) &
2232 *((lamr+fv_r)**(-cre(7)))
2233 vtnrk(k) = vtr
2234 else
2235 vtrk(k) = vtrk(k+1)
2236 vtnrk(k) = vtnrk(k+1)
2237 endif
2238
2239 if (max(vtrk(k),vtnrk(k)) .gt. 1.e-3) then
2240 ksed1(1) = max(ksed1(1), k)
2241 delta_tp = dzq(k)/(max(vtrk(k),vtnrk(k)))
2242 nstep = max(nstep, int(dt/delta_tp + 1.))
2243 endif
2244 enddo
2245 if (ksed1(1) .eq. kte) ksed1(1) = kte-1
2246 if (nstep .gt. 0) onstep(1) = 1./real(nstep)
2247 endif
2248
2249 !+---+-----------------------------------------------------------------+
2250
2251 if (any(l_qc .eqv. .true.)) then
2252 hgt_agl = 0.
2253 do_loop_hgt_agl : do k = kts, kte-1
2254 if (rc(k) .gt. r2) ksed1(5) = k
2255 hgt_agl = hgt_agl + dzq(k)
2256 if (hgt_agl .gt. 500.0) exit do_loop_hgt_agl
2257 enddo do_loop_hgt_agl
2258
2259 do k = ksed1(5), kts, -1
2260 vtc = 0.
2261 if (rc(k) .gt. r1 .and. w1d(k) .lt. 1.e-1) then
2262 if (nc(k).gt.10000.e6) then
2263 nu_c = 2
2264 elseif (nc(k).lt.100.) then
2265 nu_c = 15
2266 else
2267 nu_c = nint(nu_c_scale/nc(k)) + 2
2268 rand = 0.0
2269 if (present(rand2)) then
2270 rand = rand2
2271 endif
2272 nu_c = max(2, min(nu_c+nint(rand), 15))
2273 endif
2274 lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr
2275 ilamc = 1./lamc
2276 vtc = rhof(k)*av_c*ccg(5,nu_c)*ocg2(nu_c) * ilamc**bv_c
2277 vtck(k) = vtc
2278 vtc = rhof(k)*av_c*ccg(4,nu_c)*ocg1(nu_c) * ilamc**bv_c
2279 vtnck(k) = vtc
2280 endif
2281 enddo
2282 endif
2283
2284 !+---+-----------------------------------------------------------------+
2285
2286 if (.not. iiwarm) then
2287 if (any(l_qi .eqv. .true.)) then
2288
2289 nstep = 0
2290 do k = kte, kts, -1
2291 vti = 0.
2292
2293 if (ri(k).gt. r1) then
2294 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
2295 ilami = 1./lami
2296 vti = rhof(k)*av_i*cig(3)*oig2 * ilami**bv_i
2297 vtik(k) = vti
2298 ! First below is technically correct:
2299 ! vti = rhof(k)*av_i*cig(4)*oig1 * ilami**bv_i
2300 ! Goal: less prominent size sorting
2301 vti = rhof(k)*av_i*cig(6)/cig(7) * ilami**bv_i
2302 vtnik(k) = vti
2303 else
2304 vtik(k) = vtik(k+1)
2305 vtnik(k) = vtnik(k+1)
2306 endif
2307
2308 if (vtik(k) .gt. 1.e-3) then
2309 ksed1(2) = max(ksed1(2), k)
2310 delta_tp = dzq(k)/vtik(k)
2311 nstep = max(nstep, int(dt/delta_tp + 1.))
2312 endif
2313 enddo
2314 if (ksed1(2) .eq. kte) ksed1(2) = kte-1
2315 if (nstep .gt. 0) onstep(2) = 1./real(nstep)
2316 endif
2317
2318 !+---+-----------------------------------------------------------------+
2319
2320 if (any(l_qs .eqv. .true.)) then
2321
2322 nstep = 0
2323 do k = kte, kts, -1
2324 vts = 0.
2325
2326 if (rs(k).gt. r1) then
2327 xds = smoc(k) / smob(k)
2328 mrat = 1./xds
2329 ils1 = 1./(mrat*lam0 + fv_s)
2330 ils2 = 1./(mrat*lam1 + fv_s)
2331 t1_vts = kap0*csg(4)*ils1**cse(4)
2332 t2_vts = kap1*mrat**mu_s*csg(10)*ils2**cse(10)
2333 ils1 = 1./(mrat*lam0)
2334 ils2 = 1./(mrat*lam1)
2335 t3_vts = kap0*csg(1)*ils1**cse(1)
2336 t4_vts = kap1*mrat**mu_s*csg(7)*ils2**cse(7)
2337 vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts)
2338 if (prr_sml(k) .gt. 0.0) then
2339 sr = rs(k)/(rs(k)+rr(k))
2340 vtsk(k) = vts*sr + (1.-sr)*vtrk(k)
2341 else
2342 vtsk(k) = vts*vts_boost(k)
2343 endif
2344 else
2345 vtsk(k) = vtsk(k+1)
2346 endif
2347
2348 if (vtsk(k) .gt. 1.e-3) then
2349 ksed1(3) = max(ksed1(3), k)
2350 delta_tp = dzq(k)/vtsk(k)
2351 nstep = max(nstep, int(dt/delta_tp + 1.))
2352 endif
2353 enddo
2354 if (ksed1(3) .eq. kte) ksed1(3) = kte-1
2355 if (nstep .gt. 0) onstep(3) = 1./real(nstep)
2356 endif
2357
2358 !+---+-----------------------------------------------------------------+
2359
2360 if (any(l_qg .eqv. .true.)) then
2361 nstep = 0
2362 do k = kte, kts, -1
2363 vtg = 0.
2364
2365 if (rg(k).gt. r1) then
2366 if (configs%hail_aware) then
2367 xrho_g = max(rho_g(1),min(rg(k)/rho(k)/rb(k),rho_g(nrhg)))
2368 afall = a_coeff*((4.0*xrho_g*9.8)/(3.0*rho(k)))**b_coeff
2369 afall = afall * visco(k)**(1.0-2.0*b_coeff)
2370 else
2371 afall = av_g_old
2372 bfall = bv_g_old
2373 endif
2374 vtg = rhof(k)*afall*cgg(6,idx_bg(k))*ogg3 * ilamg(k)**bfall
2375#if defined(ccpp_default)
2376 if (temp(k).gt. t_0) then
2377 vtgk(k) = max(vtg, vtrk(k))
2378 else
2379 vtgk(k) = vtg
2380 endif
2381#else
2382 vtgk(k) = vtg
2383#endif
2384 ! Goal: less prominent size sorting
2385 ! the ELSE section below is technically (mathematically) correct:
2386 if (mu_g .eq. 0) then
2387 vtg = rhof(k)*afall*cgg(7,idx_bg(k))/cgg(12,idx_bg(k)) * ilamg(k)**bfall
2388 else
2389 vtg = rhof(k)*afall*cgg(8,idx_bg(k))*ogg2 * ilamg(k)**bfall
2390 endif
2391 vtngk(k) = vtg
2392 else
2393 vtgk(k) = vtgk(k+1)
2394 vtngk(k) = vtngk(k+1)
2395 endif
2396
2397 if (vtgk(k) .gt. 1.e-3) then
2398 ksed1(4) = max(ksed1(4), k)
2399 delta_tp = dzq(k)/vtgk(k)
2400 nstep = max(nstep, int(dt/delta_tp + 1.))
2401 endif
2402 enddo
2403 if (ksed1(4) .eq. kte) ksed1(4) = kte-1
2404 if (nstep .gt. 0) onstep(4) = 1./real(nstep)
2405 endif
2406 endif
2407
2408 !=================================================================================================================
2409 !..Sedimentation of mixing ratio is the integral of v(D)*m(D)*N(D)*dD,
2410 !.. whereas neglect m(D) term for number concentration. Therefore,
2411 !.. cloud ice has proper differential sedimentation.
2412 !.. New in v3.0+ is computing separate for rain, ice, snow, and
2413 !.. graupel species thus making code faster with credit to J. Schmidt.
2414 !.. Bug fix, 2013Nov01 to tendencies using rho(k+1) correction thanks to
2415 !.. Eric Skyllingstad.
2416 !+---+-----------------------------------------------------------------+
2417
2418 if (any(l_qr .eqv. .true.)) then
2419 nstep = nint(1./onstep(1))
2420 if(.not. sedi_semi) then
2421 do n = 1, nstep
2422 do k = kte, kts, -1
2423 sed_r(k) = vtrk(k)*rr(k)
2424 sed_n(k) = vtnrk(k)*nr(k)
2425 enddo
2426 k = kte
2427 odzq = 1./dzq(k)
2428 orho = 1./rho(k)
2429 qrten(k) = qrten(k) - sed_r(k)*odzq*onstep(1)*orho
2430 nrten(k) = nrten(k) - sed_n(k)*odzq*onstep(1)*orho
2431 rr(k) = max(r1, rr(k) - sed_r(k)*odzq*dt*onstep(1))
2432 nr(k) = max(r2, nr(k) - sed_n(k)*odzq*dt*onstep(1))
2433#if defined(ccpp_default)
2434 pfll1(k) = pfll1(k) + sed_r(k)*dt*onstep(1)
2435#endif
2436 do k = ksed1(1), kts, -1
2437 odzq = 1./dzq(k)
2438 orho = 1./rho(k)
2439 qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) &
2440 *odzq*onstep(1)*orho
2441 nrten(k) = nrten(k) + (sed_n(k+1)-sed_n(k)) &
2442 *odzq*onstep(1)*orho
2443 rr(k) = max(r1, rr(k) + (sed_r(k+1)-sed_r(k)) &
2444 *odzq*dt*onstep(1))
2445 nr(k) = max(r2, nr(k) + (sed_n(k+1)-sed_n(k)) &
2446 *odzq*dt*onstep(1))
2447#if defined(ccpp_default)
2448 pfll1(k) = pfll1(k) + sed_r(k)*dt*onstep(1)
2449#endif
2450 enddo
2451
2452 if (rr(kts).gt.r1*1000.) &
2453 pptrain = pptrain + sed_r(kts)*dt*onstep(1)
2454 enddo
2455 else !if(.not. sedi_semi)
2456 niter = 1
2457 dtcfl = dt
2458 niter = int(nstep/max(decfl_,1)) + 1
2459 dtcfl = dt/niter
2460 do n = 1, niter
2461 rr_tmp(:) = rr(:)
2462 nr_tmp(:) = nr(:)
2463 call semi_lagrange_sedim(kte,dzq,vtrk,rr,rainsfc,pfll,dtcfl,r1)
2464 call semi_lagrange_sedim(kte,dzq,vtnrk,nr,vtr,pdummy,dtcfl,r2)
2465 do k = kts, kte
2466 orhodt = 1./(rho(k)*dt)
2467 qrten(k) = qrten(k) + (rr(k) - rr_tmp(k)) * orhodt
2468 nrten(k) = nrten(k) + (nr(k) - nr_tmp(k)) * orhodt
2469 pfll1(k) = pfll1(k) + pfll(k)
2470 enddo
2471 pptrain = pptrain + rainsfc
2472
2473 do k = kte+1, kts, -1
2474 vtrk(k) = 0.
2475 vtnrk(k) = 0.
2476 enddo
2477 do k = kte, kts, -1
2478 vtr = 0.
2479 if (rr(k).gt. r1) then
2480 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
2481 vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) &
2482 *((lamr+fv_r)**(-cre(6)))
2483 vtrk(k) = vtr
2484 ! First below is technically correct:
2485 ! vtr = rhof(k)*av_r*crg(5)*org2 * lamr**cre(2) &
2486 ! *((lamr+fv_r)**(-cre(5)))
2487 ! Test: make number fall faster (but still slower than mass)
2488 ! Goal: less prominent size sorting
2489 vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) &
2490 *((lamr+fv_r)**(-cre(7)))
2491 vtnrk(k) = vtr
2492 endif
2493 enddo
2494 enddo
2495 endif! if(.not. sedi_semi)
2496 endif
2497 !+---+-----------------------------------------------------------------+
2498
2499 if (any(l_qc .eqv. .true.)) then
2500
2501 do k = kte, kts, -1
2502 sed_c(k) = vtck(k)*rc(k)
2503 sed_n(k) = vtnck(k)*nc(k)
2504 enddo
2505 do k = ksed1(5), kts, -1
2506 odzq = 1./dzq(k)
2507 orho = 1./rho(k)
2508 qcten(k) = qcten(k) + (sed_c(k+1)-sed_c(k)) *odzq*orho
2509 ncten(k) = ncten(k) + (sed_n(k+1)-sed_n(k)) *odzq*orho
2510 rc(k) = max(r1, rc(k) + (sed_c(k+1)-sed_c(k)) *odzq*dt)
2511 nc(k) = max(10., nc(k) + (sed_n(k+1)-sed_n(k)) *odzq*dt)
2512 enddo
2513 endif
2514
2515 !+---+-----------------------------------------------------------------+
2516
2517 if (any(l_qi .eqv. .true.)) then
2518
2519 nstep = nint(1./onstep(2))
2520 do n = 1, nstep
2521 do k = kte, kts, -1
2522 sed_i(k) = vtik(k)*ri(k)
2523 sed_n(k) = vtnik(k)*ni(k)
2524 enddo
2525 k = kte
2526 odzq = 1./dzq(k)
2527 orho = 1./rho(k)
2528 qiten(k) = qiten(k) - sed_i(k)*odzq*onstep(2)*orho
2529 niten(k) = niten(k) - sed_n(k)*odzq*onstep(2)*orho
2530 ri(k) = max(r1, ri(k) - sed_i(k)*odzq*dt*onstep(2))
2531 ni(k) = max(r2, ni(k) - sed_n(k)*odzq*dt*onstep(2))
2532#if defined(ccpp_default)
2533 pfil1(k) = pfil1(k) + sed_i(k)*dt*onstep(2)
2534#endif
2535 do k = ksed1(2), kts, -1
2536 odzq = 1./dzq(k)
2537 orho = 1./rho(k)
2538 qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) &
2539 *odzq*onstep(2)*orho
2540 niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) &
2541 *odzq*onstep(2)*orho
2542 ri(k) = max(r1, ri(k) + (sed_i(k+1)-sed_i(k)) &
2543 *odzq*dt*onstep(2))
2544 ni(k) = max(r2, ni(k) + (sed_n(k+1)-sed_n(k)) &
2545 *odzq*dt*onstep(2))
2546#if defined(ccpp_default)
2547 pfil1(k) = pfil1(k) + sed_i(k)*dt*onstep(2)
2548#endif
2549 enddo
2550
2551 if (ri(kts).gt.r1*1000.) &
2552 pptice = pptice + sed_i(kts)*dt*onstep(2)
2553 enddo
2554 endif
2555
2556 !+---+-----------------------------------------------------------------+
2557
2558 if (any(l_qs .eqv. .true.)) then
2559
2560 nstep = nint(1./onstep(3))
2561 do n = 1, nstep
2562 do k = kte, kts, -1
2563 sed_s(k) = vtsk(k)*rs(k)
2564 enddo
2565 k = kte
2566 odzq = 1./dzq(k)
2567 orho = 1./rho(k)
2568 qsten(k) = qsten(k) - sed_s(k)*odzq*onstep(3)*orho
2569 rs(k) = max(r1, rs(k) - sed_s(k)*odzq*dt*onstep(3))
2570#if defined(ccpp_default)
2571 pfil1(k) = pfil1(k) + sed_s(k)*dt*onstep(3)
2572#endif
2573 do k = ksed1(3), kts, -1
2574 odzq = 1./dzq(k)
2575 orho = 1./rho(k)
2576 qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) &
2577 *odzq*onstep(3)*orho
2578 rs(k) = max(r1, rs(k) + (sed_s(k+1)-sed_s(k)) &
2579 *odzq*dt*onstep(3))
2580#if defined(ccpp_default)
2581 pfil1(k) = pfil1(k) + sed_s(k)*dt*onstep(3)
2582#endif
2583 enddo
2584
2585 if (rs(kts).gt.r1*1000.) &
2586 pptsnow = pptsnow + sed_s(kts)*dt*onstep(3)
2587 enddo
2588 endif
2589
2590 !+---+-----------------------------------------------------------------+
2591
2592 if (any(l_qg .eqv. .true.)) then
2593 nstep = nint(1./onstep(4))
2594 if(.not. sedi_semi) then
2595
2596 do n = 1, nstep
2597 do k = kte, kts, -1
2598 sed_g(k) = vtgk(k)*rg(k)
2599 sed_n(k) = vtngk(k)*ng(k)
2600 sed_b(k) = vtgk(k)*rb(k)
2601 enddo
2602 k = kte
2603 odzq = 1./dzq(k)
2604 orho = 1./rho(k)
2605 qgten(k) = qgten(k) - sed_g(k)*odzq*onstep(4)*orho
2606 ngten(k) = ngten(k) - sed_n(k)*odzq*onstep(4)*orho
2607 qbten(k) = qbten(k) - sed_b(k)*odzq*onstep(4)
2608 rg(k) = max(r1, rg(k) - sed_g(k)*odzq*dt*onstep(4))
2609 ng(k) = max(r2, ng(k) - sed_n(k)*odzq*dt*onstep(4))
2610 rb(k) = max(r1/rho(k)/rho_g(nrhg), rb(k) - sed_b(k)*odzq*dt*onstep(4))
2611#if defined(ccpp_default)
2612 pfil1(k) = pfil1(k) + sed_g(k)*dt*onstep(4)
2613#endif
2614 do k = ksed1(4), kts, -1
2615 odzq = 1./dzq(k)
2616 orho = 1./rho(k)
2617 qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) &
2618 *odzq*onstep(4)*orho
2619 ngten(k) = ngten(k) + (sed_n(k+1)-sed_n(k)) &
2620 *odzq*onstep(4)*orho
2621 qbten(k) = qbten(k) + (sed_b(k+1)-sed_b(k)) &
2622 *odzq*onstep(4)
2623 rg(k) = max(r1, rg(k) + (sed_g(k+1)-sed_g(k)) &
2624 *odzq*dt*onstep(4))
2625 ng(k) = max(r2, ng(k) + (sed_n(k+1)-sed_n(k)) &
2626 *odzq*dt*onstep(4))
2627 rb(k) = max(rg(k)/rho(k)/rho_g(nrhg), rb(k) + (sed_b(k+1)-sed_b(k)) &
2628 *odzq*dt*onstep(4))
2629#if defined(ccpp_default)
2630 pfil1(k) = pfil1(k) + sed_g(k)*dt*onstep(4)
2631#endif
2632 enddo
2633
2634 if (rg(kts).gt.r1*1000.) &
2635 pptgraul = pptgraul + sed_g(kts)*dt*onstep(4)
2636 enddo
2637 else ! if(.not. sedi_semi) then
2638 niter = 1
2639 dtcfl = dt
2640 niter = int(nstep/max(decfl_,1)) + 1
2641 dtcfl = dt/niter
2642
2643 do n = 1, niter
2644 rg_tmp(:) = rg(:)
2645 call semi_lagrange_sedim(kte,dzq,vtgk,rg,graulsfc,pfil,dtcfl,r1)
2646 do k = kts, kte
2647 orhodt = 1./(rho(k)*dt)
2648 qgten(k) = qgten(k) + (rg(k) - rg_tmp(k))*orhodt
2649 pfil1(k) = pfil1(k) + pfil(k)
2650 enddo
2651 pptgraul = pptgraul + graulsfc
2652 do k = kte+1, kts, -1
2653 vtgk(k) = 0.
2654 enddo
2655 do k = kte, kts, -1
2656 vtg = 0.
2657 if (rg(k).gt. r1) then
2658 ygra1 = alog10(max(1.e-9, rg(k)))
2659 rand = 0.0
2660 if (present(rand1)) then
2661 rand = rand1
2662 endif
2663
2664 zans1 = 3.4 + 2./7.*(ygra1+8.) + rand
2665 n0_exp = 10.**(zans1)
2666 n0_exp = max(gonv_min, min(n0_exp, gonv_max))
2667 lam_exp = (n0_exp*am_g(idx_bg(k))*cgg(1,1)/rg(k))**oge1
2668 lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg
2669
2670 vtg = rhof(k)*afall*cgg(6,idx_bg(k))*ogg3 * (1./lamg)**bfall
2671 if (temp(k).gt. t_0) then
2672 vtgk(k) = max(vtg, vtrk(k))
2673 else
2674 vtgk(k) = vtg
2675 endif
2676 endif
2677 enddo
2678 enddo
2679 endif ! if(.not. sedi_semi) then
2680 endif
2681
2682 !+---+-----------------------------------------------------------------+
2683 !.. Instantly melt any cloud ice into cloud water if above 0C and
2684 !.. instantly freeze any cloud water found below HGFR.
2685 !+---+-----------------------------------------------------------------+
2686 if (.not. iiwarm) then
2687 do k = kts, kte
2688 xri = max(0.0, qi1d(k) + qiten(k)*dt)
2689 if ( (temp(k).gt. t_0) .and. (xri.gt. 0.0) ) then
2690 qcten(k) = qcten(k) + xri*odt
2691 ncten(k) = ncten(k) + ni1d(k)*odt
2692 qiten(k) = qiten(k) - xri*odt
2693 niten(k) = -ni1d(k)*odt
2694 tten(k) = tten(k) - lfus*ocp(k)*xri*odt*(1-ifdry)
2695 endif
2696
2697 xrc = max(0.0, qc1d(k) + qcten(k)*dt)
2698 if ( (temp(k).lt. hgfr) .and. (xrc.gt. 0.0) ) then
2699 lfus2 = lsub - lvap(k)
2700 xnc = nc1d(k) + ncten(k)*dt
2701 qiten(k) = qiten(k) + xrc*odt
2702 niten(k) = niten(k) + xnc*odt
2703 qcten(k) = qcten(k) - xrc*odt
2704 ncten(k) = ncten(k) - xnc*odt
2705 tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-ifdry)
2706 endif
2707 enddo
2708 endif
2709
2710 !=================================================================================================================
2711 !.. All tendencies computed, apply and pass back final values to parent.
2712 !+---+-----------------------------------------------------------------+
2713 do k = kts, kte
2714 t1d(k) = t1d(k) + tten(k)*dt
2715 qv1d(k) = max(min_qv, qv1d(k) + qvten(k)*dt)
2716 qc1d(k) = qc1d(k) + qcten(k)*dt
2717 nc1d(k) = max(2./rho(k), min(nc1d(k) + ncten(k)*dt, nt_c_max))
2718 if (configs%aerosol_aware) then
2719 nwfa1d(k) = max(nwfa_default, min(aero_max, (nwfa1d(k)+nwfaten(k)*dt)))
2720 nifa1d(k) = max(nifa_default, min(aero_max, (nifa1d(k)+nifaten(k)*dt)))
2721 endif
2722 if (qc1d(k) .le. r1) then
2723 qc1d(k) = 0.0
2724 nc1d(k) = 0.0
2725 else
2726 if (nc1d(k)*rho(k).gt.10000.e6) then
2727 nu_c = 2
2728 elseif (nc1d(k)*rho(k).lt.100.) then
2729 nu_c = 15
2730 else
2731 nu_c = nint(nu_c_scale/(nc1d(k)*rho(k))) + 2
2732 rand = 0.0
2733 if (present(rand2)) then
2734 rand = rand2
2735 endif
2736 nu_c = max(2, min(nu_c+nint(rand), 15))
2737 endif
2738
2739 lamc = (am_r*ccg(2,nu_c)*ocg1(nu_c)*nc1d(k)/qc1d(k))**obmr
2740 xdc = (bm_r + nu_c + 1.) / lamc
2741 if (xdc.lt. d0c) then
2742 lamc = cce(2,nu_c)/d0c
2743 elseif (xdc.gt. d0r*2.) then
2744 lamc = cce(2,nu_c)/(d0r*2.)
2745 endif
2746 nc1d(k) = min(ccg(1,nu_c)*ocg2(nu_c)*qc1d(k)/am_r*lamc**bm_r,&
2747 real(Nt_c_max, kind=dp)/rho(k))
2748 endif
2749
2750 qi1d(k) = qi1d(k) + qiten(k)*dt
2751 ni1d(k) = max(r2/rho(k), ni1d(k) + niten(k)*dt)
2752 if (qi1d(k) .le. r1) then
2753 qi1d(k) = 0.0
2754 ni1d(k) = 0.0
2755 else
2756 lami = (am_i*cig(2)*oig1*ni1d(k)/qi1d(k))**obmi
2757 ilami = 1./lami
2758 xdi = (bm_i + mu_i + 1.) * ilami
2759 if (xdi.lt. 5.e-6) then
2760 lami = cie(2)/5.e-6
2761 elseif (xdi.gt. 300.e-6) then
2762 lami = cie(2)/300.e-6
2763 endif
2764 ni1d(k) = min(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, max_ni/rho(k))
2765 endif
2766 qr1d(k) = qr1d(k) + qrten(k)*dt
2767 nr1d(k) = max(r2/rho(k), nr1d(k) + nrten(k)*dt)
2768 if (qr1d(k) .le. r1) then
2769 qr1d(k) = 0.0
2770 nr1d(k) = 0.0
2771 else
2772 lamr = (am_r*crg(3)*org2*nr1d(k)/qr1d(k))**obmr
2773 mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
2774 if (mvd_r(k) .gt. 2.5e-3) then
2775 mvd_r(k) = 2.5e-3
2776 elseif (mvd_r(k) .lt. d0r*0.75) then
2777 mvd_r(k) = d0r*0.75
2778 endif
2779 lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
2780 nr1d(k) = crg(2)*org3*qr1d(k)*lamr**bm_r / am_r
2781 endif
2782 qs1d(k) = qs1d(k) + qsten(k)*dt
2783 if (qs1d(k) .le. r1) qs1d(k) = 0.0
2784 qg1d(k) = qg1d(k) + qgten(k)*dt
2785 ng1d(k) = max(r2/rho(k), ng1d(k) + ngten(k)*dt)
2786 if (qg1d(k) .le. r1) then
2787 qg1d(k) = 0.0
2788 ng1d(k) = 0.0
2789 qb1d(k) = 0.0
2790 else
2791 qb1d(k) = max(qg1d(k)/rho_g(nrhg), qb1d(k) + qbten(k)*dt)
2792 qb1d(k) = min(qg1d(k)/rho_g(1), qb1d(k))
2793 idx_bg(k) = max(1,min(nint(qg1d(k)/qb1d(k) *0.01)+1,nrhg))
2794 if (.not. configs%hail_aware) idx_bg(k) = idx_bg1
2795 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng1d(k)/qg1d(k))**obmg
2796 mvd_g(k) = (3.0 + mu_g + 0.672) / lamg
2797 if (mvd_g(k) .gt. 25.4e-3) then
2798 mvd_g(k) = 25.4e-3
2799 elseif (mvd_g(k) .lt. d0r) then
2800 mvd_g(k) = d0r
2801 endif
2802 lamg = (3.0 + mu_g + 0.672) / mvd_g(k)
2803 ng1d(k) = cgg(2,1)*ogg3*qg1d(k)*lamg**bm_g / am_g(idx_bg(k))
2804 endif
2805
2806 enddo
2807
2808#if defined(ccpp_default)
2809 ! Diagnostics
2810 calculate_extended_diagnostics: if (ext_diag) then
2811 do k = kts, kte
2812 if(prw_vcd(k).gt.0)then
2813 prw_vcdc1(k) = prw_vcd(k)*dt
2814 elseif(prw_vcd(k).lt.0)then
2815 prw_vcde1(k) = -1*prw_vcd(k)*dt
2816 endif
2817 !heating/cooling diagnostics
2818 tpri_inu1(k) = pri_inu(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
2819
2820 if(pri_ide(k).gt.0)then
2821 tpri_ide1_d(k) = pri_ide(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
2822 else
2823 tpri_ide1_s(k) = -pri_ide(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
2824 endif
2825
2826 if(temp(k).lt.t_0)then
2827 tprs_ide1(k) = prs_ide(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
2828 endif
2829
2830 if(prs_sde(k).gt.0)then
2831 tprs_sde1_d(k) = prs_sde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
2832 else
2833 tprs_sde1_s(k) = -prs_sde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
2834 endif
2835
2836 if(prg_gde(k).gt.0)then
2837 tprg_gde1_d(k) = prg_gde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
2838 else
2839 tprg_gde1_s(k) = -prg_gde(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
2840 endif
2841
2842 tpri_iha1(k) = pri_iha(k)*lsub*ocp(k)*orho * (1-ifdry)*dt
2843 tpri_wfz1(k) = pri_wfz(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
2844 tpri_rfz1(k) = pri_rfz(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
2845 tprg_rfz1(k) = prg_rfz(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
2846 tprs_scw1(k) = prs_scw(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
2847 tprg_scw1(k) = prg_scw(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
2848 tprg_rcs1(k) = prg_rcs(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
2849
2850 if(temp(k).lt.t_0)then
2851 tprs_rcs1(k) = prs_rcs(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
2852 endif
2853
2854 tprr_rci1(k) = prr_rci(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
2855
2856 if(temp(k).lt.t_0)then
2857 tprg_rcg1(k) = prg_rcg(k)*lfus2*ocp(k)*orho * (1-ifdry)*dt
2858 endif
2859
2860 if(prw_vcd(k).gt.0)then
2861 tprw_vcd1_c(k) = lvap(k)*ocp(k)*prw_vcd(k)*(1-ifdry)*dt
2862 else
2863 tprw_vcd1_e(k) = -lvap(k)*ocp(k)*prw_vcd(k)*(1-ifdry)*dt
2864 endif
2865
2866 ! cooling terms
2867 tprr_sml1(k) = prr_sml(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
2868 tprr_gml1(k) = prr_gml(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
2869
2870 if(temp(k).ge.t_0)then
2871 tprr_rcg1(k) = -prr_rcg(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
2872 endif
2873
2874 if(temp(k).ge.t_0)then
2875 tprr_rcs1(k) = -prr_rcs(k)*lfus*ocp(k)*orho * (1-ifdry)*dt
2876 endif
2877
2878 tprv_rev1(k) = lvap(k)*ocp(k)*prv_rev(k)*(1-ifdry)*dt
2879 tten1(k) = tten(k)*dt
2880 qvten1(k) = qvten(k)*dt
2881 qiten1(k) = qiten(k)*dt
2882 qrten1(k) = qrten(k)*dt
2883 qsten1(k) = qsten(k)*dt
2884 qgten1(k) = qgten(k)*dt
2885 niten1(k) = niten(k)*dt
2886 nrten1(k) = nrten(k)*dt
2887 ncten1(k) = ncten(k)*dt
2888 qcten1(k) = qcten(k)*dt
2889 enddo
2890 endif calculate_extended_diagnostics
2891#endif
2892
2893 end subroutine mp_tempo_main
2894 !=================================================================================================================
2895 !..Function to compute collision efficiency of collector species (rain,
2896 !.. snow, graupel) of aerosols. Follows Wang et al, 2010, ACP, which
2897 !.. follows Slinn (1983).
2898 !+---+-----------------------------------------------------------------+
2899
2900 real function eff_aero(d, da, visc,rhoa,temp,species)
2901
2902 implicit none
2903 real(wp) :: d, da, visc, rhoa, temp
2904 character(LEN=1) :: species
2905 real(wp) :: aval, cc, diff, re, sc, st, st2, vt, eff
2906 real(wp), parameter :: boltzman = 1.3806503e-23
2907 real(wp), parameter :: meanpath = 0.0256e-6
2908
2909 vt = 1.
2910 if (species .eq. 'r') then
2911 vt = -0.1021 + 4.932e3*d - 0.9551e6*d*d &
2912 + 0.07934e9*d*d*d - 0.002362e12*d*d*d*d
2913 elseif (species .eq. 's') then
2914 vt = av_s*d**bv_s
2915 elseif (species .eq. 'g') then
2916 vt = av_g(idx_bg1)*d**bv_g(idx_bg1)
2917 endif
2918
2919 cc = 1. + 2.*meanpath/da *(1.257+0.4*exp(-0.55*da/meanpath))
2920 diff = boltzman*temp*cc/(3.*pi*visc*da)
2921
2922 re = 0.5*rhoa*d*vt/visc
2923 sc = visc/(rhoa*diff)
2924
2925 st = da*da*vt*1000./(9.*visc*d)
2926 aval = 1.+log(1.+re)
2927 st2 = (1.2 + 1./12.*aval)/(1.+aval)
2928
2929 eff = 4./(re*sc) * (1. + 0.4*sqrt(re)*sc**0.3333 &
2930 + 0.16*sqrt(re)*sqrt(sc)) &
2931 + 4.*da/d * (0.02 + da/d*(1.+2.*sqrt(re)))
2932
2933 if (st.gt.st2) eff = eff + ( (st-st2)/(st-st2+0.666667))**1.5
2934 eff_aero = max(1.e-5, min(eff, 1.0))
2935
2936 end function eff_aero
2937
2938 !=================================================================================================================
2939 !..Retrieve fraction of CCN that gets activated given the model temp,
2940 !.. vertical velocity, and available CCN concentration. The lookup
2941 !.. table (read from external file) has CCN concentration varying the
2942 !.. quickest, then updraft, then temperature, then mean aerosol radius,
2943 !.. and finally hygroscopicity, kappa.
2944 !.. TO_DO ITEM: For radiation cooling producing fog, in which case the
2945 !.. updraft velocity could easily be negative, we could use the temp
2946 !.. and its tendency to diagnose a pretend postive updraft velocity.
2947 !+---+-----------------------------------------------------------------+
2948 real function activ_ncloud(tt, ww, nccn,lsm_in)
2949
2950 implicit none
2951 REAL, INTENT(IN):: tt, ww, nccn
2952 INTEGER, INTENT(IN), optional:: lsm_in
2953
2954 REAL:: n_local, w_local
2955 INTEGER:: i, j, k, l, m, n
2956 REAL:: a, b, c, d, t, u, x1, x2, y1, y2, nx, wy, fraction
2957 REAL:: lower_lim_nuc_frac
2958
2959
2960 ! ta_Na = (/10.0, 31.6, 100.0, 316.0, 1000.0, 3160.0, 10000.0/) ntb_arc
2961 ! ta_Ww = (/0.01, 0.0316, 0.1, 0.316, 1.0, 3.16, 10.0, 31.6, 100.0/) ntb_arw
2962 ! ta_Tk = (/243.15, 253.15, 263.15, 273.15, 283.15, 293.15, 303.15/) ntb_art
2963 ! ta_Ra = (/0.01, 0.02, 0.04, 0.08, 0.16/) ntb_arr
2964 ! ta_Ka = (/0.2, 0.4, 0.6, 0.8/) ntb_ark
2965
2966 n_local = nccn * 1.e-6
2967 w_local = ww
2968
2969 if (n_local .ge. ta_na(ntb_arc)) then
2970 n_local = ta_na(ntb_arc) - 1.0
2971 elseif (n_local .le. ta_na(1)) then
2972 n_local = ta_na(1) + 1.0
2973 endif
2974 do n = 2, ntb_arc
2975 if (n_local.ge.ta_na(n-1) .and. n_local.lt.ta_na(n)) goto 8003
2976 enddo
29778003 continue
2978 i = n
2979 x1 = log(ta_na(i-1))
2980 x2 = log(ta_na(i))
2981
2982 if (w_local .ge. ta_ww(ntb_arw)) then
2983 w_local = ta_ww(ntb_arw) - 1.0
2984 elseif (w_local .le. ta_ww(1)) then
2985 w_local = ta_ww(1) + 0.001
2986 endif
2987 do n = 2, ntb_arw
2988 if (w_local.ge.ta_ww(n-1) .and. w_local.lt.ta_ww(n)) goto 8005
2989 enddo
29908005 continue
2991 j = n
2992 y1 = log(ta_ww(j-1))
2993 y2 = log(ta_ww(j))
2994
2995 k = max(1, min( nint( (tt - ta_tk(1))*0.1) + 1, ntb_art))
2996
2997 !..The next two values are indexes of mean aerosol radius and
2998 !.. hygroscopicity. Currently these are constant but a future version
2999 !.. should implement other variables to allow more freedom such as
3000 !.. at least simple separation of tiny size sulfates from larger
3001 !.. sea salts.
3002 l = 3
3003 m = 2
3004
3005 lower_lim_nuc_frac = 0.
3006 if (present(lsm_in)) then
3007 if (lsm_in .eq. 1) then ! land
3008 lower_lim_nuc_frac = 0.
3009 else if (lsm_in .eq. 0) then ! water
3010 lower_lim_nuc_frac = 0.15
3011 else
3012 lower_lim_nuc_frac = 0.15 ! catch-all for anything else
3013 endif
3014 endif
3015
3016 a = tnccn_act(i-1,j-1,k,l,m)
3017 b = tnccn_act(i,j-1,k,l,m)
3018 c = tnccn_act(i,j,k,l,m)
3019 d = tnccn_act(i-1,j,k,l,m)
3020 nx = log(n_local)
3021 wy = log(w_local)
3022
3023 t = (nx-x1)/(x2-x1)
3024 u = (wy-y1)/(y2-y1)
3025
3026 ! t = (n_local-ta(Na(i-1))/(ta_Na(i)-ta_Na(i-1))
3027 ! u = (w_local-ta_Ww(j-1))/(ta_Ww(j)-ta_Ww(j-1))
3028
3029 fraction = (1.0-t)*(1.0-u)*a + t*(1.0-u)*b + t*u*c + (1.0-t)*u*d
3030 fraction = max(fraction, lower_lim_nuc_frac)
3031
3032 ! if (NCCN*fraction .gt. 0.75*Nt_c_max) then
3033 ! write(*,*) ' DEBUG-GT ', n_local, w_local, Tt, i, j, k
3034 ! endif
3035
3036 activ_ncloud = nccn*fraction
3037
3038 end function activ_ncloud
3039
3040 !=================================================================================================================
3041
3042 !+---+-----------------------------------------------------------------+
3043 real function icedemott(tempc, qv, qvs, qvsi, rho, nifa)
3044 implicit none
3045
3046 REAL, INTENT(IN):: tempc, qv, qvs, qvsi, rho, nifa
3047
3048 !..Local vars
3049 REAL:: satw, sati, siw, p_x, si0x, dtt, dsi, dsw, dab, fc, hx
3050 REAL:: ntilde, n_in, nmax, nhat, mux, xni, nifa_cc
3051 REAL, PARAMETER:: p_c1 = 1000.
3052 REAL, PARAMETER:: p_rho_c = 0.76
3053 REAL, PARAMETER:: p_alpha = 1.0
3054 REAL, PARAMETER:: p_gam = 2.
3055 REAL, PARAMETER:: delt = 5.
3056 REAL, PARAMETER:: t0x = -40.
3057 REAL, PARAMETER:: sw0x = 0.97
3058 REAL, PARAMETER:: delsi = 0.1
3059 REAL, PARAMETER:: hdm = 0.15
3060 REAL, PARAMETER:: p_psi = 0.058707*p_gam/p_rho_c
3061 REAL, PARAMETER:: aap = 1.
3062 REAL, PARAMETER:: bbp = 0.
3063 REAL, PARAMETER:: y1p = -35.
3064 REAL, PARAMETER:: y2p = -25.
3065 REAL, PARAMETER:: rho_not0 = 101325./(287.05*273.15)
3066
3067 !+---+
3068
3069 xni = 0.0
3070 ! satw = qv/qvs
3071 ! sati = qv/qvsi
3072 ! siw = qvs/qvsi
3073 ! p_x = -1.0261+(3.1656e-3*tempc)+(5.3938e-4*(tempc*tempc)) &
3074 ! + (8.2584e-6*(tempc*tempc*tempc))
3075 ! si0x = 1.+(10.**p_x)
3076 ! if (sati.ge.si0x .and. satw.lt.0.985) then
3077 ! dtt = delta_p (tempc, T0x, T0x+delT, 1., hdm)
3078 ! dsi = delta_p (sati, Si0x, Si0x+delSi, 0., 1.)
3079 ! dsw = delta_p (satw, Sw0x, 1., 0., 1.)
3080 ! fc = dtt*dsi*0.5
3081 ! hx = min(fc+((1.-fc)*dsw), 1.)
3082 ! ntilde = p_c1*p_gam*((exp(12.96*(sati-1.1)))**0.3) / p_rho_c
3083 ! if (tempc .le. y1p) then
3084 ! n_in = ntilde
3085 ! elseif (tempc .ge. y2p) then
3086 ! n_in = p_psi*p_c1*exp(12.96*(sati-1.)-0.639)
3087 ! else
3088 ! if (tempc .le. -30.) then
3089 ! nmax = p_c1*p_gam*(exp(12.96*(siw-1.1)))**0.3/p_rho_c
3090 ! else
3091 ! nmax = p_psi*p_c1*exp(12.96*(siw-1.)-0.639)
3092 ! endif
3093 ! ntilde = MIN(ntilde, nmax)
3094 ! nhat = MIN(p_psi*p_c1*exp(12.96*(sati-1.)-0.639), nmax)
3095 ! dab = delta_p (tempc, y1p, y2p, aap, bbp)
3096 ! n_in = MIN(nhat*(ntilde/nhat)**dab, nmax)
3097 ! endif
3098 ! mux = hx*p_alpha*n_in*rho
3099 ! xni = mux*((6700.*nifa)-200.)/((6700.*5.E5)-200.)
3100 ! elseif (satw.ge.0.985 .and. tempc.gt.HGFR-273.15) then
3101 nifa_cc = max(0.5, nifa*rho_not0*1.e-6/rho)
3102 ! xni = 3.*nifa_cc**(1.25)*exp((0.46*(-tempc))-11.6) ! [DeMott, 2015]
3103 xni = (5.94e-5*(-tempc)**3.33) & ! [DeMott, 2010]
3104 * (nifa_cc**((-0.0264*(tempc))+0.0033))
3105 xni = xni*rho/rho_not0 * 1000.
3106 ! endif
3107
3108 icedemott = max(0., xni)
3109
3110 end FUNCTION icedemott
3111
3112 !+---+-----------------------------------------------------------------+
3113 !..Newer research since Koop et al (2001) suggests that the freezing
3114 !.. rate should be lower than original paper, so J_rate is reduced
3115 !.. by two orders of magnitude.
3116
3117 real function icekoop(temp, qv, qvs, naero, dt)
3118 implicit none
3119
3120 REAL, INTENT(IN):: temp, qv, qvs, naero, dt
3121 REAL:: mu_diff, a_w_i, delta_aw, log_j_rate, j_rate, prob_h, satw
3122 REAL:: xni
3123
3124 xni = 0.0
3125 satw = qv/qvs
3126 mu_diff = 210368.0 + (131.438*temp) - (3.32373e6/temp) &
3127 & - (41729.1*alog(temp))
3128 a_w_i = exp(mu_diff/(r_uni*temp))
3129 delta_aw = satw - a_w_i
3130 log_j_rate = -906.7 + (8502.0*delta_aw) &
3131 & - (26924.0*delta_aw*delta_aw) &
3132 & + (29180.0*delta_aw*delta_aw*delta_aw)
3133 log_j_rate = min(20.0, log_j_rate)
3134 j_rate = 10.**log_j_rate ! cm-3 s-1
3135 prob_h = min(1.-exp(-j_rate*ar_volume*dt), 1.)
3136 if (prob_h .gt. 0.) then
3137 xni = min(prob_h*naero, 1000.e3)
3138 endif
3139
3140 icekoop = max(0.0, xni)
3141
3142 end FUNCTION icekoop
3143
3144 !+---+-----------------------------------------------------------------+
3145 !.. Helper routine for Phillips et al (2008) ice nucleation. Trude
3146
3147 REAL function delta_p (yy, y1, y2, aa, bb)
3148 IMPLICIT NONE
3149
3150 REAL, INTENT(IN):: yy, y1, y2, aa, bb
3151 REAL:: dab, a, b, a0, a1, a2, a3
3152
3153 a = 6.*(aa-bb)/((y2-y1)*(y2-y1)*(y2-y1))
3154 b = aa+(a*y1*y1*y1/6.)-(a*y1*y1*y2*0.5)
3155 a0 = b
3156 a1 = a*y1*y2
3157 a2 = -a*(y1+y2)*0.5
3158 a3 = a/3.
3159
3160 if (yy.le.y1) then
3161 dab = aa
3162 else if (yy.ge.y2) then
3163 dab = bb
3164 else
3165 dab = a0+(a1*yy)+(a2*yy*yy)+(a3*yy*yy*yy)
3166 endif
3167
3168 if (dab.lt.aa) then
3169 dab = aa
3170 endif
3171 if (dab.gt.bb) then
3172 dab = bb
3173 endif
3174 delta_p = dab
3175
3176 END FUNCTION delta_p
3177
3178 !+---+-----------------------------------------------------------------+
3179!-------------------------------------------------------------------
3180 SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1)
3181!-------------------------------------------------------------------
3182!
3183! This routine is a semi-Lagrangain forward advection for hydrometeors
3184! with mass conservation and positive definite advection
3185! 2nd order interpolation with monotonic piecewise parabolic method is used.
3186! This routine is under assumption of decfl < 1 for semi_Lagrangian
3187!
3188! dzl depth of model layer in meter
3189! wwl terminal velocity at model layer m/s
3190! rql dry air density*mixing ratio
3191! precip precipitation at surface
3192! dt time step
3193!
3194! author: hann-ming henry juang <henry.juang@noaa.gov>
3195! implemented by song-you hong
3196! reference: Juang, H.-M., and S.-Y. Hong, 2010: Forward semi-Lagrangian advection
3197! with mass conservation and positive definiteness for falling
3198! hydrometeors. *Mon. Wea. Rev.*, *138*, 1778-1791
3199!
3200 implicit none
3201
3202 integer, intent(in) :: km
3203 real, intent(in) :: dt, R1
3204 real, intent(in) :: dzl(km),wwl(km)
3205 real, intent(out) :: precip
3206 real, intent(inout) :: rql(km)
3207 real, intent(out) :: pfsan(km)
3208 integer k,m,kk,kb,kt
3209 real tl,tl2,qql,dql,qqd
3210 real th,th2,qqh,dqh
3211 real zsum,qsum,dim,dip,con1,fa1,fa2
3212 real allold, decfl
3213 real dz(km), ww(km), qq(km)
3214 real wi(km+1), zi(km+1), za(km+2)
3215 real qn(km)
3216 real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
3217 real net_flx(km)
3218!
3219 precip = 0.0
3220 qa(:) = 0.0
3221 qq(:) = 0.0
3222 dz(:) = dzl(:)
3223 ww(:) = wwl(:)
3224 do k = 1,km
3225 if(rql(k).gt.r1) then
3226 qq(k) = rql(k)
3227 else
3228 ww(k) = 0.0
3229 endif
3230 pfsan(k) = 0.0
3231 net_flx(k) = 0.0
3232 enddo
3233! skip for no precipitation for all layers
3234 allold = 0.0
3235 do k=1,km
3236 allold = allold + qq(k)
3237 enddo
3238 if(allold.le.0.0) then
3239 return
3240 endif
3241!
3242! compute interface values
3243 zi(1)=0.0
3244 do k=1,km
3245 zi(k+1) = zi(k)+dz(k)
3246 enddo
3247! n=1
3248! plm is 2nd order, we can use 2nd order wi or 3rd order wi
3249! 2nd order interpolation to get wi
3250 wi(1) = ww(1)
3251 wi(km+1) = ww(km)
3252 do k=2,km
3253 wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
3254 enddo
3255! 3rd order interpolation to get wi
3256 fa1 = 9./16.
3257 fa2 = 1./16.
3258 wi(1) = ww(1)
3259 wi(2) = 0.5*(ww(2)+ww(1))
3260 do k=3,km-1
3261 wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
3262 enddo
3263 wi(km) = 0.5*(ww(km)+ww(km-1))
3264 wi(km+1) = ww(km)
3265
3266! terminate of top of raingroup
3267 do k=2,km
3268 if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
3269 enddo
3270
3271! diffusivity of wi
3272 con1 = 0.05
3273 do k=km,1,-1
3274 decfl = (wi(k+1)-wi(k))*dt/dz(k)
3275 if( decfl .gt. con1 ) then
3276 wi(k) = wi(k+1) - con1*dz(k)/dt
3277 endif
3278 enddo
3279! compute arrival point
3280 do k=1,km+1
3281 za(k) = zi(k) - wi(k)*dt
3282 enddo
3283 za(km+2) = zi(km+1)
3284
3285 do k=1,km+1
3286 dza(k) = za(k+1)-za(k)
3287 enddo
3288
3289! computer deformation at arrival point
3290 do k=1,km
3291 qa(k) = qq(k)*dz(k)/dza(k)
3292 enddo
3293 qa(km+1) = 0.0
3294
3295! estimate values at arrival cell interface with monotone
3296 do k=2,km
3297 dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
3298 dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
3299 if( dip*dim.le.0.0 ) then
3300 qmi(k)=qa(k)
3301 qpi(k)=qa(k)
3302 else
3303 qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
3304 qmi(k)=2.0*qa(k)-qpi(k)
3305 if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
3306 qpi(k) = qa(k)
3307 qmi(k) = qa(k)
3308 endif
3309 endif
3310 enddo
3311 qpi(1)=qa(1)
3312 qmi(1)=qa(1)
3313 qmi(km+1)=qa(km+1)
3314 qpi(km+1)=qa(km+1)
3315
3316! interpolation to regular point
3317 qn = 0.0
3318 kb=1
3319 kt=1
3320 intp : do k=1,km
3321 kb=max(kb-1,1)
3322 kt=max(kt-1,1)
3323! find kb and kt
3324 if( zi(k).ge.za(km+1) ) then
3325 exit intp
3326 else
3327 find_kb : do kk=kb,km
3328 if( zi(k).le.za(kk+1) ) then
3329 kb = kk
3330 exit find_kb
3331 else
3332 cycle find_kb
3333 endif
3334 enddo find_kb
3335 find_kt : do kk=kt,km+2
3336 if( zi(k+1).le.za(kk) ) then
3337 kt = kk
3338 exit find_kt
3339 else
3340 cycle find_kt
3341 endif
3342 enddo find_kt
3343 kt = kt - 1
3344! compute q with piecewise constant method
3345 if( kt.eq.kb ) then
3346 tl=(zi(k)-za(kb))/dza(kb)
3347 th=(zi(k+1)-za(kb))/dza(kb)
3348 tl2=tl*tl
3349 th2=th*th
3350 qqd=0.5*(qpi(kb)-qmi(kb))
3351 qqh=qqd*th2+qmi(kb)*th
3352 qql=qqd*tl2+qmi(kb)*tl
3353 qn(k) = (qqh-qql)/(th-tl)
3354 else if( kt.gt.kb ) then
3355 tl=(zi(k)-za(kb))/dza(kb)
3356 tl2=tl*tl
3357 qqd=0.5*(qpi(kb)-qmi(kb))
3358 qql=qqd*tl2+qmi(kb)*tl
3359 dql = qa(kb)-qql
3360 zsum = (1.-tl)*dza(kb)
3361 qsum = dql*dza(kb)
3362 if( kt-kb.gt.1 ) then
3363 do m=kb+1,kt-1
3364 zsum = zsum + dza(m)
3365 qsum = qsum + qa(m) * dza(m)
3366 enddo
3367 endif
3368 th=(zi(k+1)-za(kt))/dza(kt)
3369 th2=th*th
3370 qqd=0.5*(qpi(kt)-qmi(kt))
3371 dqh=qqd*th2+qmi(kt)*th
3372 zsum = zsum + th*dza(kt)
3373 qsum = qsum + dqh*dza(kt)
3374 qn(k) = qsum/zsum
3375 endif
3376 cycle intp
3377 endif
3378
3379 enddo intp
3380
3381! rain out
3382 sum_precip: do k=1,km
3383 if( za(k).lt.0.0 .and. za(k+1).le.0.0 ) then
3384 precip = precip + qa(k)*dza(k)
3385 net_flx(k) = qa(k)*dza(k)
3386 cycle sum_precip
3387 else if ( za(k).lt.0.0 .and. za(k+1).gt.0.0 ) then
3388 th = (0.0-za(k))/dza(k)
3389 th2 = th*th
3390 qqd = 0.5*(qpi(k)-qmi(k))
3391 qqh = qqd*th2+qmi(k)*th
3392 precip = precip + qqh*dza(k)
3393 net_flx(k) = qqh*dza(k)
3394 exit sum_precip
3395 endif
3396 exit sum_precip
3397 enddo sum_precip
3398
3399! calculating precipitation fluxes
3400 do k=km,1,-1
3401 if(k == km) then
3402 pfsan(k) = net_flx(k)
3403 else
3404 pfsan(k) = pfsan(k+1) + net_flx(k)
3405 end if
3406 enddo
3407!
3408! replace the new values
3409 rql(:) = max(qn(:),r1)
3410
3411 END SUBROUTINE semi_lagrange_sedim
3412
3413 !------------------
3414
3415end module module_mp_tempo_main
module_mp_tempo_main
Definition module_mp_tempo_main.F90:3
module_mp_tempo_params
Definition module_mp_tempo_params.F90:3
module_mp_tempo_utils
Definition module_mp_tempo_utils.F90:3
module_mp_tempo_params::ty_tempo_cfg
Definition module_mp_tempo_params.F90:37