program main

  use vtype_module
  use fi_module
  use ni3_module
  use netcdf

  implicit none

  character(extstr)              :: IDChar

  integer                        :: InNcID
  integer                        :: OutNcID
  character(stdstr)              :: Mode
  integer                        :: NDims
  character(stdstr), allocatable :: a_DimNames( : )
  character(stdstr)              :: Name
  character(stdstr)              :: StdName
  character(stdstr)              :: LongName
  character(stdstr)              :: Units
  character(extstr)              :: Comment
  character(extstr*4)            :: CommentStrSpeGas
  character(extstr*4)            :: CommentStrSpePtcl
  real(DP)                       :: CommentDWNGas
  real(DP)                       :: CommentDWNPtcl
  real(DP)                       :: PressForPFSort

  integer, parameter     :: NBnd = 2

  integer                :: NPress
  integer                :: NMolNum
  integer                :: NBand
  integer                :: NBin
  integer                :: NPressForPFSort

  integer                :: NBandTmp

  real(8), allocatable   :: p_Press(:)
  integer, allocatable   :: m_MolNum(:)
  integer, allocatable   :: c_BinNum(:)
  real(8), allocatable   :: c_BinWaveNum(:)  ! This can be neglected, use c_BinNum
  integer                :: a_BndNum(2)
  real(8), allocatable   :: t_Temp(:)
  real(8), allocatable   :: m_VMR(:)
  real(8), allocatable   :: u_VMRMol1  (:)
  real(8), allocatable   :: v_VMRMol2  (:)
  real(8), allocatable   :: w_VMRMol3  (:)

  integer, allocatable   :: b_BandNum(:)

  real(8), allocatable   :: ab_BandWaveNumBnds(:,:)
  integer, allocatable   :: b_BandNBin(:)
  integer, allocatable   :: ab_BandBinIndexBnds(:,:)
!!$  integer, allocatable   :: b_BandMolNum(:)
  real(8), allocatable   :: c_BinWeight(:)

!!$  real(DP), allocatable :: ct_AtmPFRatio(:,:)
  real(DP), allocatable :: pctuvw_AtmPFRatio(:,:,:,:,:,:)
  real(DP), allocatable :: cp_StrPFRatio(:,:)
  real(DP), allocatable :: bp_BandAveStrPF(:,:)
  real(DP), allocatable :: bp_RatioBandIntFluxtoTotFlux(:,:)
  real(DP), allocatable :: pmctuvw_AC(:,:,:,:,:,:,:)
  real(DP), allocatable :: b_BandAveRayScatCoefNonRadAct(:)
  real(DP), allocatable :: bm_BandAveRayScatCoef(:,:)

  integer                        :: NChar
  integer          , allocatable :: a_CharNum(:)
  integer                        :: NParticle
  integer                        :: ParticleNRadiusMax
  integer          , allocatable :: a_ParticleNum      (:)
  integer          , allocatable :: r_ParticleRadiusNum(:)
  character(stdstr), allocatable :: a_ParticleName   (:)
  integer          , allocatable :: a_ParticleNRadius(:)
  real(DP)         , allocatable :: ra_ParticleRadius(:,:)
  real(DP)         , allocatable :: bra_BandAveQExt  (:,:,:)
  real(DP)         , allocatable :: bra_BandAveSSA   (:,:,:)
  real(DP)         , allocatable :: bra_BandAveAF    (:,:,:)

  integer, allocatable   :: mc_IDAbsCoefType(:,:)

  real(8), allocatable   :: p_PressTmp(:)
  real(8), allocatable   :: p_TempTmp(:)
  real(8), allocatable   :: pm_VMRTmp(:,:)
  real(8), allocatable   :: pc_AtmPFRatioTmp(:,:)
  real(8), allocatable   :: pc_StrPFRatioTmp(:,:)
  real(DP), allocatable   :: pb_RatioBandIntFluxtoTotFluxTmp(:,:)
  real(8), allocatable   :: pmc_AbsCoef(:,:,:)

  integer              :: NMolNumTmp
  integer, allocatable :: m_MolNumTmp(:)

  character(128) :: InDirName
  character(128) :: InPtclOptPropNcFN
  character(128) :: InNcFN
  character(128) :: OutNcFN

  logical :: FlagChkConstVMR

  character(128) :: CtlFN = 'mkkdtbl.conf'
  character(extstr) :: CtlFN3 = 'many.conf'

  integer :: FileUnit
  integer :: iostat

  integer :: k
  integer :: l
  integer :: m
  integer :: Mol1Num
  integer :: Mol2Num
  integer :: Mol3Num
  integer :: iVMRMol1
  integer :: iVMRMol1s
  integer :: iVMRMol1e
  integer :: iVMRMol2
  integer :: iVMRMol2s
  integer :: iVMRMol2e
  integer :: iVMRMol3
  integer :: iVMRMol3s
  integer :: iVMRMol3e
  integer :: iTemp
  integer :: iTemps
  integer :: iTempe

  real(DP)                    :: PFTblDelTemp
  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 (:,:)

