module Spectran Bin Microphysics More...

Functions/Subroutines
subroutine, public	atmos_phy_mp_suzuki10_tracer_setup
	Config. More...

subroutine, public	atmos_phy_mp_suzuki10_setup (KA, IA, JA)
	Setup. More...

subroutine, public	atmos_phy_mp_suzuki10_tendency (KA, KS, KE, IA, IS, IE, JA, JS, JE, KIJMAX, dt, DENS, PRES, TEMP, QTRC, QDRY, CPtot, CVtot, CCN, RHOQ_t, RHOE_t, CPtot_t, CVtot_t, EVAPORATE)
	Cloud Microphysics. More...

subroutine, public	atmos_phy_mp_suzuki10_terminal_velocity (KA, vterm_o)
	get terminal velocity More...

subroutine, public	atmos_phy_mp_suzuki10_cloud_fraction (KA, KS, KE, IA, IS, IE, JA, JS, JE, QTRC0, mask_criterion, cldfrac)
	Calculate Cloud Fraction. More...

subroutine, public	atmos_phy_mp_suzuki10_effective_radius (KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS0, TEMP0, QTRC0, Re)
	Calculate Effective Radius. More...

subroutine, public	atmos_phy_mp_suzuki10_qtrc2qhyd (KA, KS, KE, IA, IS, IE, JA, JS, JE, QTRC0, Qe)
	Calculate mass ratio of each category. More...

subroutine, public	atmos_phy_mp_suzuki10_qtrc2nhyd (KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, QTRC0, Ne)
	Calculate number concentration of each category. More...

subroutine, public	atmos_phy_mp_suzuki10_qhyd2qtrc (KA, KS, KE, IA, IS, IE, JA, JS, JE, Qe, QTRC, QNUM)
	get mass ratio of each category More...

Variables
integer, public	atmos_phy_mp_suzuki10_ntracers

integer, public	atmos_phy_mp_suzuki10_nwaters

integer, public	atmos_phy_mp_suzuki10_nices

integer, public	atmos_phy_mp_suzuki10_nccn

character(len=h_short), dimension(:), allocatable, public	atmos_phy_mp_suzuki10_tracer_names

character(len=h_mid), dimension(:), allocatable, public	atmos_phy_mp_suzuki10_tracer_descriptions

character(len=h_short), dimension(:), allocatable, public	atmos_phy_mp_suzuki10_tracer_units

integer, public	nbin = 33

integer, public	nccn = 0

Detailed Description

module Spectran Bin Microphysics

Description:: This module contains subroutines for the Spectral Bin Model

Reference
- Suzuki et al., 2006 Correlation Pattern between Effective Radius and Optical Thickness of Water Clouds Simulated by a Spectral Bin Microphysics Cloud Model SOLA, 2: 116-119 doi:10.2151/sola.2006-030
- Suzuki et al., 2010 A Study of Microphysical Mechanisms for Correlation Patterns between Droplet Radius and Optical Thickness of Warm Clouds with a Spectral Bin J. Atmos. Sci., 67: 1126-1141
- Sato et al., 2009 Application of a Monte Carlo integration method to collision and coagulation growth processes of hydrometeors in a bin-type model J. Geophy. Res., 114: D09215, doi:10.1029/2008JD011247

Author: : Team SCALE

NAMELIST

PARAM_ATMOS_PHY_MP_SUZUKI10_bin

name type default value comment

NBIN integer 33 tentatively public

NCCN integer 0 tentatively public

ICEFLG integer 1

KPHASE integer 0

PARAM_ATMOS_PHY_MP_SUZUKI10

name	type	default value	comment
MP_DOAUTOCONVERSION	logical	.true.	apply collision process ?
MP_COUPLE_AEROSOL	logical	.false.	apply CCN effect?
RHO_AERO	real(RP)		— density of aerosol
R_MIN	real(RP)		— minimum radius of aerosol (um)
R_MAX	real(RP)		— maximum radius of aerosol (um)
R0_AERO	real(RP)		— center radius of aerosol (um)
S10_EMAER	real(RP)		— moleculer weight of aerosol
S10_FLAG_REGENE	logical		— flag of regeneration
S10_FLAG_NUCLEAT	logical		— flag of regeneration
S10_FLAG_ICENUCLEAT	logical		— flag of regeneration
S10_FLAG_SFAERO	logical		— flag of surface flux of aeorol
S10_RNDM_FLGP	integer		— flag of surface flux of aeorol
S10_RNDM_MSPC	integer
S10_RNDM_MBIN	integer
C_CCN	real(RP)	100.E+6_RP	N0 of Nc = N0*s^kappa
KAPPA	real(RP)	0.462_RP	kappa of Nc = N0*s^kappa
SIGMA	real(RP)	7.5E-02_RP	water surface tension N/m2
VHFCT	real(RP)	2.0_RP	van't hoff factor (i in eq.(A.11) of Suzuki (2004))

History Output: No history output

Function/Subroutine Documentation

◆ atmos_phy_mp_suzuki10_tracer_setup()

subroutine, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_tracer_setup ( )

Config.

Definition at line 287 of file scale_atmos_phy_mp_suzuki10.F90.

