module opt_prop_kd

  use vtype_module

  implicit none

  private

  public :: OptPropKDInit
  public :: OptPropKDFinalize
  public :: OptPropKDGetNWN
  public :: OptPropKDGetNPress
  public :: OptPropKDGetNMol
  public :: OptPropKDGetNBand
  public :: OptPropKDGetDelWN
  public :: OptPropKDGetBandNum
  public :: OptPropKDGetBandWNBnds
  public :: OptPropKDGetBandBinIndexBnds
  public :: OptPropKDGetStrBandAvePF
!!$  public :: OptPropKDGetTotalStrFlux
  public :: OptPropKDCalcIndexAndWeight
  public :: OptPropKDGetAbsCoefProf
  public :: OptPropKDGetAbsCoefProfWithIndex
  public :: OptPropKDGetRayScatCoef
  public :: OptPropKDGetPtclParam
  public :: OptPropKDGetPFIntedRatio
  public :: OptPropKDGetPFIntedRatioWithIndex
  public :: OptPropKDGetStrPFIntedRatio
  public :: OptPropKDPFTblGetPFDPFDT
  public :: OptPropKDPFTblChkBandIntPFRatio
  public :: OptPropKDPFTblChkBandIntStrFluxRatio

  interface OptPropKDPFTblGetPFDPFDT
    module procedure   &
      OptPropKDPFTblGetPFDPFDT2DV2, &
      OptPropKDPFTblGetPFDPFDT1DV2, &
      OptPropKDPFTblGetPFDPFDT1D, &
      OptPropKDPFTblGetPFDPFDT0D
  end interface OptPropKDPFTblGetPFDPFDT


  character(128)       , save :: OptPropNcFNSave

  integer , parameter         :: NBnd = 2

  real(DP), allocatable, save :: ab_KDTblBandWaveNumBnds    (:,:)
  integer , allocatable, save :: ab_KDTblBandBinIndexBnds   (:,:)
  !
  real(DP), allocatable, save :: c_DelWaveNum(:)
  !
  integer              , save :: KDTblIndexVMRMol1Num
  integer              , save :: KDTblIndexVMRMol2Num
  integer              , save :: KDTblIndexVMRMol3Num
  integer              , save :: KDTblNBin
  integer              , save :: KDTblNBand
  integer              , save :: KDTblNPress
  integer              , save :: KDTblNMol
  integer              , save :: KDTblNTemp
  integer              , save :: KDTblNVMRMol1
  integer              , save :: KDTblNVMRMol2
  integer              , save :: KDTblNVMRMol3
  integer , allocatable, save :: c_KDTblBinNum (:)
  real(DP), allocatable, save :: p_KDTblPress  (:)
  integer , allocatable, save :: m_KDTblMolNum (:)
  real(DP), allocatable, save :: t_KDTblTemp   (:)
  real(DP), allocatable, save :: u_KDTblVMRMol1(:)
  real(DP), allocatable, save :: v_KDTblVMRMol2(:)
  real(DP), allocatable, save :: w_KDTblVMRMol3(:)
!!$  real(DP), allocatable, save :: ct_AtmPFRatio    (:,:)
  real(DP), allocatable, save :: pctuvw_AtmPFRatio(:,:,:,:,:,:)
  real(DP), allocatable, save :: c_StrPFRatio     (:)
!!$  real(DP), allocatable, save :: b_BandAveStrPF   (:)
  real(DP), allocatable, save :: b_BandIntFluxtoTotFluxRatio(:)
!!$  real(DP), allocatable, save :: pmctuvw_AbsCoef(:,:,:,:,:,:,:)
  real(DP), allocatable, save :: pmctuvw_LogAbsCoef(:,:,:,:,:,:,:)
  real(DP), allocatable, save :: b_BandAveRayScatCoefNonRadAct(:)
  real(DP), allocatable, save :: bm_BandAveRayScatCoef(:,:)

!!$  integer                        :: NChar
!!$  integer          , allocatable :: a_CharNum(:)
  integer                       , save :: KDTblNPtcl
  integer                       , save :: KDTblPtclNRadiusMax
  integer          , allocatable, save :: a_KDTblPtclNum      (:)
  integer          , allocatable, save :: r_KDTblPtclRadiusNum(:)
  character(stdstr), allocatable, save :: a_KDTblPtclName     (:)
  integer          , allocatable, save :: a_KDTblPtclNRadius  (:)
  real(DP)         , allocatable, save :: ra_KDTblPtclRadius  (:,:)
  real(DP)         , allocatable, save :: bra_KDTblBandAveQExt(:,:,:)
  real(DP)         , allocatable, save :: bra_KDTblBandAveSSA (:,:,:)
  real(DP)         , allocatable, save :: bra_KDTblBandAveAF  (:,:,:)

  integer , allocatable, save :: mc_IDAbsCoefType(:,:)


  integer              , save :: PFTblNTemp
  integer              , save :: PFTblNBand
  real(DP), allocatable, save :: t_PFTblTemp    (:)
!!$  real(DP), allocatable, save :: b_PFTblBandNum(:)
  real(DP), allocatable, save :: tb_PFTblPF     (:,:)
  real(DP), allocatable, save :: tb_PFTblDPFDT  (:,:)

  real(DP), allocatable, save :: t_PFTblPFSum   (:)
  real(DP), allocatable, save :: tb_PFTblPFRatio(:,:)

  real(DP), allocatable, save :: tb_PFTblPFRatioCum            (:,:)
  real(DP), allocatable, save :: b_BandIntFluxtoTotFluxRatioCum(:)

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

  logical, save :: FlagInited = .false.

contains

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

  subroutine OptPropKDInit( &
    & OptPropNcFN &
    & )

    use ni3_module


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


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

    character(128) :: PressAxisName
    character(128) :: MolAxisName
    character(128) :: TempAxisName
    character(128) :: VMRMol1AxisName
    character(128) :: VMRMol2AxisName
    character(128) :: VMRMol3AxisName
    character(128) :: BinAxisName
    character(128) :: AbsCoefVarName
    character(128) :: AtmPFRatioVarName
    character(128) :: StrPFRatioVarName

    real(DP), allocatable :: cp_StrPFRatio(:,:)
!!$    real(DP), allocatable :: bp_BandAveStrPF(:,:)
    real(DP), allocatable :: bp_BandIntFluxtoTotFluxRatio(:,:)

    integer :: VMRMol1Num
    integer :: VMRMol2Num
    integer :: VMRMol3Num

    integer :: iBin
    integer :: iBand

    integer :: iTemp

    integer :: m


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


    KDTblNBin     = -1
    KDTblNBand    = -1

    KDTblNPress   = -1
    KDTblNMol     = -1
    KDTblNTemp    = -1
    KDTblNVMRMol1 = -1
    KDTblNVMRMol2 = -1
    KDTblNVMRMol3 = -1


    OptPropNcFNSave = OptPropNcFN

    PressAxisName     = "Press"
    MolAxisName       = "MolNum"
    TempAxisName      = "Temp"
    VMRMol1AxisName   = "VMRMol1"
    VMRMol2AxisName   = "VMRMol2"
    VMRMol3AxisName   = "VMRMol3"
    BinAxisName       = 'BinNum'
    AbsCoefVarName    = 'BinAbsCoef'
    AtmPFRatioVarName = 'BinAtmPFRatio'
    StrPFRatioVarName = 'BinStrPFRatio'


    Mode = "read"
    call ni3_open( OptPropNcFNSave, Mode, InNcID )

    call ni3_inq_dimlen( InNcID, BinAxisName     , KDTblNBin    )
    call ni3_inq_dimlen( InNcID, 'BandNum'       , KDTblNBand   )

    call ni3_inq_dimlen( InNcID, PressAxisName   , KDTblNPress  )
    call ni3_inq_dimlen( InNcID, MolAxisName     , KDTblNMol    )
    call ni3_inq_dimlen( InNcID, TempAxisName    , KDTblNTemp   )
    call ni3_inq_dimlen( InNcID, VMRMol1AxisName , KDTblNVMRMol1 )
    call ni3_inq_dimlen( InNcID, VMRMol2AxisName , KDTblNVMRMol2 )
    call ni3_inq_dimlen( InNcID, VMRMol3AxisName , KDTblNVMRMol3 )


    allocate( ab_KDTblBandWaveNumBnds     (NBnd,KDTblNBand) )
    allocate( ab_KDTblBandBinIndexBnds    (NBnd,KDTblNBand) )

    call ni3_get_var( InNcID, 'BandWaveNumBnds' , ab_KDTblBandWaveNumBnds  )
    call ni3_get_var( InNcID, 'BandBinIndexBnds', ab_KDTblBandBinIndexBnds )

    allocate( c_DelWaveNum(KDTblNBin) )
    call ni3_get_var( InNcID, "BinWeight", c_DelWaveNum )
    !
    call ni3_get_var( InNcID, "VMRMol1Num", VMRMol1Num )
    call ni3_get_var( InNcID, "VMRMol2Num", VMRMol2Num )
    call ni3_get_var( InNcID, "VMRMol3Num", VMRMol3Num )
    !
    allocate( c_KDTblBinNum (KDTblNBin    ) )
    allocate( p_KDTblPress  (KDTblNPress  ) )
    allocate( m_KDTblMolNum (KDTblNMol    ) )
    allocate( t_KDTblTemp   (KDTblNTemp   ) )
    allocate( u_KDTblVMRMol1(KDTblNVMRMol1) )
    allocate( v_KDTblVMRMol2(KDTblNVMRMol2) )
    allocate( w_KDTblVMRMol3(KDTblNVMRMol3) )
    !
    call ni3_get_var( InNcID, BinAxisName    , c_KDTblBinNum  )
    call ni3_get_var( InNcID, PressAxisName  , p_KDTblPress   )
    call ni3_get_var( InNcID, MolAxisName    , m_KDTblMolNum  )
    call ni3_get_var( InNcID, TempAxisName   , t_KDTblTemp    )
    call ni3_get_var( InNcID, VMRMol1AxisName, u_KDTblVMRMol1 )
    call ni3_get_var( InNcID, VMRMol2AxisName, v_KDTblVMRMol2 )
    call ni3_get_var( InNcID, VMRMol3AxisName, w_KDTblVMRMol3 )
    !
!!$    allocate( pmctuvw_AbsCoef(KDTblNPress,KDTblNMol,KDTblNBin,KDTblNTemp,KDTblNVMRMol1,KDTblNVMRMol2,KDTblNVMRMol3) )
!!$    call ni3_get_var( InNcID, AbsCoefVarName, pmctuvw_AbsCoef )
    allocate( pmctuvw_LogAbsCoef(KDTblNPress,KDTblNMol,KDTblNBin,KDTblNTemp,KDTblNVMRMol1,KDTblNVMRMol2,KDTblNVMRMol3) )
    call ni3_get_var( InNcID, AbsCoefVarName, pmctuvw_LogAbsCoef )
    pmctuvw_LogAbsCoef = log( max( pmctuvw_LogAbsCoef, 1.0d-100 ) )
    !
    allocate( mc_IDAbsCoefType(KDTblNMol,KDTblNBin) )
    call ni3_get_var( InNcID, "BinIDAbsCoefType", mc_IDAbsCoefType )
    !
!!$    allocate( ct_AtmPFRatio(KDTblNBin,KDTblNTemp) )
!!$    call ni3_get_var( InNcID, AtmPFRatioVarName, ct_AtmPFRatio )
    allocate( pctuvw_AtmPFRatio(KDTblNPress,KDTblNBin,KDTblNTemp,KDTblNVMRMol1,KDTblNVMRMol2,KDTblNVMRMol3) )
    call ni3_get_var( InNcID, AtmPFRatioVarName, pctuvw_AtmPFRatio )
    !
    allocate( cp_StrPFRatio (KDTblNBin,KDTblNPress) )
    call ni3_get_var( InNcID, StrPFRatioVarName, cp_StrPFRatio )
    allocate( c_StrPFRatio (KDTblNBin) )
    do iBin = 1, KDTblNBin
      c_StrPFRatio(iBin) = cp_StrPFRatio(iBin,1)
    end do
    deallocate( cp_StrPFRatio )
    !
!!$    allocate( bp_BandAveStrPF(KDTblNBand,KDTblNPress) )
!!$    call ni3_get_var( InNcID, StrPFRatioVarName, bp_BandAveStrPF )
!!$    allocate( b_BandAveStrPF(KDTblNBand) )
!!$    do iBand = 1, KDTblNBand
!!$      b_BandAveStrPF(iBand) = bp_BandAveStrPF(iBand,1)
!!$    end do
!!$    deallocate( bp_BandAveStrPF )
    !
    allocate( b_BandAveRayScatCoefNonRadAct(KDTblNBand) )
    allocate( bm_BandAveRayScatCoef        (KDTblNBand,KDTblNMol) )
    call ni3_get_var( InNcID, "BandAveRayScatCoefNonRadAct", b_BandAveRayScatCoefNonRadAct )
    call ni3_get_var( InNcID, "BandAveRayScatCoef"         , bm_BandAveRayScatCoef )


!!$    call ni3_inq_dimlen( InNcID, "CharNum"          , NChar              )
    call ni3_inq_dimlen( InNcID, "PtclNum"      , KDTblNPtcl          )
    call ni3_inq_dimlen( InNcID, "PtclRadiusNum", KDTblPtclNRadiusMax )
!!$    allocate( a_CharNum          (NChar             ) )
    allocate( a_KDTblPtclNum      (KDTblNPtcl         ) )
    allocate( r_KDTblPtclRadiusNum(KDTblPtclNRadiusMax) )
    allocate( a_KDTblPtclName     (KDTblNPtcl         ) )
    allocate( a_KDTblPtclNRadius  (KDTblNPtcl         ) )
    allocate( ra_KDTblPtclRadius  (KDTblPtclNRadiusMax,KDTblNPtcl) )
