Index: /LMDZ.3.3/branches/rel-LF/libf/dyn3d/vlsplt.F
===================================================================
--- /LMDZ.3.3/branches/rel-LF/libf/dyn3d/vlsplt.F	(revision 357)
+++ /LMDZ.3.3/branches/rel-LF/libf/dyn3d/vlsplt.F	(revision 358)
@@ -26,6 +26,8 @@
 c   ----------
       REAL masse(ip1jmp1,llm),pente_max
+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 w(ip1jmp1,llm),pdt
 c
@@ -33,16 +35,20 @@
 c   ---------
 c
-      INTEGER ij,l
-c
-      REAL zm(ip1jmp1,llm)
+      INTEGER i,ij,l,j,ii
+      INTEGER ijlqmin,iqmin,jqmin,lqmin
+c
+      REAL zm(ip1jmp1,llm),newmasse
       REAL mu(ip1jmp1,llm)
       REAL mv(ip1jm,llm)
       REAL mw(ip1jmp1,llm+1)
-      REAL zq(ip1jmp1,llm)
-      REAL temps1,temps2,temps3
+      REAL zq(ip1jmp1,llm),zz
+      REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
+      REAL second,temps0,temps1,temps2,temps3
+      REAL ztemps1,ztemps2,ztemps3
       REAL zzpbar, zzw
       LOGICAL testcpu
       SAVE testcpu
       SAVE temps1,temps2,temps3
+      INTEGER iminn,imaxx
 
       REAL qmin,qmax
@@ -51,5 +57,4 @@
       DATA temps1,temps2,temps3/0.,0.,0./
 
-c      PRINT*,'Debut vlsplt version debug sans vly'
 
         zzpbar = 0.5 * pdt
@@ -70,36 +75,29 @@
          mw(ij,llm+1)=0.
       ENDDO
-
+      
       CALL SCOPY(ijp1llm,q,1,zq,1)
       CALL SCOPY(ijp1llm,masse,1,zm,1)
 
-c      call minmaxq(zq,qmin,qmax,'avant vlx     ')
+	print*,'Entree vlx1'
+c	call minmaxq(zq,qmin,qmax,'avant vlx     ')
       call vlx(zq,pente_max,zm,mu)
-
-
-c     call minmaxq(zq,qmin,qmax,'avant vly     ')
-
+	print*,'Sortie vlx1'
+c	call minmaxq(zq,qmin,qmax,'apres vlx1    ')
+
+	 print*,'Entree vly1'
       call vly(zq,pente_max,zm,mv)
-
-
-c      call minmaxq(zq,qmin,qmax,'avant vlz     ')
-
+c	call minmaxq(zq,qmin,qmax,'apres vly1     ')
+	print*,'Sortie vly1'
       call vlz(zq,pente_max,zm,mw)
-
-
-c     call minmaxq(zq,qmin,qmax,'avant vly     ')
-c     call minmaxq(zm,qmin,qmax,'M avant vly     ')
+c	call minmaxq(zq,qmin,qmax,'apres vlz     ')
+
 
       call vly(zq,pente_max,zm,mv)
-
-
-c     call minmaxq(zq,qmin,qmax,'avant vlx     ')
-c     call minmaxq(zm,qmin,qmax,'M avant vlx     ')
+c	call minmaxq(zq,qmin,qmax,'apres vly     ')
+
 
       call vlx(zq,pente_max,zm,mu)
-
-c     call minmaxq(zq,qmin,qmax,'apres vlx     ')
-c     call minmaxq(zm,qmin,qmax,'M apres vlx     ')
-
+c	call minmaxq(zq,qmin,qmax,'apres vlx2    ')
+	
 
       DO l=1,llm
@@ -115,5 +113,5 @@
       END
       SUBROUTINE vlx(q,pente_max,masse,u_m)
-c
+
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
 c
@@ -137,6 +135,7 @@
 c   ----------
       REAL masse(ip1jmp1,llm),pente_max
-      REAL u_m( ip1jmp1,llm )
+      REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
       REAL q(ip1jmp1,llm)
+      REAL w(ip1jmp1,llm)
 c
 c      Local 
@@ -147,17 +146,18 @@
 c
       REAL new_m,zu_m,zdum(ip1jmp1,llm)
-      REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
+      REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
       REAL zz(ip1jmp1)
       REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
       REAL u_mq(ip1jmp1,llm)
 
-      Logical first,testcpu
+      Logical extremum,first,testcpu
       SAVE first,testcpu
 
       REAL      SSUM
       EXTERNAL  SSUM
