SCALE-RM
Functions/Subroutines | Variables
scale_atmos_phy_mp_sn14 Module Reference

module ATMOSPHERE / Physics Cloud Microphysics More...

Functions/Subroutines

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

Variables

real(rp), dimension(mp_qa), target, public atmos_phy_mp_dens
 

Detailed Description

module ATMOSPHERE / Physics Cloud Microphysics

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

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

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

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 ( character(len=*), intent(in)  MP_TYPE)

Setup Cloud Microphysics.

Definition at line 488 of file scale_atmos_phy_mp_sn14.F90.

References atmos_phy_mp_dens, scale_const::const_dice, scale_const::const_dwatr, scale_grid::grid_cdz, scale_grid_index::ia, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_l, scale_stdio::io_lnml, scale_grid_index::ja, scale_grid_index::ka, scale_process::prc_mpistop(), and scale_time::time_dtsec_atmos_phy_mp.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_setup().

488  use scale_process, only: &
489  prc_mpistop
490  use scale_grid, only: &
491  cdz => grid_cdz
492  use scale_const, only: &
493  const_dwatr, &
494  const_dice
495  use scale_time, only: &
496  time_dtsec_atmos_phy_mp
497  implicit none
498 
499  character(len=*), intent(in) :: mp_type
500 
501  namelist / param_atmos_phy_mp / &
502  mp_doautoconversion, &
503  mp_doprecipitation, &
504  mp_ssw_lim, &
505  mp_couple_aerosol, &
506  mp_ntmax_sedimentation
507 
508  real(RP), parameter :: max_term_vel = 10.0_rp !-- terminal velocity for calculate dt of sedimentation
509  integer :: nstep_max
510  integer :: ierr
511  !---------------------------------------------------------------------------
512 
513  if( io_l ) write(io_fid_log,*)
514  if( io_l ) write(io_fid_log,*) '+++ Module[Cloud Microphisics]/Categ[ATMOS]'
515  if( io_l ) write(io_fid_log,*) '*** Wrapper for SN14'
516 
517  if ( mp_type /= 'SN14' ) then
518  write(*,*) 'xxx ATMOS_PHY_MP_TYPE is not SN14. Check!'
519  call prc_mpistop
520  end if
521 
522  if ( i_qv <= 0 &
523  .OR. i_qc <= 0 &
524  .OR. i_qr <= 0 &
525  .OR. i_qi <= 0 &
526  .OR. i_qs <= 0 &
527  .OR. i_qg <= 0 &
528  .OR. i_nc <= 0 &
529  .OR. i_nr <= 0 &
530  .OR. i_ni <= 0 &
531  .OR. i_ns <= 0 &
532  .OR. i_ng <= 0 ) then
533  write(*,*) 'xxx SN14 needs QV/C/R/I/S/G and NC/R/I/S/G tracer. Check!'
534  call prc_mpistop
535  endif
536 
537  !--- read namelist
538  rewind(io_fid_conf)
539  read(io_fid_conf,nml=param_atmos_phy_mp,iostat=ierr)
540  if( ierr < 0 ) then !--- missing
541  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
542  elseif( ierr > 0 ) then !--- fatal error
543  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP. Check!'
544  call prc_mpistop
545  endif
546  if( io_lnml ) write(io_fid_log,nml=param_atmos_phy_mp)
547 
548  atmos_phy_mp_dens(i_mp_qc) = const_dwatr
549  atmos_phy_mp_dens(i_mp_qr) = const_dwatr
550  atmos_phy_mp_dens(i_mp_qi) = const_dice
551  atmos_phy_mp_dens(i_mp_qs) = const_dice
552  atmos_phy_mp_dens(i_mp_qg) = const_dice
553 
554  wlabel( 1) = "VAPOR"
555  wlabel( 2) = "CLOUD"
556  wlabel( 3) = "RAIN"
557  wlabel( 4) = "ICE"
558  wlabel( 5) = "SNOW"
559  wlabel( 6) = "GRAUPEL"
560  wlabel( 7) = "CLOUD_NUM"
561  wlabel( 8) = "RAIN_NUM"
562  wlabel( 9) = "ICE_NUM"
563  wlabel(10) = "SNOW_NUM"
564  wlabel(11) = "GRAUPEL_NUM"
565 
566  call mp_sn14_init
567 
568  allocate(nc_uplim_d(1,ia,ja))
569  nc_uplim_d(:,:,:) = 150.e6_rp
570 
571  nstep_max = int( ( time_dtsec_atmos_phy_mp * max_term_vel ) / minval( cdz ) )
572  mp_ntmax_sedimentation = max( mp_ntmax_sedimentation, nstep_max )
573 
574  mp_nstep_sedimentation = mp_ntmax_sedimentation
575  mp_rnstep_sedimentation = 1.0_rp / real(mp_ntmax_sedimentation,kind=rp)
576  mp_dtsec_sedimentation = time_dtsec_atmos_phy_mp * mp_rnstep_sedimentation
577 
578  if ( io_l ) write(io_fid_log,*)
579  if ( io_l ) write(io_fid_log,*) '*** Timestep of sedimentation is divided into : ', mp_ntmax_sedimentation, ' step'
580  if ( io_l ) write(io_fid_log,*) '*** DT of sedimentation is : ', mp_dtsec_sedimentation, '[s]'
581 
582  !--- For kij
583  allocate( gsgam2_d(ka,ia,ja) )
584  allocate( gsgam2h_d(ka,ia,ja) )
585  allocate( gam2_d(ka,ia,ja) )
586  allocate( gam2h_d(ka,ia,ja) )
587  allocate( rgsgam2_d(ka,ia,ja) )
588  allocate( rgs_d(ka,ia,ja) )
589  allocate( rgsh_d(ka,ia,ja) )
590  gsgam2_d(:,:,:) = 1.0_rp
591  gsgam2h_d(:,:,:) = 1.0_rp
592  gam2_d(:,:,:) = 1.0_rp
593  gam2h_d(:,:,:) = 1.0_rp
594  rgsgam2_d(:,:,:) = 1.0_rp
595  rgs_d(:,:,:) = 1.0_rp
596  rgsh_d(:,:,:) = 1.0_rp
597 
598  return
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:234
scale_const::const_dwatr
real(rp), parameter, public const_dwatr
density of water [kg/m3]
Definition: scale_const.F90:87
scale_time::time_dtsec_atmos_phy_mp
real(dp), public time_dtsec_atmos_phy_mp
time interval of physics(microphysics) [sec]
Definition: scale_time.F90:41
scale_const::const_dice
real(rp), parameter, public const_dice
density of ice [kg/m3]
Definition: scale_const.F90:88
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_atmos_phy_mp_dry::atmos_phy_mp_dens
real(rp), dimension(mp_qa), target, public atmos_phy_mp_dens
Definition: scale_atmos_phy_mp_dry.F90:44
scale_time
module TIME
Definition: scale_time.F90:15
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_grid
module GRID (cartesian)
Definition: scale_grid_cartesian.F90:17
scale_grid::grid_cdz
real(rp), dimension(:), allocatable, public grid_cdz
z-length of control volume [m]
Definition: scale_grid_cartesian.F90:59
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:43
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_mp_sn14()

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

Cloud Microphysics.

Definition at line 615 of file scale_atmos_phy_mp_sn14.F90.

References scale_atmos_thermodyn::aq_cp, scale_atmos_thermodyn::aq_cv, scale_atmos_phy_mp_common::atmos_phy_mp_precipitation(), aut_acc_slc_brk_kij(), debug_tem_kij(), scale_grid::dz, freezing_water_kij(), scale_grid::grid_cdz, scale_grid::grid_cz, ice_multiplication_kij(), scale_grid_index::ie, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_l, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, dc_log::log(), mixed_phase_collection_kij(), mp_negativefilter(), scale_tracer::mp_qa, nucleation_kij(), scale_process::prc_mpistop(), scale_prof::prof_rapend(), scale_prof::prof_rapstart(), scale_tracer::qa, scale_specfunc::sf_gamma(), scale_time::time_dtsec_atmos_phy_mp, and update_by_phase_change_kij().

Referenced by scale_atmos_phy_mp::atmos_phy_mp_setup().

615  use scale_grid_index
616  use scale_tracer, only: &
617  qad => qa, &
618  mp_qad => mp_qa
619  implicit none
620 
621  real(RP), intent(inout) :: dens(ka,ia,ja)
622  real(RP), intent(inout) :: momz(ka,ia,ja)
623  real(RP), intent(inout) :: momx(ka,ia,ja)
624  real(RP), intent(inout) :: momy(ka,ia,ja)
625  real(RP), intent(inout) :: rhot(ka,ia,ja)
626  real(RP), intent(inout) :: qtrc(ka,ia,ja,qad)
627  real(RP), intent(in) :: ccn(ka,ia,ja)
628  real(RP), intent(out) :: evaporate(ka,ia,ja)
629  real(RP), intent(out) :: sflx_rain(ia,ja)
630  real(RP), intent(out) :: sflx_snow(ia,ja)
631  !---------------------------------------------------------------------------
632 
633  if( io_l ) write(io_fid_log,*) '*** Physics step: Cloud microphysics(SN14)'
634 
635 #ifdef PROFILE_FIPP
636  call fipp_start()
637 #endif
638 
639  call mp_negativefilter( dens, qtrc )
640 
641  call mp_sn14( dens, & ! [INOUT]
642  momz, & ! [INOUT]
643  momx, & ! [INOUT]
644  momy, & ! [INOUT]
645  rhot, & ! [INOUT]
646  qtrc, & ! [INOUT]
647  ccn, & ! [IN]
648  evaporate, & ! [OUT]
649  sflx_rain, & ! [OUT]
650  sflx_snow ) ! [OUT]
651 
652  call mp_negativefilter( dens, qtrc )
653 
654 #ifdef PROFILE_FIPP
655  call fipp_stop()
656 #endif
657 
658  return
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:18
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_tracer::mp_qa
integer, public mp_qa
Definition: scale_tracer.F90:98
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the call graph for this function:
Here is the caller graph for this function:

◆ debug_tem_kij()

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

Definition at line 2063 of file scale_atmos_phy_mp_sn14.F90.

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

Referenced by atmos_phy_mp_sn14().

2063  use scale_process, only: &
2064  prc_myrank
2065  implicit none
2066 
2067  integer, intent(in) :: point
2068  real(RP), intent(in) :: tem(ka,ia,ja)
2069  real(RP), intent(in) :: rho(ka,ia,ja)
2070  real(RP), intent(in) :: pre(ka,ia,ja)
2071  real(RP), intent(in) :: qv (ka,ia,ja)
2072 
2073  integer :: k ,i, j
2074  !---------------------------------------------------------------------------
2075 
2076  do j = js, je
2077  do i = is, ie
2078  do k = ks, ke
2079  if ( tem(k,i,j) < tem_min &
2080  .OR. rho(k,i,j) < rho_min &
2081  .OR. pre(k,i,j) < 1.0_rp ) then
2082 
2083  if( io_l ) write(io_fid_log,'(A,I3,A,4(F16.5),3(I6))') &
2084  "*** 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
2085  endif
2086  enddo
2087  enddo
2088  enddo
2089 
2090  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_process::prc_myrank
integer, public prc_myrank
process num in local communicator
Definition: scale_process.F90:78
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the caller graph for this function:

◆ nucleation_kij()

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

Definition at line 2103 of file scale_atmos_phy_mp_sn14.F90.

References scale_atmos_saturation::atmos_saturation_dqsi_dtem_rho(), scale_grid_index::ie, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_l, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, and scale_process::prc_mpistop().

Referenced by atmos_phy_mp_sn14().

