mod_atmos_phy_rd_driver.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "inc_openmp.h"
module mod_atmos_phy_rd_driver
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_tracer
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_rd_driver_setup
  public :: atmos_phy_rd_driver_resume
  public :: atmos_phy_rd_driver

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_rd_driver_setup
    use scale_atmos_phy_rd, only: &
       atmos_phy_rd_setup
    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_downall, &
       solins       => atmos_phy_rd_solins,       &
       cossza       => atmos_phy_rd_cossza
    implicit none
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[DRIVER] / Categ[ATMOS PHY_RD] / Origin[SCALE-RM]'

    if ( atmos_sw_phy_rd ) then

       ! setup library component
       call atmos_phy_rd_setup( atmos_phy_rd_type )

    else

       if( io_l ) write(io_fid_log,*) '*** this component is never called.'
       if( io_l ) write(io_fid_log,*) '*** 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_resume
    use mod_atmos_admin, only: &
       atmos_sw_phy_rd
    implicit none

    if ( atmos_sw_phy_rd ) then

       ! run once (only for the diagnostic value)
       call prof_rapstart('ATM_Radiation', 1)
       call atmos_phy_rd_driver( update_flag = .true. )
       call prof_rapend  ('ATM_Radiation', 1)

    end if

    return
  end subroutine atmos_phy_rd_driver_resume

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_rd_driver( update_flag )
    use scale_grid_real, only: &
       real_cz,            &
       real_fz,            &
       real_lon,           &
       real_lat
    use scale_process, only: &
       prc_mpistop
    use scale_const, only: &
       pre00 => const_pre00, &
       rdry  => const_rdry,  &
       cpdry => const_cpdry
    use scale_landuse, only: &
       landuse_fact_ocean, &
       landuse_fact_land,  &
       landuse_fact_urban
    use scale_time, only: &
       dt_rd => time_dtsec_atmos_phy_rd, &
       time_nowdate,                     &
       time_offset_year
    use scale_rm_statistics, only: &
       statistics_checktotal, &
       stat_total
    use scale_history, only: &
       hist_in
    use mod_atmos_admin, only: &
       atmos_phy_rd_type
    use scale_atmos_solarins, only: &
       solarins_insolation => atmos_solarins_insolation
    use scale_atmos_phy_rd, only: &
       atmos_phy_rd
    use scale_atmos_phy_rd_mm5sw, only: &
       swrad
    use scale_atmos_phy_rd_common, only: &
       rd_heating => atmos_phy_rd_heating, &
       i_sw,     &
       i_lw,     &
       i_dn,     &
       i_up,     &
       i_direct, &
       i_diffuse
    use scale_atmos_hydrometeor, only: &
       i_qv, &
       i_qc, &
       i_qr, &
       i_qi, &
       i_qs, &
       i_qg
    use mod_atmos_vars, only: &
       temp,              &
       pres,              &
       dens   => dens_av, &
       rhot   => rhot_av, &
       qtrc   => qtrc_av, &
       rhot_t => rhot_tp
    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: &
       rhot_t_rd    => atmos_phy_rd_rhot_t,       &
       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_downall, &
       solins       => atmos_phy_rd_solins,       &
       cossza       => atmos_phy_rd_cossza
    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) :: 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)

    ! for WRF radiation scheme added by Adachi; array order is (i,k,j)
    real(RP) :: rthratensw(ia,ka,ja)
    real(RP) :: sdown3d   (ia,ka,ja)  ! downward short wave flux (W/m2)
    real(RP) :: gsw       (ia,ja)     ! net short wave flux at ground surface (W/m2)
    real(RP) :: rho3d     (ia,ka,ja)
    real(RP) :: t3d       (ia,ka,ja)
    real(RP) :: p3d       (ia,ka,ja)
    real(RP) :: pi3d      (ia,ka,ja)
    real(RP) :: dz8w      (ia,ka,ja)
    real(RP) :: qv3d      (ia,ka,ja)
    real(RP) :: qc3d      (ia,ka,ja)
    real(RP) :: qr3d      (ia,ka,ja)
    real(RP) :: qi3d      (ia,ka,ja)
    real(RP) :: qs3d      (ia,ka,ja)
    real(RP) :: qg3d      (ia,ka,ja)
    real(RP) :: flux_rad_org(ka,ia,ja,2,2,2)

    real(RP) :: total ! dummy

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

    if ( update_flag ) then

       call solarins_insolation( solins(:,:),   & ! [OUT]
                                 cossza(:,:),   & ! [OUT]
                                 real_lon(:,:),   & ! [IN]
                                 real_lat(:,:),   & ! [IN]
                                 time_nowdate(:), & ! [IN]
                                 time_offset_year ) ! [IN]


       call atmos_phy_rd( dens, rhot, qtrc,   & ! [IN]
                          real_cz, real_fz,   & ! [IN]
                          landuse_fact_ocean, & ! [IN]
                          landuse_fact_land,  & ! [IN]
                          landuse_fact_urban, & ! [IN]
                          sfc_temp,           & ! [IN]
                          sfc_albedo,         & ! [IN]
                          solins, cossza,     & ! [IN]
                          flux_rad,           & ! [OUT]
                          flux_rad_top,       & ! [OUT]
                          sflx_rad_dn         ) ! [OUT]


       ! apply radiative flux convergence -> heating rate
       call rd_heating( flux_rad(:,:,:,:,:,2), & ! [IN]
                        dens(:,:,:),       & ! [IN]
                        rhot(:,:,:),       & ! [IN]
                        qtrc(:,:,:,:),     & ! [IN]
                        real_fz(:,:,:),       & ! [IN]
                        dt_rd,                  & ! [IN]
                        temp_t(:,:,:,:),     & ! [OUT]
                        rhot_t_rd(:,:,:)        ) ! [OUT]