-      REAL temps1,temps2,temps3,temps4,temps5
-      SAVE temps1,temps2,temps3,temps4,temps5
-
+      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
+      SAVE temps0,temps1,temps2,temps3,temps4,temps5
+
+      REAL z1,z2,z3
 
       DATA first,testcpu/.true.,.false./
@@ -176,5 +176,5 @@
 
       IF (pente_max.gt.-1.e-5) THEN
-c     IF (pente_max.gt.10) THEN
+c       IF (pente_max.gt.10) THEN
 
 c   calcul des pentes avec limitation, Van Leer scheme I:
@@ -230,4 +230,5 @@
 
          ENDDO ! l=1,llm
+	print*,'Ok calcul des pentes'
 
       ELSE ! (pente_max.lt.-1.e-5)
@@ -266,5 +267,4 @@
             dxq(ij-iim,l)=dxq(ij,l)
          ENDDO
-
          DO ij=1,ip1jmp1
             iadvplus(ij,l)=0
@@ -273,4 +273,5 @@
       ENDDO
 
+	 print*,'Bouclage en iip1'
 
 c   calcul des flux a gauche et a droite
@@ -294,6 +295,9 @@
 c   on cumule le flux correspondant a toutes les mailles dont la masse
 c   au travers de la paroi pENDant le pas de temps.
+	print*,'Cumule ....'
+
       DO l=1,llm
        DO ij=iip2,ip1jm-1
+c	print*,'masse(',ij,')=',masse(ij,l)
           IF (u_m(ij,l).gt.0.) THEN
              zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
@@ -306,6 +310,7 @@
       ENDDO
 #endif
-
-
+c	stop
+
+c	go to 9999
 c   detection des points ou on advecte plus que la masse de la
 c   maille
@@ -318,4 +323,5 @@
          ENDDO
       ENDDO
+	print*,'Ok test 1'
       DO l=1,llm
        DO ij=iip1+iip1,ip1jm,iip1
@@ -323,5 +329,5 @@
        ENDDO
       ENDDO
-
+	 print*,'Ok test 2'
 
 
@@ -342,6 +348,6 @@
 
       IF(n0.gt.1) THEN
-ccc      PRINT*,'Nombre de points pour lesquels on advect plus que le'
-ccc     &       ,'contenu de la maille : ',n0
+      PRINT*,'Nombre de points pour lesquels on advect plus que le'
+     &       ,'contenu de la maille : ',n0
 
          DO l=1,llm
@@ -395,9 +401,9 @@
          ENDDO
       ENDIF  ! n0.gt.0 
-
+9999    continue
 
 
 c   bouclage en latitude
-
+	print*,'Avant bouclage en latitude'
       DO l=1,llm
         DO ij=iip1+iip1,ip1jm,iip1
@@ -423,5 +429,4 @@
          ENDDO
       ENDDO
-
 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)
@@ -456,5 +461,5 @@
       REAL masse(ip1jmp1,llm),pente_max
       REAL masse_adv_v( ip1jm,llm)
-      REAL q(ip1jmp1,llm)
+      REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
 c
 c      Local 
@@ -464,16 +469,17 @@
 c
       REAL airej2,airejjm,airescb(iim),airesch(iim)
-      REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
+      REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
       REAL adyqv(ip1jm),dyqmax(ip1jmp1)
       REAL qbyv(ip1jm,llm)
 
-      REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
+      REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
 c     REAL newq,oldmasse
-      Logical first,testcpu
-      REAL temps0,temps1,temps2,temps3,temps4,temps5
+      Logical extremum,first,testcpu
+      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
       SAVE temps0,temps1,temps2,temps3,temps4,temps5
       SAVE first,testcpu
 
       REAL convpn,convps,convmpn,convmps
+      real massepn,masseps,qpn,qps
       REAL sinlon(iip1),sinlondlon(iip1)
       REAL coslon(iip1),coslondlon(iip1)
@@ -506,5 +512,5 @@
 
 c
-
+	PRINT*,'CALCUL EN LATITUDE'
 
       DO l = 1, llm
@@ -565,4 +571,5 @@
 c   calcul des pentes limites aux poles
 
+      goto 8888
       fn=1.
       fs=1.
@@ -578,4 +585,9 @@
          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
 
@@ -682,20 +694,34 @@
          ENDDO
 c.-. ancienne version
