Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/advtrac.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/advtrac.f90	(revision 5118)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/advtrac.f90	(revision 5119)
@@ -17,4 +17,5 @@
   USE lmdz_libmath, ONLY: minmax
   USE lmdz_iniprint, ONLY: lunout, prt_level
+  USE lmdz_ssum_scopy, ONLY: scopy
 
   IMPLICIT NONE
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/caladvtrac.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/caladvtrac.F90	(revision 5118)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/caladvtrac.F90	(revision 5119)
@@ -11,4 +11,5 @@
   USE comconst_mod, ONLY: dtvr
   USE lmdz_filtreg, ONLY: filtreg
+  USE lmdz_ssum_scopy, ONLY: scopy
 
   IMPLICIT NONE
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/fluxstokenc.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/fluxstokenc.F90	(revision 5118)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/fluxstokenc.F90	(revision 5119)
@@ -7,4 +7,5 @@
   USE IOIPSL
   USE lmdz_iniprint, ONLY: lunout, prt_level
+  USE lmdz_ssum_scopy, ONLY: scopy
   !
   ! Auteur :  F. Hourdin
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/groupe.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/groupe.F90	(revision 5118)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/groupe.F90	(revision 5119)
@@ -4,4 +4,5 @@
 
   USE comconst_mod, ONLY: ngroup
+  USE lmdz_ssum_scopy, ONLY: scopy
 
   IMPLICIT NONE
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/integrd.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/integrd.F90	(revision 5118)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/integrd.F90	(revision 5119)
@@ -12,4 +12,5 @@
   USE temps_mod, ONLY: dt
   USE lmdz_iniprint, ONLY: lunout, prt_level
+  USE lmdz_ssum_scopy, ONLY: scopy
 
   IMPLICIT NONE
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/leapfrog.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/leapfrog.F90	(revision 5118)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/leapfrog.F90	(revision 5119)
@@ -27,4 +27,5 @@
   USE lmdz_description, ONLY: descript
   USE lmdz_iniprint, ONLY: lunout, prt_level
+  USE lmdz_ssum_scopy, ONLY: scopy
 
   IMPLICIT NONE
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/vlsplt.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/vlsplt.F90	(revision 5118)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/vlsplt.F90	(revision 5119)
@@ -3,6 +3,7 @@
 !
 
-SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq)
-  USE infotrac, ONLY: nqtot,tracers
+SUBROUTINE vlsplt(q, pente_max, masse, w, pbaru, pbarv, pdt, iq)
+  USE infotrac, ONLY: nqtot, tracers
+  USE lmdz_ssum_scopy, ONLY: scopy
   !
   ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
@@ -27,8 +28,8 @@
   !   Arguments:
   !   ----------
-  REAL :: masse(ip1jmp1,llm),pente_max
-  REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
-  REAL :: q(ip1jmp1,llm,nqtot)
-  REAL :: w(ip1jmp1,llm),pdt
+  REAL :: masse(ip1jmp1, llm), pente_max
+  REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
+  REAL :: q(ip1jmp1, llm, nqtot)
+  REAL :: w(ip1jmp1, llm), pdt
   INTEGER :: iq ! CRisi
   !
@@ -36,75 +37,74 @@
   !   ---------
   !
-  INTEGER :: ij,l
-  !
-  REAL :: zm(ip1jmp1,llm,nqtot)
-  REAL :: mu(ip1jmp1,llm)
-  REAL :: mv(ip1jm,llm)
-  REAL :: mw(ip1jmp1,llm+1)
-  REAL :: zq(ip1jmp1,llm,nqtot)
+  INTEGER :: ij, l
+  !
+  REAL :: zm(ip1jmp1, llm, nqtot)
+  REAL :: mu(ip1jmp1, llm)
+  REAL :: mv(ip1jm, llm)
+  REAL :: mw(ip1jmp1, llm + 1)
+  REAL :: zq(ip1jmp1, llm, nqtot)
   REAL :: zzpbar, zzw
-  INTEGER :: ifils,iq2 ! CRisi
-
-  REAL :: qmin,qmax
-  DATA qmin,qmax/0.,1.e33/
-
-    zzpbar = 0.5 * pdt
-    zzw    = pdt
-  DO l=1,llm
-    DO ij = iip2,ip1jm
-        mu(ij,l)=pbaru(ij,l) * zzpbar
-     ENDDO
-     DO ij=1,ip1jm
-        mv(ij,l)=pbarv(ij,l) * zzpbar
-     ENDDO
-     DO ij=1,ip1jmp1
-        mw(ij,l)=w(ij,l) * zzw
-     ENDDO
-  ENDDO
-
-  DO ij=1,ip1jmp1
-     mw(ij,llm+1)=0.
-  ENDDO
-
-  CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1)
-  CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
-
-  do ifils=1,tracers(iq)%nqDescen
-    iq2=tracers(iq)%iqDescen(ifils)
-    CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
-  enddo
-
-  CALL vlx(zq,pente_max,zm,mu,iq)
-  CALL vly(zq,pente_max,zm,mv,iq)
-  CALL vlz(zq,pente_max,zm,mw,iq)
-  CALL vly(zq,pente_max,zm,mv,iq)
-  CALL vlx(zq,pente_max,zm,mu,iq)
-
-  DO l=1,llm
-     DO ij=1,ip1jmp1
-       q(ij,l,iq)=zq(ij,l,iq)
-     ENDDO
-     DO ij=1,ip1jm+1,iip1
-        q(ij+iim,l,iq)=q(ij,l,iq)
-     ENDDO
+  INTEGER :: ifils, iq2 ! CRisi
+
+  REAL :: qmin, qmax
+  DATA qmin, qmax/0., 1.e33/
+
+  zzpbar = 0.5 * pdt
+  zzw = pdt
+  DO l = 1, llm
+    DO ij = iip2, ip1jm
+      mu(ij, l) = pbaru(ij, l) * zzpbar
+    ENDDO
+    DO ij = 1, ip1jm
+      mv(ij, l) = pbarv(ij, l) * zzpbar
+    ENDDO
+    DO ij = 1, ip1jmp1
+      mw(ij, l) = w(ij, l) * zzw
+    ENDDO
+  ENDDO
+
+  DO ij = 1, ip1jmp1
+    mw(ij, llm + 1) = 0.
+  ENDDO
+
+  CALL SCOPY(ijp1llm, q(1, 1, iq), 1, zq(1, 1, iq), 1)
+  CALL SCOPY(ijp1llm, masse, 1, zm(1, 1, iq), 1)
+
+  do ifils = 1, tracers(iq)%nqDescen
+    iq2 = tracers(iq)%iqDescen(ifils)
+    CALL SCOPY(ijp1llm, q(1, 1, iq2), 1, zq(1, 1, iq2), 1)
+  enddo
+
+  CALL vlx(zq, pente_max, zm, mu, iq)
+  CALL vly(zq, pente_max, zm, mv, iq)
+  CALL vlz(zq, pente_max, zm, mw, iq)
+  CALL vly(zq, pente_max, zm, mv, iq)
+  CALL vlx(zq, pente_max, zm, mu, iq)
+
+  DO l = 1, llm
+    DO ij = 1, ip1jmp1
+      q(ij, l, iq) = zq(ij, l, iq)
+    ENDDO
+    DO ij = 1, ip1jm + 1, iip1
+      q(ij + iim, l, iq) = q(ij, l, iq)
+    ENDDO
   ENDDO
   ! CRisi: aussi pour les fils
