Index: LMDZ6/branches/cirrus/libf/dyn3d/check_isotopes.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3d/check_isotopes.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3d/check_isotopes.F90	(revision 5202)
@@ -23,4 +23,5 @@
                              iso_O17, iso_HTO
    LOGICAL, SAVE :: first=.TRUE.
+   LOGICAL, PARAMETER :: tnat1=.TRUE.
 
    modname='check_isotopes'
@@ -34,5 +35,9 @@
       iso_O17 = strIdx(isoName,'H217O')
       iso_HTO = strIdx(isoName,'HTO')
-      IF(getKey('tnat', tnat)) CALL abort_gcm(modname, 'missing isotopic parameter', 1)
+      if (tnat1) then
+              tnat(:)=1.0
+      else
+         IF(getKey('tnat', tnat)) CALL abort_gcm(modname, 'missing isotopic parameter', 1)
+      endif
       first = .FALSE.
    END IF
Index: LMDZ6/branches/cirrus/libf/dyn3d/conf_gcm.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3d/conf_gcm.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3d/conf_gcm.F90	(revision 5202)
@@ -18,5 +18,5 @@
   USE logic_mod, ONLY: fxyhypb, iflag_phys, ok_etat0, ok_gradsfile, &
                        ok_guide, ok_limit, ok_strato, purmats, read_start, &
-                       ysinus, read_orop
+                       ysinus, read_orop, adv_qsat_liq
   USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy, &
                        alphax,alphay,taux,tauy
@@ -606,4 +606,13 @@
      type_trac = 'lmdz'
      CALL getin('type_trac',type_trac)
+
+
+     !Config  Key  = adv_qsat_liq
+     !Config  Desc = option for qsat calculation in the dynamics
+     !Config  Def  = n 
+     !Config  Help = controls which phase is considered for qsat calculation 
+     !Config          
+     adv_qsat_liq = .FALSE.
+     CALL getin('adv_qsat_liq',adv_qsat_liq)
 
      !Config  Key  = ok_dynzon 
@@ -672,4 +681,5 @@
      write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins 
      write(lunout,*)' ok_dyn_ave = ', ok_dyn_ave 
+     write(lunout,*)' adv_qsat_liq = ', adv_qsat_liq
   ELSE
      !Config  Key  = clon
Index: LMDZ6/branches/cirrus/libf/dyn3d/dynetat0.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3d/dynetat0.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3d/dynetat0.F90	(revision 5202)
@@ -43,4 +43,5 @@
   REAL    :: time, tnat, alpha_ideal, tab_cntrl(length)    !--- RUN PARAMS TABLE
   LOGICAL :: lSkip, ll
+  LOGICAL,PARAMETER :: tnat1=.TRUE.
 !-------------------------------------------------------------------------------
   modname="dynetat0"
@@ -155,6 +156,12 @@
       iqParent = tracers(iq)%iqParent
       IF(tracers(iq)%iso_iZone == 0) THEN
-         IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
+         if (tnat1) then
+                 tnat=1.0
+                 alpha_ideal=1.0
+                 write(*,*) 'attention dans dynetat0: les alpha_ideal sont a 1'
+         else
+          IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
             CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
