module ATMOSPHERE / Physics Cloud Microphysics More...

Functions/Subroutines
subroutine, public	atmos_phy_mp_kessler_config (MP_TYPE, QA, QS)
	Setup. More...

subroutine, public	atmos_phy_mp_kessler_setup
	Setup. More...

subroutine, public	atmos_phy_mp_kessler (DENS, MOMZ, MOMX, MOMY, RHOT, QTRC, CCN, EVAPORATE, SFLX_rain, SFLX_snow)
	Cloud Microphysics. More...

subroutine, public	atmos_phy_mp_kessler_cloudfraction (cldfrac, QTRC, mask_criterion)
	Calculate Cloud Fraction. More...

subroutine, public	atmos_phy_mp_kessler_effectiveradius (Re, QTRC0, DENS0, TEMP0)
	Calculate Effective Radius. More...

subroutine, public	atmos_phy_mp_kessler_mixingratio (Qe, QTRC0)
	Calculate mixing ratio of each category. More...

Variables
integer, parameter, public	qa_mp = 3

character(len=h_short), dimension(qa_mp), target, public	atmos_phy_mp_kessler_name

character(len=h_mid), dimension(qa_mp), target, public	atmos_phy_mp_kessler_desc

character(len=h_short), dimension(qa_mp), target, public	atmos_phy_mp_kessler_unit

real(rp), dimension(n_hyd), target, public	atmos_phy_mp_kessler_dens

Detailed Description

module ATMOSPHERE / Physics Cloud Microphysics

Description: Cloud Microphysics by Kessler-type parametarization Reference: Kessler(1969) Klemp and Wilhelmson(1978)

Author: Team SCALE

History

2012-01-14 (Y.Miyamoto) [new]
2012-03-23 (H.Yashiro) [mod] Explicit index parameter inclusion
2012-12-23 (H.Yashiro) [mod] Reconstruction
2015-09-08 (Y.Sato) [add] Add evaporated cloud number concentration

NAMELIST

PARAM_ATMOS_PHY_MP

name	type	default value	comment
MP_DOPRECIPITATION	logical	.true.	apply sedimentation (precipitation)?
MP_DONEGATIVE_FIXER	logical	.true.	apply negative fixer?
MP_LIMIT_NEGATIVE	real(RP)	1.0_RP	Abort if abs(fixed negative vaue) > abs(MP_limit_negative)
MP_NTMAX_SEDIMENTATION	integer	1	number of time step for sedimentation
MP_COUPLE_AEROSOL	logical	.false.	Consider CCN effect ?

History Output

name	description	unit	variable
Pcsat	QC production term by satadjust	kg/kg/s	QC_t_sat
Vterm_QR	terminal velocity of QR	m/s	vterm
pflux_QR	precipitation flux of QR	kg/m2/s	FLX_hydro

Function/Subroutine Documentation

◆ atmos_phy_mp_kessler_config()

subroutine, public scale_atmos_phy_mp_kessler::atmos_phy_mp_kessler_config	(	character(len=*), intent(in)	MP_TYPE,
		integer, intent(out)	QA,
		integer, intent(out)	QS
	)

Setup.

Definition at line 117 of file scale_atmos_phy_mp_kessler.F90.

