subroutine interpolateH2O_self_foreign(wn,temp,presS,presF,abcoef,firstcall,ind)

!==================================================================
!     
!     Purpose
!     -------
!     Calculates the H2O continuum self and foreign opacity
!     using lookup tables from MT_CKD version 3.3
!
!     Authors
!     -------
!     M. Turbet (2020)
!     Adapted from the file of R. Wordsworth
!     
!==================================================================

      use datafile_mod, only: datadir
      implicit none

      ! input
      double precision wn                 ! wavenumber             (cm^-1)
      double precision temp               ! temperature            (Kelvin)
      double precision presS              ! self-pressure          (Pascals)
      double precision presF              ! foreign (air) pressure (Pascals)

      ! output
      double precision abcoef             ! absorption coefficient (m^-1)

      integer nS,nT, iT
      parameter(nS=2001)    ! number of wavenumbers
      parameter(nT=39)      ! number of temperatures

      double precision kB
      parameter(kB=1.3806488e-23)

      double precision amagatS, amagatF, abcoefS, abcoefF, Nmolec
      double precision wn_arr(nS)
      double precision temp_arr(nT)
      double precision abs_arrS(nS,nT)
      double precision abs_arrF(nS,nT)
      double precision data_tmp(nT)
      
      character*43 dummy_var1
      character*46 dummy_var2

      integer nres

      double precision Ttemp

      character(LEN=*), parameter :: fmat1 = "(A43,I6,F6.1)"
      character(LEN=*), parameter :: fmat2 = "(A46,I6,F6.1)"


      integer k,ind
      logical firstcall

      save wn_arr, temp_arr, abs_arrS, abs_arrF !read by master

      character*100 dt_file
      integer strlen,ios

      amagatS=(273.15/temp)*(presS/101325.0)
      amagatF=(273.15/temp)*(presF/101325.0)

      if(firstcall)then ! called by sugas_corrk only
         print*,'----------------------------------------------------'
         print*,'Initialising H2O continuum from MT_CKD data...'

!     1.1 Open the ASCII files

!$OMP MASTER

         ! self broadening
         dt_file=TRIM(datadir)//'/continuum_data/H2O-H2O_continuum_MT_CKD3.3.cia'
         open(34,file=dt_file,form='formatted',status='old',iostat=ios)
         if (ios.ne.0) then        ! file not found
           write(*,*) 'Error from interpolateH2O_cont SELF'
           write(*,*) 'Data file ',trim(dt_file),' not found.'
           write(*,*) 'Check that your path to datagcm:',trim(datadir)
           write(*,*) ' is correct. You can change it in callphys.def with:'
           write(*,*) ' datadir = /absolute/path/to/datagcm'
           write(*,*) ' Check that there is a H2O-H2O_continuum_MT_CKD3.3.cia'
           write(*,*)'Continuum file available here:'
           write(*,*)' https://www.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/continuum_data/'
           call abort
         else
           do iT=1,nT
             read(34,fmat1) dummy_var1,nres,Ttemp
             if(nS.ne.nres)then
               print*,'Resolution given in file: ',trim(dt_file)
               print*,'is ',nres,', which does not match nS.'
               print*,'Adjust nS value in interpolateH2O_MTCKD...F90'
               stop
             endif
             temp_arr(iT)=Ttemp
	     !write(*,*) 'H2O_H2O'
             do k=1,nS
               read(34,*) wn_arr(k),abs_arrS(k,it)
	       !write(*,*) nres,temp_arr(iT),wn_arr(k),abs_arrS(k,it)
             end do
	   end do
	 endif
         close(34)
	   
         ! foreign broadening
         dt_file=TRIM(datadir)//'/continuum_data/H2O-AIR_continuum_MT_CKD3.3.cia'
         open(35,file=dt_file,form='formatted',status='old',iostat=ios)
         if (ios.ne.0) then        ! file not found
           write(*,*) 'Error from interpolateH2O_cont FOREIGN'
           write(*,*) 'Data file ',trim(dt_file),' not found.'
           write(*,*)'Check that your path to datagcm:',trim(datadir)
           write(*,*)' is correct. You can change it in callphys.def with:'
           write(*,*)' datadir = /absolute/path/to/datagcm'
           write(*,*)' Check that there is a H2O-AIR_continuum_MT_CKD3.3.cia'
           write(*,*)'Continuum file available here:'
           write(*,*)' https://www.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/continuum_data/'
           call abort
         else
           do iT=1,nT
             read(35,fmat2) dummy_var2,nres,Ttemp
             if(nS.ne.nres)then
               print*,'Resolution given in file: ',trim(dt_file)
               print*,'is ',nres,', which does not match nS.'
               print*,'Adjust nS value in interpolateH2O_MTCKD...F90'
               stop
             endif
             temp_arr(iT)=Ttemp
	     !write(*,*) 'H2O_AIR'
             do k=1,nS
               read(35,*) wn_arr(k),abs_arrF(k,it)
	       !write(*,*) nres,temp_arr(iT),wn_arr(k),abs_arrF(k,it)
             end do
	   end do
	 endif
         close(35)

         print*,'interpolateH2O_MTCKDcont: At wavenumber ',wn,' cm^-1'
         print*,'   temperature ',temp,' K'
         print*,'   H2O pressure ',presS,' Pa'
         print*,'   air pressure ',presF,' Pa'
!$OMP END MASTER
!$OMP BARRIER
	 
      endif

      call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS,ind)
      !print*,'MTCKD - self absorption is ',abcoefS,' cm^2 molecule^-1'

      call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF,ind)
      !print*,'MTCKD - foreign absorption is ',abcoefF,' cm^2 molecule^-1'

! abcoefS and abcoefF are in cm-1 amagat-2
! First we multiply by amagat^2
      abcoef = abcoefS*amagatS + abcoefF*amagatF
      abcoef = abcoef*amagatS
! Then we convert cm-1 in m-1
      abcoef = 100.*abcoef
      
      !print*,'MTCKD_v3.3 : So the total absorption is ',abcoef,' m^-1'
      !print*,'for PH2O/PN2/TEMP/wn=',presS,presF,temp,wn

      return
    end subroutine interpolateH2O_self_foreign