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  ! For the calculation of optically effective volume
402  atmos_phy_mp_dens(i_mp_qc) = dens_w
403  atmos_phy_mp_dens(i_mp_qr) = dens_w
404  atmos_phy_mp_dens(i_mp_qi) = dens_i
405  atmos_phy_mp_dens(i_mp_qs) = dens_i
406  atmos_phy_mp_dens(i_mp_qg) = dens_i
407 
408  !--- empirical coefficients A, B, C, D
409  ar = pi * dens_w / 6.0_rp
410  as = pi * dens_s / 6.0_rp
411  ag = pi * dens_g / 6.0_rp
412 
413  br = 3.0_rp
414  bs = 3.0_rp
415  bg = 3.0_rp
416 
417  cg = sqrt( ( 4.0_rp * dens_g * grav ) / ( 3.0_rp * dens00 * drag_g ) )
418 
419  dr = 0.50_rp
420  ds = 0.25_rp
421  dg = 0.50_rp
422 
423  if ( enable_roh2014 ) then ! overwrite parameters
424  sw_roh2014 = 1.0_rp
425  n0g = 4.e8_rp
426  as = 0.069_rp
427  bs = 2.0_rp
428  esi = 0.25_rp
429  endif
430 
431  gam = 1.0_rp ! =0!
432  gam_2 = 1.0_rp ! =1!
433  gam_3 = 2.0_rp ! =2!
434 
435  gam_1br = sf_gamma( 1.0_rp + br ) ! = 4!
436  gam_2br = sf_gamma( 2.0_rp + br ) ! = 5!
437  gam_3br = sf_gamma( 3.0_rp + br ) ! = 6!
438  gam_3dr = sf_gamma( 3.0_rp + dr )
439  gam_6dr = sf_gamma( 6.0_rp + dr )
440  gam_1brdr = sf_gamma( 1.0_rp + br + dr )
441  gam_5dr_h = sf_gamma( 0.5_rp * (5.0_rp+dr) )
442 
443  gam_1bs = sf_gamma( 1.0_rp + bs ) ! = 4!
444  gam_2bs = sf_gamma( 2.0_rp + bs ) ! = 5!
445  gam_3bs = sf_gamma( 3.0_rp + bs ) ! = 6!
446  gam_3ds = sf_gamma( 3.0_rp + ds )
447  gam_1bsds = sf_gamma( 1.0_rp + bs + ds )
448  gam_5ds_h = sf_gamma( 0.5_rp * (5.0_rp+ds) )
449 
450  gam_1bg = sf_gamma( 1.0_rp + bg ) ! = 4!
451  gam_3dg = sf_gamma( 3.0_rp + dg )
452  gam_1bgdg = sf_gamma( 1.0_rp + bg + dg)
453  gam_5dg_h = sf_gamma( 0.5_rp * (5.0_rp+dg) )
454 
455  do j = js, je
456  do i = is, ie
457  nc_def(i,j) = autoconv_nc
458  enddo
459  enddo
460 
461  if ( mp_doexpricit_icegen ) then
462  only_liquid = .true.
463  sw_useicegen = 1.0_rp
464  else
465  only_liquid = .false.
466  sw_useicegen = 0.0_rp
467  endif
468 
469  if ( enable_kk2000 ) then
470  sw_kk2000 = 1.0_rp
471  else
472  sw_kk2000 = 0.0_rp
473  endif
474 
475  ! detailed tendency monitor
476  w_desc(:) = w_name(:)
477  do ip = 1, w_nmax
478  call hist_reg( w_histid(ip), w_zinterp(ip), w_name(ip), w_desc(ip), w_unit(ip), ndim=3 )
479  enddo
480 
481  return
482  end subroutine atmos_phy_mp_tomita08_setup
483 
484  !-----------------------------------------------------------------------------
486  subroutine atmos_phy_mp_tomita08( &
487  DENS, &
488  MOMZ, &
489  MOMX, &
490  MOMY, &
491  RHOT, &
492  QTRC, &
493  CCN, &
494  EVAPORATE, &
495  SFLX_rain, &
496  SFLX_snow )
497  use scale_grid_index
498  use scale_const, only: &
499  dwatr => const_dwatr, &
500  pi => const_pi
501  use scale_time, only: &
502  dt => time_dtsec_atmos_phy_mp
503  use scale_history, only: &
504  hist_in
505  use scale_tracer, only: &
506  qad => qa
507  use scale_atmos_thermodyn, only: &
508  thermodyn_rhoe => atmos_thermodyn_rhoe, &
509  thermodyn_rhot => atmos_thermodyn_rhot, &
510  thermodyn_temp_pres_e => atmos_thermodyn_temp_pres_e
511  use scale_atmos_phy_mp_common, only: &
512  mp_negative_fixer => atmos_phy_mp_negative_fixer, &
513  mp_precipitation => atmos_phy_mp_precipitation, &
514  mp_saturation_adjustment => atmos_phy_mp_saturation_adjustment
515  implicit none
516 
517  real(RP), intent(inout) :: dens(ka,ia,ja)
518  real(RP), intent(inout) :: MOMZ(ka,ia,ja)
519  real(RP), intent(inout) :: MOMX(ka,ia,ja)
520  real(RP), intent(inout) :: MOMY(ka,ia,ja)
521  real(RP), intent(inout) :: RHOT(ka,ia,ja)
522  real(RP), intent(inout) :: QTRC(ka,ia,ja,qad)
523  real(RP), intent(in) :: CCN(ka,ia,ja)
524  real(RP), intent(out) :: EVAPORATE(ka,ia,ja) ! number of evaporated cloud [/m3/s]
525  real(RP), intent(out) :: SFLX_rain(ia,ja)
526  real(RP), intent(out) :: SFLX_snow(ia,ja)
527 
528  real(RP) :: RHOE_t(ka,ia,ja)
529  real(RP) :: QTRC_t(ka,ia,ja,qa)
530  real(RP) :: RHOE (ka,ia,ja)
531  real(RP) :: temp (ka,ia,ja)
532  real(RP) :: pres (ka,ia,ja)
533 
534  real(RP) :: vterm (ka,ia,ja,qa) ! terminal velocity of each tracer [m/s]
535  real(RP) :: FLX_rain (ka,ia,ja)
536  real(RP) :: FLX_snow (ka,ia,ja)
537  real(RP) :: wflux_rain(ka,ia,ja)
538  real(RP) :: wflux_snow(ka,ia,ja)
539  real(RP) :: qc_before_satadj(ka,ia,ja)
540  integer :: step
541  integer :: k, i, j
542  !---------------------------------------------------------------------------
543 
544  if( io_l ) write(io_fid_log,*) '*** Physics step: Cloud microphysics(tomita08)'
545 
546  if ( mp_donegative_fixer ) then
547  call mp_negative_fixer( dens(:,:,:), & ! [INOUT]
548  rhot(:,:,:), & ! [INOUT]
549  qtrc(:,:,:,:) ) ! [INOUT]
550  endif
551 
552  call thermodyn_rhoe( rhoe(:,:,:), & ! [OUT]
553  rhot(:,:,:), & ! [IN]
554  qtrc(:,:,:,:) ) ! [IN]
555 
556  !##### MP Main #####
557  call mp_tomita08( rhoe_t(:,:,:), & ! [OUT]
558  qtrc_t(:,:,:,:), & ! [OUT]
559  rhoe(:,:,:), & ! [INOUT]
560  qtrc(:,:,:,:), & ! [INOUT]
561  ccn(:,:,:), & ! [IN]
562  dens(:,:,:) ) ! [IN]
563 
564  flx_rain(:,:,:) = 0.0_rp
565  flx_snow(:,:,:) = 0.0_rp
566 
567  vterm(:,:,:,:) = 0.0_rp
568 
569  if ( mp_doprecipitation ) then
570 
571  do step = 1, mp_nstep_sedimentation
572 
573  call thermodyn_temp_pres_e( temp(:,:,:), & ! [OUT]
574  pres(:,:,:), & ! [OUT]
575  dens(:,:,:), & ! [IN]
576  rhoe(:,:,:), & ! [IN]
577  qtrc(:,:,:,:) ) ! [IN]
578 
579  call mp_tomita08_vterm( vterm(:,:,:,:), & ! [OUT]
580  dens(:,:,:), & ! [IN]
581  temp(:,:,:), & ! [IN]
582  qtrc(:,:,:,:) ) ! [IN]
583 
584  call mp_precipitation( wflux_rain(:,:,:), & ! [OUT]
585  wflux_snow(:,:,:), & ! [OUT]
586  dens(:,:,:), & ! [INOUT]
587  momz(:,:,:), & ! [INOUT]
588  momx(:,:,:), & ! [INOUT]
589  momy(:,:,:), & ! [INOUT]
590  rhoe(:,:,:), & ! [INOUT]
591  qtrc(:,:,:,:), & ! [INOUT]
592  vterm(:,:,:,:), & ! [IN]
593  temp(:,:,:), & ! [IN]
594  mp_dtsec_sedimentation ) ! [IN]
595 
596  do j = js, je
597  do i = is, ie
598  do k = ks-1, ke
599  flx_rain(k,i,j) = flx_rain(k,i,j) + wflux_rain(k,i,j) * mp_rnstep_sedimentation
600  flx_snow(k,i,j) = flx_snow(k,i,j) + wflux_snow(k,i,j) * mp_rnstep_sedimentation
601  enddo
602  enddo
603  enddo
604 
605  enddo
606 
607  else
608  vterm(:,:,:,:) = 0.0_rp
609  endif
610 
611  call hist_in( vterm(:,:,:,i_qr), 'Vterm_QR', 'terminal velocity of QR', 'm/s' )
612  call hist_in( vterm(:,:,:,i_qi), 'Vterm_QI', 'terminal velocity of QI', 'm/s' )
613  call hist_in( vterm(:,:,:,i_qs), 'Vterm_QS', 'terminal velocity of QS', 'm/s' )
614  call hist_in( vterm(:,:,:,i_qg), 'Vterm_QG', 'terminal velocity of QG', 'm/s' )
615 
616  do j = js, je
617  do i = is, ie
618  do k = ks, ke
619  qc_before_satadj(k,i,j) = qtrc(k,i,j,i_qc)
620  enddo
621  enddo
622  enddo
623 
624  call mp_saturation_adjustment( rhoe_t(:,:,:), & ! [INOUT]
625  qtrc_t(:,:,:,:), & ! [INOUT]
626  rhoe(:,:,:), & ! [INOUT]
627  qtrc(:,:,:,:), & ! [INOUT]
628  dens(:,:,:), & ! [IN]
629  only_liquid ) ! [IN]
630 
631  do j = js, je
632  do i = is, ie
633  do k = ks, ke
634  evaporate(k,i,j) = max( qc_before_satadj(k,i,j) - qtrc(k,i,j,i_qc), 0.0_rp ) / dt ! if negative, condensation
635  evaporate(k,i,j) = evaporate(k,i,j) * dens(k,i,j) &
636  / (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  call thermodyn_rhot( rhot(:,:,:), & ! [OUT]
644  rhoe(:,:,:), & ! [IN]
645  qtrc(:,:,:,:) ) ! [IN]
646 
647  if ( mp_donegative_fixer ) then
648  call mp_negative_fixer( dens(:,:,:), & ! [INOUT]
649  rhot(:,:,:), & ! [INOUT]
650  qtrc(:,:,:,:) ) ! [INOUT]
651  endif
652 
653  sflx_rain(:,:) = flx_rain(ks-1,:,:)
654  sflx_snow(:,:) = flx_snow(ks-1,:,:)
655 
656  return
657  end subroutine atmos_phy_mp_tomita08
658 
659  !-----------------------------------------------------------------------------
661  subroutine mp_tomita08( &
662  RHOE_t, &
663  QTRC_t, &
664  RHOE0, &
665  QTRC0, &
666  CCN, &
667  DENS0 )
668  use scale_const, only: &
669  pi => const_pi, &
670  eps => const_eps, &
671  rvap => const_rvap, &
672  cl => const_cl, &
673  tem00 => const_tem00, &
674  pre00 => const_pre00, &
675  dwatr => const_dwatr
676  use scale_time, only: &
677  dt => time_dtsec_atmos_phy_mp
678  use scale_history, only: &
679  hist_query, &
680  hist_put
681  use scale_atmos_thermodyn, only: &
682  thermodyn_temp_pres_e => atmos_thermodyn_temp_pres_e, &
683  atmos_thermodyn_templhv, &
684  atmos_thermodyn_templhs, &
685  atmos_thermodyn_templhf
686  use scale_atmos_saturation, only: &
687  saturation_dens2qsat_liq => atmos_saturation_dens2qsat_liq, &
688  saturation_dens2qsat_ice => atmos_saturation_dens2qsat_ice
689  implicit none
690 
691  real(RP), intent(out) :: RHOE_t(ka,ia,ja) ! tendency rhoe [J/m3/s]
692  real(RP), intent(out) :: QTRC_t(ka,ia,ja,qa) ! tendency tracer [kg/kg/s]
693  real(RP), intent(inout) :: RHOE0 (ka,ia,ja) ! density * internal energy [J/m3]
694  real(RP), intent(inout) :: QTRC0 (ka,ia,ja,qa) ! mass concentration [kg/kg]
695  real(RP), intent(in) :: CCN (ka,ia,ja) ! CCN number concentration [#/m3]
696  real(RP), intent(in) :: DENS0 (ka,ia,ja) ! density [kg/m3]
697 
698  ! working
699  real(RP) :: rdt
700 
701  real(RP) :: TEMP0(ka,ia,ja) ! temperature [K]
702  real(RP) :: PRES0(ka,ia,ja) ! pressure [Pa]
703  real(RP) :: QSATL(ka,ia,ja) ! saturated water vapor for liquid water [kg/kg]
704  real(RP) :: QSATI(ka,ia,ja) ! saturated water vapor for ice water [kg/kg]
705 
706  real(RP) :: dens
707  real(RP) :: rhoe
708  real(RP) :: temp
709  real(RP) :: pres
710  real(RP) :: q(qa)
711  real(RP) :: Sliq ! saturated ratio S for liquid water [0-1]
712  real(RP) :: Sice ! saturated ratio S for ice water [0-1]
713 
714  real(RP) :: Rdens
715  real(RP) :: rho_fact ! density factor
716  real(RP) :: temc ! T - T0 [K]
717 
718  real(RP) :: RLMDr, RLMDr_2, RLMDr_3
719  real(RP) :: RLMDs, RLMDs_2, RLMDs_3
720  real(RP) :: RLMDg, RLMDg_2, RLMDg_3
721  real(RP) :: RLMDr_1br, RLMDr_2br, RLMDr_3br
722  real(RP) :: RLMDs_1bs, RLMDs_2bs, RLMDs_3bs
723  real(RP) :: RLMDr_dr, RLMDr_3dr, RLMDr_5dr
724  real(RP) :: RLMDs_ds, RLMDs_3ds, RLMDs_5ds
725  real(RP) :: RLMDg_dg, RLMDg_3dg, RLMDg_5dg
726  real(RP) :: RLMDr_7
727  real(RP) :: RLMDr_6dr
728 
729  !---< Roh and Satoh (2014) >---
730  real(RP) :: tems, rhoqs, Xs2
731  real(RP) :: MOMs_0, MOMs_1, MOMs_2
732  real(RP) :: MOMs_0bs, MOMs_1bs, MOMs_2bs
733  real(RP) :: MOMs_2ds, MOMs_5ds_h, RMOMs_Vt
734  real(RP) :: coef_at(4), coef_bt(4)
735  real(RP) :: loga_, b_, nm
736 
737  real(RP) :: Vti, Vtr, Vts, Vtg
738  real(RP) :: Esi_mod, Egs_mod
739  real(RP) :: rhoqc
740  real(RP) :: Nc(ka,ia,ja)
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
748  real(RP) :: Glv, Giv, Gil
749  real(RP) :: ventr, vents, ventg
750  real(RP) :: Bergeron_sw
751  real(RP) :: a1 (ka,ia,ja)
752  real(RP) :: a2 (ka,ia,ja)
753  real(RP) :: ma2(ka,ia,ja)
754  real(RP) :: dt1
755  real(RP) :: Ni50
756 
757  ! for explicit ice generation
758  real(RP), parameter :: Di_max = 500.e-6_rp
759  real(RP), parameter :: Di_a = 11.9_rp
760  real(RP) :: sw, rhoqi, XNi, XMi, Di, Ni0, Qi0, dq, qtmp
761  real(RP) :: Pidep, Pigen
762 
763  real(RP) :: ice
764 
765  real(RP) :: net, fac, fac_sw, fac_warm, fac_ice
766  real(RP) :: w(w_nmax)
767 
768  real(RP) :: tend(i_qv:i_qg)
769 
770  real(RP) :: zerosw
771  real(RP) :: tmp
772  real(RP) :: LHVEx(ka,ia,ja)
773  real(RP) :: LHFEx(ka,ia,ja)
774  real(RP) :: LHSEx(ka,ia,ja)
775 
776  integer :: k, i, j, iq
777  !---------------------------------------------------------------------------
778 
779  call prof_rapstart('MP_tomita08', 3)
780 
781 ! RHOE_t(:,:,:) = 0.0_RP
782 ! QTRC_t(:,:,:,:) = 0.0_RP
783 
784  if ( mp_couple_aerosol ) then
785  do j = js, je
786  do i = is, ie
787  do k = ks, ke
788  nc(k,i,j) = max( ccn(k,i,j)*1.e-6_rp, nc_def(i,j) ) ! [#/m3]->[#/cc]
789  enddo
790  enddo
791  enddo
792  else
793  do j = js, je
794  do i = is, ie
795  do k = ks, ke
796  nc(k,i,j) = nc_def(i,j)
797  enddo
798  enddo
799  enddo
800  endif
801 
802  rdt = 1.0_rp / dt
803 
804  call thermodyn_temp_pres_e( temp0(:,:,:), & ! [OUT]
805  pres0(:,:,:), & ! [OUT]
806  dens0(:,:,:), & ! [IN]
807  rhoe0(:,:,:), & ! [IN]
808  qtrc0(:,:,:,:) ) ! [IN]
809 
810  call saturation_dens2qsat_liq( qsatl(:,:,:), & ! [OUT]
811  temp0(:,:,:), & ! [IN]
812  dens0(:,:,:) ) ! [IN]
813 
814  call saturation_dens2qsat_ice( qsati(:,:,:), & ! [OUT]
815  temp0(:,:,:), & ! [IN]
816  dens0(:,:,:) ) ! [IN]
817 
818  call mp_tomita08_bergeronparam( temp0(:,:,:), & ! [IN]
819  a1(:,:,:), & ! [OUT]
820  a2(:,:,:), & ! [OUT]
821  ma2(:,:,:) ) ! [OUT]
822 
823  call atmos_thermodyn_templhv( lhvex, temp0 )
824  call atmos_thermodyn_templhf( lhfex, temp0 )
825  call atmos_thermodyn_templhs( lhsex, temp0 )
826 
827 !$omp parallel do &
828 !$omp private(tend, coef_bt, coef_at, q, w) &
829 !$omp private(dens, rhoe, temp, pres) &
830 !$omp private(Sliq, Sice, Rdens, rho_fact, temc) &
831 !$omp private(Bergeron_sw, zerosw) &
832 !$omp private(RLMDr, RLMDr_dr, RLMDr_2, RLMDr_3, RLMDr_7, RLMDr_1br, RLMDr_2br, RLMDr_3br, RLMDr_3dr, RLMDr_5dr, RLMDr_6dr) &
833 !$omp private(RLMDs, RLMDs_ds, RLMDs_2, RLMDs_3, RLMDs_1bs, RLMDs_2bs, RLMDs_3bs, RLMDs_3ds, RLMDs_5ds) &
834 !$omp private(MOMs_0, MOMs_1, MOMs_2, MOMs_0bs, MOMs_1bs, MOMs_2bs, MOMs_2ds, MOMs_5ds_h, RMOMs_Vt) &
835 !$omp private(rhoqc, rhoqs, Xs2, tems, loga_, b_, nm) &
836 !$omp private(RLMDg, RLMDg_dg, RLMDg_2, RLMDg_3, RLMDg_3dg, RLMDg_5dg) &
837 !$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) &
838 !$omp private(tmp, dt1, Ni50, ice, net, fac_sw, fac, fac_warm, fac_ice) &
839 !$omp collapse(3)
840 !OCL TEMP_PRIVATE(tend, coef_bt, coef_at, q, w)
841  do j = js, je
842  do i = is, ie
843  do k = ks, ke
844 
845  ! store to work
846  dens = dens0(k,i,j)
847  rhoe = rhoe0(k,i,j)
848  temp = temp0(k,i,j)
849  pres = pres0(k,i,j)
850  do iq = i_qv, i_qg
851  q(iq) = qtrc0(k,i,j,iq)
852  enddo
853 
854  ! saturation ratio S
855  sliq = q(i_qv) / max( qsatl(k,i,j), eps )
856  sice = q(i_qv) / max( qsati(k,i,j), eps )
857 
858  rdens = 1.0_rp / dens
859  rho_fact = sqrt( dens00 * rdens )
860  temc = temp - tem00
861 
862  w(i_delta1) = ( 0.5_rp + sign(0.5_rp, q(i_qr) - 1.e-4_rp ) )
863 
864  w(i_delta2) = ( 0.5_rp + sign(0.5_rp, 1.e-4_rp - q(i_qr) ) ) &
865  * ( 0.5_rp + sign(0.5_rp, 1.e-4_rp - q(i_qs) ) )
866 
867  w(i_delta3) = 0.5_rp + sign(0.5_rp, sice - 1.0_rp )
868 
869  w(i_dqv_dt) = q(i_qv) * rdt
870  w(i_dqc_dt) = q(i_qc) * rdt
871  w(i_dqr_dt) = q(i_qr) * rdt
872  w(i_dqi_dt) = q(i_qi) * rdt
873  w(i_dqs_dt) = q(i_qs) * rdt
874  w(i_dqg_dt) = q(i_qg) * rdt
875 
876  bergeron_sw = ( 0.5_rp + sign(0.5_rp, temc + 30.0_rp ) ) &
877  * ( 0.5_rp + sign(0.5_rp, 0.0_rp - temc ) ) &
878  * ( 1.0_rp - sw_useicegen )
879 
880  ! slope parameter lambda (Rain)
881  zerosw = 0.5_rp - sign(0.5_rp, q(i_qr) - 1.e-12_rp )
882  rlmdr = sqrt(sqrt( dens * q(i_qr) / ( ar * n0r * gam_1br ) + zerosw )) * ( 1.0_rp - zerosw )
883 
884  rlmdr_dr = sqrt( rlmdr ) ! **Dr
885  rlmdr_2 = rlmdr**2
886  rlmdr_3 = rlmdr**3
887  rlmdr_7 = rlmdr**7
888  rlmdr_1br = rlmdr**4 ! (1+Br)
889  rlmdr_2br = rlmdr**5 ! (2+Br)
890  rlmdr_3br = rlmdr**6 ! (3+Br)
891  rlmdr_3dr = rlmdr**3 * rlmdr_dr
892  rlmdr_5dr = rlmdr**5 * rlmdr_dr
893  rlmdr_6dr = rlmdr**6 * rlmdr_dr
894 
895  ! slope parameter lambda (Snow)
896  zerosw = 0.5_rp - sign(0.5_rp, q(i_qs) - 1.e-12_rp )
897  rlmds = sqrt(sqrt( dens * q(i_qs) / ( as * n0s * gam_1bs ) + zerosw )) * ( 1.0_rp - zerosw )
898 
899  rlmds_ds = sqrt( sqrt(rlmds) ) ! **Ds
900  rlmds_2 = rlmds**2
901  rlmds_3 = rlmds**3
902  rlmds_1bs = rlmds**4 ! (1+Bs)
903  rlmds_2bs = rlmds**5 ! (2+Bs)
904  rlmds_3bs = rlmds**6 ! (3+Bs)
905  rlmds_3ds = rlmds**3 * rlmds_ds
906  rlmds_5ds = rlmds**5 * rlmds_ds
907 
908  moms_0 = n0s * gam * rlmds ! Ns * 0th moment
909  moms_1 = n0s * gam_2 * rlmds_2 ! Ns * 1st moment
910  moms_2 = n0s * gam_3 * rlmds_3 ! Ns * 2nd moment
911  moms_0bs = n0s * gam_1bs * rlmds_1bs ! Ns * 0+bs
912  moms_1bs = n0s * gam_2bs * rlmds_2bs ! Ns * 1+bs
913  moms_2bs = n0s * gam_3bs * rlmds_3bs ! Ns * 2+bs
914  moms_2ds = n0s * gam_3ds * rlmds_3ds ! Ns * 2+ds
915  moms_5ds_h = n0s * gam_5ds_h * sqrt(rlmds_5ds) ! Ns * (5+ds)/2
916  rmoms_vt = gam_1bsds / gam_1bs * rlmds_ds
917 
918  !---< modification by Roh and Satoh (2014) >---
919  ! bimodal size distribution of snow
920  rhoqs = dens * q(i_qs)
921  zerosw = 0.5_rp - sign(0.5_rp, rhoqs - 1.e-12_rp )
922  xs2 = rhoqs / as * ( 1.0_rp - zerosw )
923 
924  tems = min( -0.1_rp, temc )
925  coef_at(1) = coef_a( 1) + tems * ( coef_a( 2) + tems * ( coef_a( 5) + tems * coef_a( 9) ) )
926  coef_at(2) = coef_a( 3) + tems * ( coef_a( 4) + tems * coef_a( 7) )
927  coef_at(3) = coef_a( 6) + tems * coef_a( 8)
928  coef_at(4) = coef_a(10)
929  coef_bt(1) = coef_b( 1) + tems * ( coef_b( 2) + tems * ( coef_b( 5) + tems * coef_b( 9) ) )
930  coef_bt(2) = coef_b( 3) + tems * ( coef_b( 4) + tems * coef_b( 7) )
931  coef_bt(3) = coef_b( 6) + tems * coef_b( 8)
932  coef_bt(4) = coef_b(10)
933  ! 0th moment
934  loga_ = coef_at(1)
935  b_ = coef_bt(1)
936  moms_0 = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
937  + ( 1.0_rp-sw_roh2014 ) * moms_0
938  ! 1st moment
939  nm = 1.0_rp
940  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
941  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
942  moms_1 = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
943  + ( 1.0_rp-sw_roh2014 ) * moms_1
944  ! 2nd moment
945  moms_2 = ( sw_roh2014 ) * xs2 &
946  + ( 1.0_rp-sw_roh2014 ) * moms_2
947  ! 0 + Bs(=2) moment
948  nm = 2.0_rp
949  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
950  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
951  moms_0bs = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
952  + ( 1.0_rp-sw_roh2014 ) * moms_0bs
953  ! 1 + Bs(=2) moment
954  nm = 3.0_rp
955  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
956  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
957  moms_1bs = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
958  + ( 1.0_rp-sw_roh2014 ) * moms_1bs
959  ! 2 + Bs(=2) moment
960  nm = 4.0_rp
961  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
962  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
963  moms_2bs = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
964  + ( 1.0_rp-sw_roh2014 ) * moms_2bs
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 = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
970  + ( 1.0_rp-sw_roh2014 ) * moms_2ds
971  ! ( 3 + Ds(=0.25) ) / 2 moment
972  nm = 1.625_rp
973  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
974  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
975  moms_5ds_h = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ &
976  + ( 1.0_rp-sw_roh2014 ) * moms_5ds_h
977  ! Bs(=2) + Ds(=0.25) moment
978  nm = 2.25_rp
979  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
980  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
981  rmoms_vt = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ / ( moms_0bs + zerosw ) &
982  + ( 1.0_rp-sw_roh2014 ) * rmoms_vt
983 
984  ! slope parameter lambda (Graupel)
985  zerosw = 0.5_rp - sign(0.5_rp, q(i_qg) - 1.e-12_rp )
986  rlmdg = sqrt(sqrt( dens * q(i_qg) / ( ag * n0g * gam_1bg ) + zerosw )) * ( 1.0_rp - zerosw )
987 
988  rlmdg_dg = sqrt( rlmdg ) ! **Dg
989  rlmdg_2 = rlmdg**2
990  rlmdg_3 = rlmdg**3
991  rlmdg_3dg = rlmdg**3 * rlmdg_dg
992  rlmdg_5dg = rlmdg**5 * rlmdg_dg
993 
994  w(i_rlmdr) = rlmdr
995  w(i_rlmds) = rlmds
996  w(i_rlmdg) = rlmdg
997 
998  !---< terminal velocity >
999  zerosw = 0.5_rp + sign(0.5_rp, q(i_qi) - 1.e-8_rp )
1000  vti = -3.29_rp * ( dens * q(i_qi) * zerosw )**0.16_rp
1001  vtr = -cr * rho_fact * gam_1brdr / gam_1br * rlmdr_dr
1002  vts = -cs * rho_fact * rmoms_vt
1003  vtg = -cg * rho_fact * gam_1bgdg / gam_1bg * rlmdg_dg
1004 
1005  !---< Accretion >---
1006  esi_mod = min( esi, esi * exp( gamma_sacr * temc ) )
1007  egs_mod = min( egs, egs * exp( gamma_gacs * temc ) )
1008 
1009  ! [Pracw] accretion rate of cloud water by rain
1010  pracw_orig = q(i_qc) * 0.25_rp * pi * erw * n0r * cr * gam_3dr * rlmdr_3dr * rho_fact
1011 
1012  pracw_kk = 67.0_rp * ( q(i_qc) * q(i_qr) )**1.15_rp ! eq.(33) in KK(2000)
1013 
1014  ! switch orig / k-k scheme
1015  w(i_pracw) = ( 1.0_rp - sw_kk2000 ) * pracw_orig &
1016  + ( sw_kk2000 ) * pracw_kk
1017 
1018  ! [Psacw] accretion rate of cloud water by snow
1019  w(i_psacw) = q(i_qc) * 0.25_rp * pi * esw * cs * moms_2ds * rho_fact
1020 
1021  ! [Pgacw] accretion rate of cloud water by graupel
1022  w(i_pgacw) = q(i_qc) * 0.25_rp * pi * egw * n0g * cg * gam_3dg * rlmdg_3dg * rho_fact
1023 
1024  ! [Praci] accretion rate of cloud ice by rain
1025  w(i_praci) = q(i_qi) * 0.25_rp * pi * eri * n0r * cr * gam_3dr * rlmdr_3dr * rho_fact
1026 
1027  ! [Psaci] accretion rate of cloud ice by snow
1028  w(i_psaci) = q(i_qi) * 0.25_rp * pi * esi_mod * cs * moms_2ds * rho_fact
1029 
1030  ! [Pgaci] accretion rate of cloud ice by grupel
1031  w(i_pgaci) = q(i_qi) * 0.25_rp * pi * egi * n0g * cg * gam_3dg * rlmdg_3dg * rho_fact * ( 1.0_rp-sw_roh2014 )
1032 
1033  ! [Piacr] accretion rate of rain by cloud ice
1034  w(i_piacr) = q(i_qi) * ar / 4.19e-13_rp * 0.25_rp * pi * eri * n0r * cr * gam_6dr * rlmdr_6dr * rho_fact
1035 
1036  ! [Psacr] accretion rate of rain by snow
1037  w(i_psacr) = ar * 0.25_rp * pi * rdens * esr * n0r * abs(vtr-vts) &
1038  * ( gam_1br * rlmdr_1br * moms_2 &
1039  + 2.0_rp * gam_2br * rlmdr_2br * moms_1 &
1040  + gam_3br * rlmdr_3br * moms_0 )
1041 
1042  ! [Pgacr] accretion rate of rain by graupel
1043  w(i_pgacr) = ar * 0.25_rp * pi * rdens * egr * n0g * n0r * abs(vtg-vtr) &
1044  * ( gam_1br * rlmdr_1br * gam_3 * rlmdg_3 &
1045  + 2.0_rp * gam_2br * rlmdr_2br * gam_2 * rlmdg_2 &
1046  + gam_3br * rlmdr_3br * gam * rlmdg )
1047 
1048  ! [Pracs] accretion rate of snow by rain
1049  w(i_pracs) = as * 0.25_rp * pi * rdens * esr * n0r * abs(vtr-vts) &
1050  * ( moms_0bs * gam_3 * rlmdr_3 &
1051  + 2.0_rp * moms_1bs * gam_2 * rlmdr_2 &
1052  + moms_2bs * gam * rlmdr )
1053 
1054  ! [Pgacs] accretion rate of snow by graupel
1055  w(i_pgacs) = as * 0.25_rp * pi * rdens * egs_mod * n0g * abs(vtg-vts) * ( 1.0_rp-sw_roh2014 ) &
1056  * ( moms_0bs * gam_3 * rlmdg_3 &
1057  + 2.0_rp * moms_1bs * gam_2 * rlmdg_2 &
1058  + moms_2bs * gam * rlmdg )
1059 
1060  !---< Auto-conversion >---
1061  ! [Praut] auto-conversion rate from cloud water to rain
1062  rhoqc = dens * q(i_qc) * 1000.0_rp ! [g/m3]
1063  dc = 0.146_rp - 5.964e-2_rp * log( nc(k,i,j) / 2000.0_rp )
1064  praut_berry = rdens * 1.67e-5_rp * rhoqc * rhoqc / ( 5.0_rp + 3.66e-2_rp * nc(k,i,j) / ( dc * rhoqc + eps ) )
1065 
1066  praut_kk = 1350.0_rp * q(i_qc)**2.47_rp * nc(k,i,j)**(-1.79_rp) ! eq.(29) in KK(2000)
1067 
1068  ! switch berry / k-k scheme
1069  w(i_praut) = ( 1.0_rp - sw_kk2000 ) * praut_berry &
1070  + ( sw_kk2000 ) * praut_kk
1071 
1072  ! [Psaut] auto-conversion rate from cloud ice to snow
1073  betai = min( beta_saut, beta_saut * exp( gamma_saut * temc ) )
1074  w(i_psaut) = max( betai*(q(i_qi)-qicrt_saut), 0.0_rp )
1075 
1076  ! [Pgaut] auto-conversion rate from snow to graupel
1077  betas = min( beta_gaut, beta_gaut * exp( gamma_gaut * temc ) )
1078  w(i_pgaut) = max( betas*(q(i_qs)-qscrt_gaut), 0.0_rp )
1079 
1080  !---< Evaporation, Sublimation >---
1081  da = ( da0 + dda_dt * temc )
1082  kd = ( dw0 + ddw_dt * temc ) * pre00 / pres
1083  nu = ( mu0 + dmu_dt * temc ) * rdens
1084 
1085  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)) )
1086  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)) )
1087  gil = 1.0_rp / lhfex(k,i,j) * (da*temc)
1088 
1089  ! [Prevp] evaporation rate of rain
1090  ventr = f1r * gam_2 * rlmdr_2 + f2r * sqrt( cr * rho_fact / nu * rlmdr_5dr ) * gam_5dr_h
1091 
1092  w(i_prevp) = 2.0_rp * pi * rdens * n0r * ( 1.0_rp-min(sliq,1.0_rp) ) * glv * ventr
1093 
1094  ! [Psdep,Pssub] deposition/sublimation rate for snow
1095  vents = f1s * moms_1 + f2s * sqrt( cs * rho_fact / nu ) * moms_5ds_h
1096 
1097  tmp = 2.0_rp * pi * rdens * ( sice-1.0_rp ) * giv * vents
1098 
1099  w(i_psdep) = ( w(i_delta3) ) * tmp ! Sice < 1
1100  w(i_pssub) = ( w(i_delta3)-1.0_rp ) * tmp ! Sice > 1
1101 
1102  ! [Psmlt] melting rate of snow
1103  w(i_psmlt) = 2.0_rp * pi * rdens * gil * vents &
1104  + cl * temc / lhfex(k,i,j) * ( w(i_psacw) + w(i_psacr) )
1105 
1106  ! [Pgdep/pgsub] deposition/sublimation rate for graupel
1107  ventg = f1g * gam_2 * rlmdg_2 + f2g * sqrt( cg * rho_fact / nu * rlmdg_5dg ) * gam_5dg_h
1108 
1109  tmp = 2.0_rp * pi * rdens * n0g * ( sice-1.0_rp ) * giv * ventg
1110 
1111  w(i_pgdep) = ( w(i_delta3) ) * tmp ! Sice < 1
1112  w(i_pgsub) = ( w(i_delta3)-1.0_rp ) * tmp ! Sice > 1
1113 
1114  ! [Pgmlt] melting rate of graupel
1115  w(i_pgmlt) = 2.0_rp * pi * rdens * n0g * gil * ventg &
1116  + cl * temc / lhfex(k,i,j) * ( w(i_pgacw) + w(i_pgacr) )
1117 
1118  ! [Pgfrz] freezing rate of graupel
1119  w(i_pgfrz) = 2.0_rp * pi * rdens * n0r * 60.0_rp * b_frz * ar * ( exp(-a_frz*temc) - 1.0_rp ) * rlmdr_7
1120 
1121  ! [Psfw,Psfi] ( Bergeron process ) growth rate of snow by Bergeron process from cloud water/ice
1122  dt1 = ( mi50**ma2(k,i,j) - mi40**ma2(k,i,j) ) / ( a1(k,i,j) * ma2(k,i,j) )
1123  ni50 = q(i_qi) * dt / ( mi50 * dt1 )
1124 
1125  w(i_psfw) = bergeron_sw * ni50 * ( a1(k,i,j) * mi50**a2(k,i,j) &
1126  + pi * eiw * dens * q(i_qc) * ri50*ri50 * vti50 )
1127  w(i_psfi) = bergeron_sw * q(i_qi) / dt1
1128 
1129  w(i_psdep) = min( w(i_psdep), w(i_dqv_dt) )
1130  w(i_pgdep) = min( w(i_pgdep), w(i_dqv_dt) )
1131 
1132  w(i_praut) = min( w(i_praut), w(i_dqc_dt) )
1133  w(i_pracw) = min( w(i_pracw), w(i_dqc_dt) )
1134  w(i_psacw) = min( w(i_psacw), w(i_dqc_dt) )
1135  w(i_pgacw) = min( w(i_pgacw), w(i_dqc_dt) )
1136  w(i_psfw ) = min( w(i_psfw ), w(i_dqc_dt) )
1137 
1138  w(i_prevp) = min( w(i_prevp), w(i_dqr_dt) )
1139  w(i_piacr) = min( w(i_piacr), w(i_dqr_dt) )
1140  w(i_psacr) = min( w(i_psacr), w(i_dqr_dt) )
1141  w(i_pgacr) = min( w(i_pgacr), w(i_dqr_dt) )
1142  w(i_pgfrz) = min( w(i_pgfrz), w(i_dqr_dt) )
1143 
1144  w(i_psaut) = min( w(i_psaut), w(i_dqi_dt) )
1145  w(i_praci) = min( w(i_praci), w(i_dqi_dt) )
1146  w(i_psaci) = min( w(i_psaci), w(i_dqi_dt) )
1147  w(i_pgaci) = min( w(i_pgaci), w(i_dqi_dt) )
1148  w(i_psfi ) = min( w(i_psfi ), w(i_dqi_dt) )
1149 
1150  w(i_pgaut) = min( w(i_pgaut), w(i_dqs_dt) )
1151  w(i_pracs) = min( w(i_pracs), w(i_dqs_dt) )
1152  w(i_pgacs) = min( w(i_pgacs), w(i_dqs_dt) )
1153  w(i_psmlt) = max( w(i_psmlt), 0.0_rp )
1154  w(i_psmlt) = min( w(i_psmlt), w(i_dqs_dt) )
1155  w(i_pssub) = min( w(i_pssub), w(i_dqs_dt) )
1156 
1157  w(i_pgmlt) = max( w(i_pgmlt), 0.0_rp )
1158  w(i_pgmlt) = min( w(i_pgmlt), w(i_dqg_dt) )
1159  w(i_pgsub) = min( w(i_pgsub), w(i_dqg_dt) )
1160 
1161  w(i_piacr_s) = ( 1.0_rp - w(i_delta1) ) * w(i_piacr)
1162  w(i_piacr_g) = ( w(i_delta1) ) * w(i_piacr)
1163  w(i_praci_s) = ( 1.0_rp - w(i_delta1) ) * w(i_praci)
1164  w(i_praci_g) = ( w(i_delta1) ) * w(i_praci)
1165  w(i_psacr_s) = ( w(i_delta2) ) * w(i_psacr)
1166  w(i_psacr_g) = ( 1.0_rp - w(i_delta2) ) * w(i_psacr)
1167  w(i_pracs ) = ( 1.0_rp - w(i_delta2) ) * w(i_pracs)
1168 
1169  ice = 0.5_rp - sign( 0.5_rp, temp0(k,i,j)-tem00 ) ! 0: warm, 1: ice
1170 
1171  ! [QC]
1172  net = &
1173  - w(i_praut) & ! [loss] c->r
1174  - w(i_pracw) & ! [loss] c->r
1175  - w(i_psacw) & ! [loss] c->s
1176  - w(i_pgacw) & ! [loss] c->g
1177  - w(i_psfw ) * ice ! [loss] c->s
1178 
1179  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1180  fac = ( fac_sw ) &
1181  + ( 1.0_rp - fac_sw ) * min( -w(i_dqc_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1182 
1183  w(i_praut) = w(i_praut) * fac
1184  w(i_pracw) = w(i_pracw) * fac
1185  w(i_psacw) = w(i_psacw) * fac
1186  w(i_pgacw) = w(i_pgacw) * fac
1187  w(i_psfw ) = w(i_psfw ) * fac
1188 
1189  ! [QI]
1190  net = &
1191  - w(i_psaut ) & ! [loss] i->s
1192  - w(i_praci_s) & ! [loss] i->s
1193  - w(i_psaci ) & ! [loss] i->s
1194  - w(i_psfi ) & ! [loss] i->s
1195  - w(i_praci_g) & ! [loss] i->g
1196  - w(i_pgaci ) ! [loss] i->g
1197 
1198  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1199  fac = ( fac_sw ) &
1200  + ( 1.0_rp - fac_sw ) * min( -w(i_dqi_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1201 
1202  w(i_praci_s) = w(i_praci_s) * fac
1203  w(i_psaut ) = w(i_psaut ) * fac
1204  w(i_psaci ) = w(i_psaci ) * fac
1205  w(i_psfi ) = w(i_psfi ) * fac
1206  w(i_praci_g) = w(i_praci_g) * fac
1207  w(i_pgaci ) = w(i_pgaci ) * fac
1208 
1209  ! [QR]
1210  net = w(i_praut ) & ! [prod] c->r
1211  + w(i_pracw ) & ! [prod] c->r
1212  - w(i_prevp ) & ! [loss] r->v
1213  + ( &
1214  + w(i_psacw ) & ! [prod] c->s
1215  + w(i_pgacw ) & ! [prod] c->g
1216  + w(i_psmlt ) & ! [prod] s->r
1217  + w(i_pgmlt ) & ! [prod] g->r
1218  ) * ( 1.0_rp - ice ) &
1219  + ( &
1220  - w(i_piacr_s) & ! [loss] r->s
1221  - w(i_psacr_s) & ! [loss] r->s
1222  - w(i_piacr_g) & ! [loss] r->g
1223  - w(i_psacr_g) & ! [loss] r->g
1224  - w(i_pgacr ) & ! [loss] r->g
1225  - w(i_pgfrz ) & ! [loss] r->g
1226  ) * ice
1227 
1228  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1229  fac = ( fac_sw ) &
1230  + ( 1.0_rp - fac_sw ) * min( -w(i_dqr_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1231 
1232  fac_warm = fac * (1.0_rp-ice) + 1.0_rp * ice
1233  fac_ice = fac * ice + 1.0_rp * (1.0_rp-ice)
1234 
1235  w(i_praut) = w(i_praut) * fac
1236  w(i_pracw) = w(i_pracw) * fac
1237  w(i_prevp) = w(i_prevp) * fac
1238 
1239  w(i_psacw) = w(i_psacw) * fac_warm
1240  w(i_pgacw) = w(i_pgacw) * fac_warm
1241  w(i_psmlt) = w(i_psmlt) * fac_warm
1242  w(i_pgmlt) = w(i_pgmlt) * fac_warm
1243 
1244  w(i_piacr_s) = w(i_piacr_s) * fac_ice
1245  w(i_psacr_s) = w(i_psacr_s) * fac_ice
1246  w(i_piacr_g) = w(i_piacr_g) * fac_ice
1247  w(i_psacr_g) = w(i_psacr_g) * fac_ice
1248  w(i_pgacr ) = w(i_pgacr ) * fac_ice
1249  w(i_pgfrz ) = w(i_pgfrz ) * fac_ice
1250 
1251  !### in the case of Sice-1 > 0, limiter is applied to qv before qs,qg
1252  ! [QV]
1253  net = w(i_prevp) & ! [prod] r->v
1254  + w(i_pssub) & ! [prod] s->v
1255  + w(i_pgsub) & ! [prod] g->v
1256  - w(i_psdep) & ! [loss] v->s
1257  - w(i_pgdep) ! [loss] v->g
1258 
1259  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1260  fac = ( fac_sw ) &
1261  + ( 1.0_rp - fac_sw ) * min( -w(i_dqv_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1262 
1263  fac = ( w(i_delta3) ) * fac &
1264  + ( 1.0_rp - w(i_delta3) )
1265 
1266  fac = fac * ice + 1.0_rp * (1.0_rp - ice)
1267 
1268  w(i_prevp) = w(i_prevp) * fac
1269  w(i_pssub) = w(i_pssub) * fac
1270  w(i_pgsub) = w(i_pgsub) * fac
1271  w(i_psdep) = w(i_psdep) * fac
1272  w(i_pgdep) = w(i_pgdep) * fac
1273 
1274  ! [QS]
1275  net = &
1276  - w(i_pgacs ) & ! [loss] s->g
1277  + ( &
1278  - w(i_psmlt ) & ! [loss] s->r
1279  ) * ( 1.0_rp - ice ) &
1280  + ( &
1281  + w(i_psdep ) & ! [prod] v->s
1282  + w(i_psacw ) & ! [prod] c->s
1283  + w(i_psfw ) & ! [prod] c->s
1284  + w(i_piacr_s) & ! [prod] r->s
1285  + w(i_psacr_s) & ! [prod] r->s
1286  + w(i_psaut ) & ! [prod] i->s
1287  + w(i_praci_s) & ! [prod] i->s
1288  + w(i_psaci ) & ! [prod] i->s
1289  + w(i_psfi ) & ! [prod] i->s
1290  - w(i_pssub ) & ! [loss] s->v
1291  - w(i_pgaut ) & ! [loss] s->g
1292  - w(i_pracs ) & ! [loss] s->g
1293  )
1294 
1295  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1296  fac = ( fac_sw ) &
1297  + ( 1.0_rp - fac_sw ) * min( -w(i_dqs_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1298 
1299  fac_warm = fac * (1.0_rp - ice) + 1.0_rp * ice
1300  fac_ice = fac * ice + 1.0_rp * (1.0_rp-ice)
1301 
1302  w(i_pgacs) = w(i_pgacs) * fac
1303 
1304  w(i_psmlt) = w(i_psmlt) * fac_warm
1305 
1306  w(i_psdep ) = w(i_psdep ) * fac_ice
1307  w(i_psacw ) = w(i_psacw ) * fac_ice
1308  w(i_psfw ) = w(i_psfw ) * fac_ice
1309  w(i_piacr_s) = w(i_piacr_s) * fac_ice
1310  w(i_psacr_s) = w(i_psacr_s) * fac_ice
1311  w(i_psaut ) = w(i_psaut ) * fac_ice
1312  w(i_praci_s) = w(i_praci_s) * fac_ice
1313  w(i_psaci ) = w(i_psaci ) * fac_ice
1314  w(i_psfi ) = w(i_psfi ) * fac_ice
1315  w(i_pssub ) = w(i_pssub ) * fac_ice
1316  w(i_pracs ) = w(i_pracs ) * fac_ice
1317  w(i_pgaut ) = w(i_pgaut ) * fac_ice
1318 
1319  ! [QG]
1320  net = w(i_pgacs ) & ! [prod] s->g
1321  + ( &
1322  - w(i_pgmlt ) & ! [loss] g->r
1323  ) * ( 1.0_rp - ice ) &
1324  + ( &
1325  + w(i_pgdep ) & ! [prod] v->g
1326  + w(i_pgacw ) & ! [prod] c->g
1327  + w(i_piacr_g) & ! [prod] r->g
1328  + w(i_psacr_g) & ! [prod] r->g
1329  + w(i_pgacr ) & ! [prod] r->g
1330  + w(i_pgfrz ) & ! [prod] r->g
1331  + w(i_praci_g) & ! [prod] i->g
1332  + w(i_pgaci ) & ! [prod] i->g
1333  + w(i_pgaut ) & ! [prod] s->g
1334  + w(i_pracs ) & ! [prod] s->g
1335  - w(i_pgsub ) & ! [loss] g->v
1336  ) * ice
1337 
1338  fac_sw = 0.5_rp + sign( 0.5_rp, net+eps ) ! if production > loss , fac_sw=1
1339  fac = ( fac_sw ) &
1340  + ( 1.0_rp - fac_sw ) * min( -w(i_dqg_dt)/(net-fac_sw), 1.0_rp ) ! loss limiter
1341 
1342  fac_warm = fac * (1.0_rp - ice) + 1.0_rp * ice
1343  fac_ice = fac * ice + 1.0_rp * (1.0_rp-ice)
1344 
1345  w(i_pgacs) = w(i_pgacs) * fac
1346 
1347  w(i_pgmlt) = w(i_pgmlt) * fac_warm
1348 
1349  w(i_pgdep ) = w(i_pgdep ) * fac_ice
1350  w(i_pgacw ) = w(i_pgacw ) * fac_ice
1351  w(i_piacr_g) = w(i_piacr_g) * fac_ice
1352  w(i_psacr_g) = w(i_psacr_g) * fac_ice
1353  w(i_pgacr ) = w(i_pgacr ) * fac_ice
1354  w(i_pgfrz ) = w(i_pgfrz ) * fac_ice
1355  w(i_praci_g) = w(i_praci_g) * fac_ice
1356  w(i_pgaci ) = w(i_pgaci ) * fac_ice
1357  w(i_pgaut ) = w(i_pgaut ) * fac_ice
1358  w(i_pracs ) = w(i_pracs ) * fac_ice
1359  w(i_pgsub ) = w(i_pgsub ) * fac_ice
1360 
1361 
1362 
1363  ! re-calc net production & loss
1364  tend(i_qc) = &
1365  - w(i_praut ) & ! [loss] c->r
1366  - w(i_pracw ) & ! [loss] c->r
1367  - w(i_psacw ) & ! [loss] c->s
1368  - w(i_pgacw ) & ! [loss] c->g
1369  - w(i_psfw ) * ice
1370 
1371  tend(i_qr) = w(i_praut ) & ! [prod] c->r
1372  + w(i_pracw ) & ! [prod] c->r
1373  - w(i_prevp ) & ! [loss] r->v
1374  + ( &
1375  + w(i_psacw ) & ! [prod] c->s=r
1376  + w(i_pgacw ) & ! [prod] c->g=r
1377  + w(i_psmlt ) & ! [prod] s->r
1378  + w(i_pgmlt ) & ! [prod] g->r
1379  ) * ( 1.0_rp - ice ) &
1380  + ( &
1381  - w(i_piacr_s) & ! [loss] r->s
1382  - w(i_psacr_s) & ! [loss] r->s
1383  - w(i_piacr_g) & ! [loss] r->g
1384  - w(i_psacr_g) & ! [loss] r->g
1385  - w(i_pgacr ) & ! [loss] r->g
1386  - w(i_pgfrz ) & ! [loss] r->g
1387  ) * ice
1388 
1389  tend(i_qi) = ( &
1390  - w(i_psaut ) & ! [loss] i->s
1391  - w(i_praci_s) & ! [loss] i->s
1392  - w(i_psaci ) & ! [loss] i->s
1393  - w(i_psfi ) & ! [loss] i->s
1394  - w(i_praci_g) & ! [loss] i->g
1395  - w(i_pgaci ) & ! [loss] i->g
1396  ) * ice
1397 
1398  tend(i_qs) = &
1399  - w(i_pgacs ) & ! [loss] s->g
1400  + ( &
1401  - w(i_psmlt ) & ! [loss] s->r
1402  ) * ( 1.0_rp - ice ) &
1403  + ( &
1404  + w(i_psdep ) & ! [prod] v->s
1405  + w(i_psacw ) & ! [prod] c->s
1406  + w(i_psfw ) & ! [prod] c->s
1407  + w(i_piacr_s) & ! [prod] r->s
1408  + w(i_psacr_s) & ! [prod] r->s
1409  + w(i_psaut ) & ! [prod] i->s
1410  + w(i_psaci ) & ! [prod] i->s
1411  + w(i_praci_s) & ! [prod] i->s
1412  + w(i_psfi ) & ! [prod] i->s
1413  - w(i_pssub ) & ! [loss] s->v
1414  - w(i_pgaut ) & ! [loss] s->g
1415  - w(i_pracs ) & ! [loss] s->g
1416  ) * ice
1417 
1418  tend(i_qg) = w(i_pgacs ) & ! [prod] s->g
1419  + ( &
1420  - w(i_pgmlt ) & ! [loss] g->r
1421  ) * ( 1.0_rp - ice ) &
1422  + ( &
1423  + w(i_pgdep ) & ! [prod] v->g
1424  + w(i_pgacw ) & ! [prod] c->g
1425  + w(i_piacr_g) & ! [prod] r->g
1426  + w(i_psacr_g) & ! [prod] r->g
1427  + w(i_pgacr ) & ! [prod] r->g
1428  + w(i_pgfrz ) & ! [prod] r->g
1429  + w(i_praci_g) & ! [prod] i->g
1430  + w(i_pgaci ) & ! [prod] i->g
1431  + w(i_pgaut ) & ! [prod] s->g
1432  + w(i_pracs ) & ! [prod] s->g
1433  - w(i_pgsub ) & ! [loss] g->v
1434  ) * ice
1435 
1436  tend(i_qc) = max( tend(i_qc), -w(i_dqc_dt) )
1437  tend(i_qr) = max( tend(i_qr), -w(i_dqr_dt) )
1438  tend(i_qi) = max( tend(i_qi), -w(i_dqi_dt) )
1439  tend(i_qs) = max( tend(i_qs), -w(i_dqs_dt) )
1440  tend(i_qg) = max( tend(i_qg), -w(i_dqg_dt) )
1441 
1442  tend(i_qv) = - ( tend(i_qc) &
1443  + tend(i_qr) &
1444  + tend(i_qi) &
1445  + tend(i_qs) &
1446  + tend(i_qg) )
1447 
1448  qtrc_t(k,i,j,i_qv) = tend(i_qv)
1449  qtrc_t(k,i,j,i_qc) = tend(i_qc)
1450  qtrc_t(k,i,j,i_qr) = tend(i_qr)
1451  qtrc_t(k,i,j,i_qi) = tend(i_qi)
1452  qtrc_t(k,i,j,i_qs) = tend(i_qs)
1453  qtrc_t(k,i,j,i_qg) = tend(i_qg)
1454 
1455  qtrc0(k,i,j,i_qv) = qtrc0(k,i,j,i_qv) + qtrc_t(k,i,j,i_qv) * dt
1456  qtrc0(k,i,j,i_qc) = qtrc0(k,i,j,i_qc) + qtrc_t(k,i,j,i_qc) * dt
1457  qtrc0(k,i,j,i_qr) = qtrc0(k,i,j,i_qr) + qtrc_t(k,i,j,i_qr) * dt
1458  qtrc0(k,i,j,i_qi) = qtrc0(k,i,j,i_qi) + qtrc_t(k,i,j,i_qi) * dt
1459  qtrc0(k,i,j,i_qs) = qtrc0(k,i,j,i_qs) + qtrc_t(k,i,j,i_qs) * dt
1460  qtrc0(k,i,j,i_qg) = qtrc0(k,i,j,i_qg) + qtrc_t(k,i,j,i_qg) * dt
1461  enddo
1462  enddo
1463  enddo
1464 
1465  !-----< ice generation process >-----
1466  if ( mp_doexpricit_icegen ) then
1467  do j = js, je
1468  do i = is, ie
1469  do k = ks, ke
1470  ! store to work
1471  dens = dens0(k,i,j)
1472  temp = temp0(k,i,j)
1473  pres = pres0(k,i,j)
1474 
1475  ! saturation ratio S
1476  sice = qtrc0(k,i,j,i_qv) / max( qsati(k,i,j), eps )
1477 
1478  rdens = 1.0_rp / dens
1479  temc = temp - tem00
1480  rhoqi = dens * qtrc0(k,i,j,i_qi)
1481 
1482  da = ( da0 + dda_dt * temc )
1483  kd = ( dw0 + ddw_dt * temc ) * pre00 / pres
1484 
1485  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)) )
1486 
1487  ! [Pidep] deposition/sublimation : v->i or i->v
1488  sw = ( 0.5_rp + sign(0.5_rp, 0.0_rp - temc ) ) &
1489  * ( 0.5_rp + sign(0.5_rp, rhoqi ) ) ! if T < 0C & ice exists, sw=1
1490 
1491  rhoqi = max(rhoqi,eps)
1492 
1493  xni = min( max( 5.38e7_rp * rhoqi**0.75_rp, 1.e3_rp ), 1.e6_rp )
1494  xmi = rhoqi / xni
1495  di = min( di_a * sqrt(xmi), di_max )
1496 
1497  pidep = sw * 4.0_rp * di * xni * rdens * ( sice-1.0_rp ) * giv
1498 
1499  ! [Pigen] ice nucleation : v->i
1500  sw = ( 0.5_rp + sign(0.5_rp, 0.0_rp - temc ) ) &
1501  * ( 0.5_rp + sign(0.5_rp, sice - 1.0_rp ) ) ! if T < 0C & Satulated, sw=1
1502 
1503  ni0 = max( exp(-0.1_rp*temc), 1.0_rp ) * 1000.0_rp
1504  qi0 = 4.92e-11 * ni0**1.33_rp * rdens
1505 
1506  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
1507 
1508  ! update Qv, Qi with limiter
1509  dq = ( pigen + pidep ) * dt
1510  qtmp = qtrc0(k,i,j,i_qv) - dq
1511  if ( dq > 0.0_rp ) then ! v->i
1512  qtmp = max( qtmp, qsati(k,i,j) )
1513  else ! i->v
1514  qtmp = min( qtmp, qsati(k,i,j) )
1515  endif
1516  dq = qtrc0(k,i,j,i_qv) - qtmp
1517  qtmp = qtrc0(k,i,j,i_qi) + dq
1518  qtmp = max( qtmp, 0.0_rp )
1519  dq = qtmp - qtrc0(k,i,j,i_qi)
1520  qtrc0(k,i,j,i_qi) = qtrc0(k,i,j,i_qi) + dq
1521  qtrc0(k,i,j,i_qv) = qtrc0(k,i,j,i_qv) - dq
1522  qtrc_t(k,i,j,i_qi) = qtrc_t(k,i,j,i_qi) + dq /dt
1523  qtrc_t(k,i,j,i_qv) = qtrc_t(k,i,j,i_qv) - dq /dt
1524 
1525  ! [Pifrz/Pimlt] freezing and melting
1526 
1527  ! homogenious freezing : c->i
1528  sw = ( 0.5_rp - sign(0.5_rp, temc + 40.0_rp ) ) ! if T < -40C, sw=1
1529 
1530  dq = sw * qtrc0(k,i,j,i_qc)
1531  qtrc0(k,i,j,i_qi) = qtrc0(k,i,j,i_qi) + dq
1532  qtrc0(k,i,j,i_qc) = qtrc0(k,i,j,i_qc) - dq
1533  qtrc_t(k,i,j,i_qi) = qtrc_t(k,i,j,i_qi) + dq /dt
1534  qtrc_t(k,i,j,i_qc) = qtrc_t(k,i,j,i_qc) - dq /dt
1535 
1536  ! heteroginous freezing : c->i
1537  sw = ( 0.5_rp + sign(0.5_rp, temc + 40.0_rp ) ) &
1538  * ( 0.5_rp + sign(0.5_rp, 0.0_rp - temc ) ) ! if -40C < T < 0C, sw=1
1539 
1540  dq = sw * ( dens / dwatr * qtrc0(k,i,j,i_qc)**2 / (nc(k,i,j)*1.e+6) ) &
1541  * b_frz * ( exp(-a_frz*temc) - 1.0_rp ) * dt
1542  qtmp = qtrc0(k,i,j,i_qc) - dq
1543  qtmp = max( qtmp, 0.0_rp )
1544  dq = qtrc0(k,i,j,i_qc) - qtmp
1545  qtrc0(k,i,j,i_qi) = qtrc0(k,i,j,i_qi) + dq
1546  qtrc0(k,i,j,i_qc) = qtrc0(k,i,j,i_qc) - dq
1547  qtrc_t(k,i,j,i_qi) = qtrc_t(k,i,j,i_qi) + dq /dt
1548  qtrc_t(k,i,j,i_qc) = qtrc_t(k,i,j,i_qc) - dq /dt
1549 
1550  ! melting : i->c
1551  sw = ( 0.5_rp - sign(0.5_rp, 0.0_rp - temc ) ) ! if 0C < T, sw=1
1552 
1553  dq = - sw * qtrc0(k,i,j,i_qi)
1554  qtrc0(k,i,j,i_qi) = qtrc0(k,i,j,i_qi) + dq
1555  qtrc0(k,i,j,i_qc) = qtrc0(k,i,j,i_qc) - dq
1556  qtrc_t(k,i,j,i_qi) = qtrc_t(k,i,j,i_qi) + dq /dt
1557  qtrc_t(k,i,j,i_qc) = qtrc_t(k,i,j,i_qc) - dq /dt
1558 
1559  enddo
1560  enddo
1561  enddo
1562  endif
1563 
1564  do j = js, je
1565  do i = is, ie
1566  do k = ks, ke
1567  rhoe_t(k,i,j) = - dens0(k,i,j) * ( lhvex(k,i,j) * qtrc_t(k,i,j,i_qv) &
1568  - lhfex(k,i,j) * qtrc_t(k,i,j,i_qi) &
1569  - lhfex(k,i,j) * qtrc_t(k,i,j,i_qs) &
1570  - lhfex(k,i,j) * qtrc_t(k,i,j,i_qg) )
1571 
1572  rhoe0(k,i,j) = rhoe0(k,i,j) + rhoe_t(k,i,j) * dt
1573  enddo
1574  enddo
1575  enddo
1576 
1577 ! do ip = 1, w_nmax
1578 ! if ( w_histid(ip) > 0 ) then
1579 
1580 ! call HIST_query( w_histid(ip), & ! [IN]
1581 ! do_put ) ! [OUT]
1582 
1583 ! if ( do_put ) then
1584 ! do j = JS, JE
1585 ! do i = IS, IE
1586 ! do k = KS, KE
1587 ! work3D(k,i,j) = w(ip)
1588 ! enddo
1589 ! enddo
1590 ! enddo
1591 
1592 ! if( IO_L ) write(IO_FID_LOG,*) w_name(ip), "MAX/MIN:", &
1593 ! maxval(work3D(KS:KE,IS:IE,JS:JE)), minval(work3D(KS:KE,IS:IE,JS:JE))
1594 
1595 
1596 ! call HIST_put( w_histid(ip), work3D(:,:,:), w_zinterp(ip) )
1597 ! endif
1598 ! endif
1599 ! enddo
1600 
1601 ! if( IO_L ) write(IO_FID_LOG,*) "tend QV MAX/MIN:", maxval(tend(I_QV,:)), minval(tend(I_QV,:))
1602 ! if( IO_L ) write(IO_FID_LOG,*) "tend QC MAX/MIN:", maxval(tend(I_QC,:)), minval(tend(I_QC,:))
1603 ! if( IO_L ) write(IO_FID_LOG,*) "tend QR MAX/MIN:", maxval(tend(I_QR,:)), minval(tend(I_QR,:))
1604 ! if( IO_L ) write(IO_FID_LOG,*) "tend QI MAX/MIN:", maxval(tend(I_QI,:)), minval(tend(I_QI,:))
1605 ! if( IO_L ) write(IO_FID_LOG,*) "tend QS MAX/MIN:", maxval(tend(I_QS,:)), minval(tend(I_QS,:))
1606 ! if( IO_L ) write(IO_FID_LOG,*) "tend QG MAX/MIN:", maxval(tend(I_QG,:)), minval(tend(I_QG,:))
1607 !
1608 ! if( IO_L ) write(IO_FID_LOG,*) "QV MAX/MIN:", maxval(QTRC0(:,:,:,I_QV)), minval(QTRC0(:,:,:,I_QV))
1609 ! if( IO_L ) write(IO_FID_LOG,*) "QC MAX/MIN:", maxval(QTRC0(:,:,:,I_QC)), minval(QTRC0(:,:,:,I_QC))
1610 ! if( IO_L ) write(IO_FID_LOG,*) "QR MAX/MIN:", maxval(QTRC0(:,:,:,I_QR)), minval(QTRC0(:,:,:,I_QR))
1611 ! if( IO_L ) write(IO_FID_LOG,*) "QI MAX/MIN:", maxval(QTRC0(:,:,:,I_QI)), minval(QTRC0(:,:,:,I_QI))
1612 ! if( IO_L ) write(IO_FID_LOG,*) "QS MAX/MIN:", maxval(QTRC0(:,:,:,I_QS)), minval(QTRC0(:,:,:,I_QS))
1613 ! if( IO_L ) write(IO_FID_LOG,*) "QG MAX/MIN:", maxval(QTRC0(:,:,:,I_QG)), minval(QTRC0(:,:,:,I_QG))
1614 
1615 ! if( IO_L ) write(IO_FID_LOG,*) "QV_t MAX/MIN:", maxval(QTRC_t(:,:,:,I_QV)), minval(QTRC_t(:,:,:,I_QV))
1616 ! if( IO_L ) write(IO_FID_LOG,*) "QC_t MAX/MIN:", maxval(QTRC_t(:,:,:,I_QC)), minval(QTRC_t(:,:,:,I_QC))
1617 ! if( IO_L ) write(IO_FID_LOG,*) "QR_t MAX/MIN:", maxval(QTRC_t(:,:,:,I_QR)), minval(QTRC_t(:,:,:,I_QR))
1618 ! if( IO_L ) write(IO_FID_LOG,*) "QI_t MAX/MIN:", maxval(QTRC_t(:,:,:,I_QI)), minval(QTRC_t(:,:,:,I_QI))
1619 ! if( IO_L ) write(IO_FID_LOG,*) "QS_t MAX/MIN:", maxval(QTRC_t(:,:,:,I_QS)), minval(QTRC_t(:,:,:,I_QS))
1620 ! if( IO_L ) write(IO_FID_LOG,*) "QG_t MAX/MIN:", maxval(QTRC_t(:,:,:,I_QG)), minval(QTRC_t(:,:,:,I_QG))
1621 
1622 ! if( IO_L ) write(IO_FID_LOG,*) "RHOE0 MAX/MIN:", maxval(RHOE0(:,:,:)), minval(RHOE0(:,:,:))
1623 ! if( IO_L ) write(IO_FID_LOG,*) "RHOE_t MAX/MIN:", maxval(RHOE_t(:,:,:)), minval(RHOE_t(:,:,:))
1624 
1625  call prof_rapend ('MP_tomita08', 3)
1626 
1627  return
1628  end subroutine mp_tomita08
1629 
1630  !-----------------------------------------------------------------------------
1632  subroutine mp_tomita08_vterm( &
1633  vterm, &
1634  DENS0, &
1635  TEMP0, &
1636  QTRC0 )
1637  use scale_const, only: &
1638  tem00 => const_tem00
1639  implicit none
1640 
1641  real(RP), intent(out) :: vterm(ka,ia,ja,qa)
1642  real(RP), intent(in) :: DENS0(ka,ia,ja)
1643  real(RP), intent(in) :: TEMP0(ka,ia,ja)
1644  real(RP), intent(in) :: QTRC0(ka,ia,ja,qa)
1645 
1646  real(RP) :: dens
1647  real(RP) :: temc
1648  real(RP) :: q(qa)
1649 
1650  real(RP) :: rho_fact ! density factor
1651 
1652  real(RP) :: RLMDr, RLMDs, RLMDg
1653  real(RP) :: RLMDr_dr, RLMDs_ds, RLMDg_dg
1654 
1655  !---< Roh and Satoh (2014) >---
1656  real(RP) :: tems, rhoqs, Xs2
1657  real(RP) :: MOMs_0bs, RMOMs_Vt
1658  real(RP) :: coef_at(4), coef_bt(4)
1659  real(RP) :: loga_, b_, nm
1660 
1661  real(RP) :: zerosw
1662  integer :: k, i, j, iq
1663  !---------------------------------------------------------------------------
1664 
1665  do j = js, je
1666  do i = is, ie
1667  do k = ks, ke
1668  ! store to work
1669  dens = dens0(k,i,j)
1670  temc = temp0(k,i,j) - tem00
1671  do iq = i_qv, i_qg
1672  q(iq) = qtrc0(k,i,j,iq)
1673  enddo
1674 
1675  rho_fact = sqrt( dens00 / dens )
1676 
1677  ! slope parameter lambda (Rain)
1678  zerosw = 0.5_rp - sign(0.5_rp, q(i_qr) - 1.e-12_rp )
1679  rlmdr = sqrt(sqrt( dens * q(i_qr) / ( ar * n0r * gam_1br ) + zerosw )) * ( 1.0_rp - zerosw )
1680 
1681  rlmdr_dr = sqrt( rlmdr ) ! **Dr
1682 
1683  ! slope parameter lambda (Snow)
1684  zerosw = 0.5_rp - sign(0.5_rp, q(i_qs) - 1.e-12_rp )
1685  rlmds = sqrt(sqrt( dens * q(i_qs) / ( as * n0s * gam_1bs ) + zerosw )) * ( 1.0_rp - zerosw )
1686 
1687  rlmds_ds = sqrt( sqrt(rlmds) ) ! **Ds
1688  rmoms_vt = gam_1bsds / gam_1bs * rlmds_ds
1689 
1690  !---< modification by Roh and Satoh (2014) >---
1691  ! bimodal size distribution of snow
1692  rhoqs = dens * q(i_qs)
1693  zerosw = 0.5_rp - sign(0.5_rp, rhoqs - 1.e-12_rp )
1694  xs2 = rhoqs / as * ( 1.0_rp - zerosw )
1695 
1696  tems = min( -0.1_rp, temc )
1697  coef_at(1) = coef_a( 1) + tems * ( coef_a( 2) + tems * ( coef_a( 5) + tems * coef_a( 9) ) )
1698  coef_at(2) = coef_a( 3) + tems * ( coef_a( 4) + tems * coef_a( 7) )
1699  coef_at(3) = coef_a( 6) + tems * coef_a( 8)
1700  coef_at(4) = coef_a(10)
1701  coef_bt(1) = coef_b( 1) + tems * ( coef_b( 2) + tems * ( coef_b( 5) + tems * coef_b( 9) ) )
1702  coef_bt(2) = coef_b( 3) + tems * ( coef_b( 4) + tems * coef_b( 7) )
1703  coef_bt(3) = coef_b( 6) + tems * coef_b( 8)
1704  coef_bt(4) = coef_b(10)
1705  ! 0 + Bs(=2) moment
1706  nm = 2.0_rp
1707  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
1708  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
1709  moms_0bs = 10.0_rp**loga_ * xs2**b_
1710  ! Bs(=2) + Ds(=0.25) moment
1711  nm = 2.25_rp
1712  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
1713  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
1714  rmoms_vt = ( sw_roh2014 ) * 10.0_rp**loga_ * xs2**b_ / ( moms_0bs + zerosw ) &
1715  + ( 1.0_rp-sw_roh2014 ) * rmoms_vt
1716 
1717  ! slope parameter lambda (Graupel)
1718  zerosw = 0.5_rp - sign(0.5_rp, q(i_qg) - 1.e-12_rp )
1719  rlmdg = sqrt(sqrt( dens * q(i_qg) / ( ag * n0g * gam_1bg ) + zerosw )) * ( 1.0_rp - zerosw )
1720 
1721  rlmdg_dg = sqrt( rlmdg ) ! **Dg
1722 
1723  !---< terminal velocity >
1724  zerosw = 0.5_rp + sign(0.5_rp, q(i_qi) - 1.e-8_rp )
1725  vterm(k,i,j,i_qv) = 0.0_rp
1726  vterm(k,i,j,i_qc) = 0.0_rp
1727  vterm(k,i,j,i_qi) = -3.29_rp * ( dens * q(i_qi) * zerosw )**0.16_rp
1728  vterm(k,i,j,i_qr) = -cr * rho_fact * gam_1brdr / gam_1br * rlmdr_dr
1729  vterm(k,i,j,i_qs) = -cs * rho_fact * rmoms_vt
1730  vterm(k,i,j,i_qg) = -cg * rho_fact * gam_1bgdg / gam_1bg * rlmdg_dg
1731  enddo
1732  enddo
1733  enddo
1734 
1735  do j = js, je
1736  do i = is, ie
1737  vterm( 1:ks-1,i,j,:) = 0.0_rp
1738  vterm(ke+1:ka ,i,j,:) = 0.0_rp
1739  enddo
1740  enddo
1741 
1742  return
1743  end subroutine mp_tomita08_vterm
1744 
1745  !-----------------------------------------------------------------------------
1746  subroutine mp_tomita08_bergeronparam( &
1747  temp, &
1748  a1, &
1749  a2, &
1750  ma2 )
1751  use scale_const, only: &
1752  tem00 => const_tem00
1753  implicit none
1754 
1755  real(RP), intent(in) :: temp(ka,ia,ja)
1756  real(RP), intent(out) :: a1 (ka,ia,ja)
1757  real(RP), intent(out) :: a2 (ka,ia,ja)
1758  real(RP), intent(out) :: ma2 (ka,ia,ja)
1759 
1760  real(RP) :: a1_tab(32)
1761  real(RP) :: a2_tab(32)
1762 
1763  data a1_tab / 0.0001e-7_rp, 0.7939e-7_rp, 0.7841e-6_rp, 0.3369e-5_rp, 0.4336e-5_rp, &
1764  0.5285e-5_rp, 0.3728e-5_rp, 0.1852e-5_rp, 0.2991e-6_rp, 0.4248e-6_rp, &
1765  0.7434e-6_rp, 0.1812e-5_rp, 0.4394e-5_rp, 0.9145e-5_rp, 0.1725e-4_rp, &
1766  0.3348e-4_rp, 0.1725e-4_rp, 0.9175e-5_rp, 0.4412e-5_rp, 0.2252e-5_rp, &
1767  0.9115e-6_rp, 0.4876e-6_rp, 0.3473e-6_rp, 0.4758e-6_rp, 0.6306e-6_rp, &
1768  0.8573e-6_rp, 0.7868e-6_rp, 0.7192e-6_rp, 0.6513e-6_rp, 0.5956e-6_rp, &
1769  0.5333e-6_rp, 0.4834e-6_rp /
1770 
1771  data a2_tab / 0.0100_rp, 0.4006_rp, 0.4831_rp, 0.5320_rp, 0.5307_rp, &
1772  0.5319_rp, 0.5249_rp, 0.4888_rp, 0.3849_rp, 0.4047_rp, &
1773  0.4318_rp, 0.4771_rp, 0.5183_rp, 0.5463_rp, 0.5651_rp, &
1774  0.5813_rp, 0.5655_rp, 0.5478_rp, 0.5203_rp, 0.4906_rp, &
1775  0.4447_rp, 0.4126_rp, 0.3960_rp, 0.4149_rp, 0.4320_rp, &
1776  0.4506_rp, 0.4483_rp, 0.4460_rp, 0.4433_rp, 0.4413_rp, &
1777  0.4382_rp, 0.4361_rp /
1778 
1779  real(RP) :: temc
1780  integer :: itemc
1781  real(RP) :: fact
1782 
1783  integer :: k, i, j
1784  !---------------------------------------------------------------------------
1785 
1786  do j = js, je
1787  do i = is, ie
1788  do k = ks, ke
1789  temc = min( max( temp(k,i,j)-tem00, -30.99_rp ), 0.0_rp ) ! 0C <= T < 31C
1790  itemc = int( -temc ) + 1 ! 1 <= iT <= 31
1791  fact = - ( temc + real(itemc-1,kind=8) )
1792 
1793  a1(k,i,j) = ( 1.0_rp-fact ) * a1_tab(itemc ) &
1794  + ( fact ) * a1_tab(itemc+1)
1795 
1796  a2(k,i,j) = ( 1.0_rp-fact ) * a2_tab(itemc ) &
1797  + ( fact ) * a2_tab(itemc+1)
1798 
1799  ma2(k,i,j) = 1.0_rp - a2(k,i,j)
1800 
1801  a1(k,i,j) = a1(k,i,j) * 1.e-3_rp**ma2(k,i,j) ! [g->kg]
1802  enddo
1803  enddo
1804  enddo
1805 
1806  return
1807  end subroutine mp_tomita08_bergeronparam
1808 
1809  !-----------------------------------------------------------------------------
1811  subroutine atmos_phy_mp_tomita08_cloudfraction( &
1812  cldfrac, &
1813  QTRC, &
1814  mask_criterion )
1815  use scale_grid_index
1816  use scale_tracer, only: &
1817  qad => qa
1818  implicit none
1819 
1820  real(RP), intent(out) :: cldfrac(ka,ia,ja)
1821  real(RP), intent(in) :: QTRC (ka,ia,ja,qad)
1822  real(RP), intent(in) :: mask_criterion
1823 
1824  real(RP) :: qhydro
1825  integer :: k, i, j, iq
1826  !---------------------------------------------------------------------------
1827 
1828  do j = js, je
1829  do i = is, ie
1830  do k = ks, ke
1831  qhydro = 0.0_rp
1832  do iq = 1, mp_qa
1833  qhydro = qhydro + qtrc(k,i,j,i_mp2all(iq))
1834  enddo
1835  cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
1836  enddo
1837  enddo
1838  enddo
1839 
1840  return
1841  end subroutine atmos_phy_mp_tomita08_cloudfraction
1842 
1843  !-----------------------------------------------------------------------------
1845  subroutine atmos_phy_mp_tomita08_effectiveradius( &
1846  Re, &
1847  QTRC0, &
1848  DENS0, &
1849  TEMP0 )
1850  use scale_grid_index
1851  use scale_const, only: &
1852  tem00 => const_tem00
1853  use scale_tracer, only: &
1854  qad => qa, &
1855  mp_qad => mp_qa
1856  implicit none
1857 
1858  real(RP), intent(out) :: Re (ka,ia,ja,mp_qad) ! effective radius [cm]
1859  real(RP), intent(in) :: QTRC0(ka,ia,ja,qad) ! tracer mass concentration [kg/kg]
1860  real(RP), intent(in) :: DENS0(ka,ia,ja) ! density [kg/m3]
1861  real(RP), intent(in) :: TEMP0(ka,ia,ja) ! temperature [K]
1862 
1863  real(RP) :: dens
1864  real(RP) :: temc
1865  real(RP) :: q(qa)
1866  real(RP) :: RLMDr, RLMDs, RLMDg
1867 
1868  real(RP), parameter :: um2cm = 100.0_rp
1869 
1870  !---< Roh and Satoh (2014) >---
1871  real(RP) :: tems, rhoqs, Xs2
1872  real(RP) :: MOMs_2, MOMs_1bs
1873  real(RP) :: coef_at(4), coef_bt(4)
1874  real(RP) :: loga_, b_, nm
1875 
1876  real(RP) :: zerosw
1877  integer :: k, i, j, iq
1878  !---------------------------------------------------------------------------
1879 
1880  re(:,:,:,i_mp_qc) = re_qc * um2cm
1881  re(:,:,:,i_mp_qi) = re_qi * um2cm
1882 
1883  do j = js, je
1884  do i = is, ie
1885  do k = ks, ke
1886  dens = dens0(k,i,j)
1887  temc = temp0(k,i,j) - tem00
1888  do iq = i_qv, i_qg
1889  q(iq) = qtrc0(k,i,j,iq)
1890  enddo
1891 
1892  ! slope parameter lambda
1893  zerosw = 0.5_rp - sign(0.5_rp, q(i_qr) - 1.e-12_rp )
1894  rlmdr = sqrt(sqrt( dens * q(i_qr) / ( ar * n0r * gam_1br ) + zerosw )) * ( 1.0_rp - zerosw )
1895 
1896  ! Effective radius is defined by r3m/r2m = 1.5/lambda
1897  re(k,i,j,i_mp_qr) = 1.5_rp * rlmdr * um2cm
1898 
1899  zerosw = 0.5_rp - sign(0.5_rp, q(i_qs) - 1.e-12_rp )
1900  rlmds = sqrt(sqrt( dens * q(i_qs) / ( as * n0s * gam_1bs ) + zerosw )) * ( 1.0_rp - zerosw )
1901 
1902  !---< modification by Roh and Satoh (2014) >---
1903  ! bimodal size distribution of snow
1904  rhoqs = dens * q(i_qs)
1905  zerosw = 0.5_rp - sign(0.5_rp, rhoqs - 1.e-12_rp )
1906  xs2 = rhoqs / as * ( 1.0_rp - zerosw )
1907 
1908  tems = min( -0.1_rp, temc )
1909  coef_at(1) = coef_a( 1) + tems * ( coef_a( 2) + tems * ( coef_a( 5) + tems * coef_a( 9) ) )
1910  coef_at(2) = coef_a( 3) + tems * ( coef_a( 4) + tems * coef_a( 7) )
1911  coef_at(3) = coef_a( 6) + tems * coef_a( 8)
1912  coef_at(4) = coef_a(10)
1913  coef_bt(1) = coef_b( 1) + tems * ( coef_b( 2) + tems * ( coef_b( 5) + tems * coef_b( 9) ) )
1914  coef_bt(2) = coef_b( 3) + tems * ( coef_b( 4) + tems * coef_b( 7) )
1915  coef_bt(3) = coef_b( 6) + tems * coef_b( 8)
1916  coef_bt(4) = coef_b(10)
1917  ! 2nd moment
1918  moms_2 = xs2
1919  ! 1 + Bs(=2) moment
1920  nm = 3.0_rp
1921  loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
1922  b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
1923  moms_1bs = 10.0_rp**loga_ * xs2**b_
1924 
1925  re(k,i,j,i_mp_qs) = ( sw_roh2014 ) * 0.5_rp * moms_1bs / ( moms_2 + zerosw ) * um2cm &
1926  + ( 1.0_rp-sw_roh2014 ) * 1.5_rp * rlmds * um2cm
1927 
1928  zerosw = 0.5_rp - sign(0.5_rp, q(i_qg) - 1.e-12_rp )
1929  rlmdg = sqrt(sqrt( dens * q(i_qg) / ( ag * n0g * gam_1bg ) + zerosw )) * ( 1.0_rp - zerosw )
1930 
1931  re(k,i,j,i_mp_qg) = 1.5_rp * rlmdg * um2cm
1932  enddo
1933  enddo
1934  enddo
1935 
1936  return
1937  end subroutine atmos_phy_mp_tomita08_effectiveradius
1938 
1939  !-----------------------------------------------------------------------------
1941  subroutine atmos_phy_mp_tomita08_mixingratio( &
1942  Qe, &
1943  QTRC0 )
1944  use scale_grid_index
1945  use scale_tracer, only: &
1946  qad => qa, &
1947  mp_qad => mp_qa
1948  implicit none
1949 
1950  real(RP), intent(out) :: Qe (ka,ia,ja,mp_qad) ! mixing ratio of each cateory [kg/kg]
1951  real(RP), intent(in) :: QTRC0(ka,ia,ja,qad) ! tracer mass concentration [kg/kg]
1952 
1953  integer :: ihydro
1954  !---------------------------------------------------------------------------
1955 
1956  do ihydro = 1, mp_qa
1957  qe(:,:,:,ihydro) = qtrc0(:,:,:,i_mp2all(ihydro))
1958  enddo
1959 
1960  return
1961  end subroutine atmos_phy_mp_tomita08_mixingratio
1962 
1963 end module scale_atmos_phy_mp_tomita08
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:494
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_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:1944
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:602
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:497
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_atmos_phy_mp_tomita08::atmos_phy_mp_tomita08_cloudfraction
subroutine, public atmos_phy_mp_tomita08_cloudfraction(cldfrac, QTRC, mask_criterion)
Calculate Cloud Fraction.
Definition: scale_atmos_phy_mp_tomita08.F90:1815
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:59
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:43
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:1850
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