Index: trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90	(revision 2272)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90	(revision 2296)
@@ -17,5 +17,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, &
                          force_conserv_tracer
@@ -271,8 +271,11 @@
 
      !ym          call massbar_p(massem,massebx,masseby)  
-
      call vlspltgen_p( q,iadv, 2., massem, wg , &
           pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
 
+          !write(*,*) 'advtrac 162: apres appel vlspltgen_loc'
+      if (ok_iso_verif) then
+           call check_isotopes(q,1,ip1jmp1,'advtrac 162')
+      endif !if (ok_iso_verif) then
 
 !!! ATTENTION !!!! TOUT CE QUI EST ENTRE ICI ET 1234 EST OBSOLETE !!!!!!!
@@ -291,5 +294,6 @@
         if(iadv(iq).eq.10) THEN
 
-           call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
+!           call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
+           call vlsplt_p(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq) !MVals
 
            !   ----------------------------------------------------------------
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/vlsplt_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/vlsplt_p.F	(revision 2272)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/vlsplt_p.F	(revision 2296)
@@ -1,3 +1,3 @@
-      SUBROUTINE vlsplt_p(q,pente_max,masse,w,pbaru,pbarv,pdt)
+      SUBROUTINE vlsplt_p(q,pente_max,masse,w,pbaru,pbarv,pdt,iq)
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -17,4 +17,5 @@
       USE mod_hallo
       USE Vampir
+      USE infotrac, ONLY : nqtot,nqdesc,iqfils ! CRisi
       IMPLICIT NONE
 c
@@ -28,7 +29,7 @@
 c      REAL masse(iip1,jjp1,llm),pente_max
       REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
-      REAL q(ip1jmp1,llm)
-c      REAL q(iip1,jjp1,llm)
+      REAL q(ip1jmp1,llm,nqtot) !CRisi: add dimension nqtot
       REAL w(ip1jmp1,llm),pdt
+      INTEGER iq !CRisi
 c
 c      Local 
@@ -38,9 +39,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 mw(ip1jmp1,llm+1,nqtot)
+      REAL zq(ip1jmp1,llm,nqtot),zz
       REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
       REAL second,temps0,temps1,temps2,temps3
@@ -51,4 +52,5 @@
       SAVE temps1,temps2,temps3
       INTEGER iminn,imaxx
+      INTEGER ifils,iq2 ! CRisi
 
       REAL qmin,qmax
@@ -93,10 +95,10 @@
       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
 
       DO ij=ijb,ije
-         mw(ij,llm+1)=0.
+         mw(ij,llm+1,iq)=0.
       ENDDO
       
@@ -106,12 +108,12 @@
        ijb=ij_begin
        ije=ij_end
-       zq(ijb:ije,:)=q(ijb:ije,:)
-       zm(ijb:ije,:)=masse(ijb:ije,:)
+       zq(ijb:ije,:,iq)=q(ijb:ije,:,iq)
+       zm(ijb:ije,:,iq)=masse(ijb:ije,:)
       
       
 c	print*,'Entree vlx1'
 c	call minmaxq(zq,qmin,qmax,'avant vlx     ')
-      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_begin+2*iip1-1)
-      call vlx_p(zq,pente_max,zm,mu,ij_end-2*iip1+1,ij_end)
+      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_begin+2*iip1-1,iq)
+      call vlx_p(zq,pente_max,zm,mu,ij_end-2*iip1+1,ij_end,iq)
       call VTb(VTHallo)
       call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest1)
@@ -119,5 +121,5 @@
       call SendRequest(MyRequest1)
       call VTe(VTHallo)
-      call vlx_p(zq,pente_max,zm,mu,ij_begin+2*iip1,ij_end-2*iip1)
+      call vlx_p(zq,pente_max,zm,mu,ij_begin+2*iip1,ij_end-2*iip1,iq)
 c      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end)
       call VTb(VTHallo)
@@ -133,9 +135,9 @@
 c      call exchange_hallo(zm,ip1jmp1,llm,1,1)
       
-      call vly_p(zq,pente_max,zm,mv)
+      call vly_p(zq,pente_max,zm,mv,iq)
 c	call minmaxq(zq,qmin,qmax,'apres vly1     ')
 c	print*,'Sortie vly1'
-      call vlz_p(zq,pente_max,zm,mw,ij_begin,ij_begin+2*iip1-1)
-      call vlz_p(zq,pente_max,zm,mw,ij_end-2*iip1+1,ij_end)
+      call vlz_p(zq,pente_max,zm,mw,ij_begin,ij_begin+2*iip1-1,iq)
+      call vlz_p(zq,pente_max,zm,mw,ij_end-2*iip1+1,ij_end,iq)
       call VTb(VTHallo)
       call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest2)
