module band

  use dc_types

  implicit none

  private

  public :: BandInit
  public :: BandFinalize
  public :: BandGetPFInted
  public :: BandGetNBand
  public :: GetBandNum
  public :: GetSurfPFInted
  public :: r_GetPFInted
  public :: GetSurfDPFDTInted
!!$  public :: GetSolarPFInted


  integer, save              :: NBand
  integer, save              :: kmin
  integer, save              :: kmax
  real(8), allocatable, save :: a_WNBnds        (:)
  real(8), allocatable, save :: a_SurfPFInted   (:)
  real(8), allocatable, save :: ra_PFInted      (:,:)
  real(8), allocatable, save :: a_SurfDPFDTInted(:)
!!$  real(8), allocatable, save :: a_SolarFluxTOA  (:)

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

  logical, save :: FlagInited = .false.


contains

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

  subroutine BandInit( &
    & PI, Arg_kmin, Arg_kmax, SurfTemp, r_Temp, &
!!$    & SolTemp, RatioSolarConst, &
    & ArgNBand, a_ArgWNBnds &
    & )

    use planck_func_wrapper


    real(8), intent(in) :: PI
    integer, intent(in) :: Arg_kmin
    integer, intent(in) :: Arg_kmax
    real(8), intent(in) :: SurfTemp
    real(8), intent(in) :: r_Temp(Arg_kmin:Arg_kmax)
!!$    real(8), intent(in) :: SolTemp
!!$    real(8), intent(in) :: RatioSolarConst
    integer, intent(in) :: ArgNBand
    real(8), intent(in) :: a_ArgWNBnds(0:ArgNBand)

    ! Local variables
    !
    real(8) :: WNs
    real(8) :: WNe
    integer :: GQNum
    integer :: n


    kmin = Arg_kmin
    kmax = Arg_kmax

    NBand = ArgNBand
!!$    NBand = 9

    allocate( a_WNBnds      (0:NBand) )
    allocate( a_SurfPFInted (NBand) )
    allocate( ra_PFInted    (kmin:kmax,NBand) )
    allocate( a_SurfDPFDTInted(NBand) )
!!$    allocate( a_SolarFluxTOA(NBand) )

    a_WNBnds = a_ArgWNBnds

!!$    a_WNBnds( 0) =    0.0d2
!!$    a_WNBnds( 1) =  340.0d2
!!$    a_WNBnds( 2) =  540.0d2
!!$    a_WNBnds( 3) =  800.0d2
!!$    a_WNBnds( 4) =  980.0d2
!!$    a_WNBnds( 5) = 1100.0d2
!!$    a_WNBnds( 6) = 1215.0d2
!!$    a_WNBnds( 7) = 1380.0d2
!!$    a_WNBnds( 8) = 1900.0d2
!!$    a_WNBnds( 9) = 3000.0d2


    do n = 1, NBand
      WNs = a_WNBnds(n-1)
      WNe = a_WNBnds(n  )
      GQNum = int( ( WNe - WNs ) / 100.0d0 ) + 1
      call Integ_PF_GQ_Array0D( &
        & WNs, WNe, GQNum, &
        & SurfTemp, &
        & a_SurfPFInted(n) &
        & )
      a_SurfPFInted(n) = PI * a_SurfPFInted(n)
      call Integ_PF_GQ_Array1D( &
        & WNs, WNe, GQNum, &
        & 0, kmax, &
        & r_Temp, &
        & ra_PFInted(:,n) &
        & )
      ra_PFInted(:,n) = PI * ra_PFInted(:,n)
      call Integ_DPFDT_GQ_Array0D( &
        & WNs, WNe, GQNum,                 & ! (in )
        & SurfTemp,                & ! (in )
        & a_SurfDPFDTInted(n) &
        & )
      a_SurfDPFDTInted(n) = PI * a_SurfDPFDTInted(n)
    end do


    FlagInited = .true.


  end subroutine BandInit

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

  subroutine BandFinalize

    ! Local variables
    !


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


    deallocate( a_WNBnds         )
    deallocate( a_SurfPFInted    )
    deallocate( ra_PFInted       )
    deallocate( a_SurfDPFDTInted )
