module xsc_module

  use vtype_module

  implicit none

  private

  public :: XSCParam
  public :: XSCStoreDataInObj
  public :: XSCGetIMol
  public :: XSCParamFinalize
  public :: XSCReadHITRANXSCFile
  public :: XSCWriteFileNetcdf
  public :: XSCReadFileNetcdf
  public :: XSCCalcAC
  public :: XSCGetMolNumIndex


  type :: XSCParam
    integer(i4b)          :: IMol
    integer(i4b)          :: IMolOpp
    integer(i4b)          :: NWaveNum
    integer(i4b)          :: NTemp
    real(dp)              :: Press
    real(dp)    , pointer :: w_WaveNum(:)
    real(dp)    , pointer :: t_Temp   (:)
    real(dp)    , pointer :: wt_Coef  (:,:)
  end type XSCParam


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

contains

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

  subroutine XSCReadHITRANXSCFile( &
    & NFile, a_FileName, &
    & WaveNumLimS, WaveNumLimE, &
    & XSCP &
    & )

    use fi_module

    integer          , intent(in ) :: NFile
    character(extstr), intent(in ) :: a_FileName(NFile)
    real(dp)         , intent(in ) :: WaveNumLimS
    real(dp)         , intent(in ) :: WaveNumLimE
    type( XSCParam ) , intent(out) :: XSCP


    ! Local variables
    !
    character(stdstr) :: ChemSymSave

    integer :: NTemp

    character(20) :: a_Word(12)

    character(20) :: ChemSym
    real(dp)      :: WaveNumS
    real(dp)      :: WaveNumE
    integer       :: NWaveNum
    real(dp)      :: Temp
    real(dp)      :: Press
    real(dp)      :: XSCMax
!!$    real(dp)      :: InstRes
    character( 5) :: InstRes
    character(15) :: ComName
    character( 4) :: Blank
    character( 3) :: Broadener
    integer       :: Ref

    real(dp)      :: WaveNumSSave
    real(dp)      :: WaveNumESave
    integer       :: NWaveNumSave
    real(dp)      :: TempSave
    real(dp)      :: PressSave

    real(dp), allocatable :: w_WaveNum(:)
    real(dp), allocatable :: t_Temp   (:)
    real(dp), allocatable :: wt_Coef  (:,:)

    integer               :: iWaveNumS
    integer               :: iWaveNumE
    integer               :: NWaveNumBK
    real(dp), allocatable :: w_WaveNumBK(:)
    real(dp), allocatable :: wt_CoefBK  (:,:)

    integer :: IMol

    character(stdstr)  :: Mode
    integer            :: FileUnit
!!$    integer            :: iostat

    integer            :: NLoop

    integer            :: i
    integer            :: l
    integer            :: n
    integer            :: t


    ! check file
    do n = 1, NFile
      Mode = "read"
      call fi_open( a_FileName(n), Mode, FileUnit )
!!$      read( FileUnit, '(a20,2f10.3,i7,f7.3,f6.3,f10.3,f5.3,a15,a4,a3,i3)' ) &
!!$        & ChemSym, WaveNumS, WaveNumE, NWaveNum, Temp, Press, XSCMax, &
!!$        & InstRes, ComName, Blank, Broadener, Ref
      a_Word = ""
      read( FileUnit, '(a20,2a10,2a7,a6,a10,a5,a15,a4,2a3)' ) &
        & ( a_Word(l), l = 1, 12 )
      read( a_Word( 1), * ) ChemSym
      read( a_Word( 2), * ) WaveNumS
      read( a_Word( 3), * ) WaveNumE
      read( a_Word( 4), * ) NWaveNum
      read( a_Word( 5), * ) Temp
      read( a_Word( 6), * ) Press
      read( a_Word( 7), * ) XSCMax
      read( a_Word( 8), * ) InstRes
      read( a_Word( 9), * ) ComName
