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

  name	type	default value	comment
  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

  name	type	default value	comment
  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

  name	type	default value	comment
  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

  name	type	default value	comment
  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

  name	type	default value	comment
  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

  name	type	default value	comment
  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

  name	type	default value	comment
  OPT_FIX_TAUCND_C	logical	.false.
  FAC_CNDC	real(RP)	1.0_RP

  
  

History Output

name	description	unit	variable
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 714 of file scale_atmos_phy_mp_sn14.F90.

    use scale_prc, only: &
       prc_abort
    use scale_file_history, only: &
       file_history_reg
    implicit none
    integer, intent(in) :: KA
    integer, intent(in) :: IA
    integer, intent(in) :: JA
 
    namelist / param_atmos_phy_mp_sn14 / &
       mp_doautoconversion, &
       mp_ssw_lim,          &
       mp_couple_aerosol
 
    integer :: ip
    integer :: ierr
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_MP_sn14_setup",*) 'Setup'
    log_info("ATMOS_PHY_MP_sn14_setup",*) 'Seiki and Nakajima (2014) 2-moment bulk 6 category'
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_mp_sn14,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("ATMOS_PHY_MP_sn14_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_PHY_MP_sn14_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SN14. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_phy_mp_sn14)
 
    call mp_sn14_init
 
    allocate(nc_uplim_d(1,ia,ja))
    nc_uplim_d(:,:,:) = 150.e6_rp
 
    hist_max = 0
    hist_idx(:) = -1
    do ip = 1, w_nmax
       call file_history_reg( w_name(ip),  'individual tendency term in SN14', 'kg/kg/s', &
                              hist_id(ip) )
       if ( hist_id(ip) > 0 ) then
          hist_max = hist_max + 1
          hist_idx(ip) = hist_max
       end if
    enddo
    allocate( w3d(ka,ia,ja,hist_max) )
    w3d(:,:,:,:) = 0.0_rp
 
    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 772 of file scale_atmos_phy_mp_sn14.F90.

 
    deallocate(nc_uplim_d)
    deallocate(w3d)
 
    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 784 of file scale_atmos_phy_mp_sn14.F90.

    implicit none
 
    integer, intent(in) :: KA, KS, KE
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE
 
    real(RP), intent(in) :: DENS     (KA,IA,JA)
    real(RP), intent(in) :: W        (KA,IA,JA)
    real(RP), intent(in) :: QTRC     (KA,IA,JA,QA_MP)
    real(RP), intent(in) :: PRES(KA,IA,JA)
    real(RP), intent(in) :: TEMP(KA,IA,JA)
    real(RP), intent(in) :: Qdry(KA,IA,JA)
    real(RP), intent(in) :: CPtot(KA,IA,JA)
    real(RP), intent(in) :: CVtot(KA,IA,JA)
    real(RP), intent(in) :: CCN      (KA,IA,JA)
    real(DP), intent(in) :: dt
    real(RP), intent(in) :: cz(  KA,IA,JA)
    real(RP), intent(in) :: fz(0:KA,IA,JA)
 
    real(RP), intent(out) :: RHOQ_t   (KA,IA,JA,QA_MP)
    real(RP), intent(out) :: RHOE_t   (KA,IA,JA)
    real(RP), intent(out) :: CPtot_t(KA,IA,JA)
    real(RP), intent(out) :: CVtot_t(KA,IA,JA)
    real(RP), intent(out) :: EVAPORATE(KA,IA,JA)   !--- number of evaporated cloud [/m3]
 
    ! Optional for Lightning
    logical,  intent(in), optional :: flg_lt
    real(RP), intent(in), optional :: d0_crg, v0_crg
    real(RP), intent(in), optional :: dqcrg(KA,IA,JA)
    real(RP), intent(in), optional :: beta_crg(KA,IA,JA)
    real(RP), intent(in), optional :: QTRC_crg(KA,IA,JA,HYDRO_MAX)
    real(RP), intent(out), optional :: QSPLT_in(KA,IA,JA,3)
    real(RP), intent(out), optional :: Sarea(KA,IA,JA,HYDRO_MAX)
    real(RP), intent(out), optional :: RHOQcrg_t(KA,IA,JA,HYDRO_MAX)
    !---------------------------------------------------------------------------
 
    log_progress(*) 'atmosphere / physics / microphysics / SN14'
 
#ifdef PROFILE_FIPP
    call fipp_start()
#endif
 
    !##### MP Main #####
    call mp_sn14 ( &
       ka, ks, ke, ia, is, ie, ja, js, je, &
       dens(:,:,:), w(:,:,:), qtrc(:,:,:,:), pres(:,:,:), temp(:,:,:), & ! (in)
       qdry(:,:,:), cptot(:,:,:), cvtot(:,:,:), ccn(:,:,:),            & ! (in)
       real(dt,RP), cz(:,:,:), fz(:,:,:),                              & ! (in)
       RHOQ_t(:,:,:,:), RHOE_t(:,:,:), CPtot_t(:,:,:), CVtot_t(:,:,:), & ! (out)
       EVAPORATE(:,:,:),                                               & ! (out)
       flg_lt, d0_crg, v0_crg, dqcrg(:,:,:), beta_crg(:,:,:),          & ! (optional in)
       QTRC_crg(:,:,:,:),                                              & ! (optional in)
       QSPLT_in(:,:,:,:), Sarea(:,:,:,:), RHOQcrg_t(:,:,:,:)           ) ! (optional out)
 
#ifdef PROFILE_FIPP
    call fipp_stop()
#endif
 
    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 875 of file scale_atmos_phy_mp_sn14.F90.

    implicit none
    integer, intent(in) :: KA, KS, KE
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE
 
    real(RP), intent(in)  :: QTRC   (KA,IA,JA,QA_MP-1)
    real(RP), intent(in)  :: mask_criterion
 
    real(RP), intent(out) :: cldfrac(KA,IA,JA)
 
    real(RP) :: qhydro
    integer  :: k, i, j, iq
    !---------------------------------------------------------------------------
 
    !$omp parallel do &
    !$omp private(qhydro)
    do j  = js, je
    do i  = is, ie
    do k  = ks, ke
       qhydro = 0.0_rp
       do iq = i_mp_qc, i_mp_qg
          qhydro = qhydro + qtrc(k,i,j,iq)
       enddo
       cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
    enddo
    enddo
    enddo
 
    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 914 of file scale_atmos_phy_mp_sn14.F90.

    use scale_atmos_hydrometeor, only: &
       n_hyd, &
       i_hc,  &
       i_hr,  &
       i_hi,  &
       i_hs,  &
       i_hg,  &
       i_hh
    implicit none
    integer, intent(in) :: KA, KS, KE
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE
 
    real(RP), intent(in)  :: DENS0(KA,IA,JA)        ! density                   [kg/m3]
    real(RP), intent(in)  :: TEMP0(KA,IA,JA)        ! temperature               [K]
    real(RP), intent(in)  :: QTRC0(KA,IA,JA,I_QC:I_NG)     ! tracer mass concentration [kg/kg]
 
    real(RP), intent(out) :: Re   (KA,IA,JA,N_HYD) ! effective radius          [cm]
 
    ! mass concentration[kg/m3] and mean particle mass[kg]
    real(RP) :: xc(KA)
    real(RP) :: xr(KA)
    real(RP) :: xi(KA)
    real(RP) :: xs(KA)
    real(RP) :: xg(KA)
    ! diameter of average mass[kg/m3]
    real(RP) :: dc_ave(KA)
    real(RP) :: dr_ave(KA)
    ! radius of average mass
    real(RP) :: rc, rr
    ! 2nd. and 3rd. order moment of DSD
    real(RP) :: ri2m(KA), ri3m(KA)
    real(RP) :: rs2m(KA), rs3m(KA)
    real(RP) :: rg2m(KA), rg3m(KA)
 
    real(RP), parameter :: coef_Fuetal1998 = 3.0_rp / (4.0_rp*rhoi)
 
    ! r2m_min is minimum value(moment of 1 particle with 1 micron)
    real(RP), parameter :: r2m_min=1.e-12_rp
    real(RP), parameter :: um2cm = 100.0_rp
 
    real(RP) :: limitsw, zerosw
    integer :: k, i, j
    !---------------------------------------------------------------------------
 
    !$omp parallel do &
    !$omp private(xc,xr,xi,xs,xg,dc_ave,dr_ave,rc,rr,ri2m,ri3m,rs2m,rs3m,rg2m,rg3m, &
    !$omp         limitsw,zerosw)
    do j  = js, je
    do i  = is, ie
 
       ! mean particle mass[kg]
       do k  = ks, ke
          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) ) )
          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) ) )
          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) ) )
          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) ) )
          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) ) )
       enddo
 
       ! diameter of average mass : SB06 eq.(32)
       do k = ks, ke
          dc_ave(k) = a_m(i_mp_qc) * xc(k)**b_m(i_mp_qc)
          dr_ave(k) = a_m(i_mp_qr) * xr(k)**b_m(i_mp_qr)
       enddo
 
       ! cloud effective radius
       do k = ks, ke
          rc = 0.5_rp * dc_ave(k)
          limitsw = 0.5_rp + sign(0.5_rp, rc-rmin_re )
          re(k,i,j,i_hc) = coef_re(i_mp_qc) * rc * limitsw * um2cm
       enddo
 
       ! rain effective radius
       do k = ks, ke
          rr = 0.5_rp * dr_ave(k)
          limitsw = 0.5_rp + sign(0.5_rp, rr-rmin_re )
          re(k,i,j,i_hr) = coef_re(i_mp_qr) * rr * limitsw * um2cm
       enddo
 
       do k = ks, ke
          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)
          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)
          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)
       enddo
 
       ! Fu(1996), eq.(3.11) or Fu et al.(1998), eq.(2.5)
       do k = ks, ke
          ri3m(k) = coef_fuetal1998 * qtrc0(k,i,j,i_ni) * xi(k)
          rs3m(k) = coef_fuetal1998 * qtrc0(k,i,j,i_ns) * xs(k)
          rg3m(k) = coef_fuetal1998 * qtrc0(k,i,j,i_ng) * xg(k)
       enddo
 
       ! ice effective radius
       do k = ks, ke
          zerosw = 0.5_rp - sign(0.5_rp, ri2m(k) - r2m_min )
          re(k,i,j,i_hi) = ri3m(k) / ( ri2m(k) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
       enddo
 
       ! snow effective radius
       do k = ks, ke
          zerosw = 0.5_rp - sign(0.5_rp, rs2m(k) - r2m_min )
          re(k,i,j,i_hs) = rs3m(k) / ( rs2m(k) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
       enddo
 
       ! graupel effective radius
       do k = ks, ke
          zerosw = 0.5_rp - sign(0.5_rp, rg2m(k) - r2m_min )
          re(k,i,j,i_hg) = rg3m(k) / ( rg2m(k) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
       enddo
 
       do k = ks, ke
          re(k,i,j,i_hh) = 0.0_rp
       end do
 
    enddo
    enddo
 
 
    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 1043 of file scale_atmos_phy_mp_sn14.F90.

    use scale_atmos_hydrometeor, only: &
       n_hyd, &
       i_hc,  &
       i_hr,  &
       i_hi,  &
       i_hs,  &
       i_hg,  &
       i_hh
    implicit none
    integer, intent(in) :: KA, KS, KE
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE
 
    real(RP), intent(in)  :: QTRC0(KA,IA,JA,QA_MP-1) ! tracer mass concentration [kg/kg]
 
    real(RP), intent(out) :: Qe   (KA,IA,JA,N_HYD)   ! mixing ratio of each cateory [kg/kg]
 
    integer :: k, i, j
    !---------------------------------------------------------------------------
 
!OCL XFILL
    !$omp parallel do
    do j = js, je
    do i = is, ie
    do k = ks, ke
       qe(k,i,j,i_hc) = qtrc0(k,i,j,i_mp_qc)
       qe(k,i,j,i_hr) = qtrc0(k,i,j,i_mp_qr)
       qe(k,i,j,i_hi) = qtrc0(k,i,j,i_mp_qi)
       qe(k,i,j,i_hs) = qtrc0(k,i,j,i_mp_qs)
       qe(k,i,j,i_hg) = qtrc0(k,i,j,i_mp_qg)
       qe(k,i,j,i_hh) = 0.0_rp
    end do
    end do
    end do
 
    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 1086 of file scale_atmos_phy_mp_sn14.F90.

    use scale_atmos_hydrometeor, only: &
       n_hyd, &
       i_hc,  &
       i_hr,  &
       i_hi,  &
       i_hs,  &
       i_hg,  &
       i_hh
    implicit none
    integer, intent(in) :: KA, KS, KE
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE
 
    real(RP), intent(in)  :: QTRC0(KA,IA,JA,QA_MP-1) ! tracer mass concentration [kg/kg]
 
    real(RP), intent(out) :: Ne   (KA,IA,JA,N_HYD)   ! number density of each cateory [1/m3]
 
    integer :: k, i, j
    !---------------------------------------------------------------------------
 
!OCL XFILL
    !$omp parallel do
    do j = js, je
    do i = is, ie
    do k = ks, ke
       ne(k,i,j,i_hc) = qtrc0(k,i,j,i_mp_nc)
       ne(k,i,j,i_hr) = qtrc0(k,i,j,i_mp_nr)
       ne(k,i,j,i_hi) = qtrc0(k,i,j,i_mp_ni)
       ne(k,i,j,i_hs) = qtrc0(k,i,j,i_mp_ns)
       ne(k,i,j,i_hg) = qtrc0(k,i,j,i_mp_ng)
       ne(k,i,j,i_hh) = 0.0_rp
    end do
    end do
    end do
 
    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 1128 of file scale_atmos_phy_mp_sn14.F90.

    use scale_const, only: &
       pi    => const_pi, &
       undef => const_undef
    use scale_atmos_hydrometeor, only: &
       n_hyd, &
       i_hc,  &
       i_hr,  &
       i_hi,  &
       i_hs,  &
       i_hg,  &
       i_hh
    implicit none
    integer, intent(in) :: KA, KS, KE
    integer, intent(in) :: IA, IS, IE
    integer, intent(in) :: JA, JS, JE
 
    real(RP), intent(in) :: Qe(KA,IA,JA,N_HYD) ! mass ratio of each cateory [kg/kg]
 
    real(RP), intent(out) :: QTRC(KA,IA,JA,QA_MP-1) ! tracer mass concentration [kg/kg]
 
    real(RP), intent(in), optional :: QNUM(KA,IA,JA,N_HYD)
 
    real(RP), parameter :: Dc   =  20.e-6_rp ! typical particle diameter for cloud  [m]
    real(RP), parameter :: Dr   = 200.e-6_rp ! typical particle diameter for rain   [m]
    real(RP), parameter :: Di   =  80.e-6_rp ! typical particle diameter for ice    [m]
    real(RP), parameter :: Ds   =  80.e-6_rp ! typical particle diameter for snow   [m]
    real(RP), parameter :: Dg   = 200.e-6_rp ! typical particle diameter for grapel [m]
    real(RP), parameter :: b    =     3.0_rp ! assume spherical form
 
    real(RP) :: piov6
 
    integer :: k, i, j
    !---------------------------------------------------------------------------
 
 
!OCL XFILL
    !$omp parallel do
    do j = js, je
    do i = is, ie
    do k = ks, ke
       qtrc(k,i,j,i_mp_qc) = qe(k,i,j,i_hc)
    end do
    end do
    end do
 
!OCL XFILL
    !$omp parallel do
    do j = js, je
    do i = is, ie
    do k = ks, ke
       qtrc(k,i,j,i_mp_qr) = qe(k,i,j,i_hr)
    end do
    end do
    end do
 
!OCL XFILL
    !$omp parallel do
    do j = js, je
    do i = is, ie
    do k = ks, ke
       qtrc(k,i,j,i_mp_qi) = qe(k,i,j,i_hi)
    end do
    end do
    end do
 
!OCL XFILL
    !$omp parallel do
    do j = js, je
    do i = is, ie
    do k = ks, ke
       qtrc(k,i,j,i_mp_qs) = qe(k,i,j,i_hs)
    end do
    end do
    end do
 
!OCL XFILL
    !$omp parallel do
    do j = js, je
    do i = is, ie
    do k = ks, ke
       qtrc(k,i,j,i_mp_qg) = qe(k,i,j,i_hg) + qe(k,i,j,i_hh)
    end do
    end do
    end do
 
    piov6 = pi / 6.0_rp
 
    if ( present(qnum) ) then
 
!OCL XFILL
       !$omp parallel do
       do j = js, je
       do i = is, ie
       do k = ks, ke
          if ( qnum(k,i,j,i_hc) .ne. undef ) then
             qtrc(k,i,j,i_mp_nc) = qnum(k,i,j,i_hc)
          else
             qtrc(k,i,j,i_mp_nc) = qtrc(k,i,j,i_mp_qc) / ( (piov6*rhow) * dc**b )
          end if
       end do
       end do
       end do
 
!OCL XFILL
       !$omp parallel do
       do j = js, je
       do i = is, ie
       do k = ks, ke
          if ( qnum(k,i,j,i_hr) .ne. undef ) then
             qtrc(k,i,j,i_mp_nr) = qnum(k,i,j,i_hr)
          else
             qtrc(k,i,j,i_mp_nr) = qtrc(k,i,j,i_mp_qr) / ( (piov6*rhow) * dr**b )
          end if
       end do
       end do
       end do
 
!OCL XFILL
       !$omp parallel do
       do j = js, je
       do i = is, ie
       do k = ks, ke
          if ( qnum(k,i,j,i_hi) .ne. undef ) then
             qtrc(k,i,j,i_mp_ni) = qnum(k,i,j,i_hi)
          else
             qtrc(k,i,j,i_mp_ni) = qtrc(k,i,j,i_mp_qi) / ( (piov6*rhof) * di**b )
          end if
       end do
       end do
       end do
 
!OCL XFILL
       !$omp parallel do
       do j = js, je
       do i = is, ie
       do k = ks, ke
          if ( qnum(k,i,j,i_hs) .ne. undef ) then
             qtrc(k,i,j,i_mp_ns) = qnum(k,i,j,i_hs)
          else
             qtrc(k,i,j,i_mp_ns) = qtrc(k,i,j,i_mp_qs) / ( (piov6*rhof) * ds**b )
          end if
       end do
       end do
       end do
 
!OCL XFILL
       !$omp parallel do
       do j = js, je
       do i = is, ie
       do k = ks, ke
          if ( qnum(k,i,j,i_hg) .ne. undef ) then
             if ( qnum(k,i,j,i_hh) .ne. undef ) then
                qtrc(k,i,j,i_mp_ng) = qnum(k,i,j,i_hg) + qnum(k,i,j,i_hh)
             else
                qtrc(k,i,j,i_mp_ng) = qnum(k,i,j,i_hg)
             end if
          else
             qtrc(k,i,j,i_mp_ng) = qtrc(k,i,j,i_mp_qg) / ( (piov6*rhog) * dg**b )
          end if
       end do
       end do
       end do
 
    else
 
!OCL XFILL
       !$omp parallel do
       do j = js, je
       do i = is, ie
       do k = ks, ke
          qtrc(k,i,j,i_mp_nc) = qtrc(k,i,j,i_mp_qc) / ( (piov6*rhow) * dc**b )
       end do
       end do
       end do
 
!OCL XFILL
       !$omp parallel do
       do j = js, je
       do i = is, ie
       do k = ks, ke
          qtrc(k,i,j,i_mp_nr) = qtrc(k,i,j,i_mp_qr) / ( (piov6*rhow) * dr**b )
       end do
       end do
       end do
 
!OCL XFILL
       !$omp parallel do
       do j = js, je
       do i = is, ie
       do k = ks, ke
          qtrc(k,i,j,i_mp_ni) = qtrc(k,i,j,i_mp_qi) / ( (piov6*rhof) * di**b )
       end do
       end do
       end do
 
!OCL XFILL
       !$omp parallel do
       do j = js, je
       do i = is, ie
       do k = ks, ke
          qtrc(k,i,j,i_mp_ns) = qtrc(k,i,j,i_mp_qs) / ( (piov6*rhof) * ds**b )
       end do
       end do
       end do
 
!OCL XFILL
       !$omp parallel do
       do j = js, je
       do i = is, ie
       do k = ks, ke
          qtrc(k,i,j,i_mp_ng) = qtrc(k,i,j,i_mp_qg) / ( (piov6*rhog) * dg**b )
       end do
       end do
       end do
 
    end if
 
    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 1358 of file scale_atmos_phy_mp_sn14.F90.

    !$acc routine vector
    use scale_const, only: &
       const_undef
    implicit none
 
    integer, intent(in) :: KA, KS, KE
 
    real(RP), intent(in)  :: RHOQ(KA,I_QC:I_NG) ! rho * q
    real(RP), intent(in)  :: DENS(KA)    ! rho
    real(RP), intent(in)  :: TEMP(KA)    ! temperature
    real(RP), intent(in)  :: PRES(KA)    ! pressure
 
    real(RP), intent(out) :: vterm(KA,QA_MP-1) ! terminal velocity of cloud mass
 
    real(RP) :: xq, log_xq ! average mass of 1 particle( mass/number )
 
    real(RP) :: rhofac   ! density factor for terminal velocity
    real(RP) :: rhofac_q(KA), log_rhofac_q
 
    real(RP) :: rlambdar(KA) ! work for diagnosis of Rain DSD ( Seifert, 2008 )
    real(RP) :: mud_r
    real(RP) :: dq, log_dq   ! weigthed diameter.   Improved Rogers etal. (1993) formula by T.Mitsui
 
    real(RP) :: weight ! weighting coefficient for 2-branches is determined by ratio between 0.745mm and weighted diameter.  SB06 Table.1
    real(RP) :: velq_s ! terminal velocity for small branch of Rogers formula
    real(RP) :: velq_l ! terminal velocity for large branch of Rogers formula
 
    real(RP) :: tmp
    integer :: k, i, j, iq
    !---------------------------------------------------------------------------
 
    ! QC, NC
    do k = ks, ke
       rhofac = rho_0 / max( dens(k), rho_min )
 
       log_rhofac_q = log(rhofac) * gamma_v(i_mp_qc)
       log_xq = log( max( xc_min, min( xc_max, rhoq(k,i_qc) / ( rhoq(k,i_nc) + nc_min ) ) ) )
 
       vterm(k,i_mp_qc) = - exp( log_rhofac_q + log_coef_vt1(i_mp_qc,1) + log_xq * beta_v(i_mp_qc,1) )
 
       vterm(k,i_mp_nc) = - exp( log_rhofac_q + log_coef_vt0(i_mp_qc,1) + log_xq * beta_vn(i_mp_qc,1) )
    end do
 
    ! QR, NR
    mud_r = 3.0_rp * nu(i_mp_qr) + 2.0_rp
    do k = ks, ke
       rhofac = rho_0 / max( dens(k), rho_min )
       rhofac_q(k) = rhofac**gamma_v(i_mp_qr)
    end do
    do k = ks, ke
       xq = max( xr_min, min( xr_max, rhoq(k,i_qr) / ( rhoq(k,i_nr) + nr_min ) ) )
 
       rlambdar(k) = a_m(i_mp_qr) * xq**b_m(i_mp_qr) &
                   * ( (mud_r+3.0_rp) * (mud_r+2.0_rp) * (mud_r+1.0_rp) )**(-0.333333333_rp)
    end do
!OCL LOOP_FISSION_TARGET(LS)
    do k = ks, ke
       dq = ( 4.0_rp + mud_r ) * rlambdar(k) ! D^(3)+mu weighted mean diameter
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )
       velq_s = coef_vtr_ar2 * dq &
              * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar(k) )**(-5.0_rp-mud_r) )
       velq_l = coef_vtr_ar1 &
              - coef_vtr_br1 * ( 1.0_rp + coef_vtr_cr1*rlambdar(k) )**(-4.0_rp-mud_r)
       vterm(k,i_mp_qr) = -rhofac_q(k) * ( velq_l * (          weight ) &
                                         + velq_s * ( 1.0_rp - weight ) )
    end do
