!= 放射フラックス
!
!= 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


  ! モジュール引用 ; 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

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

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

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


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

  !  長波フラックス用情報
  !  Information for long wave flux
  !
  real(DP), save:: DeltaTimeLongValue = 3.0_DP
                              ! 長波フラックスを計算するステップ間隔の数値. 
                              ! Value of step interval of long wave flux calculation
  character(STRING), save:: DeltaTimeLongUnit = 'hrs.'
                              ! 長波フラックスを計算するステップ間隔の単位. 
                              ! Unit of step interval of long wave flux calculation

  integer, save:: LongBandNum = 4
                              ! 長波バンド数. 
                              ! Number of long wave band
  real(DP), save:: LongAbsorpCoeffQVap(1:MaxNmlArySize)
                              ! 水の吸収係数. デフォルト値は
                              ! (/ 8.0, 1.0, 0.1, 0.0 /). 
                              !
                              ! Absorption coefficient of water. 
                              ! Default value is 
                              ! (/ 8.0, 1.0, 0.1, 0.0 /). 
  real(DP), save:: LongAbsorpCoeffDryAir(1:MaxNmlArySize)
                              ! 空気の吸収係数. デフォルト値は
                              ! (/ 0.0, 0.0, 0.0, 5.0e-5 /)
                              !
                              ! Absorption coefficient of air. 
                              ! Default value is 
                              ! (/ 0.0, 0.0, 0.0, 5.0e-5 /)
  real(DP), save:: LongBandWeight(1:MaxNmlArySize)
                              ! バンドウェイト. デフォルト値は
                              ! (/ 0.2, 0.1, 0.1, 0.6 /)
                              !
                              ! Band weight. 
                              ! Default value is 
                              ! (/ 0.2, 0.1, 0.1, 0.6 /)
  real(DP), save:: LongPathLengthFact = 1.5_DP
                              ! 光路長のファクタ. 
                              ! 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
  !
  real(DP), save:: DeltaTimeShortValue = 1.0_DP
                              ! 短波 (日射) フラックスを計算するステップ間隔の数値. 
                              ! Value of step interval of short wave (insolation) flux calculation
  character(STRING), save:: DeltaTimeShortUnit = 'hrs.'
                              ! 短波 (日射) フラックスを計算するステップ間隔の単位. 
                              ! Unit of step interval of short wave (insolation) flux calculation

  integer, save:: ShortBandNum = 1
                              ! 短波バンド数. 
                              ! Number of short wave band
  real(DP), save:: ShortAbsorpCoeffQVap(1:MaxNmlArySize)
                              ! 水の吸収係数. デフォルト値は
                              ! (/ 0.002 /)
                              !
                              ! Absorption coefficient of water. 
                              ! Default value is 
                              ! (/ 0.002 /)
  real(DP), save:: ShortAbsorpCoeffDryAir(1:MaxNmlArySize)
                              ! 空気の吸収係数. デフォルト値は
                              ! (/ 0.0 /)
                              !
                              ! Absorption coefficient of air. 
                              ! Default value is 
                              ! (/ 0.0 /)
  real(DP), save:: ShortBandWeight(1:MaxNmlArySize)
                              ! バンドウェイト. デフォルト値は
                              ! (/ 1.0 /)
                              !
                              ! Band weight. 
                              ! Default value is 
                              ! (/ 1.0 /)
  real(DP), save:: ShortSecScat = 1.66_DP
                              ! 散乱の $ 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 = 1380.0_DP
                              ! 太陽定数. 
                              ! Solar constant
  real(DP), save:: AtmosAlbedo = 0.2_DP
                              ! 大気アルベド. 
                              ! Albedo of air

