!= dcrtm flux 計算 主プログラム
!
!= dcrtm flux calculation main program
!
! Authors::   Masanori Onishi
! Version::   $Id: main_dcrtm.f90,v 1.00 onishi$
! Tag Name::  $Name: dcrtm-20160908 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2015. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]

module flux
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !

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

  ! netcdf データ入出力
  ! netcdf data input or output
  !
  use gtool_history

  ! 格子点設定
  ! grid point setting
  !
  use gridset, only: GridsetSet, &
    &                imax,       &
    &                jmax,       &
    &                kmax

  ! 2方向近似放射モデル
  ! two stream approximation
  !
  use rad_rte_two_stream_app,                &
    &          only: RadRTETwoStreamAppInit, &
    &                RadRTETwoStreamAppLW,   &
    &                RadRTETwoStreamAppSW

  ! 物理定数設定1
  ! physical constants setteing 1
  !
  use constants0, only: PI, &
    &                   StB

  ! 物理定数設定2
  ! physical constants setteing 2
  !
  use planetconst, only: Grav, &
    &                    AvogadroNum

  ! プランク関数計算
  ! planck function calculation
  !
  use planckf, only: integral_planck

  ! 吸収散乱係数から光学的厚さ, 一次散乱アルベドへの変換
  ! convert absorption and scattering coefficient
  ! to optical depth & single scattering albedo
  !
  use convert, only: KtoTauSSA



  ! 宣言文 ; Declaration statements
  !
  implicit none

  private

  public :: calc_flux

  ! 設定変数, 他
  ! setteing parameter, values
  !
  real(8), allocatable :: r_Press(:)
                          ! 気圧 (半整数レベル) [Pa]
                          ! pressure (half level) [Pa]
  real(8), allocatable :: r_Temp(:)
                          ! 気温 (半整数レベル) [K]
                          ! temperature (half level) [K]
  real(8), allocatable :: rm_MixRatio(:,:)
                          ! 体積混合比 (半整数レベル) [1]
                          ! volume mixing ratio (half level) [1]
  real(8), allocatable :: m_MolWt(:)
                          ! 分子量 [kg mol-1]
                          ! molecular weight [kg mol-1]
  real(8), allocatable :: w_WaveNum(:)
                          ! 波数 [m-1]
                          ! wavenumber [m-1]
  real(8), allocatable :: rw_OptDepTOA(:,:)
                          ! 大気上端からの光学的厚さ (半整数レベル) [1]
                          ! optical depth from top of atmosphere (half level) [1]
  real(8), allocatable :: zw_SingleScatA(:,:)
                          ! 1次散乱アルベド (full level) [1]
                          ! single scattering albedo (full level) [1]
  real(8), allocatable :: rmw_LAbsorpK(:,:,:)
                          ! 線吸収の吸収係数 [m2 kg-1]
                          ! absorption coefficient for line absorption [m2 kg-1]
  real(8), allocatable :: rmw_CAbsorpK(:,:,:)
                          ! 連続吸収の吸収係数 [m2 kg-1]
                          ! absorption coefficient for continuum absotption [m2 kg-1]
  real(8), allocatable :: rmw_RayScatK(:,:,:)
                          ! Rayleigh 散乱の吸収係数 [m2 kg-1]
                          ! scattering coefficient for Rayleigh scattering [m2 kg-1]

  ! 計算される物理量
  ! calculated values
  !
  real(8), allocatable :: rw_RadUwFlux(:,:)
                          ! 惑星大気の射出に由来する上向きフラックス [W m-2]
                          ! upward flux originated from emittance of planetary atmospheire [W m-2]
  real(8), allocatable :: rw_RadDwFlux(:,:)
                          ! 惑星大気の射出に由来する下向きフラックス [W m-2]
                          ! downward flux originated from emittance of planetary atmospheire [W m-2]
  real(8), allocatable :: rw_SolUwFlux(:,:)
                          ! 中心星入射に由来する上向きフラックス [W m-2]
                          ! upward flux originated from emittance of central star [W m-2]
  real(8), allocatable :: rw_SolDwFlux(:,:)
                          ! 中心星入射に由来する下向きフラックス [W m-2]
                          ! downward flux originated from emittance of central star [W m-2]

  integer :: fnml = 12
                          ! 装置番号(name list file)
                          ! device number : namelist file
  integer :: wmax
                          ! 波数レベル最大値
                          ! wavenumber level maximum
  integer :: mmax
                          ! 分子数最大値
                          ! molecular number maxim
  integer :: k_w, k_r, k_m
                          ! 作業変数
                          ! work variables
  real(8) :: RatioSolarConst
                          ! 中心星の黒体放射に対する太陽定数の割合
                          ! ratio of solar constant to black body radiation
  real(8), pointer :: a_Axis(:) => null()
                          ! 配列のサイズを調べるための配列
                          ! dummy array
  integer, pointer :: a_iAxis(:) => null()
                          ! 配列のサイズを調べるための配列
                          ! dummy array

  ! for dcpam5 radiation code
  !
  !integer :: imax
  !integer :: jmax
  !integer :: kmax
  real(8), allocatable :: xyz_SSA          (:, :, :) 
                          ! single scattering albedo
  real(8), allocatable :: xyz_AF           (:, :, :) 
                          ! asynmetory factor
  real(8), allocatable :: xyr_OptDep       (:, :, :) 
                          ! optical depth from top of atmosphere
  real(8), allocatable :: xy_SurfAlbedo    (:, :) 
                          ! surface albedo
  real(8), allocatable :: xyr_PFInted      (:, :, :) 
                          ! planck function integral
  real(8), allocatable :: xy_SurfPFInted   (:, :) 
                          ! planck function in land surface
  real(8), allocatable :: xy_SurfDPFDTInted(:, :) 
                          ! planck function の温度微分
  real(8), allocatable :: xyr_RadUwFlux    (:, :, :) 
                          ! upward flux
  real(8), allocatable :: xyr_RadDwFlux    (:, :, :) 
                          ! downward flux
  real(8), allocatable :: xyra_DelRadUwFlux(:, :, :, :) 
                          ! 
  real(8), allocatable :: xyra_DelRadDwFlux(:, :, :, :)
                          !
  real(8)              :: SolarFluxTOA
                          ! solar flux of top of atmosphere
  real(8), allocatable :: xy_InAngle       (:, :)
                          ! sec (angle of incidence)

