Index: /trunk/LMDZ.GENERIC/changelog.txt
===================================================================
--- /trunk/LMDZ.GENERIC/changelog.txt	(revision 3232)
+++ /trunk/LMDZ.GENERIC/changelog.txt	(revision 3233)
@@ -1855,2 +1855,5 @@
 the reference ones from https://web.lmd.jussieu.fr/~lmdz/planets/
 
+== 23/02/2024 == GM
+Add the possibility to use a fixed vertical molar fraction profile for the 
+collision-induced absorption like the correlated-k.
Index: /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 3232)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 3233)
@@ -14,5 +14,5 @@
           int_dtaui,int_dtauv,                                 &
           tau_col,cloudfrac,totcloudfrac,                      &
-          clearsky,firstcall,lastcall)
+          clearsky,p_var,frac_var,firstcall,lastcall)
 
       use mod_phys_lmdz_para, only : is_master
@@ -40,5 +40,5 @@
                               diagdtau,kastprof,strictboundcorrk,specOLR, &
                               CLFvarying,tplanckmin,tplanckmax,global1d, &
-                              generic_condensation, aerovenus
+                              generic_condensation, aerovenus, nvarlayer, varspec
       use optcv_mod, only: optcv
       use optci_mod, only: optci
@@ -46,4 +46,5 @@
       use sfluxv_mod, only: sfluxv
       use recombin_corrk_mod, only: corrk_recombin, call_recombin
+      use pindex_mod, only: pindex
       use generic_cloud_common_h, only: Psat_generic, epsi_generic
       use generic_tracer_index_mod, only: generic_tracer_index
@@ -90,4 +91,6 @@
       REAL,INTENT(IN) :: muvar(ngrid,nlayer+1)
       REAL,INTENT(IN) :: cloudfrac(ngrid,nlayer)   ! Fraction of clouds (%).
+      REAL,INTENT(IN) :: frac_var(nvarlayer,ngasmx)! Variable molar fraction.
+      REAL,INTENT(IN) :: p_var(nvarlayer)          ! Pressure for frac_var interpolation (Pa)
       logical,intent(in) :: clearsky
       logical,intent(in) :: firstcall              ! Signals first call to physics.
@@ -227,4 +230,7 @@
 
       REAL :: maxvalue,minvalue
+      
+      real :: frac_vari(L_LEVELS)
+      real :: fracvari(ngasmx,L_LEVELS)
 
 !===============================================================
@@ -1101,4 +1107,15 @@
          endif
       enddo
+      
+!-----------------------------------------------------------------------
+!     Variation of molar fraction for CIAs
+!-----------------------------------------------------------------------
+
+      if (varspec) then
+        do k=1,ngasmx
+          call pindex(p_var,frac_var(:,k),plevrad(:),nvarlayer,L_LEVELS,frac_vari)
+          fracvari(k,:) = frac_vari
+        enddo
+      endif
 
 !=======================================================================
@@ -1133,5 +1150,5 @@
             call optcv(dtauv,tauv,taucumv,plevrad,                 &
                  qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,   &
-                 tmid,pmid,taugsurf,qvar,muvarrad)
+                 tmid,pmid,taugsurf,qvar,muvarrad,fracvari)
 
             call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
@@ -1176,5 +1193,5 @@
          call optci(plevrad,tlevrad,dtaui,taucumi,                  &
               qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,    &
-              taugsurfi,qvar,muvarrad)
+              taugsurfi,qvar,muvarrad,fracvari)
 
          call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
Index: /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 3232)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 3233)
@@ -27,5 +27,8 @@
       logical,save :: testradtimes
       logical,save :: rayleigh
-!$OMP THREADPRIVATE(newtonian,force_cpp,cpp_mugaz_mode,testradtimes,rayleigh)
+      logical,save :: varspec
+      character(64),save :: varspec_data
+      integer,save :: nvarlayer
+!$OMP THREADPRIVATE(newtonian,force_cpp,cpp_mugaz_mode,testradtimes,rayleigh,varspec,varspec_data,nvarlayer)
       logical,save :: stelbbody
       logical,save :: nearco2cond
