Index: LMDZ6/branches/cirrus/libf/dyn3dmem/check_isotopes_loc.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3dmem/check_isotopes_loc.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3dmem/check_isotopes_loc.F90	(revision 5202)
@@ -24,4 +24,5 @@
                              iso_O17, iso_HTO
    LOGICAL, SAVE :: first=.TRUE.
+   LOGICAL, PARAMETER :: tnat1=.TRUE.
 !$OMP THREADPRIVATE(first)
 
@@ -37,5 +38,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
 !$OMP END MASTER
 !$OMP BARRIER
Index: LMDZ6/branches/cirrus/libf/dyn3dmem/conf_gcm.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3dmem/conf_gcm.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3dmem/conf_gcm.F90	(revision 5202)
@@ -22,5 +22,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
@@ -660,4 +660,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 
@@ -736,4 +745,5 @@
      write(lunout,*)' ok_dyn_ave = ', ok_dyn_ave 
      write(lunout,*)' ok_dyn_xios = ', ok_dyn_xios 
+     write(lunout,*)' adv_qsat_liq = ', adv_qsat_liq
   else
      !Config  Key  = clon
Index: LMDZ6/branches/cirrus/libf/dyn3dmem/dynetat0_loc.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3dmem/dynetat0_loc.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3dmem/dynetat0_loc.F90	(revision 5202)
@@ -42,9 +42,11 @@
   INTEGER, PARAMETER :: length=100
   INTEGER :: iq, fID, vID, idecal, ierr, iqParent, iName, iZone, iPhase, ix
-  REAL    :: time, tnat, alpha_ideal, tab_cntrl(length)    !--- RUN PARAMS TABLE
+  REAL    :: time,tab_cntrl(length)    !--- RUN PARAMS TABLE
+  REAL    :: tnat, alpha_ideal
   REAL,             ALLOCATABLE :: vcov_glo(:,:),masse_glo(:,:),   ps_glo(:)
   REAL,             ALLOCATABLE :: ucov_glo(:,:),    q_glo(:,:), phis_glo(:)
   REAL,             ALLOCATABLE :: teta_glo(:,:)
   LOGICAL :: lSkip, ll
+  LOGICAL,PARAMETER :: tnat1=.TRUE.
 !-------------------------------------------------------------------------------
   modname="dynetat0_loc"
@@ -179,6 +181,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(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat*(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal-1.)
@@ -193,5 +201,5 @@
          ! remplacant 1 par izone_init dans la ligne qui suit.
          IF(tracers(iq)%iso_iZone == 1) THEN
-          q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqIsoPha(iName,iPhase))
+           q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqIsoPha(iName,iPhase))
          ELSE
            q(ijb_u:ije_u,:,iq) =  0.
Index: LMDZ6/branches/cirrus/libf/dyn3dmem/dynredem_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3dmem/dynredem_mod.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3dmem/dynredem_mod.F90	(revision 5202)
@@ -7,5 +7,5 @@
   PRIVATE
   PUBLIC :: dynredem_write_u, dynredem_write_v, dynredem_read_u, err
-  PUBLIC :: cre_var, get_var1, put_var, fil, modname, msg
+  PUBLIC :: cre_var, put_var, fil, modname, msg
   CHARACTER(LEN=256), SAVE :: fil, modname
   INTEGER,            SAVE :: nvarid
