module opt_prop

  use dc_types
  use gtool_history

  implicit none

  private

  public :: OptPropInit
  public :: OptPropFinalize
  public :: OptPropGetNWN
  public :: OptPropGetNPress
  public :: OptPropGetNMol
  public :: OptPropGetNBand
  public :: OptPropGetDelWN
  public :: OptPropGetWN
  public :: OptPropGetBandWNBnds
  public :: OptPropGetBandBinIndexBnds
  public :: OptPropGetAbsCoefProf
  public :: OptPropGetPFInted
  public :: OptPropGetPFIntedRatio


  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 :: PressAxisName
  character(128)       , save :: MolAxisName
  !
  integer              , save :: NPressSave
  integer              , save :: NMolSave
  !
  real(DP), allocatable, save :: p_PressSave(:)
  real(DP), allocatable, save :: m_MolNumSave(:)

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

  ! KDDevelop, KD
  integer , parameter         :: NBndSave = 2
  integer              , save :: NBinSave
  integer              , save :: NBandSave
  integer , allocatable, save :: a_BinNumSave(:)
  real(DP), allocatable, save :: aa_BandWaveNumBndsSave     (:,:)
!!$  integer , allocatable, save :: aa_BandWaveNumIndexBndsSave(:,:)
  integer , allocatable, save :: aa_BandBinIndexBndsSave    (:,:)
  !
  character(128), save :: WaveNumAxisName
  character(128), save :: AbsCoefVarName
  character(128), save :: PFRatioVarName
  !
  real(DP), allocatable, save :: a_DelWaveNum(:)
  !
  ! KDDevelop only
  real(DP), allocatable, save :: za_PFRatio  (:,:)
  real(DP), allocatable, save :: zma_AbsCoef (:,:,:)
  !
  ! KD only
  integer              , save :: KDTblNPress
  integer              , save :: KDTblNTemp
  integer              , save :: KDTblNVMRH2O
  real(DP), allocatable, save :: p_KDTblPress (:)
  real(DP), allocatable, save :: t_KDTblTemp  (:)
  real(DP), allocatable, save :: v_KDTblVMRH2O(:)
  real(DP), allocatable, save :: ct_PFRatio   (:,:)
  real(DP), allocatable, save :: pmctv_AbsCoef(:,:,:,:,:)


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

  logical, save :: FlagInited = .false.