!OCL LOOP_FISSION_TARGET(LS)
    do k = ks, ke
       dq = ( 1.0_rp + mud_r ) * rlambdar(k) ! D^(0)+mu weighted mean diameter
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )
       velq_s = coef_vtr_ar2 * dq &
              * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar(k) )**(-2.0_rp-mud_r) )
       velq_l = coef_vtr_ar1 &
              - coef_vtr_br1 * ( 1.0_rp + coef_vtr_cr1*rlambdar(k) )**(-1.0_rp-mud_r)
       vterm(k,i_mp_nr) = -rhofac_q(k) * ( velq_l * (          weight ) &
                                         + velq_s * ( 1.0_rp - weight ) )
    end do
 
    do k = ks, ke
       rhofac_q(k) = exp( log( pres(k)/pre0_vt ) * a_pre0_vt + log( temp(k)/tem0_vt ) * a_tem0_vt )
    end do
 
    ! QI, NI
!OCL LOOP_FISSION_TARGET(LS)
    do k = ks, ke
       log_xq = log( max( xi_min, min( xi_max, rhoq(k,i_qi) / ( rhoq(k,i_ni) + ni_min ) ) ) )
 
       tmp = log_a_m(i_mp_qi) + log_xq * b_m(i_mp_qi)
       log_dq = log_coef_dave_l(i_mp_qi) + tmp
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log_dq - log_d0_li ) ) )
 
       velq_s = exp( log_coef_vt1(i_mp_qi,1) + log_xq * beta_v(i_mp_qi,1) )
       velq_l = exp( log_coef_vt1(i_mp_qi,2) + log_xq * beta_v(i_mp_qi,2) )
       vterm(k,i_mp_qi) = - rhofac_q(k) * ( velq_l * (          weight ) &
                                          + velq_s * ( 1.0_rp - weight ) )
 
       log_dq = log_coef_dave_n(i_mp_qi) + tmp
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log_dq - log_d0_ni ) ) )
 
       velq_s = exp( log_coef_vt0(i_mp_qi,1) + log_xq * beta_vn(i_mp_qi,1) )
       velq_l = exp( log_coef_vt0(i_mp_qi,2) + log_xq * beta_vn(i_mp_qi,2) )
       vterm(k,i_mp_ni) = - rhofac_q(k) * ( velq_l * (          weight ) &
                                          + velq_s * ( 1.0_rp - weight ) )
    end do
 
    ! QS, NS
!OCL LOOP_FISSION_TARGET(LS)
    do k = ks, ke
       log_xq = log( max( xs_min, min( xs_max, rhoq(k,i_qs) / ( rhoq(k,i_ns) + ns_min ) ) ) )
 
       tmp = log_a_m(i_mp_qs) + log_xq * b_m(i_mp_qs)
       log_dq = log_coef_dave_l(i_mp_qs) + tmp
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log_dq - log_d0_ls ) ) )
 
       velq_s = exp( log_coef_vt1(i_mp_qs,1) + log_xq * beta_v(i_mp_qs,1) )
       velq_l = exp( log_coef_vt1(i_mp_qs,2) + log_xq * beta_v(i_mp_qs,2) )
       vterm(k,i_mp_qs) = - rhofac_q(k) * ( velq_l * (          weight ) &
                                          + velq_s * ( 1.0_rp - weight ) )
 
       log_dq = log_coef_dave_n(i_mp_qs) + tmp
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log_dq - log_d0_ns ) ) )
 
       velq_s = exp( log_coef_vt0(i_mp_qs,1) + log_xq * beta_vn(i_mp_qs,1) )
       velq_l = exp( log_coef_vt0(i_mp_qs,2) + log_xq * beta_vn(i_mp_qs,2) )
       vterm(k,i_mp_ns) = - rhofac_q(k) * ( velq_l * (          weight ) &
                                          + velq_s * ( 1.0_rp - weight ) )
    end do
 
    ! QG, NG
!OCL LOOP_FISSION_TARGET(LS)
    do k = ks, ke
       log_xq = log( max( xg_min, min( xg_max, rhoq(k,i_qg) / ( rhoq(k,i_ng) + ng_min ) ) ) )
 
       tmp = log_a_m(i_mp_qg) + log_xq * b_m(i_mp_qg)
       log_dq = log_coef_dave_l(i_mp_qg) + tmp
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log_dq - log_d0_lg ) ) )
 
       velq_s = exp( log_coef_vt1(i_mp_qg,1) + log_xq * beta_v(i_mp_qg,1) )
       velq_l = exp( log_coef_vt1(i_mp_qg,2) + log_xq * beta_v(i_mp_qg,2) )
       vterm(k,i_mp_qg) = - rhofac_q(k) * ( velq_l * (          weight ) &
                                          + velq_s * ( 1.0_rp - weight ) )
 
       log_dq = log_coef_dave_n(i_mp_qg) + tmp
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + log_dq - log_d0_ng ) ) )
 
       velq_s = exp( log_coef_vt0(i_mp_qg,1) + log_xq * beta_vn(i_mp_qg,1) )
       velq_l = exp( log_coef_vt0(i_mp_qg,2) + log_xq * beta_vn(i_mp_qg,2) )
       vterm(k,i_mp_ng) = - rhofac_q(k) * ( velq_l * (          weight ) &
                                          + velq_s * ( 1.0_rp - weight ) )
    enddo
 
    do iq = 1, qa_mp-1
       vterm(1:ks-2 ,iq) = const_undef
       vterm(ks-1   ,iq) = vterm(ks,iq)
       vterm(ke+1:ka,iq) = const_undef
    enddo
 
    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 3476 of file scale_atmos_phy_mp_sn14.F90.

    use scale_prc, only: &
       prc_myrank
    implicit none
    integer, intent(in) :: KA, KS, KE
 
    integer,  intent(in) :: point
    integer,  intent(in) :: i, j
    real(RP), intent(in) :: tem(KA)
    real(RP), intent(in) :: rho(KA)
    real(RP), intent(in) :: pre(KA)
    real(RP), intent(in) :: qv (KA)
 
    integer :: k
    !---------------------------------------------------------------------------
 
    do k = ks, ke
       if (      tem(k) < tem_min &
            .OR. rho(k) < rho_min &
            .OR. pre(k) < 1.0_rp    ) then
 
          log_info("ATMOS_PHY_MP_SN14_debug_tem_kij",'(A,I3,A,4(F16.5),3(I6))') &
          "point: ", point, " low tem,rho,pre:", tem(k), rho(k), pre(k), qv(k), k, i, j, prc_myrank
       endif
    enddo
 
    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 3509 of file scale_atmos_phy_mp_sn14.F90.

    use scale_prc, only: &
       prc_abort
    use scale_atmos_saturation, only: &
       moist_psat_liq       => atmos_saturation_psat_liq, &
       moist_psat_ice       => atmos_saturation_psat_ice,   &
       moist_pres2qsat_liq  => atmos_saturation_pres2qsat_liq, &
       moist_pres2qsat_ice  => atmos_saturation_pres2qsat_ice,   &
       moist_dqsi_dtem_dens => atmos_saturation_dqs_dtem_dens_liq, &
       moist_dqs_dtem_dpre_ice => atmos_saturation_dqs_dtem_dpre_ice
    use scale_atmos_hydrometeor, only: &
       cv_vapor, &
       cv_ice
    implicit none
 
    integer, intent(in) :: KA, KS, KE
 
    real(RP), intent(in)  :: cz(  KA)
    real(RP), intent(in)  :: fz(0:KA)
    real(RP), intent(in)  :: w   (KA)    ! w of full level
    real(RP), intent(in)  :: rho (KA)
    real(RP), intent(in)  :: tem (KA)
    real(RP), intent(in)  :: pre (KA)
    real(RP), intent(in)  :: qdry(KA)
    !
    real(RP), intent(in)  :: rhoq(KA,I_QV:I_NG)
    !
    real(RP), intent(in) ::  cpa(KA)
    real(RP), intent(in) :: dTdt_rad(KA)
    real(RP), intent(in) :: qke(KA)
    real(RP), intent(in) :: dt
    real(RP), intent(in) :: CCN(KA)
    real(RP), intent(in) :: nc_uplim_d
    !
    real(RP), intent(out) :: PQ(KA,PQ_MAX)
    !
    !
!    real(RP) :: c_ccn_map(1)   ! c_ccn horizontal distribution
!    real(RP) :: kappa_map(1)   ! kappa horizontal distribution
!    real(RP) :: c_in_map(1)    ! c_in  horizontal distribution ! [Add] 11/08/30 T.Mitsui
    real(RP) :: esw(KA)       ! saturation vapor pressure, water
    real(RP) :: esi(KA)       !                            ice
    real(RP) :: ssw(KA)       ! super saturation (water)
    real(RP) :: ssi(KA)       ! super saturation (ice)
!    real(RP) :: w_dsswdz(KA)  ! w*(d_ssw/ d_z) super saturation(water) flux
    real(RP) :: w_dssidz(KA)  ! w*(d_ssi/ d_z), 09/04/14 T.Mitsui
!    real(RP) :: ssw_below(KA) ! ssw(k-1)
    real(RP) :: ssi_below(KA) ! ssi(k-1), 09/04/14 T.Mitsui
    real(RP) :: z_below(KA)   ! z(k-1)
    real(RP) :: dzh           ! z(k)-z(k-1)
    real(RP) :: pv            ! vapor pressure
    ! work variables for Twomey Equation.
    real(RP) :: qsw(KA)
    real(RP) :: qsi(KA)
    real(RP) :: dqsidtem_rho(KA)
    real(RP) :: dssidt_rad(KA)
    real(RP) :: wssi, wdssi
    !
 
    real(RP) :: cva(KA)
    real(RP) :: dq_xa(KA,HYDRO_MAX)
    real(RP) :: vt_xa(KA,HYDRO_MAX,2) ! terminal velocity
    real(RP) :: dTdt_dep(KA)
    real(RP) :: PLIdep_total(KA)
    real(RP) :: wtem(KA)         ! temperature[K]
    real(RP) :: dqsidpre_tem(KA)
    real(RP) :: dqsidtem_pre(KA)
    real(RP) :: dssidt
    real(RP) :: dssidt_mp(KA)
    real(RP) :: dssidt_dyn(KA)
    !!
 
 
!    real(RP) :: xi_nuc(1)    ! xi use the value @ cloud base
!    real(RP) :: alpha_nuc(1) ! alpha_nuc
!    real(RP) :: eta_nuc(1)   ! xi use the value @ cloud base
    !
    real(RP) :: sigma_w(KA)
    real(RP) :: weff(KA)
    real(RP) :: weff_max(KA)
    real(RP) :: velz(KA)
    !
    real(RP) :: coef_ccn
    real(RP) :: slope_ccn
    real(RP) :: nc_new(KA)
    real(RP) :: nc_new_below(KA)
    real(RP) :: dnc_new
    real(RP) :: nc_new_max   ! Lohmann (2002),
    real(RP) :: a_max
    real(RP) :: b_max
    logical :: flag_nucleation(KA)
    !
    real(RP) :: r_gravity
    real(RP), parameter :: r_sqrt3=0.577350269_rp ! = sqrt(1.0/3.0)
    real(RP), parameter :: eps=1.e-30_rp
    !====> ! 09/08/18
    !
    real(RP) :: dlcdt_max, dli_max ! defined by supersaturation
    real(RP) :: dncdt_max, dni_max ! defined by supersaturation
    real(RP) :: rdt
 
    real(RP) :: tmp
    !
    integer :: k
    !
    !
 
    if( so22_het ) then
       do k = ks, ke
          ! Temperature lower limit is only used for saturation condition.
          ! On the other hand original "tem" is used for calculation of latent heat or energy equation.
          wtem(k)  = max( tem(k), tem_min )
       end do
    endif
 
 
!    c_ccn_map(1) = c_ccn
!    kappa_map(1) = kappa
!    c_in_map(1)  = c_in
    !
    !
    rdt            = 1.0_rp/dt
    r_gravity      = 1.0_rp/grav
    !
    call moist_psat_liq      ( ka, ks, ke, &
                              tem(:), esw(:) ) ! [IN], [OUT]
    call moist_psat_ice      ( ka, ks, ke, &
                              tem(:), esi(:) ) ! [IN], [OUT]
    call moist_pres2qsat_liq ( ka, ks, ke, &
                              tem(:), pre(:), qdry(:), & ! [IN]
                              qsw(:)                   ) ! [OUT]
    call moist_pres2qsat_ice ( ka, ks, ke, &
                              tem(:), pre(:), qdry(:), & ! [IN]
                              qsi(:)                   ) ! [OUT]
    call moist_dqsi_dtem_dens( ka, ks, ke, &
                               tem(:), rho(:), & ! [IN]
                               dqsidtem_rho(:) ) ! [OUT]
 
    if( so22_het ) then
       call moist_dqs_dtem_dpre_ice( ka, ks, ke, &
                                     wtem(:), pre(:), qdry(:),        & ! [IN]
                                     dqsidtem_pre(:), dqsidpre_tem(:) ) ! [OUT]
    endif
    !!
 
 
    !
    ! Lohmann (2002),JAS, eq.(1) but changing unit [cm-3] => [m-3]
    a_max = 1.e+6_rp*0.1_rp*(1.e-6_rp)**1.27_rp
    b_max = 1.27_rp
    !
    ssi_max = 1.0_rp
 
    do k = ks, ke
       pv     = rhoq(k,i_qv)*rvap*tem(k)
       ssw(k) = min( mp_ssw_lim, ( pv/esw(k)-1.0_rp ) )*100.0_rp
       ssi(k) = ( pv/esi(k) - 1.00_rp )
!       ssw_below(k+1) = ssw(k)
       ssi_below(k+1) = ssi(k)
       z_below(k+1)   = cz(k)
    end do
!    ssw_below(KS) = ssw(KS)
    ssi_below(ks) = ssi(ks)
    z_below(ks)   = cz(ks-1)
 
    ! dS/dz is evaluated by first order upstream difference
    !***  Solution for Twomey Equation ***
!    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)) * &
    coef_ccn  = 1.e+6_rp*0.88_rp*(c_ccn*1.e-6_rp)**(2.0_rp/(kappa + 2.0_rp)) &
!           * (70.0_RP)**(kappa_map(1)/(kappa_map(1) + 2.0_RP))
            * (70.0_rp)**(kappa/(kappa + 2.0_rp))
!    slope_ccn = 1.5_RP*kappa_map(1)/(kappa_map(1) + 2.0_RP)
    slope_ccn = 1.5_rp*kappa/(kappa + 2.0_rp)
    !
    do k=ks, ke
       sigma_w(k) = r_sqrt3*sqrt(max(qke(k),qke_min))
    end do
    sigma_w(ks-1) = sigma_w(ks)
    sigma_w(ke+1) = sigma_w(ke)
    ! effective vertical velocity
    do k=ks, ke
       weff(k) = w(k) - cpa(k)*r_gravity*dtdt_rad(k)
    end do
    !
    if( mp_couple_aerosol ) then
 
       do k = ks, ke
          if( ssw(k) > 1.e-10_rp .AND. pre(k) > 300.e+2_rp ) then
             nc_new(k) = max( ccn(k), c_ccn )
          else
             nc_new(k) = 0.0_rp
          endif
       enddo
 
    else
 
       if( nucl_twomey ) then
          ! diagnose cloud condensation nuclei
          do k = ks, ke
             ! effective vertical velocity (maximum vertical velocity in turbulent flow)
             weff_max(k) = weff(k) + sigma_w(k)
             ! large scale upward motion region and saturated
             if( (weff(k) > 1.e-8_rp) .AND. (ssw(k) > 1.e-10_rp)  .AND. pre(k) > 300.e+2_rp )then
                ! Lohmann (2002), eq.(1)
                nc_new_max   = coef_ccn*weff_max(k)**slope_ccn
                nc_new(k) = a_max*nc_new_max**b_max
             else
                nc_new(k) = 0.0_rp
             end if
          end do
 
       else
          ! calculate cloud condensation nuclei
          do k = ks, ke
             if( ssw(k) > 1.e-10_rp .AND. pre(k) > 300.e+2_rp ) then
                nc_new(k) = c_ccn*ssw(k)**kappa
             else
                nc_new(k) = 0.0_rp
             endif
          enddo
      endif
 
    endif
 
    do k = ks, ke
       ! nc_new is bound by upper limit
       if( nc_new(k) > nc_uplim_d ) then ! no more CCN
          flag_nucleation(k) = .false.
          nc_new_below(k+1)  = 1.e+30_rp
       else if( nc_new(k) > eps ) then ! nucleation can occur
          flag_nucleation(k) = .true.
          nc_new_below(k+1)  = nc_new(k)
       else ! nucleation cannot occur(unsaturated or negative w)
          flag_nucleation(k) = .false.
          nc_new_below(k+1)  = 0.0_rp
       end if
    end do
    nc_new_below(ks) = 0.0_rp
!    do k=KS, KE
        ! search maximum value of nc_new
!       if(  ( nc_new(k) < nc_new_below(k) ) .OR. &
!            ( nc_new_below(k) > c_ccn_map(1)*0.05_RP ) )then ! 5% of c_ccn
!            ( nc_new_below(k) > c_ccn*0.05_RP ) )then ! 5% of c_ccn
!          flag_nucleation(k) = .false.
!       end if
!    end do
 
    if( mp_couple_aerosol ) then
       do k = ks, ke
          ! nucleation occurs at only cloud base.
          ! if CCN is more than below parcel, nucleation newly occurs
          ! effective vertical velocity
          if ( flag_nucleation(k) .AND. & ! large scale upward motion region and saturated
               tem(k) > tem_ccn_low ) then
             dlcdt_max     = ( rhoq(k,i_qv) - esw(k) / ( rvap * tem(k) ) ) * rdt
             dlcdt_max     = max( dlcdt_max, 0.0_rp ) ! dlcdt_max can be artificially negative due to truncation error in floating point operation
             dncdt_max     = dlcdt_max/xc_min
