SCALE-RM
Functions/Subroutines | Variables
scale_atmos_phy_mp_sn14 Module Reference

module ATMOSPHERE / Physics Cloud Microphysics More...

Functions/Subroutines

subroutine, public atmos_phy_mp_sn14_config (MP_TYPE, QA, QS)
 Configure. More...
 
subroutine, public atmos_phy_mp_sn14_setup
 Setup Cloud Microphysics. More...
 
subroutine, public atmos_phy_mp_sn14 (DENS, MOMZ, MOMX, MOMY, RHOT, QTRC, CCN, EVAPORATE, SFLX_rain, SFLX_snow)
 Cloud Microphysics. More...
 
subroutine debug_tem_kij (point, tem, rho, pre, qv)
 
subroutine nucleation_kij (z, velz, rho, tem, pre, rhoq, PQ, cpa, dTdt_rad, qke, CCN, dt)
 
subroutine ice_multiplication_kij (PQ, Pac, tem, rhoq, xq)
 
subroutine mixed_phase_collection_kij (Pac, PQ, wtem, rhoq, xq, dq_xave, vt_xave)
 
subroutine aut_acc_slc_brk_kij (PQ, rhoq, xq, dq_xave, rho)
 
subroutine dep_vapor_melt_ice_kij (PQ, rho, tem, pre, qd, rhoq, esw, esi, xq, vt_xave, dq_xave)
 
subroutine freezing_water_kij (dt, PQ, rhoq, xq, tem)
 
subroutine update_by_phase_change_kij (ntdiv, ntmax, dt, gsgam2, z, dz, velz, dTdt_rad, rho, rhoe, rhoq, q, tem, pre, cva, esw, esi, rhoq2, PQ, qc_evaporate, sl_PLCdep, sl_PLRdep, sl_PNRdep)
 
subroutine mp_negativefilter (DENS, QTRC)
 
subroutine, public atmos_phy_mp_sn14_cloudfraction (cldfrac, QTRC, mask_criterion)
 Calculate Cloud Fraction. More...
 
subroutine, public atmos_phy_mp_sn14_effectiveradius (Re, QTRC0, DENS0, TEMP0)
 Calculate Effective Radius. More...
 
subroutine, public atmos_phy_mp_sn14_mixingratio (Qe, QTRC0)
 Calculate mixing ratio of each category. More...
 

Variables

integer, parameter, public qa_mp = 11
 
character(len=h_short), dimension(qa_mp), target, public atmos_phy_mp_sn14_name
 
character(len=h_mid), dimension(qa_mp), target, public atmos_phy_mp_sn14_desc
 
character(len=h_short), dimension(qa_mp), target, public atmos_phy_mp_sn14_unit
 
real(rp), dimension(n_hyd), target, public atmos_phy_mp_sn14_dens
 

Detailed Description

module ATMOSPHERE / Physics Cloud Microphysics

Description
Cloud Microphysics by 6 water category, double moment bulk scheme Seiki and Nakajima(2014) J. Atmos. Sci., 71, 833–853

Reference: – Journals Seifert and Beheng(2006) : Meteorol.Atmos.Phys.,vol.92,pp.45-66 Seifert and Beheng(2001) : Atmos.Res.,vol.59-60,pp.265-281 Seifert(2008) : J.Atmos.Sci.,vol.65,pp.3608-3619 Lin et al.(1983) : J.Appl.Meteor.,vol.22,pp.1065-1092 Ruttledge and Hobbs(1983) : J.Atmos.Sci.,vol.40,pp.1185-1206 Ruttledge and Hobbs(1984) : J.Atmos.Sci.,vol.40,pp.2949-2977 Cotton etal.(1986) : J.C.Appl.Meteor.,25,pp.1658-1680 Cotton and Field (2002) : QJRMS.,vol.128,pp2417-pp2437 Beard(1980) : J.Atmos.Sci.,vol.37,pp.1363-1374 [Add] 10/08/03 Berry and Reinhardt(1974a): J.Atmos.Sci.,vol.31,pp.1814-1824 Berry and Reinhardt(1974b): J.Atmos.Sci.,vol.31,pp.1825-1831 Fu(1996) : J.Climate, vol.9, pp.2058-2082 [Add] 10/08/03 Fu etal(1998) : J.Climate, vol.11, pp.2223-2237 [Add] 10/08/03 Ghan etal.(1997) : J.Geophys.Res.,vol.102,pp.21777-21794, [Add] 09/08/18 Hong et al.(2004) : Mon.Wea.Rev.,pp.103-120 Heymsfeild and Iaquinta(2000): J.Atmos.Sci., vol.57, pp.916-938 [Add] 10/08/03 Heymsfield and Kajikawa(1987): J.Atmos.Sci., vol.44, pp.1088-1099 Johnson(1981) : J.Atmos.Sci., vol.38, pp.215-218 [Add] 09/08/18 McFarquhar and Heymsfield(1996): J.Atmos.Sci.,vol.53,pp.2401-2423 Mitchell(1996) : J.Atmos.Sci., vol.53, pp.1710-1723. [Add] 10/08/03 Morrison etal.(2005) : Mon.Wea.Rev.,vol.62,pp.1665-1677, [Add] 09/08/18 Locatelli and Hobbs (1974): J.Geophys.Res., vol.70, pp.2185-2197 Lohmann(2002) : J.Atmos.Sci.,vol.59,pp.647-656 Takano and Liou(1989) : J.Atmos.Sci.,vol.46,pp.3-19 Takano and Liou(1994) : J.Atmos.Sci.,vol.52,pp.818-837 Auer and Veal(1970) : J.Atmos.Sci.,vol.27,pp.919-926 Ikawa et al.(1991) : J.M.S.J.,vol.69,pp.641-667 Murakami(1990) : J.M.S.J.,vol.68,pp.107-128 – Books Pruppacher and Klett(1997): Kluwer Academic Publishers Microphysics of Clouds and Precipitation, 2nd.edit. Seinfeld and Pandis(1998) : Wiley Interscience Atmospheric Chemistry and Physics. Jacobson(2005) : Cambridge press Fundamentals of Atmospheric Modeling, 2nd.edit.

Author
Team SCALE
History
  • 2011-10-24 (T.Seiki) [new] import from NICAM(11/08/30 ver.)
  • 2015-09-08 (Y.Sato) [add] Add evaporated cloud number concentration
NAMELIST
  • PARAM_ATMOS_PHY_MP
    nametypedefault valuecomment
    MP_DOAUTOCONVERSION logical .true.
    MP_DOPRECIPITATION logical .true.
    MP_SSW_LIM real(RP) 1.E+1_RP
    MP_COUPLE_AEROSOL logical .false. apply CCN effect?
    MP_NTMAX_SEDIMENTATION integer 1 10/08/03 [Add] T.Mitsui

  • nm_mp_sn14_init
    nametypedefault valuecomment
    OPT_DEBUG logical .false.
    OPT_DEBUG_TEM logical .false.
    OPT_DEBUG_INC logical .true.
    OPT_DEBUG_ACT logical .true.
    OPT_DEBUG_REE logical .true.
    OPT_DEBUG_BCS logical .true.
    NTMAX_PHASE_CHANGE integer 1
    NTMAX_COLLECTION integer 1

  • nm_mp_sn14_particles
    nametypedefault valuecomment
    A_M real(RP), dimension(HYDRO_MAX)
    B_M real(RP), dimension(HYDRO_MAX)
    ALPHA_V real(RP), dimension(HYDRO_MAX,2)
    BETA_V real(RP), dimension(HYDRO_MAX,2)
    GAMMA_V real(RP), dimension(HYDRO_MAX)
    ALPHA_VN real(RP), dimension(HYDRO_MAX,2)
    BETA_VN real(RP), dimension(HYDRO_MAX,2)
    A_AREA real(RP), dimension(HYDRO_MAX)
    B_AREA real(RP), dimension(HYDRO_MAX)
    CAP real(RP), dimension(HYDRO_MAX)
    NU real(RP), dimension(HYDRO_MAX)
    MU real(RP), dimension(HYDRO_MAX)
    OPT_M96_COLUMN_ICE logical .false.
    OPT_M96_ICE logical .true.
    AR_ICE_FIX real(RP) 0.7_RP

  • nm_mp_sn14_nucleation
    nametypedefault valuecomment
    IN_MAX real(RP) 1000.E+3_RP max num. of Ice-Nuclei [num/m3]
    C_CCN real(RP) 1.00E+8_RP
    KAPPA real(RP) 0.462_RP
    NM_M92 real(RP) 1.E+3_RP
    AM_M92 real(RP) -0.639_RP
    BM_M92 real(RP) 12.96_RP
    XC_CCN real(RP) 1.E-12_RP [kg]
    XI_CCN real(RP) 1.E-12_RP [kg] ! [move] 11/08/30 T.Mitsui
    TEM_CCN_LOW real(RP) 233.150_RP = -40 degC ! [Add] 10/08/03 T.Mitsui
    TEM_IN_LOW real(RP) 173.150_RP = -100 degC ! [Add] 10/08/03 T.Mitsui
    SSW_MAX real(RP) 1.1_RP [%]
    SSI_MAX real(RP) 0.60_RP
    NUCL_TWOMEY logical .false.
    INUCL_W logical .false.

  • nm_mp_sn14_collection
    nametypedefault valuecomment
    DC0 real(RP) 15.0E-6_RP lower threshold of cloud
    DC1 real(RP) 40.0E-6_RP upper threshold of cloud
    DI0 real(RP) 150.0E-6_RP lower threshold of cloud
    DS0 real(RP) 150.0E-6_RP lower threshold of cloud
    DG0 real(RP) 150.0E-6_RP lower threshold of cloud
    SIGMA_C real(RP) 0.00_RP cloud
    SIGMA_R real(RP) 0.00_RP rain
    SIGMA_I real(RP) 0.2_RP ice
    SIGMA_S real(RP) 0.2_RP snow
    SIGMA_G real(RP) 0.00_RP graupel
    OPT_STICK_KS96 logical .false.
    OPT_STICK_CO86 logical .false.
    E_IM real(RP) 0.80_RP ice max
    E_SM real(RP) 0.80_RP snow max
    E_GM real(RP) 1.00_RP graupel max
    E_IR real(RP) 1.0_RP ice x rain
    E_SR real(RP) 1.0_RP snow x rain
    E_GR real(RP) 1.0_RP graupel x rain
    E_II real(RP) 1.0_RP ice x ice
    E_SI real(RP) 1.0_RP snow x ice
    E_GI real(RP) 1.0_RP graupel x ice
    E_SS real(RP) 1.0_RP snow x snow
    E_GS real(RP) 1.0_RP graupel x snow
    E_GG real(RP) 1.0_RP graupel x graupel
    I_ICONV2G integer 1 ice => graupel
    I_SCONV2G integer 1 snow => graupel
    RHO_G real(RP) 900.0_RP [kg/m3]
    CFILL_I real(RP) 0.68_RP ice
    CFILL_S real(RP) 0.01_RP snow
    DI_CRI real(RP) 500.E-6_RP [m]

  • nm_mp_sn14_condensation
    nametypedefault valuecomment
    OPT_FIX_TAUCND_C logical .false.
    FAC_CNDC real(RP) 1.0_RP

History Output
namedescriptionunitvariable
pflux_QC precipitation flux of NC kg/m2/s FLX_hydro
pflux_QG precipitation flux of NG kg/m2/s FLX_hydro
pflux_QI precipitation flux of NI kg/m2/s FLX_hydro
pflux_QR precipitation flux of NR kg/m2/s FLX_hydro
pflux_QS precipitation flux of NS kg/m2/s FLX_hydro

Function/Subroutine Documentation

◆ atmos_phy_mp_sn14_config()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_config ( character(len=*), intent(in)  MP_TYPE,
integer, intent(out)  QA,
integer, intent(out)  QS 
)

Configure.

Definition at line 557 of file scale_atmos_phy_mp_sn14.F90.

References scale_atmos_hydrometeor::atmos_hydrometeor_regist(), atmos_phy_mp_sn14_desc, atmos_phy_mp_sn14_name, atmos_phy_mp_sn14_unit, scale_stdio::io_fid_log, scale_stdio::io_l, scale_process::prc_mpistop(), qa_mp, and scale_tracer::tracer_regist().

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config().

557  use scale_process, only: &
558  prc_mpistop
559  use scale_atmos_hydrometeor, only: &
560  atmos_hydrometeor_regist
561  use scale_tracer, only: &
562  tracer_regist
563  implicit none
564 
565  character(len=*), intent(in) :: MP_TYPE
566  integer, intent(out) :: QA
567  integer, intent(out) :: QS
568 
569  integer :: QS2
570  !---------------------------------------------------------------------------
571 
572  if( io_l ) write(io_fid_log,*)
573  if( io_l ) write(io_fid_log,*) '++++++ Module[Cloud Microphysics Tracer] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
574  if( io_l ) write(io_fid_log,*) '*** Tracers for Seiki and Nakajima (2014) 2-moment bulk 6 category'
575 
576  if ( mp_type /= 'SN14' ) then
577  write(*,*) 'xxx ATMOS_PHY_MP_TYPE is not SN14. Check!'
578  call prc_mpistop
579  end if
580 
581  call atmos_hydrometeor_regist( qs_mp, & ! [OUT]
582  1, 2, 3, & ! [IN]
583  atmos_phy_mp_sn14_name(1:6), & ! [IN]
584  atmos_phy_mp_sn14_desc(1:6), & ! [IN]
585  atmos_phy_mp_sn14_unit(1:6) ) ! [IN]
586 
587  call tracer_regist( qs2, & ! [OUT]
588  5, & ! [IN]
589  atmos_phy_mp_sn14_name(7:11), & ! [IN]
590  atmos_phy_mp_sn14_desc(7:11), & ! [IN]
591  atmos_phy_mp_sn14_unit(7:11) ) ! [IN]
592 
593  qa = qa_mp
594  qs = qs_mp
595  qe_mp = qs_mp + qa_mp - 1
596 
597  i_qv = qs
598  i_qc = qs + i_mp_qc
599  i_qr = qs + i_mp_qr
600  i_qi = qs + i_mp_qi
601  i_qs = qs + i_mp_qs
602  i_qg = qs + i_mp_qg
603  i_nc = qs + i_mp_nc
604  i_nr = qs + i_mp_nr
605  i_ni = qs + i_mp_ni
606  i_ns = qs + i_mp_ns
607  i_ng = qs + i_mp_ng
608 
609  return
scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_desc
character(len=h_mid), dimension(qa_mp), target, public atmos_phy_mp_sn14_desc
Definition: scale_atmos_phy_mp_sn14.F90:140
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_name
character(len=h_short), dimension(qa_mp), target, public atmos_phy_mp_sn14_name
Definition: scale_atmos_phy_mp_sn14.F90:139
scale_atmos_hydrometeor
module Hydrometeor
Definition: scale_atmos_sub_hydrometeor.F90:14
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_unit
character(len=h_short), dimension(qa_mp), target, public atmos_phy_mp_sn14_unit
Definition: scale_atmos_phy_mp_sn14.F90:141
scale_atmos_hydrometeor::atmos_hydrometeor_regist
subroutine, public atmos_hydrometeor_regist(Q0, NV, NL, NI, NAME, DESC, UNIT, ADVC)
Regist tracer.
Definition: scale_atmos_sub_hydrometeor.F90:196
scale_tracer::tracer_regist
subroutine, public tracer_regist(QS, NQ, NAME, DESC, UNIT, CV, CP, R, ADVC, MASS)
Regist tracer.
Definition: scale_tracer.F90:64
Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_setup()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_setup ( )

Setup Cloud Microphysics.

Definition at line 615 of file scale_atmos_phy_mp_sn14.F90.

References atmos_phy_mp_sn14_dens, scale_const::const_dice, scale_const::const_dwatr, scale_const::const_undef, scale_grid::grid_cdz, scale_atmos_hydrometeor::i_hc, scale_atmos_hydrometeor::i_hg, scale_atmos_hydrometeor::i_hi, scale_atmos_hydrometeor::i_hr, scale_atmos_hydrometeor::i_hs, scale_grid_index::ia, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_fid_nml, scale_stdio::io_l, scale_stdio::io_nml, scale_grid_index::ja, scale_grid_index::ka, scale_process::prc_mpistop(), and scale_time::time_dtsec_atmos_phy_mp.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config().

615  use scale_process, only: &
616  prc_mpistop
617  use scale_grid, only: &
618  cdz => grid_cdz
619  use scale_const, only: &
620  const_undef, &
621  const_dwatr, &
622  const_dice
623  use scale_time, only: &
624  time_dtsec_atmos_phy_mp
625  implicit none
626 
627  namelist / param_atmos_phy_mp / &
628  mp_doautoconversion, &
629  mp_doprecipitation, &
630  mp_ssw_lim, &
631  mp_couple_aerosol, &
632  mp_ntmax_sedimentation
633 
634  real(RP), parameter :: max_term_vel = 10.0_rp !-- terminal velocity for calculate dt of sedimentation
635  integer :: nstep_max
636  integer :: ierr
637  !---------------------------------------------------------------------------
638 
639  if( io_l ) write(io_fid_log,*)
640  if( io_l ) write(io_fid_log,*) '++++++ Module[Cloud Microphysics] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
641  if( io_l ) write(io_fid_log,*) '*** Seiki and Nakajima (2014) 2-moment bulk 6 category'
642 
643  !--- read namelist
644  rewind(io_fid_conf)
645  read(io_fid_conf,nml=param_atmos_phy_mp,iostat=ierr)
646  if( ierr < 0 ) then !--- missing
647  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
648  elseif( ierr > 0 ) then !--- fatal error
649  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP. Check!'
650  call prc_mpistop
651  endif
652  if( io_nml ) write(io_fid_nml,nml=param_atmos_phy_mp)
653 
654  atmos_phy_mp_sn14_dens(:) = const_undef
655  atmos_phy_mp_sn14_dens(i_hc) = const_dwatr
656  atmos_phy_mp_sn14_dens(i_hr) = const_dwatr
657  atmos_phy_mp_sn14_dens(i_hi) = const_dice
658  atmos_phy_mp_sn14_dens(i_hs) = const_dice
659  atmos_phy_mp_sn14_dens(i_hg) = const_dice
660 
661  wlabel(1) = "CLOUD"
662  wlabel(2) = "RAIN"
663  wlabel(3) = "ICE"
664  wlabel(4) = "SNOW"
665  wlabel(5) = "GRAUPEL"
666 
667  call mp_sn14_init
668 
669  allocate(nc_uplim_d(1,ia,ja))
670  nc_uplim_d(:,:,:) = 150.e6_rp
671 
672  nstep_max = int( ( time_dtsec_atmos_phy_mp * max_term_vel ) / minval( cdz ) )
673  mp_ntmax_sedimentation = max( mp_ntmax_sedimentation, nstep_max )
674 
675  mp_nstep_sedimentation = mp_ntmax_sedimentation
676  mp_rnstep_sedimentation = 1.0_rp / real(MP_ntmax_sedimentation,kind=rp)
677  mp_dtsec_sedimentation = time_dtsec_atmos_phy_mp * mp_rnstep_sedimentation
678 
679  if( io_l ) write(io_fid_log,*)
680  if( io_l ) write(io_fid_log,*) '*** Timestep of sedimentation is divided into : ', mp_ntmax_sedimentation, ' step(s)'
681  if( io_l ) write(io_fid_log,*) '*** DT of sedimentation [s] : ', mp_dtsec_sedimentation
682 
683  !--- For kij
684  allocate( gsgam2_d(ka,ia,ja) )
685  allocate( gsgam2h_d(ka,ia,ja) )
686  allocate( gam2_d(ka,ia,ja) )
687  allocate( gam2h_d(ka,ia,ja) )
688  allocate( rgsgam2_d(ka,ia,ja) )
689  allocate( rgs_d(ka,ia,ja) )
690  allocate( rgsh_d(ka,ia,ja) )
691  gsgam2_d(:,:,:) = 1.0_rp
692  gsgam2h_d(:,:,:) = 1.0_rp
693  gam2_d(:,:,:) = 1.0_rp
694  gam2h_d(:,:,:) = 1.0_rp
695  rgsgam2_d(:,:,:) = 1.0_rp
696  rgs_d(:,:,:) = 1.0_rp
697  rgsh_d(:,:,:) = 1.0_rp
698 
699  return
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_const::const_dwatr
real(rp), parameter, public const_dwatr
density of water [kg/m3]
Definition: scale_const.F90:84
scale_time::time_dtsec_atmos_phy_mp
real(dp), public time_dtsec_atmos_phy_mp
time interval of physics(microphysics) [sec]
Definition: scale_time.F90:41
scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_dens
real(rp), dimension(n_hyd), target, public atmos_phy_mp_sn14_dens
Definition: scale_atmos_phy_mp_sn14.F90:143
scale_const::const_dice
real(rp), parameter, public const_dice
density of ice [kg/m3]
Definition: scale_const.F90:85
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_time
module TIME
Definition: scale_time.F90:15
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_grid
module GRID (cartesian)
Definition: scale_grid_cartesian.F90:17
scale_grid::grid_cdz
real(rp), dimension(:), allocatable, public grid_cdz
z-length of control volume [m]
Definition: scale_grid_cartesian.F90:59
Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_mp_sn14()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14 ( real(rp), dimension (ka,ia,ja), intent(inout)  DENS,
real(rp), dimension (ka,ia,ja), intent(inout)  MOMZ,
real(rp), dimension (ka,ia,ja), intent(inout)  MOMX,
real(rp), dimension (ka,ia,ja), intent(inout)  MOMY,
real(rp), dimension (ka,ia,ja), intent(inout)  RHOT,
real(rp), dimension (ka,ia,ja,qa), intent(inout)  QTRC,
real(rp), dimension (ka,ia,ja), intent(in)  CCN,
real(rp), dimension(ka,ia,ja), intent(out)  EVAPORATE,
real(rp), dimension(ia,ja), intent(out)  SFLX_rain,
real(rp), dimension(ia,ja), intent(out)  SFLX_snow 
)

Cloud Microphysics.

Definition at line 716 of file scale_atmos_phy_mp_sn14.F90.