contains

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

  subroutine OptPropInit( &
    & Type, &
    & OptPropNcFN &
    & )

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

    if ( FlagInited ) return


    select case ( Type )
    case ( 'LBL' )
      IDMethod = IDMethodLBL
    case ( 'LBLConstPF' )
      IDMethod = IDMethodLBLConstPF
    case ( 'KDDevelop' )
      IDMethod = IDMethodKDDevelop
    case ( 'KD' )
      IDMethod = IDMethodKD
    case default
      write( 6, * ) 'In module, ', trim( ModuleName )
      write( 6, * ) 'Unexpected type: ', trim( Type )
      stop
    end select


    call OptPropInitCore( &
      & OptPropNcFN &
      & )


    FlagInited = .true.

  end subroutine OptPropInit

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

  subroutine OptPropFinalize


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


    call OptPropFinalizeCore


    FlagInited = .false.

  end subroutine OptPropFinalize

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

  function OptPropGetNWN() result( NWN )


    integer :: NWN

    ! local variables
    !


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


    NWN = OptPropGetNWNCore()


  end function OptPropGetNWN

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

  function OptPropGetNPress() result( NPress )


    integer :: NPress

    ! local variables
    !


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


    NPress = OptPropGetNPressCore()


  end function OptPropGetNPress

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

  function OptPropGetNMol() result( NMol )


    integer :: NMol

    ! local variables
    !


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


    NMol = OptPropGetNMolCore()


  end function OptPropGetNMol

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

  function OptPropGetNBand() result( NBand )


    integer :: NBand

    ! local variables
    !


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


    NBand = OptPropGetNBandCore()


  end function OptPropGetNBand

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

  function OptPropGetDelWN( iWaveNum ) result( DelWaveNum )


    integer, intent(in) :: iWaveNum


    real(DP) :: DelWaveNum

    ! local variables
    !


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


    DelWaveNum = OptPropGetDelWNCore( iWaveNum )


  end function OptPropGetDelWN

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

  function OptPropGetWN( 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 = OptPropGetWNCore( iWaveNum )


  end function OptPropGetWN

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

  function OptPropGetBandWNBnds() result( aa_BandWNBnds )


    real(DP) :: aa_BandWNBnds(2,NBandSave)


    ! local variables
    !


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


    aa_BandWNBnds = OptPropGetBandWNBndsCore()


  end function OptPropGetBandWNBnds

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

  function OptPropGetBandBinIndexBnds() result( aa_BandBinIndexBnds )


    integer :: aa_BandBinIndexBnds(2,NBandSave)


    ! local variables
    !


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


    aa_BandBinIndexBnds = OptPropGetBandBinIndexBndsCore()


  end function OptPropGetBandBinIndexBnds

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

  subroutine OptPropGetAbsCoefProf( &
    & iWaveNum, kmax, NMol, r_Press, r_Temp, r_VMRH2O, &
    & rm_AbsCoef &
    & )


    integer     , intent(in ) :: iWaveNum
    integer     , intent(in ) :: kmax
    integer     , intent(in ) :: NMol
    real(DP)    , intent(in ) :: r_Press   (1:kmax)
    real(DP)    , intent(in ) :: r_Temp    (1:kmax)
    real(DP)    , intent(in ) :: r_VMRH2O  (1:kmax)
    real(DP)    , intent(out) :: rm_AbsCoef(1:kmax,1:NMol)


    ! local variables
    !


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


    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF, IDMethodKDDevelop )
      call OptPropGetAbsCoefProfCore( &
        & iWaveNum, kmax, NMol, r_Press, &
        & rm_AbsCoef &
        & )
    case ( IDMethodKD )
      call OptPropGetAbsCoefProfKDCore( &
        & iWaveNum, kmax, NMol, r_Press, r_Temp, r_VMRH2O, &
        & rm_AbsCoef &
        & )
    case default
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Unexpected case default'
      stop
    end select


  end subroutine OptPropGetAbsCoefProf

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

  subroutine OptPropGetPFInted( &
    & iWaveNum, &
    & kmax, &
    & r_Temp, SurfTemp, &
    & r_PFInted, SurfPFInted, SurfDPFDTInted &
    & )

    integer , intent(in ) :: iWaveNum
    integer , intent(in ) :: kmax
    real(DP), intent(in ) :: r_Temp(0:kmax)
    real(DP), intent(in ) :: SurfTemp
    real(DP), intent(out) :: r_PFInted  (0:kmax)
    real(DP), intent(out) :: SurfPFInted
    real(DP), intent(out) :: SurfDPFDTInted


    ! local variables
    !


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


    call OptPropGetPFIntedCore( &
      & iWaveNum, &
      & kmax, &
      & r_Temp, SurfTemp, &
      & r_PFInted, SurfPFInted, SurfDPFDTInted &
      & )


  end subroutine OptPropGetPFInted

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

  subroutine OptPropGetPFIntedRatio( &
    & iWaveNum, kmax, r_Temp, SurfTemp, &
    & r_PFRatio, SurfPFRatio &
    & )

    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_PFRatio(1:kmax)
    real(DP)    , intent(out) :: SurfPFRatio


    ! local variables
    !


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


    select case ( IDMethod )
    case ( IDMethodLBL )
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  function OptPropGetPFIntedRatio cannot be used when IDMethod is IDMethodLBL'
      stop
    end select


    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF, IDMethodKDDevelop )
      call OptPropGetPFIntedRatioCore( &
        & iWaveNum, kmax, &
        & r_PFRatio, SurfPFRatio &
        & )
    case ( IDMethodKD )
      call OptPropGetPFIntedRatioKDCore( &
        & iWaveNum, kmax, r_Temp, SurfTemp, &
        & r_PFRatio, SurfPFRatio &
        & )
    case default
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Unexpected case default'
      stop
    end select


  end subroutine OptPropGetPFIntedRatio

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

  subroutine OptPropInitCore( &
    & OptPropNcFN &
    & )

    use rad_kd_utils, only : InqDimLength


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


    ! local variables
    !


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


    NPressSave   = -1
    NMolSave     = -1
    NWaveNumSave = -1
    NBinSave     = -1
    NBandSave    = -1

    KDTblNPress  = -1
    KDTblNTemp   = -1
    KDTblNVMRH2O = -1


    OptPropNcFNSave = OptPropNcFN

    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF, IDMethodKDDevelop )
      PressAxisName = "r_Press"
      MolAxisName   = "m_MolNum"
    case ( IDMethodKD )
      PressAxisName = "Press"
      MolAxisName   = "MolNum"
    end select
    NPressSave = InqDimLength( OptPropNcFNSave, PressAxisName )
    NMolSave   = InqDimLength( OptPropNcFNSave, MolAxisName   )
    allocate( p_PressSave (1:NPressSave) )
    allocate( m_MolNumSave(1:NMolSave  ) )
    call HistoryGet( OptPropNcFNSave, PressAxisName, p_PressSave , quiet = .true. )
    call HistoryGet( OptPropNcFNSave, MolAxisName  , m_MolNumSave, quiet = .true. )


    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF )
      WaveNumAxisName = 'w_WaveNum'
      AbsCoefVarName  = 'rmw_TAbsorpK'
      PFRatioVarName  = 'rw_PFRatio'
      !
      NWaveNumSave = InqDimLength( OptPropNcFNSave, WaveNumAxisName )
      allocate( w_WaveNumSave( 1:NWaveNumSave ) )
      call HistoryGet( OptPropNcFNSave, WaveNumAxisName, w_WaveNumSave, quiet = .true. )
    case ( IDMethodKDDevelop, IDMethodKD )
      WaveNumAxisName = 'BinNum'
      !
      NBinSave    = InqDimLength( OptPropNcFNSave, WaveNumAxisName )
      NBandSave   = InqDimLength( OptPropNcFNSave, 'BandNum' )
      allocate( a_BinNumSave               (NBinSave) )
      allocate( aa_BandWaveNumBndsSave     (NBndSave,NBandSave) )
