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_setup (KA, IA, JA)
 ATMOS_PHY_MP_sn14_setup setup. More...
 
subroutine, public atmos_phy_mp_sn14_tendency (KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, W, QTRC, PRES, TEMP, Qdry, CPtot, CVtot, CCN, dt, cz, fz, RHOQ_t, RHOE_t, CPtot_t, CVtot_t, EVAPORATE, flg_lt, d0_crg, v0_crg, dqcrg, beta_crg, QTRC_crg, QSPLT_in, Sarea, RHOQcrg_t)
 ATMOS_PHY_MP_sn14_tendency calculate tendency. More...
 
subroutine, public atmos_phy_mp_sn14_cloud_fraction (KA, KS, KE, IA, IS, IE, JA, JS, JE, QTRC, mask_criterion, cldfrac)
 ATMOS_PHY_MP_sn14_cloud_fraction Calculate Cloud Fraction. More...
 
subroutine, public atmos_phy_mp_sn14_effective_radius (KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS0, TEMP0, QTRC0, Re)
 ATMOS_PHY_MP_sn14_effective_radius Calculate Effective Radius. More...
 
subroutine, public atmos_phy_mp_sn14_qtrc2qhyd (KA, KS, KE, IA, IS, IE, JA, JS, JE, QTRC0, Qe)
 ATMOS_PHY_MP_sn14_qtrc2qhyd Calculate mass ratio of each category. More...
 
subroutine, public atmos_phy_mp_sn14_qtrc2nhyd (KA, KS, KE, IA, IS, IE, JA, JS, JE, QTRC0, Ne)
 Calculate number concentration of each category. More...
 
subroutine, public atmos_phy_mp_sn14_qhyd2qtrc (KA, KS, KE, IA, IS, IE, JA, JS, JE, Qe, QTRC, QNUM)
 
subroutine, public atmos_phy_mp_sn14_terminal_velocity (KA, KS, KE, DENS, TEMP, RHOQ, PRES, vterm)
 ATMOS_PHY_MP_sn14_terminal_velocity Calculate terminal velocity. More...
 
subroutine debug_tem (KA, KS, KE, point, i, j, tem, rho, pre, qv)
 
subroutine nucleation (KA, KS, KE, cz, fz, w, rho, tem, pre, qdry, rhoq, cpa, dTdt_rad, qke, CCN, nc_uplim_d, dt, PQ)
 
subroutine ice_multiplication (KA, KS, KE, flg_lt, Pac, tem, rhoq, rhoq_crg, xq, PQ, Pcrg1)
 
subroutine mixed_phase_collection (KA, KS, KE, flg_lt, d0_crg, v0_crg, beta_crg, dqcrg, wtem, rhoq, rhoq_crg, xq, dq_xave, vt_xave, PQ, Pcrg1, Pcrg2, Pac)
 
subroutine aut_acc_slc_brk (KA, KS, KE, flg_lt, rhoq, rhoq_crg, xq, dq_xave, rho, PQ, Pcrg)
 
subroutine freezing_water (KA, KS, KE, dt, rhoq, xq, tem, PQ)
 
subroutine cross_section (Crs, KA, KS, KE, QA_MP, QTRC0, DENS0)
 Calculate Cross Section. More...
 

Variables

integer, parameter, public qa_mp = 11
 
integer, parameter, public atmos_phy_mp_sn14_ntracers = QA_MP
 
integer, parameter, public atmos_phy_mp_sn14_nwaters = 2
 
integer, parameter, public atmos_phy_mp_sn14_nices = 3
 
character(len=h_short), dimension(qa_mp), parameter, public atmos_phy_mp_sn14_tracer_names = (/ 'QV', 'QC', 'QR', 'QI', 'QS', 'QG', 'NC', 'NR', 'NI', 'NS', 'NG' /)
 
character(len=h_mid), dimension(qa_mp), parameter, public atmos_phy_mp_sn14_tracer_descriptions = (/ 'Ratio of Water Vapor mass to total mass (Specific humidity)', 'Ratio of Cloud Water mass to total mass ', 'Ratio of Rain Water mass to total mass ', 'Ratio of Cloud Ice mass ratio to total mass ', 'Ratio of Snow miass ratio to total mass ', 'Ratio of Graupel mass ratio to total mass ', 'Cloud Water Number Density ', 'Rain Water Number Density ', 'Cloud Ice Number Density ', 'Snow Number Density ', 'Graupel Number Density '/)
 
character(len=h_short), dimension(qa_mp), parameter, public atmos_phy_mp_sn14_tracer_units = (/ 'kg/kg ', 'kg/kg ', 'kg/kg ', 'kg/kg ', 'kg/kg ', 'kg/kg ', 'num/kg', 'num/kg', 'num/kg', 'num/kg', 'num/kg' /)
 

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
NAMELIST
  • PARAM_ATMOS_PHY_MP_SN14
    nametypedefault valuecomment
    MP_DOAUTOCONVERSION logical .true.
    MP_SSW_LIM real(RP) 1.E+1_RP
    MP_COUPLE_AEROSOL logical .false. apply CCN effect?

  • PARAM_ATMOS_PHY_MP_SN14_init
    nametypedefault valuecomment
    OPT_DEBUG 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

  • PARAM_ATMOS_PHY_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

  • PARAM_ATMOS_PHY_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.

  • PARAM_ATMOS_PHY_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]

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

History Output
No history output

Function/Subroutine Documentation

◆ atmos_phy_mp_sn14_setup()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_setup ( integer, intent(in)  KA,
integer, intent(in)  IA,
integer, intent(in)  JA 
)

ATMOS_PHY_MP_sn14_setup setup.

Definition at line 597 of file scale_atmos_phy_mp_sn14.F90.

597  use scale_prc, only: &
598  prc_abort
599  implicit none
600 
601  integer, intent(in) :: KA
602  integer, intent(in) :: IA
603  integer, intent(in) :: JA
604 
605  namelist / param_atmos_phy_mp_sn14 / &
606  mp_doautoconversion, &
607  mp_ssw_lim, &
608  mp_couple_aerosol
609 
610  integer :: ierr
611  !---------------------------------------------------------------------------
612 
613  log_newline
614  log_info("ATMOS_PHY_MP_sn14_setup",*) 'Setup'
615  log_info("ATMOS_PHY_MP_sn14_setup",*) 'Seiki and Nakajima (2014) 2-moment bulk 6 category'
616 
617  !--- read namelist
618  rewind(io_fid_conf)
619  read(io_fid_conf,nml=param_atmos_phy_mp_sn14,iostat=ierr)
620  if( ierr < 0 ) then !--- missing
621  log_info("ATMOS_PHY_MP_sn14_setup",*) 'Not found namelist. Default used.'
622  elseif( ierr > 0 ) then !--- fatal error
623  log_error("ATMOS_PHY_MP_sn14_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SN14. Check!'
624  call prc_abort
625  endif
626  log_nml(param_atmos_phy_mp_sn14)
627 
628  wlabel(1) = "CLOUD"
629  wlabel(2) = "RAIN"
630  wlabel(3) = "ICE"
631  wlabel(4) = "SNOW"
632  wlabel(5) = "GRAUPEL"
633 
634  call mp_sn14_init
635 
636  allocate(nc_uplim_d(1,ia,ja))
637  nc_uplim_d(:,:,:) = 150.e6_rp
638 
639  return

References scale_io::io_fid_conf, and scale_prc::prc_abort().

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_setup().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_tendency()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_tendency ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  IS,
integer, intent(in)  IE,
integer, intent(in)  JA,
integer, intent(in)  JS,
integer, intent(in)  JE,
real(rp), dimension (ka,ia,ja), intent(in)  DENS,
real(rp), dimension (ka,ia,ja), intent(in)  W,
real(rp), dimension (ka,ia,ja,qa_mp), intent(in)  QTRC,
real(rp), dimension(ka,ia,ja), intent(in)  PRES,
real(rp), dimension(ka,ia,ja), intent(in)  TEMP,
real(rp), dimension(ka,ia,ja), intent(in)  Qdry,
real(rp), dimension(ka,ia,ja), intent(in)  CPtot,
real(rp), dimension(ka,ia,ja), intent(in)  CVtot,
real(rp), dimension (ka,ia,ja), intent(in)  CCN,
real(dp), intent(in)  dt,
real(rp), dimension( ka,ia,ja), intent(in)  cz,
real(rp), dimension(0:ka,ia,ja), intent(in)  fz,
real(rp), dimension (ka,ia,ja,qa_mp), intent(out)  RHOQ_t,
real(rp), dimension (ka,ia,ja), intent(out)  RHOE_t,
real(rp), dimension(ka,ia,ja), intent(out)  CPtot_t,
real(rp), dimension(ka,ia,ja), intent(out)  CVtot_t,
real(rp), dimension(ka,ia,ja), intent(out)  EVAPORATE,
logical, intent(in), optional  flg_lt,
real(rp), intent(in), optional  d0_crg,
real(rp), intent(in), optional  v0_crg,
real(rp), dimension(ka,ia,ja), intent(in), optional  dqcrg,
real(rp), dimension(ka,ia,ja), intent(in), optional  beta_crg,
real(rp), dimension(ka,ia,ja,hydro_max), intent(in), optional  QTRC_crg,
real(rp), dimension(ka,ia,ja,3), intent(out), optional  QSPLT_in,
real(rp), dimension(ka,ia,ja,hydro_max), intent(out), optional  Sarea,
real(rp), dimension(ka,ia,ja,hydro_max), intent(out), optional  RHOQcrg_t 
)

ATMOS_PHY_MP_sn14_tendency calculate tendency.

Definition at line 672 of file scale_atmos_phy_mp_sn14.F90.

672  implicit none
673 
674  integer, intent(in) :: KA, KS, KE
675  integer, intent(in) :: IA, IS, IE
676  integer, intent(in) :: JA, JS, JE
677 
678  real(RP), intent(in) :: DENS (KA,IA,JA)
679  real(RP), intent(in) :: W (KA,IA,JA)
680  real(RP), intent(in) :: QTRC (KA,IA,JA,QA_MP)
681  real(RP), intent(in) :: PRES(KA,IA,JA)
682  real(RP), intent(in) :: TEMP(KA,IA,JA)
683  real(RP), intent(in) :: Qdry(KA,IA,JA)
684  real(RP), intent(in) :: CPtot(KA,IA,JA)
685  real(RP), intent(in) :: CVtot(KA,IA,JA)
686  real(RP), intent(in) :: CCN (KA,IA,JA)
687  real(DP), intent(in) :: dt
688  real(RP), intent(in) :: cz( KA,IA,JA)
689  real(RP), intent(in) :: fz(0:KA,IA,JA)
690 
691  real(RP), intent(out) :: RHOQ_t (KA,IA,JA,QA_MP)
692  real(RP), intent(out) :: RHOE_t (KA,IA,JA)
693  real(RP), intent(out) :: CPtot_t(KA,IA,JA)
694  real(RP), intent(out) :: CVtot_t(KA,IA,JA)
695  real(RP), intent(out) :: EVAPORATE(KA,IA,JA) !--- number of evaporated cloud [/m3]
696 
697  ! Optional for Lightning
698  logical, intent(in), optional :: flg_lt
699  real(RP), intent(in), optional :: d0_crg, v0_crg
700  real(RP), intent(in), optional :: dqcrg(KA,IA,JA)
701  real(RP), intent(in), optional :: beta_crg(KA,IA,JA)
702  real(RP), intent(in), optional :: QTRC_crg(KA,IA,JA,HYDRO_MAX)
703  real(RP), intent(out), optional :: QSPLT_in(KA,IA,JA,3)
704  real(RP), intent(out), optional :: Sarea(KA,IA,JA,HYDRO_MAX)
705  real(RP), intent(out), optional :: RHOQcrg_t(KA,IA,JA,HYDRO_MAX)
706  !---------------------------------------------------------------------------
707 
708  log_progress(*) 'atmosphere / physics / microphysics / SN14'
709 
710 #ifdef PROFILE_FIPP
711  call fipp_start()
712 #endif
713 
714  !##### MP Main #####
715  call mp_sn14 ( &
716  ka, ks, ke, ia, is, ie, ja, js, je, &
717  dens(:,:,:), w(:,:,:), qtrc(:,:,:,:), pres(:,:,:), temp(:,:,:), & ! (in)
718  qdry(:,:,:), cptot(:,:,:), cvtot(:,:,:), ccn(:,:,:), & ! (in)
719  real(dt,RP), cz(:,:,:), fz(:,:,:), & ! (in)
720  RHOQ_t(:,:,:,:), RHOE_t(:,:,:), CPtot_t(:,:,:), CVtot_t(:,:,:), & ! (out)
721  EVAPORATE(:,:,:), & ! (out)
722  flg_lt, d0_crg, v0_crg, dqcrg(:,:,:), beta_crg(:,:,:), & ! (optional in)
723  QTRC_crg(:,:,:,:), & ! (optional in)
724  QSPLT_in(:,:,:,:), Sarea(:,:,:,:), RHOQcrg_t(:,:,:,:) ) ! (optional out)
725 
726 #ifdef PROFILE_FIPP
727  call fipp_stop()
728 #endif
729 
730  return

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_calc_tendency().

Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_cloud_fraction()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_cloud_fraction ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  IS,
integer, intent(in)  IE,
integer, intent(in)  JA,
integer, intent(in)  JS,
integer, intent(in)  JE,
real(rp), dimension (ka,ia,ja,qa_mp-1), intent(in)  QTRC,
real(rp), intent(in)  mask_criterion,
real(rp), dimension(ka,ia,ja), intent(out)  cldfrac 
)

ATMOS_PHY_MP_sn14_cloud_fraction Calculate Cloud Fraction.

Definition at line 742 of file scale_atmos_phy_mp_sn14.F90.

742  implicit none
743  integer, intent(in) :: KA, KS, KE
744  integer, intent(in) :: IA, IS, IE
745  integer, intent(in) :: JA, JS, JE
746 
747  real(RP), intent(in) :: QTRC (KA,IA,JA,QA_MP-1)
748  real(RP), intent(in) :: mask_criterion
749 
750  real(RP), intent(out) :: cldfrac(KA,IA,JA)
751 
752  real(RP) :: qhydro
753  integer :: k, i, j, iq
754  !---------------------------------------------------------------------------
755 
756  !$omp parallel do &
757  !$omp private(qhydro)
758  do j = js, je
759  do i = is, ie
760  do k = ks, ke
761  qhydro = 0.0_rp
762  do iq = 1, qa_mp-1
763  qhydro = qhydro + qtrc(k,i,j,iq)
764  enddo
765  cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
766  enddo
767  enddo
768  enddo
769 
770  return

References qa_mp.

Referenced by mod_atmos_phy_mp_vars::atmos_phy_mp_vars_get_diagnostic().

Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_effective_radius()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_effective_radius ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  IS,
integer, intent(in)  IE,
integer, intent(in)  JA,
integer, intent(in)  JS,
integer, intent(in)  JE,
real(rp), dimension(ka,ia,ja), intent(in)  DENS0,
real(rp), dimension(ka,ia,ja), intent(in)  TEMP0,
real(rp), dimension(ka,ia,ja,i_qc:i_ng), intent(in)  QTRC0,
real(rp), dimension (ka,ia,ja,n_hyd), intent(out)  Re 
)

ATMOS_PHY_MP_sn14_effective_radius Calculate Effective Radius.

Definition at line 780 of file scale_atmos_phy_mp_sn14.F90.

