scale_atmos_phy_mp_sn14 Module Reference

module ATMOSPHERE / Physics Cloud Microphysics More...

Functions/Subroutines
subroutine, public	atmos_phy_mp_sn14_setup (KA, IA, JA)
	ATMOS_PHY_MP_sn14_setup setup. More...
subroutine, public	atmos_phy_mp_sn14_tendency (KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, W, QTRC, PRES, TEMP, Qdry, CPtot, CVtot, CCN, dt, cz, fz, RHOQ_t, RHOE_t, CPtot_t, CVtot_t, EVAPORATE)
	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	debug_tem_kij (KA, KS, KE, IA, IS, IE, JA, JS, JE, point, tem, rho, pre, qv)
subroutine	nucleation_kij (KA, KS, KE, IA, IS, IE, JA, JS, JE, cz, fz, w, rho, tem, pre, qdry, rhoq, cpa, dTdt_rad, qke, CCN, dt, PQ)
subroutine	ice_multiplication_kij (KA, KS, KE, IA, IS, IE, JA, JS, JE, Pac, tem, rhoq, xq, PQ)
subroutine	mixed_phase_collection_kij (KA, KS, KE, IA, IS, IE, JA, JS, JE, wtem, rhoq, xq, dq_xave, vt_xave, PQ, Pac)
subroutine	aut_acc_slc_brk_kij (KA, KS, KE, IA, IS, IE, JA, JS, JE, rhoq, xq, dq_xave, rho, PQ)
subroutine	freezing_water_kij (KA, KS, KE, IA, IS, IE, JA, JS, JE, dt, rhoq, xq, tem, PQ)
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	update_by_phase_change_kij (KA, KS, KE, IA, IS, IE, JA, JS, JE, ntdiv, ntmax, dt, cz, fz, w, dTdt_rad, rho, qdry, esw, esi, rhoq2, pre, tem, cpa, cva, PQ, sl_PLCdep, sl_PLRdep, sl_PNRdep, RHOQ_t, RHOE_t, CPtot_t, CVtot_t, qc_evaporate)

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

Detailed Description

module ATMOSPHERE / Physics Cloud Microphysics

Description

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

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

Author

Team SCALE

NAMELIST

• PARAM_ATMOS_PHY_MP_SN14

  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_TEM	logical	.false.
  OPT_DEBUG_INC	logical	.true.
  OPT_DEBUG_ACT	logical	.true.
  OPT_DEBUG_REE	logical	.true.
  OPT_DEBUG_BCS	logical	.true.
  NTMAX_PHASE_CHANGE	integer	1
  NTMAX_COLLECTION	integer	1

  
  
• 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)	0.462_RP
  NM_M92	real(RP)	1.E+3_RP
  AM_M92	real(RP)	-0.639_RP
  BM_M92	real(RP)	12.96_RP
  XC_CCN	real(RP)	1.E-12_RP	[kg]
  XI_CCN	real(RP)	1.E-12_RP	[kg] ! [move] 11/08/30 T.Mitsui
  TEM_CCN_LOW	real(RP)	233.150_RP	= -40 degC ! [Add] 10/08/03 T.Mitsui
  TEM_IN_LOW	real(RP)	173.150_RP	= -100 degC ! [Add] 10/08/03 T.Mitsui
  SSW_MAX	real(RP)	1.1_RP	[%]
  SSI_MAX	real(RP)	0.60_RP
  NUCL_TWOMEY	logical	.false.
  INUCL_W	logical	.false.

  
  
• PARAM_ATMOS_PHY_MP_SN14_collection

  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_condensation

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

  
  

History Output

No history output

Function/Subroutine Documentation

atmos_phy_mp_sn14_setup()

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

ATMOS_PHY_MP_sn14_setup setup.

Definition at line 523 of file scale_atmos_phy_mp_sn14.F90.

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

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_setup().

    use scale_prc, only: &
       prc_abort
    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 :: 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)

    wlabel(1) = "CLOUD"
    wlabel(2) = "RAIN"
    wlabel(3) = "ICE"
    wlabel(4) = "SNOW"
    wlabel(5) = "GRAUPEL"

    call mp_sn14_init

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

    return

Here is the call graph for this function:

Here is the caller graph for this function:

atmos_phy_mp_sn14_tendency()

subroutine, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_tendency	(	integer, intent(in)	KA,
		integer, intent(in)	KS,
		integer, intent(in)	KE,
		integer, intent(in)	IA,
		integer, intent(in)	IS,
		integer, intent(in)	IE,
		integer, intent(in)	JA,
		integer, intent(in)	JS,
		integer, intent(in)	JE,
		real(rp), dimension (ka,ia,ja), intent(in)	DENS,
		real(rp), dimension (ka,ia,ja), intent(in)	W,
		real(rp), dimension (ka,ia,ja,qa_mp), intent(in)	QTRC,
		real(rp), dimension(ka,ia,ja), intent(in)	PRES,
		real(rp), dimension(ka,ia,ja), intent(in)	TEMP,
		real(rp), dimension(ka,ia,ja), intent(in)	Qdry,
		real(rp), dimension(ka,ia,ja), intent(in)	CPtot,
		real(rp), dimension(ka,ia,ja), intent(in)	CVtot,
		real(rp), dimension (ka,ia,ja), intent(in)	CCN,
		real(dp), intent(in)	dt,
		real(rp), dimension( ka,ia,ja), intent(in)	cz,
		real(rp), dimension(0:ka,ia,ja), intent(in)	fz,
		real(rp), dimension (ka,ia,ja,qa_mp), intent(out)	RHOQ_t,
		real(rp), dimension (ka,ia,ja), intent(out)	RHOE_t,
		real(rp), dimension(ka,ia,ja), intent(out)	CPtot_t,
		real(rp), dimension(ka,ia,ja), intent(out)	CVtot_t,
		real(rp), dimension(ka,ia,ja), intent(out)	EVAPORATE
	)

ATMOS_PHY_MP_sn14_tendency calculate tendency.

Definition at line 591 of file scale_atmos_phy_mp_sn14.F90.

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_calc_tendency().

    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]

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

    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)

#ifdef PROFILE_FIPP
    call fipp_stop()
#endif

    return

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

References qa_mp.

Referenced by mod_atmos_phy_mp_vars::atmos_phy_mp_vars_get_diagnostic().

    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
    !---------------------------------------------------------------------------

    do j  = js, je
    do i  = is, ie
    do k  = ks, ke
       qhydro = 0.0_rp
       do iq = 1, qa_mp-1
          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

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

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

Referenced by mod_atmos_phy_mp_vars::atmos_phy_mp_vars_get_diagnostic().

    use scale_atmos_hydrometeor, only: &
       n_hyd, &
       i_hc,  &
       i_hr,  &
       i_hi,  &
       i_hs,  &
       i_hg
    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,ia,ja)
    real(RP) :: xr(ka,ia,ja)
    real(RP) :: xi(ka,ia,ja)
    real(RP) :: xs(ka,ia,ja)
    real(RP) :: xg(ka,ia,ja)
    ! diameter of average mass[kg/m3]
    real(RP) :: dc_ave(ka,ia,ja)
    real(RP) :: dr_ave(ka,ia,ja)
    ! radius of average mass
    real(RP) :: rc, rr
    ! 2nd. and 3rd. order moment of DSD
    real(RP) :: ri2m(ka,ia,ja), ri3m(ka,ia,ja)
    real(RP) :: rs2m(ka,ia,ja), rs3m(ka,ia,ja)
    real(RP) :: rg2m(ka,ia,ja), rg3m(ka,ia,ja)

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

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

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

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

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

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

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

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

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

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

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

    re(:,:,:,i_hg+1:) = 0.0_rp

    return

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

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

Referenced by mod_atmos_phy_mp_vars::atmos_phy_mp_vars_get_diagnostic().

    use scale_atmos_hydrometeor, only: &
       n_hyd, &
       i_hc,  &
       i_hr,  &
       i_hi,  &
       i_hs,  &
       i_hg
    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]
    !---------------------------------------------------------------------------