References atmos_phy_mp_suzuki10_nccn, atmos_phy_mp_suzuki10_nices, atmos_phy_mp_suzuki10_ntracers, atmos_phy_mp_suzuki10_nwaters, atmos_phy_mp_suzuki10_tracer_descriptions, atmos_phy_mp_suzuki10_tracer_names, atmos_phy_mp_suzuki10_tracer_units, scale_io::io_fid_conf, nbin, nccn, and scale_prc::prc_abort().

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_tracer_setup().

     use scale_prc, only: &
        prc_abort
     implicit none
 
     namelist / param_atmos_phy_mp_suzuki10_bin / &
        nbin, &
        nccn, &
        iceflg, &
        kphase
 
     integer :: m, n, ierr
     !---------------------------------------------------------------------------
 
     log_newline
     log_info("ATMOS_PHY_MP_suzuki10_tracer_setup",*) 'Setup'
     log_info("ATMOS_PHY_MP_suzuki10_tracer_setup",*) 'Tracers setup for Suzuki (2010) Spectral BIN model'
 
     log_newline
     log_info("ATMOS_PHY_MP_suzuki10_tracer_setup",*) 'READ BIN NUMBER'
 
     rewind(io_fid_conf)
     read(io_fid_conf,nml=param_atmos_phy_mp_suzuki10_bin,iostat=ierr)
 
     if( ierr < 0 ) then !--- missing
       log_info("ATMOS_PHY_MP_suzuki10_tracer_setup",*)  'Not found namelist. Default used.'
     elseif( ierr > 0 ) then !--- fatal error
        log_error("ATMOS_PHY_MP_suzuki10_tracer_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SUZUKI10_bin, Check!'
        call prc_abort
     end if
 
     log_nml(param_atmos_phy_mp_suzuki10_bin)
 
     if( iceflg == 0 ) then
        nspc = 1
     elseif( iceflg == 1 ) then
        nspc = 7
     else
        log_error("ATMOS_PHY_MP_suzuki10_tracer_setup",*) "ICEFLG should be 0 (warm rain) or 1 (mixed rain) check!!"
        call prc_abort
     endif
 
     atmos_phy_mp_suzuki10_ntracers = 1 + nbin*nspc + nccn  ! number of total tracers
     atmos_phy_mp_suzuki10_nwaters  = nbin                  ! number of liquid water
     atmos_phy_mp_suzuki10_nices    = nbin * ( nspc - 1 )   ! number of ice water
     atmos_phy_mp_suzuki10_nccn     = nccn                  ! number of ccn
 
     num_start_waters = i_qv + 1
     num_end_waters   = i_qv + atmos_phy_mp_suzuki10_nwaters
     num_start_ices   = num_end_waters + 1
     num_end_ices     = num_end_waters + atmos_phy_mp_suzuki10_nices
 
     qa = atmos_phy_mp_suzuki10_ntracers
 
     allocate( atmos_phy_mp_suzuki10_tracer_names(qa) )
     allocate( atmos_phy_mp_suzuki10_tracer_descriptions(qa) )
     allocate( atmos_phy_mp_suzuki10_tracer_units(qa) )
 
     !---------------------------------------------------------------------------
     !
     !++ calculate each category and aerosol
     !
     !---------------------------------------------------------------------------
 
     do n = 1, qa
        write(atmos_phy_mp_suzuki10_tracer_units(n),'(a)')  'kg/kg'
     enddo
 
     write(atmos_phy_mp_suzuki10_tracer_names(1),'(a)') 'QV'
     write(atmos_phy_mp_suzuki10_tracer_descriptions(1),'(a)') 'Water Vapor mixing ratio'
 
     do m = 1, nspc
     do n = 1, nbin
        write(atmos_phy_mp_suzuki10_tracer_names(1+nbin*(m-1)+n),'(a,i0)') trim(namspc(m)), n
        write(atmos_phy_mp_suzuki10_tracer_descriptions(1+nbin*(m-1)+n),'(a,i0)') trim(lnamspc(m)), n
     enddo
     enddo
 
     do n = 1, nccn
        write(atmos_phy_mp_suzuki10_tracer_names(1+nbin*nspc+n),'(a,i0)') trim(namspc(8)), n
        write(atmos_phy_mp_suzuki10_tracer_descriptions(1+nbin*nspc+n),'(a,i0)') trim(lnamspc(8)), n
     enddo
 
     return

Here is the call graph for this function:

Here is the caller graph for this function:

◆ atmos_phy_mp_suzuki10_setup()

subroutine, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_setup	(	integer, intent(in)	KA,
		integer, intent(in)	IA,
		integer, intent(in)	JA
	)

Setup.

Definition at line 376 of file scale_atmos_phy_mp_suzuki10.F90.

References scale_comm_cartesc::comm_datatype, scale_comm_cartesc::comm_world, scale_const::const_dice, scale_const::const_dwatr, scale_const::const_pi, scale_atmos_hydrometeor::i_hc, scale_atmos_hydrometeor::i_hg, scale_atmos_hydrometeor::i_hh, scale_atmos_hydrometeor::i_hi, scale_atmos_hydrometeor::i_hr, scale_atmos_hydrometeor::i_hs, scale_io::io_fid_conf, scale_io::io_get_available_fid(), nbin, nccn, scale_prc::prc_abort(), scale_prc::prc_ismaster, and scale_prc::prc_masterrank.

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_setup().

     use scale_prc, only: &
        prc_abort,    &
        prc_masterrank, &
        prc_ismaster
     use scale_const, only: &
        pi => const_pi,       &
        dwatr => const_dwatr, &
        dice => const_dice
     use scale_comm_cartesc, only: &
        comm_world,    &
        comm_datatype
     use scale_atmos_hydrometeor, only: &
        i_hc, &
        i_hr, &
        i_hi, &
        i_hs, &
        i_hg, &
        i_hh
     implicit none
 
     integer, intent(in) :: ka
     integer, intent(in) :: ia
     integer, intent(in) :: ja
 
     real(RP) :: rho_aero  !--- density of aerosol
     real(RP) :: r0_aero   !--- center radius of aerosol (um)
     real(RP) :: r_min     !--- minimum radius of aerosol (um)
     real(RP) :: r_max     !--- maximum radius of aerosol (um)
     real(RP) :: s10_emaer !--- moleculer weight of aerosol
 
     logical :: s10_flag_regene  !--- flag of regeneration
     logical :: s10_flag_nucleat !--- flag of regeneration
     logical :: s10_flag_icenucleat !--- flag of regeneration
     logical :: s10_flag_sfaero  !--- flag of surface flux of aeorol
     integer :: s10_rndm_flgp  !--- flag of surface flux of aeorol
     integer :: s10_rndm_mspc
     integer :: s10_rndm_mbin
 
     namelist / param_atmos_phy_mp_suzuki10 / &
        mp_doautoconversion,    &
        mp_couple_aerosol,      &
        rho_aero,  &
        r_min, &
        r_max, &
        r0_aero,   &
        s10_emaer, &
        s10_flag_regene,  &
        s10_flag_nucleat, &
        s10_flag_icenucleat, &
        s10_flag_sfaero,  &
        s10_rndm_flgp, &
        s10_rndm_mspc, &
        s10_rndm_mbin, &
        c_ccn, kappa, &
        sigma, vhfct
 
     real(RP), parameter :: max_term_vel = 10.0_rp !-- terminal velocity for calculate dt of sedimentation
 
     integer :: nnspc, nnbin
     integer :: nn, mm, mmyu, nnyu
     integer :: myu, nyu, i, j, k, n, ierr
     !---------------------------------------------------------------------------
 
     log_newline
     log_info("ATMOS_PHY_MP_suzuki10_setup",*) 'Setup'
     log_info("ATMOS_PHY_MP_suzuki10_setup",*) 'Suzuki (2010) Spectral BIN model'
 
     !--- allocation
     allocate( xctr( nbin ) )
     allocate( xbnd( nbin+1 ) )
     allocate( radc( nbin ) )
     allocate( cctr( nbin,nspc_mk ) )
     allocate( cbnd( nbin+1,nspc_mk ) )
     allocate( ck( nspc_mk,nspc_mk,nbin,nbin ) )
     allocate( vt( nspc_mk,nbin ) )
     allocate( br( nspc_mk,nbin ) )
     allocate( ifrsl( 2,nspc_mk,nspc_mk ) )
     allocate( expxctr( nbin ) )
     allocate( expxbnd( nbin+1 ) )
     allocate( rexpxctr( nbin ) )
     allocate( rexpxbnd( nbin+1 ) )
     if ( nccn /= 0 ) then
       allocate( xactr( nccn ) )
       allocate( xabnd( nccn+1 ) )
       allocate( rada( nccn ) )
       allocate( expxactr( nccn ) )
       allocate( expxabnd( nccn+1 ) )
       allocate( rexpxactr( nccn ) )
       allocate( rexpxabnd( nccn+1 ) )
     endif
 
     mbin = nbin/2
     mspc = nspc_mk*nspc_mk
 
     rho_aero = rhoa
     s10_emaer = emaer
     r_min = rasta
     r_max = raend
     r0_aero = r0a
     s10_flag_regene = flg_regeneration
     s10_flag_nucleat = flg_nucl
     s10_flag_icenucleat = flg_icenucl
     s10_flag_sfaero = flg_sf_aero
     s10_rndm_flgp = rndm_flgp
     s10_rndm_mspc = mspc
     s10_rndm_mbin = mbin
 
     rewind(io_fid_conf)
     read(io_fid_conf,nml=param_atmos_phy_mp_suzuki10,iostat=ierr)
 
     if ( ierr < 0 ) then !--- missing
      log_info("ATMOS_PHY_MP_suzuki10_setup",*)  'Not found namelist. Default used.'
     elseif( ierr > 0 ) then !--- fatal error
      log_error("ATMOS_PHY_MP_suzuki10_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_SUZUKI10, Check!'
      call prc_abort
     endif
     log_nml(param_atmos_phy_mp_suzuki10)
 
     if ( nspc /= 1 .AND. nspc /= 7 ) then
        log_error("ATMOS_PHY_MP_suzuki10_setup",*) 'nspc should be set as 1 (warm rain) or 7 (mixed phase) check!'
        call prc_abort
     endif
 
     rhoa = rho_aero
     emaer = s10_emaer
     rasta = r_min
     raend = r_max
     r0a   = r0_aero
     flg_regeneration = s10_flag_regene
     flg_nucl = s10_flag_nucleat
     flg_icenucl = s10_flag_icenucleat
     flg_sf_aero = s10_flag_sfaero
     rndm_flgp = s10_rndm_flgp
     mspc = s10_rndm_mspc
     mbin = s10_rndm_mbin
 
     !--- read micpara.dat (microphysical parameter) and broad cast
     if ( prc_ismaster ) then
 
       fid_micpara = io_get_available_fid()
       !--- open parameter of cloud microphysics
       open ( fid_micpara, file = fname_micpara, form = 'formatted', status = 'old', iostat=ierr )
 
       !--- micpara.dat does not exist
       if ( ierr == 0 ) then
 
         read( fid_micpara,* ) nnspc, nnbin
 
         if ( nnbin /= nbin ) then
            log_error("ATMOS_PHY_MP_suzuki10_setup",*) 'nbin in namelist and nbin in micpara.dat is different check!'
            call prc_abort
         endif
 
         ! grid parameter
         log_info("ATMOS_PHY_MP_suzuki10_setup",*)  'Radius of cloud *'
         do n = 1, nbin
           read( fid_micpara,* ) nn, xctr( n ), radc( n )
           log_info("ATMOS_PHY_MP_suzuki10_setup",'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
                     "Radius of ", n, "th cloud bin (bin center)= ", radc( n ) , "[m]"
         enddo
         do n = 1, nbin+1
           read( fid_micpara,* ) nn, xbnd( n )
         enddo
         read( fid_micpara,* ) dxmic
         log_info("ATMOS_PHY_MP_suzuki10_setup",*)  'Width of Cloud SDF= ', dxmic
 
         ! capacity
         do myu = 1, nspc_mk
          do n = 1, nbin
 !          read( fid_micpara,* ) mmyu, nn, cctr( myu,n )
           read( fid_micpara,* ) mmyu, nn, cctr( n,myu )
          enddo
          do n = 1, nbin+1
 !          read( fid_micpara,* ) mmyu, nn, cbnd( myu,n )
           read( fid_micpara,* ) mmyu, nn, cbnd( n,myu )
          enddo
         enddo
 
         ! collection kernel
         do myu = 1, nspc_mk
          do nyu = 1, nspc_mk
           do i = 1, nbin
            do j = 1, nbin
             read( fid_micpara,* ) mmyu, nnyu, mm, nn, ck( myu,nyu,i,j )
            enddo
           enddo
          enddo
         enddo
 
         ! terminal velocity
         do myu = 1, nspc_mk
          do n = 1, nbin
           read( fid_micpara,* ) mmyu, nn, vt( myu,n )
          enddo
         enddo
 
         ! bulk density
         do myu = 1, nspc_mk
          do n = 1, nbin
           read( fid_micpara,* ) mmyu, nn, br( myu,n )
          enddo
         enddo
 
         close ( fid_micpara )
 
       !--- micpara.dat does not exist
       else
 
         log_info("ATMOS_PHY_MP_suzuki10_setup",*) 'micpara.dat is created'
         call mkpara
 
         fid_micpara = io_get_available_fid()
         !--- open parameter of cloud microphysics
         open ( fid_micpara, file = fname_micpara, form = 'formatted', status = 'old', iostat=ierr )
 
         read( fid_micpara,* ) nnspc, nnbin
 
         if ( nnbin /= nbin ) then
            log_error("ATMOS_PHY_MP_suzuki10_setup",*) 'nbin in inc_tracer and nbin in micpara.dat is different check!'
            call prc_abort
         endif
 
         ! grid parameter
         log_info("ATMOS_PHY_MP_suzuki10_setup",*)  'Radius of cloud *'
         do n = 1, nbin
           read( fid_micpara,* ) nn, xctr( n ), radc( n )
           log_info("ATMOS_PHY_MP_suzuki10_setup",'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
                     "Radius of ", n, "th cloud bin (bin center)= ", radc( n ) , "[m]"
         enddo
         do n = 1, nbin+1
           read( fid_micpara,* ) nn, xbnd( n )
         enddo
         read( fid_micpara,* ) dxmic
         log_info("ATMOS_PHY_MP_suzuki10_setup",*)  'Width of Cloud SDF= ', dxmic
 
         ! capacity
         do myu = 1, nspc_mk
          do n = 1, nbin
 !          read( fid_micpara,* ) mmyu, nn, cctr( myu,n )
           read( fid_micpara,* ) mmyu, nn, cctr( n,myu )
          enddo
          do n = 1, nbin+1
 !          read( fid_micpara,* ) mmyu, nn, cbnd( myu,n )
           read( fid_micpara,* ) mmyu, nn, cbnd( n,myu )
          enddo
         enddo
 
         ! collection kernel
         do myu = 1, nspc_mk
          do nyu = 1, nspc_mk
           do i = 1, nbin
            do j = 1, nbin
             read( fid_micpara,* ) mmyu, nnyu, mm, nn, ck( myu,nyu,i,j )
            enddo
           enddo
          enddo
         enddo
 
         ! terminal velocity
         do myu = 1, nspc_mk
          do n = 1, nbin
           read( fid_micpara,* ) mmyu, nn, vt( myu,n )
          enddo
         enddo
 
         ! bulk density
         do myu = 1, nspc_mk
          do n = 1, nbin
           read( fid_micpara,* ) mmyu, nn, br( myu,n )
          enddo
         enddo
 
         close ( fid_micpara )
 
       endif
 
     endif
 
     call mpi_bcast( radc, nbin,                      comm_datatype, prc_masterrank, comm_world, ierr )
     call mpi_bcast( xctr, nbin,                      comm_datatype, prc_masterrank, comm_world, ierr )
     call mpi_bcast( dxmic,1,                         comm_datatype, prc_masterrank, comm_world, ierr )
     call mpi_bcast( xbnd, nbin+1,                    comm_datatype, prc_masterrank, comm_world, ierr )
     call mpi_bcast( cctr, nbin*nspc_mk,              comm_datatype, prc_masterrank, comm_world, ierr )
     call mpi_bcast( cbnd, (nbin+1)*nspc_mk,          comm_datatype, prc_masterrank, comm_world, ierr )
     call mpi_bcast( ck,   nspc_mk*nspc_mk*nbin*nbin, comm_datatype, prc_masterrank, comm_world, ierr )
     call mpi_bcast( br, nbin*nspc_mk,                comm_datatype, prc_masterrank, comm_world, ierr )
     call mpi_bcast( vt, nbin*nspc_mk,                comm_datatype, prc_masterrank, comm_world, ierr )
 
     !--- aerosol ( CCN ) (not supported)
     if ( nccn /= 0 ) then
 
     allocate ( ncld( 1:nccn ) )
     xasta = log( rhoa*4.0_rp/3.0_rp*pi * ( rasta )**3 )
     xaend = log( rhoa*4.0_rp/3.0_rp*pi * ( raend )**3 )
 
     dxaer = ( xaend-xasta )/nccn
 
     do n = 1, nccn+1
      xabnd( n ) = xasta + dxaer*( n-1 )
     enddo
     do n = 1, nccn
      xactr( n ) = ( xabnd( n )+xabnd( n+1 ) )*0.50_rp
      rada( n )  = ( exp( xactr( n ) )*thirdovforth/pi/rhoa )**( oneovthird )
      log_info("ATMOS_PHY_MP_suzuki10_setup",'(A,1x,I3,1x,A,1x,ES15.7,1x,A)') &
           "Radius of ", n, "th aerosol bin (bin center)= ", rada( n ) , "[m]"
     enddo
 
     if ( flg_sf_aero ) then
       log_error("ATMOS_PHY_MP_suzuki10_setup",*) "flg_sf_aero=true is not supported stop!! "
       call prc_abort
 !     if ( CZ(KS) >= 10.0_RP ) then
 !          R10M1 = 10.0_RP / CZ(KS) * 0.50_RP ! scale with height
 !          R10M2 = 10.0_RP / CZ(KS) * 0.50_RP ! scale with height
 !          R10H1 = 1.0_RP * 0.50_RP
 !          R10H2 = 1.0_RP * 0.0_RP
 !          R10E1 = 1.0_RP * 0.50_RP
 !          R10E2 = 1.0_RP * 0.50_RP
 !          K10_1 = KS
 !          K10_2 = KS
 !     else
 !       k = 1
 !       do while ( CZ(k) < 10.0_RP )
 !          k = k + 1
 !          K10_1 = k
 !          K10_2 = k + 1
 !          R10M1 = ( CZ(k+1) - 10.0_RP ) / CDZ(k)
 !          R10M2 = ( 10.0_RP   - CZ(k) ) / CDZ(k)
 !          R10H1 = ( CZ(k+1) - 10.0_RP ) / CDZ(k)
 !          R10H2 = ( 10.0_RP   - CZ(k) ) / CDZ(k)
 !          R10E1 = ( CZ(k+1) - 10.0_RP ) / CDZ(k)
 !          R10E2 = ( 10.0_RP   - CZ(k) ) / CDZ(k)
 !       enddo
 !     endif
     endif
 
     endif
 
     !--- determine nbnd
     do n = 1, nbin
       if( radc( n ) > rbnd ) then
         nbnd = n
         exit
       endif
     enddo
     log_info("ATMOS_PHY_MP_suzuki10_setup",'(A,ES15.7,A)')  'Radius between cloud and rain is ', radc(nbnd), '[m]'
 
     !--- random number setup for stochastic method
     if ( rndm_flgp > 0 ) then
      call random_setup( ia*ja*ka )
     endif
 
     if ( mp_couple_aerosol .AND. nccn /=0 ) then
       log_error("ATMOS_PHY_MP_suzuki10_setup",*) 'nccn should be 0 when MP_couple_aerosol = .true. !! stop'
       call prc_abort
     endif
 
     if ( nccn /= 0 ) then
      do n = 1, nccn
       expxactr( n ) = exp( xactr( n ) )
       rexpxactr( n ) = 1.0_rp / exp( xactr( n ) )
      enddo
      do n = 1, nccn+1
       expxabnd( n ) = exp( xabnd( n ) )
       rexpxabnd( n ) = 1.0_rp / exp( xabnd( n ) )
      enddo
     endif
 
     allocate( vterm(qa-1) )
     vterm(:) = 0.0_rp
     do myu = 1, nspc
     do n = 1, nbin
       vterm((myu-1)*nbin+n) = -vt( myu,n )
     enddo
     enddo
     do n = 1, nbin
       expxctr( n ) = exp( xctr( n ) )
       rexpxctr( n ) = 1.0_rp / exp( xctr( n ) )
     enddo
     do n = 1, nbin+1
       expxbnd( n ) = exp( xbnd( n ) )
       rexpxbnd( n ) = 1.0_rp / exp( xbnd( n ) )
     enddo
 
     allocate( kindx(nbin,nbin) )
     call getrule( ifrsl,kindx )
 
     !--- determine the parameters for interpolating SDF from qxx, Nxx of parent domain
     sigma_sdf(1) = 0.2_rp
     sigma_sdf(2) = 0.35_rp
     sigma_sdf(3) = 0.35_rp
     sigma_sdf(4) = 0.35_rp
     sigma_sdf(5) = 0.35_rp
     r0_sdf(1)    = 5.e-6_rp
     r0_sdf(2)    = 2.61e-6_rp
     r0_sdf(3)    = 5.e-6_rp
     r0_sdf(4)    = 2.61e-6_rp
     r0_sdf(5)    = 2.61e-6_rp ! to be corrected
     n0_sdf(1)    = 8.0e+6_rp
     n0_sdf(2)    = 0.0_rp
     n0_sdf(3)    = 3.0e+6_rp
     n0_sdf(4)    = 4.0e+6_rp
     n0_sdf(5)    = 4.0e+6_rp ! to be corrected
     rho_sdf(1)   = dwatr
     rho_sdf(2)   = dice
     rho_sdf(3)   = 100.0_rp
     rho_sdf(4)   = 400.0_rp
     rho_sdf(5)   = 400.0_rp ! to be corrected
 
     !--- determine the parameters for interpolating SDF from qxx, Nxx of parent domain
     sigma_sdf(1) = 0.2_rp
     sigma_sdf(2) = 0.35_rp
     sigma_sdf(3) = 0.35_rp
     sigma_sdf(4) = 0.35_rp
     sigma_sdf(5) = 0.35_rp
     r0_sdf(1)    = 5.e-6_rp
     r0_sdf(2)    = 2.61e-6_rp
     r0_sdf(3)    = 5.e-6_rp
     r0_sdf(4)    = 2.61e-6_rp
     r0_sdf(5)    = 2.61e-6_rp ! to be corrected
     n0_sdf(1)    = 8.0e+6_rp
     n0_sdf(2)    = 0.0_rp
     n0_sdf(3)    = 3.0e+6_rp
     n0_sdf(4)    = 4.0e+6_rp
     n0_sdf(5)    = 4.0e+6_rp ! to be corrected
     rho_sdf(1)   = dwatr
     rho_sdf(2)   = dice
     rho_sdf(3)   = 100.0_rp
     rho_sdf(4)   = 400.0_rp
     rho_sdf(5)   = 400.0_rp ! to be corrected
 
     return

Here is the call graph for this function:

Here is the caller graph for this function:

◆ atmos_phy_mp_suzuki10_tendency()

subroutine, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_tendency	(	integer, intent(in)	KA,
		integer, intent(in)	KS,
		integer, intent(in)	KE,
		integer, intent(in)	IA,
		integer, intent(in)	IS,
		integer, intent(in)	IE,
		integer, intent(in)	JA,
		integer, intent(in)	JS,
		integer, intent(in)	JE,
		integer, intent(in)	KIJMAX,
		real(dp), intent(in)	dt,
		real(rp), dimension (ka,ia,ja), intent(in)	DENS,
		real(rp), dimension (ka,ia,ja), intent(in)	PRES,
		real(rp), dimension (ka,ia,ja), intent(in)	TEMP,
		real(rp), dimension (ka,ia,ja,qa), intent(in)	QTRC,
		real(rp), dimension (ka,ia,ja), intent(in)	QDRY,
		real(rp), dimension(ka,ia,ja), intent(in)	CPtot,
		real(rp), dimension(ka,ia,ja), intent(in)	CVtot,
		real(rp), dimension (ka,ia,ja), intent(in)	CCN,
		real(rp), dimension (ka,ia,ja,qa), intent(out)	RHOQ_t,
		real(rp), dimension (ka,ia,ja), intent(out)	RHOE_t,
		real(rp), dimension(ka,ia,ja), intent(out)	CPtot_t,
		real(rp), dimension(ka,ia,ja), intent(out)	CVtot_t,
		real(rp), dimension(ka,ia,ja), intent(out)	EVAPORATE
	)

Cloud Microphysics.

Definition at line 822 of file scale_atmos_phy_mp_suzuki10.F90.

References scale_const::const_tem00, nbin, nccn, scale_prof::prof_rapend(), and scale_prof::prof_rapstart().

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_calc_tendency().

     use scale_const, only: &
        tem00 => const_tem00
     use scale_atmos_saturation, only: &
        atmos_saturation_pres2qsat_liq, &
        atmos_saturation_pres2qsat_ice
     implicit none
 
     integer, intent(in) :: ka, ks, ke
     integer, intent(in) :: ia, is, ie
     integer, intent(in) :: ja, js, je
     integer, intent(in) :: kijmax
 
     real(DP), intent(in) :: dt
     real(RP), intent(in) :: dens (ka,ia,ja)
     real(RP), intent(in) :: pres (ka,ia,ja)
     real(RP), intent(in) :: temp (ka,ia,ja)
     real(RP), intent(in) :: qtrc (ka,ia,ja,qa)
     real(RP), intent(in) :: qdry (ka,ia,ja)
     real(RP), intent(in) :: cptot(ka,ia,ja)
     real(RP), intent(in) :: cvtot(ka,ia,ja)
     real(RP), intent(in) :: ccn  (ka,ia,ja)
 
     real(RP), intent(out) :: rhoq_t (ka,ia,ja,qa)
     real(RP), intent(out) :: rhoe_t (ka,ia,ja)
     real(RP), intent(out) :: cptot_t(ka,ia,ja)
     real(RP), intent(out) :: cvtot_t(ka,ia,ja)
     real(RP), intent(out) :: evaporate(ka,ia,ja)   !--- number of evaporated cloud [/m3]
 
     real(RP) :: qsat_l(ka,ia,ja)
     real(RP) :: qsat_i(ka,ia,ja)
     real(RP) :: ssliq(ka,ia,ja)
     real(RP) :: ssice(ka,ia,ja)
 
     integer  :: ijk_index (kijmax,3)
     integer  :: index_cld (kijmax)
     integer  :: index_cold(kijmax)
     integer  :: index_warm(kijmax)
     integer  :: ijkcount, ijkcount_cold, ijkcount_warm
     integer  :: ijk, indirect
 
     real(RP) :: dens_ijk(kijmax)
     real(RP) :: pres_ijk(kijmax)
     real(RP) :: temp_ijk(kijmax)
     real(RP) :: qdry_ijk(kijmax)
     real(RP) :: qvap_ijk(kijmax)
     real(RP) :: ccn_ijk(kijmax)
     real(RP) :: cp_ijk(kijmax)
     real(RP) :: cv_ijk(kijmax)
     real(RP) :: evaporate_ijk(kijmax)
     real(RP) :: ghyd_ijk(nbin,nspc,kijmax)
     real(RP) :: gaer_ijk(max(nccn,1),kijmax)
     real(RP) :: cldsum
     integer  :: countbin
     real(RP) :: rhoq_new
 
     integer  :: step
     integer  :: k, i, j, m, n, iq
     !---------------------------------------------------------------------------
 
     if    ( nspc == 1 ) then
        log_progress(*) 'atmosphere / physics / microphysics / SBM (Liquid water only)'
     elseif( nspc >  1 ) then
        log_progress(*) 'atmosphere / physics / microphysics / SBM (Mixed phase)'
     endif
 
     call atmos_saturation_pres2qsat_liq( ka, ks, ke, & ! [IN]
                                          ia, is, ie, & ! [IN]
                                          ja, js, je, & ! [IN]
                                          temp(:,:,:), pres(:,:,:), qdry(:,:,:), & ! [IN]
                                          qsat_l(:,:,:)                          ) ! [OUT]
 
     call atmos_saturation_pres2qsat_ice( ka, ks, ke, & ! [IN]
                                          ia, is, ie, & ! [IN]
                                          ja, js, je, & ! [IN]
                                          temp(:,:,:), pres(:,:,:), qdry(:,:,:), & ! [IN]
                                          qsat_i(:,:,:)                          ) ! [OUT]
 
     do j = js, je
     do i = is, ie
     do k = ks, ke
        ssliq(k,i,j) = qtrc(k,i,j,i_qv) / qsat_l(k,i,j) - 1.0_rp
        ssice(k,i,j) = qtrc(k,i,j,i_qv) / qsat_i(k,i,j) - 1.0_rp
     enddo
     enddo
     enddo
 
     if ( nspc == 1 ) then
        ssice(:,:,:) = 0.0_rp
     endif
 
 !--- store initial SDF of aerosol
 !--- this option is not supported
 !    if ( ofirst_sdfa ) then
 !      allocate( marate( nccn ) )
 !      do j = JS, JE
 !      do i = IS, IE
 !         do k = KS, KE
 !           sum2 = 0.0_RP
 !           do n = 1, nccn
 !             marate( n ) = gdga(k,i,j,n)*rexpxactr( n )
 !             sum2 = sum2 + gdga(k,i,j,n)*rexpxactr( n )
 !           enddo
 !         enddo
 !      enddo
 !      enddo
 !      if ( sum2 /= 0.0_RP ) then
 !        marate( 1:nccn ) = marate( 1:nccn )/sum2
 !        ofirst_sdfa = .false.
 !      endif
 !    endif
 
     !--- Arrange array for microphysics
 
     call prof_rapstart('MP_ijkconvert', 3)
 
     ijk = 0
     do j = js, je
     do i = is, ie
     do k = ks, ke
        ijk = ijk + 1
        ijk_index(ijk,1) = i
        ijk_index(ijk,2) = j
        ijk_index(ijk,3) = k
     enddo
     enddo
     enddo
 
     ijkcount = 0
     ijkcount_cold = 0
     ijkcount_warm = 0
 
     ijk = 0
     do j = js, je
     do i = is, ie
     do k = ks, ke
        ijk = ijk + 1
 
        ! calc total hydrometeors
        cldsum   = 0.0_rp
        countbin = i_qv + 1
        do m = 1, nspc
        do n = 1, nbin
          cldsum   = cldsum + qtrc(k,i,j,countbin) * dens(k,i,j) / dxmic
          countbin = countbin + 1
        enddo
        enddo
 
        if (      cldsum > cldmin       &
             .OR. ssliq(k,i,j) > 0.0_rp &
             .OR. ssice(k,i,j) > 0.0_rp ) then
 
           ijkcount = ijkcount + 1
 
           index_cld(ijkcount) = ijk
 
           dens_ijk(ijkcount) = dens(k,i,j)
           pres_ijk(ijkcount) = pres(k,i,j)
           temp_ijk(ijkcount) = temp(k,i,j)
           qdry_ijk(ijkcount) = qdry(k,i,j)
           cp_ijk(ijkcount) = cptot(k,i,j)
           cv_ijk(ijkcount) = cvtot(k,i,j)
           ccn_ijk(ijkcount)  = ccn(k,i,j)
           qvap_ijk(ijkcount) = qtrc(k,i,j,i_qv)
 
           countbin = i_qv + 1
           do m = 1, nspc
           do n = 1, nbin
              ghyd_ijk(n,m,ijkcount) = qtrc(k,i,j,countbin) * dens(k,i,j) / dxmic
              countbin = countbin + 1
           enddo
           enddo
 
           do n = 1, nccn
              gaer_ijk(n,ijkcount)   = qtrc(k,i,j,countbin) * dens(k,i,j) / dxaer
              countbin = countbin + 1
           enddo
 
           if ( temp(k,i,j) < tem00 .AND. nspc > 1 ) then ! cold
             ijkcount_cold = ijkcount_cold + 1
             index_cold(ijkcount_cold) = ijkcount
           else ! warm
             ijkcount_warm = ijkcount_warm + 1
             index_warm(ijkcount_warm) = ijkcount
           endif
 
        else
 
           ! no hudrometeors and undersaturation (no microphysical process occcurs)
           do iq = 1, qa
              rhoq_t(k,i,j,iq) = 0.0_rp
           end do
           rhoe_t(k,i,j) = 0.0_rp
           cptot_t(k,i,j) = 0.0_rp
           cvtot_t(k,i,j) = 0.0_rp
           evaporate(k,i,j) = 0.0_rp
 
        endif
 
     enddo
     enddo
     enddo
 
     call prof_rapend  ('MP_ijkconvert', 3)
 
     ! tentative timername registration
     call prof_rapstart('MP_suzuki10', 3)
     call prof_rapend  ('MP_suzuki10', 3)
     call prof_rapstart('_SBM_Nucleat',    3)
     call prof_rapend  ('_SBM_Nucleat',    3)
 !    call PROF_rapstart('_SBM_NucleatA',   3)
 !    call PROF_rapend  ('_SBM_NucleatA',   3)
     call prof_rapstart('_SBM_Liqphase',   3)
     call prof_rapend  ('_SBM_Liqphase',   3)
     call prof_rapstart('_SBM_Icephase',   3)
     call prof_rapend  ('_SBM_Icephase',   3)
     call prof_rapstart('_SBM_Mixphase',   3)
     call prof_rapend  ('_SBM_Mixphase',   3)
     call prof_rapstart('_SBM_AdvLiq',     3)
     call prof_rapend  ('_SBM_AdvLiq',     3)
     call prof_rapstart('_SBM_AdvIce',     3)
     call prof_rapend  ('_SBM_AdvIce',     3)
     call prof_rapstart('_SBM_AdvMix',     3)
     call prof_rapend  ('_SBM_AdvMix',     3)
 !    call PROF_rapstart('_SBM_FAero',      3)
 !    call PROF_rapend  ('_SBM_FAero',      3)
     call prof_rapstart('_SBM_Freezing',   3)
     call prof_rapend  ('_SBM_Freezing',   3)
     call prof_rapstart('_SBM_IceNucleat', 3)
     call prof_rapend  ('_SBM_IceNucleat', 3)
     call prof_rapstart('_SBM_Melting',    3)
     call prof_rapend  ('_SBM_Melting',    3)
     call prof_rapstart('_SBM_CollCoag',   3)
     call prof_rapend  ('_SBM_CollCoag',   3)
 !    call PROF_rapstart('_SBM_CollCoagR',  3)
 !    call PROF_rapend  ('_SBM_CollCoagR',  3)
 
     if ( ijkcount > 0 ) then
 
     call prof_rapstart('MP_suzuki10', 3)
 
     call mp_suzuki10( ka, ia, ja,                 & ! [IN]
                       ijkcount,                   & ! [IN]
                       ijkcount_cold,              & ! [IN]
                       ijkcount_warm,              & ! [IN]
                       index_cold(    1:ijkcount), & ! [IN]
                       index_warm(    1:ijkcount), & ! [IN]
                       dens_ijk(    1:ijkcount), & ! [IN]
                       pres_ijk(    1:ijkcount), & ! [IN]
                       qdry_ijk(    1:ijkcount), & ! [IN]
                       ccn_ijk(    1:ijkcount), & ! [IN]
                       temp_ijk(    1:ijkcount), & ! [INOUT]
                       qvap_ijk(    1:ijkcount), & ! [INOUT]
                       ghyd_ijk(:,:,1:ijkcount), & ! [INOUT]
                       gaer_ijk(:,  1:ijkcount), & ! [INOUT]
                       cp_ijk(    1:ijkcount), & ! [INOUT]
                       cv_ijk(    1:ijkcount), & ! [INOUT]
                       evaporate_ijk(1:ijkcount),  & ! [OUT]
                       dt                          ) ! [IN]
 
     call prof_rapend  ('MP_suzuki10', 3)
 
 !    if ( flg_sf_aero ) then
 !     do j = JS-2, JE+2
 !     do i = IS-2, IE+1
 !       VELX(i,j) = MOMX(K10_1,i,j) / ( DENS(K10_1,i+1,j)+DENS(K10_1,i,j) ) * R10M1 &
 !                 + MOMX(K10_2,i,j) / ( DENS(K10_2,i+1,j)+DENS(K10_2,i,j) ) * R10M2
 !     enddo
 !     enddo
 !
 !     do j = JS-2, JE+1
 !     do i = IS-2, IE+2
 !       VELY(i,j) = MOMY(K10_1,i,j) / ( DENS(K10_1,i,j+1)+DENS(K10_1,i,j) ) * R10M1 &
 !                 + MOMY(K10_2,i,j) / ( DENS(K10_2,i,j+1)+DENS(K10_2,i,j) ) * R10M2
 !     enddo
 !     enddo
 !    endif
 !
 !    !--- SURFACE FLUX by Monahan et al. (1986)
 !    if ( flg_sf_aero .AND. nccn /= 0 ) then
 !     do j = JS, JE
 !     do i = IS, IE
 !          ijk = ( j - JS ) * KMAX * IMAX &
 !              + ( i - IS ) * KMAX        &
 !              + ( KS - KS )              &
 !              + 1
 !       Uabs = sqrt(  ( ( VELX(i,j) + VELX(i-1,j  ) ) * 0.50_RP )**2 &
 !                   + ( ( VELY(i,j) + VELY(i  ,j-1) ) * 0.50_RP )**2 )
 !       do n = 1, nccn
 !        if ( rada( n ) <= 2.0E-5_RP .AND. rada( n ) >= 3.0E-7_RP ) then
 !         bparam = ( 0.38_RP - log( rada( n ) ) )/0.65_RP
 !         SFLX_AERO(i,j,n) = 1.373_RP * Uabs**( 3.41_RP ) * rada( n )**( -3.0_RP ) &
 !                          * ( 1.0_RP + 0.057_RP * rada( n )**( 1.05_RP ) ) &
 !                          * 10.0_RP**( 1.19_RP * exp( -bparam*bparam ) )
 !         ! convert from [#/m^2/um/s] -> [kg/m^3/unit log (m)]
 !         SFLX_AERO(i,j,n) = SFLX_AERO(i,j,n) / DENS(KS,i,j) &
 !                          / CDZ(KS) * rada( n ) / 3.0_RP * dt * expxactr( n )
 !         Gaer_ijk(n,ijk) = Gaer_ijk(n,ijk) + SFLX_AERO(i,j,n)/dxaer
 !        endif
 !       enddo
 !     enddo
 !     enddo
 !    endif
 
     call prof_rapstart('MP_ijkconvert', 3)
 
     !---- return original array
     do ijk = 1, ijkcount
        indirect = index_cld(ijk)
        i = ijk_index(indirect,1)
        j = ijk_index(indirect,2)
        k = ijk_index(indirect,3)
 
        rhoe_t(k,i,j) = ( temp_ijk(ijk) * cv_ijk(ijk) - temp(k,i,j) * cvtot(k,i,j) ) * dens(k,i,j) / dt
        cptot_t(k,i,j) = ( cp_ijk(ijk) - cptot(k,i,j) ) / dt
        cvtot_t(k,i,j) = ( cv_ijk(ijk) - cvtot(k,i,j) ) / dt
        evaporate(k,i,j) = evaporate_ijk(ijk) / dt ! [#/m3/s]
 
        rhoq_t(k,i,j,i_qv) = ( qvap_ijk(ijk) - qtrc(k,i,j,i_qv) ) * dens(k,i,j) / dt
 
        countbin = i_qv + 1
        do m = 1, nspc
        do n = 1, nbin
           rhoq_new = ghyd_ijk(n,m,ijk) * dxmic
           rhoq_t(k,i,j,countbin) = ( rhoq_new - qtrc(k,i,j,countbin)*dens(k,i,j) ) / dt
           countbin = countbin + 1
        enddo
        enddo
 
        do n = 1, nccn
           rhoq_new = gaer_ijk(n,ijk) * dxaer
           rhoq_t(k,i,j,countbin) = ( rhoq_new - qtrc(k,i,j,countbin)*dens(k,i,j) ) / dt
           countbin = countbin + 1
        enddo
     enddo
 
 
 !    if ( nccn /= 0 ) then
 !       AMR(:,:,:) = 0.0_RP
 !       do j = JS, JE
 !       do i = IS, IE
 !       do k = KS, KE
 !          do n = 1, nccn
 !             AMR(k,i,j) = AMR(k,i,j) + QTRC(k,i,j,QQE-1+n)
 !          enddo
 !       enddo
 !       enddo
 !       enddo
 !    endif
 
     call prof_rapend  ('MP_ijkconvert', 3)
 
     endif
 
     return

Here is the call graph for this function:

Here is the caller graph for this function:

◆ atmos_phy_mp_suzuki10_terminal_velocity()

subroutine, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_terminal_velocity	(	integer, intent(in)	KA,
		real(rp), dimension(ka,qa-1), intent(out)	vterm_o
	)

get terminal velocity

Definition at line 1184 of file scale_atmos_phy_mp_suzuki10.F90.

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_calc_tendency().

     implicit none
 
     integer, intent(in) :: ka
 
     real(RP), intent(out) :: vterm_o(ka,qa-1)
 
     integer :: iq
     !---------------------------------------------------------------------------
 
     do iq = 1, qa-1
        vterm_o(:,iq) = vterm(iq)
     end do
 
     return

Here is the caller graph for this function:

◆ atmos_phy_mp_suzuki10_cloud_fraction()

subroutine, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_cloud_fraction	(	integer, intent(in)	KA,
		integer, intent(in)	KS,
		integer, intent(in)	KE,
		integer, intent(in)	IA,
		integer, intent(in)	IS,
		integer, intent(in)	IE,
		integer, intent(in)	JA,
		integer, intent(in)	JS,
		integer, intent(in)	JE,
		real(rp), dimension (ka,ia,ja,nspc*nbin), intent(in)	QTRC0,
		real(rp), intent(in)	mask_criterion,
		real(rp), dimension(ka,ia,ja), intent(out)	cldfrac
	)

Calculate Cloud Fraction.

Definition at line 1209 of file scale_atmos_phy_mp_suzuki10.F90.

References nbin.

Referenced by mod_atmos_phy_mp_vars::atmos_phy_mp_vars_get_diagnostic().

     implicit none
 
     integer, intent(in) :: ka, ks, ke
     integer, intent(in) :: ia, is, ie
     integer, intent(in) :: ja, js, je
 
     real(RP), intent(in)  :: qtrc0  (ka,ia,ja,nspc*nbin)
     real(RP), intent(in)  :: mask_criterion
     real(RP), intent(out) :: cldfrac(ka,ia,ja)
 
     real(RP) :: qhydro
     integer  :: k, i, j, iq, ihydro
     !---------------------------------------------------------------------------
 
     if( nspc > 1 ) then
       do j  = js, je
       do i  = is, ie
       do k  = ks, ke
          qhydro = 0.0_rp
          do ihydro = 1, nspc
           do iq = nbin*(ihydro-1)+1, nbin*ihydro
             qhydro = qhydro + qtrc0(k,i,j,iq)
           enddo
          enddo
          cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
       enddo
       enddo
       enddo
     elseif( nspc == 1 ) then
       do j  = js, je
       do i  = is, ie
       do k  = ks, ke
          qhydro = 0.0_rp
          do ihydro = 1, i_mp_qc
           do iq = nbin*(ihydro-1)+1, nbin*ihydro
             qhydro = qhydro + qtrc0(k,i,j,iq)
           enddo
          enddo
          cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
       enddo
       enddo
       enddo
     endif
 
     return

Here is the caller graph for this function:

◆ atmos_phy_mp_suzuki10_effective_radius()

subroutine, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_effective_radius	(	integer, intent(in)	KA,
		integer, intent(in)	KS,
		integer, intent(in)	KE,
		integer, intent(in)	IA,
		integer, intent(in)	IS,
		integer, intent(in)	IE,
		integer, intent(in)	JA,
		integer, intent(in)	JS,
		integer, intent(in)	JE,
		real(rp), dimension(ka,ia,ja), intent(in)	DENS0,
		real(rp), dimension(ka,ia,ja), intent(in)	TEMP0,
		real(rp), dimension(ka,ia,ja,nspc*nbin), intent(in)	QTRC0,
		real(rp), dimension (ka,ia,ja,n_hyd), intent(out)	Re
	)

Calculate Effective Radius.

Definition at line 1266 of file scale_atmos_phy_mp_suzuki10.F90.

References scale_const::const_eps, scale_atmos_hydrometeor::i_hc, scale_atmos_hydrometeor::i_hg, scale_atmos_hydrometeor::i_hh, scale_atmos_hydrometeor::i_hi, scale_atmos_hydrometeor::i_hr, scale_atmos_hydrometeor::i_hs, scale_atmos_hydrometeor::i_qv, scale_atmos_hydrometeor::n_hyd, and nbin.

Referenced by mod_atmos_phy_mp_vars::atmos_phy_mp_vars_get_diagnostic().

     use scale_const, only: &
        eps => const_eps
     use scale_atmos_hydrometeor, only: &
        n_hyd, &
        i_qv,  &
        i_hc,  &
        i_hr,  &
        i_hi,  &
        i_hs,  &
        i_hg,  &
        i_hh
     implicit none
 
     integer, intent(in) :: ka, ks, ke
     integer, intent(in) :: ia, is, ie
     integer, intent(in) :: ja, js, je
 
     real(RP), intent(in)  :: dens0(ka,ia,ja)       ! density                   [kg/m3]
     real(RP), intent(in)  :: temp0(ka,ia,ja)       ! temperature               [K]
     real(RP), intent(in)  :: qtrc0(ka,ia,ja,nspc*nbin) ! tracer mass concentration [kg/kg]
     real(RP), intent(out) :: re   (ka,ia,ja,n_hyd) ! effective radius          [cm]
 
     real(RP), parameter :: um2cm = 100.0_rp
 
     real(RP) :: sum0(nspc), sum2, sum3, re_tmp(nspc)
     integer  :: i, j, k, iq, ihydro
     !---------------------------------------------------------------------------
 
     do k = ks, ke
     do j = js, je
     do i = is, ie
        re(k,i,j,:) = 0.0_rp
 
        ! HC
        sum3 = 0.0_rp
        sum2 = 0.0_rp
        ihydro = i_mp_qc
        do iq = 1, nbnd
           sum3 = sum3 &
                + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
                * rexpxctr( iq-nbin*(ihydro-1) ) &   !--- mass -> number
                * radc( iq-nbin*(ihydro-1) )**3.0_rp )
           sum2 = sum2 &
                + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
                * rexpxctr( iq-nbin*(ihydro-1) ) &   !--- mass -> number
                * radc( iq-nbin*(ihydro-1) )**2.0_rp )
        enddo
        sum3 = max( sum3, 0.0_rp )
        sum2 = max( sum2, 0.0_rp )
        if ( sum2 /= 0.0_rp ) then
           re(k,i,j,i_hc) = sum3 / sum2 * um2cm
        else
           re(k,i,j,i_hc) = 0.0_rp
        endif
 
        ! HR
        sum3 = 0.0_rp
        sum2 = 0.0_rp
        ihydro = i_mp_qc
        do iq = nbnd+1, nbin
           sum3 = sum3 &
                + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
                * rexpxctr( iq-nbin*(ihydro-1) ) &   !--- mass -> number
                * radc( iq-nbin*(ihydro-1) )**3.0_rp )
           sum2 = sum2 &
                + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
                * rexpxctr( iq-nbin*(ihydro-1) ) &   !--- mass -> number
                * radc( iq-nbin*(ihydro-1) )**2.0_rp )
        enddo
        sum3 = max( sum3, 0.0_rp )
        sum2 = max( sum2, 0.0_rp )
        if ( sum2 /= 0.0_rp ) then
           re(k,i,j,i_hr) = sum3 / sum2 * um2cm
        else
           re(k,i,j,i_hr) = 0.0_rp
        endif
 
     enddo
     enddo
     enddo
 
     ! other hydrometeors
     if ( nspc > 1 ) then
        do k = ks, ke
        do j = js, je
        do i = is, ie
           do ihydro = 2, nspc
              sum0(ihydro) = 0.0_rp
              sum2 = 0.0_rp
              sum3 = 0.0_rp
              do iq = nbin*(ihydro-1)+1, nbin*ihydro
                 sum0(ihydro) = sum0(ihydro) &
                      + ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) !--- [kg/kg] -> [kg/m3]
                 sum3 = sum3 &
                      + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
                      * rexpxctr( iq-nbin*(ihydro-1) ) &   !--- mass -> number
                      * radc( iq-nbin*(ihydro-1) )**3.0_rp )
                 sum2 = sum2 &
                      + ( ( qtrc0(k,i,j,iq) * dens0(k,i,j) ) & !--- [kg/kg] -> [kg/m3]
                      * rexpxctr( iq-nbin*(ihydro-1) ) &   !--- mass -> number
                      * radc( iq-nbin*(ihydro-1) )**2.0_rp )
              enddo
              sum3 = max( sum3, 0.0_rp )
              sum2 = max( sum2, 0.0_rp )
              if ( sum2 == 0.0_rp ) then
                 re_tmp(ihydro) = 0.0_rp
              else
                 re_tmp(ihydro) = sum3 / sum2 * um2cm
              end if
           end do
 
           re(k,i,j,i_hi) = ( re_tmp(i_mp_qcl) * sum0(i_mp_qcl) &
                            + re_tmp(i_mp_qp ) * sum0(i_mp_qp ) &
                            + re_tmp(i_mp_qd ) * sum0(i_mp_qd ) ) &
                          / ( sum0(i_mp_qcl) + sum0(i_mp_qp) + sum0(i_mp_qd) + eps )
           re(k,i,j,i_hs) = re_tmp(i_mp_qs)
           re(k,i,j,i_hg) = re_tmp(i_mp_qg)
           re(k,i,j,i_hh) = re_tmp(i_mp_qh)
 
        enddo
        enddo
        enddo
 
     end if
 
     return

