SCALE-RM
scale_atmos_phy_mp_tomita08.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
11 !-------------------------------------------------------------------------------
12 #include "scalelib.h"
13 module scale_atmos_phy_mp_tomita08
14  !-----------------------------------------------------------------------------
15  !
16  !++ used modules
17  !
18  use scale_precision
19  use scale_io
20  use scale_prof
21 
22  !-----------------------------------------------------------------------------
23  implicit none
24  private
25  !-----------------------------------------------------------------------------
26  !
27  !++ Public procedure
28  !
29  public :: atmos_phy_mp_tomita08_setup
30  public :: atmos_phy_mp_tomita08_finalize
31  public :: atmos_phy_mp_tomita08_putvar
32  public :: atmos_phy_mp_tomita08_getvar
33  public :: atmos_phy_mp_tomita08_adjustment
34  public :: atmos_phy_mp_tomita08_terminal_velocity
35  public :: atmos_phy_mp_tomita08_effective_radius
36  public :: atmos_phy_mp_tomita08_cloud_fraction
37  public :: atmos_phy_mp_tomita08_qtrc2qhyd
38  public :: atmos_phy_mp_tomita08_qhyd2qtrc
39 
40  !-----------------------------------------------------------------------------
41  !
42  !++ Public parameters & variables
43  !
44  integer, private, parameter :: QA_MP = 6
45 
46  integer, parameter, public :: atmos_phy_mp_tomita08_ntracers = qa_mp
47  integer, parameter, public :: atmos_phy_mp_tomita08_nwaters = 2
48  integer, parameter, public :: atmos_phy_mp_tomita08_nices = 3
49  character(len=H_SHORT), parameter, public :: atmos_phy_mp_tomita08_tracer_names(qa_mp) = (/ &
50  'QV', &
51  'QC', &
52  'QR', &
53  'QI', &
54  'QS', &
55  'QG' /)
56  character(len=H_MID) , parameter, public :: atmos_phy_mp_tomita08_tracer_descriptions(qa_mp) = (/ &
57  'Ratio of Water Vapor mass to total mass (Specific humidity)', &
58  'Ratio of Cloud Water mass to total mass ', &
59  'Ratio of Rain Water mass to total mass ', &
60  'Ratio of Cloud Ice mass ratio to total mass ', &
61  'Ratio of Snow miass ratio to total mass ', &
62  'Ratio of Graupel mass ratio to total mass '/)
63  character(len=H_SHORT), parameter, public :: atmos_phy_mp_tomita08_tracer_units(qa_mp) = (/ &
64  'kg/kg', &
65  'kg/kg', &
66  'kg/kg', &
67  'kg/kg', &
68  'kg/kg', &
69  'kg/kg' /)
70 
71  !-----------------------------------------------------------------------------
72  !
73  !++ Private procedure
74  !
75  private :: mp_tomita08
76 
77  !-----------------------------------------------------------------------------
78  !
79  !++ Private parameters & variables
80  !
81  integer, private, parameter :: i_qv = 1
82  integer, private, parameter :: i_qc = 2
83  integer, private, parameter :: i_qr = 3
84  integer, private, parameter :: i_qi = 4
85  integer, private, parameter :: i_qs = 5
86  integer, private, parameter :: i_qg = 6
87 
88  integer, private, parameter :: i_hyd_qc = 1
89  integer, private, parameter :: i_hyd_qr = 2
90  integer, private, parameter :: i_hyd_qi = 3
91  integer, private, parameter :: i_hyd_qs = 4
92  integer, private, parameter :: i_hyd_qg = 5
93 
94  logical, private :: do_couple_aerosol = .false. ! apply CCN effect?
95  logical, private :: do_explicit_icegen = .false. ! apply explicit ice generation?
96 
97  logical, private :: fixed_re = .false. ! use ice's effective radius for snow and graupel, and set rain transparent?
98  logical, private :: const_rec = .true. ! use constant effective radius for cloud water?
99  logical, private :: nofall_qr = .false. ! surpress sedimentation of rain?
100  logical, private :: nofall_qi = .false. ! surpress sedimentation of ice?
101  logical, private :: nofall_qs = .false. ! surpress sedimentation of snow?
102  logical, private :: nofall_qg = .false. ! surpress sedimentation of graupel?
103  !$acc declare create(nofall_qr, nofall_qi, nofall_qs, nofall_qg)
104 
105  real(rp), private, parameter :: dens00 = 1.28_rp
106 
107  ! Parameter for Marshall-Palmer distribution
108  real(rp), private :: n0r_def = 8.e+6_rp
109  real(rp), private :: n0s_def = 3.e+6_rp
110  real(rp), private :: n0g_def = 4.e+6_rp
111  !$acc declare create(N0r_def, N0s_def, N0g_def)
112 
113  real(rp), private :: dens_s = 100.0_rp
114  real(rp), private :: dens_g = 400.0_rp
115  ! graupel : 400
116  ! hail : 917
117  real(rp), private :: drag_g = 0.6_rp
118  real(rp), private :: re_qc = 8.e-6_rp ! effective radius for cloud water
119  real(rp), private :: re_qi = 40.e-6_rp ! effective radius for cloud ice
120  real(rp), private :: cr = 130.0_rp
121  real(rp), private :: cs = 4.84_rp
122  !$acc declare create(Cr, Cs)
123 
124  ! Empirical parameter
125  real(rp), private :: ar, as, ag
126  real(rp), private :: br, bs, bg
127  real(rp), private :: cg
128  real(rp), private :: dr, ds, dg
129  !$acc declare create(Ar, As, Ag, Cg)
130 
131  ! GAMMA function
132  real(rp), private :: gam, gam_2, gam_3
133 
134  real(rp), private :: gam_1br, gam_2br, gam_3br
135  real(rp), private :: gam_3dr
136  real(rp), private :: gam_6dr
137  real(rp), private :: gam_1brdr
138  real(rp), private :: gam_5dr_h
139 
140  real(rp), private :: gam_1bs, gam_2bs, gam_3bs
141  real(rp), private :: gam_3ds
142  real(rp), private :: gam_1bsds
143  real(rp), private :: gam_5ds_h
144 
145  real(rp), private :: gam_1bg, gam_3dg
146  real(rp), private :: gam_1bgdg
147  real(rp), private :: gam_5dg_h
148  !$acc declare create(GAM_1brdr, GAM_1bs, GAM_1bg, GAM_1br, GAM_1bsds, GAM_1bgdg)
149 
150  ! Bergeron Parameter
151  real(rp), parameter :: bergeron_a1_tab(32) = (/ &
152  0.0001e-7_rp, 0.7939e-7_rp, 0.7841e-6_rp, 0.3369e-5_rp, 0.4336e-5_rp, &
153  0.5285e-5_rp, 0.3728e-5_rp, 0.1852e-5_rp, 0.2991e-6_rp, 0.4248e-6_rp, &
154  0.7434e-6_rp, 0.1812e-5_rp, 0.4394e-5_rp, 0.9145e-5_rp, 0.1725e-4_rp, &
155  0.3348e-4_rp, 0.1725e-4_rp, 0.9175e-5_rp, 0.4412e-5_rp, 0.2252e-5_rp, &
156  0.9115e-6_rp, 0.4876e-6_rp, 0.3473e-6_rp, 0.4758e-6_rp, 0.6306e-6_rp, &
157  0.8573e-6_rp, 0.7868e-6_rp, 0.7192e-6_rp, 0.6513e-6_rp, 0.5956e-6_rp, &
158  0.5333e-6_rp, 0.4834e-6_rp /)
159  real(rp), parameter :: bergeron_a2_tab(32) = (/ &
160  0.0100_rp, 0.4006_rp, 0.4831_rp, 0.5320_rp, 0.5307_rp, &
161  0.5319_rp, 0.5249_rp, 0.4888_rp, 0.3849_rp, 0.4047_rp, &
162  0.4318_rp, 0.4771_rp, 0.5183_rp, 0.5463_rp, 0.5651_rp, &
163  0.5813_rp, 0.5655_rp, 0.5478_rp, 0.5203_rp, 0.4906_rp, &
164  0.4447_rp, 0.4126_rp, 0.3960_rp, 0.4149_rp, 0.4320_rp, &
165  0.4506_rp, 0.4483_rp, 0.4460_rp, 0.4433_rp, 0.4413_rp, &
166  0.4382_rp, 0.4361_rp /)
167  real(rp), private :: temcc, fact
168  integer, private :: itemc
169 
170  !---< Khairoutdinov and Kogan (2000) >---
171  logical, private :: enable_kk2000 = .false.
172 
173 
174  !---< Roh and Satoh (2014) >---
175  logical, private :: enable_rs2014 = .false.
176  real(rp), private, parameter :: ln10 = log(10.0_rp)
177  real(rp), private, parameter :: coef_a01 = 5.065339_rp
178  real(rp), private, parameter :: coef_a02 = -0.062659_rp
179  real(rp), private, parameter :: coef_a03 = -3.032362_rp
180  real(rp), private, parameter :: coef_a04 = 0.029469_rp
181  real(rp), private, parameter :: coef_a05 = -0.000285_rp
182  real(rp), private, parameter :: coef_a06 = 0.31255_rp
183  real(rp), private, parameter :: coef_a07 = 0.000204_rp
184  real(rp), private, parameter :: coef_a08 = 0.003199_rp
185  real(rp), private, parameter :: coef_a09 = 0.0_rp
186  real(rp), private, parameter :: coef_a10 = -0.015952_rp
187 
188  real(rp), private, parameter :: coef_b01 = 0.476221_rp
189  real(rp), private, parameter :: coef_b02 = -0.015896_rp
190  real(rp), private, parameter :: coef_b03 = 0.165977_rp
191  real(rp), private, parameter :: coef_b04 = 0.007468_rp
192  real(rp), private, parameter :: coef_b05 = -0.000141_rp
193  real(rp), private, parameter :: coef_b06 = 0.060366_rp
194  real(rp), private, parameter :: coef_b07 = 0.000079_rp
195  real(rp), private, parameter :: coef_b08 = 0.000594_rp
196  real(rp), private, parameter :: coef_b09 = 0.0_rp
197  real(rp), private, parameter :: coef_b10 = -0.003577_rp
198  !$acc declare create(enable_RS2014)
199 
200  !---< Wainwright et al. (2014) >---
201  logical, private :: enable_wdxz2014 = .false.
202  !$acc declare create(enable_WDXZ2014)
203 
204  !---< Heymsfield et al. (2007) >---
205  logical, private :: enable_hzdfhi2007 = .false. ! use scheme by Heymsfield et al. (2007)
206  character(len=12), private :: cloud_type = "synoptic" ! "synoptic" or "crystal_face"
207  real(rp), private :: coef_a0
208  real(rp), private :: coef_a1
209  real(rp), private :: coef_b0
210  real(rp), private :: coef_b1
211  !$acc declare create(enable_HZDFHI2007, coef_a0, coef_a1, coef_b0, coef_b1)
212 
213  !---< Thornberry et al. (2017) >---
214  logical, private :: enable_trawlbg2017 = .false. ! use scheme by Thornberry et al. (2017)
215 
216  ! Accretion parameter
217  real(rp), private :: eiw = 1.0_rp
218  real(rp), private :: erw = 1.0_rp
219  real(rp), private :: esw = 1.0_rp
220  real(rp), private :: egw = 1.0_rp
221  real(rp), private :: eri = 1.0_rp
222  real(rp), private :: esi = 1.0_rp
223  real(rp), private :: egi = 0.1_rp
224  real(rp), private :: esr = 1.0_rp
225  real(rp), private :: egr = 1.0_rp
226  real(rp), private :: egs = 1.0_rp
227  real(rp), private :: gamma_sacr = 25.e-3_rp
228  real(rp), private :: gamma_gacs = 90.e-3_rp
229  real(rp), private :: mi = 4.19e-13_rp
230 
231  ! Auto-conversion parameter
232  real(rp), private, parameter :: nc_lnd = 2000.0_rp
233  real(rp), private, parameter :: nc_ocn = 50.0_rp
234  real(rp), private, allocatable :: nc_def(:,:)
235  !$acc declare create(Nc_def)
236 
237  real(rp), private :: beta_saut = 6.e-3_rp
238  real(rp), private :: gamma_saut = 60.e-3_rp
239  real(rp), private :: beta_gaut = 0.0_rp
240  real(rp), private :: gamma_gaut = 90.e-3_rp
241  real(rp), private :: qicrt_saut = 0.0_rp
242  real(rp), private :: qscrt_gaut = 6.e-4_rp
243 
244  ! Evaporation, Sublimation parameter
245  real(rp), private, parameter :: da0 = 2.428e-2_rp
246  real(rp), private, parameter :: dda_dt = 7.47e-5_rp
247  real(rp), private, parameter :: dw0 = 2.222e-5_rp
248  real(rp), private, parameter :: ddw_dt = 1.37e-7_rp
249  real(rp), private, parameter :: mu0 = 1.718e-5_rp
250  real(rp), private, parameter :: dmu_dt = 5.28e-8_rp
251 
252  real(rp), private, parameter :: f1r = 0.78_rp
253  real(rp), private, parameter :: f2r = 0.27_rp
254  real(rp), private, parameter :: f1s = 0.65_rp
255  real(rp), private, parameter :: f2s = 0.39_rp
256  real(rp), private, parameter :: f1g = 0.78_rp
257  real(rp), private, parameter :: f2g = 0.27_rp
258 
259  ! Freezing parameter
260  real(rp), private, parameter :: a_frz = 0.66_rp
261  real(rp), private, parameter :: b_frz = 100.0_rp
262 
263  ! Bergeron process parameter
264  real(rp), private, parameter :: mi40 = 2.46e-10_rp
265  real(rp), private, parameter :: mi50 = 4.80e-10_rp
266  real(rp), private, parameter :: vti50 = 1.0_rp
267  real(rp), private, parameter :: ri50 = 5.e-5_rp
268 
269  ! Explicit ice generation
270  logical, private :: only_liquid = .false.
271  real(rp), private :: sw_expice = 0.0_rp
272  real(rp), private, parameter :: nc_ihtr = 300.0_rp
273  real(rp), private, parameter :: di_max = 500.e-6_rp
274  real(rp), private, parameter :: di_a = 11.9_rp
275 
276  ! For coalesence coefficient of ice and snow
277  real(rp), private :: ecoal_gi
278  real(rp), private :: ecoal_gs
279  real(rp), private :: flg_ecoali = 0.0_rp
280  real(rp), private :: flg_ecoals = 0.0_rp
281 
282  integer, private, parameter :: w_nmax = 49
283  integer, private, parameter :: i_dqv_dt = 1 !
284  integer, private, parameter :: i_dqc_dt = 2 !
285  integer, private, parameter :: i_dqr_dt = 3 !
286  integer, private, parameter :: i_dqi_dt = 4 !
287  integer, private, parameter :: i_dqs_dt = 5 !
288  integer, private, parameter :: i_dqg_dt = 6 !
289  integer, private, parameter :: i_delta1 = 7 ! separation switch for r->s,g
290  integer, private, parameter :: i_delta2 = 8 ! separation switch for s->g
291  integer, private, parameter :: i_spsati = 9 ! separation switch for ice sublimation
292  integer, private, parameter :: i_iceflg = 10 ! separation switch for T > 0
293  integer, private, parameter :: i_rlmdr = 11
294  integer, private, parameter :: i_rlmds = 12
295  integer, private, parameter :: i_rlmdg = 13
296  integer, private, parameter :: i_piacr = 14 ! r->s,g
297  integer, private, parameter :: i_psacr = 15 ! r->s,g
298  integer, private, parameter :: i_praci = 16 ! i->s,g
299  integer, private, parameter :: i_pigen = 17 ! v->i
300  integer, private, parameter :: i_pidep = 18 ! v->i
301  integer, private, parameter :: i_psdep = 19 ! v->s
302  integer, private, parameter :: i_pgdep = 20 ! v->g
303  integer, private, parameter :: i_praut = 21 ! c->r
304  integer, private, parameter :: i_pracw = 22 ! c->r
305  integer, private, parameter :: i_pihom = 23 ! c->i
306  integer, private, parameter :: i_pihtr = 24 ! c->i
307  integer, private, parameter :: i_psacw = 25 ! c->s
308  integer, private, parameter :: i_psfw = 26 ! c->s
309  integer, private, parameter :: i_pgacw = 27 ! c->g
310  integer, private, parameter :: i_prevp = 28 ! r->v
311  integer, private, parameter :: i_piacr_s = 29 ! r->s
312  integer, private, parameter :: i_psacr_s = 30 ! r->s
313  integer, private, parameter :: i_piacr_g = 31 ! r->g
314  integer, private, parameter :: i_psacr_g = 32 ! r->g
315  integer, private, parameter :: i_pgacr = 33 ! r->g
316  integer, private, parameter :: i_pgfrz = 34 ! r->g
317  integer, private, parameter :: i_pisub = 35 ! i->v
318  integer, private, parameter :: i_pimlt = 36 ! i->c
319  integer, private, parameter :: i_psaut = 37 ! i->s
320  integer, private, parameter :: i_praci_s = 38 ! i->s
321  integer, private, parameter :: i_psaci = 39 ! i->s
322  integer, private, parameter :: i_psfi = 40 ! i->s
323  integer, private, parameter :: i_praci_g = 41 ! i->g
324  integer, private, parameter :: i_pgaci = 42 ! i->g
325  integer, private, parameter :: i_pssub = 43 ! s->v
326  integer, private, parameter :: i_psmlt = 44 ! s->r
327  integer, private, parameter :: i_pgaut = 45 ! s->g
328  integer, private, parameter :: i_pracs = 46 ! s->g
329  integer, private, parameter :: i_pgacs = 47 ! s->g
330  integer, private, parameter :: i_pgsub = 48 ! g->v
331  integer, private, parameter :: i_pgmlt = 49 ! g->r
332 
333  character(len=H_SHORT), private :: w_name(w_nmax)
334 
335  data w_name / 'dqv_dt ', &
336  'dqc_dt ', &
337  'dqr_dt ', &
338  'dqi_dt ', &
339  'dqs_dt ', &
340  'dqg_dt ', &
341  'delta1 ', &
342  'delta2 ', &
343  'spsati ', &
344  'iceflg ', &
345  'RLMDr ', &
346  'RLMDs ', &
347  'RLMDg ', &
348  'Piacr ', &
349  'Psacr ', &
350  'Praci ', &
351  'Pigen ', &
352  'Pidep ', &
353  'Psdep ', &
354  'Pgdep ', &
355  'Praut ', &
356  'Pracw ', &
357  'Pihom ', &
358  'Pihtr ', &
359  'Psacw ', &
360  'Psfw ', &
361  'Pgacw ', &
362  'Prevp ', &
363  'Piacr_s', &
364  'Psacr_s', &
365  'Piacr_g', &
366  'Psacr_g', &
367  'Pgacr ', &
368  'Pgfrz ', &
369  'Pisub ', &
370  'Pimlt ', &
371  'Psaut ', &
372  'Praci_s', &
373  'Psaci ', &
374  'Psfi ', &
375  'Praci_g', &
376  'Pgaci ', &
377  'Pssub ', &
378  'Psmlt ', &
379  'Pgaut ', &
380  'Pracs ', &
381  'Pgacs ', &
382  'Pgsub ', &
383  'Pgmlt ' /
384 
385  real(rp), private, allocatable :: w3d(:,:,:,:)
386  integer, private :: hist_id(w_nmax)
387  integer, private :: hist_pcsat, hist_pisat
388  !$acc declare create(w3d)
389 
390  !-----------------------------------------------------------------------------
391 contains
392  !-----------------------------------------------------------------------------
396  subroutine atmos_phy_mp_tomita08_setup( &
397  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
398  flg_lt )
399  use scale_prc, only: &
400  prc_abort
401  use scale_const, only: &
402  pi => const_pi, &
403  grav => const_grav, &
404  dens_w => const_dwatr, &
405  dens_i => const_dice
406  use scale_specfunc, only: &
407  sf_gamma
408  use scale_file_history, only: &
409  file_history_reg
410  implicit none
411 
412  integer, intent(in) :: ka, ks, ke
413  integer, intent(in) :: ia, is, ie
414  integer, intent(in) :: ja, js, je
415 
416  logical, intent(in), optional :: flg_lt
417 
418  real(rp) :: autoconv_nc = nc_ocn
419 
420  namelist / param_atmos_phy_mp_tomita08 / &
421  do_couple_aerosol, &
422  do_explicit_icegen, &
423  autoconv_nc, &
424  enable_kk2000, &
425  enable_rs2014, &
426  enable_wdxz2014, &
427  enable_hzdfhi2007, &
428  cloud_type, &
429  enable_trawlbg2017, &
430  n0r_def, &
431  n0s_def, &
432  n0g_def, &
433  dens_s, &
434  dens_g, &
435  re_qc, &
436  re_qi, &
437  drag_g, &
438  cr, &
439  cs, &
440  eiw, &
441  erw, &
442  esw, &
443  egw, &
444  eri, &
445  esi, &
446  egi, &
447  esr, &
448  egr, &
449  egs, &
450  gamma_sacr, &
451  gamma_gacs, &
452  mi, &
453  beta_saut, &
454  gamma_saut, &
455  qicrt_saut, &
456  beta_gaut, &
457  gamma_gaut, &
458  qscrt_gaut, &
459  fixed_re, &
460  const_rec, &
461  nofall_qr, &
462  nofall_qi, &
463  nofall_qs, &
464  nofall_qg, &
465  ecoal_gi, &
466  ecoal_gs
467 
468  real(rp), parameter :: max_term_vel = 10.0_rp !-- terminal velocity for calculate dt of sedimentation
469 
470  logical :: flg_lt_
471  integer :: ierr
472  integer :: i, j, ip
473  !---------------------------------------------------------------------------
474 
475 
476  log_newline
477  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Setup'
478  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Tomita (2008) 1-moment bulk 6 category'
479 
480  allocate( w3d(ka,ia,ja,w_nmax) )
481  w3d(:,:,:,:) = 0.0_rp
482 
483  allocate( nc_def(ia,ja) )
484 
485  ecoal_gi = 0.0_rp
486  ecoal_gs = 0.0_rp
487 
488  !--- read namelist
489  rewind(io_fid_conf)
490  read(io_fid_conf,nml=param_atmos_phy_mp_tomita08,iostat=ierr)
491  if( ierr < 0 ) then !--- missing
492  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Not found namelist. Default used.'
493  elseif( ierr > 0 ) then !--- fatal error
494  log_error("ATMOS_PHY_MP_tomita08_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_TOMITA08. Check!'
495  call prc_abort
496  endif
497  log_nml(param_atmos_phy_mp_tomita08)
498 
499  log_newline
500  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'density of the snow [kg/m3] : ', dens_s
501  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'density of the graupel [kg/m3] : ', dens_g
502  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Nc for auto-conversion [num/m3]: ', autoconv_nc
503  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Use k-k scheme? : ', enable_kk2000
504  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Use Roh scheme? : ', enable_rs2014
505  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Use WDXZ scheme? : ', enable_wdxz2014
506  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Use HZDFHI scheme? : ', enable_hzdfhi2007
507  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Use TRAWLBG scheme? : ', enable_trawlbg2017
508  log_newline
509  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Use effective radius of ice for snow and graupel,'
510  log_info("ATMOS_PHY_MP_tomita08_setup",*) ' and set rain transparent? : ', fixed_re
511  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Density of the ice is used for the calculation of '
512  log_info("ATMOS_PHY_MP_tomita08_setup",*) ' optically effective volume of snow and graupel.'
513  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Surpress sedimentation of rain? : ', nofall_qr
514  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Surpress sedimentation of ice? : ', nofall_qi
515  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Surpress sedimentation of snow? : ', nofall_qs
516  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Surpress sedimentation of graupel? : ', nofall_qg
517  log_info("ATMOS_PHY_MP_tomita08_setup",*) 'Enable explicit ice generation? : ', do_explicit_icegen
518  log_newline
519 
520  do j = js, je
521  do i = is, ie
522  nc_def(i,j) = autoconv_nc
523  end do
524  end do
525 
526  !--- empirical coefficients A, B, C, D
527  ar = pi * dens_w / 6.0_rp
528  as = pi * dens_s / 6.0_rp
529  ag = pi * dens_g / 6.0_rp
530 
531  br = 3.0_rp
532  bs = 3.0_rp
533  bg = 3.0_rp
534 
535  cg = sqrt( ( 4.0_rp * dens_g * grav ) / ( 3.0_rp * dens00 * drag_g ) )
536 
537  dr = 0.50_rp
538  ds = 0.25_rp
539  dg = 0.50_rp
540 
541  if ( enable_rs2014 ) then ! overwrite parameters
542  do_explicit_icegen = .true.
543 
544  n0g_def = 4.e+8_rp
545  as = 0.069_rp
546  bs = 2.0_rp
547  esi = 0.25_rp
548  egi = 0.0_rp
549  egs = 0.0_rp
550  endif
551 
552  if ( enable_hzdfhi2007 ) then
553  select case( cloud_type )
554  case ( "synoptic" )
555  coef_a0 = 120.4_rp
556  coef_a1 = 2.21_rp
557  coef_b0 = 0.0487_rp
558  coef_b1 = 4.57e-4_rp
559  case ("crystal_face")
560  coef_a0 = 135.3_rp
561  coef_a1 = 1.93_rp
562  coef_b0 = 0.0058_rp
563  coef_b1 = -0.0024539_rp
564  case default
565  log_error("ATMOS_PHY_MP_tomita08_setup",*) 'cloud_type is invalid: ', trim(cloud_type)
566  call prc_abort
567  end select
568  end if
569 
570  if ( do_explicit_icegen ) then
571  only_liquid = .true.
572  sw_expice = 1.0_rp
573  else
574  only_liquid = .false.
575  sw_expice = 0.0_rp
576  endif
577 
578  gam = 1.0_rp ! =0!
579  gam_2 = 1.0_rp ! =1!
580  gam_3 = 2.0_rp ! =2!
581 
582  gam_1br = sf_gamma( 1.0_rp + br ) ! = 4!
583  gam_2br = sf_gamma( 2.0_rp + br ) ! = 5!
584  gam_3br = sf_gamma( 3.0_rp + br ) ! = 6!
585  gam_3dr = sf_gamma( 3.0_rp + dr )
586  gam_6dr = sf_gamma( 6.0_rp + dr )
587  gam_1brdr = sf_gamma( 1.0_rp + br + dr )
588  gam_5dr_h = sf_gamma( 0.5_rp * (5.0_rp+dr) )
589 
590  gam_1bs = sf_gamma( 1.0_rp + bs ) ! = 4!
591  gam_2bs = sf_gamma( 2.0_rp + bs ) ! = 5!
592  gam_3bs = sf_gamma( 3.0_rp + bs ) ! = 6!
593  gam_3ds = sf_gamma( 3.0_rp + ds )
594  gam_1bsds = sf_gamma( 1.0_rp + bs + ds )
595  gam_5ds_h = sf_gamma( 0.5_rp * (5.0_rp+ds) )
596 
597  gam_1bg = sf_gamma( 1.0_rp + bg ) ! = 4!
598  gam_3dg = sf_gamma( 3.0_rp + dg )
599  gam_1bgdg = sf_gamma( 1.0_rp + bg + dg)
600  gam_5dg_h = sf_gamma( 0.5_rp * (5.0_rp+dg) )
601 
602  ! history
603  do ip = 1, w_nmax
604  call file_history_reg( w_name(ip), 'individual tendency term in tomita08', 'kg/kg/s', & ! [IN]
605  hist_id(ip) ) ! [OUT]
606  end do
607  call file_history_reg( 'Pcsat', 'QC production term by satulation adjustment', 'kg/kg/s', hist_pcsat )
608  call file_history_reg( 'Pisat', 'QI production term by satulation adjustment', 'kg/kg/s', hist_pisat )
609 
610 
611  if( present(flg_lt) ) then
612  flg_lt_ = flg_lt
613  else
614  flg_lt_ = .false.
615  end if
616  if ( flg_lt_ ) then
617  if( enable_rs2014 ) then
618  write(*,*) 'xxx ROH and Satoh (2014) scheme is not supported for Lightning'
619  call prc_abort
620  endif
621 
622  if( ecoal_gi == 0.0_rp ) then
623  flg_ecoali = 0.0_rp
624  ecoal_gi = egi
625  else
626  flg_ecoali = 1.0_rp
627  endif
628  if( ecoal_gs == 0.0_rp ) then
629  flg_ecoals = 0.0_rp
630  ecoal_gs = egs
631  else
632  flg_ecoals = 1.0_rp
633  endif
634 
635  else
636  ecoal_gi = 1.0_rp
637  ecoal_gs = 1.0_rp
638  endif
639 
640  !$acc update device(Cr, Cs)
641  !$acc update device(Ar, As, Ag, Cg)
642  !$acc update device(nofall_qr, nofall_qi, nofall_qs, nofall_qg)
643  !$acc update device(N0r_def, N0s_def, N0g_def)
644  !$acc update device(Nc_def)
645  !$acc update device(GAM_1brdr, GAM_1bs, GAM_1bg, GAM_1br, GAM_1bsds, GAM_1bgdg)
646  !$acc update device(enable_RS2014)
647  !$acc update device(enable_WDXZ2014)
648  !$acc update device(enable_HZDFHI2007, coef_a0, coef_a1, coef_b0, coef_b1)
649 
650  return
651  end subroutine atmos_phy_mp_tomita08_setup
652 
653  !-----------------------------------------------------------------------------
655  subroutine atmos_phy_mp_tomita08_finalize
656 
657  deallocate( w3d )
658  deallocate( nc_def )
659 
660  return
661  end subroutine atmos_phy_mp_tomita08_finalize
662 
663  !-----------------------------------------------------------------------------
665  subroutine atmos_phy_mp_tomita08_putvar( &
666  in_drag_g, &
667  in_re_qc, &
668  in_re_qi, &
669  in_Cr, &
670  in_Cs )
671  use scale_const, only: &
672  grav => const_grav
673  implicit none
674 
675  real(rp), intent(in), optional :: in_drag_g
676  real(rp), intent(in), optional :: in_re_qc
677  real(rp), intent(in), optional :: in_re_qi
678  real(rp), intent(in), optional :: in_cr
679  real(rp), intent(in), optional :: in_cs
680  !---------------------------------------------------------------------------
681 
682  if( present(in_drag_g) ) drag_g = in_drag_g
683  if( present(in_re_qc ) ) re_qc = in_re_qc
684  if( present(in_re_qi ) ) re_qi = in_re_qi
685  if( present(in_cr ) ) cr = in_cr
686  if( present(in_cs ) ) cs = in_cs
687 
688  if( present(in_drag_g) ) then
689  cg = sqrt( ( 4.0_rp * dens_g * grav ) / ( 3.0_rp * dens00 * drag_g ) ) ! recalc. Cg
690  end if
691 
692  return
693  end subroutine atmos_phy_mp_tomita08_putvar
694 
695  !-----------------------------------------------------------------------------
697  subroutine atmos_phy_mp_tomita08_getvar( &
698  out_drag_g, &
699  out_re_qc, &
700  out_re_qi, &
701  out_Cr, &
702  out_Cs )
703  implicit none
704 
705  real(rp), intent(out), optional :: out_drag_g
706  real(rp), intent(out), optional :: out_re_qc
707  real(rp), intent(out), optional :: out_re_qi
708  real(rp), intent(out), optional :: out_cr
709  real(rp), intent(out), optional :: out_cs
710  !---------------------------------------------------------------------------
711 
712  if( present(out_drag_g) ) out_drag_g = drag_g
713  if( present(out_re_qc ) ) out_re_qc = re_qc
714  if( present(out_re_qi ) ) out_re_qi = re_qi
715  if( present(out_cr ) ) out_cr = cr
716  if( present(out_cs ) ) out_cs = cs
717 
718  return
719  end subroutine atmos_phy_mp_tomita08_getvar
720 
721  !-----------------------------------------------------------------------------
725  subroutine atmos_phy_mp_tomita08_adjustment( &
726  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
727  DENS, PRES, &
728  CCN, &
729  dt, &
730  TEMP, QTRC, CPtot, CVtot, &
731  RHOE_t, EVAPORATE, &
732  flg_lt, d0_crg, v0_crg, &
733  dqcrg, beta_crg, &
734  QSPLT_in, Sarea, QTRC_crg )
735  use scale_const, only: &
736  dwatr => const_dwatr, &
737  pi => const_pi
738  use scale_file_history, only: &
739  file_history_query, &
740  file_history_put
741  use scale_atmos_phy_mp_common, only: &
742  mp_saturation_adjustment => atmos_phy_mp_saturation_adjustment
743  implicit none
744  integer, intent(in) :: ka, ks, ke
745  integer, intent(in) :: ia, is, ie
746  integer, intent(in) :: ja, js, je
747 
748  real(rp), intent(in) :: dens (ka,ia,ja)
749  real(rp), intent(in) :: pres (ka,ia,ja)
750  real(rp), intent(in) :: ccn (ka,ia,ja)
751  real(dp), intent(in) :: dt
752 
753  real(rp), intent(inout) :: temp(ka,ia,ja)
754  real(rp), intent(inout) :: qtrc(ka,ia,ja,qa_mp)
755  real(rp), intent(inout) :: cptot(ka,ia,ja)
756  real(rp), intent(inout) :: cvtot(ka,ia,ja)
757 
758  real(rp), intent(out) :: rhoe_t (ka,ia,ja)
759  real(rp), intent(out) :: evaporate(ka,ia,ja) ! number of evaporated cloud [/m3/s]
760 
761  ! Optional for Lightning
762  logical, intent(in), optional :: flg_lt
763  real(rp), intent(in), optional :: d0_crg, v0_crg
764  real(rp), intent(in), optional :: dqcrg(ka,ia,ja)
765  real(rp), intent(in), optional :: beta_crg(ka,ia,ja)
766  real(rp), intent(out), optional :: sarea(ka,ia,ja,qa_mp-1) ! Surface area of each hydrometeor
767  real(rp), intent(out), optional :: qsplt_in(ka,ia,ja,3) ! Charge separation
768  real(rp), intent(inout), optional :: qtrc_crg(ka,ia,ja,qa_mp-1) ! Tracer for charge density
769 
770  real(rp) :: rhoe_d_sat(ka,ia,ja)
771 
772  real(rp) :: qc_t_sat(ka,ia,ja)
773  real(rp) :: qi_t_sat(ka,ia,ja)
774  logical :: hist_flag
775 
776  integer :: k, i, j
777  !---------------------------------------------------------------------------
778 
779  log_progress(*) 'atmosphere / physics / microphysics / Tomita08'
780 
781  !$acc data copy(TEMP, QTRC, CPtot, CVtot), &
782  !$acc copyin(KA, KS, KE, IA, IS, IE, JA, JS, JE, &
783  !$acc DENS, PRES, CCN, dt ) &
784  !$acc copyout(RHOE_t, EVAPORATE), &
785  !$acc create(RHOE_d_sat, QC_t_sat, QI_t_sat)
786 
787 
788  !##### MP Main #####
789  call mp_tomita08( &
790  ka, ks, ke, ia, is, ie, ja, js, je, & ! [IN]
791  dens(:,:,:), pres(:,:,:), ccn(:,:,:), & ! [IN]
792  dt, & ! [IN]
793  temp(:,:,:), qtrc(:,:,:,:), & ! [INOUT]
794  cptot(:,:,:), cvtot(:,:,:), & ! [INOUT]
795  rhoe_t(:,:,:), & ! [OUT]
796  flg_lt, d0_crg, v0_crg, & ! [IN:Optional]
797  dqcrg(:,:,:), beta_crg(:,:,:), & ! [IN:Optional]
798  qsplt_in(:,:,:,:), & ! [OUT:Optional]
799  sarea(:,:,:,:), & ! [OUT:Optional]
800  qtrc_crg(:,:,:,:) ) ! [INOUT:Optional]
801 
802  ! save value before saturation adjustment
803  !$omp parallel do collapse(2)
804  !$acc kernels
805  do j = js, je
806  do i = is, ie
807  do k = ks, ke
808  qc_t_sat(k,i,j) = qtrc(k,i,j,i_qc)
809  qi_t_sat(k,i,j) = qtrc(k,i,j,i_qi)
810  enddo
811  enddo
812  enddo
813  !$acc end kernels
814 
815  call mp_saturation_adjustment( &
816  ka, ks, ke, ia, is, ie, ja, js, je, &
817  dens(:,:,:), & ! [IN]
818  only_liquid, & ! [IN]
819  temp(:,:,:), & ! [INOUT]
820  qtrc(:,:,:,i_qv), & ! [INOUT]
821  qtrc(:,:,:,i_qc), qtrc(:,:,:,i_qi), & ! [INOUT]
822  cptot(:,:,:), cvtot(:,:,:), & ! [INOUT]
823  rhoe_d_sat(:,:,:) ) ! [OUT]
824 
825  !$omp parallel do collapse(2)
826  !$acc kernels
827  do j = js, je
828  do i = is, ie
829  do k = ks, ke
830  rhoe_t(k,i,j) = rhoe_t(k,i,j) + rhoe_d_sat(k,i,j) / dt
831  enddo
832  enddo
833  enddo
834  !$acc end kernels
835 
836  !$omp parallel do collapse(2)
837  !$acc kernels
838  do j = js, je
839  do i = is, ie
840  do k = ks, ke
841  qc_t_sat(k,i,j) = ( qtrc(k,i,j,i_qc) - qc_t_sat(k,i,j) ) / dt
842  evaporate(k,i,j) = max( -qc_t_sat(k,i,j), 0.0_rp ) & ! if negative, condensation
843  * dens(k,i,j) / (4.0_rp/3.0_rp*pi*dwatr*re_qc**3) ! mass -> number (assuming constant particle radius as re_qc)
844 
845  enddo
846  enddo
847  enddo
848  !$acc end kernels
849 
850  call file_history_query( hist_pcsat, hist_flag )
851  if ( hist_flag ) then
852  call file_history_put( hist_pcsat, qc_t_sat(:,:,:) )
853  end if
854 
855 
856  call file_history_query( hist_pisat, hist_flag )
857  if ( hist_flag ) then
858  !$omp parallel do collapse(2)
859  !$acc kernels
860  do j = js, je
861  do i = is, ie
862  do k = ks, ke
863  qi_t_sat(k,i,j) = ( qtrc(k,i,j,i_qi) - qi_t_sat(k,i,j) ) / dt
864  enddo
865  enddo
866  enddo
867  !$acc end kernels
868  call file_history_put( hist_pisat, qi_t_sat(:,:,:) )
869  end if
870 
871 
872  !$acc end data
873 
874  !##### END MP Main #####
875 
876  return
877  end subroutine atmos_phy_mp_tomita08_adjustment
878 
879  !-----------------------------------------------------------------------------
881  subroutine mp_tomita08( &
882  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
883  DENS, PRES, CCN, &
884  dt, &
885  TEMP, QTRC, &
886  CPtot0, CVtot0, &
887  RHOE_t, &
888  flg_lt, &
889  d0_crg, v0_crg, &
890  dqcrg, &
891  beta_crg, &
892  QSPLT_in, &
893  Sarea, &
894  QTRC_crg0 )
895  use scale_const, only: &
896  undef => const_undef, &
897  pi => const_pi, &
898  eps => const_eps, &
899  rvap => const_rvap, &
900  cl => const_cl, &
901  lhv0 => const_lhv0, &
902  lhs0 => const_lhs0, &
903  lhf0 => const_lhf0, &
904  tem00 => const_tem00, &
905  pre00 => const_pre00, &
906  dwatr => const_dwatr, &
907  dice => const_dice
908  use scale_file_history, only: &
909  file_history_query, &
910  file_history_put
911  use scale_atmos_hydrometeor, only: &
912  n_hyd, &
913  i_hs, &
914  lhv, &
915  lhf, &
916  cp_vapor, &
917  cp_water, &
918  cp_ice, &
919  cv_vapor, &
920  cv_water, &
921  cv_ice
922  use scale_atmos_saturation, only: &
923  saturation_dens2qsat_liq => atmos_saturation_dens2qsat_liq, &
924  saturation_dens2qsat_ice => atmos_saturation_dens2qsat_ice
925  implicit none
926  integer, intent(in) :: ka, ks, ke
927  integer, intent(in) :: ia, is, ie
928  integer, intent(in) :: ja, js, je
929 
930  real(rp), intent(in) :: dens(ka,ia,ja) ! density [km/m3]
931  real(rp), intent(in) :: pres(ka,ia,ja) ! pressure [Pa]
932  real(rp), intent(in) :: ccn (ka,ia,ja) ! [1/m3]
933  real(dp), intent(in) :: dt
934 
935  real(rp), intent(inout) :: temp (ka,ia,ja) ! temperature [K]
936  real(rp), intent(inout) :: qtrc (ka,ia,ja,qa_mp)
937  real(rp), intent(inout) :: cptot0(ka,ia,ja)
938  real(rp), intent(inout) :: cvtot0(ka,ia,ja)
939 
940  real(rp), intent(out) :: rhoe_t(ka,ia,ja)
941 
942  logical, intent(in), optional :: flg_lt
943  real(rp), intent(in), optional :: d0_crg, v0_crg
944  real(rp), intent(in), optional :: dqcrg(ka,ia,ja)
945  real(rp), intent(in), optional :: beta_crg(ka,ia,ja)
946  real(rp), intent(out), optional :: qsplt_in(ka,ia,ja,3)
947  real(rp), intent(out), optional :: sarea(ka,ia,ja,qa_mp-1)
948  real(rp), intent(inout), optional :: qtrc_crg0(ka,ia,ja,qa_mp-1)
949 
950  real(rp) :: cvtot(ks:ke)
951  real(rp) :: qv(ks:ke), qc(ks:ke), qr(ks:ke), qi(ks:ke), qs(ks:ke), qg(ks:ke)
952  real(rp) :: qv_t(ks:ke), qc_t(ks:ke), qr_t(ks:ke), qi_t(ks:ke), qs_t(ks:ke), qg_t(ks:ke)
953  real(rp) :: e_t, cp_t, cv_t
954 
955  real(rp) :: qsatl(ka) ! saturated water vapor for liquid water [kg/kg]
956  real(rp) :: qsati(ka) ! saturated water vapor for ice water [kg/kg]
957 
958  real(rp) :: sliq(ks:ke) ! saturated ratio S for liquid water (0-1)
959  real(rp) :: sice(ks:ke) ! saturated ratio S for ice water (0-1)
960 
961  real(rp) :: rho_fact(ks:ke) ! density factor
962  real(rp) :: temc(ks:ke) ! T - T0 [K]
963 
964  real(rp) :: n0r(ks:ke), n0s(ks:ke), n0g(ks:ke)
965 
966  real(rp) :: rlmdr(ks:ke), rlmdr_2(ks:ke), rlmdr_3(ks:ke)
967  real(rp) :: rlmds, rlmds_2, rlmds_3
968  real(rp) :: rlmdg(ks:ke), rlmdg_2(ks:ke), rlmdg_3(ks:ke)
969  real(rp) :: rlmdr_1br(ks:ke), rlmdr_2br(ks:ke), rlmdr_3br(ks:ke)
970  real(rp) :: rlmds_1bs, rlmds_2bs, rlmds_3bs
971  real(rp) :: rlmdr_dr(ks:ke), rlmdr_3dr(ks:ke), rlmdr_5dr(ks:ke)
972  real(rp) :: rlmds_ds, rlmds_3ds, rlmds_5ds
973  real(rp) :: rlmdg_dg(ks:ke), rlmdg_3dg(ks:ke), rlmdg_5dg(ks:ke)
974  real(rp) :: rlmdr_7(ks:ke)
975  real(rp) :: rlmdr_6dr(ks:ke)
976 
977  !---< Roh and Satoh (2014) >---
978  real(rp) :: tems, xs2
979  real(rp) :: moms_0(ks:ke), moms_1(ks:ke), moms_2(ks:ke)
980  real(rp) :: moms_0bs(ks:ke), moms_1bs(ks:ke), moms_2bs(ks:ke)
981  real(rp) :: moms_2ds(ks:ke), moms_5ds_h(ks:ke), rmoms_vt(ks:ke)
982  real(rp) :: coef_at(4), coef_bt(4)
983  real(rp) :: loga_, b_, nm
984 
985  real(rp) :: vti(ks:ke), vtr(ks:ke), vts(ks:ke), vtg(ks:ke)
986  real(rp) :: esi_mod, egs_mod(ka)
987  real(rp) :: rhoqc
988  real(rp) :: nc(ks:ke)
989  real(rp) :: pracw_orig, pracw_kk
990  real(rp) :: praut_berry, praut_kk
991  real(rp) :: dc
992  real(rp) :: betai, betas
993  real(rp) :: da
994  real(rp) :: kd
995  real(rp) :: nu(ks:ke)
996  real(rp) :: glv(ks:ke), giv(ks:ke), gil(ks:ke)
997  real(rp) :: ventr, vents, ventg
998  real(rp) :: net, fac, fack(ka), fac_sw
999  real(rp) :: zerosw, tmp
1000 
1001  !---< Bergeron process >---
1002  real(rp) :: sw_bergeron(ks:ke)
1003  real(rp) :: a1, a2
1004  real(rp) :: ma2
1005  real(rp) :: dt1
1006  real(rp) :: ni50
1007 
1008  !---< Explicit ice generation >---
1009  real(rp) :: sw, rhoqi, xni, xmi, di, nig, qig
1010 
1011  logical :: hist_sw(w_nmax), hist_flag
1012  real(rp) :: w(ks:ke,w_nmax)
1013 
1014  integer :: k, i, j, ip
1015 
1016  ! for lightning
1017  integer, parameter :: i_qgaci = 1
1018  integer, parameter :: i_qgacs = 2
1019  real(rp) :: qcrg_c(ks:ke), qcrg_r(ks:ke)
1020  real(rp) :: qcrg_i(ks:ke), qcrg_s(ks:ke), qcrg_g(ks:ke)
1021  real(rp) :: w_q(ks:ke,w_nmax)
1022  real(rp) :: w_qcrg(ks:ke,2) ! charge separation [fC/kg]
1023  real(rp) :: rdens_r, rdens_i, rdens_s, rdens_g
1024  real(rp) :: dcrg(ks:ke), beta1_crg(ks:ke), re_qs(ks:ke)
1025  real(rp) :: re(ka,ia,ja,n_hyd)
1026  real(rp) :: alpha
1027  real(rp) :: facq_qc, facq_qr, facq_qi, facq_qs, facq_qg
1028  logical :: flg_lt_l
1029  integer :: pp, qq, iq
1030  real(rp) :: qc_crg_t(ks:ke), qr_crg_t(ks:ke), qi_crg_t(ks:ke), qs_crg_t(ks:ke), qg_crg_t(ks:ke)
1031  real(rp) :: rlambda(i_qc:i_qg)
1032  !---------------------------------------------------------------------------
1033 
1034  call prof_rapstart('MP_tomita08', 3)
1035 
1036  !$acc data copy(TEMP, QTRC, CPtot0, CVtot0) &
1037  !$acc copyin(DENS, PRES, CCN) &
1038  !$acc copyout(RHOE_t) &
1039  !$acc create(HIST_sw)
1040 
1041  rdens_i = 1.0_rp / dice
1042  rdens_g = 1.0_rp / dens_g
1043  rdens_r = 1.0_rp / dwatr
1044  rdens_s = 1.0_rp / dens_s
1045  if ( present(flg_lt) ) then
1046  flg_lt_l = flg_lt
1047  else
1048  flg_lt_l = .false.
1049  end if
1050  if( flg_lt_l ) then
1051 
1052  !$acc enter data copyin(dqcrg, beta_crg, &
1053  !$acc QTRC_crg0) &
1054  !$acc create(QSPLT_in, Sarea, &
1055  !$acc Re)
1056 
1057 !OCL ZFILL
1058  !$omp parallel workshare
1059  !$acc kernels
1060  qsplt_in(:,:,:,:) = 0.0_rp
1061  !$acc end kernels
1062  !$omp end parallel workshare
1063 
1064  call atmos_phy_mp_tomita08_effective_radius( &
1065  ka, ks, ke, ia, is, ie, ja, js, je, &
1066  dens, temp, qtrc, &
1067  re )
1068  else
1069 
1070  ! dummy
1071  !$acc enter data create(dqcrg, beta_crg, &
1072  !$acc QTRC_crg0) &
1073  !$acc create(QSPLT_in, Sarea, &
1074  !$acc Re)
1075 
1076  endif
1077 
1078  hist_flag = .false.
1079  do ip = 1, w_nmax
1080  call file_history_query( hist_id(ip), hist_sw(ip) )
1081  hist_flag = hist_flag .or. hist_sw(ip)
1082  end do
1083  !$acc update device(HIST_sw)
1084 
1085  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
1086  !$omp shared(KA,KS,KE,IS,IE,JS,JE, &
1087  !$omp DENS,TEMP,PRES,QTRC,CCN,CPtot0,CVtot0,dt, &
1088  !$omp RHOE_t, &
1089  !$omp UNDEF,EPS,PRE00,LHV,LHF,LHF0,CP_VAPOR,CP_WATER,CP_ICE,CV_VAPOR,CV_WATER,CV_ICE, &
1090  !$omp do_couple_aerosol,sw_expice,enable_WDXZ2014,enable_RS2014,enable_KK2000,enable_HZDFHI2007,enable_TRAWLBG2017, &
1091  !$omp Nc_def,N0r_def,N0s_def,N0g_def, &
1092  !$omp Cr,Cs,Cg,Erw,Eri,Eiw,Esw,Esr,Esi,Egw,Egr,Egi,Egs,Ar,As,Ag, &
1093  !$omp gamma_sacr,gamma_gacs,gamma_saut,gamma_gaut,beta_saut,beta_gaut,qicrt_saut,qscrt_gaut,mi, &
1094  !$omp GAM,GAM_2,GAM_3,GAM_1br,GAM_1bs,GAM_1bsds,GAM_1bg,GAM_1bgdg,GAM_1brdr,GAM_2br,GAM_2bs,GAM_3br,GAM_3bs,GAM_3dr,GAM_3ds,GAM_3dg,GAM_5dr_h,GAM_5ds_h,GAM_5dg_h,GAM_6dr, &
1095  !$omp QTRC_crg0,QSPLT_in,Sarea,Re,beta_crg,dqcrg, &
1096  !$omp re_qc,re_qi,d0_crg,v0_crg,rdens_i,Ecoal_Gi,Ecoal_GS, &
1097  !$omp flg_lt_l,flg_ecoali,flg_ecoals, &
1098  !$omp w3d,HIST_sw,hist_flag) &
1099  !$omp private(cvtot,qv,qc,qr,qi,qs,qg,qv_t,qc_t,qr_t,qi_t,qs_t,qg_t,e_t,cp_t,cv_t, &
1100  !$omp QSATL,QSATI,Sliq,Sice,rho_fact,temc,N0r,N0s,N0g, &
1101  !$omp RLMDr,RLMDr_2,RLMDr_3,RLMDs,RLMDs_2,RLMDs_3,RLMDg,RLMDg_2,RLMDg_3, &
1102  !$omp RLMDr_1br,RLMDr_2br,RLMDr_3br,RLMDs_1bs,RLMDs_2bs,RLMDs_3bs, &
1103  !$omp RLMDr_dr,RLMDr_3dr,RLMDr_5dr,RLMDs_ds,RLMDs_3ds,RLMDs_5ds, &
1104  !$omp RLMDg_dg,RLMDg_3dg,RLMDg_5dg,RLMDr_7,RLMDr_6dr, &
1105  !$omp tems,Xs2,MOMs_0,MOMs_1,MOMs_2,MOMs_0bs,MOMs_1bs,MOMs_2bs,MOMs_2ds,MOMs_5ds_h,RMOMs_Vt, &
1106  !$omp fact,temcc,itemc, &
1107  !$omp coef_at,coef_bt,loga_,b_,nm, &
1108  !$omp coef_a0,coef_a1,coef_b0,coef_b1, &
1109  !$omp Vti,Vtr,Vts,Vtg,Esi_mod,Egs_mod,rhoqc,Nc, &
1110  !$omp Pracw_orig,Pracw_kk,Praut_berry,Praut_kk,Dc,betai,betas,Da,Kd,Nu, &
1111  !$omp Glv,Giv,Gil,ventr,vents,ventg,net,fac,fack,fac_sw,zerosw,tmp, &
1112  !$omp qc_crg_t,qr_crg_t,qi_crg_t,qs_crg_t,qg_crg_t,qcrg_c,qcrg_r,qcrg_i,qcrg_s,qcrg_g, &
1113  !$omp w_q,w_qcrg,re_qs,dcrg,beta1_crg,alpha,facq_QC,facq_QR,facq_QI,facq_QS,facq_QG,rlambda, &
1114  !$omp sw_bergeron,a1,a2,ma2,dt1,Ni50, &
1115  !$omp sw,rhoqi,XNi,XMi,Di,Nig,Qig,w)
1116  !$acc kernels
1117  do j = js, je
1118  !$acc loop private(cvtot, qv, qc, qr, qi, qs, qg, qv_t, qc_t, qr_t, qi_t, qs_t, qg_t, &
1119  !$acc QSATL, QSATI, Sliq, Sice, rho_fact, temc, N0r, N0s, N0g, &
1120  !$acc RLMDr, RLMDr_2, RLMDr_3, RLMDg, RLMDg_2, RLMDg_3, &
1121  !$acc RLMDr_1br, RLMDr_2br, RLMDr_3br, RLMDr_dr, RLMDr_3dr, RLMDr_5dr, &
1122  !$acc RLMDg_dg, RLMDg_3dg, RLMDg_5dg, RLMDr_7, RLMDr_6dr, &
1123  !$acc MOMs_0, MOMs_1, MOMs_2, MOMs_0bs, MOMs_1bs, MOMs_2bs, MOMs_2ds, MOMs_5ds_h, RMOMs_Vt, &
1124  !$acc Vti, Vtr, Vts, Vtg, Egs_mod, Nc, Nu, fack, Glv, Giv, Gil, sw_bergeron, a1, a2, ma2, w, &
1125  !$acc qcrg_c, qcrg_r, qcrg_i, qcrg_s, qcrg_g, w_q, w_qcrg, &
1126  !$acc dcrg, beta1_crg, re_qs, qc_crg_t, qr_crg_t, qi_crg_t, qs_crg_t, qg_crg_t)
1127  do i = is, ie
1128 
1129  if ( do_couple_aerosol ) then
1130  do k = ks, ke
1131  nc(k) = max( ccn(k,i,j)*1.e-6_rp, nc_def(i,j) ) ! [#/m3]->[#/cc]
1132  end do
1133  else
1134  do k = ks, ke
1135  nc(k) = nc_def(i,j)
1136  end do
1137  endif
1138 
1139  ! store to work
1140  call saturation_dens2qsat_liq( ka, ks, ke, &
1141  temp(:,i,j), dens(:,i,j), & ! [IN]
1142  qsatl(:) ) ! [OUT]
1143 
1144  call saturation_dens2qsat_ice( ka, ks, ke, &
1145  temp(:,i,j), dens(:,i,j), & ! [IN]
1146  qsati(:) ) ! [OUT]
1147 
1148  ! store to work
1149  do k = ks, ke
1150  qv(k) = max( qtrc(k,i,j,i_qv), 0.0_rp )
1151  end do
1152  do k = ks, ke
1153  qc(k) = max( qtrc(k,i,j,i_qc), 0.0_rp )
1154  end do
1155  do k = ks, ke
1156  qr(k) = max( qtrc(k,i,j,i_qr), 0.0_rp )
1157  end do
1158  do k = ks, ke
1159  qi(k) = max( qtrc(k,i,j,i_qi), 0.0_rp )
1160  end do
1161  do k = ks, ke
1162  qs(k) = max( qtrc(k,i,j,i_qs), 0.0_rp )
1163  end do
1164  do k = ks, ke
1165  qg(k) = max( qtrc(k,i,j,i_qg), 0.0_rp )
1166  end do
1167 
1168  if ( flg_lt_l ) then
1169  ! store to work
1170  do k = ks, ke
1171  qcrg_c(k) = qtrc_crg0(k,i,j,i_qc-1)
1172  enddo
1173  do k = ks, ke
1174  qcrg_r(k) = qtrc_crg0(k,i,j,i_qr-1)
1175  enddo
1176  do k = ks, ke
1177  qcrg_i(k) = qtrc_crg0(k,i,j,i_qi-1)
1178  enddo
1179  do k = ks, ke
1180  qcrg_s(k) = qtrc_crg0(k,i,j,i_qs-1)
1181  enddo
1182  do k = ks, ke
1183  qcrg_g(k) = qtrc_crg0(k,i,j,i_qg-1)
1184  enddo
1185  do k = ks, ke
1186  re_qs(k) = re(k,i,j,i_hs) * 1.0e-2_rp ! [cm] -> [m]
1187  enddo
1188  do k = ks, ke
1189  beta1_crg(k) = beta_crg(k,i,j)
1190  dcrg(k) = - dqcrg(k,i,j)
1191  end do
1192  end if
1193 
1194  ! saturation ratio S
1195  do k = ks, ke
1196  sliq(k) = qv(k) / max( qsatl(k), eps )
1197  sice(k) = qv(k) / max( qsati(k), eps )
1198 
1199  rho_fact(k) = sqrt( dens00 / dens(k,i,j) )
1200  temc(k) = temp(k,i,j) - tem00
1201  end do
1202 
1203  do k = ks, ke
1204  w(k,i_delta1) = ( 0.5_rp + sign(0.5_rp, qr(k) - 1.e-4_rp ) )
1205 
1206  w(k,i_delta2) = ( 0.5_rp + sign(0.5_rp, 1.e-4_rp - qr(k) ) ) &
1207  * ( 0.5_rp + sign(0.5_rp, 1.e-4_rp - qs(k) ) )
1208 
1209  w(k,i_spsati) = 0.5_rp + sign(0.5_rp, sice(k) - 1.0_rp )
1210 
1211  w(k,i_iceflg) = 0.5_rp - sign( 0.5_rp, temc(k) ) ! 0: warm, 1: ice
1212  end do
1213 
1214  do k = ks, ke
1215  w(k,i_dqv_dt) = qv(k) / dt
1216  w(k,i_dqc_dt) = qc(k) / dt
1217  w(k,i_dqr_dt) = qr(k) / dt
1218  w(k,i_dqi_dt) = qi(k) / dt
1219  w(k,i_dqs_dt) = qs(k) / dt
1220  w(k,i_dqg_dt) = qg(k) / dt
1221  end do
1222 
1223  do k = ks, ke
1224  sw_bergeron(k) = ( 0.5_rp + sign(0.5_rp, temc(k) + 30.0_rp ) ) &
1225  * ( 0.5_rp + sign(0.5_rp, 0.0_rp - temc(k) ) ) &
1226  * ( 1.0_rp - sw_expice )
1227  end do
1228 
1229  ! intercept parameter N0
1230  if ( enable_wdxz2014 ) then ! Wainwright et al. (2014)
1231  do k = ks, ke
1232  n0r(k) = 1.16e+5_rp * exp( log( max( dens(k,i,j)*qr(k)*1000.0_rp, 1.e-2_rp ) )*0.477_rp )
1233  n0s(k) = 4.58e+9_rp * exp( log( max( dens(k,i,j)*qs(k)*1000.0_rp, 1.e-2_rp ) )*0.788_rp )
1234  n0g(k) = 9.74e+8_rp * exp( log( max( dens(k,i,j)*qg(k)*1000.0_rp, 1.e-2_rp ) )*0.816_rp )
1235  end do
1236  else
1237  do k = ks, ke
1238  n0r(k) = n0r_def ! Marshall and Palmer (1948)
1239  n0s(k) = n0s_def ! Gunn and Marshall (1958)
1240  n0g(k) = n0g_def
1241  end do
1242  end if
1243 
1244  ! slope parameter lambda (Rain)
1245  do k = ks, ke
1246  zerosw = 0.5_rp - sign(0.5_rp, qr(k) - 1.e-12_rp )
1247  rlmdr(k) = sqrt(sqrt( dens(k,i,j) * qr(k) / ( ar * n0r(k) * gam_1br ) + zerosw )) * ( 1.0_rp-zerosw )
1248 
1249  rlmdr_dr(k) = sqrt( rlmdr(k) ) ! **Dr
1250  rlmdr_2(k) = rlmdr(k)**2
1251  rlmdr_3(k) = rlmdr(k)**3
1252  rlmdr_7(k) = rlmdr(k)**7
1253  rlmdr_1br(k) = rlmdr(k)**4 ! (1+Br)
1254  rlmdr_2br(k) = rlmdr(k)**5 ! (2+Br)
1255  rlmdr_3br(k) = rlmdr(k)**6 ! (3+Br)
1256  rlmdr_3dr(k) = rlmdr(k)**3 * rlmdr_dr(k)
1257  rlmdr_5dr(k) = rlmdr(k)**5 * rlmdr_dr(k)
1258  rlmdr_6dr(k) = rlmdr(k)**6 * rlmdr_dr(k)
1259 
1260  w(k,i_rlmdr) = rlmdr(k)
1261  end do
1262 
1263  ! slope parameter lambda (Snow)
1264  if ( enable_rs2014 ) then
1265  !$acc loop private(coef_at, coef_bt)
1266  do k = ks, ke
1267  !---< modification by Roh and Satoh (2014) >---
1268  ! bimodal size distribution of snow
1269  zerosw = 0.5_rp - sign(0.5_rp, dens(k,i,j) * qs(k) - 1.e-12_rp )
1270 
1271  xs2 = dens(k,i,j) * qs(k) / as
1272  tems = min( -0.1_rp, temc(k) )
1273  coef_at(1) = coef_a01 + tems * ( coef_a02 + tems * ( coef_a05 + tems * coef_a09 ) )
1274  coef_at(2) = coef_a03 + tems * ( coef_a04 + tems * coef_a07 )
1275  coef_at(3) = coef_a06 + tems * coef_a08
1276  coef_at(4) = coef_a10
1277  coef_bt(1) = coef_b01 + tems * ( coef_b02 + tems * ( coef_b05 + tems * coef_b09 ) )
1278  coef_bt(2) = coef_b03 + tems * ( coef_b04 + tems * coef_b07 )
1279  coef_bt(3) = coef_b06 + tems * coef_b08
1280  coef_bt(4) = coef_b10
1281  ! 0th moment
1282  loga_ = coef_at(1)
1283  b_ = coef_bt(1)
1284  moms_0(k) = exp(ln10*loga_) * exp(log(xs2+zerosw)*b_) * ( 1.0_rp-zerosw )
1285  ! 1st moment
1286  nm = 1.0_rp
1287  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
1288  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
1289  moms_1(k) = exp(ln10*loga_) * exp(log(xs2+zerosw)*b_) * ( 1.0_rp-zerosw )
1290  ! 2nd moment
1291  moms_2(k) = xs2
1292  ! 0 + Bs(=2) moment
1293  nm = 2.0_rp
1294  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
1295  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
1296  moms_0bs(k) = exp(ln10*loga_) * exp(log(xs2+zerosw)*b_) * ( 1.0_rp-zerosw )
1297  ! 1 + Bs(=2) moment
1298  nm = 3.0_rp
1299  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
1300  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
1301  moms_1bs(k) = exp(ln10*loga_) * exp(log(xs2+zerosw)*b_) * ( 1.0_rp-zerosw )
1302  ! 2 + Bs(=2) moment
1303  nm = 4.0_rp
1304  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
1305  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
1306  moms_2bs(k) = exp(ln10*loga_) * exp(log(xs2+zerosw)*b_) * ( 1.0_rp-zerosw )
1307  ! 2 + Ds(=0.25) moment
1308  nm = 2.25_rp
1309  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
1310  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
1311  moms_2ds(k) = exp(ln10*loga_) * exp(log(xs2+zerosw)*b_) * ( 1.0_rp-zerosw )
1312  ! ( 3 + Ds(=0.25) ) / 2 moment
1313  nm = 1.625_rp
1314  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
1315  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
1316  moms_5ds_h(k) = exp(ln10*loga_) * exp(log(xs2+zerosw)*b_) * ( 1.0_rp-zerosw )
1317  ! Bs(=2) + Ds(=0.25) moment
1318  nm = 2.25_rp
1319  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
1320  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
1321  rmoms_vt(k) = exp(ln10*loga_) * exp(log(xs2+zerosw)*b_) * ( 1.0_rp-zerosw ) / ( moms_0bs(k) + zerosw )
1322 
1323  w(k,i_rlmds) = undef
1324  end do
1325  else
1326  do k = ks, ke
1327  zerosw = 0.5_rp - sign(0.5_rp, dens(k,i,j) * qs(k) - 1.e-12_rp )
1328 
1329  rlmds = sqrt(sqrt( dens(k,i,j) * qs(k) / ( as * n0s(k) * gam_1bs ) + zerosw )) * ( 1.0_rp-zerosw )
1330  rlmds_ds = sqrt( sqrt(rlmds) ) ! **Ds
1331  rlmds_2 = rlmds**2
1332  rlmds_3 = rlmds_2 * rlmds
1333  rlmds_1bs = rlmds_2 * rlmds_2 ! (1+Bs)
1334  rlmds_2bs = rlmds_3 * rlmds_2 ! (2+Bs)
1335  rlmds_3bs = rlmds_3 * rlmds_3 ! (3+Bs)
1336  rlmds_3ds = rlmds_3 * rlmds_ds
1337  rlmds_5ds = rlmds_2 * rlmds_3ds
1338 
1339  moms_0(k) = n0s(k) * gam * rlmds ! Ns * 0th moment
1340  moms_1(k) = n0s(k) * gam_2 * rlmds_2 ! Ns * 1st moment
1341  moms_2(k) = n0s(k) * gam_3 * rlmds_3 ! Ns * 2nd moment
1342  moms_0bs(k) = n0s(k) * gam_1bs * rlmds_1bs ! Ns * 0+bs
1343  moms_1bs(k) = n0s(k) * gam_2bs * rlmds_2bs ! Ns * 1+bs
1344  moms_2bs(k) = n0s(k) * gam_3bs * rlmds_3bs ! Ns * 2+bs
1345  moms_2ds(k) = n0s(k) * gam_3ds * rlmds_3ds ! Ns * 2+ds
1346  moms_5ds_h(k) = n0s(k) * gam_5ds_h * sqrt(rlmds_5ds) ! Ns * (5+ds)/2
1347  rmoms_vt(k) = gam_1bsds / gam_1bs * rlmds_ds
1348 
1349  w(k,i_rlmds) = rlmds
1350  end do
1351  end if
1352 
1353  do k = ks, ke
1354  ! slope parameter lambda (Graupel)
1355  zerosw = 0.5_rp - sign(0.5_rp, qg(k) - 1.e-12_rp )
1356  rlmdg(k) = sqrt(sqrt( dens(k,i,j) * qg(k) / ( ag * n0g(k) * gam_1bg ) + zerosw )) * ( 1.0_rp-zerosw )
1357 
1358  rlmdg_dg(k) = sqrt( rlmdg(k) ) ! **Dg
1359  rlmdg_2(k) = rlmdg(k)**2
1360  rlmdg_3(k) = rlmdg(k) * rlmdg_2(k)
1361  rlmdg_3dg(k) = rlmdg_3(k) * rlmdg_dg(k)
1362  rlmdg_5dg(k) = rlmdg_2(k) * rlmdg_3dg(k)
1363 
1364  w(k,i_rlmdg) = rlmdg(k)
1365  end do
1366 
1367  !---< terminal velocity >
1368  if ( enable_hzdfhi2007 ) then
1369  do k = ks, ke
1370  zerosw = 0.5_rp - sign(0.5_rp, qi(k) - 1.e-8_rp )
1371  vti(k) = min( 0.0_rp, &
1372  - ( coef_a0 + coef_a1 * temc(k) ) &
1373  * exp( log( dens(k,i,j)*qi(k)*1000.0_rp+zerosw ) * ( coef_b0 + coef_b1 * temc(k) ) ) &
1374  * 1e-2_rp * ( 1.0_rp-zerosw ) )
1375  end do
1376  else
1377  do k = ks, ke
1378  zerosw = 0.5_rp - sign(0.5_rp, qi(k) - 1.e-8_rp )
1379  vti(k) = -3.29_rp * exp( log( dens(k,i,j)*qi(k)+zerosw )*0.16_rp ) * ( 1.0_rp-zerosw )
1380  end do
1381  end if
1382 
1383  do k = ks, ke
1384  vtr(k) = -cr * rho_fact(k) * gam_1brdr / gam_1br * rlmdr_dr(k)
1385  vts(k) = -cs * rho_fact(k) * rmoms_vt(k)
1386  vtg(k) = -cg * rho_fact(k) * gam_1bgdg / gam_1bg * rlmdg_dg(k)
1387  end do
1388 
1389 
1390  !---< Nucleation >---
1391  ! [Pigen] ice nucleation
1392  do k = ks, ke
1393  nig = max( exp(-0.1_rp*temc(k)), 1.0_rp ) * 1000.0_rp
1394  qig = 4.92e-11_rp * exp(log(nig)*1.33_rp) / dens(k,i,j)
1395  w(k,i_pigen) = max( min( qig-qi(k), qv(k)-qsati(k) ), 0.0_rp ) / dt
1396  end do
1397 
1398 
1399  !---< Accretion >---
1400 
1401  ! [Pracw] accretion rate of cloud water by rain
1402  if ( enable_kk2000 ) then
1403  do k = ks, ke
1404  zerosw = 0.5_rp - sign(0.5_rp, qc(k)*qr(k) - 1.e-12_rp )
1405  pracw_kk = 67.0_rp * exp( log( qc(k)*qr(k)+zerosw )*1.15_rp ) * ( 1.0_rp-zerosw ) ! eq.(33) in KK(2000)
1406  w(k,i_pracw) = pracw_kk
1407  end do
1408  else
1409  do k = ks, ke
1410  pracw_orig = qc(k) * 0.25_rp * pi * erw * n0r(k) * cr * gam_3dr * rlmdr_3dr(k) * rho_fact(k)
1411  w(k,i_pracw) = pracw_orig
1412  end do
1413  end if
1414 
1415  do k = ks, ke
1416  ! [Psacw] accretion rate of cloud water by snow
1417  w(k,i_psacw) = qc(k) * 0.25_rp * pi * esw * cs * moms_2ds(k) * rho_fact(k)
1418  ! [Pgacw] accretion rate of cloud water by graupel
1419  w(k,i_pgacw) = qc(k) * 0.25_rp * pi * egw * n0g(k) * cg * gam_3dg * rlmdg_3dg(k) * rho_fact(k)
1420  end do
1421 
1422  do k = ks, ke
1423  esi_mod = min( esi, esi * exp( gamma_sacr * temc(k) ) )
1424  ! [Praci] accretion rate of cloud ice by rain
1425  w(k,i_praci) = qi(k) * 0.25_rp * pi * eri * n0r(k) * cr * gam_3dr * rlmdr_3dr(k) * rho_fact(k)
1426  ! [Psaci] accretion rate of cloud ice by snow
1427  w(k,i_psaci) = qi(k) * 0.25_rp * pi * esi_mod * cs * moms_2ds(k) * rho_fact(k)
1428  ! [Pgaci] accretion rate of cloud ice by grupel
1429  w(k,i_pgaci) = qi(k) * 0.25_rp * pi * egi * n0g(k) * cg * gam_3dg * rlmdg_3dg(k) * rho_fact(k)
1430  ! [Piacr] accretion rate of rain by cloud ice
1431  w(k,i_piacr) = qi(k) * ar / mi * 0.25_rp * pi * eri * n0r(k) * cr * gam_6dr * rlmdr_6dr(k) * rho_fact(k)
1432  end do
1433 
1434  do k = ks, ke
1435  ! [Psacr] accretion rate of rain by snow
1436  w(k,i_psacr) = ar * 0.25_rp * pi / dens(k,i,j) * esr * n0r(k) * abs(vtr(k)-vts(k)) &
1437  * ( gam_1br * rlmdr_1br(k) * moms_2(k) &
1438  + 2.0_rp * gam_2br * rlmdr_2br(k) * moms_1(k) &
1439  + gam_3br * rlmdr_3br(k) * moms_0(k) )
1440  end do
1441 
1442  do k = ks, ke
1443  ! [Pgacr] accretion rate of rain by graupel
1444  w(k,i_pgacr) = ar * 0.25_rp * pi / dens(k,i,j) * egr * n0g(k) * n0r(k) * abs(vtg(k)-vtr(k)) &
1445  * ( gam_1br * rlmdr_1br(k) * gam_3 * rlmdg_3(k) &
1446  + 2.0_rp * gam_2br * rlmdr_2br(k) * gam_2 * rlmdg_2(k) &
1447  + gam_3br * rlmdr_3br(k) * gam * rlmdg(k) )
1448  end do
1449 
1450  do k = ks, ke
1451  ! [Pracs] accretion rate of snow by rain
1452  w(k,i_pracs) = as * 0.25_rp * pi / dens(k,i,j) * esr * n0r(k) * abs(vtr(k)-vts(k)) &
1453  * ( moms_0bs(k) * gam_3 * rlmdr_3(k) &
1454  + 2.0_rp * moms_1bs(k) * gam_2 * rlmdr_2(k) &
1455  + moms_2bs(k) * gam * rlmdr(k) )
1456  end do
1457 
1458  do k = ks, ke
1459  ! [Pgacs] accretion rate of snow by graupel
1460  egs_mod(k) = min( egs, egs * exp( gamma_gacs * temc(k) ) )
1461  w(k,i_pgacs) = as * 0.25_rp * pi / dens(k,i,j) * egs_mod(k) * n0g(k) * abs(vtg(k)-vts(k)) &
1462  * ( moms_0bs(k) * gam_3 * rlmdg_3(k) &
1463  + 2.0_rp * moms_1bs(k) * gam_2 * rlmdg_2(k) &
1464  + moms_2bs(k) * gam * rlmdg(k) )
1465 
1466  end do
1467 
1468 
1469  !---< Auto-conversion >---
1470 
1471  ! [Praut] auto-conversion rate from cloud water to rain
1472  if ( enable_kk2000 ) then
1473  do k = ks, ke
1474  zerosw = 0.5_rp - sign(0.5_rp, qc(k) - 1.e-12_rp )
1475  praut_kk = 1350.0_rp &
1476  * exp( log( qc(k)+zerosw )*2.47_rp + log( nc(k) )*(-1.79_rp) ) &
1477  * ( 1.0_rp-zerosw ) ! eq.(29) in KK(2000)
1478  w(k,i_praut) = praut_kk
1479  end do
1480  else
1481  do k = ks, ke
1482  rhoqc = dens(k,i,j) * qc(k) * 1000.0_rp ! [g/m3]
1483  dc = 0.146_rp - 5.964e-2_rp * log( nc(k) / 2000.0_rp )
1484  praut_berry = 1.67e-5_rp * rhoqc * rhoqc / ( ( 5.0_rp + 3.66e-2_rp * nc(k) / ( dc * rhoqc + eps ) ) * dens(k,i,j) )
1485  w(k,i_praut) = praut_berry
1486  end do
1487  end if
1488 
1489  do k = ks, ke
1490  ! [Psaut] auto-conversion rate from cloud ice to snow
1491  betai = min( beta_saut, beta_saut * exp( gamma_saut * temc(k) ) )
1492  w(k,i_psaut) = max( betai*(qi(k)-qicrt_saut), 0.0_rp )
1493  ! [Pgaut] auto-conversion rate from snow to graupel
1494  betas = min( beta_gaut, beta_gaut * exp( gamma_gaut * temc(k) ) )
1495  w(k,i_pgaut) = max( betas*(qs(k)-qscrt_gaut), 0.0_rp )
1496  end do
1497 
1498 
1499  !---< Evaporation, Sublimation, Melting, and Freezing >---
1500 
1501  do k = ks, ke
1502  da = ( da0 + dda_dt * temc(k) )
1503  kd = ( dw0 + ddw_dt * temc(k) ) * pre00 / pres(k,i,j)
1504  nu(k) = ( mu0 + dmu_dt * temc(k) ) / dens(k,i,j)
1505 
1506  glv(k) = 1.0_rp / ( lhv0/(da*temp(k,i,j)) * ( lhv0/(rvap*temp(k,i,j)) - 1.0_rp ) + 1.0_rp/(kd*dens(k,i,j)*qsatl(k)) )
1507  giv(k) = 1.0_rp / ( lhs0/(da*temp(k,i,j)) * ( lhs0/(rvap*temp(k,i,j)) - 1.0_rp ) + 1.0_rp/(kd*dens(k,i,j)*qsati(k)) )
1508  gil(k) = ( da * temc(k) ) / lhf0
1509  end do
1510 
1511  do k = ks, ke
1512  ! [Prevp] evaporation rate of rain
1513  ventr = f1r * gam_2 * rlmdr_2(k) + f2r * sqrt( cr * rho_fact(k) / nu(k) * rlmdr_5dr(k) ) * gam_5dr_h
1514  w(k,i_prevp) = 2.0_rp * pi / dens(k,i,j) * n0r(k) * ( 1.0_rp-min(sliq(k),1.0_rp) ) * glv(k) * ventr
1515  end do
1516 
1517  do k = ks, ke
1518  ! [Pidep,Pisub] deposition/sublimation rate for ice
1519  rhoqi = max(dens(k,i,j)*qi(k), eps)
1520  xni = min( max( 5.38e+7_rp * exp( log(rhoqi)*0.75_rp ), 1.e+3_rp ), 1.e+6_rp )
1521  xmi = rhoqi / xni
1522  di = min( di_a * sqrt(xmi), di_max )
1523  tmp = 4.0_rp * di * xni / dens(k,i,j) * ( sice(k)-1.0_rp ) * giv(k)
1524  w(k,i_pidep) = ( w(k,i_spsati) ) * ( tmp) ! Sice > 1
1525  w(k,i_pisub) = ( 1.0_rp-w(k,i_spsati) ) * (-tmp) ! Sice < 1
1526 
1527  ! [Pihom] homogenious freezing at T < -40C
1528  sw = ( 0.5_rp - sign(0.5_rp, temc(k) + 40.0_rp ) ) ! if T < -40C, sw=1
1529  w(k,i_pihom) = sw * qc(k) / dt
1530 
1531  ! [Pihtr] heteroginous freezing at -40C < T < 0C
1532  sw = ( 0.5_rp + sign(0.5_rp, temc(k) + 40.0_rp ) ) &
1533  * ( 0.5_rp - sign(0.5_rp, temc(k) ) ) ! if -40C < T < 0C, sw=1
1534  w(k,i_pihtr) = sw * ( dens(k,i,j) / dwatr * qc(k)**2 / ( nc_ihtr * 1.e+6_rp ) ) &
1535  * b_frz * ( exp(-a_frz*temc(k)) - 1.0_rp )
1536 
1537  ! [Pimlt] ice melting at T > 0C
1538  sw = ( 0.5_rp + sign(0.5_rp, temc(k) ) ) ! if T > 0C, sw=1
1539  w(k,i_pimlt) = sw * qi(k) / dt
1540  end do
1541 
1542  do k = ks, ke
1543  ! [Psdep,Pssub] deposition/sublimation rate for snow
1544  vents = f1s * moms_1(k) + f2s * sqrt( cs * rho_fact(k) / nu(k) ) * moms_5ds_h(k)
1545  tmp = 2.0_rp * pi / dens(k,i,j) * ( sice(k)-1.0_rp ) * giv(k) * vents
1546  w(k,i_psdep) = ( w(k,i_spsati) ) * ( tmp) ! Sice > 1
1547  w(k,i_pssub) = ( 1.0_rp-w(k,i_spsati) ) * (-tmp) ! Sice < 1
1548  ! [Psmlt] melting rate of snow
1549  w(k,i_psmlt) = 2.0_rp * pi / dens(k,i,j) * gil(k) * vents &
1550  + cl * temc(k) / lhf0 * ( w(k,i_psacw) + w(k,i_psacr) )
1551  w(k,i_psmlt) = max( w(k,i_psmlt), 0.0_rp )
1552  end do
1553 
1554  do k = ks, ke
1555  ! [Pgdep/pgsub] deposition/sublimation rate for graupel
1556  ventg = f1g * gam_2 * rlmdg_2(k) + f2g * sqrt( cg * rho_fact(k) / nu(k) * rlmdg_5dg(k) ) * gam_5dg_h
1557  tmp = 2.0_rp * pi / dens(k,i,j) * n0g(k) * ( sice(k)-1.0_rp ) * giv(k) * ventg
1558  w(k,i_pgdep) = ( w(k,i_spsati) ) * ( tmp) ! Sice > 1
1559  w(k,i_pgsub) = ( 1.0_rp-w(k,i_spsati) ) * (-tmp) ! Sice < 1
1560  ! [Pgmlt] melting rate of graupel
1561  w(k,i_pgmlt) = 2.0_rp * pi / dens(k,i,j) * n0g(k) * gil(k) * ventg &
1562  + cl * temc(k) / lhf0 * ( w(k,i_pgacw) + w(k,i_pgacr) )
1563  w(k,i_pgmlt) = max( w(k,i_pgmlt), 0.0_rp )
1564  end do
1565 
1566  ! [Pgfrz] freezing rate of graupel
1567  do k = ks, ke
1568  tmp = ( exp(-a_frz*temc(k)) - 1.0_rp ) * rlmdr_7(k) ! to avoid floating overflow
1569  w(k,i_pgfrz) = 2.0_rp * pi / dens(k,i,j) * n0r(k) * 60.0_rp * b_frz * ar * tmp
1570  end do
1571 
1572  do k = ks, ke
1573  ! [Psfw,Psfi] ( Bergeron process ) growth rate of snow by Bergeron process from cloud water/ice
1574  temcc = min( max( temc(k), -30.99_rp ), 0.0_rp ) ! 0C <= T < 31C
1575  itemc = int( -temcc ) + 1 ! 1 <= iT <= 31
1576  fact = - ( temcc + real(itemc-1,kind=8) )
1577  a1 = ( 1.0_rp-fact ) * bergeron_a1_tab(itemc ) &
1578  + ( fact ) * bergeron_a1_tab(itemc+1)
1579  a2 = ( 1.0_rp-fact ) * bergeron_a2_tab(itemc ) &
1580  + ( fact ) * bergeron_a2_tab(itemc+1)
1581  ma2 = 1.0_rp - a2
1582  a1 = a1 * exp( log(1.e-3_rp)*ma2 ) ! [g->kg]
1583  dt1 = ( exp( log(mi50)*ma2 ) &
1584  - exp( log(mi40)*ma2 ) ) / ( a1 * ma2 )
1585  ni50 = qi(k) * dt / ( mi50 * dt1 )
1586  w(k,i_psfw ) = ni50 * ( a1 * exp( log(mi50)*a2 ) &
1587  + pi * eiw * dens(k,i,j) * qc(k) * ri50*ri50 * vti50 )
1588  w(k,i_psfi ) = qi(k) / dt1
1589  end do
1590 
1591  do k = ks, ke
1592  !---< limiter >---
1593  w(k,i_pigen) = min( w(k,i_pigen), w(k,i_dqv_dt) ) * ( w(k,i_iceflg) ) * sw_expice
1594  w(k,i_pidep) = min( w(k,i_pidep), w(k,i_dqv_dt) ) * ( w(k,i_iceflg) ) * sw_expice
1595  w(k,i_psdep) = min( w(k,i_psdep), w(k,i_dqv_dt) ) * ( w(k,i_iceflg) )
1596  w(k,i_pgdep) = min( w(k,i_pgdep), w(k,i_dqv_dt) ) * ( w(k,i_iceflg) )
1597 
1598  w(k,i_pracw) = w(k,i_pracw) &
1599  + w(k,i_psacw) * ( 1.0_rp-w(k,i_iceflg) ) & ! c->r by s
1600  + w(k,i_pgacw) * ( 1.0_rp-w(k,i_iceflg) ) ! c->r by g
1601  end do
1602 
1603  do k = ks, ke
1604  w(k,i_praut) = min( w(k,i_praut), w(k,i_dqc_dt) )
1605  w(k,i_pracw) = min( w(k,i_pracw), w(k,i_dqc_dt) )
1606  w(k,i_pihom) = min( w(k,i_pihom), w(k,i_dqc_dt) ) * ( w(k,i_iceflg) ) * sw_expice
1607  w(k,i_pihtr) = min( w(k,i_pihtr), w(k,i_dqc_dt) ) * ( w(k,i_iceflg) ) * sw_expice
1608  w(k,i_psacw) = min( w(k,i_psacw), w(k,i_dqc_dt) ) * ( w(k,i_iceflg) )
1609  w(k,i_psfw ) = min( w(k,i_psfw ), w(k,i_dqc_dt) ) * ( w(k,i_iceflg) ) * sw_bergeron(k)
1610  w(k,i_pgacw) = min( w(k,i_pgacw), w(k,i_dqc_dt) ) * ( w(k,i_iceflg) )
1611  end do
1612 
1613  do k = ks, ke
1614  w(k,i_prevp) = min( w(k,i_prevp), w(k,i_dqr_dt) )
1615  w(k,i_piacr) = min( w(k,i_piacr), w(k,i_dqr_dt) ) * ( w(k,i_iceflg) )
1616  w(k,i_psacr) = min( w(k,i_psacr), w(k,i_dqr_dt) ) * ( w(k,i_iceflg) )
1617  w(k,i_pgacr) = min( w(k,i_pgacr), w(k,i_dqr_dt) ) * ( w(k,i_iceflg) )
1618  w(k,i_pgfrz) = min( w(k,i_pgfrz), w(k,i_dqr_dt) ) * ( w(k,i_iceflg) )
1619  end do
1620 
1621  do k = ks, ke
1622  w(k,i_pisub) = min( w(k,i_pisub), w(k,i_dqi_dt) ) * ( w(k,i_iceflg) ) * sw_expice
1623  w(k,i_pimlt) = min( w(k,i_pimlt), w(k,i_dqi_dt) ) * ( 1.0_rp-w(k,i_iceflg) ) * sw_expice
1624  w(k,i_psaut) = min( w(k,i_psaut), w(k,i_dqi_dt) ) * ( w(k,i_iceflg) )
1625  w(k,i_praci) = min( w(k,i_praci), w(k,i_dqi_dt) ) * ( w(k,i_iceflg) )
1626  w(k,i_psaci) = min( w(k,i_psaci), w(k,i_dqi_dt) ) * ( w(k,i_iceflg) )
1627  w(k,i_psfi ) = min( w(k,i_psfi ), w(k,i_dqi_dt) ) * ( w(k,i_iceflg) ) * sw_bergeron(k)
1628  w(k,i_pgaci) = min( w(k,i_pgaci), w(k,i_dqi_dt) ) * ( w(k,i_iceflg) )
1629  end do
1630 
1631  do k = ks, ke
1632  w(k,i_pssub) = min( w(k,i_pssub), w(k,i_dqs_dt) ) * ( w(k,i_iceflg) )
1633  w(k,i_psmlt) = min( w(k,i_psmlt), w(k,i_dqs_dt) ) * ( 1.0_rp-w(k,i_iceflg) )
1634  w(k,i_pgaut) = min( w(k,i_pgaut), w(k,i_dqs_dt) ) * ( w(k,i_iceflg) )
1635  w(k,i_pracs) = min( w(k,i_pracs), w(k,i_dqs_dt) ) * ( w(k,i_iceflg) )
1636  w(k,i_pgacs) = min( w(k,i_pgacs), w(k,i_dqs_dt) )
1637 
1638  w(k,i_pgsub) = min( w(k,i_pgsub), w(k,i_dqg_dt) ) * ( w(k,i_iceflg) )
1639  w(k,i_pgmlt) = min( w(k,i_pgmlt), w(k,i_dqg_dt) ) * ( 1.0_rp-w(k,i_iceflg) )
1640  end do
1641 
1642  do k = ks, ke
1643  w(k,i_piacr_s) = ( 1.0_rp - w(k,i_delta1) ) * w(k,i_piacr)
1644  w(k,i_piacr_g) = ( w(k,i_delta1) ) * w(k,i_piacr)
1645  w(k,i_praci_s) = ( 1.0_rp - w(k,i_delta1) ) * w(k,i_praci)
1646  w(k,i_praci_g) = ( w(k,i_delta1) ) * w(k,i_praci)
1647  w(k,i_psacr_s) = ( w(k,i_delta2) ) * w(k,i_psacr)
1648  w(k,i_psacr_g) = ( 1.0_rp - w(k,i_delta2) ) * w(k,i_psacr)
1649  w(k,i_pracs ) = ( 1.0_rp - w(k,i_delta2) ) * w(k,i_pracs)
1650  end do
1651 
1652  ! [QC]
1653  do k = ks, ke
1654  net = &
1655  + w(k,i_pimlt ) & ! [prod] i->c
1656  - w(k,i_praut ) & ! [loss] c->r
1657  - w(k,i_pracw ) & ! [loss] c->r
1658  - w(k,i_pihom ) & ! [loss] c->i
1659  - w(k,i_pihtr ) & ! [loss] c->i
1660  - w(k,i_psacw ) & ! [loss] c->s
1661  - w(k,i_psfw ) & ! [loss] c->s
1662  - w(k,i_pgacw ) ! [loss] c->g
1663 
1664  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1665  fac = ( fac_sw ) &
1666  + ( 1.0_rp - fac_sw ) * min( -w(k,i_dqc_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1667 
1668  w(k,i_pimlt ) = w(k,i_pimlt ) * fac
1669  w(k,i_praut ) = w(k,i_praut ) * fac
1670  w(k,i_pracw ) = w(k,i_pracw ) * fac
1671  w(k,i_pihom ) = w(k,i_pihom ) * fac
1672  w(k,i_pihtr ) = w(k,i_pihtr ) * fac
1673  w(k,i_psacw ) = w(k,i_psacw ) * fac
1674  w(k,i_psfw ) = w(k,i_psfw ) * fac
1675  w(k,i_pgacw ) = w(k,i_pgacw ) * fac
1676  end do
1677 
1678  ! [QI]
1679  do k = ks, ke
1680  net = &
1681  + w(k,i_pigen ) & ! [prod] v->i
1682  + w(k,i_pidep ) & ! [prod] v->i
1683  + w(k,i_pihom ) & ! [prod] c->i
1684  + w(k,i_pihtr ) & ! [prod] c->i
1685  - w(k,i_pisub ) & ! [loss] i->v
1686  - w(k,i_pimlt ) & ! [loss] i->c
1687  - w(k,i_psaut ) & ! [loss] i->s
1688  - w(k,i_praci_s) & ! [loss] i->s
1689  - w(k,i_psaci ) & ! [loss] i->s
1690  - w(k,i_psfi ) & ! [loss] i->s
1691  - w(k,i_praci_g) & ! [loss] i->g
1692  - w(k,i_pgaci ) ! [loss] i->g
1693 
1694  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1695  fac = ( fac_sw ) &
1696  + ( 1.0_rp - fac_sw ) * min( -w(k,i_dqi_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1697 
1698  w(k,i_pigen ) = w(k,i_pigen ) * fac
1699  w(k,i_pidep ) = w(k,i_pidep ) * fac
1700  w(k,i_pihom ) = w(k,i_pihom ) * fac
1701  w(k,i_pihtr ) = w(k,i_pihtr ) * fac
1702  w(k,i_pisub ) = w(k,i_pisub ) * fac
1703  w(k,i_pimlt ) = w(k,i_pimlt ) * fac
1704  w(k,i_psaut ) = w(k,i_psaut ) * fac
1705  w(k,i_praci_s) = w(k,i_praci_s) * fac
1706  w(k,i_psaci ) = w(k,i_psaci ) * fac
1707  w(k,i_psfi ) = w(k,i_psfi ) * fac
1708  w(k,i_praci_g) = w(k,i_praci_g) * fac
1709  w(k,i_pgaci ) = w(k,i_pgaci ) * fac
1710  end do
1711 
1712  ! [QR]
1713  do k = ks, ke
1714  net = &
1715  + w(k,i_praut ) & ! [prod] c->r
1716  + w(k,i_pracw ) & ! [prod] c->r
1717  + w(k,i_psmlt ) & ! [prod] s->r
1718  + w(k,i_pgmlt ) & ! [prod] g->r
1719  - w(k,i_prevp ) & ! [loss] r->v
1720  - w(k,i_piacr_s) & ! [loss] r->s
1721  - w(k,i_psacr_s) & ! [loss] r->s
1722  - w(k,i_piacr_g) & ! [loss] r->g
1723  - w(k,i_psacr_g) & ! [loss] r->g
1724  - w(k,i_pgacr ) & ! [loss] r->g
1725  - w(k,i_pgfrz ) ! [loss] r->g
1726 
1727  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1728  fac = ( fac_sw ) &
1729  + ( 1.0_rp - fac_sw ) * min( -w(k,i_dqr_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1730 
1731  w(k,i_praut ) = w(k,i_praut ) * fac
1732  w(k,i_pracw ) = w(k,i_pracw ) * fac
1733  w(k,i_psmlt ) = w(k,i_psmlt ) * fac
1734  w(k,i_pgmlt ) = w(k,i_pgmlt ) * fac
1735  w(k,i_prevp ) = w(k,i_prevp ) * fac
1736  w(k,i_piacr_s) = w(k,i_piacr_s) * fac
1737  w(k,i_psacr_s) = w(k,i_psacr_s) * fac
1738  w(k,i_piacr_g) = w(k,i_piacr_g) * fac
1739  w(k,i_psacr_g) = w(k,i_psacr_g) * fac
1740  w(k,i_pgacr ) = w(k,i_pgacr ) * fac
1741  w(k,i_pgfrz ) = w(k,i_pgfrz ) * fac
1742  end do
1743 
1744  ! [QV]
1745  do k = ks, ke
1746  net = &
1747  + w(k,i_prevp ) & ! [prod] r->v
1748  + w(k,i_pisub ) & ! [prod] i->v
1749  + w(k,i_pssub ) & ! [prod] s->v
1750  + w(k,i_pgsub ) & ! [prod] g->v
1751  - w(k,i_pigen ) & ! [loss] v->i
1752  - w(k,i_pidep ) & ! [loss] v->i
1753  - w(k,i_psdep ) & ! [loss] v->s
1754  - w(k,i_pgdep ) ! [loss] v->g
1755 
1756  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1757  fac = ( fac_sw ) &
1758  + ( 1.0_rp - fac_sw ) * min( -w(k,i_dqv_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1759 
1760  w(k,i_prevp ) = w(k,i_prevp ) * fac
1761  w(k,i_pisub ) = w(k,i_pisub ) * fac
1762  w(k,i_pssub ) = w(k,i_pssub ) * fac
1763  w(k,i_pgsub ) = w(k,i_pgsub ) * fac
1764  w(k,i_pigen ) = w(k,i_pigen ) * fac
1765  w(k,i_pidep ) = w(k,i_pidep ) * fac
1766  w(k,i_psdep ) = w(k,i_psdep ) * fac
1767  w(k,i_pgdep ) = w(k,i_pgdep ) * fac
1768  end do
1769 
1770  ! [QS]
1771  do k = ks, ke
1772  net = &
1773  + w(k,i_psdep ) & ! [prod] v->s
1774  + w(k,i_psacw ) & ! [prod] c->s
1775  + w(k,i_psfw ) & ! [prod] c->s
1776  + w(k,i_piacr_s) & ! [prod] r->s
1777  + w(k,i_psacr_s) & ! [prod] r->s
1778  + w(k,i_psaut ) & ! [prod] i->s
1779  + w(k,i_praci_s) & ! [prod] i->s
1780  + w(k,i_psaci ) & ! [prod] i->s
1781  + w(k,i_psfi ) & ! [prod] i->s
1782  - w(k,i_pssub ) & ! [loss] s->v
1783  - w(k,i_psmlt ) & ! [loss] s->r
1784  - w(k,i_pgaut ) & ! [loss] s->g
1785  - w(k,i_pracs ) & ! [loss] s->g
1786  - w(k,i_pgacs ) ! [loss] s->g
1787 
1788  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1789  fack(k) = ( fac_sw ) &
1790  + ( 1.0_rp - fac_sw ) * min( -w(k,i_dqs_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1791  end do
1792  do k = ks, ke
1793  w(k,i_psdep ) = w(k,i_psdep ) * fack(k)
1794  w(k,i_psacw ) = w(k,i_psacw ) * fack(k)
1795  w(k,i_psfw ) = w(k,i_psfw ) * fack(k)
1796  w(k,i_piacr_s) = w(k,i_piacr_s) * fack(k)
1797  w(k,i_psacr_s) = w(k,i_psacr_s) * fack(k)
1798  w(k,i_psaut ) = w(k,i_psaut ) * fack(k)
1799  w(k,i_praci_s) = w(k,i_praci_s) * fack(k)
1800  w(k,i_psaci ) = w(k,i_psaci ) * fack(k)
1801  w(k,i_psfi ) = w(k,i_psfi ) * fack(k)
1802  w(k,i_pssub ) = w(k,i_pssub ) * fack(k)
1803  w(k,i_psmlt ) = w(k,i_psmlt ) * fack(k)
1804  w(k,i_pgaut ) = w(k,i_pgaut ) * fack(k)
1805  w(k,i_pracs ) = w(k,i_pracs ) * fack(k)
1806  w(k,i_pgacs ) = w(k,i_pgacs ) * fack(k)
1807  end do
1808 
1809  ! [QG]
1810  do k = ks, ke
1811  net = &
1812  + w(k,i_pgdep ) & ! [prod] v->g
1813  + w(k,i_pgacw ) & ! [prod] c->g
1814  + w(k,i_piacr_g) & ! [prod] r->g
1815  + w(k,i_psacr_g) & ! [prod] r->g
1816  + w(k,i_pgacr ) & ! [prod] r->g
1817  + w(k,i_pgfrz ) & ! [prod] r->g
1818  + w(k,i_praci_g) & ! [prod] i->g
1819  + w(k,i_pgaci ) & ! [prod] i->g
1820  + w(k,i_pgaut ) & ! [prod] s->g
1821  + w(k,i_pracs ) & ! [prod] s->g
1822  + w(k,i_pgacs ) & ! [prod] s->g
1823  - w(k,i_pgsub ) & ! [loss] g->v
1824  - w(k,i_pgmlt ) ! [loss] g->r
1825 
1826  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1827  fack(k) = ( fac_sw ) &
1828  + ( 1.0_rp - fac_sw ) * min( -w(k,i_dqg_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1829  end do
1830  do k = ks, ke
1831  w(k,i_pgdep ) = w(k,i_pgdep ) * fack(k)
1832  w(k,i_pgacw ) = w(k,i_pgacw ) * fack(k)
1833  w(k,i_piacr_g) = w(k,i_piacr_g) * fack(k)
1834  w(k,i_psacr_g) = w(k,i_psacr_g) * fack(k)
1835  w(k,i_pgacr ) = w(k,i_pgacr ) * fack(k)
1836  w(k,i_pgfrz ) = w(k,i_pgfrz ) * fack(k)
1837  w(k,i_praci_g) = w(k,i_praci_g) * fack(k)
1838  w(k,i_pgaci ) = w(k,i_pgaci ) * fack(k)
1839  w(k,i_pgaut ) = w(k,i_pgaut ) * fack(k)
1840  w(k,i_pracs ) = w(k,i_pracs ) * fack(k)
1841  w(k,i_pgacs ) = w(k,i_pgacs ) * fack(k)
1842  w(k,i_pgsub ) = w(k,i_pgsub ) * fack(k)
1843  w(k,i_pgmlt ) = w(k,i_pgmlt ) * fack(k)
1844  end do
1845 
1846  do k = ks, ke
1847  qc_t(k) = + w(k,i_pimlt ) & ! [prod] i->c
1848  - w(k,i_praut ) & ! [loss] c->r
1849  - w(k,i_pracw ) & ! [loss] c->r
1850  - w(k,i_pihom ) & ! [loss] c->i
1851  - w(k,i_pihtr ) & ! [loss] c->i
1852  - w(k,i_psacw ) & ! [loss] c->s
1853  - w(k,i_psfw ) & ! [loss] c->s
1854  - w(k,i_pgacw ) ! [loss] c->g
1855  qc_t(k) = max( qc_t(k), -w(k,i_dqc_dt) )
1856  end do
1857 
1858  do k = ks, ke
1859  qr_t(k) = + w(k,i_praut ) & ! [prod] c->r
1860  + w(k,i_pracw ) & ! [prod] c->r
1861  + w(k,i_psmlt ) & ! [prod] s->r
1862  + w(k,i_pgmlt ) & ! [prod] g->r
1863  - w(k,i_prevp ) & ! [loss] r->v
1864  - w(k,i_piacr_s) & ! [loss] r->s
1865  - w(k,i_psacr_s) & ! [loss] r->s
1866  - w(k,i_piacr_g) & ! [loss] r->g
1867  - w(k,i_psacr_g) & ! [loss] r->g
1868  - w(k,i_pgacr ) & ! [loss] r->g
1869  - w(k,i_pgfrz ) ! [loss] r->g
1870 
1871  qr_t(k) = max( qr_t(k), -w(k,i_dqr_dt) )
1872  end do
1873 
1874  do k = ks, ke
1875  qi_t(k) = + w(k,i_pigen ) & ! [prod] v->i
1876  + w(k,i_pidep ) & ! [prod] v->i
1877  + w(k,i_pihom ) & ! [prod] c->i
1878  + w(k,i_pihtr ) & ! [prod] c->i
1879  - w(k,i_pisub ) & ! [loss] i->v
1880  - w(k,i_pimlt ) & ! [loss] i->c
1881  - w(k,i_psaut ) & ! [loss] i->s
1882  - w(k,i_praci_s) & ! [loss] i->s
1883  - w(k,i_psaci ) & ! [loss] i->s
1884  - w(k,i_psfi ) & ! [loss] i->s
1885  - w(k,i_praci_g) & ! [loss] i->g
1886  - w(k,i_pgaci ) ! [loss] i->g
1887  qi_t(k) = max( qi_t(k), -w(k,i_dqi_dt) )
1888  end do
1889 
1890  do k = ks, ke
1891  qs_t(k) = + w(k,i_psdep ) & ! [prod] v->s
1892  + w(k,i_psacw ) & ! [prod] c->s
1893  + w(k,i_psfw ) & ! [prod] c->s
1894  + w(k,i_piacr_s) & ! [prod] r->s
1895  + w(k,i_psacr_s) & ! [prod] r->s
1896  + w(k,i_psaut ) & ! [prod] i->s
1897  + w(k,i_praci_s) & ! [prod] i->s
1898  + w(k,i_psaci ) & ! [prod] i->s
1899  + w(k,i_psfi ) & ! [prod] i->s
1900  - w(k,i_pssub ) & ! [loss] s->v
1901  - w(k,i_psmlt ) & ! [loss] s->r
1902  - w(k,i_pgaut ) & ! [loss] s->g
1903  - w(k,i_pracs ) & ! [loss] s->g
1904  - w(k,i_pgacs ) ! [loss] s->g
1905  qs_t(k) = max( qs_t(k), -w(k,i_dqs_dt) )
1906  end do
1907 
1908  do k = ks, ke
1909  qg_t(k) = + w(k,i_pgdep ) & ! [prod] v->g
1910  + w(k,i_pgacw ) & ! [prod] c->g
1911  + w(k,i_piacr_g) & ! [prod] r->g
1912  + w(k,i_psacr_g) & ! [prod] r->g
1913  + w(k,i_pgacr ) & ! [prod] r->g
1914  + w(k,i_pgfrz ) & ! [prod] r->g
1915  + w(k,i_praci_g) & ! [prod] i->g
1916  + w(k,i_pgaci ) & ! [prod] i->g
1917  + w(k,i_pgaut ) & ! [prod] s->g
1918  + w(k,i_pracs ) & ! [prod] s->g
1919  + w(k,i_pgacs ) & ! [prod] s->g
1920  - w(k,i_pgsub ) & ! [loss] g->v
1921  - w(k,i_pgmlt ) ! [loss] g->r
1922  qg_t(k) = max( qg_t(k), -w(k,i_dqg_dt) )
1923  end do
1924 
1925  if ( flg_lt_l ) then
1926  do k = ks, ke
1927  ! store to work
1928  qcrg_c(k) = qtrc_crg0(k,i,j,i_qc-1)
1929  qcrg_r(k) = qtrc_crg0(k,i,j,i_qr-1)
1930  qcrg_i(k) = qtrc_crg0(k,i,j,i_qi-1)
1931  qcrg_s(k) = qtrc_crg0(k,i,j,i_qs-1)
1932  qcrg_g(k) = qtrc_crg0(k,i,j,i_qg-1)
1933  re_qs(k) = re(k,i,j,i_hs) * 1.0e-2_rp ! [cm] -> [m]
1934  beta1_crg(k) = beta_crg(k,i,j)
1935  dcrg(k) = - dqcrg(k,i,j)
1936  end do
1937 
1938  do k = ks, ke
1939  ! [Qgaci] charge separation of cloud ice
1940  alpha = 5.0_rp * ( 2.0_rp * re_qi / d0_crg )**2 * ( -vtg(k) ) / v0_crg
1941  alpha = min( alpha, 10.0_rp )
1942  w_qcrg(k,i_qgaci) = dens(k,i,j) * qi(k) * 1.5_rp &
1943  * cg * gam_3dg * rlmdg_3dg(k) * rho_fact(k) * rdens_i &
1944  * sign( min( abs(dcrg(k)*alpha),20.0_rp ),dcrg(k)*alpha ) * beta1_crg(k) &
1945  * n0g(k) / ( 2.0_rp * re_qi )**3 &
1946  * ( ( 1.0_rp - flg_ecoali ) * ( 1.0_rp-egi ) &
1947  + ( flg_ecoali ) * egi * ( 1.0_rp-ecoal_gi ) / ecoal_gi )
1948  qsplt_in(k,i,j,1) = qsplt_in(k,i,j,1) - w_qcrg(k,i_qgaci) !*dt
1949  qsplt_in(k,i,j,2) = qsplt_in(k,i,j,2) + w_qcrg(k,i_qgaci) !*dt
1950  end do
1951 
1952 !OCL LOOP_FISSION_TARGET(LS)
1953  do k = ks, ke
1954  ! [Qgacs] charge separation of snow
1955 ! w_qcrg(I_Qgacs) = 0.0_RP
1956  alpha = 5.0_rp * ( 2.0_rp * re_qs(k) / d0_crg )**2 * ( -vtg(k) ) / v0_crg
1957  alpha = min( alpha, 10.0_rp )
1958  w_qcrg(k,i_qgacs) = 0.25_rp * pi / dens(k,i,j) &
1959  * n0g(k) * n0s(k) * abs(vtg(k)-vts(k)) * beta1_crg(k) &
1960  * sign( min( abs(dcrg(k)*alpha),50.0_rp ),dcrg(k)*alpha ) &
1961  * ( gam_3 * gam * w(k,i_rlmds)**3 * rlmdg(k) &
1962  + gam_2 * gam_2 * w(k,i_rlmds)**2 * rlmdg_2(k) &
1963  + gam * gam_3 * w(k,i_rlmds) * rlmdg_3(k) ) &
1964  * ( ( 1.0_rp - flg_ecoals ) * ( 1.0_rp-egs_mod(k) ) &
1965  + ( flg_ecoals ) * egs_mod(k) * ( 1.0_rp-ecoal_gs ) / ecoal_gs )
1966  qsplt_in(k,i,j,1) = qsplt_in(k,i,j,1) - w_qcrg(k,i_qgacs) !*dt
1967  qsplt_in(k,i,j,3) = qsplt_in(k,i,j,3) + w_qcrg(k,i_qgacs) !*dt
1968  end do
1969 
1970  do k = ks, ke
1971  facq_qc = 0.5_rp + sign( 0.5_rp, qc(k) - eps )
1972  facq_qi = 0.5_rp + sign( 0.5_rp, qi(k) - eps )
1973  w_q(k,i_pimlt) = qcrg_i(k) * w(k,i_pimlt) / ( qi(k) + eps*eps ) * facq_qi
1974  w_q(k,i_praut) = qcrg_c(k) * w(k,i_praut) / ( qc(k) + eps*eps ) * facq_qc
1975  w_q(k,i_pracw) = qcrg_c(k) * w(k,i_pracw) / ( qc(k) + eps*eps ) * facq_qc
1976  w_q(k,i_pihom) = qcrg_c(k) * w(k,i_pihom) / ( qc(k) + eps*eps ) * facq_qc
1977  w_q(k,i_pihtr) = qcrg_c(k) * w(k,i_pihtr) / ( qc(k) + eps*eps ) * facq_qc
1978  w_q(k,i_psacw) = qcrg_c(k) * w(k,i_psacw) / ( qc(k) + eps*eps ) * facq_qc
1979  w_q(k,i_psfw ) = qcrg_c(k) * w(k,i_psfw ) / ( qc(k) + eps*eps ) * facq_qc
1980  w_q(k,i_pgacw) = qcrg_c(k) * w(k,i_pgacw) / ( qc(k) + eps*eps ) * facq_qc
1981  end do
1982 
1983  do k = ks, ke
1984  facq_qc = 0.5_rp + sign( 0.5_rp, qc(k) - eps )
1985  facq_qi = 0.5_rp + sign( 0.5_rp, qi(k) - eps )
1986  w_q(k,i_pigen ) = 0.0_rp
1987  w_q(k,i_pidep ) = 0.0_rp
1988  w_q(k,i_pihom ) = qcrg_c(k) * w(k,i_pihom ) / ( qc(k) + eps*eps ) * facq_qc
1989  w_q(k,i_pihtr ) = qcrg_c(k) * w(k,i_pihtr ) / ( qc(k) + eps*eps ) * facq_qc
1990  w_q(k,i_pisub ) = qcrg_i(k) * w(k,i_pisub ) / ( qi(k) + eps*eps ) * facq_qi
1991  w_q(k,i_pimlt ) = qcrg_i(k) * w(k,i_pimlt ) / ( qi(k) + eps*eps ) * facq_qi
1992  w_q(k,i_psaut ) = qcrg_i(k) * w(k,i_psaut ) / ( qi(k) + eps*eps ) * facq_qi
1993  w_q(k,i_praci_s) = qcrg_i(k) * w(k,i_praci_s) / ( qi(k) + eps*eps ) * facq_qi
1994  w_q(k,i_psaci ) = qcrg_i(k) * w(k,i_psaci ) / ( qi(k) + eps*eps ) * facq_qi
1995  w_q(k,i_psfi ) = qcrg_i(k) * w(k,i_psfi ) / ( qi(k) + eps*eps ) * facq_qi
1996  w_q(k,i_praci_g) = qcrg_i(k) * w(k,i_praci_g) / ( qi(k) + eps*eps ) * facq_qi
1997  w_q(k,i_pgaci ) = qcrg_i(k) * w(k,i_pgaci ) / ( qi(k) + eps*eps ) * facq_qi
1998  end do
1999 
2000  do k = ks, ke
2001  facq_qc = 0.5_rp + sign( 0.5_rp, qc(k) - eps )
2002  facq_qr = 0.5_rp + sign( 0.5_rp, qr(k) - eps )
2003  facq_qs = 0.5_rp + sign( 0.5_rp, qs(k) - eps )
2004  facq_qg = 0.5_rp + sign( 0.5_rp, qg(k) - eps )
2005  w_q(k,i_praut ) = qcrg_c(k) * w(k,i_praut ) / ( qc(k) + eps*eps ) * facq_qc
2006  w_q(k,i_pracw ) = qcrg_c(k) * w(k,i_pracw ) / ( qc(k) + eps*eps ) * facq_qc
2007  w_q(k,i_psmlt ) = qcrg_s(k) * w(k,i_psmlt ) / ( qs(k) + eps*eps ) * facq_qs
2008  w_q(k,i_pgmlt ) = qcrg_g(k) * w(k,i_pgmlt ) / ( qg(k) + eps*eps ) * facq_qg
2009  w_q(k,i_prevp ) = qcrg_r(k) * w(k,i_prevp ) / ( qr(k) + eps*eps ) * facq_qr
2010  w_q(k,i_piacr_s) = qcrg_r(k) * w(k,i_piacr_s) / ( qr(k) + eps*eps ) * facq_qr
2011  w_q(k,i_psacr_s) = qcrg_r(k) * w(k,i_psacr_s) / ( qr(k) + eps*eps ) * facq_qr
2012  w_q(k,i_piacr_g) = qcrg_r(k) * w(k,i_piacr_g) / ( qr(k) + eps*eps ) * facq_qr
2013  w_q(k,i_psacr_g) = qcrg_r(k) * w(k,i_psacr_g) / ( qr(k) + eps*eps ) * facq_qr
2014  w_q(k,i_pgacr ) = qcrg_r(k) * w(k,i_pgacr ) / ( qr(k) + eps*eps ) * facq_qr
2015  w_q(k,i_pgfrz ) = qcrg_r(k) * w(k,i_pgfrz ) / ( qr(k) + eps*eps ) * facq_qr
2016  end do
2017 
2018  do k = ks, ke
2019  facq_qc = 0.5_rp + sign( 0.5_rp, qc(k) - eps )
2020  facq_qr = 0.5_rp + sign( 0.5_rp, qr(k) - eps )
2021  facq_qi = 0.5_rp + sign( 0.5_rp, qi(k) - eps )
2022  facq_qs = 0.5_rp + sign( 0.5_rp, qs(k) - eps )
2023  facq_qg = 0.5_rp + sign( 0.5_rp, qg(k) - eps )
2024  w_q(k,i_psdep ) = 0.0_rp
2025  w_q(k,i_psacw ) = qcrg_c(k) * w(k,i_psacw ) / ( qc(k) + eps*eps ) * facq_qc
2026  w_q(k,i_psfw ) = qcrg_c(k) * w(k,i_psfw ) / ( qc(k) + eps*eps ) * facq_qc
2027  w_q(k,i_piacr_s) = qcrg_r(k) * w(k,i_piacr_s) / ( qr(k) + eps*eps ) * facq_qr
2028  w_q(k,i_psacr_s) = qcrg_r(k) * w(k,i_psacr_s) / ( qr(k) + eps*eps ) * facq_qr
2029  w_q(k,i_psaut ) = qcrg_i(k) * w(k,i_psaut ) / ( qi(k) + eps*eps ) * facq_qi
2030  w_q(k,i_praci_s) = qcrg_i(k) * w(k,i_praci_s) / ( qi(k) + eps*eps ) * facq_qi
2031  w_q(k,i_psaci ) = qcrg_i(k) * w(k,i_psaci ) / ( qi(k) + eps*eps ) * facq_qi
2032  w_q(k,i_psfi ) = qcrg_i(k) * w(k,i_psfi ) / ( qi(k) + eps*eps ) * facq_qi
2033  w_q(k,i_pssub ) = qcrg_s(k) * w(k,i_pssub ) / ( qs(k) + eps*eps ) * facq_qs
2034  w_q(k,i_psmlt ) = qcrg_s(k) * w(k,i_psmlt ) / ( qs(k) + eps*eps ) * facq_qs
2035  w_q(k,i_pgaut ) = qcrg_s(k) * w(k,i_pgaut ) / ( qs(k) + eps*eps ) * facq_qs
2036  w_q(k,i_pracs ) = qcrg_s(k) * w(k,i_pracs ) / ( qs(k) + eps*eps ) * facq_qs
2037  w_q(k,i_pgacs ) = qcrg_s(k) * w(k,i_pgacs ) / ( qs(k) + eps*eps ) * facq_qs
2038  end do
2039 
2040  do k = ks, ke
2041  facq_qc = 0.5_rp + sign( 0.5_rp, qc(k) - eps )
2042  facq_qr = 0.5_rp + sign( 0.5_rp, qr(k) - eps )
2043  facq_qi = 0.5_rp + sign( 0.5_rp, qi(k) - eps )
2044  facq_qs = 0.5_rp + sign( 0.5_rp, qs(k) - eps )
2045  facq_qg = 0.5_rp + sign( 0.5_rp, qg(k) - eps )
2046  w_q(k,i_pgdep ) = 0.0_rp
2047  w_q(k,i_pgacw ) = qcrg_c(k) * w(k,i_pgacw ) / ( qc(k) + eps*eps ) * facq_qc
2048  w_q(k,i_piacr_g) = qcrg_r(k) * w(k,i_piacr_g) / ( qr(k) + eps*eps ) * facq_qr
2049  w_q(k,i_psacr_g) = qcrg_r(k) * w(k,i_psacr_g) / ( qr(k) + eps*eps ) * facq_qr
2050  w_q(k,i_pgacr ) = qcrg_r(k) * w(k,i_pgacr ) / ( qr(k) + eps*eps ) * facq_qr
2051  w_q(k,i_pgfrz ) = qcrg_r(k) * w(k,i_pgfrz ) / ( qr(k) + eps*eps ) * facq_qr
2052  w_q(k,i_praci_g) = qcrg_i(k) * w(k,i_praci_g) / ( qi(k) + eps*eps ) * facq_qi
2053  w_q(k,i_pgaci ) = qcrg_i(k) * w(k,i_pgaci ) / ( qi(k) + eps*eps ) * facq_qi
2054  w_q(k,i_pgaut ) = qcrg_s(k) * w(k,i_pgaut ) / ( qs(k) + eps*eps ) * facq_qs
2055  w_q(k,i_pracs ) = qcrg_s(k) * w(k,i_pracs ) / ( qs(k) + eps*eps ) * facq_qs
2056  w_q(k,i_pgacs ) = qcrg_s(k) * w(k,i_pgacs ) / ( qs(k) + eps*eps ) * facq_qs
2057  w_q(k,i_pgsub ) = qcrg_g(k) * w(k,i_pgsub ) / ( qg(k) + eps*eps ) * facq_qg
2058  w_q(k,i_pgmlt ) = qcrg_g(k) * w(k,i_pgmlt ) / ( qg(k) + eps*eps ) * facq_qg
2059  end do
2060 
2061  do k = ks, ke
2062  qc_crg_t(k) = + w_q(k,i_pimlt ) & ! [prod] i->c
2063  - w_q(k,i_praut ) & ! [loss] c->r
2064  - w_q(k,i_pracw ) & ! [loss] c->r
2065  - w_q(k,i_pihom ) & ! [loss] c->i
2066  - w_q(k,i_pihtr ) & ! [loss] c->i
2067  - w_q(k,i_psacw ) & ! [loss] c->s
2068  - w_q(k,i_psfw ) & ! [loss] c->s
2069  - w_q(k,i_pgacw ) ! [loss] c->g
2070  end do
2071 
2072  do k = ks, ke
2073  qr_crg_t(k) = + w_q(k,i_praut ) & ! [prod] c->r
2074  + w_q(k,i_pracw ) & ! [prod] c->r
2075  + w_q(k,i_psmlt ) & ! [prod] s->r
2076  + w_q(k,i_pgmlt ) & ! [prod] g->r
2077  - w_q(k,i_prevp ) & ! [loss] r->v
2078  - w_q(k,i_piacr_s) & ! [loss] r->s
2079  - w_q(k,i_psacr_s) & ! [loss] r->s
2080  - w_q(k,i_piacr_g) & ! [loss] r->g
2081  - w_q(k,i_psacr_g) & ! [loss] r->g
2082  - w_q(k,i_pgacr ) & ! [loss] r->g
2083  - w_q(k,i_pgfrz ) ! [loss] r->g
2084  end do
2085 
2086  do k = ks, ke
2087  qi_crg_t(k) = + w_q(k,i_pigen ) & ! [prod] v->i
2088  + w_q(k,i_pidep ) & ! [prod] v->i
2089  + w_q(k,i_pihom ) & ! [prod] c->i
2090  + w_q(k,i_pihtr ) & ! [prod] c->i
2091  - w_q(k,i_pisub ) & ! [loss] i->v
2092  - w_q(k,i_pimlt ) & ! [loss] i->c
2093  - w_q(k,i_psaut ) & ! [loss] i->s
2094  - w_q(k,i_praci_s) & ! [loss] i->s
2095  - w_q(k,i_psaci ) & ! [loss] i->s
2096  - w_q(k,i_psfi ) & ! [loss] i->s
2097  - w_q(k,i_praci_g) & ! [loss] i->g
2098  - w_q(k,i_pgaci ) & ! [loss] i->g
2099  + w_qcrg(k,i_qgaci ) ! [prod] Charge Split by g-i coll.
2100  end do
2101 
2102  do k = ks, ke
2103  qs_crg_t(k) = + w_q(k,i_psdep ) & ! [prod] v->s
2104  + w_q(k,i_psacw ) & ! [prod] c->s
2105  + w_q(k,i_psfw ) & ! [prod] c->s
2106  + w_q(k,i_piacr_s) & ! [prod] r->s
2107  + w_q(k,i_psacr_s) & ! [prod] r->s
2108  + w_q(k,i_psaut ) & ! [prod] i->s
2109  + w_q(k,i_praci_s) & ! [prod] i->s
2110  + w_q(k,i_psaci ) & ! [prod] i->s
2111  + w_q(k,i_psfi ) & ! [prod] i->s
2112  - w_q(k,i_pssub ) & ! [loss] s->v
2113  - w_q(k,i_psmlt ) & ! [loss] s->r
2114  - w_q(k,i_pgaut ) & ! [loss] s->g
2115  - w_q(k,i_pracs ) & ! [loss] s->g
2116  - w_q(k,i_pgacs ) & ! [loss] s->g
2117  + w_qcrg(k,i_qgacs ) ! [prod] Charge Split by g-s coll.
2118  end do
2119 
2120  do k = ks, ke
2121  qg_crg_t(k) = + w_q(k,i_pgdep ) & ! [prod] v->g
2122  + w_q(k,i_pgacw ) & ! [prod] c->g
2123  + w_q(k,i_piacr_g) & ! [prod] r->g
2124  + w_q(k,i_psacr_g) & ! [prod] r->g
2125  + w_q(k,i_pgacr ) & ! [prod] r->g
2126  + w_q(k,i_pgfrz ) & ! [prod] r->g
2127  + w_q(k,i_praci_g) & ! [prod] i->g
2128  + w_q(k,i_pgaci ) & ! [prod] i->g
2129  + w_q(k,i_pgaut ) & ! [prod] s->g
2130  + w_q(k,i_pracs ) & ! [prod] s->g
2131  + w_q(k,i_pgacs ) & ! [prod] s->g
2132  - w_q(k,i_pgsub ) & ! [loss] g->v
2133  - w_q(k,i_pgmlt ) & ! [loss] g->r
2134  - w_qcrg(k,i_qgaci) & ! [prod] Charge Split by g-i coll.
2135  - w_qcrg(k,i_qgacs) ! [prod] Charge Split by g-s coll.
2136  end do
2137  end if
2138 
2139  do k = ks, ke
2140  qv_t(k) = - ( qc_t(k) + qr_t(k) + qi_t(k) + qs_t(k) + qg_t(k) )
2141  end do
2142 
2143  do k = ks, ke
2144  qtrc(k,i,j,i_qv) = qtrc(k,i,j,i_qv) + qv_t(k) * dt
2145  end do
2146  do k = ks, ke
2147  qtrc(k,i,j,i_qc) = qtrc(k,i,j,i_qc) + qc_t(k) * dt
2148  end do
2149  do k = ks, ke
2150  qtrc(k,i,j,i_qr) = qtrc(k,i,j,i_qr) + qr_t(k) * dt
2151  end do
2152  do k = ks, ke
2153  qtrc(k,i,j,i_qi) = qtrc(k,i,j,i_qi) + qi_t(k) * dt
2154  end do
2155  do k = ks, ke
2156  qtrc(k,i,j,i_qs) = qtrc(k,i,j,i_qs) + qs_t(k) * dt
2157  end do
2158  do k = ks, ke
2159  qtrc(k,i,j,i_qg) = qtrc(k,i,j,i_qg) + qg_t(k) * dt
2160  end do
2161  do k = ks, ke
2162  cv_t = cv_vapor * qv_t(k) &
2163  + cv_water * ( qc_t(k) + qr_t(k) ) &
2164  + cv_ice * ( qi_t(k) + qs_t(k) + qg_t(k) )
2165  cvtot(k) = cvtot0(k,i,j) + cv_t * dt
2166  end do
2167 
2168  do k = ks, ke
2169  e_t = - lhv * qv_t(k) + lhf * ( qi_t(k) + qs_t(k) + qg_t(k) ) ! internal energy change
2170  rhoe_t(k,i,j) = dens(k,i,j) * e_t
2171  temp(k,i,j) = ( temp(k,i,j) * cvtot0(k,i,j) + e_t * dt ) / cvtot(k)
2172  end do
2173 
2174  do k = ks, ke
2175  cvtot0(k,i,j) = cvtot(k)
2176  end do
2177 
2178  do k = ks, ke
2179  cp_t = cp_vapor * qv_t(k) &
2180  + cp_water * ( qc_t(k) + qr_t(k) ) &
2181  + cp_ice * ( qi_t(k) + qs_t(k) + qg_t(k) )
2182  cptot0(k,i,j) = cptot0(k,i,j) + cp_t * dt
2183  end do
2184 
2185  if (flg_lt_l) then
2186  do k = ks, ke
2187  qtrc_crg0(k,i,j,i_qc-1) = qtrc_crg0(k,i,j,i_qc-1) + qc_crg_t(k) * dt
2188  qtrc_crg0(k,i,j,i_qr-1) = qtrc_crg0(k,i,j,i_qr-1) + qr_crg_t(k) * dt
2189  qtrc_crg0(k,i,j,i_qi-1) = qtrc_crg0(k,i,j,i_qi-1) + qi_crg_t(k) * dt
2190  qtrc_crg0(k,i,j,i_qs-1) = qtrc_crg0(k,i,j,i_qs-1) + qs_crg_t(k) * dt
2191  qtrc_crg0(k,i,j,i_qg-1) = qtrc_crg0(k,i,j,i_qg-1) + qg_crg_t(k) * dt
2192  end do
2193 
2194 !OCL LOOP_FISSION_TARGET(LS)
2195  !$acc loop private(rlambda)
2196  do k = ks, ke
2197  rlambda(i_qr) = sqrt(sqrt( dens(k,i,j) * max( qtrc(k,i,j,i_qr),0.0_rp ) / ( ar * n0r(k) * gam_1br ) ))
2198  rlambda(i_qs) = sqrt(sqrt( dens(k,i,j) * max( qtrc(k,i,j,i_qs),0.0_rp ) / ( as * n0s(k) * gam_1bs ) ))
2199  rlambda(i_qg) = sqrt(sqrt( dens(k,i,j) * max( qtrc(k,i,j,i_qg),0.0_rp ) / ( ag * n0g(k) * gam_1bg ) ))
2200  sarea(k,i,j,i_qc-1) = 6.0_rp / ( re_qc*2.0_rp * dwatr ) &
2201  * max( qtrc(k,i,j,i_qc),0.0_rp ) * dens(k,i,j)
2202  sarea(k,i,j,i_qi-1) = 6.0_rp / ( re_qi*2.0_rp * dice ) &
2203  * max( qtrc(k,i,j,i_qi),0.0_rp ) * dens(k,i,j)
2204  sarea(k,i,j,i_qr-1) = pi * n0r(k) * gam_3 * rlambda(i_qr)**3
2205  sarea(k,i,j,i_qs-1) = pi * n0s(k) * gam_3 * rlambda(i_qs)**3
2206  sarea(k,i,j,i_qg-1) = pi * n0g(k) * gam_3 * rlambda(i_qg)**3
2207  end do
2208 
2209  do iq = 1, 3
2210  do k = ks, ke
2211  qsplt_in(k,i,j,iq) = qsplt_in(k,i,j,iq) * dens(k,i,j)
2212  enddo
2213  enddo
2214  end if
2215 
2216  if ( hist_flag ) then
2217  !$acc loop seq
2218  do ip = 1, w_nmax
2219  if ( hist_sw(ip) ) then
2220  do k = ks, ke
2221  w3d(k,i,j,ip) = w(k,ip)
2222  enddo
2223  end if
2224  enddo
2225  end if
2226 
2227  enddo
2228  enddo
2229  !$acc end kernels
2230 
2231  do ip = 1, w_nmax
2232  if ( hist_sw(ip) ) call file_history_put( hist_id(ip), w3d(:,:,:,ip) )
2233  enddo
2234 
2235  if( flg_lt_l ) then
2236  !$acc exit data copyout(QSPLT_in, Sarea, QTRC_crg0) &
2237  !$acc delete(Re, dqcrg, beta_crg)
2238  else
2239  !$acc exit data delete(QSPLT_in, Sarea, QTRC_crg0) &
2240  !$acc delete(Re, dqcrg, beta_crg)
2241  end if
2242 
2243  !$acc end data
2244 
2245  call prof_rapend ('MP_tomita08', 3)
2246 
2247  return
2248  end subroutine mp_tomita08
2249 
2250  !-----------------------------------------------------------------------------
2252 !OCL SERIAL
2253  subroutine atmos_phy_mp_tomita08_terminal_velocity( &
2254  KA, KS, KE, &
2255  DENS, TEMP, RHOQ, &
2256  vterm )
2257  !$acc routine vector
2258  use scale_const, only: &
2259  tem00 => const_tem00
2260  implicit none
2261  integer, intent(in), value :: ka, ks, ke
2262 
2263  real(rp), intent(in) :: dens(ka)
2264  real(rp), intent(in) :: temp(ka)
2265  real(rp), intent(in) :: rhoq(ka,qa_mp-1)
2266 
2267  real(rp), intent(out) :: vterm(ka,qa_mp-1)
2268 
2269 #ifdef _OPENACC
2270  real(rp) :: temc_tv
2271  real(rp) :: qr_tv
2272  real(rp) :: qi_tv
2273  real(rp) :: qs_tv
2274  real(rp) :: qg_tv
2275  real(rp) :: rho_fact_tv
2276 #define temc_tv(k) temc_tv
2277 #define qr_tv(k) qr_tv
2278 #define qi_tv(k) qi_tv
2279 #define qs_tv(k) qs_tv
2280 #define qg_tv(k) qg_tv
2281 #define rho_fact_tv(k) rho_fact_tv
2282 #else
2283  real(rp) :: temc_tv(ka)
2284  real(rp) :: qr_tv(ka)
2285  real(rp) :: qi_tv(ka)
2286  real(rp) :: qs_tv(ka)
2287  real(rp) :: qg_tv(ka)
2288  real(rp) :: rho_fact_tv(ka) ! density factor
2289 #endif
2290 
2291 #ifdef _OPENACC
2292  real(rp) :: n0r_tv
2293  real(rp) :: n0s_tv
2294  real(rp) :: n0g_tv
2295 #define N0r_tv(k) N0r_tv
2296 #define N0s_tv(k) N0s_tv
2297 #define N0g_tv(k) N0g_tv
2298 #else
2299  real(rp) :: n0r_tv(ka)
2300  real(rp) :: n0s_tv(ka)
2301  real(rp) :: n0g_tv(ka)
2302 #endif
2303  real(rp) :: rlmdr, rlmds, rlmdg
2304  real(rp) :: rlmdr_dr, rlmds_ds, rlmdg_dg
2305 
2306  !---< Roh and Satoh (2014) >---
2307  real(rp) :: tems, xs2
2308  real(rp) :: moms_0bs
2309 #ifdef _OPENACC
2310  real(rp) :: rmoms_vt_tv
2311 #define RMOMs_Vt_tv(k) RMOMs_Vt_tv
2312 #else
2313  real(rp) :: rmoms_vt_tv(ka)
2314 #endif
2315  real(rp) :: coef_at(4), coef_bt(4)
2316  real(rp) :: loga_, b_, nm
2317 
2318  real(rp) :: zerosw
2319  integer :: k
2320  !---------------------------------------------------------------------------
2321 
2322  !$acc loop private(coef_at, coef_bt)
2323  do k = ks, ke
2324  ! store to work
2325  temc_tv(k) = temp(k) - tem00
2326  qr_tv(k) = rhoq(k,i_hyd_qr) / dens(k)
2327  qi_tv(k) = rhoq(k,i_hyd_qi) / dens(k)
2328  qs_tv(k) = rhoq(k,i_hyd_qs) / dens(k)
2329  qg_tv(k) = rhoq(k,i_hyd_qg) / dens(k)
2330 
2331  rho_fact_tv(k) = sqrt( dens00 / dens(k) )
2332 #ifndef _OPENACC
2333  end do
2334 #endif
2335 
2336  if ( enable_wdxz2014 ) then
2337  ! Wainwright et al. (2014)
2338  ! intercept parameter N0
2339 #ifndef _OPENACC
2340  do k = ks, ke
2341 #endif
2342  n0r_tv(k) = 1.16e+5_rp * exp( log( max( dens(k)*qr_tv(k)*1000.0_rp, 1.e-2_rp ) )*0.477_rp )
2343  n0s_tv(k) = 4.58e+9_rp * exp( log( max( dens(k)*qs_tv(k)*1000.0_rp, 1.e-2_rp ) )*0.788_rp )
2344  n0g_tv(k) = 9.74e+8_rp * exp( log( max( dens(k)*qg_tv(k)*1000.0_rp, 1.e-2_rp ) )*0.816_rp )
2345 #ifndef _OPENACC
2346  end do
2347 #endif
2348  else
2349  ! intercept parameter N0
2350 #ifndef _OPENACC
2351  do k = ks, ke
2352 #endif
2353  n0r_tv(k) = n0r_def ! Marshall and Palmer (1948)
2354  n0s_tv(k) = n0s_def ! Gunn and Marshall (1958)
2355  n0g_tv(k) = n0g_def
2356 #ifndef _OPENACC
2357  end do
2358 #endif
2359  end if
2360 
2361  !---< terminal velocity >
2362  if ( enable_hzdfhi2007 ) then
2363 #ifndef _OPENACC
2364  do k = ks, ke
2365 #endif
2366  zerosw = 0.5_rp - sign(0.5_rp, qi_tv(k) - 1.e-8_rp )
2367  vterm(k,i_hyd_qi) = min( 0.0_rp, &
2368  - ( coef_a0 + coef_a1 * temc_tv(k) ) &
2369  * exp( log( dens(k)*qi_tv(k)*1000.0_rp+zerosw ) * ( coef_b0 + coef_b1 * temc_tv(k) ) ) &
2370  * 1e-2_rp * ( 1.0_rp-zerosw ) )
2371 #ifndef _OPENACC
2372  end do
2373 #endif
2374  else
2375 #ifndef _OPENACC
2376  do k = ks, ke
2377 #endif
2378  zerosw = 0.5_rp - sign(0.5_rp, qi_tv(k) - 1.e-8_rp )
2379  vterm(k,i_hyd_qi) = -3.29_rp * exp( log( dens(k)*qi_tv(k)+zerosw )*0.16_rp ) * ( 1.0_rp-zerosw )
2380 #ifndef _OPENACC
2381  end do
2382 #endif
2383  end if
2384 
2385 
2386  ! slope parameter lambda (Rain)
2387 #ifndef _OPENACC
2388  do k = ks, ke
2389 #endif
2390  zerosw = 0.5_rp - sign(0.5_rp, qr_tv(k) - 1.e-12_rp )
2391  rlmdr = sqrt(sqrt( dens(k) * qr_tv(k) / ( ar * n0r_tv(k) * gam_1br ) + zerosw )) * ( 1.0_rp-zerosw )
2392  rlmdr_dr = sqrt( rlmdr ) ! **Dr
2393  vterm(k,i_hyd_qr) = -cr * rho_fact_tv(k) * gam_1brdr / gam_1br * rlmdr_dr
2394 #ifndef _OPENACC
2395  end do
2396 #endif
2397 
2398 
2399  if ( enable_rs2014 ) then
2400  !---< modification by Roh and Satoh (2014) >---
2401  ! bimodal size distribution of snow
2402 #ifndef _OPENACC
2403  do k = ks, ke
2404 #endif
2405  zerosw = 0.5_rp - sign(0.5_rp, dens(k) * qs_tv(k) - 1.e-12_rp )
2406  xs2 = dens(k) * qs_tv(k) / as
2407 
2408  tems = min( -0.1_rp, temc_tv(k) )
2409  coef_at(1) = coef_a01 + tems * ( coef_a02 + tems * ( coef_a05 + tems * coef_a09 ) )
2410  coef_at(2) = coef_a03 + tems * ( coef_a04 + tems * coef_a07 )
2411  coef_at(3) = coef_a06 + tems * coef_a08
2412  coef_at(4) = coef_a10
2413  coef_bt(1) = coef_b01 + tems * ( coef_b02 + tems * ( coef_b05 + tems * coef_b09 ) )
2414  coef_bt(2) = coef_b03 + tems * ( coef_b04 + tems * coef_b07 )
2415  coef_bt(3) = coef_b06 + tems * coef_b08
2416  coef_bt(4) = coef_b10
2417  ! 0 + Bs(=2) moment
2418  nm = 2.0_rp
2419  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
2420  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
2421  moms_0bs = exp( ln10 * loga_ + log(xs2+zerosw) * b_ ) * ( 1.0_rp-zerosw )
2422  ! Bs(=2) + Ds(=0.25) moment
2423  nm = 2.25_rp
2424  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
2425  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
2426  rmoms_vt_tv(k) = exp( ln10 * loga_ + log(xs2+zerosw) * b_ ) * ( 1.0_rp-zerosw ) / ( moms_0bs + zerosw )
2427 #ifndef _OPENACC
2428  end do
2429 #endif
2430  else
2431  ! slope parameter lambda (Snow)
2432 #ifndef _OPENACC
2433  do k = ks, ke
2434 #endif
2435  zerosw = 0.5_rp - sign(0.5_rp, qs_tv(k) - 1.e-12_rp )
2436  rlmds = sqrt(sqrt( dens(k) * qs_tv(k) / ( as * n0s_tv(k) * gam_1bs ) + zerosw )) * ( 1.0_rp-zerosw )
2437  rlmds_ds = sqrt( sqrt(rlmds) ) ! **Ds
2438  rmoms_vt_tv(k) = gam_1bsds / gam_1bs * rlmds_ds
2439 #ifndef _OPENACC
2440  end do
2441 #endif
2442  end if
2443 
2444 #ifndef _OPENACC
2445  do k = ks, ke
2446 #endif
2447  vterm(k,i_hyd_qs) = -cs * rho_fact_tv(k) * rmoms_vt_tv(k)
2448 #ifndef _OPENACC
2449  end do
2450 
2451  do k = ks, ke
2452 #endif
2453  ! slope parameter lambda (Graupel)
2454  zerosw = 0.5_rp - sign(0.5_rp, qg_tv(k) - 1.e-12_rp )
2455  rlmdg = sqrt(sqrt( dens(k) * qg_tv(k) / ( ag * n0g_tv(k) * gam_1bg ) + zerosw )) * ( 1.0_rp-zerosw )
2456  rlmdg_dg = sqrt( rlmdg ) ! **Dg
2457  vterm(k,i_hyd_qg) = -cg * rho_fact_tv(k) * gam_1bgdg / gam_1bg * rlmdg_dg
2458 #ifndef _OPENACC
2459  enddo
2460 
2461 
2462 !OCL XFILL
2463  do k = ks, ke
2464 #endif
2465  vterm(k,i_hyd_qc) = 0.0_rp
2466 #ifndef _OPENACC
2467  end do
2468 #endif
2469 
2470  if ( nofall_qr ) then
2471 !OCL XFILL
2472 #ifndef _OPENACC
2473  do k = ks, ke
2474 #endif
2475  vterm(k,i_hyd_qr) = 0.0_rp
2476 #ifndef _OPENACC
2477  enddo
2478 #endif
2479  endif
2480 
2481  if ( nofall_qi ) then
2482 !OCL XFILL
2483 #ifndef _OPENACC
2484  do k = ks, ke
2485 #endif
2486  vterm(k,i_hyd_qi) = 0.0_rp
2487 #ifndef _OPENACC
2488  enddo
2489 #endif
2490  endif
2491 
2492  if ( nofall_qs ) then
2493 !OCL XFILL
2494 #ifndef _OPENACC
2495  do k = ks, ke
2496 #endif
2497  vterm(k,i_hyd_qs) = 0.0_rp
2498 #ifndef _OPENACC
2499  enddo
2500 #endif
2501  endif
2502 
2503  if ( nofall_qg ) then
2504 !OCL XFILL
2505 #ifndef _OPENACC
2506  do k = ks, ke
2507 #endif
2508  vterm(k,i_hyd_qg) = 0.0_rp
2509 #ifndef _OPENACC
2510  enddo
2511 #endif
2512  endif
2513 
2514 #ifdef _OPENACC
2515  end do
2516 #endif
2517 
2518  return
2519  end subroutine atmos_phy_mp_tomita08_terminal_velocity
2520 
2521  !-----------------------------------------------------------------------------
2523  subroutine atmos_phy_mp_tomita08_cloud_fraction( &
2524  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
2525  QTRC, &
2526  mask_criterion, &
2527  cldfrac )
2528  implicit none
2529  integer, intent(in) :: ka, ks, ke
2530  integer, intent(in) :: ia, is, ie
2531  integer, intent(in) :: ja, js, je
2532 
2533  real(rp), intent(in) :: qtrc(ka,ia,ja,qa_mp-1)
2534  real(rp), intent(in) :: mask_criterion
2535 
2536  real(rp), intent(out) :: cldfrac(ka,ia,ja)
2537 
2538  real(rp) :: qhydro
2539  integer :: k, i, j
2540  !---------------------------------------------------------------------------
2541 
2542  !$omp parallel do OMP_SCHEDULE_ &
2543  !$omp private(qhydro)
2544  !$acc kernels copyin(QTRC) copyout(cldfrac)
2545  do j = js, je
2546  do i = is, ie
2547  do k = ks, ke
2548  qhydro = qtrc(k,i,j,i_hyd_qc) + qtrc(k,i,j,i_hyd_qr) &
2549  + qtrc(k,i,j,i_hyd_qi) + qtrc(k,i,j,i_hyd_qs) + qtrc(k,i,j,i_hyd_qg)
2550  cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
2551  enddo
2552  enddo
2553  enddo
2554  !$acc end kernels
2555 
2556  return
2557  end subroutine atmos_phy_mp_tomita08_cloud_fraction
2558 
2559  !-----------------------------------------------------------------------------
2561  subroutine atmos_phy_mp_tomita08_effective_radius( &
2562  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
2563  DENS0, TEMP0, QTRC0, &
2564  Re )
2565  use scale_const, only: &
2566  pi => const_pi, &
2567  dens_w => const_dwatr, &
2568  dens_i => const_dice, &
2569  tem00 => const_tem00
2570  use scale_atmos_hydrometeor, only: &
2571  n_hyd, &
2572  i_hc, &
2573  i_hr, &
2574  i_hi, &
2575  i_hs, &
2576  i_hg
2577  implicit none
2578  integer, intent(in) :: ka, ks, ke
2579  integer, intent(in) :: ia, is, ie
2580  integer, intent(in) :: ja, js, je
2581 
2582  real(rp), intent(in) :: dens0(ka,ia,ja)
2583  real(rp), intent(in) :: temp0(ka,ia,ja)
2584  real(rp), intent(in) :: qtrc0(ka,ia,ja,qa_mp-1)
2585 
2586  real(rp), intent(out) :: re (ka,ia,ja,n_hyd) ! effective radius [cm]
2587  real(rp) :: dens(ka)
2588  real(rp) :: temc(ka)
2589  real(rp) :: qr(ka), qs(ka), qg(ka)
2590  real(rp) :: nc(ka)
2591  real(rp) :: n0r(ka), n0s(ka), n0g(ka)
2592  real(rp) :: rlmdr, rlmds, rlmdg
2593 
2594  real(rp), parameter :: m2cm = 100.0_rp
2595 
2596  !---< Roh and Satoh (2014) >---
2597  real(rp) :: tems, xs2
2598  real(rp) :: coef_at(4), coef_bt(4)
2599  real(rp) :: loga_, b_, nm
2600 
2601  real(rp) :: zerosw, sw
2602  integer :: k, i, j, ih
2603  !---------------------------------------------------------------------------
2604 
2605  !$acc data copyin(DENS0,TEMP0,QTRC0) copyout(Re)
2606 
2607  if ( enable_trawlbg2017 ) then
2608  !$omp parallel do OMP_SCHEDULE_ &
2609  !$omp private(sw)
2610  !$acc kernels
2611  do j = js, je
2612  do i = is, ie
2613  do k = ks, ke
2614  sw = 0.5_rp + sign(0.5_rp, temp0(k,i,j) - 192.0_rp)
2615  re(k,i,j,i_hi) = ( ( 40.0_rp + 0.53_rp * ( temp0(k,i,j) - 192.0_rp ) ) * sw &
2616  + ( 12.0_rp + 28.0_rp * exp( 0.65_rp * ( temp0(k,i,j) - 192.0_rp ) ) ) * ( 1.0_rp-sw ) &
2617  ) * 1.0e-6_rp * m2cm
2618  end do
2619  end do
2620  end do
2621  !$acc end kernels
2622  else
2623  !$omp parallel do OMP_SCHEDULE_
2624  !$acc kernels
2625  do j = js, je
2626  do i = is, ie
2627  do k = ks, ke
2628  re(k,i,j,i_hi) = re_qi * m2cm
2629  end do
2630  end do
2631  end do
2632  !$acc end kernels
2633  end if
2634  !$omp parallel do OMP_SCHEDULE_
2635  !$acc kernels
2636  do ih = i_hg+1, n_hyd
2637  do j = js, je
2638  do i = is, ie
2639  do k = ks, ke
2640  re(k,i,j,ih) = 0.0_rp
2641  end do
2642  end do
2643  end do
2644  end do
2645  !$acc end kernels
2646 
2647  if ( const_rec .or. fixed_re ) then
2648 
2649  !$omp parallel do OMP_SCHEDULE_
2650  !$acc kernels
2651  do j = js, je
2652  do i = is, ie
2653  do k = ks, ke
2654  re(k,i,j,i_hc) = re_qc * m2cm
2655  end do
2656  end do
2657  end do
2658  !$acc end kernels
2659 
2660  else
2661 
2662  !$omp parallel do OMP_SCHEDULE_ &
2663  !$omp private(Nc)
2664  !$acc kernels
2665  do j = js, je
2666  !$acc loop private(Nc)
2667  do i = is, ie
2668  if ( do_couple_aerosol ) then
2669  do k = ks, ke
2670  ! Nc(k) = max( CCN(k,i,j), Nc_def(i,j)*1.E+6_RP ) ! [#/m3] tentatively off the CCN effect
2671  nc(k) = nc_def(i,j) * 1.e+6_rp ! [#/m3]
2672  enddo
2673  else
2674  do k = ks, ke
2675  nc(k) = nc_def(i,j) * 1.e+6_rp ! [#/m3]
2676  enddo
2677  endif
2678 
2679  do k = ks, ke
2680  re(k,i,j,i_hc) = 1.1_rp &
2681  * ( dens0(k,i,j) * qtrc0(k,i,j,i_hyd_qc) / nc(k) / ( 4.0_rp / 3.0_rp * pi * dens_w ) )**(1.0_rp/3.0_rp)
2682  re(k,i,j,i_hc) = min( 1.e-3_rp, max( 1.e-6_rp, re(k,i,j,i_hc) ) ) * m2cm
2683  enddo
2684  enddo
2685  enddo
2686  !$acc end kernels
2687 
2688  endif
2689 
2690  if ( fixed_re ) then
2691 
2692  !$omp parallel do OMP_SCHEDULE_
2693  !$acc kernels
2694  do j = js, je
2695  do i = is, ie
2696  do k = ks, ke
2697  re(k,i,j,i_hr) = 10000.e-6_rp * m2cm
2698  end do
2699  end do
2700  end do
2701  !$acc end kernels
2702 
2703  !$omp parallel do OMP_SCHEDULE_
2704  !$acc kernels
2705  do j = js, je
2706  do i = is, ie
2707  do k = ks, ke
2708  re(k,i,j,i_hs) = re_qi * m2cm
2709  end do
2710  end do
2711  end do
2712  !$acc end kernels
2713 
2714  !$omp parallel do OMP_SCHEDULE_
2715  !$acc kernels
2716  do j = js, je
2717  do i = is, ie
2718  do k = ks, ke
2719  re(k,i,j,i_hg) = re_qi * m2cm
2720  end do
2721  end do
2722  end do
2723  !$acc end kernels
2724 
2725  else
2726 
2727 #ifndef __GFORTRAN__
2728  !$omp parallel do default(none) &
2729  !$omp shared(JS,JE,IS,IE,KS,KE,DENS0,TEMP0,QTRC0,enable_WDXZ2014,N0r_def) &
2730  !$omp shared(N0s_def,N0g_def,Ar,GAM_1br,Re,As,GAM_1bs) &
2731  !$omp shared(enable_RS2014,Ag,GAM_1bg) &
2732  !$omp private(i,j,k,dens,temc,qr,qs,qg,N0r,N0s,N0g,zerosw,RLMDr,RLMDs,Xs2,nm,tems,loga_) &
2733  !$omp private(b_,RLMDg,coef_at,coef_bt) &
2734  !$omp OMP_SCHEDULE_ collapse(2)
2735 #else
2736  !$omp parallel do default(shared) &
2737  !$omp private(i,j,k,dens,temc,qr,qs,qg,N0r,N0s,N0g,zerosw,RLMDr,RLMDs,Xs2,nm,tems,loga_) &
2738  !$omp private(b_,RLMDg,coef_at,coef_bt) &
2739  !$omp OMP_SCHEDULE_ collapse(2)
2740 #endif
2741  !$acc kernels
2742  do j = js, je
2743  !$acc loop private(dens,temc,qr,qs,qg,N0r,N0s,N0g)
2744  do i = is, ie
2745 
2746  do k = ks, ke
2747  dens(k) = dens0(k,i,j)
2748  temc(k) = temp0(k,i,j) - tem00
2749  qr(k) = qtrc0(k,i,j,i_hyd_qr)
2750  qs(k) = qtrc0(k,i,j,i_hyd_qs)
2751  qg(k) = qtrc0(k,i,j,i_hyd_qg)
2752  end do
2753 
2754  ! intercept parameter N0
2755  if ( enable_wdxz2014 ) then
2756  ! Wainwright et al. (2014)
2757  do k = ks, ke
2758  n0r(k) = 1.16e+5_rp * exp( log( max( dens(k)*qr(k)*1000.0_rp, 1.e-2_rp ) )*0.477_rp )
2759  n0s(k) = 4.58e+9_rp * exp( log( max( dens(k)*qs(k)*1000.0_rp, 1.e-2_rp ) )*0.788_rp )
2760  n0g(k) = 9.74e+8_rp * exp( log( max( dens(k)*qg(k)*1000.0_rp, 1.e-2_rp ) )*0.816_rp )
2761  end do
2762  else
2763  do k = ks, ke
2764  n0r(k) = n0r_def ! Marshall and Palmer (1948)
2765  n0s(k) = n0s_def ! Gunn and Marshall (1958)
2766  n0g(k) = n0g_def
2767  end do
2768  end if
2769 
2770  ! slope parameter lambda
2771  do k = ks, ke
2772  zerosw = 0.5_rp - sign(0.5_rp, qr(k) - 1.e-12_rp )
2773  rlmdr = sqrt(sqrt( dens(k) * qr(k) / ( ar * n0r(k) * gam_1br ) + zerosw )) * ( 1.0_rp-zerosw )
2774  ! Effective radius is defined by r3m/r2m = 1.5/lambda
2775  re(k,i,j,i_hr) = 1.5_rp * rlmdr * m2cm
2776  end do
2777 
2778  if ( enable_rs2014 ) then
2779  !$acc loop private(coef_at, coef_bt)
2780  do k = ks, ke
2781  zerosw = 0.5_rp - sign(0.5_rp, qs(k) - 1.e-12_rp )
2782  !---< modification by Roh and Satoh (2014) >---
2783  ! bimodal size distribution of snow
2784  zerosw = 0.5_rp - sign(0.5_rp, dens(k) * qs(k) - 1.e-12_rp )
2785  xs2 = dens(k) * qs(k) / as
2786 
2787  tems = min( -0.1_rp, temc(k) )
2788  coef_at(1) = coef_a01 + tems * ( coef_a02 + tems * ( coef_a05 + tems * coef_a09 ) )
2789  coef_at(2) = coef_a03 + tems * ( coef_a04 + tems * coef_a07 )
2790  coef_at(3) = coef_a06 + tems * coef_a08
2791  coef_at(4) = coef_a10
2792  coef_bt(1) = coef_b01 + tems * ( coef_b02 + tems * ( coef_b05 + tems * coef_b09 ) )
2793  coef_bt(2) = coef_b03 + tems * ( coef_b04 + tems * coef_b07 )
2794  coef_bt(3) = coef_b06 + tems * coef_b08
2795  coef_bt(4) = coef_b10
2796  ! 1 + Bs(=2) moment
2797  nm = 3.0_rp
2798  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
2799  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
2800 
2801  re(k,i,j,i_hs) = 0.5_rp * exp( ln10 * loga_ + log(xs2+zerosw) * b_ ) * ( 1.0_rp-zerosw ) / ( xs2+zerosw ) * m2cm
2802  end do
2803  else
2804  do k = ks, ke
2805  zerosw = 0.5_rp - sign(0.5_rp, qs(k) - 1.e-12_rp )
2806  rlmds = sqrt(sqrt( dens(k) * qs(k) / ( as * n0s(k) * gam_1bs ) + zerosw )) * ( 1.0_rp-zerosw )
2807  re(k,i,j,i_hs) = 1.5_rp * rlmds * m2cm
2808  ! Re(k,i,j,I_HS) = dens * qs / N0s / ( 2.0_RP / 3.0_RP * PI * dens_i ) / RLMDs**3 * m2cm
2809  end do
2810  end if
2811 
2812  do k = ks, ke
2813  zerosw = 0.5_rp - sign(0.5_rp, qg(k) - 1.e-12_rp )
2814  rlmdg = sqrt(sqrt( dens(k) * qg(k) / ( ag * n0g(k) * gam_1bg ) + zerosw )) * ( 1.0_rp-zerosw )
2815  re(k,i,j,i_hg) = 1.5_rp * rlmdg * m2cm
2816  ! Re(k,i,j,I_HG) = dens * qg / N0g / ( 2.0_RP / 3.0_RP * PI * dens_i ) / RLMDg**3 * m2cm
2817  enddo
2818 
2819  enddo
2820  enddo
2821  !$acc end kernels
2822 
2823  endif
2824 
2825  !$acc end data
2826 
2827  return
2828  end subroutine atmos_phy_mp_tomita08_effective_radius
2829 
2830  !-----------------------------------------------------------------------------
2832  subroutine atmos_phy_mp_tomita08_qtrc2qhyd( &
2833  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
2834  QTRC, &
2835  Qe )
2836  use scale_atmos_hydrometeor, only: &
2837  n_hyd, &
2838  i_hc, &
2839  i_hr, &
2840  i_hi, &
2841  i_hs, &
2842  i_hg
2843  implicit none
2844  integer, intent(in) :: ka, ks, ke
2845  integer, intent(in) :: ia, is, ie
2846  integer, intent(in) :: ja, js, je
2847 
2848  real(rp), intent(in) :: qtrc(ka,ia,ja,qa_mp-1)
2849 
2850  real(rp), intent(out) :: qe(ka,ia,ja,n_hyd) ! mass ratio of each cateory [kg/kg]
2851 
2852  integer :: k, i, j, ih
2853  !---------------------------------------------------------------------------
2854 
2855  !$acc data copyin(QTRC) copyout(Qe)
2856 
2857  !$omp parallel do OMP_SCHEDULE_
2858  !$acc kernels async
2859 !OCL XFILL
2860  do j = js, je
2861  do i = is, ie
2862  do k = ks, ke
2863  qe(k,i,j,i_hc) = qtrc(k,i,j,i_hyd_qc)
2864  end do
2865  end do
2866  end do
2867  !$acc end kernels
2868 
2869  !$omp parallel do OMP_SCHEDULE_
2870 !OCL XFILL
2871  !$acc kernels async
2872  do j = js, je
2873  do i = is, ie
2874  do k = ks, ke
2875  qe(k,i,j,i_hr) = qtrc(k,i,j,i_hyd_qr)
2876  end do
2877  end do
2878  end do
2879  !$acc end kernels
2880 
2881  !$omp parallel do OMP_SCHEDULE_
2882  !$acc kernels async
2883 !OCL XFILL
2884  do j = js, je
2885  do i = is, ie
2886  do k = ks, ke
2887  qe(k,i,j,i_hi) = qtrc(k,i,j,i_hyd_qi)
2888  end do
2889  end do
2890  end do
2891  !$acc end kernels
2892 
2893  !$omp parallel do OMP_SCHEDULE_
2894  !$acc kernels async
2895 !OCL XFILL
2896  do j = js, je
2897  do i = is, ie
2898  do k = ks, ke
2899  qe(k,i,j,i_hs) = qtrc(k,i,j,i_hyd_qs)
2900  end do
2901  end do
2902  end do
2903  !$acc end kernels
2904 
2905  !$omp parallel do OMP_SCHEDULE_
2906  !$acc kernels async
2907 !OCL XFILL
2908  do j = js, je
2909  do i = is, ie
2910  do k = ks, ke
2911  qe(k,i,j,i_hg) = qtrc(k,i,j,i_hyd_qg)
2912  end do
2913  end do
2914  end do
2915  !$acc end kernels
2916 
2917  !$omp parallel do OMP_SCHEDULE_
2918  !$acc kernels async
2919 !OCL XFILL
2920  do ih = i_hg+1, n_hyd
2921  do j = js, je
2922  do i = is, ie
2923  do k = ks, ke
2924  qe(k,i,j,ih) = 0.0_rp
2925  end do
2926  end do
2927  end do
2928  end do
2929  !$acc end kernels
2930 
2931  !$acc wait
2932 
2933  !$acc end data
2934 
2935  return
2936  end subroutine atmos_phy_mp_tomita08_qtrc2qhyd
2937 
2938  !-----------------------------------------------------------------------------
2940  subroutine atmos_phy_mp_tomita08_qhyd2qtrc( &
2941  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
2942  Qe, &
2943  QTRC )
2944  use scale_atmos_hydrometeor, only: &
2945  n_hyd, &
2946  i_hc, &
2947  i_hr, &
2948  i_hi, &
2949  i_hs, &
2950  i_hg, &
2951  i_hh
2952  implicit none
2953 
2954  integer, intent(in) :: ka, ks, ke
2955  integer, intent(in) :: ia, is, ie
2956  integer, intent(in) :: ja, js, je
2957  real(rp), intent(in) :: qe (ka,ia,ja,n_hyd) ! mass ratio of each cateory [kg/kg]
2958  real(rp), intent(out) :: qtrc(ka,ia,ja,qa_mp-1)
2959 
2960  integer :: k, i, j
2961  !---------------------------------------------------------------------------
2962 
2963  !$acc data copyin(Qe) copyout(QTRC)
2964 
2965  !$omp parallel do OMP_SCHEDULE_
2966  !$acc kernels async
2967 !OCL XFILL
2968  do j = js, je
2969  do i = is, ie
2970  do k = ks, ke
2971  qtrc(k,i,j,i_hyd_qc) = qe(k,i,j,i_hc)
2972  end do
2973  end do
2974  end do
2975  !$acc end kernels
2976 
2977  !$omp parallel do OMP_SCHEDULE_
2978  !$acc kernels async
2979 !OCL XFILL
2980  do j = js, je
2981  do i = is, ie
2982  do k = ks, ke
2983  qtrc(k,i,j,i_hyd_qr) = qe(k,i,j,i_hr)
2984  end do
2985  end do
2986  end do
2987  !$acc end kernels
2988 
2989  !$omp parallel do OMP_SCHEDULE_
2990  !$acc kernels async
2991 !OCL XFILL
2992  do j = js, je
2993  do i = is, ie
2994  do k = ks, ke
2995  qtrc(k,i,j,i_hyd_qi) = qe(k,i,j,i_hi)
2996  end do
2997  end do
2998  end do
2999  !$acc end kernels
3000 
3001  !$omp parallel do OMP_SCHEDULE_
3002  !$acc kernels async
3003 !OCL XFILL
3004  do j = js, je
3005  do i = is, ie
3006  do k = ks, ke
3007  qtrc(k,i,j,i_hyd_qs) = qe(k,i,j,i_hs)
3008  end do
3009  end do
3010  end do
3011  !$acc end kernels
3012 
3013  !$omp parallel do OMP_SCHEDULE_
3014  !$acc kernels async
3015 !OCL XFILL
3016  do j = js, je
3017  do i = is, ie
3018  do k = ks, ke
3019  qtrc(k,i,j,i_hyd_qg) = qe(k,i,j,i_hg) + qe(k,i,j,i_hh)
3020  end do
3021  end do
3022  end do
3023  !$acc end kernels
3024 
3025  !$acc wait
3026 
3027  !$acc end data
3028 
3029  return
3030  end subroutine atmos_phy_mp_tomita08_qhyd2qtrc
3031 
3032 end module scale_atmos_phy_mp_tomita08
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:49
scale_const::const_lhv0
real(rp), parameter, public const_lhv0
latent heat of vaporizaion at 0C [J/kg]
Definition: scale_const.F90:82
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_putvar
subroutine, public atmos_phy_mp_tomita08_putvar(in_drag_g, in_re_qc, in_re_qi, in_Cr, in_Cs)
overwrite private variables
Definition: scale_atmos_phy_mp_tomita08.F90:671
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_atmos_hydrometeor::i_hr
integer, parameter, public i_hr
liquid water rain
Definition: scale_atmos_hydrometeor.F90:98
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_qhyd2qtrc
subroutine, public atmos_phy_mp_tomita08_qhyd2qtrc(KA, KS, KE, IA, IS, IE, JA, JS, JE, Qe, QTRC)
get mass ratio of each category
Definition: scale_atmos_phy_mp_tomita08.F90:2944
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_ntracers
integer, parameter, public atmos_phy_mp_tomita08_ntracers
Definition: scale_atmos_phy_mp_tomita08.F90:46
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_terminal_velocity
subroutine, public atmos_phy_mp_tomita08_terminal_velocity(KA, KS, KE, DENS, TEMP, RHOQ, vterm)
Lin-type cold rain microphysics (terminal velocity)
Definition: scale_atmos_phy_mp_tomita08.F90:2257
scale_atmos_hydrometeor::cp_water
real(rp), public cp_water
CP for water [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:152
scale_atmos_hydrometeor::i_hs
integer, parameter, public i_hs
snow
Definition: scale_atmos_hydrometeor.F90:100
scale_atmos_phy_mp_tomita08
module atmosphere / physics / microphysics / Tomita08
Definition: scale_atmos_phy_mp_tomita08.F90:13
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_nwaters
integer, parameter, public atmos_phy_mp_tomita08_nwaters
Definition: scale_atmos_phy_mp_tomita08.F90:47
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_tracer_units
character(len=h_short), dimension(qa_mp), parameter, public atmos_phy_mp_tomita08_tracer_units
Definition: scale_atmos_phy_mp_tomita08.F90:63
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_finalize
subroutine, public atmos_phy_mp_tomita08_finalize
finalize
Definition: scale_atmos_phy_mp_tomita08.F90:656
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:68
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:35
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level, disable_barrier)
Start raptime.
Definition: scale_prof.F90:174
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_setup
subroutine, public atmos_phy_mp_tomita08_setup(KA, KS, KE, IA, IS, IE, JA, JS, JE, flg_lt)
ATMOS_PHY_MP_tomita08_setup Setup.
Definition: scale_atmos_phy_mp_tomita08.F90:399
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_getvar
subroutine, public atmos_phy_mp_tomita08_getvar(out_drag_g, out_re_qc, out_re_qi, out_Cr, out_Cs)
read private variables
Definition: scale_atmos_phy_mp_tomita08.F90:703
scale_atmos_hydrometeor::i_hh
integer, parameter, public i_hh
hail
Definition: scale_atmos_hydrometeor.F90:102
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_adjustment
subroutine, public atmos_phy_mp_tomita08_adjustment(KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, PRES, CCN, dt, TEMP, QTRC, CPtot, CVtot, RHOE_t, EVAPORATE, flg_lt, d0_crg, v0_crg, dqcrg, beta_crg, QSPLT_in, Sarea, QTRC_crg)
ATMOS_PHY_MP_tomita08_adjustment calculate state after saturation process.
Definition: scale_atmos_phy_mp_tomita08.F90:735
scale_file_history
module file_history
Definition: scale_file_history.F90:15
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_atmos_hydrometeor::cv_vapor
real(rp), public cv_vapor
CV for vapor [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:149
scale_atmos_hydrometeor::i_hi
integer, parameter, public i_hi
ice water cloud
Definition: scale_atmos_hydrometeor.F90:99
scale_io
module STDIO
Definition: scale_io.F90:10
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_qtrc2qhyd
subroutine, public atmos_phy_mp_tomita08_qtrc2qhyd(KA, KS, KE, IA, IS, IE, JA, JS, JE, QTRC, Qe)
Calculate mass ratio of each category.
Definition: scale_atmos_phy_mp_tomita08.F90:2836
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_hydrometeor::lhv
real(rp), public lhv
latent heat of vaporization for use [J/kg]
Definition: scale_atmos_hydrometeor.F90:144
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_effective_radius
subroutine, public atmos_phy_mp_tomita08_effective_radius(KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS0, TEMP0, QTRC0, Re)
Calculate Effective Radius.
Definition: scale_atmos_phy_mp_tomita08.F90:2565
scale_atmos_phy_mp_common
module ATMOSPHERE / Physics Cloud Microphysics - Common
Definition: scale_atmos_phy_mp_common.F90:13
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_tracer_names
character(len=h_short), dimension(qa_mp), parameter, public atmos_phy_mp_tomita08_tracer_names
Definition: scale_atmos_phy_mp_tomita08.F90:49
scale_specfunc
module SPECFUNC
Definition: scale_specfunc.F90:14
scale_const::const_pi
real(rp), parameter, public const_pi
pi
Definition: scale_const.F90:32
scale_atmos_hydrometeor::i_hc
integer, parameter, public i_hc
liquid water cloud
Definition: scale_atmos_hydrometeor.F90:97
scale_const::const_dwatr
real(rp), parameter, public const_dwatr
density of water [kg/m3]
Definition: scale_const.F90:89
scale_const::const_tem00
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
Definition: scale_const.F90:99
scale_const::const_cl
real(rp), parameter, public const_cl
specific heat (liquid water) [J/kg/K]
Definition: scale_const.F90:71
scale_const::const_lhf0
real(rp), public const_lhf0
latent heat of fusion at 0C [J/kg]
Definition: scale_const.F90:86
scale_atmos_hydrometeor::lhf
real(rp), public lhf
latent heat of fusion for use [J/kg]
Definition: scale_atmos_hydrometeor.F90:146
scale_file_history::file_history_reg
subroutine, public file_history_reg(name, desc, unit, itemid, standard_name, ndims, dim_type, cell_measures, fill_halo)
Register/Append variable to history file.
Definition: scale_file_history.F90:685
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level, disable_barrier)
Save raptime.
Definition: scale_prof.F90:246
scale_const::const_lhs0
real(rp), parameter, public const_lhs0
latent heat of sublimation at 0C [J/kg]
Definition: scale_const.F90:84
scale_const::const_dice
real(rp), parameter, public const_dice
density of ice [kg/m3]
Definition: scale_const.F90:90
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_nices
integer, parameter, public atmos_phy_mp_tomita08_nices
Definition: scale_atmos_phy_mp_tomita08.F90:48
scale_specfunc::sf_gamma
real(rp) function, public sf_gamma(x)
Gamma function.
Definition: scale_specfunc.F90:50
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_cloud_fraction
subroutine, public atmos_phy_mp_tomita08_cloud_fraction(KA, KS, KE, IA, IS, IE, JA, JS, JE, QTRC, mask_criterion, cldfrac)
Calculate Cloud Fraction.
Definition: scale_atmos_phy_mp_tomita08.F90:2528
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:97
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:57
scale_atmos_hydrometeor::n_hyd
integer, parameter, public n_hyd
Definition: scale_atmos_hydrometeor.F90:95
scale_atmos_hydrometeor::cp_vapor
real(rp), public cp_vapor
CP for vapor [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:150
scale_atmos_hydrometeor::cv_water
real(rp), public cv_water
CV for water [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:151
scale_atmos_hydrometeor::cv_ice
real(rp), public cv_ice
CV for ice [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:153
scale_atmos_hydrometeor::cp_ice
real(rp), public cp_ice
CP for ice [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:154
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_tracer_descriptions
character(len=h_mid), dimension(qa_mp), parameter, public atmos_phy_mp_tomita08_tracer_descriptions
Definition: scale_atmos_phy_mp_tomita08.F90:56
scale_atmos_hydrometeor::i_hg
integer, parameter, public i_hg
graupel
Definition: scale_atmos_hydrometeor.F90:101
Generated by   doxygen 1.8.17