!!$  real(DP), save:: EpsOrb = 23.5_DP
!!$                              ! 赤道傾斜角. 
!!$                              ! Inclination of equator to orbit
!!$  real(DP), save:: EqnOrb = 110.0_DP
!!$                              ! 昇降点黄経. 
!!$                              ! 
  real(DP), save:: IncomAIns   = 0.127_DP
                              ! 年平均入射の係数. 
                              ! Coefficient of annual mean incoming radiation. 
  real(DP), save:: IncomBIns   = 0.183_DP
                              ! 年平均入射の係数. AIns に同じ. 
                              ! Coefficient of annual mean incoming radiation. 
                              ! Same as "AIns". 
  real(DP), save:: IncomAZet   = 0.410_DP
                              ! 年平均入射角の係数. AIns に同じ. 
                              ! Coefficient of annual mean incoming radiation. 
                              ! Same as "AIns". 
  real(DP), save:: IncomBZet   = 0.590_DP
                              ! 年平均入射角の係数. 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)

  !  NAMELIST 変数群
  !  NAMELIST group name
  !
  namelist /radiation_nml/ &
    & DeltaTimeLongValue, DeltaTimeLongUnit, &
    & DeltaTimeShortValue, DeltaTimeShortUnit, &
!
    & LongBandNum, &
    & LongAbsorpCoeffQVap, LongAbsorpCoeffDryAir, &
    & LongBandWeight, LongPathLengthFact, &
!
    & ShortBandNum, &
    & ShortAbsorpCoeffQVap, ShortAbsorpCoeffDryAir, &
    & ShortBandWeight, ShortSecScat, &
!
    & SolarCoeff, AtmosAlbedo, &
