!= dcpam 地球大気向け短波放射モデル Ver. 1
!
!= dcpam short wave radiation model for the Earth's atmosphere Ver. 1
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: radiation_dcpam_E_SW_V1.f90,v 1.1 2010-01-11 01:28:10 yot Exp $
! Tag Name::  $Name: dcpam5-20100224 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!
module radiation_dcpam_E_SW_V1

  !
  != dcpam 地球大気向け短波放射モデル Ver. 1
  !
  != dcpam short wave radiation model for the Earth's atmosphere Ver. 1
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 地球大気向け短波放射モデル.
  !
  ! This is a short wave radiation model for the Earth's atmospehre.
  !
  ! The wavenumber range of shortwave radiation treated by this routine is 
  ! from 2600 to 66667 cm-1 (0.15 to 3.85 micron). 
  ! The wavenumber range from 40 to 2600 cm-1 is treated as routine of 
  ! long wave radiation. 
  !
  ! From 2600 to 66667 cm-1, the following effects are considered. 
  ! * 2600 to 12040 cm-1 (0.83-3.85 micron): 
  !   * absorption by H2O, 
  !     * absorption by H2O is considered by using k-distribution method 
  !       following Chou and Arking (1981), 
  !   * scattering by cloud droplets.
  ! * 12040 to 66667 cm-1 (0.15 to 0.83 micron): 
  !   * Rayleigh scattering, 
  !   * scattering by cloud droplets. 
  !
  !
  !== References
  !
  ! Chou, M.-D., and A. Arking, 
  ! An efficient method for computing the absorption of solar radiation by water vapor, 
  ! J. Atmos. Sci., 38, 798-807, 1981.
  !
  ! Mount, G. H., and  G. J. Rottman, 
  ! The solar absolute spectral irradiance 1150-3173 A: May 17, 1982, 
  ! ???, 88, 5403-5410, 1983.
  !
  !== 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#radiation_DennouAGCM_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

  real(DP), parameter :: BoltzConst = 1.3806503d-23
  real(DP), parameter :: AvogNum    = 6.0221415d23

  integer , parameter :: IDSchemeShortWave = 1
  integer , parameter :: IDSchemeLongWave  = 2

!!$  integer , parameter :: nswband = 7
!!$  integer , parameter :: nswband = 11
  integer , parameter :: nswband = 10
  real(DP), save      :: SolarSpecWavelengthLimitMicron( 1:nswband, 1:2 )
  real(DP), save      :: SolarSpecWavelength           ( 1:nswband )
  real(DP), save      :: SolarSpec                     ( 1:nswband )

  real(DP), save      :: O3AbsCrossSection             ( 1:nswband )

  real(DP), save      :: MeanMolWeight
  real(DP), save      :: O3MolWeight

  logical , save      :: FlagRayleighScattering


!!$  ! Table 2.3 of Liou (2002)
!!$  !
!!$  ! Unit of wavelength is micron meters.
!!$  !
!!$  data SolarSpecWavelength &
!!$    & / 0.2d0   ,  0.3d0   ,  0.4d0   ,  0.5d0   ,  0.6d0   ,  0.7d0   ,  0.8d0    /
!!$  !
!!$  data SolarSpecWavelengthLimitMicron &
!!$      & / 0.15d0  ,  0.25d0  ,  0.35d0  ,  0.45d0  ,  0.55d0  ,  0.65d0  ,  0.75d0,  &
!!$      &   0.25d0  ,  0.35d0  ,  0.45d0  ,  0.55d0  ,  0.65d0  ,  0.75d0  ,  0.83d0   /
!!$  !
!!$  ! Unit of solar flux is W m-2 micron-1
!!$  !
!!$  data SolarSpec &
!!$    & / 2.0832d1,  5.4765d2,  1.4042d3,  1.9619d3,  1.7632d3,  1.4300d3,  1.1257d3 /
!!$
!!$
!!$  !
!!$  ! The original data is obtained from a web page of IGACO Activity. 
!!$  ! http://igaco-o3.fmi.fi/ACSO/cross_sections.html
!!$  ! FTS-Voigt o3_203l.dat
!!$  ! Unit of solar flux is m2 molecule-1.
!!$  !
!!$!  ! o3_203l.dat
!!$!  data O3AbsCrossSection &
!!$!    & / 1.52d-22,  3.04d-22,  1.00d-26,  1.46d-25,  4.28d-25,  1.10d-25,  2.03d-26 /
!!$  ! o3_246l.dat
!!$  data O3AbsCrossSection &
!!$    & / 1.66d-22,  3.05d-22,  1.61d-27,  1.42d-25,  4.23d-25,  1.15d-25,  2.80d-26 /