!OCL XFILL
    qe(:,:,:,i_hc) = qtrc0(:,:,:,i_mp_qc)
!OCL XFILL
    qe(:,:,:,i_hr) = qtrc0(:,:,:,i_mp_qr)
!OCL XFILL
    qe(:,:,:,i_hi) = qtrc0(:,:,:,i_mp_qi)
!OCL XFILL
    qe(:,:,:,i_hs) = qtrc0(:,:,:,i_mp_qs)
!OCL XFILL
    qe(:,:,:,i_hg) = qtrc0(:,:,:,i_mp_qg)
!OCL XFILL
    qe(:,:,:,i_hg+1:) = 0.0_rp

    return

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

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

Referenced by mod_atmos_phy_mp_vars::atmos_phy_mp_vars_get_diagnostic().

    use scale_atmos_hydrometeor, only: &
       n_hyd, &
       i_hc,  &
       i_hr,  &
       i_hi,  &
       i_hs,  &
       i_hg
    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]
    !---------------------------------------------------------------------------

!OCL XFILL
    ne(:,:,:,i_hc) = qtrc0(:,:,:,i_mp_nc)
!OCL XFILL
    ne(:,:,:,i_hr) = qtrc0(:,:,:,i_mp_nr)
!OCL XFILL
    ne(:,:,:,i_hi) = qtrc0(:,:,:,i_mp_ni)
!OCL XFILL
    ne(:,:,:,i_hs) = qtrc0(:,:,:,i_mp_ns)
!OCL XFILL
    ne(:,:,:,i_hg) = qtrc0(:,:,:,i_mp_ng)
!OCL XFILL
    ne(:,:,:,i_hg+1:) = 0.0_rp

    return

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

References aut_acc_slc_brk_kij(), scale_const::const_pi, 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, scale_atmos_hydrometeor::cv_water, debug_tem_kij(), freezing_water_kij(), 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, ice_multiplication_kij(), scale_io::io_fid_conf, mixed_phase_collection_kij(), scale_atmos_hydrometeor::n_hyd, nucleation_kij(), scale_prc::prc_abort(), scale_prof::prof_rapend(), scale_prof::prof_rapstart(), qa_mp, scale_specfunc::sf_gamma(), and update_by_phase_change_kij().

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_qhyd2qtrc().

    use scale_const, only: &
       pi => const_pi
    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 :: rhow =  1000.0_rp ! typical density for warm particles   [kg/m3]
    real(RP), parameter :: rhof =   100.0_rp ! typical density for frozen particles [kg/m3]
    real(RP), parameter :: rhog =   400.0_rp ! typical density for grapel particles [kg/m3]
    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

    if ( present(qnum) ) then

!OCL XFILL
       !$omp parallel do
       do j = js, je
       do i = is, ie
       do k = ks, ke
          qtrc(k,i,j,i_mp_nc) = qnum(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_nr) = qnum(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_ni) = qnum(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_ns) = qnum(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_ng) = qnum(k,i,j,i_hg) + qnum(k,i,j,i_hh)
       end do
       end do
       end do

    else
       piov6 = pi / 6.0_rp

!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

Here is the call graph for this function:

Here is the caller graph for this function:

debug_tem_kij()

subroutine scale_atmos_phy_mp_sn14::debug_tem_kij	(	integer, intent(in)	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,
		integer, intent(in)	point,
		real(rp), dimension(ka,ia,ja), intent(in)	tem,
		real(rp), dimension(ka,ia,ja), intent(in)	rho,
		real(rp), dimension(ka,ia,ja), intent(in)	pre,
		real(rp), dimension (ka,ia,ja), intent(in)	qv
	)

Definition at line 2587 of file scale_atmos_phy_mp_sn14.F90.

References scale_prc::prc_myrank.

Referenced by atmos_phy_mp_sn14_qhyd2qtrc().

    use scale_prc, only: &
       prc_myrank
    implicit none

    integer, intent(in) :: ka, ks, ke
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je

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

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

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

          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,i,j), rho(k,i,j), pre(k,i,j), qv(k,i,j), k, i, j, prc_myrank
       endif
    enddo
    enddo
    enddo

    return

Here is the caller graph for this function:

nucleation_kij()

subroutine scale_atmos_phy_mp_sn14::nucleation_kij	(	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)	cz,
		real(rp), dimension(0:ka,ia,ja), intent(in)	fz,
		real(rp), dimension (ka,ia,ja), intent(in)	w,
		real(rp), dimension(ka,ia,ja), intent(in)	rho,
		real(rp), dimension(ka,ia,ja), intent(in)	tem,
		real(rp), dimension(ka,ia,ja), intent(in)	pre,
		real(rp), dimension(ka,ia,ja), intent(in)	qdry,
		real(rp), dimension(i_qv:i_ng,ka,ia,ja), intent(in)	rhoq,
		real(rp), dimension(ka,ia,ja), intent(in)	cpa,
		real(rp), dimension(ka,ia,ja), intent(in)	dTdt_rad,
		real(rp), dimension(ka,ia,ja), intent(in)	qke,
		real(rp), dimension(ka,ia,ja), intent(in)	CCN,
		real(rp), intent(in)	dt,
		real(rp), dimension(pq_max,ka,ia,ja), intent(out)	PQ
	)

Definition at line 2632 of file scale_atmos_phy_mp_sn14.F90.

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

Referenced by atmos_phy_mp_sn14_qhyd2qtrc().

    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
    implicit none

    integer, intent(in) :: ka, ks, ke
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je

    real(RP), intent(in)  :: cz(  ka,ia,ja)   !
    real(RP), intent(in)  :: fz(0:ka,ia,ja)   !
    real(RP), intent(in)  :: w  (ka,ia,ja)    ! w of full level
    real(RP), intent(in)  :: rho(ka,ia,ja)    ! [Add] 09/08/18 T.Mitsui
    real(RP), intent(in)  :: tem(ka,ia,ja)    ! [Add] 09/08/18 T.Mitsui
    real(RP), intent(in)  :: pre(ka,ia,ja)    ! [Add] 09/08/18 T.Mitsui
    real(RP), intent(in)  :: qdry(ka,ia,ja)
    !
    real(RP), intent(in)  :: rhoq(i_qv:i_ng,ka,ia,ja)     !
    !
    real(RP), intent(in) ::  cpa(ka,ia,ja)      ! in  09/08/18 [Add] T.Mitsui
    real(RP), intent(in)  :: dtdt_rad(ka,ia,ja) ! 09/08/18 T.Mitsui
    real(RP), intent(in)  :: qke(ka,ia,ja)      ! 09/08/18 T.Mitsui
    real(RP), intent(in)  :: dt
    real(RP), intent(in)  :: ccn(ka,ia,ja)
    !
    real(RP), intent(out) :: pq(pq_max,ka,ia,ja)
    !
    ! namelist variables
    !
    ! total aerosol number concentration [/m3]
    real(RP), parameter :: c_ccn_ocean= 1.00e+8_rp
    real(RP), parameter :: c_ccn_land = 1.26e+9_rp
    real(RP), save      :: c_ccn      = 1.00e+8_rp
    ! aerosol activation factor
    real(RP), parameter :: kappa_ocean= 0.462_rp
    real(RP), parameter :: kappa_land = 0.308_rp
    real(RP), save      :: kappa      = 0.462_rp
    real(RP), save      :: c_in       = 1.0_rp
    ! SB06 (36)
    real(RP), save :: nm_m92 = 1.e+3_rp
    real(RP), save :: am_m92 = -0.639_rp
    real(RP), save :: bm_m92 = 12.96_rp
    !
    real(RP), save :: in_max = 1000.e+3_rp ! max num. of Ice-Nuclei [num/m3]
    real(RP), save :: ssi_max= 0.60_rp
    real(RP), save :: ssw_max= 1.1_rp  ! [%]
    !
    logical, save :: flag_first = .true.
    real(RP), save :: qke_min = 0.03_rp ! sigma=0.1[m/s], 09/08/18 T.Mitsui
    real(RP), save :: tem_ccn_low=233.150_rp  ! = -40 degC  ! [Add] 10/08/03 T.Mitsui
    real(RP), save :: tem_in_low =173.150_rp  ! = -100 degC ! [Add] 10/08/03 T.Mitsui
    logical, save :: nucl_twomey = .false.
    logical, save :: inucl_w     = .false.
    !
    namelist / param_atmos_phy_mp_sn14_nucleation / &
         in_max,                     & !
         c_ccn, kappa,               & ! cloud nucleation
         nm_m92, am_m92, bm_m92,     & ! ice nucleation
         xc_ccn, xi_ccn,             &
         tem_ccn_low,                & ! [Add] 10/08/03 T.Mitsui
         tem_in_low,                 & ! [Add] 10/08/03 T.Mitsui
         ssw_max, ssi_max,           &
         nucl_twomey, inucl_w        ! [Add] 13/01/30 Y.Sato
    !