Index: /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/kcm1d.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/kcm1d.F90	(revision 3232)
+++ /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/kcm1d.F90	(revision 3233)
@@ -8,5 +8,6 @@
   use comsaison_h, only: mu0, fract, dist_star
   use planete_mod
-  use callkeys_mod, only: pceil, tstrat, tracer, global1d
+  use callkeys_mod, only: pceil, tstrat, tracer, global1d, &
+                          varspec, varspec_data, nvarlayer
   use inifis_mod, only: inifis
   use aerosol_mod, only: iniaerosol
@@ -20,4 +21,5 @@
   use geometry_mod, only: init_geometry
   use dimphy, only : init_dimphy
+  use gases_h, only: ngasmx
 
   implicit none
@@ -87,4 +89,9 @@
   real int_dtauv(1,llm,L_NSPECTV)
   real Eatmtot
+  
+  ! For fixed variable molar mass
+  real, dimension(:),allocatable,save :: p_var
+  real, dimension(:),allocatable,save :: mu_var
+  real, dimension(:,:),allocatable,save :: frac_var
 
   integer ierr
@@ -116,4 +123,32 @@
   ALLOCATE(fract(1))
   ALLOCATE(glat(1))
+
+! Initialize fixed variable mu
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
+
+  if(varspec) then
+    IF (.NOT.ALLOCATED(p_var))  ALLOCATE(p_var(nvarlayer))
+    IF (.NOT.ALLOCATED(mu_var))  ALLOCATE(mu_var(nvarlayer))
+    IF (.NOT.ALLOCATED(frac_var))  ALLOCATE(frac_var(nvarlayer,ngasmx))
+    p_var = 0.0
+    mu_var = 0.0
+    frac_var = 0.0
+    dt_file=TRIM(varspec_data)
+    open(34,file=dt_file,form='formatted',status='old',iostat=ios)
+    if (ios.ne.0) then        ! file not found
+      write(*,*) 'Error from varspec'
+      write(*,*) 'Data file ',trim(varspec_data),' not found.'
+      write(*,*) 'Check that the data is in your run repository.'
+      call abort_physic !a verifier
+    else
+      do k=1,nvarlayer
+        read(34,*) p_var(k), mu_var(k),frac_var(k,1:ngasmx)
+        !The order of columns in frac_var must correspond to order of molecules gases.def 
+        !The format of your file must be:
+        ! pressure(k) molar_mass(k), molar_fraction_of_gas_1(k), molar_fraction_of_gas_2(k),..., molar_fraction_of_gas_ngasmx(k)
+      enddo
+    endif
+    close(34)
+  endif
 
   !  Load parameters from "run.def" 
@@ -386,5 +421,5 @@
           int_dtaui,int_dtauv,                            &
           tau_col,cloudfrac,totcloudfrac,                 &
-          .false.,firstcall,lastcall)
+          .false.,p_var,frac_var,firstcall,lastcall)
 
      !write(*,*) 'BASE 3'
@@ -438,5 +473,5 @@
        int_dtaui,int_dtauv,                                    &
        tau_col,cloudfrac,totcloudfrac,                         &
-       .false.,firstcall,lastcall)
+       .false.,p_var,frac_var,firstcall,lastcall)
 
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 3232)
+++ /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 3233)
@@ -1061,4 +1061,11 @@
      call getin_p("kmixmin",kmixmin)
      if (is_master) write(*,*)trim(rname)//": kmixmin = ",kmixmin
+     
+     varspec=.false. !default value
+     call getin_p("varspec",varspec)
+     if (varspec) then
+       call getin_p("varspec_data",varspec_data)
+       call getin_p("nvarlayer",nvarlayer)
+     endif
 
      if (is_master) write(*,*)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
Index: /trunk/LMDZ.GENERIC/libf/phystd/optci.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 3232)
+++ /trunk/LMDZ.GENERIC/libf/phystd/optci.F90	(revision 3233)
@@ -7,5 +7,5 @@
 subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI,      &
      QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO,  &
-     TMID,PMID,TAUGSURF,QVAR,MUVAR)
+     TMID,PMID,TAUGSURF,QVAR,MUVAR,FRACVAR)
 
   use radinc_h, only: L_LEVELS, L_NLAYRAD, L_NSPECTI, L_NGAUSS, &