@@ -143,5 +145,5 @@
       call SendRequest(MyRequest2)
       call VTe(VTHallo)
-      call vlz_p(zq,pente_max,zm,mw,ij_begin+2*iip1,ij_end-2*iip1)
+      call vlz_p(zq,pente_max,zm,mw,ij_begin+2*iip1,ij_end-2*iip1,iq)
       call VTb(VTHallo)
       call WaitRecvRequest(MyRequest2)
@@ -154,9 +156,9 @@
       
       
-      call vly_p(zq,pente_max,zm,mv)
+      call vly_p(zq,pente_max,zm,mv,iq)
 c	call minmaxq(zq,qmin,qmax,'apres vly     ')
 
 
-      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end)
+      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end,iq)
 c	call minmaxq(zq,qmin,qmax,'apres vlx2    ')
 
@@ -167,5 +169,5 @@
       DO l=1,llm
          DO ij=ijb,ije
-           q(ij,l)=zq(ij,l)
+           q(ij,l,iq)=zq(ij,l,iq)
          ENDDO
       ENDDO
@@ -174,7 +176,22 @@
       DO l=1,llm
          DO ij=ijb,ije-iip1+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=ijb,ije
+           q(ij,l,iq2)=zq(ij,l,iq2)
+         ENDDO
+         DO ij=ijb,ije-iip1+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
 
       call WaitSendRequest(MyRequest1) 
@@ -185,5 +202,5 @@
       
       
-      SUBROUTINE vlx_p(q,pente_max,masse,u_m,ijb_x,ije_x)
+      RECURSIVE SUBROUTINE vlx_p(q,pente_max,masse,u_m,ijb_x,ije_x,iq)
 
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -197,4 +214,5 @@
 c   --------------------------------------------------------------------
       USE Parallel_lmdz
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
       IMPLICIT NONE
 c
@@ -205,8 +223,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 w(ip1jmp1,llm)
+      REAL q(ip1jmp1,llm,nqtot) ! CRisi: ajout dimension nqtot
+      !REAL w(ip1jmp1,llm,nqtot) !MVals seems not useful in this case
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -222,4 +241,7 @@
       REAL u_mq(ip1jmp1,llm)
 
+      REAL Ratio(ip1jmp1,llm,nqtot) ! CRisi,MVals: Ratio zq(fils)/zq(pere)
+      INTEGER ifils,iq2 ! CRisi,MVals: indices pour les traceurs fils
+
       Logical extremum
 
@@ -238,5 +260,5 @@
       if (pole_nord.and.ijb==1) ijb=ijb+iip1
       if (pole_sud.and.ije==ip1jmp1)  ije=ije-iip1
-         
+        
       IF (pente_max.gt.-1.e-5) THEN
 c       IF (pente_max.gt.10) THEN
@@ -250,7 +272,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
@@ -306,5 +328,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
@@ -326,5 +348,4 @@
 c$OMP END DO NOWAIT
       ENDIF ! (pente_max.lt.-1.e-5)
-
 c   bouclage de la pente en iip1:
 c   -----------------------------
@@ -348,11 +369,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)
@@ -367,11 +388,12 @@
       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
@@ -426,4 +448,5 @@
 c      PRINT*,'Nombre de points pour lesquels on advect plus que le'
 c     &       ,'contenu de la maille : ',n0
+
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO l=1,llm
@@ -450,7 +473,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
@@ -458,18 +482,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
@@ -491,18 +517,46 @@
 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
+      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
+           !MVals: veiller a ce qu'on n'ait pas de denominateur nul
+           masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),1e-16)
+           if (q(ij,l,iq).gt.1e-16) then
+             Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+           else
+             Ratio(ij,l,iq2)=0.
+           endif
+          enddo   
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         call vlx_p(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+! end CRisi
 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)+
+            !MVals: veiller a ce qu'on n'ait pas de denominateur nul
+            new_m=max(masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l),1e-16)
+            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=ijb+iip1-1,ije,iip1
-            q(ij-iim,l)=q(ij,l)
-            masse(ij-iim,l)=masse(ij,l)
+            q(ij-iim,l,iq)=q(ij,l,iq)
+            masse(ij-iim,l,iq)=masse(ij,l,iq)
          ENDDO
       ENDDO
@@ -511,4 +565,22 @@
 c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
 
+! 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
 
       RETURN
