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 precalc(band,sens_P,sens_T,sens_x,
     &     Nlnc,Nmol,i,ncl_indexes,mol_index,iso_index,
     &     mol_niso,Rgp,P,Ps,T,Tref,iso_abundance,
     &     iso_g,iso_mass,c2,
     &     sorted_index,nuc,dnuc_dP,
     &     gamma_L,dgammaL_dP,dgammaL_dT,dgammaL_dx,
     &     gamma_D,dgammaD_dT,
     &     S,dS_dT,
     &     densities,drho_dP,drho_dT,drho_dx)
      implicit none
      include 'max.inc'
      include 'formats.inc'
      include 'parameters.inc'
      include 'arrays1.inc'
c     
c     Purpose: to compute all line parameters that are needed prior to the computation of k(nu)
c
c     Inputs
      integer band
      logical sens_P,sens_T,sens_x
      integer Nlnc
      integer Nmol
      integer i
      integer ncl_indexes(1:Nline_mx)
      integer mol_index(1:Nmol_max)
      integer iso_index(1:Nmol_max,1:Niso_max)
      integer mol_niso(1:Nmol_max)
      double precision Rgp
      double precision P(1:Nmax)
      double precision Ps(1:Nmol_max)
      double precision T(1:Nmax)
      double precision Tref(1:Nmol_max)
      double precision iso_abundance(1:Nmol_max,1:Niso_max)
      integer iso_g(1:Nmol_max,1:Niso_max)
      double precision iso_mass(1:Nmol_max,1:Niso_max)
      double precision c2
c     Outputs
      integer sorted_index(1:Nline_mx)
      double precision nuc(1:Nline_mx)
      double precision dnuc_dP(1:Nline_mx)
      double precision gamma_L(1:Nline_mx)
      double precision dgammaL_dP(1:Nline_mx)
      double precision dgammaL_dT(1:Nline_mx)
      double precision dgammaL_dx(1:Nline_mx)
      double precision gamma_D(1:Nline_mx)
      double precision dgammaD_dT(1:Nline_mx)
      double precision S(1:Nline_mx)
      double precision dS_dT(1:Nline_mx)
      double precision densities(1:Nline_mx)
      double precision drho_dP(1:Nline_mx)
      double precision drho_dT(1:Nline_mx)
      double precision drho_dx(1:Nline_mx)
c     temp
      integer molec,isotop
      integer l,line,isof
      double precision Q,Qref,dQ_dT
      integer strlen
      character*(Nchar_mx) label
      label='subroutine precalc'

      do l=1,Nlnc
         line=ncl_indexes(l)
         call indexes(Nmol,mol_index,iso_index,mol_niso,
     &        index(line),isotope(line),molec,isotop,isof)
c     Debug
         if (isof.eq.0) then
            Q=0.0D+0
            dQ_dT=0.0D+0
         else
c     Debug
            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
               call Q_sensitivity(index(line),isotope(line),
     &              iso_g(molec,isotop),T(i),Q,dQ_dT)
            else
               dQ_dT=0.0D+0
            endif
c     Debug
         endif
c     Debug
         nuc(l)=nu0(line)+delta_air(line)*P(i)
         dnuc_dP(l)=delta_air(line)
c     for some reason, the BD_TIPS_2003 routine can not compute
c     partition function Q for species index 34 (Q=0. is imposed).
         if (Q.eq.0.0D+0) then  ! line intensity can not be computed
            gamma_L(l)=1.0D+0
            gamma_D(l)=1.0D+0
            S(l)=0.0D+0
            dgammaL_dP(l)=0.0D+0
            dgammaL_dT(l)=0.0D+0
            dgammaL_dx(l)=0.0D+0
            dgammaD_dT(l)=0.0D+0
            dS_dT(l)=0.0D+0
         else
c     Lorentz half-width [cm¯¹]
            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(l))
c     Debug
c            write(*,*) nuc(l),Ps(molec),gamma_L(l)
c     Debug
c     Doppler half-width [cm¯¹]
            call Doppler_linewidth(T(i),iso_mass(molec,isotop),
     &           nu0(line),gamma_D(l))
c     line intensity [cm¯¹/(molecule.cm¯²)]
            call line_intensity(Tref(index(line)),T(i),nu0(line),
     &           Sref(line),Elow(line),c2,Q,Qref,S(l))
