Index: /LMDZ5/trunk/arch/arch-X64_ADA.fcm
===================================================================
--- /LMDZ5/trunk/arch/arch-X64_ADA.fcm	(revision 2269)
+++ /LMDZ5/trunk/arch/arch-X64_ADA.fcm	(revision 2270)
@@ -5,5 +5,5 @@
 %FPP_FLAGS           -P -traditional -I/smplocal/pub/FFTW/3.3.3_dyn/include/
 %FPP_DEF             NC_DOUBLE FFT_FFTW
-%BASE_FFLAGS         -auto -mcmodel=large -integer-size 32 -real-size 64 -align all 
+%BASE_FFLAGS         -auto -recursive -mcmodel=large -integer-size 32 -real-size 64 -align all 
 %PROD_FFLAGS         -O2 -ip -fp-model strict -axAVX,SSE4.2
 %DEV_FFLAGS          -p -g -O1 -fpe0 -traceback
Index: /LMDZ5/trunk/arch/arch-X64_CURIE.fcm
===================================================================
--- /LMDZ5/trunk/arch/arch-X64_CURIE.fcm	(revision 2269)
+++ /LMDZ5/trunk/arch/arch-X64_CURIE.fcm	(revision 2270)
@@ -5,5 +5,5 @@
 %FPP_FLAGS           -P -traditional
 %FPP_DEF             NC_DOUBLE FFT_MKL
-#%BASE_FFLAGS         -xHost -i4 -r8 -auto -align all -I$(MKL_INC_DIR) -I$(MKL_INC_DIR)/intel64/lp64
+#%BASE_FFLAGS        -recursive -xHost -i4 -r8 -auto -align all -I$(MKL_INC_DIR) -I$(MKL_INC_DIR)/intel64/lp64
 %BASE_FFLAGS         -i4 -r8 -auto -align all -I$(MKL_INC_DIR) -I$(MKL_INC_DIR)/intel64/lp64 -fp-model strict
 %PROD_FFLAGS         -O2
Index: /LMDZ5/trunk/arch/arch-gfortran.fcm
===================================================================
--- /LMDZ5/trunk/arch/arch-gfortran.fcm	(revision 2269)
+++ /LMDZ5/trunk/arch/arch-gfortran.fcm	(revision 2270)
@@ -5,5 +5,5 @@
 %FPP_FLAGS           -P -traditional
 %FPP_DEF             NC_DOUBLE 
-%BASE_FFLAGS         -c -fdefault-real-8
+%BASE_FFLAGS         -c -fdefault-real-8 -frecursive
 %PROD_FFLAGS         -O3
 %DEV_FFLAGS          -O
Index: /LMDZ5/trunk/arch/arch-gfortran_CICLAD.fcm
===================================================================
--- /LMDZ5/trunk/arch/arch-gfortran_CICLAD.fcm	(revision 2269)
+++ /LMDZ5/trunk/arch/arch-gfortran_CICLAD.fcm	(revision 2270)
@@ -5,5 +5,5 @@
 %FPP_FLAGS           -P -traditional
 %FPP_DEF             NC_DOUBLE 
-%BASE_FFLAGS         -c -fdefault-real-8 -fcray-pointer
+%BASE_FFLAGS         -c -fdefault-real-8 -fcray-pointer -frecursive
 %PROD_FFLAGS         -O3
 %DEV_FFLAGS          -O -Wall -fbounds-check
Index: /LMDZ5/trunk/libf/dyn3d/advtrac.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/advtrac.F90	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3d/advtrac.F90	(revision 2270)
@@ -9,5 +9,5 @@
   !            M.A Filiberti (04/2002)
   !
-  USE infotrac, ONLY: nqtot, iadv
+  USE infotrac, ONLY: nqtot, iadv,nqperes,ok_iso_verif
   USE control_mod, ONLY: iapp_tracvl, day_step
 
@@ -223,5 +223,10 @@
      !     Appel des sous programmes d'advection
      !-----------------------------------------------------------
-     do iq=1,nqtot
+
+     if (ok_iso_verif) then
+           call check_isotopes_seq(q,1,ip1jmp1,'advtrac 162')
+     endif !if (ok_iso_verif) then
+
+     do iq=1,nqperes
         !        call clock(t_initial)
         if(iadv(iq) == 0) cycle 
@@ -230,6 +235,9 @@
         !   ----------------------------------------------------------------
         if(iadv(iq).eq.10) THEN
-           call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
-
+           ! CRisi: on fait passer tout q pour avoir acces aux fils
+           
+           write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)     
+           call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
+           
            !   ----------------------------------------------------------------
            !   Schema "pseudo amont" + test sur humidite specifique
@@ -238,6 +246,8 @@
         else if(iadv(iq).eq.14) then
            !
-           CALL vlspltqs( q(1,1,1), 2., massem, wg , &
-                pbarug,pbarvg,dtvr,p,pk,teta )
+           write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
+           CALL vlspltqs( q, 2., massem, wg , &
+                pbarug,pbarvg,dtvr,p,pk,teta,iq)
+           
            !   ----------------------------------------------------------------
            !   Schema de Frederic Hourdin
@@ -388,4 +398,7 @@
      end DO
 
+     if (ok_iso_verif) then
+           call check_isotopes_seq(q,ip1jmp1,'advtrac 397')
+     endif !if (ok_iso_verif) then
 
      !------------------------------------------------------------------
Index: /LMDZ5/trunk/libf/dyn3d/caladvtrac.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/caladvtrac.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3d/caladvtrac.F	(revision 2270)
@@ -53,5 +53,8 @@
       if (planet_type.eq."earth") then
 C initialisation
-        dq(:,:,1:2)=q(:,:,1:2)
+! CRisi: il faut gérer tous les traceurs si on veut pouvoir faire des
+! isotopes
+!        dq(:,:,1:2)=q(:,:,1:2)
+        dq(:,:,1:nqtot)=q(:,:,1:nqtot)
        
 c  test des valeurs minmax
@@ -81,6 +84,6 @@
            ENDDO
           ENDDO
-          
-          CALL qminimum( q, 2, finmasse )
+
+          CALL qminimum( q, nqtot, finmasse )
 
           CALL SCOPY   ( ip1jmp1*llm, masse, 1, finmasse,       1 )
@@ -92,5 +95,5 @@
           dtvrtrac = iapp_tracvl * dtvr
 c
-           DO iq = 1 , 2
+           DO iq = 1 , nqtot
             DO l = 1 , llm
              DO ij = 1,ip1jmp1
@@ -105,5 +108,5 @@
         if (planet_type.eq."earth") then
 ! Earth-specific treatment for the first 2 tracers (water)
-          dq(:,:,1:2)=0.
+          dq(:,:,1:nqtot)=0.
         endif ! of if (planet_type.eq."earth")
       ENDIF ! of IF( iapptrac.EQ.iapp_tracvl )
Index: /LMDZ5/trunk/libf/dyn3d/check_isotopes.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/check_isotopes.F	(revision 2270)
+++ /LMDZ5/trunk/libf/dyn3d/check_isotopes.F	(revision 2270)
@@ -0,0 +1,195 @@
+        subroutine check_isotopes_seq(q,ip1jmp1,err_msg)
+        USE infotrac
+        implicit none
+
+#include "dimensions.h"
+
+        ! inputs
+        integer ip1jmp1
+        real q(ip1jmp1,llm,nqtot)
+        character*(*) err_msg ! message d''erreur à afficher
+
+        ! locals
+        integer ixt,phase,k,i,iq,iiso,izone,ieau,iqeau
+        real xtractot,xiiso
+        real borne
+        real qmin
+        real errmax ! erreur maximale en absolu.
+        real errmaxrel ! erreur maximale en relatif autorisée
+        real deltaDmax,deltaDmin
+        real ridicule
+        parameter (borne=1e19)
+        parameter (errmax=1e-8)
+        parameter (errmaxrel=1e-3)
+        parameter (qmin=1e-11)
+        parameter (deltaDmax=200.0,deltaDmin=-999.9)
+        parameter (ridicule=1e-12)
+        real deltaD
+
+        if (ok_isotopes) then
+
+        ! verifier que rien n'est NaN
+        do ixt=1,ntraciso
+          do phase=1,nqo
+            iq=iqiso(ixt,phase)
+            do k=1,llm
+              DO i = 1,ip1jmp1
+                if ((q(i,k,iq).gt.-borne).and.
+     :            (q(i,k,iq).lt.borne)) then
+                else !if ((x(ixt,i,j).gt.-borne).and.
+                  write(*,*) 'erreur detectee par iso_verif_noNaN:'
+                  write(*,*) err_msg
+                  write(*,*) 'q,i,k,iq=',q(i,k,iq),i,k,iq
+                  write(*,*) 'borne=',borne
+                  stop
+                endif  !if ((x(ixt,i,j).gt.-borne).and.
+              enddo !DO i = 1,ip1jmp1
+            enddo !do k=1,llm
+          enddo !do phase=1,nqo
+        enddo !do ixt=1,ntraciso
+
+        ! verifier que l'eau normale est OK
+        if (use_iso(1)) then
+          ixt=indnum_fn_num(1)
+          do phase=1,nqo
+            iq=iqiso(ixt,phase)
+            do k=1,llm
+            DO i = 1,ip1jmp1  
+              if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
+     :          (abs((q(i,k,phase)-q(i,k,iq))/
+     :           max(max(abs(q(i,k,phase)),abs(q(i,k,iq))),1e-18))
+     :           .gt.errmaxrel)) then
+                  write(*,*) 'erreur detectee par iso_verif_egalite:'
+                  write(*,*) err_msg
+                  write(*,*) 'ixt,phase=',ixt,phase
+                  write(*,*) 'q,iq,i,k=',q(i,k,iq),iq,i,k
+                  write(*,*) 'q(i,k,phase)=',q(i,k,phase)
+                  stop
+              endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
+              ! bidouille pour éviter divergence:
+              q(i,k,iq)= q(i,k,phase) 
+            enddo ! DO i = 1,ip1jmp1
+            enddo !do k=1,llm
+          enddo ! do phase=1,nqo 
+        endif !if (use_iso(1)) then
+        
+        ! verifier que HDO est raisonable
+        if (use_iso(2)) then
+          ixt=indnum_fn_num(2)
+          do phase=1,nqo
+            iq=iqiso(ixt,phase)
+            do k=1,llm
+            DO i = 1,ip1jmp1
+            if (q(i,k,iq).gt.qmin) then
+             deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(2)-1)*1000
+             if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+                  write(*,*) 'erreur detectee par iso_verif_aberrant:'
+                  write(*,*) err_msg
+                  write(*,*) 'ixt,phase=',ixt,phase
+                  write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k
+                  write(*,*) 'q=',q(i,k,:)
+                  write(*,*) 'deltaD=',deltaD
+                  stop
+             endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+            endif !if (q(i,k,iq).gt.qmin) then
+            enddo !DO i = 1,ip1jmp1
+            enddo !do k=1,llm
+          enddo ! do phase=1,nqo 
+        endif !if (use_iso(2)) then
+
+        ! verifier que O18 est raisonable
+        if (use_iso(3)) then
+          ixt=indnum_fn_num(3)
+          do phase=1,nqo
+            iq=iqiso(ixt,phase)
+            do k=1,llm
+            DO i = 1,ip1jmp1
+            if (q(i,k,iq).gt.qmin) then
+             deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(3)-1)*1000
+             if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+                  write(*,*) 'erreur detectee iso_verif_aberrant O18:'
+                  write(*,*) err_msg
+                  write(*,*) 'ixt,phase=',ixt,phase
+                  write(*,*) 'q,iq,i,k,=',q(i,k,phase),iq,i,k
+                  write(*,*) 'xt=',q(i,k,:)
+                  write(*,*) 'deltaO18=',deltaD
+                  stop
+             endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+            endif !if (q(i,k,iq).gt.qmin) then
+            enddo !DO i = 1,ip1jmp1
+            enddo !do k=1,llm
+          enddo ! do phase=1,nqo 
+        endif !if (use_iso(2)) then
+
+
+        if (ok_isotrac) then
+
+          if (use_iso(2).and.use_iso(1)) then
+            do izone=1,ntraceurs_zone
+             ixt=index_trac(izone,indnum_fn_num(2))
+             ieau=index_trac(izone,indnum_fn_num(1))
+             do phase=1,nqo
+               iq=iqiso(ixt,phase)
+               iqeau=iqiso(ieau,phase)
+               do k=1,llm
+                DO i = 1,ip1jmp1
+                if (q(i,k,iq).gt.qmin) then
+                 deltaD=(q(i,k,iq)/q(i,k,iqeau)/tnat(2)-1)*1000
+                 if ((deltaD.gt.deltaDmax).or.
+     &                   (deltaD.lt.deltaDmin)) then
+                  write(*,*) 'erreur dans iso_verif_aberrant trac:'
+                  write(*,*) err_msg
+                  write(*,*) 'izone,phase=',izone,phase
+                  write(*,*) 'ixt,ieau=',ixt,ieau
+                  write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k
+                  write(*,*) 'deltaD=',deltaD
+                  stop
+                 endif !if ((deltaD.gt.deltaDmax).or.
+                endif !if (q(i,k,iq).gt.qmin) then
+                enddo !DO i = 1,ip1jmp1
+                enddo  ! do k=1,llm
+              enddo ! do phase=1,nqo    
+            enddo !do izone=1,ntraceurs_zone
+          endif !if (use_iso(2).and.use_iso(1)) then
+
+          do iiso=1,niso
+           do phase=1,nqo
+              iq=iqiso(iiso,phase)
+              do k=1,llm
+                DO i = 1,ip1jmp1
+                   xtractot=0.0
+                   xiiso=q(i,k,iq)
+                   do izone=1,ntraceurs_zone
+                      iq=iqiso(index_trac(izone,iiso),phase)
+                      xtractot=xtractot+ q(i,k,iq)
+                   enddo !do izone=1,ntraceurs_zone
+                   if ((abs(xtractot-xiiso).gt.errmax).and.
+     :                  (abs(xtractot-xiiso)/
+     :                  max(max(abs(xtractot),abs(xiiso)),1e-18)
+     :                  .gt.errmaxrel)) then
+                  write(*,*) 'erreur detectee par iso_verif_traceurs:'
+                  write(*,*) err_msg
+                  write(*,*) 'iiso,phase=',iiso,phase
+                  write(*,*) 'i,k,=',i,k
+                  write(*,*) 'q(i,k,:)=',q(i,k,:)
+                  stop
+                 endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
+                  
+                 ! bidouille pour éviter divergence:
+                 if (abs(xtractot).gt.ridicule) then
+                   do izone=1,ntraceurs_zone
+                     ixt=index_trac(izone,iiso) 
+                     q(i,k,iq)=q(i,k,iq)/xtractot*xiiso
+                   enddo !do izone=1,ntraceurs_zone                
+                  endif !if ((abs(xtractot).gt.ridicule) then
+                enddo !DO i = 1,ip1jmp1
+              enddo !do k=1,llm
+           enddo !do phase=1,nqo
+          enddo !do iiso=1,niso
+
+        endif !if (ok_isotrac) then
+
+        endif ! if (ok_isotopes)
+        
+        end
+
Index: /LMDZ5/trunk/libf/dyn3d/dynetat0.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/dynetat0.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3d/dynetat0.F	(revision 2270)
@@ -297,4 +297,19 @@
            write(lunout,*)"          Il est donc initialise a zero"
            q(:,:,:,iq)=0.
+
+           ! CRisi: pour les isotopes, on peut faire init théorique
+           ! distill de Rayleigh très simplifiée
+           if (ok_isotopes) then
+              if ((iso_num(iq).gt.0).and.(zone_num(iq).eq.0)) then
+                q(:,:,:,iq)=q(:,:,:,iqpere(iq))                         
+     &                   *tnat(iso_num(iq))                             
+     &                   *(q(:,:,:,iqpere(iq))/30.e-3)                  
+     &                   **(alpha_ideal(iso_num(iq))-1)
+              endif
+              if ((iso_num(iq).gt.0).and.(zone_num(iq).eq.1)) then
+                  q(:,:,:,iq)=q(:,:,:,iqiso(iso_indnum(iq),
+     &                   phase_num(iq)))
+              endif  
+           endif !if (ok_isotopes) then 
         ELSE
            ierr = NF90_GET_VAR(nid, nvarid, q(:,:,:,iq))
Index: /LMDZ5/trunk/libf/dyn3d/iniacademic.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/iniacademic.F90	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3d/iniacademic.F90	(revision 2270)
@@ -5,5 +5,5 @@
 
   USE filtreg_mod, ONLY: inifilr
-  USE infotrac, ONLY : nqtot
+  USE infotrac
   USE control_mod, ONLY: day_step,planet_type
 #ifdef CPP_IOIPSL
@@ -262,8 +262,27 @@
               if (i == 2) q(:,:,i)=1.e-15
               if (i.gt.2) q(:,:,i)=0.
+
+              ! CRisi: init des isotopes
+              ! distill de Rayleigh très simplifiée
+              if (ok_isotopes) then
+                if ((iso_num(i).gt.0).and.(zone_num(i).eq.0)) then          
+                   q(:,:,i)=q(:,:,iqpere(i))             &
+      &                  *tnat(iso_num(i))               &
+      &                  *(q(:,:,iqpere(i))/30.e-3)      &
+      &                  **(alpha_ideal(iso_num(i))-1)
+                endif                
+                if ((iso_num(i).gt.0).and.(zone_num(i).eq.1)) then
+                  q(:,:,i)=q(:,:,iqiso(iso_indnum(i),phase_num(i)))
+                endif
+              endif !if (ok_isotopes) then
+
            enddo
         else
            q(:,:,:)=0
         endif ! of if (planet_type=="earth")
+
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
+        endif !if (ok_iso_verif) then
 
         ! add random perturbation to temperature
Index: /LMDZ5/trunk/libf/dyn3d/leapfrog.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/leapfrog.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3d/leapfrog.F	(revision 2270)
@@ -11,5 +11,5 @@
       use IOIPSL
 #endif
-      USE infotrac, ONLY: nqtot
+      USE infotrac, ONLY: nqtot,ok_iso_verif
       USE guide_mod, ONLY : guide_main
       USE write_field, ONLY: writefield
@@ -235,4 +235,7 @@
       jH_cur = jH_cur - int(jH_cur)
 
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,ip1jmp1,'leapfrog 321')
+        endif !if (ok_iso_verif) then
 
 #ifdef CPP_IOIPSL
@@ -265,4 +268,8 @@
 !      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
 !      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
+
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,ip1jmp1,'leapfrog 400')
+        endif !if (ok_iso_verif) then
 
    2  CONTINUE ! Matsuno backward or leapfrog step begins here
@@ -305,4 +312,9 @@
       endif
 
+
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,ip1jmp1,'leapfrog 589')
+        endif !if (ok_iso_verif) then
+
 c-----------------------------------------------------------------------
 c   calcul des tendances dynamiques:
@@ -321,4 +333,9 @@
 c   calcul des tendances advection des traceurs (dont l'humidite)
 c   -------------------------------------------------------------
+
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,ip1jmp1,
+     &           'leapfrog 686: avant caladvtrac')
+        endif !if (ok_iso_verif) then
 
       IF( forward. OR . leapf )  THEN
@@ -327,4 +344,5 @@
      *        p, masse, dq,  teta,
      .        flxw, pk)
+
          
          IF (offline) THEN
@@ -346,9 +364,17 @@
 c   ----------------------------------
 
-
-       CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
+        if (ok_iso_verif) then
+           write(*,*) 'leapfrog 720' 
+           call check_isotopes_seq(q,ip1jmp1,'leapfrog 756')
+        endif !if (ok_iso_verif) then
+        
+       CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 ,
      $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
 !     $              finvmaold                                    )
 
+       if (ok_iso_verif) then
+          write(*,*) 'leapfrog 724'
+           call check_isotopes_seq(q,ip1jmp1,'leapfrog 762')
+        endif !if (ok_iso_verif) then
 
 c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
@@ -516,4 +542,7 @@
         CALL massdair(p,masse)
 
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
+        endif !if (ok_iso_verif) then
 
 c-----------------------------------------------------------------------
@@ -600,4 +629,8 @@
 c   preparation du pas d'integration suivant  ......
 
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
+        endif !if (ok_iso_verif) then
+
       IF ( .NOT.purmats ) THEN
 c       ........................................................
@@ -657,4 +690,8 @@
             ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
 
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
+        endif !if (ok_iso_verif) then
+
 c-----------------------------------------------------------------------
 c   ecriture de la bande histoire:
@@ -735,4 +772,8 @@
       ELSE ! of IF (.not.purmats)
 
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
+        endif !if (ok_iso_verif) then
+
 c       ........................................................
 c       ..............       schema  matsuno        ...............
@@ -757,4 +798,8 @@
 
             ELSE ! of IF(forward) i.e. backward step
+ 
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
+        endif !if (ok_iso_verif) then  
 
               IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
Index: /LMDZ5/trunk/libf/dyn3d/qminimum.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/qminimum.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3d/qminimum.F	(revision 2270)
@@ -2,6 +2,7 @@
 ! $Header$
 !
