module opt_prop_lblconstpf

  use vtype_module

  implicit none

  private

  public :: OptPropLBLConstPFInit
  public :: OptPropLBLConstPFFinalize
  public :: OptPropLBLConstPFGetNWN
  public :: OptPropLBLConstPFGetNPress
  public :: OptPropLBLConstPFGetNMol
  public :: OptPropLBLConstPFGetNBand
  public :: OptPropLBLConstPFGetDelWN
  public :: OptPropLBLConstPFGetWN
  public :: OptPropLBLConstPFGetBandNum
  public :: OptPropLBLConstPFGetBandWNBnds
!!$  public :: OptPropLBLConstPFGetTotalStrFlux
  public :: OptPropLBLConstPFGetAbsCoefProf
!!$  public :: OptPropLBLConstPFGetPFInted
!!$  public :: OptPropLBLConstPFGetStrPFInted
  public :: OptPropLBLConstPFGetStrBandAvePF
  public :: OptPropLBLConstPFGetPFIntedRatio
  public :: OptPropLBLConstPFGetStrPFIntedRatio


  integer, save               :: IDMethod
  integer, parameter          :: IDMethodLBL        =  1
  integer, parameter          :: IDMethodLBLConstPF =  2
  integer, parameter          :: IDMethodKDDevelop  = 10
  integer, parameter          :: IDMethodKD         = 11


  ! Common
  character(128)       , save :: OptPropNcFNSave
  character(128)       , save :: StrSpeNcFNSave
  !
  character(128)       , save :: PressAxisName
  character(128)       , save :: MolAxisName
  character(128)       , save :: PressForPFAxisName
  !
  integer              , save :: NPressSave
  integer              , save :: NMolSave
  integer              , save :: NPressForPFSave
  !
  real(DP), allocatable, save :: p_PressSave(:)
  integer , allocatable, save :: m_MolNumSave(:)

  ! LBL
  integer              , save :: NWaveNumSave
  real(DP), allocatable, save :: w_WaveNumSave(:)

!!$  real(DP), allocatable, save :: b_BandAveStrPFSave(:)
  real(DP), allocatable, save :: b_RatioBandIntFluxtoTotFluxSave(:)

  real(DP), allocatable, save :: aa_BandWaveNumBndsSave     (:,:)

  integer , parameter         :: NBndSave = 2
  integer              , save :: NBandSave
  !
  real(DP), allocatable, save :: ra_AtmPFRatio      (:,:)
  real(DP), allocatable, save :: zma_AbsCoef     (:,:,:)
  integer , allocatable, save :: ma_IDAbsCoefType(:,:)

  real(DP), allocatable, save :: a_StrPFRatio(:)

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

  logical, save :: FlagInited = .false.

contains

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

  subroutine OptPropLBLConstPFInit( &
    & OptPropNcFN, StrSpeNcFN &
    & )

    use ni3_module


    character(*), intent(in) :: OptPropNcFN
    character(*), intent(in) :: StrSpeNcFN


    ! local variables
    !
    character(stdstr) :: Mode
    integer           :: NcID

    character(128) :: WaveNumAxisName
    character(128) :: AbsCoefVarName
    character(128) :: AtmPFRatioVarName
    character(128) :: StrPFRatioVarName
    character(128) :: StrSpeVarName

    real(DP), allocatable :: ra_StrPFRatio(:,:)
