module rad_flux

  use vtype_module

  implicit none

  private

  public :: RadFluxCalcFlux
  public :: RadFluxInit
  public :: RadFluxFinalize
  public :: RadFluxOutputInit
  public :: RadFluxOutputFinalize


  logical, parameter :: FlagRayleigh = .true.
!!$  logical, parameter :: FlagRayleigh = .false.

  integer, save      :: IDMethod
  integer, parameter :: IDMethodKD         = 0
  integer, parameter :: IDMethodKDDevelop  = 1
  integer, parameter :: IDMethodLBL        = 2
  integer, parameter :: IDMethodLBLConstPF = 3
  integer, parameter :: IDMethodNoUse      = -99

  real(8), save :: Grav

  integer              , save :: NBand
  real(DP), allocatable, save :: aa_BandWNBnds      (:,:)
  integer , allocatable, save :: aa_BandBinIndexBnds(:,:)


  logical              , save :: FlagOutputInited = .false.
  integer              , save :: NcIDFlux

!!$  real(DP)             , save :: FluxRatTOAtoStelSfc

  integer              , save :: iWaveNumMaxforAlbedoLW


  character(128), save :: ModuleName = 'flux'

  logical, save :: FlagInited       = .false.


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

contains

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

  subroutine RadFluxCalcFlux( &
    & kmax, NMolNum, &
    & m_MolNum, &
    & r_Press, &
    & SurfTemp, r_Temp, rm_VMR, &
    & NPtcl, a_PtclName, za_PtclEffRad, za_DelPtclNum, &
    & AlbedoSW, AlbedoLW, SZADeg, &
    & StrFluxTOA, &
    & r_MeanMolWt, &
    & r_RadUwFlux, r_RadDwFlux, &
    & ra_DelRadUwFlux, ra_DelRadDwFlux, &
    & r_RadUwFluxLW, r_RadDwFluxLW, &
    & r_RadUwFluxSW, r_RadDwFluxSW, &
    & ra_DelRadUwFluxLW, ra_DelRadDwFluxLW, &
    & ra_DelRadUwFluxSW, ra_DelRadDwFluxSW, &
    & WNIntegStart, WNIntegEnd &
    & )

    use ni3_module

    use constants0, only : &
      & PI, &
      & AtomMassConst, &
      & AvogConst

    use rad_kd_utils, only : InqPFIntGQNum

    use rad_rte_two_stream_app_1d_wrapper, only : &
      & RadRTETwoStreamApp1DLWWrapper, &
      & RadRTETwoStreamApp1DSWWrapper, &
      & RadRTETwoStreamApp1DWrapper

    use opt_prop, only : &
      & OptPropGetNMol, &
      & OptPropGetNWN, &
      & OptPropGetNBand, &
      & OptPropGetDelWN, &
      & OptPropGetWN, &
      & OptPropGetBandNum, &
      & OptPropCalcIndexAndWeight, &
      & OptPropGetAbsCoefProf, &
      & OptPropGetAbsCoefProfForKD, &
      & OptPropGetRayScatCoef, &
      & OptPropGetPtclParam, &
      & OptPropGetPFInted, &
      & OptPropGetStrPFInted, &
      & OptPropGetStrBandAvePF, &
      & OptPropGetPFIntedRatio, &
      & OptPropGetPFIntedRatioForKD, &
      & OptPropGetStrPFIntedRatio, &
      & OptPropPFTblGetPFDPFDT, &
      & OptPropPFTblChkBandIntPFRatio, &
      & OptPropPFTblChkBandIntStrFluxRatio

      ! [DEBUG] Lines for input optical depth
!!$    use opt_prop_lbl, only : &
!!$      & OptPropLBLGetInputedOptDepProf

    use planck_func_wrapper, only : &
      & Integ_PF_GQ_Array1D   , &
      & Integ_PF_GQ_Array0D   , &
      & Integ_DPFDT_GQ_Array0D


    integer , intent(in ) :: kmax
    integer , intent(in ) :: NMolNum
    integer , intent(in ) :: m_MolNum    (1:NMolNum)
    real(DP), intent(in ) :: r_Press     (0:kmax)
    real(DP), intent(in ) :: SurfTemp
    real(DP), intent(in ) :: r_Temp      (0:kmax)
    real(DP), intent(in ) :: rm_VMR      (0:kmax,1:NMolNum)
    integer     , intent(in ) :: NPtcl
    character(*), intent(in ) :: a_PtclName   (NPtcl)
    real(8)     , intent(in ) :: za_PtclEffRad(1:kmax,NPtcl)
    real(8)     , intent(in ) :: za_DelPtclNum(1:kmax,NPtcl)
    real(DP), intent(in ) :: AlbedoSW
    real(DP), intent(in ) :: AlbedoLW
    real(DP), intent(in ) :: SZADeg
    real(DP), intent(in ) :: StrFluxTOA
    real(DP), intent(in ) :: r_MeanMolWt    (0:kmax)
    real(DP), intent(out) :: r_RadUwFlux    (0:kmax)
    real(DP), intent(out) :: r_RadDwFlux    (0:kmax)
    real(DP), intent(out) :: ra_DelRadUwFlux(0:kmax,-1:1)
    real(DP), intent(out) :: ra_DelRadDwFlux(0:kmax,-1:1)
    real(DP), intent(out) :: r_RadUwFluxLW    (0:kmax)
    real(DP), intent(out) :: r_RadDwFluxLW    (0:kmax)
    real(DP), intent(out) :: r_RadUwFluxSW    (0:kmax)
    real(DP), intent(out) :: r_RadDwFluxSW    (0:kmax)
    real(DP), intent(out) :: ra_DelRadUwFluxLW(0:kmax,-1:1)
    real(DP), intent(out) :: ra_DelRadDwFluxLW(0:kmax,-1:1)
    real(DP), intent(out) :: ra_DelRadUwFluxSW(0:kmax,-1:1)
    real(DP), intent(out) :: ra_DelRadDwFluxSW(0:kmax,-1:1)
    real(DP), intent(in ) :: WNIntegStart
    real(DP), intent(in ) :: WNIntegEnd


    ! local variables
    !
    real(DP) :: z_Press     (1:kmax)
    real(DP) :: z_Temp      (1:kmax)
    real(DP) :: zm_VMR      (1:kmax,1:NMolNum)

    real(DP), allocatable :: a_BandAveSurfPF   (:)
    real(DP), allocatable :: ra_BandAvePF      (:,:)
    real(DP), allocatable :: a_BandAveSurfDPFDT(:)
    real(DP), allocatable :: a_BandAveStrFluxTOA(:)

    integer :: NMolTmp

    integer :: NWaveNum

    integer :: iBand
    integer :: iBandStart
    integer :: iBandEnd

    integer :: iWaveNum
    integer :: iWaveNumStart
    integer :: iWaveNumEnd
    real(8) :: DelWaveNum