-      SUBROUTINE qminimum( q,nq,deltap )
+      SUBROUTINE qminimum( q,nqtot,deltap )
 
+      USE infotrac, ONLY: ok_isotopes,ntraciso,iqiso,ok_iso_verif
       IMPLICIT none
 c
@@ -13,6 +14,6 @@
 #include "comvert.h"
 c
-      INTEGER nq
-      REAL q(ip1jmp1,llm,nq), deltap(ip1jmp1,llm)
+      INTEGER nqtot
+      REAL q(ip1jmp1,llm,nqtot), deltap(ip1jmp1,llm)
 c
       INTEGER iq_vap, iq_liq
@@ -30,4 +31,7 @@
       INTEGER i, k, iq
       REAL zx_defau, zx_abc, zx_pump(ip1jmp1), pompe
+
+      real zx_defau_diag(ip1jmp1,llm,2) 
+      real q_follow(ip1jmp1,llm,2)
 c
       REAL SSUM
@@ -36,4 +40,8 @@
       SAVE imprim
       DATA imprim /0/
+      !INTEGER ijb,ije
+      !INTEGER Index_pump(ij_end-ij_begin+1)
+      !INTEGER nb_pump
+      INTEGER ixt
 c
 c Quand l'eau liquide est trop petite (ou negative), on prend
@@ -41,7 +49,20 @@
 c (sans changer la temperature !)
 c
+
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,ip1jmp1,'qminimum 52')   
+        endif !if (ok_iso_verif) then     
+
+      zx_defau_diag(:,:,:)=0.0
+      q_follow(:,:,:)=q(:,:,:)  
       DO 1000 k = 1, llm
         DO 1040 i = 1, ip1jmp1
           if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then
+
+              if (ok_isotopes) then
+                 zx_defau_diag(i,k,iq_liq)=AMAX1
+     :               ( seuil_liq - q(i,k,iq_liq), 0.0 )
+              endif !if (ok_isotopes) then
+
              q(i,k,iq_vap) = q(i,k,iq_vap) + q(i,k,iq_liq) - seuil_liq
              q(i,k,iq_liq) = seuil_liq
@@ -59,4 +80,9 @@
         DO i = 1, ip1jmp1
           if ( seuil_vap - q(i,k,iq) .gt. 0.d0 ) then
+
+            if (ok_isotopes) then
+              zx_defau_diag(i,k,iq)=AMAX1( seuil_vap - q(i,k,iq), 0.0 )
+            endif !if (ok_isotopes) then
+
             q(i,k-1,iq) =  q(i,k-1,iq) - ( seuil_vap - q(i,k,iq) ) *
      &                     deltap(i,k) / deltap(i,k-1)
@@ -83,4 +109,83 @@
          ENDDO
       ENDIF
+
+      write(*,*) 'qminimum 128'
+      if (ok_isotopes) then
+      ! CRisi: traiter de même les traceurs d'eau
+      ! Mais il faut les prendre à l'envers pour essayer de conserver la
+      ! masse.
+      ! 1) pompage dans le sol  
+      ! On suppose que ce pompage se fait sans isotopes -> on ne modifie
+      ! rien ici et on croise les doigts pour que ça ne soit pas trop
+      ! génant
+      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)
+        endif !if (zx_pump(i).gt.0.0) then
+      enddo !DO i = 1,ip1jmp1
+
+      ! 2) transfert de vap vers les couches plus hautes
+      !write(*,*) 'qminimum 139'
+      do k=2,llm
+        DO i = 1,ip1jmp1
+          if (zx_defau_diag(i,k,iq_vap).gt.0.0) then             
+              ! on ajoute la vapeur en k              
+              do ixt=1,ntraciso
+               q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap))
+     :              +zx_defau_diag(i,k,iq_vap)
+     :              *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap)
+                
+              ! et on la retranche en k-1
+               q(i,k-1,iqiso(ixt,iq_vap))=q(i,k-1,iqiso(ixt,iq_vap))
+     :              -zx_defau_diag(i,k,iq_vap)
+     :              *deltap(i,k)/deltap(i,k-1)
+     :              *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap)
+
+              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)
+     :              *deltap(i,k)/deltap(i,k-1)
+          endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then
+        enddo !DO i = 1, ip1jmp1        
+       enddo !do k=2,llm
+
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,ip1jmp1,'qminimum 168')
+        endif !if (ok_iso_verif) then
+        
+      
+        ! 3) transfert d'eau de la vapeur au liquide
+        !write(*,*) 'qminimum 164'
+        do k=1,llm
+        DO i = 1,ip1jmp1
+          if (zx_defau_diag(i,k,iq_liq).gt.0.0) then
+
+              ! on ajoute eau liquide en k en k              
+              do ixt=1,ntraciso
+               q(i,k,iqiso(ixt,iq_liq))=q(i,k,iqiso(ixt,iq_liq))
+     :              +zx_defau_diag(i,k,iq_liq)
+     :              *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap)
+              ! et on la retranche à la vapeur en k
+               q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap))
+     :              -zx_defau_diag(i,k,iq_liq)
+     :              *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap)   
+              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
+        enddo !DO i = 1, ip1jmp1
+       enddo !do k=2,llm  
+
+        if (ok_iso_verif) then
+           call check_isotopes_seq(q,ip1jmp1,'qminimum 197')
+        endif !if (ok_iso_verif) then
+
+      endif !if (ok_isotopes) then
+      !write(*,*) 'qminimum 188'
+      
 c
       RETURN
Index: /LMDZ5/trunk/libf/dyn3d/vlsplt.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/vlsplt.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3d/vlsplt.F	(revision 2270)
@@ -3,5 +3,6 @@
 c
 
-      SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt)
+      SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq)
+      USE infotrac, ONLY: nqtot,nqdesc,iqfils
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -32,7 +33,8 @@
 c      REAL masse(iip1,jjp1,llm),pente_max
       REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
-      REAL q(ip1jmp1,llm)
+      REAL q(ip1jmp1,llm,nqtot)
 c      REAL q(iip1,jjp1,llm)
       REAL w(ip1jmp1,llm),pdt
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -42,9 +44,9 @@
       INTEGER ijlqmin,iqmin,jqmin,lqmin
 c
-      REAL zm(ip1jmp1,llm),newmasse
+      REAL zm(ip1jmp1,llm,nqtot),newmasse
       REAL mu(ip1jmp1,llm)
       REAL mv(ip1jm,llm)
       REAL mw(ip1jmp1,llm+1)
-      REAL zq(ip1jmp1,llm),zz
+      REAL zq(ip1jmp1,llm,nqtot),zz
       REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
       REAL second,temps0,temps1,temps2,temps3
@@ -55,4 +57,5 @@
       SAVE temps1,temps2,temps3
       INTEGER iminn,imaxx
+      INTEGER ifils,iq2 ! CRisi
 
       REAL qmin,qmax
@@ -79,27 +82,35 @@
          mw(ij,llm+1)=0.
       ENDDO
-      
-      CALL SCOPY(ijp1llm,q,1,zq,1)
-      CALL SCOPY(ijp1llm,masse,1,zm,1)
+           
+      CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1)
+      CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
+        
+      if (nqdesc(iq).gt.0) then  
+        do ifils=1,nqdesc(iq)
+          iq2=iqfils(ifils,iq)
+          CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
+        enddo  
+      endif !if (nqfils(iq).gt.0) then
 
 cprint*,'Entree vlx1'
 c	call minmaxq(zq,qmin,qmax,'avant vlx     ')
-      call vlx(zq,pente_max,zm,mu)
+      call vlx(zq,pente_max,zm,mu,iq)
 cprint*,'Sortie vlx1'
 c	call minmaxq(zq,qmin,qmax,'apres vlx1    ')
 
 c print*,'Entree vly1'
-      call vly(zq,pente_max,zm,mv)
+
+      call vly(zq,pente_max,zm,mv,iq)
 c	call minmaxq(zq,qmin,qmax,'apres vly1     ')
 cprint*,'Sortie vly1'
-      call vlz(zq,pente_max,zm,mw)
+      call vlz(zq,pente_max,zm,mw,iq)
 c	call minmaxq(zq,qmin,qmax,'apres vlz     ')
 
 
-      call vly(zq,pente_max,zm,mv)
+      call vly(zq,pente_max,zm,mv,iq)
 c	call minmaxq(zq,qmin,qmax,'apres vly     ')
 
 
-      call vlx(zq,pente_max,zm,mu)
+      call vlx(zq,pente_max,zm,mu,iq)
 c	call minmaxq(zq,qmin,qmax,'apres vlx2    ')
 	
@@ -107,14 +118,29 @@
       DO l=1,llm
          DO ij=1,ip1jmp1
-           q(ij,l)=zq(ij,l)
+           q(ij,l,iq)=zq(ij,l,iq)
          ENDDO
          DO ij=1,ip1jm+1,iip1
-            q(ij+iim,l)=q(ij,l)
-         ENDDO
-      ENDDO
+            q(ij+iim,l,iq)=q(ij,l,iq)
+         ENDDO
+      ENDDO
+      ! CRisi: aussi pour les fils
+      if (nqdesc(iq).gt.0) then
+      do ifils=1,nqdesc(iq)
+        iq2=iqfils(ifils,iq)
+        DO l=1,llm
+         DO ij=1,ip1jmp1
+           q(ij,l,iq2)=zq(ij,l,iq2)
+         ENDDO
+         DO ij=1,ip1jm+1,iip1
+            q(ij+iim,l,iq2)=q(ij,l,iq2)
+         ENDDO
+        ENDDO
+      enddo !do ifils=1,nqdesc(iq)    
+      endif ! if (nqdesc(iq).gt.0) then    
 
       RETURN
       END
-      SUBROUTINE vlx(q,pente_max,masse,u_m)
+      SUBROUTINE vlx(q,pente_max,masse,u_m,iq)
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
 
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -139,8 +165,9 @@
 c   Arguments:
 c   ----------
-      REAL masse(ip1jmp1,llm),pente_max
+      REAL masse(ip1jmp1,llm,nqtot),pente_max
       REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
-      REAL q(ip1jmp1,llm)
+      REAL q(ip1jmp1,llm,nqtot)
       REAL w(ip1jmp1,llm)
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -155,4 +182,8 @@
       REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
       REAL u_mq(ip1jmp1,llm)
+
+      ! CRisi
+      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
+      INTEGER ifils,iq2 ! CRisi
 
       Logical extremum,first,testcpu
@@ -188,7 +219,7 @@
          DO l = 1, llm
             DO ij=iip2,ip1jm-1
-               dxqu(ij)=q(ij+1,l)-q(ij,l)
+               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
 c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
-c              sigu(ij)=u_m(ij,l)/masse(ij,l)
+c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
             ENDDO
             DO ij=iip1+iip1,ip1jm,iip1
@@ -243,5 +274,5 @@
          DO l = 1, llm
             DO ij=iip2,ip1jm-1
-               dxqu(ij)=q(ij+1,l)-q(ij,l)
+               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
             ENDDO
             DO ij=iip1+iip1,ip1jm,iip1
@@ -285,11 +316,11 @@
       DO l=1,llm
        DO ij=iip2,ip1jm-1
-          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
-     ,                     1.+u_m(ij,l)/masse(ij+1,l),
+          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
+     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
      ,                     u_m(ij,l))
           zdum(ij,l)=0.5*zdum(ij,l)
           u_mq(ij,l)=cvmgp(
-     ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
-     ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
+     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
+     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
      ,                u_m(ij,l))
           u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
@@ -303,11 +334,12 @@
       DO l=1,llm
        DO ij=iip2,ip1jm-1
-c	print*,'masse(',ij,')=',masse(ij,l)
+c	print*,'masse(',ij,')=',masse(ij,l,iq)
           IF (u_m(ij,l).gt.0.) THEN
-             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
-             u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
+             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
+             u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l))
           ELSE
-             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
-             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
+             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
+             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq)
+     &           -0.5*zdum(ij,l)*dxq(ij+1,l))
           ENDIF
        ENDDO
@@ -379,7 +411,8 @@
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(zu_m.gt.masse(ijq,l))
-                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
-                        zu_m=zu_m-masse(ijq,l)
+                     do while(zu_m.gt.masse(ijq,l,iq))
+                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
+     &                          *masse(ijq,l,iq)
+                        zu_m=zu_m-masse(ijq,l,iq)
                         i=mod(i-2+iim,iim)+1
                         ijq=(j-1)*iip1+i
@@ -387,18 +420,20 @@
 c   ajout de la maille non completement advectee
                      u_mq(ij,l)=u_mq(ij,l)+zu_m*
-     &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
+     &                  (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
+     &                  *dxq(ijq,l))
                   ELSE
                      ijq=ij+1
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(-zu_m.gt.masse(ijq,l))
-                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
-                        zu_m=zu_m+masse(ijq,l)
+                     do while(-zu_m.gt.masse(ijq,l,iq))
+                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
+     &                          *masse(ijq,l,iq)
+                        zu_m=zu_m+masse(ijq,l,iq)
                         i=mod(i,iim)+1
                         ijq=(j-1)*iip1+i
                      ENDDO
 c   ajout de la maille non completement advectee
-                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
-     &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
+                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
+     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
                   ENDIF
                ENDDO
@@ -417,4 +452,26 @@
       ENDDO
 
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+      write(*,*) 'vlsplt 326: iq,nqfils(iq)=',iq,nqfils(iq) 
+      
+      if (nqdesc(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+         DO l=1,llm
+          DO ij=iip2,ip1jm
+           ! On a besoin de q et masse seulement entre iip2 et ip1jm 
+           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+          enddo   
+         enddo
+        enddo !do ifils=1,nqdesc(iq)
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+! end CRisi
+
 
 c   calcul des tENDances
@@ -422,16 +479,34 @@
       DO l=1,llm
          DO ij=iip2+1,ip1jm
-            new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+
+            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
      &      u_mq(ij-1,l)-u_mq(ij,l))
      &      /new_m
-            masse(ij,l)=new_m
+            masse(ij,l,iq)=new_m
          ENDDO
 c   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
          DO ij=iip1+iip1,ip1jm,iip1
-            q(ij-iim,l)=q(ij,l)
-            masse(ij-iim,l)=masse(ij,l)
-         ENDDO
-      ENDDO
+            q(ij-iim,l,iq)=q(ij,l,iq)
+            masse(ij-iim,l,iq)=masse(ij,l,iq)
+         ENDDO
+      ENDDO
+
+      ! retablir les fils en rapport de melange par rapport a l'air:
+      ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio 
+      ! puis on boucle en longitude
+      if (nqdesc(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+         DO l=1,llm
+          DO ij=iip2+1,ip1jm
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+          DO ij=iip1+iip1,ip1jm,iip1
+             q(ij-iim,l,iq2)=q(ij,l,iq2)
+          enddo ! DO ij=ijb+iip1-1,ije,iip1
+         enddo !DO l=1,llm
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
+
 c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
 c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
@@ -440,5 +515,6 @@
       RETURN
       END
-      SUBROUTINE vly(q,pente_max,masse,masse_adv_v)
+      SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq)
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -464,7 +540,8 @@
 c   Arguments:
 c   ----------
-      REAL masse(ip1jmp1,llm),pente_max
+      REAL masse(ip1jmp1,llm,nqtot),pente_max
       REAL masse_adv_v( ip1jm,llm)
-      REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
+      REAL q(ip1jmp1,llm,nqtot), dq( ip1jmp1,llm)
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -491,4 +568,8 @@
       SAVE sinlon,coslon,sinlondlon,coslondlon
       SAVE airej2,airejjm
+
+      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
+      INTEGER ifils,iq2 ! CRisi
+
 c
 c
@@ -497,4 +578,6 @@
       DATA first,testcpu/.true.,.false./
       DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
+
+      write(*,*) 'vly 578: entree, iq=',iq
 
       IF(first) THEN
@@ -529,6 +612,6 @@
 
       DO i = 1, iim
-      airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
-      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
+      airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
+      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
       ENDDO
       qpns   = SSUM( iim,  airescb ,1 ) / airej2
@@ -538,5 +621,5 @@
 
       DO ij=1,ip1jm
-         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
+         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
          adyqv(ij)=abs(dyqv(ij))
       ENDDO
@@ -553,6 +636,6 @@
 
       DO ij=1,iip1
-         dyq(ij,l)=qpns-q(ij+iip1,l)
-         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
+         dyq(ij,l)=qpns-q(ij+iip1,l,iq)
+         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
       ENDDO
 
@@ -675,12 +758,15 @@
       ENDDO
 
+      write(*,*) 'vly 756'
       DO l=1,llm
        DO ij=1,ip1jm
           IF(masse_adv_v(ij,l).gt.0) THEN
-              qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
-     ,                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
+              qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)*
+     ,                   0.5*(1.-masse_adv_v(ij,l)
+     ,                   /masse(ij+iip1,l,iq))
           ELSE
-              qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
-     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
+              qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)*
+     ,                   0.5*(1.+masse_adv_v(ij,l)
+     ,                   /masse(ij,l,iq))
           ENDIF
           qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
@@ -688,12 +774,34 @@
       ENDDO
 
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+      write(*,*) 'vly 689: iq,nqfils(iq)=',iq,nqfils(iq) 
+   
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+         DO l=1,llm
+         DO ij=1,ip1jmp1
+           ! attention, chaque fils doit avoir son masseq, sinon, le 1er
+           ! fils ecrase le masseq de ses freres.
+           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
+          enddo   
+         enddo
+        enddo !do ifils=1,nqdesc(iq)
+
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         call vly(Ratio,pente_max,masseq,qbyv,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
 
       DO l=1,llm
          DO ij=iip2,ip1jm
-            newmasse=masse(ij,l)
+            newmasse=masse(ij,l,iq)
      &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
-     &         /newmasse
-            masse(ij,l)=newmasse
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
+     &         -qbyv(ij-iip1,l))/newmasse
+            masse(ij,l,iq)=newmasse
          ENDDO
 c.-. ancienne version
@@ -703,12 +811,12 @@
          convpn=SSUM(iim,qbyv(1,l),1)
          convmpn=ssum(iim,masse_adv_v(1,l),1)
-         massepn=ssum(iim,masse(1,l),1)
+         massepn=ssum(iim,masse(1,l,iq),1)
          qpn=0.
          do ij=1,iim
-            qpn=qpn+masse(ij,l)*q(ij,l)
+            qpn=qpn+masse(ij,l,iq)*q(ij,l,iq)
          enddo
          qpn=(qpn+convpn)/(massepn+convmpn)
          do ij=1,iip1
-            q(ij,l)=qpn
+            q(ij,l,iq)=qpn
          enddo
 
@@ -718,12 +826,12 @@
          convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
          convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
-         masseps=ssum(iim, masse(ip1jm+1,l),1)
+         masseps=ssum(iim, masse(ip1jm+1,l,iq),1)
          qps=0.
          do ij = ip1jm+1,ip1jmp1-1
-            qps=qps+masse(ij,l)*q(ij,l)
+            qps=qps+masse(ij,l,iq)*q(ij,l,iq)
          enddo
          qps=(qps+convps)/(masseps+convmps)
          do ij=ip1jm+1,ip1jmp1
-            q(ij,l)=qps
+            q(ij,l,iq)=qps
          enddo
 
@@ -739,5 +847,5 @@
 c        DO ij = 1,iip1
 c           q(ij,l)=newq
-c           masse(ij,l)=newmasse*aire(ij)
+c           masse(ij,l,iq)=newmasse*aire(ij)
 c        ENDDO
 c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
@@ -749,12 +857,27 @@
 c        DO ij = ip1jm+1,ip1jmp1
 c           q(ij,l)=newq
-c           masse(ij,l)=newmasse*aire(ij)
+c           masse(ij,l,iq)=newmasse*aire(ij)
 c        ENDDO
 c._. fin nouvelle version
       ENDDO
+  
+! retablir les fils en rapport de melange par rapport a l'air:
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+         DO l=1,llm
+          DO ij=1,ip1jmp1
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+         enddo
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
+
+      write(*,*) 'vly 853: sortie'
 
       RETURN
       END
-      SUBROUTINE vlz(q,pente_max,masse,w)
+      SUBROUTINE vlz(q,pente_max,masse,w,iq)
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -779,7 +902,8 @@
 c   Arguments:
 c   ----------
-      REAL masse(ip1jmp1,llm),pente_max
-      REAL q(ip1jmp1,llm)
+      REAL masse(ip1jmp1,llm,nqtot),pente_max
+      REAL q(ip1jmp1,llm,nqtot)
       REAL w(ip1jmp1,llm+1)
+      INTEGER iq
 c
 c      Local 
@@ -792,4 +916,7 @@
       REAL dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm),dzqmax
       REAL sigw
+
+      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
+      INTEGER ifils,iq2 ! CRisi
 
       LOGICAL testcpu
@@ -805,4 +932,6 @@
 c    On oriente tout dans le sens de la pression c'est a dire dans le
 c    sens de W
+
+      write(*,*) 'vlz 923: entree'
 
 #ifdef BIDON
@@ -813,5 +942,5 @@
       DO l=2,llm
          DO ij=1,ip1jmp1
