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 khi_profile(nu,nuc,T,index,
     &     usl,sl_choice,slwr,slam,slas,khi)
      implicit none
      include 'max.inc'
      double precision nu,nuc,T
      integer index
      logical usl
      integer sl_choice
      logical slwr,slam,slas
      double precision khi

      khi=1.0D+0
      if (usl) then
c     if (.not.slam), SL profiles have to be used only for CO2
c     if (slam), SL profiles have to be used for all molecules
         if (((.not.slam).and.(index.eq.2)).or.(slam)) then ! CO2
            call khi_co2(nu,nuc,T,sl_choice,slwr,slas,khi)
         endif
      endif

      return
      end



      subroutine khi_co2(nu,nuc,T,sl_choice,slwr,slas,khi)
      implicit none
      include 'max.inc'
c
c     Purpose: to compute factor KHI that accounts for the sub-Lorentzian
c     nature of line profile.
c     This routine is based on results from:
c     + Pollack et al., Icarus 1993 (vol. 103, pp. 1-42)
c     + Perrin & Hartmann, JQSRT 1989 (vol. 42, No. 4, pp. 311-317)
c     "Temperature-dependant measurements and modeling of absorption by CO2-N2 mixtures in the far
c     line-wings of the 4.3µm CO2 band"
c     + Tonkov at al. Applied Optics 1996 (vol. 35, No. 24, pp. 4863-4870):
c     "Measurements and empirical modeling of pure CO2 absorption in the 2.3 µm region at room temperature:
c     far wings, allowed and collision-induced bands"
c
c     It will compute khi for pure CO2 lines, in both 2.3 and 4.3µm bands.
c
c     Inputs:
c       + nu: wavenumber (cm-¹) where k(nu) has to be computed
c       + nuc: line center wavenumber (cm-¹)
c       + T: temperature (K)
c       + sl_choice: choice of sub-lorentzian profile
c         - 0: use whatever is relevant for current wavenumber, within validity ranges
c         - 1: use Pollack results @ 1.18µm and 2.7µm
c         - 2: use Perrin @ Hartmann @ 4.3µm
c         - 3: use Tonkov @ 2.3µm
c       + slwr: option for using CO2 SL profiles over the whole IR range
c       + slas: option to select using asymetric profile
c
c     Outputs:
c       + khi: khi factor (no dimension)
c     
      double precision nu,nuc,T
      integer sl_choice
      logical slwr,slas
      double precision khi
      double precision delta_nu
      double precision alpha(1:3),beta(1:3),epsilon(1:3),B(1:3)
      double precision sigma(1:3),zeta(1:3)
      double precision asf
      integer i
c     Debug
      integer out
c     Debug
      logical valid
      integer strlen
      character*(Nchar_mx) label
      label='subroutine ksi_co2'

      delta_nu=dabs(nu-nuc)
      khi=1.0D+0 ! default

      if (sl_choice.eq.0) then
c     Use whatever is relevant at the local wavenumber, within validity ranges
         if ((nuc.ge.7400.0D+0).and.(nuc.le.10000.0D+0)) then ! 1.18 µm band
            call khi_co2_1micron18(delta_nu,T,khi)
         else if ((nuc.ge.3800.0D+0).and.(nuc.le.4700.0D+0)) then ! 2.3µm band
            call khi_co2_2micron30(delta_nu,T,khi)
         else if ((nuc.ge.3700.0D+0).and.(nuc.le.4100.0D+0)) then ! 2.7 µm band
            call khi_co2_2micron70(delta_nu,T,khi)
         else if ((nuc.ge.2100.0D+0).and.(nuc.le.2600.0D+0)) then ! 4.3µm band
            call khi_co2_4micron30(delta_nu,T,khi)
         endif
c     
      else if (sl_choice.eq.1) then
c     Results from Pollack et al. 1993, (derived at room temperature)
         if (slwr) then
c     no respect of validity range: use 1.18 microns band
            call khi_co2_1micron18(delta_nu,T,khi)
         else                   ! inside validity range
            if ((nuc.ge.7400.0D+0).and.(nuc.le.10000.0D+0)) then ! 1.18 µm band
               call khi_co2_1micron18(delta_nu,T,khi)
            else if ((nuc.ge.3700.0D+0).and.(nuc.le.4100.0D+0)) then ! 2.7 µm band
               call khi_co2_2micron70(delta_nu,T,khi)
            endif
         endif                  ! slwr
c     
      else if (sl_choice.eq.2) then
c     Results from Perrin & Hartmann 1989 valid for 300<T<800 K
         if ((slwr).or.((.not.slwr).and.(nuc.ge.2100.0D+0)
     &        .and.(nuc.le.2600.0D+0))) then ! 4.3µm band
            call khi_co2_4micron30(delta_nu,T,khi)
         endif