!!$    real(8), allocatable :: zm_AbsCoefPerAtmMol(:,:)
    real(DP) :: zm_AbsCoefPerAtmMol(1:kmax, NMolNum)

    real(DP) :: RayScatCoefNonRadAct
    real(DP) :: m_RayScatCoef(NMolNum)

    real(DP) :: GasAbsCoef
    real(DP) :: RayScatCoef
    real(DP) :: DelGasAbsOptDep
    real(DP) :: DelRayScatOptDep
    real(DP) :: a_DelPtclExtOptDep(NPtcl)

    real(dp) :: za_PtclQExt(1:kmax, NPtcl)
    real(dp) :: za_PtclSSA (1:kmax, NPtcl)
    real(dp) :: za_PtclAF  (1:kmax, NPtcl)

    real(8) :: z_DelAtmMols   (1:kmax)
    real(8) :: DelOptDep

    real(DP) :: z_SSA          (1:kmax)
    real(DP) :: z_AF           (1:kmax)
    real(DP) :: r_OptDep       (0:kmax)
    real(DP) :: z_DelOptDep    (1:kmax)
    real(DP) :: r_PFInted      (0:kmax)
    real(DP) :: SurfPFInted
    real(DP) :: SurfDPFDTInted

    real(DP) :: r_RadUwFluxEachLW    (0:kmax)
    real(DP) :: r_RadDwFluxEachLW    (0:kmax)
    real(DP) :: ra_DelRadUwFluxEachLW(0:kmax,-1:1)
    real(DP) :: ra_DelRadDwFluxEachLW(0:kmax,-1:1)
    real(DP) :: r_RadUwFluxEachSW    (0:kmax)
    real(DP) :: r_RadDwFluxEachSW    (0:kmax)
    real(DP) :: ra_DelRadUwFluxEachSW(0:kmax,-1:1)
    real(DP) :: ra_DelRadDwFluxEachSW(0:kmax,-1:1)

    real(DP) :: r_RadUwFluxEach    (0:kmax)
    real(DP) :: r_RadDwFluxEach    (0:kmax)

    real(DP) :: z_DRadFluxDPressEach(1:kmax)

    real(DP) :: r_PFIntedRatio (0:kmax)
    real(DP) :: SurfPFIntedRatio
    real(DP) :: StrPFRatio

    real(DP) :: StrFluxTOAEach
    real(DP) :: InAngle

    real(DP) :: SurfAlbedo

!!$    real(DP) :: TempMin
    real(DP) :: VMRNonRadAct

    real(DP) :: WNs
    real(DP) :: WNe
    integer  :: GQNum

    real(dp) :: ThresholdBandIntPFRatio      = 1.0d-3
    real(dp) :: ThresholdBandIntStrFluxRatio = 1.0d-3
    logical  :: FlagBandIntPFRatio
    logical  :: FlagBandIntStrFluxRatio

    ! Variables for index and weight for KD
    integer, parameter :: NIntpol = 2
    integer            :: z_KDTblPressIndex  (1:kmax)
    integer            :: z_KDTblTempIndex   (1:kmax)
    integer            :: z_KDTblVMRMol1Index(1:kmax)
    integer            :: z_KDTblVMRMol2Index(1:kmax)
    integer            :: z_KDTblVMRMol3Index(1:kmax)
    integer            :: z_NIntpol1         (1:kmax)
    integer            :: z_NIntpol3         (1:kmax)
    integer            :: z_NIntpol4         (1:kmax)
    integer            :: z_NIntpol5         (1:kmax)
    real(DP)           :: ptuvwz_Weight(NIntpol,NIntpol,NIntpol,NIntpol,NIntpol,1:kmax)
    integer            :: r_KDTblPressIndex  (0:kmax)
    integer            :: r_KDTblTempIndex   (0:kmax)
    integer            :: r_KDTblVMRMol1Index(0:kmax)
    integer            :: r_KDTblVMRMol2Index(0:kmax)
    integer            :: r_KDTblVMRMol3Index(0:kmax)
    integer            :: r_NIntpol1         (0:kmax)
    integer            :: r_NIntpol3         (0:kmax)
    integer            :: r_NIntpol4         (0:kmax)
    integer            :: r_NIntpol5         (0:kmax)
    real(DP)           :: ptuvwr_Weight(NIntpol,NIntpol,NIntpol,NIntpol,NIntpol,0:kmax)
    integer            :: r_KDTblPressIndexSfc  (1)
    integer            :: r_KDTblTempIndexSfc   (1)
    integer            :: r_KDTblVMRMol1IndexSfc(1)
    integer            :: r_KDTblVMRMol2IndexSfc(1)
    integer            :: r_KDTblVMRMol3IndexSfc(1)
    integer            :: r_NIntpol1Sfc         (1)
    integer            :: r_NIntpol3Sfc         (1)
    integer            :: r_NIntpol4Sfc         (1)
    integer            :: r_NIntpol5SFc         (1)
    real(DP)           :: ptuvw_WeightSfc(NIntpol,NIntpol,NIntpol,NIntpol,NIntpol)

    integer :: k
    integer :: m
    integer :: n


    ! [DEBUG] Variables for input optical depth
!!$    real(DP) :: z_InputedDelOptDep(1:kmax)



    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if

    do k = 1, kmax
      z_Press(k) = ( r_Press(k-1) + r_Press(k) ) / 2.0d0
      z_Temp (k) = ( r_Temp (k-1) + r_Temp(k) ) / 2.0d0
    end do
    do m = 1, NMolNum
      do k = 1, kmax
        zm_VMR(k,m) = ( rm_VMR(k-1,m) + rm_VMR(k,m) ) / 2.0d0
      end do
    end do


    select case ( IDMethod )
