!= 放射フラックス (バンドモデル)
!
!= Radiation flux (band model)
!
! Authors::   Yasuhiro MORIKAWA, Yukiko YAMADA
! Version::   $Id: radiation_band.f90,v 1.19 2009-02-16 16:46:15 morikawa Exp $ 
! Tag Name::  $Name: dcpam5-20090306 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module radiation_band
  !
  != 放射フラックス (バンドモデル)
  !
  != Radiation flux (band model)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 温度, 比湿, 気圧から, 放射フラックスを計算する
  ! バンドモデルです. 
  !
  ! This is a band model that calculates radiation flux from
  ! temperature, specific humidity, and air pressure.
  !
  !== Procedures List
  !
  ! RadiationFlux       :: 放射フラックスの計算
  ! RadiationDTempDt    :: 放射フラックスによる温度変化の計算
  ! RadiationFluxOutput :: 放射フラックスの出力
  ! ------------        :: ------------
  ! RadiationFlux       :: Calculate radiation flux
  ! RadiationDTempDt    :: Calculate temperature tendency with radiation flux
  ! RadiationFluxOutput :: Output radiation fluxes
  !
  !== NAMELIST
  !
  ! NAMELIST#radiation_band_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

  ! 日付および時刻の取り扱い
  ! Date and time handler
  !
  use dc_date_types, only: DC_DIFFTIME
                              ! 日時の差を表現するデータ型. 
                              ! Data type for difference about date and time
  ! 種別型パラメタ
  ! 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, RadiationDTempDt, RadiationFluxOutput

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


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

  ! 長波フラックス用情報
  ! Information for long wave flux
  !
  type(DC_DIFFTIME), save:: IntTimeLong
                              ! 長波フラックスを計算する時間間隔. 
                              ! Interval time of long wave flux calculation
  type(DC_DIFFTIME), save:: PrevTimeLong
                              ! 前回長波フラックスを計算した時刻
                              ! Time when long wave flux is calculated
  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
  !
  type(DC_DIFFTIME), save:: IntTimeShort
                              ! 短波フラックスを計算する時間間隔. 
                              ! Interval time of short wave flux calculation
  type(DC_DIFFTIME), save:: PrevTimeShort
                              ! 前回短波フラックスを計算した時刻
                              ! Time when short wave flux is calculated
  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 incoming
  !
  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_band'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name: dcpam5-20090306 $' // &
    & '$Id: radiation_band.f90,v 1.19 2009-02-16 16:46:15 morikawa Exp $'
                              ! モジュールのバージョン
                              ! Module version

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

  interface RadiationDTempDt
    module procedure RadiationDTempDt
  end interface

  interface RadiationFluxOutput
    module procedure RadiationFluxOutput
  end interface

contains

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

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

    ! 短波入射 (太陽入射)
    ! Short wave (insolation) incoming
    !
    use radiation_short_income, only: ShortIncoming
