!= 放射関連ルーチン
!
!= Routines for radiation calculation
!
! Authors::   Yasuhiro MORIKAWA, Yukiko YAMADA, Satoshi NODA, Yoshiyuki O. TAKAHASHI
! Version::   $Id: rad_utils.f90,v 1.3 2011-04-06 15:28:06 yot Exp $ 
! Tag Name::  $Name: dcpam5-20110407 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module rad_utils
  !
  != 放射関連ルーチン
  !
  != Routines for radiation calculation
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 
  !
  ! 
  !
  !== Procedures List
  !
  ! RadDTempDt        :: 放射フラックスによる温度変化の計算
  ! RadFluxOutput     :: 放射フラックスの出力
  ! ------------            :: ------------
  ! RadDTempDt        :: Calculate temperature tendency with radiation flux
  ! RadFluxOutput     :: Output radiation fluxes
  !
  !== NAMELIST
  !
  ! NAMELIST#rad_utils_nml
  !

  ! モジュール引用 ; USE statements
  !

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STRING, &  ! 文字列.       Strings. 
    &                 TOKEN      ! キーワード.   Keywords. 

  ! 物理定数設定
  ! Physical constants settings
  !
  use constants, only: &
    & PI, &
                            ! $ \pi $ .
                            ! 円周率.  Circular constant
    & StB, &
                            ! $ \sigma_{SB} $ .
                            ! ステファンボルツマン定数.
                            ! Stefan-Boltzmann constant
    & Grav, &
                            ! $ g $ [m s-2]. 
                            ! 重力加速度. 
                            ! Gravitational acceleration
    & CpDry
                            ! $ C_p $ [J kg-1 K-1]. 
                            ! 乾燥大気の定圧比熱. 
                            ! Specific heat of air at constant pressure


  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: imax, & ! 経度格子点数. 
                             ! Number of grid points in longitude
    &                jmax, & ! 緯度格子点数. 
                             ! Number of grid points in latitude
    &                kmax    ! 鉛直層数. 
                             ! Number of vertical level

  ! メッセージ出力
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public :: RadRTEQNonScat
  public :: RadRTEQNonScatWrapper
  public :: RadDTempDt
  public :: RadDTempDtforNHM2DWrapper
  public :: RadFluxOutput

  ! 公開変数
  ! Public variables
  !
  logical, save, public:: rad_utils_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag


  ! 非公開変数
  ! Private variables
  !
  real(DP), save :: DiffFact


  character(*), parameter:: module_name = 'rad_utils'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name: dcpam5-20110407 $' // &
    & '$Id: rad_utils.f90,v 1.3 2011-04-06 15:28:06 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

contains

  !--------------------------------------------------------------------------------------