-  do ifils=1,tracers(iq)%nqDescen
-    iq2=tracers(iq)%iqDescen(ifils)
-    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, tracers(iq)%nqDescen
+    iq2 = tracers(iq)%iqDescen(ifils)
+    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
 
 END SUBROUTINE vlsplt
-RECURSIVE SUBROUTINE vlx(q,pente_max,masse,u_m,iq)
-  USE infotrac, ONLY: nqtot,tracers, & ! CRisi
-        min_qParent,min_qMass,min_ratio ! MVals et CRisi
+RECURSIVE SUBROUTINE vlx(q, pente_max, masse, u_m, iq)
+  USE infotrac, ONLY: nqtot, tracers, & ! CRisi
+          min_qParent, min_qMass, min_ratio ! MVals et CRisi
   USE lmdz_iniprint, ONLY: lunout, prt_level
 
@@ -126,7 +126,7 @@
   !   Arguments:
   !   ----------
-  REAL :: masse(ip1jmp1,llm,nqtot),pente_max
-  REAL :: u_m( ip1jmp1,llm )
-  REAL :: q(ip1jmp1,llm,nqtot)
+  REAL :: masse(ip1jmp1, llm, nqtot), pente_max
+  REAL :: u_m(ip1jmp1, llm)
+  REAL :: q(ip1jmp1, llm, nqtot)
   INTEGER :: iq ! CRisi
   !
@@ -134,16 +134,16 @@
   !   ---------
   !
-  INTEGER :: ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
-  INTEGER :: n0,iadvplus(ip1jmp1,llm),nl(llm)
-  !
-  REAL :: new_m,zu_m,zdum(ip1jmp1,llm)
-  REAL :: dxq(ip1jmp1,llm),dxqu(ip1jmp1)
+  INTEGER :: ij, l, j, i, iju, ijq, indu(ip1jmp1), niju
+  INTEGER :: n0, iadvplus(ip1jmp1, llm), nl(llm)
+  !
+  REAL :: new_m, zu_m, zdum(ip1jmp1, llm)
+  REAL :: dxq(ip1jmp1, llm), dxqu(ip1jmp1)
   REAL :: zz(ip1jmp1)
-  REAL :: adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
-  REAL :: u_mq(ip1jmp1,llm)
+  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
+  REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot)
+  INTEGER :: ifils, iq2 ! CRisi
 
   LOGICAL, SAVE :: first
@@ -152,80 +152,78 @@
   !   calcul de la pente a droite et a gauche de la maille
 
-
   IF (pente_max>-1.e-5) THEN
     ! IF (pente_max.gt.10) THEN
 
-  !   calcul des pentes avec limitation, Van Leer scheme I:
-  !   -----------------------------------------------------
-
-  !   calcul de la pente aux points u
-     DO l = 1, llm
-        DO ij=iip2,ip1jm-1
-           dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
-        ENDDO
-        DO ij=iip1+iip1,ip1jm,iip1
-           dxqu(ij)=dxqu(ij-iim)
-           ! sigu(ij)=sigu(ij-iim)
-        ENDDO
-
-        DO ij=iip2,ip1jm
-           adxqu(ij)=abs(dxqu(ij))
-        ENDDO
-
-  !   calcul de la pente maximum dans la maille en valeur absolue
-
-        DO ij=iip2+1,ip1jm
-           dxqmax(ij,l)=pente_max* &
-                 min(adxqu(ij-1),adxqu(ij))
-  ! limitation subtile
-  !    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
-
-
-        ENDDO
-
-        DO ij=iip1+iip1,ip1jm,iip1
-           dxqmax(ij-iim,l)=dxqmax(ij,l)
-        ENDDO
-
-        DO ij=iip2+1,ip1jm
-           IF(dxqu(ij-1)*dxqu(ij)>0) THEN
-              dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
-           ELSE
-  !   extremum local
-              dxq(ij,l)=0.
-           ENDIF
-           dxq(ij,l)=0.5*dxq(ij,l)
-           dxq(ij,l)= &
-                 sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
-        ENDDO
-
-     ENDDO ! l=1,llm
-  !print*,'Ok calcul des pentes'
+    !   calcul des pentes avec limitation, Van Leer scheme I:
+    !   -----------------------------------------------------
+
+    !   calcul de la pente aux points u
+    DO l = 1, llm
+      DO ij = iip2, ip1jm - 1
+        dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq)
+      ENDDO
+      DO ij = iip1 + iip1, ip1jm, iip1
+        dxqu(ij) = dxqu(ij - iim)
+        ! sigu(ij)=sigu(ij-iim)
+      ENDDO
+
+      DO ij = iip2, ip1jm
+        adxqu(ij) = abs(dxqu(ij))
+      ENDDO
+
+      !   calcul de la pente maximum dans la maille en valeur absolue
+
+      DO ij = iip2 + 1, ip1jm
+        dxqmax(ij, l) = pente_max * &
+                min(adxqu(ij - 1), adxqu(ij))
+        ! limitation subtile
+        !    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
+
+      ENDDO
+
+      DO ij = iip1 + iip1, ip1jm, iip1
+        dxqmax(ij - iim, l) = dxqmax(ij, l)
+      ENDDO
+
+      DO ij = iip2 + 1, ip1jm
+        IF(dxqu(ij - 1) * dxqu(ij)>0) THEN
+          dxq(ij, l) = dxqu(ij - 1) + dxqu(ij)
+        ELSE
+          !   extremum local
+          dxq(ij, l) = 0.
+        ENDIF
+        dxq(ij, l) = 0.5 * dxq(ij, l)
+        dxq(ij, l) = &
+                sign(min(abs(dxq(ij, l)), dxqmax(ij, l)), dxq(ij, l))
+      ENDDO
+
+    ENDDO ! l=1,llm
+    !print*,'Ok calcul des pentes'
 
   ELSE ! (pente_max.lt.-1.e-5)
 
