!= 放射フラックス
!
!= Radiation flux
!
! Authors::   Yasuhiro MORIKAWA, Yukiko YAMADA
! Version::   $Id: $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module radiation
  !
  != 放射フラックス
  !
  != Radiation flux
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 温度, 比湿, 気圧から, 放射フラックスを計算します. 
  !
  ! Calculate radiation flux from temperature, specific humidity, and
  ! air pressure.
  !
  !== Procedures List
  !
  ! RadiationFlux      :: 放射フラックスの計算
  ! RadiationCorrect   :: 放射フラックス補正
  ! RadiationDTempDt   :: 放射フラックスによる温度変化の計算
  ! ------------  :: ------------
  ! RadiationFlux      :: Calculate radiation flux
  ! RadiationCorrect   :: Radiation flux correction
  ! RadiationDTempDt   :: Calculate temperature tendency with radiation flux
  !
  !== NAMELIST
  !
  ! NAMELIST#radiation_nml
  !

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

  ! 格子点設定
  ! 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

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

  ! NAMELIST ファイル入力に関するユーティリティ
  ! Utilities for NAMELIST file input
  !
  use namelist_util, only: MaxNmlArySize
                              ! NAMELIST から読み込む配列の最大サイズ. 
                              ! Maximum size of arrays loaded from NAMELIST

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

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public:: RadiationFlux, RadiationCorrect, RadiationDTempDt

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


  ! 非公開変数
  ! Private variables
  !
  logical, save:: Old_Flux_saved = .false.
                              ! 一度計算したフラックスを保存したことを示すフラグ. 
                              ! Flug for saving of flux calculated once

  ! 長波フラックス用情報
  ! Information for long wave flux
  !
  integer, save:: IntStepLong
                              ! 長波フラックスを計算するステップ間隔. 
                              ! Interval step of long wave flux calculation
  real(DP), save:: DelTimeLongValue
                              ! 長波フラックスを計算する時間間隔の数値. 
                              ! Value of interval step of long wave flux calculation
  character(STRING), save:: DelTimeLongUnit
                              ! 長波フラックスを計算する時間間隔の単位. 
                              ! Unit of interval step of long wave flux calculation

  integer, save:: LongBandNum
                              ! 長波バンド数. 
                              ! Number of long wave band
  real(DP), save:: LongAbsorpCoeffQVap (1:MaxNmlArySize)
                              ! $ k_R $ . 水の吸収係数. 
                              ! Absorption coefficient of water. 
  real(DP), save:: LongAbsorpCoeffDryAir (1:MaxNmlArySize)
                              ! $ \bar{k}_R $ . 空気の吸収係数. 
                              ! Absorption coefficient of air. 
  real(DP), save:: LongBandWeight (1:MaxNmlArySize)
                              ! バンドウェイト. 
                              ! Band weight. 
  real(DP), save:: LongPathLengthFact
                              ! 光路長のファクタ. 
                              ! Factor of optical length

  real(DP), allocatable, save:: xy_TempSave (:,:)
                              ! $ T $ .     温度 (保存用). Temperature (for save)
  real(DP), allocatable, save:: xyr_RadLFluxSave (:,:,:)
                              ! 長波フラックス (保存用). 
                              ! Long wave flux (for save)
  real(DP), allocatable, save:: xyra_DelRadLFluxSave (:,:,:,:)
                              ! 長波地表温度変化 (保存用). 
                              ! Surface temperature tendency with long wave (for save)

  ! 短波フラックス用情報
  ! Information for short wave flux
  !
  integer, save:: IntStepShort
                              ! 短波フラックスを計算するステップ間隔. 
                              ! Interval step of short wave flux calculation
  real(DP), save:: DelTimeShortValue
                              ! 短波 (日射) フラックスを計算する時間間隔の数値. 
                              ! Value of interval step of short wave (insolation) flux calculation
  character(STRING), save:: DelTimeShortUnit
                              ! 短波 (日射) フラックスを計算する時間間隔の単位. 
                              ! Unit of interval step of short wave (insolation) flux calculation

  integer, save:: ShortBandNum
                              ! 短波バンド数. 
                              ! Number of short wave band
  real(DP), save:: ShortAbsorpCoeffQVap (1:MaxNmlArySize)
                              ! $ k_S $ . 水の吸収係数. 
                              ! Absorption coefficient of water. 
  real(DP), save:: ShortAbsorpCoeffDryAir (1:MaxNmlArySize)
                              ! $ \bar{k}_S $ . 空気の吸収係数. 
                              ! Absorption coefficient of air. 
  real(DP), save:: ShortBandWeight (1:MaxNmlArySize)
                              ! バンドウェイト. 
                              ! Band weight. 
  real(DP), save:: ShortSecScat
                              ! 散乱の $ \sec \zeta $ .
                              ! $ \sec \zeta $ of scattering

  real(DP), allocatable, save:: xyr_RadSFluxSave (:,:,:)
                              ! 短波 (日射) フラックス (保存用). 
                              ! Short wave (insolation) flux (for save)

  ! 短波入射用情報
  ! Information for short wave incomming
  !
  real(DP), save:: SolarCoeff
                              ! 太陽定数. 
                              ! Solar constant
  real(DP), save:: AtmosAlbedo
                              ! 大気アルベド. 
                              ! Albedo of air