!!$    deallocate( a_SolarFluxTOA   )


    FlagInited = .false.

  end subroutine BandFinalize

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

  subroutine BandGetPFInted( &
    & WaveNum, &
    & kmax, &
    & r_PFInted, SurfPFInted, SurfDPFDTInted &
    & )

    real(DP), intent(in ) :: WaveNum
    integer , intent(in ) :: kmax
    real(DP), intent(out) :: r_PFInted  (0:kmax)
    real(DP), intent(out) :: SurfPFInted
    real(DP), intent(out) :: SurfDPFDTInted


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


    if ( kmax /= size(ra_PFInted(:,1)) ) then
      write( 6, * ) trim( ModuleName )
      write( 6, * ) 'Unexpected number of layers'
      stop
    end if


    r_PFInted      = r_GetPFInted( GetBandNum( WaveNum ) )
    SurfPFInted    = GetSurfPFInted( GetBandNum( WaveNum ) )
    SurfDPFDTInted = GetSurfDPFDTInted( GetBandNum( WaveNum ) )


  end subroutine BandGetPFInted

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

  function BandGetNBand() result( NumBand )

    integer :: NumBand


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

    NumBand = NBand

  end function BandGetNBand

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

  function GetBandNum( WaveNum ) result( BandNum )

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

    integer :: BandNum

    integer :: n


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

    search_Band : do n = 1, NBand
       if ( WaveNum <= a_WNBnds(n) ) exit search_Band
    end do search_Band

    if ( n >= NBand + 1 ) then
!!$      write( 6, * ) 'Given wavenumber is out of given bands: ', &
!!$        & WaveNum, a_WNBnds(NBand)
!!$      stop
      n = -999
    end if

    BandNum = n

  end function GetBandNum

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

  function GetSurfPFInted( BandNum ) result( SurfPFInted )

    integer , intent(in) :: BandNum

    real(DP) :: SurfPFInted


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

!!$    BandNum = GetBandNum(WaveNum)
    if ( BandNum == -999 ) then
      SurfPFInted = 0.0d0
    else
      SurfPFInted = a_SurfPFInted(BandNum)
    end if

  end function GetSurfPFInted

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

  function r_GetPFInted( BandNum ) result( r_PFInted )

    integer , intent(in) :: BandNum

    real(DP) :: r_PFInted( size(ra_PFInted(:,1)) )


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

!!$    BandNum = GetBandNum(WaveNum)
    if ( BandNum == -999 ) then
      r_PFInted = 0.0d0
    else
      r_PFInted = ra_PFInted(:, BandNum)
    end if

  end function r_GetPFInted

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

  function GetSurfDPFDTInted( BandNum ) result( SurfDPFDTInted )

    integer , intent(in) :: BandNum

    real(DP) :: SurfDPFDTInted


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

!!$    BandNum = GetBandNum(WaveNum)
    if ( BandNum == -999 ) then
      SurfDPFDTInted = 0.0d0
    else
      SurfDPFDTInted = a_SurfDPFDTInted(BandNum)
    end if

  end function GetSurfDPFDTInted

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

!!$  function GetSolarPFInted( BandNum ) result( SolarPFInted )
!!$
!!$    integer , intent(in) :: BandNum
!!$
!!$    real(DP) :: SolarPFInted
!!$
!!$
!!$    if ( .not. FlagInited ) then
!!$      write( 6, * ) trim( ModuleName ), " is not initialized."
!!$      stop
!!$    end if
!!$
!!$!    BandNum = GetBandNum(WaveNum)
!!$    if ( BandNum == -999 ) then
!!$      SolarPFInted = 0.0d0
!!$    else
!!$      SolarPFInted = a_SolarFluxTOA(BandNum)
!!$    end if
!!$
!!$  end function GetSolarPFInted

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

end module band