2103  use scale_process, only: &
2104  prc_mpistop
2105  use scale_atmos_saturation, only: &
2106  moist_psat_liq => atmos_saturation_psat_liq, &
2107  moist_psat_ice => atmos_saturation_psat_ice, &
2108  moist_pres2qsat_liq => atmos_saturation_pres2qsat_liq, &
2109  moist_pres2qsat_ice => atmos_saturation_pres2qsat_ice, &
2110  moist_dqsi_dtem_rho => atmos_saturation_dqsi_dtem_rho
2111  implicit none
2112 
2113  real(RP), intent(in) :: z(ka) !
2114  real(RP), intent(in) :: velz(ka,ia,ja) ! w of half point
2115  real(RP), intent(in) :: rho(ka,ia,ja) ! [Add] 09/08/18 T.Mitsui
2116  real(RP), intent(in) :: tem(ka,ia,ja) ! [Add] 09/08/18 T.Mitsui
2117  real(RP), intent(in) :: pre(ka,ia,ja) ! [Add] 09/08/18 T.Mitsui
2118  !
2119  real(RP), intent(in) :: rhoq(qa,ka,ia,ja) !
2120  real(RP), intent(out) :: pq(pq_max,ka,ia,ja)
2121  !
2122  real(RP), intent(in) :: cpa(ka,ia,ja) ! in 09/08/18 [Add] T.Mitsui
2123  real(RP), intent(in) :: dtdt_rad(ka,ia,ja) ! 09/08/18 T.Mitsui
2124  real(RP), intent(in) :: qke(ka,ia,ja) ! 09/08/18 T.Mitsui
2125  real(DP), intent(in) :: dt
2126  real(RP), intent(in) :: ccn(ka,ia,ja)
2127  !
2128  ! namelist variables
2129  !
2130  ! total aerosol number concentration [/m3]
2131  real(RP), parameter :: c_ccn_ocean= 1.00e+8_rp
2132  real(RP), parameter :: c_ccn_land = 1.26e+9_rp
2133  real(RP), save :: c_ccn = 1.00e+8_rp
2134  ! aerosol activation factor
2135  real(RP), parameter :: kappa_ocean= 0.462_rp
2136  real(RP), parameter :: kappa_land = 0.308_rp
2137  real(RP), save :: kappa = 0.462_rp
2138  real(RP), save :: c_in = 1.0_rp
2139  ! SB06 (36)
2140  real(RP), save :: nm_m92 = 1.e+3_rp
2141  real(RP), save :: am_m92 = -0.639_rp
2142  real(RP), save :: bm_m92 = 12.96_rp
2143  !
2144  real(RP), save :: in_max = 1000.e+3_rp ! max num. of Ice-Nuclei [num/m3]
2145  real(RP), save :: ssi_max= 0.60_rp
2146  real(RP), save :: ssw_max= 1.1_rp ! [%]
2147  !
2148  logical, save :: flag_first = .true.
2149  real(RP), save :: qke_min = 0.03_rp ! sigma=0.1[m/s], 09/08/18 T.Mitsui
2150  real(RP), save :: tem_ccn_low=233.150_rp ! = -40 degC ! [Add] 10/08/03 T.Mitsui
2151  real(RP), save :: tem_in_low =173.150_rp ! = -100 degC ! [Add] 10/08/03 T.Mitsui
2152  logical, save :: nucl_twomey = .false.
2153  logical, save :: inucl_w = .false.
2154  !
2155  namelist /nm_mp_sn14_nucleation/ &
2156  in_max, & !
2157  c_ccn, kappa, & ! cloud nucleation
2158  nm_m92, am_m92, bm_m92, & ! ice nucleation
2159  xc_ccn, xi_ccn, &
2160  tem_ccn_low, & ! [Add] 10/08/03 T.Mitsui
2161  tem_in_low, & ! [Add] 10/08/03 T.Mitsui
2162  ssw_max, ssi_max, &
2163  nucl_twomey, inucl_w ! [Add] 13/01/30 Y.Sato
2164  !
2165 ! real(RP) :: c_ccn_map(1,IA,JA) ! c_ccn horizontal distribution
2166 ! real(RP) :: kappa_map(1,IA,JA) ! kappa horizontal distribution
2167 ! real(RP) :: c_in_map(1,IA,JA) ! c_in horizontal distribution ! [Add] 11/08/30 T.Mitsui
2168  real(RP) :: esw(ka,ia,ja) ! saturation vapor pressure, water
2169  real(RP) :: esi(ka,ia,ja) ! ice
2170  real(RP) :: ssw(ka,ia,ja) ! super saturation (water)
2171  real(RP) :: ssi(ka,ia,ja) ! super saturation (ice)
2172 ! real(RP) :: w_dsswdz(KA,IA,JA) ! w*(d_ssw/ d_z) super saturation(water) flux
2173  real(RP) :: w_dssidz(ka,ia,ja) ! w*(d_ssi/ d_z), 09/04/14 T.Mitsui
2174 ! real(RP) :: ssw_below(KA,IA,JA)! ssw(k-1)
2175  real(RP) :: ssi_below(ka,ia,ja)! ssi(k-1), 09/04/14 T.Mitsui
2176  real(RP) :: z_below(ka,ia,ja) ! z(k-1)
2177  real(RP) :: dz ! z(k)-z(k-1)
2178  real(RP) :: pv ! vapor pressure
2179  ! work variables for Twomey Equation.
2180  real(RP) :: qsw(ka,ia,ja)
2181  real(RP) :: qsi(ka,ia,ja)
2182  real(RP) :: dqsidtem_rho(ka,ia,ja)
2183  real(RP) :: dssidt_rad(ka,ia,ja)
2184 ! real(RP) :: dni_ratio(KA,IA,JA)
2185  real(RP) :: wssi, wdssi
2186  !
2187 ! real(RP) :: xi_nuc(1,IA,JA) ! xi use the value @ cloud base
2188 ! real(RP) :: alpha_nuc(1,IA,JA) ! alpha_nuc
2189 ! real(RP) :: eta_nuc(1,IA,JA) ! xi use the value @ cloud base
2190  !
2191  real(RP) :: sigma_w(ka,ia,ja)
2192  real(RP) :: weff(ka,ia,ja)
2193  real(RP) :: weff_max(ka,ia,ja)
2194  !
2195  real(RP) :: coef_ccn(ia,ja)
2196  real(RP) :: slope_ccn(ia,ja)
2197  real(RP) :: nc_new(ka,ia,ja)
2198  real(RP) :: nc_new_below(ka,ia,ja)
2199  real(RP) :: dnc_new
2200  real(RP) :: nc_new_max ! Lohmann (2002),
2201  real(RP) :: a_max
2202  real(RP) :: b_max
2203  logical :: flag_nucleation(ka,ia,ja)
2204  !
2205  real(RP) :: r_gravity
2206  real(RP), parameter :: r_sqrt3=0.577350269_rp ! = sqrt(1.d0/3.d0)
2207  real(RP), parameter :: eps=1.e-30_rp
2208  !====> ! 09/08/18
2209  !
2210  real(RP) :: dlcdt_max, dli_max ! defined by supersaturation
2211  real(RP) :: dncdt_max, dni_max ! defined by supersaturation
2212  real(RP) :: rdt
2213  !
2214  integer :: i, j, k
2215  !
2216  if( flag_first )then
2217  rewind(io_fid_conf)
2218  read(io_fid_conf, nml=nm_mp_sn14_nucleation, end=100)
2219 100 if( io_l ) write(io_fid_log, nml=nm_mp_sn14_nucleation)
2220  flag_first=.false.
2221  if( mp_couple_aerosol .AND. nucl_twomey ) then
2222  write(io_fid_log,*) "nucl_twomey should be false when MP_couple_aerosol is true, stop"
2223  call prc_mpistop
2224  endif
2225  endif
2226  !
2227 ! c_ccn_map(1,:,:) = c_ccn
2228 ! kappa_map(1,:,:) = kappa
2229 ! c_in_map(1,:,:) = c_in
2230  !
2231 ! nc_uplim_d(1,:,:) = c_ccn_map(1,:,:)*1.5_RP
2232  do j = js, je
2233  do i = is, ie
2234  nc_uplim_d(1,i,j) = c_ccn*1.5_rp
2235  end do
2236  end do
2237  !
2238  rdt = 1.0_rp/dt
2239  r_gravity = 1.0_rp/grav
2240  !
2241  call moist_psat_liq ( esw, tem )
2242  call moist_psat_ice ( esi, tem )
2243  call moist_pres2qsat_liq( qsw, tem, pre )
2244  call moist_pres2qsat_ice( qsi, tem, pre )
2245  call moist_dqsi_dtem_rho( dqsidtem_rho, tem, rho )
2246  !
2247  ! Lohmann (2002),JAS, eq.(1) but changing unit [cm-3] => [m-3]
2248  a_max = 1.e+6_rp*0.1_rp*(1.e-6_rp)**1.27_rp
2249  b_max = 1.27_rp
2250  !
2251  ssi_max = 1.0_rp
2252 
2253  do j = js, je
2254  do i = is, ie
2255  do k = ks, ke
2256  pv = rhoq(i_qv,k,i,j)*rvap*tem(k,i,j)
2257  ssw(k,i,j) = min( mp_ssw_lim, ( pv/esw(k,i,j)-1.0_rp ) )*100.0_rp
2258  ssi(k,i,j) = (pv/esi(k,i,j) - 1.00_rp)
2259 ! ssw_below(k+1,i,j) = ssw(k,i,j)
2260  ssi_below(k+1,i,j) = ssi(k,i,j)
2261  z_below(k+1,i,j) = z(k)
2262  end do
2263 ! ssw_below(KS,i,j) = ssw(KS,i,j)
2264  ssi_below(ks,i,j) = ssi(ks,i,j)
2265  z_below(ks,i,j) = z(ks-1)
2266 
2267  ! dS/dz is evaluated by first order upstream difference
2268  !*** Solution for Twomey Equation ***
2269 ! 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)) * &
2270  coef_ccn(i,j) = 1.e+6_rp*0.88_rp*(c_ccn*1.e-6_rp)**(2.0_rp/(kappa + 2.0_rp)) &
2271 ! * (70.0_RP)**(kappa_map(1,i,j)/(kappa_map(1,i,j) + 2.0_RP))
2272  * (70.0_rp)**(kappa/(kappa + 2.0_rp))
2273 ! slope_ccn(i,j) = 1.5_RP*kappa_map(1,i,j)/(kappa_map(1,i,j) + 2.0_RP)
2274  slope_ccn(i,j) = 1.5_rp*kappa/(kappa + 2.0_rp)
2275  !
2276  do k=ks, ke
2277  sigma_w(k,i,j) = r_sqrt3*sqrt(max(qke(k,i,j),qke_min))
2278  end do
2279  sigma_w(ks-1,i,j) = sigma_w(ks,i,j)
2280  sigma_w(ke+1,i,j) = sigma_w(ke,i,j)
2281  ! effective vertical velocity
2282  do k=ks, ke-1
2283  weff(k,i,j) = 0.5_rp*(velz(k,i,j) + velz(k+1,i,j)) - cpa(k,i,j)*r_gravity*dtdt_rad(k,i,j)
2284  end do
2285  weff(ks-1,i,j) = weff(ks,i,j)
2286  weff(ke,i,j) = weff(ke-1,i,j)
2287 
2288  end do
2289  end do
2290  !
2291  if( mp_couple_aerosol ) then
2292 
2293  do j = js, je
2294  do i = is, ie
2295  do k = ks, ke
2296  if( ssw(k,i,j) > 1.e-10_rp .AND. pre(k,i,j) > 300.e+2_rp ) then
2297  nc_new(k,i,j) = max( ccn(k,i,j), c_ccn )
2298  else
2299  nc_new(k,i,j) = 0.0_rp
2300  endif
2301  enddo
2302  enddo
2303  enddo
2304 
2305  else
2306 
2307  if( nucl_twomey ) then
2308  ! diagnose cloud condensation nuclei
2309 
2310  do j = js, je
2311  do i = is, ie
2312  do k = ks, ke
2313  ! effective vertical velocity (maximum vertical velocity in turbulent flow)
2314  weff_max(k,i,j) = weff(k,i,j) + sigma_w(k,i,j)
2315  ! large scale upward motion region and saturated
2316  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
2317  ! Lohmann (2002), eq.(1)
2318  nc_new_max = coef_ccn(i,j)*weff_max(k,i,j)**slope_ccn(i,j)
2319  nc_new(k,i,j) = a_max*nc_new_max**b_max
2320  else
2321  nc_new(k,i,j) = 0.0_rp
2322  end if
2323  end do
2324  end do
2325  end do
2326  else
2327  ! calculate cloud condensation nuclei
2328  do j = js, je
2329  do i = is, ie
2330  do k = ks, ke
2331  if( ssw(k,i,j) > 1.e-10_rp .AND. pre(k,i,j) > 300.e+2_rp ) then
2332  nc_new(k,i,j) = c_ccn*ssw(k,i,j)**kappa
2333  else
2334  nc_new(k,i,j) = 0.0_rp
2335  endif
2336  enddo
2337  enddo
2338  enddo
2339  endif
2340 
2341  endif
2342 
2343  do j = js, je
2344  do i = is, ie
2345  do k = ks, ke
2346  ! nc_new is bound by upper limit
2347  if( nc_new(k,i,j) > nc_uplim_d(1,i,j) )then ! no more CCN
2348  flag_nucleation(k,i,j) = .false.
2349  nc_new_below(k+1,i,j) = 1.e+30_rp
2350  else if( nc_new(k,i,j) > eps )then ! nucleation can occur
2351  flag_nucleation(k,i,j) = .true.
2352  nc_new_below(k+1,i,j) = nc_new(k,i,j)
2353  else ! nucleation cannot occur(unsaturated or negative w)
2354  flag_nucleation(k,i,j) = .false.
2355  nc_new_below(k+1,i,j) = 0.0_rp
2356  end if
2357  end do
2358  nc_new_below(ks,i,j) = 0.0_rp
2359 ! do k=KS, KE
2360  ! search maximum value of nc_new
2361 ! if( ( nc_new(k,i,j) < nc_new_below(k,i,j) ) .OR. &
2362 ! ( nc_new_below(k,i,j) > c_ccn_map(1,i,j)*0.05_RP ) )then ! 5% of c_ccn
2363 ! ( nc_new_below(k,i,j) > c_ccn*0.05_RP ) )then ! 5% of c_ccn
2364 ! flag_nucleation(k,i,j) = .false.
2365 ! end if
2366 ! end do
2367 
2368  end do
2369  end do
2370 
2371  if( mp_couple_aerosol ) then
2372 
2373  do j = js, je
2374  do i = is, ie
2375  do k = ks, ke
2376  ! nucleation occurs at only cloud base.
2377  ! if CCN is more than below parcel, nucleation newly occurs
2378  ! effective vertical velocity
2379  if ( flag_nucleation(k,i,j) .AND. & ! large scale upward motion region and saturated
2380  tem(k,i,j) > tem_ccn_low ) then
2381  dlcdt_max = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
2382  dncdt_max = dlcdt_max/xc_min
2383 ! dnc_new = nc_new(k,i,j)-rhoq(I_NC,k,i,j)
2384  dnc_new = nc_new(k,i,j)
2385  pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
2386  pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
2387  else
2388  pq(i_ncccn,k,i,j) = 0.0_rp
2389  pq(i_lcccn,k,i,j) = 0.0_rp
2390  end if
2391  end do
2392  end do
2393  end do
2394 
2395  else
2396 
2397  if( nucl_twomey ) then
2398  do j = js, je
2399  do i = is, ie
2400  do k = ks, ke
2401  ! nucleation occurs at only cloud base.
2402  ! if CCN is more than below parcel, nucleation newly occurs
2403  ! effective vertical velocity
2404  if ( flag_nucleation(k,i,j) .AND. & ! large scale upward motion region and saturated
2405  tem(k,i,j) > tem_ccn_low .AND. &
2406  nc_new(k,i,j) > rhoq(i_nc,k,i,j) ) then
2407  dlcdt_max = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
2408  dncdt_max = dlcdt_max/xc_min
2409  dnc_new = nc_new(k,i,j)-rhoq(i_nc,k,i,j)
2410  pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
2411  pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
2412  else
2413  pq(i_ncccn,k,i,j) = 0.0_rp
2414  pq(i_lcccn,k,i,j) = 0.0_rp
2415  end if
2416  end do
2417  end do
2418  end do
2419  else
2420  do j = js, je
2421  do i = is, ie
2422  do k = ks, ke
2423  ! effective vertical velocity
2424  if( tem(k,i,j) > tem_ccn_low .AND. &
2425  nc_new(k,i,j) > rhoq(i_nc,k,i,j) ) then
2426  dlcdt_max = (rhoq(i_qv,k,i,j) - esw(k,i,j)/(rvap*tem(k,i,j)))*rdt
2427  dncdt_max = dlcdt_max/xc_min
2428  dnc_new = nc_new(k,i,j)-rhoq(i_nc,k,i,j)
2429  pq(i_ncccn,k,i,j) = min( dncdt_max, dnc_new*rdt )
2430  pq(i_lcccn,k,i,j) = min( dlcdt_max, xc_min*pq(i_ncccn,k,i,j) )
2431  else
2432  pq(i_ncccn,k,i,j) = 0.0_rp
2433  pq(i_lcccn,k,i,j) = 0.0_rp
2434  end if
2435  end do
2436  end do
2437  end do
2438  endif
2439  endif
2440 
2441  !
2442  ! ice nucleation
2443  !
2444  ! +++ NOTE ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2445  ! Based on Phillips etal.(2006).
2446  ! However this approach doesn't diagnose Ni itself but diagnose tendency.
2447  ! Original approach adjust Ni instantaneously .
2448  ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2449  do j = js, je
2450  do i = is, ie
2451  do k = ks, ke
2452  dz = z(k) - z_below(k,i,j)
2453  w_dssidz(k,i,j) = velz(k,i,j)*(ssi(k,i,j) - ssi_below(k,i,j))/dz ! 09/04/14 [Add] T.Mitsui
2454  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)
2455  dli_max = (rhoq(i_qv,k,i,j) - esi(k,i,j)/(rvap*tem(k,i,j)))*rdt
2456  dni_max = min( dli_max/xi_ccn, (in_max-rhoq(i_ni,k,i,j))*rdt )
2457  wdssi = min( w_dssidz(k,i,j)+dssidt_rad(k,i,j), 0.01_rp)
2458  wssi = min( ssi(k,i,j), ssi_max)
2459  ! SB06(34),(35)
2460  if( ( wdssi > eps ) .AND. & !
2461  (tem(k,i,j) < 273.15_rp ) .AND. & !
2462  (rhoq(i_ni,k,i,j) < in_max ) .AND. &
2463  (wssi >= eps ) )then !
2464 ! 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)
2465  if( inucl_w ) then
2466  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)
2467  else
2468  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 )
2469  endif
2470  pq(i_liccn,k,i,j) = min(dli_max, pq(i_niccn,k,i,j)*xi_ccn )
2471  ! only for output
2472 ! dni_ratio(k,i,j) = dssidt_rad(k,i,j)/( w_dssidz(k,i,j)+dssidt_rad(k,i,j) )
2473  else
2474  pq(i_niccn,k,i,j) = 0.0_rp
2475  pq(i_liccn,k,i,j) = 0.0_rp
2476  end if
2477  end do
2478  end do
2479  end do
2480 
2481  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_atmos_saturation
module ATMOSPHERE / Saturation adjustment
Definition: scale_atmos_sub_saturation.F90:16
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:234
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_atmos_saturation::atmos_saturation_dqsi_dtem_rho
subroutine, public atmos_saturation_dqsi_dtem_rho(dqsdtem, temp, dens)
Definition: scale_atmos_sub_saturation.F90:1285
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ice_multiplication_kij()

subroutine scale_atmos_phy_mp_sn14::ice_multiplication_kij ( real(rp), dimension(pq_max,ka,ia,ja), intent(out)  PQ,
real(rp), dimension(pac_max,ka,ia,ja), intent(in)  Pac,
real(rp), dimension(ka,ia,ja), intent(in)  tem,
real(rp), dimension(qa,ka,ia,ja), intent(in)  rhoq,
real(rp), dimension(5,ka,ia,ja), intent(in)  xq 
)

Definition at line 2488 of file scale_atmos_phy_mp_sn14.F90.

References scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, dc_log::log(), and scale_specfunc::sf_gamma().

Referenced by atmos_phy_mp_sn14().