780  use scale_atmos_hydrometeor, only: &
781  n_hyd, &
782  i_hc, &
783  i_hr, &
784  i_hi, &
785  i_hs, &
786  i_hg, &
787  i_hh
788  implicit none
789  integer, intent(in) :: KA, KS, KE
790  integer, intent(in) :: IA, IS, IE
791  integer, intent(in) :: JA, JS, JE
792 
793  real(RP), intent(in) :: DENS0(KA,IA,JA) ! density [kg/m3]
794  real(RP), intent(in) :: TEMP0(KA,IA,JA) ! temperature [K]
795  real(RP), intent(in) :: QTRC0(KA,IA,JA,I_QC:I_NG) ! tracer mass concentration [kg/kg]
796 
797  real(RP), intent(out) :: Re (KA,IA,JA,N_HYD) ! effective radius [cm]
798 
799  ! mass concentration[kg/m3] and mean particle mass[kg]
800  real(RP) :: xc(KA)
801  real(RP) :: xr(KA)
802  real(RP) :: xi(KA)
803  real(RP) :: xs(KA)
804  real(RP) :: xg(KA)
805  ! diameter of average mass[kg/m3]
806  real(RP) :: dc_ave(KA)
807  real(RP) :: dr_ave(KA)
808  ! radius of average mass
809  real(RP) :: rc, rr
810  ! 2nd. and 3rd. order moment of DSD
811  real(RP) :: ri2m(KA), ri3m(KA)
812  real(RP) :: rs2m(KA), rs3m(KA)
813  real(RP) :: rg2m(KA), rg3m(KA)
814 
815  real(RP), parameter :: coef_Fuetal1998 = 3.0_rp / (4.0_rp*rhoi)
816 
817  ! r2m_min is minimum value(moment of 1 particle with 1 micron)
818  real(RP), parameter :: r2m_min=1.e-12_rp
819  real(RP), parameter :: um2cm = 100.0_rp
820 
821  real(RP) :: limitsw, zerosw
822  integer :: k, i, j
823  !---------------------------------------------------------------------------
824 
825  !$omp parallel do &
826  !$omp private(xc,xr,xi,xs,xg,dc_ave,dr_ave,rc,rr,ri2m,ri3m,rs2m,rs3m,rg2m,rg3m, &
827  !$omp limitsw,zerosw)
828  do j = js, je
829  do i = is, ie
830 
831  ! mean particle mass[kg]
832  do k = ks, ke
833  xc(k) = min( xc_max, max( xc_min, dens0(k,i,j)*qtrc0(k,i,j,i_qc)/(qtrc0(k,i,j,i_nc)+nc_min) ) )
834  xr(k) = min( xr_max, max( xr_min, dens0(k,i,j)*qtrc0(k,i,j,i_qr)/(qtrc0(k,i,j,i_nr)+nr_min) ) )
835  xi(k) = min( xi_max, max( xi_min, dens0(k,i,j)*qtrc0(k,i,j,i_qi)/(qtrc0(k,i,j,i_ni)+ni_min) ) )
836  xs(k) = min( xs_max, max( xs_min, dens0(k,i,j)*qtrc0(k,i,j,i_qs)/(qtrc0(k,i,j,i_ns)+ns_min) ) )
837  xg(k) = min( xg_max, max( xg_min, dens0(k,i,j)*qtrc0(k,i,j,i_qg)/(qtrc0(k,i,j,i_ng)+ng_min) ) )
838  enddo
839 
840  ! diameter of average mass : SB06 eq.(32)
841  do k = ks, ke
842  dc_ave(k) = a_m(i_mp_qc) * xc(k)**b_m(i_mp_qc)
843  dr_ave(k) = a_m(i_mp_qr) * xr(k)**b_m(i_mp_qr)
844  enddo
845 
846  ! cloud effective radius
847  do k = ks, ke
848  rc = 0.5_rp * dc_ave(k)
849  limitsw = 0.5_rp + sign(0.5_rp, rc-rmin_re )
850  re(k,i,j,i_hc) = coef_re(i_mp_qc) * rc * limitsw * um2cm
851  enddo
852 
853  ! rain effective radius
854  do k = ks, ke
855  rr = 0.5_rp * dr_ave(k)
856  limitsw = 0.5_rp + sign(0.5_rp, rr-rmin_re )
857  re(k,i,j,i_hr) = coef_re(i_mp_qr) * rr * limitsw * um2cm
858  enddo
859 
860  do k = ks, ke
861  ri2m(k) = pi * coef_rea2(i_mp_qi) * qtrc0(k,i,j,i_ni) * a_rea2(i_mp_qi) * xi(k)**b_rea2(i_mp_qi)
862  rs2m(k) = pi * coef_rea2(i_mp_qs) * qtrc0(k,i,j,i_ns) * a_rea2(i_mp_qs) * xs(k)**b_rea2(i_mp_qs)
863  rg2m(k) = pi * coef_rea2(i_mp_qg) * qtrc0(k,i,j,i_ng) * a_rea2(i_mp_qg) * xg(k)**b_rea2(i_mp_qg)
864  enddo
865 
866  ! Fu(1996), eq.(3.11) or Fu et al.(1998), eq.(2.5)
867  do k = ks, ke
868  ri3m(k) = coef_fuetal1998 * qtrc0(k,i,j,i_ni) * xi(k)
869  rs3m(k) = coef_fuetal1998 * qtrc0(k,i,j,i_ns) * xs(k)
870  rg3m(k) = coef_fuetal1998 * qtrc0(k,i,j,i_ng) * xg(k)
871  enddo
872 
873  ! ice effective radius
874  do k = ks, ke
875  zerosw = 0.5_rp - sign(0.5_rp, ri2m(k) - r2m_min )
876  re(k,i,j,i_hi) = ri3m(k) / ( ri2m(k) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
877  enddo
878 
879  ! snow effective radius
880  do k = ks, ke
881  zerosw = 0.5_rp - sign(0.5_rp, rs2m(k) - r2m_min )
882  re(k,i,j,i_hs) = rs3m(k) / ( rs2m(k) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
883  enddo
884 
885  ! graupel effective radius
886  do k = ks, ke
887  zerosw = 0.5_rp - sign(0.5_rp, rg2m(k) - r2m_min )
888  re(k,i,j,i_hg) = rg3m(k) / ( rg2m(k) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
889  enddo
890 
891  do k = ks, ke
892  re(k,i,j,i_hh) = 0.0_rp
893  end do
894 
895  enddo
896  enddo
897 
898 
899  return

References scale_atmos_hydrometeor::i_hc, scale_atmos_hydrometeor::i_hg, scale_atmos_hydrometeor::i_hh, scale_atmos_hydrometeor::i_hi, scale_atmos_hydrometeor::i_hr, scale_atmos_hydrometeor::i_hs, and scale_atmos_hydrometeor::n_hyd.

Referenced by mod_atmos_phy_mp_vars::atmos_phy_mp_vars_get_diagnostic().

Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_qtrc2qhyd()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_qtrc2qhyd ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  IS,
integer, intent(in)  IE,
integer, intent(in)  JA,
integer, intent(in)  JS,
integer, intent(in)  JE,
real(rp), dimension(ka,ia,ja,qa_mp-1), intent(in)  QTRC0,
real(rp), dimension (ka,ia,ja,n_hyd), intent(out)  Qe 
)

ATMOS_PHY_MP_sn14_qtrc2qhyd Calculate mass ratio of each category.

Definition at line 909 of file scale_atmos_phy_mp_sn14.F90.

909  use scale_atmos_hydrometeor, only: &
910  n_hyd, &
911  i_hc, &
912  i_hr, &
913  i_hi, &
914  i_hs, &
915  i_hg, &
916  i_hh
917  implicit none
918  integer, intent(in) :: KA, KS, KE
919  integer, intent(in) :: IA, IS, IE
920  integer, intent(in) :: JA, JS, JE
921 
922  real(RP), intent(in) :: QTRC0(KA,IA,JA,QA_MP-1) ! tracer mass concentration [kg/kg]
923 
924  real(RP), intent(out) :: Qe (KA,IA,JA,N_HYD) ! mixing ratio of each cateory [kg/kg]
925 
926  integer :: k, i, j
927  !---------------------------------------------------------------------------
928 
929 !OCL XFILL
930  !$omp parallel do
931  do j = js, je
932  do i = is, ie
933  do k = ks, ke
934  qe(k,i,j,i_hc) = qtrc0(k,i,j,i_mp_qc)
935  qe(k,i,j,i_hr) = qtrc0(k,i,j,i_mp_qr)
936  qe(k,i,j,i_hi) = qtrc0(k,i,j,i_mp_qi)
937  qe(k,i,j,i_hs) = qtrc0(k,i,j,i_mp_qs)
938  qe(k,i,j,i_hg) = qtrc0(k,i,j,i_mp_qg)
939  qe(k,i,j,i_hh) = 0.0_rp
940  end do
941  end do
942  end do
943 
944  return

References scale_atmos_hydrometeor::i_hc, scale_atmos_hydrometeor::i_hg, scale_atmos_hydrometeor::i_hh, scale_atmos_hydrometeor::i_hi, scale_atmos_hydrometeor::i_hr, scale_atmos_hydrometeor::i_hs, and scale_atmos_hydrometeor::n_hyd.

Referenced by mod_atmos_phy_mp_vars::atmos_phy_mp_vars_get_diagnostic().

Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_qtrc2nhyd()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_qtrc2nhyd ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  IS,
integer, intent(in)  IE,
integer, intent(in)  JA,
integer, intent(in)  JS,
integer, intent(in)  JE,
real(rp), dimension(ka,ia,ja,qa_mp-1), intent(in)  QTRC0,
real(rp), dimension (ka,ia,ja,n_hyd), intent(out)  Ne 
)

Calculate number concentration of each category.

Definition at line 952 of file scale_atmos_phy_mp_sn14.F90.

952  use scale_atmos_hydrometeor, only: &
953  n_hyd, &
954  i_hc, &
955  i_hr, &
956  i_hi, &
957  i_hs, &
958  i_hg, &
959  i_hh
960  implicit none
961  integer, intent(in) :: KA, KS, KE
962  integer, intent(in) :: IA, IS, IE
963  integer, intent(in) :: JA, JS, JE
964 
965  real(RP), intent(in) :: QTRC0(KA,IA,JA,QA_MP-1) ! tracer mass concentration [kg/kg]
966 
967  real(RP), intent(out) :: Ne (KA,IA,JA,N_HYD) ! number density of each cateory [1/m3]
968 
969  integer :: k, i, j
970  !---------------------------------------------------------------------------
971 
972 !OCL XFILL
973  !$omp parallel do
974  do j = js, je
975  do i = is, ie
976  do k = ks, ke
977  ne(k,i,j,i_hc) = qtrc0(k,i,j,i_mp_nc)
978  ne(k,i,j,i_hr) = qtrc0(k,i,j,i_mp_nr)
979  ne(k,i,j,i_hi) = qtrc0(k,i,j,i_mp_ni)
980  ne(k,i,j,i_hs) = qtrc0(k,i,j,i_mp_ns)
981  ne(k,i,j,i_hg) = qtrc0(k,i,j,i_mp_ng)
982  ne(k,i,j,i_hh) = 0.0_rp
983  end do
984  end do
985  end do
986 
987  return

References scale_atmos_hydrometeor::i_hc, scale_atmos_hydrometeor::i_hg, scale_atmos_hydrometeor::i_hh, scale_atmos_hydrometeor::i_hi, scale_atmos_hydrometeor::i_hr, scale_atmos_hydrometeor::i_hs, and scale_atmos_hydrometeor::n_hyd.

Referenced by mod_atmos_phy_mp_vars::atmos_phy_mp_vars_get_diagnostic().

Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_qhyd2qtrc()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_qhyd2qtrc ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  IS,
integer, intent(in)  IE,
integer, intent(in)  JA,
integer, intent(in)  JS,
integer, intent(in)  JE,
real(rp), dimension(ka,ia,ja,n_hyd), intent(in)  Qe,
real(rp), dimension(ka,ia,ja,qa_mp-1), intent(out)  QTRC,
real(rp), dimension(ka,ia,ja,n_hyd), intent(in), optional  QNUM 
)

Definition at line 995 of file scale_atmos_phy_mp_sn14.F90.

995  use scale_const, only: &
996  pi => const_pi, &
997  undef => const_undef
998  use scale_atmos_hydrometeor, only: &
999  n_hyd, &
1000  i_hc, &
1001  i_hr, &
1002  i_hi, &
1003  i_hs, &
1004  i_hg, &
1005  i_hh
1006  implicit none
1007  integer, intent(in) :: KA, KS, KE
1008  integer, intent(in) :: IA, IS, IE
1009  integer, intent(in) :: JA, JS, JE
1010 
1011  real(RP), intent(in) :: Qe(KA,IA,JA,N_HYD) ! mass ratio of each cateory [kg/kg]
1012 
1013  real(RP), intent(out) :: QTRC(KA,IA,JA,QA_MP-1) ! tracer mass concentration [kg/kg]
1014 
1015  real(RP), intent(in), optional :: QNUM(KA,IA,JA,N_HYD)
1016 
1017  real(RP), parameter :: Dc = 20.e-6_rp ! typical particle diameter for cloud [m]
1018  real(RP), parameter :: Dr = 200.e-6_rp ! typical particle diameter for rain [m]
1019  real(RP), parameter :: Di = 80.e-6_rp ! typical particle diameter for ice [m]
1020  real(RP), parameter :: Ds = 80.e-6_rp ! typical particle diameter for snow [m]
1021  real(RP), parameter :: Dg = 200.e-6_rp ! typical particle diameter for grapel [m]
1022  real(RP), parameter :: RHOw = 1000.0_rp ! typical density for warm particles [kg/m3]
1023  real(RP), parameter :: RHOf = 100.0_rp ! typical density for frozen particles [kg/m3]
1024  real(RP), parameter :: RHOg = 400.0_rp ! typical density for grapel particles [kg/m3]
1025  real(RP), parameter :: b = 3.0_rp ! assume spherical form
1026 
1027  real(RP) :: piov6
1028 
1029  integer :: k, i, j
1030  !---------------------------------------------------------------------------
1031 
1032 
1033 !OCL XFILL
1034  !$omp parallel do
1035  do j = js, je
1036  do i = is, ie
1037  do k = ks, ke
1038  qtrc(k,i,j,i_mp_qc) = qe(k,i,j,i_hc)
1039  end do
1040  end do
1041  end do
1042 
1043 !OCL XFILL
1044  !$omp parallel do
1045  do j = js, je
1046  do i = is, ie
1047  do k = ks, ke
1048  qtrc(k,i,j,i_mp_qr) = qe(k,i,j,i_hr)
1049  end do
1050  end do
1051  end do
1052 
1053 !OCL XFILL
1054  !$omp parallel do
1055  do j = js, je
1056  do i = is, ie
1057  do k = ks, ke
1058  qtrc(k,i,j,i_mp_qi) = qe(k,i,j,i_hi)
1059  end do
1060  end do
1061  end do
1062 
1063 !OCL XFILL
1064  !$omp parallel do
1065  do j = js, je
1066  do i = is, ie
1067  do k = ks, ke
1068  qtrc(k,i,j,i_mp_qs) = qe(k,i,j,i_hs)
1069  end do
1070  end do
1071  end do
1072 
1073 !OCL XFILL
1074  !$omp parallel do
1075  do j = js, je
1076  do i = is, ie
1077  do k = ks, ke
1078  qtrc(k,i,j,i_mp_qg) = qe(k,i,j,i_hg) + qe(k,i,j,i_hh)
1079  end do
1080  end do
1081  end do
1082 
1083 
1084  piov6 = pi / 6.0_rp
1085 
1086  if ( present(qnum) ) then
1087 
1088 !OCL XFILL
1089  !$omp parallel do
1090  do j = js, je
1091  do i = is, ie
1092  do k = ks, ke
1093  if ( qnum(k,i,j,i_hc) .ne. undef ) then
1094  qtrc(k,i,j,i_mp_nc) = qnum(k,i,j,i_hc)
1095  else
1096  qtrc(k,i,j,i_mp_nc) = qtrc(k,i,j,i_mp_qc) / ( (piov6*rhow) * dc**b )
1097  end if
1098  end do
1099  end do
1100  end do
1101 
1102 !OCL XFILL
1103  !$omp parallel do
1104  do j = js, je
1105  do i = is, ie
1106  do k = ks, ke
1107  if ( qnum(k,i,j,i_hr) .ne. undef ) then
1108  qtrc(k,i,j,i_mp_nr) = qnum(k,i,j,i_hr)
1109  else
1110  qtrc(k,i,j,i_mp_nr) = qtrc(k,i,j,i_mp_qr) / ( (piov6*rhow) * dr**b )
1111  end if
1112  end do
1113  end do
1114  end do
1115 
1116 !OCL XFILL
1117  !$omp parallel do
1118  do j = js, je
1119  do i = is, ie
1120  do k = ks, ke
1121  if ( qnum(k,i,j,i_hi) .ne. undef ) then
1122  qtrc(k,i,j,i_mp_ni) = qnum(k,i,j,i_hi)
1123  else
1124  qtrc(k,i,j,i_mp_ni) = qtrc(k,i,j,i_mp_qi) / ( (piov6*rhof) * di**b )
1125  end if
1126  end do
1127  end do
1128  end do
1129 
1130 !OCL XFILL
1131  !$omp parallel do
1132  do j = js, je
1133  do i = is, ie
1134  do k = ks, ke
1135  if ( qnum(k,i,j,i_hs) .ne. undef ) then
1136  qtrc(k,i,j,i_mp_ns) = qnum(k,i,j,i_hs)
1137  else
1138  qtrc(k,i,j,i_mp_ns) = qtrc(k,i,j,i_mp_qs) / ( (piov6*rhof) * ds**b )
1139  end if
1140  end do
1141  end do
1142  end do
1143 
1144 !OCL XFILL
1145  !$omp parallel do
1146  do j = js, je
1147  do i = is, ie
1148  do k = ks, ke
1149  if ( qnum(k,i,j,i_hg) .ne. undef ) then
1150  if ( qnum(k,i,j,i_hh) .ne. undef ) then
1151  qtrc(k,i,j,i_mp_ng) = qnum(k,i,j,i_hg) + qnum(k,i,j,i_hh)
1152  else
1153  qtrc(k,i,j,i_mp_ng) = qnum(k,i,j,i_hg)
1154  end if
1155  else
1156  qtrc(k,i,j,i_mp_ng) = qtrc(k,i,j,i_mp_qg) / ( (piov6*rhog) * dg**b )
1157  end if
1158  end do
1159  end do
1160  end do
1161 
1162  else
1163 
1164 !OCL XFILL
1165  !$omp parallel do
1166  do j = js, je
1167  do i = is, ie
1168  do k = ks, ke
1169  qtrc(k,i,j,i_mp_nc) = qtrc(k,i,j,i_mp_qc) / ( (piov6*rhow) * dc**b )
1170  end do
1171  end do
1172  end do
1173 
1174 !OCL XFILL
1175  !$omp parallel do
1176  do j = js, je
1177  do i = is, ie
1178  do k = ks, ke
1179  qtrc(k,i,j,i_mp_nr) = qtrc(k,i,j,i_mp_qr) / ( (piov6*rhow) * dr**b )
1180  end do
1181  end do
1182  end do
1183 
1184 !OCL XFILL
1185  !$omp parallel do
1186  do j = js, je
1187  do i = is, ie
1188  do k = ks, ke
1189  qtrc(k,i,j,i_mp_ni) = qtrc(k,i,j,i_mp_qi) / ( (piov6*rhof) * di**b )
1190  end do
1191  end do
1192  end do
1193 
1194 !OCL XFILL
1195  !$omp parallel do
1196  do j = js, je
1197  do i = is, ie
1198  do k = ks, ke
1199  qtrc(k,i,j,i_mp_ns) = qtrc(k,i,j,i_mp_qs) / ( (piov6*rhof) * ds**b )
1200  end do
1201  end do
1202  end do
1203 
1204 !OCL XFILL
1205  !$omp parallel do
1206  do j = js, je
1207  do i = is, ie
1208  do k = ks, ke
1209  qtrc(k,i,j,i_mp_ng) = qtrc(k,i,j,i_mp_qg) / ( (piov6*rhog) * dg**b )
1210  end do
1211  end do
1212  end do
1213 
1214  end if
1215 
1216  return

