module rad_correlatedk_read_opacity_tables_mod

implicit none

contains

      subroutine rad_correlatedk_read_opacity_tables 

!==================================================================
!
!     Purpose
!     -------
!     Set up gaseous absorption parameters used by the radiation code.
!     This subroutine is a replacement for the old 'setrad', which contained
!     both absorption and scattering data.
!
!     Authors
!     -------
!     Adapted and generalised from the NASA Ames code by Robin Wordsworth (2009)
!     Added double gray case by Jeremy Leconte (2012)
!     New HITRAN continuum data section by RW (2012)
!     Modern traceur.def & corrk recombing scheme by J.Vatant d'Ollone (2020)
!
!     Summary
!     -------
!
!==================================================================

      use radinc_h, only: corrkdir, banddir, L_NSPECTI, L_NSPECTV, &
                          L_NGAUSS, L_NPREF, L_NTREF, L_REFVAR, L_PINT
      use radcommon_h, only : pgasref,pfgasref,pgasmin,pgasmax
      use radcommon_h, only : tgasref,tgasmin,tgasmax
      use radcommon_h, only : gasv,gasi,FZEROI,FZEROV,gweight
      use radcommon_h, only : wrefvar,WNOI,WNOV
      use datafile_mod, only: datadir
      use comcstfi_mod, only: mugaz
      use gases_h, only: gnom, ngasmx, &
                         igas_H2, igas_H2O, igas_He, igas_N2, igas_CH4, &
                         igas_CO2, igas_O2
      use ioipsl_getin_p_mod, only: getin_p
      use callkeys_mod, only: varactive,varfixed,graybody,callgasvis,&
		continuum
      use tracer_h, only : nqtot, moderntracdef, is_recomb, noms
      use rad_correlatedk_online_recombination_mod, only: rad_correlatedk_recombination_setup,        &
                corrk_recombin, use_premix, nrecomb_tot, rcb2tot_idx
      use rad_correlatedk_continuum_interpolation_mod, only: rad_correlatedk_continuum_interpolation
      use util_lagrange_interpolation_mod, only: lagrange4
      use mod_phys_lmdz_para, only : is_master, bcast
      implicit none

!==================================================================

      logical file_ok

      integer n, nt, np, nh, ng, nw, m, i

      character(len=200) :: file_id
      character(len=500) :: file_path

      ! ALLOCATABLE ARRAYS -- AS 12/2011
      REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi8, gasv8 	!read by master
      character*20,allocatable,DIMENSION(:) :: gastype ! for check with gnom, read by master

      real*8 x, xi(4), yi(4), ans, p
!     For gray case (JL12)
      real kappa_IR, kappa_VI, IR_VI_wnlimit
      integer nVI_limit,nIR_limit

      integer ngas, igas, jgas

      double precision testcont ! for continuum absorption initialisation

      if (.not. moderntracdef) use_premix=.true. ! Added by JVO for compatibility with 'old' traceur.def
      
      if (use_premix) then ! use_premix flag added by JVO, thus if pure recombining then premix is skipped

       if (is_master) then ! only the master needs read the file