References scale_atmos_phy_mp_common::atmos_phy_mp_precipitation(), aut_acc_slc_brk_kij(), debug_tem_kij(), dep_vapor_melt_ice_kij(), scale_grid::dz, freezing_water_kij(), scale_grid::grid_cdz, scale_grid::grid_cz, scale_atmos_hydrometeor::i_nc, scale_atmos_hydrometeor::i_ng, scale_atmos_hydrometeor::i_ni, scale_atmos_hydrometeor::i_nr, scale_atmos_hydrometeor::i_ns, scale_atmos_hydrometeor::i_qc, scale_atmos_hydrometeor::i_qg, scale_atmos_hydrometeor::i_qi, scale_atmos_hydrometeor::i_qr, scale_atmos_hydrometeor::i_qs, scale_atmos_hydrometeor::i_qv, scale_grid_index::ia, ice_multiplication_kij(), scale_grid_index::ie, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_fid_nml, scale_stdio::io_l, scale_stdio::io_nml, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, mixed_phase_collection_kij(), mp_negativefilter(), nucleation_kij(), scale_process::prc_mpistop(), scale_prof::prof_rapend(), scale_prof::prof_rapstart(), scale_tracer::qa, qa_mp, scale_specfunc::sf_gamma(), scale_time::time_dtsec_atmos_phy_mp, scale_tracer::tracer_cv, scale_tracer::tracer_mass, scale_tracer::tracer_r, and update_by_phase_change_kij().

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config().

716  use scale_grid_index
717  use scale_tracer, only: &
718  qa
719  implicit none
720 
721  real(RP), intent(inout) :: DENS (KA,IA,JA)
722  real(RP), intent(inout) :: MOMZ (KA,IA,JA)
723  real(RP), intent(inout) :: MOMX (KA,IA,JA)
724  real(RP), intent(inout) :: MOMY (KA,IA,JA)
725  real(RP), intent(inout) :: RHOT (KA,IA,JA)
726  real(RP), intent(inout) :: QTRC (KA,IA,JA,QA)
727  real(RP), intent(in) :: CCN (KA,IA,JA)
728  real(RP), intent(out) :: EVAPORATE(KA,IA,JA) !--- number of evaporated cloud [/m3]
729  real(RP), intent(out) :: SFLX_rain(IA,JA)
730  real(RP), intent(out) :: SFLX_snow(IA,JA)
731  !---------------------------------------------------------------------------
732 
733  if( io_l ) write(io_fid_log,*) '*** Atmos physics step: Cloud microphysics(SN14)'
734 
735 #ifdef PROFILE_FIPP
736  call fipp_start()
737 #endif
738 
739  call mp_negativefilter( dens, qtrc )
740 
741  call mp_sn14( dens, & ! [INOUT]
742  momz, & ! [INOUT]
743  momx, & ! [INOUT]
744  momy, & ! [INOUT]
745  rhot, & ! [INOUT]
746  qtrc, & ! [INOUT]
747  ccn, & ! [IN]
748  evaporate, & ! [OUT]
749  sflx_rain, & ! [OUT]
750  sflx_snow ) ! [OUT]
751 
752  call mp_negativefilter( dens, qtrc )
753 
754 #ifdef PROFILE_FIPP
755  call fipp_stop()
756 #endif
757 
758  return
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
Here is the call graph for this function:
Here is the caller graph for this function:

◆ debug_tem_kij()

subroutine scale_atmos_phy_mp_sn14::debug_tem_kij ( integer, intent(in)  point,
real(rp), dimension(ka,ia,ja), intent(in)  tem,
real(rp), dimension(ka,ia,ja), intent(in)  rho,
real(rp), dimension(ka,ia,ja), intent(in)  pre,
real(rp), dimension (ka,ia,ja), intent(in)  qv 
)

Definition at line 2174 of file scale_atmos_phy_mp_sn14.F90.

References scale_grid_index::ie, scale_stdio::io_fid_log, scale_stdio::io_l, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, and scale_process::prc_myrank.

Referenced by atmos_phy_mp_sn14().

2174  use scale_process, only: &
2175  prc_myrank
2176  implicit none
2177 
2178  integer, intent(in) :: point
2179  real(RP), intent(in) :: tem(KA,IA,JA)
2180  real(RP), intent(in) :: rho(KA,IA,JA)
2181  real(RP), intent(in) :: pre(KA,IA,JA)
2182  real(RP), intent(in) :: qv (KA,IA,JA)
2183 
2184  integer :: k ,i, j
2185  !---------------------------------------------------------------------------
2186 
2187  do j = js, je
2188  do i = is, ie
2189  do k = ks, ke
2190  if ( tem(k,i,j) < tem_min &
2191  .OR. rho(k,i,j) < rho_min &
2192  .OR. pre(k,i,j) < 1.0_rp ) then
2193 
2194  if( io_l ) write(io_fid_log,'(A,I3,A,4(F16.5),3(I6))') &
2195  "*** 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
2196  endif
2197  enddo
2198  enddo
2199  enddo
2200 
2201  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_process::prc_myrank
integer, public prc_myrank
process num in local communicator
Definition: scale_process.F90:79
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
Here is the caller graph for this function:

◆ nucleation_kij()

subroutine scale_atmos_phy_mp_sn14::nucleation_kij ( real(rp), dimension(ka), intent(in)  z,
real(rp), dimension(ka,ia,ja), intent(in)  velz,
real(rp), dimension(ka,ia,ja), intent(in)  rho,
real(rp), dimension(ka,ia,ja), intent(in)  tem,
real(rp), dimension(ka,ia,ja), intent(in)  pre,
real(rp), dimension(i_qv:i_ng,ka,ia,ja), intent(in)  rhoq,
real(rp), dimension(pq_max,ka,ia,ja), intent(out)  PQ,
real(rp), dimension(ka,ia,ja), intent(in)  cpa,
real(rp), dimension(ka,ia,ja), intent(in)  dTdt_rad,
real(rp), dimension(ka,ia,ja), intent(in)  qke,
real(rp), dimension(ka,ia,ja), intent(in)  CCN,
real(rp), intent(in)  dt 
)

Definition at line 2214 of file scale_atmos_phy_mp_sn14.F90.

References scale_atmos_saturation::atmos_saturation_dqsi_dtem_rho(), scale_atmos_hydrometeor::i_nc, scale_atmos_hydrometeor::i_ni, scale_grid_index::ie, scale_stdio::io_fid_conf, scale_stdio::io_fid_nml, scale_stdio::io_nml, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, scale_process::prc_mpistop(), and scale_tracer::qa.

Referenced by atmos_phy_mp_sn14().

2214  use scale_process, only: &
2215  prc_mpistop
2216  use scale_tracer, only: &
2217  qa
2218  use scale_atmos_saturation, only: &
2219  moist_psat_liq => atmos_saturation_psat_liq, &
2220  moist_psat_ice => atmos_saturation_psat_ice, &
2221  moist_pres2qsat_liq => atmos_saturation_pres2qsat_liq, &
2222  moist_pres2qsat_ice => atmos_saturation_pres2qsat_ice, &
2223  moist_dqsi_dtem_rho => atmos_saturation_dqsi_dtem_rho
2224  implicit none
2225 
2226  real(RP), intent(in) :: z(KA) !
2227  real(RP), intent(in) :: velz(KA,IA,JA) ! w of half point
2228  real(RP), intent(in) :: rho(KA,IA,JA) ! [Add] 09/08/18 T.Mitsui
2229  real(RP), intent(in) :: tem(KA,IA,JA) ! [Add] 09/08/18 T.Mitsui
2230  real(RP), intent(in) :: pre(KA,IA,JA) ! [Add] 09/08/18 T.Mitsui
2231  !
2232  real(RP), intent(in) :: rhoq(I_QV:I_NG,KA,IA,JA) !
2233  real(RP), intent(out) :: PQ(PQ_MAX,KA,IA,JA)
2234  !
2235  real(RP), intent(in) :: cpa(KA,IA,JA) ! in 09/08/18 [Add] T.Mitsui
2236  real(RP), intent(in) :: dTdt_rad(KA,IA,JA) ! 09/08/18 T.Mitsui
2237  real(RP), intent(in) :: qke(KA,IA,JA) ! 09/08/18 T.Mitsui
2238  real(RP), intent(in) :: dt
2239  real(RP), intent(in) :: CCN(KA,IA,JA)
2240  !
2241  ! namelist variables
2242  !
2243  ! total aerosol number concentration [/m3]
2244  real(RP), parameter :: c_ccn_ocean= 1.00e+8_rp
2245  real(RP), parameter :: c_ccn_land = 1.26e+9_rp
2246  real(RP), save :: c_ccn = 1.00e+8_rp
2247  ! aerosol activation factor
2248  real(RP), parameter :: kappa_ocean= 0.462_rp
2249  real(RP), parameter :: kappa_land = 0.308_rp
2250  real(RP), save :: kappa = 0.462_rp
2251  real(RP), save :: c_in = 1.0_rp
2252  ! SB06 (36)
2253  real(RP), save :: nm_M92 = 1.e+3_rp
2254  real(RP), save :: am_M92 = -0.639_rp
2255  real(RP), save :: bm_M92 = 12.96_rp
2256  !
2257  real(RP), save :: in_max = 1000.e+3_rp ! max num. of Ice-Nuclei [num/m3]
2258  real(RP), save :: ssi_max= 0.60_rp
2259  real(RP), save :: ssw_max= 1.1_rp ! [%]
2260  !
2261  logical, save :: flag_first = .true.
2262  real(RP), save :: qke_min = 0.03_rp ! sigma=0.1[m/s], 09/08/18 T.Mitsui
2263  real(RP), save :: tem_ccn_low=233.150_rp ! = -40 degC ! [Add] 10/08/03 T.Mitsui
2264  real(RP), save :: tem_in_low =173.150_rp ! = -100 degC ! [Add] 10/08/03 T.Mitsui
2265  logical, save :: nucl_twomey = .false.
2266  logical, save :: inucl_w = .false.
2267  !
2268  namelist /nm_mp_sn14_nucleation/ &
2269  in_max, & !
2270  c_ccn, kappa, & ! cloud nucleation
2271  nm_m92, am_m92, bm_m92, & ! ice nucleation
2272  xc_ccn, xi_ccn, &
2273  tem_ccn_low, & ! [Add] 10/08/03 T.Mitsui
2274  tem_in_low, & ! [Add] 10/08/03 T.Mitsui
2275  ssw_max, ssi_max, &
2276  nucl_twomey, inucl_w ! [Add] 13/01/30 Y.Sato
2277  !
2278 ! real(RP) :: c_ccn_map(1,IA,JA) ! c_ccn horizontal distribution
2279 ! real(RP) :: kappa_map(1,IA,JA) ! kappa horizontal distribution
2280 ! real(RP) :: c_in_map(1,IA,JA) ! c_in horizontal distribution ! [Add] 11/08/30 T.Mitsui
2281  real(RP) :: esw(KA,IA,JA) ! saturation vapor pressure, water
2282  real(RP) :: esi(KA,IA,JA) ! ice
2283  real(RP) :: ssw(KA,IA,JA) ! super saturation (water)
2284  real(RP) :: ssi(KA,IA,JA) ! super saturation (ice)
2285 ! real(RP) :: w_dsswdz(KA,IA,JA) ! w*(d_ssw/ d_z) super saturation(water) flux
2286  real(RP) :: w_dssidz(KA,IA,JA) ! w*(d_ssi/ d_z), 09/04/14 T.Mitsui
2287 ! real(RP) :: ssw_below(KA,IA,JA)! ssw(k-1)
2288  real(RP) :: ssi_below(KA,IA,JA)! ssi(k-1), 09/04/14 T.Mitsui
2289  real(RP) :: z_below(KA,IA,JA) ! z(k-1)
2290  real(RP) :: dz ! z(k)-z(k-1)
2291  real(RP) :: pv ! vapor pressure
2292  ! work variables for Twomey Equation.
2293  real(RP) :: qsw(KA,IA,JA)
2294  real(RP) :: qsi(KA,IA,JA)
2295  real(RP) :: dqsidtem_rho(KA,IA,JA)
2296  real(RP) :: dssidt_rad(KA,IA,JA)
2297 ! real(RP) :: dni_ratio(KA,IA,JA)
2298  real(RP) :: wssi, wdssi
2299  !
2300 ! real(RP) :: xi_nuc(1,IA,JA) ! xi use the value @ cloud base
2301 ! real(RP) :: alpha_nuc(1,IA,JA) ! alpha_nuc
2302 ! real(RP) :: eta_nuc(1,IA,JA) ! xi use the value @ cloud base
2303  !
2304  real(RP) :: sigma_w(KA,IA,JA)
2305  real(RP) :: weff(KA,IA,JA)
2306  real(RP) :: weff_max(KA,IA,JA)
2307  !
2308  real(RP) :: coef_ccn(IA,JA)
2309  real(RP) :: slope_ccn(IA,JA)
2310  real(RP) :: nc_new(KA,IA,JA)
2311  real(RP) :: nc_new_below(KA,IA,JA)
2312  real(RP) :: dnc_new
2313  real(RP) :: nc_new_max ! Lohmann (2002),
2314  real(RP) :: a_max
2315  real(RP) :: b_max
2316  logical :: flag_nucleation(KA,IA,JA)
2317  !
2318  real(RP) :: r_gravity
2319  real(RP), parameter :: r_sqrt3=0.577350269_rp ! = sqrt(1.d0/3.d0)
2320  real(RP), parameter :: eps=1.e-30_rp
2321  !====> ! 09/08/18
2322  !
2323  real(RP) :: dlcdt_max, dli_max ! defined by supersaturation
2324  real(RP) :: dncdt_max, dni_max ! defined by supersaturation
2325  real(RP) :: rdt
2326  !
2327  integer :: i, j, k
2328  !
2329  if( flag_first )then
2330  rewind(io_fid_conf)
2331  read(io_fid_conf, nml=nm_mp_sn14_nucleation, end=100)
2332 100 if( io_nml ) write(io_fid_nml,nml=nm_mp_sn14_nucleation)
2333  flag_first=.false.
2334 
2335  if ( mp_couple_aerosol .AND. nucl_twomey ) then
2336  write(*,*) "xxx [mp_sn14/nucleation] nucl_twomey should be false when MP_couple_aerosol is true, stop"
2337  call prc_mpistop
2338  endif
2339  endif
2340  !
2341 ! c_ccn_map(1,:,:) = c_ccn
2342 ! kappa_map(1,:,:) = kappa
2343 ! c_in_map(1,:,:) = c_in
2344  !
2345 ! nc_uplim_d(1,:,:) = c_ccn_map(1,:,:)*1.5_RP
2346  do j = js, je
2347  do i = is, ie
2348  nc_uplim_d(1,i,j) = c_ccn*1.5_rp
2349  end do
2350  end do
2351  !
2352  rdt = 1.0_rp/dt
2353  r_gravity = 1.0_rp/grav
2354  !
2355  call moist_psat_liq ( esw, tem )
2356  call moist_psat_ice ( esi, tem )
2357  call moist_pres2qsat_liq( qsw, tem, pre )
2358  call moist_pres2qsat_ice( qsi, tem, pre )
2359  call moist_dqsi_dtem_rho( dqsidtem_rho, tem, rho )
2360  !
2361  ! Lohmann (2002),JAS, eq.(1) but changing unit [cm-3] => [m-3]
2362  a_max = 1.e+6_rp*0.1_rp*(1.e-6_rp)**1.27_rp
2363  b_max = 1.27_rp
2364  !
2365  ssi_max = 1.0_rp
2366 
2367  do j = js, je
2368  do i = is, ie
2369  do k = ks, ke
2370  pv = rhoq(i_qv,k,i,j)*rvap*tem(k,i,j)
2371  ssw(k,i,j) = min( mp_ssw_lim, ( pv/esw(k,i,j)-1.0_rp ) )*100.0_rp
2372  ssi(k,i,j) = (pv/esi(k,i,j) - 1.00_rp)
2373 ! ssw_below(k+1,i,j) = ssw(k,i,j)
2374  ssi_below(k+1,i,j) = ssi(k,i,j)
2375  z_below(k+1,i,j) = z(k)
2376  end do
2377 ! ssw_below(KS,i,j) = ssw(KS,i,j)
2378  ssi_below(ks,i,j) = ssi(ks,i,j)
2379  z_below(ks,i,j) = z(ks-1)
2380 
2381  ! dS/dz is evaluated by first order upstream difference
2382  !*** Solution for Twomey Equation ***
2383 ! 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)) * &
2384  coef_ccn(i,j) = 1.e+6_rp*0.88_rp*(c_ccn*1.e-6_rp)**(2.0_rp/(kappa + 2.0_rp)) &
2385 ! * (70.0_RP)**(kappa_map(1,i,j)/(kappa_map(1,i,j) + 2.0_RP))
2386  * (70.0_rp)**(kappa/(kappa + 2.0_rp))
2387 ! slope_ccn(i,j) = 1.5_RP*kappa_map(1,i,j)/(kappa_map(1,i,j) + 2.0_RP)
2388  slope_ccn(i,j) = 1.5_rp*kappa/(kappa + 2.0_rp)
2389  !
2390  do k=ks, ke
2391  sigma_w(k,i,j) = r_sqrt3*sqrt(max(qke(k,i,j),qke_min))
2392  end do
2393  sigma_w(ks-1,i,j) = sigma_w(ks,i,j)
2394  sigma_w(ke+1,i,j) = sigma_w(ke,i,j)
2395  ! effective vertical velocity
2396  do k=ks, ke
2397  weff(k,i,j) = 0.5_rp*(velz(k-1,i,j) + velz(k,i,j)) - cpa(k,i,j)*r_gravity*dtdt_rad(k,i,j)
2398  end do
2399 
2400  end do
2401  end do
2402  !
2403  if( mp_couple_aerosol ) then
2404 
2405  do j = js, je
2406  do i = is, ie
2407  do k = ks, ke
2408  if( ssw(k,i,j) > 1.e-10_rp .AND. pre(k,i,j) > 300.e+2_rp ) then
2409  nc_new(k,i,j) = max( ccn(k,i,j), c_ccn )
2410  else
2411  nc_new(k,i,j) = 0.0_rp
2412  endif
2413  enddo
2414  enddo
2415  enddo
2416 
2417  else
2418 
2419  if( nucl_twomey ) then
2420  ! diagnose cloud condensation nuclei
2421 
2422  do j = js, je
2423  do i = is, ie
2424  do k = ks, ke
2425  ! effective vertical velocity (maximum vertical velocity in turbulent flow)
2426  weff_max(k,i,j) = weff(k,i,j) + sigma_w(k,i,j)
2427  ! large scale upward motion region and saturated
2428  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
2429  ! Lohmann (2002), eq.(1)
2430  nc_new_max = coef_ccn(i,j)*weff_max(k,i,j)**slope_ccn(i,j)
2431  nc_new(k,i,j) = a_max*nc_new_max**b_max
2432  else
2433  nc_new(k,i,j) = 0.0_rp
2434  end if
2435  end do
2436  end do
2437  end do
2438  else
2439  ! calculate cloud condensation nuclei
2440  do j = js, je
2441  do i = is, ie
2442  do k = ks, ke
2443  if( ssw(k,i,j) > 1.e-10_rp .AND. pre(k,i,j) > 300.e+2_rp ) then
2444  nc_new(k,i,j) = c_ccn*ssw(k,i,j)**kappa
2445  else
2446  nc_new(k,i,j) = 0.0_rp
2447  endif
2448  enddo
2449  enddo
2450  enddo
2451  endif
2452 
2453  endif
2454 
2455  do j = js, je
2456  do i = is, ie
2457  do k = ks, ke
2458  ! nc_new is bound by upper limit
2459  if( nc_new(k,i,j) > nc_uplim_d(1,i,j) )then ! no more CCN
2460  flag_nucleation(k,i,j) = .false.
2461  nc_new_below(k+1,i,j) = 1.e+30_rp
2462  else if( nc_new(k,i,j) > eps )then ! nucleation can occur
2463  flag_nucleation(k,i,j) = .true.
2464  nc_new_below(k+1,i,j) = nc_new(k,i,j)
2465  else ! nucleation cannot occur(unsaturated or negative w)
2466  flag_nucleation(k,i,j) = .false.
2467  nc_new_below(k+1,i,j) = 0.0_rp
2468  end if
2469  end do
2470  nc_new_below(ks,i,j) = 0.0_rp
2471 ! do k=KS, KE
2472  ! search maximum value of nc_new
2473 ! if( ( nc_new(k,i,j) < nc_new_below(k,i,j) ) .OR. &
2474 ! ( nc_new_below(k,i,j) > c_ccn_map(1,i,j)*0.05_RP ) )then ! 5% of c_ccn
2475 ! ( nc_new_below(k,i,j) > c_ccn*0.05_RP ) )then ! 5% of c_ccn
2476 ! flag_nucleation(k,i,j) = .false.
2477 ! end if
2478 ! end do
2479 
2480  end do
2481  end do
2482 
2483  if( mp_couple_aerosol ) then
2484 
2485  do j = js, je
2486  do i = is, ie
2487  do k = ks, ke
2488  ! nucleation occurs at only cloud base.
2489  ! if CCN is more than below parcel, nucleation newly occurs
2490  ! effective vertical velocity
2491  if ( flag_nucleation(k,i,j) .AND. & ! large scale upward motion region and saturated
2492  tem(k,i,j) > tem_ccn_low ) then
2493  dlcdt_max = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
2494  dncdt_max = dlcdt_max/xc_min
2495 ! dnc_new = nc_new(k,i,j)-rhoq(I_NC,k,i,j)
2496  dnc_new = nc_new(k,i,j)
2497  pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
2498  pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
2499  else
2500  pq(i_ncccn,k,i,j) = 0.0_rp
2501  pq(i_lcccn,k,i,j) = 0.0_rp
2502  end if
2503  end do
2504  end do
2505  end do
2506 
2507  else
2508 
2509  if( nucl_twomey ) then
2510  do j = js, je
2511  do i = is, ie
2512  do k = ks, ke
2513  ! nucleation occurs at only cloud base.
2514  ! if CCN is more than below parcel, nucleation newly occurs
2515  ! effective vertical velocity
2516  if ( flag_nucleation(k,i,j) .AND. & ! large scale upward motion region and saturated
2517  tem(k,i,j) > tem_ccn_low .AND. &
2518  nc_new(k,i,j) > rhoq(i_nc,k,i,j) ) then
2519  dlcdt_max = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
2520  dncdt_max = dlcdt_max/xc_min
2521  dnc_new = nc_new(k,i,j)-rhoq(i_nc,k,i,j)
2522  pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
2523  pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
2524  else
2525  pq(i_ncccn,k,i,j) = 0.0_rp
2526  pq(i_lcccn,k,i,j) = 0.0_rp
2527  end if
2528  end do
2529  end do
2530  end do
2531  else
2532  do j = js, je
2533  do i = is, ie
2534  do k = ks, ke
2535  ! effective vertical velocity
2536  if( tem(k,i,j) > tem_ccn_low .AND. &
2537  nc_new(k,i,j) > rhoq(i_nc,k,i,j) ) then
2538  dlcdt_max = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
2539  dncdt_max = dlcdt_max/xc_min
2540  dnc_new = nc_new(k,i,j)-rhoq(i_nc,k,i,j)
2541  pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
2542  pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
2543  else
2544  pq(i_ncccn,k,i,j) = 0.0_rp
2545  pq(i_lcccn,k,i,j) = 0.0_rp
2546  end if
2547  end do
2548  end do
2549  end do
2550  endif
2551  endif
2552 
2553  !
2554  ! ice nucleation
2555  !
2556  ! +++ NOTE ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2557  ! Based on Phillips etal.(2006).
2558  ! However this approach doesn't diagnose Ni itself but diagnose tendency.
2559  ! Original approach adjust Ni instantaneously .
2560  ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2561  do j = js, je
2562  do i = is, ie
2563  do k = ks, ke
2564  dz = z(k) - z_below(k,i,j)
2565  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
2566  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)
2567  dli_max = (rhoq(i_qv,k,i,j) - esi(k,i,j)/(rvap*tem(k,i,j)))*rdt
2568  dni_max = min( dli_max/xi_ccn, (in_max-rhoq(i_ni,k,i,j))*rdt )
2569  wdssi = min( w_dssidz(k,i,j)+dssidt_rad(k,i,j), 0.01_rp)
2570  wssi = min( ssi(k,i,j), ssi_max)
2571  ! SB06(34),(35)
2572  if( ( wdssi > eps ) .AND. & !
2573  (tem(k,i,j) < 273.15_rp ) .AND. & !
2574  (rhoq(i_ni,k,i,j) < in_max ) .AND. &
2575  (wssi >= eps ) )then !
2576 ! 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)
2577  if( inucl_w ) then
2578  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)
2579  else
2580  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 )
2581  endif
2582  pq(i_liccn,k,i,j) = min(dli_max, pq(i_niccn,k,i,j)*xi_ccn )
2583  ! only for output
2584 ! dni_ratio(k,i,j) = dssidt_rad(k,i,j)/( w_dssidz(k,i,j)+dssidt_rad(k,i,j) )
2585  else
2586  pq(i_niccn,k,i,j) = 0.0_rp
2587  pq(i_liccn,k,i,j) = 0.0_rp
2588  end if
2589  end do
2590  end do
2591  end do
2592 
2593  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_atmos_saturation
module ATMOSPHERE / Saturation adjustment
Definition: scale_atmos_sub_saturation.F90:16
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_atmos_saturation::atmos_saturation_dqsi_dtem_rho
subroutine, public atmos_saturation_dqsi_dtem_rho(dqsdtem, temp, dens)
Definition: scale_atmos_sub_saturation.F90:1219
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ice_multiplication_kij()

