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_finalize
 finalize 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, cva, dTdt_rad, qke, CCN, nc_uplim_d, dt, dq_xa, vt_xa, PQ)
 
subroutine nucleation_ice_hom (KA, KS, KE, tem, pre, rho, qd, rhoq_qv, cva, cpa, w, dTdt_rad, dTdt_dep, PLIdep, dt, PLIhom, PNIhom)
 
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 mixed_phase_collection_bin (KA, KS, KE, flg_lt, d0_crg, v0_crg, beta_crg, dqcrg, wtem, rhoq, rhoq_crg, xq, dq_xave, vt_xave, rho, PQ, Pcrg1, Pcrg2, Pac)
 
subroutine aut_acc_slc_brk (KA, KS, KE, flg_lt, rhoq, rhoq_crg, xq, dq_xave, rho, PQ, Pcrg)
 
subroutine dep_vapor_ice_wrk (KA, KS, KE, PLIdep_total, rho, tem, pre, qd, esi, qsi, rhoq, vt_xave, dq_xave, dt)
 
subroutine freezing_water (KA, KS, KE, dt, rhoq, xq, tem, PQ)
 
subroutine cross_section (KA, KS, KE, QA_MP, QTRC0, DENS0, Crs)
 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 mass 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.
    OPT_COLLECTION_BIN logical .false. SO22
    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)
    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.
    SO22_HET logical .false. SO22
    OPT_NUCLEATION_ICE_HOM logical .false. SO22

  • 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_collection_bin
    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
    OPT_STICK_KS96 logical .false.
    OPT_STICK_CO86 logical .false.
    TEM_MIN_ESTICK real(RP) 253.0_RP
    OPT_STICK_RHH57 logical .false.
    OPT_STICK_RHKS96 logical .false.
    OPT_STICK_C12 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
namedescriptionunitvariable
I_CGNGacNI2NG individual tendency term in SN14 kg/kg/s I_CGNGacNI2NG

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

716  use scale_prc, only: &
717  prc_abort
718  use scale_file_history, only: &
719  file_history_reg
720  implicit none
721  integer, intent(in) :: KA
722  integer, intent(in) :: IA
723  integer, intent(in) :: JA
724 
725  namelist / param_atmos_phy_mp_sn14 / &
726  mp_doautoconversion, &
727  mp_ssw_lim, &
728  mp_couple_aerosol
729 
730  integer :: ip
731  integer :: ierr
732  !---------------------------------------------------------------------------
733 
734  log_newline
735  log_info("ATMOS_PHY_MP_sn14_setup",*) 'Setup'
736  log_info("ATMOS_PHY_MP_sn14_setup",*) 'Seiki and Nakajima (2014) 2-moment bulk 6 category'
737 
738  !--- read namelist
739  rewind(io_fid_conf)
740  read(io_fid_conf,nml=param_atmos_phy_mp_sn14,iostat=ierr)
741  if( ierr < 0 ) then !--- missing
742  log_info("ATMOS_PHY_MP_sn14_setup",*) 'Not found namelist. Default used.'
743  elseif( ierr > 0 ) then !--- fatal error
744  log_error("ATMOS_PHY_MP_sn14_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SN14. Check!'
745  call prc_abort
746  endif
747  log_nml(param_atmos_phy_mp_sn14)
748 
749  call mp_sn14_init
750 
751  allocate(nc_uplim_d(1,ia,ja))
752  nc_uplim_d(:,:,:) = 150.e6_rp
753 
754  hist_max = 0
755  hist_idx(:) = -1
756  do ip = 1, w_nmax
757  call file_history_reg( w_name(ip), 'individual tendency term in SN14', 'kg/kg/s', &
758  hist_id(ip) )
759  if ( hist_id(ip) > 0 ) then
760  hist_max = hist_max + 1
761  hist_idx(ip) = hist_max
762  end if
763  enddo
764  allocate( w3d(ka,ia,ja,hist_max) )
765  w3d(:,:,:,:) = 0.0_rp
766 
767  return

References scale_file_history::file_history_reg(), 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_finalize()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_finalize

finalize

Definition at line 773 of file scale_atmos_phy_mp_sn14.F90.

773 
774  deallocate(nc_uplim_d)
775 
776  return

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_finalize().

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

809  implicit none
810 
811  integer, intent(in) :: KA, KS, KE
812  integer, intent(in) :: IA, IS, IE
813  integer, intent(in) :: JA, JS, JE
814 
815  real(RP), intent(in) :: DENS (KA,IA,JA)
816  real(RP), intent(in) :: W (KA,IA,JA)
817  real(RP), intent(in) :: QTRC (KA,IA,JA,QA_MP)
818  real(RP), intent(in) :: PRES(KA,IA,JA)
819  real(RP), intent(in) :: TEMP(KA,IA,JA)
820  real(RP), intent(in) :: Qdry(KA,IA,JA)
821  real(RP), intent(in) :: CPtot(KA,IA,JA)
822  real(RP), intent(in) :: CVtot(KA,IA,JA)
823  real(RP), intent(in) :: CCN (KA,IA,JA)
824  real(DP), intent(in) :: dt
825  real(RP), intent(in) :: cz( KA,IA,JA)
826  real(RP), intent(in) :: fz(0:KA,IA,JA)
827 
828  real(RP), intent(out) :: RHOQ_t (KA,IA,JA,QA_MP)
829  real(RP), intent(out) :: RHOE_t (KA,IA,JA)
830  real(RP), intent(out) :: CPtot_t(KA,IA,JA)
831  real(RP), intent(out) :: CVtot_t(KA,IA,JA)
832  real(RP), intent(out) :: EVAPORATE(KA,IA,JA) !--- number of evaporated cloud [/m3]
833 
834  ! Optional for Lightning
835  logical, intent(in), optional :: flg_lt
836  real(RP), intent(in), optional :: d0_crg, v0_crg
837  real(RP), intent(in), optional :: dqcrg(KA,IA,JA)
838  real(RP), intent(in), optional :: beta_crg(KA,IA,JA)
839  real(RP), intent(in), optional :: QTRC_crg(KA,IA,JA,HYDRO_MAX)
840  real(RP), intent(out), optional :: QSPLT_in(KA,IA,JA,3)
841  real(RP), intent(out), optional :: Sarea(KA,IA,JA,HYDRO_MAX)
842  real(RP), intent(out), optional :: RHOQcrg_t(KA,IA,JA,HYDRO_MAX)
843  !---------------------------------------------------------------------------
844 
845  log_progress(*) 'atmosphere / physics / microphysics / SN14'
846 
847 #ifdef PROFILE_FIPP
848  call fipp_start()
849 #endif
850 
851  !##### MP Main #####
852  call mp_sn14 ( &
853  ka, ks, ke, ia, is, ie, ja, js, je, &
854  dens(:,:,:), w(:,:,:), qtrc(:,:,:,:), pres(:,:,:), temp(:,:,:), & ! (in)
855  qdry(:,:,:), cptot(:,:,:), cvtot(:,:,:), ccn(:,:,:), & ! (in)
856  real(dt,RP), cz(:,:,:), fz(:,:,:), & ! (in)
857  RHOQ_t(:,:,:,:), RHOE_t(:,:,:), CPtot_t(:,:,:), CVtot_t(:,:,:), & ! (out)
858  EVAPORATE(:,:,:), & ! (out)
859  flg_lt, d0_crg, v0_crg, dqcrg(:,:,:), beta_crg(:,:,:), & ! (optional in)
860  QTRC_crg(:,:,:,:), & ! (optional in)
861  QSPLT_in(:,:,:,:), Sarea(:,:,:,:), RHOQcrg_t(:,:,:,:) ) ! (optional out)
862 
863 #ifdef PROFILE_FIPP
864  call fipp_stop()
865 #endif
866 
867  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 879 of file scale_atmos_phy_mp_sn14.F90.

879  implicit none
880  integer, intent(in) :: KA, KS, KE
881  integer, intent(in) :: IA, IS, IE
882  integer, intent(in) :: JA, JS, JE
883 
884  real(RP), intent(in) :: QTRC (KA,IA,JA,QA_MP-1)
885  real(RP), intent(in) :: mask_criterion
886 
887  real(RP), intent(out) :: cldfrac(KA,IA,JA)
888 
889  real(RP) :: qhydro
890  integer :: k, i, j, iq
891  !---------------------------------------------------------------------------
892 
893  !$omp parallel do &
894  !$omp private(qhydro)
895  do j = js, je
896  do i = is, ie
897  do k = ks, ke
898  qhydro = 0.0_rp
899  do iq = i_mp_qc, i_mp_qg
900  qhydro = qhydro + qtrc(k,i,j,iq)
901  enddo
902  cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
903  enddo
904  enddo
905  enddo
906 
907  return

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

917  use scale_atmos_hydrometeor, only: &
918  n_hyd, &
919  i_hc, &
920  i_hr, &
921  i_hi, &
922  i_hs, &
923  i_hg, &
924  i_hh
925  implicit none
926  integer, intent(in) :: KA, KS, KE
927  integer, intent(in) :: IA, IS, IE
928  integer, intent(in) :: JA, JS, JE
929 
930  real(RP), intent(in) :: DENS0(KA,IA,JA) ! density [kg/m3]
931  real(RP), intent(in) :: TEMP0(KA,IA,JA) ! temperature [K]
932  real(RP), intent(in) :: QTRC0(KA,IA,JA,I_QC:I_NG) ! tracer mass concentration [kg/kg]
933 
934  real(RP), intent(out) :: Re (KA,IA,JA,N_HYD) ! effective radius [cm]
935 
936  ! mass concentration[kg/m3] and mean particle mass[kg]
937  real(RP) :: xc(KA)
938  real(RP) :: xr(KA)
939  real(RP) :: xi(KA)
940  real(RP) :: xs(KA)
941  real(RP) :: xg(KA)
942  ! diameter of average mass[kg/m3]
943  real(RP) :: dc_ave(KA)
944  real(RP) :: dr_ave(KA)
945  ! radius of average mass
946  real(RP) :: rc, rr
947  ! 2nd. and 3rd. order moment of DSD
948  real(RP) :: ri2m(KA), ri3m(KA)
949  real(RP) :: rs2m(KA), rs3m(KA)
950  real(RP) :: rg2m(KA), rg3m(KA)
951 
952  real(RP), parameter :: coef_Fuetal1998 = 3.0_rp / (4.0_rp*rhoi)
953 
954  ! r2m_min is minimum value(moment of 1 particle with 1 micron)
955  real(RP), parameter :: r2m_min=1.e-12_rp
956  real(RP), parameter :: um2cm = 100.0_rp
957 
958  real(RP) :: limitsw, zerosw
959  integer :: k, i, j
960  !---------------------------------------------------------------------------
961 
962  !$omp parallel do &
963  !$omp private(xc,xr,xi,xs,xg,dc_ave,dr_ave,rc,rr,ri2m,ri3m,rs2m,rs3m,rg2m,rg3m, &
964  !$omp limitsw,zerosw)
965  do j = js, je
966  do i = is, ie
967 
968  ! mean particle mass[kg]
969  do k = ks, ke
970  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) ) )
971  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) ) )
972  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) ) )
973  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) ) )
974  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) ) )
975  enddo
976 
977  ! diameter of average mass : SB06 eq.(32)
978  do k = ks, ke
979  dc_ave(k) = a_m(i_mp_qc) * xc(k)**b_m(i_mp_qc)
980  dr_ave(k) = a_m(i_mp_qr) * xr(k)**b_m(i_mp_qr)
981  enddo
982 
983  ! cloud effective radius
984  do k = ks, ke
985  rc = 0.5_rp * dc_ave(k)
986  limitsw = 0.5_rp + sign(0.5_rp, rc-rmin_re )
987  re(k,i,j,i_hc) = coef_re(i_mp_qc) * rc * limitsw * um2cm
988  enddo
989 
990  ! rain effective radius
991  do k = ks, ke
992  rr = 0.5_rp * dr_ave(k)
993  limitsw = 0.5_rp + sign(0.5_rp, rr-rmin_re )
994  re(k,i,j,i_hr) = coef_re(i_mp_qr) * rr * limitsw * um2cm
995  enddo
996 
997  do k = ks, ke
998  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)
999  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)
1000  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)
1001  enddo
1002 
1003  ! Fu(1996), eq.(3.11) or Fu et al.(1998), eq.(2.5)
1004  do k = ks, ke
1005  ri3m(k) = coef_fuetal1998 * qtrc0(k,i,j,i_ni) * xi(k)
1006  rs3m(k) = coef_fuetal1998 * qtrc0(k,i,j,i_ns) * xs(k)
1007  rg3m(k) = coef_fuetal1998 * qtrc0(k,i,j,i_ng) * xg(k)
1008  enddo
1009 
1010  ! ice effective radius
1011  do k = ks, ke
1012  zerosw = 0.5_rp - sign(0.5_rp, ri2m(k) - r2m_min )
1013  re(k,i,j,i_hi) = ri3m(k) / ( ri2m(k) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
1014  enddo
1015 
1016  ! snow effective radius
1017  do k = ks, ke
1018  zerosw = 0.5_rp - sign(0.5_rp, rs2m(k) - r2m_min )
1019  re(k,i,j,i_hs) = rs3m(k) / ( rs2m(k) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
1020  enddo
1021 
1022  ! graupel effective radius
1023  do k = ks, ke
1024  zerosw = 0.5_rp - sign(0.5_rp, rg2m(k) - r2m_min )
1025  re(k,i,j,i_hg) = rg3m(k) / ( rg2m(k) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
1026  enddo
1027 
1028  do k = ks, ke
1029  re(k,i,j,i_hh) = 0.0_rp
1030  end do
1031 
1032  enddo
1033  enddo
1034 
1035 
1036  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 1046 of file scale_atmos_phy_mp_sn14.F90.

1046  use scale_atmos_hydrometeor, only: &
1047  n_hyd, &
1048  i_hc, &
1049  i_hr, &
1050  i_hi, &
1051  i_hs, &
1052  i_hg, &
1053  i_hh
1054  implicit none
1055  integer, intent(in) :: KA, KS, KE
1056  integer, intent(in) :: IA, IS, IE
1057  integer, intent(in) :: JA, JS, JE
1058 
1059  real(RP), intent(in) :: QTRC0(KA,IA,JA,QA_MP-1) ! tracer mass concentration [kg/kg]
1060 
1061  real(RP), intent(out) :: Qe (KA,IA,JA,N_HYD) ! mixing ratio of each cateory [kg/kg]
1062 
1063  integer :: k, i, j
1064  !---------------------------------------------------------------------------
1065 
1066 !OCL XFILL
1067  !$omp parallel do
1068  do j = js, je
1069  do i = is, ie
1070  do k = ks, ke
1071  qe(k,i,j,i_hc) = qtrc0(k,i,j,i_mp_qc)
1072  qe(k,i,j,i_hr) = qtrc0(k,i,j,i_mp_qr)
1073  qe(k,i,j,i_hi) = qtrc0(k,i,j,i_mp_qi)
1074  qe(k,i,j,i_hs) = qtrc0(k,i,j,i_mp_qs)
1075  qe(k,i,j,i_hg) = qtrc0(k,i,j,i_mp_qg)
1076  qe(k,i,j,i_hh) = 0.0_rp
1077  end do
1078  end do
1079  end do
1080 
1081  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 1089 of file scale_atmos_phy_mp_sn14.F90.

1089  use scale_atmos_hydrometeor, only: &
1090  n_hyd, &
1091  i_hc, &
1092  i_hr, &
1093  i_hi, &
1094  i_hs, &
1095  i_hg, &
1096  i_hh
1097  implicit none
1098  integer, intent(in) :: KA, KS, KE
1099  integer, intent(in) :: IA, IS, IE
1100  integer, intent(in) :: JA, JS, JE
1101 
1102  real(RP), intent(in) :: QTRC0(KA,IA,JA,QA_MP-1) ! tracer mass concentration [kg/kg]
1103 
1104  real(RP), intent(out) :: Ne (KA,IA,JA,N_HYD) ! number density of each cateory [1/m3]
1105 
1106  integer :: k, i, j
1107  !---------------------------------------------------------------------------
1108 
1109 !OCL XFILL
1110  !$omp parallel do
1111  do j = js, je
1112  do i = is, ie
1113  do k = ks, ke
1114  ne(k,i,j,i_hc) = qtrc0(k,i,j,i_mp_nc)
1115  ne(k,i,j,i_hr) = qtrc0(k,i,j,i_mp_nr)
1116  ne(k,i,j,i_hi) = qtrc0(k,i,j,i_mp_ni)
1117  ne(k,i,j,i_hs) = qtrc0(k,i,j,i_mp_ns)
1118  ne(k,i,j,i_hg) = qtrc0(k,i,j,i_mp_ng)
1119  ne(k,i,j,i_hh) = 0.0_rp
1120  end do
1121  end do
1122  end do
1123 
1124  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 1132 of file scale_atmos_phy_mp_sn14.F90.

1132  use scale_const, only: &
1133  pi => const_pi, &
1134  undef => const_undef
1135  use scale_atmos_hydrometeor, only: &
1136  n_hyd, &
1137  i_hc, &
1138  i_hr, &
1139  i_hi, &
1140  i_hs, &
1141  i_hg, &
1142  i_hh
1143  implicit none
1144  integer, intent(in) :: KA, KS, KE
1145  integer, intent(in) :: IA, IS, IE
1146  integer, intent(in) :: JA, JS, JE
1147 
1148  real(RP), intent(in) :: Qe(KA,IA,JA,N_HYD) ! mass ratio of each cateory [kg/kg]
1149 
1150  real(RP), intent(out) :: QTRC(KA,IA,JA,QA_MP-1) ! tracer mass concentration [kg/kg]
1151 
1152  real(RP), intent(in), optional :: QNUM(KA,IA,JA,N_HYD)
1153 
1154  real(RP), parameter :: Dc = 20.e-6_rp ! typical particle diameter for cloud [m]
1155  real(RP), parameter :: Dr = 200.e-6_rp ! typical particle diameter for rain [m]
1156  real(RP), parameter :: Di = 80.e-6_rp ! typical particle diameter for ice [m]
1157  real(RP), parameter :: Ds = 80.e-6_rp ! typical particle diameter for snow [m]
1158  real(RP), parameter :: Dg = 200.e-6_rp ! typical particle diameter for grapel [m]
1159  real(RP), parameter :: b = 3.0_rp ! assume spherical form
1160 
1161  real(RP) :: piov6
1162 
1163  integer :: k, i, j
1164  !---------------------------------------------------------------------------
1165 
1166 
1167 !OCL XFILL
1168  !$omp parallel do
1169  do j = js, je
1170  do i = is, ie
1171  do k = ks, ke
1172  qtrc(k,i,j,i_mp_qc) = qe(k,i,j,i_hc)
1173  end do
1174  end do
1175  end do
1176 
1177 !OCL XFILL
1178  !$omp parallel do
1179  do j = js, je
1180  do i = is, ie
1181  do k = ks, ke
1182  qtrc(k,i,j,i_mp_qr) = qe(k,i,j,i_hr)
1183  end do
1184  end do
1185  end do
1186 
1187 !OCL XFILL
1188  !$omp parallel do
1189  do j = js, je
1190  do i = is, ie
1191  do k = ks, ke
1192  qtrc(k,i,j,i_mp_qi) = qe(k,i,j,i_hi)
1193  end do
1194  end do
1195  end do
1196 
1197 !OCL XFILL
1198  !$omp parallel do
1199  do j = js, je
1200  do i = is, ie
1201  do k = ks, ke
1202  qtrc(k,i,j,i_mp_qs) = qe(k,i,j,i_hs)
1203  end do
1204  end do
1205  end do
1206 
1207 !OCL XFILL
1208  !$omp parallel do
1209  do j = js, je
1210  do i = is, ie
1211  do k = ks, ke
1212  qtrc(k,i,j,i_mp_qg) = qe(k,i,j,i_hg) + qe(k,i,j,i_hh)
1213  end do
1214  end do
1215  end do
1216 
1217  piov6 = pi / 6.0_rp
1218 
1219  if ( present(qnum) ) then
1220 
1221 !OCL XFILL
1222  !$omp parallel do
1223  do j = js, je
1224  do i = is, ie
1225  do k = ks, ke
1226  if ( qnum(k,i,j,i_hc) .ne. undef ) then
1227  qtrc(k,i,j,i_mp_nc) = qnum(k,i,j,i_hc)
1228  else
1229  qtrc(k,i,j,i_mp_nc) = qtrc(k,i,j,i_mp_qc) / ( (piov6*rhow) * dc**b )
1230  end if
1231  end do
1232  end do
1233  end do
1234 
1235 !OCL XFILL
1236  !$omp parallel do
1237  do j = js, je
1238  do i = is, ie
1239  do k = ks, ke
1240  if ( qnum(k,i,j,i_hr) .ne. undef ) then
1241  qtrc(k,i,j,i_mp_nr) = qnum(k,i,j,i_hr)
1242  else
1243  qtrc(k,i,j,i_mp_nr) = qtrc(k,i,j,i_mp_qr) / ( (piov6*rhow) * dr**b )
1244  end if
1245  end do
1246  end do
1247  end do
1248 
1249 !OCL XFILL
1250  !$omp parallel do
1251  do j = js, je
1252  do i = is, ie
1253  do k = ks, ke
1254  if ( qnum(k,i,j,i_hi) .ne. undef ) then
1255  qtrc(k,i,j,i_mp_ni) = qnum(k,i,j,i_hi)
1256  else
1257  qtrc(k,i,j,i_mp_ni) = qtrc(k,i,j,i_mp_qi) / ( (piov6*rhof) * di**b )
1258  end if
1259  end do
1260  end do
1261  end do
1262 
1263 !OCL XFILL
1264  !$omp parallel do
1265  do j = js, je
1266  do i = is, ie
1267  do k = ks, ke
1268  if ( qnum(k,i,j,i_hs) .ne. undef ) then
1269  qtrc(k,i,j,i_mp_ns) = qnum(k,i,j,i_hs)
1270  else
1271  qtrc(k,i,j,i_mp_ns) = qtrc(k,i,j,i_mp_qs) / ( (piov6*rhof) * ds**b )
1272  end if
1273  end do
1274  end do
1275  end do
1276 
1277 !OCL XFILL
1278  !$omp parallel do
1279  do j = js, je
1280  do i = is, ie
1281  do k = ks, ke
1282  if ( qnum(k,i,j,i_hg) .ne. undef ) then
1283  if ( qnum(k,i,j,i_hh) .ne. undef ) then
1284  qtrc(k,i,j,i_mp_ng) = qnum(k,i,j,i_hg) + qnum(k,i,j,i_hh)
1285  else
1286  qtrc(k,i,j,i_mp_ng) = qnum(k,i,j,i_hg)
1287  end if
1288  else
1289  qtrc(k,i,j,i_mp_ng) = qtrc(k,i,j,i_mp_qg) / ( (piov6*rhog) * dg**b )
1290  end if
1291  end do
1292  end do
1293  end do
1294 
1295  else
1296 
1297 !OCL XFILL
1298  !$omp parallel do
1299  do j = js, je
1300  do i = is, ie
1301  do k = ks, ke
1302  qtrc(k,i,j,i_mp_nc) = qtrc(k,i,j,i_mp_qc) / ( (piov6*rhow) * dc**b )
1303  end do
1304  end do
1305  end do
1306 
1307 !OCL XFILL
1308  !$omp parallel do
1309  do j = js, je
1310  do i = is, ie
1311  do k = ks, ke
1312  qtrc(k,i,j,i_mp_nr) = qtrc(k,i,j,i_mp_qr) / ( (piov6*rhow) * dr**b )
1313  end do
1314  end do
1315  end do
1316 
1317 !OCL XFILL
1318  !$omp parallel do
1319  do j = js, je
1320  do i = is, ie
1321  do k = ks, ke
1322  qtrc(k,i,j,i_mp_ni) = qtrc(k,i,j,i_mp_qi) / ( (piov6*rhof) * di**b )
1323  end do
1324  end do
1325  end do
1326 
1327 !OCL XFILL
1328  !$omp parallel do
1329  do j = js, je
1330  do i = is, ie
1331  do k = ks, ke
1332  qtrc(k,i,j,i_mp_ns) = qtrc(k,i,j,i_mp_qs) / ( (piov6*rhof) * ds**b )
1333  end do
1334  end do
1335  end do
1336 
1337 !OCL XFILL
1338  !$omp parallel do
1339  do j = js, je
1340  do i = is, ie
1341  do k = ks, ke
1342  qtrc(k,i,j,i_mp_ng) = qtrc(k,i,j,i_mp_qg) / ( (piov6*rhog) * dg**b )
1343  end do
1344  end do
1345  end do
1346 
1347  end if
1348 
1349  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 1364 of file scale_atmos_phy_mp_sn14.F90.

1364  !$acc routine vector
1365  use scale_const, only: &
1366  const_undef
1367  implicit none
1368 
1369  integer, intent(in) :: KA, KS, KE
1370 
1371  real(RP), intent(in) :: RHOQ(KA,I_QC:I_NG) ! rho * q
1372  real(RP), intent(in) :: DENS(KA) ! rho
1373  real(RP), intent(in) :: TEMP(KA) ! temperature
1374  real(RP), intent(in) :: PRES(KA) ! pressure
1375 
1376  real(RP), intent(out) :: vterm(KA,QA_MP-1) ! terminal velocity of cloud mass
1377 
1378  real(RP) :: xq, log_xq ! average mass of 1 particle( mass/number )
1379 
1380  real(RP) :: rhofac ! density factor for terminal velocity
1381  real(RP) :: rhofac_q(KA), log_rhofac_q
1382 
1383  real(RP) :: rlambdar(KA) ! work for diagnosis of Rain DSD ( Seifert, 2008 )
1384  real(RP) :: mud_r
1385  real(RP) :: dq, log_dq ! weigthed diameter. Improved Rogers etal. (1993) formula by T.Mitsui
1386 
1387  real(RP) :: weight ! weighting coefficient for 2-branches is determined by ratio between 0.745mm and weighted diameter. SB06 Table.1
1388  real(RP) :: velq_s ! terminal velocity for small branch of Rogers formula
1389  real(RP) :: velq_l ! terminal velocity for large branch of Rogers formula
1390 
1391  real(RP) :: tmp
1392  integer :: k, i, j, iq
1393  !---------------------------------------------------------------------------
1394 
1395  ! QC, NC
1396  do k = ks, ke
1397  rhofac = rho_0 / max( dens(k), rho_min )
1398 
1399  log_rhofac_q = log(rhofac) * gamma_v(i_mp_qc)
1400  log_xq = log( max( xc_min, min( xc_max, rhoq(k,i_qc) / ( rhoq(k,i_nc) + nc_min ) ) ) )
1401 
1402  vterm(k,i_mp_qc) = - exp( log_rhofac_q + log_coef_vt1(i_mp_qc,1) + log_xq * beta_v(i_mp_qc,1) )
1403 
1404  vterm(k,i_mp_nc) = - exp( log_rhofac_q + log_coef_vt0(i_mp_qc,1) + log_xq * beta_vn(i_mp_qc,1) )
1405  end do
1406 
1407  ! QR, NR
1408  mud_r = 3.0_rp * nu(i_mp_qr) + 2.0_rp
1409  do k = ks, ke
1410  rhofac = rho_0 / max( dens(k), rho_min )
1411  rhofac_q(k) = rhofac**gamma_v(i_mp_qr)
1412  end do
1413  do k = ks, ke
1414  xq = max( xr_min, min( xr_max, rhoq(k,i_qr) / ( rhoq(k,i_nr) + nr_min ) ) )
1415 
1416  rlambdar(k) = a_m(i_mp_qr) * xq**b_m(i_mp_qr) &
1417  * ( (mud_r+3.0_rp) * (mud_r+2.0_rp) * (mud_r+1.0_rp) )**(-0.333333333_rp)
1418  end do
1419 !OCL LOOP_FISSION_TARGET(LS)
1420  do k = ks, ke
1421  dq = ( 4.0_rp + mud_r ) * rlambdar(k) ! D^(3)+mu weighted mean diameter
1422  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )
1423  velq_s = coef_vtr_ar2 * dq &
1424  * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar(k) )**(-5.0_rp-mud_r) )
1425  velq_l = coef_vtr_ar1 &
1426  - coef_vtr_br1 * ( 1.0_rp + coef_vtr_cr1*rlambdar(k) )**(-4.0_rp-mud_r)
1427  vterm(k,i_mp_qr) = -rhofac_q(k) * ( velq_l * ( weight ) &
1428  + velq_s * ( 1.0_rp - weight ) )
1429  end do
1430 !OCL LOOP_FISSION_TARGET(LS)
1431  do k = ks, ke
1432  dq = ( 1.0_rp + mud_r ) * rlambdar(k) ! D^(0)+mu weighted mean diameter
1433  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )
1434  velq_s = coef_vtr_ar2 * dq &
1435  * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar(k) )**(-2.0_rp-mud_r) )
1436  velq_l = coef_vtr_ar1 &
1437  - coef_vtr_br1 * ( 1.0_rp + coef_vtr_cr1*rlambdar(k) )**(-1.0_rp-mud_r)
1438  vterm(k,i_mp_nr) = -rhofac_q(k) * ( velq_l * ( weight ) &
1439  + velq_s * ( 1.0_rp - weight ) )
1440  end do
1441 
1442  do k = ks, ke
1443  rhofac_q(k) = exp( log( pres(k)/pre0_vt ) * a_pre0_vt + log( temp(k)/tem0_vt ) * a_tem0_vt )
1444  end do
1445 
1446  ! QI, NI
1447 !OCL LOOP_FISSION_TARGET(LS)
1448  do k = ks, ke
1449  log_xq = log( max( xi_min, min( xi_max, rhoq(k,i_qi) / ( rhoq(k,i_ni) + ni_min ) ) ) )
1450 
1451  tmp = log_a_m(i_mp_qi) + log_xq * b_m(i_mp_qi)
1452  log_dq = log_coef_dave_l(i_mp_qi) + tmp
1453  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log_dq - log_d0_li ) ) )
1454 
1455  velq_s = exp( log_coef_vt1(i_mp_qi,1) + log_xq * beta_v(i_mp_qi,1) )
1456  velq_l = exp( log_coef_vt1(i_mp_qi,2) + log_xq * beta_v(i_mp_qi,2) )
1457  vterm(k,i_mp_qi) = - rhofac_q(k) * ( velq_l * ( weight ) &
1458  + velq_s * ( 1.0_rp - weight ) )
1459 
1460  log_dq = log_coef_dave_n(i_mp_qi) + tmp
1461  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log_dq - log_d0_ni ) ) )
1462 
1463  velq_s = exp( log_coef_vt0(i_mp_qi,1) + log_xq * beta_vn(i_mp_qi,1) )
1464  velq_l = exp( log_coef_vt0(i_mp_qi,2) + log_xq * beta_vn(i_mp_qi,2) )
1465  vterm(k,i_mp_ni) = - rhofac_q(k) * ( velq_l * ( weight ) &
1466  + velq_s * ( 1.0_rp - weight ) )
1467  end do
1468 
1469  ! QS, NS
1470 !OCL LOOP_FISSION_TARGET(LS)
1471  do k = ks, ke
1472  log_xq = log( max( xs_min, min( xs_max, rhoq(k,i_qs) / ( rhoq(k,i_ns) + ns_min ) ) ) )
1473 
1474  tmp = log_a_m(i_mp_qs) + log_xq * b_m(i_mp_qs)
1475  log_dq = log_coef_dave_l(i_mp_qs) + tmp
1476  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log_dq - log_d0_ls ) ) )
1477 
1478  velq_s = exp( log_coef_vt1(i_mp_qs,1) + log_xq * beta_v(i_mp_qs,1) )
1479  velq_l = exp( log_coef_vt1(i_mp_qs,2) + log_xq * beta_v(i_mp_qs,2) )
1480  vterm(k,i_mp_qs) = - rhofac_q(k) * ( velq_l * ( weight ) &
1481  + velq_s * ( 1.0_rp - weight ) )
1482 
1483  log_dq = log_coef_dave_n(i_mp_qs) + tmp
1484  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log_dq - log_d0_ns ) ) )
1485 
1486  velq_s = exp( log_coef_vt0(i_mp_qs,1) + log_xq * beta_vn(i_mp_qs,1) )
1487  velq_l = exp( log_coef_vt0(i_mp_qs,2) + log_xq * beta_vn(i_mp_qs,2) )
1488  vterm(k,i_mp_ns) = - rhofac_q(k) * ( velq_l * ( weight ) &
1489  + velq_s * ( 1.0_rp - weight ) )
1490  end do
1491 
1492  ! QG, NG
1493 !OCL LOOP_FISSION_TARGET(LS)
1494  do k = ks, ke
1495  log_xq = log( max( xg_min, min( xg_max, rhoq(k,i_qg) / ( rhoq(k,i_ng) + ng_min ) ) ) )
1496 
1497  tmp = log_a_m(i_mp_qg) + log_xq * b_m(i_mp_qg)
1498  log_dq = log_coef_dave_l(i_mp_qg) + tmp
1499  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log_dq - log_d0_lg ) ) )
1500 
1501  velq_s = exp( log_coef_vt1(i_mp_qg,1) + log_xq * beta_v(i_mp_qg,1) )
1502  velq_l = exp( log_coef_vt1(i_mp_qg,2) + log_xq * beta_v(i_mp_qg,2) )
1503  vterm(k,i_mp_qg) = - rhofac_q(k) * ( velq_l * ( weight ) &
1504  + velq_s * ( 1.0_rp - weight ) )
1505 
1506  log_dq = log_coef_dave_n(i_mp_qg) + tmp
1507  weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log_dq - log_d0_ng ) ) )
1508 
1509  velq_s = exp( log_coef_vt0(i_mp_qg,1) + log_xq * beta_vn(i_mp_qg,1) )
1510  velq_l = exp( log_coef_vt0(i_mp_qg,2) + log_xq * beta_vn(i_mp_qg,2) )
1511  vterm(k,i_mp_ng) = - rhofac_q(k) * ( velq_l * ( weight ) &
1512  + velq_s * ( 1.0_rp - weight ) )
1513  enddo
1514 
1515  do iq = 1, qa_mp-1
1516  vterm(1:ks-2 ,iq) = const_undef
1517  vterm(ks-1 ,iq) = vterm(ks,iq)
1518  vterm(ke+1:ka,iq) = const_undef
1519  enddo
1520 
1521  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(), mixed_phase_collection_bin(), 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 3479 of file scale_atmos_phy_mp_sn14.F90.

