c     Copyright (C) 2008-2014 Vincent Eymet
c
c     KDISTRIBUTION is free software; you can redistribute it and/or modify
c     it under the terms of the GNU General Public License as published by
c     the Free Software Foundation; either version 3, or (at your option)
c     any later version.

c     KDISTRIBUTION is distributed in the hope that it will be useful,
c     but WITHOUT ANY WARRANTY; without even the implied warranty of
c     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c     GNU General Public License for more details.
c
c     You should have received a copy of the GNU General Public License
c     along with KDISTRIBUTION; if not, see <http://www.gnu.org/licenses/>
c
      subroutine read_ka(kspectrum,nl,imin,n,nu_min,nu_max,
     &     dkdP_spectrum,dkdT_spectrum,dkdx_spectrum,Nmol,
     &     nu,k,dk_dP,dk_dT,dk_dx)
      implicit none
      include 'max.inc'
c
c     Purpose: to retrieve the k-spectrum in the required spectral interval
c     
c     Inputs:
c       + kspectrum: full nu/k spectrum
c       + nl: number of points in "kspectrum"
c       + imin: index of first value that has to be taken from "kspectrum"
c       + n: number of values that have to be taken from "kspectrum"
c       + nu_min: minimum value of "nu" for the spectral interval
c       + nu_max: maximum value of "nu" for the spectral interval
c       + dkdP_spectrum: values of d(k)/d(P) for the same values of "nu"
c       + dkdT_spectrum: values of d(k)/d(T) for the same values of "nu"
c       + dkdx_spectrum: values of d(k)/d(x) for the same values of "nu"
c       + Nmol: number of molecules in the dkdx_spectrum array
c
c     Outputs:
c       + nu: values of nu for the spectral interval
c       + k: associated values of k
c       + dk_dP: values of d(k)/d(P)
c       + dk_dT: values of d(k)/d(T)
c       + dk_dx: values of d(k)/d(x)
c
      integer i,imin,n,nl,k1i,kni
      double precision kspectrum(1:kmax,2)
      double precision dkdP_spectrum(1:kmax)
      double precision dkdT_spectrum(1:kmax)
      double precision dkdx_spectrum(1:kmax,1:Nmol_mx)
      integer Nmol,mol
      double precision nu_min,nu_max
      double precision nu_inf,nu_sup,k_inf,k_sup
      double precision nu(1:kmax),k(1:kmax)
      double precision dk_dP(1:kmax)
      double precision dk_dT(1:kmax)
      double precision dk_dx(1:kmax,1:Nmol_mx)
      double precision dkdP_inf,dkdP_sup
      double precision dkdT_inf,dkdT_sup
      double precision dkdx_inf(1:Nmol_mx),dkdx_sup(1:Nmol_mx)
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine read_ka'

c     Debug
c      write(*,*) 'imin=',imin,' n=',n
c      write(*,*) 'k(',imin,')=',kspectrum(imin,1),kspectrum(imin,2)
c      write(*,*) 'k(',imin+n-1,')=',kspectrum(imin+n-1,1),
c     &     kspectrum(imin+n-1,2)
c      write(*,*) 'dk_dT(',imin,')=',dkdT_spectrum(imin)
c     Debug
      k1i=0
      nu(1)=nu_min
      k(1)=kspectrum(imin,2)
      dk_dP(1)=dkdP_spectrum(imin)
      dk_dT(1)=dkdT_spectrum(imin)
      do mol=1,Nmol
         dk_dx(1,mol)=dkdx_spectrum(imin,mol)
      enddo ! mol
      if (imin.gt.1) then
         if ((nu_min.gt.kspectrum(imin-1,1))
     &        .and.(nu_min.lt.kspectrum(imin,1))) then
            k1i=1
            nu_inf=kspectrum(imin-1,1)
            k_inf=kspectrum(imin-1,2)
            dkdP_inf=dkdP_spectrum(imin-1)
            dkdT_inf=dkdT_spectrum(imin-1)
            do mol=1,Nmol
               dkdx_inf(mol)=dkdx_spectrum(imin-1,mol)
            enddo
            nu(2)=kspectrum(imin,1)
            k(2)=kspectrum(imin,2)
            dk_dP(2)=dkdP_spectrum(imin)
            dk_dT(2)=dkdT_spectrum(imin)
            do mol=1,Nmol
               dk_dx(2,mol)=dkdx_spectrum(imin,mol)
            enddo ! mol
            k(1)=k_inf+(k(2)-k_inf)/(nu(2)-nu_inf)*(nu_min-nu_inf)
            dk_dP(1)=dkdP_inf+
     &           (dk_dP(2)-dkdP_inf)
     &           /(nu(2)-nu_inf)*(nu_min-nu_inf)
            dk_dT(1)=dkdT_inf+
     &           (dk_dT(2)-dkdT_inf)
     &           /(nu(2)-nu_inf)*(nu_min-nu_inf)
            do mol=1,Nmol
               dk_dx(1,mol)=dkdx_inf(mol)+
     &              (dk_dx(2,mol)-dkdx_inf(mol))
     &              /(nu(2)-nu_inf)*(nu_min-nu_inf)
            enddo ! mol
         endif
      endif