!             dnc_new       = nc_new(k)-rhoq(k,I_NC)
             dnc_new       = nc_new(k)
             pq(k,i_ncccn) = min( dncdt_max, dnc_new*rdt )
             pq(k,i_lcccn) = min( dlcdt_max, xc_min*pq(k,i_ncccn) )
          else
             pq(k,i_ncccn) = 0.0_rp
             pq(k,i_lcccn) = 0.0_rp
          end if
       end do
    else
 
       if( nucl_twomey ) then
          do k = ks, ke
             ! nucleation occurs at only cloud base.
             ! if CCN is more than below parcel, nucleation newly occurs
             ! effective vertical velocity
             if ( flag_nucleation(k) .AND. & ! large scale upward motion region and saturated
                  tem(k) > tem_ccn_low .AND. &
                  nc_new(k) > rhoq(k,i_nc) ) then
                dlcdt_max     = ( rhoq(k,i_qv) - esw(k) / ( rvap * tem(k) ) ) * rdt
                dlcdt_max     = max( dlcdt_max, 0.0_rp ) ! dlcdt_max can be artificially negative due to truncation error in floating point operation
                dncdt_max     = dlcdt_max/xc_min
                dnc_new       = nc_new(k)-rhoq(k,i_nc)
                pq(k,i_ncccn) = min( dncdt_max, dnc_new*rdt )
                pq(k,i_lcccn) = min( dlcdt_max, xc_min*pq(k,i_ncccn) )
             else
                pq(k,i_ncccn) = 0.0_rp
                pq(k,i_lcccn) = 0.0_rp
             end if
          end do
       else
          do k = ks, ke
             ! effective vertical velocity
             if(  tem(k) > tem_ccn_low .AND. &
                  nc_new(k) > rhoq(k,i_nc) ) then
                dlcdt_max     = ( rhoq(k,i_qv) - esw(k) / ( rvap * tem(k) ) ) * rdt
                dlcdt_max     = max( dlcdt_max, 0.0_rp ) ! dlcdt_max can be artificially negative due to truncation error in floating point operation
                dncdt_max     = dlcdt_max/xc_min
                dnc_new       = nc_new(k)-rhoq(k,i_nc)
                pq(k,i_ncccn) = min( dncdt_max, dnc_new*rdt )
                pq(k,i_lcccn) = min( dlcdt_max, xc_min*pq(k,i_ncccn) )
             else
                pq(k,i_ncccn) = 0.0_rp
                pq(k,i_lcccn) = 0.0_rp
             end if
          end do
       endif
 
    endif
 
    !
    ! ice nucleation
    !
    ! +++ NOTE ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    ! Based on Phillips etal.(2006).
    ! However this approach doesn't diagnose Ni itself but diagnose tendency.
    ! Original approach adjust Ni instantaneously .
    ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
    if( so22_het .or. opt_nucleation_ice_hom ) then
       call dep_vapor_ice_wrk(         & ! in
            ka, ks, ke,                & ! in
            plidep_total(:),           & ! out for ice nucleation
            rho(:), tem(:), pre(:),    & ! in
            qdry(:), esi(:), qsi(:),   & ! in
            rhoq(:,:),                 & ! in
            vt_xa, dq_xa,              & ! in
            dt                         ) ! in
 
       do k = ks, ke
          dtdt_dep(k) = (lhs0+(cv_vapor-cv_ice)*tem(k))*plidep_total(k)/(rho(k)*cva(k))
       enddo
    else
       dtdt_dep(:) = 0.0_rp
    endif
 
    do k = ks, ke-1
       velz(k)           = ( w(k) * ( cz(k+1) - fz(k) ) + w(k+1) * ( fz(k) - cz(k) ) ) / ( cz(k+1) - cz(k) ) ! @ half level
    end do
    velz(ke) = 0.0_rp
    do k = ks, ke
       dzh           = cz(k) - z_below(k)
       w_dssidz(k)   = velz(k) * (ssi(k) - ssi_below(k))/dzh
       dssidt_rad(k) = -rhoq(k,i_qv)/(rho(k)*qsi(k)*qsi(k))*dqsidtem_rho(k)*dtdt_rad(k)
       dli_max       = ( rhoq(k,i_qv) - esi(k) / ( rvap * tem(k) ) ) * rdt
       dli_max       = max( dli_max, 0.0_rp ) ! dli_max can be artificially negative due to truncation error in floating point operation
       dni_max       = min( dli_max/xi_ccn, (in_max-rhoq(k,i_ni))*rdt )
       wdssi         = min( w_dssidz(k)+dssidt_rad(k), 0.01_rp)
       wssi          = min( ssi(k), ssi_max)
 
       !! Seiki and Ohno 2022
       if( so22_het ) then
          dssidt_mp(k)  = -plidep_total(k)/(rho(k)*qsi(k))
          !!  PLIdep(k) -> PQ(k,I_LIdep)
          dssidt_rad(k) = -rhoq(k,i_qv)/(rho(k)*qsi(k)*qsi(k))*dqsidtem_rho(k)*dtdt_rad(k)
          !!
          dssidt_dyn(k) = +rhoq(k,i_qv)/(rho(k)*qsi(k)*qsi(k))&
               * velz(k)*grav*(dqsidtem_pre(k)/cpa(k)+dqsidpre_tem(k)*rho(k))
          !    * w(k)*CNST_GRAV*(dqsidtem_pre(ij,k)/cpa(ij,k)+dqsidpre_tem(ij,k)*rho(ij,k))
          dssidt        = dssidt_mp(k) + dssidt_rad(k) + dssidt_dyn(k)
       endif
 
       ! SB06(34),(35)
!#       if(  ( wdssi       > eps         ) .AND. & !
!#            (tem(k)       < 273.15_RP   ) .AND. & !
!#            (rhoq(k,I_NI) < in_max      ) .AND. &
!#            (wssi      >= eps       ) )then   !
!#          tmp = c_in * nm_M92 * exp( 0.3_RP * bm_M92 * ( wssi - 0.1_RP ) )
!#          if( inucl_w ) then
!#             tmp = bm_M92 * 0.3_RP * tmp * wdssi
!#          else
!#             tmp = max( tmp - rhoq(k,I_NI), 0.0_RP ) * rdt
!#          endif
!#          PQ(k,I_NIccn) = min(dni_max, tmp)
!#          PQ(k,I_LIccn) = min(dli_max, PQ(k,I_NIccn)*xi_ccn )
!#       else
!#          PQ(k,I_NIccn) = 0.0_RP
!#          PQ(k,I_LIccn) = 0.0_RP
!#       end if
       if(  (tem(k)     < 273.15_rp   ) .AND. & !
            (rhoq(k,i_ni) < in_max      ) .AND. &
            (wssi        >= eps       ) )then   !
          tmp = c_in * nm_m92 * exp( 0.3_rp * bm_m92 * ( wssi - 0.1_rp ) )
          if( inucl_w .and. wdssi > eps ) then
             tmp = bm_m92 * 0.3_rp * tmp * wdssi
          elseif( so22_het .and. dssidt > eps ) then
             tmp = bm_m92 * 0.3_rp * tmp * dssidt
          else
             tmp = max( tmp - rhoq(k,i_ni), 0.0_rp ) * rdt
          endif
          pq(k,i_niccn) = min(dni_max, tmp)
          pq(k,i_liccn) = min(dli_max, pq(k,i_niccn)*xi_ccn )
       else
          pq(k,i_niccn) = 0.0_rp
          pq(k,i_liccn) = 0.0_rp
       end if
 
 
    end do
 
    if( opt_nucleation_ice_hom ) then
       call nucleation_ice_hom(   &
            ka, ks, ke,           & !in
            tem, pre, rho,        & !in
            qdry, rhoq(:,i_qv),   & !in
            cva, cpa,             & !in
            w,                    & !in
            dtdt_rad,             & !in
            dtdt_dep,             & !in
            plidep_total, dt,     & !in
            pq(:,i_lihom),        & !out
            pq(:,i_nihom)         ) !out
    else
       pq(:,i_lihom) = 0.0_rp
       pq(:,i_nihom) = 0.0_rp
    endif
 
 
    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 3939 of file scale_atmos_phy_mp_sn14.F90.

    use scale_prc, only: &
         prc_abort
    use scale_const, only: &
         psat0  => const_psat0,  &
         t00    => const_tem00,  &
         pstd   => const_pstd,   &
         cpvap  => const_cpvap,  &
         cvvap  => const_cvvap,  &
         ci     => const_ci,     &
         cl     => const_cl,     &
         rvap   => const_rvap,   &
         rdry   => const_rdry,   &
         lhs00  => const_lhs00,  &
         lhs0   => const_lhs0,   &
         lhv00  => const_lhv00,  &
         lhf00  => const_lhf00,  &
         epsvap => const_epsvap, &
         grav   => const_grav,   &
         pi     => const_pi
    implicit none
 
    integer, intent(in)  :: KA, KS, KE
 
    real(RP), intent(in) :: tem(KA)
    real(RP), intent(in) :: pre(KA)
    real(RP), intent(in) :: rho(KA)
    real(RP), intent(in) :: qd(KA)
    real(RP), intent(in) :: rhoq_qv(KA)
    real(RP), intent(in) :: cpa(KA)        ! specific heat @ cnst. pressure
    real(RP), intent(in) :: cva(KA)        ! specific heat @ cnst. volume
    ! balance equation @ homogeneous ice nucleation
    real(RP), intent(in) :: w(KA)          ! adiabatic ascending
    real(RP), intent(in) :: dTdt_rad(KA)  ! effect of radiative heating
    real(RP), intent(in) :: dTdt_dep(KA)  ! effect of radiative heating
    real(RP), intent(in) :: PLIdep(KA)
!    real(RP), intent(out) :: PQ(KA,PQ_MAX)     ! competition effect by vapor consumption
    real(RP), intent(out):: PLIhom(KA)
    real(RP), intent(out):: PNIhom(KA)
    real(RP), intent(in) :: dt
    real(RP), parameter :: rhoi=916.0_rp           ! ice density [kg/m3]
    real(RP), parameter :: rrhoi=1.0_rp/rhoi       !
    real(RP), parameter :: Mw=18.01528_rp         ! Water Molar Weight
    real(RP), parameter :: Nav=6.0221415e+23_rp      ! Avogadro Number
    real(RP), parameter :: vw=(mw*1.e-3_rp/nav)/rhoi ! volume of a water molecule in ice phase
    !                     18nm*exp( 3.0*log(1.5)**2 )=29.4760367
    real(RP), parameter :: r0=29.5e-9_rp          ! aerosol mass(volume) mode radius
    real(RP), parameter :: c_gf =  1.01187_rp     ! dAlmeida
    real(RP), parameter :: g_gf = -0.206449_rp    ! dAlmeida
    real(RP), parameter :: rho_min=1.e-5_rp         ! 3.e-3 is lower limit recognized in many experiments.
    real(RP), parameter :: tem_min=150.0_rp
    real(RP), parameter :: ni_max =300.e+6_rp
!    real(RP) :: dTdt_dep(KA)  ! effect of deposition heating
    real(RP) :: wtem
    real(RP) :: esi, esw
    real(RP) :: qsi
    real(RP) :: lhs
    real(RP) :: dqsidtem
    real(RP) :: den1, den2
    real(RP) :: dqsidt_pre
    real(RP) :: dqsidp_tem
    real(RP) :: rw
    real(RP) :: temc_lim
    real(RP) :: rho_lim
    real(RP) :: pre_lim
    real(RP) :: Dw
    real(RP) :: si, sw
    real(RP) :: Scr
    real(RP) :: dsidt_mp
    real(RP) :: dsidt_rd
    real(RP) :: wp
    real(RP) :: a1,a2,a3
    real(RP) :: b1,b2
    real(RP) :: dlogJdT
    real(RP) :: dtemdt_dyn
    real(RP) :: rtau ! 1/tau
    real(RP) :: delta
    real(RP) :: kappa
    real(RP) :: Rim_w
    real(RP) :: ri_wrk
    real(RP) :: ri
    !
    real(RP), parameter :: r2pi    = 0.5_rp/pi ! 1/2pi
    real(RP), parameter :: sqrt_pi = sqrt(pi)
    real(RP), parameter :: coef_mi = 4.0_rp/3.0_rp*pi*rhoi
    real(RP), parameter :: eps     = 1.e-30_rp
    real(RP) :: rdt
    integer :: ierr
    integer :: ij,k
    !
 
    rdt     = 1.0_rp/dt
    plihom(:)     = 0.0_rp
    pnihom(:)     = 0.0_rp
!    PQ(:,I_NIhom)     = 0.0_RP
!    PQ(:,I_LIhom)     = 0.0_RP
!    PQ(1:KS,I_NIhom)     = 0.0_RP
!    PQ(1:KS,I_LIhom)     = 0.0_RP
!    PQ(KE:KA,I_NIhom)     = 0.0_RP
!    PQ(KE:KA,I_LIhom)     = 0.0_RP
    do k = ks, ke
          wtem= max(tem(k), tem_min)
          esi = min( psat0 &
               * ( wtem / t00 ) ** ( ( cpvap - ci ) / rvap ) &
               * exp( lhs00 / rvap &
               * ( 1.0_rp / t00 - 1.0_rp / wtem ) ), pre(k))
          esw = psat0 &
               * ( wtem / t00 ) ** ( ( cpvap - cl ) / rvap ) &
               * exp( lhv00 / rvap &
               * ( 1.0_rp / t00 - 1.0_rp / wtem ) )
          ! (10) in Ren and MacKenzie (2005)
          scr       = 2.349_rp - wtem/259.0_rp
          qsi       = epsvap * esi / ( pre(k) - ( 1.0_rp - epsvap ) * esi )
          si        = rhoq_qv(k)*rvap*wtem/esi ! rho(k)*qv(k)
          ! sw must be less than 1 to calculate hygroscopic growth (otherwise, already activated)
          sw        = min(rho(k)*rhoq_qv(k)*rvap*wtem/esw,0.999_rp)
          if ( si < scr ) then
             ! No nucleation occurs
             qsi       = 0.0_rp
             lhs       = lhs00 + (cpvap - ci )*(wtem-t00)
             dqsidtem  = 0.0_rp
             si        = 0.0_rp
             sw        = 0.0_rp
             scr       = 100.0_rp
             a1        = 0.0_rp
             wp        = -1.0_rp
             rtau      = -1.0_rp
          else
             ! (dsi/dt)_MP
             lhs       = lhs0 + (cpvap - ci)*(wtem-t00)
             dqsidtem  = esi/(rho(k)*rvap*wtem*wtem)&
                  * (lhs/(rvap*wtem)-1.0_rp)
             dsidt_mp  = -1.0_rp/(rho(k)*qsi) &
                  * (1.0_rp+si*(lhv00+lhf00+(cvvap-ci)*wtem)/cva(k)*dqsidtem)&
                  * plidep(k)
             dsidt_rd  = -si/qsi*dqsidtem*dtdt_rad(k)
             ! This is an simple formulation from original paper
             a1        = lhs0*grav/(cpa(k)*rvap*tem(k)*tem(k))&
                  - grav/(rdry*tem(k))
!!$          ! alternative formulation for NICAM-EXACT
!!$          den1      = (pre(k)-(1.0_RP-EPSvap)*esi)*(pre(k)-(1.0_RP-EPSvap)*esi)
!!$          den2      = den1*Rvap*wtem*wtem
!!$          dqsidp_tem=-EPSvap*              esi/den1
!!$          dqsidt_pre= EPSvap*pre(k)*lhs*esi/den2
!!$          a1        = GRAV/qsi*(dqsidt_pre/cpa(k)+dqsidp_tem*rho(k))
             wp        = w(k) + 1.0_rp/(a1*si)*(dsidt_mp+dsidt_rd)
             ! (21) in RM05
             dlogjdt   = -0.004_rp*wtem*wtem + 2.0_rp*wtem - 304.4_rp
             ! (22) in RM05
             dtemdt_dyn= - grav*w(k)/cpa(k)
             rtau      =  dlogjdt*(dtdt_rad(k)+dtdt_dep(k)+dtemdt_dyn )
          endif
          rw = r0*c_gf*(1-sw)**g_gf
 
          if ( wp > eps .AND. rtau > eps ) then
             ! a1,a2,a3,b1,b2 are given by Appendix B in RM05
             temc_lim= max(tem(k)-t00, temc_lim_diff )
             rho_lim = max(rho(k),rho_min)             !
             pre_lim = rho_lim*(qd(k)*rdry + rhoq_qv(k)*rvap)*(temc_lim+t00)! only for Dw
             dw      = 0.211e-4_rp* (((temc_lim+t00)/t00)**1.94_rp) *(pstd/pre_lim)
             ! v_th = sqrt(8*Rv*T/pi) in the paragraph after eq(2) in Karcher etal.(2006)
             b2      = 0.5_rp/dw*sqrt(rvap*wtem*r2pi)
             b1      = (si-1.0_rp)*0.5_rp*rrhoi*esi/sqrt(2.0_rp*pi*rvap*wtem)
             a2      = mw*rvap*wtem/(nav*esi)
             a3      = epsvap*mw*lhs0*lhs0/(nav*cpa(k)*pre(k)*wtem)
             delta   = b2*rw
             kappa   = 2.0_rp*b1*b2/( rtau*(1.0_rp+delta)*(1.0_rp+delta) )
             ! (A9) in RM05
             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))) &
                  +                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 )
             pnihom(k)  = min( scr/(scr-1.0_rp)*a1*wp/(rim_w*4.0_rp*pi*dw/b2), ni_max )* rdt
             ri_wrk        = 1.0_rp+b2*rw
             ri            = ( sqrt(ri_wrk*ri_wrk + 2.0_rp*b1*b2*dt )-1.0_rp )/b2
             plihom(k)  = coef_mi*ri*ri*ri*pnihom(k)
          else
             plihom(k)  = 0.0_rp
             pnihom(k)  = 0.0_rp
          endif
    enddo
    !
    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 4128 of file scale_atmos_phy_mp_sn14.F90.

 
    ! ice multiplication by splintering
    ! we consider Hallet-Mossop process
    use scale_specfunc, only: &
       gammafunc => sf_gamma
    implicit none
    integer, intent(in) :: KA, KS, KE
 
    real(RP), intent(in) :: Pac(KA,Pac_MAX)
    real(RP), intent(in) :: tem(KA)
    real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
    real(RP), intent(in) :: xq(KA,HYDRO_MAX)
    !
    real(RP), intent(inout):: PQ(KA,PQ_MAX)
    ! for lightning
    logical,  intent(in) :: flg_lt
    real(RP), intent(in) :: rhoq_crg(KA,I_QC:I_QG)
    real(RP), intent(inout):: Pcrg1(KA,PQ_MAX)
    !
    ! production of (350.d3 ice particles)/(cloud riming [g]) => 350*d6 [/kg]
    real(RP), parameter :: pice = 350.0e+6_rp
    ! production of (1 ice particle)/(250 cloud particles riming)
    real(RP), parameter :: pnc  = 250.0_rp
    ! temperature factor
    real(RP) :: fp
 
    real(RP) :: igm            ! in complete gamma(x,alpha)
    real(RP) :: x
    ! coefficient of expansion using in calculation of igm
    real(RP) :: a0,a1,a2,a3,a4,a5
    real(RP) :: a6,a7,a8,a9,a10
    real(RP) :: an1,an2,b0,b1,b2,c0,c1,c2
    real(RP) :: d0,d1,d2,e1,e2,h0,h1,h2
    real(RP), parameter :: eps=1.0e-30_rp
    ! number of cloud droplets larger than 12 micron(radius).
    real(RP) :: n12
    !
    real(RP) :: wn, wni, wns, wng
    integer :: k, iq
    real(RP) :: sw1
    !
    !