!!$    case ( IDMethodLBLConstPF, IDMethodKDDevelop, IDMethodKD )
    case ( IDMethodLBLConstPF, IDMethodKDDevelop )

      allocate( a_BandAveSurfPF    (NBand) )
      allocate( ra_BandAvePF       (0:kmax,NBand) )
      allocate( a_BandAveSurfDPFDT (NBand) )
      allocate( a_BandAveStrFluxTOA(NBand) )

      do n = 1, NBand
        WNs = aa_BandWNBnds(1,n)
        WNe = aa_BandWNBnds(2,n)
        !
        GQNum = InqPFIntGQNum( WNe - WNs )
        !
        call Integ_PF_GQ_Array0D( &
          & WNs, WNe, GQNum, &
          & SurfTemp, &
          & a_BandAveSurfPF(n) &
          & )
        a_BandAveSurfPF(n) = PI * a_BandAveSurfPF(n) / ( WNe - WNs )
        call Integ_PF_GQ_Array1D( &
          & WNs, WNe, GQNum, &
          & 0, kmax, &
          & r_Temp, &
          & ra_BandAvePF(:,n) &
          & )
        ra_BandAvePF(:,n) = PI * ra_BandAvePF(:,n) / ( WNe - WNs )
        call Integ_DPFDT_GQ_Array0D( &
          & WNs, WNe, GQNum,                 & ! (in )
          & SurfTemp,                & ! (in )
          & a_BandAveSurfDPFDT(n) &
          & )
        a_BandAveSurfDPFDT(n) = PI * a_BandAveSurfDPFDT(n) / ( WNe - WNs )
      end do

      do n = 1, NBand
        call OptPropGetStrBandAvePF( &
          & n, StrFluxTOA, &
          & a_BandAveStrFluxTOA(n) &
          & )
      end do

!!$    case ( IDMethodNoUse )
    case ( IDMethodKD )

      allocate( a_BandAveSurfPF    (NBand) )
      allocate( ra_BandAvePF       (0:kmax,NBand) )
      allocate( a_BandAveSurfDPFDT (NBand) )
      allocate( a_BandAveStrFluxTOA(NBand) )

      do n = 1, NBand
        call OptPropPFTblGetPFDPFDT( &
          & n, SurfTemp, &
          & a_BandAveSurfPF(n), &
          & .false. &
          & )
        call OptPropPFTblGetPFDPFDT( &
          & n, kmax+1, r_Temp, &
          & ra_BandAvePF(:,n), &
          & .false. &
          & )

        call OptPropPFTblGetPFDPFDT( &
          & n, SurfTemp, &
          & a_BandAveSurfDPFDT(n), &
          & .true. &
          & )
      end do

      do n = 1, NBand
        call OptPropGetStrBandAvePF( &
          & n, StrFluxTOA, &
          & a_BandAveStrFluxTOA(n) &
          & )
      end do

    end select



    do k = 1, kmax
!!$      z_DelAtmMols(k) = &
!!$        &   ( r_Press(k-1) - r_Press(k) ) / Grav
      z_DelAtmMols(k) = &
        &   ( r_Press(k-1) - r_Press(k) ) / Grav &
        & / ( ( r_MeanMolWt(k-1) + r_MeanMolWt(k) ) / 2.0d0 &
        &     / AvogConst )
    end do


    NMolTmp    = OptPropGetNMol()
    if ( NMolTmp /= NMolNum ) then
      write( 6, * ) 'NMolNum in argument does not agree with NMol in data.'
      write( 6, * ) NMolNum, NMolTmp
      stop
    end if
    NWaveNum   = OptPropGetNWN()

!!$    allocate( zm_AbsCoefPerAtmMol(1:kmax, NMolNum) )


!!$    r_RadUwFlux     = 0.0d0
!!$    r_RadDwFlux     = 0.0d0
!!$    ra_DelRadUwFlux = 0.0d0
!!$    ra_DelRadDwFlux = 0.0d0

    r_RadUwFluxLW     = 0.0d0
    r_RadDwFluxLW     = 0.0d0
    ra_DelRadUwFluxLW = 0.0d0
    ra_DelRadDwFluxLW = 0.0d0
    r_RadUwFluxSW     = 0.0d0
    r_RadDwFluxSW     = 0.0d0
    ra_DelRadUwFluxSW = 0.0d0
    ra_DelRadDwFluxSW = 0.0d0

    select case ( IDMethod )
    case ( IDMethodKDDevelop, IDMethodKD )
      if ( ( WNIntegStart < 0.0d0 ) .or. ( WNIntegEnd < 0.0d0 ) ) then
        iWaveNumStart = 1
        iWaveNumEnd   = NWaveNum
      else
        do iBand = 1, NBand
          if ( aa_BandWNBnds(1,iBand) >= WNIntegStart ) exit
        end do
        iBandStart = iBand
        iBandStart = min( iBandStart, NBand )
        do iBand = iBandStart, NBand
          if ( aa_BandWNBnds(2,iBand) >= WNIntegEnd ) exit
        end do
        iBandEnd = min( iBand, NBand )
        iWaveNumStart = aa_BandBinIndexBnds(1,iBandStart)
        iWaveNumEnd   = aa_BandBinIndexBnds(2,iBandEnd  )
      end if
    case ( IDMethodLBL, IDMethodLBLConstPF )
      if ( ( WNIntegStart < 0.0d0 ) .or. ( WNIntegEnd < 0.0d0 ) ) then
        iWaveNumStart = 1
        iWaveNumEnd   = NWaveNum
      else
        loop_search_wnstart: do iWaveNum = 1, NWaveNum
          if ( OptPropGetWN( iWaveNum ) >= WNIntegStart ) exit
        end do loop_search_wnstart
        iWaveNumStart = iWaveNum
        iWaveNumStart = min( iWaveNumStart, NWaveNum )
        loop_search_wnend: do iWaveNum = iWaveNumStart, NWaveNum
          if ( OptPropGetWN( iWaveNum ) > WNIntegEnd ) exit
        end do loop_search_wnend
        iWaveNumEnd = iWaveNum - 1
      end if
    case default
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Unexpected case default.'
      stop
    end select



