Index: /trunk/LMDZ.VENUS/libf/phyvenus/interpolate_continuum.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/interpolate_continuum.F90	(revision 3794)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/interpolate_continuum.F90	(revision 3794)
@@ -0,0 +1,833 @@
+module interpolate_continuum_mod
+
+implicit none
+
+contains
+
+     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/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,intent(in) :: ind_WN            ! wavenumber index 
+      integer,intent(in) :: igas_X            ! index of molecule X
+      integer,intent(in) :: igas_Y            ! index of molecule Y
+      double precision,intent(in) :: temp     ! temperature (Kelvin)
+      double precision,intent(in) :: pres_X   ! partial pressure of molecule X (Pascals)
+      double precision,intent(in) :: pres_Y   ! partial pressure of molecule Y (Pascals)
+      character(len=*),intent(in) :: filename ! name of the lookup table
+      character(len=2),intent(in) :: c_WN     ! wavelength chanel: infrared (IR) or visible (VI)
+      logical,intent(in) :: firstcall
+
+      ! output
+      double precision,intent(out) :: 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,save :: num_T_N2N2
+      double precision,save,dimension(:),allocatable :: temp_arr_N2N2
+      double precision,save,dimension(:,:),allocatable :: abs_arr_N2N2_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_N2N2_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair O2-O2
+      integer,save :: num_T_O2O2
+      double precision,save,dimension(:),allocatable :: temp_arr_O2O2
+      double precision,save,dimension(:,:),allocatable :: abs_arr_O2O2_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_O2O2_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair H2-H2
+      integer,save :: num_T_H2H2
+      double precision,save,dimension(:),allocatable :: temp_arr_H2H2
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2H2_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2H2_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair CO2-CO2
+      integer,save :: num_T_CO2CO2
+      double precision,save,dimension(:),allocatable :: temp_arr_CO2CO2
+      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CO2_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CO2_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair CH4-CH4
+      integer,save :: num_T_CH4CH4
+      double precision,save,dimension(:),allocatable :: temp_arr_CH4CH4
+      double precision,save,dimension(:,:),allocatable :: abs_arr_CH4CH4_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_CH4CH4_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair H2O-H2O
+      integer,save :: num_T_H2OH2O
+      double precision,save,dimension(:),allocatable :: temp_arr_H2OH2O
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2OH2O_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2OH2O_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair H2-He
+      integer,save :: num_T_H2He
+      double precision,save,dimension(:),allocatable :: temp_arr_H2He
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2He_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2He_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair H2-CH4
+      integer,save :: num_T_H2CH4
+      double precision,save,dimension(:),allocatable :: temp_arr_H2CH4
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2CH4_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2CH4_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair CO2-H2
+      integer,save :: num_T_CO2H2
+      double precision,save,dimension(:),allocatable :: temp_arr_CO2H2
+      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2H2_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2H2_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair CO2-CH4
+      integer,save :: num_T_CO2CH4
+      double precision,save,dimension(:),allocatable :: temp_arr_CO2CH4
+      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CH4_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CH4_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair N2-H2
+      integer,save :: num_T_N2H2
+      double precision,save,dimension(:),allocatable :: temp_arr_N2H2
+      double precision,save,dimension(:,:),allocatable :: abs_arr_N2H2_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_N2H2_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair N2-CH4
+      integer,save :: num_T_N2CH4
+      double precision,save,dimension(:),allocatable :: temp_arr_N2CH4
+      double precision,save,dimension(:,:),allocatable :: abs_arr_N2CH4_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_N2CH4_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair CO2-O2
+      integer,save :: num_T_CO2O2
+      double precision,save,dimension(:),allocatable :: temp_arr_CO2O2
+      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2O2_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_CO2O2_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair N2-O2
+      integer,save :: num_T_N2O2
+      double precision,save,dimension(:), allocatable :: temp_arr_N2O2
+      double precision,save,dimension(:,:), allocatable :: abs_arr_N2O2_IR
+      double precision,save,dimension(:,:), allocatable :: abs_arr_N2O2_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair H2O-N2
+      integer,save :: num_T_H2ON2
+      double precision,save,dimension(:),allocatable :: temp_arr_H2ON2
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2ON2_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2ON2_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair H2O-O2
+      integer,save :: num_T_H2OO2
+      double precision,save,dimension(:),allocatable :: temp_arr_H2OO2
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2OO2_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2OO2_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+      ! Temperature array, continuum absorption grid for the pair H2O-CO2
+      integer,save :: num_T_H2OCO2
+      double precision,save,dimension(:),allocatable :: temp_arr_H2OCO2
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2OCO2_IR
+      double precision,save,dimension(:,:),allocatable :: abs_arr_H2OCO2_VI
+! None of these saved variables are THREADPRIVATE because read by master
+! and then only accessed but never modified and thus can be shared
+
+
+      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/generic/datagcm/continuum_data/'
+          endif
+          call abort_physic("interpolate_continuum","missing input file",1)
+        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)
+
+        allocate(temp_arr(num_T)) 
+        allocate(wn_arr(num_wn)) 
+        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
+          allocate(temp_arr_CO2CO2(num_T_CO2CO2)) 
+          allocate(abs_arr_CO2CO2_VI(L_NSPECTV,num_T_CO2CO2)) 
+          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
+          allocate(temp_arr_N2N2(num_T_N2N2)) 
+          allocate(abs_arr_N2N2_VI(L_NSPECTV,num_T_N2N2)) 
+          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
+          allocate(temp_arr_O2O2(num_T_O2O2)) 
+          allocate(abs_arr_O2O2_VI(L_NSPECTV,num_T_O2O2)) 
+          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
+          allocate(temp_arr_CH4CH4(num_T_CH4CH4)) 
+          allocate(abs_arr_CH4CH4_VI(L_NSPECTV,num_T_CH4CH4)) 
+          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
+          allocate(temp_arr_H2H2(num_T_H2H2)) 
+          allocate(abs_arr_H2H2_VI(L_NSPECTV,num_T_H2H2)) 
+          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
+          allocate(temp_arr_H2OH2O(num_T_H2OH2O)) 
+          allocate(abs_arr_H2OH2O_VI(L_NSPECTV,num_T_H2OH2O)) 
+          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
+          allocate(temp_arr_N2H2(num_T_N2H2)) 
+          allocate(abs_arr_N2H2_VI(L_NSPECTV,num_T_N2H2)) 
+          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
+          allocate(temp_arr_N2O2(num_T_N2O2)) 
+          allocate(abs_arr_N2O2_VI(L_NSPECTV,num_T_N2O2)) 
+          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
+          allocate(temp_arr_N2CH4(num_T_N2CH4)) 
+          allocate(abs_arr_N2CH4_VI(L_NSPECTV,num_T_N2CH4)) 
+          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
+          allocate(temp_arr_CO2O2(num_T_CO2O2)) 
+          allocate(abs_arr_CO2O2_VI(L_NSPECTV,num_T_CO2O2)) 
+          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
+          allocate(temp_arr_H2CH4(num_T_H2CH4)) 
+          allocate(abs_arr_H2CH4_VI(L_NSPECTV,num_T_H2CH4)) 
+          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
+          allocate(temp_arr_H2He(num_T_H2He)) 
+          allocate(abs_arr_H2He_VI(L_NSPECTV,num_T_H2He)) 
+          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
+          allocate(temp_arr_H2ON2(num_T_H2ON2)) 
+          allocate(abs_arr_H2ON2_VI(L_NSPECTV,num_T_H2ON2)) 
+          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
+          allocate(temp_arr_H2OO2(num_T_H2OO2)) 
+          allocate(abs_arr_H2OO2_VI(L_NSPECTV,num_T_H2OO2)) 
+          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
+          allocate(temp_arr_H2OCO2(num_T_H2OCO2)) 
+          allocate(abs_arr_H2OCO2_VI(L_NSPECTV,num_T_H2OCO2)) 
+          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
+          allocate(temp_arr_CO2CO2(num_T_CO2CO2)) 
+          allocate(abs_arr_CO2CO2_VI(L_NSPECTV,num_T_CO2CO2)) 
+          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
+          allocate(temp_arr_CO2H2(num_T_CO2H2)) 
+          allocate(abs_arr_CO2H2_VI(L_NSPECTV,num_T_CO2H2)) 
+          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
+          allocate(temp_arr_CO2CH4(num_T_CO2CH4)) 
+          allocate(abs_arr_CO2CH4_VI(L_NSPECTV,num_T_CO2CH4)) 
+          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
+
+      !!!! TEST !!!
+      !abs_coef=abs_coef*30
+
+      !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'
+      
+    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,intent(in) :: num_T
+      integer,intent(in) :: num_wn
+      double precision,intent(in) :: wn_arr(num_wn)
+      double precision,intent(in) :: abs_arr_in(num_wn,num_T)
+      double precision,intent(out) :: abs_arr_out_IR(L_NSPECTI,num_T)
+      double precision,intent(out) :: 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 infrared (IR) 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
+	
+    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,intent(out) :: z_temp
+      double precision,intent(in) :: temp
+      integer,intent(in) :: num_T
+      double precision,intent(in) :: temp_arr(num_T)
+      
+      character(len=22) :: rname = "T_boundaries_continuum"
+      
+      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
+          call abort_physic(rname,"temperature too low",1)
+!        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
+          call abort_physic(rname,"temperature too low",1)
+!        else
+!          z_temp=maxval(temp_arr)
+!        endif
+      endif
+      
+    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'
+
+      
+    end subroutine interpolate_T_abs_coeff
+
+end module interpolate_continuum_mod
Index: /trunk/LMDZ.VENUS/libf/phyvenus/optcv.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/optcv.F90	(revision 3793)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/optcv.F90	(revision 3794)
@@ -11,5 +11,6 @@
   use radinc_h, only: L_NLAYRAD, L_NLEVRAD, L_LEVELS, L_NSPECTV, L_NGAUSS, L_REFVAR, NAERKIND
   use radcommon_h, only: gasv, tlimit, tgasref, pfgasref,wnov,scalep,indv
