module opt_prop_kd

  use dc_types
  use gtool_history

  implicit none

  private

  public :: OptPropKDInit
  public :: OptPropKDFinalize
  public :: OptPropKDGetNWN
  public :: OptPropKDGetNPress
  public :: OptPropKDGetNMol
  public :: OptPropKDGetNBand
  public :: OptPropKDGetDelWN
  public :: OptPropKDGetBandWNBnds
  public :: OptPropKDGetBandBinIndexBnds
  public :: OptPropKDGetAbsCoefProf
  public :: OptPropKDGetPFIntedRatio


  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 :: KDTblNBin
  integer              , save :: KDTblNBand
  integer              , save :: KDTblNPress
  integer              , save :: KDTblNMol
  integer              , save :: KDTblNTemp
  integer              , save :: KDTblNVMRH2O
  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 :: v_KDTblVMRH2O(:)
  real(DP), allocatable, save :: ct_PFRatio   (:,:)
  real(DP), allocatable, save :: pmctv_AbsCoef(:,:,:,:,:)


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

  logical, save :: FlagInited = .false.

contains

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

  subroutine OptPropKDInit( &
    & OptPropNcFN &
    & )

    use rad_kd_utils, only : InqDimLength


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


    ! local variables
    !
    character(128) :: PressAxisName
    character(128) :: MolAxisName
    character(128) :: TempAxisName
    character(128) :: VMRH2OAxisName
    character(128) :: BinAxisName
    character(128) :: AbsCoefVarName
    character(128) :: PFRatioVarName


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


    KDTblNBin    = -1
    KDTblNBand   = -1

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


    OptPropNcFNSave = OptPropNcFN

    PressAxisName   = "Press"
    MolAxisName     = "MolNum"
    TempAxisName    = "Temp"
    VMRH2OAxisName  = "VMRH2O"
    BinAxisName     = 'BinNum'
    AbsCoefVarName  = 'AbsCoef'
    PFRatioVarName  = 'PFRatio'


    KDTblNBin    = InqDimLength( OptPropNcFNSave, BinAxisName     )
    KDTblNBand   = InqDimLength( OptPropNcFNSave, 'BandNum'       )

    KDTblNPress  = InqDimLength( OptPropNcFNSave, PressAxisName   )
    KDTblNMol    = InqDimLength( OptPropNcFNSave, MolAxisName     )
    KDTblNTemp   = InqDimLength( OptPropNcFNSave, TempAxisName    )
    KDTblNVMRH2O = InqDimLength( OptPropNcFNSave, VMRH2OAxisName  )


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

    call HistoryGet( OptPropNcFNSave, 'BandWaveNumBnds', ab_KDTblBandWaveNumBnds, &
      & quiet = .true. )
    call HistoryGet( OptPropNcFNSave, 'BandBinIndexBnds', ab_KDTblBandBinIndexBnds, &
      & quiet = .true. )

    allocate( c_DelWaveNum(KDTblNBin) )
    call HistoryGet( OptPropNcFNSave, "BinWeight", c_DelWaveNum, &
      & quiet = .true. )
    !
    allocate( c_KDTblBinNum(KDTblNBin   ) )
    allocate( p_KDTblPress (KDTblNPress ) )
    allocate( m_KDTblMolNum(KDTblNMol   ) )
    allocate( t_KDTblTemp  (KDTblNTemp  ) )
    allocate( v_KDTblVMRH2O(KDTblNVMRH2O) )
    !
    call HistoryGet( OptPropNcFNSave, BinAxisName   , c_KDTblBinNum, &
      & quiet = .true. )
    call HistoryGet( OptPropNcFNSave, PressAxisName , p_KDTblPress, &
      & quiet = .true. )
    call HistoryGet( OptPropNcFNSave, MolAxisName   , m_KDTblMolNum, &
      & quiet = .true. )
    call HistoryGet( OptPropNcFNSave, TempAxisName  , t_KDTblTemp, &
      & quiet = .true. )
    call HistoryGet( OptPropNcFNSave, VMRH2OAxisName, v_KDTblVMRH2O, &
      & quiet = .true. )
    !
    allocate( pmctv_AbsCoef(KDTblNPress,KDTblNMol,KDTblNBin,KDTblNTemp,KDTblNVMRH2O) )
    call HistoryGet( OptPropNcFNSave, AbsCoefVarName, pmctv_AbsCoef, &
      & quiet = .true. )
    !
    allocate( ct_PFRatio(KDTblNBin,KDTblNTemp) )
    call HistoryGet( OptPropNcFNSave, PFRatioVarName, ct_PFRatio, &
      & quiet = .true. )


    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
    KDTblNVMRH2O = -1


    deallocate( ab_KDTblBandWaveNumBnds  )
    deallocate( ab_KDTblBandBinIndexBnds )
    !
    deallocate( c_DelWaveNum )
    !
    deallocate( c_KDTblBinNum )
    deallocate( p_KDTblPress  )
    deallocate( m_KDTblMolNum )
    deallocate( t_KDTblTemp   )
    deallocate( v_KDTblVMRH2O )
    !
    deallocate( pmctv_AbsCoef )
    deallocate( ct_PFRatio )


    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

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

  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 OptPropKDGetAbsCoefProf( &
    & iWaveNum, kmax, NMol, r_Press, r_Temp, r_VMRH2O, &
    & rm_AbsCoef &
    & )


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


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


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



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

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

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


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


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

      do m = 1, NMol
        do n3 = 1, 2
          do n2 = 1, 2
            do n1 = 1, 2
              l1 = iPress -2+n1
              l2 = iTemp  -2+n2
              l3 = iVMRH2O-2+n3
!!$              aaam_Val(n1,n2,n3,m) = pmctv_AbsCoef(l1,m,iWaveNum,l2,l3)
              aaam_Val(n1,n2,n3,m) = &
                & log( max( pmctv_AbsCoef(l1,m,iWaveNum,l2,l3), 1.0d-100 ) )
              if ( ( iWaveNum == 1 ) .and. ( m == 3 ) ) then
                write( 71, * ) m, pmctv_AbsCoef(l1,m,iWaveNum,l2,l3)
              end if
            end do
          end do
        end do
      end do
      if ( iWaveNum == 1 ) then
        write( 71, * ) '---'
      end if


      do m = 1, NMol

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


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

!!$        rm_AbsCoef(k,m) = Val
        rm_AbsCoef(k,m) = exp( Val )
      end do
    end do

    if ( iWaveNum == 1 ) then
      m = 3
      do k = 1, kmax
        write( 6, * ) k, rm_AbsCoef(k,m), r_Temp(k), r_VMRH2O(k), r_Press(k)
      end do
    end if

  end subroutine OptPropKDGetAbsCoefProf

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

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


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


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


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


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

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

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

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

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

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

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


  end subroutine OptPropKDGetPFIntedRatio

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

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


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


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


    ! check
    if ( a_TblGrid(1) > a_TblGrid(2) ) then
      write( 6, * ) 'Unexpected order of axis ', trim( AxisName ), &
        & ', i.e., descending order.'
      stop
    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
    end do


  end subroutine FindTblIndex

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

end module opt_prop_kd