!OCL LOOP_FISSION_TARGET(LS)
    do k = ks, ke
       ! Here we assume particle temperature is same as environment temperature.
       ! If you want to treat in a better manner,
       ! you can diagnose with eq.(64) in CT(86)
       if     (tem(k) >  270.16_rp)then
          fp   = 0.0_rp
       else if(tem(k) >= 268.16_rp)then
          fp   = (270.16_rp-tem(k))*0.5_rp
       else if(tem(k) >= 265.16_rp)then
          fp   = (tem(k)-265.16_rp)*0.333333333_rp
       else
          fp   = 0.0_rp
       end if
       ! Approximation of Incomplete Gamma function
       ! Here we calculate with algorithm by Numerical Recipes.
       ! This approach is based on eq.(78) in Cotton etal.(1986),
       ! but more accurate and expanded for Generalized Gamma Distribution.
       x       = coef_lambda(i_mp_qc)*(xc_cr/xq(k,i_mp_qc))**mu(i_mp_qc)
       !
       if(x<1.e-2_rp*alpha)then        ! negligible
          igm  = 0.0_rp
       else if(x<alpha+1.0_rp)then    ! series expansion
          ! 10th-truncation is enough for cloud droplet.
          a0   = 1.0_rp/alpha         ! n=0
          a1   = a0*x/(alpha+1.0_rp)  ! n=1
          a2   = a1*x/(alpha+2.0_rp)  ! n=2
          a3   = a2*x/(alpha+3.0_rp)  ! n=3
          a4   = a3*x/(alpha+4.0_rp)  ! n=4
          a5   = a4*x/(alpha+5.0_rp)  ! n=5
          a6   = a5*x/(alpha+6.0_rp)  ! n=6
          a7   = a6*x/(alpha+7.0_rp)  ! n=7
          a8   = a7*x/(alpha+8.0_rp)  ! n=8
          a9   = a8*x/(alpha+9.0_rp)  ! n=9
          a10  = a9*x/(alpha+10.0_rp) ! n=10
          igm  = (a0+a1+a2+a3+a4+a5+a6+a7+a8+a9+a10)*exp( -x + alpha*log(x) - lgm )
       else if(x<alpha*100.0_rp) then ! continued fraction expansion
          ! 2nd-truncation is enough for cloud droplet.
          ! setup
          b0   = x+1.0_rp-alpha
          c0   = 1.0_rp/eps
          d0   = 1.0_rp/b0
          h0   = d0
          ! n=1
          an1  = -(1.0_rp-alpha)
          b1   = b0 + 2.0_rp
          d1   = 1.0_rp/(an1*d0+b1)
          c1   = b1+an1/c0
          e1   = d1*c1
          h1   = h0*e1
          ! n=2
          an2  = -2.0_rp*(2.0_rp-alpha)
          b2   = b1 + 2.0_rp
          d2   = 1.0_rp/(an2*d1+b2)
          c2   = b2+an2/c1
          e2   = d2*c2
          h2   = h1*e2
          !
          igm  = 1.0_rp - exp( -x + alpha*log(x) - lgm )*h2
       else                       ! negligible
          igm  = 1.0_rp
       end if
       ! n12 is number of cloud droplets larger than 12 micron.
       n12 = rhoq(k,i_nc)*(1.0_rp-igm)
       ! eq.(82) CT(86)
       wn            = (pice + n12/((rhoq(k,i_qc)+xc_min)*pnc) )*fp ! filtered by xc_min
       wni           = wn*(-pac(k,i_liaclc2li)  ) ! riming production rate is all negative
       wns           = wn*(-pac(k,i_lsaclc2ls)  )
       wng           = wn*(-pac(k,i_lgaclc2lg)  )
       pq(k,i_nispl) = wni+wns+wng
       !
       pq(k,i_lsspl) = - wns*xq(k,i_mp_qi) ! snow    => ice
       pq(k,i_lgspl) = - wng*xq(k,i_mp_qi) ! graupel => ice
       if (flg_lt) then
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NC is small,ignore charge transfer
          pcrg1(k,i_nsspl) = - wns*(1.0_rp-sw1) &
                                / (rhoq(k,i_ns)+sw1)*rhoq_crg(k,i_qs)
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small ) !--- if NC is small,ignore charge transfer
          pcrg1(k,i_ngspl) = - wng*(1.0_rp-sw1) &
                                / (rhoq(k,i_ng)+sw1)*rhoq_crg(k,i_qg)
          pcrg1(k,i_nispl) = - ( pcrg1(k,i_nsspl) + pcrg1(k,i_ngspl) )
       end if
       !
    end do
    !
    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 4267 of file scale_atmos_phy_mp_sn14.F90.

    use scale_atmos_saturation, only: &
       moist_psat_ice => atmos_saturation_psat_ice
    implicit none
 
    integer, intent(in) :: KA, KS, KE
 
    !--- mixed-phase collection process
    !                  And all we set all production term as a negative sign to avoid confusion.
    !
    real(RP), intent(in) :: wtem(KA)
    !--- mass/number concentration[kg/m3]
    real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
    ! necessary ?
    real(RP), intent(in) :: xq(KA,HYDRO_MAX)
    !--- diameter of averaged mass( D(ave x) )
    real(RP), intent(in) :: dq_xave(KA,HYDRO_MAX)
    !--- terminal velocity of averaged mass( vt(ave x) )
    real(RP), intent(in) :: vt_xave(KA,HYDRO_MAX,2)
    ! [Add] 11/08/30 T.Mitsui, for autoconversion of ice
    !    real(RP), intent(in) :: rho(KA)
    !--- partial conversion
    real(RP), intent(inout):: PQ(KA,PQ_MAX)
    !
    real(RP), intent(out):: Pac(KA,Pac_MAX)
    !--- for lightning component
    logical,  intent(in) :: flg_lt
    real(RP), intent(in) :: beta_crg(KA)
    real(RP), intent(in) :: dqcrg(KA)
    real(RP), intent(in) :: d0_crg, v0_crg
    real(RP), intent(in) :: rhoq_crg(KA,I_QC:I_QG)
    real(RP), intent(inout):: Pcrg1(KA,PQ_MAX)
    real(RP), intent(inout):: Pcrg2(KA,Pcrg_MAX)
 
    real(RP), parameter :: a_dec = 0.883_rp
    real(RP), parameter :: b_dec = 0.093_rp
    real(RP), parameter :: c_dec = 0.00348_rp
    real(RP), parameter :: d_dec = 4.5185e-5_rp
    !
    !
    real(RP) :: tem(KA)
    !
    !--- collection efficency of each specie
    real(RP) :: E_c(KA), E_r, E_i, E_s, E_g
    real(RP) :: E_ic, E_sc, E_gc
    !--- sticking efficiency
    real(RP) :: E_stick(KA)
    ! [Add] 10/08/03 T.Mitsui
    real(RP) :: temc, temc2, temc3
    real(RP) :: E_dec
    real(RP) :: esi_rat
    real(RP) :: esi(KA)
    !
    real(RP) :: temc_p, temc_m             ! celcius tem.
!    real(RP) :: ci_aut(KA)
!    real(RP) :: taui_aut(KA)
!    real(RP) :: tau_sce(KA)
    !--- DSD averaged diameter for each species
    real(RP) :: ave_dc                     ! cloud
!    real(RP) :: ave_dr                     ! rain
    real(RP) :: ave_di(KA)                 ! ice
    real(RP) :: ave_ds(KA)                 ! snow
    real(RP) :: ave_dg                     ! graupel
    !--- coefficient of collection equations(L:mass, N:number)
    real(RP) :: coef_acc_LCI,   coef_acc_NCI   ! cloud     - cloud ice
    real(RP) :: coef_acc_LCS,   coef_acc_NCS   ! cloud     - snow
    !
    real(RP) :: coef_acc_LCG,   coef_acc_NCG   ! cloud     - graupel
    real(RP) :: coef_acc_LRI_I, coef_acc_NRI_I ! rain      - cloud ice
    real(RP) :: coef_acc_LRI_R, coef_acc_NRI_R ! rain      - cloud ice
    real(RP) :: coef_acc_LRS_S, coef_acc_NRS_S ! rain      - snow
    real(RP) :: coef_acc_LRS_R, coef_acc_NRS_R ! rain      - snow
    real(RP) :: coef_acc_LRG,   coef_acc_NRG   ! rain      - graupel
    real(RP) :: coef_acc_LII,   coef_acc_NII   ! cloud ice - cloud ice
    real(RP) :: coef_acc_LIS,   coef_acc_NIS   ! cloud ice - snow
    real(RP) ::                 coef_acc_NSS   ! snow      - snow
    real(RP) ::                 coef_acc_NGG   ! grauepl   - graupel
    real(RP) :: coef_acc_LSG,   coef_acc_NSG(KA) ! snow      - graupel
    !--- (diameter) x (diameter)
    real(RP) :: dcdc(KA), dcdi, dcds, dcdg
    real(RP) :: drdr(KA), drdi(KA), drds(KA), drdg
    real(RP) :: didi(KA), dids, didg
    real(RP) :: dsds(KA), dsdg
    real(RP) :: dgdg(KA)
    !--- (terminal velocity) x (terminal velocity)
    real(RP) :: vcvc(KA), vcvi, vcvs, vcvg
    real(RP) :: vrvr(KA), vrvi(KA), vrvs(KA), vrvg
    real(RP) :: vivi(KA), vivs, vivg
    real(RP) :: vsvs(KA), vsvg
    real(RP) :: vgvg(KA)
    !
    real(RP) :: wx_cri, wx_crs
    real(RP) :: coef_emelt
    real(RP) :: w1
 
    real(RP) :: sw, sw1, sw2
    real(RP) :: alpha_lt
    !
    integer :: k, iqw
    !
    !
    do k = ks, ke
       tem(k) = max( wtem(k), tem_min ) ! 11/08/30 T.Mitsui
    end do
 
    call moist_psat_ice( ka, ks, ke, &
                         tem(:), esi(:)  ) ! [IN], [OUT]
 
    if( opt_stick_ks96 )then
       do k = ks, ke
          ! Khain and Sednev (1996), eq.(3.15)
          temc       = tem(k) - t00
          temc2      = temc*temc
          temc3      = temc2*temc
          e_dec      = max(0.0_rp, a_dec + b_dec*temc + c_dec*temc2 + d_dec*temc3 )
          esi_rat    = rhoq(k,i_qv)*rvap*tem(k)/esi(k)
          e_stick(k) = min(1.0_rp, e_dec*esi_rat)
       end do
    else if( opt_stick_co86 )then
       do k = ks, ke
          ! [Add] 11/08/30 T.Mitsui, Cotton et al. (1986)
          temc       = min(tem(k) - t00,0.0_rp)
          w1         = 0.035_rp*temc-0.7_rp
          e_stick(k) = 10._rp**w1
       end do
    else
       do k = ks, ke
          ! Lin et al. (1983)
          temc_m     = min(tem(k) - t00,0.0_rp) ! T < 273.15
          e_stick(k) = exp(0.09_rp*temc_m)
       end do
    end if
 
    do k = ks, ke
       ! averaged diameter using SB06(82)
       ave_dc = coef_d(i_mp_qc)*xq(k,i_mp_qc)**b_m(i_mp_qc)
       !------------------------------------------------------------------------
       ! coellection efficiency are given as follows
       e_c(k) = max(0.0_rp, min(1.0_rp, (ave_dc-dc0)/(dc1-dc0) ))
    end do
 
    !------------------------------------------------------------------------
    ! Collection:  a collects b ( assuming particle size a>b )
    do k = ks, ke
       dcdc(k) = dq_xave(k,i_mp_qc) * dq_xave(k,i_mp_qc)
       drdr(k) = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qr)
       didi(k) = dq_xave(k,i_mp_qi) * dq_xave(k,i_mp_qi)
       dsds(k) = dq_xave(k,i_mp_qs) * dq_xave(k,i_mp_qs)
       dgdg(k) = dq_xave(k,i_mp_qg) * dq_xave(k,i_mp_qg)
       drdi(k) = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qi)
       drds(k) = dq_xave(k,i_mp_qr) * dq_xave(k,i_mp_qs)
    end do
    do k = ks, ke
       vcvc(k) = vt_xave(k,i_mp_qc,2) * vt_xave(k,i_mp_qc,2)
       vrvr(k) = vt_xave(k,i_mp_qr,2) * vt_xave(k,i_mp_qr,2)
       vivi(k) = vt_xave(k,i_mp_qi,2) * vt_xave(k,i_mp_qi,2)
       vsvs(k) = vt_xave(k,i_mp_qs,2) * vt_xave(k,i_mp_qs,2)
       vgvg(k) = vt_xave(k,i_mp_qg,2) * vt_xave(k,i_mp_qg,2)
       vrvi(k) = vt_xave(k,i_mp_qr,2) * vt_xave(k,i_mp_qi,2)
       vrvs(k) = vt_xave(k,i_mp_qr,2) * vt_xave(k,i_mp_qs,2)
    end do
 
    do k = ks, ke
       ave_di(k) = coef_d(i_mp_qi)*xq(k,i_mp_qi)**b_m(i_mp_qi)
       ave_ds(k) = coef_d(i_mp_qs)*xq(k,i_mp_qs)**b_m(i_mp_qs)
    end do
 
    !------------------------------------------------------------------------
    !
    !+++ pattern 1: a + b => a  (a>b)
    !                           (i-c, s-c, g-c, s-i, g-r, s-g)
    !------------------------------------------------------------------------
 
    ! cloud-ice => ice
    ! reduction term of cloud
    do k = ks, ke
       dcdi = dq_xave(k,i_mp_qc)   * dq_xave(k,i_mp_qi)
       vcvi = vt_xave(k,i_mp_qc,2) * vt_xave(k,i_mp_qi,2)
       sw = 0.5_rp - sign(0.5_rp, di0-ave_di(k)) ! if(ave_di>di0)then sw=1
       e_i = e_im * sw
       e_ic = e_i*e_c(k)
       coef_acc_lci = &
              ( delta_b1(i_mp_qc)*dcdc(k) + delta_ab1(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi(k) ) &
        * sqrt( theta_b1(i_mp_qc)*vcvc(k) - theta_ab1(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi(k) &
            +  sigma_i + sigma_c )
       coef_acc_nci = &
              ( delta_b0(i_mp_qc)*dcdc(k) + delta_ab0(i_mp_qi,i_mp_qc)*dcdi + delta_b0(i_mp_qi)*didi(k) ) &
        * sqrt( theta_b0(i_mp_qc)*vcvc(k) - theta_ab0(i_mp_qi,i_mp_qc)*vcvi + theta_b0(i_mp_qi)*vivi(k) &
            +  sigma_i + sigma_c )
       pac(k,i_liaclc2li)= -0.25_rp*pi*e_ic*rhoq(k,i_ni)*rhoq(k,i_qc)*coef_acc_lci
       pac(k,i_niacnc2ni)= -0.25_rp*pi*e_ic*rhoq(k,i_ni)*rhoq(k,i_nc)*coef_acc_nci
    end do
 
    ! cloud-snow => snow
    ! reduction term of cloud
    do k = ks, ke
       dcds = dq_xave(k,i_mp_qc)   * dq_xave(k,i_mp_qs)
       vcvs = vt_xave(k,i_mp_qc,2) * vt_xave(k,i_mp_qs,2)
       sw = 0.5_rp - sign(0.5_rp, ds0-ave_ds(k)) ! if(ave_ds>ds0)then sw=1
       e_s = e_sm * sw
       e_sc = e_s*e_c(k)
       coef_acc_lcs = &
              ( delta_b1(i_mp_qc)*dcdc(k) + delta_ab1(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds(k) ) &
        * sqrt( theta_b1(i_mp_qc)*vcvc(k) - theta_ab1(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs(k) &
            +  sigma_s + sigma_c )
       coef_acc_ncs = &
              ( delta_b0(i_mp_qc)*dcdc(k) + delta_ab0(i_mp_qs,i_mp_qc)*dcds + delta_b0(i_mp_qs)*dsds(k) ) &
        * sqrt( theta_b0(i_mp_qc)*vcvc(k) - theta_ab0(i_mp_qs,i_mp_qc)*vcvs + theta_b0(i_mp_qs)*vsvs(k) &
            +  sigma_s + sigma_c )
       pac(k,i_lsaclc2ls)= -0.25_rp*pi*e_sc*rhoq(k,i_ns)*rhoq(k,i_qc)*coef_acc_lcs
       pac(k,i_nsacnc2ns)= -0.25_rp*pi*e_sc*rhoq(k,i_ns)*rhoq(k,i_nc)*coef_acc_ncs
    end do
 
    ! cloud-graupel => graupel
    ! reduction term of cloud
    do k = ks, ke
       dcdg = dq_xave(k,i_mp_qc)   * dq_xave(k,i_mp_qg)
       vcvg = vt_xave(k,i_mp_qc,2) * vt_xave(k,i_mp_qg,2)
       ave_dg = coef_d(i_mp_qg)*xq(k,i_mp_qg)**b_m(i_mp_qg)
       sw = 0.5_rp - sign(0.5_rp, dg0-ave_dg) ! if(ave_dg>dg0)then sw=1
       e_g = e_gm * sw
       e_gc = e_g*e_c(k)
       coef_acc_lcg = &
              ( delta_b1(i_mp_qc)*dcdc(k) + delta_ab1(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg(k) ) &
        * sqrt( theta_b1(i_mp_qc)*vcvc(k) - theta_ab1(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg(k) &
            +  sigma_g + sigma_c )
       coef_acc_ncg = &
              ( delta_b0(i_mp_qc)*dcdc(k) + delta_ab0(i_mp_qg,i_mp_qc)*dcdg + delta_b0(i_mp_qg)*dgdg(k) ) &
        * sqrt( theta_b0(i_mp_qc)*vcvc(k) - theta_ab0(i_mp_qg,i_mp_qc)*vcvg + theta_b0(i_mp_qg)*vgvg(k) &
            +  sigma_g + sigma_c )
       pac(k,i_lgaclc2lg)= -0.25_rp*pi*e_gc*rhoq(k,i_ng)*rhoq(k,i_qc)*coef_acc_lcg
       pac(k,i_ngacnc2ng)= -0.25_rp*pi*e_gc*rhoq(k,i_ng)*rhoq(k,i_nc)*coef_acc_ncg
    end do
 
    ! snow-graupel => graupel
    do k = ks, ke
       dsdg = dq_xave(k,i_mp_qs)   * dq_xave(k,i_mp_qg)
       vsvg = vt_xave(k,i_mp_qs,2) * vt_xave(k,i_mp_qg,2)
       coef_acc_lsg = &
              ( delta_b1(i_mp_qs)*dsds(k) + delta_ab1(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg(k) ) &
        * sqrt( theta_b1(i_mp_qs)*vsvs(k) - theta_ab1(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg(k) &
            +  sigma_g + sigma_s )
       coef_acc_nsg(k) = &
              ( delta_b0(i_mp_qs)*dsds(k) + delta_ab0(i_mp_qg,i_mp_qs)*dsdg + delta_b0(i_mp_qg)*dgdg(k) ) &
            ! [fix] T.Mitsui 08/05/08
        * sqrt( theta_b0(i_mp_qs)*vsvs(k) - theta_ab0(i_mp_qg,i_mp_qs)*vsvg + theta_b0(i_mp_qg)*vgvg(k) &
            +  sigma_g + sigma_s )
       pac(k,i_lgacls2lg)= -0.25_rp*pi*e_stick(k)*e_gs*rhoq(k,i_ng)*rhoq(k,i_qs)*coef_acc_lsg
       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)
    end do
 
    !-----------------
    !  (start) Y.Sato added on 2018/8/31
    !--- ice-graupel => graupel
!!$    do k = KS, KE
!!$       didg = dq_xave(k,I_mp_QI)   * dq_xave(k,I_mp_QG)
!!$       vivg = vt_xave(k,I_mp_QI,2) * vt_xave(k,I_mp_QG,2)
!!$       coef_acc_LIG = &
!!$              ( delta_b1(I_QI)*didi(k) + delta_ab1(I_QG,I_QI)*didg + delta_b0(I_QG)*dgdg(k) ) &
!!$        * sqrt( theta_b1(I_QI)*vivi(k) - theta_ab1(I_QG,I_QI)*vivg + theta_b0(I_QG)*vgvg(k)   &
!!$            +  sigma_g + sigma_i )
!!$       coef_acc_NIG(k) = &
!!$              ( delta_b0(I_QI)*didi(k) + delta_ab0(I_QG,I_QI)*didg + delta_b0(I_QG)*dgdg(k) ) &
!!$        * sqrt( theta_b0(I_QI)*vivi(k) - theta_ab0(I_QG,I_QI)*vivg + theta_b0(I_QG)*vgvg(k) &
!!$            +  sigma_g + sigma_i )
!!$       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
!!$       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
!!$    end do
 
    !------------------------------------------------------------------------
    ! ice-snow => snow
    ! reduction term of ice
    do k = ks, ke
       dids = dq_xave(k,i_mp_qi)   * dq_xave(k,i_mp_qs)
       vivs = vt_xave(k,i_mp_qi,2) * vt_xave(k,i_mp_qs,2)
       coef_acc_lis = &
              ( delta_b1(i_mp_qi)*didi(k) + delta_ab1(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds(k) ) &
        * sqrt( theta_b1(i_mp_qi)*vivi(k) - theta_ab1(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs(k) &
            +  sigma_i + sigma_s )
       coef_acc_nis = &
              ( delta_b0(i_mp_qi)*didi(k) + delta_ab0(i_mp_qs,i_mp_qi)*dids + delta_b0(i_mp_qs)*dsds(k) ) &
        * sqrt( theta_b0(i_mp_qi)*vivi(k) - theta_ab0(i_mp_qs,i_mp_qi)*vivs + theta_b0(i_mp_qs)*vsvs(k) &
            +  sigma_i + sigma_s )
       pac(k,i_liacls2ls)= -0.25_rp*pi*e_stick(k)*e_si*rhoq(k,i_ns)*rhoq(k,i_qi)*coef_acc_lis
       pac(k,i_niacns2ns)= -0.25_rp*pi*e_stick(k)*e_si*rhoq(k,i_ns)*rhoq(k,i_ni)*coef_acc_nis
    end do
 
    do k = ks, ke
       drdg = dq_xave(k,i_mp_qr)   * dq_xave(k,i_mp_qg)
       vrvg = vt_xave(k,i_mp_qr,2) * vt_xave(k,i_mp_qg,2)
       sw = sign(0.5_rp, t00-tem(k)) + 0.5_rp
       ! if ( tem(k) <= T00 )then
          ! rain-graupel => graupel
          ! reduction term of rain
          ! sw = 1
       ! else
          ! rain-graupel => rain
          ! reduction term of graupel
          ! sw = 0
       coef_acc_lrg = &
              ( ( delta_b1(i_mp_qr)*drdr(k) + delta_ab1(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg(k) ) * sw &
              + ( 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) ) &
            * 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 &
                  + ( 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) &
                  + sigma_r + sigma_g )
       pac(k,i_lraclg2lg) = -0.25_rp*pi*e_gr*coef_acc_lrg * ( rhoq(k,i_ng)*rhoq(k,i_qr) * sw &
                                                            + rhoq(k,i_nr)*rhoq(k,i_qg) * (1.0_rp-sw) )
       coef_acc_nrg = &
              ( delta_b0(i_mp_qr)*drdr(k) + delta_ab0(i_mp_qg,i_mp_qr)*drdg + delta_b0(i_mp_qg)*dgdg(k) ) &
        * sqrt( theta_b0(i_mp_qr)*vrvr(k) - theta_ab0(i_mp_qg,i_mp_qr)*vrvg + theta_b0(i_mp_qg)*vgvg(k) &
            +  sigma_r + sigma_g )
       pac(k,i_nracng2ng) = -0.25_rp*pi*e_gr*rhoq(k,i_ng)*rhoq(k,i_nr)*coef_acc_nrg
    end do
 
    !------------------------------------------------------------------------
    !
    !+++ pattern 2: a + b => c  (a>b)
    !                           (r-i,r-s)
    !------------------------------------------------------------------------
 
    ! rain-ice => graupel
    ! reduction term of ice
    do k = ks, ke
       coef_acc_lri_i = &
              ( delta_b1(i_mp_qi)*didi(k) + delta_ab1(i_mp_qr,i_mp_qi)*drdi(k) + delta_b0(i_mp_qr)*drdr(k) ) &
        * 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) &
            +  sigma_r + sigma_i )
       coef_acc_nri_i = &
              ( delta_b0(i_mp_qi)*didi(k) + delta_ab0(i_mp_qr,i_mp_qi)*drdi(k) + delta_b0(i_mp_qr)*drdr(k) ) &
        * 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) &
            +  sigma_r + sigma_i )
       pac(k,i_lracli2lg_i)= -0.25_rp*pi*e_ir*rhoq(k,i_nr)*rhoq(k,i_qi)*coef_acc_lri_i
       pac(k,i_nracni2ng_i)= -0.25_rp*pi*e_ir*rhoq(k,i_nr)*rhoq(k,i_ni)*coef_acc_nri_i
    end do
 
    ! reduction term of rain
    do k = ks, ke
       coef_acc_lri_r = &
              ( delta_b1(i_mp_qr)*drdr(k) + delta_ab1(i_mp_qi,i_mp_qr)*drdi(k) + delta_b0(i_mp_qi)*didi(k) ) &
        * 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) &
            +  sigma_r + sigma_i )
       coef_acc_nri_r = &
              ( delta_b0(i_mp_qr)*drdr(k) + delta_ab0(i_mp_qi,i_mp_qr)*drdi(k) + delta_b0(i_mp_qi)*didi(k) ) &
        * 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) &
            +  sigma_r + sigma_i )
       pac(k,i_lracli2lg_r)= -0.25_rp*pi*e_ir*rhoq(k,i_ni)*rhoq(k,i_qr)*coef_acc_lri_r
       pac(k,i_nracni2ng_r)= -0.25_rp*pi*e_ir*rhoq(k,i_ni)*rhoq(k,i_nr)*coef_acc_nri_r
    end do
 
    ! rain-snow => graupel
    ! reduction term of snow
    do k = ks, ke
       coef_acc_lrs_s = &
              ( delta_b1(i_mp_qs)*dsds(k) + delta_ab1(i_mp_qr,i_mp_qs)*drds(k) + delta_b0(i_mp_qr)*drdr(k) ) &
        * 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) &
            +  sigma_r + sigma_s )
       coef_acc_nrs_s = &
              ( delta_b0(i_mp_qs)*dsds(k) + delta_ab0(i_mp_qr,i_mp_qs)*drds(k) + delta_b0(i_mp_qr)*drdr(k) ) &
        * 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) &
            +  sigma_r + sigma_s )
       pac(k,i_lracls2lg_s)= -0.25_rp*pi*e_sr*rhoq(k,i_nr)*rhoq(k,i_qs)*coef_acc_lrs_s
       pac(k,i_nracns2ng_s)= -0.25_rp*pi*e_sr*rhoq(k,i_nr)*rhoq(k,i_ns)*coef_acc_nrs_s
    end do
 
    ! reduction term of rain
    do k = ks, ke
       coef_acc_lrs_r = &
              ( delta_b1(i_mp_qr)*drdr(k) + delta_ab1(i_mp_qs,i_mp_qr)*drds(k) + delta_b0(i_mp_qs)*dsds(k) ) &
        * 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) &
            +  sigma_r + sigma_s )
       coef_acc_nrs_r = &
              ( delta_b0(i_mp_qr)*drdr(k) + delta_ab0(i_mp_qs,i_mp_qr)*drds(k) + delta_b0(i_mp_qs)*dsds(k) ) &
        * 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) &
            +  sigma_r + sigma_s )
       pac(k,i_lracls2lg_r)= -0.25_rp*pi*e_sr*rhoq(k,i_ns)*rhoq(k,i_qr)*coef_acc_lrs_r
       pac(k,i_nracns2ng_r)= -0.25_rp*pi*e_sr*rhoq(k,i_ns)*rhoq(k,i_nr)*coef_acc_nrs_r
    end do
 
    !------------------------------------------------------------------------
    !
    !+++ pattern 3: a + a => b  (i-i)
    !
    !------------------------------------------------------------------------
 
    ! ice-ice ( reduction is double, but includes double-count)
    do k = ks, ke
       coef_acc_lii = &
              ( delta_b0(i_mp_qi)*didi(k) + delta_ab1(i_mp_qi,i_mp_qi)*didi(k) + delta_b1(i_mp_qi)*didi(k) ) &
        * 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) &
            +  sigma_i + sigma_i )
       coef_acc_nii = &
              ( delta_b0(i_mp_qi)*didi(k) + delta_ab0(i_mp_qi,i_mp_qi)*didi(k) + delta_b0(i_mp_qi)*didi(k) ) &
        * 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) &
            +  sigma_i + sigma_i )
       pac(k,i_liacli2ls)= -0.25_rp*pi*e_stick(k)*e_ii*rhoq(k,i_ni)*rhoq(k,i_qi)*coef_acc_lii
       pac(k,i_niacni2ns)= -0.25_rp*pi*e_stick(k)*e_ii*rhoq(k,i_ni)*rhoq(k,i_ni)*coef_acc_nii
 