-  !   Pentes produits:
-  !   ----------------
-
-     DO l = 1, llm
-        DO ij=iip2,ip1jm-1
-           dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
-        ENDDO
-        DO ij=iip1+iip1,ip1jm,iip1
-           dxqu(ij)=dxqu(ij-iim)
-        ENDDO
-
-        DO ij=iip2+1,ip1jm
-           zz(ij)=dxqu(ij-1)*dxqu(ij)
-           zz(ij)=zz(ij)+zz(ij)
-           IF(zz(ij)>0) THEN
-              dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
-           ELSE
-  !   extremum local
-              dxq(ij,l)=0.
-           ENDIF
-        ENDDO
-
-     ENDDO
+    !   Pentes produits:
+    !   ----------------
+
+    DO l = 1, llm
+      DO ij = iip2, ip1jm - 1
+        dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq)
+      ENDDO
+      DO ij = iip1 + iip1, ip1jm, iip1
+        dxqu(ij) = dxqu(ij - iim)
+      ENDDO
+
+      DO ij = iip2 + 1, ip1jm
+        zz(ij) = dxqu(ij - 1) * dxqu(ij)
+        zz(ij) = zz(ij) + zz(ij)
+        IF(zz(ij)>0) THEN
+          dxq(ij, l) = zz(ij) / (dxqu(ij - 1) + dxqu(ij))
+        ELSE
+          !   extremum local
+          dxq(ij, l) = 0.
+        ENDIF
+      ENDDO
+
+    ENDDO
 
   ENDIF ! (pente_max.lt.-1.e-5)
@@ -234,11 +232,11 @@
   !   -----------------------------
 
-  DO l=1,llm
-     DO ij=iip1+iip1,ip1jm,iip1
-        dxq(ij-iim,l)=dxq(ij,l)
-     ENDDO
-     DO ij=1,ip1jmp1
-        iadvplus(ij,l)=0
-     ENDDO
+  DO l = 1, llm
+    DO ij = iip1 + iip1, ip1jm, iip1
+      dxq(ij - iim, l) = dxq(ij, l)
+    ENDDO
+    DO ij = 1, ip1jmp1
+      iadvplus(ij, l) = 0
+    ENDDO
 
   ENDDO
@@ -247,31 +245,31 @@
   !   on cumule le flux correspondant a toutes les mailles dont la masse
   !   au travers de la paroi pENDant le pas de temps.
-  DO l=1,llm
-   DO ij=iip2,ip1jm-1
-      IF (u_m(ij,l)>0.) THEN
-         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))
+  DO l = 1, llm
+    DO ij = iip2, ip1jm - 1
+      IF (u_m(ij, l)>0.) THEN
+        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,iq)
-         u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq) &
-               -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
+    ENDDO
   ENDDO
 
   !   detection des points ou on advecte plus que la masse de la
   !   maille
-  DO l=1,llm
-     DO ij=iip2,ip1jm-1
-        IF(zdum(ij,l)<0) THEN
-           iadvplus(ij,l)=1
-           u_mq(ij,l)=0.
-        ENDIF
-     ENDDO
-  ENDDO
-  DO l=1,llm
-   DO ij=iip1+iip1,ip1jm,iip1
-      iadvplus(ij,l)=iadvplus(ij-iim,l)
-   ENDDO
+  DO l = 1, llm
+    DO ij = iip2, ip1jm - 1
+      IF(zdum(ij, l)<0) THEN
+        iadvplus(ij, l) = 1
+        u_mq(ij, l) = 0.
+      ENDIF
+    ENDDO
+  ENDDO
+  DO l = 1, llm
+    DO ij = iip1 + iip1, ip1jm, iip1
+      iadvplus(ij, l) = iadvplus(ij - iim, l)
+    ENDDO
   ENDDO
 
@@ -283,69 +281,69 @@
   !  calcul du nombre de maille sur lequel on advecte plus que la maille.
 
-  n0=0
-  DO l=1,llm
-     nl(l)=0
-     DO ij=iip2,ip1jm
-        nl(l)=nl(l)+iadvplus(ij,l)
-     ENDDO
-     n0=n0+nl(l)
+  n0 = 0
+  DO l = 1, llm
+    nl(l) = 0
+    DO ij = iip2, ip1jm
+      nl(l) = nl(l) + iadvplus(ij, l)
+    ENDDO
+    n0 = n0 + nl(l)
   ENDDO
 
   IF(n0>0) THEN