References scale_const::const_pi, scale_const::const_undef, scale_atmos_hydrometeor::i_hc, scale_atmos_hydrometeor::i_hg, scale_atmos_hydrometeor::i_hh, scale_atmos_hydrometeor::i_hi, scale_atmos_hydrometeor::i_hr, scale_atmos_hydrometeor::i_hs, and scale_atmos_hydrometeor::n_hyd.

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_qhyd2qtrc().

Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_terminal_velocity()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_terminal_velocity ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
real(rp), dimension(ka), intent(in)  DENS,
real(rp), dimension(ka), intent(in)  TEMP,
real(rp), dimension(ka,i_qc:i_ng), intent(in)  RHOQ,
real(rp), dimension(ka), intent(in)  PRES,
real(rp), dimension(ka,qa_mp-1), intent(out)  vterm 
)

ATMOS_PHY_MP_sn14_terminal_velocity Calculate terminal velocity.

Definition at line 1231 of file scale_atmos_phy_mp_sn14.F90.

1231  use scale_const, only: &
1232  const_undef
1233  implicit none
1234 
1235  integer, intent(in) :: KA, KS, KE
1236 
1237  real(RP), intent(in) :: RHOQ(KA,I_QC:I_NG) ! rho * q
1238  real(RP), intent(in) :: DENS(KA) ! rho
1239  real(RP), intent(in) :: TEMP(KA) ! temperature
1240  real(RP), intent(in) :: PRES(KA) ! pressure
1241 
1242  real(RP), intent(out) :: vterm(KA,QA_MP-1) ! terminal velocity of cloud mass
1243 
1244  real(RP) :: xq ! average mass of 1 particle( mass/number )
1245 
1246  real(RP) :: rhofac ! density factor for terminal velocity( air friction )
1247  real(RP) :: rhofac_q(HYDRO_MAX)
1248 
1249  real(RP) :: rlambdar ! work for diagnosis of Rain DSD ( Seifert, 2008 )
1250  real(RP) :: mud_r
1251  real(RP) :: dq, dql ! weigthed diameter. Improved Rogers etal. (1993) formula by T.Mitsui
1252 
1253 
1254  real(RP) :: weight ! weighting coefficient for 2-branches is determined by ratio between 0.745mm and weighted diameter. SB06 Table.1
1255  real(RP) :: velq_s ! terminal velocity for small branch of Rogers formula
1256  real(RP) :: velq_l ! terminal velocity for large branch of Rogers formula
1257 
1258  real(RP) :: tmp
1259  integer :: k, i, j, iq
1260  !---------------------------------------------------------------------------
1261 
1262  mud_r = 3.0_rp * nu(i_mp_qr) + 2.0_rp
1263 
1264 
1265  do k = ks, ke
1266  rhofac = rho_0 / max( dens(k), rho_min )
1267 
1268  ! QC
1269  rhofac_q(i_mp_qc) = rhofac ** gamma_v(i_mp_qc)
1270  xq = max( xqmin(i_mp_qc), min( xqmax(i_mp_qc), rhoq(k,i_qc) / ( rhoq(k,i_nc) + nqmin(i_mp_qc) ) ) )
1271 
1272  vterm(k,i_mp_qc) = -rhofac_q(i_mp_qc) * coef_vt1(i_mp_qc,1) * xq**beta_v(i_mp_qc,1)
1273  ! NC
1274  vterm(k,i_mp_nc) = -rhofac_q(i_mp_qc) * coef_vt0(i_mp_qc,1) * xq**beta_vn(i_mp_qc,1)
1275 
1276  ! QR
1277  rhofac_q(i_mp_qr) = rhofac ** gamma_v(i_mp_qr)
1278  xq = max( xqmin(i_mp_qr), min( xqmax(i_mp_qr), rhoq(k,i_qr) / ( rhoq(k,i_nr) + nqmin(i_mp_qr) ) ) )
1279 
1280  rlambdar = a_m(i_mp_qr) * xq**b_m(i_mp_qr) &
1281  * ( (mud_r+3.0_rp) * (mud_r+2.0_rp) * (mud_r+1.0_rp) )**(-0.333333333_rp)
1282  dq = ( 4.0_rp + mud_r ) * rlambdar ! D^(3)+mu weighted mean diameter
1283  dql = dq
1284  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )
1285 
1286  velq_s = coef_vtr_ar2 * dq &
1287  * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar )**(-5.0_rp-mud_r) )
1288  velq_l = coef_vtr_ar1 - coef_vtr_br1 &
1289  * ( 1.0_rp + coef_vtr_cr1*rlambdar )**(-4.0_rp-mud_r)
1290  vterm(k,i_mp_qr) = -rhofac_q(i_mp_qr) &
1291  * ( velq_l * ( weight ) &
1292  + velq_s * ( 1.0_rp - weight ) )
1293  ! NR
1294  dq = ( 1.0_rp + mud_r ) * rlambdar
1295  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )
1296 
1297  velq_s = coef_vtr_ar2 * dql &
1298  * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar )**(-2.0_rp-mud_r) )
1299  velq_l = coef_vtr_ar1 - coef_vtr_br1 &
1300  * ( 1.0_rp + coef_vtr_cr1*rlambdar )**(-1.0_rp-mud_r)
1301  vterm(k,i_mp_nr) = -rhofac_q(i_mp_qr) &
1302  * ( velq_l * ( weight ) &
1303  + velq_s * ( 1.0_rp - weight ) )
1304 
1305  ! QI
1306  rhofac_q(i_mp_qi) = ( pres(k)/pre0_vt )**a_pre0_vt * ( temp(k)/tem0_vt )**a_tem0_vt
1307  xq = max( xqmin(i_mp_qi), min( xqmax(i_mp_qi), rhoq(k,i_qi) / ( rhoq(k,i_ni) + nqmin(i_mp_qi) ) ) )
1308 
1309  tmp = a_m(i_mp_qi) * xq**b_m(i_mp_qi)
1310  dq = coef_dave_l(i_mp_qi) * tmp
1311  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_li ) ) ) )
1312 
1313  velq_s = coef_vt1(i_mp_qi,1) * xq**beta_v(i_mp_qi,1)
1314  velq_l = coef_vt1(i_mp_qi,2) * xq**beta_v(i_mp_qi,2)
1315  vterm(k,i_mp_qi) = -rhofac_q(i_mp_qi) &
1316  * ( velq_l * ( weight ) &
1317  + velq_s * ( 1.0_rp - weight ) )
1318  ! NI
1319  dq = coef_dave_n(i_mp_qi) * tmp
1320  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ni ) ) ) )
1321 
1322  velq_s = coef_vt0(i_mp_qi,1) * xq**beta_vn(i_mp_qi,1)
1323  velq_l = coef_vt0(i_mp_qi,2) * xq**beta_vn(i_mp_qi,2)
1324  vterm(k,i_mp_ni) = -rhofac_q(i_mp_qi) &
1325  * ( velq_l * ( weight ) &
1326  + velq_s * ( 1.0_rp - weight ) )
1327 
1328  ! QS
1329  rhofac_q(i_mp_qs) = rhofac_q(i_mp_qi)
1330  xq = max( xqmin(i_mp_qs), min( xqmax(i_mp_qs), rhoq(k,i_qs) / ( rhoq(k,i_ns) + nqmin(i_mp_qs) ) ) )
1331 
1332  tmp = a_m(i_mp_qs) * xq**b_m(i_mp_qs)
1333  dq = coef_dave_l(i_mp_qs) * tmp
1334  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ls ) ) ) )
1335 
1336  velq_s = coef_vt1(i_mp_qs,1) * xq**beta_v(i_mp_qs,1)
1337  velq_l = coef_vt1(i_mp_qs,2) * xq**beta_v(i_mp_qs,2)
1338  vterm(k,i_mp_qs) = -rhofac_q(i_mp_qs) &
1339  * ( velq_l * ( weight ) &
1340  + velq_s * ( 1.0_rp - weight ) )
1341  ! NS
1342  dq = coef_dave_n(i_mp_qs) * tmp
1343  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ns ) ) ) )
1344 
1345  velq_s = coef_vt0(i_mp_qs,1) * xq**beta_vn(i_mp_qs,1)
1346  velq_l = coef_vt0(i_mp_qs,2) * xq**beta_vn(i_mp_qs,2)
1347  vterm(k,i_mp_ns) = -rhofac_q(i_mp_qs) &
1348  * ( velq_l * ( weight ) &
1349  + velq_s * ( 1.0_rp - weight ) )
1350 
1351  ! QG
1352  rhofac_q(i_mp_qg) = rhofac_q(i_mp_qi)
1353  xq = max( xqmin(i_mp_qg), min( xqmax(i_mp_qg), rhoq(k,i_qg) / ( rhoq(k,i_ng) + nqmin(i_mp_qg) ) ) )
1354 
1355  tmp = a_m(i_mp_qg) * xq**b_m(i_mp_qg)
1356  dq = coef_dave_l(i_mp_qg) * tmp
1357  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_lg ) ) ) )
1358 
1359  velq_s = coef_vt1(i_mp_qg,1) * xq**beta_v(i_mp_qg,1)
1360  velq_l = coef_vt1(i_mp_qg,2) * xq**beta_v(i_mp_qg,2)
1361  vterm(k,i_mp_qg) = -rhofac_q(i_mp_qg) &
1362  * ( velq_l * ( weight ) &
1363  + velq_s * ( 1.0_rp - weight ) )
1364  ! NG
1365  dq = coef_dave_n(i_mp_qg) * tmp
1366  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log( dq/d0_ng ) ) ) )
1367 
1368  velq_s = coef_vt0(i_mp_qg,1) * xq**beta_vn(i_mp_qg,1)
1369  velq_l = coef_vt0(i_mp_qg,2) * xq**beta_vn(i_mp_qg,2)
1370  vterm(k,i_mp_ng) = -rhofac_q(i_mp_qg) &
1371  * ( velq_l * ( weight ) &
1372  + velq_s * ( 1.0_rp - weight ) )
1373  enddo
1374 
1375  do iq = 1, qa_mp-1
1376  vterm(1:ks-2,iq) = const_undef
1377  vterm(ks-1,iq) = vterm(ks,iq)
1378  vterm(ke+1:ka,iq) = const_undef
1379  enddo
1380 
1381  return

References aut_acc_slc_brk(), scale_const::const_undef, scale_atmos_hydrometeor::cp_ice, scale_atmos_hydrometeor::cp_vapor, scale_atmos_hydrometeor::cp_water, cross_section(), scale_atmos_hydrometeor::cv_ice, scale_atmos_hydrometeor::cv_vapor, scale_atmos_hydrometeor::cv_water, debug_tem(), freezing_water(), ice_multiplication(), scale_io::io_fid_conf, mixed_phase_collection(), nucleation(), scale_prc::prc_abort(), qa_mp, and scale_specfunc::sf_gamma().

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_calc_tendency().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ debug_tem()

subroutine scale_atmos_phy_mp_sn14::debug_tem ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  point,
integer, intent(in)  i,
integer, intent(in)  j,
real(rp), dimension(ka), intent(in)  tem,
real(rp), dimension(ka), intent(in)  rho,
real(rp), dimension(ka), intent(in)  pre,
real(rp), dimension (ka), intent(in)  qv 
)

Definition at line 3208 of file scale_atmos_phy_mp_sn14.F90.

3208  use scale_prc, only: &
3209  prc_myrank
3210  implicit none
3211  integer, intent(in) :: KA, KS, KE
3212 
3213  integer, intent(in) :: point
3214  integer, intent(in) :: i, j
3215  real(RP), intent(in) :: tem(KA)
3216  real(RP), intent(in) :: rho(KA)
3217  real(RP), intent(in) :: pre(KA)
3218  real(RP), intent(in) :: qv (KA)
3219 
3220  integer :: k
3221  !---------------------------------------------------------------------------
3222 
3223  do k = ks, ke
3224  if ( tem(k) < tem_min &
3225  .OR. rho(k) < rho_min &
3226  .OR. pre(k) < 1.0_rp ) then
3227 
3228  log_info("ATMOS_PHY_MP_SN14_debug_tem_kij",'(A,I3,A,4(F16.5),3(I6))') &
3229  "point: ", point, " low tem,rho,pre:", tem(k), rho(k), pre(k), qv(k), k, i, j, prc_myrank
3230  endif
3231  enddo
3232 
3233  return

References scale_prc::prc_myrank.

Referenced by atmos_phy_mp_sn14_terminal_velocity().

Here is the caller graph for this function:

◆ nucleation()

subroutine scale_atmos_phy_mp_sn14::nucleation ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
real(rp), dimension( ka), intent(in)  cz,
real(rp), dimension(0:ka), intent(in)  fz,
real(rp), dimension (ka), intent(in)  w,
real(rp), dimension (ka), intent(in)  rho,
real(rp), dimension (ka), intent(in)  tem,
real(rp), dimension (ka), intent(in)  pre,
real(rp), dimension(ka), intent(in)  qdry,
real(rp), dimension(ka,i_qv:i_ng), intent(in)  rhoq,
real(rp), dimension(ka), intent(in)  cpa,
real(rp), dimension(ka), intent(in)  dTdt_rad,
real(rp), dimension(ka), intent(in)  qke,
real(rp), dimension(ka), intent(in)  CCN,
real(rp), intent(in)  nc_uplim_d,
real(rp), intent(in)  dt,
real(rp), dimension(ka,pq_max), intent(out)  PQ 
)

Definition at line 3248 of file scale_atmos_phy_mp_sn14.F90.