contains

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

  subroutine calc_flux( &
    & FlagCKD, FlagConstPF, &
    & Flag_InputFormat, InFileName, OutFileFlux, &
    & SolTemp, InAngle, SurfAlbedo, RatioSolarConst, &
    & NVLev, &
    & BandWNs, BandWNe, &
    & r_ArgPress, r_ArgRadUwFlux, r_ArgRadDwFlux, &
    & FlagOutput &
    & )

    use band_module, only: GetBandNum, GetSurfPFInted, GetPFInted, GetSolarPFInted

    logical, intent(in) :: FlagCKD
    logical, intent(in) :: FlagConstPF
    integer, intent(in) :: Flag_InputFormat
                          ! 入力物理量の指定フラグ
                          ! flag of input value 
                          ! 1: optical depth and single scattering albedo
                          ! 2: absorption, scattering coefficient
    character(*), intent(in) :: InFileName
                          ! 入力 ファイル名
                          ! input file name
    character(*), intent(in) :: OutFileFlux
                          ! 出力ファイル名
                          ! output file name
    real(8), intent(in) :: SolTemp
                          ! 中心星表面温度 [K]
                          ! surface temperature of central star [K]
    real(8), intent(in) :: InAngle
                          ! sec (太陽光入射角) [1]
                          ! sec (angle of incidence) [1]
    real(8), intent(in) :: SurfAlbedo
                          ! 地表面アルベド [1]
                          ! surface albedo [1]
    real(8), intent(in) :: RatioSolarConst

    integer, intent(in) :: NVLev

    real(8), intent(in) :: BandWNs
    real(8), intent(in) :: BandWNe

    real(8), intent(out) :: r_ArgPress(0:NVLev)
    real(8), intent(out) :: r_ArgRadUwFlux(0:NVLev)
    real(8), intent(out) :: r_ArgRadDwFlux(0:NVLev)

    logical, intent(in ), optional :: FlagOutput

    integer              :: NBand
    integer              :: NBin
    integer, allocatable :: a_BandNBin ( : )
    real(8), allocatable :: a_BinWeight( : )
    integer, allocatable :: aa_BandBinIndexBnds( :, : )
    real(8), allocatable :: aa_BandWaveNumBnds( :, : )
    real(8), allocatable :: ra_PFRatio        ( :, : )
    integer              :: iBand
    integer              :: iBin

    integer :: BandWNIndexs
    integer :: BandWNIndexe

    integer :: l


