Index: LMDZ5/branches/testing/libf/dyn3dmem/vlsplt_loc.F
===================================================================
--- LMDZ5/branches/testing/libf/dyn3dmem/vlsplt_loc.F	(revision 1910)
+++ LMDZ5/branches/testing/libf/dyn3dmem/vlsplt_loc.F	(revision 2298)
@@ -2,5 +2,5 @@
 ! $Id$
 !
-      SUBROUTINE vlx_loc(q,pente_max,masse,u_m,ijb_x,ije_x)
+      RECURSIVE 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 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
+           masse(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,masse,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)
+      RECURSIVE 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 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
+           masse(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,masse,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,5 +827,5 @@
 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
@@ -720,4 +833,24 @@
 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
       END
@@ -725,5 +858,5 @@
       
       
-      SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x)
+      RECURSIVE 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,37 +982,86 @@
 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
+           masse(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) ! 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_loc(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,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
 
       RETURN