Here is the caller graph for this function:

◆ atmos_phy_mp_suzuki10_qtrc2qhyd()

subroutine, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_qtrc2qhyd	(	integer, intent(in)	KA,
		integer, intent(in)	KS,
		integer, intent(in)	KE,
		integer, intent(in)	IA,
		integer, intent(in)	IS,
		integer, intent(in)	IE,
		integer, intent(in)	JA,
		integer, intent(in)	JS,
		integer, intent(in)	JE,
		real(rp), dimension(ka,ia,ja,nbin*nspc), intent(in)	QTRC0,
		real(rp), dimension (ka,ia,ja,n_hyd), intent(out)	Qe
	)

Calculate mass ratio of each category.

Definition at line 1402 of file scale_atmos_phy_mp_suzuki10.F90.

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

Referenced by mod_atmos_phy_mp_vars::atmos_phy_mp_vars_get_diagnostic().

     use scale_atmos_hydrometeor, only: &
        n_hyd, &
        i_hc,  &
        i_hr,  &
        i_hi,  &
        i_hs,  &
        i_hg,  &
        i_hh
     implicit none
 
     integer,  intent(in)  :: ka, ks, ke
     integer,  intent(in)  :: ia, is, ie
     integer,  intent(in)  :: ja, js, je
     real(RP), intent(in)  :: qtrc0(ka,ia,ja,nbin*nspc) ! tracer mass concentration [kg/kg]
     real(RP), intent(out) :: qe   (ka,ia,ja,n_hyd)     ! mixing ratio of each cateory [kg/kg]
 
     integer :: ihydro, ibin, iq, icateg
     integer :: k, i, j
     !---------------------------------------------------------------------------
 
 !OCL XFILL
     qe(:,:,:,:) = 0.0_rp
 
     do ihydro = 1, nspc
     do ibin   = 1, nbin
        iq = nbin*(ihydro-1) + ibin
 
        if    ( iq > 0                 .AND. iq <= nbin*(i_mp_qc -1) ) then ! liquid
           if    ( iq > 0    .AND. iq <= nbnd ) then ! cloud
              icateg = i_hc
           elseif( iq > nbnd .AND. iq <= nbin ) then ! rain
              icateg = i_hr
           endif
        elseif( iq > nbin*(i_mp_qc -1) .AND. iq <= nbin*(i_mp_qcl-1) ) then ! ice (column)
           icateg = i_hi
        elseif( iq > nbin*(i_mp_qcl-1) .AND. iq <= nbin*(i_mp_qp -1) ) then ! ice (plate)
           icateg = i_hi
        elseif( iq > nbin*(i_mp_qp -1) .AND. iq <= nbin*(i_mp_qs -1) ) then ! ice (dendrite)
           icateg = i_hi
        elseif( iq > nbin*(i_mp_qs -1) .AND. iq <= nbin*(i_mp_qg -1) ) then ! snow
           icateg = i_hs
        elseif( iq > nbin*(i_mp_qg -1) .AND. iq <= nbin*(i_mp_qh -1) ) then ! graupel
           icateg = i_hg
        elseif( iq > nbin*(i_mp_qh -1) .AND. iq <= nbin*nspc         ) then ! hail
           icateg = i_hh
        endif
 
        do j = js, je
        do i = is, ie
        do k = ks, ke
           qe(k,i,j,icateg) = qe(k,i,j,icateg) + qtrc0(k,i,j,iq)
        enddo
        enddo
        enddo
     enddo
     enddo
 
     return