subroutine scale_atmos_phy_mp_sn14::ice_multiplication_kij ( real(rp), dimension(pq_max,ka,ia,ja), intent(inout)  PQ,
real(rp), dimension(pac_max,ka,ia,ja), intent(in)  Pac,
real(rp), dimension(ka,ia,ja), intent(in)  tem,
real(rp), dimension(i_qv:i_ng,ka,ia,ja), intent(in)  rhoq,
real(rp), dimension(hydro_max,ka,ia,ja), intent(in)  xq 
)

Definition at line 2600 of file scale_atmos_phy_mp_sn14.F90.

References scale_atmos_hydrometeor::i_nc, scale_atmos_hydrometeor::i_qc, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, scale_tracer::qa, and scale_specfunc::sf_gamma().

Referenced by atmos_phy_mp_sn14().

2600 
2601  ! ice multiplication by splintering
2602  ! we consider Hallet-Mossop process
2603  use scale_specfunc, only: &
2604  gammafunc => sf_gamma
2605  use scale_tracer, only: &
2606  qa
2607  implicit none
2608  real(RP), intent(inout):: PQ(PQ_MAX,KA,IA,JA)
2609  !
2610  real(RP), intent(in) :: Pac(Pac_MAX,KA,IA,JA)
2611  real(RP), intent(in) :: tem(KA,IA,JA)
2612  real(RP), intent(in) :: rhoq(I_QV:I_NG,KA,IA,JA)
2613  real(RP), intent(in) :: xq(HYDRO_MAX,KA,IA,JA)
2614  !
2615  logical, save :: flag_first = .true.
2616  ! production of (350.d3 ice particles)/(cloud riming [g]) => 350*d6 [/kg]
2617  real(RP), parameter :: pice = 350.0e+6_rp
2618  ! production of (1 ice particle)/(250 cloud particles riming)
2619  real(RP), parameter :: pnc = 250.0_rp
2620  real(RP), parameter :: rc_cr= 12.e-6_rp ! critical size[micron]
2621  ! temperature factor
2622  real(RP) :: fp
2623  ! work for incomplete gamma function
2624  real(RP), save :: xc_cr ! mass[kg] of cloud with r=critical size[micron]
2625  real(RP), save :: alpha ! slope parameter of gamma function
2626  real(RP), save :: gm, lgm ! gamma(alpha), log(gamma(alpha))
2627  real(RP) :: igm ! in complete gamma(x,alpha)
2628  real(RP) :: x
2629  ! coefficient of expansion using in calculation of igm
2630  real(RP) :: a0,a1,a2,a3,a4,a5
2631  real(RP) :: a6,a7,a8,a9,a10
2632  real(RP) :: an1,an2,b0,b1,b2,c0,c1,c2
2633  real(RP) :: d0,d1,d2,e1,e2,h0,h1,h2
2634  real(RP), parameter :: eps=1.0e-30_rp
2635  ! number of cloud droplets larger than 12 micron(radius).
2636  real(RP) :: n12
2637  !
2638  real(RP) :: wn, wni, wns, wng
2639  integer :: i, j, k
2640  !
2641  if( flag_first )then
2642  flag_first = .false.
2643  ! work for Incomplete Gamma function
2644  xc_cr = (2.0_rp*rc_cr/a_m(i_mp_qc))**(1.0_rp/b_m(i_mp_qc))
2645  alpha = (nu(i_mp_qc)+1.0_rp)/mu(i_mp_qc)
2646  gm = gammafunc(alpha)
2647  lgm = log(gm)
2648  end if
2649  !
2650  do j = js, je
2651  do i = is, ie
2652  do k = ks, ke
2653  ! Here we assume particle temperature is same as environment temperature.
2654  ! If you want to treat in a better manner,
2655  ! you can diagnose with eq.(64) in CT(86)
2656  if (tem(k,i,j) > 270.16_rp)then
2657  fp = 0.0_rp
2658  else if(tem(k,i,j) >= 268.16_rp)then
2659  fp = (270.16_rp-tem(k,i,j))*0.5_rp
2660  else if(tem(k,i,j) >= 265.16_rp)then
2661  fp = (tem(k,i,j)-265.16_rp)*0.333333333_rp
2662  else
2663  fp = 0.0_rp
2664  end if
2665  ! Approximation of Incomplete Gamma function
2666  ! Here we calculate with algorithm by Numerical Recipes.
2667  ! This approach is based on eq.(78) in Cotton etal.(1986),
2668  ! but more accurate and expanded for Generalized Gamma Distribution.
2669  x = coef_lambda(i_mp_qc)*(xc_cr/xq(i_mp_qc,k,i,j))**mu(i_mp_qc)
2670  !
2671  if(x<1.e-2_rp*alpha)then ! negligible
2672  igm = 0.0_rp
2673  else if(x<alpha+1.0_rp)then ! series expansion
2674  ! 10th-truncation is enough for cloud droplet.
2675  a0 = 1.0_rp/alpha ! n=0
2676  a1 = a0*x/(alpha+1.0_rp) ! n=1
2677  a2 = a1*x/(alpha+2.0_rp) ! n=2
2678  a3 = a2*x/(alpha+3.0_rp) ! n=3
2679  a4 = a3*x/(alpha+4.0_rp) ! n=4
2680  a5 = a4*x/(alpha+5.0_rp) ! n=5
2681  a6 = a5*x/(alpha+6.0_rp) ! n=6
2682  a7 = a6*x/(alpha+7.0_rp) ! n=7
2683  a8 = a7*x/(alpha+8.0_rp) ! n=8
2684  a9 = a8*x/(alpha+9.0_rp) ! n=9
2685  a10 = a9*x/(alpha+10.0_rp) ! n=10
2686  igm = (a0+a1+a2+a3+a4+a5+a6+a7+a8+a9+a10)*exp( -x + alpha*log(x) - lgm )
2687  else if(x<alpha*1.d2) then ! continued fraction expansion
2688  ! 2nd-truncation is enough for cloud droplet.
2689  ! setup
2690  b0 = x+1.0_rp-alpha
2691  c0 = 1.0_rp/eps
2692  d0 = 1.0_rp/b0
2693  h0 = d0
2694  ! n=1
2695  an1 = -(1.0_rp-alpha)
2696  b1 = b0 + 2.0_rp
2697  d1 = 1.0_rp/(an1*d0+b1)
2698  c1 = b1+an1/c0
2699  e1 = d1*c1
2700  h1 = h0*e1
2701  ! n=2
2702  an2 = -2.0_rp*(2.0_rp-alpha)
2703  b2 = b1 + 2.0_rp
2704  d2 = 1.0_rp/(an2*d1+b2)
2705  c2 = b2+an2/c1
2706  e2 = d2*c2
2707  h2 = h1*e2
2708  !
2709  igm = 1.0_rp - exp( -x + alpha*log(x) - lgm )*h2
2710  else ! negligible
2711  igm = 1.0_rp
2712  end if
2713  ! n12 is number of cloud droplets larger than 12 micron.
2714  n12 = rhoq(i_nc,k,i,j)*(1.0_rp-igm)
2715  ! eq.(82) CT(86)
2716  wn = (pice + n12/((rhoq(i_qc,k,i,j)+xc_min)*pnc) )*fp ! filtered by xc_min
2717  wni = wn*(-pac(i_liaclc2li,k,i,j) ) ! riming production rate is all negative
2718  wns = wn*(-pac(i_lsaclc2ls,k,i,j) )
2719  wng = wn*(-pac(i_lgaclc2lg,k,i,j) )
2720  pq(i_nispl,k,i,j) = wni+wns+wng
2721  !
2722  pq(i_lsspl,k,i,j) = - wns*xq(i_mp_qi,k,i,j) ! snow => ice
2723  pq(i_lgspl,k,i,j) = - wng*xq(i_mp_qi,k,i,j) ! graupel => ice
2724  !
2725  end do
2726  end do
2727  end do
2728  !
2729  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_specfunc::sf_gamma
real(rp) function, public sf_gamma(x)
Gamma function.
Definition: scale_specfunc.F90:48
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_specfunc
module SPECFUNC
Definition: scale_specfunc.F90:13
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
Here is the call graph for this function:
Here is the caller graph for this function:

◆ mixed_phase_collection_kij()

subroutine scale_atmos_phy_mp_sn14::mixed_phase_collection_kij ( real(rp), dimension(pac_max,ka,ia,ja), intent(out)  Pac,
real(rp), dimension(pq_max,ka,ia,ja), intent(inout)  PQ,
real(rp), dimension(ka,ia,ja), intent(in)  wtem,
real(rp), dimension(i_qv:i_ng,ka,ia,ja), intent(in)  rhoq,
real(rp), dimension(hydro_max,ka,ia,ja), intent(in)  xq,
real(rp), dimension(hydro_max,ka,ia,ja), intent(in)  dq_xave,
real(rp), dimension(hydro_max,2,ka,ia,ja), intent(in)  vt_xave 
)

Definition at line 2737 of file scale_atmos_phy_mp_sn14.F90.

References scale_atmos_hydrometeor::i_nc, scale_atmos_hydrometeor::i_ni, scale_atmos_hydrometeor::i_nr, scale_atmos_hydrometeor::i_ns, scale_atmos_hydrometeor::i_qc, scale_atmos_hydrometeor::i_qg, scale_atmos_hydrometeor::i_qi, scale_atmos_hydrometeor::i_qr, scale_atmos_hydrometeor::i_qs, scale_grid_index::ie, scale_stdio::io_fid_conf, scale_stdio::io_fid_nml, scale_stdio::io_nml, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, and scale_tracer::qa.

Referenced by atmos_phy_mp_sn14().

