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
      use mod_phys_lmdz_para, only : is_master

      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
         if (is_master) print*,'----------------------------------------------------'
         if (is_master) 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
           if (is_master) then 
             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/'
           endif
           call abort
         else
           do iT=1,nT
             read(34,fmat1) dummy_var1,nres,Ttemp
             if(nS.ne.nres)then
               if (is_master) 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'
               endif
               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
           if (is_master) then
             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/'
           endif
           call abort
         else
           do iT=1,nT
             read(35,fmat2) dummy_var2,nres,Ttemp
             if(nS.ne.nres)then
               if (is_master) 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'
               endif
               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)
         if (is_master) then
           print*,'interpolateH2O_MTCKDcont: At wavenumber ',wn,' cm^-1'
           print*,'   temperature ',temp,' K'
           print*,'   H2O pressure ',presS,' Pa'
           print*,'   air pressure ',presF,' Pa'
         endif
!$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