c     
      else if (sl_choice.eq.3) then
c     Results from Tonkov et al. valid @ T=296K
         if ((slwr).or.((.not.slwr).and.(nuc.ge.3800.0D+0)
     &        .and.(nuc.le.4700.0D+0))) then ! 2.3µm band
            call khi_co2_2micron30(delta_nu,T,khi)
         endif
      else                      ! value of sl_choice
         call error(label)
         write(*,*) 'sl_choice=',sl_choice
         stop
      endif                     ! sl_choice

      if (slas) then
         call asymetry(nu,nuc,asf)
         khi=khi*asf
      endif                     ! slas

      return
      end



      subroutine khi_co2_1micron18(delta_nu,T,khi)
      implicit none
      include 'max.inc'
c
c     Purpose: to compute khi factor based on results from Pollack et al. in the 1.18 micron band of CO2
c
c     Input:
c       + delta_nu: distance from line center (inv. cm)
c       + T: temperature (K)
c
c     Output:
c       khi: khi factor
c
c     I/O
      double precision delta_nu
      double precision T
      double precision khi
c     temp
      double precision alpha(1:3)
      double precision beta(1:3)
      double precision sigma(1:3)
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine khi_co2_1micron18'

      alpha(1)=1.350D+0
      alpha(2)=0.6140D+0
      beta(1)=1.0D+1
      beta(2)=9.00D+1
      sigma(1)=3.0D+0           ! cm-¹
      sigma(2)=1.0D+1           ! cm-¹
      sigma(3)=1.60D+2          ! cm-¹

      if ((delta_nu.ge.0.0D+0).and.
     &     (delta_nu.lt.sigma(1))) then
         khi=1.0D+0
      else if ((delta_nu.ge.sigma(1)).and.
     &        (delta_nu.lt.sigma(2))) then
         khi=alpha(1)*dexp(-delta_nu/beta(1))
      else if ((delta_nu.ge.sigma(2)).and.
     &        (delta_nu.lt.sigma(3))) then
         khi=alpha(2)*dexp(-delta_nu/beta(2))
      else if (delta_nu.ge.sigma(3)) then
         khi=0.0D+0
      endif

      return
      end



      subroutine khi_co2_2micron30(delta_nu,T,khi)
      implicit none
      include 'max.inc'
c
c     Purpose: to compute khi factor based on results from Tonkov et al. in the 2.30 micron band of CO2
c
c     Input:
c       + delta_nu: distance from line center (inv. cm)
c       + T: temperature (K)
c
c     Output:
c       khi: khi factor
c
c     I/O
      double precision delta_nu
      double precision T
      double precision khi
c     temp
      double precision sigma(1:3)
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine khi_co2_2micron30'

      sigma(1)=3.0D+0            ! cm-¹
      sigma(2)=1.5D+2            ! cm-¹
      sigma(3)=3.0D+2            ! cm-¹

      if ((delta_nu.ge.0.0D+0).and.
     &     (delta_nu.lt.sigma(1))) then
         khi=1.0D+0
      else if ((delta_nu.ge.sigma(1)).and.
     &        (delta_nu.lt.sigma(2))) then
         khi=1.084D+0*dexp(-2.7D-2*delta_nu)
      else if ((delta_nu.ge.sigma(2)).and.
     &        (delta_nu.lt.sigma(3))) then
         khi=2.08D-1*dexp(-1.60D-2*delta_nu)
      else if (delta_nu.gt.sigma(3)) then
         khi=2.5D-2*dexp(-9.0D-3*delta_nu)
      endif

      return
      end



      subroutine khi_co2_2micron70(delta_nu,T,khi)
      implicit none
      include 'max.inc'
c
c     Purpose: to compute khi factor based on results from Pollack et al. in the 2.70 micron band of CO2
c
c     Input:
c       + delta_nu: distance from line center (inv. cm)
c       + T: temperature (K)
c
c     Output:
c       khi: khi factor
c
c     I/O
      double precision delta_nu
      double precision T
      double precision khi
c     temp
      double precision alpha(1:3)
      double precision beta(1:3)
      double precision sigma(1:3)
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine khi_co2_2micron70'

      alpha(1)=1.350D+0
      alpha(2)=0.6140D+0
      beta(1)=1.0D+1
      beta(2)=4.70D+1
      sigma(1)=3.0D+0           ! cm-¹
      sigma(2)=1.0D+1           ! cm-¹
      sigma(3)=1.20D+2          ! cm-¹

      if ((delta_nu.ge.0.0D+0).and.
     &     (delta_nu.lt.sigma(1))) then
         khi=1.0D+0
      else if ((delta_nu.ge.sigma(1)).and.
     &        (delta_nu.lt.sigma(2))) then
         khi=alpha(1)*dexp(-delta_nu/beta(1))
      else if ((delta_nu.ge.sigma(2)).and.
     &        (delta_nu.lt.sigma(3))) then
         khi=alpha(2)*dexp(-delta_nu/beta(2))
      else if (delta_nu.ge.sigma(3)) then
         khi=0.0D+0
      endif

      return
      end



      subroutine khi_co2_4micron30(delta_nu,T,khi)
      implicit none
      include 'max.inc'