!!$    call ni3_get_var( InNcID, "CharNum"          , a_CharNum           )
    call ni3_get_var( InNcID, "PtclNum"      , a_KDTblPtclNum       )
    call ni3_get_var( InNcID, "PtclRadiusNum", r_KDTblPtclRadiusNum )
    call ni3_get_var( InNcID, "PtclName"     , a_KDTblPtclName      )
    call ni3_get_var( InNcID, "PtclNRadius"  , a_KDTblPtclNRadius   )
    call ni3_get_var( InNcID, "PtclRadius"   , ra_KDTblPtclRadius   )
    allocate( bra_KDTblBandAveQExt  (KDTblNBand,KDTblPtclNRadiusMax,KDTblNPtcl) )
    allocate( bra_KDTblBandAveSSA   (KDTblNBand,KDTblPtclNRadiusMax,KDTblNPtcl) )
    allocate( bra_KDTblBandAveAF    (KDTblNBand,KDTblPtclNRadiusMax,KDTblNPtcl) )
    call ni3_get_var( InNcID, "BandAveQExt"      , bra_KDTblBandAveQExt     )
    call ni3_get_var( InNcID, "BandAveSSA"       , bra_KDTblBandAveSSA      )
    call ni3_get_var( InNcID, "BandAveAF"        , bra_KDTblBandAveAF       )


    allocate( bp_BandIntFluxtoTotFluxRatio(KDTblNBand,KDTblNPress) )

!!$    call ni3_get_var( InNcID, "RatioBandIntFluxtoTotFlux", pb_RatioBandIntFluxtoTotFlux )
!!$    allocate( b_RatioBandIntFluxtoTotFlux(KDTblNBand) )
!!$    do iBand = 1, KDTblNBand
!!$      b_RatioBandIntFluxtoTotFlux(iBand) = bp_RatioBandIntFluxtoTotFlux(iBand,1)
!!$    end do
!!$    deallocate( bp_RatioBandIntFluxtoTotFlux )

!!$    call ni3_close( InNcID )


!!$    Mode = "read"
!!$    call ni3_open( StrSpeNcFN, Mode, InNcID )
!!$    call ni3_inq_dimlen( InNcID, "KDTblNBand", TmpNum )
!!$    if ( TmpNum /= KDTblNBand ) then
!!$      write( 6, * ) 'Array size for stellar spectrum is inappropriate: ', &
!!$        TmpNum, KDTblNBand
!!$      stop
!!$    end if
    call ni3_get_var( InNcID, "RatioBandIntFluxtoTotFlux", bp_BandIntFluxtoTotFluxRatio )
    allocate( b_BandIntFluxtoTotFluxRatio   (KDTblNBand) )
    allocate( b_BandIntFluxtoTotFluxRatioCum(KDTblNBand) )
    do iBand = 1, KDTblNBand
      b_BandIntFluxtoTotFluxRatio(iBand) = bp_BandIntFluxtoTotFluxRatio(iBand,1)
    end do
    b_BandIntFluxtoTotFluxRatioCum(1) = b_BandIntFluxtoTotFluxRatio(1)
    do iBand = 1+1, KDTblNBand
      b_BandIntFluxtoTotFluxRatioCum(iBand) = &
        &   b_BandIntFluxtoTotFluxRatioCum(iBand-1) &
        & + b_BandIntFluxtoTotFluxRatio   (iBand  )
    end do
    deallocate( bp_BandIntFluxtoTotFluxRatio )


    write( 6, * ) 'Cumulative ratio of band integrated steller flux to total'
    do iBand = 1, KDTblNBand
      write( 6, * ) &
        & iBand, &
        & int( ab_KDTblBandWaveNumBnds(1,iBand) * 1.0d-2 ), '-', &
        & int( ab_KDTblBandWaveNumBnds(2,iBand) * 1.0d-2 ), ':', &
!!$        & b_BandIntFluxtoTotFluxRatio(iBand)
        & b_BandIntFluxtoTotFluxRatioCum(iBand)
    end do

    ! Read Planck function table
    !
    call ni3_inq_dimlen( InNcID, "PFTblTemp", PFTblNTemp )
    PFTblNBand = KDTblNBand
    allocate( t_PFTblTemp     (PFTblNTemp) )
!!$    allocate( b_PFTblBandNum  (PFTblNBand) )
    allocate( tb_PFTblPF      (PFTblNTemp, KDTblNBand) )
    allocate( tb_PFTblDPFDT   (PFTblNTemp, KDTblNBand) )
    call ni3_get_var( InNcID, "PFTblTemp"   , t_PFTblTemp    )
!!$    call ni3_get_var( InNcID, "PFTblBandNum", b_PFTblBandNum )
    call ni3_get_var( InNcID, "PFTblPF"     , tb_PFTblPF     )
    call ni3_get_var( InNcID, "PFTblDPFDT"  , tb_PFTblDPFDT  )


    call ni3_close( InNcID )


    ! Table for checking ratio of integrated Planck function to total
    allocate( t_PFTblPFSum      (PFTblNTemp) )
    allocate( tb_PFTblPFRatio   (PFTblNTemp, KDTblNBand) )
    allocate( tb_PFTblPFRatioCum(PFTblNTemp, KDTblNBand) )
    t_PFTblPFSum = 0.0d0
    do iBand = 1, KDTblNBand
      do iTemp = 1, PFTblNTemp
        t_PFTblPFSum(iTemp) = t_PFTblPFSum(iTemp) + tb_PFTblPF(iTemp,iBand)
      end do
    end do
    do iBand = 1, KDTblNBand
      tb_PFTblPFRatio(:,iBand) = tb_PFTblPF(:,iBand) / t_PFTblPFSum
    end do
    tb_PFTblPFRatioCum(:,KDTblNBand) = tb_PFTblPFRatio(:,KDTblNBand)
    do iBand = KDTblNBand-1, 1, -1
      tb_PFTblPFRatioCum(:,iBand) = &
        &   tb_PFTblPFRatioCum(:,iBand+1) &
        & + tb_PFTblPFRatio   (:,iBand  )
    end do
    !   check
    write( 6, * ) 'Cumulative ratio of band integrated Planck function to total'
    do iBand = 1, KDTblNBand
      write( 6, * ) &
        & iBand, &
        & int( ab_KDTblBandWaveNumBnds(1,iBand) * 1.0d-2 ), '-', &
        & int( ab_KDTblBandWaveNumBnds(2,iBand) * 1.0d-2 ), ':', &
!!$        & tb_PFTblPFRatio(PFTblNTemp,iBand)
        & tb_PFTblPFRatioCum(PFTblNTemp,iBand)
    end do



    do m = 1, KDTblNMol
      if ( m_KDTblMolNum(m) == VMRMol1Num ) exit
    end do
    if ( m > KDTblNMol ) then
      stop 'Unexpected events when searching index for VMRMol1Num.'
    end if
    KDTblIndexVMRMol1Num = m
    !
    do m = 1, KDTblNMol
      if ( m_KDTblMolNum(m) == VMRMol2Num ) exit
    end do
    if ( m > KDTblNMol ) then
      stop 'Unexpected events when searching index for VMRMol2Num.'
    end if
    KDTblIndexVMRMol2Num = m
    !
    do m = 1, KDTblNMol
      if ( m_KDTblMolNum(m) == VMRMol3Num ) exit
    end do
    if ( m > KDTblNMol ) then
      stop 'Unexpected events when searching index for VMRMol3Num.'
    end if
    KDTblIndexVMRMol3Num = m


    ! check key and non-key molecules
    if ( mc_IDAbsCoefType(KDTblIndexVMRMol1Num,1) == 0 ) then
      if ( KDTblNVMRMol1 > 1 ) then
        write( 6, * ) &
          & 'VMR Axis 1 has an unexpected size nevertheless it is non-key molecule: ', &
          & KDTblNVMRMol1
        stop
      end if
    end if
    if ( mc_IDAbsCoefType(KDTblIndexVMRMol2Num,2) == 0 ) then
      if ( KDTblNVMRMol2 > 1 ) then
        write( 6, * ) &
          & 'VMR Axis 2 has an unexpected size nevertheless it is non-key molecule: ', &
          & KDTblNVMRMol2
        stop
      end if
    end if
    if ( mc_IDAbsCoefType(KDTblIndexVMRMol3Num,1) == 0 ) then
      if ( KDTblNVMRMol3 > 1 ) then
        write( 6, * ) &
          & 'VMR Axis 3 has an unexpected size nevertheless it is non-key molecule: ', &
          & KDTblNVMRMol3
        stop
      end if
    end if


    ! Preparation of tables for Planck function
    !
!!$    call OptPropKDPFTblInit


    FlagInited = .true.


  end subroutine OptPropKDInit

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

  subroutine OptPropKDFinalize


    ! local variables
    !


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


    KDTblNBin     = -1
    KDTblNBand    = -1

    KDTblNPress   = -1
    KDTblNMol     = -1
    KDTblNTemp    = -1
    KDTblNVMRMol1 = -1
    KDTblNVMRMol2 = -1
    KDTblNVMRMol3 = -1


    deallocate( ab_KDTblBandWaveNumBnds  )
    deallocate( ab_KDTblBandBinIndexBnds )
    !
    deallocate( c_DelWaveNum )
    !
    deallocate( c_KDTblBinNum  )
    deallocate( p_KDTblPress   )
    deallocate( m_KDTblMolNum  )
    deallocate( t_KDTblTemp    )
    deallocate( u_KDTblVMRMol1 )
    deallocate( v_KDTblVMRMol2 )
    deallocate( w_KDTblVMRMol3 )
    !
!!$    deallocate( pmctuvw_AbsCoef )
    deallocate( pmctuvw_LogAbsCoef )
!!$    deallocate( ct_AtmPFRatio )
    deallocate( pctuvw_AtmPFRatio )
    deallocate( c_StrPFRatio )
!!$    deallocate( b_BandAveStrPF )
    deallocate( b_BandIntFluxtoTotFluxRatio )

    deallocate( b_BandAveRayScatCoefNonRadAct )
    deallocate( bm_BandAveRayScatCoef         )