!    real(RP) :: c_ccn_map(1,IA,JA)   ! c_ccn horizontal distribution
!    real(RP) :: kappa_map(1,IA,JA)   ! kappa horizontal distribution
!    real(RP) :: c_in_map(1,IA,JA)    ! c_in  horizontal distribution ! [Add] 11/08/30 T.Mitsui
    real(RP) :: esw(ka,ia,ja)      ! saturation vapor pressure, water
    real(RP) :: esi(ka,ia,ja)      !                            ice
    real(RP) :: ssw(ka,ia,ja)      ! super saturation (water)
    real(RP) :: ssi(ka,ia,ja)      ! super saturation (ice)
!    real(RP) :: w_dsswdz(KA,IA,JA) ! w*(d_ssw/ d_z) super saturation(water) flux
    real(RP) :: w_dssidz(ka,ia,ja) ! w*(d_ssi/ d_z), 09/04/14 T.Mitsui
!    real(RP) :: ssw_below(KA,IA,JA)! ssw(k-1)
    real(RP) :: ssi_below(ka,ia,ja)! ssi(k-1), 09/04/14 T.Mitsui
    real(RP) :: z_below(ka,ia,ja)  ! z(k-1)
    real(RP) :: dzh                ! z(k)-z(k-1)
    real(RP) :: pv                 ! vapor pressure
    ! work variables for Twomey Equation.
    real(RP) :: qsw(ka,ia,ja)
    real(RP) :: qsi(ka,ia,ja)
    real(RP) :: dqsidtem_rho(ka,ia,ja)
    real(RP) :: dssidt_rad(ka,ia,ja)
!    real(RP) :: dni_ratio(KA,IA,JA)
    real(RP) :: wssi, wdssi
    !
!    real(RP) :: xi_nuc(1,IA,JA)    ! xi use the value @ cloud base
!    real(RP) :: alpha_nuc(1,IA,JA) ! alpha_nuc
!    real(RP) :: eta_nuc(1,IA,JA)   ! xi use the value @ cloud base
    !
    real(RP) :: sigma_w(ka,ia,ja)
    real(RP) :: weff(ka,ia,ja)
    real(RP) :: weff_max(ka,ia,ja)
    real(RP) :: velz(ka)
    !
    real(RP) :: coef_ccn(ia,ja)
    real(RP) :: slope_ccn(ia,ja)
    real(RP) :: nc_new(ka,ia,ja)
    real(RP) :: nc_new_below(ka,ia,ja)
    real(RP) :: dnc_new
    real(RP) :: nc_new_max   ! Lohmann (2002),
    real(RP) :: a_max
    real(RP) :: b_max
    logical :: flag_nucleation(ka,ia,ja)
    !
    real(RP) :: r_gravity
    real(RP), parameter :: r_sqrt3=0.577350269_rp ! = sqrt(1.d0/3.d0)
    real(RP), parameter :: eps=1.e-30_rp
    !====> ! 09/08/18
    !
    real(RP) :: dlcdt_max, dli_max ! defined by supersaturation
    real(RP) :: dncdt_max, dni_max ! defined by supersaturation
    real(RP) :: rdt
    !
    integer :: i, j, k
    !
    if( flag_first )then
       rewind(io_fid_conf)
       read(io_fid_conf, nml=param_atmos_phy_mp_sn14_nucleation, end=100)
100    log_nml(param_atmos_phy_mp_sn14_nucleation)
       flag_first=.false.

       if ( mp_couple_aerosol .AND. nucl_twomey ) then
          log_error("ATMOS_PHY_MP_SN14_nucleation_kij",*) "nucl_twomey should be false when MP_couple_aerosol is true, stop"
          call prc_abort
       endif
    endif
    !
!    c_ccn_map(1,:,:) = c_ccn
!    kappa_map(1,:,:) = kappa
!    c_in_map(1,:,:)  = c_in
    !
!    nc_uplim_d(1,:,:)  = c_ccn_map(1,:,:)*1.5_RP
    do j = js, je
    do i = is, ie
       nc_uplim_d(1,i,j)  = c_ccn*1.5_rp
    end do
    end do
    !
    rdt            = 1.0_rp/dt
    r_gravity      = 1.0_rp/grav
    !
    call moist_psat_liq      ( ka, ks, ke, ia, is, ie, ja, js, je, &
                              tem(:,:,:), esw(:,:,:) ) ! [IN], [OUT]
    call moist_psat_ice      ( ka, ks, ke, ia, is, ie, ja, js, je, &
                              tem(:,:,:), esi(:,:,:) ) ! [IN], [OUT]
    call moist_pres2qsat_liq ( ka, ks, ke, ia, is, ie, ja, js, je, &
                              tem(:,:,:), pre(:,:,:), qdry(:,:,:), & ! [IN]
                              qsw(:,:,:)                           ) ! [OUT]
    call moist_pres2qsat_ice ( ka, ks, ke, ia, is, ie, ja, js, je, &
                              tem(:,:,:), pre(:,:,:), qdry(:,:,:), & ! [IN]
                              qsi(:,:,:)                           ) ! [OUT]
    call moist_dqsi_dtem_dens( ka, ks, ke, ia, is, ie, ja, js, je, &
                               tem(:,:,:), rho(:,:,:), & ! [IN]
                               dqsidtem_rho(:,:,:)     ) ! [OUT]
    !
    ! Lohmann (2002),JAS, eq.(1) but changing unit [cm-3] => [m-3]
    a_max = 1.e+6_rp*0.1_rp*(1.e-6_rp)**1.27_rp
    b_max = 1.27_rp
    !
    ssi_max = 1.0_rp

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

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

    end do
    end do
    !
    if( mp_couple_aerosol ) then

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

    else

      if( nucl_twomey ) then
        ! diagnose cloud condensation nuclei

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

    endif

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

    end do
    end do

    if( mp_couple_aerosol ) then

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

    else

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

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

    return

Here is the call graph for this function:

Here is the caller graph for this function:

ice_multiplication_kij()

subroutine scale_atmos_phy_mp_sn14::ice_multiplication_kij	(	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(pac_max,ka,ia,ja), intent(in)	Pac,
		real(rp), dimension(ka,ia,ja), intent(in)	tem,
		real(rp), dimension(i_qv:i_ng,ka,ia,ja), intent(in)	rhoq,
		real(rp), dimension(hydro_max,ka,ia,ja), intent(in)	xq,
		real(rp), dimension(pq_max,ka,ia,ja), intent(inout)	PQ
	)

Definition at line 3037 of file scale_atmos_phy_mp_sn14.F90.

References scale_specfunc::sf_gamma().