-  IF (prt_level > 2) PRINT *, &
-        'Nombre de points pour lesquels on advect plus que le' &
-        ,'contenu de la maille : ',n0
-
-     DO l=1,llm
-        IF(nl(l)>0) THEN
-           iju=0
-  !   indicage des mailles concernees par le traitement special
-           DO ij=iip2,ip1jm
-              IF(iadvplus(ij,l)==1.AND.mod(ij,iip1)/=0) THEN
-                 iju=iju+1
-                 indu(iju)=ij
-              ENDIF
-           ENDDO
-           niju=iju
-
-  !  traitement des mailles
-           DO iju=1,niju
-              ij=indu(iju)
-              j=(ij-1)/iip1+1
-              zu_m=u_m(ij,l)
-              u_mq(ij,l)=0.
-              IF(zu_m>0.) THEN
-                 ijq=ij
-                 i=ijq-(j-1)*iip1
-  !   accumulation pour les mailles completements advectees
-                 do while(zu_m>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
-  !   ajout de la maille non completement advectee
-                 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
-  !   accumulation pour les mailles completements advectees
-                 do while(-zu_m>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
-  !   ajout de la maille non completement advectee
-                 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
-        ENDIF
-     ENDDO
+    IF (prt_level > 2) PRINT *, &
+            'Nombre de points pour lesquels on advect plus que le' &
+            , 'contenu de la maille : ', n0
+
+    DO l = 1, llm
+      IF(nl(l)>0) THEN
+        iju = 0
+        !   indicage des mailles concernees par le traitement special
+        DO ij = iip2, ip1jm
+          IF(iadvplus(ij, l)==1.AND.mod(ij, iip1)/=0) THEN
+            iju = iju + 1
+            indu(iju) = ij
+          ENDIF
+        ENDDO
+        niju = iju
+
+        !  traitement des mailles
+        DO iju = 1, niju
+          ij = indu(iju)
+          j = (ij - 1) / iip1 + 1
+          zu_m = u_m(ij, l)
+          u_mq(ij, l) = 0.
+          IF(zu_m>0.) THEN
+            ijq = ij
+            i = ijq - (j - 1) * iip1
+            !   accumulation pour les mailles completements advectees
+            do while(zu_m>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
+            !   ajout de la maille non completement advectee
+            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
+            !   accumulation pour les mailles completements advectees
+            do while(-zu_m>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
+            !   ajout de la maille non completement advectee
+            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
+      ENDIF
+    ENDDO
   ENDIF  ! n0.gt.0
 
@@ -353,7 +351,7 @@
   !   bouclage en latitude
   !print*,'cvant bouclage en latitude'
-  DO l=1,llm
-    DO ij=iip1+iip1,ip1jm,iip1
-       u_mq(ij,l)=u_mq(ij-iim,l)
+  DO l = 1, llm
+    DO ij = iip1 + iip1, ip1jm, iip1
+      u_mq(ij, l) = u_mq(ij - iim, l)
     ENDDO
   ENDDO
@@ -363,24 +361,24 @@
   !WRITE(*,*) 'vlsplt 326: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen
 
-  do ifils=1,tracers(iq)%nqDescen
-    iq2=tracers(iq)%iqDescen(ifils)
-    DO l=1,llm
-      DO ij=iip2,ip1jm
+  do ifils = 1, tracers(iq)%nqDescen
+    iq2 = tracers(iq)%iqDescen(ifils)
+    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)
+        !           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
         !Mvals: veiller a ce qu'on n'ait pas de denominateur nul
-        masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
-        IF (q(ij,l,iq)>min_qParent) THEN
-          Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+        masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
+        IF (q(ij, l, iq)>min_qParent) THEN
+          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
         else
-          Ratio(ij,l,iq2)=min_ratio
+          Ratio(ij, l, iq2) = min_ratio
         endif
       enddo
     enddo
   enddo
-  do ifils=1,tracers(iq)%nqChildren
-    iq2=tracers(iq)%iqDescen(ifils)
-    CALL vlx(Ratio,pente_max,masseq,u_mq,iq2)
+  do ifils = 1, tracers(iq)%nqChildren
+    iq2 = tracers(iq)%iqDescen(ifils)
+    CALL vlx(Ratio, pente_max, masseq, u_mq, iq2)
   enddo
   ! end CRisi
@@ -389,17 +387,17 @@
   !   calcul des tENDances
 
-  DO l=1,llm
-     DO ij=iip2+1,ip1jm
-        !MVals: veiller a ce qu'on ait pas de denominateur nul
-        new_m=max(masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l),min_qMass)
-        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
-     DO ij=iip1+iip1,ip1jm,iip1
-        q(ij-iim,l,iq)=q(ij,l,iq)
-        masse(ij-iim,l,iq)=masse(ij,l,iq)
-     ENDDO
+  DO l = 1, llm
+    DO ij = iip2 + 1, ip1jm
+      !MVals: veiller a ce qu'on ait pas de denominateur nul
+      new_m = max(masse(ij, l, iq) + u_m(ij - 1, l) - u_m(ij, l), min_qMass)
+      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
+    DO ij = iip1 + iip1, ip1jm, iip1
+      q(ij - iim, l, iq) = q(ij, l, iq)
+      masse(ij - iim, l, iq) = masse(ij, l, iq)
+    ENDDO
   ENDDO
 
@@ -407,21 +405,20 @@
   ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
   ! puis on boucle en longitude
-  do ifils=1,tracers(iq)%nqDescen
-    iq2=tracers(iq)%iqDescen(ifils)
-    DO l=1,llm
-      DO ij=iip2+1,ip1jm
-        q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
+  do ifils = 1, tracers(iq)%nqDescen
+    iq2 = tracers(iq)%iqDescen(ifils)
+    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)
+      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
 
-
 END SUBROUTINE vlx
-RECURSIVE SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq)
-  USE infotrac, ONLY: nqtot,tracers, & ! CRisi
-        min_qParent,min_qMass,min_ratio ! MVals et CRisi
+RECURSIVE SUBROUTINE vly(q, pente_max, masse, masse_adv_v, iq)
+  USE infotrac, ONLY: nqtot, tracers, & ! CRisi
+          min_qParent, min_qMass, min_ratio ! MVals et CRisi
   !
   ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
@@ -445,7 +442,7 @@
   !   Arguments:
   !   ----------
-  REAL :: masse(ip1jmp1,llm,nqtot),pente_max
-  REAL :: masse_adv_v( ip1jm,llm)
-  REAL :: q(ip1jmp1,llm,nqtot)
+  REAL :: masse(ip1jmp1, llm, nqtot), pente_max
+  REAL :: masse_adv_v(ip1jm, llm)
+  REAL :: q(ip1jmp1, llm, nqtot)
   INTEGER :: iq ! CRisi
   !
@@ -453,23 +450,23 @@
   !   ---------
   !
-  INTEGER :: i,ij,l
-  !
-  REAL :: airej2,airejjm,airescb(iim),airesch(iim)
-  REAL :: dyq(ip1jmp1,llm),dyqv(ip1jm)
-  REAL :: adyqv(ip1jm),dyqmax(ip1jmp1)
-  REAL :: qbyv(ip1jm,llm)
-
-  REAL :: qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
+  INTEGER :: i, ij, l
+  !
+  REAL :: airej2, airejjm, airescb(iim), airesch(iim)
+  REAL :: dyq(ip1jmp1, llm), dyqv(ip1jm)
+  REAL :: adyqv(ip1jm), dyqmax(ip1jmp1)
+  REAL :: qbyv(ip1jm, llm)
+
+  REAL :: qpns, qpsn, dyn1, dys1, dyn2, dys2, newmasse, fn, fs
   LOGICAL, SAVE :: first
 