!!$    do iWaveNum = 1, NWaveNum
    do iWaveNum = iWaveNumStart, iWaveNumEnd

      DelWaveNum = OptPropGetDelWN( iWaveNum )


      select case ( IDMethod )
      case ( IDMethodKD )
        if ( iWaveNum == iWaveNumStart ) then
          ! NOTE:
          !   r_KDTbl*IndexSfc, r_NIntpol?Sfc, and ptuwv_WeightSfc will be 
          !   overwritten below.
          call OptPropCalcIndexAndWeight( &
            & kmax, NMolNum, m_MolNum, z_Press, z_Temp, zm_VMR, SurfTemp, &
            & NIntpol, &
            & z_KDTblPressIndex, z_KDTblTempIndex, &
            & z_KDTblVMRMol1Index, z_KDTblVMRMol2Index, z_KDTblVMRMol3Index, &
            & z_NIntpol1, z_NIntpol3, z_NIntpol4, z_NIntpol5, &
            & ptuvwz_Weight, &
            & r_KDTblPressIndexSfc, r_KDTblTempIndexSfc, &
            & r_KDTblVMRMol1IndexSfc, r_KDTblVMRMol2IndexSfc, r_KDTblVMRMol3IndexSfc, &
            & r_NIntpol1Sfc, r_NIntpol3Sfc, r_NIntpol4Sfc, r_NIntpol5SFc, &
            & ptuvw_WeightSfc &
            & )
          call OptPropCalcIndexAndWeight( &
            & kmax+1, NMolNum, m_MolNum, r_Press, r_Temp, rm_VMR, SurfTemp, &
            & NIntpol, &
            & r_KDTblPressIndex, r_KDTblTempIndex, &
            & r_KDTblVMRMol1Index, r_KDTblVMRMol2Index, r_KDTblVMRMol3Index, &
            & r_NIntpol1, r_NIntpol3, r_NIntpol4, r_NIntpol5, &
            & ptuvwr_Weight, &
            & r_KDTblPressIndexSfc, r_KDTblTempIndexSfc, &
            & r_KDTblVMRMol1IndexSfc, r_KDTblVMRMol2IndexSfc, r_KDTblVMRMol3IndexSfc, &
            & r_NIntpol1Sfc, r_NIntpol3Sfc, r_NIntpol4Sfc, r_NIntpol5SFc, &
            & ptuvw_WeightSfc &
            & )
        end if
        call OptPropGetAbsCoefProfForKD( &
          & iWaveNum, kmax, NMolNum, m_MolNum, z_Press, z_Temp, zm_VMR, &
          & NIntpol, &
          & z_KDTblPressIndex, z_KDTblTempIndex, &
          & z_KDTblVMRMol1Index, z_KDTblVMRMol2Index, z_KDTblVMRMol3Index, &
          & z_NIntpol1, z_NIntpol3, z_NIntpol4, z_NIntpol5, &
          & ptuvwz_Weight, &
          & zm_AbsCoefPerAtmMol &
          & )
      case default
        call OptPropGetAbsCoefProf( &
          & iWaveNum, kmax, NMolNum, m_MolNum, z_Press, z_Temp, zm_VMR, &
          & zm_AbsCoefPerAtmMol &
          & )
      end select


      if ( FlagRayleigh ) then
        call OptPropGetRayScatCoef( &
          & iWaveNum, &
          & NMolNum, m_MolNum, &
          & RayScatCoefNonRadAct, m_RayScatCoef &
          & )
      else
        RayScatCoefNonRadAct = 0.0d0
        m_RayScatCoef        = 0.0d0
      end if

      call OptPropGetPtclParam( &
        & iWaveNum, &
        & kmax, NPtcl, a_PtclName, za_PtclEffRad, &
        & za_PtclQExt, za_PtclSSA, za_PtclAF &
        & )


      ! [DEBUG] Lines for input optical depth
!!$      call OptPropLBLGetInputedOptDepProf( &
!!$        & iWaveNum, kmax, &
!!$        & z_InputedDelOptDep &
!!$        & )


      r_OptDep(kmax) = 0.0d0
      do k = kmax, 1, -1
        GasAbsCoef = 0.0d0
        do m = 1, NMolNum
          GasAbsCoef = GasAbsCoef + zm_AbsCoefPerAtmMol(k,m)
        end do
        DelGasAbsOptDep  = GasAbsCoef * z_DelAtmMols(k)
        !
        !
        ! [DEBUG] Lines for input optical depth
!!$        DelGasAbsOptDep  = z_InputedDelOptDep(k)
        !
        !
        VMRNonRadAct = 1.0d0
        do m = 1, NMolNum
          VMRNonRadAct = VMRNonRadAct - zm_VMR(k,m)
        end do
        VMRNonRadAct = max( VMRNonRadAct, 0.0d0 )
        RayScatCoef = RayScatCoefNonRadAct * VMRNonRadAct
        do m = 1, NMolNum
          RayScatCoef = RayScatCoef &
            & + m_RayScatCoef(m) * zm_VMR(k,m)
        end do
        DelRayScatOptDep = RayScatCoef * z_DelAtmMols(k)
        !
        do m = 1, NPtcl
          a_DelPtclExtOptDep(m) = &
            & za_PtclQExt(k,m) * PI * za_PtclEffRad(k,m)**2 &
            &   * za_DelPtclNum(k,m)
        end do

        DelOptDep = DelGasAbsOptDep + DelRayScatOptDep
        do m = 1, NPtcl
          DelOptDep = DelOptDep + a_DelPtclExtOptDep(m)
        end do


        r_OptDep(k-1) = r_OptDep(k) + DelOptDep
        !
        if ( DelOptDep == 0.0d0 ) then
          z_SSA(k) = 0.0d0
        else
          z_SSA(k) = 0.0d0 * DelGasAbsOptDep + 1.0d0 * DelRayScatOptDep
          do m = 1, NPtcl
            z_SSA(k) = z_SSA(k) + za_PtclSSA(k,m) * a_DelPtclExtOptDep(m)
          end do
          z_SSA(k) = z_SSA(k) / DelOptDep
        end if
        !
        if ( z_SSA(k) == 0.0d0 ) then
          z_AF(k) = 0.0d0
        else
          z_AF(k) = &
            &   0.0d0 * 0.0d0 * DelGasAbsOptDep  &
            & + 0.0d0 * 1.0d0 * DelRayScatOptDep
          do m = 1, NPtcl
            z_AF(k) = z_AF(k) &
              & + za_PtclAF(k,m) * za_PtclSSA(k,m) * a_DelPtclExtOptDep(m)
          end do
          z_AF(k) = z_AF(k) / ( z_SSA(k) * DelOptDep )
        end if

      end do

      z_SSA = min( z_SSA, 1.0d0 - 1.0d-10 )
      z_AF  = min( z_AF , 1.0d0           )


      select case ( IDMethod )
      case ( IDMethodLBL )
        call OptPropGetPFInted( &
          & iWaveNum, &
          & kmax+1, &
          & r_Temp, SurfTemp, &
          & r_PFInted, SurfPFInted, SurfDPFDTInted &
          & )
        !
        call OptPropGetStrPFInted( &
          & iWaveNum, StrFluxTOA, &
          & StrFluxTOAEach &
          & )
      case ( IDMethodLBLConstPF, IDMethodKDDevelop, IDMethodKD )

        n = OptPropGetBandNum( iWaveNum )

        if ( n >= NBand + 1 ) then
          write( 6, * ) 'In ', trim( ModuleName )
          write( 6, * ) '  Unable to find appropriate band number: ', n
          stop
        end if

        r_PFInted      = ra_BandAvePF(:,n)
        SurfPFInted    = a_BandAveSurfPF(n)
        SurfDPFDTInted = a_BandAveSurfDPFDT(n)
        !
        select case ( IDMethod )
        case ( IDMethodLBLConstPF, IDMethodKDDevelop )
          call OptPropGetPFIntedRatio( &
            & iWaveNum, kmax+1, r_Temp, SurfTemp, &
            & r_PFIntedRatio, SurfPFIntedRatio &
            & )
        case ( IDMethodKD )