Here is the caller graph for this function:

◆ atmos_phy_mp_suzuki10_qtrc2nhyd()

subroutine, public scale_atmos_phy_mp_suzuki10::atmos_phy_mp_suzuki10_qtrc2nhyd	(	integer, intent(in)	KA,
		integer, intent(in)	KS,
		integer, intent(in)	KE,
		integer, intent(in)	IA,
		integer, intent(in)	IS,
		integer, intent(in)	IE,
		integer, intent(in)	JA,
		integer, intent(in)	JS,
		integer, intent(in)	JE,
		real(rp), dimension (ka,ia,ja), intent(in)	DENS,
		real(rp), dimension(ka,ia,ja,nbin*nspc), intent(in)	QTRC0,
		real(rp), dimension (ka,ia,ja,n_hyd), intent(out)	Ne
	)

Calculate number concentration of each category.

Definition at line 1470 of file scale_atmos_phy_mp_suzuki10.F90.

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

Referenced by mod_atmos_phy_mp_vars::atmos_phy_mp_vars_get_diagnostic().

     use scale_atmos_hydrometeor, only: &
        n_hyd, &
        i_hc,  &
        i_hr,  &
        i_hi,  &
        i_hs,  &
        i_hg,  &
        i_hh
     implicit none
 
     integer,  intent(in)  :: ka, ks, ke
     integer,  intent(in)  :: ia, is, ie
     integer,  intent(in)  :: ja, js, je
     real(RP), intent(in)  :: dens (ka,ia,ja)           ! density [kg/m3]
     real(RP), intent(in)  :: qtrc0(ka,ia,ja,nbin*nspc) ! tracer mass concentration [kg/kg]
     real(RP), intent(out) :: ne   (ka,ia,ja,n_hyd)     ! number concentration of each cateory [1/m3]
 
     integer :: ihydro, ibin, iq, icateg
     integer :: k, i, j
     !---------------------------------------------------------------------------
 
 !OCL XFILL
     ne(:,:,:,:) = 0.0_rp
 
     do ihydro = 1, nspc
     do ibin   = 1, nbin
        iq = nbin*(ihydro-1) + ibin
 
        if    ( iq > 0                 .AND. iq <= nbin*(i_mp_qc -1) ) then ! liquid
           if    ( iq > 0    .AND. iq <= nbnd ) then ! cloud
              icateg = i_hc
           elseif( iq > nbnd .AND. iq <= nbin ) then ! rain
              icateg = i_hr
           endif
        elseif( iq > nbin*(i_mp_qc -1) .AND. iq <= nbin*(i_mp_qcl-1) ) then ! ice (column)
           icateg = i_hi
        elseif( iq > nbin*(i_mp_qcl-1) .AND. iq <= nbin*(i_mp_qp -1) ) then ! ice (plate)
           icateg = i_hi
        elseif( iq > nbin*(i_mp_qp -1) .AND. iq <= nbin*(i_mp_qs -1) ) then ! ice (dendrite)
           icateg = i_hi
        elseif( iq > nbin*(i_mp_qs -1) .AND. iq <= nbin*(i_mp_qg -1) ) then ! snow
           icateg = i_hs
        elseif( iq > nbin*(i_mp_qg -1) .AND. iq <= nbin*(i_mp_qh -1) ) then ! graupel
           icateg = i_hg
        elseif( iq > nbin*(i_mp_qh -1) .AND. iq <= nbin*nspc         ) then ! hail
           icateg = i_hh
        endif
 
        do j = js, je
        do i = is, ie
        do k = ks, ke
           ne(k,i,j,icateg) = ne(k,i,j,icateg) + dens(k,i,j) * qtrc0(k,i,j,iq) * rexpxctr(ibin)
        enddo
        enddo
        enddo
     enddo
     enddo
 
     return

