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