!=======================================================================
!       Load variable species data, exit if we have wrong database
        file_id='/corrk_data/' // TRIM(corrkdir) // '/Q.dat'
        file_path=TRIM(datadir)//TRIM(file_id)

        ! check that the file exists
        inquire(FILE=file_path,EXIST=file_ok)
        if(.not.file_ok) then
         write(*,*)'The file ',TRIM(file_path)
         write(*,*)'was not found by rad_correlatedk_read_opacity_tables .F90, exiting.'
         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(*,*)'Also check that the corrkdir you chose in callphys.def exists.'
         call abort_physic("rad_correlatedk_read_opacity_tables ", "Unable to read file", 1)
        endif

        ! check that database matches varactive toggle
        open(111,file=TRIM(file_path),form='formatted')
        read(111,*) ngas

        if(moderntracdef) then
           ! JVO 20 - TODO : Sanity check with nspcrad !
           print*, 'Warning : Sanity check on # of gases still not implemented !!'
        else
         if(ngas.ne.ngasmx)then
           print*,'Number of gases in radiative transfer data (',ngas,') does not', &
                  'match that in gases.def (',ngasmx,'), exiting.'
           call abort_physic("rad_correlatedk_read_opacity_tables ", "Number of gases in radiative transfer data does not match that in gases.def", 1)
         endif 
        endif ! of if (moderntracdef)

        if(ngas.eq.1 .and. (varactive.or.varfixed))then
         print*,'You have varactive/fixed=.true. but the database [',  &
                     corrkdir(1:LEN_TRIM(corrkdir)),           &
                '] has no variable species, exiting.'
         call abort_physic("rad_correlatedk_read_opacity_tables ", "You have varactive/fixed=.true. but the database has no variable species",1)
        elseif(ngas.gt.5 .or. ngas.lt.1)then
         print*,ngas,' species in database [',               &
                     corrkdir(1:LEN_TRIM(corrkdir)),           &
                '], radiative code cannot handle this.'
         call abort_physic("rad_correlatedk_read_opacity_tables ", "No gas or too many gases for radiative code", 1)
        endif 

        ! allocate gastype and read from Q.dat
        ALLOCATE( gastype( ngas ) )

        do igas=1,ngas
         read(111,*) gastype(igas)
         if(corrk_recombin) then
            print*,'Premix gas ',igas,' is ',gastype(igas)
         else
           print*,'Gas ',igas,' is ',gastype(igas)
         endif
        enddo ! of do igas=1,ngas

        ! get array size, load the coefficients
        open(111,file=TRIM(file_path),form='formatted')
        read(111,*) L_REFVAR
        ALLOCATE( WREFVAR(L_REFVAR) )
        read(111,*) wrefvar
        close(111)
      endif ! of if (is_master)
      
      ! share the information with all cores
      call bcast(ngas)
      call bcast(L_REFVAR)
      if (.not.is_master) ALLOCATE(WREFVAR(L_REFVAR))
      call bcast(WREFVAR)
      

      if(L_REFVAR.gt.1 .and. (.not.varactive) .and. (.not.varfixed))then
         print*,'You have varactive and varfixed=.false. and the database [', &
                     corrkdir(1:LEN_TRIM(corrkdir)),           &
                '] has a variable species.'
         call abort_physic("rad_correlatedk_read_opacity_tables ", "You have varactive and varfixed=.false. and the database has a variable species",1)
      endif

      if (moderntracdef) then
          ! JVO 20 - TODO : Sanity check with nspcrad !
          print*, 'Warning : Sanity check on name of gases still not implemented !!'
      else
        ! Check that gastype and gnom match
        if (is_master) then
         do igas=1,ngas
           print*,'Gas ',igas,' is ',trim(gnom(igas))
           if (trim(gnom(igas)).ne.trim(gastype(igas))) then
              print*,'Name of a gas in radiative transfer data (',trim(gastype(igas)),') does not ', &
                   'match that in gases.def (',trim(gnom(igas)),'), exiting. You should compare ', &
                   'gases.def with Q.dat in your radiative transfer directory.' 
              call abort_physic("rad_correlatedk_read_opacity_tables ", "Name of a gas in radiative transfer data does not match that in gases.def",1)
           endif
         enddo
        endif ! of if (is_master)
        print*,'Confirmed gas match in radiative transfer and gases.def!'
      endif

      ! display the values
      print*,'Variable gas volume mixing ratios:'
      do n=1,L_REFVAR
         !print*,n,'.',wrefvar(n),' kg/kg' ! pay attention!
         print*,n,'.',wrefvar(n),' mol/mol'
      end do
      print*,''
      
      else
        L_REFVAR = 1
      endif ! of if (use_premix)

!=======================================================================
!     Set the weighting in g-space

      if (is_master) then ! only the master needs read the file
        file_id='/corrk_data/' // TRIM(corrkdir) // '/g.dat'
        file_path=TRIM(datadir)//TRIM(file_id)

        ! check that the file exists
        inquire(FILE=file_path,EXIST=file_ok)
        if(.not.file_ok) then
         write(*,*)'The file ',TRIM(file_path)
         write(*,*)'was not found by rad_correlatedk_read_opacity_tables .F90, exiting.'
         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(*,*)'Also check that the corrkdir you chose in callphys.def exists.'
         call abort_physic("rad_correlatedk_read_opacity_tables ", "Unable to read file", 1)
        endif
      
        ! check the array size is correct, load the coefficients
        open(111,file=TRIM(file_path),form='formatted')
        read(111,*) L_NGAUSS
        ALLOCATE( GWEIGHT(L_NGAUSS) )
        read(111,*) gweight
        close(111)
 
        ! display the values
        print*,'Correlated-k g-space grid:'
        do n=1,L_NGAUSS
         print*,n,'.',gweight(n)
        end do
        print*,''
      endif ! of if (is_master)

      ! share the information with all cores
      call bcast(L_NGAUSS)
      if (.not.is_master) ALLOCATE(GWEIGHT(L_NGAUSS))
      call bcast(GWEIGHT)