2737 ! rho ) ! [Add] 11/08/30
2738  use scale_tracer, only: &
2739  qa
2740  use scale_atmos_saturation, only: &
2741  moist_psat_ice => atmos_saturation_psat_ice
2742  implicit none
2743 
2744  !--- mixed-phase collection process
2745  ! And all we set all production term as a negative sign to avoid confusion.
2746  !
2747  real(RP), intent(out):: Pac(Pac_MAX,KA,IA,JA)
2748  !--- partial conversion
2749  real(RP), intent(inout):: PQ(PQ_MAX,KA,IA,JA)
2750  !
2751  real(RP), intent(in) :: wtem(KA,IA,JA)
2752  !--- mass/number concentration[kg/m3]
2753  real(RP), intent(in) :: rhoq(I_QV:I_NG,KA,IA,JA)
2754  ! necessary ?
2755  real(RP), intent(in) :: xq(HYDRO_MAX,KA,IA,JA)
2756  !--- diameter of averaged mass( D(ave x) )
2757  real(RP), intent(in) :: dq_xave(HYDRO_MAX,KA,IA,JA)
2758  !--- terminal velocity of averaged mass( vt(ave x) )
2759  real(RP), intent(in) :: vt_xave(HYDRO_MAX,2,KA,IA,JA)
2760  ! [Add] 11/08/30 T.Mitsui, for autoconversion of ice
2761 ! real(RP), intent(in) :: rho(KA,IA,JA)
2762  !
2763  ! namelist variables
2764  !=== for collection
2765  !--- threshold of diameters to collide with others
2766  real(RP), save :: dc0 = 15.0e-6_rp ! lower threshold of cloud
2767  real(RP), save :: dc1 = 40.0e-6_rp ! upper threshold of cloud
2768  real(RP), save :: di0 = 150.0e-6_rp ! lower threshold of cloud
2769  real(RP), save :: ds0 = 150.0e-6_rp ! lower threshold of cloud
2770  real(RP), save :: dg0 = 150.0e-6_rp ! lower threshold of cloud
2771  !--- standard deviation of terminal velocity[m/s]
2772  real(RP), save :: sigma_c=0.00_rp ! cloud
2773  real(RP), save :: sigma_r=0.00_rp ! rain
2774  real(RP), save :: sigma_i=0.2_rp ! ice
2775  real(RP), save :: sigma_s=0.2_rp ! snow
2776  real(RP), save :: sigma_g=0.00_rp ! graupel
2777  !--- max collection efficiency for cloud
2778  real(RP), save :: E_im = 0.80_rp ! ice max
2779  real(RP), save :: E_sm = 0.80_rp ! snow max
2780  real(RP), save :: E_gm = 1.00_rp ! graupel max
2781  !--- collection efficiency between 2 species
2782  real(RP), save :: E_ir=1.0_rp ! ice x rain
2783  real(RP), save :: E_sr=1.0_rp ! snow x rain
2784  real(RP), save :: E_gr=1.0_rp ! graupel x rain
2785  real(RP), save :: E_ii=1.0_rp ! ice x ice
2786  real(RP), save :: E_si=1.0_rp ! snow x ice
2787  real(RP), save :: E_gi=1.0_rp ! graupel x ice
2788  real(RP), save :: E_ss=1.0_rp ! snow x snow
2789  real(RP), save :: E_gs=1.0_rp ! graupel x snow
2790  real(RP), save :: E_gg=1.0_rp ! graupel x graupel
2791  !=== for partial conversion
2792  !--- flag: 1=> partial conversion to graupel, 0=> no conversion
2793  integer, save :: i_iconv2g=1 ! ice => graupel
2794  integer, save :: i_sconv2g=1 ! snow => graupel
2795  !--- bulk density of graupel
2796  real(RP), save :: rho_g = 900.0_rp ! [kg/m3]
2797  !--- space filling coefficient [%]
2798  real(RP), save :: cfill_i = 0.68_rp ! ice
2799  real(RP), save :: cfill_s = 0.01_rp ! snow
2800  !--- critical diameter for ice conversion
2801  real(RP), save :: di_cri = 500.e-6_rp ! [m]
2802  ! [Add] 10/08/03 T.Mitsui
2803  logical, save :: opt_stick_KS96=.false.
2804  logical, save :: opt_stick_CO86=.false.
2805  real(RP), parameter :: a_dec = 0.883_rp
2806  real(RP), parameter :: b_dec = 0.093_rp
2807  real(RP), parameter :: c_dec = 0.00348_rp
2808  real(RP), parameter :: d_dec = 4.5185e-5_rp
2809  !
2810  logical, save :: flag_first = .true.
2811  namelist /nm_mp_sn14_collection/ &
2812  dc0, dc1, di0, ds0, dg0, &
2813  sigma_c, sigma_r, sigma_i, sigma_s, sigma_g, &
2814  opt_stick_ks96, &
2815  opt_stick_co86, &
2816  e_im, e_sm, e_gm, &
2817  e_ir, e_sr, e_gr, e_ii, e_si, e_gi, e_ss, e_gs, e_gg, &
2818  i_iconv2g, i_sconv2g, rho_g, cfill_i, cfill_s, di_cri
2819  !
2820  real(RP) :: tem(KA,IA,JA)
2821  !
2822  !--- collection efficency of each specie
2823  real(RP) :: E_c, E_r, E_i, E_s, E_g !
2824  real(RP) :: E_ic, E_sc, E_gc !
2825  !--- sticking efficiency
2826  real(RP) :: E_stick(KA,IA,JA)
2827  ! [Add] 10/08/03 T.Mitsui
2828  real(RP) :: temc, temc2, temc3
2829  real(RP) :: E_dec
2830  real(RP) :: esi_rat
2831  real(RP) :: esi(KA,IA,JA)
2832  !
2833  real(RP) :: temc_p, temc_m ! celcius tem.
2834  ! [Add] 11/08/30 T.Mitsui, estimation of autoconversion time
2835 ! real(RP) :: ci_aut(KA,IA,JA)
2836 ! real(RP) :: taui_aut(KA,IA,JA)
2837 ! real(RP) :: tau_sce(KA,IA,JA)
2838  !--- DSD averaged diameter for each species
2839  real(RP) :: ave_dc ! cloud
2840 ! real(RP) :: ave_dr ! rain
2841  real(RP) :: ave_di ! ice
2842  real(RP) :: ave_ds ! snow
2843  real(RP) :: ave_dg ! graupel
2844  !--- coefficient of collection equations(L:mass, N:number)
2845  real(RP) :: coef_acc_LCI, coef_acc_NCI ! cloud - cloud ice
2846  real(RP) :: coef_acc_LCS, coef_acc_NCS ! cloud - snow
2847  !
2848  real(RP) :: coef_acc_LCG, coef_acc_NCG ! cloud - graupel
2849  real(RP) :: coef_acc_LRI_I, coef_acc_NRI_I ! rain - cloud ice
2850  real(RP) :: coef_acc_LRI_R, coef_acc_NRI_R ! rain - cloud ice
2851  real(RP) :: coef_acc_LRS_S, coef_acc_NRS_S ! rain - snow
2852  real(RP) :: coef_acc_LRS_R, coef_acc_NRS_R ! rain - snow
2853  real(RP) :: coef_acc_LRG, coef_acc_NRG ! rain - graupel
2854  real(RP) :: coef_acc_LII, coef_acc_NII ! cloud ice - cloud ice
2855  real(RP) :: coef_acc_LIS, coef_acc_NIS ! cloud ice - snow
2856  real(RP) :: coef_acc_NSS ! snow - snow
2857  real(RP) :: coef_acc_NGG ! grauepl - graupel
2858  real(RP) :: coef_acc_LSG, coef_acc_NSG ! snow - graupel
2859  !--- (diameter) x (diameter)
2860  real(RP) :: dcdc, dcdi, dcds, dcdg
2861  real(RP) :: drdr, drdi, drds, drdg
2862  real(RP) :: didi, dids, didg
2863  real(RP) :: dsds, dsdg
2864  real(RP) :: dgdg
2865  !--- (terminal velocity) x (terminal velocity)
2866  real(RP) :: vcvc, vcvi, vcvs, vcvg
2867  real(RP) :: vrvr, vrvi, vrvs, vrvg
2868  real(RP) :: vivi, vivs, vivg
2869  real(RP) :: vsvs, vsvg
2870  real(RP) :: vgvg
2871  !
2872  real(RP) :: wx_cri, wx_crs
2873  real(RP) :: coef_emelt
2874  real(RP) :: w1
2875 
2876  real(RP) :: sw
2877  !
2878  integer :: i, j, k
2879  !
2880  if( flag_first )then
2881  rewind( io_fid_conf )
2882  read( io_fid_conf, nml=nm_mp_sn14_collection, end=100 )
2883 100 if( io_nml ) write(io_fid_nml,nml=nm_mp_sn14_collection)
2884  flag_first = .false.
2885  end if
2886  !
2887  ! [Add] 10/08/03 T.Mitsui
2888  do j = js, je
2889  do i = is, ie
2890  do k = ks, ke
2891  tem(k,i,j) = max( wtem(k,i,j), tem_min ) ! 11/08/30 T.Mitsui
2892  end do
2893  end do
2894  end do
2895 
2896  call moist_psat_ice( esi, tem )
2897 
2898  if( opt_stick_ks96 )then
2899  do j = js, je
2900  do i = is, ie
2901  do k = ks, ke
2902  ! Khain and Sednev (1996), eq.(3.15)
2903  temc = tem(k,i,j) - t00
2904  temc2 = temc*temc
2905  temc3 = temc2*temc
2906  e_dec = max(0.0_rp, a_dec + b_dec*temc + c_dec*temc2 + d_dec*temc3 )
2907  esi_rat = rhoq(i_qv,k,i,j)*rvap*tem(k,i,j)/esi(k,i,j)
2908  e_stick(k,i,j) = min(1.0_rp, e_dec*esi_rat)
2909  end do
2910  end do
2911  end do
2912  else if( opt_stick_co86 )then
2913  do j = js, je
2914  do i = is, ie
2915  do k = ks, ke
2916  ! [Add] 11/08/30 T.Mitsui, Cotton et al. (1986)
2917  temc = min(tem(k,i,j) - t00,0.0_rp)
2918  w1 = 0.035_rp*temc-0.7_rp
2919  e_stick(k,i,j) = 10._rp**w1
2920  end do
2921  end do
2922  end do
2923  else
2924  do j = js, je
2925  do i = is, ie
2926  do k = ks, ke
2927  ! Lin et al. (1983)
2928  temc_m = min(tem(k,i,j) - t00,0.0_rp) ! T < 273.15
2929  e_stick(k,i,j) = exp(0.09_rp*temc_m)
2930  end do
2931  end do
2932  end do
2933  end if
2934  !
2935  profile_start("sn14_collection")
2936 !OCL NOSIMD
2937  do j = js, je
2938  do i = is, ie
2939  do k = ks, ke
2940  !
2941 ! temc_m = min(tem(k,i,j) - T00,0.0_RP) ! T < 273.15
2942  temc_p = max(tem(k,i,j) - t00,0.0_rp) ! T > 273.15
2943  ! averaged diameter using SB06(82)
2944  ave_dc = coef_d(i_mp_qc)*xq(i_mp_qc,k,i,j)**b_m(i_mp_qc)
2945  ave_di = coef_d(i_mp_qi)*xq(i_mp_qi,k,i,j)**b_m(i_mp_qi)
2946  ave_ds = coef_d(i_mp_qs)*xq(i_mp_qs,k,i,j)**b_m(i_mp_qs)
2947  ave_dg = coef_d(i_mp_qg)*xq(i_mp_qg,k,i,j)**b_m(i_mp_qg)
2948  !------------------------------------------------------------------------
2949  ! coellection efficiency are given as follows
2950  e_c = max(0.0_rp, min(1.0_rp, (ave_dc-dc0)/(dc1-dc0) ))
2951  sw = 0.5_rp - sign(0.5_rp, di0-ave_di) ! if(ave_di>di0)then sw=1
2952  e_i = e_im * sw
2953  sw = 0.5_rp - sign(0.5_rp, ds0-ave_ds) ! if(ave_ds>ds0)then sw=1
2954  e_s = e_sm * sw
2955  sw = 0.5_rp - sign(0.5_rp, dg0-ave_dg) ! if(ave_dg>dg0)then sw=1
2956  e_g = e_gm * sw
2957  e_ic = e_i*e_c
2958  e_sc = e_s*e_c
2959  e_gc = e_g*e_c
2960  !------------------------------------------------------------------------
2961  ! Collection: a collects b ( assuming particle size a>b )
2962  dcdc = dq_xave(i_mp_qc,k,i,j) * dq_xave(i_mp_qc,k,i,j)
2963  drdr = dq_xave(i_mp_qr,k,i,j) * dq_xave(i_mp_qr,k,i,j)
2964  didi = dq_xave(i_mp_qi,k,i,j) * dq_xave(i_mp_qi,k,i,j)
2965  dsds = dq_xave(i_mp_qs,k,i,j) * dq_xave(i_mp_qs,k,i,j)
2966  dgdg = dq_xave(i_mp_qg,k,i,j) * dq_xave(i_mp_qg,k,i,j)
2967  dcdi = dq_xave(i_mp_qc,k,i,j) * dq_xave(i_mp_qi,k,i,j)
2968  dcds = dq_xave(i_mp_qc,k,i,j) * dq_xave(i_mp_qs,k,i,j)
2969  dcdg = dq_xave(i_mp_qc,k,i,j) * dq_xave(i_mp_qg,k,i,j)
2970  drdi = dq_xave(i_mp_qr,k,i,j) * dq_xave(i_mp_qi,k,i,j)
2971  drds = dq_xave(i_mp_qr,k,i,j) * dq_xave(i_mp_qs,k,i,j)
2972  drdg = dq_xave(i_mp_qr,k,i,j) * dq_xave(i_mp_qg,k,i,j)
2973  dids = dq_xave(i_mp_qi,k,i,j) * dq_xave(i_mp_qs,k,i,j)
2974 ! didg = dq_xave(I_mp_QI,k,i,j) * dq_xave(I_mp_QG,k,i,j)
2975  dsdg = dq_xave(i_mp_qs,k,i,j) * dq_xave(i_mp_qg,k,i,j)
2976  !
2977  vcvc = vt_xave(i_mp_qc,2,k,i,j)* vt_xave(i_mp_qc,2,k,i,j)
2978  vrvr = vt_xave(i_mp_qr,2,k,i,j)* vt_xave(i_mp_qr,2,k,i,j)
2979  vivi = vt_xave(i_mp_qi,2,k,i,j)* vt_xave(i_mp_qi,2,k,i,j)
2980  vsvs = vt_xave(i_mp_qs,2,k,i,j)* vt_xave(i_mp_qs,2,k,i,j)
2981  vgvg = vt_xave(i_mp_qg,2,k,i,j)* vt_xave(i_mp_qg,2,k,i,j)
2982  vcvi = vt_xave(i_mp_qc,2,k,i,j)* vt_xave(i_mp_qi,2,k,i,j)
2983  vcvs = vt_xave(i_mp_qc,2,k,i,j)* vt_xave(i_mp_qs,2,k,i,j)
2984  vcvg = vt_xave(i_mp_qc,2,k,i,j)* vt_xave(i_mp_qg,2,k,i,j)
2985  vrvi = vt_xave(i_mp_qr,2,k,i,j)* vt_xave(i_mp_qi,2,k,i,j)
2986  vrvs = vt_xave(i_mp_qr,2,k,i,j)* vt_xave(i_mp_qs,2,k,i,j)
2987  vrvg = vt_xave(i_mp_qr,2,k,i,j)* vt_xave(i_mp_qg,2,k,i,j)
2988  vivs = vt_xave(i_mp_qi,2,k,i,j)* vt_xave(i_mp_qs,2,k,i,j)
2989 ! vivg = vt_xave(I_mp_QI,2,k,i,j)* vt_xave(I_mp_QG,2,k,i,j)
2990  vsvg = vt_xave(i_mp_qs,2,k,i,j)* vt_xave(i_mp_qg,2,k,i,j)
2991  !------------------------------------------------------------------------
2992  !
2993  !+++ pattern 1: a + b => a (a>b)
2994  ! (i-c, s-c, g-c, s-i, g-r, s-g)
2995  !------------------------------------------------------------------------
2996  ! cloud-ice => ice
2997  ! reduction term of cloud
2998  coef_acc_lci = &
2999  ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi ) &
3000  * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi &
3001  + sigma_i + sigma_c )**0.5_rp
3002  coef_acc_nci = &
3003  ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi ) &
3004  * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi &
3005  + sigma_i + sigma_c )**0.5_rp
3006  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
3007  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
3008  ! cloud-snow => snow
3009  ! reduction term of cloud
3010  coef_acc_lcs = &
3011  ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds ) &
3012  * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs &
3013  + sigma_s + sigma_c )**0.5_rp
3014  coef_acc_ncs = &
3015  ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds ) &
3016  * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs &
3017  + sigma_s + sigma_c )**0.5_rp
3018  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
3019  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
3020  ! cloud-graupel => graupel
3021  ! reduction term of cloud
3022  coef_acc_lcg = &
3023  ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg ) &
3024  * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg &
3025  + sigma_g + sigma_c )**0.5_rp
3026  coef_acc_ncg = &
3027  ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg ) &
3028  * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg &
3029  + sigma_g + sigma_c )**0.5_rp
3030  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
3031  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
3032  ! snow-graupel => graupel
3033  coef_acc_lsg = &
3034  ( delta_b1(i_mp_qs)*dsds + delta_ab1(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg ) &
3035  * ( theta_b1(i_mp_qs)*vsvs - theta_ab1(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg &
3036  + sigma_g + sigma_s )**0.5_rp
3037  coef_acc_nsg = &
3038  ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg ) &
3039  ! [fix] T.Mitsui 08/05/08
3040  * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg &
3041  + sigma_g + sigma_s )**0.5_rp
3042  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
3043  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
3044  !------------------------------------------------------------------------
3045  ! ice-snow => snow
3046  ! reduction term of ice
3047  coef_acc_lis = &
3048  ( delta_b1(i_mp_qi)*didi + delta_ab1(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds ) &
3049  * ( theta_b1(i_mp_qi)*vivi - theta_ab1(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs &
3050  + sigma_i + sigma_s )**0.5_rp
3051  coef_acc_nis = &
3052  ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds ) &
3053  * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs &
3054  + sigma_i + sigma_s )**0.5_rp
3055  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
3056  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
3057  !
3058  sw = sign(0.5_rp, t00-tem(k,i,j)) + 0.5_rp
3059  ! if ( tem(k,i,j) <= T00 )then
3060  ! rain-graupel => graupel
3061  ! reduction term of rain
3062  ! sw = 1
3063  ! else
3064  ! rain-graupel => rain
3065  ! reduction term of graupel
3066  ! sw = 0
3067  coef_acc_lrg = &
3068  ( ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg ) * sw &
3069  + ( delta_b1(i_mp_qg)*dgdg + delta_ab1(i_mp_qr,i_mp_qg)*drdg + delta_b0(i_mp_qr)*drdr ) * (1.0_rp-sw) ) &
3070  * sqrt( ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qg,i_mp_qr)*vrvg + theta_b0(i_mp_qg)*vgvg ) * sw &
3071  + ( theta_b1(i_mp_qg)*vgvg - theta_ab1(i_mp_qr,i_mp_qg)*vrvg + theta_b0(i_mp_qr)*vrvr ) * (1.0_rp-sw) &
3072  + sigma_r + sigma_g )
3073  pac(i_lraclg2lg,k,i,j) = -0.25_rp*pi*e_gr*coef_acc_lrg &
3074  * ( rhoq(i_ng,k,i,j)*rhoq(i_qr,k,i,j) * sw &
3075  + rhoq(i_nr,k,i,j)*rhoq(i_qg,k,i,j) * (1.0_rp-sw) )
3076  coef_acc_nrg = &
3077  ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg ) &
3078  * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qg,i_mp_qr)*vrvg + theta_b0(i_mp_qg)*vgvg &
3079  + sigma_r + sigma_g )**0.5_rp
3080  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
3081  !
3082  !------------------------------------------------------------------------
3083  !
3084  !+++ pattern 2: a + b => c (a>b)
3085  ! (r-i,r-s)
3086  !------------------------------------------------------------------------
3087  ! rain-ice => graupel
3088  ! reduction term of ice
3089  coef_acc_lri_i = &
3090  ( delta_b1(i_mp_qi)*didi + delta_ab1(i_mp_qr,i_mp_qi)*drdi + delta_b0(i_mp_qr)*drdr ) &
3091  * ( theta_b1(i_mp_qi)*vivi - theta_ab1(i_mp_qr,i_mp_qi)*vrvi + theta_b0(i_mp_qr)*vrvr &
3092  + sigma_r + sigma_i )**0.5_rp
3093  coef_acc_nri_i = &
3094  ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qr,i_mp_qi)*drdi + delta_b0(i_mp_qr)*drdr ) &
3095  * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qr,i_mp_qi)*vrvi + theta_b0(i_mp_qr)*vrvr &
3096  + sigma_r + sigma_i )**0.5_rp
3097  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
3098  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
3099  ! reduction term of rain
3100  coef_acc_lri_r = &
3101  ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qi,i_mp_qr)*drdi + delta_b0(i_mp_qi)*didi ) &
3102  * ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qi,i_mp_qr)*vrvi + theta_b0(i_mp_qi)*vivi &
3103  + sigma_r + sigma_i )**0.5_rp
3104  coef_acc_nri_r = &
3105  ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qi,i_mp_qr)*drdi + delta_b0(i_mp_qi)*didi ) &
3106  * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qi,i_mp_qr)*vrvi + theta_b0(i_mp_qi)*vivi &
3107  + sigma_r + sigma_i )**0.5_rp
3108  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
3109  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
3110  ! rain-snow => graupel
3111  ! reduction term of snow
3112  coef_acc_lrs_s = &
3113  ( delta_b1(i_mp_qs)*dsds + delta_ab1(i_mp_qr,i_mp_qs)*drds + delta_b0(i_mp_qr)*drdr ) &
3114  * ( theta_b1(i_mp_qs)*vsvs - theta_ab1(i_mp_qr,i_mp_qs)*vrvs + theta_b0(i_mp_qr)*vrvr &
3115  + sigma_r + sigma_s )**0.5_rp
3116  coef_acc_nrs_s = &
3117  ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qr,i_mp_qs)*drds + delta_b0(i_mp_qr)*drdr ) &
3118  * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qr,i_mp_qs)*vrvs + theta_b0(i_mp_qr)*vrvr &
3119  + sigma_r + sigma_s )**0.5_rp
3120  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
3121  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
3122  ! reduction term of rain
3123  coef_acc_lrs_r = &
3124  ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qs,i_mp_qr)*drds + delta_b0(i_mp_qs)*dsds ) &
3125  * ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qs,i_mp_qr)*vrvs + theta_b0(i_mp_qs)*vsvs &
3126  + sigma_r + sigma_s )**0.5_rp
3127  coef_acc_nrs_r = &
3128  ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qs,i_mp_qr)*drds + delta_b0(i_mp_qs)*dsds ) &
3129  * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qs,i_mp_qr)*vrvs + theta_b0(i_mp_qs)*vsvs &
3130  + sigma_r + sigma_s )**0.5_rp
3131  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
3132  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
3133  !------------------------------------------------------------------------
3134  !
3135  !+++ pattern 3: a + a => b (i-i)
3136  !
3137  !------------------------------------------------------------------------
3138  ! ice-ice ( reduction is double, but includes double-count)
3139  coef_acc_lii = &
3140  ( delta_b0(i_mp_qi)*didi + delta_ab1(i_mp_qi,i_mp_qi)*didi + delta_b1(i_mp_qi)*didi ) &
3141  * ( theta_b0(i_mp_qi)*vivi - theta_ab1(i_mp_qi,i_mp_qi)*vivi + theta_b1(i_mp_qi)*vivi &
3142  + sigma_i + sigma_i )**0.5_rp
3143  coef_acc_nii = &
3144  ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qi,i_mp_qi)*didi + delta_b0(i_mp_qi)*didi ) &
3145  * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qi,i_mp_qi)*vivi + theta_b0(i_mp_qi)*vivi &
3146  + sigma_i + sigma_i )**0.5_rp
3147  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
3148  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
3149  !
3150 ! ci_aut(k,i,j) = 0.25_RP*pi*E_ii*rhoq(I_NI,k,i,j)*coef_acc_LII
3151 ! taui_aut(k,i,j) = 1._RP/max(E_stick(k,i,j)*ci_aut(k,i,j),1.E-10_RP)
3152 ! 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)
3153  !------------------------------------------------------------------------
3154  !
3155  !+++ pattern 4: a + a => a (s-s)
3156  !
3157  !------------------------------------------------------------------------
3158  ! snow-snow => snow
3159  coef_acc_nss = &
3160  ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qs,i_mp_qs)*dsds + delta_b0(i_mp_qs)*dsds ) &
3161  * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qs,i_mp_qs)*vsvs + theta_b0(i_mp_qs)*vsvs &
3162  + sigma_s + sigma_s )**0.5_rp
3163  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
3164  !
3165  ! graupel-grauple => graupel
3166  coef_acc_ngg = &
3167  ( delta_b0(i_mp_qg)*dgdg + delta_ab0(i_mp_qg,i_mp_qg)*dgdg + delta_b0(i_mp_qg)*dgdg ) &
3168  * ( theta_b0(i_mp_qg)*vgvg - theta_ab0(i_mp_qg,i_mp_qg)*vgvg + theta_b0(i_mp_qg)*vgvg &
3169  + sigma_g + sigma_g )**0.5_rp
3170  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
3171  !
3172  !------------------------------------------------------------------------
3173  !--- Partial conversion
3174  ! SB06(70),(71)
3175  ! i_iconv2g: option whether partial conversions work or not
3176  ! ice-cloud => graupel
3177  sw = 0.5_rp - sign(0.5_rp,di_cri-ave_di) ! if( ave_di > di_cri )then sw=1
3178  wx_cri = cfill_i*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_di*ave_di*ave_di/xq(i_mp_qi,k,i,j) - 1.0_rp ) * sw
3179  pq(i_licon,k,i,j) = i_iconv2g * pac(i_liaclc2li,k,i,j)/max(1.0_rp, wx_cri) * sw
3180  pq(i_nicon,k,i,j) = i_iconv2g * pq(i_licon,k,i,j)/xq(i_mp_qi,k,i,j) * sw
3181 
3182  ! snow-cloud => graupel
3183  wx_crs = cfill_s*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_ds*ave_ds*ave_ds/xq(i_mp_qs,k,i,j) - 1.0_rp )
3184  pq(i_lscon,k,i,j) = i_sconv2g * (pac(i_lsaclc2ls,k,i,j))/max(1.0_rp, wx_crs)
3185  pq(i_nscon,k,i,j) = i_sconv2g * pq(i_lscon,k,i,j)/xq(i_mp_qs,k,i,j)
3186  !------------------------------------------------------------------------
3187  !--- enhanced melting( due to collection-freezing of water droplets )
3188  ! originally from Rutledge and Hobbs(1984). eq.(A.21)
3189  ! if T > 273.15 then temc_p=T-273.15, else temc_p=0
3190  ! 08/05/08 [fix] T.Mitsui LHF00 => LHF0
3191  ! melting occurs around T=273K, so LHF0 is suitable both SIMPLE and EXACT,
3192  ! otherwise LHF can have sign both negative(EXACT) and positive(SIMPLE).
3193 !!$ coef_emelt = -CL/LHF00*temc_p
3194  coef_emelt = cl/lhf0*temc_p
3195  ! cloud-graupel
3196  pq(i_lgacm,k,i,j) = coef_emelt*pac(i_lgaclc2lg,k,i,j)
3197  pq(i_ngacm,k,i,j) = pq(i_lgacm,k,i,j)/xq(i_mp_qg,k,i,j)
3198  ! rain-graupel
3199  pq(i_lgarm,k,i,j) = coef_emelt*pac(i_lraclg2lg,k,i,j)
3200  pq(i_ngarm,k,i,j) = pq(i_lgarm,k,i,j)/xq(i_mp_qg,k,i,j)
3201  ! cloud-snow
3202  pq(i_lsacm,k,i,j) = coef_emelt*(pac(i_lsaclc2ls,k,i,j))
3203  pq(i_nsacm,k,i,j) = pq(i_lsacm,k,i,j)/xq(i_mp_qs,k,i,j)
3204  ! rain-snow
3205  pq(i_lsarm,k,i,j) = coef_emelt*(pac(i_lracls2lg_r,k,i,j)+pac(i_lracls2lg_s,k,i,j))
3206  pq(i_nsarm,k,i,j) = pq(i_lsarm,k,i,j)/xq(i_mp_qg,k,i,j) ! collect? might be I_mp_QS
3207  ! cloud-ice
3208  pq(i_liacm,k,i,j) = coef_emelt*pac(i_liaclc2li,k,i,j)
3209  pq(i_niacm,k,i,j) = pq(i_liacm,k,i,j)/xq(i_mp_qi,k,i,j)
3210  ! rain-ice
3211  pq(i_liarm,k,i,j) = coef_emelt*(pac(i_lracli2lg_r,k,i,j)+pac(i_lracli2lg_i,k,i,j))
3212  pq(i_niarm,k,i,j) = pq(i_liarm,k,i,j)/xq(i_mp_qg,k,i,j) ! collect? might be I_mp_QI
3213  end do
3214  end do
3215  end do
3216  profile_stop("sn14_collection")
3217 
3218  !
3219  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_atmos_saturation
module ATMOSPHERE / Saturation adjustment
Definition: scale_atmos_sub_saturation.F90:16
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
Here is the caller graph for this function:

◆ aut_acc_slc_brk_kij()

subroutine scale_atmos_phy_mp_sn14::aut_acc_slc_brk_kij ( real(rp), dimension(pq_max,ka,ia,ja), intent(inout)  PQ,
real(rp), dimension(i_qv:i_ng,ka,ia,ja), intent(in)  rhoq,
real(rp), dimension(hydro_max,ka,ia,ja), intent(in)  xq,
real(rp), dimension(hydro_max,ka,ia,ja), intent(in)  dq_xave,
real(rp), dimension(ka,ia,ja), intent(in)  rho 
)

Definition at line 3227 of file scale_atmos_phy_mp_sn14.F90.

References scale_atmos_hydrometeor::i_nc, scale_atmos_hydrometeor::i_nr, scale_atmos_hydrometeor::i_qc, scale_atmos_hydrometeor::i_qr, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, and scale_tracer::qa.

Referenced by atmos_phy_mp_sn14().