3248  use scale_prc, only: &
3249  prc_abort
3250  use scale_atmos_saturation, only: &
3251  moist_psat_liq => atmos_saturation_psat_liq, &
3252  moist_psat_ice => atmos_saturation_psat_ice, &
3253  moist_pres2qsat_liq => atmos_saturation_pres2qsat_liq, &
3254  moist_pres2qsat_ice => atmos_saturation_pres2qsat_ice, &
3255  moist_dqsi_dtem_dens => atmos_saturation_dqs_dtem_dens_liq
3256  implicit none
3257 
3258  integer, intent(in) :: KA, KS, KE
3259 
3260  real(RP), intent(in) :: cz( KA)
3261  real(RP), intent(in) :: fz(0:KA)
3262  real(RP), intent(in) :: w (KA) ! w of full level
3263  real(RP), intent(in) :: rho (KA)
3264  real(RP), intent(in) :: tem (KA)
3265  real(RP), intent(in) :: pre (KA)
3266  real(RP), intent(in) :: qdry(KA)
3267  !
3268  real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
3269  !
3270  real(RP), intent(in) :: cpa(KA)
3271  real(RP), intent(in) :: dTdt_rad(KA)
3272  real(RP), intent(in) :: qke(KA)
3273  real(RP), intent(in) :: dt
3274  real(RP), intent(in) :: CCN(KA)
3275  real(RP), intent(in) :: nc_uplim_d
3276  !
3277  real(RP), intent(out) :: PQ(KA,PQ_MAX)
3278  !
3279  !
3280 ! real(RP) :: c_ccn_map(1) ! c_ccn horizontal distribution
3281 ! real(RP) :: kappa_map(1) ! kappa horizontal distribution
3282 ! real(RP) :: c_in_map(1) ! c_in horizontal distribution ! [Add] 11/08/30 T.Mitsui
3283  real(RP) :: esw(KA) ! saturation vapor pressure, water
3284  real(RP) :: esi(KA) ! ice
3285  real(RP) :: ssw(KA) ! super saturation (water)
3286  real(RP) :: ssi(KA) ! super saturation (ice)
3287 ! real(RP) :: w_dsswdz(KA) ! w*(d_ssw/ d_z) super saturation(water) flux
3288  real(RP) :: w_dssidz(KA) ! w*(d_ssi/ d_z), 09/04/14 T.Mitsui
3289 ! real(RP) :: ssw_below(KA) ! ssw(k-1)
3290  real(RP) :: ssi_below(KA) ! ssi(k-1), 09/04/14 T.Mitsui
3291  real(RP) :: z_below(KA) ! z(k-1)
3292  real(RP) :: dzh ! z(k)-z(k-1)
3293  real(RP) :: pv ! vapor pressure
3294  ! work variables for Twomey Equation.
3295  real(RP) :: qsw(KA)
3296  real(RP) :: qsi(KA)
3297  real(RP) :: dqsidtem_rho(KA)
3298  real(RP) :: dssidt_rad(KA)
3299  real(RP) :: wssi, wdssi
3300  !
3301 ! real(RP) :: xi_nuc(1) ! xi use the value @ cloud base
3302 ! real(RP) :: alpha_nuc(1) ! alpha_nuc
3303 ! real(RP) :: eta_nuc(1) ! xi use the value @ cloud base
3304  !
3305  real(RP) :: sigma_w(KA)
3306  real(RP) :: weff(KA)
3307  real(RP) :: weff_max(KA)
3308  real(RP) :: velz(KA)
3309  !
3310  real(RP) :: coef_ccn
3311  real(RP) :: slope_ccn
3312  real(RP) :: nc_new(KA)
3313  real(RP) :: nc_new_below(KA)
3314  real(RP) :: dnc_new
3315  real(RP) :: nc_new_max ! Lohmann (2002),
3316  real(RP) :: a_max
3317  real(RP) :: b_max
3318  logical :: flag_nucleation(KA)
3319  !
3320  real(RP) :: r_gravity
3321  real(RP), parameter :: r_sqrt3=0.577350269_rp ! = sqrt(1.d0/3.d0)
3322  real(RP), parameter :: eps=1.e-30_rp
3323  !====> ! 09/08/18
3324  !
3325  real(RP) :: dlcdt_max, dli_max ! defined by supersaturation
3326  real(RP) :: dncdt_max, dni_max ! defined by supersaturation
3327  real(RP) :: rdt
3328 
3329  real(RP) :: tmp
3330  !
3331  integer :: k
3332  !
3333  !
3334 ! c_ccn_map(1) = c_ccn
3335 ! kappa_map(1) = kappa
3336 ! c_in_map(1) = c_in
3337  !
3338  !
3339  rdt = 1.0_rp/dt
3340  r_gravity = 1.0_rp/grav
3341  !
3342  call moist_psat_liq ( ka, ks, ke, &
3343  tem(:), esw(:) ) ! [IN], [OUT]
3344  call moist_psat_ice ( ka, ks, ke, &
3345  tem(:), esi(:) ) ! [IN], [OUT]
3346  call moist_pres2qsat_liq ( ka, ks, ke, &
3347  tem(:), pre(:), qdry(:), & ! [IN]
3348  qsw(:) ) ! [OUT]
3349  call moist_pres2qsat_ice ( ka, ks, ke, &
3350  tem(:), pre(:), qdry(:), & ! [IN]
3351  qsi(:) ) ! [OUT]
3352  call moist_dqsi_dtem_dens( ka, ks, ke, &
3353  tem(:), rho(:), & ! [IN]
3354  dqsidtem_rho(:) ) ! [OUT]
3355  !
3356  ! Lohmann (2002),JAS, eq.(1) but changing unit [cm-3] => [m-3]
3357  a_max = 1.e+6_rp*0.1_rp*(1.e-6_rp)**1.27_rp
3358  b_max = 1.27_rp
3359  !
3360  ssi_max = 1.0_rp
3361 
3362  do k = ks, ke
3363  pv = rhoq(k,i_qv)*rvap*tem(k)
3364  ssw(k) = min( mp_ssw_lim, ( pv/esw(k)-1.0_rp ) )*100.0_rp
3365  ssi(k) = ( pv/esi(k) - 1.00_rp )
3366 ! ssw_below(k+1) = ssw(k)
3367  ssi_below(k+1) = ssi(k)
3368  z_below(k+1) = cz(k)
3369  end do
3370 ! ssw_below(KS) = ssw(KS)
3371  ssi_below(ks) = ssi(ks)
3372  z_below(ks) = cz(ks-1)
3373 
3374  ! dS/dz is evaluated by first order upstream difference
3375  !*** Solution for Twomey Equation ***
3376 ! coef_ccn = 1.E+6_RP*0.88_RP*(c_ccn_map(1)*1.E-6_RP)**(2.0_RP/(kappa_map(1) + 2.0_RP)) * &
3377  coef_ccn = 1.e+6_rp*0.88_rp*(c_ccn*1.e-6_rp)**(2.0_rp/(kappa + 2.0_rp)) &
3378 ! * (70.0_RP)**(kappa_map(1)/(kappa_map(1) + 2.0_RP))
3379  * (70.0_rp)**(kappa/(kappa + 2.0_rp))
3380 ! slope_ccn = 1.5_RP*kappa_map(1)/(kappa_map(1) + 2.0_RP)
3381  slope_ccn = 1.5_rp*kappa/(kappa + 2.0_rp)
3382  !
3383  do k=ks, ke
3384  sigma_w(k) = r_sqrt3*sqrt(max(qke(k),qke_min))
3385  end do
3386  sigma_w(ks-1) = sigma_w(ks)
3387  sigma_w(ke+1) = sigma_w(ke)
3388  ! effective vertical velocity
3389  do k=ks, ke
3390  weff(k) = w(k) - cpa(k)*r_gravity*dtdt_rad(k)
3391  end do
3392  !
3393  if( mp_couple_aerosol ) then
3394 
3395  do k = ks, ke
3396  if( ssw(k) > 1.e-10_rp .AND. pre(k) > 300.e+2_rp ) then
3397  nc_new(k) = max( ccn(k), c_ccn )
3398  else
3399  nc_new(k) = 0.0_rp
3400  endif
3401  enddo
3402 
3403  else
3404 
3405  if( nucl_twomey ) then
3406  ! diagnose cloud condensation nuclei
3407  do k = ks, ke
3408  ! effective vertical velocity (maximum vertical velocity in turbulent flow)
3409  weff_max(k) = weff(k) + sigma_w(k)
3410  ! large scale upward motion region and saturated
3411  if( (weff(k) > 1.e-8_rp) .AND. (ssw(k) > 1.e-10_rp) .AND. pre(k) > 300.e+2_rp )then
3412  ! Lohmann (2002), eq.(1)
3413  nc_new_max = coef_ccn*weff_max(k)**slope_ccn
3414  nc_new(k) = a_max*nc_new_max**b_max
3415  else
3416  nc_new(k) = 0.0_rp
3417  end if
3418  end do
3419 
3420  else
3421  ! calculate cloud condensation nuclei
3422  do k = ks, ke
3423  if( ssw(k) > 1.e-10_rp .AND. pre(k) > 300.e+2_rp ) then
3424  nc_new(k) = c_ccn*ssw(k)**kappa
3425  else
3426  nc_new(k) = 0.0_rp
3427  endif
3428  enddo
3429  endif
3430 
3431  endif
3432 
3433  do k = ks, ke
3434  ! nc_new is bound by upper limit
3435  if( nc_new(k) > nc_uplim_d ) then ! no more CCN
3436  flag_nucleation(k) = .false.
3437  nc_new_below(k+1) = 1.e+30_rp
3438  else if( nc_new(k) > eps ) then ! nucleation can occur
3439  flag_nucleation(k) = .true.
3440  nc_new_below(k+1) = nc_new(k)
3441  else ! nucleation cannot occur(unsaturated or negative w)
3442  flag_nucleation(k) = .false.
3443  nc_new_below(k+1) = 0.0_rp
3444  end if
3445  end do
3446  nc_new_below(ks) = 0.0_rp
3447 ! do k=KS, KE
3448  ! search maximum value of nc_new
3449 ! if( ( nc_new(k) < nc_new_below(k) ) .OR. &
3450 ! ( nc_new_below(k) > c_ccn_map(1)*0.05_RP ) )then ! 5% of c_ccn
3451 ! ( nc_new_below(k) > c_ccn*0.05_RP ) )then ! 5% of c_ccn
3452 ! flag_nucleation(k) = .false.
3453 ! end if
3454 ! end do
3455 
3456  if( mp_couple_aerosol ) then
3457  do k = ks, ke
3458  ! nucleation occurs at only cloud base.
3459  ! if CCN is more than below parcel, nucleation newly occurs
3460  ! effective vertical velocity
3461  if ( flag_nucleation(k) .AND. & ! large scale upward motion region and saturated
3462  tem(k) > tem_ccn_low ) then
3463  dlcdt_max = ( rhoq(k,i_qv) - esw(k) / ( rvap * tem(k) ) ) * rdt
3464  dlcdt_max = max( dlcdt_max, 0.0_rp ) ! dlcdt_max can be artificially negative due to truncation error in floating point operation
3465  dncdt_max = dlcdt_max/xc_min
3466 ! dnc_new = nc_new(k)-rhoq(k,I_NC)
3467  dnc_new = nc_new(k)
3468  pq(k,i_ncccn) = min( dncdt_max, dnc_new*rdt )
3469  pq(k,i_lcccn) = min( dlcdt_max, xc_min*pq(k,i_ncccn) )
3470  else
3471  pq(k,i_ncccn) = 0.0_rp
3472  pq(k,i_lcccn) = 0.0_rp
3473  end if
3474  end do
3475  else
3476 
3477  if( nucl_twomey ) then
3478  do k = ks, ke
3479  ! nucleation occurs at only cloud base.
3480  ! if CCN is more than below parcel, nucleation newly occurs
3481  ! effective vertical velocity
3482  if ( flag_nucleation(k) .AND. & ! large scale upward motion region and saturated
3483  tem(k) > tem_ccn_low .AND. &
3484  nc_new(k) > rhoq(k,i_nc) ) then
3485  dlcdt_max = ( rhoq(k,i_qv) - esw(k) / ( rvap * tem(k) ) ) * rdt
3486  dlcdt_max = max( dlcdt_max, 0.0_rp ) ! dlcdt_max can be artificially negative due to truncation error in floating point operation
3487  dncdt_max = dlcdt_max/xc_min
3488  dnc_new = nc_new(k)-rhoq(k,i_nc)
3489  pq(k,i_ncccn) = min( dncdt_max, dnc_new*rdt )
3490  pq(k,i_lcccn) = min( dlcdt_max, xc_min*pq(k,i_ncccn) )
3491  else
3492  pq(k,i_ncccn) = 0.0_rp
3493  pq(k,i_lcccn) = 0.0_rp
3494  end if
3495  end do
3496  else
3497  do k = ks, ke
3498  ! effective vertical velocity
3499  if( tem(k) > tem_ccn_low .AND. &
3500  nc_new(k) > rhoq(k,i_nc) ) then
3501  dlcdt_max = ( rhoq(k,i_qv) - esw(k) / ( rvap * tem(k) ) ) * rdt
3502  dlcdt_max = max( dlcdt_max, 0.0_rp ) ! dlcdt_max can be artificially negative due to truncation error in floating point operation
3503  dncdt_max = dlcdt_max/xc_min
3504  dnc_new = nc_new(k)-rhoq(k,i_nc)
3505  pq(k,i_ncccn) = min( dncdt_max, dnc_new*rdt )
3506  pq(k,i_lcccn) = min( dlcdt_max, xc_min*pq(k,i_ncccn) )
3507  else
3508  pq(k,i_ncccn) = 0.0_rp
3509  pq(k,i_lcccn) = 0.0_rp
3510  end if
3511  end do
3512  endif
3513 
3514  endif
3515 
3516  !
3517  ! ice nucleation
3518  !
3519  ! +++ NOTE ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3520  ! Based on Phillips etal.(2006).
3521  ! However this approach doesn't diagnose Ni itself but diagnose tendency.
3522  ! Original approach adjust Ni instantaneously .
3523  ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3524  do k = ks, ke-1
3525  velz(k) = ( w(k) * ( cz(k+1) - fz(k) ) + w(k+1) * ( fz(k) - cz(k) ) ) / ( cz(k+1) - cz(k) ) ! @ half level
3526  end do
3527  velz(ke) = 0.0_rp
3528  do k = ks, ke
3529  dzh = cz(k) - z_below(k)
3530  w_dssidz(k) = velz(k) * (ssi(k) - ssi_below(k))/dzh
3531  dssidt_rad(k) = -rhoq(k,i_qv)/(rho(k)*qsi(k)*qsi(k))*dqsidtem_rho(k)*dtdt_rad(k)
3532  dli_max = ( rhoq(k,i_qv) - esi(k) / ( rvap * tem(k) ) ) * rdt
3533  dli_max = max( dli_max, 0.0_rp ) ! dli_max can be artificially negative due to truncation error in floating point operation
3534  dni_max = min( dli_max/xi_ccn, (in_max-rhoq(k,i_ni))*rdt )
3535  wdssi = min( w_dssidz(k)+dssidt_rad(k), 0.01_rp)
3536  wssi = min( ssi(k), ssi_max)
3537  ! SB06(34),(35)
3538  if( ( wdssi > eps ) .AND. & !
3539  (tem(k) < 273.15_rp ) .AND. & !
3540  (rhoq(k,i_ni) < in_max ) .AND. &
3541  (wssi >= eps ) )then !
3542  tmp = c_in * nm_m92 * exp( 0.3_rp * bm_m92 * ( wssi - 0.1_rp ) )
3543  if( inucl_w ) then
3544  tmp = bm_m92 * 0.3_rp * tmp * wdssi
3545  else
3546  tmp = max( tmp - rhoq(k,i_ni), 0.0_rp ) * rdt
3547  endif
3548  pq(k,i_niccn) = min(dni_max, tmp)
3549  pq(k,i_liccn) = min(dli_max, pq(k,i_niccn)*xi_ccn )
3550  else
3551  pq(k,i_niccn) = 0.0_rp
3552  pq(k,i_liccn) = 0.0_rp
3553  end if
3554  end do
3555 
3556  return

References scale_prc::prc_abort().

Referenced by atmos_phy_mp_sn14_terminal_velocity().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ice_multiplication()

subroutine scale_atmos_phy_mp_sn14::ice_multiplication ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
logical, intent(in)  flg_lt,
real(rp), dimension(ka,pac_max), intent(in)  Pac,
real(rp), dimension(ka), intent(in)  tem,
real(rp), dimension(ka,i_qv:i_ng), intent(in)  rhoq,
real(rp), dimension(i_qc:i_qg,ka), intent(in)  rhoq_crg,
real(rp), dimension(ka,hydro_max), intent(in)  xq,
real(rp), dimension(ka,pq_max), intent(inout)  PQ,
real(rp), dimension(pq_max,ka), intent(inout)  Pcrg1 
)

Definition at line 3567 of file scale_atmos_phy_mp_sn14.F90.