Index: LMDZ6/branches/cirrus/libf/dyn3dmem/filtreg_p.F
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3dmem/filtreg_p.F	(revision 4945)
+++ 	(revision )
@@ -1,432 +1,0 @@
-
-
-      SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv, 
-     &     ifiltre, iaire, griscal ,iter)
-      USE parallel_lmdz, only : OMP_CHUNK
-      USE mod_filtre_fft
-      USE timer_filtre
-      
-      USE filtreg_mod
-      
-      IMPLICIT NONE
-      
-c=======================================================================
-c
-c   Auteur: P. Le Van        07/10/97
-c   ------
-c
-c   Objet: filtre matriciel longitudinal ,avec les matrices precalculees
-c                     pour l'operateur  Filtre    .
-c   ------
-c
-c   Arguments:
-c   ----------
-c
-c      
-c      ibeg..iend            lattitude a filtrer
-c      nlat                  nombre de latitudes du champ
-c      nbniv                 nombre de niveaux verticaux a filtrer
-c      champ(iip1,nblat,nbniv)  en entree : champ a filtrer
-c                            en sortie : champ filtre
-c      ifiltre               +1  Transformee directe
-c                            -1  Transformee inverse
-c                            +2  Filtre directe
-c                            -2  Filtre inverse
-c
-c      iaire                 1   si champ intensif
-c                            2   si champ extensif (pondere par les aires)
-c
-c      iter                  1   filtre simple
-c
-c=======================================================================
-c
-c
-c                      Variable Intensive
-c                ifiltre = 1     filtre directe
-c                ifiltre =-1     filtre inverse
-c
-c                      Variable Extensive
-c                ifiltre = 2     filtre directe
-c                ifiltre =-2     filtre inverse
-c
-c
-      INCLUDE "dimensions.h"
-      INCLUDE "paramet.h"
-      INCLUDE "coefils.h"
-c
-      INTEGER ibeg,iend,nlat,nbniv,ifiltre,iter
-      INTEGER i,j,l,k
-      INTEGER iim2,immjm
-      INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
-      
-      REAL  champ( iip1,nlat,nbniv)
-      
-      LOGICAL    griscal
-      INTEGER    hemisph, iaire
-      
-      REAL :: champ_fft(iip1,nlat,nbniv)
-      REAL :: champ_in(iip1,nlat,nbniv)
-      
-      LOGICAL,SAVE     :: first=.TRUE.
-c$OMP THREADPRIVATE(first) 
-
-      REAL, DIMENSION(iip1,nlat,nbniv) :: champ_loc
-      INTEGER :: ll_nb, nbniv_loc
-      REAL, SAVE :: sdd12(iim,4)
-c$OMP THREADPRIVATE(sdd12) 
-
-      INTEGER, PARAMETER :: type_sddu=1
-      INTEGER, PARAMETER :: type_sddv=2
-      INTEGER, PARAMETER :: type_unsddu=3
-      INTEGER, PARAMETER :: type_unsddv=4
-
-      INTEGER :: sdd1_type, sdd2_type
-
-      IF (first) THEN
-         sdd12(1:iim,type_sddu) = sddu(1:iim)
-         sdd12(1:iim,type_sddv) = sddv(1:iim)
-         sdd12(1:iim,type_unsddu) = unsddu(1:iim)
-         sdd12(1:iim,type_unsddv) = unsddv(1:iim)
-
-         CALL Init_timer
-         first=.FALSE.
-      ENDIF
-
-c$OMP MASTER      
-      CALL start_timer
-c$OMP END MASTER
-
-c-------------------------------------------------------c
-
-      IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
-     & CALL abort_gcm("fitreg_p","Pas de transformee simple 
-     &dans cette version",1)
-      
-      IF( iter.EQ. 2 )  THEN
-         PRINT *,' Pas d iteration du filtre dans cette version !'
-     &        , ' Utiliser old_filtreg et repasser !'
-         CALL abort_gcm("fitreg_p","stopped",1)
-      ENDIF
-
-      IF( ifiltre.EQ. -2 .AND..NOT.griscal )     THEN
-         PRINT *,' Cette routine ne calcule le filtre inverse que '
-     &        , ' sur la grille des scalaires !'
-         CALL abort_gcm("fitreg_p","stopped",1)
-      ENDIF
-
-      IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 )  THEN
-         PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
-     &        , ' corriger et repasser !'
-         CALL abort_gcm("fitreg_p","stopped",1)
-       ENDIF
-c
-
-      iim2   = iim * iim
-      immjm  = iim * jjm
-c
-c
-      IF( griscal )   THEN
-         IF( nlat. NE. jjp1 )  THEN
-            CALL abort_gcm("fitreg_p","nlat. NE. jjp1",1)
-         ELSE
-c     
-            IF( iaire.EQ.1 )  THEN
-               sdd1_type = type_sddv
-               sdd2_type = type_unsddv
-            ELSE
-               sdd1_type = type_unsddv
-               sdd2_type = type_sddv
-            ENDIF
-c
-            jdfil1 = 2
-            jffil1 = jfiltnu
-            jdfil2 = jfiltsu
-            jffil2 = jjm
-         ENDIF
-      ELSE
-         IF( nlat.NE.jjm )  THEN
-            CALL abort_gcm("fitreg_p","nlat. NE. jjm",1)
-         ELSE
-c
-            IF( iaire.EQ.1 )  THEN
-               sdd1_type = type_sddu
-               sdd2_type = type_unsddu
-            ELSE
-               sdd1_type = type_unsddu
-               sdd2_type = type_sddu
-            ENDIF
-c     
-            jdfil1 = 1
-            jffil1 = jfiltnv
-            jdfil2 = jfiltsv
-            jffil2 = jjm
-         ENDIF
-      ENDIF
-c      
-      DO hemisph = 1, 2
-c     
-         IF ( hemisph.EQ.1 )  THEN
-cym
-            jdfil = max(jdfil1,ibeg)
-            jffil = min(jffil1,iend)
-         ELSE
-cym
-            jdfil = max(jdfil2,ibeg)
-            jffil = min(jffil2,iend)
-         ENDIF
-
-
-cccccccccccccccccccccccccccccccccccccccccccc
-c Utilisation du filtre classique
-cccccccccccccccccccccccccccccccccccccccccccc
-
-         IF (.NOT. use_filtre_fft) THEN
-      
-c     !---------------------------------!
-c     ! Agregation des niveau verticaux !
-c     ! uniquement necessaire pour une  !
-c     ! execution OpenMP                !
-c     !---------------------------------!
-            ll_nb = 0
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-            DO l = 1, nbniv
-               ll_nb = ll_nb+1
-               DO j = jdfil,jffil
-                  DO i = 1, iim
-                     champ_loc(i,j,ll_nb) = 
-     &                    champ(i,j,l) * sdd12(i,sdd1_type)
-                  ENDDO
-               ENDDO
-            ENDDO
-c$OMP END DO NOWAIT
-
-            nbniv_loc = ll_nb
-
-            IF( hemisph.EQ.1 )      THEN
-               
-               IF( ifiltre.EQ.-2 )   THEN
-                  DO j = jdfil,jffil
-#ifdef BLAS
-                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
-     &                    matrinvn(1,1,j), iim, 
-     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
-     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
-#else
-                     champ_fft(:iim,j-jdfil+1,:)
-     &                    =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
-#endif
-                  ENDDO
-                  
-               ELSE IF ( griscal )     THEN
-                  DO j = jdfil,jffil
-#ifdef BLAS
-                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
-     &                    matriceun(1,1,j), iim, 
-     &                    champ_loc(1,j,1), iip1*nlat, 0.0,
-     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
-#else
-                     champ_fft(:iim,j-jdfil+1,:)
-     &                    =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
-#endif
-                  ENDDO
-                  
-               ELSE 
-                  DO j = jdfil,jffil
-#ifdef BLAS
-                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
-     &                    matricevn(1,1,j), iim, 
-     &                    champ_loc(1,j,1), iip1*nlat, 0.0, 
-     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
-#else
-                     champ_fft(:iim,j-jdfil+1,:)
-     &                    =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
-#endif
-                  ENDDO
-                  
-               ENDIF
-               
-            ELSE
-               
-               IF( ifiltre.EQ.-2 )   THEN
-                  DO j = jdfil,jffil
-#ifdef BLAS
-                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
-     &                    matrinvs(1,1,j-jfiltsu+1), iim, 
-     &                    champ_loc(1,j,1), iip1*nlat, 0.0, 
-     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
-#else
-                     champ_fft(:iim,j-jdfil+1,:)
-     &                    =matmul(matrinvs(:,:,j-jfiltsu+1),
-     &                            champ_loc(:iim,j,:))
-#endif
-                  ENDDO
-                  
-               ELSE IF ( griscal )     THEN
-                  
-                  DO j = jdfil,jffil
-#ifdef BLAS
-                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
-     &                    matriceus(1,1,j-jfiltsu+1), iim, 
-     &                    champ_loc(1,j,1), iip1*nlat, 0.0, 
-     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
-#else
-                     champ_fft(:iim,j-jdfil+1,:)
-     &                    =matmul(matriceus(:,:,j-jfiltsu+1),
-     &                            champ_loc(:iim,j,:))
-#endif
-                  ENDDO
-                  
-               ELSE 
-                  
-                  DO j = jdfil,jffil
-#ifdef BLAS
-                     CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, 
-     &                    matricevs(1,1,j-jfiltsv+1), iim, 
-     &                    champ_loc(1,j,1), iip1*nlat, 0.0, 
-     &                    champ_fft(1,j-jdfil+1,1), iip1*nlat)
-#else
-                     champ_fft(:iim,j-jdfil+1,:)
-     &                    =matmul(matricevs(:,:,j-jfiltsv+1),
-     &                            champ_loc(:iim,j,:))
-#endif
-                  ENDDO
-                  
-               ENDIF
-               
-            ENDIF
-!     c     
-            IF( ifiltre.EQ.2 )  THEN
-               
-c     !-------------------------------------!
-c     ! Dés-agregation des niveau verticaux !
-c     ! uniquement necessaire pour une      !
-c     ! execution OpenMP                    !
-c     !-------------------------------------!
-               ll_nb = 0
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-               DO l = 1, nbniv
-                  ll_nb = ll_nb + 1
-                  DO j = jdfil,jffil
-                     DO i = 1, iim
-                        champ( i,j,l ) = (champ_loc(i,j,ll_nb) 
-     &                       + champ_fft(i,j-jdfil+1,ll_nb))
-     &                       * sdd12(i,sdd2_type)
-                     ENDDO
-                  ENDDO
-               ENDDO
-c$OMP END DO NOWAIT
-               
-            ELSE
-               
-               ll_nb = 0
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-               DO l = 1, nbniv_loc
-                  ll_nb = ll_nb + 1
-                  DO j = jdfil,jffil
-                     DO i = 1, iim
-                        champ( i,j,l ) = (champ_loc(i,j,ll_nb) 
-     &                       - champ_fft(i,j-jdfil+1,ll_nb))
-     &                       * sdd12(i,sdd2_type)
-                     ENDDO
-                  ENDDO
-               ENDDO
-c$OMP END DO NOWAIT
-               
-            ENDIF
-            
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-            DO l = 1, nbniv
-               DO j = jdfil,jffil
-                  champ( iip1,j,l ) = champ( 1,j,l )
-               ENDDO
-            ENDDO
-c$OMP END DO NOWAIT
-            
-ccccccccccccccccccccccccccccccccccccccccccccc
-c Utilisation du filtre FFT
-ccccccccccccccccccccccccccccccccccccccccccccc
-        
-         ELSE
-       
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-            DO l=1,nbniv
-               DO j=jdfil,jffil
-                  DO  i = 1, iim
-                     champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
-                     champ_fft( i,j,l) = champ(i,j,l)
-                  ENDDO
-               ENDDO
-            ENDDO
-c$OMP END DO NOWAIT
-
-            IF (jdfil<=jffil) THEN
-               IF( ifiltre. EQ. -2 )   THEN
-                  CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv) 
-               ELSE IF ( griscal )     THEN
-                  CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
-               ELSE
-                  CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
-               ENDIF
-            ENDIF
-
-
-            IF( ifiltre.EQ. 2 )  THEN
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)          
-               DO l=1,nbniv
-                  DO j=jdfil,jffil
-                     DO  i = 1, iim
-                        champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
-     &                       *sdd12(i,sdd2_type)
-                     ENDDO
-                  ENDDO
-               ENDDO
-c$OMP END DO NOWAIT	  
-            ELSE
-        
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)  
-               DO l=1,nbniv
-                  DO j=jdfil,jffil
-                     DO  i = 1, iim
-                        champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
-     &                       *sdd12(i,sdd2_type)
-                     ENDDO
-                  ENDDO
-               ENDDO
-c$OMP END DO NOWAIT          
-            ENDIF
-c
-c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
-            DO l=1,nbniv
-               DO j=jdfil,jffil
-!            champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
-                  champ( iip1,j,l ) = champ( 1,j,l )
-               ENDDO
-            ENDDO
-c$OMP END DO NOWAIT          	
-         ENDIF 
-c Fin de la zone de filtrage
-
-	
-      ENDDO
-
-!      DO j=1,nlat
-!     
-!          PRINT *,"check FFT ----> Delta(",j,")=",
-!     &            sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
-!     &            sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:)) 
-!      ENDDO
-      
-!          PRINT *,"check FFT ----> Delta(",j,")=",
-!     &            sum(champ-champ_fft)/sum(champ)
-!      
-      
-c
- 1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a 
-     &     filtrer, sur la grille des scalaires'/)
- 2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
-     &     ltrer, sur la grille de V ou de Z'/)
-c$OMP MASTER      
-      CALL stop_timer
-c$OMP END MASTER
-      RETURN
-      END
Index: LMDZ6/branches/cirrus/libf/dyn3dmem/gcm.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3dmem/gcm.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3dmem/gcm.F90	(revision 5202)
@@ -480,5 +480,5 @@
   !$OMP COPYIN(saison,ecripar,fxyhypb,ysinus,read_start,ok_guide) &
   !$OMP COPYIN(ok_strato,ok_gradsfile,ok_limit,ok_etat0) &