!!$      read( a_Word(10), * ) Blank
!!$      read( a_Word(11), * ) Broadener
      read( a_Word(12), * ) Ref
      close( FileUnit )
      if ( n == 1 ) then
        ChemSymSave  = ChemSym
        WaveNumSSave = WaveNumS
        WaveNumESave = WaveNumE
        NWaveNumSave = NWaveNum
        TempSave     = Temp
        PressSave    = Press
      else
        if ( ChemSymSave /= ChemSym ) then
          write( 6, * ) 'ChemSym is not the same for all files.'
          write( 6, * ) ChemSymSave, ChemSym
          stop
        end if
        if ( WaveNumSSave /= WaveNumS ) then
          write( 6, * ) 'WaveNumS is not the same for all files.'
          write( 6, * ) WaveNumSSave, WaveNumS
          stop
        end if
        if ( WaveNumESave /= WaveNumE ) then
          write( 6, * ) 'WaveNumE is not the same for all files.'
          write( 6, * ) WaveNumESave, WaveNumE
          stop
        end if
        if ( NWaveNumSave /= NWaveNum ) then
          write( 6, * ) 'NWaveNum is not the same for all files.'
          write( 6, * ) NWaveNumSave, NWaveNum
          stop
        end if
!!$        if ( TempSave == Temp ) then
!!$          write( 6, * ) 'Temp is same with that in another file.'
!!$          write( 6, * ) TempSave, Temp
!!$          stop
!!$        end if
        if ( PressSave /= Press ) then
          write( 6, * ) 'Press is not the same for all files.'
          write( 6, * ) PressSave, Press
          stop
        end if
      end if
    end do


    NTemp        = NFile

    if ( NTemp <= 1 ) then
      write( 6, * ) 'ERROR: Number of temperature has to be greater than 1.'
      stop
    end if

    allocate( w_WaveNum( NWaveNum ) )
    allocate( t_Temp   ( NTemp    ) )
    allocate( wt_Coef  ( NWaveNum, NTemp ) )

    WaveNumS = WaveNumS * 1.0d2
    WaveNumE = WaveNumE * 1.0d2
    do l = 1, NWaveNum
      w_WaveNum(l) = ( WaveNumE - WaveNumS ) / ( NWaveNum-1 ) * ( l - 1 ) + WaveNumS
    end do

    do n = 1, NFile
      Mode = "read"
      call fi_open( a_FileName(n), Mode, FileUnit )
      !
!!$      read( FileUnit, '(a20,2f10.3,i7,f7.3,f6.3,f10.3,f5.3,a15,a4,a3,i3)' ) &
!!$        & ChemSym, WaveNumS, WaveNumE, NWaveNum, Temp, Press, XSCMax, &
!!$        & InstRes, ComName, Blank, Broadener, Ref
      a_Word = ""
      read( FileUnit, '(a20,2a10,a7,a7,a6,a10,a5,a15,a4,a3,a3)' ) &
        & ( a_Word(l), l = 1, 12 )
      read( a_Word( 1), * ) ChemSym
      read( a_Word( 2), * ) WaveNumS
      read( a_Word( 3), * ) WaveNumE
      read( a_Word( 4), * ) NWaveNum
      read( a_Word( 5), * ) Temp
      read( a_Word( 6), * ) Press
      read( a_Word( 7), * ) XSCMax
      read( a_Word( 8), * ) InstRes
      read( a_Word( 9), * ) ComName