3567 
3568  ! ice multiplication by splintering
3569  ! we consider Hallet-Mossop process
3570  use scale_specfunc, only: &
3571  gammafunc => sf_gamma
3572  implicit none
3573 
3574  integer, intent(in) :: KA, KS, KE
3575  !
3576  real(RP), intent(in) :: Pac(KA,Pac_MAX)
3577  real(RP), intent(in) :: tem(KA)
3578  real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
3579  real(RP), intent(in) :: xq(KA,HYDRO_MAX)
3580  !
3581  real(RP), intent(inout):: PQ(KA,PQ_MAX)
3582  ! for lightning
3583  logical, intent(in) :: flg_lt
3584  real(RP), intent(in) :: rhoq_crg(I_QC:I_QG,KA)
3585  real(RP), intent(inout):: Pcrg1(PQ_MAX,KA)
3586  !
3587  ! production of (350.d3 ice particles)/(cloud riming [g]) => 350*d6 [/kg]
3588  real(RP), parameter :: pice = 350.0e+6_rp
3589  ! production of (1 ice particle)/(250 cloud particles riming)
3590  real(RP), parameter :: pnc = 250.0_rp
3591  ! temperature factor
3592  real(RP) :: fp
3593 
3594  real(RP) :: igm ! in complete gamma(x,alpha)
3595  real(RP) :: x
3596  ! coefficient of expansion using in calculation of igm
3597  real(RP) :: a0,a1,a2,a3,a4,a5
3598  real(RP) :: a6,a7,a8,a9,a10
3599  real(RP) :: an1,an2,b0,b1,b2,c0,c1,c2
3600  real(RP) :: d0,d1,d2,e1,e2,h0,h1,h2
3601  real(RP), parameter :: eps=1.0e-30_rp
3602  ! number of cloud droplets larger than 12 micron(radius).
3603  real(RP) :: n12
3604  !
3605  real(RP) :: wn, wni, wns, wng
3606  integer :: k, iq
3607  real(RP) :: sw1
3608  !
3609  !
3610  do k = ks, ke
3611  ! Here we assume particle temperature is same as environment temperature.
3612  ! If you want to treat in a better manner,
3613  ! you can diagnose with eq.(64) in CT(86)
3614  if (tem(k) > 270.16_rp)then
3615  fp = 0.0_rp
3616  else if(tem(k) >= 268.16_rp)then
3617  fp = (270.16_rp-tem(k))*0.5_rp
3618  else if(tem(k) >= 265.16_rp)then
3619  fp = (tem(k)-265.16_rp)*0.333333333_rp
3620  else
3621  fp = 0.0_rp
3622  end if
3623  ! Approximation of Incomplete Gamma function
3624  ! Here we calculate with algorithm by Numerical Recipes.
3625  ! This approach is based on eq.(78) in Cotton etal.(1986),
3626  ! but more accurate and expanded for Generalized Gamma Distribution.
3627  x = coef_lambda(i_mp_qc)*(xc_cr/xq(k,i_mp_qc))**mu(i_mp_qc)
3628  !
3629  if(x<1.e-2_rp*alpha)then ! negligible
3630  igm = 0.0_rp
3631  else if(x<alpha+1.0_rp)then ! series expansion
3632  ! 10th-truncation is enough for cloud droplet.
3633  a0 = 1.0_rp/alpha ! n=0
3634  a1 = a0*x/(alpha+1.0_rp) ! n=1
3635  a2 = a1*x/(alpha+2.0_rp) ! n=2
3636  a3 = a2*x/(alpha+3.0_rp) ! n=3
3637  a4 = a3*x/(alpha+4.0_rp) ! n=4
3638  a5 = a4*x/(alpha+5.0_rp) ! n=5
3639  a6 = a5*x/(alpha+6.0_rp) ! n=6
3640  a7 = a6*x/(alpha+7.0_rp) ! n=7
3641  a8 = a7*x/(alpha+8.0_rp) ! n=8
3642  a9 = a8*x/(alpha+9.0_rp) ! n=9
3643  a10 = a9*x/(alpha+10.0_rp) ! n=10
3644  igm = (a0+a1+a2+a3+a4+a5+a6+a7+a8+a9+a10)*exp( -x + alpha*log(x) - lgm )
3645  else if(x<alpha*1.d2) then ! continued fraction expansion
3646  ! 2nd-truncation is enough for cloud droplet.
3647  ! setup
3648  b0 = x+1.0_rp-alpha
3649  c0 = 1.0_rp/eps
3650  d0 = 1.0_rp/b0
3651  h0 = d0
3652  ! n=1
3653  an1 = -(1.0_rp-alpha)
3654  b1 = b0 + 2.0_rp
3655  d1 = 1.0_rp/(an1*d0+b1)
3656  c1 = b1+an1/c0
3657  e1 = d1*c1
3658  h1 = h0*e1
3659  ! n=2
3660  an2 = -2.0_rp*(2.0_rp-alpha)
3661  b2 = b1 + 2.0_rp
3662  d2 = 1.0_rp/(an2*d1+b2)
3663  c2 = b2+an2/c1
3664  e2 = d2*c2
3665  h2 = h1*e2
3666  !
3667  igm = 1.0_rp - exp( -x + alpha*log(x) - lgm )*h2
3668  else ! negligible
3669  igm = 1.0_rp
3670  end if
3671  ! n12 is number of cloud droplets larger than 12 micron.
3672  n12 = rhoq(k,i_nc)*(1.0_rp-igm)
3673  ! eq.(82) CT(86)
3674  wn = (pice + n12/((rhoq(k,i_qc)+xc_min)*pnc) )*fp ! filtered by xc_min
3675  wni = wn*(-pac(k,i_liaclc2li) ) ! riming production rate is all negative
3676  wns = wn*(-pac(k,i_lsaclc2ls) )
3677  wng = wn*(-pac(k,i_lgaclc2lg) )
3678  pq(k,i_nispl) = wni+wns+wng
3679  !
3680  pq(k,i_lsspl) = - wns*xq(k,i_mp_qi) ! snow => ice
3681  pq(k,i_lgspl) = - wng*xq(k,i_mp_qi) ! graupel => ice
3682  if (flg_lt) then
3683  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NC is small,ignore charge transfer
3684  pcrg1(i_nsspl,k) = - wns*(1.0_rp-sw1) &
3685  / (rhoq(k,i_ns)+sw1)*rhoq_crg(i_qs,k)
3686  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small ) !--- if NC is small,ignore charge transfer
3687  pcrg1(i_ngspl,k) = - wng*(1.0_rp-sw1) &
3688  / (rhoq(k,i_ng)+sw1)*rhoq_crg(i_qg,k)
3689  pcrg1(i_nispl,k) = - ( pcrg1(i_nsspl,k) + pcrg1(i_ngspl,k) )
3690  end if
3691  !
3692  end do
3693  !
3694  return

References scale_specfunc::sf_gamma().

Referenced by atmos_phy_mp_sn14_terminal_velocity().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ mixed_phase_collection()

subroutine scale_atmos_phy_mp_sn14::mixed_phase_collection ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
logical, intent(in)  flg_lt,
real(rp), intent(in)  d0_crg,
real(rp), intent(in)  v0_crg,
real(rp), dimension(ka), intent(in)  beta_crg,
real(rp), dimension(ka), intent(in)  dqcrg,
real(rp), dimension(ka), intent(in)  wtem,
real(rp), dimension(ka,i_qv:i_ng), intent(in)  rhoq,
real(rp), dimension(i_qc:i_qg,ka), intent(in)  rhoq_crg,
real(rp), dimension(ka,hydro_max), intent(in)  xq,
real(rp), dimension(ka,hydro_max), intent(in)  dq_xave,
real(rp), dimension(ka,hydro_max,2), intent(in)  vt_xave,
real(rp), dimension(ka,pq_max), intent(inout)  PQ,
real(rp), dimension(pq_max,ka), intent(inout)  Pcrg1,
real(rp), dimension(pcrg_max,ka), intent(inout)  Pcrg2,
real(rp), dimension(ka,pac_max), intent(out)  Pac 
)

Definition at line 3710 of file scale_atmos_phy_mp_sn14.F90.