@@ -516,5 +588,5 @@
 
 
-      SUBROUTINE vly_p(q,pente_max,masse,masse_adv_v)
+      RECURSIVE SUBROUTINE vly_p(q,pente_max,masse,masse_adv_v,iq)
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -529,4 +601,5 @@
 c   --------------------------------------------------------------------
       USE parallel_lmdz
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
       USE comconst_mod, ONLY: pi
 
@@ -540,7 +613,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) ! CRisi: ajout dimension nqtot
+      INTEGER iq ! CRisi
 c
 c      Local 
@@ -571,4 +645,7 @@
       SAVE airej2,airejjm
 c$OMP THREADPRIVATE(airej2,airejjm)
+
+      REAL Ratio(ip1jmp1,llm,nqtot) ! CRisi,MVals: Ratio zq(fils)/zq(pere)
+      INTEGER ifils,iq2 ! CRisi,MVals: indices pour les traceurs fils
 c
 c
@@ -613,5 +690,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
@@ -620,5 +697,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
@@ -635,5 +712,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
@@ -654,5 +731,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
         
@@ -676,5 +753,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
 
@@ -812,9 +889,9 @@
        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)
@@ -822,4 +899,40 @@
       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
+           !masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+           !MVals: veiller a ce qu'on n'ait pas de denominateur nul
+           masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),1e-16)
+           if (q(ij,l,iq).gt.1e-16) then
+             Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+           else
+             Ratio(ij,l,iq2)=0.
+           endif   
+          enddo   
+         enddo
+c$OMP END DO NOWAIT
+        enddo !do ifils=1,nqdesc(iq)
+
+        do ifils=1,nqfils(iq)
+         iq2=iqfils(ifils,iq)
+         call vly_p(Ratio,pente_max,masse,qbyv,iq2)
+        enddo !do ifils=1,nqfils(iq)
+      endif !if (nqfils(iq).gt.0) then
+! end CRisi
       
       ijb=ij_begin
@@ -831,10 +944,10 @@
       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)-qbyv(ij-iip1,l))/newmasse
+            masse(ij,l,iq)=newmasse
          ENDDO
 c.-. ancienne version
@@ -844,12 +957,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
@@ -862,12 +975,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
@@ -899,4 +1012,23 @@
 c$OMP END DO NOWAIT
 
+! CRisi: 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
       END
@@ -904,5 +1036,5 @@
       
       
-      SUBROUTINE vlz_p(q,pente_max,masse,w,ijb_x,ije_x)
+      RECURSIVE SUBROUTINE vlz_p(q,pente_max,masse,w,ijb_x,ije_x,iq)
 c
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
@@ -917,4 +1049,5 @@
 c   --------------------------------------------------------------------
       USE Parallel_lmdz
+      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
       IMPLICIT NONE
 c
@@ -925,7 +1058,8 @@
 c   Arguments:
 c   ----------
-      REAL masse(ip1jmp1,llm),pente_max
-      REAL q(ip1jmp1,llm)
-      REAL w(ip1jmp1,llm+1)
+      REAL masse(ip1jmp1,llm,nqtot),pente_max
+      REAL q(ip1jmp1,llm,nqtot)
+      REAL w(ip1jmp1,llm+1,nqtot)
+      INTEGER iq
 c
 c      Local 
@@ -934,5 +1068,7 @@
       INTEGER i,ij,l,j,ii
 c
-      REAL,SAVE :: wq(ip1jmp1,llm+1)
+!      REAL,SAVE :: wq(ip1jmp1,llm+1)
+      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: wq !MVals
+      
       REAL newmasse
 
@@ -956,4 +1092,21 @@
 c    On oriente tout dans le sens de la pression c'est a dire dans le
 c    sens de W
+      !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. 
+      !REAL Ratio(ip1jmp1,llm,nqtot) ! CRisi 
+      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: Ratio ! CRisi 
+      INTEGER ifils,iq2 ! CRisi
+      LOGICAL, SAVE :: first=.TRUE.!MVals
+!$OMP THREADPRIVATE(first)
+
+         IF (first) THEN !MVals
+            first=.FALSE.
+!$OMP MASTER
+            ALLOCATE(wq(ip1jmp1,llm+1,nqtot)) !MVals
+            ALLOCATE(Ratio(ip1jmp1,llm,nqtot)) !MVals
+!$OMP END MASTER
+!$OMP BARRIER
+         END IF !MVals
 
 #ifdef BIDON
@@ -969,5 +1122,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
@@ -999,4 +1152,5 @@
          dzq(ij,llm)=0.
       ENDDO
+!      ALLOCATE(wq(ip1jmp1,llm+1,nqtot)) !MVals
 c$OMP END MASTER
 c$OMP BARRIER
@@ -1015,10 +1169,12 @@
        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