!          ci_aut(k)   =  0.25_RP*pi*E_ii*rhoq(k,I_NI)*coef_acc_LII
!          taui_aut(k) = 1._RP/max(E_stick(k)*ci_aut(k),1.E-10_RP)
!          tau_sce(k)  = rhoq(k,I_QI)/max(rhoq(k,I_QIj)+rhoq(k,I_QS),1.E-10_RP)
    end do
 
    !------------------------------------------------------------------------
    !
    !+++ pattern 4: a + a => a  (s-s)
    !
    !------------------------------------------------------------------------
 
    ! snow-snow => snow
    do k = ks, ke
       coef_acc_nss = &
              ( delta_b0(i_mp_qs)*dsds(k) + delta_ab0(i_mp_qs,i_mp_qs)*dsds(k) + delta_b0(i_mp_qs)*dsds(k) ) &
        * 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) &
            +  sigma_s + sigma_s )
       pac(k,i_nsacns2ns)= -0.125_rp*pi*e_stick(k)*e_ss*rhoq(k,i_ns)*rhoq(k,i_ns)*coef_acc_nss
    end do
 
    ! graupel-grauple => graupel
    do k = ks, ke
       coef_acc_ngg = &
              ( delta_b0(i_mp_qg)*dgdg(k) + delta_ab0(i_mp_qg,i_mp_qg)*dgdg(k) + delta_b0(i_mp_qg)*dgdg(k) ) &
        * 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) &
            +  sigma_g + sigma_g )
       pac(k,i_ngacng2ng)= -0.125_rp*pi*e_stick(k)*e_gg*rhoq(k,i_ng)*rhoq(k,i_ng)*coef_acc_ngg
    end do
 
    !------------------------------------------------------------------------
    !--- Partial conversion
    ! SB06(70),(71)
    ! i_iconv2g: option whether partial conversions work or not
 
    ! ice-cloud => graupel
    do k = ks, ke
       sw = 0.5_rp - sign(0.5_rp,di_cri-ave_di(k)) ! if( ave_di > di_cri )then sw=1
       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
       pq(k,i_licon) = i_iconv2g * pac(k,i_liaclc2li)/max(1.0_rp, wx_cri) * sw
       pq(k,i_nicon) = i_iconv2g * pq(k,i_licon)/xq(k,i_mp_qi) * sw
    end do
 
    ! snow-cloud => graupel
    do k = ks, ke
       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 )
       pq(k,i_lscon) = i_sconv2g * (pac(k,i_lsaclc2ls))/max(1.0_rp, wx_crs)
       pq(k,i_nscon) = i_sconv2g * pq(k,i_lscon)/xq(k,i_mp_qs)
    end do
 
    !--- enhanced melting( due to collection-freezing of water droplets )
    !    originally from Rutledge and Hobbs(1984). eq.(A.21)
    ! if T > 273.15 then temc_p=T-273.15, else temc_p=0
    ! 08/05/08 [fix] T.Mitsui LHF00 => LHF0
    ! melting occurs around T=273K, so LHF0 is suitable both SIMPLE and EXACT,
    ! otherwise LHF can have sign both negative(EXACT) and positive(SIMPLE).
    do k = ks, ke
!       temc_m = min(tem(k) - T00, 0.0_RP) ! T < 273.15
       temc_p = max(tem(k) - t00, 0.0_rp) ! T > 273.15
!!$       coef_emelt   = -CL/LHF00*temc_p
       coef_emelt   =  cl/lhf0*temc_p
       ! cloud-graupel
       pq(k,i_lgacm) =  coef_emelt*pac(k,i_lgaclc2lg)
       pq(k,i_ngacm) =  pq(k,i_lgacm)/xq(k,i_mp_qg)
       ! rain-graupel
       pq(k,i_lgarm) =  coef_emelt*pac(k,i_lraclg2lg)
       pq(k,i_ngarm) =  pq(k,i_lgarm)/xq(k,i_mp_qg)
       ! cloud-snow
       pq(k,i_lsacm) =  coef_emelt*(pac(k,i_lsaclc2ls))
       pq(k,i_nsacm) =  pq(k,i_lsacm)/xq(k,i_mp_qs)
       ! rain-snow
       pq(k,i_lsarm) =  coef_emelt*(pac(k,i_lracls2lg_r)+pac(k,i_lracls2lg_s))
       pq(k,i_nsarm) =  pq(k,i_lsarm)/xq(k,i_mp_qg)
       ! cloud-ice
       pq(k,i_liacm) =  coef_emelt*pac(k,i_liaclc2li)
       pq(k,i_niacm) =  pq(k,i_liacm)/xq(k,i_mp_qi)
       ! rain-ice
       pq(k,i_liarm) =  coef_emelt*(pac(k,i_lracli2lg_r)+pac(k,i_lracli2lg_i))
       pq(k,i_niarm) =  pq(k,i_liarm)/xq(k,i_mp_qg)
    end do
 
 
    !---- for charge density
    if ( flg_lt ) then
       !--- C + I -> I (decrease from cloud charge)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small ) !--- if NC is small,  ignore charge transfer
          pcrg2(k,i_niacnc2ni) = pac(k,i_niacnc2ni)*(1.0_rp-sw1) / (rhoq(k,i_nc)+sw1) * rhoq_crg(k,i_qc)
       end do
 
       !--- C + S -> S (decrease from cloud charge)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small ) !--- if NC is small,  ignore charge transfer
          pcrg2(k,i_nsacnc2ns) = pac(k,i_nsacnc2ns)*(1.0_rp-sw1) / (rhoq(k,i_nc)+sw1) * rhoq_crg(k,i_qc)
       end do
 
       !--- C + G -> G (decrease from cloud charge)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small ) !--- if NC is small,  ignore charge transfer
          pcrg2(k,i_ngacnc2ng) = pac(k,i_ngacnc2ng)*(1.0_rp-sw1) / (rhoq(k,i_nc)+sw1) * rhoq_crg(k,i_qc)
       end do
 
       !--- S + G -> G (decrease from snow charge)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NS is small,  ignore charge transfer
          pcrg2(k,i_ngacns2ng) = pac(k,i_ngacns2ng)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
       end do
 
       !--- Charge separation by Snow-Graupel rebound--------------------------------
       do k = ks, ke
          alpha_lt = 5.0_rp * ( dq_xave(k,i_mp_qs) / d0_crg )**2 * vt_xave(k,i_mp_qg,2) / v0_crg
          alpha_lt = min( alpha_lt, 10.0_rp )
          pcrg2(k,i_cgngacns2ng)= 0.25_rp*pi*( 1.0_rp - e_stick(k) )*e_gs &
                                * rhoq(k,i_ng)*rhoq(k,i_ns)*coef_acc_nsg(k) &
                                * ( dqcrg(k)*alpha_lt ) &
                                * beta_crg(k)
       end do
 
       !--- I + G -> G (decrease from snow charge)
!!$       do k = KS, KE
!!$         sw1 = 0.5_RP - sign( 0.5_RP, rhoq(k,I_NI)-SMALL ) !--- if NS is small,  ignore charge transfer
!!$         Pcrg2(k,I_NGacNI2NG) = Pac(k,I_NGacNI2NG)*(1.0_RP-sw1) / (rhoq(k,I_NI)+sw1) * rhoq_crg(k,I_QI) * flg_igcol
!!$       !--- Charge separation by Ice-Graupel rebound--------------------------------
!!$         alpha_lt = 5.0_RP * ( dq_xave(k,I_mp_QI) / d0_crg )**2 * vt_xave(k,I_mp_QG,2) / v0_crg
!!$         alpha_lt = min( alpha_lt, 10.0_RP )
!!$         Pcrg2(k,I_CGNGacNI2NG)= 0.25_RP*pi*( 1.0_RP - E_stick(k) )*E_gi &
!!$                                  * rhoq(k,I_NG)*rhoq(k,I_NI)*coef_acc_NIG(k) &
!!$                                  * ( dqcrg(k)*alpha_lt ) &
!!$                                  * beta_crg(k) * flg_igcol
!!$       end do
 
       !--- I + S -> S (decrease from ice charge)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small,  ignore charge transfer
          pcrg2(k,i_niacns2ns) = pac(k,i_niacns2ns)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
       end do
 
       !--- R+G->R (T>T00 sw=0, decrease from graupel charge), ->G(T<=T00 sw=1, dcrerase from rain charge)
       do k = ks, ke
          sw = 0.5_rp + sign( 0.5_rp, t00-tem(k) )
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small ) !--- if NR is small,  ignore charge transfer
          sw2 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small ) !--- if NG is small,  ignore charge transfer
          pcrg2(k,i_nracng2ng) = pac(k,i_nracng2ng)*(1.0_rp-sw1)/(rhoq(k,i_nr)+sw1) * rhoq_crg(k,i_qr) * sw &
                               + pac(k,i_nracng2ng)*(1.0_rp-sw2)/(rhoq(k,i_ng)+sw2) * rhoq_crg(k,i_qg) * (1.0_rp-sw)
       end do
 
       !--- R + I -> G (decrease from both ice and rain charge, but only ice charge at this part)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small,  ignore charge transfer
          pcrg2(k,i_nracni2ng_i) = pac(k,i_nracni2ng_i)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
       end do
 
       !--- R + I -> G (decrease from both ice and rain charge, but only rain charge at this part)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small ) !--- if NR is small,  ignore charge transfer
          pcrg2(k,i_nracni2ng_r) = pac(k,i_nracni2ng_r)*(1.0_rp-sw1) / (rhoq(k,i_nr)+sw1) * rhoq_crg(k,i_qr)
       end do
 
       !--- R + S -> G (decrease from both snow and rain charge, but only snow charge at this part)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NS is small,  ignore charge transfer
          pcrg2(k,i_nracns2ng_s) = pac(k,i_nracns2ng_s)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
       end do
 
       !--- R + S -> G (decrease from both snow and rain charge, but only rain charge at this part)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small ) !--- if NR is small,  ignore charge transfer
          pcrg2(k,i_nracns2ng_r) = pac(k,i_nracns2ng_r)*(1.0_rp-sw1) / (rhoq(k,i_nr)+sw1) * rhoq_crg(k,i_qr)
       end do
 
       !--- I + I -> S (decrease from ice charge)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small,  ignore charge transfer
          pcrg2(k,i_niacni2ns) = pac(k,i_niacni2ns)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
       end do
 
       !--- I + C -> G (decrease from ice charge)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small,  ignore charge transfer
          pcrg1(k,i_nicon) = i_iconv2g * pq(k,i_nicon)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
       end do
 
       !--- S + C -> G (decrease from snow charge)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NS is small,  ignore charge transfer
          pcrg1(k,i_nscon) = i_sconv2g * pq(k,i_nscon)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
       end do
 
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small ) !--- if NG is small,  ignore charge transfer
          pcrg1(k,i_ngacm) = pq(k,i_ngacm)*(1.0_rp-sw1) / (rhoq(k,i_ng)+sw1) * rhoq_crg(k,i_qg)
       end do
 
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small ) !--- if NG is small,  ignore charge transfer
          pcrg1(k,i_ngarm) = pq(k,i_ngarm)*(1.0_rp-sw1) / (rhoq(k,i_ng)+sw1) * rhoq_crg(k,i_qg)
       end do
 
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NS is small,  ignore charge transfer
          pcrg1(k,i_nsacm) = pq(k,i_nsacm)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
       end do
 
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small ) !--- if NS is small,  ignore charge transfer
          pcrg1(k,i_nsarm) = pq(k,i_nsarm)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
       end do
 
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small,  ignore charge transfer
          pcrg1(k,i_niacm) = pq(k,i_niacm)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
       end do
 
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small ) !--- if NI is small,  ignore charge transfer
          pcrg1(k,i_niarm) = pq(k,i_niarm)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
       end do
 
    end if
 
    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 4898 of file scale_atmos_phy_mp_sn14.F90.

    use scale_atmos_saturation, only: &
       moist_psat_ice => atmos_saturation_psat_ice
    use scale_prc, only: &
       prc_abort
    implicit none
 
    integer, intent(in) :: KA, KS, KE
 
    !--- mixed-phase collection process
    !                  And all we set all production term as a negative sign to avoid confusion.
    !
    real(RP), intent(in) :: wtem(KA)
    !--- mass/number concentration[kg/m3]
    real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
    real(RP), intent(in) :: rho(KA) ! air density
    ! necessary ?
    real(RP), intent(in) :: xq(KA,HYDRO_MAX)
    !--- diameter of averaged mass( D(ave x) )
    real(RP), intent(in) :: dq_xave(KA,HYDRO_MAX)
    !--- terminal velocity of averaged mass( vt(ave x) )
    real(RP), intent(in) :: vt_xave(KA,HYDRO_MAX,2)
    ! [Add] 11/08/30 T.Mitsui, for autoconversion of ice
    !    real(RP), intent(in) :: rho(KA)
    !--- partial conversion
    real(RP), intent(inout):: PQ(KA,PQ_MAX)
    !
    real(RP), intent(out):: Pac(KA,Pac_MAX)
    !--- for lightning component
    logical,  intent(in) :: flg_lt
    real(RP), intent(in) :: beta_crg(KA)
    real(RP), intent(in) :: dqcrg(KA)
    real(RP), intent(in) :: d0_crg, v0_crg
    real(RP), intent(in) :: rhoq_crg(KA,I_QC:I_QG)
    real(RP), intent(inout):: Pcrg1(KA,PQ_MAX)
    real(RP), intent(inout):: Pcrg2(KA,Pcrg_MAX)
 
    real(RP), parameter :: a_dec = 0.883_rp
    real(RP), parameter :: b_dec = 0.093_rp
    real(RP), parameter :: c_dec = 0.00348_rp
    real(RP), parameter :: d_dec = 4.5185e-5_rp
    !
    !!!!from default SN14
 
    real(RP), parameter :: E_C12(6)=(/&
         0.010_rp, 0.080_rp, 0.10_rp, 0.60_rp, 0.20_rp, 0.10_rp/)
 
    real(RP) :: tem(KA)
    !
    !--- collection efficency of each specie
    real(RP) :: E_c, E_r, E_i, E_s, E_g
!    real(RP) ::  E_ic, E_sc, E_gc
    !--- sticking efficiency
    real(RP):: E_stick(KA)
    ! [Add] 10/08/03 T.Mitsui
    real(RP) :: temc, temc2, temc3
    real(RP) :: E_dec
    real(RP) :: esi_rat(KA)
    real(RP) :: esi(KA)
    !
    real(RP) :: temc_p, temc_m             ! celcius tem.
!    real(RP) :: ci_aut(KA)
!    real(RP) :: taui_aut(KA)
!    real(RP) :: tau_sce(KA)
    !--- DSD averaged diameter for each species
    real(RP) :: ave_dc                     ! cloud
!    real(RP) :: ave_dr                     ! rain
    real(RP) :: ave_di(KA)                 ! ice
    real(RP) :: ave_ds(KA)                 ! snow
    real(RP) :: ave_dg                     ! graupel
    !--- coefficient of collection equations(L:mass, N:number)
    real(RP) :: coef_acc_LCI,   coef_acc_NCI   ! cloud     - cloud ice
    real(RP) :: coef_acc_LCS,   coef_acc_NCS   ! cloud     - snow
    !
    real(RP) :: coef_acc_LCG,   coef_acc_NCG   ! cloud     - graupel
    real(RP) :: coef_acc_LRI_I, coef_acc_NRI_I ! rain      - cloud ice
    real(RP) :: coef_acc_LRI_R, coef_acc_NRI_R ! rain      - cloud ice
    real(RP) :: coef_acc_LRS_S, coef_acc_NRS_S ! rain      - snow
    real(RP) :: coef_acc_LRS_R, coef_acc_NRS_R ! rain      - snow
    real(RP) :: coef_acc_LRG,   coef_acc_NRG   ! rain      - graupel
    real(RP) :: coef_acc_LII,   coef_acc_NII   ! cloud ice - cloud ice
    real(RP) :: coef_acc_LIS,   coef_acc_NIS   ! cloud ice - snow
    real(RP) ::                 coef_acc_NSS   ! snow      - snow
    real(RP) ::                 coef_acc_NGG   ! grauepl   - graupel
    real(RP) :: coef_acc_LSG,   coef_acc_NSG(KA) ! snow      - graupel
    !--- (diameter) x (diameter)
    real(RP) :: dcdc(KA), dcdi, dcds, dcdg
    real(RP) :: drdr(KA), drdi(KA), drds(KA), drdg
    real(RP) :: didi(KA), dids, didg
    real(RP) :: dsds(KA), dsdg
    real(RP) :: dgdg(KA)