3710  use scale_atmos_saturation, only: &
3711  moist_psat_ice => atmos_saturation_psat_ice
3712  implicit none
3713 
3714  integer, intent(in) :: KA, KS, KE
3715 
3716  !--- mixed-phase collection process
3717  ! And all we set all production term as a negative sign to avoid confusion.
3718  !
3719  real(RP), intent(in) :: wtem(KA)
3720  !--- mass/number concentration[kg/m3]
3721  real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
3722  ! necessary ?
3723  real(RP), intent(in) :: xq(KA,HYDRO_MAX)
3724  !--- diameter of averaged mass( D(ave x) )
3725  real(RP), intent(in) :: dq_xave(KA,HYDRO_MAX)
3726  !--- terminal velocity of averaged mass( vt(ave x) )
3727  real(RP), intent(in) :: vt_xave(KA,HYDRO_MAX,2)
3728  ! [Add] 11/08/30 T.Mitsui, for autoconversion of ice
3729  ! real(RP), intent(in) :: rho(KA)
3730  !--- partial conversion
3731  real(RP), intent(inout):: PQ(KA,PQ_MAX)
3732  !
3733  real(RP), intent(out):: Pac(KA,Pac_MAX)
3734  !--- for lightning component
3735  logical, intent(in) :: flg_lt
3736  real(RP), intent(in) :: beta_crg(KA)
3737  real(RP), intent(in) :: dqcrg(KA)
3738  real(RP), intent(in) :: d0_crg, v0_crg
3739  real(RP), intent(in) :: rhoq_crg(I_QC:I_QG,KA)
3740  real(RP), intent(inout):: Pcrg1(PQ_MAX,KA)
3741  real(RP), intent(inout):: Pcrg2(Pcrg_MAX,KA)
3742 
3743  real(RP), parameter :: a_dec = 0.883_rp
3744  real(RP), parameter :: b_dec = 0.093_rp
3745  real(RP), parameter :: c_dec = 0.00348_rp
3746  real(RP), parameter :: d_dec = 4.5185e-5_rp
3747  !
3748  !
3749  real(RP) :: tem(KA)
3750  !
3751  !--- collection efficency of each specie
3752  real(RP) :: E_c, E_r, E_i, E_s, E_g !
3753  real(RP) :: E_ic, E_sc, E_gc !
3754  !--- sticking efficiency
3755  real(RP) :: E_stick(KA)
3756  ! [Add] 10/08/03 T.Mitsui
3757  real(RP) :: temc, temc2, temc3
3758  real(RP) :: E_dec
3759  real(RP) :: esi_rat
3760  real(RP) :: esi(KA)
3761  !
3762  real(RP) :: temc_p, temc_m ! celcius tem.
3763 ! real(RP) :: ci_aut(KA)
3764 ! real(RP) :: taui_aut(KA)
3765 ! real(RP) :: tau_sce(KA)
3766  !--- DSD averaged diameter for each species
3767  real(RP) :: ave_dc ! cloud
3768 ! real(RP) :: ave_dr ! rain
3769  real(RP) :: ave_di ! ice
3770  real(RP) :: ave_ds ! snow
3771  real(RP) :: ave_dg ! graupel
3772  !--- coefficient of collection equations(L:mass, N:number)
3773  real(RP) :: coef_acc_LCI, coef_acc_NCI ! cloud - cloud ice
3774  real(RP) :: coef_acc_LCS, coef_acc_NCS ! cloud - snow
3775  !
3776  real(RP) :: coef_acc_LCG, coef_acc_NCG ! cloud - graupel
3777  real(RP) :: coef_acc_LRI_I, coef_acc_NRI_I ! rain - cloud ice
3778  real(RP) :: coef_acc_LRI_R, coef_acc_NRI_R ! rain - cloud ice
3779  real(RP) :: coef_acc_LRS_S, coef_acc_NRS_S ! rain - snow
3780  real(RP) :: coef_acc_LRS_R, coef_acc_NRS_R ! rain - snow
3781  real(RP) :: coef_acc_LRG, coef_acc_NRG ! rain - graupel
3782  real(RP) :: coef_acc_LII, coef_acc_NII ! cloud ice - cloud ice
3783  real(RP) :: coef_acc_LIS, coef_acc_NIS ! cloud ice - snow
3784  real(RP) :: coef_acc_NSS ! snow - snow
3785  real(RP) :: coef_acc_NGG ! grauepl - graupel
3786  real(RP) :: coef_acc_LSG, coef_acc_NSG ! snow - graupel
3787  !--- (diameter) x (diameter)
3788  real(RP) :: dcdc, dcdi, dcds, dcdg
3789  real(RP) :: drdr, drdi, drds, drdg
3790  real(RP) :: didi, dids, didg
3791  real(RP) :: dsds, dsdg
3792  real(RP) :: dgdg
3793  !--- (terminal velocity) x (terminal velocity)
3794  real(RP) :: vcvc, vcvi, vcvs, vcvg
3795  real(RP) :: vrvr, vrvi, vrvs, vrvg
3796  real(RP) :: vivi, vivs, vivg
3797  real(RP) :: vsvs, vsvg
3798  real(RP) :: vgvg
3799  !
3800  real(RP) :: wx_cri, wx_crs
3801  real(RP) :: coef_emelt
3802  real(RP) :: w1
3803 
3804  real(RP) :: sw, sw1, sw2
3805  real(RP) :: alpha_lt
3806  !
3807  integer :: k, iqw
3808  !
3809  !
3810  do k = ks, ke
3811  tem(k) = max( wtem(k), tem_min ) ! 11/08/30 T.Mitsui
3812  end do
3813 
3814  call moist_psat_ice( ka, ks, ke, &
3815  tem(:), esi(:) ) ! [IN], [OUT]
3816 
3817  if( opt_stick_ks96 )then
3818  do k = ks, ke
3819  ! Khain and Sednev (1996), eq.(3.15)
3820  temc = tem(k) - t00
3821  temc2 = temc*temc
3822  temc3 = temc2*temc
3823  e_dec = max(0.0_rp, a_dec + b_dec*temc + c_dec*temc2 + d_dec*temc3 )
3824  esi_rat = rhoq(k,i_qv)*rvap*tem(k)/esi(k)
3825  e_stick(k) = min(1.0_rp, e_dec*esi_rat)
3826  end do
3827  else if( opt_stick_co86 )then
3828  do k = ks, ke
3829  ! [Add] 11/08/30 T.Mitsui, Cotton et al. (1986)
3830  temc = min(tem(k) - t00,0.0_rp)
3831  w1 = 0.035_rp*temc-0.7_rp
3832  e_stick(k) = 10._rp**w1
3833  end do
3834  else
3835  do k = ks, ke
3836  ! Lin et al. (1983)
3837  temc_m = min(tem(k) - t00,0.0_rp) ! T < 273.15
3838  e_stick(k) = exp(0.09_rp*temc_m)
3839  end do
3840  end if
3841  !
3842 !OCL NOSIMD
3843  do k = ks, ke
3844  !
3845 ! temc_m = min(tem(k) - T00,0.0_RP) ! T < 273.15
3846  temc_p = max(tem(k) - t00,0.0_rp) ! T > 273.15
3847  ! averaged diameter using SB06(82)
3848  ave_dc = coef_d(i_mp_qc)*xq(k,i_mp_qc)**b_m(i_mp_qc)
3849  ave_di = coef_d(i_mp_qi)*xq(k,i_mp_qi)**b_m(i_mp_qi)
3850  ave_ds = coef_d(i_mp_qs)*xq(k,i_mp_qs)**b_m(i_mp_qs)
3851  ave_dg = coef_d(i_mp_qg)*xq(k,i_mp_qg)**b_m(i_mp_qg)
3852  !------------------------------------------------------------------------
3853  ! coellection efficiency are given as follows
3854  e_c = max(0.0_rp, min(1.0_rp, (ave_dc-dc0)/(dc1-dc0) ))
3855  sw = 0.5_rp - sign(0.5_rp, di0-ave_di) ! if(ave_di>di0)then sw=1
3856  e_i = e_im * sw
3857  sw = 0.5_rp - sign(0.5_rp, ds0-ave_ds) ! if(ave_ds>ds0)then sw=1
3858  e_s = e_sm * sw
3859  sw = 0.5_rp - sign(0.5_rp, dg0-ave_dg) ! if(ave_dg>dg0)then sw=1
3860  e_g = e_gm * sw
3861  e_ic = e_i*e_c
3862  e_sc = e_s*e_c
3863  e_gc = e_g*e_c
3864  !------------------------------------------------------------------------
3865  ! Collection: a collects b ( assuming particle size a>b )
3866  dcdc = dq_xave(k,i_mp_qc) * dq_xave(k,i_mp_qc)
3867  drdr = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qr)
3868  didi = dq_xave(k,i_mp_qi) * dq_xave(k,i_mp_qi)
3869  dsds = dq_xave(k,i_mp_qs) * dq_xave(k,i_mp_qs)
3870  dgdg = dq_xave(k,i_mp_qg) * dq_xave(k,i_mp_qg)
3871  dcdi = dq_xave(k,i_mp_qc) * dq_xave(k,i_mp_qi)
3872  dcds = dq_xave(k,i_mp_qc) * dq_xave(k,i_mp_qs)
3873  dcdg = dq_xave(k,i_mp_qc) * dq_xave(k,i_mp_qg)
3874  drdi = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qi)
3875  drds = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qs)
3876  drdg = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qg)
3877  dids = dq_xave(k,i_mp_qi) * dq_xave(k,i_mp_qs)
3878 ! didg = dq_xave(k,I_mp_QI) * dq_xave(k,I_mp_QG)
3879  dsdg = dq_xave(k,i_mp_qs) * dq_xave(k,i_mp_qg)
3880  !
3881  vcvc = vt_xave(k,i_mp_qc,2)* vt_xave(k,i_mp_qc,2)
3882  vrvr = vt_xave(k,i_mp_qr,2)* vt_xave(k,i_mp_qr,2)
3883  vivi = vt_xave(k,i_mp_qi,2)* vt_xave(k,i_mp_qi,2)
3884  vsvs = vt_xave(k,i_mp_qs,2)* vt_xave(k,i_mp_qs,2)
3885  vgvg = vt_xave(k,i_mp_qg,2)* vt_xave(k,i_mp_qg,2)
3886  vcvi = vt_xave(k,i_mp_qc,2)* vt_xave(k,i_mp_qi,2)
3887  vcvs = vt_xave(k,i_mp_qc,2)* vt_xave(k,i_mp_qs,2)
3888  vcvg = vt_xave(k,i_mp_qc,2)* vt_xave(k,i_mp_qg,2)
3889  vrvi = vt_xave(k,i_mp_qr,2)* vt_xave(k,i_mp_qi,2)
3890  vrvs = vt_xave(k,i_mp_qr,2)* vt_xave(k,i_mp_qs,2)
3891  vrvg = vt_xave(k,i_mp_qr,2)* vt_xave(k,i_mp_qg,2)
3892  vivs = vt_xave(k,i_mp_qi,2)* vt_xave(k,i_mp_qs,2)
3893 ! vivg = vt_xave(k,I_mp_QI,2)* vt_xave(k,I_mp_QG,2)
3894  vsvg = vt_xave(k,i_mp_qs,2)* vt_xave(k,i_mp_qg,2)
3895  !------------------------------------------------------------------------
3896  !
3897  !+++ pattern 1: a + b => a (a>b)
3898  ! (i-c, s-c, g-c, s-i, g-r, s-g)
3899  !------------------------------------------------------------------------
3900  ! cloud-ice => ice
3901  ! reduction term of cloud
3902  coef_acc_lci = &
3903  ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi ) &
3904  * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi &
3905  + sigma_i + sigma_c )**0.5_rp
3906  coef_acc_nci = &
3907  ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi ) &
3908  * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi &
3909  + sigma_i + sigma_c )**0.5_rp
3910  pac(k,i_liaclc2li)= -0.25_rp*pi*e_ic*rhoq(k,i_ni)*rhoq(k,i_qc)*coef_acc_lci
3911  pac(k,i_niacnc2ni)= -0.25_rp*pi*e_ic*rhoq(k,i_ni)*rhoq(k,i_nc)*coef_acc_nci
3912  !---- for charge density
3913  if (flg_lt) then !--- C + I -> I (decrease from cloud chgarge)
3914  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small ) !--- if NC is small, ignore charge transfer
3915  pcrg2(i_niacnc2ni,k) = pac(k,i_niacnc2ni)*(1.0_rp-sw1) &
3916  / (rhoq(k,i_nc)+sw1)*rhoq_crg(i_qc,k)
3917  end if
3918  ! cloud-snow => snow
3919  ! reduction term of cloud
3920  coef_acc_lcs = &
3921  ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds ) &
3922  * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs &
3923  + sigma_s + sigma_c )**0.5_rp
3924  coef_acc_ncs = &
3925  ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds ) &
3926  * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs &
3927  + sigma_s + sigma_c )**0.5_rp
3928  pac(k,i_lsaclc2ls)= -0.25_rp*pi*e_sc*rhoq(k,i_ns)*rhoq(k,i_qc)*coef_acc_lcs
3929  pac(k,i_nsacnc2ns)= -0.25_rp*pi*e_sc*rhoq(k,i_ns)*rhoq(k,i_nc)*coef_acc_ncs
3930  !---- for charge density
3931  if (flg_lt) then !--- C + S -> S (decrease from cloud charge)
3932  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small ) !--- if NC is small, ignore charge transfer
3933  pcrg2(i_nsacnc2ns,k) = pac(k,i_nsacnc2ns)*(1.0_rp-sw1) &
3934  / (rhoq(k,i_nc)+sw1)*rhoq_crg(i_qc,k)
3935  end if
3936  ! cloud-graupel => graupel
3937  ! reduction term of cloud
3938  coef_acc_lcg = &
3939  ( delta_b1(i_mp_qc)*dcdc + delta_ab1(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg ) &
3940  * ( theta_b1(i_mp_qc)*vcvc - theta_ab1(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg &
3941  + sigma_g + sigma_c )**0.5_rp
3942  coef_acc_ncg = &
3943  ( delta_b0(i_mp_qc)*dcdc + delta_ab0(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg ) &
3944  * ( theta_b0(i_mp_qc)*vcvc - theta_ab0(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg &
3945  + sigma_g + sigma_c )**0.5_rp
3946  pac(k,i_lgaclc2lg)= -0.25_rp*pi*e_gc*rhoq(k,i_ng)*rhoq(k,i_qc)*coef_acc_lcg
3947  pac(k,i_ngacnc2ng)= -0.25_rp*pi*e_gc*rhoq(k,i_ng)*rhoq(k,i_nc)*coef_acc_ncg
3948  !---- for charge density
3949  if (flg_lt) then !--- C + G -> G (decrease from cloud charge)
3950  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small ) !--- if NC is small, ignore charge transfer
3951  pcrg2(i_ngacnc2ng,k) = pac(k,i_ngacnc2ng)*(1.0_rp-sw1) &
3952  / (rhoq(k,i_nc)+sw1)*rhoq_crg(i_qc,k)
3953  end if
3954  ! snow-graupel => graupel
3955  coef_acc_lsg = &
3956  ( delta_b1(i_mp_qs)*dsds + delta_ab1(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg ) &
3957  * ( theta_b1(i_mp_qs)*vsvs - theta_ab1(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg &
3958  + sigma_g + sigma_s )**0.5_rp
3959  coef_acc_nsg = &
3960  ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg ) &
3961  ! [fix] T.Mitsui 08/05/08
3962  * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg &
3963  + sigma_g + sigma_s )**0.5_rp
3964  pac(k,i_lgacls2lg)= -0.25_rp*pi*e_stick(k)*e_gs*rhoq(k,i_ng)*rhoq(k,i_qs)*coef_acc_lsg
3965  pac(k,i_ngacns2ng)= -0.25_rp*pi*e_stick(k)*e_gs*rhoq(k,i_ng)*rhoq(k,i_ns)*coef_acc_nsg
3966  !---- for charge density
3967  if (flg_lt) then !--- S + G -> G (decrease from snow charge)
3968  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NS is small, ignore charge transfer
3969  pcrg2(i_ngacns2ng,k) = pac(k,i_ngacns2ng)*(1.0_rp-sw1) &
3970  / (rhoq(k,i_ns)+sw1)*rhoq_crg(i_qs,k)
3971  !--- Charge split by Snow-Graupel rebound--------------------------------
3972  alpha_lt = 5.0_rp * ( dq_xave(k,i_mp_qs) / d0_crg )**2 * vt_xave(k,i_mp_qg,2) / v0_crg
3973  alpha_lt = min( alpha_lt, 10.0_rp )
3974  pcrg2(i_cgngacns2ng,k)= 0.25_rp*pi*( 1.0_rp - e_stick(k) )*e_gs &
3975  * rhoq(k,i_ng)*rhoq(k,i_ns)*coef_acc_nsg &
3976  * ( dqcrg(k)*alpha_lt ) &
3977  * beta_crg(k)
3978  end if
3979 !! !-----------------
3980 !! ! (start) Y.Sato added on 2018/8/31
3981 !! !--- ice-graupel => graupel
3982 !! coef_acc_LIG = &
3983 !! ( delta_b1(I_QI)*didi + delta_ab1(I_QG,I_QI)*didg + delta_b0(I_QG)*dgdg ) &
3984 !! * ( theta_b1(I_QI)*vivi - theta_ab1(I_QG,I_QI)*vivg + theta_b0(I_QG)*vgvg &
3985 !! + sigma_g + sigma_i )**0.5_RP
3986 !! coef_acc_NIG = &
3987 !! ( delta_b0(I_QI)*didi + delta_ab0(I_QG,I_QI)*didg + delta_b0(I_QG)*dgdg ) &
3988 !! ! [fix] T.Mitsui 08/05/08
3989 !! * ( theta_b0(I_QI)*vivi - theta_ab0(I_QG,I_QI)*vivg + theta_b0(I_QG)*vgvg &
3990 !! + sigma_g + sigma_i )**0.5_RP
3991 !! Pac(k,I_LGacLI2LG)= -0.25_RP*pi*E_stick(k)*E_gi*rhoq(k,I_NG)*rhoq(k,I_QI)*coef_acc_LIG*flg_igcol
3992 !! Pac(k,I_NGacNI2NG)= -0.25_RP*pi*E_stick(k)*E_gi*rhoq(k,I_NG)*rhoq(k,I_NI)*coef_acc_NIG*flg_igcol
3993 !! !---- for charge density
3994 !! if (flg_lt) then !--- I + G -> G (decrease from snow charge)
3995 !! sw1 = 0.5_RP - sign( 0.5_RP, rhoq(k,I_NI)-SMALL ) !--- if NS is small, ignore charge transfer
3996 !! Pcrg2(I_NGacNI2NG,k) = Pac(k,I_NGacNI2NG)*(1.0_RP-sw1) &
3997 !! / (rhoq(k,I_NI)+sw1)*rhoq_crg(I_QI,k)*flg_lt*flg_igcol
3998 !! !--- Charge split by Ice-Graupel rebound--------------------------------
3999 !! alpha_lt = 5.0_RP * ( dq_xave(k,I_mp_QI) / d0_crg )**2 * vt_xave(k,I_mp_QG,2) / v0_crg
4000 !! alpha_lt = min( alpha_lt, 10.0_RP )
4001 !! Pcrg2(I_CGNGacNI2NG,k)= 0.25_RP*pi*( 1.0_RP - E_stick(k) )*E_gi &
4002 !! * rhoq(k,I_NG)*rhoq(k,I_NI)*coef_acc_NIG &
4003 !! * ( dqcrg(k)*alpha_lt ) &
4004 !! * beta_crg(k) * flg_lt * flg_igcol
4005 !! end if
4006 !! ! (end) Y.Sato added on 2018/8/31
4007 !! !------------------
4008  !------------------------------------------------------------------------
4009  ! ice-snow => snow
4010  ! reduction term of ice
4011  coef_acc_lis = &
4012  ( delta_b1(i_mp_qi)*didi + delta_ab1(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds ) &
4013  * ( theta_b1(i_mp_qi)*vivi - theta_ab1(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs &
4014  + sigma_i + sigma_s )**0.5_rp
4015  coef_acc_nis = &
4016  ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds ) &
4017  * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs &
4018  + sigma_i + sigma_s )**0.5_rp
4019  pac(k,i_liacls2ls)= -0.25_rp*pi*e_stick(k)*e_si*rhoq(k,i_ns)*rhoq(k,i_qi)*coef_acc_lis
4020  pac(k,i_niacns2ns)= -0.25_rp*pi*e_stick(k)*e_si*rhoq(k,i_ns)*rhoq(k,i_ni)*coef_acc_nis
4021  !
4022  !---- for charge density
4023  if (flg_lt) then !--- I + S -> S (decrease from ice charge)
4024  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small, ignore charge transfer
4025  pcrg2(i_niacns2ns,k) = pac(k,i_niacns2ns)*(1.0_rp-sw1) &
4026  / (rhoq(k,i_ni)+sw1)*rhoq_crg(i_qi,k)
4027  end if
4028  sw = sign(0.5_rp, t00-tem(k)) + 0.5_rp
4029  ! if ( tem(k) <= T00 )then
4030  ! rain-graupel => graupel
4031  ! reduction term of rain
4032  ! sw = 1
4033  ! else
4034  ! rain-graupel => rain
4035  ! reduction term of graupel
4036  ! sw = 0
4037  coef_acc_lrg = &
4038  ( ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg ) * sw &
4039  + ( 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) ) &
4040  * sqrt( ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qg,i_mp_qr)*vrvg + theta_b0(i_mp_qg)*vgvg ) * sw &
4041  + ( 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) &
4042  + sigma_r + sigma_g )
4043  pac(k,i_lraclg2lg) = -0.25_rp*pi*e_gr*coef_acc_lrg &
4044  * ( rhoq(k,i_ng)*rhoq(k,i_qr) * sw &
4045  + rhoq(k,i_nr)*rhoq(k,i_qg) * (1.0_rp-sw) )
4046  coef_acc_nrg = &
4047  ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg ) &
4048  * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qg,i_mp_qr)*vrvg + theta_b0(i_mp_qg)*vgvg &
4049  + sigma_r + sigma_g )**0.5_rp
4050  pac(k,i_nracng2ng) = -0.25_rp*pi*e_gr*rhoq(k,i_ng)*rhoq(k,i_nr)*coef_acc_nrg
4051  !---- for charge density
4052  if (flg_lt) then !--- R+G->R (T>T00 sw=0, decrease from graupel charge), ->G(T<=T00 sw=1, dcrerase from rain charge)
4053  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small ) !--- if NR is small, ignore charge transfer
4054  sw2 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small ) !--- if NG is small, ignore charge transfer
4055  pcrg2(i_nracng2ng,k) = pac(k,i_nracng2ng)*(1.0_rp-sw1)/(rhoq(k,i_nr)+sw1) &
4056  *rhoq_crg(i_qr,k) * sw &
4057  + pac(k,i_nracng2ng)*(1.0_rp-sw2)/(rhoq(k,i_ng)+sw2) &
4058  *rhoq_crg(i_qg,k) * (1.0_rp-sw)
4059  end if
4060  !
4061  !------------------------------------------------------------------------
4062  !
4063  !+++ pattern 2: a + b => c (a>b)
4064  ! (r-i,r-s)
4065  !------------------------------------------------------------------------
4066  ! rain-ice => graupel
4067  ! reduction term of ice
4068  coef_acc_lri_i = &
4069  ( delta_b1(i_mp_qi)*didi + delta_ab1(i_mp_qr,i_mp_qi)*drdi + delta_b0(i_mp_qr)*drdr ) &
4070  * ( theta_b1(i_mp_qi)*vivi - theta_ab1(i_mp_qr,i_mp_qi)*vrvi + theta_b0(i_mp_qr)*vrvr &
4071  + sigma_r + sigma_i )**0.5_rp
4072  coef_acc_nri_i = &
4073  ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qr,i_mp_qi)*drdi + delta_b0(i_mp_qr)*drdr ) &
4074  * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qr,i_mp_qi)*vrvi + theta_b0(i_mp_qr)*vrvr &
4075  + sigma_r + sigma_i )**0.5_rp
4076  pac(k,i_lracli2lg_i)= -0.25_rp*pi*e_ir*rhoq(k,i_nr)*rhoq(k,i_qi)*coef_acc_lri_i
4077  pac(k,i_nracni2ng_i)= -0.25_rp*pi*e_ir*rhoq(k,i_nr)*rhoq(k,i_ni)*coef_acc_nri_i
4078  !---- for charge density
4079  if (flg_lt) then !--- R + I -> G (decrease from both ice and rain charge, but only ice charge at this part)
4080  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NG is small, ignore charge transfer
4081  pcrg2(i_nracni2ng_i,k) = pac(k,i_nracni2ng_i)*(1.0_rp-sw1) &
4082  / (rhoq(k,i_ni)+sw1)*rhoq_crg(i_qi,k)
4083  end if
4084  ! reduction term of rain
4085  coef_acc_lri_r = &
4086  ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qi,i_mp_qr)*drdi + delta_b0(i_mp_qi)*didi ) &
4087  * ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qi,i_mp_qr)*vrvi + theta_b0(i_mp_qi)*vivi &
4088  + sigma_r + sigma_i )**0.5_rp
4089  coef_acc_nri_r = &
4090  ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qi,i_mp_qr)*drdi + delta_b0(i_mp_qi)*didi ) &
4091  * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qi,i_mp_qr)*vrvi + theta_b0(i_mp_qi)*vivi &
4092  + sigma_r + sigma_i )**0.5_rp
4093  pac(k,i_lracli2lg_r)= -0.25_rp*pi*e_ir*rhoq(k,i_ni)*rhoq(k,i_qr)*coef_acc_lri_r
4094  pac(k,i_nracni2ng_r)= -0.25_rp*pi*e_ir*rhoq(k,i_ni)*rhoq(k,i_nr)*coef_acc_nri_r
4095  !---- for charge density
4096  if (flg_lt) then !--- R + I -> G (decrease from both ice and rain charge, but only rain charge at this part)
4097  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small ) !--- if NG is small, ignore charge transfer
4098  pcrg2(i_nracni2ng_r,k) = pac(k,i_nracni2ng_r)*(1.0_rp-sw1)&
4099  / (rhoq(k,i_nr)+sw1)*rhoq_crg(i_qr,k)
4100  end if
4101  ! rain-snow => graupel
4102  ! reduction term of snow
4103  coef_acc_lrs_s = &
4104  ( delta_b1(i_mp_qs)*dsds + delta_ab1(i_mp_qr,i_mp_qs)*drds + delta_b0(i_mp_qr)*drdr ) &
4105  * ( theta_b1(i_mp_qs)*vsvs - theta_ab1(i_mp_qr,i_mp_qs)*vrvs + theta_b0(i_mp_qr)*vrvr &
4106  + sigma_r + sigma_s )**0.5_rp
4107  coef_acc_nrs_s = &
4108  ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qr,i_mp_qs)*drds + delta_b0(i_mp_qr)*drdr ) &
4109  * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qr,i_mp_qs)*vrvs + theta_b0(i_mp_qr)*vrvr &
4110  + sigma_r + sigma_s )**0.5_rp
4111  pac(k,i_lracls2lg_s)= -0.25_rp*pi*e_sr*rhoq(k,i_nr)*rhoq(k,i_qs)*coef_acc_lrs_s
4112  pac(k,i_nracns2ng_s)= -0.25_rp*pi*e_sr*rhoq(k,i_nr)*rhoq(k,i_ns)*coef_acc_nrs_s
4113  !---- for charge density
4114  if (flg_lt) then !--- R + S -> G (decrease from both snow and rain charge, but only snow charge at this part)
4115  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NS is small, ignore charge transfer
4116  pcrg2(i_nracns2ng_s,k) = pac(k,i_nracns2ng_s)*(1.0_rp-sw1)&
4117  / (rhoq(k,i_ns)+sw1)*rhoq_crg(i_qs,k)
4118  end if
4119  ! reduction term of rain
4120  coef_acc_lrs_r = &
4121  ( delta_b1(i_mp_qr)*drdr + delta_ab1(i_mp_qs,i_mp_qr)*drds + delta_b0(i_mp_qs)*dsds ) &
4122  * ( theta_b1(i_mp_qr)*vrvr - theta_ab1(i_mp_qs,i_mp_qr)*vrvs + theta_b0(i_mp_qs)*vsvs &
4123  + sigma_r + sigma_s )**0.5_rp
4124  coef_acc_nrs_r = &
4125  ( delta_b0(i_mp_qr)*drdr + delta_ab0(i_mp_qs,i_mp_qr)*drds + delta_b0(i_mp_qs)*dsds ) &
4126  * ( theta_b0(i_mp_qr)*vrvr - theta_ab0(i_mp_qs,i_mp_qr)*vrvs + theta_b0(i_mp_qs)*vsvs &
4127  + sigma_r + sigma_s )**0.5_rp
4128  pac(k,i_lracls2lg_r)= -0.25_rp*pi*e_sr*rhoq(k,i_ns)*rhoq(k,i_qr)*coef_acc_lrs_r
4129  pac(k,i_nracns2ng_r)= -0.25_rp*pi*e_sr*rhoq(k,i_ns)*rhoq(k,i_nr)*coef_acc_nrs_r
4130  !---- for charge density
4131  if (flg_lt) then !--- R + S -> G (decrease from both snow and rain charge, but only rain charge at this part)
4132  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small ) !--- if NR is small, ignore charge transfer
4133  pcrg2(i_nracns2ng_r,k) = pac(k,i_nracns2ng_r)*(1.0_rp-sw1) &
4134  / (rhoq(k,i_nr)+sw1)*rhoq_crg(i_qr,k)
4135  end if
4136  !------------------------------------------------------------------------
4137  !
4138  !+++ pattern 3: a + a => b (i-i)
4139  !
4140  !------------------------------------------------------------------------
4141  ! ice-ice ( reduction is double, but includes double-count)
4142  coef_acc_lii = &
4143  ( delta_b0(i_mp_qi)*didi + delta_ab1(i_mp_qi,i_mp_qi)*didi + delta_b1(i_mp_qi)*didi ) &
4144  * ( theta_b0(i_mp_qi)*vivi - theta_ab1(i_mp_qi,i_mp_qi)*vivi + theta_b1(i_mp_qi)*vivi &
4145  + sigma_i + sigma_i )**0.5_rp
4146  coef_acc_nii = &
4147  ( delta_b0(i_mp_qi)*didi + delta_ab0(i_mp_qi,i_mp_qi)*didi + delta_b0(i_mp_qi)*didi ) &
4148  * ( theta_b0(i_mp_qi)*vivi - theta_ab0(i_mp_qi,i_mp_qi)*vivi + theta_b0(i_mp_qi)*vivi &
4149  + sigma_i + sigma_i )**0.5_rp
4150  pac(k,i_liacli2ls)= -0.25_rp*pi*e_stick(k)*e_ii*rhoq(k,i_ni)*rhoq(k,i_qi)*coef_acc_lii
4151  pac(k,i_niacni2ns)= -0.25_rp*pi*e_stick(k)*e_ii*rhoq(k,i_ni)*rhoq(k,i_ni)*coef_acc_nii
4152  !---- for charge density
4153  if (flg_lt) then !--- I + I -> S (decrease from ice charge)
4154  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small, ignore charge transfer
4155  pcrg2(i_niacni2ns,k) = pac(k,i_niacni2ns)*(1.0_rp-sw1) &
4156  / (rhoq(k,i_ni)+sw1)*rhoq_crg(i_qi,k)
4157  end if
4158  !
4159 ! ci_aut(k) = 0.25_RP*pi*E_ii*rhoq(k,I_NI)*coef_acc_LII
4160 ! taui_aut(k) = 1._RP/max(E_stick(k)*ci_aut(k),1.E-10_RP)
4161 ! tau_sce(k) = rhoq(k,I_QI)/max(rhoq(k,I_QIj)+rhoq(k,I_QS),1.E-10_RP)
4162  !------------------------------------------------------------------------
4163  !
4164  !+++ pattern 4: a + a => a (s-s)
4165  !
4166  !------------------------------------------------------------------------
4167  ! snow-snow => snow
4168  coef_acc_nss = &
4169  ( delta_b0(i_mp_qs)*dsds + delta_ab0(i_mp_qs,i_mp_qs)*dsds + delta_b0(i_mp_qs)*dsds ) &
4170  * ( theta_b0(i_mp_qs)*vsvs - theta_ab0(i_mp_qs,i_mp_qs)*vsvs + theta_b0(i_mp_qs)*vsvs &
4171  + sigma_s + sigma_s )**0.5_rp
4172  pac(k,i_nsacns2ns)= -0.125_rp*pi*e_stick(k)*e_ss*rhoq(k,i_ns)*rhoq(k,i_ns)*coef_acc_nss
4173  !
4174  ! graupel-grauple => graupel
4175  coef_acc_ngg = &
4176  ( delta_b0(i_mp_qg)*dgdg + delta_ab0(i_mp_qg,i_mp_qg)*dgdg + delta_b0(i_mp_qg)*dgdg ) &
4177  * ( theta_b0(i_mp_qg)*vgvg - theta_ab0(i_mp_qg,i_mp_qg)*vgvg + theta_b0(i_mp_qg)*vgvg &
4178  + sigma_g + sigma_g )**0.5_rp
4179  pac(k,i_ngacng2ng)= -0.125_rp*pi*e_stick(k)*e_gg*rhoq(k,i_ng)*rhoq(k,i_ng)*coef_acc_ngg
4180  !
4181  !------------------------------------------------------------------------
4182  !--- Partial conversion
4183  ! SB06(70),(71)
4184  ! i_iconv2g: option whether partial conversions work or not
4185  ! ice-cloud => graupel
4186  sw = 0.5_rp - sign(0.5_rp,di_cri-ave_di) ! if( ave_di > di_cri )then sw=1
4187  wx_cri = cfill_i*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_di*ave_di*ave_di/xq(k,i_mp_qi) - 1.0_rp ) * sw
4188  pq(k,i_licon) = i_iconv2g * pac(k,i_liaclc2li)/max(1.0_rp, wx_cri) * sw
4189  pq(k,i_nicon) = i_iconv2g * pq(k,i_licon)/xq(k,i_mp_qi) * sw
4190  if (flg_lt) then !--- I + C -> G (decrease from ice charge)
4191  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small, ignore charge transfer
4192  pcrg1(i_nicon,k) = i_iconv2g * pq(k,i_nicon)*(1.0_rp-sw) &
4193  / (rhoq(k,i_ni)+sw1)*rhoq_crg(i_qi,k)
4194  end if
4195 
4196  ! snow-cloud => graupel
4197  wx_crs = cfill_s*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_ds*ave_ds*ave_ds/xq(k,i_mp_qs) - 1.0_rp )
4198  pq(k,i_lscon) = i_sconv2g * (pac(k,i_lsaclc2ls))/max(1.0_rp, wx_crs)
4199  pq(k,i_nscon) = i_sconv2g * pq(k,i_lscon)/xq(k,i_mp_qs)
4200  if (flg_lt) then !--- S + C -> G (decrease from snow charge)
4201  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NS is small, ignore charge transfer
4202  pcrg1(i_nscon,k) = i_sconv2g * pq(k,i_nscon)*(1.0_rp-sw) &
4203  / (rhoq(k,i_ns)+sw1)*rhoq_crg(i_qs,k)
4204  end if
4205  !------------------------------------------------------------------------
4206  !--- enhanced melting( due to collection-freezing of water droplets )
4207  ! originally from Rutledge and Hobbs(1984). eq.(A.21)
4208  ! if T > 273.15 then temc_p=T-273.15, else temc_p=0
4209  ! 08/05/08 [fix] T.Mitsui LHF00 => LHF0
4210  ! melting occurs around T=273K, so LHF0 is suitable both SIMPLE and EXACT,
4211  ! otherwise LHF can have sign both negative(EXACT) and positive(SIMPLE).
4212 !!$ coef_emelt = -CL/LHF00*temc_p
4213  coef_emelt = cl/lhf0*temc_p
4214  ! cloud-graupel
4215  pq(k,i_lgacm) = coef_emelt*pac(k,i_lgaclc2lg)
4216  pq(k,i_ngacm) = pq(k,i_lgacm)/xq(k,i_mp_qg)
4217  if (flg_lt) then
4218  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small ) !--- if NG is small, ignore charge transfer
4219  pcrg1(i_ngacm,k) = pq(k,i_ngacm)*(1.0_rp-sw1) &
4220  / (rhoq(k,i_ng)+sw1)*rhoq_crg(i_qg,k)
4221  end if
4222  ! rain-graupel
4223  pq(k,i_lgarm) = coef_emelt*pac(k,i_lraclg2lg)
4224  pq(k,i_ngarm) = pq(k,i_lgarm)/xq(k,i_mp_qg)
4225  if (flg_lt) then
4226  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small ) !--- if NG is small, ignore charge transfer
4227  pcrg1(i_ngarm,k) = pq(k,i_ngarm)*(1.0_rp-sw1) &
4228  / (rhoq(k,i_ng)+sw1)*rhoq_crg(i_qg,k)
4229  end if
4230  ! cloud-snow
4231  pq(k,i_lsacm) = coef_emelt*(pac(k,i_lsaclc2ls))
4232  pq(k,i_nsacm) = pq(k,i_lsacm)/xq(k,i_mp_qs)
4233  if (flg_lt) then
4234  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NG is small, ignore charge transfer
4235  pcrg1(i_nsacm,k) = pq(k,i_nsacm)*(1.0_rp-sw1) &
4236  / (rhoq(k,i_ns)+sw1)*rhoq_crg(i_qs,k)
4237  end if
4238  ! rain-snow
4239  pq(k,i_lsarm) = coef_emelt*(pac(k,i_lracls2lg_r)+pac(k,i_lracls2lg_s))
4240  pq(k,i_nsarm) = pq(k,i_lsarm)/xq(k,i_mp_qg)
4241  if (flg_lt) then
4242  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NG is small, ignore charge transfer
4243  pcrg1(i_nsarm,k) = pq(k,i_nsarm)*(1.0_rp-sw1) &
4244  / (rhoq(k,i_ns)+sw1)*rhoq_crg(i_qs,k)
4245  end if
4246  ! cloud-ice
4247  pq(k,i_liacm) = coef_emelt*pac(k,i_liaclc2li)
4248  pq(k,i_niacm) = pq(k,i_liacm)/xq(k,i_mp_qi)
4249  if (flg_lt) then
4250  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NG is small, ignore charge transfer
4251  pcrg1(i_niacm,k) = pq(k,i_niacm)*(1.0_rp-sw1) &
4252  / (rhoq(k,i_ni)+sw1)*rhoq_crg(i_qi,k)
4253  end if
4254  ! rain-ice
4255  pq(k,i_liarm) = coef_emelt*(pac(k,i_lracli2lg_r)+pac(k,i_lracli2lg_i))
4256  pq(k,i_niarm) = pq(k,i_liarm)/xq(k,i_mp_qg)
4257  if (flg_lt) then
4258  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small ) !--- if NG is small, ignore charge transfer
4259  pcrg1(i_niarm,k) = pq(k,i_niarm)*(1.0_rp-sw1) &
4260  / (rhoq(k,i_nr)+sw1)*rhoq_crg(i_qi,k)
4261  end if
4262  end do
4263  !
4264  return