!!$      read( a_Word(10), * ) Blank
!!$      read( a_Word(11), * ) Broadener
      read( a_Word(12), * ) Ref
      !
      t_Temp(n) = Temp
      !
      NLoop = NWaveNum/10
      if ( mod( NWaveNum, 10 ) >= 1 ) NLoop = NLoop + 1
      do i = 1, NLoop
        read( FileUnit, * ) &
          & ( wt_Coef(l,n), l = (i-1)*10+1, min((i-1)*10+10,NWaveNum) )
      end do
      !
      close( FileUnit )
    end do

    if ( ( NTemp         == 2             ) .and. &
      &  ( t_Temp(1)     == t_Temp(2)     ) .and. &
      &  ( a_FileName(1) == a_FileName(2) ) ) then
      write( 6, * ) 'NOTE: Temperature is modified: '
      write( 6, * ) '  Original values: ', t_Temp
      t_Temp(2) = t_Temp(2) + 1.0d0
      write( 6, * ) '  Modified valies: ', t_Temp
    else
      do t = 1+1, NTemp
        if ( t_Temp(t-1) >= t_Temp(t) ) then
          write( 6, * ) 'Order of file (temperature) is inappropriate.'
          write( 6, * ) t_Temp
          stop
        end if
      end do
    end if

    ! unit of coefficient is changed from cm2 molecule-1 to m2 molecule-1
    wt_Coef = wt_Coef * 1.0d-4

    !
    ! Restrict wavenumber range
    !
    NWaveNumBK = NWaveNum
    allocate( w_WaveNumBK(NWaveNumBK)       )
    allocate( wt_CoefBK  (NWaveNumBK,NTemp) )
    w_WaveNumBK = w_WaveNum
    wt_CoefBK   = wt_Coef
    deallocate( w_WaveNum )
    deallocate( wt_Coef   )
    !
    NWaveNum = 0
    do l = 1, NWaveNumBK
      if ( WaveNumLimS <= w_WaveNumBK(l) ) exit
    end do
    iWaveNumS = l
    do l = NWaveNumBK, 1, -1
      if ( w_WaveNumBK(l) <= WaveNumLimE ) exit
    end do
    iWaveNumE = l
    !
    NWaveNum = iWaveNumE - iWaveNumS + 1
    allocate( w_WaveNum(NWaveNum)       )
    allocate( wt_Coef  (NWaveNum,NTemp) )
    do t = 1, NTemp
      do l = 1, NWaveNum
        w_WaveNum(l)   = w_WaveNumBK(l+iWaveNumS-1)
        wt_Coef  (l,t) = wt_CoefBK  (l+iWaveNumS-1,t)
      end do
    end do
    deallocate( w_WaveNumBK )
    deallocate( wt_CoefBK   )


    call XSCGetIMol( &
      & ChemSym, &
      & IMol &
      & )

    call XSCStoreDataInObj( &
      & IMol, &
      & NWaveNum, NTemp, Press, w_WaveNum, t_Temp, wt_Coef, &
      & XSCP  &
      & )

    deallocate( w_WaveNum )
    deallocate( t_Temp    )
    deallocate( wt_Coef   )


  end subroutine XSCReadHITRANXSCFile

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

  subroutine XSCStoreDataInObj( &
    & IMol, &
    & NWaveNum, NTemp, Press, w_WaveNum, t_Temp, wt_Coef, &
    & XSCP  &
    & )

    integer , intent(in ) :: IMol
    integer , intent(in ) :: NWaveNum
    integer , intent(in ) :: NTemp
    real(dp), intent(in ) :: Press
    real(dp), intent(in ) :: w_WaveNum(NWaveNum)
    real(dp), intent(in ) :: t_Temp   (NTemp   )
    real(dp), intent(in ) :: wt_Coef  (NWaveNum,NTemp)
    type( XSCParam ), intent(out) :: XSCP


    ! Local variables
    !


    XSCP % IMol     = IMol
    XSCP % NWaveNum = NWaveNum
    XSCP % NTemp    = NTemp

    allocate( XSCP % w_WaveNum( XSCP % NWaveNum ) )
    allocate( XSCP % t_Temp   ( XSCP % NTemp    ) )
    allocate( XSCP % wt_Coef  ( XSCP % NWaveNum, XSCP % NTemp ) )

    XSCP % Press     = Press
    XSCP % w_WaveNum = w_WaveNum
    XSCP % t_Temp    = t_Temp
    XSCP % wt_Coef   = wt_Coef


  end subroutine XSCStoreDataInObj

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

  subroutine XSCGetIMol( &
    & ChemSym, &
    & IMol &
    & )

    use hitranconst, only : &
      & NHITRANMolNum, &
      & m_HITRANMolName

    character(*), intent(in ) :: ChemSym
    integer     , intent(out) :: IMol


    ! Local variables
    !
    character(stdstr) :: CMol
    integer :: lm
    integer :: ls
    integer :: le

    integer :: l
    integer :: m


    CMol  = trim( ChemSym )
    lm = len_trim( CMol )
    do l = 1, lm
      if ( CMol(l:l) /= ' ' ) exit
    end do
    ls = l
    do l = ls+1, lm
      if ( CMol(l:l) == '-' ) exit
    end do
    le = l-1
    CMol = CMol(ls:le)


    do m = 1, NHITRANMolNum
      if ( CMol == m_HITRANMolName(m) ) exit
    end do
    if ( m > NHITRANMolNum ) then
      write( 6, * ) 'Unable to find molecular name, ', trim( CMol ), ', in HITRAN.'
      stop
    end if
    IMol = m


  end subroutine XSCGetIMol

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

  subroutine XSCParamFinalize( XSCP )

    type( XSCParam ), intent(inout) :: XSCP


    XSCP % IMol     = -1
    XSCP % IMolOpp  = -1
    XSCP % NWaveNum = -1
    XSCP % NTemp    = -1
    XSCP % Press    = -1.0d0
    deallocate( XSCP % w_WaveNum )
    deallocate( XSCP % t_Temp    )
    deallocate( XSCP % wt_Coef   )


  end subroutine XSCParamFinalize

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

  subroutine XSCWriteFileNetcdf( fn, Src, Com, Ref, XSCP )

    use ni3_module

    character(len=*), intent(in) :: fn
    character(len=*), intent(in) :: Src
    character(len=*), intent(in) :: Com
    character(len=*), intent(in) :: Ref
    type( XSCParam ), intent(in) :: XSCP


    !
    ! local variables
    !
    integer(i4b)                   :: NDims
    character(stdstr), allocatable :: a_DimNames(:)
    integer(i4b)                   :: NcID
    character(stdstr)              :: Mode
    character(stdstr)              :: Name
    character(stdstr)              :: StdName
    character(extstr)              :: LongName
    character(extstr)              :: Unit