2488 
2489  ! ice multiplication by splintering
2490  ! we consider Hallet-Mossop process
2491  use scale_specfunc, only: &
2492  gammafunc => sf_gamma
2493  implicit none
2494  real(RP), intent(out):: pq(pq_max,ka,ia,ja)
2495  !
2496  real(RP), intent(in) :: pac(pac_max,ka,ia,ja)
2497  real(RP), intent(in) :: tem(ka,ia,ja)
2498  real(RP), intent(in) :: rhoq(qa,ka,ia,ja)
2499  real(RP), intent(in) :: xq(5,ka,ia,ja)
2500  !
2501  logical, save :: flag_first = .true.
2502  ! production of (350.d3 ice particles)/(cloud riming [g]) => 350*d6 [/kg]
2503  real(RP), parameter :: pice = 350.0e+6_rp
2504  ! production of (1 ice particle)/(250 cloud particles riming)
2505  real(RP), parameter :: pnc = 250.0_rp
2506  real(RP), parameter :: rc_cr= 12.e-6_rp ! critical size[micron]
2507  ! temperature factor
2508  real(RP) :: fp
2509  ! work for incomplete gamma function
2510  real(RP), save :: xc_cr ! mass[kg] of cloud with r=critical size[micron]
2511  real(RP), save :: alpha ! slope parameter of gamma function
2512  real(RP), save :: gm, lgm ! gamma(alpha), log(gamma(alpha))
2513  real(RP) :: igm ! in complete gamma(x,alpha)
2514  real(RP) :: x
2515  ! coefficient of expansion using in calculation of igm
2516  real(RP) :: a0,a1,a2,a3,a4,a5
2517  real(RP) :: a6,a7,a8,a9,a10
2518  real(RP) :: an1,an2,b0,b1,b2,c0,c1,c2
2519  real(RP) :: d0,d1,d2,e1,e2,h0,h1,h2
2520  real(RP), parameter :: eps=1.0e-30_rp
2521  ! number of cloud droplets larger than 12 micron(radius).
2522  real(RP) :: n12
2523  !
2524  real(RP) :: wn, wni, wns, wng
2525  integer :: i, j, k
2526  !
2527  if( flag_first )then
2528  flag_first = .false.
2529  ! work for Incomplete Gamma function
2530  xc_cr = (2.0_rp*rc_cr/a_m(i_qc))**(1.0_rp/b_m(i_qc))
2531  alpha = (nu(i_qc)+1.0_rp)/mu(i_qc)
2532  gm = gammafunc(alpha)
2533  lgm = log(gm)
2534  end if
2535  !
2536  do j = js, je
2537  do i = is, ie
2538  do k = ks, ke
2539  ! Here we assume particle temperature is same as environment temperature.
2540  ! If you want to treat in a better manner,
2541  ! you can diagnose with eq.(64) in CT(86)
2542  if (tem(k,i,j) > 270.16_rp)then
2543  fp = 0.0_rp
2544  else if(tem(k,i,j) >= 268.16_rp)then
2545  fp = (270.16_rp-tem(k,i,j))*0.5_rp
2546  else if(tem(k,i,j) >= 265.16_rp)then
2547  fp = (tem(k,i,j)-265.16_rp)*0.333333333_rp
2548  else
2549  fp = 0.0_rp
2550  end if
2551  ! Approximation of Incomplete Gamma function
2552  ! Here we calculate with algorithm by Numerical Recipes.
2553  ! This approach is based on eq.(78) in Cotton etal.(1986),
2554  ! but more accurate and expanded for Generalized Gamma Distribution.
2555  x = coef_lambda(i_qc)*(xc_cr/xq(i_c,k,i,j))**mu(i_qc)
2556  !
2557  if(x<1.e-2_rp*alpha)then ! negligible
2558  igm = 0.0_rp
2559  else if(x<alpha+1.0_rp)then ! series expansion
2560  ! 10th-truncation is enough for cloud droplet.
2561  a0 = 1.0_rp/alpha ! n=0
2562  a1 = a0*x/(alpha+1.0_rp) ! n=1
2563  a2 = a1*x/(alpha+2.0_rp) ! n=2
2564  a3 = a2*x/(alpha+3.0_rp) ! n=3
2565  a4 = a3*x/(alpha+4.0_rp) ! n=4
2566  a5 = a4*x/(alpha+5.0_rp) ! n=5
2567  a6 = a5*x/(alpha+6.0_rp) ! n=6
2568  a7 = a6*x/(alpha+7.0_rp) ! n=7
2569  a8 = a7*x/(alpha+8.0_rp) ! n=8
2570  a9 = a8*x/(alpha+9.0_rp) ! n=9
2571  a10 = a9*x/(alpha+10.0_rp) ! n=10
2572  igm = (a0+a1+a2+a3+a4+a5+a6+a7+a8+a9+a10)*exp( -x + alpha*log(x) - lgm )
2573  else if(x<alpha*1.d2) then ! continued fraction expansion
2574  ! 2nd-truncation is enough for cloud droplet.
2575  ! setup
2576  b0 = x+1.0_rp-alpha
2577  c0 = 1.0_rp/eps
2578  d0 = 1.0_rp/b0
2579  h0 = d0
2580  ! n=1
2581  an1 = -(1.0_rp-alpha)
2582  b1 = b0 + 2.0_rp
2583  d1 = 1.0_rp/(an1*d0+b1)
2584  c1 = b1+an1/c0
2585  e1 = d1*c1
2586  h1 = h0*e1
2587  ! n=2
2588  an2 = -2.0_rp*(2.0_rp-alpha)
2589  b2 = b1 + 2.0_rp
2590  d2 = 1.0_rp/(an2*d1+b2)
2591  c2 = b2+an2/c1
2592  e2 = d2*c2
2593  h2 = h1*e2
2594  !
2595  igm = 1.0_rp - exp( -x + alpha*log(x) - lgm )*h2
2596  else ! negligible
2597  igm = 1.0_rp
2598  end if
2599  ! n12 is number of cloud droplets larger than 12 micron.
2600  n12 = rhoq(i_nc,k,i,j)*(1.0_rp-igm)
2601  ! eq.(82) CT(86)
2602  wn = (pice + n12/((rhoq(i_qc,k,i,j)+xc_min)*pnc) )*fp ! filtered by xc_min
2603  wni = wn*(-pac(i_liaclc2li,k,i,j) ) ! riming production rate is all negative
2604  wns = wn*(-pac(i_lsaclc2ls,k,i,j) )
2605  wng = wn*(-pac(i_lgaclc2lg,k,i,j) )
2606  pq(i_nispl,k,i,j) = wni+wns+wng
2607  !
2608  pq(i_lsspl,k,i,j) = - wns*xq(i_i,k,i,j) ! snow => ice
2609  pq(i_lgspl,k,i,j) = - wng*xq(i_i,k,i,j) ! graupel => ice
2610  !
2611  end do
2612  end do
2613  end do
2614  !
2615  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_specfunc::sf_gamma
real(rp) function, public sf_gamma(x)
Gamma function.
Definition: scale_specfunc.F90:48
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_specfunc
module SPECFUNC
Definition: scale_specfunc.F90:13
dc_log::log
subroutine, public log(type, message)
Definition: dc_log.f90:133
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the call graph for this function:
Here is the caller graph for this function:

◆ mixed_phase_collection_kij()

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

Definition at line 2623 of file scale_atmos_phy_mp_sn14.F90.

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

Referenced by atmos_phy_mp_sn14().

