subroutine RadiationSL09Flux( xyr_Press, xyz_Temp, xyr_RadSFlux, xyr_RadLFlux, xyra_DelRadLFlux )


    ! USE statements
    !

    !
    ! Physical constants settings
    !
    use constants, only: PI       ! $ \pi $ .
                                  ! Circular constant

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

    !
    ! Solve radiative transfer equation in two stream approximation
    !
    use radiation_two_stream_app, only: RadiationTwoStreamApp, IDSchemeShortWave, RadiationTwoStreamAppHomogAtm

    ! 放射関連ルーチン
    ! Routines for radiation calculation
    !
    use radiation_utils, only : RadiationRTEQNonScat

    real(DP), intent(in ) :: xyr_Press    ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(in ) :: xyz_Temp     ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(out) :: xyr_RadSFlux ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP), intent(out) :: xyr_RadLFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyra_DelRadLFlux(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_SurfTemp  (0:imax-1, 1:jmax)
    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) :: xy_SurfEmis  (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. radiation_SL09_inited ) call RadiationSL09Init


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

    SolarFluxTOA  = SolarConst / PI


    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


    do j = 1, jmax
      xy_CosSZA(:,j) = cos( y_Lat(j) )
    end do
    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 RadiationTwoStreamAppHomogAtm( xy_SurfAlbedo, SolarFluxTOA, xy_InAngle, SSA, AF, xyr_OptDep, xyr_RadSFlux, FlagSemiInfAtm = .true., FlagSL09 = .true. )



    ! Long wave radiation
    !

    !   Although the surface temperature and surface emissivity are set, but are not used.
    !
    xy_SurfTemp = 1.0d100
    xy_SurfEmis = 1.0d100

    xyr_OptDep = LWOptDepAtRefPress * ( xyr_Press / LWRefPress )**LWOrd


    call RadiationRTEQNonScat( xyz_Temp, xy_SurfTemp, xy_SurfEmis, xyr_OptDep, xyr_RadLFlux, xyra_DelRadLFlux, xy_SurfUpRadLFluxBase = xyr_RadSFlux(:,:,0) )


  end subroutine RadiationSL09Flux

  subroutine RadiationSL09Init

    ! ファイル入出力補助
    ! 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


    ! 宣言文 ; 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 /radiation_SL09_nml/ SWOptDepAtRefPress, SWRefPress, SWOrd, LWOptDepAtRefPress, LWRefPress, LWOrd, SolarConst
          !
          ! デフォルト値については初期化手続 "radiation_SL09#RadiationSL09Init"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "radiation_SL09#RadiationSL09Init" for the default values.
          !


    ! デフォルト値の設定
    ! Default values settings
    !
    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, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = radiation_SL09_nml, iostat = iostat_nml )             ! (out)
      close( unit_nml )

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



    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    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) )


    radiation_SL09_inited = .true.

  end subroutine RadiationSL09Init

  subroutine RadiationSL09SWFlux( 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 RadiationSL09SWFlux