mod_atmos_phy_rd_driver.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module mod_atmos_phy_rd_driver
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_atmos_grid_cartesc_index
  use scale_tracer
  use scale_cpl_sfc_index
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_rd_driver_setup
  public :: atmos_phy_rd_driver_calc_tendency

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_rd_driver_setup
    use scale_prc, only: &
       prc_abort
    use scale_atmos_phy_rd_mstrnx, only: &
       atmos_phy_rd_mstrnx_setup
    use scale_atmos_phy_rd_offline, only: &
       atmos_phy_rd_offline_setup
    use scale_atmos_grid_cartesc, only: &
       cz => atmos_grid_cartesc_cz, &
       fz => atmos_grid_cartesc_fz
    use mod_atmos_admin, only: &
       atmos_phy_rd_type, &
       atmos_sw_phy_rd
    use mod_atmos_phy_rd_vars, only: &
       sfcflx_lw_up => atmos_phy_rd_sflx_lw_up,   &
       sfcflx_lw_dn => atmos_phy_rd_sflx_lw_dn,   &
       sfcflx_sw_up => atmos_phy_rd_sflx_sw_up,   &
       sfcflx_sw_dn => atmos_phy_rd_sflx_sw_dn,   &
       toaflx_lw_up => atmos_phy_rd_toaflx_lw_up, &
       toaflx_lw_dn => atmos_phy_rd_toaflx_lw_dn, &
       toaflx_sw_up => atmos_phy_rd_toaflx_sw_up, &
       toaflx_sw_dn => atmos_phy_rd_toaflx_sw_dn, &
       sflx_rad_dn  => atmos_phy_rd_sflx_down,    &
       solins       => atmos_phy_rd_solins,       &
       cossza       => atmos_phy_rd_cossza
    implicit none
    !---------------------------------------------------------------------------

    log_newline
    log_info("ATMOS_PHY_RD_driver_setup",*) 'Setup'

    if ( atmos_sw_phy_rd ) then

       select case ( atmos_phy_rd_type )
       case ( "MSTRNX" )
          call atmos_phy_rd_mstrnx_setup( ka, ks, ke, cz(:), fz(:) )
       case ( "OFFLINE" )
          call atmos_phy_rd_offline_setup
       case default
          log_error("ATMOS_PHY_RD_driver_setup",*) 'invalid Radiation type(', trim(atmos_phy_rd_type), '). CHECK!'
          call prc_abort
       end select

    else

       log_info("ATMOS_PHY_RD_driver_setup",*) 'this component is never called.'
       log_info("ATMOS_PHY_RD_driver_setup",*) 'radiation fluxes are set to zero.'
       sfcflx_lw_up(:,:)     = 0.0_rp
       sfcflx_lw_dn(:,:)     = 0.0_rp
       sfcflx_sw_up(:,:)     = 0.0_rp
       sfcflx_sw_dn(:,:)     = 0.0_rp
       toaflx_lw_up(:,:)     = 0.0_rp
       toaflx_lw_dn(:,:)     = 0.0_rp
       toaflx_sw_up(:,:)     = 0.0_rp
       toaflx_sw_dn(:,:)     = 0.0_rp
       sflx_rad_dn(:,:,:,:) = 0.0_rp
       solins(:,:)     = 0.0_rp
       cossza(:,:)     = 0.0_rp

    endif

    return
  end subroutine atmos_phy_rd_driver_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_rd_driver_calc_tendency( update_flag )
    use scale_atmos_grid_cartesc_real, only: &
       real_cz  => atmos_grid_cartesc_real_cz,  &
       real_fz  => atmos_grid_cartesc_real_fz,  &
       real_lon => atmos_grid_cartesc_real_lon, &
       real_lat => atmos_grid_cartesc_real_lat, &
       atmos_grid_cartesc_real_vol, &
       atmos_grid_cartesc_real_totvol
    use scale_landuse, only: &
       fact_ocean => landuse_fact_ocean, &
       fact_land  => landuse_fact_land,  &
       fact_urban => landuse_fact_urban
    use scale_time, only: &
       time_nowdate,                     &
       time_nowdaysec,                   &
       time_offset_year
    use scale_statistics, only: &
       statistics_checktotal, &
       statistics_total
    use scale_file_history, only: &
       file_history_in
    use scale_atmos_hydrometeor, only: &
       n_hyd
    use scale_atmos_aerosol, only: &
       n_ae
    use mod_atmos_admin, only: &
       atmos_phy_rd_type
    use scale_atmos_solarins, only: &
       solarins_insolation => atmos_solarins_insolation
    use scale_atmos_phy_rd_mstrnx, only: &
       atmos_phy_rd_mstrnx_flux
    use scale_atmos_phy_rd_offline, only: &
       atmos_phy_rd_offline_flux
    use scale_atmos_phy_rd_common, only: &
       atmos_phy_rd_calc_heating, &
       i_sw, &
       i_lw, &
       i_dn, &
       i_up
    use mod_atmos_vars, only: &
       temp,              &
       pres,              &
       qv,                &
       cvtot,             &
       dens   => dens_av, &
       qtrc   => qtrc_av, &
       rhoh   => rhoh_p
    use mod_atmos_phy_sf_vars, only: &
       sfc_temp    => atmos_phy_sf_sfc_temp,   &
       sfc_albedo  => atmos_phy_sf_sfc_albedo
    use mod_atmos_phy_rd_vars, only: &
       rhoh_rd      => atmos_phy_rd_rhoh,         &
       sfcflx_lw_up => atmos_phy_rd_sflx_lw_up,   &
       sfcflx_lw_dn => atmos_phy_rd_sflx_lw_dn,   &
       sfcflx_sw_up => atmos_phy_rd_sflx_sw_up,   &
       sfcflx_sw_dn => atmos_phy_rd_sflx_sw_dn,   &
       toaflx_lw_up => atmos_phy_rd_toaflx_lw_up, &
       toaflx_lw_dn => atmos_phy_rd_toaflx_lw_dn, &
       toaflx_sw_up => atmos_phy_rd_toaflx_sw_up, &
       toaflx_sw_dn => atmos_phy_rd_toaflx_sw_dn, &
       sflx_rad_dn  => atmos_phy_rd_sflx_down,    &
       solins       => atmos_phy_rd_solins,       &
       cossza       => atmos_phy_rd_cossza
    use mod_atmos_vars, only: &
       atmos_vars_get_diagnostic
    use mod_atmos_phy_mp_vars, only: &
       atmos_phy_mp_vars_get_diagnostic
    use mod_atmos_phy_ae_vars, only: &
       atmos_phy_ae_vars_get_diagnostic
    implicit none

    logical, intent(in) :: update_flag

    real(RP) :: TEMP_t      (ka,ia,ja,3)
    real(RP) :: flux_rad    (ka,ia,ja,2,2,2)
    real(RP) :: flux_rad_top(   ia,ja,2,2,2)

    real(RP) :: dtau_s      (ka,ia,ja) ! 0.67 micron cloud optical depth
    real(RP) :: dem_s       (ka,ia,ja) ! 10.5 micron cloud emissivity

    real(RP) :: flux_up     (ka,ia,ja,2)
    real(RP) :: flux_dn     (ka,ia,ja,2)
    real(RP) :: flux_net    (ka,ia,ja,2)
    real(RP) :: flux_net_toa(   ia,ja,2)
    real(RP) :: flux_net_sfc(   ia,ja,2)

    real(RP) :: SFCFLX_LW_up_c(ia,ja)
    real(RP) :: SFCFLX_LW_dn_c(ia,ja)
    real(RP) :: SFCFLX_SW_up_c(ia,ja)
    real(RP) :: SFCFLX_SW_dn_c(ia,ja)
    real(RP) :: TOAFLX_LW_up_c(ia,ja)
    real(RP) :: TOAFLX_LW_dn_c(ia,ja)
    real(RP) :: TOAFLX_SW_up_c(ia,ja)
    real(RP) :: TOAFLX_SW_dn_c(ia,ja)

    real(RP) :: CLDFRAC(ka,ia,ja)
    real(RP) :: MP_Re  (ka,ia,ja,n_hyd)
    real(RP) :: MP_Qe  (ka,ia,ja,n_hyd)
    real(RP) :: AE_Re  (ka,ia,ja,n_ae)
    real(RP) :: AE_Qe  (ka,ia,ja,n_ae)

    real(RP) :: RH(ka,ia,ja)

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

    if ( update_flag ) then

       call solarins_insolation( ia, is, ie, ja, js, je, &
                                 real_lon(:,:), real_lat(:,:),      & ! [IN]
                                 time_nowdate(:), time_offset_year, & ! [IN]
                                 solins(:,:), cossza(:,:)         ) ! [OUT]

       call atmos_phy_mp_vars_get_diagnostic( &
            dens(:,:,:), temp(:,:,:), qtrc(:,:,:,:), & ! [IN]
            cldfrac=cldfrac, re=mp_re, qe=mp_qe      ) ! [IN]

       call atmos_vars_get_diagnostic( "RH", rh )
       call atmos_phy_ae_vars_get_diagnostic( &
            qtrc(:,:,:,:), rh(:,:,:), & ! [IN]
            re=ae_re                  ) ! [IN]