!!$  !
!!$  ! Unit of wavelength is micron meters.
!!$  !
!!$  data SolarSpecWavelength &
!!$    & / 0.2d0   ,  0.26d0  ,  0.28d0  ,  0.30d0  ,  0.32d0  ,  0.34d0  ,  0.4d0   ,  0.5d0   ,  0.6d0   ,  0.7d0   ,  0.8d0    /
!!$  data SolarSpecWavelengthLimitMicron &
!!$    & / 0.15d0  ,  0.25d0  ,  0.27d0  ,  0.29d0  ,  0.31d0  ,  0.33d0  , 0.35d0  ,  0.45d0  ,  0.55d0  ,  0.65d0  ,  0.75d0,  &
!!$    &   0.25d0  ,  0.27d0  ,  0.29d0  ,  0.31d0  ,  0.33d0  ,  0.35d0  , 0.45d0  ,  0.55d0  ,  0.65d0  ,  0.75d0  ,  0.83d0   /
!!$  !
!!$  ! Unit of solar flux is W m-2 micron-1
!!$  !
!!$  data SolarSpec &
!!$    & / 2.0832d1,  5.4765d2,  5.4765d2,  5.4765d2,  5.4765d2,  5.4765d2,  1.4042d3,  1.9619d3,  1.7632d3,  1.4300d3,  1.1257d3 /
!!$
!!$
!!$  !
!!$  ! The original data is obtained from a web page of IGACO Activity. 
!!$  ! http://igaco-o3.fmi.fi/ACSO/cross_sections.html
!!$  ! FTS-Voigt o3_203l.dat
!!$  ! Unit of solar flux is m2 molecule-1.
!!$  !
!!$!  ! o3_203l.dat
!!$!  data O3AbsCrossSection &
!!$!    & / 1.52d-22,  3.04d-22,  1.00d-26,  1.46d-25,  4.28d-25,  1.10d-25,  2.03d-26 /
!!$  ! o3_246l.dat
!!$  data O3AbsCrossSection &
!!$    & / 1.66d-22,  1.05d-21,  4.16d-22,  4.92d-23,  3.11d-24, 1.53d-25,  1.61d-27,  1.42d-25,  4.23d-25,  1.15d-25,  2.80d-26 /




  !
  ! Solar spectrum from 2600 to 12040 cm-1 (0.83-3.85 micron) is originally obtained 
  ! from table 2.3 of Liou (2002), then subivided in the wavelength range from 0.25 
  ! to 0.35 micron by using the solar spectrum data by Mount and Rottman (1983).
  !

  !
  ! Unit of wavelength is micron meters.
  !
  data SolarSpecWavelength &
    & / 0.2d0   ,  0.26d0  ,  0.28d0  ,  0.30d0  ,  0.33d0  ,  0.4d0   ,  0.5d0   ,  0.6d0   ,  0.7d0   ,  0.8d0    /
  data SolarSpecWavelengthLimitMicron &
    & / 0.15d0  ,  0.25d0  ,  0.27d0  ,  0.29d0  ,  0.31d0  ,  0.35d0  ,  0.45d0  ,  0.55d0  ,  0.65d0  ,  0.75d0,  &
    &   0.25d0  ,  0.27d0  ,  0.29d0  ,  0.31d0  ,  0.35d0  ,  0.45d0  ,  0.55d0  ,  0.65d0  ,  0.75d0  ,  0.83d0   /
  !
  ! Unit of solar flux is W m-2 micron-1
  !
  data SolarSpec &
    & / 2.0832d1,  1.4088d2,  2.4832d2,  5.4289d2,  1.2585d3,  1.4042d3,  1.9619d3,  1.7632d3,  1.4300d3,  1.1257d3 /

  !
  ! The original data is obtained from a web page of IGACO Activity. 
  ! http://igaco-o3.fmi.fi/ACSO/cross_sections.html
  ! FTS-Voigt o3_246l.dat
  ! O3 absorption cross-section at 246 K and at 100 hPa is used. 

  !
  ! Unit of solar flux is m2 molecule-1.
  !