-  use gases_h, only: gfrac, ngasmx, igas_H2, igas_H2O, igas_He, igas_N2
+  use gases_h, only: gfrac, ngasmx, igas_H2O, igas_CO2
+  use interpolate_continuum_mod, only: interpolate_continuum
 
   implicit none
@@ -72,5 +73,4 @@
   real*8 DRAYAER
   double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
-  double precision p_cross
 
   ! variable species mixing ratio variables
@@ -166,5 +166,4 @@
 !        if(continuum)then
            ! include continua if necessary
-           wn_cont = dble(wnov(nw))
            T_cont  = dble(TMID(k))
            do igas=1,ngasmx
@@ -178,17 +177,12 @@
               dtemp=0.0
 
-! For Venus: only H2O, using CKD
-              if(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
-
-                 p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
-!                 if(H2Ocont_simple)then
-!                    call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
-!                 else
-                    interm = indv(nw,igas,igas)
-                    call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm)
-                    indv(nw,igas,igas) = interm
-!                 endif
-
-              endif
+! For Venus: only H2O and CO2
+              do jgas=1,ngasmx
+                if ( ((igas .eq. igas_H2O) .and. (jgas .eq. igas_H2O)) .or.  &
+                     ((igas .eq. igas_H2O) .and. (jgas .eq. igas_CO2)) .or.   &
+                     ((igas .eq. igas_CO2) .and. (jgas .eq. igas_CO2)) ) then
+                  call interpolate_continuum('',igas,jgas,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
+                endif
+              enddo
 
               DCONT = DCONT + dtemp
Index: /trunk/LMDZ.VENUS/libf/phyvenus/radcommon_h.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/radcommon_h.F90	(revision 3793)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/radcommon_h.F90	(revision 3794)
@@ -1,4 +1,4 @@
 module radcommon_h
-      use radinc_h, only: L_NSPECTV, NTstart, NTstop, &
+      use radinc_h, only: L_NSPECTI, L_NSPECTV, NTstart, NTstop, &
                           naerkind, nsizemax
       implicit none
@@ -13,4 +13,13 @@
 !  some or all of this common data set
 !
+!     WNOI       - Array of wavenumbers at the spectral interval
+!                  centers for the infrared.  Array is NSPECTI
+!                  elements long.
+!     DWNI       - Array of "delta wavenumber", i.e., the width,
+!                  in wavenumbers (cm^-1) of each IR spectral
+!                  interval.  NSPECTI elements long.
+!     WAVEI      - Array (NSPECTI elements long) of the wavelenght
+!                  (in microns) at the center of each IR spectral
+!                  interval.
 !     WNOV       - Array of wavenumbers at the spectral interval
 !                  center for the VISUAL.  Array is NSPECTV
@@ -29,4 +38,6 @@
 !     TAURAYVAR  - Array (NSPECTV elements) of the pressure-independent
 !                  part of Rayleigh scattering optical depth for the variable gas.
+!     FZEROI     - Fraction of zeros in the IR CO2 k-coefficients, for
+!                  each temperature, pressure, and spectral interval
 !     FZEROV     - Fraction of zeros in the VISUAL CO2 k-coefficients, for
 !                  each temperature, pressure, and spectral interval
@@ -42,30 +53,51 @@
 ! gvis       :  assymetry factor
 ! 
+!   Longwave
+!   ~~~~~~~~
+! 
+! For the "naerkind" kind of aerosol radiative properties : 
+! QIRsQREF :  Qext / Qext("longrefir")
+! omegaIR  :  mean single scattering albedo
+! gIR      :  mean assymetry factor
+!
+!
+! Note - QIRsQREF in the martian model was scaled to longrefvis,
+! here it is scaled to longrefir, which is actually a dummy parameter,
+! as we do not output scaled aerosol opacity ...
+!
 
+      REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) !BWNI read by master in setspi
       REAL*8 BWNV(L_NSPECTV+1), WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV) !BWNV read by master in setspv
       REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV), TAURAYVAR(L_NSPECTV)
