module band_module

  use dc_types

  implicit none

  private

  public :: BandInit
  public :: BandFinish
  public :: GetNBand
  public :: GetBandNum
  public :: GetSurfPFInted
  public :: GetPFInted
  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_SolarFluxTOA(:)


contains

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

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


    real(8), intent(in) :: PI
    integer, intent(in) :: Arg_kmin
    integer, intent(in) :: Arg_kmax
    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
    !
    integer :: k
    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_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
       a_SurfPFInted(n) = PI * integral_planck(a_WNBnds(n-1), a_WNBnds(n), r_Temp(0))
       do k = kmin, kmax
          ra_PFInted(k,n) = PI * integral_planck(a_WNBnds(n-1), a_WNBnds(n), r_Temp(k))
       end do
       a_SolarFluxTOA(n) = PI * RatioSolarConst * integral_planck(a_WNBnds(n-1), a_WNBnds(n), SolTemp)
    end do

    do n = 1, NBand
!!$      a_SurfPFInted(n)  = a_SurfPFInted(n)  * DelWN / ( a_WNBnds(n) - a_WNBnds(n-1) )
!!$      ra_PFInted(:,n)   = ra_PFInted(:,n)   * DelWN / ( a_WNBnds(n) - a_WNBnds(n-1) )
!!$      a_SolarFluxTOA(n) = a_SolarFluxTOA(n) * DelWN / ( a_WNBnds(n) - a_WNBnds(n-1) )
      a_SurfPFInted(n)  = a_SurfPFInted(n)  / ( a_WNBnds(n) - a_WNBnds(n-1) )
      ra_PFInted(:,n)   = ra_PFInted(:,n)   / ( a_WNBnds(n) - a_WNBnds(n-1) )
      a_SolarFluxTOA(n) = a_SolarFluxTOA(n) / ( a_WNBnds(n) - a_WNBnds(n-1) )
    end do


  end subroutine BandInit

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

  subroutine BandFinish

    ! Local variables
    !


    deallocate( a_WNBnds       )
    deallocate( a_SurfPFInted  )
    deallocate( ra_PFInted     )
    deallocate( a_SolarFluxTOA )


  end subroutine BandFinish

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

  function GetNBand() result( NumBand )

    integer :: NumBand

    NumBand = NBand

  end function GetNBand

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

  function GetBandNum( WaveNum ) result( BandNum )

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

    integer :: BandNum

    integer :: n

    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


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

  end function GetSurfPFInted

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

  function GetPFInted( k, BandNum ) result( PFInted )

    integer , intent(in) :: k
    integer , intent(in) :: BandNum

    real(DP) :: PFInted


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

  end function GetPFInted

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

  function GetSolarPFInted( BandNum ) result( SolarPFInted )

    integer , intent(in) :: BandNum

    real(DP) :: SolarPFInted


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

  end function GetSolarPFInted

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

end module band_module