!!$    real(DP), allocatable :: pb_BandAveStrPF(:,:)
    real(DP), allocatable :: rb_RatioBandIntFluxtoTotFlux(:,:)

    integer :: iWaveNum
    integer :: iBand

    integer :: TmpNum


    NPressSave      = -1
    NMolSave        = -1
    NWaveNumSave    = -1
    NBandSave       = -1
    NPressForPFSave = -1


    OptPropNcFNSave = OptPropNcFN
    StrSpeNcFNSave  = StrSpeNcFN

    PressAxisName      = "Press"
    MolAxisName        = "MolNum"
    PressForPFAxisName = "PressForPF"

    Mode = "read"
    call ni3_open( OptPropNcFNSave, Mode, NcID )
    call ni3_inq_dimlen( NcID, PressAxisName     , NPressSave      )
    call ni3_inq_dimlen( NcID, MolAxisName       , NMolSave        )
    call ni3_inq_dimlen( NcID, PressForPFAxisName, NPressForPFSave )
    allocate( p_PressSave (1:NPressSave) )
    allocate( m_MolNumSave(1:NMolSave  ) )
    call ni3_get_var( NcID, PressAxisName, p_PressSave  )
    call ni3_get_var( NcID, MolAxisName  , m_MolNumSave )


    WaveNumAxisName = 'WaveNum'
    AbsCoefVarName  = 'AbsCoef'
    AtmPFRatioVarName = 'AtmPFRatio'
    StrPFRatioVarName = 'StrPFRatio'
    !
    call ni3_inq_dimlen( NcID, WaveNumAxisName, NWaveNumSave )
    allocate( w_WaveNumSave( 1:NWaveNumSave ) )
    call ni3_get_var( NcID, WaveNumAxisName, w_WaveNumSave )

    allocate( zma_AbsCoef (NPressSave,NMolSave,NWaveNumSave) )
    call ni3_get_var( NcID, AbsCoefVarName, zma_AbsCoef )

    allocate( ma_IDAbsCoefType(NMolSave,NWaveNumSave) )
    call ni3_get_var( NcID, "IDAbsCoefType", ma_IDAbsCoefType )

    allocate( ra_AtmPFRatio(NPressForPFSave,NWaveNumSave) )
    call ni3_get_var( NcID, AtmPFRatioVarName, ra_AtmPFRatio )
    allocate( ra_StrPFRatio(NPressForPFSave,NWaveNumSave) )
    call ni3_get_var( NcID, StrPFRatioVarName, ra_StrPFRatio )
    allocate( a_StrPFRatio(NWaveNumSave) )
    do iWaveNum = 1, NWaveNumSave
      a_StrPFRatio(iWaveNum) = ra_StrPFRatio(1,iWaveNum)
    end do
    deallocate( ra_StrPFRatio )

    call ni3_inq_dimlen( NcID, "BandNum", NBandSave )
    allocate( aa_BandWaveNumBndsSave     (NBndSave,NBandSave) )
    call ni3_get_var( NcID, 'BandWaveNumBnds' , aa_BandWaveNumBndsSave  )

    call ni3_close( NcID )


!!$    allocate( b_BandAveStrPFSave( NBandSave ) )
!!$    allocate( pb_BandAveStrPF   ( NPressSave, NBandSave ) )
    allocate( b_RatioBandIntFluxtoTotFluxSave( NBandSave ) )
    allocate( rb_RatioBandIntFluxtoTotFlux( NPressForPFSave, NBandSave ) )

!!$    StrSpeVarName   = 'BandAveStrPF'
    StrSpeVarName   = 'RatioBandIntFluxtoTotFlux'
    !
    Mode = "read"
    call ni3_open( StrSpeNcFNSave, Mode, NcID )
    call ni3_inq_dimlen( NcID, "BandNum", TmpNum )
    if ( TmpNum /= NBandSave ) then
      write( 6, * ) 'Array size for stellar spectrum is inappropriate: '
      write( 6, * ) '  ', trim( StrSpeNcFN  ), ': ', TmpNum
      write( 6, * ) '  ', trim( OptPropNcFN ), ': ', NWaveNumSave
      stop
    end if
!!$    call ni3_get_var( NcID, "BandAveStrPF", pb_BandAveStrPF )
    call ni3_get_var( NcID, "RatioBandIntFluxtoTotFlux", rb_RatioBandIntFluxtoTotFlux )
    call ni3_close( NcID )
    !
    do iBand = 1, NBandSave
!!$      b_BandAveStrPFSave(iBand) = pb_BandAveStrPF(1,iBand)
      b_RatioBandIntFluxtoTotFluxSave(iBand) = rb_RatioBandIntFluxtoTotFlux(1,iBand)
    enddo
    !