c     Debug
c      write(*,*) 'nu(1)=',nu(1)
c      write(*,*) 'k(1)=',k(1)
c     Debug
      do i=2+k1i,n+1
         nu(i)=kspectrum(imin+i-1-k1i,1)
         k(i)=kspectrum(imin+i-1-k1i,2)
         dk_dP(i)=dkdP_spectrum(imin+i-1-k1i)
         dk_dT(i)=dkdT_spectrum(imin+i-1-k1i)
         do mol=1,Nmol
            dk_dx(i,mol)=dkdx_spectrum(imin+i-1-k1i,mol)
         enddo ! mol
      enddo
      kni=0
      if ((imin+n-1.lt.nl).or.(n.eq.0)) then
         if ((nu_max.gt.kspectrum(imin+n-1,1))
     &        .and.(nu_max.lt.kspectrum(imin+n,1))) then
            kni=1
            nu_sup=kspectrum(imin+n,1)
            k_sup=kspectrum(imin+n,2)
            dkdP_sup=dkdP_spectrum(imin+n)
            dkdT_sup=dkdT_spectrum(imin+n)
            do mol=1,Nmol
               dkdx_sup(mol)=dkdx_spectrum(imin+n,mol)
            enddo ! mol
            nu(n+2)=nu_max
            k(n+2)=k(n+1)+(k_sup-k(n+1))/(nu_sup-nu(n+1))
     &           *(nu_max-nu(n+1))
            dk_dP(n+2)=dk_dP(n+1)+
     &           (dkdP_sup-dk_dP(n+1))
     &           /(nu_sup-nu(n+1))*(nu_max-nu(n+1))
            dk_dT(n+2)=dk_dT(n+1)+
     &           (dkdT_sup-dk_dT(n+1))
     &           /(nu_sup-nu(n+1))*(nu_max-nu(n+1))
            do mol=1,Nmol
               dk_dx(n+2,mol)=dk_dx(n+1,mol)+
     &              (dkdx_sup(mol)-dk_dx(n+1,mol))
     &              /(nu_sup-nu(n+1))*(nu_max-nu(n+1))
            enddo ! mol
c     Debug
c     write(*,*) 'nu(',n+2,')=',nu(n+2)
c     write(*,*) 'k(',n+2,')=',k(n+2)
c     Debug
            n=n+1               ! because nu(n+2) is interpolated
         endif
      endif
      if (k1i.eq.1) then
         n=n+1                  ! because k(1) is interpolated
      endif

c     Debug
c      write(*,*) 'n=',n
c     Debug

c     Debug
      do i=2,n
         if (nu(i).eq.nu(i-1)) then
            call error(label)
            write(*,*) 'nu(',i-1,')=',nu(i-1)
            write(*,*) 'nu(',i,')=',nu(i)
            stop
         endif
      enddo
c     Debug

      return
      end



      subroutine read_single_kspectrum(level,band,
     &     Nlines,Nvalues,Nvalues_read,
     &     nu_start,nu_end,pindex,
     &     nu_prev,k_prev,
     &     nu_tmp,k_tmp,
     &     Nk,nu,k,Nmol,k_is_null)
      implicit none
      include 'max.inc'
      include 'formats.inc'