!!$    & EpsOrb, EqnOrb, &
    & IncomAIns, IncomBIns, IncomAZet, IncomBZet


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

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

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


    character(*), parameter:: subname = 'Radiation'

    ! 実行文 ; Executable statement

    if ( .not. initialized ) call RadiationInit

    ! 光学的厚さの計算
    ! Calculate optical depth
    !
    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

    ! current_time / delta_time_RadLong == 0, または
    ! delta_time_RadLong 未設定, または
    ! まだ一度も計算していない場合には計算を行う. 
    !
    ! If current_time / delta_time_RadLong == 0, or
    ! "delta_time_RadLong" is not configured, or
    ! calculation was never performed, 
    ! calculation is performed. 
    !
    if ( mod( phy_radflx % current_time, &
      &       phy_radflx % delta_time_RadLong) == 0 .or. &
      &  EvalSec( phy_radflx % delta_time_RadLong ) < 0.0_DP .or. &
      &  .not. save_flag                                            ) 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)

    ! current_time / delta_time_RadShort == 0, または
    ! delta_time_RadShort 未設定, または
    ! まだ一度も計算していない場合には計算を行う. 
    !
    ! If current_time / delta_time_RadShort == 0, or
    ! "delta_time_RadShort" is not configured, or
    ! calculation was never performed, 
    ! calculation is performed. 
    ! 
    if ( mod( phy_radflx % current_time, &
      &       phy_radflx % delta_time_RadShort) == 0 .or. &
      &  EvalSec( phy_radflx % delta_time_RadShort ) < 0.0_DP .or. &
      &  .not. save_flag                                             ) then

      call DbgMessage( 'ShortFlux is calculated' )

      ! 短波フラックスの計算
      ! Calculate short wave (insolation) flux
      !
      call ShortFlux( &
        & xyr_TauQVap, xyr_TauDryAir, & ! (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

    save_flag = .true.

  end subroutine RadiationFlux

  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

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

    ! 宣言文 ; 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)
                              ! 光学的厚さ：水
                              ! Optical depth of water
    real(DP), intent(in):: xyr_TauDryAir (0:imax-1, 1:jmax, 0:kmax)
                              ! 光学的厚さ：空気
                              ! Optical depth 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 , bn      ! DO ループ用作業変数
                              ! Work variables for DO loop

    ! 実行文 ; Executable statement

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

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

    do bn = 1, LongBandNumber
      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, LongBandNumber
        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

  subroutine ShortFlux( &
    & xyr_TauQVap, xyr_TauDryAir, & ! (in)
    & xyr_RadSFlux )                ! (inout)
    !
    ! 短波フラックスを計算します.
    !
    ! Calculate short wave flux. 
    !
    real(DP), intent(in):: xyr_TauQVap (0:imax-1, 1:jmax, 0:kmax)
                              ! 光学的厚さ：水
                              ! Optical depth of water
    real(DP), intent(in):: xyr_TauDryAir (0:imax-1, 1:jmax, 0:kmax)
                              ! 光学的厚さ：空気
                              ! Optical depth of air
    real(DP), intent(inout):: xyr_RadSFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 短波 (日射) フラックス. 
                              ! Shortwave (insolation) flux

    ! 宣言文 ; Declaration statements
    !


    ! 作業変数
    ! Work variables
    !
    real(DP):: xy_InAngle (0:imax-1, 1:jmax)
                              ! sec (入射角). 
                              ! sec (angle of incidence)

    real(DP):: xy_SurfAlbedo (0:imax-1, 1:jmax)
                              ! 地表アルベド. 
                              ! Surface albedo

    real(DP):: BandWeightSum  ! バンドウェイトの和
                              ! Sum of band weights
    real(DP):: SecScat        ! 散乱の $ sec \zeta $
                              ! $ sec \zeta $ of scattering
    integer:: k, bn           ! DO ループ用作業変数
                              ! Work variables for DO loop

    ! 実行文 ; Executable statement

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

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

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

    do bn = 1, ShortBandNumber
      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( - SecScat &
          &            * (   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 ( ShortBandNumber == 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 RadiationInit

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

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

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

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

    implicit none

    integer:: valid_count     ! NAMELIST から読み込んだ配列の有効値の数. 
                              ! Number of valid values in an array loaded from NAMELIST
    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                              ! IOSTAT of NAMELIST read

    character(*), parameter:: subname = 'RadiationInit'

    ! 実行文 ; Executable statement

    if ( initialized ) return

    call InitCheck

    ! デフォルト値の設定 (この他は宣言文にて設定)
    ! Default values settings (Others are specified in declaration statements)
    !
    LongAbsorpCoeffQVap    = -999.9_DP
    LongAbsorpCoeffDryAir  = -999.9_DP
    LongBandWeight         = -999.9_DP

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

    ShortAbsorpCoeffQVap   = -999.9_DP
    ShortAbsorpCoeffDryAir = -999.9_DP
    ShortBandWeight        = -999.9_DP

    ShortAbsorpCoeffQVap   (1:ShortBandNum) = (/ 0.002 /)
    ShortAbsorpCoeffDryAir (1:ShortBandNum) = (/ 0.0 /)
    ShortBandWeight        (1:ShortBandNum) = (/ 1.0 /)


    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    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)

    ! ステップ間隔の処理
    ! Handle step interval
    !

    ! バンド数, 吸収係数, バンドウェイトのチェック
    ! Check number of band, absorption coefficients, band weight
    !

    if ( any( LongAbsorpCoeffQVap(1:LongBandNum) < 0.0_DP ) ) then
      valid_count = count( LongAbsorpCoeffQVap(1:LongBandNum) > -1.0_DP )
      call MessageNotify( 'E', module_name, &
        &   'number of valid values of LongAbsorpCoeffQVap=<%*f> is %d. ' \\ &
        &   'That is less than LongBandNum=<%d>.', &
        &   d = LongAbsorpCoeffQVap(1:valid_count), n = (/ valid_count /), &
        &   i = (/ valid_count, LongBandNum /) )
    end if

    if ( any( LongAbsorpCoeffDryAir(1:LongBandNum) < 0.0_DP ) ) then
      valid_count = count( LongAbsorpCoeffDryAir(1:LongBandNum) > -1.0_DP )
      call MessageNotify( 'E', module_name, &
        &   'number of valid values of LongAbsorpCoeffDryAir=<%*f> is %d. ' \\ &
        &   'That is less than LongBandNum=<%d>.', &
        &   d = LongAbsorpCoeffDryAir(1:valid_count), n = (/ valid_count /), &
        &   i = (/ valid_count, LongBandNum /) )
    end if

    if ( any( LongBandWeight(1:LongBandNum) < 0.0_DP ) ) then
      valid_count = count( LongBandWeight(1:LongBandNum) > -1.0_DP )
      call MessageNotify( 'E', module_name, &
        &   'number of valid values of LongBandWeight=<%*f> is %d. ' \\ &
        &   'That is less than LongBandNum=<%d>.', &
        &   d = LongBandWeight(1:valid_count), n = (/ valid_count /), &
        &   i = (/ valid_count, LongBandNum /) )
    end if

    if ( any( ShortAbsorpCoeffQVap(1:ShortBandNum) < 0.0_DP ) ) then
      valid_count = count( ShortAbsorpCoeffQVap(1:ShortBandNum) > -1.0_DP )
      call MessageNotify( 'E', module_name, &
        &   'number of valid values of ShortAbsorpCoeffQVap=<%*f> is %d. ' \\ &
        &   'That is less than ShortBandNum=<%d>.', &
        &   d = ShortAbsorpCoeffQVap(1:valid_count), n = (/ valid_count /), &
        &   i = (/ valid_count, ShortBandNum /) )
    end if

    if ( any( ShortAbsorpCoeffDryAir(1:ShortBandNum) < 0.0_DP ) ) then
      valid_count = count( ShortAbsorpCoeffDryAir(1:ShortBandNum) > -1.0_DP )
      call MessageNotify( 'E', module_name, &
        &   'number of valid values of ShortAbsorpCoeffDryAir=<%*f> is %d. ' \\ &
        &   'That is less than ShortBandNum=<%d>.', &
        &   d = ShortAbsorpCoeffDryAir(1:valid_count), n = (/ valid_count /), &
        &   i = (/ valid_count, ShortBandNum /) )
    end if

    if ( any( ShortBandWeight(1:ShortBandNum) < 0.0_DP ) ) then
      valid_count = count( ShortBandWeight(1:ShortBandNum) > -1.0_DP )
      call MessageNotify( 'E', module_name, &
        &   'number of valid values of ShortBandWeight=<%*f> is %d. ' \\ &
        &   'That is less than ShortBandNum=<%d>.', &
        &   d = ShortBandWeight(1:valid_count), n = (/ valid_count /), &
        &   i = (/ valid_count, ShortBandNum /) )
    end if

    ! 印字 ; Print
    !
!!$    call MessageNotify( 'M', subname, 'imax = %d', i = (/ imax /) )
!!$    call MessageNotify( 'M', subname, 'jmax = %d', i = (/ jmax /) )
!!$    call MessageNotify( 'M', subname, 'kmax = %d', i = (/ kmax /) )
!!$    call MessageNotify( 'M', subname, 'r_Sigma = %*f', &
!!$      &                 d = (/ r_Sigma(0:kmax) /), n = (/ kmax + 1 /) )


    initialized = .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_init => initialized

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: gridset_init => initialized

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: constants_init => initialized

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

    ! 実行文 ; Executable statement

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

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

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

  end subroutine InitCheck

!!$  subroutine InitDepmod
!!$    !
!!$    ! 依存モジュールの初期化
!!$    !
!!$    ! Initialize dependency modules
!!$
!!$    ! NAMELIST ファイル入力に関するユーティリティ
!!$    ! Utilities for NAMELIST file input
!!$    !
!!$    use namelist_util, only: NmlutilInit
!!$
!!$    ! 格子点設定
!!$    ! Grid points settings
!!$    !
!!$    use gridset, only: GridsetInit
!!$
!!$    ! 物理定数設定
!!$    ! Physical constants settings
!!$    !
!!$    use constants, only: ConstantsInit
!!$
!!$    ! 実行文 ; Executable statement
!!$
!!$    call NmlutilInit
!!$    call GridsetInit
!!$    call ConstantsInit
!!$
!!$  end subroutine InitDepmod

end module radiation