!!$    deallocate( pb_BandAveStrPF )
    deallocate( rb_RatioBandIntFluxtoTotFlux )







    FlagInited = .true.


  end subroutine OptPropLBLConstPFInit

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

  subroutine OptPropLBLConstPFFinalize


    ! local variables
    !


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


    NPressSave      = -1
    NMolSave        = -1
    NWaveNumSave    = -1
    NBandSave       = -1
    NPressForPFSave = -1


    deallocate( p_PressSave  )
    deallocate( m_MolNumSave )

    NWaveNumSave    = 0
    deallocate( w_WaveNumSave )
    select case ( IDMethod )
    case ( IDMethodLBLConstPF )
      deallocate( aa_BandWaveNumBndsSave )
    end select

    deallocate( zma_AbsCoef      )

    deallocate( ma_IDAbsCoefType )

    deallocate( aa_BandWaveNumBndsSave )

    deallocate( ra_AtmPFRatio )
    deallocate( a_StrPFRatio )

!!$    deallocate( b_BandAveStrPFSave )
    deallocate( b_RatioBandIntFluxtoTotFluxSave )


  end subroutine OptPropLBLConstPFFinalize

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

  function OptPropLBLConstPFGetNWN() result( NWN )


    integer :: NWN

    ! local variables
    !


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


    NWN = NWaveNumSave


  end function OptPropLBLConstPFGetNWN

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

  function OptPropLBLConstPFGetNPress() result( NPress )


    integer :: NPress

    ! local variables
    !


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


    NPress = NPressSave


  end function OptPropLBLConstPFGetNPress

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

  function OptPropLBLConstPFGetNMol() result( NMol )


    integer :: NMol

    ! local variables
    !


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


    NMol = NMolSave


  end function OptPropLBLConstPFGetNMol

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

  function OptPropLBLConstPFGetNBand() result( NBand )


    integer :: NBand

    ! local variables
    !


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


    NBand = NBandSave


  end function OptPropLBLConstPFGetNBand

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

  function OptPropLBLConstPFGetDelWN() result( DelWaveNum )


    real(DP) :: DelWaveNum

    ! local variables
    !


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


    DelWaveNum = w_WaveNumSave( 2 ) - w_WaveNumSave( 1 )


  end function OptPropLBLConstPFGetDelWN

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

  function OptPropLBLConstPFGetWN( iWaveNum ) result( WaveNum )


    integer, intent(in) :: iWaveNum

    real(DP) :: WaveNum

    ! local variables
    !


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


    WaveNum = w_WaveNumSave( iWaveNum )


  end function OptPropLBLConstPFGetWN

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

  function OptPropLBLConstPFGetBandNum( iWaveNum ) result( BandNum )


    integer, intent(in) :: iWaveNum

    integer :: BandNum


    ! local variables
    !
    real(Dp) :: WaveNum
    integer  :: n


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


    WaveNum = w_WaveNumSave( iWaveNum )
    do n = 1, NBandSave
      if ( WaveNum <= aa_BandWaveNumBndsSave(2,n) ) exit
    end do

    BandNum = n


  end function OptPropLBLConstPFGetBandNum

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

  subroutine OptPropLBLConstPFGetBandWNBnds( &
    & NBand, &
    & aa_BandWNBnds &
    & )


    integer , intent(in ) :: NBand
    real(DP), intent(out) :: aa_BandWNBnds(2,NBand)

    ! local variables
    !


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


    ! check
    if ( NBand /= NBandSave ) then
      write( 6, * ) 'Number of bands is inappropriate: ', NBand, NBandSave
      stop
    end if
    aa_BandWNBnds = aa_BandWaveNumBndsSave


  end subroutine OptPropLBLConstPFGetBandWNBnds

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

!!$  subroutine OptPropLBLConstPFGetTotalStrFlux( &
!!$    & TotalStrFlux &
!!$    & )
!!$
!!$
!!$    real(DP), intent(out) :: TotalStrFlux
!!$
!!$
!!$    ! local variables
!!$    !
!!$    real(DP) :: DelWaveNum
!!$
!!$    integer  :: l
!!$
!!$
!!$    TotalStrFlux = 0.0d0
!!$    do l = 1, NBandSave
!!$      DelWaveNum = aa_BandWaveNumBndsSave(2,l) - aa_BandWaveNumBndsSave(1,l)
!!$      TotalStrFlux = TotalStrFlux &
!!$        & + b_BandAveStrPFSave( l ) * DelWaveNum
!!$    end do
!!$
!!$
!!$  end subroutine OptPropLBLConstPFGetTotalStrFlux

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

  subroutine OptPropLBLConstPFGetAbsCoefProf( &
    & iWaveNum, kmax, NMol, m_MolNum, pm_VMR, &
    & pm_AbsCoef &
    & )


    integer     , intent(in ) :: iWaveNum
    integer     , intent(in ) :: kmax
    integer     , intent(in ) :: NMol
    integer     , intent(in ) :: m_MolNum  (1:NMol)
    real(DP)    , intent(in ) :: pm_VMR    (1:kmax,1:NMol)
    real(DP)    , intent(out) :: pm_AbsCoef(1:kmax,1:NMol)


    ! local variables
    !