-         convpn=SSUM(iim,qbyv(1,l),1)/apoln
-         convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
-         DO ij = 1,iip1
-            newmasse=masse(ij,l)+convmpn*aire(ij)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
-     &               newmasse
-            masse(ij,l)=newmasse
-         ENDDO
-         convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
-         convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
-         DO ij = ip1jm+1,ip1jmp1
-            newmasse=masse(ij,l)+convmps*aire(ij)
-            q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
-     &               newmasse
-            masse(ij,l)=newmasse
-         ENDDO
+c        convpn=SSUM(iim,qbyv(1,l),1)/apoln
+c        convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
+
+         convpn=SSUM(iim,qbyv(1,l),1)
+         convmpn=ssum(iim,masse_adv_v(1,l),1)
+         massepn=ssum(iim,masse(1,l),1)
+         qpn=0.
+         do ij=1,iim
+            qpn=qpn+masse(ij,l)*q(ij,l)
+         enddo
+         qpn=(qpn+convpn)/(massepn+convmpn)
+         do ij=1,iip1
+            q(ij,l)=qpn
+         enddo
+
+c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
+c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
+
+         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)
+         qps=0.
+         do ij = ip1jm+1,ip1jmp1-1
+            qps=qps+masse(ij,l)*q(ij,l)
+         enddo
+         qps=(qps+convps)/(masseps+convmps)
+         do ij=ip1jm+1,ip1jmp1
+            q(ij,l)=qps
+         enddo
+
 c.-. fin ancienne version
 
@@ -852,32 +878,76 @@
       RETURN
       END
-      SUBROUTINE minmaxq(zq,qmin,qmax,comment)
+c      SUBROUTINE minmaxq(zq,qmin,qmax,comment)
+c
+c#include "dimensions.h"
+c#include "paramet.h"
+
+c      CHARACTER*(*) comment
+c      real qmin,qmax
+c      real zq(ip1jmp1,llm)
+
+c      INTEGER jadrs(ip1jmp1), jbad, k, i
+
+
+c      DO k = 1, llm
+c         jbad = 0
+c         DO i = 1, ip1jmp1
+c         IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
+c            jbad = jbad + 1
+c            jadrs(jbad) = i
+c         ENDIF
+c         ENDDO
+c         IF (jbad.GT.0) THEN
+c         PRINT*, comment
+c         DO i = 1, jbad
+cc            PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
+c         ENDDO
+c         ENDIF
+c      ENDDO
+
+c      return
+c      end
+      subroutine minmaxq(zq,qmin,qmax,comment)
 
 #include "dimensions.h"
 #include "paramet.h"
 
-      CHARACTER*(*) comment
+      character*20 comment
       real qmin,qmax
       real zq(ip1jmp1,llm)
-
-      INTEGER jadrs(ip1jmp1), jbad, k, i
-
-      DO k = 1, llm
-         jbad = 0
-         DO i = 1, ip1jmp1
-         IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
-            jbad = jbad + 1
-            jadrs(jbad) = i
-         ENDIF
-         ENDDO
-         IF (jbad.GT.0) THEN
-         PRINT*, comment
-         DO i = 1, jbad
-cc            PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
-         ENDDO
-         ENDIF
-      ENDDO
+      real zzq(iip1,jjp1,llm)
+
+      integer imin,jmin,lmin,ijlmin
+      integer imax,jmax,lmax,ijlmax
+
+      integer ismin,ismax
+
+      call scopy (ip1jmp1*llm,zq,1,zzq,1)
+
+      ijlmin=ismin(ijp1llm,zq,1)
+      lmin=(ijlmin-1)/ip1jmp1+1
+      ijlmin=ijlmin-(lmin-1.)*ip1jmp1
+      jmin=(ijlmin-1)/iip1+1
+      imin=ijlmin-(jmin-1.)*iip1
+      zqmin=zq(ijlmin,lmin)
+
+      ijlmax=ismax(ijp1llm,zq,1)
+      lmax=(ijlmax-1)/ip1jmp1+1
+      ijlmax=ijlmax-(lmax-1.)*ip1jmp1
+      jmax=(ijlmax-1)/iip1+1
+      imax=ijlmax-(jmax-1.)*iip1
+      zqmax=zq(ijlmax,lmax)
+
+       if(zqmin.lt.qmin) 
+c     s     write(*,9999) comment,
+     s     write(*,*) comment,
+     s     imin,jmin,lmin,zqmin,zzq(imin,jmin,lmin)
+       if(zqmax.gt.qmax) 
+c     s     write(*,9999) comment,
+     s     write(*,*) comment,
+     s     imax,jmax,lmax,zqmax,zzq(imax,jmax,lmax)
 
       return
+9999  format(a20,'  q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
       end
 