!            Re=AE_Re, Qe=AE_Qe      ) ! [IN]


       select case ( atmos_phy_rd_type )
       case ( "MSTRNX" )

          call atmos_phy_rd_mstrnx_flux( &
               ka, ks, ke, ia, is, ie, ja, js, je, &
               dens(:,:,:), temp(:,:,:), pres(:,:,:), qv(:,:,:), & ! [IN]
               real_cz(:,:,:), real_fz(:,:,:),                   & ! [IN]
               fact_ocean(:,:), fact_land(:,:), fact_urban(:,:), & ! [IN]
               sfc_temp(:,:), sfc_albedo(:,:,:,:),               & ! [IN]
               solins(:,:), cossza(:,:),                         & ! [IN]
               cldfrac(:,:,:), mp_re(:,:,:,:), mp_qe(:,:,:,:),   & ! [IN]
               ae_re(:,:,:,:),                                   & ! [IN]
!               AE_Re(:,:,:,:), AE_Qe(:,:,:,:),                   & ! [IN]
               flux_rad(:,:,:,:,:,:),                            & ! [OUT]
               flux_rad_top(:,:,:,:,:), sflx_rad_dn(:,:,:,:),    & ! [OUT]
               dtau_s = dtau_s(:,:,:), dem_s = dem_s(:,:,:)      ) ! [OUT]

       case ( "OFFLINE" )

          call atmos_phy_rd_offline_flux( &
               ka, ks, ke, ia, is, ie, ja, js, je, &
               time_nowdaysec,        & ! [IN]
               flux_rad(:,:,:,:,:,2), & ! [OUT]
               sflx_rad_dn(:,:,:,:)   ) ! [OUT]
          flux_rad(:,:,:,:,:,1)   = 0.0_rp ! clear sky
          flux_rad_top(:,:,:,:,:) = 0.0_rp
          dtau_s(:,:,:) = 0.0_rp
          dem_s(:,:,:) = 0.0_rp

       end select


