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 CIA(T,nu,Nmol,P,Ps,mol_index,ucia,ciac,ciawr,ciatr,
     &     sens_P,sens_T,sens_x,
     &     kCIA,dkCIA_dP,dkCIA_dT,dkCIA_dxco2)
      implicit none
      include 'max.inc'
c     
c     Purpose: to compute collision-induced absorption (CIA)
c
c     Inputs:
c       + T: temperature (K)
c       + nu: wavenumber (cm-¹)
c       + Nmol: number of radiatively active species
c       + P: total pressure (atm)
c       + Ps: partial pressures (atm)
c       + mol_index: indexes of the Nmol radiatively active species
c       + ucia: option for using CIA for CO2
c       + ciac: option for CIA source
c       + ciawr: option for using CIA outside its wavenumber validity range
c       + ciatr: option for using CIA outside its temperature validity range
c       + sens_*: options for sensitivities computation
c
c     Outputs:
c       + kCIA: collision-induced absorption coefficient (m-¹)
c       + dkCIA_dP: sensitivity of kCIA to total pressure (m-¹/atm)
c       + dkCIA_dT: sensitivity of kCIA to temperature (m-¹/K)
c       + dkCIA_dxco2: sensitivity of kCIA to the concentration of CO2 (m-¹)
c

c     I/O
      double precision T,nu
      integer Nmol
      double precision P,Ps(1:Nmol_max)
      integer mol_index(1:Nmol_max)
      logical ucia
      integer ciac
      logical ciawr,ciatr
      logical sens_P,sens_T,sens_x
      double precision kCIA
      double precision dkCIA_dP
      double precision dkCIA_dT
      double precision dkCIA_dxco2
c     temp
      integer imol,molf,mol
      double precision k1,k2,k,amg
      double precision damg_dP,damg_dT,damg_dxco2
      double precision xco2,dT,T2,kT2,k1T2,k2T2,dkdT
      integer strlen
      character*(Nchar_mx) label
      label='subroutine CIA'


      kCIA=0.0D+0
      if (ucia) then
c     ----------------------------------------------------------------------------------------------
c     computation of the CO2 number density (amagats)
c     ----------------------------------------------------------------------------------------------

c     looking for partial pressure of CO2
         molf=0
         do mol=1,Nmol
            if (mol_index(mol).eq.2) then
               molf=1
               imol=mol
               goto 123
            endif
         enddo
c     imol should be the index of CO2: we are looking at Ps(imol)
 123     continue
         if (molf.eq.0) then
            write(*,*) 'Error from routine CIA:'
            write(*,*) 'CO2 could not be identified:'
            do mol=1,Nmol
               write(*,*) 'mol_index(',mol,')=',mol_index(mol),
     &              ' Ps(',mol,')=',Ps(mol)
            enddo
            stop
         endif ! molf=0
         amg=273.15D+0/T*Ps(imol) ! amagats of CO2
c     no division by standard pressure (1013.0D+0 Pa) is needed since Ps is given in atm
c     ----------------------------------------------------------------------------------------------
         if (((.not.ciawr).and.(nu.ge.0.0D+0).and.(nu.le.250.0D+0))
     &        .or.(ciawr)) then
            if (((.not.ciatr).and.(T.ge.200.0D+0).and.(T.le.800.0D+0))
     &           .or.(ciatr)) then
               if (ciac.eq.1) then ! use CIA program from Gruszka
                  call getspc(T,nu,k) ! k is in cm^-1.amagat^-2
c     Debug
c                  write(*,*) 'using getspc'
c     Debug
               else if (ciac.eq.2) then ! use CIA data from Baranov 2004
                  call baranov(T,nu,k) ! k is in cm^-1.amagat^-2
               else ! use both
                  call getspc(T,nu,k1)
                  call baranov(T,nu,k2)
                  k=k1+k2 ! k is in cm^-1.amagat^-2
               endif
               kCIA=k*1.0D+2*amg**2.0D+0 ! m^-1
c     -------------------------------------------------------------------
c     sensitivities 
                  if (sens_P) then
                     damg_dP=amg/P
                     dkCIA_dP=1.0D+2*k*2.0D+0*amg*damg_dP
                  endif
                  if (sens_T) then
                     damg_dT=-amg/T
                     dT=1.0D+0 ! K
                     T2=T+dT
                     if (ciac.eq.1) then ! use CIA program from Gruszka
                        call getspc(T2,nu,kT2) ! kT2 is in cm^-1.amagat^2
                     else if (ciac.eq.2) then ! use CIA data from Baranov 2004
                        call baranov(T2,nu,kT2) ! k is in cm^-1.amagat^2
                     else       ! use both
                        call getspc(T2,nu,k1T2)
                        call baranov(T2,nu,k2T2)
                        kT2=k1T2+k2T2 ! k is in cm^-1.amagat^2
                     endif
                     dkdT=(kT2-k)/dT
                     dkCIA_dT=1.0D+2*
     &                    (dkdT*amg**2.0D+0+k*2.0D+0*amg*damg_dT)
                  endif
                  if (sens_x) then
                     xco2=Ps(imol)/P
                     damg_dxco2=amg/xco2
                     dkCIA_dxco2=1.0D+2*k*2.0D+0*amg*damg_dxco2
                  endif
c     -------------------------------------------------------------------
            endif
         endif
      endif
      
      return
      end