@@ -1028,21 +1184,73 @@
 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)
+         !print*,"margaux:vlsplt,PERE",iq
+         !print*,"margaux:vlsplt,FILS",iq2
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+         DO l=1,llm
+          DO ij=ijb,ije
+           !masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
+           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+           !MVals: veiller a ce qu'on n'ait pas de denominateur nul
+           masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),1e-16)
+           if (q(ij,l,iq).gt.1e-16) then
+             Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+           else 
+             Ratio(ij,l,iq2)=0.
+           endif
+           !wq(ij,l,iq2)=wq(ij,l,iq) ! correction bug le 15mai2015
+           w(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_p(Ratio,pente_max,masse,w,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
+
+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)=newmasse
-         ENDDO
-      ENDDO
-c$OMP END DO NOWAIT
-
+            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,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
 
       RETURN
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/vlspltgen_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/vlspltgen_p.F	(revision 2272)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/vlspltgen_p.F	(revision 2296)
@@ -26,5 +26,6 @@
       USE Write_Field_p
       USE VAMPIR
-      USE infotrac, ONLY : nqtot
+      USE infotrac, ONLY : nqtot,nqperes,nqdesc,nqfils,iqfils, 
+     &    ok_iso_verif
       USE comconst_mod, ONLY: cpp
       IMPLICIT NONE
@@ -53,5 +54,6 @@
       REAL,SAVE :: mu(ip1jmp1,llm)
       REAL,SAVE :: mv(ip1jm,llm)
-      REAL,SAVE :: mw(ip1jmp1,llm+1)
+      !REAL,SAVE :: mw(ip1jmp1,llm+1)
+      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: mw !MVals
       REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: zq
       REAL zzpbar, zzw
@@ -65,5 +67,5 @@
       REAL ptarg,pdelarg,foeew,zdelta
       REAL tempe(ip1jmp1)
-      INTEGER ijb,ije,iq
+      INTEGER ijb,ije,iq,iq2,ifils
       LOGICAL, SAVE :: firstcall=.TRUE.
 !$OMP THREADPRIVATE(firstcall)
@@ -91,4 +93,5 @@
 !$OMP MASTER
             ALLOCATE(zm(ip1jmp1,llm,nqtot))
+            ALLOCATE(mw(ip1jmp1,llm+1,nqtot)) !MVals
             ALLOCATE(zq(ip1jmp1,llm,nqtot))
 !$OMP END MASTER
@@ -155,17 +158,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.
-      ENDDO
-c$OMP END MASTER
+         mw(ij,llm+1,iq)=0.
+      ENDDO
+c$OMP END MASTER
+      ENDDO
 
 c      CALL SCOPY(ijp1llm,q,1,zq,1)
@@ -184,8 +191,14 @@
       ENDDO
 
+      ! verif temporaire
+      ijb=ij_begin
+      ije=ij_end
+      if (ok_iso_verif) then
+        call check_isotopes(zq,ijb,ije,'vlspltgen_loc 197')
+      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
         if(iadv(iq) == 0) then
 	
@@ -193,13 +206,12 @@
 	
 	else if (iadv(iq)==10) then
-
-#ifdef _ADV_HALO        
-	  call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
-     &	             ij_begin,ij_begin+2*iip1-1)
-          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
-     &               ij_end-2*iip1+1,ij_end)
+#ifdef _ADV_HALO     
+	  call vlx_p(zq,pente_max,zm,mu,
+     &	             ij_begin,ij_begin+2*iip1-1,iq)
+          call vlx_p(zq,pente_max,zm,mu,
+     &               ij_end-2*iip1+1,ij_end,iq)
 #else
-	  call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
-     &	             ij_begin,ij_end)
+	  call vlx_p(zq,pente_max,zm,mu,
+     &	             ij_begin,ij_end,iq)
 #endif
 
@@ -209,5 +221,10 @@
           call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
           call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
-
+! CRisi
+          do ifils=1,nqdesc(iq)
+            iq2=iqfils(ifils,iq)
+         call Register_Hallo(zq(1,1,iq2),ip1jmp1,llm,2,2,2,2,MyRequest1)
+         call Register_Hallo(zm(1,1,iq2),ip1jmp1,llm,1,1,1,1,MyRequest1)
+          enddo
 c$OMP MASTER
           call VTe(VTHallo)
@@ -232,5 +249,10 @@
           call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
           call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
-
+!CRisi
+          do ifils=1,nqdesc(iq)
+            iq2=iqfils(ifils,iq)
+         call Register_Hallo(zq(1,1,iq2),ip1jmp1,llm,2,2,2,2,MyRequest1)
+         call Register_Hallo(zm(1,1,iq2),ip1jmp1,llm,1,1,1,1,MyRequest1)
+          enddo
 c$OMP MASTER
           call VTe(VTHallo)