!!$      allocate( aa_BandWaveNumIndexBndsSave(NBndSave,NBandSave) )
      allocate( aa_BandBinIndexBndsSave    (NBndSave,NBandSave) )

      call HistoryGet( OptPropNcFNSave, WaveNumAxisName       , a_BinNumSave, quiet = .true. )
      call HistoryGet( OptPropNcFNSave, 'BandWaveNumBnds'     , aa_BandWaveNumBndsSave, quiet = .true. )
!!$      call HistoryGet( OptPropNcFNSave, 'BandWaveNumIndexBnds', aa_BandWaveNumIndexBndsSave, quiet = .true. )
      call HistoryGet( OptPropNcFNSave, 'BandBinIndexBnds'    , aa_BandBinIndexBndsSave, quiet = .true. )
      !
      select case ( IDMethod )
      case ( IDMethodKDDevelop )
        AbsCoefVarName  = 'rma_TAbsorpK'
        PFRatioVarName  = 'ra_PFRatio'
      case ( IDMethodKD )
        AbsCoefVarName  = 'AbsCoef'
        PFRatioVarName  = 'PFRatio'
      end select
    end select

    select case ( IDMethod )
    case ( IDMethodKDDevelop, IDMethodKD )
      allocate( a_DelWaveNum(NBinSave) )
      call HistoryGet( OptPropNcFNSave, "BinWeight", a_DelWaveNum, &
        & quiet = .true. )
    end select
    !
    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF )
      allocate( zma_AbsCoef (NPressSave,NMolSave,NWaveNumSave) )
      call HistoryGet( OptPropNcFNSave, AbsCoefVarName, zma_AbsCoef, &
        & quiet = .true. )
    case ( IDMethodKDDevelop )
      allocate( zma_AbsCoef (NPressSave,NMolSave,NBinSave) )
      call HistoryGet( OptPropNcFNSave, AbsCoefVarName, zma_AbsCoef, &
        & quiet = .true. )
    case ( IDMethodKD )
      KDTblNPress  = NPressSave
      KDTblNTemp   = InqDimLength( OptPropNcFNSave, 'Temp'   )
      KDTblNVMRH2O = InqDimLength( OptPropNcFNSave, 'VMRH2O' )
      allocate( p_KDTblPress (KDTblNPress ) )
      allocate( t_KDTblTemp  (KDTblNTemp  ) )
      allocate( v_KDTblVMRH2O(KDTblNVMRH2O) )
      allocate( pmctv_AbsCoef(KDTblNPress,NMolSave,NBinSave,KDTblNTemp,KDTblNVMRH2O) )
      call HistoryGet( OptPropNcFNSave, 'Press' , p_KDTblPress, &
        & quiet = .true. )
      call HistoryGet( OptPropNcFNSave, 'Temp'  , t_KDTblTemp, &
        & quiet = .true. )
      call HistoryGet( OptPropNcFNSave, 'VMRH2O', v_KDTblVMRH2O, &
        & quiet = .true. )
      call HistoryGet( OptPropNcFNSave, AbsCoefVarName, pmctv_AbsCoef, &
        & quiet = .true. )
      !
      select case ( IDMethod )
      case ( IDMethodLBL )
      case ( IDMethodLBLConstPF )
        allocate( za_PFRatio(NPressSave,NWaveNumSave) )
        call HistoryGet( OptPropNcFNSave, PFRatioVarName, za_PFRatio, &
          & quiet = .true. )
      case ( IDMethodKDDevelop )
        allocate( za_PFRatio(NPressSave,NBinSave) )
        call HistoryGet( OptPropNcFNSave, PFRatioVarName, za_PFRatio, &
          & quiet = .true. )
      case ( IDMethodKD )
        allocate( ct_PFRatio(NBinSave,KDTblNTemp) )
        call HistoryGet( OptPropNcFNSave, PFRatioVarName, ct_PFRatio, &
          & quiet = .true. )
      end select
    end select


  end subroutine OptPropInitCore

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

  subroutine OptPropFinalizeCore


    ! local variables
    !


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


    NPressSave   = -1
    NMolSave     = -1
    NWaveNumSave = -1
    NBinSave     = -1
    NBandSave    = -1

    KDTblNPress  = -1
    KDTblNTemp   = -1
    KDTblNVMRH2O = -1


    deallocate( p_PressSave  )
    deallocate( m_MolNumSave )

    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF )
      NWaveNumSave    = 0
      deallocate( w_WaveNumSave )
    case ( IDMethodKDDevelop, IDMethodKD )
      NBinSave    = 0
      deallocate( a_BinNumSave                )
      deallocate( aa_BandWaveNumBndsSave      )