2623 ! rho ) ! [Add] 11/08/30
2624  use scale_atmos_saturation, only: &
2625  moist_psat_ice => atmos_saturation_psat_ice
2626  implicit none
2627 
2628  !--- mixed-phase collection process
2629  ! And all we set all production term as a negative sign to avoid confusion.
2630  !
2631  real(RP), intent(out):: pac(pac_max,ka,ia,ja)
2632  !--- partial conversion
2633  real(RP), intent(out):: pq(pq_max,ka,ia,ja)
2634  !
2635  real(RP), intent(in) :: wtem(ka,ia,ja)
2636  !--- mass/number concentration[kg/m3]
2637  real(RP), intent(in) :: rhoq(qa,ka,ia,ja)
2638  ! necessary ?
2639  real(RP), intent(in) :: xq(5,ka,ia,ja)
2640  !--- diameter of averaged mass( D(ave x) )
2641  real(RP), intent(in) :: dq_xave(5,ka,ia,ja)
2642  !--- terminal velocity of averaged mass( vt(ave x) )
2643  real(RP), intent(in) :: vt_xave(5,2,ka,ia,ja)
2644  ! [Add] 11/08/30 T.Mitsui, for autoconversion of ice
2645 ! real(RP), intent(in) :: rho(KA,IA,JA)
2646  !
2647  ! namelist variables
2648  !=== for collection
2649  !--- threshold of diameters to collide with others
2650  real(RP), save :: dc0 = 15.0e-6_rp ! lower threshold of cloud
2651  real(RP), save :: dc1 = 40.0e-6_rp ! upper threshold of cloud
2652  real(RP), save :: di0 = 150.0e-6_rp ! lower threshold of cloud
2653  real(RP), save :: ds0 = 150.0e-6_rp ! lower threshold of cloud
2654  real(RP), save :: dg0 = 150.0e-6_rp ! lower threshold of cloud
2655  !--- standard deviation of terminal velocity[m/s]
2656  real(RP), save :: sigma_c=0.00_rp ! cloud
2657  real(RP), save :: sigma_r=0.00_rp ! rain
2658  real(RP), save :: sigma_i=0.2_rp ! ice
2659  real(RP), save :: sigma_s=0.2_rp ! snow
2660  real(RP), save :: sigma_g=0.00_rp ! graupel
2661  !--- max collection efficiency for cloud
2662  real(RP), save :: e_im = 0.80_rp ! ice max
2663  real(RP), save :: e_sm = 0.80_rp ! snow max
2664  real(RP), save :: e_gm = 1.00_rp ! graupel max
2665  !--- collection efficiency between 2 species
2666  real(RP), save :: e_ir=1.0_rp ! ice x rain
2667  real(RP), save :: e_sr=1.0_rp ! snow x rain
2668  real(RP), save :: e_gr=1.0_rp ! graupel x rain
2669  real(RP), save :: e_ii=1.0_rp ! ice x ice
2670  real(RP), save :: e_si=1.0_rp ! snow x ice
2671  real(RP), save :: e_gi=1.0_rp ! graupel x ice
2672  real(RP), save :: e_ss=1.0_rp ! snow x snow
2673  real(RP), save :: e_gs=1.0_rp ! graupel x snow
2674  real(RP), save :: e_gg=1.0_rp ! graupel x graupel
2675  !=== for partial conversion
2676  !--- flag: 1=> partial conversion to graupel, 0=> no conversion
2677  integer, save :: i_iconv2g=1 ! ice => graupel
2678  integer, save :: i_sconv2g=1 ! snow => graupel
2679  !--- bulk density of graupel
2680  real(RP), save :: rho_g = 900.0_rp ! [kg/m3]
2681  !--- space filling coefficient [%]
2682  real(RP), save :: cfill_i = 0.68_rp ! ice
2683  real(RP), save :: cfill_s = 0.01_rp ! snow
2684  !--- critical diameter for ice conversion
2685  real(RP), save :: di_cri = 500.e-6_rp ! [m]
2686  ! [Add] 10/08/03 T.Mitsui
2687  logical, save :: opt_stick_ks96=.false.
2688  logical, save :: opt_stick_co86=.false.
2689  real(RP), parameter :: a_dec = 0.883_rp
2690  real(RP), parameter :: b_dec = 0.093_rp
2691  real(RP), parameter :: c_dec = 0.00348_rp
2692  real(RP), parameter :: d_dec = 4.5185e-5_rp
2693  !
2694  logical, save :: flag_first = .true.
2695  namelist /nm_mp_sn14_collection/ &
2696  dc0, dc1, di0, ds0, dg0, &
2697  sigma_c, sigma_r, sigma_i, sigma_s, sigma_g, &
2698  opt_stick_ks96, &
2699  opt_stick_co86, &
2700  e_im, e_sm, e_gm, &
2701  e_ir, e_sr, e_gr, e_ii, e_si, e_gi, e_ss, e_gs, e_gg, &
2702  i_iconv2g, i_sconv2g, rho_g, cfill_i, cfill_s, di_cri
2703  !
2704  real(RP) :: tem(ka,ia,ja)
2705  !
2706  !--- collection efficency of each specie
2707  real(RP) :: e_c, e_r, e_i, e_s, e_g !
2708  real(RP) :: e_ic, e_sc, e_gc !
2709  !--- sticking efficiency
2710  real(RP) :: e_stick(ka,ia,ja)
2711  ! [Add] 10/08/03 T.Mitsui
2712  real(RP) :: temc, temc2, temc3
2713  real(RP) :: e_dec
2714  real(RP) :: esi_rat
2715  real(RP) :: esi(ka,ia,ja)
2716  !
2717  real(RP) :: temc_p, temc_m ! celcius tem.
2718  ! [Add] 11/08/30 T.Mitsui, estimation of autoconversion time
2719 ! real(RP) :: ci_aut(KA,IA,JA)
2720 ! real(RP) :: taui_aut(KA,IA,JA)
2721 ! real(RP) :: tau_sce(KA,IA,JA)
2722  !--- DSD averaged diameter for each species
2723  real(RP) :: ave_dc ! cloud
2724 ! real(RP) :: ave_dr ! rain
2725  real(RP) :: ave_di ! ice
2726  real(RP) :: ave_ds ! snow
2727  real(RP) :: ave_dg ! graupel
2728  !--- coefficient of collection equations(L:mass, N:number)
2729  real(RP) :: coef_acc_lci, coef_acc_nci ! cloud - cloud ice
2730  real(RP) :: coef_acc_lcs, coef_acc_ncs ! cloud - snow
2731  !
2732  real(RP) :: coef_acc_lcg, coef_acc_ncg ! cloud - graupel
2733  real(RP) :: coef_acc_lri_i, coef_acc_nri_i ! rain - cloud ice
2734  real(RP) :: coef_acc_lri_r, coef_acc_nri_r ! rain - cloud ice
2735  real(RP) :: coef_acc_lrs_s, coef_acc_nrs_s ! rain - snow
2736  real(RP) :: coef_acc_lrs_r, coef_acc_nrs_r ! rain - snow
2737  real(RP) :: coef_acc_lrg, coef_acc_nrg ! rain - graupel
2738  real(RP) :: coef_acc_lii, coef_acc_nii ! cloud ice - cloud ice
2739  real(RP) :: coef_acc_lis, coef_acc_nis ! cloud ice - snow
2740  real(RP) :: coef_acc_nss ! snow - snow
2741  real(RP) :: coef_acc_ngg ! grauepl - graupel
2742  real(RP) :: coef_acc_lsg, coef_acc_nsg ! snow - graupel
2743  !--- (diameter) x (diameter)
2744  real(RP) :: dcdc, dcdi, dcds, dcdg
2745  real(RP) :: drdr, drdi, drds, drdg
2746  real(RP) :: didi, dids, didg
2747  real(RP) :: dsds, dsdg
2748  real(RP) :: dgdg
2749  !--- (terminal velocity) x (terminal velocity)
2750  real(RP) :: vcvc, vcvi, vcvs, vcvg
2751  real(RP) :: vrvr, vrvi, vrvs, vrvg
2752  real(RP) :: vivi, vivs, vivg
2753  real(RP) :: vsvs, vsvg
2754  real(RP) :: vgvg
2755  !
2756  real(RP) :: wx_cri, wx_crs
2757  real(RP) :: coef_emelt
2758  real(RP) :: w1
2759 
2760  real(RP) :: sw
2761  !
2762  integer :: i, j, k
2763  !
2764  if( flag_first )then
2765  rewind( io_fid_conf )
2766  read( io_fid_conf, nml=nm_mp_sn14_collection, end=100 )
2767 100 if( io_l ) write( io_fid_log, nml=nm_mp_sn14_collection )
2768  flag_first = .false.
2769  end if
2770  !
2771  ! [Add] 10/08/03 T.Mitsui
2772  do j = js, je
2773  do i = is, ie
2774  do k = ks, ke
2775  tem(k,i,j) = max( wtem(k,i,j), tem_min ) ! 11/08/30 T.Mitsui
2776  end do
2777  end do
2778  end do
2779 
2780  call moist_psat_ice( esi, tem )
2781 
2782  if( opt_stick_ks96 )then
2783  do j = js, je
2784  do i = is, ie
2785  do k = ks, ke
2786  ! Khain and Sednev (1996), eq.(3.15)
2787  temc = tem(k,i,j) - t00
2788  temc2 = temc*temc
2789  temc3 = temc2*temc
2790  e_dec = max(0.0_rp, a_dec + b_dec*temc + c_dec*temc2 + d_dec*temc3 )
2791  esi_rat = rhoq(i_qv,k,i,j)*rvap*tem(k,i,j)/esi(k,i,j)
2792  e_stick(k,i,j) = min(1.0_rp, e_dec*esi_rat)
2793  end do
2794  end do
2795  end do
2796  else if( opt_stick_co86 )then
2797  do j = js, je
2798  do i = is, ie
2799  do k = ks, ke
2800  ! [Add] 11/08/30 T.Mitsui, Cotton et al. (1986)
2801  temc = min(tem(k,i,j) - t00,0.0_rp)
2802  w1 = 0.035_rp*temc-0.7_rp
2803  e_stick(k,i,j) = 10._rp**w1
2804  end do
2805  end do
2806  end do
2807  else
2808  do j = js, je
2809  do i = is, ie
2810  do k = ks, ke
2811  ! Lin et al. (1983)
2812  temc_m = min(tem(k,i,j) - t00,0.0_rp) ! T < 273.15
2813  e_stick(k,i,j) = exp(0.09_rp*temc_m)
2814  end do
2815  end do
2816  end do
2817  end if
2818  !
2819  profile_start("sn14_collection")
2820 !OCL NOSIMD
2821  do j = js, je
2822  do i = is, ie
2823  do k = ks, ke
2824  !
2825 ! temc_m = min(tem(k,i,j) - T00,0.0_RP) ! T < 273.15
2826  temc_p = max(tem(k,i,j) - t00,0.0_rp) ! T > 273.15
2827  ! averaged diameter using SB06(82)
2828  ave_dc = coef_d(i_qc)*xq(i_c,k,i,j)**b_m(i_qc)
2829  ave_di = coef_d(i_qi)*xq(i_i,k,i,j)**b_m(i_qi)
2830  ave_ds = coef_d(i_qs)*xq(i_s,k,i,j)**b_m(i_qs)
2831  ave_dg = coef_d(i_qg)*xq(i_g,k,i,j)**b_m(i_qg)
2832  !------------------------------------------------------------------------
2833  ! coellection efficiency are given as follows
2834  e_c = max(0.0_rp, min(1.0_rp, (ave_dc-dc0)/(dc1-dc0) ))
2835  sw = 0.5_rp - sign(0.5_rp, di0-ave_di) ! if(ave_di>di0)then sw=1
2836  e_i = e_im * sw
2837  sw = 0.5_rp - sign(0.5_rp, ds0-ave_ds) ! if(ave_ds>ds0)then sw=1
2838  e_s = e_sm * sw
2839  sw = 0.5_rp - sign(0.5_rp, dg0-ave_dg) ! if(ave_dg>dg0)then sw=1
2840  e_g = e_gm * sw
2841  e_ic = e_i*e_c
2842  e_sc = e_s*e_c
2843  e_gc = e_g*e_c
2844  !------------------------------------------------------------------------
2845  ! Collection: a collects b ( assuming particle size a>b )
2846  dcdc = dq_xave(i_c,k,i,j) * dq_xave(i_c,k,i,j)
2847  drdr = dq_xave(i_r,k,i,j) * dq_xave(i_r,k,i,j)
2848  didi = dq_xave(i_i,k,i,j) * dq_xave(i_i,k,i,j)
2849  dsds = dq_xave(i_s,k,i,j) * dq_xave(i_s,k,i,j)
2850  dgdg = dq_xave(i_g,k,i,j) * dq_xave(i_g,k,i,j)
2851  dcdi = dq_xave(i_c,k,i,j) * dq_xave(i_i,k,i,j)
2852  dcds = dq_xave(i_c,k,i,j) * dq_xave(i_s,k,i,j)
2853  dcdg = dq_xave(i_c,k,i,j) * dq_xave(i_g,k,i,j)
2854  drdi = dq_xave(i_r,k,i,j) * dq_xave(i_i,k,i,j)
2855  drds = dq_xave(i_r,k,i,j) * dq_xave(i_s,k,i,j)
2856  drdg = dq_xave(i_r,k,i,j) * dq_xave(i_g,k,i,j)
2857  dids = dq_xave(i_i,k,i,j) * dq_xave(i_s,k,i,j)
2858 ! didg = dq_xave(I_I,k,i,j) * dq_xave(I_G,k,i,j)
2859  dsdg = dq_xave(i_s,k,i,j) * dq_xave(i_g,k,i,j)
2860  !
2861  vcvc = vt_xave(i_c,2,k,i,j)* vt_xave(i_c,2,k,i,j)
2862  vrvr = vt_xave(i_r,2,k,i,j)* vt_xave(i_r,2,k,i,j)
2863  vivi = vt_xave(i_i,2,k,i,j)* vt_xave(i_i,2,k,i,j)
2864  vsvs = vt_xave(i_s,2,k,i,j)* vt_xave(i_s,2,k,i,j)
2865  vgvg = vt_xave(i_g,2,k,i,j)* vt_xave(i_g,2,k,i,j)
2866  vcvi = vt_xave(i_c,2,k,i,j)* vt_xave(i_i,2,k,i,j)
2867  vcvs = vt_xave(i_c,2,k,i,j)* vt_xave(i_s,2,k,i,j)
2868  vcvg = vt_xave(i_c,2,k,i,j)* vt_xave(i_g,2,k,i,j)
2869  vrvi = vt_xave(i_r,2,k,i,j)* vt_xave(i_i,2,k,i,j)
2870  vrvs = vt_xave(i_r,2,k,i,j)* vt_xave(i_s,2,k,i,j)
2871  vrvg = vt_xave(i_r,2,k,i,j)* vt_xave(i_g,2,k,i,j)
2872  vivs = vt_xave(i_i,2,k,i,j)* vt_xave(i_s,2,k,i,j)
2873 ! vivg = vt_xave(I_I,2,k,i,j)* vt_xave(I_G,2,k,i,j)
2874  vsvg = vt_xave(i_s,2,k,i,j)* vt_xave(i_g,2,k,i,j)
2875  !------------------------------------------------------------------------
2876  !
2877  !+++ pattern 1: a + b => a (a>b)
2878  ! (i-c, s-c, g-c, s-i, g-r, s-g)
2879  !------------------------------------------------------------------------
2880  ! cloud-ice => ice
2881  ! reduction term of cloud
2882  coef_acc_lci = &
2883  ( delta_b1(i_qc)*dcdc + delta_ab1(i_qi,i_qc)*dcdi + delta_b0(i_qi)*didi ) &
2884  * ( theta_b1(i_qc)*vcvc - theta_ab1(i_qi,i_qc)*vcvi + theta_b0(i_qi)*vivi &
2885  + sigma_i + sigma_c )**0.5_rp
2886  coef_acc_nci = &
2887  ( delta_b0(i_qc)*dcdc + delta_ab0(i_qi,i_qc)*dcdi + delta_b0(i_qi)*didi ) &
2888  * ( theta_b0(i_qc)*vcvc - theta_ab0(i_qi,i_qc)*vcvi + theta_b0(i_qi)*vivi &
2889  + sigma_i + sigma_c )**0.5_rp
2890  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
2891  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
2892  ! cloud-snow => snow
2893  ! reduction term of cloud
2894  coef_acc_lcs = &
2895  ( delta_b1(i_qc)*dcdc + delta_ab1(i_qs,i_qc)*dcds + delta_b0(i_qs)*dsds ) &
2896  * ( theta_b1(i_qc)*vcvc - theta_ab1(i_qs,i_qc)*vcvs + theta_b0(i_qs)*vsvs &
2897  + sigma_s + sigma_c )**0.5_rp
2898  coef_acc_ncs = &
2899  ( delta_b0(i_qc)*dcdc + delta_ab0(i_qs,i_qc)*dcds + delta_b0(i_qs)*dsds ) &
2900  * ( theta_b0(i_qc)*vcvc - theta_ab0(i_qs,i_qc)*vcvs + theta_b0(i_qs)*vsvs &
2901  + sigma_s + sigma_c )**0.5_rp
2902  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
2903  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
2904  ! cloud-graupel => graupel
2905  ! reduction term of cloud
2906  coef_acc_lcg = &
2907  ( delta_b1(i_qc)*dcdc + delta_ab1(i_qg,i_qc)*dcdg + delta_b0(i_qg)*dgdg ) &
2908  * ( theta_b1(i_qc)*vcvc - theta_ab1(i_qg,i_qc)*vcvg + theta_b0(i_qg)*vgvg &
2909  + sigma_g + sigma_c )**0.5_rp
2910  coef_acc_ncg = &
2911  ( delta_b0(i_qc)*dcdc + delta_ab0(i_qg,i_qc)*dcdg + delta_b0(i_qg)*dgdg ) &
2912  * ( theta_b0(i_qc)*vcvc - theta_ab0(i_qg,i_qc)*vcvg + theta_b0(i_qg)*vgvg &
2913  + sigma_g + sigma_c )**0.5_rp
2914  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
2915  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
2916  ! snow-graupel => graupel
2917  coef_acc_lsg = &
2918  ( delta_b1(i_qs)*dsds + delta_ab1(i_qg,i_qs)*dsdg + delta_b0(i_qg)*dgdg ) &
2919  * ( theta_b1(i_qs)*vsvs - theta_ab1(i_qg,i_qs)*vsvg + theta_b0(i_qg)*vgvg &
2920  + sigma_g + sigma_s )**0.5_rp
2921  coef_acc_nsg = &
2922  ( delta_b0(i_qs)*dsds + delta_ab0(i_qg,i_qs)*dsdg + delta_b0(i_qg)*dgdg ) &
2923  ! [fix] T.Mitsui 08/05/08
2924  * ( theta_b0(i_qs)*vsvs - theta_ab0(i_qg,i_qs)*vsvg + theta_b0(i_qg)*vgvg &
2925  + sigma_g + sigma_s )**0.5_rp
2926  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
2927  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
2928  !------------------------------------------------------------------------
2929  ! ice-snow => snow
2930  ! reduction term of ice
2931  coef_acc_lis = &
2932  ( delta_b1(i_qi)*didi + delta_ab1(i_qs,i_qi)*dids + delta_b0(i_qs)*dsds ) &
2933  * ( theta_b1(i_qi)*vivi - theta_ab1(i_qs,i_qi)*vivs + theta_b0(i_qs)*vsvs &
2934  + sigma_i + sigma_s )**0.5_rp
2935  coef_acc_nis = &
2936  ( delta_b0(i_qi)*didi + delta_ab0(i_qs,i_qi)*dids + delta_b0(i_qs)*dsds ) &
2937  * ( theta_b0(i_qi)*vivi - theta_ab0(i_qs,i_qi)*vivs + theta_b0(i_qs)*vsvs &
2938  + sigma_i + sigma_s )**0.5_rp
2939  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
2940  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
2941  !
2942  sw = sign(0.5_rp, t00-tem(k,i,j)) + 0.5_rp
2943  ! if ( tem(k,i,j) <= T00 )then
2944  ! rain-graupel => graupel
2945  ! reduction term of rain
2946  ! sw = 1
2947  ! else
2948  ! rain-graupel => rain
2949  ! reduction term of graupel
2950  ! sw = 0
2951  coef_acc_lrg = &
2952  ( ( delta_b1(i_qr)*drdr + delta_ab1(i_qg,i_qr)*drdg + delta_b0(i_qg)*dgdg ) * sw &
2953  + ( delta_b1(i_qg)*dgdg + delta_ab1(i_qr,i_qg)*drdg + delta_b0(i_qr)*drdr ) * (1.0_rp-sw) ) &
2954  * sqrt( ( theta_b1(i_qr)*vrvr - theta_ab1(i_qg,i_qr)*vrvg + theta_b0(i_qg)*vgvg ) * sw &
2955  + ( theta_b1(i_qg)*vgvg - theta_ab1(i_qr,i_qg)*vrvg + theta_b0(i_qr)*vrvr ) * (1.0_rp-sw) &
2956  + sigma_r + sigma_g )
2957  pac(i_lraclg2lg,k,i,j) = -0.25_rp*pi*e_gr*coef_acc_lrg &
2958  * ( rhoq(i_ng,k,i,j)*rhoq(i_qr,k,i,j) * sw &
2959  + rhoq(i_nr,k,i,j)*rhoq(i_qg,k,i,j) * (1.0_rp-sw) )
2960  coef_acc_nrg = &
2961  ( delta_b0(i_qr)*drdr + delta_ab0(i_qg,i_qr)*drdg + delta_b0(i_qg)*dgdg ) &
2962  * ( theta_b0(i_qr)*vrvr - theta_ab0(i_qg,i_qr)*vrvg + theta_b0(i_qg)*vgvg &
2963  + sigma_r + sigma_g )**0.5_rp
2964  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
2965  !
2966  !------------------------------------------------------------------------
2967  !
2968  !+++ pattern 2: a + b => c (a>b)
2969  ! (r-i,r-s)
2970  !------------------------------------------------------------------------
2971  ! rain-ice => graupel
2972  ! reduction term of ice
2973  coef_acc_lri_i = &
2974  ( delta_b1(i_qi)*didi + delta_ab1(i_qr,i_qi)*drdi + delta_b0(i_qr)*drdr ) &
2975  * ( theta_b1(i_qi)*vivi - theta_ab1(i_qr,i_qi)*vrvi + theta_b0(i_qr)*vrvr &
2976  + sigma_r + sigma_i )**0.5_rp
2977  coef_acc_nri_i = &
2978  ( delta_b0(i_qi)*didi + delta_ab0(i_qr,i_qi)*drdi + delta_b0(i_qr)*drdr ) &
2979  * ( theta_b0(i_qi)*vivi - theta_ab0(i_qr,i_qi)*vrvi + theta_b0(i_qr)*vrvr &
2980  + sigma_r + sigma_i )**0.5_rp
2981  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
2982  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
2983  ! reduction term of rain
2984  coef_acc_lri_r = &
2985  ( delta_b1(i_qr)*drdr + delta_ab1(i_qi,i_qr)*drdi + delta_b0(i_qi)*didi ) &
2986  * ( theta_b1(i_qr)*vrvr - theta_ab1(i_qi,i_qr)*vrvi + theta_b0(i_qi)*vivi &
2987  + sigma_r + sigma_i )**0.5_rp
2988  coef_acc_nri_r = &
2989  ( delta_b0(i_qr)*drdr + delta_ab0(i_qi,i_qr)*drdi + delta_b0(i_qi)*didi ) &
2990  * ( theta_b0(i_qr)*vrvr - theta_ab0(i_qi,i_qr)*vrvi + theta_b0(i_qi)*vivi &
2991  + sigma_r + sigma_i )**0.5_rp
2992  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
2993  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
2994  ! rain-snow => graupel
2995  ! reduction term of snow
2996  coef_acc_lrs_s = &
2997  ( delta_b1(i_qs)*dsds + delta_ab1(i_qr,i_qs)*drds + delta_b0(i_qr)*drdr ) &
2998  * ( theta_b1(i_qs)*vsvs - theta_ab1(i_qr,i_qs)*vrvs + theta_b0(i_qr)*vrvr &
2999  + sigma_r + sigma_s )**0.5_rp
3000  coef_acc_nrs_s = &
3001  ( delta_b0(i_qs)*dsds + delta_ab0(i_qr,i_qs)*drds + delta_b0(i_qr)*drdr ) &
3002  * ( theta_b0(i_qs)*vsvs - theta_ab0(i_qr,i_qs)*vrvs + theta_b0(i_qr)*vrvr &
3003  + sigma_r + sigma_s )**0.5_rp
3004  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
3005  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
3006  ! reduction term of rain
3007  coef_acc_lrs_r = &
3008  ( delta_b1(i_qr)*drdr + delta_ab1(i_qs,i_qr)*drds + delta_b0(i_qs)*dsds ) &
3009  * ( theta_b1(i_qr)*vrvr - theta_ab1(i_qs,i_qr)*vrvs + theta_b0(i_qs)*vsvs &
3010  + sigma_r + sigma_s )**0.5_rp
3011  coef_acc_nrs_r = &
3012  ( delta_b0(i_qr)*drdr + delta_ab0(i_qs,i_qr)*drds + delta_b0(i_qs)*dsds ) &
3013  * ( theta_b0(i_qr)*vrvr - theta_ab0(i_qs,i_qr)*vrvs + theta_b0(i_qs)*vsvs &
3014  + sigma_r + sigma_s )**0.5_rp
3015  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
3016  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
3017  !------------------------------------------------------------------------
3018  !
3019  !+++ pattern 3: a + a => b (i-i)
3020  !
3021  !------------------------------------------------------------------------
3022  ! ice-ice ( reduction is double, but includes double-count)
3023  coef_acc_lii = &
3024  ( delta_b0(i_qi)*didi + delta_ab1(i_qi,i_qi)*didi + delta_b1(i_qi)*didi ) &
3025  * ( theta_b0(i_qi)*vivi - theta_ab1(i_qi,i_qi)*vivi + theta_b1(i_qi)*vivi &
3026  + sigma_i + sigma_i )**0.5_rp
3027  coef_acc_nii = &
3028  ( delta_b0(i_qi)*didi + delta_ab0(i_qi,i_qi)*didi + delta_b0(i_qi)*didi ) &
3029  * ( theta_b0(i_qi)*vivi - theta_ab0(i_qi,i_qi)*vivi + theta_b0(i_qi)*vivi &
3030  + sigma_i + sigma_i )**0.5_rp
3031  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
3032  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
3033  !
3034 ! ci_aut(k,i,j) = 0.25_RP*pi*E_ii*rhoq(I_NI,k,i,j)*coef_acc_LII
3035 ! taui_aut(k,i,j) = 1._RP/max(E_stick(k,i,j)*ci_aut(k,i,j),1.E-10_RP)
3036 ! 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)
3037  !------------------------------------------------------------------------
3038  !
3039  !+++ pattern 4: a + a => a (s-s)
3040  !
3041  !------------------------------------------------------------------------
3042  ! snow-snow => snow
3043  coef_acc_nss = &
3044  ( delta_b0(i_qs)*dsds + delta_ab0(i_qs,i_qs)*dsds + delta_b0(i_qs)*dsds ) &
3045  * ( theta_b0(i_qs)*vsvs - theta_ab0(i_qs,i_qs)*vsvs + theta_b0(i_qs)*vsvs &
3046  + sigma_s + sigma_s )**0.5_rp
3047  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
3048  !
3049  ! graupel-grauple => graupel
3050  coef_acc_ngg = &
3051  ( delta_b0(i_qg)*dgdg + delta_ab0(i_qg,i_qg)*dgdg + delta_b0(i_qg)*dgdg ) &
3052  * ( theta_b0(i_qg)*vgvg - theta_ab0(i_qg,i_qg)*vgvg + theta_b0(i_qg)*vgvg &
3053  + sigma_g + sigma_g )**0.5_rp
3054  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
3055  !
3056  !------------------------------------------------------------------------
3057  !--- Partial conversion
3058  ! SB06(70),(71)
3059  ! i_iconv2g: option whether partial conversions work or not
3060  ! ice-cloud => graupel
3061  sw = 0.5_rp - sign(0.5_rp,di_cri-ave_di) ! if( ave_di > di_cri )then sw=1
3062  wx_cri = cfill_i*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_di*ave_di*ave_di/xq(i_i,k,i,j) - 1.0_rp ) * sw
3063  pq(i_licon,k,i,j) = i_iconv2g * pac(i_liaclc2li,k,i,j)/max(1.0_rp, wx_cri) * sw
3064  pq(i_nicon,k,i,j) = i_iconv2g * pq(i_licon,k,i,j)/xq(i_i,k,i,j) * sw
3065 
3066  ! snow-cloud => graupel
3067  wx_crs = cfill_s*dwatr/rho_g*( pi/6.0_rp*rho_g*ave_ds*ave_ds*ave_ds/xq(i_s,k,i,j) - 1.0_rp )
3068  pq(i_lscon,k,i,j) = i_sconv2g * (pac(i_lsaclc2ls,k,i,j))/max(1.0_rp, wx_crs)
3069  pq(i_nscon,k,i,j) = i_sconv2g * pq(i_lscon,k,i,j)/xq(i_s,k,i,j)
3070  !------------------------------------------------------------------------
3071  !--- enhanced melting( due to collection-freezing of water droplets )
3072  ! originally from Rutledge and Hobbs(1984). eq.(A.21)
3073  ! if T > 273.15 then temc_p=T-273.15, else temc_p=0
3074  ! 08/05/08 [fix] T.Mitsui LHF00 => LHF0
3075  ! melting occurs around T=273K, so LHF0 is suitable both SIMPLE and EXACT,
3076  ! otherwise LHF can have sign both negative(EXACT) and positive(SIMPLE).
3077 !!$ coef_emelt = -CL/LHF00*temc_p
3078  coef_emelt = cl/lhf0*temc_p
3079  ! cloud-graupel
3080  pq(i_lgacm,k,i,j) = coef_emelt*pac(i_lgaclc2lg,k,i,j)
3081  pq(i_ngacm,k,i,j) = pq(i_lgacm,k,i,j)/xq(i_g,k,i,j)
3082  ! rain-graupel
3083  pq(i_lgarm,k,i,j) = coef_emelt*pac(i_lraclg2lg,k,i,j)
3084  pq(i_ngarm,k,i,j) = pq(i_lgarm,k,i,j)/xq(i_g,k,i,j)
3085  ! cloud-snow
3086  pq(i_lsacm,k,i,j) = coef_emelt*(pac(i_lsaclc2ls,k,i,j))
3087  pq(i_nsacm,k,i,j) = pq(i_lsacm,k,i,j)/xq(i_s,k,i,j)
3088  ! rain-snow
3089  pq(i_lsarm,k,i,j) = coef_emelt*(pac(i_lracls2lg_r,k,i,j)+pac(i_lracls2lg_s,k,i,j))
3090  pq(i_nsarm,k,i,j) = pq(i_lsarm,k,i,j)/xq(i_g,k,i,j) ! collect? might be I_S
3091  ! cloud-ice
3092  pq(i_liacm,k,i,j) = coef_emelt*pac(i_liaclc2li,k,i,j)
3093  pq(i_niacm,k,i,j) = pq(i_liacm,k,i,j)/xq(i_i,k,i,j)
3094  ! rain-ice
3095  pq(i_liarm,k,i,j) = coef_emelt*(pac(i_lracli2lg_r,k,i,j)+pac(i_lracli2lg_i,k,i,j))
3096  pq(i_niarm,k,i,j) = pq(i_liarm,k,i,j)/xq(i_g,k,i,j) ! collect? might be I_I
3097  end do
3098  end do
3099  end do
3100  profile_stop("sn14_collection")
3101 
3102  !
3103  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_atmos_saturation
module ATMOSPHERE / Saturation adjustment
Definition: scale_atmos_sub_saturation.F90:16
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the caller graph for this function:

◆ aut_acc_slc_brk_kij()

subroutine scale_atmos_phy_mp_sn14::aut_acc_slc_brk_kij ( real(rp), dimension(pq_max,ka,ia,ja), intent(out)  PQ,
real(rp), dimension(qa,ka,ia,ja), intent(in)  rhoq,
real(rp), dimension(5,ka,ia,ja), intent(in)  xq,
real(rp), dimension(5,ka,ia,ja), intent(in)  dq_xave,
real(rp), dimension(ka,ia,ja), intent(in)  rho 
)

Definition at line 3111 of file scale_atmos_phy_mp_sn14.F90.

References scale_const::const_eps, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, and scale_grid_index::ks.

Referenced by atmos_phy_mp_sn14().

3111  implicit none
3112 
3113  real(RP), intent(out) :: pq(pq_max,ka,ia,ja)
3114  !
3115  real(RP), intent(in) :: rhoq(qa,ka,ia,ja)
3116  real(RP), intent(in) :: xq(5,ka,ia,ja)
3117  real(RP), intent(in) :: dq_xave(5,ka,ia,ja)
3118  real(RP), intent(in) :: rho(ka,ia,ja)
3119  !
3120  ! parameter for autoconversion
3121  real(RP), parameter :: kcc = 4.44e+9_rp ! collision efficiency [m3/kg2/sec]
3122  real(RP), parameter :: tau_min = 1.e-20_rp ! empirical filter by T.Mitsui
3123  real(RP), parameter :: rx_sep = 1.0_rp/x_sep ! 1/x_sep, 10/08/03 [Add] T.Mitsui
3124  !
3125  ! parameter for accretion
3126  real(RP), parameter :: kcr = 5.8_rp ! collision efficiency [m3/kg2/sec]
3127  real(RP), parameter :: thr_acc = 5.e-5_rp ! threshold for universal function original
3128  !
3129  ! parameter for self collection and collison break-up
3130  real(RP), parameter :: krr = 4.33_rp ! k_rr, S08 (35)
3131  real(RP), parameter :: kaprr = 60.7_rp ! kappa_rr, SB06(11)
3132  real(RP), parameter :: kbr = 1000._rp ! k_br, SB06(14)
3133  real(RP), parameter :: kapbr = 2.3e+3_rp ! kappa_br, SB06(15)
3134  real(RP), parameter :: dr_min = 0.35e-3_rp ! minimum diameter, SB06(13)-(15)
3135  !
3136  ! work variables
3137  real(RP) :: coef_nuc0 ! coefficient of number for Auto-conversion
3138  real(RP) :: coef_nuc1 ! mass
3139  real(RP) :: coef_aut0 ! number
3140  real(RP) :: coef_aut1 ! mass
3141  real(RP) :: lwc ! lc+lr
3142  real(RP) :: tau ! conversion ratio: qr/(qc+qr) ranges [0:1]
3143  real(RP) :: rho_fac ! factor of air density
3144  real(RP) :: psi_aut ! Universal function of Auto-conversion
3145  real(RP) :: psi_acc ! Universal function of Accretion
3146  real(RP) :: psi_brk ! Universal function of Breakup
3147  real(RP) :: ddr ! diameter difference from equilibrium
3148  !
3149  integer :: i, j, k
3150  !
3151  coef_nuc0 = (nu(i_qc)+2.0_rp)/(nu(i_qc)+1.0_rp)
3152  coef_nuc1 = (nu(i_qc)+2.0_rp)*(nu(i_qc)+4.0_rp)/(nu(i_qc)+1.0_rp)/(nu(i_qc)+1.0_rp)
3153  coef_aut0 = -kcc*coef_nuc0
3154  coef_aut1 = -kcc/x_sep/20._rp*coef_nuc1
3155  !
3156  do j = js, je
3157  do i = is, ie
3158  do k = ks, ke
3159  lwc = rhoq(i_qr,k,i,j)+rhoq(i_qc,k,i,j)
3160  if( lwc > xc_min )then
3161  tau = max(tau_min, rhoq(i_qr,k,i,j)/lwc)
3162  else
3163  tau = tau_min
3164  end if
3165  rho_fac = sqrt(rho_0/max(rho(k,i,j),rho_min))
3166  !
3167  ! Auto-conversion ( cloud-cloud => rain )
3168  psi_aut = 400._rp*(tau**0.7_rp)*(1.0_rp - (tau**0.7_rp))**3 ! (6) SB06
3169  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
3170  ! lc = lwc*(1-tau), lr = lwc*tau
3171  pq(i_lcaut,k,i,j) = coef_aut1*lwc*lwc*xq(i_c,k,i,j)*xq(i_c,k,i,j) & ! (4) SB06
3172  *((1.0_rp-tau)*(1.0_rp-tau) + psi_aut)*rho_fac*rho_fac !
3173  pq(i_nraut,k,i,j) = -rx_sep*pq(i_lcaut,k,i,j) ! (A7) SB01
3174  !
3175  ! Accretion ( cloud-rain => rain )
3176  psi_acc =(tau/(tau+thr_acc))**4 ! (8) SB06
3177  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
3178  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
3179  !
3180  ! Self-collection ( rain-rain => rain )
3181  pq(i_nrslc,k,i,j) = -krr*rhoq(i_nr,k,i,j)*rhoq(i_qr,k,i,j)*rho_fac ! (A.8) SB01
3182  !
3183  ! Collisional breakup of rain
3184  ddr = min(1.e-3_rp, dq_xave(i_r,k,i,j) - dr_eq )
3185  if (dq_xave(i_r,k,i,j) < dr_min )then ! negligible
3186  psi_brk = -1.0_rp
3187  pq(i_nrbrk,k,i,j) = 0.0_rp
3188  else if (dq_xave(i_r,k,i,j) <= dr_eq )then
3189  psi_brk = kbr*ddr + 1.0_rp ! (14) SB06 (+1 is necessary)
3190  pq(i_nrbrk,k,i,j) = - (psi_brk + 1.0_rp)*pq(i_nrslc,k,i,j) ! (13) SB06
3191  else
3192  psi_brk = 2.0_rp*exp(kapbr*ddr) - 1.0_rp ! (15) SB06
3193  pq(i_nrbrk,k,i,j) = - (psi_brk + 1.0_rp)*pq(i_nrslc,k,i,j) ! (13) SB06
3194  end if
3195  !
3196  end do
3197  end do
3198  end do
3199  !
3200  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the caller graph for this function:

◆ freezing_water_kij()

subroutine scale_atmos_phy_mp_sn14::freezing_water_kij ( real(dp), intent(in)  dt,
real(rp), dimension(pq_max,ka,ia,ja), intent(out)  PQ,
real(rp), dimension(qa,ka,ia,ja), intent(in)  rhoq,
real(rp), dimension(5,ka,ia,ja), intent(in)  xq,
real(rp), dimension(ka,ia,ja), intent(in)  tem 
)

Definition at line 3446 of file scale_atmos_phy_mp_sn14.F90.

References scale_const::const_undef, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ka, scale_grid_index::ke, scale_grid_index::ks, and dc_log::log().

Referenced by atmos_phy_mp_sn14().

3446  !
3447  ! In this subroutine,
3448  ! We assumed surface temperature of droplets are same as environment.
3449  implicit none
3450 
3451  real(DP), intent(in) :: dt
3452  real(RP), intent(out):: pq(pq_max,ka,ia,ja)
3453  !
3454  real(RP), intent(in) :: tem(ka,ia,ja)
3455  !
3456  real(RP), intent(in) :: rhoq(qa,ka,ia,ja)
3457  real(RP), intent(in) :: xq(5,ka,ia,ja)
3458  !
3459  real(RP), parameter :: temc_min = -65.0_rp
3460  real(RP), parameter :: a_het = 0.2_rp ! SB06 (44)
3461  real(RP), parameter :: b_het = 0.65_rp ! SB06 (44)
3462  !
3463  real(RP) :: coef_m2_c
3464  real(RP) :: coef_m2_r
3465  ! temperature [celsius]
3466  real(RP) :: temc, temc2, temc3, temc4
3467  ! temperature function of homegenous/heterogenous freezing
3468  real(RP) :: jhom, jhet
3469  real(RP) :: rdt
3470  !
3471  integer :: i,j,k
3472  !
3473  rdt = 1.0_rp/dt
3474  !
3475  coef_m2_c = coef_m2(i_qc)
3476  coef_m2_r = coef_m2(i_qr)
3477  !
3478  profile_start("sn14_freezing")
3479  do j = js, je
3480  do i = is, ie
3481  do k = ks, ke
3482  temc = max( tem(k,i,j) - t00, temc_min )
3483  ! These cause from aerosol-droplet interaction.
3484  ! Bigg(1953) formula, Khain etal.(2000) eq.(4.5), Pruppacher and Klett(1997) eq.(9-48)
3485  jhet = a_het*exp( -b_het*temc - 1.0_rp )
3486  ! These cause in nature.
3487  ! Cotton and Field 2002, QJRMS. (12)
3488  if( temc < -65.0_rp )then
3489  jhom = 10.0_rp**(24.37236_rp)*1.e+3_rp
3490  jhet = a_het*exp( 65.0_rp*b_het - 1.0_rp ) ! 09/04/14 [Add], fixer T.Mitsui
3491  else if( temc < -30.0_rp ) then
3492  temc2 = temc*temc
3493  temc3 = temc*temc2
3494  temc4 = temc2*temc2
3495  jhom = 10.0_rp**(&
3496  - 243.40_rp - 14.75_rp*temc - 0.307_rp*temc2 &
3497  - 0.00287_rp*temc3 - 0.0000102_rp*temc4 ) *1.e+3_rp
3498  else if( temc < 0.0_rp) then
3499  jhom = 10._rp**(-7.63_rp-2.996_rp*(temc+30.0_rp))*1.e+3_rp
3500  else
3501  jhom = 0.0_rp
3502  jhet = 0.0_rp
3503  end if
3504  ! Note, xc should be limited in range[xc_min:xc_max].
3505  ! and PNChom need to be calculated by NC
3506  ! because reduction rate of Nc need to be bound by NC.
3507  ! For the same reason PLChom also bound by LC and xc.
3508  ! Basically L and N should be independent
3509  ! but average particle mass x should be in suitable range.
3510  ! Homogenous Freezing
3511  pq(i_lchom,k,i,j) = 0.0_rp
3512  pq(i_nchom,k,i,j) = 0.0_rp
3513  ! Heterogenous Freezing
3514  pq(i_lchet,k,i,j) = -rdt*rhoq(i_qc,k,i,j)*( 1.0_rp - exp( -coef_m2_c*xq(i_c,k,i,j)*(jhet+jhom)*dt ) )
3515  pq(i_nchet,k,i,j) = -rdt*rhoq(i_nc,k,i,j)*( 1.0_rp - exp( - xq(i_c,k,i,j)*(jhet+jhom)*dt ) )
3516  pq(i_lrhet,k,i,j) = -rdt*rhoq(i_qr,k,i,j)*( 1.0_rp - exp( -coef_m2_r*xq(i_r,k,i,j)*(jhet+jhom)*dt ) )
3517  pq(i_nrhet,k,i,j) = -rdt*rhoq(i_nr,k,i,j)*( 1.0_rp - exp( - xq(i_r,k,i,j)*(jhet+jhom)*dt ) )
3518  end do
3519  end do
3520  end do
3521  profile_stop("sn14_freezing")
3522  !
3523  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the call graph for this function:
Here is the caller graph for this function:

◆ update_by_phase_change_kij()

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

Definition at line 3720 of file scale_atmos_phy_mp_sn14.F90.

References scale_atmos_thermodyn::aq_cp, scale_atmos_thermodyn::aq_cv, scale_atmos_saturation::atmos_saturation_dqsi_dtem_dpre(), scale_atmos_saturation::atmos_saturation_dqsi_dtem_rho(), scale_atmos_saturation::atmos_saturation_dqsw_dtem_dpre(), scale_atmos_saturation::atmos_saturation_dqsw_dtem_rho(), scale_grid_index::ie, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_l, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, and scale_grid_index::ks.

Referenced by atmos_phy_mp_sn14().

