!= Schneider and Liu (2009) の放射モデル
!
!= Radiation model by Schneider and Liu (2009)
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: rad_SL09.f90,v 1.6 2013/05/25 06:49:44 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!
module rad_SL09
  !
  != Schneider and Liu (2009) の放射モデル
  !
  != Radiation model by Schneider and Liu (2009)
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 
  !
  ! This is a radiation model described by Schneider and Liu (2009). 
  !
  !== References
  !
  !  Schneider, T. and J. Liu, 
  !    Formation of jets and equatorial superrotation on Jupiter, 
  !    J. Atmos. Sci., 69, 579, 2009.
  !
  !== Procedures List
  !
!!$  ! RadiationFluxDennouAGCM :: 放射フラックスの計算
!!$  ! RadiationDTempDt        :: 放射フラックスによる温度変化の計算
!!$  ! RadiationFluxOutput     :: 放射フラックスの出力
!!$  ! RadiationFinalize       :: 終了処理 (モジュール内部の変数の割り付け解除)
!!$  ! ------------            :: ------------
!!$  ! RadiationFluxDennouAGCM :: Calculate radiation flux
!!$  ! RadiationDTempDt        :: Calculate temperature tendency with radiation flux
!!$  ! RadiationFluxOutput     :: Output radiation fluxes
!!$  ! RadiationFinalize       :: Termination (deallocate variables in this module)
  !
  !== NAMELIST
  !
  ! NAMELIST#rad_SL09_nml
  !

  ! USE statements
  !

  !
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! Double precision.
    &                 STRING, &  ! Strings.
    &                 TOKEN      ! Keywords.

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

  implicit none

  private


  ! Private variables
  !
  logical, save :: FlagGMIns             ! 全球平均 radiation スイッチ
                                         ! Global Average Radiation Switch
  real(DP), save :: SWOptDepAtRefPress   ! 参照圧力での短波光学的厚さ
                                         ! SW optical thickness at reference pressure
  real(DP), save :: SWRefPress           ! 短波参照圧力
                                         ! Reference pressure for SW
  real(DP), save :: SWOrd                ! 短波光学的厚さの圧力べき
                                         ! Power of press. for SW optical thickness
  real(DP), save :: LWOptDepAtRefPress   ! 参照圧力での長波光学的厚さ
                                         ! LW optical thickness at reference pressure
  real(DP), save :: LWRefPress           ! 長波参照圧力              
                                         ! Reference pressure for LW
  real(DP), save :: LWOrd                ! 長波光学的厚さの圧力べき  
                                         ! Power of press. for LW optical thickness
  real(DP), save :: SolarConst           ! 太陽定数
                                         ! Solar constant


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

  public :: RadSL09Init
  public :: RadSL09Flux

  character(*), parameter:: module_name = 'rad_SL09'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: rad_SL09.f90,v 1.6 2013/05/25 06:49:44 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

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

contains

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

  subroutine RadSL09Flux(                     &
    & xyr_Press, xyz_Press, xyz_Temp,         & ! (in)
    & xyr_RadSUwFlux, xyr_RadSDwFlux,         & ! (out)
    & xyr_RadLUwFlux, xyr_RadLDwFlux,         & ! (out)
    & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux  & ! (out)
    & )


    ! USE statements
    !

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

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & PI                    ! $ \pi $ .
                              ! 円周率.  Circular constant

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

    !
    ! Solve radiative transfer equation in two stream approximation
    !
!!$    use rad_rte_two_stream_app, only: OLD_RadRTETwoStreamAppHomogAtm
    use rad_rte_two_stream_app, only: RadRTETwoStreamAppHomogAtm

    real(DP), intent(in ) :: xyr_Press         (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ) :: xyz_Press         (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_Temp          (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out) :: xyr_RadSUwFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyr_RadSDwFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyr_RadLUwFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyr_RadLDwFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyra_DelRadLUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
    real(DP), intent(out) :: xyra_DelRadLDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)


    ! Work variables
    !
    real(DP) :: SolarFluxTOA