!!$  real(DP), save:: EpsOrb
!!$                              ! 赤道傾斜角. 
!!$                              ! Inclination of equator to orbit
!!$  real(DP), save:: EqnOrb
!!$                              ! 昇降点黄経. 
!!$                              ! 
  real(DP), save:: IncomAIns
                              ! $ A_{ins} $ . 年平均入射の係数. 
                              ! Coefficient of annual mean incoming radiation. 
  real(DP), save:: IncomBIns
                              ! $ B_{ins} $ . 年平均入射の係数. AIns に同じ. 
                              ! Coefficient of annual mean incoming radiation. 
                              ! Same as "AIns". 
  real(DP), save:: IncomAZet
                              ! $ A_{\zeta} $ . 年平均入射角の係数. AIns に同じ. 
                              ! Coefficient of annual mean incoming radiation. 
                              ! Same as "AIns". 
  real(DP), save:: IncomBZet
                              ! $ B_{\zeta} $ . 年平均入射角の係数. AIns に同じ. 
                              ! Coefficient of annual mean incoming radiation. 
                              ! Same as "AIns". 


  real(DP), allocatable, save:: xy_IncomRadSFlux (:,:)
                              ! 短波 (日射) フラックス. 
                              ! Short wave (insolation) flux
  real(DP), allocatable, save:: xy_InAngle (:,:)
                              ! sec (入射角). 
                              ! sec (angle of incidence)


  character(*), parameter:: module_name = 'radiation'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id:  $'
                              ! モジュールのバージョン
                              ! Module version

  ! INTERFACE 文 ; INTERFACE statements
  !
  interface RadiationFlux
    module procedure RadiationFlux
  end interface

  interface RadiationCorrect
    module procedure RadiationCorrect
  end interface

  interface RadiationDTempDt
    module procedure RadiationDTempDt
  end interface

