Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 2859)
+++ /trunk/LMDZ.GENERIC/README	(revision 2860)
@@ -1764,2 +1764,7 @@
 aerovenus3, aerovenusUV flags which may be specified in callphys.def).
 
+== 04/01/2023 == MT
+- Add the possibility to include CO2-H2 CIA opacity (based on Turbet et al. 2020, Icarus, Volume 346, article id. 113762)
+- Add interpolateCO2H2 routine + updated associated RT routines
+(Be careful, if you use the CO2-H2 CIA, to use the associated correlated-k table, with CO2 and H2, and not just CO2)
+
Index: /trunk/LMDZ.GENERIC/libf/phystd/interpolateCO2H2.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/interpolateCO2H2.F90	(revision 2860)
+++ /trunk/LMDZ.GENERIC/libf/phystd/interpolateCO2H2.F90	(revision 2860)
@@ -0,0 +1,140 @@
+subroutine interpolateCO2H2(wn,temp,presCO2,presH2,abcoef,firstcall,ind)
+
+  !==================================================================
+  !     
+  !     Purpose
+  !     -------
+  !     Calculates the CO2-H2 CIA opacity, using a lookup table from
+  !     Turbet et al. 2020, Icarus, Volume 346, article id. 113762
+  !
+  !     Authors
+  !     -------
+  !     M. Turbet (2011)
+  !     
+  !==================================================================
+
+  use datafile_mod, only: datadir
+  implicit none
+
+  ! input
+  double precision wn                 ! wavenumber             (cm^-1)
+  double precision temp               ! temperature            (Kelvin)
+  double precision presCO2            ! CO2 partial pressure    (Pascals)
+  double precision presH2             ! H2 partial pressure    (Pascals)
+
+  ! output
+  double precision abcoef             ! absorption coefficient (m^-1)
+
+  integer nS,nT
+  parameter(nS=300)
+  parameter(nT=6)
+  double precision, parameter :: losch = 2.6867774e19
+  ! Loschmit's number (molecule cm^-3 at STP) 
+  ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2
+  ! see Richard et al. 2011, JQSRT for details
+
+  double precision amagatCO2
+  double precision amagatH2
+  double precision wn_arr(nS)
+  double precision temp_arr(nT)
+  double precision abs_arr(nS,nT)
+
+  integer k,iT
+  logical firstcall
+
+  save wn_arr, temp_arr, abs_arr !read by master
+
+  character*100 dt_file
+  integer strlen,ios
+
+  character(LEN=*), parameter :: fmat1 = "(A21,F10.3,F10.3,I7,F7.1,E10.3,F5.3)"
+
+  character*21 bleh
+  double precision blah, Ttemp
+  integer nres
+  integer ind
+
+  if(temp.gt.600)then
+     print*,'Your temperatures are too high for this CO2-H2 CIA dataset.'
+     print*,'Please run mixed CO2-H2 atmospheres below T = 600 K.'      
+     stop
+  endif
+  
+  if(temp.lt.100)then
+     print*,'Your temperatures are too low for this CO2-H2 CIA dataset.'
+     print*,'Please run mixed CO2-H2 atmospheres above T = 100 K.'      
+     stop
+  endif
+
+  amagatCO2 = (273.15/temp)*(presCO2/101325.0)
+  amagatH2 = (273.15/temp)*(presH2/101325.0)
+
+  if(firstcall)then ! called by sugas_corrk only
+     print*,'----------------------------------------------------'
+     print*,'Initialising CO2-H2 continuum from HITRAN database...'
+
+     !     1.1 Open the ASCII files
+     dt_file=TRIM(datadir)//'/continuum_data/CO2-H2_TURBET2020.cia'
+
+
+!$OMP MASTER
+     open(33,file=dt_file,form='formatted',status='old',iostat=ios)
+     if (ios.ne.0) then        ! file not found
+        write(*,*) 'Error from interpolateCO2H2'
+        write(*,*) 'Data file ',trim(dt_file),' not found.'
+        write(*,*) 'Check that your path to datagcm:',trim(datadir)
+        write(*,*) 'is correct. You can change it in callphys.def with:'
+        write(*,*) 'datadir = /absolute/path/to/datagcm'
+        write(*,*) 'Also check that the continuum data continuum_data/CO2-H2_TURBET2020.cia is there.'       
+        call abort
+     else
+
+        do iT=1,nT
+
+           read(33,fmat1) bleh,blah,blah,nres,Ttemp
+           if(nS.ne.nres)then
+              print*,'Resolution given in file: ',trim(dt_file)
+              print*,'is ',nres,', which does not match nS.'
+              print*,'Please adjust nS value in interpolateCO2H2.F90'
+              stop
+           endif
+           temp_arr(iT)=Ttemp
+
+           do k=1,nS
+              read(33,*) wn_arr(k),abs_arr(k,it)
+              write(*,*) 'for k = ', k, ' we have ', wn_arr(k),abs_arr(k,it)
+           end do
+
+        end do
+
+     endif
+     close(33)
+!$OMP END MASTER
+!$OMP BARRIER
+
+     print*,'interpolateCO2H2: At wavenumber ',wn,' cm^-1'
+     print*,'   temperature                 ',temp,' K'
+     print*,'   CO2 partial pressure         ',presCO2,' Pa'
+     print*,'   and H2 partial pressure     ',presH2,' Pa'
+
+  endif
+
+  call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind)
+
+!     print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
+!     print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
+
+  abcoef=abcoef*100.0*amagatCO2*amagatH2 ! convert to m^-1
+
+!     print*,'We have ',amagatCO2,' amagats of CO2'
+!     print*,'and     ',amagatH2,' amagats of H2'
+!     print*,'So the absorption is ',abcoef,' m^-1'
+
+
+!  unlike for Rayleigh scattering, we do not currently weight by the BB function
+!  however our bands are normally thin, so this is no big deal.
+
+
+  return
+end subroutine interpolateCO2H2
+
Index: /trunk/LMDZ.GENERIC/libf/phystd/optci.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 2859)
+++ /trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 2860)
@@ -200,4 +200,8 @@
                        call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
                        indi(nw,igas,jgas) = interm
+                    elseif(jgas.eq.igas_CO2)then 
+                       interm = indi(nw,igas,jgas)
+                       call interpolateCO2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
+                       indi(nw,igas,jgas) = interm
                     elseif(jgas.eq.igas_He)then 
                        interm = indi(nw,igas,jgas)
Index: /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 2859)
+++ /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 2860)
@@ -208,4 +208,9 @@
                        indv(nw,igas,jgas) = interm
                        ! should be irrelevant in the visible
+                    elseif(jgas.eq.igas_CO2)then 
+                       interm = indv(nw,igas,jgas)
+                       call interpolateCO2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
+                       indv(nw,igas,jgas) = interm
+                       ! might not be relevant in the visible
                     elseif(jgas.eq.igas_He)then 
                        interm = indv(nw,igas,jgas)
Index: /trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 2859)
+++ /trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90	(revision 2860)
@@ -732,4 +732,7 @@
                   dummy = -9999
                   call interpolateN2H2(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.,dummy)
+               elseif (jgas .eq. igas_CO2) then
+                  dummy = -9999
+                  call interpolateCO2H2(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.,dummy)
                elseif (jgas .eq. igas_He) then
                   dummy = -9999
