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_finalize
  public :: atmos_phy_rd_driver_calc_tendency
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  logical, private :: RD_use_PBL_cloud = .false.
 
  !-----------------------------------------------------------------------------
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,   &
       tomflx_lw_up => atmos_phy_rd_tomflx_lw_up, &
       tomflx_lw_dn => atmos_phy_rd_tomflx_lw_dn, &
       tomflx_sw_up => atmos_phy_rd_tomflx_sw_up, &
       tomflx_sw_dn => atmos_phy_rd_tomflx_sw_dn, &
       sflx_rad_dn  => atmos_phy_rd_sflx_down,    &
       solins       => atmos_phy_rd_solins,       &
       cossza       => atmos_phy_rd_cossza
    implicit none
    
    namelist / param_atmos_phy_rd / &
         rd_use_pbl_cloud
 
    integer :: ierr
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_RD_driver_setup",*) 'Setup'
 
    if ( atmos_sw_phy_rd ) then
 
       !--- read namelist
       rewind(io_fid_conf)
       read(io_fid_conf,nml=param_atmos_phy_rd,iostat=ierr)
       if( ierr < 0 ) then !--- missing
          log_info("ATMOS_PHY_RD_driver_setup",*) 'Not found namelist. Default used.'
       elseif( ierr > 0 ) then !--- fatal error
          log_error("ATMOS_PHY_RD_driver_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_RD. Check!'
          call prc_abort
       endif
       log_nml(param_atmos_phy_rd)
 
       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.'
       !$acc kernels
       sfcflx_lw_up(:,:)     = 0.0_rp
       sfcflx_lw_dn(:,:)     = 0.0_rp
       sfcflx_sw_up(:,:)     = 0.0_rp
       sfcflx_sw_dn(:,:)     = 0.0_rp
       tomflx_lw_up(:,:)     = 0.0_rp
       tomflx_lw_dn(:,:)     = 0.0_rp
       tomflx_sw_up(:,:)     = 0.0_rp
       tomflx_sw_dn(:,:)     = 0.0_rp
       sflx_rad_dn(:,:,:,:) = 0.0_rp
       solins(:,:)     = 0.0_rp
       cossza(:,:)     = 0.0_rp
       !$acc end kernels
 
    endif
 
    return
  end subroutine atmos_phy_rd_driver_setup
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_rd_driver_finalize
    use scale_atmos_phy_rd_mstrnx, only: &
       atmos_phy_rd_mstrnx_finalize
    use scale_atmos_phy_rd_profile, only: &
       atmos_phy_rd_profile_finalize
    use mod_atmos_admin, only: &
       atmos_phy_rd_type, &
       atmos_sw_phy_rd
    implicit none
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_RD_driver_finalize",*) 'Finalize'
 
    if ( atmos_sw_phy_rd ) then
       select case ( atmos_phy_rd_type )
       case ( "MSTRNX" )
          call atmos_phy_rd_mstrnx_finalize
       case ( "OFFLINE" )
       end select
    end if
 
    call atmos_phy_rd_profile_finalize
 
    return
  end subroutine atmos_phy_rd_driver_finalize
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_rd_driver_calc_tendency( update_flag )
    use scale_const, only: &
       eps => const_eps
    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, &
       i_hc, &
       i_hi
    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,   &
       tomflx_lw_up => atmos_phy_rd_tomflx_lw_up, &
       tomflx_lw_dn => atmos_phy_rd_tomflx_lw_dn, &
       tomflx_sw_up => atmos_phy_rd_tomflx_sw_up, &
       tomflx_sw_dn => atmos_phy_rd_tomflx_sw_dn, &
       sflx_rad_dn  => atmos_phy_rd_sflx_down,    &
       solins       => atmos_phy_rd_solins,       &
       cossza       => atmos_phy_rd_cossza
    use mod_atmos_phy_bl_vars, only: &
       atmos_phy_bl_zi,     &
       atmos_phy_bl_ql,     &
       atmos_phy_bl_cldfrac
    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_sfc(   ia,ja,2)
    real(rp) :: flux_net_toa(   ia,ja,2)
    real(rp) :: flux_net_tom(   ia,ja,2)
 
    real(rp) :: toaflx_lw_up(ia,ja)
    real(rp) :: toaflx_lw_dn(ia,ja)
    real(rp) :: toaflx_sw_up(ia,ja)
    real(rp) :: toaflx_sw_dn(ia,ja)
 
    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) :: tomflx_lw_up_c(ia,ja)
    real(rp) :: tomflx_lw_dn_c(ia,ja)
    real(rp) :: tomflx_sw_up_c(ia,ja)
    real(rp) :: tomflx_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
 
       !$acc data create(TEMP_t,flux_rad,flux_rad_top,dtau_s,dem_s,flux_up,flux_dn,flux_net,flux_net_sfc,flux_net_toa,flux_net_tom,TOAFLX_LW_up,TOAFLX_LW_dn,TOAFLX_SW_up,TOAFLX_SW_dn,SFCFLX_LW_up_c,SFCFLX_LW_dn_c,SFCFLX_SW_up_c,SFCFLX_SW_dn_c,TOAFLX_LW_up_c,TOAFLX_LW_dn_c,TOAFLX_SW_up_c,TOAFLX_SW_dn_c,TOMFLX_LW_up_c,TOMFLX_LW_dn_c,TOMFLX_SW_up_c,TOMFLX_SW_dn_c,CLDFRAC,MP_Re,MP_Qe,AE_Re,AE_Qe,RH)
 
       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      ) ! [OUT]
 
       if ( rd_use_pbl_cloud ) then
          !$omp parallel do
          !$acc kernels
          do j = js, je
          do i = is, ie
          do k = ks, ke
             if ( real_cz(k,i,j) < atmos_phy_bl_zi(i,j) + real_fz(ks-1,i,j) ) then
                cldfrac(k,i,j) = atmos_phy_bl_cldfrac(k,i,j)
                if ( atmos_phy_bl_cldfrac(k,i,j) > eps ) then
                   mp_qe(k,i,j,i_hc) = atmos_phy_bl_ql(k,i,j) / atmos_phy_bl_cldfrac(k,i,j)
                else
                   mp_qe(k,i,j,i_hc) = 0.0_rp
                end if
                mp_qe(k,i,j,i_hi) = 0.0_rp
                mp_re(k,i,j,i_hc) = 8.0e-4
             end if
          end do
          end do
          end do
          !$acc end kernels
          
       end if
 
       call atmos_vars_get_diagnostic( "RH", rh )
       call atmos_phy_ae_vars_get_diagnostic( &
            qtrc(:,:,:,:), rh(:,:,:), & ! [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(:,:,:,:), 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]
          !$acc kernels
          flux_rad(:,:,:,:,:,1)   = 0.0_rp ! clear sky
          flux_rad_top(:,:,:,:,:) = 0.0_rp
          dtau_s(:,:,:) = 0.0_rp
          dem_s(:,:,:) = 0.0_rp
          !$acc end kernels
 
       end select
 
 
       ! surface
!OCL XFILL
       !$omp parallel do
       !$acc kernels
       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
       !$acc end kernels
 
       ! top of the atmosphere
!OCL XFILL
       !$omp parallel do
       !$acc kernels
       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
       !$acc end kernels
 
       ! top of the model
!OCL XFILL
       !$omp parallel do
       !$acc kernels
       do j = js, je
       do i = is, ie
          ! for clear-sky
          tomflx_lw_up_c(i,j)    = flux_rad(ke,i,j,i_lw,i_up,1)
          tomflx_lw_dn_c(i,j)    = flux_rad(ke,i,j,i_lw,i_dn,1)
          tomflx_sw_up_c(i,j)    = flux_rad(ke,i,j,i_sw,i_up,1)
          tomflx_sw_dn_c(i,j)    = flux_rad(ke,i,j,i_sw,i_dn,1)
          ! for all-sky
          tomflx_lw_up(i,j)    = flux_rad(ke,i,j,i_lw,i_up,2)
          tomflx_lw_dn(i,j)    = flux_rad(ke,i,j,i_lw,i_dn,2)
          tomflx_sw_up(i,j)    = flux_rad(ke,i,j,i_sw,i_up,2)
          tomflx_sw_dn(i,j)    = flux_rad(ke,i,j,i_sw,i_dn,2)
 
          flux_net_tom(i,j,i_lw) = tomflx_lw_up(i,j) - tomflx_lw_dn(i,j)
          flux_net_tom(i,j,i_sw) = tomflx_sw_up(i,j) - tomflx_sw_dn(i,j)
       enddo
       enddo
       !$acc end kernels
 
!OCL XFILL
       !$omp parallel do collapse(2)
       !$acc kernels
       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
       !$acc end kernels
 
       ! 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( tomflx_lw_up_c(:,:),      'TOMFLX_LW_up_c', 'TOM upward   longwave  radiation flux (clr)', 'W/m2', fill_halo=.true. )
       call file_history_in( tomflx_lw_dn_c(:,:),      'TOMFLX_LW_dn_c', 'TOM downward longwave  radiation flux (clr)', 'W/m2', fill_halo=.true. )
       call file_history_in( tomflx_sw_up_c(:,:),      'TOMFLX_SW_up_c', 'TOM upward   shortwave radiation flux (clr)', 'W/m2', fill_halo=.true. )
       call file_history_in( tomflx_sw_dn_c(:,:),      'TOMFLX_SW_dn_c', 'TOM downward shortwave radiation flux (clr)', 'W/m2', fill_halo=.true. )
 
       call file_history_in( tomflx_lw_up(:,:),      'TOMFLX_LW_up',   'TOM upward   longwave  radiation flux',       'W/m2', fill_halo=.true. )
       call file_history_in( tomflx_lw_dn(:,:),      'TOMFLX_LW_dn',   'TOM downward longwave  radiation flux',       'W/m2', fill_halo=.true. )
       call file_history_in( tomflx_sw_up(:,:),      'TOMFLX_SW_up',   'TOM upward   shortwave radiation flux',       'W/m2', fill_halo=.true. )
       call file_history_in( tomflx_sw_dn(:,:),      'TOMFLX_SW_dn',   'TOM downward shortwave radiation flux',       'W/m2', fill_halo=.true. )
 
       call file_history_in( flux_net_tom(:,:,i_lw), 'TOMFLX_LW_net',  'TOM net      longwave  radiation flux',       'W/m2', fill_halo=.true. )
       call file_history_in( flux_net_tom(:,:,i_sw), 'TOMFLX_SW_net',  'TOM 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. )
 
       !$acc end data
 
    endif
 
    !$acc kernels
    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
    !$acc end kernels
 
    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