!  ! o3_203l.dat
!  data O3AbsCrossSection &
!    & / 1.52d-22,  3.04d-22,  1.00d-26,  1.46d-25,  4.28d-25,  1.10d-25,  2.03d-26 /
  ! o3_246l.dat
  data O3AbsCrossSection &
    & / 1.66d-22,  1.05d-21,  4.18d-22,  4.92d-23,  1.63d-24, 1.61d-27,  1.42d-25,  4.23d-25,  1.15d-25,  2.80d-26 /



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

  public :: RadiationDcpamESWV1Flux

  character(*), parameter:: module_name = 'radiation_dcpam_E_SW_V1'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name: dcpam5-20100224 $' // &
    & '$Id: radiation_dcpam_E_SW_V1.f90,v 1.1 2010-01-11 01:28:10 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

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

contains

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

  subroutine RadiationDcpamESWV1Flux(                                 &
    & xy_SurfAlbedo,                                                  &
    & xyz_Press, xyr_Press, xyz_Temp, xyz_QVap, xyz_QO3, xyz_Height,  &
    & xyr_RadSFlux                                                    &
    & )

    ! USE statements
    !

    !
    ! Physical constants settings
    !
    use constants, only: Grav,  & ! $ g $ [m s-2].
                                  !
                                  ! Gravitational acceleration
      &                  PI       ! $ \pi $ .
                                  ! Circular constant

    use radiation_two_stream_app, only: RadiationTwoStreamApp

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN, &              ! ステップ $ t $ の時刻.
                              ! Time of step $ t $.
      & EndTime, &            ! 計算終了時刻.
                              ! End time of calculation
      & TimesetClockStart, TimesetClockStop

    ! Cho and Arking (1981) による短波放射モデル
    ! Short wave radiation model described by Lacis and Hansen (1974)
    !
    use radiation_CA81, only: RadiationCA81Flux, RadiationCA81NumKDFBin, RadiationCA81H2ODelOptDep

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

    use set_cloud, only : SetCloudSW


    real(DP), intent(in ) :: xy_SurfAlbedo( 0:imax-1, 1:jmax )
    real(DP), intent(in ) :: xyz_Press    ( 0:imax-1, 1:jmax, 1:kmax )
    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(in ) :: xyz_QVap     ( 0:imax-1, 1:jmax, 1:kmax )
                              ! $ q $ .   混合比. Mass mixing ratio of constituents (1)
    real(DP), intent(in ) :: xyz_QO3      ( 0:imax-1, 1:jmax, 1:kmax )
                              ! O3 分布 (1)
                              ! O3 distribution (1)
    real(DP), intent(in ) :: xyz_Height   ( 0:imax-1, 1:jmax, 1:kmax )
    real(DP), intent(out) :: xyr_RadSFlux ( 0:imax-1, 1:jmax, 0:kmax )


    real(DP):: xy_IncomRadSFlux(0:imax-1, 1:jmax)
                              ! 短波 (日射) フラックス.
                              ! Short wave (insolation) flux
    real(DP) :: xy_InAngle    (0:imax-1, 1:jmax)
                              ! sec (入射角).
                              ! sec (angle of incidence)
    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) :: SolarFluxTOA
    integer  :: IDScheme

    real(DP), parameter :: CloudSinScatAlb   = 1.0d0 - 1.0d-10
    real(DP), parameter :: CloudAsymFact     = 0.85d0
    real(DP), parameter :: RayScatSinScatAlb = 1.0d0 - 1.0d-10
    real(DP), parameter :: RayScatAsymFact   = 0.0d0

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

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

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

    real(DP)            :: xyz_TotDelOptDep( 0:imax-1, 1:jmax, 1:kmax )
    real(DP)            :: xyr_TotOptDep   ( 0:imax-1, 1:jmax, 0:kmax )
    real(DP)            :: xyr_RadFlux     ( 0:imax-1, 1:jmax, 0:kmax )

    real(DP) :: Wavelength
    real(DP) :: NumDenSTP
    real(DP) :: RefIndexReal
    real(DP) :: CorFactorDelta
    real(DP) :: CorFactor
    real(DP) :: RayScatCrossSec

    integer  :: nCA81KDFBin
    real(DP) :: KDFWeight
    real(DP) :: xyz_H2ODelOptDep( 0:imax-1, 1:jmax, 1:kmax )

    logical  :: flag_cloud

    integer  :: i
    integer  :: j
    integer  :: k
    integer  :: l


    ! 初期化
    ! Initialization
    !
    if ( .not. radiation_dcpam_E_SW_V1_inited ) call RadiationDcpamESWV1Init


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