!!$          call OptPropGetPFIntedRatio( &
!!$            & iWaveNum, kmax+1, NMolNum, m_MolNum, r_Press, r_Temp, SurfTemp, rm_VMR, &
!!$            & r_PFIntedRatio, SurfPFIntedRatio &
!!$            & )
          call OptPropGetPFIntedRatioForKD( &
            & iWaveNum, kmax+1, NMolNum, m_MolNum, r_Press, r_Temp, SurfTemp, rm_VMR, &
            & r_PFIntedRatio, SurfPFIntedRatio, &
            & NIntpol, &
            & r_KDTblPressIndex, r_KDTblTempIndex, &
            & r_KDTblVMRMol1Index, r_KDTblVMRMol2Index, r_KDTblVMRMol3Index, &
            & r_NIntpol1, r_NIntpol3, r_NIntpol4, r_NIntpol5, &
            & ptuvwr_Weight, &
            & r_KDTblPressIndexSfc, r_KDTblTempIndexSfc, &
            & r_KDTblVMRMol1IndexSfc, r_KDTblVMRMol2IndexSfc, r_KDTblVMRMol3IndexSfc, &
            & r_NIntpol1Sfc, r_NIntpol3Sfc, r_NIntpol4Sfc, r_NIntpol5SFc, &
            & ptuvw_WeightSfc &
            & )
        end select
        r_PFInted      = r_PFInted * r_PFIntedRatio
        SurfPFInted    = SurfPFInted * SurfPFIntedRatio
        SurfDPFDTInted = SurfDPFDTInted * SurfPFIntedRatio

        StrFluxTOAEach = a_BandAveStrFluxTOA(n)
        !
        call OptPropGetStrPFIntedRatio( &
          & iWaveNum, &
          & StrPFRatio &
          & )
        StrFluxTOAEach = StrFluxTOAEach * StrPFRatio

        r_PFInted      = r_PFInted      * DelWaveNum
        SurfPFInted    = SurfPFInted    * DelWaveNum
        SurfDPFDTInted = SurfDPFDTInted * DelWaveNum
        !
        StrFluxTOAEach = StrFluxTOAEach * DelWaveNum

      case default
        write( 6, * ) 'In ', trim( ModuleName )
        write( 6, * ) '  Unexpected case default'
        stop
      end select

      if ( iWaveNum <= iWaveNumMaxforAlbedoLW ) then
        SurfAlbedo         = AlbedoLW
      else
        SurfAlbedo         = AlbedoSW
      end if


!!$      write( 6, * ) iWaveNum, OptPropGetWN( iWaveNum ), SurfAlbedo
!!$      write( 6, * ) iWaveNum, SurfAlbedo, iWaveNumMaxforAlbedoLW


      select case ( IDMethod )
      case ( IDMethodKD )
        iBand = OptPropGetBandNum( iWaveNum )
        !
        call OptPropPFTblChkBandIntPFRatio( &
          & iBand, ThresholdBandIntPFRatio, &
          & FlagBandIntPFRatio &
          & )
        call OptPropPFTblChkBandIntStrFluxRatio( &
          & iBand, ThresholdBandIntStrFluxRatio, &
          & FlagBandIntStrFluxRatio &
          & )
      case default
        FlagBandIntPFRatio      = .true.
        FlagBandIntStrFluxRatio = .true.
      end select



      if ( FlagBandIntPFRatio ) then
        call RadRTETwoStreamApp1DLWWrapper(            &
          & kmax,                                    & ! (in)
          & z_SSA, z_AF,                             & ! (in)
          & r_OptDep,                                & ! (in)
          & SurfAlbedo,                              & ! (in)
          & r_PFInted, SurfPFInted, SurfDPFDTInted,  & ! (in)
!        & r_RadUwFluxEach, r_RadDwFluxEach,        & ! (out)
!        & ra_DelRadUwFluxEach, ra_DelRadDwFluxEach & ! (out)
          & r_RadUwFluxEachLW, r_RadDwFluxEachLW,        & ! (out)
          & ra_DelRadUwFluxEachLW, ra_DelRadDwFluxEachLW & ! (out)
          & )
      else
        r_RadUwFluxEachLW     = 0.0d0
        r_RadDwFluxEachLW     = 0.0d0
        ra_DelRadUwFluxEachLW = 0.0d0
        ra_DelRadDwFluxEachLW = 0.0d0
      end if


      if ( FlagBandIntStrFluxRatio ) then

        InAngle        = 1.0d0 / cos( SZADeg * PI / 180.0d0 )

        call RadRTETwoStreamApp1DSWWrapper(      &
          & kmax,                                    & ! (in)
          & z_SSA, z_AF,                  & ! (in)
          & r_OptDep,                     & ! (in)
          & SurfAlbedo,                   & ! (in)
          & StrFluxTOAEach, InAngle,        & ! (in)
!        & r_RadUwFluxEach, r_RadDwFluxEach      & ! (out)
          & r_RadUwFluxEachSW, r_RadDwFluxEachSW      & ! (out)
          & )
