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 line_intensity(Tref,T,nu0,Sref,Elow,c2,Q,Qref,S)
      implicit none
      include 'max.inc'
c
c     Purpose: to compute the intensity of the given spectral line
c
c     Inputs:
c       + Tref: reference temperature
c       + T: temperature [K]
c       + nu0: vacuum wavenumber [cm¯¹]
c       + Sref: line intensity [cm¯¹/molecule.cm¯²] at standard Tref
c       + Elow: lower-state energy [cm¯¹]
c       + c2:  second radiation constant [K.cm]
c       + Q: total internal partition sum
c       + Qref: total internal partition sum at reference temperature Tref
c
c     Outputs:
c       + S: line intensity [cm¯¹/molecule.cm¯²]
c
      include 'parameters.inc'
c     I/O
      double precision Tref
      double precision T
      double precision nu0
      double precision Sref
      double precision Elow
      double precision c2
      double precision Q
      double precision Qref
      double precision S
c     temp
      integer strlen
      character*(Nchar_mx) label
c     Debug
      integer out1,out2
c     Debug
      label='subroutine line_intensity'

      S=Sref
     &     *Qref/Q
     &     *dexp(-c2*Elow/T)/dexp(-c2*Elow/Tref)
     &     *(1.0D+0-dexp(-c2*nu0/T))/(1.0D+0-dexp(-c2*nu0/Tref)) ! [cm¯¹/(molecule.cm¯²)]
c     Debug
c$$$      write(*,*) 'This is routine line_intensity'
c$$$      write(*,*) 'S=',S
c$$$      write(*,*) 'Tref=',Tref
c$$$      write(*,*) 'T=',T
c$$$      write(*,*) 'nu0=',nu0
c$$$      write(*,*) 'Sref=',Sref
c$$$      write(*,*) 'Elow=',Elow
c$$$      write(*,*) 'c2=',c2
c$$$      write(*,*) 'Q=',Q
c$$$      write(*,*) 'Qref=',Qref
c$$$      write(*,*) dexp(-c2*Elow/T)
c$$$      write(*,*) dexp(-c2*Elow/Tref)
c$$$      write(*,*) (1.0D+0-dexp(-c2*nu0/T))
c$$$      write(*,*) (1.0D+0-dexp(-c2*nu0/Tref))
c$$$      stop
c     Debug
c     Debug
      call test_nan(S,out1)
      call test_inf(S,out2)
      if ((out1.eq.1).or.(out2.eq.1)) then
         write(*,*) 'This is routine line_intensity'
         write(*,*) 'S=',S
         write(*,*) 'Tref=',Tref
         write(*,*) 'T=',T
         write(*,*) 'nu0=',nu0
         write(*,*) 'Sref=',Sref
         write(*,*) 'Elow=',Elow
         write(*,*) 'c2=',c2
         write(*,*) 'Q=',Q
         write(*,*) 'Qref=',Qref
         write(*,*) dexp(-c2*Elow/T)
         write(*,*) dexp(-c2*Elow/Tref)
         write(*,*) (1.0D+0-dexp(-c2*nu0/T))
         write(*,*) (1.0D+0-dexp(-c2*nu0/Tref))
         stop
      endif
c     Debug

      return
      end