CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
module_mp_tempo.F90
1! 3D TEMPO Driver for MPAS
2!=================================================================================================================
3module module_mp_tempo
4
5 use machine, only: wp => kind_phys, sp => kind_sngl_prec, dp => kind_dbl_prec
6 use module_mp_tempo_params
7 use module_mp_tempo_utils, only : create_bins, table_efrw, table_efsw, table_dropevap, &
8 calc_refl10cm, calc_effectrad
9 use module_mp_tempo_main, only : mp_tempo_main
10 ! use mp_radar
11
12 implicit none
13
14contains
15 !=================================================================================================================
16 ! This subroutine handles initialzation of the microphysics scheme including building of lookup tables,
17 ! allocating arrays for the microphysics scheme, and defining gamma function variables.
18
19 ! Input:
20 ! l_mp_tables = .false. to build lookup tables. If l_mp_tables = .true., lookup tables are built.
21
22 ! AAJ No support yet for hail_aware in microphysics driver
23 subroutine tempo_init(l_mp_tables, hail_aware_flag, aerosol_aware_flag)
24
25 ! Input arguments:
26 logical, intent(in) :: l_mp_tables
27 logical, intent(in), optional :: aerosol_aware_flag, hail_aware_flag
28
29 integer, parameter :: open_OK = 0
30 integer, parameter :: num_records = 5
31 integer :: qr_acr_qg_filesize, qr_acr_qg_check, qr_acr_qg_dim1size, qr_acr_qg_dim9size
32 logical :: qr_acr_qg_exists, qr_acr_qg_hailaware_exists
33 integer :: i, j, k, l, m, n
34 integer :: istat
35 logical :: micro_init
36 integer :: mp_unit
37 character(len=132) :: message
38
39 ! Initialize physical constants
40 call mp_tempo_params_init()
41
42 if (present(hail_aware_flag)) then
43 configs%hail_aware = hail_aware_flag
44 else
45 configs%hail_aware = .false.
46 endif
47
48 ! If lookup tables are already built
49 if (l_mp_tables) then
50 ! configs%hail_aware = hail_aware_flag
51 ! write(message, '(L1)') configs%hail_aware
52
53 if (present(aerosol_aware_flag)) then
54 configs%aerosol_aware = aerosol_aware_flag
55 ! write(message, '(L1)') configs%aerosol_aware
56 endif
57 endif
58
59 if (configs%hail_aware) then
60 dimnrhg = nrhg
61 else
62 av_g(idx_bg1) = av_g_old
63 bv_g(idx_bg1) = bv_g_old
64 dimnrhg = nrhg1
65 endif
66
67 micro_init = .false.
68
69 !=================================================================================================================
70 ! Check the qr_acr_qg lookup table to make sure it is compatible with runtime options
71
72 ! If lookup tables are already built
73 ! if (l_mp_tables) then
74
75 ! inquire(file='MP_TEMPO_QRacrQG_DATA.DBL', exist=qr_acr_qg_exists)
76 ! if (qr_acr_qg_exists) then ! Check again that file exists
77
78 ! ! Add check on qr_ac_qg filesize to determine if table includes hail-awareness (dimNRHG=9)
79 ! qr_acr_qg_check = dp * num_records * (dimNRHG * ntb_g1 * ntb_g * ntb_r1 * ntb_r + 1)
80 ! qr_acr_qg_dim1size = dp * num_records * (NRHG1 * ntb_g1 * ntb_g * ntb_r1 * ntb_r + 1)
81 ! qr_acr_qg_dim9size = dp * num_records * (NRHG * ntb_g1 * ntb_g * ntb_r1 * ntb_r + 1)
82
83 ! inquire(file='MP_TEMPO_QRacrQG_DATA.DBL', size=qr_acr_qg_filesize)
84
85 ! if (qr_acr_qg_filesize == qr_acr_qg_dim1size) then
86 ! using_hail_aware_table = .false.
87 ! call physics_message('--- tempo_init() ' // &
88 ! 'Lookup table for qr_acr_qg is not hail aware.')
89 ! dimNRHG = NRHG1
90 ! if (hail_aware_flag) then
91 ! call physics_error_fatal('--- tempo_init() Cannot use hail-aware microphysics ' // &
92 ! 'with non hail-aware qr_acr_qg lookup table. ' // &
93 ! 'Please rebuild table with parameter build_hail_aware_table set to true.')
94 ! endif
95 ! elseif (qr_acr_qg_filesize == qr_acr_qg_dim9size) then
96 ! using_hail_aware_table = .true.
97 ! call physics_message('--- tempo_init() ' // &
98 ! 'Lookup table for qr_acr_qg is hail aware.')
99 ! dimNRHG = NRHG
100 ! else
101 ! using_hail_aware_table = .false.
102 ! if (hail_aware_flag) using_hail_aware_table = .true.
103 ! call physics_message('--- tempo_init() ' // &
104 ! 'Could not determine if lookup table for qr_acr_qg is hail aware based on file size.')
105 ! endif
106 ! endif
107 ! endif
108
109 !=================================================================================================================
110 ! Allocate space for lookup tables (J. Michalakes 2009Jun08).
111 if (.not. allocated(tcg_racg)) then
112 allocate(tcg_racg(ntb_g1,ntb_g,dimnrhg,ntb_r1,ntb_r))
113 micro_init = .true.
114 endif
115
116 ! Rain-graupel (including array above tcg_racg)
117 if (.not. allocated(tmr_racg)) allocate(tmr_racg(ntb_g1,ntb_g,dimnrhg,ntb_r1,ntb_r))
118 if (.not. allocated(tcr_gacr)) allocate(tcr_gacr(ntb_g1,ntb_g,dimnrhg,ntb_r1,ntb_r))
119 if (.not. allocated(tnr_racg)) allocate(tnr_racg(ntb_g1,ntb_g,dimnrhg,ntb_r1,ntb_r))
120 if (.not. allocated(tnr_gacr)) allocate(tnr_gacr(ntb_g1,ntb_g,dimnrhg,ntb_r1,ntb_r))
121
122 ! Rain-snow
123 if (.not. allocated(tcs_racs1)) allocate(tcs_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
124 if (.not. allocated(tmr_racs1)) allocate(tmr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
125 if (.not. allocated(tcs_racs2)) allocate(tcs_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
126 if (.not. allocated(tmr_racs2)) allocate(tmr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
127 if (.not. allocated(tcr_sacr1)) allocate(tcr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
128 if (.not. allocated(tms_sacr1)) allocate(tms_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
129 if (.not. allocated(tcr_sacr2)) allocate(tcr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
130 if (.not. allocated(tms_sacr2)) allocate(tms_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
131 if (.not. allocated(tnr_racs1)) allocate(tnr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
132 if (.not. allocated(tnr_racs2)) allocate(tnr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
133 if (.not. allocated(tnr_sacr1)) allocate(tnr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
134 if (.not. allocated(tnr_sacr2)) allocate(tnr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
135
136 ! Cloud water freezing
137 if (.not. allocated(tpi_qcfz)) allocate(tpi_qcfz(ntb_c,nbc,ntb_t1,ntb_in))
138 if (.not. allocated(tni_qcfz)) allocate(tni_qcfz(ntb_c,nbc,ntb_t1,ntb_in))
139
140 ! Rain freezing
141 if (.not. allocated(tpi_qrfz)) allocate(tpi_qrfz(ntb_r,ntb_r1,ntb_t1,ntb_in))
142 if (.not. allocated(tpg_qrfz)) allocate(tpg_qrfz(ntb_r,ntb_r1,ntb_t1,ntb_in))
143 if (.not. allocated(tni_qrfz)) allocate(tni_qrfz(ntb_r,ntb_r1,ntb_t1,ntb_in))
144 if (.not. allocated(tnr_qrfz)) allocate(tnr_qrfz(ntb_r,ntb_r1,ntb_t1,ntb_in))
145
146 ! Ice growth and conversion to snow
147 if (.not. allocated(tps_iaus)) allocate(tps_iaus(ntb_i,ntb_i1))
148 if (.not. allocated(tni_iaus)) allocate(tni_iaus(ntb_i,ntb_i1))
149 if (.not. allocated(tpi_ide)) allocate(tpi_ide(ntb_i,ntb_i1))
150
151 ! Collision efficiencies
152 if (.not. allocated(t_efrw)) allocate(t_efrw(nbr,nbc))
153 if (.not. allocated(t_efsw)) allocate(t_efsw(nbs,nbc))
154
155 ! Cloud water
156 if (.not. allocated(tnr_rev)) allocate(tnr_rev(nbr,ntb_r1,ntb_r))
157 if (.not. allocated(tpc_wev)) allocate(tpc_wev(nbc,ntb_c,nbc))
158 if (.not. allocated(tnc_wev)) allocate(tnc_wev(nbc,ntb_c,nbc))
159
160 ! CCN
161 if (.not. allocated(tnccn_act)) allocate(tnccn_act(ntb_arc,ntb_arw,ntb_art,ntb_arr,ntb_ark))
162
163 !=================================================================================================================
164 if (micro_init) then
165
166 ! Schmidt number to one-third used numerous times.
167 sc3 = sc**(1./3.)
168
169 ! Compute minimum ice diameter from mass and minimum snow/graupel mass from diameter
170 d0i = (xm0i/am_i)**(1.0/bm_i)
171 xm0s = am_s * d0s**bm_s
172 xm0g = am_g(nrhg) * d0g**bm_g
173
174 ! These constants various exponents and gamma() assoc with cloud, rain, snow, and graupel.
175 do n = 1, 15
176 cce(1,n) = n + 1.
177 cce(2,n) = bm_r + n + 1.
178 cce(3,n) = bm_r + n + 4.
179 cce(4,n) = n + bv_c + 1.
180 cce(5,n) = bm_r + n + bv_c + 1.
181 ccg(1,n) = gamma(cce(1,n))
182 ccg(2,n) = gamma(cce(2,n))
183 ccg(3,n) = gamma(cce(3,n))
184 ccg(4,n) = gamma(cce(4,n))
185 ccg(5,n) = gamma(cce(5,n))
186 ocg1(n) = 1.0 / ccg(1,n)
187 ocg2(n) = 1.0 / ccg(2,n)
188 enddo
189
190 cie(1) = mu_i + 1.
191 cie(2) = bm_i + mu_i + 1.
192 cie(3) = bm_i + mu_i + bv_i + 1.
193 cie(4) = mu_i + bv_i + 1.
194 cie(5) = mu_i + 2.
195 cie(6) = bm_i*0.5 + mu_i + bv_i + 1.
196 cie(7) = bm_i*0.5 + mu_i + 1.
197 cig(1) = gamma(cie(1))
198 cig(2) = gamma(cie(2))
199 cig(3) = gamma(cie(3))
200 cig(4) = gamma(cie(4))
201 cig(5) = gamma(cie(5))
202 cig(6) = gamma(cie(6))
203 cig(7) = gamma(cie(7))
204 oig1 = 1.0 / cig(1)
205 oig2 = 1.0 / cig(2)
206 obmi = 1.0 / bm_i
207
208 cre(1) = bm_r + 1.
209 cre(2) = mu_r + 1.
210 cre(3) = bm_r + mu_r + 1.
211 cre(4) = bm_r*2. + mu_r + 1.
212 cre(5) = mu_r + bv_r + 1.
213 cre(6) = bm_r + mu_r + bv_r + 1.
214 cre(7) = bm_r*0.5 + mu_r + bv_r + 1.
215 cre(8) = bm_r + mu_r + bv_r + 3.
216 cre(9) = mu_r + bv_r + 3.
217 cre(10) = mu_r + 2.
218 cre(11) = 0.5*(bv_r + 5. + 2.*mu_r)
219 cre(12) = bm_r*0.5 + mu_r + 1.
220 cre(13) = bm_r*2. + mu_r + bv_r + 1.
221
222 do n = 1, 13
223 crg(n) = gamma(cre(n))
224 enddo
225
226 obmr = 1.0 / bm_r
227 ore1 = 1.0 / cre(1)
228 org1 = 1.0 / crg(1)
229 org2 = 1.0 / crg(2)
230 org3 = 1.0 / crg(3)
231
232 cse(1) = bm_s + 1.
233 cse(2) = bm_s + 2.
234 cse(3) = bm_s*2.
235 cse(4) = bm_s + bv_s + 1.
236 cse(5) = bm_s*2. + bv_s + 1.
237 cse(6) = bm_s*2. + 1.
238 cse(7) = bm_s + mu_s + 1.
239 cse(8) = bm_s + mu_s + 2.
240 cse(9) = bm_s + mu_s + 3.
241 cse(10) = bm_s + mu_s + bv_s + 1.
242 cse(11) = bm_s*2. + mu_s + bv_s + 1.
243 cse(12) = bm_s*2. + mu_s + 1.
244 cse(13) = bv_s + 2.
245 cse(14) = bm_s + bv_s
246 cse(15) = mu_s + 1.
247 cse(16) = 1.0 + (1.0 + bv_s)/2.
248 cse(17) = bm_s + bv_s + 2.
249
250 do n = 1, 17
251 csg(n) = gamma(cse(n))
252 enddo
253
254 oams = 1.0 / am_s
255 obms = 1.0 / bm_s
256 ocms = oams**obms
257
258 cge(1,:) = bm_g + 1.
259 cge(2,:) = mu_g + 1.
260 cge(3,:) = bm_g + mu_g + 1.
261 cge(4,:) = bm_g*2. + mu_g + 1.
262 cge(10,:) = mu_g + 2.
263 cge(12,:) = bm_g*0.5 + mu_g + 1.
264
265 do m = 1, nrhg
266 cge(5,m) = bm_g*2. + mu_g + bv_g(m) + 1.
267 cge(6,m) = bm_g + mu_g + bv_g(m) + 1.
268 cge(7,m) = bm_g*0.5 + mu_g + bv_g(m) + 1.
269 cge(8,m) = mu_g + bv_g(m) + 1. ! not used
270 cge(9,m) = mu_g + bv_g(m) + 3.
271 cge(11,m) = 0.5*(bv_g(m) + 5. + 2.*mu_g)
272 enddo
273
274 do m = 1, nrhg
275 do n = 1, 12
276 cgg(n,m) = gamma(cge(n,m))
277 enddo
278 enddo
279
280 oamg = 1.0 / am_g
281 obmg = 1.0 / bm_g
282
283 do m = 1, nrhg
284 oamg(m) = 1.0 / am_g(m)
285 ocmg(m) = oamg(m)**obmg
286 enddo
287
288 oge1 = 1.0 / cge(1,1)
289 ogg1 = 1.0 / cgg(1,1)
290 ogg2 = 1.0 / cgg(2,1)
291 ogg3 = 1.0 / cgg(3,1)
292
293 !=================================================================================================================
294 ! Simplify various rate eqns the best we can now.
295
296 ! Rain collecting cloud water and cloud ice
297 t1_qr_qc = pi * 0.25 * av_r * crg(9)
298 t1_qr_qi = pi * 0.25 * av_r * crg(9)
299 t2_qr_qi = pi * 0.25 * am_r*av_r * crg(8)
300
301 ! Graupel collecting cloud water
302 ! t1_qg_qc = PI*.25*av_g * cgg(9)
303
304 ! Snow collecting cloud water
305 t1_qs_qc = pi * 0.25 * av_s
306
307 ! Snow collecting cloud ice
308 t1_qs_qi = pi * 0.25 * av_s
309
310 ! Evaporation of rain; ignore depositional growth of rain.
311 t1_qr_ev = 0.78 * crg(10)
312 t2_qr_ev = 0.308 * sc3 * sqrt(av_r) * crg(11)
313
314 ! Sublimation/depositional growth of snow
315 t1_qs_sd = 0.86
316 t2_qs_sd = 0.28 * sc3 * sqrt(av_s)
317
318 ! Melting of snow
319 t1_qs_me = pi * 4. *c_sqrd * olfus * 0.86
320 t2_qs_me = pi * 4. *c_sqrd * olfus * 0.28 * sc3 * sqrt(av_s)
321
322 ! Sublimation/depositional growth of graupel
323 t1_qg_sd = 0.86 * cgg(10,1)
324 ! t2_qg_sd = 0.28*Sc3*SQRT(av_g) * cgg(11)
325
326 ! Melting of graupel
327 t1_qg_me = pi * 4. * c_cube * olfus * 0.86 * cgg(10,1)
328 ! t2_qg_me = PI*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g) * cgg(11)
329
330
331 ! Constants for helping find lookup table indexes.
332 nic2 = nint(log10(r_c(1)))
333 nii2 = nint(log10(r_i(1)))
334 nii3 = nint(log10(nt_i(1)))
335 nir2 = nint(log10(r_r(1)))
336 nir3 = nint(log10(n0r_exp(1)))
337 nis2 = nint(log10(r_s(1)))
338 nig2 = nint(log10(r_g(1)))
339 nig3 = nint(log10(n0g_exp(1)))
340 niin2 = nint(log10(nt_in(1)))
341
342 ! Create bins of cloud water (from minimum diameter to 100 microns).
343 dc(1) = d0c*1.0_dp
344 dtc(1) = d0c*1.0_dp
345 do n = 2, nbc
346 dc(n) = dc(n-1) + 1.0e-6_dp
347 dtc(n) = (dc(n) - dc(n-1))
348 enddo
349
350 ! Create bins of cloud ice (from min diameter up to 2x min snow size).
351 call create_bins(numbins=nbi, lowbin=d0i*1.0_dp, highbin=d0s*2.0_dp, &
352 bins=di, deltabins=dti)
353
354 ! Create bins of rain (from min diameter up to 5 mm).
355 call create_bins(numbins=nbr, lowbin=d0r*1.0_dp, highbin=0.005_dp, &
356 bins=dr, deltabins=dtr)
357
358 ! Create bins of snow (from min diameter up to 2 cm).
359 call create_bins(numbins=nbs, lowbin=d0s*1.0_dp, highbin=0.02_dp, &
360 bins=ds, deltabins=dts)
361
362 ! Create bins of graupel (from min diameter up to 5 cm).
363 call create_bins(numbins=nbg, lowbin=d0g*1.0_dp, highbin=0.05_dp, &
364 bins=dg, deltabins=dtg)
365
366 ! Create bins of cloud droplet number concentration (1 to 3000 per cc).
367 call create_bins(numbins=nbc, lowbin=1.0_dp, highbin=3000.0_dp, &
368 bins=t_nc)
369 t_nc = t_nc * 1.0e6_dp
370 nic1 = log(t_nc(nbc)/t_nc(1))
371
372 !=================================================================================================================
373 ! Create lookup tables for most costly calculations.
374
375 do m = 1, ntb_r
376 do k = 1, ntb_r1
377 do n = 1, dimnrhg
378 do j = 1, ntb_g
379 do i = 1, ntb_g1
380 tcg_racg(i,j,n,k,m) = 0.0_dp
381 tmr_racg(i,j,n,k,m) = 0.0_dp
382 tcr_gacr(i,j,n,k,m) = 0.0_dp
383 tnr_racg(i,j,n,k,m) = 0.0_dp
384 tnr_gacr(i,j,n,k,m) = 0.0_dp
385 enddo
386 enddo
387 enddo
388 enddo
389 enddo
390
391 do m = 1, ntb_r
392 do k = 1, ntb_r1
393 do j = 1, ntb_t
394 do i = 1, ntb_s
395 tcs_racs1(i,j,k,m) = 0.0_dp
396 tmr_racs1(i,j,k,m) = 0.0_dp
397 tcs_racs2(i,j,k,m) = 0.0_dp
398 tmr_racs2(i,j,k,m) = 0.0_dp
399 tcr_sacr1(i,j,k,m) = 0.0_dp
400 tms_sacr1(i,j,k,m) = 0.0_dp
401 tcr_sacr2(i,j,k,m) = 0.0_dp
402 tms_sacr2(i,j,k,m) = 0.0_dp
403 tnr_racs1(i,j,k,m) = 0.0_dp
404 tnr_racs2(i,j,k,m) = 0.0_dp
405 tnr_sacr1(i,j,k,m) = 0.0_dp
406 tnr_sacr2(i,j,k,m) = 0.0_dp
407 enddo
408 enddo
409 enddo
410 enddo
411
412 do m = 1, ntb_in
413 do k = 1, ntb_t1
414 do j = 1, ntb_r1
415 do i = 1, ntb_r
416 tpi_qrfz(i,j,k,m) = 0.0_dp
417 tni_qrfz(i,j,k,m) = 0.0_dp
418 tpg_qrfz(i,j,k,m) = 0.0_dp
419 tnr_qrfz(i,j,k,m) = 0.0_dp
420 enddo
421 enddo
422 do j = 1, nbc
423 do i = 1, ntb_c
424 tpi_qcfz(i,j,k,m) = 0.0_dp
425 tni_qcfz(i,j,k,m) = 0.0_dp
426 enddo
427 enddo
428 enddo
429 enddo
430
431 do j = 1, ntb_i1
432 do i = 1, ntb_i
433 tps_iaus(i,j) = 0.0_dp
434 tni_iaus(i,j) = 0.0_dp
435 tpi_ide(i,j) = 0.0_dp
436 enddo
437 enddo
438
439 do j = 1, nbc
440 do i = 1, nbr
441 t_efrw(i,j) = 0.0
442 enddo
443 do i = 1, nbs
444 t_efsw(i,j) = 0.0
445 enddo
446 enddo
447
448 do k = 1, ntb_r
449 do j = 1, ntb_r1
450 do i = 1, nbr
451 tnr_rev(i,j,k) = 0.0_dp
452 enddo
453 enddo
454 enddo
455
456 do k = 1, nbc
457 do j = 1, ntb_c
458 do i = 1, nbc
459 tpc_wev(i,j,k) = 0.0_dp
460 tnc_wev(i,j,k) = 0.0_dp
461 enddo
462 enddo
463 enddo
464
465 do m = 1, ntb_ark
466 do l = 1, ntb_arr
467 do k = 1, ntb_art
468 do j = 1, ntb_arw
469 do i = 1, ntb_arc
470 tnccn_act(i,j,k,l,m) = 1.0
471 enddo
472 enddo
473 enddo
474 enddo
475 enddo
476
477 !=================================================================================================================
478 ! Check that the look-up tables are available.
479 if (.not. l_mp_tables) return
480
481 ! Collision efficiency between rain/snow and cloud water.
482 call table_efrw ! => fills t_Efrw
483 call table_efsw ! => fills t_Efsw
484
485 ! Drop evaporation
486 call table_dropevap ! => fills tpc_wev and tnc_wev
487
488 ! Read a static file containing CCN activation of aerosols. The data were created from a parcel model
489 ! by Feingold & Heymsfield with further changes by Eidhammer and Kriedenweis.
490 ! call mpas_new_unit(mp_unit, unformatted = .true.)
491 mp_unit = 11
492 ! open(unit=mp_unit,file='CCN_ACTIVATE.BIN',form='unformatted',status='old',action='read',iostat=istat)
493 ! if (istat /= open_OK) then
494 ! call physics_error_fatal('--- tempo_init() failure opening CCN_ACTIVATE.BIN')
495 ! endif
496 ! read(mp_unit) tnccn_act
497 ! close(unit=mp_unit)
498
499 ! Rain collecting graupel & graupel collecting rain.
500 if (configs%hail_aware) then
501 using_hail_aware_table = .true.
502 open(unit=mp_unit,file='MP_TEMPO_HAILAWARE_QRacrQG_DATA.DBL',form='unformatted',status='old',action='read', &
503 iostat=istat, convert='big_endian')
504! if (istat /= open_OK) then
505! call physics_error_fatal('--- tempo_init() failure opening MP_TEMPO_HAILAWARE_QRacrQG.DBL')
506! endif
507 read(mp_unit) tcg_racg
508 read(mp_unit) tmr_racg
509 read(mp_unit) tcr_gacr
510 read(mp_unit) tnr_racg
511 read(mp_unit) tnr_gacr
512 close(unit=mp_unit)
513 else
514 inquire(file='MP_TEMPO_HAILAWARE_QRacrQG_DATA.DBL', exist=qr_acr_qg_hailaware_exists)
515 inquire(file='MP_TEMPO_QRacrQG_DATA.DBL', exist=qr_acr_qg_exists)
516
517 if (qr_acr_qg_hailaware_exists) then
518 using_hail_aware_table = .true.
519 open(unit=mp_unit,file='MP_TEMPO_HAILAWARE_QRacrQG_DATA.DBL',form='unformatted',status='old', &
520 action='read',iostat=istat, convert='big_endian')
521! if (istat /= open_OK) then
522! call physics_error_fatal('--- tempo_init() failure opening MP_TEMPO_HAILAWARE_QRacrQG.DBL')
523! endif
524 elseif (qr_acr_qg_exists) then
525 using_hail_aware_table = .false.
526 open(unit=mp_unit,file='MP_TEMPO_QRacrQG_DATA.DBL',form='unformatted',status='old', &
527 action='read',iostat=istat, convert='big_endian')
528! if (istat /= open_OK) then
529! call physics_error_fatal('--- tempo_init() failure opening MP_TEMPO_QRacrQG.DBL')
530! endif
531 else
532! call physics_error_fatal('--- tempo_init() could not find file to read QRacrQG data.')
533 endif
534 read(mp_unit) tcg_racg
535 read(mp_unit) tmr_racg
536 read(mp_unit) tcr_gacr
537 read(mp_unit) tnr_racg
538 read(mp_unit) tnr_gacr
539 close(unit=mp_unit)
540 endif
541
542 ! Rain collecting snow & snow collecting rain.
543 open(unit=mp_unit,file='MP_TEMPO_QRacrQS_DATA.DBL',form='unformatted',status='old',action='read', &
544 iostat=istat, convert='big_endian')
545! if (istat /= open_OK) then
546! call physics_error_fatal('--- tempo_init() failure opening MP_TEMPO_QRacrQS.DBL')
547! endif
548 read(mp_unit) tcs_racs1
549 read(mp_unit) tmr_racs1
550 read(mp_unit) tcs_racs2
551 read(mp_unit) tmr_racs2
552 read(mp_unit) tcr_sacr1
553 read(mp_unit) tms_sacr1
554 read(mp_unit) tcr_sacr2
555 read(mp_unit) tms_sacr2
556 read(mp_unit) tnr_racs1
557 read(mp_unit) tnr_racs2
558 read(mp_unit) tnr_sacr1
559 read(mp_unit) tnr_sacr2
560 close(unit=mp_unit)
561
562 ! Cloud water and rain freezing (Bigg, 1953).
563 open(unit=mp_unit,file='MP_TEMPO_freezeH2O_DATA.DBL',form='unformatted',status='old',action='read', &
564 iostat=istat, convert='big_endian')
565! if (istat /= open_OK) then
566! call physics_error_fatal('--- tempo_init() failure opening MP_TEMPO_freezeH2O.DBL')
567! endif
568 read(mp_unit) tpi_qrfz
569 read(mp_unit) tni_qrfz
570 read(mp_unit) tpg_qrfz
571 read(mp_unit) tnr_qrfz
572 read(mp_unit) tpi_qcfz
573 read(mp_unit) tni_qcfz
574 close(unit=mp_unit)
575
576 ! Conversion of some ice mass into snow category.
577 open(unit=mp_unit,file='MP_TEMPO_QIautQS_DATA.DBL',form='unformatted',status='old',action='read', &
578 iostat=istat, convert='big_endian')
579! if (istat /= open_OK) then
580! call physics_error_fatal('--- tempo_init() failure opening MP_TEMPO_QIautQS.DBL')
581! endif
582 read(mp_unit) tpi_ide
583 read(mp_unit) tps_iaus
584 read(mp_unit) tni_iaus
585 close(unit=mp_unit)
586 ! call mpas_release_unit(mp_unit)
587
588 ! Initialize various constants for computing radar reflectivity.
589! xam_r = am_r
590! xbm_r = bm_r
591! xmu_r = mu_r
592! xam_s = am_s
593! xbm_s = bm_s
594! xmu_s = mu_s
595! xam_g = am_g(idx_bg1)
596! xbm_g = bm_g
597! xmu_g = mu_g
598! call radar_init
599
600 endif ! micro_init
601
602 end subroutine tempo_init
603
604 !=================================================================================================================
605 ! This is a wrapper routine designed to transfer values from 3D to 1D.
606 ! Required microphysics variables are qv, qc, qr, nr, qi, ni, qs, qg
607 ! Optional microphysics variables are aerosol aware (nc, nwfa, nifa, nwfa2d, nifa2d), and hail aware (ng, qg)
608
609 subroutine tempo_3d_to_1d_driver(qv, qc, qr, qi, qs, qg, qb, ni, nr, nc, ng, &
610 nwfa, nifa, nwfa2d, nifa2d, th, pii, p, w, dz, dt_in, itimestep, &
611 rainnc, rainncv, snownc, snowncv, graupelnc, graupelncv, sr, frainnc, &
612 refl_10cm, diagflag, do_radar_ref, re_cloud, re_ice, re_snow, &
613 has_reqc, has_reqi, has_reqs, ntc, muc, &
614 ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte)
615
616 ! Subroutine (3D) arguments
617 integer, intent(in) :: ids,ide, jds,jde, kds,kde, ims,ime, jms,jme, kms,kme, its,ite, jts,jte, kts,kte
618 real, dimension(ims:ime, kms:kme, jms:jme), intent(inout) :: qv, qc, qr, qi, qs, qg, ni, nr, th
619 real, dimension(ims:ime, kms:kme, jms:jme), intent(inout) :: re_cloud, re_ice, re_snow
620 integer, intent(in) :: has_reqc, has_reqi, has_reqs
621 real, dimension(ims:ime, kms:kme, jms:jme), intent(in) :: pii, p, w, dz
622 real, dimension(ims:ime, jms:jme), intent(inout) :: rainnc, rainncv, sr
623 real, optional, dimension(:,:), intent(inout) :: frainnc
624 real, dimension(ims:ime, jms:jme), intent(in), optional :: ntc, muc
625 real, dimension(ims:ime, kms:kme, jms:jme), intent(inout), optional :: nc, nwfa, nifa, qb, ng
626 real, dimension(ims:ime, jms:jme), intent(in), optional :: nwfa2d, nifa2d
627 real, dimension(ims:ime, kms:kme, jms:jme), intent(inout), optional :: refl_10cm
628 real, dimension(ims:ime, jms:jme), intent(inout), optional :: snownc, snowncv, graupelnc, graupelncv
629 real, intent(in) :: dt_in
630 integer, intent(in) :: itimestep
631
632 ! Local (1d) variables
633 real, dimension(kts:kte) :: qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, ni1d, nr1d, nc1d, ng1d, &
634 nwfa1d, nifa1d, t1d, p1d, w1d, dz1d, rho, dbz
635 real, dimension(kts:kte) :: re_qc1d, re_qi1d, re_qs1d
636 real, dimension(its:ite, jts:jte) :: pcp_ra, pcp_sn, pcp_gr, pcp_ic, frain
637 real :: dt, pptrain, pptsnow, pptgraul, pptice
638 real :: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max
639 real :: nwfa1
640 real :: ygra1, zans1
641 real :: graupel_vol
642 double precision :: lamg, lam_exp, lamr, n0_min, n0_exp
643 integer :: i, j, k
644 integer :: imax_qc, imax_qr, imax_qi, imax_qs, imax_qg, imax_ni, imax_nr
645 integer :: jmax_qc, jmax_qr, jmax_qi, jmax_qs, jmax_qg, jmax_ni, jmax_nr
646 integer :: kmax_qc, kmax_qr, kmax_qi, kmax_qs, kmax_qg, kmax_ni, kmax_nr
647 integer :: i_start, j_start, i_end, j_end
648 logical, optional, intent(in) :: diagflag
649 integer, optional, intent(in) :: do_radar_ref
650 character(len=132) :: message
651
652 !=================================================================================================================
653 i_start = its
654 j_start = jts
655 i_end = min(ite, ide)
656 j_end = min(jte, jde)
657 dt = dt_in
658
659 qc_max = 0.0
660 qr_max = 0.0
661 qs_max = 0.0
662 qi_max = 0.0
663 qg_max = 0.0
664 ni_max = 0.0
665 nr_max = 0.0
666 imax_qc = 0
667 imax_qr = 0
668 imax_qi = 0
669 imax_qs = 0
670 imax_qg = 0
671 imax_ni = 0
672 imax_nr = 0
673 jmax_qc = 0
674 jmax_qr = 0
675 jmax_qi = 0
676 jmax_qs = 0
677 jmax_qg = 0
678 jmax_ni = 0
679 jmax_nr = 0
680 kmax_qc = 0
681 kmax_qr = 0
682 kmax_qi = 0
683 kmax_qs = 0
684 kmax_qg = 0
685 kmax_ni = 0
686 kmax_nr = 0
687
688 !=================================================================================================================
689 j_loop: do j = j_start, j_end
690 i_loop: do i = i_start, i_end
691 pptrain = 0.0
692 pptsnow = 0.0
693 pptgraul = 0.0
694 pptice = 0.0
695 rainncv(i,j) = 0.0
696 if (present(snowncv)) then
697 snowncv(i,j) = 0.0
698 endif
699 if (present(graupelncv)) then
700 graupelncv(i,j) = 0.0
701 endif
702 sr(i,j) = 0.0
703
704 ! ntc and muc are defined in mpas submodule based on landmask
705 if (present(ntc)) then
706 nt_c = ntc(i,j)
707 mu_c = muc(i,j)
708 else
709 nt_c = nt_c_o
710 mu_c = 4
711 endif
712
713 !=================================================================================================================
714 ! Begin k loop
715 do k = kts, kte
716 t1d(k) = th(i,k,j) * pii(i,k,j)
717 p1d(k) = p(i,k,j)
718 w1d(k) = w(i,k,j)
719 dz1d(k) = dz(i,k,j)
720 qv1d(k) = qv(i,k,j)
721 qc1d(k) = qc(i,k,j)
722 qi1d(k) = qi(i,k,j)
723 qr1d(k) = qr(i,k,j)
724 qs1d(k) = qs(i,k,j)
725 qg1d(k) = qg(i,k,j)
726 ni1d(k) = ni(i,k,j)
727 nr1d(k) = nr(i,k,j)
728 rho(k) = roverrv * p1d(k) / (r * t1d(k) * (qv1d(k)+roverrv))
729
730 ! nwfa, nifa, and nc are optional aerosol-aware variables
731 if (present(nwfa)) then
732 if (present(nwfa2d)) then
733 if (k == kts) then
734 nwfa(i,k,j) = nwfa(i,k,j) + nwfa2d(i,j) * dt
735 endif
736 endif
737 nwfa(i,k,j) = max(nwfa_default, min(aero_max, nwfa(i,k,j)))
738 nwfa1d(k) = nwfa(i,k,j)
739 else
740 nwfa1d(k) = nwfa_default / rho(k)
741 configs%aerosol_aware = .false.
742 endif
743
744 if (present(nifa)) then
745 nifa1d(k) = nifa(i,k,j)
746 else
747 nifa1d(k) = nifa_default / rho(k)
748 configs%aerosol_aware = .false.
749 endif
750
751 if (present(nc)) then
752 nc1d(k) = nc(i,k,j)
753 else
754 nc1d(k) = nt_c / rho(k)
755 configs%aerosol_aware = .false.
756 endif
757 enddo
758
759 ! ng and qb are optional hail-aware variables
760 if ((present(ng)) .and. (present(qb))) then
761 configs%hail_aware = .true.
762 do k = kts, kte
763 ng1d(k) = ng(i,k,j)
764 qb1d(k) = qb(i,k,j)
765 enddo
766 else
767 do k = kte, kts, -1
768 ! This is the one-moment graupel formulation
769 if (qg1d(k) > r1) then
770 ygra1 = log10(max(1.e-9, qg1d(k)*rho(k)))
771 zans1 = 3.0 + 2.0/7.0*(ygra1+8.0)
772 zans1 = max(2.0, min(zans1, 6.0))
773 n0_exp = 10.0**(zans1)
774 lam_exp = (n0_exp*am_g(idx_bg1)*cgg(1,1) / (rho(k)*qg1d(k)))**oge1
775 lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg
776 ng1d(k) = cgg(2,1) * ogg3*rho(k) * qg1d(k) * lamg**bm_g / am_g(idx_bg1)
777 ng1d(k) = max(r2, (ng1d(k)/rho(k)))
778 qb1d(k) = qg1d(k) / rho_g(idx_bg1)
779 else
780 ng1d(k) = 0
781 qb1d(k) = 0
782 endif
783 enddo
784 endif
785
786 !if (itimestep == 1) then
787 ! call physics_message('--- tempo_3d_to_1d_driver() configuration...')
788 ! write(message, '(L1)') configs%hail_aware
789 ! call physics_message(' hail_aware_flag = ' // trim(message))
790 ! write(message, '(L1)') configs%aerosol_aware
791 ! call physics_message(' aerosol_aware_flag = ' // trim(message))
792 ! call physics_message('calling mp_tempo_main() at itimestep = 1')
793 !endif
794 !=================================================================================================================
795 ! Main call to the 1D microphysics
796 call mp_tempo_main(qv1d=qv1d, qc1d=qc1d, qi1d=qi1d, qr1d=qr1d, qs1d=qs1d, qg1d=qg1d, qb1d=qb1d, &
797 ni1d=ni1d, nr1d=nr1d, nc1d=nc1d, ng1d=ng1d, nwfa1d=nwfa1d, nifa1d=nifa1d, t1d=t1d, p1d=p1d, &
798 w1d=w1d, dzq=dz1d, pptrain=pptrain, pptsnow=pptsnow, pptgraul=pptgraul, pptice=pptice, &
799 kts=kts, kte=kte, dt=dt, ii=i, jj=j, configs=configs)
800
801 !=================================================================================================================
802 ! Compute diagnostics and return output to 3D
803 pcp_ra(i,j) = pptrain
804 pcp_sn(i,j) = pptsnow
805 pcp_gr(i,j) = pptgraul
806 pcp_ic(i,j) = pptice
807 rainncv(i,j) = pptrain + pptsnow + pptgraul + pptice
808 rainnc(i,j) = rainnc(i,j) + pptrain + pptsnow + pptgraul + pptice
809 if (present(snowncv) .and. present(snownc)) then
810 snowncv(i,j) = pptsnow + pptice
811 snownc(i,j) = snownc(i,j) + pptsnow + pptice
812 endif
813 if (present(graupelncv) .and. present(graupelnc)) then
814 graupelncv(i,j) = pptgraul
815 graupelnc(i,j) = graupelnc(i,j) + pptgraul
816 endif
817 if (present(frainnc)) then
818 frain(i,j) = 0.
819 if(t1d(1) <= 273.) then
820 frain(i,j) = pcp_ra(i,j)
821 endif
822 frainnc(i,j) = frainnc(i,j) + frain(i,j)
823 endif
824
825 sr(i,j) = (pptsnow + pptgraul + pptice) / (rainncv(i,j) + r1)
826
827 ! ng and qb are optional hail-aware variables
828 if ((present(ng)) .and. (present(qb))) then
829 do k = kts, kte
830 ng(i,k,j) = ng1d(k)
831 qb(i,k,j) = qb1d(k)
832 enddo
833 else
834 do k = kte, kts, -1
835 ! This is the one-moment graupel formulation
836 if (qg1d(k) > r1) then
837 ygra1 = log10(max(1.e-9, qg1d(k)*rho(k)))
838 zans1 = 3.0 + 2.0/7.0*(ygra1+8.0)
839 zans1 = max(2.0, min(zans1, 6.0))
840 n0_exp = 10.0**(zans1)
841 lam_exp = (n0_exp*am_g(idx_bg1)*cgg(1,1) / (rho(k)*qg1d(k)))**oge1
842 lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg
843 ng1d(k) = cgg(2,1) * ogg3*rho(k) * qg1d(k) * lamg**bm_g / am_g(idx_bg1)
844 ng1d(k) = max(r2, (ng1d(k)/rho(k)))
845 qb1d(k) = qg1d(k) / rho_g(idx_bg1)
846 else
847 ng1d(k) = 0
848 qb1d(k) = 0
849 endif
850 enddo
851 endif
852
853 do k = kts, kte
854 if (present(nc)) nc(i,k,j) = nc1d(k)
855 if (present(nwfa)) nwfa(i,k,j) = nwfa1d(k)
856 if (present(nifa)) nifa(i,k,j) = nifa1d(k)
857 qv(i,k,j) = qv1d(k)
858 qc(i,k,j) = qc1d(k)
859 qi(i,k,j) = qi1d(k)
860 qr(i,k,j) = qr1d(k)
861 qs(i,k,j) = qs1d(k)
862 qg(i,k,j) = qg1d(k)
863 ni(i,k,j) = ni1d(k)
864 nr(i,k,j) = nr1d(k)
865 th(i,k,j) = t1d(k) / pii(i,k,j)
866 enddo
867
868 !=================================================================================================================
869 ! Reflectivity
870 call calc_refl10cm (qv1d=qv1d, qc1d=qc1d, qr1d=qr1d, nr1d=nr1d, qs1d=qs1d, qg1d=qg1d, ng1d=ng1d, qb1d=qb1d, &
871 t1d=t1d, p1d=p1d, dbz=dbz, kts=kts, kte=kte, ii=i, jj=j, configs=configs)
872 do k = kts, kte
873 refl_10cm(i,k,j) = max(-35.0_wp, dbz(k))
874 enddo
875
876 ! Cloud, ice, and snow effective radius
877 if (has_reqc /= 0 .and. has_reqi /= 0 .and. has_reqs /= 0) then
878 do k = kts, kte
879 re_qc1d(k) = 2.49e-6
880 re_qi1d(k) = 4.99e-6
881 re_qs1d(k) = 9.99e-6
882 enddo
883 call calc_effectrad (t1d=t1d, p1d=p1d, qv1d=qv1d, qc1d=qc1d, nc1d=nc1d, qi1d=qi1d, &
884 ni1d=ni1d, qs1d=qs1d, re_qc1d=re_qc1d, re_qi1d=re_qi1d, re_qs1d=re_qs1d, &
885 kts=kts, kte=kte, configs=configs)
886 do k = kts, kte
887 re_cloud(i,k,j) = max(2.49e-6, min(re_qc1d(k), 50.e-6))
888 re_ice(i,k,j) = max(4.99e-6, min(re_qi1d(k), 125.e-6))
889 re_snow(i,k,j) = max(9.99e-6, min(re_qs1d(k), 999.e-6))
890 enddo
891 endif
892
893 enddo i_loop
894 enddo j_loop
895
896 end subroutine tempo_3d_to_1d_driver
897 !=================================================================================================================
898
899end module module_mp_tempo
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
Definition module_mp_tempo.F90:3