!#    !--- (terminal velocity) x (terminal velocity)
!#    real(RP) :: vcvc(KA), vcvi, vcvs, vcvg
!#    real(RP) :: vrvr(KA), vrvi(KA), vrvs(KA), vrvg
!#    real(RP) :: vivi(KA), vivs, vivg
!#    real(RP) :: vsvs(KA), vsvg
!#    real(RP) :: vgvg(KA)
    !
!    real(RP) :: wx_cri, wx_crs
!    real(RP) :: coef_emelt
!    real(RP) :: w1
 
    real(RP) :: sw, sw1, sw2, tmp
    real(RP) :: alpha_lt
    !
    integer  :: k, iqw
    !
 
    real(RP) :: tem_e(KA) ! [Add] 15/05/19 T.Seiki
 
 
    !
    ! work for binary collision
    !
    ! work for Gauss-Legendre quadrature
    integer, parameter :: ngmax=4
 
    real(RP) :: lambdac(KA), lambdar(KA), lambdai(KA), lambdas(KA), lambdag(KA)
    real(RP) :: A_dsdc(KA), A_dsdr(KA), A_dsdi(KA), A_dsds(KA), A_dsdg(KA)
    real(RP) :: dNdx
    real(RP) :: dxdd
    real(RP) :: dNdD
    real(RP) :: dNc_glx(KA,ngmax), dNr_glx(KA,ngmax), dNi_glx(KA,ngmax), dNs_glx(KA,ngmax), dNg_glx(KA,ngmax)
    real(RP) :: dNc_gly, dNr_gly(KA,ngmax), dNi_gly(KA,ngmax), dNs_gly(KA,ngmax), dNg_gly(KA,ngmax)
    !
    real(RP) :: dc_glx, dr_glx, di_glx, ds_glx, dg_glx
    real(RP) :: dc_gly, dr_gly, di_gly, ds_gly, dg_gly
    real(RP) :: xc_glx(KA,ngmax), xr_glx(KA,ngmax), xi_glx(KA,ngmax), xs_glx(KA,ngmax), xg_glx(KA,ngmax)
    real(RP) :: xc_gly, xr_gly(KA,ngmax), xi_gly, xs_gly, xg_gly
 
    real(RP) :: vtc_glx(KA,ngmax), vtr_glx(KA,ngmax), vti_glx(KA,ngmax), vts_glx(KA,ngmax), vtg_glx(KA,ngmax)
    real(RP) :: vtc_gly, vtr_gly(KA,ngmax), vti_gly(KA,ngmax), vts_gly(KA,ngmax), vtg_gly(KA,ngmax)
    real(RP) :: dac_glx(KA,ngmax), dar_glx(KA,ngmax), dai_glx(KA,ngmax), das_glx(KA,ngmax), dag_glx(KA,ngmax)
    real(RP) :: dac_gly, dar_gly(KA,ngmax), dai_gly(KA,ngmax), das_gly(KA,ngmax), dag_gly(KA,ngmax)
    !
    integer :: ngx, ngy
    !
    real(RP) :: E_ic(KA), E_sc(KA), E_gc(KA)
    !
    ! Parameters to calculate terminal velocity formulated by Mitchell (1996)
    !
    real(RP) :: acx, bcx, gcx, scx ! geometric parameters of hexagonal column ice defined by Mitchell (1996)
    real(RP) :: acy, bcy, gcy, scy ! geometric parameters of hexagonal column ice defined by Mitchell (1996)
    real(RP), parameter :: as = 0.59452551_rp
    real(RP), parameter :: bs = 2.4490_rp
    real(RP), parameter :: gs = 0.131488_rp
    real(RP), parameter :: ss = 1.880000_rp
    real(RP), parameter :: ag = 19.5072514_rp !0.049d0*1.d-3*(100.0d0**2.8d0)
    real(RP), parameter :: bg = 2.8_rp
    real(RP), parameter :: gg = 0.5_rp
    real(RP), parameter :: sg = 2.0_rp
    real(RP) :: num_Besti_glx, num_Bests_glx, num_Bestg_glx
    real(RP) :: num_Besti_gly, num_Bests_gly, num_Bestg_gly
    real(RP) :: num_Rei_glx, num_Res_glx, num_Reg_glx
    real(RP) :: num_Rei_gly, num_Res_gly, num_Reg_gly
    real(RP), parameter :: c0=0.6_rp  ! Bohm (1989)
    real(RP), parameter :: d0=5.83_rp ! Bohm (1989)
    !
    real(RP) :: mua(KA), nua(KA)
    !--- Dynamic viscosity
    real(RP), parameter :: mua0 = 1.718e-5_rp
    
    real(RP), parameter :: dmua_dT = 5.28e-8_rp
    
    !======  mua = mua0 + temc*dmua_dT
    !
    ! collection Kernel
    !
    real(RP) :: kernel_cg, kernel_cs, kernel_ci
    real(RP) :: kernel_rg, kernel_rs, kernel_ri
    real(RP) :: kernel_ig, kernel_is, kernel_ii
    real(RP) :: kernel_sg, kernel_ss
    real(RP) :: kernel_gg
    real(RP) :: kernel_sg_reb
    !
    !--------------------------
    !
    ! accurate for PSDs with optimization
    !
    !-------------------------------
    ! cloud
    ! gauss_range = 2.0_RP
    ! ngmax = 4
    ! rain,
    ! gauss_range = xxx
    ! ngmax = 4
    ! ice, snow, graupel
    ! gauss_range = 5.0_RP
    ! ngmax = 4
    !
    real(RP), parameter :: gauss_rangec=2.0_rp
    real(RP), parameter :: wc_gl(ngmax)=(/&
         0.2411146051511425e+00_rp, &
         0.4520325754088027e+00_rp, &
         0.4520325754088027e+00_rp, &
         0.2411146051511425e+00_rp &
         /)
    real(RP), parameter :: coefc_d_gl(ngmax)=(/&
         0.5505187813766612e+00_rp, &
         0.7900516927471093e+00_rp, &
         0.1265739962562290e+01_rp, &
         0.1816468454535444e+01_rp &
         /)
    real(RP), parameter :: gauss_ranger=8.0_rp
    real(RP), parameter :: wr_gl(ngmax)=(/&
         0.723343815453428e+00_rp,&
         1.35609772622641e+00_rp, &
         1.35609772622641e+00_rp, &
         0.723343815453428e+00_rp &
         /)
    real(RP), parameter :: coefr_d_gl(ngmax)=(/&
         0.166846238310235e+00_rp, &
         0.493135790663523e+00_rp, &
         2.02783902311062e+00_rp, &
         5.99354237846583e+00_rp &
         /)
    real(RP), parameter :: gauss_range=5.0_rp
    real(RP), parameter :: w_gl(ngmax)=(/&
         0.559850775788111e+00_rp, &
         1.04958713664599e+00_rp, &
         1.04958713664599e+00_rp, &
         0.559850775788111e+00_rp &
         /)
    real(RP), parameter :: coef_d_gl(ngmax)=(/&
         0.2500872485877803e+00_rp, &
         0.5785800417604080e+00_rp, &
         0.1728369331505741e+01_rp, &
         0.3998604509613778e+01_rp &
         /)
    !
    real(RP) :: wx_cri, wx_crs
    real(RP) :: coef_emelt
    real(RP) :: w1
!    integer :: ij, k
    integer :: ierr
 
    !---------------------------------------------------------------------------
    !
    !
 
    do k = ks, ke
       tem(k) = max( wtem(k), tem_min ) ! 11/08/30 T.Mitsui
    end do
 
    call moist_psat_ice( ka, ks, ke, &
                         tem(:), esi(:)  ) ! [IN], [OUT]
 
    if( opt_stick_ks96 )then
       do k = ks, ke
          ! Khain and Sednev (1996), eq.(3.15)
          temc       = tem(k) - t00
          temc2      = temc*temc
          temc3      = temc2*temc
          e_dec      = max(0.0_rp, a_dec + b_dec*temc + c_dec*temc2 + d_dec*temc3 )
          esi_rat(k)    = rhoq(k,i_qv)*rvap*tem(k)/esi(k)
          e_stick(k) = min(1.0_rp, e_dec*esi_rat(k))
       end do
    else if( opt_stick_co86 )then
       do k = ks, ke
          ! [Add] 11/08/30 T.Mitsui, Cotton et al. (1986)
          temc       = min(tem(k) - t00,0.0_rp)
          w1         = 0.035_rp*temc-0.7_rp
          e_stick(k) = 10._rp**w1
       end do
    else
       do k = ks, ke
          ! Lin et al. (1983)
          temc_m     = min(tem(k) - t00,0.0_rp) ! T < 273.15
          e_stick(k) = exp(0.09_rp*temc_m)
       end do
    end if
 
    !------------------------------------------------------------------------
 
 
    pac(:,:) = 0.0_rp
 
    tem(:) = max(wtem(:), tem_min )
    tem_e(:) = max(wtem(:), tem_min_estick )
 
    call moist_psat_ice( ka, ks, ke, tem(:), esi(:) )
 
 
    if ( opt_stick_ks96 ) then
       do k=ks, ke
          ! Khain and Sednev (1996), eq.(3.15)
          temc          = tem_e(k) - t00 ![Mod] 15/05/19 T.Seiki
          temc2         = temc*temc
          temc3         = temc2*temc
          e_dec         = max(0.0_rp, a_dec + b_dec*temc + c_dec*temc2 + d_dec*temc3 )
          esi_rat(k) = rhoq(k,i_qv)*rvap*tem(k)/esi(k)
          e_stick(k) = min(1.0_rp, e_dec*esi_rat(k))
       enddo
    elseif( opt_stick_co86 ) then
       do k=ks, ke
          temc          = min(tem_e(k) - t00,0.0_rp)
          w1            = 0.035_rp*temc-0.7_rp
          e_stick(k) = 10.0_rp**w1
       enddo
    elseif( opt_stick_c12 ) then
       do k=ks, ke
          if     (tem_e(k)>273.15_rp) then
             e_stick(k)=1.0_rp
          elseif(tem_e(k)<243.15_rp) then  !-30degC
             e_stick(k)=e_c12(1)
          elseif(tem_e(k)<248.15_rp) then  !-25degC
             e_stick(k)=e_c12(1)+0.2_rp*(e_c12(2)-e_c12(1))*(tem_e(k)-243.15_rp)
          elseif(tem_e(k)<253.15_rp) then  !-20degC
             e_stick(k)=e_c12(2)+0.2_rp*(e_c12(3)-e_c12(2))*(tem_e(k)-248.15_rp)
          elseif(tem_e(k)<258.15_rp) then  !-15degC
             e_stick(k)=e_c12(3)+0.2_rp*(e_c12(4)-e_c12(3))*(tem_e(k)-253.15_rp)
          elseif(tem_e(k)<263.15_rp) then  !-10degC
             e_stick(k)=e_c12(4)+0.2_rp*(e_c12(5)-e_c12(4))*(tem_e(k)-258.15_rp)
          elseif(tem_e(k)<268.15_rp) then  !- 5degC
             e_stick(k)=e_c12(5)+0.2_rp*(e_c12(6)-e_c12(5))*(tem_e(k)-263.15_rp)
          else                               !  0degC
             e_stick(k)=e_c12(6)
          endif
       enddo
    else
       if ( opt_stick_rhh57 ) then
          do k=ks, ke
             if ( tem_e(k) < 270.0_rp .AND. rhoq(k,i_qv)*rvap*tem(k) < esi(k) ) then
                esi_rat(k) = 0.0_rp
             else
                esi_rat(k) = 1.0_rp
             endif
          enddo
       elseif( opt_stick_rhks96 ) then
          do k=ks, ke
             esi_rat(k) = min( rhoq(k,i_qv)*rvap*tem(k)/esi(k),1.0_rp )
          enddo
       else
          esi_rat(:) = 1.0_rp
       endif
       !
       do k=ks, ke
          ! Lin et al. (1983)
          temc_m        = min(tem_e(k) - t00,0.0_rp) ! T < 273.15
          e_stick(k) = exp(0.09_rp*temc_m)*esi_rat(k)
       enddo
    endif
 
    !
    ! Integration Start
    !
 
    do k = ks, ke
       mua(k) = mua0 + dmua_dt*(tem(k)-273.15_rp)
       nua(k) = mua(k)/rho(k)        ! [m2/s]
    end do
    do k = ks, ke
       lambdac(k) = xq(k,i_mp_qc)**(-mu(i_mp_qc))*coef_lambda(i_mp_qc)
       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))
    end do
    do k = ks, ke
       lambdar(k) = xq(k,i_mp_qr)**(-mu(i_mp_qr))*coef_lambda(i_mp_qr)
       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))
    end do
    do k = ks, ke
       lambdai(k) = xq(k,i_mp_qi)**(-mu(i_mp_qi))*coef_lambda(i_mp_qi)
       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))
    end do
    do k = ks, ke
       lambdas(k) = xq(k,i_mp_qs)**(-mu(i_mp_qs))*coef_lambda(i_mp_qs)
       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))
    end do
    do k = ks, ke
       lambdag(k) = xq(k,i_mp_qg)**(-mu(i_mp_qg))*coef_lambda(i_mp_qg)
       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))
    end do
 
    !
    ! Particle X
    !
    do ngx=1, ngmax ! X < Y
!OCL LOOP_FISSION_TARGET(LS)
       do k = ks, ke
          dc_glx         = dq_xave(k,i_mp_qc)*coefc_d_gl(ngx)
          xc_glx(k,ngx) = ( (dc_glx/a_m(i_mp_qc)) )**(1.0_rp/b_m(i_mp_qc))
          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))&
               *(xc_glx(k,ngx)/(b_m(i_mp_qc)*dc_glx))*dc_glx*wc_gl(ngx) ! dNdlogD*weight
          !
          dr_glx         = dq_xave(k,i_mp_qr)*coefr_d_gl(ngx)
          xr_glx(k,ngx) = ( (dr_glx/a_m(i_mp_qr)) )**(1.0_rp/b_m(i_mp_qr))
          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))&
               *(xr_glx(k,ngx)/(b_m(i_mp_qr)*dr_glx))*dr_glx*wr_gl(ngx) ! dNdlogD*weight
          !
          di_glx         = dq_xave(k,i_mp_qi)*coef_d_gl(ngx)
          xi_glx(k,ngx) = ( (di_glx/a_m(i_mp_qi)) )**(1.0_rp/b_m(i_mp_qi))
          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))&
               *(xi_glx(k,ngx)/(b_m(i_mp_qi)*di_glx))*di_glx*w_gl(ngx) ! dNdlogD*weight
          !
          ds_glx         = dq_xave(k,i_mp_qs)*coef_d_gl(ngx)
          xs_glx(k,ngx) = ( (ds_glx/a_m(i_mp_qs)) )**(1.0_rp/b_m(i_mp_qs))
          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))&
               *(xs_glx(k,ngx)/(b_m(i_mp_qs)*ds_glx))*ds_glx*w_gl(ngx)! dNdlogD*weight
          !
          dg_glx         = dq_xave(k,i_mp_qg)*coef_d_gl(ngx)
          xg_glx(k,ngx) = ( (dg_glx/a_m(i_mp_qg)) )**(1.0_rp/b_m(i_mp_qg))
          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))&
               *(xg_glx(k,ngx)/(b_m(i_mp_qg)*dg_glx))*dg_glx*w_gl(ngx)! dNdlogD*weight
          !
          ! Hexagonal Columns
!!$          if( di_glx <= 100.e-6_RP )then
!!$             acx = 0.1677_RP*1.e-3_RP*(100.0_RP**2.91_RP)
!!$             bcx = 2.91_RP
!!$             gcx = (0.684_RP*1.e-4_RP)*10.0_RP**(2.0_RP*2.0_RP)
!!$             scx = 2.0_RP
!!$          else
!!$             acx = 0.00166_RP*1.e-3_RP*(100.0_RP**1.91_RP)
!!$             bcx = 1.91_RP
!!$             gcx = (0.0696_RP*1.e-4_RP)*10.0_RP**(2.0_RP*1.5_RP)
!!$             scx = 1.5_RP
!!$          end if
          sw = 0.5_rp + sign(0.5_rp, di_glx - 100.e-6_rp )
          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 )
          bcx = 2.91_rp*(1.0_rp-sw) + 1.91_rp*sw
          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 )
          scx = 2.0_rp*(1.0_rp-sw) + 1.5_rp*sw
          num_besti_glx = 2.0_rp*acx*grav*rho(k)*di_glx**(bcx+2.0_rp-scx)/(gcx*mua(k)**2)
          num_bests_glx = 2.0_rp*as *grav*rho(k)*ds_glx**(bs +2.0_rp-ss )/(gs *mua(k)**2)
          num_bestg_glx = 2.0_rp*ag *grav*rho(k)*dg_glx**(bg +2.0_rp-sg )/(gg *mua(k)**2)
          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
          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
          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
          !
          vtc_glx(k,ngx) = coef_vtr_ar2*dc_glx*(1.0_rp-exp(-coef_vtr_br2*dc_glx))
          !
!!$          if( dr_glx < d_vtr_branch )then
!!$             vtr_glx = coef_vtr_ar2*dr_glx*(1.0_RP-exp(-coef_vtr_br2*dr_glx))
!!$          else
!!$             vtr_glx = coef_vtr_ar1-coef_vtr_br1*exp(-coef_vtr_cr1*dr_glx)
!!$          end if
          sw = 0.5_rp + sign( 0.5_rp, dr_glx - d_vtr_branch )
          tmp = exp( - ( coef_vtr_br2*(1.0_rp-sw) + coef_vtr_cr1*sw ) * dr_glx )
          vtr_glx(k,ngx) = coef_vtr_ar2 * dr_glx * ( 1.0_rp - tmp ) * (1.0_rp-sw) &
                         + ( coef_vtr_ar1 - coef_vtr_br1 * tmp ) * sw
          !
          vti_glx(k,ngx) = num_rei_glx*nua(k)/di_glx
          vts_glx(k,ngx) = num_res_glx*nua(k)/ds_glx
          vtg_glx(k,ngx) = num_reg_glx*nua(k)/dg_glx
          ! equivalent area diameter
          dac_glx(k,ngx) = dc_glx
          dar_glx(k,ngx) = dr_glx
          dai_glx(k,ngx) = 2.0_rp*sqrt( (gcx*di_glx**scx)/pi )
          das_glx(k,ngx) = 2.0_rp*sqrt( (gs *ds_glx**ss )/pi )
          dag_glx(k,ngx) = 2.0_rp*sqrt( (gg *dg_glx**sg )/pi )
       end do
    end do
 
    !
    ! Particle Y
    !
    do ngy=1, ngmax
!OCL LOOP_FISSION_TARGET(LS)
       do k = ks, ke
          dr_gly         = dq_xave(k,i_mp_qr)*coefr_d_gl(ngy)
          xr_gly(k,ngy) = ( (dr_gly/a_m(i_mp_qr)) )**(1.0_rp/b_m(i_mp_qr))
          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))&
               *(xr_gly(k,ngy)/(b_m(i_mp_qr)*dr_gly))*dr_gly*wr_gl(ngy) ! dNdlogD*weight
          !
          di_gly         = dq_xave(k,i_mp_qi)*coef_d_gl(ngy)
          xi_gly         = ( (di_gly/a_m(i_mp_qi)) )**(1.0_rp/b_m(i_mp_qi))
          dni_gly(k,ngy) = a_dsdi(k)*(xi_gly**nu(i_mp_qi)) * exp(-lambdai(k)*xi_gly**mu(i_mp_qi))&
               *(xi_gly/(b_m(i_mp_qi)*di_gly))*di_gly*w_gl(ngy) ! dNdlogD*weight
          !
          ds_gly         = dq_xave(k,i_mp_qs)*coef_d_gl(ngy)
          xs_gly         = ( (ds_gly/a_m(i_mp_qs)) )**(1.0_rp/b_m(i_mp_qs))
          dns_gly(k,ngy) = a_dsds(k)*(xs_gly**nu(i_mp_qs)) * exp(-lambdas(k)*xs_gly**mu(i_mp_qs))&
               *(xs_gly/(b_m(i_mp_qs)*ds_gly))*ds_gly*w_gl(ngy)! dNdlogD*weight
          !
          dg_gly         = dq_xave(k,i_mp_qg)*coef_d_gl(ngy)
          xg_gly         = ( (dg_gly/a_m(i_mp_qg)) )**(1.0_rp/b_m(i_mp_qg))
          dng_gly(k,ngy) = a_dsdg(k)*(xg_gly**nu(i_mp_qg)) * exp(-lambdag(k)*xg_gly**mu(i_mp_qg))&
               *(xg_gly/(b_m(i_mp_qg)*dg_gly))*dg_gly*w_gl(ngy)! dNdlogD*weight
          !
          ! Hexagonal Columns