-!$OMP THREADPRIVATE(WNOV,DWNV,WAVEV,&
+!$OMP THREADPRIVATE(WNOI,DWNI,WAVEI,&
+	!$OMP WNOV,DWNV,WAVEV,&
 	!$OMP STELLARF,TAURAY,TAURAYVAR)
 
+      REAL*8 blami(L_NSPECTI+1)
       REAL*8 blamv(L_NSPECTV+1) ! these are needed by suaer.F90
-!$OMP THREADPRIVATE(blamv)
+!$OMP THREADPRIVATE(blami,blamv)
 
       !! AS: introduced to avoid doing same computations again for continuum
+      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: indi
       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: indv
-!$OMP THREADPRIVATE(indv)
+!$OMP THREADPRIVATE(indi,indv)
 
       !!! ALLOCATABLE STUFF SO THAT DIMENSIONS ARE READ in *.dat FILES -- AS 12/2011  
-      REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasv
+      REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi, gasv
       REAL*8, DIMENSION(:), ALLOCATABLE :: PGASREF, TGASREF, PFGASREF, GWEIGHT
+      real*8 FZEROI(L_NSPECTI)
       real*8 FZEROV(L_NSPECTV)
       real*8 pgasmin, pgasmax
       real*8 tgasmin, tgasmax
-!$OMP THREADPRIVATE(gasv,&  !pgasref,tgasref,pfgasref read by master in sugas_corrk
-	!$OMP FZEROV)           !pgasmin,pgasmax,tgasmin,tgasmax read by master in sugas_corrk
+!$OMP THREADPRIVATE(gasi,gasv,&  !pgasref,tgasref,pfgasref read by master in sugas_corrk
+	!$OMP FZEROI,FZEROV)           !pgasmin,pgasmax,tgasmin,tgasmax read by master in sugas_corrk
 
       real QVISsQREF(L_NSPECTV,naerkind,nsizemax)
       real omegavis(L_NSPECTV,naerkind,nsizemax)
       real gvis(L_NSPECTV,naerkind,nsizemax)
-!$OMP THREADPRIVATE(QVISsQREF,omegavis,gvis)
+      real QIRsQREF(L_NSPECTV,naerkind,nsizemax)
+      real omegair(L_NSPECTV,naerkind,nsizemax)
+      real gir(L_NSPECTV,naerkind,nsizemax)
+!$OMP THREADPRIVATE(QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir)
 
 
@@ -73,5 +105,5 @@
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-      REAL lamrefvis(naerkind)
+      REAL lamrefir(naerkind),lamrefvis(naerkind)
 
 ! Actual number of grain size classes in each domain for a
@@ -87,16 +119,18 @@
 
 ! Extinction coefficient at reference wavelengths;
-!   These wavelengths are defined in aeroptproperties, and called longrefvis
+!   These wavelengths are defined in aeroptproperties, and called longrefvis and longrefir
 
       REAL :: QREFvis(naerkind,nsizemax)
+      REAL :: QREFir(naerkind,nsizemax)
+      REAL :: omegaREFir(naerkind,nsizemax)
 
       REAL,SAVE :: tstellar ! Stellar brightness temperature (SW)
-
+      REAL*8,SAVE :: planckir(naerkind,nsizemax)
       real*8,save :: PTOP
 
       real*8,parameter :: UBARI = 0.5D0
 
-!$OMP THREADPRIVATE(QREFvis,& 	! gweight read by master in sugas_corrk
-		!$OMP tstellar,PTOP)
+!$OMP THREADPRIVATE(QREFvis,QREFir,omegaREFir,& 	! gweight read by master in sugas_corrk
+		!$OMP tstellar,planckir,PTOP)
 
 !     If the gas optical depth (top to the surface) is less than
Index: /trunk/LMDZ.VENUS/libf/phyvenus/radinc_h.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/radinc_h.F90	(revision 3793)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/radinc_h.F90	(revision 3794)
@@ -79,4 +79,5 @@
 
       integer, parameter :: L_NSPECTV = NBvisible
+      integer, parameter :: L_NSPECTI = NBinfrared
 
       integer, parameter :: NAERKIND  = 5
Index: /trunk/LMDZ.VENUS/libf/phyvenus/sfluxv.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/sfluxv.F	(revision 3793)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/sfluxv.F	(revision 3794)
@@ -1,5 +1,5 @@
       SUBROUTINE SFLUXV(DTAUV,TAUV,TAUCUMV,RSFV,DWNV,WBARV,COSBV,
      *                  UBAR0,STEL,NFLUXTOPV,FLUXTOPVDN,
-     *                  NFLUXOUTV_nu,NFLUXGNDV_nu,
+     *                  NFLUXV_nu,NFLUXGNDV_nu,
      *                  FMNETV,FLUXUPV,FLUXDNV,FZEROV,taugsurf)
 
@@ -19,5 +19,5 @@
       real*8 FLUXUPV(L_NLAYRAD), FLUXDNV(L_NLAYRAD)
       real*8 NFLUXTOPV, FLUXUP, FLUXDN,FLUXTOPVDN
-      real*8 NFLUXOUTV_nu(L_NSPECTV)
+      real*8 NFLUXV_nu(L_NLAYRAD,L_NSPECTV)
       real*8 NFLUXGNDV_nu(L_NSPECTV)
 
@@ -40,5 +40,5 @@
 
       DO NW=1,L_NSPECTV
-         NFLUXOUTV_nu(NW)=0.0
+         NFLUXV_nu(:,NW)=0.0
          NFLUXGNDV_nu(NW)=0.0
       END DO
@@ -115,7 +115,7 @@
           END DO
 
-c     band-resolved flux leaving TOA (RDW)
-          NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW)
-     *      +FLUXUP*GWEIGHT(NG)*(1.0-FZEROV(NW))
+c     band-resolved net flux 
+          NFLUXV_nu(:,NW) = NFLUXV_nu(:,NW)
+     *      +(FMUPV-FMDV)*GWEIGHT(NG)*(1.0-FZEROV(NW))
 
 c     band-resolved flux at ground (RDW)
