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_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 mod_atmos_vars, only: &
       temp,              &
       pres,              &
       dens,              &
       rhot,              &
       qtrc,              &
       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) = flux_rad(ks-1,i,j,i_lw,i_up,2) - flux_rad(ks-1,i,j,i_lw,i_dn,2)
          flux_net_sfc(i,j,i_sw) = flux_rad(ks-1,i,j,i_sw,i_up,2) - flux_rad(ks-1,i,j,i_sw,i_dn,2)
       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) = flux_rad_top(i,j,i_lw,i_up,2) - flux_rad_top(i,j,i_lw,i_dn,2)
          flux_net_toa(i,j,i_sw) = flux_rad_top(i,j,i_sw,i_up,2) - flux_rad_top(i,j,i_sw,i_dn,2)
       enddo
       enddo

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

          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_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_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_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

       if ( atmos_phy_rd_type == 'WRF' ) then

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

          do j = 1, ja
          do i = 1, ia
          do k = 1, ka
             qv3d(i,k,j) = qtrc(k,i,j,i_qv)
             qc3d(i,k,j) = qtrc(k,i,j,i_qc)
             qr3d(i,k,j) = qtrc(k,i,j,i_qr)
             qi3d(i,k,j) = qtrc(k,i,j,i_qi)
             qs3d(i,k,j) = qtrc(k,i,j,i_qs)
             qg3d(i,k,j) = qtrc(k,i,j,i_qg)
             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)', '0-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( 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( 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( sflx_rad_dn(:,:,i_sw,i_direct), 'SFLX_SW_dn_direct',  'SFC down. shortwave ratio direct',  'W/m2', nohalo=.true. )
       call hist_in( sflx_rad_dn(:,:,i_sw,i_diffuse),'SFLX_SW_dn_diffuse', 'SFC down. shortwave ratio diffuse', '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( 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. )

       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