!!$      deallocate( aa_BandWaveNumIndexBndsSave )
      deallocate( aa_BandBinIndexBndsSave     )
    end select
    !
    select case ( IDMethod )
    case ( IDMethodKDDevelop, IDMethodKD )
      deallocate( a_DelWaveNum )
    end select
    !
    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF )
      deallocate( zma_AbsCoef )
    case ( IDMethodKDDevelop )
      deallocate( zma_AbsCoef )
    case ( IDMethodKD )
      deallocate( p_KDTblPress  )
      deallocate( t_KDTblTemp   )
      deallocate( v_KDTblVMRH2O )
      deallocate( pmctv_AbsCoef )
      !
      select case ( IDMethod )
      case ( IDMethodLBL )
      case ( IDMethodLBLConstPF )
        deallocate( za_PFRatio )
      case ( IDMethodKDDevelop )
        deallocate( za_PFRatio )
      case ( IDMethodKD )
        deallocate( ct_PFRatio )
      end select
    end select


  end subroutine OptPropFinalizeCore

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

  function OptPropGetNWNCore() result( NWN )


    integer :: NWN

    ! local variables
    !


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


    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF )
      NWN = NWaveNumSave
    case ( IDMethodKDDevelop, IDMethodKD )
      NWN = NBinSave
    case default
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Unexpected case default'
      stop
    end select

  end function OptPropGetNWNCore

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

  function OptPropGetNPressCore() result( NPress )


    integer :: NPress

    ! local variables
    !


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


    NPress = NPressSave


  end function OptPropGetNPressCore

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

  function OptPropGetNMolCore() result( NMol )


    integer :: NMol

    ! local variables
    !


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


    NMol = NMolSave


  end function OptPropGetNMolCore

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

  function OptPropGetNBandCore() result( NBand )


    integer :: NBand

    ! local variables
    !


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


    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF )
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  function OptPropGetNBandCore cannot be used when IDMethod is IDMethodLBL or IDMethodLBLConstPF'
      stop
    case ( IDMethodKDDevelop, IDMethodKD )
      NBand = NBandSave
    case default
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Unexpected case default'
      stop
    end select


  end function OptPropGetNBandCore

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

  function OptPropGetDelWNCore( iWaveNum ) result( DelWaveNum )

    use rad_kd_utils, only : MkRangeString

    integer, intent(in) :: iWaveNum

    real(DP) :: DelWaveNum

    ! local variables
    !
    character(128) :: String


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


    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF )
      DelWaveNum = w_WaveNumSave( 2 ) - w_WaveNumSave( 1 )
    case ( IDMethodKDDevelop, IDMethodKD )
      DelWaveNum = a_DelWaveNum(iWaveNum)
    end select


  end function OptPropGetDelWNCore

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

  function OptPropGetWNCore( 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


    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF )
      WaveNum = w_WaveNumSave( iWaveNum )
    case ( IDMethodKDDevelop, IDMethodKD )
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  function OptPropGetWNCore cannot be used when IDMethod is IDMethodKDDevelop'
      stop
    case default
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Unexpected case default'
      stop
    end select


  end function OptPropGetWNCore

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

  function OptPropGetBandWNBndsCore( ) result( aa_BandWNBnds )


    real(DP) :: aa_BandWNBnds(2,NBandSave)

    ! local variables
    !


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


    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF )
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  function OptPropGetBandWNBndsCore cannot be used when IDMethod is IDMethodLBL'
      stop
    case ( IDMethodKDDevelop, IDMethodKD )
      aa_BandWNBnds = aa_BandWaveNumBndsSave
    case default
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Unexpected case default'
      stop
    end select


  end function OptPropGetBandWNBndsCore

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

  function OptPropGetBandBinIndexBndsCore( ) result( aa_BandBinIndexBnds )


    integer :: aa_BandBinIndexBnds(2,NBandSave)


    ! local variables
    !


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


    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF )
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  function OptPropGetBandBinIndexBndsCore cannot be used when IDMethod is IDMethodLBL or IDMethodLBLConstPF'
      stop
    case ( IDMethodKDDevelop, IDMethodKD )
      aa_BandBinIndexBnds = aa_BandBinIndexBndsSave
    end select


  end function OptPropGetBandBinIndexBndsCore

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

  subroutine OptPropGetAbsCoefProfCore( &
    & iWaveNum, kmax, NMol, r_Press, &
    & pm_AbsCoef &
    & )

    use rad_kd_utils, only : MkRangeString


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


    ! local variables
    !
    character(128) :: String