-  !$OMP COPYIN(iflag_phys,iflag_trac)
+  !$OMP COPYIN(iflag_phys,iflag_trac,adv_qsat_liq)
   CALL leapfrog_loc(ucov,vcov,teta,ps,masse,phis,q,time_0)
   !$OMP END PARALLEL
Index: LMDZ6/branches/cirrus/libf/dyn3dmem/iniacademic_loc.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3dmem/iniacademic_loc.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3dmem/iniacademic_loc.F90	(revision 5202)
@@ -85,4 +85,5 @@
 
   REAL zdtvr, tnat, alpha_ideal
+  LOGICAL,PARAMETER :: tnat1=.true.
   
   character(len=*),parameter :: modname="iniacademic"
@@ -323,10 +324,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(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqParent)*tnat*(q(ijb_u:ije_u,:,iqParent)/30.e-3)**(alpha_ideal-1.)
-              ELSE
-                 q(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqIsoPha(iName,iPhase))
-              END IF
+              ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
+                 IF(tracers(iq)%iso_iZone == 1) THEN ! a verifier.
+                    ! 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(ijb_u:ije_u,:,iq) = q(ijb_u:ije_u,:,iqIsoPha(iName,iPhase))
+                 else !IF(tracers(iq)%iso_iZone == 1) THEN
+                    q(ijb_u:ije_u,:,iq) = 0.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/dyn3dmem/logic_mod.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3dmem/logic_mod.F90	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3dmem/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
@@ -37,5 +39,5 @@
 !$OMP     apdiss,apdelq,saison,ecripar,fxyhypb,ysinus, &
 !$OMP     read_start,ok_guide,ok_strato,ok_gradsfile, &