!=======================================================================
!     Set the reference pressure and temperature arrays.  These are
!     the pressures and temperatures at which we have k-coefficients.

!-----------------------------------------------------------------------
! pressure

      if (is_master) then ! only the master needs read the file
        file_id='/corrk_data/' // TRIM(corrkdir) // '/p.dat'
        file_path=TRIM(datadir)//TRIM(file_id)

        ! check that the file exists
        inquire(FILE=file_path,EXIST=file_ok)
        if(.not.file_ok) then
         write(*,*)'The file ',TRIM(file_path)
         write(*,*)'was not found by rad_correlatedk_read_opacity_tables .F90, exiting.'
         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(*,*)'Also check that the corrkdir you chose in callphys.def exists.'
         call abort_physic("rad_correlatedk_read_opacity_tables ", "Unable to read file", 1)
        endif
     
        ! get array size, load the coefficients
        open(111,file=TRIM(file_path),form='formatted')
        read(111,*) L_NPREF
        ALLOCATE( PGASREF(L_NPREF) )
        read(111,*) pgasref
        close(111)
        L_PINT = (L_NPREF-1)*5+1
        ALLOCATE( PFGASREF(L_PINT) )

        ! display the values
        print*,'Correlated-k pressure grid (mBar):'
        do n=1,L_NPREF
         print*,n,'. 1 x 10^',pgasref(n),' mBar'
        end do
        print*,''
      endif ! of if (is_master)

      ! share the information with all cores
      call bcast(L_NPREF)
      if (.not.is_master) ALLOCATE(PGASREF(L_NPREF))
      call bcast(PGASREF)
      call bcast(L_PINT)
      if (.not.is_master) ALLOCATE(PFGASREF(L_PINT)) ! values computed below

      ! save the min / max matrix values
      pgasmin = 10.0**pgasref(1)
      pgasmax = 10.0**pgasref(L_NPREF)

      ! interpolate to finer grid, adapted to uneven grids
      do n=1,L_NPREF-1
         do m=1,5
            pfgasref((n-1)*5+m) = pgasref(n)+(m-1)*(pgasref(n+1) - pgasref(n))/5.
         end do
      end do
      pfgasref(L_PINT) = pgasref(L_NPREF)

!-----------------------------------------------------------------------
! temperature

      if (is_master) then ! only the master needs read the file
        file_id='/corrk_data/' // TRIM(corrkdir) // '/T.dat'
        file_path=TRIM(datadir)//TRIM(file_id)

        ! check that the file exists
        inquire(FILE=file_path,EXIST=file_ok)
        if(.not.file_ok) then
         write(*,*)'The file ',TRIM(file_path)
         write(*,*)'was not found by rad_correlatedk_read_opacity_tables .F90, exiting.'
         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(*,*)'Also check that the corrkdir you chose in callphys.def exists.'
         call abort_physic("rad_correlatedk_read_opacity_tables ", "Unable to read file",1)
        endif

        ! get array size, load the coefficients
        open(111,file=TRIM(file_path),form='formatted')
        read(111,*) L_NTREF
        ALLOCATE( TGASREF(L_NTREF) )
        read(111,*) tgasref
        close(111)

        ! display the values
        print*,'Correlated-k temperature grid:'
        do n=1,L_NTREF
         print*,n,'.',tgasref(n),' K'
        end do
      endif ! of if (is_master)

      ! share the information with all cores
      call bcast(L_NTREF)
      if (.not.is_master) ALLOCATE(TGASREF(L_NTREF))
      call bcast(TGASREF)
      
      ! save the min / max matrix values
      tgasmin = tgasref(1)
      tgasmax = tgasref(L_NTREF)

      if(corrk_recombin) then ! even if no premix we keep a L_REFVAR=1 to store precombined at firstcall (see
         ALLOCATE(gasi8(L_NTREF,L_NPREF,L_REFVAR+nrecomb_tot,L_NSPECTI,L_NGAUSS))
         ALLOCATE(gasv8(L_NTREF,L_NPREF,L_REFVAR+nrecomb_tot,L_NSPECTV,L_NGAUSS))
      else
         ALLOCATE(gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS))
         ALLOCATE(gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS))
      endif

! JVO 20 : In these gasi/gasi8 matrix  we stack in same dim. variable specie and species to recombine (to keep code small)

