c     kspectrum (http://www.meso-star.com/en_Products.html) - This file is part of kspectrum
c     Copyright (C) 2008-2015 - Méso-Star - Vincent Eymet
c     
c     This file must be used under the terms of the CeCILL license.
c     This source file is licensed as described in the file COPYING, which
c     you should have received as part of this distribution.  The terms
c     are also available at
c     http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
c     
      subroutine constantk(nu_min,nu_max,Nlines,Tref,T,P,Ps,i,
     &     Nmol,mol_index,mol_niso,iso_index,
     &     iso_abundance,iso_mass,
     &     c2,Rgp,iso_g,k_mxep,
     &     line_profile,usl,sl_choice,slwr,slam,slas,
     &     sens_T,sens_P,sens_x,
     &     ltr,dorg,tr_dist,ngamma,
     &     c,sigmac,kc,
     &     dsigmac_dP,dsigmac_dT,dsigmac_dx,
     &     dkc_dP,dkc_dT,dkc_dx,
     &     Nlc,Nlnc,ncl_indexes)
      implicit none
      include 'max.inc'
      include 'parameters.inc'
      include 'formats.inc'
      include 'arrays1.inc'
c
c     Purpose: determine which lines can be considered as constant
c     over the specified narrowband index
c
c     Inputs:
c       + nu_min: minimum value of nu for the given narrowband interval
c       + nu_max: maximum value of nu for the given narrowband interval
c       + Nlines: number of spectral lines in the HITRAN database
c       + Tref: array of reference temperatures
c       + T: temperature, at each atmospheric level
c       + P: the total pressure, at each atmospheric level
c       + Ps: partial pressure of each species, at each level
c       + m : number of atmospheric levels
c       + i: index of the current atmospheric level
c       + index: number of the HITRAN molecule
c       + isotope: number of the HITRAN isotope
c       + Nmol: number of molecules used in the HITRAN database
c       + mol_index: indexes of the HITRAN molecules
c       + mol_niso: number of isotopes for each molecule
c       + iso_index: indexes of the HITRAN isotopes
c       + iso_abundance: abundances of isotopes for each molecule
c       + iso_mass(mol,iso) are molar masses [g/mol]
c       + nu0: vacuum wavenumber [cm¯¹]
c       + Sref: line intensities [cm¯¹/molecule.cm¯²] at standard Tref
c       + gamma_air: air-broadened half-width [cm¯¹.atm¯¹]
c       + gamma_self: self-broadened half-width [cm¯¹.atm¯¹]
c       + Elow: lower-state energy [cm¯¹]
c       + n_air: temperature-dependance expnent for gamma_air [unitless]
c       + delta_air: air pressure-induced line shift [cm¯¹.atm¯¹]
c       + n_self: temperature-dependance expnent for gamma_self [unitless]
c       + c2:  second radiation constant [K.cm]
c       + Rgp: perfect gas constant [atm.m³/(mol.K)]
c       + iso_Qref: total internal partition sum at reference temperature (Tref)
c       + iso_g(mol,iso) are state independent degeneracy factors
c       + k_mxep: maximum allowed error percentage of k
c       + line_profile: 1 for a Lorentz profile, 2 for a Voigt profile
c       + usl: options for using sub-lorentzian profiles
c       + sl_choice: choice of sub-lorentzian profile
c       + slwr: option for using CO2 SL profiles over the whole IR range
c       + slam: option for using CO2 SL profiles for all molecules
c       + slas: option for taking into account SL profile asymetry
c       + sens_T: option for computing sensitivities to temperature
c       + sens_P: option for computing sensitivities to pressure
c       + sens_x: option for computing sensitivities to concentrations
c       + ltr: true if line truncation is required
c       + dorg: (1) use constant truncation distance or (2) multiple of line width
c       + tr_dist: truncation distance when truncation is required
c       + ngamma: multiple of line width for truncation
c
c     Output:
c       + c: 0 if line does not have a contant value; yes otherwise
c       + sigmac: contant value of the absorption cross-section if c=1 [cm²/molecule]
c       + kc: contant value of k over the narrowband, if c=1
c       + dsigmac_dP: sensitivity of sigmac to total pressure [cm²/molecule/atm]
c       + dsigmac_dT: sensitivity of sigmac to temperature [cm²/molecule/K]
c       + dsigmac_dx: sensitivity of sigmac to the concentration of chemical species (array) [cm²/molecule]
c       + dkc_dP: sensitivity of kc to total pressure [m¯¹/atm] 
c       + dkc_dT: sensitivity of kc to temperature [m¯¹/K]
c       + dkc_dx: sensitivity of kc to the concentration of chemical species (array) [m¯¹]
c       + Nlc: number of lines that have a constant value
c       + Nlnc: number of lines that do not have a constant value
c       + ncl_indexes: indexes of lines that do not have a constant value
c
c     Inputs
      double precision nu_min,nu_max
      integer Nlines
      double precision Tref(1:Nmol_max)
      double precision T(1:Nmax)
      double precision P(1:Nmax)
      double precision Ps(1:Nmol_max)
c      integer m
      integer i
      double precision c2,Rgp
c      double precision iso_Qref(1:Nmol_max,1:Niso_max)
      integer iso_g(1:Nmol_max,1:Niso_max)
      integer Nmol
      integer mol_index(1:Nmol_max)
      integer mol_niso(1:Nmol_max)
      integer iso_index(1:Nmol_max,1:Niso_max)
      double precision iso_abundance(1:Nmol_max,1:Niso_max)
      double precision iso_mass(1:Nmol_max,1:Niso_max)
      double precision k_mxep
      integer line_profile
      logical usl
      integer sl_choice
      logical slwr
      logical slam
      logical slas
      logical sens_T,sens_P,sens_x
      logical ltr
      integer dorg
      double precision tr_dist
      integer ngamma
c     Outputs
      integer c(1:Nline_mx),Nlc,Nlnc
      integer ncl_indexes(1:Nline_mx)
      double precision sigmac,kc
      double precision dsigmac_dP
      double precision dsigmac_dT
      double precision dsigmac_dx(1:Nmol_max)
      double precision dkc_dP
      double precision dkc_dT
      double precision dkc_dx(1:Nmol_max)
c     temp
      integer line,molec,isotop,isof,mol
      double precision Q,Qref,dQ_dT
      double precision dens,drho_dP,drho_dT,drho_dx
      double precision gamma_L
      double precision sigma1,sigma2,nuc
      double precision k1,k2,kmin
      double precision dsigma1_dP,dsigma1_dT,dsigma1_dx
      double precision dk1_dP,dk1_dT,dk1_dx
      double precision dsigma2_dP,dsigma2_dT,dsigma2_dx
      double precision dk2_dP,dk2_dT,dk2_dx
      double precision eps
      logical truncate
      integer strlen
      character*(Nchar_mx) label
      label='subroutine constantk'

c     Debug
c      open(60,file='./constant_lines.txt')
c     Debug
      Nlc=0
      Nlnc=0
      sigmac=0.0D+0
      kc=0.0D+0
      dsigmac_dP=0.0D+0
      dsigmac_dT=0.0D+0
      do mol=1,Nmol
         dsigmac_dx(mol)=0.0D+0
      enddo
      dkc_dP=0.0D+0
      dkc_dT=0.0D+0
      do mol=1,Nmol
         dkc_dx(mol)=0.0D+0
      enddo
      write(*,40) 'Line classification for spectral interval [',
     &     nu_min,'-',nu_max,'] cm^{-1}...'
c     Debug
      if (Nlines.gt.Nline_mx) then
         call error(label)
         write(*,*) 'Nlines=',Nlines
         write(*,*) 'while Nline_mx=',Nline_mx
         stop
      endif
c     Debug
      do line=1,Nlines
c     Debug
c         write(*,*) 'line=',line,' / ',Nlines
c     Debug
         nuc=nu0(line)+delta_air(line)*P(i)
c     lines whose center is within [nu_min,nu_max] are not examined
         if ((nuc.gt.nu_min).and.(nuc.lt.nu_max)) then
            c(line)=0
            Nlnc=Nlnc+1
            ncl_indexes(Nlnc)=line
            goto 123
         endif
c     get "molec" and "isotop", indexes (in arrays iso_*) of molecular species
c     described by index "line"
c     Debug
         if (index(line).eq.0) then
            call error(label)
            write(*,*) 'index(',line,')=',index(line)
            stop
         endif
c     Debug
         call indexes(Nmol,mol_index,iso_index,mol_niso,
     &        index(line),isotope(line),molec,isotop,isof)
c     Debug
         if (isof.eq.0) then
            write(*,*) 'error from routine constantk:'
            write(*,*) 'isof=',isof
            stop
         endif
c     Debug
         call density(Rgp,P(i),Ps(molec),T(i),
     &        iso_abundance(molec,isotop),
     &        dens,drho_dP,drho_dT,drho_dx)
         call Lorentz_linewidth(Tref(index(line)),T(i),P(i),Ps(molec),
     &        n_air(line),n_self(line),gamma_air(line),gamma_self(line),
     &        gamma_L)
c     2013-10-15: line truncation has been displaced in this routine, in order
c     to provide only the indexes of lines whose contribution should be considered
         call line_truncation(ltr,dorg,tr_dist,ngamma,
     &        nu_min,nu_max,nuc,gamma_L,usl,index(line),
     &        truncate)
         if (truncate) then
            c(line)=1
            Nlc=Nlc+1
            goto 123
         endif
c     2013-10-15
         call TIPS(index(line),Tref(index(line)),isotope(line),
     &        iso_g(molec,isotop),Qref)
         call TIPS(index(line),T(i),isotope(line),
     &        iso_g(molec,isotop),Q)
         if (sens_T) then
c     Debug
c            write(*,*) 'ok0'
c     Debug
            call Q_sensitivity(index(line),isotope(line),
     &           iso_g(molec,isotop),T(i),Q,dQ_dT)
c     Debug
c            write(*,*) 'ok1'
c     Debug
         else
            dQ_dT=0.0D+0
         endif
         if (Q.eq.0.0D+0) then
            goto 123
         endif
c     Debug
         if (nu0(line).lt.1.0D-10) then
            write(*,*) 'nu0(',line,')=',nu0(line)
            stop
         endif
c     Debug
         call calc_k(nu_min,Tref(index(line)),T(i),index(line),
     &        P(i),Ps(molec),
     &        nu0(line),iso_mass(molec,isotop),Sref(line),
     &        gamma_air(line),gamma_self(line),Elow(line),
     &        n_air(line),delta_air(line),n_self(line),
     &        c2,dens,drho_dP,drho_dT,drho_dx,Q,Qref,dQ_dT,
     &        line_profile,usl,sl_choice,slwr,slam,slas,
     &        sens_T,sens_P,sens_x,
     &        sigma1,k1,
     &        dsigma1_dP,dsigma1_dT,dsigma1_dx,dk1_dP,dk1_dT,dk1_dx)
         call calc_k(nu_max,Tref(index(line)),T(i),index(line),
     &        P(i),Ps(molec),
     &        nu0(line),iso_mass(molec,isotop),Sref(line),
     &        gamma_air(line),gamma_self(line),Elow(line),
     &        n_air(line),delta_air(line),n_self(line),
     &        c2,dens,drho_dP,drho_dT,drho_dx,Q,Qref,dQ_dT,
     &        line_profile,usl,sl_choice,slwr,slam,slas,
     &        sens_T,sens_P,sens_x,
     &        sigma2,k2,
     &        dsigma2_dP,dsigma2_dT,dsigma2_dx,dk2_dP,dk2_dT,dk2_dx)
         kmin=k1
         if (k2.lt.kmin) then
            kmin=k2
         endif
         eps=dabs(k1-k2)/kmin
         if (eps.lt.k_mxep) then
            c(line)=1
            Nlc=Nlc+1
            sigmac=sigmac+(sigma1+sigma2)/2.0D+0
            kc=kc+(k1+k2)/2.0D+0
c            kc=kc+k2
c     Debug
c            if (index(line).eq.2) then
c               write(60,*) line,nu0(line),(k1+k2)/2.0D+0
c            endif
c     Debug
            if (sens_P) then
               dsigmac_dP=dsigmac_dP+(dsigma1_dP+dsigma2_dP)/2.0D+0
               dkc_dP=dkc_dP+(dk1_dP+dk2_dP)/2.0D+0
            endif
            if (sens_T) then
               dsigmac_dT=dsigmac_dT+(dsigma1_dT+dsigma2_dT)/2.0D+0
               dkc_dT=dkc_dT+(dk1_dT+dk2_dT)/2.0D+0
            endif
            if (sens_x) then
               dsigmac_dx(molec)=dsigmac_dx(molec)
     &              +(dsigma1_dx+dsigma2_dx)/2.0D+0
               dkc_dx(molec)=dkc_dx(molec)
     &              +(dk1_dx+dk2_dx)/2.0D+0
            endif
         else
            c(line)=0
            Nlnc=Nlnc+1
            ncl_indexes(Nlnc)=line
         endif
 123     continue
      enddo ! lines
c      call save_temp(Nlnc,ncl_indexes)
      write(*,*) '... done'
c     Debug
c      close(60)
c     Debug

      return
      end