-!$OMP     ok_limit,ok_etat0,hybrid)
+!$OMP     ok_limit,ok_etat0,hybrid, adv_qsat_liq)
 !$OMP THREADPRIVATE(iflag_phys,iflag_trac)
 
Index: LMDZ6/branches/cirrus/libf/dyn3dmem/qminimum_loc.F
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3dmem/qminimum_loc.F	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3dmem/qminimum_loc.F	(revision 5202)
@@ -31,4 +31,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(ijb_u:ije_u), pompe
@@ -49,9 +52,4 @@
       INTEGER ixt
       INTEGER iso_verif_noNaN_nostop
-c
-c Quand l'eau liquide est trop petite (ou negative), on prend
-c l'eau vapeur de la meme couche et la convertit en eau liquide
-c (sans changer la temperature !)
-c
 
 c$OMP BARRIER
@@ -63,4 +61,10 @@
          first = .FALSE.
       END IF
+c
+c Quand l'eau liquide est trop petite (ou negative), on prend
+c l'eau vapeur de la meme couche et la convertit en eau liquide
+c (sans changer la temperature !)
+c
+
       call check_isotopes(q,ij_begin,ij_end,'qminimum 52')   
 
@@ -73,6 +77,6 @@
           zx_defau_diag(i,k,1)=0.0
           zx_defau_diag(i,k,2)=0.0