c     
c     Purpose: to read the spectrum for a single file
c
c     Input:
c       + level: level index
c       + band: band index
c       + Nlines: number of lines in main high-resolution spectrum data file
c       + Nvalues: total number of (nu/k) couples in the file
c       + Nvalues_read: number of (nu/k) couples read from file
c       + nu_start: starting wavenumber
c       + nu_end: ending wavenumber
c       + pindex: index of the process that is running this routine
c       + nu_prev, k_prev: last-1 (nu/k) couple read from file
c       + nu_tmp, k_tmp: last (nu/k) couple read from file
c     
c     Output:
c       + Nk: number of nu/k(nu) values that have been read
c       + nu: values of nu for the spectral interval
c       + k: associated values of k
c       + Nmol: number of molecular species
c       + k_is_null: true if every value of k is null
c       + Nvalues_read: number of (nu/k) couples read from file
c       + nu_prev, k_prev: last-1 (nu/k) couple read from file
c       + nu_tmp, k_tmp: last (nu/k) couple read from file
c

c     I/O
      integer level,band
      integer Nlines,Nvalues,Nvalues_read
      double precision nu_start,nu_end
      integer pindex
      integer Nk
      double precision nu(1:kmax)
      double precision k(1:kmax)
      integer Nmol
      logical k_is_null
      double precision nu_prev,k_prev
      double precision nu_tmp,k_tmp
c     temp
      logical band_done,nu_start_found,nu_end_found
      double precision raf
      integer ik
      double precision epsilon_k
      parameter(epsilon_k=1.0D-10)
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine read_single_kspectrum'

c     Debug
c      write(*,*) 'nu_start=',nu_start
c      write(*,*) 'nu_end=',nu_end
c     Debug

      if (nu_end.le.nu_start) then
         call error(label)
         write(*,*) 'Inconsistent input parameters:'
         write(*,*) 'nu_start=',nu_start
         write(*,*) 'nu_end=',nu_end
         stop
      endif
      
      Nk=0
      band_done=.false.
      nu_start_found=.false.
      nu_end_found=.false.

      do while (.not.band_done)

         if ((.not.nu_start_found).and.
     &        (nu_prev.lt.nu_start).and.
     &        (nu_tmp.ge.nu_start)) then
            nu_start_found=.true.
            Nk=Nk+1
            if (Nk.gt.kmax) then
               call error(label)
               write(*,*) 'Nk=',Nk
               write(*,*) 'while kmax=',kmax
               write(*,*) 'when nu_start has been found'
               stop
            endif
c     Debug
c            write(*,*) 'Nvalues_read=',Nvalues_read
c     Debug
            if (Nvalues_read.eq.1) then
               nu(Nk)=nu_tmp
               k(Nk)=k_tmp
            else
c     Debug
c               write(*,*) 'nu_start=',nu_start
c     Debug
               nu(Nk)=nu_start
               call linear_interpolation(nu_prev,nu_tmp,k_prev,k_tmp,
     &              nu(Nk),k(Nk))
               call check_k(k(Nk),epsilon_k)
            endif
c     Debug
c            write(*,*) 'nu_start was identified:'
c            write(*,*) 'nu_prev=',nu_prev,' k_prev=',k_prev
c            write(*,*) 'nu_tmp=',nu_tmp,' k_tmp=',k_tmp
c            write(*,*) 'nu(',Nk,')=',nu(Nk),' k(',Nk,')=',k(Nk)
c            write(*,*) 'Nvalues_read=',Nvalues_read
c            stop
c     Debug
         endif

         if ((.not.nu_end_found).and.
     &        (nu_prev.lt.nu_end).and.
     &        (nu_tmp.ge.nu_end)) then
            nu_end_found=.true.
            Nk=Nk+1
            if (Nk.gt.kmax) then
               call error(label)
               write(*,*) 'Nk=',Nk
               write(*,*) 'while kmax=',kmax
               write(*,*) 'when nu_end has been found'
               stop
            endif
            nu(Nk)=nu_end
            call linear_interpolation(nu_prev,nu_tmp,k_prev,k_tmp,
     &           nu(Nk),k(Nk))
            call check_k(k(Nk),epsilon_k)