Referenced by atmos_phy_mp_sn14_qhyd2qtrc().

    ! ice multiplication by splintering
    ! we consider Hallet-Mossop process
    use scale_specfunc, only: &
       gammafunc => sf_gamma
    implicit none

    integer, intent(in) :: ka, ks, ke
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je
    !
    real(RP), intent(in) :: pac(pac_max,ka,ia,ja)
    real(RP), intent(in) :: tem(ka,ia,ja)
    real(RP), intent(in) :: rhoq(i_qv:i_ng,ka,ia,ja)
    real(RP), intent(in) :: xq(hydro_max,ka,ia,ja)
    !
    real(RP), intent(inout):: pq(pq_max,ka,ia,ja)
    !
    logical, save :: flag_first      = .true.
    ! production of (350.d3 ice particles)/(cloud riming [g]) => 350*d6 [/kg]
    real(RP), parameter :: pice = 350.0e+6_rp
    ! production of (1 ice particle)/(250 cloud particles riming)
    real(RP), parameter :: pnc  = 250.0_rp
    real(RP), parameter :: rc_cr= 12.e-6_rp ! critical size[micron]
    ! temperature factor
    real(RP) :: fp
    ! work for incomplete gamma function
    real(RP), save :: xc_cr    ! mass[kg] of cloud with r=critical size[micron]
    real(RP), save :: alpha    ! slope parameter of gamma function
    real(RP), save :: gm, lgm  ! gamma(alpha), log(gamma(alpha))
    real(RP) :: igm            ! in complete gamma(x,alpha)
    real(RP) :: x
    ! coefficient of expansion using in calculation of igm
    real(RP) :: a0,a1,a2,a3,a4,a5
    real(RP) :: a6,a7,a8,a9,a10
    real(RP) :: an1,an2,b0,b1,b2,c0,c1,c2
    real(RP) :: d0,d1,d2,e1,e2,h0,h1,h2
    real(RP), parameter :: eps=1.0e-30_rp
    ! number of cloud droplets larger than 12 micron(radius).
    real(RP) :: n12
    !
    real(RP) :: wn, wni, wns, wng
    integer :: i, j, k
    !
    if( flag_first )then
       flag_first = .false.
       ! work for Incomplete Gamma function
       xc_cr = (2.0_rp*rc_cr/a_m(i_mp_qc))**(1.0_rp/b_m(i_mp_qc))
       alpha = (nu(i_mp_qc)+1.0_rp)/mu(i_mp_qc)
       gm    = gammafunc(alpha)
       lgm   = log(gm)
    end if
    !
    do j = js, je
    do i = is, ie
    do k = ks, ke
       ! Here we assume particle temperature is same as environment temperature.
       ! If you want to treat in a better manner,
       ! you can diagnose with eq.(64) in CT(86)
       if     (tem(k,i,j) >  270.16_rp)then
          fp   = 0.0_rp
       else if(tem(k,i,j) >= 268.16_rp)then
          fp   = (270.16_rp-tem(k,i,j))*0.5_rp
       else if(tem(k,i,j) >= 265.16_rp)then
          fp   = (tem(k,i,j)-265.16_rp)*0.333333333_rp
       else
          fp   = 0.0_rp
       end if
       ! Approximation of Incomplete Gamma function
       ! Here we calculate with algorithm by Numerical Recipes.
       ! This approach is based on eq.(78) in Cotton etal.(1986),
       ! but more accurate and expanded for Generalized Gamma Distribution.
       x       = coef_lambda(i_mp_qc)*(xc_cr/xq(i_mp_qc,k,i,j))**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*1.d2) then ! continued fraction expansion
          ! 2nd-truncation is enough for cloud droplet.
          ! setup
          b0   = x+1.0_rp-alpha
          c0   = 1.0_rp/eps
          d0   = 1.0_rp/b0
          h0   = d0
          ! n=1
          an1  = -(1.0_rp-alpha)
          b1   = b0 + 2.0_rp
          d1   = 1.0_rp/(an1*d0+b1)
          c1   = b1+an1/c0
          e1   = d1*c1
          h1   = h0*e1
          ! n=2
          an2  = -2.0_rp*(2.0_rp-alpha)
          b2   = b1 + 2.0_rp
          d2   = 1.0_rp/(an2*d1+b2)
          c2   = b2+an2/c1
          e2   = d2*c2
          h2   = h1*e2
          !
          igm  = 1.0_rp - exp( -x + alpha*log(x) - lgm )*h2
       else                       ! negligible
          igm  = 1.0_rp
       end if
       ! n12 is number of cloud droplets larger than 12 micron.
       n12 = rhoq(i_nc,k,i,j)*(1.0_rp-igm)
       ! eq.(82) CT(86)
       wn           = (pice + n12/((rhoq(i_qc,k,i,j)+xc_min)*pnc) )*fp ! filtered by xc_min
       wni          = wn*(-pac(i_liaclc2li,k,i,j)  ) ! riming production rate is all negative
       wns          = wn*(-pac(i_lsaclc2ls,k,i,j)  )
       wng          = wn*(-pac(i_lgaclc2lg,k,i,j)  )
       pq(i_nispl,k,i,j) = wni+wns+wng
       !
       pq(i_lsspl,k,i,j) = - wns*xq(i_mp_qi,k,i,j) ! snow    => ice
       pq(i_lgspl,k,i,j) = - wng*xq(i_mp_qi,k,i,j) ! graupel => ice
       !
    end do
    end do
    end do
    !
    return

Here is the call graph for this function:

Here is the caller graph for this function:

mixed_phase_collection_kij()

subroutine scale_atmos_phy_mp_sn14::mixed_phase_collection_kij	(	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)	wtem,
		real(rp), dimension(i_qv:i_ng,ka,ia,ja), intent(in)	rhoq,
		real(rp), dimension(hydro_max,ka,ia,ja), intent(in)	xq,
		real(rp), dimension(hydro_max,ka,ia,ja), intent(in)	dq_xave,
		real(rp), dimension(hydro_max,2,ka,ia,ja), intent(in)	vt_xave,
		real(rp), dimension(pq_max,ka,ia,ja), intent(inout)	PQ,
		real(rp), dimension(pac_max,ka,ia,ja), intent(out)	Pac
	)

Definition at line 3180 of file scale_atmos_phy_mp_sn14.F90.

References scale_io::io_fid_conf.