-          q_follow(i,k,1)=q(i,k,1)
-          q_follow(i,k,2)=q(i,k,2)
+          q_follow(i,k,1)=q(i,k,iq_vap)
+          q_follow(i,k,2)=q(i,k,iq_liq)
         ENDDO
 c$OMP END DO NOWAIT
@@ -80,18 +84,18 @@
 
       !write(lunout,*) 'qminimum 57'
-      DO 1000 k = 1, llm
+      DO k = 1, llm
 c$OMP DO SCHEDULE(STATIC)       
-      DO 1040 i = ijb, ije
-            if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then
-
-              if (niso > 0) zx_defau_diag(i,k,iq_liq)=AMAX1
+        DO i = ijb, ije
+          if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then
+
+            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
-c$OMP END DO NOWAIT
- 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
+        END DO
+c$OMP END DO NOWAIT
+      END DO
 
 c
@@ -100,22 +104,20 @@
 c
       !write(lunout,*) 'qminimum 81'
-      iq = iq_vap
-c
       DO k = llm, 2, -1
 ccc      zx_abc = dpres(k) / dpres(k-1)
 c$OMP DO SCHEDULE(STATIC)
-      DO i = ijb, ije
-
-         if ( seuil_vap - q(i,k,iq) .gt. 0.d0 ) then
-
-            if (niso > 0)
-     &        zx_defau_diag(i,k,iq)=AMAX1( seuil_vap - q(i,k,iq), 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  
-
-         endif
-      ENDDO
+        DO i = ijb, ije
+
+          if ( seuil_vap - q(i,k,iq_vap) .gt. 0.d0 ) then
+
+            if (niso > 0) zx_defau_diag(i,k,1)
+     &           = AMAX1( seuil_vap - q(i,k,iq_vap), 0.0 )
+
+            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
 c$OMP END DO NOWAIT
       ENDDO
@@ -129,6 +131,6 @@
 c$OMP DO SCHEDULE(STATIC)
       DO i = ijb, ije
-         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 )
          IF (zx_pump(i) > 0.0) THEN
             nb_pump = nb_pump+1