3479  use scale_prc, only: &
3480  prc_myrank
3481  implicit none
3482  integer, intent(in) :: KA, KS, KE
3483 
3484  integer, intent(in) :: point
3485  integer, intent(in) :: i, j
3486  real(RP), intent(in) :: tem(KA)
3487  real(RP), intent(in) :: rho(KA)
3488  real(RP), intent(in) :: pre(KA)
3489  real(RP), intent(in) :: qv (KA)
3490 
3491  integer :: k
3492  !---------------------------------------------------------------------------
3493 
3494  do k = ks, ke
3495  if ( tem(k) < tem_min &
3496  .OR. rho(k) < rho_min &
3497  .OR. pre(k) < 1.0_rp ) then
3498 
3499  log_info("ATMOS_PHY_MP_SN14_debug_tem_kij",'(A,I3,A,4(F16.5),3(I6))') &
3500  "point: ", point, " low tem,rho,pre:", tem(k), rho(k), pre(k), qv(k), k, i, j, prc_myrank
3501  endif
3502  enddo
3503 
3504  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)  cva,
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,hydro_max)  dq_xa,
real(rp), dimension(ka,hydro_max,2)  vt_xa,
real(rp), dimension(ka,pq_max), intent(out)  PQ 
)

Definition at line 3519 of file scale_atmos_phy_mp_sn14.F90.

3519  use scale_prc, only: &
3520  prc_abort
3521  use scale_atmos_saturation, only: &
3522  moist_psat_liq => atmos_saturation_psat_liq, &
3523  moist_psat_ice => atmos_saturation_psat_ice, &
3524  moist_pres2qsat_liq => atmos_saturation_pres2qsat_liq, &
3525  moist_pres2qsat_ice => atmos_saturation_pres2qsat_ice, &
3526  moist_dqsi_dtem_dens => atmos_saturation_dqs_dtem_dens_liq, &
3527  moist_dqs_dtem_dpre_ice => atmos_saturation_dqs_dtem_dpre_ice
3528  use scale_atmos_hydrometeor, only: &
3529  cv_vapor, &
3530  cv_ice
3531  implicit none
3532 
3533  integer, intent(in) :: KA, KS, KE
3534 
3535  real(RP), intent(in) :: cz( KA)
3536  real(RP), intent(in) :: fz(0:KA)
3537  real(RP), intent(in) :: w (KA) ! w of full level
3538  real(RP), intent(in) :: rho (KA)
3539  real(RP), intent(in) :: tem (KA)
3540  real(RP), intent(in) :: pre (KA)
3541  real(RP), intent(in) :: qdry(KA)
3542  !
3543  real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
3544  !
3545  real(RP), intent(in) :: cpa(KA)
3546  real(RP), intent(in) :: dTdt_rad(KA)
3547  real(RP), intent(in) :: qke(KA)
3548  real(RP), intent(in) :: dt
3549  real(RP), intent(in) :: CCN(KA)
3550  real(RP), intent(in) :: nc_uplim_d
3551  !
3552  real(RP), intent(out) :: PQ(KA,PQ_MAX)
3553  !
3554  !
3555 ! real(RP) :: c_ccn_map(1) ! c_ccn horizontal distribution
3556 ! real(RP) :: kappa_map(1) ! kappa horizontal distribution
3557 ! real(RP) :: c_in_map(1) ! c_in horizontal distribution ! [Add] 11/08/30 T.Mitsui
3558  real(RP) :: esw(KA) ! saturation vapor pressure, water
3559  real(RP) :: esi(KA) ! ice
3560  real(RP) :: ssw(KA) ! super saturation (water)
3561  real(RP) :: ssi(KA) ! super saturation (ice)
3562 ! real(RP) :: w_dsswdz(KA) ! w*(d_ssw/ d_z) super saturation(water) flux
3563  real(RP) :: w_dssidz(KA) ! w*(d_ssi/ d_z), 09/04/14 T.Mitsui
3564 ! real(RP) :: ssw_below(KA) ! ssw(k-1)
3565  real(RP) :: ssi_below(KA) ! ssi(k-1), 09/04/14 T.Mitsui
3566  real(RP) :: z_below(KA) ! z(k-1)
3567  real(RP) :: dzh ! z(k)-z(k-1)
3568  real(RP) :: pv ! vapor pressure
3569  ! work variables for Twomey Equation.
3570  real(RP) :: qsw(KA)
3571  real(RP) :: qsi(KA)
3572  real(RP) :: dqsidtem_rho(KA)
3573  real(RP) :: dssidt_rad(KA)
3574  real(RP) :: wssi, wdssi
3575  !
3576 
3577  real(RP) :: cva(KA)
3578  real(RP) :: dq_xa(KA,HYDRO_MAX)
3579  real(RP) :: vt_xa(KA,HYDRO_MAX,2) ! terminal velocity
3580  real(RP) :: dTdt_dep(KA)
3581  real(RP) :: PLIdep_total(KA)
3582  real(RP) :: wtem(KA) ! temperature[K]
3583  real(RP) :: dqsidpre_tem(KA)
3584  real(RP) :: dqsidtem_pre(KA)
3585  real(RP) :: dssidt
3586  real(RP) :: dssidt_mp(KA)
3587  real(RP) :: dssidt_dyn(KA)
3588  !!
3589 
3590 
3591 ! real(RP) :: xi_nuc(1) ! xi use the value @ cloud base
3592 ! real(RP) :: alpha_nuc(1) ! alpha_nuc
3593 ! real(RP) :: eta_nuc(1) ! xi use the value @ cloud base
3594  !
3595  real(RP) :: sigma_w(KA)
3596  real(RP) :: weff(KA)
3597  real(RP) :: weff_max(KA)
3598  real(RP) :: velz(KA)
3599  !
3600  real(RP) :: coef_ccn
3601  real(RP) :: slope_ccn
3602  real(RP) :: nc_new(KA)
3603  real(RP) :: nc_new_below(KA)
3604  real(RP) :: dnc_new
3605  real(RP) :: nc_new_max ! Lohmann (2002),
3606  real(RP) :: a_max
3607  real(RP) :: b_max
3608  logical :: flag_nucleation(KA)
3609  !
3610  real(RP) :: r_gravity
3611  real(RP), parameter :: r_sqrt3=0.577350269_rp ! = sqrt(1.0/3.0)
3612  real(RP), parameter :: eps=1.e-30_rp
3613  !====> ! 09/08/18
3614  !
3615  real(RP) :: dlcdt_max, dli_max ! defined by supersaturation
3616  real(RP) :: dncdt_max, dni_max ! defined by supersaturation
3617  real(RP) :: rdt
3618 
3619  real(RP) :: tmp
3620  !
3621  integer :: k
3622  !
3623  !
3624 
3625  if( so22_het ) then
3626  do k = ks, ke
3627  ! Temperature lower limit is only used for saturation condition.
3628  ! On the other hand original "tem" is used for calculation of latent heat or energy equation.
3629  wtem(k) = max( tem(k), tem_min )
3630  end do
3631  endif
3632 
3633 
3634 ! c_ccn_map(1) = c_ccn
3635 ! kappa_map(1) = kappa
3636 ! c_in_map(1) = c_in
3637  !
3638  !
3639  rdt = 1.0_rp/dt
3640  r_gravity = 1.0_rp/grav
3641  !
3642  call moist_psat_liq ( ka, ks, ke, &
3643  tem(:), esw(:) ) ! [IN], [OUT]
3644  call moist_psat_ice ( ka, ks, ke, &
3645  tem(:), esi(:) ) ! [IN], [OUT]
3646  call moist_pres2qsat_liq ( ka, ks, ke, &
3647  tem(:), pre(:), qdry(:), & ! [IN]
3648  qsw(:) ) ! [OUT]
3649  call moist_pres2qsat_ice ( ka, ks, ke, &
3650  tem(:), pre(:), qdry(:), & ! [IN]
3651  qsi(:) ) ! [OUT]
3652  call moist_dqsi_dtem_dens( ka, ks, ke, &
3653  tem(:), rho(:), & ! [IN]
3654  dqsidtem_rho(:) ) ! [OUT]
3655 
3656  if( so22_het ) then
3657  call moist_dqs_dtem_dpre_ice( ka, ks, ke, &
3658  wtem(:), pre(:), qdry(:), & ! [IN]
3659  dqsidtem_pre(:), dqsidpre_tem(:) ) ! [OUT]
3660  endif
3661  !!
3662 
3663 
3664  !
3665  ! Lohmann (2002),JAS, eq.(1) but changing unit [cm-3] => [m-3]
3666  a_max = 1.e+6_rp*0.1_rp*(1.e-6_rp)**1.27_rp
3667  b_max = 1.27_rp
3668  !
3669  ssi_max = 1.0_rp
3670 
3671  do k = ks, ke
3672  pv = rhoq(k,i_qv)*rvap*tem(k)
3673  ssw(k) = min( mp_ssw_lim, ( pv/esw(k)-1.0_rp ) )*100.0_rp
3674  ssi(k) = ( pv/esi(k) - 1.00_rp )
3675 ! ssw_below(k+1) = ssw(k)
3676  ssi_below(k+1) = ssi(k)
3677  z_below(k+1) = cz(k)
3678  end do
3679 ! ssw_below(KS) = ssw(KS)
3680  ssi_below(ks) = ssi(ks)
3681  z_below(ks) = cz(ks-1)
3682 
3683  ! dS/dz is evaluated by first order upstream difference
3684  !*** Solution for Twomey Equation ***
3685 ! 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)) * &
3686  coef_ccn = 1.e+6_rp*0.88_rp*(c_ccn*1.e-6_rp)**(2.0_rp/(kappa + 2.0_rp)) &
3687 ! * (70.0_RP)**(kappa_map(1)/(kappa_map(1) + 2.0_RP))
3688  * (70.0_rp)**(kappa/(kappa + 2.0_rp))
3689 ! slope_ccn = 1.5_RP*kappa_map(1)/(kappa_map(1) + 2.0_RP)
3690  slope_ccn = 1.5_rp*kappa/(kappa + 2.0_rp)
3691  !
3692  do k=ks, ke
3693  sigma_w(k) = r_sqrt3*sqrt(max(qke(k),qke_min))
3694  end do
3695  sigma_w(ks-1) = sigma_w(ks)
3696  sigma_w(ke+1) = sigma_w(ke)
3697  ! effective vertical velocity
3698  do k=ks, ke
3699  weff(k) = w(k) - cpa(k)*r_gravity*dtdt_rad(k)
3700  end do
3701  !
3702  if( mp_couple_aerosol ) then
3703 
3704  do k = ks, ke
3705  if( ssw(k) > 1.e-10_rp .AND. pre(k) > 300.e+2_rp ) then
3706  nc_new(k) = max( ccn(k), c_ccn )
3707  else
3708  nc_new(k) = 0.0_rp
3709  endif
3710  enddo
3711 
3712  else
3713 
3714  if( nucl_twomey ) then
3715  ! diagnose cloud condensation nuclei
3716  do k = ks, ke
3717  ! effective vertical velocity (maximum vertical velocity in turbulent flow)
3718  weff_max(k) = weff(k) + sigma_w(k)
3719  ! large scale upward motion region and saturated
3720  if( (weff(k) > 1.e-8_rp) .AND. (ssw(k) > 1.e-10_rp) .AND. pre(k) > 300.e+2_rp )then
3721  ! Lohmann (2002), eq.(1)
3722  nc_new_max = coef_ccn*weff_max(k)**slope_ccn
3723  nc_new(k) = a_max*nc_new_max**b_max
3724  else
3725  nc_new(k) = 0.0_rp
3726  end if
3727  end do
3728 
3729  else
3730  ! calculate cloud condensation nuclei
3731  do k = ks, ke
3732  if( ssw(k) > 1.e-10_rp .AND. pre(k) > 300.e+2_rp ) then
3733  nc_new(k) = c_ccn*ssw(k)**kappa
3734  else
3735  nc_new(k) = 0.0_rp
3736  endif
3737  enddo
3738  endif
3739 
3740  endif
3741 
3742  do k = ks, ke
3743  ! nc_new is bound by upper limit
3744  if( nc_new(k) > nc_uplim_d ) then ! no more CCN
3745  flag_nucleation(k) = .false.
3746  nc_new_below(k+1) = 1.e+30_rp
3747  else if( nc_new(k) > eps ) then ! nucleation can occur
3748  flag_nucleation(k) = .true.
3749  nc_new_below(k+1) = nc_new(k)
3750  else ! nucleation cannot occur(unsaturated or negative w)
3751  flag_nucleation(k) = .false.
3752  nc_new_below(k+1) = 0.0_rp
3753  end if
3754  end do
3755  nc_new_below(ks) = 0.0_rp
3756 ! do k=KS, KE
3757  ! search maximum value of nc_new
3758 ! if( ( nc_new(k) < nc_new_below(k) ) .OR. &
3759 ! ( nc_new_below(k) > c_ccn_map(1)*0.05_RP ) )then ! 5% of c_ccn
3760 ! ( nc_new_below(k) > c_ccn*0.05_RP ) )then ! 5% of c_ccn
3761 ! flag_nucleation(k) = .false.
3762 ! end if
3763 ! end do
3764 
3765  if( mp_couple_aerosol ) then
3766  do k = ks, ke
3767  ! nucleation occurs at only cloud base.
3768  ! if CCN is more than below parcel, nucleation newly occurs
3769  ! effective vertical velocity
3770  if ( flag_nucleation(k) .AND. & ! large scale upward motion region and saturated
3771  tem(k) > tem_ccn_low ) then
3772  dlcdt_max = ( rhoq(k,i_qv) - esw(k) / ( rvap * tem(k) ) ) * rdt
3773  dlcdt_max = max( dlcdt_max, 0.0_rp ) ! dlcdt_max can be artificially negative due to truncation error in floating point operation
3774  dncdt_max = dlcdt_max/xc_min
3775 ! dnc_new = nc_new(k)-rhoq(k,I_NC)
3776  dnc_new = nc_new(k)
3777  pq(k,i_ncccn) = min( dncdt_max, dnc_new*rdt )
3778  pq(k,i_lcccn) = min( dlcdt_max, xc_min*pq(k,i_ncccn) )
3779  else
3780  pq(k,i_ncccn) = 0.0_rp
3781  pq(k,i_lcccn) = 0.0_rp
3782  end if
3783  end do
3784  else
3785 
3786  if( nucl_twomey ) then
3787  do k = ks, ke
3788  ! nucleation occurs at only cloud base.
3789  ! if CCN is more than below parcel, nucleation newly occurs
3790  ! effective vertical velocity
3791  if ( flag_nucleation(k) .AND. & ! large scale upward motion region and saturated
3792  tem(k) > tem_ccn_low .AND. &
3793  nc_new(k) > rhoq(k,i_nc) ) then
3794  dlcdt_max = ( rhoq(k,i_qv) - esw(k) / ( rvap * tem(k) ) ) * rdt
3795  dlcdt_max = max( dlcdt_max, 0.0_rp ) ! dlcdt_max can be artificially negative due to truncation error in floating point operation
3796  dncdt_max = dlcdt_max/xc_min
3797  dnc_new = nc_new(k)-rhoq(k,i_nc)
3798  pq(k,i_ncccn) = min( dncdt_max, dnc_new*rdt )
3799  pq(k,i_lcccn) = min( dlcdt_max, xc_min*pq(k,i_ncccn) )
3800  else
3801  pq(k,i_ncccn) = 0.0_rp
3802  pq(k,i_lcccn) = 0.0_rp
3803  end if
3804  end do
3805  else
3806  do k = ks, ke
3807  ! effective vertical velocity
3808  if( tem(k) > tem_ccn_low .AND. &
3809  nc_new(k) > rhoq(k,i_nc) ) then
3810  dlcdt_max = ( rhoq(k,i_qv) - esw(k) / ( rvap * tem(k) ) ) * rdt
3811  dlcdt_max = max( dlcdt_max, 0.0_rp ) ! dlcdt_max can be artificially negative due to truncation error in floating point operation
3812  dncdt_max = dlcdt_max/xc_min
3813  dnc_new = nc_new(k)-rhoq(k,i_nc)
3814  pq(k,i_ncccn) = min( dncdt_max, dnc_new*rdt )
3815  pq(k,i_lcccn) = min( dlcdt_max, xc_min*pq(k,i_ncccn) )
3816  else
3817  pq(k,i_ncccn) = 0.0_rp
3818  pq(k,i_lcccn) = 0.0_rp
3819  end if
3820  end do
3821  endif
3822 
3823  endif
3824 
3825  !
3826  ! ice nucleation
3827  !
3828  ! +++ NOTE ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3829  ! Based on Phillips etal.(2006).
3830  ! However this approach doesn't diagnose Ni itself but diagnose tendency.
3831  ! Original approach adjust Ni instantaneously .
3832  ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3833 
3834  if( so22_het .or. opt_nucleation_ice_hom ) then
3835  call dep_vapor_ice_wrk( & ! in
3836  ka, ks, ke, & ! in
3837  plidep_total(:), & ! out for ice nucleation
3838  rho(:), tem(:), pre(:), & ! in
3839  qdry(:), esi(:), qsi(:), & ! in
3840  rhoq(:,:), & ! in
3841  vt_xa, dq_xa, & ! in
3842  dt ) ! in
3843 
3844  do k = ks, ke
3845  dtdt_dep(k) = (lhs0+(cv_vapor-cv_ice)*tem(k))*plidep_total(k)/(rho(k)*cva(k))
3846  enddo
3847  else
3848  dtdt_dep(:) = 0.0_rp
3849  endif
3850 
3851  do k = ks, ke-1
3852  velz(k) = ( w(k) * ( cz(k+1) - fz(k) ) + w(k+1) * ( fz(k) - cz(k) ) ) / ( cz(k+1) - cz(k) ) ! @ half level
3853  end do
3854  velz(ke) = 0.0_rp
3855  do k = ks, ke
3856  dzh = cz(k) - z_below(k)
3857  w_dssidz(k) = velz(k) * (ssi(k) - ssi_below(k))/dzh
3858  dssidt_rad(k) = -rhoq(k,i_qv)/(rho(k)*qsi(k)*qsi(k))*dqsidtem_rho(k)*dtdt_rad(k)
3859  dli_max = ( rhoq(k,i_qv) - esi(k) / ( rvap * tem(k) ) ) * rdt
3860  dli_max = max( dli_max, 0.0_rp ) ! dli_max can be artificially negative due to truncation error in floating point operation
3861  dni_max = min( dli_max/xi_ccn, (in_max-rhoq(k,i_ni))*rdt )
3862  wdssi = min( w_dssidz(k)+dssidt_rad(k), 0.01_rp)
3863  wssi = min( ssi(k), ssi_max)
3864 
3865  !! Seiki and Ohno 2022
3866  if( so22_het ) then
3867  dssidt_mp(k) = -plidep_total(k)/(rho(k)*qsi(k))
3868  !! PLIdep(k) -> PQ(k,I_LIdep)
3869  dssidt_rad(k) = -rhoq(k,i_qv)/(rho(k)*qsi(k)*qsi(k))*dqsidtem_rho(k)*dtdt_rad(k)
3870  !!
3871  dssidt_dyn(k) = +rhoq(k,i_qv)/(rho(k)*qsi(k)*qsi(k))&
3872  * velz(k)*grav*(dqsidtem_pre(k)/cpa(k)+dqsidpre_tem(k)*rho(k))
3873  ! * w(k)*CNST_GRAV*(dqsidtem_pre(ij,k)/cpa(ij,k)+dqsidpre_tem(ij,k)*rho(ij,k))
3874  dssidt = dssidt_mp(k) + dssidt_rad(k) + dssidt_dyn(k)
3875  endif
3876 
3877  ! SB06(34),(35)
3878 !# if( ( wdssi > eps ) .AND. & !
3879 !# (tem(k) < 273.15_RP ) .AND. & !
3880 !# (rhoq(k,I_NI) < in_max ) .AND. &
3881 !# (wssi >= eps ) )then !
3882 !# tmp = c_in * nm_M92 * exp( 0.3_RP * bm_M92 * ( wssi - 0.1_RP ) )
3883 !# if( inucl_w ) then
3884 !# tmp = bm_M92 * 0.3_RP * tmp * wdssi
3885 !# else
3886 !# tmp = max( tmp - rhoq(k,I_NI), 0.0_RP ) * rdt
3887 !# endif
3888 !# PQ(k,I_NIccn) = min(dni_max, tmp)
3889 !# PQ(k,I_LIccn) = min(dli_max, PQ(k,I_NIccn)*xi_ccn )
3890 !# else
3891 !# PQ(k,I_NIccn) = 0.0_RP
3892 !# PQ(k,I_LIccn) = 0.0_RP
3893 !# end if
3894  if( (tem(k) < 273.15_rp ) .AND. & !
3895  (rhoq(k,i_ni) < in_max ) .AND. &
3896  (wssi >= eps ) )then !
3897  tmp = c_in * nm_m92 * exp( 0.3_rp * bm_m92 * ( wssi - 0.1_rp ) )
3898  if( inucl_w .and. wdssi > eps ) then
3899  tmp = bm_m92 * 0.3_rp * tmp * wdssi
3900  elseif( so22_het .and. dssidt > eps ) then
3901  tmp = bm_m92 * 0.3_rp * tmp * dssidt
3902  else
3903  tmp = max( tmp - rhoq(k,i_ni), 0.0_rp ) * rdt
3904  endif
3905  pq(k,i_niccn) = min(dni_max, tmp)
3906  pq(k,i_liccn) = min(dli_max, pq(k,i_niccn)*xi_ccn )
3907  else
3908  pq(k,i_niccn) = 0.0_rp
3909  pq(k,i_liccn) = 0.0_rp
3910  end if
3911 
3912 
3913  end do
3914 
3915  if( opt_nucleation_ice_hom ) then
3916  call nucleation_ice_hom( &
3917  ka, ks, ke, & !in
3918  tem, pre, rho, & !in
3919  qdry, rhoq(:,i_qv), & !in
3920  cva, cpa, & !in
3921  w, & !in
3922  dtdt_rad, & !in
3923  dtdt_dep, & !in
3924  plidep_total, dt, & !in
3925  pq(:,i_lihom), & !out
3926  pq(:,i_nihom) ) !out
3927  else
3928  pq(:,i_lihom) = 0.0_rp
3929  pq(:,i_nihom) = 0.0_rp
3930  endif
3931 
3932 
3933  return

