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, PI       ! $ \pi $ .
                                  ! Circular constant

    use radiation_two_stream_app, only: RadiationTwoStreamApp

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimeN, EndTime, 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 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 )

!!$    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, xyz_H2ODelOptDep, KDFWeight )


        ! 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, xyr_RadFlux )

      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, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = radiation_dcpam_E_SW_V1_nml, 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