!!$    integer        :: k
    integer        :: m


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


    ! check
    if ( kmax /= NPressSave ) then
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Number of layers does not match with number of pressure levels in optical properties: ', &
        & kmax, NPressSave
      stop
    end if


    if ( NMol /= NMolSave ) then
      write( 6, * ) 'NMol in arguments does not match with NMol in absorption coefficient file.'
      write( 6, * ) NMol, NMolSave
      stop
    end if
    do m = 1, NMol
      if ( m_MolNum(m) /= m_MolNumSave(m) ) then
        write( 6, * ) 'Molecular number is inappropriate: ', &
          & m, m_MolNum(m), m_MolNumSave(m)
        stop
      end if
    end do


    pm_AbsCoef = zma_AbsCoef(:,:,iWaveNum)

    do m = 1, NMolSave
      select case ( ma_IDAbsCoefType(m,iWaveNum) )
      case ( -1 )
        pm_AbsCoef(:,m) = 0.0d0
      case ( 0 )
        pm_AbsCoef(:,m) = pm_AbsCoef(:,m) * pm_VMR(:,m)
!!$          case ( 1 )
!!$            pm_AbsCoef(:,m) = pm_AbsCoef(:,m)
      end select
    end do


  end subroutine OptPropLBLConstPFGetAbsCoefProf

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

!!$  subroutine OptPropLBLConstPFGetPFInted( &
!!$    & iWaveNum, &
!!$    & kmax, &
!!$    & r_Temp, SurfTemp, &
!!$    & r_AtmPFInted, SurfPFInted, SurfDPFDTInted &
!!$    & )
!!$
!!$    use constants0, only : PI
!!$
!!$    use planck_func_wrapper
!!$
!!$    use rad_kd_utils, only : InqPFIntGQNum
!!$
!!$
!!$    integer , intent(in ) :: iWaveNum
!!$    integer , intent(in ) :: kmax
!!$    real(DP), intent(in ) :: r_Temp   (1:kmax)
!!$    real(DP), intent(in ) :: SurfTemp
!!$    real(DP), intent(out) :: r_AtmPFInted(1:kmax)
!!$    real(DP), intent(out) :: SurfPFInted
!!$    real(DP), intent(out) :: SurfDPFDTInted
!!$
!!$
!!$    ! local variables
!!$    !
!!$    real(DP) :: WaveNum
!!$    real(DP) :: DelWaveNum
!!$    real(DP) :: WNs
!!$    real(DP) :: WNe
!!$    integer  :: GQNum
!!$
!!$
!!$    WaveNum    = OptPropLBLConstPFGetWN( iWaveNum )
!!$    DelWaveNum = OptPropLBLConstPFGetDelWN()
!!$    WNs = WaveNum - DelWaveNum / 2.0d0
!!$    WNe = WaveNum + DelWaveNum / 2.0d0
!!$
!!$!    GQNum = 1
!!$!    GQNum = int( DelWaveNum / 10.0d0 ) + 1
!!$    GQNum = InqPFIntGQNum( DelWaveNum )
!!$
!!$    call Integ_PF_GQ_Array1D( &
!!$      & WNs, WNe, GQNum, &
!!$      & 1, kmax, &
!!$      & r_Temp, &
!!$      & r_AtmPFInted &
!!$      & )
!!$    r_AtmPFInted = PI * r_AtmPFInted
!!$    call Integ_PF_GQ_Array0D( &
!!$      & WNs, WNe, GQNum, &
!!$      & SurfTemp, &
!!$      & SurfPFInted &
!!$      & )
!!$    SurfPFInted = PI * SurfPFInted
!!$    call Integ_DPFDT_GQ_Array0D( &
!!$      & WNs, WNe, GQNum,                 & ! (in )
!!$      & SurfTemp,                & ! (in )
!!$      & SurfDPFDTInted                  & ! (out)
!!$      & )
!!$    SurfDPFDTInted = PI * SurfDPFDTInted
!!$
!!$
!!$  end subroutine OptPropLBLConstPFGetPFInted

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