References scale_atmos_hydrometeor::cv_ice, scale_atmos_hydrometeor::cv_vapor, dep_vapor_ice_wrk(), nucleation_ice_hom(), and 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:

◆ nucleation_ice_hom()

subroutine scale_atmos_phy_mp_sn14::nucleation_ice_hom ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
real(rp), dimension(ka), intent(in)  tem,
real(rp), dimension(ka), intent(in)  pre,
real(rp), dimension(ka), intent(in)  rho,
real(rp), dimension(ka), intent(in)  qd,
real(rp), dimension(ka), intent(in)  rhoq_qv,
real(rp), dimension(ka), intent(in)  cva,
real(rp), dimension(ka), intent(in)  cpa,
real(rp), dimension(ka), intent(in)  w,
real(rp), dimension(ka), intent(in)  dTdt_rad,
real(rp), dimension(ka), intent(in)  dTdt_dep,
real(rp), dimension(ka), intent(in)  PLIdep,
real(rp), intent(in)  dt,
real(rp), dimension(ka), intent(out)  PLIhom,
real(rp), dimension(ka), intent(out)  PNIhom 
)

Definition at line 3944 of file scale_atmos_phy_mp_sn14.F90.

3944  use scale_prc, only: &
3945  prc_abort
3946  use scale_const, only: &
3947  psat0 => const_psat0, &
3948  t00 => const_tem00, &
3949  pstd => const_pstd, &
3950  cpvap => const_cpvap, &
3951  cvvap => const_cvvap, &
3952  ci => const_ci, &
3953  cl => const_cl, &
3954  rvap => const_rvap, &
3955  rdry => const_rdry, &
3956  lhs00 => const_lhs00, &
3957  lhs0 => const_lhs0, &
3958  lhv00 => const_lhv00, &
3959  lhf00 => const_lhf00, &
3960  epsvap => const_epsvap, &
3961  grav => const_grav, &
3962  pi => const_pi
3963  implicit none
3964 
3965  integer, intent(in) :: KA, KS, KE
3966 
3967  real(RP), intent(in) :: tem(KA)
3968  real(RP), intent(in) :: pre(KA)
3969  real(RP), intent(in) :: rho(KA)
3970  real(RP), intent(in) :: qd(KA)
3971  real(RP), intent(in) :: rhoq_qv(KA)
3972  real(RP), intent(in) :: cpa(KA) ! specific heat @ cnst. pressure
3973  real(RP), intent(in) :: cva(KA) ! specific heat @ cnst. volume
3974  ! balance equation @ homogeneous ice nucleation
3975  real(RP), intent(in) :: w(KA) ! adiabatic ascending
3976  real(RP), intent(in) :: dTdt_rad(KA) ! effect of radiative heating
3977  real(RP), intent(in) :: dTdt_dep(KA) ! effect of radiative heating
3978  real(RP), intent(in) :: PLIdep(KA)
3979 ! real(RP), intent(out) :: PQ(KA,PQ_MAX) ! competition effect by vapor consumption
3980  real(RP), intent(out):: PLIhom(KA)
3981  real(RP), intent(out):: PNIhom(KA)
3982  real(RP), intent(in) :: dt
3983  real(RP), parameter :: rhoi=916.0_rp ! ice density [kg/m3]
3984  real(RP), parameter :: rrhoi=1.0_rp/rhoi !
3985  real(RP), parameter :: Mw=18.01528_rp ! Water Molar Weight
3986  real(RP), parameter :: Nav=6.0221415e+23_rp ! Avogadro Number
3987  real(RP), parameter :: vw=(mw*1.e-3_rp/nav)/rhoi ! volume of a water molecule in ice phase
3988  ! 18nm*exp( 3.0*log(1.5)**2 )=29.4760367
3989  real(RP), parameter :: r0=29.5e-9_rp ! aerosol mass(volume) mode radius
3990  real(RP), parameter :: c_gf = 1.01187_rp ! dAlmeida
3991  real(RP), parameter :: g_gf = -0.206449_rp ! dAlmeida
3992  real(RP), parameter :: rho_min=1.e-5_rp ! 3.e-3 is lower limit recognized in many experiments.
3993  real(RP), parameter :: tem_min=150.0_rp
3994  real(RP), parameter :: ni_max =300.e+6_rp
3995 ! real(RP) :: dTdt_dep(KA) ! effect of deposition heating
3996  real(RP) :: wtem
3997  real(RP) :: esi, esw
3998  real(RP) :: qsi
3999  real(RP) :: lhs
4000  real(RP) :: dqsidtem
4001  real(RP) :: den1, den2
4002  real(RP) :: dqsidt_pre
4003  real(RP) :: dqsidp_tem
4004  real(RP) :: rw
4005  real(RP) :: temc_lim
4006  real(RP) :: rho_lim
4007  real(RP) :: pre_lim
4008  real(RP) :: Dw
4009  real(RP) :: si, sw
4010  real(RP) :: Scr
4011  real(RP) :: dsidt_mp
4012  real(RP) :: dsidt_rd
4013  real(RP) :: wp
4014  real(RP) :: a1,a2,a3
4015  real(RP) :: b1,b2
4016  real(RP) :: dlogJdT
4017  real(RP) :: dtemdt_dyn
4018  real(RP) :: rtau ! 1/tau
4019  real(RP) :: delta
4020  real(RP) :: kappa
4021  real(RP) :: Rim_w
4022  real(RP) :: ri_wrk
4023  real(RP) :: ri
4024  !
4025  real(RP), parameter :: r2pi = 0.5_rp/pi ! 1/2pi
4026  real(RP), parameter :: sqrt_pi = sqrt(pi)
4027  real(RP), parameter :: coef_mi = 4.0_rp/3.0_rp*pi*rhoi
4028  real(RP), parameter :: eps = 1.e-30_rp
4029  real(RP) :: rdt
4030  integer :: ierr
4031  integer :: ij,k
4032  !
4033 
4034  rdt = 1.0_rp/dt
4035  plihom(:) = 0.0_rp
4036  pnihom(:) = 0.0_rp
4037 ! PQ(:,I_NIhom) = 0.0_RP
4038 ! PQ(:,I_LIhom) = 0.0_RP
4039 ! PQ(1:KS,I_NIhom) = 0.0_RP
4040 ! PQ(1:KS,I_LIhom) = 0.0_RP
4041 ! PQ(KE:KA,I_NIhom) = 0.0_RP
4042 ! PQ(KE:KA,I_LIhom) = 0.0_RP
4043  do k = ks, ke
4044  wtem= max(tem(k), tem_min)
4045  esi = min( psat0 &
4046  * ( wtem / t00 ) ** ( ( cpvap - ci ) / rvap ) &
4047  * exp( lhs00 / rvap &
4048  * ( 1.0_rp / t00 - 1.0_rp / wtem ) ), pre(k))
4049  esw = psat0 &
4050  * ( wtem / t00 ) ** ( ( cpvap - cl ) / rvap ) &
4051  * exp( lhv00 / rvap &
4052  * ( 1.0_rp / t00 - 1.0_rp / wtem ) )
4053  ! (10) in Ren and MacKenzie (2005)
4054  scr = 2.349_rp - wtem/259.0_rp
4055  qsi = epsvap * esi / ( pre(k) - ( 1.0_rp - epsvap ) * esi )
4056  si = rhoq_qv(k)*rvap*wtem/esi ! rho(k)*qv(k)
4057  ! sw must be less than 1 to calculate hygroscopic growth (otherwise, already activated)
4058  sw = min(rho(k)*rhoq_qv(k)*rvap*wtem/esw,0.999_rp)
4059  if ( si < scr ) then
4060  ! No nucleation occurs
4061  qsi = 0.0_rp
4062  lhs = lhs00 + (cpvap - ci )*(wtem-t00)
4063  dqsidtem = 0.0_rp
4064  si = 0.0_rp
4065  sw = 0.0_rp
4066  scr = 100.0_rp
4067  a1 = 0.0_rp
4068  wp = -1.0_rp
4069  rtau = -1.0_rp
4070  else
4071  ! (dsi/dt)_MP
4072  lhs = lhs0 + (cpvap - ci)*(wtem-t00)
4073  dqsidtem = esi/(rho(k)*rvap*wtem*wtem)&
4074  * (lhs/(rvap*wtem)-1.0_rp)
4075  dsidt_mp = -1.0_rp/(rho(k)*qsi) &
4076  * (1.0_rp+si*(lhv00+lhf00+(cvvap-ci)*wtem)/cva(k)*dqsidtem)&
4077  * plidep(k)
4078  dsidt_rd = -si/qsi*dqsidtem*dtdt_rad(k)
4079  ! This is an simple formulation from original paper
4080  a1 = lhs0*grav/(cpa(k)*rvap*tem(k)*tem(k))&
4081  - grav/(rdry*tem(k))
4082 !!$ ! alternative formulation for NICAM-EXACT
4083 !!$ den1 = (pre(k)-(1.0_RP-EPSvap)*esi)*(pre(k)-(1.0_RP-EPSvap)*esi)
4084 !!$ den2 = den1*Rvap*wtem*wtem
4085 !!$ dqsidp_tem=-EPSvap* esi/den1
4086 !!$ dqsidt_pre= EPSvap*pre(k)*lhs*esi/den2
4087 !!$ a1 = GRAV/qsi*(dqsidt_pre/cpa(k)+dqsidp_tem*rho(k))
4088  wp = w(k) + 1.0_rp/(a1*si)*(dsidt_mp+dsidt_rd)
4089  ! (21) in RM05
4090  dlogjdt = -0.004_rp*wtem*wtem + 2.0_rp*wtem - 304.4_rp
4091  ! (22) in RM05
4092  dtemdt_dyn= - grav*w(k)/cpa(k)
4093  rtau = dlogjdt*(dtdt_rad(k)+dtdt_dep(k)+dtemdt_dyn )
4094  endif
4095  rw = r0*c_gf*(1-sw)**g_gf
4096 
4097  if ( wp > eps .AND. rtau > eps ) then
4098  ! a1,a2,a3,b1,b2 are given by Appendix B in RM05
4099  temc_lim= max(tem(k)-t00, temc_lim_diff )
4100  rho_lim = max(rho(k),rho_min) !
4101  pre_lim = rho_lim*(qd(k)*rdry + rhoq_qv(k)*rvap)*(temc_lim+t00)! only for Dw
4102  dw = 0.211e-4_rp* (((temc_lim+t00)/t00)**1.94_rp) *(pstd/pre_lim)
4103  ! v_th = sqrt(8*Rv*T/pi) in the paragraph after eq(2) in Karcher etal.(2006)
4104  b2 = 0.5_rp/dw*sqrt(rvap*wtem*r2pi)
4105  b1 = (si-1.0_rp)*0.5_rp*rrhoi*esi/sqrt(2.0_rp*pi*rvap*wtem)
4106  a2 = mw*rvap*wtem/(nav*esi)
4107  a3 = epsvap*mw*lhs0*lhs0/(nav*cpa(k)*pre(k)*wtem)
4108  delta = b2*rw
4109  kappa = 2.0_rp*b1*b2/( rtau*(1.0_rp+delta)*(1.0_rp+delta) )
4110  ! (A9) in RM05
4111  rim_w = max(1.e-20_rp, 0.5_rp*(1.0_rp+delta)*(3.0_rp*kappa/(2.0_rp+sqrt(1.0_rp+9.0_rp/pi*kappa))) &
4112  + 1.0_rp/(1.0_rp+delta)*(3.0_rp /(2.0_rp+sqrt(1.0_rp+9.0_rp/pi*kappa)))+delta-1.0_rp )
4113  pnihom(k) = min( scr/(scr-1.0_rp)*a1*wp/(rim_w*4.0_rp*pi*dw/b2), ni_max )* rdt
4114  ri_wrk = 1.0_rp+b2*rw
4115  ri = ( sqrt(ri_wrk*ri_wrk + 2.0_rp*b1*b2*dt )-1.0_rp )/b2
4116  plihom(k) = coef_mi*ri*ri*ri*pnihom(k)
4117  else
4118  plihom(k) = 0.0_rp
4119  pnihom(k) = 0.0_rp
4120  endif
4121  enddo
4122  !
4123  return

References scale_const::const_ci, scale_const::const_cl, scale_const::const_cpvap, scale_const::const_cvvap, scale_const::const_epsvap, scale_const::const_grav, scale_const::const_lhf00, scale_const::const_lhs0, scale_const::const_lhs00, scale_const::const_lhv00, scale_const::const_pi, scale_const::const_psat0, scale_const::const_pstd, scale_const::const_rdry, scale_const::const_rvap, scale_const::const_tem00, and scale_prc::prc_abort().

Referenced by nucleation().

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(ka,i_qc:i_qg), 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(ka,pq_max), intent(inout)  Pcrg1 
)

Definition at line 4134 of file scale_atmos_phy_mp_sn14.F90.

4134 
4135  ! ice multiplication by splintering
4136  ! we consider Hallet-Mossop process
4137  use scale_specfunc, only: &
4138  gammafunc => sf_gamma
4139  implicit none
4140  integer, intent(in) :: KA, KS, KE
4141 
4142  real(RP), intent(in) :: Pac(KA,Pac_MAX)
4143  real(RP), intent(in) :: tem(KA)
4144  real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
4145  real(RP), intent(in) :: xq(KA,HYDRO_MAX)
4146  !
4147  real(RP), intent(inout):: PQ(KA,PQ_MAX)
4148  ! for lightning
4149  logical, intent(in) :: flg_lt
4150  real(RP), intent(in) :: rhoq_crg(KA,I_QC:I_QG)
4151  real(RP), intent(inout):: Pcrg1(KA,PQ_MAX)
4152  !
4153  ! production of (350.d3 ice particles)/(cloud riming [g]) => 350*d6 [/kg]
4154  real(RP), parameter :: pice = 350.0e+6_rp
4155  ! production of (1 ice particle)/(250 cloud particles riming)
4156  real(RP), parameter :: pnc = 250.0_rp
4157  ! temperature factor
4158  real(RP) :: fp
4159 
4160  real(RP) :: igm ! in complete gamma(x,alpha)
4161  real(RP) :: x
4162  ! coefficient of expansion using in calculation of igm
4163  real(RP) :: a0,a1,a2,a3,a4,a5
4164  real(RP) :: a6,a7,a8,a9,a10
4165  real(RP) :: an1,an2,b0,b1,b2,c0,c1,c2
4166  real(RP) :: d0,d1,d2,e1,e2,h0,h1,h2
4167  real(RP), parameter :: eps=1.0e-30_rp
4168  ! number of cloud droplets larger than 12 micron(radius).
4169  real(RP) :: n12
4170  !
4171  real(RP) :: wn, wni, wns, wng
4172  integer :: k, iq
4173  real(RP) :: sw1
4174  !
4175  !
4176 !OCL LOOP_FISSION_TARGET(LS)
4177  do k = ks, ke
4178  ! Here we assume particle temperature is same as environment temperature.
4179  ! If you want to treat in a better manner,
4180  ! you can diagnose with eq.(64) in CT(86)
4181  if (tem(k) > 270.16_rp)then
4182  fp = 0.0_rp
4183  else if(tem(k) >= 268.16_rp)then
4184  fp = (270.16_rp-tem(k))*0.5_rp
4185  else if(tem(k) >= 265.16_rp)then
4186  fp = (tem(k)-265.16_rp)*0.333333333_rp
4187  else
4188  fp = 0.0_rp
4189  end if
4190  ! Approximation of Incomplete Gamma function
4191  ! Here we calculate with algorithm by Numerical Recipes.
4192  ! This approach is based on eq.(78) in Cotton etal.(1986),
4193  ! but more accurate and expanded for Generalized Gamma Distribution.
4194  x = coef_lambda(i_mp_qc)*(xc_cr/xq(k,i_mp_qc))**mu(i_mp_qc)
4195  !
4196  if(x<1.e-2_rp*alpha)then ! negligible
4197  igm = 0.0_rp
4198  else if(x<alpha+1.0_rp)then ! series expansion
4199  ! 10th-truncation is enough for cloud droplet.
4200  a0 = 1.0_rp/alpha ! n=0
4201  a1 = a0*x/(alpha+1.0_rp) ! n=1
4202  a2 = a1*x/(alpha+2.0_rp) ! n=2
4203  a3 = a2*x/(alpha+3.0_rp) ! n=3
4204  a4 = a3*x/(alpha+4.0_rp) ! n=4
4205  a5 = a4*x/(alpha+5.0_rp) ! n=5
4206  a6 = a5*x/(alpha+6.0_rp) ! n=6
4207  a7 = a6*x/(alpha+7.0_rp) ! n=7
4208  a8 = a7*x/(alpha+8.0_rp) ! n=8
4209  a9 = a8*x/(alpha+9.0_rp) ! n=9
4210  a10 = a9*x/(alpha+10.0_rp) ! n=10
4211  igm = (a0+a1+a2+a3+a4+a5+a6+a7+a8+a9+a10)*exp( -x + alpha*log(x) - lgm )
4212  else if(x<alpha*100.0_rp) then ! continued fraction expansion
4213  ! 2nd-truncation is enough for cloud droplet.
4214  ! setup
4215  b0 = x+1.0_rp-alpha
4216  c0 = 1.0_rp/eps
4217  d0 = 1.0_rp/b0
4218  h0 = d0
4219  ! n=1
4220  an1 = -(1.0_rp-alpha)
4221  b1 = b0 + 2.0_rp
4222  d1 = 1.0_rp/(an1*d0+b1)
4223  c1 = b1+an1/c0
4224  e1 = d1*c1
4225  h1 = h0*e1
4226  ! n=2
4227  an2 = -2.0_rp*(2.0_rp-alpha)
4228  b2 = b1 + 2.0_rp
4229  d2 = 1.0_rp/(an2*d1+b2)
4230  c2 = b2+an2/c1
4231  e2 = d2*c2
4232  h2 = h1*e2
4233  !
4234  igm = 1.0_rp - exp( -x + alpha*log(x) - lgm )*h2
4235  else ! negligible
4236  igm = 1.0_rp
4237  end if
4238  ! n12 is number of cloud droplets larger than 12 micron.
4239  n12 = rhoq(k,i_nc)*(1.0_rp-igm)
4240  ! eq.(82) CT(86)
4241  wn = (pice + n12/((rhoq(k,i_qc)+xc_min)*pnc) )*fp ! filtered by xc_min
4242  wni = wn*(-pac(k,i_liaclc2li) ) ! riming production rate is all negative
4243  wns = wn*(-pac(k,i_lsaclc2ls) )
4244  wng = wn*(-pac(k,i_lgaclc2lg) )
4245  pq(k,i_nispl) = wni+wns+wng
4246  !
4247  pq(k,i_lsspl) = - wns*xq(k,i_mp_qi) ! snow => ice
4248  pq(k,i_lgspl) = - wng*xq(k,i_mp_qi) ! graupel => ice
4249  if (flg_lt) then
4250  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NC is small,ignore charge transfer
4251  pcrg1(k,i_nsspl) = - wns*(1.0_rp-sw1) &
4252  / (rhoq(k,i_ns)+sw1)*rhoq_crg(k,i_qs)
4253  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small ) !--- if NC is small,ignore charge transfer
4254  pcrg1(k,i_ngspl) = - wng*(1.0_rp-sw1) &
4255  / (rhoq(k,i_ng)+sw1)*rhoq_crg(k,i_qg)
4256  pcrg1(k,i_nispl) = - ( pcrg1(k,i_nsspl) + pcrg1(k,i_ngspl) )
4257  end if
4258  !
4259  end do
4260  !
4261  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(ka,i_qc:i_qg), 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(ka,pq_max), intent(inout)  Pcrg1,
real(rp), dimension(ka,pcrg_max), intent(inout)  Pcrg2,
real(rp), dimension(ka,pac_max), intent(out)  Pac 
)

Definition at line 4278 of file scale_atmos_phy_mp_sn14.F90.