@@ -165,5 +167,5 @@
       DO i = ijb, ije
         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 = ijb, ije  
@@ -175,13 +177,13 @@
 c$OMP DO SCHEDULE(STATIC)      
         DO i = ijb, ije
-          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     
-!              write(lunout,*) 'i,k,q_follow(i,k-1,iq_vap)=',
-!     :                 i,k,q_follow(i,k-1,iq_vap)          
-              if (q_follow(i,k-1,iq_vap).lt.min_qParent) then
+!              write(lunout,*) 'i,k,q_follow(i,k-1,ivap)=',
+!     :                 i,k,q_follow(i,k-1,1)          
+              if (q_follow(i,k-1,1).lt.min_qParent) then
                 write(lunout,*) 'tmp qmin: on stoppe'
                 write(lunout,*) 'zx_pump(i)=',zx_pump(i)
-                write(lunout,*) 'q_follow(i,:,iq_vap)=',
-     :                   q_follow(i,:,iq_vap)
+                write(lunout,*) 'q_follow(i,:,ivap)=',
+     :                   q_follow(i,:,1)
                 write(lunout,*) 'k=',k
                 call abort_gcm("qminimum","not enough vapor",1)
@@ -189,14 +191,14 @@
             do ixt=1,ntiso
 !                write(lunout,*) 'qmin 168: ixt=',ixt
-!                write(lunout,*) 'q(i,k,iqIsoPha(ixt,iq_vap)=',
+!                write(lunout,*) 'q(i,k,iqIsoPha(ixt,iq_vap))=',
 !     :             q(i,k,iqIsoPha(ixt,iq_vap))
-!                write(lunout,*) 'zx_defau_diag(i,k,iq_vap)=',
-!     :                  zx_defau_diag(i,k,iq_vap)
-!                write(lunout,*) 'q(i,k-1,iqIsoPha(ixt,iq_vap)=',
+!                write(lunout,*) 'zx_defau_diag(i,k,ivap)=',
+!     :                  zx_defau_diag(i,k,1)
+!                write(lunout,*) 'q(i,k-1,iqIsoPha(ixt,iq_vap))=',
 !     :                   q(i,k-1,iqIsoPha(ixt,iq_vap))     
 
                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)
                 
               if (isoCheck) then