!    use radiation_short_income_sr, only: ShortIncoming

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

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

    ! 日付および時刻の取り扱い
    ! Date and time handler
    !
    use dc_date, only: operator(-), operator(>=), operator(+), toChar

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

    ! 宣言文 ; 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):: 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
    character(STRING):: str_debug
                              ! デバッグ用変数
                              ! Variable for debug

    ! 実行文 ; Executable statement
    !

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

    ! 初期化
    ! Initialization
    !
    if ( .not. radiation_band_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 ( TimeN - PrevTimeLong >= IntTimeLong .or. .not. Old_Flux_saved ) then

      if ( .not. Old_Flux_saved ) then
        PrevTimeLong = TimeN
      else
        PrevTimeLong = PrevTimeLong + IntTimeLong
      end if

      if ( Debug() ) then
        str_debug = toChar(TimeN)
        call DbgMessage( '%c: LongFlux is calculated at %c', &
          & c1 = module_name, c2 = trim(str_debug) )
      end if

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

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

      if ( Debug() ) then
        str_debug = toChar(TimeN)
        call DbgMessage( &
          & '%c: LongFlux is not calculated at %c. Save values are used.', &
          & c1 = module_name, c2 = trim(str_debug) )
      end if

      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 short wave (insolation)
    !
    if ( TimeN - PrevTimeShort >= IntTimeShort .or. .not. Old_Flux_saved ) then

      if ( .not. Old_Flux_saved ) then
        PrevTimeShort = TimeN
      else
        PrevTimeShort = PrevTimeShort + IntTimeShort
      end if

      if ( Debug() ) then
        str_debug = toChar(TimeN)
        call DbgMessage( &
          & '%c: ShortFlux is calculated at %c', &
          & c1 = module_name, c2 = trim(str_debug) )
      end if

      ! 短波入射の計算
      ! Calculate short wave (insolation) incoming radiation
      !
      call ShortIncoming( &
        & xy_IncomRadSFlux, xy_InAngle ) ! (out)

      ! 短波フラックスの計算
      ! 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

      if ( Debug() ) then
        str_debug = toChar(TimeN)
        call DbgMessage( &
          & '%c: ShortFlux is not calculated at %c. Save values are used.', &
          & c1 = module_name, c2 = trim(str_debug) )
      end if

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

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

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

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

  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 RadiationFluxOutput( &
    & xyr_RadSFlux, xyr_RadLFlux, &      ! (in)
    & xyra_DelRadLFlux, xy_DSurfTempDt & ! (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

    ! 出力のための作業変数
    ! 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. radiation_band_inited ) call RadiationInit

    ! 長波フラックスの補正 ( 地表フラックス分の補正 )
    ! Correct longwave flux ( amount of surface flux )
    !
    do k = 0, kmax
      xyr_RadLFluxCor (:,:,k) = &
        &     xyr_RadLFlux (:,:,k) &
        &   + xyra_DelRadLFlux(:,:,k,0) * xy_DSurfTempDt (:,:) * 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) )

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

  end subroutine RadiationFluxOutput

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

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

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & StartTime             ! 計算開始時刻. 

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

    ! 日付および時刻の取り扱い
    ! Date and time handler
    !
    use dc_date, only: DCDiffTimeCreate

    ! ファイル入出力補助
    ! 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 gtool_historyauto, only: HistoryAutoAddVariable

    ! 宣言文 ; 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_band_nml/ &
      & DelTimeLongValue, DelTimeLongUnit, &
      & DelTimeShortValue, DelTimeShortUnit, &
!
      & LongBandNum, &
      & LongAbsorpCoeffQVap, LongAbsorpCoeffDryAir, &
      & LongBandWeight, LongPathLengthFact, &
!
      & ShortBandNum, &
      & ShortAbsorpCoeffQVap, ShortAbsorpCoeffDryAir, &
      & ShortBandWeight, ShortSecScat
          !
          ! デフォルト値については初期化手続 "radiation_band#RadiationInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "radiation_band#RadiationInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    if ( radiation_band_inited ) return
    call InitCheck

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

    ! 長波フラックス用情報
    ! Information for long wave flux
    !
    PrevTimeLong = StartTime

    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
    !
    PrevTimeShort = StartTime

    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

    ! 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_band_nml, &  ! (out)
        & iostat = iostat_nml )        ! (out)
      close( unit_nml )

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

    ! 時間間隔の処理
    ! Handle interval time
    !
    call DCDiffTimeCreate( IntTimeLong, &     ! (out)
      & DelTimeLongValue, DelTimeLongUnit )   ! (in)
    call DCDiffTimeCreate( IntTimeShort, &    ! (out)
      & DelTimeShortValue, DelTimeShortUnit ) ! (in)

    ! バンド数, 吸収係数, バンドウェイトのチェック
    ! 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)

    ! 短波入射用変数の割付
    ! Allocate variables for short wave (insolation) incoming radiation
    !
    allocate( xy_IncomRadSFlux (0:imax-1, 1:jmax) )
    allocate( xy_InAngle (0:imax-1, 1:jmax) )

    ! 保存用の変数の割り付け
    ! 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 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' )

    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'DelTime:' )
    call MessageNotify( 'M', module_name, '  DelTimeLong  = %f [%c]', &
      & d = (/ DelTimeLongValue /), c1 = trim( DelTimeLongUnit ) )
    call MessageNotify( 'M', module_name, '  DelTimeShort = %f [%c]', &
      & d = (/ DelTimeShortValue /), 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, '-- version = %c', c1 = trim(version) )

    radiation_band_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_band