Referenced by atmos_phy_mp_sn14_terminal_velocity().

Here is the caller graph for this function:

◆ aut_acc_slc_brk()

subroutine scale_atmos_phy_mp_sn14::aut_acc_slc_brk ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
logical, intent(in)  flg_lt,
real(rp), dimension(ka,i_qv:i_ng), intent(in)  rhoq,
real(rp), dimension(i_qc:i_qg,ka), intent(in)  rhoq_crg,
real(rp), dimension(ka,hydro_max), intent(in)  xq,
real(rp), dimension(ka,hydro_max), intent(in)  dq_xave,
real(rp), dimension(ka), intent(in)  rho,
real(rp), dimension(ka,pq_max), intent(inout)  PQ,
real(rp), dimension(pq_max,ka), intent(inout)  Pcrg 
)

Definition at line 4276 of file scale_atmos_phy_mp_sn14.F90.

4276  implicit none
4277 
4278  integer, intent(in) :: KA, KS, KE
4279  !
4280  real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
4281  real(RP), intent(in) :: rhoq_crg(I_QC:I_QG,KA)
4282  logical, intent(in) :: flg_lt
4283  real(RP), intent(in) :: xq(KA,HYDRO_MAX)
4284  real(RP), intent(in) :: dq_xave(KA,HYDRO_MAX)
4285  real(RP), intent(in) :: rho(KA)
4286  !
4287  real(RP), intent(inout) :: PQ(KA,PQ_MAX)
4288  real(RP), intent(inout) :: Pcrg(PQ_MAX,KA)
4289  !
4290  ! parameter for autoconversion
4291  real(RP), parameter :: kcc = 4.44e+9_rp ! collision efficiency [m3/kg2/sec]
4292  real(RP), parameter :: tau_min = 1.e-20_rp ! empirical filter by T.Mitsui
4293  real(RP), parameter :: rx_sep = 1.0_rp/x_sep ! 1/x_sep, 10/08/03 [Add] T.Mitsui
4294  !
4295  ! parameter for accretion
4296  real(RP), parameter :: kcr = 5.8_rp ! collision efficiency [m3/kg2/sec]
4297  real(RP), parameter :: thr_acc = 5.e-5_rp ! threshold for universal function original
4298  !
4299  ! parameter for self collection and collison break-up
4300  real(RP), parameter :: krr = 4.33_rp ! k_rr, S08 (35)
4301  real(RP), parameter :: kaprr = 60.7_rp ! kappa_rr, SB06(11)
4302  real(RP), parameter :: kbr = 1000._rp ! k_br, SB06(14)
4303  real(RP), parameter :: kapbr = 2.3e+3_rp ! kappa_br, SB06(15)
4304  real(RP), parameter :: dr_min = 0.35e-3_rp ! minimum diameter, SB06(13)-(15)
4305  !
4306  ! work variables
4307  real(RP) :: coef_nuc0 ! coefficient of number for Auto-conversion
4308  real(RP) :: coef_nuc1 ! mass
4309  real(RP) :: coef_aut0 ! number
4310  real(RP) :: coef_aut1 ! mass
4311  real(RP) :: lwc ! lc+lr
4312  real(RP) :: tau ! conversion ratio: qr/(qc+qr) ranges [0:1]
4313  real(RP) :: rho_fac ! factor of air density
4314  real(RP) :: psi_aut ! Universal function of Auto-conversion
4315  real(RP) :: psi_acc ! Universal function of Accretion
4316  real(RP) :: psi_brk ! Universal function of Breakup
4317  real(RP) :: ddr ! diameter difference from equilibrium
4318  !
4319  integer :: k, iqw
4320  real(RP) :: sw
4321  !
4322  coef_nuc0 = (nu(i_mp_qc)+2.0_rp)/(nu(i_mp_qc)+1.0_rp)
4323  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)
4324  coef_aut0 = -kcc*coef_nuc0
4325  coef_aut1 = -kcc/x_sep/20._rp*coef_nuc1
4326  !
4327  do k = ks, ke
4328  lwc = rhoq(k,i_qr) + rhoq(k,i_qc)
4329  if( lwc > xc_min )then
4330  tau = max(tau_min, rhoq(k,i_qr)/lwc)
4331  else
4332  tau = tau_min
4333  end if
4334  rho_fac = sqrt(rho_0/max(rho(k),rho_min))
4335  !
4336  ! Auto-conversion ( cloud-cloud => rain )
4337  psi_aut = 400._rp*(tau**0.7_rp)*(1.0_rp - (tau**0.7_rp))**3 ! (6) SB06
4338  pq(k,i_ncaut) = coef_aut0*rhoq(k,i_qc)*rhoq(k,i_qc)*rho_fac*rho_fac ! (9) SB06 sc+aut
4339  ! lc = lwc*(1-tau), lr = lwc*tau
4340  pq(k,i_lcaut) = coef_aut1*lwc*lwc*xq(k,i_mp_qc)*xq(k,i_mp_qc) & ! (4) SB06
4341  *((1.0_rp-tau)*(1.0_rp-tau) + psi_aut)*rho_fac*rho_fac !
4342  pq(k,i_nraut) = -rx_sep*pq(k,i_lcaut) ! (A7) SB01
4343  !--- for Charge density
4344  if (flg_lt) then
4345  sw = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
4346  pcrg(i_ncaut,k) = pq(k,i_ncaut)*(1.0_rp-sw)/(rhoq(k,i_nc)+sw)*rhoq_crg(i_qc,k)
4347  pcrg(i_nraut,k) = -pcrg(i_ncaut,k)
4348  end if
4349  !
4350  ! Accretion ( cloud-rain => rain )
4351  psi_acc =(tau/(tau+thr_acc))**4 ! (8) SB06
4352  pq(k,i_lcacc) = -kcr*rhoq(k,i_qc)*rhoq(k,i_qr)*rho_fac*psi_acc ! (7) SB06
4353  pq(k,i_ncacc) = -kcr*rhoq(k,i_nc)*rhoq(k,i_qr)*rho_fac*psi_acc ! (A6) SB01
4354  !--- for Charge density
4355  if (flg_lt) then
4356  sw = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
4357  pcrg(i_ncacc,k) = pq(k,i_ncacc)*(1.0_rp-sw)/(rhoq(k,i_nc)+sw)*rhoq(k,i_qc)
4358  end if
4359  !
4360  ! Self-collection ( rain-rain => rain )
4361  pq(k,i_nrslc) = -krr*rhoq(k,i_nr)*rhoq(k,i_qr)*rho_fac ! (A.8) SB01
4362  !
4363  ! Collisional breakup of rain
4364  ddr = min(1.e-3_rp, dq_xave(k,i_mp_qr) - dr_eq )
4365  if ( dq_xave(k,i_mp_qr) < dr_min ) then ! negligible
4366  psi_brk = -1.0_rp
4367  else if ( dq_xave(k,i_mp_qr) <= dr_eq ) then
4368  psi_brk = kbr*ddr ! (14) SB06
4369  else
4370  psi_brk = exp(kapbr*ddr) - 1.0_rp ! (15) SB06 (SB06 has a typo)
4371  end if
4372  pq(k,i_nrbrk) = - (psi_brk + 1.0_rp)*pq(k,i_nrslc) ! (13) SB06
4373  !
4374  end do
4375  !
4376  return

