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