References scale_atmos_hydrometeor::atmos_hydrometeor_regist(), atmos_phy_mp_kessler_desc, atmos_phy_mp_kessler_name, atmos_phy_mp_kessler_unit, scale_stdio::io_fid_log, scale_stdio::io_l, scale_process::prc_mpistop(), and qa_mp.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config().

     use scale_process, only: &
        prc_mpistop
     use scale_atmos_hydrometeor, only: &
        atmos_hydrometeor_regist
     implicit none
 
     character(len=*), intent(in)  :: MP_TYPE
     integer,          intent(out) :: QA
     integer,          intent(out) :: QS
     !---------------------------------------------------------------------------
 
     if( io_l ) write(io_fid_log,*)
     if( io_l ) write(io_fid_log,*) '++++++ Module[Cloud Microphysics Tracer] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
     if( io_l ) write(io_fid_log,*) '*** Tracers for KESSLER-type 1-moment bulk 3 category'
 
     if ( mp_type /= 'KESSLER' ) then
        write(*,*) 'xxx ATMOS_PHY_MP_TYPE is not KESSLER. Check!'
        call prc_mpistop
     endif
 
     call atmos_hydrometeor_regist( qs_mp,                     & ! [OUT]
                                    1, 2, 0,                   & ! [IN]
                                    atmos_phy_mp_kessler_name, & ! [IN]
                                    atmos_phy_mp_kessler_desc, & ! [IN]
                                    atmos_phy_mp_kessler_unit  ) ! [IN]
 
     qa = qa_mp
     qs = qs_mp
     qe_mp = qs_mp + qa_mp - 1
 
     i_qv = qs
     i_qc = qs + i_mp_qc
     i_qr = qs + i_mp_qr
 
     return

Here is the call graph for this function:

Here is the caller graph for this function:

◆ atmos_phy_mp_kessler_setup()

subroutine, public scale_atmos_phy_mp_kessler::atmos_phy_mp_kessler_setup ( )

Setup.

Definition at line 157 of file scale_atmos_phy_mp_kessler.F90.

References atmos_phy_mp_kessler_dens, scale_const::const_dwatr, scale_const::const_undef, scale_grid::grid_cdz, scale_atmos_hydrometeor::i_hc, scale_atmos_hydrometeor::i_hr, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_fid_nml, scale_stdio::io_l, scale_stdio::io_nml, 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_config().

     use scale_process, only: &
        prc_mpistop
     use scale_const, only: &
        const_undef, &
        const_dwatr
     use scale_time, only: &
        time_dtsec_atmos_phy_mp
     use scale_grid, only: &
        cdz => grid_cdz
     implicit none
 
     real(RP), parameter :: max_term_vel = 10.0_rp  !-- terminal velocity for calculate dt of sedimentation
     integer :: nstep_max
 
     namelist / param_atmos_phy_mp / &
        mp_doprecipitation,     &
        mp_donegative_fixer,    &
        mp_limit_negative,      &
        mp_ntmax_sedimentation, &
        mp_couple_aerosol
 
     integer :: ierr
 
     !---------------------------------------------------------------------------
 
     if( io_l ) write(io_fid_log,*)
     if( io_l ) write(io_fid_log,*) '++++++ Module[Cloud Microphysics] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
     if( io_l ) write(io_fid_log,*) '*** KESSLER-type 1-moment bulk 3 category'
 
     !--- read namelist
     rewind(io_fid_conf)
     read(io_fid_conf,nml=param_atmos_phy_mp,iostat=ierr)
     if( ierr < 0 ) then !--- missing
        if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
     elseif( ierr > 0 ) then !--- fatal error
        write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_MP. Check!'
        call prc_mpistop
     endif
     if( io_nml ) write(io_fid_nml,nml=param_atmos_phy_mp)
 
     if( mp_couple_aerosol ) then
        write(*,*) 'xxx MP_aerosol_couple should be .false. for KESSLER type MP!'
        call prc_mpistop
     endif
 
     if( io_l ) write(io_fid_log,*)
     if( io_l ) write(io_fid_log,*) '*** Enable negative fixer?                    : ', mp_donegative_fixer
     if( io_l ) write(io_fid_log,*) '*** Value limit of negative fixer (abs)       : ', abs(mp_limit_negative)
     if( io_l ) write(io_fid_log,*) '*** Enable sedimentation (precipitation)?     : ', mp_doprecipitation
 
     nstep_max = int( ( time_dtsec_atmos_phy_mp * max_term_vel ) / minval( cdz ) )
     mp_ntmax_sedimentation = max( mp_ntmax_sedimentation, nstep_max )
 
     mp_nstep_sedimentation  = mp_ntmax_sedimentation
     mp_rnstep_sedimentation = 1.0_rp / real(MP_ntmax_sedimentation,kind=rp)
     mp_dtsec_sedimentation  = time_dtsec_atmos_phy_mp * mp_rnstep_sedimentation
 
     if( io_l ) write(io_fid_log,*) '*** Timestep of sedimentation is divided into : ', mp_ntmax_sedimentation, 'step'
     if( io_l ) write(io_fid_log,*) '*** DT of sedimentation                       : ', mp_dtsec_sedimentation, '[s]'
 
     atmos_phy_mp_kessler_dens(:)    = const_undef
     atmos_phy_mp_kessler_dens(i_hc) = const_dwatr
     atmos_phy_mp_kessler_dens(i_hr) = const_dwatr
 
     allocate( factor_vterm(ka) )
 
     return