!!$  subroutine RadiationRTEQNonScat(                    &
!!$    & xyz_Temp, xy_SurfTemp, xy_SurfEmis, xyr_OptDep, & ! (in )
!!$    & xyr_RadLFlux, xyra_DelRadLFlux,                 & ! (out)
!!$    & WNs, WNe, NumGaussNode,                         & ! (in ) optional
!!$    & xyrr_Trans                                      & ! (in ) optional
!!$    & )
!!$    !
!!$    ! 散乱なしの場合の放射伝達方程式の計算
!!$    !
!!$    ! Integrate radiative transfer equation without scattering
!!$    !
!!$
!!$    ! モジュール引用 ; USE statements
!!$    !
!!$
!!$    ! 物理定数設定
!!$    ! Physical constants settings
!!$    !
!!$    use constants, only: &
!!$      & PI, &
!!$                              ! $ \pi $ .
!!$                              ! 円周率.  Circular constant
!!$      & StB
!!$                              ! $ \sigma_{SB} $ .
!!$                              ! ステファンボルツマン定数.
!!$                              ! Stefan-Boltzmann constant
!!$
!!$
!!$    ! プランク関数の計算
!!$    ! Calculate Planck function
!!$    !
!!$    use planck_func, only :                            &
!!$      & Integ_PF_GQ_Array3D   , Integ_PF_GQ_Array2D,   &
!!$      & Integ_DPFDT_GQ_Array3D, Integ_DPFDT_GQ_Array2D
!!$
!!$    ! 宣言文 ; Declaration statements
!!$    !
!!$
!!$    real(DP), intent(in ):: xyz_Temp    (0:imax-1, 1:jmax, 1:kmax)
!!$                              ! $ T $ .     温度. Temperature
!!$    real(DP), intent(in ):: xy_SurfTemp (0:imax-1, 1:jmax)
!!$                              ! 地表面温度. 
!!$                              ! Surface temperature
!!$    real(DP), intent(in ):: xy_SurfEmis (0:imax-1, 1:jmax)
!!$                              ! 惑星表面射出率. 
!!$                              ! Surface emissivity
!!$    real(DP), intent(in ):: xyr_OptDep  (0:imax-1, 1:jmax, 0:kmax)
!!$                              ! Optical depth
!!$    real(DP), intent(out):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
!!$                              ! 長波フラックス. 
!!$                              ! Longwave flux
!!$    real(DP), intent(out):: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
!!$                              ! 長波地表温度変化. 
!!$                              ! Surface temperature tendency with longwave
!!$    real(DP), intent(in ), optional :: WNs
!!$    real(DP), intent(in ), optional :: WNe
!!$    integer , intent(in ), optional :: NumGaussNode
!!$    real(DP), intent(in ), optional :: xyrr_Trans   (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
!!$                              ! 透過係数. 
!!$                              ! Transmission coefficient
!!$
!!$
!!$    ! 作業変数
!!$    ! Work variables
!!$    !
!!$    real(DP):: xyr_RadLDoFlux (0:imax-1, 1:jmax, 0:kmax)
!!$    real(DP):: xyr_RadLUpFlux (0:imax-1, 1:jmax, 0:kmax)
!!$
!!$    real(DP):: xyz_DelTrans (0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP):: xyrr_TransLV (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
!!$                              ! 透過係数. 
!!$                              ! Transmission coefficient
!!$    real(DP):: xyz_IntPF        (0:imax-1, 1:jmax, 1:kmax)
!!$                              ! Integrated Planck function
!!$    real(DP):: xy_SurfIntPF     (0:imax-1, 1:jmax)
!!$                              ! Integrated Planck function with surface temperature
!!$    real(DP):: xyz_IntDPFDT     (0:imax-1, 1:jmax, 1:kmax)
!!$                              ! Integrated temperature derivative of Planck function
!!$    real(DP):: xy_SurfIntDPFDT  (0:imax-1, 1:jmax)
!!$                              ! Integrated temperature derivative of Planck function
!!$                              ! with surface temperature
!!$
!!$    integer:: k, kk           ! 鉛直方向に回る DO ループ用作業変数
!!$                              ! Work variables for DO loop in vertical direction
!!$
!!$    ! 実行文 ; Executable statement
!!$    !
!!$
!!$    ! 初期化
!!$    ! Initialization
!!$    !
!!$    if ( .not. radiation_utils_inited ) call RadiationInit
!!$
!!$    ! Check arguments
!!$    !
!!$    if ( present( WNs ) ) then
!!$      if ( ( .not. present( WNe ) ) .or. ( .not. present( NumGaussNode ) ) ) then
!!$        call MessageNotify( 'E', module_name, &
!!$          & 'All of WNs, WNe, and NumGaussNode have to be present.' )
!!$      end if
!!$    else
!!$      if ( present( WNe ) .or. present( NumGaussNode ) ) then
!!$        call MessageNotify( 'E', module_name, &
!!$          & 'All of WNs, WNe, and NumGaussNode have to be present.' )
!!$      end if
!!$    end if
!!$
!!$    if ( present( WNs ) ) then
!!$      ! Case for non-grey atmosphere
!!$      !
!!$
!!$      ! Integrate Planck function and temperature derivative of it
!!$      !
!!$      call Integ_PF_GQ_Array3D( &
!!$        & WNs, WNe, NumGaussNode, &
!!$        & 0, imax-1, 1, jmax, 1, kmax, &
!!$        & xyz_Temp, &
!!$        & xyz_IntPF &
!!$        & )
!!$      call Integ_PF_GQ_Array2D( &
!!$        & WNs, WNe, NumGaussNode, &
!!$        & 0, imax-1, 1, jmax, &
!!$        & xy_SurfTemp, &
!!$        & xy_SurfIntPF &
!!$        & )
!!$      call Integ_DPFDT_GQ_Array3D( &
!!$        & 0, imax-1, 1, jmax, 1, kmax, &
!!$        & WNs, WNe, NumGaussNode, xyz_Temp, & ! (in )
!!$        & xyz_IntDPFDT          & ! (out)
!!$        & )
!!$      call Integ_DPFDT_GQ_Array2D( &
!!$        & 0, imax-1, 1, jmax, &
!!$        & WNs, WNe, NumGaussNode, xy_SurfTemp, & ! (in )
!!$        & xy_SurfIntDPFDT           & ! (out)
!!$        & )
!!$
!!$      xyz_IntPF       = PI * xyz_IntPF
!!$      xy_SurfIntPF    = PI * xy_SurfIntPF
!!$      xyz_IntDPFDT    = PI * xyz_IntDPFDT
!!$      xy_SurfIntDPFDT = PI * xy_SurfIntDPFDT
!!$
!!$    else
!!$
!!$      ! Case for grey atmosphere
!!$      !
!!$      xyz_IntPF       = StB * xyz_Temp**4
!!$      xy_SurfIntPF    = StB * xy_SurfTemp**4
!!$      xyz_IntDPFDT    = 4.0_DP * StB * xyz_Temp**3
!!$      xy_SurfIntDPFDT = 4.0_DP * StB * xy_SurfTemp**3
!!$
!!$    end if
!!$
!!$
!!$    if ( present( xyrr_Trans ) ) then
!!$      do k = 0, kmax
!!$        do kk = k, k
!!$          xyrr_TransLV(:,:,k,kk) = 1.0_DP
!!$        end do
!!$        do kk = k+1, kmax
!!$          xyrr_TransLV(:,:,k,kk) = xyrr_Trans(:,:,k,kk)
!!$        end do
!!$      end do
!!$    else
!!$      ! 透過関数計算
!!$      ! Calculate transmission functions
!!$      !
!!$      do k = 1, kmax
!!$        xyz_DelTrans(:,:,k) = &
!!$          & exp( - DiffFact * ( xyr_OptDep(:,:,k-1) - xyr_OptDep(:,:,k) ) )
!!$      end do
!!$      !
!!$      do k = 0, kmax
!!$        do kk = k, k
!!$          xyrr_TransLV(:,:,k,kk) = 1.0_DP
!!$        end do
!!$        do kk = k+1, kmax
!!$          xyrr_TransLV(:,:,k,kk) = xyrr_TransLV(:,:,k,kk-1) * xyz_DelTrans(:,:,kk)
!!$        end do
!!$      end do
!!$    end if
!!$    do k = 0, kmax
!!$      do kk = 0, k-1
!!$        xyrr_TransLV(:,:,k,kk) = xyrr_TransLV(:,:,kk,k)
!!$      end do
!!$    end do
!!$
!!$
!!$    ! 放射フラックス計算
!!$    ! Calculate radiation flux
!!$    !
!!$!    do k = 0, kmax
!!$!
!!$!      xyr_RadLFlux(:,:,k) = xy_SurfEmis * PI * xy_SurfIntPF * xyrr_Trans(:,:,k,0)
!!$!
!!$!      do kk = 0, kmax-1
!!$!        xyr_RadLFlux(:,:,k) = xyr_RadLFlux(:,:,k)                &
!!$!          & - PI * xyz_IntPF(:,:,kk+1)                           &
!!$!          & * ( xyrr_Trans(:,:,k,kk) - xyrr_Trans(:,:,k,kk+1) )
!!$!      end do
!!$!
!!$!    end do
!!$
!!$
!!$    !   Initialization
!!$    !
!!$    xyr_RadLDoFlux = 0.0_DP
!!$    xyr_RadLUpFlux = 0.0_DP
!!$    !
!!$    !   Downward flux
!!$    !
!!$    do k = kmax, 0, -1
!!$
!!$      do kk = kmax, k+1, -1
!!$        xyr_RadLDoFlux(:,:,k) = xyr_RadLDoFlux(:,:,k)          &
!!$          & + xyz_IntPF(:,:,kk)                                &
!!$          & * ( xyrr_TransLV(:,:,k,kk-1) - xyrr_TransLV(:,:,k,kk) )
!!$      end do
!!$
!!$    end do
!!$    !
!!$    !   Upward flux
!!$    !
!!$    !     Set upward flux
!!$    !
!!$    do k = 0, kmax
!!$
!!$      xyr_RadLUpFlux(:,:,k) = xy_SurfEmis * xy_SurfIntPF * xyrr_TransLV(:,:,k,0)
!!$
!!$      do kk = 1, k
!!$        xyr_RadLUpFlux(:,:,k) = xyr_RadLUpFlux(:,:,k)          &
!!$          & - xyz_IntPF(:,:,kk)                                &
!!$          & * ( xyrr_TransLV(:,:,k,kk-1) - xyrr_TransLV(:,:,k,kk) )
!!$      end do
!!$
!!$    end do
!!$
!!$    xyr_RadLFlux = xyr_RadLUpFlux - xyr_RadLDoFlux
!!$
!!$
!!$    ! 放射フラックスの変化率の計算
!!$    ! Calculate rate of change of radiative flux
!!$    !
!!$    do k = 0, kmax
!!$      xyra_DelRadLFlux(:,:,k,0) =                           &
!!$        & xy_SurfEmis * xy_SurfIntDPFDT * xyrr_TransLV(:,:,k,0)
!!$
!!$      xyra_DelRadLFlux(:,:,k,1) =                           &
!!$        & - xyz_IntDPFDT(:,:,1)                             &
!!$        &   * ( xyrr_TransLV(:,:,k,0) - xyrr_TransLV(:,:,k,1) )
!!$    end do
!!$
!!$
!!$  end subroutine RadiationRTEQNonScat

  !--------------------------------------------------------------------------------------
!!$
!!$  subroutine RadiationRTEQNonScat(                    &
!!$    & xyz_Temp, xy_SurfTemp, xy_SurfEmis, xyrr_Trans, & ! (in )
!!$    & xyr_RadLFlux, xyra_DelRadLFlux,                 & ! (out)
!!$    & WNs, WNe, NumGaussNode                          & ! (in ) optional
!!$    & )
!!$    !
!!$    ! 散乱なしの場合の放射伝達方程式の計算
!!$    !
!!$    ! Integrate radiative transfer equation without scattering
!!$    !
!!$
!!$    ! モジュール引用 ; USE statements
!!$    !
!!$
!!$    ! 物理定数設定
!!$    ! Physical constants settings
!!$    !
!!$    use constants, only: &
!!$      & PI, &
!!$                              ! $ \pi $ .
!!$                              ! 円周率.  Circular constant
!!$      & StB
!!$                              ! $ \sigma_{SB} $ .
!!$                              ! ステファンボルツマン定数.
!!$                              ! Stefan-Boltzmann constant
!!$
!!$
!!$    ! プランク関数の計算
!!$    ! Calculate Planck function
!!$    !
!!$    use planck_func, only :                            &
!!$      & Integ_PF_GQ_Array3D   , Integ_PF_GQ_Array2D,   &
!!$      & Integ_DPFDT_GQ_Array3D, Integ_DPFDT_GQ_Array2D
!!$
!!$    ! 宣言文 ; Declaration statements
!!$    !
!!$
!!$    real(DP), intent(in ) :: xyz_Temp    (0:imax-1, 1:jmax, 1:kmax)
!!$                              ! $ T $ .     温度. Temperature
!!$    real(DP), intent(in ) :: xy_SurfTemp (0:imax-1, 1:jmax)
!!$                              ! 地表面温度. 
!!$                              ! Surface temperature
!!$    real(DP), intent(in ) :: xy_SurfEmis (0:imax-1, 1:jmax)
!!$                              ! 惑星表面射出率. 
!!$                              ! Surface emissivity
!!$    real(DP), intent(in ) :: xyrr_Trans   (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
!!$                              ! 透過係数. 
!!$                              ! Transmission coefficient
!!$    real(DP), intent(out) :: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
!!$                              ! 長波フラックス. 
!!$                              ! Longwave flux
!!$    real(DP), intent(out) :: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
!!$                              ! 長波地表温度変化. 
!!$                              ! Surface temperature tendency with longwave
!!$    real(DP), intent(in ), optional :: WNs
!!$    real(DP), intent(in ), optional :: WNe
!!$    integer , intent(in ), optional :: NumGaussNode
!!$
!!$
!!$    ! 作業変数
!!$    ! Work variables
!!$    !
!!$    real(DP):: xyr_RadLDoFlux (0:imax-1, 1:jmax, 0:kmax)
!!$    real(DP):: xyr_RadLUpFlux (0:imax-1, 1:jmax, 0:kmax)
!!$
!!$    real(DP):: xyz_DelTrans (0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP):: xyz_IntPF        (0:imax-1, 1:jmax, 1:kmax)
!!$                              ! Integrated Planck function
!!$    real(DP):: xy_SurfIntPF     (0:imax-1, 1:jmax)
!!$                              ! Integrated Planck function with surface temperature
!!$    real(DP):: xyz_IntDPFDT     (0:imax-1, 1:jmax, 1:kmax)
!!$                              ! Integrated temperature derivative of Planck function
!!$    real(DP):: xy_SurfIntDPFDT  (0:imax-1, 1:jmax)
!!$                              ! Integrated temperature derivative of Planck function
!!$                              ! with surface temperature
!!$
!!$    integer:: k, kk           ! 鉛直方向に回る DO ループ用作業変数
!!$                              ! Work variables for DO loop in vertical direction
!!$
!!$    ! 実行文 ; Executable statement
!!$    !
!!$
!!$    ! 初期化
!!$    ! Initialization
!!$    !
!!$    if ( .not. radiation_utils_inited ) call RadiationInit
!!$
!!$    ! Check arguments
!!$    !
!!$    if ( present( WNs ) ) then
!!$      if ( ( .not. present( WNe ) ) .or. ( .not. present( NumGaussNode ) ) ) then
!!$        call MessageNotify( 'E', module_name, &
!!$          & 'All of WNs, WNe, and NumGaussNode have to be present.' )
!!$      end if
!!$    else
!!$      if ( present( WNe ) .or. present( NumGaussNode ) ) then
!!$        call MessageNotify( 'E', module_name, &
!!$          & 'All of WNs, WNe, and NumGaussNode have to be present.' )
!!$      end if
!!$    end if
!!$
!!$    if ( present( WNs ) ) then
!!$      ! Case for non-grey atmosphere
!!$      !
!!$
!!$      ! Integrate Planck function and temperature derivative of it
!!$      !
!!$      call Integ_PF_GQ_Array3D( &
!!$        & WNs, WNe, NumGaussNode, &
!!$        & 0, imax-1, 1, jmax, 1, kmax, &
!!$        & xyz_Temp, &
!!$        & xyz_IntPF &
!!$        & )
!!$      call Integ_PF_GQ_Array2D( &
!!$        & WNs, WNe, NumGaussNode, &
!!$        & 0, imax-1, 1, jmax, &
!!$        & xy_SurfTemp, &
!!$        & xy_SurfIntPF &
!!$        & )
!!$      call Integ_DPFDT_GQ_Array3D( &
!!$        & 0, imax-1, 1, jmax, 1, kmax, &
!!$        & WNs, WNe, NumGaussNode, xyz_Temp, & ! (in )
!!$        & xyz_IntDPFDT          & ! (out)
!!$        & )
!!$      call Integ_DPFDT_GQ_Array2D( &
!!$        & 0, imax-1, 1, jmax, &
!!$        & WNs, WNe, NumGaussNode, xy_SurfTemp, & ! (in )
!!$        & xy_SurfIntDPFDT           & ! (out)
!!$        & )
!!$
!!$      xyz_IntPF       = PI * xyz_IntPF
!!$      xy_SurfIntPF    = PI * xy_SurfIntPF
!!$      xyz_IntDPFDT    = PI * xyz_IntDPFDT
!!$      xy_SurfIntDPFDT = PI * xy_SurfIntDPFDT
!!$
!!$    else
!!$
!!$      ! Case for grey atmosphere
!!$      !
!!$      xyz_IntPF       = StB * xyz_Temp**4
!!$      xy_SurfIntPF    = StB * xy_SurfTemp**4
!!$      xyz_IntDPFDT    = 4.0_DP * StB * xyz_Temp**3
!!$      xy_SurfIntDPFDT = 4.0_DP * StB * xy_SurfTemp**3
!!$
!!$    end if
!!$
!!$
!!$
!!$    ! 放射フラックス計算
!!$    ! Calculate radiation flux
!!$    !
!!$!    do k = 0, kmax
!!$!
!!$!      xyr_RadLFlux(:,:,k) = xy_SurfEmis * PI * xy_SurfIntPF * xyrr_Trans(:,:,k,0)
!!$!
!!$!      do kk = 0, kmax-1
!!$!        xyr_RadLFlux(:,:,k) = xyr_RadLFlux(:,:,k)                &
!!$!          & - PI * xyz_IntPF(:,:,kk+1)                           &
!!$!          & * ( xyrr_Trans(:,:,k,kk) - xyrr_Trans(:,:,k,kk+1) )
!!$!      end do
!!$!
!!$!    end do
!!$
!!$
!!$    !   Initialization
!!$    !
!!$    xyr_RadLDoFlux = 0.0_DP
!!$    xyr_RadLUpFlux = 0.0_DP
!!$    !
!!$    !   Downward flux
!!$    !
!!$    do k = kmax, 0, -1
!!$
!!$      do kk = kmax, k+1, -1
!!$        xyr_RadLDoFlux(:,:,k) = xyr_RadLDoFlux(:,:,k)          &
!!$          & + xyz_IntPF(:,:,kk)                                &
!!$          & * ( xyrr_Trans(:,:,k,kk-1) - xyrr_Trans(:,:,k,kk) )
!!$      end do
!!$
!!$    end do
!!$    !
!!$    !   Upward flux
!!$    !
!!$    !     Set upward flux
!!$    !
!!$    do k = 0, kmax
!!$
!!$      xyr_RadLUpFlux(:,:,k) = xy_SurfEmis * xy_SurfIntPF * xyrr_Trans(:,:,k,0)
!!$
!!$      do kk = 1, k
!!$        xyr_RadLUpFlux(:,:,k) = xyr_RadLUpFlux(:,:,k)          &
!!$          & - xyz_IntPF(:,:,kk)                                &
!!$          & * ( xyrr_Trans(:,:,k,kk-1) - xyrr_Trans(:,:,k,kk) )
!!$      end do
!!$
!!$    end do
!!$
!!$    xyr_RadLFlux = xyr_RadLUpFlux - xyr_RadLDoFlux
!!$
!!$
!!$    ! 放射フラックスの変化率の計算
!!$    ! Calculate rate of change of radiative flux
!!$    !
!!$    do k = 0, kmax
!!$      xyra_DelRadLFlux(:,:,k,0) =                           &
!!$        & xy_SurfEmis * xy_SurfIntDPFDT * xyrr_Trans(:,:,k,0)
!!$
!!$      xyra_DelRadLFlux(:,:,k,1) =                           &
!!$        & - xyz_IntDPFDT(:,:,1)                             &
!!$        &   * ( xyrr_Trans(:,:,k,0) - xyrr_Trans(:,:,k,1) )
!!$    end do
!!$
!!$
!!$  end subroutine RadiationRTEQNonScat
!!$
  !--------------------------------------------------------------------------------------

  subroutine RadRTEQNonScat(                                  &
    & xyz_IntPF, xy_SurfIntPF, xy_IntDPFDT1, xy_SurfIntDPFDT, & ! (in)
    & xyrr_Trans,                                             & ! (in)
    & xyr_RadLFlux, xyra_DelRadLFlux                          & ! (out)
    & )
    !
    ! 散乱なしの場合の放射伝達方程式の計算
    !
    ! Integrate radiative transfer equation without scattering
    !

    ! モジュール引用 ; USE statements
    !


    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in ) :: xyz_IntPF        (0:imax-1, 1:jmax, 1:kmax)
                              ! Integrated Planck function
    real(DP), intent(in ) :: xy_SurfIntPF     (0:imax-1, 1:jmax)
                              ! Integrated Planck function with surface temperature
    real(DP), intent(in ) :: xy_IntDPFDT1     (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
    real(DP), intent(in ) :: xy_SurfIntDPFDT  (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
                              ! with surface temperature
    real(DP), intent(in ) :: xyrr_Trans   (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
                              ! 透過係数. 
                              ! Transmission coefficient
    real(DP), intent(out) :: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux
    real(DP), intent(out) :: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with longwave


    ! 作業変数
    ! Work variables
    !
    real(DP):: xyr_RadLDoFlux (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_RadLUpFlux (0:imax-1, 1:jmax, 0:kmax)

    integer:: k, kk           ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction

    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. rad_utils_inited ) call RadInit



    ! 放射フラックス計算
    ! Calculate radiation flux
    !
!!$    do k = 0, kmax
!!$
!!$      xyr_RadLFlux(:,:,k) = xy_SurfEmis * PI * xy_SurfIntPF * xyrr_Trans(:,:,k,0)
!!$
!!$      do kk = 0, kmax-1
!!$        xyr_RadLFlux(:,:,k) = xyr_RadLFlux(:,:,k)                &
!!$          & - PI * xyz_IntPF(:,:,kk+1)                           &
!!$          & * ( xyrr_Trans(:,:,k,kk) - xyrr_Trans(:,:,k,kk+1) )
!!$      end do
!!$
!!$    end do


    !   Initialization
    !
    xyr_RadLDoFlux = 0.0_DP
    xyr_RadLUpFlux = 0.0_DP
    !
    !   Downward flux
    !
    do k = kmax, 0, -1

      do kk = kmax, k+1, -1
        xyr_RadLDoFlux(:,:,k) = xyr_RadLDoFlux(:,:,k)          &
          & + xyz_IntPF(:,:,kk)                                &
          & * ( xyrr_Trans(:,:,k,kk-1) - xyrr_Trans(:,:,k,kk) )
      end do

    end do
    !
    !   Upward flux
    !
    !     Set upward flux
    !
    do k = 0, kmax

      xyr_RadLUpFlux(:,:,k) = xy_SurfIntPF * xyrr_Trans(:,:,k,0)

      do kk = 1, k
        xyr_RadLUpFlux(:,:,k) = xyr_RadLUpFlux(:,:,k)          &
          & - xyz_IntPF(:,:,kk)                                &
          & * ( xyrr_Trans(:,:,k,kk-1) - xyrr_Trans(:,:,k,kk) )
      end do

    end do

    xyr_RadLFlux = xyr_RadLUpFlux - xyr_RadLDoFlux


    ! 放射フラックスの変化率の計算
    ! Calculate rate of change of radiative flux
    !
    do k = 0, kmax
      xyra_DelRadLFlux(:,:,k,0) =                           &
        & xy_SurfIntDPFDT * xyrr_Trans(:,:,k,0)

      xyra_DelRadLFlux(:,:,k,1) =                           &
        & - xy_IntDPFDT1(:,:)                               &
        &   * ( xyrr_Trans(:,:,k,0) - xyrr_Trans(:,:,k,1) )
    end do


  end subroutine RadRTEQNonScat

  !--------------------------------------------------------------------------------------

  subroutine RadRTEQNonScatWrapper(                   &
    & xyz_Temp, xy_SurfTemp, xy_SurfEmis, xyrr_Trans, & ! (in )
    & xyr_RadLFlux, xyra_DelRadLFlux,                 & ! (out)
    & WNs, WNe, NumGaussNode                          & ! (in ) optional
    & )
    !
    ! 散乱なしの場合の放射伝達方程式の計算
    !
    ! Integrate radiative transfer equation without scattering
    !

    ! モジュール引用 ; USE statements
    !


    ! プランク関数の計算
    ! Calculate Planck function
    !
    use planck_func, only :                            &
      & Integ_PF_GQ_Array3D   , Integ_PF_GQ_Array2D,   &
      & Integ_DPFDT_GQ_Array3D, Integ_DPFDT_GQ_Array2D

    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in ) :: xyz_Temp    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(in ) :: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度. 
                              ! Surface temperature
    real(DP), intent(in ) :: xy_SurfEmis (0:imax-1, 1:jmax)
                              ! 惑星表面射出率. 
                              ! Surface emissivity
    real(DP), intent(in ) :: xyrr_Trans   (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
                              ! 透過係数. 
                              ! Transmission coefficient
    real(DP), intent(out) :: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux
    real(DP), intent(out) :: xyra_DelRadLFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with longwave
    real(DP), intent(in ), optional :: WNs
    real(DP), intent(in ), optional :: WNe
    integer , intent(in ), optional :: NumGaussNode


    ! 作業変数
    ! Work variables
    !
    real(DP):: xyz_IntPF        (0:imax-1, 1:jmax, 1:kmax)
                              ! Integrated Planck function
    real(DP):: xy_SurfIntPF     (0:imax-1, 1:jmax)
                              ! Integrated Planck function with surface temperature
    real(DP):: xyz_IntDPFDT     (0:imax-1, 1:jmax, 1:kmax)
                              ! Integrated temperature derivative of Planck function
    real(DP):: xy_SurfIntDPFDT  (0:imax-1, 1:jmax)
                              ! Integrated temperature derivative of Planck function
                              ! with surface temperature


    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. rad_utils_inited ) call RadInit


    ! Check arguments
    !
    if ( present( WNs ) .or. present( WNe ) .or. present( NumGaussNode ) ) then
      if ( .not. ( present( WNs ) .and. present( WNe ) .and. present( NumGaussNode ) ) ) then
        call MessageNotify( 'E', module_name, &
          & 'All of WNs, WNe, and NumGaussNode have to be present.' )
      end if
    end if


    if ( present( WNs ) ) then
      ! Case for non-grey atmosphere
      !

      ! Integrate Planck function and temperature derivative of it
      !
      call Integ_PF_GQ_Array3D( &
        & WNs, WNe, NumGaussNode, &
        & 0, imax-1, 1, jmax, 1, kmax, &
        & xyz_Temp, &
        & xyz_IntPF &
        & )
      call Integ_PF_GQ_Array2D( &
        & WNs, WNe, NumGaussNode, &
        & 0, imax-1, 1, jmax, &
        & xy_SurfTemp, &
        & xy_SurfIntPF &
        & )
      call Integ_DPFDT_GQ_Array3D( &
        & 0, imax-1, 1, jmax, 1, kmax, &
        & WNs, WNe, NumGaussNode, xyz_Temp, & ! (in )
        & xyz_IntDPFDT          & ! (out)
        & )
      call Integ_DPFDT_GQ_Array2D( &
        & 0, imax-1, 1, jmax, &
        & WNs, WNe, NumGaussNode, xy_SurfTemp, & ! (in )
        & xy_SurfIntDPFDT           & ! (out)
        & )

      xyz_IntPF       = PI * xyz_IntPF
      xy_SurfIntPF    = PI * xy_SurfIntPF
      xyz_IntDPFDT    = PI * xyz_IntDPFDT
      xy_SurfIntDPFDT = PI * xy_SurfIntDPFDT

    else

      ! Case for grey atmosphere
      !
      xyz_IntPF       =                        StB * xyz_Temp**4
      xy_SurfIntPF    = xy_SurfEmis          * StB * xy_SurfTemp**4
      xyz_IntDPFDT    =               4.0_DP * StB * xyz_Temp**3
      xy_SurfIntDPFDT = xy_SurfEmis * 4.0_DP * StB * xy_SurfTemp**3

    end if


    call RadRTEQNonScat(                                               &
      & xyz_IntPF, xy_SurfIntPF, xyz_IntDPFDT(:,:,1), xy_SurfIntDPFDT, & ! (in)
      & xyrr_Trans,                                                    & ! (in)
      & xyr_RadLFlux, xyra_DelRadLFlux                                 & ! (out)
      & )


  end subroutine RadRTEQNonScatWrapper

  !--------------------------------------------------------------------------------------

  subroutine RadDTempDt( &
    & xyr_RadLFlux, xyr_RadSFlux, xyr_Press, & ! (in)
    & xyz_DTempDtRadL, xyz_DTempDtRadS &       ! (out)
    & )
    !
    ! 放射による温度変化率を計算します. 
    ! 
    ! Temperature tendency with radiation is calculated. 
    !

    ! モジュール引用 ; USE statements
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut


    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux
    real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux
    real(DP), intent(in):: xyr_Press    (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(out):: xyz_DTempDtRadL (0:imax-1, 1:jmax, 1:kmax)
                              ! 長波加熱率. 
                              ! Temperature tendency with longwave
    real(DP), intent(out):: xyz_DTempDtRadS (0:imax-1, 1:jmax, 1:kmax)
                              ! 短波加熱率. 
                              ! Temperature tendency with shortwave

    ! 作業変数
    ! Work variables
    !
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction

    ! 実行文 ; Executable statement
    !

    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )

    ! 初期化
    ! Initialization
    !
    if ( .not. rad_utils_inited ) call RadInit

    ! 放射冷却率の演算
    ! Calculate radiation cooling rate
    !
    do k = 1, kmax
      xyz_DTempDtRadL(:,:,k) = &
        & (     xyr_RadLFlux(:,:,k-1) - xyr_RadLFlux(:,:,k) ) &
        &   / ( xyr_Press(:,:,k-1)    - xyr_Press(:,:,k) ) &
        &   / CpDry * Grav

      xyz_DTempDtRadS(:,:,k) = &
        & (     xyr_RadSFlux(:,:,k-1) - xyr_RadSFlux(:,:,k) ) &
        &   / ( xyr_Press(:,:,k-1)    - xyr_Press(:,:,k) ) &
        &   / CpDry * Grav
    end do


    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'DTempDtRadL' , xyz_DTempDtRadL )
    call HistoryAutoPut( TimeN, 'DTempDtRadS' , xyz_DTempDtRadS )


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine RadDTempDt

  !--------------------------------------------------------------------------------------

  subroutine RadDTempDtforNHM(                        &
    & xyr_RadLFlux, xyr_RadSFlux, r_Height, xyz_Dens, & ! (in)
    & xyz_DTempDtRadL, xyz_DTempDtRadS                & ! (out)
    & )
    !
    ! 放射による温度変化率を計算します. 
    ! 
    ! Temperature tendency with radiation is calculated. 
    !

    ! モジュール引用 ; USE statements
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut


    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(in):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux
    real(DP), intent(in):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux
    real(DP), intent(in):: r_Height                       (0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: xyz_Dens     (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out):: xyz_DTempDtRadL (0:imax-1, 1:jmax, 1:kmax)
                              ! 長波加熱率. 
                              ! Temperature tendency with longwave
    real(DP), intent(out):: xyz_DTempDtRadS (0:imax-1, 1:jmax, 1:kmax)
                              ! 短波加熱率. 
                              ! Temperature tendency with shortwave

    ! 作業変数
    ! Work variables
    !
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction

    ! 実行文 ; Executable statement
    !

    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )

    ! 初期化
    ! Initialization
    !
    if ( .not. rad_utils_inited ) call RadInit

    ! 放射冷却率の演算
    ! Calculate radiation cooling rate
    !
    do k = 1, kmax
      xyz_DTempDtRadL(:,:,k) = &
        &   - ( xyr_RadLFlux(:,:,k) - xyr_RadLFlux(:,:,k-1) ) &
        &   / (   r_Height      (k) -   r_Height      (k-1) ) &
        &   / ( CpDry * xyz_Dens(:,:,k) )

      xyz_DTempDtRadS(:,:,k) = &
        &   - ( xyr_RadSFlux(:,:,k) - xyr_RadSFlux(:,:,k-1) ) &
        &   / (   r_Height      (k) -   r_Height      (k-1) ) &
        &   / ( CpDry * xyz_Dens(:,:,k) )
    end do


    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'DTempDtRadL' , xyz_DTempDtRadL )
    call HistoryAutoPut( TimeN, 'DTempDtRadS' , xyz_DTempDtRadS )


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine RadDTempDtforNHM

  !--------------------------------------------------------------------------------------

  subroutine RadDTempDtforNHM2DWrapper(             &
    & xr_RadLFlux, xr_RadSFlux, r_Height, xz_Dens,  & ! (in)
    & xz_DTempDtRadL, xz_DTempDtRadS                & ! (out)
    & )
    !
    ! 放射による温度変化率を計算します. 
    ! 
    ! Temperature tendency with radiation is calculated. 
    !

    ! モジュール引用 ; USE statements
    !

    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(in):: xr_RadLFlux (0:imax-1, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux
    real(DP), intent(in):: xr_RadSFlux (0:imax-1, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux
    real(DP), intent(in):: r_Height    (0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: xz_Dens     (0:imax-1, 1:kmax)
    real(DP), intent(out):: xz_DTempDtRadL (0:imax-1, 1:kmax)
                              ! 長波加熱率. 
                              ! Temperature tendency with longwave
    real(DP), intent(out):: xz_DTempDtRadS (0:imax-1, 1:kmax)
                              ! 短波加熱率. 
                              ! Temperature tendency with shortwave

    ! 作業変数
    ! Work variables
    !
    real(DP) :: xyr_RadLFlux   (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_RadSFlux   (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyz_Dens       (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DTempDtRadL(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DTempDtRadS(0:imax-1, 1:jmax, 1:kmax)


    ! 実行文 ; Executable statement
    !

    xyr_RadLFlux(:,1,:) = xr_RadLFlux(:,:)
    xyr_RadSFlux(:,1,:) = xr_RadSFlux(:,:)
    xyz_Dens    (:,1,:) = xz_Dens    (:,:)

    call RadDTempDtforNHM(                              &
      & xyr_RadLFlux, xyr_RadSFlux, r_Height, xyz_Dens, & ! (in)
      & xyz_DTempDtRadL, xyz_DTempDtRadS                & ! (out)
      & )

    xz_DTempDtRadL(:,:) = xyz_DTempDtRadL(:,1,:)
    xz_DTempDtRadS(:,:) = xyz_DTempDtRadS(:,1,:)


  end subroutine RadDTempDtforNHM2DWrapper

  !-------------------------------------------------------------------

  subroutine RadFluxOutput(                          &
    & xyr_RadSFlux, xyr_RadLFlux,                    & ! (in)
    & xyra_DelRadLFlux, xy_DSurfTempDt, xyz_DTempDt  & ! (in)
    & )
    !
    ! 放射フラックス (xyr_RadSFlux, xyr_RadLFlux) 
    ! について, その他の引数を用いて補正し, 出力を行う. 
    !
    ! Radiation fluxes (xyr_RadSFlux, xyr_RadLFlux) are
    ! corrected by using other arguments, and the corrected values are output.
    !

    ! モジュール引用 ; USE statements
    !

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & DelTime, &            ! $ \Delta t $ [s]
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(in):: xyr_RadSFlux    (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux
    real(DP), intent(in):: xyr_RadLFlux    (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux
    real(DP), intent(in):: xyra_DelRadLFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with longwave
    real(DP), intent(in):: xy_DSurfTempDt  (0:imax-1, 1:jmax)
                              ! 地表面温度変化率. 
                              ! Surface temperature tendency
    real(DP), intent(in):: xyz_DTempDt     (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \DP{T}{t} $ . 温度変化 (K s-1)
                              ! Temperature tendency (K s-1)

    ! 出力のための作業変数
    ! Work variables for output
    !
    real(DP):: xyr_RadLFluxCor (0:imax-1, 1:jmax, 0:kmax)
                              ! 補正された長波フラックス. 
                              ! Corrected longwave flux

    ! 作業変数
    ! Work variables
    !
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    ! 実行文 ; Executable statement
    !

    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )

    ! 初期化
    ! Initialization
    !
    if ( .not. rad_utils_inited ) call RadInit

    ! 長波フラックスの補正 ( 地表フラックス分の補正 )
    ! Correct longwave flux ( amount of surface flux )
    !
    ! Lines commented out below will be deleted soon (yot, 2010/10/31).
!!$    do k = 0, kmax
!!$      xyr_RadLFluxCor (:,:,k) = &
!!$        &     xyr_RadLFlux (:,:,k) &
!!$        &   + xyra_DelRadLFlux(:,:,k,0) * xy_DSurfTempDt (:,:) * DelTime
!!$    end do
    do k = 0, kmax
      xyr_RadLFluxCor(:,:,k) = &
        &   xyr_RadLFlux(:,:,k) &
        & + (   xy_DSurfTempDt     * xyra_DelRadLFlux(:,:,k,0)   &
        &     + xyz_DTempDt(:,:,1) * xyra_DelRadLFlux(:,:,k,1) ) &
        &   * DelTime
    end do

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'OLR' , xyr_RadLFluxCor(:,:,kmax) )
    call HistoryAutoPut( TimeN, 'SLR' , xyr_RadLFluxCor(:,:,0)    )
    call HistoryAutoPut( TimeN, 'OSR' , xyr_RadSFlux   (:,:,kmax) )
    call HistoryAutoPut( TimeN, 'SSR' , xyr_RadSFlux   (:,:,0)    )


    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'OLRB', xyr_RadLFlux   (:,:,kmax) )
    call HistoryAutoPut( TimeN, 'SLRB', xyr_RadLFlux   (:,:,0)    )
    call HistoryAutoPut( TimeN, 'OSRB', xyr_RadSFlux   (:,:,kmax) )
    call HistoryAutoPut( TimeN, 'SSRB', xyr_RadSFlux   (:,:,0)    )


    ! 長波フラックスの補正 ( 地表フラックス分の補正 )
    ! Correct longwave flux ( amount of surface flux )
    !
    ! Lines commented out below will be deleted soon (yot, 2010/10/31).
!!$    do k = 0, kmax
!!$      xyr_RadLFluxCor (:,:,k) = &
!!$        &     xyr_RadLFlux (:,:,k) &
!!$        &   + xyra_DelRadLFlux(:,:,k,0) * xy_DSurfTempDt (:,:) * 2.0d0 * DelTime
!!$    end do
    do k = 0, kmax
      xyr_RadLFluxCor(:,:,k) = &
        &   xyr_RadLFlux(:,:,k) &
        & + (   xy_DSurfTempDt     * xyra_DelRadLFlux(:,:,k,0)   &
        &     + xyz_DTempDt(:,:,1) * xyra_DelRadLFlux(:,:,k,1) ) &
        &   * 2.0_DP * DelTime
    end do

    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'OLRA', xyr_RadLFluxCor(:,:,kmax) )
    call HistoryAutoPut( TimeN, 'SLRA', xyr_RadLFluxCor(:,:,0)    )
    call HistoryAutoPut( TimeN, 'OSRA', xyr_RadSFlux   (:,:,kmax) )
    call HistoryAutoPut( TimeN, 'SSRA', xyr_RadSFlux   (:,:,0)    )


    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )

  end subroutine RadFluxOutput

  !-------------------------------------------------------------------

  subroutine RadInit
    !
    ! rad_utils モジュールの初期化を行います. 
    ! NAMELIST#rad_utils_nml の読み込みはこの手続きで行われます. 
    !
    ! "rad_utils" module is initialized. 
    ! "NAMELIST#rad_utils_nml" is loaded in this procedure. 
    !

    ! モジュール引用 ; USE statements
    !

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable

    ! 文字列操作
    ! Character handling
    !
    use dc_string, only: toChar

    ! 宣言文 ; Declaration statements
    !
    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                              ! IOSTAT of NAMELIST read


    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /rad_utils_nml/ &
      & DiffFact
!!$      & DelTimeLongValue, DelTimeLongUnit, &
!!$      & DelTimeShortValue, DelTimeShortUnit, &
!!$!
!!$      & LongBandNum, &
!!$      & LongAbsorpCoefQVap, LongAbsorpCoefDryAir, &
!!$      & LongBandWeight, LongPathLengthFact, &
!!$!
!!$      & ShortBandNum, &
!!$      & ShortAbsorpCoefQVap, ShortAbsorpCoefDryAir, &
!!$      & ShortBandWeight, ShortSecScat, &
!!$      & ShortAtmosAlbedo, &
!!$!
!!$      & RstInputFile, RstOutputFile
          !
          ! デフォルト値については初期化手続 "rad_utils#RadInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "rad_utils#RadInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( rad_utils_inited ) return

    ! デフォルト値の設定
    ! Default values settings
    !
    DiffFact = 1.66_DP

    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, &          ! (out)
        & namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml,                     & ! (in)
        & nml = rad_utils_nml,            & ! (out)
        & iostat = iostat_nml )             ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if


    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'OLR', &
      & (/ 'lon ', 'lat ', 'time' /), &
      & 'outgoing longwave', 'W m-2' )

    call HistoryAutoAddVariable( 'SLR', &
      & (/ 'lon ', 'lat ', 'time' /), &
      & 'surface longwave', 'W m-2' )

    call HistoryAutoAddVariable( 'OSR', &
      & (/ 'lon ', 'lat ', 'time' /), &
      & 'outgoing shortwave', 'W m-2' )

    call HistoryAutoAddVariable( 'SSR', &
      & (/ 'lon ', 'lat ', 'time' /), &
      & 'surface shortwave', 'W m-2' )

    call HistoryAutoAddVariable( 'OLRB', &
      & (/ 'lon ', 'lat ', 'time' /), &
      & 'outgoing longwave', 'W m-2' )

    call HistoryAutoAddVariable( 'SLRB', &
      & (/ 'lon ', 'lat ', 'time' /), &
      & 'surface longwave', 'W m-2' )

    call HistoryAutoAddVariable( 'OSRB', &
      & (/ 'lon ', 'lat ', 'time' /), &
      & 'outgoing shortwave', 'W m-2' )

    call HistoryAutoAddVariable( 'SSRB', &
      & (/ 'lon ', 'lat ', 'time' /), &
      & 'surface shortwave', 'W m-2' )

    call HistoryAutoAddVariable( 'OLRA', &
      & (/ 'lon ', 'lat ', 'time' /), &
      & 'outgoing longwave', 'W m-2' )

    call HistoryAutoAddVariable( 'SLRA', &
      & (/ 'lon ', 'lat ', 'time' /), &
      & 'surface longwave', 'W m-2' )

    call HistoryAutoAddVariable( 'OSRA', &
      & (/ 'lon ', 'lat ', 'time' /), &
      & 'outgoing shortwave', 'W m-2' )

    call HistoryAutoAddVariable( 'SSRA', &
      & (/ 'lon ', 'lat ', 'time' /), &
      & 'surface shortwave', 'W m-2' )

    call HistoryAutoAddVariable( 'DTempDtRadL', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'long wave radiative heating rate', 'K s-1' )

    call HistoryAutoAddVariable( 'DTempDtRadS', &
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &
      & 'short wave radiative heating rate', 'K s-1' )


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'DiffFact = %f', d = (/ DiffFact /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    rad_utils_inited = .true.

  end subroutine RadInit

  !-------------------------------------------------------------------

end module rad_utils