!      ra_DelRadUwFluxEach = 0.0d0
        !      ra_DelRadDwFluxEach = 0.0d0
        ra_DelRadUwFluxEachSW = 0.0d0
        ra_DelRadDwFluxEachSW = 0.0d0
      else
        r_RadUwFluxEachSW     = 0.0d0
        r_RadDwFluxEachSW     = 0.0d0
        ra_DelRadUwFluxEachSW = 0.0d0
        ra_DelRadDwFluxEachSW = 0.0d0
      end if


      r_RadUwFluxLW     = r_RadUwFluxLW     + r_RadUwFluxEachLW
      r_RadDwFluxLW     = r_RadDwFluxLW     + r_RadDwFluxEachLW
      ra_DelRadUwFluxLW = ra_DelRadUwFluxLW + ra_DelRadUwFluxEachLW
      ra_DelRadDwFluxLW = ra_DelRadDwFluxLW + ra_DelRadDwFluxEachLW
      r_RadUwFluxSW     = r_RadUwFluxSW     + r_RadUwFluxEachSW
      r_RadDwFluxSW     = r_RadDwFluxSW     + r_RadDwFluxEachSW
      ra_DelRadUwFluxSW = ra_DelRadUwFluxSW + ra_DelRadUwFluxEachSW
      ra_DelRadDwFluxSW = ra_DelRadDwFluxSW + ra_DelRadDwFluxEachSW

!!$      InAngle        = 1.0d0 / cos( SZADeg * PI / 180.0d0 )
!!$
!!$      call RadRTETwoStreamApp1DWrapper(            &
!!$        & IDScatApprox, FlagSrcFuncTech,           & ! (in)
!!$        & kmax,                                    & ! (in)
!!$        & z_SSA, z_AF,                             & ! (in)
!!$        & r_OptDep,                                & ! (in)
!!$        & SurfAlbedo,                              & ! (in)
!!$!        & r_RadUwFlux, r_RadDwFlux,                & ! (out)
!!$        & r_RadUwFluxEach, r_RadDwFluxEach,        & ! (out)
!!$!        & SolarFluxTOA, InAngle,                   & ! (in)
!!$        & StrFluxTOAEach, InAngle,                   & ! (in)
!!$        & r_PFInted, SurfPFInted, SurfDPFDTInted,  & ! (in)
!!$!        & ra_DelRadUwFlux, ra_DelRadDwFlux         & ! (out)
!!$        & ra_DelRadUwFluxEach, ra_DelRadDwFluxEach & ! (out)
!!$        & )


!!$      if ( iWaveNumEnd < 10000 ) then
!!$      select case ( IDMethod )
!!$      case ( IDMethodLBLConstPF, IDMethodKDDevelop )
!!$        write( 60, * ) iWaveNum, r_OptDep(0), StrFluxTOAEach, r_RadDwFluxEach(0), pm_AbsCoefPerAtmMol(0,1), StrFluxTOAEach, InAngle
!!$        write( 60, * ) iWaveNum, r_PFInted(0), SurfPFInted, SurfDPFDTInted, &
!!$          & IDScatApprox, FlagSrcFuncTech
!!$      end select
!!$      end if

      if ( FlagOutputInited ) then
        call ni3_put_varss( NcIDFlux, "RadUwFluxLW", iWaveNum, r_RadUwFluxEachLW/DelWaveNum )
        call ni3_put_varss( NcIDFlux, "RadDwFluxLW", iWaveNum, r_RadDwFluxEachLW/DelWaveNum )
        call ni3_put_varss( NcIDFlux, "RadUwFluxSW", iWaveNum, r_RadUwFluxEachSW/DelWaveNum )
        call ni3_put_varss( NcIDFlux, "RadDwFluxSW", iWaveNum, r_RadDwFluxEachSW/DelWaveNum )
        r_RadUwFluxEach = r_RadUwFluxEachLW + r_RadUwFluxEachSW
        r_RadDwFluxEach = r_RadDwFluxEachLW + r_RadDwFluxEachSW
        call ni3_put_varss( NcIDFlux, "RadUwFlux", iWaveNum, r_RadUwFluxEach/DelWaveNum )
        call ni3_put_varss( NcIDFlux, "RadDwFlux", iWaveNum, r_RadDwFluxEach/DelWaveNum )
        call ni3_put_varss( NcIDFlux, "OptDep"   , iWaveNum, r_OptDep )
        do k = 1, kmax
          z_DelOptDep(k) = r_OptDep(k-1) - r_OptDep(k)
        end do
        call ni3_put_varss( NcIDFlux, "DelOptDep", iWaveNum, z_DelOptDep )
        do k = 1, kmax
          z_DRadFluxDPressEach(k) = &
            &   (   ( r_RadUwFluxEach(k  ) - r_RadDwFluxEach(k  ) ) &
            &     - ( r_RadUwFluxEach(k-1) - r_RadDwFluxEach(k-1) ) ) &
            & / ( r_Press(k) - r_Press(k-1) )
        end do
        z_DRadFluxDPressEach = z_DRadFluxDPressEach / DelWaveNum
        call ni3_put_varss( NcIDFlux, "DRadFluxDPress", iWaveNum, z_DRadFluxDPressEach )
      end if

    end do

    r_RadUwFlux     = r_RadUwFluxLW     + r_RadUwFluxSW
    r_RadDwFlux     = r_RadDwFluxLW     + r_RadDwFluxSW
    ra_DelRadUwFlux = ra_DelRadUwFluxLW + ra_DelRadUwFluxSW
    ra_DelRadDwFlux = ra_DelRadDwFluxLW + ra_DelRadDwFluxSW

!!$    deallocate( zm_AbsCoefPerAtmMol )


    select case ( IDMethod )
    case ( IDMethodLBLConstPF, IDMethodKDDevelop, IDMethodKD )
      deallocate( a_BandAveSurfPF     )
      deallocate( ra_BandAvePF        )
      deallocate( a_BandAveSurfDPFDT  )
      deallocate( a_BandAveStrFluxTOA )
    end select


  end subroutine RadFluxCalcFlux

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

  subroutine RadFluxInit( &
    & ArgGrav, &
    & WaveNumAlbedoSwitch, &
    & OptPropNcFN, &
    & ArgRayScatCoefNcFN, ArgStrSpeNcFN, &
    & Method, &
    & a_ArgBandWNBnds &
    & )

    use constants0, only : StB

    use opt_prop, only : &
      & OptPropInit, &
      & OptPropGetNBand, &
      & OptPropGetBandWNBnds, &
      & OptPropGetBandBinIndexBnds