Referenced by atmos_phy_mp_sn14_qhyd2qtrc().

    use scale_atmos_saturation, only: &
       moist_psat_ice => atmos_saturation_psat_ice
    implicit none

    integer, intent(in) :: ka, ks, ke
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je

    !--- mixed-phase collection process
    !                  And all we set all production term as a negative sign to avoid confusion.
    !
    real(RP), intent(in) :: wtem(ka,ia,ja)
    !--- mass/number concentration[kg/m3]
    real(RP), intent(in) :: rhoq(i_qv:i_ng,ka,ia,ja)
    ! necessary ?
    real(RP), intent(in) :: xq(hydro_max,ka,ia,ja)
    !--- diameter of averaged mass( D(ave x) )
    real(RP), intent(in) :: dq_xave(hydro_max,ka,ia,ja)
    !--- terminal velocity of averaged mass( vt(ave x) )
    real(RP), intent(in) :: vt_xave(hydro_max,2,ka,ia,ja)
    ! [Add] 11/08/30 T.Mitsui, for autoconversion of ice
    !    real(RP), intent(in) :: rho(KA,IA,JA)
    !--- partial conversion
    real(RP), intent(inout):: pq(pq_max,ka,ia,ja)
    !
    real(RP), intent(out):: pac(pac_max,ka,ia,ja)
    !
    ! namelist variables
    !=== for collection
    !--- threshold of diameters to collide with others
    real(RP), save :: dc0 =  15.0e-6_rp       ! lower threshold of cloud
    real(RP), save :: dc1 =  40.0e-6_rp       ! upper threshold of cloud
    real(RP), save :: di0 = 150.0e-6_rp       ! lower threshold of cloud
    real(RP), save :: ds0 = 150.0e-6_rp       ! lower threshold of cloud
    real(RP), save :: dg0 = 150.0e-6_rp       ! lower threshold of cloud
    !--- standard deviation of terminal velocity[m/s]
    real(RP), save :: sigma_c=0.00_rp        ! cloud
    real(RP), save :: sigma_r=0.00_rp        ! rain
    real(RP), save :: sigma_i=0.2_rp        ! ice
    real(RP), save :: sigma_s=0.2_rp        ! snow
    real(RP), save :: sigma_g=0.00_rp        ! graupel
    !--- max collection efficiency for cloud
    real(RP), save :: e_im = 0.80_rp        ! ice max
    real(RP), save :: e_sm = 0.80_rp        ! snow max
    real(RP), save :: e_gm = 1.00_rp        ! graupel max
    !--- collection efficiency between 2 species
    real(RP), save :: e_ir=1.0_rp            ! ice     x rain
    real(RP), save :: e_sr=1.0_rp            ! snow    x rain
    real(RP), save :: e_gr=1.0_rp            ! graupel x rain
    real(RP), save :: e_ii=1.0_rp            ! ice     x ice
    real(RP), save :: e_si=1.0_rp            ! snow    x ice
    real(RP), save :: e_gi=1.0_rp            ! graupel x ice
    real(RP), save :: e_ss=1.0_rp            ! snow    x snow
    real(RP), save :: e_gs=1.0_rp            ! graupel x snow
    real(RP), save :: e_gg=1.0_rp            ! graupel x graupel
    !=== for partial conversion
    !--- flag: 1=> partial conversion to graupel, 0=> no conversion
    integer, save :: i_iconv2g=1          ! ice  => graupel
    integer, save :: i_sconv2g=1          ! snow => graupel
    !--- bulk density of graupel
    real(RP), save :: rho_g   = 900.0_rp     ! [kg/m3]
    !--- space filling coefficient [%]
    real(RP), save :: cfill_i = 0.68_rp     ! ice
    real(RP), save :: cfill_s = 0.01_rp     ! snow
    !--- critical diameter for ice conversion
    real(RP), save :: di_cri  = 500.e-6_rp    ! [m]
    ! [Add] 10/08/03 T.Mitsui
    logical, save :: opt_stick_ks96=.false.
    logical, save :: opt_stick_co86=.false.
    real(RP), parameter :: a_dec = 0.883_rp
    real(RP), parameter :: b_dec = 0.093_rp
    real(RP), parameter :: c_dec = 0.00348_rp
    real(RP), parameter :: d_dec = 4.5185e-5_rp
    !
    logical, save :: flag_first = .true.
    namelist / param_atmos_phy_mp_sn14_collection / &
         dc0, dc1, di0, ds0, dg0,    &
         sigma_c, sigma_r, sigma_i, sigma_s, sigma_g, &
         opt_stick_ks96,   &
         opt_stick_co86,   &
         e_im, e_sm, e_gm, &
         e_ir, e_sr, e_gr, e_ii, e_si, e_gi, e_ss, e_gs, e_gg, &
         i_iconv2g, i_sconv2g, rho_g, cfill_i, cfill_s, di_cri
    !
    real(RP) :: tem(ka,ia,ja)
    !
    !--- collection efficency of each specie
    real(RP) :: e_c, e_r, e_i, e_s, e_g    !
    real(RP) :: e_ic, e_sc, e_gc           !
    !--- sticking efficiency
    real(RP) :: e_stick(ka,ia,ja)
    ! [Add] 10/08/03 T.Mitsui
    real(RP) :: temc, temc2, temc3
    real(RP) :: e_dec
    real(RP) :: esi_rat
    real(RP) :: esi(ka,ia,ja)
    !
    real(RP) :: temc_p, temc_m             ! celcius tem.
    ! [Add] 11/08/30 T.Mitsui, estimation of autoconversion time
!    real(RP) :: ci_aut(KA,IA,JA)
!    real(RP) :: taui_aut(KA,IA,JA)
!    real(RP) :: tau_sce(KA,IA,JA)
    !--- DSD averaged diameter for each species
    real(RP) :: ave_dc                     ! cloud
!    real(RP) :: ave_dr                     ! rain
    real(RP) :: ave_di                     ! ice
    real(RP) :: ave_ds                     ! snow
    real(RP) :: ave_dg                     ! graupel
    !--- coefficient of collection equations(L:mass, N:number)
    real(RP) :: coef_acc_lci,   coef_acc_nci   ! cloud     - cloud ice
    real(RP) :: coef_acc_lcs,   coef_acc_ncs   ! cloud     - snow
    !
    real(RP) :: coef_acc_lcg,   coef_acc_ncg   ! cloud     - graupel
    real(RP) :: coef_acc_lri_i, coef_acc_nri_i ! rain      - cloud ice
    real(RP) :: coef_acc_lri_r, coef_acc_nri_r ! rain      - cloud ice
    real(RP) :: coef_acc_lrs_s, coef_acc_nrs_s ! rain      - snow
    real(RP) :: coef_acc_lrs_r, coef_acc_nrs_r ! rain      - snow
    real(RP) :: coef_acc_lrg,   coef_acc_nrg   ! rain      - graupel
    real(RP) :: coef_acc_lii,   coef_acc_nii   ! cloud ice - cloud ice
    real(RP) :: coef_acc_lis,   coef_acc_nis   ! cloud ice - snow
    real(RP) ::                 coef_acc_nss   ! snow      - snow
    real(RP) ::                 coef_acc_ngg   ! grauepl   - graupel
    real(RP) :: coef_acc_lsg,   coef_acc_nsg   ! snow      - graupel
    !--- (diameter) x (diameter)
    real(RP) :: dcdc, dcdi, dcds, dcdg
    real(RP) :: drdr, drdi, drds, drdg
    real(RP) :: didi, dids, didg
    real(RP) :: dsds, dsdg
    real(RP) :: dgdg
    !--- (terminal velocity) x (terminal velocity)
    real(RP) :: vcvc, vcvi, vcvs, vcvg
    real(RP) :: vrvr, vrvi, vrvs, vrvg
    real(RP) :: vivi, vivs, vivg
    real(RP) :: vsvs, vsvg
    real(RP) :: vgvg
    !
    real(RP) :: wx_cri, wx_crs
    real(RP) :: coef_emelt
    real(RP) :: w1

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

    call moist_psat_ice( ka, ks, ke, ia, is, ie, ja, js, je, &
                         tem(:,:,:), esi(:,:,:)  ) ! [IN], [OUT]

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

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

    !
    return

Here is the caller graph for this function:

aut_acc_slc_brk_kij()

subroutine scale_atmos_phy_mp_sn14::aut_acc_slc_brk_kij	(	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(i_qv:i_ng,ka,ia,ja), intent(in)	rhoq,
		real(rp), dimension(hydro_max,ka,ia,ja), intent(in)	xq,
		real(rp), dimension(hydro_max,ka,ia,ja), intent(in)	dq_xave,
		real(rp), dimension(ka,ia,ja), intent(in)	rho,
		real(rp), dimension(pq_max,ka,ia,ja), intent(inout)	PQ
	)

Definition at line 3673 of file scale_atmos_phy_mp_sn14.F90.

References scale_const::const_eps.

Referenced by atmos_phy_mp_sn14_qhyd2qtrc().

    implicit none

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

Here is the caller graph for this function:

freezing_water_kij()

subroutine scale_atmos_phy_mp_sn14::freezing_water_kij	(	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), intent(in)	dt,
		real(rp), dimension(i_qv:i_ng,ka,ia,ja), intent(in)	rhoq,
		real(rp), dimension(hydro_max,ka,ia,ja), intent(in)	xq,
		real(rp), dimension(ka,ia,ja), intent(in)	tem,
		real(rp), dimension(pq_max,ka,ia,ja), intent(inout)	PQ
	)

Definition at line 4016 of file scale_atmos_phy_mp_sn14.F90.

