Index: trunk/LMDZ.GENERIC/changelog.txt
===================================================================
--- trunk/LMDZ.GENERIC/changelog.txt	(revision 3621)
+++ trunk/LMDZ.GENERIC/changelog.txt	(revision 3641)
@@ -1994,2 +1994,6 @@
 Follow-up of recent changes in iostart (making 1D restarts) from r3562,
 Correspondingly adapt newstart.
+
+== 23/02/2025 == MT
+Implement the new generic continuum database
+More details here: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Continuum_Database
Index: trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90	(revision 3621)
+++ trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90	(revision 3641)
@@ -45,4 +45,6 @@
             elseif(igas.eq.igas_N2)then
                mugaz_c = mugaz_c + 28.01*gfrac(igas)
+	    elseif(igas.eq.igas_O2)then
+               mugaz_c = mugaz_c + 31.999*gfrac(igas)
             elseif(igas.eq.igas_H2)then
                mugaz_c = mugaz_c + 2.01*gfrac(igas)
@@ -98,4 +100,6 @@
             elseif(igas.eq.igas_N2)then
                cpp_c   = cpp_c   + 1.040*gfrac(igas)*28.01/mugaz_c
+            elseif(igas.eq.igas_O2)then
+               cpp_c   = cpp_c   + 0.918*gfrac(igas)*31.999/mugaz_c
             elseif(igas.eq.igas_H2)then
                cpp_c   = cpp_c   + 14.31*gfrac(igas)*2.01/mugaz_c
Index: trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90	(revision 3621)
+++ trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90	(revision 3641)
@@ -29,5 +29,6 @@
       use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, taurayvar, scalep
       use gases_h, only: ngasmx, vgas, gnom, gfrac, igas_CO2, igas_H2, &
-                           igas_H2O, igas_He, igas_N2, igas_CH4, igas_CO
+                           igas_H2O, igas_He, igas_N2, igas_CH4, igas_CO, &
+			   igas_O2
       use comcstfi_mod, only: g, mugaz, pi
 
Index: trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 3621)
+++ trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 3641)
@@ -13,6 +13,6 @@
       logical,save :: season,diurnal,tlocked,rings_shadow,lwrite
 !$OMP THREADPRIVATE(season,diurnal,tlocked,rings_shadow,lwrite)
-      logical,save :: callgasvis,continuum,graybody
-!$OMP THREADPRIVATE(callgasvis,continuum,graybody)
+      logical,save :: callgasvis,continuum,generic_continuum_database,graybody
+!$OMP THREADPRIVATE(callgasvis,continuum,generic_continuum_database,graybody)
       logical,save :: strictboundcorrk                                     
 !$OMP THREADPRIVATE(strictboundcorrk)
Index: trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 3621)
+++ trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 3641)
@@ -341,4 +341,12 @@
      call getin_p("continuum",continuum)
      if (is_master) write(*,*) trim(rname)//": continuum = ",continuum
+     
+     if (is_master) write(*,*) trim(rname)//&
+       ": use generic continuum opacity database ?"//&
+       " (matters only if callrad=T)"
+     generic_continuum_database=.true. ! default value
+     call getin_p("generic_continuum_database",generic_continuum_database)
+     if (is_master) write(*,*) trim(rname)//": generic_continuum_database = ", &
+     generic_continuum_database
  
      if (is_master) write(*,*) trim(rname)//": version for H2H2 CIA file ?"