contains

  subroutine RadiationFlux( &
    & xyz_Temp, xyz_QVap, xyr_Press, &         ! (in)
    & xy_SurfTemp, xy_SurfAlbedo, &            ! (in)
    & xyr_RadLFlux, xyr_RadSFlux, &            ! (out)
    & xya_SurfRadLMtx, xyra_DelRadLFlux &      ! (out)
    & )
    !
    ! 温度, 比湿, 気圧から, 放射フラックスを計算します. 
    !
    ! Calculate radiation flux from temperature, specific humidity, and 
    ! air pressure. 
    !

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration

    ! 時刻管理
    ! Time control
    !
    use timeset, only: Cstep, & ! 現在のステップ数. Current steps
      & TimesetClockStart, TimesetClockStop

    ! デバッグ用ユーティリティ
    ! Utilities for debug
    !
    use dc_trace, only: DbgMessage, BeginSub, EndSub

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: xyz_Temp  (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(in):: xyz_QVap  (0:imax-1, 1:jmax, 1:kmax)
                              ! $ q $ .     比湿. Specific humidity
    real(DP), intent(in):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \hat{p} $ . 気圧 (半整数レベル). 
                              ! Air pressure (half level)
    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 1:jmax)
                              ! 地表面温度. 
                              ! Surface temperature
    real(DP), intent(in):: xy_SurfAlbedo (0:imax-1, 1:jmax)
                              ! 地表アルベド. 
                              ! Surface albedo

    real(DP), intent(out):: xyr_RadLFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux
    real(DP), intent(out):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux
    real(DP), intent(out):: xya_SurfRadLMtx (0:imax-1, 1:jmax, -1:1)
                              ! $ T $ .  陰解行列: 地表. 
                              ! implicit matrix: surface
    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_TauQVap (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \int_z^{\infty} \rho q \, dz $ . 
                              ! 水に関する光学的厚さ 
                              ! $ \tau = \int_z^{\infty} k \rho q \, dz $
                              ! を吸収係数 $ k $ (高度 $ z $ に拠らず一定) 
                              ! で割ったもの.
                              ! 
                              ! Optical depth of water 
                              ! $ \tau = \int_z^{\infty} k \rho q \, dz $
                              ! divided by absorption coefficient $ k $ 
                              ! (this does not depend to height $ z $
                              ! and constant) of water
                              ! 
    real(DP):: xyr_TauDryAir (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \int_z^{\infty} \rho \, dz $ . 
                              ! 空気に関する光学的厚さ 
                              ! $ \tau = \int_z^{\infty} k \rho \, dz $
                              ! を吸収係数 $ k $ (高度 $ z $ に拠らず一定) 
                              ! で割ったもの.
                              ! 
                              ! Optical depth of air 
                              ! $ \tau = \int_z^{\infty} k \rho \, dz $
                              ! divided by absorption coefficient $ k $ 
                              ! (this does not depend to height $ z $
                              ! and constant) of air
                              ! 

    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. radiation_inited ) call RadiationInit

    ! 光学的厚さを吸収係数で割ったものの計算
    ! Calculate optical depth divided by absorption coefficient
    !
    xyr_TauQVap   = 0.0_DP
    xyr_TauDryAir = 0.0_DP

    do k = kmax-1, 0, -1
      xyr_TauQVap(:,:,k) = &
        &   xyr_TauQVap(:,:,k+1) &
        & + xyz_QVap(:,:,k+1) &
        &   * ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) / Grav

      xyr_TauDryAir(:,:,k) = &
        &   xyr_TauDryAir(:,:,k+1) &
        & + ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) / Grav
    end do

    ! 長波フラックスの算出
    ! Calculate long wave flux
    !
    if ( mod( Cstep, IntStepLong ) == 0 .or. .not. Old_Flux_saved ) then

      call DbgMessage( 'LongFlux is calculated' )

      call LongFlux( &
        & xyz_Temp, xy_SurfTemp, xyr_TauQVap, xyr_TauDryAir, & ! (in)
        & xyr_RadLFlux, xyra_DelRadLFlux )                     ! (out)

    ! 前回の値を利用
    ! Use values in last time
    !
    else

      call DbgMessage( &
        & 'LongFlux is not calculated. Save values are used' )

      xyr_RadLFlux = xyr_RadLFluxSave
      xyra_DelRadLFlux = xyra_DelRadLFluxSave
       
      do k = 0, kmax
        xyr_RadLFlux(:,:,k) = &
          &   xyr_RadLFlux(:,:,k) &
          & + xyra_DelRadLFlux(:,:,k,1) &
          &     * ( xyz_Temp(:,:,1) - xy_TempSave )

        xyra_DelRadLFlux(:,:,k,1) = &
          &   xyra_DelRadLFlux(:,:,k,1) &
          &     / ( xy_TempSave**3 ) * ( xyz_Temp(:,:,1)**3 )
      end do
    end if

    ! 長波陰解用行列の算出
    ! Calculate long wave implicit matrix
    !
    xya_SurfRadLMtx(:,:,0)  = xyra_DelRadLFlux(:,:,0,0)
    xya_SurfRadLMtx(:,:,1)  = xyra_DelRadLFlux(:,:,0,1)
    xya_SurfRadLMtx(:,:,-1) = 0.0_DP


    ! 短波 (日射) フラックスの算出
    ! Calculate short wave (insolation)
    !
    if ( mod( Cstep, IntStepShort ) == 0 .or. .not. Old_Flux_saved ) then

      call DbgMessage( 'ShortFlux is calculated' )

      ! 短波フラックスの計算
      ! Calculate short wave (insolation) flux
      !
      xyr_RadSFlux = 0.0_DP
      xyr_RadSFlux(:,:,kmax) = xy_IncomRadSFlux

      call ShortFlux( &
        & xyr_TauQVap, xyr_TauDryAir, xy_SurfAlbedo, & ! (in)
        & xyr_RadSFlux )                               ! (inout)

    ! 前回の値を利用
    ! Use values in last time
    ! 
    else

      call DbgMessage( &
        & 'ShortFlux is not calculated. Save values are used' )

      xyr_RadSFlux = xyr_RadSFluxSave
    end if


    ! 今回計算した値を保存
    ! Save calculated values in this time
    !
    xy_TempSave          = xyz_Temp (:,:,1)
    xyr_RadLFluxSave     = xyr_RadLFlux
    xyr_RadSFluxSave     = xyr_RadSFlux
    xyra_DelRadLFluxSave = xyra_DelRadLFlux

    Old_Flux_saved = .true.

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

  end subroutine RadiationFlux



  subroutine RadiationCorrect( &
    & xyz_DTempDt, xy_DSurfTempDt, xyra_DelRadLFlux, & ! (in)
    & xyr_RadLFlux &                                   ! (inout)
    & )
    !
    ! 放射フラックスの補正を行います. 
    !
    ! Radiation flux is corrected. 
    !

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime   ! $ \Delta t $ [s]

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

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

    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. radiation_inited ) call RadiationInit

    ! 放射フラックスの補正
    ! Correct radiation flux
    !
    do k = 0, kmax
      xyr_RadLFlux(:,:,k) = &
        &   xyr_RadLFlux(:,:,k) &
        & + (   xy_DSurfTempDt     * xyra_DelRadLFlux(:,:,k,0)   &
        &     + xyz_DTempDt(:,:,1) * xyra_DelRadLFlux(:,:,k,1) ) &
        &   * 2.0_DP * DelTime
    end do

  end subroutine RadiationCorrect



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

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, &  ! $ g $ [m s-2]. 
                                  ! 重力加速度. 
                                  ! Gravitational acceleration
      &                  CpDry    ! $ C_p $ [J kg-1 K-1]. 
                                  ! 乾燥大気の定圧比熱. 
                                  ! Specific heat of air at constant pressure

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

    ! 宣言文 ; 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
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. radiation_inited ) call RadiationInit

    ! 放射冷却率の演算
    ! 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 HistoryFileOutput( 'RadLFlux',    xyr_RadLFlux )
    call HistoryFileOutput( 'RadSFlux',    xyr_RadSFlux )
    call HistoryFileOutput( 'DTempDtRadL', xyz_DTempDtRadL )
    call HistoryFileOutput( 'DTempDtRadS', xyz_DTempDtRadS )

  end subroutine RadiationDTempDt



  subroutine LongFlux( &
    & xyz_Temp, xy_SurfTemp, xyr_TauQVap, xyr_TauDryAir, & ! (in)
    & xyr_RadLFlux, xyra_DelRadLFlux &                     ! (out)
    & )
    !
    ! 長波フラックスの計算
    !
    ! Calculate long wave flux
    !

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: StB  ! $ \sigma_{SB} $ . 
                              ! ステファンボルツマン定数. 
                              ! Stefan-Boltzmann constant

    ! 宣言文 ; Declaration statements
    !
    implicit none
    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):: xyr_TauQVap (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \int_z^{\infty} \rho q \, dz $ . 
                              ! 水に関する光学的厚さ 
                              ! $ \tau = \int_z^{\infty} k \rho q \, dz $
                              ! を吸収係数 $ k $ (高度 $ z $ に拠らず一定) 
                              ! で割ったもの.
                              ! 
                              ! Optical depth of water 
                              ! $ \tau = \int_z^{\infty} k \rho q \, dz $
                              ! divided by absorption coefficient $ k $ 
                              ! (this does not depend to height $ z $
                              ! and constant) of water
                              ! 
    real(DP), intent(in):: xyr_TauDryAir (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \int_z^{\infty} \rho \, dz $ . 
                              ! 空気に関する光学的厚さ 
                              ! $ \tau = \int_z^{\infty} k \rho \, dz $
                              ! を吸収係数 $ k $ (高度 $ z $ に拠らず一定) 
                              ! で割ったもの.
                              ! 
                              ! Optical depth of air 
                              ! $ \tau = \int_z^{\infty} k \rho \, dz $
                              ! divided by absorption coefficient $ k $ 
                              ! (this does not depend to height $ z $
                              ! and constant) of air
                              ! 
    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_Trans (0:imax-1, 1:jmax, 0:kmax)
                              ! 透過係数. 
                              ! Transmission coefficient
    real(DP):: xyr_Trans1 (0:imax-1, 1:jmax, 0:kmax)
                              ! 1/2 レベルからの透過係数. 
                              ! Transmission coefficient above 1/2 level
    real(DP):: xyr_Trans2 (0:imax-1, 1:jmax, 0:kmax)
                              ! 3/2 レベルからの透過係数. 
                              ! Transmission coefficient above 3/2 level
    real(DP):: xyz_PiB (0:imax-1, 1:jmax, 1:kmax)
                              ! $ \pi B = \sigma T^{4} $
    real(DP):: xy_SurfPiB (0:imax-1, 1:jmax)
                              ! 地表面の $ \pi B $ . 
                              ! $ \pi B $ on surface

    real(DP):: BandWeightSum  ! バンドウェイトの和
                              ! Sum of band weights

    integer:: k, kk           ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: bn              ! 波長について回る DO ループ用作業変数
                              ! Work variables for DO loop in wavenumber bands

    ! 実行文 ; Executable statement
    !

    ! バンドウェイトの設定
    ! Configure band weight
    !
    BandWeightSum = 0.0_DP

    do bn = 1, LongBandNum
      BandWeightSum =  BandWeightSum + LongBandWeight(bn)
    end do

    do bn = 1, LongBandNum
      LongBandWeight(bn) = LongBandWeight(bn) / BandWeightSum
    end do

    ! $ \pi B $ の計算
    ! Calculate $ \pi B $
    !
    xyz_PiB    = StB * ( xyz_Temp**4 )
    xy_SurfPiB = StB * ( xy_SurfTemp**4 )

    do k = 0, kmax

      ! 透過関数計算
      ! Calculate transmission functions
      !
      xyr_Trans = 0.0_DP
      
      do bn = 1, LongBandNum
        do kk = 0, kmax
          xyr_Trans(:,:,kk) = &
            &   xyr_Trans(:,:,kk) &
            & + LongBandWeight(bn) &
            &   * exp( - LongPathLengthFact &
            &            * (   LongAbsorpCoeffQVap(bn) &
            &                  * abs(   xyr_TauQVap(:,:,kk) &
            &                         - xyr_TauQVap(:,:,k)  ) &
            &                + LongAbsorpCoeffDryAir(bn) &
            &                  * abs(   xyr_TauDryAir(:,:,kk) &
            &                         - xyr_TauDryAir(:,:,k)  ) ) )
        end do
      end do

      ! 放射フラックス計算
      ! Calculate radiation flux
      !
      xyr_RadLFlux(:,:,k) = xy_SurfPiB * xyr_Trans(:,:,0)

      do kk = 0, kmax-1
        xyr_RadLFlux(:,:,k) = &
          &   xyr_RadLFlux(:,:,k) &
          & - xyz_PiB(:,:,kk+1) * ( xyr_Trans(:,:,kk) - xyr_Trans(:,:,kk+1) )
      end do

      ! 補正項計算用透過関数計算
      ! Calculate transmission functions for correction terms
      !
      xyr_Trans1(:,:,k) = xyr_Trans(:,:,0)
      xyr_Trans2(:,:,k) = xyr_Trans(:,:,1)

    end do

    ! 長波地表温度変化の計算
    ! Calclate surface temperature tendency with long wave
    !
    do k = 0, kmax
      xyra_DelRadLFlux(:,:,k,0) = &
        & 4.0_DP * xy_SurfPiB / xy_SurfTemp * xyr_Trans1(:,:,k)

      xyra_DelRadLFlux(:,:,k,1) = &
        & 4.0_DP * xyz_PiB(:,:,1) / xyz_Temp(:,:,1) &
        &   * ( xyr_Trans2(:,:,k) - xyr_Trans1(:,:,k) )
    end do

  end subroutine LongFlux



  subroutine ShortFlux( &
    & xyr_TauQVap, xyr_TauDryAir, xy_SurfAlbedo, & ! (in)
    & xyr_RadSFlux &                               ! (inout)
    & )
    !
    ! 短波フラックスを計算します.
    !
    ! Calculate short wave flux. 
    !

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

    ! 宣言文 ; Declaration statements
    !
    implicit none

    real(DP), intent(in):: xyr_TauQVap (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \int_z^{\infty} \rho q \, dz $ . 
                              ! 水に関する光学的厚さ 
                              ! $ \tau = \int_z^{\infty} k \rho q \, dz $
                              ! を吸収係数 $ k $ (高度 $ z $ に拠らず一定) 
                              ! で割ったもの.
                              ! 
                              ! Optical depth of water 
                              ! $ \tau = \int_z^{\infty} k \rho q \, dz $
                              ! divided by absorption coefficient $ k $ 
                              ! (this does not depend to height $ z $
                              ! and constant) of water
                              ! 
    real(DP), intent(in):: xyr_TauDryAir (0:imax-1, 1:jmax, 0:kmax)
                              ! $ \int_z^{\infty} \rho \, dz $ . 
                              ! 空気に関する光学的厚さ 
                              ! $ \tau = \int_z^{\infty} k \rho \, dz $
                              ! を吸収係数 $ k $ (高度 $ z $ に拠らず一定) 
                              ! で割ったもの.
                              ! 
                              ! Optical depth of air 
                              ! $ \tau = \int_z^{\infty} k \rho \, dz $
                              ! divided by absorption coefficient $ k $ 
                              ! (this does not depend to height $ z $
                              ! and constant) of air
                              ! 
    real(DP), intent(in):: xy_SurfAlbedo (0:imax-1, 1:jmax)
                              ! 地表アルベド. 
                              ! Surface albedo
    real(DP), intent(inout):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux

    ! 作業変数
    ! Work variables
    !
    real(DP):: BandWeightSum  ! バンドウェイトの和
                              ! Sum of band weights
    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
    integer:: bn              ! 波長について回る DO ループ用作業変数
                              ! Work variables for DO loop in wavenumber bands

    ! 実行文 ; Executable statement
    !

    ! バンドウェイトの設定
    ! Configure band weight
    !
    BandWeightSum = 0.0_DP

    do bn = 1, ShortBandNum
      BandWeightSum = BandWeightSum + ShortBandWeight(bn)
    end do

    do bn = 1, ShortBandNum
      ShortBandWeight(bn) = ShortBandWeight(bn) / BandWeightSum
    end do

    do bn = 1, ShortBandNum
      do k = 0, kmax

        ! 各レベルでの下向き透過
        ! Downward transmission on each level
        !
        if ( k /= kmax ) then
          xyr_RadSFlux(:,:,k) = &
            &   xyr_RadSFlux(:,:,k) & 
            & + ShortBandWeight(bn) * xyr_RadSFlux(:,:,kmax) &
            &   * exp( - xy_InAngle(:,:) &
            &            * (   ShortAbsorpCoeffQVap(bn) * xyr_TauQVap(:,:,k) &
            &                + ShortAbsorpCoeffDryAir(bn) * xyr_TauDryAir(:,:,k) ) &
            &       )
        end if
        
        ! 各レベルでの上向き透過
        ! Upward transmission on each level
        !
        xyr_RadSFlux(:,:,k) = &
          &   xyr_RadSFlux(:,:,k) & 
          & - ShortBandWeight(bn) * xyr_RadSFlux(:,:,kmax) &
          &   * exp( - xy_InAngle(:,:) &
          &            * (   ShortAbsorpCoeffQVap(bn) * xyr_TauQVap(:,:,0) &
          &                + ShortAbsorpCoeffDryAir(bn) * xyr_TauDryAir(:,:,0) ) &
          &       ) &
          &   * xy_SurfAlbedo &
          &   * exp( - ShortSecScat &
          &            * (   ShortAbsorpCoeffQVap(bn) &
          &                  * ( xyr_TauQVap(:,:,0) - xyr_TauQVap(:,:,k) ) &
          &                + ShortAbsorpCoeffDryAir(bn) &
          &                  * ( xyr_TauDryAir(:,:,0) - xyr_TauDryAir(:,:,k) ) &
          &              ) &
          &        )
      end do
    end do

    ! 吸収なしの場合
    ! In the case of no absorption
    !
    if ( ShortBandNum == 0 ) then
      do k = 0, kmax
        xyr_RadSFlux(:,:,k) = &
          & ( 1.0_DP - xy_SurfAlbedo ) * xyr_RadSFlux(:,:,kmax)
      end do
    end if

  end subroutine ShortFlux



  subroutine ShortIncoming( &
    & xy_IncomRadSFlux, xy_InAngle & ! (out)
    & )
    !
    ! 短波入射を計算します.
    !
    ! Calculate short wave (insolation) incoming radiation. 
    !

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: y_Lat  ! $ \varphi $ [rad.] . 緯度. Latitude

    ! 宣言文 ; Declaration statements
    !
    implicit none
    real(DP), intent(out):: xy_IncomRadSFlux (0:imax-1, 1:jmax)
                              ! 短波 (日射) フラックス. 
                              ! Short wave (insolation) flux
    real(DP), intent(out):: xy_InAngle (0:imax-1, 1:jmax)
                              ! sec (入射角). 
                              ! sec (angle of incidence)

    ! 作業変数
    ! Work variables
    !
    integer:: i               ! 経度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in longitude
    integer:: j               ! 緯度方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in latitude

    ! 実行文 ; Executable statement
    !

    ! 年, 日平均日射の計算
    ! Calculate annual mean, daily mean insolation
    !
    do i = 0, imax - 1
      do j = 1, jmax
        xy_IncomRadSFlux(i,j) = &
          &  - SolarCoeff * (1.0_DP - AtmosAlbedo ) &
          &    * ( IncomAIns + IncomBIns * cos( y_Lat(j) )**2 )

        if ( xy_IncomRadSFlux(i,j) < 0.0_DP ) then
          xy_InAngle(i,j) = 1.0_DP / ( IncomAZet + IncomBZet * cos( y_Lat(j) )**2 )
        else
          xy_IncomRadSFlux(i,j) = 0.0_DP
          xy_InAngle(i,j) = 0.0_DP
        end if
      end do
    end do


    ! 季節変化, 日変化がある場合の計算
    ! Calculate with seasonal change and diurnal change
    !

    ! AGCM5 から移植予定. 
    ! Importation from AGCM5 is planed

  end subroutine ShortIncoming



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

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimesetIntStepEval

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

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

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output

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

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

    ! 宣言文 ; Declaration statements
    !
    implicit none

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

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /radiation_nml/ &
      & DelTimeLongValue, DelTimeLongUnit, &
      & DelTimeShortValue, DelTimeShortUnit, &
!
      & LongBandNum, &
      & LongAbsorpCoeffQVap, LongAbsorpCoeffDryAir, &
      & LongBandWeight, LongPathLengthFact, &
!
      & ShortBandNum, &
      & ShortAbsorpCoeffQVap, ShortAbsorpCoeffDryAir, &
      & ShortBandWeight, ShortSecScat, &
!
      & SolarCoeff, AtmosAlbedo, &
!!!$      & EpsOrb, EqnOrb, &
      & IncomAIns, IncomBIns, IncomAZet, IncomBZet
          !
          ! デフォルト値については初期化手続 "radiation#RadiationInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "radiation#RadiationInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( radiation_inited ) return
    call InitCheck

    ! デフォルト値の設定
    ! Default values settings
    !

    ! 長波フラックス用情報
    ! Information for long wave flux
    !
    IntStepLong      = 1
    DelTimeLongValue = 3.0_DP
    DelTimeLongUnit  = 'hrs.'

    LongBandNum      = 4
    LongAbsorpCoeffQVap    = -999.9_DP
    LongAbsorpCoeffDryAir  = -999.9_DP
    LongBandWeight         = -999.9_DP
    LongAbsorpCoeffQVap   (1:LongBandNum) = (/ 8.0_DP, 1.0_DP, 0.1_DP, 0.0_DP /)
    LongAbsorpCoeffDryAir (1:LongBandNum) = (/ 0.0_DP, 0.0_DP, 0.0_DP, 5.0e-5_DP /)
    LongBandWeight        (1:LongBandNum) = (/ 0.2_DP, 0.1_DP, 0.1_DP, 0.6_DP /)
    LongPathLengthFact = 1.5_DP

    ! 短波フラックス用情報
    ! Information for short wave flux
    !
    IntStepShort      = 1
    DelTimeShortValue = 1.0_DP
    DelTimeShortUnit  = 'hrs.'

    ShortBandNum = 1
    ShortAbsorpCoeffQVap   = -999.9_DP
    ShortAbsorpCoeffDryAir = -999.9_DP
    ShortBandWeight        = -999.9_DP
    ShortAbsorpCoeffQVap   (1:ShortBandNum) = (/ 0.002_DP /)
    ShortAbsorpCoeffDryAir (1:ShortBandNum) = (/ 0.0_DP /)
    ShortBandWeight        (1:ShortBandNum) = (/ 1.0_DP /)
    ShortSecScat = 1.66_DP

    ! 短波入射用情報
    ! Information for short wave incomming
    !
    SolarCoeff  = 1380.0_DP
    AtmosAlbedo = 0.2_DP
!!$    EpsOrb = 23.5_DP
!!$    EqnOrb = 110.0_DP
    IncomAIns   = 0.127_DP
    IncomBIns   = 0.183_DP
    IncomAZet   = 0.410_DP
    IncomBZet   = 0.590_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 = radiation_nml, &  ! (out)
        & iostat = iostat_nml )   ! (out)
      close( unit_nml )

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

    ! ステップ間隔の処理
    ! Handle interval step
    !
    call TimesetIntStepEval( &
      & DelTimeLongValue, DelTimeLongUnit, & ! (in)
      & IntStepLong )                        ! (out)

    call TimesetIntStepEval( &
      & DelTimeShortValue, DelTimeShortUnit, & ! (in)
      & IntStepShort )                         ! (out)

    ! バンド数, 吸収係数, バンドウェイトのチェック
    ! Check number of band, absorption coefficients, band weight
    !
    call NmlutilAryValid( module_name, &              ! (in)
      & LongAbsorpCoeffQVap, 'LongAbsorpCoeffQVap', & ! (in)
      & LongBandNum,         'LongBandNum' )          ! (in)

    call NmlutilAryValid( module_name, &                  ! (in)
      & LongAbsorpCoeffDryAir, 'LongAbsorpCoeffDryAir', & ! (in)
      & LongBandNum,           'LongBandNum' )            ! (in)

    call NmlutilAryValid( module_name, &     ! (in)
      & LongBandWeight, 'LongBandWeight', &  ! (in)
      & LongBandNum,    'LongBandNum' )      ! (in)

    call NmlutilAryValid( module_name, &                ! (in)
      & ShortAbsorpCoeffQVap, 'ShortAbsorpCoeffQVap', & ! (in)
      & ShortBandNum,         'ShortBandNum' )          ! (in)

    call NmlutilAryValid( module_name, &                    ! (in)
      & ShortAbsorpCoeffDryAir, 'ShortAbsorpCoeffDryAir', & ! (in)
      & ShortBandNum,           'ShortBandNum' )            ! (in)

    call NmlutilAryValid( module_name, &       ! (in)
      & ShortBandWeight, 'ShortBandWeight', &  ! (in)
      & ShortBandNum,    'ShortBandNum' )      ! (in)

    ! 短波入射の計算
    ! Calculate short wave (insolation) incoming radiation
    !
    allocate( xy_IncomRadSFlux (0:imax-1, 1:jmax) )
    allocate( xy_InAngle (0:imax-1, 1:jmax) )

    call ShortIncoming( &
      & xy_IncomRadSFlux, xy_InAngle ) ! (out)

    ! 保存用の変数の割り付け
    ! Allocate variables for saving
    !
    allocate( xy_TempSave          (0:imax-1, 1:jmax) )
    allocate( xyr_RadLFluxSave     (0:imax-1, 1:jmax, 0:kmax) )
    allocate( xyr_RadSFluxSave     (0:imax-1, 1:jmax, 0:kmax) )
    allocate( xyra_DelRadLFluxSave (0:imax-1, 1:jmax, 0:kmax, 0:1) )

    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryFileRegister( 'RadLFlux' ,   'longwave flux', &
      & 'W m-2', StoA('lon', 'lat', 'sigm', 'time') )
    call HistoryFileRegister( 'RadSFlux' ,    'shortwave flux', &
      & 'W m-2', StoA('lon', 'lat', 'sigm', 'time') )
    call HistoryFileRegister( 'DTempDtRadL' , 'longwave heating', &
      & 'K s-1', StoA('lon', 'lat', 'sig', 'time') )
    call HistoryFileRegister( 'DTempDtRadS' , 'shortwave heating', &
      & 'K s-1', StoA('lon', 'lat', 'sig', 'time') )

    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'IntStep:' )
    call MessageNotify( 'M', module_name, '  IntStepLong  = %d', i = (/ IntStepLong  /) )
    call MessageNotify( 'M', module_name, '  IntStepShort = %d', i = (/ IntStepShort  /) )
    call MessageNotify( 'M', module_name, 'DelTime:' )
    call MessageNotify( 'M', module_name, '  DelTimeLongValue  = %f', d = (/ DelTimeLongValue  /) )
    call MessageNotify( 'M', module_name, '  DelTimeLongUnit   = %c', c1 = trim( DelTimeLongUnit   ) )
    call MessageNotify( 'M', module_name, '  DelTimeShortValue = %f', d = (/ DelTimeShortValue /) )
    call MessageNotify( 'M', module_name, '  DelTimeShortUnit  = %c', c1 = trim( DelTimeShortUnit  ) )
!
    call MessageNotify( 'M', module_name, 'LongFlux:' )
    call MessageNotify( 'M', module_name, '  LongBandNum            = %d', i = (/ LongBandNum /) )
    call MessageNotify( 'M', module_name, '  LongAbsorpCoeffQVap    = (/ %*r /)', &
      & r = real( LongAbsorpCoeffQVap(1:LongBandNum) ), n = (/ LongBandNum /) )
    call MessageNotify( 'M', module_name, '  LongAbsorpCoeffDryAir  = (/ %*r /)', &
      & r = real( LongAbsorpCoeffDryAir(1:LongBandNum) ), n = (/ LongBandNum /) )
    call MessageNotify( 'M', module_name, '  LongBandWeight         = (/ %*r /)', &
      & r = real( LongBandWeight(1:LongBandNum) ), n = (/ LongBandNum /) )
    call MessageNotify( 'M', module_name, '  LongPathLengthFact     = %f', d = (/ LongPathLengthFact /) )
!
    call MessageNotify( 'M', module_name, 'ShortFlux:' )
    call MessageNotify( 'M', module_name, '  ShortBandNum           = %d', i = (/ ShortBandNum /) )
    call MessageNotify( 'M', module_name, '  ShortAbsorpCoeffQVap   = (/ %*r /)', &
      & r = real( ShortAbsorpCoeffQVap(1:ShortBandNum) ), n = (/ ShortBandNum /) )
    call MessageNotify( 'M', module_name, '  ShortAbsorpCoeffDryAir = (/ %*r /)', &
      & r = real( ShortAbsorpCoeffDryAir(1:ShortBandNum) ), n = (/ ShortBandNum /) )
    call MessageNotify( 'M', module_name, '  ShortBandWeight        = (/ %*r /)', &
      & r = real( ShortBandWeight(1:ShortBandNum) ), n = (/ ShortBandNum /) )
    call MessageNotify( 'M', module_name, '  ShortSecScat           = %f', d = (/ ShortSecScat /) )
!
    call MessageNotify( 'M', module_name, 'ShortIncomming:' )
    call MessageNotify( 'M', module_name, '  SolarCoeff  = %f', d = (/ SolarCoeff  /) )
    call MessageNotify( 'M', module_name, '  AtmosAlbedo = %f', d = (/ AtmosAlbedo /) )
!!$    call MessageNotify( 'M', module_name, '  EpsOrb      = %f', d = (/ EpsOrb      /) )
!!$    call MessageNotify( 'M', module_name, '  EqnOrb      = %f', d = (/ EqnOrb      /) )
    call MessageNotify( 'M', module_name, '  IncomAIns   = %f', d = (/ IncomAIns   /) )
    call MessageNotify( 'M', module_name, '  IncomBIns   = %f', d = (/ IncomBIns   /) )
    call MessageNotify( 'M', module_name, '  IncomAZet   = %f', d = (/ IncomAZet   /) )
    call MessageNotify( 'M', module_name, '  IncomBZet   = %f', d = (/ IncomBZet   /) )

    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    radiation_inited = .true.
  end subroutine RadiationInit



  subroutine InitCheck
    !
    ! 依存モジュールの初期化チェック
    !
    ! Check initialization of dependency modules

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

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

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: gridset_inited

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: constants_inited

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: axesset_inited

    ! 時刻管理
    ! Time control
    !
    use timeset, only: timeset_inited

    ! 実行文 ; Executable statement
    !

    if ( .not. namelist_util_inited ) &
      & call MessageNotify( 'E', module_name, '"namelist_util" module is not initialized.' )

    if ( .not. gridset_inited ) &
      & call MessageNotify( 'E', module_name, '"gridset" module is not initialized.' )

    if ( .not. constants_inited ) &
      & call MessageNotify( 'E', module_name, '"constants" module is not initialized.' )

    if ( .not. axesset_inited ) &
      & call MessageNotify( 'E', module_name, '"axesset" module is not initialized.' )

    if ( .not. timeset_inited ) &
      & call MessageNotify( 'E', module_name, '"timeset" module is not initialized.' )

  end subroutine InitCheck

end module radiation