!!$    integer        :: k


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


    ! check
    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF, IDMethodKDDevelop )
    case default
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Unexpected case default'
      stop
    end select


    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF, IDMethodKDDevelop )
      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
!!$      do k = 1, kmax
!!$        if ( abs( p_Press(k) - p_PressSave(k) ) > p_PressSave(k) * 1.0d-6 ) then
!!$          write( 6, * ) 'In ', trim( ModuleName )
!!$          write( 6, * ) '  Pressure is inappropriate: ', p_PressSave(k), p_Press(k)
!!$          stop
!!$        end if
!!$      end do
    end select


    pm_AbsCoef = zma_AbsCoef(:,:,iWaveNum)


  end subroutine OptPropGetAbsCoefProfCore

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

  subroutine OptPropGetAbsCoefProfKDCore( &
    & iWaveNum, kmax, NMol, r_Press, r_Temp, r_VMRH2O, &
    & rm_AbsCoef &
    & )


    integer     , intent(in ) :: iWaveNum
    integer     , intent(in ) :: kmax
    integer     , intent(in ) :: NMol
    real(DP)    , intent(in ) :: r_Press   (1:kmax)
    real(DP)    , intent(in ) :: r_Temp    (1:kmax)
    real(DP)    , intent(in ) :: r_VMRH2O  (1:kmax)
    real(DP)    , intent(out) :: rm_AbsCoef(1:kmax,1:NMol)


    ! local variables
    !
    integer        :: r_KDTblPressIndex (1:kmax)
    integer        :: r_KDTblTempIndex  (1:kmax)
    integer        :: r_KDTblVMRH2OIndex(1:kmax)
    real(DP)       :: a_Grid   (2)
    real(DP)       :: a_Weight1(2)
    real(DP)       :: a_Weight2(2)
    real(DP)       :: a_Weight3(2)
    real(DP)       :: a_Grid1(2)
    real(DP)       :: a_Grid2(2)
    real(DP)       :: a_Grid3(2)
    real(DP)       :: aaa_Weight(2,2,2)
    real(DP)       :: aaam_Val  (2,2,2,NMol)
    real(DP)       :: Val
    real(DP)       :: Press
    real(DP)       :: Temp
    real(DP)       :: VMRH2O
    integer        :: iPress
    integer        :: iTemp
    integer        :: iVMRH2O
    integer        :: k
    integer        :: l
    integer        :: m
    integer        :: n
    integer        :: l1
    integer        :: l2
    integer        :: l3
    integer        :: n1
    integer        :: n2
    integer        :: n3


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


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



    call FindTblIndex( &
      & 'Press', &
      & KDTblNPress, p_KDTblPress, kmax, r_Press, &
      & r_KDTblPressIndex &
      & )
    call FindTblIndex( &
      & 'Temp', &
      & KDTblNTemp, t_KDTblTemp, kmax, r_Temp, &
      & r_KDTblTempIndex &
      & )
    call FindTblIndex( &
      & 'VMRH2O', &
      & KDTblNVMRH2O, v_KDTblVMRH2O, kmax, r_VMRH2O, &
      & r_KDTblVMRH2OIndex &
      & )

    do k = 1, kmax
      Press  = r_Press (k)
      Temp   = r_Temp  (k)
      VMRH2O = r_VMRH2O(k)

      iPress  = r_KDTblPressIndex(k)
      iTemp   = r_KDTblTempIndex(k)
      iVMRH2O = r_KDTblVMRH2OIndex(k)


      ! l=3-n
      !   when n=1, l=2
      !   when n=2, l=1
      !
      if ( iPress == 2 ) then
        a_Weight1(1) = 1.0d0
        a_Weight1(2) = 0.0d0
      else
        a_Grid(1) = p_KDTblPress (iPress-1 )
        a_Grid(2) = p_KDTblPress (iPress   )
        do n = 1, 2
          l = 3-n
          a_Weight1(n) = ( Press    - a_Grid(l) ) / ( a_Grid(2) - a_Grid(1) )
        end do
      end if
      !
      a_Grid(1) = t_KDTblTemp  (iTemp-1  )
      a_Grid(2) = t_KDTblTemp  (iTemp    )
      do n = 1, 2
        l = 3-n
        a_Weight2(n) = ( Temp     - a_Grid(l) ) / ( a_Grid(2) - a_Grid(1) )
      end do
      !
      if ( iVMRH2O == 2 ) then
        a_Weight3(1) = 1.0d0
        a_Weight3(2) = 0.0d0
      else
        a_Grid(1) = v_KDTblVMRH2O(iVMRH2O-1)
        a_Grid(2) = v_KDTblVMRH2O(iVMRH2O  )
        do n = 1, 2
          l = 3-n
          a_Weight3(n) = ( VMRH2O   - a_Grid(l) ) / ( a_Grid(2) - a_Grid(1) )
        end do
      end if
      !
      aaa_Weight = 1.0d0
      do n3 = 1, 2
        do n2 = 1, 2
          do n1 = 1, 2
            aaa_Weight(n1,n2,n3) = aaa_Weight(n1,n2,n3) &
              & * a_Weight1(n1) * a_Weight2(n2) * a_Weight3(n3)
          end do
        end do
      end do
      aaa_Weight = abs( aaa_Weight )