-            dzqw(ij,l)=q(ij,l-1)-q(ij,l)
+            dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq)
             adzqw(ij,l)=abs(dzqw(ij,l))
          ENDDO
@@ -835,4 +964,5 @@
       ENDDO
 
+      write(*,*) 'vlz 954'
       DO ij=1,ip1jmp1
          dzq(ij,1)=0.
@@ -851,12 +981,14 @@
 c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
 
+       write(*,*) 'vlz 969'
        DO l = 1,llm-1
          do  ij = 1,ip1jmp1
           IF(w(ij,l+1).gt.0.) THEN
-             sigw=w(ij,l+1)/masse(ij,l+1)
-             wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
+             sigw=w(ij,l+1)/masse(ij,l+1,iq)
+             wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1,iq)
+     &           +0.5*(1.-sigw)*dzq(ij,l+1))
           ELSE
-             sigw=w(ij,l+1)/masse(ij,l)
-             wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
+             sigw=w(ij,l+1)/masse(ij,l,iq)
+             wq(ij,l+1)=w(ij,l+1)*(q(ij,l,iq)-0.5*(1.+sigw)*dzq(ij,l))
           ENDIF
          ENDDO
@@ -868,13 +1000,46 @@
        ENDDO
 
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+      write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq)
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+         DO l=1,llm
+          DO ij=1,ip1jmp1
+           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)        
+          enddo   
+         enddo
+        enddo !do ifils=1,nqdesc(iq)
+        
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)         
+         call vlz(Ratio,pente_max,masseq,wq,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+! end CRisi  
+
       DO l=1,llm
          DO ij=1,ip1jmp1
-            newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
+            newmasse=masse(ij,l,iq)+w(ij,l+1)-w(ij,l)
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+wq(ij,l+1)-wq(ij,l))
      &         /newmasse
-            masse(ij,l)=newmasse
-         ENDDO
-      ENDDO
-
+            masse(ij,l,iq)=newmasse
+         ENDDO
+      ENDDO
+
+! retablir les fils en rapport de melange par rapport a l'air:
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+         DO l=1,llm
+          DO ij=1,ip1jmp1
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+         enddo
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
+      write(*,*) 'vlsplt 1032'
 
       RETURN
Index: /LMDZ5/trunk/libf/dyn3d/vlspltqs.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d/vlspltqs.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3d/vlspltqs.F	(revision 2270)
@@ -3,5 +3,6 @@
 c
        SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
-     ,                                  p,pk,teta                 )
+     ,                                  p,pk,teta,iq             )
+       USE infotrac, ONLY: nqtot,nqdesc,iqfils
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron 
@@ -35,7 +36,8 @@
       REAL masse(ip1jmp1,llm),pente_max
       REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
-      REAL q(ip1jmp1,llm)
+      REAL q(ip1jmp1,llm,nqtot)
       REAL w(ip1jmp1,llm),pdt
       REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -43,11 +45,12 @@
 c
       INTEGER i,ij,l,j,ii
+      INTEGER ifils,iq2 ! CRisi
 c
       REAL qsat(ip1jmp1,llm)
-      REAL zm(ip1jmp1,llm)
+      REAL zm(ip1jmp1,llm,nqtot)
       REAL mu(ip1jmp1,llm)
       REAL mv(ip1jm,llm)
       REAL mw(ip1jmp1,llm+1)
-      REAL zq(ip1jmp1,llm)
+      REAL zq(ip1jmp1,llm,nqtot)
       REAL temps1,temps2,temps3
       REAL zzpbar, zzw
@@ -116,31 +119,33 @@
       ENDDO
 
-      CALL SCOPY(ijp1llm,q,1,zq,1)
-      CALL SCOPY(ijp1llm,masse,1,zm,1)
+      CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1)
+      CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
+      if (nqdesc(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+        iq2=iqfils(ifils,iq)
+        CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
+       enddo  
+      endif !if (nqfils(iq).gt.0) then
 
 c      call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
-      call vlxqs(zq,pente_max,zm,mu,qsat)
-
+      call vlxqs(zq,pente_max,zm,mu,qsat,iq)
 
 c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
 
-      call vlyqs(zq,pente_max,zm,mv,qsat)
-
+      call vlyqs(zq,pente_max,zm,mv,qsat,iq)
 
 c      call minmaxq(zq,qmin,qmax,'avant vlz     ')
 
-      call vlz(zq,pente_max,zm,mw)
-
+      call vlz(zq,pente_max,zm,mw,iq)
 
 c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
 c     call minmaxq(zm,qmin,qmax,'M avant vlyqs     ')
 
-      call vlyqs(zq,pente_max,zm,mv,qsat)
-
+      call vlyqs(zq,pente_max,zm,mv,qsat,iq)
 
 c     call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
 c     call minmaxq(zm,qmin,qmax,'M avant vlxqs     ')
 
-      call vlxqs(zq,pente_max,zm,mu,qsat)
+      call vlxqs(zq,pente_max,zm,mu,qsat,iq)
 
 c     call minmaxq(zq,qmin,qmax,'apres vlxqs     ')
@@ -150,14 +155,31 @@
       DO l=1,llm
          DO ij=1,ip1jmp1
-           q(ij,l)=zq(ij,l)
+           q(ij,l,iq)=zq(ij,l,iq)
          ENDDO
          DO ij=1,ip1jm+1,iip1
-            q(ij+iim,l)=q(ij,l)
-         ENDDO
-      ENDDO
+            q(ij+iim,l,iq)=q(ij,l,iq)
+         ENDDO
+      ENDDO
+      ! CRisi: aussi pour les fils
+      if (nqdesc(iq).gt.0) then
+      do ifils=1,nqdesc(iq)
+        iq2=iqfils(ifils,iq)
+        DO l=1,llm
+         DO ij=1,ip1jmp1
+           q(ij,l,iq2)=zq(ij,l,iq2)
+         ENDDO
+         DO ij=1,ip1jm+1,iip1
+            q(ij+iim,l,iq2)=q(ij,l,iq2)
+         ENDDO
+        ENDDO
+      enddo !do ifils=1,nqdesc(iq)  
+      endif ! if (nqfils(iq).gt.0) then
+      write(*,*) 'vlspltqs 183: fin de la routine'
 
       RETURN
       END
-      SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat)
+      SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat,iq)
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
+
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -179,8 +201,9 @@
 c   Arguments:
 c   ----------
-      REAL masse(ip1jmp1,llm),pente_max
+      REAL masse(ip1jmp1,llm,nqtot),pente_max
       REAL u_m( ip1jmp1,llm )
-      REAL q(ip1jmp1,llm)
+      REAL q(ip1jmp1,llm,nqtot)
       REAL qsat(ip1jmp1,llm)
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -195,4 +218,8 @@
       REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
       REAL u_mq(ip1jmp1,llm)
+
+      ! CRisi
+      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
+      INTEGER ifils,iq2 ! CRisi
 
       Logical first,testcpu
@@ -227,7 +254,7 @@
          DO l = 1, llm
             DO ij=iip2,ip1jm-1
-               dxqu(ij)=q(ij+1,l)-q(ij,l)
+               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
 c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
-c              sigu(ij)=u_m(ij,l)/masse(ij,l)
+c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
             ENDDO
             DO ij=iip1+iip1,ip1jm,iip1
@@ -281,5 +308,5 @@
          DO l = 1, llm
             DO ij=iip2,ip1jm-1
-               dxqu(ij)=q(ij+1,l)-q(ij,l)
+               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
             ENDDO
             DO ij=iip1+iip1,ip1jm,iip1
@@ -323,11 +350,11 @@
       DO l=1,llm
        DO ij=iip2,ip1jm-1
-          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
-     ,                     1.+u_m(ij,l)/masse(ij+1,l),
+          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
+     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
      ,                     u_m(ij,l))
           zdum(ij,l)=0.5*zdum(ij,l)
           u_mq(ij,l)=cvmgp(
-     ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
-     ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
+     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
+     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
      ,                u_m(ij,l))
           u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
@@ -341,11 +368,11 @@
        DO ij=iip2,ip1jm-1
           IF (u_m(ij,l).gt.0.) THEN
-             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
+             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
              u_mq(ij,l)=u_m(ij,l)*
-     $         min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
+     $         min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
           ELSE
-             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
+             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
              u_mq(ij,l)=u_m(ij,l)*
-     $         min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
+     $         min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
           ENDIF
        ENDDO
@@ -416,7 +443,8 @@
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(zu_m.gt.masse(ijq,l))
-                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
-                        zu_m=zu_m-masse(ijq,l)
+                     do while(zu_m.gt.masse(ijq,l,iq))
+                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
+     &                          *masse(ijq,l,iq)
+                        zu_m=zu_m-masse(ijq,l,iq)
                         i=mod(i-2+iim,iim)+1
                         ijq=(j-1)*iip1+i
@@ -424,18 +452,20 @@
 c   ajout de la maille non completement advectee
                      u_mq(ij,l)=u_mq(ij,l)+zu_m*
-     &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
+     &                  (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
+     &                  *dxq(ijq,l))
                   ELSE
                      ijq=ij+1
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(-zu_m.gt.masse(ijq,l))
-                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
-                        zu_m=zu_m+masse(ijq,l)
+                     do while(-zu_m.gt.masse(ijq,l,iq))
+                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
+     &                          *masse(ijq,l,iq)
+                        zu_m=zu_m+masse(ijq,l,iq)
                         i=mod(i,iim)+1
                         ijq=(j-1)*iip1+i
                      ENDDO
 c   ajout de la maille non completement advectee
-                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
-     &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
+                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
+     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
                   ENDIF
                ENDDO
@@ -454,4 +484,25 @@
       ENDDO
 
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+      write(*,*) 'vlspltqs 326: iq,nqfils(iq)=',iq,nqfils(iq) 
+      
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+         DO l=1,llm
+          DO ij=iip2,ip1jm
+           ! On a besoin de q et masse seulement entre iip2 et ip1jm 
+           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+          enddo   
+         enddo
+        enddo !do ifils=1,nqdesc(iq)
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+! end CRisi
 
 c   calcul des tendances
@@ -459,16 +510,33 @@
       DO l=1,llm
          DO ij=iip2+1,ip1jm
-            new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+
+            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
      &      u_mq(ij-1,l)-u_mq(ij,l))
      &      /new_m
-            masse(ij,l)=new_m
+            masse(ij,l,iq)=new_m
          ENDDO
 c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
          DO ij=iip1+iip1,ip1jm,iip1
-            q(ij-iim,l)=q(ij,l)
-            masse(ij-iim,l)=masse(ij,l)
-         ENDDO
-      ENDDO
+            q(ij-iim,l,iq)=q(ij,l,iq)
+            masse(ij-iim,l,iq)=masse(ij,l,iq)
+         ENDDO
+      ENDDO
+
+      ! retablir les fils en rapport de melange par rapport a l'air:
+      ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio 
+      ! puis on boucle en longitude
+      if (nqdesc(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+         DO l=1,llm
+          DO ij=iip2+1,ip1jm
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+          DO ij=iip1+iip1,ip1jm,iip1
+             q(ij-iim,l,iq2)=q(ij,l,iq2)
+          enddo ! DO ij=ijb+iip1-1,ije,iip1
+         enddo !DO l=1,llm
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
 
 c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
@@ -478,5 +546,6 @@
       RETURN
       END
-      SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat)
+      SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq)
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -502,8 +571,9 @@
 c   Arguments:
 c   ----------
-      REAL masse(ip1jmp1,llm),pente_max
+      REAL masse(ip1jmp1,llm,nqtot),pente_max
       REAL masse_adv_v( ip1jm,llm)
-      REAL q(ip1jmp1,llm)
+      REAL q(ip1jmp1,llm,nqtot)
       REAL qsat(ip1jmp1,llm)
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -529,4 +599,7 @@
       SAVE sinlon,coslon,sinlondlon,coslondlon
       SAVE airej2,airejjm
+
+      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
+      INTEGER ifils,iq2 ! CRisi
 c
 c
@@ -567,6 +640,6 @@
 
       DO i = 1, iim
-      airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
-      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
+      airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
+      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
       ENDDO
       qpns   = SSUM( iim,  airescb ,1 ) / airej2
@@ -576,5 +649,5 @@
 
       DO ij=1,ip1jm
-         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
+         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
          adyqv(ij)=abs(dyqv(ij))
       ENDDO
@@ -591,6 +664,6 @@
 
       DO ij=1,iip1
-         dyq(ij,l)=qpns-q(ij+iip1,l)
-         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
+         dyq(ij,l)=qpns-q(ij+iip1,l,iq)
+         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
       ENDDO
 
@@ -710,9 +783,10 @@
        DO ij=1,ip1jm
          IF( masse_adv_v(ij,l).GT.0. ) THEN
-           qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l )  +
-     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
+           qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq )  +
+     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)
+     ,      /masse(ij+iip1,l,iq)))
          ELSE
-              qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
-     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
+              qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) *
+     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) )
          ENDIF
           qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
@@ -721,11 +795,33 @@
 
 
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+      write(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq) 
+   
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+         DO l=1,llm
+         DO ij=1,ip1jmp1
+           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
+          enddo   
+         enddo
+        enddo !do ifils=1,nqdesc(iq)
+
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2
+         call vly(Ratio,pente_max,masseq,qbyv,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+
       DO l=1,llm
          DO ij=iip2,ip1jm
-            newmasse=masse(ij,l)
+            newmasse=masse(ij,l,iq)
      &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
-     &         /newmasse
-            masse(ij,l)=newmasse
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
+     &         -qbyv(ij-iip1,l))/newmasse
+            masse(ij,l,iq)=newmasse
          ENDDO
 c.-. ancienne version
@@ -733,16 +829,16 @@
          convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
          DO ij = 1,iip1
-            newmasse=masse(ij,l)+convmpn*aire(ij)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
+            newmasse=masse(ij,l,iq)+convmpn*aire(ij)
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convpn*aire(ij))/
      &               newmasse
-            masse(ij,l)=newmasse
+            masse(ij,l,iq)=newmasse
          ENDDO
          convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
          convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
          DO ij = ip1jm+1,ip1jmp1
-            newmasse=masse(ij,l)+convmps*aire(ij)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
+            newmasse=masse(ij,l,iq)+convmps*aire(ij)
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convps*aire(ij))/
      &               newmasse
-            masse(ij,l)=newmasse
+            masse(ij,l,iq)=newmasse
          ENDDO
 c.-. fin ancienne version
@@ -757,5 +853,5 @@
 c        DO ij = 1,iip1
 c           q(ij,l)=newq
-c           masse(ij,l)=newmasse*aire(ij)
+c           masse(ij,l,iq)=newmasse*aire(ij)
 c        ENDDO
 c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
@@ -767,9 +863,24 @@
 c        DO ij = ip1jm+1,ip1jmp1
 c           q(ij,l)=newq
-c           masse(ij,l)=newmasse*aire(ij)
+c           masse(ij,l,iq)=newmasse*aire(ij)
 c        ENDDO
 c._. fin nouvelle version
       ENDDO
 
+      write(*,*) 'vly 866'
+
+! retablir les fils en rapport de melange par rapport a l'air:
+      if (nqdesc(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+         DO l=1,llm
+          DO ij=1,ip1jmp1
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+         enddo
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
+      write(*,*) 'vly 879'
+
       RETURN
       END
Index: /LMDZ5/trunk/libf/dyn3d_common/infotrac.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3d_common/infotrac.F90	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3d_common/infotrac.F90	(revision 2270)
@@ -12,4 +12,7 @@
   INTEGER, SAVE :: nbtr
 
+! CRisi: nb traceurs pères= directement advectés par l'air
+  INTEGER, SAVE :: nqperes
+
 ! Name variables
   CHARACTER(len=20), ALLOCATABLE, DIMENSION(:), SAVE :: tname ! tracer short name for restart and diagnostics
@@ -22,4 +25,11 @@
 !         dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code. 
   INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: niadv ! equivalent dyn / physique
+
+! CRisi: tableaux de fils
+  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: nqfils
+  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: nqdesc ! nombres de fils + nombre de tous les petits fils sur toutes les générations
+  INTEGER, SAVE :: nqdesc_tot
+  INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE    :: iqfils
+  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: iqpere
 
 ! conv_flg(it)=0 : convection desactivated for tracer number it 
@@ -30,4 +40,19 @@
   CHARACTER(len=4),SAVE :: type_trac
   CHARACTER(len=8),DIMENSION(:),ALLOCATABLE, SAVE :: solsym
+   
+    ! CRisi: cas particulier des isotopes
+    LOGICAL,SAVE :: ok_isotopes,ok_iso_verif,ok_isotrac,ok_init_iso
+    INTEGER :: niso_possibles   
+    PARAMETER ( niso_possibles=5)
+    real, DIMENSION (niso_possibles),SAVE :: tnat,alpha_ideal
+    LOGICAL, DIMENSION(niso_possibles),SAVE ::  use_iso
+    INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE ::  iqiso ! donne indice iq en fn de (ixt,phase) 
+    INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  iso_num ! donne numéro iso entre 1 et niso_possibles en fn de nqtot
+    INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  iso_indnum ! donne numéro iso entre 1 et niso effectif en fn de nqtot
+    INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  zone_num ! donne numéro de la zone de tracage en fn de nqtot
+    INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  phase_num ! donne numéro de la zone de tracage en fn de nqtot
+    INTEGER, DIMENSION(niso_possibles), SAVE :: indnum_fn_num ! donne indice entre entre 1 et niso en fonction du numéro d isotope entre 1 et niso_possibles
+    INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE ::  index_trac ! numéro ixt en fn izone, indnum entre 1 et niso
+    INTEGER,SAVE :: niso,ntraceurs_zone,ntraciso
  
 CONTAINS
@@ -63,4 +88,5 @@
 
     CHARACTER(len=15), ALLOCATABLE, DIMENSION(:) :: tnom_0  ! tracer short name
+    CHARACTER(len=15), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name: CRisi
     CHARACTER(len=3), DIMENSION(30) :: descrq
     CHARACTER(len=1), DIMENSION(3)  :: txts
@@ -70,4 +96,8 @@
     INTEGER :: nqtrue  ! number of tracers read from tracer.def, without higer order of moment
     INTEGER :: iq, new_iq, iiq, jq, ierr
+    INTEGER :: ifils,ipere,generation ! CRisi
+    LOGICAL :: continu,nouveau_traceurdef
+    INTEGER :: IOstatus ! gestion de la retrocompatibilite de traceur.def
+    CHARACTER(len=15) :: tchaine    
 
     character(len=*),parameter :: modname="infotrac_init"
@@ -134,4 +164,5 @@
           WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
           READ(90,*) nqtrue
+          write(lunout,*) 'nqtrue=',nqtrue
        ELSE 
           WRITE(lunout,*) trim(modname),': Problem in opening traceur.def'
@@ -192,5 +223,5 @@
 ! Allocate variables depending on nqtrue
 !
-    ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue))
+    ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue),tnom_transp(nqtrue))
 !
 !jyg<
@@ -230,7 +261,47 @@
           ! Continue to read tracer.def
           DO iq=1,nqtrue
-             READ(90,*) hadv(iq),vadv(iq),tnom_0(iq)
-          END DO
+
+             write(*,*) 'infotrac 237: iq=',iq
+             ! CRisi: ajout du nom du fluide transporteur
+             ! mais rester retro compatible
+             READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine
+             write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq)
+             write(lunout,*) 'tchaine=',trim(tchaine)
+             write(*,*) 'infotrac 238: IOstatus=',IOstatus
+             if (IOstatus.ne.0) then
+                CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1)
+             endif
+             ! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un
+             ! espace ou pas au milieu de la chaine.
+             continu=1
+             nouveau_traceurdef=0
+             iiq=1
+             do while (continu)
+                if (tchaine(iiq:iiq).eq.' ') then
+                  nouveau_traceurdef=1
+                  continu=0
+                else if (iiq.lt.LEN_TRIM(tchaine)) then
+                  iiq=iiq+1
+                else
+                  continu=0      
+                endif
+             enddo
+             write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef
+             if (nouveau_traceurdef) then
+                write(lunout,*) 'C''est la nouvelle version de traceur.def'
+                tnom_0(iq)=tchaine(1:iiq-1)
+                tnom_transp(iq)=tchaine(iiq+1:15)
+             else
+                write(lunout,*) 'C''est l''ancienne version de traceur.def'
+                write(lunout,*) 'On suppose que les traceurs sont tous d''air'
+                tnom_0(iq)=tchaine
+                tnom_transp(iq) = 'air'
+             endif
+             write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>'
+             write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>'
+
+          END DO !DO iq=1,nqtrue
           CLOSE(90)  
+
        ELSE ! Without tracer.def, set default values 
          if (planet_type=="earth") then
@@ -239,17 +310,22 @@
           vadv(1) = 14
           tnom_0(1) = 'H2Ov'
+          tnom_transp(1) = 'air'
           hadv(2) = 10
           vadv(2) = 10
           tnom_0(2) = 'H2Ol'
+          tnom_transp(2) = 'air'
           hadv(3) = 10
           vadv(3) = 10
           tnom_0(3) = 'RN'