!!$    write( 6, * ) '************************'
!!$    write( 6, * ) '************************'
!!$    write( 6, * ) '************************'
!!$    write( 6, * ) '************************'
!!$    write( 6, * ) '************************'
!!$    write( 6, * ) 'xy_In_Angle is set'
!!$    xy_InAngle = 1.0d0
!!$    xy_InAngle = 1.0d0 / cos( 60.0d0 * 3.141592d0 / 180.0d0 )
!!$    write( 6, * ) '************************'
!!$    write( 6, * ) '************************'
!!$    write( 6, * ) '************************'
!!$    write( 6, * ) '************************'
!!$    write( 6, * ) '************************'




    ! Cloud optical depth
    !

    call SetCloudSW(                                       &
      & xyz_Press, xyr_Press, xyz_Temp, xyz_QVap, xyz_Height,  &
      & xyz_CloudDelOptDep                                     &
      & )


    xyr_RadSFlux = 0.0d0


    ! * 12040 to 66667 cm-1 (0.15 to 0.83 micron): 
    !   * Rayleigh scattering, 
    !   * scattering by cloud droplets. 
    !
    do l = 1, nswband

      IDScheme       = IDSchemeShortWave

      QeRatio        = 1.0d0

      SolarFluxTOA   = SolarSpec(l)
      Wavelength     = SolarSpecWavelength(l)


!!$      xyz_ssa        = 1.0d0 - 1.0d-10
!!$      xyz_af         = 0.0d0



      ! Rayleigh scattering
      !
      if ( FlagRayleighScattering ) then

        ! (3.3.17) in p.93
        RefIndexReal =                                             &
          &      6432.8d0                                          &
          & + 2949810.0d0 / ( 146.0d0 - (Wavelength*1.0d6)**(-2) ) &
          & +   25540.0d0 / (  41.0d0 - (Wavelength*1.0d6)**(-2) )
        RefIndexReal = 1.0d0 + RefIndexReal * 1.0d-8
        !
        ! equation in text in p.93
        CorFactorDelta = 0.035d0
        CorFactor      =                         &
          &   ( 6.0d0 + 3.0d0 * CorFactorDelta ) &
          & / ( 6.0d0 - 7.0d0 * CorFactorDelta )
        !
        ! (3.3.19) in p.93
        NumDenSTP = 101325.0d0 / ( BoltzConst * 273.15d0 )
        RayScatCrossSec =                                     &
          &   8.0d0 * PI**3 * ( RefIndexReal**2 - 1.0d0 )**2  &
          & / ( 3.0d0 * Wavelength**4 * NumDenSTP**2 )        &
          & * CorFactor
        !
        do k = 1, kmax
          xyz_RayScatDelOptDep(:,:,k) =                            &
            & + ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k ) ) / Grav  &
            & / ( MeanMolWeight / AvogNum )                        &
            & * RayScatCrossSec
        end do

      else

        xyz_RayScatDelOptDep = 0.0d0

      end if



      ! O3 absorption
      !
      do k = 1, kmax
        xyz_O3AbsDelOptDep(:,:,k) =                              &
          & + ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k ) ) / Grav  &
          & / ( O3MolWeight / AvogNum )                          &
          & * xyz_QO3(:,:,k)                                     &
          & * O3AbsCrossSection(l)
      end do



      ! Total optical parameter
      !
      xyz_TotDelOptDep = xyz_CloudDelOptDep + xyz_RayScatDelOptDep + xyz_O3AbsDelOptDep
      !
      xyr_TotOptDep(:,:,kmax) = 0.0d0
      do k = kmax-1, 0, -1
        xyr_TotOptDep(:,:,k) = xyr_TotOptDep(:,:,k+1) + xyz_TotDelOptDep(:,:,k+1)
      end do
      !
      xyz_ssa =                                           &
        &   ( CloudSinScatAlb   * xyz_CloudDelOptDep      &
        &   + RayScatSinScatAlb * xyz_RayScatDelOptDep    &
        &   + 0.0d0             * xyz_O3AbsDelOptDep )    &
        & / ( xyz_TotDelOptDep + 1.0d-100 )
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if ( xyz_ssa(i,j,k) >= 1.0d0 ) then
              xyz_ssa(i,j,k) = 1.0d0 - 1.0d-10
            end if
          end do
        end do
      end do
      xyz_af  =                                                               &
        &   ( CloudAsymFact     * CloudSinScatAlb   * xyz_CloudDelOptDep      &
        &   + RayScatAsymFact   * RayScatSinScatAlb * xyz_RayScatDelOptDep    &
        &   + 0.0d0             * 0.0d0             * xyz_O3AbsDelOptDep   )  &
        & / ( xyz_ssa * xyz_TotDelOptDep + 1.0d-100 )