@@ -173,7 +173,7 @@
         END DO
 
-c     band-resolved flux leaving TOA (RDW)
-          NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW)
-     *      +FLUXUP*FZERO
+c     band-resolved net flux 
+          NFLUXV_nu(:,NW) = NFLUXV_nu(:,NW)
+     *      +(FMUPV-FMDV)*FZERO
 
 c     band-resolved flux at ground (RDW)
Index: /trunk/LMDZ.VENUS/libf/phyvenus/su_gases.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/su_gases.F90	(revision 3793)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/su_gases.F90	(revision 3794)
@@ -39,4 +39,12 @@
    igas_OCS=5
 
+  ! We don't care about the others
+  ! but initialized in order to avoid problems in interpolate_continuum.F90
+  igas_H2=99
+  igas_He=99
+  igas_N2=99
+  igas_CH4=99
+  igas_O2=99
+
   vgas=0
   if(.not.allocated(gfrac)) allocate(gfrac(ngasmx))
Index: /trunk/LMDZ.VENUS/libf/phyvenus/sugas_corrk.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/sugas_corrk.F90	(revision 3793)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/sugas_corrk.F90	(revision 3794)
@@ -27,5 +27,6 @@
       use radcommon_h, only : WNOV
       use datafile_mod, only: datadir,corrkdir,banddir