4278  use scale_atmos_saturation, only: &
4279  moist_psat_ice => atmos_saturation_psat_ice
4280  implicit none
4281 
4282  integer, intent(in) :: KA, KS, KE
4283 
4284  !--- mixed-phase collection process
4285  ! And all we set all production term as a negative sign to avoid confusion.
4286  !
4287  real(RP), intent(in) :: wtem(KA)
4288  !--- mass/number concentration[kg/m3]
4289  real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
4290  ! necessary ?
4291  real(RP), intent(in) :: xq(KA,HYDRO_MAX)
4292  !--- diameter of averaged mass( D(ave x) )
4293  real(RP), intent(in) :: dq_xave(KA,HYDRO_MAX)
4294  !--- terminal velocity of averaged mass( vt(ave x) )
4295  real(RP), intent(in) :: vt_xave(KA,HYDRO_MAX,2)
4296  ! [Add] 11/08/30 T.Mitsui, for autoconversion of ice
4297  ! real(RP), intent(in) :: rho(KA)
4298  !--- partial conversion
4299  real(RP), intent(inout):: PQ(KA,PQ_MAX)
4300  !
4301  real(RP), intent(out):: Pac(KA,Pac_MAX)
4302  !--- for lightning component
4303  logical, intent(in) :: flg_lt
4304  real(RP), intent(in) :: beta_crg(KA)
4305  real(RP), intent(in) :: dqcrg(KA)
4306  real(RP), intent(in) :: d0_crg, v0_crg
4307  real(RP), intent(in) :: rhoq_crg(KA,I_QC:I_QG)
4308  real(RP), intent(inout):: Pcrg1(KA,PQ_MAX)
4309  real(RP), intent(inout):: Pcrg2(KA,Pcrg_MAX)
4310 
4311  real(RP), parameter :: a_dec = 0.883_rp
4312  real(RP), parameter :: b_dec = 0.093_rp
4313  real(RP), parameter :: c_dec = 0.00348_rp
4314  real(RP), parameter :: d_dec = 4.5185e-5_rp
4315  !
4316  !
4317  real(RP) :: tem(KA)
4318  !
4319  !--- collection efficency of each specie
4320  real(RP) :: E_c(KA), E_r, E_i, E_s, E_g
4321  real(RP) :: E_ic, E_sc, E_gc
4322  !--- sticking efficiency
4323  real(RP) :: E_stick(KA)
4324  ! [Add] 10/08/03 T.Mitsui
4325  real(RP) :: temc, temc2, temc3
4326  real(RP) :: E_dec
4327  real(RP) :: esi_rat
4328  real(RP) :: esi(KA)
4329  !
4330  real(RP) :: temc_p, temc_m ! celcius tem.
4331 ! real(RP) :: ci_aut(KA)
4332 ! real(RP) :: taui_aut(KA)
4333 ! real(RP) :: tau_sce(KA)
4334  !--- DSD averaged diameter for each species
4335  real(RP) :: ave_dc ! cloud
4336 ! real(RP) :: ave_dr ! rain
4337  real(RP) :: ave_di(KA) ! ice
4338  real(RP) :: ave_ds(KA) ! snow
4339  real(RP) :: ave_dg ! graupel
4340  !--- coefficient of collection equations(L:mass, N:number)
4341  real(RP) :: coef_acc_LCI, coef_acc_NCI ! cloud - cloud ice
4342  real(RP) :: coef_acc_LCS, coef_acc_NCS ! cloud - snow
4343  !
4344  real(RP) :: coef_acc_LCG, coef_acc_NCG ! cloud - graupel
4345  real(RP) :: coef_acc_LRI_I, coef_acc_NRI_I ! rain - cloud ice
4346  real(RP) :: coef_acc_LRI_R, coef_acc_NRI_R ! rain - cloud ice
4347  real(RP) :: coef_acc_LRS_S, coef_acc_NRS_S ! rain - snow
4348  real(RP) :: coef_acc_LRS_R, coef_acc_NRS_R ! rain - snow
4349  real(RP) :: coef_acc_LRG, coef_acc_NRG ! rain - graupel
4350  real(RP) :: coef_acc_LII, coef_acc_NII ! cloud ice - cloud ice
4351  real(RP) :: coef_acc_LIS, coef_acc_NIS ! cloud ice - snow
4352  real(RP) :: coef_acc_NSS ! snow - snow
4353  real(RP) :: coef_acc_NGG ! grauepl - graupel
4354  real(RP) :: coef_acc_LSG, coef_acc_NSG(KA) ! snow - graupel
4355  !--- (diameter) x (diameter)
4356  real(RP) :: dcdc(KA), dcdi, dcds, dcdg
4357  real(RP) :: drdr(KA), drdi(KA), drds(KA), drdg
4358  real(RP) :: didi(KA), dids, didg
4359  real(RP) :: dsds(KA), dsdg
4360  real(RP) :: dgdg(KA)
4361  !--- (terminal velocity) x (terminal velocity)
4362  real(RP) :: vcvc(KA), vcvi, vcvs, vcvg
4363  real(RP) :: vrvr(KA), vrvi(KA), vrvs(KA), vrvg
4364  real(RP) :: vivi(KA), vivs, vivg
4365  real(RP) :: vsvs(KA), vsvg
4366  real(RP) :: vgvg(KA)
4367  !
4368  real(RP) :: wx_cri, wx_crs
4369  real(RP) :: coef_emelt
4370  real(RP) :: w1
4371 
4372  real(RP) :: sw, sw1, sw2
4373  real(RP) :: alpha_lt
4374  !
4375  integer :: k, iqw
4376  !
4377  !
4378  do k = ks, ke
4379  tem(k) = max( wtem(k), tem_min ) ! 11/08/30 T.Mitsui
4380  end do
4381 
4382  call moist_psat_ice( ka, ks, ke, &
4383  tem(:), esi(:) ) ! [IN], [OUT]
4384 
4385  if( opt_stick_ks96 )then
4386  do k = ks, ke
4387  ! Khain and Sednev (1996), eq.(3.15)
4388  temc = tem(k) - t00
4389  temc2 = temc*temc
4390  temc3 = temc2*temc
4391  e_dec = max(0.0_rp, a_dec + b_dec*temc + c_dec*temc2 + d_dec*temc3 )
4392  esi_rat = rhoq(k,i_qv)*rvap*tem(k)/esi(k)
4393  e_stick(k) = min(1.0_rp, e_dec*esi_rat)
4394  end do
4395  else if( opt_stick_co86 )then
4396  do k = ks, ke
4397  ! [Add] 11/08/30 T.Mitsui, Cotton et al. (1986)
4398  temc = min(tem(k) - t00,0.0_rp)
4399  w1 = 0.035_rp*temc-0.7_rp
4400  e_stick(k) = 10._rp**w1
4401  end do
4402  else
4403  do k = ks, ke
4404  ! Lin et al. (1983)
4405  temc_m = min(tem(k) - t00,0.0_rp) ! T < 273.15
4406  e_stick(k) = exp(0.09_rp*temc_m)
4407  end do
4408  end if
4409 
4410  do k = ks, ke
4411  ! averaged diameter using SB06(82)
4412  ave_dc = coef_d(i_mp_qc)*xq(k,i_mp_qc)**b_m(i_mp_qc)
4413  !------------------------------------------------------------------------
4414  ! coellection efficiency are given as follows
4415  e_c(k) = max(0.0_rp, min(1.0_rp, (ave_dc-dc0)/(dc1-dc0) ))
4416  end do
4417 
4418  !------------------------------------------------------------------------
4419  ! Collection: a collects b ( assuming particle size a>b )
4420  do k = ks, ke
4421  dcdc(k) = dq_xave(k,i_mp_qc) * dq_xave(k,i_mp_qc)
4422  drdr(k) = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qr)
4423  didi(k) = dq_xave(k,i_mp_qi) * dq_xave(k,i_mp_qi)
4424  dsds(k) = dq_xave(k,i_mp_qs) * dq_xave(k,i_mp_qs)
4425  dgdg(k) = dq_xave(k,i_mp_qg) * dq_xave(k,i_mp_qg)
4426  drdi(k) = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qi)
4427  drds(k) = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qs)
4428  end do
4429  do k = ks, ke
4430  vcvc(k) = vt_xave(k,i_mp_qc,2) * vt_xave(k,i_mp_qc,2)
4431  vrvr(k) = vt_xave(k,i_mp_qr,2) * vt_xave(k,i_mp_qr,2)
4432  vivi(k) = vt_xave(k,i_mp_qi,2) * vt_xave(k,i_mp_qi,2)
4433  vsvs(k) = vt_xave(k,i_mp_qs,2) * vt_xave(k,i_mp_qs,2)
4434  vgvg(k) = vt_xave(k,i_mp_qg,2) * vt_xave(k,i_mp_qg,2)
4435  vrvi(k) = vt_xave(k,i_mp_qr,2) * vt_xave(k,i_mp_qi,2)
4436  vrvs(k) = vt_xave(k,i_mp_qr,2) * vt_xave(k,i_mp_qs,2)
4437  end do
4438 
4439  do k = ks, ke
4440  ave_di(k) = coef_d(i_mp_qi)*xq(k,i_mp_qi)**b_m(i_mp_qi)
4441  ave_ds(k) = coef_d(i_mp_qs)*xq(k,i_mp_qs)**b_m(i_mp_qs)
4442  end do
4443 
4444  !------------------------------------------------------------------------
4445  !
4446  !+++ pattern 1: a + b => a (a>b)
4447  ! (i-c, s-c, g-c, s-i, g-r, s-g)
4448  !------------------------------------------------------------------------
4449 
4450  ! cloud-ice => ice
4451  ! reduction term of cloud
4452  do k = ks, ke
4453  dcdi = dq_xave(k,i_mp_qc) * dq_xave(k,i_mp_qi)
4454  vcvi = vt_xave(k,i_mp_qc,2) * vt_xave(k,i_mp_qi,2)
4455  sw = 0.5_rp - sign(0.5_rp, di0-ave_di(k)) ! if(ave_di>di0)then sw=1
4456  e_i = e_im * sw
4457  e_ic = e_i*e_c(k)
4458  coef_acc_lci = &
4459  ( delta_b1(i_mp_qc)*dcdc(k) + delta_ab1(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi(k) ) &
4460  * sqrt( theta_b1(i_mp_qc)*vcvc(k) - theta_ab1(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi(k) &
4461  + sigma_i + sigma_c )
4462  coef_acc_nci = &
4463  ( delta_b0(i_mp_qc)*dcdc(k) + delta_ab0(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi(k) ) &
4464  * sqrt( theta_b0(i_mp_qc)*vcvc(k) - theta_ab0(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi(k) &
4465  + sigma_i + sigma_c )
4466  pac(k,i_liaclc2li)= -0.25_rp*pi*e_ic*rhoq(k,i_ni)*rhoq(k,i_qc)*coef_acc_lci
4467  pac(k,i_niacnc2ni)= -0.25_rp*pi*e_ic*rhoq(k,i_ni)*rhoq(k,i_nc)*coef_acc_nci
4468  end do
4469 
4470  ! cloud-snow => snow
4471  ! reduction term of cloud
4472  do k = ks, ke
4473  dcds = dq_xave(k,i_mp_qc) * dq_xave(k,i_mp_qs)
4474  vcvs = vt_xave(k,i_mp_qc,2) * vt_xave(k,i_mp_qs,2)
4475  sw = 0.5_rp - sign(0.5_rp, ds0-ave_ds(k)) ! if(ave_ds>ds0)then sw=1
4476  e_s = e_sm * sw
4477  e_sc = e_s*e_c(k)
4478  coef_acc_lcs = &
4479  ( delta_b1(i_mp_qc)*dcdc(k) + delta_ab1(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds(k) ) &
4480  * sqrt( theta_b1(i_mp_qc)*vcvc(k) - theta_ab1(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs(k) &
4481  + sigma_s + sigma_c )
4482  coef_acc_ncs = &
4483  ( delta_b0(i_mp_qc)*dcdc(k) + delta_ab0(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds(k) ) &
4484  * sqrt( theta_b0(i_mp_qc)*vcvc(k) - theta_ab0(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs(k) &
4485  + sigma_s + sigma_c )
4486  pac(k,i_lsaclc2ls)= -0.25_rp*pi*e_sc*rhoq(k,i_ns)*rhoq(k,i_qc)*coef_acc_lcs
4487  pac(k,i_nsacnc2ns)= -0.25_rp*pi*e_sc*rhoq(k,i_ns)*rhoq(k,i_nc)*coef_acc_ncs
4488  end do
4489 
4490  ! cloud-graupel => graupel
4491  ! reduction term of cloud
4492  do k = ks, ke
4493  dcdg = dq_xave(k,i_mp_qc) * dq_xave(k,i_mp_qg)
4494  vcvg = vt_xave(k,i_mp_qc,2) * vt_xave(k,i_mp_qg,2)
4495  ave_dg = coef_d(i_mp_qg)*xq(k,i_mp_qg)**b_m(i_mp_qg)
4496  sw = 0.5_rp - sign(0.5_rp, dg0-ave_dg) ! if(ave_dg>dg0)then sw=1
4497  e_g = e_gm * sw
4498  e_gc = e_g*e_c(k)
4499  coef_acc_lcg = &
4500  ( delta_b1(i_mp_qc)*dcdc(k) + delta_ab1(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg(k) ) &
4501  * sqrt( theta_b1(i_mp_qc)*vcvc(k) - theta_ab1(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg(k) &
4502  + sigma_g + sigma_c )
4503  coef_acc_ncg = &
4504  ( delta_b0(i_mp_qc)*dcdc(k) + delta_ab0(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg(k) ) &
4505  * sqrt( theta_b0(i_mp_qc)*vcvc(k) - theta_ab0(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg(k) &
4506  + sigma_g + sigma_c )
4507  pac(k,i_lgaclc2lg)= -0.25_rp*pi*e_gc*rhoq(k,i_ng)*rhoq(k,i_qc)*coef_acc_lcg
4508  pac(k,i_ngacnc2ng)= -0.25_rp*pi*e_gc*rhoq(k,i_ng)*rhoq(k,i_nc)*coef_acc_ncg
4509  end do
4510 
4511  ! snow-graupel => graupel
4512  do k = ks, ke
4513  dsdg = dq_xave(k,i_mp_qs) * dq_xave(k,i_mp_qg)
4514  vsvg = vt_xave(k,i_mp_qs,2) * vt_xave(k,i_mp_qg,2)
4515  coef_acc_lsg = &
4516  ( delta_b1(i_mp_qs)*dsds(k) + delta_ab1(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg(k) ) &
4517  * sqrt( theta_b1(i_mp_qs)*vsvs(k) - theta_ab1(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg(k) &
4518  + sigma_g + sigma_s )
4519  coef_acc_nsg(k) = &
4520  ( delta_b0(i_mp_qs)*dsds(k) + delta_ab0(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg(k) ) &
4521  ! [fix] T.Mitsui 08/05/08
4522  * sqrt( theta_b0(i_mp_qs)*vsvs(k) - theta_ab0(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg(k) &
4523  + sigma_g + sigma_s )
4524  pac(k,i_lgacls2lg)= -0.25_rp*pi*e_stick(k)*e_gs*rhoq(k,i_ng)*rhoq(k,i_qs)*coef_acc_lsg
4525  pac(k,i_ngacns2ng)= -0.25_rp*pi*e_stick(k)*e_gs*rhoq(k,i_ng)*rhoq(k,i_ns)*coef_acc_nsg(k)
4526  end do
4527 
4528  !-----------------
4529  ! (start) Y.Sato added on 2018/8/31
4530  !--- ice-graupel => graupel
4531 !!$ do k = KS, KE
4532 !!$ didg = dq_xave(k,I_mp_QI) * dq_xave(k,I_mp_QG)
4533 !!$ vivg = vt_xave(k,I_mp_QI,2) * vt_xave(k,I_mp_QG,2)
4534 !!$ coef_acc_LIG = &
4535 !!$ ( delta_b1(I_QI)*didi(k) + delta_ab1(I_QG,I_QI)*didg + delta_b0(I_QG)*dgdg(k) ) &
4536 !!$ * sqrt( theta_b1(I_QI)*vivi(k) - theta_ab1(I_QG,I_QI)*vivg + theta_b0(I_QG)*vgvg(k) &
4537 !!$ + sigma_g + sigma_i )
4538 !!$ coef_acc_NIG(k) = &
4539 !!$ ( delta_b0(I_QI)*didi(k) + delta_ab0(I_QG,I_QI)*didg + delta_b0(I_QG)*dgdg(k) ) &
4540 !!$ * sqrt( theta_b0(I_QI)*vivi(k) - theta_ab0(I_QG,I_QI)*vivg + theta_b0(I_QG)*vgvg(k) &
4541 !!$ + sigma_g + sigma_i )
4542 !!$ 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
4543 !!$ Pac(k,I_NGacNI2NG)= -0.25_RP*pi*E_stick(k)*E_gi*rhoq(k,I_NG)*rhoq(k,I_NI)*coef_acc_NIG(k)*flg_igcol
4544 !!$ end do
4545 
4546  !------------------------------------------------------------------------
4547  ! ice-snow => snow
4548  ! reduction term of ice
4549  do k = ks, ke
4550  dids = dq_xave(k,i_mp_qi) * dq_xave(k,i_mp_qs)
4551  vivs = vt_xave(k,i_mp_qi,2) * vt_xave(k,i_mp_qs,2)
4552  coef_acc_lis = &
4553  ( delta_b1(i_mp_qi)*didi(k) + delta_ab1(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds(k) ) &
4554  * sqrt( theta_b1(i_mp_qi)*vivi(k) - theta_ab1(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs(k) &
4555  + sigma_i + sigma_s )
4556  coef_acc_nis = &
4557  ( delta_b0(i_mp_qi)*didi(k) + delta_ab0(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds(k) ) &
4558  * sqrt( theta_b0(i_mp_qi)*vivi(k) - theta_ab0(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs(k) &
4559  + sigma_i + sigma_s )
4560  pac(k,i_liacls2ls)= -0.25_rp*pi*e_stick(k)*e_si*rhoq(k,i_ns)*rhoq(k,i_qi)*coef_acc_lis
4561  pac(k,i_niacns2ns)= -0.25_rp*pi*e_stick(k)*e_si*rhoq(k,i_ns)*rhoq(k,i_ni)*coef_acc_nis
4562  end do
4563 
4564  do k = ks, ke
4565  drdg = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qg)
4566  vrvg = vt_xave(k,i_mp_qr,2) * vt_xave(k,i_mp_qg,2)
4567  sw = sign(0.5_rp, t00-tem(k)) + 0.5_rp
4568  ! if ( tem(k) <= T00 )then
4569  ! rain-graupel => graupel
4570  ! reduction term of rain
4571  ! sw = 1
4572  ! else
4573  ! rain-graupel => rain
4574  ! reduction term of graupel
4575  ! sw = 0
4576  coef_acc_lrg = &
4577  ( ( delta_b1(i_mp_qr)*drdr(k) + delta_ab1(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg(k) ) * sw &
4578  + ( delta_b1(i_mp_qg)*dgdg(k) + delta_ab1(i_mp_qr,i_mp_qg)*drdg + delta_b0(i_mp_qr)*drdr(k) ) * (1.0_rp-sw) ) &
4579  * sqrt( ( theta_b1(i_mp_qr)*vrvr(k) - theta_ab1(i_mp_qg,i_mp_qr)*vrvg + theta_b0(i_mp_qg)*vgvg(k) ) * sw &
4580  + ( theta_b1(i_mp_qg)*vgvg(k) - theta_ab1(i_mp_qr,i_mp_qg)*vrvg + theta_b0(i_mp_qr)*vrvr(k) ) * (1.0_rp-sw) &
4581  + sigma_r + sigma_g )
4582  pac(k,i_lraclg2lg) = -0.25_rp*pi*e_gr*coef_acc_lrg * ( rhoq(k,i_ng)*rhoq(k,i_qr) * sw &
4583  + rhoq(k,i_nr)*rhoq(k,i_qg) * (1.0_rp-sw) )
4584  coef_acc_nrg = &
4585  ( delta_b0(i_mp_qr)*drdr(k) + delta_ab0(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg(k) ) &
4586  * sqrt( theta_b0(i_mp_qr)*vrvr(k) - theta_ab0(i_mp_qg,i_mp_qr)*vrvg + theta_b0(i_mp_qg)*vgvg(k) &
4587  + sigma_r + sigma_g )
4588  pac(k,i_nracng2ng) = -0.25_rp*pi*e_gr*rhoq(k,i_ng)*rhoq(k,i_nr)*coef_acc_nrg
4589  end do
4590 
4591  !------------------------------------------------------------------------
4592  !
4593  !+++ pattern 2: a + b => c (a>b)
4594  ! (r-i,r-s)
4595  !------------------------------------------------------------------------
4596 
4597  ! rain-ice => graupel
4598  ! reduction term of ice
4599  do k = ks, ke
4600  coef_acc_lri_i = &
4601  ( delta_b1(i_mp_qi)*didi(k) + delta_ab1(i_mp_qr,i_mp_qi)*drdi(k) + delta_b0(i_mp_qr)*drdr(k) ) &
4602  * sqrt( theta_b1(i_mp_qi)*vivi(k) - theta_ab1(i_mp_qr,i_mp_qi)*vrvi(k) + theta_b0(i_mp_qr)*vrvr(k) &
4603  + sigma_r + sigma_i )
4604  coef_acc_nri_i = &
4605  ( delta_b0(i_mp_qi)*didi(k) + delta_ab0(i_mp_qr,i_mp_qi)*drdi(k) + delta_b0(i_mp_qr)*drdr(k) ) &
4606  * sqrt( theta_b0(i_mp_qi)*vivi(k) - theta_ab0(i_mp_qr,i_mp_qi)*vrvi(k) + theta_b0(i_mp_qr)*vrvr(k) &
4607  + sigma_r + sigma_i )
4608  pac(k,i_lracli2lg_i)= -0.25_rp*pi*e_ir*rhoq(k,i_nr)*rhoq(k,i_qi)*coef_acc_lri_i
4609  pac(k,i_nracni2ng_i)= -0.25_rp*pi*e_ir*rhoq(k,i_nr)*rhoq(k,i_ni)*coef_acc_nri_i
4610  end do
4611 
4612  ! reduction term of rain
4613  do k = ks, ke
4614  coef_acc_lri_r = &
4615  ( delta_b1(i_mp_qr)*drdr(k) + delta_ab1(i_mp_qi,i_mp_qr)*drdi(k) + delta_b0(i_mp_qi)*didi(k) ) &
4616  * sqrt( theta_b1(i_mp_qr)*vrvr(k) - theta_ab1(i_mp_qi,i_mp_qr)*vrvi(k) + theta_b0(i_mp_qi)*vivi(k) &
4617  + sigma_r + sigma_i )
4618  coef_acc_nri_r = &
4619  ( delta_b0(i_mp_qr)*drdr(k) + delta_ab0(i_mp_qi,i_mp_qr)*drdi(k) + delta_b0(i_mp_qi)*didi(k) ) &
4620  * sqrt( theta_b0(i_mp_qr)*vrvr(k) - theta_ab0(i_mp_qi,i_mp_qr)*vrvi(k) + theta_b0(i_mp_qi)*vivi(k) &
4621  + sigma_r + sigma_i )
4622  pac(k,i_lracli2lg_r)= -0.25_rp*pi*e_ir*rhoq(k,i_ni)*rhoq(k,i_qr)*coef_acc_lri_r
4623  pac(k,i_nracni2ng_r)= -0.25_rp*pi*e_ir*rhoq(k,i_ni)*rhoq(k,i_nr)*coef_acc_nri_r
4624  end do
4625 
4626  ! rain-snow => graupel
4627  ! reduction term of snow
4628  do k = ks, ke
4629  coef_acc_lrs_s = &
4630  ( delta_b1(i_mp_qs)*dsds(k) + delta_ab1(i_mp_qr,i_mp_qs)*drds(k) + delta_b0(i_mp_qr)*drdr(k) ) &
4631  * sqrt( theta_b1(i_mp_qs)*vsvs(k) - theta_ab1(i_mp_qr,i_mp_qs)*vrvs(k) + theta_b0(i_mp_qr)*vrvr(k) &
4632  + sigma_r + sigma_s )
4633  coef_acc_nrs_s = &
4634  ( delta_b0(i_mp_qs)*dsds(k) + delta_ab0(i_mp_qr,i_mp_qs)*drds(k) + delta_b0(i_mp_qr)*drdr(k) ) &
4635  * sqrt( theta_b0(i_mp_qs)*vsvs(k) - theta_ab0(i_mp_qr,i_mp_qs)*vrvs(k) + theta_b0(i_mp_qr)*vrvr(k) &
4636  + sigma_r + sigma_s )
4637  pac(k,i_lracls2lg_s)= -0.25_rp*pi*e_sr*rhoq(k,i_nr)*rhoq(k,i_qs)*coef_acc_lrs_s
4638  pac(k,i_nracns2ng_s)= -0.25_rp*pi*e_sr*rhoq(k,i_nr)*rhoq(k,i_ns)*coef_acc_nrs_s
4639  end do
4640 
4641  ! reduction term of rain
4642  do k = ks, ke
4643  coef_acc_lrs_r = &
4644  ( delta_b1(i_mp_qr)*drdr(k) + delta_ab1(i_mp_qs,i_mp_qr)*drds(k) + delta_b0(i_mp_qs)*dsds(k) ) &
4645  * sqrt( theta_b1(i_mp_qr)*vrvr(k) - theta_ab1(i_mp_qs,i_mp_qr)*vrvs(k) + theta_b0(i_mp_qs)*vsvs(k) &
4646  + sigma_r + sigma_s )
4647  coef_acc_nrs_r = &
4648  ( delta_b0(i_mp_qr)*drdr(k) + delta_ab0(i_mp_qs,i_mp_qr)*drds(k) + delta_b0(i_mp_qs)*dsds(k) ) &
4649  * sqrt( theta_b0(i_mp_qr)*vrvr(k) - theta_ab0(i_mp_qs,i_mp_qr)*vrvs(k) + theta_b0(i_mp_qs)*vsvs(k) &
4650  + sigma_r + sigma_s )
4651  pac(k,i_lracls2lg_r)= -0.25_rp*pi*e_sr*rhoq(k,i_ns)*rhoq(k,i_qr)*coef_acc_lrs_r
4652  pac(k,i_nracns2ng_r)= -0.25_rp*pi*e_sr*rhoq(k,i_ns)*rhoq(k,i_nr)*coef_acc_nrs_r
4653  end do
4654 
4655  !------------------------------------------------------------------------
4656  !
4657  !+++ pattern 3: a + a => b (i-i)
4658  !
4659  !------------------------------------------------------------------------
4660 
4661  ! ice-ice ( reduction is double, but includes double-count)
4662  do k = ks, ke
4663  coef_acc_lii = &
4664  ( delta_b0(i_mp_qi)*didi(k) + delta_ab1(i_mp_qi,i_mp_qi)*didi(k) + delta_b1(i_mp_qi)*didi(k) ) &
4665  * sqrt( theta_b0(i_mp_qi)*vivi(k) - theta_ab1(i_mp_qi,i_mp_qi)*vivi(k) + theta_b1(i_mp_qi)*vivi(k) &
4666  + sigma_i + sigma_i )
4667  coef_acc_nii = &
4668  ( delta_b0(i_mp_qi)*didi(k) + delta_ab0(i_mp_qi,i_mp_qi)*didi(k) + delta_b0(i_mp_qi)*didi(k) ) &
4669  * sqrt( theta_b0(i_mp_qi)*vivi(k) - theta_ab0(i_mp_qi,i_mp_qi)*vivi(k) + theta_b0(i_mp_qi)*vivi(k) &
4670  + sigma_i + sigma_i )
4671  pac(k,i_liacli2ls)= -0.25_rp*pi*e_stick(k)*e_ii*rhoq(k,i_ni)*rhoq(k,i_qi)*coef_acc_lii
4672  pac(k,i_niacni2ns)= -0.25_rp*pi*e_stick(k)*e_ii*rhoq(k,i_ni)*rhoq(k,i_ni)*coef_acc_nii
4673 
4674 ! ci_aut(k) = 0.25_RP*pi*E_ii*rhoq(k,I_NI)*coef_acc_LII
4675 ! taui_aut(k) = 1._RP/max(E_stick(k)*ci_aut(k),1.E-10_RP)
4676 ! tau_sce(k) = rhoq(k,I_QI)/max(rhoq(k,I_QIj)+rhoq(k,I_QS),1.E-10_RP)
4677  end do
4678 
4679  !------------------------------------------------------------------------
4680  !
4681  !+++ pattern 4: a + a => a (s-s)
4682  !
4683  !------------------------------------------------------------------------
4684 
4685  ! snow-snow => snow
4686  do k = ks, ke
4687  coef_acc_nss = &
4688  ( delta_b0(i_mp_qs)*dsds(k) + delta_ab0(i_mp_qs,i_mp_qs)*dsds(k) + delta_b0(i_mp_qs)*dsds(k) ) &
4689  * sqrt( theta_b0(i_mp_qs)*vsvs(k) - theta_ab0(i_mp_qs,i_mp_qs)*vsvs(k) + theta_b0(i_mp_qs)*vsvs(k) &
4690  + sigma_s + sigma_s )
4691  pac(k,i_nsacns2ns)= -0.125_rp*pi*e_stick(k)*e_ss*rhoq(k,i_ns)*rhoq(k,i_ns)*coef_acc_nss
4692  end do
4693 
4694  ! graupel-grauple => graupel
4695  do k = ks, ke
4696  coef_acc_ngg = &
4697  ( delta_b0(i_mp_qg)*dgdg(k) + delta_ab0(i_mp_qg,i_mp_qg)*dgdg(k) + delta_b0(i_mp_qg)*dgdg(k) ) &
4698  * sqrt( theta_b0(i_mp_qg)*vgvg(k) - theta_ab0(i_mp_qg,i_mp_qg)*vgvg(k) + theta_b0(i_mp_qg)*vgvg(k) &
4699  + sigma_g + sigma_g )
4700  pac(k,i_ngacng2ng)= -0.125_rp*pi*e_stick(k)*e_gg*rhoq(k,i_ng)*rhoq(k,i_ng)*coef_acc_ngg
4701  end do
4702 
4703  !------------------------------------------------------------------------
4704  !--- Partial conversion
4705  ! SB06(70),(71)
4706  ! i_iconv2g: option whether partial conversions work or not
4707 
4708  ! ice-cloud => graupel
4709  do k = ks, ke
4710  sw = 0.5_rp - sign(0.5_rp,di_cri-ave_di(k)) ! if( ave_di > di_cri )then sw=1
4711  wx_cri = cfill_i*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_di(k)**3/xq(k,i_mp_qi) - 1.0_rp ) * sw
4712  pq(k,i_licon) = i_iconv2g * pac(k,i_liaclc2li)/max(1.0_rp, wx_cri) * sw
4713  pq(k,i_nicon) = i_iconv2g * pq(k,i_licon)/xq(k,i_mp_qi) * sw
4714  end do
4715 
4716  ! snow-cloud => graupel
4717  do k = ks, ke
4718  wx_crs = cfill_s*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_ds(k)**3/xq(k,i_mp_qs) - 1.0_rp )
4719  pq(k,i_lscon) = i_sconv2g * (pac(k,i_lsaclc2ls))/max(1.0_rp, wx_crs)
4720  pq(k,i_nscon) = i_sconv2g * pq(k,i_lscon)/xq(k,i_mp_qs)
4721  end do
4722 
4723  !--- enhanced melting( due to collection-freezing of water droplets )
4724  ! originally from Rutledge and Hobbs(1984). eq.(A.21)
4725  ! if T > 273.15 then temc_p=T-273.15, else temc_p=0
4726  ! 08/05/08 [fix] T.Mitsui LHF00 => LHF0
4727  ! melting occurs around T=273K, so LHF0 is suitable both SIMPLE and EXACT,
4728  ! otherwise LHF can have sign both negative(EXACT) and positive(SIMPLE).
4729  do k = ks, ke
4730 ! temc_m = min(tem(k) - T00, 0.0_RP) ! T < 273.15
4731  temc_p = max(tem(k) - t00, 0.0_rp) ! T > 273.15
4732 !!$ coef_emelt = -CL/LHF00*temc_p
4733  coef_emelt = cl/lhf0*temc_p
4734  ! cloud-graupel
4735  pq(k,i_lgacm) = coef_emelt*pac(k,i_lgaclc2lg)
4736  pq(k,i_ngacm) = pq(k,i_lgacm)/xq(k,i_mp_qg)
4737  ! rain-graupel
4738  pq(k,i_lgarm) = coef_emelt*pac(k,i_lraclg2lg)
4739  pq(k,i_ngarm) = pq(k,i_lgarm)/xq(k,i_mp_qg)
4740  ! cloud-snow
4741  pq(k,i_lsacm) = coef_emelt*(pac(k,i_lsaclc2ls))
4742  pq(k,i_nsacm) = pq(k,i_lsacm)/xq(k,i_mp_qs)
4743  ! rain-snow
4744  pq(k,i_lsarm) = coef_emelt*(pac(k,i_lracls2lg_r)+pac(k,i_lracls2lg_s))
4745  pq(k,i_nsarm) = pq(k,i_lsarm)/xq(k,i_mp_qg)
4746  ! cloud-ice
4747  pq(k,i_liacm) = coef_emelt*pac(k,i_liaclc2li)
4748  pq(k,i_niacm) = pq(k,i_liacm)/xq(k,i_mp_qi)
4749  ! rain-ice
4750  pq(k,i_liarm) = coef_emelt*(pac(k,i_lracli2lg_r)+pac(k,i_lracli2lg_i))
4751  pq(k,i_niarm) = pq(k,i_liarm)/xq(k,i_mp_qg)
4752  end do
4753 
4754 
4755  !---- for charge density
4756  if ( flg_lt ) then
4757  !--- C + I -> I (decrease from cloud charge)
4758  do k = ks, ke
4759  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small ) !--- if NC is small, ignore charge transfer
4760  pcrg2(k,i_niacnc2ni) = pac(k,i_niacnc2ni)*(1.0_rp-sw1) / (rhoq(k,i_nc)+sw1) * rhoq_crg(k,i_qc)
4761  end do
4762 
4763  !--- C + S -> S (decrease from cloud charge)
4764  do k = ks, ke
4765  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small ) !--- if NC is small, ignore charge transfer
4766  pcrg2(k,i_nsacnc2ns) = pac(k,i_nsacnc2ns)*(1.0_rp-sw1) / (rhoq(k,i_nc)+sw1) * rhoq_crg(k,i_qc)
4767  end do
4768 
4769  !--- C + G -> G (decrease from cloud charge)
4770  do k = ks, ke
4771  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small ) !--- if NC is small, ignore charge transfer
4772  pcrg2(k,i_ngacnc2ng) = pac(k,i_ngacnc2ng)*(1.0_rp-sw1) / (rhoq(k,i_nc)+sw1) * rhoq_crg(k,i_qc)
4773  end do
4774 
4775  !--- S + G -> G (decrease from snow charge)
4776  do k = ks, ke
4777  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NS is small, ignore charge transfer
4778  pcrg2(k,i_ngacns2ng) = pac(k,i_ngacns2ng)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
4779  end do
4780 
4781  !--- Charge separation by Snow-Graupel rebound--------------------------------
4782  do k = ks, ke
4783  alpha_lt = 5.0_rp * ( dq_xave(k,i_mp_qs) / d0_crg )**2 * vt_xave(k,i_mp_qg,2) / v0_crg
4784  alpha_lt = min( alpha_lt, 10.0_rp )
4785  pcrg2(k,i_cgngacns2ng)= 0.25_rp*pi*( 1.0_rp - e_stick(k) )*e_gs &
4786  * rhoq(k,i_ng)*rhoq(k,i_ns)*coef_acc_nsg(k) &
4787  * ( dqcrg(k)*alpha_lt ) &
4788  * beta_crg(k)
4789  end do
4790 
4791  !--- I + G -> G (decrease from snow charge)
4792 !!$ do k = KS, KE
4793 !!$ sw1 = 0.5_RP - sign( 0.5_RP, rhoq(k,I_NI)-SMALL ) !--- if NS is small, ignore charge transfer
4794 !!$ Pcrg2(k,I_NGacNI2NG) = Pac(k,I_NGacNI2NG)*(1.0_RP-sw1) / (rhoq(k,I_NI)+sw1) * rhoq_crg(k,I_QI) * flg_igcol
4795 !!$ !--- Charge separation by Ice-Graupel rebound--------------------------------
4796 !!$ alpha_lt = 5.0_RP * ( dq_xave(k,I_mp_QI) / d0_crg )**2 * vt_xave(k,I_mp_QG,2) / v0_crg
4797 !!$ alpha_lt = min( alpha_lt, 10.0_RP )
4798 !!$ Pcrg2(k,I_CGNGacNI2NG)= 0.25_RP*pi*( 1.0_RP - E_stick(k) )*E_gi &
4799 !!$ * rhoq(k,I_NG)*rhoq(k,I_NI)*coef_acc_NIG(k) &
4800 !!$ * ( dqcrg(k)*alpha_lt ) &
4801 !!$ * beta_crg(k) * flg_igcol
4802 !!$ end do
4803 
4804  !--- I + S -> S (decrease from ice charge)
4805  do k = ks, ke
4806  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small, ignore charge transfer
4807  pcrg2(k,i_niacns2ns) = pac(k,i_niacns2ns)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
4808  end do
4809 
4810  !--- R+G->R (T>T00 sw=0, decrease from graupel charge), ->G(T<=T00 sw=1, dcrerase from rain charge)
4811  do k = ks, ke
4812  sw = 0.5_rp + sign( 0.5_rp, t00-tem(k) )
4813  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small ) !--- if NR is small, ignore charge transfer
4814  sw2 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small ) !--- if NG is small, ignore charge transfer
4815  pcrg2(k,i_nracng2ng) = pac(k,i_nracng2ng)*(1.0_rp-sw1)/(rhoq(k,i_nr)+sw1) * rhoq_crg(k,i_qr) * sw &
4816  + pac(k,i_nracng2ng)*(1.0_rp-sw2)/(rhoq(k,i_ng)+sw2) * rhoq_crg(k,i_qg) * (1.0_rp-sw)
4817  end do
4818 
4819  !--- R + I -> G (decrease from both ice and rain charge, but only ice charge at this part)
4820  do k = ks, ke
4821  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small, ignore charge transfer
4822  pcrg2(k,i_nracni2ng_i) = pac(k,i_nracni2ng_i)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
4823  end do
4824 
4825  !--- R + I -> G (decrease from both ice and rain charge, but only rain charge at this part)
4826  do k = ks, ke
4827  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small ) !--- if NR is small, ignore charge transfer
4828  pcrg2(k,i_nracni2ng_r) = pac(k,i_nracni2ng_r)*(1.0_rp-sw1) / (rhoq(k,i_nr)+sw1) * rhoq_crg(k,i_qr)
4829  end do
4830 
4831  !--- R + S -> G (decrease from both snow and rain charge, but only snow charge at this part)
4832  do k = ks, ke
4833  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NS is small, ignore charge transfer
4834  pcrg2(k,i_nracns2ng_s) = pac(k,i_nracns2ng_s)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
4835  end do
4836 
4837  !--- R + S -> G (decrease from both snow and rain charge, but only rain charge at this part)
4838  do k = ks, ke
4839  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small ) !--- if NR is small, ignore charge transfer
4840  pcrg2(k,i_nracns2ng_r) = pac(k,i_nracns2ng_r)*(1.0_rp-sw1) / (rhoq(k,i_nr)+sw1) * rhoq_crg(k,i_qr)
4841  end do
4842 
4843  !--- I + I -> S (decrease from ice charge)
4844  do k = ks, ke
4845  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small, ignore charge transfer
4846  pcrg2(k,i_niacni2ns) = pac(k,i_niacni2ns)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
4847  end do
4848 
4849  !--- I + C -> G (decrease from ice charge)
4850  do k = ks, ke
4851  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small, ignore charge transfer
4852  pcrg1(k,i_nicon) = i_iconv2g * pq(k,i_nicon)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
4853  end do
4854 
4855  !--- S + C -> G (decrease from snow charge)
4856  do k = ks, ke
4857  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NS is small, ignore charge transfer
4858  pcrg1(k,i_nscon) = i_sconv2g * pq(k,i_nscon)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
4859  end do
4860 
4861  do k = ks, ke
4862  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small ) !--- if NG is small, ignore charge transfer
4863  pcrg1(k,i_ngacm) = pq(k,i_ngacm)*(1.0_rp-sw1) / (rhoq(k,i_ng)+sw1) * rhoq_crg(k,i_qg)
4864  end do
4865 
4866  do k = ks, ke
4867  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small ) !--- if NG is small, ignore charge transfer
4868  pcrg1(k,i_ngarm) = pq(k,i_ngarm)*(1.0_rp-sw1) / (rhoq(k,i_ng)+sw1) * rhoq_crg(k,i_qg)
4869  end do
4870 
4871  do k = ks, ke
4872  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NS is small, ignore charge transfer
4873  pcrg1(k,i_nsacm) = pq(k,i_nsacm)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
4874  end do
4875 
4876  do k = ks, ke
4877  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NS is small, ignore charge transfer
4878  pcrg1(k,i_nsarm) = pq(k,i_nsarm)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
4879  end do
4880 
4881  do k = ks, ke
4882  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small, ignore charge transfer
4883  pcrg1(k,i_niacm) = pq(k,i_niacm)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
4884  end do
4885 
4886  do k = ks, ke
4887  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small, ignore charge transfer
4888  pcrg1(k,i_niarm) = pq(k,i_niarm)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
4889  end do
4890 
4891  end if
4892 
4893  return

Referenced by atmos_phy_mp_sn14_terminal_velocity().

Here is the caller graph for this function:

◆ mixed_phase_collection_bin()

subroutine scale_atmos_phy_mp_sn14::mixed_phase_collection_bin ( 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(ka,i_qc:i_qg), 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), intent(in)  rho,
real(rp), dimension(ka,pq_max), intent(inout)  PQ,
real(rp), dimension(ka,pq_max), intent(inout)  Pcrg1,
real(rp), dimension(ka,pcrg_max), intent(inout)  Pcrg2,
real(rp), dimension(ka,pac_max), intent(out)  Pac 
)

0degC

Definition at line 4909 of file scale_atmos_phy_mp_sn14.F90.

4909  use scale_atmos_saturation, only: &
4910  moist_psat_ice => atmos_saturation_psat_ice
4911  use scale_prc, only: &
4912  prc_abort
4913  implicit none
4914 
4915  integer, intent(in) :: KA, KS, KE
4916 
4917  !--- mixed-phase collection process
4918  ! And all we set all production term as a negative sign to avoid confusion.
4919  !
4920  real(RP), intent(in) :: wtem(KA)
4921  !--- mass/number concentration[kg/m3]
4922  real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
4923  real(RP), intent(in) :: rho(KA) ! air density
4924  ! necessary ?
4925  real(RP), intent(in) :: xq(KA,HYDRO_MAX)
4926  !--- diameter of averaged mass( D(ave x) )
4927  real(RP), intent(in) :: dq_xave(KA,HYDRO_MAX)
4928  !--- terminal velocity of averaged mass( vt(ave x) )
4929  real(RP), intent(in) :: vt_xave(KA,HYDRO_MAX,2)
4930  ! [Add] 11/08/30 T.Mitsui, for autoconversion of ice
4931  ! real(RP), intent(in) :: rho(KA)
4932  !--- partial conversion
4933  real(RP), intent(inout):: PQ(KA,PQ_MAX)
4934  !
4935  real(RP), intent(out):: Pac(KA,Pac_MAX)
4936  !--- for lightning component
4937  logical, intent(in) :: flg_lt
4938  real(RP), intent(in) :: beta_crg(KA)
4939  real(RP), intent(in) :: dqcrg(KA)
4940  real(RP), intent(in) :: d0_crg, v0_crg
4941  real(RP), intent(in) :: rhoq_crg(KA,I_QC:I_QG)
4942  real(RP), intent(inout):: Pcrg1(KA,PQ_MAX)
4943  real(RP), intent(inout):: Pcrg2(KA,Pcrg_MAX)
4944 
4945  real(RP), parameter :: a_dec = 0.883_rp
4946  real(RP), parameter :: b_dec = 0.093_rp
4947  real(RP), parameter :: c_dec = 0.00348_rp
4948  real(RP), parameter :: d_dec = 4.5185e-5_rp
4949  !
4950  !!!!from default SN14
4951 
4952  real(RP), parameter :: E_C12(6)=(/&
4953  0.010_rp, 0.080_rp, 0.10_rp, 0.60_rp, 0.20_rp, 0.10_rp/)
4954 
4955  real(RP) :: tem(KA)
4956  !
4957  !--- collection efficency of each specie
4958  real(RP) :: E_c, E_r, E_i, E_s, E_g
4959 ! real(RP) :: E_ic, E_sc, E_gc
4960  !--- sticking efficiency
4961  real(RP):: E_stick(KA)
4962  ! [Add] 10/08/03 T.Mitsui
4963  real(RP) :: temc, temc2, temc3
4964  real(RP) :: E_dec
4965  real(RP) :: esi_rat(KA)
4966  real(RP) :: esi(KA)
4967  !
4968  real(RP) :: temc_p, temc_m ! celcius tem.
4969 ! real(RP) :: ci_aut(KA)
4970 ! real(RP) :: taui_aut(KA)
4971 ! real(RP) :: tau_sce(KA)
4972  !--- DSD averaged diameter for each species
4973  real(RP) :: ave_dc ! cloud
4974 ! real(RP) :: ave_dr ! rain
4975  real(RP) :: ave_di(KA) ! ice
4976  real(RP) :: ave_ds(KA) ! snow
4977  real(RP) :: ave_dg ! graupel
4978  !--- coefficient of collection equations(L:mass, N:number)
4979  real(RP) :: coef_acc_LCI, coef_acc_NCI ! cloud - cloud ice
4980  real(RP) :: coef_acc_LCS, coef_acc_NCS ! cloud - snow
4981  !
4982  real(RP) :: coef_acc_LCG, coef_acc_NCG ! cloud - graupel
4983  real(RP) :: coef_acc_LRI_I, coef_acc_NRI_I ! rain - cloud ice
4984  real(RP) :: coef_acc_LRI_R, coef_acc_NRI_R ! rain - cloud ice
4985  real(RP) :: coef_acc_LRS_S, coef_acc_NRS_S ! rain - snow
4986  real(RP) :: coef_acc_LRS_R, coef_acc_NRS_R ! rain - snow
4987  real(RP) :: coef_acc_LRG, coef_acc_NRG ! rain - graupel
4988  real(RP) :: coef_acc_LII, coef_acc_NII ! cloud ice - cloud ice
4989  real(RP) :: coef_acc_LIS, coef_acc_NIS ! cloud ice - snow
4990  real(RP) :: coef_acc_NSS ! snow - snow
4991  real(RP) :: coef_acc_NGG ! grauepl - graupel
4992  real(RP) :: coef_acc_LSG, coef_acc_NSG(KA) ! snow - graupel
4993  !--- (diameter) x (diameter)
4994  real(RP) :: dcdc(KA), dcdi, dcds, dcdg
4995  real(RP) :: drdr(KA), drdi(KA), drds(KA), drdg
4996  real(RP) :: didi(KA), dids, didg
4997  real(RP) :: dsds(KA), dsdg
4998  real(RP) :: dgdg(KA)
4999 !# !--- (terminal velocity) x (terminal velocity)
5000 !# real(RP) :: vcvc(KA), vcvi, vcvs, vcvg
5001 !# real(RP) :: vrvr(KA), vrvi(KA), vrvs(KA), vrvg
5002 !# real(RP) :: vivi(KA), vivs, vivg
5003 !# real(RP) :: vsvs(KA), vsvg
5004 !# real(RP) :: vgvg(KA)
5005  !
5006 ! real(RP) :: wx_cri, wx_crs
5007 ! real(RP) :: coef_emelt
5008 ! real(RP) :: w1
5009 
5010  real(RP) :: sw, sw1, sw2, tmp
5011  real(RP) :: alpha_lt
5012  !
5013  integer :: k, iqw
5014  !
5015 
5016  real(RP) :: tem_e(KA) ! [Add] 15/05/19 T.Seiki
5017 
5018 
5019  !
5020  ! work for binary collision
5021  !
5022  ! work for Gauss-Legendre quadrature
5023  integer, parameter :: ngmax=4
5024 
5025  real(RP) :: lambdac(KA), lambdar(KA), lambdai(KA), lambdas(KA), lambdag(KA)
5026  real(RP) :: A_dsdc(KA), A_dsdr(KA), A_dsdi(KA), A_dsds(KA), A_dsdg(KA)
5027  real(RP) :: dNdx
5028  real(RP) :: dxdd
5029  real(RP) :: dNdD
5030  real(RP) :: dNc_glx(KA,ngmax), dNr_glx(KA,ngmax), dNi_glx(KA,ngmax), dNs_glx(KA,ngmax), dNg_glx(KA,ngmax)
5031  real(RP) :: dNc_gly, dNr_gly(KA,ngmax), dNi_gly(KA,ngmax), dNs_gly(KA,ngmax), dNg_gly(KA,ngmax)
5032  !
5033  real(RP) :: dc_glx, dr_glx, di_glx, ds_glx, dg_glx
5034  real(RP) :: dc_gly, dr_gly, di_gly, ds_gly, dg_gly
5035  real(RP) :: xc_glx(KA,ngmax), xr_glx(KA,ngmax), xi_glx(KA,ngmax), xs_glx(KA,ngmax), xg_glx(KA,ngmax)
5036  real(RP) :: xc_gly, xr_gly(KA,ngmax), xi_gly, xs_gly, xg_gly
5037 
5038  real(RP) :: vtc_glx(KA,ngmax), vtr_glx(KA,ngmax), vti_glx(KA,ngmax), vts_glx(KA,ngmax), vtg_glx(KA,ngmax)
5039  real(RP) :: vtc_gly, vtr_gly(KA,ngmax), vti_gly(KA,ngmax), vts_gly(KA,ngmax), vtg_gly(KA,ngmax)
5040  real(RP) :: dac_glx(KA,ngmax), dar_glx(KA,ngmax), dai_glx(KA,ngmax), das_glx(KA,ngmax), dag_glx(KA,ngmax)
5041  real(RP) :: dac_gly, dar_gly(KA,ngmax), dai_gly(KA,ngmax), das_gly(KA,ngmax), dag_gly(KA,ngmax)
5042  !
5043  integer :: ngx, ngy
5044  !
5045  real(RP) :: E_ic(KA), E_sc(KA), E_gc(KA)
5046  !
5047  ! Parameters to calculate terminal velocity formulated by Mitchell (1996)
5048  !
5049  real(RP) :: acx, bcx, gcx, scx ! geometric parameters of hexagonal column ice defined by Mitchell (1996)
5050  real(RP) :: acy, bcy, gcy, scy ! geometric parameters of hexagonal column ice defined by Mitchell (1996)
5051  real(RP), parameter :: as = 0.59452551_rp
5052  real(RP), parameter :: bs = 2.4490_rp
5053  real(RP), parameter :: gs = 0.131488_rp
5054  real(RP), parameter :: ss = 1.880000_rp
5055  real(RP), parameter :: ag = 19.5072514_rp !0.049d0*1.d-3*(100.0d0**2.8d0)
5056  real(RP), parameter :: bg = 2.8_rp
5057  real(RP), parameter :: gg = 0.5_rp
5058  real(RP), parameter :: sg = 2.0_rp
5059  real(RP) :: num_Besti_glx, num_Bests_glx, num_Bestg_glx
5060  real(RP) :: num_Besti_gly, num_Bests_gly, num_Bestg_gly
5061  real(RP) :: num_Rei_glx, num_Res_glx, num_Reg_glx
5062  real(RP) :: num_Rei_gly, num_Res_gly, num_Reg_gly
5063  real(RP), parameter :: c0=0.6_rp ! Bohm (1989)
5064  real(RP), parameter :: d0=5.83_rp ! Bohm (1989)
5065  !
5066  real(RP) :: mua(KA), nua(KA)
5067  !--- Dynamic viscosity
5068  real(RP), parameter :: mua0 = 1.718e-5_rp
5069 
5070  real(RP), parameter :: dmua_dT = 5.28e-8_rp
5071 
5072  !====== mua = mua0 + temc*dmua_dT
5073  !
5074  ! collection Kernel
5075  !
5076  real(RP) :: kernel_cg, kernel_cs, kernel_ci
5077  real(RP) :: kernel_rg, kernel_rs, kernel_ri
5078  real(RP) :: kernel_ig, kernel_is, kernel_ii
5079  real(RP) :: kernel_sg, kernel_ss
5080  real(RP) :: kernel_gg
5081  real(RP) :: kernel_sg_reb
5082  !
5083  !--------------------------
5084  !
5085  ! accurate for PSDs with optimization
5086  !
5087  !-------------------------------
5088  ! cloud
5089  ! gauss_range = 2.0_RP
5090  ! ngmax = 4
5091  ! rain,
5092  ! gauss_range = xxx
5093  ! ngmax = 4
5094  ! ice, snow, graupel
5095  ! gauss_range = 5.0_RP
5096  ! ngmax = 4
5097  !
5098  real(RP), parameter :: gauss_rangec=2.0_rp
5099  real(RP), parameter :: wc_gl(ngmax)=(/&
5100  0.2411146051511425e+00_rp, &
5101  0.4520325754088027e+00_rp, &
5102  0.4520325754088027e+00_rp, &
5103  0.2411146051511425e+00_rp &
5104  /)
5105  real(RP), parameter :: coefc_d_gl(ngmax)=(/&
5106  0.5505187813766612e+00_rp, &
5107  0.7900516927471093e+00_rp, &
5108  0.1265739962562290e+01_rp, &
5109  0.1816468454535444e+01_rp &
5110  /)
5111  real(RP), parameter :: gauss_ranger=8.0_rp
5112  real(RP), parameter :: wr_gl(ngmax)=(/&
5113  0.723343815453428e+00_rp,&
5114  1.35609772622641e+00_rp, &
5115  1.35609772622641e+00_rp, &
5116  0.723343815453428e+00_rp &
5117  /)
5118  real(RP), parameter :: coefr_d_gl(ngmax)=(/&
5119  0.166846238310235e+00_rp, &
5120  0.493135790663523e+00_rp, &
5121  2.02783902311062e+00_rp, &
5122  5.99354237846583e+00_rp &
5123  /)
5124  real(RP), parameter :: gauss_range=5.0_rp
5125  real(RP), parameter :: w_gl(ngmax)=(/&
5126  0.559850775788111e+00_rp, &
5127  1.04958713664599e+00_rp, &
5128  1.04958713664599e+00_rp, &
5129  0.559850775788111e+00_rp &
5130  /)
5131  real(RP), parameter :: coef_d_gl(ngmax)=(/&
5132  0.2500872485877803e+00_rp, &
5133  0.5785800417604080e+00_rp, &
5134  0.1728369331505741e+01_rp, &
5135  0.3998604509613778e+01_rp &
5136  /)
5137  !
5138  real(RP) :: wx_cri, wx_crs
5139  real(RP) :: coef_emelt
5140  real(RP) :: w1
5141 ! integer :: ij, k
5142  integer :: ierr
5143 
5144  !---------------------------------------------------------------------------
5145  !
5146  !
5147 
5148  do k = ks, ke
5149  tem(k) = max( wtem(k), tem_min ) ! 11/08/30 T.Mitsui
5150  end do
5151 
5152  call moist_psat_ice( ka, ks, ke, &
5153  tem(:), esi(:) ) ! [IN], [OUT]
5154 
5155  if( opt_stick_ks96 )then
5156  do k = ks, ke
5157  ! Khain and Sednev (1996), eq.(3.15)
5158  temc = tem(k) - t00
5159  temc2 = temc*temc
5160  temc3 = temc2*temc
5161  e_dec = max(0.0_rp, a_dec + b_dec*temc + c_dec*temc2 + d_dec*temc3 )
5162  esi_rat(k) = rhoq(k,i_qv)*rvap*tem(k)/esi(k)
5163  e_stick(k) = min(1.0_rp, e_dec*esi_rat(k))
5164  end do
5165  else if( opt_stick_co86 )then
5166  do k = ks, ke
5167  ! [Add] 11/08/30 T.Mitsui, Cotton et al. (1986)
5168  temc = min(tem(k) - t00,0.0_rp)
5169  w1 = 0.035_rp*temc-0.7_rp
5170  e_stick(k) = 10._rp**w1
5171  end do
5172  else
5173  do k = ks, ke
5174  ! Lin et al. (1983)
5175  temc_m = min(tem(k) - t00,0.0_rp) ! T < 273.15
5176  e_stick(k) = exp(0.09_rp*temc_m)
5177  end do
5178  end if
5179 
5180  !------------------------------------------------------------------------
5181 
5182 
5183  pac(:,:) = 0.0_rp
5184 
5185  tem(:) = max(wtem(:), tem_min )
5186  tem_e(:) = max(wtem(:), tem_min_estick )
5187 
5188  call moist_psat_ice( ka, ks, ke, tem(:), esi(:) )
5189 
5190 
5191  if ( opt_stick_ks96 ) then
5192  do k=ks, ke
5193  ! Khain and Sednev (1996), eq.(3.15)
5194  temc = tem_e(k) - t00 ![Mod] 15/05/19 T.Seiki
5195  temc2 = temc*temc
5196  temc3 = temc2*temc
5197  e_dec = max(0.0_rp, a_dec + b_dec*temc + c_dec*temc2 + d_dec*temc3 )
5198  esi_rat(k) = rhoq(k,i_qv)*rvap*tem(k)/esi(k)
5199  e_stick(k) = min(1.0_rp, e_dec*esi_rat(k))
5200  enddo
5201  elseif( opt_stick_co86 ) then
5202  do k=ks, ke
5203  temc = min(tem_e(k) - t00,0.0_rp)
5204  w1 = 0.035_rp*temc-0.7_rp
5205  e_stick(k) = 10.0_rp**w1
5206  enddo
5207  elseif( opt_stick_c12 ) then
5208  do k=ks, ke
5209  if (tem_e(k)>273.15_rp) then
5210  e_stick(k)=1.0_rp
5211  elseif(tem_e(k)<243.15_rp) then !-30degC
5212  e_stick(k)=e_c12(1)
5213  elseif(tem_e(k)<248.15_rp) then !-25degC
5214  e_stick(k)=e_c12(1)+0.2_rp*(e_c12(2)-e_c12(1))*(tem_e(k)-243.15_rp)
5215  elseif(tem_e(k)<253.15_rp) then !-20degC
5216  e_stick(k)=e_c12(2)+0.2_rp*(e_c12(3)-e_c12(2))*(tem_e(k)-248.15_rp)
5217  elseif(tem_e(k)<258.15_rp) then !-15degC
5218  e_stick(k)=e_c12(3)+0.2_rp*(e_c12(4)-e_c12(3))*(tem_e(k)-253.15_rp)
5219  elseif(tem_e(k)<263.15_rp) then !-10degC
5220  e_stick(k)=e_c12(4)+0.2_rp*(e_c12(5)-e_c12(4))*(tem_e(k)-258.15_rp)
5221  elseif(tem_e(k)<268.15_rp) then !- 5degC
5222  e_stick(k)=e_c12(5)+0.2_rp*(e_c12(6)-e_c12(5))*(tem_e(k)-263.15_rp)
5223  else ! 0degC
5224  e_stick(k)=e_c12(6)
5225  endif
5226  enddo
5227  else
5228  if ( opt_stick_rhh57 ) then
5229  do k=ks, ke
5230  if ( tem_e(k) < 270.0_rp .AND. rhoq(k,i_qv)*rvap*tem(k) < esi(k) ) then
5231  esi_rat(k) = 0.0_rp
5232  else
5233  esi_rat(k) = 1.0_rp
5234  endif
5235  enddo
5236  elseif( opt_stick_rhks96 ) then
5237  do k=ks, ke
5238  esi_rat(k) = min( rhoq(k,i_qv)*rvap*tem(k)/esi(k),1.0_rp )
5239  enddo
5240  else
5241  esi_rat(:) = 1.0_rp
5242  endif
5243  !
5244  do k=ks, ke
5245  ! Lin et al. (1983)
5246  temc_m = min(tem_e(k) - t00,0.0_rp) ! T < 273.15
5247  e_stick(k) = exp(0.09_rp*temc_m)*esi_rat(k)
5248  enddo
5249  endif
5250 
5251  !
5252  ! Integration Start
5253  !
5254 
5255  do k = ks, ke
5256  mua(k) = mua0 + dmua_dt*(tem(k)-273.15_rp)
5257  nua(k) = mua(k)/rho(k) ! [m2/s]
5258  end do
5259  do k = ks, ke
5260  lambdac(k) = xq(k,i_mp_qc)**(-mu(i_mp_qc))*coef_lambda(i_mp_qc)
5261  a_dsdc(k) = rhoq(k,i_nc)*coef_a(i_mp_qc)*lambdac(k)**((nu(i_mp_qc)+1.0_rp)/mu(i_mp_qc))
5262  end do
5263  do k = ks, ke
5264  lambdar(k) = xq(k,i_mp_qr)**(-mu(i_mp_qr))*coef_lambda(i_mp_qr)
5265  a_dsdr(k) = rhoq(k,i_nr)*coef_a(i_mp_qr)*lambdar(k)**((nu(i_mp_qr)+1.0_rp)/mu(i_mp_qr))
5266  end do
5267  do k = ks, ke
5268  lambdai(k) = xq(k,i_mp_qi)**(-mu(i_mp_qi))*coef_lambda(i_mp_qi)
5269  a_dsdi(k) = rhoq(k,i_ni)*coef_a(i_mp_qi)*lambdai(k)**((nu(i_mp_qi)+1.0_rp)/mu(i_mp_qi))
5270  end do
5271  do k = ks, ke
5272  lambdas(k) = xq(k,i_mp_qs)**(-mu(i_mp_qs))*coef_lambda(i_mp_qs)
5273  a_dsds(k) = rhoq(k,i_ns)*coef_a(i_mp_qs)*lambdas(k)**((nu(i_mp_qs)+1.0_rp)/mu(i_mp_qs))
5274  end do
5275  do k = ks, ke
5276  lambdag(k) = xq(k,i_mp_qg)**(-mu(i_mp_qg))*coef_lambda(i_mp_qg)
5277  a_dsdg(k) = rhoq(k,i_ng)*coef_a(i_mp_qg)*lambdag(k)**((nu(i_mp_qg)+1.0_rp)/mu(i_mp_qg))
5278  end do
5279 
5280  !
5281  ! Particle X
5282  !
5283  do ngx=1, ngmax ! X < Y
5284 !OCL LOOP_FISSION_TARGET(LS)
5285  do k = ks, ke
5286  dc_glx = dq_xave(k,i_mp_qc)*coefc_d_gl(ngx)
5287  xc_glx(k,ngx) = ( (dc_glx/a_m(i_mp_qc)) )**(1.0_rp/b_m(i_mp_qc))
5288  dnc_glx(k,ngx) = a_dsdc(k)*(xc_glx(k,ngx)**nu(i_mp_qc)) * exp(-lambdac(k)*xc_glx(k,ngx)**mu(i_mp_qc))&
5289  *(xc_glx(k,ngx)/(b_m(i_mp_qc)*dc_glx))*dc_glx*wc_gl(ngx) ! dNdlogD*weight
5290  !
5291  dr_glx = dq_xave(k,i_mp_qr)*coefr_d_gl(ngx)
5292  xr_glx(k,ngx) = ( (dr_glx/a_m(i_mp_qr)) )**(1.0_rp/b_m(i_mp_qr))
5293  dnr_glx(k,ngx) = a_dsdr(k)*(xr_glx(k,ngx)**nu(i_mp_qr)) * exp(-lambdar(k)*xr_glx(k,ngx)**mu(i_mp_qr))&
5294  *(xr_glx(k,ngx)/(b_m(i_mp_qr)*dr_glx))*dr_glx*wr_gl(ngx) ! dNdlogD*weight
5295  !
5296  di_glx = dq_xave(k,i_mp_qi)*coef_d_gl(ngx)
5297  xi_glx(k,ngx) = ( (di_glx/a_m(i_mp_qi)) )**(1.0_rp/b_m(i_mp_qi))
5298  dni_glx(k,ngx) = a_dsdi(k)*(xi_glx(k,ngx)**nu(i_mp_qi)) * exp(-lambdai(k)*xi_glx(k,ngx)**mu(i_mp_qi))&
5299  *(xi_glx(k,ngx)/(b_m(i_mp_qi)*di_glx))*di_glx*w_gl(ngx) ! dNdlogD*weight
5300  !
5301  ds_glx = dq_xave(k,i_mp_qs)*coef_d_gl(ngx)
5302  xs_glx(k,ngx) = ( (ds_glx/a_m(i_mp_qs)) )**(1.0_rp/b_m(i_mp_qs))
5303  dns_glx(k,ngx) = a_dsds(k)*(xs_glx(k,ngx)**nu(i_mp_qs)) * exp(-lambdas(k)*xs_glx(k,ngx)**mu(i_mp_qs))&
5304  *(xs_glx(k,ngx)/(b_m(i_mp_qs)*ds_glx))*ds_glx*w_gl(ngx)! dNdlogD*weight
5305  !
5306  dg_glx = dq_xave(k,i_mp_qg)*coef_d_gl(ngx)
5307  xg_glx(k,ngx) = ( (dg_glx/a_m(i_mp_qg)) )**(1.0_rp/b_m(i_mp_qg))
5308  dng_glx(k,ngx) = a_dsdg(k)*(xg_glx(k,ngx)**nu(i_mp_qg)) * exp(-lambdag(k)*xg_glx(k,ngx)**mu(i_mp_qg))&
5309  *(xg_glx(k,ngx)/(b_m(i_mp_qg)*dg_glx))*dg_glx*w_gl(ngx)! dNdlogD*weight
5310  !
5311  ! Hexagonal Columns
5312 !!$ if( di_glx <= 100.e-6_RP )then
5313 !!$ acx = 0.1677_RP*1.e-3_RP*(100.0_RP**2.91_RP)
5314 !!$ bcx = 2.91_RP
5315 !!$ gcx = (0.684_RP*1.e-4_RP)*10.0_RP**(2.0_RP*2.0_RP)
5316 !!$ scx = 2.0_RP
5317 !!$ else
5318 !!$ acx = 0.00166_RP*1.e-3_RP*(100.0_RP**1.91_RP)
5319 !!$ bcx = 1.91_RP
5320 !!$ gcx = (0.0696_RP*1.e-4_RP)*10.0_RP**(2.0_RP*1.5_RP)
5321 !!$ scx = 1.5_RP
5322 !!$ end if
5323  sw = 0.5_rp + sign(0.5_rp, di_glx - 100.e-6_rp )
5324  acx = ( 0.1677_rp*(1.0_rp-sw) + 0.00166_rp*sw ) * 1.e-3_rp * 100.0_rp**( 2.91_rp*(1.0_rp-sw) + 1.91_rp*sw )
5325  bcx = 2.91_rp*(1.0_rp-sw) + 1.91_rp*sw
5326  gcx = (0.684_rp*(1.0_rp-sw) + 0.0696_rp*sw ) * 1.e-4_rp * 10.0_rp**( 4.0_rp*(1.0_rp-sw) + 3.0_rp*sw )
5327  scx = 2.0_rp*(1.0_rp-sw) + 1.5_rp*sw
5328  num_besti_glx = 2.0_rp*acx*grav*rho(k)*di_glx**(bcx+2.0_rp-scx)/(gcx*mua(k)**2)
5329  num_bests_glx = 2.0_rp*as *grav*rho(k)*ds_glx**(bs +2.0_rp-ss )/(gs *mua(k)**2)
5330  num_bestg_glx = 2.0_rp*ag *grav*rho(k)*dg_glx**(bg +2.0_rp-sg )/(gg *mua(k)**2)
5331  num_rei_glx = 0.25_rp*d0*d0*( sqrt(1.0_rp+4.0_rp*sqrt(num_besti_glx)/(d0*d0*sqrt(c0)))-1.0_rp )**2
5332  num_res_glx = 0.25_rp*d0*d0*( sqrt(1.0_rp+4.0_rp*sqrt(num_bests_glx)/(d0*d0*sqrt(c0)))-1.0_rp )**2
5333  num_reg_glx = 0.25_rp*d0*d0*( sqrt(1.0_rp+4.0_rp*sqrt(num_bestg_glx)/(d0*d0*sqrt(c0)))-1.0_rp )**2
5334  !
5335  vtc_glx(k,ngx) = coef_vtr_ar2*dc_glx*(1.0_rp-exp(-coef_vtr_br2*dc_glx))
5336  !
5337 !!$ if( dr_glx < d_vtr_branch )then
5338 !!$ vtr_glx = coef_vtr_ar2*dr_glx*(1.0_RP-exp(-coef_vtr_br2*dr_glx))
5339 !!$ else
5340 !!$ vtr_glx = coef_vtr_ar1-coef_vtr_br1*exp(-coef_vtr_cr1*dr_glx)
5341 !!$ end if
5342  sw = 0.5_rp + sign( 0.5_rp, dr_glx - d_vtr_branch )
5343  tmp = exp( - ( coef_vtr_br2*(1.0_rp-sw) + coef_vtr_cr1*sw ) * dr_glx )
5344  vtr_glx(k,ngx) = coef_vtr_ar2 * dr_glx * ( 1.0_rp - tmp ) * (1.0_rp-sw) &
5345  + ( coef_vtr_ar1 - coef_vtr_br1 * tmp ) * sw
5346  !
5347  vti_glx(k,ngx) = num_rei_glx*nua(k)/di_glx
5348  vts_glx(k,ngx) = num_res_glx*nua(k)/ds_glx
5349  vtg_glx(k,ngx) = num_reg_glx*nua(k)/dg_glx
5350  ! equivalent area diameter
5351  dac_glx(k,ngx) = dc_glx
5352  dar_glx(k,ngx) = dr_glx
5353  dai_glx(k,ngx) = 2.0_rp*sqrt( (gcx*di_glx**scx)/pi )
5354  das_glx(k,ngx) = 2.0_rp*sqrt( (gs *ds_glx**ss )/pi )
5355  dag_glx(k,ngx) = 2.0_rp*sqrt( (gg *dg_glx**sg )/pi )
5356  end do
5357  end do
5358 
5359  !
5360  ! Particle Y
5361  !
5362  do ngy=1, ngmax
5363 !OCL LOOP_FISSION_TARGET(LS)
5364  do k = ks, ke
5365  dr_gly = dq_xave(k,i_mp_qr)*coefr_d_gl(ngy)
5366  xr_gly(k,ngy) = ( (dr_gly/a_m(i_mp_qr)) )**(1.0_rp/b_m(i_mp_qr))
5367  dnr_gly(k,ngy) = a_dsdr(k)*(xr_gly(k,ngy)**nu(i_mp_qr)) * exp(-lambdar(k)*xr_gly(k,ngy)**mu(i_mp_qr))&
5368  *(xr_gly(k,ngy)/(b_m(i_mp_qr)*dr_gly))*dr_gly*wr_gl(ngy) ! dNdlogD*weight
5369  !
5370  di_gly = dq_xave(k,i_mp_qi)*coef_d_gl(ngy)
5371  xi_gly = ( (di_gly/a_m(i_mp_qi)) )**(1.0_rp/b_m(i_mp_qi))
5372  dni_gly(k,ngy) = a_dsdi(k)*(xi_gly**nu(i_mp_qi)) * exp(-lambdai(k)*xi_gly**mu(i_mp_qi))&
5373  *(xi_gly/(b_m(i_mp_qi)*di_gly))*di_gly*w_gl(ngy) ! dNdlogD*weight
5374  !
5375  ds_gly = dq_xave(k,i_mp_qs)*coef_d_gl(ngy)
5376  xs_gly = ( (ds_gly/a_m(i_mp_qs)) )**(1.0_rp/b_m(i_mp_qs))
5377  dns_gly(k,ngy) = a_dsds(k)*(xs_gly**nu(i_mp_qs)) * exp(-lambdas(k)*xs_gly**mu(i_mp_qs))&
5378  *(xs_gly/(b_m(i_mp_qs)*ds_gly))*ds_gly*w_gl(ngy)! dNdlogD*weight
5379  !
5380  dg_gly = dq_xave(k,i_mp_qg)*coef_d_gl(ngy)
5381  xg_gly = ( (dg_gly/a_m(i_mp_qg)) )**(1.0_rp/b_m(i_mp_qg))
5382  dng_gly(k,ngy) = a_dsdg(k)*(xg_gly**nu(i_mp_qg)) * exp(-lambdag(k)*xg_gly**mu(i_mp_qg))&
5383  *(xg_gly/(b_m(i_mp_qg)*dg_gly))*dg_gly*w_gl(ngy)! dNdlogD*weight
5384  !
5385  ! Hexagonal Columns
5386 !!$ if( di_gly <= 100.e-6_RP )then
5387 !!$ acy = 0.1677_RP*1.d-3*(100.0_RP**2.91_RP)
5388 !!$ bcy = 2.91_RP
5389 !!$ gcy = (0.684_RP*1.e-4_RP)*10.0_RP**(2.0_RP*2.0_RP)
5390 !!$ scy = 2.0_RP
5391 !!$ else
5392 !!$ acy = 0.00166_RP*1.e-3_RP*(100.0_RP**1.91_RP)
5393 !!$ bcy = 1.91_RP
5394 !!$ gcy = (0.0696_RP*1.e-4_RP)*10.0_RP**(2.0_RP*1.5_RP)
5395 !!$ scy = 1.5_RP
5396 !!$ end if
5397  sw = 0.5_rp + sign( 0.5_rp, di_gly - 100.e-6_rp )
5398  acy = ( 0.1677_rp*(1.0_rp-sw) + 0.00166_rp*sw ) * 1.e-3_rp * 100.0_rp**( 2.91_rp*(1.0_rp-sw) + 1.91_rp*sw )
5399  bcy = 2.91_rp*(1.0_rp-sw) + 1.91_rp*sw
5400  gcy = ( 0.684_rp*(1.0_rp-sw) + 0.0696_rp*sw ) * 1.e-4_rp * 10.0_rp**( 4.0_rp*(1.0_rp-sw) + 3.0_rp*sw )
5401  scy = 2.0_rp*(1.0_rp-sw) + 1.5_rp*sw
5402  num_besti_gly = 2.0_rp*acy*grav*rho(k)*di_gly**(bcy+2.0_rp-scy)/(gcy*mua(k)**2)
5403  num_bests_gly = 2.0_rp*as *grav*rho(k)*ds_gly**(bs +2.0_rp-ss )/(gs *mua(k)**2)
5404  num_bestg_gly = 2.0_rp*ag *grav*rho(k)*dg_gly**(bg +2.0_rp-sg )/(gg *mua(k)**2)
5405  num_rei_gly = 0.25_rp*d0*d0*( sqrt(1.0_rp+4.0_rp*sqrt(num_besti_gly)/(d0*d0*sqrt(c0)))-1.0_rp )**2
5406  num_res_gly = 0.25_rp*d0*d0*( sqrt(1.0_rp+4.0_rp*sqrt(num_bests_gly)/(d0*d0*sqrt(c0)))-1.0_rp )**2
5407  num_reg_gly = 0.25_rp*d0*d0*( sqrt(1.0_rp+4.0_rp*sqrt(num_bestg_gly)/(d0*d0*sqrt(c0)))-1.0_rp )**2
5408  !
5409 !!$ if( dr_gly < d_vtr_branch )then
5410 !!$ vtr_gly(k,ngy) = coef_vtr_ar2*dr_gly*(1.0_RP-exp(-coef_vtr_br2*dr_gly))
5411 !!$ else
5412 !!$ vtr_gly(k,ngy) = coef_vtr_ar1-coef_vtr_br1*exp(-coef_vtr_cr1*dr_gly)
5413 !!$ end if
5414  sw = 0.5_rp + sign( 0.5_rp, dr_gly - d_vtr_branch )
5415  tmp = exp( - ( coef_vtr_br2*(1.0_rp-sw) + coef_vtr_cr1*sw ) * dr_gly )
5416  vtr_gly(k,ngy) = coef_vtr_ar2 * dr_gly * ( 1.0_rp - tmp ) * (1.0_rp-sw) &
5417  + ( coef_vtr_ar1 - coef_vtr_br1 * tmp ) * sw
5418  !
5419  vti_gly(k,ngy) = num_rei_gly*nua(k)/di_gly
5420  vts_gly(k,ngy) = num_res_gly*nua(k)/ds_gly
5421  vtg_gly(k,ngy) = num_reg_gly*nua(k)/dg_gly
5422  ! equivalent area diameter
5423  dar_gly(k,ngy) = dr_gly
5424  dai_gly(k,ngy) = 2.0_rp*sqrt( (gcy*di_gly**scy)/pi )
5425  das_gly(k,ngy) = 2.0_rp*sqrt( (gs *ds_gly**ss )/pi )
5426  dag_gly(k,ngy) = 2.0_rp*sqrt( (gg *dg_gly**sg )/pi )
5427  end do
5428  end do
5429 
5430  !
5431  ! BULK collection efficiency are given as follows
5432  !
5433  do k = ks, ke
5434  e_c = max(0.0_rp, min(1.0_rp, (dq_xave(k,i_mp_qc)-dc0)/(dc1-dc0) ))
5435  !
5436 !!$ if (dq_xave(k,I_mp_QI)>di0) then
5437 !!$ E_i = E_im
5438 !!$ else
5439 !!$ E_i = 0.0_RP
5440 !!$ endif
5441  sw = 0.5_rp + sign( 0.5_rp, dq_xave(k,i_mp_qi)-di0 )
5442  e_i = e_im * sw
5443 !!$ if (dq_xave(k,I_mp_QS)>ds0) then
5444 !!$ E_s = E_sm
5445 !!$ else
5446 !!$ E_s = 0.0_RP
5447 !!$ endif
5448  sw = 0.5_rp + sign( 0.5_rp, dq_xave(k,i_mp_qs)-ds0 )
5449  e_s = e_sm * sw
5450 !!$ if (dq_xave(k,I_mp_QG)>dg0) then
5451 !!$ E_g = E_gm
5452 !!$ else
5453 !!$ E_g = 0.0_RP
5454 !!$ endif
5455  sw = 0.5_rp + sign( 0.5_rp, dq_xave(k,i_mp_qg)-dg0 )
5456  e_g = e_gm * sw
5457 
5458  e_ic(k) = e_i*e_c
5459  e_sc(k) = e_s*e_c
5460  e_gc(k) = e_g*e_c
5461  end do
5462 
5463 
5464  !=========================================================================================
5465  ! collection equation
5466  !=========================================================================================
5467  do ngx=1, ngmax ! X < Y
5468  do ngy=1, ngmax
5469 
5470  do k = ks, ke
5471  !
5472  ! 1.c-g (X=Cloud, Y=Graupel)
5473  !
5474 ! kernel_cg = 0.25_RP*pi*(dag_gly(k,ngy)+dac_glx(k,ngx))*(dag_gly(k,ngy)+dac_glx(k,ngx))*sqrt((vtg_gly(k,ngy)-vtc_glx(k,ngx))*(vtg_gly(k,ngy)-vtc_glx(k,ngx))) * E_gc(k)
5475  kernel_cg = 0.25_rp * pi * (dag_gly(k,ngy)+dac_glx(k,ngx))**2 * abs(vtg_gly(k,ngy)-vtc_glx(k,ngx)) * e_gc(k)
5476  pac(k,i_ngacnc2ng) = pac(k,i_ngacnc2ng) - kernel_cg *dnc_glx(k,ngx)*dng_gly(k,ngy)
5477  pac(k,i_lgaclc2lg) = pac(k,i_lgaclc2lg) - kernel_cg*xc_glx(k,ngx)*dnc_glx(k,ngx)*dng_gly(k,ngy)
5478  end do
5479  do k = ks, ke
5480  !
5481  ! 2.c-s (X=Cloud, Y=Snow)
5482  !
5483 ! kernel_cs = 0.25_RP*pi*(das_gly(k,ngy)+dac_glx(k,ngx))*(das_gly(k,ngy)+dac_glx(k,ngx))*sqrt((vts_gly(k,ngy)-vtc_glx(k,ngx))*(vts_gly(k,ngy)-vtc_glx(k,ngx))) * E_sc(k)
5484  kernel_cs = 0.25_rp * pi * (das_gly(k,ngy)+dac_glx(k,ngx))**2 * abs(vts_gly(k,ngy)-vtc_glx(k,ngx)) * e_sc(k)
5485  pac(k,i_nsacnc2ns) = pac(k,i_nsacnc2ns) - kernel_cs *dnc_glx(k,ngx)*dns_gly(k,ngy)
5486  pac(k,i_lsaclc2ls) = pac(k,i_lsaclc2ls) - kernel_cs*xc_glx(k,ngx)*dnc_glx(k,ngx)*dns_gly(k,ngy)
5487  end do
5488  do k = ks, ke
5489  !
5490  ! 3.c-i (X=Cloud, Y=Cloud Ice)
5491  !
5492 ! kernel_ci = 0.25_RP*pi*(dai_gly(k,ngy)+dac_glx(k,ngx))*(dai_gly(k,ngy)+dac_glx(k,ngx))*sqrt((vti_gly(k,ngy)-vtc_glx(k,ngx))*(vti_gly(k,ngy)-vtc_glx(k,ngx))) * E_ic(k)
5493  kernel_ci = 0.25_rp * pi * (dai_gly(k,ngy)+dac_glx(k,ngx))**2 * abs(vti_gly(k,ngy)-vtc_glx(k,ngx)) * e_ic(k)
5494  pac(k,i_niacnc2ni) = pac(k,i_niacnc2ni) - kernel_ci *dnc_glx(k,ngx)*dni_gly(k,ngy)
5495  pac(k,i_liaclc2li) = pac(k,i_liaclc2li) - kernel_ci*xc_glx(k,ngx)*dnc_glx(k,ngx)*dni_gly(k,ngy)
5496  end do
5497  do k = ks, ke
5498  !
5499  ! 4.r-g
5500  !
5501 ! kernel_rg = 0.25_RP*pi*(dag_gly(k,ngy)+dar_glx(k,ngx))*(dag_gly(k,ngy)+dar_glx(k,ngx))*sqrt((vtg_gly(k,ngy)-vtr_glx(k,ngx))*(vtg_gly(k,ngy)-vtr_glx(k,ngx))) * E_gr
5502  kernel_rg = 0.25_rp * pi * (dag_gly(k,ngy)+dar_glx(k,ngx))**2 * abs(vtg_gly(k,ngy)-vtr_glx(k,ngx)) * e_gr
5503  ! T < 273K (X=Rain , Y=Graupel)
5504  pac(k,i_nracng2ng) = pac(k,i_nracng2ng) - kernel_rg *dnr_glx(k,ngx)*dng_gly(k,ngy)
5505  pac(k,i_lraclg2lg) = pac(k,i_lraclg2lg) - kernel_rg*xr_glx(k,ngx)*dnr_glx(k,ngx)*dng_gly(k,ngy)
5506  ! T > 273K (X=Graupel, Y=Rain )
5507  pac(k,i_nracng2nr) = pac(k,i_nracng2nr) - kernel_rg *dng_glx(k,ngx)*dnr_gly(k,ngy)
5508  pac(k,i_lraclg2lr) = pac(k,i_lraclg2lr) - kernel_rg*xg_glx(k,ngx)*dng_glx(k,ngx)*dnr_gly(k,ngy)
5509  end do
5510  do k = ks, ke
5511  !
5512  ! 5.r-s
5513  !
5514 ! kernel_rs = 0.25_RP*pi*(das_glx(k,ngx)+dar_gly(k,ngy))*(das_glx(k,ngx)+dar_gly(k,ngy))*sqrt((vts_glx(k,ngx)-vtr_gly(k,ngy))*(vts_glx(k,ngx)-vtr_gly(k,ngy))) * E_sr
5515  kernel_rs = 0.25_rp * pi * (das_glx(k,ngx)+dar_gly(k,ngy))**2 * abs(vts_glx(k,ngx)-vtr_gly(k,ngy)) * e_sr
5516  ! (X=Snow, Y=Rain)
5517  pac(k,i_nracns2ng_r) = pac(k,i_nracns2ng_r) - kernel_rs *dns_glx(k,ngx)*dnr_gly(k,ngy)
5518  pac(k,i_lracls2lg_r) = pac(k,i_lracls2lg_r) - kernel_rs*xr_gly(k,ngy)*dns_glx(k,ngx)*dnr_gly(k,ngy)
5519  !
5520  pac(k,i_nracns2ng_s) = pac(k,i_nracns2ng_s) - kernel_rs *dns_glx(k,ngx)*dnr_gly(k,ngy)
5521  pac(k,i_lracls2lg_s) = pac(k,i_lracls2lg_s) - kernel_rs*xs_glx(k,ngx)*dns_glx(k,ngx)*dnr_gly(k,ngy)
5522  end do
5523  do k = ks, ke
5524  !
5525  ! 6.r-i
5526  !
5527 ! kernel_ri = 0.25_RP*pi*(dai_glx(k,ngx)+dar_gly(k,ngy))*(dai_glx(k,ngx)+dar_gly(k,ngy))*sqrt((vti_glx(k,ngx)-vtr_gly(k,ngy))*(vti_glx(k,ngx)-vtr_gly(k,ngy))) * E_ir
5528  kernel_ri = 0.25_rp * pi * (dai_glx(k,ngx)+dar_gly(k,ngy))**2 * abs(vti_glx(k,ngx)-vtr_gly(k,ngy)) * e_ir
5529  ! (X=Cloud Ice, Y=Rain)
5530  pac(k,i_nracni2ng_r) = pac(k,i_nracni2ng_r) - kernel_ri *dni_glx(k,ngx)*dnr_gly(k,ngy)
5531  pac(k,i_lracli2lg_r) = pac(k,i_lracli2lg_r) - kernel_ri*xr_gly(k,ngy)*dni_glx(k,ngx)*dnr_gly(k,ngy)
5532  !
5533  pac(k,i_nracni2ng_i) = pac(k,i_nracni2ng_i) - kernel_ri *dni_glx(k,ngx)*dnr_gly(k,ngy)
5534  pac(k,i_lracli2lg_i) = pac(k,i_lracli2lg_i) - kernel_ri*xi_glx(k,ngx)*dni_glx(k,ngx)*dnr_gly(k,ngy)
5535  end do
5536  do k = ks, ke
5537  !
5538  ! 7.i-g
5539  !
5540 ! kernel_ig = 0.25_RP*pi*(dai_glx(k,ngx)+dag_gly(k,ngy))*(dai_glx(k,ngx)+dag_gly(k,ngy))*sqrt((vti_glx(k,ngx)-vtg_gly(k,ngy))*(vti_glx(k,ngx)-vtg_gly(k,ngy))) * E_stick(k) * E_gi
5541  kernel_ig = 0.25_rp * pi * (dai_glx(k,ngx)+dag_gly(k,ngy))**2 * abs(vti_glx(k,ngx)-vtg_gly(k,ngy)) * e_stick(k) * e_gi
5542  pac(k,i_niacng2ng) = pac(k,i_niacng2ng) - kernel_ig *dni_glx(k,ngx)*dng_gly(k,ngy)
5543  pac(k,i_liaclg2lg) = pac(k,i_liaclg2lg) - kernel_ig*xi_glx(k,ngx)*dni_glx(k,ngx)*dng_gly(k,ngy)
5544  end do
5545  do k = ks, ke
5546  !
5547  ! 8.i-s
5548  !
5549 ! kernel_is = 0.25_RP*pi*(dai_glx(k,ngx)+das_gly(k,ngy))*(dai_glx(k,ngx)+das_gly(k,ngy))*sqrt((vti_glx(k,ngx)-vts_gly(k,ngy))*(vti_glx(k,ngx)-vts_gly(k,ngy))) * E_stick(k) * E_si
5550  kernel_is = 0.25_rp * pi * (dai_glx(k,ngx)+das_gly(k,ngy))**2 * abs(vti_glx(k,ngx)-vts_gly(k,ngy)) * e_stick(k) * e_si
5551  pac(k,i_niacns2ns) = pac(k,i_niacns2ns) - kernel_is *dni_glx(k,ngx)*dns_gly(k,ngy)
5552  pac(k,i_liacls2ls) = pac(k,i_liacls2ls) - kernel_is*xi_glx(k,ngx)*dni_glx(k,ngx)*dns_gly(k,ngy)
5553  end do
5554  do k = ks, ke
5555  !
5556  ! 9.i-i
5557  !
5558 ! kernel_ii = 0.25_RP*pi*(dai_glx(k,ngx)+dai_gly(k,ngy))*(dai_glx(k,ngx)+dai_gly(k,ngy))*sqrt((vti_glx(k,ngx)-vti_gly(k,ngy))*(vti_glx(k,ngx)-vti_gly(k,ngy))) * E_stick(k) * E_ii
5559  kernel_ii = 0.25_rp * pi * (dai_glx(k,ngx)+dai_gly(k,ngy))**2 * abs(vti_glx(k,ngx)-vti_gly(k,ngy)) * e_stick(k) * e_ii
5560  pac(k,i_niacni2ns) = pac(k,i_niacni2ns) - kernel_ii *dni_glx(k,ngx)*dni_gly(k,ngy)
5561  pac(k,i_liacli2ls) = pac(k,i_liacli2ls) - kernel_ii*xi_glx(k,ngx)*dni_glx(k,ngx)*dni_gly(k,ngy)
5562  end do
5563  do k = ks, ke
5564  !
5565  ! 10.s-g
5566  !
5567 ! kernel_sg = 0.25_RP*pi*(das_glx(k,ngx)+dag_gly(k,ngy))*(das_glx(k,ngx)+dag_gly(k,ngy))*sqrt((vts_glx(k,ngx)-vtg_gly(k,ngy))*(vts_glx(k,ngx)-vtg_gly(k,ngy))) * E_stick(k) * E_gs
5568  kernel_sg = 0.25_rp * pi * (das_glx(k,ngx)+dag_gly(k,ngy))**2 * abs(vts_glx(k,ngx)-vtg_gly(k,ngy)) * e_stick(k) * e_gs
5569  pac(k,i_ngacns2ng) = pac(k,i_ngacns2ng) - kernel_sg *dns_glx(k,ngx)*dng_gly(k,ngy)
5570  pac(k,i_lgacls2lg) = pac(k,i_lgacls2lg) - kernel_sg*xs_glx(k,ngx)*dns_glx(k,ngx)*dng_gly(k,ngy)
5571  end do
5572  do k = ks, ke
5573  !
5574  ! 11.s-s
5575  !
5576 ! kernel_ss = 0.125_RP*pi*(das_glx(k,ngx)+das_gly(k,ngy))*(das_glx(k,ngx)+das_gly(k,ngy))*sqrt((vts_glx(k,ngx)-vts_gly(k,ngy))*(vts_glx(k,ngx)-vts_gly(k,ngy))) * E_stick(k) * E_ss
5577  kernel_ss = 0.125_rp * pi * (das_glx(k,ngx)+das_gly(k,ngy))**2 * abs(vts_glx(k,ngx)-vts_gly(k,ngy)) * e_stick(k) * e_ss
5578  pac(k,i_nsacns2ns) = pac(k,i_nsacns2ns) - kernel_ss*dns_glx(k,ngx)*dns_gly(k,ngy)
5579  end do
5580  do k = ks, ke
5581  !
5582  ! 12.g-g
5583  !
5584 ! kernel_gg = 0.125_RP*pi*(dag_glx(k,ngx)+dag_gly(k,ngy))*(dag_glx(k,ngx)+dag_gly(k,ngy))*sqrt((vtg_glx(k,ngx)-vtg_gly(k,ngy))*(vtg_glx(k,ngx)-vtg_gly(k,ngy))) * E_stick(k) * E_gg
5585  kernel_gg = 0.125_rp * pi * (dag_glx(k,ngx)+dag_gly(k,ngy))**2 * abs(vtg_glx(k,ngx)-vtg_gly(k,ngy)) * e_stick(k) * e_gg
5586  pac(k,i_ngacng2ng) = pac(k,i_ngacng2ng) - kernel_gg*dng_glx(k,ngx)*dng_gly(k,ngy)
5587 
5588  end do
5589  end do
5590  end do
5591 
5592  !--- Charge separation by Snow-Graupel rebound--------------------------------
5593  if ( flg_lt ) then
5594  do k = ks, ke
5595 !OCL UNROLL('full')
5596  do ngy = 1, ngmax
5597 !OCL UNROLL('full')
5598  do ngx = 1, ngmax
5599  alpha_lt = 5.0_rp * ( das_glx(k,ngx) / d0_crg )**2*vtg_gly(k,ngy)/v0_crg
5600  alpha_lt = min( alpha_lt, 10.0_rp )
5601 ! kernel_sg_reb = 0.25_RP*pi*(das_glx(k,ngx)+dag_gly(k,ngy))*(das_glx(k,ngx)+dag_gly(k,ngy))*sqrt((vts_glx(k,ngx)-vtg_gly(k,ngy))*(vts_glx(k,ngx)-vtg_gly(k,ngy))) &
5602  kernel_sg_reb = 0.25_rp * pi * (das_glx(k,ngx)+dag_gly(k,ngy))**2 * abs(vts_glx(k,ngx)-vtg_gly(k,ngy)) &
5603  * ( 1.0_rp - e_stick(k) ) * e_gs
5604  pcrg2(k,i_cgngacns2ng) = pcrg2(k,i_cgngacns2ng) + kernel_sg_reb*dns_glx(k,ngx)*dng_gly(k,ngy)*dqcrg(k)*alpha_lt*beta_crg(k)
5605  end do
5606  end do
5607  end do
5608  end if
5609 
5610  !
5611  !
5612  do k=ks, ke
5613  temc_p = max(tem(k) - t00,0.0_rp) ! T > 273.15
5614  !------------------------------------------------------------------------
5615  !--- Partial conversion
5616  ! SB06(70),(71)
5617  ! i_iconv2g: option whether partial conversions work or not
5618  ! ice-cloud => graupel
5619  if ( dq_xave(k,i_mp_qi) > di_cri ) then
5620  wx_cri = cfill_i*rhow/rho_g*( pi/6.0_rp*rho_g*dq_xave(k,i_mp_qi)*dq_xave(k,i_mp_qi)*dq_xave(k,i_mp_qi)/xq(k,i_mp_qi) - 1.0_rp )
5621  pq(k,i_licon) = i_iconv2g* pac(k,i_liaclc2li)/max(1.0_rp, wx_cri)
5622  pq(k,i_nicon) = i_iconv2g* pq(k,i_licon)/xq(k,i_mp_qi)
5623  else
5624  wx_cri = 0.0_rp
5625  pq(k,i_licon) = 0.0_rp
5626  pq(k,i_nicon) = 0.0_rp
5627  endif
5628  ! snow-cloud => graupel
5629  wx_crs = cfill_s*rhow/rho_g*( pi/6.0_rp*rho_g*dq_xave(k,i_mp_qs)*dq_xave(k,i_mp_qs)*dq_xave(k,i_mp_qs)/xq(k,i_mp_qs) - 1.0_rp )
5630  pq(k,i_lscon) = i_sconv2g* (pac(k,i_lsaclc2ls))/max(1.0_rp, wx_crs)
5631  pq(k,i_nscon) = i_sconv2g* pq(k,i_lscon)/xq(k,i_mp_qs)
5632  !------------------------------------------------------------------------
5633  !--- enhanced melting( due to collection-freezing of water droplets )
5634  ! originally from Rutledge and Hobbs(1984). eq.(A.21)
5635  ! if T > 273.15 then temc_p=T-273.15, else temc_p=0
5636  ! 08/05/08 [fix] T.Mitsui LHF00 => LHF0
5637  ! melting occurs around T=273K, so LHF0 is suitable both SIMPLE and EXACT,
5638  ! otherwise LHF can have sign both negative(EXACT) and positive(SIMPLE).
5639  coef_emelt = cl/lhf0*temc_p
5640  ! cloud-graupel
5641  pq(k,i_lgacm) = coef_emelt*pac(k,i_lgaclc2lg)
5642  pq(k,i_ngacm) = pq(k,i_lgacm)/xq(k,i_mp_qg)
5643  ! rain-graupel
5644  pq(k,i_lgarm) = coef_emelt*pac(k,i_lraclg2lg)
5645  pq(k,i_ngarm) = pq(k,i_lgarm)/xq(k,i_mp_qg)
5646  ! cloud-snow
5647  pq(k,i_lsacm) = coef_emelt*(pac(k,i_lsaclc2ls))
5648  pq(k,i_nsacm) = pq(k,i_lsacm)/xq(k,i_mp_qs)
5649  ! rain-snow
5650  pq(k,i_lsarm) = coef_emelt*(pac(k,i_lracls2lg_r)+pac(k,i_lracls2lg_s))
5651  pq(k,i_nsarm) = pq(k,i_lsarm)/xq(k,i_mp_qg)
5652  ! cloud-ice
5653  pq(k,i_liacm) = coef_emelt*pac(k,i_liaclc2li)
5654  pq(k,i_niacm) = pq(k,i_liacm)/xq(k,i_mp_qi)
5655  ! rain-ice
5656  pq(k,i_liarm) = coef_emelt*(pac(k,i_lracli2lg_r)+pac(k,i_lracli2lg_i))
5657  pq(k,i_niarm) = pq(k,i_liarm)/xq(k,i_mp_qg)
5658  enddo
5659 
5660  !---- for charge density
5661  if ( flg_lt ) then
5662 
5663  ! 1.c-g (X=Cloud, Y=Graupel)
5664  do k = ks, ke
5665  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
5666  pcrg2(k,i_ngacnc2ng) = pac(k,i_ngacnc2ng)*(1.0_rp-sw1) / (rhoq(k,i_nc)+sw1) * rhoq_crg(k,i_qc)
5667  end do
5668 
5669  ! 2.c-s (X=Cloud, Y=Snow)
5670  do k = ks, ke
5671  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
5672  pcrg2(k,i_nsacnc2ns) = pac(k,i_nsacnc2ns)*(1.0_rp-sw1) / (rhoq(k,i_nc)+sw1) * rhoq_crg(k,i_qc)
5673  end do
5674 
5675  ! 3.c-i (X=Cloud, Y=Cloud Ice)
5676  do k = ks, ke
5677  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
5678  pcrg2(k,i_niacnc2ni) = pac(k,i_niacnc2ni)*(1.0_rp-sw1) / (rhoq(k,i_nc)+sw1) * rhoq_crg(k,i_qc)
5679  end do
5680 
5681  ! 4.r-g
5682  do k = ks, ke
5683  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small )
5684  sw2 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small )
5685  pcrg2(k,i_nracng2ng) = pac(k,i_nracng2ng)*(1.0_rp-sw1) / (rhoq(k,i_nr)+sw1) * rhoq_crg(k,i_qr)
5686  pcrg2(k,i_nracng2nr) = pac(k,i_nracng2nr)*(1.0_rp-sw2) / (rhoq(k,i_ng)+sw2) * rhoq_crg(k,i_qg)
5687  end do
5688 
5689  ! 5.r-s
5690  do k = ks, ke
5691  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small )
5692  sw2 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
5693  pcrg2(k,i_nracns2ng_r) = pac(k,i_nracns2ng_r)*(1.0_rp-sw1) / (rhoq(k,i_nr)+sw1) * rhoq_crg(k,i_qr)
5694  pcrg2(k,i_nracns2ng_s) = pac(k,i_nracns2ng_s)*(1.0_rp-sw2) / (rhoq(k,i_ns)+sw2) * rhoq_crg(k,i_qs)
5695  end do
5696 
5697 
5698  ! 6.r-i
5699  do k = ks, ke
5700  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small )
5701  sw2 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
5702  pcrg2(k,i_nracni2ng_r) = pac(k,i_nracni2ng_r)*(1.0_rp-sw1) / (rhoq(k,i_nr)+sw1) * rhoq_crg(k,i_qr)
5703  pcrg2(k,i_nracni2ng_i) = pac(k,i_nracni2ng_i)*(1.0_rp-sw2) / (rhoq(k,i_ni)+sw2) * rhoq_crg(k,i_qi)
5704  end do
5705 
5706  ! 7.i-g
5707  do k = ks, ke
5708  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
5709  pcrg2(k,i_niacng2ng) = pac(k,i_niacng2ng)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
5710  end do
5711 
5712  ! 8.i-s
5713  do k = ks, ke
5714  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
5715  pcrg2(k,i_niacns2ns) = pac(k,i_niacns2ns)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
5716  end do
5717 
5718  ! 9.i-i
5719  do k = ks, ke
5720  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
5721  pcrg2(k,i_niacni2ns) = pac(k,i_niacni2ns)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
5722  end do
5723 
5724  ! 10.s-g
5725  do k = ks, ke
5726  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
5727  pcrg2(k,i_ngacns2ng) = pac(k,i_ngacns2ng)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
5728  end do
5729 
5730  ! 11.s-s
5731  do k = ks, ke
5732  pcrg2(k,i_nsacns2ns) = 0.0_rp ! no charge transfer between category due to the collection of the same category (snow-snow)
5733  end do
5734 
5735  ! 12.g-g
5736  do k = ks, ke
5737  pcrg2(k,i_ngacng2ng) = 0.0_rp ! no charge transfer between category due to the collection of the same category (graupel-graupel)
5738  end do
5739 
5740  !--- Partial conversion
5741  ! ice-cloud => graupel
5742  do k = ks, ke
5743  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
5744  pcrg1(k,i_nicon) = pq(k,i_nicon)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
5745  end do
5746 
5747  ! snow-cloud => graupel
5748  do k = ks, ke
5749  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
5750  pcrg1(k,i_nscon) = pq(k,i_nscon)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
5751  end do
5752 
5753  !--- enhanced melting( due to collection-freezing of water droplets )
5754  ! cloud-graupel
5755  do k = ks, ke
5756  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small )
5757  pcrg1(k,i_ngacm) = pq(k,i_ngacm)*(1.0_rp-sw1) / (rhoq(k,i_ng)+sw1) * rhoq_crg(k,i_qg)
5758  end do
5759 
5760  ! rain-graupel
5761  do k = ks, ke
5762  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small )
5763  pcrg1(k,i_ngarm) = pq(k,i_ngarm)*(1.0_rp-sw1) / (rhoq(k,i_ng)+sw1) * rhoq_crg(k,i_qg)
5764  end do
5765 
5766  ! cloud-snow
5767  do k = ks, ke
5768  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
5769  pcrg1(k,i_nsacm) = pq(k,i_nsacm)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
5770  end do
5771 
5772  ! rain-snow
5773  do k = ks, ke
5774  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
5775  pcrg1(k,i_nsarm) = pq(k,i_nsarm)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
5776  end do
5777 
5778  ! cloud-ice
5779  do k = ks, ke
5780  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
5781  pcrg1(k,i_niacm) = pq(k,i_niacm)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
5782  end do
5783 
5784  ! rain-ice
5785  do k = ks, ke
5786  sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
5787  pcrg1(k,i_niarm) = pq(k,i_niarm)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
5788  end do
5789 
5790  endif
5791 
5792  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:

◆ 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(ka,i_qc:i_qg), 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(ka,pq_max), intent(inout)  Pcrg 
)

Definition at line 5805 of file scale_atmos_phy_mp_sn14.F90.

5805  implicit none
5806 
5807  integer, intent(in) :: KA, KS, KE
5808  !
5809  real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
5810  real(RP), intent(in) :: rhoq_crg(KA,I_QC:I_QG)
5811  logical, intent(in) :: flg_lt
5812  real(RP), intent(in) :: xq(KA,HYDRO_MAX)
5813  real(RP), intent(in) :: dq_xave(KA,HYDRO_MAX)
5814  real(RP), intent(in) :: rho(KA)
5815  !
5816  real(RP), intent(inout) :: PQ(KA,PQ_MAX)
5817  real(RP), intent(inout) :: Pcrg(KA,PQ_MAX)
5818  !
5819  ! parameter for autoconversion
5820  real(RP), parameter :: kcc = 4.44e+9_rp ! collision efficiency [m3/kg2/sec]
5821  real(RP), parameter :: tau_min = 1.e-20_rp ! empirical filter by T.Mitsui
5822  real(RP), parameter :: rx_sep = 1.0_rp/x_sep ! 1/x_sep, 10/08/03 [Add] T.Mitsui
5823  !
5824  ! parameter for accretion
5825  real(RP), parameter :: kcr = 5.8_rp ! collision efficiency [m3/kg2/sec]
5826  real(RP), parameter :: thr_acc = 5.e-5_rp ! threshold for universal function original
5827  !
5828  ! parameter for self collection and collison break-up
5829  real(RP), parameter :: krr = 4.33_rp ! k_rr, S08 (35)
5830  real(RP), parameter :: kaprr = 60.7_rp ! kappa_rr, SB06(11)
5831  real(RP), parameter :: kbr = 1000._rp ! k_br, SB06(14)
5832  real(RP), parameter :: kapbr = 2.3e+3_rp ! kappa_br, SB06(15)
5833  real(RP), parameter :: dr_min = 0.35e-3_rp ! minimum diameter, SB06(13)-(15)
5834  !
5835  ! work variables
5836  real(RP) :: coef_nuc0 ! coefficient of number for Auto-conversion
5837  real(RP) :: coef_nuc1 ! mass
5838  real(RP) :: coef_aut0 ! number
5839  real(RP) :: coef_aut1 ! mass
5840  real(RP) :: lwc ! lc+lr
5841  real(RP) :: tau ! conversion ratio: qr/(qc+qr) ranges [0:1]
5842  real(RP) :: rho_fac ! factor of air density
5843  real(RP) :: psi_aut ! Universal function of Auto-conversion
5844  real(RP) :: psi_acc ! Universal function of Accretion
5845  real(RP) :: psi_brk ! Universal function of Breakup
5846  real(RP) :: ddr ! diameter difference from equilibrium
5847  !
5848  integer :: k, iqw
5849  real(RP) :: sw
5850  !
5851  coef_nuc0 = (nu(i_mp_qc)+2.0_rp)/(nu(i_mp_qc)+1.0_rp)
5852  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)
5853  coef_aut0 = -kcc*coef_nuc0
5854  coef_aut1 = -kcc/x_sep/20._rp*coef_nuc1
5855  !
5856 !OCL LOOP_FISSION_TARGET(LS)
5857  do k = ks, ke
5858  lwc = rhoq(k,i_qr) + rhoq(k,i_qc)
5859  if( lwc > xc_min )then
5860  tau = max(tau_min, rhoq(k,i_qr)/lwc)
5861  else
5862  tau = tau_min
5863  end if
5864  rho_fac = sqrt(rho_0/max(rho(k),rho_min))
5865  !
5866  ! Auto-conversion ( cloud-cloud => rain )
5867  psi_aut = 400._rp*(tau**0.7_rp)*(1.0_rp - (tau**0.7_rp))**3 ! (6) SB06
5868  pq(k,i_ncaut) = coef_aut0*rhoq(k,i_qc)*rhoq(k,i_qc)*rho_fac*rho_fac ! (9) SB06 sc+aut
5869  ! lc = lwc*(1-tau), lr = lwc*tau
5870  pq(k,i_lcaut) = coef_aut1*lwc*lwc*xq(k,i_mp_qc)*xq(k,i_mp_qc) & ! (4) SB06
5871  *((1.0_rp-tau)*(1.0_rp-tau) + psi_aut)*rho_fac*rho_fac !
5872  pq(k,i_nraut) = -rx_sep*pq(k,i_lcaut) ! (A7) SB01
5873  !--- for Charge density
5874  if (flg_lt) then
5875  sw = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
5876  pcrg(k,i_ncaut) = pq(k,i_ncaut)*(1.0_rp-sw)/(rhoq(k,i_nc)+sw)*rhoq_crg(k,i_qc)
5877  pcrg(k,i_nraut) = -pcrg(k,i_ncaut)
5878  end if
5879  !
5880  ! Accretion ( cloud-rain => rain )
5881  psi_acc =(tau/(tau+thr_acc))**4 ! (8) SB06
5882  pq(k,i_lcacc) = -kcr*rhoq(k,i_qc)*rhoq(k,i_qr)*rho_fac*psi_acc ! (7) SB06
5883  pq(k,i_ncacc) = -kcr*rhoq(k,i_nc)*rhoq(k,i_qr)*rho_fac*psi_acc ! (A6) SB01
5884  !--- for Charge density
5885  if (flg_lt) then
5886  sw = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
5887  pcrg(k,i_ncacc) = pq(k,i_ncacc)*(1.0_rp-sw)/(rhoq(k,i_nc)+sw)*rhoq(k,i_qc)
5888  end if
5889  !
5890  ! Self-collection ( rain-rain => rain )
5891  pq(k,i_nrslc) = -krr*rhoq(k,i_nr)*rhoq(k,i_qr)*rho_fac ! (A.8) SB01
5892  !
5893  ! Collisional breakup of rain
5894  ddr = min(1.e-3_rp, dq_xave(k,i_mp_qr) - dr_eq )
5895  if ( dq_xave(k,i_mp_qr) < dr_min ) then ! negligible
5896  psi_brk = -1.0_rp
5897  else if ( dq_xave(k,i_mp_qr) <= dr_eq ) then
5898  psi_brk = kbr*ddr ! (14) SB06
5899  else
5900  psi_brk = exp(kapbr*ddr) - 1.0_rp ! (15) SB06 (SB06 has a typo)
5901  end if
5902  pq(k,i_nrbrk) = - (psi_brk + 1.0_rp)*pq(k,i_nrslc) ! (13) SB06
5903  !
5904  end do
5905  !
5906  return

References scale_const::const_eps.

Referenced by atmos_phy_mp_sn14_terminal_velocity().

Here is the caller graph for this function:

◆ dep_vapor_ice_wrk()

subroutine scale_atmos_phy_mp_sn14::dep_vapor_ice_wrk ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
real(rp), dimension(ka), intent(out)  PLIdep_total,
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)  qd,
real(rp), dimension(ka), intent(in)  esi,
real(rp), dimension(ka), intent(in)  qsi,
real(rp), dimension(ka,i_qv:i_ng), intent(in)  rhoq,
real(rp), dimension(ka,hydro_max,1:2), intent(in)  vt_xave,
real(rp), dimension(ka,hydro_max), intent(in)  dq_xave,
real(rp), intent(in)  dt 
)