!!$    real(8) :: ResWN


    ! 配列サイズを調べる
    ! check array size
    !
!!$    call HistoryGetPointer(''//InFileName//'', 'r_Press', a_Axis )
!!$    kmax = size( a_Axis ) - 1
!!$    Nullify( a_Axis )


    if ( ( FlagCKD ) .and. ( .not. FlagConstPF ) ) then
      stop 'FlagConstPF has to be .true. when FlagCKD is .true.'
    end if

    kmax = NVLev

    !print *, kmax

    if ( FlagCKD ) then
      call HistoryGetPointer(''//InFileName//'', 'BinNum', a_iAxis, quiet = .true. )
      wmax = size( a_iAxis ) - 1
      Nullify( a_iAxis )
    else
      call HistoryGetPointer(''//InFileName//'', 'w_WaveNum', a_Axis, quiet = .true. )
      wmax = size( a_Axis ) - 1
      Nullify( a_Axis )
    end if

!!$    ResWN = a_Axis(2) - a_Axis(1)


    !print *, wmax
    !print *, ResWN
    !wmax = abs((MaxWN - MinWN)/ResWN)



    call RadRTETwoStreamAppInit

    !== set1: grid
    imax = 1
    jmax = 1
    !kmax = 32

    call GridsetSet(imax, jmax, kmax)


    allocate (xyz_SSA          (0:imax-1, 1:jmax, 1:kmax) )
    allocate (xyz_AF           (0:imax-1, 1:jmax, 1:kmax) )
    allocate (xyr_OptDep       (0:imax-1, 1:jmax, 0:kmax) )
    allocate (xy_SurfAlbedo    (0:imax-1, 1:jmax)         )
    allocate (xyr_PFInted      (0:imax-1, 1:jmax, 0:kmax) )
    allocate (xy_SurfPFInted   (0:imax-1, 1:jmax)         )
    allocate (xy_SurfDPFDTInted(0:imax-1, 1:jmax)         )
    allocate (xyr_RadUwFlux    (0:imax-1, 1:jmax, 0:kmax) )
    allocate (xyr_RadDwFlux    (0:imax-1, 1:jmax, 0:kmax) )
    allocate (xyra_DelRadUwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) )
    allocate (xyra_DelRadDwFlux(0:imax-1, 1:jmax, 0:kmax, 0:1) )
    allocate (xy_InAngle       (0:imax-1, 1:jmax)         ) 


    allocate( r_Press       (0:kmax        ) )
    allocate( r_Temp        (0:kmax        ) )
    allocate( w_WaveNum     (        0:wmax) )
    allocate( rw_RadUwFlux  (0:kmax, 0:wmax) )
    allocate( rw_RadDwFlux  (0:kmax, 0:wmax) )
    allocate( rw_SolUwFlux  (0:kmax, 0:wmax) )
    allocate( rw_SolDwFlux  (0:kmax, 0:wmax) )
    allocate( rw_OptDepTOA  (0:kmax, 0:wmax) )
    allocate( zw_SingleScatA(1:kmax, 0:wmax) )


    if ( FlagCKD ) then
      w_WaveNum    = -1.0d10
    else
      call HistoryGet(''//InFileName//'', 'w_WaveNum', w_WaveNum, quiet = .true.)
      do k_w = 0, wmax
        if( w_WaveNum(k_w) == 0.0_8 ) then
          write( 6, * ) 'STOP: Wavenumber is zero at k = ', k_w, '.'
          stop
        end if
      end do
    end if


    ! 鉛直プロファイルの読み込み
    ! read profile data
    !
    if(Flag_InputFormat == 1) then

      call HistoryGet(''//InFileName//'', 'r_Press', r_Press, quiet = .true.)
      call HistoryGet(''//InFileName//'', 'r_Temp', r_Temp, quiet = .true.)
      call HistoryGet(''//InFileName//'', 'rw_OptDep', rw_OptDepTOA, quiet = .true.)
      call HistoryGet(''//InFileName//'', 'zw_SingleScatA', zw_SingleScatA, quiet = .true.)

    elseif(Flag_InputFormat == 2) then

      ! 配列サイズを調べる
      ! check array size
      !
      call HistoryGetPointer(''//InFileName//'', 'm_MolWt', a_Axis, quiet = .true.)

      mmax = size( a_Axis ) - 1

      Nullify( a_Axis )

      !print *, mmax

      allocate( m_MolWt       (        0:mmax        ) )
      allocate( rm_MixRatio   (0:kmax, 0:mmax        ) )
      allocate( rmw_LAbsorpK  (0:kmax, 0:mmax, 0:wmax) )
      allocate( rmw_CAbsorpK  (0:kmax, 0:mmax, 0:wmax) )
      allocate( rmw_RayScatK  (0:kmax, 0:mmax, 0:wmax) )

      call HistoryGet(''//InFileName//'', 'r_Press', r_Press, quiet = .true.)
      call HistoryGet(''//InFileName//'', 'r_Temp', r_Temp, quiet = .true.)
      call HistoryGet(''//InFileName//'', 'm_MolWt', m_MolWt, quiet = .true.)
      call HistoryGet(''//InFileName//'', 'rm_MixRatio', rm_MixRatio, quiet = .true.)
      if ( FlagCKD ) then
        call HistoryGet(''//InFileName//'', 'rma_TAbsorpK', rmw_LAbsorpK, quiet = .true.)
        rmw_CAbsorpK = 0.0d0
        rmw_RayScatK = 0.0d0
      else
!!$        call HistoryGet(''//InFileName//'', 'rmw_LAbsorpK', rmw_LAbsorpK, quiet = .true.)
!!$        call HistoryGet(''//InFileName//'', 'rmw_CAbsorpK', rmw_CAbsorpK, quiet = .true.)
!!$        call HistoryGet(''//InFileName//'', 'rmw_RayScatK', rmw_RayScatK, quiet = .true.)
        call HistoryGet(''//InFileName//'', 'rmw_TAbsorpK', rmw_LAbsorpK, quiet = .true.)
        rmw_CAbsorpK = 0.0d0
        rmw_RayScatK = 0.0d0
      end if

      call KtoTauSSA (                                          & ! 2016/08/05
        & kmax, mmax, wmax,                                     &
        & r_Press, r_Temp, rm_MixRatio(0:kmax,0:mmax), m_MolWt, &
        & rmw_RayScatK, rmw_LAbsorpK, rmw_CAbsorpK,             &
        & rw_OptDepTOA, zw_SingleScatA                          )

    else

      print *, 'invalid value (Flag_InputFormat):', Flag_InputFormat
      print *, 'stop'
      stop
 
    end if


    if ( FlagCKD ) then
      call HistoryGetPointer(''//InFileName//'', 'BandNum', a_iAxis, quiet = .true. )
      NBand = size( a_iAxis )
      Nullify( a_iAxis )
      call HistoryGetPointer(''//InFileName//'', 'BinNum', a_iAxis, quiet = .true. )
      NBin = size( a_iAxis )
      Nullify( a_iAxis )

      allocate( a_BandNBin ( NBand ) )
      allocate( a_BinWeight( NBin  ) )
      allocate( aa_BandBinIndexBnds( 2, NBand ) )
      allocate( aa_BandWaveNumBnds ( 2, NBand ) )

      call HistoryGet(''//InFileName//'', 'BandNBin' , a_BandNBin , quiet = .true.)
      call HistoryGet(''//InFileName//'', 'BinWeight', a_BinWeight, quiet = .true.)
      call HistoryGet(''//InFileName//'', 'BandBinIndexBnds', aa_BandBinIndexBnds, quiet = .true.)
      call HistoryGet(''//InFileName//'', 'BandWaveNumBnds', aa_BandWaveNumBnds, quiet = .true.)
    end if
    if ( FlagCKD .or. FlagConstPF ) then
      call HistoryGetPointer(''//InFileName//'', 'BinNum', a_iAxis, quiet = .true. )
      NBin = size( a_iAxis )
      Nullify( a_iAxis )
      allocate( ra_PFRatio         ( 0:kmax, NBin ) )
      call HistoryGet(''//InFileName//'', 'ra_PFRatio'     , ra_PFRatio, quiet = .true.)
    end if


    ! フラックスの計算
    ! calculate flux
    !

    ! dcpam5
    !xyz_SSA          (:, :, :) 
    xyz_AF           (:, :, :) = 0.0_8
    !xyr_OptDep       (:, :, :) 
    xy_SurfAlbedo    (:, :)    = SurfAlbedo
    !xyr_PFInted      (:, :, :)
    !xy_SurfPFInted   (:, :)
    xy_SurfDPFDTInted(:, :)    = 0.0_8
    xy_InAngle       (:, :)    = InAngle




    ! dcpam5 radiation code
    !
    do k_w = 0, wmax

!      if( w_WaveNum(k_w) <= 300000.0_8 ) then
!        xy_SurfAlbedo(:,:) = 0.0_8
!      else
      xy_SurfAlbedo(:,:) = SurfAlbedo
!      endif

      do k_r = 1, kmax
        xyz_SSA(:,:,k_r) = zw_SingleScatA(k_r,k_w)

        if(xyz_SSA(0,1,k_r) >= 1.0_8) then ! 2014/04/21
          xyz_SSA(:,:,k_r) = 1.0_8 - 10.0**(-16) !2014/04/22 !0.999999_8
        end if

      end do

      do k_r = 0, kmax
        xyr_OptDep(:,:,k_r) = rw_OptDepTOA(k_r, k_w)
      end do

      if ( FlagCKD ) then
        iBin = k_w + 1
        !
        search_bandnum : do iBand = 1, NBand
          if ( ( aa_BandBinIndexBnds(1,iBand) <= iBin ) .and. &
            &  ( iBin <= aa_BandBinIndexBnds(2,iBand) ) ) exit search_bandnum
        end do search_bandnum
        if ( ( iBand < 1 ) .or. ( iBand > NBand ) ) then
          write( 6, * ) 'iBand is inappropriate: ', iBand, NBand
          stop
        end if
        !
        xy_SurfPFInted(:, :) = GetSurfPFInted( iBand ) &
          & * ra_PFRatio( 0, iBin )
        do k_r = 0, kmax
          xyr_PFInted(:, :,k_r) = GetPFInted( k_r, iBand ) &
            & * ra_PFRatio( k_r, iBin )
        end do
        ! This may need to correct.
!!$        SolarFluxTOA = GetSolarPFInted( iBand ) &
!!$          & * ra_PFRatio( k_r, iBin )
        SolarFluxTOA = GetSolarPFInted( iBand ) &
          & * 1.0d0

      else
        if ( FlagConstPF ) then
          iBin = k_w + 1
          xy_SurfPFInted(:, :) = GetSurfPFInted( GetBandNum(w_WaveNum(k_w)) ) &
            & * ra_PFRatio( 0, iBin )
          do k_r = 0, kmax
            xyr_PFInted(:, :,k_r) = GetPFInted( k_r, GetBandNum(w_WaveNum(k_w)) ) &
              & * ra_PFRatio( k_r, iBin )
          end do
          SolarFluxTOA = GetSolarPFInted( GetBandNum(w_WaveNum(k_w)) ) &
            & * 1.0d0
        else
          xy_SurfPFInted(:, :) = PI * integral_planck(w_WaveNum(k_w) - 0.5_8, w_WaveNum(k_w) + 0.5_8, r_Temp(0))
          do k_r = 0, kmax
            xyr_PFInted(:, :,k_r) = PI * integral_planck(w_WaveNum(k_w) - 0.5_8, w_WaveNum(k_w) + 0.5_8, r_Temp(k_r))
          end do
          SolarFluxTOA = PI * RatioSolarConst * integral_planck(w_WaveNum(k_w) - 0.5_8, w_WaveNum(k_w) + 0.5_8, SolTemp)
        end if
      end if


        ! dcpam5 radiation code

      call RadRTETwoStreamAppSW(            &
        & xyz_SSA, xyz_AF,                  & ! (in)
        & xyr_OptDep,                       & ! (in)
        & xy_SurfAlbedo,                    & ! (in)
        & SolarFluxTOA, xy_InAngle,         & ! (in)
        & xyr_RadUwFlux, xyr_RadDwFlux      & ! (out)
        & )

      rw_SolUwFlux(0:kmax,k_w) = xyr_RadUwFlux(0,1,0:kmax)
      rw_SolDwFlux(0:kmax,k_w) = xyr_RadDwFlux(0,1,0:kmax)

      call RadRTETwoStreamAppLW(                           &
        & xyz_SSA, xyz_AF,                                 & ! (in)
        & xyr_OptDep,                                      & ! (in)
        & xy_SurfAlbedo,                                   & ! (in)
        & xyr_PFInted, xy_SurfPFInted, xy_SurfDPFDTInted,  & ! (in)
        & xyr_RadUwFlux, xyr_RadDwFlux,                    & ! (out)
        & xyra_DelRadUwFlux, xyra_DelRadDwFlux             & ! (out)
        & )

      rw_RadUwFlux(0:kmax,k_w) = xyr_RadUwFlux(0,1,0:kmax)
      rw_RadDwFlux(0:kmax,k_w) = xyr_RadDwFlux(0,1,0:kmax)

    end do

    ! 出力ファイルの設定
    ! setting of output file
    !
    call HistoryCreate(                                    &
      & quiet = .true., &
      & file=''//OutFileFlux//'', &
      & title='Flux; '//OutFileFlux//'', &
      & source='main_flux.f90;'//InFileName//'', &
      & institution='DCRTM',       &
      & dims=(/'r_Press  ', 'w_WaveNum'/), dimsizes=(/kmax+1,wmax+1/),      &
      & longnames=(/'pressure  ', 'wavenumber'/),       &
      & units=(/'Pa ', 'm-1'/)                    )

    call HistoryPut('r_Press  ',r_Press)
    call HistoryPut('w_WaveNum',w_WaveNum)

    call HistoryAddVariable( &
      & varname='r_Temp', dims=(/'r_Press'/), &
      & longname='temperature', units='K', xtype='double')
    call HistoryPut('r_Temp',r_Temp)

    call HistoryAddVariable( &
      & varname='rw_RadFluxUp', dims=(/'r_Press  ', 'w_WaveNum'/), &
      & longname='upward flux(planet)', units='W m-2 m', xtype='double')
    call HistoryPut('rw_RadFluxUp',rw_RadUwFlux)

    call HistoryAddVariable( &
      & varname='rw_RadFluxDn', dims=(/'r_Press  ', 'w_WaveNum'/), &
      & longname='downward flux(planet)', units='W m-2 m', xtype='double')
    call HistoryPut('rw_RadFluxDn',rw_RadDwFlux)

    call HistoryAddVariable( &
      & varname='rw_SolFluxUp', dims=(/'r_Press  ', 'w_WaveNum'/), &
      & longname='upward flux(solar)', units='W m-2 m', xtype='double')
    call HistoryPut('rw_SolFluxUp',rw_SolUwFlux)

    call HistoryAddVariable( &
      & varname='rw_SolFluxDn', dims=(/'r_Press  ', 'w_WaveNum'/), &
      & longname='downward flux(solar)', units='W m-2 m', xtype='double')
    call HistoryPut('rw_SolFluxDn',rw_SolDwFlux)

    call HistoryClose( quiet = .true. )


    r_ArgPress = r_Press

    if ( FlagCKD ) then
      r_ArgRadUwFlux = 0.0d0
      r_ArgRadDwFlux = 0.0d0
      do iBand = 1, NBand
        if ( ( BandWNs <= aa_BandWaveNumBnds(1,iBand) ) .and. &
          &  ( aa_BandWaveNumBnds(2,iBand) <= BandWNe ) ) then
          do iBin = aa_BandBinIndexBnds(1,iBand), aa_BandBinIndexBnds(2,iBand)
            k_w = iBin - 1
!!$          r_ArgRadUwFlux = r_ArgRadUwFlux &
!!$            & + rw_RadUwFlux(:,k_w) * a_BinWeight(iBin) &
!!$            & * ( aa_BandWaveNumBnds(2,iBand) - aa_BandWaveNumBnds(1,iBand) )
!!$          r_ArgRadDwFlux = r_ArgRadDwFlux &
!!$            & + rw_RadDwFlux(:,k_w) * a_BinWeight(iBin) &
!!$            & * ( aa_BandWaveNumBnds(2,iBand) - aa_BandWaveNumBnds(1,iBand) )
            r_ArgRadUwFlux = r_ArgRadUwFlux &
              & + rw_RadUwFlux(:,k_w) * a_BinWeight(iBin)
            r_ArgRadDwFlux = r_ArgRadDwFlux &
              & + rw_RadDwFlux(:,k_w) * a_BinWeight(iBin)
          end do
        end if
      end do

    else

      ! Look for bound index
      loop_search_bound_index_1 : do l = 0+1, wmax
        if ( w_WaveNum(l) > BandWNs ) exit loop_search_bound_index_1
      end do loop_search_bound_index_1
      BandWNIndexs = l-1
      !
      loop_search_bound_index : do l = BandWNIndexs+1, wmax
        if ( w_WaveNum(l) > BandWNe ) exit loop_search_bound_index
      end do loop_search_bound_index
      BandWNIndexe = l-1
      !
      r_ArgRadUwFlux = 0.0d0
      do k_w = BandWNIndexs, BandWNIndexe
        r_ArgRadUwFlux = r_ArgRadUwFlux &
          & + rw_RadUwFlux(:,k_w) * ( w_WaveNum(1) - w_WaveNum(0) )
      end do
      r_ArgRadDwFlux = 0.0d0
      do k_w = BandWNIndexs, BandWNIndexe
        r_ArgRadDwFlux = r_ArgRadDwFlux &
          & + rw_RadDwFlux(:,k_w) * ( w_WaveNum(1) - w_WaveNum(0) )
      end do

    end if


    if ( FlagCKD ) then
      deallocate( a_BandNBin )
      deallocate( a_BinWeight )
      deallocate( aa_BandBinIndexBnds )
      deallocate( aa_BandWaveNumBnds  )
      deallocate( ra_PFRatio          )
    end if


    if(Flag_InputFormat == 1) then

    elseif(Flag_InputFormat == 2) then

      deallocate( m_MolWt      )
      deallocate( rm_MixRatio  )
      deallocate( rmw_LAbsorpK )
      deallocate( rmw_CAbsorpK )
      deallocate( rmw_RayScatK )
    end if

    deallocate( r_Press           )
    deallocate( r_Temp            )
    deallocate( w_WaveNum         )
    deallocate( rw_RadUwFlux      )
    deallocate( rw_RadDwFlux      )
    deallocate( rw_SolUwFlux      )
    deallocate( rw_SolDwFlux      )
    deallocate( rw_OptDepTOA      )
    deallocate( zw_SingleScatA    )

    deallocate (xyz_SSA           )
    deallocate (xyz_AF            )
    deallocate (xyr_OptDep        )
    deallocate (xy_SurfAlbedo     )
    deallocate (xyr_PFInted       )
    deallocate (xy_SurfPFInted    )
    deallocate (xy_SurfDPFDTInted )
    deallocate (xyr_RadUwFlux     )
    deallocate (xyr_RadDwFlux     )
    deallocate (xyra_DelRadUwFlux )
    deallocate (xyra_DelRadDwFlux )
    deallocate (xy_InAngle        )


  end subroutine calc_flux

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

end module flux