!OCL XFILL
       do j = js, je
       do i = is, ie
          ! for clear-sky
          sfcflx_lw_up_c(i,j)    = flux_rad(ks-1,i,j,i_lw,i_up,1)
          sfcflx_lw_dn_c(i,j)    = flux_rad(ks-1,i,j,i_lw,i_dn,1)
          sfcflx_sw_up_c(i,j)    = flux_rad(ks-1,i,j,i_sw,i_up,1)
          sfcflx_sw_dn_c(i,j)    = flux_rad(ks-1,i,j,i_sw,i_dn,1)
          ! for all-sky
          sfcflx_lw_up(i,j)    = flux_rad(ks-1,i,j,i_lw,i_up,2)
          sfcflx_lw_dn(i,j)    = flux_rad(ks-1,i,j,i_lw,i_dn,2)
          sfcflx_sw_up(i,j)    = flux_rad(ks-1,i,j,i_sw,i_up,2)
          sfcflx_sw_dn(i,j)    = flux_rad(ks-1,i,j,i_sw,i_dn,2)

          flux_net_sfc(i,j,i_lw) = sfcflx_lw_up(i,j) - sfcflx_lw_dn(i,j)
          flux_net_sfc(i,j,i_sw) = sfcflx_sw_up(i,j) - sfcflx_sw_dn(i,j)
       enddo
       enddo