Definition at line 6159 of file scale_atmos_phy_mp_sn14.F90.

6159  use scale_const, only: &
6160  rvap => const_rvap, &
6161  rdry => const_rdry, &
6162  t00 => const_tem00, &
6163  lhs0 => const_lhs0, &
6164  lhf00 => const_lhf00, &
6165  pi => const_pi, &
6166  cpdry => const_cpdry, &
6167  pstd => const_pstd, &
6168  psat0 => const_psat0
6169  implicit none
6170 
6171  integer, intent(in) :: KA
6172  integer, intent(in) :: KS
6173  integer, intent(in) :: KE
6174  ! Diffusion growth or Evaporation, Sublimation
6175  real(RP), intent(out) :: PLIdep_total(KA) ! mass for cloud ice
6176  real(RP), intent(in) :: rho(KA) ! air density
6177  real(RP), intent(in) :: tem(KA) ! air temperature
6178  real(RP), intent(in) :: pre(KA) ! air pressure
6179  real(RP), intent(in) :: qd(KA) ! mixing ratio of dry air
6180  real(RP), intent(in) :: esi(KA) ! saturation vapor pressure(solid water)
6181  real(RP), intent(in) :: qsi(KA) ! saturation vapor pressure(solid water)
6182  real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
6183  ! Notice following values differ from mean terminal velocity or diameter.
6184  ! mean(vt(x)) /= vt(mean(x)) and mean(D(x)) /= D(mean(x))
6185  ! Following ones are vt(mean(x)) and D(mean(x)).
6186  real(RP), intent(in) :: vt_xave(KA,HYDRO_MAX,1:2) ! terminal velocity of mean cloud
6187  real(RP), intent(in) :: dq_xave(KA,HYDRO_MAX) !
6188  real(RP), intent(in) :: dt
6189  !
6190 ! real(RP), intent(in) :: dtime
6191 ! real(RP) :: dtime
6192  !
6193  real(RP) :: rho_lim ! limited density
6194  real(RP) :: temc_lim ! limited temperature[celsius]
6195  real(RP) :: pre_lim ! limited density
6196  real(RP) :: temc ! temperature[celsius]
6197  real(RP) :: pv ! vapor pressure
6198  real(RP) :: qv ! vapor pressure
6199  real(RP) :: ssi ! super saturation ratio(ice water)
6200  real(RP) :: nua, r_nua ! kinematic viscosity of air
6201  real(RP) :: mua ! viscosity of air
6202  real(RP) :: Kat ! thermal conductance
6203  real(RP) :: Dw ! diffusivity of water vapor
6204  real(RP) :: Dht ! diffusivity of heat
6205  real(RP) :: Gi ! diffusion factor by balance between heat and vapor
6206  real(RP) :: Gii, Gis, Gig ! for rain, ice, snow and graupel.
6207  real(RP) :: Gm ! melting factor by balance between heat and vapor
6208  real(RP) :: Nsc_r3 !
6209  !
6210  real(RP) :: Nreis_r2 !
6211  real(RP) :: Nress_r2, Nregs_r2 !
6212  real(RP) :: Nreil_r2 !
6213  real(RP) :: Nresl_r2, Nregl_r2 !
6214  real(RP) :: NscNrei_s, NscNrei_l
6215  real(RP) :: NscNres_s, NscNres_l
6216  real(RP) :: NscNreg_s, NscNreg_l
6217  real(RP) :: ventLI_s, ventLI_l
6218  real(RP) :: ventLS_s, ventLS_l
6219  real(RP) :: ventLG_s, ventLG_l
6220  real(RP) :: wti, wts, wtg
6221  real(RP), parameter :: r_14=1.0_rp/1.4_rp
6222  real(RP), parameter :: r_15=1.0_rp/1.5_rp
6223  real(RP) :: ventLI !
6224  real(RP) :: ventLS !
6225  real(RP) :: ventLG !
6226  !
6227  real(RP) :: total_dep
6228  real(RP) :: PLIdep_wrk
6229  real(RP) :: PLSdep_wrk
6230  real(RP) :: PLGdep_wrk
6231  real(RP) :: dep_limiter
6232  !
6233  real(RP), parameter :: Re_max=1.e3_rp
6234  real(RP), parameter :: Re_min=1.e-4_rp
6235  integer :: k
6236  !
6237 ! PLIdep_total(1:KS)=0.0_RP
6238 ! PLIdep_total(KE:KA)=0.0_RP
6239  plidep_total(:)=0.0_rp !! 22/10/14
6240 
6241  !
6242  do k=ks, ke
6243  temc = tem(k) - t00 ! degC
6244  temc_lim= max(temc, temc_lim_diff) !
6245  rho_lim = max(rho(k),rho_min) !
6246  qv = rhoq(k,i_qv)/rho_lim
6247  pre_lim = rho_lim*(qd(k)*rdry + qv*rvap)*(temc_lim+t00)
6248  !--------------------------------------------------------------------
6249  ! Diffusion growth part is described in detail
6250  ! by Pruppacher and Klett (1997) Sec. 13.2(liquid) and 13.3(solid)
6251  ! G:factor of thermal diffusion(1st.term) and vapor diffusion(2nd. term)
6252  ! SB06(23),(38), Lin et al(31),(52) or others
6253  ! Dw is introduced by Pruppacher and Klett(1997),(13-3)
6254  dw = 0.211e-4_rp* (((temc_lim+t00)/t00)**1.94_rp) *(pstd/pre_lim)
6255  kat = ka0 + temc_lim*dka_dt
6256  mua = mua0 + temc_lim*dmua_dt
6257  nua = mua/rho_lim
6258  r_nua = 1.0_rp/nua
6259  gi = (lhs0/kat/tem(k))*(lhs0/rvap/tem(k)-1.0_rp)+(rvap*tem(k)/dw/esi(k))
6260  ! capacities account for their surface geometries
6261  gii = 4.0_rp*pi/cap(i_mp_qi)/gi
6262  gis = 4.0_rp*pi/cap(i_mp_qs)/gi
6263  gig = 4.0_rp*pi/cap(i_mp_qg)/gi
6264  ! vent: ventilation effect( asymmetry vapor field around particles due to aerodynamic )
6265  ! SB06 (30),(31) and each coefficient is by (88),(89)
6266  nsc_r3 = (nua/dw)**(0.33333333_rp) ! (Schmidt number )^(1/3)
6267  ! Beard and Pruppacher(1971) had performed in the range [0<Re<=320],
6268  ! So here we limit Re in the range Re_max=1000, Re_min=0.0001.
6269  nreis_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qi,1)*dq_xave(k,i_mp_qi)*r_nua))) ! (Reynolds number)^(1/2) cloud ice
6270  nress_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qs,1)*dq_xave(k,i_mp_qs)*r_nua))) ! (Reynolds number)^(1/2) snow
6271  nregs_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qg,1)*dq_xave(k,i_mp_qg)*r_nua))) ! (Reynolds number)^(1/2) graupel
6272  nreil_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qi,2)*dq_xave(k,i_mp_qi)*r_nua))) ! (Reynolds number)^(1/2) cloud ice
6273  nresl_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qs,2)*dq_xave(k,i_mp_qs)*r_nua))) ! (Reynolds number)^(1/2) snow
6274  nregl_r2 = sqrt(max(re_min,min(re_max,vt_xave(k,i_mp_qg,2)*dq_xave(k,i_mp_qg)*r_nua))) ! (Reynolds number)^(1/2) graupel
6275  !
6276  nscnrei_s=nsc_r3*nreis_r2 ! small ice
6277  nscnrei_l=nsc_r3*nreil_r2 ! large ice
6278  nscnres_s=nsc_r3*nress_r2 ! small snow
6279  nscnres_l=nsc_r3*nresl_r2 ! large snow
6280  nscnreg_s=nsc_r3*nregs_r2 ! small snow
6281  nscnreg_l=nsc_r3*nregl_r2 ! large snow
6282  !
6283  ventli_s = ah_vent1(i_mp_qi,1) + bh_vent1(i_mp_qi,1)*nscnrei_s
6284  ventli_l = ah_vent1(i_mp_qi,2) + bh_vent1(i_mp_qi,2)*nscnrei_l
6285  ventls_s = ah_vent1(i_mp_qs,1) + bh_vent1(i_mp_qs,1)*nscnres_s
6286  ventls_l = ah_vent1(i_mp_qs,2) + bh_vent1(i_mp_qs,2)*nscnres_l
6287  ventlg_s = ah_vent1(i_mp_qg,1) + bh_vent1(i_mp_qg,1)*nscnreg_s
6288  ventlg_l = ah_vent1(i_mp_qg,2) + bh_vent1(i_mp_qg,2)*nscnreg_l
6289  ! branch is 1.4 for rain, snow, graupel; is 1.0 for ice (PK97, 13-60,-61,-88,-89).
6290  wti = ( min(max( nscnrei_s , 0.5_rp), 2.0_rp) -0.5_rp )*r_15 ! weighting between 1.0*0.5 and 1.0*2
6291  wts = ( min(max( nscnres_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15 ! weighting between 1.4*0.5 and 1.4*2
6292  wtg = ( min(max( nscnreg_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15 ! weighting between 1.4*0.5 and 1.4*2
6293  ! interpolation between two branches
6294  ventli = (1.0_rp-wti)*ventli_s + wti*ventli_l
6295  ventls = (1.0_rp-wts)*ventls_s + wts*ventls_l
6296  ventlg = (1.0_rp-wtg)*ventlg_s + wtg*ventlg_l
6297  ! SB06(29)
6298  ssi = qv/qsi(k) - 1.0_rp ! supersaturation
6299  plidep_wrk = gii*ssi*max(rhoq(k,i_ni),0.0_rp)*dq_xave(k,i_mp_qi)*ventli
6300  plsdep_wrk = gis*ssi*max(rhoq(k,i_ns),0.0_rp)*dq_xave(k,i_mp_qs)*ventls
6301  plgdep_wrk = gig*ssi*max(rhoq(k,i_ng),0.0_rp)*dq_xave(k,i_mp_qg)*ventlg
6302  !
6303  dep_limiter = rho(k)*(qv-qsi(k))/dt
6304  if (ssi < -1.e-30_rp)then ! unsaturated
6305  plidep_total(k) = max(plidep_wrk+plsdep_wrk+plgdep_wrk, dep_limiter)
6306  else if (ssi > 1.e-30_rp)then ! saturated
6307  plidep_total(k) = min(plidep_wrk+plsdep_wrk+plgdep_wrk, dep_limiter)
6308  else
6309  plidep_total(k) = 0.0_rp
6310  end if
6311  enddo
6312 
6313  return

References scale_const::const_cpdry, scale_const::const_lhf00, scale_const::const_lhs0, scale_const::const_pi, scale_const::const_psat0, scale_const::const_pstd, scale_const::const_rdry, scale_const::const_rvap, and scale_const::const_tem00.

Referenced by nucleation().

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

6323  implicit none
6324  !
6325  ! In this subroutine,
6326  ! We assumed surface temperature of droplets are same as environment.
6327 
6328  integer, intent(in) :: KA, KS, KE
6329 
6330  real(RP), intent(in) :: dt
6331  !
6332  real(RP), intent(in) :: tem(KA)
6333  !
6334  real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
6335  real(RP), intent(in) :: xq(KA,HYDRO_MAX)
6336  !
6337  real(RP), intent(inout):: PQ(KA,PQ_MAX)
6338  !
6339  real(RP), parameter :: temc_min = -65.0_rp
6340  real(RP), parameter :: a_het = 0.2_rp ! SB06 (44)
6341  real(RP), parameter :: b_het = 0.65_rp ! SB06 (44)
6342  !
6343  real(RP) :: coef_m2_c
6344  real(RP) :: coef_m2_r
6345  ! temperature [celsius]
6346  real(RP) :: temc, temc2, temc3, temc4
6347  ! temperature function of homegenous/heterogenous freezing
6348  real(RP) :: Jhom, Jhet, Jh(KA)
6349  real(RP) :: rdt
6350  real(RP) :: tmp
6351  !
6352  integer :: k
6353  !
6354  rdt = 1.0_rp/dt
6355  !
6356  coef_m2_c = coef_m2(i_mp_qc)
6357  coef_m2_r = coef_m2(i_mp_qr)
6358  !
6359 
6360  ! Note, xc should be limited in range[xc_min:xc_max].
6361  ! and PNChom need to be calculated by NC
6362  ! because reduction rate of Nc need to be bound by NC.
6363  ! For the same reason PLChom also bound by LC and xc.
6364  ! Basically L and N should be independent
6365  ! but average particle mass x should be in suitable range.
6366 
6367  ! Homogenous Freezing
6368  do k = ks, ke
6369  pq(k,i_lchom) = 0.0_rp
6370  pq(k,i_nchom) = 0.0_rp
6371  end do
6372 
6373  ! Heterogenous Freezing
6374  do k = ks, ke
6375  temc = max( tem(k) - t00, temc_min )
6376  ! These cause from aerosol-droplet interaction.
6377  ! Bigg(1953) formula, Khain etal.(2000) eq.(4.5), Pruppacher and Klett(1997) eq.(9-48)
6378  jhet = a_het*( exp( -b_het*temc ) - 1.0_rp )
6379  ! These cause in nature.
6380  ! Cotton and Field 2002, QJRMS. (12)
6381  if( temc < -30.0_rp ) then
6382  temc2 = temc*temc
6383  temc3 = temc*temc2
6384  temc4 = temc2*temc2
6385  jhom = 10.0_rp**(&
6386  - 243.40_rp - 14.75_rp*temc - 0.307_rp*temc2 &
6387  - 0.00287_rp*temc3 - 0.0000102_rp*temc4 ) *1.e+3_rp
6388  else if( temc < 0.0_rp) then
6389  jhom = 10._rp**(-7.63_rp-2.996_rp*(temc+30.0_rp))*1.e+3_rp
6390  else
6391  jhom = 0.0_rp
6392  jhet = 0.0_rp
6393  end if
6394  jh(k) = ( jhet + jhom ) * dt
6395  end do
6396 
6397  do k = ks, ke
6398 #if defined(NVIDIA) || defined(SX)
6399  tmp = min( xq(k,i_mp_qc)*jh(k), 1.e+3_rp) ! apply exp limiter
6400  pq(k,i_lchet) = -rdt*rhoq(k,i_qc)*( 1.0_rp - exp( -coef_m2_c*tmp ) )
6401  pq(k,i_nchet) = -rdt*rhoq(k,i_nc)*( 1.0_rp - exp( - tmp ) )
6402 
6403  tmp = min( xq(k,i_mp_qr)*jh(k), 1.e+3_rp) ! apply exp limiter
6404  pq(k,i_lrhet) = -rdt*rhoq(k,i_qr)*( 1.0_rp - exp( -coef_m2_r*tmp ) )
6405  pq(k,i_nrhet) = -rdt*rhoq(k,i_nr)*( 1.0_rp - exp( - tmp ) )
6406 #else
6407  pq(k,i_lchet) = -rdt*rhoq(k,i_qc)*( 1.0_rp - exp( -coef_m2_c*xq(k,i_mp_qc)*jh(k) ) )
6408  pq(k,i_nchet) = -rdt*rhoq(k,i_nc)*( 1.0_rp - exp( - xq(k,i_mp_qc)*jh(k) ) )
6409  pq(k,i_lrhet) = -rdt*rhoq(k,i_qr)*( 1.0_rp - exp( -coef_m2_r*xq(k,i_mp_qr)*jh(k) ) )
6410  pq(k,i_nrhet) = -rdt*rhoq(k,i_nr)*( 1.0_rp - exp( - xq(k,i_mp_qr)*jh(k) ) )
6411 #endif
6412  end do
6413 
6414  !
6415  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 ( 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,
real(rp), dimension(ka,hydro_max), intent(out)  Crs 
)

Calculate Cross Section.

Definition at line 7143 of file scale_atmos_phy_mp_sn14.F90.

7143  implicit none
7144 
7145  integer, intent(in) :: KA, KS, KE
7146  integer, intent(in) :: QA_MP
7147  real(RP), intent(in) :: QTRC0(KA,QA_MP) ! tracer mass concentration [kg/kg]
7148  real(RP), intent(in) :: DENS0(KA) ! density [kg/m3]
7149  real(RP), intent(out) :: Crs(KA,HYDRO_MAX)! Cross section [cm]
7150 
7151  ! mass concentration[kg/m3] and mean particle mass[kg]
7152  real(RP) :: xc(KA)
7153  real(RP) :: xr(KA)
7154  real(RP) :: xi(KA)
7155  real(RP) :: xs(KA)
7156  real(RP) :: xg(KA)
7157  ! radius of average mass
7158  real(RP) :: rc, rr
7159 
7160  real(RP) :: coef_Fuetal1998
7161  ! r2m_min is minimum value(moment of 1 particle with 1 micron)
7162  real(RP), parameter :: r2m_min=1.e-12_rp
7163  real(RP), parameter :: um2cm = 100.0_rp
7164 
7165  real(RP) :: limitsw, zerosw
7166  integer :: k
7167  !---------------------------------------------------------------------------
7168 
7169  ! mean particle mass[kg]
7170  do k = ks, ke
7171  xc(k) = min( xc_max, max( xc_min, dens0(k)*qtrc0(k,i_qc)/(qtrc0(k,i_nc)+nc_min) ) )
7172  xr(k) = min( xr_max, max( xr_min, dens0(k)*qtrc0(k,i_qr)/(qtrc0(k,i_nr)+nr_min) ) )
7173  xi(k) = min( xi_max, max( xi_min, dens0(k)*qtrc0(k,i_qi)/(qtrc0(k,i_ni)+ni_min) ) )
7174  xs(k) = min( xs_max, max( xs_min, dens0(k)*qtrc0(k,i_qs)/(qtrc0(k,i_ns)+ns_min) ) )
7175  xg(k) = min( xg_max, max( xg_min, dens0(k)*qtrc0(k,i_qg)/(qtrc0(k,i_ng)+ng_min) ) )
7176  enddo
7177 
7178 
7179  do k = ks, ke
7180  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)
7181  enddo
7182 
7183  do k = ks, ke
7184  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)
7185  enddo
7186 
7187  do k = ks, ke
7188  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)
7189  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)
7190  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)
7191  enddo
7192 
7193  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 110 of file scale_atmos_phy_mp_sn14.F90.

110  integer, public, parameter :: QA_MP = 11

Referenced by 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 112 of file scale_atmos_phy_mp_sn14.F90.

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

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

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

115  character(len=H_SHORT), parameter, public :: ATMOS_PHY_MP_SN14_tracer_names(QA_MP) = (/ &
116  'QV', &
117  'QC', &
118  'QR', &
119  'QI', &
120  'QS', &
121  'QG', &
122  'NC', &
123  'NR', &
124  'NI', &
125  'NS', &
126  '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 mass 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 127 of file scale_atmos_phy_mp_sn14.F90.

127  character(len=H_MID) , parameter, public :: ATMOS_PHY_MP_SN14_tracer_descriptions(QA_MP) = (/ &
128  'Ratio of Water Vapor mass to total mass (Specific humidity)', &
129  'Ratio of Cloud Water mass to total mass ', &
130  'Ratio of Rain Water mass to total mass ', &
131  'Ratio of Cloud Ice mass ratio to total mass ', &
132  'Ratio of Snow mass ratio to total mass ', &
133  'Ratio of Graupel mass ratio to total mass ', &
134  'Cloud Water Number Density ', &
135  'Rain Water Number Density ', &
136  'Cloud Ice Number Density ', &
137  'Snow Number Density ', &
138  '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 139 of file scale_atmos_phy_mp_sn14.F90.

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

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:49
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_atmos_hydrometeor::i_hr
integer, parameter, public i_hr
liquid water rain
Definition: scale_atmos_hydrometeor.F90:98
scale_atmos_hydrometeor::i_hs
integer, parameter, public i_hs
snow
Definition: scale_atmos_hydrometeor.F90:100
scale_const::const_lhs00
real(rp), public const_lhs00
latent heat of sublimation at 0K [J/kg]
Definition: scale_const.F90:85
scale_const::const_epsvap
real(rp), public const_epsvap
Rdry / Rvap.
Definition: scale_const.F90:75
scale_const::const_lhv00
real(rp), public const_lhv00
latent heat of vaporizaion at 0K [J/kg]
Definition: scale_const.F90:83
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:68
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:91
scale_atmos_hydrometeor::i_hh
integer, parameter, public i_hh
hail
Definition: scale_atmos_hydrometeor.F90:102
scale_file_history
module file_history
Definition: scale_file_history.F90:15
scale_const::const_cpvap
real(rp), parameter, public const_cpvap
specific heat (water vapor, constant pressure) [J/kg/K]
Definition: scale_const.F90:69
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_atmos_hydrometeor::cv_vapor
real(rp), public cv_vapor
CV for vapor [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:149
scale_atmos_hydrometeor::i_hi
integer, parameter, public i_hi
ice water cloud
Definition: scale_atmos_hydrometeor.F90:99
scale_const::const_pstd
real(rp), public const_pstd
standard pressure [Pa]
Definition: scale_const.F90:96
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:60
scale_specfunc
module SPECFUNC
Definition: scale_specfunc.F90:14
scale_const::const_psat0
real(rp), parameter, public const_psat0
saturate pressure of water vapor at 0C [Pa]
Definition: scale_const.F90:88
scale_const::const_pi
real(rp), parameter, public const_pi
pi
Definition: scale_const.F90:32
scale_atmos_hydrometeor::i_hc
integer, parameter, public i_hc
liquid water cloud
Definition: scale_atmos_hydrometeor.F90:97
scale_const::const_ci
real(rp), parameter, public const_ci
specific heat (ice) [J/kg/K]
Definition: scale_const.F90:72
scale_const::const_tem00
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
Definition: scale_const.F90:99
scale_const::const_cl
real(rp), parameter, public const_cl
specific heat (liquid water) [J/kg/K]
Definition: scale_const.F90:71
scale_const::const_cvvap
real(rp), public const_cvvap
specific heat (water vapor, constant volume) [J/kg/K]
Definition: scale_const.F90:70
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:59
scale_file_history::file_history_reg
subroutine, public file_history_reg(name, desc, unit, itemid, standard_name, ndims, dim_type, cell_measures, fill_halo)
Register/Append variable to history file.
Definition: scale_file_history.F90:685
scale_const::const_lhs0
real(rp), parameter, public const_lhs0
latent heat of sublimation at 0C [J/kg]
Definition: scale_const.F90:84
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:43
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
scale_atmos_hydrometeor::cv_ice
real(rp), public cv_ice
CV for ice [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:153
scale_const::const_lhf00
real(rp), public const_lhf00
latent heat of fusion at 0K [J/kg]
Definition: scale_const.F90:87
scale_atmos_hydrometeor::i_hg
integer, parameter, public i_hg
graupel
Definition: scale_atmos_hydrometeor.F90:101
Generated by   doxygen 1.8.17