Here is the call graph for this function:

Here is the caller graph for this function:

◆ atmos_phy_mp_kessler()

subroutine, public scale_atmos_phy_mp_kessler::atmos_phy_mp_kessler	(	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,qa), 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 239 of file scale_atmos_phy_mp_kessler.F90.

References scale_atmos_phy_mp_common::atmos_phy_mp_negative_fixer(), scale_atmos_phy_mp_common::atmos_phy_mp_precipitation(), scale_atmos_phy_mp_common::atmos_phy_mp_saturation_adjustment(), scale_atmos_refstate::atmos_refstate_dens, scale_comm::comm_horizontal_mean(), scale_const::const_dwatr, scale_const::const_pi, scale_atmos_hydrometeor::i_qc, scale_atmos_hydrometeor::i_qr, scale_atmos_hydrometeor::i_qv, 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, scale_atmos_hydrometeor::lhv, scale_prof::prof_rapend(), scale_prof::prof_rapstart(), scale_tracer::qa, qa_mp, scale_time::time_dtsec_atmos_phy_mp, scale_tracer::tracer_cv, scale_tracer::tracer_mass, and scale_tracer::tracer_r.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config().

     use scale_grid_index
     use scale_const, only: &
        dwatr => const_dwatr, &
        pi    => const_pi
     use scale_comm, only: &
        comm_horizontal_mean
     use scale_time, only: &
        dt_mp => time_dtsec_atmos_phy_mp
     use scale_history, only: &
        hist_in
     use scale_tracer, only: &
        qa,         &
        tracer_r,   &
        tracer_cv,  &
        tracer_mass
     use scale_atmos_thermodyn, only: &
        thermodyn_rhoe        => atmos_thermodyn_rhoe,       &
        thermodyn_rhot        => atmos_thermodyn_rhot,       &
        thermodyn_temp_pres_e => atmos_thermodyn_temp_pres_e
     use scale_atmos_phy_mp_common, only: &
        mp_negative_fixer        => atmos_phy_mp_negative_fixer,       &
        mp_precipitation         => atmos_phy_mp_precipitation,        &
        mp_saturation_adjustment => atmos_phy_mp_saturation_adjustment
     implicit none
 
     real(RP), intent(inout) :: DENS     (KA,IA,JA)
     real(RP), intent(inout) :: MOMZ     (KA,IA,JA)
     real(RP), intent(inout) :: MOMX     (KA,IA,JA)
     real(RP), intent(inout) :: MOMY     (KA,IA,JA)
     real(RP), intent(inout) :: RHOT     (KA,IA,JA)
     real(RP), intent(inout) :: QTRC     (KA,IA,JA,QA)
     real(RP), intent(in)    :: CCN      (KA,IA,JA)
     real(RP), intent(out)   :: EVAPORATE(KA,IA,JA)   !--- number of evaporated cloud [/m3]
     real(RP), intent(out)   :: SFLX_rain(IA,JA)
     real(RP), intent(out)   :: SFLX_snow(IA,JA)
 
     real(RP) :: RHOE_t   (KA,IA,JA)
     real(RP) :: QTRC_t   (KA,IA,JA,QA)
     real(RP) :: RHOE     (KA,IA,JA)
     real(RP) :: temp     (KA,IA,JA)
     real(RP) :: pres     (KA,IA,JA)
 
     real(RP) :: vterm    (KA,IA,JA,QA_MP-1) ! terminal velocity of each tracer  [m/s]
     real(RP) :: pflux    (KA,IA,JA,QA_MP-1) ! precipitation flux of each tracer [kg/m2/s]
     real(RP) :: FLX_hydro(KA,IA,JA,QA_MP-1)
     real(RP) :: QC_t_sat (KA,IA,JA)
     real(RP) :: rho_prof (KA)               ! averaged profile of rho
 
     integer  :: step
     integer  :: k, i, j
     !---------------------------------------------------------------------------
 
     if( io_l ) write(io_fid_log,*) '*** Atmos physics  step: Cloud microphysics(kessler)'
 
     if ( first ) then
        ! Calculate collection factor for terminal velocity of QR
        call comm_horizontal_mean( rho_prof(:), dens(:,:,:) )
        rho_prof(:) = rho_prof(:) * 1.e-3_rp ! [kg/m3]->[g/cc]
 
        do k = ks, ke
           factor_vterm(k) = sqrt( rho_prof(ks)/rho_prof(k) )
        enddo
        first = .false.
     endif
 
     if ( mp_donegative_fixer ) then
        call mp_negative_fixer( dens(:,:,:),      & ! [INOUT]
                                rhot(:,:,:),      & ! [INOUT]
                                qtrc(:,:,:,:),    & ! [INOUT]
                                i_qv,             & ! [IN]
                                mp_limit_negative ) ! [IN]
     endif
 
     call thermodyn_rhoe( rhoe(:,:,:),   & ! [OUT]
                          rhot(:,:,:),   & ! [IN]
                          qtrc(:,:,:,:), & ! [IN]
                          tracer_cv(:),       & ! [IN]
                          tracer_r(:),       & ! [IN]
                          tracer_mass(:)        ) ! [IN]
 
     !##### MP Main #####
     call mp_kessler( rhoe_t(:,:,:),   & ! [OUT]
                      qtrc_t(:,:,:,:), & ! [OUT]
                      rhoe(:,:,:),   & ! [INOUT]
                      qtrc(:,:,:,:), & ! [INOUT]
                      dens(:,:,:)    ) ! [IN]
 
     vterm(:,:,:,:) = 0.0_rp
     flx_hydro(:,:,:,:) = 0.0_rp
 
     if ( mp_doprecipitation ) then
 
        do step = 1, mp_nstep_sedimentation
 
           call thermodyn_temp_pres_e( temp(:,:,:),   & ! [OUT]
                                       pres(:,:,:),   & ! [OUT]
                                       dens(:,:,:),   & ! [IN]
                                       rhoe(:,:,:),   & ! [IN]
                                       qtrc(:,:,:,:), & ! [IN]
                                       tracer_cv(:),       & ! [IN]
                                       tracer_r(:),       & ! [IN]
                                       tracer_mass(:)        ) ! [IN]
 
           call mp_kessler_vterm( vterm(:,:,:,:), & ! [OUT]
                                  dens(:,:,:),   & ! [IN]
                                  qtrc(:,:,:,:)  ) ! [IN]
 
           call mp_precipitation( qa_mp,                 & ! [IN]
                                  qs_mp,                 & ! [IN]
                                  pflux(:,:,:,:),    & ! [OUT]
                                  vterm(:,:,:,:),    & ! [INOUT]
                                  dens(:,:,:),      & ! [INOUT]
                                  momz(:,:,:),      & ! [INOUT]
                                  momx(:,:,:),      & ! [INOUT]
                                  momy(:,:,:),      & ! [INOUT]
                                  rhoe(:,:,:),      & ! [INOUT]
                                  qtrc(:,:,:,:),    & ! [INOUT]
                                  temp(:,:,:),      & ! [IN]
                                  tracer_cv(:),          & ! [IN]
                                  mp_dtsec_sedimentation ) ! [IN]
 
           do j  = js, je
           do i  = is, ie
           do k  = ks-1, ke-1
              flx_hydro(k,i,j,i_mp_qr) = flx_hydro(k,i,j,i_mp_qr) + pflux(k,i,j,i_mp_qr) * mp_rnstep_sedimentation
           enddo
           enddo
           enddo
 
        enddo
 
     endif
 
     call hist_in( flx_hydro(:,:,:,i_mp_qr), 'pflux_QR', 'precipitation flux of QR', 'kg/m2/s', nohalo=.true. )
 
     call hist_in( vterm(:,:,:,i_mp_qr), 'Vterm_QR', 'terminal velocity of QR', 'm/s' )
 
     do j = js, je
     do i = is, ie
     do k = ks, ke
        qc_t_sat(k,i,j) = qtrc(k,i,j,i_qc)
     enddo
     enddo
     enddo
 
     call mp_saturation_adjustment( rhoe_t(:,:,:),       & ! [INOUT]
                                    qtrc_t(:,:,:,:),     & ! [INOUT]
                                    rhoe(:,:,:),       & ! [INOUT]
                                    qtrc(:,:,:,:),     & ! [INOUT]
                                    dens(:,:,:),       & ! [IN]
                                    i_qv,                & ! [IN]
                                    i_qc,                & ! [IN]
                                    -1,                  & ! [IN]
                                    flag_liquid = .true. ) ! [IN]
 
     do j = js, je
     do i = is, ie
     do k = ks, ke
        qc_t_sat(k,i,j) = ( qtrc(k,i,j,i_qc) - qc_t_sat(k,i,j) ) / dt_mp
     enddo
     enddo
     enddo
 
     call hist_in( qc_t_sat(:,:,:), 'Pcsat', 'QC production term by satadjust', 'kg/kg/s' )
 
     do j = js, je
     do i = is, ie
     do k = ks, ke
        evaporate(k,i,j) = max( -qc_t_sat(k,i,j), 0.0_rp ) ! if negative, condensation
        evaporate(k,i,j) = evaporate(k,i,j) * dens(k,i,j) / (4.0_rp/3.0_rp*pi*dwatr*re_qc**3) ! mass -> number (assuming constant particle radius as re_qc)
     enddo
     enddo
     enddo
 
     !##### END MP Main #####
 
     call thermodyn_rhot( rhot(:,:,:),   & ! [OUT]
                          rhoe(:,:,:),   & ! [IN]
                          qtrc(:,:,:,:), & ! [IN]
                          tracer_cv(:),  & ! [IN]
                          tracer_r(:),   & ! [IN]
                          tracer_mass(:) ) ! [IN]
 
     if ( mp_donegative_fixer ) then
        call mp_negative_fixer( dens(:,:,:),      & ! [INOUT]
                                rhot(:,:,:),      & ! [INOUT]
                                qtrc(:,:,:,:),    & ! [INOUT]
                                i_qv,             & ! [IN]
                                mp_limit_negative ) ! [IN]
     endif
 
     do j = js, je
     do i = is, ie
        sflx_rain(i,j) = - flx_hydro(ks-1,i,j,i_mp_qr)
        sflx_snow(i,j) = 0.0_rp
     enddo
     enddo
 
     return