Index: trunk/LMDZ.GENERIC/libf/phystd/interpolate_continuum.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/interpolate_continuum.F90	(revision 3641)
+++ trunk/LMDZ.GENERIC/libf/phystd/interpolate_continuum.F90	(revision 3641)
@@ -0,0 +1,809 @@
+
+     subroutine interpolate_continuum(filename,igas_X,igas_Y,c_WN,ind_WN,temp,pres_X,pres_Y,abs_coef,firstcall)
+
+!==================================================================
+!     
+!     Purpose
+!     -------
+!     Generic routine to calculate continuum opacities, using lookup tables provided here: https://web.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/continuum_data/
+!     More information on the data here: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Continuum_Database
+!
+!     Author
+!     -------
+!     M. Turbet (2025)
+!
+!==================================================================
+
+      use datafile_mod, only: datadir
+      use mod_phys_lmdz_para, only : is_master
+
+      use gases_h, only: ngasmx, gnom, &
+                         igas_H2, igas_H2O, igas_He, igas_N2, &
+                         igas_CH4, igas_CO2, igas_O2
+
+      use radinc_h, only: L_NSPECTI, L_NSPECTV
+
+      use radcommon_h, only : BWNV,BWNI,WNOI,WNOV
+
+
+      implicit none
+
+      ! input
+      integer ind_WN                      ! wavenumber index 
+      integer igas_X                      ! index of molecule X
+      integer igas_Y                      ! index of molecule Y
+      double precision temp               ! temperature (Kelvin)
+      double precision pres_X             ! partial pressure of molecule X (Pascals)
+      double precision pres_Y             ! partial pressure of molecule Y (Pascals)
+      character*200 filename              ! name of the lookup table
+      character*2 c_WN                    ! wavelength chanel: infrared (IR) or visible (VI)
+      logical firstcall
+
+      ! output
+      double precision abs_coef           ! absorption coefficient (m^-1)
+
+      ! intermediate variables
+      double precision amagat_X           ! density of molecule X (in amagat units)
+      double precision amagat_Y           ! density of molecule Y (in amagat units)
+
+      character(len=512) :: line
+
+      integer i, pos, iT, iW, iB, count_norm, igas
+
+      double precision temp_value, temp_abs, temp_wn
+
+      double precision z_temp
+
+      integer num_wn, num_T
+
+      double precision, dimension(:), allocatable :: temp_arr
+      double precision, dimension(:),   allocatable :: wn_arr
+      double precision, dimension(:,:), allocatable :: abs_arr
+
+      integer ios
+
+      ! Temperature array, continuum absorption grid for the pair N2-N2
+      integer num_T_N2N2
+      double precision, dimension(:), allocatable :: temp_arr_N2N2
+      double precision, dimension(:,:), allocatable :: abs_arr_N2N2_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_N2N2_VI
+      save num_T_N2N2,temp_arr_N2N2,abs_arr_N2N2_IR,abs_arr_N2N2_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair O2-O2
+      integer num_T_O2O2
+      double precision, dimension(:), allocatable :: temp_arr_O2O2
+      double precision, dimension(:,:), allocatable :: abs_arr_O2O2_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_O2O2_VI
+      save num_T_O2O2,temp_arr_O2O2,abs_arr_O2O2_IR,abs_arr_O2O2_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair H2-H2
+      integer num_T_H2H2
+      double precision, dimension(:), allocatable :: temp_arr_H2H2
+      double precision, dimension(:,:), allocatable :: abs_arr_H2H2_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_H2H2_VI
+      save num_T_H2H2,temp_arr_H2H2,abs_arr_H2H2_IR,abs_arr_H2H2_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair CO2-CO2
+      integer num_T_CO2CO2
+      double precision, dimension(:), allocatable :: temp_arr_CO2CO2
+      double precision, dimension(:,:), allocatable :: abs_arr_CO2CO2_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_CO2CO2_VI
+      save num_T_CO2CO2,temp_arr_CO2CO2,abs_arr_CO2CO2_IR,abs_arr_CO2CO2_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair CH4-CH4
+      integer num_T_CH4CH4
+      double precision, dimension(:), allocatable :: temp_arr_CH4CH4
+      double precision, dimension(:,:), allocatable :: abs_arr_CH4CH4_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_CH4CH4_VI
+      save num_T_CH4CH4,temp_arr_CH4CH4,abs_arr_CH4CH4_IR,abs_arr_CH4CH4_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair H2O-H2O
+      integer num_T_H2OH2O
+      double precision, dimension(:), allocatable :: temp_arr_H2OH2O
+      double precision, dimension(:,:), allocatable :: abs_arr_H2OH2O_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_H2OH2O_VI
+      save num_T_H2OH2O,temp_arr_H2OH2O,abs_arr_H2OH2O_IR,abs_arr_H2OH2O_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair H2-He
+      integer num_T_H2He
+      double precision, dimension(:), allocatable :: temp_arr_H2He
+      double precision, dimension(:,:), allocatable :: abs_arr_H2He_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_H2He_VI
+      save num_T_H2He,temp_arr_H2He,abs_arr_H2He_IR,abs_arr_H2He_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair H2-CH4
+      integer num_T_H2CH4
+      double precision, dimension(:), allocatable :: temp_arr_H2CH4
+      double precision, dimension(:,:), allocatable :: abs_arr_H2CH4_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_H2CH4_VI
+      save num_T_H2CH4,temp_arr_H2CH4,abs_arr_H2CH4_IR,abs_arr_H2CH4_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair CO2-H2
+      integer num_T_CO2H2
+      double precision, dimension(:), allocatable :: temp_arr_CO2H2
+      double precision, dimension(:,:), allocatable :: abs_arr_CO2H2_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_CO2H2_VI
+      save num_T_CO2H2,temp_arr_CO2H2,abs_arr_CO2H2_IR,abs_arr_CO2H2_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair CO2-CH4
+      integer num_T_CO2CH4
+      double precision, dimension(:), allocatable :: temp_arr_CO2CH4
+      double precision, dimension(:,:), allocatable :: abs_arr_CO2CH4_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_CO2CH4_VI
+      save num_T_CO2CH4,temp_arr_CO2CH4,abs_arr_CO2CH4_IR,abs_arr_CO2CH4_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair N2-H2
+      integer num_T_N2H2
+      double precision, dimension(:), allocatable :: temp_arr_N2H2
+      double precision, dimension(:,:), allocatable :: abs_arr_N2H2_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_N2H2_VI
+      save num_T_N2H2,temp_arr_N2H2,abs_arr_N2H2_IR,abs_arr_N2H2_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair N2-CH4
+      integer num_T_N2CH4
+      double precision, dimension(:), allocatable :: temp_arr_N2CH4
+      double precision, dimension(:,:), allocatable :: abs_arr_N2CH4_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_N2CH4_VI
+      save num_T_N2CH4,temp_arr_N2CH4,abs_arr_N2CH4_IR,abs_arr_N2CH4_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair CO2-O2
+      integer num_T_CO2O2
+      double precision, dimension(:), allocatable :: temp_arr_CO2O2
+      double precision, dimension(:,:), allocatable :: abs_arr_CO2O2_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_CO2O2_VI
+      save num_T_CO2O2,temp_arr_CO2O2,abs_arr_CO2O2_IR,abs_arr_CO2O2_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair N2-O2
+      integer num_T_N2O2
+      double precision, dimension(:), allocatable :: temp_arr_N2O2
+      double precision, dimension(:,:), allocatable :: abs_arr_N2O2_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_N2O2_VI
+      save num_T_N2O2,temp_arr_N2O2,abs_arr_N2O2_IR,abs_arr_N2O2_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair H2O-N2
+      integer num_T_H2ON2
+      double precision, dimension(:), allocatable :: temp_arr_H2ON2
+      double precision, dimension(:,:), allocatable :: abs_arr_H2ON2_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_H2ON2_VI
+      save num_T_H2ON2,temp_arr_H2ON2,abs_arr_H2ON2_IR,abs_arr_H2ON2_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair H2O-O2
+      integer num_T_H2OO2
+      double precision, dimension(:), allocatable :: temp_arr_H2OO2
+      double precision, dimension(:,:), allocatable :: abs_arr_H2OO2_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_H2OO2_VI
+      save num_T_H2OO2,temp_arr_H2OO2,abs_arr_H2OO2_IR,abs_arr_H2OO2_VI !read by master
+
+      ! Temperature array, continuum absorption grid for the pair H2O-CO2
+      integer num_T_H2OCO2
+      double precision, dimension(:), allocatable :: temp_arr_H2OCO2
+      double precision, dimension(:,:), allocatable :: abs_arr_H2OCO2_IR
+      double precision, dimension(:,:), allocatable :: abs_arr_H2OCO2_VI
+      save num_T_H2OCO2,temp_arr_H2OCO2,abs_arr_H2OCO2_IR,abs_arr_H2OCO2_VI !read by master
+
+
+
+      if(firstcall)then ! called by sugas_corrk only
+        if (is_master) print*,'----------------------------------------------------'
+        if (is_master) print*,'Initialising continuum (interpolate_continuum routine) from ', trim(filename)
+
+!$OMP MASTER
+
+        open(unit=33, file=trim(filename), status="old", action="read",iostat=ios)
+
+        if (ios.ne.0) then        ! file not found
+          if (is_master) then
+            write(*,*) 'Error from interpolate_continuum routine'
+            write(*,*) 'Data file ',trim(filename),' 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(*,*) 'Also check that the continuum data is there.'
+            write(*,*) 'Latest continuum data can be downloaded here:'
+            write(*,*) 'https://web.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/continuum_data/'
+          endif
+          call abort
+        endif
+
+        ! We read the first line of the file to get the number of temperatures provided in the data file
+        read(33, '(A)') line
+
+        i = 1
+        iT = 0
+
+        do while (i .lt. len_trim(line))
+          pos = index(line(i:), 'T=')
+	  if (pos == 0) exit
+          i = i + pos
+          iT = iT + 1
+          read(line(i+2:i+10), '(E9.2)') temp_value
+        end do
+
+        num_T=iT ! num_T is the number of temperatures provided in the data file
+	
+	! We read all the remaining lines of the file to get the number of wavenumbers provided in the data file
+        iW = 0
+        do
+          read(33,*, end=501) line
+          iW = iW + 1
+        end do
+	
+501 continue
+	
+        num_wn=iW ! num_wn is the number of wavenumbers provided in the data file
+	
+        close(33)
+
+        if(.not.allocated(temp_arr)) allocate(temp_arr(num_T)) 
+        if(.not.allocated(wn_arr))  allocate(wn_arr(num_wn)) 
+        if(.not.allocated(abs_arr)) allocate(abs_arr(num_wn,num_T)) 
+	
+	! We now open and read the file a second time to extract the temperature array, wavenumber array and continuum absorption data
+
+        open(unit=33, file=trim(filename), status="old", action="read")
+	
+	! We extract the temperature array (temp_arr)
+	
+        read(33, '(A)') line
+
+        i = 1
+        iT = 0
+
+        do while (i .lt. len_trim(line))
+          pos = index(line(i:), 'T=')
+	  if (pos == 0) exit
+          i = i + pos
+          iT = iT + 1
+          read(line(i+2:i+10), '(E9.2)') temp_arr(iT)
+        end do
+	
+	! We extract the wavenumber array (wn_arr) and continuum absorption (abs_arr)
+
+        do iW=1,num_wn
+          read(33,*) wn_arr(iW), (abs_arr(iW, iT), iT=1,num_T)
+        end do
+
+        close(33)
+	
+        print*,'We read continuum absorption data for the pair ', trim(gnom(igas_X)),'-',trim(gnom(igas_Y))
+	print*,'Temperature grid of the dataset: ', temp_arr(:)
+	
+	! We loop on all molecular pairs with available continuum data and fill the corresponding array 
+	
+        if ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CO2)) then
+          num_T_CO2CO2=num_T
+          if(.not.allocated(temp_arr_CO2CO2)) allocate(temp_arr_CO2CO2(num_T_CO2CO2)) 
+          if(.not.allocated(abs_arr_CO2CO2_VI)) allocate(abs_arr_CO2CO2_VI(L_NSPECTV,num_T_CO2CO2)) 
+          if(.not.allocated(abs_arr_CO2CO2_IR)) allocate(abs_arr_CO2CO2_IR(L_NSPECTI,num_T_CO2CO2)) 
+          temp_arr_CO2CO2(:)=temp_arr(:)
+          abs_arr_CO2CO2_VI(:,:)=0.
+          abs_arr_CO2CO2_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2CO2_VI,abs_arr_CO2CO2_IR,num_T_CO2CO2)
+	elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_N2)) then
+          num_T_N2N2=num_T
+          if(.not.allocated(temp_arr_N2N2)) allocate(temp_arr_N2N2(num_T_N2N2)) 
+          if(.not.allocated(abs_arr_N2N2_VI)) allocate(abs_arr_N2N2_VI(L_NSPECTV,num_T_N2N2)) 
+          if(.not.allocated(abs_arr_N2N2_IR)) allocate(abs_arr_N2N2_IR(L_NSPECTI,num_T_N2N2)) 
+          temp_arr_N2N2(:)=temp_arr(:)
+          abs_arr_N2N2_VI(:,:)=0.
+          abs_arr_N2N2_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_N2N2_VI,abs_arr_N2N2_IR,num_T_N2N2)
+	elseif ((igas_X .eq. igas_O2) .and. (igas_Y .eq. igas_O2)) then
+          num_T_O2O2=num_T
+          if(.not.allocated(temp_arr_O2O2)) allocate(temp_arr_O2O2(num_T_O2O2)) 
+          if(.not.allocated(abs_arr_O2O2_VI)) allocate(abs_arr_O2O2_VI(L_NSPECTV,num_T_O2O2)) 
+          if(.not.allocated(abs_arr_O2O2_IR)) allocate(abs_arr_O2O2_IR(L_NSPECTI,num_T_O2O2)) 
+          temp_arr_O2O2(:)=temp_arr(:)
+          abs_arr_O2O2_VI(:,:)=0.
+          abs_arr_O2O2_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_O2O2_VI,abs_arr_O2O2_IR,num_T_O2O2)
+	elseif ((igas_X .eq. igas_CH4) .and. (igas_Y .eq. igas_CH4)) then
+          num_T_CH4CH4=num_T
+          if(.not.allocated(temp_arr_CH4CH4)) allocate(temp_arr_CH4CH4(num_T_CH4CH4)) 
+          if(.not.allocated(abs_arr_CH4CH4_VI)) allocate(abs_arr_CH4CH4_VI(L_NSPECTV,num_T_CH4CH4)) 
+          if(.not.allocated(abs_arr_CH4CH4_IR)) allocate(abs_arr_CH4CH4_IR(L_NSPECTI,num_T_CH4CH4)) 
+          temp_arr_CH4CH4(:)=temp_arr(:)
+          abs_arr_CH4CH4_VI(:,:)=0.
+          abs_arr_CH4CH4_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CH4CH4_VI,abs_arr_CH4CH4_IR,num_T_CH4CH4)
+	elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_H2)) then
+          num_T_H2H2=num_T
+          if(.not.allocated(temp_arr_H2H2)) allocate(temp_arr_H2H2(num_T_H2H2)) 
+          if(.not.allocated(abs_arr_H2H2_VI)) allocate(abs_arr_H2H2_VI(L_NSPECTV,num_T_H2H2)) 
+          if(.not.allocated(abs_arr_H2H2_IR)) allocate(abs_arr_H2H2_IR(L_NSPECTI,num_T_H2H2)) 
+          temp_arr_H2H2(:)=temp_arr(:)
+          abs_arr_H2H2_VI(:,:)=0.
+          abs_arr_H2H2_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2H2_VI,abs_arr_H2H2_IR,num_T_H2H2)
+	elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_H2O)) then
+          num_T_H2OH2O=num_T
+          if(.not.allocated(temp_arr_H2OH2O)) allocate(temp_arr_H2OH2O(num_T_H2OH2O)) 
+          if(.not.allocated(abs_arr_H2OH2O_VI)) allocate(abs_arr_H2OH2O_VI(L_NSPECTV,num_T_H2OH2O)) 
+          if(.not.allocated(abs_arr_H2OH2O_IR)) allocate(abs_arr_H2OH2O_IR(L_NSPECTI,num_T_H2OH2O)) 
+          temp_arr_H2OH2O(:)=temp_arr(:)
+          abs_arr_H2OH2O_VI(:,:)=0.
+          abs_arr_H2OH2O_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2OH2O_VI,abs_arr_H2OH2O_IR,num_T_H2OH2O)
+	elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_H2)) then
+          num_T_N2H2=num_T
+          if(.not.allocated(temp_arr_N2H2)) allocate(temp_arr_N2H2(num_T_N2H2)) 
+          if(.not.allocated(abs_arr_N2H2_VI)) allocate(abs_arr_N2H2_VI(L_NSPECTV,num_T_N2H2)) 
+          if(.not.allocated(abs_arr_N2H2_IR)) allocate(abs_arr_N2H2_IR(L_NSPECTI,num_T_N2H2)) 
+          temp_arr_N2H2(:)=temp_arr(:)
+          abs_arr_N2H2_VI(:,:)=0.
+          abs_arr_N2H2_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_N2H2_VI,abs_arr_N2H2_IR,num_T_N2H2)
+	elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_O2)) then
+          num_T_N2O2=num_T
+          if(.not.allocated(temp_arr_N2O2)) allocate(temp_arr_N2O2(num_T_N2O2)) 
+          if(.not.allocated(abs_arr_N2O2_VI)) allocate(abs_arr_N2O2_VI(L_NSPECTV,num_T_N2O2)) 
+          if(.not.allocated(abs_arr_N2O2_IR)) allocate(abs_arr_N2O2_IR(L_NSPECTI,num_T_N2O2)) 
+          temp_arr_N2O2(:)=temp_arr(:)
+          abs_arr_N2O2_VI(:,:)=0.
+          abs_arr_N2O2_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_N2O2_VI,abs_arr_N2O2_IR,num_T_N2O2)
+	elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_CH4)) then
+          num_T_N2CH4=num_T
+          if(.not.allocated(temp_arr_N2CH4)) allocate(temp_arr_N2CH4(num_T_N2CH4)) 
+          if(.not.allocated(abs_arr_N2CH4_VI)) allocate(abs_arr_N2CH4_VI(L_NSPECTV,num_T_N2CH4)) 
+          if(.not.allocated(abs_arr_N2CH4_IR)) allocate(abs_arr_N2CH4_IR(L_NSPECTI,num_T_N2CH4)) 
+          temp_arr_N2CH4(:)=temp_arr(:)
+          abs_arr_N2CH4_VI(:,:)=0.
+          abs_arr_N2CH4_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_N2CH4_VI,abs_arr_N2CH4_IR,num_T_N2CH4)
+	elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_O2)) then
+          num_T_CO2O2=num_T
+          if(.not.allocated(temp_arr_CO2O2)) allocate(temp_arr_CO2O2(num_T_CO2O2)) 
+          if(.not.allocated(abs_arr_CO2O2_VI)) allocate(abs_arr_CO2O2_VI(L_NSPECTV,num_T_CO2O2)) 
+          if(.not.allocated(abs_arr_CO2O2_IR)) allocate(abs_arr_CO2O2_IR(L_NSPECTI,num_T_CO2O2)) 
+          temp_arr_CO2O2(:)=temp_arr(:)
+          abs_arr_CO2O2_VI(:,:)=0.
+          abs_arr_CO2O2_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2O2_VI,abs_arr_CO2O2_IR,num_T_CO2O2)
+	elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_CH4)) then
+          num_T_H2CH4=num_T
+          if(.not.allocated(temp_arr_H2CH4)) allocate(temp_arr_H2CH4(num_T_H2CH4)) 
+          if(.not.allocated(abs_arr_H2CH4_VI)) allocate(abs_arr_H2CH4_VI(L_NSPECTV,num_T_H2CH4)) 
+          if(.not.allocated(abs_arr_H2CH4_IR)) allocate(abs_arr_H2CH4_IR(L_NSPECTI,num_T_H2CH4)) 
+          temp_arr_H2CH4(:)=temp_arr(:)
+          abs_arr_H2CH4_VI(:,:)=0.
+          abs_arr_H2CH4_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2CH4_VI,abs_arr_H2CH4_IR,num_T_H2CH4)
+	elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_He)) then
+          num_T_H2He=num_T
+          if(.not.allocated(temp_arr_H2He)) allocate(temp_arr_H2He(num_T_H2He)) 
+          if(.not.allocated(abs_arr_H2He_VI)) allocate(abs_arr_H2He_VI(L_NSPECTV,num_T_H2He)) 
+          if(.not.allocated(abs_arr_H2He_IR)) allocate(abs_arr_H2He_IR(L_NSPECTI,num_T_H2He)) 
+          temp_arr_H2He(:)=temp_arr(:)
+          abs_arr_H2He_VI(:,:)=0.
+          abs_arr_H2He_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2He_VI,abs_arr_H2He_IR,num_T_H2He)
+	elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_N2)) then
+          num_T_H2ON2=num_T
+          if(.not.allocated(temp_arr_H2ON2)) allocate(temp_arr_H2ON2(num_T_H2ON2)) 
+          if(.not.allocated(abs_arr_H2ON2_VI)) allocate(abs_arr_H2ON2_VI(L_NSPECTV,num_T_H2ON2)) 
+          if(.not.allocated(abs_arr_H2ON2_IR)) allocate(abs_arr_H2ON2_IR(L_NSPECTI,num_T_H2ON2)) 
+          temp_arr_H2ON2(:)=temp_arr(:)
+          abs_arr_H2ON2_VI(:,:)=0.
+          abs_arr_H2ON2_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2ON2_VI,abs_arr_H2ON2_IR,num_T_H2ON2)
+	elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_O2)) then
+          num_T_H2OO2=num_T
+          if(.not.allocated(temp_arr_H2OO2)) allocate(temp_arr_H2OO2(num_T_H2OO2)) 
+          if(.not.allocated(abs_arr_H2OO2_VI)) allocate(abs_arr_H2OO2_VI(L_NSPECTV,num_T_H2OO2)) 
+          if(.not.allocated(abs_arr_H2OO2_IR)) allocate(abs_arr_H2OO2_IR(L_NSPECTI,num_T_H2OO2)) 
+          temp_arr_H2OO2(:)=temp_arr(:)
+          abs_arr_H2OO2_VI(:,:)=0.
+          abs_arr_H2OO2_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2OO2_VI,abs_arr_H2OO2_IR,num_T_H2OO2)
+	elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_CO2)) then
+          num_T_H2OCO2=num_T
+          if(.not.allocated(temp_arr_H2OCO2)) allocate(temp_arr_H2OCO2(num_T_H2OCO2)) 
+          if(.not.allocated(abs_arr_H2OCO2_VI)) allocate(abs_arr_H2OCO2_VI(L_NSPECTV,num_T_H2OCO2)) 
+          if(.not.allocated(abs_arr_H2OCO2_IR)) allocate(abs_arr_H2OCO2_IR(L_NSPECTI,num_T_H2OCO2)) 
+          temp_arr_H2OCO2(:)=temp_arr(:)
+          abs_arr_H2OCO2_VI(:,:)=0.
+          abs_arr_H2OCO2_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2OCO2_VI,abs_arr_H2OCO2_IR,num_T_H2OCO2)
+	elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CO2)) then
+          num_T_CO2CO2=num_T
+          if(.not.allocated(temp_arr_CO2CO2)) allocate(temp_arr_CO2CO2(num_T_CO2CO2)) 
+          if(.not.allocated(abs_arr_CO2CO2_VI)) allocate(abs_arr_CO2CO2_VI(L_NSPECTV,num_T_CO2CO2)) 
+          if(.not.allocated(abs_arr_CO2CO2_IR)) allocate(abs_arr_CO2CO2_IR(L_NSPECTI,num_T_CO2CO2)) 
+          temp_arr_CO2CO2(:)=temp_arr(:)
+          abs_arr_CO2CO2_VI(:,:)=0.
+          abs_arr_CO2CO2_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2CO2_VI,abs_arr_CO2CO2_IR,num_T_CO2CO2)
+	elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_H2)) then
+          num_T_CO2H2=num_T
+          if(.not.allocated(temp_arr_CO2H2)) allocate(temp_arr_CO2H2(num_T_CO2H2)) 
+          if(.not.allocated(abs_arr_CO2H2_VI)) allocate(abs_arr_CO2H2_VI(L_NSPECTV,num_T_CO2H2)) 
+          if(.not.allocated(abs_arr_CO2H2_IR)) allocate(abs_arr_CO2H2_IR(L_NSPECTI,num_T_CO2H2)) 
+          temp_arr_CO2H2(:)=temp_arr(:)
+          abs_arr_CO2H2_VI(:,:)=0.
+          abs_arr_CO2H2_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2H2_VI,abs_arr_CO2H2_IR,num_T_CO2H2)
+	elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CH4)) then
+          num_T_CO2CH4=num_T
+          if(.not.allocated(temp_arr_CO2CH4)) allocate(temp_arr_CO2CH4(num_T_CO2CH4)) 
+          if(.not.allocated(abs_arr_CO2CH4_VI)) allocate(abs_arr_CO2CH4_VI(L_NSPECTV,num_T_CO2CH4)) 
+          if(.not.allocated(abs_arr_CO2CH4_IR)) allocate(abs_arr_CO2CH4_IR(L_NSPECTI,num_T_CO2CH4)) 
+          temp_arr_CO2CH4(:)=temp_arr(:)
+          abs_arr_CO2CH4_VI(:,:)=0.
+          abs_arr_CO2CH4_IR(:,:)=0.
+	  call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2CH4_VI,abs_arr_CO2CH4_IR,num_T_CO2CH4)  
+        endif ! igas_X / igas_Y condition
+	
+
+!$OMP END MASTER
+!$OMP BARRIER
+
+
+      endif ! firstcall
+
+      ! We loop on all molecular pairs with available continuum data and interpolate in the temperature field
+      ! Two options: we call visible (VI) or infrared (IR) tables, depending on the value of c_WN 
+      
+      if ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CO2)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_CO2CO2,num_T_CO2CO2)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_CO2CO2,num_T_CO2CO2,abs_coef,abs_arr_CO2CO2_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_CO2CO2,num_T_CO2CO2,abs_coef,abs_arr_CO2CO2_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif
+      elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_N2)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_N2N2,num_T_N2N2)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_N2N2,num_T_N2N2,abs_coef,abs_arr_N2N2_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_N2N2,num_T_N2N2,abs_coef,abs_arr_N2N2_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif
+      elseif ((igas_X .eq. igas_O2) .and. (igas_Y .eq. igas_O2)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_O2O2,num_T_O2O2)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_O2O2,num_T_O2O2,abs_coef,abs_arr_O2O2_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_O2O2,num_T_O2O2,abs_coef,abs_arr_O2O2_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif
+      elseif ((igas_X .eq. igas_CH4) .and. (igas_Y .eq. igas_CH4)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_CH4CH4,num_T_CH4CH4)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_CH4CH4,num_T_CH4CH4,abs_coef,abs_arr_CH4CH4_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_CH4CH4,num_T_CH4CH4,abs_coef,abs_arr_CH4CH4_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif	
+      elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_H2)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_H2H2,num_T_H2H2)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_H2H2,num_T_H2H2,abs_coef,abs_arr_H2H2_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_H2H2,num_T_H2H2,abs_coef,abs_arr_H2H2_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif
+      elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_H2O)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_H2OH2O,num_T_H2OH2O)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_H2OH2O,num_T_H2OH2O,abs_coef,abs_arr_H2OH2O_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_H2OH2O,num_T_H2OH2O,abs_coef,abs_arr_H2OH2O_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif
+      elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_H2)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_N2H2,num_T_N2H2)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_N2H2,num_T_N2H2,abs_coef,abs_arr_N2H2_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_N2H2,num_T_N2H2,abs_coef,abs_arr_N2H2_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif
+      elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_O2)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_N2O2,num_T_N2O2)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_N2O2,num_T_N2O2,abs_coef,abs_arr_N2O2_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_N2O2,num_T_N2O2,abs_coef,abs_arr_N2O2_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif
+      elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_CH4)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_N2CH4,num_T_N2CH4)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_N2CH4,num_T_N2CH4,abs_coef,abs_arr_N2CH4_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_N2CH4,num_T_N2CH4,abs_coef,abs_arr_N2CH4_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif
+      elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_O2)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_CO2O2,num_T_CO2O2)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_CO2O2,num_T_CO2O2,abs_coef,abs_arr_CO2O2_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_CO2O2,num_T_CO2O2,abs_coef,abs_arr_CO2O2_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif
+      elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_CH4)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_H2CH4,num_T_H2CH4)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_H2CH4,num_T_H2CH4,abs_coef,abs_arr_H2CH4_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_H2CH4,num_T_H2CH4,abs_coef,abs_arr_H2CH4_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif
+      elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_He)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_H2He,num_T_H2He)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_H2He,num_T_H2He,abs_coef,abs_arr_H2He_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_H2He,num_T_H2He,abs_coef,abs_arr_H2He_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif	
+      elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_N2)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_H2ON2,num_T_H2ON2)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_H2ON2,num_T_H2ON2,abs_coef,abs_arr_H2ON2_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_H2ON2,num_T_H2ON2,abs_coef,abs_arr_H2ON2_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif	
+      elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_O2)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_H2OO2,num_T_H2OO2)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_H2OO2,num_T_H2OO2,abs_coef,abs_arr_H2OO2_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_H2OO2,num_T_H2OO2,abs_coef,abs_arr_H2OO2_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif	
+      elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_CO2)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_H2OCO2,num_T_H2OCO2)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_H2OCO2,num_T_H2OCO2,abs_coef,abs_arr_H2OCO2_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_H2OCO2,num_T_H2OCO2,abs_coef,abs_arr_H2OCO2_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif
+      elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_H2)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_CO2H2,num_T_CO2H2)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_CO2H2,num_T_CO2H2,abs_coef,abs_arr_CO2H2_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_CO2H2,num_T_CO2H2,abs_coef,abs_arr_CO2H2_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif	
+      elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CH4)) then
+        call T_boundaries_continuum(z_temp,temp,temp_arr_CO2CH4,num_T_CO2CH4)
+	if(c_WN .eq. 'IR') then
+          call interpolate_T_abs_coeff(z_temp,temp_arr_CO2CH4,num_T_CO2CH4,abs_coef,abs_arr_CO2CH4_IR(ind_WN,:))
+	elseif(c_WN .eq. 'VI') then
+	  call interpolate_T_abs_coeff(z_temp,temp_arr_CO2CH4,num_T_CO2CH4,abs_coef,abs_arr_CO2CH4_VI(ind_WN,:))
+	else
+	  print*,'You must select visible (VI) or infrared (IR) canal.'
+          stop
+	endif									
+      endif ! igas_X / igas_Y condition
+      
+      ! We compute the values of amagat for molecules X and Y
+      amagat_X = (273.15/temp)*(pres_X/101325.0)
+      amagat_Y = (273.15/temp)*(pres_Y/101325.0)
+
+      ! We convert the absorption coefficient from cm^-1 amagat^-2 into m^-1
+      abs_coef=abs_coef*100.0*amagat_X*amagat_Y
+
+      !print*,'We have ',amagat_X,' amagats of molecule ', trim(gnom(igas_X))
+      !print*,'We have ',amagat_X,' amagats of molecule ', trim(gnom(igas_Y))
+      !print*,'So the absorption is ',abs_coef,' m^-1'
+      
+      return
+    end subroutine interpolate_continuum
+    
+    
+    subroutine interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr_in,abs_arr_out_VI,abs_arr_out_IR,num_T)
+    
+!==================================================================
+!     
+!     Purpose
+!     -------
+!     Interpolate the continuum data into the visible (VI) and infrared (IR) spectral chanels.
+!
+!     Author
+!     -------
+!     M. Turbet (2025)
+!
+!==================================================================
+
+      use radcommon_h, only : BWNV,BWNI,WNOI,WNOV
+      use radinc_h, only: L_NSPECTI, L_NSPECTV
+      use mod_phys_lmdz_para, only : is_master
+
+      implicit none
+            
+      integer iW, iB, count_norm
+      integer num_T
+      integer num_wn
+      double precision wn_arr(num_wn)
+      double precision abs_arr_in(num_wn,num_T)
+      double precision abs_arr_out_IR(L_NSPECTI,num_T)
+      double precision abs_arr_out_VI(L_NSPECTV,num_T)
+
+      ! First visible (VI) chanel
+
+      ! We get read of all the wavenumbers lower than the minimum wavenumber in the visible wavenumber grid
+      iW=1
+      do while((wn_arr(iW) .lt. BWNV(1)) .and. (iW .lt. num_wn))
+        iW=iW+1
+      enddo
+      
+      ! We compute the mean of the continuum absorption inside each wavenumber visible (VI) chanel      
+      do iB = 1, L_NSPECTV
+        count_norm=0
+        do while((wn_arr(iW) .lt. BWNV(iB+1)) .and. (iW .lt. num_wn))
+          abs_arr_out_VI(iB,:)=abs_arr_out_VI(iB,:)+abs_arr_in(iW,:)
+          count_norm=count_norm+1
+          iW=iW+1
+        enddo
+        if(count_norm .ge. 1) abs_arr_out_VI(iB,:)=abs_arr_out_VI(iB,:)/count_norm
+      end do
+      
+      ! Then infrared (IR) chanel
+      
+      ! We get read of all the wavenumbers lower than the minimum wavenumber in the infrared wavenumber grid
+      iW=1
+      do while((wn_arr(iW) .lt. BWNI(1)) .and. (iW .lt. num_wn))
+        iW=iW+1
+      enddo
+
+      ! We compute the mean of the continuum absorption inside each wavenumber visible (VI) chanel      
+      do iB = 1, L_NSPECTI
+        count_norm=0
+        do while((wn_arr(iW) .lt. BWNI(iB+1)) .and. (iW .lt. num_wn))
+          abs_arr_out_IR(iB,:)=abs_arr_out_IR(iB,:)+abs_arr_in(iW,:)
+          count_norm=count_norm+1
+          iW=iW+1
+        enddo
+        if(count_norm .ge. 1) abs_arr_out_IR(iB,:)=abs_arr_out_IR(iB,:)/count_norm
+      end do
+
+      if (is_master) then
+        print*, 'Continuum absorption, first temperature, visible (VI):'
+        do iB = 1, L_NSPECTV
+          print*,WNOV(iB),' cm-1',abs_arr_out_VI(iB,1), ' cm-1 amagat-2'
+        end do
+
+        print*, 'Continuum absorption, first temperature, infrared (IR):'
+        do iB = 1, L_NSPECTI
+          print*,WNOI(iB),' cm-1',abs_arr_out_IR(iB,1), ' cm-1 amagat-2'
+        end do
+      endif
+	
+      return
+    end subroutine interpolate_wn_abs_coeff
+
+
+    subroutine T_boundaries_continuum(z_temp,temp,temp_arr,num_T)
+    
+!==================================================================
+!     
+!     Purpose
+!     -------
+!     Check if the temperature is outside the boundaries of the continuum data temperatures.
+!
+!     Author
+!     -------
+!     M. Turbet (2025)
+!
+!==================================================================
+    
+      use callkeys_mod, only: strictboundcia
+      use mod_phys_lmdz_para, only : is_master
+
+      implicit none
+      
+      double precision z_temp
+      double precision temp
+      integer num_T
+      double precision temp_arr(num_T)
+      
+      z_temp=temp
+      
+      if(z_temp .lt. minval(temp_arr)) then
+        if (strictboundcia) then
+          if (is_master) then
+            print*,'Your temperatures are too low for this continuum dataset'
+            print*, 'Minimum temperature is ', minval(temp_arr), ' K'
+          endif
+          stop
+        else
+          z_temp=minval(temp_arr)
+        endif
+      elseif(z_temp .gt. maxval(temp_arr)) then
+        if (strictboundcia) then
+          if (is_master) then
+            print*,'Your temperatures are too high for this continuum dataset'
+            print*, 'Maximum temperature is ', maxval(temp_arr), ' K'
+          endif
+          stop
+        else
+          z_temp=maxval(temp_arr)
+        endif
+      endif
+      
+      return
+    end subroutine T_boundaries_continuum
+
+
+    subroutine interpolate_T_abs_coeff(z_temp,temp_arr,num_T,abs_coef,abs_arr)
+
+!==================================================================
+!     
+!     Purpose
+!     -------
+!     Interpolate in the continuum data using the temperature field
+!
+!     Author
+!     -------
+!     M. Turbet (2025)
+!
+!==================================================================
+
+      implicit none
+      
+      integer iT
+      double precision z_temp
+      integer num_T
+      double precision temp_arr(num_T)
+      
+      double precision abs_coef
+      double precision abs_arr(num_T)
+      
+      ! Check where to interpolate
+      iT=1
+      do while ( z_temp .gt. temp_arr(iT) )
+        iT=iT+1
+      end do
+      
+      ! We proceed to a simple linear interpolation using the two most nearby temperatures
+      if(iT .lt. num_T) then
+        abs_coef=abs_arr(iT-1)+(abs_arr(iT)-abs_arr(iT-1))*(z_temp-temp_arr(iT-1))/(temp_arr(iT)-temp_arr(iT-1))
+      else
+        abs_coef=abs_arr(iT)
+      endif
+      
+      !print*,'the absorption is ',abs_coef,' cm^-1 amagat^-2'
+
+      
+      return
+    end subroutine interpolate_T_abs_coeff
Index: trunk/LMDZ.GENERIC/libf/phystd/optci.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 3621)
+++ trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 3641)
@@ -13,7 +13,8 @@
   use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig
   use gases_h, only: gfrac, ngasmx, igas_N2, igas_He, igas_H2O, igas_H2, &