!!$      aaa_Weight = 1.0d0
!!$      do n3 = 1, 2
!!$        do n2 = 1, 2
!!$          do n1 = 1, 2
!!$            ! when n1=1, l1=2
!!$            ! when n1=2, l1=1
!!$            l1 = 3-n1
!!$            l2 = 3-n2
!!$            l3 = 3-n3
!!$            aaa_Weight(n1,n2,n3) = aaa_Weight(n1,n2,n3) &
!!$              & * ( Press    - a_Grid1(l1) ) / ( a_Grid1(2) - a_Grid1(1) ) &
!!$              & * ( Temp     - a_Grid2(l2) ) / ( a_Grid2(2) - a_Grid2(1) ) &
!!$              & * ( VMRH2O   - a_Grid3(l3) ) / ( a_Grid3(2) - a_Grid3(1) )
!!$          end do
!!$        end do
!!$      end do
!!$      aaa_Weight = abs( aaa_Weight )

      do m = 1, NMol
        do n3 = 1, 2
          do n2 = 1, 2
            do n1 = 1, 2
              l1 = iPress -2+n1
              l2 = iTemp  -2+n2
              l3 = iVMRH2O-2+n3
              aaam_Val(n1,n2,n3,m) = pmctv_AbsCoef(l1,m,iWaveNum,l2,l3)
            end do
          end do
        end do
      end do


      do m = 1, NMol

        Val = 0.0d0
        do n3 = 1, 2
          do n2 = 1, 2
            do n1 = 1, 2
              Val = Val + aaam_Val(n1,n2,n3,m) * aaa_Weight(n1,n2,n3)
            end do
          end do
        end do