!-----------------------------------------------------------------------
! allocate the multidimensional arrays in radcommon_h
     if(corrk_recombin) then
        ALLOCATE(gasi(L_NTREF,L_PINT,L_REFVAR+nrecomb_tot,L_NSPECTI,L_NGAUSS))
        ALLOCATE(gasv(L_NTREF,L_PINT,L_REFVAR+nrecomb_tot,L_NSPECTV,L_NGAUSS))
        ! display the values
        print*,''
        print*,'Correlated-k matrix size:'
        print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR+nrecomb_tot,' (',L_REFVAR,'+',nrecomb_tot,'),',L_NGAUSS,']' 
      else
        ALLOCATE(gasi(L_NTREF,L_PINT,L_REFVAR,L_NSPECTI,L_NGAUSS))
        ALLOCATE(gasv(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS))
        ! display the values
        print*,''
        print*,'Correlated-k matrix size:' 
        print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR,',',L_NGAUSS,']' 
      endif

!=======================================================================
!     Get gaseous k-coefficients and interpolate onto finer pressure grid


!        wavelength used to separate IR from VI in graybody. We will need that anyway
         IR_VI_wnlimit=3000.
         write(*,*)"graybody: Visible / Infrared separation set at",10000./IR_VI_wnlimit,"um"
	 
	 nVI_limit=0
	 Do nw=1,L_NSPECTV
	    if ((WNOV(nw).gt.IR_VI_wnlimit).and.(L_NSPECTV.gt.1)) then
	       nVI_limit=nw-1
	       exit
	    endif
	 End do
	 nIR_limit=L_NSPECTI
	 Do nw=1,L_NSPECTI
	    if ((WNOI(nw).gt.IR_VI_wnlimit).and.(L_NSPECTI.gt.1)) then
	       nIR_limit=nw-1
	       exit
	    endif
	 End do

      if (graybody) then
!        constant absorption coefficient in visible
         write(*,*)"graybody: constant absorption coefficient in visible:"
         kappa_VI=-100000.
         call getin_p("kappa_VI",kappa_VI)
         write(*,*)" kappa_VI = ",kappa_VI
	 kappa_VI=kappa_VI*1.e4* mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule         
      