!OCL XFILL
       do j = js, je
       do i = is, ie
          ! for clear-sky
          toaflx_lw_up_c(i,j)    = flux_rad_top(i,j,i_lw,i_up,1)
          toaflx_lw_dn_c(i,j)    = flux_rad_top(i,j,i_lw,i_dn,1)
          toaflx_sw_up_c(i,j)    = flux_rad_top(i,j,i_sw,i_up,1)
          toaflx_sw_dn_c(i,j)    = flux_rad_top(i,j,i_sw,i_dn,1)
          ! for all-sky
          toaflx_lw_up(i,j)    = flux_rad_top(i,j,i_lw,i_up,2)
          toaflx_lw_dn(i,j)    = flux_rad_top(i,j,i_lw,i_dn,2)
          toaflx_sw_up(i,j)    = flux_rad_top(i,j,i_sw,i_up,2)
          toaflx_sw_dn(i,j)    = flux_rad_top(i,j,i_sw,i_dn,2)

          flux_net_toa(i,j,i_lw) = toaflx_lw_up(i,j) - toaflx_lw_dn(i,j)
          flux_net_toa(i,j,i_sw) = toaflx_sw_up(i,j) - toaflx_sw_dn(i,j)
       enddo
       enddo

!OCL XFILL
       do j = js, je
       do i = is, ie
       do k = ks, ke
          flux_up(k,i,j,i_lw) = 0.5_rp * ( flux_rad(k-1,i,j,i_lw,i_up,2) + flux_rad(k,i,j,i_lw,i_up,2) )
          flux_dn(k,i,j,i_lw) = 0.5_rp * ( flux_rad(k-1,i,j,i_lw,i_dn,2) + flux_rad(k,i,j,i_lw,i_dn,2) )
          flux_up(k,i,j,i_sw) = 0.5_rp * ( flux_rad(k-1,i,j,i_sw,i_up,2) + flux_rad(k,i,j,i_sw,i_up,2) )
          flux_dn(k,i,j,i_sw) = 0.5_rp * ( flux_rad(k-1,i,j,i_sw,i_dn,2) + flux_rad(k,i,j,i_sw,i_dn,2) )

          flux_net(k,i,j,i_lw) = flux_up(k,i,j,i_lw) - flux_dn(k,i,j,i_lw)
          flux_net(k,i,j,i_sw) = flux_up(k,i,j,i_sw) - flux_dn(k,i,j,i_sw)
       enddo
       enddo
       enddo

       ! apply radiative flux convergence -> heating rate
       call atmos_phy_rd_calc_heating( &
            ka, ks, ke, ia, is, ie, ja, js, je, &
            flux_rad(:,:,:,:,:,2),     & ! [IN]
            dens(:,:,:), temp(:,:,:),  & ! [IN]
            cvtot(:,:,:),              & ! [IN]
            real_fz(:,:,:),            & ! [IN]
            rhoh_rd(:,:,:),            & ! [OUT]
            temp_t = temp_t(:,:,:,:)   ) ! [OUT]



       call file_history_in( solins(:,:), 'SOLINS', 'solar insolation',        'W/m2', fill_halo=.true. )
       call file_history_in( cossza(:,:), 'COSZ',   'cos(solar zenith angle)', '1',    fill_halo=.true. )

       call file_history_in( sfcflx_lw_up_c(:,:),      'SFLX_LW_up_c',   'SFC upward   longwave  radiation flux (clr)', 'W/m2', fill_halo=.true. )
       call file_history_in( sfcflx_lw_dn_c(:,:),      'SFLX_LW_dn_c',   'SFC downward longwave  radiation flux (clr)', 'W/m2', fill_halo=.true. )
       call file_history_in( sfcflx_sw_up_c(:,:),      'SFLX_SW_up_c',   'SFC upward   shortwave radiation flux (clr)', 'W/m2', fill_halo=.true. )
       call file_history_in( sfcflx_sw_dn_c(:,:),      'SFLX_SW_dn_c',   'SFC downward shortwave radiation flux (clr)', 'W/m2', fill_halo=.true. )

       call file_history_in( sfcflx_lw_up(:,:),      'SFLX_LW_up',     'SFC upward   longwave  radiation flux',       'W/m2', fill_halo=.true. )
       call file_history_in( sfcflx_lw_dn(:,:),      'SFLX_LW_dn',     'SFC downward longwave  radiation flux',       'W/m2', fill_halo=.true. )
       call file_history_in( sfcflx_sw_up(:,:),      'SFLX_SW_up',     'SFC upward   shortwave radiation flux',       'W/m2', fill_halo=.true. )
       call file_history_in( sfcflx_sw_dn(:,:),      'SFLX_SW_dn',     'SFC downward shortwave radiation flux',       'W/m2', fill_halo=.true. )

       call file_history_in( flux_net_sfc(:,:,i_lw), 'SFLX_LW_net',    'SFC net      longwave  radiation flux',       'W/m2', fill_halo=.true. )
       call file_history_in( flux_net_sfc(:,:,i_sw), 'SFLX_SW_net',    'SFC net      shortwave radiation flux',       'W/m2', fill_halo=.true. )

       call file_history_in( toaflx_lw_up_c(:,:),      'TOAFLX_LW_up_c', 'TOA upward   longwave  radiation flux (clr)', 'W/m2', fill_halo=.true. )
       call file_history_in( toaflx_lw_dn_c(:,:),      'TOAFLX_LW_dn_c', 'TOA downward longwave  radiation flux (clr)', 'W/m2', fill_halo=.true. )
       call file_history_in( toaflx_sw_up_c(:,:),      'TOAFLX_SW_up_c', 'TOA upward   shortwave radiation flux (clr)', 'W/m2', fill_halo=.true. )
       call file_history_in( toaflx_sw_dn_c(:,:),      'TOAFLX_SW_dn_c', 'TOA downward shortwave radiation flux (clr)', 'W/m2', fill_halo=.true. )

       call file_history_in( toaflx_lw_up(:,:),      'TOAFLX_LW_up',   'TOA upward   longwave  radiation flux',       'W/m2', fill_halo=.true. )
       call file_history_in( toaflx_lw_dn(:,:),      'TOAFLX_LW_dn',   'TOA downward longwave  radiation flux',       'W/m2', fill_halo=.true. )
       call file_history_in( toaflx_sw_up(:,:),      'TOAFLX_SW_up',   'TOA upward   shortwave radiation flux',       'W/m2', fill_halo=.true. )
       call file_history_in( toaflx_sw_dn(:,:),      'TOAFLX_SW_dn',   'TOA downward shortwave radiation flux',       'W/m2', fill_halo=.true. )

       call file_history_in( flux_net_toa(:,:,i_lw), 'TOAFLX_LW_net',  'TOA net      longwave  radiation flux',       'W/m2', fill_halo=.true. )
       call file_history_in( flux_net_toa(:,:,i_sw), 'TOAFLX_SW_net',  'TOA net      shortwave radiation flux',       'W/m2', fill_halo=.true. )

       call file_history_in( flux_net_sfc(:,:,i_lw),   'SLR',          'SFC net longwave  radiation flux',  'W/m2', fill_halo=.true. )
       call file_history_in( flux_net_sfc(:,:,i_sw),   'SSR',          'SFC net shortwave radiation flux',  'W/m2', fill_halo=.true. )
       call file_history_in( toaflx_lw_up(:,:),        'OLR',          'outgoing longwave  radiation flux', 'W/m2', fill_halo=.true. )
       call file_history_in( toaflx_sw_up(:,:),        'OSR',          'outgoing shortwave radiation flux', 'W/m2', fill_halo=.true. )

       call file_history_in( flux_up(:,:,:,i_lw), 'RADFLUX_LWUP', 'upward   longwave  radiation flux', 'W/m2', fill_halo=.true. )
       call file_history_in( flux_dn(:,:,:,i_lw), 'RADFLUX_LWDN', 'downward longwave  radiation flux', 'W/m2', fill_halo=.true. )
       call file_history_in( flux_net(:,:,:,i_lw), 'RADFLUX_LW',   'net      longwave  radiation flux', 'W/m2', fill_halo=.true. )
       call file_history_in( flux_up(:,:,:,i_sw), 'RADFLUX_SWUP', 'upward   shortwave radiation flux', 'W/m2', fill_halo=.true. )
       call file_history_in( flux_dn(:,:,:,i_sw), 'RADFLUX_SWDN', 'downward shortwave radiation flux', 'W/m2', fill_halo=.true. )
       call file_history_in( flux_net(:,:,:,i_sw), 'RADFLUX_SW',   'net      shortwave radiation flux', 'W/m2', fill_halo=.true. )

       call file_history_in( sflx_rad_dn(:,:,i_r_direct ,i_r_ir ), 'SFLX_IR_dn_dir',  'SFC downward radiation flux (direct; IR)',   'W/m2', fill_halo=.true. )
       call file_history_in( sflx_rad_dn(:,:,i_r_diffuse,i_r_ir ), 'SFLX_IR_dn_dif',  'SFC downward radiation flux (diffuse; IR)',  'W/m2', fill_halo=.true. )
       call file_history_in( sflx_rad_dn(:,:,i_r_direct ,i_r_nir), 'SFLX_NIR_dn_dir', 'SFC downward radiation flux (direct; NIR)',  'W/m2', fill_halo=.true. )
       call file_history_in( sflx_rad_dn(:,:,i_r_diffuse,i_r_nir), 'SFLX_NIR_dn_dif', 'SFC downward radiation flux (diffuse; NIR)', 'W/m2', fill_halo=.true. )
       call file_history_in( sflx_rad_dn(:,:,i_r_direct ,i_r_vis), 'SFLX_VIS_dn_dir', 'SFC downward radiation flux (direct; VIS)',  'W/m2', fill_halo=.true. )
       call file_history_in( sflx_rad_dn(:,:,i_r_diffuse,i_r_vis), 'SFLX_VIS_dn_dif', 'SFC downward radiation flux (diffuse; VIS)', 'W/m2', fill_halo=.true. )

       call file_history_in( temp_t(:,:,:,i_lw), 'TEMP_t_rd_LW', 'tendency of temp in rd(LW)',  'K/day',     fill_halo=.true. )
       call file_history_in( temp_t(:,:,:,i_sw), 'TEMP_t_rd_SW', 'tendency of temp in rd(SW)',  'K/day',     fill_halo=.true. )
       call file_history_in( temp_t(:,:,:,3   ), 'TEMP_t_rd',    'tendency of temp in rd',      'K/day',     fill_halo=.true. )
       call file_history_in( rhoh_rd(:,:,:),      'RHOH_RD',      'diabatic heating rate in rd', 'J/m3/s',    fill_halo=.true. )

       ! output of raw data, for offline output
       call file_history_in( flux_rad(:,:,:,i_lw,i_up,2), 'RFLX_LW_up', 'upward   longwave  radiation flux (cell face)', 'W/m2', fill_halo=.true. )
       call file_history_in( flux_rad(:,:,:,i_lw,i_dn,2), 'RFLX_LW_dn', 'downward longwave  radiation flux (cell face)', 'W/m2', fill_halo=.true. )
       call file_history_in( flux_rad(:,:,:,i_sw,i_up,2), 'RFLX_SW_up', 'upward   shortwave radiation flux (cell face)', 'W/m2', fill_halo=.true. )
       call file_history_in( flux_rad(:,:,:,i_sw,i_dn,2), 'RFLX_SW_dn', 'downward shortwave radiation flux (cell face)', 'W/m2', fill_halo=.true. )

       call file_history_in( dtau_s(:,:,:), 'dtau_s', '0.67 micron cloud optical depth', '1', fill_halo=.true. )
       call file_history_in( dem_s(:,:,:), 'dem_s',  '10.5 micron cloud emissivity',    '1', fill_halo=.true. )

    endif

    do j = js, je
    do i = is, ie
    do k = ks, ke
       rhoh(k,i,j) = rhoh(k,i,j) + rhoh_rd(k,i,j)
    enddo
    enddo
    enddo

    if ( statistics_checktotal ) then
       call statistics_total( ka, ks, ke, ia, is, ie, ja, js, je, &
                              rhoh_rd(:,:,:), 'RHOH_RD',          &
                              atmos_grid_cartesc_real_vol(:,:,:), &
                              atmos_grid_cartesc_real_totvol      )
    endif

    return
  end subroutine atmos_phy_rd_driver_calc_tendency

end module mod_atmos_phy_rd_driver