!!$        aaa_Val(1,1,1) = pmctv_AbsCoef(iPress1,iMol,iWaveNum,iTemp1,iVMRH2O1)
!!$        aaa_Val(2,1,1) = pmctv_AbsCoef(iPress2,iMol,iWaveNum,iTemp1,iVMRH2O1)
!!$        aaa_Val(1,2,1) = pmctv_AbsCoef(iPress1,iMol,iWaveNum,iTemp2,iVMRH2O1)
!!$        aaa_Val(2,2,1) = pmctv_AbsCoef(iPress2,iMol,iWaveNum,iTemp2,iVMRH2O1)
!!$        aaa_Val(1,1,2) = pmctv_AbsCoef(iPress1,iMol,iWaveNum,iTemp1,iVMRH2O2)
!!$        aaa_Val(2,1,2) = pmctv_AbsCoef(iPress2,iMol,iWaveNum,iTemp1,iVMRH2O2)
!!$        aaa_Val(1,2,2) = pmctv_AbsCoef(iPress1,iMol,iWaveNum,iTemp2,iVMRH2O2)
!!$        aaa_Val(2,2,2) = pmctv_AbsCoef(iPress2,iMol,iWaveNum,iTemp2,iVMRH2O2)
!!$
!!$        do n2 = 1, 2
!!$          do n1 = 1, 2
!!$            aa_Val(n1,n2) = &
!!$              &   ( aaa_Val(n1,n2,1) - aaa_Val(n1,n2,2) ) &
!!$              & / ( a_Grid3(1)       - a_Grid3(2)     ) &
!!$              & * ( VMRH2O           - a_Grid3(2)     ) + aaa_Val(n1,n2,2)
!!$          end do
!!$        end do
!!$        do n1 = 1, 2
!!$          a_Val(n1) = &
!!$            &   ( aa_Val(n1,1) - aa_Val(n1,2) ) &
!!$            & / ( a_Grid2(1)   - a_Grid2(2)     ) &
!!$            & * ( Temp         - a_Grid2(2)     ) + aa_Val(n1,2)
!!$        end do
!!$        Val = &
!!$          &   ( a_Val(1)   - a_Val(2)   ) &
!!$          & / ( a_Grid1(1) - a_Grid1(2) ) &
!!$          & * ( Press      - a_Grid1(2) ) + a_Val(2)

        rm_AbsCoef(k,m) = Val
      end do
    end do


  end subroutine OptPropGetAbsCoefProfKDCore

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

  subroutine FindTblIndex( &
    & AxisName, &
    & NTblGrid, a_TblGrid, NGrid, a_Grid, &
    & a_Index &
    & )


    character(*), intent(in ) :: AxisName
    integer     , intent(in ) :: NTblGrid
    real(DP)    , intent(in ) :: a_TblGrid (1:NTblGrid)
    integer     , intent(in ) :: NGrid
    real(DP)    , intent(in ) :: a_Grid    (1:NGrid)
    integer     , intent(out) :: a_Index   (1:NGrid)


    ! local variables
    !
    integer :: i
    integer :: j


    if ( a_TblGrid(2) > a_TblGrid(1) ) then

      do i = 1, NGrid
        if ( ( a_Grid(i)           < a_TblGrid(1) ) .or. &
          &  ( a_TblGrid(NTblGrid) < a_Grid(i)    ) ) then
          write( 6, * ) trim( AxisName ), " grid is out of range of a table."
          write( 6, * ) a_TblGrid(1), a_Grid(i), a_TblGrid(NTblGrid)
          stop
        end if
        do j = 1+1, NTblGrid
          if ( a_TblGrid(j) >= a_Grid(i) ) exit
        end do
        a_Index(i) = j
      end do

    else

      do i = 1, NGrid
        if ( ( a_Grid(i)    < a_TblGrid(NTblGrid) ) .or. &
          &  ( a_TblGrid(1) < a_Grid(i)           ) ) then
          write( 6, * ) trim( AxisName ), " grid is out of range of a table."
          write( 6, * ) a_TblGrid(NTblGrid), a_Grid(i), a_TblGrid(1)
          stop
        end if
        do j = 1+1, NTblGrid
          if ( a_TblGrid(j) <= a_Grid(i) ) exit
        end do
        a_Index(i) = j
      end do

    end if


  end subroutine FindTblIndex

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

  subroutine OptPropGetPFIntedCore( &
    & iWaveNum, &
    & kmax, &
    & r_Temp, SurfTemp, &
    & r_PFInted, SurfPFInted, SurfDPFDTInted &
    & )

    use constants0, only : PI

    use planck_func_wrapper

    integer , intent(in ) :: iWaveNum
    integer , intent(in ) :: kmax
    real(DP), intent(in ) :: r_Temp(0:kmax)
    real(DP), intent(in ) :: SurfTemp
    real(DP), intent(out) :: r_PFInted  (0: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


    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF )
      WaveNum    = OptPropGetWNCore( iWaveNum )
      DelWaveNum = OptPropGetDelWNCore( iWaveNum )
      WNs = WaveNum - DelWaveNum / 2.0d0
      WNe = WaveNum + DelWaveNum / 2.0d0
    case ( IDMethodKDDevelop, IDMethodKD )
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  function OptPropGetPFIntedCore cannot be used when IDMethod is IDMethodKDDevelop, IDMethodKD'
      stop
    case default
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Unexpected case default'
      stop
    end select