!!$      ! Total optical parameter
!!$      !
!!$      xyz_TotDelOptDep = xyz_O3AbsDelOptDep
!!$
!!$      xyr_TotOptDep(:,:,kmax) = 0.0d0
!!$      do k = kmax-1, 0, -1
!!$        xyr_TotOptDep(:,:,k) = xyr_TotOptDep(:,:,k+1) + xyz_TotDelOptDep(:,:,k+1)
!!$      end do
!!$
!!$      xyz_ssa = 0.0d0
!!$      xyz_af  = 0.0d0



      call RadiationTwoStreamApp( &
        & QeRatio, xyz_ssa, xyz_af, &
        & SolarFluxTOA, &
        & xy_SurfAlbedo, &
        & IDScheme, &
        & xy_InAngle, &
        & xyr_TotOptDep, &
        & xyr_RadFlux &
        & )

      xyr_RadSFlux = xyr_RadSFlux + xyr_RadFlux

!!$      i = imax / 2
!!$      j = jmax / 2
!!$      do k = 0, kmax
!!$        write( 91, * ) l, Wavelength, xyr_RadFlux(i,j,k), xyr_Press(i,j,k), xyr_OptDepBase(i,j,k)
!!$      end do
!!$      write( 91, * )
!!$      write( 91, * )
!!$      call flush( 91 )

    end do

!!$    i = imax / 2
!!$    j = jmax / 2
!!$    j = jmax * 3 / 4
!!$    do k = 0, kmax
!!$      write( 92, * ) xyr_RadSFlux(i,j,k), xyr_Press(i,j,k)
!!$    end do
!!$    call flush( 92 )
!!$
!!$    do k = 1, kmax
!!$      write( 82, * ) &
!!$        & + (     xyr_RadSFlux(i,j,k-1) - xyr_RadSFlux(i,j,k) )  &
!!$        &     / ( xyr_Press(i,j,k-1)    - xyr_Press(i,j,k) )     &
!!$        &     / 1004.6 * Grav, &
!!$        & xyz_Press(i,j,k)
!!$    end do
!!$    call flush( 82 )
!!$
!!$    write( 6, * ) '########## ', acos( 1.0d0 / xy_InAngle(i,j) ) * 180.0d0 / 3.141592d0


    ! Check for dry atmosphere
    !
    if ( all( xyz_CloudDelOptDep <= 0.0d0 ) ) then
      flag_cloud = .false.