-  REAL :: convpn,convps,convmpn,convmps
-  REAL :: massepn,masseps,qpn,qps
-  REAL :: sinlon(iip1),sinlondlon(iip1)
-  REAL :: coslon(iip1),coslondlon(iip1)
-  SAVE sinlon,coslon,sinlondlon,coslondlon
-  SAVE airej2,airejjm
-
-  REAL :: masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
-  INTEGER :: ifils,iq2 ! CRisi
+  REAL :: convpn, convps, convmpn, convmps
+  REAL :: massepn, masseps, qpn, qps
+  REAL :: sinlon(iip1), sinlondlon(iip1)
+  REAL :: coslon(iip1), coslondlon(iip1)
+  SAVE sinlon, coslon, sinlondlon, coslondlon
+  SAVE airej2, airejjm
+
+  REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot) ! CRisi
+  INTEGER :: ifils, iq2 ! CRisi
 
   !
@@ -482,18 +479,18 @@
 
   IF(first) THEN
-     PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
-     first=.FALSE.
-     do i=2,iip1
-        coslon(i)=cos(rlonv(i))
-        sinlon(i)=sin(rlonv(i))
-        coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
-        sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
-     ENDDO
-     coslon(1)=coslon(iip1)
-     coslondlon(1)=coslondlon(iip1)
-     sinlon(1)=sinlon(iip1)
-     sinlondlon(1)=sinlondlon(iip1)
-     airej2 = SSUM( iim, aire(iip2), 1 )
-     airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
+    PRINT*, 'Shema  Amont nouveau  appele dans  Vanleer   '
+    first = .FALSE.
+    do i = 2, iip1
+      coslon(i) = cos(rlonv(i))
+      sinlon(i) = sin(rlonv(i))
+      coslondlon(i) = coslon(i) * (rlonu(i) - rlonu(i - 1)) / pi
+      sinlondlon(i) = sinlon(i) * (rlonu(i) - rlonu(i - 1)) / pi
+    ENDDO
+    coslon(1) = coslon(iip1)
+    coslondlon(1) = coslondlon(iip1)
+    sinlon(1) = sinlon(iip1)
+    sinlondlon(1) = sinlondlon(iip1)
+    airej2 = SSUM(iim, aire(iip2), 1)
+    airejjm = SSUM(iim, aire(ip1jm - iim), 1)
   ENDIF
 
@@ -502,194 +499,194 @@
 
   DO l = 1, llm