!!$    real(DP) :: QeRatio
!!$    real(DP) :: xyz_SSA      (0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP) :: xyz_AF       (0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP) :: xy_SurfAlbedo(0:imax-1, 1:jmax)
!!$    real(DP) :: xy_InAngle   (0:imax-1, 1:jmax)
    real(DP) :: xy_CosSZA    (0:imax-1, 1:jmax)
    real(DP) :: xyr_OptDep   (0:imax-1, 1:jmax, 0:kmax)

    real(DP) :: SSA
    real(DP) :: AF

!!$    integer  :: i
    integer  :: j
!!$    integer  :: k


    ! 初期化
    ! Initialization
    !
    if ( .not. rad_SL09_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! Short wave radiation
    !
    xyr_OptDep = SWOptDepAtRefPress * ( xyr_Press / SWRefPress )**SWOrd


    SSA = 0.8_DP
    AF  = 0.204_DP
    !   Af = 0 may be much better than 0.204 when Eddington approximation is used.
!!$    AF         = 0.0_DP


    if ( FlagGMIns ) then
      SolarFluxTOA = SolarConst / 4.0_DP
      xy_CosSZA    = 1.0_DP
    else
      SolarFluxTOA = SolarConst / PI
      do j = 1, jmax
        xy_CosSZA(:,j) = cos( y_Lat(j) )
      end do
    end if
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        if ( xy_CosSZA(i,j) > 0.0_DP ) then
!!$          xy_InAngle(i,j) = 1.0_DP / xy_CosSZA(i,j)
!!$        else
!!$          xy_InAngle(i,j) = 0.0_DP
!!$        end if
!!$      end do
!!$    end do
!!$
!!$    !   Unused variable but this is required as an argument
!!$    !
!!$    xy_SurfAlbedo = 1.0d100
!!$
!!$    call OLD_RadRTETwoStreamAppHomogAtm(                                   &
!!$      & xy_SurfAlbedo, SolarFluxTOA, xy_InAngle, SSA, AF, xyr_OptDep,  & ! (in )
!!$      & xyr_RadSFlux,                                                  & ! (out)
!!$      & FlagSemiInfAtm = .true., FlagSL09 = .true.                     & ! (in ) optional
!!$      & )
!!$
!!$    call RadRTETwoStreamAppHomogAtm(                               &
!!$      & xy_SurfAlbedo, SolarFluxTOA, xy_InAngle, SSA, AF, xyr_OptDep,  & ! (in )
!!$      & xyr_RadSUwFlux, xyr_RadSDwFlux,                                & ! (out)
!!$      & FlagSemiInfAtm = .true., FlagSL09 = .true                      & ! (in ) optional
!!$      & )

    call RadSL09SWFlux(                                     &
      & SolarFluxTOA, xy_CosSZA, SSA, AF, xyr_OptDep,       &
      & xyr_RadSUwFlux                                      &
      & )
    xyr_RadSDwFlux = 0.0_DP



    ! Long wave radiation
    !

    !   Although the surface temperature and surface emissivity are set, but are not used.
    !
    xyr_OptDep = LWOptDepAtRefPress * ( xyr_Press / LWRefPress )**LWOrd


!!$    call RadiationRTEQNonScat(                    &
!!$      & xyz_Temp, xy_SurfTemp, xy_SurfEmis, xyr_OptDep, & ! (in )
!!$      & xyr_RadLFlux, xyra_DelRadLFlux,                 & ! (out)
!!$      & xy_SurfUpRadLFluxBase = xyr_RadSFlux(:,:,0)     & ! (in ) optional
!!$      & )
!!$    call RadSL09LWFlux(                           &
!!$      & xyz_Temp, xyr_OptDep,                     & ! (in )
!!$      & xyr_RadSFlux(:,:,0),                      & ! (in )
!!$      & xyr_RadLFlux, xyra_DelRadLFlux            & ! (out)
!!$      & )
    call RadSL09LWFlux(                       &
      & xyr_Press, xyz_Press, xyz_Temp, xyr_OptDep,  & ! (in )
      & xyr_RadSUwFlux(:,:,0)-xyr_RadSDwFlux(:,:,0), & ! (in )
      & xyr_RadLUwFlux, xyr_RadLDwFlux,              & ! (out)
      & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux       & ! (out)
      & )


  end subroutine RadSL09Flux

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

  subroutine OLD_RadSL09LWFlux(                           &
    & xyz_Temp, xyr_OptDep,                           & ! (in )
    & xy_SurfUpRadLFluxBase,                          & ! (in )
    & xyr_RadLFlux, xyra_DelRadLFlux                  & ! (out)
    & )
    !
    ! 散乱なしの場合の放射伝達方程式の計算
    !
    ! Integrate radiative transfer equation without scattering
    !

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

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

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & PI, &
                              ! $ \pi $ .
                              ! 円周率.  Circular constant
      & StB
                              ! $ \sigma_{SB} $ .
                              ! ステファンボルツマン定数.
                              ! Stefan-Boltzmann constant

    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in ) :: xyz_Temp    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(in ) :: xyr_OptDep  (0:imax-1, 1:jmax, 0:kmax)
                              ! Optical depth
    real(DP), intent(in ) :: xy_SurfUpRadLFluxBase(0:imax-1, 1:jmax)
    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), parameter :: DiffFact = 1.66_DP

    real(DP):: xyr_RadLDoFlux (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_RadLUpFlux (0:imax-1, 1:jmax, 0:kmax)

    real(DP):: xyz_DelTrans (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyrr_Trans   (0:imax-1, 1:jmax, 0:kmax, 0:kmax)
                              ! 透過係数. 
                              ! Transmission coefficient
    real(DP):: xyz_IntPF        (0:imax-1, 1:jmax, 1:kmax)
                              ! Integrated Planck function
    real(DP):: xyz_IntDPFDT     (0:imax-1, 1:jmax, 1:kmax)
                              ! Integrated temperature derivative of Planck function

    real(DP):: xy_SurfUpRadLFlux(0:imax-1, 1:jmax)

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

    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. rad_SL09_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    ! Case for grey atmosphere
    !
    xyz_IntPF       = StB * xyz_Temp**4
    xyz_IntDPFDT    = 4.0_DP * StB * xyz_Temp**3


    ! 透過関数計算
    ! Calculate transmission functions
    !
    do k = 1, kmax
      xyz_DelTrans(:,:,k) = &
        & exp( - DiffFact * ( xyr_OptDep(:,:,k-1) - xyr_OptDep(:,:,k) ) )
    end do
    !
    do k = 0, kmax
      do kk = k, k
        xyrr_Trans(:,:,k,kk) = 1.0_DP
      end do
      do kk = k+1, kmax
        xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,k,kk-1) * xyz_DelTrans(:,:,kk)
      end do
    end do
    do k = 0, kmax
      do kk = 0, k-1
        xyrr_Trans(:,:,k,kk) = xyrr_Trans(:,:,kk,k)
      end do
    end do


    ! 放射フラックス計算
    ! Calculate radiation flux
    !

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

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

    end do
    !
    !   Upward flux
    !
    !     Set upward flux
    !
    xy_SurfUpRadLFlux = xyr_RadLDoFlux(:,:,0) - xy_SurfUpRadLFluxBase
    !
    do k = 0, kmax

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

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

    end do

    xyr_RadLFlux = xyr_RadLUpFlux - xyr_RadLDoFlux


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

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

  end subroutine OLD_RadSL09LWFlux

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

  subroutine RadSL09LWFlux(                       &
    & xyr_Press, xyz_Press, xyz_Temp, xyr_OptDep, & ! (in )
    & xy_SurfRadSFlux,                            & ! (in )
    & xyr_RadLUwFlux, xyr_RadLDwFlux,             & ! (out)
    & xyra_DelRadLUwFlux, xyra_DelRadLDwFlux      & ! (out)
    & )
    !
    ! 散乱なしの場合の放射伝達方程式の計算
    !
    ! Integrate radiative transfer equation without scattering
    !

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

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

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

    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in ) :: xyr_Press   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ) :: xyz_Press   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_Temp    (0:imax-1, 1:jmax, 1:kmax)
                              ! $ T $ .     温度. Temperature
    real(DP), intent(in ) :: xyr_OptDep  (0:imax-1, 1:jmax, 0:kmax)
                              ! Optical depth
    real(DP), intent(in ) :: xy_SurfRadSFlux(0:imax-1, 1:jmax)
    real(DP), intent(out) :: xyr_RadLUwFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux
    real(DP), intent(out) :: xyr_RadLDwFlux (0:imax-1, 1:jmax, 0:kmax)
                              ! 長波フラックス. 
                              ! Longwave flux
    real(DP), intent(out) :: xyra_DelRadLUwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with longwave
    real(DP), intent(out) :: xyra_DelRadLDwFlux (0:imax-1, 1:jmax, 0:kmax, 0:1)
                              ! 長波地表温度変化. 
                              ! Surface temperature tendency with longwave


    ! 作業変数
    ! Work variables
    !
    real(DP):: xyr_Temp       (0:imax-1, 1:jmax, 0:kmax)
    real(DP):: xyr_IntPF      (0:imax-1, 1:jmax, 0:kmax)
                              ! Integrated Planck function
    real(DP):: xyz_DPFDOptDep (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_DelOptDep     (0:imax-1, 1:jmax, 1:kmax)
    real(DP):: xyz_TransEachLayer(0:imax-1, 1:jmax, 1:kmax)
                              ! 透過係数. 
                              ! Transmission coefficient

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

    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. rad_SL09_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    k = 0
    xyr_Temp(:,:,k) =                                    &
      &      ( xyz_Temp (:,:,k+2) - xyz_Temp (:,:,k+1) ) &
      & / log( xyz_Press(:,:,k+2) / xyz_Press(:,:,k+1) ) &
      & * log( xyr_Press(:,:,k  ) / xyz_Press(:,:,k+1) ) &
      & + xyz_Temp(:,:,k+1)
    do k = 1, kmax-1
      xyr_Temp(:,:,k) =                                &
        &      ( xyz_Temp (:,:,k+1) - xyz_Temp (:,:,k) ) &
        & / log( xyz_Press(:,:,k+1) / xyz_Press(:,:,k) ) &
        & * log( xyr_Press(:,:,k  ) / xyz_Press(:,:,k) ) &
        & + xyz_Temp(:,:,k)
    end do
    k = kmax
    xyr_Temp(:,:,k) = xyz_Temp(:,:,k)


    xyr_IntPF       = StB * xyr_Temp**4
    !
    do k = 1, kmax
      xyz_DelOptDep(:,:,k) = xyr_OptDep(:,:,k-1) - xyr_OptDep(:,:,k)
    end do
    !
    xyz_TransEachLayer = exp( - xyz_DelOptDep )
    !
    do k = 1, kmax
      xyz_DPFDOptDep(:,:,k) =                      &
        &   ( xyr_IntPF(:,:,k-1) - xyr_IntPF(:,:,k) ) &
        & / max( xyz_DelOptDep(:,:,k), 1.0d-100 )
    end do


    ! 放射フラックス計算
    ! Calculate radiation flux
    !

    !
    !   Downward flux
    !
    k = kmax
    xyr_RadLDwFlux(:,:,k) = 0.0_DP
    do k = kmax-1, 0, -1
      xyr_RadLDwFlux(:,:,k) =                                &
        &   ( xyr_RadLDwFlux(:,:,k+1) - xyr_IntPF(:,:,k+1) ) &
        & * xyz_TransEachLayer(:,:,k+1)                      &
        & + xyr_IntPF(:,:,k)                                 &
        & - xyz_DPFDOptDep(:,:,k+1)                          &
        &   * ( 1.0_DP - xyz_TransEachLayer(:,:,k+1) )
    end do
    !
    !   Upward flux
    !
    !     Set upward flux
    k = 0
    xyr_RadLUwFlux(:,:,k) = xyr_RadLDwFlux(:,:,0) - xy_SurfRadSFlux
    !
    do k = 1, kmax
      xyr_RadLUwFlux(:,:,k) =                                &
        &   ( xyr_RadLUwFlux(:,:,k-1) - xyr_IntPF(:,:,k-1) ) &
        & * xyz_TransEachLayer(:,:,k)                        &
        & + xyr_IntPF(:,:,k)                                 &
        & + xyz_DPFDOptDep(:,:,k)                            &
        &   * ( 1.0_DP - xyz_TransEachLayer(:,:,k) )
    end do


    ! 放射フラックスの変化率の計算
    ! Calculate rate of change of radiative flux
    !
    xyra_DelRadLUwFlux = 0.0_DP
    xyra_DelRadLDwFlux = 0.0_DP


  end subroutine RadSL09LWFlux

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

  subroutine RadSL09SWFlux(                                     &
    & SolarFluxTOA, xy_CosSZA, SSA, AF, xyr_OptDep,             &
    & xyr_RadSFlux                                              &
    & )

    real(DP), intent(in ) :: SolarFluxTOA
    real(DP), intent(in ) :: xy_CosSZA(0:imax-1, 1:jmax)
    real(DP), intent(in ) :: SSA
    real(DP), intent(in ) :: AF
    real(DP), intent(in ) :: xyr_OptDep   ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(out) :: xyr_RadSFlux ( 0:imax-1, 1:jmax, 0:kmax )

    ! Work variables
    !
    real(DP) :: BondAlbedo
    real(DP) :: Gamma
    integer  :: j, k


    BondAlbedo = &
      &   ( sqrt( 1.0_DP - SSA * AF ) - sqrt( 1.0_DP - SSA ) ) &
      & / ( sqrt( 1.0_DP - SSA * AF ) + sqrt( 1.0_DP - SSA ) )

    Gamma = 2.0_DP * sqrt( 1.0_DP - SSA ) * sqrt( 1.0_DP - SSA * AF )

    do k = 0, kmax
      do j = 1, jmax
        xyr_RadSFlux(:,j,k) =                   &
          & - SolarFluxTOA * xy_CosSZA(:,j) &
          & * ( 1.0_DP - BondAlbedo ) * exp( -Gamma * xyr_OptDep(:,j,k) )
      end do
    end do


  end subroutine RadSL09SWFlux

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

  subroutine RadSL09Init

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

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

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

    !
    ! Solve radiative transfer equation in two stream approximation
    !
    use rad_rte_two_stream_app, only: RadRTETwoStreamAppInit


    ! 宣言文 ; Declaration statements
    !

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

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /rad_SL09_nml/ &
      & FlagGMIns,                             &
      & SWOptDepAtRefPress, SWRefPress, SWOrd, &
      & LWOptDepAtRefPress, LWRefPress, LWOrd, &
      & SolarConst
          !
          ! デフォルト値については初期化手続 "rad_SL09#RadSL09Init"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "rad_SL09#RadSL09Init" for the default values.
          !


    if ( rad_SL09_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !
    FlagGMIns          = .false.

    SWOptDepAtRefPress =  3.0_DP
    SWRefPress         =  3.0d5
    SWOrd              =  1.0_DP

    LWOptDepAtRefPress = 80.0_DP
    LWRefPress         =  3.0d5
    LWOrd              =  2.0_DP

    SolarConst         = 50.7_DP




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

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

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


    ! Initialization of modules used in this module
    !

    !
    ! Solve radiative transfer equation in two stream approximation
    !
    call RadRTETwoStreamAppInit


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'FlagGMIns          = %b', l = (/ FlagGMIns /) )
    call MessageNotify( 'M', module_name, 'SWOptDepAtRefPress = %f', d = (/ SWOptDepAtRefPress /) )
    call MessageNotify( 'M', module_name, 'SWRefPress         = %f', d = (/ SWRefPress /) )
    call MessageNotify( 'M', module_name, 'SWOrd              = %f', d = (/ SWOrd /) )
    call MessageNotify( 'M', module_name, 'LWOptDepAtRefPress = %f', d = (/ LWOptDepAtRefPress /) )
    call MessageNotify( 'M', module_name, 'LWRefPress         = %f', d = (/ LWRefPress /) )
    call MessageNotify( 'M', module_name, 'LWOrd              = %f', d = (/ LWOrd /) )
    call MessageNotify( 'M', module_name, 'SolarConst         = %f', d = (/ SolarConst /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    rad_SL09_inited = .true.

  end subroutine RadSL09Init

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

end module rad_SL09