-      use gases_h, only: gnom, ngasmx, igas_H2O
+      use gases_h, only: gnom, ngasmx, igas_H2O, igas_CO2
+      use interpolate_continuum_mod, only: interpolate_continuum
       implicit none
 #include "YOMCST.h"
@@ -124,4 +125,6 @@
                  'match that in gases.def (',trim(gnom(igas)),'). You should compare ', &
                  'gases.def with Q.dat in your radiative transfer directory.' 
+            print*,'For Venus, no gases.def, but su_gases.F90 hardcoded... Should be identical ', &
+                 'to Q.dat in your radiative transfer directory.'
             message='exiting.'
             call abort_physic(subname,message,1)
@@ -452,22 +455,30 @@
 !=======================================================================
 !     Initialise the continuum absorption data
-!      if(continuum)then
-      do igas=1,ngasmx
-
-! For Venus: only H2O, using CKD
-         if (igas .eq. igas_H2O) then
-
-            ! H2O is special 
-!            if(H2Ocont_simple)then
-!               call interpolateH2Ocont_PPC(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
-!            else
-               dummy = -9999
-               call interpolateH2Ocont_CKD(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.,dummy)
-!            endif
-
-         endif
-
-      enddo
-!      endif
+!     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/
+
+! For Venus: only H2O and CO2
+        if (igas .eq. igas_H2O) then
+         file_id='/continuum_data/' // 'H2O-H2O_2022.cia'
+         ! H2O-H2O_2022.cia = 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_CO2) then
+           file_id='/continuum_data/' // 'H2O-CO2_2022.cia'
+           ! H2O-CO2_2022.cia = 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_2025.cia'
+         ! CO2-CO2_2025.cia = 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.)
+        endif
+       enddo ! igas=1,ngasmx
+!     endif
 
 !      print*,'----------------------------------------------------'
