SCALE-RM
scale_atmos_phy_mp_tomita08.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
15 !-------------------------------------------------------------------------------
16 module scale_atmos_phy_mp_tomita08
17  !-----------------------------------------------------------------------------
18  !
19  !++ used modules
20  !
21  use scale_precision
22  use scale_stdio
23  use scale_prof
24  use scale_grid_index
25 
26  use scale_tracer_tomita08
27  use scale_tracer, only: qa
28  !-----------------------------------------------------------------------------
29  implicit none
30  private
31  !-----------------------------------------------------------------------------
32  !
33  !++ Public procedure
34  !
35  public :: atmos_phy_mp_tomita08_setup
36  public :: atmos_phy_mp_tomita08
37  public :: atmos_phy_mp_tomita08_cloudfraction
38  public :: atmos_phy_mp_tomita08_effectiveradius
39  public :: atmos_phy_mp_tomita08_mixingratio
40 
41  !-----------------------------------------------------------------------------
42  !
43  !++ Public parameters & variables
44  !
45  real(RP), public, target :: atmos_phy_mp_dens(mp_qa) ! hydrometeor density [kg/m3]=[g/L]
46 
47  !-----------------------------------------------------------------------------
48  !
49  !++ Private procedure
50  !
51  private :: mp_tomita08
52  private :: mp_tomita08_vterm
53  private :: mp_tomita08_bergeronparam
54 
55  !-----------------------------------------------------------------------------
56  !
57  !++ Private parameters & variables
58  !
59  logical, private :: mp_donegative_fixer = .true. ! apply negative fixer?
60  logical, private :: mp_doprecipitation = .true. ! apply sedimentation (precipitation)?
61  logical, private :: mp_couple_aerosol = .false. ! apply CCN effect?
62 
63  logical, private :: mp_doexpricit_icegen = .false. ! apply explicit ice generation?
64  logical, private :: only_liquid = .false.
65  real(RP), private :: sw_useicegen = 0.0_rp
66 
67  real(RP), private :: dens00 = 1.28_rp
68 
69  ! Parameter for Marshall-Palmer distribution
70  real(RP), private :: n0r = 8.e6_rp
71  real(RP), private :: n0s = 3.e6_rp
72  real(RP), private :: n0g = 4.e6_rp
73 
74  real(RP), private :: dens_s = 100.0_rp
75  real(RP), private :: dens_g = 400.0_rp
76  ! graupel : 400
77  ! hail : 917
78  real(RP), private :: drag_g = 0.6_rp
79  real(RP), private :: re_qc = 8.e-6_rp ! effective radius for cloud water
80  real(RP), private :: re_qi = 40.e-6_rp ! effective radius for cloud ice
81  real(RP), private :: cr = 130.0_rp
82  real(RP), private :: cs = 4.84_rp
83 
84  ! Empirical parameter
85  real(RP), private :: ar, as, ag
86  real(RP), private :: br, bs, bg
87  real(RP), private :: cg
88  real(RP), private :: dr, ds, dg
89 
90  ! GAMMA function
91  real(RP), private :: gam, gam_2, gam_3
92 
93  real(RP), private :: gam_1br, gam_2br, gam_3br
94  real(RP), private :: gam_3dr
95  real(RP), private :: gam_6dr
96  real(RP), private :: gam_1brdr
97  real(RP), private :: gam_5dr_h
98 
99  real(RP), private :: gam_1bs, gam_2bs, gam_3bs
100  real(RP), private :: gam_3ds
101  real(RP), private :: gam_1bsds
102  real(RP), private :: gam_5ds_h
103 
104  real(RP), private :: gam_1bg, gam_3dg
105  real(RP), private :: gam_1bgdg
106  real(RP), private :: gam_5dg_h
107 
108  !---< Khairoutdinov and Kogan (2000) >---
109  real(RP), private :: sw_kk2000 = 0.0_rp
110 
111  !---< Roh and Satoh (2014) >---
112  real(RP), private :: sw_roh2014 = 0.0_rp
113  real(RP), private, parameter :: coef_a(10) = (/ 5.065339_rp, -0.062659_rp, -3.032362_rp, 0.029469_rp, -0.000285_rp, &
114  0.31255_rp, 0.000204_rp, 0.003199_rp, 0.0_rp, -0.015952_rp /)
115  real(RP), private, parameter :: coef_b(10) = (/ 0.476221_rp, -0.015896_rp, 0.165977_rp, 0.007468_rp, -0.000141_rp, &
116  0.060366_rp, 0.000079_rp, 0.000594_rp, 0.0_rp, -0.003577_rp /)
117 
118  ! Accretion parameter
119  real(RP), private :: eiw = 1.0_rp
120  real(RP), private :: erw = 1.0_rp
121  real(RP), private :: esw = 1.0_rp
122  real(RP), private :: egw = 1.0_rp
123  real(RP), private :: eri = 1.0_rp
124  real(RP), private :: esi = 1.0_rp
125  real(RP), private :: egi = 0.1_rp
126  real(RP), private :: esr = 1.0_rp
127  real(RP), private :: egr = 1.0_rp
128  real(RP), private :: egs = 1.0_rp
129  real(RP), private :: gamma_sacr = 25.e-3_rp
130  real(RP), private :: gamma_gacs = 90.e-3_rp
131 
132  ! Auto-conversion parameter
133  real(RP), private, parameter :: nc_lnd = 2000.0_rp
134  real(RP), private, parameter :: nc_ocn = 50.0_rp
135  real(RP), private, allocatable :: nc_def(:,:)
136 
137  real(RP), private :: beta_saut = 6.e-3_rp
138  real(RP), private :: gamma_saut = 60.e-3_rp
139  real(RP), private :: beta_gaut = 1.e-3_rp
140  real(RP), private :: gamma_gaut = 90.e-3_rp
141  real(RP), private :: qicrt_saut = 0.0_rp
142  real(RP), private :: qscrt_gaut = 6.e-4_rp
143 
144  ! Evaporation, Sublimation parameter
145  real(RP), private, parameter :: da0 = 2.428e-2_rp
146  real(RP), private, parameter :: dda_dt = 7.47e-5_rp
147  real(RP), private, parameter :: dw0 = 2.222e-5_rp
148  real(RP), private, parameter :: ddw_dt = 1.37e-7_rp
149  real(RP), private, parameter :: mu0 = 1.718e-5_rp
150  real(RP), private, parameter :: dmu_dt = 5.28e-8_rp
151 
152  real(RP), private :: f1r = 0.78_rp
153  real(RP), private :: f2r = 0.27_rp
154  real(RP), private :: f1s = 0.65_rp
155  real(RP), private :: f2s = 0.39_rp
156  real(RP), private :: f1g = 0.78_rp
157  real(RP), private :: f2g = 0.27_rp
158 
159  ! Freezing parameter
160  real(RP), private :: a_frz = 0.66_rp
161  real(RP), private :: b_frz = 100.0_rp
162 
163  ! Bergeron process parameter
164  real(RP), private :: mi40 = 2.46e-10_rp
165  real(RP), private :: mi50 = 4.80e-10_rp
166  real(RP), private :: vti50 = 1.0_rp
167  real(RP), private :: ri50 = 5.e-5_rp
168 
169  integer, private, parameter :: w_nmax = 42
170  integer, private, parameter :: i_dqv_dt = 1 !
171  integer, private, parameter :: i_dqc_dt = 2 !
172  integer, private, parameter :: i_dqr_dt = 3 !
173  integer, private, parameter :: i_dqi_dt = 4 !
174  integer, private, parameter :: i_dqs_dt = 5 !
175  integer, private, parameter :: i_dqg_dt = 6 !
176  integer, private, parameter :: i_delta1 = 7 ! separation switch for r->s,g
177  integer, private, parameter :: i_delta2 = 8 ! separation switch for s->g
178  integer, private, parameter :: i_delta3 = 9 ! separation switch for ice sublimation
179  integer, private, parameter :: i_rlmdr = 10
180  integer, private, parameter :: i_rlmds = 11
181  integer, private, parameter :: i_rlmdg = 12
182  integer, private, parameter :: i_piacr = 13 ! r->s,g
183  integer, private, parameter :: i_psacr = 14 ! r->s,g
184  integer, private, parameter :: i_praci = 15 ! i->s,g
185  integer, private, parameter :: i_psmlt = 16 ! s->r
186  integer, private, parameter :: i_pgmlt = 17 ! g->r
187  integer, private, parameter :: i_praut = 18 ! c->r
188  integer, private, parameter :: i_pracw = 19 ! c->r
189  integer, private, parameter :: i_psacw = 20 ! c->s
190  integer, private, parameter :: i_psfw = 21 ! c->s
191  integer, private, parameter :: i_pgacw = 22 ! c->g
192  integer, private, parameter :: i_prevp = 23 ! r->v
193  integer, private, parameter :: i_piacr_s = 24 ! r->s
194  integer, private, parameter :: i_psacr_s = 25 ! r->s
195  integer, private, parameter :: i_piacr_g = 26 ! r->g
196  integer, private, parameter :: i_psacr_g = 27 ! r->g
197  integer, private, parameter :: i_pgacr = 28 ! r->g
198  integer, private, parameter :: i_pgfrz = 29 ! r->g
199  integer, private, parameter :: i_psaut = 30 ! i->s
200  integer, private, parameter :: i_praci_s = 31 ! i->s
201  integer, private, parameter :: i_psaci = 32 ! i->s
202  integer, private, parameter :: i_psfi = 33 ! i->s
203  integer, private, parameter :: i_praci_g = 34 ! i->g
204  integer, private, parameter :: i_pgaci = 35 ! i->g
205  integer, private, parameter :: i_psdep = 36 ! v->s
206  integer, private, parameter :: i_pssub = 37 ! s->v
207  integer, private, parameter :: i_pgaut = 38 ! s->g
208  integer, private, parameter :: i_pracs = 39 ! s->g
209  integer, private, parameter :: i_pgacs = 40 ! s->g
210  integer, private, parameter :: i_pgdep = 41 ! v->g
211  integer, private, parameter :: i_pgsub = 42 ! g->v
212 
213  integer, private :: w_histid (w_nmax) = -999
214  logical, private :: w_zinterp(w_nmax) = .false.
215  character(len=H_SHORT), private :: w_name (w_nmax)
216  character(len=H_MID), private :: w_desc (w_nmax) = ''
217  character(len=H_SHORT), private :: w_unit (w_nmax) = 'kg/kg/s'
218 
219  data w_name / 'dqv_dt ', &
220  'dqc_dt ', &
221  'dqr_dt ', &
222  'dqi_dt ', &
223  'dqs_dt ', &
224  'dqg_dt ', &
225  'delta1 ', &
226  'delta2 ', &
227  'delta3 ', &
228  'RLMDr ', &
229  'RLMDs ', &
230  'RLMDg ', &
231  'Piacr ', &
232  'Psacr ', &
233  'Praci ', &
234  'Psmlt ', &
235  'Pgmlt ', &
236  'Praut ', &
237  'Pracw ', &
238  'Psacw ', &
239  'Psfw ', &
240  'Pgacw ', &
241  'Prevp ', &
242  'Piacr_s', &
243  'Psacr_s', &
244  'Piacr_g', &
245  'Psacr_g', &
246  'Pgacr ', &
247  'Pgfrz ', &
248  'Psaut ', &
249  'Praci_s', &
250  'Psaci ', &
251  'Psfi ', &
252  'Praci_g', &
253  'Pgaci ', &
254  'Psdep ', &
255  'Pssub ', &
256  'Pgaut ', &
257  'Pracs ', &
258  'Pgacs ', &
259  'Pgdep ', &
260  'Pgsub ' /
261 
262  real(RP), private, allocatable :: work3d(:,:,:)
263 
264  integer, private :: mp_ntmax_sedimentation = 1 ! number of time step for sedimentation
265 
266  integer, private :: mp_nstep_sedimentation
267  real(RP), private :: mp_rnstep_sedimentation
268  real(DP), private :: mp_dtsec_sedimentation
269 
270  logical, private :: debug
271 
272  !-----------------------------------------------------------------------------
273 contains
274  !-----------------------------------------------------------------------------
276  subroutine atmos_phy_mp_tomita08_setup( MP_TYPE )
277  use scale_process, only: &
278  prc_mpistop
279  use scale_const, only: &
280  pi => const_pi, &
281  grav => const_grav, &
282  dens_w => const_dwatr, &
283  dens_i => const_dice
284  use scale_specfunc, only: &
285  sf_gamma
286  use scale_time, only: &
287  time_dtsec_atmos_phy_mp
288  use scale_grid, only: &
289  cdz => grid_cdz
290  use scale_history, only: &
291  hist_reg
292  implicit none
293 
294  character(len=*), intent(in) :: MP_TYPE
295 
296  real(RP) :: autoconv_nc = nc_ocn
297  logical :: enable_kk2000 = .false.
298  logical :: enable_roh2014 = .false.
299 
300  namelist / param_atmos_phy_mp / &
301  mp_doprecipitation, &
302  mp_donegative_fixer, &
303  mp_doexpricit_icegen, &
304  mp_ntmax_sedimentation, &
305  mp_couple_aerosol
306 
307  namelist / param_atmos_phy_mp_tomita08 / &
308  autoconv_nc, &
309  enable_kk2000, &
310  enable_roh2014, &
311  dens_s, &
312  dens_g, &
313  re_qc, &
314  re_qi, &
315  cr, &
316  cs, &
317  drag_g, &
318  beta_saut, &
319  gamma_saut, &
320  gamma_sacr, &
321  n0s, &
322  n0g, &
323  debug
324 
325  real(RP), parameter :: max_term_vel = 10.0_rp !-- terminal velocity for calculate dt of sedimentation
326  integer :: nstep_max
327 
328  integer :: ierr
329  integer :: i, j, ip
330  !---------------------------------------------------------------------------
331 
332  if( io_l ) write(io_fid_log,*)
333  if( io_l ) write(io_fid_log,*) '++++++ Module[Cloud Microphysics] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
334  if( io_l ) write(io_fid_log,*) '*** TOMITA08: 1-moment bulk 6 category'
335 
336  if ( mp_type /= 'TOMITA08' ) then
337  write(*,*) 'xxx ATMOS_PHY_MP_TYPE is not TOMITA08. Check!'
338  call prc_mpistop
339  endif
340 
341  if ( i_qv <= 0 &
342  .OR. i_qc <= 0 &
343  .OR. i_qr <= 0 &
344  .OR. i_qi <= 0 &
345  .OR. i_qs <= 0 &
346  .OR. i_qg <= 0 ) then
347  write(*,*) 'xxx TOMITA08 needs QV/C/R/I/S/G tracer. Check!'
348  call prc_mpistop
349  endif
350 
351  allocate( work3d(ka,ia,ja) )
352  work3d(:,:,:) = 0.0_rp
353 
354  allocate( nc_def(ia,ja) )
355 
356  !--- read namelist
357  rewind(io_fid_conf)
358  read(io_fid_conf,nml=param_atmos_phy_mp,iostat=ierr)
359  if( ierr < 0 ) then !--- missing
360  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
361  elseif( ierr > 0 ) then !--- fatal error
362  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP. Check!'
363  call prc_mpistop
364  endif
365  if( io_lnml ) write(io_fid_log,nml=param_atmos_phy_mp)
366 
367  if ( io_l ) write(io_fid_log,*)
368  if ( io_l ) write(io_fid_log,*) '*** Enable negative fixer? : ', mp_donegative_fixer
369  if ( io_l ) write(io_fid_log,*) '*** Enable sedimentation (precipitation)? : ', mp_doprecipitation
370  if ( io_l ) write(io_fid_log,*) '*** Enable explicit ice generation? : ', mp_doexpricit_icegen
371 
372  nstep_max = ceiling( ( time_dtsec_atmos_phy_mp * max_term_vel ) / minval( cdz(:) ) )
373  mp_ntmax_sedimentation = max( mp_ntmax_sedimentation, nstep_max )
374 
375  mp_nstep_sedimentation = mp_ntmax_sedimentation
376  mp_rnstep_sedimentation = 1.0_rp / real(mp_ntmax_sedimentation,kind=rp)
377  mp_dtsec_sedimentation = time_dtsec_atmos_phy_mp * mp_rnstep_sedimentation
378 
379  if ( io_l ) write(io_fid_log,*) '*** Timestep of sedimentation is divided into : ', mp_ntmax_sedimentation, 'step'
380  if ( io_l ) write(io_fid_log,*) '*** DT of sedimentation : ', mp_dtsec_sedimentation, '[s]'
381 
382  !--- read namelist
383  rewind(io_fid_conf)
384  read(io_fid_conf,nml=param_atmos_phy_mp_tomita08,iostat=ierr)
385  if( ierr < 0 ) then !--- missing
386  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
387  elseif( ierr > 0 ) then !--- fatal error
388  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP_TOMITA08. Check!'
389  call prc_mpistop
390  endif
391  if( io_lnml ) write(io_fid_log,nml=param_atmos_phy_mp_tomita08)
392 
393  if ( io_l ) write(io_fid_log,*)
394  if ( io_l ) write(io_fid_log,*) '*** density of the snow [kg/m3] : ', dens_s
395  if ( io_l ) write(io_fid_log,*) '*** density of the graupel [kg/m3] : ', dens_g
396  if ( io_l ) write(io_fid_log,*) '*** Nc for auto-conversion [num/m3]: ', autoconv_nc
397  if ( io_l ) write(io_fid_log,*) '*** Use k-k scheme? : ', enable_kk2000
398  if ( io_l ) write(io_fid_log,*) '*** Use Roh scheme? : ', enable_roh2014
399  if ( io_l ) write(io_fid_log,*)
400 
401  atmos_phy_mp_dens(i_mp_qc) = dens_w
402  atmos_phy_mp_dens(i_mp_qr) = dens_w
403  atmos_phy_mp_dens(i_mp_qi) = dens_i
404  atmos_phy_mp_dens(i_mp_qs) = dens_s
405  atmos_phy_mp_dens(i_mp_qg) = dens_g
406 
407  !--- empirical coefficients A, B, C, D
408  ar = pi * dens_w / 6.0_rp
409  as = pi * dens_s / 6.0_rp
410  ag = pi * dens_g / 6.0_rp
411 
412  br = 3.0_rp
413  bs = 3.0_rp
414  bg = 3.0_rp
415 
416  cg = sqrt( ( 4.0_rp * dens_g * grav ) / ( 3.0_rp * dens00 * drag_g ) )
417 
418  dr = 0.50_rp
419  ds = 0.25_rp
420  dg = 0.50_rp
421 
422  if ( enable_roh2014 ) then ! overwrite parameters
423  sw_roh2014 = 1.0_rp
424  n0g = 4.e8_rp
425  as = 0.069_rp
426  bs = 2.0_rp
427  esi = 0.25_rp
428  endif
429 
430  gam = 1.0_rp ! =0!
431  gam_2 = 1.0_rp ! =1!
432  gam_3 = 2.0_rp ! =2!
433 
434  gam_1br = sf_gamma( 1.0_rp + br ) ! = 4!
435  gam_2br = sf_gamma( 2.0_rp + br ) ! = 5!
436  gam_3br = sf_gamma( 3.0_rp + br ) ! = 6!
437  gam_3dr = sf_gamma( 3.0_rp + dr )
438  gam_6dr = sf_gamma( 6.0_rp + dr )
439  gam_1brdr = sf_gamma( 1.0_rp + br + dr )
440  gam_5dr_h = sf_gamma( 0.5_rp * (5.0_rp+dr) )
441 
442  gam_1bs = sf_gamma( 1.0_rp + bs ) ! = 4!
443  gam_2bs = sf_gamma( 2.0_rp + bs ) ! = 5!
444  gam_3bs = sf_gamma( 3.0_rp + bs ) ! = 6!
445  gam_3ds = sf_gamma( 3.0_rp + ds )
446  gam_1bsds = sf_gamma( 1.0_rp + bs + ds )
447  gam_5ds_h = sf_gamma( 0.5_rp * (5.0_rp+ds) )
448 
449  gam_1bg = sf_gamma( 1.0_rp + bg ) ! = 4!
450  gam_3dg = sf_gamma( 3.0_rp + dg )
451  gam_1bgdg = sf_gamma( 1.0_rp + bg + dg)
452  gam_5dg_h = sf_gamma( 0.5_rp * (5.0_rp+dg) )
453 
454  do j = js, je
455  do i = is, ie
456  nc_def(i,j) = autoconv_nc
457  enddo
458  enddo
459 
460  if ( mp_doexpricit_icegen ) then
461  only_liquid = .true.
462  sw_useicegen = 1.0_rp
463  else
464  only_liquid = .false.
465  sw_useicegen = 0.0_rp
466  endif
467 
468  if ( enable_kk2000 ) then
469  sw_kk2000 = 1.0_rp
470  else
471  sw_kk2000 = 0.0_rp
472  endif
473 
474  ! detailed tendency monitor
475  w_desc(:) = w_name(:)
476  do ip = 1, w_nmax
477  call hist_reg( w_histid(ip), w_zinterp(ip), w_name(ip), w_desc(ip), w_unit(ip), ndim=3 )
478  enddo
479 
480  return
481  end subroutine atmos_phy_mp_tomita08_setup
482 
483  !-----------------------------------------------------------------------------
485  subroutine atmos_phy_mp_tomita08( &
486  DENS, &
487  MOMZ, &
488  MOMX, &
489  MOMY, &
490  RHOT, &
491  QTRC, &
492  CCN, &
493  EVAPORATE, &
494  SFLX_rain, &
495  SFLX_snow )
496  use scale_grid_index
497  use scale_const, only: &
498  dwatr => const_dwatr, &
499  pi => const_pi
500  use scale_time, only: &
501  dt => time_dtsec_atmos_phy_mp
502  use scale_history, only: &
503  hist_in
504  use scale_tracer, only: &
505  qad => qa
506  use scale_atmos_thermodyn, only: &
507  thermodyn_rhoe => atmos_thermodyn_rhoe, &
508  thermodyn_rhot => atmos_thermodyn_rhot, &
509  thermodyn_temp_pres_e => atmos_thermodyn_temp_pres_e
510  use scale_atmos_phy_mp_common, only: &
511  mp_negative_fixer => atmos_phy_mp_negative_fixer, &
512  mp_precipitation => atmos_phy_mp_precipitation, &
513  mp_saturation_adjustment => atmos_phy_mp_saturation_adjustment
514  implicit none
515 
516  real(RP), intent(inout) :: dens(ka,ia,ja)
517  real(RP), intent(inout) :: MOMZ(ka,ia,ja)
518  real(RP), intent(inout) :: MOMX(ka,ia,ja)
519  real(RP), intent(inout) :: MOMY(ka,ia,ja)
520  real(RP), intent(inout) :: RHOT(ka,ia,ja)
521  real(RP), intent(inout) :: QTRC(ka,ia,ja,qad)
522  real(RP), intent(in) :: CCN(ka,ia,ja)
523  real(RP), intent(out) :: EVAPORATE(ka,ia,ja) ! number of evaporated cloud [/m3/s]
524  real(RP), intent(out) :: SFLX_rain(ia,ja)
525  real(RP), intent(out) :: SFLX_snow(ia,ja)
526 
527  real(RP) :: RHOE_t(ka,ia,ja)
528  real(RP) :: QTRC_t(ka,ia,ja,qa)
529  real(RP) :: RHOE (ka,ia,ja)
530  real(RP) :: temp (ka,ia,ja)
531  real(RP) :: pres (ka,ia,ja)
532 
533  real(RP) :: vterm (ka,ia,ja,qa) ! terminal velocity of each tracer [m/s]
534  real(RP) :: FLX_rain (ka,ia,ja)
535  real(RP) :: FLX_snow (ka,ia,ja)
536  real(RP) :: wflux_rain(ka,ia,ja)
537  real(RP) :: wflux_snow(ka,ia,ja)
538  real(RP) :: qc_before_satadj(ka,ia,ja)
539  integer :: step
540  integer :: k, i, j
541  !---------------------------------------------------------------------------
542 
543  if( io_l ) write(io_fid_log,*) '*** Physics step: Cloud microphysics(tomita08)'
544 
545  if ( mp_donegative_fixer ) then
546  call mp_negative_fixer( dens(:,:,:), & ! [INOUT]
547  rhot(:,:,:), & ! [INOUT]
548  qtrc(:,:,:,:) ) ! [INOUT]
549  endif
550 
551  call thermodyn_rhoe( rhoe(:,:,:), & ! [OUT]
552  rhot(:,:,:), & ! [IN]
553  qtrc(:,:,:,:) ) ! [IN]
554 
555  !##### MP Main #####
556  call mp_tomita08( rhoe_t(:,:,:), & ! [OUT]
557  qtrc_t(:,:,:,:), & ! [OUT]
558  rhoe(:,:,:), & ! [INOUT]
559  qtrc(:,:,:,:), & ! [INOUT]
560  ccn(:,:,:), & ! [IN]
561  dens(:,:,:) ) ! [IN]
562 
563  flx_rain(:,:,:) = 0.0_rp
564  flx_snow(:,:,:) = 0.0_rp
565 
566  vterm(:,:,:,:) = 0.0_rp
567 
568  if ( mp_doprecipitation ) then
569 
570  do step = 1, mp_nstep_sedimentation
571 
572  call thermodyn_temp_pres_e( temp(:,:,:), & ! [OUT]
573  pres(:,:,:), & ! [OUT]
574  dens(:,:,:), & ! [IN]
575  rhoe(:,:,:), & ! [IN]
576  qtrc(:,:,:,:) ) ! [IN]
577 
578  call mp_tomita08_vterm( vterm(:,:,:,:), & ! [OUT]
579  dens(:,:,:), & ! [IN]
580  temp(:,:,:), & ! [IN]
581  qtrc(:,:,:,:) ) ! [IN]
582 
583  call mp_precipitation( wflux_rain(:,:,:), & ! [OUT]
584  wflux_snow(:,:,:), & ! [OUT]
585  dens(:,:,:), & ! [INOUT]
586  momz(:,:,:), & ! [INOUT]
587  momx(:,:,:), & ! [INOUT]
588  momy(:,:,:), & ! [INOUT]
589  rhoe(:,:,:), & ! [INOUT]
590  qtrc(:,:,:,:), & ! [INOUT]
591  vterm(:,:,:,:), & ! [IN]
592  temp(:,:,:), & ! [IN]
593  mp_dtsec_sedimentation ) ! [IN]
594 
595  do j = js, je
596  do i = is, ie
597  do k = ks-1, ke
598  flx_rain(k,i,j) = flx_rain(k,i,j) + wflux_rain(k,i,j) * mp_rnstep_sedimentation
599  flx_snow(k,i,j) = flx_snow(k,i,j) + wflux_snow(k,i,j) * mp_rnstep_sedimentation
600  enddo
601  enddo
602  enddo
603 
604  enddo
605 
606  else
607  vterm(:,:,:,:) = 0.0_rp
608  endif
609 
610  call hist_in( vterm(:,:,:,i_qr), 'Vterm_QR', 'terminal velocity of QR', 'm/s' )
611  call hist_in( vterm(:,:,:,i_qi), 'Vterm_QI', 'terminal velocity of QI', 'm/s' )
612  call hist_in( vterm(:,:,:,i_qs), 'Vterm_QS', 'terminal velocity of QS', 'm/s' )
613  call hist_in( vterm(:,:,:,i_qg), 'Vterm_QG', 'terminal velocity of QG', 'm/s' )
614 
615  do j = js, je
616  do i = is, ie
617  do k = ks, ke
618  qc_before_satadj(k,i,j) = qtrc(k,i,j,i_qc)
619  enddo
620  enddo
621  enddo
622 
623  call mp_saturation_adjustment( rhoe_t(:,:,:), & ! [INOUT]
624  qtrc_t(:,:,:,:), & ! [INOUT]
625  rhoe(:,:,:), & ! [INOUT]
626  qtrc(:,:,:,:), & ! [INOUT]
627  dens(:,:,:), & ! [IN]
628  only_liquid ) ! [IN]
629 
630  do j = js, je
631  do i = is, ie
632  do k = ks, ke
633  evaporate(k,i,j) = max( qc_before_satadj(k,i,j) - qtrc(k,i,j,i_qc), 0.0_rp ) / dt ! if negative, condensation
634  evaporate(k,i,j) = evaporate(k,i,j) * dens(k,i,j) &
635  / (4.0_rp/3.0_rp*pi*dwatr*re_qc**3) ! mass -> number (assuming constant particle radius as re_qc)
636  enddo
637  enddo
638  enddo
639 
640  !##### END MP Main #####
641 
642  call thermodyn_rhot( rhot(:,:,:), & ! [OUT]
643  rhoe(:,:,:), & ! [IN]
644  qtrc(:,:,:,:) ) ! [IN]
645 
646  if ( mp_donegative_fixer ) then
647  call mp_negative_fixer( dens(:,:,:), & ! [INOUT]
648  rhot(:,:,:), & ! [INOUT]
649  qtrc(:,:,:,:) ) ! [INOUT]
650  endif
651 
652  sflx_rain(:,:) = flx_rain(ks-1,:,:)
653  sflx_snow(:,:) = flx_snow(ks-1,:,:)
654 
655  return
656  end subroutine atmos_phy_mp_tomita08
657 
658  !-----------------------------------------------------------------------------
660  subroutine mp_tomita08( &
661  RHOE_t, &
662  QTRC_t, &
663  RHOE0, &
664  QTRC0, &
665  CCN, &
666  DENS0 )
667  use scale_const, only: &
668  pi => const_pi, &
669  eps => const_eps, &
670  rvap => const_rvap, &
671  cl => const_cl, &
672  tem00 => const_tem00, &
673  pre00 => const_pre00, &
674  dwatr => const_dwatr
675  use scale_time, only: &
676  dt => time_dtsec_atmos_phy_mp
677  use scale_history, only: &
678  hist_query, &
679  hist_put
680  use scale_atmos_thermodyn, only: &
681  thermodyn_temp_pres_e => atmos_thermodyn_temp_pres_e, &
682  atmos_thermodyn_templhv, &
683  atmos_thermodyn_templhs, &
684  atmos_thermodyn_templhf
685  use scale_atmos_saturation, only: &
686  saturation_dens2qsat_liq => atmos_saturation_dens2qsat_liq, &
687  saturation_dens2qsat_ice => atmos_saturation_dens2qsat_ice
688  implicit none
689 
690  real(RP), intent(out) :: RHOE_t(ka,ia,ja) ! tendency rhoe [J/m3/s]
691  real(RP), intent(out) :: QTRC_t(ka,ia,ja,qa) ! tendency tracer [kg/kg/s]
692  real(RP), intent(inout) :: RHOE0 (ka,ia,ja) ! density * internal energy [J/m3]
693  real(RP), intent(inout) :: QTRC0 (ka,ia,ja,qa) ! mass concentration [kg/kg]
694  real(RP), intent(in) :: CCN (ka,ia,ja) ! CCN number concentration [#/m3]
695  real(RP), intent(in) :: DENS0 (ka,ia,ja) ! density [kg/m3]
696 
697  ! working
698  real(RP) :: rdt
699 
700  real(RP) :: TEMP0(ka,ia,ja) ! temperature [K]
701  real(RP) :: PRES0(ka,ia,ja) ! pressure [Pa]
702  real(RP) :: QSATL(ka,ia,ja) ! saturated water vapor for liquid water [kg/kg]
703  real(RP) :: QSATI(ka,ia,ja) ! saturated water vapor for ice water [kg/kg]
704 
705  real(RP) :: dens
706  real(RP) :: rhoe
707  real(RP) :: temp
708  real(RP) :: pres
709  real(RP) :: q(qa)
710  real(RP) :: Sliq ! saturated ratio S for liquid water [0-1]
711  real(RP) :: Sice ! saturated ratio S for ice water [0-1]
712 
713  real(RP) :: Rdens
714  real(RP) :: rho_fact ! density factor
715  real(RP) :: temc ! T - T0 [K]
716 
717  real(RP) :: RLMDr, RLMDr_2, RLMDr_3
718  real(RP) :: RLMDs, RLMDs_2, RLMDs_3
719  real(RP) :: RLMDg, RLMDg_2, RLMDg_3
720  real(RP) :: RLMDr_1br, RLMDr_2br, RLMDr_3br
721  real(RP) :: RLMDs_1bs, RLMDs_2bs, RLMDs_3bs
722  real(RP) :: RLMDr_dr, RLMDr_3dr, RLMDr_5dr
723  real(RP) :: RLMDs_ds, RLMDs_3ds, RLMDs_5ds
724  real(RP) :: RLMDg_dg, RLMDg_3dg, RLMDg_5dg
725  real(RP) :: RLMDr_7
726  real(RP) :: RLMDr_6dr
727 
728  !---< Roh and Satoh (2014) >---
729  real(RP) :: tems, rhoqs, Xs2
730  real(RP) :: MOMs_0, MOMs_1, MOMs_2
731  real(RP) :: MOMs_0bs, MOMs_1bs, MOMs_2bs
732  real(RP) :: MOMs_2ds, MOMs_5ds_h, RMOMs_Vt
733  real(RP) :: coef_at(4), coef_bt(4)
734  real(RP) :: loga_, b_, nm
735 
736  real(RP) :: Vti, Vtr, Vts, Vtg
737  real(RP) :: Esi_mod, Egs_mod
738  real(RP) :: rhoqc
739  real(RP) :: Nc(ka,ia,ja)
740  real(RP) :: Pracw_orig, Pracw_kk
741  real(RP) :: Praut_berry, Praut_kk
742  real(RP) :: Dc
743  real(RP) :: betai, betas
744  real(RP) :: Da
745  real(RP) :: Kd
746  real(RP) :: Nu
747  real(RP) :: Glv, Giv, Gil
748  real(RP) :: ventr, vents, ventg
749  real(RP) :: Bergeron_sw
750  real(RP) :: a1 (ka,ia,ja)
751  real(RP) :: a2 (ka,ia,ja)
752  real(RP) :: ma2(ka,ia,ja)
753  real(RP) :: dt1
754  real(RP) :: Ni50
755 
756  ! for explicit ice generation
757  real(RP), parameter :: Di_max = 500.e-6_rp
758  real(RP), parameter :: Di_a = 11.9_rp
759  real(RP) :: sw, rhoqi, XNi, XMi, Di, Ni0, Qi0, dq, qtmp
760  real(RP) :: Pidep, Pigen
761 
762  real(RP) :: ice
763 
764  real(RP) :: net, fac, fac_sw, fac_warm, fac_ice
765  real(RP) :: w(w_nmax)
766 
767  real(RP) :: tend(i_qv:i_qg)
768 
769  real(RP) :: zerosw
770  real(RP) :: tmp
771  real(RP) :: LHVEx(ka,ia,ja)
772  real(RP) :: LHFEx(ka,ia,ja)
773  real(RP) :: LHSEx(ka,ia,ja)
774 
775  integer :: k, i, j, iq
776  !---------------------------------------------------------------------------
777 
778  call prof_rapstart('MP_tomita08', 3)
779 
780 ! RHOE_t(:,:,:) = 0.0_RP
781 ! QTRC_t(:,:,:,:) = 0.0_RP
782 
783  if ( mp_couple_aerosol ) then
784  do j = js, je
785  do i = is, ie
786  do k = ks, ke
787  nc(k,i,j) = max( ccn(k,i,j)*1.e-6_rp, nc_def(i,j) ) ! [#/m3]->[#/cc]
788  enddo
789  enddo
790  enddo
791  else
792  do j = js, je
793  do i = is, ie
794  do k = ks, ke
795  nc(k,i,j) = nc_def(i,j)
796  enddo
797  enddo
798  enddo
799  endif
800 
801  rdt = 1.0_rp / dt
802 
803  call thermodyn_temp_pres_e( temp0(:,:,:), & ! [OUT]
804  pres0(:,:,:), & ! [OUT]
805  dens0(:,:,:), & ! [IN]
806  rhoe0(:,:,:), & ! [IN]
807  qtrc0(:,:,:,:) ) ! [IN]
808 
809  call saturation_dens2qsat_liq( qsatl(:,:,:), & ! [OUT]
810  temp0(:,:,:), & ! [IN]
811  dens0(:,:,:) ) ! [IN]
812 
813  call saturation_dens2qsat_ice( qsati(:,:,:), & ! [OUT]
814  temp0(:,:,:), & ! [IN]
815  dens0(:,:,:) ) ! [IN]
816 
817  call mp_tomita08_bergeronparam( temp0(:,:,:), & ! [IN]
818  a1(:,:,:), & ! [OUT]
819  a2(:,:,:), & ! [OUT]
820  ma2(:,:,:) ) ! [OUT]
821 
822  call atmos_thermodyn_templhv( lhvex, temp0 )
823  call atmos_thermodyn_templhf( lhfex, temp0 )
824  call atmos_thermodyn_templhs( lhsex, temp0 )
825 
826 !$omp parallel do &
827 !$omp private(tend, coef_bt, coef_at, q, w) &
828 !$omp private(dens, rhoe, temp, pres) &
829 !$omp private(Sliq, Sice, Rdens, rho_fact, temc) &
830 !$omp private(Bergeron_sw, zerosw) &
831 !$omp private(RLMDr, RLMDr_dr, RLMDr_2, RLMDr_3, RLMDr_7, RLMDr_1br, RLMDr_2br, RLMDr_3br, RLMDr_3dr, RLMDr_5dr, RLMDr_6dr) &
832 !$omp private(RLMDs, RLMDs_ds, RLMDs_2, RLMDs_3, RLMDs_1bs, RLMDs_2bs, RLMDs_3bs, RLMDs_3ds, RLMDs_5ds) &
833 !$omp private(MOMs_0, MOMs_1, MOMs_2, MOMs_0bs, MOMs_1bs, MOMs_2bs, MOMs_2ds, MOMs_5ds_h, RMOMs_Vt) &
834 !$omp private(rhoqc, rhoqs, Xs2, tems, loga_, b_, nm) &
835 !$omp private(RLMDg, RLMDg_dg, RLMDg_2, RLMDg_3, RLMDg_3dg, RLMDg_5dg) &
836 !$omp private(Vti, Vtr, Vts, Vtg, Esi_mod, Egs_mod, Dc, Praut_berry, Praut_kk, betai, betas, Da, Kd, NU, Glv, Giv, Gil, ventr, vents, ventg) &
837 !$omp private(tmp, dt1, Ni50, ice, net, fac_sw, fac, fac_warm, fac_ice) &
838 !$omp collapse(3)
839 !OCL TEMP_PRIVATE(tend, coef_bt, coef_at, q, w)
840  do j = js, je
841  do i = is, ie
842  do k = ks, ke
843 
844  ! store to work
845  dens = dens0(k,i,j)
846  rhoe = rhoe0(k,i,j)
847  temp = temp0(k,i,j)
848  pres = pres0(k,i,j)
849  do iq = i_qv, i_qg
850  q(iq) = qtrc0(k,i,j,iq)
851  enddo
852 
853  ! saturation ratio S
854  sliq = q(i_qv) / max( qsatl(k,i,j), eps )
855  sice = q(i_qv) / max( qsati(k,i,j), eps )
856 
857  rdens = 1.0_rp / dens
858  rho_fact = sqrt( dens00 * rdens )
859  temc = temp - tem00
860 
861  w(i_delta1) = ( 0.5_rp + sign(0.5_rp, q(i_qr) - 1.e-4_rp ) )
862 
863  w(i_delta2) = ( 0.5_rp + sign(0.5_rp, 1.e-4_rp - q(i_qr) ) ) &
864  * ( 0.5_rp + sign(0.5_rp, 1.e-4_rp - q(i_qs) ) )
865 
866  w(i_delta3) = 0.5_rp + sign(0.5_rp, sice - 1.0_rp )
867 
868  w(i_dqv_dt) = q(i_qv) * rdt
869  w(i_dqc_dt) = q(i_qc) * rdt
870  w(i_dqr_dt) = q(i_qr) * rdt
871  w(i_dqi_dt) = q(i_qi) * rdt
872  w(i_dqs_dt) = q(i_qs) * rdt
873  w(i_dqg_dt) = q(i_qg) * rdt
874 
875  bergeron_sw = ( 0.5_rp + sign(0.5_rp, temc + 30.0_rp ) ) &
876  * ( 0.5_rp + sign(0.5_rp, 0.0_rp - temc ) ) &
877  * ( 1.0_rp - sw_useicegen )
878 
879  ! slope parameter lambda (Rain)
880  zerosw = 0.5_rp - sign(0.5_rp, q(i_qr) - 1.e-12_rp )
881  rlmdr = sqrt(sqrt( dens * q(i_qr) / ( ar * n0r * gam_1br ) + zerosw )) * ( 1.0_rp - zerosw )
882 
883  rlmdr_dr = sqrt( rlmdr ) ! **Dr
884  rlmdr_2 = rlmdr**2
885  rlmdr_3 = rlmdr**3
886  rlmdr_7 = rlmdr**7
887  rlmdr_1br = rlmdr**4 ! (1+Br)
888  rlmdr_2br = rlmdr**5 ! (2+Br)
889  rlmdr_3br = rlmdr**6 ! (3+Br)
890  rlmdr_3dr = rlmdr**3 * rlmdr_dr
891  rlmdr_5dr = rlmdr**5 * rlmdr_dr
892  rlmdr_6dr = rlmdr**6 * rlmdr_dr
893 
894  ! slope parameter lambda (Snow)
895  zerosw = 0.5_rp - sign(0.5_rp, q(i_qs) - 1.e-12_rp )
896  rlmds = sqrt(sqrt( dens * q(i_qs) / ( as * n0s * gam_1bs ) + zerosw )) * ( 1.0_rp - zerosw )
897 
898  rlmds_ds = sqrt( sqrt(rlmds) ) ! **Ds
899  rlmds_2 = rlmds**2
900  rlmds_3 = rlmds**3
901  rlmds_1bs = rlmds**4 ! (1+Bs)
902  rlmds_2bs = rlmds**5 ! (2+Bs)
903  rlmds_3bs = rlmds**6 ! (3+Bs)
904  rlmds_3ds = rlmds**3 * rlmds_ds
905  rlmds_5ds = rlmds**5 * rlmds_ds
906 
907  moms_0 = n0s * gam * rlmds ! Ns * 0th moment
908  moms_1 = n0s * gam_2 * rlmds_2 ! Ns * 1st moment
909  moms_2 = n0s * gam_3 * rlmds_3 ! Ns * 2nd moment
910  moms_0bs = n0s * gam_1bs * rlmds_1bs ! Ns * 0+bs
911  moms_1bs = n0s * gam_2bs * rlmds_2bs ! Ns * 1+bs
912  moms_2bs = n0s * gam_3bs * rlmds_3bs ! Ns * 2+bs
913  moms_2ds = n0s * gam_3ds * rlmds_3ds ! Ns * 2+ds
914  moms_5ds_h = n0s * gam_5ds_h * sqrt(rlmds_5ds) ! Ns * (5+ds)/2
915  rmoms_vt = gam_1bsds / gam_1bs * rlmds_ds
916 
917  !---< modification by Roh and Satoh (2014) >---
918  ! bimodal size distribution of snow
919  rhoqs = dens * q(i_qs)
920  zerosw = 0.5_rp - sign(0.5_rp, rhoqs - 1.e-12_rp )
921  xs2 = rhoqs / as * ( 1.0_rp - zerosw )
922 
923  tems = min( -0.1_rp, temc )
924  coef_at(1) = coef_a( 1) + tems * ( coef_a( 2) + tems * ( coef_a( 5) + tems * coef_a( 9) ) )
925  coef_at(2) = coef_a( 3) + tems * ( coef_a( 4) + tems * coef_a( 7) )
926  coef_at(3) = coef_a( 6) + tems * coef_a( 8)
927  coef_at(4) = coef_a(10)
928  coef_bt(1) = coef_b( 1) + tems * ( coef_b( 2) + tems * ( coef_b( 5) + tems * coef_b( 9) ) )
929  coef_bt(2) = coef_b( 3) + tems * ( coef_b( 4) + tems * coef_b( 7) )
930  coef_bt(3) = coef_b( 6) + tems * coef_b( 8)
931  coef_bt(4) = coef_b(10)
932  ! 0th moment
933  loga_ = coef_at(1)
934  b_ = coef_bt(1)
935  moms_0 = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
936  + ( 1.0_rp-sw_roh2014 ) * moms_0
937  ! 1st moment
938  nm = 1.0_rp
939  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
940  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
941  moms_1 = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
942  + ( 1.0_rp-sw_roh2014 ) * moms_1
943  ! 2nd moment
944  moms_2 = ( sw_roh2014 ) * xs2 &
945  + ( 1.0_rp-sw_roh2014 ) * moms_2
946  ! 0 + Bs(=2) moment
947  nm = 2.0_rp
948  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
949  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
950  moms_0bs = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
951  + ( 1.0_rp-sw_roh2014 ) * moms_0bs
952  ! 1 + Bs(=2) moment
953  nm = 3.0_rp
954  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
955  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
956  moms_1bs = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
957  + ( 1.0_rp-sw_roh2014 ) * moms_1bs
958  ! 2 + Bs(=2) moment
959  nm = 4.0_rp
960  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
961  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
962  moms_2bs = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
963  + ( 1.0_rp-sw_roh2014 ) * moms_2bs
964  ! 2 + Ds(=0.25) moment
965  nm = 2.25_rp
966  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
967  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
968  moms_2ds = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
969  + ( 1.0_rp-sw_roh2014 ) * moms_2ds
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 = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
975  + ( 1.0_rp-sw_roh2014 ) * moms_5ds_h
976  ! Bs(=2) + Ds(=0.25) moment
977  nm = 2.25_rp
978  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
979  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
980  rmoms_vt = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ / ( moms_0bs + zerosw ) &
981  + ( 1.0_rp-sw_roh2014 ) * rmoms_vt
982 
983  ! slope parameter lambda (Graupel)
984  zerosw = 0.5_rp - sign(0.5_rp, q(i_qg) - 1.e-12_rp )
985  rlmdg = sqrt(sqrt( dens * q(i_qg) / ( ag * n0g * gam_1bg ) + zerosw )) * ( 1.0_rp - zerosw )
986 
987  rlmdg_dg = sqrt( rlmdg ) ! **Dg
988  rlmdg_2 = rlmdg**2
989  rlmdg_3 = rlmdg**3
990  rlmdg_3dg = rlmdg**3 * rlmdg_dg
991  rlmdg_5dg = rlmdg**5 * rlmdg_dg
992 
993  w(i_rlmdr) = rlmdr
994  w(i_rlmds) = rlmds
995  w(i_rlmdg) = rlmdg
996 
997  !---< terminal velocity >
998  zerosw = 0.5_rp + sign(0.5_rp, q(i_qi) - 1.e-8_rp )
999  vti = -3.29_rp * ( dens * q(i_qi) * zerosw )**0.16_rp
1000  vtr = -cr * rho_fact * gam_1brdr / gam_1br * rlmdr_dr
1001  vts = -cs * rho_fact * rmoms_vt
1002  vtg = -cg * rho_fact * gam_1bgdg / gam_1bg * rlmdg_dg
1003 
1004  !---< Accretion >---
1005  esi_mod = min( esi, esi * exp( gamma_sacr * temc ) )
1006  egs_mod = min( egs, egs * exp( gamma_gacs * temc ) )
1007 
1008  ! [Pracw] accretion rate of cloud water by rain
1009  pracw_orig = q(i_qc) * 0.25_rp * pi * erw * n0r * cr * gam_3dr * rlmdr_3dr * rho_fact
1010 
1011  pracw_kk = 67.0_rp * ( q(i_qc) * q(i_qr) )**1.15_rp ! eq.(33) in KK(2000)
1012 
1013  ! switch orig / k-k scheme
1014  w(i_pracw) = ( 1.0_rp - sw_kk2000 ) * pracw_orig &
1015  + ( sw_kk2000 ) * pracw_kk
1016 
1017  ! [Psacw] accretion rate of cloud water by snow
1018  w(i_psacw) = q(i_qc) * 0.25_rp * pi * esw * cs * moms_2ds * rho_fact
1019 
1020  ! [Pgacw] accretion rate of cloud water by graupel
1021  w(i_pgacw) = q(i_qc) * 0.25_rp * pi * egw * n0g * cg * gam_3dg * rlmdg_3dg * rho_fact
1022 
1023  ! [Praci] accretion rate of cloud ice by rain
1024  w(i_praci) = q(i_qi) * 0.25_rp * pi * eri * n0r * cr * gam_3dr * rlmdr_3dr * rho_fact
1025 
1026  ! [Psaci] accretion rate of cloud ice by snow
1027  w(i_psaci) = q(i_qi) * 0.25_rp * pi * esi_mod * cs * moms_2ds * rho_fact
1028 
1029  ! [Pgaci] accretion rate of cloud ice by grupel
1030  w(i_pgaci) = q(i_qi) * 0.25_rp * pi * egi * n0g * cg * gam_3dg * rlmdg_3dg * rho_fact * ( 1.0_rp-sw_roh2014 )
1031 
1032  ! [Piacr] accretion rate of rain by cloud ice
1033  w(i_piacr) = q(i_qi) * ar / 4.19e-13_rp * 0.25_rp * pi * eri * n0r * cr * gam_6dr * rlmdr_6dr * rho_fact
1034 
1035  ! [Psacr] accretion rate of rain by snow
1036  w(i_psacr) = ar * 0.25_rp * pi * rdens * esr * n0r * abs(vtr-vts) &
1037  * ( gam_1br * rlmdr_1br * moms_2 &
1038  + 2.0_rp * gam_2br * rlmdr_2br * moms_1 &
1039  + gam_3br * rlmdr_3br * moms_0 )
1040 
1041  ! [Pgacr] accretion rate of rain by graupel
1042  w(i_pgacr) = ar * 0.25_rp * pi * rdens * egr * n0g * n0r * abs(vtg-vtr) &
1043  * ( gam_1br * rlmdr_1br * gam_3 * rlmdg_3 &
1044  + 2.0_rp * gam_2br * rlmdr_2br * gam_2 * rlmdg_2 &
1045  + gam_3br * rlmdr_3br * gam * rlmdg )
1046 
1047  ! [Pracs] accretion rate of snow by rain
1048  w(i_pracs) = as * 0.25_rp * pi * rdens * esr * n0r * abs(vtr-vts) &
1049  * ( moms_0bs * gam_3 * rlmdr_3 &
1050  + 2.0_rp * moms_1bs * gam_2 * rlmdr_2 &
1051  + moms_2bs * gam * rlmdr )
1052 
1053  ! [Pgacs] accretion rate of snow by graupel
1054  w(i_pgacs) = as * 0.25_rp * pi * rdens * egs_mod * n0g * abs(vtg-vts) * ( 1.0_rp-sw_roh2014 ) &
1055  * ( moms_0bs * gam_3 * rlmdg_3 &
1056  + 2.0_rp * moms_1bs * gam_2 * rlmdg_2 &
1057  + moms_2bs * gam * rlmdg )
1058 
1059  !---< Auto-conversion >---
1060  ! [Praut] auto-conversion rate from cloud water to rain
1061  rhoqc = dens * q(i_qc) * 1000.0_rp ! [g/m3]
1062  dc = 0.146_rp - 5.964e-2_rp * log( nc(k,i,j) / 2000.0_rp )
1063  praut_berry = rdens * 1.67e-5_rp * rhoqc * rhoqc / ( 5.0_rp + 3.66e-2_rp * nc(k,i,j) / ( dc * rhoqc + eps ) )
1064 
1065  praut_kk = 1350.0_rp * q(i_qc)**2.47_rp * nc(k,i,j)**(-1.79_rp) ! eq.(29) in KK(2000)
1066 
1067  ! switch berry / k-k scheme
1068  w(i_praut) = ( 1.0_rp - sw_kk2000 ) * praut_berry &
1069  + ( sw_kk2000 ) * praut_kk
1070 
1071  ! [Psaut] auto-conversion rate from cloud ice to snow
1072  betai = min( beta_saut, beta_saut * exp( gamma_saut * temc ) )
1073  w(i_psaut) = max( betai*(q(i_qi)-qicrt_saut), 0.0_rp )
1074 
1075  ! [Pgaut] auto-conversion rate from snow to graupel
1076  betas = min( beta_gaut, beta_gaut * exp( gamma_gaut * temc ) )
1077  w(i_pgaut) = max( betas*(q(i_qs)-qscrt_gaut), 0.0_rp )
1078 
1079  !---< Evaporation, Sublimation >---
1080  da = ( da0 + dda_dt * temc )
1081  kd = ( dw0 + ddw_dt * temc ) * pre00 / pres
1082  nu = ( mu0 + dmu_dt * temc ) * rdens
1083 
1084  glv = 1.0_rp / ( lhvex(k,i,j)/(da*temp) * ( lhvex(k,i,j)/(rvap*temp) - 1.0_rp ) + 1.0_rp/(kd*dens*qsatl(k,i,j)) )
1085  giv = 1.0_rp / ( lhsex(k,i,j)/(da*temp) * ( lhsex(k,i,j)/(rvap*temp) - 1.0_rp ) + 1.0_rp/(kd*dens*qsati(k,i,j)) )
1086  gil = 1.0_rp / lhfex(k,i,j) * (da*temc)
1087 
1088  ! [Prevp] evaporation rate of rain
1089  ventr = f1r * gam_2 * rlmdr_2 + f2r * sqrt( cr * rho_fact / nu * rlmdr_5dr ) * gam_5dr_h
1090 
1091  w(i_prevp) = 2.0_rp * pi * rdens * n0r * ( 1.0_rp-min(sliq,1.0_rp) ) * glv * ventr
1092 
1093  ! [Psdep,Pssub] deposition/sublimation rate for snow
1094  vents = f1s * moms_1 + f2s * sqrt( cs * rho_fact / nu ) * moms_5ds_h
1095 
1096  tmp = 2.0_rp * pi * rdens * ( sice-1.0_rp ) * giv * vents
1097 
1098  w(i_psdep) = ( w(i_delta3) ) * tmp ! Sice < 1
1099  w(i_pssub) = ( w(i_delta3)-1.0_rp ) * tmp ! Sice > 1
1100 
1101  ! [Psmlt] melting rate of snow
1102  w(i_psmlt) = 2.0_rp * pi * rdens * gil * vents &
1103  + cl * temc / lhfex(k,i,j) * ( w(i_psacw) + w(i_psacr) )
1104 
1105  ! [Pgdep/pgsub] deposition/sublimation rate for graupel
1106  ventg = f1g * gam_2 * rlmdg_2 + f2g * sqrt( cg * rho_fact / nu * rlmdg_5dg ) * gam_5dg_h
1107 
1108  tmp = 2.0_rp * pi * rdens * n0g * ( sice-1.0_rp ) * giv * ventg
1109 
1110  w(i_pgdep) = ( w(i_delta3) ) * tmp ! Sice < 1
1111  w(i_pgsub) = ( w(i_delta3)-1.0_rp ) * tmp ! Sice > 1
1112 
1113  ! [Pgmlt] melting rate of graupel
1114  w(i_pgmlt) = 2.0_rp * pi * rdens * n0g * gil * ventg &
1115  + cl * temc / lhfex(k,i,j) * ( w(i_pgacw) + w(i_pgacr) )
1116 
1117  ! [Pgfrz] freezing rate of graupel
1118  w(i_pgfrz) = 2.0_rp * pi * rdens * n0r * 60.0_rp * b_frz * ar * ( exp(-a_frz*temc) - 1.0_rp ) * rlmdr_7
1119 
1120  ! [Psfw,Psfi] ( Bergeron process ) growth rate of snow by Bergeron process from cloud water/ice
1121  dt1 = ( mi50**ma2(k,i,j) - mi40**ma2(k,i,j) ) / ( a1(k,i,j) * ma2(k,i,j) )
1122  ni50 = q(i_qi) * dt / ( mi50 * dt1 )
1123 
1124  w(i_psfw) = bergeron_sw * ni50 * ( a1(k,i,j) * mi50**a2(k,i,j) &
1125  + pi * eiw * dens * q(i_qc) * ri50*ri50 * vti50 )
1126  w(i_psfi) = bergeron_sw * q(i_qi) / dt1
1127 
1128  w(i_psdep) = min( w(i_psdep), w(i_dqv_dt) )
1129  w(i_pgdep) = min( w(i_pgdep), w(i_dqv_dt) )
1130 
1131  w(i_praut) = min( w(i_praut), w(i_dqc_dt) )
1132  w(i_pracw) = min( w(i_pracw), w(i_dqc_dt) )
1133  w(i_psacw) = min( w(i_psacw), w(i_dqc_dt) )
1134  w(i_pgacw) = min( w(i_pgacw), w(i_dqc_dt) )
1135  w(i_psfw ) = min( w(i_psfw ), w(i_dqc_dt) )
1136 
1137  w(i_prevp) = min( w(i_prevp), w(i_dqr_dt) )
1138  w(i_piacr) = min( w(i_piacr), w(i_dqr_dt) )
1139  w(i_psacr) = min( w(i_psacr), w(i_dqr_dt) )
1140  w(i_pgacr) = min( w(i_pgacr), w(i_dqr_dt) )
1141  w(i_pgfrz) = min( w(i_pgfrz), w(i_dqr_dt) )
1142 
1143  w(i_psaut) = min( w(i_psaut), w(i_dqi_dt) )
1144  w(i_praci) = min( w(i_praci), w(i_dqi_dt) )
1145  w(i_psaci) = min( w(i_psaci), w(i_dqi_dt) )
1146  w(i_pgaci) = min( w(i_pgaci), w(i_dqi_dt) )
1147  w(i_psfi ) = min( w(i_psfi ), w(i_dqi_dt) )
1148 
1149  w(i_pgaut) = min( w(i_pgaut), w(i_dqs_dt) )
1150  w(i_pracs) = min( w(i_pracs), w(i_dqs_dt) )
1151  w(i_pgacs) = min( w(i_pgacs), w(i_dqs_dt) )
1152  w(i_psmlt) = max( w(i_psmlt), 0.0_rp )
1153  w(i_psmlt) = min( w(i_psmlt), w(i_dqs_dt) )
1154  w(i_pssub) = min( w(i_pssub), w(i_dqs_dt) )
1155 
1156  w(i_pgmlt) = max( w(i_pgmlt), 0.0_rp )
1157  w(i_pgmlt) = min( w(i_pgmlt), w(i_dqg_dt) )
1158  w(i_pgsub) = min( w(i_pgsub), w(i_dqg_dt) )
1159 
1160  w(i_piacr_s) = ( 1.0_rp - w(i_delta1) ) * w(i_piacr)
1161  w(i_piacr_g) = ( w(i_delta1) ) * w(i_piacr)
1162  w(i_praci_s) = ( 1.0_rp - w(i_delta1) ) * w(i_praci)
1163  w(i_praci_g) = ( w(i_delta1) ) * w(i_praci)
1164  w(i_psacr_s) = ( w(i_delta2) ) * w(i_psacr)
1165  w(i_psacr_g) = ( 1.0_rp - w(i_delta2) ) * w(i_psacr)
1166  w(i_pracs ) = ( 1.0_rp - w(i_delta2) ) * w(i_pracs)
1167 
1168  ice = 0.5_rp - sign( 0.5_rp, temp0(k,i,j)-tem00 ) ! 0: warm, 1: ice
1169 
1170  ! [QC]
1171  net = &
1172  - w(i_praut) & ! [loss] c->r
1173  - w(i_pracw) & ! [loss] c->r
1174  - w(i_psacw) & ! [loss] c->s
1175  - w(i_pgacw) & ! [loss] c->g
1176  - w(i_psfw ) * ice ! [loss] c->s
1177 
1178  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1179  fac = ( fac_sw ) &
1180  + ( 1.0_rp - fac_sw ) * min( -w(i_dqc_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1181 
1182  w(i_praut) = w(i_praut) * fac
1183  w(i_pracw) = w(i_pracw) * fac
1184  w(i_psacw) = w(i_psacw) * fac
1185  w(i_pgacw) = w(i_pgacw) * fac
1186  w(i_psfw ) = w(i_psfw ) * fac
1187 
1188  ! [QI]
1189  net = &
1190  - w(i_psaut ) & ! [loss] i->s
1191  - w(i_praci_s) & ! [loss] i->s
1192  - w(i_psaci ) & ! [loss] i->s
1193  - w(i_psfi ) & ! [loss] i->s
1194  - w(i_praci_g) & ! [loss] i->g
1195  - w(i_pgaci ) ! [loss] i->g
1196 
1197  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1198  fac = ( fac_sw ) &
1199  + ( 1.0_rp - fac_sw ) * min( -w(i_dqi_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1200 
1201  w(i_praci_s) = w(i_praci_s) * fac
1202  w(i_psaut ) = w(i_psaut ) * fac
1203  w(i_psaci ) = w(i_psaci ) * fac
1204  w(i_psfi ) = w(i_psfi ) * fac
1205  w(i_praci_g) = w(i_praci_g) * fac
1206  w(i_pgaci ) = w(i_pgaci ) * fac
1207 
1208  ! [QR]
1209  net = w(i_praut ) & ! [prod] c->r
1210  + w(i_pracw ) & ! [prod] c->r
1211  - w(i_prevp ) & ! [loss] r->v
1212  + ( &
1213  + w(i_psacw ) & ! [prod] c->s
1214  + w(i_pgacw ) & ! [prod] c->g
1215  + w(i_psmlt ) & ! [prod] s->r
1216  + w(i_pgmlt ) & ! [prod] g->r
1217  ) * ( 1.0_rp - ice ) &
1218  + ( &
1219  - w(i_piacr_s) & ! [loss] r->s
1220  - w(i_psacr_s) & ! [loss] r->s
1221  - w(i_piacr_g) & ! [loss] r->g
1222  - w(i_psacr_g) & ! [loss] r->g
1223  - w(i_pgacr ) & ! [loss] r->g
1224  - w(i_pgfrz ) & ! [loss] r->g
1225  ) * ice
1226 
1227  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1228  fac = ( fac_sw ) &
1229  + ( 1.0_rp - fac_sw ) * min( -w(i_dqr_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1230 
1231  fac_warm = fac * (1.0_rp-ice) + 1.0_rp * ice
1232  fac_ice = fac * ice + 1.0_rp * (1.0_rp-ice)
1233 
1234  w(i_praut) = w(i_praut) * fac
1235  w(i_pracw) = w(i_pracw) * fac
1236  w(i_prevp) = w(i_prevp) * fac
1237 
1238  w(i_psacw) = w(i_psacw) * fac_warm
1239  w(i_pgacw) = w(i_pgacw) * fac_warm
1240  w(i_psmlt) = w(i_psmlt) * fac_warm
1241  w(i_pgmlt) = w(i_pgmlt) * fac_warm
1242 
1243  w(i_piacr_s) = w(i_piacr_s) * fac_ice
1244  w(i_psacr_s) = w(i_psacr_s) * fac_ice
1245  w(i_piacr_g) = w(i_piacr_g) * fac_ice
1246  w(i_psacr_g) = w(i_psacr_g) * fac_ice
1247  w(i_pgacr ) = w(i_pgacr ) * fac_ice
1248  w(i_pgfrz ) = w(i_pgfrz ) * fac_ice
1249 
1250  !### in the case of Sice-1 > 0, limiter is applied to qv before qs,qg
1251  ! [QV]
1252  net = w(i_prevp) & ! [prod] r->v
1253  + w(i_pssub) & ! [prod] s->v
1254  + w(i_pgsub) & ! [prod] g->v
1255  - w(i_psdep) & ! [loss] v->s
1256  - w(i_pgdep) ! [loss] v->g
1257 
1258  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1259  fac = ( fac_sw ) &
1260  + ( 1.0_rp - fac_sw ) * min( -w(i_dqv_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1261 
1262  fac = ( w(i_delta3) ) * fac &
1263  + ( 1.0_rp - w(i_delta3) )
1264 
1265  fac = fac * ice + 1.0_rp * (1.0_rp - ice)
1266 
1267  w(i_prevp) = w(i_prevp) * fac
1268  w(i_pssub) = w(i_pssub) * fac
1269  w(i_pgsub) = w(i_pgsub) * fac
1270  w(i_psdep) = w(i_psdep) * fac
1271  w(i_pgdep) = w(i_pgdep) * fac
1272 
1273  ! [QS]
1274  net = &
1275  - w(i_pgacs ) & ! [loss] s->g
1276  + ( &
1277  - w(i_psmlt ) & ! [loss] s->r
1278  ) * ( 1.0_rp - ice ) &
1279  + ( &
1280  + w(i_psdep ) & ! [prod] v->s
1281  + w(i_psacw ) & ! [prod] c->s
1282  + w(i_psfw ) & ! [prod] c->s
1283  + w(i_piacr_s) & ! [prod] r->s
1284  + w(i_psacr_s) & ! [prod] r->s
1285  + w(i_psaut ) & ! [prod] i->s
1286  + w(i_praci_s) & ! [prod] i->s
1287  + w(i_psaci ) & ! [prod] i->s
1288  + w(i_psfi ) & ! [prod] i->s
1289  - w(i_pssub ) & ! [loss] s->v
1290  - w(i_pgaut ) & ! [loss] s->g
1291  - w(i_pracs ) & ! [loss] s->g
1292  )
1293 
1294  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1295  fac = ( fac_sw ) &
1296  + ( 1.0_rp - fac_sw ) * min( -w(i_dqs_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1297 
1298  fac_warm = fac * (1.0_rp - ice) + 1.0_rp * ice
1299  fac_ice = fac * ice + 1.0_rp * (1.0_rp-ice)
1300 
1301  w(i_pgacs) = w(i_pgacs) * fac
1302 
1303  w(i_psmlt) = w(i_psmlt) * fac_warm
1304 
1305  w(i_psdep ) = w(i_psdep ) * fac_ice
1306  w(i_psacw ) = w(i_psacw ) * fac_ice
1307  w(i_psfw ) = w(i_psfw ) * fac_ice
1308  w(i_piacr_s) = w(i_piacr_s) * fac_ice
1309  w(i_psacr_s) = w(i_psacr_s) * fac_ice
1310  w(i_psaut ) = w(i_psaut ) * fac_ice
1311  w(i_praci_s) = w(i_praci_s) * fac_ice
1312  w(i_psaci ) = w(i_psaci ) * fac_ice
1313  w(i_psfi ) = w(i_psfi ) * fac_ice
1314  w(i_pssub ) = w(i_pssub ) * fac_ice
1315  w(i_pracs ) = w(i_pracs ) * fac_ice
1316  w(i_pgaut ) = w(i_pgaut ) * fac_ice
1317 
1318  ! [QG]
1319  net = w(i_pgacs ) & ! [prod] s->g
1320  + ( &
1321  - w(i_pgmlt ) & ! [loss] g->r
1322  ) * ( 1.0_rp - ice ) &
1323  + ( &
1324  + w(i_pgdep ) & ! [prod] v->g
1325  + w(i_pgacw ) & ! [prod] c->g
1326  + w(i_piacr_g) & ! [prod] r->g
1327  + w(i_psacr_g) & ! [prod] r->g
1328  + w(i_pgacr ) & ! [prod] r->g
1329  + w(i_pgfrz ) & ! [prod] r->g
1330  + w(i_praci_g) & ! [prod] i->g
1331  + w(i_pgaci ) & ! [prod] i->g
1332  + w(i_pgaut ) & ! [prod] s->g
1333  + w(i_pracs ) & ! [prod] s->g
1334  - w(i_pgsub ) & ! [loss] g->v
1335  ) * ice
1336 
1337  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1338  fac = ( fac_sw ) &
1339  + ( 1.0_rp - fac_sw ) * min( -w(i_dqg_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1340 
1341  fac_warm = fac * (1.0_rp - ice) + 1.0_rp * ice
1342  fac_ice = fac * ice + 1.0_rp * (1.0_rp-ice)
1343 
1344  w(i_pgacs) = w(i_pgacs) * fac
1345 
1346  w(i_pgmlt) = w(i_pgmlt) * fac_warm
1347 
1348  w(i_pgdep ) = w(i_pgdep ) * fac_ice
1349  w(i_pgacw ) = w(i_pgacw ) * fac_ice
1350  w(i_piacr_g) = w(i_piacr_g) * fac_ice
1351  w(i_psacr_g) = w(i_psacr_g) * fac_ice
1352  w(i_pgacr ) = w(i_pgacr ) * fac_ice
1353  w(i_pgfrz ) = w(i_pgfrz ) * fac_ice
1354  w(i_praci_g) = w(i_praci_g) * fac_ice
1355  w(i_pgaci ) = w(i_pgaci ) * fac_ice
1356  w(i_pgaut ) = w(i_pgaut ) * fac_ice
1357  w(i_pracs ) = w(i_pracs ) * fac_ice
1358  w(i_pgsub ) = w(i_pgsub ) * fac_ice
1359 
1360 
1361 
1362  ! re-calc net production & loss
1363  tend(i_qc) = &
1364  - w(i_praut ) & ! [loss] c->r
1365  - w(i_pracw ) & ! [loss] c->r
1366  - w(i_psacw ) & ! [loss] c->s
1367  - w(i_pgacw ) & ! [loss] c->g
1368  - w(i_psfw ) * ice
1369 
1370  tend(i_qr) = w(i_praut ) & ! [prod] c->r
1371  + w(i_pracw ) & ! [prod] c->r
1372  - w(i_prevp ) & ! [loss] r->v
1373  + ( &
1374  + w(i_psacw ) & ! [prod] c->s=r
1375  + w(i_pgacw ) & ! [prod] c->g=r
1376  + w(i_psmlt ) & ! [prod] s->r
1377  + w(i_pgmlt ) & ! [prod] g->r
1378  ) * ( 1.0_rp - ice ) &
1379  + ( &
1380  - w(i_piacr_s) & ! [loss] r->s
1381  - w(i_psacr_s) & ! [loss] r->s
1382  - w(i_piacr_g) & ! [loss] r->g
1383  - w(i_psacr_g) & ! [loss] r->g
1384  - w(i_pgacr ) & ! [loss] r->g
1385  - w(i_pgfrz ) & ! [loss] r->g
1386  ) * ice
1387 
1388  tend(i_qi) = ( &
1389  - w(i_psaut ) & ! [loss] i->s
1390  - w(i_praci_s) & ! [loss] i->s
1391  - w(i_psaci ) & ! [loss] i->s
1392  - w(i_psfi ) & ! [loss] i->s
1393  - w(i_praci_g) & ! [loss] i->g
1394  - w(i_pgaci ) & ! [loss] i->g
1395  ) * ice
1396 
1397  tend(i_qs) = &
1398  - w(i_pgacs ) & ! [loss] s->g
1399  + ( &
1400  - w(i_psmlt ) & ! [loss] s->r
1401  ) * ( 1.0_rp - ice ) &
1402  + ( &
1403  + w(i_psdep ) & ! [prod] v->s
1404  + w(i_psacw ) & ! [prod] c->s
1405  + w(i_psfw ) & ! [prod] c->s
1406  + w(i_piacr_s) & ! [prod] r->s
1407  + w(i_psacr_s) & ! [prod] r->s
1408  + w(i_psaut ) & ! [prod] i->s
1409  + w(i_psaci ) & ! [prod] i->s
1410  + w(i_praci_s) & ! [prod] i->s
1411  + w(i_psfi ) & ! [prod] i->s
1412  - w(i_pssub ) & ! [loss] s->v
1413  - w(i_pgaut ) & ! [loss] s->g
1414  - w(i_pracs ) & ! [loss] s->g
1415  ) * ice
1416 
1417  tend(i_qg) = w(i_pgacs ) & ! [prod] s->g
1418  + ( &
1419  - w(i_pgmlt ) & ! [loss] g->r
1420  ) * ( 1.0_rp - ice ) &
1421  + ( &
1422  + w(i_pgdep ) & ! [prod] v->g
1423  + w(i_pgacw ) & ! [prod] c->g
1424  + w(i_piacr_g) & ! [prod] r->g
1425  + w(i_psacr_g) & ! [prod] r->g
1426  + w(i_pgacr ) & ! [prod] r->g
1427  + w(i_pgfrz ) & ! [prod] r->g
1428  + w(i_praci_g) & ! [prod] i->g
1429  + w(i_pgaci ) & ! [prod] i->g
1430  + w(i_pgaut ) & ! [prod] s->g
1431  + w(i_pracs ) & ! [prod] s->g
1432  - w(i_pgsub ) & ! [loss] g->v
1433  ) * ice
1434 
1435  tend(i_qc) = max( tend(i_qc), -w(i_dqc_dt) )
1436  tend(i_qr) = max( tend(i_qr), -w(i_dqr_dt) )
1437  tend(i_qi) = max( tend(i_qi), -w(i_dqi_dt) )
1438  tend(i_qs) = max( tend(i_qs), -w(i_dqs_dt) )
1439  tend(i_qg) = max( tend(i_qg), -w(i_dqg_dt) )
1440 
1441  tend(i_qv) = - ( tend(i_qc) &
1442  + tend(i_qr) &
1443  + tend(i_qi) &
1444  + tend(i_qs) &
1445  + tend(i_qg) )
1446 
1447  qtrc_t(k,i,j,i_qv) = tend(i_qv)
1448  qtrc_t(k,i,j,i_qc) = tend(i_qc)
1449  qtrc_t(k,i,j,i_qr) = tend(i_qr)
1450  qtrc_t(k,i,j,i_qi) = tend(i_qi)
1451  qtrc_t(k,i,j,i_qs) = tend(i_qs)
1452  qtrc_t(k,i,j,i_qg) = tend(i_qg)
1453 
1454  qtrc0(k,i,j,i_qv) = qtrc0(k,i,j,i_qv) + qtrc_t(k,i,j,i_qv) * dt
1455  qtrc0(k,i,j,i_qc) = qtrc0(k,i,j,i_qc) + qtrc_t(k,i,j,i_qc) * dt
1456  qtrc0(k,i,j,i_qr) = qtrc0(k,i,j,i_qr) + qtrc_t(k,i,j,i_qr) * dt
1457  qtrc0(k,i,j,i_qi) = qtrc0(k,i,j,i_qi) + qtrc_t(k,i,j,i_qi) * dt
1458  qtrc0(k,i,j,i_qs) = qtrc0(k,i,j,i_qs) + qtrc_t(k,i,j,i_qs) * dt
1459  qtrc0(k,i,j,i_qg) = qtrc0(k,i,j,i_qg) + qtrc_t(k,i,j,i_qg) * dt
1460  enddo
1461  enddo
1462  enddo
1463 
1464  !-----< ice generation process >-----
1465  if ( mp_doexpricit_icegen ) then
1466  do j = js, je
1467  do i = is, ie
1468  do k = ks, ke
1469  ! store to work
1470  dens = dens0(k,i,j)
1471  temp = temp0(k,i,j)
1472  pres = pres0(k,i,j)
1473 
1474  ! saturation ratio S
1475  sice = qtrc0(k,i,j,i_qv) / max( qsati(k,i,j), eps )
1476 
1477  rdens = 1.0_rp / dens
1478  temc = temp - tem00
1479  rhoqi = dens * qtrc0(k,i,j,i_qi)
1480 
1481  da = ( da0 + dda_dt * temc )
1482  kd = ( dw0 + ddw_dt * temc ) * pre00 / pres
1483 
1484  giv = 1.0_rp / ( lhsex(k,i,j)/(da*temp) * ( lhsex(k,i,j)/(rvap*temp) - 1.0_rp ) + 1.0_rp/(kd*dens*qsati(k,i,j)) )
1485 
1486  ! [Pidep] deposition/sublimation : v->i or i->v
1487  sw = ( 0.5_rp + sign(0.5_rp, 0.0_rp - temc ) ) &
1488  * ( 0.5_rp + sign(0.5_rp, rhoqi ) ) ! if T < 0C & ice exists, sw=1
1489 
1490  rhoqi = max(rhoqi,eps)
1491 
1492  xni = min( max( 5.38e7_rp * rhoqi**0.75_rp, 1.e3_rp ), 1.e6_rp )
1493  xmi = rhoqi / xni
1494  di = min( di_a * sqrt(xmi), di_max )
1495 
1496  pidep = sw * 4.0_rp * di * xni * rdens * ( sice-1.0_rp ) * giv
1497 
1498  ! [Pigen] ice nucleation : v->i
1499  sw = ( 0.5_rp + sign(0.5_rp, 0.0_rp - temc ) ) &
1500  * ( 0.5_rp + sign(0.5_rp, sice - 1.0_rp ) ) ! if T < 0C & Satulated, sw=1
1501 
1502  ni0 = max( exp(-0.1_rp*temc), 1.0_rp ) * 1000.0_rp
1503  qi0 = 4.92e-11 * ni0**1.33_rp * rdens
1504 
1505  pigen = sw * max( min( qi0-qtrc0(k,i,j,i_qi), qtrc0(k,i,j,i_qv)-qsati(k,i,j) ), 0.0_rp ) / dt
1506 
1507  ! update Qv, Qi with limiter
1508  dq = ( pigen + pidep ) * dt
1509  qtmp = qtrc0(k,i,j,i_qv) - dq
1510  if ( dq > 0.0_rp ) then ! v->i
1511  qtmp = max( qtmp, qsati(k,i,j) )
1512  else ! i->v
1513  qtmp = min( qtmp, qsati(k,i,j) )
1514  endif
1515  dq = qtrc0(k,i,j,i_qv) - qtmp
1516  qtmp = qtrc0(k,i,j,i_qi) + dq
1517  qtmp = max( qtmp, 0.0_rp )
1518  dq = qtmp - qtrc0(k,i,j,i_qi)
1519  qtrc0(k,i,j,i_qi) = qtrc0(k,i,j,i_qi) + dq
1520  qtrc0(k,i,j,i_qv) = qtrc0(k,i,j,i_qv) - dq
1521  qtrc_t(k,i,j,i_qi) = qtrc_t(k,i,j,i_qi) + dq /dt
1522  qtrc_t(k,i,j,i_qv) = qtrc_t(k,i,j,i_qv) - dq /dt
1523 
1524  ! [Pifrz/Pimlt] freezing and melting
1525 
1526  ! homogenious freezing : c->i
1527  sw = ( 0.5_rp - sign(0.5_rp, temc + 40.0_rp ) ) ! if T < -40C, sw=1
1528 
1529  dq = sw * qtrc0(k,i,j,i_qc)
1530  qtrc0(k,i,j,i_qi) = qtrc0(k,i,j,i_qi) + dq
1531  qtrc0(k,i,j,i_qc) = qtrc0(k,i,j,i_qc) - dq
1532  qtrc_t(k,i,j,i_qi) = qtrc_t(k,i,j,i_qi) + dq /dt
1533  qtrc_t(k,i,j,i_qc) = qtrc_t(k,i,j,i_qc) - dq /dt
1534 
1535  ! heteroginous freezing : c->i
1536  sw = ( 0.5_rp + sign(0.5_rp, temc + 40.0_rp ) ) &
1537  * ( 0.5_rp + sign(0.5_rp, 0.0_rp - temc ) ) ! if -40C < T < 0C, sw=1
1538 
1539  dq = sw * ( dens / dwatr * qtrc0(k,i,j,i_qc)**2 / (nc(k,i,j)*1.e+6) ) &
1540  * b_frz * ( exp(-a_frz*temc) - 1.0_rp ) * dt
1541  qtmp = qtrc0(k,i,j,i_qc) - dq
1542  qtmp = max( qtmp, 0.0_rp )
1543  dq = qtrc0(k,i,j,i_qc) - qtmp
1544  qtrc0(k,i,j,i_qi) = qtrc0(k,i,j,i_qi) + dq
1545  qtrc0(k,i,j,i_qc) = qtrc0(k,i,j,i_qc) - dq
1546  qtrc_t(k,i,j,i_qi) = qtrc_t(k,i,j,i_qi) + dq /dt
1547  qtrc_t(k,i,j,i_qc) = qtrc_t(k,i,j,i_qc) - dq /dt
1548 
1549  ! melting : i->c
1550  sw = ( 0.5_rp - sign(0.5_rp, 0.0_rp - temc ) ) ! if 0C < T, sw=1
1551 
1552  dq = - sw * qtrc0(k,i,j,i_qi)
1553  qtrc0(k,i,j,i_qi) = qtrc0(k,i,j,i_qi) + dq
1554  qtrc0(k,i,j,i_qc) = qtrc0(k,i,j,i_qc) - dq
1555  qtrc_t(k,i,j,i_qi) = qtrc_t(k,i,j,i_qi) + dq /dt
1556  qtrc_t(k,i,j,i_qc) = qtrc_t(k,i,j,i_qc) - dq /dt
1557 
1558  enddo
1559  enddo
1560  enddo
1561  endif
1562 
1563  do j = js, je
1564  do i = is, ie
1565  do k = ks, ke
1566  rhoe_t(k,i,j) = - dens0(k,i,j) * ( lhvex(k,i,j) * qtrc_t(k,i,j,i_qv) &
1567  - lhfex(k,i,j) * qtrc_t(k,i,j,i_qi) &
1568  - lhfex(k,i,j) * qtrc_t(k,i,j,i_qs) &
1569  - lhfex(k,i,j) * qtrc_t(k,i,j,i_qg) )
1570 
1571  rhoe0(k,i,j) = rhoe0(k,i,j) + rhoe_t(k,i,j) * dt
1572  enddo
1573  enddo
1574  enddo
1575 
1576 ! do ip = 1, w_nmax
1577 ! if ( w_histid(ip) > 0 ) then
1578 
1579 ! call HIST_query( w_histid(ip), & ! [IN]
1580 ! do_put ) ! [OUT]
1581 
1582 ! if ( do_put ) then
1583 ! do j = JS, JE
1584 ! do i = IS, IE
1585 ! do k = KS, KE
1586 ! work3D(k,i,j) = w(ip)
1587 ! enddo
1588 ! enddo
1589 ! enddo
1590 
1591 ! if( IO_L ) write(IO_FID_LOG,*) w_name(ip), "MAX/MIN:", &
1592 ! maxval(work3D(KS:KE,IS:IE,JS:JE)), minval(work3D(KS:KE,IS:IE,JS:JE))
1593 
1594 
1595 ! call HIST_put( w_histid(ip), work3D(:,:,:), w_zinterp(ip) )
1596 ! endif
1597 ! endif
1598 ! enddo
1599 
1600 ! if( IO_L ) write(IO_FID_LOG,*) "tend QV MAX/MIN:", maxval(tend(I_QV,:)), minval(tend(I_QV,:))
1601 ! if( IO_L ) write(IO_FID_LOG,*) "tend QC MAX/MIN:", maxval(tend(I_QC,:)), minval(tend(I_QC,:))
1602 ! if( IO_L ) write(IO_FID_LOG,*) "tend QR MAX/MIN:", maxval(tend(I_QR,:)), minval(tend(I_QR,:))
1603 ! if( IO_L ) write(IO_FID_LOG,*) "tend QI MAX/MIN:", maxval(tend(I_QI,:)), minval(tend(I_QI,:))
1604 ! if( IO_L ) write(IO_FID_LOG,*) "tend QS MAX/MIN:", maxval(tend(I_QS,:)), minval(tend(I_QS,:))
1605 ! if( IO_L ) write(IO_FID_LOG,*) "tend QG MAX/MIN:", maxval(tend(I_QG,:)), minval(tend(I_QG,:))
1606 !
1607 ! if( IO_L ) write(IO_FID_LOG,*) "QV MAX/MIN:", maxval(QTRC0(:,:,:,I_QV)), minval(QTRC0(:,:,:,I_QV))
1608 ! if( IO_L ) write(IO_FID_LOG,*) "QC MAX/MIN:", maxval(QTRC0(:,:,:,I_QC)), minval(QTRC0(:,:,:,I_QC))
1609 ! if( IO_L ) write(IO_FID_LOG,*) "QR MAX/MIN:", maxval(QTRC0(:,:,:,I_QR)), minval(QTRC0(:,:,:,I_QR))
1610 ! if( IO_L ) write(IO_FID_LOG,*) "QI MAX/MIN:", maxval(QTRC0(:,:,:,I_QI)), minval(QTRC0(:,:,:,I_QI))
1611 ! if( IO_L ) write(IO_FID_LOG,*) "QS MAX/MIN:", maxval(QTRC0(:,:,:,I_QS)), minval(QTRC0(:,:,:,I_QS))
1612 ! if( IO_L ) write(IO_FID_LOG,*) "QG MAX/MIN:", maxval(QTRC0(:,:,:,I_QG)), minval(QTRC0(:,:,:,I_QG))
1613 
1614 ! if( IO_L ) write(IO_FID_LOG,*) "QV_t MAX/MIN:", maxval(QTRC_t(:,:,:,I_QV)), minval(QTRC_t(:,:,:,I_QV))
1615 ! if( IO_L ) write(IO_FID_LOG,*) "QC_t MAX/MIN:", maxval(QTRC_t(:,:,:,I_QC)), minval(QTRC_t(:,:,:,I_QC))
1616 ! if( IO_L ) write(IO_FID_LOG,*) "QR_t MAX/MIN:", maxval(QTRC_t(:,:,:,I_QR)), minval(QTRC_t(:,:,:,I_QR))
1617 ! if( IO_L ) write(IO_FID_LOG,*) "QI_t MAX/MIN:", maxval(QTRC_t(:,:,:,I_QI)), minval(QTRC_t(:,:,:,I_QI))
1618 ! if( IO_L ) write(IO_FID_LOG,*) "QS_t MAX/MIN:", maxval(QTRC_t(:,:,:,I_QS)), minval(QTRC_t(:,:,:,I_QS))
1619 ! if( IO_L ) write(IO_FID_LOG,*) "QG_t MAX/MIN:", maxval(QTRC_t(:,:,:,I_QG)), minval(QTRC_t(:,:,:,I_QG))
1620 
1621 ! if( IO_L ) write(IO_FID_LOG,*) "RHOE0 MAX/MIN:", maxval(RHOE0(:,:,:)), minval(RHOE0(:,:,:))
1622 ! if( IO_L ) write(IO_FID_LOG,*) "RHOE_t MAX/MIN:", maxval(RHOE_t(:,:,:)), minval(RHOE_t(:,:,:))
1623 
1624  call prof_rapend ('MP_tomita08', 3)
1625 
1626  return
1627  end subroutine mp_tomita08
1628 
1629  !-----------------------------------------------------------------------------
1631  subroutine mp_tomita08_vterm( &
1632  vterm, &
1633  DENS0, &
1634  TEMP0, &
1635  QTRC0 )
1636  use scale_const, only: &
1637  tem00 => const_tem00
1638  implicit none
1639 
1640  real(RP), intent(out) :: vterm(ka,ia,ja,qa)
1641  real(RP), intent(in) :: DENS0(ka,ia,ja)
1642  real(RP), intent(in) :: TEMP0(ka,ia,ja)
1643  real(RP), intent(in) :: QTRC0(ka,ia,ja,qa)
1644 
1645  real(RP) :: dens
1646  real(RP) :: temc
1647  real(RP) :: q(qa)
1648 
1649  real(RP) :: rho_fact ! density factor
1650 
1651  real(RP) :: RLMDr, RLMDs, RLMDg
1652  real(RP) :: RLMDr_dr, RLMDs_ds, RLMDg_dg
1653 
1654  !---< Roh and Satoh (2014) >---
1655  real(RP) :: tems, rhoqs, Xs2
1656  real(RP) :: MOMs_0bs, RMOMs_Vt
1657  real(RP) :: coef_at(4), coef_bt(4)
1658  real(RP) :: loga_, b_, nm
1659 
1660  real(RP) :: zerosw
1661  integer :: k, i, j, iq
1662  !---------------------------------------------------------------------------
1663 
1664  do j = js, je
1665  do i = is, ie
1666  do k = ks, ke
1667  ! store to work
1668  dens = dens0(k,i,j)
1669  temc = temp0(k,i,j) - tem00
1670  do iq = i_qv, i_qg
1671  q(iq) = qtrc0(k,i,j,iq)
1672  enddo
1673 
1674  rho_fact = sqrt( dens00 / dens )
1675 
1676  ! slope parameter lambda (Rain)
1677  zerosw = 0.5_rp - sign(0.5_rp, q(i_qr) - 1.e-12_rp )
1678  rlmdr = sqrt(sqrt( dens * q(i_qr) / ( ar * n0r * gam_1br ) + zerosw )) * ( 1.0_rp - zerosw )
1679 
1680  rlmdr_dr = sqrt( rlmdr ) ! **Dr
1681 
1682  ! slope parameter lambda (Snow)
1683  zerosw = 0.5_rp - sign(0.5_rp, q(i_qs) - 1.e-12_rp )
1684  rlmds = sqrt(sqrt( dens * q(i_qs) / ( as * n0s * gam_1bs ) + zerosw )) * ( 1.0_rp - zerosw )
1685 
1686  rlmds_ds = sqrt( sqrt(rlmds) ) ! **Ds
1687  rmoms_vt = gam_1bsds / gam_1bs * rlmds_ds
1688 
1689  !---< modification by Roh and Satoh (2014) >---
1690  ! bimodal size distribution of snow
1691  rhoqs = dens * q(i_qs)
1692  zerosw = 0.5_rp - sign(0.5_rp, rhoqs - 1.e-12_rp )
1693  xs2 = rhoqs / as * ( 1.0_rp - zerosw )
1694 
1695  tems = min( -0.1_rp, temc )
1696  coef_at(1) = coef_a( 1) + tems * ( coef_a( 2) + tems * ( coef_a( 5) + tems * coef_a( 9) ) )
1697  coef_at(2) = coef_a( 3) + tems * ( coef_a( 4) + tems * coef_a( 7) )
1698  coef_at(3) = coef_a( 6) + tems * coef_a( 8)
1699  coef_at(4) = coef_a(10)
1700  coef_bt(1) = coef_b( 1) + tems * ( coef_b( 2) + tems * ( coef_b( 5) + tems * coef_b( 9) ) )
1701  coef_bt(2) = coef_b( 3) + tems * ( coef_b( 4) + tems * coef_b( 7) )
1702  coef_bt(3) = coef_b( 6) + tems * coef_b( 8)
1703  coef_bt(4) = coef_b(10)
1704  ! 0 + Bs(=2) moment
1705  nm = 2.0_rp
1706  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
1707  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
1708  moms_0bs = 10.0_rp**loga_ * xs2**b_
1709  ! Bs(=2) + Ds(=0.25) moment
1710  nm = 2.25_rp
1711  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
1712  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
1713  rmoms_vt = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ / ( moms_0bs + zerosw ) &
1714  + ( 1.0_rp-sw_roh2014 ) * rmoms_vt
1715 
1716  ! slope parameter lambda (Graupel)
1717  zerosw = 0.5_rp - sign(0.5_rp, q(i_qg) - 1.e-12_rp )
1718  rlmdg = sqrt(sqrt( dens * q(i_qg) / ( ag * n0g * gam_1bg ) + zerosw )) * ( 1.0_rp - zerosw )
1719 
1720  rlmdg_dg = sqrt( rlmdg ) ! **Dg
1721 
1722  !---< terminal velocity >
1723  zerosw = 0.5_rp + sign(0.5_rp, q(i_qi) - 1.e-8_rp )
1724  vterm(k,i,j,i_qv) = 0.0_rp
1725  vterm(k,i,j,i_qc) = 0.0_rp
1726  vterm(k,i,j,i_qi) = -3.29_rp * ( dens * q(i_qi) * zerosw )**0.16_rp
1727  vterm(k,i,j,i_qr) = -cr * rho_fact * gam_1brdr / gam_1br * rlmdr_dr
1728  vterm(k,i,j,i_qs) = -cs * rho_fact * rmoms_vt
1729  vterm(k,i,j,i_qg) = -cg * rho_fact * gam_1bgdg / gam_1bg * rlmdg_dg
1730  enddo
1731  enddo
1732  enddo
1733 
1734  do j = js, je
1735  do i = is, ie
1736  vterm( 1:ks-1,i,j,:) = 0.0_rp
1737  vterm(ke+1:ka ,i,j,:) = 0.0_rp
1738  enddo
1739  enddo
1740 
1741  return
1742  end subroutine mp_tomita08_vterm
1743 
1744  !-----------------------------------------------------------------------------
1745  subroutine mp_tomita08_bergeronparam( &
1746  temp, &
1747  a1, &
1748  a2, &
1749  ma2 )
1750  use scale_const, only: &
1751  tem00 => const_tem00
1752  implicit none
1753 
1754  real(RP), intent(in) :: temp(ka,ia,ja)
1755  real(RP), intent(out) :: a1 (ka,ia,ja)
1756  real(RP), intent(out) :: a2 (ka,ia,ja)
1757  real(RP), intent(out) :: ma2 (ka,ia,ja)
1758 
1759  real(RP) :: a1_tab(32)
1760  real(RP) :: a2_tab(32)
1761 
1762  data a1_tab / 0.0001e-7_rp, 0.7939e-7_rp, 0.7841e-6_rp, 0.3369e-5_rp, 0.4336e-5_rp, &
1763  0.5285e-5_rp, 0.3728e-5_rp, 0.1852e-5_rp, 0.2991e-6_rp, 0.4248e-6_rp, &
1764  0.7434e-6_rp, 0.1812e-5_rp, 0.4394e-5_rp, 0.9145e-5_rp, 0.1725e-4_rp, &
1765  0.3348e-4_rp, 0.1725e-4_rp, 0.9175e-5_rp, 0.4412e-5_rp, 0.2252e-5_rp, &
1766  0.9115e-6_rp, 0.4876e-6_rp, 0.3473e-6_rp, 0.4758e-6_rp, 0.6306e-6_rp, &
1767  0.8573e-6_rp, 0.7868e-6_rp, 0.7192e-6_rp, 0.6513e-6_rp, 0.5956e-6_rp, &
1768  0.5333e-6_rp, 0.4834e-6_rp /
1769 
1770  data a2_tab / 0.0100_rp, 0.4006_rp, 0.4831_rp, 0.5320_rp, 0.5307_rp, &
1771  0.5319_rp, 0.5249_rp, 0.4888_rp, 0.3849_rp, 0.4047_rp, &
1772  0.4318_rp, 0.4771_rp, 0.5183_rp, 0.5463_rp, 0.5651_rp, &
1773  0.5813_rp, 0.5655_rp, 0.5478_rp, 0.5203_rp, 0.4906_rp, &
1774  0.4447_rp, 0.4126_rp, 0.3960_rp, 0.4149_rp, 0.4320_rp, &
1775  0.4506_rp, 0.4483_rp, 0.4460_rp, 0.4433_rp, 0.4413_rp, &
1776  0.4382_rp, 0.4361_rp /
1777 
1778  real(RP) :: temc
1779  integer :: itemc
1780  real(RP) :: fact
1781 
1782  integer :: k, i, j
1783  !---------------------------------------------------------------------------
1784 
1785  do j = js, je
1786  do i = is, ie
1787  do k = ks, ke
1788  temc = min( max( temp(k,i,j)-tem00, -30.99_rp ), 0.0_rp ) ! 0C <= T < 31C
1789  itemc = int( -temc ) + 1 ! 1 <= iT <= 31
1790  fact = - ( temc + real(itemc-1,kind=8) )
1791 
1792  a1(k,i,j) = ( 1.0_rp-fact ) * a1_tab(itemc ) &
1793  + ( fact ) * a1_tab(itemc+1)
1794 
1795  a2(k,i,j) = ( 1.0_rp-fact ) * a2_tab(itemc ) &
1796  + ( fact ) * a2_tab(itemc+1)
1797 
1798  ma2(k,i,j) = 1.0_rp - a2(k,i,j)
1799 
1800  a1(k,i,j) = a1(k,i,j) * 1.e-3_rp**ma2(k,i,j) ! [g->kg]
1801  enddo
1802  enddo
1803  enddo
1804 
1805  return
1806  end subroutine mp_tomita08_bergeronparam
1807 
1808  !-----------------------------------------------------------------------------
1810  subroutine atmos_phy_mp_tomita08_cloudfraction( &
1811  cldfrac, &
1812  QTRC )
1813  use scale_grid_index
1814  use scale_tracer, only: &
1815  qad => qa
1816  implicit none
1817 
1818  real(RP), intent(out) :: cldfrac(ka,ia,ja)
1819  real(RP), intent(in) :: QTRC (ka,ia,ja,qad)
1820 
1821  real(RP) :: qcriteria = 0.005e-3_rp ! 0.005g/kg, Tompkins & Craig
1822 
1823  real(RP) :: qhydro
1824  integer :: k, i, j, iq
1825  !---------------------------------------------------------------------------
1826 
1827  do j = js, je
1828  do i = is, ie
1829  do k = ks, ke
1830  qhydro = 0.0_rp
1831  do iq = 1, mp_qa
1832  qhydro = qhydro + qtrc(k,i,j,i_mp2all(iq))
1833  enddo
1834  cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-qcriteria)
1835  enddo
1836  enddo
1837  enddo
1838 
1839  return
1840  end subroutine atmos_phy_mp_tomita08_cloudfraction
1841 
1842  !-----------------------------------------------------------------------------
1844  subroutine atmos_phy_mp_tomita08_effectiveradius( &
1845  Re, &
1846  QTRC0, &
1847  DENS0, &
1848  TEMP0 )
1849  use scale_grid_index
1850  use scale_const, only: &
1851  tem00 => const_tem00
1852  use scale_tracer, only: &
1853  qad => qa, &
1854  mp_qad => mp_qa
1855  implicit none
1856 
1857  real(RP), intent(out) :: Re (ka,ia,ja,mp_qad) ! effective radius [cm]
1858  real(RP), intent(in) :: QTRC0(ka,ia,ja,qad) ! tracer mass concentration [kg/kg]
1859  real(RP), intent(in) :: DENS0(ka,ia,ja) ! density [kg/m3]
1860  real(RP), intent(in) :: TEMP0(ka,ia,ja) ! temperature [K]
1861 
1862  real(RP) :: dens
1863  real(RP) :: temc
1864  real(RP) :: q(qa)
1865  real(RP) :: RLMDr, RLMDs, RLMDg
1866 
1867  real(RP), parameter :: um2cm = 100.0_rp
1868 
1869  !---< Roh and Satoh (2014) >---
1870  real(RP) :: tems, rhoqs, Xs2
1871  real(RP) :: MOMs_2, MOMs_1bs
1872  real(RP) :: coef_at(4), coef_bt(4)
1873  real(RP) :: loga_, b_, nm
1874 
1875  real(RP) :: zerosw
1876  integer :: k, i, j, iq
1877  !---------------------------------------------------------------------------
1878 
1879  re(:,:,:,i_mp_qc) = re_qc * um2cm
1880  re(:,:,:,i_mp_qi) = re_qi * um2cm
1881 
1882  do j = js, je
1883  do i = is, ie
1884  do k = ks, ke
1885  dens = dens0(k,i,j)
1886  temc = temp0(k,i,j) - tem00
1887  do iq = i_qv, i_qg
1888  q(iq) = qtrc0(k,i,j,iq)
1889  enddo
1890 
1891  ! slope parameter lambda
1892  zerosw = 0.5_rp - sign(0.5_rp, q(i_qr) - 1.e-12_rp )
1893  rlmdr = sqrt(sqrt( dens * q(i_qr) / ( ar * n0r * gam_1br ) + zerosw )) * ( 1.0_rp - zerosw )
1894 
1895  ! Effective radius is defined by r3m/r2m = 1.5/lambda
1896  re(k,i,j,i_mp_qr) = 1.5_rp * rlmdr * um2cm
1897 
1898  zerosw = 0.5_rp - sign(0.5_rp, q(i_qs) - 1.e-12_rp )
1899  rlmds = sqrt(sqrt( dens * q(i_qs) / ( as * n0s * gam_1bs ) + zerosw )) * ( 1.0_rp - zerosw )
1900 
1901  !---< modification by Roh and Satoh (2014) >---
1902  ! bimodal size distribution of snow
1903  rhoqs = dens * q(i_qs)
1904  zerosw = 0.5_rp - sign(0.5_rp, rhoqs - 1.e-12_rp )
1905  xs2 = rhoqs / as * ( 1.0_rp - zerosw )
1906 
1907  tems = min( -0.1_rp, temc )
1908  coef_at(1) = coef_a( 1) + tems * ( coef_a( 2) + tems * ( coef_a( 5) + tems * coef_a( 9) ) )
1909  coef_at(2) = coef_a( 3) + tems * ( coef_a( 4) + tems * coef_a( 7) )
1910  coef_at(3) = coef_a( 6) + tems * coef_a( 8)
1911  coef_at(4) = coef_a(10)
1912  coef_bt(1) = coef_b( 1) + tems * ( coef_b( 2) + tems * ( coef_b( 5) + tems * coef_b( 9) ) )
1913  coef_bt(2) = coef_b( 3) + tems * ( coef_b( 4) + tems * coef_b( 7) )
1914  coef_bt(3) = coef_b( 6) + tems * coef_b( 8)
1915  coef_bt(4) = coef_b(10)
1916  ! 2nd moment
1917  moms_2 = xs2
1918  ! 1 + Bs(=2) moment
1919  nm = 3.0_rp
1920  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
1921  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
1922  moms_1bs = 10.0_rp**loga_ * xs2**b_
1923 
1924  re(k,i,j,i_mp_qs) = ( sw_roh2014 ) * 0.5_rp * moms_1bs / ( moms_2 + zerosw ) * um2cm &
1925  + ( 1.0_rp-sw_roh2014 ) * 1.5_rp * rlmds * um2cm
1926 
1927  zerosw = 0.5_rp - sign(0.5_rp, q(i_qg) - 1.e-12_rp )
1928  rlmdg = sqrt(sqrt( dens * q(i_qg) / ( ag * n0g * gam_1bg ) + zerosw )) * ( 1.0_rp - zerosw )
1929 
1930  re(k,i,j,i_mp_qg) = 1.5_rp * rlmdg * um2cm
1931  enddo
1932  enddo
1933  enddo
1934 
1935  return
1936  end subroutine atmos_phy_mp_tomita08_effectiveradius
1937 
1938  !-----------------------------------------------------------------------------
1940  subroutine atmos_phy_mp_tomita08_mixingratio( &
1941  Qe, &
1942  QTRC0 )
1943  use scale_grid_index
1944  use scale_tracer, only: &
1945  qad => qa, &
1946  mp_qad => mp_qa
1947  implicit none
1948 
1949  real(RP), intent(out) :: Qe (ka,ia,ja,mp_qad) ! mixing ratio of each cateory [kg/kg]
1950  real(RP), intent(in) :: QTRC0(ka,ia,ja,qad) ! tracer mass concentration [kg/kg]
1951 
1952  integer :: ihydro
1953  !---------------------------------------------------------------------------
1954 
1955  do ihydro = 1, mp_qa
1956  qe(:,:,:,ihydro) = qtrc0(:,:,:,i_mp2all(ihydro))
1957  enddo
1958 
1959  return
1960  end subroutine atmos_phy_mp_tomita08_mixingratio
1961 
1962 end module scale_atmos_phy_mp_tomita08
1963 !-------------------------------------------------------------------------------
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_history::hist_reg
subroutine, public hist_reg(itemid, zinterp, item, desc, unit, ndim, xdim, ydim, zdim)
Register/Append variable to history file.
Definition: scale_history.F90:496
scale_atmos_saturation
module ATMOSPHERE / Saturation adjustment
Definition: scale_atmos_sub_saturation.F90:16
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:234
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_cloudfraction
subroutine, public atmos_phy_mp_tomita08_cloudfraction(cldfrac, QTRC)
Calculate Cloud Fraction.
Definition: scale_atmos_phy_mp_tomita08.F90:1813
scale_atmos_phy_mp_common::atmos_phy_mp_precipitation
subroutine, public atmos_phy_mp_precipitation(flux_rain, flux_snow, DENS, MOMZ, MOMX, MOMY, RHOE, QTRC, vterm, temp, dt)
precipitation transport
Definition: scale_atmos_phy_mp_common.F90:657
scale_const::const_dwatr
real(rp), parameter, public const_dwatr
density of water [kg/m3]
Definition: scale_const.F90:87
scale_time::time_dtsec_atmos_phy_mp
real(dp), public time_dtsec_atmos_phy_mp
time interval of physics(microphysics) [sec]
Definition: scale_time.F90:41
scale_stdio::io_l
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:59
scale_specfunc::sf_gamma
real(rp) function, public sf_gamma(x)
Gamma function.
Definition: scale_specfunc.F90:48
scale_const::const_cl
real(rp), parameter, public const_cl
specific heat (liquid water) [J/kg/K]
Definition: scale_const.F90:68
scale_atmos_phy_mp_tomita08::atmos_phy_mp_dens
real(rp), dimension(mp_qa), target, public atmos_phy_mp_dens
Definition: scale_atmos_phy_mp_tomita08.F90:45
scale_stdio
module STDIO
Definition: scale_stdio.F90:12
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_const::const_tem00
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
Definition: scale_const.F90:95
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_mixingratio
subroutine, public atmos_phy_mp_tomita08_mixingratio(Qe, QTRC0)
Calculate mixing ratio of each category.
Definition: scale_atmos_phy_mp_tomita08.F90:1943
scale_const::const_dice
real(rp), parameter, public const_dice
density of ice [kg/m3]
Definition: scale_const.F90:88
scale_atmos_phy_mp_common
module ATMOSPHERE / Physics Cloud Microphysics - Common
Definition: scale_atmos_phy_mp_common.F90:16
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_history::hist_query
subroutine, public hist_query(itemid, answer)
Check time to putting data.
Definition: scale_history.F90:604
scale_tracer
module TRACER
Definition: scale_tracer.F90:18
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_tracer::mp_qa
integer, public mp_qa
Definition: scale_tracer.F90:98
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:93
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:48
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_time
module TIME
Definition: scale_time.F90:15
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_specfunc
module SPECFUNC
Definition: scale_specfunc.F90:13
dc_log::log
subroutine, public log(type, message)
Definition: dc_log.f90:133
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:65
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_atmos_phy_mp_common::atmos_phy_mp_negative_fixer
subroutine, public atmos_phy_mp_negative_fixer(DENS, RHOT, QTRC)
Negative fixer.
Definition: scale_atmos_phy_mp_common.F90:64
scale_grid
module GRID (cartesian)
Definition: scale_grid_cartesian.F90:17
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level)
Start raptime.
Definition: scale_prof.F90:132
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08
subroutine, public atmos_phy_mp_tomita08(DENS, MOMZ, MOMX, MOMY, RHOT, QTRC, CCN, EVAPORATE, SFLX_rain, SFLX_snow)
Cloud Microphysics.
Definition: scale_atmos_phy_mp_tomita08.F90:496
scale_prof
module profiler
Definition: scale_prof.F90:10
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:36
scale_atmos_thermodyn
module ATMOSPHERE / Thermodynamics
Definition: scale_atmos_sub_thermodyn.F90:19
scale_stdio::io_lnml
logical, public io_lnml
output log or not? (for namelist, this process)
Definition: scale_stdio.F90:60
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_grid::grid_cdz
real(rp), dimension(:), allocatable, public grid_cdz
z-length of control volume [m]
Definition: scale_grid_cartesian.F90:56
scale_history
module HISTORY
Definition: scale_history.F90:16
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:34
scale_stdio::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_stdio.F90:55
scale_atmos_phy_mp_common::atmos_phy_mp_saturation_adjustment
subroutine, public atmos_phy_mp_saturation_adjustment(RHOE_t, QTRC_t, RHOE0, QTRC0, DENS0, flag_liquid)
Saturation adjustment.
Definition: scale_atmos_phy_mp_common.F90:129
scale_stdio::io_fid_log
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
scale_tracer_tomita08
module TRACER / tomita08
Definition: scale_tracer_tomita08.F90:14
scale_atmos_phy_mp_tomita08
module ATMOSPHERE / Physics Cloud Microphysics
Definition: scale_atmos_phy_mp_tomita08.F90:16
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level)
Save raptime.
Definition: scale_prof.F90:178
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:38
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_setup
subroutine, public atmos_phy_mp_tomita08_setup(MP_TYPE)
Setup.
Definition: scale_atmos_phy_mp_tomita08.F90:277
scale_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_effectiveradius
subroutine, public atmos_phy_mp_tomita08_effectiveradius(Re, QTRC0, DENS0, TEMP0)
Calculate Effective Radius.
Definition: scale_atmos_phy_mp_tomita08.F90:1849
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Generated by   doxygen 1.8.13