@@ -15,5 +15,5 @@
                      igas_CH4, igas_CO2
   use comcstfi_mod, only: g, r, mugaz
-  use callkeys_mod, only: kastprof,continuum,graybody
+  use callkeys_mod, only: kastprof,continuum,graybody,varspec
   use recombin_corrk_mod, only: corrk_recombin, gasi_recomb
   use tpindex_mod, only: tpindex
@@ -78,4 +78,5 @@
   real*8,intent(in) :: QVAR(L_LEVELS)
   real*8,intent(in) :: MUVAR(L_LEVELS)
+  real*8,intent(in) ::  FRACVAR(ngasmx,L_LEVELS)
   real*8  WRATIO(L_LEVELS)
   real*8  KCOEF(4)
@@ -193,4 +194,6 @@
               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)))
@@ -213,5 +216,9 @@
                  ! then cross-interactions with other gases
                  do jgas=1,ngasmx
-                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
+                    if(varspec) then
+                      p_cross = dble(PMID(k)*scalep*FRACVAR(jgas,k)*(1.-QVAR(k)))
+                    else
+                      p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
+                    endif
                     dtempc  = 0.0d0
                     if(jgas.eq.igas_N2)then 
@@ -240,5 +247,9 @@
                  ! then cross-interactions with other gases
                  do jgas=1,ngasmx
-                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
+                    if(varspec) then
+                      p_cross = dble(PMID(k)*scalep*FRACVAR(jgas,k)*(1.-QVAR(k)))
+                    else
+                      p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
+                    endif
                     dtempc  = 0.0d0
                     if(jgas.eq.igas_H2)then 
Index: /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 3232)
+++ /trunk/LMDZ.GENERIC/libf/phystd/optcv.F90	(revision 3233)
@@ -7,5 +7,5 @@
 SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,  &
      QXVAER,QSVAER,GVAER,WBARV,COSBV,       &
-     TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
+     TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR,FRACVAR)
 
   use radinc_h, only: L_NLAYRAD, L_NLEVRAD, L_LEVELS, L_NSPECTV, L_NGAUSS, L_REFVAR, NAERKIND
@@ -14,5 +14,5 @@
                      igas_CH4, igas_CO2
   use comcstfi_mod, only: g, r, mugaz
-  use callkeys_mod, only: kastprof,continuum,graybody,callgasvis
+  use callkeys_mod, only: kastprof,continuum,graybody,callgasvis,varspec
   use recombin_corrk_mod, only: corrk_recombin, gasv_recomb
   use tpindex_mod, only: tpindex
@@ -84,4 +84,5 @@
   real*8,intent(in) :: QVAR(L_LEVELS)
   real*8,intent(in) :: MUVAR(L_LEVELS)
+  real*8,intent(in) :: FRACVAR(ngasmx,L_LEVELS)
   real*8 :: WRATIO(L_LEVELS)
   real*8  KCOEF(4)
@@ -198,4 +199,6 @@
               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)))
@@ -219,5 +222,9 @@
                  ! then cross-interactions with other gases
                  do jgas=1,ngasmx
-                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
+                    if(varspec) then
+                      p_cross = dble(PMID(k)*scalep*FRACVAR(jgas,k)*(1.-QVAR(k)))
+                    else
+                      p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
+                    endif
                     dtempc  = 0.0
                     if(jgas.eq.igas_N2)then 
@@ -248,5 +255,9 @@
                  ! then cross-interactions with other gases
                  do jgas=1,ngasmx
-                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
+                    if(varspec) then
+                      p_cross = dble(PMID(k)*scalep*FRACVAR(jgas,k)*(1.-QVAR(k)))
+                    else
+                      p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
+                    endif
                     dtempc  = 0.0d0
                     if(jgas.eq.igas_H2)then 
Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 3232)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 3233)
@@ -21,5 +21,5 @@
       use generic_cloud_common_h, only : epsi_generic, Psat_generic
       use thermcell_mod, only: init_thermcell_mod
-      use gases_h, only: gnom, gfrac
+      use gases_h, only: gnom, gfrac, ngasmx
       use radcommon_h, only: sigma, glat, grav, BWNV, WNOI, DWNI, DWNV, WNOV
       use suaer_corrk_mod, only: suaer_corrk