-  !
-  !   --------------------------------
-  !  CALCUL EN LATITUDE
-  !   --------------------------------
-
-  !   On commence par calculer la valeur du traceur moyenne sur le premier cercle
-  !   de latitude autour du pole (qpns pour le pole nord et qpsn pour
-  !    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
-
-  DO i = 1, iim
-  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
-  qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
-
-  !   calcul des pentes aux points v
-
-  DO ij=1,ip1jm
-     dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
-     adyqv(ij)=abs(dyqv(ij))
-  ENDDO
-
-  !   calcul des pentes aux points scalaires
-
-  DO ij=iip2,ip1jm
-     dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
-     dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
-     dyqmax(ij)=pente_max*dyqmax(ij)
-  ENDDO
-
-  !   calcul des pentes aux poles
-
-  DO ij=1,iip1
-     dyq(ij,l)=qpns-q(ij+iip1,l,iq)
-     dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
-  ENDDO
-
-  !   filtrage de la derivee
-  dyn1=0.
-  dys1=0.
-  dyn2=0.
-  dys2=0.
-  DO ij=1,iim
-     dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
-     dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
-     dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
-     dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
-  ENDDO
-  DO ij=1,iip1
-     dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
-     dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
-  ENDDO
-
-  !   calcul des pentes limites aux poles
-
-  goto 8888
-  fn=1.
-  fs=1.
-  DO ij=1,iim
-     IF(pente_max*adyqv(ij)<abs(dyq(ij,l))) THEN
-        fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
-     ENDIF
-  IF(pente_max*adyqv(ij+ip1jm-iip1)<abs(dyq(ij+ip1jm,l))) THEN
-     fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
-     ENDIF
-  ENDDO
-  DO ij=1,iip1
-     dyq(ij,l)=fn*dyq(ij,l)
-     dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
-  ENDDO
-8888   continue
-  DO ij=1,iip1
-     dyq(ij,l)=0.
-     dyq(ip1jm+ij,l)=0.
-  ENDDO
-
-  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-  !  En memoire de dIFferents tests sur la
-  !  limitation des pentes aux poles.
-  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-  ! PRINT*,dyq(1)
-  ! PRINT*,dyqv(iip1+1)
-  ! appn=abs(dyq(1)/dyqv(iip1+1))
-  ! PRINT*,dyq(ip1jm+1)
-  ! PRINT*,dyqv(ip1jm-iip1+1)
-  ! apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
-  ! DO ij=2,iim
-  !    appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
-  !    apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
-  ! ENDDO
-  ! appn=min(pente_max/appn,1.)
-  ! apps=min(pente_max/apps,1.)
-  !
-  !
-  !   cas ou on a un extremum au pole
-  !
-  ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
-  !    &   appn=0.
-  ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
-  !    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
-  !    &   apps=0.
-  !
-  !   limitation des pentes aux poles
-  ! DO ij=1,iip1
-  !    dyq(ij)=appn*dyq(ij)
-  !    dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
-  ! ENDDO
-  !
-  !   test
-  !  DO ij=1,iip1
-  !     dyq(iip1+ij)=0.
-  !     dyq(ip1jm+ij-iip1)=0.
-  !  ENDDO
-  !  DO ij=1,ip1jmp1
-  !     dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
-  !  ENDDO
-  !
-  ! changement 10 07 96
-  ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
-  !    &   THEN
-  !    DO ij=1,iip1
-  !       dyqmax(ij)=0.
-  !    ENDDO
-  ! ELSE
-  !    DO ij=1,iip1
-  !       dyqmax(ij)=pente_max*abs(dyqv(ij))
-  !    ENDDO
-  ! ENDIF
-  !
-  ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
-  !    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
-  !    &THEN
-  !    DO ij=ip1jm+1,ip1jmp1
-  !       dyqmax(ij)=0.
-  !    ENDDO
-  ! ELSE
-  !    DO ij=ip1jm+1,ip1jmp1
-  !       dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
-  !    ENDDO
-  ! ENDIF
-  !   fin changement 10 07 96
-  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-
-  !   calcul des pentes limitees
-
-  DO ij=iip2,ip1jm
-     IF(dyqv(ij)*dyqv(ij-iip1)>0.) THEN
-        dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
-     ELSE
-        dyq(ij,l)=0.
-     ENDIF
-  ENDDO
+    !
+    !   --------------------------------
+    !  CALCUL EN LATITUDE
+    !   --------------------------------
+
+    !   On commence par calculer la valeur du traceur moyenne sur le premier cercle
+    !   de latitude autour du pole (qpns pour le pole nord et qpsn pour
+    !    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
+
+    DO i = 1, iim
+      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
+    qpsn = SSUM(iim, airesch, 1) / airejjm
+
+    !   calcul des pentes aux points v
+
+    DO ij = 1, ip1jm
+      dyqv(ij) = q(ij, l, iq) - q(ij + iip1, l, iq)
+      adyqv(ij) = abs(dyqv(ij))
+    ENDDO
+
+    !   calcul des pentes aux points scalaires
+
+    DO ij = iip2, ip1jm
+      dyq(ij, l) = .5 * (dyqv(ij - iip1) + dyqv(ij))
+      dyqmax(ij) = min(adyqv(ij - iip1), adyqv(ij))
+      dyqmax(ij) = pente_max * dyqmax(ij)
+    ENDDO
+
+    !   calcul des pentes aux poles
+
+    DO ij = 1, iip1
+      dyq(ij, l) = qpns - q(ij + iip1, l, iq)
+      dyq(ip1jm + ij, l) = q(ip1jm + ij - iip1, l, iq) - qpsn
+    ENDDO
+
+    !   filtrage de la derivee
+    dyn1 = 0.
+    dys1 = 0.
+    dyn2 = 0.
+    dys2 = 0.
+    DO ij = 1, iim
+      dyn1 = dyn1 + sinlondlon(ij) * dyq(ij, l)
+      dys1 = dys1 + sinlondlon(ij) * dyq(ip1jm + ij, l)
+      dyn2 = dyn2 + coslondlon(ij) * dyq(ij, l)
+      dys2 = dys2 + coslondlon(ij) * dyq(ip1jm + ij, l)
+    ENDDO
+    DO ij = 1, iip1
+      dyq(ij, l) = dyn1 * sinlon(ij) + dyn2 * coslon(ij)
+      dyq(ip1jm + ij, l) = dys1 * sinlon(ij) + dys2 * coslon(ij)
+    ENDDO
+
+    !   calcul des pentes limites aux poles
+
+    goto 8888
+    fn = 1.
+    fs = 1.
+    DO ij = 1, iim
+      IF(pente_max * adyqv(ij)<abs(dyq(ij, l))) THEN
+        fn = min(pente_max * adyqv(ij) / abs(dyq(ij, l)), fn)
+      ENDIF
+      IF(pente_max * adyqv(ij + ip1jm - iip1)<abs(dyq(ij + ip1jm, l))) THEN
+        fs = min(pente_max * adyqv(ij + ip1jm - iip1) / abs(dyq(ij + ip1jm, l)), fs)
+      ENDIF
+    ENDDO
+    DO ij = 1, iip1
+      dyq(ij, l) = fn * dyq(ij, l)
+      dyq(ip1jm + ij, l) = fs * dyq(ip1jm + ij, l)
+    ENDDO
+    8888   continue
+    DO ij = 1, iip1
+      dyq(ij, l) = 0.
+      dyq(ip1jm + ij, l) = 0.
+    ENDDO
+
+    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+    !  En memoire de dIFferents tests sur la
+    !  limitation des pentes aux poles.
+    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+    ! PRINT*,dyq(1)
+    ! PRINT*,dyqv(iip1+1)
+    ! appn=abs(dyq(1)/dyqv(iip1+1))
+    ! PRINT*,dyq(ip1jm+1)
+    ! PRINT*,dyqv(ip1jm-iip1+1)
+    ! apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
+    ! DO ij=2,iim
+    !    appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
+    !    apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
+    ! ENDDO
+    ! appn=min(pente_max/appn,1.)
+    ! apps=min(pente_max/apps,1.)
+    !
+    !
+    !   cas ou on a un extremum au pole
+    !
+    ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
+    !    &   appn=0.
+    ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
+    !    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
+    !    &   apps=0.
+    !
+    !   limitation des pentes aux poles
+    ! DO ij=1,iip1
+    !    dyq(ij)=appn*dyq(ij)
+    !    dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
+    ! ENDDO
+    !
+    !   test
+    !  DO ij=1,iip1
+    !     dyq(iip1+ij)=0.
+    !     dyq(ip1jm+ij-iip1)=0.
+    !  ENDDO
+    !  DO ij=1,ip1jmp1
+    !     dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
+    !  ENDDO
+    !
+    ! changement 10 07 96
+    ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
+    !    &   THEN
+    !    DO ij=1,iip1
+    !       dyqmax(ij)=0.
+    !    ENDDO
+    ! ELSE
+    !    DO ij=1,iip1
+    !       dyqmax(ij)=pente_max*abs(dyqv(ij))
+    !    ENDDO
+    ! ENDIF
+    !
+    ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
+    !    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
+    !    &THEN
+    !    DO ij=ip1jm+1,ip1jmp1
+    !       dyqmax(ij)=0.
+    !    ENDDO
+    ! ELSE
+    !    DO ij=ip1jm+1,ip1jmp1
+    !       dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
+    !    ENDDO
+    ! ENDIF
+    !   fin changement 10 07 96
+    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+
+    !   calcul des pentes limitees
+
+    DO ij = iip2, ip1jm
+      IF(dyqv(ij) * dyqv(ij - iip1)>0.) THEN
+        dyq(ij, l) = sign(min(abs(dyq(ij, l)), dyqmax(ij)), dyq(ij, l))
+      ELSE
+        dyq(ij, l) = 0.
+      ENDIF
+    ENDDO
 
   ENDDO
 
   ! !WRITE(*,*) 'vly 756'
-  DO l=1,llm
-   DO ij=1,ip1jm
-      IF(masse_adv_v(ij,l)>0) THEN
-          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))
+  DO l = 1, llm
+    DO ij = 1, ip1jm
+      IF(masse_adv_v(ij, l)>0) THEN
+        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,iq)-dyq(ij,l)* &
-                0.5*(1.+masse_adv_v(ij,l) &
-                /masse(ij,l,iq))
+        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)
-   ENDDO
+      qbyv(ij, l) = masse_adv_v(ij, l) * qbyv(ij, l)
+    ENDDO
   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,nqDesc(iq)=',iq,tracers(iq)%nqDescen