@@ -242,5 +264,5 @@
         endif
       
-      enddo
+      enddo !DO iq=1,nqperes
       
       
@@ -256,5 +278,16 @@
 c$OMP END MASTER       
 c$OMP BARRIER
-      do iq=1,nqtot
+
+!      ! 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 293')
+      endif !if (ok_iso_verif) then
+
+!      do iq=1,nqtot
+      do iq=1,nqperes !CRisi
 
         if(iadv(iq) == 0) then
@@ -265,6 +298,6 @@
 
 #ifdef _ADV_HALLO
-          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
-     &               ij_begin+2*iip1,ij_end-2*iip1)
+          call vlx_p(zq,pente_max,zm,mu,
+     &               ij_begin+2*iip1,ij_end-2*iip1,iq)
 #endif        
 	else if (iadv(iq)==14) then
@@ -295,6 +328,18 @@
 c$OMP END MASTER
 c$OMP BARRIER
- 
-      do iq=1,nqtot
+
+      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,nqtot
+      do iq=1,nqperes !CRisi
 
         if(iadv(iq) == 0) then
@@ -304,5 +349,5 @@
 	else if (iadv(iq)==10) then
         
-          call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
+          call vly_p(zq,pente_max,zm,mv,iq)
   
 	else if (iadv(iq)==14) then
@@ -318,6 +363,10 @@
        enddo
 
-
-      do iq=1,nqtot
+      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,nqtot
+      do iq=1,nqperes !CRisi
 
         if(iadv(iq) == 0) then 
@@ -329,11 +378,11 @@
 c$OMP BARRIER        
 #ifdef _ADV_HALLO
-          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
-     &               ij_begin,ij_begin+2*iip1-1)
-          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
-     &               ij_end-2*iip1+1,ij_end)
+          call vlz_p(zq,pente_max,zm,mw,
+     &               ij_begin,ij_begin+2*iip1-1,iq)
+          call vlz_p(zq,pente_max,zm,mw,
+     &               ij_end-2*iip1+1,ij_end,iq)
 #else
-          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
-     &               ij_begin,ij_end)
+          call vlz_p(zq,pente_max,zm,mw,
+     &               ij_begin,ij_end,iq)
 #endif
 c$OMP BARRIER
@@ -345,5 +394,10 @@
           call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest2)
           call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest2)
-
+          ! CRisi
+          do ifils=1,nqdesc(iq)
+            iq2=iqfils(ifils,iq)
+         call Register_Hallo(zq(1,1,iq2),ip1jmp1,llm,2,2,2,2,MyRequest2)
+         call Register_Hallo(zm(1,1,iq2),ip1jmp1,llm,1,1,1,1,MyRequest2)
+          enddo 
 c$OMP MASTER
           call VTe(VTHallo)
@@ -369,6 +423,11 @@
 c$OMP END MASTER	
 
-c$OMP BARRIER
-      do iq=1,nqtot
+      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,nqtot
+      do iq=1,nqperes !CRisi
 
         if(iadv(iq) == 0) then
@@ -380,6 +439,6 @@
 
 #ifdef _ADV_HALLO
-          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
-     &               ij_begin+2*iip1,ij_end-2*iip1)
+          call vlz_p(zq,pente_max,zm,mw,
+     &               ij_begin+2*iip1,ij_end-2*iip1,iq)
 #endif
 
@@ -408,6 +467,11 @@
 c$OMP BARRIER
 
-
-      do iq=1,nqtot
+      !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,nqtot
+      do iq=1,nqperes !CRisi
 
         if(iadv(iq) == 0) then
@@ -417,5 +481,5 @@
 	else if (iadv(iq)==10) then
         
-          call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
+          call vly_p(zq,pente_max,zm,mv,iq)
   
 	else if (iadv(iq)==14) then
@@ -429,7 +493,12 @@
         endif
        
-       enddo
-
-      do iq=1,nqtot
+       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,nqtot
+      do iq=1,nqperes !CRisi
 
         if(iadv(iq) == 0) then 
@@ -439,6 +508,6 @@
 	else if (iadv(iq)==10) then
         
-          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
-     &               ij_begin,ij_end)
+          call vlx_p(zq,pente_max,zm,mu,
+     &               ij_begin,ij_end,iq)
   
 	else if (iadv(iq)==14) then
@@ -453,6 +522,10 @@
         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
@@ -483,5 +556,8 @@
       ENDDO
         
-      
+      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
 