+          tnom_transp(3) = 'air'
           hadv(4) = 10
           vadv(4) = 10
           tnom_0(4) = 'PB'
+          tnom_transp(4) = 'air'
          else ! default for other planets
           hadv(1) = 10
           vadv(1) = 10
           tnom_0(1) = 'dummy'
+          tnom_transp(1) = 'dummy'
          endif ! of if (planet_type=="earth")
        END IF
@@ -258,5 +334,5 @@
        WRITE(lunout,*) trim(modname),': nombre de traceurs ',nqtrue
        DO iq=1,nqtrue
-          WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq)
+          WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq),tnom_transp(iq)
        END DO
 
@@ -447,11 +523,301 @@
     END DO
 
+
+! CRisi: quels sont les traceurs fils et les traceurs pères.
+! initialiser tous les tableaux d'indices liés aux traceurs familiaux
+! + vérifier que tous les pères sont écrits en premières positions
+    ALLOCATE(nqfils(nqtot),nqdesc(nqtot))    
+    ALLOCATE(iqfils(nqtot,nqtot))   
+    ALLOCATE(iqpere(nqtot))
+    nqperes=0
+    nqfils(:)=0
+    nqdesc(:)=0
+    iqfils(:,:)=0
+    iqpere(:)=0
+    nqdesc_tot=0    
+    DO iq=1,nqtot
+      if (tnom_transp(iq) == 'air') then
+        ! ceci est un traceur père
+        WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un pere'
+        nqperes=nqperes+1
+        iqpere(iq)=0
+      else !if (tnom_transp(iq) == 'air') then
+        ! ceci est un fils. Qui est son père?
+        WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un fils'
+        continu=.true.
+        ipere=1
+        do while (continu)           
+          if (tnom_transp(iq) == tnom_0(ipere)) then
+            ! Son père est ipere
+            WRITE(lunout,*) 'Le traceur',iq,'appele ', &
+      &          trim(tnom_0(iq)),' est le fils de ',ipere,'appele ',trim(tnom_0(ipere))
+            nqfils(ipere)=nqfils(ipere)+1  
+            iqfils(nqfils(ipere),ipere)=iq
+            iqpere(iq)=ipere         
+            continu=.false.
+          else !if (tnom_transp(iq) == tnom_0(ipere)) then
+            ipere=ipere+1 
+            if (ipere.gt.nqtot) then
+                WRITE(lunout,*) 'Le traceur',iq,'appele ', &
+      &          trim(tnom_0(iq)),', est orpelin.'
+                CALL abort_gcm('infotrac_init','Un traceur est orphelin',1)
+            endif !if (ipere.gt.nqtot) then
+          endif !if (tnom_transp(iq) == tnom_0(ipere)) then
+        enddo !do while (continu)
+      endif !if (tnom_transp(iq) == 'air') then
+    enddo !DO iq=1,nqtot
+    WRITE(lunout,*) 'infotrac: nqperes=',nqperes    
+    WRITE(lunout,*) 'nqfils=',nqfils
+    WRITE(lunout,*) 'iqpere=',iqpere
+    WRITE(lunout,*) 'iqfils=',iqfils
+
+! Calculer le nombre de descendants à partir de iqfils et de nbfils
+    DO iq=1,nqtot    
+      generation=0
+      continu=.true.
+      ifils=iq
+      do while (continu)
+        ipere=iqpere(ifils)
+        if (ipere.gt.0) then
+         nqdesc(ipere)=nqdesc(ipere)+1   
+         nqdesc_tot=nqdesc_tot+1      
+         iqfils(nqdesc(ipere),ipere)=iq
+         ifils=ipere
+         generation=generation+1
+        else !if (ipere.gt.0) then
+         continu=.false.
+        endif !if (ipere.gt.0) then
+      enddo !do while (continu)    
+      WRITE(lunout,*) 'Le traceur ',iq,', appele ',trim(tnom_0(iq)),' est un traceur de generation: ',generation
+    enddo !DO iq=1,nqtot
+    WRITE(lunout,*) 'infotrac: nqdesc=',nqdesc
+    WRITE(lunout,*) 'iqfils=',iqfils
+    WRITE(lunout,*) 'nqdesc_tot=',nqdesc_tot
+
+! Interdire autres schémas que 10 pour les traceurs fils, et autres schémas
+! que 10 et 14 si des pères ont des fils
+    do iq=1,nqtot
+      if (iqpere(iq).gt.0) then
+        ! ce traceur a un père qui n'est pas l'air
+        ! Seul le schéma 10 est autorisé
+        if (iadv(iq)/=10) then
+           WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for sons'
+          CALL abort_gcm('infotrac_init','Sons should be advected by scheme 10',1)
+        endif
+        ! Le traceur père ne peut être advecté que par schéma 10 ou 14:
+        IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN
+          WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for fathers'
+          CALL abort_gcm('infotrac_init','Fathers should be advected by scheme 10 ou 14',1)
+        endif !IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN
+     endif !if (iqpere(iq).gt.0) the
+    enddo !do iq=1,nqtot
+
+
+! detecter quels sont les traceurs isotopiques parmi des traceurs
+    call infotrac_isoinit(tnom_0,nqtrue)
+        
 !-----------------------------------------------------------------------
 ! Finalize :
 !
-    DEALLOCATE(tnom_0, hadv, vadv)
+    DEALLOCATE(tnom_0, hadv, vadv,tnom_transp)
 
 
   END SUBROUTINE infotrac_init
 
+  SUBROUTINE infotrac_isoinit(tnom_0,nqtrue)
+
+#ifdef CPP_IOIPSL
+  use IOIPSL
+#else
+  ! if not using IOIPSL, we still need to use (a local version of) getin
+  use ioipsl_getincom
+#endif
+  implicit none
+ 
+    ! inputs
+    INTEGER nqtrue
+    CHARACTER(len=15) tnom_0(nqtrue)
+    
+    ! locals    
+    CHARACTER(len=3), DIMENSION(niso_possibles) :: tnom_iso
+    INTEGER, ALLOCATABLE,DIMENSION(:,:) :: nb_iso,nb_traciso
+    INTEGER, ALLOCATABLE,DIMENSION(:) :: nb_isoind
+    INTEGER :: ntraceurs_zone_prec,iq,phase,ixt,iiso,izone
+    CHARACTER(len=19) :: tnom_trac
+    INCLUDE "iniprint.h"
+
+    tnom_iso=(/'eau','HDO','O18','O17','HTO'/)
+
+    ALLOCATE(nb_iso(niso_possibles,nqo))
+    ALLOCATE(nb_isoind(nqo))
+    ALLOCATE(nb_traciso(niso_possibles,nqo))
+    ALLOCATE(iso_num(nqtot))
+    ALLOCATE(iso_indnum(nqtot))
+    ALLOCATE(zone_num(nqtot))
+    ALLOCATE(phase_num(nqtot))
+     
+    iso_num(:)=0
+    iso_indnum(:)=0
+    zone_num(:)=0
+    phase_num(:)=0
+    indnum_fn_num(:)=0
+    use_iso(:)=.false.  
+    nb_iso(:,:)=0  
+    nb_isoind(:)=0     
+    nb_traciso(:,:)=0
+    niso=0
+    ntraceurs_zone=0  
+    ntraceurs_zone_prec=0
+    ntraciso=0
+
+    do iq=nqo+1,nqtot
+       write(lunout,*) 'infotrac 569: iq,tnom_0(iq)=',iq,tnom_0(iq)
+       do phase=1,nqo    
+        do ixt= 1,niso_possibles   
+         tnom_trac=trim(tnom_0(phase))//'_'
+         tnom_trac=trim(tnom_trac)//trim(tnom_iso(ixt))
+         write(*,*) 'phase,ixt,tnom_trac=',phase,ixt,tnom_trac      
+         IF (tnom_0(iq) == tnom_trac) then
+          write(lunout,*) 'Ce traceur est un isotope'
+          nb_iso(ixt,phase)=nb_iso(ixt,phase)+1   
+          nb_isoind(phase)=nb_isoind(phase)+1   
+          iso_num(iq)=ixt
+          iso_indnum(iq)=nb_isoind(phase)
+          indnum_fn_num(ixt)=iso_indnum(iq)
+          phase_num(iq)=phase
+          write(lunout,*) 'iso_num(iq)=',iso_num(iq)
+          write(lunout,*) 'iso_indnum(iq)=',iso_indnum(iq)
+          write(lunout,*) 'indnum_fn_num(ixt)=',indnum_fn_num(ixt)
+          write(lunout,*) 'phase_num(iq)=',phase_num(iq)
+          goto 20
+         else if (iqpere(iq).gt.0) then          
+          if (tnom_0(iqpere(iq)) == tnom_trac) then
+           write(lunout,*) 'Ce traceur est le fils d''un isotope'
+           ! c'est un traceur d'isotope
+           nb_traciso(ixt,phase)=nb_traciso(ixt,phase)+1
+           iso_num(iq)=ixt
+           iso_indnum(iq)=indnum_fn_num(ixt)
+           zone_num(iq)=nb_traciso(ixt,phase)
+           phase_num(iq)=phase
+           write(lunout,*) 'iso_num(iq)=',iso_num(iq)
+           write(lunout,*) 'phase_num(iq)=',phase_num(iq)
+           write(lunout,*) 'zone_num(iq)=',zone_num(iq)
+           goto 20
+          endif !if (tnom_0(iqpere(iq)) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then
+         endif !IF (tnom_0(iq) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then
+        enddo !do ixt= niso_possibles
+       enddo !do phase=1,nqo
+  20   continue
+      enddo !do iq=1,nqtot
+
+      write(lunout,*) 'iso_num=',iso_num
+      write(lunout,*) 'iso_indnum=',iso_indnum
+      write(lunout,*) 'zone_num=',zone_num  
+      write(lunout,*) 'phase_num=',phase_num
+      write(lunout,*) 'indnum_fn_num=',indnum_fn_num
+
+      do ixt= 1,niso_possibles  
+
+        if (nb_iso(ixt,1).eq.1) then
+          ! on vérifie que toutes les phases ont le même nombre de
+          ! traceurs
+          do phase=2,nqo
+            if (nb_iso(ixt,phase).ne.nb_iso(ixt,1)) then
+              write(lunout,*) 'ixt,phase,nb_iso=',ixt,phase,nb_iso(ixt,phase)
+              CALL abort_gcm('infotrac_init','Phases must have same number of isotopes',1)
+            endif
+          enddo !do phase=2,nqo
+
+          niso=niso+1
+          use_iso(ixt)=.true.
+          ntraceurs_zone=nb_traciso(ixt,1)
+
+          ! on vérifie que toutes les phases ont le même nombre de
+          ! traceurs
+          do phase=2,nqo
+            if (nb_traciso(ixt,phase).ne.ntraceurs_zone) then
+              write(lunout,*) 'ixt,phase,nb_traciso=',ixt,phase,nb_traciso(ixt,phase)
+              write(lunout,*) 'ntraceurs_zone=',ntraceurs_zone
+              CALL abort_gcm('infotrac_init','Phases must have same number of tracers',1)
+            endif  
+          enddo  !do phase=2,nqo
+          ! on vérifie que tous les isotopes ont le même nombre de
+          ! traceurs
+          if (ntraceurs_zone_prec.gt.0) then               
+            if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
+              ntraceurs_zone_prec=ntraceurs_zone
+            else !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
+              write(*,*) 'ntraceurs_zone_prec,ntraceurs_zone=',ntraceurs_zone_prec,ntraceurs_zone   
+              CALL abort_gcm('infotrac_init', &
+               &'Isotope tracers are not well defined in traceur.def',1)           
+            endif !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
+           endif !if (ntraceurs_zone_prec.gt.0) then
+
+        else if (nb_iso(ixt,1).ne.0) then
+           WRITE(lunout,*) 'nqo,ixt=',nqo,ixt
+           WRITE(lunout,*) 'nb_iso(ixt,1)=',nb_iso(ixt,1)    
+           CALL abort_gcm('infotrac_init','Isotopes are not well defined in traceur.def',1)     
+        endif   !if (nb_iso(ixt,1).eq.1) then       
+    enddo ! do ixt= niso_possibles
+
+    ! dimensions isotopique:
+    ntraciso=niso*(ntraceurs_zone+1)
+    WRITE(lunout,*) 'niso=',niso
+    WRITE(lunout,*) 'ntraceurs_zone,ntraciso=',ntraceurs_zone,ntraciso    
+ 
+    ! flags isotopiques:
+    if (niso.gt.0) then
+        ok_isotopes=.true.
+    else
+        ok_isotopes=.false.
+    endif 
+    WRITE(lunout,*) 'ok_isotopes=',ok_isotopes
+ 
+    if (ok_isotopes) then
+        ok_iso_verif=.false.
+        call getin('ok_iso_verif',ok_iso_verif)
+        ok_init_iso=.false.
+        call getin('ok_init_iso',ok_init_iso)
+        tnat=(/1.0,155.76e-6,2005.2e-6,0.004/100.,0.0/)
+        alpha_ideal=(/1.0,1.01,1.006,1.003,1.0/)
+    endif !if (ok_isotopes) then  
+    WRITE(lunout,*) 'ok_iso_verif=',ok_iso_verif
+    WRITE(lunout,*) 'ok_init_iso=',ok_init_iso
+
+    if (ntraceurs_zone.gt.0) then
+        ok_isotrac=.true.
+    else
+        ok_isotrac=.false.
+    endif   
+    WRITE(lunout,*) 'ok_isotrac=',ok_isotrac
+
+    ! remplissage du tableau iqiso(ntraciso,phase)
+    ALLOCATE(iqiso(ntraciso,nqo))   
+    iqiso(:,:)=0     
+    do iq=1,nqtot
+        if (iso_num(iq).gt.0) then
+          ixt=iso_indnum(iq)+zone_num(iq)*niso
+          iqiso(ixt,phase_num(iq))=iq
+        endif
+    enddo
+    WRITE(lunout,*) 'iqiso=',iqiso
+
+    ! replissage du tableau index_trac(ntraceurs_zone,niso)
+    ALLOCATE(index_trac(ntraceurs_zone,niso))  
+    if (ok_isotrac) then
+        do iiso=1,niso
+          do izone=1,ntraceurs_zone 
+             index_trac(izone,iiso)=iiso+izone*niso
+          enddo
+        enddo
+    else !if (ok_isotrac) then      
+        index_trac(:,:)=0.0
+    endif !if (ok_isotrac) then
+    write(lunout,*) 'index_trac=',index_trac    
+
+! Finalize :
+    DEALLOCATE(nb_iso)
+
+  END SUBROUTINE infotrac_isoinit
+
 END MODULE infotrac
Index: /LMDZ5/trunk/libf/dyn3d_common/iso_verif_dyn.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3d_common/iso_verif_dyn.F	(revision 2270)
+++ /LMDZ5/trunk/libf/dyn3d_common/iso_verif_dyn.F	(revision 2270)
@@ -0,0 +1,101 @@
+        function iso_verif_noNaN_nostop(x,err_msg)
+        implicit none
+        ! si x est NaN, on affiche message
+        ! d'erreur et return 1 si erreur
+
+        ! input:
+        real x
+        character*(*) err_msg ! message d''erreur à afficher
+
+        ! output 
+        real borne
+        parameter (borne=1e19)
+        integer iso_verif_noNaN_nostop
+
+        if ((x.gt.-borne).and.(x.lt.borne)) then
+                iso_verif_noNAN_nostop=0
+        else
+            write(*,*) 'erreur detectee par iso_verif_nonNaN:'
+            write(*,*) err_msg
+            write(*,*) 'x=',x
+            iso_verif_noNaN_nostop=1
+        endif      
+
+        return
+        end
+
+        function iso_verif_egalite_nostop
+     :           (a,b,err_msg)
+        implicit none
+        ! compare a et b. Si pas egal, on affiche message
+        ! d'erreur et stoppe
+        ! pour egalite, on verifie erreur absolue et arreur relative
+
+        ! input:
+        real a, b
+        character*(*) err_msg ! message d''erreur à afficher
+
+        ! locals
+        real errmax ! erreur maximale en absolu.
+        real errmaxrel ! erreur maximale en relatif autorisée
+        parameter (errmax=1e-8)
+        parameter (errmaxrel=1e-3)
+
+        ! output
+        integer iso_verif_egalite_nostop
+
+        iso_verif_egalite_nostop=0
+
+        if (abs(a-b).gt.errmax) then
+          if (abs((a-b)/max(max(abs(b),abs(a)),1e-18))
+     :            .gt.errmaxrel) then
+            write(*,*) 'erreur detectee par iso_verif_egalite:'
+            write(*,*) err_msg
+            write(*,*) 'a=',a
+            write(*,*) 'b=',b
+            iso_verif_egalite_nostop=1
+          endif
+        endif     
+        
+        return
+        end        
+
+
+        function iso_verif_aberrant_nostop
+     :           (x,iso,q,err_msg)
+        USE infotrac
+        implicit none
+        
+        ! input:
+        real x,q
+        integer iso ! 2=HDO, 1=O18
+        character*(*) err_msg ! message d''erreur à afficher
+
+        ! locals
+        real qmin,deltaD
+        real deltaDmax,deltaDmin
+        parameter (qmin=1e-11)
+        parameter (deltaDmax=200.0,deltaDmin=-999.9)
+
+        ! output
+        integer iso_verif_aberrant_nostop
+
+        iso_verif_aberrant_nostop=0
+
+        ! verifier que HDO est raisonable
+         if (q.gt.qmin) then
+             deltaD=(x/q/tnat(iso)-1)*1000
+             if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+                  write(*,*) 'erreur detectee par iso_verif_aberrant:'
+                  write(*,*) err_msg
+                  write(*,*) 'q=',q
+                  write(*,*) 'deltaD=',deltaD
+                  write(*,*) 'iso=',iso
+                  iso_verif_aberrant_nostop=1
+             endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+          endif !if (q(i,k,iq).gt.qmin) then
+
+       
+        return
+        end        
+
Index: /LMDZ5/trunk/libf/dyn3dmem/advtrac_loc.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dmem/advtrac_loc.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3dmem/advtrac_loc.F	(revision 2270)
@@ -24,5 +24,5 @@
       USE Vampir
       USE times
-      USE infotrac, ONLY: nqtot, iadv
+      USE infotrac, ONLY: nqtot, iadv, ok_iso_verif
       USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
       USE advtrac_mod, ONLY: finmasse
@@ -82,5 +82,5 @@
 !$OMP THREADPRIVATE(testRequest)
 
-c  test sur l'eventuelle creation de valeurs negatives de la masse
+c  test sur l''eventuelle creation de valeurs negatives de la masse
          ijb=ij_begin
          ije=ij_end
@@ -155,7 +155,13 @@
 c$OMP BARRIER
                  
+          write(*,*) 'advtrac 157: appel de vlspltgen_loc'
           call vlspltgen_loc( q,iadv, 2., massem, wg ,
      *                        pbarug,pbarvg,dtvr,p,
      *                        pk,teta )
+
+          write(*,*) 'advtrac 162: apres appel vlspltgen_loc'
+      if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'advtrac 162')
+      endif !if (ok_iso_verif) then
 
 #ifdef DEBUG_IO     
@@ -356,5 +362,6 @@
 c$OMP END DO
 
-       CALL qminimum_loc( q, 2, finmasse )
+        ! CRisi: on passe nqtot et non nq
+       CALL qminimum_loc( q, nqtot, finmasse )
 
       endif ! of if (planet_type=="earth")
Index: /LMDZ5/trunk/libf/dyn3dmem/caladvtrac_loc.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dmem/caladvtrac_loc.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3dmem/caladvtrac_loc.F	(revision 2270)
@@ -56,5 +56,5 @@
 !$OMP THREADPRIVATE(Request_vanleer)
 
-           
+      write(*,*) 'caladvtrac 58: entree'     
       ijbu=ij_begin
       ijeu=ij_end
@@ -109,4 +109,5 @@
 
       IF ( iadvtr.EQ.iapp_tracvl ) THEN
+      write(*,*) 'caladvtrac 133'
 c$OMP MASTER
         call suspend_timer(timer_caldyn)
@@ -183,7 +184,9 @@
          CALL WriteField_u('wg1',wg_adv)
 #endif        
+      write(*,*) 'caladvtrac 185'  
       CALL advtrac_loc( pbarug_adv,pbarvg_adv,wg_adv, 
      *             p_adv,  massem_adv,q_adv, teta_adv,
-     .             pk_adv)
+     .             pk_adv)     
+      write(*,*) 'caladvtrac 189'
 
 
Index: /LMDZ5/trunk/libf/dyn3dmem/check_isotopes_loc.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dmem/check_isotopes_loc.F	(revision 2270)
+++ /LMDZ5/trunk/libf/dyn3dmem/check_isotopes_loc.F	(revision 2270)
@@ -0,0 +1,197 @@
+        subroutine check_isotopes(q,ijb,ije,err_msg)
+        USE infotrac
+        USE parallel_lmdz
+        implicit none
+
+#include "dimensions.h"
+
+        ! inputs
+        integer ijb,ije ! peut être local et différent de ijb_u,ije_u, ex: dans qminimum
+        real q(ijb_u:ije_u,llm,nqtot)
+        character*(*) err_msg ! message d''erreur à afficher
+
+        ! locals
+        integer ixt,phase,k,i,iq,iiso,izone,ieau,iqeau
+        real xtractot,xiiso
+        real borne
+        real qmin
+        real errmax ! erreur maximale en absolu.
+        real errmaxrel ! erreur maximale en relatif autorisée
+        real deltaDmax,deltaDmin
+        real ridicule
+        parameter (borne=1e19)
+        parameter (errmax=1e-8)
+        parameter (errmaxrel=1e-3)
+        parameter (qmin=1e-11)
+        parameter (deltaDmax=200.0,deltaDmin=-999.9)
+        parameter (ridicule=1e-12)
+        real deltaD
+
+        if (ok_isotopes) then
+
+        ! verifier que rien n'est NaN
+        do ixt=1,ntraciso
+          do phase=1,nqo
+            iq=iqiso(ixt,phase)
+            do k=1,llm
+              DO i = ijb,ije
+                if ((q(i,k,iq).gt.-borne).and.
+     :            (q(i,k,iq).lt.borne)) then
+                else !if ((x(ixt,i,j).gt.-borne).and.
+                  write(*,*) 'erreur detectee par iso_verif_noNaN:'
+                  write(*,*) err_msg
+                  write(*,*) 'q,i,k,iq=',q(i,k,iq),i,k,iq
+                  write(*,*) 'borne=',borne
+                  stop
+                endif  !if ((x(ixt,i,j).gt.-borne).and.
+              enddo !DO i = ijb,ije
+            enddo !do k=1,llm
+          enddo !do phase=1,nqo
+        enddo !do ixt=1,ntraciso
+
+        ! verifier que l'eau normale est OK
+        if (use_iso(1)) then
+          ixt=indnum_fn_num(1)
+          do phase=1,nqo
+            iq=iqiso(ixt,phase)
+            do k=1,llm
+            DO i = ijb,ije  
+              if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
+     :          (abs((q(i,k,phase)-q(i,k,iq))/
+     :           max(max(abs(q(i,k,phase)),abs(q(i,k,iq))),1e-18))
+     :           .gt.errmaxrel)) then
+                  write(*,*) 'erreur detectee par iso_verif_egalite:'
+                  write(*,*) err_msg
+                  write(*,*) 'ixt,phase=',ixt,phase
+                  write(*,*) 'q,iq,i,k=',q(i,k,iq),iq,i,k
+                  write(*,*) 'q(i,k,phase)=',q(i,k,phase)
+                  stop
+              endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
+              ! bidouille pour éviter divergence:
+              q(i,k,iq)= q(i,k,phase) 
+            enddo ! DO i = ijb,ije
+            enddo !do k=1,llm
+          enddo ! do phase=1,nqo 
+        endif !if (use_iso(1)) then
+        
+        ! verifier que HDO est raisonable
+        if (use_iso(2)) then
+          ixt=indnum_fn_num(2)
+          do phase=1,nqo
+            iq=iqiso(ixt,phase)
+            do k=1,llm
+            DO i = ijb,ije
+            if (q(i,k,iq).gt.qmin) then
+             deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(2)-1)*1000
+             if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+                  write(*,*) 'erreur detectee par iso_verif_aberrant:'
+                  write(*,*) err_msg
+                  write(*,*) 'ixt,phase=',ixt,phase
+                  write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k
+                  write(*,*) 'q=',q(i,k,:)
+                  write(*,*) 'deltaD=',deltaD
+                  stop
+             endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+            endif !if (q(i,k,iq).gt.qmin) then
+            enddo !DO i = ijb,ije
+            enddo !do k=1,llm
+          enddo ! do phase=1,nqo 
+        endif !if (use_iso(2)) then
+
+        ! verifier que O18 est raisonable
+        if (use_iso(3)) then
+          ixt=indnum_fn_num(3)
+          do phase=1,nqo
+            iq=iqiso(ixt,phase)
+            do k=1,llm
+            DO i = ijb,ije
+            if (q(i,k,iq).gt.qmin) then
+             deltaD=(q(i,k,iq)/q(i,k,phase)/tnat(3)-1)*1000
+             if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+                  write(*,*) 'erreur detectee iso_verif_aberrant O18:'
+                  write(*,*) err_msg
+                  write(*,*) 'ixt,phase=',ixt,phase
+                  write(*,*) 'q,iq,i,k,=',q(i,k,phase),iq,i,k
+                  write(*,*) 'xt=',q(i,k,:)
+                  write(*,*) 'deltaO18=',deltaD
+                  stop
+             endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
+            endif !if (q(i,k,iq).gt.qmin) then
+            enddo !DO i = ijb,ije
+            enddo !do k=1,llm
+          enddo ! do phase=1,nqo 
+        endif !if (use_iso(2)) then
+
+
+        if (ok_isotrac) then
+
+          if (use_iso(2).and.use_iso(1)) then
+            do izone=1,ntraceurs_zone
+             ixt=index_trac(izone,indnum_fn_num(2))
+             ieau=index_trac(izone,indnum_fn_num(1))
+             do phase=1,nqo
+               iq=iqiso(ixt,phase)
+               iqeau=iqiso(ieau,phase)
+               do k=1,llm
+                DO i = ijb,ije
+                if (q(i,k,iq).gt.qmin) then
+                 deltaD=(q(i,k,iq)/q(i,k,iqeau)/tnat(2)-1)*1000
+                 if ((deltaD.gt.deltaDmax).or.
+     &                   (deltaD.lt.deltaDmin)) then
+                  write(*,*) 'erreur dans iso_verif_aberrant trac:'
+                  write(*,*) err_msg
+                  write(*,*) 'izone,phase=',izone,phase
+                  write(*,*) 'ixt,ieau=',ixt,ieau
+                  write(*,*) 'q,iq,i,k,=',q(i,k,iq),iq,i,k
+                  write(*,*) 'deltaD=',deltaD
+                  stop
+                 endif !if ((deltaD.gt.deltaDmax).or.
+                endif !if (q(i,k,iq).gt.qmin) then
+                enddo !DO i = ijb,ije
+                enddo  ! do k=1,llm
+              enddo ! do phase=1,nqo    
+            enddo !do izone=1,ntraceurs_zone
+          endif !if (use_iso(2).and.use_iso(1)) then
+
+          do iiso=1,niso
+           do phase=1,nqo
+              iq=iqiso(iiso,phase)
+              do k=1,llm
+                DO i = ijb,ije
+                   xtractot=0.0
+                   xiiso=q(i,k,iq)
+                   do izone=1,ntraceurs_zone
+                      iq=iqiso(index_trac(izone,iiso),phase)
+                      xtractot=xtractot+ q(i,k,iq)
+                   enddo !do izone=1,ntraceurs_zone
+                   if ((abs(xtractot-xiiso).gt.errmax).and.
+     :                  (abs(xtractot-xiiso)/
+     :                  max(max(abs(xtractot),abs(xiiso)),1e-18)
+     :                  .gt.errmaxrel)) then
+                  write(*,*) 'erreur detectee par iso_verif_traceurs:'
+                  write(*,*) err_msg
+                  write(*,*) 'iiso,phase=',iiso,phase
+                  write(*,*) 'i,k,=',i,k
+                  write(*,*) 'q(i,k,:)=',q(i,k,:)
+                  stop
+                 endif !if ((abs(q(i,k,phase)-q(i,k,iq)).gt.errmax).and.
+                  
+                 ! bidouille pour éviter divergence:
+                 if (abs(xtractot).gt.ridicule) then
+                   do izone=1,ntraceurs_zone
+                     ixt=index_trac(izone,iiso) 
+                     q(i,k,iq)=q(i,k,iq)/xtractot*xiiso
+                   enddo !do izone=1,ntraceurs_zone                
+                  endif !if ((abs(xtractot).gt.ridicule) then
+                enddo !DO i = ijb,ije
+              enddo !do k=1,llm
+           enddo !do phase=1,nqo
+          enddo !do iiso=1,niso
+
+        endif !if (ok_isotrac) then
+
+        endif ! if (ok_isotopes)
+        
+        end
+
+
Index: /LMDZ5/trunk/libf/dyn3dmem/dynetat0_loc.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dmem/dynetat0_loc.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3dmem/dynetat0_loc.F	(revision 2270)
@@ -366,4 +366,19 @@
            write(lunout,*)"Il est donc initialise a zero"
            q(:,:,iq)=0.
+
+           ! CRisi: pour les isotopes, on peut faire init théorique
+           ! distill de Rayleigh très simplifiée
+           if (ok_isotopes) then
+              if ((iso_num(iq).gt.0).and.(zone_num(iq).eq.0)) then
+                q(:,:,iq)=q(:,:,iqpere(iq))                             &
+     &                   *tnat(iso_num(iq))                             &
+     &                   *(q(:,:,iqpere(iq))/30.e-3)                    &
+     &                   **(alpha_ideal(iso_num(iq))-1)
+              endif
+              if ((iso_num(iq).gt.0).and.(zone_num(iq).eq.1)) then
+                  q(:,:,iq)=q(:,:,iqiso(iso_indnum(iq),phase_num(iq)))
+              endif  
+           endif !if (ok_isotopes) then        
+
         ELSE
 #ifdef NC_DOUBLE
@@ -380,5 +395,9 @@
 
         ENDIF
-      ENDDO
+      ENDDO !DO iq=1,nqtot
+
+      if (ok_iso_verif) then
+         call check_isotopes(q,ijb_u,ije_u,'dynetat0_loc')
+      endif !if (ok_iso_verif) then
 
       DEALLOCATE(q_glo)
Index: /LMDZ5/trunk/libf/dyn3dmem/iniacademic_loc.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3dmem/iniacademic_loc.F90	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3dmem/iniacademic_loc.F90	(revision 2270)
@@ -7,5 +7,6 @@
   use exner_hyb_m, only: exner_hyb
   use exner_milieu_m, only: exner_milieu
-  USE infotrac, ONLY : nqtot
+  USE infotrac, ONLY: nqtot,niso_possibles,ok_isotopes,iqpere,ok_iso_verif,tnat,alpha_ideal, &
+        & iqiso,phase_num,iso_indnum,iso_num,zone_num
   USE control_mod, ONLY: day_step,planet_type
   USE parallel_lmdz, ONLY: ijb_u, ije_u, ijb_v, ije_v
@@ -110,5 +111,5 @@
   ztot0      = 0.
   stot0      = 0.
-  ang0       = 0.
+  ang0       = 0.      
 
   if (llm == 1) then
@@ -269,12 +270,32 @@
         if (planet_type=="earth") then
            ! Earth: first two tracers will be water
+
            do i=1,nqtot
               if (i == 1) q(ijb_u:ije_u,:,i)=1.e-10
               if (i == 2) q(ijb_u:ije_u,:,i)=1.e-15
               if (i.gt.2) q(ijb_u:ije_u,:,i)=0.
+
+              ! CRisi: init des isotopes
+              ! distill de Rayleigh très simplifiée
+              if (ok_isotopes) then
+                if ((iso_num(i).gt.0).and.(zone_num(i).eq.0)) then          
+                   q(ijb_u:ije_u,:,i)=q(ijb_u:ije_u,:,iqpere(i))       &
+      &                  *tnat(iso_num(i))                             &
+      &                  *(q(ijb_u:ije_u,:,iqpere(i))/30.e-3)                              &
+     &                   **(alpha_ideal(iso_num(i))-1)
+                endif                
+                if ((iso_num(i).gt.0).and.(zone_num(i).eq.1)) then
+                  q(ijb_u:ije_u,:,i)=q(ijb_u:ije_u,:,iqiso(iso_indnum(i),phase_num(i)))
+                endif
+              endif !if (ok_isotopes) then
+
            enddo
         else
            q(ijb_u:ije_u,:,:)=0
         endif ! of if (planet_type=="earth")
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'iniacademic_loc')
+        endif !if (ok_iso_verif) then
 
         ! add random perturbation to temperature
Index: /LMDZ5/trunk/libf/dyn3dmem/integrd_loc.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dmem/integrd_loc.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3dmem/integrd_loc.F	(revision 2270)
@@ -11,4 +11,5 @@
       USE write_field
       USE integrd_mod
+      USE infotrac, ONLY: ok_iso_verif ! ajout CRisi
       IMPLICIT NONE
 
@@ -86,5 +87,11 @@
       INTEGER :: ierr
 
+      !write(*,*) 'integrd 88: entree, nq=',nq
 c-----------------------------------------------------------------------
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb,ije,'integrd 92')
+        endif !if (ok_iso_verif) then
+
 c$OMP BARRIER     
       if (pole_nord) THEN
@@ -125,6 +132,8 @@
       DO 2 ij = ijb,ije
        pscr (ij)    = ps0(ij)
-       ps (ij)      = psm1(ij) + dt * dp(ij)
+       ps (ij)      = psm1(ij) + dt * dp(ij)     
+
    2  CONTINUE
+
 c$OMP END DO  
 c$OMP BARRIER
@@ -214,4 +223,5 @@
 c$OMP END MASTER
 c$OMP BARRIER
+      !write(*,*) 'integrd 217'  
 c
 c  ... Calcul  de la nouvelle masse d'air au dernier temps integre t+1 ...
@@ -219,4 +229,5 @@
 
       CALL pression_loc ( ip1jmp1, ap, bp, ps, p )
+
 c$OMP BARRIER
       CALL massdair_loc (     p  , masse         )
@@ -334,8 +345,19 @@
            ENDDO
           ENDDO
+          
 c$OMP END DO NOWAIT
 c$OMP BARRIER
 
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb,ije,'integrd 342')
+        endif !if (ok_iso_verif) then
+
+          !write(*,*) 'integrd 341'
           CALL qminimum_loc( q, nq, deltap )