Here is the call graph for this function:

Here is the caller graph for this function:

◆ atmos_phy_mp_kessler_cloudfraction()

subroutine, public scale_atmos_phy_mp_kessler::atmos_phy_mp_kessler_cloudfraction	(	real(rp), dimension(ka,ia,ja), intent(out)	cldfrac,
		real(rp), dimension (ka,ia,ja,qa), intent(in)	QTRC,
		real(rp), intent(in)	mask_criterion
	)

Calculate Cloud Fraction.

Definition at line 639 of file scale_atmos_phy_mp_kessler.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_config().

     use scale_grid_index
     use scale_tracer, only: &
        qa
     implicit none
 
     real(RP), intent(out) :: cldfrac(KA,IA,JA)
     real(RP), intent(in)  :: QTRC   (KA,IA,JA,QA)
     real(RP), intent(in)  :: mask_criterion
 
     real(RP) :: qhydro
     integer  :: k, i, j, iq
     !---------------------------------------------------------------------------
 
     do j  = js, je
     do i  = is, ie
     do k  = ks, ke
        qhydro = 0.0_rp
        do iq = qs_mp+1, qe_mp
           qhydro = qhydro + qtrc(k,i,j,iq)
        enddo
        cldfrac(k,i,j) = 0.5_rp + sign(0.5_rp,qhydro-mask_criterion)
     enddo
     enddo
     enddo
 
     return