!!$          if( di_gly <= 100.e-6_RP )then
!!$             acy = 0.1677_RP*1.d-3*(100.0_RP**2.91_RP)
!!$             bcy = 2.91_RP
!!$             gcy = (0.684_RP*1.e-4_RP)*10.0_RP**(2.0_RP*2.0_RP)
!!$             scy = 2.0_RP
!!$          else
!!$             acy = 0.00166_RP*1.e-3_RP*(100.0_RP**1.91_RP)
!!$             bcy = 1.91_RP
!!$             gcy = (0.0696_RP*1.e-4_RP)*10.0_RP**(2.0_RP*1.5_RP)
!!$             scy = 1.5_RP
!!$          end if
          sw = 0.5_rp + sign( 0.5_rp, di_gly - 100.e-6_rp )
          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 )
          bcy = 2.91_rp*(1.0_rp-sw) + 1.91_rp*sw
          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 )
          scy = 2.0_rp*(1.0_rp-sw) + 1.5_rp*sw
          num_besti_gly = 2.0_rp*acy*grav*rho(k)*di_gly**(bcy+2.0_rp-scy)/(gcy*mua(k)**2)
          num_bests_gly = 2.0_rp*as *grav*rho(k)*ds_gly**(bs +2.0_rp-ss )/(gs *mua(k)**2)
          num_bestg_gly = 2.0_rp*ag *grav*rho(k)*dg_gly**(bg +2.0_rp-sg )/(gg *mua(k)**2)
          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
          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
          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
          !
!!$          if( dr_gly < d_vtr_branch )then
!!$             vtr_gly(k,ngy) = coef_vtr_ar2*dr_gly*(1.0_RP-exp(-coef_vtr_br2*dr_gly))
!!$          else
!!$             vtr_gly(k,ngy) = coef_vtr_ar1-coef_vtr_br1*exp(-coef_vtr_cr1*dr_gly)
!!$          end if
          sw = 0.5_rp + sign( 0.5_rp, dr_gly - d_vtr_branch )
          tmp = exp( - ( coef_vtr_br2*(1.0_rp-sw) + coef_vtr_cr1*sw ) * dr_gly )
          vtr_gly(k,ngy) = coef_vtr_ar2 * dr_gly * ( 1.0_rp - tmp ) * (1.0_rp-sw) &
                         + ( coef_vtr_ar1 - coef_vtr_br1 * tmp ) * sw
          !
          vti_gly(k,ngy) = num_rei_gly*nua(k)/di_gly
          vts_gly(k,ngy) = num_res_gly*nua(k)/ds_gly
          vtg_gly(k,ngy) = num_reg_gly*nua(k)/dg_gly
          ! equivalent area diameter
          dar_gly(k,ngy) = dr_gly
          dai_gly(k,ngy) = 2.0_rp*sqrt( (gcy*di_gly**scy)/pi )
          das_gly(k,ngy) = 2.0_rp*sqrt( (gs *ds_gly**ss )/pi )
          dag_gly(k,ngy) = 2.0_rp*sqrt( (gg *dg_gly**sg )/pi )
       end do
    end do
 
    !
    ! BULK collection efficiency are given as follows
    !
    do k = ks, ke
       e_c = max(0.0_rp, min(1.0_rp, (dq_xave(k,i_mp_qc)-dc0)/(dc1-dc0) ))
       !
!!$       if (dq_xave(k,I_mp_QI)>di0) then
!!$          E_i = E_im
!!$       else
!!$          E_i = 0.0_RP
!!$       endif
       sw = 0.5_rp + sign( 0.5_rp, dq_xave(k,i_mp_qi)-di0 )
       e_i = e_im * sw
!!$       if (dq_xave(k,I_mp_QS)>ds0) then
!!$          E_s = E_sm
!!$       else
!!$          E_s = 0.0_RP
!!$       endif
       sw = 0.5_rp + sign( 0.5_rp, dq_xave(k,i_mp_qs)-ds0 )
       e_s = e_sm * sw
!!$       if (dq_xave(k,I_mp_QG)>dg0) then
!!$          E_g = E_gm
!!$       else
!!$          E_g = 0.0_RP
!!$       endif
       sw = 0.5_rp + sign( 0.5_rp, dq_xave(k,i_mp_qg)-dg0 )
       e_g = e_gm * sw
 
       e_ic(k) = e_i*e_c
       e_sc(k) = e_s*e_c
       e_gc(k) = e_g*e_c
    end do
 
 
    !=========================================================================================
    ! collection equation
    !=========================================================================================
    do ngx=1, ngmax ! X < Y
       do ngy=1, ngmax
 
          do k = ks, ke
             !
             ! 1.c-g (X=Cloud, Y=Graupel)
             !
!             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)
             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)
             pac(k,i_ngacnc2ng) = pac(k,i_ngacnc2ng) - kernel_cg              *dnc_glx(k,ngx)*dng_gly(k,ngy)
             pac(k,i_lgaclc2lg) = pac(k,i_lgaclc2lg) - kernel_cg*xc_glx(k,ngx)*dnc_glx(k,ngx)*dng_gly(k,ngy)
          end do
          do k = ks, ke
             !
             ! 2.c-s (X=Cloud, Y=Snow)
             !
!             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)
             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)
             pac(k,i_nsacnc2ns) = pac(k,i_nsacnc2ns) - kernel_cs              *dnc_glx(k,ngx)*dns_gly(k,ngy)
             pac(k,i_lsaclc2ls) = pac(k,i_lsaclc2ls) - kernel_cs*xc_glx(k,ngx)*dnc_glx(k,ngx)*dns_gly(k,ngy)
          end do
          do k = ks, ke
             !
             ! 3.c-i (X=Cloud, Y=Cloud Ice)
             !
!             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)
             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)
             pac(k,i_niacnc2ni) = pac(k,i_niacnc2ni) - kernel_ci              *dnc_glx(k,ngx)*dni_gly(k,ngy)
             pac(k,i_liaclc2li) = pac(k,i_liaclc2li) - kernel_ci*xc_glx(k,ngx)*dnc_glx(k,ngx)*dni_gly(k,ngy)
          end do
          do k = ks, ke
             !
             ! 4.r-g
             !
!             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
             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
             ! T < 273K  (X=Rain   , Y=Graupel)
             pac(k,i_nracng2ng) = pac(k,i_nracng2ng) - kernel_rg              *dnr_glx(k,ngx)*dng_gly(k,ngy)
             pac(k,i_lraclg2lg) = pac(k,i_lraclg2lg) - kernel_rg*xr_glx(k,ngx)*dnr_glx(k,ngx)*dng_gly(k,ngy)
             ! T > 273K  (X=Graupel, Y=Rain   )
             pac(k,i_nracng2nr) = pac(k,i_nracng2nr) - kernel_rg              *dng_glx(k,ngx)*dnr_gly(k,ngy)
             pac(k,i_lraclg2lr) = pac(k,i_lraclg2lr) - kernel_rg*xg_glx(k,ngx)*dng_glx(k,ngx)*dnr_gly(k,ngy)
          end do
          do k = ks, ke
             !
             ! 5.r-s
             !
!             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
             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
             ! (X=Snow, Y=Rain)
             pac(k,i_nracns2ng_r) = pac(k,i_nracns2ng_r) - kernel_rs              *dns_glx(k,ngx)*dnr_gly(k,ngy)
             pac(k,i_lracls2lg_r) = pac(k,i_lracls2lg_r) - kernel_rs*xr_gly(k,ngy)*dns_glx(k,ngx)*dnr_gly(k,ngy)
             !
             pac(k,i_nracns2ng_s) = pac(k,i_nracns2ng_s) - kernel_rs              *dns_glx(k,ngx)*dnr_gly(k,ngy)
             pac(k,i_lracls2lg_s) = pac(k,i_lracls2lg_s) - kernel_rs*xs_glx(k,ngx)*dns_glx(k,ngx)*dnr_gly(k,ngy)
          end do
          do k = ks, ke
             !
             ! 6.r-i
             !
!             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
             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
             ! (X=Cloud Ice, Y=Rain)
             pac(k,i_nracni2ng_r) = pac(k,i_nracni2ng_r) - kernel_ri              *dni_glx(k,ngx)*dnr_gly(k,ngy)
             pac(k,i_lracli2lg_r) = pac(k,i_lracli2lg_r) - kernel_ri*xr_gly(k,ngy)*dni_glx(k,ngx)*dnr_gly(k,ngy)
             !
             pac(k,i_nracni2ng_i) = pac(k,i_nracni2ng_i) - kernel_ri              *dni_glx(k,ngx)*dnr_gly(k,ngy)
             pac(k,i_lracli2lg_i) = pac(k,i_lracli2lg_i) - kernel_ri*xi_glx(k,ngx)*dni_glx(k,ngx)*dnr_gly(k,ngy)
          end do
          do k = ks, ke
             !
             ! 7.i-g
             !
!             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
             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
             pac(k,i_niacng2ng) = pac(k,i_niacng2ng) - kernel_ig              *dni_glx(k,ngx)*dng_gly(k,ngy)
             pac(k,i_liaclg2lg) = pac(k,i_liaclg2lg) - kernel_ig*xi_glx(k,ngx)*dni_glx(k,ngx)*dng_gly(k,ngy)
          end do
          do k = ks, ke
             !
             ! 8.i-s
             !
!             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
             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
             pac(k,i_niacns2ns) = pac(k,i_niacns2ns) - kernel_is              *dni_glx(k,ngx)*dns_gly(k,ngy)
             pac(k,i_liacls2ls) = pac(k,i_liacls2ls) - kernel_is*xi_glx(k,ngx)*dni_glx(k,ngx)*dns_gly(k,ngy)
          end do
          do k = ks, ke
             !
             ! 9.i-i
             !
!             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
             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
             pac(k,i_niacni2ns) = pac(k,i_niacni2ns) - kernel_ii              *dni_glx(k,ngx)*dni_gly(k,ngy)
             pac(k,i_liacli2ls) = pac(k,i_liacli2ls) - kernel_ii*xi_glx(k,ngx)*dni_glx(k,ngx)*dni_gly(k,ngy)
          end do
          do k = ks, ke
             !
             ! 10.s-g
             !
!             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
             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
             pac(k,i_ngacns2ng) = pac(k,i_ngacns2ng) - kernel_sg              *dns_glx(k,ngx)*dng_gly(k,ngy)
             pac(k,i_lgacls2lg) = pac(k,i_lgacls2lg) - kernel_sg*xs_glx(k,ngx)*dns_glx(k,ngx)*dng_gly(k,ngy)
          end do
          do k = ks, ke
             !
             ! 11.s-s
             !
!             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
             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
             pac(k,i_nsacns2ns) = pac(k,i_nsacns2ns) - kernel_ss*dns_glx(k,ngx)*dns_gly(k,ngy)
          end do
          do k = ks, ke
             !
             ! 12.g-g
             !
!             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
             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
             pac(k,i_ngacng2ng) = pac(k,i_ngacng2ng) - kernel_gg*dng_glx(k,ngx)*dng_gly(k,ngy)
 
          end do
       end do
    end do
 
    !--- Charge separation by Snow-Graupel rebound--------------------------------
    if ( flg_lt ) then
       do k = ks, ke
!OCL UNROLL('full')
          do ngy = 1, ngmax
!OCL UNROLL('full')
             do ngx = 1, ngmax
                alpha_lt = 5.0_rp * ( das_glx(k,ngx) / d0_crg )**2*vtg_gly(k,ngy)/v0_crg
                alpha_lt = min( alpha_lt, 10.0_rp )
!                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)))  &
                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)) &
                              * ( 1.0_rp - e_stick(k) ) * e_gs
                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)
             end do
          end do
       end do
    end if
 
    !
    !
    do k=ks, ke
       temc_p = max(tem(k) - t00,0.0_rp) ! T > 273.15
       !------------------------------------------------------------------------
       !--- Partial conversion
       ! SB06(70),(71)
       ! i_iconv2g: option whether partial conversions work or not
       ! ice-cloud => graupel
       if ( dq_xave(k,i_mp_qi) > di_cri ) then
          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 )
          pq(k,i_licon) = i_iconv2g*  pac(k,i_liaclc2li)/max(1.0_rp, wx_cri)
          pq(k,i_nicon) = i_iconv2g*  pq(k,i_licon)/xq(k,i_mp_qi)
       else
          wx_cri       = 0.0_rp
          pq(k,i_licon) = 0.0_rp
          pq(k,i_nicon) = 0.0_rp
       endif
       ! snow-cloud => graupel
       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 )
       pq(k,i_lscon) = i_sconv2g*  (pac(k,i_lsaclc2ls))/max(1.0_rp, wx_crs)
       pq(k,i_nscon) = i_sconv2g*  pq(k,i_lscon)/xq(k,i_mp_qs)
       !------------------------------------------------------------------------
       !--- enhanced melting( due to collection-freezing of water droplets )
       !    originally from Rutledge and Hobbs(1984). eq.(A.21)
       ! if T > 273.15 then temc_p=T-273.15, else temc_p=0
       ! 08/05/08 [fix] T.Mitsui LHF00 => LHF0
       ! melting occurs around T=273K, so LHF0 is suitable both SIMPLE and EXACT,
       ! otherwise LHF can have sign both negative(EXACT) and positive(SIMPLE).
       coef_emelt   =  cl/lhf0*temc_p
       ! cloud-graupel
       pq(k,i_lgacm) =  coef_emelt*pac(k,i_lgaclc2lg)
       pq(k,i_ngacm) =  pq(k,i_lgacm)/xq(k,i_mp_qg)
       ! rain-graupel
       pq(k,i_lgarm) =  coef_emelt*pac(k,i_lraclg2lg)
       pq(k,i_ngarm) =  pq(k,i_lgarm)/xq(k,i_mp_qg)
       ! cloud-snow
       pq(k,i_lsacm) =  coef_emelt*(pac(k,i_lsaclc2ls))
       pq(k,i_nsacm) =  pq(k,i_lsacm)/xq(k,i_mp_qs)
       ! rain-snow
       pq(k,i_lsarm) =  coef_emelt*(pac(k,i_lracls2lg_r)+pac(k,i_lracls2lg_s))
       pq(k,i_nsarm) =  pq(k,i_lsarm)/xq(k,i_mp_qg)
       ! cloud-ice
       pq(k,i_liacm) =  coef_emelt*pac(k,i_liaclc2li)
       pq(k,i_niacm) =  pq(k,i_liacm)/xq(k,i_mp_qi)
       ! rain-ice
       pq(k,i_liarm) =  coef_emelt*(pac(k,i_lracli2lg_r)+pac(k,i_lracli2lg_i))
       pq(k,i_niarm) =  pq(k,i_liarm)/xq(k,i_mp_qg)
    enddo
 
    !---- for charge density
    if ( flg_lt ) then
 
       ! 1.c-g (X=Cloud, Y=Graupel)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
          pcrg2(k,i_ngacnc2ng) = pac(k,i_ngacnc2ng)*(1.0_rp-sw1) / (rhoq(k,i_nc)+sw1) * rhoq_crg(k,i_qc)
       end do
 
       ! 2.c-s (X=Cloud, Y=Snow)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
          pcrg2(k,i_nsacnc2ns) = pac(k,i_nsacnc2ns)*(1.0_rp-sw1) / (rhoq(k,i_nc)+sw1) * rhoq_crg(k,i_qc)
       end do
 
       ! 3.c-i (X=Cloud, Y=Cloud Ice)
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
          pcrg2(k,i_niacnc2ni) = pac(k,i_niacnc2ni)*(1.0_rp-sw1) / (rhoq(k,i_nc)+sw1) * rhoq_crg(k,i_qc)
       end do
 
       ! 4.r-g
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small )
          sw2 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small )
          pcrg2(k,i_nracng2ng) = pac(k,i_nracng2ng)*(1.0_rp-sw1) / (rhoq(k,i_nr)+sw1) * rhoq_crg(k,i_qr)
          pcrg2(k,i_nracng2nr) = pac(k,i_nracng2nr)*(1.0_rp-sw2) / (rhoq(k,i_ng)+sw2) * rhoq_crg(k,i_qg)
       end do
 
       ! 5.r-s
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small )
          sw2 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
          pcrg2(k,i_nracns2ng_r) = pac(k,i_nracns2ng_r)*(1.0_rp-sw1) / (rhoq(k,i_nr)+sw1) * rhoq_crg(k,i_qr)
          pcrg2(k,i_nracns2ng_s) = pac(k,i_nracns2ng_s)*(1.0_rp-sw2) / (rhoq(k,i_ns)+sw2) * rhoq_crg(k,i_qs)
       end do
 
 
       ! 6.r-i
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nr)-small )
          sw2 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
          pcrg2(k,i_nracni2ng_r) = pac(k,i_nracni2ng_r)*(1.0_rp-sw1) / (rhoq(k,i_nr)+sw1) * rhoq_crg(k,i_qr)
          pcrg2(k,i_nracni2ng_i) = pac(k,i_nracni2ng_i)*(1.0_rp-sw2) / (rhoq(k,i_ni)+sw2) * rhoq_crg(k,i_qi)
       end do
 
       ! 7.i-g
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
          pcrg2(k,i_niacng2ng) = pac(k,i_niacng2ng)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
       end do
 
       ! 8.i-s
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
          pcrg2(k,i_niacns2ns) = pac(k,i_niacns2ns)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
       end do
 
       ! 9.i-i
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
          pcrg2(k,i_niacni2ns) = pac(k,i_niacni2ns)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
       end do
 
       ! 10.s-g
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
          pcrg2(k,i_ngacns2ng) = pac(k,i_ngacns2ng)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
       end do
 
       ! 11.s-s
       do k = ks, ke
          pcrg2(k,i_nsacns2ns) = 0.0_rp ! no charge transfer between category due to the collection of the same category (snow-snow)
       end do
 
       ! 12.g-g
       do k = ks, ke
          pcrg2(k,i_ngacng2ng) = 0.0_rp ! no charge transfer between category due to the collection of the same category (graupel-graupel)
       end do
 
       !--- Partial conversion
       ! ice-cloud => graupel
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
          pcrg1(k,i_nicon) = pq(k,i_nicon)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
       end do
 
       ! snow-cloud => graupel
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
          pcrg1(k,i_nscon) = pq(k,i_nscon)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
       end do
 
       !--- enhanced melting( due to collection-freezing of water droplets )
       ! cloud-graupel
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small )
          pcrg1(k,i_ngacm) = pq(k,i_ngacm)*(1.0_rp-sw1) / (rhoq(k,i_ng)+sw1) * rhoq_crg(k,i_qg)
       end do
 
       ! rain-graupel
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ng)-small )
          pcrg1(k,i_ngarm) = pq(k,i_ngarm)*(1.0_rp-sw1) / (rhoq(k,i_ng)+sw1) * rhoq_crg(k,i_qg)
       end do
 
       ! cloud-snow
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
          pcrg1(k,i_nsacm) = pq(k,i_nsacm)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
       end do
 
       ! rain-snow
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ns)-small )
          pcrg1(k,i_nsarm) = pq(k,i_nsarm)*(1.0_rp-sw1) / (rhoq(k,i_ns)+sw1) * rhoq_crg(k,i_qs)
       end do
 
       ! cloud-ice
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
          pcrg1(k,i_niacm) = pq(k,i_niacm)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
       end do
 
       ! rain-ice
       do k = ks, ke
          sw1 = 0.5_rp - sign( 0.5_rp, rhoq(k,i_ni)-small )
          pcrg1(k,i_niarm) = pq(k,i_niarm)*(1.0_rp-sw1) / (rhoq(k,i_ni)+sw1) * rhoq_crg(k,i_qi)
       end do
 
    endif
 
    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 5799 of file scale_atmos_phy_mp_sn14.F90.

    implicit none
 
    integer, intent(in) :: KA, KS, KE
    !
    real(RP), intent(in)  :: rhoq(KA,I_QV:I_NG)
    real(RP), intent(in)  :: rhoq_crg(KA,I_QC:I_QG)
    logical,  intent(in)  :: flg_lt
    real(RP), intent(in)  :: xq(KA,HYDRO_MAX)
    real(RP), intent(in)  :: dq_xave(KA,HYDRO_MAX)
    real(RP), intent(in)  :: rho(KA)
    !
    real(RP), intent(inout) :: PQ(KA,PQ_MAX)
    real(RP), intent(inout) :: Pcrg(KA,PQ_MAX)
    !
    ! parameter for autoconversion
    real(RP), parameter :: kcc     = 4.44e+9_rp  ! collision efficiency [m3/kg2/sec]
    real(RP), parameter :: tau_min = 1.e-20_rp  ! empirical filter by T.Mitsui
    real(RP), parameter :: rx_sep  = 1.0_rp/x_sep ! 1/x_sep, 10/08/03 [Add] T.Mitsui
    !
    ! parameter for accretion
    real(RP), parameter :: kcr     = 5.8_rp     ! collision efficiency [m3/kg2/sec]
    real(RP), parameter :: thr_acc = 5.e-5_rp   ! threshold for universal function original
    !
    ! parameter for self collection and collison break-up
    real(RP), parameter :: krr     = 4.33_rp ! k_rr,      S08 (35)
    real(RP), parameter :: kaprr   = 60.7_rp  ! kappa_rr,  SB06(11)
    real(RP), parameter :: kbr     = 1000._rp ! k_br,      SB06(14)
    real(RP), parameter :: kapbr   = 2.3e+3_rp   ! kappa_br,  SB06(15)
    real(RP), parameter :: dr_min  = 0.35e-3_rp ! minimum diameter, SB06(13)-(15)
    !
    ! work variables
    real(RP) :: coef_nuc0 ! coefficient of number for Auto-conversion
    real(RP) :: coef_nuc1 !                mass
    real(RP) :: coef_aut0 !                number
    real(RP) :: coef_aut1 !                mass
    real(RP) :: lwc       ! lc+lr
    real(RP) :: tau       ! conversion ratio: qr/(qc+qr) ranges [0:1]
    real(RP) :: rho_fac   ! factor of air density
    real(RP) :: psi_aut   ! Universal function of Auto-conversion
    real(RP) :: psi_acc   ! Universal function of Accretion
    real(RP) :: psi_brk   ! Universal function of Breakup
    real(RP) :: ddr       ! diameter difference from equilibrium
    !
    integer :: k, iqw
    real(RP) :: sw
    !
    coef_nuc0 = (nu(i_mp_qc)+2.0_rp)/(nu(i_mp_qc)+1.0_rp)
    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)
    coef_aut0 =  -kcc*coef_nuc0
    coef_aut1 =  -kcc/x_sep/20._rp*coef_nuc1
    !