c     Debug
c            write(*,*) 'S(',l,')=',S(l)
c            if (index(line).eq.2) then
c               write(61,*) line,nu0(line),gamma_L(l),gamma_D(l),S(l)
c            endif
c     Debug
c     ----------------------------------
c     sensitivities
            if ((sens_P).or.(sens_T).or.(sens_x)) then
               call gammaL_sensitivities(Tref(index(line)),
     &              gamma_air(line),gamma_self(line),
     &              n_air(line),n_self(line),
     &              P(i),T(i),Ps(molec),gamma_L(l),
     &              dgammaL_dP(l),dgammaL_dT(l),dgammaL_dx(l))
               call gammaD_sensitivity(iso_mass(molec,isotop),
     &              nu0(line),T(i),gamma_D(l),dgammaD_dT(l))
               call S_sensitivity(Tref(index(line)),T(i),nu0(line),
     &              Sref(line),Elow(line),c2,Q,dQ_dT,Qref,S(l),
     &              dS_dT(l))
            else
               dgammaL_dP(l)=0.0D+0
               dgammaL_dT(l)=0.0D+0
               dgammaL_dx(l)=0.0D+0
               dgammaD_dT(l)=0.0D+0
               dS_dT(l)=0.0D+0
            endif
c     ----------------------------------
         endif
         call density(Rgp,P(i),Ps(molec),T(i),
     &        iso_abundance(molec,isotop),
     &        densities(l),drho_dP(l),drho_dT(l),drho_dx(l))

c     Debug
c         if (S(l).gt.Smax) then
c            Smax=S(l)
c            lmax=l
c         endif
c         if ((l.eq.79182).or.(l.eq.79221)) then
c            write(*,*) Sref(line)
c            write(*,*) Qref,Q,Qref/Q
c            write(*,*) c2,Elow(line),T(i),dexp(-c2*Elow(line)/T(i))
c            write(*,*) c2,Elow(line),Tref(index(line)),
c     &           dexp(-c2*Elow(line)/Tref(index(line)))
c            write(*,*) 1.0D+0-dexp(-c2*nu0(line)/T(i))
c            write(*,*) 1.0D+0-dexp(-c2*nu0(line)/Tref(index(line)))
c            write(*,*) Sref(line)*Qref/Q*
c     &           dexp(-c2*Elow(line)/T(i))/
c     &           dexp(-c2*Elow(line)/Tref(index(line)))*
c     &           (1.0D+0-dexp(-c2*nu0(line)/T(i)))/
c     &           (1.0D+0-dexp(-c2*nu0(line)/Tref(index(line))))
c            write(*,*) 'S(',l,')=',S(l),nu0(line),nuc(l)

c            write(*,*) Tref(index(line)),T(i),n_air(line),
c     &           (Tref(index(line))/T(i))**n_air(line)
c            write(*,*) gamma_air(line),P(i),Ps(molec),
c     &           gamma_air(line)*(P(i)-Ps(molec))
c            write(*,*) gamma_self(line),Ps(molec),
c     &           gamma_self(line)*Ps(molec)
c            write(*,*) (Tref(index(line))/T(i))**n_air(line)*(
c     &           gamma_air(line)*(P(i)-Ps(molec))+
c     &           gamma_self(line)*Ps(molec))
c            write(*,*) 'gamma_L(',l,')=',gamma_L(l),nu0(line),nuc(l)
 
c            write(*,*) T(i),iso_mass(molec,isotop)
c            write(*,*) nu0(line)*cdop*
c     &           dsqrt(T(i)/iso_mass(molec,isotop)*1.0D+3)
c           write(*,*) 'gamma_D(',l,')=',gamma_D(l),nu0(line),nuc(l)
c
c            write(*,*) 'densities(',l,')=',densities(l)
c         endif
c     Debug
      enddo ! l

c     Debug
c      write(*,*) 'Smax=',Smax,' nuc=',nuc(lmax)
c      do l=1,Nlnc
c         line=ncl_indexes(l)
c         if (S(l).eq.Smax) then
c            write(*,*) 'line=',line,' nu0=',nu0(line),' nuc=',nuc(l)
c         endif
c      enddo 
c     Debug

c     Debug
c      close(61)
c     Debug

c     Sortng output arrays by descending order of line intensity
c      do l=1,Nlnc
c         sorted_index(l)=index(l)
c      enddo ! l
      call sort_precalc_data(Nlnc,index,
     &     sorted_index,nuc,dnuc_dP,
     &     gamma_L,dgammaL_dP,dgammaL_dT,dgammaL_dx,
     &     gamma_D,dgammaD_dT,
     &     S,dS_dT,densities,
     &     drho_dP,drho_dT,drho_dx)

      return
      end