+          !write(*,*) 'integrd 343'
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb,ije,'integrd 346')
+        endif !if (ok_iso_verif) then
 c
 c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
@@ -387,4 +409,8 @@
       
       ENDIF
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb,ije,'integrd 409')
+        endif !if (ok_iso_verif) then
       
 ! Ehouarn: forget about finvmaold
@@ -404,4 +430,5 @@
 
 15    continue
+          !write(*,*) 'integrd 410'
 
 c$OMP DO SCHEDULE(STATIC)
Index: /LMDZ5/trunk/libf/dyn3dmem/leapfrog_loc.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dmem/leapfrog_loc.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3dmem/leapfrog_loc.F	(revision 2270)
@@ -200,4 +200,8 @@
       LOGICAL, SAVE :: firstcall=.TRUE.
       TYPE(distrib),SAVE :: new_dist
+
+      if (ok_iso_verif) then
+         call check_isotopes(q0,ijb_u,ije_u,'leapfrog204: debut')
+      endif !if (ok_iso_verif) then
       
 c$OMP MASTER
@@ -219,4 +223,8 @@
       itaufinp1 = itaufin +1
 
+      if (ok_iso_verif) then
+        call check_isotopes(q0,ijb_u,ije_u,'leapfrog 226')
+      endif !if (ok_iso_verif) then
+
       itau = 0
       physic=.true.
@@ -231,4 +239,8 @@
       phis=phis0
       q=q0
+
+      if (ok_iso_verif) then
+        call check_isotopes(q,ijb_u,ije_u,'leapfrog 239')
+      endif !if (ok_iso_verif) then
       
 !      iday = day_ini+itau/day_step
@@ -296,5 +308,5 @@
 
    1  CONTINUE ! Matsuno Forward step begins here
-
+      write(*,*) 'leapfrog 298: itau=',itau
       jD_cur = jD_ref + day_ini - day_ref +                             &
      &          itau/day_step
@@ -306,4 +318,7 @@
       endif
 
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 321')
+        endif !if (ok_iso_verif) then
 
 #ifdef CPP_IOIPSL
@@ -384,5 +399,15 @@
 cym      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
 
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 400')
+        endif !if (ok_iso_verif) then
+
    2  CONTINUE ! Matsuno backward or leapfrog step begins here
+
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 402')
+        endif !if (ok_iso_verif) then
 
 c$OMP MASTER
@@ -455,4 +480,8 @@
 c$OMP END MASTER      
 
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 471')
+        endif !if (ok_iso_verif) then
 
 !ym  PAS D'AJUSTEMENT POUR LE MOMENT     
@@ -574,4 +603,7 @@
      
       
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 589')
+        endif !if (ok_iso_verif) then
       
 c-----------------------------------------------------------------------
@@ -635,5 +667,8 @@
       ! compute geopotential phi()
       CALL geopot_loc  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
-
+        
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 651')
+        endif !if (ok_iso_verif) then
       
       call VTb(VTcaldyn)
@@ -644,4 +679,5 @@
 !      CALL FTRACE_REGION_BEGIN("caldyn")
       time = jD_cur + jH_cur 
+
       CALL caldyn_loc 
      $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
@@ -670,10 +706,21 @@
 c   -------------------------------------------------------------
 
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,
+     &           'leapfrog 686: avant caladvtrac')
+        endif !if (ok_iso_verif) then
       
       IF( forward. OR . leapf )  THEN
 ! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
+        write(*,*) 'leapfrog 679: avant CALL caladvtrac_loc'
          CALL caladvtrac_loc(q,pbaru,pbarv,
      *        p, masse, dq,  teta,
      .        flxw,pk, iapptrac)
+
+         write(*,*) 'leapfrog 719'
+         if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,
+     &           'leapfrog 698: apres caladvtrac')
+         endif !if (ok_iso_verif) then
 
 !      do j=1,nqtot
@@ -708,8 +755,19 @@
 !       CALL FTRACE_REGION_BEGIN("integrd")
 
-       CALL integrd_loc ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
+       write(*,*) 'leapfrog 720' 
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 756')
+        endif !if (ok_iso_verif) then
+
+       ! CRisi: pourquoi aller jusqu'à 2 et non pas jusqu'à nqtot?? 
+       CALL integrd_loc ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 ,
      $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis)
 !     $              finvmaold                                    )
 
+       write(*,*) 'leapfrog 724'       
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 762')
+        endif !if (ok_iso_verif) then 
+ 
 !       CALL FTRACE_REGION_END("integrd")
 c$OMP BARRIER
@@ -724,4 +782,9 @@
       call WriteField_u('ps_int',ps)
 #endif    
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 775')
+        endif !if (ok_iso_verif) then
+
 c      do j=1,nqtot
 c        call WriteField_p('q'//trim(int2str(j)),
@@ -1082,4 +1145,9 @@
        ENDIF ! of IF( apphys )
 
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1132')
+        endif !if (ok_iso_verif) then
+        write(*,*) 'leapfrog 1134: iflag_phys=',iflag_phys
+
       IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
 c$OMP MASTER
@@ -1146,9 +1214,12 @@
 
 cc$OMP END PARALLEL
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1196')
+        endif !if (ok_iso_verif) then
 
 c-----------------------------------------------------------------------
 c   dissipation horizontale et verticale  des petites echelles:
 c   ----------------------------------------------------------
-
+      write(*,*) 'leapfrog 1163: apdiss=',apdiss
       IF(apdiss) THEN
       
@@ -1379,5 +1450,9 @@
 c                call abort_gcm(modname,abort_message,0)
 c              ENDIF
-        
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1430')
+        endif !if (ok_iso_verif) then     
+ 
 c   ********************************************************************
 c   ********************************************************************
@@ -1455,4 +1530,8 @@
       ENDIF
       
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1509')
+        endif !if (ok_iso_verif) then
+
       IF ( .NOT.purmats ) THEN
 c       ........................................................
@@ -1526,4 +1605,8 @@
             ENDIF
 
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1584')
+        endif !if (ok_iso_verif) then
+
 c-----------------------------------------------------------------------
 c   ecriture de la bande histoire:
@@ -1562,4 +1645,8 @@
             ENDIF ! of IF (itau.EQ.itaufin)
 
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1624')
+        endif !if (ok_iso_verif) then
+
 c-----------------------------------------------------------------------
 c   gestion de l'integration temporelle:
@@ -1596,4 +1683,9 @@
 
       ELSE ! of IF (.not.purmats)
+
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1664')
+        endif !if (ok_iso_verif) then
 
 c       ........................................................
@@ -1631,4 +1723,9 @@
 
             ELSE ! of IF(forward) i.e. backward step
+
+              
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1698')
+        endif !if (ok_iso_verif) then  
 
               IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
@@ -1683,4 +1780,9 @@
             ENDIF ! of IF (forward)
 
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb_u,ije_u,'leapfrog 1750')
+        endif !if (ok_iso_verif) then
+
       END IF ! of IF(.not.purmats)
 c$OMP MASTER
Index: /LMDZ5/trunk/libf/dyn3dmem/qminimum_loc.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dmem/qminimum_loc.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3dmem/qminimum_loc.F	(revision 2270)
@@ -1,4 +1,5 @@
-      SUBROUTINE qminimum_loc( q,nq,deltap )
+      SUBROUTINE qminimum_loc( q,nqtot,deltap )
       USE parallel_lmdz
+      USE infotrac, ONLY: ok_isotopes,ntraciso,iqiso,ok_iso_verif
       IMPLICIT none
 c