References scale_const::const_eps.

Referenced by atmos_phy_mp_sn14_terminal_velocity().

Here is the caller graph for this function:

◆ freezing_water()

subroutine scale_atmos_phy_mp_sn14::freezing_water ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
real(rp), intent(in)  dt,
real(rp), dimension(ka,i_qv:i_ng), intent(in)  rhoq,
real(rp), dimension(ka,hydro_max), intent(in)  xq,
real(rp), dimension(ka), intent(in)  tem,
real(rp), dimension(ka,pq_max), intent(inout)  PQ 
)

Definition at line 4619 of file scale_atmos_phy_mp_sn14.F90.

4619  implicit none
4620  !
4621  ! In this subroutine,
4622  ! We assumed surface temperature of droplets are same as environment.
4623 
4624  integer, intent(in) :: KA, KS, KE
4625 
4626  real(RP), intent(in) :: dt
4627  !
4628  real(RP), intent(in) :: tem(KA)
4629  !
4630  real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
4631  real(RP), intent(in) :: xq(KA,HYDRO_MAX)
4632  !
4633  real(RP), intent(inout):: PQ(KA,PQ_MAX)
4634  !
4635  real(RP), parameter :: temc_min = -65.0_rp
4636  real(RP), parameter :: a_het = 0.2_rp ! SB06 (44)
4637  real(RP), parameter :: b_het = 0.65_rp ! SB06 (44)
4638  !
4639  real(RP) :: coef_m2_c
4640  real(RP) :: coef_m2_r
4641  ! temperature [celsius]
4642  real(RP) :: temc, temc2, temc3, temc4
4643  ! temperature function of homegenous/heterogenous freezing
4644  real(RP) :: Jhom, Jhet
4645  real(RP) :: rdt
4646  real(RP) :: tmp
4647  !
4648  integer :: k
4649  !
4650  rdt = 1.0_rp/dt
4651  !
4652  coef_m2_c = coef_m2(i_mp_qc)
4653  coef_m2_r = coef_m2(i_mp_qr)
4654  !
4655 
4656  do k = ks, ke
4657  temc = max( tem(k) - t00, temc_min )
4658  ! These cause from aerosol-droplet interaction.
4659  ! Bigg(1953) formula, Khain etal.(2000) eq.(4.5), Pruppacher and Klett(1997) eq.(9-48)
4660  jhet = a_het*(exp( -b_het*temc) - 1.0_rp )
4661  ! These cause in nature.
4662  ! Cotton and Field 2002, QJRMS. (12)
4663  if( temc < -65.0_rp )then
4664  jhom = 10.0_rp**(24.37236_rp)*1.e+3_rp
4665  jhet = a_het*(exp( 65.0_rp*b_het) - 1.0_rp )
4666  else if( temc < -30.0_rp ) then
4667  temc2 = temc*temc
4668  temc3 = temc*temc2
4669  temc4 = temc2*temc2
4670  jhom = 10.0_rp**(&
4671  - 243.40_rp - 14.75_rp*temc - 0.307_rp*temc2 &
4672  - 0.00287_rp*temc3 - 0.0000102_rp*temc4 ) *1.e+3_rp
4673  else if( temc < 0.0_rp) then
4674  jhom = 10._rp**(-7.63_rp-2.996_rp*(temc+30.0_rp))*1.e+3_rp
4675  else
4676  jhom = 0.0_rp
4677  jhet = 0.0_rp
4678  end if
4679  ! Note, xc should be limited in range[xc_min:xc_max].
4680  ! and PNChom need to be calculated by NC
4681  ! because reduction rate of Nc need to be bound by NC.
4682  ! For the same reason PLChom also bound by LC and xc.
4683  ! Basically L and N should be independent
4684  ! but average particle mass x should be in suitable range.
4685  ! Homogenous Freezing
4686  pq(k,i_lchom) = 0.0_rp
4687  pq(k,i_nchom) = 0.0_rp
4688  ! Heterogenous Freezing
4689 #if defined(PGI) || defined(SX)
4690  tmp = min( xq(k,i_mp_qc)*(jhet+jhom)*dt, 1.e+3_rp) ! apply exp limiter
4691  pq(k,i_lchet) = -rdt*rhoq(k,i_qc)*( 1.0_rp - exp( -coef_m2_c*tmp ) )
4692  pq(k,i_nchet) = -rdt*rhoq(k,i_nc)*( 1.0_rp - exp( - tmp ) )
4693 
4694  tmp = min( xq(k,i_mp_qr)*(jhet+jhom)*dt, 1.e+3_rp) ! apply exp limiter
4695  pq(k,i_lrhet) = -rdt*rhoq(k,i_qr)*( 1.0_rp - exp( -coef_m2_r*tmp ) )
4696  pq(k,i_nrhet) = -rdt*rhoq(k,i_nr)*( 1.0_rp - exp( - tmp ) )
4697 #else
4698  pq(k,i_lchet) = -rdt*rhoq(k,i_qc)*( 1.0_rp - exp( -coef_m2_c*xq(k,i_mp_qc)*(jhet+jhom)*dt ) )
4699  pq(k,i_nchet) = -rdt*rhoq(k,i_nc)*( 1.0_rp - exp( - xq(k,i_mp_qc)*(jhet+jhom)*dt ) )
4700  pq(k,i_lrhet) = -rdt*rhoq(k,i_qr)*( 1.0_rp - exp( -coef_m2_r*xq(k,i_mp_qr)*(jhet+jhom)*dt ) )
4701  pq(k,i_nrhet) = -rdt*rhoq(k,i_nr)*( 1.0_rp - exp( - xq(k,i_mp_qr)*(jhet+jhom)*dt ) )
4702 #endif
4703  end do
4704 
4705  !
4706  return

References scale_atmos_hydrometeor::cp_ice, scale_atmos_hydrometeor::cp_vapor, scale_atmos_hydrometeor::cp_water, scale_atmos_hydrometeor::cv_ice, scale_atmos_hydrometeor::cv_vapor, and scale_atmos_hydrometeor::cv_water.

Referenced by atmos_phy_mp_sn14_terminal_velocity().

Here is the caller graph for this function:

◆ cross_section()

subroutine scale_atmos_phy_mp_sn14::cross_section ( real(rp), dimension (ka,hydro_max), intent(out)  Crs,
integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  QA_MP,
real(rp), dimension(ka,qa_mp), intent(in)  QTRC0,
real(rp), dimension(ka), intent(in)  DENS0 
)

Calculate Cross Section.

Definition at line 5367 of file scale_atmos_phy_mp_sn14.F90.

5367  implicit none
5368 
5369  integer, intent(in) :: KA, KS, KE
5370  integer, intent(in) :: QA_MP
5371  real(RP), intent(out) :: Crs (KA,HYDRO_MAX)! Cross section [cm]
5372  real(RP), intent(in) :: QTRC0(KA,QA_MP) ! tracer mass concentration [kg/kg]
5373  real(RP), intent(in) :: DENS0(KA) ! density [kg/m3]
5374 
5375  ! mass concentration[kg/m3] and mean particle mass[kg]
5376  real(RP) :: xc(KA)
5377  real(RP) :: xr(KA)
5378  real(RP) :: xi(KA)
5379  real(RP) :: xs(KA)
5380  real(RP) :: xg(KA)
5381  ! diameter of average mass[kg/m3]
5382  real(RP) :: dc_ave(KA)
5383  real(RP) :: dr_ave(KA)
5384  ! radius of average mass
5385  real(RP) :: rc, rr
5386 
5387  real(RP) :: coef_Fuetal1998
5388  ! r2m_min is minimum value(moment of 1 particle with 1 micron)
5389  real(RP), parameter :: r2m_min=1.e-12_rp
5390  real(RP), parameter :: um2cm = 100.0_rp
5391 
5392  real(RP) :: limitsw, zerosw
5393  integer :: k
5394  !---------------------------------------------------------------------------
5395 
5396  ! mean particle mass[kg]
5397  do k = ks, ke
5398  xc(k) = min( xc_max, max( xc_min, dens0(k)*qtrc0(k,i_qc)/(qtrc0(k,i_nc)+nc_min) ) )
5399  xr(k) = min( xr_max, max( xr_min, dens0(k)*qtrc0(k,i_qr)/(qtrc0(k,i_nr)+nr_min) ) )
5400  xi(k) = min( xi_max, max( xi_min, dens0(k)*qtrc0(k,i_qi)/(qtrc0(k,i_ni)+ni_min) ) )
5401  xs(k) = min( xs_max, max( xs_min, dens0(k)*qtrc0(k,i_qs)/(qtrc0(k,i_ns)+ns_min) ) )
5402  xg(k) = min( xg_max, max( xg_min, dens0(k)*qtrc0(k,i_qg)/(qtrc0(k,i_ng)+ng_min) ) )
5403  enddo
5404 
5405  ! diameter of average mass : SB06 eq.(32)
5406  do k = ks, ke
5407  dc_ave(k) = a_m(i_qc) * xc(k)**b_m(i_qc)
5408  dr_ave(k) = a_m(i_qr) * xr(k)**b_m(i_qr)
5409  enddo
5410 
5411  do k = ks, ke
5412  crs(k,i_mp_qc) = pi * coef_r2(i_mp_qc) * qtrc0(k,i_nc) * a_rea2(i_mp_qc) * xc(k)**b_rea2(i_mp_qc)
5413  enddo
5414 
5415  do k = ks, ke
5416  crs(k,i_mp_qr) = pi * coef_r2(i_mp_qr) * qtrc0(k,i_nr) * a_rea2(i_mp_qr) * xr(k)**b_rea2(i_mp_qr)
5417  enddo
5418 
5419  do k = ks, ke
5420  crs(k,i_mp_qi) = pi * coef_rea2(i_mp_qi) * qtrc0(k,i_ni) * a_rea2(i_mp_qi) * xi(k)**b_rea2(i_mp_qi)
5421  crs(k,i_mp_qs) = pi * coef_rea2(i_mp_qs) * qtrc0(k,i_ns) * a_rea2(i_mp_qs) * xs(k)**b_rea2(i_mp_qs)
5422  crs(k,i_mp_qg) = pi * coef_rea2(i_mp_qg) * qtrc0(k,i_ng) * a_rea2(i_mp_qg) * xg(k)**b_rea2(i_mp_qg)
5423  enddo
5424 
5425  return

Referenced by atmos_phy_mp_sn14_terminal_velocity().

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 109 of file scale_atmos_phy_mp_sn14.F90.

109  integer, public, parameter :: QA_MP = 11

Referenced by atmos_phy_mp_sn14_cloud_fraction(), and atmos_phy_mp_sn14_terminal_velocity().

◆ atmos_phy_mp_sn14_ntracers

integer, parameter, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_ntracers = QA_MP

Definition at line 111 of file scale_atmos_phy_mp_sn14.F90.

111  integer, parameter, public :: ATMOS_PHY_MP_SN14_ntracers = qa_mp

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

◆ atmos_phy_mp_sn14_nwaters

integer, parameter, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_nwaters = 2

Definition at line 112 of file scale_atmos_phy_mp_sn14.F90.

112  integer, parameter, public :: ATMOS_PHY_MP_SN14_nwaters = 2

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

◆ atmos_phy_mp_sn14_nices

integer, parameter, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_nices = 3

Definition at line 113 of file scale_atmos_phy_mp_sn14.F90.

113  integer, parameter, public :: ATMOS_PHY_MP_SN14_nices = 3

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

◆ atmos_phy_mp_sn14_tracer_names

character(len=h_short), dimension(qa_mp), parameter, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_tracer_names = (/ 'QV', 'QC', 'QR', 'QI', 'QS', 'QG', 'NC', 'NR', 'NI', 'NS', 'NG' /)

Definition at line 114 of file scale_atmos_phy_mp_sn14.F90.

114  character(len=H_SHORT), parameter, public :: ATMOS_PHY_MP_SN14_tracer_names(QA_MP) = (/ &
115  'QV', &
116  'QC', &
117  'QR', &
118  'QI', &
119  'QS', &
120  'QG', &
121  'NC', &
122  'NR', &
123  'NI', &
124  'NS', &
125  'NG' /)

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

◆ atmos_phy_mp_sn14_tracer_descriptions

character(len=h_mid), dimension(qa_mp), parameter, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_tracer_descriptions = (/ 'Ratio of Water Vapor mass to total mass (Specific humidity)', 'Ratio of Cloud Water mass to total mass ', 'Ratio of Rain Water mass to total mass ', 'Ratio of Cloud Ice mass ratio to total mass ', 'Ratio of Snow miass ratio to total mass ', 'Ratio of Graupel mass ratio to total mass ', 'Cloud Water Number Density ', 'Rain Water Number Density ', 'Cloud Ice Number Density ', 'Snow Number Density ', 'Graupel Number Density '/)

Definition at line 126 of file scale_atmos_phy_mp_sn14.F90.

126  character(len=H_MID) , parameter, public :: ATMOS_PHY_MP_SN14_tracer_descriptions(QA_MP) = (/ &
127  'Ratio of Water Vapor mass to total mass (Specific humidity)', &
128  'Ratio of Cloud Water mass to total mass ', &
129  'Ratio of Rain Water mass to total mass ', &
130  'Ratio of Cloud Ice mass ratio to total mass ', &
131  'Ratio of Snow miass ratio to total mass ', &
132  'Ratio of Graupel mass ratio to total mass ', &
133  'Cloud Water Number Density ', &
134  'Rain Water Number Density ', &
135  'Cloud Ice Number Density ', &
136  'Snow Number Density ', &
137  'Graupel Number Density '/)

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

◆ atmos_phy_mp_sn14_tracer_units

character(len=h_short), dimension(qa_mp), parameter, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_tracer_units = (/ 'kg/kg ', 'kg/kg ', 'kg/kg ', 'kg/kg ', 'kg/kg ', 'kg/kg ', 'num/kg', 'num/kg', 'num/kg', 'num/kg', 'num/kg' /)

Definition at line 138 of file scale_atmos_phy_mp_sn14.F90.

138  character(len=H_SHORT), parameter, public :: ATMOS_PHY_MP_SN14_tracer_units(QA_MP) = (/ &
139  'kg/kg ', &
140  'kg/kg ', &
141  'kg/kg ', &
142  'kg/kg ', &
143  'kg/kg ', &
144  'kg/kg ', &
145  'num/kg', &
146  'num/kg', &
147  'num/kg', &
148  'num/kg', &
149  'num/kg' /)

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
scale_atmos_hydrometeor::i_hr
integer, parameter, public i_hr
liquid water rain
Definition: scale_atmos_hydrometeor.F90:82
scale_atmos_hydrometeor::i_hs
integer, parameter, public i_hs
snow
Definition: scale_atmos_hydrometeor.F90:84
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_prc::prc_myrank
integer, public prc_myrank
process num in local communicator
Definition: scale_prc.F90:90
scale_atmos_hydrometeor::i_hh
integer, parameter, public i_hh
hail
Definition: scale_atmos_hydrometeor.F90:86
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:31
scale_atmos_hydrometeor::i_hi
integer, parameter, public i_hi
ice water cloud
Definition: scale_atmos_hydrometeor.F90:83
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_specfunc
module SPECFUNC
Definition: scale_specfunc.F90:14
scale_atmos_hydrometeor::i_hc
integer, parameter, public i_hc
liquid water cloud
Definition: scale_atmos_hydrometeor.F90:81
scale_specfunc::sf_gamma
real(rp) function, public sf_gamma(x)
Gamma function.
Definition: scale_specfunc.F90:50
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:41
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
scale_atmos_hydrometeor::i_hg
integer, parameter, public i_hg
graupel
Definition: scale_atmos_hydrometeor.F90:85
Generated by   doxygen 1.8.17