module cia_module

  use vtype_module

  implicit none

  private

  public :: CIAParam
  public :: CIAReadFile
  public :: CIAParamFinalize
  public :: CIAWriteFileNetcdf
  public :: CIAReadFileNetcdf
  public :: CIACalcAC
  public :: CIAGetMolNumIndex


    ! IMol            : molecular number for specie
    ! nline           : number of lines
    ! iiso            : isotopomer index
    ! wnlc0(nline)    : line center wave number (m^-1)
    ! stren0(nline)   : line strength (m^-1/(molecule m^-2) @ 296K)
    ! lowene(nline)   : lower state energy (m^-1)
    ! tdep_abw(nline) : coefficient of temperature dependence of
    !                 : air-broadened halfwidth
    ! tdep_sbw(nline) : coefficient of temperature dependence of
    !                 : self-broadened halfwidth
    ! abwid0(nline)   : air-broadened half width (m^-1 @ 296K)
    ! sbwid0(nline)   : self-broadened half width (m^-1 @ 296K)
    ! apshift(nline)  : air pressure-induced line shift (m^-1 atm^-1 @ 296K)

    ! wnlc1(linen)    : wave number at line center after temperature/pressure
    !                 : correction (m^-1) 
    ! stren1(nline)     : line strength after temperature/pressure correction
    !                   : (m^-1 / molecule m^-2)
    ! pbwid1(nline)     : pressure broadened (lorentz) halfwidth
    !                   : after temperature/pressure correction (m^-1)
    ! dowid1(nline)     : doppler halfwidth
    !                   : after temperature/pressure correction (m^-1)
    ! abwid1(nline)     : air-broadened half width
    !                   : after temperature/pressure correction (m^-1)
    ! sbwid1(nline)     : self-broadened half width
    !                   : after temperature/pressure correction (m^-1)

    ! wnco(nline)     : cutoff wave number (m^-1)
    ! renfac(nline)   : renormalization factor

  type :: CIAParam
    integer(i4b)          :: IMol
    integer(i4b)          :: IMolOpp
    integer(i4b)          :: NWaveNum
    integer(i4b)          :: NTemp
    real(dp)    , pointer :: w_WaveNum(:)
    real(dp)    , pointer :: t_Temp   (:)
    real(dp)    , pointer :: wt_Coef  (:,:)
!!$    integer(i4b)          :: km
!!$    real(dp)    , pointer :: wz_Coef(:,:)
  end type CIAParam


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