-                     igas_CH4, igas_CO2
+                     igas_CH4, igas_CO2, igas_O2
   use comcstfi_mod, only: g, r, mugaz
-  use callkeys_mod, only: kastprof,continuum,graybody,varspec
+  use callkeys_mod, only: kastprof,continuum,graybody,varspec, &
+                          generic_continuum_database
   use recombin_corrk_mod, only: corrk_recombin, gasi_recomb
   use tpindex_mod, only: tpindex
@@ -188,7 +189,77 @@
         if(continuum.and.(.not.graybody))then
            ! include continua if necessary
-           wn_cont = dble(wnoi(nw))
-           T_cont  = dble(TMID(k))
-           do igas=1,ngasmx
+	   
+	   if(generic_continuum_database)then
+	     T_cont  = dble(TMID(k))
+	     do igas=1,ngasmx
+	     
+               if(gfrac(igas).eq.-1)then ! variable
+                 p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
+               elseif(varspec) then
+                 p_cont  = dble(PMID(k)*scalep*FRACVAR(igas,k)*(1.-QVAR(k)))
+               else
+                 p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
+               endif
+	     
+	       dtemp=0.0
+	       
+	       if (igas .eq. igas_N2) then
+	        call interpolate_continuum('',igas_N2,igas_N2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                do jgas=1,ngasmx
+                 if (jgas .eq. igas_H2) then
+		  call interpolate_continuum('',igas_N2,igas_H2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                 elseif (jgas .eq. igas_O2) then
+		  call interpolate_continuum('',igas_N2,igas_O2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                 elseif (jgas .eq. igas_CH4) then
+		  call interpolate_continuum('',igas_N2,igas_CH4,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                 endif
+                enddo
+               elseif (igas .eq. igas_O2) then
+	        call interpolate_continuum('',igas_O2,igas_O2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+	        do jgas=1,ngasmx
+                 if (jgas .eq. igas_CO2) then
+		  call interpolate_continuum('',igas_CO2,igas_O2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+	         endif
+                enddo
+               elseif (igas .eq. igas_H2) then
+	        call interpolate_continuum('',igas_H2,igas_H2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                do jgas=1,ngasmx
+                 if (jgas .eq. igas_CH4) then
+		  call interpolate_continuum('',igas_H2,igas_CH4,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                 elseif (jgas .eq. igas_He) then
+		  call interpolate_continuum('',igas_H2,igas_He,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                 endif
+                enddo	 
+               elseif (igas .eq. igas_CH4) then
+	        call interpolate_continuum('',igas_CH4,igas_CH4,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+               elseif (igas .eq. igas_H2O) then
+	        call interpolate_continuum('',igas_H2O,igas_H2O,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                do jgas=1,ngasmx
+                 if (jgas .eq. igas_N2) then
+		  call interpolate_continuum('',igas_H2O,igas_N2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                 elseif (jgas .eq. igas_O2) then
+		  call interpolate_continuum('',igas_H2O,igas_O2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                 elseif (jgas .eq. igas_CO2) then
+		  call interpolate_continuum('',igas_H2O,igas_CO2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                 endif
+                enddo	 
+               elseif (igas .eq. igas_CO2) then
+	        call interpolate_continuum('',igas_CO2,igas_CO2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+	        do jgas=1,ngasmx
+                 if (jgas .eq. igas_H2) then
+		  call interpolate_continuum('',igas_CO2,igas_H2,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+	         elseif (jgas .eq. igas_CH4) then
+		  call interpolate_continuum('',igas_CO2,igas_CH4,'IR',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                 endif
+                enddo
+               endif
+	       
+	       DCONT = DCONT + dtemp
+	       
+	     enddo ! igas=1,ngasmx
+	   else ! generic_continuum_database
+            wn_cont = dble(wnoi(nw))
+            T_cont  = dble(TMID(k))
+            do igas=1,ngasmx
 
               if(gfrac(igas).eq.-1)then ! variable
@@ -295,5 +366,7 @@
            DCONT = DCONT*dz(k)
 
-        endif
+        endif ! generic_continuum_database
+	
+	endif ! continuum
 
         do ng=1,L_NGAUSS-1
Index: trunk/LMDZ.GENERIC/libf/phystd/optcv.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 3621)
+++ trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 3641)
@@ -12,7 +12,8 @@
   use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig
   use gases_h, only: gfrac, ngasmx, igas_H2, igas_H2O, igas_He, igas_N2, &
-                     igas_CH4, igas_CO2
+                     igas_CH4, igas_CO2, igas_O2
   use comcstfi_mod, only: g, r, mugaz
-  use callkeys_mod, only: kastprof,continuum,graybody,callgasvis,varspec
+  use callkeys_mod, only: kastprof,continuum,graybody,callgasvis,varspec, &
+                          generic_continuum_database
   use recombin_corrk_mod, only: corrk_recombin, gasv_recomb
   use tpindex_mod, only: tpindex
@@ -193,4 +194,77 @@
         if(continuum.and.(.not.graybody).and.callgasvis)then
            ! include continua if necessary
+	   
+	  if(generic_continuum_database)then
+	    T_cont  = dble(TMID(k))
+	    do igas=1,ngasmx
+	     
+              if(gfrac(igas).eq.-1)then ! variable
+                p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
+              elseif(varspec) then
+                p_cont  = dble(PMID(k)*scalep*FRACVAR(igas,k)*(1.-QVAR(k)))
+              else
+                p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
+              endif
+	     
+	      dtemp=0.0
+	      
+	      if (igas .eq. igas_N2) then
+	       call interpolate_continuum('',igas_N2,igas_N2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+               do jgas=1,ngasmx
+                if (jgas .eq. igas_H2) then
+		 call interpolate_continuum('',igas_N2,igas_H2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                elseif (jgas .eq. igas_O2) then
+		 call interpolate_continuum('',igas_N2,igas_O2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                elseif (jgas .eq. igas_CH4) then
+		 call interpolate_continuum('',igas_N2,igas_CH4,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                endif
+               enddo
+              elseif (igas .eq. igas_O2) then
+	       call interpolate_continuum('',igas_O2,igas_O2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+	       do jgas=1,ngasmx
+                if (jgas .eq. igas_CO2) then
+		 call interpolate_continuum('',igas_CO2,igas_O2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+	        endif
+               enddo
+              elseif (igas .eq. igas_H2) then
+	       call interpolate_continuum('',igas_H2,igas_H2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+               do jgas=1,ngasmx
+                if (jgas .eq. igas_CH4) then
+		 call interpolate_continuum('',igas_H2,igas_CH4,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                elseif (jgas .eq. igas_He) then
+		 call interpolate_continuum('',igas_H2,igas_He,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                endif
+               enddo	 
+              elseif (igas .eq. igas_CH4) then
+	       call interpolate_continuum('',igas_CH4,igas_CH4,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+              elseif (igas .eq. igas_H2O) then
+	       call interpolate_continuum('',igas_H2O,igas_H2O,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+               do jgas=1,ngasmx
+                if (jgas .eq. igas_N2) then
+		 call interpolate_continuum('',igas_H2O,igas_N2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                elseif (jgas .eq. igas_O2) then
+		 call interpolate_continuum('',igas_H2O,igas_O2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                elseif (jgas .eq. igas_CO2) then
+		 call interpolate_continuum('',igas_H2O,igas_CO2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                endif
+               enddo	 
+              elseif (igas .eq. igas_CO2) then
+	       call interpolate_continuum('',igas_CO2,igas_CO2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+	       do jgas=1,ngasmx
+                if (jgas .eq. igas_H2) then
+		 call interpolate_continuum('',igas_CO2,igas_H2,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+	        elseif (jgas .eq. igas_CH4) then
+		 call interpolate_continuum('',igas_CO2,igas_CH4,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                endif
+               enddo
+              endif
+	       
+	      DCONT = DCONT + dtemp
+	       
+	    enddo ! igas=1,ngasmx
+	    
+	  else ! generic_continuum_database
+	   
+	   
            wn_cont = dble(wnov(nw))
            T_cont  = dble(TMID(k))
@@ -293,6 +367,7 @@
            DCONT = DCONT*dz(k)
 
-        endif
-
+        endif ! generic_continuum_database
+        endif ! continuum
+	
         do ng=1,L_NGAUSS-1
 
Index: trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 3621)
+++ trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 3641)
@@ -31,8 +31,8 @@
       use gases_h, only: gnom, ngasmx, &
                          igas_H2, igas_H2O, igas_He, igas_N2, igas_CH4, &
-                         igas_CO2
+                         igas_CO2, igas_O2
       use ioipsl_getin_p_mod, only: getin_p
       use callkeys_mod, only: varactive,varfixed,graybody,callgasvis,&
-		continuum
+		continuum,generic_continuum_database
       use tracer_h, only : nqtot, moderntracdef, is_recomb, noms
       use recombin_corrk_mod, only: su_recombin,        &
@@ -716,6 +716,95 @@
 !     Initialise the continuum absorption data
       if(continuum)then
-      do igas=1,ngasmx
-
+      if(generic_continuum_database)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/LMDZ.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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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 interpolate_continuum(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
+       
+      else ! generic_continuum_database tag
+        do igas=1,ngasmx
          if (igas .eq. igas_N2) then
 
@@ -771,6 +860,7 @@
          endif
 
-      enddo
-      endif
+      enddo ! igas=1,ngasmx
+      endif ! generic_continuum_database tag
+      endif ! continuum flag
 
       print*,'----------------------------------------------------'