c
c     Purpose: to compute khi factor based on results from Perrin & Hartmann in the 4.30 micron band of CO2
c
c     Input:
c       + delta_nu: distance from line center (inv. cm)
c       + T: temperature (K)
c
c     Output:
c       khi: khi factor
c
c     I/O
      double precision delta_nu
      double precision T
      double precision khi
c     temp
      double precision alpha(1:3)
      double precision beta(1:3)
      double precision epsilon(1:3)
      double precision sigma(1:3)
      double precision B(1:3)
      integer i
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine khi_co2_4micron30'

      alpha(1)=8.88D-2
      alpha(2)=0.0D+0
      alpha(3)=2.32D-2
      beta(1)=-1.60D-1
      beta(2)=5.26D-2
      beta(3)=0.0D+0
      epsilon(1)=4.10D-3
      epsilon(2)=1.52D-3
      epsilon(3)=0.0D+0
      sigma(1)=3.0D+0           ! cm-¹
      sigma(2)=3.0D+1           ! cm-¹
      sigma(3)=1.2D+2           ! cm-¹
      do i=1,3
         B(i)=alpha(i)+beta(i)*dexp(-epsilon(i)*T)
      enddo ! i

      if ((delta_nu.ge.0.0D+0).and.
     &     (delta_nu.lt.sigma(1))) then
         khi=1.0D+0
      else if ((delta_nu.ge.sigma(1)).and.
     &        (delta_nu.lt.sigma(2))) then
         khi=dexp(-B(1)*(delta_nu-sigma(1)))
      else if ((delta_nu.ge.sigma(2)).and.
     &        (delta_nu.lt.sigma(3))) then
         khi=dexp(-B(1)*(sigma(2)-sigma(1))
     &        -B(2)*(delta_nu-sigma(2)))
      else if (delta_nu.ge.sigma(3)) then
         khi=dexp(-B(1)*(sigma(2)-sigma(1))
     &        -B(2)*(sigma(3)-sigma(2))
     &        -B(3)*(delta_nu-sigma(3)))
      endif

      return
      end


      
      subroutine asymetry(nu,nuc,asf)
      implicit none
      include 'max.inc'
c     
c     Purpose: to compute asymetry factor of sub-lorentzian profile khi(nu)
c     
c     Inputs:
c       + nu: wavenumber (cm-¹) at which asymetry factor has to be computed
c       + nuc: center wavenumber (cm-¹) of the line
c
c     Outputs:
c       + asf: asymetry factor
c
c     Asymetry factors are taken from Menoux, LeDoucen et Boulet, 1987:
c     "Line shape in the low-frequency wing of self-broadened CO2 lines"
c     Appled Optics, Vol. 26, No. 3, pp. 554-562
c
      double precision nu
      double precision nuc
      double precision asf
      double precision f_asymetry
      double precision sigma,sigma0
      double precision alpha1,alpha2
      double precision beta1,beta2
      double precision f1,f2
c     Debug
      integer out
c     Debug
      integer strlen
      character*(Nchar_mx) label
      label='subroutine asymetry'

      alpha1=0.889D+0
      alpha2=0.220D+0
      beta1=31.627D+0
      beta2=50.063D+0

      sigma0=1.20D+2 ! cm-¹
      sigma=dabs(nu-nuc)
c     for distances to nuc lower that sigma0, there is no asymetry
      if (sigma.le.sigma0) then
         asf=1.0D+0
c     for distances to nuc greater than sigma0, asymetry factor has to be computed
      else
         f1=f_asymetry(alpha1,beta1,sigma)
         f2=f_asymetry(alpha2,beta2,sigma)
         if (f1+f2.eq.0.0D+0) then
            asf=1.0D+0
         else
            if (nu.gt.nuc) then
               asf=2.0D+0*f1/(f1+f2)
            else
               asf=2.0D+0*f2/(f1+f2)
            endif
         endif
      endif

      return
      end




      double precision function f_asymetry(alpha,beta,sigma)
      implicit none
      include 'max.inc'
      double precision alpha,beta,sigma
      integer strlen
      character*(Nchar_mx) label
      label='subroutine f_asymetry'

      f_asymetry=alpha*dexp(-sigma/beta)

      return
      end