!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

       if ( atmos_phy_rd_type == 'WRF' ) then

          flux_rad_org(:,:,:,:,:,:) = flux_rad(:,:,:,:,:,:)
          rthratensw = 0.0_rp
          sdown3d    = 0.0_rp
          gsw        = 0.0_rp

          if ( i_qv > 0 ) then
             qv3d(:,:,:) = qtrc(:,:,:,i_qv)
          else
             qv3d = 0.0_rp
          end if
          if ( i_qc > 0 ) then
             qc3d(:,:,:) = qtrc(:,:,:,i_qc)
          else
             qc3d = 0.0_rp
          end if
          if ( i_qr > 0 ) then
             qr3d(:,:,:) = qtrc(:,:,:,i_qr)
          else
             qr3d = 0.0_rp
          end if
          if ( i_qi > 0 ) then
             qi3d(:,:,:) = qtrc(:,:,:,i_qi)
          else
             qi3d = 0.0_rp
          end if
          if ( i_qs > 0 ) then
             qs3d(:,:,:) = qtrc(:,:,:,i_qs)
          else
             qs3d = 0.0_rp
          end if
          if ( i_qg > 0 ) then
             qg3d(:,:,:) = qtrc(:,:,:,i_qg)
          else
             qg3d = 0.0_rp
          end if

          do j = 1, ja
          do i = 1, ia
          do k = 1, ka
             t3d(i,k,j) = temp(k,i,j)                       ! temperature (K)
             rho3d(i,k,j) = dens(k,i,j)                       ! density (kg/m^3)
             p3d(i,k,j) = pres(k,i,j)                       ! pressure (Pa)
             pi3d(i,k,j) = (pres(k,i,j)/pre00)**(rdry/cpdry) ! exner function (dimensionless)
             dz8w(i,k,j) = real_fz(k,i,j)-real_fz(k-1,i,j)   ! dz between full levels(m)
          enddo
          enddo
          enddo

          call swrad( dt_rd,                & ! [IN]
                      rthratensw,           & ! [INOUT]
                      sdown3d,              & ! [INOUT]
                      gsw,                  & ! [INOUT]
                      real_lat,             & ! [IN]
                      real_lon,             & ! [IN]
                      sfc_albedo(:,:,i_sw), & ! [IN]
                      rho3d,                & ! [IN]
                      t3d,                  & ! [IN]
                      p3d,                  & ! [IN]
                      pi3d,                 & ! [IN]
                      dz8w,                 & ! [IN]
                      solins(:,:),          & ! [IN]
                      cossza(:,:),          & ! [IN]
                      qv3d,                 & ! [IN]
                      qc3d,                 & ! [IN]
                      qr3d,                 & ! [IN]
                      qi3d,                 & ! [IN]
                      qs3d,                 & ! [IN]
                      qg3d,                 & ! [IN]
                      f_qv      = .true.,   & ! [IN]
                      f_qc      = .true.,   & ! [IN]
                      f_qr      = .true.,   & ! [IN]
                      f_qi      = .true.,   & ! [IN]
                      f_qs      = .true.,   & ! [IN]
                      f_qg      = .true.,   & ! [IN]
                      icloud    = 1,        & ! [IN]
                      warm_rain = .true.    ) ! [IN]

          do j = js, je
          do i = is, ie
             flux_net_sfc(i,j,i_sw) = gsw(i,j)
             do k = ks-1, ke
                flux_rad(k,i,j,i_sw,i_up,2) = 0.0_rp
                flux_rad(k,i,j,i_sw,i_dn,2) = sdown3d(i,k,j)
             enddo

             do k = 1, ks-2
                flux_rad(k,i,j,i_sw,i_dn,2) = sdown3d(i,ks-1,j)
             enddo

             do k = ke+1, ka
                flux_rad(k,i,j,i_sw,i_dn,2) = sdown3d(i,ke,j)
             enddo
          enddo
          enddo

          do j = js, je
          do i = is, ie
             sfcflx_sw_up(i,j) = flux_rad(ks-1,i,j,i_sw,i_dn,2) * sfc_albedo(i,j,i_sw)
             sfcflx_sw_dn(i,j) = flux_rad(ks-1,i,j,i_sw,i_dn,2)

             toaflx_sw_up(i,j) = 0.0_rp
             toaflx_sw_dn(i,j) = flux_rad(ke,i,j,i_sw,i_dn,2) ! sometimes TOA altitude is very low
          enddo
          enddo

          do j = js, je
          do i = is, ie
             flux_net_toa(i,j,i_sw) = 0.0_rp
             do k = ks, ke
                flux_net(k,i,j,i_sw) = 0.0_rp
                flux_up(k,i,j,i_sw) = 0.0_rp
                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) )
             enddo
          enddo
          enddo

          do j = js, je
          do i = is, ie
             flux_net_sfc(i,j,i_sw) = flux_rad(ks-1,i,j,i_sw,i_dn,2)*sfc_albedo(i,j,i_sw)-flux_rad(ks-1,i,j,i_sw,i_dn,2)
          enddo
          enddo

       endif

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

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

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

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

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

       call hist_in( flux_net_sfc(:,:,i_lw),   'SLR',          'SFC net longwave  radiation flux',  'W/m2', nohalo=.true. )
       call hist_in( flux_net_sfc(:,:,i_sw),   'SSR',          'SFC net shortwave radiation flux',  'W/m2', nohalo=.true. )
       call hist_in( flux_net_toa(:,:,i_lw),   'OLR',          'TOA net longwave  radiation flux',  'W/m2', nohalo=.true. )
       call hist_in( flux_net_toa(:,:,i_sw),   'OSR',          'TOA net shortwave radiation flux',  'W/m2', nohalo=.true. )

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

       call hist_in( sflx_rad_dn(:,:,i_lw,i_direct ), 'SFLX_LW_dn_dir', 'SFC downward longwave  flux (direct )', 'W/m2', nohalo=.true. )
       call hist_in( sflx_rad_dn(:,:,i_lw,i_diffuse), 'SFLX_LW_dn_dif', 'SFC downward longwave  flux (diffuse)', 'W/m2', nohalo=.true. )
       call hist_in( sflx_rad_dn(:,:,i_sw,i_direct ), 'SFLX_SW_dn_dir', 'SFC downward shortwave flux (direct )', 'W/m2', nohalo=.true. )
       call hist_in( sflx_rad_dn(:,:,i_sw,i_diffuse), 'SFLX_SW_dn_dif', 'SFC downward shortwave flux (diffuse)', 'W/m2', nohalo=.true. )

       call hist_in( temp_t(:,:,:,i_lw), 'TEMP_t_rd_LW', 'tendency of temp in rd(LW)', 'K/day',     nohalo=.true. )
       call hist_in( temp_t(:,:,:,i_sw), 'TEMP_t_rd_SW', 'tendency of temp in rd(SW)', 'K/day',     nohalo=.true. )
       call hist_in( temp_t(:,:,:,3   ), 'TEMP_t_rd',    'tendency of temp in rd',     'K/day',     nohalo=.true. )
       call hist_in( rhot_t_rd(:,:,:),      'RHOT_t_RD',    'tendency of RHOT in rd',     'K.kg/m3/s', nohalo=.true. )

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

       if ( atmos_phy_rd_type ==  'WRF' ) then
          ! revert all radiation flux from MM5 scheme to default
          flux_rad(:,:,:,:,:,:) = flux_rad_org(:,:,:,:,:,:)

          do j = js, je
          do i = is, ie
             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)

             toaflx_sw_up(i,j) = flux_rad_top(i,j,i_sw,i_up,2) ! mstrnx
             toaflx_sw_dn(i,j) = flux_rad_top(i,j,i_sw,i_dn,2) ! mstrnx
          enddo
          enddo

       endif

    endif

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

    if ( statistics_checktotal ) then
       call stat_total( total, rhot_t_rd(:,:,:), 'RHOT_t_RD' )
    endif

    return
  end subroutine atmos_phy_rd_driver

end module mod_atmos_phy_rd_driver