!        constant absorption coefficient in IR
         write(*,*)"graybody: constant absorption coefficient in InfraRed:"
         kappa_IR=-100000.
         call getin_p("kappa_IR",kappa_IR)
         write(*,*)" kappa_IR = ",kappa_IR	 
         kappa_IR=kappa_IR*1.e4* mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule 

         write(*,*)"graybody: Visible / Infrared separation set at band: IR=",nIR_limit,", VI=",nVI_limit
	       
      Else
         kappa_VI=1.e-30      
         kappa_IR=1.e-30        
      End if

      if (is_master) then ! only the master needs to fill gasv8

       ! VISIBLE
       if (callgasvis) then
            
        ! Looking for premixed corrk files corrk_gcm_VI.dat if needed
        if (use_premix) then
           ! print*,corrkdir(1:4)
           if ((corrkdir(1:4).eq.'null'))then   !(TRIM(corrkdir).eq.'null_LowTeffStar')) then
              gasv8(:,:,:,:,:)=0.0
              print*,'using no corrk data'
              print*,'Visible corrk gaseous absorption is set to zero if graybody=F'
           else
              file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat' 
              file_path=TRIM(datadir)//TRIM(file_id)
              
              ! check that the file exists
              inquire(FILE=file_path,EXIST=file_ok)
              if(.not.file_ok) then
                 write(*,*)'The file ',TRIM(file_path)
                 write(*,*)'was not found by rad_correlatedk_read_opacity_tables .F90.'
                 write(*,*)'Are you sure you have absorption data for these bands?'
                 call abort_physic("rad_correlatedk_read_opacity_tables ", "No absorption data found", 1)
              endif
           
              open(111,file=TRIM(file_path),form='formatted')
              read(111,*) gasv8(:,:,:L_REFVAR,:,:)
              close(111)
	   end if
        else
           gasv8(:,:,1,:,:)=0.0 ! dummy but needs to be initialized
        endif ! use_premix
        
        ! Looking for others files corrk_gcm_VI_XXX.dat if needed
        if (corrk_recombin) then
          do igas=1,nrecomb_tot

            file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI_'//trim(adjustl(noms(rcb2tot_idx(igas))))//'.dat' 
            file_path=TRIM(datadir)//TRIM(file_id)
             
            ! check that the file exists
            inquire(FILE=file_path,EXIST=file_ok)
            if(.not.file_ok) then
               write(*,*)'The file ',TRIM(file_path)
               write(*,*)'was not found by rad_correlatedk_read_opacity_tables .F90.'
               write(*,*)'Are you sure you have absorption data for this specie at these bands?'
               call abort_physic("rad_correlatedk_read_opacity_tables ", "No absorption data found", 1)
            endif
         
            open(111,file=TRIM(file_path),form='formatted')
            read(111,*) gasv8(:,:,L_REFVAR+igas,:,:)
            close(111)
          enddo                           
        endif ! corrk_recombin

        if(nVI_limit.eq.0) then
         gasv8(:,:,:,:,:)= gasv8(:,:,:,:,:)+kappa_VI
        else if (nVI_limit.eq.L_NSPECTV) then
         gasv8(:,:,:,:,:)= gasv8(:,:,:,:,:)+kappa_IR
        else
         gasv8(:,:,:,1:nVI_limit,:)= gasv8(:,:,:,1:nVI_limit,:)+kappa_IR
         gasv8(:,:,:,nVI_limit+1:L_NSPECTV,:)= gasv8(:,:,:,nVI_limit+1:L_NSPECTV,:)+kappa_VI
        end if
           
       else 
            print*,'Visible corrk gaseous absorption is set to zero.'
            gasv8(:,:,:,:,:)=0.0
       endif ! callgasvis
     
     endif ! of (is_master)
     
     ! share the information with all cores
     call bcast(gasv8)
     
      if (is_master) then ! only the master needs to fill gasi8
      ! INFRA-RED
      
      ! Looking for premixed corrk files corrk_gcm_IR.dat if needed
      if (use_premix) then
        if ((corrkdir(1:4).eq.'null'))then       !.or.(TRIM(corrkdir).eq.'null_LowTeffStar')) then
           print*,'Infrared corrk gaseous absorption is set to zero if graybody=F'
           gasi8(:,:,:,:,:)=0.0
        else 
           file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat' 
           file_path=TRIM(datadir)//TRIM(file_id)
        
           ! check that the file exists
           inquire(FILE=file_path,EXIST=file_ok)
           if(.not.file_ok) then
              write(*,*)'The file ',TRIM(file_path)
              write(*,*)'was not found by rad_correlatedk_read_opacity_tables .F90.'
              write(*,*)'Are you sure you have absorption data for these bands?'
              call abort_physic("rad_correlatedk_read_opacity_tables ", "No absorption data found", 1)
           endif
         
           open(111,file=TRIM(file_path),form='formatted')
           read(111,*) gasi8(:,:,:L_REFVAR,:,:)
           close(111)

           ! 'fzero' is a currently unused feature that allows optimisation
           ! of the radiative transfer by neglecting bands where absorption
           ! is close to zero. As it could be useful in the future, this 
           ! section of the code has been kept commented and not erased.
           ! RW 7/3/12.

           do nw=1,L_NSPECTI
              fzeroi(nw) = 0.d0