Referenced by atmos_phy_mp_sn14_qhyd2qtrc().

    implicit none
    !
    ! In this subroutine,
    ! We assumed surface temperature of droplets are same as environment.

    integer, intent(in) :: ka, ks, ke
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je

    real(RP), intent(in) :: dt
    !
    real(RP), intent(in) :: tem(ka,ia,ja)
    !
    real(RP), intent(in) :: rhoq(i_qv:i_ng,ka,ia,ja)
    real(RP), intent(in) :: xq(hydro_max,ka,ia,ja)
    !
    real(RP), intent(inout):: pq(pq_max,ka,ia,ja)
    !
    real(RP), parameter :: temc_min = -65.0_rp
    real(RP), parameter :: a_het = 0.2_rp  ! SB06 (44)
    real(RP), parameter :: b_het = 0.65_rp ! SB06 (44)
    !
    real(RP) :: coef_m2_c
    real(RP) :: coef_m2_r
    ! temperature [celsius]
    real(RP) :: temc, temc2, temc3, temc4
    ! temperature function of homegenous/heterogenous freezing
    real(RP) :: jhom, jhet
    real(RP) :: rdt
    real(RP) :: tmp
    !
    integer :: i,j,k
    !
    rdt = 1.0_rp/dt
    !
    coef_m2_c =   coef_m2(i_mp_qc)
    coef_m2_r =   coef_m2(i_mp_qr)
    !
    profile_start("sn14_freezing")
    do j = js, je
    do i = is, ie
       do k = ks, ke
          temc = max( tem(k,i,j) - t00, temc_min )
          ! These cause from aerosol-droplet interaction.
          ! Bigg(1953) formula, Khain etal.(2000) eq.(4.5), Pruppacher and Klett(1997) eq.(9-48)
          jhet =  a_het*exp( -b_het*temc - 1.0_rp )
          ! These cause in nature.
          ! Cotton and Field 2002, QJRMS. (12)
          if( temc < -65.0_rp )then
             jhom = 10.0_rp**(24.37236_rp)*1.e+3_rp
             jhet =  a_het*exp( 65.0_rp*b_het - 1.0_rp ) ! 09/04/14 [Add], fixer T.Mitsui
          else if( temc < -30.0_rp ) then
             temc2 = temc*temc
             temc3 = temc*temc2
             temc4 = temc2*temc2
             jhom = 10.0_rp**(&
                  - 243.40_rp - 14.75_rp*temc - 0.307_rp*temc2 &
                  - 0.00287_rp*temc3 - 0.0000102_rp*temc4 ) *1.e+3_rp
          else if( temc < 0.0_rp) then
             jhom = 10._rp**(-7.63_rp-2.996_rp*(temc+30.0_rp))*1.e+3_rp
          else
             jhom = 0.0_rp
             jhet = 0.0_rp
          end if
          ! Note, xc should be limited in range[xc_min:xc_max].
          ! and PNChom need to be calculated by NC
          ! because reduction rate of Nc need to be bound by NC.
          ! For the same reason PLChom also bound by LC and xc.
          ! Basically L and N should be independent
          ! but average particle mass x should be in suitable range.
          ! Homogenous Freezing
          pq(i_lchom,k,i,j) = 0.0_rp
          pq(i_nchom,k,i,j) = 0.0_rp
          ! Heterogenous Freezing
#if defined(PGI) || defined(SX)
          tmp = min( xq(i_mp_qc,k,i,j)*(jhet+jhom)*dt, 1.e+3_rp) ! apply exp limiter
          pq(i_lchet,k,i,j) = -rdt*rhoq(i_qc,k,i,j)*( 1.0_rp - exp( -coef_m2_c*tmp ) )
          pq(i_nchet,k,i,j) = -rdt*rhoq(i_nc,k,i,j)*( 1.0_rp - exp( -          tmp ) )

          tmp = min( xq(i_mp_qr,k,i,j)*(jhet+jhom)*dt, 1.e+3_rp) ! apply exp limiter
          pq(i_lrhet,k,i,j) = -rdt*rhoq(i_qr,k,i,j)*( 1.0_rp - exp( -coef_m2_r*tmp ) )
          pq(i_nrhet,k,i,j) = -rdt*rhoq(i_nr,k,i,j)*( 1.0_rp - exp( -          tmp ) )
#else
          pq(i_lchet,k,i,j) = -rdt*rhoq(i_qc,k,i,j)*( 1.0_rp - exp( -coef_m2_c*xq(i_mp_qc,k,i,j)*(jhet+jhom)*dt ) )
          pq(i_nchet,k,i,j) = -rdt*rhoq(i_nc,k,i,j)*( 1.0_rp - exp( -          xq(i_mp_qc,k,i,j)*(jhet+jhom)*dt ) )
          pq(i_lrhet,k,i,j) = -rdt*rhoq(i_qr,k,i,j)*( 1.0_rp - exp( -coef_m2_r*xq(i_mp_qr,k,i,j)*(jhet+jhom)*dt ) )
          pq(i_nrhet,k,i,j) = -rdt*rhoq(i_nr,k,i,j)*( 1.0_rp - exp( -          xq(i_mp_qr,k,i,j)*(jhet+jhom)*dt ) )
#endif
       end do
    end do
    end do
    profile_stop("sn14_freezing")
    !
    return

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

References scale_const::const_undef, and qa_mp.

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_calc_tendency().

    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       ! average mass of 1 particle( mass/number )

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

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


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

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

    profile_start("sn14_terminal_vel")

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


    do k = ks, ke
       rhofac = rho_0 / max( dens(k), rho_min )

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

       vterm(k,i_mp_qc) = -rhofac_q(i_mp_qc) * coef_vt1(i_mp_qc,1) * xq**beta_v(i_mp_qc,1)
       ! NC
       vterm(k,i_mp_nc) = -rhofac_q(i_mp_qc) * coef_vt0(i_mp_qc,1) * xq**beta_vn(i_mp_qc,1)

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

       rlambdar = 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)
       dq = ( 4.0_rp + mud_r ) * rlambdar ! D^(3)+mu weighted mean diameter
       dql = dq
       weight = min( 1.0_rp, max( 0.0_rp, 0.5_rp * ( 1.0_rp + tanh( pi * log( dq/d_vtr_branch ) ) ) ) )

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

       velq_s = coef_vtr_ar2 * dql &
            * ( 1.0_rp - ( 1.0_rp + coef_vtr_br2*rlambdar )**(-2.0_rp-mud_r) )
       velq_l = coef_vtr_ar1 - coef_vtr_br1 &
            * ( 1.0_rp + coef_vtr_cr1*rlambdar )**(-1.0_rp-mud_r)
       vterm(k,i_mp_nr) = -rhofac_q(i_mp_qr) &
            * ( velq_l * (          weight ) &
              + velq_s * ( 1.0_rp - weight ) )

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

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

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

       velq_s = coef_vt0(i_mp_qi,1) * xq**beta_vn(i_mp_qi,1)
       velq_l = coef_vt0(i_mp_qi,2) * xq**beta_vn(i_mp_qi,2)
       vterm(k,i_mp_ni) = -rhofac_q(i_mp_qi) &
            * ( velq_l * (          weight ) &
              + velq_s * ( 1.0_rp - weight ) )

       ! QS
       rhofac_q(i_mp_qs) = rhofac_q(i_mp_qi)
       xq = max( xqmin(i_mp_qs), min( xqmax(i_mp_qs), rhoq(k,i_qs) / ( rhoq(k,i_ns) + nqmin(i_mp_qs) ) ) )

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

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

       velq_s = coef_vt0(i_mp_qs,1) * xq**beta_vn(i_mp_qs,1)
       velq_l = coef_vt0(i_mp_qs,2) * xq**beta_vn(i_mp_qs,2)
       vterm(k,i_mp_ns) = -rhofac_q(i_mp_qs) &
            * ( velq_l * (          weight ) &
              + velq_s * ( 1.0_rp - weight ) )

       ! QG
       rhofac_q(i_mp_qg) = rhofac_q(i_mp_qi)
       xq = max( xqmin(i_mp_qg), min( xqmax(i_mp_qg), rhoq(k,i_qg) / ( rhoq(k,i_ng) + nqmin(i_mp_qg) ) ) )

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

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

       velq_s = coef_vt0(i_mp_qg,1) * xq**beta_vn(i_mp_qg,1)
       velq_l = coef_vt0(i_mp_qg,2) * xq**beta_vn(i_mp_qg,2)
       vterm(k,i_mp_ng) = -rhofac_q(i_mp_qg) &
            * ( 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

    profile_stop("sn14_terminal_vel")

    return

Here is the caller graph for this function:

update_by_phase_change_kij()

subroutine scale_atmos_phy_mp_sn14::update_by_phase_change_kij	(	integer, intent(in)	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,
		integer, intent(in)	ntdiv,
		integer, intent(in)	ntmax,
		real(rp), intent(in)	dt,
		real(rp), dimension(ka,ia,ja), intent(in)	cz,
		real(rp), dimension(ka,ia,ja), intent(in)	fz,
		real(rp), dimension(ka,ia,ja), intent(in)	w,
		real(rp), dimension(ka,ia,ja), intent(in)	dTdt_rad,
		real(rp), dimension(ka,ia,ja), intent(in)	rho,
		real(rp), dimension(ka,ia,ja), intent(in)	qdry,
		real(rp), dimension(ka,ia,ja), intent(in)	esw,
		real(rp), dimension(ka,ia,ja), intent(in)	esi,
		real(rp), dimension(i_qv:i_ng,ka,ia,ja), intent(in)	rhoq2,
		real(rp), dimension(ka,ia,ja), intent(in)	pre,
		real(rp), dimension(ka,ia,ja), intent(in)	tem,
		real(rp), dimension(ka,ia,ja), intent(in)	cpa,
		real(rp), dimension(ka,ia,ja), intent(in)	cva,
		real(rp), dimension(pq_max,ka,ia,ja), intent(inout)	PQ,
		real(rp), dimension(ia,ja), intent(inout)	sl_PLCdep,
		real(rp), dimension(ia,ja), intent(inout)	sl_PLRdep,
		real(rp), dimension(ia,ja), intent(inout)	sl_PNRdep,
		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)	qc_evaporate
	)

