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