3720  use scale_atmos_thermodyn, only: &
3721  aq_cv, &
3722  aq_cp
3723  use scale_atmos_saturation, only: &
3724  moist_pres2qsat_liq => atmos_saturation_pres2qsat_liq, &
3725  moist_pres2qsat_ice => atmos_saturation_pres2qsat_ice, &
3726  moist_dqsw_dtem_rho => atmos_saturation_dqsw_dtem_rho, &
3727  moist_dqsi_dtem_rho => atmos_saturation_dqsi_dtem_rho, &
3728  moist_dqsw_dtem_dpre => atmos_saturation_dqsw_dtem_dpre, &
3729  moist_dqsi_dtem_dpre => atmos_saturation_dqsi_dtem_dpre
3730  implicit none
3731 
3732  integer, intent(in) :: ntdiv ! [Add] 10/08/03
3733  integer, intent(in) :: ntmax ! [Add] 10/08/03
3734  !
3735  real(DP), intent(in) :: dt ! time step[s]
3736  real(RP), intent(in) :: gsgam2(ka,ia,ja) ! metric
3737  real(RP), intent(in) :: z(ka) ! altitude [m]
3738  real(RP), intent(in) :: dz(ka) ! altitude [m]
3739  real(RP), intent(in) :: velz(ka,ia,ja) ! vertical velocity @ half point[m/s]
3740  real(RP), intent(in) :: dtdt_rad(ka,ia,ja) ! temperture tendency by radiation[K/s]
3741  real(RP), intent(in) :: rho(ka,ia,ja) ! density[kg/m3]
3742  real(RP), intent(inout) :: rhoe(ka,ia,ja) ! internal energy[J/m3]
3743  real(RP), intent(inout) :: rhoq(qa,ka,ia,ja) ! tracers[kg/m3]
3744  real(RP), intent(inout) :: q(ka,ia,ja,qa) ! tracers mixing ratio[kg/kg]
3745  real(RP), intent(inout) :: tem(ka,ia,ja) ! temperature[K]
3746  real(RP), intent(inout) :: pre(ka,ia,ja) ! pressure[Pa]
3747  real(RP), intent(out) :: cva(ka,ia,ja) ! specific heat at constant volume
3748  real(RP), intent(in) :: esw(ka,ia,ja) ! saturated vapor pressure for liquid
3749  real(RP), intent(in) :: esi(ka,ia,ja) ! for ice
3750  real(RP), intent(in) :: rhoq2(qa,ka,ia,ja)
3751  !+++ tendency[kg/m3/s]
3752  real(RP), intent(inout) :: pq(pq_max,ka,ia,ja)
3753  real(RP), intent(out) :: qc_evaporate(ka,ia,ja)
3754  !+++ Column integrated tendency[kg/m2/s]
3755  real(RP), intent(inout) :: sl_plcdep(ia,ja)
3756  real(RP), intent(inout) :: sl_plrdep(ia,ja), sl_pnrdep(ia,ja)
3757  !
3758  real(RP) :: rmoist
3759  !
3760  real(RP) :: xi ! mean mass of ice particles
3761  real(RP) :: rrho ! 1/rho
3762  real(RP) :: wtem(ka,ia,ja) ! temperature[K]
3763  real(RP) :: qdry ! mixing ratio of dry air
3764  !
3765  real(RP) :: r_cva ! specific heat at constant volume
3766  real(RP) :: r_cpa ! specific heat at constant pressure
3767  real(RP) :: qsw(ka,ia,ja) ! saturated mixing ratio for liquid
3768  real(RP) :: qsi(ka,ia,ja) ! saturated mixing ratio for solid
3769  real(RP) :: dqswdtem_rho(ka,ia,ja) ! (dqsw/dtem)_rho
3770  real(RP) :: dqsidtem_rho(ka,ia,ja) ! (dqsi/dtem)_rho
3771  real(RP) :: dqswdtem_pre(ka,ia,ja) ! (dqsw/dtem)_pre
3772  real(RP) :: dqsidtem_pre(ka,ia,ja) ! (dqsi/dtem)_pre
3773  real(RP) :: dqswdpre_tem(ka,ia,ja) ! (dqsw/dpre)_tem
3774  real(RP) :: dqsidpre_tem(ka,ia,ja) ! (dqsi/dpre)_tem
3775  !
3776  real(RP) :: w ! vetical velocity[m/s]
3777  real(RP) :: acnd ! Pdynliq + Bergeron-Findeisen
3778  real(RP) :: adep ! Pdyndep + Bergeron-Findeisen
3779  real(RP) :: aliqliq, asolliq
3780  real(RP) :: aliqsol, asolsol
3781  real(RP) :: pdynliq ! production term of ssw by vertical motion
3782  real(RP) :: pdynsol ! production term of ssi by vertical motion
3783  real(RP) :: pradliq ! production term of ssw by radiation
3784  real(RP) :: pradsol ! production term of ssi by radiation
3785  real(RP) :: taucnd, r_taucnd ! time scale of ssw change by MP
3786  real(RP) :: taudep, r_taudep ! time scale of ssi change by MP
3787  real(RP) :: taucnd_c(ka,ia,ja), r_taucnd_c ! by cloud
3788  real(RP) :: taucnd_r(ka,ia,ja), r_taucnd_r ! by rain
3789  real(RP) :: taudep_i(ka,ia,ja), r_taudep_i ! by ice
3790  real(RP) :: taudep_s(ka,ia,ja), r_taudep_s ! by snow
3791  real(RP) :: taudep_g(ka,ia,ja), r_taudep_g ! by graupel
3792  ! alternative tendency through changing ssw and ssi
3793  real(RP) :: pncdep ! [Add] 11/08/30 T.Mitsui
3794  real(RP) :: plr2nr, pli2ni, pls2ns, plg2ng
3795  real(RP) :: coef_a_cnd, coef_b_cnd
3796  real(RP) :: coef_a_dep, coef_b_dep
3797  !
3798  real(RP) :: frz_dqc
3799  real(RP) :: frz_dnc
3800  real(RP) :: frz_dqr
3801  real(RP) :: frz_dnr
3802  real(RP) :: mlt_dqi
3803  real(RP) :: mlt_dni
3804  real(RP) :: mlt_dqs
3805  real(RP) :: mlt_dns
3806  real(RP) :: mlt_dqg
3807  real(RP) :: mlt_dng
3808  real(RP) :: dep_dqi
3809  real(RP) :: dep_dni
3810  real(RP) :: dep_dqs
3811  real(RP) :: dep_dns
3812  real(RP) :: dep_dqg
3813  real(RP) :: dep_dng
3814  real(RP) :: dep_dqr
3815  real(RP) :: dep_dnr
3816  real(RP) :: dep_dqc
3817  real(RP) :: dep_dnc ! 11/08/30 [Add] T.Mitsui, dep_dnc
3818  real(RP) :: r_xc_ccn, r_xi_ccn ! 11/08/30 [Add] T.Mitsui
3819  !
3820  real(RP) :: drhoqv
3821  real(RP) :: drhoqc, drhoqr, drhoqi, drhoqs, drhoqg
3822  real(RP) :: drhonc, drhonr, drhoni, drhons, drhong
3823  !
3824  real(RP) :: fac1, fac2, fac3, fac4, fac5, fac6
3825  real(RP) :: r_rvaptem ! 1/(Rvap*tem)
3826  real(RP) :: pv ! vapor pressure
3827  real(RP) :: lvsw, lvsi ! saturated vapor density
3828  real(RP) :: dlvsw, dlvsi !
3829  ! [Add] 11/08/30 T.Mitsui
3830  real(RP) :: dcnd, ddep ! total cndensation/deposition
3831  real(RP) :: uplim_cnd ! upper limit of condensation
3832  real(RP) :: lowlim_cnd ! lower limit of evaporation
3833  ! [Add] 11/08/30 T.Mitsui
3834  real(RP) :: uplim_dep ! upper limit of condensation
3835  real(RP) :: lowlim_dep ! lower limit of evaporation
3836  real(RP) :: ssw, ssi ! supersaturation ratio
3837  real(RP) :: r_esw, r_esi ! 1/esw, 1/esi
3838  real(RP) :: r_lvsw, r_lvsi ! 1/(lvsw*ssw), 1/(lvsi*ssi)
3839  real(RP) :: r_dt ! 1/dt
3840  real(RP) :: ssw_o, ssi_o
3841 ! real(RP) :: dt_dyn
3842 ! real(RP) :: dt_mp
3843  !
3844 ! real(RP) :: tem_lh(KA,IA,JA)
3845 ! real(RP) :: dtemdt_lh(KA,IA,JA)
3846  real(RP), save :: fac_cndc = 1.0_rp
3847  logical, save :: opt_fix_taucnd_c=.false.
3848  logical, save :: flag_first =.true.
3849  !
3850  namelist /nm_mp_sn14_condensation/ &
3851  opt_fix_taucnd_c, fac_cndc
3852 
3853  real(RP) :: fac_cndc_wrk
3854  !
3855  real(RP), parameter :: tau100day = 1.e+7_rp
3856  real(RP), parameter :: r_tau100day = 1.e-7_rp
3857  real(RP), parameter :: eps=1.e-30_rp
3858  !
3859  integer :: i,j,k,iqw
3860  real(RP) :: sw
3861  !
3862 
3863  ! [Add] 11/08/30 T.Mitsui
3864  if( flag_first )then
3865  flag_first = .false.
3866  rewind(io_fid_conf)
3867  read (io_fid_conf,nml=nm_mp_sn14_condensation, end=100)
3868 100 if( io_l ) write (io_fid_log,nml=nm_mp_sn14_condensation)
3869  end if
3870  !
3871 ! dt_dyn = dt*ntmax
3872 ! dt_mp = dt*(ntdiv-1)
3873  !
3874  r_dt = 1.0_rp/dt
3875  !
3876  r_xc_ccn=1.0_rp/xc_ccn
3877 ! r_xi_ccn=1.0_RP/xi_ccn
3878  !
3879  if( opt_fix_taucnd_c )then
3880  fac_cndc_wrk = fac_cndc**(1.0_rp-b_m(i_qc))
3881  do j = js, je
3882  do i = is, ie
3883  do k = ks, ke
3884  pq(i_lcdep,k,i,j) = pq(i_lcdep,k,i,j)*fac_cndc_wrk
3885  end do
3886  end do
3887  end do
3888  if( io_l ) write(io_fid_log,*) "taucnd:fac_cndc_wrk=",fac_cndc_wrk
3889  end if
3890 
3891 !OCL XFILL
3892  do j = js, je
3893  do i = is, ie
3894  do k = ks, ke
3895  ! Temperature lower limit is only used for saturation condition.
3896  ! On the other hand original "tem" is used for calculation of latent heat or energy equation.
3897  wtem(k,i,j) = max( tem(k,i,j), tem_min )
3898  end do
3899  end do
3900  end do
3901 
3902  call moist_pres2qsat_liq ( qsw, wtem, pre )
3903  call moist_pres2qsat_ice ( qsi, wtem, pre )
3904  call moist_dqsw_dtem_rho ( dqswdtem_rho, wtem, rho )
3905  call moist_dqsi_dtem_rho ( dqsidtem_rho, wtem, rho )
3906  call moist_dqsw_dtem_dpre( dqswdtem_pre, dqswdpre_tem, wtem, pre )
3907  call moist_dqsi_dtem_dpre( dqsidtem_pre, dqsidpre_tem, wtem, pre )
3908 
3909  profile_start("sn14_update")
3910  do j = js, je
3911  do i = is, ie
3912  do k = ks, ke
3913  if( z(k) <= 25000.0_rp )then
3914  w = 0.5_rp*(velz(k,i,j) + velz(k+1,i,j))
3915  else
3916  w = 0.0_rp
3917  end if
3918  if( pre(k,i,j) < esw(k,i,j)+1.e-10_rp )then
3919  qsw(k,i,j) = 1.0_rp
3920  dqswdtem_rho(k,i,j) = 0.0_rp
3921  dqswdtem_pre(k,i,j) = 0.0_rp
3922  dqswdpre_tem(k,i,j) = 0.0_rp
3923  end if
3924  if( pre(k,i,j) < esi(k,i,j)+1.e-10_rp )then
3925  qsi(k,i,j) = 1.0_rp
3926  dqsidtem_rho(k,i,j) = 0.0_rp
3927  dqsidtem_pre(k,i,j) = 0.0_rp
3928  dqsidpre_tem(k,i,j) = 0.0_rp
3929  end if
3930 
3931  r_rvaptem = 1.0_rp/(rvap*wtem(k,i,j))
3932  lvsw = esw(k,i,j)*r_rvaptem ! rho=p/(Rv*T)
3933  lvsi = esi(k,i,j)*r_rvaptem !
3934  pv = rhoq2(i_qv,k,i,j)*rvap*tem(k,i,j)
3935  r_esw = 1.0_rp/esw(k,i,j)
3936  r_esi = 1.0_rp/esi(k,i,j)
3937  ssw = min( mp_ssw_lim, ( pv*r_esw-1.0_rp ) )
3938  ssi = pv*r_esi - 1.0_rp
3939  r_lvsw = 1.0_rp/lvsw
3940  r_lvsi = 1.0_rp/lvsi
3941  r_taucnd_c = pq(i_lcdep,k,i,j)*r_lvsw
3942  r_taucnd_r = pq(i_lrdep,k,i,j)*r_lvsw
3943  r_taudep_i = pq(i_lidep,k,i,j)*r_lvsi
3944  r_taudep_s = pq(i_lsdep,k,i,j)*r_lvsi
3945  r_taudep_g = pq(i_lgdep,k,i,j)*r_lvsi
3946 ! taucnd_c(k,i,j) = 1.0_RP/(r_taucnd_c+r_tau100day)
3947 ! taucnd_r(k,i,j) = 1.0_RP/(r_taucnd_r+r_tau100day)
3948 ! taudep_i(k,i,j) = 1.0_RP/(r_taudep_i+r_tau100day)
3949 ! taudep_s(k,i,j) = 1.0_RP/(r_taudep_s+r_tau100day)
3950 ! taudep_g(k,i,j) = 1.0_RP/(r_taudep_g+r_tau100day)
3951 
3952  calc_qdry( qdry, q, k, i, j, iqw )
3953  calc_cv( cva(k,i,j), qdry, q, k, i, j, iqw, cvdry, aq_cv )
3954  calc_r( rmoist, q(k,i,j,i_qv), qdry, rdry, rvap )
3955  r_cva = 1.0_rp / cva(k,i,j)
3956  r_cpa = 1.0_rp / (cva(k,i,j) + rmoist)
3957 
3958  ! Coefficient of latent heat release for ssw change by PLCdep and PLRdep
3959  aliqliq = 1.0_rp &
3960  + r_cva*( lhv00 + (cvvap-cl)*tem(k,i,j) )*dqswdtem_rho(k,i,j)
3961  ! Coefficient of latent heat release for ssw change by PLIdep, PLSdep and PLGdep
3962  asolliq = 1.0_rp &
3963  + r_cva*( lhv00 + lhf00 + (cvvap-ci)*tem(k,i,j) )*dqswdtem_rho(k,i,j)
3964  ! Coefficient of latent heat release for ssi change by PLCdep and PLRdep
3965  aliqsol = 1.0_rp &
3966  + r_cva*( lhv00 + (cvvap-cl)*tem(k,i,j) )*dqsidtem_rho(k,i,j)
3967  ! Coefficient of latent heat release for ssi change by PLIdep, PLSdep and PLGdep
3968  asolsol = 1.0_rp &
3969  + r_cva*( lhv00 + lhf00 + (cvvap-ci)*tem(k,i,j) )*dqsidtem_rho(k,i,j)
3970  pdynliq = w * grav * ( r_cpa*dqswdtem_pre(k,i,j) + rho(k,i,j)*dqswdpre_tem(k,i,j) )
3971  pdynsol = w * grav * ( r_cpa*dqsidtem_pre(k,i,j) + rho(k,i,j)*dqsidpre_tem(k,i,j) )
3972  pradliq = -dtdt_rad(k,i,j) * dqswdtem_rho(k,i,j)
3973  pradsol = -dtdt_rad(k,i,j) * dqsidtem_rho(k,i,j)
3974 
3975  ssw_o = ssw
3976  ssi_o = ssi
3977 ! ssw_o = ssw - Pdynliq*(dt_dyn-dt_mp)/qsw(k,i,j) + Pradliq*r_qsw*dt_mp
3978 ! ssi_o = ssi - Pdynsol*(dt_dyn-dt_mp)/qsi(k,i,j) + Pradsol*r_qsi*dt_mp
3979 
3980  acnd = pdynliq + pradliq &
3981  - ( r_taudep_i+r_taudep_s+r_taudep_g ) * ( qsw(k,i,j) - qsi(k,i,j) )
3982  adep = pdynsol + pradsol &
3983  + ( r_taucnd_c+r_taucnd_r ) * ( qsw(k,i,j) - qsi(k,i,j) )
3984  r_taucnd = &
3985  + aliqliq*( r_taucnd_c+r_taucnd_r ) &
3986  + asolliq*( r_taudep_i+r_taudep_s+r_taudep_g )
3987  r_taudep = &
3988  + aliqsol*( r_taucnd_c+r_taucnd_r )&
3989  + asolsol*( r_taudep_i+r_taudep_s+r_taudep_g )
3990 
3991  uplim_cnd = max( rho(k,i,j)*ssw_o*qsw(k,i,j)*r_dt, 0.0_rp )
3992  lowlim_cnd = min( rho(k,i,j)*ssw_o*qsw(k,i,j)*r_dt, 0.0_rp )
3993  if( r_taucnd < r_tau100day )then
3994 ! taucnd = tau100day
3995  pq(i_lcdep,k,i,j) = max(lowlim_cnd, min(uplim_cnd, pq(i_lcdep,k,i,j)*ssw_o ))
3996  pq(i_lrdep,k,i,j) = max(lowlim_cnd, min(uplim_cnd, pq(i_lrdep,k,i,j)*ssw_o ))
3997  pq(i_nrdep,k,i,j) = min(0.0_rp, pq(i_nrdep,k,i,j)*ssw_o )
3998 ! PLR2NR = 0.0_RP
3999  else
4000  taucnd = 1.0_rp/r_taucnd
4001  ! Production term for liquid water content
4002  coef_a_cnd = rho(k,i,j)*acnd*taucnd
4003  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 )
4004  pq(i_lcdep,k,i,j) = coef_a_cnd*r_taucnd_c - coef_b_cnd*r_taucnd_c
4005  plr2nr = pq(i_nrdep,k,i,j)/(pq(i_lrdep,k,i,j)+1.e-30_rp)
4006  pq(i_lrdep,k,i,j) = coef_a_cnd*r_taucnd_r - coef_b_cnd*r_taucnd_r
4007  pq(i_nrdep,k,i,j) = min(0.0_rp, pq(i_lrdep,k,i,j)*plr2nr )
4008  end if
4009 
4010  uplim_dep = max( rho(k,i,j)*ssi_o*qsi(k,i,j)*r_dt, 0.0_rp )
4011  lowlim_dep = min( rho(k,i,j)*ssi_o*qsi(k,i,j)*r_dt, 0.0_rp )
4012  if( r_taudep < r_tau100day )then
4013 ! taudep = tau100day
4014  pq(i_lidep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lidep,k,i,j)*ssi_o ))
4015  pq(i_lsdep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lsdep,k,i,j)*ssi_o ))
4016  pq(i_lgdep,k,i,j) = max(lowlim_dep, min(uplim_dep, pq(i_lgdep,k,i,j)*ssi_o ))
4017  pq(i_nidep,k,i,j) = min(0.0_rp, pq(i_nidep,k,i,j)*ssi_o )
4018  pq(i_nsdep,k,i,j) = min(0.0_rp, pq(i_nsdep,k,i,j)*ssi_o )
4019  pq(i_ngdep,k,i,j) = min(0.0_rp, pq(i_ngdep,k,i,j)*ssi_o )
4020  else
4021  taudep = 1.0_rp/r_taudep
4022  ! Production term for ice water content
4023  coef_a_dep = rho(k,i,j)*adep*taudep
4024  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 )
4025  pli2ni = pq(i_nidep,k,i,j)/max(pq(i_lidep,k,i,j),1.e-30_rp)
4026  pls2ns = pq(i_nsdep,k,i,j)/max(pq(i_lsdep,k,i,j),1.e-30_rp)
4027  plg2ng = pq(i_ngdep,k,i,j)/max(pq(i_lgdep,k,i,j),1.e-30_rp)
4028  pq(i_lidep,k,i,j) = coef_a_dep*r_taudep_i - coef_b_dep*r_taudep_i
4029  pq(i_lsdep,k,i,j) = coef_a_dep*r_taudep_s - coef_b_dep*r_taudep_s
4030  pq(i_lgdep,k,i,j) = coef_a_dep*r_taudep_g - coef_b_dep*r_taudep_g
4031  pq(i_nidep,k,i,j) = min(0.0_rp, pq(i_lidep,k,i,j)*pli2ni )
4032  pq(i_nsdep,k,i,j) = min(0.0_rp, pq(i_lsdep,k,i,j)*pls2ns )
4033  pq(i_ngdep,k,i,j) = min(0.0_rp, pq(i_lgdep,k,i,j)*plg2ng )
4034  end if
4035 
4036  sw = 0.5_rp - sign(0.5_rp, pq(i_lcdep,k,i,j)+eps) != 1 for PLCdep<=-eps
4037  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
4038 ! if( PQ(I_LCdep,k,i,j) < -eps )then
4039 ! 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 )
4040 ! else
4041 ! PNCdep = 0.0_RP
4042 ! end if
4043 ! if( PQ(I_LIdep,k,i,j) < -eps )then
4044 ! 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 )
4045 ! else
4046 ! PQ(I_NIdep,k,i,j) = 0.0_RP
4047 ! end if
4048 
4049  !--- evaporation/condensation
4050  r_rvaptem = 1.0_rp/(rvap*wtem(k,i,j))
4051  lvsw = esw(k,i,j)*r_rvaptem
4052  dlvsw = rhoq2(i_qv,k,i,j)-lvsw
4053  dcnd = dt*(pq(i_lcdep,k,i,j)+pq(i_lrdep,k,i,j))
4054 
4055  sw = ( sign(0.5_rp,dcnd) + sign(0.5_rp,dlvsw) ) &
4056  * ( 0.5_rp + sign(0.5_rp,abs(dcnd)-eps) ) ! to avoid zero division
4057  ! sw= 1: always supersaturated
4058  ! sw=-1: always unsaturated
4059  ! sw= 0: partially unsaturated during timestep
4060  fac1 = min(dlvsw*sw,dcnd*sw)*sw / (abs(sw)-1.0_rp+dcnd) & ! sw=1,-1
4061  + 1.0_rp - abs(sw) ! sw=0
4062  dep_dqc = max( dt*pq(i_lcdep,k,i,j)*fac1, &
4063  -rhoq2(i_qc,k,i,j) - 1e30_rp*(sw+1.0_rp) )*abs(sw) != -lc for sw=-1, -inf for sw=1
4064  dep_dqr = max( dt*pq(i_lrdep,k,i,j)*fac1, &
4065  -rhoq2(i_qr,k,i,j) - 1e30_rp*(sw+1.0_rp) )*abs(sw) != -lr for sw=-1, -inf for sw=1
4066 ! if ( (dcnd > eps) .AND. (dlvsw > eps) )then
4067 ! ! always supersaturated
4068 ! fac1 = min(dlvsw,dcnd)/dcnd
4069 ! dep_dqc = dt*PQ(I_LCdep,k,i,j)*fac1
4070 ! dep_dqr = dt*PQ(I_LRdep,k,i,j)*fac1
4071 ! else if( (dcnd < -eps) .AND. (dlvsw < -eps) )then
4072 ! ! always unsaturated
4073 ! fac1 = max( dlvsw,dcnd )/dcnd
4074 ! dep_dqc = max( dt*PQ(I_LCdep,k,i,j)*fac1, -rhoq2(I_QC,k,i,j) )
4075 ! dep_dqr = max( dt*PQ(I_LRdep,k,i,j)*fac1, -rhoq2(I_QR,k,i,j) )
4076 ! else
4077 ! ! partially unsaturated during timestep
4078 ! fac1 = 1.0_RP
4079 ! dep_dqc = 0.0_RP
4080 ! dep_dqr = 0.0_RP
4081 ! end if
4082 
4083  ! evaporation always lose number(always negative).
4084  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
4085  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
4086 
4087  qc_evaporate(k,i,j) = - dep_dnc ! [Add] Y.Sato 15/09/08
4088 
4089  !--- deposition/sublimation
4090  lvsi = esi(k,i,j)*r_rvaptem
4091  ddep = dt*(pq(i_lidep,k,i,j)+pq(i_lsdep,k,i,j)+pq(i_lgdep,k,i,j))
4092  dlvsi = rhoq2(i_qv,k,i,j)-lvsi ! limiter for esi>1.d0
4093 
4094  sw = ( sign(0.5_rp,ddep) + sign(0.5_rp,dlvsi) ) &
4095  * ( 0.5_rp + sign(0.5_rp,abs(ddep)-eps) ) ! to avoid zero division
4096  ! sw= 1: always supersaturated
4097  ! sw=-1: always unsaturated
4098  ! sw= 0: partially unsaturated during timestep
4099  fac2 = min(dlvsi*sw,ddep*sw)*sw / (abs(sw)-1.0_rp+ddep) & ! sw=1,-1
4100  + 1.0_rp - abs(sw) ! sw=0
4101  dep_dqi = max( dt*pq(i_lidep,k,i,j) &
4102  * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), & != fac2 for sw=-1,1, 1 for sw=0
4103  -rhoq2(i_qi,k,i,j) - 1e30_rp*(sw+1.0_rp) ) != -li for sw=-1, -inf for sw=0,1
4104  dep_dqs = max( dt*pq(i_lsdep,k,i,j) &
4105  * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), & != fac2 for sw=-1,1, 1 for sw=0
4106  -rhoq2(i_qs,k,i,j) - 1e30_rp*(sw+1.0_rp) ) != -ls for sw=-1, -inf for sw=0,1
4107  dep_dqg = max( dt*pq(i_lgdep,k,i,j) &
4108  * ( 1.0_rp-abs(sw) + fac2*abs(sw) ), & != fac2 for sw=-1,1, 1 for sw=0
4109  -rhoq2(i_qg,k,i,j) - 1e30_rp*(sw+1.0_rp) ) != -lg for sw=-1, -inf for sw=0,1
4110 ! if ( (ddep > eps) .AND. (dlvsi > eps) )then
4111 ! ! always supersaturated
4112 ! fac2 = min(dlvsi,ddep)/ddep
4113 ! dep_dqi = dt*PQ(I_LIdep,k,i,j)*fac2
4114 ! dep_dqs = dt*PQ(I_LSdep,k,i,j)*fac2
4115 ! dep_dqg = dt*PQ(I_LGdep,k,i,j)*fac2
4116 ! else if ( (ddep < -eps) .AND. (dlvsi < -eps) )then
4117 ! ! always unsaturated
4118 ! fac2 = max(dlvsi,ddep)/ddep
4119 ! dep_dqi = max(dt*PQ(I_LIdep,k,i,j)*fac2, -rhoq2(I_QI,k,i,j) )
4120 ! dep_dqs = max(dt*PQ(I_LSdep,k,i,j)*fac2, -rhoq2(I_QS,k,i,j) )
4121 ! dep_dqg = max(dt*PQ(I_LGdep,k,i,j)*fac2, -rhoq2(I_QG,k,i,j) )
4122 ! else
4123 ! ! partially unsaturated during timestep
4124 ! fac2 = 1.0_RP
4125 ! dep_dqi = dt*PQ(I_LIdep,k,i,j)
4126 ! dep_dqs = dt*PQ(I_LSdep,k,i,j)
4127 ! dep_dqg = dt*PQ(I_LGdep,k,i,j)
4128 ! end if
4129 
4130  ! evaporation always lose number(always negative).
4131  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
4132  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
4133  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
4134 
4135  !--- freezing of cloud drop
4136  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
4137  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
4138  fac3 = ( frz_dqc-eps )/( dt*(pq(i_lchom,k,i,j)+pq(i_lchet,k,i,j))-eps )
4139  fac4 = ( frz_dnc-eps )/( dt*(pq(i_nchom,k,i,j)+pq(i_nchet,k,i,j))-eps )
4140  pq(i_lchom,k,i,j) = fac3*pq(i_lchom,k,i,j)
4141  pq(i_lchet,k,i,j) = fac3*pq(i_lchet,k,i,j)
4142  pq(i_nchom,k,i,j) = fac4*pq(i_nchom,k,i,j)
4143  pq(i_nchet,k,i,j) = fac4*pq(i_nchet,k,i,j)
4144 
4145  !--- melting
4146  ! ice change
4147  mlt_dqi = max( dt*pq(i_limlt,k,i,j), -rhoq2(i_qi,k,i,j)-dep_dqi ) ! negative value
4148  mlt_dni = max( dt*pq(i_nimlt,k,i,j), -rhoq2(i_ni,k,i,j)-dep_dni ) ! negative value
4149  ! snow change
4150  mlt_dqs = max( dt*pq(i_lsmlt,k,i,j), -rhoq2(i_qs,k,i,j)-dep_dqs ) ! negative value
4151  mlt_dns = max( dt*pq(i_nsmlt,k,i,j), -rhoq2(i_ns,k,i,j)-dep_dns ) ! negative value
4152  ! graupel change
4153  mlt_dqg = max( dt*pq(i_lgmlt,k,i,j), -rhoq2(i_qg,k,i,j)-dep_dqg ) ! negative value
4154  mlt_dng = max( dt*pq(i_ngmlt,k,i,j), -rhoq2(i_ng,k,i,j)-dep_dng ) ! negative value
4155 
4156  !--- freezing of larger droplets
4157  frz_dqr = max( dt*(pq(i_lrhet,k,i,j)), min(0.0_rp, -rhoq2(i_qr,k,i,j)-dep_dqr) ) ! negative value
4158  frz_dnr = max( dt*(pq(i_nrhet,k,i,j)), min(0.0_rp, -rhoq2(i_nr,k,i,j)-dep_dnr) ) ! negative value
4159 
4160  fac5 = ( frz_dqr-eps )/( dt*pq(i_lrhet,k,i,j)-eps )
4161  pq(i_lrhet,k,i,j) = fac5*pq(i_lrhet,k,i,j)
4162  fac6 = ( frz_dnr-eps )/( dt*pq(i_nrhet,k,i,j)-eps )
4163  pq(i_nrhet,k,i,j) = fac6*pq(i_nrhet,k,i,j)
4164 
4165  ! water vapor change
4166  drhoqv = -(dep_dqc+dep_dqi+dep_dqs+dep_dqg+dep_dqr)
4167 
4168  rhoq(i_qv,k,i,j) = max(0.0_rp, rhoq(i_qv,k,i,j) + drhoqv )
4169 
4170  rhoe(k,i,j) = rhoe(k,i,j) - lhv * drhoqv
4171 
4172  xi = min(xi_max, max(xi_min, rhoq2(i_qi,k,i,j)/(rhoq2(i_ni,k,i,j)+ni_min) ))
4173  sw = 0.5_rp + sign(0.5_rp,xi-x_sep) ! if (xi>=x_sep) then sw=1 else sw=0
4174  ! sw=1: large ice crystals turn into rain by melting
4175 
4176  ! total cloud change
4177  drhoqc = ( frz_dqc - mlt_dqi*(1.0_rp-sw) + dep_dqc )
4178  drhonc = ( frz_dnc - mlt_dni*(1.0_rp-sw) + dep_dnc )
4179  ! total rain change
4180  drhoqr = ( frz_dqr - mlt_dqg - mlt_dqs - mlt_dqi*sw + dep_dqr )
4181  drhonr = ( frz_dnr - mlt_dng - mlt_dns - mlt_dni*sw + dep_dnr )
4182 
4183  rhoq(i_qc,k,i,j) = max(0.0_rp, rhoq(i_qc,k,i,j) + drhoqc )
4184  rhoq(i_nc,k,i,j) = max(0.0_rp, rhoq(i_nc,k,i,j) + drhonc )
4185  rhoq(i_qr,k,i,j) = max(0.0_rp, rhoq(i_qr,k,i,j) + drhoqr )
4186  rhoq(i_nr,k,i,j) = max(0.0_rp, rhoq(i_nr,k,i,j) + drhonr )
4187 
4188  ! total ice change
4189  drhoqi = (-frz_dqc + mlt_dqi + dep_dqi )
4190  drhoni = (-frz_dnc + mlt_dni + dep_dni )
4191 
4192  rhoq(i_qi,k,i,j) = max(0.0_rp, rhoq(i_qi,k,i,j) + drhoqi )
4193  rhoq(i_ni,k,i,j) = max(0.0_rp, rhoq(i_ni,k,i,j) + drhoni )
4194 
4195  rhoe(k,i,j) = rhoe(k,i,j) + lhf * drhoqi
4196 
4197  ! total snow change
4198  drhoqs = ( mlt_dqs + dep_dqs )
4199  drhons = ( mlt_dns + dep_dns )
4200 
4201  rhoq(i_qs,k,i,j) = max(0.0_rp, rhoq(i_qs,k,i,j) + drhoqs )
4202  rhoq(i_ns,k,i,j) = max(0.0_rp, rhoq(i_ns,k,i,j) + drhons )
4203 
4204  rhoe(k,i,j) = rhoe(k,i,j) + lhf * drhoqs
4205 
4206  ! total graupel change
4207  drhoqg = (-frz_dqr + mlt_dqg + dep_dqg )
4208  drhong = (-frz_dnr + mlt_dng + dep_dng )
4209 
4210  rhoq(i_qg,k,i,j) = max(0.0_rp, rhoq(i_qg,k,i,j) + drhoqg )
4211  rhoq(i_ng,k,i,j) = max(0.0_rp, rhoq(i_ng,k,i,j) + drhong )
4212 
4213  rhoe(k,i,j) = rhoe(k,i,j) + lhf * drhoqg
4214 
4215  !--- update mixing ratio
4216  rrho = 1.0_rp/rho(k,i,j)
4217 
4218  q(k,i,j,i_qv) = rhoq(i_qv,k,i,j) * rrho
4219  q(k,i,j,i_qc) = rhoq(i_qc,k,i,j) * rrho
4220  q(k,i,j,i_qr) = rhoq(i_qr,k,i,j) * rrho
4221  q(k,i,j,i_qi) = rhoq(i_qi,k,i,j) * rrho
4222  q(k,i,j,i_qs) = rhoq(i_qs,k,i,j) * rrho
4223  q(k,i,j,i_qg) = rhoq(i_qg,k,i,j) * rrho
4224  q(k,i,j,i_nc) = rhoq(i_nc,k,i,j) * rrho
4225  q(k,i,j,i_nr) = rhoq(i_nr,k,i,j) * rrho
4226  q(k,i,j,i_ni) = rhoq(i_ni,k,i,j) * rrho
4227  q(k,i,j,i_ns) = rhoq(i_ns,k,i,j) * rrho
4228  q(k,i,j,i_ng) = rhoq(i_ng,k,i,j) * rrho
4229 
4230  calc_qdry( qdry, q, k, i, j, iqw )
4231  calc_cv( cva(k,i,j), qdry, q, k, i, j, iqw, cvdry, aq_cv )
4232  calc_r( rmoist, q(k,i,j,i_qv), qdry, rdry, rvap )
4233  tem(k,i,j) = rhoe(k,i,j) / ( rho(k,i,j) * cva(k,i,j) )
4234  pre(k,i,j) = rho(k,i,j) * rmoist * tem(k,i,j)
4235 
4236  sl_plcdep(i,j) = sl_plcdep(i,j) + dep_dqc*dz(k)*gsgam2(k,i,j)
4237  sl_plrdep(i,j) = sl_plrdep(i,j) + dep_dqr*dz(k)*gsgam2(k,i,j)
4238  sl_pnrdep(i,j) = sl_pnrdep(i,j) + dep_dnr*dz(k)*gsgam2(k,i,j)
4239  end do
4240  end do
4241  end do
4242  profile_stop("sn14_update")
4243 
4244  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_atmos_saturation
module ATMOSPHERE / Saturation adjustment
Definition: scale_atmos_sub_saturation.F90:16
scale_atmos_thermodyn::aq_cp
real(rp), dimension(:), allocatable, public aq_cp
CP for each hydrometeors [J/kg/K].
Definition: scale_atmos_sub_thermodyn.F90:142
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_atmos_saturation::atmos_saturation_dqsw_dtem_rho
subroutine, public atmos_saturation_dqsw_dtem_rho(dqsdtem, temp, dens)
Definition: scale_atmos_sub_saturation.F90:1242
scale_atmos_saturation::atmos_saturation_dqsw_dtem_dpre
subroutine, public atmos_saturation_dqsw_dtem_dpre(dqsdtem, dqsdpre, temp, pres)
Definition: scale_atmos_sub_saturation.F90:1327
scale_atmos_saturation::atmos_saturation_dqsi_dtem_rho
subroutine, public atmos_saturation_dqsi_dtem_rho(dqsdtem, temp, dens)
Definition: scale_atmos_sub_saturation.F90:1285
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_atmos_saturation::atmos_saturation_dqsi_dtem_dpre
subroutine, public atmos_saturation_dqsi_dtem_dpre(dqsdtem, dqsdpre, temp, pres)
Definition: scale_atmos_sub_saturation.F90:1375
scale_atmos_thermodyn::aq_cv
real(rp), dimension(:), allocatable, public aq_cv
CV for each hydrometeors [J/kg/K].
Definition: scale_atmos_sub_thermodyn.F90:143
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_atmos_thermodyn
module ATMOSPHERE / Thermodynamics
Definition: scale_atmos_sub_thermodyn.F90:19
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the call graph for this function:
Here is the caller graph for this function:

◆ mp_negativefilter()

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

Definition at line 4250 of file scale_atmos_phy_mp_sn14.F90.

References scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, scale_prof::prof_rapend(), scale_prof::prof_rapstart(), and scale_tracer::qa.

Referenced by atmos_phy_mp_sn14().

4250  implicit none
4251  real(RP), intent(inout) :: dens(ka,ia,ja)
4252  real(RP), intent(inout) :: qtrc(ka,ia,ja,qa)
4253 
4254  real(RP) :: diffq(ka,ia,ja)
4255  real(RP) :: r_xmin
4256 
4257  integer :: k, i, j, iq
4258  !---------------------------------------------------------------------------
4259 
4260  call prof_rapstart('MP_filter', 3)
4261 
4262  r_xmin = 1.0_rp / xmin_filter
4263 
4264  ! total hydrometeor (before correction)
4265  do j = js, je
4266  do i = is, ie
4267 
4268  do k = ks, ke
4269  diffq(k,i,j) = qtrc(k,i,j,i_qv) &
4270  + qtrc(k,i,j,i_qc) &
4271  + qtrc(k,i,j,i_qr) &
4272  + qtrc(k,i,j,i_qi) &
4273  + qtrc(k,i,j,i_qs) &
4274  + qtrc(k,i,j,i_qg)
4275  enddo
4276 
4277  ! remove negative value of hydrometeor (mass, number)
4278  do iq = 1, qa
4279  do k = ks, ke
4280  qtrc(k,i,j,iq) = max(0.0_rp, qtrc(k,i,j,iq))
4281  enddo
4282  enddo
4283 
4284  ! apply correction of hydrometeor to total density
4285  ! [note] mass conservation is broken here to fill rounding error.
4286  do k = ks, ke
4287  dens(k,i,j) = dens(k,i,j) &
4288  * ( 1.0_rp &
4289  + qtrc(k,i,j,i_qv) &
4290  + qtrc(k,i,j,i_qc) &
4291  + qtrc(k,i,j,i_qr) &
4292  + qtrc(k,i,j,i_qi) &
4293  + qtrc(k,i,j,i_qs) &
4294  + qtrc(k,i,j,i_qg) &
4295  - diffq(k,i,j) ) ! after-before
4296  enddo
4297 
4298  ! avoid unrealistical value of number concentration
4299  ! due to numerical diffusion in advection
4300 
4301  do k = ks, ke
4302  if ( qtrc(k,i,j,i_nc) > qtrc(k,i,j,i_qc)*r_xmin ) then
4303  qtrc(k,i,j,i_nc) = qtrc(k,i,j,i_qc)*r_xmin
4304  endif
4305  enddo
4306  do k = ks, ke
4307  if ( qtrc(k,i,j,i_nr) > qtrc(k,i,j,i_qr)*r_xmin ) then
4308  qtrc(k,i,j,i_nr) = qtrc(k,i,j,i_qr)*r_xmin
4309  endif
4310  enddo
4311  do k = ks, ke
4312  if ( qtrc(k,i,j,i_ni) > qtrc(k,i,j,i_qi)*r_xmin ) then
4313  qtrc(k,i,j,i_ni) = qtrc(k,i,j,i_qi)*r_xmin
4314  endif
4315  enddo
4316  do k = ks, ke
4317  if ( qtrc(k,i,j,i_ns) > qtrc(k,i,j,i_qs)*r_xmin ) then
4318  qtrc(k,i,j,i_ns) = qtrc(k,i,j,i_qs)*r_xmin
4319  endif
4320  enddo
4321  do k = ks, ke
4322  if ( qtrc(k,i,j,i_ng) > qtrc(k,i,j,i_qg)*r_xmin ) then
4323  qtrc(k,i,j,i_ng) = qtrc(k,i,j,i_qg)*r_xmin
4324  endif
4325  enddo
4326 
4327  enddo
4328  enddo
4329 
4330  call prof_rapend('MP_filter', 3)
4331 
4332  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_cloudfraction()

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