!OCL LOOP_FISSION_TARGET(LS)
    do k = ks, ke
       lwc = rhoq(k,i_qr) + rhoq(k,i_qc)
       if( lwc > xc_min )then
          tau  = max(tau_min, rhoq(k,i_qr)/lwc)
       else
          tau  = tau_min
       end if
       rho_fac = sqrt(rho_0/max(rho(k),rho_min))
       !
       ! Auto-conversion ( cloud-cloud => rain )
       psi_aut       = 400._rp*(tau**0.7_rp)*(1.0_rp - (tau**0.7_rp))**3   ! (6) SB06
       pq(k,i_ncaut)  = coef_aut0*rhoq(k,i_qc)*rhoq(k,i_qc)*rho_fac*rho_fac    ! (9) SB06 sc+aut
       ! lc = lwc*(1-tau), lr = lwc*tau
       pq(k,i_lcaut)  = coef_aut1*lwc*lwc*xq(k,i_mp_qc)*xq(k,i_mp_qc) &          ! (4) SB06
            *((1.0_rp-tau)*(1.0_rp-tau) + psi_aut)*rho_fac*rho_fac        !
       pq(k,i_nraut)  = -rx_sep*pq(k,i_lcaut)                           ! (A7) SB01
       !--- for Charge density
       if (flg_lt) then
          sw = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
          pcrg(k,i_ncaut) = pq(k,i_ncaut)*(1.0_rp-sw)/(rhoq(k,i_nc)+sw)*rhoq_crg(k,i_qc)
          pcrg(k,i_nraut) = -pcrg(k,i_ncaut)
       end if
       !
       ! Accretion ( cloud-rain => rain )
       psi_acc       =(tau/(tau+thr_acc))**4                          ! (8) SB06
       pq(k,i_lcacc)  = -kcr*rhoq(k,i_qc)*rhoq(k,i_qr)*rho_fac*psi_acc         ! (7) SB06
       pq(k,i_ncacc)  = -kcr*rhoq(k,i_nc)*rhoq(k,i_qr)*rho_fac*psi_acc         ! (A6) SB01
       !--- for Charge density
       if (flg_lt) then
          sw = 0.5_rp - sign( 0.5_rp, rhoq(k,i_nc)-small )
          pcrg(k,i_ncacc)  = pq(k,i_ncacc)*(1.0_rp-sw)/(rhoq(k,i_nc)+sw)*rhoq(k,i_qc)
       end if
       !
       ! Self-collection ( rain-rain => rain )
       pq(k,i_nrslc)  = -krr*rhoq(k,i_nr)*rhoq(k,i_qr)*rho_fac                 ! (A.8) SB01
       !
       ! Collisional breakup of rain
       ddr           = min(1.e-3_rp, dq_xave(k,i_mp_qr) - dr_eq )
       if      ( dq_xave(k,i_mp_qr) < dr_min ) then       ! negligible
          psi_brk      = -1.0_rp
       else if ( dq_xave(k,i_mp_qr) <= dr_eq  ) then
          psi_brk      = kbr*ddr                          ! (14) SB06
       else
          psi_brk      = exp(kapbr*ddr) - 1.0_rp          ! (15) SB06 (SB06 has a typo)
       end if
       pq(k,i_nrbrk) = - (psi_brk + 1.0_rp)*pq(k,i_nrslc) ! (13) SB06
       !
    end do
    !
    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 6151 of file scale_atmos_phy_mp_sn14.F90.

    use scale_const, only: &
         rvap    =>  const_rvap,      &
         rdry    =>  const_rdry,      &
         t00     =>  const_tem00,     &
         lhs0    =>  const_lhs0,      &
         lhf00   =>  const_lhf00,     &
         pi      =>  const_pi,        &
         cpdry   =>  const_cpdry,     &
         pstd    =>  const_pstd,      &
         psat0   =>  const_psat0
    implicit none
 
    integer, intent(in)  :: KA
    integer, intent(in)  :: KS
    integer, intent(in)  :: KE
    ! Diffusion growth or Evaporation, Sublimation
    real(RP), intent(out) :: PLIdep_total(KA)  ! mass          for cloud ice
    real(RP), intent(in)  :: rho(KA)     ! air density
    real(RP), intent(in)  :: tem(KA)     ! air temperature
    real(RP), intent(in)  :: pre(KA)     ! air pressure
    real(RP), intent(in)  :: qd(KA)      ! mixing ratio of dry air
    real(RP), intent(in)  :: esi(KA)     ! saturation vapor pressure(solid water)
    real(RP), intent(in)  :: qsi(KA)     ! saturation vapor pressure(solid water)
    real(RP), intent(in)  :: rhoq(KA,I_QV:I_NG)
    ! Notice following values differ from mean terminal velocity or diameter.
    ! mean(vt(x)) /= vt(mean(x)) and mean(D(x)) /= D(mean(x))
    ! Following ones are vt(mean(x)) and D(mean(x)).
    real(RP), intent(in)  :: vt_xave(KA,HYDRO_MAX,1:2) ! terminal velocity of mean cloud
    real(RP), intent(in)  :: dq_xave(KA,HYDRO_MAX) !
    real(RP), intent(in)  :: dt
    !
!    real(RP), intent(in)  :: dtime
!    real(RP) :: dtime
    !
    real(RP) :: rho_lim            ! limited density
    real(RP) :: temc_lim           ! limited temperature[celsius]
    real(RP) :: pre_lim            ! limited density
    real(RP) :: temc               ! temperature[celsius]
    real(RP) :: pv                 ! vapor pressure
    real(RP) :: qv                 ! vapor pressure
    real(RP) :: ssi                ! super saturation ratio(ice water)
    real(RP) :: nua, r_nua         ! kinematic viscosity of air
    real(RP) :: mua                ! viscosity of air
    real(RP) :: Kat                 ! thermal conductance
    real(RP) :: Dw                 ! diffusivity of water vapor
    real(RP) :: Dht                 ! diffusivity of heat
    real(RP) :: Gi                 ! diffusion factor by balance between heat and vapor
    real(RP) :: Gii, Gis, Gig      ! for rain, ice, snow and graupel.
    real(RP) :: Gm                 ! melting factor by balance between heat and vapor
    real(RP) :: Nsc_r3             !
    !
    real(RP) :: Nreis_r2  !
    real(RP) :: Nress_r2, Nregs_r2  !
    real(RP) :: Nreil_r2  !
    real(RP) :: Nresl_r2, Nregl_r2  !
    real(RP) :: NscNrei_s, NscNrei_l
    real(RP) :: NscNres_s, NscNres_l
    real(RP) :: NscNreg_s, NscNreg_l
    real(RP) :: ventLI_s, ventLI_l
    real(RP) :: ventLS_s, ventLS_l
    real(RP) :: ventLG_s, ventLG_l
    real(RP) :: wti, wts, wtg
    real(RP), parameter :: r_14=1.0_rp/1.4_rp
    real(RP), parameter :: r_15=1.0_rp/1.5_rp
    real(RP) :: ventLI !
    real(RP) :: ventLS !
    real(RP) :: ventLG !
    !
    real(RP) :: total_dep
    real(RP) :: PLIdep_wrk
    real(RP) :: PLSdep_wrk
    real(RP) :: PLGdep_wrk
    real(RP) :: dep_limiter
    !
    real(RP), parameter :: Re_max=1.e3_rp
    real(RP), parameter :: Re_min=1.e-4_rp
    integer :: k
    !
!    PLIdep_total(1:KS)=0.0_RP
!    PLIdep_total(KE:KA)=0.0_RP
    plidep_total(:)=0.0_rp  !! 22/10/14
 
    !
    do k=ks, ke
       temc    = tem(k) - t00   ! degC
       temc_lim= max(temc, temc_lim_diff) !
       rho_lim = max(rho(k),rho_min)   !
       qv      = rhoq(k,i_qv)/rho_lim
       pre_lim = rho_lim*(qd(k)*rdry + qv*rvap)*(temc_lim+t00)
       !--------------------------------------------------------------------
       ! Diffusion growth part is described in detail
       ! by Pruppacher and Klett (1997) Sec. 13.2(liquid) and 13.3(solid)
       ! G:factor of thermal diffusion(1st.term) and vapor diffusion(2nd. term)
       ! SB06(23),(38), Lin et al(31),(52) or others
       ! Dw is introduced by Pruppacher and Klett(1997),(13-3)
       dw      = 0.211e-4_rp* (((temc_lim+t00)/t00)**1.94_rp) *(pstd/pre_lim)
       kat     = ka0  + temc_lim*dka_dt
       mua     = mua0 + temc_lim*dmua_dt
       nua     = mua/rho_lim
       r_nua   = 1.0_rp/nua
       gi      = (lhs0/kat/tem(k))*(lhs0/rvap/tem(k)-1.0_rp)+(rvap*tem(k)/dw/esi(k))
       ! capacities account for their surface geometries
       gii     = 4.0_rp*pi/cap(i_mp_qi)/gi
       gis     = 4.0_rp*pi/cap(i_mp_qs)/gi
       gig     = 4.0_rp*pi/cap(i_mp_qg)/gi
       ! vent: ventilation effect( asymmetry vapor field around particles due to aerodynamic )
       ! SB06 (30),(31) and each coefficient is by (88),(89)
       nsc_r3  = (nua/dw)**(0.33333333_rp)                    ! (Schmidt number )^(1/3)
       ! Beard and Pruppacher(1971) had performed in the range [0<Re<=320],
       ! So here we limit Re in the range Re_max=1000, Re_min=0.0001.
       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
       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
       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
       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
       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
       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
       !
       nscnrei_s=nsc_r3*nreis_r2 ! small ice
       nscnrei_l=nsc_r3*nreil_r2 ! large ice
       nscnres_s=nsc_r3*nress_r2 ! small snow
       nscnres_l=nsc_r3*nresl_r2 ! large snow
       nscnreg_s=nsc_r3*nregs_r2 ! small snow
       nscnreg_l=nsc_r3*nregl_r2 ! large snow
       !
       ventli_s = ah_vent1(i_mp_qi,1) + bh_vent1(i_mp_qi,1)*nscnrei_s
       ventli_l = ah_vent1(i_mp_qi,2) + bh_vent1(i_mp_qi,2)*nscnrei_l
       ventls_s = ah_vent1(i_mp_qs,1) + bh_vent1(i_mp_qs,1)*nscnres_s
       ventls_l = ah_vent1(i_mp_qs,2) + bh_vent1(i_mp_qs,2)*nscnres_l
       ventlg_s = ah_vent1(i_mp_qg,1) + bh_vent1(i_mp_qg,1)*nscnreg_s
       ventlg_l = ah_vent1(i_mp_qg,2) + bh_vent1(i_mp_qg,2)*nscnreg_l
       ! branch is 1.4 for rain, snow, graupel; is 1.0 for ice (PK97, 13-60,-61,-88,-89).
       wti     = ( min(max( nscnrei_s     , 0.5_rp), 2.0_rp) -0.5_rp )*r_15 ! weighting between 1.0*0.5 and 1.0*2
       wts     = ( min(max( nscnres_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15 ! weighting between 1.4*0.5 and 1.4*2
       wtg     = ( min(max( nscnreg_s*r_14, 0.5_rp), 2.0_rp) -0.5_rp )*r_15 ! weighting between 1.4*0.5 and 1.4*2
       ! interpolation between two branches
       ventli  = (1.0_rp-wti)*ventli_s + wti*ventli_l
       ventls  = (1.0_rp-wts)*ventls_s + wts*ventls_l
       ventlg  = (1.0_rp-wtg)*ventlg_s + wtg*ventlg_l
       ! SB06(29)
       ssi         = qv/qsi(k) - 1.0_rp  ! supersaturation
       plidep_wrk  = gii*ssi*max(rhoq(k,i_ni),0.0_rp)*dq_xave(k,i_mp_qi)*ventli
       plsdep_wrk  = gis*ssi*max(rhoq(k,i_ns),0.0_rp)*dq_xave(k,i_mp_qs)*ventls
       plgdep_wrk  = gig*ssi*max(rhoq(k,i_ng),0.0_rp)*dq_xave(k,i_mp_qg)*ventlg
       !
       dep_limiter = rho(k)*(qv-qsi(k))/dt
       if      (ssi < -1.e-30_rp)then ! unsaturated
          plidep_total(k) = max(plidep_wrk+plsdep_wrk+plgdep_wrk, dep_limiter)
       else if (ssi >  1.e-30_rp)then !  saturated
          plidep_total(k) = min(plidep_wrk+plsdep_wrk+plgdep_wrk, dep_limiter)
       else
          plidep_total(k) = 0.0_rp
       end if
    enddo
 
    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 6319 of file scale_atmos_phy_mp_sn14.F90.

    implicit none
    !
    ! In this subroutine,
    ! We assumed surface temperature of droplets are same as environment.
 
    integer, intent(in) :: KA, KS, KE
 
    real(RP), intent(in) :: dt
    !
    real(RP), intent(in) :: tem(KA)
    !
    real(RP), intent(in) :: rhoq(KA,I_QV:I_NG)
    real(RP), intent(in) :: xq(KA,HYDRO_MAX)
    !
    real(RP), intent(inout):: PQ(KA,PQ_MAX)
    !
    real(RP), parameter :: temc_min = -65.0_rp
    real(RP), parameter :: a_het = 0.2_rp  ! SB06 (44)
    real(RP), parameter :: b_het = 0.65_rp ! SB06 (44)
    !
    real(RP) :: coef_m2_c
    real(RP) :: coef_m2_r
    ! temperature [celsius]
    real(RP) :: temc, temc2, temc3, temc4
    ! temperature function of homegenous/heterogenous freezing
    real(RP) :: Jhom, Jhet, Jh(KA)
    real(RP) :: rdt
    real(RP) :: tmp
    !
    integer :: k
    !
    rdt = 1.0_rp/dt
    !
    coef_m2_c =   coef_m2(i_mp_qc)
    coef_m2_r =   coef_m2(i_mp_qr)
    !
 
    ! Note, xc should be limited in range[xc_min:xc_max].
    ! and PNChom need to be calculated by NC
    ! because reduction rate of Nc need to be bound by NC.
    ! For the same reason PLChom also bound by LC and xc.
    ! Basically L and N should be independent
    ! but average particle mass x should be in suitable range.
 
    ! Homogenous Freezing
    do k = ks, ke
       pq(k,i_lchom) = 0.0_rp
       pq(k,i_nchom) = 0.0_rp
    end do
 
    ! Heterogenous Freezing
    do k = ks, ke
       temc = max( tem(k) - t00, temc_min )
       ! These cause from aerosol-droplet interaction.
       ! Bigg(1953) formula, Khain etal.(2000) eq.(4.5), Pruppacher and Klett(1997) eq.(9-48)
       jhet =  a_het*( exp( -b_het*temc ) - 1.0_rp )
       ! These cause in nature.
       ! Cotton and Field 2002, QJRMS. (12)
       if( temc < -30.0_rp ) then
          temc2 = temc*temc
          temc3 = temc*temc2
          temc4 = temc2*temc2
          jhom = 10.0_rp**(&
               - 243.40_rp - 14.75_rp*temc - 0.307_rp*temc2 &
               - 0.00287_rp*temc3 - 0.0000102_rp*temc4 ) *1.e+3_rp
       else if( temc < 0.0_rp) then
          jhom = 10._rp**(-7.63_rp-2.996_rp*(temc+30.0_rp))*1.e+3_rp
       else
          jhom = 0.0_rp
          jhet = 0.0_rp
       end if
       jh(k) = ( jhet + jhom ) * dt
    end do
 
    do k = ks, ke
#if defined(NVIDIA) || defined(SX)
       tmp = min( xq(k,i_mp_qc)*jh(k), 1.e+3_rp) ! apply exp limiter
       pq(k,i_lchet) = -rdt*rhoq(k,i_qc)*( 1.0_rp - exp( -coef_m2_c*tmp ) )
       pq(k,i_nchet) = -rdt*rhoq(k,i_nc)*( 1.0_rp - exp( -          tmp ) )
 
       tmp = min( xq(k,i_mp_qr)*jh(k), 1.e+3_rp) ! apply exp limiter
       pq(k,i_lrhet) = -rdt*rhoq(k,i_qr)*( 1.0_rp - exp( -coef_m2_r*tmp ) )
       pq(k,i_nrhet) = -rdt*rhoq(k,i_nr)*( 1.0_rp - exp( -          tmp ) )
#else
       pq(k,i_lchet) = -rdt*rhoq(k,i_qc)*( 1.0_rp - exp( -coef_m2_c*xq(k,i_mp_qc)*jh(k) ) )
       pq(k,i_nchet) = -rdt*rhoq(k,i_nc)*( 1.0_rp - exp( -          xq(k,i_mp_qc)*jh(k) ) )
       pq(k,i_lrhet) = -rdt*rhoq(k,i_qr)*( 1.0_rp - exp( -coef_m2_r*xq(k,i_mp_qr)*jh(k) ) )
       pq(k,i_nrhet) = -rdt*rhoq(k,i_nr)*( 1.0_rp - exp( -          xq(k,i_mp_qr)*jh(k) ) )
#endif
    end do
 
    !
    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 7138 of file scale_atmos_phy_mp_sn14.F90.

    implicit none
 
    integer,  intent(in)  :: KA, KS, KE
    integer,  intent(in)  :: QA_MP
    real(RP), intent(in)  :: QTRC0(KA,QA_MP)  ! tracer mass concentration [kg/kg]
    real(RP), intent(in)  :: DENS0(KA)        ! density                   [kg/m3]
    real(RP), intent(out) :: Crs(KA,HYDRO_MAX)! Cross section             [cm]
 
    ! mass concentration[kg/m3] and mean particle mass[kg]
    real(RP) :: xc(KA)
    real(RP) :: xr(KA)
    real(RP) :: xi(KA)
    real(RP) :: xs(KA)
    real(RP) :: xg(KA)
    ! radius of average mass
    real(RP) :: rc, rr
 
    real(RP) :: coef_Fuetal1998
    ! r2m_min is minimum value(moment of 1 particle with 1 micron)
    real(RP), parameter :: r2m_min=1.e-12_rp
    real(RP), parameter :: um2cm = 100.0_rp
 
    real(RP) :: limitsw, zerosw
    integer :: k
    !---------------------------------------------------------------------------
 
    ! mean particle mass[kg]
    do k  = ks, ke
       xc(k) = min( xc_max, max( xc_min, dens0(k)*qtrc0(k,i_qc)/(qtrc0(k,i_nc)+nc_min) ) )
       xr(k) = min( xr_max, max( xr_min, dens0(k)*qtrc0(k,i_qr)/(qtrc0(k,i_nr)+nr_min) ) )
       xi(k) = min( xi_max, max( xi_min, dens0(k)*qtrc0(k,i_qi)/(qtrc0(k,i_ni)+ni_min) ) )
       xs(k) = min( xs_max, max( xs_min, dens0(k)*qtrc0(k,i_qs)/(qtrc0(k,i_ns)+ns_min) ) )
       xg(k) = min( xg_max, max( xg_min, dens0(k)*qtrc0(k,i_qg)/(qtrc0(k,i_ng)+ng_min) ) )
    enddo
 
 
    do k = ks, ke
       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)
    enddo
 
    do k = ks, ke
       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)
    enddo
 
    do k = ks, ke
       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)
       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)
       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)
    enddo
 
    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.

  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.

  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.

  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.

  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.

  character(len=H_SHORT), parameter, public :: ATMOS_PHY_MP_SN14_tracer_names(QA_MP) = (/ &
       'QV', &
       'QC', &
       'QR', &
       'QI', &
       'QS', &
       'QG', &
       'NC', &
       'NR', &
       'NI', &
       'NS', &
       '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.

  character(len=H_MID)  , parameter, public :: ATMOS_PHY_MP_SN14_tracer_descriptions(QA_MP) = (/ &
       '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                                     '/)

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.

  character(len=H_SHORT), parameter, public :: ATMOS_PHY_MP_SN14_tracer_units(QA_MP) = (/ &
       'kg/kg ', &
       'kg/kg ', &
       'kg/kg ', &
       'kg/kg ', &
       'kg/kg ', &
       'kg/kg ', &
       'num/kg', &
       'num/kg', &
       'num/kg', &
       'num/kg', &
       'num/kg'  /)

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().