c     Debug
c            write(*,*) 'nu_end was identified:'
c            write(*,*) 'nu_prev=',nu_prev,' k_prev=',k_prev
c            write(*,*) 'nu_tmp=',nu_tmp,' k_tmp=',k_tmp
c            write(*,*) 'nu(',Nk,')=',nu(Nk),' k(',Nk,')=',k(Nk)
c     Debug
         endif

c     Debug
c         write(*,*) 'nu_start_found=',nu_start_found
c         write(*,*) 'nu_end_found=',nu_end_found
c     Debug

         if ((nu_start_found).and.(.not.nu_end_found)) then
            if (nu_tmp.gt.nu(Nk)) then
               Nk=Nk+1
               if (Nk.gt.kmax) then
                  call error(label)
                  write(*,*) 'Nk=',Nk
                  write(*,*) 'while kmax=',kmax
                  write(*,*) 'while reading band'
                  write(*,*) 'level:',level
                  write(*,*) 'band=',band
                  write(*,*) 'nu_start=',nu_start
                  write(*,*) 'nu_end=',nu_end
                  stop
               endif
               nu(Nk)=nu_tmp
               k(Nk)=k_tmp
            endif
         endif

         if (Nvalues_read.eq.Nvalues) then
            nu_end_found=.true.
         endif

         if (nu_start_found.and.nu_end_found) then ! stop reading the file
            band_done=.true.
         else ! store current value and read another value from file
c     Debug
c            write(*,*) 'Nvalues_read=',Nvalues_read,' Nvalues=',Nvalues
c     Debug
            if (Nvalues_read.lt.Nvalues) then
               nu_prev=nu_tmp
               k_prev=k_tmp
               read(13,52) nu_tmp,raf,k_tmp
               call check_k(k_tmp,epsilon_k)
               Nvalues_read=Nvalues_read+1
c     Debug
c            write(*,*) 'read: nu=',nu_tmp,' k=',k_tmp
c            if (Nvalues_read.eq.3) then
c               stop
c            endif
c     Debug
            endif
         endif

         if (Nvalues_read.eq.1) then
            if (.not.nu_start_found) then
               if ((nu_tmp.ge.nu_start).and.(nu_tmp.lt.nu_end)) then
                  nu_start_found=.true.
                  Nk=Nk+1
                  if (Nk.gt.kmax) then
                     call error(label)
                     write(*,*) 'Nk=',Nk
                     write(*,*) 'while kmax=',kmax
                     write(*,*) 'when nu_start has been found'
                     stop
                  endif
                  nu(Nk)=nu_tmp
                  k(Nk)=k_tmp
               endif
            endif ! nu_start_found=.false.
         endif ! Nvalues_read=1

      enddo ! while (band_done=.false.)

      if (Nk.le.0) then
         call error(label)
         write(*,*) 'Nk=',Nk
         stop
      else if (Nk.eq.1) then
         call error(label)
         write(*,*) 'band=',band
         write(*,*) 'nu_start=',nu_start,' nu_end=',nu_end
         write(*,*) 'Nk=',Nk
         write(*,*) 'Probable issue with k-file level:',level
         stop
      endif

      k_is_null=.true.
      do ik=1,Nk
         if (k(ik).gt.0.0D+0) then
            k_is_null=.false.
            goto 123
         endif
      enddo ! ik

 123  continue
c     double-check every k(ik) for negative values
      do ik=1,Nk
         call check_k(k(ik),epsilon_k)
      enddo ! ik
c     Debug
c      write(*,*) 'Nk=',Nk
c     Debug

 666  continue
      return
      end



      subroutine check_k(k,epsilon_k)
      implicit none
      include 'max.inc'
c     
c     Purpose: to check if the provided value of k is positive
c     and raise a flag if it is negative
c
c     Input:
c       + k: value of k
c       + epsilon_k: maximum (negative) value of k allowed
c
c     Output:
c       + k: can be modified
c
c     I/O
      double precision k,epsilon_k
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine check_k'

      if ((k.lt.0.0D+0).and.(dabs(k).le.epsilon_k)) then
         k=0.0D+0
      else if ((k.lt.0.0D+0).and.((dabs(k).gt.epsilon_k))) then
         call error(label)
         write(*,*) 'k=',k
         stop
      endif

      return
      end