!!$      write( 6, * ) 'SHORTWAVE TEST: No cloud'
    else
      flag_cloud = .true.
    end if


    ! * 2600 to 12040 cm-1 (0.83-3.85 micron): 
    !   * absorption by H2O, 
    !     * absorption by H2O is considered by using k-distribution method 
    !       following Chou and Arking (1981), 
    !   * scattering by cloud droplets.
    !
    if ( flag_cloud ) then

      call RadiationCA81NumKDFBin( nCA81KDFBin )

      do l = 1, nCA81KDFBin

        call RadiationCA81H2ODelOptDep(        &
          & xyz_QVap, xyr_Press, xyz_Press, l, & ! (in )
          & xyz_H2ODelOptDep, KDFWeight              & ! (out)
          & )


        ! Total optical parameter
        !
        xyz_TotDelOptDep = xyz_CloudDelOptDep + xyz_H2ODelOptDep

        xyr_TotOptDep(:,:,kmax) = 0.0d0
        do k = kmax-1, 0, -1
          xyr_TotOptDep(:,:,k) = xyr_TotOptDep(:,:,k+1) + xyz_TotDelOptDep(:,:,k+1)
        end do

        xyz_ssa =                                           &
          &   ( CloudSinScatAlb   * xyz_CloudDelOptDep      &
          &   + 0.0d0             * xyz_H2ODelOptDep     )  &
          & / ( xyz_TotDelOptDep + 1.0d-100 )
        do k = 1, kmax
          do j = 1, jmax
            do i = 0, imax-1
              if ( xyz_ssa(i,j,k) >= 1.0d0 ) then
                xyz_ssa(i,j,k) = 1.0d0 - 1.0d-10
              end if
            end do
          end do
        end do
        xyz_af  =                                                               &
          &   ( CloudAsymFact     * CloudSinScatAlb   * xyz_CloudDelOptDep      &
          &   + 0.0d0             * 0.0d0             * xyz_H2ODelOptDep     )  &
          & / ( xyz_ssa * xyz_TotDelOptDep + 1.0d-100 )



        IDScheme     = IDSchemeShortWave

        QeRatio      = 1.0d0

        SolarFluxTOA = 1.0d0

        call RadiationTwoStreamApp( &
          & QeRatio, xyz_ssa, xyz_af, &
          & SolarFluxTOA, &
          & xy_SurfAlbedo, &
          & IDScheme, &
          & xy_InAngle, &
          & xyr_TotOptDep, &
          & xyr_RadFlux &
          & )
        xyr_RadFlux = xyr_RadFlux * KDFWeight

        xyr_RadSFlux = xyr_RadSFlux + xyr_RadFlux

      end do

    else

      ! Cho and Arking (1981) による短波放射モデル
      ! !!!Short wave radiation model described by Lacis and Hansen (1981)
      !
      call RadiationCA81Flux( &
        & xy_SurfAlbedo, xy_InAngle, xyz_QVap, xyr_Press, xyz_Press, & ! (in )
        & xyr_RadFlux                                                & ! (out)
        & )

      xyr_RadSFlux = xyr_RadSFlux + xyr_RadFlux

    end if






!!$    i = imax / 2
!!$    j = jmax / 2
!!$    write( 93, * ) '# ', xy_SurfAlbedo(i,j)
!!$    do k = 0, kmax
!!$      write( 93, * ) xyr_RadSFlux(i,j,k), xyr_Press(i,j,k)
!!$    end do
!!$    call flush( 93 )
!!$    stop



  end subroutine RadiationDcpamESWV1Flux

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

  subroutine RadiationDcpamESWV1Init

    ! ファイル入出力補助
    ! 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_dcpam_E_SW_V1_nml/ &
      & MeanMolWeight, &
      & O3MolWeight, &
      & FlagRayleighScattering
          !
          ! デフォルト値については初期化手続 "radiation_dcpam_SWEV1#RadiationDcpamSWEV1Init"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "radiation_dcpam_SWEV1#RadiationDcpamSWEV1Init" for the default values.
          !


    ! デフォルト値の設定
    ! Default values settings
    !
    MeanMolWeight          = 28.0d-3

    O3MolWeight            = 48.0d-3

    FlagRayleighScattering = .true.


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

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



    ! Unit of wavelength is changed from micron meter to meter.
    !
    SolarSpecWavelength = SolarSpecWavelength * 1.0d-6

    ! Unit of solar flux is changed from W m-2 micron-1 to W m-2 
    ! (multiply wavelength width)
    !
    SolarSpec           = SolarSpec &
      & * ( SolarSpecWavelengthLimitMicron(:,2) - SolarSpecWavelengthLimitMicron(:,1) )
!!$    !
!!$    ! Value of solar flux integrated from 0.75 to 0.85 micron is scaled to 
!!$    ! that integrated from 0.75 to 0.83 micron, because the wavelength range 
!!$    ! from 0.83 to 0.85 micron is treated by other band. 
!!$    SolarSpec(nswband)  = SolarSpec(nswband) * (0.83d0-0.75d0) / (0.85d0-0.75d0)



    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'MeanMolWeight          = %f', d = (/ MeanMolWeight /) )
    call MessageNotify( 'M', module_name, 'O3MolWeight            = %f', d = (/ O3MolWeight /) )
    call MessageNotify( 'M', module_name, 'FlagRayleighScattering = %b', l = (/ FlagRayleighScattering /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    radiation_dcpam_E_SW_V1_inited = .true.

  end subroutine RadiationDcpamESWV1Init

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

end module radiation_dcpam_E_SW_V1