!!$  integer :: kIndexForPF


  namelist /id_nml/ &
    & IDChar

  namelist /file_nml/ &
    & InDirName, InPtclOptPropNcFN, OutNcFN

  namelist /flag_nml/ &
    & FlagChkConstVMR

  namelist /many_nml/ &
    & iTemps, iTempe, &
    & Mol1Num, iVMRMol1s, iVMRMol1e, &
    & Mol2Num, iVMRMol2s, iVMRMol2e, &
    & Mol3Num, iVMRMol3s, iVMRMol3e

  namelist /pftbl_nml/ &
    & PFTblDelTemp


  IDChar            = ""

  FlagChkConstVMR   = .true.

  InDirName         = "indir"
  InPtclOptPropNcFN = "rsc.nc"
  OutNcFN           = "out/kdactbl.nc"

  PFTblDelTemp      = 1.0d0


  call fi_open( CtlFN, "read", FileUnit )
  !
  rewind( FileUnit )
  read( FileUnit, nml = file_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist, file_nml, is something wrong.'
  end if
  write( 6, nml = file_nml )
  !
  rewind( FileUnit )
  read( FileUnit, nml = id_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist, id_nml, is something wrong.'
  end if
  write( 6, nml = id_nml )
  !
  rewind( FileUnit )
  read( FileUnit, nml = flag_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist, flag_nml, is something wrong.'
  end if
  write( 6, nml = flag_nml )
  !
  rewind( FileUnit )
  read( FileUnit, nml = pftbl_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist, pftbl_nml, is something wrong.'
  end if
  write( 6, nml = pftbl_nml )
  !
  close( FileUnit )

  call fi_open( CtlFN3, "read", FileUnit )
  !
  rewind( FileUnit )
  read( FileUnit, nml = many_nml, iostat = iostat )
  if ( iostat > 0 ) then
    stop 'Namelist in many_nml is something wrong.'
  end if
  write( 6, nml = many_nml )
  !
  close( FileUnit )


!!$  NVMRH2O = iVMRH2Oe - iVMRH2Os + 1


  do iVMRMol3 = iVMRMol3s, iVMRMol3e
    do iVMRMol2 = iVMRMol2s, iVMRMol2e
      do iVMRMol1 = iVMRMol1s, iVMRMol1e

        do iTemp = iTemps, iTempe

          write( InNcFN, '(a,i3.3,3(a,i2.2,a,i3.3),a)' ) &
            & trim( InDirName ) // "/" &
            & // 'ac_T-', iTemp, &
            & '_M', Mol1Num, '-', iVMRMol1, &
            & '_M', Mol2Num, '-', iVMRMol2, &
            & '_M', Mol3Num, '-', iVMRMol3, &
            & '-00-10.nc'

          Mode = "read"
          call ni3_open( InNcFN, Mode, InNcID )
          !
          call ni3_inq_dimlen( InNcID, "Press" , NPress  )
          call ni3_inq_dimlen( InNcID, "MolNum", NMolNum )

          call ni3_inq_dimlen( InNcID, "BandNum", NBand )
          call ni3_inq_dimlen( InNcID, "BinNum" , NBin  )


          if ( ( iVMRMol1 == iVMRMol1s ) .and. &
            &  ( iVMRMol2 == iVMRMol2s ) .and. &
            &  ( iVMRMol3 == iVMRMol3s ) .and. &
            &  ( iTemp == iTemps ) ) then

            allocate( p_Press(NPress) )
            allocate( m_MolNum(NMolNum) )

            allocate( t_Temp(iTemps:iTempe) )
            allocate( m_VMR(NMolNum) )
            allocate( u_VMRMol1(iVMRMol1s:iVMRMol1e) )
            allocate( v_VMRMol2(iVMRMol2s:iVMRMol2e) )
            allocate( w_VMRMol3(iVMRMol3s:iVMRMol3e) )


            allocate( b_BandNum(NBand) )
            allocate( c_BinNum    (NBin) )
            allocate( c_BinWaveNum(NBin) )
            allocate( ab_BandWaveNumBnds(NBnd,NBand) )
            allocate( b_BandNBin(NBand) )
            allocate( ab_BandBinIndexBnds(NBnd,NBand) )
!!$          allocate( b_BandMolNum(NBand) )
            allocate( c_BinWeight(NBin) )

!!$            allocate( ct_AtmPFRatio(NBin,iTemps:iTempe) )
            allocate( pctuvw_AtmPFRatio &
              & (NPress,NBin, &
              &  iTemps:iTempe, &
              &  iVMRMol1s:iVMRMol1e, &
              &  iVMRMol2s:iVMRMol2e, &
              &  iVMRMol3s:iVMRMol3e) )
            allocate( cp_StrPFRatio(NBin,NPress) )
            allocate( bp_BandAveStrPF(NBand,NPress) )
            allocate( bp_RatioBandIntFluxtoTotFlux(NBand,NPress) )
            allocate( pmctuvw_AC &
              & (NPress,NMolNum,NBin, &
              &  iTemps:iTempe, &
              &  iVMRMol1s:iVMRMol1e, &
              &  iVMRMol2s:iVMRMol2e, &
              &  iVMRMol3s:iVMRMol3e) )
            allocate( mc_IDAbsCoefType(NMolNum,NBin) )

            allocate( p_TempTmp(NPress) )
            allocate( pm_VMRTmp(NPress,NMolNum) )
            allocate( pc_AtmPFRatioTmp(NPress,NBin) )
            allocate( pc_StrPFRatioTmp(NPress,NBin) )
            allocate( pb_RatioBandIntFluxtoTotFluxTmp(NPress,NBand) )
          end if


          if ( ( iVMRMol1 == iVMRMol1s ) .and. &
            &  ( iVMRMol2 == iVMRMol2s ) .and. &
            &  ( iVMRMol3 == iVMRMol3s ) .and. &
            &  ( iTemp == iTemps ) ) then

            call ni3_get_att( InNcID, '_global_', 'comment_dwn_m-1', CommentDWNGas )
            call ni3_get_att( InNcID, '_global_', 'comment_strspe', CommentStrSpeGas )
!!$            call ni3_get_att( InNcID, '_global_', 'pressure_for_P.F._sort_Pa', PressForPFSort )

            call ni3_get_var( InNcID, "BndNum", a_BndNum )
            call ni3_get_var( InNcID, "Press" , p_Press  )

            allocate( p_PressTmp(NPress) )
            p_PressTmp = p_Press
            do k = 1, NPress
              p_Press(k) = p_PressTmp(NPress-(k-1))
            end do
            deallocate( p_PressTmp )

            call ni3_get_var( InNcID, "MolNum"          , m_MolNum            )

            call ni3_get_var( InNcID, "BandNum"         , b_BandNum           )
            call ni3_get_var( InNcID, "BinNum"          , c_BinNum            )
            call ni3_get_var( InNcID, "BinWaveNum"      , c_BinWaveNum        )
            call ni3_get_var( InNcID, "BandWaveNumBnds" , ab_BandWaveNumBnds  )
            call ni3_get_var( InNcID, "BandNBin"        , b_BandNBin          )
            call ni3_get_var( InNcID, "BandBinIndexBnds", ab_BandBinIndexBnds )
!!$          call ni3_get_var( InNcID, "BandMolNum"   , b_BandMolNum     )
            call ni3_get_var( InNcID, "BinWeight"       , c_BinWeight         )

            call ni3_get_var( InNcID, "VMR", pm_VMRTmp )
            do m = 1, NMolNum
              m_VMR(m) = pm_VMRTmp(1,m)
            end do

          end if

          call ni3_get_var( InNcID, "Temp", p_TempTmp )
          do k = 1+1, NPress
            if ( p_TempTmp(k) /= p_TempTmp(1) ) then
              write( 6, * ) "Error: Temperature is not uniform in vertical direction."
              write( 6, * ) p_TempTmp
              stop
            end if
          end do
          if ( ( iVMRMol1 == iVMRMol1s ) .and. &
            &  ( iVMRMol2 == iVMRMol2s ) .and. &
            &  ( iVMRMol3 == iVMRMol3s ) ) then
            t_Temp(iTemp) = p_TempTmp(1)
          end if

          call ni3_get_var( InNcID, "VMR", pm_VMRTmp )
          if ( FlagChkConstVMR ) then
            do m = 1, NMolNum
              do k = 1, NPress
                if ( pm_VMRTmp(1,m) /= pm_VMRTmp(k,m) ) then
                  write( 6, * ) "Error: Mixing ratio is not uniform in vertical direction."
                  write( 6, * ) "  ", pm_VMRTmp
                  stop
                end if
              end do
            end do
          end if

          ! check mixing ratio of non-key molecule
          do m = 1, NMolNum
            if ( ( m_MolNum(m) == Mol1Num ) .or. &
              &  ( m_MolNum(m) == Mol2Num ) .or. &
              &  ( m_MolNum(m) == Mol3Num ) ) &
              & cycle
            if ( m_VMR(m) /= pm_VMRTmp(1,m) ) then
              write( 6, * ) "Error: Unexpected mixing ratio: ", m_VMR(m), pm_VMRTmp(1,m)
              stop
            end if
          end do
          !
          if ( iTemp == iTemps ) then
            ! set volume mixing ratio of key molecule 1
            if ( ( iVMRMol2 == iVMRMol2s ) .and. &
              &  ( iVMRMol3 == iVMRMol3s ) ) then
              do m = 1, NMolNum
                if ( m_MolNum(m) == Mol1Num ) exit
              end do
              if ( m > NMolNum ) then
                write( 6, * ) "Error: Unexpected molecular number for key molecule 1"
                stop
              end if
              u_VMRMol1(iVMRMol1) = pm_VMRTmp(1,m)
            end if
            !
            ! set volume mixing ratio of key molecule 2
            if ( ( iVMRMol1 == iVMRMol1s ) .and. &
              &  ( iVMRMol3 == iVMRMol3s ) ) then
              do m = 1, NMolNum
                if ( m_MolNum(m) == Mol2Num ) exit
              end do
              if ( m > NMolNum ) then
                write( 6, * ) "Error: Unexpected molecular number for key molecule 2"
                stop
              end if
              v_VMRMol2(iVMRMol2) = pm_VMRTmp(1,m)
            end if
            !
            ! set volume mixing ratio of key molecule 3
            if ( ( iVMRMol1 == iVMRMol1s ) .and. &
              &  ( iVMRMol2 == iVMRMol2s ) ) then
              do m = 1, NMolNum
                if ( m_MolNum(m) == Mol3Num ) exit
              end do
              if ( m > NMolNum ) then
                write( 6, * ) "Error: Unexpected molecular number for key molecule 3"
                stop
              end if
              w_VMRMol3(iVMRMol3) = pm_VMRTmp(1,m)
            end if
          end if

          call ni3_inq_dimlen( InNcID, "PressForPF" , NPressForPFSort  )
          if ( NPressForPFSort /= NPress ) then
            write( 6, * ) 'NPressForPF does not match with NPress.'
            write( 6, * ) NPressForPFSort, NPress
            stop
          end if

          call ni3_get_var( InNcID, "BinAtmPFRatio", pc_AtmPFRatioTmp )
!!$          kIndexForPF = 1
!!$          do k = 1, NPress
!!$            if ( abs( p_Press(k          ) - PressForPFSort ) <        &
!!$              &  abs( p_Press(kIndexForPF) - PressForPFSort )   ) then
!!$              kIndexForPF = k
!!$            end if
!!$          end do
!!$          write( 6, * ) &
!!$            & 'Pressure for Planck function sorting : ', &
!!$            & p_Press(kIndexForPF), ' Pa.'
!!$          write( 6, * ) &
!!$            & '  Specified pressure                 : ', &
!!$            & PressForPFSort, ' Pa.'
!!$          do l = 1, NBin
!!$            ! It should be understood that order of p_Press has been reversed
!!$            ! above.
!!$            k = kIndexForPF
!!$            ct_AtmPFRatio(l,iTemp) = pc_AtmPFRatioTmp(NPress-(k-1),l)
!!$          end do
          do l = 1, NBin
            do m = 1, NMolNum
              do k = 1, NPress
                pctuvw_AtmPFRatio(k,l,iTemp,iVMRMol1,iVMRMol2,iVMRMol3) = &
                  & pc_AtmPFRatioTmp(NPress-(k-1),l)
              end do
            end do
          end do

          if ( ( iVMRMol1 == iVMRMol1s ) .and. &
            &  ( iVMRMol2 == iVMRMol2s ) .and. &
            &  ( iVMRMol3 == iVMRMol3s ) .and. &
            &  ( iTemp == iTemps ) ) then
            call ni3_get_var( InNcID, "BinStrPFRatio", pc_StrPFRatioTmp )
            do l = 1, NBin
              do k = 1, NPress
!!$                cp_StrPFRatio(l,k) = pc_StrPFRatioTmp(k,l)
                cp_StrPFRatio(l,k) = pc_StrPFRatioTmp(NPress-(k-1),l)
              end do
            end do
            !
            call ni3_get_var( InNcID, "RatioBandIntFluxtoTotFlux", pb_RatioBandIntFluxtoTotFluxTmp )
            do l = 1, NBand
              do k = 1, NPress
!!$                bp_RatioBandIntFluxtoTotFlux(l,k) = pb_RatioBandIntFluxtoTotFluxTmp(k,l)
                bp_RatioBandIntFluxtoTotFlux(l,k) = pb_RatioBandIntFluxtoTotFluxTmp(NPress-(k-1),l)
              end do
            end do
          end if

!!$      call HistoryGet( InNcFN, "rma_TAbsorpK", pmctv_AC(:,:,:,iTemp,iVMRH2O), quiet = .true. )
          allocate( pmc_AbsCoef(NPress,NMolNum,NBin) )
          call ni3_get_var( InNcID, "BinAbsCoef", pmc_AbsCoef )
          do l = 1, NBin
            do m = 1, NMolNum
              do k = 1, NPress
                pmctuvw_AC(k,m,l,iTemp,iVMRMol1,iVMRMol2,iVMRMol3) = &
                  & pmc_AbsCoef(NPress-(k-1),m,l)
              end do
            end do
          end do
          deallocate( pmc_AbsCoef )

          call ni3_get_var( InNcID, "BinIDAbsCoefType", mc_IDAbsCoefType )

          call ni3_close( InNcID )

        end do
      end do
    end do
  end do

  ! Input Rayleigh scattering coefficient
  allocate( b_BandAveRayScatCoefNonRadAct(NBand) )
  allocate( bm_BandAveRayScatCoef(NBand,NMolNum) )
  !
  Mode = "read"
  InNcFN = InPtclOptPropNcFN
  call ni3_open( InNcFN, Mode, InNcID )
  !
  call ni3_inq_dimlen( InNcID, 'BandNum', NBandTmp )
  if ( NBand /= NBandTmp ) then
    write( 6, * ) 'Length of BandNum dimension in a file for particle is different from that in a file for gas.'
    write( 6, * ) NBandTmp, NBand
    stop
  end if
  !
  call ni3_get_att( InNcID, '_global_', 'comment_dwn_m-1', CommentDWNPtcl )
  call ni3_get_att( InNcID, '_global_', 'comment_strspe', CommentStrSpePtcl )
  !
  call ni3_inq_dimlen( InNcID, "MolNum", NMolNumTmp )
  if ( NMolNumTmp /= NMolNum ) then
    write( 6, * ) 'Length of MolNum dimension in a file for particle is different from that in a file for gas.'
    write( 6, * ) NMolNumTmp, NMolNum
    stop
  end if
  allocate( m_MolNumTmp(1:NMolNum) )
  call ni3_get_var( InNcID, "MolNum", m_MolNumTmp )
  do m = 1, NMolNum
    if ( m_MolNumTmp(m) /= m_MolNum(m) ) then
      write( 6, * ) 'Molecular number ', m, ' is different between files for gas and particles.'
      write( 6, * ) '  ', trim( InNcFN ), ': ', m_MolNumTmp
      write( 6, * ) '  ', trim( 'Gas ' ), ': ', m_MolNum
      stop
    end if
  end do
  deallocate( m_MolNumTmp )
  !
  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"      , NParticle          )
  call ni3_inq_dimlen( InNcID, "PtclRadiusNum", ParticleNRadiusMax )
  allocate( a_CharNum          (NChar             ) )
  allocate( a_ParticleNum      (NParticle         ) )
  allocate( r_ParticleRadiusNum(ParticleNRadiusMax) )
  allocate( a_ParticleName     (NParticle         ) )
  allocate( a_ParticleNRadius  (NParticle         ) )
  allocate( ra_ParticleRadius(ParticleNRadiusMax,NParticle) )
  call ni3_get_var( InNcID, "CharNum"      , a_CharNum           )
  call ni3_get_var( InNcID, "PtclNum"      , a_ParticleNum       )
  call ni3_get_var( InNcID, "PtclRadiusNum", r_ParticleRadiusNum )
  call ni3_get_var( InNcID, "PtclName"     , a_ParticleName      )
  call ni3_get_var( InNcID, "PtclNRadius"  , a_ParticleNRadius   )
  call ni3_get_var( InNcID, "PtclRadius"   , ra_ParticleRadius   )
  allocate( bra_BandAveQExt  (NBand,ParticleNRadiusMax,NParticle) )
  allocate( bra_BandAveSSA   (NBand,ParticleNRadiusMax,NParticle) )
  allocate( bra_BandAveAF    (NBand,ParticleNRadiusMax,NParticle) )
  call ni3_get_var( InNcID, "BandAveQExt"      , bra_BandAveQExt     )
  call ni3_get_var( InNcID, "BandAveSSA"       , bra_BandAveSSA      )
  call ni3_get_var( InNcID, "BandAveAF"        , bra_BandAveAF       )



  ! Table of Planck function
  !
  call PFTblInit


  ! Output

  Mode = "new"
  call ni3_open( OutNcFN, Mode, OutNcID )

  call ni3_put_att( OutNcID, '_global_', 'id', IDChar )

  write( 6, * ) 'Wavenumber resolution'
  write( 6, * ) '  gas abs. coef.     : ', CommentDWNGas
  write( 6, * ) '  particle opt. prop.: ', CommentDWNPtcl
  call ni3_put_att( OutNcID, '_global_', 'comment_dwn_m-1_gas' , CommentDWNGas  )
  call ni3_put_att( OutNcID, '_global_', 'comment_dwn_m-1_ptcl', CommentDWNPtcl )

  write( 6, * ) 'Stellar flux'
  write( 6, * ) '  gas abs. coef.     : ', trim( CommentStrSpeGas  )
  write( 6, * ) '  particle opt. prop.: ', trim( CommentStrSpePtcl )
  call ni3_put_att( OutNcID, '_global_', 'comment_strspe_gas' , CommentStrSpeGas  )
  call ni3_put_att( OutNcID, '_global_', 'comment_strspe_ptcl', CommentStrSpePtcl )

  call ni3_put_att( OutNcID, '_global_', 'pressure_for_P.F._sort', PressForPFSort )

  Name     = "Press"
  StdName  = "plev"
  LongName = "pressure"
  Units    = "Pa"
  call ni3_set_dim( OutNcID, Name, NF90_DOUBLE, p_Press  , &
    & stdname = StdName, longname = LongName, units  = Units )
  !
  Name     = "MolNum"
  StdName  = "MolNum"
  LongName = "molecular number"
  Units    = "1"
  call ni3_set_dim( OutNcID, Name, NF90_INT   , m_MolNum, &
    & stdname = StdName, longname = LongName, units  = Units )
  !
  Name     = "BandNum"
  StdName  = "BandNum"
  LongName = "band number"
  Units    = "1"
  call ni3_set_dim( OutNcID, Name, NF90_INT   , b_BandNum, &
    & stdname = StdName, longname = LongName, units  = Units )
  !
  Name     = "BinNum"
  StdName  = "BinNum"
  LongName = "bin number"
  Units    = "1"
  call ni3_set_dim( OutNcID, Name, NF90_INT   , c_BinNum, &
    & stdname = StdName, longname = LongName, units  = Units )
  !
  Name     = "BinWaveNum"
  StdName  = "BinWaveNum"
  LongName = "bin wavenumber"
  Units    = "m-1"
  call ni3_set_dim( OutNcID, Name, NF90_DOUBLE, c_BinWaveNum, &
    & stdname = StdName, longname = LongName, units  = Units )
  Comment = "This dimension is present for visualization purpose. This can be replaced with BinNum dimension."
  call ni3_put_att( OutNcID, Name, "comment", Comment )
  !
  Name     = "BndNum"
  StdName  = "BndNum"
  LongName = "bound number"
  Units    = "1"
  call ni3_set_dim( OutNcID, Name, NF90_INT   , a_BndNum, &
    & stdname = StdName, longname = LongName, units  = Units )
  !
  Name     = "Temp"
  StdName  = Name
  LongName = "temperature"
  Units    = "K"
  call ni3_set_dim( OutNcID, Name, NF90_DOUBLE, t_Temp  , &
    & stdname = StdName, longname = LongName, units  = Units )
  !
  Name     = "VMRMol1"
  StdName  = Name
  write( LongName, '(a,i2.2)' ) &
    & "volume mixing ratio of molecular 1, number = ", Mol1Num
  Units    = "1"
  call ni3_set_dim( OutNcID, Name, NF90_DOUBLE, u_VMRMol1, &
    & stdname = StdName, longname = LongName, units  = Units )
  !
  Name     = "VMRMol2"
  StdName  = Name
  write( LongName, '(a,i2.2)' ) &
    & "volume mixing ratio of molecular 2, number = ", Mol2Num
  Units    = "1"
  call ni3_set_dim( OutNcID, Name, NF90_DOUBLE, v_VMRMol2, &
    & stdname = StdName, longname = LongName, units  = Units )
  !
  Name     = "VMRMol3"
  StdName  = Name
  write( LongName, '(a,i2.2)' ) &
    & "volume mixing ratio of molecular 3, number = ", Mol3Num
  Units    = "1"
  call ni3_set_dim( OutNcID, Name, NF90_DOUBLE, w_VMRMol3, &
    & stdname = StdName, longname = LongName, units  = Units )
  !
  Name     = "PFTblTemp"
  StdName  = Name
  LongName = "temperature for Planck function table"
  Units    = "K"
  call ni3_set_dim( OutNcID, Name, NF90_DOUBLE, t_PFTblTemp, &
    & stdname = StdName, longname = LongName, units  = Units )
  !
  !
  !
  NDims = 0
  allocate( a_DimNames( NDims ) )
  Name     = "VMRMol1Num"
  StdName  = Name
  LongName = "molecular number of volume mixing ratio dimension 1"
  Units    = "1"
  call ni3_def_var( OutNcID, Name, NF90_INT, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, Mol1Num )
  deallocate( a_DimNames )
  !
  NDims = 0
  allocate( a_DimNames( NDims ) )
  Name     = "VMRMol2Num"
  StdName  = Name
  LongName = "molecular number of volume mixing ratio dimension 2"
  Units    = "1"
  call ni3_def_var( OutNcID, Name, NF90_INT, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, Mol2Num )
  deallocate( a_DimNames )
  !
  NDims = 0
  allocate( a_DimNames( NDims ) )
  Name     = "VMRMol3Num"
  StdName  = Name
  LongName = "molecular number of volume mixing ratio dimension 3"
  Units    = "1"
  call ni3_def_var( OutNcID, Name, NF90_INT, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, Mol3Num )
  deallocate( a_DimNames )
  !
  NDims = 1
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "MolNum"
  Name     = "VMR"
  StdName  = Name
  LongName = "volume mixing ratio"
  Units    = "1"
  call ni3_def_var( OutNcID, Name, NF90_DOUBLE, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, m_VMR )
  deallocate( a_DimNames )
  !
  NDims = 2
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "BndNum"
  a_DimNames( 2 ) = "BandNum"
  Name     = "BandWaveNumBnds"
  StdName  = Name
  LongName = "band wavenumber bound"
  Units    = "m-1"
  call ni3_def_var( OutNcID, Name, NF90_DOUBLE, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, ab_BandWaveNumBnds )
  deallocate( a_DimNames )

  NDims = 1
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "BandNum"
  Name     = "BandNBin"
  StdName  = Name
  LongName = "number of bins in a band"
  Units    = "1"
  call ni3_def_var( OutNcID, Name, NF90_INT, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, b_BandNBin )
  deallocate( a_DimNames )


  NDims = 2
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "BndNum"
  a_DimNames( 2 ) = "BandNum"
  Name     = "BandBinIndexBnds"
  StdName  = Name
  LongName = "index bound of bins in a band"
  Units    = "1"
  call ni3_def_var( OutNcID, Name, NF90_INT, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, ab_BandBinIndexBnds )
  deallocate( a_DimNames )


!!$  NDims = 1
!!$  allocate( a_DimNames( NDims ) )
!!$  a_DimNames( 1 ) = "BandNum"
!!$  Name     = "BandMolNum"
!!$  StdName  = Name
!!$  LongName = "band key molecular number"
!!$  Units    = "1"
!!$  call ni3_def_var( OutNcID, Name, NF90_INT, NDims, a_DimNames, &
!!$    & stdname = StdName, longname = LongName, units = Units )
!!$  call ni3_put_var( OutNcID, Name, b_BandMolNum )
!!$  deallocate( a_DimNames )


  NDims = 1
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "BinNum"
  Name     = "BinWeight"
  StdName  = Name
  LongName = "bin weight"
  Units    = "m-1"
  call ni3_def_var( OutNcID, Name, NF90_DOUBLE, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, c_BinWeight )
  deallocate( a_DimNames )

  NDims = 7
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "Press"
  a_DimNames( 2 ) = "MolNum"
  a_DimNames( 3 ) = "BinWaveNum"
  a_DimNames( 4 ) = "Temp"
  a_DimNames( 5 ) = "VMRMol1"
  a_DimNames( 6 ) = "VMRMol2"
  a_DimNames( 7 ) = "VMRMol3"
  Name     = "BinAbsCoef"
  StdName  = Name
  LongName = "absorption_coefficient"
!!$!          Units    = "(kg m-2)-1"
  Units    = "(molecule m-2)-1"
  call ni3_def_var( OutNcID, name, NF90_DOUBLE, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, pmctuvw_AC )
  deallocate( a_DimNames )
  !
  NDims = 2
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "MolNum"
  a_DimNames( 2 ) = "BinWaveNum"
  Name     = "BinIDAbsCoefType"
  StdName  = Name
  LongName = "absorption type ID"
  Units    = "1"
  call ni3_def_var( OutNcID, Name, NF90_INT, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, mc_IDAbsCoefType )
  Comment = "-1: zero abs. coef., 0: abs. coef. per unit abs. molecule, 1: abs. coef. per unit atm. molecule"
  call ni3_put_att( OutNcID, Name, "comment", Comment )
  deallocate( a_DimNames )
  !
!!$  NDims = 2
!!$  allocate( a_DimNames( NDims ) )
!!$  a_DimNames( 1 ) = "BinWaveNum"
!!$  a_DimNames( 2 ) = "Temp"
!!$  Name     = "BinAtmPFRatio"
!!$  StdName  = Name
!!$  LongName = "Planck function ratio"
!!$  Units    = "1"
!!$  call ni3_def_var( OutNcID, Name, NF90_DOUBLE, NDims, a_DimNames, &
!!$    & stdname = StdName, longname = LongName, units = Units )
!!$  call ni3_put_var( OutNcID, Name, ct_AtmPFRatio )
!!$  deallocate( a_DimNames )
  !
  NDims = 6
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "Press"
  a_DimNames( 2 ) = "BinWaveNum"
  a_DimNames( 3 ) = "Temp"
  a_DimNames( 4 ) = "VMRMol1"
  a_DimNames( 5 ) = "VMRMol2"
  a_DimNames( 6 ) = "VMRMol3"
  Name     = "BinAtmPFRatio"
  StdName  = Name
  LongName = "Planck function ratio"
  Units    = "1"
  call ni3_def_var( OutNcID, Name, NF90_DOUBLE, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, pctuvw_AtmPFRatio )
  deallocate( a_DimNames )
  !
  NDims = 2
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "BinWaveNum"
  a_DimNames( 2 ) = "Press"
  Name     = "BinStrPFRatio"
  StdName  = Name
  LongName = "Star Planck function Ratio"
  Units    = "1"
  call ni3_def_var( OutNcID, Name, NF90_DOUBLE, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, cp_StrPFRatio )
  deallocate( a_DimNames )
  !
!!$  NDims = 2
!!$  allocate( a_DimNames( NDims ) )
!!$  a_DimNames( 1 ) = "BandNum"
!!$  a_DimNames( 2 ) = "Press"
!!$  Name     = "BandAveStrPF"
!!$  StdName  = Name
!!$  LongName = "Band averaged stellar flux (Planck function)"
!!$  Units    = "1"
!!$  call ni3_def_var( OutNcID, Name, NF90_DOUBLE, NDims, a_DimNames, &
!!$    & stdname = StdName, longname = LongName, units = Units )
!!$  call ni3_put_var( OutNcID, Name, bp_BandAveStrPF )
!!$  deallocate( a_DimNames )
  !
  NDims = 2
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "BandNum"
  a_DimNames( 2 ) = "Press"
  Name     = "RatioBandIntFluxtoTotFlux"
  StdName  = Name
  LongName = "Ratio of band integrated flux to total flux (Planck function)"
  Units    = "1"
  call ni3_def_var( OutNcID, Name, NF90_DOUBLE, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, bp_RatioBandIntFluxtoTotFlux )
  deallocate( a_DimNames )


  Name     = "PtclNum"
  call ni3_set_dim( OutNcID, Name, NI3_INT, a_ParticleNum )
  call ni3_cp_atts( InNcID, OutNcID, Name, Name )
  !
  Name     = "PtclRadiusNum"
  call ni3_set_dim( OutNcID, Name, NI3_INT, r_ParticleRadiusNum )
  call ni3_cp_atts( InNcID, OutNcID, Name, Name )
  !
  name    = "CharNum"
  call ni3_set_dim( OutNcID, name, NI3_INT, a_CharNum )
  call ni3_cp_atts( InNcID, OutNcID, Name, Name )

  !
  NDims = 1
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "BandNum"
  Name     = "BandAveRayScatCoefNonRadAct"
  call ni3_def_var( OutNcID, Name, NF90_DOUBLE, NDims, a_DimNames )
  call ni3_cp_atts( InNcID, OutNcID, Name, Name )
  call ni3_put_var( OutNcID, Name, b_BandAveRayScatCoefNonRadAct )
  deallocate( a_DimNames )
  !
  NDims = 2
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "BandNum"
  a_DimNames( 2 ) = "MolNum"
  Name     = "BandAveRayScatCoef"
  call ni3_def_var( OutNcID, Name, NF90_DOUBLE, NDims, a_DimNames )
  call ni3_cp_atts( InNcID, OutNcID, Name, Name )
  call ni3_put_var( OutNcID, Name, bm_BandAveRayScatCoef )
  deallocate( a_DimNames )
  !
  NDims = 2
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "CharNum"
  a_DimNames( 2 ) = "PtclNum"
  name    = "PtclName"
  call ni3_def_var( OutNcID, name, NI3_CHAR, NDims, a_DimNames )
  call ni3_cp_atts( InNcID, OutNcID, Name, Name )
  call ni3_put_var( OutNcID, name, a_ParticleName )
  !
  deallocate( a_DimNames )
  !
  NDims = 1
  allocate( a_DimNames( NDims ) )
  !
  a_DimNames( 1 ) = "PtclNum"
  name    = "PtclNRadius"
  call ni3_def_var( OutNcID, name, NI3_DOUBLE, NDims, a_DimNames )
  call ni3_cp_atts( InNcID, OutNcID, Name, Name )
  call ni3_put_var( OutNcID, name, a_ParticleNRadius )
  !
  deallocate( a_DimNames )
  !
  NDims = 2
  allocate( a_DimNames( NDims ) )
  !
  a_DimNames( 1 ) = "PtclRadiusNum"
  a_DimNames( 2 ) = "PtclNum"
  name    = "PtclRadius"
  call ni3_def_var( OutNcID, name, NI3_DOUBLE, NDims, a_DimNames )
  call ni3_cp_atts( InNcID, OutNcID, Name, Name )
  call ni3_put_var( OutNcID, name, ra_ParticleRadius )
  !
  deallocate( a_DimNames )

  NDims = 3
  allocate( a_DimNames( NDims ) )
  !
  a_DimNames( 1 ) = "BandNum"
!!$  a_DimNames( 1 ) = "BandWaveNum"
  a_DimNames( 2 ) = "PtclRadiusNum"
  a_DimNames( 3 ) = "PtclNum"
  !
  Name    = "BandAveQExt"
  call ni3_def_var( OutNcID, Name, NI3_DOUBLE, NDims, a_DimNames )
  call ni3_cp_atts( InNcID, OutNcID, Name, Name )
  call ni3_put_var( OutNcID, Name, bra_BandAveQExt )
  !
  Name    = "BandAveSSA"
  call ni3_def_var( OutNcID, Name, NI3_DOUBLE, NDims, a_DimNames )
  call ni3_cp_atts( InNcID, OutNcID, Name, Name )
  call ni3_put_var( OutNcID, Name, bra_BandAveSSA  )
  !
  Name    = "BandAveAF"
  call ni3_def_var( OutNcID, Name, NI3_DOUBLE, NDims, a_DimNames )
  call ni3_cp_atts( InNcID, OutNcID, Name, Name )
  call ni3_put_var( OutNcID, Name, bra_BandAveAF   )
  !
  deallocate( a_DimNames )


  NDims = 2
  allocate( a_DimNames( NDims ) )
  a_DimNames( 1 ) = "PFTblTemp"
  a_DimNames( 2 ) = "BandNum"
  !
  Name     = "PFTblPF"
  StdName  = Name
  LongName = "table of Planck function"
  Units    = "W m-2"
  call ni3_def_var( OutNcID, Name, NF90_DOUBLE, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, tb_PFTblPF )
  !
  Name     = "PFTblDPFDT"
  StdName  = Name
  LongName = "table of derivative of Planck function"
  Units    = "W m-2 K-1"
  call ni3_def_var( OutNcID, Name, NF90_DOUBLE, NDims, a_DimNames, &
    & stdname = StdName, longname = LongName, units = Units )
  call ni3_put_var( OutNcID, Name, tb_PFTblDPFDT )
  !
  deallocate( a_DimNames )


  call ni3_close( OutNcID )

  call ni3_close( InNcID )


  ! Deallocate table of Planck function
  !
  call PFTblFinalize


  write( 6, * ) 'Finish'

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

contains

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

  subroutine PFTblInit

    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), allocatable :: a_PFTblChkTemp    (:)
    real(dp), allocatable :: a_PFTblChkVal     (:)
    real(dp), allocatable :: ab_PFTblChkPF     (:,:)
    real(dp), allocatable :: ab_PFTblChkDPFDT  (:,:)
    real(dp), allocatable :: a_PFTblChkPFSum   (:)
    real(dp), allocatable :: a_PFTblChkDPFDTSum(:)

    real(dp), parameter :: PFTblThreshold = 1.0d-5
    real(dp), parameter :: PFTblRelError  = 1.0d-3

    real(dp) :: WNs
    real(dp) :: WNe

    integer :: GQNum

    integer :: iBand
    integer :: iTemp

    integer :: NLev


    write( 6, * ) 'Preparing Planck function table ...'

    ! Last +1 is preparation for quadratic interpolation.
    PFTblNTemp = int( ( t_Temp(iTempe) - t_Temp(iTemps) ) / PFTblDelTemp ) + 1 + 1
    PFTblNBand = NBand

    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_Temp(iTemps) + PFTblDelTemp * (iTemp-1)
    end do

    do iBand = 1, PFTblNBand
      WNs = ab_BandWaveNumBnds(1,iBand)
      WNe = ab_BandWaveNumBnds(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


    NLev = PFTblNTemp - 1
    allocate( a_PFTblChkTemp    ( NLev ) )
    allocate( a_PFTblChkVal     ( NLev ) )
    allocate( ab_PFTblChkPF     ( NLev, PFTblNBand ) )
    allocate( ab_PFTblChkDPFDT  ( NLev, PFTblNBand ) )
    allocate( a_PFTblChkPFSum   ( NLev ) )
    allocate( a_PFTblChkDPFDTSum( NLev ) )


    do iTemp = 1, NLev
      a_PFTblChkTemp(iTemp) = ( t_PFTblTemp(iTemp) + t_PFTblTemp(iTemp+1) ) / 2.0d0
    end do

    do iBand = 1, PFTblNBand

      WNs = ab_BandWaveNumBnds(1,iBand)
      WNe = ab_BandWaveNumBnds(2,iBand)
      !
      GQNum = InqPFIntGQNum( WNe - WNs )

      call Integ_PF_GQ_Array1D( &
        & WNs, WNe, GQNum, &
        & 1, NLev, &
        & a_PFTblChkTemp, &
        & ab_PFTblChkPF(:,iBand) &
        & )
      ab_PFTblChkPF(:,iBand) = PI * ab_PFTblChkPF(:,iBand) / ( WNe - WNs )
      !
      call Integ_DPFDT_GQ_Array1D( &
        & WNs, WNe, GQNum,         & ! (in )
        & 1, NLev,           &
        & a_PFTblChkTemp,             &
        & ab_PFTblChkDPFDT(:,iBand) &
        & )
      ab_PFTblChkDPFDT(:,iBand) = PI * ab_PFTblChkDPFDT(:,iBand) / ( WNe - WNs )

    end do

    a_PFTblChkPFSum    = 0.0d0
    a_PFTblChkDPFDTSum = 0.0d0
    !
    do iBand = 1, PFTblNBand
      WNs = ab_BandWaveNumBnds(1,iBand)
      WNe = ab_BandWaveNumBnds(2,iBand)
      !
      a_PFTblChkPFSum    = a_PFTblChkPFSum    + ab_PFTblChkPF   (:,iBand) * ( WNe - WNs )
      a_PFTblChkDPFDTSum = a_PFTblChkDPFDTSum + ab_PFTblChkDPFDT(:,iBand) * ( WNe - WNs )
    end do
    WNs = ab_BandWaveNumBnds(1,1         )
    WNe = ab_BandWaveNumBnds(2,PFTblNBand)
    a_PFTblChkPFSum    = a_PFTblChkPFSum    / ( WNe - WNs )
    a_PFTblChkDPFDTSum = a_PFTblChkDPFDTSum / ( WNe - WNs )

    do iBand = 1, PFTblNBand

      call PFTblInterpolatePFDPFDT1D( &
        & PFTblNTemp, PFTblNBand, t_PFTblTemp, tb_PFTblPF, &
        & iBand, NLev, a_PFTblChkTemp, &
        & a_PFTblChkVal &
        & )
      !
      do k = 1, NLev
        if ( ( ab_PFTblChkPF(k,iBand) / a_PFTblChkPFSum(k) ) > PFTblThreshold ) then
          if ( ( abs( ab_PFTblChkPF(k,iBand) - a_PFTblChkVal(k) ) / ab_PFTblChkPF(k,iBand) ) > PFTblRelError ) then
            write( 6, * ) 'Checking PF table: ', iBand, a_PFTblChkTemp(k)
            write( 6, * ) '   Values         : ', &
              & ab_PFTblChkPF(k,iBand), a_PFTblChkVal(k)
            write( 6, * ) '   Relative error : ', &
              & abs( ab_PFTblChkPF(k,iBand) - a_PFTblChkVal(k) ) / a_PFTblChkVal(k)
            stop
          end if
        end if
      end do

      call PFTblInterpolatePFDPFDT1D( &
        & PFTblNTemp, PFTblNBand, t_PFTblTemp, tb_PFTblDPFDT, &
        & iBand, NLev, a_PFTblChkTemp, &
        & a_PFTblChkVal &
        & )
      !
      do k = 1, NLev
        if ( ( ab_PFTblChkDPFDT(k,iBand) / a_PFTblChkDPFDTSum(k) ) > PFTblThreshold ) then
          if ( ( abs( ab_PFTblChkDPFDT(k,iBand) - a_PFTblChkVal(k) ) / ab_PFTblChkDPFDT(k,iBand) ) > PFTblRelError ) then
            write( 6, * ) 'Checking DPFDT table: ', iBand, a_PFTblChkTemp(k)
            write( 6, * ) '   Values         : ', &
              & ab_PFTblChkDPFDT(k,iBand), a_PFTblChkVal(k)
            write( 6, * ) '   Relative error : ', &
              & abs( ab_PFTblChkDPFDT(k,iBand) - a_PFTblChkVal(k) ) / ab_PFTblChkDPFDT(k,iBand)
            stop
          end if
        end if
      end do

    end do

    deallocate( a_PFTblChkTemp     )
    deallocate( a_PFTblChkVal      )
    deallocate( ab_PFTblChkPF      )
    deallocate( ab_PFTblChkDPFDT   )
    deallocate( a_PFTblChkPFSum    )
    deallocate( a_PFTblChkDPFDTSum )


  end subroutine PFTblInit

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

  subroutine PFTblFinalize

    ! Local variables
    !


    deallocate( t_PFTblTemp    )
    deallocate( b_PFTblBandNum )
    deallocate( tb_PFTblPF     )
    deallocate( tb_PFTblDPFDT  )


  end subroutine PFTblFinalize

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

  subroutine PFTblInterpolatePFDPFDT1D( &
    & PFTblNTemp, PFTblNBand, t_PFTblTemp, tb_PFTblVal, &
    & iBand, NLev, a_Temp, &
    & a_Val &
    & )

    integer , intent(in ) :: PFTblNTemp
    integer , intent(in ) :: PFTblNBand
    real(dp), intent(in ) :: t_PFTblTemp(PFTblNTemp)
    real(dp), intent(in ) :: tb_PFTblVal(PFTblNTemp, PFTblNBand)
    integer , intent(in ) :: iBand
    integer , intent(in ) :: NLev
    real(dp), intent(in ) :: a_Temp(NLev)
    real(dp), intent(out) :: a_Val (NLev)


    ! 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

    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_PFTblVal(iTemp-1,iBand)
      a_Y(2) = tb_PFTblVal(iTemp  ,iBand)
      a_Y(3) = tb_PFTblVal(iTemp+1,iBand)
      a_Val(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 subroutine PFTblInterpolatePFDPFDT1D

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

end program main