@@ -204,10 +206,10 @@
      :                   'qminimum 155').eq.1) then
                    write(*,*) 'i,k,ixt=',i,k,ixt
-                   write(*,*) 'q_follow(i,k-1,iq_vap)=',
-     :                   q_follow(i,k-1,iq_vap)
+                   write(*,*) 'q_follow(i,k-1,ivap)=',
+     :                   q_follow(i,k-1,1)
                    write(*,*) 'q(i,k,iqIsoPha(ixt,iq_vap))=',
      :                   q(i,k,iqIsoPha(ixt,iq_vap))
-                   write(*,*) 'zx_defau_diag(i,k,iq_vap)=',
-     :                   zx_defau_diag(i,k,iq_vap)
+                   write(*,*) 'zx_defau_diag(i,k,ivap)=',
+     :                   zx_defau_diag(i,k,1)
                    write(*,*) 'q(i,k-1,iqIsoPha(ixt,iq_vap))=',
      :                   q(i,k-1,iqIsoPha(ixt,iq_vap))
@@ -219,8 +221,8 @@
                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)
 
                if (isoCheck) then
@@ -229,10 +231,10 @@
      :                   'qminimum 175').eq.1) then
                    write(*,*) 'k,i,ixt=',k,i,ixt
-                   write(*,*) 'q_follow(i,k-1,iq_vap)=',
-     :                   q_follow(i,k-1,iq_vap)
+                   write(*,*) 'q_follow(i,k-1,ivap)=',
+     :                   q_follow(i,k-1,1)
                    write(*,*) 'q(i,k,iqIsoPha(ixt,iq_vap))=',
      :                   q(i,k,iqIsoPha(ixt,iq_vap))
-                   write(*,*) 'zx_defau_diag(i,k,iq_vap)=',
-     :                   zx_defau_diag(i,k,iq_vap)
+                   write(*,*) 'zx_defau_diag(i,k,ivap)=',
+     :                   zx_defau_diag(i,k,1)
                    write(*,*) 'q(i,k-1,iqIsoPha(ixt,iq_vap))=',
      :                   q(i,k-1,iqIsoPha(ixt,iq_vap))
@@ -242,10 +244,10 @@
 
               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        
 c$OMP END DO NOWAIT
@@ -260,21 +262,21 @@
 c$OMP DO SCHEDULE(STATIC)
         DO i = ijb, ije
-          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 = ijb, ije
 c$OMP END DO NOWAIT       
Index: LMDZ6/branches/cirrus/libf/dyn3dmem/vlspltgen_loc.F
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3dmem/vlspltgen_loc.F	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3dmem/vlspltgen_loc.F	(revision 5202)
@@ -10,5 +10,5 @@
 c
 c    ********************************************************************
-c          Shema  d'advection " pseudo amont " .
+c          Schema  d'advection " pseudo amont " .
 c      + test sur humidite specifique: Q advecte< Qsat aval
 c                   (F. Codron, 10/99)
@@ -32,4 +32,5 @@
       USE vlspltgen_mod
       USE comconst_mod, ONLY: cpp
+      USE logic_mod, ONLY: adv_qsat_liq
       IMPLICIT NONE
 
@@ -108,5 +109,9 @@
          ENDDO
          DO ij = ijb, ije
-          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 )
Index: LMDZ6/branches/cirrus/libf/dyn3dmem/vlspltqs_loc.F
===================================================================
--- LMDZ6/branches/cirrus/libf/dyn3dmem/vlspltqs_loc.F	(revision 4945)
+++ LMDZ6/branches/cirrus/libf/dyn3dmem/vlspltqs_loc.F	(revision 5202)
@@ -806,5 +806,5 @@
          IF (pole_sud) THEN
          
-           convps  = -SSUM(iim,qbyv(ip1jm-iim,l,iq),iq,1)/apols
+           convps  = -SSUM(iim,qbyv(ip1jm-iim,l,iq),1)/apols
            convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
            DO ij = ip1jm+1,ip1jmp1
