scale_atmos_phy_mp_sn14.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------

#ifdef PROFILE_FAPP
#define PROFILE_START(name) call fapp_start(name, 1, 1)
#define PROFILE_STOP(name)  call fapp_stop (name, 1, 1)
#elif defined(PROFILE_FINEPA)
#define PROFILE_START(name) call start_collection(name)
#define PROFILE_STOP(name)  call stop_collection (name)
#else
#define PROFILE_START(name)
#define PROFILE_STOP(name)
#endif

#include "macro_thermodyn.h"
module scale_atmos_phy_mp_sn14
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index

  use scale_tracer_sn14
  use scale_tracer, only: qa
  use scale_const, only: &
     grav   => const_grav,   &
     pi     => const_pi,     &
     undef8 => const_undef8, &
     rdry   => const_rdry,   &
     cpdry  => const_cpdry,  &
     cvdry  => const_cvdry,  &
     p00    => const_pre00,  &
     t00    => const_tem00,  &
     rvap   => const_rvap,   &
     cpvap  => const_cpvap,  &
     cvvap  => const_cvvap,  &
     cl     => const_cl,     &
     ci     => const_ci,     &
     lhv    => const_lhv00,  &
     lhf    => const_lhf00,  &
     lhv0   => const_lhv0,   &
     lhf0   => const_lhf0,   &
     lhs0   => const_lhs0,   &
     lhv00  => const_lhv00,  &
     lhf00  => const_lhf00,  &
     psat0  => const_psat0,  &
     emelt  => const_emelt,  &
     dwatr  => const_dwatr
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_mp_sn14_setup
  public :: atmos_phy_mp_sn14
  public :: atmos_phy_mp_sn14_cloudfraction
  public :: atmos_phy_mp_sn14_effectiveradius
  public :: atmos_phy_mp_sn14_mixingratio

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  real(RP), public, target :: atmos_phy_mp_dens(mp_qa) ! hydrometeor density [kg/m3]=[g/L]

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: mp_sn14_init
  private :: mp_sn14
  private :: mp_terminal_velocity

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters
  !
  integer, parameter :: hydro_max  = 6    ! total number of mixing ratio of water
  integer, parameter :: i_c = 1
  integer, parameter :: i_r = 2
  integer, parameter :: i_i = 3
  integer, parameter :: i_s = 4
  integer, parameter :: i_g = 5

  ! production rate
  !   nucleation
  integer, parameter :: i_lcccn =  1
  integer, parameter :: i_ncccn =  2
  integer, parameter :: i_liccn =  3
  integer, parameter :: i_niccn =  4
  !   freezing
  integer, parameter :: i_lchom =  5
  integer, parameter :: i_nchom =  6
  integer, parameter :: i_lchet =  7
  integer, parameter :: i_nchet =  8
  integer, parameter :: i_lrhet =  9
  integer, parameter :: i_nrhet = 10
  !   melting
  integer, parameter :: i_limlt = 11
  integer, parameter :: i_nimlt = 12
  integer, parameter :: i_lsmlt = 13
  integer, parameter :: i_nsmlt = 14
  integer, parameter :: i_lgmlt = 15
  integer, parameter :: i_ngmlt = 16
  !   vapor deposition
  integer, parameter :: i_lrdep = 17
  integer, parameter :: i_nrdep = 18
  integer, parameter :: i_lidep = 19
  integer, parameter :: i_nidep = 20
  integer, parameter :: i_lsdep = 21
  integer, parameter :: i_nsdep = 22
  integer, parameter :: i_lgdep = 23
  integer, parameter :: i_ngdep = 24
  integer, parameter :: i_lcdep = 25
!    integer, parameter :: I_NCdep = 26
  !   warm collection process
  !     auto-conversion
  integer, parameter :: i_lcaut = 26
  integer, parameter :: i_ncaut = 27
  integer, parameter :: i_nraut = 28
  !     accretion
  integer, parameter :: i_lcacc = 29
  integer, parameter :: i_ncacc = 30
  !     self-colletion, break-up
  integer, parameter :: i_nrslc = 31
  integer, parameter :: i_nrbrk = 32

  !   partial conversion(ice, snow => graupel)
  integer, parameter :: i_licon = 33
  integer, parameter :: i_nicon = 34
  integer, parameter :: i_lscon = 35
  integer, parameter :: i_nscon = 36
  !   enhanced melting due to
  integer, parameter :: i_liacm = 37 ! ice-cloud
  integer, parameter :: i_niacm = 38
  integer, parameter :: i_liarm = 39 ! ice-rain
  integer, parameter :: i_niarm = 40
  integer, parameter :: i_lsacm = 41 ! snow-cloud
  integer, parameter :: i_nsacm = 42
  integer, parameter :: i_lsarm = 43 ! snow-rain
  integer, parameter :: i_nsarm = 44
  integer, parameter :: i_lgacm = 45 ! graupel-cloud
  integer, parameter :: i_ngacm = 46
  integer, parameter :: i_lgarm = 47 ! graupel-rain
  integer, parameter :: i_ngarm = 48
  !   ice multiplication by splintering
  integer, parameter :: i_lgspl = 49
  integer, parameter :: i_lsspl = 50
  integer, parameter :: i_nispl = 51

  integer, parameter :: pq_max = 51

  ! production rate of mixed-phase collection process
  ! PXXacYY2ZZ means XX collect YY produce ZZ
  integer, parameter :: i_liaclc2li   =  1 ! cloud-ice
  integer, parameter :: i_niacnc2ni   =  2
  integer, parameter :: i_lsaclc2ls   =  3 ! cloud-snow(cloud change)
  integer, parameter :: i_nsacnc2ns   =  4
  integer, parameter :: i_lgaclc2lg   =  5 ! cloud-graupel
  integer, parameter :: i_ngacnc2ng   =  6
  integer, parameter :: i_lracli2lg_i =  7 ! rain-ice(ice change)
  integer, parameter :: i_nracni2ng_i =  8
  integer, parameter :: i_lracli2lg_r =  9 ! rain-ice(rain change)
  integer, parameter :: i_nracni2ng_r = 10
  integer, parameter :: i_lracls2lg_s = 11 ! rain-snow(snow change)
  integer, parameter :: i_nracns2ng_s = 12
  integer, parameter :: i_lracls2lg_r = 13 ! rain-snow(rain change)
  integer, parameter :: i_nracns2ng_r = 14
  integer, parameter :: i_lraclg2lg   = 15 ! rain-graupel(rain change)
  integer, parameter :: i_nracng2ng   = 16
  integer, parameter :: i_liacli2ls   = 17 ! ice-ice
  integer, parameter :: i_niacni2ns   = 18
  integer, parameter :: i_liacls2ls   = 19 ! ice-snow(ice change)
  integer, parameter :: i_niacns2ns   = 20
  integer, parameter :: i_nsacns2ns   = 21 ! snow-snow
  integer, parameter :: i_ngacng2ng   = 22 ! graupel-graupel
  integer, parameter :: i_lgacls2lg   = 23 ! snow-graupel
  integer, parameter :: i_ngacns2ng   = 24

  integer, parameter :: pac_max       = 24



  character(len=H_SHORT), save :: wlabel(11)

  ! empirical value from Meyers etal.(1991), 1[/liter] = 1.d3[/m3]
  real(RP), private, parameter :: nqmin(6) = (/ 0.0_rp, 1.e+4_rp, 1.0_rp, 1.0_rp, 1.e-4_rp, 1.e-4_rp /) ! [1/m3]
  ! refer to Seifert(2002) (Dr. Thesis, Table.5.1)
  ! max mass, for D_min=79um, 2mm, 5mm, 1cm, 1cm
  real(RP), private, parameter :: xqmax(6) = (/ 0.0_rp, 2.6e-10_rp, 5.0e-6_rp, 1.377e-6_rp, 7.519e-6_rp, 4.90e-5_rp /)! [kg]
  ! SB06, Table 1.
  ! min mass, for D_min=2um, 79um, 10um, 20um, 100um
  real(RP), private, parameter :: xqmin(6) = (/ 0.0_rp, 4.20e-15_rp, 2.60e-10_rp, 3.382e-13_rp, 1.847e-12_rp, 1.230e-10_rp /)! [kg]



  ! for all processes
  ! SB06, Table 1.
  real(RP), private, parameter :: xc_min = 4.20e-15_rp   ! [kg] : min mass, D_min=2um
  real(RP), private, parameter :: xr_min = 2.60e-10_rp   ! [kg] : min mass, D_min=79um
  real(RP), private, parameter :: xi_min = 3.382e-13_rp  ! [kg] : min mass, D_min=10um
  real(RP), private, parameter :: xs_min = 1.847e-12_rp  ! [kg] : min mass, D_min=20um
  real(RP), private, parameter :: xg_min = 1.230e-10_rp  ! [kg] : min mass, D_min=100um
  ! refer to Seifert(2002) (Dr. Thesis, Table.5.1)
  real(RP), private, parameter :: xc_max = 2.6e-10_rp    ! [kg] : max, D_max=79um
  real(RP), private, parameter :: xr_max = 5.00e-6_rp    ! [kg] : max, D_max=2mm
  real(RP), private, parameter :: xi_max = 1.377e-6_rp   ! [kg] : max, D_max=5mm
  real(RP), private, parameter :: xs_max = 7.519e-6_rp   ! [kg] : max, D_max=1cm
  real(RP), private, parameter :: xg_max = 4.900e-5_rp   ! [kg] : max, D_max=1cm
  ! filter similar to Ikawa et al.(1991) sec.3.5
  real(RP), private, parameter :: xmin_filter= xc_min
  ! filter of effective radius(1 micron)
  real(RP), private, parameter :: rmin_re= 1.e-6_rp
  !
  ! SB06(95),(96)
  real(RP), private, parameter :: n0r_min= 2.5e+5_rp    ! [m-4]: min intercept parameter of rain
  real(RP), private, parameter :: n0r_max= 2.0e+7_rp    ! [m-4]: max
  real(RP), private, parameter :: lambdar_min= 1.e+3_rp ! [m-1]: min slope parameter of rain
  real(RP), private, parameter :: lambdar_max= 1.e+4_rp ! [m-1]: max
  ! empirical value from Meyers etal.(1991), 1[/liter] = 1.d3[/m3]
  real(RP), private, parameter :: nc_min = 1.e+4_rp     ! [m-3] empirical T.Mitsui
  real(RP), private, parameter :: nr_min = 1.0_rp     ! [m-3] 1/1000 [/liter]
  real(RP), private, parameter :: ni_min = 1.0_rp     ! [m-3]
  real(RP), private, parameter :: ns_min = 1.e-4_rp    ! [m-3]
  real(RP), private, parameter :: ng_min = 1.e-4_rp    ! [m-3]
  ! empirical filter
  real(RP), private, parameter :: lc_min = xc_min*nc_min
  real(RP), private, parameter :: lr_min = xr_min*nr_min
  real(RP), private, parameter :: li_min = xi_min*ni_min
  real(RP), private, parameter :: ls_min = xs_min*ns_min
  real(RP), private, parameter :: lg_min = xg_min*ng_min
  !
  real(RP), private, parameter :: x_sep   = 2.6e-10_rp ! boundary mass between cloud and rain
  !
  real(RP), private, parameter :: tem_min=100.0_rp
  real(RP), private, parameter :: rho_min=1.e-5_rp     ! 3.e-3 is lower limit recognized in many experiments.
  real(RP), private, parameter :: rhoi   = 916.70_rp
  ! for Seifert(2008)
  ! work parameter for gamma function, imported from MISC_gammafunc
  real(RP), private, parameter :: gfac_coef(6)=(/&
       +76.180091729471460_rp, -86.505320329416770_rp, +24.014098240830910_rp,&
       -1.2317395724501550_rp, +0.1208650973866179e-2_rp, -0.5395239384953e-5_rp /)
  real(RP), private, parameter :: gfac_ser0=1.000000000190015_rp
  !
  integer, private, save :: ntmax_phase_change = 1
  integer, private, save :: ntmax_collection   = 1
  integer, private, save :: mp_ntmax_sedimentation= 1 ! 10/08/03 [Add] T.Mitsui
  !
  !--- standard density
  real(RP), private, parameter :: rho_0 = 1.280_rp
  !--- max number of Nc( activatable aerosol number concentration )
  real(RP), allocatable, private, save :: nc_uplim_d(:,:,:)
  !
  !--- thermal conductivity of air
  real(RP), private, parameter :: ka0  = 2.428e-2_rp
  
  real(RP), private, parameter :: dka_dt = 7.47e-5_rp
  
  !====== Ka = Ka0 + temc*dKa_dT
  !
  !--- Dynamic viscosity
  real(RP), private, parameter :: mua0 = 1.718e-5_rp
  
  real(RP), private, parameter :: dmua_dt = 5.28e-8_rp
  
  !======  mua = mua0 + temc*dmua_dT
  !
  real(RP), private, save :: xc_ccn = 1.e-12_rp ! [kg]
  real(RP), private, save :: xi_ccn = 1.e-12_rp ! [kg] ! [move] 11/08/30 T.Mitsui
  !
  ! capacity of diffusional growth
  ! ( dependent of their geometries )
  real(RP), private, save :: cap(hydro_max)
  !
  ! constants for Diameter-Mass relation
  ! D = a * x^b
  real(RP), private, save :: a_m(hydro_max)
  real(RP), private, save :: b_m(hydro_max)
  ! constants for Terminal velocity-Mass relation
  ! vt = alpha * x^beta * f
  real(RP), private, save :: alpha_v(hydro_max,2)
  real(RP), private, save :: beta_v(hydro_max,2)
  real(RP), private, save :: alpha_vn(hydro_max,2)   !
  real(RP), private, save :: beta_vn(hydro_max,2)    !
  real(RP), private, save :: gamma_v(hydro_max)
  !  Aerodynamical factor for correction of terminal velocity.(Heymsfield and Iaquinta, 2000)
  !  vt(tem,pre) = vt0 * (pre/pre0)**a_pre0 * (tem/tem0)**a_tem0
  real(RP), private, parameter :: pre0_vt   = 300.e+2_rp  ! 300hPa
  real(RP), private, parameter :: tem0_vt   = 233.0_rp  ! -40degC
  real(RP), private, parameter :: a_pre0_vt = -0.1780_rp
  real(RP), private, parameter :: a_tem0_vt = -0.3940_rp
  ! Parameters to determine Droplet Size Distribution
  ! as a General Gamma Distribution
  ! f(x) = A x^nu exp(-lambda x^mu )
  ! for Marshall Palmer Distribution ( popular for rain )
  !     mu=1/3, nu=-2/3
  ! for Gamma Distribution ( popular for cloud )
  !     mu=1
  real(RP), private, save :: nu(hydro_max)
  real(RP), private, save :: mu(hydro_max)
  ! Mitchell(1996), JAS, vol.53, No.12, pp.1710-1723
  !  area = a_area*D^b_area
  !  area = ax_area*x^bx_area
  ! Auer and Veal(1970), JAS, vol.27, pp.919-pp.926
  ! height = a_h*x^b_h( based on h=a_ar*D^b_ar,  ar:aspect ratio)
  real(RP), private, save :: a_area(hydro_max)       !
  real(RP), private, save :: b_area(hydro_max)       !
  real(RP), private, save :: ax_area(hydro_max)      !
  real(RP), private, save :: bx_area(hydro_max)      !
  ! parameters for radius of equivalent area
  ! r_ea = a_rea*x**b_rea
  real(RP), private, save :: a_rea(hydro_max)        !
  real(RP), private, save :: b_rea(hydro_max)        !
  real(RP), private, save :: a_rea2(hydro_max)       !
  real(RP), private, save :: b_rea2(hydro_max)       !
  real(RP), private, save :: a_rea3(hydro_max)       !
  real(RP), private, save :: b_rea3(hydro_max)       !
  !
  real(RP), private, save :: a_d2vt(hydro_max)       !
  real(RP), private, save :: b_d2vt(hydro_max)       !
  ! coefficient of x^2 moment of DSD
  ! Z = integral x*x*f(x) dx
  !   = coef_m2*N*(L/N)^2
  real(RP), private, save :: coef_m2(hydro_max)
  ! radar reflectivity coefficient defined by diameter
  real(RP), private, save :: coef_d6(hydro_max)      !
  ! volume coefficient defined by diameter
  real(RP), private, save :: coef_d3(hydro_max)      !
  ! coefficient of weighted mean diameter
  real(RP), private, save :: coef_d(hydro_max)
  ! coefficient of weighted mean d*d*v
  real(RP), private, save :: coef_d2v(hydro_max)     !
  ! coefficient of moment of d*d*v
  real(RP), private, save :: coef_md2v(hydro_max)    !
  !
  ! for effective radius(spherical particle)
  real(RP), private, save :: coef_r2(hydro_max)
  real(RP), private, save :: coef_r3(hydro_max)
  real(RP), private, save :: coef_re(hydro_max)
  ! for effective radius(hexagonal plate)
  real(RP), private, save :: coef_rea2(hydro_max)    !
  real(RP), private, save :: coef_rea3(hydro_max)    !
  logical, private, save :: opt_m96_ice=.true.               !
  logical, private, save :: opt_m96_column_ice=.false.       !
  !
  ! coefficeint of weighted mean terminal velocity
  ! vt0 is number weighted and
  ! vt1 is mass   weighted.
  real(RP), private, save :: coef_vt0(hydro_max,2)
  real(RP), private, save :: coef_vt1(hydro_max,2)
  real(RP), private, save :: coef_deplc
  real(RP), private, save :: coef_dave_n(hydro_max) !
  real(RP), private, save :: coef_dave_l(hydro_max) !
  ! diameter of terminal velocity branch
  !
  real(RP), private, save :: d0_ni=261.76e-6_rp   !
  real(RP), private, save :: d0_li=398.54e-6_rp
  real(RP), private, parameter     :: d0_ns=270.03e-6_rp   !
  real(RP), private, parameter     :: d0_ls=397.47e-6_rp   !
  real(RP), private, parameter     :: d0_ng=269.08e-6_rp   !
  real(RP), private, parameter     :: d0_lg=376.36e-6_rp   !
  !
  real(RP), private, parameter :: coef_vtr_ar1=9.65_rp    ! coef. for large branch
  ! original parameter of Rogers etal.(1993)
  real(RP), private, parameter :: coef_vtr_br1=10.43_rp    ! ...
  real(RP), private, parameter :: coef_vtr_cr1=600.0_rp    ! ...
  real(RP), private, parameter :: coef_vtr_ar2=4.e+3_rp     ! coef. for small branch
  real(RP), private, parameter :: coef_vtr_br2=12.e+3_rp     ! ...
  real(RP), private, parameter :: d_vtr_branch=0.745e-3_rp  ! 0.745 mm (diameter dividing 2-branches)
  ! equilibrium diameter of rain break-up
  real(RP), private, parameter :: dr_eq   = 1.10e-3_rp ! eqilibrium diameter, Seifert 2008(36)
  ! coefficient of General Gamma.
  !  f(x)  = A x^nu exp(-lambda x^mu )
  ! lambda = coef_lambda * (L/N)^{-mu}
  !     A  = coef_A*N*lambda^slope_A
  real(RP), private, save :: coef_a(hydro_max)
  real(RP), private, save :: coef_lambda(hydro_max)
!  real(RP), private, save :: slope_A(HYDRO_MAX)
  ! coefficeint of weighted ventilation effect.
  ! large, and small branch is by PK97(13-60),(13-61),(13-88),(13-89)
  real(RP), private, save :: ah_vent  (hydro_max,2) !
  real(RP), private, save :: bh_vent  (hydro_max,2) !
  real(RP), private, save :: ah_vent0 (hydro_max,2) !
  real(RP), private, save :: bh_vent0 (hydro_max,2) !
  real(RP), private, save :: ah_vent1 (hydro_max,2) !
  real(RP), private, save :: bh_vent1 (hydro_max,2) !
  ! coefficient of collision growth
  real(RP), private, save :: delta_b0 (hydro_max)
  real(RP), private, save :: delta_b1 (hydro_max)
  real(RP), private, save :: delta_ab0(hydro_max,hydro_max)
  real(RP), private, save :: delta_ab1(hydro_max,hydro_max)
  !
  real(RP), private, save :: theta_b0 (hydro_max)
  real(RP), private, save :: theta_b1 (hydro_max)
  real(RP), private, save :: theta_ab0(hydro_max,hydro_max)
  real(RP), private, save :: theta_ab1(hydro_max,hydro_max)
  !
  logical, private, save :: opt_debug=.false.
  !
  logical, private, save :: opt_debug_tem=.false.
  logical, private, save :: opt_debug_inc=.true.
  logical, private, save :: opt_debug_act=.true.
  logical, private, save :: opt_debug_ree=.true.
  logical, private, save :: opt_debug_bcs=.true.

  integer, private, save :: mp_nstep_sedimentation
  real(RP), private, save :: mp_rnstep_sedimentation
  real(DP), private, save :: mp_dtsec_sedimentation

  !
  ! metrics of vertical coordinate
  ! not used in SCALE-RM
  !
  real(RP), private, allocatable, save :: gsgam2_d (:,:,:)
  real(RP), private, allocatable, save :: gsgam2h_d(:,:,:)
  real(RP), private, allocatable, save :: gam2_d   (:,:,:)
  real(RP), private, allocatable, save :: gam2h_d  (:,:,:)
  real(RP), private, allocatable, save :: rgsgam2_d(:,:,:)
  real(RP), private, allocatable, save :: rgs_d    (:,:,:)
  real(RP), private, allocatable, save :: rgsh_d   (:,:,:)

  logical, private, save :: mp_doautoconversion = .true.
  logical, private, save :: mp_doprecipitation  = .true.
  logical, private, save :: mp_couple_aerosol   = .false. ! apply CCN effect?
  real(RP), private, save :: mp_ssw_lim = 1.e+1_rp


  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_sn14_setup( MP_TYPE )
    use scale_process, only: &
       prc_mpistop
    use scale_grid, only: &
       cdz => grid_cdz
    use scale_const, only: &
       const_dwatr, &
       const_dice
    use scale_time, only: &
       time_dtsec_atmos_phy_mp
    implicit none

    character(len=*), intent(in) :: MP_TYPE

    namelist / param_atmos_phy_mp / &
       mp_doautoconversion, &
       mp_doprecipitation,  &
       mp_ssw_lim,          &
       mp_couple_aerosol,   &
       mp_ntmax_sedimentation

    real(RP), parameter :: max_term_vel = 10.0_rp  !-- terminal velocity for calculate dt of sedimentation
    integer :: nstep_max
    integer :: ierr
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '+++ Module[Cloud Microphisics]/Categ[ATMOS]'
    if( io_l ) write(io_fid_log,*) '*** Wrapper for SN14'

    if ( mp_type /= 'SN14' ) then
       write(*,*) 'xxx ATMOS_PHY_MP_TYPE is not SN14. Check!'
       call prc_mpistop
    end if

    if (      i_qv <= 0 &
         .OR. i_qc <= 0 &
         .OR. i_qr <= 0 &
         .OR. i_qi <= 0 &
         .OR. i_qs <= 0 &
         .OR. i_qg <= 0 &
         .OR. i_nc <= 0 &
         .OR. i_nr <= 0 &
         .OR. i_ni <= 0 &
         .OR. i_ns <= 0 &
         .OR. i_ng <= 0 ) then
       write(*,*) 'xxx SN14 needs QV/C/R/I/S/G and NC/R/I/S/G tracer. Check!'
       call prc_mpistop
    endif

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_mp,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP. Check!'
       call prc_mpistop
    endif
    if( io_lnml ) write(io_fid_log,nml=param_atmos_phy_mp)

    atmos_phy_mp_dens(i_mp_qc) = const_dwatr
    atmos_phy_mp_dens(i_mp_qr) = const_dwatr
    atmos_phy_mp_dens(i_mp_qi) = const_dice
    atmos_phy_mp_dens(i_mp_qs) = const_dice
    atmos_phy_mp_dens(i_mp_qg) = const_dice

    wlabel( 1) = "VAPOR"
    wlabel( 2) = "CLOUD"
    wlabel( 3) = "RAIN"
    wlabel( 4) = "ICE"
    wlabel( 5) = "SNOW"
    wlabel( 6) = "GRAUPEL"
    wlabel( 7) = "CLOUD_NUM"
    wlabel( 8) = "RAIN_NUM"
    wlabel( 9) = "ICE_NUM"
    wlabel(10) = "SNOW_NUM"
    wlabel(11) = "GRAUPEL_NUM"

    call mp_sn14_init

    allocate(nc_uplim_d(1,ia,ja))
    nc_uplim_d(:,:,:) = 150.e6_rp

    nstep_max = int( ( time_dtsec_atmos_phy_mp * max_term_vel ) / minval( cdz ) )
    mp_ntmax_sedimentation = max( mp_ntmax_sedimentation, nstep_max )

    mp_nstep_sedimentation  = mp_ntmax_sedimentation
    mp_rnstep_sedimentation = 1.0_rp / real(mp_ntmax_sedimentation,kind=rp)
    mp_dtsec_sedimentation  = time_dtsec_atmos_phy_mp * mp_rnstep_sedimentation

    if ( io_l ) write(io_fid_log,*)
    if ( io_l ) write(io_fid_log,*) '*** Timestep of sedimentation is divided into : ', mp_ntmax_sedimentation, ' step'
    if ( io_l ) write(io_fid_log,*) '*** DT of sedimentation is : ', mp_dtsec_sedimentation, '[s]'

    !--- For kij
    allocate( gsgam2_d(ka,ia,ja) )
    allocate( gsgam2h_d(ka,ia,ja) )
    allocate( gam2_d(ka,ia,ja) )
    allocate( gam2h_d(ka,ia,ja) )
    allocate( rgsgam2_d(ka,ia,ja) )
    allocate( rgs_d(ka,ia,ja) )
    allocate( rgsh_d(ka,ia,ja) )
    gsgam2_d(:,:,:) = 1.0_rp
    gsgam2h_d(:,:,:) = 1.0_rp
    gam2_d(:,:,:) = 1.0_rp
    gam2h_d(:,:,:) = 1.0_rp
    rgsgam2_d(:,:,:) = 1.0_rp
    rgs_d(:,:,:) = 1.0_rp
    rgsh_d(:,:,:) = 1.0_rp

    return
  end subroutine atmos_phy_mp_sn14_setup

  !-----------------------------------------------------------------------------
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_sn14( &
       DENS,      &
       MOMZ,      &
       MOMX,      &
       MOMY,      &
       RHOT,      &
       QTRC,      &
       CCN,       &
       EVAPORATE, &
       SFLX_rain, &
       SFLX_snow  )
    use scale_grid_index
    use scale_tracer, only: &
       qad => qa, &
       mp_qad => mp_qa
    implicit none

    real(RP), intent(inout) :: dens(ka,ia,ja)
    real(RP), intent(inout) :: MOMZ(ka,ia,ja)
    real(RP), intent(inout) :: MOMX(ka,ia,ja)
    real(RP), intent(inout) :: MOMY(ka,ia,ja)
    real(RP), intent(inout) :: RHOT(ka,ia,ja)
    real(RP), intent(inout) :: QTRC(ka,ia,ja,qad)
    real(RP), intent(in)    :: CCN(ka,ia,ja)
    real(RP), intent(out)   :: EVAPORATE(ka,ia,ja)
    real(RP), intent(out)   :: SFLX_rain(ia,ja)
    real(RP), intent(out)   :: SFLX_snow(ia,ja)
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*) '*** Physics step: Cloud microphysics(SN14)'

#ifdef PROFILE_FIPP
    call fipp_start()
#endif

    call mp_negativefilter( dens, qtrc )

    call mp_sn14( dens,      & ! [INOUT]
                  momz,      & ! [INOUT]
                  momx,      & ! [INOUT]
                  momy,      & ! [INOUT]
                  rhot,      & ! [INOUT]
                  qtrc,      & ! [INOUT]
                  ccn,       & ! [IN]
                  evaporate, & ! [OUT]
                  sflx_rain, & ! [OUT]
                  sflx_snow  ) ! [OUT]

    call mp_negativefilter( dens, qtrc )

#ifdef PROFILE_FIPP
    call fipp_stop()
#endif

    return
  end subroutine atmos_phy_mp_sn14

  !-----------------------------------------------------------------------------
  subroutine mp_sn14_init
    use scale_process, only: &
       prc_mpistop
    use scale_specfunc, only: &
        gammafunc => sf_gamma
    implicit none

    real(RP), allocatable :: w1(:),w2(:),w3(:),w4(:),w5(:),w6(:),w7(:),w8(:)
    ! work for calculation of capacity, Mitchell and Arnott (1994) , eq.(9)
    real(RP) :: ar_ice_fix = 0.7_rp
    real(RP) :: wcap1, wcap2
    ! work for ventilation coefficient
    logical :: flag_vent0(hydro_max), flag_vent1(hydro_max)
    integer :: ierr
    integer :: iw, ia, ib
    integer :: n
    !
    namelist /nm_mp_sn14_init/       &
         opt_debug,                  &
         opt_debug_tem,              &
         opt_debug_inc,              &
         opt_debug_act,              &
         opt_debug_ree,              &
         opt_debug_bcs,              &
         ntmax_phase_change,         &
         ntmax_collection
    !
    namelist /nm_mp_sn14_particles/ &
         a_m, b_m, alpha_v, beta_v, gamma_v, &
         alpha_vn, beta_vn,    &
         a_area, b_area, cap,  &
         nu, mu,               &
         opt_m96_column_ice,   &
         opt_m96_ice,          &
         ar_ice_fix
    real(RP), parameter :: eps_gamma=1.e-30_rp

    a_m(:)         = undef8
    b_m(:)         = undef8
    alpha_v(:,:)   = undef8
    beta_v(:,:)    = undef8
    alpha_vn(:,:)  = undef8
    beta_vn(:,:)   = undef8
    gamma_v(:)     = undef8
    a_d2vt(:)      = undef8
    b_d2vt(:)      = undef8
    a_area(:)      = undef8
    b_area(:)      = undef8
    ax_area(:)     = undef8
    bx_area(:)     = undef8
    a_rea(:)       = undef8
    b_rea(:)       = undef8
    a_rea2(:)      = undef8
    b_rea2(:)      = undef8
    a_rea3(:)      = undef8
    b_rea3(:)      = undef8
    nu(:)          = undef8
    mu(:)          = undef8
    cap(:)         = undef8
    coef_m2(:)     = undef8
    coef_dave_n(:) = undef8
    coef_dave_l(:) = undef8
    coef_d(:)      = undef8
    coef_d3(:)     = undef8
    coef_d6(:)     = undef8
    coef_d2v(:)    = undef8
    coef_md2v(:)   = undef8
    coef_r2(:)     = undef8
    coef_r3(:)     = undef8
    coef_re(:)     = undef8
    coef_rea2(:)   = undef8
    coef_rea3(:)   = undef8
    coef_a(:)      = undef8