@@ -10,6 +11,6 @@
 #include "comvert.h"
 c
-      INTEGER nq
-      REAL q(ijb_u:ije_u,llm,nq), deltap(ijb_u:ije_u,llm)
+      INTEGER nqtot ! CRisi: on remplace nq par nqtot
+      REAL q(ijb_u:ije_u,llm,nqtot), deltap(ijb_u:ije_u,llm)
 c
       INTEGER iq_vap, iq_liq
@@ -27,4 +28,7 @@
       INTEGER i, k, iq
       REAL zx_defau, zx_abc, zx_pump(ijb_u:ije_u), pompe
+
+      real zx_defau_diag(ijb_u:ije_u,llm,2) 
+      real q_follow(ijb_u:ije_u,llm,2) 
 c
       REAL SSUM
@@ -38,4 +42,6 @@
       INTEGER Index_pump(ij_end-ij_begin+1)
       INTEGER nb_pump
+      INTEGER ixt
+      INTEGER iso_verif_noNaN_nostop
 c
 c Quand l'eau liquide est trop petite (ou negative), on prend
@@ -44,11 +50,26 @@
 c
 
+        write(*,*) 'qminimum 52: entree'
+        if (ok_iso_verif) then
+           call check_isotopes(q,ij_begin,ij_end,'qminimum 52')   
+        endif !if (ok_iso_verif) then     
+
       ijb=ij_begin
       ije=ij_end
 
+      zx_defau_diag(ijb:ije,:,:)=0.0
+      q_follow(ijb:ije,:,:)=q(ijb:ije,:,:)  
+
+      !write(*,*) 'qminimum 57'
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
       DO 1000 k = 1, llm
       DO 1040 i = ijb, ije
             if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then
+
+              if (ok_isotopes) then
+                 zx_defau_diag(i,k,iq_liq)=AMAX1
+     :               ( seuil_liq - q(i,k,iq_liq), 0.0 )
+              endif !if (ok_isotopes) then
+
                q(i,k,iq_vap) = q(i,k,iq_vap) + q(i,k,iq_liq) - seuil_liq
                q(i,k,iq_liq) = seuil_liq
@@ -60,8 +81,10 @@
 c --->  SYNCHRO OPENMP ICI
 
+
 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
+      !write(*,*) 'qminimum 81'
       iq = iq_vap
 c
@@ -70,8 +93,15 @@
 c$OMP DO SCHEDULE(STATIC)
       DO i = ijb, ije
+
          if ( seuil_vap - q(i,k,iq) .gt. 0.d0 ) then
+
+            if (ok_isotopes) then
+              zx_defau_diag(i,k,iq)=AMAX1( seuil_vap - q(i,k,iq), 0.0 )
+            endif !if (ok_isotopes) then
+
             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
@@ -79,8 +109,10 @@
       ENDDO
 c$OMP BARRIER
+
 c
 c Quand il s'agit de la premiere couche au-dessus du sol, on
 c doit imprimer un message d'avertissement (saturation possible).
 c
+      !write(*,*) 'qminimum 106'
       nb_pump=0
 c$OMP DO SCHEDULE(STATIC)
@@ -103,4 +135,114 @@
          ENDDO
       ENDIF
+
+      !write(*,*) 'qminimum 128'
+      if (ok_isotopes) then
+      ! CRisi: traiter de même les traceurs d'eau
+      ! Mais il faut les prendre à l'envers pour essayer de conserver la
+      ! masse.
+      ! 1) pompage dans le sol  
+      ! On suppose que ce pompage se fait sans isotopes -> on ne modifie
+      ! rien ici et on croise les doigts pour que ça ne soit pas trop
+      ! génant
+      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)
+        endif !if (zx_pump(i).gt.0.0) then
+      enddo !DO i = ijb, ije  
+
+      ! 2) transfert de vap vers les couches plus hautes
+      !write(*,*) 'qminimum 139'
+      do k=2,llm
+        DO i = ijb, ije
+          if (zx_defau_diag(i,k,iq_vap).gt.0.0) then             
+              ! on ajoute la vapeur en k              
+              do ixt=1,ntraciso
+               q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap))
+     :              +zx_defau_diag(i,k,iq_vap)
+     :              *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap)
+                
+              if (ok_iso_verif) then
+                if (iso_verif_noNaN_nostop(q(i,k,iqiso(ixt,iq_vap)),
+     :                   '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(i,k,iqiso(ixt,iq_vap))=',
+     :                   q(i,k,iqiso(ixt,iq_vap))
+                   write(*,*) 'zx_defau_diag(i,k,iq_vap)=',
+     :                   zx_defau_diag(i,k,iq_vap)
+                   write(*,*) 'q(i,k-1,iqiso(ixt,iq_vap))=',
+     :                   q(i,k-1,iqiso(ixt,iq_vap))
+                   stop
+                endif
+              endif
+
+              ! et on la retranche en k-1
+               q(i,k-1,iqiso(ixt,iq_vap))=q(i,k-1,iqiso(ixt,iq_vap))
+     :              -zx_defau_diag(i,k,iq_vap)
+     :              *deltap(i,k)/deltap(i,k-1)
+     :              *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap)
+
+               if (ok_iso_verif) then
+                if (iso_verif_noNaN_nostop(q(i,k-1,iqiso(ixt,iq_vap)),
+     :                   '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(i,k,iqiso(ixt,iq_vap))=',
+     :                   q(i,k,iqiso(ixt,iq_vap))
+                   write(*,*) 'zx_defau_diag(i,k,iq_vap)=',
+     :                   zx_defau_diag(i,k,iq_vap)
+                   write(*,*) 'q(i,k-1,iqiso(ixt,iq_vap))=',
+     :                   q(i,k-1,iqiso(ixt,iq_vap))
+                   stop
+                endif
+              endif 
+
+              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)
+     :              *deltap(i,k)/deltap(i,k-1)
+          endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then
+        enddo !DO i = 1, ip1jmp1        
+       enddo !do k=2,llm
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb,ije,'qminimum 168')
+        endif !if (ok_iso_verif) then
+        
+      
+        ! 3) transfert d'eau de la vapeur au liquide
+        !write(*,*) 'qminimum 164'
+        do k=1,llm
+        DO i = ijb, ije
+          if (zx_defau_diag(i,k,iq_liq).gt.0.0) then
+
+              ! on ajoute eau liquide en k en k              
+              do ixt=1,ntraciso
+               q(i,k,iqiso(ixt,iq_liq))=q(i,k,iqiso(ixt,iq_liq))
+     :              +zx_defau_diag(i,k,iq_liq)
+     :              *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap)
+              ! et on la retranche à la vapeur en k
+               q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap))
+     :              -zx_defau_diag(i,k,iq_liq)
+     :              *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap)   
+              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
+        enddo !DO i = 1, ip1jmp1
+       enddo !do k=2,llm  
+
+        if (ok_iso_verif) then
+           call check_isotopes(q,ijb,ije,'qminimum 197')
+        endif !if (ok_iso_verif) then
+
+      endif !if (ok_isotopes) then
+      !write(*,*) 'qminimum 188'
 c
       RETURN
Index: /LMDZ5/trunk/libf/dyn3dmem/vlsplt_loc.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dmem/vlsplt_loc.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3dmem/vlsplt_loc.F	(revision 2270)
@@ -2,5 +2,5 @@
 ! $Id$
 !
-      SUBROUTINE vlx_loc(q,pente_max,masse,u_m,ijb_x,ije_x)
+      SUBROUTINE vlx_loc(q,pente_max,masse,u_m,ijb_x,ije_x,iq)
 
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -14,4 +14,5 @@
 c   --------------------------------------------------------------------
       USE parallel_lmdz
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
       IMPLICIT NONE
 c
@@ -25,8 +26,9 @@
 c   Arguments:
 c   ----------
-      REAL masse(ijb_u:ije_u,llm),pente_max
-      REAL u_m( ijb_u:ije_u,llm ),pbarv( iip1,jjb_v:jje_v,llm)
-      REAL q(ijb_u:ije_u,llm)
-      REAL w(ijb_u:ije_u,llm)
+      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
+      REAL u_m( ijb_u:ije_u,llm),pbarv( iip1,jjb_v:jje_v,llm)
+      REAL q(ijb_u:ije_u,llm,nqtot) ! CRisi: ajout dimension nqtot
+      REAL w(ijb_u:ije_u,llm) 
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -42,4 +44,7 @@
       REAL u_mq(ijb_u:ije_u,llm)
 
+      REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
+      INTEGER ifils,iq2 ! CRisi
+
       Logical extremum
 
@@ -51,4 +56,6 @@
       INTEGER ijb,ije,ijb_x,ije_x
       
+      write(*,*) 'vlsplt 58: entree dans vlx_loc, iq,ijb_x=',
+     &   iq,ijb_x
 c   calcul de la pente a droite et a gauche de la maille
 
@@ -64,5 +71,5 @@
 c   calcul des pentes avec limitation, Van Leer scheme I:
 c   -----------------------------------------------------
-
+      ! on a besoin de q entre ijb et ije
 c   calcul de la pente aux points u
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
@@ -70,7 +77,7 @@
             
             DO ij=ijb,ije-1
-               dxqu(ij)=q(ij+1,l)-q(ij,l)
+               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
 c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
-c              sigu(ij)=u_m(ij,l)/masse(ij,l)
+c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
             ENDDO
             DO ij=ijb+iip1-1,ije,iip1
@@ -126,5 +133,5 @@
          DO l = 1, llm
             DO ij=ijb,ije-1
-               dxqu(ij)=q(ij+1,l)-q(ij,l)
+               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
             ENDDO
             DO ij=ijb+iip1-1,ije,iip1
@@ -147,4 +154,6 @@
       ENDIF ! (pente_max.lt.-1.e-5)
 
+      write(*,*) 'vlx 156: iq,ijb_x=',iq,ijb_x
+
 c   bouclage de la pente en iip1:
 c   -----------------------------
@@ -168,11 +177,11 @@
       DO l=1,llm
        DO ij=ijb,ije-1
-          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
-     ,                     1.+u_m(ij,l)/masse(ij+1,l),
-     ,                     u_m(ij,l))
+          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
+     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
+     ,                     u_m(ij,l,iq))
           zdum(ij,l)=0.5*zdum(ij,l)
           u_mq(ij,l)=cvmgp(
-     ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
-     ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
+     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
+     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
      ,                u_m(ij,l))
           u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
@@ -185,13 +194,16 @@
 c	print*,'Cumule ....'
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        ! on a besoin de masse entre ijb et ije
       DO l=1,llm
        DO ij=ijb,ije-1
-c	print*,'masse(',ij,')=',masse(ij,l)
+c	print*,'masse(',ij,')=',masse(ij,l,iq)
           IF (u_m(ij,l).gt.0.) THEN
-             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
-             u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
+             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
+             u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq)
+     :           +0.5*zdum(ij,l)*dxq(ij,l))
           ELSE
-             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
-             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
+             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
+             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq)
+     :           -0.5*zdum(ij,l)*dxq(ij+1,l))
           ENDIF
        ENDDO
@@ -215,4 +227,5 @@
 c$OMP END DO NOWAIT
 c	print*,'Ok test 1'
+
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
@@ -223,5 +236,5 @@
 c$OMP END DO NOWAIT
 c	 print*,'Ok test 2'
-
+        
 
 c   traitement special pour le cas ou on advecte en longitude plus que le
@@ -247,4 +260,6 @@
 c     &       ,'contenu de la maille : ',n0
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+
+
          DO l=1,llm
             IF(nl(l).gt.0) THEN
@@ -258,5 +273,5 @@
                ENDDO
                niju=iju
-c              PRINT*,'niju,nl',niju,nl(l)
+               PRINT*,'vlx 278, niju,nl',niju,nl(l)
 
 c  traitement des mailles
@@ -270,7 +285,8 @@
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(zu_m.gt.masse(ijq,l))
-                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
-                        zu_m=zu_m-masse(ijq,l)
+                     do while(zu_m.gt.masse(ijq,l,iq))
+                        u_mq(ij,l)=u_mq(ij,l)
+     &                          +q(ijq,l,iq)*masse(ijq,l,iq)
+                        zu_m=zu_m-masse(ijq,l,iq)
                         i=mod(i-2+iim,iim)+1
                         ijq=(j-1)*iip1+i
@@ -278,18 +294,20 @@
 c   ajout de la maille non completement advectee
                      u_mq(ij,l)=u_mq(ij,l)+zu_m*
-     &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
+     &               (q(ijq,l,iq)+0.5*
+     &               (1.-zu_m/masse(ijq,l,iq))*dxq(ijq,l))
                   ELSE
                      ijq=ij+1
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(-zu_m.gt.masse(ijq,l))
-                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
-                        zu_m=zu_m+masse(ijq,l)
+                     do while(-zu_m.gt.masse(ijq,l,iq))
+                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
+     &                           *masse(ijq,l,iq)
+                        zu_m=zu_m+masse(ijq,l,iq)
                         i=mod(i,iim)+1
                         ijq=(j-1)*iip1+i
                      ENDDO
 c   ajout de la maille non completement advectee
-                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
-     &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
+                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
+     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
                   ENDIF
                ENDDO
@@ -299,5 +317,4 @@
 cym      ENDIF  ! n0.gt.0 
 9999    continue
-
 
 c   bouclage en latitude
@@ -311,21 +328,71 @@
 c$OMP END DO NOWAIT
 
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+
+      write(*,*) 'vlsplt 326: iq,ijb_x,nqfils(iq)=',iq,ijb_x,nqfils(iq) 
+
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+         DO l=1,llm
+          DO ij=ijb,ije
+           ! On a besoin de q et masse seulement entre ijb et ije. On ne
+           ! les calcule donc que de ijb à ije
+           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+          enddo   
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         call vlx_loc(Ratio,pente_max,masseq,u_mq,ijb_x,ije_x,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+! end CRisi
+
+      write(*,*) 'vlsplt 360: iq,ijb_x=',iq,ijb_x
+
 c   calcul des tENDances
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
          DO ij=ijb+1,ije
-            new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+
-     &      u_mq(ij-1,l)-u_mq(ij,l))
-     &      /new_m
-            masse(ij,l)=new_m
+            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
+     &        u_mq(ij-1,l)-u_mq(ij,l))
+     &        /new_m
+            masse(ij,l,iq)=new_m
          ENDDO
 c   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
          DO ij=ijb+iip1-1,ije,iip1
-            q(ij-iim,l)=q(ij,l)
-            masse(ij-iim,l)=masse(ij,l)
-         ENDDO
-      ENDDO
-c$OMP END DO NOWAIT
+            q(ij-iim,l,iq)=q(ij,l,iq)
+            masse(ij-iim,l,iq)=masse(ij,l,iq)
+         ENDDO
+      ENDDO
+c$OMP END DO NOWAIT
+      write(*,*) 'vlsplt 380: iq,ijb_x=',iq,ijb_x
+
+! retablir les fils en rapport de melange par rapport a l'air:
+      ! On calcule q entre ijb+1 et ije -> on fait pareil pour ratio 
+      ! puis on boucle en longitude
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb+1,ije
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+          DO ij=ijb+iip1-1,ije,iip1
+             q(ij-iim,l,iq2)=q(ij,l,iq2)
+          enddo ! DO ij=ijb+iip1-1,ije,iip1
+         enddo !DO l=1,llm
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
+
+      write(*,*) 'vlsplt 399: iq,ijb_x=',iq,ijb_x
 c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
 c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
@@ -336,5 +403,5 @@
 
 
-      SUBROUTINE vly_loc(q,pente_max,masse,masse_adv_v)
+      SUBROUTINE vly_loc(q,pente_max,masse,masse_adv_v,iq)
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -349,4 +416,5 @@
 c   --------------------------------------------------------------------
       USE parallel_lmdz
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
       IMPLICIT NONE
 c
@@ -361,7 +429,8 @@
 c   Arguments:
 c   ----------
-      REAL masse(ijb_u:ije_u,llm),pente_max
+      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
       REAL masse_adv_v( ijb_v:ije_v,llm)
-      REAL q(ijb_u:ije_u,llm), dq( ijb_u:ije_u,llm)
+      REAL q(ijb_u:ije_u,llm,nqtot), dq( ijb_u:ije_u,llm)
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -392,4 +461,7 @@
       SAVE airej2,airejjm
 c$OMP THREADPRIVATE(airej2,airejjm)
+
+      REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
+      INTEGER ifils,iq2 ! CRisi
 c
 c
@@ -401,6 +473,11 @@
       INTEGER ijb,ije
 
+      ijb=ij_begin-2*iip1
+      ije=ij_end+2*iip1  
+      if (pole_nord) ijb=ij_begin
+      if (pole_sud)  ije=ij_end
+
       IF(first) THEN
-c         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
+         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
          first=.false.
          do i=2,iip1
@@ -434,5 +511,5 @@
       if (pole_nord) then
         DO i = 1, iim
-          airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
+          airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
         ENDDO
         qpns   = SSUM( iim,  airescb ,1 ) / airej2
@@ -441,11 +518,9 @@
       if (pole_sud) then
         DO i = 1, iim
-          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
+          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
         ENDDO
         qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
       endif
       
-      
-
 c   calcul des pentes aux points v
 
@@ -455,8 +530,12 @@
       if (pole_sud)  ije=ij_end-iip1
       
+      ! on a besoin de q entre ij_begin-2*iip1 et ij_end+2*iip1
+      ! Si pole sud, entre ij_begin-2*iip1 et ij_end
+      ! Si pole Nord, entre ij_begin et ij_end+2*iip1
       DO ij=ijb,ije
-         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
+         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
          adyqv(ij)=abs(dyqv(ij))
       ENDDO
+ 
 
 c   calcul des pentes aux points scalaires
@@ -475,5 +554,5 @@
       IF (pole_nord) THEN
         DO ij=1,iip1
-           dyq(ij,l)=qpns-q(ij+iip1,l)
+           dyq(ij,l)=qpns-q(ij+iip1,l,iq)
         ENDDO
         
@@ -497,5 +576,5 @@
 
         DO ij=1,iip1
-           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
+           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
         ENDDO
 
@@ -633,9 +712,10 @@
        DO ij=ijb,ije
           IF(masse_adv_v(ij,l).gt.0) THEN
-              qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
-     ,                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
+              qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)*
+     ,                   0.5*(1.-masse_adv_v(ij,l)
+     ,                   /masse(ij+iip1,l,iq))
           ELSE
-              qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
-     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
+              qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)*
+     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq))
           ENDIF
           qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
@@ -643,4 +723,33 @@
       ENDDO
 c$OMP END DO NOWAIT
+
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+      write(*,*) 'vly 689: iq,nqfils(iq)=',iq,nqfils(iq) 
+
+      ijb=ij_begin-2*iip1
+      ije=ij_end+2*iip1
+      if (pole_nord) ijb=ij_begin
+      if (pole_sud)  ije=ij_end
+   
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
+         DO l=1,llm
+         DO ij=ijb,ije
+           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
+          enddo   
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         call vly_loc(Ratio,pente_max,masseq,qbyv,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+! end CRisi
       
       ijb=ij_begin
@@ -652,11 +761,15 @@
       DO l=1,llm
          DO ij=ijb,ije
-            newmasse=masse(ij,l)
-     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
-     
-            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
-     &         /newmasse
-            masse(ij,l)=newmasse
-         ENDDO
+            newmasse=masse(ij,l,iq)
+     &         +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
+
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
+     &         -qbyv(ij-iip1,l))/newmasse
+
+            masse(ij,l,iq)=newmasse
+
+         ENDDO
+
+
 c.-. ancienne version
 c        convpn=SSUM(iim,qbyv(1,l),1)/apoln
@@ -665,12 +778,12 @@
            convpn=SSUM(iim,qbyv(1,l),1)
            convmpn=ssum(iim,masse_adv_v(1,l),1)
-           massepn=ssum(iim,masse(1,l),1)
+           massepn=ssum(iim,masse(1,l,iq),1)
            qpn=0.
            do ij=1,iim
-              qpn=qpn+masse(ij,l)*q(ij,l)
+              qpn=qpn+masse(ij,l,iq)*q(ij,l,iq)
            enddo
            qpn=(qpn+convpn)/(massepn+convmpn)
            do ij=1,iip1
-              q(ij,l)=qpn
+              q(ij,l,iq)=qpn
            enddo
          endif
@@ -683,12 +796,12 @@
            convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
            convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
-           masseps=ssum(iim, masse(ip1jm+1,l),1)
+           masseps=ssum(iim, masse(ip1jm+1,l,iq),1)
            qps=0.
            do ij = ip1jm+1,ip1jmp1-1
-              qps=qps+masse(ij,l)*q(ij,l)
+              qps=qps+masse(ij,l,iq)*q(ij,l,iq)
            enddo
            qps=(qps+convps)/(masseps+convmps)
            do ij=ip1jm+1,ip1jmp1
-              q(ij,l)=qps
+              q(ij,l,iq)=qps
            enddo
          endif
@@ -704,5 +817,5 @@
 c        DO ij = 1,iip1
 c           q(ij,l)=newq
-c           masse(ij,l)=newmasse*aire(ij)
+c           masse(ij,l,iq)=newmasse*aire(ij)
 c        ENDDO
 c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