Definition at line 4303 of file scale_atmos_phy_mp_sn14.F90.

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, scale_atmos_hydrometeor::cv_water, and scale_io::io_fid_conf.

Referenced by atmos_phy_mp_sn14_qhyd2qtrc().

    use scale_atmos_hydrometeor, only: &
       cp_vapor, &
       cp_water, &
       cp_ice,   &
       cv_vapor, &
       cv_water, &
       cv_ice
    use scale_atmos_saturation, only: &
       moist_pres2qsat_liq  => atmos_saturation_pres2qsat_liq,  &
       moist_pres2qsat_ice  => atmos_saturation_pres2qsat_ice,  &
       moist_dqs_dtem_dens_liq  => atmos_saturation_dqs_dtem_dens_liq,  &
       moist_dqs_dtem_dens_ice  => atmos_saturation_dqs_dtem_dens_ice,  &
       moist_dqs_dtem_dpre_liq => atmos_saturation_dqs_dtem_dpre_liq, &
       moist_dqs_dtem_dpre_ice => atmos_saturation_dqs_dtem_dpre_ice
    implicit none

    integer, intent(in) :: ka, ks, ke
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je

    integer, intent(in)    :: ntdiv               ! [Add] 10/08/03
    integer, intent(in)    :: ntmax               ! [Add] 10/08/03
    !
    real(RP), intent(in)    :: dt                 ! time step[s]
    !real(RP), intent(in)    :: gsgam2(KA,IA,JA)   ! metric
    real(RP), intent(in)    :: cz(ka,ia,ja)        ! altitude [m]
    real(RP), intent(in)    :: fz(ka,ia,ja)       ! altitude difference [m]
    real(RP), intent(in)    :: w(ka,ia,ja)        ! vertical velocity @ full level [m/s]
    real(RP), intent(in)    :: dtdt_rad(ka,ia,ja) ! temperture tendency by radiation[K/s]
    real(RP), intent(in)    :: rho(ka,ia,ja)      ! density[kg/m3]
    real(RP), intent(in)    :: qdry(ka,ia,ja)     ! dry air mass ratio [kg/kg]
    real(RP), intent(in)    :: esw(ka,ia,ja)      ! saturated vapor pressure for liquid
    real(RP), intent(in)    :: esi(ka,ia,ja)      !                          for ice
    real(RP), intent(in)    :: rhoq2(i_qv:i_ng,ka,ia,ja)

    real(RP), intent(in) :: tem(ka,ia,ja)      ! temperature[K]
    real(RP), intent(in) :: pre(ka,ia,ja)      ! pressure[Pa]
    real(RP), intent(in)   :: cva(ka,ia,ja)      ! specific heat at constant volume
    real(RP), intent(in)   :: cpa(ka,ia,ja)      !

    !+++ tendency[kg/m3/s]
    real(RP), intent(inout) :: pq(pq_max,ka,ia,ja)
    !+++ Column integrated tendency[kg/m2/s]
    real(RP), intent(inout) :: sl_plcdep(ia,ja)
    real(RP), intent(inout) :: sl_plrdep(ia,ja), sl_pnrdep(ia,ja)

    real(RP),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)

    !+++ tendency[kg/m3/s]
    real(RP), intent(out)   :: qc_evaporate(ka,ia,ja)

    real(RP) :: xi                     ! mean mass of ice particles
    real(RP) :: rrho                   ! 1/rho
    real(RP) :: wtem(ka,ia,ja)         ! temperature[K]
    !
    real(RP) :: r_cva                  ! specific heat at constant volume
    !real(RP) :: cpa                    ! specific heat at constant pressure
    real(RP) :: r_cpa                  ! specific heat at constant pressure
    real(RP) :: qsw(ka,ia,ja)          ! saturated mixing ratio for liquid
    real(RP) :: qsi(ka,ia,ja)          ! saturated mixing ratio for solid
    real(RP) :: dqswdtem_rho(ka,ia,ja) ! (dqsw/dtem)_rho
    real(RP) :: dqsidtem_rho(ka,ia,ja) ! (dqsi/dtem)_rho
    real(RP) :: dqswdtem_pre(ka,ia,ja) ! (dqsw/dtem)_pre
    real(RP) :: dqsidtem_pre(ka,ia,ja) ! (dqsi/dtem)_pre
    real(RP) :: dqswdpre_tem(ka,ia,ja) ! (dqsw/dpre)_tem
    real(RP) :: dqsidpre_tem(ka,ia,ja) ! (dqsi/dpre)_tem
    !
    real(RP) :: w2                     ! vetical velocity[m/s]
    real(RP) :: acnd                   ! Pdynliq + Bergeron-Findeisen
    real(RP) :: adep                   ! Pdyndep + Bergeron-Findeisen
    real(RP) :: aliqliq, asolliq
    real(RP) :: aliqsol, asolsol
    real(RP) :: pdynliq                ! production term of ssw by vertical motion
    real(RP) :: pdynsol                ! production term of ssi by vertical motion
    real(RP) :: pradliq                ! production term of ssw by radiation
    real(RP) :: pradsol                ! production term of ssi by radiation
    real(RP) :: taucnd, r_taucnd       ! time scale of ssw change by MP
    real(RP) :: taudep, r_taudep       ! time scale of ssi change by MP
    real(RP) :: taucnd_c(ka,ia,ja), r_taucnd_c  ! by cloud
    real(RP) :: taucnd_r(ka,ia,ja), r_taucnd_r  ! by rain
    real(RP) :: taudep_i(ka,ia,ja), r_taudep_i  ! by ice
    real(RP) :: taudep_s(ka,ia,ja), r_taudep_s  ! by snow
    real(RP) :: taudep_g(ka,ia,ja), r_taudep_g  ! by graupel
    ! alternative tendency through changing ssw and ssi
    real(RP) :: pncdep ! [Add] 11/08/30 T.Mitsui
    real(RP) :: plr2nr, pli2ni, pls2ns, plg2ng
    real(RP) :: coef_a_cnd, coef_b_cnd
    real(RP) :: coef_a_dep, coef_b_dep
    !
    real(RP) :: frz_dqc
    real(RP) :: frz_dnc
    real(RP) :: frz_dqr
    real(RP) :: frz_dnr
    real(RP) :: mlt_dqi
    real(RP) :: mlt_dni
    real(RP) :: mlt_dqs
    real(RP) :: mlt_dns
    real(RP) :: mlt_dqg
    real(RP) :: mlt_dng
    real(RP) :: dep_dqi
    real(RP) :: dep_dni
    real(RP) :: dep_dqs
    real(RP) :: dep_dns
    real(RP) :: dep_dqg
    real(RP) :: dep_dng
    real(RP) :: dep_dqr
    real(RP) :: dep_dnr
    real(RP) :: dep_dqc
    real(RP) :: dep_dnc   ! 11/08/30 [Add] T.Mitsui, dep_dnc
    real(RP) :: r_xc_ccn, r_xi_ccn ! 11/08/30 [Add] T.Mitsui
    !
    real(RP) :: drhoqv
    real(RP) :: drhoqc, drhoqr, drhoqi, drhoqs, drhoqg
    real(RP) :: drhonc, drhonr, drhoni, drhons, drhong
    !
    real(RP) :: fac1, fac2, fac3, fac4, fac5, fac6
    real(RP) :: r_rvaptem        ! 1/(Rvap*tem)
    real(RP) :: pv               ! vapor pressure
    real(RP) :: lvsw, lvsi       ! saturated vapor density
    real(RP) :: dlvsw, dlvsi     !
    ! [Add] 11/08/30 T.Mitsui
    real(RP) :: dcnd, ddep       ! total cndensation/deposition
    real(RP) :: uplim_cnd        ! upper limit of condensation
    real(RP) :: lowlim_cnd       ! lower limit of evaporation
    ! [Add] 11/08/30 T.Mitsui
    real(RP) :: uplim_dep        ! upper limit of condensation
    real(RP) :: lowlim_dep       ! lower limit of evaporation
    real(RP) :: ssw, ssi         ! supersaturation ratio
    real(RP) :: r_esw, r_esi     ! 1/esw, 1/esi
    real(RP) :: r_lvsw, r_lvsi   ! 1/(lvsw*ssw), 1/(lvsi*ssi)
    real(RP) :: r_dt             ! 1/dt
    real(RP) :: ssw_o, ssi_o