+         endif
          CALL msg('Tracer <'//TRIM(var)//'> is missing => initialized with a simplified Rayleigh distillation law.', modname)
          q(:,:,:,iq) = q(:,:,:,iqParent)*tnat*(q(:,:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
Index: LMDZ6/branches/cirrus/libf/dyn3d/dynredem_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3d/dynredem_mod.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3d/dynredem_mod.F90	(revision 5202)
@@ -4,5 +4,5 @@
   PRIVATE
   PUBLIC :: dynredem_write_u, dynredem_write_v, dynredem_read_u, err
-  PUBLIC :: cre_var, get_var1, put_var1, put_var2, fil, modname, msg
+  PUBLIC :: cre_var, put_var1, put_var2, fil, modname, msg
   include "dimensions.h"
   include "paramet.h"
Index: LMDZ6/branches/cirrus/libf/dyn3d/iniacademic.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3d/iniacademic.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3d/iniacademic.F90	(revision 5202)
@@ -80,4 +80,5 @@
 
   REAL zdtvr, tnat, alpha_ideal
+  LOGICAL,PARAMETER :: tnat1=.true.
   
   character(len=*),parameter :: modname="iniacademic"
@@ -321,10 +322,22 @@
               iqParent = tracers(iq)%iqParent
               IF(tracers(iq)%iso_iZone == 0) THEN
-                 IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
+                 if (tnat1) then
+                         tnat=1.0
+                         alpha_ideal=1.0
+                         write(*,*) 'Attention dans iniacademic: alpha_ideal=1'
+                 else
+                    IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
                     CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
+                 endif
                  q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
-              ELSE
-                 q(:,:,iq) = q(:,:,iqIsoPha(iName,iPhase))
-              END IF
+              ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
+                 IF(tracers(iq)%iso_iZone == 1) THEN
+                    ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
+                    ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
+                    q(:,:,iq) = q(:,:,iqIsoPha(iName,iPhase))
+                 else !IF(tracers(iq)%iso_iZone == 1) THEN
+                    q(:,:,iq) = 0.
+                 endif !IF(tracers(iq)%iso_iZone == 1) THEN
+              END IF !IF(tracers(iq)%iso_iZone == 0) THEN
            enddo
         else
Index: LMDZ6/branches/cirrus/libf/dyn3d/logic_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3d/logic_mod.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3d/logic_mod.F90	(revision 5202)
@@ -30,4 +30,6 @@
   LOGICAL hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
                  ! (only used if disvert_type==2)
+  LOGICAL adv_qsat_liq ! true if qsat is calculated alwats wrt liquid for
+                       ! adapted Van Leer advection scheme
   INTEGER iflag_phys ! type of physics to call: 0 none, 1: phy*** package,
                      ! 2: Held & Suarez, 101-200: aquaplanets & terraplanets
Index: LMDZ6/branches/cirrus/libf/dyn3d/qminimum.F
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3d/qminimum.F	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3d/qminimum.F	(revision 5202)
@@ -28,4 +28,7 @@
 c     .................................................................
 c
+cDC iq_val and iq_liq are usable for q only, NOT for q_follow
+c   and zx_defau_diag (crash if iq_val/liq==3) => vapor/liquid
+c   water at hardcoded indices 1/2 in these variables
       INTEGER i, k, iq
       REAL zx_defau, zx_abc, zx_pump(ip1jmp1), pompe
@@ -58,37 +61,38 @@
 
       zx_defau_diag(:,:,:)=0.0
-      q_follow(:,:,1:2)=q(:,:,1:2)  
-      DO 1000 k = 1, llm
-        DO 1040 i = 1, ip1jmp1
+      q_follow(:,:,1)=q(:,:,iq_vap)  
+      q_follow(:,:,2)=q(:,:,iq_liq)  
+      DO k = 1, llm
+        DO i = 1, ip1jmp1
           if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then
 
-              if (niso > 0) zx_defau_diag(i,k,iq_liq)=AMAX1
+            if (niso > 0) zx_defau_diag(i,k,2)=AMAX1
      :               ( seuil_liq - q(i,k,iq_liq), 0.0 )
 
-             q(i,k,iq_vap) = q(i,k,iq_vap) + q(i,k,iq_liq) - seuil_liq
-             q(i,k,iq_liq) = seuil_liq
-           endif
- 1040   CONTINUE
- 1000 CONTINUE
+            q(i,k,iq_vap) = q(i,k,iq_vap) + q(i,k,iq_liq) - seuil_liq
+            q(i,k,iq_liq) = seuil_liq
+          endif
+        ENDDO
+      ENDDO
 c
 c Quand l'eau vapeur est trop faible (ou negative), on complete
 c le defaut en prennant de l'eau vapeur de la couche au-dessous.
 c
-      iq = iq_vap
-c
       DO k = llm, 2, -1
 ccc      zx_abc = dpres(k) / dpres(k-1)
         DO i = 1, ip1jmp1
-          if ( seuil_vap - q(i,k,iq) .gt. 0.d0 ) then
+          if ( seuil_vap - q(i,k,iq_vap) .gt. 0.d0 ) then
 
-            if (niso > 0)
-     &        zx_defau_diag(i,k,iq)=AMAX1( seuil_vap - q(i,k,iq), 0.0 )
+            if (niso > 0) zx_defau_diag(i,k,1)
+     &           = AMAX1( seuil_vap - q(i,k,iq_vap), 0.0 )
 
-            q(i,k-1,iq) =  q(i,k-1,iq) - ( seuil_vap - q(i,k,iq) ) *
-     &                     deltap(i,k) / deltap(i,k-1)
-            q(i,k,iq)   =  seuil_vap  
+            q(i,k-1,iq_vap) = q(i,k-1,iq_vap) - (seuil_vap
+     &           -q(i,k,iq_vap)) * deltap(i,k)/deltap(i,k-1)
+            q(i,k,iq_vap)   =  seuil_vap  
+
           endif
         ENDDO
       ENDDO
+
 c
 c Quand il s'agit de la premiere couche au-dessus du sol, on
@@ -96,6 +100,6 @@
 c
       DO i = 1, ip1jmp1
-         zx_pump(i) = AMAX1( 0.0, seuil_vap - q(i,1,iq) )
-         q(i,1,iq)  = AMAX1( q(i,1,iq), seuil_vap )
+         zx_pump(i) = AMAX1( 0.0, seuil_vap - q(i,1,iq_vap) )
+         q(i,1,iq_vap)  = AMAX1( q(i,1,iq_vap), seuil_vap )
       ENDDO
       pompe = SSUM(ip1jmp1,zx_pump,1)
@@ -121,5 +125,5 @@
       DO i = 1,ip1jmp1
         if (zx_pump(i).gt.0.0) then
-          q_follow(i,1,iq_vap)=q_follow(i,1,iq_vap)+zx_pump(i)
+          q_follow(i,1,1)=q_follow(i,1,1)+zx_pump(i)
         endif !if (zx_pump(i).gt.0.0) then
       enddo !DO i = 1,ip1jmp1
@@ -129,26 +133,26 @@
       do k=2,llm
         DO i = 1,ip1jmp1
-          if (zx_defau_diag(i,k,iq_vap).gt.0.0) then             
+          if (zx_defau_diag(i,k,1).gt.0.0) then             
               ! on ajoute la vapeur en k              
               do ixt=1,ntiso
                q(i,k,iqIsoPha(ixt,iq_vap))=q(i,k,iqIsoPha(ixt,iq_vap))
-     :           +zx_defau_diag(i,k,iq_vap)
-     :           *q(i,k-1,iqIsoPha(ixt,iq_vap))/q_follow(i,k-1,iq_vap)
+     :           +zx_defau_diag(i,k,1)
+     :           *q(i,k-1,iqIsoPha(ixt,iq_vap))/q_follow(i,k-1,1)
                 
               ! et on la retranche en k-1
                q(i,k-1,iqIsoPha(ixt,iq_vap))=
      :            q(i,k-1,iqIsoPha(ixt,iq_vap))
-     :              -zx_defau_diag(i,k,iq_vap)
+     :              -zx_defau_diag(i,k,1)
      :              *deltap(i,k)/deltap(i,k-1)
      :              *q(i,k-1,iqIsoPha(ixt,iq_vap))
-     :              /q_follow(i,k-1,iq_vap)
+     :              /q_follow(i,k-1,1)
 
               enddo !do ixt=1,niso
-              q_follow(i,k,iq_vap)=   q_follow(i,k,iq_vap)
-     :               +zx_defau_diag(i,k,iq_vap)
-              q_follow(i,k-1,iq_vap)=   q_follow(i,k-1,iq_vap)
-     :               -zx_defau_diag(i,k,iq_vap)
+              q_follow(i,k,1)=   q_follow(i,k,1)
+     :               +zx_defau_diag(i,k,1)
+              q_follow(i,k-1,1)=   q_follow(i,k-1,1)
+     :               -zx_defau_diag(i,k,1)
      :              *deltap(i,k)/deltap(i,k-1)
-          endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then
+          endif !if (zx_defau_diag(i,k,1).gt.0.0) then
         enddo !DO i = 1, ip1jmp1        
        enddo !do k=2,llm
@@ -161,21 +165,21 @@
         do k=1,llm
         DO i = 1,ip1jmp1
-          if (zx_defau_diag(i,k,iq_liq).gt.0.0) then
+          if (zx_defau_diag(i,k,2).gt.0.0) then
 
               ! on ajoute eau liquide en k en k              
               do ixt=1,ntiso
                q(i,k,iqIsoPha(ixt,iq_liq))=q(i,k,iqIsoPha(ixt,iq_liq))
-     :              +zx_defau_diag(i,k,iq_liq)
-     :              *q(i,k,iqIsoPha(ixt,iq_vap))/q_follow(i,k,iq_vap)
+     :              +zx_defau_diag(i,k,2)
+     :              *q(i,k,iqIsoPha(ixt,iq_vap))/q_follow(i,k,1)
               ! et on la retranche à la vapeur en k
                q(i,k,iqIsoPha(ixt,iq_vap))=q(i,k,iqIsoPha(ixt,iq_vap))
-     :              -zx_defau_diag(i,k,iq_liq)
-     :              *q(i,k,iqIsoPha(ixt,iq_vap))/q_follow(i,k,iq_vap)   
+     :              -zx_defau_diag(i,k,2)
+     :              *q(i,k,iqIsoPha(ixt,iq_vap))/q_follow(i,k,1)   
               enddo !do ixt=1,niso
-              q_follow(i,k,iq_liq)=   q_follow(i,k,iq_liq)
-     :               +zx_defau_diag(i,k,iq_liq)
-              q_follow(i,k,iq_vap)=   q_follow(i,k,iq_vap)
-     :               -zx_defau_diag(i,k,iq_liq)
-          endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then
+              q_follow(i,k,2)=   q_follow(i,k,2)
+     :               +zx_defau_diag(i,k,2)
+              q_follow(i,k,1)=   q_follow(i,k,1)
+     :               -zx_defau_diag(i,k,2)
+          endif !if (zx_defau_diag(i,k,1).gt.0.0) then
         enddo !DO i = 1, ip1jmp1
        enddo !do k=2,llm  
Index: LMDZ6/branches/cirrus/libf/dyn3d/vlspltqs.F
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3d/vlspltqs.F	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3d/vlspltqs.F	(revision 5202)
@@ -25,5 +25,5 @@
       
       USE comconst_mod, ONLY: cpp
-      
+      USE logic_mod, ONLY: adv_qsat_liq 
       IMPLICIT NONE
 c
@@ -92,5 +92,9 @@
          ENDDO
          DO ij = 1, ip1jmp1
-          zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
+          IF (adv_qsat_liq) THEN
+             zdelta = 0.
+          ELSE
+             zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
+          ENDIF
           play   = 0.5*(p(ij,l)+p(ij,l+1))
           qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