Here is the caller graph for this function:

◆ atmos_phy_mp_suzuki10_qhyd2qtrc()

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

get mass ratio of each category

Definition at line 1538 of file scale_atmos_phy_mp_suzuki10.F90.

Referenced by mod_atmos_phy_mp_driver::atmos_phy_mp_driver_qhyd2qtrc().

     use scale_const, only: &
        pi => const_pi, &
        eps => const_eps
     use scale_atmos_hydrometeor, only: &
        n_hyd, &
        i_hc, &
        i_hr, &
        i_hi, &
        i_hs, &
        i_hg, &
        i_hh
     use scale_specfunc, only: &
        sf_gamma
     implicit none
     integer, intent(in) :: ka, ks, ke
     integer, intent(in) :: ia, is, ie
     integer, intent(in) :: ja, js, je
 
     real(RP), intent(in) :: qe(ka,ia,ja,n_hyd) ! mass ratio of each cateory [kg/kg]
 
     real(RP), intent(out)  :: qtrc(ka,ia,ja,qa-1)
 
     real(RP), intent(in), optional :: qnum(ka,ia,ja,n_hyd) ! number concentratio
 
     real(RP) :: coef0, coef1, coef2
     real(RP) :: dummy(nbin)
     real(RP) :: tmp_hyd, num_hyd, lambda_hyd
 
     integer :: k, i, j, iq
 
     if ( present(qnum) ) then
        log_warn("ATMOS_PHY_MP_suzuki10_qhyd2qtrc",*) 'At this moment, number concentratio is ignored'
     end if
 
     !--- define coefficients
     coef0 = 4.0_rp/3.0_rp*pi
     coef1 = 4.0_rp/3.0_rp*sqrt(pi/2.0_rp)
 
     if( nspc == 1 ) then !--- put all hydrometeors to liquid (warm bin)
 
        do j = js, je
        do i = is, ie
        do k = ks, ke
 
           tmp_hyd = 0.0_rp
           do iq = 1, nbin
              dummy(iq) = coef1 / sigma_sdf(1) * rho_sdf(1) * radc( iq )**3 &
                   * exp(  &
                   - ( log( radc(iq) )-log( r0_sdf(1) ) )**2*0.5_rp &
                   / sigma_sdf(1) / sigma_sdf(1) &
                   )
              tmp_hyd = tmp_hyd + dummy(iq)
           enddo
 
           coef2 = ( qe(k,i,j,i_hc) + qe(k,i,j,i_hr) &
                   + qe(k,i,j,i_hi) + qe(k,i,j,i_hs) + qe(k,i,j,i_hg) ) &
                 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
 
           do iq = 1, nbin
              qtrc(k+2,i,j,(il-1)*nbin+iq) = coef2 * dummy(iq)
           enddo
 
        enddo
        enddo
        enddo
 
     elseif( nspc > 1 ) then  !--- put each hydrometeor to each category (ice bin)
 
        do j = js, je
        do i = is, ie
        do k = ks, ke
 
           !--- Rain and Cloud put into liquid bin (log-normal)
           tmp_hyd = 0.0_rp
           do iq = 1, nbin
              dummy(iq) = coef1 / sigma_sdf(1) * rho_sdf(1) * radc( iq )**3 &
                   * exp(  &
                   - ( log( radc(iq) )-log( r0_sdf(1) ) )**2*0.5_rp &
                   / sigma_sdf(1) / sigma_sdf(1) &
                   )
              tmp_hyd = tmp_hyd + dummy(iq)
           enddo
 
           coef2 = ( qe(k,i,j,i_hc) + qe(k,i,j,i_hr) ) &
                 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
 
           do iq = 1, nbin
              qtrc(k,i,j,(il-1)*nbin+iq) = coef2 * dummy(iq)
           enddo
 
           !--- Ice put into plate bin (log-normal)
           tmp_hyd = 0.0_rp
           do iq = 1, nbin
              dummy(iq) = coef1 / sigma_sdf(2) * rho_sdf(2) * radc( iq )**3 &
                   * exp(  &
                   - ( log( radc(iq) )-log( r0_sdf(2) ) )**2*0.5_rp &
                   / sigma_sdf(2) / sigma_sdf(2) &
                   )
              tmp_hyd = tmp_hyd + dummy(iq)
           enddo
 
           coef2 = ( qe(k,i,j,i_hi) ) &
                 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
 
           do iq = 1, nbin
              qtrc(k,i,j,(ip-1)*nbin+iq) = coef2 * dummy(iq)
           enddo
 
           !--- Snow put into snow bin (gamma)
           num_hyd = coef0 * n0_sdf(3) * rho_sdf(3)
           lambda_hyd = ( pi * rho_sdf(3) / 6.0_rp *n0_sdf(3) * sf_gamma(4.0_rp) &
                / ( qe(k,i,j,i_hs) &
                + (0.50_rp-sign(0.50_rp,qe(k,i,j,i_hs)-eps)) &
                ) )**(0.25_rp)
 
           tmp_hyd = 0.0_rp
           do iq = 1, nbin
              dummy(iq) = num_hyd * radc( iq )**3 &
                   * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
              tmp_hyd = tmp_hyd + dummy(iq)
           enddo
 
           coef2 = ( qe(k,i,j,i_hs) ) &
                 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
 
           do iq = 1, nbin
              qtrc(k,i,j,(iss-1)*nbin+iq) = coef2 * dummy(iq)
           enddo
 
           !--- Graupel put into Graupel bin (gamma)
           num_hyd = coef0 * n0_sdf(4) * rho_sdf(4)
           lambda_hyd = ( pi * rho_sdf(4) / 6.0_rp *n0_sdf(4) * sf_gamma(4.0_rp) &
                / ( qe(k,i,j,i_hg) &
                + (0.50_rp-sign(0.50_rp,qe(k,i,j,i_hg)-eps)) &
                ) )**(0.25_rp)
 
           tmp_hyd = 0.0_rp
           do iq = 1, nbin
              dummy(iq) = num_hyd * radc( iq )**3 &
                   * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
              tmp_hyd = tmp_hyd + dummy(iq)
           enddo
 
           coef2 = ( qe(k,i,j,i_hg) ) &
                 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
 
           do iq = 1, nbin
              qtrc(k,i,j,(ig-1)*nbin+iq) = coef2 * dummy(iq)
           enddo
 
           !--- Hail put into Hail bin (gamma)
           num_hyd = coef0 * n0_sdf(5) * rho_sdf(5)
           lambda_hyd = ( pi * rho_sdf(5) / 6.0_rp *n0_sdf(5) * sf_gamma(4.0_rp) &
                / ( qe(k,i,j,i_hh) &
                + (0.50_rp-sign(0.50_rp,qe(k,i,j,i_hh)-eps)) &
                ) )**(0.25_rp)
 
           tmp_hyd = 0.0_rp
           do iq = 1, nbin
              dummy(iq) = num_hyd * radc( iq )**3 &
                   * exp( -lambda_hyd * 0.5_rp * radc( iq ) )
              tmp_hyd = tmp_hyd + dummy(iq)
           enddo
 
           coef2 = ( qe(k,i,j,i_hh) ) &
                 / ( tmp_hyd + ( 0.50_rp - sign(0.50_rp,tmp_hyd-eps) ) )
 
           do iq = 1, nbin
              qtrc(k,i,j,(ih-1)*nbin+iq) = coef2 * dummy(iq)
           enddo
 
        enddo
        enddo
        enddo
 
     endif
 
     do iq = nbin*nspc+1, qa-1
     do j = js, je
     do i = is, ie
     do k = ks, ke
        qtrc(k,i,j,iq) = 0.0_rp
     end do
     end do
     end do
     end do
 
     return