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