!    real(RP) :: dt_dyn
!    real(RP) :: dt_mp
    !
!    real(RP)       :: tem_lh(KA,IA,JA)
!    real(RP)       :: dtemdt_lh(KA,IA,JA)
    real(RP), save :: fac_cndc        = 1.0_rp
    logical, save :: opt_fix_taucnd_c=.false.
    logical, save :: flag_first      =.true.
    !
    namelist / param_atmos_phy_mp_sn14_condensation / &
         opt_fix_taucnd_c, fac_cndc

    real(RP) :: fac_cndc_wrk
    !
    real(RP), parameter :: tau100day   = 1.e+7_rp
    real(RP), parameter :: r_tau100day = 1.e-7_rp
    real(RP), parameter :: eps=1.e-30_rp
    !
    real(RP) :: dz
    !
    integer :: i,j,k,iqw
    real(RP) :: sw
    real(RP) :: dqv, dqc, dqr, dqi, dqs, dqg, dcv, dcp
    !

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

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

    call moist_pres2qsat_liq( ka, ks, ke, ia, is, ie, ja, js, je, &
                              wtem(:,:,:), pre(:,:,:), qdry(:,:,:), & ! [IN]
                              qsw(:,:,:)                            ) ! [OUT]
    call moist_pres2qsat_ice( ka, ks, ke, ia, is, ie, ja, js, je, &
                              wtem(:,:,:), pre(:,:,:), qdry(:,:,:), & ! [IN]
                              qsi(:,:,:)                            ) ! [OUT]
    call moist_dqs_dtem_dens_liq( ka, ks, ke, ia, is, ie, ja, js, je, &
                                  wtem(:,:,:), rho(:,:,:), & ! [IN]
                                  dqswdtem_rho(:,:,:)      ) ! [OUT]
    call moist_dqs_dtem_dens_ice( ka, ks, ke, ia, is, ie, ja, js, je, &
                                  wtem(:,:,:), rho(:,:,:), & ! [IN]
                                  dqsidtem_rho(:,:,:)      ) ! [OUT]
    call moist_dqs_dtem_dpre_liq( ka, ks, ke, ia, is, ie, ja, js, je, &
                                  wtem(:,:,:), pre(:,:,:), qdry(:,:,:),    & ! [IN]
                                  dqswdtem_pre(:,:,:), dqswdpre_tem(:,:,:) ) ! [OUT]
    call moist_dqs_dtem_dpre_ice( ka, ks, ke, ia, is, ie, ja, js, je, &
                                  wtem(:,:,:), pre(:,:,:), qdry(:,:,:),    & ! [IN]
                                  dqsidtem_pre(:,:,:), dqsidpre_tem(:,:,:) ) ! [OUT]

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

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

       r_cva = 1.0_rp / cva(k,i,j)
       r_cpa = 1.0_rp / cpa(k,i,j)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

       ! tendency
       rhoq_t(k,i,j,i_qv) = drhoqv / dt
       rhoq_t(k,i,j,i_qc) = drhoqc / dt
       rhoq_t(k,i,j,i_nc) = drhonc / dt
       rhoq_t(k,i,j,i_qr) = drhoqr / dt
       rhoq_t(k,i,j,i_nr) = drhonr / dt
       rhoq_t(k,i,j,i_qi) = drhoqi / dt
       rhoq_t(k,i,j,i_ni) = drhoni / dt
       rhoq_t(k,i,j,i_qs) = drhoqs / dt
       rhoq_t(k,i,j,i_ns) = drhons / dt
       rhoq_t(k,i,j,i_qg) = drhoqg / dt
       rhoq_t(k,i,j,i_ng) = drhong / dt

       rhoe_t(k,i,j) = ( - lhv * drhoqv + lhf * ( drhoqi+ drhoqs + drhoqg ) ) / dt

       rrho = 1.0_rp/rho(k,i,j)
       dqv = drhoqv * rrho
       dqc = drhoqc * rrho
       dqr = drhoqr * rrho
       dqi = drhoqi * rrho
       dqs = drhoqs * rrho
       dqg = drhoqg * rrho

       dcv = cv_vapor * dqv + cv_water * ( dqc + dqr ) + cv_ice * ( dqi + dqs + dqg )
       dcp = cp_vapor * dqv + cp_water * ( dqc + dqr ) + cp_ice * ( dqi + dqs + dqg )

       cvtot_t(k,i,j) = dcv/dt
       cptot_t(k,i,j) = dcp/dt

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

    return

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

Referenced by atmos_phy_mp_sn14_cloud_fraction(), atmos_phy_mp_sn14_qhyd2qtrc(), and atmos_phy_mp_sn14_terminal_velocity().

  integer, public, parameter :: qa_mp  = 11

atmos_phy_mp_sn14_ntracers

integer, parameter, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_ntracers = QA_MP

Definition at line 121 of file scale_atmos_phy_mp_sn14.F90.

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

  integer,                parameter, public :: atmos_phy_mp_sn14_ntracers = qa_mp

atmos_phy_mp_sn14_nwaters

integer, parameter, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_nwaters = 2

Definition at line 122 of file scale_atmos_phy_mp_sn14.F90.

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

  integer,                parameter, public :: atmos_phy_mp_sn14_nwaters = 2

atmos_phy_mp_sn14_nices

integer, parameter, public scale_atmos_phy_mp_sn14::atmos_phy_mp_sn14_nices = 3

Definition at line 123 of file scale_atmos_phy_mp_sn14.F90.

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

  integer,                parameter, public :: atmos_phy_mp_sn14_nices = 3

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

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

  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'  /)

atmos_phy_mp_sn14_tracer_descriptions

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

Definition at line 136 of file scale_atmos_phy_mp_sn14.F90.

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

  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 miass ratio to total mass                    ', &
       'Ratio of Graupel mass ratio to total mass                  ', &
       'Cloud Water Number Density                                 ', &
       'Rain Water Number Density                                  ', &
       'Cloud Ice Number Density                                   ', &
       'Snow Number Density                                        ', &
       'Graupel Number Density                                     '/)

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

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

  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'  /)