3227  use scale_tracer, only: &
3228  qa
3229  implicit none
3230 
3231  real(RP), intent(inout) :: PQ(PQ_MAX,KA,IA,JA)
3232  !
3233  real(RP), intent(in) :: rhoq(I_QV:I_NG,KA,IA,JA)
3234  real(RP), intent(in) :: xq(HYDRO_MAX,KA,IA,JA)
3235  real(RP), intent(in) :: dq_xave(HYDRO_MAX,KA,IA,JA)
3236  real(RP), intent(in) :: rho(KA,IA,JA)
3237  !
3238  ! parameter for autoconversion
3239  real(RP), parameter :: kcc = 4.44e+9_rp ! collision efficiency [m3/kg2/sec]
3240  real(RP), parameter :: tau_min = 1.e-20_rp ! empirical filter by T.Mitsui
3241  real(RP), parameter :: rx_sep = 1.0_rp/x_sep ! 1/x_sep, 10/08/03 [Add] T.Mitsui
3242  !
3243  ! parameter for accretion
3244  real(RP), parameter :: kcr = 5.8_rp ! collision efficiency [m3/kg2/sec]
3245  real(RP), parameter :: thr_acc = 5.e-5_rp ! threshold for universal function original
3246  !
3247  ! parameter for self collection and collison break-up
3248  real(RP), parameter :: krr = 4.33_rp ! k_rr, S08 (35)
3249  real(RP), parameter :: kaprr = 60.7_rp ! kappa_rr, SB06(11)
3250  real(RP), parameter :: kbr = 1000._rp ! k_br, SB06(14)
3251  real(RP), parameter :: kapbr = 2.3e+3_rp ! kappa_br, SB06(15)
3252  real(RP), parameter :: dr_min = 0.35e-3_rp ! minimum diameter, SB06(13)-(15)
3253  !
3254  ! work variables
3255  real(RP) :: coef_nuc0 ! coefficient of number for Auto-conversion
3256  real(RP) :: coef_nuc1 ! mass
3257  real(RP) :: coef_aut0 ! number
3258  real(RP) :: coef_aut1 ! mass
3259  real(RP) :: lwc ! lc+lr
3260  real(RP) :: tau ! conversion ratio: qr/(qc+qr) ranges [0:1]
3261  real(RP) :: rho_fac ! factor of air density
3262  real(RP) :: psi_aut ! Universal function of Auto-conversion
3263  real(RP) :: psi_acc ! Universal function of Accretion
3264  real(RP) :: psi_brk ! Universal function of Breakup
3265  real(RP) :: ddr ! diameter difference from equilibrium
3266  !
3267  integer :: i, j, k
3268  !
3269  coef_nuc0 = (nu(i_mp_qc)+2.0_rp)/(nu(i_mp_qc)+1.0_rp)
3270  coef_nuc1 = (nu(i_mp_qc)+2.0_rp)*(nu(i_mp_qc)+4.0_rp)/(nu(i_mp_qc)+1.0_rp)/(nu(i_mp_qc)+1.0_rp)
3271  coef_aut0 = -kcc*coef_nuc0
3272  coef_aut1 = -kcc/x_sep/20._rp*coef_nuc1
3273  !
3274  do j = js, je
3275  do i = is, ie
3276  do k = ks, ke
3277  lwc = rhoq(i_qr,k,i,j)+rhoq(i_qc,k,i,j)
3278  if( lwc > xc_min )then
3279  tau = max(tau_min, rhoq(i_qr,k,i,j)/lwc)
3280  else
3281  tau = tau_min
3282  end if
3283  rho_fac = sqrt(rho_0/max(rho(k,i,j),rho_min))
3284  !
3285  ! Auto-conversion ( cloud-cloud => rain )
3286  psi_aut = 400._rp*(tau**0.7_rp)*(1.0_rp - (tau**0.7_rp))**3 ! (6) SB06
3287  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
3288  ! lc = lwc*(1-tau), lr = lwc*tau
3289  pq(i_lcaut,k,i,j) = coef_aut1*lwc*lwc*xq(i_mp_qc,k,i,j)*xq(i_mp_qc,k,i,j) & ! (4) SB06
3290  *((1.0_rp-tau)*(1.0_rp-tau) + psi_aut)*rho_fac*rho_fac !
3291  pq(i_nraut,k,i,j) = -rx_sep*pq(i_lcaut,k,i,j) ! (A7) SB01
3292  !
3293  ! Accretion ( cloud-rain => rain )
3294  psi_acc =(tau/(tau+thr_acc))**4 ! (8) SB06
3295  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
3296  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
3297  !
3298  ! Self-collection ( rain-rain => rain )
3299  pq(i_nrslc,k,i,j) = -krr*rhoq(i_nr,k,i,j)*rhoq(i_qr,k,i,j)*rho_fac ! (A.8) SB01
3300  !
3301  ! Collisional breakup of rain
3302  ddr = min(1.e-3_rp, dq_xave(i_mp_qr,k,i,j) - dr_eq )
3303  if (dq_xave(i_mp_qr,k,i,j) < dr_min )then ! negligible
3304  psi_brk = -1.0_rp
3305  pq(i_nrbrk,k,i,j) = 0.0_rp
3306  else if (dq_xave(i_mp_qr,k,i,j) <= dr_eq )then
3307  psi_brk = kbr*ddr + 1.0_rp ! (14) SB06 (+1 is necessary)
3308  pq(i_nrbrk,k,i,j) = - (psi_brk + 1.0_rp)*pq(i_nrslc,k,i,j) ! (13) SB06
3309  else
3310  psi_brk = 2.0_rp*exp(kapbr*ddr) - 1.0_rp ! (15) SB06
3311  pq(i_nrbrk,k,i,j) = - (psi_brk + 1.0_rp)*pq(i_nrslc,k,i,j) ! (13) SB06
3312  end if
3313  !
3314  end do
3315  end do
3316  end do
3317  !
3318  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
Here is the caller graph for this function:

◆ dep_vapor_melt_ice_kij()

subroutine scale_atmos_phy_mp_sn14::dep_vapor_melt_ice_kij ( real(rp), dimension(pq_max,ka,ia,ja), intent(inout)  PQ,
real(rp), dimension(ka,ia,ja), intent(in)  rho,
real(rp), dimension(ka,ia,ja), intent(in)  tem,
real(rp), dimension(ka,ia,ja), intent(in)  pre,
real(rp), dimension(ka,ia,ja), intent(in)  qd,
real(rp), dimension(i_qv:i_ng,ka,ia,ja), intent(in)  rhoq,
real(rp), dimension(ka,ia,ja), intent(in)  esw,
real(rp), dimension(ka,ia,ja), intent(in)  esi,
real(rp), dimension(hydro_max,ka,ia,ja), intent(in)  xq,
real(rp), dimension(hydro_max,2,ka,ia,ja), intent(in)  vt_xave,
real(rp), dimension(hydro_max,ka,ia,ja), intent(in)  dq_xave 
)

Definition at line 3327 of file scale_atmos_phy_mp_sn14.F90.

References scale_const::const_eps, scale_atmos_hydrometeor::i_nc, scale_atmos_hydrometeor::i_ni, scale_atmos_hydrometeor::i_nr, scale_atmos_hydrometeor::i_ns, scale_atmos_hydrometeor::i_qg, scale_atmos_hydrometeor::i_qi, scale_atmos_hydrometeor::i_qs, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, and scale_tracer::qa.

Referenced by atmos_phy_mp_sn14().

3327  use scale_const, only: &
3328  eps => const_eps
3329  use scale_tracer, only: &
3330  qa
3331  implicit none
3332 
3333  ! Diffusion growth or Evaporation, Sublimation
3334  real(RP), intent(inout) :: PQ(PQ_MAX,KA,IA,JA) ! mass change for cloud, [Add] 09/08/18 T.Mitsui
3335 
3336  real(RP), intent(in) :: rho(KA,IA,JA) ! air density
3337  real(RP), intent(in) :: tem(KA,IA,JA) ! air temperature
3338  real(RP), intent(in) :: pre(KA,IA,JA) ! air pressure
3339  real(RP), intent(in) :: qd(KA,IA,JA) ! mixing ratio of dry air
3340  real(RP), intent(in) :: esw(KA,IA,JA) ! saturation vapor pressure(liquid water)
3341  real(RP), intent(in) :: esi(KA,IA,JA) ! saturation vapor pressure(solid water)
3342  real(RP), intent(in) :: rhoq(I_QV:I_NG,KA,IA,JA)
3343  real(RP), intent(in) :: xq(HYDRO_MAX,KA,IA,JA) ! mean mass
3344  ! Notice following values differ from mean terminal velocity or diameter.
3345  ! mean(vt(x)) /= vt(mean(x)) and mean(D(x)) /= D(mean(x))
3346  ! Following ones are vt(mean(x)) and D(mean(x)).
3347  real(RP), intent(in) :: vt_xave(HYDRO_MAX,2,KA,IA,JA) ! terminal velocity of mean cloud 09/08/18 [Add], T.Mitsui
3348  !
3349  real(RP), intent(in) :: dq_xave(HYDRO_MAX,KA,IA,JA) ! diameter
3350  !
3351  real(RP) :: rho_lim ! limited density 09/08/18 T.Mitsui
3352  real(RP) :: temc_lim ! limited temperature[celsius] 09/08/18 T.Mitsui
3353  real(RP) :: pre_lim ! limited density 09/08/18 T.Mitsui
3354  real(RP) :: temc ! temperature[celsius]
3355 ! real(RP) :: pv ! vapor pressure
3356  real(RP) :: qv ! mixing ratio of water vapor [Add] 09/08/18
3357 ! real(RP) :: ssw ! super saturation ratio(liquid water)
3358 ! real(RP) :: ssi ! super saturation ratio(ice water)
3359  real(RP) :: nua, r_nua ! kinematic viscosity of air
3360  real(RP) :: mua ! viscosity of air
3361  real(RP) :: Kalfa ! thermal conductance
3362  real(RP) :: Dw ! diffusivity of water vapor
3363  real(RP) :: Dt ! diffusivity of heat
3364  real(RP) :: Gw, Gi ! diffusion factor by balance between heat and vapor
3365  real(RP) :: Gwr, Gii, Gis, Gig ! for rain, ice, snow and graupel.
3366  real(RP) :: Gm ! melting factor by balance between heat and vapor
3367  real(RP) :: Nsc_r3 !
3368  ! [Mod] 11/08/30 T.Mitsui, considering large and small branches
3369 ! real(RP) :: Nrecs_r2 ! 09/08/18 [Add] T.Mitsui
3370  real(RP) :: Nrers_r2, Nreis_r2 !
3371  real(RP) :: Nress_r2, Nregs_r2 !
3372 ! real(RP) :: Nrecl_r2 ! 09/08/18 [Add] T.Mitsui
3373  real(RP) :: Nrerl_r2, Nreil_r2 !
3374  real(RP) :: Nresl_r2, Nregl_r2 !
3375  real(RP) :: NscNrer_s, NscNrer_l
3376  real(RP) :: NscNrei_s, NscNrei_l
3377  real(RP) :: NscNres_s, NscNres_l
3378  real(RP) :: NscNreg_s, NscNreg_l
3379  real(RP) :: ventLR_s, ventLR_l
3380  real(RP) :: ventNI_s, ventNI_l, ventLI_s, ventLI_l
3381  real(RP) :: ventNS_s, ventNS_l, ventLS_s, ventLS_l
3382  real(RP) :: ventNG_s, ventNG_l, ventLG_s, ventLG_l
3383  !
3384  real(RP) :: wtr, wti, wts, wtg
3385  real(RP), parameter :: r_14=1.0_rp/1.4_rp
3386  real(RP), parameter :: r_15=1.0_rp/1.5_rp
3387  !
3388  real(RP) :: ventLR
3389  real(RP) :: ventNI, ventLI
3390  real(RP) :: ventNS, ventLS
3391  real(RP) :: ventNG, ventLG
3392  !
3393  real(RP), parameter :: Re_max=1.e+3_rp
3394  real(RP), parameter :: Re_min=1.e-4_rp
3395 
3396  real(RP) :: sw
3397  !
3398  integer :: i, j, k
3399  !
3400  ! Notice,T.Mitsui
3401  ! Vapor deposition and melting would not be solved iteratively to reach equilibrium.
3402  ! Because following phenomena are not adjustment but transition.
3403  ! Just time-scales differ among them.
3404  ! If we would treat more appropreately, there would be time-splitting method to solve each ones.
3405 
3406  profile_start("sn14_dep_vapor")
3407  do j = js, je
3408  do i = is, ie
3409  do k = ks, ke
3410  temc = tem(k,i,j) - t00 ! degC
3411  temc_lim= max(temc, -40._rp ) ! [Add] 09/08/18 T.Mitsui, Pruppacher and Klett(1997),(13-3)
3412  rho_lim = max(rho(k,i,j),rho_min) ! [Add] 09/08/18 T.Mitsui
3413  qv = rhoq(i_qv,k,i,j)/rho_lim
3414  pre_lim = rho_lim*(qd(k,i,j)*rdry + qv*rvap)*(temc_lim+t00) ![Add] 09/08/18 T.Mitsui
3415  !--------------------------------------------------------------------
3416  ! Diffusion growth part is described in detail
3417  ! by Pruppacher and Klett (1997) Sec. 13.2(liquid) and 13.3(solid)
3418  !
3419  ! G:factor of thermal diffusion(1st.term) and vapor diffusion(2nd. term)
3420  ! SB06(23),(38), Lin et al(31),(52) or others
3421  ! Dw is introduced by Pruppacher and Klett(1997),(13-3)
3422  dw = 0.211e-4_rp* (((temc_lim+t00)/t00)**1.94_rp) *(p00/pre_lim)
3423  kalfa = ka0 + temc_lim*dka_dt
3424  mua = mua0 + temc_lim*dmua_dt
3425  nua = mua/rho_lim
3426  r_nua = 1.0_rp/nua
3427  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))
3428  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))
3429  ! capacities account for their surface geometries
3430  gwr = 4.0_rp*pi/cap(i_mp_qr)/gw
3431  gii = 4.0_rp*pi/cap(i_mp_qi)/gi
3432  gis = 4.0_rp*pi/cap(i_mp_qs)/gi
3433  gig = 4.0_rp*pi/cap(i_mp_qg)/gi
3434  ! vent: ventilation effect( asymmetry vapor field around particles due to aerodynamic )
3435  ! SB06 (30),(31) and each coefficient is by (88),(89)
3436  nsc_r3 = (nua/dw)**(0.33333333_rp) ! (Schmidt number )^(1/3)
3437  !
3438 ! Nrecs_r2 = sqrt(max(Re_min,min(Re_max,vt_xave(I_mp_QC,1,k,i,j)*dq_xave(I_mp_QC,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) cloud
3439  nrers_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qr,1,k,i,j)*dq_xave(i_mp_qr,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) rain
3440  nreis_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qi,1,k,i,j)*dq_xave(i_mp_qi,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) cloud ice
3441  nress_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qs,1,k,i,j)*dq_xave(i_mp_qs,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) snow
3442  nregs_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qg,1,k,i,j)*dq_xave(i_mp_qg,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) graupel
3443  !
3444 ! Nrecl_r2 = sqrt(max(Re_min,min(Re_max,vt_xave(I_mp_QC,2,k,i,j)*dq_xave(I_mp_QC,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) cloud
3445  nrerl_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qr,2,k,i,j)*dq_xave(i_mp_qr,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) rain
3446  nreil_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qi,2,k,i,j)*dq_xave(i_mp_qi,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) cloud ice
3447  nresl_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qs,2,k,i,j)*dq_xave(i_mp_qs,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) snow
3448  nregl_r2 = sqrt(max(re_min,min(re_max,vt_xave(i_mp_qg,2,k,i,j)*dq_xave(i_mp_qg,k,i,j)*r_nua))) ! (Reynolds number)^(1/2) graupel
3449  nscnrer_s=nsc_r3*nrers_r2 ! small rain
3450  nscnrer_l=nsc_r3*nrerl_r2 ! large rain
3451  !
3452  nscnrei_s=nsc_r3*nreis_r2 ! small ice
3453  nscnrei_l=nsc_r3*nreil_r2 ! large ice
3454  !
3455  nscnres_s=nsc_r3*nress_r2 ! small snow
3456  nscnres_l=nsc_r3*nresl_r2 ! large snow
3457  !
3458  nscnreg_s=nsc_r3*nregs_r2 ! small snow
3459  nscnreg_l=nsc_r3*nregl_r2 ! large snow
3460  !
3461  ventlr_s = ah_vent1(i_mp_qr,1) + bh_vent1(i_mp_qr,1)*nscnrer_s
3462  ventlr_l = ah_vent1(i_mp_qr,2) + bh_vent1(i_mp_qr,2)*nscnrer_l
3463  !
3464  ventni_s = ah_vent0(i_mp_qi,1) + bh_vent0(i_mp_qi,1)*nscnrei_s
3465  ventni_l = ah_vent0(i_mp_qi,2) + bh_vent0(i_mp_qi,2)*nscnrei_l
3466  ventli_s = ah_vent1(i_mp_qi,1) + bh_vent1(i_mp_qi,1)*nscnrei_s
3467  ventli_l = ah_vent1(i_mp_qi,2) + bh_vent1(i_mp_qi,2)*nscnrei_l
3468  !
3469  ventns_s = ah_vent0(i_mp_qs,1) + bh_vent0(i_mp_qs,1)*nscnres_s
3470  ventns_l = ah_vent0(i_mp_qs,2) + bh_vent0(i_mp_qs,2)*nscnres_l
3471  ventls_s = ah_vent1(i_mp_qs,1) + bh_vent1(i_mp_qs,1)*nscnres_s
3472  ventls_l = ah_vent1(i_mp_qs,2) + bh_vent1(i_mp_qs,2)*nscnres_l
3473  !
3474  ventng_s = ah_vent0(i_mp_qg,1) + bh_vent0(i_mp_qg,1)*nscnreg_s
3475  ventng_l = ah_vent0(i_mp_qg,2) + bh_vent0(i_mp_qg,2)*nscnreg_l
3476  ventlg_s = ah_vent1(i_mp_qg,1) + bh_vent1(i_mp_qg,1)*nscnreg_s
3477  ventlg_l = ah_vent1(i_mp_qg,2) + bh_vent1(i_mp_qg,2)*nscnreg_l
3478  !
3479  ! branch is 1.4 for rain, snow, graupel; is 1.0 for ice (PK97, 13-60,-61,-88,-89).
3480  !
3481  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
3482  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
3483  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
3484  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
3485  ! interpolation between two branches
3486  ventni = (1.0_rp-wti)*ventni_s + wti*ventni_l
3487  ventns = (1.0_rp-wts)*ventns_s + wts*ventns_l
3488  ventng = (1.0_rp-wtg)*ventng_s + wtg*ventng_l
3489  !
3490  ventlr = (1.0_rp-wtr)*ventlr_s + wtr*ventlr_l
3491  ventli = (1.0_rp-wti)*ventli_s + wti*ventli_l
3492  ventls = (1.0_rp-wts)*ventls_s + wts*ventls_l
3493  ventlg = (1.0_rp-wtg)*ventlg_s + wtg*ventlg_l
3494  !
3495  ! SB06(29)
3496  ! [Mod] 08/05/08 T.Mitsui, recover PNXdep, and rain is only evaporation.
3497  ! Ni, Ns, Ng should decrease in nature so we add this term.
3498  ! And vapor deposition never occur unless number exist.
3499  ! [Add comment] 09/08/18
3500  ! recover condensation/evaporation of rain,
3501  ! and ventilation effects are not taken into account for cloud.
3502  !
3503 !!$***************************************************************************
3504 !!$ NOTICE:
3505 !!$ 09/08/18 [Mod] Hereafter PLxdep means inverse of timescale.
3506 !!$***************************************************************************
3507 !!$ PQ(I_LCdep,k,i,j) = Gwr*ssw*rhoq(I_NC,k,i,j)*dq_xave(I_mp_QC,k,i,j)*coef_deplc
3508 !!$ PQ(I_LRdep,k,i,j) = Gwr*ssw*rhoq(I_NR,k,i,j)*dq_xave(I_mp_QR,k,i,j)*ventLR
3509 !!$ PQ(I_LIdep,k,i,j) = Gii*ssi*rhoq(I_NI,k,i,j)*dq_xave(I_mp_QI,k,i,j)*ventLI
3510 !!$ PQ(I_LSdep,k,i,j) = Gis*ssi*rhoq(I_NS,k,i,j)*dq_xave(I_mp_QS,k,i,j)*ventLS
3511 !!$ PQ(I_LGdep,k,i,j) = Gig*ssi*rhoq(I_NG,k,i,j)*dq_xave(I_mp_QG,k,i,j)*ventLG
3512  pq(i_lcdep,k,i,j) = gwr*rhoq(i_nc,k,i,j)*dq_xave(i_mp_qc,k,i,j)*coef_deplc
3513  pq(i_lrdep,k,i,j) = gwr*rhoq(i_nr,k,i,j)*dq_xave(i_mp_qr,k,i,j)*ventlr
3514  pq(i_lidep,k,i,j) = gii*rhoq(i_ni,k,i,j)*dq_xave(i_mp_qi,k,i,j)*ventli
3515  pq(i_lsdep,k,i,j) = gis*rhoq(i_ns,k,i,j)*dq_xave(i_mp_qs,k,i,j)*ventls
3516  pq(i_lgdep,k,i,j) = gig*rhoq(i_ng,k,i,j)*dq_xave(i_mp_qg,k,i,j)*ventlg
3517  pq(i_nrdep,k,i,j) = pq(i_lrdep,k,i,j)/xq(i_mp_qr,k,i,j)
3518  pq(i_nidep,k,i,j) = 0.0_rp
3519  pq(i_nsdep,k,i,j) = pq(i_lsdep,k,i,j)/xq(i_mp_qs,k,i,j)
3520  pq(i_ngdep,k,i,j) = pq(i_lgdep,k,i,j)/xq(i_mp_qg,k,i,j)
3521  !
3522  !------------------------------------------------------------------------
3523  ! Melting part is described by Pruppacher and Klett (1997) Sec.16.3.1
3524  ! Here we omit "Shedding" of snow-flakes and ice-particles.
3525  ! "Shedding" may be applicative if you refer
3526  ! eq.(38) in Cotton etal.(1986) Jour. Clim. Appl. Meteor. p.1658-1680.
3527  ! SB06(73)
3528  dt = kalfa/(cpvap*rho_0)
3529  ! Gm: factor caused by balance between
3530  ! "water evaporation cooling(1st.)" and "fusion heating(2nd.)"
3531  ! SB06(76)
3532  ! [fix] 08/05/08 T.Mitsui LHF00 => EMELT and esw => PSAT0
3533  ! LHS0 is more suitable than LHS because melting occurs around 273.15 K.
3534  gm = 2.0_rp*pi/emelt&
3535  * ( (kalfa*dt/dw)*(temc) + (dw*lhs0/rvap)*(esi(k,i,j)/tem(k,i,j)-psat0/t00) )
3536  ! SB06(76)
3537  ! Notice! melting only occurs where T > 273.15 K else doesn't.
3538  ! [fix] 08/05/08 T.Mitsui, Gm could be both positive and negative value.
3539  ! See Pruppacher and Klett(1997) eq.(16-79) or Rasmussen and Pruppacher(1982)
3540  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
3541  ! if Gm==0 then rh and tem is critical value for melting process.
3542  ! 08/05/16 [Mod] T.Mitsui, change term of PLimlt. N_i => L_i/ (limited x_i)
3543  ! because melting never occur when N_i=0.
3544  pq(i_limlt,k,i,j) = - gm * rhoq(i_qi,k,i,j)*dq_xave(i_mp_qi,k,i,j)*ventli/xq(i_mp_qi,k,i,j) * sw
3545  ! [Mod] 08/08/23 T.Mitsui for Seifert(2008)
3546  pq(i_nimlt,k,i,j) = - gm * rhoq(i_ni,k,i,j)*dq_xave(i_mp_qi,k,i,j)*ventni/xq(i_mp_qi,k,i,j) * sw ! 09/08/18 [Mod] recover, T.Mitsui
3547  pq(i_lsmlt,k,i,j) = - gm * rhoq(i_qs,k,i,j)*dq_xave(i_mp_qs,k,i,j)*ventls/xq(i_mp_qs,k,i,j) * sw
3548  ! [Mod] 08/08/23 T.Mitsui for Seifert(2008)
3549  pq(i_nsmlt,k,i,j) = - gm * rhoq(i_ns,k,i,j)*dq_xave(i_mp_qs,k,i,j)*ventns/xq(i_mp_qs,k,i,j) * sw ! 09/08/18 [Mod] recover, T.Mitsui
3550  pq(i_lgmlt,k,i,j) = - gm * rhoq(i_qg,k,i,j)*dq_xave(i_mp_qg,k,i,j)*ventlg/xq(i_mp_qg,k,i,j) * sw
3551  ! [Mod] 08/08/23 T.Mitsui for Seifert(2008)
3552  pq(i_ngmlt,k,i,j) = - gm * rhoq(i_ng,k,i,j)*dq_xave(i_mp_qg,k,i,j)*ventng/xq(i_mp_qg,k,i,j) * sw ! 09/08/18 [Mod] recover, T.Mitsui
3553 
3554  end do
3555  end do
3556  end do
3557  profile_stop("sn14_dep_vapor")
3558  !
3559  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:36
Here is the caller graph for this function:

◆ freezing_water_kij()

subroutine scale_atmos_phy_mp_sn14::freezing_water_kij ( real(rp), intent(in)  dt,
real(rp), dimension(pq_max,ka,ia,ja), intent(inout)  PQ,
real(rp), dimension(i_qv:i_ng,ka,ia,ja), intent(in)  rhoq,
real(rp), dimension(hydro_max,ka,ia,ja), intent(in)  xq,
real(rp), dimension(ka,ia,ja), intent(in)  tem 
)

Definition at line 3566 of file scale_atmos_phy_mp_sn14.F90.

References scale_const::const_undef, scale_atmos_hydrometeor::i_nc, scale_atmos_hydrometeor::i_ni, scale_atmos_hydrometeor::i_nr, scale_atmos_hydrometeor::i_ns, scale_atmos_hydrometeor::i_qc, scale_atmos_hydrometeor::i_qg, scale_atmos_hydrometeor::i_qi, scale_atmos_hydrometeor::i_qr, scale_atmos_hydrometeor::i_qs, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, and scale_tracer::qa.

Referenced by atmos_phy_mp_sn14().

3566  use scale_tracer, only: &
3567  qa
3568  implicit none
3569  !
3570  ! In this subroutine,
3571  ! We assumed surface temperature of droplets are same as environment.
3572 
3573  real(RP), intent(in) :: dt
3574  real(RP), intent(inout):: PQ(PQ_MAX,KA,IA,JA)
3575  !
3576  real(RP), intent(in) :: tem(KA,IA,JA)
3577  !
3578  real(RP), intent(in) :: rhoq(I_QV:I_NG,KA,IA,JA)
3579  real(RP), intent(in) :: xq(HYDRO_MAX,KA,IA,JA)
3580  !
3581  real(RP), parameter :: temc_min = -65.0_rp
3582  real(RP), parameter :: a_het = 0.2_rp ! SB06 (44)
3583  real(RP), parameter :: b_het = 0.65_rp ! SB06 (44)
3584  !
3585  real(RP) :: coef_m2_c
3586  real(RP) :: coef_m2_r
3587  ! temperature [celsius]
3588  real(RP) :: temc, temc2, temc3, temc4
3589  ! temperature function of homegenous/heterogenous freezing
3590  real(RP) :: Jhom, Jhet
3591  real(RP) :: rdt
3592  real(RP) :: tmp
3593  !
3594  integer :: i,j,k
3595  !
3596  rdt = 1.0_rp/dt
3597  !
3598  coef_m2_c = coef_m2(i_mp_qc)
3599  coef_m2_r = coef_m2(i_mp_qr)
3600  !
3601  profile_start("sn14_freezing")
3602  do j = js, je
3603  do i = is, ie
3604  do k = ks, ke
3605  temc = max( tem(k,i,j) - t00, temc_min )
3606  ! These cause from aerosol-droplet interaction.
3607  ! Bigg(1953) formula, Khain etal.(2000) eq.(4.5), Pruppacher and Klett(1997) eq.(9-48)
3608  jhet = a_het*exp( -b_het*temc - 1.0_rp )
3609  ! These cause in nature.
3610  ! Cotton and Field 2002, QJRMS. (12)
3611  if( temc < -65.0_rp )then
3612  jhom = 10.0_rp**(24.37236_rp)*1.e+3_rp
3613  jhet = a_het*exp( 65.0_rp*b_het - 1.0_rp ) ! 09/04/14 [Add], fixer T.Mitsui
3614  else if( temc < -30.0_rp ) then
3615  temc2 = temc*temc
3616  temc3 = temc*temc2
3617  temc4 = temc2*temc2
3618  jhom = 10.0_rp**(&
3619  - 243.40_rp - 14.75_rp*temc - 0.307_rp*temc2 &
3620  - 0.00287_rp*temc3 - 0.0000102_rp*temc4 ) *1.e+3_rp
3621  else if( temc < 0.0_rp) then
3622  jhom = 10._rp**(-7.63_rp-2.996_rp*(temc+30.0_rp))*1.e+3_rp
3623  else
3624  jhom = 0.0_rp
3625  jhet = 0.0_rp
3626  end if
3627  ! Note, xc should be limited in range[xc_min:xc_max].
3628  ! and PNChom need to be calculated by NC
3629  ! because reduction rate of Nc need to be bound by NC.
3630  ! For the same reason PLChom also bound by LC and xc.
3631  ! Basically L and N should be independent
3632  ! but average particle mass x should be in suitable range.
3633  ! Homogenous Freezing
3634  pq(i_lchom,k,i,j) = 0.0_rp
3635  pq(i_nchom,k,i,j) = 0.0_rp
3636  ! Heterogenous Freezing
3637 #if defined(PGI) || defined(SX)
3638  tmp = min( xq(i_mp_qc,k,i,j)*(jhet+jhom)*dt, 1.e+3_rp) ! apply exp limiter
3639  pq(i_lchet,k,i,j) = -rdt*rhoq(i_qc,k,i,j)*( 1.0_rp - exp( -coef_m2_c*tmp ) )
3640  pq(i_nchet,k,i,j) = -rdt*rhoq(i_nc,k,i,j)*( 1.0_rp - exp( - tmp ) )
3641 
3642  tmp = min( xq(i_mp_qr,k,i,j)*(jhet+jhom)*dt, 1.e+3_rp) ! apply exp limiter
3643  pq(i_lrhet,k,i,j) = -rdt*rhoq(i_qr,k,i,j)*( 1.0_rp - exp( -coef_m2_r*tmp ) )
3644  pq(i_nrhet,k,i,j) = -rdt*rhoq(i_nr,k,i,j)*( 1.0_rp - exp( - tmp ) )
3645 #else
3646  pq(i_lchet,k,i,j) = -rdt*rhoq(i_qc,k,i,j)*( 1.0_rp - exp( -coef_m2_c*xq(i_mp_qc,k,i,j)*(jhet+jhom)*dt ) )
3647  pq(i_nchet,k,i,j) = -rdt*rhoq(i_nc,k,i,j)*( 1.0_rp - exp( - xq(i_mp_qc,k,i,j)*(jhet+jhom)*dt ) )
3648  pq(i_lrhet,k,i,j) = -rdt*rhoq(i_qr,k,i,j)*( 1.0_rp - exp( -coef_m2_r*xq(i_mp_qr,k,i,j)*(jhet+jhom)*dt ) )
3649  pq(i_nrhet,k,i,j) = -rdt*rhoq(i_nr,k,i,j)*( 1.0_rp - exp( - xq(i_mp_qr,k,i,j)*(jhet+jhom)*dt ) )
3650 #endif
3651  end do
3652  end do
3653  end do
3654  profile_stop("sn14_freezing")
3655  !
3656  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
Here is the caller graph for this function:

◆ update_by_phase_change_kij()

subroutine scale_atmos_phy_mp_sn14::update_by_phase_change_kij ( integer, intent(in)  ntdiv,
integer, intent(in)  ntmax,
real(rp), intent(in)  dt,
real(rp), dimension(ka,ia,ja), intent(in)  gsgam2,
real(rp), dimension(ka), intent(in)  z,
real(rp), dimension(ka), intent(in)  dz,
real(rp), dimension(ka,ia,ja), intent(in)  velz,
real(rp), dimension(ka,ia,ja), intent(in)  dTdt_rad,
real(rp), dimension(ka,ia,ja), intent(in)  rho,
real(rp), dimension(ka,ia,ja), intent(inout)  rhoe,
real(rp), dimension(i_qv:i_ng,ka,ia,ja), intent(inout)  rhoq,
real(rp), dimension(ka,ia,ja,qa), intent(inout)  q,
real(rp), dimension(ka,ia,ja), intent(inout)  tem,
real(rp), dimension(ka,ia,ja), intent(inout)  pre,
real(rp), dimension(ka,ia,ja), intent(out)  cva,
real(rp), dimension(ka,ia,ja), intent(in)  esw,
real(rp), dimension(ka,ia,ja), intent(in)  esi,
real(rp), dimension(i_qv:i_ng,ka,ia,ja), intent(in)  rhoq2,
real(rp), dimension(pq_max,ka,ia,ja), intent(inout)  PQ,
real(rp), dimension(ka,ia,ja), intent(out)  qc_evaporate,
real(rp), dimension(ia,ja), intent(inout)  sl_PLCdep,
real(rp), dimension(ia,ja), intent(inout)  sl_PLRdep,
real(rp), dimension(ia,ja), intent(inout)  sl_PNRdep 
)

Definition at line 3847 of file scale_atmos_phy_mp_sn14.F90.

References scale_atmos_saturation::atmos_saturation_dqsi_dtem_dpre(), scale_atmos_saturation::atmos_saturation_dqsi_dtem_rho(), scale_atmos_saturation::atmos_saturation_dqsw_dtem_dpre(), scale_atmos_saturation::atmos_saturation_dqsw_dtem_rho(), scale_atmos_hydrometeor::i_nc, scale_atmos_hydrometeor::i_ni, scale_atmos_hydrometeor::i_nr, scale_atmos_hydrometeor::i_ns, scale_atmos_hydrometeor::i_qc, scale_atmos_hydrometeor::i_qg, scale_atmos_hydrometeor::i_qi, scale_atmos_hydrometeor::i_qr, scale_atmos_hydrometeor::i_qs, scale_grid_index::ie, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_fid_nml, scale_stdio::io_l, scale_stdio::io_nml, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, scale_tracer::qa, scale_tracer::tracer_cp, scale_tracer::tracer_cv, scale_tracer::tracer_mass, and scale_tracer::tracer_r.

Referenced by atmos_phy_mp_sn14().