@@ -714,9 +827,29 @@
 c        DO ij = ip1jm+1,ip1jmp1
 c           q(ij,l)=newq
-c           masse(ij,l)=newmasse*aire(ij)
+c           masse(ij,l,iq)=newmasse*aire(ij)
 c        ENDDO
 c._. fin nouvelle version
       ENDDO
 c$OMP END DO NOWAIT
+
+! retablir les fils en rapport de melange par rapport a l'air:
+      ijb=ij_begin
+      ije=ij_end
+!      if (pole_nord) ijb=ij_begin
+!      if (pole_sud)  ije=ij_end
+
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb,ije
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
+
 
       RETURN
@@ -725,5 +858,5 @@
       
       
-      SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x)
+      SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x,iq)
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -739,4 +872,5 @@
       USE parallel_lmdz
       USE vlz_mod
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi  
       IMPLICIT NONE
 c
@@ -750,7 +884,8 @@
 c   Arguments:
 c   ----------
-      REAL masse(ijb_u:ije_u,llm),pente_max
-      REAL q(ijb_u:ije_u,llm)
-      REAL w(ijb_u:ije_u,llm+1)
+      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
+      REAL q(ijb_u:ije_u,llm,nqtot)
+      REAL w(ijb_u:ije_u,llm+1,nqtot)
+      INTEGER iq
 c
 c      Local 
@@ -779,5 +914,9 @@
       LOGICAL,SAVE :: first=.TRUE.
 !$OMP THREADPRIVATE(first)
-      
+
+      !REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
+      ! Ces varibles doivent être déclarées en pointer et en save dans
+      ! vlz_loc si on veut qu'elles soient vues par tous les threads.  
+      INTEGER ifils,iq2 ! CRisi
 
       IF (first) THEN
@@ -787,4 +926,5 @@
 c    sens de W
 
+      write(*,*) 'vlsplt 926: entree dans vlz_loc, iq=',iq
 #ifdef BIDON
       IF(testcpu) THEN
@@ -799,5 +939,5 @@
       DO l=2,llm
          DO ij=ijb,ije
-            dzqw(ij,l)=q(ij,l-1)-q(ij,l)
+            dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq)
             adzqw(ij,l)=abs(dzqw(ij,l))
          ENDDO
@@ -842,36 +982,84 @@
 c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
 
+       write(*,*) 'vlz 982,ijb,ije=',ijb,ije
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        DO l = 1,llm-1
          do  ij = ijb,ije
-          IF(w(ij,l+1).gt.0.) THEN
-             sigw=w(ij,l+1)/masse(ij,l+1)
-             wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
+          IF(w(ij,l+1,iq).gt.0.) THEN
+             sigw=w(ij,l+1,iq)/masse(ij,l+1,iq)
+             wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l+1,iq)
+     :           +0.5*(1.-sigw)*dzq(ij,l+1))
           ELSE
-             sigw=w(ij,l+1)/masse(ij,l)
-             wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
+             sigw=w(ij,l+1,iq)/masse(ij,l,iq)
+             wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l,iq)
+     :           -0.5*(1.+sigw)*dzq(ij,l))
           ENDIF
          ENDDO
        ENDDO
-c$OMP END DO NOWAIT
+c$OMP END DO NOWAIT    
+       write(*,*) 'vlz 1001'    
 
 c$OMP MASTER
        DO ij=ijb,ije
-          wq(ij,llm+1)=0.
-          wq(ij,1)=0.
+          wq(ij,llm+1,iq)=0.
+          wq(ij,1,iq)=0.
        ENDDO
 c$OMP END MASTER
 c$OMP BARRIER
 
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+      write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq)
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+         DO l=1,llm
+          DO ij=ijb,ije
+           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+           wq(ij,l,iq2)=wq(ij,l,iq)
+          enddo   
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+c$OMP BARRIER
+        
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         call vlz_loc(Ratio,pente_max,masseq,wq,ijb_x,ije_x,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+! end CRisi  
+
+! CRisi: On rajoute ici une barrière car on veut être sur que tous les
+! wq soient synchronisés
+      write(*,*) 'vlz 1032'  
+c$OMP BARRIER
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
          DO ij=ijb,ije
-            newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
+            newmasse=masse(ij,l,iq)+w(ij,l+1,iq)-w(ij,l,iq)
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)
+     &         +wq(ij,l+1,iq)-wq(ij,l,iq))
      &         /newmasse
-            masse(ij,l)=newmasse
-         ENDDO
-      ENDDO
-c$OMP END DO NOWAIT
+            masse(ij,l,iq)=newmasse
+         ENDDO
+      ENDDO
+c$OMP END DO NOWAIT
+
+! retablir les fils en rapport de melange par rapport a l'air:
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb,ije
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
 
 
Index: /LMDZ5/trunk/libf/dyn3dmem/vlspltgen_loc.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dmem/vlspltgen_loc.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3dmem/vlspltgen_loc.F	(revision 2270)
@@ -27,5 +27,7 @@
       USE Write_Field_loc
       USE VAMPIR
-      USE infotrac, ONLY : nqtot
+      ! CRisi: on rajoute variables utiles d'infotrac  
+      USE infotrac, ONLY : nqtot,nqperes,nqdesc,nqfils,iqfils, 
+     &    ok_iso_verif
       USE vlspltgen_mod
       IMPLICIT NONE
@@ -64,5 +66,5 @@
       REAL ptarg,pdelarg,foeew,zdelta
       REAL tempe(ijb_u:ije_u)
-      INTEGER ijb,ije,iq
+      INTEGER ijb,ije,iq,iq2,ifils
       LOGICAL, SAVE :: firstcall=.TRUE.
 !$OMP THREADPRIVATE(firstcall)
@@ -150,17 +152,21 @@
       ije=ij_end
 
+      DO iq=1,nqtot
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
       DO l=1,llm
          DO ij=ijb,ije
-            mw(ij,l)=w(ij,l) * zzw
+            mw(ij,l,iq)=w(ij,l) * zzw
          ENDDO
       ENDDO
 c$OMP END DO NOWAIT
-
+      ENDDO
+
+      DO iq=1,nqtot  
 c$OMP MASTER
       DO ij=ijb,ije
-         mw(ij,llm+1)=0.
+         mw(ij,llm+1,iq)=0.
       ENDDO
 c$OMP END MASTER
+      ENDDO
 
 c      CALL SCOPY(ijp1llm,q,1,zq,1)
@@ -170,5 +176,5 @@
        ije=ij_end
 
-      DO iq=1,nqtot
+      DO iq=1,nqtot       
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
         DO l=1,llm
@@ -179,5 +185,5 @@
       ENDDO
 
-#ifdef DEBUG_IO    
+#ifdef DEBUG_IO     
        CALL WriteField_u('mu',mu)
        CALL WriteField_v('mv',mv)
@@ -186,11 +192,171 @@
 #endif
 
+      ! verif temporaire
+      ijb=ij_begin-2*iip1
+      ije=ij_end+2*iip1  
+      if (pole_nord) ijb=ij_begin
+      if (pole_sud)  ije=ij_end  
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ijb,ije,'vlspltgen_loc 191')
+      endif !if (ok_iso_verif) then   
+
 c$OMP BARRIER           
-      DO iq=1,nqtot
-
+!      DO iq=1,nqtot
+      DO iq=1,nqperes ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air
+       write(*,*) 'vlspltgen 192: iq,iadv=',iq,iadv(iq)
+#ifdef DEBUG_IO    
+       CALL WriteField_u('zq',zq(:,:,iq))
+       CALL WriteField_u('zm',zm(:,:,iq)) 
+#endif
+        if(iadv(iq) == 0) then
+	
+	  cycle 
+	
+	else if (iadv(iq)==10) then
+
+#ifdef _ADV_HALO        
+! CRisi: on ajoute les nombres de fils et tableaux des fils
+! On suppose qu'on ne peut advecter les fils que par le schéma 10.  
+	  call vlx_loc(zq,pente_max,zm,mu,
+     &	             ij_begin,ij_begin+2*iip1-1,iq)
+          call vlx_loc(zq,pente_max,zm,mu,
+     &               ij_end-2*iip1+1,ij_end,iq)
+#else
+	  call vlx_loc(zq,pente_max,zm,mu,
+     &	             ij_begin,ij_end,iq)
+#endif
+
+c$OMP MASTER
+          call VTb(VTHallo)
+c$OMP END MASTER
+          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
+          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
+! CRisi
+          do ifils=1,nqdesc(iq)
+            iq2=iqfils(ifils,iq)
+            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1)
+            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1)
+          enddo
+
+c$OMP MASTER
+          call VTe(VTHallo)
+c$OMP END MASTER
+	else if (iadv(iq)==14) then
+
+#ifdef _ADV_HALO           
+          call vlxqs_loc(zq,pente_max,zm,mu,
+     &                   qsat,ij_begin,ij_begin+2*iip1-1,iq)
+          call vlxqs_loc(zq,pente_max,zm,mu,
+     &                   qsat,ij_end-2*iip1+1,ij_end,iq)
+#else
+          call vlxqs_loc(zq,pente_max,zm,mu,
+     &                   qsat,ij_begin,ij_end,iq)
+#endif
+
+c$OMP MASTER
+          call VTb(VTHallo)
+c$OMP END MASTER
+
+          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
+          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
+          do ifils=1,nqdesc(iq)
+            iq2=iqfils(ifils,iq)
+            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1)
+            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1)
+          enddo
+
+c$OMP MASTER
+          call VTe(VTHallo)
+c$OMP END MASTER 
+        else
+	
+	  stop 'vlspltgen_p : schema non parallelise'
+      
+        endif
+      
+      enddo !DO iq=1,nqperes
+      
+      
+c$OMP BARRIER      
+c$OMP MASTER      
+      call VTb(VTHallo)
+c$OMP END MASTER
+
+      call SendRequest(MyRequest1)
+
+c$OMP MASTER
+      call VTe(VTHallo)
+c$OMP END MASTER       
+c$OMP BARRIER
+
+      ! verif temporaire
+      ijb=ij_begin-2*iip1
+      ije=ij_end+2*iip1  
+      if (pole_nord) ijb=ij_begin
+      if (pole_sud)  ije=ij_end  
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 280')
+      endif !if (ok_iso_verif) then
+
+      do iq=1,nqperes
+        write(*,*) 'vlspltgen 279: iq=',iq
+
+        if(iadv(iq) == 0) then
+	
+	  cycle 
+	
+	else if (iadv(iq)==10) then
+
+#ifdef _ADV_HALLO
+          call vlx_loc(zq,pente_max,zm,mu,
+     &                 ij_begin+2*iip1,ij_end-2*iip1,iq)
+#endif        
+	else if (iadv(iq)==14) then
+#ifdef _ADV_HALLO
+          call vlxqs_loc(zq,pente_max,zm,mu,
+     &                    qsat,ij_begin+2*iip1,ij_end-2*iip1,iq)
+#endif    
+        else
+	
+	  stop 'vlspltgen_p : schema non parallelise'
+      
+        endif
+      
+      enddo
+c$OMP BARRIER      
+c$OMP MASTER
+      call VTb(VTHallo)
+c$OMP END MASTER
+
+!      call WaitRecvRequest(MyRequest1)
+!      call WaitSendRequest(MyRequest1)
+c$OMP BARRIER
+       call WaitRequest(MyRequest1)
+
+
+c$OMP MASTER
+      call VTe(VTHallo)
+c$OMP END MASTER
+c$OMP BARRIER
+
+      
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 326')
+      endif !if (ok_iso_verif) then        
+      if (ok_iso_verif) then
+           ijb=ij_begin-2*iip1
+           ije=ij_end+2*iip1
+           if (pole_nord) ijb=ij_begin
+           if (pole_sud)  ije=ij_end
+           call check_isotopes(zq,ijb,ije,'vlspltgen_loc 336')
+      endif !if (ok_iso_verif) then  
+
+      do iq=1,nqperes
+       write(*,*) 'vlspltgen 321: iq=',iq
 #ifdef DEBUG_IO    
        CALL WriteField_u('zq',zq(:,:,iq))
        CALL WriteField_u('zm',zm(:,:,iq))
 #endif
+
         if(iadv(iq) == 0) then
 	
@@ -198,47 +364,12 @@
 	
 	else if (iadv(iq)==10) then
-
-#ifdef _ADV_HALO        
-	  call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &	             ij_begin,ij_begin+2*iip1-1)
-          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &               ij_end-2*iip1+1,ij_end)
-#else
-	  call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &	             ij_begin,ij_end)
-#endif
-
-c$OMP MASTER
-          call VTb(VTHallo)
-c$OMP END MASTER
-          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
-          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
-
-c$OMP MASTER
-          call VTe(VTHallo)
-c$OMP END MASTER
+        
+          call vly_loc(zq,pente_max,zm,mv,iq)
+  
 	else if (iadv(iq)==14) then
-
-#ifdef _ADV_HALO           
-          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &                   qsat,ij_begin,ij_begin+2*iip1-1)
-          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &                   qsat,ij_end-2*iip1+1,ij_end)
-#else
-
-          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &                   qsat,ij_begin,ij_end)
-#endif
-
-c$OMP MASTER
-          call VTb(VTHallo)
-c$OMP END MASTER
-
-          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
-          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
-
-c$OMP MASTER
-          call VTe(VTHallo)
-c$OMP END MASTER 
+      
+          call vlyqs_loc(zq,pente_max,zm,mv,
+     &                   qsat,iq)
+ 
         else
 	
@@ -246,78 +377,94 @@
       
         endif
-      
-      enddo
-      
-      
-c$OMP BARRIER      
-c$OMP MASTER      
-      call VTb(VTHallo)
-c$OMP END MASTER
-
-      call SendRequest(MyRequest1)
-
-c$OMP MASTER
-      call VTe(VTHallo)
-c$OMP END MASTER       
-c$OMP BARRIER
-      do iq=1,nqtot
-
-        if(iadv(iq) == 0) then
-	
-	  cycle 
-	
-	else if (iadv(iq)==10) then
-
-#ifdef _ADV_HALLO
-          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &                 ij_begin+2*iip1,ij_end-2*iip1)
-#endif        
-	else if (iadv(iq)==14) then
-#ifdef _ADV_HALLO
-          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &                    qsat,ij_begin+2*iip1,ij_end-2*iip1)
-#endif    
-        else
-	
-	  stop 'vlspltgen_p : schema non parallelise'
-      
-        endif
-      
-      enddo
-c$OMP BARRIER      
-c$OMP MASTER
-      call VTb(VTHallo)
-c$OMP END MASTER
-
-!      call WaitRecvRequest(MyRequest1)
-!      call WaitSendRequest(MyRequest1)
-c$OMP BARRIER
-       call WaitRequest(MyRequest1)
-
-
-c$OMP MASTER
-      call VTe(VTHallo)
-c$OMP END MASTER
-c$OMP BARRIER
-
-
-      do iq=1,nqtot
+       
+       enddo
+
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 357')
+      endif !if (ok_iso_verif) then
+
+      do iq=1,nqperes
+      write(*,*) 'vlspltgen 349: iq=',iq
 #ifdef DEBUG_IO    
        CALL WriteField_u('zq',zq(:,:,iq))
        CALL WriteField_u('zm',zm(:,:,iq))
 #endif
+        if(iadv(iq) == 0) then 
+	  
+	  cycle 
+	
+	else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
+
+c$OMP BARRIER        
+#ifdef _ADV_HALLO
+          call vlz_loc(zq,pente_max,zm,mw,
+     &               ij_begin,ij_begin+2*iip1-1,iq)
+          call vlz_loc(zq,pente_max,zm,mw,
+     &               ij_end-2*iip1+1,ij_end,iq)
+#else
+          call vlz_loc(zq,pente_max,zm,mw,
+     &               ij_begin,ij_end,iq)
+#endif
+c$OMP BARRIER
+
+c$OMP MASTER
+          call VTb(VTHallo)
+c$OMP END MASTER
+
+          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest2)
+          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest2)
+          ! CRisi
+          do ifils=1,nqdesc(iq)
+            iq2=iqfils(ifils,iq)
+            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest2)
+            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest2)
+          enddo      
+c$OMP MASTER
+          call VTe(VTHallo)
+c$OMP END MASTER	
+c$OMP BARRIER
+        else
+	
+	  stop 'vlspltgen_p : schema non parallelise'
+      
+        endif
+      
+      enddo
+c$OMP BARRIER      
+
+c$OMP MASTER        
+      call VTb(VTHallo)
+c$OMP END MASTER
+
+      call SendRequest(MyRequest2)
+
+c$OMP MASTER
+      call VTe(VTHallo)
+c$OMP END MASTER	
+
+
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 429')
+      endif !if (ok_iso_verif) then 
+
+c$OMP BARRIER
+      do iq=1,nqperes
+      write(*,*) 'vlspltgen 409: iq=',iq
+
         if(iadv(iq) == 0) then
-	
+	  
 	  cycle 
 	
-	else if (iadv(iq)==10) then
-        
-          call vly_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv)
-  
-	else if (iadv(iq)==14) then
-      
-          call vlyqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv,
-     &                   qsat)
- 
+	else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
+c$OMP BARRIER        
+
+          write(*,*) 'vlspltgen_loc 461'
+#ifdef _ADV_HALLO
+          write(*,*) 'vlspltgen_loc 462'
+          call vlz_loc(zq,pente_max,zm,mw,
+     &               ij_begin+2*iip1,ij_end-2*iip1,iq)
+#endif
+
+c$OMP BARRIER        
         else
 	
@@ -325,162 +472,101 @@
       
         endif
-       
-       enddo
-
-
-      do iq=1,nqtot
+      
+      enddo
+      write(*,*) 'vlspltgen_loc 476'
+
+c$OMP BARRIER
+c$OMP MASTER
+      call VTb(VTHallo)
+c$OMP END MASTER
+
+!      call WaitRecvRequest(MyRequest2)
+!      call WaitSendRequest(MyRequest2)
+c$OMP BARRIER
+       CALL WaitRequest(MyRequest2)
+
+c$OMP MASTER
+      call VTe(VTHallo)
+c$OMP END MASTER
+c$OMP BARRIER
+
+
+      write(*,*) 'vlspltgen_loc 494'
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 461')
+      endif !if (ok_iso_verif) then
+
+      do iq=1,nqperes
+      write(*,*) 'vlspltgen 449: iq=',iq
 #ifdef DEBUG_IO    
        CALL WriteField_u('zq',zq(:,:,iq))
        CALL WriteField_u('zm',zm(:,:,iq))
 #endif
+        if(iadv(iq) == 0) then
+	
+	  cycle 
+	
+	else if (iadv(iq)==10) then
+        
+          call vly_loc(zq,pente_max,zm,mv,iq)
+  
+	else if (iadv(iq)==14) then
+      
+          call vlyqs_loc(zq,pente_max,zm,mv,
+     &                   qsat,iq)
+ 
+        else
+	
+	  stop 'vlspltgen_p : schema non parallelise'
+      
+        endif
+       
+       enddo !do iq=1,nqperes
+
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 493')
+      endif !if (ok_iso_verif) then
+
+      do iq=1,nqperes
+      write(*,*) 'vlspltgen 477: iq=',iq
+#ifdef DEBUG_IO    
+       CALL WriteField_u('zq',zq(:,:,iq))
+       CALL WriteField_u('zm',zm(:,:,iq))
+#endif
         if(iadv(iq) == 0) then 
 	  
 	  cycle 
 	
-	else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
-
-c$OMP BARRIER        
-#ifdef _ADV_HALLO
-          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
-     &               ij_begin,ij_begin+2*iip1-1)
-          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
-     &               ij_end-2*iip1+1,ij_end)
-#else
-          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
-     &               ij_begin,ij_end)
-#endif
-c$OMP BARRIER
-
-c$OMP MASTER
-          call VTb(VTHallo)
-c$OMP END MASTER
-
-          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest2)
-          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest2)
-
-c$OMP MASTER
-          call VTe(VTHallo)
-c$OMP END MASTER	
-c$OMP BARRIER
-        else
-	
-	  stop 'vlspltgen_p : schema non parallelise'
-      
-        endif
-      
-      enddo
-c$OMP BARRIER      
-
-c$OMP MASTER        
-      call VTb(VTHallo)
-c$OMP END MASTER
-
-      call SendRequest(MyRequest2)
-
-c$OMP MASTER
-      call VTe(VTHallo)
-c$OMP END MASTER	
-
-c$OMP BARRIER
-      do iq=1,nqtot
-
-        if(iadv(iq) == 0) then
-	  
-	  cycle 
-	
-	else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
-c$OMP BARRIER        
-
-#ifdef _ADV_HALLO
-          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
-     &               ij_begin+2*iip1,ij_end-2*iip1)
-#endif
-
-c$OMP BARRIER        
-        else
-	
-	  stop 'vlspltgen_p : schema non parallelise'
-      
-        endif
-      
-      enddo
-
-c$OMP BARRIER
-c$OMP MASTER
-      call VTb(VTHallo)
-c$OMP END MASTER
-
-!      call WaitRecvRequest(MyRequest2)
-!      call WaitSendRequest(MyRequest2)
-c$OMP BARRIER
-       CALL WaitRequest(MyRequest2)
-
-c$OMP MASTER
-      call VTe(VTHallo)
-c$OMP END MASTER
-c$OMP BARRIER
-
-
-      do iq=1,nqtot
-#ifdef DEBUG_IO    
-       CALL WriteField_u('zq',zq(:,:,iq))
-       CALL WriteField_u('zm',zm(:,:,iq))
-#endif
-        if(iadv(iq) == 0) then
-	
-	  cycle 
-	
 	else if (iadv(iq)==10) then
         