Index: /trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_corrk.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_corrk.F90	(revision 3793)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_corrk.F90	(revision 3794)
@@ -9,5 +9,5 @@
                              gweight, pfgasref, pgasmax, pgasmin, &
                              pgasref, tgasmax, tgasmin, tgasref, scalep, &
-                             stellarf, dwnv, tauray
+                             stellarf, dwnv, tauray, wavev, wnov
       use datafile_mod, only: datadir,banddir,corrkdir
       use ioipsl_getin_p_mod, only: getin_p
@@ -99,5 +99,6 @@
       REAL*8 tauaero(L_LEVELS,naerkind)
       REAL*8 nfluxtopv,fluxtopvdn
-      REAL*8 nfluxoutv_nu(L_NSPECTV)                 ! Outgoing band-resolved VI flux at TOA (W/m2).
+      REAL*8 nfluxv_nu(L_NLAYRAD,L_NSPECTV)  ! vertical profil, band-resolved VI net flux (W/m2)
+      REAL*8 heatrate_nu(L_NLAYRAD,L_NSPECTV)  ! vertical profil, band-resolved VI heating rate (K/s/micron)
       REAL*8 fmnetv(L_NLAYRAD)
       REAL*8 fluxupv(L_NLAYRAD)
@@ -504,5 +505,5 @@
             call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
                  acosz,stel_fract,                                 &
-                 nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu,   &
+                 nfluxtopv,fluxtopvdn,nfluxv_nu,nfluxgndv_nu,      &
                  fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
 
@@ -510,5 +511,5 @@
             nfluxtopv       = 0.0d0
             fluxtopvdn      = 0.0d0
-            nfluxoutv_nu(:) = 0.0d0
+            nfluxv_nu(:,:)  = 0.0d0
             nfluxgndv_nu(:) = 0.0d0
             do l=1,L_NLAYRAD
@@ -564,4 +565,5 @@
              *RG/(cpdet(tmid(1))*scalep*(plevrad(3)-plevrad(2)))
 
+
 !-----------------------------------------------------------------------    
 !-----------------------------------------------------------------------    
@@ -570,4 +572,21 @@
 !-----------------------------------------------------------------------
 
+! Diagnostic in 1D: 
+!     output the vertical profil of band-resolved heating rates, in K/s/microns
+      if ((1.eq.1).and.is_master.and.(ngrid.eq.1).and.firstcall) then
+         open(unit=11,file="nfluxnu.dat")
+         write(11,*) "lambda(microns)",wavev(:)
+         write(11,*) "dlbda(microns)",1.e4*dwnv(:)/wnov(:)**2
+         write(11,*) "pressure(Pa) dT_nu(K/s/micron)"
+         DO l=2,L_NLAYRAD
+            heatrate_nu(L_NLAYRAD+1-l,:)=(nfluxv_nu(l,:)-nfluxv_nu(l-1,:))  &
+                *RG/(cpdet(tmid(l))*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))  &
+            ! dlbda(microns)=-dnu(cm-1)*1E4/nu(cm-1)^2
+                / (1.e4*dwnv(:)/wnov(:)**2)
+            write(11,*) pplay(1,L_NLAYRAD+1-l),heatrate_nu(L_NLAYRAD+1-l,:)
+         END DO
+         close(11)
+      endif
+         
     end subroutine sw_venus_corrk
 