Here is the caller graph for this function:

◆ atmos_phy_mp_kessler_effectiveradius()

subroutine, public scale_atmos_phy_mp_kessler::atmos_phy_mp_kessler_effectiveradius	(	real(rp), dimension (ka,ia,ja,n_hyd), intent(out)	Re,
		real(rp), dimension(ka,ia,ja,qa), 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 674 of file scale_atmos_phy_mp_kessler.F90.

References scale_atmos_hydrometeor::n_hyd, and scale_tracer::qa.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config().

     use scale_grid_index
     use scale_tracer, only: &
        qa
     use scale_atmos_hydrometeor, only: &
        n_hyd
     implicit none
     real(RP), intent(out) :: Re   (KA,IA,JA,N_HYD) ! effective radius          [cm]
     real(RP), intent(in)  :: QTRC0(KA,IA,JA,QA)    ! tracer mass concentration [kg/kg]
     real(RP), intent(in)  :: DENS0(KA,IA,JA)       ! density                   [kg/m3]
     real(RP), intent(in)  :: TEMP0(KA,IA,JA)       ! temperature               [K]
 
     real(RP), parameter :: um2cm = 100.0_rp
 
     integer :: ih
     !---------------------------------------------------------------------------
 
     do ih = 1, n_hyd
        if ( ih == i_hc ) then
           re(:,:,:,ih) =   8.e-6_rp * um2cm
        else if ( ih == i_hr ) then
           re(:,:,:,ih) = 100.e-6_rp * um2cm
        else
           re(:,:,:,ih) = 0.0_rp
        end if
     end do
 
     return

Here is the caller graph for this function:

◆ atmos_phy_mp_kessler_mixingratio()

subroutine, public scale_atmos_phy_mp_kessler::atmos_phy_mp_kessler_mixingratio	(	real(rp), dimension (ka,ia,ja,n_hyd), intent(out)	Qe,
		real(rp), dimension(ka,ia,ja,qa), intent(in)	QTRC0
	)

Calculate mixing ratio of each category.

Definition at line 708 of file scale_atmos_phy_mp_kessler.F90.

References scale_atmos_hydrometeor::n_hyd, and scale_tracer::qa.

Referenced by scale_atmos_phy_mp::atmos_phy_mp_config().

     use scale_grid_index
     use scale_tracer, only: &
        qa
     use scale_atmos_hydrometeor, only: &
        n_hyd
     implicit none
     real(RP), intent(out) :: Qe   (KA,IA,JA,N_HYD) ! mixing ratio of each cateory [kg/kg]
     real(RP), intent(in)  :: QTRC0(KA,IA,JA,QA)    ! tracer mass concentration [kg/kg]
 
     integer :: ih
     !---------------------------------------------------------------------------
 
     do ih = 1, n_hyd
        if ( ih == i_hc ) then
           qe(:,:,:,ih) = qtrc0(:,:,:,i_qc)
        else if ( ih == i_hr ) then
           qe(:,:,:,ih) = qtrc0(:,:,:,i_qr)
        else
           qe(:,:,:,ih) = 0.0_rp
        end if
     end do
 
     return