!!$    integer :: t



    Mode = "new"
    call ni3_open( fn, Mode, NcID )

    call ni3_set_ga( NcID, src = Src, com = Com, ref = Ref )

    call ni3_put_att( NcID, "_global_", "imol", XSCP % IMol )

    Name     = "WaveNum"
    StdName  = Name
    LongName = "wavenumber"
    Unit     = "m-1"
    call ni3_set_dim( NcID, Name, NI3_DOUBLE, XSCP % w_WaveNum, &
      & StdName, LongName, Unit )
!      call ni3_def_dim( NcID, Name, NI3_INT, NI3_UNLIMITED, &
!        & StdName, LongName, Unit )

    Name     = "Temp"
    StdName  = Name
    LongName = "temperature"
    Unit     = "K"
    call ni3_set_dim( NcID, Name, NI3_DOUBLE, XSCP % t_Temp, &
      & StdName, LongName, Unit )
!!$    call ni3_def_dim( NcID, Name, NI3_DOUBLE, NI3_UNLIMITED, &
!!$      & StdName, LongName, Unit )
!!$    do t = 1, CIAP % NTemp
!!$      call ni3_put_varss( NcID, Name, t, CIAP % t_Temp(t) )
!!$    end do

    NDims = 2
    allocate( a_DimNames( NDims ) )
    a_DimNames(1) = "WaveNum"
    a_DimNames(2) = "Temp"

    Name     = "Coef"
    StdName  = Name
    LongName = "intensity coefficient"
    Unit     = "m2 (molecule)-1"
    call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
      & StdName, LongName, Unit )
    call ni3_put_var( NcID, Name, XSCP % wt_Coef )