!    slope_A(:)     = UNDEF8
    coef_lambda(:) = undef8
    coef_vt0(:,:)  = undef8
    coef_vt1(:,:)  = undef8
    delta_b0(:)    = undef8
    delta_b1(:)    = undef8
    delta_ab0(:,:) = undef8
    delta_ab1(:,:) = undef8
    theta_b0(:)    = undef8
    theta_b1(:)    = undef8
    theta_ab0(:,:) = undef8
    theta_ab1(:,:) = undef8
    !
    ah_vent(:,:)   = undef8
    ah_vent0(:,:)  = undef8
    ah_vent1(:,:)  = undef8
    bh_vent(:,:)   = undef8
    bh_vent0(:,:)  = undef8
    bh_vent1(:,:)  = undef8
    !

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '+++ Module[SN14]/Categ[ATMOS]'

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=nm_mp_sn14_init,iostat=ierr)

    if( ierr < 0 ) then !--- missing
       if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist nm_mp_sn14_init. Check!'
       call prc_mpistop
    endif
    if( io_l ) write(io_fid_log,nml=nm_mp_sn14_init)

    !
    ! default setting
    !
    ! Area parameters with mks unit originated by Mitchell(1996)
    a_area(i_qc) = pi/4.0_rp ! sphere
    a_area(i_qr) = pi/4.0_rp ! sphere
    a_area(i_qi) = 0.65_rp*1.e-4_rp*100.0_rp**(2.00_rp)   ! Mitchell(1996), Hexagonal Plate
    a_area(i_qs) = 0.2285_rp*1.e-4_rp*100.0_rp**(1.88_rp) ! Mitchell(1996), Aggregates
    a_area(i_qg) = 0.50_rp*1.e-4_rp*100.0_rp**(2.0_rp)    ! Mitchell(1996), Lump Graupel
    b_area(i_qc) = 2.0_rp
    b_area(i_qr) = 2.0_rp
    b_area(i_qi) = 2.0_rp
    b_area(i_qs) = 1.88_rp
    b_area(i_qg) = 2.0_rp
    !
    ! Seifert and Beheng(2006), Table. 1 or List of symbols
    !----------------------------------------------------------
    ! Diameter-Mass relationship
    ! D = a * x^b
    a_m(i_qc) = 0.124_rp
    a_m(i_qr) = 0.124_rp
    a_m(i_qi) = 0.217_rp
    a_m(i_qs) = 8.156_rp
    a_m(i_qg) = 0.190_rp
    b_m(i_qc) = 1.0_rp/3.0_rp
    b_m(i_qr) = 1.0_rp/3.0_rp
    b_m(i_qi) = 0.302_rp
    b_m(i_qs) = 0.526_rp
    b_m(i_qg) = 0.323_rp
    !----------------------------------------------------------
    ! Terminal velocity-Mass relationship
    ! vt = alpha * x^beta * (rho0/rho)^gamma
    alpha_v(i_qc,:)= 3.75e+5_rp
    alpha_v(i_qr,:)= 159.0_rp ! not for sedimantation
    alpha_v(i_qi,:)= 317.0_rp
    alpha_v(i_qs,:)= 27.70_rp
    alpha_v(i_qg,:)= 40.0_rp
    beta_v(i_qc,:) = 2.0_rp/3.0_rp
    beta_v(i_qr,:) = 0.266_rp ! not for sedimantation
    beta_v(i_qi,:) = 0.363_rp
    beta_v(i_qs,:) = 0.216_rp
    beta_v(i_qg,:) = 0.230_rp
    gamma_v(i_qc)  = 1.0_rp
    ! This is high Reynolds number limit(Beard 1980)
    gamma_v(i_qr)  = 1.0_rp/2.0_rp
    gamma_v(i_qi)  = 1.0_rp/2.0_rp
    gamma_v(i_qs)  = 1.0_rp/2.0_rp
    gamma_v(i_qg)  = 1.0_rp/2.0_rp
    !----------------------------------------------------------
    ! DSD parameters
    ! f(x) = A x^nu exp( -lambda x^mu )
    ! Gamma Disribution           : mu=1  , nu:arbitrary
    ! Marshall-Palmer Distribution: mu=1/3, nu:-2/3
    ! In the case of MP, f(D) dD = f(x)dx
    !                    f(x)    = c * f(D)/D^2 (c:coefficient)
    nu(i_qc) =  1.0_rp          ! arbitrary for Gamma
    nu(i_qr) = -1.0_rp/3.0_rp     ! nu(diameter)=1, equilibrium condition.
    nu(i_qi) =  1.0_rp          !
    nu(i_qs) =  1.0_rp          !
    nu(i_qg) =  1.0_rp          !
    !
    mu(i_qc) = 1.0_rp           ! Gamma
    mu(i_qr) = 1.0_rp/3.0_rp      ! Marshall Palmer
    mu(i_qi) = 1.0_rp/3.0_rp      !
    mu(i_qs) = 1.0_rp/3.0_rp      !
    mu(i_qg) = 1.0_rp/3.0_rp      !
    !----------------------------------------------------------
    ! Geomeries for diffusion growth
    ! Pruppacher and Klett(1997), (13-77)-(13-80) and
    ! originally derived by McDonald(1963b)
    ! sphere: cap=2
    ! plate : cap=pi
    ! needle with aspect ratio a/b
    !       : cap=log(2*a/b)
    cap(i_qc) = 2.0_rp    ! sphere
    cap(i_qr) = 2.0_rp    ! sphere
    cap(i_qi) = pi ! hexagonal plate
    cap(i_qs) = 2.0_rp    ! mix aggregates
    cap(i_qg) = 2.0_rp    ! lump
    !
    alpha_vn(:,:) = alpha_v(:,:)
    beta_vn(:,:) = beta_v(:,:)
    !------------------------------------------------------------------------
    !
    ! additional setting
    !

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=nm_mp_sn14_particles,iostat=ierr)

    if( ierr < 0 ) then !--- missing
       if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist nm_mp_sn14_particles. Check!'
       call prc_mpistop
    endif
    if( io_l ) write(io_fid_log,nml=nm_mp_sn14_particles)

    ! [Add] 10/08/03 T.Mitsui
    ! particles shapes are
    if( opt_m96_ice ) then
       ! ice is randomly oriented Hexagonal plate (Auer and Veal 1970, Takano and Liou 1995, Mitchell 1996)
       ! snow is assemblages of planar polycrystals(Mitchell 1996)
       ! graupel is Lump graupel(R4b) is assumed(Mitchell 1996)
       a_area(i_qi)   = 0.120284936_rp
       a_area(i_qs)   = 0.131488_rp
       a_area(i_qg)   = 0.5_rp
       b_area(i_qi)   = 1.850000_rp
       b_area(i_qs)   = 1.880000_rp
       b_area(i_qg)   = 2.0_rp
       a_m(i_qi)      = 1.23655360084766_rp
       a_m(i_qs)      = a_m(i_qi)
       a_m(i_qg)      = 0.346111225718402_rp
       b_m(i_qi)      = 0.408329930583912_rp
       b_m(i_qs)      = b_m(i_qi)
       b_m(i_qg)      = 0.357142857142857_rp
       !
       if( opt_m96_column_ice )then
          d0_ni=240.49e-6_rp   ! this is column
          d0_li=330.09e-6_rp   ! this is column
          a_area(i_qi)= (0.684_rp*1.e-4_rp)*10.0_rp**(2.0_rp*2.00_rp)
          b_area(i_qi)= 2.0_rp
          a_m(i_qi)   = 0.19834046116844_rp
          b_m(i_qi)   = 0.343642611683849_rp
          ! [Add] 11/08/30 T.Mitsui
          ! approximated by the capacity for prolate spheroid with constant aspect ratio
          wcap1       = sqrt(1.0_rp-ar_ice_fix**2)
          wcap2       = log( (1.0_rp+wcap1)/ar_ice_fix )
          cap(i_qi)   = 2.0_rp*wcap2/wcap1
          !
       end if
       !
       ! These value are derived by least-square fitting in the range
       ! qi [100um:1000um] in diameter
       ! qs [100um:1000um] in diameter
       ! qg [200um:2000um] in diameter
       !                      small branch          , large branch
       alpha_v(i_qi,:) =(/ 5798.60107421875_rp, 167.347076416016_rp/)
       alpha_vn(i_qi,:) =(/ 12408.177734375_rp, 421.799865722656_rp/)
       if( opt_m96_column_ice )then
          alpha_v(i_qi,:) = (/2901.0_rp, 32.20_rp/)
          alpha_vn(i_qi,:) = (/9675.2_rp, 64.16_rp/)
       end if
       alpha_v(i_qs,:)    =(/ 15173.3916015625_rp, 305.678619384766_rp/)
       alpha_vn(i_qs,:)    =(/ 29257.1601562500_rp, 817.985717773438_rp/)
       alpha_v(i_qg,:)    =(/ 15481.6904296875_rp, 311.642242431641_rp/)
       alpha_vn(i_qg,:)    =(/ 27574.6562500000_rp, 697.536132812500_rp/)
       !
       beta_v(i_qi,:)  =(/ 0.504873454570770_rp, 0.324817866086960_rp/)
       beta_vn(i_qi,:)  =(/ 0.548495233058929_rp, 0.385287821292877_rp/)
       if( opt_m96_column_ice )then
          beta_v(i_qi,:)  =(/ 0.465552181005478_rp, 0.223826110363007_rp/)
          beta_vn(i_qi,:)  =(/ 0.530453503131866_rp, 0.273761242628098_rp/)
       end if
       beta_v(i_qs,:)     =(/ 0.528109610080719_rp, 0.329863965511322_rp/)
       beta_vn(i_qs,:)     =(/ 0.567154467105865_rp, 0.393876969814301_rp/)
       beta_v(i_qg,:)     =(/ 0.534656763076782_rp, 0.330253750085831_rp/)
       beta_vn(i_qg,:)     =(/ 0.570551633834839_rp, 0.387124240398407_rp/)
    end if
    !
    ! area-diameter relation => area-mass relation
    ax_area(i_qc:i_qg) = a_area(i_qc:i_qg)*a_m(i_qc:i_qg)**b_area(i_qc:i_qg)
    bx_area(i_qc:i_qg) = b_area(i_qc:i_qg)*b_m(i_qc:i_qg)
    !
    ! radius of equivalent area - m ass relation
    ! pi*rea**2 = ax_area*x**bx_area
    a_rea(i_qc:i_qg)   = sqrt(ax_area(i_qc:i_qg)/pi)
    b_rea(i_qc:i_qg)   = bx_area(i_qc:i_qg)/2.0_rp
    a_rea2(i_qc:i_qg)  = a_rea(i_qc:i_qg)**2
    b_rea2(i_qc:i_qg)  = b_rea(i_qc:i_qg)*2.0_rp
    a_rea3(i_qc:i_qg)  = a_rea(i_qc:i_qg)**3
    b_rea3(i_qc:i_qg)  = b_rea(i_qc:i_qg)*3.0_rp
    !
    a_d2vt(i_qc:i_qg)=alpha_v(i_qc:i_qg,2)*(1.0_rp/alpha_v(i_qc:i_qg,2))**(beta_v(i_qc:i_qg,2)/b_m(i_qc:i_qg))
    b_d2vt(i_qc:i_qg)=(beta_v(i_qc:i_qg,2)/b_m(i_qc:i_qg))
    !
    ! Calculation of Moment Coefficient
    !
    allocate( w1(i_qc:i_qg), w2(i_qc:i_qg), w3(i_qc:i_qg), w4(i_qc:i_qg) )
    allocate( w5(i_qc:i_qg), w6(i_qc:i_qg), w7(i_qc:i_qg), w8(i_qc:i_qg) )
    w1(:) = 0.0_rp
    w2(:) = 0.0_rp
    w3(:) = 0.0_rp
    w4(:) = 0.0_rp
    w5(:) = 0.0_rp
    w6(:) = 0.0_rp
    w7(:) = 0.0_rp
    w8(:) = 0.0_rp
    !-------------------------------------------------------
    ! moment coefficient
    ! SB06 (82)
    ! M^n = coef_mn * N * (L/N)**n
    ! M^2 = Z = coef_m2 * N *(L/N)**2
    ! a*M^b = a*integral x^b f(x) dx = ave D
    do iw=i_qc,i_qg
       n = 2
       w1(iw) = gammafunc( (n+nu(iw)+1.0_rp)/mu(iw) )
       w2(iw) = gammafunc( (nu(iw)+1.0_rp)/mu(iw) )
       w3(iw) = gammafunc( (nu(iw)+2.0_rp)/mu(iw) )
       coef_m2(iw) = w1(iw)/w2(iw)*(  w2(iw)/w3(iw)  )**n
       !
       w4(iw) = gammafunc( (b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
       coef_d(iw) = a_m(iw) * w4(iw)/w2(iw)*(  w2(iw)/w3(iw)  )**b_m(iw)
       w5(iw) = gammafunc( (2.0_rp*b_m(iw)+beta_v(iw,2)+nu(iw)+1.0_rp)/mu(iw) )
       w6(iw) = gammafunc( (3.0_rp*b_m(iw)+beta_v(iw,2)+nu(iw)+1.0_rp)/mu(iw) )
       coef_d2v(iw) = a_m(iw) * w6(iw)/w5(iw)* ( w2(iw)/w3(iw) )**b_m(iw)
       coef_md2v(iw)=           w5(iw)/w2(iw)* ( w2(iw)/w3(iw) )**(2.0_rp*b_m(iw)+beta_v(iw,2))
       ! 09/04/14 [Add] T.Mitsui, volume and radar reflectivity
       w7(iw) = gammafunc( (3.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
       coef_d3(iw)  = a_m(iw)**3 * w7(iw)/w2(iw)*(  w2(iw)/w3(iw)  )**(3.0_rp*b_m(iw))
       w8(iw) = gammafunc( (6.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
       coef_d6(iw)  = a_m(iw)**6 * w8(iw)/w2(iw)*(  w2(iw)/w3(iw)  )**(6.0_rp*b_m(iw))
    end do
    !
    coef_deplc = coef_d(i_qc)/a_m(i_qc)
    !-------------------------------------------------------
    ! coefficient of 2nd and 3rd moments for effective radius
    ! for spherical particle
    do iw=i_qc,i_qg
       ! integ r^2 f(x)dx
       w1(iw) = gammafunc( (2.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
       w2(iw) = gammafunc( (nu(iw)+1.0_rp)/mu(iw) )
       w3(iw) = gammafunc( (nu(iw)+2.0_rp)/mu(iw) )
       ! integ r^3 f(x)dx
       w4(iw) = gammafunc( (3.0_rp*b_m(iw)+nu(iw)+1.0_rp)/mu(iw) )
       !
       coef_r2(iw)=w1(iw)/w2(iw)*( w2(iw)/w3(iw) )**(2.0_rp*b_m(iw))
       coef_r3(iw)=w4(iw)/w2(iw)*( w2(iw)/w3(iw) )**(3.0_rp*b_m(iw))
       coef_re(iw)=coef_r3(iw)/coef_r2(iw)
       !
    end do
    !-------------------------------------------------------
    ! coefficient for effective radius of equivalent area and
    ! coefficient for volume of equivalent area
    do iw=i_qc,i_qg
       w1(iw) = gammafunc( (nu(iw)+1.0_rp)/mu(iw) )
       w2(iw) = gammafunc( (nu(iw)+2.0_rp)/mu(iw) )
       w3(iw) = gammafunc( (b_rea2(iw)+nu(iw)+1.0_rp)/mu(iw) )
       w4(iw) = gammafunc( (b_rea3(iw)+nu(iw)+1.0_rp)/mu(iw) )
       !
       coef_rea2(iw) = w3(iw)/w1(iw)*( w1(iw)/w2(iw) )**b_rea2(iw)
       coef_rea3(iw) = w4(iw)/w1(iw)*( w1(iw)/w2(iw) )**b_rea3(iw)
    end do
    !-------------------------------------------------------
    ! coefficient of gamma-distribution
    ! SB06(80)
    do iw=i_qc,i_qg
       w1(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
       w2(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
       coef_a(iw)      = mu(iw)/w1(iw)
!       slope_A(iw)     = w1(iw)
       coef_lambda(iw) = (w1(iw)/w2(iw))**(-mu(iw))
    end do
    !-------------------------------------------------------
    ! coefficient for terminal velocity in sedimentation
    ! SB06(78)
    do ia=1,2
       do iw=i_qc,i_qg
          n = 0
          w1(iw) = gammafunc( (beta_vn(iw,ia) + nu(iw) + 1.0_rp + n)/mu(iw) )
          w2(iw) = gammafunc( (                 nu(iw) + 1.0_rp + n)/mu(iw) )
          w3(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
          w4(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
          ! coefficient of terminal velocity for number
          coef_vt0(iw,ia)=alpha_vn(iw,ia)*w1(iw)/w2(iw)*(w3(iw)/w4(iw))**beta_vn(iw,ia)
          n = 1
          w1(iw) = gammafunc( (beta_v(iw,ia) + nu(iw) + 1.0_rp + n)/mu(iw) )
          w2(iw) = gammafunc( (                nu(iw) + 1.0_rp + n)/mu(iw) )
          ! coefficient of terminal velocity for mass
          coef_vt1(iw,ia)=alpha_v(iw,ia)*w1(iw)/w2(iw)*(w3(iw)/w4(iw))**beta_v(iw,ia)
       end do
    end do
    ! coefficient for weighted diameter used in calculation of terminal velocity
    do iw=i_qc,i_qg
       w1(iw) = gammafunc( (       b_m(iw) + nu(iw) + 1.0_rp)/mu(iw) )
       w2(iw) = gammafunc( (1.0_rp + b_m(iw) + nu(iw) + 1.0_rp)/mu(iw) )
       w3(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
       w4(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
       coef_dave_n(iw) =  (w1(iw)/w3(iw))*(w3(iw)/w4(iw))**      b_m(iw)
       coef_dave_l(iw) =  (w2(iw)/w3(iw))*(w3(iw)/w4(iw))**(1.0_rp+b_m(iw))
    end do
    !-------------------------------------------------------
    !
    ah_vent(i_qc,1:2) = (/1.0000_rp,1.0000_rp/) ! no effect
    ah_vent(i_qr,1:2) = (/1.0000_rp,0.780_rp/)
    ah_vent(i_qi,1:2) = (/1.0000_rp,0.860_rp/)
    ah_vent(i_qs,1:2) = (/1.0000_rp,0.780_rp/)
    ah_vent(i_qg,1:2) = (/1.0000_rp,0.780_rp/)
    bh_vent(i_qc,1:2) = (/0.0000_rp,0.0000_rp/)
    bh_vent(i_qr,1:2) = (/0.108_rp,0.308_rp/)
    bh_vent(i_qi,1:2) = (/0.140_rp,0.280_rp/)
    bh_vent(i_qs,1:2) = (/0.108_rp,0.308_rp/)
    bh_vent(i_qg,1:2) = (/0.108_rp,0.308_rp/)
    !
    do iw=i_qc,i_qg
       n = 0
       if( (nu(iw) + b_m(iw) + n) > eps_gamma  )then
          w1(iw) = gammafunc( (nu(iw) + b_m(iw) + n)/mu(iw) )
          w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
          w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
          ah_vent0(iw,1)= ah_vent(iw,1)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
          ah_vent0(iw,2)= ah_vent(iw,2)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
          flag_vent0(iw)=.true.
       else
          ah_vent0(iw,1)= 1.0_rp
          ah_vent0(iw,2)= 1.0_rp
          flag_vent0(iw)=.false.
       end if
       n = 1
       if( (nu(iw) + b_m(iw) + n) > eps_gamma  )then
          w1(iw) = gammafunc( (nu(iw) + b_m(iw) + n)/mu(iw) )
          w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
          w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
          ah_vent1(iw,1)= ah_vent(iw,1)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
          ah_vent1(iw,2)= ah_vent(iw,2)*(w1(iw)/w2(iw))*(w2(iw)/w3(iw))**(b_m(iw)+n-1.0_rp)
          flag_vent1(iw)=.true.
       else
          ah_vent1(iw,1)= 1.0_rp
          ah_vent1(iw,2)= 1.0_rp
          flag_vent1(iw)=.true.
       end if
    end do
    do iw=i_qc,i_qg
       n = 0
       if( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n) < eps_gamma )then
          flag_vent0(iw)=.false.
       end if
       if(flag_vent0(iw))then
          w1(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n)/mu(iw) )
          w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
          w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
          ! [Add] 11/08/30 T.Mitsui
          w4(iw) = gammafunc( (nu(iw) + 2.0_rp*b_m(iw) + beta_v(iw,1) + n)/mu(iw) )
          bh_vent0(iw,1)=bh_vent(iw,1)*(w4(iw)/w2(iw))*(w2(iw)/w3(iw))**(2.00_rp*b_m(iw)+beta_v(iw,1)+n-1.0_rp)
          w5(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,2) + n)/mu(iw) )
          bh_vent0(iw,2)=bh_vent(iw,2)*(w5(iw)/w2(iw))*(w2(iw)/w3(iw))**(1.5_rp*b_m(iw)+0.5_rp*beta_v(iw,2)+n-1.0_rp)
       else
          bh_vent0(iw,1) = 0.0_rp
          bh_vent0(iw,2) = 0.0_rp
       end if
       !
       n = 1
       if( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n) < eps_gamma )then
          flag_vent1(iw)=.false.
       end if
       if(flag_vent1(iw))then
          w1(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,1) + n)/mu(iw) )
          w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
          w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
          ! [Add] 11/08/30 T.Mitsui
          w4(iw) = gammafunc( (nu(iw) + 2.0_rp*b_m(iw) + beta_v(iw,1) + n)/mu(iw) )
          bh_vent1(iw,1)=bh_vent(iw,1)*(w4(iw)/w2(iw))*(w2(iw)/w3(iw))**(2.00_rp*b_m(iw)+beta_v(iw,1)+n-1.0_rp)
          !
          w5(iw) = gammafunc( (nu(iw) + 1.5_rp*b_m(iw) + 0.5_rp*beta_v(iw,2) + n)/mu(iw) )
          bh_vent1(iw,2)=bh_vent(iw,2)*(w5(iw)/w2(iw))*(w2(iw)/w3(iw))**(1.5_rp*b_m(iw)+0.5_rp*beta_v(iw,2)+n-1.0_rp)
       else
          bh_vent1(iw,1) = 0.0_rp
          bh_vent1(iw,2) = 0.0_rp
       end if
    end do
    !-------------------------------------------------------
    ! coefficient for collision process
    ! stochastic coefficient for collision cross section
    ! sb06 (90) -- self collection
    do iw=i_qc,i_qg
       n = 0
       w1(iw) = gammafunc( (2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
       w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
       w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
       delta_b0(iw) = w1(iw)/w2(iw) &
            *( w2(iw)/w3(iw) )**(2.0_rp*b_rea(iw) + n)
       n = 1
       w1(iw) = gammafunc( (2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
       delta_b1(iw) = w1(iw)/w2(iw) &
            *( w2(iw)/w3(iw) )**(2.0_rp*b_rea(iw) + n)
    end do
    ! stochastic coefficient for collision cross section
    ! sb06(91) -- riming( collide with others )
    do iw=i_qc,i_qg
       n = 0
       w1(iw) = gammafunc( (b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
       w2(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
       w3(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
       w4(iw) = gammafunc( (b_rea(iw) + nu(iw) + 1.0_rp    )/mu(iw) )
       n = 1
       w5(iw) = gammafunc( (b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
    end do
    ! ia > ib ( larger particles "a" catch smaller particles "b" )
    do ia=i_qc,i_qg
       do ib=i_qc,i_qg
          n=0 !
          ! NOTE, collected  particle has a moment of n.
          !       collecting particle has only number(n=0).
          delta_ab0(ia,ib) = 2.0_rp*(w1(ib)/w2(ib))*(w4(ia)/w2(ia)) &
               * ( w2(ib)/w3(ib) )**(b_rea(ib)+n) &
               * ( w2(ia)/w3(ia) )**(b_rea(ia)  )
          n=1 !
          delta_ab1(ia,ib) = 2.0_rp*(w5(ib)/w2(ib))*(w4(ia)/w2(ia)) &
               * ( w2(ib)/w3(ib) )**(b_rea(ib)+n) &
               * ( w2(ia)/w3(ia) )**(b_rea(ia)  )
       end do
    end do
    ! stochastic coefficient for terminal velocity
    ! sb06(92) -- self collection
    ! assuming equivalent area circle.
    do iw=i_qc,i_qg
       n = 0
       w1(iw) = gammafunc( (2.0_rp*beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
       w2(iw) = gammafunc( (                        2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
       w3(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
       w4(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
       theta_b0(iw) = w1(iw)/w2(iw) * ( w3(iw)/w4(iw) )**(2.0_rp*beta_v(iw,2))
       n = 1
       w1(iw) = gammafunc( (2.0_rp*beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
       w2(iw) = gammafunc( (                        2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
       theta_b1(iw) = w1(iw)/w2(iw) * ( w3(iw)/w4(iw) )**(2.0_rp*beta_v(iw,2))
    end do
    !
    ! stochastic coefficient for terminal velocity
    ! sb06(93) -- riming( collide with others )
    do iw=i_qc,i_qg
       n = 0
       w1(iw) = gammafunc( (beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
       w2(iw) = gammafunc( (                   2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
       w3(iw) = gammafunc( (beta_v(iw,2) + 2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp    )/mu(iw) )
       w4(iw) = gammafunc( (                   2.0_rp*b_rea(iw) + nu(iw) + 1.0_rp    )/mu(iw) )
       !
       w5(iw) = gammafunc( (nu(iw) + 1.0_rp)/mu(iw) )
       w6(iw) = gammafunc( (nu(iw) + 2.0_rp)/mu(iw) )
       n = 1
       w7(iw) = gammafunc( (beta_v(iw,2) + b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
       w8(iw) = gammafunc( (                   b_rea(iw) + nu(iw) + 1.0_rp + n)/mu(iw) )
    end do
    ! ia > ib ( larger particles "a" catch smaller particles "b" )
    do ia=i_qc,i_qg
       do ib=i_qc,i_qg
          theta_ab0(ia,ib) = 2.0_rp * (w1(ib)/w2(ib))*(w3(ia)/w4(ia)) &
               * (w5(ia)/w6(ia))**beta_v(ia,2) &
               * (w5(ib)/w6(ib))**beta_v(ib,2)
          theta_ab1(ia,ib) = 2.0_rp * (w7(ib)/w8(ib))*(w3(ia)/w4(ia)) &
               * (w5(ia)/w6(ia))**beta_v(ia,2) &
               * (w5(ib)/w6(ib))**beta_v(ib,2)
       end do
    end do

    deallocate(w1,w2,w3,w4,w5,w6,w7,w8)

    if( io_l ) write(io_fid_log,'(100a16)')     "LABEL       ",wlabel(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "capacity    ",cap(:) ! [Add] 11/08/30 T.Mitsui
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_m2     ",coef_m2(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_d      ",coef_d(:)
    !
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_d3     ",coef_d3(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_d6     ",coef_d6(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_d2v    ",coef_d2v(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_md2v   ",coef_md2v(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "a_d2vt      ",a_d2vt(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "b_d2vt      ",b_d2vt(:)
    !
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_r2     ",coef_r2(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_r3     ",coef_r3(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_re     ",coef_re(:)
    !
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "a_area      ",a_area(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "b_area      ",b_area(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "ax_area     ",ax_area(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "bx_area     ",bx_area(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "a_rea       ",a_rea(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "b_rea       ",b_rea(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "a_rea3      ",a_rea3(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "b_rea3      ",b_rea3(:)
    !
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_rea2   ",coef_rea2(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_rea3   ",coef_rea3(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_vt0    ",coef_vt0(:,1)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_vt1    ",coef_vt1(:,1)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_A      ",coef_a(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "coef_lambda ",coef_lambda(:)

    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "ah_vent0 sml",ah_vent0(:,1)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "ah_vent0 lrg",ah_vent0(:,2)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "ah_vent1 sml",ah_vent1(:,1)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "ah_vent1 lrg",ah_vent1(:,2)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "bh_vent0 sml",bh_vent0(:,1)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "bh_vent0 lrg",bh_vent0(:,2)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "bh_vent1 sml",bh_vent1(:,1)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "bh_vent1 lrg",bh_vent1(:,2)

    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "delta_b0    ",delta_b0(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "delta_b1    ",delta_b1(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "theta_b0    ",theta_b0(:)
    if( io_l ) write(io_fid_log,'(a,100ES16.6)') "theta_b1    ",theta_b1(:)

    do ia=qqs,qqe
       if( io_l ) write(io_fid_log,'(a,a10,a,100ES16.6)') "delta0(a,b)=(",trim(wlabel(ia)),",b)=",(delta_ab0(ia,ib),ib=qqs,qqe)
    enddo
    do ia=qqs,qqe
       if( io_l ) write(io_fid_log,'(a,a10,a,100ES16.6)') "delta1(a,b)=(",trim(wlabel(ia)),",b)=",(delta_ab1(ia,ib),ib=qqs,qqe)
    enddo
    do ia=qqs,qqe
       if( io_l ) write(io_fid_log,'(a,a10,a,100ES16.6)') "theta0(a,b)=(",trim(wlabel(ia)),",b)=",(theta_ab0(ia,ib),ib=qqs,qqe)
    enddo
    do ia=qqs,qqe
       if( io_l ) write(io_fid_log,'(a,a10,a,100ES16.6)') "theta1(a,b)=(",trim(wlabel(ia)),",b)=",(theta_ab1(ia,ib),ib=qqs,qqe)
    enddo

    return
  end subroutine mp_sn14_init
  !-----------------------------------------------------------------------------
  subroutine mp_sn14 ( &
       DENS,      &
       MOMZ,      &
       MOMX,      &
       MOMY,      &
       RHOT,      &
       QTRC,      &
       CCN,       &
       EVAPORATE, &
       SFLX_rain, &
       SFLX_snow  )
    use scale_time, only: &
       dt => time_dtsec_atmos_phy_mp
    use scale_grid, only: &
       z    => grid_cz, &
       dz   => grid_cdz
    use scale_atmos_phy_mp_common, only: &
         mp_precipitation => atmos_phy_mp_precipitation
    use scale_atmos_thermodyn, only: &
       aq_cv, &
       aq_cp
    use scale_atmos_saturation, only: &
       moist_psat_liq      => atmos_saturation_psat_liq,   &
       moist_psat_ice      => atmos_saturation_psat_ice
    implicit none

    real(RP), intent(inout) :: DENS(ka,ia,ja)
    real(RP), intent(inout) :: MOMZ(ka,ia,ja)
    real(RP), intent(inout) :: MOMX(ka,ia,ja)
    real(RP), intent(inout) :: MOMY(ka,ia,ja)
    real(RP), intent(inout) :: RHOT(ka,ia,ja)
    real(RP), intent(inout) :: QTRC(ka,ia,ja,qa)
    real(RP), intent(in)    :: CCN(ka,ia,ja)
    real(RP), intent(out)   :: EVAPORATE(ka,ia,ja)
    real(RP), intent(out)   :: SFLX_rain(ia,ja)
    real(RP), intent(out)   :: SFLX_snow(ia,ja)

    !
    ! primary variables
    !
    real(RP) :: rhoq(qa,ka,ia,ja)
    real(RP) :: pott   (ka,ia,ja)
    !
    ! diagnostic variables
    !
    real(RP) :: qdry(ka,ia,ja)
    real(RP) :: temp(ka,ia,ja)
    real(RP) :: pres(ka,ia,ja)
    real(RP) :: rhoe(ka,ia,ja)
    real(RP) :: velz(ka,ia,ja)
    !
    real(RP) :: rhoq2(qa,ka,ia,ja)
    !
    real(RP) :: xq(5,ka,ia,ja)
    !
    real(RP) :: dq_xa(5,ka,ia,ja)
    real(RP) :: vt_xa(5,2,ka,ia,ja) ! terminal velocity

    real(RP) :: wtemp(ka,ia,ja) ! filtered temperature
    real(RP) :: esw(ka,ia,ja) ! saturated vapor pressure(water)
    real(RP) :: esi(ka,ia,ja) ! saturated vapor pressure(ice)
    !
    real(RP) :: rho_fac
    real(RP) :: rho_fac_q(5,ka,ia,ja) ! factor for tracers, 1:cloud, 2:rain, 3:ice, 4: snow, 5:graupel
    real(RP) :: cva(ka,ia,ja)    !
    real(RP) :: cpa(ka,ia,ja)       ! [Add] 09/08/18 T.Mitsui
    !
    real(RP) :: drhogqv               ! d (rho*qv*gsgam2)
    real(RP) :: drhogqc, drhognc      !        qc, nc
    real(RP) :: drhogqr, drhognr      !        qr, nr
    real(RP) :: drhogqi, drhogni      !        qi, ni
    real(RP) :: drhogqs, drhogns      !        qs, ns
    real(RP) :: drhogqg, drhogng      !        qg, ng

    ! production rate
    real(RP) :: PQ(pq_max,ka,ia,ja)

    real(RP) :: wrm_dqc, wrm_dnc
    real(RP) :: wrm_dqr, wrm_dnr

    ! production rate of mixed-phase collection process
    real(RP) :: Pac(pac_max,ka,ia,ja)

    real(RP) :: gc_dqc, gc_dnc
    real(RP) :: sc_dqc, sc_dnc
    real(RP) :: ic_dqc, ic_dnc
    real(RP) :: rg_dqg, rg_dng
    real(RP) :: rg_dqr, rg_dnr
    real(RP) :: rs_dqr, rs_dnr, rs_dqs, rs_dns
    real(RP) :: ri_dqr, ri_dnr
    real(RP) :: ri_dqi, ri_dni
    real(RP) :: ii_dqi, ii_dni
    real(RP) :: is_dqi, is_dni, ss_dns
    real(RP) :: gs_dqs, gs_dns, gg_dng
    ! mixed-phase collection process total plus(clp_), total minus(clm_)
    real(RP) :: clp_dqc, clp_dnc, clm_dqc, clm_dnc
    real(RP) :: clp_dqr, clp_dnr, clm_dqr, clm_dnr
    real(RP) :: clp_dqi, clp_dni, clm_dqi, clm_dni
    real(RP) :: clp_dqs, clp_dns, clm_dqs, clm_dns
    real(RP) :: clp_dqg, clp_dng, clm_dqg, clm_dng
    real(RP) :: fac1, fac3, fac4, fac6, fac7, fac9, fac10
    ! production rate of partial conversion(ice, snow => graupel)
    real(RP) :: pco_dqi, pco_dni
    real(RP) :: pco_dqs, pco_dns
    real(RP) :: pco_dqg, pco_dng
    ! production rate of enhanced melting due to
    real(RP) :: eml_dqc, eml_dnc
    real(RP) :: eml_dqr, eml_dnr
    real(RP) :: eml_dqi, eml_dni
    real(RP) :: eml_dqs, eml_dns
    real(RP) :: eml_dqg, eml_dng
    ! production rate of ice multiplication by splintering
    real(RP) :: spl_dqi, spl_dni
    real(RP) :: spl_dqg, spl_dqs

    real(RP) :: rrho(ka,ia,ja)

    !-----------------------------------------------
    ! work for explicit supersaturation modeling
    !-----------------------------------------------
    real(RP) :: dTdt_equiv_d(ka,ia,ja)    !
    !--------------------------------------------------
    !
    ! variables for output
    !
    !--------------------------------------------------
    ! work for column production term
    real(RP) :: sl_PLCdep(ia,ja)
    real(RP) :: sl_PLRdep(ia,ja), sl_PNRdep(ia,ja) !
    !--------------------------------------------------
    real(RP) :: qke_d(ka,ia,ja)

    real(RP), parameter :: eps       = 1.e-19_rp
    real(RP), parameter :: eps_qv    = 1.e-19_rp
    real(RP), parameter :: eps_rhoge = 1.e-19_rp
    real(RP), parameter :: eps_rhog  = 1.e-19_rp
    integer :: ntdiv

    real(RP) :: Rmoist

    real(RP) :: velw(ka,ia,ja,qa)
    real(RP) :: FLX_rain (ka,ia,ja)
    real(RP) :: FLX_snow (ka,ia,ja)
    real(RP) :: FLX_tot  (ka,ia,ja)
    real(RP) :: wflux_rain(ka,ia,ja)
    real(RP) :: wflux_snow(ka,ia,ja)
    integer :: step

    real(RP) :: sw

    integer :: k, i, j, iq
    !---------------------------------------------------------------------------

    !============================================================================
    !
    !--  Each process is integrated sequentially.
    !    1. Nucleation and filter
    !    2. Phase change
    !    3. Collection
    !    4. Saturation adjustment( only for qc, nc )
    !    5. Sedimentation
    !    6. filter( keep non-negative value for radiation scheme )
    !    0. calculation of optical moments
    !
    !============================================================================
    !----------------------------------------------------------------------------
    !
    ! 1.Nucleation of cloud water and cloud ice
    !
    !----------------------------------------------------------------------------
    call prof_rapstart('MP_Preprocess', 3)

    do j = js, je
    do i = is, ie
    do k = ks, ke
       do iq = 1, qa
          rhoq(iq,k,i,j) = dens(k,i,j) * qtrc(k,i,j,iq)
       enddo
       rhoq2(i_qv,k,i,j) = dens(k,i,j)*qtrc(k,i,j,i_qv)
       rhoq2(i_ni,k,i,j) = max( 0.0_rp, dens(k,i,j)*qtrc(k,i,j,i_ni) )
       rhoq2(i_nc,k,i,j) = max( 0.0_rp, dens(k,i,j)*qtrc(k,i,j,i_nc) )
    enddo
    enddo
    enddo

    do j = js, je
    do i = is, ie
       velz(ks-1,i,j) = 0.0_rp
       do k = ks, ke-1
          velz(k,i,j) = momz(k,i,j) / ( dens(k,i,j) + dens(k+1,i,j) ) * 2.0_rp
       enddo
       velz(ke,i,j) = 0.0_rp
    end do
    end do

    do j = js, je
    do i = is, ie
    do k = ks, ke
       rrho(k,i,j) = 1.0_rp / dens(k,i,j)
       pott(k,i,j) = rhot(k,i,j) * rrho(k,i,j)
       calc_qdry( qdry(k,i,j), qtrc, k, i, j, iq )
       calc_cv( cva(k,i,j), qdry(k,i,j), qtrc, k, i, j, iq, cvdry, aq_cv )
       calc_r( rmoist, qtrc(k,i,j,i_qv), qdry(k,i,j), rdry, rvap )
       cpa(k,i,j) = cva(k,i,j) + rmoist
       calc_pre( pres(k,i,j), dens(k,i,j), pott(k,i,j), rmoist, cpa(k,i,j), p00 )
       temp(k,i,j) = pres(k,i,j) / ( dens(k,i,j) * rmoist )
       rhoe(k,i,j) = dens(k,i,j) * temp(k,i,j) * cva(k,i,j)
       wtemp(k,i,j) = max(temp(k,i,j), tem_min)
    enddo
    enddo
    enddo

    if( opt_debug_tem ) call debug_tem_kij( 1, temp(:,:,:), dens(:,:,:), pres(:,:,:), qtrc(:,:,:,i_qv) )

    do j = js, je
    do i = is, ie
    do k = ks, ke
       rho_fac              = rho_0 / max(dens(k,i,j),rho_min)
       rho_fac_q(i_c,k,i,j) = rho_fac**gamma_v(i_qc)
       rho_fac_q(i_r,k,i,j) = rho_fac**gamma_v(i_qr)
       rho_fac_q(i_i,k,i,j) = (pres(k,i,j)/pre0_vt)**a_pre0_vt * (temp(k,i,j)/tem0_vt)**a_tem0_vt
       rho_fac_q(i_s,k,i,j) = rho_fac_q(i_i,k,i,j)
       rho_fac_q(i_g,k,i,j) = rho_fac_q(i_i,k,i,j)
    enddo
    enddo
    enddo

!OCL XFILL
    do j = js, je
    do i = is, ie
       sl_plcdep(i,j) = 0.0_rp
       sl_plrdep(i,j) = 0.0_rp
       sl_pnrdep(i,j) = 0.0_rp
    end do
    end do

!OCL XFILL
    do j = js, je
    do i = is, ie
    do k = ks, ke
       qke_d(k,i,j) = 0.0_rp ! 2*TKE
    enddo
    enddo
    enddo

!OCL XFILL
    do j = js, je
    do i = is, ie
    do k = ks, ke
       dtdt_equiv_d(k,i,j) = 0.0_rp
    enddo
    enddo
    enddo

    call prof_rapend  ('MP_Preprocess', 3)

    call prof_rapstart('MP_Nucleation', 3)

    call nucleation_kij(    &
         z, velz,           & ! in
         dens, wtemp, pres, & ! in
         rhoq2,             & ! (in)
         pq,                & ! out
         cpa,               & ! in
         dtdt_equiv_d,      & ! in
         qke_d,             & ! in
         ccn,               & ! in
         dt                 ) ! in

    do j = js, je
    do i = is, ie
    do k = ks, ke
       ! nucleation
       drhogqc = dt * pq(i_lcccn,k,i,j)
       drhognc = dt * pq(i_ncccn,k,i,j)
       drhogqi = dt * pq(i_liccn,k,i,j)
       drhogni = dt * pq(i_niccn,k,i,j)
       drhogqv = max( -rhoq(i_qv,k,i,j), -drhogqc-drhogqi )
       fac1    = drhogqv / min( -drhogqc-drhogqi, -eps ) ! limiting coefficient

       rhoq(i_qv,k,i,j) = rhoq(i_qv,k,i,j) + drhogqv
       rhoq(i_qc,k,i,j) = max(0.0_rp, rhoq(i_qc,k,i,j) + drhogqc*fac1)
       rhoq(i_qi,k,i,j) = max(0.0_rp, rhoq(i_qi,k,i,j) + drhogqi*fac1)
       rhoq(i_nc,k,i,j) = max(0.0_rp, rhoq(i_nc,k,i,j) + drhognc)
       rhoq(i_ni,k,i,j) = max(0.0_rp, rhoq(i_ni,k,i,j) + drhogni)

       ! cloud number concentration filter
       rhoq(i_nc,k,i,j) = min( rhoq(i_nc,k,i,j), nc_uplim_d(1,i,j) )

       rhoe(k,i,j) = rhoe(k,i,j) - lhv * drhogqv + lhf * drhogqi*fac1

       qtrc(k,i,j,i_qv) = rhoq(i_qv,k,i,j) * rrho(k,i,j)
       qtrc(k,i,j,i_qc) = rhoq(i_qc,k,i,j) * rrho(k,i,j)
       qtrc(k,i,j,i_qi) = rhoq(i_qi,k,i,j) * rrho(k,i,j)
       qtrc(k,i,j,i_nc) = rhoq(i_nc,k,i,j) * rrho(k,i,j)
       qtrc(k,i,j,i_ni) = rhoq(i_ni,k,i,j) * rrho(k,i,j)

       calc_qdry( qdry(k,i,j), qtrc, k, i, j, iq )
       calc_cv( cva(k,i,j), qdry(k,i,j), qtrc, k, i, j, iq, cvdry, aq_cv )
       calc_r( rmoist, qtrc(k,i,j,i_qv), qdry(k,i,j), rdry, rvap )
       temp(k,i,j) = rhoe(k,i,j) / ( dens(k,i,j) * cva(k,i,j) )
       pres(k,i,j) = dens(k,i,j) * rmoist * temp(k,i,j)
       wtemp(k,i,j) = max( temp(k,i,j), tem_min )
    enddo
    enddo
    enddo

!    if( opt_debug )     call debugreport_nucleation
    if( opt_debug_tem ) call debug_tem_kij( 2, temp(:,:,:), dens(:,:,:), pres(:,:,:), qtrc(:,:,:,i_qv) )


    call prof_rapend  ('MP_Nucleation', 3)
    !----------------------------------------------------------------------------
    !
    ! 2.Phase change: Freezing, Melting, Vapor deposition
    !
    !----------------------------------------------------------------------------
    call prof_rapstart('MP_Phase_change', 3)

    do j = js, je
    do i = is, ie
    do k = ks, ke
       rhoq2(i_qr,k,i,j)     = rhoq(i_qr,k,i,j)
       rhoq2(i_nr,k,i,j)     = rhoq(i_nr,k,i,j)
       xq(i_r,k,i,j)     = max(xr_min,  min(xr_max, rhoq2(i_qr,k,i,j)/(rhoq2(i_nr,k,i,j)+nr_min) ))

       dq_xa(i_r,k,i,j)  = a_m(i_qr)*xq(i_r,k,i,j)**b_m(i_qr)
       vt_xa(i_r,1,k,i,j) = alpha_v(i_qr,1)*(xq(i_r,k,i,j)**beta_v(i_qr,1))*rho_fac_q(i_r,k,i,j)
       vt_xa(i_r,2,k,i,j) = vt_xa(i_r,1,k,i,j)

       !! Following values shoud be already filtered to be non-zero before sbroutine was called.
       ! Mass concentration [kg/m3]
       rhoq2(i_qv,k,i,j) = rhoq(i_qv,k,i,j)
       rhoq2(i_qc,k,i,j) = rhoq(i_qc,k,i,j)
       rhoq2(i_qi,k,i,j) = rhoq(i_qi,k,i,j)
       rhoq2(i_qs,k,i,j) = rhoq(i_qs,k,i,j)
       rhoq2(i_qg,k,i,j) = rhoq(i_qg,k,i,j)
       ! Number concentration[/m3] (should be filtered to avoid zero division.)
       rhoq2(i_nc,k,i,j) = rhoq(i_nc,k,i,j)
       rhoq2(i_ni,k,i,j) = rhoq(i_ni,k,i,j)
       rhoq2(i_ns,k,i,j) = rhoq(i_ns,k,i,j)
       rhoq2(i_ng,k,i,j) = rhoq(i_ng,k,i,j)

       ! Mass of mean particle [kg]
       ! SB06(94)
       !
       xq(i_c,k,i,j)     = min(xc_max, max(xc_min, rhoq2(i_qc,k,i,j)/(rhoq2(i_nc,k,i,j)+nc_min) ))
       xq(i_i,k,i,j)     = min(xi_max, max(xi_min, rhoq2(i_qi,k,i,j)/(rhoq2(i_ni,k,i,j)+ni_min) ))
       xq(i_s,k,i,j)     = min(xs_max, max(xs_min, rhoq2(i_qs,k,i,j)/(rhoq2(i_ns,k,i,j)+ns_min) ))
       xq(i_g,k,i,j)     = min(xg_max, max(xg_min, rhoq2(i_qg,k,i,j)/(rhoq2(i_ng,k,i,j)+ng_min) ))
       ! diamter of average mass
       ! SB06(32)
       dq_xa(i_c,k,i,j)  = a_m(i_qc)*xq(i_c,k,i,j)**b_m(i_qc)
       dq_xa(i_i,k,i,j)  = a_m(i_qi)*xq(i_i,k,i,j)**b_m(i_qi)
       dq_xa(i_s,k,i,j)  = a_m(i_qs)*xq(i_s,k,i,j)**b_m(i_qs)
       dq_xa(i_g,k,i,j)  = a_m(i_qg)*xq(i_g,k,i,j)**b_m(i_qg)

       ! terminal velocity of average mass
       vt_xa(i_c,1,k,i,j) = alpha_v(i_qc,1)*(xq(i_c,k,i,j)**beta_v(i_qc,1))*rho_fac_q(i_c,k,i,j)
       vt_xa(i_i,1,k,i,j) = alpha_v(i_qi,1)*(xq(i_i,k,i,j)**beta_v(i_qi,1))*rho_fac_q(i_i,k,i,j)
       vt_xa(i_s,1,k,i,j) = alpha_v(i_qs,1)*(xq(i_s,k,i,j)**beta_v(i_qs,1))*rho_fac_q(i_s,k,i,j)
       vt_xa(i_g,1,k,i,j) = alpha_v(i_qg,1)*(xq(i_g,k,i,j)**beta_v(i_qg,1))*rho_fac_q(i_g,k,i,j)
       vt_xa(i_c,2,k,i,j) = alpha_v(i_qc,2)*(xq(i_c,k,i,j)**beta_v(i_qc,2))*rho_fac_q(i_c,k,i,j)
       vt_xa(i_i,2,k,i,j) = alpha_v(i_qi,2)*(xq(i_i,k,i,j)**beta_v(i_qi,2))*rho_fac_q(i_i,k,i,j)
       vt_xa(i_s,2,k,i,j) = alpha_v(i_qs,2)*(xq(i_s,k,i,j)**beta_v(i_qs,2))*rho_fac_q(i_s,k,i,j)
       vt_xa(i_g,2,k,i,j) = alpha_v(i_qg,2)*(xq(i_g,k,i,j)**beta_v(i_qg,2))*rho_fac_q(i_g,k,i,j)

    end do
    end do
    end do

    call moist_psat_liq( esw, wtemp )
    call moist_psat_ice( esi, wtemp )

    call freezing_water_kij( &
         dt,             & ! in
         pq,             & ! out
         rhoq2, xq, temp ) ! in

    call dep_vapor_melt_ice_kij( &
         pq,                 & ! out
         dens, wtemp, pres, qdry, & ! in
         rhoq2,               & ! in
         esw, esi,           & ! in
         xq,                 & ! in
         vt_xa,              & ! in
         dq_xa               ) ! in

    !
    ! update subroutine
    !
    call update_by_phase_change_kij( &
         ntdiv, ntmax_phase_change,  & ! in
         dt,                         & ! in
         gsgam2_d,                   & ! in
         z,                          & ! in
         dz,                         & ! in
         velz,                       & ! in
         dtdt_equiv_d,               & ! in
         dens,                       & ! in
         rhoe,                       & ! inout
         rhoq, qtrc,                 & ! inout
         temp, pres,                 & ! inout
         cva,                        & ! out
         esw, esi, rhoq2,            & ! in
         pq,                         & ! inout
         evaporate,                  & ! out
         sl_plcdep,                  & ! inout
         sl_plrdep, sl_pnrdep        ) ! inout

!    if( opt_debug )     call debugreport_phasechange
    if( opt_debug_tem ) call debug_tem_kij( 3, temp(:,:,:), dens(:,:,:), pres(:,:,:), qtrc(:,:,:,i_qv) )

    call prof_rapend  ('MP_Phase_change', 3)

    !---------------------------------------------------------------------------
    !
    ! 3.Collection process
    !
    !---------------------------------------------------------------------------
    call prof_rapstart('MP_Collection', 3)

    ! parameter setting
    do j = js, je
    do i = is, ie
    do k = ks, ke
       ! Mass concentration [kg/m3]
       rhoq2(i_qv,k,i,j) = rhoq(i_qv,k,i,j)
       rhoq2(i_qc,k,i,j) = rhoq(i_qc,k,i,j)
       rhoq2(i_qr,k,i,j) = rhoq(i_qr,k,i,j)
       rhoq2(i_qi,k,i,j) = rhoq(i_qi,k,i,j)
       rhoq2(i_qs,k,i,j) = rhoq(i_qs,k,i,j)
       rhoq2(i_qg,k,i,j) = rhoq(i_qg,k,i,j)
       ! Number concentration[/m3]
       rhoq2(i_nc,k,i,j) = rhoq(i_nc,k,i,j)
       rhoq2(i_nr,k,i,j) = rhoq(i_nr,k,i,j)
       rhoq2(i_ni,k,i,j) = rhoq(i_ni,k,i,j)
       rhoq2(i_ns,k,i,j) = rhoq(i_ns,k,i,j)
       rhoq2(i_ng,k,i,j) = rhoq(i_ng,k,i,j)

       ! Mass of mean particle [kg]
       xq(i_c,k,i,j) = min(xc_max, max(xc_min, rhoq2(i_qc,k,i,j)/(rhoq2(i_nc,k,i,j)+nc_min) ) )
       xq(i_r,k,i,j) = min(xr_max, max(xr_min, rhoq2(i_qr,k,i,j)/(rhoq2(i_nr,k,i,j)+nr_min) ) )
       xq(i_i,k,i,j) = min(xi_max, max(xi_min, rhoq2(i_qi,k,i,j)/(rhoq2(i_ni,k,i,j)+ni_min) ) )
       xq(i_s,k,i,j) = min(xs_max, max(xs_min, rhoq2(i_qs,k,i,j)/(rhoq2(i_ns,k,i,j)+ns_min) ) )
       xq(i_g,k,i,j) = min(xg_max, max(xg_min, rhoq2(i_qg,k,i,j)/(rhoq2(i_ng,k,i,j)+ng_min) ) )

       ! effective cross section is assume as area equivalent circle
       dq_xa(i_c,k,i,j) = 2.0_rp*a_rea(i_qc)*xq(i_c,k,i,j)**b_rea(i_qc)
       dq_xa(i_r,k,i,j) = 2.0_rp*a_rea(i_qr)*xq(i_r,k,i,j)**b_rea(i_qr)
       dq_xa(i_i,k,i,j) = 2.0_rp*a_rea(i_qi)*xq(i_i,k,i,j)**b_rea(i_qi)
       dq_xa(i_s,k,i,j) = 2.0_rp*a_rea(i_qs)*xq(i_s,k,i,j)**b_rea(i_qs)
       dq_xa(i_g,k,i,j) = 2.0_rp*a_rea(i_qg)*xq(i_g,k,i,j)**b_rea(i_qg)

       ! terminal velocity of average mass
       ! SB06(33)
       vt_xa(i_c,2,k,i,j) = alpha_v(i_qc,2)*(xq(i_c,k,i,j)**beta_v(i_qc,2))*rho_fac_q(i_c,k,i,j)
       vt_xa(i_r,2,k,i,j) = alpha_v(i_qr,2)*(xq(i_r,k,i,j)**beta_v(i_qr,2))*rho_fac_q(i_r,k,i,j)
       vt_xa(i_i,2,k,i,j) = alpha_v(i_qi,2)*(xq(i_i,k,i,j)**beta_v(i_qi,2))*rho_fac_q(i_i,k,i,j)
       vt_xa(i_s,2,k,i,j) = alpha_v(i_qs,2)*(xq(i_s,k,i,j)**beta_v(i_qs,2))*rho_fac_q(i_s,k,i,j)
       vt_xa(i_g,2,k,i,j) = alpha_v(i_qg,2)*(xq(i_g,k,i,j)**beta_v(i_qg,2))*rho_fac_q(i_g,k,i,j)
    enddo
    enddo
    enddo

    ! Auto-conversion, Accretion, Self-collection, Break-up
    ! [Mod] T.Seiki
    if ( mp_doautoconversion ) then
       call aut_acc_slc_brk_kij(  &
            pq, &
            rhoq2, xq, dq_xa, &
            dens               )
    else
!OCL XFILL
       do j = js, je
       do i = is, ie
       do k = ks, ke
          pq(i_lcaut,k,i,j) = 0.0_rp
          pq(i_ncaut,k,i,j) = 0.0_rp
          pq(i_nraut,k,i,j) = 0.0_rp
          pq(i_lcacc,k,i,j) = 0.0_rp
          pq(i_ncacc,k,i,j) = 0.0_rp
          pq(i_nrslc,k,i,j) = 0.0_rp
          pq(i_nrbrk,k,i,j) = 0.0_rp
       end do
       end do
       end do
    endif

    call mixed_phase_collection_kij( &
         ! collection process
         pac, pq,                    & ! out
         temp, rhoq2,                & ! in
         xq, dq_xa, vt_xa            ) ! in
!         DENS(:,:,:),                ) ! in

    call ice_multiplication_kij( &
         pq,                     & ! out
         pac,                    & ! in
         temp, rhoq2, xq         ) ! in

    !
    ! update
    ! rhogq = l*gsgam
    !
    profile_start("sn14_update_rhoq")
    do j = js, je
    do i = is, ie
    do k = ks, ke

       ! warm collection process
       wrm_dqc = max( dt*( pq(i_lcaut,k,i,j)+pq(i_lcacc,k,i,j) ), -rhoq2(i_qc,k,i,j)  )
       wrm_dnc = max( dt*( pq(i_ncaut,k,i,j)+pq(i_ncacc,k,i,j) ), -rhoq2(i_nc,k,i,j)  )
       wrm_dnr = max( dt*( pq(i_nraut,k,i,j)+pq(i_nrslc,k,i,j)+pq(i_nrbrk,k,i,j) ), -rhoq2(i_nr,k,i,j) )
       wrm_dqr = -wrm_dqc
       ! mixed phase collection
       ! Pxxacyy2zz xx and yy decrease and zz increase .
       !
       ! At first fixer is applied to decreasing particles.
       ! order of fixer: graupel-cloud, snow-cloud, ice-cloud, graupel-rain, snow-rain, ice-rain,
       !                 snow-ice,  ice-ice, graupel-snow, snow-snow
       ! cloud mass decrease
       gc_dqc  = max( dt*pac(i_lgaclc2lg,k,i,j)  , min(0.0_rp, -rhoq2(i_qc,k,i,j)-wrm_dqc               )) ! => dqg
       sc_dqc  = max( dt*pac(i_lsaclc2ls,k,i,j)  , min(0.0_rp, -rhoq2(i_qc,k,i,j)-wrm_dqc-gc_dqc        )) ! => dqs
       ic_dqc  = max( dt*pac(i_liaclc2li,k,i,j)  , min(0.0_rp, -rhoq2(i_qc,k,i,j)-wrm_dqc-gc_dqc-sc_dqc )) ! => dqi
       ! cloud num. decrease
       gc_dnc  = max( dt*pac(i_ngacnc2ng,k,i,j)  , min(0.0_rp, -rhoq2(i_nc,k,i,j)-wrm_dnc               )) ! => dnc:minus
       sc_dnc  = max( dt*pac(i_nsacnc2ns,k,i,j)  , min(0.0_rp, -rhoq2(i_nc,k,i,j)-wrm_dnc-gc_dnc        )) ! => dnc:minus
       ic_dnc  = max( dt*pac(i_niacnc2ni,k,i,j)  , min(0.0_rp, -rhoq2(i_nc,k,i,j)-wrm_dnc-gc_dnc-sc_dnc )) ! => dnc:minus

       ! rain mass decrease ( tem < 273.15K)
       sw = sign(0.5_rp, t00-temp(k,i,j)) + 0.5_rp ! if( temp(k,i,j) <= T00 )then sw=1, else sw=0
       rg_dqr  = max( dt*pac(i_lraclg2lg,  k,i,j), min(0.0_rp, -rhoq2(i_qr,k,i,j)-wrm_dqr               )) * sw
       rg_dqg  = max( dt*pac(i_lraclg2lg,  k,i,j), min(0.0_rp, -rhoq2(i_qg,k,i,j)                       )) * ( 1.0_rp - sw )
       rs_dqr  = max( dt*pac(i_lracls2lg_r,k,i,j), min(0.0_rp, -rhoq2(i_qr,k,i,j)-wrm_dqr-rg_dqr        )) * sw
       ri_dqr  = max( dt*pac(i_lracli2lg_r,k,i,j), min(0.0_rp, -rhoq2(i_qr,k,i,j)-wrm_dqr-rg_dqr-rs_dqr )) * sw
       ! rain num. decrease
       rg_dnr  = max( dt*pac(i_nracng2ng,  k,i,j), min(0.0_rp, -rhoq2(i_nr,k,i,j)-wrm_dnr               )) * sw
       rg_dng  = max( dt*pac(i_nracng2ng,  k,i,j), min(0.0_rp, -rhoq2(i_ng,k,i,j)                       )) * ( 1.0_rp - sw )
       rs_dnr  = max( dt*pac(i_nracns2ng_r,k,i,j), min(0.0_rp, -rhoq2(i_nr,k,i,j)-wrm_dnr-rg_dnr        )) * sw
       ri_dnr  = max( dt*pac(i_nracni2ng_r,k,i,j), min(0.0_rp, -rhoq2(i_nr,k,i,j)-wrm_dnr-rg_dnr-rs_dnr )) * sw

       ! ice mass decrease
       fac1    = (ri_dqr-eps)/ (dt*pac(i_lracli2lg_r,k,i,j)-eps) ! suppress factor by filter of rain
       ri_dqi  = max( dt*pac(i_lracli2lg_i,k,i,j)*fac1, min(0.0_rp, -rhoq2(i_qi,k,i,j)+ic_dqc               )) ! => dqg
       ii_dqi  = max( dt*pac(i_liacli2ls,k,i,j)       , min(0.0_rp, -rhoq2(i_qi,k,i,j)+ic_dqc-ri_dqi        )) ! => dqs
       is_dqi  = max( dt*pac(i_liacls2ls,k,i,j)       , min(0.0_rp, -rhoq2(i_qi,k,i,j)+ic_dqc-ri_dqi-ii_dqi )) ! => dqs
       ! ice num. decrease
       fac4    = (ri_dnr-eps)/ (dt*pac(i_nracni2ng_r,k,i,j)-eps) ! suppress factor by filter of rain
       ri_dni  = max( dt*pac(i_nracni2ng_i,k,i,j)*fac4, min(0.0_rp, -rhoq2(i_ni,k,i,j)               )) ! => dni:minus
       ii_dni  = max( dt*pac(i_niacni2ns,k,i,j)       , min(0.0_rp, -rhoq2(i_ni,k,i,j)-ri_dni        )) ! => dni:minus,dns:plus(*0.5)
       is_dni  = max( dt*pac(i_niacns2ns,k,i,j)       , min(0.0_rp, -rhoq2(i_ni,k,i,j)-ri_dni-ii_dni )) ! => dni:minus,dns:plus
       ! snow mass decrease
       fac3    = (rs_dqr-eps)/(dt*pac(i_lracls2lg_r,k,i,j)-eps) ! suppress factor by filter of rain
       rs_dqs  = max( dt*pac(i_lracls2lg_s,k,i,j)*fac3, min(0.0_rp, -rhoq2(i_qs,k,i,j)+sc_dqc+ii_dqi+is_dqi        )) ! => dqg
       gs_dqs  = max( dt*pac(i_lgacls2lg,k,i,j)       , min(0.0_rp, -rhoq2(i_qs,k,i,j)+sc_dqc+ii_dqi+is_dqi-rs_dqs )) ! => dqg
       ! snow num. decrease
       fac6    = (rs_dnr-eps)/(dt*pac(i_nracns2ng_r,k,i,j)-eps) ! suppress factor by filter of rain
!       fac7    = (is_dni-eps)/(dt*Pac(I_NIacNS2NS,  k,i,j)-eps) ! suppress factor by filter of ice
       rs_dns  = max( dt*pac(i_nracns2ng_s,k,i,j)*fac6, min(0.0_rp, -rhoq2(i_ns,k,i,j)+0.50_rp*ii_dni+is_dni       )) ! => dns:minus
       gs_dns  = max( dt*pac(i_ngacns2ng,k,i,j)       , min(0.0_rp, -rhoq2(i_ns,k,i,j)+0.50_rp*ii_dni+is_dni-rs_dns )) ! => dns:minus
       ss_dns  = max( dt*pac(i_nsacns2ns,k,i,j)       , min(0.0_rp, -rhoq2(i_ns,k,i,j)+0.50_rp*ii_dni+is_dni-rs_dns-gs_dns ))
       !
       gg_dng  = max( dt*pac(i_ngacng2ng,k,i,j)       , min(0.0_rp, -rhoq2(i_ng,k,i,j) ))
       !
       ! total plus in mixed phase collection(clp_)
       ! mass
       ! if( temp(k,i,j) <= T00 )then sw=1, else sw=0
       clp_dqc =  0.0_rp
       clp_dqr = (-rg_dqg-rs_dqs-ri_dqi) * (1.0_rp-sw)
       clp_dqi = -ic_dqc
       clp_dqs = -sc_dqc-ii_dqi-is_dqi
       clp_dqg = -gc_dqc -gs_dqs  + (-rg_dqr-rs_dqr-rs_dqs-ri_dqr-ri_dqi) * sw
       ! num.( number only increase when a+b=>c,  dnc=-dna)
       clp_dnc = 0.0_rp
       clp_dnr = 0.0_rp
       clp_dni = 0.0_rp
       clp_dns = -ii_dni*0.5_rp
       clp_dng = (-rs_dnr-ri_dnr) * sw
       ! total minus in mixed phase collection(clm_)
       ! mass
       clm_dqc = gc_dqc+sc_dqc+ic_dqc
       clm_dqr = (rg_dqr+rs_dqr+ri_dqr) * sw
       clm_dqi = ri_dqi+ii_dqi+is_dqi
       clm_dqs = rs_dqs+gs_dqs
       clm_dqg = rg_dqg * (1.0_rp-sw)
       ! num.
       clm_dnc = gc_dnc+sc_dnc+ic_dnc
       clm_dnr = (rg_dnr+rs_dnr+ri_dnr) * sw
       clm_dni = ri_dni+ii_dni+is_dni
       clm_dns = rs_dns+ss_dns+gs_dns
       clm_dng = gg_dng + rg_dng * (1.0_rp-sw)

       ! partial conversion
       ! 08/05/08 [Mod] T.Mitsui
       pco_dqi = max( dt*pq(i_licon,k,i,j), -clp_dqi )
       pco_dqs = max( dt*pq(i_lscon,k,i,j), -clp_dqs )
       pco_dqg = -pco_dqi-pco_dqs
       ! 08/05/08 [Mod] T.Mitsui
       pco_dni = max( dt*pq(i_nicon,k,i,j), -clp_dni )
       pco_dns = max( dt*pq(i_nscon,k,i,j), -clp_dns )
       pco_dng = -pco_dni-pco_dns
       ! enhanced melting ( always negative value )
       ! ice-cloud melting produces cloud, others produce rain
       eml_dqi =  max( dt*pq(i_liacm,k,i,j), min(0.0_rp, -rhoq2(i_qi,k,i,j)-(clp_dqi+clm_dqi)-pco_dqi ))
       eml_dqs =  max( dt*pq(i_lsacm,k,i,j), min(0.0_rp, -rhoq2(i_qs,k,i,j)-(clp_dqs+clm_dqs)-pco_dqs ))
       eml_dqg =  max( dt*(pq(i_lgacm,k,i,j)+pq(i_lgarm,k,i,j)+pq(i_lsarm,k,i,j)+pq(i_liarm,k,i,j)), &
                  min(0.0_rp, -rhoq2(i_qg,k,i,j)-(clp_dqg+clm_dqg)-pco_dqg ))
       eml_dqc = -eml_dqi
       eml_dqr = -eml_dqs-eml_dqg
       !
       eml_dni =  max( dt*pq(i_niacm,k,i,j), min(0.0_rp, -rhoq2(i_ni,k,i,j)-(clp_dni+clm_dni)-pco_dni ))
       eml_dns =  max( dt*pq(i_nsacm,k,i,j), min(0.0_rp, -rhoq2(i_ns,k,i,j)-(clp_dns+clm_dns)-pco_dns ))
       eml_dng =  max( dt*(pq(i_ngacm,k,i,j)+pq(i_ngarm,k,i,j)+pq(i_nsarm,k,i,j)+pq(i_niarm,k,i,j)), &
                  min(0.0_rp, -rhoq2(i_ng,k,i,j)-(clp_dng+clm_dng)-pco_dng ))
       eml_dnc = -eml_dni
       eml_dnr = -eml_dns-eml_dng
       !
       ! ice multiplication
       spl_dqg = max( dt*pq(i_lgspl,k,i,j), min(0.0_rp, -rhoq2(i_qg,k,i,j)-(clp_dqg+clm_dqg)-pco_dqg-eml_dqg ))
       spl_dqs = max( dt*pq(i_lsspl,k,i,j), min(0.0_rp, -rhoq2(i_qs,k,i,j)-(clp_dqs+clm_dqs)-pco_dqs-eml_dqs ))
       spl_dqi = -spl_dqg-spl_dqs
       fac9    = (spl_dqg-eps)/(dt*pq(i_lgspl,k,i,j)-eps)
       fac10   = (spl_dqs-eps)/(dt*pq(i_lsspl,k,i,j)-eps)
       spl_dni = dt*pq(i_nispl,k,i,j)*fac9*fac10
       !
       ! total cloud change
       drhogqc = (wrm_dqc + clp_dqc + clm_dqc           + eml_dqc )
       drhognc = (wrm_dnc + clp_dnc + clm_dnc           + eml_dnc )
       ! total rain change
       drhogqr = (wrm_dqr + clp_dqr + clm_dqr           + eml_dqr )
       drhognr = (wrm_dnr + clp_dnr + clm_dnr           + eml_dnr )
       ! total ice change
       drhogqi = (          clp_dqi + clm_dqi + pco_dqi + eml_dqi + spl_dqi)
       drhogni = (          clp_dni + clm_dni + pco_dni + eml_dni + spl_dni)
       ! total snow change
       drhogqs = (          clp_dqs + clm_dqs + pco_dqs + eml_dqs + spl_dqs)
       drhogns = (          clp_dns + clm_dns + pco_dns + eml_dns )
       ! total graupel change
       drhogqg = (          clp_dqg + clm_dqg + pco_dqg + eml_dqg + spl_dqg)
       drhogng = (          clp_dng + clm_dng + pco_dng + eml_dng )
       !
       !--- update
       !
       rhoq(i_qc,k,i,j) = max(0.0_rp, rhoq(i_qc,k,i,j) + drhogqc )
       rhoq(i_nc,k,i,j) = max(0.0_rp, rhoq(i_nc,k,i,j) + drhognc )
       rhoq(i_qr,k,i,j) = max(0.0_rp, rhoq(i_qr,k,i,j) + drhogqr )
       rhoq(i_nr,k,i,j) = max(0.0_rp, rhoq(i_nr,k,i,j) + drhognr )
       rhoq(i_qi,k,i,j) = max(0.0_rp, rhoq(i_qi,k,i,j) + drhogqi )
       rhoq(i_ni,k,i,j) = max(0.0_rp, rhoq(i_ni,k,i,j) + drhogni )
       rhoq(i_qs,k,i,j) = max(0.0_rp, rhoq(i_qs,k,i,j) + drhogqs )
       rhoq(i_ns,k,i,j) = max(0.0_rp, rhoq(i_ns,k,i,j) + drhogns )
       rhoq(i_qg,k,i,j) = max(0.0_rp, rhoq(i_qg,k,i,j) + drhogqg )
       rhoq(i_ng,k,i,j) = max(0.0_rp, rhoq(i_ng,k,i,j) + drhogng )
       !
       ! update
       ! rhogq = l*gsgam
       rhoe(k,i,j) = rhoe(k,i,j) + lhf * ( drhogqi + drhogqs + drhogqg )
    enddo
    enddo
    enddo
    profile_stop("sn14_update_rhoq")

    call prof_rapend  ('MP_Collection', 3)

    call prof_rapstart('MP_Postprocess', 3)

    !--- update mixing ratio
    do j  = js, je
    do i  = is, ie
    do k  = ks, ke
       do iq = 1,  qa
          qtrc(k,i,j,iq) = rhoq(iq,k,i,j) * rrho(k,i,j)
       enddo

       calc_qdry( qdry(k,i,j), qtrc, k, i, j, iq )
       calc_cv( cva(k,i,j), qdry(k,i,j), qtrc, k, i, j, iq, cvdry, aq_cv )
       calc_r( rmoist, qtrc(k,i,j,i_qv), qdry(k,i,j), rdry, rvap )
       cpa(k,i,j) = cva(k,i,j) + rmoist
       temp(k,i,j) = rhoe(k,i,j) / ( dens(k,i,j) * cva(k,i,j) )
       pres(k,i,j) = dens(k,i,j) * rmoist * temp(k,i,j)
       rhot(k,i,j) = temp(k,i,j) * ( p00 / pres(k,i,j) )**(rmoist/cpa(k,i,j)) &
               * dens(k,i,j)
    enddo
    enddo
    enddo

!    if( opt_debug )     call debugreport_collection
    if( opt_debug_tem ) call debug_tem_kij( 4, temp(:,:,:), dens(:,:,:), pres(:,:,:), qtrc(:,:,:,i_qv) )

    do iq = 1, qa
       do j = js, je
       do i = is, ie
       do k = ks, ke
          qtrc(k,i,j,iq) = rhoq(iq,k,i,j) * rrho(k,i,j)
       enddo
       enddo
       enddo
    enddo

    call prof_rapend  ('MP_Postprocess', 3)

    !----------------------------------------------------------------------------
    !
    ! 4.Saturation adjustment
    !
    !----------------------------------------------------------------------------
    call prof_rapstart('MP_Saturation_adjustment', 3)
    ! nothing to do
    call prof_rapend  ('MP_Saturation_adjustment', 3)
    !----------------------------------------------------------------------------
    !
    ! 5. Sedimentation ( terminal velocity must be negative )
    !
    !----------------------------------------------------------------------------
    call prof_rapstart('MP_Sedimentation', 3)

    if ( mp_doprecipitation ) then

    do j = js, je
    do i = is, ie
    do k = ks-1, ke
       flx_rain(k,i,j) = 0.0_rp
       flx_snow(k,i,j) = 0.0_rp
    enddo
    enddo
    enddo

    velw(:,:,:,:) = 0.0_rp

    do step = 1, mp_nstep_sedimentation

       call mp_terminal_velocity( velw(:,:,:,:), & ! [OUT]
                                  rhoq(:,:,:,:), & ! [IN]
                                  dens(:,:,:),   & ! [IN]
                                  temp(:,:,:),   & ! [IN]
                                  pres(:,:,:)    ) ! [IN]

       call mp_precipitation( wflux_rain(:,:,:),     & ! [OUT]
                              wflux_snow(:,:,:),     & ! [OUT]
                              dens(:,:,:),     & ! [INOUT]
                              momz(:,:,:),     & ! [INOUT]
                              momx(:,:,:),     & ! [INOUT]
                              momy(:,:,:),     & ! [INOUT]
                              rhoe(:,:,:),     & ! [INOUT]
                              qtrc(:,:,:,:),   & ! [INOUT]
                              velw(:,:,:,:),   & ! [IN]
                              temp(:,:,:),     & ! [IN]
                              mp_dtsec_sedimentation ) ! [IN]

       do j = js, je
       do i = is, ie
       do k = ks-1, ke
          flx_rain(k,i,j) = flx_rain(k,i,j) + wflux_rain(k,i,j) * mp_rnstep_sedimentation
          flx_snow(k,i,j) = flx_snow(k,i,j) + wflux_snow(k,i,j) * mp_rnstep_sedimentation
       enddo
       enddo
       enddo

    enddo

    endif

    do j = js, je
    do i = is, ie
       sflx_rain(i,j) = flx_rain(ks-1,i,j)
       sflx_snow(i,j) = flx_snow(ks-1,i,j)
    end do
    end do

    call prof_rapend  ('MP_Sedimentation', 3)

    return
  end subroutine mp_sn14

  !-----------------------------------------------------------------------------
  subroutine debug_tem_kij( &
      point,    &
      tem,      &
      rho,      &
      pre,      &
      qv        )
    use scale_process, only: &
       prc_myrank
    implicit none

    integer, intent(in) :: point
    real(RP), intent(in) :: tem(ka,ia,ja)
    real(RP), intent(in) :: rho(ka,ia,ja)
    real(RP), intent(in) :: pre(ka,ia,ja)
    real(RP), intent(in) :: qv (ka,ia,ja)

    integer :: k ,i, j
    !---------------------------------------------------------------------------

    do j = js, je
    do i = is, ie
    do k = ks, ke
       if (      tem(k,i,j) < tem_min &
            .OR. rho(k,i,j) < rho_min &
            .OR. pre(k,i,j) < 1.0_rp    ) then

          if( io_l ) write(io_fid_log,'(A,I3,A,4(F16.5),3(I6))') &
          "*** point: ", point, " low tem,rho,pre:", tem(k,i,j), rho(k,i,j), pre(k,i,j), qv(k,i,j), k, i, j, prc_myrank
       endif
    enddo
    enddo
    enddo

    return
  end subroutine debug_tem_kij

  subroutine nucleation_kij( &
       z, velz,          &
       rho, tem, pre,    &
       rhoq,             &
       PQ,               &
       cpa,              & ! in
       dTdt_rad,         & ! in
       qke,              & ! in
       CCN,               & ! in
       dt                ) ! in
    use scale_process, only: &
       prc_mpistop
    use scale_atmos_saturation, only: &
       moist_psat_liq       => atmos_saturation_psat_liq, &
       moist_psat_ice       => atmos_saturation_psat_ice,   &
       moist_pres2qsat_liq  => atmos_saturation_pres2qsat_liq, &
       moist_pres2qsat_ice  => atmos_saturation_pres2qsat_ice,   &
       moist_dqsi_dtem_rho  => atmos_saturation_dqsi_dtem_rho
    implicit none

    real(RP), intent(in)  :: z(ka)      !
    real(RP), intent(in)  :: velz(ka,ia,ja)   ! w of half point
    real(RP), intent(in)  :: rho(ka,ia,ja)    ! [Add] 09/08/18 T.Mitsui
    real(RP), intent(in)  :: tem(ka,ia,ja)    ! [Add] 09/08/18 T.Mitsui
    real(RP), intent(in)  :: pre(ka,ia,ja)    ! [Add] 09/08/18 T.Mitsui
    !
    real(RP), intent(in)  :: rhoq(qa,ka,ia,ja)     !
    real(RP), intent(out) :: PQ(pq_max,ka,ia,ja)
    !
    real(RP), intent(in) ::  cpa(ka,ia,ja)      ! in  09/08/18 [Add] T.Mitsui
    real(RP), intent(in)  :: dTdt_rad(ka,ia,ja) ! 09/08/18 T.Mitsui
    real(RP), intent(in)  :: qke(ka,ia,ja)      ! 09/08/18 T.Mitsui
    real(DP), intent(in)  :: dt
    real(RP), intent(in)  :: CCN(ka,ia,ja)
    !
    ! namelist variables
    !
    ! total aerosol number concentration [/m3]
    real(RP), parameter :: c_ccn_ocean= 1.00e+8_rp
    real(RP), parameter :: c_ccn_land = 1.26e+9_rp
    real(RP), save      :: c_ccn      = 1.00e+8_rp
    ! aerosol activation factor
    real(RP), parameter :: kappa_ocean= 0.462_rp
    real(RP), parameter :: kappa_land = 0.308_rp
    real(RP), save      :: kappa      = 0.462_rp
    real(RP), save      :: c_in       = 1.0_rp
    ! SB06 (36)
    real(RP), save :: nm_M92 = 1.e+3_rp
    real(RP), save :: am_M92 = -0.639_rp
    real(RP), save :: bm_M92 = 12.96_rp
    !
    real(RP), save :: in_max = 1000.e+3_rp ! max num. of Ice-Nuclei [num/m3]
    real(RP), save :: ssi_max= 0.60_rp
    real(RP), save :: ssw_max= 1.1_rp  ! [%]
    !
    logical, save :: flag_first = .true.
    real(RP), save :: qke_min = 0.03_rp ! sigma=0.1[m/s], 09/08/18 T.Mitsui
    real(RP), save :: tem_ccn_low=233.150_rp  ! = -40 degC  ! [Add] 10/08/03 T.Mitsui
    real(RP), save :: tem_in_low =173.150_rp  ! = -100 degC ! [Add] 10/08/03 T.Mitsui
    logical, save :: nucl_twomey = .false.
    logical, save :: inucl_w     = .false.
    !
    namelist /nm_mp_sn14_nucleation/ &
         in_max,                     & !
         c_ccn, kappa,               & ! cloud nucleation
         nm_m92, am_m92, bm_m92,     & ! ice nucleation
         xc_ccn, xi_ccn,             &
         tem_ccn_low,                & ! [Add] 10/08/03 T.Mitsui
         tem_in_low,                 & ! [Add] 10/08/03 T.Mitsui
         ssw_max, ssi_max,           &
         nucl_twomey, inucl_w        ! [Add] 13/01/30 Y.Sato
    !
!    real(RP) :: c_ccn_map(1,IA,JA)   ! c_ccn horizontal distribution
!    real(RP) :: kappa_map(1,IA,JA)   ! kappa horizontal distribution
!    real(RP) :: c_in_map(1,IA,JA)    ! c_in  horizontal distribution ! [Add] 11/08/30 T.Mitsui
    real(RP) :: esw(ka,ia,ja)      ! saturation vapor pressure, water
    real(RP) :: esi(ka,ia,ja)      !                            ice
    real(RP) :: ssw(ka,ia,ja)      ! super saturation (water)
    real(RP) :: ssi(ka,ia,ja)      ! super saturation (ice)
!    real(RP) :: w_dsswdz(KA,IA,JA) ! w*(d_ssw/ d_z) super saturation(water) flux
    real(RP) :: w_dssidz(ka,ia,ja) ! w*(d_ssi/ d_z), 09/04/14 T.Mitsui
!    real(RP) :: ssw_below(KA,IA,JA)! ssw(k-1)
    real(RP) :: ssi_below(ka,ia,ja)! ssi(k-1), 09/04/14 T.Mitsui
    real(RP) :: z_below(ka,ia,ja)  ! z(k-1)
    real(RP) :: dz                   ! z(k)-z(k-1)
    real(RP) :: pv                   ! vapor pressure
    ! work variables for Twomey Equation.
    real(RP) :: qsw(ka,ia,ja)
    real(RP) :: qsi(ka,ia,ja)
    real(RP) :: dqsidtem_rho(ka,ia,ja)
    real(RP) :: dssidt_rad(ka,ia,ja)
!    real(RP) :: dni_ratio(KA,IA,JA)
    real(RP) :: wssi, wdssi
    !
!    real(RP) :: xi_nuc(1,IA,JA)    ! xi use the value @ cloud base
!    real(RP) :: alpha_nuc(1,IA,JA) ! alpha_nuc
!    real(RP) :: eta_nuc(1,IA,JA)   ! xi use the value @ cloud base
    !
    real(RP) :: sigma_w(ka,ia,ja)
    real(RP) :: weff(ka,ia,ja)
    real(RP) :: weff_max(ka,ia,ja)
    !
    real(RP) :: coef_ccn(ia,ja)
    real(RP) :: slope_ccn(ia,ja)
    real(RP) :: nc_new(ka,ia,ja)
    real(RP) :: nc_new_below(ka,ia,ja)
    real(RP) :: dnc_new
    real(RP) :: nc_new_max   ! Lohmann (2002),
    real(RP) :: a_max
    real(RP) :: b_max
    logical :: flag_nucleation(ka,ia,ja)
    !
    real(RP) :: r_gravity
    real(RP), parameter :: r_sqrt3=0.577350269_rp ! = sqrt(1.d0/3.d0)
    real(RP), parameter :: eps=1.e-30_rp
    !====> ! 09/08/18
    !
    real(RP) :: dlcdt_max, dli_max ! defined by supersaturation
    real(RP) :: dncdt_max, dni_max ! defined by supersaturation
    real(RP) :: rdt
    !
    integer :: i, j, k
    !
    if( flag_first )then
       rewind(io_fid_conf)
       read(io_fid_conf, nml=nm_mp_sn14_nucleation, end=100)
100    if( io_l ) write(io_fid_log, nml=nm_mp_sn14_nucleation)
       flag_first=.false.
       if( mp_couple_aerosol .AND. nucl_twomey ) then
        write(io_fid_log,*) "nucl_twomey should be false when MP_couple_aerosol is true, stop"
        call prc_mpistop
       endif
    endif
    !
!    c_ccn_map(1,:,:) = c_ccn
!    kappa_map(1,:,:) = kappa
!    c_in_map(1,:,:)  = c_in
    !
!    nc_uplim_d(1,:,:)  = c_ccn_map(1,:,:)*1.5_RP
    do j = js, je
    do i = is, ie
       nc_uplim_d(1,i,j)  = c_ccn*1.5_rp
    end do
    end do
    !
    rdt            = 1.0_rp/dt
    r_gravity      = 1.0_rp/grav
    !
    call moist_psat_liq     ( esw, tem )
    call moist_psat_ice     ( esi, tem )
    call moist_pres2qsat_liq( qsw, tem, pre )
    call moist_pres2qsat_ice( qsi, tem, pre )
    call moist_dqsi_dtem_rho( dqsidtem_rho, tem, rho )
    !
    ! Lohmann (2002),JAS, eq.(1) but changing unit [cm-3] => [m-3]
    a_max = 1.e+6_rp*0.1_rp*(1.e-6_rp)**1.27_rp
    b_max = 1.27_rp
    !
    ssi_max = 1.0_rp

    do j = js, je
    do i = is, ie
       do k = ks, ke
          pv        = rhoq(i_qv,k,i,j)*rvap*tem(k,i,j)
          ssw(k,i,j) = min( mp_ssw_lim, ( pv/esw(k,i,j)-1.0_rp ) )*100.0_rp
          ssi(k,i,j) = (pv/esi(k,i,j) - 1.00_rp)
!          ssw_below(k+1,i,j) = ssw(k,i,j)
          ssi_below(k+1,i,j) = ssi(k,i,j)
          z_below(k+1,i,j)   = z(k)
       end do
!       ssw_below(KS,i,j) = ssw(KS,i,j)
       ssi_below(ks,i,j) = ssi(ks,i,j)
       z_below(ks,i,j)   = z(ks-1)

       ! dS/dz is evaluated by first order upstream difference
       !***  Solution for Twomey Equation ***
!       coef_ccn(i,j)  = 1.E+6_RP*0.88_RP*(c_ccn_map(1,i,j)*1.E-6_RP)**(2.0_RP/(kappa_map(1,i,j) + 2.0_RP)) * &
       coef_ccn(i,j)  = 1.e+6_rp*0.88_rp*(c_ccn*1.e-6_rp)**(2.0_rp/(kappa + 2.0_rp)) &
!           * (70.0_RP)**(kappa_map(1,i,j)/(kappa_map(1,i,j) + 2.0_RP))
            * (70.0_rp)**(kappa/(kappa + 2.0_rp))
!       slope_ccn(i,j) = 1.5_RP*kappa_map(1,i,j)/(kappa_map(1,i,j) + 2.0_RP)
       slope_ccn(i,j) = 1.5_rp*kappa/(kappa + 2.0_rp)
       !
       do k=ks, ke
          sigma_w(k,i,j) = r_sqrt3*sqrt(max(qke(k,i,j),qke_min))
       end do
       sigma_w(ks-1,i,j) = sigma_w(ks,i,j)
       sigma_w(ke+1,i,j) = sigma_w(ke,i,j)
       ! effective vertical velocity
       do k=ks, ke-1
          weff(k,i,j) = 0.5_rp*(velz(k,i,j) + velz(k+1,i,j)) - cpa(k,i,j)*r_gravity*dtdt_rad(k,i,j)
       end do
       weff(ks-1,i,j) = weff(ks,i,j)
       weff(ke,i,j)   = weff(ke-1,i,j)

    end do
    end do
    !
    if( mp_couple_aerosol ) then

         do j = js, je
         do i = is, ie
         do k = ks, ke
            if( ssw(k,i,j) > 1.e-10_rp .AND. pre(k,i,j) > 300.e+2_rp ) then
               nc_new(k,i,j) = max( ccn(k,i,j), c_ccn )
            else
               nc_new(k,i,j) = 0.0_rp
            endif
         enddo
         enddo
         enddo

    else

      if( nucl_twomey ) then
        ! diagnose cloud condensation nuclei

        do j = js, je
        do i = is, ie
        do k = ks, ke
           ! effective vertical velocity (maximum vertical velocity in turbulent flow)
           weff_max(k,i,j) = weff(k,i,j) + sigma_w(k,i,j)
           ! large scale upward motion region and saturated
           if( (weff(k,i,j) > 1.e-8_rp) .AND. (ssw(k,i,j) > 1.e-10_rp)  .AND. pre(k,i,j) > 300.e+2_rp )then
              ! Lohmann (2002), eq.(1)
              nc_new_max   = coef_ccn(i,j)*weff_max(k,i,j)**slope_ccn(i,j)
              nc_new(k,i,j) = a_max*nc_new_max**b_max
           else
              nc_new(k,i,j) = 0.0_rp
           end if
        end do
        end do
        end do
      else
        ! calculate cloud condensation nuclei
          do j = js, je
          do i = is, ie
          do k = ks, ke
             if( ssw(k,i,j) > 1.e-10_rp .AND. pre(k,i,j) > 300.e+2_rp ) then
                nc_new(k,i,j) = c_ccn*ssw(k,i,j)**kappa
             else
                nc_new(k,i,j) = 0.0_rp
             endif
          enddo
          enddo
          enddo
      endif

    endif

    do j = js, je
    do i = is, ie
       do k = ks, ke
          ! nc_new is bound by upper limit
          if( nc_new(k,i,j) > nc_uplim_d(1,i,j) )then ! no more CCN
             flag_nucleation(k,i,j) = .false.
             nc_new_below(k+1,i,j)  = 1.e+30_rp
          else if( nc_new(k,i,j) > eps )then ! nucleation can occur
             flag_nucleation(k,i,j) = .true.
             nc_new_below(k+1,i,j)  = nc_new(k,i,j)
          else ! nucleation cannot occur(unsaturated or negative w)
             flag_nucleation(k,i,j) = .false.
             nc_new_below(k+1,i,j)  = 0.0_rp
          end if
       end do
       nc_new_below(ks,i,j) = 0.0_rp
!       do k=KS, KE
           ! search maximum value of nc_new
!          if(  ( nc_new(k,i,j) < nc_new_below(k,i,j) ) .OR. &
!               ( nc_new_below(k,i,j) > c_ccn_map(1,i,j)*0.05_RP ) )then ! 5% of c_ccn
!               ( nc_new_below(k,i,j) > c_ccn*0.05_RP ) )then ! 5% of c_ccn
!             flag_nucleation(k,i,j) = .false.
!          end if
!       end do

    end do
    end do

    if( mp_couple_aerosol ) then

         do j = js, je
         do i = is, ie
         do k = ks, ke
            ! nucleation occurs at only cloud base.
            ! if CCN is more than below parcel, nucleation newly occurs
            ! effective vertical velocity
            if ( flag_nucleation(k,i,j) .AND. & ! large scale upward motion region and saturated
                 tem(k,i,j) > tem_ccn_low ) then
               dlcdt_max    = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
               dncdt_max    = dlcdt_max/xc_min
!               dnc_new      = nc_new(k,i,j)-rhoq(I_NC,k,i,j)
               dnc_new      = nc_new(k,i,j)
               pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
               pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
            else
               pq(i_ncccn,k,i,j) = 0.0_rp
               pq(i_lcccn,k,i,j) = 0.0_rp
            end if
         end do
         end do
         end do

    else

      if( nucl_twomey ) then
         do j = js, je
         do i = is, ie
         do k = ks, ke
            ! nucleation occurs at only cloud base.
            ! if CCN is more than below parcel, nucleation newly occurs
            ! effective vertical velocity
            if ( flag_nucleation(k,i,j) .AND. & ! large scale upward motion region and saturated
                 tem(k,i,j) > tem_ccn_low .AND. &
                 nc_new(k,i,j) > rhoq(i_nc,k,i,j) ) then
               dlcdt_max    = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
               dncdt_max    = dlcdt_max/xc_min
               dnc_new      = nc_new(k,i,j)-rhoq(i_nc,k,i,j)
               pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
               pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
            else
               pq(i_ncccn,k,i,j) = 0.0_rp
               pq(i_lcccn,k,i,j) = 0.0_rp
            end if
         end do
         end do
         end do
      else
         do j = js, je
         do i = is, ie
         do k = ks, ke
            ! effective vertical velocity
            if(  tem(k,i,j) > tem_ccn_low .AND. &
                 nc_new(k,i,j) > rhoq(i_nc,k,i,j) ) then
               dlcdt_max    = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
               dncdt_max    = dlcdt_max/xc_min
               dnc_new      = nc_new(k,i,j)-rhoq(i_nc,k,i,j)
               pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
               pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
            else
               pq(i_ncccn,k,i,j) = 0.0_rp
               pq(i_lcccn,k,i,j) = 0.0_rp
            end if
         end do
         end do
         end do
      endif
    endif

    !
    ! ice nucleation
    !
    ! +++ NOTE ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    ! Based on Phillips etal.(2006).
    ! However this approach doesn't diagnose Ni itself but diagnose tendency.
    ! Original approach adjust Ni instantaneously .
    ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    do j = js, je
    do i = is, ie
    do k = ks, ke
       dz             = z(k) - z_below(k,i,j)
       w_dssidz(k,i,j) = velz(k,i,j)*(ssi(k,i,j) - ssi_below(k,i,j))/dz ! 09/04/14 [Add] T.Mitsui
       dssidt_rad(k,i,j) = -rhoq(i_qv,k,i,j)/(rho(k,i,j)*qsi(k,i,j)*qsi(k,i,j))*dqsidtem_rho(k,i,j)*dtdt_rad(k,i,j)
       dli_max        = (rhoq(i_qv,k,i,j) - esi(k,i,j)/(rvap*tem(k,i,j)))*rdt
       dni_max        = min( dli_max/xi_ccn, (in_max-rhoq(i_ni,k,i,j))*rdt )
       wdssi          = min( w_dssidz(k,i,j)+dssidt_rad(k,i,j), 0.01_rp)
       wssi           = min( ssi(k,i,j), ssi_max)
       ! SB06(34),(35)
       if(  ( wdssi    > eps        ) .AND. & !
            (tem(k,i,j) < 273.15_rp   ) .AND. & !
            (rhoq(i_ni,k,i,j)  < in_max     ) .AND. &
            (wssi      >= eps       ) )then   !
!             PNIccn(k,i,j) = min(dni_max, c_in_map(1,i,j)*bm_M92*nm_M92*0.3_RP*exp(0.3_RP*bm_M92*(wssi-0.1_RP))*wdssi)
          if( inucl_w ) then
             pq(i_niccn,k,i,j) = min(dni_max, c_in*bm_m92*nm_m92*0.3_rp*exp(0.3_rp*bm_m92*(wssi-0.1_rp))*wdssi)
          else
             pq(i_niccn,k,i,j) = min(dni_max, max(c_in*nm_m92*exp(0.3_rp*bm_m92*(wssi-0.1_rp) )-rhoq(i_ni,k,i,j),0.0_rp )*rdt )
          endif
          pq(i_liccn,k,i,j) = min(dli_max, pq(i_niccn,k,i,j)*xi_ccn )
          ! only for output
!             dni_ratio(k,i,j) = dssidt_rad(k,i,j)/( w_dssidz(k,i,j)+dssidt_rad(k,i,j) )
       else
          pq(i_niccn,k,i,j) = 0.0_rp
          pq(i_liccn,k,i,j) = 0.0_rp
       end if
    end do
    end do
    end do

    return
  end subroutine nucleation_kij
  !----------------------------
  subroutine ice_multiplication_kij( &
       PQ,           & ! out
       Pac,          & ! in
       tem, rhoq, xq ) ! in

    ! ice multiplication by splintering
    ! we consider Hallet-Mossop process
    use scale_specfunc, only: &
       gammafunc => sf_gamma
    implicit none
    real(RP), intent(out):: PQ(pq_max,ka,ia,ja)
    !
    real(RP), intent(in) :: Pac(pac_max,ka,ia,ja)
    real(RP), intent(in) :: tem(ka,ia,ja)
    real(RP), intent(in) :: rhoq(qa,ka,ia,ja)
    real(RP), intent(in) :: xq(5,ka,ia,ja)
    !
    logical, save :: flag_first      = .true.
    ! production of (350.d3 ice particles)/(cloud riming [g]) => 350*d6 [/kg]
    real(RP), parameter :: pice = 350.0e+6_rp
    ! production of (1 ice particle)/(250 cloud particles riming)
    real(RP), parameter :: pnc  = 250.0_rp
    real(RP), parameter :: rc_cr= 12.e-6_rp ! critical size[micron]
    ! temperature factor
    real(RP) :: fp
    ! work for incomplete gamma function
    real(RP), save :: xc_cr    ! mass[kg] of cloud with r=critical size[micron]
    real(RP), save :: alpha    ! slope parameter of gamma function
    real(RP), save :: gm, lgm  ! gamma(alpha), log(gamma(alpha))
    real(RP) :: igm            ! in complete gamma(x,alpha)
    real(RP) :: x
    ! coefficient of expansion using in calculation of igm
    real(RP) :: a0,a1,a2,a3,a4,a5
    real(RP) :: a6,a7,a8,a9,a10
    real(RP) :: an1,an2,b0,b1,b2,c0,c1,c2
    real(RP) :: d0,d1,d2,e1,e2,h0,h1,h2
    real(RP), parameter :: eps=1.0e-30_rp
    ! number of cloud droplets larger than 12 micron(radius).
    real(RP) :: n12
    !
    real(RP) :: wn, wni, wns, wng
    integer :: i, j, k
    !
    if( flag_first )then
       flag_first = .false.
       ! work for Incomplete Gamma function
       xc_cr = (2.0_rp*rc_cr/a_m(i_qc))**(1.0_rp/b_m(i_qc))
       alpha = (nu(i_qc)+1.0_rp)/mu(i_qc)
       gm    = gammafunc(alpha)
       lgm   = log(gm)
    end if
    !
    do j = js, je
    do i = is, ie
    do k = ks, ke
       ! Here we assume particle temperature is same as environment temperature.
       ! If you want to treat in a better manner,
       ! you can diagnose with eq.(64) in CT(86)
       if     (tem(k,i,j) >  270.16_rp)then
          fp   = 0.0_rp
       else if(tem(k,i,j) >= 268.16_rp)then
          fp   = (270.16_rp-tem(k,i,j))*0.5_rp
       else if(tem(k,i,j) >= 265.16_rp)then
          fp   = (tem(k,i,j)-265.16_rp)*0.333333333_rp
       else
          fp   = 0.0_rp
       end if
       ! Approximation of Incomplete Gamma function
       ! Here we calculate with algorithm by Numerical Recipes.
       ! This approach is based on eq.(78) in Cotton etal.(1986),
       ! but more accurate and expanded for Generalized Gamma Distribution.
       x       = coef_lambda(i_qc)*(xc_cr/xq(i_c,k,i,j))**mu(i_qc)
       !
       if(x<1.e-2_rp*alpha)then        ! negligible
          igm  = 0.0_rp
       else if(x<alpha+1.0_rp)then    ! series expansion
          ! 10th-truncation is enough for cloud droplet.
          a0   = 1.0_rp/alpha         ! n=0
          a1   = a0*x/(alpha+1.0_rp)  ! n=1
          a2   = a1*x/(alpha+2.0_rp)  ! n=2
          a3   = a2*x/(alpha+3.0_rp)  ! n=3
          a4   = a3*x/(alpha+4.0_rp)  ! n=4
          a5   = a4*x/(alpha+5.0_rp)  ! n=5
          a6   = a5*x/(alpha+6.0_rp)  ! n=6
          a7   = a6*x/(alpha+7.0_rp)  ! n=7
          a8   = a7*x/(alpha+8.0_rp)  ! n=8
          a9   = a8*x/(alpha+9.0_rp)  ! n=9
          a10  = a9*x/(alpha+10.0_rp) ! n=10
          igm  = (a0+a1+a2+a3+a4+a5+a6+a7+a8+a9+a10)*exp( -x + alpha*log(x) - lgm )
       else if(x<alpha*1.d2) then ! continued fraction expansion
          ! 2nd-truncation is enough for cloud droplet.
          ! setup
          b0   = x+1.0_rp-alpha
          c0   = 1.0_rp/eps
          d0   = 1.0_rp/b0
          h0   = d0
          ! n=1
          an1  = -(1.0_rp-alpha)
          b1   = b0 + 2.0_rp
          d1   = 1.0_rp/(an1*d0+b1)
          c1   = b1+an1/c0
          e1   = d1*c1
          h1   = h0*e1
          ! n=2
          an2  = -2.0_rp*(2.0_rp-alpha)
          b2   = b1 + 2.0_rp
          d2   = 1.0_rp/(an2*d1+b2)
          c2   = b2+an2/c1
          e2   = d2*c2
          h2   = h1*e2
          !
          igm  = 1.0_rp - exp( -x + alpha*log(x) - lgm )*h2
       else                       ! negligible
          igm  = 1.0_rp
       end if
       ! n12 is number of cloud droplets larger than 12 micron.
       n12 = rhoq(i_nc,k,i,j)*(1.0_rp-igm)
       ! eq.(82) CT(86)
       wn           = (pice + n12/((rhoq(i_qc,k,i,j)+xc_min)*pnc) )*fp ! filtered by xc_min
       wni          = wn*(-pac(i_liaclc2li,k,i,j)  ) ! riming production rate is all negative
       wns          = wn*(-pac(i_lsaclc2ls,k,i,j)  )
       wng          = wn*(-pac(i_lgaclc2lg,k,i,j)  )
       pq(i_nispl,k,i,j) = wni+wns+wng
       !
       pq(i_lsspl,k,i,j) = - wns*xq(i_i,k,i,j) ! snow    => ice
       pq(i_lgspl,k,i,j) = - wng*xq(i_i,k,i,j) ! graupel => ice
       !
    end do
    end do
    end do
    !
    return
  end subroutine ice_multiplication_kij
  !----------------------------
  subroutine mixed_phase_collection_kij( &
       ! collection process
       Pac, PQ,                          & ! out
       wtem, rhoq,                       & ! in
       xq, dq_xave,  vt_xave             ) ! in
!       rho                               ) ! [Add] 11/08/30
    use scale_atmos_saturation, only: &
       moist_psat_ice => atmos_saturation_psat_ice
    implicit none

    !--- mixed-phase collection process
    !                  And all we set all production term as a negative sign to avoid confusion.
    !
    real(RP), intent(out):: Pac(pac_max,ka,ia,ja)
    !--- partial conversion
    real(RP), intent(out):: PQ(pq_max,ka,ia,ja)
    !
    real(RP), intent(in) :: wtem(ka,ia,ja)
    !--- mass/number concentration[kg/m3]
    real(RP), intent(in) :: rhoq(qa,ka,ia,ja)
    ! necessary ?
    real(RP), intent(in) :: xq(5,ka,ia,ja)
    !--- diameter of averaged mass( D(ave x) )
    real(RP), intent(in) :: dq_xave(5,ka,ia,ja)
    !--- terminal velocity of averaged mass( vt(ave x) )
    real(RP), intent(in) :: vt_xave(5,2,ka,ia,ja)
    ! [Add] 11/08/30 T.Mitsui, for autoconversion of ice
!    real(RP), intent(in) :: rho(KA,IA,JA)
    !
    ! namelist variables
    !=== for collection
    !--- threshold of diameters to collide with others
    real(RP), save :: dc0 =  15.0e-6_rp       ! lower threshold of cloud
    real(RP), save :: dc1 =  40.0e-6_rp       ! upper threshold of cloud
    real(RP), save :: di0 = 150.0e-6_rp       ! lower threshold of cloud
    real(RP), save :: ds0 = 150.0e-6_rp       ! lower threshold of cloud
    real(RP), save :: dg0 = 150.0e-6_rp       ! lower threshold of cloud
    !--- standard deviation of terminal velocity[m/s]
    real(RP), save :: sigma_c=0.00_rp        ! cloud
    real(RP), save :: sigma_r=0.00_rp        ! rain
    real(RP), save :: sigma_i=0.2_rp        ! ice
    real(RP), save :: sigma_s=0.2_rp        ! snow
    real(RP), save :: sigma_g=0.00_rp        ! graupel
    !--- max collection efficiency for cloud
    real(RP), save :: E_im = 0.80_rp        ! ice max
    real(RP), save :: E_sm = 0.80_rp        ! snow max
    real(RP), save :: E_gm = 1.00_rp        ! graupel max
    !--- collection efficiency between 2 species
    real(RP), save :: E_ir=1.0_rp            ! ice     x rain
    real(RP), save :: E_sr=1.0_rp            ! snow    x rain
    real(RP), save :: E_gr=1.0_rp            ! graupel x rain
    real(RP), save :: E_ii=1.0_rp            ! ice     x ice
    real(RP), save :: E_si=1.0_rp            ! snow    x ice
    real(RP), save :: E_gi=1.0_rp            ! graupel x ice
    real(RP), save :: E_ss=1.0_rp            ! snow    x snow
    real(RP), save :: E_gs=1.0_rp            ! graupel x snow
    real(RP), save :: E_gg=1.0_rp            ! graupel x graupel
    !=== for partial conversion
    !--- flag: 1=> partial conversion to graupel, 0=> no conversion
    integer, save :: i_iconv2g=1          ! ice  => graupel
    integer, save :: i_sconv2g=1          ! snow => graupel
    !--- bulk density of graupel
    real(RP), save :: rho_g   = 900.0_rp     ! [kg/m3]
    !--- space filling coefficient [%]
    real(RP), save :: cfill_i = 0.68_rp     ! ice
    real(RP), save :: cfill_s = 0.01_rp     ! snow
    !--- critical diameter for ice conversion
    real(RP), save :: di_cri  = 500.e-6_rp    ! [m]
    ! [Add] 10/08/03 T.Mitsui
    logical, save :: opt_stick_KS96=.false.
    logical, save :: opt_stick_CO86=.false.
    real(RP), parameter :: a_dec = 0.883_rp
    real(RP), parameter :: b_dec = 0.093_rp
    real(RP), parameter :: c_dec = 0.00348_rp
    real(RP), parameter :: d_dec = 4.5185e-5_rp
    !
    logical, save :: flag_first = .true.
    namelist /nm_mp_sn14_collection/ &
         dc0, dc1, di0, ds0, dg0,    &
         sigma_c, sigma_r, sigma_i, sigma_s, sigma_g, &
         opt_stick_ks96,   &
         opt_stick_co86,   &
         e_im, e_sm, e_gm, &
         e_ir, e_sr, e_gr, e_ii, e_si, e_gi, e_ss, e_gs, e_gg, &
         i_iconv2g, i_sconv2g, rho_g, cfill_i, cfill_s, di_cri
    !
    real(RP) :: tem(ka,ia,ja)
    !
    !--- collection efficency of each specie
    real(RP) :: E_c, E_r, E_i, E_s, E_g    !
    real(RP) :: E_ic, E_sc, E_gc           !
    !--- sticking efficiency
    real(RP) :: E_stick(ka,ia,ja)
    ! [Add] 10/08/03 T.Mitsui
    real(RP) :: temc, temc2, temc3
    real(RP) :: E_dec
    real(RP) :: esi_rat
    real(RP) :: esi(ka,ia,ja)
    !
    real(RP) :: temc_p, temc_m             ! celcius tem.
    ! [Add] 11/08/30 T.Mitsui, estimation of autoconversion time
!    real(RP) :: ci_aut(KA,IA,JA)
!    real(RP) :: taui_aut(KA,IA,JA)
!    real(RP) :: tau_sce(KA,IA,JA)
    !--- DSD averaged diameter for each species
    real(RP) :: ave_dc                     ! cloud
!    real(RP) :: ave_dr                     ! rain
    real(RP) :: ave_di                     ! ice
    real(RP) :: ave_ds                     ! snow
    real(RP) :: ave_dg                     ! graupel
    !--- coefficient of collection equations(L:mass, N:number)
    real(RP) :: coef_acc_LCI,   coef_acc_NCI   ! cloud     - cloud ice
    real(RP) :: coef_acc_LCS,   coef_acc_NCS   ! cloud     - snow
    !
    real(RP) :: coef_acc_LCG,   coef_acc_NCG   ! cloud     - graupel
    real(RP) :: coef_acc_LRI_I, coef_acc_NRI_I ! rain      - cloud ice
    real(RP) :: coef_acc_LRI_R, coef_acc_NRI_R ! rain      - cloud ice
    real(RP) :: coef_acc_LRS_S, coef_acc_NRS_S ! rain      - snow
    real(RP) :: coef_acc_LRS_R, coef_acc_NRS_R ! rain      - snow
    real(RP) :: coef_acc_LRG,   coef_acc_NRG   ! rain      - graupel
    real(RP) :: coef_acc_LII,   coef_acc_NII   ! cloud ice - cloud ice
    real(RP) :: coef_acc_LIS,   coef_acc_NIS   ! cloud ice - snow
    real(RP) ::                 coef_acc_NSS   ! snow      - snow
    real(RP) ::                 coef_acc_NGG   ! grauepl   - graupel
    real(RP) :: coef_acc_LSG,   coef_acc_NSG   ! snow      - graupel
    !--- (diameter) x (diameter)
    real(RP) :: dcdc, dcdi, dcds, dcdg
    real(RP) :: drdr, drdi, drds, drdg
    real(RP) :: didi, dids, didg
    real(RP) :: dsds, dsdg
    real(RP) :: dgdg
    !--- (terminal velocity) x (terminal velocity)
    real(RP) :: vcvc, vcvi, vcvs, vcvg
    real(RP) :: vrvr, vrvi, vrvs, vrvg
    real(RP) :: vivi, vivs, vivg
    real(RP) :: vsvs, vsvg
    real(RP) :: vgvg
    !
    real(RP) :: wx_cri, wx_crs
    real(RP) :: coef_emelt
    real(RP) :: w1

    real(RP) :: sw
    !
    integer :: i, j, k
    !
    if( flag_first )then
       rewind( io_fid_conf )
       read( io_fid_conf, nml=nm_mp_sn14_collection, end=100 )
100    if( io_l ) write( io_fid_log, nml=nm_mp_sn14_collection )
       flag_first = .false.
    end if
    !
    ! [Add] 10/08/03 T.Mitsui
    do j = js, je
    do i = is, ie
    do k = ks, ke
       tem(k,i,j) = max( wtem(k,i,j), tem_min ) ! 11/08/30 T.Mitsui
    end do
    end do
    end do

    call moist_psat_ice( esi, tem )

    if( opt_stick_ks96 )then
       do j = js, je
       do i = is, ie
       do k = ks, ke
          ! Khain and Sednev (1996), eq.(3.15)
          temc          = tem(k,i,j) - t00
          temc2         = temc*temc
          temc3         = temc2*temc
          e_dec         = max(0.0_rp, a_dec + b_dec*temc + c_dec*temc2 + d_dec*temc3 )
          esi_rat       = rhoq(i_qv,k,i,j)*rvap*tem(k,i,j)/esi(k,i,j)
          e_stick(k,i,j) = min(1.0_rp, e_dec*esi_rat)
       end do
       end do
       end do
    else if( opt_stick_co86 )then
       do j = js, je
       do i = is, ie
       do k = ks, ke
          ! [Add] 11/08/30 T.Mitsui, Cotton et al. (1986)
          temc          = min(tem(k,i,j) - t00,0.0_rp)
          w1            = 0.035_rp*temc-0.7_rp
          e_stick(k,i,j) = 10._rp**w1
       end do
       end do
       end do
    else
       do j = js, je
       do i = is, ie
       do k = ks, ke
          ! Lin et al. (1983)
          temc_m        = min(tem(k,i,j) - t00,0.0_rp) ! T < 273.15
          e_stick(k,i,j) = exp(0.09_rp*temc_m)
       end do
       end do
       end do
    end if
    !
    profile_start("sn14_collection")
!OCL NOSIMD
    do j = js, je
    do i = is, ie
    do k = ks, ke
       !
!          temc_m = min(tem(k,i,j) - T00,0.0_RP) ! T < 273.15
       temc_p = max(tem(k,i,j) - t00,0.0_rp) ! T > 273.15
       ! averaged diameter using SB06(82)
       ave_dc = coef_d(i_qc)*xq(i_c,k,i,j)**b_m(i_qc)
       ave_di = coef_d(i_qi)*xq(i_i,k,i,j)**b_m(i_qi)
       ave_ds = coef_d(i_qs)*xq(i_s,k,i,j)**b_m(i_qs)
       ave_dg = coef_d(i_qg)*xq(i_g,k,i,j)**b_m(i_qg)
       !------------------------------------------------------------------------
       ! coellection efficiency are given as follows
       e_c = max(0.0_rp, min(1.0_rp, (ave_dc-dc0)/(dc1-dc0) ))
       sw = 0.5_rp - sign(0.5_rp, di0-ave_di) ! if(ave_di>di0)then sw=1
       e_i = e_im * sw
       sw = 0.5_rp - sign(0.5_rp, ds0-ave_ds) ! if(ave_ds>ds0)then sw=1
       e_s = e_sm * sw
       sw = 0.5_rp - sign(0.5_rp, dg0-ave_dg) ! if(ave_dg>dg0)then sw=1
       e_g = e_gm * sw
       e_ic = e_i*e_c
       e_sc = e_s*e_c
       e_gc = e_g*e_c
       !------------------------------------------------------------------------
       ! Collection:  a collects b ( assuming particle size a>b )
       dcdc = dq_xave(i_c,k,i,j) * dq_xave(i_c,k,i,j)
       drdr = dq_xave(i_r,k,i,j) * dq_xave(i_r,k,i,j)
       didi = dq_xave(i_i,k,i,j) * dq_xave(i_i,k,i,j)
       dsds = dq_xave(i_s,k,i,j) * dq_xave(i_s,k,i,j)
       dgdg = dq_xave(i_g,k,i,j) * dq_xave(i_g,k,i,j)
       dcdi = dq_xave(i_c,k,i,j) * dq_xave(i_i,k,i,j)
       dcds = dq_xave(i_c,k,i,j) * dq_xave(i_s,k,i,j)
       dcdg = dq_xave(i_c,k,i,j) * dq_xave(i_g,k,i,j)
       drdi = dq_xave(i_r,k,i,j) * dq_xave(i_i,k,i,j)
       drds = dq_xave(i_r,k,i,j) * dq_xave(i_s,k,i,j)
       drdg = dq_xave(i_r,k,i,j) * dq_xave(i_g,k,i,j)
       dids = dq_xave(i_i,k,i,j) * dq_xave(i_s,k,i,j)
!       didg = dq_xave(I_I,k,i,j) * dq_xave(I_G,k,i,j)
       dsdg = dq_xave(i_s,k,i,j) * dq_xave(i_g,k,i,j)
       !
       vcvc = vt_xave(i_c,2,k,i,j)* vt_xave(i_c,2,k,i,j)
       vrvr = vt_xave(i_r,2,k,i,j)* vt_xave(i_r,2,k,i,j)
       vivi = vt_xave(i_i,2,k,i,j)* vt_xave(i_i,2,k,i,j)
       vsvs = vt_xave(i_s,2,k,i,j)* vt_xave(i_s,2,k,i,j)
       vgvg = vt_xave(i_g,2,k,i,j)* vt_xave(i_g,2,k,i,j)
       vcvi = vt_xave(i_c,2,k,i,j)* vt_xave(i_i,2,k,i,j)
       vcvs = vt_xave(i_c,2,k,i,j)* vt_xave(i_s,2,k,i,j)
       vcvg = vt_xave(i_c,2,k,i,j)* vt_xave(i_g,2,k,i,j)
       vrvi = vt_xave(i_r,2,k,i,j)* vt_xave(i_i,2,k,i,j)
       vrvs = vt_xave(i_r,2,k,i,j)* vt_xave(i_s,2,k,i,j)
       vrvg = vt_xave(i_r,2,k,i,j)* vt_xave(i_g,2,k,i,j)
       vivs = vt_xave(i_i,2,k,i,j)* vt_xave(i_s,2,k,i,j)
!       vivg = vt_xave(I_I,2,k,i,j)* vt_xave(I_G,2,k,i,j)
       vsvg = vt_xave(i_s,2,k,i,j)* vt_xave(i_g,2,k,i,j)
       !------------------------------------------------------------------------
       !
       !+++ pattern 1: a + b => a  (a>b)
       !                           (i-c, s-c, g-c, s-i, g-r, s-g)
       !------------------------------------------------------------------------
       ! cloud-ice => ice
       ! reduction term of cloud
       coef_acc_lci = &
              ( delta_b1(i_qc)*dcdc + delta_ab1(i_qi,i_qc)*dcdi + delta_b0(i_qi)*didi ) &
            * ( theta_b1(i_qc)*vcvc - theta_ab1(i_qi,i_qc)*vcvi + theta_b0(i_qi)*vivi   &
            +  sigma_i + sigma_c )**0.5_rp
       coef_acc_nci = &
              ( delta_b0(i_qc)*dcdc + delta_ab0(i_qi,i_qc)*dcdi + delta_b0(i_qi)*didi ) &
            * ( theta_b0(i_qc)*vcvc - theta_ab0(i_qi,i_qc)*vcvi + theta_b0(i_qi)*vivi   &
            +  sigma_i + sigma_c )**0.5_rp
       pac(i_liaclc2li,k,i,j)= -0.25_rp*pi*e_ic*rhoq(i_ni,k,i,j)*rhoq(i_qc,k,i,j)*coef_acc_lci
       pac(i_niacnc2ni,k,i,j)= -0.25_rp*pi*e_ic*rhoq(i_ni,k,i,j)*rhoq(i_nc,k,i,j)*coef_acc_nci
       ! cloud-snow => snow
       ! reduction term of cloud
       coef_acc_lcs = &
              ( delta_b1(i_qc)*dcdc + delta_ab1(i_qs,i_qc)*dcds + delta_b0(i_qs)*dsds ) &
            * ( theta_b1(i_qc)*vcvc - theta_ab1(i_qs,i_qc)*vcvs + theta_b0(i_qs)*vsvs   &
            +  sigma_s + sigma_c )**0.5_rp
       coef_acc_ncs = &
              ( delta_b0(i_qc)*dcdc + delta_ab0(i_qs,i_qc)*dcds + delta_b0(i_qs)*dsds ) &
            * ( theta_b0(i_qc)*vcvc - theta_ab0(i_qs,i_qc)*vcvs + theta_b0(i_qs)*vsvs   &
            +  sigma_s + sigma_c )**0.5_rp
       pac(i_lsaclc2ls,k,i,j)= -0.25_rp*pi*e_sc*rhoq(i_ns,k,i,j)*rhoq(i_qc,k,i,j)*coef_acc_lcs
       pac(i_nsacnc2ns,k,i,j)= -0.25_rp*pi*e_sc*rhoq(i_ns,k,i,j)*rhoq(i_nc,k,i,j)*coef_acc_ncs
       ! cloud-graupel => graupel
       ! reduction term of cloud
       coef_acc_lcg = &
              ( delta_b1(i_qc)*dcdc + delta_ab1(i_qg,i_qc)*dcdg + delta_b0(i_qg)*dgdg ) &
            * ( theta_b1(i_qc)*vcvc - theta_ab1(i_qg,i_qc)*vcvg + theta_b0(i_qg)*vgvg   &
            +  sigma_g + sigma_c )**0.5_rp
       coef_acc_ncg = &
              ( delta_b0(i_qc)*dcdc + delta_ab0(i_qg,i_qc)*dcdg + delta_b0(i_qg)*dgdg ) &
            * ( theta_b0(i_qc)*vcvc - theta_ab0(i_qg,i_qc)*vcvg + theta_b0(i_qg)*vgvg   &
            +  sigma_g + sigma_c )**0.5_rp
       pac(i_lgaclc2lg,k,i,j)= -0.25_rp*pi*e_gc*rhoq(i_ng,k,i,j)*rhoq(i_qc,k,i,j)*coef_acc_lcg
       pac(i_ngacnc2ng,k,i,j)= -0.25_rp*pi*e_gc*rhoq(i_ng,k,i,j)*rhoq(i_nc,k,i,j)*coef_acc_ncg
       ! snow-graupel => graupel
       coef_acc_lsg = &
              ( delta_b1(i_qs)*dsds + delta_ab1(i_qg,i_qs)*dsdg + delta_b0(i_qg)*dgdg ) &
            * ( theta_b1(i_qs)*vsvs - theta_ab1(i_qg,i_qs)*vsvg + theta_b0(i_qg)*vgvg   &
            +  sigma_g + sigma_s )**0.5_rp
       coef_acc_nsg = &
              ( delta_b0(i_qs)*dsds + delta_ab0(i_qg,i_qs)*dsdg + delta_b0(i_qg)*dgdg ) &
            ! [fix] T.Mitsui 08/05/08
            * ( theta_b0(i_qs)*vsvs - theta_ab0(i_qg,i_qs)*vsvg + theta_b0(i_qg)*vgvg   &
            +  sigma_g + sigma_s )**0.5_rp
       pac(i_lgacls2lg,k,i,j)= -0.25_rp*pi*e_stick(k,i,j)*e_gs*rhoq(i_ng,k,i,j)*rhoq(i_qs,k,i,j)*coef_acc_lsg
       pac(i_ngacns2ng,k,i,j)= -0.25_rp*pi*e_stick(k,i,j)*e_gs*rhoq(i_ng,k,i,j)*rhoq(i_ns,k,i,j)*coef_acc_nsg
       !------------------------------------------------------------------------
       ! ice-snow => snow
       ! reduction term of ice
       coef_acc_lis = &
              ( delta_b1(i_qi)*didi + delta_ab1(i_qs,i_qi)*dids + delta_b0(i_qs)*dsds ) &
            * ( theta_b1(i_qi)*vivi - theta_ab1(i_qs,i_qi)*vivs + theta_b0(i_qs)*vsvs   &
            +  sigma_i + sigma_s )**0.5_rp
       coef_acc_nis = &
              ( delta_b0(i_qi)*didi + delta_ab0(i_qs,i_qi)*dids + delta_b0(i_qs)*dsds ) &
            * ( theta_b0(i_qi)*vivi - theta_ab0(i_qs,i_qi)*vivs + theta_b0(i_qs)*vsvs   &
            +  sigma_i + sigma_s )**0.5_rp
       pac(i_liacls2ls,k,i,j)= -0.25_rp*pi*e_stick(k,i,j)*e_si*rhoq(i_ns,k,i,j)*rhoq(i_qi,k,i,j)*coef_acc_lis
       pac(i_niacns2ns,k,i,j)= -0.25_rp*pi*e_stick(k,i,j)*e_si*rhoq(i_ns,k,i,j)*rhoq(i_ni,k,i,j)*coef_acc_nis
       !
       sw = sign(0.5_rp, t00-tem(k,i,j)) + 0.5_rp
       ! if ( tem(k,i,j) <= T00 )then
          ! rain-graupel => graupel
          ! reduction term of rain
          ! sw = 1
       ! else
          ! rain-graupel => rain
          ! reduction term of graupel
          ! sw = 0
       coef_acc_lrg = &
              ( ( delta_b1(i_qr)*drdr + delta_ab1(i_qg,i_qr)*drdg + delta_b0(i_qg)*dgdg ) * sw &
              + ( delta_b1(i_qg)*dgdg + delta_ab1(i_qr,i_qg)*drdg + delta_b0(i_qr)*drdr ) * (1.0_rp-sw) ) &
            * sqrt( ( theta_b1(i_qr)*vrvr - theta_ab1(i_qg,i_qr)*vrvg + theta_b0(i_qg)*vgvg ) * sw &
                  + ( theta_b1(i_qg)*vgvg - theta_ab1(i_qr,i_qg)*vrvg + theta_b0(i_qr)*vrvr ) * (1.0_rp-sw) &
                  + sigma_r + sigma_g )
       pac(i_lraclg2lg,k,i,j) = -0.25_rp*pi*e_gr*coef_acc_lrg &
            * ( rhoq(i_ng,k,i,j)*rhoq(i_qr,k,i,j) * sw &
              + rhoq(i_nr,k,i,j)*rhoq(i_qg,k,i,j) * (1.0_rp-sw) )
       coef_acc_nrg = &
              ( delta_b0(i_qr)*drdr + delta_ab0(i_qg,i_qr)*drdg + delta_b0(i_qg)*dgdg ) &
            * ( theta_b0(i_qr)*vrvr - theta_ab0(i_qg,i_qr)*vrvg + theta_b0(i_qg)*vgvg   &
            +  sigma_r + sigma_g )**0.5_rp
       pac(i_nracng2ng,k,i,j) = -0.25_rp*pi*e_gr*rhoq(i_ng,k,i,j)*rhoq(i_nr,k,i,j)*coef_acc_nrg
       !
       !------------------------------------------------------------------------
       !
       !+++ pattern 2: a + b => c  (a>b)
       !                           (r-i,r-s)
       !------------------------------------------------------------------------
       ! rain-ice => graupel
       ! reduction term of ice
       coef_acc_lri_i = &
              ( delta_b1(i_qi)*didi + delta_ab1(i_qr,i_qi)*drdi + delta_b0(i_qr)*drdr ) &
            * ( theta_b1(i_qi)*vivi - theta_ab1(i_qr,i_qi)*vrvi + theta_b0(i_qr)*vrvr   &
            +  sigma_r + sigma_i )**0.5_rp
       coef_acc_nri_i = &
              ( delta_b0(i_qi)*didi + delta_ab0(i_qr,i_qi)*drdi + delta_b0(i_qr)*drdr ) &
            * ( theta_b0(i_qi)*vivi - theta_ab0(i_qr,i_qi)*vrvi + theta_b0(i_qr)*vrvr   &
            +  sigma_r + sigma_i )**0.5_rp
       pac(i_lracli2lg_i,k,i,j)= -0.25_rp*pi*e_ir*rhoq(i_nr,k,i,j)*rhoq(i_qi,k,i,j)*coef_acc_lri_i
       pac(i_nracni2ng_i,k,i,j)= -0.25_rp*pi*e_ir*rhoq(i_nr,k,i,j)*rhoq(i_ni,k,i,j)*coef_acc_nri_i
       ! reduction term of rain
       coef_acc_lri_r = &
              ( delta_b1(i_qr)*drdr + delta_ab1(i_qi,i_qr)*drdi + delta_b0(i_qi)*didi ) &
            * ( theta_b1(i_qr)*vrvr - theta_ab1(i_qi,i_qr)*vrvi + theta_b0(i_qi)*vivi   &
            +  sigma_r + sigma_i )**0.5_rp
       coef_acc_nri_r = &
              ( delta_b0(i_qr)*drdr + delta_ab0(i_qi,i_qr)*drdi + delta_b0(i_qi)*didi ) &
            * ( theta_b0(i_qr)*vrvr - theta_ab0(i_qi,i_qr)*vrvi + theta_b0(i_qi)*vivi   &
            +  sigma_r + sigma_i )**0.5_rp
       pac(i_lracli2lg_r,k,i,j)= -0.25_rp*pi*e_ir*rhoq(i_ni,k,i,j)*rhoq(i_qr,k,i,j)*coef_acc_lri_r
       pac(i_nracni2ng_r,k,i,j)= -0.25_rp*pi*e_ir*rhoq(i_ni,k,i,j)*rhoq(i_nr,k,i,j)*coef_acc_nri_r
       ! rain-snow => graupel
       ! reduction term of snow
       coef_acc_lrs_s = &
              ( delta_b1(i_qs)*dsds + delta_ab1(i_qr,i_qs)*drds + delta_b0(i_qr)*drdr ) &
            * ( theta_b1(i_qs)*vsvs - theta_ab1(i_qr,i_qs)*vrvs + theta_b0(i_qr)*vrvr   &
            +  sigma_r + sigma_s )**0.5_rp
       coef_acc_nrs_s = &
              ( delta_b0(i_qs)*dsds + delta_ab0(i_qr,i_qs)*drds + delta_b0(i_qr)*drdr ) &
            * ( theta_b0(i_qs)*vsvs - theta_ab0(i_qr,i_qs)*vrvs + theta_b0(i_qr)*vrvr   &
            +  sigma_r + sigma_s )**0.5_rp
       pac(i_lracls2lg_s,k,i,j)= -0.25_rp*pi*e_sr*rhoq(i_nr,k,i,j)*rhoq(i_qs,k,i,j)*coef_acc_lrs_s
       pac(i_nracns2ng_s,k,i,j)= -0.25_rp*pi*e_sr*rhoq(i_nr,k,i,j)*rhoq(i_ns,k,i,j)*coef_acc_nrs_s
       ! reduction term of rain
       coef_acc_lrs_r = &
              ( delta_b1(i_qr)*drdr + delta_ab1(i_qs,i_qr)*drds + delta_b0(i_qs)*dsds ) &
            * ( theta_b1(i_qr)*vrvr - theta_ab1(i_qs,i_qr)*vrvs + theta_b0(i_qs)*vsvs   &
            +  sigma_r + sigma_s )**0.5_rp
       coef_acc_nrs_r = &
              ( delta_b0(i_qr)*drdr + delta_ab0(i_qs,i_qr)*drds + delta_b0(i_qs)*dsds ) &
            * ( theta_b0(i_qr)*vrvr - theta_ab0(i_qs,i_qr)*vrvs + theta_b0(i_qs)*vsvs   &
            +  sigma_r + sigma_s )**0.5_rp
       pac(i_lracls2lg_r,k,i,j)= -0.25_rp*pi*e_sr*rhoq(i_ns,k,i,j)*rhoq(i_qr,k,i,j)*coef_acc_lrs_r
       pac(i_nracns2ng_r,k,i,j)= -0.25_rp*pi*e_sr*rhoq(i_ns,k,i,j)*rhoq(i_nr,k,i,j)*coef_acc_nrs_r
       !------------------------------------------------------------------------
       !
       !+++ pattern 3: a + a => b  (i-i)
       !
       !------------------------------------------------------------------------
       ! ice-ice ( reduction is double, but includes double-count)
       coef_acc_lii = &
              ( delta_b0(i_qi)*didi + delta_ab1(i_qi,i_qi)*didi + delta_b1(i_qi)*didi ) &
            * ( theta_b0(i_qi)*vivi - theta_ab1(i_qi,i_qi)*vivi + theta_b1(i_qi)*vivi   &
            +  sigma_i + sigma_i )**0.5_rp
       coef_acc_nii = &
              ( delta_b0(i_qi)*didi + delta_ab0(i_qi,i_qi)*didi + delta_b0(i_qi)*didi ) &
            * ( theta_b0(i_qi)*vivi - theta_ab0(i_qi,i_qi)*vivi + theta_b0(i_qi)*vivi   &
            +  sigma_i + sigma_i )**0.5_rp
       pac(i_liacli2ls,k,i,j)= -0.25_rp*pi*e_stick(k,i,j)*e_ii*rhoq(i_ni,k,i,j)*rhoq(i_qi,k,i,j)*coef_acc_lii
       pac(i_niacni2ns,k,i,j)= -0.25_rp*pi*e_stick(k,i,j)*e_ii*rhoq(i_ni,k,i,j)*rhoq(i_ni,k,i,j)*coef_acc_nii
       !
!          ci_aut(k,i,j)   =  0.25_RP*pi*E_ii*rhoq(I_NI,k,i,j)*coef_acc_LII
!          taui_aut(k,i,j) = 1._RP/max(E_stick(k,i,j)*ci_aut(k,i,j),1.E-10_RP)
!          tau_sce(k,i,j)  = rhoq(I_QI,k,i,j)/max(rhoq(I_QI,k,i,j)+rhoq(I_QS,k,i,j),1.E-10_RP)
       !------------------------------------------------------------------------
       !
       !+++ pattern 4: a + a => a  (s-s)
       !
       !------------------------------------------------------------------------
       ! snow-snow => snow
       coef_acc_nss = &
              ( delta_b0(i_qs)*dsds + delta_ab0(i_qs,i_qs)*dsds + delta_b0(i_qs)*dsds ) &
            * ( theta_b0(i_qs)*vsvs - theta_ab0(i_qs,i_qs)*vsvs + theta_b0(i_qs)*vsvs   &
            +  sigma_s + sigma_s )**0.5_rp
       pac(i_nsacns2ns,k,i,j)= -0.125_rp*pi*e_stick(k,i,j)*e_ss*rhoq(i_ns,k,i,j)*rhoq(i_ns,k,i,j)*coef_acc_nss
       !
       ! graupel-grauple => graupel
       coef_acc_ngg = &
              ( delta_b0(i_qg)*dgdg + delta_ab0(i_qg,i_qg)*dgdg + delta_b0(i_qg)*dgdg ) &
            * ( theta_b0(i_qg)*vgvg - theta_ab0(i_qg,i_qg)*vgvg + theta_b0(i_qg)*vgvg   &
            +  sigma_g + sigma_g )**0.5_rp
       pac(i_ngacng2ng,k,i,j)= -0.125_rp*pi*e_stick(k,i,j)*e_gg*rhoq(i_ng,k,i,j)*rhoq(i_ng,k,i,j)*coef_acc_ngg
       !
       !------------------------------------------------------------------------
       !--- Partial conversion
       ! SB06(70),(71)
       ! i_iconv2g: option whether partial conversions work or not
       ! ice-cloud => graupel
       sw = 0.5_rp - sign(0.5_rp,di_cri-ave_di) ! if( ave_di > di_cri )then sw=1
       wx_cri = cfill_i*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_di*ave_di*ave_di/xq(i_i,k,i,j) - 1.0_rp ) * sw
       pq(i_licon,k,i,j) = i_iconv2g * pac(i_liaclc2li,k,i,j)/max(1.0_rp, wx_cri) * sw
       pq(i_nicon,k,i,j) = i_iconv2g * pq(i_licon,k,i,j)/xq(i_i,k,i,j) * sw

       ! snow-cloud => graupel
       wx_crs = cfill_s*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_ds*ave_ds*ave_ds/xq(i_s,k,i,j) - 1.0_rp )
       pq(i_lscon,k,i,j) = i_sconv2g * (pac(i_lsaclc2ls,k,i,j))/max(1.0_rp, wx_crs)
       pq(i_nscon,k,i,j) = i_sconv2g * pq(i_lscon,k,i,j)/xq(i_s,k,i,j)
       !------------------------------------------------------------------------
       !--- enhanced melting( due to collection-freezing of water droplets )
       !    originally from Rutledge and Hobbs(1984). eq.(A.21)
       ! if T > 273.15 then temc_p=T-273.15, else temc_p=0
       ! 08/05/08 [fix] T.Mitsui LHF00 => LHF0
       ! melting occurs around T=273K, so LHF0 is suitable both SIMPLE and EXACT,
       ! otherwise LHF can have sign both negative(EXACT) and positive(SIMPLE).
!!$       coef_emelt   = -CL/LHF00*temc_p
       coef_emelt   =  cl/lhf0*temc_p
       ! cloud-graupel
       pq(i_lgacm,k,i,j) =  coef_emelt*pac(i_lgaclc2lg,k,i,j)
       pq(i_ngacm,k,i,j) =  pq(i_lgacm,k,i,j)/xq(i_g,k,i,j)
       ! rain-graupel
       pq(i_lgarm,k,i,j) =  coef_emelt*pac(i_lraclg2lg,k,i,j)
       pq(i_ngarm,k,i,j) =  pq(i_lgarm,k,i,j)/xq(i_g,k,i,j)
       ! cloud-snow
       pq(i_lsacm,k,i,j) =  coef_emelt*(pac(i_lsaclc2ls,k,i,j))
       pq(i_nsacm,k,i,j) =  pq(i_lsacm,k,i,j)/xq(i_s,k,i,j)
       ! rain-snow
       pq(i_lsarm,k,i,j) =  coef_emelt*(pac(i_lracls2lg_r,k,i,j)+pac(i_lracls2lg_s,k,i,j))
       pq(i_nsarm,k,i,j) =  pq(i_lsarm,k,i,j)/xq(i_g,k,i,j) ! collect? might be I_S
       ! cloud-ice
       pq(i_liacm,k,i,j) =  coef_emelt*pac(i_liaclc2li,k,i,j)
       pq(i_niacm,k,i,j) =  pq(i_liacm,k,i,j)/xq(i_i,k,i,j)
       ! rain-ice
       pq(i_liarm,k,i,j) =  coef_emelt*(pac(i_lracli2lg_r,k,i,j)+pac(i_lracli2lg_i,k,i,j))
       pq(i_niarm,k,i,j) =  pq(i_liarm,k,i,j)/xq(i_g,k,i,j) ! collect? might be I_I
    end do
    end do
    end do
    profile_stop("sn14_collection")

    !
    return
  end subroutine mixed_phase_collection_kij
  !----------------------------
  ! Auto-conversion, Accretion, Self-collection, Break-up
  subroutine aut_acc_slc_brk_kij(  &
       PQ,                     &
       rhoq, xq, dq_xave,      &
       rho                     )
    implicit none

    real(RP), intent(out) :: PQ(pq_max,ka,ia,ja)
    !
    real(RP), intent(in)  :: rhoq(qa,ka,ia,ja)
    real(RP), intent(in)  :: xq(5,ka,ia,ja)
    real(RP), intent(in)  :: dq_xave(5,ka,ia,ja)
    real(RP), intent(in)  :: rho(ka,ia,ja)
    !
    ! parameter for autoconversion
    real(RP), parameter :: kcc     = 4.44e+9_rp  ! collision efficiency [m3/kg2/sec]
    real(RP), parameter :: tau_min = 1.e-20_rp  ! empirical filter by T.Mitsui
    real(RP), parameter :: rx_sep  = 1.0_rp/x_sep ! 1/x_sep, 10/08/03 [Add] T.Mitsui
    !
    ! parameter for accretion
    real(RP), parameter :: kcr     = 5.8_rp     ! collision efficiency [m3/kg2/sec]
    real(RP), parameter :: thr_acc = 5.e-5_rp   ! threshold for universal function original
    !
    ! parameter for self collection and collison break-up
    real(RP), parameter :: krr     = 4.33_rp ! k_rr,      S08 (35)
    real(RP), parameter :: kaprr   = 60.7_rp  ! kappa_rr,  SB06(11)
    real(RP), parameter :: kbr     = 1000._rp ! k_br,      SB06(14)
    real(RP), parameter :: kapbr   = 2.3e+3_rp   ! kappa_br,  SB06(15)
    real(RP), parameter :: dr_min  = 0.35e-3_rp ! minimum diameter, SB06(13)-(15)
    !
    ! work variables
    real(RP) :: coef_nuc0 ! coefficient of number for Auto-conversion
    real(RP) :: coef_nuc1 !                mass
    real(RP) :: coef_aut0 !                number
    real(RP) :: coef_aut1 !                mass
    real(RP) :: lwc       ! lc+lr
    real(RP) :: tau       ! conversion ratio: qr/(qc+qr) ranges [0:1]
    real(RP) :: rho_fac   ! factor of air density
    real(RP) :: psi_aut   ! Universal function of Auto-conversion
    real(RP) :: psi_acc   ! Universal function of Accretion
    real(RP) :: psi_brk   ! Universal function of Breakup
    real(RP) :: ddr       ! diameter difference from equilibrium
    !
    integer :: i, j, k
    !
    coef_nuc0 = (nu(i_qc)+2.0_rp)/(nu(i_qc)+1.0_rp)
    coef_nuc1 = (nu(i_qc)+2.0_rp)*(nu(i_qc)+4.0_rp)/(nu(i_qc)+1.0_rp)/(nu(i_qc)+1.0_rp)
    coef_aut0 =  -kcc*coef_nuc0
    coef_aut1 =  -kcc/x_sep/20._rp*coef_nuc1
    !
    do j = js, je
    do i = is, ie
    do k = ks, ke
       lwc = rhoq(i_qr,k,i,j)+rhoq(i_qc,k,i,j)
       if( lwc > xc_min )then
          tau  = max(tau_min, rhoq(i_qr,k,i,j)/lwc)
       else
          tau  = tau_min
       end if
       rho_fac = sqrt(rho_0/max(rho(k,i,j),rho_min))
       !
       ! Auto-conversion ( cloud-cloud => rain )
       psi_aut       = 400._rp*(tau**0.7_rp)*(1.0_rp - (tau**0.7_rp))**3   ! (6) SB06
       pq(i_ncaut,k,i,j)  = coef_aut0*rhoq(i_qc,k,i,j)*rhoq(i_qc,k,i,j)*rho_fac*rho_fac    ! (9) SB06 sc+aut
       ! lc = lwc*(1-tau), lr = lwc*tau
       pq(i_lcaut,k,i,j)  = coef_aut1*lwc*lwc*xq(i_c,k,i,j)*xq(i_c,k,i,j) &          ! (4) SB06
            *((1.0_rp-tau)*(1.0_rp-tau) + psi_aut)*rho_fac*rho_fac        !
       pq(i_nraut,k,i,j)  = -rx_sep*pq(i_lcaut,k,i,j)                           ! (A7) SB01
       !
       ! Accretion ( cloud-rain => rain )
       psi_acc       =(tau/(tau+thr_acc))**4                          ! (8) SB06
       pq(i_lcacc,k,i,j)  = -kcr*rhoq(i_qc,k,i,j)*rhoq(i_qr,k,i,j)*rho_fac*psi_acc         ! (7) SB06
       pq(i_ncacc,k,i,j)  = -kcr*rhoq(i_nc,k,i,j)*rhoq(i_qr,k,i,j)*rho_fac*psi_acc         ! (A6) SB01
       !
       ! Self-collection ( rain-rain => rain )
       pq(i_nrslc,k,i,j)  = -krr*rhoq(i_nr,k,i,j)*rhoq(i_qr,k,i,j)*rho_fac                 ! (A.8) SB01
       !
       ! Collisional breakup of rain
       ddr           = min(1.e-3_rp, dq_xave(i_r,k,i,j) - dr_eq )
       if      (dq_xave(i_r,k,i,j) < dr_min )then                          ! negligible
          psi_brk      = -1.0_rp
          pq(i_nrbrk,k,i,j) = 0.0_rp
       else if (dq_xave(i_r,k,i,j) <= dr_eq  )then
          psi_brk      = kbr*ddr + 1.0_rp                               ! (14) SB06 (+1 is necessary)
          pq(i_nrbrk,k,i,j) = - (psi_brk + 1.0_rp)*pq(i_nrslc,k,i,j)              ! (13) SB06
       else
          psi_brk      = 2.0_rp*exp(kapbr*ddr) - 1.0_rp                   ! (15) SB06
          pq(i_nrbrk,k,i,j) = - (psi_brk + 1.0_rp)*pq(i_nrslc,k,i,j)              ! (13) SB06
       end if
       !
    end do
    end do
    end do
    !
    return
  end subroutine aut_acc_slc_brk_kij
  ! Vapor Deposition, Ice Melting
  subroutine dep_vapor_melt_ice_kij( &
       PQ,                   & ! out
       rho, tem, pre, qd,  & ! in
       rhoq,                   & ! in
       esw, esi,             & ! in
       xq, vt_xave, dq_xave  ) ! in
    use scale_const, only: &
       eps => const_eps
    implicit none

    ! Diffusion growth or Evaporation, Sublimation
    real(RP), intent(out) :: PQ(pq_max,ka,ia,ja)  ! mass change   for cloud, [Add]  09/08/18 T.Mitsui

    real(RP), intent(in)  :: rho(ka,ia,ja)     ! air density
    real(RP), intent(in)  :: tem(ka,ia,ja)     ! air temperature
    real(RP), intent(in)  :: pre(ka,ia,ja)     ! air pressure
    real(RP), intent(in)  :: qd(ka,ia,ja)      ! mixing ratio of dry air
    real(RP), intent(in)  :: esw(ka,ia,ja)     ! saturation vapor pressure(liquid water)
    real(RP), intent(in)  :: esi(ka,ia,ja)     ! saturation vapor pressure(solid water)
    real(RP), intent(in)  :: rhoq(qa,ka,ia,ja)
    real(RP), intent(in)  :: xq(5,ka,ia,ja)    ! mean mass
    ! Notice following values differ from mean terminal velocity or diameter.
    ! mean(vt(x)) /= vt(mean(x)) and mean(D(x)) /= D(mean(x))
    ! Following ones are vt(mean(x)) and D(mean(x)).
    real(RP), intent(in)  :: vt_xave(5,2,ka,ia,ja) ! terminal velocity of mean cloud 09/08/18 [Add], T.Mitsui
    !
    real(RP), intent(in)  :: dq_xave(5,ka,ia,ja) ! diameter
    !
    real(RP) :: rho_lim            ! limited density              09/08/18 T.Mitsui
    real(RP) :: temc_lim           ! limited temperature[celsius] 09/08/18 T.Mitsui
    real(RP) :: pre_lim            ! limited density              09/08/18 T.Mitsui
    real(RP) :: temc               ! temperature[celsius]
!    real(RP) :: pv                 ! vapor pressure
    real(RP) :: qv                 ! mixing ratio of water vapor [Add] 09/08/18
!    real(RP) :: ssw                ! super saturation ratio(liquid water)
!    real(RP) :: ssi                ! super saturation ratio(ice water)
    real(RP) :: nua, r_nua         ! kinematic viscosity of air
    real(RP) :: mua                ! viscosity of air
    real(RP) :: Kalfa              ! thermal conductance
    real(RP) :: Dw                 ! diffusivity of water vapor
    real(RP) :: Dt                 ! diffusivity of heat
    real(RP) :: Gw, Gi             ! diffusion factor by balance between heat and vapor
    real(RP) :: Gwr, Gii, Gis, Gig ! for rain, ice, snow and graupel.
    real(RP) :: Gm                 ! melting factor by balance between heat and vapor
    real(RP) :: Nsc_r3             !
    ! [Mod] 11/08/30 T.Mitsui, considering large and small branches
!    real(RP) :: Nrecs_r2            ! 09/08/18 [Add] T.Mitsui
    real(RP) :: Nrers_r2, Nreis_r2  !
    real(RP) :: Nress_r2, Nregs_r2  !
!    real(RP) :: Nrecl_r2            ! 09/08/18 [Add] T.Mitsui
    real(RP) :: Nrerl_r2, Nreil_r2  !
    real(RP) :: Nresl_r2, Nregl_r2  !
    real(RP) :: NscNrer_s, NscNrer_l
    real(RP) :: NscNrei_s, NscNrei_l
    real(RP) :: NscNres_s, NscNres_l
    real(RP) :: NscNreg_s, NscNreg_l
    real(RP) :: ventLR_s, ventLR_l
    real(RP) :: ventNI_s, ventNI_l, ventLI_s, ventLI_l
    real(RP) :: ventNS_s, ventNS_l, ventLS_s, ventLS_l
    real(RP) :: ventNG_s, ventNG_l, ventLG_s, ventLG_l
    !
    real(RP) :: wtr, wti, wts, wtg
    real(RP), parameter :: r_14=1.0_rp/1.4_rp
    real(RP), parameter :: r_15=1.0_rp/1.5_rp
    !
    real(RP) ::         ventLR
    real(RP) :: ventNI, ventLI
    real(RP) :: ventNS, ventLS
    real(RP) :: ventNG, ventLG
    !
    real(RP), parameter :: Re_max=1.e+3_rp
    real(RP), parameter :: Re_min=1.e-4_rp

    real(RP) :: sw
    !
    integer :: i, j, k
    !
    ! Notice,T.Mitsui
    ! Vapor deposition and melting would not be solved iteratively to reach equilibrium.
    ! Because following phenomena are not adjustment but transition.
    ! Just time-scales differ among them.
    ! If we would treat more appropreately, there would be time-splitting method to solve each ones.

    profile_start("sn14_dep_vapor")
    do j = js, je
    do i = is, ie
    do k = ks, ke
       temc    = tem(k,i,j) - t00   ! degC
       temc_lim= max(temc, -40._rp )       ! [Add] 09/08/18 T.Mitsui, Pruppacher and Klett(1997),(13-3)
       rho_lim = max(rho(k,i,j),rho_min)   ! [Add] 09/08/18 T.Mitsui
       qv      = rhoq(i_qv,k,i,j)/rho_lim
       pre_lim = rho_lim*(qd(k,i,j)*rdry + qv*rvap)*(temc_lim+t00) ![Add] 09/08/18 T.Mitsui
       !--------------------------------------------------------------------
       ! Diffusion growth part is described in detail
       ! by Pruppacher and Klett (1997) Sec. 13.2(liquid) and 13.3(solid)
       !
       ! G:factor of thermal diffusion(1st.term) and vapor diffusion(2nd. term)
       ! SB06(23),(38), Lin et al(31),(52) or others
       ! Dw is introduced by Pruppacher and Klett(1997),(13-3)
       dw      = 0.211e-4_rp* (((temc_lim+t00)/t00)**1.94_rp) *(p00/pre_lim)
       kalfa      = ka0  + temc_lim*dka_dt
       mua     = mua0 + temc_lim*dmua_dt
       nua     = mua/rho_lim
       r_nua   = 1.0_rp/nua
       gw      = (lhv0/kalfa/tem(k,i,j))*(lhv0/rvap/tem(k,i,j)-1.0_rp)+(rvap*tem(k,i,j)/dw/esw(k,i,j))
       gi      = (lhs0/kalfa/tem(k,i,j))*(lhs0/rvap/tem(k,i,j)-1.0_rp)+(rvap*tem(k,i,j)/dw/esi(k,i,j))
       ! capacities account for their surface geometries
       gwr     = 4.0_rp*pi/cap(i_qr)/gw
       gii     = 4.0_rp*pi/cap(i_qi)/gi
       gis     = 4.0_rp*pi/cap(i_qs)/gi
       gig     = 4.0_rp*pi/cap(i_qg)/gi
       ! vent: ventilation effect( asymmetry vapor field around particles due to aerodynamic )
       ! SB06 (30),(31) and each coefficient is by (88),(89)
       nsc_r3  = (nua/dw)**(0.33333333_rp)                    ! (Schmidt number )^(1/3)
       !
!       Nrecs_r2 = sqrt(max(Re_min,min(Re_max,vt_xave(I_C,1,k,i,j)*dq_xave(I_C,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) cloud
       nrers_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_r,1,k,i,j)*dq_xave(i_r,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) rain
       nreis_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_i,1,k,i,j)*dq_xave(i_i,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) cloud ice
       nress_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_s,1,k,i,j)*dq_xave(i_s,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) snow
       nregs_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_g,1,k,i,j)*dq_xave(i_g,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) graupel
       !
!       Nrecl_r2 = sqrt(max(Re_min,min(Re_max,vt_xave(I_C,2,k,i,j)*dq_xave(I_C,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) cloud
       nrerl_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_r,2,k,i,j)*dq_xave(i_r,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) rain
       nreil_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_i,2,k,i,j)*dq_xave(i_i,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) cloud ice
       nresl_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_s,2,k,i,j)*dq_xave(i_s,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) snow
       nregl_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_g,2,k,i,j)*dq_xave(i_g,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) graupel
       nscnrer_s=nsc_r3*nrers_r2 ! small rain
       nscnrer_l=nsc_r3*nrerl_r2 ! large rain
       !
       nscnrei_s=nsc_r3*nreis_r2 ! small ice
       nscnrei_l=nsc_r3*nreil_r2 ! large ice
       !
       nscnres_s=nsc_r3*nress_r2 ! small snow
       nscnres_l=nsc_r3*nresl_r2 ! large snow
       !
       nscnreg_s=nsc_r3*nregs_r2 ! small snow
       nscnreg_l=nsc_r3*nregl_r2 ! large snow
       !
       ventlr_s = ah_vent1(i_qr,1) + bh_vent1(i_qr,1)*nscnrer_s
       ventlr_l = ah_vent1(i_qr,2) + bh_vent1(i_qr,2)*nscnrer_l
       !
       ventni_s = ah_vent0(i_qi,1) + bh_vent0(i_qi,1)*nscnrei_s
       ventni_l = ah_vent0(i_qi,2) + bh_vent0(i_qi,2)*nscnrei_l
       ventli_s = ah_vent1(i_qi,1) + bh_vent1(i_qi,1)*nscnrei_s
       ventli_l = ah_vent1(i_qi,2) + bh_vent1(i_qi,2)*nscnrei_l
       !
       ventns_s = ah_vent0(i_qs,1) + bh_vent0(i_qs,1)*nscnres_s
       ventns_l = ah_vent0(i_qs,2) + bh_vent0(i_qs,2)*nscnres_l
       ventls_s = ah_vent1(i_qs,1) + bh_vent1(i_qs,1)*nscnres_s
       ventls_l = ah_vent1(i_qs,2) + bh_vent1(i_qs,2)*nscnres_l
       !
       ventng_s = ah_vent0(i_qg,1) + bh_vent0(i_qg,1)*nscnreg_s
       ventng_l = ah_vent0(i_qg,2) + bh_vent0(i_qg,2)*nscnreg_l
       ventlg_s = ah_vent1(i_qg,1) + bh_vent1(i_qg,1)*nscnreg_s
       ventlg_l = ah_vent1(i_qg,2) + bh_vent1(i_qg,2)*nscnreg_l
       !
       ! branch is 1.4 for rain, snow, graupel; is 1.0 for ice (PK97, 13-60,-61,-88,-89).
       !
       wtr     = ( min(max( nscnrer_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15 ! weighting between 1.4*0.5 and 1.4*2
       wti     = ( min(max( nscnrei_s     , 0.5_rp), 2.0_rp) -0.5_rp )*r_15 ! weighting between 1.0*0.5 and 1.0*2
       wts     = ( min(max( nscnres_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15 ! weighting between 1.4*0.5 and 1.4*2
       wtg     = ( min(max( nscnreg_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15 ! weighting between 1.4*0.5 and 1.4*2
       ! interpolation between two branches
       ventni  = (1.0_rp-wti)*ventni_s + wti*ventni_l
       ventns  = (1.0_rp-wts)*ventns_s + wts*ventns_l
       ventng  = (1.0_rp-wtg)*ventng_s + wtg*ventng_l
       !
       ventlr  = (1.0_rp-wtr)*ventlr_s + wtr*ventlr_l
       ventli  = (1.0_rp-wti)*ventli_s + wti*ventli_l
       ventls  = (1.0_rp-wts)*ventls_s + wts*ventls_l
       ventlg  = (1.0_rp-wtg)*ventlg_s + wtg*ventlg_l
       !
       ! SB06(29)
       ! [Mod] 08/05/08 T.Mitsui, recover PNXdep, and rain is only evaporation.
       ! Ni, Ns, Ng should decrease in nature so we add this term.
       ! And vapor deposition never occur unless number exist.
       ! [Add comment]  09/08/18
       ! recover condensation/evaporation of rain,
       ! and ventilation effects are not taken into account for cloud.
       !
!!$***************************************************************************
!!$        NOTICE:
!!$         09/08/18 [Mod] Hereafter PLxdep means inverse of timescale.
!!$***************************************************************************
!!$     PQ(I_LCdep,k,i,j) = Gwr*ssw*rhoq(I_NC,k,i,j)*dq_xave(I_C,k,i,j)*coef_deplc
!!$     PQ(I_LRdep,k,i,j) = Gwr*ssw*rhoq(I_NR,k,i,j)*dq_xave(I_R,k,i,j)*ventLR
!!$     PQ(I_LIdep,k,i,j) = Gii*ssi*rhoq(I_NI,k,i,j)*dq_xave(I_I,k,i,j)*ventLI
!!$     PQ(I_LSdep,k,i,j) = Gis*ssi*rhoq(I_NS,k,i,j)*dq_xave(I_S,k,i,j)*ventLS
!!$     PQ(I_LGdep,k,i,j) = Gig*ssi*rhoq(I_NG,k,i,j)*dq_xave(I_G,k,i,j)*ventLG
       pq(i_lcdep,k,i,j) = gwr*rhoq(i_nc,k,i,j)*dq_xave(i_c,k,i,j)*coef_deplc
       pq(i_lrdep,k,i,j) = gwr*rhoq(i_nr,k,i,j)*dq_xave(i_r,k,i,j)*ventlr
       pq(i_lidep,k,i,j) = gii*rhoq(i_ni,k,i,j)*dq_xave(i_i,k,i,j)*ventli
       pq(i_lsdep,k,i,j) = gis*rhoq(i_ns,k,i,j)*dq_xave(i_s,k,i,j)*ventls
       pq(i_lgdep,k,i,j) = gig*rhoq(i_ng,k,i,j)*dq_xave(i_g,k,i,j)*ventlg
       pq(i_nrdep,k,i,j) = pq(i_lrdep,k,i,j)/xq(i_r,k,i,j)
       pq(i_nidep,k,i,j) = 0.0_rp
       pq(i_nsdep,k,i,j) = pq(i_lsdep,k,i,j)/xq(i_s,k,i,j)
       pq(i_ngdep,k,i,j) = pq(i_lgdep,k,i,j)/xq(i_g,k,i,j)
       !
       !------------------------------------------------------------------------
       ! Melting part is described by Pruppacher and Klett (1997) Sec.16.3.1
       ! Here we omit "Shedding" of snow-flakes and ice-particles.
       ! "Shedding" may be applicative if you refer
       ! eq.(38) in Cotton etal.(1986) Jour. Clim. Appl. Meteor. p.1658-1680.
       ! SB06(73)
       dt      = kalfa/(cpvap*rho_0)
       ! Gm: factor caused by balance between
       !     "water evaporation cooling(1st.)" and "fusion heating(2nd.)"
       ! SB06(76)
       ! [fix] 08/05/08 T.Mitsui  LHF00 => EMELT  and  esw => PSAT0
       ! LHS0 is more suitable than LHS because melting occurs around 273.15 K.
       gm      = 2.0_rp*pi/emelt&
               * ( (kalfa*dt/dw)*(temc) + (dw*lhs0/rvap)*(esi(k,i,j)/tem(k,i,j)-psat0/t00) )
       ! SB06(76)
       ! Notice! melting only occurs where T > 273.15 K else doesn't.
       ! [fix] 08/05/08 T.Mitsui, Gm could be both positive and negative value.
       !       See Pruppacher and Klett(1997) eq.(16-79) or Rasmussen and Pruppacher(1982)
       sw = ( sign(0.5_rp,temc) + 0.5_rp ) * ( sign(0.5_rp,gm-eps) + 0.5_rp ) ! sw = 1 if( (temc>=0.0_RP) .AND. (Gm>0.0_RP) ), otherwise sw = 0
       !  if Gm==0 then rh and tem is critical value for melting process.
       ! 08/05/16 [Mod] T.Mitsui, change term of PLimlt. N_i => L_i/ (limited x_i)
       ! because melting never occur when N_i=0.
       pq(i_limlt,k,i,j) = - gm * rhoq(i_qi,k,i,j)*dq_xave(i_i,k,i,j)*ventli/xq(i_i,k,i,j) * sw
       ! [Mod] 08/08/23 T.Mitsui for Seifert(2008)
       pq(i_nimlt,k,i,j) = - gm * rhoq(i_ni,k,i,j)*dq_xave(i_i,k,i,j)*ventni/xq(i_i,k,i,j) * sw ! 09/08/18 [Mod] recover, T.Mitsui
       pq(i_lsmlt,k,i,j) = - gm * rhoq(i_qs,k,i,j)*dq_xave(i_s,k,i,j)*ventls/xq(i_s,k,i,j) * sw
       ! [Mod] 08/08/23 T.Mitsui for Seifert(2008)
       pq(i_nsmlt,k,i,j) = - gm * rhoq(i_ns,k,i,j)*dq_xave(i_s,k,i,j)*ventns/xq(i_s,k,i,j) * sw ! 09/08/18 [Mod] recover, T.Mitsui
       pq(i_lgmlt,k,i,j) = - gm * rhoq(i_qg,k,i,j)*dq_xave(i_g,k,i,j)*ventlg/xq(i_g,k,i,j) * sw
       ! [Mod] 08/08/23 T.Mitsui for Seifert(2008)
       pq(i_ngmlt,k,i,j) = - gm * rhoq(i_ng,k,i,j)*dq_xave(i_g,k,i,j)*ventng/xq(i_g,k,i,j) * sw ! 09/08/18 [Mod] recover, T.Mitsui

    end do
    end do
    end do
    profile_stop("sn14_dep_vapor")
    !
    return
  end subroutine dep_vapor_melt_ice_kij
  !-----------------------------------------------------------------------------
  subroutine freezing_water_kij( &
       dt,           &
       PQ,           &
       rhoq, xq, tem )
    !
    ! In this subroutine,
    ! We assumed surface temperature of droplets are same as environment.
    implicit none

    real(DP), intent(in) :: dt
    real(RP), intent(out):: PQ(pq_max,ka,ia,ja)
    !
    real(RP), intent(in) :: tem(ka,ia,ja)
    !
    real(RP), intent(in) :: rhoq(qa,ka,ia,ja)
    real(RP), intent(in) :: xq(5,ka,ia,ja)
    !
    real(RP), parameter :: temc_min = -65.0_rp
    real(RP), parameter :: a_het = 0.2_rp  ! SB06 (44)
    real(RP), parameter :: b_het = 0.65_rp ! SB06 (44)
    !
    real(RP) :: coef_m2_c
    real(RP) :: coef_m2_r
    ! temperature [celsius]
    real(RP) :: temc, temc2, temc3, temc4
    ! temperature function of homegenous/heterogenous freezing
    real(RP) :: Jhom, Jhet
    real(RP) :: rdt
    !
    integer :: i,j,k
    !
    rdt = 1.0_rp/dt
    !
    coef_m2_c =   coef_m2(i_qc)
    coef_m2_r =   coef_m2(i_qr)
    !
    profile_start("sn14_freezing")
    do j = js, je
    do i = is, ie
       do k = ks, ke
          temc = max( tem(k,i,j) - t00, temc_min )
          ! These cause from aerosol-droplet interaction.
          ! Bigg(1953) formula, Khain etal.(2000) eq.(4.5), Pruppacher and Klett(1997) eq.(9-48)
          jhet =  a_het*exp( -b_het*temc - 1.0_rp )
          ! These cause in nature.
          ! Cotton and Field 2002, QJRMS. (12)
          if( temc < -65.0_rp )then
             jhom = 10.0_rp**(24.37236_rp)*1.e+3_rp
             jhet =  a_het*exp( 65.0_rp*b_het - 1.0_rp ) ! 09/04/14 [Add], fixer T.Mitsui
          else if( temc < -30.0_rp ) then
             temc2 = temc*temc
             temc3 = temc*temc2
             temc4 = temc2*temc2
             jhom = 10.0_rp**(&
                  - 243.40_rp - 14.75_rp*temc - 0.307_rp*temc2 &
                  - 0.00287_rp*temc3 - 0.0000102_rp*temc4 ) *1.e+3_rp
          else if( temc < 0.0_rp) then
             jhom = 10._rp**(-7.63_rp-2.996_rp*(temc+30.0_rp))*1.e+3_rp
          else
             jhom = 0.0_rp
             jhet = 0.0_rp
          end if
          ! Note, xc should be limited in range[xc_min:xc_max].
          ! and PNChom need to be calculated by NC
          ! because reduction rate of Nc need to be bound by NC.
          ! For the same reason PLChom also bound by LC and xc.
          ! Basically L and N should be independent
          ! but average particle mass x should be in suitable range.
          ! Homogenous Freezing
          pq(i_lchom,k,i,j) = 0.0_rp
          pq(i_nchom,k,i,j) = 0.0_rp
          ! Heterogenous Freezing
          pq(i_lchet,k,i,j) = -rdt*rhoq(i_qc,k,i,j)*( 1.0_rp - exp( -coef_m2_c*xq(i_c,k,i,j)*(jhet+jhom)*dt ) )
          pq(i_nchet,k,i,j) = -rdt*rhoq(i_nc,k,i,j)*( 1.0_rp - exp( -          xq(i_c,k,i,j)*(jhet+jhom)*dt ) )
          pq(i_lrhet,k,i,j) = -rdt*rhoq(i_qr,k,i,j)*( 1.0_rp - exp( -coef_m2_r*xq(i_r,k,i,j)*(jhet+jhom)*dt ) )
          pq(i_nrhet,k,i,j) = -rdt*rhoq(i_nr,k,i,j)*( 1.0_rp - exp( -          xq(i_r,k,i,j)*(jhet+jhom)*dt ) )
       end do
    end do
    end do
    profile_stop("sn14_freezing")
    !
    return
  end subroutine freezing_water_kij
  !-----------------------------------------------------------------------------
  subroutine mp_terminal_velocity( &
      velw, &
      rhoq, &
      DENS, &
      temp, &
      pres  )
    use scale_const, only: &
       const_undef
    implicit none

    real(RP), intent(out) :: velw(ka,ia,ja,qa) ! terminal velocity of cloud mass
    real(RP), intent(in)  :: rhoq(qa,ka,ia,ja) ! rho * q
    real(RP), intent(in)  :: DENS(ka,ia,ja)    ! rho
    real(RP), intent(in)  :: temp(ka,ia,ja)    ! temperature
    real(RP), intent(in)  :: pres(ka,ia,ja)    ! pressure

    real(RP) :: xq       ! average mass of 1 particle( mass/number )

    real(RP) :: rhofac   ! density factor for terminal velocity( air friction )
    real(RP) :: rhofac_q(5)

    real(RP) :: rlambdar ! work for diagnosis of Rain DSD ( Seifert, 2008 )
    real(RP) :: mud_r
    real(RP) :: dq, dql  ! weigthed diameter.   Improved Rogers etal. (1993) formula by T.Mitsui


    real(RP) :: weight ! weighting coefficient for 2-branches is determined by ratio between 0.745mm and weighted diameter.  SB06 Table.1
    real(RP) :: velq_s ! terminal velocity for small branch of Rogers formula
    real(RP) :: velq_l ! terminal velocity for large branch of Rogers formula

    real(RP) :: tmp
    integer :: k, i, j, iq
    !---------------------------------------------------------------------------

    profile_start("sn14_terminal_vel")

    mud_r = 3.0_rp * nu(i_qr) + 2.0_rp

!OCL XFILL
    do j = js, je
    do i = is, ie
    do k = 1, ka
       velw(k,i,j,i_qv) = const_undef
    end do
    end do
    end do

    do j = js, je
    do i = is, ie
    do k = ks, ke
       rhofac = rho_0 / max( dens(k,i,j), rho_min )

       ! QC
       rhofac_q(i_c) = rhofac ** gamma_v(i_qc)
       xq = max( xqmin(i_qc), min( xqmax(i_qc), rhoq(i_qc,k,i,j) / ( rhoq(i_nc,k,i,j) + nqmin(i_qc) ) ) )

       velw(k,i,j,i_qc) = -rhofac_q(i_c) * coef_vt1(i_qc,1) * xq**beta_v(i_qc,1)
       ! NC
       velw(k,i,j,i_nc) = -rhofac_q(i_c) * coef_vt0(i_qc,1) * xq**beta_vn(i_qc,1)

       ! QR
       rhofac_q(i_r) = rhofac ** gamma_v(i_qr)
       xq = max( xqmin(i_qr), min( xqmax(i_qr), rhoq(i_qr,k,i,j) / ( rhoq(i_nr,k,i,j) + nqmin(i_qr) ) ) )

       rlambdar = a_m(i_qr) * xq**b_m(i_qr) &
                   * ( (mud_r+3.0_rp) * (mud_r+2.0_rp) * (mud_r+1.0_rp) )**(-0.333333333_rp)
       dq = ( 4.0_rp + mud_r ) * rlambdar ! D^(3)+mu weighted mean diameter
       dql = dq
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )

       velq_s = coef_vtr_ar2 * dq &
            * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar )**(-5.0_rp-mud_r) )
       velq_l = coef_vtr_ar1 - coef_vtr_br1 &
            * ( 1.0_rp + coef_vtr_cr1*rlambdar )**(-4.0_rp-mud_r)
       velw(k,i,j,i_qr) = -rhofac_q(i_r) &
            * ( velq_l * (          weight ) &
              + velq_s * ( 1.0_rp - weight ) )
       ! NR
       dq = ( 1.0_rp + mud_r ) * rlambdar
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )

       velq_s = coef_vtr_ar2 * dql &
            * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar )**(-2.0_rp-mud_r) )
       velq_l = coef_vtr_ar1 - coef_vtr_br1 &
            * ( 1.0_rp + coef_vtr_cr1*rlambdar )**(-1.0_rp-mud_r)
       velw(k,i,j,i_nr) = -rhofac_q(i_r) &
            * ( velq_l * (          weight ) &
              + velq_s * ( 1.0_rp - weight ) )

       ! QI
       rhofac_q(i_i) = ( pres(k,i,j)/pre0_vt )**a_pre0_vt * ( temp(k,i,j)/tem0_vt )**a_tem0_vt
       xq = max( xqmin(i_qi), min( xqmax(i_qi), rhoq(i_qi,k,i,j) / ( rhoq(i_ni,k,i,j) + nqmin(i_qi) ) ) )

       tmp = a_m(i_qi) * xq**b_m(i_qi)
       dq = coef_dave_l(i_qi) * tmp
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_li ) ) ) )

       velq_s = coef_vt1(i_qi,1) * xq**beta_v(i_qi,1)
       velq_l = coef_vt1(i_qi,2) * xq**beta_v(i_qi,2)
       velw(k,i,j,i_qi) = -rhofac_q(i_i) &
            * ( velq_l * (          weight ) &
              + velq_s * ( 1.0_rp - weight ) )
       ! NI
       dq = coef_dave_n(i_qi) * tmp
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ni ) ) ) )

       velq_s = coef_vt0(i_qi,1) * xq**beta_vn(i_qi,1)
       velq_l = coef_vt0(i_qi,2) * xq**beta_vn(i_qi,2)
       velw(k,i,j,i_ni) = -rhofac_q(i_i) &
            * ( velq_l * (          weight ) &
              + velq_s * ( 1.0_rp - weight ) )

       ! QS
       rhofac_q(i_s) = rhofac_q(i_i)
       xq = max( xqmin(i_qs), min( xqmax(i_qs), rhoq(i_qs,k,i,j) / ( rhoq(i_ns,k,i,j) + nqmin(i_qs) ) ) )

       tmp = a_m(i_qs) * xq**b_m(i_qs)
       dq = coef_dave_l(i_qs) * tmp
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ls ) ) ) )

       velq_s = coef_vt1(i_qs,1) * xq**beta_v(i_qs,1)
       velq_l = coef_vt1(i_qs,2) * xq**beta_v(i_qs,2)
       velw(k,i,j,i_qs) = -rhofac_q(i_s) &
            * ( velq_l * (          weight ) &
              + velq_s * ( 1.0_rp - weight ) )
       ! NS
       dq = coef_dave_n(i_qs) * tmp
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ns ) ) ) )

       velq_s = coef_vt0(i_qs,1) * xq**beta_vn(i_qs,1)
       velq_l = coef_vt0(i_qs,2) * xq**beta_vn(i_qs,2)
       velw(k,i,j,i_ns) = -rhofac_q(i_s) &
            * ( velq_l * (          weight ) &
              + velq_s * ( 1.0_rp - weight ) )

       ! QG
       rhofac_q(i_g) = rhofac_q(i_i)
       xq = max( xqmin(i_qg), min( xqmax(i_qg), rhoq(i_qg,k,i,j) / ( rhoq(i_ng,k,i,j) + nqmin(i_qg) ) ) )

       tmp = a_m(i_qg) * xq**b_m(i_qg)
       dq = coef_dave_l(i_qg) * tmp
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_lg ) ) ) )

       velq_s = coef_vt1(i_qg,1) * xq**beta_v(i_qg,1)
       velq_l = coef_vt1(i_qg,2) * xq**beta_v(i_qg,2)
       velw(k,i,j,i_qg) = -rhofac_q(i_g) &
            * ( velq_l * (          weight ) &
              + velq_s * ( 1.0_rp - weight ) )
       ! NG
       dq = coef_dave_n(i_qg) * tmp
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ng ) ) ) )

       velq_s = coef_vt0(i_qg,1) * xq**beta_vn(i_qg,1)
       velq_l = coef_vt0(i_qg,2) * xq**beta_vn(i_qg,2)
       velw(k,i,j,i_ng) = -rhofac_q(i_g) &
            * ( velq_l * (          weight ) &
              + velq_s * ( 1.0_rp - weight ) )
    enddo
    enddo
    enddo

    do iq = i_qc, i_ng
    do j = js, je
    do i = is, ie
       velw(1:ks-2,i,j,iq) = const_undef
       velw(ks-1,i,j,iq) = velw(ks,i,j,iq)
       velw(ke+1:ka,i,j,iq) = const_undef
    enddo
    enddo
    enddo

    profile_stop("sn14_terminal_vel")

    return
  end subroutine mp_terminal_velocity
  !----------------------------------------------------------------
  subroutine update_by_phase_change_kij(   &
       ntdiv, ntmax,        & ! in [Add] 10/08/03
       dt,                  & ! in
       gsgam2,              & ! in
       z,                   & ! in
       dz,                  & ! in
       velz,                & ! in
       dTdt_rad,            & ! in
       rho,                 & ! in
       rhoe,                & ! inout
       rhoq, q,             & ! inout
       tem, pre,            & ! inout
       cva,                 & ! out
       esw, esi, rhoq2,     & ! in
       PQ,                  & ! in
       qc_evaporate,        & ! in
       sl_PLCdep,           &
       sl_PLRdep, sl_PNRdep )
    use scale_atmos_thermodyn, only: &
       aq_cv, &
       aq_cp
    use scale_atmos_saturation, only: &
       moist_pres2qsat_liq  => atmos_saturation_pres2qsat_liq,  &
       moist_pres2qsat_ice  => atmos_saturation_pres2qsat_ice,  &
       moist_dqsw_dtem_rho  => atmos_saturation_dqsw_dtem_rho,  &
       moist_dqsi_dtem_rho  => atmos_saturation_dqsi_dtem_rho,  &
       moist_dqsw_dtem_dpre => atmos_saturation_dqsw_dtem_dpre, &
       moist_dqsi_dtem_dpre => atmos_saturation_dqsi_dtem_dpre
    implicit none

    integer, intent(in)    :: ntdiv               ! [Add] 10/08/03
    integer, intent(in)    :: ntmax               ! [Add] 10/08/03
    !
    real(DP), intent(in)    :: dt                 ! time step[s]
    real(RP), intent(in)    :: gsgam2(ka,ia,ja)   ! metric
    real(RP), intent(in)    :: z(ka)              ! altitude [m]
    real(RP), intent(in)    :: dz(ka)             ! altitude [m]
    real(RP), intent(in)    :: velz(ka,ia,ja)     ! vertical velocity @ half point[m/s]
    real(RP), intent(in)    :: dTdt_rad(ka,ia,ja) ! temperture tendency by radiation[K/s]
    real(RP), intent(in)    :: rho(ka,ia,ja)      ! density[kg/m3]
    real(RP), intent(inout) :: rhoe(ka,ia,ja)     ! internal energy[J/m3]
    real(RP), intent(inout) :: rhoq(qa,ka,ia,ja)  ! tracers[kg/m3]
    real(RP), intent(inout) :: q(ka,ia,ja,qa)     ! tracers mixing ratio[kg/kg]
    real(RP), intent(inout) :: tem(ka,ia,ja)      ! temperature[K]
    real(RP), intent(inout) :: pre(ka,ia,ja)      ! pressure[Pa]
    real(RP), intent(out)   :: cva(ka,ia,ja)      ! specific heat at constant volume
    real(RP), intent(in)    :: esw(ka,ia,ja)      ! saturated vapor pressure for liquid
    real(RP), intent(in)    :: esi(ka,ia,ja)      !                          for ice
    real(RP), intent(in)    :: rhoq2(qa,ka,ia,ja)
    !+++ tendency[kg/m3/s]
    real(RP), intent(inout) :: PQ(pq_max,ka,ia,ja)
    real(RP), intent(out)   :: qc_evaporate(ka,ia,ja)
    !+++ Column integrated tendency[kg/m2/s]
    real(RP), intent(inout) :: sl_PLCdep(ia,ja)
    real(RP), intent(inout) :: sl_PLRdep(ia,ja), sl_PNRdep(ia,ja)
    !
    real(RP) :: Rmoist
    !
    real(RP) :: xi                     ! mean mass of ice particles
    real(RP) :: rrho                   ! 1/rho
    real(RP) :: wtem(ka,ia,ja)         ! temperature[K]
    real(RP) :: qdry                   ! mixing ratio of dry air
    !
    real(RP) :: r_cva                  ! specific heat at constant volume
    real(RP) :: r_cpa                  ! specific heat at constant pressure
    real(RP) :: qsw(ka,ia,ja)          ! saturated mixing ratio for liquid
    real(RP) :: qsi(ka,ia,ja)          ! saturated mixing ratio for solid
    real(RP) :: dqswdtem_rho(ka,ia,ja) ! (dqsw/dtem)_rho
    real(RP) :: dqsidtem_rho(ka,ia,ja) ! (dqsi/dtem)_rho
    real(RP) :: dqswdtem_pre(ka,ia,ja) ! (dqsw/dtem)_pre
    real(RP) :: dqsidtem_pre(ka,ia,ja) ! (dqsi/dtem)_pre
    real(RP) :: dqswdpre_tem(ka,ia,ja) ! (dqsw/dpre)_tem
    real(RP) :: dqsidpre_tem(ka,ia,ja) ! (dqsi/dpre)_tem
    !
    real(RP) :: w                      ! vetical velocity[m/s]
    real(RP) :: Acnd                   ! Pdynliq + Bergeron-Findeisen
    real(RP) :: Adep                   ! Pdyndep + Bergeron-Findeisen
    real(RP) :: aliqliq, asolliq
    real(RP) :: aliqsol, asolsol
    real(RP) :: Pdynliq                ! production term of ssw by vertical motion
    real(RP) :: Pdynsol                ! production term of ssi by vertical motion
    real(RP) :: Pradliq                ! production term of ssw by radiation
    real(RP) :: Pradsol                ! production term of ssi by radiation
    real(RP) :: taucnd, r_taucnd       ! time scale of ssw change by MP
    real(RP) :: taudep, r_taudep       ! time scale of ssi change by MP
    real(RP) :: taucnd_c(ka,ia,ja), r_taucnd_c  ! by cloud
    real(RP) :: taucnd_r(ka,ia,ja), r_taucnd_r  ! by rain
    real(RP) :: taudep_i(ka,ia,ja), r_taudep_i  ! by ice
    real(RP) :: taudep_s(ka,ia,ja), r_taudep_s  ! by snow
    real(RP) :: taudep_g(ka,ia,ja), r_taudep_g  ! by graupel
    ! alternative tendency through changing ssw and ssi
    real(RP) :: PNCdep ! [Add] 11/08/30 T.Mitsui
    real(RP) :: PLR2NR, PLI2NI, PLS2NS, PLG2NG
    real(RP) :: coef_a_cnd, coef_b_cnd
    real(RP) :: coef_a_dep, coef_b_dep
    !
    real(RP) :: frz_dqc
    real(RP) :: frz_dnc
    real(RP) :: frz_dqr
    real(RP) :: frz_dnr
    real(RP) :: mlt_dqi
    real(RP) :: mlt_dni
    real(RP) :: mlt_dqs
    real(RP) :: mlt_dns
    real(RP) :: mlt_dqg
    real(RP) :: mlt_dng
    real(RP) :: dep_dqi
    real(RP) :: dep_dni
    real(RP) :: dep_dqs
    real(RP) :: dep_dns
    real(RP) :: dep_dqg
    real(RP) :: dep_dng
    real(RP) :: dep_dqr
    real(RP) :: dep_dnr
    real(RP) :: dep_dqc
    real(RP) :: dep_dnc   ! 11/08/30 [Add] T.Mitsui, dep_dnc
    real(RP) :: r_xc_ccn, r_xi_ccn ! 11/08/30 [Add] T.Mitsui
    !
    real(RP) :: drhoqv
    real(RP) :: drhoqc, drhoqr, drhoqi, drhoqs, drhoqg
    real(RP) :: drhonc, drhonr, drhoni, drhons, drhong
    !
    real(RP) :: fac1, fac2, fac3, fac4, fac5, fac6
    real(RP) :: r_rvaptem        ! 1/(Rvap*tem)
    real(RP) :: pv               ! vapor pressure
    real(RP) :: lvsw, lvsi       ! saturated vapor density
    real(RP) :: dlvsw, dlvsi     !
    ! [Add] 11/08/30 T.Mitsui
    real(RP) :: dcnd, ddep       ! total cndensation/deposition
    real(RP) :: uplim_cnd        ! upper limit of condensation
    real(RP) :: lowlim_cnd       ! lower limit of evaporation
    ! [Add] 11/08/30 T.Mitsui
    real(RP) :: uplim_dep        ! upper limit of condensation
    real(RP) :: lowlim_dep       ! lower limit of evaporation
    real(RP) :: ssw, ssi         ! supersaturation ratio
    real(RP) :: r_esw, r_esi     ! 1/esw, 1/esi
    real(RP) :: r_lvsw, r_lvsi   ! 1/(lvsw*ssw), 1/(lvsi*ssi)
    real(RP) :: r_dt             ! 1/dt
    real(RP) :: ssw_o, ssi_o
!    real(RP) :: dt_dyn
!    real(RP) :: dt_mp
    !
!    real(RP)       :: tem_lh(KA,IA,JA)
!    real(RP)       :: dtemdt_lh(KA,IA,JA)
    real(RP), save :: fac_cndc        = 1.0_rp
    logical, save :: opt_fix_taucnd_c=.false.
    logical, save :: flag_first      =.true.
    !
    namelist /nm_mp_sn14_condensation/ &
         opt_fix_taucnd_c, fac_cndc

    real(RP) :: fac_cndc_wrk
    !
    real(RP), parameter :: tau100day   = 1.e+7_rp
    real(RP), parameter :: r_tau100day = 1.e-7_rp
    real(RP), parameter :: eps=1.e-30_rp
    !
    integer :: i,j,k,iqw
    real(RP) :: sw
    !

    ! [Add] 11/08/30 T.Mitsui
    if( flag_first )then
       flag_first = .false.
       rewind(io_fid_conf)
       read  (io_fid_conf,nml=nm_mp_sn14_condensation, end=100)
100    if( io_l ) write (io_fid_log,nml=nm_mp_sn14_condensation)
    end if
    !
!    dt_dyn     = dt*ntmax
!    dt_mp      = dt*(ntdiv-1)
    !
    r_dt       = 1.0_rp/dt
    !
    r_xc_ccn=1.0_rp/xc_ccn
!    r_xi_ccn=1.0_RP/xi_ccn
    !
    if( opt_fix_taucnd_c )then
       fac_cndc_wrk = fac_cndc**(1.0_rp-b_m(i_qc))
       do j = js, je
       do i = is, ie
       do k = ks, ke
          pq(i_lcdep,k,i,j)  = pq(i_lcdep,k,i,j)*fac_cndc_wrk
       end do
       end do
       end do
       if( io_l ) write(io_fid_log,*) "taucnd:fac_cndc_wrk=",fac_cndc_wrk
    end if

!OCL XFILL
    do j = js, je
    do i = is, ie
    do k = ks, ke
       ! Temperature lower limit is only used for saturation condition.
       ! On the other hand original "tem" is used for calculation of latent heat or energy equation.
       wtem(k,i,j)  = max( tem(k,i,j), tem_min )
    end do
    end do
    end do

    call moist_pres2qsat_liq ( qsw, wtem, pre )
    call moist_pres2qsat_ice ( qsi, wtem, pre )
    call moist_dqsw_dtem_rho ( dqswdtem_rho, wtem, rho )
    call moist_dqsi_dtem_rho ( dqsidtem_rho, wtem, rho )
    call moist_dqsw_dtem_dpre( dqswdtem_pre, dqswdpre_tem, wtem, pre )
    call moist_dqsi_dtem_dpre( dqsidtem_pre, dqsidpre_tem, wtem, pre )

    profile_start("sn14_update")
    do j = js, je
    do i = is, ie
    do k = ks, ke
       if( z(k) <= 25000.0_rp )then
          w = 0.5_rp*(velz(k,i,j) + velz(k+1,i,j))
       else
          w = 0.0_rp
       end if
       if( pre(k,i,j) < esw(k,i,j)+1.e-10_rp )then
          qsw(k,i,j) = 1.0_rp
          dqswdtem_rho(k,i,j) = 0.0_rp
          dqswdtem_pre(k,i,j) = 0.0_rp
          dqswdpre_tem(k,i,j) = 0.0_rp
       end if
       if( pre(k,i,j) < esi(k,i,j)+1.e-10_rp )then
          qsi(k,i,j) = 1.0_rp
          dqsidtem_rho(k,i,j) = 0.0_rp
          dqsidtem_pre(k,i,j) = 0.0_rp
          dqsidpre_tem(k,i,j) = 0.0_rp
       end if

       r_rvaptem        = 1.0_rp/(rvap*wtem(k,i,j))
       lvsw             = esw(k,i,j)*r_rvaptem        ! rho=p/(Rv*T)
       lvsi             = esi(k,i,j)*r_rvaptem        !
       pv               = rhoq2(i_qv,k,i,j)*rvap*tem(k,i,j)
       r_esw            = 1.0_rp/esw(k,i,j)
       r_esi            = 1.0_rp/esi(k,i,j)
       ssw              = min( mp_ssw_lim, ( pv*r_esw-1.0_rp ) )
       ssi              = pv*r_esi - 1.0_rp
       r_lvsw           = 1.0_rp/lvsw
       r_lvsi           = 1.0_rp/lvsi
       r_taucnd_c       = pq(i_lcdep,k,i,j)*r_lvsw
       r_taucnd_r       = pq(i_lrdep,k,i,j)*r_lvsw
       r_taudep_i       = pq(i_lidep,k,i,j)*r_lvsi
       r_taudep_s       = pq(i_lsdep,k,i,j)*r_lvsi
       r_taudep_g       = pq(i_lgdep,k,i,j)*r_lvsi
!       taucnd_c(k,i,j)   = 1.0_RP/(r_taucnd_c+r_tau100day)
!       taucnd_r(k,i,j)   = 1.0_RP/(r_taucnd_r+r_tau100day)
!       taudep_i(k,i,j)   = 1.0_RP/(r_taudep_i+r_tau100day)
!       taudep_s(k,i,j)   = 1.0_RP/(r_taudep_s+r_tau100day)
!       taudep_g(k,i,j)   = 1.0_RP/(r_taudep_g+r_tau100day)

       calc_qdry( qdry, q, k, i, j, iqw )
       calc_cv( cva(k,i,j), qdry, q, k, i, j, iqw, cvdry, aq_cv )
       calc_r( rmoist, q(k,i,j,i_qv), qdry, rdry, rvap )
       r_cva = 1.0_rp / cva(k,i,j)
       r_cpa = 1.0_rp / (cva(k,i,j) + rmoist)

       ! Coefficient of latent heat release for ssw change by PLCdep and PLRdep
       aliqliq          = 1.0_rp &
               + r_cva*( lhv00              + (cvvap-cl)*tem(k,i,j) )*dqswdtem_rho(k,i,j)
       ! Coefficient of latent heat release for ssw change by PLIdep, PLSdep and PLGdep
       asolliq          = 1.0_rp &
               + r_cva*( lhv00 + lhf00 + (cvvap-ci)*tem(k,i,j) )*dqswdtem_rho(k,i,j)
       ! Coefficient of latent heat release for ssi change by PLCdep and PLRdep
       aliqsol          = 1.0_rp &
               + r_cva*( lhv00              + (cvvap-cl)*tem(k,i,j) )*dqsidtem_rho(k,i,j)
       ! Coefficient of latent heat release for ssi change by PLIdep, PLSdep and PLGdep
       asolsol          = 1.0_rp &
               + r_cva*( lhv00 + lhf00 + (cvvap-ci)*tem(k,i,j) )*dqsidtem_rho(k,i,j)
       pdynliq          = w * grav * ( r_cpa*dqswdtem_pre(k,i,j) + rho(k,i,j)*dqswdpre_tem(k,i,j) )
       pdynsol          = w * grav * ( r_cpa*dqsidtem_pre(k,i,j) + rho(k,i,j)*dqsidpre_tem(k,i,j) )
       pradliq          = -dtdt_rad(k,i,j)    * dqswdtem_rho(k,i,j)
       pradsol          = -dtdt_rad(k,i,j)    * dqsidtem_rho(k,i,j)

       ssw_o            = ssw
       ssi_o            = ssi
!       ssw_o            = ssw - Pdynliq*(dt_dyn-dt_mp)/qsw(k,i,j) + Pradliq*r_qsw*dt_mp
!       ssi_o            = ssi - Pdynsol*(dt_dyn-dt_mp)/qsi(k,i,j) + Pradsol*r_qsi*dt_mp

       acnd             = pdynliq + pradliq &
               - ( r_taudep_i+r_taudep_s+r_taudep_g ) * ( qsw(k,i,j) - qsi(k,i,j) )
       adep             = pdynsol + pradsol &
               + ( r_taucnd_c+r_taucnd_r )            * ( qsw(k,i,j) - qsi(k,i,j) )
       r_taucnd         = &
               + aliqliq*( r_taucnd_c+r_taucnd_r ) &
               + asolliq*( r_taudep_i+r_taudep_s+r_taudep_g )
       r_taudep         = &
               + aliqsol*( r_taucnd_c+r_taucnd_r )&
               + asolsol*( r_taudep_i+r_taudep_s+r_taudep_g )

       uplim_cnd        = max( rho(k,i,j)*ssw_o*qsw(k,i,j)*r_dt, 0.0_rp )
       lowlim_cnd       = min( rho(k,i,j)*ssw_o*qsw(k,i,j)*r_dt, 0.0_rp )
       if( r_taucnd < r_tau100day )then
!          taucnd            = tau100day
          pq(i_lcdep,k,i,j) = max(lowlim_cnd, min(uplim_cnd, pq(i_lcdep,k,i,j)*ssw_o ))
          pq(i_lrdep,k,i,j) = max(lowlim_cnd, min(uplim_cnd, pq(i_lrdep,k,i,j)*ssw_o ))
          pq(i_nrdep,k,i,j) = min(0.0_rp, pq(i_nrdep,k,i,j)*ssw_o )
!          PLR2NR           = 0.0_RP
       else
          taucnd     = 1.0_rp/r_taucnd
          ! Production term for liquid water content
          coef_a_cnd = rho(k,i,j)*acnd*taucnd
          coef_b_cnd = rho(k,i,j)*taucnd*r_dt*(ssw_o*qsw(k,i,j)-acnd*taucnd) * ( exp(-dt*r_taucnd) - 1.0_rp )
          pq(i_lcdep,k,i,j) = coef_a_cnd*r_taucnd_c - coef_b_cnd*r_taucnd_c
          plr2nr            = pq(i_nrdep,k,i,j)/(pq(i_lrdep,k,i,j)+1.e-30_rp)
          pq(i_lrdep,k,i,j) = coef_a_cnd*r_taucnd_r - coef_b_cnd*r_taucnd_r
          pq(i_nrdep,k,i,j) = min(0.0_rp, pq(i_lrdep,k,i,j)*plr2nr )
       end if

       uplim_dep        = max( rho(k,i,j)*ssi_o*qsi(k,i,j)*r_dt, 0.0_rp )
       lowlim_dep       = min( rho(k,i,j)*ssi_o*qsi(k,i,j)*r_dt, 0.0_rp )
       if( r_taudep < r_tau100day )then
!          taudep            = tau100day
          pq(i_lidep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lidep,k,i,j)*ssi_o ))
          pq(i_lsdep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lsdep,k,i,j)*ssi_o ))
          pq(i_lgdep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lgdep,k,i,j)*ssi_o ))
          pq(i_nidep,k,i,j) = min(0.0_rp, pq(i_nidep,k,i,j)*ssi_o )
          pq(i_nsdep,k,i,j) = min(0.0_rp, pq(i_nsdep,k,i,j)*ssi_o )
          pq(i_ngdep,k,i,j) = min(0.0_rp, pq(i_ngdep,k,i,j)*ssi_o )
       else
          taudep     = 1.0_rp/r_taudep
          ! Production term for ice water content
          coef_a_dep = rho(k,i,j)*adep*taudep
          coef_b_dep = rho(k,i,j)*taudep*r_dt*(ssi_o*qsi(k,i,j)-adep*taudep) * ( exp(-dt*r_taudep) - 1.0_rp )
          pli2ni           = pq(i_nidep,k,i,j)/max(pq(i_lidep,k,i,j),1.e-30_rp)
          pls2ns           = pq(i_nsdep,k,i,j)/max(pq(i_lsdep,k,i,j),1.e-30_rp)
          plg2ng           = pq(i_ngdep,k,i,j)/max(pq(i_lgdep,k,i,j),1.e-30_rp)
          pq(i_lidep,k,i,j) =  coef_a_dep*r_taudep_i - coef_b_dep*r_taudep_i
          pq(i_lsdep,k,i,j) =  coef_a_dep*r_taudep_s - coef_b_dep*r_taudep_s
          pq(i_lgdep,k,i,j) =  coef_a_dep*r_taudep_g - coef_b_dep*r_taudep_g
          pq(i_nidep,k,i,j) = min(0.0_rp, pq(i_lidep,k,i,j)*pli2ni )
          pq(i_nsdep,k,i,j) = min(0.0_rp, pq(i_lsdep,k,i,j)*pls2ns )
          pq(i_ngdep,k,i,j) = min(0.0_rp, pq(i_lgdep,k,i,j)*plg2ng )
       end if

       sw = 0.5_rp - sign(0.5_rp, pq(i_lcdep,k,i,j)+eps) != 1 for PLCdep<=-eps
       pncdep = min(0.0_rp, ((rhoq2(i_qc,k,i,j)+pq(i_lcdep,k,i,j)*dt)*r_xc_ccn - rhoq2(i_nc,k,i,j))*r_dt ) * sw
!       if( PQ(I_LCdep,k,i,j) < -eps )then
!          PNCdep = min(0.0_RP, ((rhoq2(I_QC,k,i,j)+PQ(I_LCdep,k,i,j)*dt)*r_xc_ccn - rhoq2(I_NC,k,i,j))*r_dt )
!       else
!          PNCdep = 0.0_RP
!       end if
!       if( PQ(I_LIdep,k,i,j) < -eps )then
!          PQ(I_NIdep,k,i,j) = min(0.0_RP, ((li(k,i,j)+PQ(I_LIdep,k,i,j)*dt)*r_xi_ccn - rhoq2(I_NI,k,i,j))*r_dt )
!       else
!          PQ(I_NIdep,k,i,j) = 0.0_RP
!       end if

       !--- evaporation/condensation
       r_rvaptem = 1.0_rp/(rvap*wtem(k,i,j))
       lvsw    = esw(k,i,j)*r_rvaptem
       dlvsw   = rhoq2(i_qv,k,i,j)-lvsw
       dcnd    = dt*(pq(i_lcdep,k,i,j)+pq(i_lrdep,k,i,j))

       sw = ( sign(0.5_rp,dcnd) + sign(0.5_rp,dlvsw) ) &
          * ( 0.5_rp + sign(0.5_rp,abs(dcnd)-eps) ) ! to avoid zero division
       ! sw= 1: always supersaturated
       ! sw=-1: always unsaturated
       ! sw= 0: partially unsaturated during timestep
       fac1 = min(dlvsw*sw,dcnd*sw)*sw / (abs(sw)-1.0_rp+dcnd) & ! sw=1,-1
            + 1.0_rp - abs(sw)                                   ! sw=0
       dep_dqc = max( dt*pq(i_lcdep,k,i,j)*fac1, &
                     -rhoq2(i_qc,k,i,j) - 1e30_rp*(sw+1.0_rp) )*abs(sw) != -lc for sw=-1, -inf for sw=1
       dep_dqr = max( dt*pq(i_lrdep,k,i,j)*fac1, &
                     -rhoq2(i_qr,k,i,j) - 1e30_rp*(sw+1.0_rp) )*abs(sw) != -lr for sw=-1, -inf for sw=1
!       if     ( (dcnd >  eps) .AND. (dlvsw > eps) )then
!          ! always supersaturated
!          fac1    = min(dlvsw,dcnd)/dcnd
!          dep_dqc =  dt*PQ(I_LCdep,k,i,j)*fac1
!          dep_dqr =  dt*PQ(I_LRdep,k,i,j)*fac1
!       else if( (dcnd < -eps) .AND. (dlvsw < -eps) )then
!          ! always unsaturated
!          fac1    = max( dlvsw,dcnd )/dcnd
!          dep_dqc = max( dt*PQ(I_LCdep,k,i,j)*fac1, -rhoq2(I_QC,k,i,j) )
!          dep_dqr = max( dt*PQ(I_LRdep,k,i,j)*fac1, -rhoq2(I_QR,k,i,j) )
!       else
!          ! partially unsaturated during timestep
!          fac1    = 1.0_RP
!          dep_dqc = 0.0_RP
!          dep_dqr = 0.0_RP
!       end if

       ! evaporation always lose number(always negative).
       dep_dnc = max( dt*pncdep*fac1, -rhoq2(i_nc,k,i,j) ) ! ss>0 dep=0, ss<0 dep<0 ! [Add] 11/08/30 T.Mitsui
       dep_dnr = max( dt*pq(i_nrdep,k,i,j)*fac1, -rhoq2(i_nr,k,i,j) ) ! ss>0 dep=0, ss<0 dep<0

       qc_evaporate(k,i,j) = - dep_dnc ! [Add] Y.Sato 15/09/08

       !--- deposition/sublimation
       lvsi    = esi(k,i,j)*r_rvaptem
       ddep    = dt*(pq(i_lidep,k,i,j)+pq(i_lsdep,k,i,j)+pq(i_lgdep,k,i,j))
       dlvsi   = rhoq2(i_qv,k,i,j)-lvsi  ! limiter for esi>1.d0

       sw = ( sign(0.5_rp,ddep) + sign(0.5_rp,dlvsi) ) &
          * ( 0.5_rp + sign(0.5_rp,abs(ddep)-eps) ) ! to avoid zero division
       ! sw= 1: always supersaturated
       ! sw=-1: always unsaturated
       ! sw= 0: partially unsaturated during timestep
       fac2 = min(dlvsi*sw,ddep*sw)*sw / (abs(sw)-1.0_rp+ddep) & ! sw=1,-1
            + 1.0_rp - abs(sw)                                   ! sw=0
       dep_dqi = max( dt*pq(i_lidep,k,i,j) &
                      * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), & != fac2 for sw=-1,1, 1 for sw=0
                     -rhoq2(i_qi,k,i,j) - 1e30_rp*(sw+1.0_rp) )      != -li for sw=-1, -inf for sw=0,1
       dep_dqs = max( dt*pq(i_lsdep,k,i,j) &
                      * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), & != fac2 for sw=-1,1, 1 for sw=0
                     -rhoq2(i_qs,k,i,j) - 1e30_rp*(sw+1.0_rp) )      != -ls for sw=-1, -inf for sw=0,1
       dep_dqg = max( dt*pq(i_lgdep,k,i,j) &
                      * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), & != fac2 for sw=-1,1, 1 for sw=0
                     -rhoq2(i_qg,k,i,j) - 1e30_rp*(sw+1.0_rp) )      != -lg for sw=-1, -inf for sw=0,1
!       if      ( (ddep >  eps) .AND. (dlvsi > eps) )then
!          ! always supersaturated
!          fac2    = min(dlvsi,ddep)/ddep
!          dep_dqi = dt*PQ(I_LIdep,k,i,j)*fac2
!          dep_dqs = dt*PQ(I_LSdep,k,i,j)*fac2
!          dep_dqg = dt*PQ(I_LGdep,k,i,j)*fac2
!       else if ( (ddep < -eps) .AND. (dlvsi < -eps) )then
!          ! always unsaturated
!          fac2    = max(dlvsi,ddep)/ddep
!          dep_dqi = max(dt*PQ(I_LIdep,k,i,j)*fac2, -rhoq2(I_QI,k,i,j) )
!          dep_dqs = max(dt*PQ(I_LSdep,k,i,j)*fac2, -rhoq2(I_QS,k,i,j) )
!          dep_dqg = max(dt*PQ(I_LGdep,k,i,j)*fac2, -rhoq2(I_QG,k,i,j) )
!       else
!          ! partially unsaturated during timestep
!          fac2    = 1.0_RP
!          dep_dqi = dt*PQ(I_LIdep,k,i,j)
!          dep_dqs = dt*PQ(I_LSdep,k,i,j)
!          dep_dqg = dt*PQ(I_LGdep,k,i,j)
!       end if

       ! evaporation always lose number(always negative).
       dep_dni = max( dt*pq(i_nidep,k,i,j)*fac2, -rhoq2(i_ni,k,i,j) ) ! ss>0 dep=0, ss<0 dep<0
       dep_dns = max( dt*pq(i_nsdep,k,i,j)*fac2, -rhoq2(i_ns,k,i,j) ) ! ss>0 dep=0, ss<0 dep<0
       dep_dng = max( dt*pq(i_ngdep,k,i,j)*fac2, -rhoq2(i_ng,k,i,j) ) ! ss>0 dep=0, ss<0 dep<0

       !--- freezing of cloud drop
       frz_dqc = max( dt*(pq(i_lchom,k,i,j)+pq(i_lchet,k,i,j)), -rhoq2(i_qc,k,i,j)-dep_dqc ) ! negative value
       frz_dnc = max( dt*(pq(i_nchom,k,i,j)+pq(i_nchet,k,i,j)), -rhoq2(i_nc,k,i,j)-dep_dnc ) ! negative value
       fac3    = ( frz_dqc-eps )/( dt*(pq(i_lchom,k,i,j)+pq(i_lchet,k,i,j))-eps )
       fac4    = ( frz_dnc-eps )/( dt*(pq(i_nchom,k,i,j)+pq(i_nchet,k,i,j))-eps )
       pq(i_lchom,k,i,j) = fac3*pq(i_lchom,k,i,j)
       pq(i_lchet,k,i,j) = fac3*pq(i_lchet,k,i,j)
       pq(i_nchom,k,i,j) = fac4*pq(i_nchom,k,i,j)
       pq(i_nchet,k,i,j) = fac4*pq(i_nchet,k,i,j)

       !--- melting
       ! ice change
       mlt_dqi = max( dt*pq(i_limlt,k,i,j), -rhoq2(i_qi,k,i,j)-dep_dqi )  ! negative value
       mlt_dni = max( dt*pq(i_nimlt,k,i,j), -rhoq2(i_ni,k,i,j)-dep_dni )  ! negative value
       ! snow change
       mlt_dqs = max( dt*pq(i_lsmlt,k,i,j), -rhoq2(i_qs,k,i,j)-dep_dqs )  ! negative value
       mlt_dns = max( dt*pq(i_nsmlt,k,i,j), -rhoq2(i_ns,k,i,j)-dep_dns )  ! negative value
       ! graupel change
       mlt_dqg = max( dt*pq(i_lgmlt,k,i,j), -rhoq2(i_qg,k,i,j)-dep_dqg )  ! negative value
       mlt_dng = max( dt*pq(i_ngmlt,k,i,j), -rhoq2(i_ng,k,i,j)-dep_dng )  ! negative value

       !--- freezing of larger droplets
       frz_dqr = max( dt*(pq(i_lrhet,k,i,j)), min(0.0_rp, -rhoq2(i_qr,k,i,j)-dep_dqr) ) ! negative value
       frz_dnr = max( dt*(pq(i_nrhet,k,i,j)), min(0.0_rp, -rhoq2(i_nr,k,i,j)-dep_dnr) ) ! negative value

       fac5         = ( frz_dqr-eps )/( dt*pq(i_lrhet,k,i,j)-eps )
       pq(i_lrhet,k,i,j) = fac5*pq(i_lrhet,k,i,j)
       fac6         = ( frz_dnr-eps )/( dt*pq(i_nrhet,k,i,j)-eps )
       pq(i_nrhet,k,i,j) = fac6*pq(i_nrhet,k,i,j)

       ! water vapor change
       drhoqv = -(dep_dqc+dep_dqi+dep_dqs+dep_dqg+dep_dqr)

       rhoq(i_qv,k,i,j) = max(0.0_rp, rhoq(i_qv,k,i,j) + drhoqv )

       rhoe(k,i,j) = rhoe(k,i,j) - lhv * drhoqv

       xi = min(xi_max, max(xi_min, rhoq2(i_qi,k,i,j)/(rhoq2(i_ni,k,i,j)+ni_min) ))
       sw = 0.5_rp + sign(0.5_rp,xi-x_sep) ! if (xi>=x_sep) then sw=1 else sw=0
                                                  ! sw=1: large ice crystals turn into rain by melting

       ! total cloud change
       drhoqc = ( frz_dqc - mlt_dqi*(1.0_rp-sw) + dep_dqc )
       drhonc = ( frz_dnc - mlt_dni*(1.0_rp-sw) + dep_dnc )
       ! total rain change
       drhoqr = ( frz_dqr - mlt_dqg - mlt_dqs - mlt_dqi*sw + dep_dqr )
       drhonr = ( frz_dnr - mlt_dng - mlt_dns - mlt_dni*sw + dep_dnr )

       rhoq(i_qc,k,i,j) = max(0.0_rp, rhoq(i_qc,k,i,j) + drhoqc )
       rhoq(i_nc,k,i,j) = max(0.0_rp, rhoq(i_nc,k,i,j) + drhonc )
       rhoq(i_qr,k,i,j) = max(0.0_rp, rhoq(i_qr,k,i,j) + drhoqr )
       rhoq(i_nr,k,i,j) = max(0.0_rp, rhoq(i_nr,k,i,j) + drhonr )

       ! total ice change
       drhoqi = (-frz_dqc + mlt_dqi             + dep_dqi )
       drhoni = (-frz_dnc + mlt_dni             + dep_dni )

       rhoq(i_qi,k,i,j) = max(0.0_rp, rhoq(i_qi,k,i,j) + drhoqi )
       rhoq(i_ni,k,i,j) = max(0.0_rp, rhoq(i_ni,k,i,j) + drhoni )

       rhoe(k,i,j) = rhoe(k,i,j) + lhf * drhoqi

       ! total snow change
       drhoqs = (           mlt_dqs             + dep_dqs )
       drhons = (           mlt_dns             + dep_dns )

       rhoq(i_qs,k,i,j) = max(0.0_rp, rhoq(i_qs,k,i,j) + drhoqs )
       rhoq(i_ns,k,i,j) = max(0.0_rp, rhoq(i_ns,k,i,j) + drhons )

       rhoe(k,i,j) = rhoe(k,i,j) + lhf * drhoqs

       ! total graupel change
       drhoqg = (-frz_dqr + mlt_dqg             + dep_dqg )
       drhong = (-frz_dnr + mlt_dng             + dep_dng )

       rhoq(i_qg,k,i,j) = max(0.0_rp, rhoq(i_qg,k,i,j) + drhoqg )
       rhoq(i_ng,k,i,j) = max(0.0_rp, rhoq(i_ng,k,i,j) + drhong )

       rhoe(k,i,j) = rhoe(k,i,j) + lhf * drhoqg

       !--- update mixing ratio
       rrho = 1.0_rp/rho(k,i,j)

       q(k,i,j,i_qv) = rhoq(i_qv,k,i,j) * rrho
       q(k,i,j,i_qc) = rhoq(i_qc,k,i,j) * rrho
       q(k,i,j,i_qr) = rhoq(i_qr,k,i,j) * rrho
       q(k,i,j,i_qi) = rhoq(i_qi,k,i,j) * rrho
       q(k,i,j,i_qs) = rhoq(i_qs,k,i,j) * rrho
       q(k,i,j,i_qg) = rhoq(i_qg,k,i,j) * rrho
       q(k,i,j,i_nc) = rhoq(i_nc,k,i,j) * rrho
       q(k,i,j,i_nr) = rhoq(i_nr,k,i,j) * rrho
       q(k,i,j,i_ni) = rhoq(i_ni,k,i,j) * rrho
       q(k,i,j,i_ns) = rhoq(i_ns,k,i,j) * rrho
       q(k,i,j,i_ng) = rhoq(i_ng,k,i,j) * rrho

       calc_qdry( qdry, q, k, i, j, iqw )
       calc_cv( cva(k,i,j), qdry, q, k, i, j, iqw, cvdry, aq_cv )
       calc_r( rmoist, q(k,i,j,i_qv), qdry, rdry, rvap )
       tem(k,i,j) = rhoe(k,i,j) / ( rho(k,i,j) * cva(k,i,j) )
       pre(k,i,j) = rho(k,i,j) * rmoist * tem(k,i,j)

       sl_plcdep(i,j) = sl_plcdep(i,j) + dep_dqc*dz(k)*gsgam2(k,i,j)
       sl_plrdep(i,j) = sl_plrdep(i,j) + dep_dqr*dz(k)*gsgam2(k,i,j)
       sl_pnrdep(i,j) = sl_pnrdep(i,j) + dep_dnr*dz(k)*gsgam2(k,i,j)
    end do
    end do
    end do
    profile_stop("sn14_update")

    return
  end subroutine update_by_phase_change_kij
  !-------------------------------------------------------------------------------
  subroutine mp_negativefilter( &
       DENS, &
       QTRC  )
    implicit none
    real(RP), intent(inout) :: DENS(ka,ia,ja)
    real(RP), intent(inout) :: QTRC(ka,ia,ja,qa)

    real(RP) :: diffq(ka,ia,ja)
    real(RP) :: r_xmin

    integer :: k, i, j, iq
    !---------------------------------------------------------------------------

    call prof_rapstart('MP_filter', 3)

    r_xmin = 1.0_rp / xmin_filter

    ! total hydrometeor (before correction)
    do j = js, je
    do i = is, ie

    do k = ks, ke
       diffq(k,i,j) = qtrc(k,i,j,i_qv) &
                    + qtrc(k,i,j,i_qc) &
                    + qtrc(k,i,j,i_qr) &
                    + qtrc(k,i,j,i_qi) &
                    + qtrc(k,i,j,i_qs) &
                    + qtrc(k,i,j,i_qg)
    enddo

    ! remove negative value of hydrometeor (mass, number)
    do iq = 1, qa
    do k  = ks, ke
       qtrc(k,i,j,iq) = max(0.0_rp, qtrc(k,i,j,iq))
    enddo
    enddo

    ! apply correction of hydrometeor to total density
    ! [note] mass conservation is broken here to fill rounding error.
    do k  = ks, ke
       dens(k,i,j) = dens(k,i,j)        &
                   * ( 1.0_rp           &
                     + qtrc(k,i,j,i_qv) &
                     + qtrc(k,i,j,i_qc) &
                     + qtrc(k,i,j,i_qr) &
                     + qtrc(k,i,j,i_qi) &
                     + qtrc(k,i,j,i_qs) &
                     + qtrc(k,i,j,i_qg) &
                     - diffq(k,i,j)      ) ! after-before
    enddo

    ! avoid unrealistical value of number concentration
    ! due to numerical diffusion in advection

    do k  = ks, ke
       if ( qtrc(k,i,j,i_nc) > qtrc(k,i,j,i_qc)*r_xmin ) then
          qtrc(k,i,j,i_nc) = qtrc(k,i,j,i_qc)*r_xmin
       endif
    enddo
    do k  = ks, ke
       if ( qtrc(k,i,j,i_nr) > qtrc(k,i,j,i_qr)*r_xmin ) then
          qtrc(k,i,j,i_nr) = qtrc(k,i,j,i_qr)*r_xmin
       endif
    enddo
    do k  = ks, ke
       if ( qtrc(k,i,j,i_ni) > qtrc(k,i,j,i_qi)*r_xmin ) then
          qtrc(k,i,j,i_ni) = qtrc(k,i,j,i_qi)*r_xmin
       endif
    enddo
    do k  = ks, ke
       if ( qtrc(k,i,j,i_ns) > qtrc(k,i,j,i_qs)*r_xmin ) then
          qtrc(k,i,j,i_ns) = qtrc(k,i,j,i_qs)*r_xmin
       endif
    enddo
    do k  = ks, ke
       if ( qtrc(k,i,j,i_ng) > qtrc(k,i,j,i_qg)*r_xmin ) then
          qtrc(k,i,j,i_ng) = qtrc(k,i,j,i_qg)*r_xmin
       endif
    enddo

    enddo
    enddo

    call prof_rapend('MP_filter', 3)

    return
  end subroutine mp_negativefilter
  !-------------------------------------------------------------------------------

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_sn14_cloudfraction( &
       cldfrac,       &
       QTRC,          &
       mask_criterion )
    use scale_grid_index
    use scale_tracer, only: &
       qad => qa
    implicit none

    real(RP), intent(out) :: cldfrac(ka,ia,ja)
    real(RP), intent(in)  :: QTRC   (ka,ia,ja,qad)
    real(RP), intent(in)  :: mask_criterion

    real(RP) :: qhydro
    integer  :: k, i, j, iq
    !---------------------------------------------------------------------------

    do j  = js, je
    do i  = is, ie
    do k  = ks, ke
       qhydro = 0.0_rp
       do iq = 1, mp_qa
          qhydro = qhydro + qtrc(k,i,j,i_mp2all(iq))
       enddo
       cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
    enddo
    enddo
    enddo

    return
  end subroutine atmos_phy_mp_sn14_cloudfraction

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_sn14_effectiveradius( &
       Re,    &
       QTRC0, &
       DENS0, &
       TEMP0  )
    use scale_grid_index
    use scale_tracer, only: &
       qad => qa, &
       mp_qad => mp_qa
    implicit none

    real(RP), intent(out) :: Re   (ka,ia,ja,mp_qad) ! effective radius          [cm]
    real(RP), intent(in)  :: QTRC0(ka,ia,ja,qad)    ! tracer mass concentration [kg/kg]
    real(RP), intent(in)  :: DENS0(ka,ia,ja)        ! density                   [kg/m3]
    real(RP), intent(in)  :: TEMP0(ka,ia,ja)        ! temperature               [K]

    ! mass concentration[kg/m3] and mean particle mass[kg]
    real(RP) :: xc(ka,ia,ja)
    real(RP) :: xr(ka,ia,ja)
    real(RP) :: xi(ka,ia,ja)
    real(RP) :: xs(ka,ia,ja)
    real(RP) :: xg(ka,ia,ja)
    ! diameter of average mass[kg/m3]
    real(RP) :: dc_ave(ka,ia,ja)
    real(RP) :: dr_ave(ka,ia,ja)
    ! radius of average mass
    real(RP) :: rc, rr
    ! 2nd. and 3rd. order moment of DSD
    real(RP) :: ri2m(ka,ia,ja), ri3m(ka,ia,ja)
    real(RP) :: rs2m(ka,ia,ja), rs3m(ka,ia,ja)
    real(RP) :: rg2m(ka,ia,ja), rg3m(ka,ia,ja)

    real(RP) :: coef_Fuetal1998
    ! r2m_min is minimum value(moment of 1 particle with 1 micron)
    real(RP), parameter :: r2m_min=1.e-12_rp
    real(RP), parameter :: um2cm = 100.0_rp

    real(RP) :: limitsw, zerosw
    integer :: k, i, j
    !---------------------------------------------------------------------------

    ! mean particle mass[kg]
    do j  = js, je
    do i  = is, ie
    do k  = ks, ke
       xc(k,i,j) = min( xc_max, max( xc_min, dens0(k,i,j)*qtrc0(k,i,j,i_qc)/(qtrc0(k,i,j,i_nc)+nc_min) ) )
       xr(k,i,j) = min( xr_max, max( xr_min, dens0(k,i,j)*qtrc0(k,i,j,i_qr)/(qtrc0(k,i,j,i_nr)+nr_min) ) )
       xi(k,i,j) = min( xi_max, max( xi_min, dens0(k,i,j)*qtrc0(k,i,j,i_qi)/(qtrc0(k,i,j,i_ni)+ni_min) ) )
       xs(k,i,j) = min( xs_max, max( xs_min, dens0(k,i,j)*qtrc0(k,i,j,i_qs)/(qtrc0(k,i,j,i_ns)+ns_min) ) )
       xg(k,i,j) = min( xg_max, max( xg_min, dens0(k,i,j)*qtrc0(k,i,j,i_qg)/(qtrc0(k,i,j,i_ng)+ng_min) ) )
    enddo
    enddo
    enddo

    ! diameter of average mass : SB06 eq.(32)
    do j = js, je
    do i = is, ie
    do k = ks, ke
       dc_ave(k,i,j) = a_m(i_qc) * xc(k,i,j)**b_m(i_qc)
       dr_ave(k,i,j) = a_m(i_qr) * xr(k,i,j)**b_m(i_qr)
    enddo
    enddo
    enddo

    ! cloud effective radius
    do j = js, je
    do i = is, ie
    do k = ks, ke
       rc = 0.5_rp * dc_ave(k,i,j)
       limitsw = 0.5_rp + sign(0.5_rp, rc-rmin_re )
       re(k,i,j,i_mp_qc) = coef_re(i_qc) * rc * limitsw * um2cm
    enddo
    enddo
    enddo

    ! rain effective radius
    do j = js, je
    do i = is, ie
    do k = ks, ke
       rr = 0.5_rp * dr_ave(k,i,j)
       limitsw = 0.5_rp + sign(0.5_rp, rr-rmin_re )
       re(k,i,j,i_mp_qr) = coef_re(i_qr) * rr * limitsw * um2cm
    enddo
    enddo
    enddo

    do j = js, je
    do i = is, ie
    do k = ks, ke
       ri2m(k,i,j) = pi * coef_rea2(i_qi) * qtrc0(k,i,j,i_ni) * a_rea2(i_qi) * xi(k,i,j)**b_rea2(i_qi)
       rs2m(k,i,j) = pi * coef_rea2(i_qs) * qtrc0(k,i,j,i_ns) * a_rea2(i_qs) * xs(k,i,j)**b_rea2(i_qs)
       rg2m(k,i,j) = pi * coef_rea2(i_qg) * qtrc0(k,i,j,i_ng) * a_rea2(i_qg) * xg(k,i,j)**b_rea2(i_qg)
    enddo
    enddo
    enddo

    ! Fu(1996), eq.(3.11) or Fu et al.(1998), eq.(2.5)
    coef_fuetal1998 = 3.0_rp / (4.0_rp*rhoi)
    do j = js, je
    do i = is, ie
    do k = ks, ke
       ri3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ni) * xi(k,i,j)
       rs3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ns) * xs(k,i,j)
       rg3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ng) * xg(k,i,j)
    enddo
    enddo
    enddo

    ! ice effective radius
    do j = js, je
    do i = is, ie
    do k = ks, ke
       zerosw = 0.5_rp - sign(0.5_rp, ri2m(k,i,j) - r2m_min )
       re(k,i,j,i_mp_qi) = ri3m(k,i,j) / ( ri2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
    enddo
    enddo
    enddo

    ! snow effective radius
    do j = js, je
    do i = is, ie
    do k = ks, ke
       zerosw = 0.5_rp - sign(0.5_rp, rs2m(k,i,j) - r2m_min )
       re(k,i,j,i_mp_qs) = rs3m(k,i,j) / ( rs2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
    enddo
    enddo
    enddo

    ! graupel effective radius
    do j = js, je
    do i = is, ie
    do k = ks, ke
       zerosw = 0.5_rp - sign(0.5_rp, rg2m(k,i,j) - r2m_min )
       re(k,i,j,i_mp_qg) = rg3m(k,i,j) / ( rg2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
    enddo
    enddo
    enddo

    return
  end subroutine atmos_phy_mp_sn14_effectiveradius
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_sn14_mixingratio( &
       Qe,    &
       QTRC0  )
    use scale_grid_index
    use scale_tracer, only: &
       qad => qa, &
       mp_qad => mp_qa
    implicit none

    real(RP), intent(out) :: Qe   (ka,ia,ja,mp_qad) ! mixing ratio of each cateory [kg/kg]
    real(RP), intent(in)  :: QTRC0(ka,ia,ja,qad)    ! tracer mass concentration [kg/kg]

    integer  :: ihydro
    !---------------------------------------------------------------------------

    do ihydro = 1, mp_qa
       qe(:,:,:,ihydro) = qtrc0(:,:,:,i_mp2all(ihydro))
    enddo

    return
  end subroutine atmos_phy_mp_sn14_mixingratio

end module scale_atmos_phy_mp_sn14