Calculate Cloud Fraction.

Definition at line 4342 of file scale_atmos_phy_mp_sn14.F90.

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

Referenced by scale_atmos_phy_mp::atmos_phy_mp_setup().

4342  use scale_grid_index
4343  use scale_tracer, only: &
4344  qad => qa
4345  implicit none
4346 
4347  real(RP), intent(out) :: cldfrac(ka,ia,ja)
4348  real(RP), intent(in) :: qtrc (ka,ia,ja,qad)
4349  real(RP), intent(in) :: mask_criterion
4350 
4351  real(RP) :: qhydro
4352  integer :: k, i, j, iq
4353  !---------------------------------------------------------------------------
4354 
4355  do j = js, je
4356  do i = is, ie
4357  do k = ks, ke
4358  qhydro = 0.0_rp
4359  do iq = 1, mp_qa
4360  qhydro = qhydro + qtrc(k,i,j,i_mp2all(iq))
4361  enddo
4362  cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
4363  enddo
4364  enddo
4365  enddo
4366 
4367  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:18
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_effectiveradius()

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

Calculate Effective Radius.

Definition at line 4377 of file scale_atmos_phy_mp_sn14.F90.

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

Referenced by scale_atmos_phy_mp::atmos_phy_mp_setup().

4377  use scale_grid_index
4378  use scale_tracer, only: &
4379  qad => qa, &
4380  mp_qad => mp_qa
4381  implicit none
4382 
4383  real(RP), intent(out) :: re (ka,ia,ja,mp_qad) ! effective radius [cm]
4384  real(RP), intent(in) :: qtrc0(ka,ia,ja,qad) ! tracer mass concentration [kg/kg]
4385  real(RP), intent(in) :: dens0(ka,ia,ja) ! density [kg/m3]
4386  real(RP), intent(in) :: temp0(ka,ia,ja) ! temperature [K]
4387 
4388  ! mass concentration[kg/m3] and mean particle mass[kg]
4389  real(RP) :: xc(ka,ia,ja)
4390  real(RP) :: xr(ka,ia,ja)
4391  real(RP) :: xi(ka,ia,ja)
4392  real(RP) :: xs(ka,ia,ja)
4393  real(RP) :: xg(ka,ia,ja)
4394  ! diameter of average mass[kg/m3]
4395  real(RP) :: dc_ave(ka,ia,ja)
4396  real(RP) :: dr_ave(ka,ia,ja)
4397  ! radius of average mass
4398  real(RP) :: rc, rr
4399  ! 2nd. and 3rd. order moment of DSD
4400  real(RP) :: ri2m(ka,ia,ja), ri3m(ka,ia,ja)
4401  real(RP) :: rs2m(ka,ia,ja), rs3m(ka,ia,ja)
4402  real(RP) :: rg2m(ka,ia,ja), rg3m(ka,ia,ja)
4403 
4404  real(RP) :: coef_fuetal1998
4405  ! r2m_min is minimum value(moment of 1 particle with 1 micron)
4406  real(RP), parameter :: r2m_min=1.e-12_rp
4407  real(RP), parameter :: um2cm = 100.0_rp
4408 
4409  real(RP) :: limitsw, zerosw
4410  integer :: k, i, j
4411  !---------------------------------------------------------------------------
4412 
4413  ! mean particle mass[kg]
4414  do j = js, je
4415  do i = is, ie
4416  do k = ks, ke
4417  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) ) )
4418  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) ) )
4419  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) ) )
4420  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) ) )
4421  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) ) )
4422  enddo
4423  enddo
4424  enddo
4425 
4426  ! diameter of average mass : SB06 eq.(32)
4427  do j = js, je
4428  do i = is, ie
4429  do k = ks, ke
4430  dc_ave(k,i,j) = a_m(i_qc) * xc(k,i,j)**b_m(i_qc)
4431  dr_ave(k,i,j) = a_m(i_qr) * xr(k,i,j)**b_m(i_qr)
4432  enddo
4433  enddo
4434  enddo
4435 
4436  ! cloud effective radius
4437  do j = js, je
4438  do i = is, ie
4439  do k = ks, ke
4440  rc = 0.5_rp * dc_ave(k,i,j)
4441  limitsw = 0.5_rp + sign(0.5_rp, rc-rmin_re )
4442  re(k,i,j,i_mp_qc) = coef_re(i_qc) * rc * limitsw * um2cm
4443  enddo
4444  enddo
4445  enddo
4446 
4447  ! rain effective radius
4448  do j = js, je
4449  do i = is, ie
4450  do k = ks, ke
4451  rr = 0.5_rp * dr_ave(k,i,j)
4452  limitsw = 0.5_rp + sign(0.5_rp, rr-rmin_re )
4453  re(k,i,j,i_mp_qr) = coef_re(i_qr) * rr * limitsw * um2cm
4454  enddo
4455  enddo
4456  enddo
4457 
4458  do j = js, je
4459  do i = is, ie
4460  do k = ks, ke
4461  ri2m(k,i,j) = pi * coef_rea2(i_qi) * qtrc0(k,i,j,i_ni) * a_rea2(i_qi) * xi(k,i,j)**b_rea2(i_qi)
4462  rs2m(k,i,j) = pi * coef_rea2(i_qs) * qtrc0(k,i,j,i_ns) * a_rea2(i_qs) * xs(k,i,j)**b_rea2(i_qs)
4463  rg2m(k,i,j) = pi * coef_rea2(i_qg) * qtrc0(k,i,j,i_ng) * a_rea2(i_qg) * xg(k,i,j)**b_rea2(i_qg)
4464  enddo
4465  enddo
4466  enddo
4467 
4468  ! Fu(1996), eq.(3.11) or Fu et al.(1998), eq.(2.5)
4469  coef_fuetal1998 = 3.0_rp / (4.0_rp*rhoi)
4470  do j = js, je
4471  do i = is, ie
4472  do k = ks, ke
4473  ri3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ni) * xi(k,i,j)
4474  rs3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ns) * xs(k,i,j)
4475  rg3m(k,i,j) = coef_fuetal1998 * qtrc0(k,i,j,i_ng) * xg(k,i,j)
4476  enddo
4477  enddo
4478  enddo
4479 
4480  ! ice effective radius
4481  do j = js, je
4482  do i = is, ie
4483  do k = ks, ke
4484  zerosw = 0.5_rp - sign(0.5_rp, ri2m(k,i,j) - r2m_min )
4485  re(k,i,j,i_mp_qi) = ri3m(k,i,j) / ( ri2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
4486  enddo
4487  enddo
4488  enddo
4489 
4490  ! snow effective radius
4491  do j = js, je
4492  do i = is, ie
4493  do k = ks, ke
4494  zerosw = 0.5_rp - sign(0.5_rp, rs2m(k,i,j) - r2m_min )
4495  re(k,i,j,i_mp_qs) = rs3m(k,i,j) / ( rs2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
4496  enddo
4497  enddo
4498  enddo
4499 
4500  ! graupel effective radius
4501  do j = js, je
4502  do i = is, ie
4503  do k = ks, ke
4504  zerosw = 0.5_rp - sign(0.5_rp, rg2m(k,i,j) - r2m_min )
4505  re(k,i,j,i_mp_qg) = rg3m(k,i,j) / ( rg2m(k,i,j) + zerosw ) * ( 1.0_rp - zerosw ) * um2cm
4506  enddo
4507  enddo
4508  enddo
4509 
4510  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:18
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_tracer::mp_qa
integer, public mp_qa
Definition: scale_tracer.F90:98
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the caller graph for this function:

◆ atmos_phy_mp_sn14_mixingratio()

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

Calculate mixing ratio of each category.

Definition at line 4517 of file scale_atmos_phy_mp_sn14.F90.

References scale_tracer::mp_qa, and scale_tracer::qa.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_setup().

4517  use scale_grid_index
4518  use scale_tracer, only: &
4519  qad => qa, &
4520  mp_qad => mp_qa
4521  implicit none
4522 
4523  real(RP), intent(out) :: qe (ka,ia,ja,mp_qad) ! mixing ratio of each cateory [kg/kg]
4524  real(RP), intent(in) :: qtrc0(ka,ia,ja,qad) ! tracer mass concentration [kg/kg]
4525 
4526  integer :: ihydro
4527  !---------------------------------------------------------------------------
4528 
4529  do ihydro = 1, mp_qa
4530  qe(:,:,:,ihydro) = qtrc0(:,:,:,i_mp2all(ihydro))
4531  enddo
4532 
4533  return
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:44
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:18
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_tracer::mp_qa
integer, public mp_qa
Definition: scale_tracer.F90:98
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the caller graph for this function:

Variable Documentation

◆ atmos_phy_mp_dens

real(rp), dimension(mp_qa), target, public scale_atmos_phy_mp_sn14::atmos_phy_mp_dens

Definition at line 119 of file scale_atmos_phy_mp_sn14.F90.

Referenced by atmos_phy_mp_sn14_setup().

119  real(RP), public, target :: atmos_phy_mp_dens(mp_qa) ! hydrometeor density [kg/m3]=[g/L]
scale_tracer::mp_qa
integer, public mp_qa
Definition: scale_tracer.F90:98
scale_atmos_phy_mp_dry::atmos_phy_mp_dens
real(rp), dimension(mp_qa), target, public atmos_phy_mp_dens
Definition: scale_atmos_phy_mp_dry.F90:44
Generated by   doxygen 1.8.13