!!$    GQNum = 1
    GQNum = int( DelWaveNum / 10.0d0 ) + 1

    call Integ_PF_GQ_Array1D( &
      & WNs, WNe, GQNum, &
      & 0, kmax, &
      & r_Temp, &
      & r_PFInted &
      & )
    r_PFInted = PI * r_PFInted
    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 OptPropGetPFIntedCore

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

  subroutine OptPropGetPFIntedRatioCore( &
    & iWaveNum, kmax, &
    & r_PFRatio, SurfPFRatio &
    & )

    use rad_kd_utils, only : MkRangeString


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


    ! local variables
    !
    character(128) :: String

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


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


    ! check
    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF, IDMethodKDDevelop )
    case default
      write( 6, * ) 'In ', trim( ModuleName )
      write( 6, * ) '  Unexpected case default'
      stop
    end select
    select case ( IDMethod )
    case ( IDMethodLBL, IDMethodLBLConstPF, IDMethodKDDevelop )
      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
    end select


    r_PFRatio = za_PFRatio(:,iWaveNum)


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


    SurfPFRatio = r_PFRatio(1)


  end subroutine OptPropGetPFIntedRatioCore

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

  subroutine OptPropGetPFIntedRatioKDCore( &
    & iWaveNum, kmax, r_Temp, SurfTemp, &
    & r_PFRatio, SurfPFRatio &
    & )


    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_PFRatio(1:kmax)
    real(DP)    , intent(out) :: SurfPFRatio


    ! local variables
    !
    integer        :: r_KDTblTempIndex  (1:kmax)
    real(DP)       :: a_Grid(2)
    real(DP)       :: a_Val (2)
    real(DP)       :: Val
    real(DP)       :: a_Weight(2)
    real(DP)       :: a_SurfTemp          (1)
    integer        :: a_KDTblSurfTempIndex(1)
    real(DP)       :: Temp
    integer        :: iTemp
    integer        :: k
    integer        :: l
    integer        :: n


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


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


    call FindTblIndex( &
      & 'Temp', &
      & KDTblNTemp, t_KDTblTemp, kmax, r_Temp, &
      & r_KDTblTempIndex &
      & )
    a_SurfTemp = SurfTemp
    call FindTblIndex( &
      & 'SurfTemp', &
      & KDTblNTemp, t_KDTblTemp, 1, a_SurfTemp, &
      & a_KDTblSurfTempIndex &
      & )

    do k = 0, kmax
      if ( k == 0 ) then
        Temp   = a_SurfTemp          (1)
        iTemp  = a_KDTblSurfTempIndex(1)
      else
        Temp   = r_Temp          (k)
        iTemp  = r_KDTblTempIndex(k)
      end if

      a_Grid(1) = t_KDTblTemp  (iTemp-1  )
      a_Grid(2) = t_KDTblTemp  (iTemp    )

      ! l = 3-n
      !   when n=1, l=2
      !   when n=2, l=1
      !
      do n = 1, 2
        l = 3-n
        a_Weight(n) = ( Temp     - a_Grid(l) ) / ( a_Grid(2) - a_Grid(1) )
      end do
      a_Weight = abs( a_Weight )

      do n = 1, 2
        l = iTemp  -2+n
        a_Val(n) = ct_PFRatio(iWaveNum,l)
      end do

      Val = 0.0d0
      do n = 1, 2
        Val = Val + a_Val(n) * a_Weight(n)
      end do

      if ( k == 0 ) then
        SurfPFRatio  = Val
      else
        r_PFRatio(k) = Val
      end if
    end do


  end subroutine OptPropGetPFIntedRatioKDCore

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

end module opt_prop