!!$    deallocate( a_CharNum           )
    deallocate( a_KDTblPtclNum       )
    deallocate( r_KDTblPtclRadiusNum )
    deallocate( a_KDTblPtclName      )
    deallocate( a_KDTblPtclNRadius   )
    deallocate( ra_KDTblPtclRadius   )


    ! Finalize tables for Planck function
    !
    call OptPropKDPFTblFinalize


    FlagInited = .false.


  end subroutine OptPropKDFinalize

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

  function OptPropKDGetNWN() result( NWN )


    integer :: NWN

    ! local variables
    !


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


    NWN = KDTblNBin


  end function OptPropKDGetNWN

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

  function OptPropKDGetNPress() result( NPress )


    integer :: NPress

    ! local variables
    !


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


    NPress = KDTblNPress


  end function OptPropKDGetNPress

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

  function OptPropKDGetNMol() result( NMol )


    integer :: NMol

    ! local variables
    !


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


    NMol = KDTblNMol


  end function OptPropKDGetNMol

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

  function OptPropKDGetNBand() result( NBand )


    integer :: NBand

    ! local variables
    !


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


    NBand = KDTblNBand


  end function OptPropKDGetNBand

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

  function OptPropKDGetDelWN( 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 = c_DelWaveNum(iWaveNum)


  end function OptPropKDGetDelWN

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

  function OptPropKDGetBandNum( iWaveNum ) result( BandNum )


    integer, intent(in) :: iWaveNum

    integer :: BandNum


    ! local variables
    !
    integer :: n


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


    do n = 1, KDTblNBand
      if ( ( ab_KDTblBandBinIndexBnds(1,n) <= iWaveNum ) .and. &
        &  ( iWaveNum <= ab_KDTblBandBinIndexBnds(2,n) ) ) &
        exit
    end do

    BandNum = n


  end function OptPropKDGetBandNum

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

  subroutine OptPropKDGetBandWNBnds( &
    & 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 /= KDTblNBand ) then
      write( 6, * ) 'Number of bands is inappropriate: ', NBand, KDTblNBand
      stop
    end if

    aa_BandWNBnds = ab_KDTblBandWaveNumBnds


  end subroutine OptPropKDGetBandWNBnds

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

  subroutine OptPropKDGetBandBinIndexBnds( &
    & NBand, &
    & ab_BandBinIndexBnds &
    & )

    integer, intent(in ) :: NBand
    integer, intent(out) :: ab_BandBinIndexBnds(2,NBand)


    ! local variables
    !


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


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

    ab_BandBinIndexBnds = ab_KDTblBandBinIndexBnds


  end subroutine OptPropKDGetBandBinIndexBnds

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

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

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


    ! local variables
    !


!!$    StrBandAvePF = b_BandAveStrPF( iBand )
    StrBandAvePF = b_BandIntFluxtoTotFluxRatio( iBand )
    StrBandAvePF = StrBandAvePF * StrFluxTOA
    StrBandAvePF = StrBandAvePF &
      & / ( ab_KDTblBandWaveNumBnds(2,iBand) - ab_KDTblBandWaveNumBnds(1,iBand) )


  end subroutine OptPropKDGetStrBandAvePF

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

!!$  subroutine OptPropKDGetTotalStrFlux( &
!!$    & TotalStrFlux &
!!$    & )
!!$
!!$    real(DP), intent(out) :: TotalStrFlux
!!$
!!$
!!$    ! local variables
!!$    !
!!$    real(DP) :: DelWaveNum
!!$
!!$    integer  :: l
!!$
!!$
!!$    TotalStrFlux = 0.0d0
!!$    do l = 1, KDTblNBand
!!$      DelWaveNum = ab_KDTblBandWaveNumBnds(2,l) - ab_KDTblBandWaveNumBnds(1,l)
!!$      TotalStrFlux = TotalStrFlux &
!!$        & + b_BandAveStrPF( l ) * DelWaveNum
!!$    end do
!!$
!!$
!!$  end subroutine OptPropKDGetTotalStrFlux

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

  subroutine OptPropKDCalcIndexAndWeight( &
    & kmax, NMol, 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 &
    & )


    integer     , intent(in ) :: kmax
    integer     , intent(in ) :: NMol
    integer     , intent(in ) :: m_MolNum  (1:NMol)
    real(DP)    , intent(in ) :: r_Press   (1:kmax)
    real(DP)    , intent(in ) :: r_Temp    (1:kmax)
    real(DP)    , intent(in ) :: rm_VMR    (1:kmax,1:NMol)
    real(DP)    , intent(in ) :: SurfTemp
    integer     , intent(in ) :: NIntpol
    integer     , intent(out) :: r_KDTblPressIndex  (1:kmax)
    integer     , intent(out) :: r_KDTblTempIndex   (1:kmax)
    integer     , intent(out) :: r_KDTblVMRMol1Index(1:kmax)
    integer     , intent(out) :: r_KDTblVMRMol2Index(1:kmax)
    integer     , intent(out) :: r_KDTblVMRMol3Index(1:kmax)
    integer     , intent(out) :: r_NIntpol1         (1:kmax)
    integer     , intent(out) :: r_NIntpol3         (1:kmax)
    integer     , intent(out) :: r_NIntpol4         (1:kmax)
    integer     , intent(out) :: r_NIntpol5         (1:kmax)
    real(DP)    , intent(out) :: ptuvwr_Weight(NIntpol,NIntpol,NIntpol,NIntpol,NIntpol,1:kmax)
    integer     , intent(out) :: r_KDTblPressIndexSfc  (1)
    integer     , intent(out) :: r_KDTblTempIndexSfc   (1)
    integer     , intent(out) :: r_KDTblVMRMol1IndexSfc(1)
    integer     , intent(out) :: r_KDTblVMRMol2IndexSfc(1)
    integer     , intent(out) :: r_KDTblVMRMol3IndexSfc(1)
    integer     , intent(out) :: r_NIntpol1Sfc         (1)
    integer     , intent(out) :: r_NIntpol3Sfc         (1)
    integer     , intent(out) :: r_NIntpol4Sfc         (1)
    integer     , intent(out) :: r_NIntpol5SFc         (1)
    real(DP)    , intent(out) :: ptuvw_WeightSfc(NIntpol,NIntpol,NIntpol,NIntpol,NIntpol)


    ! local variables
    !
    !       It should be noticed that NIntpol is defined twice in this module!
!!$    integer , parameter :: NIntpol = 2
!!$    integer , parameter :: NIntpol = 3
    real(DP)            :: a_Grid   (NIntpol)
    real(DP)            :: a_Weight (NIntpol)
    real(DP)            :: p_Weight1(NIntpol)
    real(DP)            :: t_Weight2(NIntpol)
    real(DP)            :: u_Weight3(NIntpol)
    real(DP)            :: v_Weight4(NIntpol)
    real(DP)            :: w_Weight5(NIntpol)
!!$    integer             :: NIntpolMax
    integer             :: NIntpol1
    integer             :: NIntpol3
    integer             :: NIntpol4
    integer             :: NIntpol5
    integer             :: IndexVMRMol1Num
    integer             :: IndexVMRMol2Num
    integer             :: IndexVMRMol3Num
    real(DP)       :: Val
    real(DP)       :: Press
    real(DP)       :: Temp
    real(DP)       :: VMRMol1
    real(DP)       :: VMRMol2
    real(DP)       :: VMRMol3

    integer        :: iPress
    integer        :: iTemp
    integer        :: iVMRMol1
    integer        :: iVMRMol2
    integer        :: iVMRMol3
    integer        :: k
    integer        :: l
    integer        :: m
    integer        :: n
    integer        :: n1
    integer        :: n2
    integer        :: n3
    integer        :: n4
    integer        :: n5


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


    ! check size of VMR of argument and number of molecules in table
    if ( NMol < KDTblNMol ) then
      write( 6, * ) 'Array size for VMR is not enough: ', NMol, KDTblNMol
      stop
    end if
    do m = 1, NMol
      if ( m_MolNum(m) /= m_KDTblMolNum(m) ) then
        write( 6, * ) 'Given molecular number, ', m_MolNum(m), ',  does not match with that in table.'
        write( 6, * ) '  Argument: ', m_MolNum
        write( 6, * ) '  Table   : ', m_KDTblMolNum
        stop
      end if
    end do


    if ( KDTblNPress < NIntpol ) then
      write( 6, * ) 'Array size for Press is not enough for interpolation: ', &
        & KDTblNPress, NIntpol
      stop
    end if
    if ( KDTblNTemp < NIntpol ) then
      write( 6, * ) 'Array size for Temp is not enough for interpolation: ', &
        & KDTblNTemp, NIntpol
      stop
    end if
    if ( KDTblNVMRMol1 > 1 ) then
      if ( KDTblNVMRMol1 < NIntpol ) then
        write( 6, * ) 'Array size for VMRMol1 is not enough for interpolation: ', &
          & KDTblNVMRMol1, NIntpol
        stop
      end if
    end if
    if ( KDTblNVMRMol2 > 1 ) then
      if ( KDTblNVMRMol2 < NIntpol ) then
        write( 6, * ) 'Array size for VMRMol2 is not enough for interpolation: ', &
          & KDTblNVMRMol2, NIntpol
        stop
      end if
    end if
    if ( KDTblNVMRMol3 > 1 ) then
      if ( KDTblNVMRMol3 < NIntpol ) then
        write( 6, * ) 'Array size for VMRMol3 is not enough for interpolation: ', &
          & KDTblNVMRMol3, NIntpol
        stop
      end if
    end if


    do m = 1, NMol
      if ( m_MolNum(m) == m_KDTblMolNum(KDTblIndexVMRMol1Num) ) exit
    end do
    if ( m > NMol ) then
      write( 6, * ) &
        & 'Unable to find a molecular number, ', &
        & m_KDTblMolNum(KDTblIndexVMRMol1Num), ', in MolNum in argument.'
      stop
    end if
    IndexVMRMol1Num = m
    !
    do m = 1, NMol
      if ( m_MolNum(m) == m_KDTblMolNum(KDTblIndexVMRMol2Num) ) exit
    end do
    if ( m > NMol ) then
      write( 6, * ) &
        & 'Unable to find a molecular number, ', &
        & m_KDTblMolNum(KDTblIndexVMRMol2Num), ', in MolNum in argument.'
      stop
    end if
    IndexVMRMol2Num = m
    !
    do m = 1, NMol
      if ( m_MolNum(m) == m_KDTblMolNum(KDTblIndexVMRMol3Num) ) exit
    end do
    if ( m > NMol ) then
      write( 6, * ) &
        & 'Unable to find a molecular number, ', &
        & m_KDTblMolNum(KDTblIndexVMRMol3Num), ', in MolNum in argument.'
      stop
    end if
    IndexVMRMol3Num = m


    call FindTblIndex( &
      & 'Press', &
      & KDTblNPress, p_KDTblPress, NIntpol, kmax, r_Press, &
      & r_KDTblPressIndex &
      & )
    call FindTblIndex( &
      & 'Temp', &
      & KDTblNTemp, t_KDTblTemp, NIntpol, kmax, r_Temp, &
      & r_KDTblTempIndex &
      & )
    if ( KDTblNVMRMol1 > 1 ) then
      call FindTblIndex( &
        & 'VMRMol1', &
        & KDTblNVMRMol1, u_KDTblVMRMol1, NIntpol, kmax, rm_VMR(:,IndexVMRMol1Num), &
        & r_KDTblVMRMol1Index &
        & )
    else
      r_KDTblVMRMol1Index = 1
    end if
    if ( KDTblNVMRMol2 > 1 ) then
      call FindTblIndex( &
        & 'VMRMol2', &
        & KDTblNVMRMol2, v_KDTblVMRMol2, NIntpol, kmax, rm_VMR(:,IndexVMRMol2Num), &
        & r_KDTblVMRMol2Index &
        & )
    else
      r_KDTblVMRMol2Index = 1
    end if
    if ( KDTblNVMRMol3 > 1 ) then
      call FindTblIndex( &
        & 'VMRMol3', &
        & KDTblNVMRMol3, w_KDTblVMRMol3, NIntpol, kmax, rm_VMR(:,IndexVMRMol3Num), &
        & r_KDTblVMRMol3Index &
        & )
    else
      r_KDTblVMRMol3Index = 1
    end if


    do k = 1, kmax
      Press   = r_Press (k)
      Temp    = r_Temp  (k)
      VMRMol1 = rm_VMR  (k,IndexVMRMol1Num)
      VMRMol2 = rm_VMR  (k,IndexVMRMol2Num)
      VMRMol3 = rm_VMR  (k,IndexVMRMol3Num)

      iPress   = r_KDTblPressIndex  (k)
      iTemp    = r_KDTblTempIndex   (k)
      iVMRMol1 = r_KDTblVMRMol1Index(k)
      iVMRMol2 = r_KDTblVMRMol2Index(k)
      iVMRMol3 = r_KDTblVMRMol3Index(k)


      call CalcWeight( &
        & KDTblNPress, p_KDTblPress, iPress, &
        & Press, &
        & NIntpol, &
        & p_Weight1, NIntpol1 &
        & )
      !
      do n = 1, NIntpol
        a_Grid(n) = t_KDTblTemp(iTemp+n-1)
      end do
      Val = Temp
      a_Weight = 1.0d0
      do n = 1, NIntpol
        do l = 1, NIntpol
          if ( l == n ) cycle
          a_Weight(n) = a_Weight(n) &
            & * ( Val - a_Grid(l) ) / ( a_Grid(n) - a_Grid(l) )
        end do
      end do
      t_Weight2 = a_Weight
      !
      call CalcWeight( &
        & KDTblNVMRMol1, u_KDTblVMRMol1, iVMRMol1, &
        & VMRMol1, &
        & NIntpol, &
        & u_Weight3, NIntpol3 &
        & )
      !
      call CalcWeight( &
        & KDTblNVMRMol2, v_KDTblVMRMol2, iVMRMol2, &
        & VMRMol2, &
        & NIntpol, &
        & v_Weight4, NIntpol4 &
        & )
      !
      call CalcWeight( &
        & KDTblNVMRMol3, w_KDTblVMRMol3, iVMRMol3, &
        & VMRMol3, &
        & NIntpol, &
        & w_Weight5, NIntpol5 &
        & )
      !
      r_NIntpol1(k) = NIntpol1
      r_NIntpol3(k) = NIntpol3
      r_NIntpol4(k) = NIntpol4
      r_NIntpol5(k) = NIntpol5
      !
      ptuvwr_Weight(:,:,:,:,:,k) = 1.0d0
      do n5 = 1, NIntpol5
        do n4 = 1, NIntpol4
          do n3 = 1, NIntpol3
            do n2 = 1, NIntpol
              do n1 = 1, NIntpol1
                ptuvwr_Weight(n1,n2,n3,n4,n5,k) = ptuvwr_Weight(n1,n2,n3,n4,n5,k) &
                  & * p_Weight1(n1) * t_Weight2(n2) &
                  & * u_Weight3(n3) * v_Weight4(n4) * w_Weight5(n5)
              end do
            end do
          end do
        end do
      end do

    end do


    r_KDTblPressIndexSfc  (1) = r_KDTblPressIndex  (1)
    call FindTblIndex( &
      & 'Temp', &
      & KDTblNTemp, t_KDTblTemp, NIntpol, 1, (/SurfTemp/), &
      & r_KDTblTempIndexSfc(1:1) &
      & )
    r_KDTblVMRMol1IndexSfc(1) = r_KDTblVMRMol1Index(1)
    r_KDTblVMRMol2IndexSfc(1) = r_KDTblVMRMol2Index(1)
    r_KDTblVMRMol3IndexSfc(1) = r_KDTblVMRMol3Index(1)
    do k = 1, 1
      Press   = r_Press (k)
      Temp    = SurfTemp
      VMRMol1 = rm_VMR  (k,IndexVMRMol1Num)
      VMRMol2 = rm_VMR  (k,IndexVMRMol2Num)
      VMRMol3 = rm_VMR  (k,IndexVMRMol3Num)

      iPress   = r_KDTblPressIndexSfc  (k)
      iTemp    = r_KDTblTempIndexSfc   (k)
      iVMRMol1 = r_KDTblVMRMol1IndexSfc(k)
      iVMRMol2 = r_KDTblVMRMol2IndexSfc(k)
      iVMRMol3 = r_KDTblVMRMol3IndexSfc(k)


      call CalcWeight( &
        & KDTblNPress, p_KDTblPress, iPress, &
        & Press, &
        & NIntpol, &
        & p_Weight1, NIntpol1 &
        & )
      !
      do n = 1, NIntpol
        a_Grid(n) = t_KDTblTemp(iTemp+n-1)
      end do
      Val = Temp
      a_Weight = 1.0d0
      do n = 1, NIntpol
        do l = 1, NIntpol
          if ( l == n ) cycle
          a_Weight(n) = a_Weight(n) &
            & * ( Val - a_Grid(l) ) / ( a_Grid(n) - a_Grid(l) )
        end do
      end do
      t_Weight2 = a_Weight
      !
      call CalcWeight( &
        & KDTblNVMRMol1, u_KDTblVMRMol1, iVMRMol1, &
        & VMRMol1, &
        & NIntpol, &
        & u_Weight3, NIntpol3 &
        & )
      !
      call CalcWeight( &
        & KDTblNVMRMol2, v_KDTblVMRMol2, iVMRMol2, &
        & VMRMol2, &
        & NIntpol, &
        & v_Weight4, NIntpol4 &
        & )
      !
      call CalcWeight( &
        & KDTblNVMRMol3, w_KDTblVMRMol3, iVMRMol3, &
        & VMRMol3, &
        & NIntpol, &
        & w_Weight5, NIntpol5 &
        & )
      !
      r_NIntpol1Sfc(k) = NIntpol1
      r_NIntpol3Sfc(k) = NIntpol3
      r_NIntpol4Sfc(k) = NIntpol4
      r_NIntpol5Sfc(k) = NIntpol5
      !
      ptuvw_WeightSfc = 1.0d0
      do n5 = 1, NIntpol5
        do n4 = 1, NIntpol4
          do n3 = 1, NIntpol3
            do n2 = 1, NIntpol
!!$              do n1 = 1, NIntpol
              do n1 = 1, NIntpol1
                ptuvw_WeightSfc(n1,n2,n3,n4,n5) = ptuvw_WeightSfc(n1,n2,n3,n4,n5) &
                  & * p_Weight1(n1) * t_Weight2(n2) &
                  & * u_Weight3(n3) * v_Weight4(n4) * w_Weight5(n5)
              end do
            end do
          end do
        end do
      end do

    end do


  end subroutine OptPropKDCalcIndexAndWeight

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

  subroutine OptPropKDGetAbsCoefProf( &
    & iWaveNum, kmax, NMol, m_MolNum, r_Press, r_Temp, rm_VMR, &
    & rm_AbsCoef &
    & )


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


    ! local variables
    !
    integer        :: r_KDTblPressIndex  (1:kmax)
    integer        :: r_KDTblTempIndex   (1:kmax)
    integer        :: r_KDTblVMRMol1Index(1:kmax)
    integer        :: r_KDTblVMRMol2Index(1:kmax)
    integer        :: r_KDTblVMRMol3Index(1:kmax)
    !       It should be noticed that NIntpol is defined twice in this module!
    integer , parameter :: NIntpol = 2
!!$    integer , parameter :: NIntpol = 3
    real(DP)            :: a_Grid   (NIntpol)
    real(DP)            :: a_Weight (NIntpol)
    real(DP)            :: p_Weight1(NIntpol)
    real(DP)            :: t_Weight2(NIntpol)
    real(DP)            :: u_Weight3(NIntpol)
    real(DP)            :: v_Weight4(NIntpol)
    real(DP)            :: w_Weight5(NIntpol)
    real(DP)            :: ptuvw_Weight(NIntpol,NIntpol,NIntpol,NIntpol,NIntpol)
    real(DP)            :: ptuvw_Val   (NIntpol,NIntpol,NIntpol,NIntpol,NIntpol)
    real(DP)            :: m_Val       (KDTblNMol)
!!$    integer             :: NIntpolMax
    integer             :: NIntpol1
    integer             :: NIntpol3
    integer             :: NIntpol4
    integer             :: NIntpol5
    integer             :: IndexVMRMol1Num
    integer             :: IndexVMRMol2Num
    integer             :: IndexVMRMol3Num
    real(DP)       :: Val
    real(DP)       :: Press
    real(DP)       :: Temp
    real(DP)       :: VMRMol1
    real(DP)       :: VMRMol2
    real(DP)       :: VMRMol3

    integer        :: LoopLimit1
    integer        :: LoopLimit2

    integer        :: iPress
    integer        :: iTemp
    integer        :: iVMRMol1
    integer        :: iVMRMol2
    integer        :: iVMRMol3
    integer        :: k
    integer        :: l
    integer        :: m
    integer        :: n
    integer        :: l1
    integer        :: l2
    integer        :: l3
    integer        :: l4
    integer        :: l5
    integer        :: n1
    integer        :: n2
    integer        :: n3
    integer        :: n4
    integer        :: n5


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


    ! check size of VMR of argument and number of molecules in table
    if ( NMol < KDTblNMol ) then
      write( 6, * ) 'Array size for VMR is not enough: ', NMol, KDTblNMol
      stop
    end if
    do m = 1, NMol
      if ( m_MolNum(m) /= m_KDTblMolNum(m) ) then
        write( 6, * ) 'Given molecular number, ', m_MolNum(m), ',  does not match with that in table.'
        write( 6, * ) '  Argument: ', m_MolNum
        write( 6, * ) '  Table   : ', m_KDTblMolNum
        stop
      end if
    end do


    if ( KDTblNPress < NIntpol ) then
      write( 6, * ) 'Array size for Press is not enough for interpolation: ', &
        & KDTblNPress, NIntpol
      stop
    end if
    if ( KDTblNTemp < NIntpol ) then
      write( 6, * ) 'Array size for Temp is not enough for interpolation: ', &
        & KDTblNTemp, NIntpol
      stop
    end if
    if ( KDTblNVMRMol1 > 1 ) then
      if ( KDTblNVMRMol1 < NIntpol ) then
        write( 6, * ) 'Array size for VMRMol1 is not enough for interpolation: ', &
          & KDTblNVMRMol1, NIntpol
        stop
      end if
    end if
    if ( KDTblNVMRMol2 > 1 ) then
      if ( KDTblNVMRMol2 < NIntpol ) then
        write( 6, * ) 'Array size for VMRMol2 is not enough for interpolation: ', &
          & KDTblNVMRMol2, NIntpol
        stop
      end if
    end if
    if ( KDTblNVMRMol3 > 1 ) then
      if ( KDTblNVMRMol3 < NIntpol ) then
        write( 6, * ) 'Array size for VMRMol3 is not enough for interpolation: ', &
          & KDTblNVMRMol3, NIntpol
        stop
      end if
    end if


    do m = 1, NMol
      if ( m_MolNum(m) == m_KDTblMolNum(KDTblIndexVMRMol1Num) ) exit
    end do
    if ( m > NMol ) then
      write( 6, * ) &
        & 'Unable to find a molecular number, ', &
        & m_KDTblMolNum(KDTblIndexVMRMol1Num), ', in MolNum in argument.'
      stop
    end if
    IndexVMRMol1Num = m
    !
    do m = 1, NMol
      if ( m_MolNum(m) == m_KDTblMolNum(KDTblIndexVMRMol2Num) ) exit
    end do
    if ( m > NMol ) then
      write( 6, * ) &
        & 'Unable to find a molecular number, ', &
        & m_KDTblMolNum(KDTblIndexVMRMol2Num), ', in MolNum in argument.'
      stop
    end if
    IndexVMRMol2Num = m
    !
    do m = 1, NMol
      if ( m_MolNum(m) == m_KDTblMolNum(KDTblIndexVMRMol3Num) ) exit
    end do
    if ( m > NMol ) then
      write( 6, * ) &
        & 'Unable to find a molecular number, ', &
        & m_KDTblMolNum(KDTblIndexVMRMol3Num), ', in MolNum in argument.'
      stop
    end if
    IndexVMRMol3Num = m


    call FindTblIndex( &
      & 'Press', &
      & KDTblNPress, p_KDTblPress, NIntpol, kmax, r_Press, &
      & r_KDTblPressIndex &
      & )
    call FindTblIndex( &
      & 'Temp', &
      & KDTblNTemp, t_KDTblTemp, NIntpol, kmax, r_Temp, &
      & r_KDTblTempIndex &
      & )
    if ( KDTblNVMRMol1 > 1 ) then
      call FindTblIndex( &
        & 'VMRMol1', &
        & KDTblNVMRMol1, u_KDTblVMRMol1, NIntpol, kmax, rm_VMR(:,IndexVMRMol1Num), &
        & r_KDTblVMRMol1Index &
        & )
    else
      r_KDTblVMRMol1Index = 1
    end if
    if ( KDTblNVMRMol2 > 1 ) then
      call FindTblIndex( &
        & 'VMRMol2', &
        & KDTblNVMRMol2, v_KDTblVMRMol2, NIntpol, kmax, rm_VMR(:,IndexVMRMol2Num), &
        & r_KDTblVMRMol2Index &
        & )
    else
      r_KDTblVMRMol2Index = 1
    end if
    if ( KDTblNVMRMol3 > 1 ) then
      call FindTblIndex( &
        & 'VMRMol3', &
        & KDTblNVMRMol3, w_KDTblVMRMol3, NIntpol, kmax, rm_VMR(:,IndexVMRMol3Num), &
        & r_KDTblVMRMol3Index &
        & )
    else
      r_KDTblVMRMol3Index = 1
    end if


    LoopLimit1 = min( NMol, KDTblNMol )
    LoopLimit2 = max( NMol, KDTblNMol )

    do k = 1, kmax
      Press   = r_Press (k)
      Temp    = r_Temp  (k)
      VMRMol1 = rm_VMR  (k,IndexVMRMol1Num)
      VMRMol2 = rm_VMR  (k,IndexVMRMol2Num)
      VMRMol3 = rm_VMR  (k,IndexVMRMol3Num)

      iPress   = r_KDTblPressIndex  (k)
      iTemp    = r_KDTblTempIndex   (k)
      iVMRMol1 = r_KDTblVMRMol1Index(k)
      iVMRMol2 = r_KDTblVMRMol2Index(k)
      iVMRMol3 = r_KDTblVMRMol3Index(k)


      call CalcWeight( &
        & KDTblNPress, p_KDTblPress, iPress, &
        & Press, &
        & NIntpol, &
        & p_Weight1, NIntpol1 &
        & )
      !
      do n = 1, NIntpol
        a_Grid(n) = t_KDTblTemp(iTemp+n-1)
      end do
      Val = Temp
      a_Weight = 1.0d0
      do n = 1, NIntpol
        do l = 1, NIntpol
          if ( l == n ) cycle
          a_Weight(n) = a_Weight(n) &
            & * ( Val - a_Grid(l) ) / ( a_Grid(n) - a_Grid(l) )
        end do
      end do
      t_Weight2 = a_Weight
      !
      call CalcWeight( &
        & KDTblNVMRMol1, u_KDTblVMRMol1, iVMRMol1, &
        & VMRMol1, &
        & NIntpol, &
        & u_Weight3, NIntpol3 &
        & )
      !
      call CalcWeight( &
        & KDTblNVMRMol2, v_KDTblVMRMol2, iVMRMol2, &
        & VMRMol2, &
        & NIntpol, &
        & v_Weight4, NIntpol4 &
        & )
      !
      call CalcWeight( &
        & KDTblNVMRMol3, w_KDTblVMRMol3, iVMRMol3, &
        & VMRMol3, &
        & NIntpol, &
        & w_Weight5, NIntpol5 &
        & )
      !
      ptuvw_Weight = 1.0d0
      do n5 = 1, NIntpol5
        do n4 = 1, NIntpol4
          do n3 = 1, NIntpol3
            do n2 = 1, NIntpol
              do n1 = 1, NIntpol1
                ptuvw_Weight(n1,n2,n3,n4,n5) = ptuvw_Weight(n1,n2,n3,n4,n5) &
                  & * p_Weight1(n1) * t_Weight2(n2) &
                  & * u_Weight3(n3) * v_Weight4(n4) * w_Weight5(n5)
              end do
            end do
          end do
        end do
      end do


      do m = 1, KDTblNMol

        if ( mc_IDAbsCoefType(m,iWaveNum) == -1 ) then
          ! zero abs. coef.

          m_Val(m) = 0.0d0

        else
          ! else

          do n5 = 1, NIntpol5
            do n4 = 1, NIntpol4
              do n3 = 1, NIntpol3
                do n2 = 1, NIntpol
                  do n1 = 1, NIntpol
                    l1 = iPress  +n1-1
                    l2 = iTemp   +n2-1
                    l3 = iVMRMol1+n3-1
                    l4 = iVMRMol2+n4-1
                    l5 = iVMRMol3+n5-1
!!$                aaaam_Val(n1,n2,n3,n4,m) = pmctvv_AbsCoef(l1,m,iWaveNum,l2,l3,l4)
!!$                  ptuvwm_Val(n1,n2,n3,n4,n5,m) = &
!!$                    & log( max( pmctuvw_AbsCoef(l1,m,iWaveNum,l2,l3,l4,l5), 1.0d-100 ) )
                    ptuvw_Val(n1,n2,n3,n4,n5) = &
                      & pmctuvw_LogAbsCoef(l1,m,iWaveNum,l2,l3,l4,l5)
                  end do
                end do
              end do
            end do
          end do

          m_Val(m) = 0.0d0
          do n5 = 1, NIntpol5
            do n4 = 1, NIntpol4
              do n3 = 1, NIntpol3
                do n2 = 1, NIntpol
                  do n1 = 1, NIntpol
                    m_Val(m) = m_Val(m) &
                      & + ptuvw_Val(n1,n2,n3,n4,n5) * ptuvw_Weight(n1,n2,n3,n4,n5)
                  end do
                end do
              end do
            end do
          end do

          m_Val(m) = exp( m_Val(m) )

        end if

      end do

      do m = 1, LoopLimit1
        rm_AbsCoef(k,m) = m_Val(m)
      end do
      do m = LoopLimit1+1, LoopLimit2
        rm_AbsCoef(k,m) = 0.0d0
      end do

    end do

    do m = 1, NMol
      select case ( mc_IDAbsCoefType(m,iWaveNum) )
      case ( -1 ) ! zero abs. coef.
        rm_AbsCoef(:,m) = 0.0d0
      case (  0 ) ! abs. coef. per unit abs. molecule
        rm_AbsCoef(:,m) = rm_AbsCoef(:,m) * rm_VMR(:,m)
!!$      case (  1 ) ! abs. coef. per unit atm. molecule
!!$        rm_AbsCoef(:,m) = rm_AbsCoef(:,m)
      end select
    end do

  end subroutine OptPropKDGetAbsCoefProf

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

  subroutine OptPropKDGetAbsCoefProfWithIndex( &
    & iWaveNum, kmax, NMol, m_MolNum, r_Press, r_Temp, rm_VMR, &
    & rm_AbsCoef, &
    & NIntpol, &
    & r_KDTblPressIndex, r_KDTblTempIndex, &
    & r_KDTblVMRMol1Index, r_KDTblVMRMol2Index, r_KDTblVMRMol3Index, &
    & r_NIntpol1, r_NIntpol3, r_NIntpol4, r_NIntpol5, &
    & ptuvwr_Weight &
    & )


    integer     , intent(in ) :: iWaveNum
    integer     , intent(in ) :: kmax
    integer     , intent(in ) :: NMol
    integer     , intent(in ) :: m_MolNum  (1:NMol)
    real(DP)    , intent(in ) :: r_Press   (1:kmax)
    real(DP)    , intent(in ) :: r_Temp    (1:kmax)
    real(DP)    , intent(in ) :: rm_VMR    (1:kmax,1:NMol)
    real(DP)    , intent(out) :: rm_AbsCoef(1:kmax,1:NMol)
    integer     , intent(in ) :: NIntpol
    integer     , intent(in ) :: r_KDTblPressIndex  (1:kmax)
    integer     , intent(in ) :: r_KDTblTempIndex   (1:kmax)
    integer     , intent(in ) :: r_KDTblVMRMol1Index(1:kmax)
    integer     , intent(in ) :: r_KDTblVMRMol2Index(1:kmax)
    integer     , intent(in ) :: r_KDTblVMRMol3Index(1:kmax)
    integer     , intent(in ) :: r_NIntpol1         (1:kmax)
    integer     , intent(in ) :: r_NIntpol3         (1:kmax)
    integer     , intent(in ) :: r_NIntpol4         (1:kmax)
    integer     , intent(in ) :: r_NIntpol5         (1:kmax)
    real(DP)    , intent(in ) :: ptuvwr_Weight(NIntpol,NIntpol,NIntpol,NIntpol,NIntpol,1:kmax)


    ! local variables
    !
    real(DP)            :: ptuvw_Val   (NIntpol,NIntpol,NIntpol,NIntpol,NIntpol)
    real(DP)            :: m_Val       (KDTblNMol)

    integer             :: NIntpol1
    integer             :: NIntpol3
    integer             :: NIntpol4
    integer             :: NIntpol5

    integer        :: LoopLimit1
    integer        :: LoopLimit2

    integer        :: iPress
    integer        :: iTemp
    integer        :: iVMRMol1
    integer        :: iVMRMol2
    integer        :: iVMRMol3
    integer        :: k
    integer        :: m
    integer        :: l1
    integer        :: l2
    integer        :: l3
    integer        :: l4
    integer        :: l5
    integer        :: n1
    integer        :: n2
    integer        :: n3
    integer        :: n4
    integer        :: n5


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



    LoopLimit1 = min( NMol, KDTblNMol )
    LoopLimit2 = max( NMol, KDTblNMol )

    do k = 1, kmax
      iPress   = r_KDTblPressIndex  (k)
      iTemp    = r_KDTblTempIndex   (k)
      iVMRMol1 = r_KDTblVMRMol1Index(k)
      iVMRMol2 = r_KDTblVMRMol2Index(k)
      iVMRMol3 = r_KDTblVMRMol3Index(k)

      NIntpol1 = r_NIntpol1         (k)
      NIntpol3 = r_NIntpol3         (k)
      NIntpol4 = r_NIntpol4         (k)
      NIntpol5 = r_NIntpol5         (k)

      do m = 1, KDTblNMol

        if ( mc_IDAbsCoefType(m,iWaveNum) == -1 ) then
          ! zero abs. coef.

          m_Val(m) = 0.0d0

        else
          ! else

          do n5 = 1, NIntpol5
            do n4 = 1, NIntpol4
              do n3 = 1, NIntpol3
                do n2 = 1, NIntpol
                  do n1 = 1, NIntpol1
                    l1 = iPress  +n1-1
                    l2 = iTemp   +n2-1
                    l3 = iVMRMol1+n3-1
                    l4 = iVMRMol2+n4-1
                    l5 = iVMRMol3+n5-1
!!$                aaaam_Val(n1,n2,n3,n4,m) = pmctvv_AbsCoef(l1,m,iWaveNum,l2,l3,l4)
!!$                  ptuvwm_Val(n1,n2,n3,n4,n5,m) = &
!!$                    & log( max( pmctuvw_AbsCoef(l1,m,iWaveNum,l2,l3,l4,l5), 1.0d-100 ) )
                    ptuvw_Val(n1,n2,n3,n4,n5) = &
                      & pmctuvw_LogAbsCoef(l1,m,iWaveNum,l2,l3,l4,l5)
                  end do
                end do
              end do
            end do
          end do

          m_Val(m) = 0.0d0
          do n5 = 1, NIntpol5
            do n4 = 1, NIntpol4
              do n3 = 1, NIntpol3
                do n2 = 1, NIntpol
                  do n1 = 1, NIntpol1
                    m_Val(m) = m_Val(m) &
!!$                      & + ptuvw_Val(n1,n2,n3,n4,n5) * ptuvw_Weight(n1,n2,n3,n4,n5)
                      & + ptuvw_Val(n1,n2,n3,n4,n5) * ptuvwr_Weight(n1,n2,n3,n4,n5,k)
                  end do
                end do
              end do
            end do
          end do

          m_Val(m) = exp( m_Val(m) )

        end if

      end do

      do m = 1, LoopLimit1
        rm_AbsCoef(k,m) = m_Val(m)
      end do
      do m = LoopLimit1+1, LoopLimit2
        rm_AbsCoef(k,m) = 0.0d0
      end do

    end do

    do m = 1, NMol
      select case ( mc_IDAbsCoefType(m,iWaveNum) )
      case ( -1 ) ! zero abs. coef.
        rm_AbsCoef(:,m) = 0.0d0
      case (  0 ) ! abs. coef. per unit abs. molecule
        rm_AbsCoef(:,m) = rm_AbsCoef(:,m) * rm_VMR(:,m)
!!$      case (  1 ) ! abs. coef. per unit atm. molecule
!!$        rm_AbsCoef(:,m) = rm_AbsCoef(:,m)
      end select
    end do

  end subroutine OptPropKDGetAbsCoefProfWithIndex

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

  subroutine CalcWeight( &
    & KDTblNAxis, a_KDTblAxis, iAxis, &
    & TargetValue, &
    & NIntpol, &
    & a_Weight, NIntpolMax &
    & )

    integer, intent(in ) :: KDTblNAxis
    real(8), intent(in ) :: a_KDTblAxis(KDTblNAxis)
    integer, intent(in ) :: iAxis
    real(8), intent(in ) :: TargetValue
    integer, intent(in ) :: NIntpol
    real(8), intent(out) :: a_Weight(NIntpol)
    integer, intent(out) :: NIntpolMax


    ! Local variables
    !
    real(8) :: a_Grid(NIntpol)
    real(8) :: Val
    integer :: l
    integer :: n


    if ( KDTblNAxis > 1 ) then
      do n = 1, NIntpol
        a_Grid(n) = a_KDTblAxis(iAxis+n-1)
      end do
      Val = TargetValue
      if ( a_Grid(1) <= 0.0d0 ) then
        a_Weight = 1.0d0
        do n = 1, NIntpol
          do l = 1, NIntpol
            if ( l == n ) cycle
            a_Weight(n) = a_Weight(n) &
              & * ( Val - a_Grid(l) ) / ( a_Grid(n) - a_Grid(l) )
          end do
        end do
      else
        a_Weight = 1.0d0
        do n = 1, NIntpol
          do l = 1, NIntpol
            if ( l == n ) cycle
            a_Weight(n) = a_Weight(n) &
              & * log( Val / a_Grid(l) ) / log( a_Grid(n) / a_Grid(l) )
          end do
        end do
      end if
      NIntpolMax  = NIntpol
    else
      a_Weight    = 0.0d0
      a_Weight(1) = 1.0d0
      NIntpolMax  = 1
    end if


  end subroutine CalcWeight

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

  subroutine OptPropKDGetRayScatCoef( &
    & iWaveNum, &
    & NMol, m_MolNum, &
    & RayScatCoefNonRadAct, m_RayScatCoef &
    & )


    integer     , intent(in ) :: iWaveNum
    integer     , intent(in ) :: NMol
    integer     , intent(in ) :: m_MolNum(NMol)
    real(DP)    , intent(out) :: RayScatCoefNonRadAct
    real(DP)    , intent(out) :: m_RayScatCoef(NMol)


    ! local variables
    !
    integer :: iBand

    integer :: m


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


    ! check size of VMR of argument and number of molecules in table
    if ( NMol < KDTblNMol ) then
      write( 6, * ) 'Array size for VMR is not enough: ', NMol, KDTblNMol
      stop
    end if
    do m = 1, NMol
      if ( m_MolNum(m) /= m_KDTblMolNum(m) ) then
        write( 6, * ) 'Given molecular number, ', m_MolNum(m), ',  does not match with that in table.'
        write( 6, * ) '  Argument: ', m_MolNum
        write( 6, * ) '  Table   : ', m_KDTblMolNum
        stop
      end if
    end do


    iBand = OptPropKDGetBandNum( iWaveNum )
    RayScatCoefNonRadAct = b_BandAveRayScatCoefNonRadAct( iBand )
    do m = 1, NMol
      m_RayScatCoef(m) = bm_BandAveRayScatCoef( iBand, m )
    end do


  end subroutine OptPropKDGetRayScatCoef

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

  subroutine OptPropKDGetPtclParam( &
    & iWaveNum, &
    & kmax, NPtcl, a_PtclName, za_PtclRadius, &
    & za_QExt, za_SSA, za_AF &
    & )


    integer     , intent(in ) :: iWaveNum
    integer     , intent(in ) :: kmax
    integer     , intent(in ) :: NPtcl
    character(*), intent(in ) :: a_PtclName(NPtcl)
    real(DP)    , intent(in ) :: za_PtclRadius(kmax,NPtcl)
    real(DP)    , intent(out) :: za_QExt(kmax,NPtcl)
    real(DP)    , intent(out) :: za_SSA (kmax,NPtcl)
    real(DP)    , intent(out) :: za_AF  (kmax,NPtcl)


    ! local variables
    !
    integer  :: NRadius
    real(DP) :: Radius
    real(DP) :: a_X(2)
    real(DP) :: a_Val(2)

    integer  :: iBand
    integer  :: iPtcl
    integer  :: i
    integer  :: k
    integer  :: l
    integer  :: n


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


    iBand = OptPropKDGetBandNum( iWaveNum )

    i = iBand

    do n = 1, NPtcl

      do iPtcl = 1, KDTblNPtcl
        if ( trim( a_PtclName(n) ) == trim( a_KDTblPtclName(iPtcl) ) ) exit
      end do
      if ( iPtcl > KDTblNPtcl ) then
        write( 6, * ) trim( a_PtclName(n) ) // ' is not found.'
        stop
      end if

      NRadius = a_KDTblPtclNRadius(iPtcl)

      if ( NRadius == 1 ) then
        l = 1
        do k = 1, kmax
          za_QExt(k,n) = bra_KDTblBandAveQExt(i,l,iPtcl)
          za_SSA (k,n) = bra_KDTblBandAveSSA (i,l,iPtcl)
          za_AF  (k,n) = bra_KDTblBandAveAF  (i,l,iPtcl)
        end do
      else
        do k = 1, kmax
          Radius = za_PtclRadius(k,n)
          if ( Radius <= ra_KDTblPtclRadius(1,n) ) then
            l = 1
            za_QExt(k,n) = bra_KDTblBandAveQExt(i,l,iPtcl)
            za_SSA (k,n) = bra_KDTblBandAveSSA (i,l,iPtcl)
            za_AF  (k,n) = bra_KDTblBandAveAF  (i,l,iPtcl)
          else if ( ra_KDTblPtclRadius(NRadius,n) <= Radius ) then
            l = NRadius
            za_QExt(k,n) = bra_KDTblBandAveQExt(i,l,iPtcl)
            za_SSA (k,n) = bra_KDTblBandAveSSA (i,l,iPtcl)
            za_AF  (k,n) = bra_KDTblBandAveAF  (i,l,iPtcl)
          else
            do l = 1+1, NRadius
              if ( Radius <= ra_KDTblPtclRadius(l,n) ) exit
            end do
            a_X(1)   = ra_KDTblPtclRadius(l-1,iPtcl)
            a_X(2)   = ra_KDTblPtclRadius(l  ,iPtcl)
            !
            a_Val(1) = bra_KDTblBandAveQExt(i,l-1,iPtcl)
            a_Val(2) = bra_KDTblBandAveQExt(i,l  ,iPtcl)
            za_QExt(k,n) = &
              &   ( a_Val(2) - a_Val(1) ) &
              & / ( a_X  (2) - a_X  (1) ) &
              & * ( Radius   - a_X  (1) ) &
              & + a_Val(1)
            !
            a_Val(1) = bra_KDTblBandAveSSA (i,l-1,iPtcl)
            a_Val(2) = bra_KDTblBandAveSSA (i,l  ,iPtcl)
            za_SSA (k,n) = &
              &   ( a_Val(2) - a_Val(1) ) &
              & / ( a_X  (2) - a_X  (1) ) &
              & * ( Radius   - a_X  (1) ) &
              & + a_Val(1)
            !
            a_Val(1) = bra_KDTblBandAveAF  (i,l-1,iPtcl)
            a_Val(2) = bra_KDTblBandAveAF  (i,l  ,iPtcl)
            za_AF  (k,n) = &
              &   ( a_Val(2) - a_Val(1) ) &
              & / ( a_X  (2) - a_X  (1) ) &
              & * ( Radius   - a_X  (1) ) &
              & + a_Val(1)
          end if
        end do
      end if
    end do


  end subroutine OptPropKDGetPtclParam

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

!!$  subroutine OLDOptPropKDGetPFIntedRatio( &
!!$    & iWaveNum, kmax, r_Temp, SurfTemp, &
!!$    & r_AtmPFRatio, 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_AtmPFRatio(1:kmax)
!!$    real(DP)    , intent(out) :: SurfPFRatio
!!$
!!$
!!$    ! local variables
!!$    !
!!$    integer        :: r_KDTblTempIndex  (1:kmax)
!!$    !       It should be noticed that NIntpol is defined twice in this module!
!!$    integer , parameter :: NIntpol = 2
!!$!    integer , parameter :: NIntpol = 3
!!$    real(DP)            :: a_Grid  (NIntpol)
!!$    real(DP)            :: a_Val   (NIntpol)
!!$    real(DP)            :: a_Weight(NIntpol)
!!$    real(DP)            :: Val
!!$    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
!!$
!!$
!!$    call FindTblIndex( &
!!$      & 'Temp', &
!!$      & KDTblNTemp, t_KDTblTemp, NIntpol, kmax, r_Temp, &
!!$      & r_KDTblTempIndex &
!!$      & )
!!$    a_SurfTemp = SurfTemp
!!$    call FindTblIndex( &
!!$      & 'SurfTemp', &
!!$      & KDTblNTemp, t_KDTblTemp, NIntpol, 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
!!$
!!$
!!$      do n = 1, NIntpol
!!$        a_Grid(n) = t_KDTblTemp(iTemp+n-1)
!!$      end do
!!$      a_Weight = 1.0d0
!!$      do n = 1, NIntpol
!!$        do l = 1, NIntpol
!!$          if ( l == n ) cycle
!!$          a_Weight(n) = a_Weight(n) &
!!$            & * ( Temp - a_Grid(l) ) / ( a_Grid(n) - a_Grid(l) )
!!$        end do
!!$      end do
!!$
!!$      do n = 1, NIntpol
!!$        l = iTemp+n-1
!!$        a_Val(n) = ct_AtmPFRatio(iWaveNum,l)
!!$      end do
!!$
!!$      Val = 0.0d0
!!$      do n = 1, NIntpol
!!$        Val = Val + a_Val(n) * a_Weight(n)
!!$      end do
!!$
!!$      if ( k == 0 ) then
!!$        SurfPFRatio  = Val
!!$      else
!!$        r_AtmPFRatio(k) = Val
!!$      end if
!!$    end do
!!$
!!$
!!$  end subroutine OLDOptPropKDGetPFIntedRatio

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

  subroutine OptPropKDGetPFIntedRatio( &
    & iWaveNum, kmax, NMol, m_MolNum, r_Press, r_Temp, SurfTemp, rm_VMR, &
    & r_AtmPFRatio, SurfPFRatio &
    & )

    integer     , intent(in ) :: iWaveNum
    integer     , intent(in ) :: kmax
    integer     , intent(in ) :: NMol
    integer     , intent(in ) :: m_MolNum(NMol)
    real(DP)    , intent(in ) :: r_Press  (1:kmax)
    real(DP)    , intent(in ) :: r_Temp   (1:kmax)
    real(DP)    , intent(in ) :: SurfTemp
    real(DP)    , intent(in ) :: rm_VMR    (1:kmax,1:NMol)
    real(DP)    , intent(out) :: r_AtmPFRatio(1:kmax)
    real(DP)    , intent(out) :: SurfPFRatio



    ! local variables
    !
    integer        :: r_KDTblPressIndex  (1:kmax)
    integer        :: r_KDTblTempIndex   (1:kmax)
    integer        :: r_KDTblVMRMol1Index(1:kmax)
    integer        :: r_KDTblVMRMol2Index(1:kmax)
    integer        :: r_KDTblVMRMol3Index(1:kmax)
    !       It should be noticed that NIntpol is defined twice in this module!
    integer , parameter :: NIntpol = 2
!!$    integer , parameter :: NIntpol = 3
    real(DP)            :: a_Grid   (NIntpol)
    real(DP)            :: a_Weight (NIntpol)
    real(DP)            :: p_Weight1(NIntpol)
    real(DP)            :: t_Weight2(NIntpol)
    real(DP)            :: u_Weight3(NIntpol)
    real(DP)            :: v_Weight4(NIntpol)
    real(DP)            :: w_Weight5(NIntpol)
    real(DP)            :: ptuvw_Weight(NIntpol,NIntpol,NIntpol,NIntpol,NIntpol)
    real(DP)            :: ptuvw_Val  (NIntpol,NIntpol,NIntpol,NIntpol,NIntpol)
    integer             :: NIntpol1
    integer             :: NIntpol3
    integer             :: NIntpol4
    integer             :: NIntpol5
    integer             :: IndexVMRMol1Num
    integer             :: IndexVMRMol2Num
    integer             :: IndexVMRMol3Num
    real(DP)       :: Val
    real(DP)       :: Press
    real(DP)       :: Temp
    real(DP)       :: VMRMol1
    real(DP)       :: VMRMol2
    real(DP)       :: VMRMol3
    integer        :: iPress
    integer        :: iTemp
    integer        :: iVMRMol1
    integer        :: iVMRMol2
    integer        :: iVMRMol3
    integer        :: k
    integer        :: l
    integer        :: m
    integer        :: n
    integer        :: l1
    integer        :: l2
    integer        :: l3
    integer        :: l4
    integer        :: l5
    integer        :: n1
    integer        :: n2
    integer        :: n3
    integer        :: n4
    integer        :: n5


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


    ! check size of VMR of argument and number of molecules in table
    if ( NMol < KDTblNMol ) then
      write( 6, * ) 'Array size for VMR is not enough: ', NMol, KDTblNMol
      stop
    end if
    do m = 1, NMol
      if ( m_MolNum(m) /= m_KDTblMolNum(m) ) then
        write( 6, * ) 'Given molecular number, ', m_MolNum(m), ',  does not match with that in table.'
        write( 6, * ) '  Argument: ', m_MolNum
        write( 6, * ) '  Table   : ', m_KDTblMolNum
        stop
      end if
    end do


    if ( KDTblNPress < NIntpol ) then
      write( 6, * ) 'Array size for Press is not enough for interpolation: ', &
        & KDTblNPress, NIntpol
      stop
    end if
    if ( KDTblNTemp < NIntpol ) then
      write( 6, * ) 'Array size for Temp is not enough for interpolation: ', &
        & KDTblNTemp, NIntpol
      stop
    end if
    if ( KDTblNVMRMol1 > 1 ) then
      if ( KDTblNVMRMol1 < NIntpol ) then
        write( 6, * ) 'Array size for VMRMol1 is not enough for interpolation: ', &
          & KDTblNVMRMol1, NIntpol
        stop
      end if
    end if
    if ( KDTblNVMRMol2 > 1 ) then
      if ( KDTblNVMRMol2 < NIntpol ) then
        write( 6, * ) 'Array size for VMRMol2 is not enough for interpolation: ', &
          & KDTblNVMRMol2, NIntpol
        stop
      end if
    end if
    if ( KDTblNVMRMol3 > 1 ) then
      if ( KDTblNVMRMol3 < NIntpol ) then
        write( 6, * ) 'Array size for VMRMol3 is not enough for interpolation: ', &
          & KDTblNVMRMol3, NIntpol
        stop
      end if
    end if


    do m = 1, NMol
      if ( m_MolNum(m) == m_KDTblMolNum(KDTblIndexVMRMol1Num) ) exit
    end do
    if ( m > NMol ) then
      write( 6, * ) &
        & 'Unable to find a molecular number, ', &
        & m_KDTblMolNum(KDTblIndexVMRMol1Num), ', in MolNum in argument.'
      stop
    end if
    IndexVMRMol1Num = m
    !
    do m = 1, NMol
      if ( m_MolNum(m) == m_KDTblMolNum(KDTblIndexVMRMol2Num) ) exit
    end do
    if ( m > NMol ) then
      write( 6, * ) &
        & 'Unable to find a molecular number, ', &
        & m_KDTblMolNum(KDTblIndexVMRMol2Num), ', in MolNum in argument.'
      stop
    end if
    IndexVMRMol2Num = m
    !
    do m = 1, NMol
      if ( m_MolNum(m) == m_KDTblMolNum(KDTblIndexVMRMol3Num) ) exit
    end do
    if ( m > NMol ) then
      write( 6, * ) &
        & 'Unable to find a molecular number, ', &
        & m_KDTblMolNum(KDTblIndexVMRMol3Num), ', in MolNum in argument.'
      stop
    end if
    IndexVMRMol3Num = m


    call FindTblIndex( &
      & 'Press', &
      & KDTblNPress, p_KDTblPress, NIntpol, kmax, r_Press, &
      & r_KDTblPressIndex &
      & )
    call FindTblIndex( &
      & 'Temp', &
      & KDTblNTemp, t_KDTblTemp, NIntpol, kmax, r_Temp, &
      & r_KDTblTempIndex &
      & )
    if ( KDTblNVMRMol1 > 1 ) then
      call FindTblIndex( &
        & 'VMRMol1', &
        & KDTblNVMRMol1, u_KDTblVMRMol1, NIntpol, kmax, rm_VMR(:,IndexVMRMol1Num), &
        & r_KDTblVMRMol1Index &
        & )
    else
      r_KDTblVMRMol1Index = 1
    end if
    if ( KDTblNVMRMol2 > 1 ) then
      call FindTblIndex( &
        & 'VMRMol2', &
        & KDTblNVMRMol2, v_KDTblVMRMol2, NIntpol, kmax, rm_VMR(:,IndexVMRMol2Num), &
        & r_KDTblVMRMol2Index &
        & )
    else
      r_KDTblVMRMol2Index = 1
    end if
    if ( KDTblNVMRMol3 > 1 ) then
      call FindTblIndex( &
        & 'VMRMol3', &
        & KDTblNVMRMol3, w_KDTblVMRMol3, NIntpol, kmax, rm_VMR(:,IndexVMRMol3Num), &
        & r_KDTblVMRMol3Index &
        & )
    else
      r_KDTblVMRMol3Index = 1
    end if


    do k = 1, kmax
      Press   = r_Press (k)
      Temp    = r_Temp  (k)
      VMRMol1 = rm_VMR  (k,IndexVMRMol1Num)
      VMRMol2 = rm_VMR  (k,IndexVMRMol2Num)
      VMRMol3 = rm_VMR  (k,IndexVMRMol3Num)

      iPress   = r_KDTblPressIndex  (k)
      iTemp    = r_KDTblTempIndex   (k)
      iVMRMol1 = r_KDTblVMRMol1Index(k)
      iVMRMol2 = r_KDTblVMRMol2Index(k)
      iVMRMol3 = r_KDTblVMRMol3Index(k)


      call CalcWeight( &
        & KDTblNPress, p_KDTblPress, iPress, &
        & Press, &
        & NIntpol, &
        & p_Weight1, NIntpol1 &
        & )
      !
      do n = 1, NIntpol
        a_Grid(n) = t_KDTblTemp(iTemp+n-1)
      end do
      Val = Temp
      a_Weight = 1.0d0
      do n = 1, NIntpol
        do l = 1, NIntpol
          if ( l == n ) cycle
          a_Weight(n) = a_Weight(n) &
            & * ( Val - a_Grid(l) ) / ( a_Grid(n) - a_Grid(l) )
        end do
      end do
      t_Weight2 = a_Weight
      !
      call CalcWeight( &
        & KDTblNVMRMol1, u_KDTblVMRMol1, iVMRMol1, &
        & VMRMol1, &
        & NIntpol, &
        & u_Weight3, NIntpol3 &
        & )
      !
      call CalcWeight( &
        & KDTblNVMRMol2, v_KDTblVMRMol2, iVMRMol2, &
        & VMRMol2, &
        & NIntpol, &
        & v_Weight4, NIntpol4 &
        & )
      !
      call CalcWeight( &
        & KDTblNVMRMol3, w_KDTblVMRMol3, iVMRMol3, &
        & VMRMol3, &
        & NIntpol, &
        & w_Weight5, NIntpol5 &
        & )
      !
      ptuvw_Weight = 1.0d0
      do n5 = 1, NIntpol5
        do n4 = 1, NIntpol4
          do n3 = 1, NIntpol3
            do n2 = 1, NIntpol
!!$              do n1 = 1, NIntpol
              do n1 = 1, NIntpol1
                ptuvw_Weight(n1,n2,n3,n4,n5) = ptuvw_Weight(n1,n2,n3,n4,n5) &
                  & * p_Weight1(n1) * t_Weight2(n2) &
                  & * u_Weight3(n3) * v_Weight4(n4) * w_Weight5(n5)
              end do
            end do
          end do
        end do
      end do

      do n5 = 1, NIntpol5
        do n4 = 1, NIntpol4
          do n3 = 1, NIntpol3
            do n2 = 1, NIntpol
!!$              do n1 = 1, NIntpol
              do n1 = 1, NIntpol1
                l1 = iPress  +n1-1
                l2 = iTemp   +n2-1
                l3 = iVMRMol1+n3-1
                l4 = iVMRMol2+n4-1
                l5 = iVMRMol3+n5-1
                ptuvw_Val(n1,n2,n3,n4,n5) = &
                  & pctuvw_AtmPFRatio(l1,iWaveNum,l2,l3,l4,l5)
              end do
            end do
          end do
        end do
      end do

      Val = 0.0d0
      do n5 = 1, NIntpol5
        do n4 = 1, NIntpol4
          do n3 = 1, NIntpol3
            do n2 = 1, NIntpol
!!$              do n1 = 1, NIntpol
              do n1 = 1, NIntpol1
                Val = Val &
                  & + ptuvw_Val(n1,n2,n3,n4,n5) * ptuvw_Weight(n1,n2,n3,n4,n5)
              end do
            end do
          end do
        end do

      end do

      r_AtmPFRatio(k) = Val

    end do



    call FindTblIndex( &
      & 'Temp', &
      & KDTblNTemp, t_KDTblTemp, NIntpol, 1, (/SurfTemp/), &
      & r_KDTblTempIndex(1:1) &
      & )
    do k = 1, 1
      Press   = r_Press (k)
      Temp    = SurfTemp
      VMRMol1 = rm_VMR  (k,IndexVMRMol1Num)
      VMRMol2 = rm_VMR  (k,IndexVMRMol2Num)
      VMRMol3 = rm_VMR  (k,IndexVMRMol3Num)

      iPress   = r_KDTblPressIndex  (k)
      iTemp    = r_KDTblTempIndex   (k)
      iVMRMol1 = r_KDTblVMRMol1Index(k)
      iVMRMol2 = r_KDTblVMRMol2Index(k)
      iVMRMol3 = r_KDTblVMRMol3Index(k)


      call CalcWeight( &
        & KDTblNPress, p_KDTblPress, iPress, &
        & Press, &
        & NIntpol, &
        & p_Weight1, NIntpol1 &
        & )
      !
      do n = 1, NIntpol
        a_Grid(n) = t_KDTblTemp(iTemp+n-1)
      end do
      Val = Temp
      a_Weight = 1.0d0
      do n = 1, NIntpol
        do l = 1, NIntpol
          if ( l == n ) cycle
          a_Weight(n) = a_Weight(n) &
            & * ( Val - a_Grid(l) ) / ( a_Grid(n) - a_Grid(l) )
        end do
      end do
      t_Weight2 = a_Weight
      !
      call CalcWeight( &
        & KDTblNVMRMol1, u_KDTblVMRMol1, iVMRMol1, &
        & VMRMol1, &
        & NIntpol, &
        & u_Weight3, NIntpol3 &
        & )
      !
      call CalcWeight( &
        & KDTblNVMRMol2, v_KDTblVMRMol2, iVMRMol2, &
        & VMRMol2, &
        & NIntpol, &
        & v_Weight4, NIntpol4 &
        & )
      !
      call CalcWeight( &
        & KDTblNVMRMol3, w_KDTblVMRMol3, iVMRMol3, &
        & VMRMol3, &
        & NIntpol, &
        & w_Weight5, NIntpol5 &
        & )
      !
      ptuvw_Weight = 1.0d0
      do n5 = 1, NIntpol5
        do n4 = 1, NIntpol4
          do n3 = 1, NIntpol3
            do n2 = 1, NIntpol
!!$              do n1 = 1, NIntpol
              do n1 = 1, NIntpol1
                ptuvw_Weight(n1,n2,n3,n4,n5) = ptuvw_Weight(n1,n2,n3,n4,n5) &
                  & * p_Weight1(n1) * t_Weight2(n2) &
                  & * u_Weight3(n3) * v_Weight4(n4) * w_Weight5(n5)
              end do
            end do
          end do
        end do
      end do

      do n5 = 1, NIntpol5
        do n4 = 1, NIntpol4
          do n3 = 1, NIntpol3
            do n2 = 1, NIntpol
!!$              do n1 = 1, NIntpol
              do n1 = 1, NIntpol1
                l1 = iPress  +n1-1
                l2 = iTemp   +n2-1
                l3 = iVMRMol1+n3-1
                l4 = iVMRMol2+n4-1
                l5 = iVMRMol3+n5-1
                ptuvw_Val(n1,n2,n3,n4,n5) = &
                  & pctuvw_AtmPFRatio(l1,iWaveNum,l2,l3,l4,l5)
              end do
            end do
          end do
        end do
      end do

      Val = 0.0d0
      do n5 = 1, NIntpol5
        do n4 = 1, NIntpol4
          do n3 = 1, NIntpol3
            do n2 = 1, NIntpol
!!$              do n1 = 1, NIntpol
              do n1 = 1, NIntpol1
                Val = Val &
                  & + ptuvw_Val(n1,n2,n3,n4,n5) * ptuvw_Weight(n1,n2,n3,n4,n5)
              end do
            end do
          end do
        end do

      end do

      SurfPFRatio = Val

    end do


  end subroutine OptPropKDGetPFIntedRatio

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

  subroutine OptPropKDGetPFIntedRatioWithIndex( &
    & iWaveNum, kmax, NMol, m_MolNum, r_Press, r_Temp, SurfTemp, rm_VMR, &
    & r_AtmPFRatio, SurfPFRatio, &
    & 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 &
    & )

    integer     , intent(in ) :: iWaveNum
    integer     , intent(in ) :: kmax
    integer     , intent(in ) :: NMol
    integer     , intent(in ) :: m_MolNum(NMol)
    real(DP)    , intent(in ) :: r_Press  (1:kmax)
    real(DP)    , intent(in ) :: r_Temp   (1:kmax)
    real(DP)    , intent(in ) :: SurfTemp
    real(DP)    , intent(in ) :: rm_VMR    (1:kmax,1:NMol)
    real(DP)    , intent(out) :: r_AtmPFRatio(1:kmax)
    real(DP)    , intent(out) :: SurfPFRatio
    integer     , intent(in ) :: NIntpol
    integer     , intent(in ) :: r_KDTblPressIndex  (1:kmax)
    integer     , intent(in ) :: r_KDTblTempIndex   (1:kmax)
    integer     , intent(in ) :: r_KDTblVMRMol1Index(1:kmax)
    integer     , intent(in ) :: r_KDTblVMRMol2Index(1:kmax)
    integer     , intent(in ) :: r_KDTblVMRMol3Index(1:kmax)
    integer     , intent(in ) :: r_NIntpol1         (1:kmax)
    integer     , intent(in ) :: r_NIntpol3         (1:kmax)
    integer     , intent(in ) :: r_NIntpol4         (1:kmax)
    integer     , intent(in ) :: r_NIntpol5         (1:kmax)
    real(DP)    , intent(in ) :: ptuvwr_Weight(NIntpol,NIntpol,NIntpol,NIntpol,NIntpol,1:kmax)
    integer     , intent(in ) :: r_KDTblPressIndexSfc  (1)
    integer     , intent(in ) :: r_KDTblTempIndexSfc   (1)
    integer     , intent(in ) :: r_KDTblVMRMol1IndexSfc(1)
    integer     , intent(in ) :: r_KDTblVMRMol2IndexSfc(1)
    integer     , intent(in ) :: r_KDTblVMRMol3IndexSfc(1)
    integer     , intent(in ) :: r_NIntpol1Sfc         (1)
    integer     , intent(in ) :: r_NIntpol3Sfc         (1)
    integer     , intent(in ) :: r_NIntpol4Sfc         (1)
    integer     , intent(in ) :: r_NIntpol5SFc         (1)
    real(DP)    , intent(in ) :: ptuvw_WeightSfc(NIntpol,NIntpol,NIntpol,NIntpol,NIntpol)



    ! local variables
    !
    integer             :: r_KDTblTempIndex1(1)
    real(DP)            :: ptuvw_Weight(NIntpol,NIntpol,NIntpol,NIntpol,NIntpol)
    real(DP)            :: ptuvw_Val  (NIntpol,NIntpol,NIntpol,NIntpol,NIntpol)
    integer             :: NIntpol1
    integer             :: NIntpol3
    integer             :: NIntpol4
    integer             :: NIntpol5
    integer             :: IndexVMRMol1Num
    integer             :: IndexVMRMol2Num
    integer             :: IndexVMRMol3Num
    real(DP)       :: Val
    real(DP)       :: Press
    real(DP)       :: Temp
    real(DP)       :: VMRMol1
    real(DP)       :: VMRMol2
    real(DP)       :: VMRMol3
    integer        :: iPress
    integer        :: iTemp
    integer        :: iVMRMol1
    integer        :: iVMRMol2
    integer        :: iVMRMol3
    integer        :: k
    integer        :: l
    integer        :: m
    integer        :: n
    integer        :: l1
    integer        :: l2
    integer        :: l3
    integer        :: l4
    integer        :: l5
    integer        :: n1
    integer        :: n2
    integer        :: n3
    integer        :: n4
    integer        :: n5


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



    do k = 1, kmax

      iPress   = r_KDTblPressIndex  (k)
      iTemp    = r_KDTblTempIndex   (k)
      iVMRMol1 = r_KDTblVMRMol1Index(k)
      iVMRMol2 = r_KDTblVMRMol2Index(k)
      iVMRMol3 = r_KDTblVMRMol3Index(k)

      NIntpol1 = r_NIntpol1         (k)
      NIntpol3 = r_NIntpol3         (k)
      NIntpol4 = r_NIntpol4         (k)
      NIntpol5 = r_NIntpol5         (k)


      do n5 = 1, NIntpol5
        do n4 = 1, NIntpol4
          do n3 = 1, NIntpol3
            do n2 = 1, NIntpol
!!$              do n1 = 1, NIntpol
              do n1 = 1, NIntpol1
                l1 = iPress  +n1-1
                l2 = iTemp   +n2-1
                l3 = iVMRMol1+n3-1
                l4 = iVMRMol2+n4-1
                l5 = iVMRMol3+n5-1
                ptuvw_Val(n1,n2,n3,n4,n5) = &
                  & pctuvw_AtmPFRatio(l1,iWaveNum,l2,l3,l4,l5)
              end do
            end do
          end do
        end do
      end do

      Val = 0.0d0
      do n5 = 1, NIntpol5
        do n4 = 1, NIntpol4
          do n3 = 1, NIntpol3
            do n2 = 1, NIntpol
!!$              do n1 = 1, NIntpol
              do n1 = 1, NIntpol1
                Val = Val &
                  & + ptuvw_Val(n1,n2,n3,n4,n5) * ptuvwr_Weight(n1,n2,n3,n4,n5,k)
              end do
            end do
          end do
        end do

      end do

      r_AtmPFRatio(k) = Val

    end do



    do k = 1, 1
      iPress   = r_KDTblPressIndexSfc  (k)
      iTemp    = r_KDTblTempIndexSfc   (k)
      iVMRMol1 = r_KDTblVMRMol1IndexSfc(k)
      iVMRMol2 = r_KDTblVMRMol2IndexSfc(k)
      iVMRMol3 = r_KDTblVMRMol3IndexSfc(k)

      NIntpol1 = r_NIntpol1Sfc         (k)
      NIntpol3 = r_NIntpol3Sfc         (k)
      NIntpol4 = r_NIntpol4Sfc         (k)
      NIntpol5 = r_NIntpol5Sfc         (k)


      do n5 = 1, NIntpol5
        do n4 = 1, NIntpol4
          do n3 = 1, NIntpol3
            do n2 = 1, NIntpol
!!$              do n1 = 1, NIntpol
              do n1 = 1, NIntpol1
                l1 = iPress  +n1-1
                l2 = iTemp   +n2-1
                l3 = iVMRMol1+n3-1
                l4 = iVMRMol2+n4-1
                l5 = iVMRMol3+n5-1
                ptuvw_Val(n1,n2,n3,n4,n5) = &
                  & pctuvw_AtmPFRatio(l1,iWaveNum,l2,l3,l4,l5)
              end do
            end do
          end do
        end do
      end do

      Val = 0.0d0
      do n5 = 1, NIntpol5
        do n4 = 1, NIntpol4
          do n3 = 1, NIntpol3
            do n2 = 1, NIntpol
!!$              do n1 = 1, NIntpol
              do n1 = 1, NIntpol1
                Val = Val &
                  & + ptuvw_Val(n1,n2,n3,n4,n5) * ptuvw_WeightSfc(n1,n2,n3,n4,n5)
              end do
            end do
          end do
        end do

      end do

      SurfPFRatio = Val

    end do


  end subroutine OptPropKDGetPFIntedRatioWithIndex

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

  subroutine OptPropKDGetStrPFIntedRatio( &
    & iWaveNum, &
    & StrPFRatio &
    & )


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


    ! local variables
    !

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

    StrPFRatio = c_StrPFRatio(iWaveNum)


  end subroutine OptPropKDGetStrPFIntedRatio

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

!!$  subroutine OptPropKDPFTblInit
!!$
!!$    use constants0, only : PI
!!$    use rad_kd_utils, only : InqPFIntGQNum
!!$
!!$    use planck_func_wrapper, only : &
!!$      & Integ_PF_GQ_Array1D   , &
!!$      & Integ_DPFDT_GQ_Array1D
!!$
!!$    ! Local variables
!!$    !
!!$    real(dp), parameter :: DelTemp = 1.0d0
!!$
!!$    real(dp) :: WNs
!!$    real(dp) :: WNe
!!$
!!$    integer :: GQNum
!!$
!!$    integer :: iBand
!!$    integer :: iTemp
!!$
!!$
!!$    write( 6, * ) 'Preparing Planck function table ...'
!!$
!!$    ! Last +1 is preparation for quadratic interpolation.
!!$    PFTblNTemp = int( ( t_KDTblTemp(KDTblNTemp) - t_KDTblTemp(1) ) / DelTemp ) + 1 + 1
!!$    PFTblNBand = KDTblNBand
!!$
!!$    allocate( t_PFTblTemp   (PFTblNTemp) )
!!$!    allocate( b_PFTblBandNum(PFTblNBand) )
!!$    allocate( tb_PFTblPF    (PFTblNTemp, PFTblNBand) )
!!$    allocate( tb_PFTblDPFDT (PFTblNTemp, PFTblNBand) )
!!$
!!$    do iTemp = 1, PFTblNTemp
!!$      t_PFTblTemp(iTemp) = t_KDTblTemp(1) + DelTemp * (iTemp-1)
!!$    end do
!!$
!!$    do iBand = 1, PFTblNBand
!!$      WNs = ab_KDTblBandWaveNumBnds(1,iBand)
!!$      WNe = ab_KDTblBandWaveNumBnds(2,iBand)
!!$      !
!!$      GQNum = InqPFIntGQNum( WNe - WNs )
!!$      !
!!$      call Integ_PF_GQ_Array1D( &
!!$        & WNs, WNe, GQNum, &
!!$        & 1, PFTblNTemp, &
!!$        & t_PFTblTemp, &
!!$        & tb_PFTblPF(:,iBand) &
!!$        & )
!!$      tb_PFTblPF(:,iBand) = PI * tb_PFTblPF(:,iBand) / ( WNe - WNs )
!!$      call Integ_DPFDT_GQ_Array1D( &
!!$        & WNs, WNe, GQNum,         & ! (in )
!!$        & 1, PFTblNTemp,           &
!!$        & t_PFTblTemp,             &
!!$        & tb_PFTblDPFDT(:,iBand) &
!!$        & )
!!$      tb_PFTblDPFDT(:,iBand) = PI * tb_PFTblDPFDT(:,iBand) / ( WNe - WNs )
!!$    end do
!!$
!!$
!!$  end subroutine OptPropKDPFTblInit

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

  subroutine OptPropKDPFTblFinalize

    ! Local variables
    !


    deallocate( t_PFTblTemp    )
!!$    deallocate( b_PFTblBandNum )
    deallocate( tb_PFTblPF     )
    deallocate( tb_PFTblDPFDT  )

    deallocate( t_PFTblPFSum    )
    deallocate( tb_PFTblPFRatio )

  end subroutine OptPropKDPFTblFinalize

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

  subroutine OptPropKDPFTblGetPFDPFDT2DV2( &
    & iBand, &
    & imax, &
    & is, ie, &
    & NLev, xa_Temp, &
    & xa_PF, &
    & FlagDPFDT &
    & )

    integer , intent(in ) :: iBand
    integer , intent(in ) :: imax
    integer , intent(in ) :: is
    integer , intent(in ) :: ie
    integer , intent(in ) :: NLev
    real(dp), intent(in ) :: xa_Temp(imax,NLev)
    real(dp), intent(out) :: xa_PF  (imax,NLev)
    logical , intent(in ), optional :: FlagDPFDT


    ! Local variables
    !
    real(dp) :: a_Temp(NLev)
    real(dp) :: a_PF  (NLev)

    real(dp) :: DelTemp
    integer  :: iTemp
    integer  :: a_iTemp(NLev)

    real(dp) :: XVal
    real(dp) :: a_X(3)
    real(dp) :: a_Y(3)

    integer  :: i
    integer  :: k


    do i = is, ie
      do k = 1, NLev
        a_Temp(k) = xa_Temp(i,k)
      end do


      DelTemp = t_PFTblTemp(2) - t_PFTblTemp(1)

      do k = 1, NLev
        if (      ( a_Temp(k) < t_PFTblTemp(1)          ) &
          &  .or. ( t_PFTblTemp(PFTblNTemp) < a_Temp(k) ) ) then
          write( 6, * ) 'PF Interpolation: Temperature is out of range:'
          write( 6, * ) '     ', k, a_Temp(k)
          write( 6, * ) '     ', t_PFTblTemp(1), t_PFTblTemp(PFTblNTemp)
          stop
        end if
      end do

      iTemp = int( ( a_Temp(k) - t_PFTblTemp(1) ) / DelTemp ) + 1 + 1
      iTemp = max( iTemp, 2            )
      iTemp = min( iTemp, PFTblNTemp-1 )

      if (      ( a_Temp(k) < t_PFTblTemp(iTemp-1) ) &
        &  .or. ( t_PFTblTemp(iTemp+1) < a_Temp(k) ) ) then
        write( 6, * ) 'PF Interpolation: Deduced position is out of range:'
        write( 6, * ) '     ', k, a_Temp(k)
        write( 6, * ) '     ', t_PFTblTemp(iTemp-1), t_PFTblTemp(iTemp+1)
        stop
      end if

      a_iTemp(k) = iTemp


      if ( present( FlagDPFDT ) .and. FlagDPFDT ) then
        do k = 1, NLev
          iTemp = a_iTemp(k)

!!$        a_PF(k) = &
!!$          &   ( tb_PFTblDPFDT(iTemp,iBand) - tb_PFTblDPFDT(iTemp-1,iBand) ) &
!!$          & / ( t_PFTblTemp(iTemp)         - t_PFTblTemp(iTemp-1)         ) &
!!$          & * ( a_Temp(k)                  - t_PFTblTemp(iTemp-1)         ) &
!!$          & + tb_PFTblDPFDT(iTemp-1,iBand)

          XVal   = a_Temp(k)
          a_X(1) = t_PFTblTemp(iTemp-1)
          a_X(2) = t_PFTblTemp(iTemp  )
          a_X(3) = t_PFTblTemp(iTemp+1)
          a_Y(1) = tb_PFTblDPFDT(iTemp-1,iBand)
          a_Y(2) = tb_PFTblDPFDT(iTemp  ,iBand)
          a_Y(3) = tb_PFTblDPFDT(iTemp+1,iBand)
          a_PF(k) = &
            &   ( ( XVal   - a_X(2) ) * ( XVal   - a_X(3) ) ) &
            & / ( ( a_X(1) - a_X(2) ) * ( a_X(1) - a_X(3) ) ) &
            & * a_Y(1)                                        &
            & + ( ( XVal   - a_X(3) ) * ( XVal   - a_X(1) ) ) &
            & / ( ( a_X(2) - a_X(3) ) * ( a_X(2) - a_X(1) ) ) &
            & * a_Y(2)                                        &
            & + ( ( XVal   - a_X(1) ) * ( XVal   - a_X(2) ) ) &
            & / ( ( a_X(3) - a_X(1) ) * ( a_X(3) - a_X(2) ) ) &
            & * a_Y(3)
        end do
      else
        do k = 1, NLev
          iTemp = a_iTemp(k)

!!$        a_PF(k) = &
!!$          &   ( tb_PFTblPF(iTemp,iBand) - tb_PFTblPF(iTemp-1,iBand) ) &
!!$          & / ( t_PFTblTemp(iTemp)      - t_PFTblTemp(iTemp-1)      ) &
!!$          & * ( a_Temp(k)               - t_PFTblTemp(iTemp-1)      ) &
!!$          & + tb_PFTblPF(iTemp-1,iBand)

          XVal   = a_Temp(k)
          a_X(1) = t_PFTblTemp(iTemp-1)
          a_X(2) = t_PFTblTemp(iTemp  )
          a_X(3) = t_PFTblTemp(iTemp+1)
          a_Y(1) = tb_PFTblPF(iTemp-1,iBand)
          a_Y(2) = tb_PFTblPF(iTemp  ,iBand)
          a_Y(3) = tb_PFTblPF(iTemp+1,iBand)
          a_PF(k) = &
            &   ( ( XVal   - a_X(2) ) * ( XVal   - a_X(3) ) ) &
            & / ( ( a_X(1) - a_X(2) ) * ( a_X(1) - a_X(3) ) ) &
            & * a_Y(1)                                        &
            & + ( ( XVal   - a_X(3) ) * ( XVal   - a_X(1) ) ) &
            & / ( ( a_X(2) - a_X(3) ) * ( a_X(2) - a_X(1) ) ) &
            & * a_Y(2)                                        &
            & + ( ( XVal   - a_X(1) ) * ( XVal   - a_X(2) ) ) &
            & / ( ( a_X(3) - a_X(1) ) * ( a_X(3) - a_X(2) ) ) &
            & * a_Y(3)
        end do
      end if

      do k = 1, NLev
        xa_PF(i,k) = a_PF(k)
      end do

    end do


  end subroutine OptPropKDPFTblGetPFDPFDT2DV2

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

  subroutine OptPropKDPFTblGetPFDPFDT1DV2( &
    & iBand, &
    & imax, &
    & is, ie, &
    & x_Temp, &
    & x_PF, &
    & FlagDPFDT &
    & )

    integer , intent(in ) :: iBand
    integer , intent(in ) :: imax
    integer , intent(in ) :: is
    integer , intent(in ) :: ie
    real(dp), intent(in ) :: x_Temp(imax)
    real(dp), intent(out) :: x_PF  (imax)
    logical , intent(in ), optional :: FlagDPFDT


    ! Local variables
    !
    real(dp) :: xa_Temp(imax,1)
    real(dp) :: xa_PF  (imax,1)

    integer  :: i


    do i = is, ie
      xa_Temp(i,1) = x_Temp(i)
    end do

    call OptPropKDPFTblGetPFDPFDT( &
      & iBand, &
      & imax, &
      & is, ie, &
      & 1, xa_Temp, &
      & xa_PF, &
      & FlagDPFDT &
      & )

    do i = is, ie
      x_PF(i) = xa_PF(i,1)
    end do


  end subroutine OptPropKDPFTblGetPFDPFDT1DV2

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

  subroutine OptPropKDPFTblGetPFDPFDT1D( &
    & iBand, NLev, a_Temp, &
    & a_PF, &
    & FlagDPFDT &
    & )

    integer , intent(in ) :: iBand
    integer , intent(in ) :: NLev
    real(dp), intent(in ) :: a_Temp(NLev)
    real(dp), intent(out) :: a_PF  (NLev)
    logical , intent(in ), optional :: FlagDPFDT


    ! Local variables
    !
    real(dp) :: DelTemp
    integer  :: iTemp
    integer  :: a_iTemp(NLev)

    real(dp) :: XVal
    real(dp) :: a_X(3)
    real(dp) :: a_Y(3)

    integer  :: k


    DelTemp = t_PFTblTemp(2) - t_PFTblTemp(1)

    do k = 1, NLev

      if (      ( a_Temp(k) < t_PFTblTemp(1)          ) &
        &  .or. ( t_PFTblTemp(PFTblNTemp) < a_Temp(k) ) ) then
        write( 6, * ) 'PF Interpolation: Temperature is out of range:'
        write( 6, * ) '     ', k, a_Temp(k)
        write( 6, * ) '     ', t_PFTblTemp(1), t_PFTblTemp(PFTblNTemp)
        stop
      end if

      iTemp = int( ( a_Temp(k) - t_PFTblTemp(1) ) / DelTemp ) + 1 + 1
      iTemp = max( iTemp, 2            )
      iTemp = min( iTemp, PFTblNTemp-1 )

      if (      ( a_Temp(k) < t_PFTblTemp(iTemp-1) ) &
        &  .or. ( t_PFTblTemp(iTemp+1) < a_Temp(k) ) ) then
        write( 6, * ) 'PF Interpolation: Deduced position is out of range:'
        write( 6, * ) '     ', k, a_Temp(k)
        write( 6, * ) '     ', t_PFTblTemp(iTemp-1), t_PFTblTemp(iTemp+1)
        stop
      end if

      a_iTemp(k) = iTemp

    end do

    if ( present( FlagDPFDT ) .and. FlagDPFDT ) then
      do k = 1, NLev
        iTemp = a_iTemp(k)

!!$        a_PF(k) = &
!!$          &   ( tb_PFTblDPFDT(iTemp,iBand) - tb_PFTblDPFDT(iTemp-1,iBand) ) &
!!$          & / ( t_PFTblTemp(iTemp)         - t_PFTblTemp(iTemp-1)         ) &
!!$          & * ( a_Temp(k)                  - t_PFTblTemp(iTemp-1)         ) &
!!$          & + tb_PFTblDPFDT(iTemp-1,iBand)

        XVal   = a_Temp(k)
        a_X(1) = t_PFTblTemp(iTemp-1)
        a_X(2) = t_PFTblTemp(iTemp  )
        a_X(3) = t_PFTblTemp(iTemp+1)
        a_Y(1) = tb_PFTblDPFDT(iTemp-1,iBand)
        a_Y(2) = tb_PFTblDPFDT(iTemp  ,iBand)
        a_Y(3) = tb_PFTblDPFDT(iTemp+1,iBand)
        a_PF(k) = &
          &   ( ( XVal   - a_X(2) ) * ( XVal   - a_X(3) ) ) &
          & / ( ( a_X(1) - a_X(2) ) * ( a_X(1) - a_X(3) ) ) &
          & * a_Y(1)                                        &
          & + ( ( XVal   - a_X(3) ) * ( XVal   - a_X(1) ) ) &
          & / ( ( a_X(2) - a_X(3) ) * ( a_X(2) - a_X(1) ) ) &
          & * a_Y(2)                                        &
          & + ( ( XVal   - a_X(1) ) * ( XVal   - a_X(2) ) ) &
          & / ( ( a_X(3) - a_X(1) ) * ( a_X(3) - a_X(2) ) ) &
          & * a_Y(3)
      end do
    else
      do k = 1, NLev
        iTemp = a_iTemp(k)

!!$        a_PF(k) = &
!!$          &   ( tb_PFTblPF(iTemp,iBand) - tb_PFTblPF(iTemp-1,iBand) ) &
!!$          & / ( t_PFTblTemp(iTemp)      - t_PFTblTemp(iTemp-1)      ) &
!!$          & * ( a_Temp(k)               - t_PFTblTemp(iTemp-1)      ) &
!!$          & + tb_PFTblPF(iTemp-1,iBand)

        XVal   = a_Temp(k)
        a_X(1) = t_PFTblTemp(iTemp-1)
        a_X(2) = t_PFTblTemp(iTemp  )
        a_X(3) = t_PFTblTemp(iTemp+1)
        a_Y(1) = tb_PFTblPF(iTemp-1,iBand)
        a_Y(2) = tb_PFTblPF(iTemp  ,iBand)
        a_Y(3) = tb_PFTblPF(iTemp+1,iBand)
        a_PF(k) = &
          &   ( ( XVal   - a_X(2) ) * ( XVal   - a_X(3) ) ) &
          & / ( ( a_X(1) - a_X(2) ) * ( a_X(1) - a_X(3) ) ) &
          & * a_Y(1)                                        &
          & + ( ( XVal   - a_X(3) ) * ( XVal   - a_X(1) ) ) &
          & / ( ( a_X(2) - a_X(3) ) * ( a_X(2) - a_X(1) ) ) &
          & * a_Y(2)                                        &
          & + ( ( XVal   - a_X(1) ) * ( XVal   - a_X(2) ) ) &
          & / ( ( a_X(3) - a_X(1) ) * ( a_X(3) - a_X(2) ) ) &
          & * a_Y(3)
      end do
    end if


  end subroutine OptPropKDPFTblGetPFDPFDT1D

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

  subroutine OptPropKDPFTblGetPFDPFDT0D( &
    & iBand, Temp, &
    & PF, &
    & FlagDPFDT &
    & )

    integer , intent(in ) :: iBand
    real(dp), intent(in ) :: Temp
    real(dp), intent(out) :: PF
    logical , intent(in ), optional :: FlagDPFDT


    ! Local variables
    !
    integer , parameter :: NLev = 1
    real(dp)            :: a_Temp(NLev)
    real(dp)            :: a_PF  (NLev)


    a_Temp = Temp

    call OptPropKDPFTblGetPFDPFDT1D( &
      & iBand, NLev, a_Temp, &
      & a_PF, &
      & FlagDPFDT &
      & )

    PF = a_PF(1)


  end subroutine OptPropKDPFTblGetPFDPFDT0D

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

  subroutine OptPropKDPFTblChkBandIntPFRatio( &
    & iBand, ThresholdBandIntPFRatio, &
    & FlagBandIntPFRatio &
    & )

    integer , intent(in ) :: iBand
    real(dp), intent(in ) :: ThresholdBandIntPFRatio
    logical , intent(out) :: FlagBandIntPFRatio


    ! Local variables
    !
    integer :: iTemp


    iTemp = PFTblNTemp
!!$    if ( tb_PFTblPFRatio(iTemp,iBand) >= ThresholdBandIntPFRatio ) then
    if ( tb_PFTblPFRatioCum(iTemp,iBand) >= ThresholdBandIntPFRatio ) then
      FlagBandIntPFRatio = .true.
    else
      FlagBandIntPFRatio = .false.
    end if


  end subroutine OptPropKDPFTblChkBandIntPFRatio

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

  subroutine OptPropKDPFTblChkBandIntStrFluxRatio( &
    & iBand, ThresholdBandIntStrFluxRatio, &
    & FlagBandIntStrFluxRatio &
    & )

    integer , intent(in ) :: iBand
    real(dp), intent(in ) :: ThresholdBandIntStrFluxRatio
    logical , intent(out) :: FlagBandIntStrFluxRatio


    ! Local variables
    !


!!$    if ( b_BandIntFluxtoTotFluxRatio(iBand) >= ThresholdBandIntStrFluxRatio ) then
    if ( b_BandIntFluxtoTotFluxRatioCum(iBand) >= ThresholdBandIntStrFluxRatio ) then
      FlagBandIntStrFluxRatio = .true.
    else
      FlagBandIntStrFluxRatio = .true.
    end if


  end subroutine OptPropKDPFTblChkBandIntStrFluxRatio

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

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


    character(*), intent(in ) :: AxisName
    integer     , intent(in ) :: NTblGrid
    real(DP)    , intent(in ) :: a_TblGrid (1:NTblGrid)
    integer     , intent(in ) :: NIntpol
    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


    ! check
    if ( a_TblGrid(1) > a_TblGrid(2) ) then
      write( 6, * ) 'Unexpected order of axis ', trim( AxisName ), &
        & ', i.e., descending order.'
      stop
    end if

    if ( NTblGrid-NIntpol+1 <= 0 ) then
      write( 6, * ) &
        & 'ERROR: Number of points in table axis is insufficient: ', &
        & NTblGrid, NIntpol
    end if

    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
      a_Index(i) = a_Index(i) - 1
      a_Index(i) = min( a_Index(i), NTblGrid-NIntpol+1 )
    end do


  end subroutine FindTblIndex

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

end module opt_prop_kd