-          call vly_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv)
+          call vlx_loc(zq,pente_max,zm,mu,
+     &               ij_begin,ij_end,iq)
   
 	else if (iadv(iq)==14) then
       
-          call vlyqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv,
-     &                   qsat)
+          call vlxqs_loc(zq,pente_max,zm,mu,
+     &                 qsat, ij_begin,ij_end,iq)
  
         else
 	
-	  stop 'vlspltgen_p : schema non parallelise'
+          stop 'vlspltgen_p : schema non parallelise'
       
         endif
        
-       enddo
-
-
-      do iq=1,nqtot
-#ifdef DEBUG_IO    
-       CALL WriteField_u('zq',zq(:,:,iq))
-       CALL WriteField_u('zm',zm(:,:,iq))
-#endif
-        if(iadv(iq) == 0) then 
-	  
-	  cycle 
-	
-	else if (iadv(iq)==10) then
-        
-          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &               ij_begin,ij_end)
-  
-	else if (iadv(iq)==14) then
-      
-          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
-     &                 qsat, ij_begin,ij_end)
- 
-        else
-	
-          stop 'vlspltgen_p : schema non parallelise'
-      
-        endif
-       
-       enddo
-
+       enddo !do iq=1,nqperes
+
+      write(*,*) 'vlspltgen 550: apres derniere serie de call vlx'
+      if (ok_iso_verif) then
+           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 521')
+      endif !if (ok_iso_verif) then
      
       ijb=ij_begin
       ije=ij_end
+      write(*,*) 'vlspltgen_loc 557'
 c$OMP BARRIER      
 
-
+      write(*,*) 'vlspltgen_loc 559'  
       DO iq=1,nqtot
+       write(*,*) 'vlspltgen_loc 561, iq=',iq  
 #ifdef DEBUG_IO    
        CALL WriteField_u('zq',zq(:,:,iq))
@@ -495,5 +581,6 @@
            ENDDO
         ENDDO
-c$OMP END DO NOWAIT          
+c$OMP END DO NOWAIT    
+      write(*,*) 'vlspltgen_loc 575'      
 
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
@@ -504,8 +591,11 @@
         ENDDO
 c$OMP END DO NOWAIT  
-
-      ENDDO
+      write(*,*) 'vlspltgen_loc 583'  
+      ENDDO !DO iq=1,nqtot
         
-      
+      if (ok_iso_verif) then
+           call check_isotopes(q,ij_begin,ij_end,'vlspltgen_loc 557')
+      endif !if (ok_iso_verif) then
+
 c$OMP BARRIER
 
@@ -516,4 +606,5 @@
 cc$OMP BARRIER
 
+      write(*,*) 'vlspltgen 597: sortie'  
       RETURN
       END
Index: /LMDZ5/trunk/libf/dyn3dmem/vlspltgen_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3dmem/vlspltgen_mod.F90	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3dmem/vlspltgen_mod.F90	(revision 2270)
@@ -2,7 +2,7 @@
 
   REAL,POINTER,SAVE :: qsat(:,:)
-  REAL,POINTER,SAVE :: mu(:,:)
+  REAL,POINTER,SAVE :: mu(:,:) ! CRisi: on ajoute une dimension
   REAL,POINTER,SAVE :: mv(:,:)
-  REAL,POINTER,SAVE :: mw(:,:)
+  REAL,POINTER,SAVE :: mw(:,:,:)
   REAL,POINTER,SAVE :: zm(:,:,:)
   REAL,POINTER,SAVE :: zq(:,:,:)
@@ -25,5 +25,5 @@
     CALL allocate_u(mu,llm,d)
     CALL allocate_v(mv,llm,d)
-    CALL allocate_u(mw,llm+1,d)
+    CALL allocate_u(mw,llm+1,nqtot,d)
     CALL allocate_u(zm,llm,nqtot,d)
     CALL allocate_u(zq,llm,nqtot,d)
Index: /LMDZ5/trunk/libf/dyn3dmem/vlspltqs_loc.F
===================================================================
--- /LMDZ5/trunk/libf/dyn3dmem/vlspltqs_loc.F	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3dmem/vlspltqs_loc.F	(revision 2270)
@@ -1,12 +1,13 @@
-      SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x)
+      SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x,iq)
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
 c
 c    ********************************************************************
-c     Shema  d'advection " pseudo amont " .
+c     Shema  d''advection " pseudo amont " .
 c    ********************************************************************
 c
 c   --------------------------------------------------------------------
       USE parallel_lmdz
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi  
       IMPLICIT NONE
 c
@@ -20,8 +21,9 @@
 c   Arguments:
 c   ----------
-      REAL masse(ijb_u:ije_u,llm),pente_max
+      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
       REAL u_m( ijb_u:ije_u,llm )
-      REAL q(ijb_u:ije_u,llm)
+      REAL q(ijb_u:ije_u,llm,nqtot)
       REAL qsat(ijb_u:ije_u,llm)
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -36,4 +38,7 @@
       REAL adxqu(ijb_u:ije_u),dxqmax(ijb_u:ije_u,llm)
       REAL u_mq(ijb_u:ije_u,llm)
+      REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
+      INTEGER ifils,iq2 ! CRisi
+
 
       REAL      SSUM
@@ -42,4 +47,6 @@
       INTEGER ijb,ije,ijb_x,ije_x
       
+      write(*,*) 'vlspltqs 58: entree vlxqs_loc, iq,ijb_x=',
+     &   iq,ijb_x
 
 c   calcul de la pente a droite et a gauche de la maille
@@ -65,7 +72,7 @@
          DO l = 1, llm
             DO ij=ijb,ije-1
-               dxqu(ij)=q(ij+1,l)-q(ij,l)
+               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
 c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
-c              sigu(ij)=u_m(ij,l)/masse(ij,l)
+c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
             ENDDO
             DO ij=ijb+iip1-1,ije,iip1
@@ -120,5 +127,5 @@
          DO l = 1, llm
             DO ij=ijb,ije-1
-               dxqu(ij)=q(ij+1,l)-q(ij,l)
+               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
             ENDDO
             DO ij=ijb+iip1-1,ije,iip1
@@ -179,11 +186,11 @@
       DO l=1,llm
        DO ij=ijb,ije-1
-          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
-     ,                     1.+u_m(ij,l)/masse(ij+1,l),
+          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
+     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
      ,                     u_m(ij,l))
           zdum(ij,l)=0.5*zdum(ij,l)
           u_mq(ij,l)=cvmgp(
-     ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
-     ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
+     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
+     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
      ,                u_m(ij,l))
           u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
@@ -195,16 +202,16 @@
 c   on cumule le flux correspondant a toutes les mailles dont la masse
 c   au travers de la paroi pENDant le pas de temps.
-c   le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
+c   le rapport de melange de l''air advecte est min(q_vanleer, Qsat_downwind)
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
        DO ij=ijb,ije-1
           IF (u_m(ij,l).gt.0.) THEN
-             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
+             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
              u_mq(ij,l)=u_m(ij,l)*
-     $         min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
+     $         min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
           ELSE
-             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
+             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
              u_mq(ij,l)=u_m(ij,l)*
-     $         min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
+     $         min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
           ENDIF
        ENDDO
@@ -273,5 +280,5 @@
                ENDDO
                niju=iju
-c              PRINT*,'niju,nl',niju,nl(l)
+               PRINT*,'vlxqs 280: niju,nl',niju,nl(l)
 
 c  traitement des mailles
@@ -285,26 +292,28 @@
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(zu_m.gt.masse(ijq,l))
-                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
-                        zu_m=zu_m-masse(ijq,l)
+                     do while(zu_m.gt.masse(ijq,l,iq))
+                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
+     &                     *masse(ijq,l,iq) 
+                        zu_m=zu_m-masse(ijq,l,iq)
                         i=mod(i-2+iim,iim)+1
                         ijq=(j-1)*iip1+i
                      ENDDO
 c   ajout de la maille non completement advectee
-                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
-     &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
+                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)
+     &                 +0.5*(1.-zu_m/masse(ijq,l,iq))*dxq(ijq,l))
                   ELSE
                      ijq=ij+1
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(-zu_m.gt.masse(ijq,l))
-                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
-                        zu_m=zu_m+masse(ijq,l)
+                     do while(-zu_m.gt.masse(ijq,l,iq))
+                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
+     &                   *masse(ijq,l,iq)
+                        zu_m=zu_m+masse(ijq,l,iq)
                         i=mod(i,iim)+1
                         ijq=(j-1)*iip1+i
                      ENDDO
 c   ajout de la maille non completement advectee
-                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
-     &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
+                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
+     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
                   ENDIF
                ENDDO
@@ -325,28 +334,77 @@
 c$OMP END DO NOWAIT
 
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+      write(*,*) 'vlspltqs 336: iq,ijb_x,nqfils(iq)=',
+     &     iq,ijb_x,nqfils(iq)  
+
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb,ije
+           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)      
+          enddo   
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqfils(iq)
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         write(*,*) 'vlxqs 349: on appelle vlx pour iq2=',iq2
+         call vlx_loc(Ratio,pente_max,masseq,u_mq,ijb_x,ije_x,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+! end CRisi
+
+      write(*,*) 'vlspltqs 360: iq,ijb_x=',iq,ijb_x   
+
 c   calcul des tendances
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
          DO ij=ijb+1,ije
-            new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+
+            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
      &      u_mq(ij-1,l)-u_mq(ij,l))
      &      /new_m
-            masse(ij,l)=new_m
-         ENDDO
-c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
+            masse(ij,l,iq)=new_m
+         ENDDO
+c   Modif Fred 22 03 96 correction d''un bug (les scopy ci-dessous)
          DO ij=ijb+iip1-1,ije,iip1
-            q(ij-iim,l)=q(ij,l)
-            masse(ij-iim,l)=masse(ij,l)
-         ENDDO
-      ENDDO
-c$OMP END DO NOWAIT
+            q(ij-iim,l,iq)=q(ij,l,iq)
+            masse(ij-iim,l,iq)=masse(ij,l,iq)
+         ENDDO
+      ENDDO
+c$OMP END DO NOWAIT
+
+      write(*,*) 'vlspltqs 380: iq,ijb_x=',iq,ijb_x
+
+! retablir les fils en rapport de melange par rapport a l'air:
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb+1,ije
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+          DO ij=ijb+iip1-1,ije,iip1
+             q(ij-iim,l,iq2)=q(ij,l,iq2)
+          enddo ! DO ij=ijb+iip1-1,ije,iip1
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
+
+      write(*,*) 'vlspltqs 399: iq,ijb_x=',iq,ijb_x
+
 c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
-c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
+c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1,iq),iip1,masse(iip2,1,iq),iip1)
 
 
       RETURN
       END
-      SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat)
+      SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat,iq)
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -361,4 +419,5 @@
 c   --------------------------------------------------------------------
       USE parallel_lmdz
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi  
       IMPLICIT NONE
 c
@@ -373,8 +432,9 @@
 c   Arguments:
 c   ----------
-      REAL masse(ijb_u:ije_u,llm),pente_max
+      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
       REAL masse_adv_v( ijb_v:ije_v,llm)
-      REAL q(ijb_u:ije_u,llm)
+      REAL q(ijb_u:ije_u,llm,nqtot)
       REAL qsat(ijb_u:ije_u,llm)
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -386,5 +446,5 @@
       REAL dyq(ijb_u:ije_u,llm),dyqv(ijb_v:ije_v)
       REAL adyqv(ijb_v:ije_v),dyqmax(ijb_u:ije_u)
-      REAL qbyv(ijb_v:ije_v,llm)
+      REAL qbyv(ijb_v:ije_v,llm,nqtot)
 
       REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
@@ -402,4 +462,7 @@
 c
 c
+      REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
+      INTEGER ifils,iq2 ! CRisi
+
       REAL      SSUM
 
@@ -407,6 +470,18 @@
       INTEGER ijb,ije
 
+      ijb=ij_begin-2*iip1
+      ije=ij_end+2*iip1  
+      if (pole_nord) ijb=ij_begin
+      if (pole_sud)  ije=ij_end
+      ij=3525
+      l=3
+      if ((ij.ge.ijb).and.(ij.le.ije)) then
+        write(*,*) 'vlyqs 480: ij,l,iq,ijb,q(ij,l,:)=',
+     &             ij,l,iq,ijb,q(ij,l,:)
+      endif  
+
       IF(first) THEN
          PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
+         PRINT*,'vlyqs_loc, iq=',iq
          first=.false.
          do i=2,iip1
@@ -439,5 +514,5 @@
       if (pole_nord) then
         DO i = 1, iim
-          airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
+          airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
         ENDDO
         qpns   = SSUM( iim,  airescb ,1 ) / airej2
@@ -446,5 +521,5 @@
       if (pole_sud) then
         DO i = 1, iim
-          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
+          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
         ENDDO
         qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
@@ -460,5 +535,5 @@
       
       DO ij=ijb,ije
-         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
+         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
          adyqv(ij)=abs(dyqv(ij))
       ENDDO
@@ -482,5 +557,5 @@
 c   calcul des pentes aux poles
         DO ij=1,iip1
-           dyq(ij,l)=qpns-q(ij+iip1,l)
+           dyq(ij,l)=qpns-q(ij+iip1,l,iq)
         ENDDO
 
@@ -513,5 +588,5 @@
 
         DO ij=1,iip1
-           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
+           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
         ENDDO
 
@@ -636,14 +711,45 @@
        DO ij=ijb,ije
          IF( masse_adv_v(ij,l).GT.0. ) THEN
-           qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l )  +
-     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
+           qbyv(ij,l,iq)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq )  +
+     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)
+     ,      /masse(ij+iip1,l,iq)))
          ELSE
-              qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
-     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
+              qbyv(ij,l,iq)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) *
+     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) )
          ENDIF
-          qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
+          qbyv(ij,l,iq) = masse_adv_v(ij,l)*qbyv(ij,l,iq)
        ENDDO
       ENDDO
 c$OMP END DO NOWAIT
+
+! CRisi: appel récursif de l'advection sur les fils.
+! Il faut faire ça avant d'avoir mis à jour q et masse
+      write(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq)
+      
+      ijb=ij_begin-2*iip1
+      ije=ij_end+2*iip1
+      if (pole_nord) ijb=ij_begin
+      if (pole_sud)  ije=ij_end  
+
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb,ije
+           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)   
+          enddo   
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         call vly_loc(Ratio,pente_max,masseq,qbyv,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+
+        
+! end CRisi
 
       ijb=ij_begin
@@ -655,9 +761,9 @@
       DO l=1,llm
          DO ij=ijb,ije
-            newmasse=masse(ij,l)
+            newmasse=masse(ij,l,iq)
      &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
-     &         /newmasse
-            masse(ij,l)=newmasse
+            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l,iq)
+     &         -qbyv(ij-iip1,l,iq))/newmasse
+            masse(ij,l,iq)=newmasse
          ENDDO
 c.-. ancienne version
@@ -665,11 +771,11 @@
          IF (pole_nord) THEN
 
-           convpn=SSUM(iim,qbyv(1,l),1)/apoln
+           convpn=SSUM(iim,qbyv(1,l,iq),1)/apoln
            convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
            DO ij = 1,iip1
-              newmasse=masse(ij,l)+convmpn*aire(ij)
-              q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
+              newmasse=masse(ij,l,iq)+convmpn*aire(ij)
+              q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convpn*aire(ij))/
      &                 newmasse
-              masse(ij,l)=newmasse
+              masse(ij,l,iq)=newmasse
            ENDDO
 	 
@@ -678,11 +784,11 @@
 	 IF (pole_sud) THEN
 	 
-	   convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
+	   convps  = -SSUM(iim,qbyv(ip1jm-iim,l,iq),iq,1)/apols
            convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
            DO ij = ip1jm+1,ip1jmp1
-              newmasse=masse(ij,l)+convmps*aire(ij)
-              q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
+              newmasse=masse(ij,l,iq)+convmps*aire(ij)
+              q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convps*aire(ij))/
      &                 newmasse
-              masse(ij,l)=newmasse
+              masse(ij,l,iq)=newmasse
            ENDDO
 	 
@@ -691,27 +797,48 @@
 
 c._. nouvelle version
-c        convpn=SSUM(iim,qbyv(1,l),1)
+c        convpn=SSUM(iim,qbyv(1,l,iq),1)
 c        convmpn=ssum(iim,masse_adv_v(1,l),1)
-c        oldmasse=ssum(iim,masse(1,l),1)
+c        oldmasse=ssum(iim,masse(1,l,iq),1)
 c        newmasse=oldmasse+convmpn
-c        newq=(q(1,l)*oldmasse+convpn)/newmasse
+c        newq=(q(1,l,iq)*oldmasse+convpn)/newmasse
 c        newmasse=newmasse/apoln
 c        DO ij = 1,iip1
-c           q(ij,l)=newq
-c           masse(ij,l)=newmasse*aire(ij)
+c           q(ij,l,iq)=newq
+c           masse(ij,l,iq)=newmasse*aire(ij)
 c        ENDDO
-c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
+c        convps=-SSUM(iim,qbyv(ip1jm-iim,l,iq),1)
 c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
-c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
+c        oldmasse=ssum(iim,masse(ip1jm-iim,l,iq),1)
 c        newmasse=oldmasse+convmps
-c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
+c        newq=(q(ip1jmp1,l,iq)*oldmasse+convps)/newmasse
 c        newmasse=newmasse/apols
 c        DO ij = ip1jm+1,ip1jmp1
-c           q(ij,l)=newq
-c           masse(ij,l)=newmasse*aire(ij)
+c           q(ij,l,iq)=newq
+c           masse(ij,l,iq)=newmasse*aire(ij)
 c        ENDDO
 c._. fin nouvelle version
       ENDDO
 c$OMP END DO NOWAIT
+
+! retablir les fils en rapport de melange par rapport a l'air:
+      ijb=ij_begin
+      ije=ij_end
+!      if (pole_nord) ijb=ij_begin+iip1
+!      if (pole_sud)  ije=ij_end-iip1
+  
+      if (nqfils(iq).gt.0) then  
+       do ifils=1,nqdesc(iq)
+         iq2=iqfils(ifils,iq)  
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
+         DO l=1,llm
+          DO ij=ijb,ije
+            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)            
+          enddo
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+      endif !if (nqfils(iq).gt.0) then
+
+
       RETURN
       END
Index: /LMDZ5/trunk/libf/dyn3dmem/vlz_mod.F90
===================================================================
--- /LMDZ5/trunk/libf/dyn3dmem/vlz_mod.F90	(revision 2269)
+++ /LMDZ5/trunk/libf/dyn3dmem/vlz_mod.F90	(revision 2270)
@@ -1,8 +1,11 @@
 MODULE vlz_mod
 
-  REAL,POINTER,SAVE :: wq(:,:)
+  REAL,POINTER,SAVE :: wq(:,:,:)
   REAL,POINTER,SAVE :: dzq(:,:)
   REAL,POINTER,SAVE :: dzqw(:,:)
   REAL,POINTER,SAVE :: adzqw(:,:)
+  ! CRisi: pour les traceurs:  
+  REAL,POINTER,SAVE :: masseq(:,:,:)
+  REAL,POINTER,SAVE :: Ratio(:,:,:)
   
 CONTAINS
@@ -18,8 +21,12 @@
     
     d=>distrib_vanleer
-    CALL allocate_u(wq,llm+1,d)
+    CALL allocate_u(wq,llm+1,nqtot,d)
     CALL allocate_u(dzq,llm,d)
     CALL allocate_u(dzqw,llm,d)
     CALL allocate_u(adzqw,llm,d)
+    if (nqdesc_tot.gt.0) then
+    CALL allocate_u(masseq,llm,nqtot,d)
+    CALL allocate_u(Ratio,llm,nqtot,d)
+    endif !if (nqdesc_tot.gt.0) then
 
   END SUBROUTINE vlz_allocate
@@ -29,4 +36,5 @@
   USE bands
   USE parallel_lmdz
+  USE infotrac
   IMPLICIT NONE
     TYPE(distrib),INTENT(IN) :: dist
@@ -36,4 +44,9 @@
     CALL switch_u(dzqw,distrib_vanleer,dist)
     CALL switch_u(adzqw,distrib_vanleer,dist)
+    ! CRisi:
+    if (nqdesc_tot.gt.0) then    
+    CALL switch_u(masseq,distrib_vanleer,dist)
+    CALL switch_u(Ratio,distrib_vanleer,dist)
+    endif !if (nqdesc_tot.gt.0) then     
 
   END SUBROUTINE vlz_switch_vanleer  