-
-  do ifils=1,tracers(iq)%nqDescen
-    iq2=tracers(iq)%iqDescen(ifils)
-    DO l=1,llm
-      DO ij=1,ip1jmp1
+  ! WRITE(*,*) 'vly 689: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen
+
+  do ifils = 1, tracers(iq)%nqDescen
+    iq2 = tracers(iq)%iqDescen(ifils)
+    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)
+        !           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
         ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul
-        masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
-        IF (q(ij,l,iq)>min_qParent) THEN
-          Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+        masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
+        IF (q(ij, l, iq)>min_qParent) THEN
+          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
         else
-          Ratio(ij,l,iq2)=min_ratio
+          Ratio(ij, l, iq2) = min_ratio
         endif
       enddo
@@ -697,47 +694,47 @@
   enddo
 
-  do ifils=1,tracers(iq)%nqDescen
-    iq2=tracers(iq)%iqDescen(ifils)
-    CALL vly(Ratio,pente_max,masseq,qbyv,iq2)
-  enddo
-
-  DO l=1,llm
-     DO ij=iip2,ip1jm
-        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
-     convpn=SSUM(iim,qbyv(1,l),1)
-     convmpn=ssum(iim,masse_adv_v(1,l),1)
-     massepn=ssum(iim,masse(1,l,iq),1)
-     qpn=0.
-     do ij=1,iim
-        qpn=qpn+masse(ij,l,iq)*q(ij,l,iq)
-     enddo
-     qpn=(qpn+convpn)/(massepn+convmpn)
-     do ij=1,iip1
-        q(ij,l,iq)=qpn
-     enddo
-     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,iq),1)
-     qps=0.
-     do ij = ip1jm+1,ip1jmp1-1
-        qps=qps+masse(ij,l,iq)*q(ij,l,iq)
-     enddo
-     qps=(qps+convps)/(masseps+convmps)
-     do ij=ip1jm+1,ip1jmp1
-        q(ij,l,iq)=qps
-     enddo
+  do ifils = 1, tracers(iq)%nqDescen
+    iq2 = tracers(iq)%iqDescen(ifils)
+    CALL vly(Ratio, pente_max, masseq, qbyv, iq2)
+  enddo
+
+  DO l = 1, llm
+    DO ij = iip2, ip1jm
+      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
+    convpn = SSUM(iim, qbyv(1, l), 1)
+    convmpn = ssum(iim, masse_adv_v(1, l), 1)
+    massepn = ssum(iim, masse(1, l, iq), 1)
+    qpn = 0.
+    do ij = 1, iim
+      qpn = qpn + masse(ij, l, iq) * q(ij, l, iq)
+    enddo
+    qpn = (qpn + convpn) / (massepn + convmpn)
+    do ij = 1, iip1
+      q(ij, l, iq) = qpn
+    enddo
+    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, iq), 1)
+    qps = 0.
+    do ij = ip1jm + 1, ip1jmp1 - 1
+      qps = qps + masse(ij, l, iq) * q(ij, l, iq)
+    enddo
+    qps = (qps + convps) / (masseps + convmps)
+    do ij = ip1jm + 1, ip1jmp1
+      q(ij, l, iq) = qps
+    enddo
   ENDDO
 
   ! retablir les fils en rapport de melange par rapport a l'air:
-  do ifils=1,tracers(iq)%nqDescen
-    iq2=tracers(iq)%iqDescen(ifils)
-    DO l=1,llm
-      DO ij=1,ip1jmp1
-        q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
+  do ifils = 1, tracers(iq)%nqDescen
+    iq2 = tracers(iq)%iqDescen(ifils)
+    DO l = 1, llm
+      DO ij = 1, ip1jmp1
+        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
       enddo
     enddo
@@ -746,9 +743,8 @@
   ! !WRITE(*,*) 'vly 853: sortie'
 
-
 END SUBROUTINE vly
-RECURSIVE SUBROUTINE vlz(q,pente_max,masse,w,iq)
-  USE infotrac, ONLY: nqtot,tracers, & ! CRisi
-        min_qParent,min_qMass,min_ratio ! MVals et CRisi
+RECURSIVE SUBROUTINE vlz(q, pente_max, masse, w, iq)
+  USE infotrac, ONLY: nqtot, tracers, & ! CRisi
+          min_qParent, min_qMass, min_ratio ! MVals et CRisi
   !
   ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
@@ -768,7 +764,7 @@
   !   Arguments:
   !   ----------
-  REAL :: masse(ip1jmp1,llm,nqtot),pente_max
-  REAL :: q(ip1jmp1,llm,nqtot)
-  REAL :: w(ip1jmp1,llm+1)
+  REAL :: masse(ip1jmp1, llm, nqtot), pente_max
+  REAL :: q(ip1jmp1, llm, nqtot)
+  REAL :: w(ip1jmp1, llm + 1)
   INTEGER :: iq
   !
@@ -776,13 +772,13 @@
   !   ---------
   !
-  INTEGER :: ij,l
-  !
-  REAL :: wq(ip1jmp1,llm+1),newmasse
-
-  REAL :: dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm),dzqmax
+  INTEGER :: ij, l
+  !
+  REAL :: wq(ip1jmp1, llm + 1), newmasse
+
+  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
+  REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot) ! CRisi
+  INTEGER :: ifils, iq2 ! CRisi
 
 #ifdef BIDON
