CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
module_mp_tempo_utils.F90
1! Utilities for TEMPO Microphysics
2!=================================================================================================================
3module module_mp_tempo_utils
4
5 use module_mp_tempo_params
6
7#if defined(mpas)
8 use mpas_kind_types, only: wp => rkind, sp => r4kind, dp => r8kind
9 use mp_radar
10#elif defined(standalone)
11 use machine, only: wp => kind_phys, sp => kind_sngl_prec, dp => kind_dbl_prec
12 use module_mp_radar
13#else
14 use machine, only: wp => kind_phys, sp => kind_sngl_prec, dp => kind_dbl_prec
15 use module_mp_radar
16#define ccpp_default 1
17#endif
18
19#if defined(ccpp_default) && defined(MPI)
20 use mpi_f08
21#endif
22
23 implicit none
24
25contains
26 !=================================================================================================================
27 ! Normalized lower gamma function calculated either with a series expansion or continued fraction method
28 ! Input:
29 ! a = gamma function argument, x = upper limit of integration
30 ! Output:
31 ! gamma_p = gamma(a, x) / Gamma(a)
32
33 real function gamma_p(a, x)
34
35 real(wp), intent(in) :: a, x
36
37 if ((x < 0.0) .or. (a <= 0.0)) stop "Invalid arguments for function gamma_p"
38
39 if (x < (a+1.0)) then
40 gamma_p = gamma_series(a, x)
41 else
42 ! gammma_cf computes the upper series
43 gamma_p = 1.0 - gamma_cf(a, x)
44 endif
45
46 end function gamma_p
47
48 !=================================================================================================================
49 ! https://dlmf.nist.gov/8.7 (Equation 8.7.1)
50 ! Solves the normalized lower gamma function: gamma(a,x) / Gamma(a) = x**a * gamma*(a,x)
51 ! Where gamma*(a,x) = exp(-x) * sum (x**k / Gamma(a+k+1))
52 ! Input:
53 ! a = gamma function argument, x = upper limit of integration
54 ! Output:
55 ! Normalized LOWER gamma function: gamma(a, x) / Gamma(a)
56
57 real function gamma_series(a, x)
58
59 real(wp), intent(in) :: a, x
60 integer :: k
61 integer, parameter :: it_max = 100
62 real(wp), parameter :: smallvalue = 1.e-7
63 real(wp) :: ap1, sum_term, sum
64
65 if (x <= 0.0) stop "Invalid arguments for function gamma_series"
66
67 ! k = 0 summation term is 1 / Gamma(a+1)
68 ap1 = a
69 sum_term = 1.0 / gamma(ap1+1.0)
70 sum = sum_term
71 do k = 1, it_max
72 ap1 = ap1 + 1.0
73 sum_term = sum_term * x / ap1
74 sum = sum + sum_term
75 if (abs(sum_term) < (abs(sum) * smallvalue)) exit
76 enddo
77 if (k == it_max) stop "gamma_series solution did not converge"
78
79 gamma_series = sum * x**a * exp(-x)
80
81 end function gamma_series
82
83 !=================================================================================================================
84 ! http://functions.wolfram.com/06.06.10.0003.01
85 ! Solves the normalized upper gamma function: gamma(a,x) / Gamma(a)
86 ! Using a continued fractions method (modifed Lentz Algorithm)
87 ! Input:
88 ! a = gamma function argument, x = lower limit of integration
89 ! Output:
90 ! Normalized UPPER gamma function: gamma(a, x) / Gamma(a)
91
92 real function gamma_cf(a, x)
93
94 real(wp), intent(in) :: a, x
95 integer :: k
96 integer, parameter :: it_max = 100
97 real(wp), parameter :: smallvalue = 1.e-7
98 real(wp), parameter :: offset = 1.e-30
99 real(wp) :: b, d, h0, c, delta, h, aj
100
101 b = 1.0 - a + x
102 d = 1.0 / b
103 h0 = offset
104 c = b + (1.0/offset)
105 delta = c * d
106 h = h0 * delta
107
108 do k = 1, it_max
109 aj = k * (a-k)
110 b = b + 2.0
111 d = b + aj*d
112 if(abs(d) < offset) d = offset
113 c = b + aj/c
114 if(abs(c) < offset) c = offset
115 d = 1.0 / d
116 delta = c * d
117 h = h * delta
118 if (abs(delta-1.0) < smallvalue) exit
119 enddo
120 if (k == it_max) stop "gamma_cf solution did not converge"
121
122 gamma_cf = exp(-x+a*log(x)) * h / gamma(a)
123
124 end function gamma_cf
125
126 !=================================================================================================================
127 ! Calculates log-spaced bins for hydrometer sizes to simplify calculations later
128 ! Input:
129 ! numbins, lowbin, highbin
130 ! Output:
131 ! bins, deltabins
132
133 subroutine create_bins(numbins, lowbin, highbin, bins, deltabins)
134
135 integer, intent(in) :: numbins
136 real(dp), intent(in) :: lowbin, highbin
137
138 integer :: n
139 real(dp), dimension(numbins+1) :: xDx
140 real(dp), dimension(:), intent(out) :: bins
141 real(dp), dimension(:), intent(out), optional :: deltabins
142
143 xdx(1) = lowbin
144 xdx(numbins+1) = highbin
145
146 do n = 2, numbins
147 xdx(n) = exp(real(n-1, kind=dp)/real(numbins, kind=dp) * log(xdx(numbins+1)/xdx(1)) + log(xdx(1)))
148 enddo
149
150 do n = 1, numbins
151 bins(n) = sqrt(xdx(n)*xdx(n+1))
152 enddo
153
154 if (present(deltabins)) then
155 do n = 1, numbins
156 deltabins(n) = xdx(n+1) - xdx(n)
157 enddo
158 endif
159
160 end subroutine create_bins
161
162 !=================================================================================================================
163 ! Variable collision efficiency for rain collecting cloud water using method of Beard and Grover, 1974
164 ! if a/A less than 0.25; otherwise uses polynomials to get close match of Pruppacher & Klett Fig 14-9.
165 ! Output:
166 ! t_Efrw
167
168 subroutine table_efrw
169
170 ! Local variables
171 real(dp) :: vtr, stokes, reynolds, Ef_rw
172 real(dp) :: p, yc0, F, G, H, z, K0, X
173 integer :: i, j
174
175 do j = 1, nbc
176 do i = 1, nbr
177 ef_rw = 0.0
178 p = dc(j) / dr(i)
179 if (dr(i) < 50.e-6 .or. dc(j) < 3.e-6) then
180 t_efrw(i,j) = 0.0
181 elseif (p > 0.25) then
182 x = dc(j) * 1.e6_dp
183 if (dr(i) < 75.e-6) then
184 ef_rw = 0.026794*x - 0.20604
185 elseif (dr(i) < 125.e-6) then
186 ef_rw = -0.00066842*x*x + 0.061542*x - 0.37089
187 elseif (dr(i) < 175.e-6) then
188 ef_rw = 4.091e-06*x*x*x*x - 0.00030908*x*x*x + 0.0066237*x*x - 0.0013687*x - 0.073022
189 elseif (dr(i) < 250.e-6) then
190 ef_rw = 9.6719e-5*x*x*x - 0.0068901*x*x + 0.17305*x - 0.65988
191 elseif (dr(i) < 350.e-6) then
192 ef_rw = 9.0488e-5*x*x*x - 0.006585*x*x + 0.16606*x - 0.56125
193 else
194 ef_rw = 0.00010721*x*x*x - 0.0072962*x*x + 0.1704*x - 0.46929
195 endif
196 else
197 vtr = -0.1021 + 4.932e3*dr(i) - 0.9551e6*dr(i)*dr(i) + 0.07934e9*dr(i)*dr(i)*dr(i) &
198 - 0.002362e12*dr(i)*dr(i)*dr(i)*dr(i)
199 stokes = dc(j) * dc(j) * vtr * rho_w2 / (9.*1.718e-5*dr(i))
200 reynolds = 9. * stokes / (p*p*rho_w2)
201
202 f = log(reynolds)
203 g = -0.1007_dp - 0.358_dp*f + 0.0261_dp*f*f
204 k0 = exp(g)
205 z = log(stokes / (k0+1.e-15_dp))
206 h = 0.1465_dp + 1.302_dp*z - 0.607_dp*z*z + 0.293_dp*z*z*z
207 yc0 = 2.0_dp / pi * atan(h)
208 ef_rw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p))
209
210 endif
211 t_efrw(i,j) = max(0.0, min(real(ef_rw, kind=wp), 0.95))
212 enddo
213 enddo
214
215 end subroutine table_efrw
216
217 !=================================================================================================================
218 ! Variable collision efficiency for snow collecting cloud water using method of Wang and Ji, 2000 except
219 ! equate melted snow diameter to their "effective collision cross-section."
220 ! Output:
221 ! t_Efsw
222
223 subroutine table_efsw
224
225 ! Local variables
226 real(dp) :: Ds_m, vts, vtc, stokes, reynolds, Ef_sw
227 real(dp) :: p, yc0, F, G, H, z, K0
228 integer :: i, j
229
230 do j = 1, nbc
231 vtc = 1.19e4_dp * (1.0e4_dp*dc(j)*dc(j)*0.25_dp)
232 do i = 1, nbs
233 vts = av_s*ds(i)**bv_s * exp(-fv_s*ds(i)) - vtc
234 ds_m = (am_s*ds(i)**bm_s / am_r)**obmr
235 p = dc(j) / ds_m
236 if (p > 0.25 .or. ds(i) < d0s .or. dc(j) < 6.e-6 .or. vts < 1.e-3) then
237 t_efsw(i,j) = 0.0
238 else
239 stokes = dc(j) * dc(j) * vts * rho_w2 / (9.*1.718e-5*ds_m)
240 reynolds = 9. * stokes / (p*p*rho_w2)
241
242 f = log(reynolds)
243 g = -0.1007_dp - 0.358_dp*f + 0.0261_dp*f*f
244 k0 = exp(g)
245 z = log(stokes / (k0+1.e-15_dp))
246 h = 0.1465_dp + 1.302_dp*z - 0.607_dp*z*z + 0.293_dp*z*z*z
247 yc0 = 2.0_dp / pi * atan(h)
248 ef_sw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p))
249
250 t_efsw(i,j) = max(0.0, min(real(ef_sw, kind=wp), 0.95))
251 endif
252 enddo
253 enddo
254
255 end subroutine table_efsw
256
257 !=================================================================================================================
258 ! Droplet evaporation
259 ! Output:
260 ! tpc_wev, tnc_wev
261
262 subroutine table_dropevap
263
264 ! Local variables
265 integer :: i, j, k, n
266 real(dp), dimension(nbc) :: N_c, massc
267 real(dp) :: summ, summ2, lamc, N0_c
268 integer :: nu_c
269
270 do n = 1, nbc
271 massc(n) = am_r*dc(n)**bm_r
272 enddo
273
274 do k = 1, nbc
275 nu_c = min(nu_c_max, nint(nu_c_scale/t_nc(k)) + nu_c_min)
276 do j = 1, ntb_c
277 lamc = (t_nc(k)*am_r* ccg(2,nu_c)*ocg1(nu_c) / r_c(j))**obmr
278 n0_c = t_nc(k)*ocg1(nu_c) * lamc**cce(1,nu_c)
279 do i = 1, nbc
280 n_c(i) = n0_c* dc(i)**nu_c*exp(-lamc*dc(i))*dtc(i)
281 summ = 0.
282 summ2 = 0.
283 do n = 1, i
284 summ = summ + massc(n)*n_c(n)
285 summ2 = summ2 + n_c(n)
286 enddo
287 tpc_wev(i,j,k) = summ
288 tnc_wev(i,j,k) = summ2
289 enddo
290 enddo
291 enddo
292
293 end subroutine table_dropevap
294
295 !=================================================================================================================
296 ! Rain collecting graupel (and inverse). Explicit CE integration.
297
298 subroutine qr_acr_qg(NRHGtable)
299 implicit none
300
301 INTEGER, INTENT(IN) ::NRHGtable
302
303 !..Local variables
304 INTEGER:: i, j, k, m, n, n2, n3, idx_bg
305 INTEGER:: km, km_s, km_e
306 DOUBLE PRECISION, DIMENSION(nbg):: N_g
307 DOUBLE PRECISION, DIMENSION(nbg,NRHGtable):: vg
308 DOUBLE PRECISION, DIMENSION(nbr):: vr, N_r
309 DOUBLE PRECISION:: N0_r, N0_g, lam_exp, lamg, lamr
310 DOUBLE PRECISION:: massg, massr, dvg, dvr, t1, t2, z1, z2, y1, y2
311
312 !+---+
313
314 do n2 = 1, nbr
315 ! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
316 vr(n2) = -0.1021 + 4.932e3*dr(n2) - 0.9551e6*dr(n2)*dr(n2) &
317 + 0.07934e9*dr(n2)*dr(n2)*dr(n2) &
318 - 0.002362e12*dr(n2)*dr(n2)*dr(n2)*dr(n2)
319 enddo
320
321 do n3 = 1, nrhgtable
322 do n = 1, nbg
323 if (nrhgtable == nrhg) idx_bg = n3
324 if (nrhgtable == nrhg1) idx_bg = idx_bg1
325 vg(n,n3) = av_g(idx_bg)*dg(n)**bv_g(idx_bg)
326 enddo
327 enddo
328
329 km_s = 0
330 km_e = ntb_r*ntb_r1 - 1
331
332 do km = km_s, km_e
333 m = km / ntb_r1 + 1
334 k = mod( km , ntb_r1 ) + 1
335
336 lam_exp = (n0r_exp(k)*am_r*crg(1)/r_r(m))**ore1
337 lamr = lam_exp * (crg(3)*org2*org1)**obmr
338 n0_r = n0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2)
339 do n2 = 1, nbr
340 n_r(n2) = n0_r*dr(n2)**mu_r *dexp(-lamr*dr(n2))*dtr(n2)
341 enddo
342
343 do n3 = 1, nrhgtable
344 if (nrhgtable == nrhg) idx_bg = n3
345 if (nrhgtable == nrhg1) idx_bg = idx_bg1
346
347 do j = 1, ntb_g
348 do i = 1, ntb_g1
349 lam_exp = (n0g_exp(i)*am_g(idx_bg)*cgg(1,1)/r_g(j))**oge1
350 lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg
351 n0_g = n0g_exp(i)/(cgg(2,1)*lam_exp) * lamg**cge(2,1)
352 do n = 1, nbg
353 n_g(n) = n0_g*dg(n)**mu_g * dexp(-lamg*dg(n))*dtg(n)
354 enddo
355
356 t1 = 0.0d0
357 t2 = 0.0d0
358 z1 = 0.0d0
359 z2 = 0.0d0
360 y1 = 0.0d0
361 y2 = 0.0d0
362 do n2 = 1, nbr
363 massr = am_r * dr(n2)**bm_r
364 do n = 1, nbg
365 massg = am_g(idx_bg) * dg(n)**bm_g
366
367 dvg = 0.5d0*((vr(n2) - vg(n,n3)) + dabs(vr(n2)-vg(n,n3)))
368 dvr = 0.5d0*((vg(n,n3) - vr(n2)) + dabs(vg(n,n3)-vr(n2)))
369
370 t1 = t1+ pi*.25*ef_rg*(dg(n)+dr(n2))*(dg(n)+dr(n2)) &
371 *dvg*massg * n_g(n)* n_r(n2)
372 z1 = z1+ pi*.25*ef_rg*(dg(n)+dr(n2))*(dg(n)+dr(n2)) &
373 *dvg*massr * n_g(n)* n_r(n2)
374 y1 = y1+ pi*.25*ef_rg*(dg(n)+dr(n2))*(dg(n)+dr(n2)) &
375 *dvg * n_g(n)* n_r(n2)
376
377 t2 = t2+ pi*.25*ef_rg*(dg(n)+dr(n2))*(dg(n)+dr(n2)) &
378 *dvr*massr * n_g(n)* n_r(n2)
379 y2 = y2+ pi*.25*ef_rg*(dg(n)+dr(n2))*(dg(n)+dr(n2)) &
380 *dvr * n_g(n)* n_r(n2)
381 z2 = z2+ pi*.25*ef_rg*(dg(n)+dr(n2))*(dg(n)+dr(n2)) &
382 *dvr*massg * n_g(n)* n_r(n2)
383 enddo
38497 continue
385 enddo
386 tcg_racg(i,j,n3,k,m) = t1
387 tmr_racg(i,j,n3,k,m) = dmin1(z1, r_r(m)*1.0d0)
388 tcr_gacr(i,j,n3,k,m) = t2
389 tnr_racg(i,j,n3,k,m) = y1
390 tnr_gacr(i,j,n3,k,m) = y2
391 enddo
392 enddo
393 enddo
394 enddo
395
396 end subroutine qr_acr_qg
397
398 !=================================================================================================================
399 ! Rain collecting snow (and inverse). Explicit CE integration.
400
401 subroutine qr_acr_qs
402
403 implicit none
404
405 !..Local variables
406 INTEGER:: i, j, k, m, n, n2
407 INTEGER:: km, km_s, km_e
408 DOUBLE PRECISION, DIMENSION(nbr):: vr, D1, N_r
409 DOUBLE PRECISION, DIMENSION(nbs):: vs, N_s
410 DOUBLE PRECISION:: loga_, a_, b_, second, M0, M2, M3, Mrat, oM3
411 DOUBLE PRECISION:: N0_r, lam_exp, lamr, slam1, slam2
412 DOUBLE PRECISION:: dvs, dvr, masss, massr
413 DOUBLE PRECISION:: t1, t2, t3, t4, z1, z2, z3, z4
414 DOUBLE PRECISION:: y1, y2, y3, y4
415
416 !+---+
417
418 do n2 = 1, nbr
419 ! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
420 vr(n2) = -0.1021 + 4.932e3*dr(n2) - 0.9551e6*dr(n2)*dr(n2) &
421 + 0.07934e9*dr(n2)*dr(n2)*dr(n2) &
422 - 0.002362e12*dr(n2)*dr(n2)*dr(n2)*dr(n2)
423 d1(n2) = (vr(n2)/av_s)**(1./bv_s)
424 enddo
425 do n = 1, nbs
426 vs(n) = 1.5*av_s*ds(n)**bv_s * dexp(-fv_s*ds(n))
427 enddo
428
429 km_s = 0
430 km_e = ntb_r*ntb_r1 - 1
431
432 do km = km_s, km_e
433 m = km / ntb_r1 + 1
434 k = mod( km , ntb_r1 ) + 1
435
436 lam_exp = (n0r_exp(k)*am_r*crg(1)/r_r(m))**ore1
437 lamr = lam_exp * (crg(3)*org2*org1)**obmr
438 n0_r = n0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2)
439 do n2 = 1, nbr
440 n_r(n2) = n0_r*dr(n2)**mu_r * dexp(-lamr*dr(n2))*dtr(n2)
441 enddo
442
443 do j = 1, ntb_t
444 do i = 1, ntb_s
445
446 !..From the bm_s moment, compute plus one moment. If we are not
447 !.. using bm_s=2, then we must transform to the pure 2nd moment
448 !.. (variable called "second") and then to the bm_s+1 moment.
449
450 m2 = r_s(i)*oams *1.0d0
451 if (bm_s.gt.2.0-1.e-3 .and. bm_s.lt.2.0+1.e-3) then
452 loga_ = sa(1) + sa(2)*tc(j) + sa(3)*bm_s &
453 + sa(4)*tc(j)*bm_s + sa(5)*tc(j)*tc(j) &
454 + sa(6)*bm_s*bm_s + sa(7)*tc(j)*tc(j)*bm_s &
455 + sa(8)*tc(j)*bm_s*bm_s + sa(9)*tc(j)*tc(j)*tc(j) &
456 + sa(10)*bm_s*bm_s*bm_s
457 a_ = 10.0**loga_
458 b_ = sb(1) + sb(2)*tc(j) + sb(3)*bm_s &
459 + sb(4)*tc(j)*bm_s + sb(5)*tc(j)*tc(j) &
460 + sb(6)*bm_s*bm_s + sb(7)*tc(j)*tc(j)*bm_s &
461 + sb(8)*tc(j)*bm_s*bm_s + sb(9)*tc(j)*tc(j)*tc(j) &
462 + sb(10)*bm_s*bm_s*bm_s
463 second = (m2/a_)**(1./b_)
464 else
465 second = m2
466 endif
467
468 loga_ = sa(1) + sa(2)*tc(j) + sa(3)*cse(1) &
469 + sa(4)*tc(j)*cse(1) + sa(5)*tc(j)*tc(j) &
470 + sa(6)*cse(1)*cse(1) + sa(7)*tc(j)*tc(j)*cse(1) &
471 + sa(8)*tc(j)*cse(1)*cse(1) + sa(9)*tc(j)*tc(j)*tc(j) &
472 + sa(10)*cse(1)*cse(1)*cse(1)
473 a_ = 10.0**loga_
474 b_ = sb(1)+sb(2)*tc(j)+sb(3)*cse(1) + sb(4)*tc(j)*cse(1) &
475 + sb(5)*tc(j)*tc(j) + sb(6)*cse(1)*cse(1) &
476 + sb(7)*tc(j)*tc(j)*cse(1) + sb(8)*tc(j)*cse(1)*cse(1) &
477 + sb(9)*tc(j)*tc(j)*tc(j)+sb(10)*cse(1)*cse(1)*cse(1)
478 m3 = a_ * second**b_
479
480 om3 = 1./m3
481 mrat = m2*(m2*om3)*(m2*om3)*(m2*om3)
482 m0 = (m2*om3)**mu_s
483 slam1 = m2 * om3 * lam0
484 slam2 = m2 * om3 * lam1
485
486 do n = 1, nbs
487 n_s(n) = mrat*(kap0*dexp(-slam1*ds(n)) &
488 + kap1*m0*ds(n)**mu_s * dexp(-slam2*ds(n)))*dts(n)
489 enddo
490
491 t1 = 0.0d0
492 t2 = 0.0d0
493 t3 = 0.0d0
494 t4 = 0.0d0
495 z1 = 0.0d0
496 z2 = 0.0d0
497 z3 = 0.0d0
498 z4 = 0.0d0
499 y1 = 0.0d0
500 y2 = 0.0d0
501 y3 = 0.0d0
502 y4 = 0.0d0
503 do n2 = 1, nbr
504 massr = am_r * dr(n2)**bm_r
505 do n = 1, nbs
506 masss = am_s * ds(n)**bm_s
507
508 dvs = 0.5d0*((vr(n2) - vs(n)) + dabs(vr(n2)-vs(n)))
509 dvr = 0.5d0*((vs(n) - vr(n2)) + dabs(vs(n)-vr(n2)))
510 if (massr .gt. 1.5*masss) then
511 t1 = t1+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
512 *dvs*masss * n_s(n)* n_r(n2)
513 z1 = z1+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
514 *dvs*massr * n_s(n)* n_r(n2)
515 y1 = y1+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
516 *dvs * n_s(n)* n_r(n2)
517 else
518 t3 = t3+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
519 *dvs*masss * n_s(n)* n_r(n2)
520 z3 = z3+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
521 *dvs*massr * n_s(n)* n_r(n2)
522 y3 = y3+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
523 *dvs * n_s(n)* n_r(n2)
524 endif
525
526 if (massr .gt. 1.5*masss) then
527 t2 = t2+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
528 *dvr*massr * n_s(n)* n_r(n2)
529 y2 = y2+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
530 *dvr * n_s(n)* n_r(n2)
531 z2 = z2+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
532 *dvr*masss * n_s(n)* n_r(n2)
533 else
534 t4 = t4+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
535 *dvr*massr * n_s(n)* n_r(n2)
536 y4 = y4+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
537 *dvr * n_s(n)* n_r(n2)
538 z4 = z4+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
539 *dvr*masss * n_s(n)* n_r(n2)
540 endif
541
542 enddo
543 enddo
544 tcs_racs1(i,j,k,m) = t1
545 tmr_racs1(i,j,k,m) = dmin1(z1, r_r(m)*1.0d0)
546 tcs_racs2(i,j,k,m) = t3
547 tmr_racs2(i,j,k,m) = z3
548 tcr_sacr1(i,j,k,m) = t2
549 tms_sacr1(i,j,k,m) = z2
550 tcr_sacr2(i,j,k,m) = t4
551 tms_sacr2(i,j,k,m) = z4
552 tnr_racs1(i,j,k,m) = y1
553 tnr_racs2(i,j,k,m) = y3
554 tnr_sacr1(i,j,k,m) = y2
555 tnr_sacr2(i,j,k,m) = y4
556 enddo
557 enddo
558 enddo
559
560 end subroutine qr_acr_qs
561
562 !=================================================================================================================
563 ! This is a literal adaptation of Bigg (1954) probability of drops of a particular volume freezing.
564 ! Given this probability, simply freeze the proportion of drops summing their masses.
565
566 subroutine freezeh2o
567
568 implicit none
569
570 !..Local variables
571 INTEGER:: i, j, k, m, n, n2
572 DOUBLE PRECISION :: N_r, N_c
573 DOUBLE PRECISION, DIMENSION(nbr):: massr
574 DOUBLE PRECISION, DIMENSION(nbc):: massc
575 DOUBLE PRECISION:: sum1, sum2, sumn1, sumn2, &
576 prob, vol, Texp, orho_w, &
577 lam_exp, lamr, N0_r, lamc, N0_c, y
578 INTEGER:: nu_c
579 REAL:: T_adjust
580
581 orho_w = 1./rho_w2
582
583 do n2 = 1, nbr
584 massr(n2) = am_r*dr(n2)**bm_r
585 enddo
586 do n = 1, nbc
587 massc(n) = am_r*dc(n)**bm_r
588 enddo
589
590 !..Freeze water (smallest drops become cloud ice, otherwise graupel).
591 do m = 1, ntb_in
592 t_adjust = max(-3.0, min(3.0 - alog10(nt_in(m)), 3.0))
593 do k = 1, 45
594 ! print*, ' Freezing water for temp = ', -k
595 texp = dexp( real(k,kind=dp) - t_adjust*1.0d0 ) - 1.0d0
596 do j = 1, ntb_r1
597 do i = 1, ntb_r
598 lam_exp = (n0r_exp(j)*am_r*crg(1)/r_r(i))**ore1
599 lamr = lam_exp * (crg(3)*org2*org1)**obmr
600 n0_r = n0r_exp(j)/(crg(2)*lam_exp) * lamr**cre(2)
601 sum1 = 0.0d0
602 sum2 = 0.0d0
603 sumn1 = 0.0d0
604 sumn2 = 0.0d0
605 do n2 = nbr, 1, -1
606 n_r = n0_r*dr(n2)**mu_r*dexp(-lamr*dr(n2))*dtr(n2)
607 vol = massr(n2)*orho_w
608 prob = max(0.0d0, 1.0d0 - dexp(-120.0d0*vol*5.2d-4 * texp))
609 if (massr(n2) .lt. xm0g) then
610 sumn1 = sumn1 + prob*n_r
611 sum1 = sum1 + prob*n_r*massr(n2)
612 else
613 sumn2 = sumn2 + prob*n_r
614 sum2 = sum2 + prob*n_r*massr(n2)
615 endif
616 if ((sum1+sum2).ge.r_r(i)) EXIT
617 enddo
618 tpi_qrfz(i,j,k,m) = sum1
619 tni_qrfz(i,j,k,m) = sumn1
620 tpg_qrfz(i,j,k,m) = sum2
621 tnr_qrfz(i,j,k,m) = sumn2
622 enddo
623 enddo
624
625 do j = 1, nbc
626 nu_c = min(15, nint(1000.e6/t_nc(j)) + 2)
627 do i = 1, ntb_c
628 lamc = (t_nc(j)*am_r* ccg(2,nu_c) * ocg1(nu_c) / r_c(i))**obmr
629 n0_c = t_nc(j)*ocg1(nu_c) * lamc**cce(1,nu_c)
630 sum1 = 0.0d0
631 sumn2 = 0.0d0
632 do n = nbc, 1, -1
633 vol = massc(n)*orho_w
634 prob = max(0.0d0, 1.0d0 - dexp(-120.0d0*vol*5.2d-4 * texp))
635 n_c = n0_c*dc(n)**nu_c*exp(-lamc*dc(n))*dtc(n)
636 sumn2 = min(t_nc(j), sumn2 + prob*n_c)
637 sum1 = sum1 + prob*n_c*massc(n)
638 if (sum1 .ge. r_c(i)) EXIT
639 enddo
640 tpi_qcfz(i,j,k,m) = sum1
641 tni_qcfz(i,j,k,m) = sumn2
642 enddo
643 enddo
644 enddo
645 enddo
646
647 end subroutine freezeh2o
648
649 !=================================================================================================================
650 ! Cloud ice converting to snow since portion greater than min snow size. Given cloud ice content (kg/m**3),
651 ! number concentration (#/m**3) and gamma shape parameter, mu_i, break the distrib into bins and figure out
652 ! the mass/number of ice with sizes larger than D0s. Also, compute incomplete gamma function for the
653 ! integration of ice depositional growth from diameter=0 to D0s. Amount of ice depositional growth is this
654 ! portion of distrib while larger diameters contribute to snow growth (as in Harrington et al. 1995).
655
656 subroutine qi_aut_qs
657
658 implicit none
659
660 !..Local variables
661 INTEGER:: i, j, n2
662 DOUBLE PRECISION, DIMENSION(nbi):: N_i
663 DOUBLE PRECISION:: N0_i, lami, Di_mean, t1, t2
664 REAL:: xlimit_intg
665
666 do j = 1, ntb_i1
667 do i = 1, ntb_i
668 lami = (am_i*cig(2)*oig1*nt_i(j)/r_i(i))**obmi
669 di_mean = (bm_i + mu_i + 1.) / lami
670 n0_i = nt_i(j)*oig1 * lami**cie(1)
671 t1 = 0.0d0
672 t2 = 0.0d0
673 if (sngl(di_mean) .gt. 5.*d0s) then
674 t1 = r_i(i)
675 t2 = nt_i(j)
676 tpi_ide(i,j) = 0.0d0
677 elseif (sngl(di_mean) .lt. d0i) then
678 t1 = 0.0d0
679 t2 = 0.0d0
680 tpi_ide(i,j) = 1.0d0
681 else
682 xlimit_intg = lami*d0s
683 tpi_ide(i,j) = gamma_p(mu_i+2.0, xlimit_intg) * 1.0d0
684 do n2 = 1, nbi
685 n_i(n2) = n0_i*di(n2)**mu_i * dexp(-lami*di(n2))*dti(n2)
686 if (di(n2).ge.d0s) then
687 t1 = t1 + n_i(n2) * am_i*di(n2)**bm_i
688 t2 = t2 + n_i(n2)
689 endif
690 enddo
691 endif
692 tps_iaus(i,j) = t1
693 tni_iaus(i,j) = t2
694 enddo
695 enddo
696
697 end subroutine qi_aut_qs
698
699#if defined (ccpp_default)
700 !=================================================================================================================
701 !..Fill the table of CCN activation data created from parcel model run
702 !.. by Trude Eidhammer with inputs of aerosol number concentration,
703 !.. vertical velocity, temperature, lognormal mean aerosol radius, and
704 !.. hygroscopicity, kappa. The data are read from external file and
705 !.. contain activated fraction of CCN for given conditions.
706 !+---+-----------------------------------------------------------------+
707 !+---+-----------------------------------------------------------------+
714
715 subroutine table_ccnact(errmess,errflag)
716
717 implicit none
718
719 !..Error handling variables
720 CHARACTER(len=*), INTENT(INOUT) :: errmess
721 INTEGER, INTENT(INOUT) :: errflag
722
723 !..Local variables
724 INTEGER:: iunit_mp_th1, i
725 LOGICAL:: opened
726
727 iunit_mp_th1 = -1
728 DO i = 20,99
729 INQUIRE ( i , opened = opened )
730 IF ( .NOT. opened ) THEN
731 iunit_mp_th1 = i
732 GOTO 2010
733 ENDIF
734 ENDDO
7352010 CONTINUE
736 IF ( iunit_mp_th1 < 0 ) THEN
737 write(0,*) 'module_mp_tempo: table_ccnAct: '// &
738 'Can not find unused fortran unit to read in lookup table.'
739 return
740 ENDIF
741
742 !WRITE(*, '(A,I2)') 'module_mp_tempo: opening CCN_ACTIVATE.BIN on unit ',iunit_mp_th1
743 OPEN(iunit_mp_th1,file='CCN_ACTIVATE.BIN', &
744 form='UNFORMATTED',status='OLD',convert='BIG_ENDIAN',err=9009)
745
746 !sms$serial begin
747 READ(iunit_mp_th1,err=9010) tnccn_act
748 !sms$serial end
749
750 RETURN
7519009 CONTINUE
752 WRITE( errmess , '(A,I2)' ) 'module_mp_tempo: error opening CCN_ACTIVATE.BIN on unit ',iunit_mp_th1
753 errflag = 1
754 RETURN
7559010 CONTINUE
756 WRITE( errmess , '(A,I2)' ) 'module_mp_tempo: error reading CCN_ACTIVATE.BIN on unit ',iunit_mp_th1
757 errflag = 1
758 RETURN
759
760 end subroutine table_ccnact
761
762 !=================================================================================================================
763 ! Rain collecting graupel (and inverse). Explicit CE integration.
764 subroutine qr_acr_qg_par(NRHGtable, qrqg_file)
765
766 implicit none
767
768 INTEGER, INTENT(IN) ::NRHGtable
769
770 !..Local variables
771 INTEGER:: i, j, k, m, n, n2, n3, idx_bg
772 INTEGER:: km, km_s, km_e
773 DOUBLE PRECISION, DIMENSION(nbg):: N_g
774 DOUBLE PRECISION, DIMENSION(nbg,NRHGtable):: vg
775 DOUBLE PRECISION, DIMENSION(nbr):: vr, N_r
776 DOUBLE PRECISION:: N0_r, N0_g, lam_exp, lamg, lamr
777 DOUBLE PRECISION:: massg, massr, dvg, dvr, t1, t2, z1, z2, y1, y2
778 LOGICAL force_read_thompson, write_thompson_tables
779 LOGICAL lexist,lopen
780 INTEGER good,ierr
781 character(len=*), intent(in) :: qrqg_file
782
783 force_read_thompson = .false.
784 write_thompson_tables = .false.
785 !+---+
786
787 good = 0
788 INQUIRE(file=qrqg_file, exist=lexist)
789#ifdef MPI
790 call mpi_barrier(mpi_communicator,ierr)
791#endif
792 IF ( lexist ) THEN
793 OPEN(63,file=qrqg_file,form="unformatted",err=1234)
794 !sms$serial begin
795 READ(63,err=1234) tcg_racg
796 READ(63,err=1234) tmr_racg
797 READ(63,err=1234) tcr_gacr
798!! READ(63,err=1234) tmg_gacr
799 READ(63,err=1234) tnr_racg
800 READ(63,err=1234) tnr_gacr
801 !sms$serial end
802 good = 1
8031234 CONTINUE
804 IF ( good .NE. 1 ) THEN
805 INQUIRE(63,opened=lopen)
806 IF (lopen) THEN
807 IF( force_read_thompson ) THEN
808 write(0,*) "Error reading "//qrqg_file//" Aborting because force_read_thompson is .true."
809 return
810 ENDIF
811 CLOSE(63)
812 ELSE
813 IF( force_read_thompson ) THEN
814 write(0,*) "Error opening "//qrqg_file//" Aborting because force_read_thompson is .true."
815 return
816 ENDIF
817 ENDIF
818 ELSE
819 INQUIRE(63,opened=lopen)
820 IF (lopen) THEN
821 CLOSE(63)
822 ENDIF
823 ENDIF
824 ELSE
825 IF( force_read_thompson ) THEN
826 write(0,*) "Non-existent "//qrqg_file//" Aborting because force_read_thompson is .true."
827 return
828 ENDIF
829 ENDIF
830
831 IF (.NOT. good .EQ. 1 ) THEN
832 if (thompson_table_writer) then
833 write_thompson_tables = .true.
834 write(0,*) "ThompMP: computing qr_acr_qg"
835 endif
836 do n2 = 1, nbr
837 ! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
838 vr(n2) = -0.1021 + 4.932e3*dr(n2) - 0.9551e6*dr(n2)*dr(n2) &
839 + 0.07934e9*dr(n2)*dr(n2)*dr(n2) &
840 - 0.002362e12*dr(n2)*dr(n2)*dr(n2)*dr(n2)
841 enddo
842 ! do n = 1, nbg
843 ! vg(n) = av_g*Dg(n)**bv_g
844 ! enddo
845
846 do n3 = 1, nrhgtable
847 do n = 1, nbg
848 if (nrhgtable == nrhg) idx_bg = n3
849 if (nrhgtable == nrhg1) idx_bg = idx_bg1
850 vg(n,n3) = av_g(idx_bg)*dg(n)**bv_g(idx_bg)
851 enddo
852 enddo
853
854 !..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for
855 !.. fortran indices. J. Michalakes, 2009Oct30.
856
857#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
858 CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e )
859#else
860 km_s = 0
861 km_e = ntb_r*ntb_r1 - 1
862#endif
863
864 do km = km_s, km_e
865 m = km / ntb_r1 + 1
866 k = mod( km , ntb_r1 ) + 1
867
868 lam_exp = (n0r_exp(k)*am_r*crg(1)/r_r(m))**ore1
869 lamr = lam_exp * (crg(3)*org2*org1)**obmr
870 n0_r = n0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2)
871 do n2 = 1, nbr
872 n_r(n2) = n0_r*dr(n2)**mu_r *dexp(-lamr*dr(n2))*dtr(n2)
873 enddo
874
875 do n3 = 1, nrhgtable
876 if (nrhgtable == nrhg) idx_bg = n3
877 if (nrhgtable == nrhg1) idx_bg = idx_bg1
878
879 do j = 1, ntb_g
880 do i = 1, ntb_g1
881 lam_exp = (n0g_exp(i)*am_g(idx_bg)*cgg(1,1)/r_g(j))**oge1
882 lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg
883 n0_g = n0g_exp(i)/(cgg(2,1)*lam_exp) * lamg**cge(2,1)
884 do n = 1, nbg
885 n_g(n) = n0_g*dg(n)**mu_g * dexp(-lamg*dg(n))*dtg(n)
886 enddo
887
888 t1 = 0.0d0
889 t2 = 0.0d0
890 z1 = 0.0d0
891 z2 = 0.0d0
892 y1 = 0.0d0
893 y2 = 0.0d0
894 do n2 = 1, nbr
895 massr = am_r * dr(n2)**bm_r
896 do n = 1, nbg
897 massg = am_g(idx_bg) * dg(n)**bm_g
898
899 dvg = 0.5d0*((vr(n2) - vg(n,n3)) + dabs(vr(n2)-vg(n,n3)))
900 dvr = 0.5d0*((vg(n,n3) - vr(n2)) + dabs(vg(n,n3)-vr(n2)))
901
902 t1 = t1+ pi*.25*ef_rg*(dg(n)+dr(n2))*(dg(n)+dr(n2)) &
903 *dvg*massg * n_g(n)* n_r(n2)
904 z1 = z1+ pi*.25*ef_rg*(dg(n)+dr(n2))*(dg(n)+dr(n2)) &
905 *dvg*massr * n_g(n)* n_r(n2)
906 y1 = y1+ pi*.25*ef_rg*(dg(n)+dr(n2))*(dg(n)+dr(n2)) &
907 *dvg * n_g(n)* n_r(n2)
908
909 t2 = t2+ pi*.25*ef_rg*(dg(n)+dr(n2))*(dg(n)+dr(n2)) &
910 *dvr*massr * n_g(n)* n_r(n2)
911 y2 = y2+ pi*.25*ef_rg*(dg(n)+dr(n2))*(dg(n)+dr(n2)) &
912 *dvr * n_g(n)* n_r(n2)
913 z2 = z2+ pi*.25*ef_rg*(dg(n)+dr(n2))*(dg(n)+dr(n2)) &
914 *dvr*massg * n_g(n)* n_r(n2)
915 enddo
91697 continue
917 enddo
918 tcg_racg(i,j,n3,k,m) = t1
919 tmr_racg(i,j,n3,k,m) = dmin1(z1, r_r(m)*1.0d0)
920 tcr_gacr(i,j,n3,k,m) = t2
921!! tmg_gacr(i,j,k,m) = DMIN1(z2, r_g(j)*1.0d0)
922 tnr_racg(i,j,n3,k,m) = y1
923 tnr_gacr(i,j,n3,k,m) = y2
924 enddo
925 enddo
926 enddo
927 enddo
928 IF ( write_thompson_tables ) THEN
929 write(0,*) "Writing "//qrqg_file//" in Tempo MP init"
930 OPEN(63,file=qrqg_file,form="unformatted",err=9234)
931 WRITE(63,err=9234) tcg_racg
932 WRITE(63,err=9234) tmr_racg
933 WRITE(63,err=9234) tcr_gacr
934 WRITE(63,err=9234) tnr_racg
935 WRITE(63,err=9234) tnr_gacr
936 CLOSE(63)
937 RETURN ! ----- RETURN
9389234 CONTINUE
939 write(0,*) "Error writing "//qrqg_file
940 return
941 ENDIF
942 ENDIF
943
944 end subroutine qr_acr_qg_par
945
946 !=================================================================================================================
947 ! Rain collecting snow (and inverse). Explicit CE integration.
948
949 subroutine qr_acr_qs_par
950
951 implicit none
952
953 !..Local variables
954 INTEGER:: i, j, k, m, n, n2
955 INTEGER:: km, km_s, km_e
956 DOUBLE PRECISION, DIMENSION(nbr):: vr, D1, N_r
957 DOUBLE PRECISION, DIMENSION(nbs):: vs, N_s
958 DOUBLE PRECISION:: loga_, a_, b_, second, M0, M2, M3, Mrat, oM3
959 DOUBLE PRECISION:: N0_r, lam_exp, lamr, slam1, slam2
960 DOUBLE PRECISION:: dvs, dvr, masss, massr
961 DOUBLE PRECISION:: t1, t2, t3, t4, z1, z2, z3, z4
962 DOUBLE PRECISION:: y1, y2, y3, y4
963 LOGICAL force_read_thompson, write_thompson_tables
964 LOGICAL lexist,lopen
965 INTEGER good,ierr
966
967 !+---+
968
969 force_read_thompson = .false.
970 write_thompson_tables = .false.
971
972 good = 0
973 INQUIRE(file=qr_acr_qs_file, exist=lexist)
974#ifdef MPI
975 call mpi_barrier(mpi_communicator,ierr)
976#endif
977 IF ( lexist ) THEN
978 !write(0,*) "ThompMP: read "//qr_acr_qs_file//" instead of computing"
979 OPEN(63,file=qr_acr_qs_file,form="unformatted",err=1234)
980 !sms$serial begin
981 READ(63,err=1234)tcs_racs1
982 READ(63,err=1234)tmr_racs1
983 READ(63,err=1234)tcs_racs2
984 READ(63,err=1234)tmr_racs2
985 READ(63,err=1234)tcr_sacr1
986 READ(63,err=1234)tms_sacr1
987 READ(63,err=1234)tcr_sacr2
988 READ(63,err=1234)tms_sacr2
989 READ(63,err=1234)tnr_racs1
990 READ(63,err=1234)tnr_racs2
991 READ(63,err=1234)tnr_sacr1
992 READ(63,err=1234)tnr_sacr2
993 !sms$serial end
994 good = 1
9951234 CONTINUE
996 IF ( good .NE. 1 ) THEN
997 INQUIRE(63,opened=lopen)
998 IF (lopen) THEN
999 IF( force_read_thompson ) THEN
1000 write(0,*) "Error reading "//qr_acr_qs_file//" Aborting because force_read_thompson is .true."
1001 return
1002 ENDIF
1003 CLOSE(63)
1004 ELSE
1005 IF( force_read_thompson ) THEN
1006 write(0,*) "Error opening "//qr_acr_qs_file//" Aborting because force_read_thompson is .true."
1007 return
1008 ENDIF
1009 ENDIF
1010 ELSE
1011 INQUIRE(63,opened=lopen)
1012 IF (lopen) THEN
1013 CLOSE(63)
1014 ENDIF
1015 ENDIF
1016 ELSE
1017 IF( force_read_thompson ) THEN
1018 write(0,*) "Non-existent "//qr_acr_qs_file//" Aborting because force_read_thompson is .true."
1019 return
1020 ENDIF
1021 ENDIF
1022
1023 IF (.NOT. good .EQ. 1 ) THEN
1024 if (thompson_table_writer) then
1025 write_thompson_tables = .true.
1026 write(0,*) "ThompMP: computing qr_acr_qs"
1027 endif
1028 do n2 = 1, nbr
1029 ! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2))
1030 vr(n2) = -0.1021 + 4.932e3*dr(n2) - 0.9551e6*dr(n2)*dr(n2) &
1031 + 0.07934e9*dr(n2)*dr(n2)*dr(n2) &
1032 - 0.002362e12*dr(n2)*dr(n2)*dr(n2)*dr(n2)
1033 d1(n2) = (vr(n2)/av_s)**(1./bv_s)
1034 enddo
1035 do n = 1, nbs
1036 vs(n) = 1.5*av_s*ds(n)**bv_s * dexp(-fv_s*ds(n))
1037 enddo
1038
1039 !..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for
1040 !.. fortran indices. J. Michalakes, 2009Oct30.
1041
1042#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1043 CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e )
1044#else
1045 km_s = 0
1046 km_e = ntb_r*ntb_r1 - 1
1047#endif
1048
1049 do km = km_s, km_e
1050 m = km / ntb_r1 + 1
1051 k = mod( km , ntb_r1 ) + 1
1052
1053 lam_exp = (n0r_exp(k)*am_r*crg(1)/r_r(m))**ore1
1054 lamr = lam_exp * (crg(3)*org2*org1)**obmr
1055 n0_r = n0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2)
1056 do n2 = 1, nbr
1057 n_r(n2) = n0_r*dr(n2)**mu_r * dexp(-lamr*dr(n2))*dtr(n2)
1058 enddo
1059
1060 do j = 1, ntb_t
1061 do i = 1, ntb_s
1062
1063 !..From the bm_s moment, compute plus one moment. If we are not
1064 !.. using bm_s=2, then we must transform to the pure 2nd moment
1065 !.. (variable called "second") and then to the bm_s+1 moment.
1066
1067 m2 = r_s(i)*oams *1.0d0
1068 if (bm_s.gt.2.0-1.e-3 .and. bm_s.lt.2.0+1.e-3) then
1069 loga_ = sa(1) + sa(2)*tc(j) + sa(3)*bm_s &
1070 + sa(4)*tc(j)*bm_s + sa(5)*tc(j)*tc(j) &
1071 + sa(6)*bm_s*bm_s + sa(7)*tc(j)*tc(j)*bm_s &
1072 + sa(8)*tc(j)*bm_s*bm_s + sa(9)*tc(j)*tc(j)*tc(j) &
1073 + sa(10)*bm_s*bm_s*bm_s
1074 a_ = 10.0**loga_
1075 b_ = sb(1) + sb(2)*tc(j) + sb(3)*bm_s &
1076 + sb(4)*tc(j)*bm_s + sb(5)*tc(j)*tc(j) &
1077 + sb(6)*bm_s*bm_s + sb(7)*tc(j)*tc(j)*bm_s &
1078 + sb(8)*tc(j)*bm_s*bm_s + sb(9)*tc(j)*tc(j)*tc(j) &
1079 + sb(10)*bm_s*bm_s*bm_s
1080 second = (m2/a_)**(1./b_)
1081 else
1082 second = m2
1083 endif
1084
1085 loga_ = sa(1) + sa(2)*tc(j) + sa(3)*cse(1) &
1086 + sa(4)*tc(j)*cse(1) + sa(5)*tc(j)*tc(j) &
1087 + sa(6)*cse(1)*cse(1) + sa(7)*tc(j)*tc(j)*cse(1) &
1088 + sa(8)*tc(j)*cse(1)*cse(1) + sa(9)*tc(j)*tc(j)*tc(j) &
1089 + sa(10)*cse(1)*cse(1)*cse(1)
1090 a_ = 10.0**loga_
1091 b_ = sb(1)+sb(2)*tc(j)+sb(3)*cse(1) + sb(4)*tc(j)*cse(1) &
1092 + sb(5)*tc(j)*tc(j) + sb(6)*cse(1)*cse(1) &
1093 + sb(7)*tc(j)*tc(j)*cse(1) + sb(8)*tc(j)*cse(1)*cse(1) &
1094 + sb(9)*tc(j)*tc(j)*tc(j)+sb(10)*cse(1)*cse(1)*cse(1)
1095 m3 = a_ * second**b_
1096
1097 om3 = 1./m3
1098 mrat = m2*(m2*om3)*(m2*om3)*(m2*om3)
1099 m0 = (m2*om3)**mu_s
1100 slam1 = m2 * om3 * lam0
1101 slam2 = m2 * om3 * lam1
1102
1103 do n = 1, nbs
1104 n_s(n) = mrat*(kap0*dexp(-slam1*ds(n)) &
1105 + kap1*m0*ds(n)**mu_s * dexp(-slam2*ds(n)))*dts(n)
1106 enddo
1107
1108 t1 = 0.0d0
1109 t2 = 0.0d0
1110 t3 = 0.0d0
1111 t4 = 0.0d0
1112 z1 = 0.0d0
1113 z2 = 0.0d0
1114 z3 = 0.0d0
1115 z4 = 0.0d0
1116 y1 = 0.0d0
1117 y2 = 0.0d0
1118 y3 = 0.0d0
1119 y4 = 0.0d0
1120 do n2 = 1, nbr
1121 massr = am_r * dr(n2)**bm_r
1122 do n = 1, nbs
1123 masss = am_s * ds(n)**bm_s
1124
1125 dvs = 0.5d0*((vr(n2) - vs(n)) + dabs(vr(n2)-vs(n)))
1126 dvr = 0.5d0*((vs(n) - vr(n2)) + dabs(vs(n)-vr(n2)))
1127
1128 if (massr .gt. 1.5*masss) then
1129 t1 = t1+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
1130 *dvs*masss * n_s(n)* n_r(n2)
1131 z1 = z1+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
1132 *dvs*massr * n_s(n)* n_r(n2)
1133 y1 = y1+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
1134 *dvs * n_s(n)* n_r(n2)
1135 else
1136 t3 = t3+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
1137 *dvs*masss * n_s(n)* n_r(n2)
1138 z3 = z3+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
1139 *dvs*massr * n_s(n)* n_r(n2)
1140 y3 = y3+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
1141 *dvs * n_s(n)* n_r(n2)
1142 endif
1143
1144 if (massr .gt. 1.5*masss) then
1145 t2 = t2+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
1146 *dvr*massr * n_s(n)* n_r(n2)
1147 y2 = y2+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
1148 *dvr * n_s(n)* n_r(n2)
1149 z2 = z2+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
1150 *dvr*masss * n_s(n)* n_r(n2)
1151 else
1152 t4 = t4+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
1153 *dvr*massr * n_s(n)* n_r(n2)
1154 y4 = y4+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
1155 *dvr * n_s(n)* n_r(n2)
1156 z4 = z4+ pi*.25*ef_rs*(ds(n)+dr(n2))*(ds(n)+dr(n2)) &
1157 *dvr*masss * n_s(n)* n_r(n2)
1158 endif
1159
1160 enddo
1161 enddo
1162 tcs_racs1(i,j,k,m) = t1
1163 tmr_racs1(i,j,k,m) = dmin1(z1, r_r(m)*1.0d0)
1164 tcs_racs2(i,j,k,m) = t3
1165 tmr_racs2(i,j,k,m) = z3
1166 tcr_sacr1(i,j,k,m) = t2
1167 tms_sacr1(i,j,k,m) = z2
1168 tcr_sacr2(i,j,k,m) = t4
1169 tms_sacr2(i,j,k,m) = z4
1170 tnr_racs1(i,j,k,m) = y1
1171 tnr_racs2(i,j,k,m) = y3
1172 tnr_sacr1(i,j,k,m) = y2
1173 tnr_sacr2(i,j,k,m) = y4
1174 enddo
1175 enddo
1176 enddo
1177
1178 IF ( write_thompson_tables ) THEN
1179 write(0,*) "Writing "//qr_acr_qs_file//" in Tempo MP init"
1180 OPEN(63,file=qr_acr_qs_file,form="unformatted",err=9234)
1181 WRITE(63,err=9234)tcs_racs1
1182 WRITE(63,err=9234)tmr_racs1
1183 WRITE(63,err=9234)tcs_racs2
1184 WRITE(63,err=9234)tmr_racs2
1185 WRITE(63,err=9234)tcr_sacr1
1186 WRITE(63,err=9234)tms_sacr1
1187 WRITE(63,err=9234)tcr_sacr2
1188 WRITE(63,err=9234)tms_sacr2
1189 WRITE(63,err=9234)tnr_racs1
1190 WRITE(63,err=9234)tnr_racs2
1191 WRITE(63,err=9234)tnr_sacr1
1192 WRITE(63,err=9234)tnr_sacr2
1193 CLOSE(63)
1194 RETURN ! ----- RETURN
11959234 CONTINUE
1196 write(0,*) "Error writing "//qr_acr_qs_file
1197 ENDIF
1198 ENDIF
1199
1200 end subroutine qr_acr_qs_par
1201
1202 !=================================================================================================================
1203 !! This is a literal adaptation of Bigg (1954) probability of drops of
1204 !! a particular volume freezing. Given this probability, simply freeze
1205 !! the proportion of drops summing their masses.
1206
1207 subroutine freezeh2o_par(threads)
1208
1209 implicit none
1210
1211 !..Interface variables
1212 INTEGER, INTENT(IN):: threads
1213
1214 !..Local variables
1215 INTEGER:: i, j, k, m, n, n2
1216 DOUBLE PRECISION:: N_r, N_c
1217 DOUBLE PRECISION, DIMENSION(nbr):: massr
1218 DOUBLE PRECISION, DIMENSION(nbc):: massc
1219 DOUBLE PRECISION:: sum1, sum2, sumn1, sumn2, &
1220 prob, vol, Texp, orho_w, &
1221 lam_exp, lamr, N0_r, lamc, N0_c, y
1222 INTEGER:: nu_c
1223 REAL:: T_adjust
1224 LOGICAL force_read_thompson, write_thompson_tables
1225 LOGICAL lexist,lopen
1226 INTEGER good,ierr
1227
1228 !+---+
1229 force_read_thompson = .false.
1230 write_thompson_tables = .false.
1231
1232 good = 0
1233 INQUIRE(file=freeze_h2o_file,exist=lexist)
1234#ifdef MPI
1235 call mpi_barrier(mpi_communicator,ierr)
1236#endif
1237 IF ( lexist ) THEN
1238 !write(0,*) "ThompMP: read "//freeze_h2o_file//" instead of computing"
1239 OPEN(63,file=freeze_h2o_file,form="unformatted",err=1234)
1240 !sms$serial begin
1241 READ(63,err=1234)tpi_qrfz
1242 READ(63,err=1234)tni_qrfz
1243 READ(63,err=1234)tpg_qrfz
1244 READ(63,err=1234)tnr_qrfz
1245 READ(63,err=1234)tpi_qcfz
1246 READ(63,err=1234)tni_qcfz
1247 !sms$serial end
1248 good = 1
12491234 CONTINUE
1250 IF ( good .NE. 1 ) THEN
1251 INQUIRE(63,opened=lopen)
1252 IF (lopen) THEN
1253 IF( force_read_thompson ) THEN
1254 write(0,*) "Error reading "//freeze_h2o_file//" Aborting because force_read_thompson is .true."
1255 return
1256 ENDIF
1257 CLOSE(63)
1258 ELSE
1259 IF( force_read_thompson ) THEN
1260 write(0,*) "Error opening "//freeze_h2o_file//" Aborting because force_read_thompson is .true."
1261 return
1262 ENDIF
1263 ENDIF
1264 ELSE
1265 INQUIRE(63,opened=lopen)
1266 IF (lopen) THEN
1267 CLOSE(63)
1268 ENDIF
1269 ENDIF
1270 ELSE
1271 IF( force_read_thompson ) THEN
1272 write(0,*) "Non-existent "//freeze_h2o_file//" Aborting because force_read_thompson is .true."
1273 return
1274 ENDIF
1275 ENDIF
1276
1277 IF (.NOT. good .EQ. 1 ) THEN
1278 if (thompson_table_writer) then
1279 write_thompson_tables = .true.
1280 write(0,*) "ThompMP: computing freezeH2O"
1281 endif
1282
1283 orho_w = 1./rho_w2
1284
1285 do n2 = 1, nbr
1286 massr(n2) = am_r*dr(n2)**bm_r
1287 enddo
1288 do n = 1, nbc
1289 massc(n) = am_r*dc(n)**bm_r
1290 enddo
1291
1292 !..Freeze water (smallest drops become cloud ice, otherwise graupel).
1293 do m = 1, ntb_in
1294 t_adjust = max(-3.0, min(3.0 - alog10(nt_in(m)), 3.0))
1295 do k = 1, 45
1296 ! print*, ' Freezing water for temp = ', -k
1297 texp = dexp( dfloat(k) - t_adjust*1.0d0 ) - 1.0d0
1298 !$OMP PARALLEL DO SCHEDULE(dynamic) num_threads(threads) &
1299 !$OMP PRIVATE(j,i,lam_exp,lamr,N0_r,sum1,sum2,sumn1,sumn2,n2,N_r,vol,prob)
1300 do j = 1, ntb_r1
1301 do i = 1, ntb_r
1302 lam_exp = (n0r_exp(j)*am_r*crg(1)/r_r(i))**ore1
1303 lamr = lam_exp * (crg(3)*org2*org1)**obmr
1304 n0_r = n0r_exp(j)/(crg(2)*lam_exp) * lamr**cre(2)
1305 sum1 = 0.0d0
1306 sum2 = 0.0d0
1307 sumn1 = 0.0d0
1308 sumn2 = 0.0d0
1309 do n2 = nbr, 1, -1
1310 n_r = n0_r*dr(n2)**mu_r*dexp(-lamr*dr(n2))*dtr(n2)
1311 vol = massr(n2)*orho_w
1312 prob = max(0.0d0, 1.0d0 - dexp(-120.0d0*vol*5.2d-4 * texp))
1313 if (massr(n2) .lt. xm0g) then
1314 sumn1 = sumn1 + prob*n_r
1315 sum1 = sum1 + prob*n_r*massr(n2)
1316 else
1317 sumn2 = sumn2 + prob*n_r
1318 sum2 = sum2 + prob*n_r*massr(n2)
1319 endif
1320 if ((sum1+sum2).ge.r_r(i)) EXIT
1321 enddo
1322 tpi_qrfz(i,j,k,m) = sum1
1323 tni_qrfz(i,j,k,m) = sumn1
1324 tpg_qrfz(i,j,k,m) = sum2
1325 tnr_qrfz(i,j,k,m) = sumn2
1326 enddo
1327 enddo
1328 !$OMP END PARALLEL DO
1329
1330 !$OMP PARALLEL DO SCHEDULE(dynamic) num_threads(threads) &
1331 !$OMP PRIVATE(j,i,nu_c,lamc,N0_c,sum1,sumn2,vol,prob,N_c)
1332 do j = 1, nbc
1333 nu_c = min(15, nint(1000.e6/t_nc(j)) + 2)
1334 do i = 1, ntb_c
1335 lamc = (t_nc(j)*am_r* ccg(2,nu_c) * ocg1(nu_c) / r_c(i))**obmr
1336 n0_c = t_nc(j)*ocg1(nu_c) * lamc**cce(1,nu_c)
1337 sum1 = 0.0d0
1338 sumn2 = 0.0d0
1339 do n = nbc, 1, -1
1340 vol = massc(n)*orho_w
1341 prob = max(0.0d0, 1.0d0 - dexp(-120.0d0*vol*5.2d-4 * texp))
1342 n_c = n0_c*dc(n)**nu_c*exp(-lamc*dc(n))*dtc(n)
1343 sumn2 = min(t_nc(j), sumn2 + prob*n_c)
1344 sum1 = sum1 + prob*n_c*massc(n)
1345 if (sum1 .ge. r_c(i)) EXIT
1346 enddo
1347 tpi_qcfz(i,j,k,m) = sum1
1348 tni_qcfz(i,j,k,m) = sumn2
1349 enddo
1350 enddo
1351 !$OMP END PARALLEL DO
1352 enddo
1353 enddo
1354
1355 IF ( write_thompson_tables ) THEN
1356 write(0,*) "Writing "//freeze_h2o_file//" in Tempo MP init"
1357 OPEN(63,file=freeze_h2o_file,form="unformatted",err=9234)
1358 WRITE(63,err=9234)tpi_qrfz
1359 WRITE(63,err=9234)tni_qrfz
1360 WRITE(63,err=9234)tpg_qrfz
1361 WRITE(63,err=9234)tnr_qrfz
1362 WRITE(63,err=9234)tpi_qcfz
1363 WRITE(63,err=9234)tni_qcfz
1364 CLOSE(63)
1365 RETURN ! ----- RETURN
13669234 CONTINUE
1367 write(0,*) "Error writing "//freeze_h2o_file
1368 return
1369 ENDIF
1370 ENDIF
1371
1372 end subroutine freezeh2o_par
1373
1374#endif
1375
1376 !=================================================================================================================
1377 !..Compute _radiation_ effective radii of cloud water, ice, and snow.
1378 !.. These are entirely consistent with microphysics assumptions, not
1379 !.. constant or otherwise ad hoc as is internal to most radiation
1380 !.. schemes. Since only the smallest snowflakes should impact
1381 !.. radiation, compute from first portion of complicated Field number
1382 !.. distribution, not the second part, which is the larger sizes.
1383
1384 subroutine calc_effectrad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, &
1385 & re_qc1d, re_qi1d, re_qs1d, kts, kte, lsml, configs)
1386
1387 IMPLICIT NONE
1388
1389 !..Sub arguments
1390 INTEGER, INTENT(IN):: kts, kte
1391 REAL, DIMENSION(kts:kte), INTENT(IN):: &
1392 & t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d
1393 REAL, DIMENSION(kts:kte), INTENT(INOUT):: re_qc1d, re_qi1d, re_qs1d
1394 type(ty_tempo_cfg), intent(in) :: configs
1395 integer, intent(in), optional :: lsml
1396
1397 !..Local variables
1398 INTEGER:: k
1399 REAL, DIMENSION(kts:kte):: rho, rc, nc, ri, ni, rs
1400 REAL:: smo2, smob, smoc
1401 REAL:: tc0, loga_, a_, b_
1402 DOUBLE PRECISION:: lamc, lami
1403 LOGICAL:: has_qc, has_qi, has_qs
1404 INTEGER:: inu_c
1405 real, dimension(15), parameter:: g_ratio = (/24,60,120,210,336, &
1406 & 504,720,990,1320,1716,2184,2730,3360,4080,4896/)
1407
1408 has_qc = .false.
1409 has_qi = .false.
1410 has_qs = .false.
1411
1412 re_qc1d = 0.
1413 re_qi1d = 0.
1414 re_qs1d = 0.
1415
1416 do k = kts, kte
1417 rho(k) = 0.622*p1d(k)/(r*t1d(k)*(qv1d(k)+0.622))
1418 rc(k) = max(r1, qc1d(k)*rho(k))
1419 nc(k) = max(2., min(nc1d(k)*rho(k), nt_c_max))
1420 if (.not. (configs%aerosol_aware .or. merra2_aerosol_aware)) then
1421 nc(k) = nt_c
1422 if (present(lsml)) then
1423 if( lsml == 1) then
1424 nc(k) = nt_c_l
1425 else
1426 nc(k) = nt_c_o
1427 endif
1428 endif
1429 endif
1430 if (rc(k).gt.r1 .and. nc(k).gt.r2) has_qc = .true.
1431 ri(k) = max(r1, qi1d(k)*rho(k))
1432 ni(k) = max(r2, ni1d(k)*rho(k))
1433 if (ri(k).gt.r1 .and. ni(k).gt.r2) has_qi = .true.
1434 rs(k) = max(r1, qs1d(k)*rho(k))
1435 if (rs(k).gt.r1) has_qs = .true.
1436 enddo
1437
1438 if (has_qc) then
1439 do k = kts, kte
1440 if (rc(k).le.r1 .or. nc(k).le.r2) cycle
1441 if (nc(k).lt.100) then
1442 inu_c = 15
1443 elseif (nc(k).gt.1.e10) then
1444 inu_c = 2
1445 else
1446 inu_c = min(15, nint(1000.e6/nc(k)) + 2)
1447 endif
1448 lamc = (nc(k)*am_r*g_ratio(inu_c)/rc(k))**obmr
1449 re_qc1d(k) = max(2.51e-6, min(sngl(0.5d0 * dble(3.+inu_c)/lamc), 50.e-6))
1450 enddo
1451 endif
1452
1453 if (has_qi) then
1454 do k = kts, kte
1455 if (ri(k).le.r1 .or. ni(k).le.r2) cycle
1456 lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi
1457 re_qi1d(k) = max(2.51e-6, min(sngl(0.5d0 * dble(3.+mu_i)/lami), 125.e-6))
1458 enddo
1459 endif
1460
1461 if (has_qs) then
1462 do k = kts, kte
1463 if (rs(k).le.r1) cycle
1464 tc0 = min(-0.1, t1d(k)-273.15)
1465 smob = rs(k)*oams
1466
1467 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
1468 !.. then we must compute actual 2nd moment and use as reference.
1469 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
1470 smo2 = smob
1471 else
1472 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
1473 & + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
1474 & + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
1475 & + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
1476 & + sa(10)*bm_s*bm_s*bm_s
1477 a_ = 10.0**loga_
1478 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
1479 & + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
1480 & + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
1481 & + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
1482 & + sb(10)*bm_s*bm_s*bm_s
1483 smo2 = (smob/a_)**(1./b_)
1484 endif
1485 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
1486 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
1487 & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
1488 & + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
1489 & + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
1490 & + sa(10)*cse(1)*cse(1)*cse(1)
1491 a_ = 10.0**loga_
1492 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
1493 & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
1494 & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
1495 & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
1496 smoc = a_ * smo2**b_
1497 re_qs1d(k) = max(5.01e-6, min(0.5*(smoc/smob), 999.e-6))
1498 enddo
1499 endif
1500
1501 end subroutine calc_effectrad
1502
1503 !=================================================================================================================
1504 !..Compute radar reflectivity assuming 10 cm wavelength radar and using
1505 !.. Rayleigh approximation. Only complication is melted snow/graupel
1506 !.. which we treat as water-coated ice spheres and use Uli Blahak's
1507 !.. library of routines. The meltwater fraction is simply the amount
1508 !.. of frozen species remaining from what initially existed at the
1509 !.. melting level interface.
1510
1511 subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, ng1d, qb1d, &
1512 t1d, p1d, dBZ, kts, kte, ii, jj, configs, rand1, melti, &
1513 vt_dBZ, first_time_step)
1514
1515 IMPLICIT NONE
1516
1517 !..Sub arguments
1518 INTEGER, INTENT(IN):: kts, kte, ii, jj
1519 REAL, OPTIONAL, INTENT(IN):: rand1
1520 REAL, DIMENSION(kts:kte), INTENT(IN):: &
1521 qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, ng1d, qb1d, t1d, p1d
1522 REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
1523 REAL, DIMENSION(kts:kte), OPTIONAL, INTENT(INOUT):: vt_dBZ
1524 LOGICAL, OPTIONAL, INTENT(IN) :: first_time_step
1525
1526 type(ty_tempo_cfg), intent(in) :: configs
1527
1528 !..Local variables
1529 LOGICAL :: do_vt_dBZ
1530 LOGICAL :: allow_wet_graupel
1531 LOGICAL :: allow_wet_snow
1532 REAL, DIMENSION(kts:kte):: temp, pres, qv, rho, rhof
1533 REAL, DIMENSION(kts:kte):: rc, rr, nr, rs, rg, ng, rb
1534 INTEGER, DIMENSION(kts:kte):: idx_bg
1535
1536 DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g
1537 REAL, DIMENSION(kts:kte):: mvd_r
1538 REAL, DIMENSION(kts:kte):: smob, smo2, smoc, smoz
1539 REAL:: oM3, M0, Mrat, slam1, slam2, xDs
1540 REAL:: ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts
1541 REAL:: vtr_dbz_wt, vts_dbz_wt, vtg_dbz_wt
1542
1543 REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
1544
1545 DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamr, lamg
1546 REAL:: a_, b_, loga_, tc0, SR
1547 DOUBLE PRECISION:: fmelt_s, fmelt_g
1548
1549 INTEGER:: i, k, k_0, kbot, n
1550 LOGICAL, OPTIONAL, INTENT(IN):: melti
1551 LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
1552
1553 DOUBLE PRECISION:: cback, x, eta, f_d
1554 REAL:: xslw1, ygra1, zans1
1555
1556 !+---+
1557 if (present(vt_dbz) .and. present(first_time_step)) then
1558 do_vt_dbz = .true.
1559 if (first_time_step) then
1560 ! no bright banding, to be consistent with hydrometeor retrieval in GSI
1561 allow_wet_snow = .false.
1562 else
1563 allow_wet_snow = .true.
1564 endif
1565 allow_wet_graupel = .false.
1566 else
1567 do_vt_dbz = .false.
1568 allow_wet_snow = .true.
1569 allow_wet_graupel = .false.
1570 endif
1571
1572 do k = kts, kte
1573 dbz(k) = -35.0
1574 enddo
1575
1576 !+---+-----------------------------------------------------------------+
1577 !..Put column of data into local arrays.
1578 !+---+-----------------------------------------------------------------+
1579 do k = kts, kte
1580 temp(k) = t1d(k)
1581 qv(k) = max(1.e-10, qv1d(k))
1582 pres(k) = p1d(k)
1583 rho(k) = 0.622*pres(k)/(r*temp(k)*(qv(k)+0.622))
1584 rhof(k) = sqrt(rho_not/rho(k))
1585 rc(k) = max(r1, qc1d(k)*rho(k))
1586 if (qr1d(k) .gt. r1) then
1587 rr(k) = qr1d(k)*rho(k)
1588 nr(k) = max(r2, nr1d(k)*rho(k))
1589 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
1590 ilamr(k) = 1./lamr
1591 n0_r(k) = nr(k)*org2*lamr**cre(2)
1592 mvd_r(k) = (3.0 + mu_r + 0.672) * ilamr(k)
1593 l_qr(k) = .true.
1594 else
1595 rr(k) = r1
1596 nr(k) = r1
1597 mvd_r(k) = 50.e-6
1598 l_qr(k) = .false.
1599 endif
1600 if (qs1d(k) .gt. r2) then
1601 rs(k) = qs1d(k)*rho(k)
1602 l_qs(k) = .true.
1603 else
1604 rs(k) = r1
1605 l_qs(k) = .false.
1606 endif
1607 if (qg1d(k) .gt. r2) then
1608 rg(k) = qg1d(k)*rho(k)
1609 ng(k) = max(r2, ng1d(k)*rho(k))
1610 rb(k) = max(qg1d(k)/rho_g(nrhg), qb1d(k))
1611 rb(k) = min(qg1d(k)/rho_g(1), rb(k))
1612 idx_bg(k) = max(1,min(nint(qg1d(k)/rb(k) *0.01)+1,nrhg))
1613 if (.not. configs%hail_aware) idx_bg(k) = idx_bg1
1614 l_qg(k) = .true.
1615 else
1616 rg(k) = r1
1617 ng(k) = r2
1618 idx_bg(k) = idx_bg1
1619 l_qg(k) = .false.
1620 endif
1621 enddo
1622
1623 !+---+-----------------------------------------------------------------+
1624 !..Calculate y-intercept, slope, and useful moments for snow.
1625 !+---+-----------------------------------------------------------------+
1626 do k = kts, kte
1627 smo2(k) = 0.
1628 smob(k) = 0.
1629 smoc(k) = 0.
1630 smoz(k) = 0.
1631 enddo
1632 if (any(l_qs .eqv. .true.)) then
1633 do k = kts, kte
1634 if (.not. l_qs(k)) cycle
1635 tc0 = min(-0.1, temp(k)-273.15)
1636 smob(k) = rs(k)*oams
1637
1638 !..All other moments based on reference, 2nd moment. If bm_s.ne.2,
1639 !.. then we must compute actual 2nd moment and use as reference.
1640 if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then
1641 smo2(k) = smob(k)
1642 else
1643 loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s &
1644 & + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 &
1645 & + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s &
1646 & + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 &
1647 & + sa(10)*bm_s*bm_s*bm_s
1648 a_ = 10.0**loga_
1649 b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s &
1650 & + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 &
1651 & + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s &
1652 & + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 &
1653 & + sb(10)*bm_s*bm_s*bm_s
1654 smo2(k) = (smob(k)/a_)**(1./b_)
1655 endif
1656
1657 !..Calculate bm_s+1 (th) moment. Useful for diameter calcs.
1658 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
1659 & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
1660 & + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
1661 & + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
1662 & + sa(10)*cse(1)*cse(1)*cse(1)
1663 a_ = 10.0**loga_
1664 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
1665 & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
1666 & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
1667 & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
1668 smoc(k) = a_ * smo2(k)**b_
1669
1670 !..Calculate bm_s*2 (th) moment. Useful for reflectivity.
1671 loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(3) &
1672 & + sa(4)*tc0*cse(3) + sa(5)*tc0*tc0 &
1673 & + sa(6)*cse(3)*cse(3) + sa(7)*tc0*tc0*cse(3) &
1674 & + sa(8)*tc0*cse(3)*cse(3) + sa(9)*tc0*tc0*tc0 &
1675 & + sa(10)*cse(3)*cse(3)*cse(3)
1676 a_ = 10.0**loga_
1677 b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(3) + sb(4)*tc0*cse(3) &
1678 & + sb(5)*tc0*tc0 + sb(6)*cse(3)*cse(3) &
1679 & + sb(7)*tc0*tc0*cse(3) + sb(8)*tc0*cse(3)*cse(3) &
1680 & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(3)*cse(3)*cse(3)
1681 smoz(k) = a_ * smo2(k)**b_
1682 enddo
1683 endif
1684
1685 !+---+-----------------------------------------------------------------+
1686 !..Calculate y-intercept, slope values for graupel.
1687 !+---+-----------------------------------------------------------------+
1688 if (any(l_qg .eqv. .true.)) then
1689 do k = kte, kts, -1
1690 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg
1691 ilamg(k) = 1./lamg
1692 n0_g(k) = ng(k)*ogg2*lamg**cge(2,1)
1693 enddo
1694 else
1695 ilamg(:) = 0.
1696 n0_g(:) = 0.
1697 endif
1698
1699 !+---+-----------------------------------------------------------------+
1700 !..Locate K-level of start of melting (k_0 is level above).
1701 !+---+-----------------------------------------------------------------+
1702 k_0 = kts
1703 if (present(melti)) then
1704 if ( melti ) then
1705 k_loop:do k = kte-1, kts, -1
1706 if ((temp(k).gt.273.15) .and. l_qr(k) &
1707 & .and. (l_qs(k+1).or.l_qg(k+1)) ) then
1708 k_0 = max(k+1, k_0)
1709 EXIT k_loop
1710 endif
1711 enddo k_loop
1712 endif
1713 endif
1714 !+---+-----------------------------------------------------------------+
1715 !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
1716 !.. and non-water-coated snow and graupel when below freezing are
1717 !.. simple. Integrations of m(D)*m(D)*N(D)*dD.
1718 !+---+-----------------------------------------------------------------+
1719
1720 do k = kts, kte
1721 ze_rain(k) = 1.e-22
1722 ze_snow(k) = 1.e-22
1723 ze_graupel(k) = 1.e-22
1724 if (l_qr(k)) ze_rain(k) = n0_r(k)*crg(4)*ilamr(k)**cre(4)
1725 if (l_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/pi)*(6.0/pi) &
1726 & * (am_s/900.0)*(am_s/900.0)*smoz(k)
1727 if (l_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/pi)*(6.0/pi) &
1728 & * (am_g(idx_bg(k))/900.0)*(am_g(idx_bg(k))/900.0) &
1729 & * n0_g(k)*cgg(4,1)*ilamg(k)**cge(4,1)
1730 enddo
1731
1732 !+---+-----------------------------------------------------------------+
1733 !..Special case of melting ice (snow/graupel) particles. Assume the
1734 !.. ice is surrounded by the liquid water. Fraction of meltwater is
1735 !.. extremely simple based on amount found above the melting level.
1736 !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
1737 !.. routines).
1738 !+---+-----------------------------------------------------------------+
1739 if (present(melti)) then
1740 if (.not. iiwarm .and. melti .and. k_0.ge.2) then
1741 do k = k_0-1, kts, -1
1742
1743 !..Reflectivity contributed by melting snow
1744 if (allow_wet_snow .and. l_qs(k) .and. l_qs(k_0) ) then
1745 sr = max(0.01, min(1.0 - rs(k)/(rs(k) + rr(k)), 0.99))
1746 fmelt_s = dble(sr*sr)
1747 eta = 0.d0
1748 om3 = 1./smoc(k)
1749 m0 = (smob(k)*om3)
1750 mrat = smob(k)*m0*m0*m0
1751 slam1 = m0 * lam0
1752 slam2 = m0 * lam1
1753 do n = 1, nrbins
1754 x = am_s * xxds(n)**bm_s
1755 call rayleigh_soak_wetgraupel (x, dble(ocms), dble(obms), &
1756 & fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
1757 & cback, mixingrulestring_s, matrixstring_s, &
1758 & inclusionstring_s, hoststring_s, &
1759 & hostmatrixstring_s, hostinclusionstring_s)
1760 f_d = mrat*(kap0*dexp(-slam1*xxds(n)) &
1761 & + kap1*(m0*xxds(n))**mu_s * dexp(-slam2*xxds(n)))
1762 eta = eta + f_d * cback * simpson(n) * xdts(n)
1763 enddo
1764 ze_snow(k) = sngl(lamda4 / (pi5 * k_w) * eta)
1765 endif
1766
1767 !..Reflectivity contributed by melting graupel
1768 if (allow_wet_graupel .and. l_qg(k) .and. l_qg(k_0) ) then
1769 sr = max(0.01, min(1.0 - rg(k)/(rg(k) + rr(k)), 0.99))
1770 fmelt_g = dble(sr*sr)
1771 eta = 0.d0
1772 lamg = 1./ilamg(k)
1773 do n = 1, nrbins
1774 x = am_g(idx_bg(k)) * xxdg(n)**bm_g
1775 call rayleigh_soak_wetgraupel (x, dble(ocmg(idx_bg(k))), dble(obmg), &
1776 & fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
1777 & cback, mixingrulestring_g, matrixstring_g, &
1778 & inclusionstring_g, hoststring_g, &
1779 & hostmatrixstring_g, hostinclusionstring_g)
1780 f_d = n0_g(k)*xxdg(n)**mu_g * dexp(-lamg*xxdg(n))
1781 eta = eta + f_d * cback * simpson(n) * xdtg(n)
1782 enddo
1783 ze_graupel(k) = sngl(lamda4 / (pi5 * k_w) * eta)
1784 endif
1785
1786 enddo
1787 endif
1788 endif
1789
1790 do k = kte, kts, -1
1791 dbz(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
1792 enddo
1793
1794 !..Reflectivity-weighted terminal velocity (snow, rain, graupel, mix).
1795 if (do_vt_dbz) then
1796 do k = kte, kts, -1
1797 vt_dbz(k) = 1.e-3
1798 if (rs(k).gt.r2) then
1799 mrat = smob(k) / smoc(k)
1800 ils1 = 1./(mrat*lam0 + fv_s)
1801 ils2 = 1./(mrat*lam1 + fv_s)
1802 t1_vts = kap0*csg(5)*ils1**cse(5)
1803 t2_vts = kap1*mrat**mu_s*csg(11)*ils2**cse(11)
1804 ils1 = 1./(mrat*lam0)
1805 ils2 = 1./(mrat*lam1)
1806 t3_vts = kap0*csg(6)*ils1**cse(6)
1807 t4_vts = kap1*mrat**mu_s*csg(12)*ils2**cse(12)
1808 vts_dbz_wt = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts)
1809 if (temp(k).ge.273.15 .and. temp(k).lt.275.15) then
1810 vts_dbz_wt = vts_dbz_wt*1.5
1811 elseif (temp(k).ge.275.15) then
1812 vts_dbz_wt = vts_dbz_wt*2.0
1813 endif
1814 else
1815 vts_dbz_wt = 1.e-3
1816 endif
1817
1818 if (rr(k).gt.r1) then
1819 lamr = 1./ilamr(k)
1820 vtr_dbz_wt = rhof(k)*av_r*crg(13)*(lamr+fv_r)**(-cre(13)) &
1821 / (crg(4)*lamr**(-cre(4)))
1822 else
1823 vtr_dbz_wt = 1.e-3
1824 endif
1825
1826 if (rg(k).gt.r2) then
1827 lamg = 1./ilamg(k)
1828 vtg_dbz_wt = rhof(k)*av_g(idx_bg(k))*cgg(5,idx_bg(k))*lamg**(-cge(5,idx_bg(k))) &
1829 & / (cgg(4,1)*lamg**(-cge(4,1)))
1830 else
1831 vtg_dbz_wt = 1.e-3
1832 endif
1833
1834 ! if (rg(k).gt.R2) then
1835 ! lamg = 1./ilamg(k)
1836 ! vtg_dbz_wt = rhof(k)*av_g*cgg(5)*lamg**(-cge(5)) &
1837 ! / (cgg(4)*lamg**(-cge(4)))
1838 ! else
1839 ! vtg_dbz_wt = 1.E-3
1840 ! endif
1841
1842 vt_dbz(k) = (vts_dbz_wt*ze_snow(k) + vtr_dbz_wt*ze_rain(k) &
1843 + vtg_dbz_wt*ze_graupel(k)) &
1844 / (ze_rain(k)+ze_snow(k)+ze_graupel(k))
1845 enddo
1846 endif
1847
1848 end subroutine calc_refl10cm
1849
1850 !=================================================================================================================
1851 !..Compute max hail size aloft and at the ground (both 2D fields)
1852
1853 subroutine hail_size_diagnostics(kts, kte, qg1d, ng1d, qb1d, t1d, p1d, qv1d, qg_max_diam1d, configs)
1854
1855 implicit none
1856
1857 integer, intent(in) :: kts, kte
1858 real(wp), dimension(kts:kte), intent(in) :: qg1d, ng1d, qb1d, t1d, p1d, qv1d
1859 real(wp), dimension(kts:kte), intent(out) :: qg_max_diam1d
1860 type(ty_tempo_cfg), intent(in) :: configs
1861
1862 ! local variables
1863 real(wp), dimension(kts:kte) :: rho, rg, ng, rb
1864 integer, dimension(kts:kte) :: idx_bg
1865 real(dp) :: lamg, N0_g, f_d, sum_nh, sum_t, hail_max
1866 integer :: k, n
1867 integer, parameter :: nhbins = 50
1868 real(dp), dimension(nhbins):: hbins, dhbins
1869 real(dp), parameter :: lowbin = 500.e-6
1870 real(dp), parameter :: highbin = 0.075
1871 real(dp), parameter :: threshold_conc = 0.0005
1872
1873 ! Mass distribution (faster calculation, more spread, noisier field)
1874 ! do k = kts, kte
1875 ! qg_max_diam1d(k) = 0.
1876 ! if(qg1d(k) >= 1.e-6) then
1877 ! rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622))
1878 ! rg(k) = qg1d(k)*rho(k)
1879 ! ng(k) = max(R2, ng1d(k)*rho(k))
1880 ! rb(k) = max(qg1d(k)/rho_g(nrhg), qb1d(k))
1881 ! rb(k) = min(qg1d(k)/rho_g(1), rb(k))
1882 ! idx_bg(k) = max(1,min(nint(qg1d(k)/rb(k) *0.01)+1,nrhg))
1883 ! if (.not. configs%hail_aware) idx_bg(k) = idx_bg1
1884 ! if (rho_g(idx_bg(k)) < 550.) cycle
1885 ! lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg
1886 ! qg_max_diam1d(k) = 1000. * 10.05 / lamg ! For graupel, use the 99th percent of mass distribution
1887 ! endif
1888 ! enddo
1889
1890 ! Binned number distribution method
1891 call create_bins(numbins=nhbins, lowbin=lowbin, highbin=highbin, bins=hbins, deltabins=dhbins)
1892
1893 do k = kts, kte
1894 qg_max_diam1d(k) = 0.
1895 if(qg1d(k) >= 1.e-6) then
1896 rho(k) = 0.622*p1d(k)/(r*t1d(k)*(qv1d(k)+0.622))
1897 rg(k) = qg1d(k)*rho(k)
1898 ng(k) = max(r2, ng1d(k)*rho(k))
1899 rb(k) = max(qg1d(k)/rho_g(nrhg), qb1d(k))
1900 rb(k) = min(qg1d(k)/rho_g(1), rb(k))
1901 idx_bg(k) = max(1,min(nint(qg1d(k)/rb(k) *0.01)+1,nrhg))
1902 if (.not. configs%hail_aware) idx_bg(k) = idx_bg1
1903 if (rho_g(idx_bg(k)) < 350.) cycle
1904
1905 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg
1906 n0_g = ng(k)*ogg2*lamg**cge(2,1)
1907
1908 sum_nh = 0.
1909 sum_t = 0.
1910 do n = nhbins, 1, -1
1911 f_d = n0_g*hbins(n)**mu_g * exp(-lamg*hbins(n)) * dhbins(n)
1912 sum_nh = sum_nh + f_d
1913 if (sum_nh > threshold_conc) exit
1914 sum_t = sum_nh
1915 enddo
1916
1917 if (n >= nhbins) then
1918 hail_max = hbins(nhbins)
1919 elseif (hbins(n+1) .gt. 1.e-3) then
1920 hail_max = hbins(n) - (sum_nh-threshold_conc)/(sum_nh-sum_t) * (hbins(n)-hbins(n+1))
1921 else
1922 hail_max = 1.e-4
1923 endif
1924 qg_max_diam1d(k) = 1000. * hail_max ! convert to mm
1925 endif
1926 enddo
1927
1928 end subroutine hail_size_diagnostics
1929
1930 !=================================================================================================================
1931
1932 elemental subroutine make_hydrometeor_number_concentrations(qc, qr, qi, nwfa, temp, rhoa, nc, nr, ni)
1933 implicit none
1934
1935 real, intent(in) :: qc, qr, qi, nwfa, temp, rhoa
1936 real, intent(inout) :: nc, nr, ni
1937
1938 end subroutine make_hydrometeor_number_concentrations
1939
1940 !=================================================================================================================
1941
1947 elemental real function make_icenumber (Q_ice, temp)
1948
1949 !IMPLICIT NONE
1950 REAL, PARAMETER:: ice_density = 890.0
1951 REAL, PARAMETER:: pi = 3.1415926536
1952 real, intent(in):: q_ice, temp
1953 integer idx_rei
1954 real corr, reice, deice
1955 double precision lambda
1956
1957 !+---+-----------------------------------------------------------------+
1958 !..Table of lookup values of radiative effective radius of ice crystals
1959 !.. as a function of Temperature from -94C to 0C. Taken from WRF RRTMG
1960 !.. radiation code where it is attributed to Jon Egill Kristjansson
1961 !.. and coauthors.
1962 !+---+-----------------------------------------------------------------+
1963
1964 !real retab(95)
1965 !data retab / &
1966 ! 5.92779, 6.26422, 6.61973, 6.99539, 7.39234, &
1967 ! 7.81177, 8.25496, 8.72323, 9.21800, 9.74075, 10.2930, &
1968 ! 10.8765, 11.4929, 12.1440, 12.8317, 13.5581, 14.2319, &
1969 ! 15.0351, 15.8799, 16.7674, 17.6986, 18.6744, 19.6955, &
1970 ! 20.7623, 21.8757, 23.0364, 24.2452, 25.5034, 26.8125, &
1971 ! 27.7895, 28.6450, 29.4167, 30.1088, 30.7306, 31.2943, &
1972 ! 31.8151, 32.3077, 32.7870, 33.2657, 33.7540, 34.2601, &
1973 ! 34.7892, 35.3442, 35.9255, 36.5316, 37.1602, 37.8078, &
1974 ! 38.4720, 39.1508, 39.8442, 40.5552, 41.2912, 42.0635, &
1975 ! 42.8876, 43.7863, 44.7853, 45.9170, 47.2165, 48.7221, &
1976 ! 50.4710, 52.4980, 54.8315, 57.4898, 60.4785, 63.7898, &
1977 ! 65.5604, 71.2885, 75.4113, 79.7368, 84.2351, 88.8833, &
1978 ! 93.6658, 98.5739, 103.603, 108.752, 114.025, 119.424, &
1979 ! 124.954, 130.630, 136.457, 142.446, 148.608, 154.956, &
1980 ! 161.503, 168.262, 175.248, 182.473, 189.952, 197.699, &
1981 ! 205.728, 214.055, 222.694, 231.661, 240.971, 250.639/
1982 real, dimension(95), parameter:: retab = (/ &
1983 5.92779, 6.26422, 6.61973, 6.99539, 7.39234, &
1984 7.81177, 8.25496, 8.72323, 9.21800, 9.74075, 10.2930, &
1985 10.8765, 11.4929, 12.1440, 12.8317, 13.5581, 14.2319, &
1986 15.0351, 15.8799, 16.7674, 17.6986, 18.6744, 19.6955, &
1987 20.7623, 21.8757, 23.0364, 24.2452, 25.5034, 26.8125, &
1988 27.7895, 28.6450, 29.4167, 30.1088, 30.7306, 31.2943, &
1989 31.8151, 32.3077, 32.7870, 33.2657, 33.7540, 34.2601, &
1990 34.7892, 35.3442, 35.9255, 36.5316, 37.1602, 37.8078, &
1991 38.4720, 39.1508, 39.8442, 40.5552, 41.2912, 42.0635, &
1992 42.8876, 43.7863, 44.7853, 45.9170, 47.2165, 48.7221, &
1993 50.4710, 52.4980, 54.8315, 57.4898, 60.4785, 63.7898, &
1994 65.5604, 71.2885, 75.4113, 79.7368, 84.2351, 88.8833, &
1995 93.6658, 98.5739, 103.603, 108.752, 114.025, 119.424, &
1996 124.954, 130.630, 136.457, 142.446, 148.608, 154.956, &
1997 161.503, 168.262, 175.248, 182.473, 189.952, 197.699, &
1998 205.728, 214.055, 222.694, 231.661, 240.971, 250.639 /)
1999
2000 if (q_ice == 0) then
2001 make_icenumber = 0
2002 return
2003 end if
2004
2005 !+---+-----------------------------------------------------------------+
2006 !..From the model 3D temperature field, subtract 179K for which
2007 !.. index value of retab as a start. Value of corr is for
2008 !.. interpolating between neighboring values in the table.
2009 !+---+-----------------------------------------------------------------+
2010
2011 idx_rei = int(temp-179.)
2012 idx_rei = min(max(idx_rei,1),94)
2013 corr = temp - int(temp)
2014 reice = retab(idx_rei)*(1.-corr) + retab(idx_rei+1)*corr
2015 deice = 2.*reice * 1.e-6
2016
2017 !+---+-----------------------------------------------------------------+
2018 !..Now we have the final radiative effective size of ice (as function
2019 !.. of temperature only). This size represents 3rd moment divided by
2020 !.. second moment of the ice size distribution, so we can compute a
2021 !.. number concentration from the mean size and mass mixing ratio.
2022 !.. The mean (radiative effective) diameter is 3./Slope for an inverse
2023 !.. exponential size distribution. So, starting with slope, work
2024 !.. backwords to get number concentration.
2025 !+---+-----------------------------------------------------------------+
2026
2027 lambda = 3.0 / deice
2028 make_icenumber = q_ice * lambda*lambda*lambda / (pi*ice_density)
2029
2030 !+---+-----------------------------------------------------------------+
2031 !..Example1: Common ice size coming from Thompson scheme is about 30 microns.
2032 !.. An example ice mixing ratio could be 0.001 g/kg for a temperature of -50C.
2033 !.. Remember to convert both into MKS units. This gives N_ice=357652 per kg.
2034 !..Example2: Lower in atmosphere at T=-10C matching ~162 microns in retab,
2035 !.. and assuming we have 0.1 g/kg mixing ratio, then N_ice=28122 per kg,
2036 !.. which is 28 crystals per liter of air if the air density is 1.0.
2037 !+---+-----------------------------------------------------------------+
2038
2039 return
2040 end function make_icenumber
2041
2042 !=================================================================================================================
2043 elemental real function make_dropletnumber (Q_cloud, qnwfa)
2044
2045 !IMPLICIT NONE
2046
2047 real, intent(in):: q_cloud, qnwfa
2048
2049 real, parameter:: pi = 3.1415926536
2050 real, parameter:: am_r = pi*1000./6.
2051 real, dimension(15), parameter:: g_ratio = (/24,60,120,210,336, &
2052 & 504,720,990,1320,1716,2184,2730,3360,4080,4896/)
2053 double precision:: lambda, qnc
2054 real:: q_nwfa, x1, xdc
2055 integer:: nu_c
2056
2057 if (q_cloud == 0) then
2058 make_dropletnumber = 0
2059 return
2060 end if
2061
2062 !+---+
2063
2064 q_nwfa = max(99.e6, min(qnwfa,5.e10))
2065 nu_c = max(2, min(nint(2.5e10/q_nwfa), 15))
2066
2067 x1 = max(1., min(q_nwfa*1.e-9, 10.)) - 1.
2068 xdc = (30. - x1*20./9.) * 1.e-6
2069
2070 lambda = (4.0d0 + nu_c) / xdc
2071 qnc = q_cloud / g_ratio(nu_c) * lambda*lambda*lambda / am_r
2072 make_dropletnumber = sngl(qnc)
2073
2074 return
2075 end function make_dropletnumber
2076
2077 !=================================================================================================================
2078 elemental real function make_rainnumber (Q_rain, temp)
2079
2080 IMPLICIT NONE
2081
2082 real, intent(in):: q_rain, temp
2083 double precision:: lambda, n0, qnr
2084 real, parameter:: pi = 3.1415926536
2085 real, parameter:: am_r = pi*1000./6.
2086
2087 if (q_rain == 0) then
2088 make_rainnumber = 0
2089 return
2090 end if
2091
2092 !+---+-----------------------------------------------------------------+
2093 !.. Not thrilled with it, but set Y-intercept parameter to Marshal-Palmer value
2094 !.. that basically assumes melting snow becomes typical rain. However, for
2095 !.. -2C < T < 0C, make linear increase in exponent to attempt to keep
2096 !.. supercooled collision-coalescence (warm-rain) similar to drizzle rather
2097 !.. than bigger rain drops. While this could also exist at T>0C, it is
2098 !.. more difficult to assume it directly from having mass and not number.
2099 !+---+-----------------------------------------------------------------+
2100
2101 n0 = 8.e6
2102
2103 if (temp .le. 271.15) then
2104 n0 = 8.e8
2105 elseif (temp .gt. 271.15 .and. temp.lt.273.15) then
2106 n0 = 8. * 10**(279.15-temp)
2107 endif
2108
2109 lambda = sqrt(sqrt(n0*am_r*6.0/q_rain))
2110 qnr = q_rain / 6.0 * lambda*lambda*lambda / am_r
2111 make_rainnumber = sngl(qnr)
2112
2113 return
2114 end function make_rainnumber
2115
2116!=================================================================================================================
2117 !+---+-----------------------------------------------------------------+
2118 ! THIS FUNCTION CALCULATES THE LIQUID SATURATION VAPOR MIXING RATIO AS
2119 ! A FUNCTION OF TEMPERATURE AND PRESSURE
2120 !
2121 REAL function rslf(p,t)
2122
2123 IMPLICIT NONE
2124 REAL, INTENT(IN):: p, t
2125 REAL:: esl,x
2126 REAL, PARAMETER:: c0= .611583699e03
2127 REAL, PARAMETER:: c1= .444606896e02
2128 REAL, PARAMETER:: c2= .143177157e01
2129 REAL, PARAMETER:: c3= .264224321e-1
2130 REAL, PARAMETER:: c4= .299291081e-3
2131 REAL, PARAMETER:: c5= .203154182e-5
2132 REAL, PARAMETER:: c6= .702620698e-8
2133 REAL, PARAMETER:: c7= .379534310e-11
2134 REAL, PARAMETER:: c8=-.321582393e-13
2135
2136 x=max(-80.,t-273.16)
2137
2138 ! ESL=612.2*EXP(17.67*X/(T-29.65))
2139 esl=c0+x*(c1+x*(c2+x*(c3+x*(c4+x*(c5+x*(c6+x*(c7+x*c8)))))))
2140 esl=min(esl, p*0.15) ! Even with P=1050mb and T=55C, the sat. vap. pres only contributes to ~15% of total pres.
2141 rslf=.622*esl/(p-esl)
2142
2143 ! ALTERNATIVE
2144 ! ; Source: Murphy and Koop, Review of the vapour pressure of ice and
2145 ! supercooled water for atmospheric applications, Q. J. R.
2146 ! Meteorol. Soc (2005), 131, pp. 1539-1565.
2147 ! ESL = EXP(54.842763 - 6763.22 / T - 4.210 * ALOG(T) + 0.000367 * T
2148 ! + TANH(0.0415 * (T - 218.8)) * (53.878 - 1331.22
2149 ! / T - 9.44523 * ALOG(T) + 0.014025 * T))
2150
2151 END FUNCTION rslf
2152 !+---+-----------------------------------------------------------------+
2153 ! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A
2154 ! FUNCTION OF TEMPERATURE AND PRESSURE
2155 !
2156 REAL function rsif(p,t)
2157
2158 IMPLICIT NONE
2159 REAL, INTENT(IN):: p, t
2160 REAL:: esi,x
2161 REAL, PARAMETER:: c0= .609868993e03
2162 REAL, PARAMETER:: c1= .499320233e02
2163 REAL, PARAMETER:: c2= .184672631e01
2164 REAL, PARAMETER:: c3= .402737184e-1
2165 REAL, PARAMETER:: c4= .565392987e-3
2166 REAL, PARAMETER:: c5= .521693933e-5
2167 REAL, PARAMETER:: c6= .307839583e-7
2168 REAL, PARAMETER:: c7= .105785160e-9
2169 REAL, PARAMETER:: c8= .161444444e-12
2170
2171 x=max(-80.,t-273.16)
2172 esi=c0+x*(c1+x*(c2+x*(c3+x*(c4+x*(c5+x*(c6+x*(c7+x*c8)))))))
2173 esi=min(esi, p*0.15)
2174 rsif=.622*esi/max(1.e-4,(p-esi))
2175
2176 ! ALTERNATIVE
2177 ! ; Source: Murphy and Koop, Review of the vapour pressure of ice and
2178 ! supercooled water for atmospheric applications, Q. J. R.
2179 ! Meteorol. Soc (2005), 131, pp. 1539-1565.
2180 ! ESI = EXP(9.550426 - 5723.265/T + 3.53068*ALOG(T) - 0.00728332*T)
2181
2182 END FUNCTION rsif
2183 !=================================================================================================================
2184
2185end module module_mp_tempo_utils
module_mp_tempo_utils::table_ccnact
subroutine table_ccnact(errmess, errflag)
Fill the table of CCN activation data created from parcel model run by Trude Eidhammer with inputs of...
Definition module_mp_tempo_utils.F90:716
module_mp_tempo_utils::make_icenumber
elemental real function make_icenumber(q_ice, temp)
Table of lookup values of radiative effective radius of ice crystals as a function of Temperature fro...
Definition module_mp_tempo_utils.F90:1948
module_mp_radar
This module is more library code whereas the individual microphysics schemes contains specific detail...
Definition module_mp_radar.F90:26
module_mp_tempo_params
Definition module_mp_tempo_params.F90:3
module_mp_tempo_utils
Definition module_mp_tempo_utils.F90:3