!!$      & OptPropGetTotalStrFlux

    use rad_rte_two_stream_app_1d_wrapper, only : RadRTETwoStreamApp1DWrapperInit

    real(8)     , intent(in )           :: ArgGrav
    real(DP)    , intent(in )           :: WaveNumAlbedoSwitch
    character(*), intent(in )           :: OptPropNcFN
    character(*), intent(in ), optional :: ArgRayScatCoefNcFN
    character(*), intent(in ), optional :: ArgStrSpeNcFN
    character(*), intent(in ), optional :: Method
    real(DP)    , intent(in ), optional :: a_ArgBandWNBnds(:)


    ! local variables
    !
    character(extstr) :: MethodLocal
    character(extstr) :: RayScatCoefNcFN
    character(extstr) :: StrSpeNcFN

    integer  :: n


    if ( FlagInited ) return


    if ( present( Method ) ) then
      select case ( Method )
      case ( 'KD' )
        IDMethod = IDMethodKD
      case ( 'KDDevelop' )
        IDMethod = IDMethodKDDevelop
      case ( 'LBL' )
        IDMethod = IDMethodLBL
      case ( 'LBLConstPF' )
        IDMethod = IDMethodLBLConstPF
      case default
        write( 6, * ) "In module, ", trim( ModuleName )
        write( 6, * ) "  Unexpected method: ", trim( Method )
        stop
      end select
      MethodLocal = Method
    else
      MethodLocal = 'KD'
      IDMethod    = IDMethodKD
    end if


    ! Check arguments
    !
    select case ( IDMethod )
    case ( IDMethodLBLConstPF )
      if ( present( a_ArgBandWNBnds ) ) then
      else
        write( 6, * ) "In module, ", trim( ModuleName )
        write( 6, * ) 'a_BandWNBnds has to be present when IDMethod is IDMethodLBLConstPF.'
        stop
      end if
    end select


    Grav      = ArgGrav

    if ( present( ArgRayScatCoefNcFN ) ) then
      RayScatCoefNcFN = ArgRayScatCoefNcFN
    else
      RayScatCoefNcFN = ""
    end if
    if ( present( ArgStrSpeNcFN ) ) then
      StrSpeNcFN = ArgStrSpeNcFN
    else
      StrSpeNcFN = ""
    end if



    call OptPropInit( &
      & MethodLocal, &
      & OptPropNcFN, RayScatCoefNcFN, StrSpeNcFN &
      & )

    call RadRTETwoStreamApp1DWrapperInit


    select case ( IDMethod )
    case ( IDMethodKDDevelop, IDMethodKD )

      NBand = OptPropGetNBand()
      allocate( aa_BandWNBnds      (2,NBand) )
      allocate( aa_BandBinIndexBnds(2,NBand) )
      call OptPropGetBandWNBnds( &
        & NBand, &
        & aa_BandWNBnds &
        & )
      call OptPropGetBandBinIndexBnds( &
        & NBand, &
        & aa_BandBinIndexBnds &
        & )

    case ( IDMethodLBLConstPF )

      NBand = size( a_ArgBandWNBnds ) - 1
      allocate( aa_BandWNBnds(2,NBand) )

      do n = 1, NBand
        aa_BandWNBnds(1,n) = a_ArgBandWNBnds(n  )
        aa_BandWNBnds(2,n) = a_ArgBandWNBnds(n+1)
      end do

    end select


    ! Check index for albedo switch from LW to SW
    iWaveNumMaxforAlbedoLW = SearchBoundaryIndex( WaveNumAlbedoSwitch )