3847  use scale_tracer, only: &
3848  qa, &
3849  tracer_r, &
3850  tracer_cv, &
3851  tracer_cp, &
3852  tracer_mass
3853  use scale_atmos_saturation, only: &
3854  moist_pres2qsat_liq => atmos_saturation_pres2qsat_liq, &
3855  moist_pres2qsat_ice => atmos_saturation_pres2qsat_ice, &
3856  moist_dqsw_dtem_rho => atmos_saturation_dqsw_dtem_rho, &
3857  moist_dqsi_dtem_rho => atmos_saturation_dqsi_dtem_rho, &
3858  moist_dqsw_dtem_dpre => atmos_saturation_dqsw_dtem_dpre, &
3859  moist_dqsi_dtem_dpre => atmos_saturation_dqsi_dtem_dpre
3860  implicit none
3861 
3862  integer, intent(in) :: ntdiv ! [Add] 10/08/03
3863  integer, intent(in) :: ntmax ! [Add] 10/08/03
3864  !
3865  real(RP), intent(in) :: dt ! time step[s]
3866  real(RP), intent(in) :: gsgam2(KA,IA,JA) ! metric
3867  real(RP), intent(in) :: z(KA) ! altitude [m]
3868  real(RP), intent(in) :: dz(KA) ! altitude [m]
3869  real(RP), intent(in) :: velz(KA,IA,JA) ! vertical velocity @ half point[m/s]
3870  real(RP), intent(in) :: dTdt_rad(KA,IA,JA) ! temperture tendency by radiation[K/s]
3871  real(RP), intent(in) :: rho(KA,IA,JA) ! density[kg/m3]
3872  real(RP), intent(inout) :: rhoe(KA,IA,JA) ! internal energy[J/m3]
3873  real(RP), intent(inout) :: rhoq(I_QV:I_NG,KA,IA,JA) ! tracers[kg/m3]
3874  real(RP), intent(inout) :: q(KA,IA,JA,QA) ! tracers mixing ratio[kg/kg]
3875  real(RP), intent(inout) :: tem(KA,IA,JA) ! temperature[K]
3876  real(RP), intent(inout) :: pre(KA,IA,JA) ! pressure[Pa]
3877  real(RP), intent(out) :: cva(KA,IA,JA) ! specific heat at constant volume
3878  real(RP), intent(in) :: esw(KA,IA,JA) ! saturated vapor pressure for liquid
3879  real(RP), intent(in) :: esi(KA,IA,JA) ! for ice
3880  real(RP), intent(in) :: rhoq2(I_QV:I_NG,KA,IA,JA)
3881  !+++ tendency[kg/m3/s]
3882  real(RP), intent(inout) :: PQ(PQ_MAX,KA,IA,JA)
3883  real(RP), intent(out) :: qc_evaporate(KA,IA,JA)
3884  !+++ Column integrated tendency[kg/m2/s]
3885  real(RP), intent(inout) :: sl_PLCdep(IA,JA)
3886  real(RP), intent(inout) :: sl_PLRdep(IA,JA), sl_PNRdep(IA,JA)
3887  !
3888  real(RP) :: Rmoist
3889  !
3890  real(RP) :: xi ! mean mass of ice particles
3891  real(RP) :: rrho ! 1/rho
3892  real(RP) :: wtem(KA,IA,JA) ! temperature[K]
3893  real(RP) :: qdry ! mixing ratio of dry air
3894  !
3895  real(RP) :: r_cva ! specific heat at constant volume
3896  real(RP) :: cpa ! specific heat at constant pressure
3897  real(RP) :: r_cpa ! specific heat at constant pressure
3898  real(RP) :: qsw(KA,IA,JA) ! saturated mixing ratio for liquid
3899  real(RP) :: qsi(KA,IA,JA) ! saturated mixing ratio for solid
3900  real(RP) :: dqswdtem_rho(KA,IA,JA) ! (dqsw/dtem)_rho
3901  real(RP) :: dqsidtem_rho(KA,IA,JA) ! (dqsi/dtem)_rho
3902  real(RP) :: dqswdtem_pre(KA,IA,JA) ! (dqsw/dtem)_pre
3903  real(RP) :: dqsidtem_pre(KA,IA,JA) ! (dqsi/dtem)_pre
3904  real(RP) :: dqswdpre_tem(KA,IA,JA) ! (dqsw/dpre)_tem
3905  real(RP) :: dqsidpre_tem(KA,IA,JA) ! (dqsi/dpre)_tem
3906  !
3907  real(RP) :: w ! vetical velocity[m/s]
3908  real(RP) :: Acnd ! Pdynliq + Bergeron-Findeisen
3909  real(RP) :: Adep ! Pdyndep + Bergeron-Findeisen
3910  real(RP) :: aliqliq, asolliq
3911  real(RP) :: aliqsol, asolsol
3912  real(RP) :: Pdynliq ! production term of ssw by vertical motion
3913  real(RP) :: Pdynsol ! production term of ssi by vertical motion
3914  real(RP) :: Pradliq ! production term of ssw by radiation
3915  real(RP) :: Pradsol ! production term of ssi by radiation
3916  real(RP) :: taucnd, r_taucnd ! time scale of ssw change by MP
3917  real(RP) :: taudep, r_taudep ! time scale of ssi change by MP
3918  real(RP) :: taucnd_c(KA,IA,JA), r_taucnd_c ! by cloud
3919  real(RP) :: taucnd_r(KA,IA,JA), r_taucnd_r ! by rain
3920  real(RP) :: taudep_i(KA,IA,JA), r_taudep_i ! by ice
3921  real(RP) :: taudep_s(KA,IA,JA), r_taudep_s ! by snow
3922  real(RP) :: taudep_g(KA,IA,JA), r_taudep_g ! by graupel
3923  ! alternative tendency through changing ssw and ssi
3924  real(RP) :: PNCdep ! [Add] 11/08/30 T.Mitsui
3925  real(RP) :: PLR2NR, PLI2NI, PLS2NS, PLG2NG
3926  real(RP) :: coef_a_cnd, coef_b_cnd
3927  real(RP) :: coef_a_dep, coef_b_dep
3928  !
3929  real(RP) :: frz_dqc
3930  real(RP) :: frz_dnc
3931  real(RP) :: frz_dqr
3932  real(RP) :: frz_dnr
3933  real(RP) :: mlt_dqi
3934  real(RP) :: mlt_dni
3935  real(RP) :: mlt_dqs
3936  real(RP) :: mlt_dns
3937  real(RP) :: mlt_dqg
3938  real(RP) :: mlt_dng
3939  real(RP) :: dep_dqi
3940  real(RP) :: dep_dni
3941  real(RP) :: dep_dqs
3942  real(RP) :: dep_dns
3943  real(RP) :: dep_dqg
3944  real(RP) :: dep_dng
3945  real(RP) :: dep_dqr
3946  real(RP) :: dep_dnr
3947  real(RP) :: dep_dqc
3948  real(RP) :: dep_dnc ! 11/08/30 [Add] T.Mitsui, dep_dnc
3949  real(RP) :: r_xc_ccn, r_xi_ccn ! 11/08/30 [Add] T.Mitsui
3950  !
3951  real(RP) :: drhoqv
3952  real(RP) :: drhoqc, drhoqr, drhoqi, drhoqs, drhoqg
3953  real(RP) :: drhonc, drhonr, drhoni, drhons, drhong
3954  !
3955  real(RP) :: fac1, fac2, fac3, fac4, fac5, fac6
3956  real(RP) :: r_rvaptem ! 1/(Rvap*tem)
3957  real(RP) :: pv ! vapor pressure
3958  real(RP) :: lvsw, lvsi ! saturated vapor density
3959  real(RP) :: dlvsw, dlvsi !
3960  ! [Add] 11/08/30 T.Mitsui
3961  real(RP) :: dcnd, ddep ! total cndensation/deposition
3962  real(RP) :: uplim_cnd ! upper limit of condensation
3963  real(RP) :: lowlim_cnd ! lower limit of evaporation
3964  ! [Add] 11/08/30 T.Mitsui
3965  real(RP) :: uplim_dep ! upper limit of condensation
3966  real(RP) :: lowlim_dep ! lower limit of evaporation
3967  real(RP) :: ssw, ssi ! supersaturation ratio
3968  real(RP) :: r_esw, r_esi ! 1/esw, 1/esi
3969  real(RP) :: r_lvsw, r_lvsi ! 1/(lvsw*ssw), 1/(lvsi*ssi)
3970  real(RP) :: r_dt ! 1/dt
3971  real(RP) :: ssw_o, ssi_o
3972 ! real(RP) :: dt_dyn
3973 ! real(RP) :: dt_mp
3974  !
3975 ! real(RP) :: tem_lh(KA,IA,JA)
3976 ! real(RP) :: dtemdt_lh(KA,IA,JA)
3977  real(RP), save :: fac_cndc = 1.0_rp
3978  logical, save :: opt_fix_taucnd_c=.false.
3979  logical, save :: flag_first =.true.
3980  !
3981  namelist /nm_mp_sn14_condensation/ &
3982  opt_fix_taucnd_c, fac_cndc
3983 
3984  real(RP) :: fac_cndc_wrk
3985  !
3986  real(RP), parameter :: tau100day = 1.e+7_rp
3987  real(RP), parameter :: r_tau100day = 1.e-7_rp
3988  real(RP), parameter :: eps=1.e-30_rp
3989  !
3990  integer :: i,j,k,iqw
3991  real(RP) :: sw
3992  !
3993 
3994  ! [Add] 11/08/30 T.Mitsui
3995  if( flag_first )then
3996  flag_first = .false.
3997  rewind(io_fid_conf)
3998  read (io_fid_conf,nml=nm_mp_sn14_condensation, end=100)
3999 100 if( io_nml ) write(io_fid_nml,nml=nm_mp_sn14_condensation)
4000  end if
4001  !
4002 ! dt_dyn = dt*ntmax
4003 ! dt_mp = dt*(ntdiv-1)
4004  !
4005  r_dt = 1.0_rp/dt
4006  !
4007  r_xc_ccn=1.0_rp/xc_ccn
4008 ! r_xi_ccn=1.0_RP/xi_ccn
4009  !
4010  if( opt_fix_taucnd_c )then
4011  fac_cndc_wrk = fac_cndc**(1.0_rp-b_m(i_mp_qc))
4012  do j = js, je
4013  do i = is, ie
4014  do k = ks, ke
4015  pq(i_lcdep,k,i,j) = pq(i_lcdep,k,i,j)*fac_cndc_wrk
4016  end do
4017  end do
4018  end do
4019  if( io_l ) write(io_fid_log,*) "taucnd:fac_cndc_wrk=",fac_cndc_wrk
4020  end if
4021 
4022 !OCL XFILL
4023  do j = js, je
4024  do i = is, ie
4025  do k = ks, ke
4026  ! Temperature lower limit is only used for saturation condition.
4027  ! On the other hand original "tem" is used for calculation of latent heat or energy equation.
4028  wtem(k,i,j) = max( tem(k,i,j), tem_min )
4029  end do
4030  end do
4031  end do
4032 
4033  call moist_pres2qsat_liq ( qsw, wtem, pre )
4034  call moist_pres2qsat_ice ( qsi, wtem, pre )
4035  call moist_dqsw_dtem_rho ( dqswdtem_rho, wtem, rho )
4036  call moist_dqsi_dtem_rho ( dqsidtem_rho, wtem, rho )
4037  call moist_dqsw_dtem_dpre( dqswdtem_pre, dqswdpre_tem, wtem, pre )
4038  call moist_dqsi_dtem_dpre( dqsidtem_pre, dqsidpre_tem, wtem, pre )
4039 
4040  profile_start("sn14_update")
4041  do j = js, je
4042  do i = is, ie
4043  do k = ks, ke
4044  if( z(k) <= 25000.0_rp )then
4045  w = 0.5_rp*(velz(k-1,i,j) + velz(k,i,j))
4046  else
4047  w = 0.0_rp
4048  end if
4049  if( pre(k,i,j) < esw(k,i,j)+1.e-10_rp )then
4050  qsw(k,i,j) = 1.0_rp
4051  dqswdtem_rho(k,i,j) = 0.0_rp
4052  dqswdtem_pre(k,i,j) = 0.0_rp
4053  dqswdpre_tem(k,i,j) = 0.0_rp
4054  end if
4055  if( pre(k,i,j) < esi(k,i,j)+1.e-10_rp )then
4056  qsi(k,i,j) = 1.0_rp
4057  dqsidtem_rho(k,i,j) = 0.0_rp
4058  dqsidtem_pre(k,i,j) = 0.0_rp
4059  dqsidpre_tem(k,i,j) = 0.0_rp
4060  end if
4061 
4062  r_rvaptem = 1.0_rp/(rvap*wtem(k,i,j))
4063  lvsw = esw(k,i,j)*r_rvaptem ! rho=p/(Rv*T)
4064  lvsi = esi(k,i,j)*r_rvaptem !
4065  pv = rhoq2(i_qv,k,i,j)*rvap*tem(k,i,j)
4066  r_esw = 1.0_rp/esw(k,i,j)
4067  r_esi = 1.0_rp/esi(k,i,j)
4068  ssw = min( mp_ssw_lim, ( pv*r_esw-1.0_rp ) )
4069  ssi = pv*r_esi - 1.0_rp
4070  r_lvsw = 1.0_rp/lvsw
4071  r_lvsi = 1.0_rp/lvsi
4072  r_taucnd_c = pq(i_lcdep,k,i,j)*r_lvsw
4073  r_taucnd_r = pq(i_lrdep,k,i,j)*r_lvsw
4074  r_taudep_i = pq(i_lidep,k,i,j)*r_lvsi
4075  r_taudep_s = pq(i_lsdep,k,i,j)*r_lvsi
4076  r_taudep_g = pq(i_lgdep,k,i,j)*r_lvsi
4077 ! taucnd_c(k,i,j) = 1.0_RP/(r_taucnd_c+r_tau100day)
4078 ! taucnd_r(k,i,j) = 1.0_RP/(r_taucnd_r+r_tau100day)
4079 ! taudep_i(k,i,j) = 1.0_RP/(r_taudep_i+r_tau100day)
4080 ! taudep_s(k,i,j) = 1.0_RP/(r_taudep_s+r_tau100day)
4081 ! taudep_g(k,i,j) = 1.0_RP/(r_taudep_g+r_tau100day)
4082 
4083  calc_qdry( qdry, q, tracer_mass, k, i, j, iqw )
4084  calc_cv( cva(k,i,j), qdry, q, k, i, j, iqw, cvdry, tracer_cv )
4085  calc_cp( cpa, qdry, q, k, i, j, iqw, cpdry, tracer_cp )
4086  calc_r( rmoist, qdry, q, k, i, j, iqw, rdry, tracer_r )
4087  r_cva = 1.0_rp / cva(k,i,j)
4088  r_cpa = 1.0_rp / cpa
4089 
4090  ! Coefficient of latent heat release for ssw change by PLCdep and PLRdep
4091  aliqliq = 1.0_rp &
4092  + r_cva*( lhv00 + (cvvap-cl)*tem(k,i,j) )*dqswdtem_rho(k,i,j)
4093  ! Coefficient of latent heat release for ssw change by PLIdep, PLSdep and PLGdep
4094  asolliq = 1.0_rp &
4095  + r_cva*( lhv00 + lhf00 + (cvvap-ci)*tem(k,i,j) )*dqswdtem_rho(k,i,j)
4096  ! Coefficient of latent heat release for ssi change by PLCdep and PLRdep
4097  aliqsol = 1.0_rp &
4098  + r_cva*( lhv00 + (cvvap-cl)*tem(k,i,j) )*dqsidtem_rho(k,i,j)
4099  ! Coefficient of latent heat release for ssi change by PLIdep, PLSdep and PLGdep
4100  asolsol = 1.0_rp &
4101  + r_cva*( lhv00 + lhf00 + (cvvap-ci)*tem(k,i,j) )*dqsidtem_rho(k,i,j)
4102  pdynliq = w * grav * ( r_cpa*dqswdtem_pre(k,i,j) + rho(k,i,j)*dqswdpre_tem(k,i,j) )
4103  pdynsol = w * grav * ( r_cpa*dqsidtem_pre(k,i,j) + rho(k,i,j)*dqsidpre_tem(k,i,j) )
4104  pradliq = -dtdt_rad(k,i,j) * dqswdtem_rho(k,i,j)
4105  pradsol = -dtdt_rad(k,i,j) * dqsidtem_rho(k,i,j)
4106 
4107  ssw_o = ssw
4108  ssi_o = ssi
4109 ! ssw_o = ssw - Pdynliq*(dt_dyn-dt_mp)/qsw(k,i,j) + Pradliq*r_qsw*dt_mp
4110 ! ssi_o = ssi - Pdynsol*(dt_dyn-dt_mp)/qsi(k,i,j) + Pradsol*r_qsi*dt_mp
4111 
4112  acnd = pdynliq + pradliq &
4113  - ( r_taudep_i+r_taudep_s+r_taudep_g ) * ( qsw(k,i,j) - qsi(k,i,j) )
4114  adep = pdynsol + pradsol &
4115  + ( r_taucnd_c+r_taucnd_r ) * ( qsw(k,i,j) - qsi(k,i,j) )
4116  r_taucnd = &
4117  + aliqliq*( r_taucnd_c+r_taucnd_r ) &
4118  + asolliq*( r_taudep_i+r_taudep_s+r_taudep_g )
4119  r_taudep = &
4120  + aliqsol*( r_taucnd_c+r_taucnd_r )&
4121  + asolsol*( r_taudep_i+r_taudep_s+r_taudep_g )
4122 
4123  uplim_cnd = max( rho(k,i,j)*ssw_o*qsw(k,i,j)*r_dt, 0.0_rp )
4124  lowlim_cnd = min( rho(k,i,j)*ssw_o*qsw(k,i,j)*r_dt, 0.0_rp )
4125  if( r_taucnd < r_tau100day )then
4126 ! taucnd = tau100day
4127  pq(i_lcdep,k,i,j) = max(lowlim_cnd, min(uplim_cnd, pq(i_lcdep,k,i,j)*ssw_o ))
4128  pq(i_lrdep,k,i,j) = max(lowlim_cnd, min(uplim_cnd, pq(i_lrdep,k,i,j)*ssw_o ))
4129  pq(i_nrdep,k,i,j) = min(0.0_rp, pq(i_nrdep,k,i,j)*ssw_o )
4130 ! PLR2NR = 0.0_RP
4131  else
4132  taucnd = 1.0_rp/r_taucnd
4133  ! Production term for liquid water content
4134  coef_a_cnd = rho(k,i,j)*acnd*taucnd
4135  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 )
4136  pq(i_lcdep,k,i,j) = coef_a_cnd*r_taucnd_c - coef_b_cnd*r_taucnd_c
4137  plr2nr = pq(i_nrdep,k,i,j)/(pq(i_lrdep,k,i,j)+1.e-30_rp)
4138  pq(i_lrdep,k,i,j) = coef_a_cnd*r_taucnd_r - coef_b_cnd*r_taucnd_r
4139  pq(i_nrdep,k,i,j) = min(0.0_rp, pq(i_lrdep,k,i,j)*plr2nr )
4140  end if
4141 
4142  uplim_dep = max( rho(k,i,j)*ssi_o*qsi(k,i,j)*r_dt, 0.0_rp )
4143  lowlim_dep = min( rho(k,i,j)*ssi_o*qsi(k,i,j)*r_dt, 0.0_rp )
4144  if( r_taudep < r_tau100day )then
4145 ! taudep = tau100day
4146  pq(i_lidep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lidep,k,i,j)*ssi_o ))
4147  pq(i_lsdep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lsdep,k,i,j)*ssi_o ))
4148  pq(i_lgdep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lgdep,k,i,j)*ssi_o ))
4149  pq(i_nidep,k,i,j) = min(0.0_rp, pq(i_nidep,k,i,j)*ssi_o )
4150  pq(i_nsdep,k,i,j) = min(0.0_rp, pq(i_nsdep,k,i,j)*ssi_o )
4151  pq(i_ngdep,k,i,j) = min(0.0_rp, pq(i_ngdep,k,i,j)*ssi_o )
4152  else
4153  taudep = 1.0_rp/r_taudep
4154  ! Production term for ice water content
4155  coef_a_dep = rho(k,i,j)*adep*taudep
4156  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 )
4157  pli2ni = pq(i_nidep,k,i,j)/max(pq(i_lidep,k,i,j),1.e-30_rp)
4158  pls2ns = pq(i_nsdep,k,i,j)/max(pq(i_lsdep,k,i,j),1.e-30_rp)
4159  plg2ng = pq(i_ngdep,k,i,j)/max(pq(i_lgdep,k,i,j),1.e-30_rp)
4160  pq(i_lidep,k,i,j) = coef_a_dep*r_taudep_i - coef_b_dep*r_taudep_i
4161  pq(i_lsdep,k,i,j) = coef_a_dep*r_taudep_s - coef_b_dep*r_taudep_s
4162  pq(i_lgdep,k,i,j) = coef_a_dep*r_taudep_g - coef_b_dep*r_taudep_g
4163  pq(i_nidep,k,i,j) = min(0.0_rp, pq(i_lidep,k,i,j)*pli2ni )
4164  pq(i_nsdep,k,i,j) = min(0.0_rp, pq(i_lsdep,k,i,j)*pls2ns )
4165  pq(i_ngdep,k,i,j) = min(0.0_rp, pq(i_lgdep,k,i,j)*plg2ng )
4166  end if
4167 
4168  sw = 0.5_rp - sign(0.5_rp, pq(i_lcdep,k,i,j)+eps) != 1 for PLCdep<=-eps
4169  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
4170 ! if( PQ(I_LCdep,k,i,j) < -eps )then
4171 ! 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 )
4172 ! else
4173 ! PNCdep = 0.0_RP
4174 ! end if
4175 ! if( PQ(I_LIdep,k,i,j) < -eps )then
4176 ! 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 )
4177 ! else
4178 ! PQ(I_NIdep,k,i,j) = 0.0_RP
4179 ! end if
4180 
4181  !--- evaporation/condensation
4182  r_rvaptem = 1.0_rp/(rvap*wtem(k,i,j))
4183  lvsw = esw(k,i,j)*r_rvaptem
4184  dlvsw = rhoq2(i_qv,k,i,j)-lvsw
4185  dcnd = dt*(pq(i_lcdep,k,i,j)+pq(i_lrdep,k,i,j))
4186 
4187  sw = ( sign(0.5_rp,dcnd) + sign(0.5_rp,dlvsw) ) &
4188  * ( 0.5_rp + sign(0.5_rp,abs(dcnd)-eps) ) ! to avoid zero division
4189  ! sw= 1: always supersaturated
4190  ! sw=-1: always unsaturated
4191  ! sw= 0: partially unsaturated during timestep
4192  fac1 = min(dlvsw*sw,dcnd*sw)*sw / (abs(sw)-1.0_rp+dcnd) & ! sw=1,-1
4193  + 1.0_rp - abs(sw) ! sw=0
4194  dep_dqc = max( dt*pq(i_lcdep,k,i,j)*fac1, &
4195  -rhoq2(i_qc,k,i,j) - 1e30_rp*(sw+1.0_rp) )*abs(sw) != -lc for sw=-1, -inf for sw=1
4196  dep_dqr = max( dt*pq(i_lrdep,k,i,j)*fac1, &
4197  -rhoq2(i_qr,k,i,j) - 1e30_rp*(sw+1.0_rp) )*abs(sw) != -lr for sw=-1, -inf for sw=1
4198 ! if ( (dcnd > eps) .AND. (dlvsw > eps) )then
4199 ! ! always supersaturated
4200 ! fac1 = min(dlvsw,dcnd)/dcnd
4201 ! dep_dqc = dt*PQ(I_LCdep,k,i,j)*fac1
4202 ! dep_dqr = dt*PQ(I_LRdep,k,i,j)*fac1
4203 ! else if( (dcnd < -eps) .AND. (dlvsw < -eps) )then
4204 ! ! always unsaturated
4205 ! fac1 = max( dlvsw,dcnd )/dcnd
4206 ! dep_dqc = max( dt*PQ(I_LCdep,k,i,j)*fac1, -rhoq2(I_QC,k,i,j) )
4207 ! dep_dqr = max( dt*PQ(I_LRdep,k,i,j)*fac1, -rhoq2(I_QR,k,i,j) )
4208 ! else
4209 ! ! partially unsaturated during timestep
4210 ! fac1 = 1.0_RP
4211 ! dep_dqc = 0.0_RP
4212 ! dep_dqr = 0.0_RP
4213 ! end if
4214 
4215  ! evaporation always lose number(always negative).
4216  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
4217  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
4218 
4219  qc_evaporate(k,i,j) = - dep_dnc ! [Add] Y.Sato 15/09/08
4220 
4221  !--- deposition/sublimation
4222  lvsi = esi(k,i,j)*r_rvaptem
4223  ddep = dt*(pq(i_lidep,k,i,j)+pq(i_lsdep,k,i,j)+pq(i_lgdep,k,i,j))
4224  dlvsi = rhoq2(i_qv,k,i,j)-lvsi ! limiter for esi>1.d0
4225 
4226  sw = ( sign(0.5_rp,ddep) + sign(0.5_rp,dlvsi) ) &
4227  * ( 0.5_rp + sign(0.5_rp,abs(ddep)-eps) ) ! to avoid zero division
4228  ! sw= 1: always supersaturated
4229  ! sw=-1: always unsaturated
4230  ! sw= 0: partially unsaturated during timestep
4231  fac2 = min(dlvsi*sw,ddep*sw)*sw / (abs(sw)-1.0_rp+ddep) & ! sw=1,-1
4232  + 1.0_rp - abs(sw) ! sw=0
4233  dep_dqi = max( dt*pq(i_lidep,k,i,j) &
4234  * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), & != fac2 for sw=-1,1, 1 for sw=0
4235  -rhoq2(i_qi,k,i,j) - 1e30_rp*(sw+1.0_rp) ) != -li for sw=-1, -inf for sw=0,1
4236  dep_dqs = max( dt*pq(i_lsdep,k,i,j) &
4237  * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), & != fac2 for sw=-1,1, 1 for sw=0
4238  -rhoq2(i_qs,k,i,j) - 1e30_rp*(sw+1.0_rp) ) != -ls for sw=-1, -inf for sw=0,1
4239  dep_dqg = max( dt*pq(i_lgdep,k,i,j) &
4240  * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), & != fac2 for sw=-1,1, 1 for sw=0
4241  -rhoq2(i_qg,k,i,j) - 1e30_rp*(sw+1.0_rp) ) != -lg for sw=-1, -inf for sw=0,1
4242 ! if ( (ddep > eps) .AND. (dlvsi > eps) )then
4243 ! ! always supersaturated
4244 ! fac2 = min(dlvsi,ddep)/ddep
4245 ! dep_dqi = dt*PQ(I_LIdep,k,i,j)*fac2
4246 ! dep_dqs = dt*PQ(I_LSdep,k,i,j)*fac2
4247 ! dep_dqg = dt*PQ(I_LGdep,k,i,j)*fac2
4248 ! else if ( (ddep < -eps) .AND. (dlvsi < -eps) )then
4249 ! ! always unsaturated
4250 ! fac2 = max(dlvsi,ddep)/ddep
4251 ! dep_dqi = max(dt*PQ(I_LIdep,k,i,j)*fac2, -rhoq2(I_QI,k,i,j) )
4252 ! dep_dqs = max(dt*PQ(I_LSdep,k,i,j)*fac2, -rhoq2(I_QS,k,i,j) )
4253 ! dep_dqg = max(dt*PQ(I_LGdep,k,i,j)*fac2, -rhoq2(I_QG,k,i,j) )
4254 ! else
4255 ! ! partially unsaturated during timestep
4256 ! fac2 = 1.0_RP
4257 ! dep_dqi = dt*PQ(I_LIdep,k,i,j)
4258 ! dep_dqs = dt*PQ(I_LSdep,k,i,j)
4259 ! dep_dqg = dt*PQ(I_LGdep,k,i,j)
4260 ! end if
4261 
4262  ! evaporation always lose number(always negative).
4263  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
4264  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
4265  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
4266 
4267  !--- freezing of cloud drop
4268  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
4269  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
4270  fac3 = ( frz_dqc-eps )/( dt*(pq(i_lchom,k,i,j)+pq(i_lchet,k,i,j))-eps )
4271  fac4 = ( frz_dnc-eps )/( dt*(pq(i_nchom,k,i,j)+pq(i_nchet,k,i,j))-eps )
4272  pq(i_lchom,k,i,j) = fac3*pq(i_lchom,k,i,j)
4273  pq(i_lchet,k,i,j) = fac3*pq(i_lchet,k,i,j)
4274  pq(i_nchom,k,i,j) = fac4*pq(i_nchom,k,i,j)
4275  pq(i_nchet,k,i,j) = fac4*pq(i_nchet,k,i,j)
4276 
4277  !--- melting
4278  ! ice change
4279  mlt_dqi = max( dt*pq(i_limlt,k,i,j), -rhoq2(i_qi,k,i,j)-dep_dqi ) ! negative value
4280  mlt_dni = max( dt*pq(i_nimlt,k,i,j), -rhoq2(i_ni,k,i,j)-dep_dni ) ! negative value
4281  ! snow change
4282  mlt_dqs = max( dt*pq(i_lsmlt,k,i,j), -rhoq2(i_qs,k,i,j)-dep_dqs ) ! negative value
4283  mlt_dns = max( dt*pq(i_nsmlt,k,i,j), -rhoq2(i_ns,k,i,j)-dep_dns ) ! negative value
4284  ! graupel change
4285  mlt_dqg = max( dt*pq(i_lgmlt,k,i,j), -rhoq2(i_qg,k,i,j)-dep_dqg ) ! negative value
4286  mlt_dng = max( dt*pq(i_ngmlt,k,i,j), -rhoq2(i_ng,k,i,j)-dep_dng ) ! negative value
4287 
4288  !--- freezing of larger droplets
4289  frz_dqr = max( dt*(pq(i_lrhet,k,i,j)), min(0.0_rp, -rhoq2(i_qr,k,i,j)-dep_dqr) ) ! negative value
4290  frz_dnr = max( dt*(pq(i_nrhet,k,i,j)), min(0.0_rp, -rhoq2(i_nr,k,i,j)-dep_dnr) ) ! negative value
4291 
4292  fac5 = ( frz_dqr-eps )/( dt*pq(i_lrhet,k,i,j)-eps )
4293  pq(i_lrhet,k,i,j) = fac5*pq(i_lrhet,k,i,j)
4294  fac6 = ( frz_dnr-eps )/( dt*pq(i_nrhet,k,i,j)-eps )
4295  pq(i_nrhet,k,i,j) = fac6*pq(i_nrhet,k,i,j)
4296 
4297  ! water vapor change
4298  drhoqv = -(dep_dqc+dep_dqi+dep_dqs+dep_dqg+dep_dqr)
4299 
4300  rhoq(i_qv,k,i,j) = max(0.0_rp, rhoq(i_qv,k,i,j) + drhoqv )
4301 
4302  rhoe(k,i,j) = rhoe(k,i,j) - lhv * drhoqv
4303 
4304  xi = min(xi_max, max(xi_min, rhoq2(i_qi,k,i,j)/(rhoq2(i_ni,k,i,j)+ni_min) ))
4305  sw = 0.5_rp + sign(0.5_rp,xi-x_sep) ! if (xi>=x_sep) then sw=1 else sw=0
4306  ! sw=1: large ice crystals turn into rain by melting
4307 
4308  ! total cloud change
4309  drhoqc = ( frz_dqc - mlt_dqi*(1.0_rp-sw) + dep_dqc )
4310  drhonc = ( frz_dnc - mlt_dni*(1.0_rp-sw) + dep_dnc )
4311  ! total rain change
4312  drhoqr = ( frz_dqr - mlt_dqg - mlt_dqs - mlt_dqi*sw + dep_dqr )
4313  drhonr = ( frz_dnr - mlt_dng - mlt_dns - mlt_dni*sw + dep_dnr )
4314 
4315  rhoq(i_qc,k,i,j) = max(0.0_rp, rhoq(i_qc,k,i,j) + drhoqc )
4316  rhoq(i_nc,k,i,j) = max(0.0_rp, rhoq(i_nc,k,i,j) + drhonc )
4317  rhoq(i_qr,k,i,j) = max(0.0_rp, rhoq(i_qr,k,i,j) + drhoqr )
4318  rhoq(i_nr,k,i,j) = max(0.0_rp, rhoq(i_nr,k,i,j) + drhonr )
4319 
4320  ! total ice change
4321  drhoqi = (-frz_dqc + mlt_dqi + dep_dqi )
4322  drhoni = (-frz_dnc + mlt_dni + dep_dni )
4323 
4324  rhoq(i_qi,k,i,j) = max(0.0_rp, rhoq(i_qi,k,i,j) + drhoqi )
4325  rhoq(i_ni,k,i,j) = max(0.0_rp, rhoq(i_ni,k,i,j) + drhoni )
4326 
4327  rhoe(k,i,j) = rhoe(k,i,j) + lhf * drhoqi
4328 
4329  ! total snow change
4330  drhoqs = ( mlt_dqs + dep_dqs )
4331  drhons = ( mlt_dns + dep_dns )
4332 
4333  rhoq(i_qs,k,i,j) = max(0.0_rp, rhoq(i_qs,k,i,j) + drhoqs )
4334  rhoq(i_ns,k,i,j) = max(0.0_rp, rhoq(i_ns,k,i,j) + drhons )
4335 
4336  rhoe(k,i,j) = rhoe(k,i,j) + lhf * drhoqs
4337 
4338  ! total graupel change
4339  drhoqg = (-frz_dqr + mlt_dqg + dep_dqg )
4340  drhong = (-frz_dnr + mlt_dng + dep_dng )
4341 
4342  rhoq(i_qg,k,i,j) = max(0.0_rp, rhoq(i_qg,k,i,j) + drhoqg )
4343  rhoq(i_ng,k,i,j) = max(0.0_rp, rhoq(i_ng,k,i,j) + drhong )
4344 
4345  rhoe(k,i,j) = rhoe(k,i,j) + lhf * drhoqg
4346 
4347  !--- update mixing ratio
4348  rrho = 1.0_rp/rho(k,i,j)
4349 
4350  q(k,i,j,i_qv) = rhoq(i_qv,k,i,j) * rrho
4351  q(k,i,j,i_qc) = rhoq(i_qc,k,i,j) * rrho
4352  q(k,i,j,i_qr) = rhoq(i_qr,k,i,j) * rrho
4353  q(k,i,j,i_qi) = rhoq(i_qi,k,i,j) * rrho
4354  q(k,i,j,i_qs) = rhoq(i_qs,k,i,j) * rrho
4355  q(k,i,j,i_qg) = rhoq(i_qg,k,i,j) * rrho
4356  q(k,i,j,i_nc) = rhoq(i_nc,k,i,j) * rrho
4357  q(k,i,j,i_nr) = rhoq(i_nr,k,i,j) * rrho
4358  q(k,i,j,i_ni) = rhoq(i_ni,k,i,j) * rrho
4359  q(k,i,j,i_ns) = rhoq(i_ns,k,i,j) * rrho
4360  q(k,i,j,i_ng) = rhoq(i_ng,k,i,j) * rrho
4361 
4362  calc_qdry( qdry, q, tracer_mass, k, i, j, iqw )
4363  calc_cv( cva(k,i,j), qdry, q, k, i, j, iqw, cvdry, tracer_cv )
4364  calc_r( rmoist, qdry, q, k, i, j, iqw, rdry, tracer_r )
4365  tem(k,i,j) = rhoe(k,i,j) / ( rho(k,i,j) * cva(k,i,j) )
4366  pre(k,i,j) = rho(k,i,j) * rmoist * tem(k,i,j)
4367 
4368  sl_plcdep(i,j) = sl_plcdep(i,j) + dep_dqc*dz(k)*gsgam2(k,i,j)
4369  sl_plrdep(i,j) = sl_plrdep(i,j) + dep_dqr*dz(k)*gsgam2(k,i,j)
4370  sl_pnrdep(i,j) = sl_pnrdep(i,j) + dep_dnr*dz(k)*gsgam2(k,i,j)
4371  end do
4372  end do
4373  end do
4374  profile_stop("sn14_update")
4375 
4376  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_atmos_saturation
module ATMOSPHERE / Saturation adjustment
Definition: scale_atmos_sub_saturation.F90:16
scale_tracer::tracer_r
real(rp), dimension(qa_max), public tracer_r
Definition: scale_tracer.F90:41
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer::tracer_cv
real(rp), dimension(qa_max), public tracer_cv
Definition: scale_tracer.F90:40
scale_atmos_saturation::atmos_saturation_dqsw_dtem_rho
subroutine, public atmos_saturation_dqsw_dtem_rho(dqsdtem, temp, dens)
Definition: scale_atmos_sub_saturation.F90:1178
scale_atmos_saturation::atmos_saturation_dqsw_dtem_dpre
subroutine, public atmos_saturation_dqsw_dtem_dpre(dqsdtem, dqsdpre, temp, pres)
Definition: scale_atmos_sub_saturation.F90:1261
scale_atmos_saturation::atmos_saturation_dqsi_dtem_rho
subroutine, public atmos_saturation_dqsi_dtem_rho(dqsdtem, temp, dens)
Definition: scale_atmos_sub_saturation.F90:1219
scale_tracer::tracer_cp
real(rp), dimension(qa_max), public tracer_cp
Definition: scale_tracer.F90:39
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_atmos_saturation::atmos_saturation_dqsi_dtem_dpre
subroutine, public atmos_saturation_dqsi_dtem_dpre(dqsdtem, dqsdpre, temp, pres)
Definition: scale_atmos_sub_saturation.F90:1310
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_tracer::tracer_mass
real(rp), dimension(qa_max), public tracer_mass
Definition: scale_tracer.F90:42
Here is the call graph for this function:
Here is the caller graph for this function:

◆ mp_negativefilter()

subroutine scale_atmos_phy_mp_sn14::mp_negativefilter ( real(rp), dimension(ka,ia,ja), intent(inout)  DENS,
real(rp), dimension(ka,ia,ja,qa), intent(inout)  QTRC 
)

Definition at line 4382 of file scale_atmos_phy_mp_sn14.F90.

References scale_atmos_hydrometeor::i_nc, scale_atmos_hydrometeor::i_ng, scale_atmos_hydrometeor::i_ni, scale_atmos_hydrometeor::i_nr, scale_atmos_hydrometeor::i_ns, scale_atmos_hydrometeor::i_qc, scale_atmos_hydrometeor::i_qg, scale_atmos_hydrometeor::i_qi, scale_atmos_hydrometeor::i_qr, scale_atmos_hydrometeor::i_qs, scale_atmos_hydrometeor::i_qv, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, scale_prof::prof_rapend(), scale_prof::prof_rapstart(), and scale_tracer::qa.

Referenced by atmos_phy_mp_sn14().

4382  use scale_tracer, only: &
4383  qa
4384  implicit none
4385  real(RP), intent(inout) :: DENS(KA,IA,JA)
4386  real(RP), intent(inout) :: QTRC(KA,IA,JA,QA)
4387 
4388  real(RP) :: diffq(KA,IA,JA)
4389  real(RP) :: r_xmin
4390 
4391  integer :: k, i, j, iq
4392  !---------------------------------------------------------------------------
4393 
4394  call prof_rapstart('MP_filter', 3)
4395 
4396  r_xmin = 1.0_rp / xmin_filter
4397 
4398  ! total hydrometeor (before correction)
4399  do j = js, je
4400  do i = is, ie
4401 
4402  do k = ks, ke
4403  diffq(k,i,j) = qtrc(k,i,j,i_qv) &
4404  + qtrc(k,i,j,i_qc) &
4405  + qtrc(k,i,j,i_qr) &
4406  + qtrc(k,i,j,i_qi) &
4407  + qtrc(k,i,j,i_qs) &
4408  + qtrc(k,i,j,i_qg)
4409  enddo
4410 
4411  ! remove negative value of hydrometeor (mass, number)
4412  do iq = i_qc, i_ng
4413  do k = ks, ke
4414  qtrc(k,i,j,iq) = max(0.0_rp, qtrc(k,i,j,iq))
4415  enddo
4416  enddo
4417 
4418  ! apply correction of hydrometeor to total density
4419  ! [note] mass conservation is broken here to fill rounding error.
4420  do k = ks, ke
4421  dens(k,i,j) = dens(k,i,j) &
4422  * ( 1.0_rp &
4423  + qtrc(k,i,j,i_qv) &
4424  + qtrc(k,i,j,i_qc) &
4425  + qtrc(k,i,j,i_qr) &
4426  + qtrc(k,i,j,i_qi) &
4427  + qtrc(k,i,j,i_qs) &
4428  + qtrc(k,i,j,i_qg) &
4429  - diffq(k,i,j) ) ! after-before
4430  enddo
4431 
4432  ! avoid unrealistical value of number concentration
4433  ! due to numerical diffusion in advection
4434 
4435  do k = ks, ke
4436  if ( qtrc(k,i,j,i_nc) > qtrc(k,i,j,i_qc)*r_xmin ) then
4437  qtrc(k,i,j,i_nc) = qtrc(k,i,j,i_qc)*r_xmin
4438  endif
4439  enddo
4440  do k = ks, ke
4441  if ( qtrc(k,i,j,i_nr) > qtrc(k,i,j,i_qr)*r_xmin ) then
4442  qtrc(k,i,j,i_nr) = qtrc(k,i,j,i_qr)*r_xmin
4443  endif
4444  enddo
4445  do k = ks, ke
4446  if ( qtrc(k,i,j,i_ni) > qtrc(k,i,j,i_qi)*r_xmin ) then
4447  qtrc(k,i,j,i_ni) = qtrc(k,i,j,i_qi)*r_xmin
4448  endif
4449  enddo
4450  do k = ks, ke
4451  if ( qtrc(k,i,j,i_ns) > qtrc(k,i,j,i_qs)*r_xmin ) then
4452  qtrc(k,i,j,i_ns) = qtrc(k,i,j,i_qs)*r_xmin
4453  endif
4454  enddo
4455  do k = ks, ke
4456  if ( qtrc(k,i,j,i_ng) > qtrc(k,i,j,i_qg)*r_xmin ) then
4457  qtrc(k,i,j,i_ng) = qtrc(k,i,j,i_qg)*r_xmin
4458  endif
4459  enddo
4460 
4461  enddo
4462  enddo
4463 
4464  call prof_rapend('MP_filter', 3)
4465 
4466  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_cloudfraction()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_cloudfraction ( real(rp), dimension(ka,ia,ja), intent(out)  cldfrac,
real(rp), dimension (ka,ia,ja,qa), intent(in)  QTRC,
real(rp), intent(in)  mask_criterion 
)

Calculate Cloud Fraction.

Definition at line 4476 of file scale_atmos_phy_mp_sn14.F90.

References scale_atmos_hydrometeor::i_qc, scale_atmos_hydrometeor::i_qg, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, and scale_tracer::qa.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config().

4476  use scale_grid_index
4477  use scale_tracer, only: &
4478  qa
4479  implicit none
4480 
4481  real(RP), intent(out) :: cldfrac(KA,IA,JA)
4482  real(RP), intent(in) :: QTRC (KA,IA,JA,QA)
4483  real(RP), intent(in) :: mask_criterion
4484 
4485  real(RP) :: qhydro
4486  integer :: k, i, j, iq
4487  !---------------------------------------------------------------------------
4488 
4489  do j = js, je
4490  do i = is, ie
4491  do k = ks, ke
4492  qhydro = 0.0_rp
4493  do iq = i_qc, i_qg
4494  qhydro = qhydro + qtrc(k,i,j,iq)
4495  enddo
4496  cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
4497  enddo
4498  enddo
4499  enddo
4500 
4501  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_effectiveradius()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_effectiveradius ( real(rp), dimension (ka,ia,ja,n_hyd), intent(out)  Re,
real(rp), dimension(ka,ia,ja,qa), intent(in)  QTRC0,
real(rp), dimension(ka,ia,ja), intent(in)  DENS0,
real(rp), dimension(ka,ia,ja), intent(in)  TEMP0 
)

Calculate Effective Radius.

Definition at line 4511 of file scale_atmos_phy_mp_sn14.F90.

References scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, scale_atmos_hydrometeor::n_hyd, and scale_tracer::qa.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config().

4511  use scale_grid_index
4512  use scale_tracer, only: &
4513  qa
4514  use scale_atmos_hydrometeor, only: &
4515  n_hyd
4516  implicit none
4517 
4518  real(RP), intent(out) :: Re (KA,IA,JA,N_HYD) ! effective radius [cm]
4519  real(RP), intent(in) :: QTRC0(KA,IA,JA,QA) ! tracer mass concentration [kg/kg]
4520  real(RP), intent(in) :: DENS0(KA,IA,JA) ! density [kg/m3]
4521  real(RP), intent(in) :: TEMP0(KA,IA,JA) ! temperature [K]
4522 
4523  ! mass concentration[kg/m3] and mean particle mass[kg]
4524  real(RP) :: xc(KA,IA,JA)
4525  real(RP) :: xr(KA,IA,JA)
4526  real(RP) :: xi(KA,IA,JA)
4527  real(RP) :: xs(KA,IA,JA)
4528  real(RP) :: xg(KA,IA,JA)
4529  ! diameter of average mass[kg/m3]
4530  real(RP) :: dc_ave(KA,IA,JA)
4531  real(RP) :: dr_ave(KA,IA,JA)
4532  ! radius of average mass
4533  real(RP) :: rc, rr
4534  ! 2nd. and 3rd. order moment of DSD
4535  real(RP) :: ri2m(KA,IA,JA), ri3m(KA,IA,JA)
4536  real(RP) :: rs2m(KA,IA,JA), rs3m(KA,IA,JA)
4537  real(RP) :: rg2m(KA,IA,JA), rg3m(KA,IA,JA)
4538 
4539  real(RP) :: coef_Fuetal1998
4540  ! r2m_min is minimum value(moment of 1 particle with 1 micron)
4541  real(RP), parameter :: r2m_min=1.e-12_rp
4542  real(RP), parameter :: um2cm = 100.0_rp
4543 
4544  real(RP) :: limitsw, zerosw
4545  integer :: k, i, j
4546  !---------------------------------------------------------------------------
4547 
4548  ! mean particle mass[kg]
4549  do j = js, je
4550  do i = is, ie
4551  do k = ks, ke
4552  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) ) )
4553  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) ) )
4554  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) ) )
4555  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) ) )
4556  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) ) )
4557  enddo
4558  enddo
4559  enddo
4560 
4561  ! diameter of average mass : SB06 eq.(32)
4562  do j = js, je
4563  do i = is, ie
4564  do k = ks, ke
4565  dc_ave(k,i,j) = a_m(i_mp_qc) * xc(k,i,j)**b_m(i_mp_qc)
4566  dr_ave(k,i,j) = a_m(i_mp_qr) * xr(k,i,j)**b_m(i_mp_qr)
4567  enddo
4568  enddo
4569  enddo
4570 
4571  ! cloud effective radius
4572  do j = js, je
4573  do i = is, ie
4574  do k = ks, ke
4575  rc = 0.5_rp * dc_ave(k,i,j)
4576  limitsw = 0.5_rp + sign(0.5_rp, rc-rmin_re )
4577  re(k,i,j,i_hc) = coef_re(i_mp_qc) * rc * limitsw * um2cm
4578  enddo
4579  enddo
4580  enddo
4581 
4582  ! rain effective radius
4583  do j = js, je
4584  do i = is, ie
4585  do k = ks, ke
4586  rr = 0.5_rp * dr_ave(k,i,j)
4587  limitsw = 0.5_rp + sign(0.5_rp, rr-rmin_re )
4588  re(k,i,j,i_hr) = coef_re(i_mp_qr) * rr * limitsw * um2cm
4589  enddo
4590  enddo
4591  enddo
4592 
4593  do j = js, je
4594  do i = is, ie
4595  do k = ks, ke
4596  ri2m(k,i,j) = pi * coef_rea2(i_mp_qi) * qtrc0(k,i,j,i_ni) * a_rea2(i_mp_qi) * xi(k,i,j)**b_rea2(i_mp_qi)
4597  rs2m(k,i,j) = pi * coef_rea2(i_mp_qs) * qtrc0(k,i,j,i_ns) * a_rea2(i_mp_qs) * xs(k,i,j)**b_rea2(i_mp_qs)
4598  rg2m(k,i,j) = pi * coef_rea2(i_mp_qg) * qtrc0(k,i,j,i_ng) * a_rea2(i_mp_qg) * xg(k,i,j)**b_rea2(i_mp_qg)
4599  enddo
4600  enddo
4601  enddo
4602 
4603  ! Fu(1996), eq.(3.11) or Fu et al.(1998), eq.(2.5)
4604  coef_fuetal1998 = 3.0_rp / (4.0_rp*rhoi)
4605  do j = js, je
4606  do i = is, ie
4607  do k = ks, ke
4608  ri3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ni) * xi(k,i,j)
4609  rs3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ns) * xs(k,i,j)
4610  rg3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ng) * xg(k,i,j)
4611  enddo
4612  enddo
4613  enddo
4614 
4615  ! ice effective radius
4616  do j = js, je
4617  do i = is, ie
4618  do k = ks, ke
4619  zerosw = 0.5_rp - sign(0.5_rp, ri2m(k,i,j) - r2m_min )
4620  re(k,i,j,i_hi) = ri3m(k,i,j) / ( ri2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
4621  enddo
4622  enddo
4623  enddo
4624 
4625  ! snow effective radius
4626  do j = js, je
4627  do i = is, ie
4628  do k = ks, ke
4629  zerosw = 0.5_rp - sign(0.5_rp, rs2m(k,i,j) - r2m_min )
4630  re(k,i,j,i_hs) = rs3m(k,i,j) / ( rs2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
4631  enddo
4632  enddo
4633  enddo
4634 
4635  ! graupel effective radius
4636  do j = js, je
4637  do i = is, ie
4638  do k = ks, ke
4639  zerosw = 0.5_rp - sign(0.5_rp, rg2m(k,i,j) - r2m_min )
4640  re(k,i,j,i_hg) = rg3m(k,i,j) / ( rg2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
4641  enddo
4642  enddo
4643  enddo
4644 
4645  re(:,:,:,i_hg+1:) = 0.0_rp
4646 
4647  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_atmos_hydrometeor
module Hydrometeor
Definition: scale_atmos_sub_hydrometeor.F90:14
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_mixingratio()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_mixingratio ( real(rp), dimension (ka,ia,ja,n_hyd), intent(out)  Qe,
real(rp), dimension(ka,ia,ja,qa), intent(in)  QTRC0 
)

Calculate mixing ratio of each category.

Definition at line 4654 of file scale_atmos_phy_mp_sn14.F90.

References scale_atmos_hydrometeor::n_hyd, and scale_tracer::qa.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config().

4654  use scale_grid_index
4655  use scale_tracer, only: &
4656  qa
4657  use scale_atmos_hydrometeor, only: &
4658  n_hyd
4659  implicit none
4660 
4661  real(RP), intent(out) :: Qe (KA,IA,JA,N_HYD) ! mixing ratio of each cateory [kg/kg]
4662  real(RP), intent(in) :: QTRC0(KA,IA,JA,QA) ! tracer mass concentration [kg/kg]
4663 
4664  integer :: ihydro, iqa
4665  !---------------------------------------------------------------------------
4666 
4667 
4668  qe(:,:,:,i_hc) = qtrc0(:,:,:,i_qc)
4669  qe(:,:,:,i_hr) = qtrc0(:,:,:,i_qr)
4670  qe(:,:,:,i_hi) = qtrc0(:,:,:,i_qi)
4671  qe(:,:,:,i_hs) = qtrc0(:,:,:,i_qs)
4672  qe(:,:,:,i_hg) = qtrc0(:,:,:,i_qg)
4673  qe(:,:,:,i_hg+1:) = 0.0_rp
4674 
4675  return
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_atmos_hydrometeor
module Hydrometeor
Definition: scale_atmos_sub_hydrometeor.F90:14
Here is the caller graph for this function:

Variable Documentation

◆ qa_mp

integer, parameter, public scale_atmos_phy_mp_sn14::qa_mp = 11

Definition at line 137 of file scale_atmos_phy_mp_sn14.F90.

Referenced by atmos_phy_mp_sn14(), and atmos_phy_mp_sn14_config().

137  integer, public, parameter :: QA_MP = 11

◆ atmos_phy_mp_sn14_name

character(len=h_short), dimension(qa_mp), target, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_name

Definition at line 139 of file scale_atmos_phy_mp_sn14.F90.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config(), and atmos_phy_mp_sn14_config().

139  character(len=H_SHORT), public, target :: ATMOS_PHY_MP_sn14_NAME(QA_MP)

◆ atmos_phy_mp_sn14_desc

character(len=h_mid), dimension(qa_mp), target, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_desc

Definition at line 140 of file scale_atmos_phy_mp_sn14.F90.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config(), and atmos_phy_mp_sn14_config().

140  character(len=H_MID) , public, target :: ATMOS_PHY_MP_sn14_DESC(QA_MP)

◆ atmos_phy_mp_sn14_unit

character(len=h_short), dimension(qa_mp), target, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_unit

Definition at line 141 of file scale_atmos_phy_mp_sn14.F90.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config(), and atmos_phy_mp_sn14_config().

141  character(len=H_SHORT), public, target :: ATMOS_PHY_MP_sn14_UNIT(QA_MP)

◆ atmos_phy_mp_sn14_dens

real(rp), dimension(n_hyd), target, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_dens

Definition at line 143 of file scale_atmos_phy_mp_sn14.F90.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config(), and atmos_phy_mp_sn14_setup().

143  real(RP), public, target :: ATMOS_PHY_MP_sn14_DENS(N_HYD) ! hydrometeor density [kg/m3]=[g/L]
Generated by   doxygen 1.8.14