CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
module_gfdl_cloud_microphys.F90
1
7!***********************************************************************
8!* GNU Lesser General Public License
9!*
10!* This file is part of the GFDL Cloud Microphysics.
11!*
12!* The GFDL Cloud Microphysics is free software: you can
13!* redistribute it and/or modify it under the terms of the
14!* GNU Lesser General Public License as published by the
15!* Free Software Foundation, either version 3 of the License, or
16!* (at your option) any later version.
17!*
18!* The GFDL Cloud Microphysics is distributed in the hope it will be
19!* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
20!* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
21!* See the GNU General Public License for more details.
22!*
23!* You should have received a copy of the GNU Lesser General Public
24!* License along with the GFDL Cloud Microphysics.
25!* If not, see <http://www.gnu.org/licenses/>.
26!***********************************************************************
27! =======================================================================
32
33 ! use mpp_mod, only: stdlog, mpp_pe, mpp_root_pe, mpp_clock_id, &
34 ! mpp_clock_begin, mpp_clock_end, clock_routine, &
35 ! input_nml_file
36 ! use diag_manager_mod, only: register_diag_field, send_data
37 ! use time_manager_mod, only: time_type, get_time
38 ! use constants_mod, only: grav, rdgas, rvgas, cp_air, hlv, hlf, pi => pi_8
39 ! use fms_mod, only: write_version_number, open_namelist_file, &
40 ! check_nml_error, file_exist, close_file
41
43
44 implicit none
45
46 private
47
50! public wqs1, wqs2, qs_blend, wqsat_moist, wqsat2_moist
51! public qsmith_init, qsmith, es2_table1d, es3_table1d, esw_table1d
52! public setup_con, wet_bulb
53
54 real :: missing_value = - 1.e10
55
56 logical :: module_is_initialized = .false.
57 logical :: qsmith_tables_initialized = .false.
58
59 character (len = 17) :: mod_name = 'gfdl_cloud_microphys'
60
61 real, parameter :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6
62 real, parameter :: rhos = 0.1e3, rhog = 0.4e3
63 real, parameter :: grav = 9.80665
64 real, parameter :: rdgas = 287.05
65 real, parameter :: rvgas = 461.50
66 real, parameter :: cp_air = 1004.6
67 real, parameter :: hlv = 2.5e6
68 real, parameter :: hlf = 3.3358e5
69 real, parameter :: pi = 3.1415926535897931
70
71 ! real, parameter :: rdgas = 287.04 !< gfdl: gas constant for dry air
72
73 ! real, parameter :: cp_air = rdgas * 7. / 2. ! 1004.675, heat capacity of dry air at constant pressure
74 real, parameter :: cp_vap = 4.0 * rvgas
75 ! real, parameter :: cv_air = 717.56 ! satoh value
76 real, parameter :: cv_air = cp_air - rdgas
77 ! real, parameter :: cv_vap = 1410.0 ! emanuel value
78 real, parameter :: cv_vap = 3.0 * rvgas
79
80 ! the following two are from emanuel's book "atmospheric convection"
81 ! real, parameter :: c_ice = 2106.0 ! heat capacity of ice at 0 deg c: c = c_ice + 7.3 * (t - tice)
82 ! real, parameter :: c_liq = 4190.0 ! heat capacity of water at 0 deg c
83
84 real, parameter :: c_ice = 1972.0
85 real, parameter :: c_liq = 4185.5
86 ! real, parameter :: c_liq = 4218.0 ! ifs: heat capacity of liquid at 0 deg c
87
88 real, parameter :: eps = rdgas / rvgas
89 real, parameter :: zvir = rvgas / rdgas - 1.
90
91 real, parameter :: t_ice = 273.16
92 real, parameter :: table_ice = 273.16
93
94 ! real, parameter :: e00 = 610.71 ! gfdl: saturation vapor pressure at 0 deg c
95 real, parameter :: e00 = 611.21
96
97 real, parameter :: dc_vap = cp_vap - c_liq
98 real, parameter :: dc_ice = c_liq - c_ice
99
100 real, parameter :: hlv0 = hlv
101 ! real, parameter :: hlv0 = 2.501e6 ! emanuel appendix - 2
102 real, parameter :: hlf0 = hlf
103 ! real, parameter :: hlf0 = 3.337e5 ! emanuel
104
105 real, parameter :: lv0 = hlv0 - dc_vap * t_ice
106 real, parameter :: li00 = hlf0 - dc_ice * t_ice
107
108 real, parameter :: d2ice = dc_vap + dc_ice
109 real, parameter :: li2 = lv0 + li00
110
111 real, parameter :: qrmin = 1.e-8
112 real, parameter :: qvmin = 1.e-20
113 real, parameter :: qcmin = 1.e-12
114
115 real, parameter :: vr_min = 1.e-3
116 real, parameter :: vf_min = 1.e-5
117
118 real, parameter :: dz_min = 1.e-2
119
120 real, parameter :: sfcrho = 1.2
121 real, parameter :: rhor = 1.e3
122 ! intercept parameters
123
124 real, parameter :: rnzr = 8.0e6 ! lin83
125 real, parameter :: rnzs = 3.0e6 ! lin83
126 real, parameter :: rnzg = 4.0e6 ! rh84
127 real, parameter :: rnzh = 4.0e4 ! lin83 --- lmh 29 sep 17
128
129 ! density parameters
130
131 real, parameter :: rhoh = 0.917e3 ! lin83 --- lmh 29 sep 17
132
133 public rhor, rhos, rhog, rhoh, rnzr, rnzs, rnzg, rnzh
134 real :: cracs, csacr, cgacr, cgacs, csacw, craci, csaci, cgacw, cgaci, cracw
135 real :: acco (3, 4)
136 real :: cssub (5), cgsub (5), crevp (5), cgfr (2), csmlt (5), cgmlt (5)
137
138 real :: es0, ces0
139 real :: pie, rgrav, fac_rc
140 real :: c_air, c_vap
141
142 real :: lati, latv, lats, lat2, lcp, icp, tcp
143
144 real :: d0_vap
145 real :: lv00
146
147 ! cloud microphysics switchers
148
149 integer :: icloud_f = 0
150 integer :: irain_f = 0
151
152 logical :: de_ice = .false.
153 logical :: sedi_transport = .true.
154 logical :: do_sedi_w = .false.
155 logical :: do_sedi_heat = .true.
156 logical :: prog_ccn = .false.
157 logical :: do_qa = .true.
158 logical :: rad_snow = .true.
159 logical :: rad_graupel = .true.
160 logical :: rad_rain = .true.
161 logical :: fix_negative = .false.
162 logical :: do_setup = .true.
163 logical :: p_nonhydro = .false.
164
165 real, allocatable :: table (:), table2 (:), table3 (:), tablew (:)
166 real, allocatable :: des (:), des2 (:), des3 (:), desw (:)
167
168 logical :: tables_are_initialized = .false.
169
170 ! logical :: master
171 ! integer :: id_rh, id_vtr, id_vts, id_vtg, id_vti, id_rain, id_snow, id_graupel, &
172 ! id_ice, id_prec, id_cond, id_var, id_droplets
173 real, parameter :: dt_fr = 8.
174 ! minimum temperature water can exist (moore & molinero nov. 2011, nature)
175 ! dt_fr can be considered as the error bar
176
177 real :: p_min = 100.
178
179 ! slj, the following parameters are for cloud - resolving resolution: 1 - 5 km
180
181 ! qi0_crt = 0.8e-4
182 ! qs0_crt = 0.6e-3
183 ! c_psaci = 0.1
184 ! c_pgacs = 0.1
185
186 ! -----------------------------------------------------------------------
187 ! namelist parameters
188 ! -----------------------------------------------------------------------
189
190 real :: cld_min = 0.05
191 real :: tice = 273.16
192
193 real :: t_min = 178.
194 real :: t_sub = 184.
195 real :: mp_time = 150.
196
197 ! relative humidity increment
198
199 real :: rh_inc = 0.25
200 real :: rh_inr = 0.25
201 real :: rh_ins = 0.25
202
203 ! conversion time scale
204
205 real :: tau_r2g = 900.
206 real :: tau_smlt = 900.
207 real :: tau_g2r = 600.
208 real :: tau_imlt = 600.
209 real :: tau_i2s = 1000.
210 real :: tau_l2r = 900.
211 real :: tau_v2l = 150.
212 real :: tau_l2v = 300.
213 real :: tau_g2v = 900.
214 real :: tau_v2g = 21600.
215
216 ! horizontal subgrid variability
217
218 real :: dw_land = 0.20
219 real :: dw_ocean = 0.10
220
221 ! prescribed ccn
222
223 real :: ccn_o = 90.
224 real :: ccn_l = 270.
225
226 real :: rthresh = 10.0e-6
227
228 ! -----------------------------------------------------------------------
229 ! wrf / wsm6 scheme: qi_gen = 4.92e-11 * (1.e3 * exp (0.1 * tmp)) ** 1.33
230 ! optimized: qi_gen = 4.92e-11 * exp (1.33 * log (1.e3 * exp (0.1 * tmp)))
231 ! qi_gen ~ 4.808e-7 at 0 c; 1.818e-6 at - 10 c, 9.82679e-5 at - 40c
232 ! the following value is constructed such that qc_crt = 0 at zero c and @ - 10c matches
233 ! wrf / wsm6 ice initiation scheme; qi_crt = qi_gen * min (qi_lim, 0.1 * tmp) / den
234 ! -----------------------------------------------------------------------
235
236 real :: sat_adj0 = 0.90
237
238 real :: qc_crt = 5.0e-8
239
240 real :: qi_lim = 1.
241
242 real :: ql_mlt = 2.0e-3
243 real :: qs_mlt = 1.0e-6
244
245 real :: ql_gen = 1.0e-3
246 real :: qi_gen = 1.82e-6
247
248 ! cloud condensate upper bounds: "safety valves" for ql & qi
249
250 real :: ql0_max = 2.0e-3
251 real :: qi0_max = 1.0e-4
252
253 real :: qi0_crt = 1.0e-4
255 real :: qr0_crt = 1.0e-4
256 ! lfo used * mixing ratio * = 1.e-4 (hail in lfo)
257 real :: qs0_crt = 1.0e-3
258
259 real :: c_paut = 0.55
260 real :: c_psaci = 0.02
261 real :: c_piacr = 5.0
262 real :: c_cracw = 0.9
263 real :: c_pgacs = 2.0e-3
264
265 ! decreasing clin to reduce csacw (so as to reduce cloud water --- > snow)
266
267 real :: alin = 842.0
268 real :: clin = 4.8
269
270 ! fall velocity tuning constants:
271
272 logical :: const_vi = .false.
273 logical :: const_vs = .false.
274 logical :: const_vg = .false.
275 logical :: const_vr = .false.
276
277 ! good values:
278
279 real :: vi_fac = 1.
280 real :: vs_fac = 1.
281 real :: vg_fac = 1.
282 real :: vr_fac = 1.
283
284 ! upper bounds of fall speed (with variable speed option)
285
286 real :: vi_max = 0.5
287 real :: vs_max = 5.0
288 real :: vg_max = 8.0
289 real :: vr_max = 12.
290
291 ! cloud microphysics switchers
292
293 logical :: fast_sat_adj = .false.
294 logical :: z_slope_liq = .true.
295 logical :: z_slope_ice = .false.
296 logical :: use_ccn = .false.
297 logical :: use_ppm = .false.
298 logical :: mono_prof = .true.
299 logical :: mp_print = .false.
300 logical :: do_hail = .false.
301
302 ! real :: global_area = - 1.
303
304 real :: log_10, tice0, t_wfr
305
306 integer :: reiflag = 1
307 ! 1: Heymsfield and Mcfarquhar, 1996
308 ! 2: Wyser, 1998
309
310 logical :: tintqs = .false.
311
312 real :: rewmin = 5.0, rewmax = 10.0
313 real :: reimin = 10.0, reimax = 150.0
314 real :: rermin = 10.0, rermax = 10000.0
315 real :: resmin = 150.0, resmax = 10000.0
316 real :: regmin = 300.0, regmax = 10000.0
317
318 ! -----------------------------------------------------------------------
319 ! namelist
320 ! -----------------------------------------------------------------------
321
322 namelist / gfdl_cloud_microphysics_nml / &
323 mp_time, t_min, t_sub, tau_r2g, tau_smlt, tau_g2r, dw_land, dw_ocean, &
324 vi_fac, vr_fac, vs_fac, vg_fac, ql_mlt, do_qa, fix_negative, vi_max, &
325 vs_max, vg_max, vr_max, qs_mlt, qs0_crt, qi_gen, ql0_max, qi0_max, &
326 qi0_crt, qr0_crt, fast_sat_adj, rh_inc, rh_ins, rh_inr, const_vi, &
327 const_vs, const_vg, const_vr, use_ccn, rthresh, ccn_l, ccn_o, qc_crt, &
328 tau_g2v, tau_v2g, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v, &
329 tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs, &
330 z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice, &
331 rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof, &
332 do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, &
333 mp_print, reiflag, rewmin, rewmax, reimin, reimax, rermin, rermax, &
334 resmin, resmax, regmin, regmax, tintqs, do_hail
335
336 public &
337 mp_time, t_min, t_sub, tau_r2g, tau_smlt, tau_g2r, dw_land, dw_ocean, &
338 vi_fac, vr_fac, vs_fac, vg_fac, ql_mlt, do_qa, fix_negative, vi_max, &
339 vs_max, vg_max, vr_max, qs_mlt, qs0_crt, qi_gen, ql0_max, qi0_max, &
340 qi0_crt, qr0_crt, fast_sat_adj, rh_inc, rh_ins, rh_inr, const_vi, &
341 const_vs, const_vg, const_vr, use_ccn, rthresh, ccn_l, ccn_o, qc_crt, &
342 tau_g2v, tau_v2g, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v, &
343 tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs, &
344 z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice, &
345 rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof, &
346 do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, &
347 mp_print, reiflag, rewmin, rewmax, reimin, reimax, rermin, rermax, &
348 resmin, resmax, regmin, regmax, tintqs, do_hail
349
350contains
351
352! -----------------------------------------------------------------------
353! the driver of the gfdl cloud microphysics
354! -----------------------------------------------------------------------
355
359 iis, iie, jjs, jje, kks, kke, ktop, kbot, &
360 qv, ql, qr, qi, qs, qg, qa, qn, &
361 qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, pt_dt, pt, w, &
362 uin, vin, udt, vdt, dz, delp, area, dt_in, land, &
363 rain, snow, ice, graupel, hydrostatic, phys_hydrostatic, &
364 p, lradar, refl_10cm, reset, pfils, pflls)
365
366 implicit none
367
368 ! Interface variables
369 integer, intent (in) :: iis, iie, jjs, jje ! physics window
370 integer, intent (in) :: kks, kke ! vertical dimension
371 integer, intent (in) :: ktop, kbot ! vertical compute domain
372
373 real, intent (in) :: dt_in ! physics time step
374
375 real, intent (in), dimension (iis:iie, jjs:jje) :: area ! cell area
376 real, intent (in), dimension (iis:iie, jjs:jje) :: land ! land fraction
377
378 real, intent (in), dimension (iis:iie, jjs:jje, kks:kke) :: delp, dz, uin, vin
379 real, intent (in), dimension (iis:iie, jjs:jje, kks:kke) :: pt, qv, ql, qr, qg, qa, qn
380
381 real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: qi, qs
382 real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: pt_dt, qa_dt, udt, vdt, w
383 real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: qv_dt, ql_dt, qr_dt
384 real, intent (inout), dimension (iis:iie, jjs:jje, kks:kke) :: qi_dt, qs_dt, qg_dt
385
386 real, intent (out), dimension (iis:iie, jjs:jje) :: rain, snow, ice, graupel
387
388 logical, intent (in) :: hydrostatic, phys_hydrostatic
389
390 !integer, intent (in) :: seconds
391 real, intent (in), dimension (iis:iie, jjs:jje, kks:kke) :: p
392 logical, intent (in) :: lradar
393 real, intent (out), dimension (iis:iie, jjs:jje, kks:kke) :: refl_10cm
394 logical, intent (in) :: reset
395 real, intent (out), dimension (iis:iie, jjs:jje, kks:kke) :: pfils, pflls
396
397 ! Local variables
398 logical :: melti = .false.
399
400 real :: mpdt, rdt, dts, convt, tot_prec
401
402 integer :: i, j, k
403 integer :: is, ie, js, je ! physics window
404 integer :: ks, ke ! vertical dimension
405 integer :: days, ntimes, kflip
406
407 real, dimension (iie-iis+1, jje-jjs+1) :: prec_mp, prec1, cond, w_var, rh0
408
409 real, dimension (iie-iis+1, jje-jjs+1,kke-kks+1) :: vt_r, vt_s, vt_g, vt_i, qn2
410
411 real, dimension (size(pt,1), size(pt,3)) :: m2_rain, m2_sol
412
413 real :: allmax
414!+---+-----------------------------------------------------------------+
415!For 3D reflectivity calculations
416 real, dimension(ktop:kbot):: qv1d, t1d, p1d, qr1d, qs1d, qg1d, dbz
417!+---+-----------------------------------------------------------------+
418
419 is = 1
420 js = 1
421 ks = 1
422 ie = iie - iis + 1
423 je = jje - jjs + 1
424 ke = kke - kks + 1
425 ! call mpp_clock_begin (gfdl_mp_clock)
426
427 ! -----------------------------------------------------------------------
428 ! define heat capacity of dry air and water vapor based on hydrostatical property
429 ! -----------------------------------------------------------------------
430
431 if (phys_hydrostatic .or. hydrostatic) then
432 c_air = cp_air
433 c_vap = cp_vap
434 p_nonhydro = .false.
435 else
436 c_air = cv_air
437 c_vap = cv_vap
438 p_nonhydro = .true.
439 endif
440 d0_vap = c_vap - c_liq
441 lv00 = hlv0 - d0_vap * t_ice
442
443 if (hydrostatic) do_sedi_w = .false.
444
445 ! -----------------------------------------------------------------------
446 ! define latent heat coefficient used in wet bulb and Bigg mechanism
447 ! -----------------------------------------------------------------------
448
449 latv = hlv
450 lati = hlf
451 lats = latv + lati
452 lat2 = lats * lats
453
454 lcp = latv / cp_air
455 icp = lati / cp_air
456 tcp = (latv + lati) / cp_air
457
458 ! tendency zero out for am moist processes should be done outside the driver
459
460 ! -----------------------------------------------------------------------
461 ! define cloud microphysics sub time step
462 ! -----------------------------------------------------------------------
463
464 mpdt = min(dt_in, mp_time)
465 rdt = 1. / dt_in
466 ntimes = nint(dt_in / mpdt)
467
468 ! small time step:
469 dts = dt_in / real(ntimes)
470
471 ! call get_time (time, seconds, days)
472
473 ! -----------------------------------------------------------------------
474 ! initialize precipitation
475 ! -----------------------------------------------------------------------
476
477 do j = js, je
478 do i = is, ie
479 graupel(i, j) = 0.
480 rain(i, j) = 0.
481 snow(i, j) = 0.
482 ice(i, j) = 0.
483 cond(i, j) = 0.
484 enddo
485 enddo
486
487 pfils = 0.
488 pflls = 0.
489
490 ! -----------------------------------------------------------------------
491 ! major cloud microphysics
492 ! -----------------------------------------------------------------------
493
494 do j = js, je
495 call mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, qg,&
496 qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes, &
497 rain(:, j), snow(:, j), graupel(:, j), ice(:, j), m2_rain, &
498 m2_sol, cond(:, j), area(:, j), land(:, j), udt, vdt, pt_dt, &
499 qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, w_var, vt_r, &
500 vt_s, vt_g, vt_i, qn2)
501 do k = ktop, kbot
502 do i = is, ie
503 pfils(i, j, k) = m2_sol(i, k)
504 pflls(i, j, k) = m2_rain(i, k)
505 enddo
506 enddo
507 enddo
508
509 ! -----------------------------------------------------------------------
510 ! no clouds allowed above ktop
511 ! -----------------------------------------------------------------------
512
513 if (ks < ktop) then
514 do k = ks, ktop
515 if (do_qa) then
516 do j = js, je
517 do i = is, ie
518 qa_dt(i, j, k) = 0.
519 enddo
520 enddo
521 else
522 do j = js, je
523 do i = is, ie
524 ! qa_dt (i, j, k) = - qa (i, j, k) * rdt
525 qa_dt(i, j, k) = 0. ! gfs
526 enddo
527 enddo
528 endif
529 enddo
530 endif
531
532 ! -----------------------------------------------------------------------
533 ! diagnostic output
534 ! -----------------------------------------------------------------------
535
536 ! if (id_vtr > 0) then
537 ! used = send_data (id_vtr, vt_r, time, is_in = iis, js_in = jjs)
538 ! endif
539
540 ! if (id_vts > 0) then
541 ! used = send_data (id_vts, vt_s, time, is_in = iis, js_in = jjs)
542 ! endif
543
544 ! if (id_vtg > 0) then
545 ! used = send_data (id_vtg, vt_g, time, is_in = iis, js_in = jjs)
546 ! endif
547
548 ! if (id_vti > 0) then
549 ! used = send_data (id_vti, vt_i, time, is_in = iis, js_in = jjs)
550 ! endif
551
552 ! if (id_droplets > 0) then
553 ! used = send_data (id_droplets, qn2, time, is_in = iis, js_in = jjs)
554 ! endif
555
556 ! if (id_var > 0) then
557 ! used = send_data (id_var, w_var, time, is_in = iis, js_in = jjs)
558 ! endif
559
560 ! convert to mm / day
561
562 convt = 86400. * rdt * rgrav
563 do j = js, je
564 do i = is, ie
565 rain(i, j) = rain(i, j) * convt
566 snow(i, j) = snow(i, j) * convt
567 ice(i, j) = ice(i, j) * convt
568 graupel(i, j) = graupel(i, j) * convt
569 prec_mp(i, j) = rain(i, j) + snow(i, j) + ice(i, j) + graupel(i, j)
570 enddo
571 enddo
572
573 ! if (id_cond > 0) then
574 ! do j = js, je
575 ! do i = is, ie
576 ! cond (i, j) = cond (i, j) * rgrav
577 ! enddo
578 ! enddo
579 ! used = send_data (id_cond, cond, time, is_in = iis, js_in = jjs)
580 ! endif
581
582 ! if (id_snow > 0) then
583 ! used = send_data (id_snow, snow, time, iis, jjs)
584 ! used = send_data (id_snow, snow, time, is_in = iis, js_in = jjs)
585 ! if (mp_print .and. seconds == 0) then
586 ! tot_prec = g_sum (snow, is, ie, js, je, area, 1)
587 ! if (master) write (*, *) 'mean snow = ', tot_prec
588 ! endif
589 ! endif
590 !
591 ! if (id_graupel > 0) then
592 ! used = send_data (id_graupel, graupel, time, iis, jjs)
593 ! used = send_data (id_graupel, graupel, time, is_in = iis, js_in = jjs)
594 ! if (mp_print .and. seconds == 0) then
595 ! tot_prec = g_sum (graupel, is, ie, js, je, area, 1)
596 ! if (master) write (*, *) 'mean graupel = ', tot_prec
597 ! endif
598 ! endif
599 !
600 ! if (id_ice > 0) then
601 ! used = send_data (id_ice, ice, time, iis, jjs)
602 ! used = send_data (id_ice, ice, time, is_in = iis, js_in = jjs)
603 ! if (mp_print .and. seconds == 0) then
604 ! tot_prec = g_sum (ice, is, ie, js, je, area, 1)
605 ! if (master) write (*, *) 'mean ice_mp = ', tot_prec
606 ! endif
607 ! endif
608 !
609 ! if (id_rain > 0) then
610 ! used = send_data (id_rain, rain, time, iis, jjs)
611 ! used = send_data (id_rain, rain, time, is_in = iis, js_in = jjs)
612 ! if (mp_print .and. seconds == 0) then
613 ! tot_prec = g_sum (rain, is, ie, js, je, area, 1)
614 ! if (master) write (*, *) 'mean rain = ', tot_prec
615 ! endif
616 ! endif
617 !
618 ! if (id_rh > 0) then !not used?
619 ! used = send_data (id_rh, rh0, time, iis, jjs)
620 ! used = send_data (id_rh, rh0, time, is_in = iis, js_in = jjs)
621 ! endif
622 !
623 !
624 ! if (id_prec > 0) then
625 ! used = send_data (id_prec, prec_mp, time, iis, jjs)
626 ! used = send_data (id_prec, prec_mp, time, is_in = iis, js_in = jjs)
627 ! endif
628
629 ! if (mp_print) then
630 ! prec1 (:, :) = prec1 (:, :) + prec_mp (:, :)
631 ! if (seconds == 0) then
632 ! prec1 (:, :) = prec1 (:, :) * dt_in / 86400.
633 ! tot_prec = g_sum (prec1, is, ie, js, je, area, 1)
634 ! if (master) write (*, *) 'daily prec_mp = ', tot_prec
635 ! prec1 (:, :) = 0.
636 ! endif
637 ! endif
638
639 ! call mpp_clock_end (gfdl_mp_clock)
640 if(lradar) then
641 ! Only set melti to true at the output times
642 if (reset) then
643 melti=.true.
644 else
645 melti=.false.
646 endif
647 do j = js, je
648 do i = is, ie
649 do k = ktop,kbot
650 kflip = kbot-ktop+1-k+1
651 t1d(k) = pt(i,j,kflip)
652 p1d(k) = p(i,j,kflip)
653 qv1d(k) = qv(i,j,kflip)/(1-qv(i,j,kflip))
654 qr1d(k) = qr(i,j,kflip)
655 qs1d(k) = qs(i,j,kflip)
656 qg1d(k) = qg(i,j,kflip)
657 enddo
658 call refl10cm_gfdl (qv1d, qr1d, qs1d, qg1d, &
659 t1d, p1d, dbz, ktop, kbot, i, j, melti)
660 do k = ktop,kbot
661 kflip = kbot-ktop+1-k+1
662 refl_10cm(i,j,kflip) = max(-35., dbz(k))
663 enddo
664 enddo
665 enddo
666 endif
667
669
670! -----------------------------------------------------------------------
677subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, &
678 qg, qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes, &
679 rain, snow, graupel, ice, m2_rain, m2_sol, cond, area1, land, &
680 u_dt, v_dt, pt_dt, qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, &
681 w_var, vt_r, vt_s, vt_g, vt_i, qn2)
682
683 implicit none
684
685 logical, intent (in) :: hydrostatic
686
687 integer, intent (in) :: j, is, ie, js, je, ks, ke
688 integer, intent (in) :: ntimes, ktop, kbot
689
690 real, intent (in) :: dt_in
691
692 real, intent (in), dimension (is:) :: area1, land
693
694 real, intent (in), dimension (is:, js:, ks:) :: uin, vin, delp, pt, dz
695 real, intent (in), dimension (is:, js:, ks:) :: qv, ql, qr, qg, qa, qn
696
697 real, intent (inout), dimension (is:, js:, ks:) :: qi, qs
698 real, intent (inout), dimension (is:, js:, ks:) :: u_dt, v_dt, w, pt_dt, qa_dt
699 real, intent (inout), dimension (is:, js:, ks:) :: qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt
700
701 real, intent (inout), dimension (is:) :: rain, snow, ice, graupel, cond
702
703 real, intent (out), dimension (is:, js:) :: w_var
704
705 real, intent (out), dimension (is:, js:, ks:) :: vt_r, vt_s, vt_g, vt_i, qn2
706
707 real, intent (out), dimension (is:, ks:) :: m2_rain, m2_sol
708
709 real, dimension (ktop:kbot) :: qvz, qlz, qrz, qiz, qsz, qgz, qaz
710 real, dimension (ktop:kbot) :: vtiz, vtsz, vtgz, vtrz
711 real, dimension (ktop:kbot) :: dp0, dp1, dz0, dz1
712 real, dimension (ktop:kbot) :: qv0, ql0, qr0, qi0, qs0, qg0, qa0
713 real, dimension (ktop:kbot) :: t0, den, den0, tz, p1, denfac
714 real, dimension (ktop:kbot) :: ccn, c_praut, m1_rain, m1_sol, m1
715 real, dimension (ktop:kbot) :: u0, v0, u1, v1, w1
716
717 real :: cpaut, rh_adj, rh_rain
718 real :: r1, s1, i1, g1, rdt, ccn0
719 real :: dt_rain, dts
720 real :: s_leng, t_land, t_ocean, h_var
721 real :: cvm, tmp, omq
722 real :: dqi, qio, qin
723
724 integer :: i, k, n
725
726 dts = dt_in / real(ntimes)
727 dt_rain = dts * 0.5
728 rdt = 1. / dt_in
729
730 ! -----------------------------------------------------------------------
731 ! use local variables
732 ! -----------------------------------------------------------------------
733
734 !GJF: assign values to intent(out) variables that are commented out
735 w_var = 0.0
736 vt_r = 0.0
737 vt_s = 0.0
738 vt_g = 0.0
739 vt_i = 0.0
740 qn2 = 0.0
741 !GJF
742
743 do i = is, ie
744
745 do k = ktop, kbot
746 qiz(k) = qi(i, j, k)
747 qsz(k) = qs(i, j, k)
748 enddo
749
750 ! -----------------------------------------------------------------------
752 ! -----------------------------------------------------------------------
753
754 if (de_ice) then
755 do k = ktop, kbot
756 qio = qiz(k) - dt_in * qi_dt(i, j, k) ! original qi before phys
757 qin = max(qio, qi0_max) ! adjusted value
758 if (qiz(k) > qin) then
759 qsz(k) = qsz(k) + qiz(k) - qin
760 qiz(k) = qin
761 dqi = (qin - qio) * rdt ! modified qi tendency
762 qs_dt(i, j, k) = qs_dt(i, j, k) + qi_dt(i, j, k) - dqi
763 qi_dt(i, j, k) = dqi
764 qi(i, j, k) = qiz(k)
765 qs(i, j, k) = qsz(k)
766 endif
767 enddo
768 endif
769
770 do k = ktop, kbot
771
772 t0(k) = pt(i, j, k)
773 tz(k) = t0(k)
774 dp1(k) = delp(i, j, k)
775 dp0(k) = dp1(k) ! moist air mass * grav
776
777 ! -----------------------------------------------------------------------
779 ! -----------------------------------------------------------------------
780
781 qvz(k) = qv(i, j, k)
782 qlz(k) = ql(i, j, k)
783 qrz(k) = qr(i, j, k)
784 qgz(k) = qg(i, j, k)
785
786 ! dp1: dry air_mass
787 ! dp1 (k) = dp1 (k) * (1. - (qvz (k) + qlz (k) + qrz (k) + qiz (k) + qsz (k) + qgz (k)))
788 dp1(k) = dp1(k) * (1. - qvz(k)) ! gfs
789 omq = dp0(k) / dp1(k)
790
791 qvz(k) = qvz(k) * omq
792 qlz(k) = qlz(k) * omq
793 qrz(k) = qrz(k) * omq
794 qiz(k) = qiz(k) * omq
795 qsz(k) = qsz(k) * omq
796 qgz(k) = qgz(k) * omq
797
798 qa0(k) = qa(i, j, k)
799 qaz(k) = 0.
800 dz0(k) = dz(i, j, k)
801
802 den0(k) = - dp1(k) / (grav * dz0(k)) ! density of dry air
803 p1(k) = den0(k) * rdgas * t0(k) ! dry air pressure
804
805 ! -----------------------------------------------------------------------
806 ! save a copy of old value for computing tendencies
807 ! -----------------------------------------------------------------------
808
809 qv0(k) = qvz(k)
810 ql0(k) = qlz(k)
811 qr0(k) = qrz(k)
812 qi0(k) = qiz(k)
813 qs0(k) = qsz(k)
814 qg0(k) = qgz(k)
815
816 ! -----------------------------------------------------------------------
817 ! for sedi_momentum
818 ! -----------------------------------------------------------------------
819
820 m1(k) = 0.
821 u0(k) = uin(i, j, k)
822 v0(k) = vin(i, j, k)
823 u1(k) = u0(k)
824 v1(k) = v0(k)
825
826 enddo
827
828 if (do_sedi_w) then
829 do k = ktop, kbot
830 w1(k) = w(i, j, k)
831 enddo
832 ! Set to -999 if not used
833 else
834 w1(:) = -999.0
835 endif
836
837 ! -----------------------------------------------------------------------
840 ! -----------------------------------------------------------------------
841
842 cpaut = c_paut * 0.104 * grav / 1.717e-5
843
844 if (prog_ccn) then
845 do k = ktop, kbot
846 ! convert # / cc to # / m^3
847 ccn(k) = qn(i, j, k) * 1.e6
848 c_praut(k) = cpaut * (ccn(k) * rhor) ** (- 1. / 3.)
849 enddo
850 use_ccn = .false.
851 else
852 ccn0 = (ccn_l * land(i) + ccn_o * (1. - land(i))) * 1.e6
853 if (use_ccn) then
854 ! -----------------------------------------------------------------------
855 ! ccn is formulted as ccn = ccn_surface * (den / den_surface)
856 ! -----------------------------------------------------------------------
857 ccn0 = ccn0 * rdgas * tz(kbot) / p1(kbot)
858 endif
859 tmp = cpaut * (ccn0 * rhor) ** (- 1. / 3.)
860 do k = ktop, kbot
861 c_praut(k) = tmp
862 ccn(k) = ccn0
863 enddo
864 endif
865
866 ! -----------------------------------------------------------------------
888 ! total water subgrid deviation in horizontal direction
889 ! default area dependent form: use dx ~ 100 km as the base
890 ! -----------------------------------------------------------------------
891
892 s_leng = sqrt(sqrt(area1(i) / 1.e10))
893 t_land = dw_land * s_leng
894 t_ocean = dw_ocean * s_leng
895 h_var = t_land * land(i) + t_ocean * (1. - land(i))
896 h_var = min(0.20, max(0.01, h_var))
897 ! if (id_var > 0) w_var (i, j) = h_var
898
899 ! -----------------------------------------------------------------------
901 ! -----------------------------------------------------------------------
902
903 rh_adj = 1. - h_var - rh_inc
904 rh_rain = max(0.35, rh_adj - rh_inr) ! rh_inr = 0.25
905
906 ! -----------------------------------------------------------------------
908 ! -----------------------------------------------------------------------
909
910 if (fix_negative) &
911 call neg_adj (ktop, kbot, tz, dp1, qvz, qlz, qrz, qiz, qsz, qgz)
912
913 m2_rain(i, :) = 0.
914 m2_sol(i, :) = 0.
915
917 do n = 1, ntimes
918
919 ! -----------------------------------------------------------------------
921 ! -----------------------------------------------------------------------
922
923 if (p_nonhydro) then
924 do k = ktop, kbot
925 dz1(k) = dz0(k)
926 den(k) = den0(k) ! dry air density remains the same
927 denfac(k) = sqrt(sfcrho / den(k))
928 enddo
929 else
930 do k = ktop, kbot
931 dz1(k) = dz0(k) * tz(k) / t0(k) ! hydrostatic balance
932 den(k) = den0(k) * dz0(k) / dz1(k)
933 denfac(k) = sqrt(sfcrho / den(k))
934 enddo
935 endif
936
937 ! -----------------------------------------------------------------------
939 ! -----------------------------------------------------------------------
940
941 call warm_rain (dt_rain, ktop, kbot, dp1, dz1, tz, qvz, qlz, qrz, qiz, qsz, &
942 qgz, den, denfac, ccn, c_praut, rh_rain, vtrz, r1, m1_rain, w1, h_var)
943
944 rain(i) = rain(i) + r1
945
946 do k = ktop, kbot
947 m2_rain(i, k) = m2_rain(i, k) + m1_rain(k)
948 m1(k) = m1(k) + m1_rain(k)
949 enddo
950
951 ! -----------------------------------------------------------------------
953 ! -----------------------------------------------------------------------
954
957 call fall_speed (ktop, kbot, den, qsz, qiz, qgz, qlz, tz, vtsz, vtiz, vtgz)
958
960 call terminal_fall (dts, ktop, kbot, tz, qvz, qlz, qrz, qgz, qsz, qiz, &
961 dz1, dp1, den, vtgz, vtsz, vtiz, r1, g1, s1, i1, m1_sol, w1)
962
963 rain(i) = rain(i) + r1 ! from melted snow & ice that reached the ground
964 snow(i) = snow(i) + s1
965 graupel(i) = graupel(i) + g1
966 ice(i) = ice(i) + i1
967
968 ! -----------------------------------------------------------------------
970 ! -----------------------------------------------------------------------
971
972 if (do_sedi_heat) &
973 call sedi_heat (ktop, kbot, dp1, m1_sol, dz1, tz, qvz, qlz, qrz, qiz, &
974 qsz, qgz, c_ice)
975
976 ! -----------------------------------------------------------------------
978 ! -----------------------------------------------------------------------
979
980 call warm_rain (dt_rain, ktop, kbot, dp1, dz1, tz, qvz, qlz, qrz, qiz, qsz, &
981 qgz, den, denfac, ccn, c_praut, rh_rain, vtrz, r1, m1_rain, w1, h_var)
982
983 rain(i) = rain(i) + r1
984
985 do k = ktop, kbot
986 m2_rain(i, k) = m2_rain(i, k) + m1_rain(k)
987 m2_sol(i, k) = m2_sol(i, k) + m1_sol(k)
988 m1(k) = m1(k) + m1_rain(k) + m1_sol(k)
989 enddo
990
991 ! -----------------------------------------------------------------------
993 ! -----------------------------------------------------------------------
994
995 call icloud (ktop, kbot, tz, p1, qvz, qlz, qrz, qiz, qsz, qgz, dp1, den, &
996 denfac, vtsz, vtgz, vtrz, qaz, rh_adj, rh_rain, dts, h_var)
997
998 enddo
999
1000 ! convert units from Pa*kg/kg to kg/m^2/s
1001 m2_rain(i, :) = m2_rain(i, :) * rdt * rgrav
1002 m2_sol(i, :) = m2_sol(i, :) * rdt * rgrav
1003
1004 ! -----------------------------------------------------------------------
1006 ! note: dp1 is dry mass; dp0 is the old moist (total) mass
1007 ! -----------------------------------------------------------------------
1008
1009 if (sedi_transport) then
1010 do k = ktop + 1, kbot
1011 u1(k) = (dp0(k) * u1(k) + m1(k - 1) * u1(k - 1)) / (dp0(k) + m1(k - 1))
1012 v1(k) = (dp0(k) * v1(k) + m1(k - 1) * v1(k - 1)) / (dp0(k) + m1(k - 1))
1013 u_dt(i, j, k) = u_dt(i, j, k) + (u1(k) - u0(k)) * rdt
1014 v_dt(i, j, k) = v_dt(i, j, k) + (v1(k) - v0(k)) * rdt
1015 enddo
1016 endif
1017
1018 if (do_sedi_w) then
1019 do k = ktop, kbot
1020 w(i, j, k) = w1(k)
1021 enddo
1022 endif
1023
1024 ! -----------------------------------------------------------------------
1026 ! convert to dry mixing ratios
1027 ! -----------------------------------------------------------------------
1028
1029 do k = ktop, kbot
1030 omq = dp1(k) / dp0(k)
1031 qv_dt(i, j, k) = qv_dt(i, j, k) + rdt * (qvz(k) - qv0(k)) * omq
1032 ql_dt(i, j, k) = ql_dt(i, j, k) + rdt * (qlz(k) - ql0(k)) * omq
1033 qr_dt(i, j, k) = qr_dt(i, j, k) + rdt * (qrz(k) - qr0(k)) * omq
1034 qi_dt(i, j, k) = qi_dt(i, j, k) + rdt * (qiz(k) - qi0(k)) * omq
1035 qs_dt(i, j, k) = qs_dt(i, j, k) + rdt * (qsz(k) - qs0(k)) * omq
1036 qg_dt(i, j, k) = qg_dt(i, j, k) + rdt * (qgz(k) - qg0(k)) * omq
1037 cvm = c_air + qvz(k) * c_vap + (qrz(k) + qlz(k)) * c_liq + (qiz(k) + qsz(k) + qgz(k)) * c_ice
1038 pt_dt(i, j, k) = pt_dt(i, j, k) + rdt * (tz(k) - t0(k)) * cvm / cp_air
1039 enddo
1040
1041 ! -----------------------------------------------------------------------
1043 ! -----------------------------------------------------------------------
1044
1045 do k = ktop, kbot
1046 if (do_qa) then
1047 qa_dt(i, j, k) = 0.
1048 else
1049 qa_dt(i, j, k) = qa_dt(i, j, k) + rdt * (qaz(k) / real(ntimes) - qa0(k))
1050 endif
1051 enddo
1052
1053 ! -----------------------------------------------------------------------
1054 ! fms diagnostics:
1055 ! -----------------------------------------------------------------------
1056
1057 ! if (id_cond > 0) then
1058 ! do k = ktop, kbot ! total condensate
1059 ! cond (i) = cond (i) + dp1 (k) * (qlz (k) + qrz (k) + qsz (k) + qiz (k) + qgz (k))
1060 ! enddo
1061 ! endif
1062 !
1063 ! if (id_vtr > 0) then
1064 ! do k = ktop, kbot
1065 ! vt_r (i, j, k) = vtrz (k)
1066 ! enddo
1067 ! endif
1068 !
1069 ! if (id_vts > 0) then
1070 ! do k = ktop, kbot
1071 ! vt_s (i, j, k) = vtsz (k)
1072 ! enddo
1073 ! endif
1074 !
1075 ! if (id_vtg > 0) then
1076 ! do k = ktop, kbot
1077 ! vt_g (i, j, k) = vtgz (k)
1078 ! enddo
1079 ! endif
1080 !
1081 ! if (id_vts > 0) then
1082 ! do k = ktop, kbot
1083 ! vt_i (i, j, k) = vtiz (k)
1084 ! enddo
1085 ! endif
1086 !
1087 ! if (id_droplets > 0) then
1088 ! do k = ktop, kbot
1089 ! qn2 (i, j, k) = ccn (k)
1090 ! enddo
1091 ! endif
1092
1093 enddo
1094
1095end subroutine mpdrv
1096
1097! -----------------------------------------------------------------------
1100subroutine sedi_heat (ktop, kbot, dm, m1, dz, tz, qv, ql, qr, qi, qs, qg, cw)
1101
1102 implicit none
1103
1104 ! input q fields are dry mixing ratios, and dm is dry air mass
1105
1106 integer, intent (in) :: ktop, kbot
1107
1108 real, intent (in), dimension (ktop:kbot) :: dm, m1, dz, qv, ql, qr, qi, qs, qg
1109
1110 real, intent (inout), dimension (ktop:kbot) :: tz
1111
1112 real, intent (in) :: cw ! heat capacity
1113
1114 real, dimension (ktop:kbot) :: dgz, cvn
1115
1116 real :: tmp
1117
1118 integer :: k
1119
1120 do k = ktop, kbot
1121 dgz(k) = - 0.5 * grav * dz(k) ! > 0
1122 cvn(k) = dm(k) * (cv_air + qv(k) * cv_vap + (qr(k) + ql(k)) * &
1123 c_liq + (qi(k) + qs(k) + qg(k)) * c_ice)
1124 enddo
1125
1126 ! -----------------------------------------------------------------------
1127 ! sjl, july 2014
1128 ! assumption: the ke in the falling condensates is negligible compared to the potential energy
1129 ! that was unaccounted for. local thermal equilibrium is assumed, and the loss in pe is transformed
1130 ! into internal energy (to heat the whole grid box)
1131 ! backward time - implicit upwind transport scheme:
1132 ! dm here is dry air mass
1133 ! -----------------------------------------------------------------------
1134
1135 k = ktop
1136 tmp = cvn(k) + m1(k) * cw
1137 tz(k) = (tmp * tz(k) + m1(k) * dgz(k)) / tmp
1138
1139 ! -----------------------------------------------------------------------
1140 ! implicit algorithm: can't be vectorized
1141 ! needs an inner i - loop for vectorization
1142 ! -----------------------------------------------------------------------
1143
1144 do k = ktop + 1, kbot
1145 tz(k) = ((cvn(k) + cw * (m1(k) - m1(k - 1))) * tz(k) + m1(k - 1) * &
1146 cw * tz(k - 1) + dgz(k) * (m1(k - 1) + m1(k))) / (cvn(k) + cw * m1(k))
1147 enddo
1148
1149end subroutine sedi_heat
1150
1151! -----------------------------------------------------------------------
1155subroutine warm_rain (dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, &
1156 den, denfac, ccn, c_praut, rh_rain, vtr, r1, m1_rain, w1, h_var)
1157
1158 implicit none
1159
1160 integer, intent (in) :: ktop, kbot
1161
1162 real, intent (in) :: dt ! time step (s)
1163 real, intent (in) :: rh_rain, h_var
1164
1165 real, intent (in), dimension (ktop:kbot) :: dp, dz, den
1166 real, intent (in), dimension (ktop:kbot) :: denfac, ccn, c_praut
1167
1168 real, intent (inout), dimension (ktop:kbot) :: tz, vtr
1169 real, intent (inout), dimension (ktop:kbot) :: qv, ql, qr, qi, qs, qg
1170 real, intent (inout), dimension (ktop:kbot) :: m1_rain, w1
1171
1172 real, intent (out) :: r1
1173
1174 real, parameter :: so3 = 7. / 3.
1175
1176 real, dimension (ktop:kbot) :: dl, dm
1177 real, dimension (ktop:kbot + 1) :: ze, zt
1178
1179 real :: sink, dq, qc0, qc
1180 real :: qden
1181 real :: zs = 0.
1182 real :: dt5
1183
1184 integer :: k
1185
1186 ! fall velocity constants:
1187
1188 real, parameter :: vconr = 2503.23638966667
1189 real, parameter :: normr = 25132741228.7183
1190 real, parameter :: thr = 1.e-8
1191
1192 logical :: no_fall
1193
1194 dt5 = 0.5 * dt
1195
1196 ! -----------------------------------------------------------------------
1198 ! -----------------------------------------------------------------------
1199
1200 m1_rain(:) = 0.
1201
1202 call check_column (ktop, kbot, qr, no_fall)
1203
1204 if (no_fall) then
1205 vtr(:) = vf_min
1206 r1 = 0.
1207 else
1208
1209 ! -----------------------------------------------------------------------
1211 ! -----------------------------------------------------------------------
1212
1213 if (const_vr) then
1214 vtr(:) = vr_fac ! ifs_2016: 4.0
1215 else
1216 do k = ktop, kbot
1217 qden = qr(k) * den(k)
1218 if (qr(k) < thr) then
1219 vtr(k) = vr_min
1220 else
1221 vtr(k) = vr_fac * vconr * sqrt(min(10., sfcrho / den(k))) * &
1222 exp(0.2 * log(qden / normr))
1223 vtr(k) = min(vr_max, max(vr_min, vtr(k)))
1224 endif
1225 enddo
1226 endif
1227
1228 ze(kbot + 1) = zs
1229 do k = kbot, ktop, - 1
1230 ze(k) = ze(k + 1) - dz(k) ! dz < 0
1231 enddo
1232
1233 ! -----------------------------------------------------------------------
1236 ! -----------------------------------------------------------------------
1237
1238 ! if (.not. fast_sat_adj) &
1239 call revap_racc (ktop, kbot, dt5, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
1240
1241 if (do_sedi_w) then
1242 do k = ktop, kbot
1243 dm(k) = dp(k) * (1. + qv(k) + ql(k) + qr(k) + qi(k) + qs(k) + qg(k))
1244 enddo
1245 endif
1246
1247 ! -----------------------------------------------------------------------
1250 ! -----------------------------------------------------------------------
1251
1252 if (use_ppm) then
1253 zt(ktop) = ze(ktop)
1254 do k = ktop + 1, kbot
1255 zt(k) = ze(k) - dt5 * (vtr(k - 1) + vtr(k))
1256 enddo
1257 zt(kbot + 1) = zs - dt * vtr(kbot)
1258
1259 do k = ktop, kbot
1260 if (zt(k + 1) >= zt(k)) zt(k + 1) = zt(k) - dz_min
1261 enddo
1262 call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qr, r1, m1_rain, mono_prof)
1263 else
1264 call implicit_fall (dt, ktop, kbot, ze, vtr, dp, qr, r1, m1_rain)
1265 endif
1266
1267 ! -----------------------------------------------------------------------
1268 ! - Calculate vertical velocity transportation during sedimentation.
1269 ! (do_sedi_w =.true. to turn on vertical motion tranport during sedimentation
1270 ! .false. by default)
1271 ! -----------------------------------------------------------------------
1272
1273 if (do_sedi_w) then
1274 w1(ktop) = (dm(ktop) * w1(ktop) + m1_rain(ktop) * vtr(ktop)) / (dm(ktop) - m1_rain(ktop))
1275 do k = ktop + 1, kbot
1276 w1(k) = (dm(k) * w1(k) - m1_rain(k - 1) * vtr(k - 1) + m1_rain(k) * vtr(k)) &
1277 / (dm(k) + m1_rain(k - 1) - m1_rain(k))
1278 enddo
1279 endif
1280
1281 ! -----------------------------------------------------------------------
1283 ! -----------------------------------------------------------------------
1284
1285 if (do_sedi_heat) &
1286 call sedi_heat (ktop, kbot, dp, m1_rain, dz, tz, qv, ql, qr, qi, qs, qg, c_liq)
1287
1288 ! -----------------------------------------------------------------------
1291 ! -----------------------------------------------------------------------
1292
1293 call revap_racc (ktop, kbot, dt5, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
1294
1295 endif
1296
1297 ! -----------------------------------------------------------------------
1301 ! -----------------------------------------------------------------------
1302
1303 if (irain_f /= 0) then
1304
1305 ! -----------------------------------------------------------------------
1306 ! no subgrid varaibility
1307 ! -----------------------------------------------------------------------
1308
1309 do k = ktop, kbot
1310 qc0 = fac_rc * ccn(k)
1311 if (tz(k) > t_wfr) then
1312 if (use_ccn) then
1313 ! -----------------------------------------------------------------------
1314 ! ccn is formulted as ccn = ccn_surface * (den / den_surface)
1315 ! -----------------------------------------------------------------------
1316 qc = qc0
1317 else
1318 qc = qc0 / den(k)
1319 endif
1320 dq = ql(k) - qc
1321 if (dq > 0.) then
1322 sink = min(dq, dt * c_praut(k) * den(k) * exp(so3 * log(ql(k))))
1323 ql(k) = ql(k) - sink
1324 qr(k) = qr(k) + sink
1325 endif
1326 endif
1327 enddo
1328
1329 else
1330
1331 ! -----------------------------------------------------------------------
1333 ! -----------------------------------------------------------------------
1334
1335 call linear_prof (kbot - ktop + 1, ql(ktop), dl(ktop), z_slope_liq, h_var)
1336
1337 do k = ktop, kbot
1338 qc0 = fac_rc * ccn(k)
1339 if (tz(k) > t_wfr + dt_fr) then
1340 dl(k) = min(max(1.e-6, dl(k)), 0.5 * ql(k))
1341 ! --------------------------------------------------------------------
1342 ! as in klein's gfdl am2 stratiform scheme (with subgrid variations)
1343 ! --------------------------------------------------------------------
1344 if (use_ccn) then
1345 ! --------------------------------------------------------------------
1346 ! ccn is formulted as ccn = ccn_surface * (den / den_surface)
1347 ! --------------------------------------------------------------------
1348 qc = qc0
1349 else
1350 qc = qc0 / den(k)
1351 endif
1352 dq = 0.5 * (ql(k) + dl(k) - qc)
1353 ! --------------------------------------------------------------------
1354 ! dq = dl if qc == q_minus = ql - dl
1355 ! dq = 0 if qc == q_plus = ql + dl
1356 ! --------------------------------------------------------------------
1357 if (dq > 0.) then ! q_plus > qc
1358 ! --------------------------------------------------------------------
1360 ! --------------------------------------------------------------------
1361 sink = min(1., dq / dl(k)) * dt * c_praut(k) * den(k) * exp(so3 * log(ql(k)))
1362 ql(k) = ql(k) - sink
1363 qr(k) = qr(k) + sink
1364 endif
1365 endif
1366 enddo
1367 endif
1368
1369end subroutine warm_rain
1370
1371! -----------------------------------------------------------------------
1375subroutine revap_racc (ktop, kbot, dt, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
1376
1377 implicit none
1378
1379 integer, intent (in) :: ktop, kbot
1380
1381 real, intent (in) :: dt ! time step (s)
1382 real, intent (in) :: rh_rain, h_var
1383
1384 real, intent (in), dimension (ktop:kbot) :: den, denfac
1385
1386 real, intent (inout), dimension (ktop:kbot) :: tz, qv, qr, ql, qi, qs, qg
1387
1388 real, dimension (ktop:kbot) :: lhl, cvm, q_liq, q_sol, lcpk
1389
1390 real :: dqv, qsat, dqsdt, evap, t2, qden, q_plus, q_minus, sink
1391 real :: qpz, dq, dqh, tin
1392
1393 integer :: k
1394
1395 do k = ktop, kbot
1396
1397 if (tz(k) > t_wfr .and. qr(k) > qrmin) then
1398
1399 ! -----------------------------------------------------------------------
1401 ! -----------------------------------------------------------------------
1402
1403 lhl(k) = lv00 + d0_vap * tz(k)
1404 q_liq(k) = ql(k) + qr(k)
1405 q_sol(k) = qi(k) + qs(k) + qg(k)
1406 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1407 lcpk(k) = lhl(k) / cvm(k)
1408
1409 tin = tz(k) - lcpk(k) * ql(k) ! presence of clouds suppresses the rain evap
1410 qpz = qv(k) + ql(k)
1411 qsat = wqs2(tin, den(k), dqsdt)
1412 dqh = max(ql(k), h_var * max(qpz, qcmin))
1413 dqh = min(dqh, 0.2 * qpz) ! new limiter
1414 dqv = qsat - qv(k) ! use this to prevent super - sat the gird box
1415 q_minus = qpz - dqh
1416 q_plus = qpz + dqh
1417
1418 ! -----------------------------------------------------------------------
1421 ! -----------------------------------------------------------------------
1422
1423 ! -----------------------------------------------------------------------
1425 ! -----------------------------------------------------------------------
1426
1427 if (dqv > qvmin .and. qsat > q_minus) then
1428 if (qsat > q_plus) then
1429 dq = qsat - qpz
1430 else
1431 ! -----------------------------------------------------------------------
1432 ! q_minus < qsat < q_plus
1433 ! dq == dqh if qsat == q_minus
1434 ! -----------------------------------------------------------------------
1435 dq = 0.25 * (q_minus - qsat) ** 2 / dqh
1436 endif
1437 qden = qr(k) * den(k)
1438 t2 = tin * tin
1439 evap = crevp(1) * t2 * dq * (crevp(2) * sqrt(qden) + crevp(3) * &
1440 exp(0.725 * log(qden))) / (crevp(4) * t2 + crevp(5) * qsat * den(k))
1441 evap = min(qr(k), dt * evap, dqv / (1. + lcpk(k) * dqsdt))
1442 ! -----------------------------------------------------------------------
1443 ! alternative minimum evap in dry environmental air
1444 ! sink = min (qr (k), dim (rh_rain * qsat, qv (k)) / (1. + lcpk (k) * dqsdt))
1445 ! evap = max (evap, sink)
1446 ! -----------------------------------------------------------------------
1447 qr(k) = qr(k) - evap
1448 qv(k) = qv(k) + evap
1449 q_liq(k) = q_liq(k) - evap
1450 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1451 tz(k) = tz(k) - evap * lhl(k) / cvm(k)
1452 endif
1453
1454 ! -----------------------------------------------------------------------
1456 ! -----------------------------------------------------------------------
1457
1458 ! if (qr (k) > qrmin .and. ql (k) > 1.e-7 .and. qsat < q_plus) then
1459 if (qr(k) > qrmin .and. ql(k) > 1.e-6 .and. qsat < q_minus) then
1460 sink = dt * denfac(k) * cracw * exp(0.95 * log(qr(k) * den(k)))
1461 sink = sink / (1. + sink) * ql(k)
1462 ql(k) = ql(k) - sink
1463 qr(k) = qr(k) + sink
1464 endif
1465
1466 endif ! warm - rain
1467 enddo
1468
1469end subroutine revap_racc
1470
1471! -----------------------------------------------------------------------
1475! qi -- > ql & ql -- > qr
1476! edges: qe == qbar + / - dm
1478subroutine linear_prof (km, q, dm, z_var, h_var)
1479
1480 implicit none
1481
1482 integer, intent (in) :: km
1483
1484 real, intent (in) :: q (km), h_var
1485
1486 real, intent (out) :: dm (km)
1487
1488 logical, intent (in) :: z_var
1489
1490 real :: dq (km)
1491
1492 integer :: k
1493
1494 if (z_var) then
1495 do k = 2, km
1496 dq(k) = 0.5 * (q(k) - q(k - 1))
1497 enddo
1498 dm(1) = 0.
1499
1500 ! -----------------------------------------------------------------------
1503 ! -----------------------------------------------------------------------
1504
1505 do k = 2, km - 1
1506 dm(k) = 0.5 * min(abs(dq(k) + dq(k + 1)), 0.5 * q(k))
1507 if (dq(k) * dq(k + 1) <= 0.) then
1508 if (dq(k) > 0.) then ! local max
1509 dm(k) = min(dm(k), dq(k), - dq(k + 1))
1510 else
1511 dm(k) = 0.
1512 endif
1513 endif
1514 enddo
1515 dm(km) = 0.
1516
1517 ! -----------------------------------------------------------------------
1520 ! -----------------------------------------------------------------------
1521
1522 do k = 1, km
1523 dm(k) = max(dm(k), qvmin, h_var * q(k))
1524 enddo
1525 else
1526 do k = 1, km
1527 dm(k) = max(qvmin, h_var * q(k))
1528 enddo
1529 endif
1530
1531end subroutine linear_prof
1532
1533! =======================================================================
1543subroutine icloud (ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, &
1544 den, denfac, vts, vtg, vtr, qak, rh_adj, rh_rain, dts, h_var)
1545
1546 implicit none
1547
1548 integer, intent (in) :: ktop, kbot
1549
1550 real, intent (in), dimension (ktop:kbot) :: p1, dp1, den, denfac, vts, vtg, vtr
1551
1552 real, intent (inout), dimension (ktop:kbot) :: tzk, qvk, qlk, qrk, qik, qsk, qgk, qak
1553
1554 real, intent (in) :: rh_adj, rh_rain, dts, h_var
1555
1556 real, dimension (ktop:kbot) :: lcpk, icpk, tcpk, di, lhl, lhi
1557 real, dimension (ktop:kbot) :: cvm, q_liq, q_sol
1558
1559 real :: rdts, fac_g2v, fac_v2g, fac_i2s, fac_imlt
1560 real :: tz, qv, ql, qr, qi, qs, qg, melt
1561 real :: pracs, psacw, pgacw, psacr, pgacr, pgaci, praci, psaci
1562 real :: pgmlt, psmlt, pgfr, pgaut, psaut, pgsub
1563 real :: tc, tsq, dqs0, qden, qim, qsm
1564 real :: dt5, factor, sink, qi_crt
1565 real :: tmp, qsw, qsi, dqsdt, dq
1566 real :: dtmp, qc, q_plus, q_minus
1567
1568 integer :: k
1569
1570 dt5 = 0.5 * dts
1571
1572 rdts = 1. / dts
1573
1574 ! -----------------------------------------------------------------------
1576 ! -----------------------------------------------------------------------
1577
1578 fac_i2s = 1. - exp(- dts / tau_i2s)
1579 fac_g2v = 1. - exp(- dts / tau_g2v)
1580 fac_v2g = 1. - exp(- dts / tau_v2g)
1581
1582 fac_imlt = 1. - exp(- dt5 / tau_imlt)
1583
1584 ! -----------------------------------------------------------------------
1586 ! -----------------------------------------------------------------------
1587
1588 do k = ktop, kbot
1589 lhi(k) = li00 + dc_ice * tzk(k)
1590 q_liq(k) = qlk(k) + qrk(k)
1591 q_sol(k) = qik(k) + qsk(k) + qgk(k)
1592 cvm(k) = c_air + qvk(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1593 icpk(k) = lhi(k) / cvm(k)
1594 enddo
1595
1596 ! -----------------------------------------------------------------------
1597 ! sources of cloud ice: pihom, cold rain, and the sat_adj
1598 ! (initiation plus deposition)
1599 ! sources of snow: cold rain, auto conversion + accretion (from cloud ice)
1600 ! sat_adj (deposition; requires pre - existing snow) ; initial snow comes from auto conversion
1601 ! -----------------------------------------------------------------------
1602
1603 do k = ktop, kbot
1604 if (tzk(k) > tice .and. qik(k) > qcmin) then
1605
1606 ! -----------------------------------------------------------------------
1608 ! -----------------------------------------------------------------------
1609
1610 melt = min(qik(k), fac_imlt * (tzk(k) - tice) / icpk(k))
1611 tmp = min(melt, dim(ql_mlt, qlk(k))) ! max ql amount
1612 qlk(k) = qlk(k) + tmp
1613 qrk(k) = qrk(k) + melt - tmp
1614 qik(k) = qik(k) - melt
1615 q_liq(k) = q_liq(k) + melt
1616 q_sol(k) = q_sol(k) - melt
1617 cvm(k) = c_air + qvk(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1618 tzk(k) = tzk(k) - melt * lhi(k) / cvm(k)
1619
1620 elseif (tzk(k) < t_wfr .and. qlk(k) > qcmin) then
1621
1622 ! -----------------------------------------------------------------------
1625 ! -----------------------------------------------------------------------
1626
1627 dtmp = t_wfr - tzk(k)
1628 factor = min(1., dtmp / dt_fr)
1629 sink = min(qlk(k) * factor, dtmp / icpk(k))
1630 qi_crt = qi_gen * min(qi_lim, 0.1 * (tice - tzk(k))) / den(k)
1631 tmp = min(sink, dim(qi_crt, qik(k)))
1632 qlk(k) = qlk(k) - sink
1633 qsk(k) = qsk(k) + sink - tmp
1634 qik(k) = qik(k) + tmp
1635 q_liq(k) = q_liq(k) - sink
1636 q_sol(k) = q_sol(k) + sink
1637 cvm(k) = c_air + qvk(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1638 tzk(k) = tzk(k) + sink * lhi(k) / cvm(k)
1639
1640 endif
1641 enddo
1642
1643 ! -----------------------------------------------------------------------
1645 ! -----------------------------------------------------------------------
1646
1647 call linear_prof (kbot - ktop + 1, qik(ktop), di(ktop), z_slope_ice, h_var)
1648
1649 ! -----------------------------------------------------------------------
1651 ! -----------------------------------------------------------------------
1652
1653 do k = ktop, kbot
1654 lhl(k) = lv00 + d0_vap * tzk(k)
1655 lhi(k) = li00 + dc_ice * tzk(k)
1656 lcpk(k) = lhl(k) / cvm(k)
1657 icpk(k) = lhi(k) / cvm(k)
1658 tcpk(k) = lcpk(k) + icpk(k)
1659 enddo
1660
1661 do k = ktop, kbot
1662
1663 ! -----------------------------------------------------------------------
1664 ! do nothing above p_min
1665 ! -----------------------------------------------------------------------
1666
1667 if (p1(k) < p_min) cycle
1668
1669 tz = tzk(k)
1670 qv = qvk(k)
1671 ql = qlk(k)
1672 qi = qik(k)
1673 qr = qrk(k)
1674 qs = qsk(k)
1675 qg = qgk(k)
1676
1677 pgacr = 0.
1678 pgacw = 0.
1679 tc = tz - tice
1680
1681 if (tc .ge. 0.) then
1682
1683 ! -----------------------------------------------------------------------
1685 ! -----------------------------------------------------------------------
1686
1687 dqs0 = ces0 / p1(k) - qv
1688
1689 if (qs > qcmin) then
1690
1691 ! -----------------------------------------------------------------------
1693 ! only rate is used (for snow melt) since tc > 0.
1694 ! -----------------------------------------------------------------------
1695
1696 if (ql > qrmin) then
1697 factor = denfac(k) * csacw * exp(0.8125 * log(qs * den(k)))
1698 psacw = factor / (1. + dts * factor) * ql ! rate
1699 else
1700 psacw = 0.
1701 endif
1702
1703 ! -----------------------------------------------------------------------
1706 ! -----------------------------------------------------------------------
1707
1708 if (qr > qrmin) then
1709 psacr = min(acr3d(vts(k), vtr(k), qr, qs, csacr, acco(1, 2), &
1710 den(k)), qr * rdts)
1711 pracs = acr3d(vtr(k), vts(k), qs, qr, cracs, acco(1, 1), den(k))
1712 else
1713 psacr = 0.
1714 pracs = 0.
1715 endif
1716
1717 ! -----------------------------------------------------------------------
1720 ! -----------------------------------------------------------------------
1721
1722 psmlt = max(0., smlt(tc, dqs0, qs * den(k), psacw, psacr, csmlt, &
1723 den(k), denfac(k)))
1724 sink = min(qs, dts * (psmlt + pracs), tc / icpk(k))
1725 qs = qs - sink
1726 ! sjl, 20170321:
1727 tmp = min(sink, dim(qs_mlt, ql)) ! max ql due to snow melt
1728 ql = ql + tmp
1729 qr = qr + sink - tmp
1730 ! qr = qr + sink
1731 ! sjl, 20170321:
1732 q_liq(k) = q_liq(k) + sink
1733 q_sol(k) = q_sol(k) - sink
1734 cvm(k) = c_air + qv * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1735 tz = tz - sink * lhi(k) / cvm(k)
1736 tc = tz - tice
1737
1738 endif
1739
1740 ! -----------------------------------------------------------------------
1742 ! -----------------------------------------------------------------------
1743
1744 lhi(k) = li00 + dc_ice * tz
1745 icpk(k) = lhi(k) / cvm(k)
1746
1747 ! -----------------------------------------------------------------------
1749 ! -----------------------------------------------------------------------
1750
1751 if (qg > qcmin .and. tc > 0.) then
1752
1753 ! -----------------------------------------------------------------------
1755 ! -----------------------------------------------------------------------
1756
1757 if (qr > qrmin) &
1758 pgacr = min(acr3d(vtg(k), vtr(k), qr, qg, cgacr, acco(1, 3), &
1759 den(k)), rdts * qr)
1760
1761 ! -----------------------------------------------------------------------
1763 ! -----------------------------------------------------------------------
1764
1765 qden = qg * den(k)
1766 if (ql > qrmin) then
1767 factor = cgacw * qden / sqrt(den(k) * sqrt(sqrt(qden)))
1768 pgacw = factor / (1. + dts * factor) * ql ! rate
1769 endif
1770
1771 ! -----------------------------------------------------------------------
1773 ! -----------------------------------------------------------------------
1774
1775 pgmlt = dts * gmlt(tc, dqs0, qden, pgacw, pgacr, cgmlt, den(k))
1776 pgmlt = min(max(0., pgmlt), qg, tc / icpk(k))
1777 qg = qg - pgmlt
1778 qr = qr + pgmlt
1779 q_liq(k) = q_liq(k) + pgmlt
1780 q_sol(k) = q_sol(k) - pgmlt
1781 cvm(k) = c_air + qv * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1782 tz = tz - pgmlt * lhi(k) / cvm(k)
1783
1784 endif
1785
1786 else
1787
1788 ! -----------------------------------------------------------------------
1790 ! -----------------------------------------------------------------------
1791
1792 ! -----------------------------------------------------------------------
1794 ! -----------------------------------------------------------------------
1795
1796 if (qi > 3.e-7) then ! cloud ice sink terms
1797
1798 if (qs > 1.e-7) then
1799 ! -----------------------------------------------------------------------
1800 ! sjl added (following lin eq. 23) the temperature dependency
1801 ! to reduce accretion, use esi = exp (0.05 * tc) as in hong et al 2004
1802 ! -----------------------------------------------------------------------
1803 factor = dts * denfac(k) * csaci * exp(0.05 * tc + 0.8125 * log(qs * den(k)))
1804 psaci = factor / (1. + factor) * qi
1805 else
1806 psaci = 0.
1807 endif
1808
1809 ! -----------------------------------------------------------------------
1814 ! -----------------------------------------------------------------------
1815 if (qi0_crt < 0.) then
1816 qim = - qi0_crt
1817 else
1818 qim = qi0_crt / den(k)
1819 endif
1820 ! -----------------------------------------------------------------------
1821 ! assuming linear subgrid vertical distribution of cloud ice
1822 ! the mismatch computation following lin et al. 1994, mwr
1823 ! -----------------------------------------------------------------------
1824
1825 if (const_vi) then
1826 tmp = fac_i2s
1827 else
1828 tmp = fac_i2s * exp(0.025 * tc)
1829 endif
1830
1831 di(k) = max(di(k), qrmin)
1832 q_plus = qi + di(k)
1833 if (q_plus > (qim + qrmin)) then
1834 if (qim > (qi - di(k))) then
1835 dq = (0.25 * (q_plus - qim) ** 2) / di(k)
1836 else
1837 dq = qi - qim
1838 endif
1839 psaut = tmp * dq
1840 else
1841 psaut = 0.
1842 endif
1843 ! -----------------------------------------------------------------------
1844 ! sink is no greater than 75% of qi
1845 ! -----------------------------------------------------------------------
1846 sink = min(0.75 * qi, psaci + psaut)
1847 qi = qi - sink
1848 qs = qs + sink
1849
1850 ! -----------------------------------------------------------------------
1852 ! -----------------------------------------------------------------------
1853
1854 if (qg > 1.e-6) then
1855 ! -----------------------------------------------------------------------
1856 ! factor = dts * cgaci / sqrt (den (k)) * exp (0.05 * tc + 0.875 * log (qg * den (k)))
1857 ! simplified form: remove temp dependency & set the exponent "0.875" -- > 1
1858 ! -----------------------------------------------------------------------
1859 factor = dts * cgaci * sqrt(den(k)) * qg
1860 pgaci = factor / (1. + factor) * qi
1861 qi = qi - pgaci
1862 qg = qg + pgaci
1863 endif
1864
1865 endif
1866
1867 ! -----------------------------------------------------------------------
1869 ! -----------------------------------------------------------------------
1870
1871 ! -----------------------------------------------------------------------
1872 ! rain to ice, snow, graupel processes:
1873 ! -----------------------------------------------------------------------
1874
1875 tc = tz - tice
1876
1877 if (qr > 1.e-7 .and. tc < 0.) then
1878
1879 ! -----------------------------------------------------------------------
1880 ! * sink * terms to qr: psacr + pgfr
1881 ! source terms to qs: psacr
1882 ! source terms to qg: pgfr
1883 ! -----------------------------------------------------------------------
1884
1885 ! -----------------------------------------------------------------------
1887 ! -----------------------------------------------------------------------
1888
1889 if (qs > 1.e-7) then ! if snow exists
1890 psacr = dts * acr3d(vts(k), vtr(k), qr, qs, csacr, acco(1, 2), den(k))
1891 else
1892 psacr = 0.
1893 endif
1894
1895 ! -----------------------------------------------------------------------
1897 ! -----------------------------------------------------------------------
1898
1899 pgfr = dts * cgfr(1) / den(k) * (exp(- cgfr(2) * tc) - 1.) * &
1900 exp(1.75 * log(qr * den(k)))
1901
1902 ! -----------------------------------------------------------------------
1907 ! -----------------------------------------------------------------------
1908
1909 sink = psacr + pgfr
1910 factor = min(sink, qr, - tc / icpk(k)) / max(sink, qrmin)
1911
1912 psacr = factor * psacr
1913 pgfr = factor * pgfr
1914
1915 sink = psacr + pgfr
1916 qr = qr - sink
1917 qs = qs + psacr
1918 qg = qg + pgfr
1919 q_liq(k) = q_liq(k) - sink
1920 q_sol(k) = q_sol(k) + sink
1921 cvm(k) = c_air + qv * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
1922 tz = tz + sink * lhi(k) / cvm(k)
1923
1924 endif
1925
1926 ! -----------------------------------------------------------------------
1928 ! -----------------------------------------------------------------------
1929
1930 lhi(k) = li00 + dc_ice * tz
1931 icpk(k) = lhi(k) / cvm(k)
1932
1933 ! -----------------------------------------------------------------------
1935 ! -----------------------------------------------------------------------
1936
1937 if (qs > 1.e-7) then
1938
1939 ! -----------------------------------------------------------------------
1941 ! -----------------------------------------------------------------------
1942
1943 if (qg > qrmin) then
1944 sink = dts * acr3d(vtg(k), vts(k), qs, qg, cgacs, acco(1, 4), den(k))
1945 else
1946 sink = 0.
1947 endif
1948
1949 ! -----------------------------------------------------------------------
1951 ! -----------------------------------------------------------------------
1952
1953 qsm = qs0_crt / den(k)
1954 if (qs > qsm) then
1955 factor = dts * 1.e-3 * exp(0.09 * (tz - tice))
1956 sink = sink + factor / (1. + factor) * (qs - qsm)
1957 endif
1958 sink = min(qs, sink)
1959 qs = qs - sink
1960 qg = qg + sink
1961
1962 endif ! snow existed
1963
1964 if (qg > 1.e-7 .and. tz < tice0) then
1965
1966 ! -----------------------------------------------------------------------
1968 ! -----------------------------------------------------------------------
1969
1970 if (ql > 1.e-6) then
1971 qden = qg * den(k)
1972 factor = dts * cgacw * qden / sqrt(den(k) * sqrt(sqrt(qden)))
1973 pgacw = factor / (1. + factor) * ql
1974 else
1975 pgacw = 0.
1976 endif
1977
1978 ! -----------------------------------------------------------------------
1980 ! -----------------------------------------------------------------------
1981
1982 if (qr > 1.e-6) then
1983 pgacr = min(dts * acr3d(vtg(k), vtr(k), qr, qg, cgacr, acco(1, 3), &
1984 den(k)), qr)
1985 else
1986 pgacr = 0.
1987 endif
1988
1989 sink = pgacr + pgacw
1990 factor = min(sink, dim(tice, tz) / icpk(k)) / max(sink, qrmin)
1991 pgacr = factor * pgacr
1992 pgacw = factor * pgacw
1993
1994 sink = pgacr + pgacw
1995 qg = qg + sink
1996 qr = qr - pgacr
1997 ql = ql - pgacw
1998 q_liq(k) = q_liq(k) - sink
1999 q_sol(k) = q_sol(k) + sink
2000 cvm(k) = c_air + qv * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2001 tz = tz + sink * lhi(k) / cvm(k)
2002
2003 endif
2004
2005 endif
2006
2007 tzk(k) = tz
2008 qvk(k) = qv
2009 qlk(k) = ql
2010 qik(k) = qi
2011 qrk(k) = qr
2012 qsk(k) = qs
2013 qgk(k) = qg
2014
2015 enddo
2016
2017 ! -----------------------------------------------------------------------
2019 ! -----------------------------------------------------------------------
2020
2021 call subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, rh_adj, tzk, qvk, &
2022 qlk, qrk, qik, qsk, qgk, qak, h_var, rh_rain)
2023
2024end subroutine icloud
2025
2026! =======================================================================
2030subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, rh_adj, tz, qv, &
2031 ql, qr, qi, qs, qg, qa, h_var, rh_rain)
2032
2033 implicit none
2034
2035 integer, intent (in) :: ktop, kbot
2036
2037 real, intent (in), dimension (ktop:kbot) :: p1, den, denfac
2038
2039 real, intent (in) :: dts, rh_adj, h_var, rh_rain
2040
2041 real, intent (inout), dimension (ktop:kbot) :: tz, qv, ql, qr, qi, qs, qg, qa
2042
2043 real, dimension (ktop:kbot) :: lcpk, icpk, tcpk, tcp3, lhl, lhi
2044 real, dimension (ktop:kbot) :: cvm, q_liq, q_sol, q_cond
2045
2046 real :: fac_v2l, fac_l2v
2047
2048 real :: pidep, qi_crt
2049
2050 ! -----------------------------------------------------------------------
2051 ! qstar over water may be accurate only down to - 80 deg c with ~10% uncertainty
2052 ! must not be too large to allow psc
2053 ! -----------------------------------------------------------------------
2054
2055 real :: rh, rqi, tin, qsw, qsi, qpz, qstar
2056 real :: dqsdt, dwsdt, dq, dq0, factor, tmp
2057 real :: q_plus, q_minus, dt_evap, dt_pisub
2058 real :: evap, sink, tc, pisub, q_adj, dtmp
2059 real :: pssub, pgsub, tsq, qden, fac_g2v, fac_v2g
2060
2061 integer :: k
2062
2063 if (fast_sat_adj) then
2064 dt_evap = 0.5 * dts
2065 else
2066 dt_evap = dts
2067 endif
2068
2069 ! -----------------------------------------------------------------------
2071 ! -----------------------------------------------------------------------
2072
2073 fac_v2l = 1. - exp(- dt_evap / tau_v2l)
2074 fac_l2v = 1. - exp(- dt_evap / tau_l2v)
2075
2076 fac_g2v = 1. - exp(- dts / tau_g2v)
2077 fac_v2g = 1. - exp(- dts / tau_v2g)
2078
2079 ! -----------------------------------------------------------------------
2081 ! -----------------------------------------------------------------------
2082
2083 do k = ktop, kbot
2084 lhl(k) = lv00 + d0_vap * tz(k)
2085 lhi(k) = li00 + dc_ice * tz(k)
2086 q_liq(k) = ql(k) + qr(k)
2087 q_sol(k) = qi(k) + qs(k) + qg(k)
2088 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2089 lcpk(k) = lhl(k) / cvm(k)
2090 icpk(k) = lhi(k) / cvm(k)
2091 tcpk(k) = lcpk(k) + icpk(k)
2092 tcp3(k) = lcpk(k) + icpk(k) * min(1., dim(tice, tz(k)) / (tice - t_wfr))
2093 enddo
2094
2095 do k = ktop, kbot
2096
2097 if (p1(k) < p_min) cycle
2098
2099 ! -----------------------------------------------------------------------
2101 ! -----------------------------------------------------------------------
2102
2103 if (tz(k) < t_min) then
2104 sink = dim(qv(k), 1.e-7)
2105 qv(k) = qv(k) - sink
2106 qi(k) = qi(k) + sink
2107 q_sol(k) = q_sol(k) + sink
2108 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2109 tz(k) = tz(k) + sink * (lhl(k) + lhi(k)) / cvm(k)
2110 if (.not. do_qa) qa(k) = qa(k) + 1. ! air fully saturated; 100 % cloud cover
2111 cycle
2112 endif
2113
2114 ! -----------------------------------------------------------------------
2116 ! -----------------------------------------------------------------------
2117
2118 lhl(k) = lv00 + d0_vap * tz(k)
2119 lhi(k) = li00 + dc_ice * tz(k)
2120 lcpk(k) = lhl(k) / cvm(k)
2121 icpk(k) = lhi(k) / cvm(k)
2122 tcpk(k) = lcpk(k) + icpk(k)
2123 tcp3(k) = lcpk(k) + icpk(k) * min(1., dim(tice, tz(k)) / (tice - t_wfr))
2124
2125 ! -----------------------------------------------------------------------
2127 ! -----------------------------------------------------------------------
2128
2129 qpz = qv(k) + ql(k) + qi(k)
2130 tin = tz(k) - (lhl(k) * (ql(k) + qi(k)) + lhi(k) * qi(k)) / (c_air + &
2131 qpz * c_vap + qr(k) * c_liq + (qs(k) + qg(k)) * c_ice)
2132 if (tin > t_sub + 6.) then
2133 rh = qpz / iqs1(tin, den(k))
2134 if (rh < rh_adj) then ! qpz / rh_adj < qs
2135 tz(k) = tin
2136 qv(k) = qpz
2137 ql(k) = 0.
2138 qi(k) = 0.
2139 cycle ! cloud free
2140 endif
2141 endif
2142
2143 ! -----------------------------------------------------------------------
2145 ! -----------------------------------------------------------------------
2146
2147 qsw = wqs2(tz(k), den(k), dwsdt)
2148 dq0 = qsw - qv(k)
2149 if (dq0 > 0.) then
2150 ! SJL 20170703 added ql factor to prevent the situation of high ql and low RH
2151 ! factor = min (1., fac_l2v * sqrt (max (0., ql (k)) / 1.e-5) * 10. * dq0 / qsw)
2152 ! factor = fac_l2v
2153 ! factor = 1
2154 factor = min(1., fac_l2v * (10. * dq0 / qsw)) ! the rh dependent factor = 1 at 90%
2155 evap = min(ql(k), factor * dq0 / (1. + tcp3(k) * dwsdt))
2156 else ! condensate all excess vapor into cloud water
2157 ! -----------------------------------------------------------------------
2158 ! evap = fac_v2l * dq0 / (1. + tcp3 (k) * dwsdt)
2159 ! sjl, 20161108
2160 ! -----------------------------------------------------------------------
2161 evap = dq0 / (1. + tcp3(k) * dwsdt)
2162 endif
2163 qv(k) = qv(k) + evap
2164 ql(k) = ql(k) - evap
2165 q_liq(k) = q_liq(k) - evap
2166 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2167 tz(k) = tz(k) - evap * lhl(k) / cvm(k)
2168
2169 ! -----------------------------------------------------------------------
2171 ! -----------------------------------------------------------------------
2172
2173 lhi(k) = li00 + dc_ice * tz(k)
2174 icpk(k) = lhi(k) / cvm(k)
2175
2176 ! -----------------------------------------------------------------------
2178 ! -----------------------------------------------------------------------
2179
2180 dtmp = t_wfr - tz(k) ! [ - 40, - 48]
2181 if (dtmp > 0. .and. ql(k) > qcmin) then
2182 sink = min(ql(k), ql(k) * dtmp * 0.125, dtmp / icpk(k))
2183 ql(k) = ql(k) - sink
2184 qi(k) = qi(k) + sink
2185 q_liq(k) = q_liq(k) - sink
2186 q_sol(k) = q_sol(k) + sink
2187 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2188 tz(k) = tz(k) + sink * lhi(k) / cvm(k)
2189 endif
2190
2191 ! -----------------------------------------------------------------------
2193 ! -----------------------------------------------------------------------
2194
2195 lhi(k) = li00 + dc_ice * tz(k)
2196 icpk(k) = lhi(k) / cvm(k)
2197
2198 ! -----------------------------------------------------------------------
2200 ! -----------------------------------------------------------------------
2201
2202 if (fast_sat_adj) then
2203 dt_pisub = 0.5 * dts
2204 else
2205 dt_pisub = dts
2206 tc = tice - tz(k)
2207 if (ql(k) > qrmin .and. tc > 0.) then
2208 sink = 3.3333e-10 * dts * (exp(0.66 * tc) - 1.) * den(k) * ql(k) * ql(k)
2209 sink = min(ql(k), tc / icpk(k), sink)
2210 ql(k) = ql(k) - sink
2211 qi(k) = qi(k) + sink
2212 q_liq(k) = q_liq(k) - sink
2213 q_sol(k) = q_sol(k) + sink
2214 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2215 tz(k) = tz(k) + sink * lhi(k) / cvm(k)
2216 endif ! significant ql existed
2217 endif
2218
2219 ! -----------------------------------------------------------------------
2221 ! -----------------------------------------------------------------------
2222
2223 lhl(k) = lv00 + d0_vap * tz(k)
2224 lhi(k) = li00 + dc_ice * tz(k)
2225 lcpk(k) = lhl(k) / cvm(k)
2226 icpk(k) = lhi(k) / cvm(k)
2227 tcpk(k) = lcpk(k) + icpk(k)
2228
2229 ! -----------------------------------------------------------------------
2231 ! -----------------------------------------------------------------------
2232
2233 if (tz(k) < tice) then
2234 qsi = iqs2(tz(k), den(k), dqsdt)
2235 dq = qv(k) - qsi
2236 sink = dq / (1. + tcpk(k) * dqsdt)
2237 if (qi(k) > qrmin) then
2238 ! eq 9, hong et al. 2004, mwr
2239 ! for a and b, see dudhia 1989: page 3103 eq (b7) and (b8)
2240 pidep = dt_pisub * dq * 349138.78 * exp(0.875 * log(qi(k) * den(k))) &
2241 / (qsi * den(k) * lat2 / (0.0243 * rvgas * tz(k) ** 2) + 4.42478e4)
2242 else
2243 pidep = 0.
2244 endif
2245 if (dq > 0.) then ! vapor - > ice
2246 tmp = tice - tz(k)
2247 ! 20160912: the following should produce more ice at higher altitude
2248 ! qi_crt = 4.92e-11 * exp (1.33 * log (1.e3 * exp (0.1 * tmp))) / den (k)
2249 qi_crt = qi_gen * min(qi_lim, 0.1 * tmp) / den(k)
2250 sink = min(sink, max(qi_crt - qi(k), pidep), tmp / tcpk(k))
2251 else ! ice -- > vapor
2252 pidep = pidep * min(1., dim(tz(k), t_sub) * 0.2)
2253 sink = max(pidep, sink, - qi(k))
2254 endif
2255 qv(k) = qv(k) - sink
2256 qi(k) = qi(k) + sink
2257 q_sol(k) = q_sol(k) + sink
2258 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2259 tz(k) = tz(k) + sink * (lhl(k) + lhi(k)) / cvm(k)
2260 endif
2261
2262 ! -----------------------------------------------------------------------
2264 ! -----------------------------------------------------------------------
2265
2266 lhl(k) = lv00 + d0_vap * tz(k)
2267 lhi(k) = li00 + dc_ice * tz(k)
2268 lcpk(k) = lhl(k) / cvm(k)
2269 icpk(k) = lhi(k) / cvm(k)
2270 tcpk(k) = lcpk(k) + icpk(k)
2271
2272 ! -----------------------------------------------------------------------
2274 ! this process happens for all temp rage
2275 ! -----------------------------------------------------------------------
2276
2277 if (qs(k) > qrmin) then
2278 qsi = iqs2(tz(k), den(k), dqsdt)
2279 qden = qs(k) * den(k)
2280 tmp = exp(0.65625 * log(qden))
2281 tsq = tz(k) * tz(k)
2282 dq = (qsi - qv(k)) / (1. + tcpk(k) * dqsdt)
2283 pssub = cssub(1) * tsq * (cssub(2) * sqrt(qden) + cssub(3) * tmp * &
2284 sqrt(denfac(k))) / (cssub(4) * tsq + cssub(5) * qsi * den(k))
2285 pssub = (qsi - qv(k)) * dts * pssub
2286 if (pssub > 0.) then ! qs -- > qv, sublimation
2287 pssub = min(pssub * min(1., dim(tz(k), t_sub) * 0.2), qs(k))
2288 else
2289 if (tz(k) > tice) then
2290 pssub = 0. ! no deposition
2291 else
2292 pssub = max(pssub, dq, (tz(k) - tice) / tcpk(k))
2293 endif
2294 endif
2295 qs(k) = qs(k) - pssub
2296 qv(k) = qv(k) + pssub
2297 q_sol(k) = q_sol(k) - pssub
2298 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2299 tz(k) = tz(k) - pssub * (lhl(k) + lhi(k)) / cvm(k)
2300 endif
2301
2302 ! -----------------------------------------------------------------------
2304 ! -----------------------------------------------------------------------
2305
2306 lhl(k) = lv00 + d0_vap * tz(k)
2307 lhi(k) = li00 + dc_ice * tz(k)
2308 lcpk(k) = lhl(k) / cvm(k)
2309 icpk(k) = lhi(k) / cvm(k)
2310 tcpk(k) = lcpk(k) + icpk(k)
2311
2312 ! -----------------------------------------------------------------------
2314 ! -----------------------------------------------------------------------
2315
2316 if (qg(k) > qrmin) then
2317 qsi = iqs2(tz(k), den(k), dqsdt)
2318 dq = (qv(k) - qsi) / (1. + tcpk(k) * dqsdt)
2319 pgsub = (qv(k) / qsi - 1.) * qg(k)
2320 if (pgsub > 0.) then ! deposition
2321 if (tz(k) > tice) then
2322 pgsub = 0. ! no deposition
2323 else
2324 pgsub = min(fac_v2g * pgsub, 0.2 * dq, ql(k) + qr(k), &
2325 (tice - tz(k)) / tcpk(k))
2326 endif
2327 else ! submilation
2328 pgsub = max(fac_g2v * pgsub, dq) * min(1., dim(tz(k), t_sub) * 0.1)
2329 endif
2330 qg(k) = qg(k) + pgsub
2331 qv(k) = qv(k) - pgsub
2332 q_sol(k) = q_sol(k) + pgsub
2333 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2334 tz(k) = tz(k) + pgsub * (lhl(k) + lhi(k)) / cvm(k)
2335 endif
2336
2337#ifdef USE_MIN_EVAP
2338 ! -----------------------------------------------------------------------
2340 ! -----------------------------------------------------------------------
2341
2342 lhl(k) = lv00 + d0_vap * tz(k)
2343 lcpk(k) = lhl(k) / cvm(k)
2344
2345 ! -----------------------------------------------------------------------
2347 ! -----------------------------------------------------------------------
2348
2349 if (qr(k) > qcmin) then
2350 qsw = wqs2(tz(k), den(k), dqsdt)
2351 sink = min(qr(k), dim(rh_rain * qsw, qv(k)) / (1. + lcpk(k) * dqsdt))
2352 qv(k) = qv(k) + sink
2353 qr(k) = qr(k) - sink
2354 q_liq(k) = q_liq(k) - sink
2355 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2356 tz(k) = tz(k) - sink * lhl(k) / cvm(k)
2357 endif
2358#endif
2359
2360 ! -----------------------------------------------------------------------
2362 ! -----------------------------------------------------------------------
2363
2364 lhl(k) = lv00 + d0_vap * tz(k)
2365 cvm(k) = c_air + (qv(k) + q_liq(k) + q_sol(k)) * c_vap
2366 lcpk(k) = lhl(k) / cvm(k)
2367
2368 ! -----------------------------------------------------------------------
2369 ! compute cloud fraction
2370 ! -----------------------------------------------------------------------
2371
2372 ! -----------------------------------------------------------------------
2374 ! -----------------------------------------------------------------------
2375
2376 if (do_qa) cycle
2377
2378 if (rad_snow) then
2379 q_sol(k) = qi(k) + qs(k)
2380 else
2381 q_sol(k) = qi(k)
2382 endif
2383 if (rad_rain) then
2384 q_liq(k) = ql(k) + qr(k)
2385 else
2386 q_liq(k) = ql(k)
2387 endif
2388 q_cond(k) = q_liq(k) + q_sol(k)
2389
2390 qpz = qv(k) + q_cond(k) ! qpz is conserved
2391
2392 ! -----------------------------------------------------------------------
2394 ! -----------------------------------------------------------------------
2395
2396 tin = tz(k) - (lcpk(k) * q_cond(k) + icpk(k) * q_sol(k)) ! minimum temperature
2397 ! tin = tz (k) - ((lv00 + d0_vap * tz (k)) * q_cond (k) + &
2398 ! (li00 + dc_ice * tz (k)) * q_sol (k)) / (c_air + qpz * c_vap)
2399
2400 ! -----------------------------------------------------------------------
2401 ! determine saturated specific humidity
2402 ! -----------------------------------------------------------------------
2403
2404 if (tin <= t_wfr) then
2405 ! ice phase:
2406 qstar = iqs1(tin, den(k))
2407 elseif (tin >= tice) then
2408 ! liquid phase:
2409 qstar = wqs1(tin, den(k))
2410 else
2411 ! mixed phase:
2412 qsi = iqs1(tin, den(k))
2413 qsw = wqs1(tin, den(k))
2414 if (q_cond(k) > 3.e-6) then
2415 rqi = q_sol(k) / q_cond(k)
2416 else
2417 ! -----------------------------------------------------------------------
2418 ! mostly liquid water q_cond (k) at initial cloud development stage
2419 ! -----------------------------------------------------------------------
2420 rqi = (tice - tin) / (tice - t_wfr)
2421 endif
2422 qstar = rqi * qsi + (1. - rqi) * qsw
2423 endif
2424
2425 ! -----------------------------------------------------------------------
2428 ! -----------------------------------------------------------------------
2429
2430 if (qpz > qrmin) then
2431 ! partial cloudiness by pdf:
2432 dq = max(qcmin, h_var * qpz)
2433 q_plus = qpz + dq ! cloud free if qstar > q_plus
2434 q_minus = qpz - dq
2435 if (qstar < q_minus) then
2436 qa(k) = qa(k) + 1. ! air fully saturated; 100 % cloud cover
2437 elseif (qstar < q_plus .and. q_cond(k) > qc_crt) then
2438 qa(k) = qa(k) + (q_plus - qstar) / (dq + dq) ! partial cloud cover
2439 ! qa (k) = sqrt (qa (k) + (q_plus - qstar) / (dq + dq))
2440 endif
2441 endif
2442
2443 enddo
2444
2445end subroutine subgrid_z_proc
2446
2447! =======================================================================
2450subroutine revap_rac1 (hydrostatic, is, ie, dt, tz, qv, ql, qr, qi, qs, qg, den, hvar)
2451
2452 implicit none
2453
2454 logical, intent (in) :: hydrostatic
2455
2456 integer, intent (in) :: is, ie
2457
2458 real, intent (in) :: dt ! time step (s)
2459
2460 real, intent (in), dimension (is:ie) :: den, hvar, qi, qs, qg
2461
2462 real, intent (inout), dimension (is:ie) :: tz, qv, qr, ql
2463
2464 real, dimension (is:ie) :: lcp2, denfac, q_liq, q_sol, cvm, lhl
2465
2466 real :: dqv, qsat, dqsdt, evap, qden, q_plus, q_minus, sink
2467 real :: tin, t2, qpz, dq, dqh
2468
2469 integer :: i
2470
2471 ! -----------------------------------------------------------------------
2472 ! define latend heat coefficient
2473 ! -----------------------------------------------------------------------
2474
2475 do i = is, ie
2476 lhl(i) = lv00 + d0_vap * tz(i)
2477 q_liq(i) = ql(i) + qr(i)
2478 q_sol(i) = qi(i) + qs(i) + qg(i)
2479 cvm(i) = c_air + qv(i) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
2480 lcp2(i) = lhl(i) / cvm(i)
2481 ! denfac (i) = sqrt (sfcrho / den (i))
2482 enddo
2483
2484 do i = is, ie
2485 if (qr(i) > qrmin .and. tz(i) > t_wfr) then
2486 qpz = qv(i) + ql(i)
2487 tin = tz(i) - lcp2(i) * ql(i) ! presence of clouds suppresses the rain evap
2488 qsat = wqs2(tin, den(i), dqsdt)
2489 dqh = max(ql(i), hvar(i) * max(qpz, qcmin))
2490 dqv = qsat - qv(i)
2491 q_minus = qpz - dqh
2492 q_plus = qpz + dqh
2493
2494 ! -----------------------------------------------------------------------
2495 ! qsat must be > q_minus to activate evaporation
2496 ! qsat must be < q_plus to activate accretion
2497 ! -----------------------------------------------------------------------
2498
2499 ! -----------------------------------------------------------------------
2500 ! rain evaporation
2501 ! -----------------------------------------------------------------------
2502
2503 if (dqv > qvmin .and. qsat > q_minus) then
2504 if (qsat > q_plus) then
2505 dq = qsat - qpz
2506 else
2507 ! q_minus < qsat < q_plus
2508 ! dq == dqh if qsat == q_minus
2509 dq = 0.25 * (q_minus - qsat) ** 2 / dqh
2510 endif
2511 qden = qr(i) * den(i)
2512 t2 = tin * tin
2513 evap = crevp(1) * t2 * dq * (crevp(2) * sqrt(qden) + crevp(3) * exp(0.725 * log(qden))) &
2514 / (crevp(4) * t2 + crevp(5) * qsat * den(i))
2515 evap = min(qr(i), dt * evap, dqv / (1. + lcp2(i) * dqsdt))
2516 qr(i) = qr(i) - evap
2517 qv(i) = qv(i) + evap
2518 q_liq(i) = q_liq(i) - evap
2519 cvm(i) = c_air + qv(i) * c_vap + q_liq(i) * c_liq + q_sol(i) * c_ice
2520 tz(i) = tz(i) - evap * lhl(i) / cvm(i)
2521 endif
2522
2523 ! -----------------------------------------------------------------------
2524 ! accretion: pracc
2525 ! -----------------------------------------------------------------------
2526
2527 if (qr(i) > qrmin .and. ql(i) > 1.e-8 .and. qsat < q_plus) then
2528 denfac(i) = sqrt(sfcrho / den(i))
2529 sink = dt * denfac(i) * cracw * exp(0.95 * log(qr(i) * den(i)))
2530 sink = sink / (1. + sink) * ql(i)
2531 ql(i) = ql(i) - sink
2532 qr(i) = qr(i) + sink
2533 endif
2534 endif
2535 enddo
2536
2537end subroutine revap_rac1
2538
2539! =======================================================================
2543subroutine terminal_fall (dtm, ktop, kbot, tz, qv, ql, qr, qg, qs, qi, dz, dp, &
2544 den, vtg, vts, vti, r1, g1, s1, i1, m1_sol, w1)
2545
2546 implicit none
2547
2548 integer, intent (in) :: ktop, kbot
2549
2550 real, intent (in) :: dtm ! time step (s)
2551
2552 real, intent (in), dimension (ktop:kbot) :: vtg, vts, vti, den, dp, dz
2553
2554 real, intent (inout), dimension (ktop:kbot) :: qv, ql, qr, qg, qs, qi, tz, m1_sol, w1
2555
2556 real, intent (out) :: r1, g1, s1, i1
2557
2558 real, dimension (ktop:kbot + 1) :: ze, zt
2559
2560 real :: qsat, dqsdt, dt5, evap, dtime
2561 real :: factor, frac
2562 real :: tmp, precip, tc, sink
2563
2564 real, dimension (ktop:kbot) :: lcpk, icpk, cvm, q_liq, q_sol, lhl, lhi
2565 real, dimension (ktop:kbot) :: m1, dm
2566
2567 real :: zs = 0.
2568 real :: fac_imlt
2569
2570 integer :: k, k0, m
2571
2572 logical :: no_fall
2573
2574 dt5 = 0.5 * dtm
2575 fac_imlt = 1. - exp(- dt5 / tau_imlt)
2576
2577 ! -----------------------------------------------------------------------
2578 ! define heat capacity and latend heat coefficient
2579 ! -----------------------------------------------------------------------
2580
2581 do k = ktop, kbot
2582 m1_sol(k) = 0.
2583 lhl(k) = lv00 + d0_vap * tz(k)
2584 lhi(k) = li00 + dc_ice * tz(k)
2585 q_liq(k) = ql(k) + qr(k)
2586 q_sol(k) = qi(k) + qs(k) + qg(k)
2587 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2588 lcpk(k) = lhl(k) / cvm(k)
2589 icpk(k) = lhi(k) / cvm(k)
2590 enddo
2591
2592 ! -----------------------------------------------------------------------
2593 ! find significant melting level
2594 ! -----------------------------------------------------------------------
2595
2596 k0 = kbot
2597 do k = ktop, kbot - 1
2598 if (tz(k) > tice) then
2599 k0 = k
2600 exit
2601 endif
2602 enddo
2603
2604 ! -----------------------------------------------------------------------
2605 ! melting of cloud_ice (before fall) :
2606 ! -----------------------------------------------------------------------
2607
2608 do k = k0, kbot
2609 tc = tz(k) - tice
2610 if (qi(k) > qcmin .and. tc > 0.) then
2611 sink = min(qi(k), fac_imlt * tc / icpk(k))
2612 tmp = min(sink, dim(ql_mlt, ql(k)))
2613 ql(k) = ql(k) + tmp
2614 qr(k) = qr(k) + sink - tmp
2615 qi(k) = qi(k) - sink
2616 q_liq(k) = q_liq(k) + sink
2617 q_sol(k) = q_sol(k) - sink
2618 cvm(k) = c_air + qv(k) * c_vap + q_liq(k) * c_liq + q_sol(k) * c_ice
2619 tz(k) = tz(k) - sink * lhi(k) / cvm(k)
2620 tc = tz(k) - tice
2621 endif
2622 enddo
2623
2624 ! -----------------------------------------------------------------------
2625 ! turn off melting when cloud microphysics time step is small
2626 ! -----------------------------------------------------------------------
2627
2628 if (dtm < 60.) k0 = kbot
2629
2630 ! sjl, turn off melting of falling cloud ice, snow and graupel
2631 k0 = kbot
2632 ! sjl, turn off melting of falling cloud ice, snow and graupel
2633
2634 ze(kbot + 1) = zs
2635 do k = kbot, ktop, - 1
2636 ze(k) = ze(k + 1) - dz(k) ! dz < 0
2637 enddo
2638
2639 zt(ktop) = ze(ktop)
2640
2641 ! -----------------------------------------------------------------------
2642 ! update capacity heat and latend heat coefficient
2643 ! -----------------------------------------------------------------------
2644
2645 do k = k0, kbot
2646 lhi(k) = li00 + dc_ice * tz(k)
2647 icpk(k) = lhi(k) / cvm(k)
2648 enddo
2649
2650 ! -----------------------------------------------------------------------
2651 ! melting of falling cloud ice into rain
2652 ! -----------------------------------------------------------------------
2653
2654 call check_column (ktop, kbot, qi, no_fall)
2655
2656 if (vi_fac < 1.e-5 .or. no_fall) then
2657 i1 = 0.
2658 else
2659
2660 do k = ktop + 1, kbot
2661 zt(k) = ze(k) - dt5 * (vti(k - 1) + vti(k))
2662 enddo
2663 zt(kbot + 1) = zs - dtm * vti(kbot)
2664
2665 do k = ktop, kbot
2666 if (zt(k + 1) >= zt(k)) zt(k + 1) = zt(k) - dz_min
2667 enddo
2668
2669 if (k0 < kbot) then
2670 do k = kbot - 1, k0, - 1
2671 if (qi(k) > qrmin) then
2672 do m = k + 1, kbot
2673 if (zt(k + 1) >= ze(m)) exit
2674 if (zt(k) < ze(m + 1) .and. tz(m) > tice) then
2675 dtime = min(1.0, (ze(m) - ze(m + 1)) / (max(vr_min, vti(k)) * tau_imlt))
2676 sink = min(qi(k) * dp(k) / dp(m), dtime * (tz(m) - tice) / icpk(m))
2677 tmp = min(sink, dim(ql_mlt, ql(m)))
2678 ql(m) = ql(m) + tmp
2679 qr(m) = qr(m) - tmp + sink
2680 tz(m) = tz(m) - sink * icpk(m)
2681 qi(k) = qi(k) - sink * dp(m) / dp(k)
2682 endif
2683 enddo
2684 endif
2685 enddo
2686 endif
2687
2688 if (do_sedi_w) then
2689 do k = ktop, kbot
2690 dm(k) = dp(k) * (1. + qv(k) + ql(k) + qr(k) + qi(k) + qs(k) + qg(k))
2691 enddo
2692 endif
2693
2694 if (use_ppm) then
2695 call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qi, i1, m1_sol, mono_prof)
2696 else
2697 call implicit_fall (dtm, ktop, kbot, ze, vti, dp, qi, i1, m1_sol)
2698 endif
2699
2700 if (do_sedi_w) then
2701 w1(ktop) = (dm(ktop) * w1(ktop) + m1_sol(ktop) * vti(ktop)) / (dm(ktop) - m1_sol(ktop))
2702 do k = ktop + 1, kbot
2703 w1(k) = (dm(k) * w1(k) - m1_sol(k - 1) * vti(k - 1) + m1_sol(k) * vti(k)) &
2704 / (dm(k) + m1_sol(k - 1) - m1_sol(k))
2705 enddo
2706 endif
2707
2708 endif
2709
2710 ! -----------------------------------------------------------------------
2711 ! melting of falling snow into rain
2712 ! -----------------------------------------------------------------------
2713
2714 r1 = 0.
2715
2716 call check_column (ktop, kbot, qs, no_fall)
2717
2718 if (no_fall) then
2719 s1 = 0.
2720 else
2721
2722 do k = ktop + 1, kbot
2723 zt(k) = ze(k) - dt5 * (vts(k - 1) + vts(k))
2724 enddo
2725 zt(kbot + 1) = zs - dtm * vts(kbot)
2726
2727 do k = ktop, kbot
2728 if (zt(k + 1) >= zt(k)) zt(k + 1) = zt(k) - dz_min
2729 enddo
2730
2731 if (k0 < kbot) then
2732 do k = kbot - 1, k0, - 1
2733 if (qs(k) > qrmin) then
2734 do m = k + 1, kbot
2735 if (zt(k + 1) >= ze(m)) exit
2736 dtime = min(dtm, (ze(m) - ze(m + 1)) / (vr_min + vts(k)))
2737 if (zt(k) < ze(m + 1) .and. tz(m) > tice) then
2738 dtime = min(1.0, dtime / tau_smlt)
2739 sink = min(qs(k) * dp(k) / dp(m), dtime * (tz(m) - tice) / icpk(m))
2740 tz(m) = tz(m) - sink * icpk(m)
2741 qs(k) = qs(k) - sink * dp(m) / dp(k)
2742 if (zt(k) < zs) then
2743 r1 = r1 + sink * dp(m) ! precip as rain
2744 else
2745 ! qr source here will fall next time step (therefore, can evap)
2746 qr(m) = qr(m) + sink
2747 endif
2748 endif
2749 if (qs(k) < qrmin) exit
2750 enddo
2751 endif
2752 enddo
2753 endif
2754
2755 if (do_sedi_w) then
2756 do k = ktop, kbot
2757 dm(k) = dp(k) * (1. + qv(k) + ql(k) + qr(k) + qi(k) + qs(k) + qg(k))
2758 enddo
2759 endif
2760
2761 if (use_ppm) then
2762 call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qs, s1, m1, mono_prof)
2763 else
2764 call implicit_fall (dtm, ktop, kbot, ze, vts, dp, qs, s1, m1)
2765 endif
2766
2767 do k = ktop, kbot
2768 m1_sol(k) = m1_sol(k) + m1(k)
2769 enddo
2770
2771 if (do_sedi_w) then
2772 w1(ktop) = (dm(ktop) * w1(ktop) + m1(ktop) * vts(ktop)) / (dm(ktop) - m1(ktop))
2773 do k = ktop + 1, kbot
2774 w1(k) = (dm(k) * w1(k) - m1(k - 1) * vts(k - 1) + m1(k) * vts(k)) &
2775 / (dm(k) + m1(k - 1) - m1(k))
2776 enddo
2777 endif
2778
2779 endif
2780
2781 ! ----------------------------------------------
2782 ! melting of falling graupel into rain
2783 ! ----------------------------------------------
2784
2785 call check_column (ktop, kbot, qg, no_fall)
2786
2787 if (no_fall) then
2788 g1 = 0.
2789 else
2790
2791 do k = ktop + 1, kbot
2792 zt(k) = ze(k) - dt5 * (vtg(k - 1) + vtg(k))
2793 enddo
2794 zt(kbot + 1) = zs - dtm * vtg(kbot)
2795
2796 do k = ktop, kbot
2797 if (zt(k + 1) >= zt(k)) zt(k + 1) = zt(k) - dz_min
2798 enddo
2799
2800 if (k0 < kbot) then
2801 do k = kbot - 1, k0, - 1
2802 if (qg(k) > qrmin) then
2803 do m = k + 1, kbot
2804 if (zt(k + 1) >= ze(m)) exit
2805 dtime = min(dtm, (ze(m) - ze(m + 1)) / vtg(k))
2806 if (zt(k) < ze(m + 1) .and. tz(m) > tice) then
2807 dtime = min(1., dtime / tau_g2r)
2808 sink = min(qg(k) * dp(k) / dp(m), dtime * (tz(m) - tice) / icpk(m))
2809 tz(m) = tz(m) - sink * icpk(m)
2810 qg(k) = qg(k) - sink * dp(m) / dp(k)
2811 if (zt(k) < zs) then
2812 r1 = r1 + sink * dp(m)
2813 else
2814 qr(m) = qr(m) + sink
2815 endif
2816 endif
2817 if (qg(k) < qrmin) exit
2818 enddo
2819 endif
2820 enddo
2821 endif
2822
2823 if (do_sedi_w) then
2824 do k = ktop, kbot
2825 dm(k) = dp(k) * (1. + qv(k) + ql(k) + qr(k) + qi(k) + qs(k) + qg(k))
2826 enddo
2827 endif
2828
2829 if (use_ppm) then
2830 call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qg, g1, m1, mono_prof)
2831 else
2832 call implicit_fall (dtm, ktop, kbot, ze, vtg, dp, qg, g1, m1)
2833 endif
2834
2835 do k = ktop, kbot
2836 m1_sol(k) = m1_sol(k) + m1(k)
2837 enddo
2838
2839 if (do_sedi_w) then
2840 w1(ktop) = (dm(ktop) * w1(ktop) + m1(ktop) * vtg(ktop)) / (dm(ktop) - m1(ktop))
2841 do k = ktop + 1, kbot
2842 w1(k) = (dm(k) * w1(k) - m1(k - 1) * vtg(k - 1) + m1(k) * vtg(k)) &
2843 / (dm(k) + m1(k - 1) - m1(k))
2844 enddo
2845 endif
2846
2847 endif
2848
2849end subroutine terminal_fall
2850
2851! =======================================================================
2855subroutine check_column (ktop, kbot, q, no_fall)
2856
2857 implicit none
2858
2859 integer, intent (in) :: ktop, kbot
2860
2861 real, intent (in) :: q (ktop:kbot)
2862
2863 logical, intent (out) :: no_fall
2864
2865 integer :: k
2866
2867 no_fall = .true.
2868
2869 do k = ktop, kbot
2870 if (q(k) > qrmin) then
2871 no_fall = .false.
2872 exit
2873 endif
2874 enddo
2875
2876end subroutine check_column
2877
2878! =======================================================================
2883subroutine implicit_fall (dt, ktop, kbot, ze, vt, dp, q, precip, m1)
2884
2885 implicit none
2886
2887 integer, intent (in) :: ktop, kbot
2888
2889 real, intent (in) :: dt
2890
2891 real, intent (in), dimension (ktop:kbot + 1) :: ze
2892
2893 real, intent (in), dimension (ktop:kbot) :: vt, dp
2894
2895 real, intent (inout), dimension (ktop:kbot) :: q
2896
2897 real, intent (out), dimension (ktop:kbot) :: m1
2898
2899 real, intent (out) :: precip
2900
2901 real, dimension (ktop:kbot) :: dz, qm, dd
2902
2903 integer :: k
2904
2905 do k = ktop, kbot
2906 dz(k) = ze(k) - ze(k + 1)
2907 dd(k) = dt * vt(k)
2908 q(k) = q(k) * dp(k)
2909 enddo
2910
2911 ! -----------------------------------------------------------------------
2912 ! sedimentation: non - vectorizable loop
2913 ! -----------------------------------------------------------------------
2914
2915 qm(ktop) = q(ktop) / (dz(ktop) + dd(ktop))
2916 do k = ktop + 1, kbot
2917 qm(k) = (q(k) + dd(k - 1) * qm(k - 1)) / (dz(k) + dd(k))
2918 enddo
2919
2920 ! -----------------------------------------------------------------------
2921 ! qm is density at this stage
2922 ! -----------------------------------------------------------------------
2923
2924 do k = ktop, kbot
2925 qm(k) = qm(k) * dz(k)
2926 enddo
2927
2928 ! -----------------------------------------------------------------------
2929 ! output mass fluxes: non - vectorizable loop
2930 ! -----------------------------------------------------------------------
2931
2932 m1(ktop) = q(ktop) - qm(ktop)
2933 do k = ktop + 1, kbot
2934 m1(k) = m1(k - 1) + q(k) - qm(k)
2935 enddo
2936 precip = m1(kbot)
2937
2938 ! -----------------------------------------------------------------------
2939 ! update:
2940 ! -----------------------------------------------------------------------
2941
2942 do k = ktop, kbot
2943 q(k) = qm(k) / dp(k)
2944 enddo
2945
2946end subroutine implicit_fall
2947
2948! =======================================================================
2952subroutine lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, q, precip, m1, mono)
2953
2954 implicit none
2955
2956 integer, intent (in) :: ktop, kbot
2957
2958 real, intent (in) :: zs
2959
2960 logical, intent (in) :: mono
2961
2962 real, intent (in), dimension (ktop:kbot + 1) :: ze, zt
2963
2964 real, intent (in), dimension (ktop:kbot) :: dp
2965
2966 ! m1: flux
2967 real, intent (inout), dimension (ktop:kbot) :: q, m1
2968
2969 real, intent (out) :: precip
2970
2971 real, dimension (ktop:kbot) :: qm, dz
2972
2973 real :: a4 (4, ktop:kbot)
2974
2975 real :: pl, pr, delz, esl
2976
2977 integer :: k, k0, n, m
2978
2979 real, parameter :: r3 = 1. / 3., r23 = 2. / 3.
2980
2981 ! -----------------------------------------------------------------------
2982 ! density:
2983 ! -----------------------------------------------------------------------
2984
2985 do k = ktop, kbot
2986 dz(k) = zt(k) - zt(k + 1) ! note: dz is positive
2987 q(k) = q(k) * dp(k)
2988 a4(1, k) = q(k) / dz(k)
2989 qm(k) = 0.
2990 enddo
2991
2992 ! -----------------------------------------------------------------------
2993 ! construct vertical profile with zt as coordinate
2994 ! -----------------------------------------------------------------------
2995
2996 call cs_profile (a4(1, ktop), dz(ktop), kbot - ktop + 1, mono)
2997
2998 k0 = ktop
2999 do k = ktop, kbot
3000 do n = k0, kbot
3001 if (ze(k) <= zt(n) .and. ze(k) >= zt(n + 1)) then
3002 pl = (zt(n) - ze(k)) / dz(n)
3003 if (zt(n + 1) <= ze(k + 1)) then
3004 ! entire new grid is within the original grid
3005 pr = (zt(n) - ze(k + 1)) / dz(n)
3006 qm(k) = a4(2, n) + 0.5 * (a4(4, n) + a4(3, n) - a4(2, n)) * (pr + pl) - &
3007 a4(4, n) * r3 * (pr * (pr + pl) + pl ** 2)
3008 qm(k) = qm(k) * (ze(k) - ze(k + 1))
3009 k0 = n
3010 goto 555
3011 else
3012 qm(k) = (ze(k) - zt(n + 1)) * (a4(2, n) + 0.5 * (a4(4, n) + &
3013 a4(3, n) - a4(2, n)) * (1. + pl) - a4(4, n) * (r3 * (1. + pl * (1. + pl))))
3014 if (n < kbot) then
3015 do m = n + 1, kbot
3016 ! locate the bottom edge: ze (k + 1)
3017 if (ze(k + 1) < zt(m + 1)) then
3018 qm(k) = qm(k) + q(m)
3019 else
3020 delz = zt(m) - ze(k + 1)
3021 esl = delz / dz(m)
3022 qm(k) = qm(k) + delz * (a4(2, m) + 0.5 * esl * &
3023 (a4(3, m) - a4(2, m) + a4(4, m) * (1. - r23 * esl)))
3024 k0 = m
3025 goto 555
3026 endif
3027 enddo
3028 endif
3029 goto 555
3030 endif
3031 endif
3032 enddo
3033 555 continue
3034 enddo
3035
3036 m1(ktop) = q(ktop) - qm(ktop)
3037 do k = ktop + 1, kbot
3038 m1(k) = m1(k - 1) + q(k) - qm(k)
3039 enddo
3040 precip = m1(kbot)
3041
3042 ! convert back to * dry * mixing ratio:
3043 ! dp must be dry air_mass (because moist air mass will be changed due to terminal fall) .
3044
3045 do k = ktop, kbot
3046 q(k) = qm(k) / dp(k)
3047 enddo
3048
3049end subroutine lagrangian_fall_ppm
3050
3052subroutine cs_profile (a4, del, km, do_mono)
3053
3054 implicit none
3055
3056 integer, intent (in) :: km ! vertical dimension
3057
3058 real, intent (in) :: del (km)
3059
3060 logical, intent (in) :: do_mono
3061
3062 real, intent (inout) :: a4 (4, km)
3063
3064 real, parameter :: qp_min = 1.e-6
3065
3066 real :: gam (km)
3067 real :: q (km + 1)
3068 real :: d4, bet, a_bot, grat, pmp, lac
3069 real :: pmp_1, lac_1, pmp_2, lac_2
3070 real :: da1, da2, a6da
3071
3072 integer :: k
3073
3074 logical extm (km)
3075
3076 grat = del(2) / del(1) ! grid ratio
3077 bet = grat * (grat + 0.5)
3078 q(1) = (2. * grat * (grat + 1.) * a4(1, 1) + a4(1, 2)) / bet
3079 gam(1) = (1. + grat * (grat + 1.5)) / bet
3080
3081 do k = 2, km
3082 d4 = del(k - 1) / del(k)
3083 bet = 2. + 2. * d4 - gam(k - 1)
3084 q(k) = (3. * (a4(1, k - 1) + d4 * a4(1, k)) - q(k - 1)) / bet
3085 gam(k) = d4 / bet
3086 enddo
3087
3088 a_bot = 1. + d4 * (d4 + 1.5)
3089 q(km + 1) = (2. * d4 * (d4 + 1.) * a4(1, km) + a4(1, km - 1) - a_bot * q(km)) &
3090 / (d4 * (d4 + 0.5) - a_bot * gam(km))
3091
3092 do k = km, 1, - 1
3093 q(k) = q(k) - gam(k) * q(k + 1)
3094 enddo
3095
3096 ! -----------------------------------------------------------------------
3097 ! apply constraints
3098 ! -----------------------------------------------------------------------
3099
3100 do k = 2, km
3101 gam(k) = a4(1, k) - a4(1, k - 1)
3102 enddo
3103
3104 ! -----------------------------------------------------------------------
3105 ! apply large - scale constraints to all fields if not local max / min
3106 ! -----------------------------------------------------------------------
3107
3108 ! -----------------------------------------------------------------------
3109 ! top:
3110 ! -----------------------------------------------------------------------
3111
3112 q(1) = max(q(1), 0.)
3113 q(2) = min(q(2), max(a4(1, 1), a4(1, 2)))
3114 q(2) = max(q(2), min(a4(1, 1), a4(1, 2)), 0.)
3115
3116 ! -----------------------------------------------------------------------
3117 ! interior:
3118 ! -----------------------------------------------------------------------
3119
3120 do k = 3, km - 1
3121 if (gam(k - 1) * gam(k + 1) > 0.) then
3122 q(k) = min(q(k), max(a4(1, k - 1), a4(1, k)))
3123 q(k) = max(q(k), min(a4(1, k - 1), a4(1, k)))
3124 else
3125 if (gam(k - 1) > 0.) then
3126 ! there exists a local max
3127 q(k) = max(q(k), min(a4(1, k - 1), a4(1, k)))
3128 else
3129 ! there exists a local min
3130 q(k) = min(q(k), max(a4(1, k - 1), a4(1, k)))
3131 q(k) = max(q(k), 0.0)
3132 endif
3133 endif
3134 enddo
3135
3136 ! -----------------------------------------------------------------------
3137 ! bottom :
3138 ! -----------------------------------------------------------------------
3139
3140 q(km) = min(q(km), max(a4(1, km - 1), a4(1, km)))
3141 q(km) = max(q(km), min(a4(1, km - 1), a4(1, km)), 0.)
3142 ! q (km + 1) = max (q (km + 1), 0.)
3143
3144 ! -----------------------------------------------------------------------
3145 ! f (s) = al + s * [ (ar - al) + a6 * (1 - s) ] (0 <= s <= 1)
3146 ! -----------------------------------------------------------------------
3147
3148 do k = 1, km - 1
3149 a4(2, k) = q(k)
3150 a4(3, k) = q(k + 1)
3151 enddo
3152
3153 do k = 2, km - 1
3154 if (gam(k) * gam(k + 1) > 0.0) then
3155 extm(k) = .false.
3156 else
3157 extm(k) = .true.
3158 endif
3159 enddo
3160
3161 if (do_mono) then
3162 do k = 3, km - 2
3163 if (extm(k)) then
3164 ! positive definite constraint only if true local extrema
3165 if (a4(1, k) < qp_min .or. extm(k - 1) .or. extm(k + 1)) then
3166 a4(2, k) = a4(1, k)
3167 a4(3, k) = a4(1, k)
3168 endif
3169 else
3170 a4(4, k) = 6. * a4(1, k) - 3. * (a4(2, k) + a4(3, k))
3171 if (abs(a4(4, k)) > abs(a4(2, k) - a4(3, k))) then
3172 ! check within the smooth region if subgrid profile is non - monotonic
3173 pmp_1 = a4(1, k) - 2.0 * gam(k + 1)
3174 lac_1 = pmp_1 + 1.5 * gam(k + 2)
3175 a4(2, k) = min(max(a4(2, k), min(a4(1, k), pmp_1, lac_1)), &
3176 max(a4(1, k), pmp_1, lac_1))
3177 pmp_2 = a4(1, k) + 2.0 * gam(k)
3178 lac_2 = pmp_2 - 1.5 * gam(k - 1)
3179 a4(3, k) = min(max(a4(3, k), min(a4(1, k), pmp_2, lac_2)), &
3180 max(a4(1, k), pmp_2, lac_2))
3181 endif
3182 endif
3183 enddo
3184 else
3185 do k = 3, km - 2
3186 if (extm(k)) then
3187 if (a4(1, k) < qp_min .or. extm(k - 1) .or. extm(k + 1)) then
3188 a4(2, k) = a4(1, k)
3189 a4(3, k) = a4(1, k)
3190 endif
3191 endif
3192 enddo
3193 endif
3194
3195 do k = 1, km - 1
3196 a4(4, k) = 6. * a4(1, k) - 3. * (a4(2, k) + a4(3, k))
3197 enddo
3198
3199 k = km - 1
3200 if (extm(k)) then
3201 a4(2, k) = a4(1, k)
3202 a4(3, k) = a4(1, k)
3203 a4(4, k) = 0.
3204 else
3205 da1 = a4(3, k) - a4(2, k)
3206 da2 = da1 ** 2
3207 a6da = a4(4, k) * da1
3208 if (a6da < - da2) then
3209 a4(4, k) = 3. * (a4(2, k) - a4(1, k))
3210 a4(3, k) = a4(2, k) - a4(4, k)
3211 elseif (a6da > da2) then
3212 a4(4, k) = 3. * (a4(3, k) - a4(1, k))
3213 a4(2, k) = a4(3, k) - a4(4, k)
3214 endif
3215 endif
3216
3217 call cs_limiters (km - 1, a4)
3218
3219 ! -----------------------------------------------------------------------
3220 ! bottom layer:
3221 ! -----------------------------------------------------------------------
3222
3223 a4(2, km) = a4(1, km)
3224 a4(3, km) = a4(1, km)
3225 a4(4, km) = 0.
3226
3227end subroutine cs_profile
3228
3231subroutine cs_limiters (km, a4)
3232
3233 implicit none
3234
3235 integer, intent (in) :: km
3236
3237 real, intent (inout) :: a4 (4, km) ! ppm array
3238
3239 real, parameter :: r12 = 1. / 12.
3240
3241 integer :: k
3242
3243 ! -----------------------------------------------------------------------
3244 ! positive definite constraint
3245 ! -----------------------------------------------------------------------
3246
3247 do k = 1, km
3248 if (abs(a4(3, k) - a4(2, k)) < - a4(4, k)) then
3249 if ((a4(1, k) + 0.25 * (a4(3, k) - a4(2, k)) ** 2 / a4(4, k) + a4(4, k) * r12) < 0.) then
3250 if (a4(1, k) < a4(3, k) .and. a4(1, k) < a4(2, k)) then
3251 a4(3, k) = a4(1, k)
3252 a4(2, k) = a4(1, k)
3253 a4(4, k) = 0.
3254 elseif (a4(3, k) > a4(2, k)) then
3255 a4(4, k) = 3. * (a4(2, k) - a4(1, k))
3256 a4(3, k) = a4(2, k) - a4(4, k)
3257 else
3258 a4(4, k) = 3. * (a4(3, k) - a4(1, k))
3259 a4(2, k) = a4(3, k) - a4(4, k)
3260 endif
3261 endif
3262 endif
3263 enddo
3264
3265end subroutine cs_limiters
3266
3267! =======================================================================
3270subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg)
3271
3272 implicit none
3273
3274 integer, intent (in) :: ktop, kbot
3275
3276 real, intent (in), dimension (ktop:kbot) :: den, qs, qi, qg, ql, tk
3277 real, intent (out), dimension (ktop:kbot) :: vts, vti, vtg
3278
3279 ! fall velocity constants:
3280
3281 real, parameter :: thi = 1.0e-8
3282 real, parameter :: thg = 1.0e-8
3283 real, parameter :: ths = 1.0e-8
3284
3285 ! coefficient for the parameterization of mass weighted fall velocity
3286 ! as a function of temperature and IWC.
3287 ! Table 1 in Deng and Mace (2008) \cite deng_and_mace_2008.
3288 real, parameter :: aa = - 4.14122e-5
3289 real, parameter :: bb = - 0.00538922
3290 real, parameter :: cc = - 0.0516344
3291 real, parameter :: dd = 0.00216078
3292 real, parameter :: ee = 1.9714
3293
3294 ! marshall - palmer constants
3295
3296 real, parameter :: vcons = 6.6280504
3297 real, parameter :: vcong = 87.2382675
3298 real, parameter :: vconh = vcong * sqrt(rhoh / rhog)
3299 real, parameter :: norms = 942477796.076938
3300 real, parameter :: normg = 5026548245.74367
3301 real, parameter :: normh = pi * rhoh * rnzh
3302
3303 real, dimension (ktop:kbot) :: qden, tc, rhof
3304
3305 real :: vi0
3306
3307 integer :: k
3308
3309 ! -----------------------------------------------------------------------
3310 ! marshall - palmer formula
3311 ! -----------------------------------------------------------------------
3312
3313 ! -----------------------------------------------------------------------
3314 ! try the local air density -- for global model; the true value could be
3315 ! much smaller than sfcrho over high mountains
3316 ! -----------------------------------------------------------------------
3317
3318 do k = ktop, kbot
3319 rhof(k) = sqrt(min(10., sfcrho / den(k)))
3320 enddo
3321
3322 ! -----------------------------------------------------------------------
3325 ! -----------------------------------------------------------------------
3326
3327 if (const_vi) then
3328 vti(:) = vi_fac
3329 else
3330 ! -----------------------------------------------------------------------
3331 ! use deng and mace (2008, grl), which gives smaller fall speed than hd90 formula
3332 ! -----------------------------------------------------------------------
3333 vi0 = 0.01 * vi_fac
3334 do k = ktop, kbot
3335 if (qi(k) < thi) then ! this is needed as the fall - speed maybe problematic for small qi
3336 vti(k) = vf_min
3337 else
3338 tc(k) = tk(k) - tice
3339 vti(k) = (3. + log10(qi(k) * den(k))) * (tc(k) * (aa * tc(k) + bb) + cc) + dd * tc(k) + ee
3340 vti(k) = vi0 * exp(log_10 * vti(k)) * 0.9
3341 vti(k) = min(vi_max, max(vf_min, vti(k)))
3342 endif
3343 enddo
3344 endif
3345
3346 ! -----------------------------------------------------------------------
3348 ! -----------------------------------------------------------------------
3349
3350 if (const_vs) then
3351 vts(:) = vs_fac ! 1. ifs_2016
3352 else
3353 do k = ktop, kbot
3354 if (qs(k) < ths) then
3355 vts(k) = vf_min
3356 else
3357 vts(k) = vs_fac * vcons * rhof(k) * exp(0.0625 * log(qs(k) * den(k) / norms))
3358 vts(k) = min(vs_max, max(vf_min, vts(k)))
3359 endif
3360 enddo
3361 endif
3362
3363 ! -----------------------------------------------------------------------
3365 ! -----------------------------------------------------------------------
3366 if (const_vg) then
3367 vtg(:) = vg_fac ! 2.
3368 else
3369 if (do_hail) then
3370 do k = ktop, kbot
3371 if (qg(k) < thg) then
3372 vtg(k) = vf_min
3373 else
3374 vtg(k) = vg_fac * vconh * rhof(k) * sqrt(sqrt(sqrt(qg(k) * den(k) / normh)))
3375 vtg(k) = min(vg_max, max(vf_min, vtg(k)))
3376 endif
3377 enddo
3378 else
3379 do k = ktop, kbot
3380 if (qg(k) < thg) then
3381 vtg(k) = vf_min
3382 else
3383 vtg(k) = vg_fac * vcong * rhof(k) * sqrt(sqrt(sqrt(qg(k) * den(k) / normg)))
3384 vtg(k) = min(vg_max, max(vf_min, vtg(k)))
3385 endif
3386 enddo
3387 endif
3388 endif
3389
3390end subroutine fall_speed
3391
3392! =======================================================================
3395subroutine setupm
3396
3397 implicit none
3398
3399 real :: gcon, cd, scm3, pisq, act (8)
3400 real :: vdifu, tcond
3401 real :: visk
3402 real :: ch2o, hltf
3403 real :: hlts, hltc, ri50
3404
3405 real, parameter :: gam263 = 1.456943, gam275 = 1.608355, gam290 = 1.827363, &
3406 gam325 = 2.54925, gam350 = 3.323363, gam380 = 4.694155, &
3407 gam425 = 8.285063, gam450 = 11.631769, gam480 = 17.837789, &
3408 gam625 = 184.860962, gam680 = 496.604067
3409
3410 ! intercept parameters
3411
3412! real, parameter :: rnzr = 8.0e6 ! lin83
3413! real, parameter :: rnzs = 3.0e6 ! lin83
3414! real, parameter :: rnzg = 4.0e6 ! rh84
3415
3416 ! density parameters
3417
3418! real, parameter :: rhos = 0.1e3 !< lin83 (snow density; 1 / 10 of water)
3419! real, parameter :: rhog = 0.4e3 !< rh84 (graupel density)
3420 real, parameter :: acc (3) = (/ 5.0, 2.0, 0.5 /)
3421
3422 real den_rc
3423
3424 integer :: i, k
3425
3426 pie = 4. * atan(1.0)
3427
3428 ! s. klein's formular (eq 16) from am2
3429
3430 fac_rc = (4. / 3.) * pie * rhor * rthresh ** 3
3431
3432 if (prog_ccn) then
3433 ! if (master) write (*, *) 'prog_ccn option is .t.'
3434 else
3435 den_rc = fac_rc * ccn_o * 1.e6
3436 ! if (master) write (*, *) 'mp: for ccn_o = ', ccn_o, 'ql_rc = ', den_rc
3437 den_rc = fac_rc * ccn_l * 1.e6
3438 ! if (master) write (*, *) 'mp: for ccn_l = ', ccn_l, 'ql_rc = ', den_rc
3439 endif
3440
3441 vdifu = 2.11e-5
3442 tcond = 2.36e-2
3443
3444 visk = 1.259e-5
3445 hlts = 2.8336e6
3446 hltc = 2.5e6
3447 hltf = 3.336e5
3448
3449 ch2o = 4.1855e3
3450 ri50 = 1.e-4
3451
3452 pisq = pie * pie
3453 scm3 = (visk / vdifu) ** (1. / 3.)
3454
3455 cracs = pisq * rnzr * rnzs * rhos
3456 csacr = pisq * rnzr * rnzs * rhor
3457 if (do_hail) then
3458 cgacr = pisq * rnzr * rnzh * rhor
3459 cgacs = pisq * rnzh * rnzs * rhos
3460 else
3461 cgacr = pisq * rnzr * rnzg * rhor
3462 cgacs = pisq * rnzg * rnzs * rhos
3463 endif
3464 cgacs = cgacs * c_pgacs
3465
3466 ! act: 1 - 2:racs (s - r) ; 3 - 4:sacr (r - s) ;
3467 ! 5 - 6:gacr (r - g) ; 7 - 8:gacs (s - g)
3468
3469 act(1) = pie * rnzs * rhos
3470 act(2) = pie * rnzr * rhor
3471 if (do_hail) then
3472 act(6) = pie * rnzh * rhoh
3473 else
3474 act(6) = pie * rnzg * rhog
3475 endif
3476 act(3) = act(2)
3477 act(4) = act(1)
3478 act(5) = act(2)
3479 act(7) = act(1)
3480 act(8) = act(6)
3481
3482 do i = 1, 3
3483 do k = 1, 4
3484 acco(i, k) = acc(i) / (act(2 * k - 1) ** ((7 - i) * 0.25) * act(2 * k) ** (i * 0.25))
3485 enddo
3486 enddo
3487
3488 gcon = 40.74 * sqrt(sfcrho) ! 44.628
3489
3490 csacw = pie * rnzs * clin * gam325 / (4. * act(1) ** 0.8125)
3491 ! decreasing csacw to reduce cloud water --- > snow
3492
3493 craci = pie * rnzr * alin * gam380 / (4. * act(2) ** 0.95)
3494 csaci = csacw * c_psaci
3495
3496 if (do_hail) then
3497 cgacw = pie * rnzh * gam350 * gcon / (4. * act(6) ** 0.875)
3498 else
3499 cgacw = pie * rnzg * gam350 * gcon / (4. * act(6) ** 0.875)
3500 endif
3501 ! cgaci = cgacw * 0.1
3502
3503 ! sjl, may 28, 2012
3504 cgaci = cgacw * 0.05
3505 ! sjl, may 28, 2012
3506
3507 cracw = craci ! cracw = 3.27206196043822
3508 cracw = c_cracw * cracw
3509
3510 ! subl and revp: five constants for three separate processes
3511
3512 cssub(1) = 2. * pie * vdifu * tcond * rvgas * rnzs
3513 if (do_hail) then
3514 cgsub(1) = 2. * pie * vdifu * tcond * rvgas * rnzh
3515 else
3516 cgsub(1) = 2. * pie * vdifu * tcond * rvgas * rnzg
3517 endif
3518 crevp(1) = 2. * pie * vdifu * tcond * rvgas * rnzr
3519 cssub(2) = 0.78 / sqrt(act(1))
3520 cgsub(2) = 0.78 / sqrt(act(6))
3521 crevp(2) = 0.78 / sqrt(act(2))
3522 cssub(3) = 0.31 * scm3 * gam263 * sqrt(clin / visk) / act(1) ** 0.65625
3523 cgsub(3) = 0.31 * scm3 * gam275 * sqrt(gcon / visk) / act(6) ** 0.6875
3524 crevp(3) = 0.31 * scm3 * gam290 * sqrt(alin / visk) / act(2) ** 0.725
3525 cssub(4) = tcond * rvgas
3526 cssub(5) = hlts ** 2 * vdifu
3527 cgsub(4) = cssub(4)
3528 crevp(4) = cssub(4)
3529 cgsub(5) = cssub(5)
3530 crevp(5) = hltc ** 2 * vdifu
3531
3532 cgfr(1) = 20.e2 * pisq * rnzr * rhor / act(2) ** 1.75
3533 cgfr(2) = 0.66
3534
3535 ! smlt: five constants (lin et al. 1983)
3536
3537 csmlt(1) = 2. * pie * tcond * rnzs / hltf
3538 csmlt(2) = 2. * pie * vdifu * rnzs * hltc / hltf
3539 csmlt(3) = cssub(2)
3540 csmlt(4) = cssub(3)
3541 csmlt(5) = ch2o / hltf
3542
3543 ! gmlt: five constants
3544
3545 if (do_hail) then
3546 cgmlt(1) = 2. * pie * tcond * rnzh / hltf
3547 cgmlt(2) = 2. * pie * vdifu * rnzh * hltc / hltf
3548 else
3549 cgmlt(1) = 2. * pie * tcond * rnzg / hltf
3550 cgmlt(2) = 2. * pie * vdifu * rnzg * hltc / hltf
3551 endif
3552 cgmlt(3) = cgsub(2)
3553 cgmlt(4) = cgsub(3)
3554 cgmlt(5) = ch2o / hltf
3555
3556 es0 = 6.107799961e2 ! ~6.1 mb
3557 ces0 = eps * es0
3558
3559end subroutine setupm
3560
3561! =======================================================================
3562! initialization of gfdl cloud microphysics
3566subroutine gfdl_cloud_microphys_mod_init (me, master, nlunit, input_nml_file, logunit, &
3567 fn_nml, errmsg, errflg)
3568
3569 implicit none
3570
3571 integer, intent (in) :: me
3572 integer, intent (in) :: master
3573 integer, intent (in) :: nlunit
3574 integer, intent (in) :: logunit
3575
3576 character (len = 64), intent (in) :: fn_nml
3577 character (len = *), intent (in) :: input_nml_file(:)
3578 character(len=*), intent(out) :: errmsg
3579 integer, intent(out) :: errflg
3580
3581 integer :: ios
3582 logical :: exists
3583
3584 ! integer, intent (in) :: id, jd, kd
3585 ! integer, intent (in) :: axes (4)
3586 ! type (time_type), intent (in) :: time
3587
3588 ! integer :: unit, io, ierr, k, logunit
3589 ! logical :: flag
3590 ! real :: tmp, q1, q2
3591
3592 ! master = (mpp_pe () .eq.mpp_root_pe ())
3593
3594 ! Initialize CCPP error-handling
3595 errflg = 0
3596 errmsg = ''
3597
3598#ifdef INTERNAL_FILE_NML
3599 read (input_nml_file, nml = gfdl_cloud_microphysics_nml)
3600#else
3601 inquire (file = trim(fn_nml), exist = exists)
3602 if (.not. exists) then
3603 write (6, *) 'gfdl - mp :: namelist file: ', trim(fn_nml), ' does not exist'
3604 errflg = 1
3605 errmsg = 'ERROR(gfdl_cloud_microphys_mod_init): namelist file '//trim(fn_nml)//' does not exist'
3606 return
3607 else
3608 open (unit = nlunit, file = fn_nml, action = 'read' , status = 'old', iostat = ios)
3609 endif
3610 rewind(nlunit)
3611 read (nlunit, nml = gfdl_cloud_microphysics_nml)
3612 close (nlunit)
3613#endif
3614
3615 ! write version number and namelist to log file
3616 if (me == master) then
3617 write (logunit, *) " ================================================================== "
3618 write (logunit, *) "gfdl_cloud_microphys_mod"
3619 write (logunit, nml = gfdl_cloud_microphysics_nml)
3620 endif
3621
3622 if (do_setup) then
3623 call setup_con
3624 call setupm
3625 do_setup = .false.
3626 endif
3627
3628 log_10 = log(10.)
3629
3630 tice0 = tice - 0.01
3631 t_wfr = tice - 40.0 ! supercooled water can exist down to - 48 c, which is the "absolute"
3632
3633 ! if (master) write (logunit, nml = gfdl_cloud_microphys_nml)
3634 !
3635 ! id_vtr = register_diag_field (mod_name, 'vt_r', axes (1:3), time, &
3636 ! 'rain fall speed', 'm / s', missing_value = missing_value)
3637 ! id_vts = register_diag_field (mod_name, 'vt_s', axes (1:3), time, &
3638 ! 'snow fall speed', 'm / s', missing_value = missing_value)
3639 ! id_vtg = register_diag_field (mod_name, 'vt_g', axes (1:3), time, &
3640 ! 'graupel fall speed', 'm / s', missing_value = missing_value)
3641 ! id_vti = register_diag_field (mod_name, 'vt_i', axes (1:3), time, &
3642 ! 'ice fall speed', 'm / s', missing_value = missing_value)
3643
3644 ! id_droplets = register_diag_field (mod_name, 'droplets', axes (1:3), time, &
3645 ! 'droplet number concentration', '# / m3', missing_value = missing_value)
3646 ! id_rh = register_diag_field (mod_name, 'rh_lin', axes (1:2), time, &
3647 ! 'relative humidity', 'n / a', missing_value = missing_value)
3648
3649 ! id_rain = register_diag_field (mod_name, 'rain_lin', axes (1:2), time, &
3650 ! 'rain_lin', 'mm / day', missing_value = missing_value)
3651 ! id_snow = register_diag_field (mod_name, 'snow_lin', axes (1:2), time, &
3652 ! 'snow_lin', 'mm / day', missing_value = missing_value)
3653 ! id_graupel = register_diag_field (mod_name, 'graupel_lin', axes (1:2), time, &
3654 ! 'graupel_lin', 'mm / day', missing_value = missing_value)
3655 ! id_ice = register_diag_field (mod_name, 'ice_lin', axes (1:2), time, &
3656 ! 'ice_lin', 'mm / day', missing_value = missing_value)
3657 ! id_prec = register_diag_field (mod_name, 'prec_lin', axes (1:2), time, &
3658 ! 'prec_lin', 'mm / day', missing_value = missing_value)
3659
3660 ! if (master) write (*, *) 'prec_lin diagnostics initialized.', id_prec
3661
3662 ! id_cond = register_diag_field (mod_name, 'cond_lin', axes (1:2), time, &
3663 ! 'total condensate', 'kg / m ** 2', missing_value = missing_value)
3664 ! id_var = register_diag_field (mod_name, 'var_lin', axes (1:2), time, &
3665 ! 'subgrid variance', 'n / a', missing_value = missing_value)
3666
3667 ! call qsmith_init
3668
3669 ! testing the water vapor tables
3670
3671 ! if (mp_debug .and. master) then
3672 ! write (*, *) 'testing water vapor tables in gfdl_cloud_microphys'
3673 ! tmp = tice - 90.
3674 ! do k = 1, 25
3675 ! q1 = wqsat_moist (tmp, 0., 1.e5)
3676 ! q2 = qs1d_m (tmp, 0., 1.e5)
3677 ! write (*, *) nint (tmp - tice), q1, q2, 'dq = ', q1 - q2
3678 ! tmp = tmp + 5.
3679 ! enddo
3680 ! endif
3681
3682 ! if (master) write (*, *) 'gfdl_cloud_micrphys diagnostics initialized.'
3683
3684 module_is_initialized = .true.
3685
3686!+---+-----------------------------------------------------------------+
3687!..Set these variables needed for computing radar reflectivity. These
3688!.. get used within radar_init to create other variables used in the
3689!.. radar module.
3690
3691 xam_r = pi*rhor/6.
3692 xbm_r = 3.
3693 xmu_r = 0.
3694 xam_s = pi*rhos/6.
3695 xbm_s = 3.
3696 xmu_s = 0.
3697 xam_g = pi*rhog/6.
3698 xbm_g = 3.
3699 xmu_g = 0.
3700
3701 call radar_init
3702
3703end subroutine gfdl_cloud_microphys_mod_init
3704
3705! =======================================================================
3706! end of gfdl cloud microphysics
3711
3712 implicit none
3713
3714 deallocate (table)
3715 deallocate (table2)
3716 deallocate (table3)
3717 deallocate (tablew)
3718 deallocate (des)
3719 deallocate (des2)
3720 deallocate (des3)
3721 deallocate (desw)
3722
3723 tables_are_initialized = .false.
3724
3725end subroutine gfdl_cloud_microphys_mod_end
3726
3727! =======================================================================
3728! qsmith table initialization
3731subroutine setup_con
3732
3733 implicit none
3734
3735 ! master = (mpp_pe () .eq.mpp_root_pe ())
3736
3737 rgrav = 1. / grav
3738
3739 if (.not. qsmith_tables_initialized) call qsmith_init
3740
3741 qsmith_tables_initialized = .true.
3742
3743end subroutine setup_con
3744
3745! =======================================================================
3749! =======================================================================
3750
3751real function acr3d (v1, v2, q1, q2, c, cac, rho)
3752
3753 implicit none
3754
3755 real, intent (in) :: v1, v2, c, rho
3756 real, intent (in) :: q1, q2 ! mixing ratio!!!
3757 real, intent (in) :: cac (3)
3758
3759 real :: t1, s1, s2
3760
3761 ! integer :: k
3762 !
3763 ! real :: a
3764 !
3765 ! a = 0.0
3766 ! do k = 1, 3
3767 ! a = a + cac (k) * ((q1 * rho) ** ((7 - k) * 0.25) * (q2 * rho) ** (k * 0.25))
3768 ! enddo
3769 ! acr3d = c * abs (v1 - v2) * a / rho
3770
3771 ! optimized
3772
3773 t1 = sqrt(q1 * rho)
3774 s1 = sqrt(q2 * rho)
3775 s2 = sqrt(s1) ! s1 = s2 ** 2
3776 acr3d = c * abs(v1 - v2) * q1 * s2 * (cac(1) * t1 + cac(2) * sqrt(t1) * s2 + cac(3) * s1)
3777
3778end function acr3d
3779
3780! =======================================================================
3784! =======================================================================
3785
3786real function smlt (tc, dqs, qsrho, psacw, psacr, c, rho, rhofac)
3787
3788 implicit none
3789
3790 real, intent (in) :: tc, dqs, qsrho, psacw, psacr, c (5), rho, rhofac
3791
3792 smlt = (c(1) * tc / rho - c(2) * dqs) * (c(3) * sqrt(qsrho) + &
3793 c(4) * qsrho ** 0.65625 * sqrt(rhofac)) + c(5) * tc * (psacw + psacr)
3794
3795end function smlt
3796
3797! =======================================================================
3800! =======================================================================
3801
3802real function gmlt (tc, dqs, qgrho, pgacw, pgacr, c, rho)
3803
3804 implicit none
3805
3806 real, intent (in) :: tc, dqs, qgrho, pgacw, pgacr, c (5), rho
3807
3808 gmlt = (c(1) * tc / rho - c(2) * dqs) * (c(3) * sqrt(qgrho) + &
3809 c(4) * qgrho ** 0.6875 / rho ** 0.25) + c(5) * tc * (pgacw + pgacr)
3810
3811end function gmlt
3812
3813! =======================================================================
3814! initialization
3815! prepare saturation water vapor pressure tables
3816! =======================================================================
3823! =======================================================================
3824subroutine qsmith_init
3825
3826 implicit none
3827
3828 integer, parameter :: length = 2621
3829
3830 integer :: i
3831
3832 if (.not. tables_are_initialized) then
3833
3834 ! master = (mpp_pe () .eq. mpp_root_pe ())
3835 ! if (master) print *, ' gfdl mp: initializing qs tables'
3836
3837 ! debug code
3838 ! print *, mpp_pe (), allocated (table), allocated (table2), &
3839 ! allocated (table3), allocated (tablew), allocated (des), &
3840 ! allocated (des2), allocated (des3), allocated (desw)
3841 ! end debug code
3842
3843 ! generate es table (dt = 0.1 deg. c)
3844
3845 allocate (table(length))
3846 allocate (table2(length))
3847 allocate (table3(length))
3848 allocate (tablew(length))
3849 allocate (des(length))
3850 allocate (des2(length))
3851 allocate (des3(length))
3852 allocate (desw(length))
3853
3854 call qs_table (length)
3855 call qs_table2 (length)
3856 call qs_table3 (length)
3857 call qs_tablew (length)
3858
3859 do i = 1, length - 1
3860 des(i) = max(0., table(i + 1) - table(i))
3861 des2(i) = max(0., table2(i + 1) - table2(i))
3862 des3(i) = max(0., table3(i + 1) - table3(i))
3863 desw(i) = max(0., tablew(i + 1) - tablew(i))
3864 enddo
3865 des(length) = des(length - 1)
3866 des2(length) = des2(length - 1)
3867 des3(length) = des3(length - 1)
3868 desw(length) = desw(length - 1)
3869
3870 tables_are_initialized = .true.
3871
3872 endif
3873
3874end subroutine qsmith_init
3875
3876! =======================================================================
3877! compute the saturated specific humidity for table ii
3881real function wqs1 (ta, den)
3882
3883 implicit none
3884
3887
3888 real, intent (in) :: ta, den
3889
3890 real :: es, ap1, tmin
3891
3892 integer :: it
3893
3894 tmin = table_ice - 160.
3895 ap1 = 10. * dim(ta, tmin) + 1.
3896 ap1 = min(2621., ap1)
3897 it = ap1
3898 es = tablew(it) + (ap1 - it) * desw(it)
3899 wqs1 = es / (rvgas * ta * den)
3900
3901end function wqs1
3902
3903! =======================================================================
3904! compute the gradient of saturated specific humidity for table ii
3908! =======================================================================
3909
3910real function wqs2 (ta, den, dqdt)
3911
3912 implicit none
3913
3914 ! pure water phase; universal dry / moist formular using air density
3915 ! input "den" can be either dry or moist air density
3916
3917 real, intent (in) :: ta, den
3918
3919 real, intent (out) :: dqdt
3920
3921 real :: es, ap1, tmin
3922
3923 integer :: it
3924
3925 tmin = table_ice - 160.
3926
3927 if (.not. tables_are_initialized) call qsmith_init
3928
3929 ap1 = 10. * dim(ta, tmin) + 1.
3930 ap1 = min(2621., ap1)
3931 it = ap1
3932 es = tablew(it) + (ap1 - it) * desw(it)
3933 wqs2 = es / (rvgas * ta * den)
3934 it = ap1 - 0.5
3935 ! finite diff, del_t = 0.1:
3936 dqdt = 10. * (desw(it) + (ap1 - it) * (desw(it + 1) - desw(it))) / (rvgas * ta * den)
3937
3938end function wqs2
3939
3940! =======================================================================
3941! compute wet buld temperature
3944! =======================================================================
3945
3946real function wet_bulb (q, t, den)
3947
3948 implicit none
3949
3950 real, intent (in) :: t, q, den
3951
3952 real :: qs, tp, dqdt
3953
3954 wet_bulb = t
3955 qs = wqs2(wet_bulb, den, dqdt)
3956 tp = 0.5 * (qs - q) / (1. + lcp * dqdt) * lcp
3957 wet_bulb = wet_bulb - tp
3958
3959 ! tp is negative if super - saturated
3960 if (tp > 0.01) then
3961 qs = wqs2(wet_bulb, den, dqdt)
3962 tp = (qs - q) / (1. + lcp * dqdt) * lcp
3963 wet_bulb = wet_bulb - tp
3964 endif
3965
3966end function wet_bulb
3967
3968! =======================================================================
3971! =======================================================================
3972
3973real function iqs1 (ta, den)
3974
3975 implicit none
3976
3979
3980 real, intent (in) :: ta, den
3981
3982 real :: es, ap1, tmin
3983
3984 integer :: it
3985
3986 tmin = table_ice - 160.
3987 ap1 = 10. * dim(ta, tmin) + 1.
3988 ap1 = min(2621., ap1)
3989 it = ap1
3990 es = table2(it) + (ap1 - it) * des2(it)
3991 iqs1 = es / (rvgas * ta * den)
3992
3993end function iqs1
3994
3995! =======================================================================
3998! =======================================================================
3999
4000real function iqs2 (ta, den, dqdt)
4001
4002 implicit none
4003
4004 ! water - ice phase; universal dry / moist formular using air density
4005 ! input "den" can be either dry or moist air density
4006
4007 real, intent (in) :: ta, den
4008
4009 real, intent (out) :: dqdt
4010
4011 real :: es, ap1, tmin
4012
4013 integer :: it
4014
4015 tmin = table_ice - 160.
4016 ap1 = 10. * dim(ta, tmin) + 1.
4017 ap1 = min(2621., ap1)
4018 it = ap1
4019 es = table2(it) + (ap1 - it) * des2(it)
4020 iqs2 = es / (rvgas * ta * den)
4021 it = ap1 - 0.5
4022 dqdt = 10. * (des2(it) + (ap1 - it) * (des2(it + 1) - des2(it))) / (rvgas * ta * den)
4023
4024end function iqs2
4025
4026! =======================================================================
4029! =======================================================================
4030
4031real function qs1d_moist (ta, qv, pa, dqdt)
4032
4033 implicit none
4034
4035 real, intent (in) :: ta, pa, qv
4036
4037 real, intent (out) :: dqdt
4038
4039 real :: es, ap1, tmin, eps10
4040
4041 integer :: it
4042
4043 tmin = table_ice - 160.
4044 eps10 = 10. * eps
4045 ap1 = 10. * dim(ta, tmin) + 1.
4046 ap1 = min(2621., ap1)
4047 it = ap1
4048 es = table2(it) + (ap1 - it) * des2(it)
4049 qs1d_moist = eps * es * (1. + zvir * qv) / pa
4050 it = ap1 - 0.5
4051 dqdt = eps10 * (des2(it) + (ap1 - it) * (des2(it + 1) - des2(it))) * (1. + zvir * qv) / pa
4052
4053end function qs1d_moist
4054
4055! =======================================================================
4056! compute the gradient of saturated specific humidity for table ii
4059! =======================================================================
4060
4061real function wqsat2_moist (ta, qv, pa, dqdt)
4062
4063 implicit none
4064
4065 real, intent (in) :: ta, pa, qv
4066
4067 real, intent (out) :: dqdt
4068
4069 real :: es, ap1, tmin, eps10
4070
4071 integer :: it
4072
4073 tmin = table_ice - 160.
4074 eps10 = 10. * eps
4075 ap1 = 10. * dim(ta, tmin) + 1.
4076 ap1 = min(2621., ap1)
4077 it = ap1
4078 es = tablew(it) + (ap1 - it) * desw(it)
4079 wqsat2_moist = eps * es * (1. + zvir * qv) / pa
4080 it = ap1 - 0.5
4081 dqdt = eps10 * (desw(it) + (ap1 - it) * (desw(it + 1) - desw(it))) * (1. + zvir * qv) / pa
4082
4083end function wqsat2_moist
4084
4085! =======================================================================
4086! compute the saturated specific humidity for table ii
4089! =======================================================================
4090
4091real function wqsat_moist (ta, qv, pa)
4092
4093 implicit none
4094
4095 real, intent (in) :: ta, pa, qv
4096
4097 real :: es, ap1, tmin
4098
4099 integer :: it
4100
4101 tmin = table_ice - 160.
4102 ap1 = 10. * dim(ta, tmin) + 1.
4103 ap1 = min(2621., ap1)
4104 it = ap1
4105 es = tablew(it) + (ap1 - it) * desw(it)
4106 wqsat_moist = eps * es * (1. + zvir * qv) / pa
4107
4108end function wqsat_moist
4109
4110! =======================================================================
4113! =======================================================================
4114
4115real function qs1d_m (ta, qv, pa)
4116
4117 implicit none
4118
4119 real, intent (in) :: ta, pa, qv
4120
4121 real :: es, ap1, tmin
4122
4123 integer :: it
4124
4125 tmin = table_ice - 160.
4126 ap1 = 10. * dim(ta, tmin) + 1.
4127 ap1 = min(2621., ap1)
4128 it = ap1
4129 es = table2(it) + (ap1 - it) * des2(it)
4130 qs1d_m = eps * es * (1. + zvir * qv) / pa
4131
4132end function qs1d_m
4133
4134! =======================================================================
4137! =======================================================================
4138
4139real function d_sat (ta, den)
4140
4141 implicit none
4142
4143 real, intent (in) :: ta, den
4144
4145 real :: es_w, es_i, ap1, tmin
4146
4147 integer :: it
4148
4149 tmin = table_ice - 160.
4150 ap1 = 10. * dim(ta, tmin) + 1.
4151 ap1 = min(2621., ap1)
4152 it = ap1
4153 es_w = tablew(it) + (ap1 - it) * desw(it)
4154 es_i = table2(it) + (ap1 - it) * des2(it)
4155 d_sat = dim(es_w, es_i) / (rvgas * ta * den) ! take positive difference
4156
4157end function d_sat
4158
4159! =======================================================================
4162! =======================================================================
4163
4164real function esw_table (ta)
4165
4166 implicit none
4167
4168 real, intent (in) :: ta
4169
4170 real :: ap1, tmin
4171
4172 integer :: it
4173
4174 tmin = table_ice - 160.
4175 ap1 = 10. * dim(ta, tmin) + 1.
4176 ap1 = min(2621., ap1)
4177 it = ap1
4178 esw_table = tablew(it) + (ap1 - it) * desw(it)
4179
4180end function esw_table
4181
4182! =======================================================================
4185! =======================================================================
4186
4187real function es2_table (ta)
4188
4189 implicit none
4190
4191 real, intent (in) :: ta
4192
4193 real :: ap1, tmin
4194
4195 integer :: it
4196
4197 tmin = table_ice - 160.
4198 ap1 = 10. * dim(ta, tmin) + 1.
4199 ap1 = min(2621., ap1)
4200 it = ap1
4201 es2_table = table2(it) + (ap1 - it) * des2(it)
4202
4203end function es2_table
4204
4205! =======================================================================
4209subroutine esw_table1d (ta, es, n)
4210
4211 implicit none
4212
4213 integer, intent (in) :: n
4214
4215 real, intent (in) :: ta (n)
4216
4217 real, intent (out) :: es (n)
4218
4219 real :: ap1, tmin
4220
4221 integer :: i, it
4222
4223 tmin = table_ice - 160.
4224
4225 do i = 1, n
4226 ap1 = 10. * dim(ta(i), tmin) + 1.
4227 ap1 = min(2621., ap1)
4228 it = ap1
4229 es(i) = tablew(it) + (ap1 - it) * desw(it)
4230 enddo
4231
4232end subroutine esw_table1d
4233
4234! =======================================================================
4238subroutine es2_table1d (ta, es, n)
4239
4240 implicit none
4241
4242 integer, intent (in) :: n
4243
4244 real, intent (in) :: ta (n)
4245
4246 real, intent (out) :: es (n)
4247
4248 real :: ap1, tmin
4249
4250 integer :: i, it
4251
4252 tmin = table_ice - 160.
4253
4254 do i = 1, n
4255 ap1 = 10. * dim(ta(i), tmin) + 1.
4256 ap1 = min(2621., ap1)
4257 it = ap1
4258 es(i) = table2(it) + (ap1 - it) * des2(it)
4259 enddo
4260
4261end subroutine es2_table1d
4262
4263! =======================================================================
4267subroutine es3_table1d (ta, es, n)
4268
4269 implicit none
4270
4271 integer, intent (in) :: n
4272
4273 real, intent (in) :: ta (n)
4274
4275 real, intent (out) :: es (n)
4276
4277 real :: ap1, tmin
4278
4279 integer :: i, it
4280
4281 tmin = table_ice - 160.
4282
4283 do i = 1, n
4284 ap1 = 10. * dim(ta(i), tmin) + 1.
4285 ap1 = min(2621., ap1)
4286 it = ap1
4287 es(i) = table3(it) + (ap1 - it) * des3(it)
4288 enddo
4289
4290end subroutine es3_table1d
4291
4292! =======================================================================
4295! 1 - phase table
4296subroutine qs_tablew (n)
4297
4298 implicit none
4299
4300 integer, intent (in) :: n
4301
4302 real :: delt = 0.1
4303 real :: tmin, tem, fac0, fac1, fac2
4304
4305 integer :: i
4306
4307 tmin = table_ice - 160.
4308
4309 ! -----------------------------------------------------------------------
4310 ! compute es over water
4311 ! -----------------------------------------------------------------------
4312
4313 do i = 1, n
4314 tem = tmin + delt * real(i - 1)
4315 fac0 = (tem - t_ice) / (tem * t_ice)
4316 fac1 = fac0 * lv0
4317 fac2 = (dc_vap * log(tem / t_ice) + fac1) / rvgas
4318 tablew(i) = e00 * exp(fac2)
4319 enddo
4320
4321end subroutine qs_tablew
4322
4323! =======================================================================
4326! 2 - phase table
4327subroutine qs_table2 (n)
4328
4329 implicit none
4330
4331 integer, intent (in) :: n
4332
4333 real :: delt = 0.1
4334 real :: tmin, tem0, tem1, fac0, fac1, fac2
4335
4336 integer :: i, i0, i1
4337
4338 tmin = table_ice - 160.
4339
4340 do i = 1, n
4341 tem0 = tmin + delt * real(i - 1)
4342 fac0 = (tem0 - t_ice) / (tem0 * t_ice)
4343 if (i <= 1600) then
4344 ! -----------------------------------------------------------------------
4345 ! compute es over ice between - 160 deg c and 0 deg c.
4346 ! -----------------------------------------------------------------------
4347 fac1 = fac0 * li2
4348 fac2 = (d2ice * log(tem0 / t_ice) + fac1) / rvgas
4349 else
4350 ! -----------------------------------------------------------------------
4351 ! compute es over water between 0 deg c and 102 deg c.
4352 ! -----------------------------------------------------------------------
4353 fac1 = fac0 * lv0
4354 fac2 = (dc_vap * log(tem0 / t_ice) + fac1) / rvgas
4355 endif
4356 table2(i) = e00 * exp(fac2)
4357 enddo
4358
4359 ! -----------------------------------------------------------------------
4360 ! smoother around 0 deg c
4361 ! -----------------------------------------------------------------------
4362
4363 i0 = 1600
4364 i1 = 1601
4365 tem0 = 0.25 * (table2(i0 - 1) + 2. * table(i0) + table2(i0 + 1))
4366 tem1 = 0.25 * (table2(i1 - 1) + 2. * table(i1) + table2(i1 + 1))
4367 table2(i0) = tem0
4368 table2(i1) = tem1
4369
4370end subroutine qs_table2
4371
4372! =======================================================================
4375! 2 - phase table with " - 2 c" as the transition point
4376subroutine qs_table3 (n)
4377
4378 implicit none
4379
4380 integer, intent (in) :: n
4381
4382 real :: delt = 0.1
4383 real :: esbasw, tbasw, esbasi, tmin, tem, aa, b, c, d, e
4384 real :: tem0, tem1
4385
4386 integer :: i, i0, i1
4387
4388 esbasw = 1013246.0
4389 tbasw = table_ice + 100.
4390 esbasi = 6107.1
4391 tmin = table_ice - 160.
4392
4393 do i = 1, n
4394 tem = tmin + delt * real(i - 1)
4395 ! if (i <= 1600) then
4396 if (i <= 1580) then ! change to - 2 c
4397 ! -----------------------------------------------------------------------
4398 ! compute es over ice between - 160 deg c and 0 deg c.
4399 ! see smithsonian meteorological tables page 350.
4400 ! -----------------------------------------------------------------------
4401 aa = - 9.09718 * (table_ice / tem - 1.)
4402 b = - 3.56654 * alog10(table_ice / tem)
4403 c = 0.876793 * (1. - tem / table_ice)
4404 e = alog10(esbasi)
4405 table3(i) = 0.1 * 10 ** (aa + b + c + e)
4406 else
4407 ! -----------------------------------------------------------------------
4408 ! compute es over water between - 2 deg c and 102 deg c.
4409 ! see smithsonian meteorological tables page 350.
4410 ! -----------------------------------------------------------------------
4411 aa = - 7.90298 * (tbasw / tem - 1.)
4412 b = 5.02808 * alog10(tbasw / tem)
4413 c = - 1.3816e-7 * (10 ** ((1. - tem / tbasw) * 11.344) - 1.)
4414 d = 8.1328e-3 * (10 ** ((tbasw / tem - 1.) * (- 3.49149)) - 1.)
4415 e = alog10(esbasw)
4416 table3(i) = 0.1 * 10 ** (aa + b + c + d + e)
4417 endif
4418 enddo
4419
4420 ! -----------------------------------------------------------------------
4421 ! smoother around - 2 deg c
4422 ! -----------------------------------------------------------------------
4423
4424 i0 = 1580
4425 i1 = 1581
4426 tem0 = 0.25 * (table3(i0 - 1) + 2. * table(i0) + table3(i0 + 1))
4427 tem1 = 0.25 * (table3(i1 - 1) + 2. * table(i1) + table3(i1 + 1))
4428 table3(i0) = tem0
4429 table3(i1) = tem1
4430
4431end subroutine qs_table3
4432
4433! =======================================================================
4434! compute the saturated specific humidity for table
4435! note: this routine is based on "moist" mixing ratio
4439real function qs_blend (t, p, q)
4440
4441 implicit none
4442
4443 real, intent (in) :: t, p, q
4444
4445 real :: es, ap1, tmin
4446
4447 integer :: it
4448
4449 tmin = table_ice - 160.
4450 ap1 = 10. * dim(t, tmin) + 1.
4451 ap1 = min(2621., ap1)
4452 it = ap1
4453 es = table(it) + (ap1 - it) * des(it)
4454 qs_blend = eps * es * (1. + zvir * q) / p
4455
4456end function qs_blend
4457
4458! =======================================================================
4461! 3 - phase table
4462subroutine qs_table (n)
4463
4464 implicit none
4465
4466 integer, intent (in) :: n
4467
4468 real :: delt = 0.1
4469 real :: tmin, tem, esh20
4470 real :: wice, wh2o, fac0, fac1, fac2
4471 real :: esupc (200)
4472
4473 integer :: i
4474
4475 tmin = table_ice - 160.
4476
4477 ! -----------------------------------------------------------------------
4478 ! compute es over ice between - 160 deg c and 0 deg c.
4479 ! -----------------------------------------------------------------------
4480
4481 do i = 1, 1600
4482 tem = tmin + delt * real(i - 1)
4483 fac0 = (tem - t_ice) / (tem * t_ice)
4484 fac1 = fac0 * li2
4485 fac2 = (d2ice * log(tem / t_ice) + fac1) / rvgas
4486 table(i) = e00 * exp(fac2)
4487 enddo
4488
4489 ! -----------------------------------------------------------------------
4490 ! compute es over water between - 20 deg c and 102 deg c.
4491 ! -----------------------------------------------------------------------
4492
4493 do i = 1, 1221
4494 tem = 253.16 + delt * real(i - 1)
4495 fac0 = (tem - t_ice) / (tem * t_ice)
4496 fac1 = fac0 * lv0
4497 fac2 = (dc_vap * log(tem / t_ice) + fac1) / rvgas
4498 esh20 = e00 * exp(fac2)
4499 if (i <= 200) then
4500 esupc(i) = esh20
4501 else
4502 table(i + 1400) = esh20
4503 endif
4504 enddo
4505
4506 ! -----------------------------------------------------------------------
4507 ! derive blended es over ice and supercooled water between - 20 deg c and 0 deg c
4508 ! -----------------------------------------------------------------------
4509
4510 do i = 1, 200
4511 tem = 253.16 + delt * real(i - 1)
4512 wice = 0.05 * (table_ice - tem)
4513 wh2o = 0.05 * (tem - 253.16)
4514 table(i + 1400) = wice * table(i + 1400) + wh2o * esupc(i)
4515 enddo
4516
4517end subroutine qs_table
4518
4519! =======================================================================
4520! compute the saturated specific humidity and the gradient of saturated specific humidity
4521! input t in deg k, p in pa; p = rho rdry tv, moist pressure
4525!@details It als oincludes the option for computing des/dT.
4526subroutine qsmith (im, km, ks, t, p, q, qs, dqdt)
4527
4528 implicit none
4529
4530 integer, intent (in) :: im, km, ks
4531
4532 real, intent (in), dimension (im, km) :: t, p, q
4533
4534 real, intent (out), dimension (im, km) :: qs
4535
4536 real, intent (out), dimension (im, km), optional :: dqdt
4537
4538 real :: eps10, ap1, tmin
4539
4540 real, dimension (im, km) :: es
4541
4542 integer :: i, k, it
4543
4544 tmin = table_ice - 160.
4545 eps10 = 10. * eps
4546
4547 if (.not. tables_are_initialized) then
4548 call qsmith_init
4549 endif
4550
4551 do k = ks, km
4552 do i = 1, im
4553 ap1 = 10. * dim(t(i, k), tmin) + 1.
4554 ap1 = min(2621., ap1)
4555 it = ap1
4556 es(i, k) = table(it) + (ap1 - it) * des(it)
4557 qs(i, k) = eps * es(i, k) * (1. + zvir * q(i, k)) / p(i, k)
4558 enddo
4559 enddo
4560
4561 if (present (dqdt)) then
4562 do k = ks, km
4563 do i = 1, im
4564 ap1 = 10. * dim(t(i, k), tmin) + 1.
4565 ap1 = min(2621., ap1) - 0.5
4566 it = ap1
4567 dqdt(i, k) = eps10 * (des(it) + (ap1 - it) * (des(it + 1) - des(it))) * (1. + zvir * q(i, k)) / p(i, k)
4568 enddo
4569 enddo
4570 endif
4571
4572end subroutine qsmith
4573
4574! =======================================================================
4578subroutine neg_adj (ktop, kbot, pt, dp, qv, ql, qr, qi, qs, qg)
4579
4580 implicit none
4581
4582 integer, intent (in) :: ktop, kbot
4583
4584 real, intent (in), dimension (ktop:kbot) :: dp
4585
4586 real, intent (inout), dimension (ktop:kbot) :: pt, qv, ql, qr, qi, qs, qg
4587
4588 real, dimension (ktop:kbot) :: lcpk, icpk
4589
4590 real :: dq, cvm
4591
4592 integer :: k
4593
4594 ! -----------------------------------------------------------------------
4595 ! define heat capacity and latent heat coefficient
4596 ! -----------------------------------------------------------------------
4597
4598 do k = ktop, kbot
4599 cvm = c_air + qv(k) * c_vap + (qr(k) + ql(k)) * c_liq + (qi(k) + qs(k) + qg(k)) * c_ice
4600 lcpk(k) = (lv00 + d0_vap * pt(k)) / cvm
4601 icpk(k) = (li00 + dc_ice * pt(k)) / cvm
4602 enddo
4603
4604 do k = ktop, kbot
4605
4606 ! -----------------------------------------------------------------------
4607 ! ice phase:
4608 ! -----------------------------------------------------------------------
4609
4610 ! if cloud ice < 0, borrow from snow
4611 if (qi(k) < 0.) then
4612 qs(k) = qs(k) + qi(k)
4613 qi(k) = 0.
4614 endif
4615 ! if snow < 0, borrow from graupel
4616 if (qs(k) < 0.) then
4617 qg(k) = qg(k) + qs(k)
4618 qs(k) = 0.
4619 endif
4620 ! if graupel < 0, borrow from rain
4621 if (qg(k) < 0.) then
4622 qr(k) = qr(k) + qg(k)
4623 pt(k) = pt(k) - qg(k) * icpk(k) ! heating
4624 qg(k) = 0.
4625 endif
4626
4627 ! -----------------------------------------------------------------------
4628 ! liquid phase:
4629 ! -----------------------------------------------------------------------
4630
4631 ! if rain < 0, borrow from cloud water
4632 if (qr(k) < 0.) then
4633 ql(k) = ql(k) + qr(k)
4634 qr(k) = 0.
4635 endif
4636 ! if cloud water < 0, borrow from water vapor
4637 if (ql(k) < 0.) then
4638 qv(k) = qv(k) + ql(k)
4639 pt(k) = pt(k) - ql(k) * lcpk(k) ! heating
4640 ql(k) = 0.
4641 endif
4642
4643 enddo
4644
4645 ! -----------------------------------------------------------------------
4646 ! fix water vapor; borrow from below
4647 ! -----------------------------------------------------------------------
4648
4649 do k = ktop, kbot - 1
4650 if (qv(k) < 0.) then
4651 qv(k + 1) = qv(k + 1) + qv(k) * dp(k) / dp(k + 1)
4652 qv(k) = 0.
4653 endif
4654 enddo
4655
4656 ! -----------------------------------------------------------------------
4657 ! bottom layer; borrow from above
4658 ! -----------------------------------------------------------------------
4659
4660 if (qv(kbot) < 0. .and. qv(kbot - 1) > 0.) then
4661 dq = min(- qv(kbot) * dp(kbot), qv(kbot - 1) * dp(kbot - 1))
4662 qv(kbot - 1) = qv(kbot - 1) - dq / dp(kbot - 1)
4663 qv(kbot) = qv(kbot) + dq / dp(kbot)
4664 endif
4665
4666end subroutine neg_adj
4667
4668! =======================================================================
4669! compute global sum
4671! =======================================================================
4672
4673!real function g_sum (p, ifirst, ilast, jfirst, jlast, area, mode)
4674!
4675! use mpp_mod, only: mpp_sum
4676!
4677! implicit none
4678!
4679! integer, intent (in) :: ifirst, ilast, jfirst, jlast
4680! integer, intent (in) :: mode ! if == 1 divided by area
4681!
4682! real, intent (in), dimension (ifirst:ilast, jfirst:jlast) :: p, area
4683!
4684! integer :: i, j
4685!
4686! real :: gsum
4687!
4688! if (global_area < 0.) then
4689! global_area = 0.
4690! do j = jfirst, jlast
4691! do i = ifirst, ilast
4692! global_area = global_area + area (i, j)
4693! enddo
4694! enddo
4695! call mpp_sum (global_area)
4696! endif
4697!
4698! gsum = 0.
4699! do j = jfirst, jlast
4700! do i = ifirst, ilast
4701! gsum = gsum + p (i, j) * area (i, j)
4702! enddo
4703! enddo
4704! call mpp_sum (gsum)
4705!
4706! if (mode == 1) then
4707! g_sum = gsum / global_area
4708! else
4709! g_sum = gsum
4710! endif
4711!
4712!end function g_sum
4713
4714! ==========================================================================
4717subroutine interpolate_z (is, ie, js, je, km, zl, hgt, a3, a2)
4718
4719 implicit none
4720
4721 integer, intent (in) :: is, ie, js, je, km
4722
4723 real, intent (in), dimension (is:ie, js:je, km) :: a3
4724
4725 real, intent (in), dimension (is:ie, js:je, km + 1) :: hgt ! hgt (k) > hgt (k + 1)
4726
4727 real, intent (in) :: zl
4728
4729 real, intent (out), dimension (is:ie, js:je) :: a2
4730
4731 real, dimension (km) :: zm
4732
4733 integer :: i, j, k
4734
4735 !$omp parallel do default (none) shared (is, ie, js, je, km, hgt, zl, a2, a3) private (zm)
4736
4737 do j = js, je
4738 do i = is, ie
4739 do k = 1, km
4740 zm(k) = 0.5 * (hgt(i, j, k) + hgt(i, j, k + 1))
4741 enddo
4742 if (zl >= zm(1)) then
4743 a2(i, j) = a3(i, j, 1)
4744 elseif (zl <= zm(km)) then
4745 a2(i, j) = a3(i, j, km)
4746 else
4747 do k = 1, km - 1
4748 if (zl <= zm(k) .and. zl >= zm(k + 1)) then
4749 a2(i, j) = a3(i, j, k) + (a3(i, j, k + 1) - a3(i, j, k)) * (zm(k) - zl) / (zm(k) - zm(k + 1))
4750 exit
4751 endif
4752 enddo
4753 endif
4754 enddo
4755 enddo
4756
4757end subroutine interpolate_z
4758
4759! =======================================================================
4764! =======================================================================
4765subroutine cloud_diagnosis (is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms, qmg, t, &
4766 rew, rei, rer, res, reg)
4767
4768 implicit none
4769
4770 integer, intent (in) :: is, ie, ks, ke
4771 integer, intent (in), dimension (is:ie) :: lsm ! land sea mask, 0: ocean, 1: land, 2: sea ice
4772
4773 real, intent (in), dimension (is:ie, ks:ke) :: den, delp, t
4774 real, intent (in), dimension (is:ie, ks:ke) :: qmw, qmi, qmr, qms, qmg
4775
4776 real, intent (out), dimension (is:ie, ks:ke) :: rew, rei, rer, res, reg
4777
4778 real, dimension (is:ie, ks:ke) :: qcw, qci, qcr, qcs, qcg
4779
4780 integer :: i, k
4781
4782 real :: lambdar, lambdas, lambdag
4783 real :: dpg, rei_fac, mask, ccn, bw
4784 real, parameter :: rho_0 = 50.e-3
4785
4786 real :: rhow = 1.0e3, rhor = 1.0e3, rhos = 1.0e2, rhog = 4.0e2
4787 real :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6
4788 real :: alphar = 0.8, alphas = 0.25, alphag = 0.5
4789 real :: gammar = 17.837789, gammas = 8.2850630, gammag = 11.631769
4790 real :: qmin = 1.0e-12, beta = 1.22, qmin1 = 9.e-6
4791
4792 do k = ks, ke
4793 do i = is, ie
4794
4795 dpg = abs(delp(i, k)) / grav
4796 mask = min(max(real(lsm(i)), 0.0), 2.0)
4797
4798 ! -----------------------------------------------------------------------
4799 ! cloud water (Martin et al., 1994)
4800 ! -----------------------------------------------------------------------
4801
4802 ccn = 0.80 * (- 1.15e-3 * (ccn_o ** 2) + 0.963 * ccn_o + 5.30) * abs(mask - 1.0) + &
4803 0.67 * (- 2.10e-4 * (ccn_l ** 2) + 0.568 * ccn_l - 27.9) * (1.0 - abs(mask - 1.0))
4804
4805 if (qmw(i, k) .gt. qmin) then
4806 qcw(i, k) = dpg * qmw(i, k) * 1.0e3
4807 rew(i, k) = exp(1.0 / 3.0 * log((3.0 * den(i, k) * qmw(i, k)) / (4.0 * pi * rhow * ccn))) * 1.0e4
4808 rew(i, k) = max(rewmin, min(rewmax, rew(i, k)))
4809 else
4810 qcw(i, k) = 0.0
4811 rew(i, k) = rewmin
4812 endif
4813
4814 if (reiflag .eq. 1) then
4815
4816 ! -----------------------------------------------------------------------
4817 ! cloud ice (Heymsfield and Mcfarquhar, 1996)
4818 ! -----------------------------------------------------------------------
4819
4820 if (qmi(i, k) .gt. qmin1) then
4821 qci(i, k) = dpg * qmi(i, k) * 1.0e3
4822 rei_fac = log(1.0e3 * qmi(i, k) * den(i, k))
4823 if (t(i, k) - tice .lt. - 50) then
4824 rei(i, k) = beta / 9.917 * exp(0.109 * rei_fac) * 1.0e3
4825 elseif (t(i, k) - tice .lt. - 40) then
4826 rei(i, k) = beta / 9.337 * exp(0.080 * rei_fac) * 1.0e3
4827 elseif (t(i, k) - tice .lt. - 30) then
4828 rei(i, k) = beta / 9.208 * exp(0.055 * rei_fac) * 1.0e3
4829 else
4830 rei(i, k) = beta / 9.387 * exp(0.031 * rei_fac) * 1.0e3
4831 endif
4832 rei(i, k) = max(reimin, min(reimax, rei(i, k)))
4833 else
4834 qci(i, k) = 0.0
4835 rei(i, k) = reimin
4836 endif
4837
4838 endif
4839
4840 if (reiflag .eq. 2) then
4841
4842 ! -----------------------------------------------------------------------
4843 ! cloud ice (Wyser, 1998)
4844 ! -----------------------------------------------------------------------
4845
4846 if (qmi(i, k) .gt. qmin1) then
4847 qci(i, k) = dpg * qmi(i, k) * 1.0e3
4848 bw = - 2. + 1.e-3 * log10(den(i, k) * qmi(i, k) / rho_0) * max(0.0, tice - t(i, k)) ** 1.5
4849 rei(i, k) = 377.4 + bw * (203.3 + bw * (37.91 + 2.3696 * bw))
4850 rei(i, k) = max(reimin, min(reimax, rei(i, k)))
4851 else
4852 qci(i, k) = 0.0
4853 rei(i, k) = reimin
4854 endif
4855
4856 endif
4857
4858 ! -----------------------------------------------------------------------
4859 ! rain (Lin et al., 1983)
4860 ! -----------------------------------------------------------------------
4861
4862 if (qmr(i, k) .gt. qmin) then
4863 qcr(i, k) = dpg * qmr(i, k) * 1.0e3
4864 lambdar = exp(0.25 * log(pi * rhor * n0r / qmr(i, k) / den(i, k)))
4865 rer(i, k) = 0.5 * exp(log(gammar / 6) / alphar) / lambdar * 1.0e6
4866 rer(i, k) = max(rermin, min(rermax, rer(i, k)))
4867 else
4868 qcr(i, k) = 0.0
4869 rer(i, k) = rermin
4870 endif
4871
4872 ! -----------------------------------------------------------------------
4873 ! snow (Lin et al., 1983)
4874 ! -----------------------------------------------------------------------
4875
4876 if (qms(i, k) .gt. qmin1) then
4877 qcs(i, k) = dpg * qms(i, k) * 1.0e3
4878 lambdas = exp(0.25 * log(pi * rhos * n0s / qms(i, k) / den(i, k)))
4879 res(i, k) = 0.5 * exp(log(gammas / 6) / alphas) / lambdas * 1.0e6
4880 res(i, k) = max(resmin, min(resmax, res(i, k)))
4881 else
4882 qcs(i, k) = 0.0
4883 res(i, k) = resmin
4884 endif
4885
4886 ! -----------------------------------------------------------------------
4887 ! graupel (Lin et al., 1983)
4888 ! -----------------------------------------------------------------------
4889
4890 if (qmg(i, k) .gt. qmin) then
4891 qcg(i, k) = dpg * qmg(i, k) * 1.0e3
4892 lambdag = exp(0.25 * log(pi * rhog * n0g / qmg(i, k) / den(i, k)))
4893 reg(i, k) = 0.5 * exp(log(gammag / 6) / alphag) / lambdag * 1.0e6
4894 reg(i, k) = max(regmin, min(regmax, reg(i, k)))
4895 else
4896 qcg(i, k) = 0.0
4897 reg(i, k) = regmin
4898 endif
4899
4900 enddo
4901 enddo
4902
4903end subroutine cloud_diagnosis
4904
4905!+---+-----------------------------------------------------------------+
4908 subroutine refl10cm_gfdl (qv1d, qr1d, qs1d, qg1d, &
4909 t1d, p1d, dBZ, kts, kte, ii, jj, melti)
4910
4911 IMPLICIT NONE
4912
4913!..Sub arguments
4914 INTEGER, INTENT(IN):: kts, kte, ii,jj
4915 REAL, DIMENSION(kts:kte), INTENT(IN):: &
4916 qv1d, qr1d, qs1d, qg1d, t1d, p1d
4917 REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
4918
4919!..Local variables
4920 REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
4921 REAL, DIMENSION(kts:kte):: rr, rs, rg
4922! REAL:: temp_C
4923
4924 DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams, ilamg
4925 DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s, N0_g
4926 DOUBLE PRECISION:: lamr, lams, lamg
4927 LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
4928
4929 REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
4930 DOUBLE PRECISION:: fmelt_s, fmelt_g
4931
4932 INTEGER:: i, k, k_0, kbot, n
4933 LOGICAL, INTENT(IN):: melti
4934 DOUBLE PRECISION:: cback, x, eta, f_d
4935!+---+
4936
4937 do k = kts, kte
4938 dbz(k) = -35.0
4939 enddo
4940
4941!+---+-----------------------------------------------------------------+
4942!..Put column of data into local arrays.
4943!+---+-----------------------------------------------------------------+
4944 do k = kts, kte
4945 temp(k) = t1d(k)
4946! temp_C = min(-0.001, temp(K)-273.15)
4947 qv(k) = max(1.e-10, qv1d(k))
4948 pres(k) = p1d(k)
4949 rho(k) = 0.622*pres(k)/(rdgas*temp(k)*(qv(k)+0.622))
4950
4951 if (qr1d(k) .gt. 1.e-9) then
4952 rr(k) = qr1d(k)*rho(k)
4953 n0_r(k) = n0r
4954 lamr = (xam_r*xcrg(3)*n0_r(k)/rr(k))**(1./xcre(1))
4955 ilamr(k) = 1./lamr
4956 l_qr(k) = .true.
4957 else
4958 rr(k) = 1.e-12
4959 l_qr(k) = .false.
4960 endif
4961
4962 if (qs1d(k) .gt. 1.e-9) then
4963 rs(k) = qs1d(k)*rho(k)
4964 n0_s(k) = n0s
4965 lams = (xam_s*xcsg(3)*n0_s(k)/rs(k))**(1./xcse(1))
4966 ilams(k) = 1./lams
4967 l_qs(k) = .true.
4968 else
4969 rs(k) = 1.e-12
4970 l_qs(k) = .false.
4971 endif
4972
4973 if (qg1d(k) .gt. 1.e-9) then
4974 rg(k) = qg1d(k)*rho(k)
4975 n0_g(k) = n0g
4976 lamg = (xam_g*xcgg(3)*n0_g(k)/rg(k))**(1./xcge(1))
4977 ilamg(k) = 1./lamg
4978 l_qg(k) = .true.
4979 else
4980 rg(k) = 1.e-12
4981 l_qg(k) = .false.
4982 endif
4983 enddo
4984
4985!+---+-----------------------------------------------------------------+
4986!..Locate K-level of start of melting (k_0 is level above).
4987!+---+-----------------------------------------------------------------+
4988 k_0 = kts
4989 k_loop:do k = kte-1, kts, -1
4990 if ( melti .and. (temp(k).gt.273.15) .and. l_qr(k) &
4991 .and. (l_qs(k+1).or.l_qg(k+1)) ) then
4992 k_0 = max(k+1, k_0)
4993 EXIT k_loop
4994 endif
4995 enddo k_loop
4996!+---+-----------------------------------------------------------------+
4997!..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
4998!.. and non-water-coated snow and graupel when below freezing are
4999!.. simple. Integrations of m(D)*m(D)*N(D)*dD.
5000!+---+-----------------------------------------------------------------+
5001 do k = kts, kte
5002 ze_rain(k) = 1.e-22
5003 ze_snow(k) = 1.e-22
5004 ze_graupel(k) = 1.e-22
5005 if (l_qr(k)) ze_rain(k) = n0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
5006 if (l_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/pi)*(6.0/pi) &
5007 * (xam_s/900.0)*(xam_s/900.0) &
5008 * n0_s(k)*xcsg(4)*ilams(k)**xcse(4)
5009 if (l_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/pi)*(6.0/pi) &
5010 * (xam_g/900.0)*(xam_g/900.0) &
5011 * n0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
5012 enddo
5013
5014
5015!+---+-----------------------------------------------------------------+
5016!..Special case of melting ice (snow/graupel) particles. Assume the
5017!.. ice is surrounded by the liquid water. Fraction of meltwater is
5018!.. extremely simple based on amount found above the melting level.
5019!.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
5020!.. routines).
5021!+---+-----------------------------------------------------------------+
5022
5023 if (melti .and. k_0.ge.kts+1) then
5024 do k = k_0-1, kts, -1
5025
5026!..Reflectivity contributed by melting snow
5027 if (l_qs(k) .and. l_qs(k_0) ) then
5028 fmelt_s = max(0.005d0, min(1.0d0-rs(k)/rs(k_0), 0.99d0))
5029 eta = 0.d0
5030 lams = 1./ilams(k)
5031 do n = 1, nrbins
5032 x = xam_s * xxds(n)**xbm_s
5033 call rayleigh_soak_wetgraupel (x,dble(xocms),dble(xobms), &
5034 fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
5035 cback, mixingrulestring_s, matrixstring_s, &
5036 inclusionstring_s, hoststring_s, &
5037 hostmatrixstring_s, hostinclusionstring_s)
5038 f_d = n0_s(k)*xxds(n)**xmu_s * dexp(-lams*xxds(n))
5039 eta = eta + f_d * cback * simpson(n) * xdts(n)
5040 enddo
5041 ze_snow(k) = sngl(lamda4 / (pi5 * k_w) * eta)
5042 endif
5043
5044
5045!..Reflectivity contributed by melting graupel
5046
5047 if (l_qg(k) .and. l_qg(k_0) ) then
5048 fmelt_g = max(0.005d0, min(1.0d0-rg(k)/rg(k_0), 0.99d0))
5049 eta = 0.d0
5050 lamg = 1./ilamg(k)
5051 do n = 1, nrbins
5052 x = xam_g * xxdg(n)**xbm_g
5053 call rayleigh_soak_wetgraupel (x,dble(xocmg),dble(xobmg), &
5054 fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
5055 cback, mixingrulestring_g, matrixstring_g, &
5056 inclusionstring_g, hoststring_g, &
5057 hostmatrixstring_g, hostinclusionstring_g)
5058 f_d = n0_g(k)*xxdg(n)**xmu_g * dexp(-lamg*xxdg(n))
5059 eta = eta + f_d * cback * simpson(n) * xdtg(n)
5060 enddo
5061 ze_graupel(k) = sngl(lamda4 / (pi5 * k_w) * eta)
5062 endif
5063
5064 enddo
5065 endif
5066
5067 do k = kte, kts, -1
5068 dbz(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
5069 enddo
5070
5071
5072 end subroutine refl10cm_gfdl
5073!+---+-----------------------------------------------------------------+
5074
5075end module gfdl_cloud_microphys_mod
real function acr3d(v1, v2, q1, q2, c, cac, rho)
The function is an accretion function (Lin et al.(1983) lin_et_al_1983 )
real function smlt(tc, dqs, qsrho, psacw, psacr, c, rho, rhofac)
Melting of snow function (Lin et al.(1983) lin_et_al_1983) note: psacw and psacr must be calc before ...
subroutine esw_table1d(ta, es, n)
The subroutine 'esw_table1d' computes the saturated water vapor pressure for table ii.
subroutine interpolate_z(is, ie, js, je, km, zl, hgt, a3, a2)
quick local sum algorithm
subroutine qs_table(n)
saturation water vapor pressure table i
subroutine qs_table3(n)
saturation water vapor pressure table iv
subroutine refl10cm_gfdl(qv1d, qr1d, qs1d, qg1d, t1d, p1d, dbz, kts, kte, ii, jj, melti)
This subroutine calculates radar reflectivity.
subroutine implicit_fall(dt, ktop, kbot, ze, vt, dp, q, precip, m1)
The subroutine computes the time-implicit monotonic fall scheme.
subroutine qsmith(im, km, ks, t, p, q, qs, dqdt)
The function 'qsmith' computes the saturated specific humidity with a blend of water and ice dependin...
subroutine, public cloud_diagnosis(is, ie, ks, ke, den, delp, lsm, qmw, qmi, qmr, qms, qmg, t, rew, rei, rer, res, reg)
The subroutine 'cloud_diagnosis' diagnoses the radius of cloud species.
subroutine cs_limiters(km, a4)
This subroutine perform positive definite constraint.
subroutine linear_prof(km, q, dm, z_var, h_var)
Definition of vertical subgrid variability used for cloud ice and cloud water autoconversion.
subroutine qs_table2(n)
saturation water vapor pressure table iii
subroutine fall_speed(ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg)
The subroutine calculates vertical fall speed of snow/ice/graupel.
subroutine check_column(ktop, kbot, q, no_fall)
The subroutine 'check_column' checks if the water species is large enough to fall.
subroutine revap_racc(ktop, kbot, dt, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
This subroutine calculates evaporation of rain and accretion of rain.
subroutine icloud(ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, den, denfac, vts, vtg, vtr, qak, rh_adj, rh_rain, dts, h_var)
This subroutine includes cloud ice microphysics processes.
real function wqs1(ta, den)
The function 'wqs1' returns the saturation vapor pressure over pure liquid water for a given temperat...
subroutine revap_rac1(hydrostatic, is, ie, dt, tz, qv, ql, qr, qi, qs, qg, den, hvar)
This subroutine calculates rain evaporation.
subroutine terminal_fall(dtm, ktop, kbot, tz, qv, ql, qr, qg, qs, qi, dz, dp, den, vtg, vts, vti, r1, g1, s1, i1, m1_sol, w1)
The subroutine 'terminal_fall' computes terminal fall speed.
subroutine neg_adj(ktop, kbot, pt, dp, qv, ql, qr, qi, qs, qg)
The subroutine 'neg_adj' fixes negative water species.
subroutine mpdrv(hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, qg, qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes, rain, snow, graupel, ice, m2_rain, m2_sol, cond, area1, land, u_dt, v_dt, pt_dt, qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, w_var, vt_r, vt_s, vt_g, vt_i, qn2)
GFDL cloud microphysics, major program, and is based on Lin et al.(1983) lin_et_al_1983 and Rutledge ...
subroutine setupm
The subroutine sets up gfdl cloud microphysics parameters.
subroutine es3_table1d(ta, es, n)
The subroutine 'es3_table1d' computes the saturated water vapor pressure for table iv.
subroutine, public gfdl_cloud_microphys_mod_end()
The subroutine 'gfdl_cloud_microphys_init' terminates the GFDL cloud microphysics.
subroutine, public gfdl_cloud_microphys_mod_driver(iis, iie, jjs, jje, kks, kke, ktop, kbot, qv, ql, qr, qi, qs, qg, qa, qn, qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, pt_dt, pt, w, uin, vin, udt, vdt, dz, delp, area, dt_in, land, rain, snow, ice, graupel, hydrostatic, phys_hydrostatic, p, lradar, refl_10cm, reset, pfils, pflls)
This subroutine is the driver of the GFDL cloud microphysics.
subroutine lagrangian_fall_ppm(ktop, kbot, zs, ze, zt, dp, q, precip, m1, mono)
Lagrangian scheme.
subroutine qsmith_init
The subroutine 'qsmith_init' initializes lookup tables for saturation water vapor pressure for the fo...
subroutine, public gfdl_cloud_microphys_mod_init(me, master, nlunit, input_nml_file, logunit, fn_nml, errmsg, errflg)
The subroutine 'gfdl_cloud_microphys_init' initializes the GFDL cloud microphysics.
subroutine warm_rain(dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, den, denfac, ccn, c_praut, rh_rain, vtr, r1, m1_rain, w1, h_var)
This subroutine includes warm rain cloud microphysics.
subroutine qs_tablew(n)
saturation water vapor pressure table ii
subroutine sedi_heat(ktop, kbot, dm, m1, dz, tz, qv, ql, qr, qi, qs, qg, cw)
This subroutine calculates sedimentation of heat.
real function qs_blend(t, p, q)
The function 'qs_blend' computes the saturated specific humidity with a blend of water and ice depend...
subroutine setup_con
The subroutine 'setup_con' sets up constants and calls 'qsmith_init'.
subroutine cs_profile(a4, del, km, do_mono)
subroutine subgrid_z_proc(ktop, kbot, p1, den, denfac, dts, rh_adj, tz, qv, ql, qr, qi, qs, qg, qa, h_var, rh_rain)
This subroutine calculates temperature sentive high vertical resolution processes.
subroutine es2_table1d(ta, es, n)
The subroutine 'es3_table1d' computes the saturated water vapor pressure for table iii.
This module is more library code whereas the individual microphysics schemes contains specific detail...