!!$    call OptPropGetTotalStrFlux( &
!!$      & TotalStrFlux &
!!$      & )
!!$    FluxRatTOAtoStelSfc = StrFluxTOA / TotalStrFlux


    FlagInited = .true.


  end subroutine RadFluxInit

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

  subroutine RadFluxFinalize

    use opt_prop, only : OptPropFinalize


    ! local variables
    !


    call OptPropFinalize

    if ( FlagOutputInited ) then
      call RadFluxOutputFinalize
    end if


    select case ( IDMethod )
    case ( IDMethodKDDevelop )
      deallocate( aa_BandWNBnds       )
      deallocate( aa_BandBinIndexBnds )
    case ( IDMethodLBLConstPF )
      deallocate( aa_BandWNBnds )
    end select


    FlagInited = .false.

  end subroutine RadFluxFinalize

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

  function SearchBoundaryIndex( WaveNumSwitch ) result( iWaveNum )

    use opt_prop, only : &
      & OptPropGetNWN, &
      & OptPropGetWN

    real(DP), intent(in ) :: WaveNumSwitch

    integer               :: iWaveNum


    ! Local variables
    !
    integer :: NWaveNum
    integer :: iBand


    NWaveNum   = OptPropGetNWN()

    ! Check index for switching from LW to SW
    select case ( IDMethod )
    case ( IDMethodKDDevelop, IDMethodKD )
      if ( WaveNumSwitch < aa_BandWNBnds(1,1) ) then
        ! ... for short wave is used in all bands.
        iWaveNum = 0
      else if ( aa_BandWNBnds(2,NBand) < WaveNumSwitch ) then
        ! ... for long wave is used in all bands.
        iWaveNum = NWaveNum
      else
        do iBand = 1, NBand
          if ( ( aa_BandWNBnds(1,iBand) <= WaveNumSwitch          ) .and. &
            &  ( WaveNumSwitch          <= aa_BandWNBnds(2,iBand) ) ) then
            exit
          end if
        end do
        if (   ( WaveNumSwitch          - aa_BandWNBnds(1,iBand) ) &
          &  < ( aa_BandWNBnds(2,iBand) - WaveNumSwitch           ) ) then
          iBand = iBand - 1
        end if
        ! obtain max bin number for ... for long wave
        if ( iBand >= 1 ) then
          iWaveNum = aa_BandBinIndexBnds(2,iBand)
        else
          iWaveNum = 0
        end if
      end if
    case ( IDMethodLBL, IDMethodLBLConstPF )
      if ( WaveNumSwitch < OptPropGetWN( 1 ) ) then
        ! ... for short wave is used in all bands.
        iWaveNum = 0
      else if ( OptPropGetWN( NWaveNum ) < WaveNumSwitch ) then
        ! ... for long wave is used in all bands.
        iWaveNum = NWaveNum
      else
        do iWaveNum = 1, NWaveNum
          if ( WaveNumSwitch <= OptPropGetWN( iWaveNum ) ) exit
        end do
        iWaveNum = iWaveNum - 1
      end if
    case default
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Unexpected case default.'
      stop
    end select

  end function SearchBoundaryIndex

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

  subroutine RadFluxOutputInit( &
    & OptPropNcFN, &
    & FluxOutNcFn, &
    & kmax, r_Press &
    & )

    use ni3_module
    use netcdf

    character(*), intent(in) :: OptPropNcFN
    character(*), intent(in) :: FluxOutNcFn
    integer     , intent(in) :: kmax
    real(8)     , intent(in) :: r_Press    (0:kmax)


    ! local variables
    !
    real(8)           :: z_Press    (1:kmax)
    integer           :: NcID
    character(stdstr) :: Mode
    character(stdstr) :: Name
    character(stdstr) :: StdName
    character(stdstr) :: LongName
    character(stdstr) :: Units
    character(stdstr) :: aa_DimNames(2)

    integer :: LBLNWN
    real(DP), allocatable :: w_WaveNum(:)

    integer :: k


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    do k = 1, kmax
      z_Press(k) = ( r_Press(k-1) + r_Press(k) ) / 2.0d0
    end do


    Mode = "read"
    call ni3_open( OptPropNcFN, Mode, NcID )
    call ni3_inq_dimlen( NcID, 'WaveNum', LBLNWN )
    allocate( w_WaveNum( 1:LBLNWN ) )
    call ni3_get_var( NcID, 'WaveNum', w_WaveNum )
    call ni3_close( NcID )


    ! Output

    Mode = "new"
    call ni3_open( FluxOutNcFN, Mode, NcIDFlux )

    Name    = "Press"
    StdName = "plev"
    Units   = "Pa"
    call ni3_set_dim( NcIDFlux, Name, NF90_DOUBLE, r_Press  , &
      & stdname = stdname, units  = units )
    Name    = "WaveNum"
    StdName = "wavenumber"
    Units   = "m-1"
    call ni3_set_dim( NcIDFlux, Name, NF90_DOUBLE, w_WaveNum , &
      & stdname = stdname, units  = units )
    Name    = "PressLC"
    StdName = "plev"
    Units   = "Pa"
    call ni3_set_dim( NcIDFlux, Name, NF90_DOUBLE, z_Press  , &
      & stdname = stdname, units  = units )


    aa_DimNames(1) = "Press"
    aa_DimNames(2) = "WaveNum"

    Name     = "RadUwFlux"
    StdName  = Name
    LongName = "radiative upward flux"
    Units   = "W m-2 (m-1)-1"
    call ni3_def_var( NcIDFlux, Name, NF90_DOUBLE, 2, aa_DimNames, &
      & stdname = StdName, longname = LongName, Units = Units )

    Name     = "RadDwFlux"
    StdName  = Name
    LongName = "radiative downward flux"
    Units   = "W m-2 (m-1)-1"
    call ni3_def_var( NcIDFlux, Name, NF90_DOUBLE, 2, aa_DimNames, &
      & stdname = StdName, longname = LongName, Units = Units )

    Name     = "RadUwFluxLW"
    StdName  = Name
    LongName = "long wave radiative upward flux"
    Units   = "W m-2 (m-1)-1"
    call ni3_def_var( NcIDFlux, Name, NF90_DOUBLE, 2, aa_DimNames, &
      & stdname = StdName, longname = LongName, Units = Units )

    Name     = "RadDwFluxLW"
    StdName  = Name
    LongName = "long wave radiative downward flux"
    Units   = "W m-2 (m-1)-1"
    call ni3_def_var( NcIDFlux, Name, NF90_DOUBLE, 2, aa_DimNames, &
      & stdname = StdName, longname = LongName, Units = Units )

    Name     = "RadUwFluxSW"
    StdName  = Name
    LongName = "short wave radiative upward flux"
    Units   = "W m-2 (m-1)-1"
    call ni3_def_var( NcIDFlux, Name, NF90_DOUBLE, 2, aa_DimNames, &
      & stdname = StdName, longname = LongName, Units = Units )

    Name     = "RadDwFluxSW"
    StdName  = Name
    LongName = "short wave radiative downward flux"
    Units   = "W m-2 (m-1)-1"
    call ni3_def_var( NcIDFlux, Name, NF90_DOUBLE, 2, aa_DimNames, &
      & stdname = StdName, longname = LongName, Units = Units )

    Name     = "OptDep"
    StdName  = Name
    LongName = "optical depth"
    Units   = "1"
    call ni3_def_var( NcIDFlux, Name, NF90_DOUBLE, 2, aa_DimNames, &
      & stdname = StdName, longname = LongName, Units = Units )


    aa_DimNames(1) = "PressLC"
    aa_DimNames(2) = "WaveNum"

    Name     = "DRadFluxDPress"
    StdName  = Name
    LongName = "flux divergence"
    Units    = "W m-2 Pa-1 (m-1)-1"
    call ni3_def_var( NcIDFlux, Name, NF90_DOUBLE, 2, aa_DimNames, &
      & stdname = StdName, longname = LongName, Units = Units )

    Name     = "DelOptDep"
    StdName  = Name
    LongName = "optical depth increment"
    Units   = "1"
    call ni3_def_var( NcIDFlux, Name, NF90_DOUBLE, 2, aa_DimNames, &
      & stdname = StdName, longname = LongName, Units = Units )

    deallocate( w_WaveNum )


    FlagOutputInited = .true.

  end subroutine RadFluxOutputInit

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

  subroutine RadFluxOutputFinalize

    use ni3_module

    ! local variables
    !


    if ( .not. FlagInited ) then
      write( 6, * ) trim( ModuleName ), " is not initialized."
      stop
    end if


    if ( FlagOutputInited ) then
      call ni3_close( NcIDFlux )
      FlagOutputInited = .false.
    end if


  end subroutine RadFluxOutputFinalize

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

end module rad_flux