!!$  subroutine OptPropLBLConstPFGetStrPFInted( &
!!$    & iWaveNum, &
!!$    & StrPFInted &
!!$    & )
!!$
!!$    use constants0, only : PI
!!$
!!$    use planck_func_wrapper
!!$
!!$    use rad_kd_utils, only : InqPFIntGQNum
!!$
!!$    integer , intent(in ) :: iWaveNum
!!$    real(DP), intent(out) :: StrPFInted
!!$
!!$
!!$    ! local variables
!!$    !
!!$    real(DP) :: WaveNum
!!$    real(DP) :: DelWaveNum
!!$    integer  :: iBand
!!$
!!$
!!$    WaveNum = w_WaveNumSave( iWaveNum )
!!$    do iBand = 1, NBandSave
!!$      if ( WaveNum <= aa_BandWaveNumBndsSave(2,iBand) ) then
!!$        exit
!!$      end if
!!$    end do
!!$    StrPFInted = b_BandAveStrPFSave( iBand )
!!$    DelWaveNum = aa_BandWaveNumBndsSave(2,iBand) - aa_BandWaveNumBndsSave(1,iBand)
!!$    StrPFInted = StrPFInted * DelWaveNum
!!$
!!$
!!$  end subroutine OptPropLBLConstPFGetStrPFInted

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

  subroutine OptPropLBLConstPFGetStrBandAvePF( &
    & iBand, StrFluxTOA, &
    & StrBandAvePF &
    & )

    integer , intent(in ) :: iBand
    real(DP), intent(in ) :: StrFluxTOA
    real(DP), intent(out) :: StrBandAvePF


    ! local variables
    !


!!$    StrBandAvePF = b_BandAveStrPFSave( iBand )
    StrBandAvePF = b_RatioBandIntFluxtoTotFluxSave( iBand )
    StrBandAvePF = StrBandAvePF * StrFluxTOA
    StrBandAvePF = StrBandAvePF &
      & / ( aa_BandWaveNumBndsSave(2,iBand) - aa_BandWaveNumBndsSave(1,iBand) )


  end subroutine OptPropLBLConstPFGetStrBandAvePF

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

  subroutine OptPropLBLConstPFGetPFIntedRatio( &
    & iWaveNum, kmax, &
    & r_AtmPFRatio, SurfPFRatio &
    & )


    integer     , intent(in ) :: iWaveNum
    integer     , intent(in ) :: kmax
    real(DP)    , intent(out) :: r_AtmPFRatio(1:kmax)
    real(DP)    , intent(out) :: SurfPFRatio


    ! local variables
    !

!!$    real(DP), allocatable :: za_PFRatio(:,:)


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


    ! check
    if ( kmax /= NPressForPFSave ) then
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Number of layers does not match with number of pressure levels in optical properties: ', &
        & kmax, NPressForPFSave
      stop
    end if


    r_AtmPFRatio = ra_AtmPFRatio(:,iWaveNum)


!!$    allocate( za_PFRatio( NPressSave, NWaveNumSave ) )
!!$    call HistoryGet( OptPropNcFNSave, PFRatioVarName, za_PFRatio, &
!!$      & quiet = .true. )
!!$    z_PFRatio = za_PFRatio(:,iWaveNum)
!!$    deallocate( za_PFRatio )


    SurfPFRatio = r_AtmPFRatio(1)


  end subroutine OptPropLBLConstPFGetPFIntedRatio

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

  subroutine OptPropLBLConstPFGetStrPFIntedRatio( &
    & iWaveNum, &
    & StrPFRatio &
    & )


    integer     , intent(in ) :: iWaveNum
    real(DP)    , intent(out) :: StrPFRatio


    ! local variables
    !

!!$    real(DP), allocatable :: za_PFRatio(:,:)


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


    StrPFRatio = a_StrPFRatio(iWaveNum)


  end subroutine OptPropLBLConstPFGetStrPFIntedRatio

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

end module opt_prop_lblconstpf