!              do nt=1,L_NTREF
!                 do np=1,L_NPREF
!                    do nh=1,L_REFVAR
!                       do ng = 1,L_NGAUSS
!                          if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
!                             fzeroi(nw)=fzeroi(nw)+1.d0
!                          endif
!                       end do
!                    end do
!                 end do
!              end do
!              fzeroi(nw)=fzeroi(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
           end do

           do nw=1,L_NSPECTV
              fzerov(nw) = 0.d0
!              do nt=1,L_NTREF
!                 do np=1,L_NPREF
!                    do nh=1,L_REFVAR
!                       do ng = 1,L_NGAUSS
!                          if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
!                             fzerov(nw)=fzerov(nw)+1.d0
!                          endif
!                       end do
!                    end do
!                 end do
!              end do
!              fzerov(nw)=fzerov(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
           end do

        endif ! if corrkdir=null

      else
         gasi8(:,:,1,:,:)=0.0 ! dummy but needs to be initialized
      endif ! use_premix
    
      ! Looking for others files corrk_gcm_IR_XXX.dat if needed
      if (corrk_recombin) then
        do igas=1,nrecomb_tot

          file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR_'//trim(adjustl(noms(rcb2tot_idx(igas))))//'.dat' 
          file_path=TRIM(datadir)//TRIM(file_id)
           
          ! check that the file exists
          inquire(FILE=file_path,EXIST=file_ok)
          if(.not.file_ok) then
             write(*,*)'The file ',TRIM(file_path)
             write(*,*)'was not found by rad_correlatedk_read_opacity_tables .F90.'
             write(*,*)'Are you sure you have absorption data for this specie at these bands?'
             call abort_physic("rad_correlatedk_read_opacity_tables ", "No absorption data found", 1)
          endif

          open(111,file=TRIM(file_path),form='formatted')
          read(111,*) gasi8(:,:,L_REFVAR+igas,:,:)
          close(111)
        enddo                           
      endif ! corrk_recombin

      endif ! of if (is_master)

      ! share the information with all cores
      call bcast(gasi8)
      call bcast(fzeroi)
      call bcast(fzerov)
      
      if(nIR_limit.eq.0) then
         gasi8(:,:,:,:,:)= gasi8(:,:,:,:,:)+kappa_VI
      else if (nIR_limit.eq.L_NSPECTI) then
	 gasi8(:,:,:,:,:)= gasi8(:,:,:,:,:)+kappa_IR
      else
	 gasi8(:,:,:,1:nIR_limit,:)= gasi8(:,:,:,1:nIR_limit,:)+kappa_IR
	 gasi8(:,:,:,nIR_limit+1:L_NSPECTI,:)= gasi8(:,:,:,nIR_limit+1:L_NSPECTI,:)+kappa_VI
      end if


!     Take log10 of the values - this is what we will interpolate.
!     Smallest value is 1.0E-200.

      do nt=1,L_NTREF
         do np=1,L_NPREF
            do nh=1,L_REFVAR+nrecomb_tot
               do ng = 1,L_NGAUSS

                  do nw=1,L_NSPECTV
                     if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then
                        gasv8(nt,np,nh,nw,ng) = log10(gasv8(nt,np,nh,nw,ng))
                     else
                        gasv8(nt,np,nh,nw,ng) = -200.0
                     end if
                  end do

                  do nw=1,L_NSPECTI
                     if(gasi8(nt,np,nh,nw,ng).gt.1.0d-200) then
                        gasi8(nt,np,nh,nw,ng) = log10(gasi8(nt,np,nh,nw,ng))
                     else
                        gasi8(nt,np,nh,nw,ng) = -200.0
                     end if
                  end do
                  
               end do
            end do
         end do
      end do

!     Interpolate the values:  first the longwave

      do nt=1,L_NTREF
         do nh=1,L_REFVAR+nrecomb_tot
            do nw=1,L_NSPECTI
               do ng=1,L_NGAUSS

!     First, the initial interval

                  n = 1 
                  do m=1,5
                     x     = pfgasref(m)
                     xi(1) = pgasref(n)
                     xi(2) = pgasref(n+1)
                     xi(3) = pgasref(n+2)
                     xi(4) = pgasref(n+3)
                     yi(1) = gasi8(nt,n,nh,nw,ng)
                     yi(2) = gasi8(nt,n+1,nh,nw,ng)
                     yi(3) = gasi8(nt,n+2,nh,nw,ng)
                     yi(4) = gasi8(nt,n+3,nh,nw,ng)
                     call lagrange4(x,xi,yi,ans)
                     gasi(nt,m,nh,nw,ng) = 10.0**ans
                  end do 
                  
                  do n=2,L_NPREF-2
                     do m=1,5
                        i     = (n-1)*5+m
                        x     = pfgasref(i)
                        xi(1) = pgasref(n-1)
                        xi(2) = pgasref(n)
                        xi(3) = pgasref(n+1)
                        xi(4) = pgasref(n+2)
                        yi(1) = gasi8(nt,n-1,nh,nw,ng)
                        yi(2) = gasi8(nt,n,nh,nw,ng)
                        yi(3) = gasi8(nt,n+1,nh,nw,ng)
                        yi(4) = gasi8(nt,n+2,nh,nw,ng)
                        call lagrange4(x,xi,yi,ans)
                        gasi(nt,i,nh,nw,ng) = 10.0**ans
                     end do 
                  end do

!     Now, get the last interval

                  n = L_NPREF-1                 
                  do m=1,5
                     i     = (n-1)*5+m
                     x     = pfgasref(i)
                     xi(1) = pgasref(n-2)
                     xi(2) = pgasref(n-1)
                     xi(3) = pgasref(n)
                     xi(4) = pgasref(n+1)
                     yi(1) = gasi8(nt,n-2,nh,nw,ng)
                     yi(2) = gasi8(nt,n-1,nh,nw,ng)
                     yi(3) = gasi8(nt,n,nh,nw,ng)
                     yi(4) = gasi8(nt,n+1,nh,nw,ng)
                     call lagrange4(x,xi,yi,ans)
                     gasi(nt,i,nh,nw,ng) = 10.0**ans
                  end do  

!     Fill the last pressure point

                  gasi(nt,L_PINT,nh,nw,ng) = &
                       10.0**gasi8(nt,L_NPREF,nh,nw,ng)

               end do
            end do
         end do
      end do

!     Interpolate the values:  now the shortwave

      do nt=1,L_NTREF
         do nh=1,L_REFVAR+nrecomb_tot
            do nw=1,L_NSPECTV
               do ng=1,L_NGAUSS

!     First, the initial interval

                  n = 1 
                  do m=1,5
                     x     = pfgasref(m)
                     xi(1) = pgasref(n)
                     xi(2) = pgasref(n+1)
                     xi(3) = pgasref(n+2)
                     xi(4) = pgasref(n+3)
                     yi(1) = gasv8(nt,n,nh,nw,ng)
                     yi(2) = gasv8(nt,n+1,nh,nw,ng)
                     yi(3) = gasv8(nt,n+2,nh,nw,ng)
                     yi(4) = gasv8(nt,n+3,nh,nw,ng)
                     call lagrange4(x,xi,yi,ans)
                     gasv(nt,m,nh,nw,ng) = 10.0**ans
                  end do 
                  
                  do n=2,L_NPREF-2
                     do m=1,5
                        i     = (n-1)*5+m
                        x     = pfgasref(i)
                        xi(1) = pgasref(n-1)
                        xi(2) = pgasref(n)
                        xi(3) = pgasref(n+1)
                        xi(4) = pgasref(n+2)
                        yi(1) = gasv8(nt,n-1,nh,nw,ng)
                        yi(2) = gasv8(nt,n,nh,nw,ng)
                        yi(3) = gasv8(nt,n+1,nh,nw,ng)
                        yi(4) = gasv8(nt,n+2,nh,nw,ng)
                        call lagrange4(x,xi,yi,ans)
                        gasv(nt,i,nh,nw,ng) = 10.0**ans
                     end do 
                  end do

!     Now, get the last interval

                  n = L_NPREF-1
                  do m=1,5
                     i     = (n-1)*5+m
                     x     = pfgasref(i)
                     xi(1) = pgasref(n-2)
                     xi(2) = pgasref(n-1)
                     xi(3) = pgasref(n)
                     xi(4) = pgasref(n+1)
                     yi(1) = gasv8(nt,n-2,nh,nw,ng)
                     yi(2) = gasv8(nt,n-1,nh,nw,ng)
                     yi(3) = gasv8(nt,n,nh,nw,ng)
                     yi(4) = gasv8(nt,n+1,nh,nw,ng)
                     call lagrange4(x,xi,yi,ans)
                     gasv(nt,i,nh,nw,ng) = 10.0**ans
                  end do  

!     Fill the last pressure point

                  gasv(nt,L_PINT,nh,nw,ng) = &
                      10.0**gasv8(nt,L_NPREF,nh,nw,ng)
                  
               end do
            end do
         end do
      end do

! Allocate and initialise arrays for corrk_recombin
      if (corrk_recombin) call rad_correlatedk_recombination_setup

!=======================================================================
!     Initialise the continuum absorption data
      if(continuum)then
      
       do igas=1,ngasmx ! we loop on all pairs of molecules that have data available
       ! data can be downloaded from https://web.lmd.jussieu.fr/~lmdz/planets/generic/datagcm/continuum_data/
        if (igas .eq. igas_N2) then
         file_id='/continuum_data/' // 'N2-N2_continuum_70-500K_2025.dat'
         file_path=TRIM(datadir)//TRIM(file_id)
         call rad_correlatedk_continuum_interpolation(file_path,igas_N2,igas_N2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
         do jgas=1,ngasmx
          if (jgas .eq. igas_H2) then
           file_id='/continuum_data/' // 'H2-N2_continuum_40-600K_2025.dat'
           file_path=TRIM(datadir)//TRIM(file_id)
           call rad_correlatedk_continuum_interpolation(file_path,igas_N2,igas_H2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
          elseif (jgas .eq. igas_O2) then
           file_id='/continuum_data/' // 'O2-N2_continuum_100-500K_2025.dat'
           file_path=TRIM(datadir)//TRIM(file_id)
           call rad_correlatedk_continuum_interpolation(file_path,igas_N2,igas_O2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
          elseif (jgas .eq. igas_CH4) then
           file_id='/continuum_data/' // 'CH4-N2_continuum_40-600K_2025.dat'
           file_path=TRIM(datadir)//TRIM(file_id)
           call rad_correlatedk_continuum_interpolation(file_path,igas_N2,igas_CH4,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
          endif
         enddo
        elseif (igas .eq. igas_O2) then
         file_id='/continuum_data/' // 'O2-O2_continuum_100-400K_2025.dat'
         file_path=TRIM(datadir)//TRIM(file_id)
	 call rad_correlatedk_continuum_interpolation(file_path,igas_O2,igas_O2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
	 do jgas=1,ngasmx
          if (jgas .eq. igas_CO2) then
           file_id='/continuum_data/' // 'CO2-O2_continuum_150-600K_2025.dat'
           file_path=TRIM(datadir)//TRIM(file_id)
	   call rad_correlatedk_continuum_interpolation(file_path,igas_CO2,igas_O2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
	  endif
         enddo
        elseif (igas .eq. igas_H2) then
         file_id='/continuum_data/' // 'H2-H2_continuum_40-7000K_2025.dat'
         file_path=TRIM(datadir)//TRIM(file_id)
         call rad_correlatedk_continuum_interpolation(file_path,igas_H2,igas_H2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
         do jgas=1,ngasmx
          if (jgas .eq. igas_CH4) then
           file_id='/continuum_data/' // 'H2-CH4_continuum_40-600K_2025.dat'
           file_path=TRIM(datadir)//TRIM(file_id)
           call rad_correlatedk_continuum_interpolation(file_path,igas_H2,igas_CH4,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
          elseif (jgas .eq. igas_He) then
           file_id='/continuum_data/' // 'H2-He_continuum_40-5500K_2025.dat'
           file_path=TRIM(datadir)//TRIM(file_id)
	   call rad_correlatedk_continuum_interpolation(file_path,igas_H2,igas_He,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
          endif
         enddo	 
        elseif (igas .eq. igas_CH4) then
         file_id='/continuum_data/' // 'CH4-CH4_continuum_40-500K_2025.dat'
         file_path=TRIM(datadir)//TRIM(file_id)
         call rad_correlatedk_continuum_interpolation(file_path,igas_CH4,igas_CH4,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
        elseif (igas .eq. igas_H2O) then
         file_id='/continuum_data/' // 'H2O-H2O_continuum_100-2000K_2025.dat'
         file_path=TRIM(datadir)//TRIM(file_id)
         call rad_correlatedk_continuum_interpolation(file_path,igas_H2O,igas_H2O,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
         do jgas=1,ngasmx
          if (jgas .eq. igas_N2) then
           file_id='/continuum_data/' // 'H2O-N2_continuum_100-2000K_2025.dat'
           file_path=TRIM(datadir)//TRIM(file_id)
           call rad_correlatedk_continuum_interpolation(file_path,igas_H2O,igas_N2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
          elseif (jgas .eq. igas_O2) then
           file_id='/continuum_data/' // 'H2O-O2_continuum_100-2000K_2025.dat'
           file_path=TRIM(datadir)//TRIM(file_id)
           call rad_correlatedk_continuum_interpolation(file_path,igas_H2O,igas_O2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
          elseif (jgas .eq. igas_CO2) then
           file_id='/continuum_data/' // 'H2O-CO2_continuum_100-800K_2025.dat'
           file_path=TRIM(datadir)//TRIM(file_id)
           call rad_correlatedk_continuum_interpolation(file_path,igas_H2O,igas_CO2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
          endif
         enddo	 
        elseif (igas .eq. igas_CO2) then
         file_id='/continuum_data/' // 'CO2-CO2_continuum_100-800K_2025.dat'
         file_path=TRIM(datadir)//TRIM(file_id)
         call rad_correlatedk_continuum_interpolation(file_path,igas_CO2,igas_CO2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
	 do jgas=1,ngasmx
          if (jgas .eq. igas_H2) then
           file_id='/continuum_data/' // 'CO2-H2_continuum_100-800K_2025.dat'
           file_path=TRIM(datadir)//TRIM(file_id)
	   call rad_correlatedk_continuum_interpolation(file_path,igas_CO2,igas_H2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
	  elseif (jgas .eq. igas_CH4) then
           file_id='/continuum_data/' // 'CO2-CH4_continuum_100-800K_2025.dat'
           file_path=TRIM(datadir)//TRIM(file_id)
           call rad_correlatedk_continuum_interpolation(file_path,igas_CO2,igas_CH4,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
          endif
         enddo
        endif
       enddo ! igas=1,ngasmx
       
      endif ! continuum flag

      print*,'----------------------------------------------------'
      print*,'And that`s all we have. It`s possible that other'
      print*,'continuum absorption may be present, but if it is we'
      print*,'don`t yet have data for it...'
      print*,''

!     Deallocate local arrays
      DEALLOCATE( gasi8 )
      DEALLOCATE( gasv8 )
      if (is_master) DEALLOCATE( gastype )

    end subroutine rad_correlatedk_read_opacity_tables 

end module rad_correlatedk_read_opacity_tables_mod