@@ -69,5 +69,6 @@
                               startphy_file, testradtimes, tlocked, &
                               tracer, UseTurbDiff, water, watercond, &
-                              waterrain, generic_rain, global1d, moistadjustment, szangle
+                              waterrain, generic_rain, global1d, szangle, &
+                              moistadjustment, varspec, varspec_data, nvarlayer
       use generic_tracer_index_mod, only: generic_tracer_index
       use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
@@ -76,4 +77,5 @@
       use phys_state_var_mod
       use callcorrk_mod, only: callcorrk
+      use pindex_mod, only: pindex
       use vdifc_mod, only: vdifc
       use turbdiff_mod, only: turbdiff
@@ -450,4 +452,14 @@
 
       real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
+      
+      ! For fixed variable molar mass
+      real, dimension(:),allocatable,save :: p_var
+      real, dimension(:),allocatable,save :: mu_var
+      real, dimension(:,:),allocatable,save :: frac_var
+!$OMP THREADPRIVATE(p_var,mu_var,frac_var)
+      real :: mu_vari(nlayer)
+      real :: muvari(ngrid,nlayer)
+      integer ios,k
+      character*100 dt_file
 
       real,save,allocatable :: reffcol(:,:)
@@ -506,4 +518,33 @@
          zdtsw(:,:) = 0.0
          zdtlw(:,:) = 0.0
+         
+!        Initialize fixed variable mu
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
+
+         if(varspec) then
+           IF (.NOT.ALLOCATED(p_var))  ALLOCATE(p_var(nvarlayer))
+           IF (.NOT.ALLOCATED(mu_var))  ALLOCATE(mu_var(nvarlayer))
+           IF (.NOT.ALLOCATED(frac_var))  ALLOCATE(frac_var(nvarlayer,ngasmx))
+           p_var = 0.0
+           mu_var = 0.0
+           frac_var = 0.0
+           muvari = 0.0 
+           dt_file=TRIM(varspec_data)
+           open(33,file=dt_file,form='formatted',status='old',iostat=ios)
+           if (ios.ne.0) then        ! file not found
+             write(*,*) 'Error from varspec'
+             write(*,*) 'Data file ',trim(varspec_data),' not found.'
+             write(*,*) 'Check that the data is in your run repository.'
+             call abort_physic !a verifier
+           else
+             do k=1,nvarlayer
+               read(33,*) p_var(k), mu_var(k),frac_var(k,1:ngasmx)
+               !The order of columns in frac_var must correspond to order of molecules gases.def 
+               !The format of your file must be:
+               ! pressure(k) molar_mass(k), molar_fraction_of_gas_1(k), molar_fraction_of_gas_2(k),..., molar_fraction_of_gas_ngasmx(k)
+             enddo
+           endif
+           close(33)
+         endif
 
 !        Initialize tracer names, indexes and properties.
@@ -952,8 +993,9 @@
                endif
                if(water) then
+                  ! take into account water effect on mean molecular weight
                   muvar(1:ngrid,1:nlayer)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,1:nlayer,igcm_h2o_vap))
                   muvar(1:ngrid,nlayer+1)=mugaz/(1.e0+(1.e0/epsi-1.e0)*pq(1:ngrid,nlayer,igcm_h2o_vap))
-                  ! take into account water effect on mean molecular weight
                elseif(generic_condensation) then
+                  ! take into account generic condensable specie (GCS) effect on mean molecular weight
                   do iq=1,nq
 
@@ -969,5 +1011,13 @@
                      endif
                   end do ! do iq=1,nq loop on tracers
-               ! take into account generic condensable specie (GCS) effect on mean molecular weight
+               elseif(varspec) then
+                  !take into account fixed variable mean molecular weight
+                  do ig=1,ngrid
+                    call pindex(p_var,mu_var,pplay(ig,:),nvarlayer,nlayer,mu_vari)
+                    muvari(ig,:) = mu_vari
+                  enddo
+
+                  muvar(1:ngrid,1:nlayer) = muvari(1:ngrid,1:nlayer)
+                  muvar(1:ngrid,nlayer+1) = muvari(1:ngrid,nlayer)
 
                else