!!$    do t = 1, CIAP % NTemp
!!$      write( 6, * ) t, CIAP % wt_Coef(:,t)
!!$      call ni3_put_varss( NcID, Name, t, CIAP % wt_Coef(:,t) )
!!$    end do


    deallocate( a_DimNames )

    call ni3_close( NcID )


  end subroutine XSCWriteFileNetcdf

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

  subroutine XSCReadFileNetcdf( fn, XSCP )

    use ni3_module

    character(len=*), intent(in ) :: fn
    type( XSCParam ), intent(out) :: XSCP


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


    Mode = "read"
    call ni3_open( fn, Mode, NcID )

    call ni3_get_att( NcID, "_global_", "imol", XSCP % IMol )

    Name     = "WaveNum"
    call ni3_inq_dimlen( NcID, Name, XSCP % NWaveNum )
    Name     = "Temp"
    call ni3_inq_dimlen( NcID, Name, XSCP % NTemp    )

    allocate( XSCP % w_WaveNum( XSCP % NWaveNum ) )
    allocate( XSCP % t_Temp   ( XSCP % NTemp    ) )
    allocate( XSCP % wt_Coef  ( XSCP % NWaveNum, XSCP % NTemp ) )

    Name     = "WaveNum"
    call ni3_get_var( NcID, Name, XSCP % w_WaveNum )
    Name     = "Temp"
    call ni3_get_var( NcID, Name, XSCP % t_Temp    )
    Name     = "Coef"
    call ni3_get_var( NcID, Name, XSCP % wt_Coef   )

    call ni3_close( NcID )


  end subroutine XSCReadFileNetcdf

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

  function XSCGetMolNumIndex( NMolNum, m_MolNum, IMol ) result( Index )

    integer, intent(in ) :: NMolNum
    integer, intent(in ) :: m_MolNum( NMolNum )
    integer, intent(in ) :: IMol

    integer :: Index

    ! Local variables
    !
    integer :: m


    do m = 1, NMolNum
      if ( IMol == m_MolNum(m) ) exit
    end do
    if ( m > NMolNum ) then
      write( 6, * ) 'Unable to find molecular number ', IMol, ' in ', m_MolNum
      stop
    end if

    Index = m


  end function XSCGetMolNumIndex

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

  subroutine XSCCalcAC( &
    & XSCP, Temp, &
    & NWaveNum, w_WaveNum, &
    & w_AbsCoef &
    & )

    type(XSCParam), intent(in ) :: XSCP
    real(8)       , intent(in ) :: Temp
    integer       , intent(in ) :: NWaveNum
    real(8)       , intent(in ) :: w_WaveNum(NWaveNum)
    real(8)       , intent(out) :: w_AbsCoef(NWaveNum)


    ! Local variables
    !
    real(8) :: a_WeightWaveNum(2)
    real(8) :: a_WeightTemp   (2)
    integer :: iTemp
    integer :: l
    integer :: l2
    integer :: lxsc


    if ( Temp <= XSCP%t_Temp(1) ) then
      iTemp           = 2
      a_WeightTemp(1) = 1.0d0
      a_WeightTemp(2) = 0.0d0
    else if ( XSCP%t_Temp( XSCP%NTemp ) <= Temp ) then
      iTemp           = XSCP%NTemp
      a_WeightTemp(1) = 0.0d0
      a_WeightTemp(2) = 1.0d0
    else
      do iTemp = 1+1, XSCP%NTemp
        if ( Temp <= XSCP%t_Temp(iTemp) ) exit
      end do
      if ( iTemp > XSCP%NTemp ) then
        stop 'Unexpected error when searching index for temperature.'
      end if
      a_WeightTemp(1) = &
        & ( XSCP%t_Temp(iTemp  ) - Temp ) / ( XSCP%t_Temp(iTemp) - XSCP%t_Temp(iTemp-1) )
      a_WeightTemp(2) = &
        & ( Temp - XSCP%t_Temp(iTemp-1) ) / ( XSCP%t_Temp(iTemp) - XSCP%t_Temp(iTemp-1) )
    end if


    lxsc = 1+1
    do l = 1, NWaveNum

      if ( w_WaveNum(l) < XSCP%w_WaveNum(1) ) then
        w_AbsCoef(l) = 0.0d0
      else if ( XSCP%w_WaveNum(XSCP%NWaveNum) < w_WaveNum(l) ) then
        w_AbsCoef(l) = 0.0d0
      else
        do l2 = lxsc, XSCP%NWaveNum
          if ( w_WaveNum(l) <= XSCP%w_WaveNum(l2) ) exit
        end do
        lxsc = l2
        a_WeightWaveNum(1) = &
          &   ( XSCP%w_WaveNum(lxsc  ) - w_WaveNum(l)           ) &
          & / ( XSCP%w_WaveNum(lxsc  ) - XSCP%w_WaveNum(lxsc-1) )
        a_WeightWaveNum(2) = &
          &   ( w_WaveNum(l)           - XSCP%w_WaveNum(lxsc-1) ) &
          & / ( XSCP%w_WaveNum(lxsc  ) - XSCP%w_WaveNum(lxsc-1) )

        w_AbsCoef(l) = &
          &   XSCP%wt_Coef(lxsc-1,iTemp-1)           &
          &   * a_WeightWaveNum(1) * a_WeightTemp(1) &
          & + XSCP%wt_Coef(lxsc  ,iTemp-1)           &
          &   * a_WeightWaveNum(2) * a_WeightTemp(1) &
          & + XSCP%wt_Coef(lxsc-1,iTemp  )           &
          &   * a_WeightWaveNum(1) * a_WeightTemp(2) &
          & + XSCP%wt_Coef(lxsc  ,iTemp  )           &
          &   * a_WeightWaveNum(2) * a_WeightTemp(2)
      end if

    end do


  end subroutine XSCCalcAC

  !****************************************************************************

end module xsc_module