contains

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

  subroutine CIAReadFile( FileName, ciap )

    use fi_module

    character(*)    , intent(in ) :: FileName
    type( CIAParam ), intent(out) :: ciap


    ! Local variables
    !
    character(stdstr) :: ChemSymSave
    integer :: NWaveNumSave
    integer :: NTemp

    character(stdstr) :: ChemSym
    real(dp) :: WaveNumS
    real(dp) :: WaveNumE
    integer  :: NWaveNum
    real(dp) :: Temp
    real(dp) :: Press
    real(dp) :: MaxCrossSec

    real(dp), allocatable :: w_WaveNum(:)
    real(dp)              :: WaveNum
    real(dp)              :: CIACoef

    integer :: IMol1
    integer :: IMol2

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

    integer            :: l
    integer            :: t


    NWaveNumSave = 0
    NTemp        = 0

    Mode = "read"
    call fi_open( FileName, Mode, FileUnit )
    do
      read( FileUnit, '(a20,2f10.3,i7,f7.1,e10.3,f10.3)', iostat = iostat ) &
        & ChemSym, WaveNumS, WaveNumE, NWaveNum, Temp, Press, MaxCrossSec
      if( iostat .ne. 0 ) exit
      NTemp      = NTemp + 1
      ! check
      if ( NTemp == 1 ) then
        ChemSymSave = ChemSym
        allocate( w_WaveNum(NWaveNum) )
        NWaveNumSave = NWaveNum
      else
        if ( ChemSymSave /= ChemSym ) then
          write( 6, * ) 'Molecular symbol is inappropriate: ', ChemSym, ChemSymSave
          stop
        end if
        if ( NWaveNumSave /= NWaveNum ) then
          write( 6, * ) 'Number of wavenumber is inappropriate: ', NWaveNum, NWaveNumSave
          stop
        end if
      end if
      do l = 1, NWaveNum
        read( FileUnit, * ) WaveNum, CIACoef
        if ( NTemp == 1 ) then
          w_WaveNum(l) = WaveNum
        else
          if ( w_WaveNum(l) /= WaveNum ) then
            write( 6, * ) 'Wavenumber is inappropriate: ', w_WaveNum(l), WaveNum
            stop
          end if
        end if
      end do
    end do
    close( FileUnit )

    deallocate( w_WaveNum )


    call CIAGetIMols( &
      & ChemSym, &
      & IMol1, IMol2 &
      & )

    CIAP % IMol     = IMol1
    CIAP % IMolOpp  = IMol2
    CIAP % NWaveNum = NWaveNum
    CIAP % NTemp    = NTemp
    allocate( CIAP % w_WaveNum( CIAP % NWaveNum ) )
    allocate( CIAP % t_Temp   ( CIAP % NTemp    ) )
    allocate( CIAP % wt_Coef  ( CIAP % NWaveNum, CIAP % NTemp ) )


    Mode = "read"
    call fi_open( FileName, Mode, FileUnit )
    do t = 1, CIAP % NTemp
      read( FileUnit, '(a20,2f10.3,i7,f7.1,e10.3,f10.3)', iostat = iostat ) &
        & ChemSym, WaveNumS, WaveNumE, NWaveNum, Temp, Press, MaxCrossSec
      CIAP % t_Temp(t) = Temp
      do l = 1, CIAP % NWaveNum
        read( FileUnit, * ) WaveNum, CIACoef
        ! unit of wavenumber is changed from cm-1 to m-1
        ! unit of coefficient is changed from cm5 molecule-2 to m5 molecule-2
        WaveNum = WaveNum * 1.0d2
        CIACoef = CIACoef * 1.0d-10
        CIAP % w_WaveNum(l) = WaveNum
        CIAP % wt_Coef(l,t) = CIACoef
      end do
    end do
    close( FileUnit )


  end subroutine CIAReadFile

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

  subroutine CIAStoreDataInObj( FileName, ciap )

    use fi_module

    character(*)    , intent(in ) :: FileName
    type( CIAParam ), intent(out) :: ciap


    ! Local variables
    !
    character(stdstr) :: ChemSymSave
    integer :: NWaveNumSave
    integer :: NTemp

    character(stdstr) :: ChemSym
    real(dp) :: WaveNumS
    real(dp) :: WaveNumE
    integer  :: NWaveNum
    real(dp) :: Temp
    real(dp) :: Press
    real(dp) :: MaxCrossSec

    real(dp), allocatable :: w_WaveNum(:)
    real(dp)              :: WaveNum
    real(dp)              :: CIACoef

    integer :: IMol1
    integer :: IMol2

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

    integer            :: l
    integer            :: t


    NWaveNumSave = 0
    NTemp        = 0

    Mode = "read"
    call fi_open( FileName, Mode, FileUnit )
    do
      read( FileUnit, '(a20,2f10.3,i7,f7.1,e10.3,f10.3)', iostat = iostat ) &
        & ChemSym, WaveNumS, WaveNumE, NWaveNum, Temp, Press, MaxCrossSec
      if( iostat .ne. 0 ) exit
      NTemp      = NTemp + 1
      ! check
      if ( NTemp == 1 ) then
        ChemSymSave = ChemSym
        allocate( w_WaveNum(NWaveNum) )
        NWaveNumSave = NWaveNum
      else
        if ( ChemSymSave /= ChemSym ) then
          write( 6, * ) 'Molecular symbol is inappropriate: ', ChemSym, ChemSymSave
          stop
        end if
        if ( NWaveNumSave /= NWaveNum ) then
          write( 6, * ) 'Number of wavenumber is inappropriate: ', NWaveNum, NWaveNumSave
          stop
        end if
      end if
      do l = 1, NWaveNum
        read( FileUnit, * ) WaveNum, CIACoef
        if ( NTemp == 1 ) then
          w_WaveNum(l) = WaveNum
        else
          if ( w_WaveNum(l) /= WaveNum ) then
            write( 6, * ) 'Wavenumber is inappropriate: ', w_WaveNum(l), WaveNum
            stop
          end if
        end if
      end do
    end do
    close( FileUnit )

    deallocate( w_WaveNum )


    call CIAGetIMols( &
      & ChemSym, &
      & IMol1, IMol2 &
      & )

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

    CIAP % IMol     = IMol1
    CIAP % IMolOpp  = IMol2
    CIAP % NWaveNum = NWaveNum
    CIAP % NTemp    = NTemp
    allocate( CIAP % w_WaveNum( CIAP % NWaveNum ) )
    allocate( CIAP % t_Temp   ( CIAP % NTemp    ) )
    allocate( CIAP % wt_Coef  ( CIAP % NWaveNum, CIAP % NTemp ) )


    Mode = "read"
    call fi_open( FileName, Mode, FileUnit )
    do t = 1, CIAP % NTemp
      read( FileUnit, '(a20,2f10.3,i7,f7.1,e10.3,f10.3)', iostat = iostat ) &
        & ChemSym, WaveNumS, WaveNumE, NWaveNum, Temp, Press, MaxCrossSec
      CIAP % t_Temp(t) = Temp
      do l = 1, CIAP % NWaveNum
        read( FileUnit, * ) WaveNum, CIACoef
        ! unit of wavenumber is changed from cm-1 to m-1
        ! unit of coefficient is changed from cm5 molecule-2 to m5 molecule-2
        WaveNum = WaveNum * 1.0d2
        CIACoef = CIACoef * 1.0d-10
        CIAP % w_WaveNum(l) = WaveNum
        CIAP % wt_Coef(l,t) = CIACoef
      end do
    end do
    close( FileUnit )


  end subroutine CIAStoreDataInObj

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

  subroutine CIAGetIMols( &
    & ChemSym, &
    & IMol1, IMol2 &
    & )

    use hitranconst

    character(*), intent(in ) :: ChemSym
    integer     , intent(out) :: IMol1
    integer     , intent(out) :: IMol2


    ! Local variables
    !
    character(stdstr) :: CMol1
    character(stdstr) :: CMol2
    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
    CMol1 = CMol(ls:le)

    ls = le+1+1
    le = lm
    CMol2 = CMol(ls:le)


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

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


  end subroutine CIAGetIMols

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

  subroutine CIAParamFinalize( CIAP )

    type( CIAParam ), intent(inout) :: CIAP


    CIAP % IMol     = -1
    CIAP % IMolOpp  = -1
    CIAP % NWaveNum = -1
    CIAP % NTemp    = -1
    deallocate( CIAP % w_WaveNum )
    deallocate( CIAP % t_Temp    )
    deallocate( CIAP % wt_Coef   )


  end subroutine CIAParamFinalize

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

  subroutine CIAWriteFileNetcdf( fn, Src, Ref, CIAP )

    use ni3_module

    character(len=*), intent(in) :: fn
    character(len=*), intent(in) :: Src
    character(len=*), intent(in) :: Ref
    type( CIAParam ), intent(in) :: CIAP


    !
    ! 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, ref = Ref )

    call ni3_put_att( NcID, "_global_", "imol"   , CIAP % IMol    )
    call ni3_put_att( NcID, "_global_", "imolopp", CIAP % IMolOpp )

    Name     = "WaveNum"
    StdName  = Name
    LongName = "wavenumber"
    Unit     = "m-1"
    call ni3_set_dim( NcID, Name, NI3_DOUBLE, CIAP % 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, CIAP % 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     = "m5 (molecule)-2"
    call ni3_def_var( NcID, Name, NI3_DOUBLE, NDims, a_DimNames, &
      & StdName, LongName, Unit )
    call ni3_put_var( NcID, Name, CIAP % 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 CIAWriteFileNetcdf

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

  subroutine CIAReadFileNetcdf( fn, CIAP )

    use ni3_module

    character(len=*), intent(in ) :: fn
    type( CIAParam ), intent(out) :: CIAP


    !
    ! 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"   , CIAP % IMol    )
    call ni3_get_att( NcID, "_global_", "imolopp", CIAP % IMolOpp )

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

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

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

    call ni3_close( NcID )


  end subroutine CIAReadFileNetcdf

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

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

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

  subroutine CIACalcAC( &
    & CIAP, VMROpp, NumDen, Temp, &
    & NWaveNum, w_WaveNum, &
    & w_AbsCoef &
    & )

    type(CIAParam), intent(in ) :: CIAP
    real(8)       , intent(in ) :: VMROpp
    real(8)       , intent(in ) :: NumDen
    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 :: lcia


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


    lcia = 1+1
    do l = 1, NWaveNum

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

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

    end do

    w_AbsCoef = w_AbsCoef * VMROpp * NumDen


  end subroutine CIACalcAC

  !****************************************************************************
!!$
!!$    subroutine lbl_readfile_netcdf( fn, km, lp )
!!$
!!$      use ni3_module
!!$
!!$      character(len=*) , intent(in ) :: fn
!!$      integer          , intent(in ) :: km
!!$      type( LineParam ), intent(out) :: lp
!!$
!!$
!!$      !
!!$      ! local variables
!!$      !
!!$      integer(i4b)      :: NcID
!!$      character(stdstr) :: Mode
!!$      character(stdstr) :: Name
!!$      integer(i4b)      :: IMol
!!$      integer(i4b)      :: NLine
!!$      integer(i4b)      :: l
!!$
!!$
!!$      Mode = "read"
!!$      call ni3_open( fn, Mode, NcID )
!!$
!!$      ! IMol            : molecular number for specie
!!$      Name     = "imol"
!!$      call ni3_get_var( NcID, Name, IMol )
!!$
!!$      Name     = "Line"
!!$      call ni3_inq_dimlen( NcID, Name, NLine )
!!$
!!$      call lbl_lineparam_init( IMol, NLine, km, lp )
!!$
!!$      ! iiso            : isotopomer index
!!$      Name     = "iiso"
!!$      call ni3_get_var( NcID, Name, lp%iiso )
!!$
!!$      ! wnlc0(nline)    : line center wave number (m^-1)
!!$      Name     = "wnlc0"
!!$      call ni3_get_var( NcID, Name, lp%wnlc0 )
!!$
!!$      ! stren0(nline)   : line strength (m^-1/(molecule m^-2) @ 296K)
!!$      Name     = "stren0"
!!$      call ni3_get_var( NcID, Name, lp%stren0 )
!!$
!!$      ! lowene(nline)   : lower state energy (m^-1)
!!$      Name     = "lowene0"
!!$      call ni3_get_var( NcID, Name, lp%lowene )
!!$
!!$      ! tdep_abw(nline) : coefficient of temperature dependence of
!!$      !                 : air-broadened halfwidth
!!$      Name     = "tdep_abw"
!!$      call ni3_get_var( NcID, Name, lp%tdep_abw )
!!$
!!$      ! tdep_sbw(nline) : coefficient of temperature dependence of
!!$      !                 : self-broadened halfwidth
!!$      Name     = "tdep_sbw"
!!$      call ni3_get_var( NcID, Name, lp%tdep_sbw )
!!$
!!$      ! abwid0(nline)   : air-broadened half width (m^-1 @ 296K)
!!$      Name     = "abwid0"
!!$      call ni3_get_var( NcID, Name, lp%abwid0 )
!!$
!!$      ! sbwid0(nline)   : self-broadened half width (m^-1 @ 296K)
!!$      Name     = "sbwid0"
!!$      call ni3_get_var( NcID, Name, lp%sbwid0 )
!!$
!!$      ! apshift(nline)  : air pressure-induced line shift (m^-1 atm^-1 @ 296K)
!!$      Name     = "apshift"
!!$      call ni3_get_var( NcID, Name, lp%apshift )
!!$
!!$      call ni3_close( NcID )
!!$
!!$
!!$    end subroutine lbl_readfile_netcdf

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

end module cia_module