@@ -795,27 +791,27 @@
   !    On oriente tout dans le sens de la pression c'est a dire dans le
   !    sens de W
-  DO l=2,llm
-     DO ij=1,ip1jmp1
-        dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq)
-        adzqw(ij,l)=abs(dzqw(ij,l))
-     ENDDO
-  ENDDO
-
-  DO l=2,llm-1
-     DO ij=1,ip1jmp1
-        IF(dzqw(ij,l)*dzqw(ij,l+1)>0.) THEN
-            dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
-        ELSE
-            dzq(ij,l)=0.
-        ENDIF
-        dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
-        dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
-     ENDDO
+  DO l = 2, llm
+    DO ij = 1, ip1jmp1
+      dzqw(ij, l) = q(ij, l - 1, iq) - q(ij, l, iq)
+      adzqw(ij, l) = abs(dzqw(ij, l))
+    ENDDO
+  ENDDO
+
+  DO l = 2, llm - 1
+    DO ij = 1, ip1jmp1
+      IF(dzqw(ij, l) * dzqw(ij, l + 1)>0.) THEN
+        dzq(ij, l) = 0.5 * (dzqw(ij, l) + dzqw(ij, l + 1))
+      ELSE
+        dzq(ij, l) = 0.
+      ENDIF
+      dzqmax = pente_max * min(adzqw(ij, l), adzqw(ij, l + 1))
+      dzq(ij, l) = sign(min(abs(dzq(ij, l)), dzqmax), dzq(ij, l))
+    ENDDO
   ENDDO
 
   ! !WRITE(*,*) 'vlz 954'
-  DO ij=1,ip1jmp1
-     dzq(ij,1)=0.
-     dzq(ij,llm)=0.
+  DO ij = 1, ip1jmp1
+    dzq(ij, 1) = 0.
+    dzq(ij, llm) = 0.
   ENDDO
 
@@ -826,37 +822,37 @@
   ! calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
 
-   DO l = 1,llm-1
-     do  ij = 1,ip1jmp1
-      IF(w(ij,l+1)>0.) THEN
-         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))
+  DO l = 1, llm - 1
+    do  ij = 1, ip1jmp1
+      IF(w(ij, l + 1)>0.) THEN
+        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,iq)
-         wq(ij,l+1)=w(ij,l+1)*(q(ij,l,iq)-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
-   ENDDO
-
-   DO ij=1,ip1jmp1
-      wq(ij,llm+1)=0.
-      wq(ij,1)=0.
-   ENDDO
+    ENDDO
+  ENDDO
+
+  DO ij = 1, ip1jmp1
+    wq(ij, llm + 1) = 0.
+    wq(ij, 1) = 0.
+  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,nqChildren(iq)=',iq,nqChildren(iq)
-  do ifils=1,tracers(iq)%nqDescen
-    iq2=tracers(iq)%iqDescen(ifils)
-    DO l=1,llm
-      DO ij=1,ip1jmp1
+  do ifils = 1, tracers(iq)%nqDescen
+    iq2 = tracers(iq)%iqDescen(ifils)
+    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)
+        !           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
         ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul
-        masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
-        IF (q(ij,l,iq)>min_qParent) THEN
-          Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
+        masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
+        IF (q(ij, l, iq)>min_qParent) THEN
+          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
         else
-          Ratio(ij,l,iq2)=min_ratio
+          Ratio(ij, l, iq2) = min_ratio
         endif
       enddo
@@ -864,41 +860,40 @@
   enddo
 
-  do ifils=1,tracers(iq)%nqChildren
-    iq2=tracers(iq)%iqDescen(ifils)
-    CALL vlz(Ratio,pente_max,masseq,wq,iq2)
+  do ifils = 1, tracers(iq)%nqChildren
+    iq2 = tracers(iq)%iqDescen(ifils)
+    CALL vlz(Ratio, pente_max, masseq, wq, iq2)
   enddo
   ! end CRisi
 
-  DO l=1,llm
-     DO ij=1,ip1jmp1
-        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,iq)=newmasse
-     ENDDO
+  DO l = 1, llm
+    DO ij = 1, ip1jmp1
+      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, iq) = newmasse
+    ENDDO
   ENDDO
 
   ! retablir les fils en rapport de melange par rapport a l'air:
-  do ifils=1,tracers(iq)%nqDescen
-    iq2=tracers(iq)%iqDescen(ifils)
-    DO l=1,llm
-      DO ij=1,ip1jmp1
-        q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
+  do ifils = 1, tracers(iq)%nqDescen
+    iq2 = tracers(iq)%iqDescen(ifils)
+    DO l = 1, llm
+      DO ij = 1, ip1jmp1
+        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
       enddo
     enddo
   enddo
 
-
 END SUBROUTINE vlz
 
-SUBROUTINE minmaxq(zq,qmin,qmax,comment)
+SUBROUTINE minmaxq(zq, qmin, qmax, comment)
 
   INCLUDE "dimensions.h"
   INCLUDE "paramet.h"
 
-  CHARACTER(LEN=20) :: comment
-  REAL :: qmin,qmax
-  REAL :: zq(ip1jmp1,llm)
-  REAL :: zzq(iip1,jjp1,llm)
+  CHARACTER(LEN = 20) :: comment
+  REAL :: qmin, qmax
+  REAL :: zq(ip1jmp1, llm)
+  REAL :: zzq(iip1, jjp1, llm)
 
 END SUBROUTINE  minmaxq
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/vlspltqs.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/vlspltqs.F90	(revision 5118)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/vlspltqs.F90	(revision 5119)
@@ -26,4 +26,5 @@
   USE comconst_mod, ONLY: cpp
   USE logic_mod, ONLY: adv_qsat_liq
+  USE lmdz_ssum_scopy, ONLY: scopy
   IMPLICIT NONE
   !
@@ -172,5 +173,4 @@
   enddo
   !WRITE(*,*) 'vlspltqs 183: fin de la routine'
-
 
 END SUBROUTINE vlspltqs
@@ -505,5 +505,4 @@
   ! CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
   ! CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
-
 
 END SUBROUTINE vlxqs
@@ -840,4 +839,3 @@
   ! !WRITE(*,*) 'vly 879'
 
-
 END SUBROUTINE vlyqs