@@ -995,5 +1045,5 @@
                               int_dtaui,int_dtauv,                                &
                               tau_col,cloudfrac,totcloudfrac,                     &
-                              clearsky,firstcall,lastcall)
+                              clearsky,p_var,frac_var,firstcall,lastcall)
 
                 !GG (feb2021): Option to "artificially" decrease the raditive time scale in
@@ -1026,5 +1076,5 @@
                                  int_dtaui1,int_dtauv1,                              &
                                  tau_col1,cloudfrac,totcloudfrac,                    &
-                                 clearsky,firstcall,lastcall)
+                                 clearsky,p_var,frac_var,firstcall,lastcall)
                   clearsky = .false.  ! just in case.
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/pindex.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/pindex.F90	(revision 3233)
+++ /trunk/LMDZ.GENERIC/libf/phystd/pindex.F90	(revision 3233)
@@ -0,0 +1,87 @@
+module pindex_mod
+
+implicit none
+
+contains
+
+subroutine pindex(pti,vti,pref,nvarlev,nlev,nv)
+
+INTEGER, INTENT(IN) :: nvarlev
+INTEGER, INTENT(IN) :: nlev
+REAL, INTENT(IN) :: pti(nvarlev)
+REAL, INTENT(IN) :: vti(nvarlev)
+REAL, INTENT(IN) :: pref(nlev)
+
+REAL, INTENT(OUT) :: nv(nlev)
+
+REAL :: a,b
+INTEGER :: n,l
+
+!==================================================================
+!     
+!     Purpose
+!     -------
+!     Interpolate only data to the given P values.
+!
+!     Inputs
+!     ------
+!     pti                - The pressure to interpolate to
+!     vti                - The variable to interpolate to
+!     pref               - The tabulated pressure grid
+!     nvarlev            - Number of levels of pti
+!     nlev               - Number of levels of pref
+!   
+!     Outputs
+!     -------
+!     nv                 - Interpolated variable
+!
+!     Local variables
+!     -------
+!     n,l                - integer
+!     m                  - coefficient for interpolation
+!     k                  - coefficient for interpolation
+!                        
+!==================================================================
+
+DO n=1,nlev
+  IF(pti(1) .LT. pti(nvarlev)) THEN
+    IF(pref(n).LE.pti(1)) THEN
+      nv(n)=vti(1)
+    ELSEIF(pref(n).GE.pti(nvarlev)) THEN
+      nv(n)=vti(nvarlev)
+    ELSE
+      DO l=2,nvarlev
+        IF(pti(l-1).LE.pref(n) .AND. pref(n).LE.pti(l)) THEN
+          IF(abs(vti(l-1)-vti(l)).LT.1e-6) THEN
+            nv(n) = (vti(l-1) + vti(l))/2.
+          ELSE
+            a = (pti(l)-pti(l-1))/(vti(l)-vti(l-1))
+            b = (vti(l)*pti(l-1)-vti(l-1)*pti(l))/(vti(l)-vti(l-1))
+            nv(n) = (pref(n)-b)/a
+          ENDIF
+        ENDIF
+      ENDDO
+    END IF
+  ELSE
+    IF(pref(n).LE.pti(nvarlev)) THEN
+      nv(n)=vti(nvarlev)
+    ELSEIF(pref(n).GE.pti(1)) THEN
+      nv(n)=vti(1)
+    ELSE
+      DO l=2,nvarlev
+        IF(pti(l-1).GE.pref(n) .AND. pref(n).GE.pti(l)) THEN
+          IF(abs(vti(l-1)-vti(l)).LT.1e-6) THEN
+            nv(n) = (vti(l-1) + vti(l))/2.
+          ELSE
+            a = (pti(l)-pti(l-1))/(vti(l)-vti(l-1))
+            b = (vti(l)*pti(l-1)-vti(l-1)*pti(l))/(vti(l)-vti(l-1))
+            nv(n) = (pref(n)-b)/a
+          ENDIF
+        ENDIF
+      ENDDO
+    END IF
+  END IF
+ENDDO
+
+end subroutine pindex
+end module pindex_mod
