Index: LMDZ6/trunk/libf/dyn3d_common/advy.f90
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/advy.F	(revision 5084)
+++ LMDZ6/trunk/libf/dyn3d_common/advy.f90	(revision 5246)
@@ -2,419 +2,422 @@
 ! $Header$
 !
-      SUBROUTINE advy(limit,dty,pbarv,sm,s0,sx,sy,sz)
-      IMPLICIT NONE
-
-CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                                                C
-C  first-order moments (SOM) advection of tracer in Y direction  C
-C                                                                C
-C  Source : Pascal Simon ( Meteo, CNRM )                         C
-C  Adaptation : A.A. (LGGE)                                      C
-C  Derniere Modif : 15/12/94 LAST
-C                                                                C
-C  sont les arguments d'entree pour le s-pg                      C
-C                                                                C
-C  argument de sortie du s-pg                                    C
-C                                                                C
-CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C
-C  Rem : Probleme aux poles il faut reecrire ce cas specifique
-C        Attention au sens de l'indexation 
-C
-C  parametres principaux du modele
-C
-C
-      include "dimensions.h"
-      include "paramet.h"
-      include "comgeom2.h"
- 
-C  Arguments :
-C  ----------
-C  dty : frequence fictive d'appel du transport
-C  parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
-
-      INTEGER lon,lat,niv
-      INTEGER i,j,jv,k,kp,l
-      INTEGER ntra
-      PARAMETER (ntra = 1)
-
-      REAL dty
-      REAL pbarv ( iip1,jjm, llm )
-
-C  moments: SM  total mass in each grid box
-C           S0  mass of tracer in each grid box
-C           Si  1rst order moment in i direction
-C
-      REAL SM(iip1,jjp1,llm)
-     +    ,S0(iip1,jjp1,llm,ntra)
-      REAL sx(iip1,jjp1,llm,ntra)
-     +    ,sy(iip1,jjp1,llm,ntra)
-     +    ,sz(iip1,jjp1,llm,ntra)
-
-
-C  Local :
-C  ------- 
-
-C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
-C  mass fluxes in kg
-C  declaration :
-
-      REAL VGRI(iip1,0:jjp1,llm)
-
-C  Rem : UGRI et WGRI ne sont pas utilises dans 
-C  cette subroutine ( advection en y uniquement )
-C  Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
-C
-C  the moments F are similarly defined and used as temporary
-C  storage for portions of the grid boxes in transit
-C
-      REAL F0(iim,0:jjp1,ntra),FM(iim,0:jjp1)
-      REAL FX(iim,jjm,ntra),FY(iim,jjm,ntra)
-      REAL FZ(iim,jjm,ntra)
-      REAL S00(ntra)
-      REAL SM0             ! Just temporal variable
-C
-C  work arrays
-C
-      REAL ALF(iim,0:jjp1),ALF1(iim,0:jjp1)
-      REAL ALFQ(iim,0:jjp1),ALF1Q(iim,0:jjp1)
-      REAL TEMPTM          ! Just temporal variable
-c
-C  Special pour poles 
-c
-      REAL sbms,sfms,sfzs,sbmn,sfmn,sfzn
-      REAL sns0(ntra),snsz(ntra),snsm
-      REAL s1v(llm),slatv(llm)
-      REAL qy1(iim,llm,ntra),qylat(iim,llm,ntra)
-      REAL cx1(llm,ntra), cxLAT(llm,ntra)
-      REAL cy1(llm,ntra), cyLAT(llm,ntra)
-      REAL z1(iim), zcos(iim), zsin(iim)
-      real smpn,smps,s0pn,s0ps
-      REAL SSUM
-      EXTERNAL SSUM
-C
-      REAL sqi,sqf
-      LOGICAL LIMIT
-
-      lon = iim         ! rem : Il est possible qu'un pbl. arrive ici
-      lat = jjp1        ! a cause des dim. differentes entre les
-      niv=llm
-
-C
-C  the moments Fi are used as temporary storage for
-C  portions of the grid boxes in transit at the current level
-C
-C  work arrays
-C
-
-      DO l = 1,llm
-         DO j = 1,jjm
-            DO i = 1,iip1  
-            vgri (i,j,llm+1-l)=-1.*pbarv(i,j,l)  
-            enddo
-         enddo
-         do i=1,iip1
-             vgri(i,0,l) = 0.
-             vgri(i,jjp1,l) = 0.
-         enddo
-      enddo
-
-      DO 1 L=1,NIV
-C
-C  place limits on appropriate moments before transport
-C      (if flux-limiting is to be applied)
-C
-      IF(.NOT.LIMIT) GO TO 11
-C
-      DO 10 JV=1,NTRA
-      DO 10 K=1,LAT
-      DO 100 I=1,LON
-         sy(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.),
-     +                           ABS(sy(I,K,L,JV))),sy(I,K,L,JV))
- 100  CONTINUE
- 10   CONTINUE
-C
+SUBROUTINE advy(limit,dty,pbarv,sm,s0,sx,sy,sz)
+  IMPLICIT NONE
+
+  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+  !                                                                C
+  !  first-order moments (SOM) advection of tracer in Y direction  C
+  !                                                                C
+  !  Source : Pascal Simon ( Meteo, CNRM )                         C
+  !  Adaptation : A.A. (LGGE)                                      C
+  !  Derniere Modif : 15/12/94 LAST
+                                                             ! C
+  !  sont les arguments d'entree pour le s-pg                      C
+  !                                                                C
+  !  argument de sortie du s-pg                                    C
+  !                                                                C
+  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+  !
+  !  Rem : Probleme aux poles il faut reecrire ce cas specifique
+  !    Attention au sens de l'indexation
+  !
+  !  parametres principaux du modele
+  !
+  !
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom2.h"
+
+  !  Arguments :
+  !  ----------
+  !  dty : frequence fictive d'appel du transport
+  !  parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
+
+  INTEGER :: lon,lat,niv
+  INTEGER :: i,j,jv,k,kp,l
+  INTEGER :: ntra
+  PARAMETER (ntra = 1)
+
+  REAL :: dty
+  REAL :: pbarv ( iip1,jjm, llm )
+
+  !  moments: SM  total mass in each grid box
+        ! S0  mass of tracer in each grid box
+        ! Si  1rst order moment in i direction
+  !
+  REAL :: SM(iip1,jjp1,llm) &
+        ,S0(iip1,jjp1,llm,ntra)
+  REAL :: sx(iip1,jjp1,llm,ntra) &
+        ,sy(iip1,jjp1,llm,ntra) &
+        ,sz(iip1,jjp1,llm,ntra)
+
+
+  !  Local :
+  !  -------
+
+  !  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
+  !  mass fluxes in kg
+  !  declaration :
+
+  REAL :: VGRI(iip1,0:jjp1,llm)
+
+  !  Rem : UGRI et WGRI ne sont pas utilises dans
+  !  cette subroutine ( advection en y uniquement )
+  !  Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
+  !
+  !  the moments F are similarly defined and used as temporary
+  !  storage for portions of the grid boxes in transit
+  !
+  REAL :: F0(iim,0:jjp1,ntra),FM(iim,0:jjp1)
+  REAL :: FX(iim,jjm,ntra),FY(iim,jjm,ntra)
+  REAL :: FZ(iim,jjm,ntra)
+  REAL :: S00(ntra)
+  REAL :: SM0             ! Just temporal variable
+  !
+  !  work arrays
+  !
+  REAL :: ALF(iim,0:jjp1),ALF1(iim,0:jjp1)
+  REAL :: ALFQ(iim,0:jjp1),ALF1Q(iim,0:jjp1)
+  REAL :: TEMPTM          ! Just temporal variable
+  !
+  !  Special pour poles
+  !
+  REAL :: sbms,sfms,sfzs,sbmn,sfmn,sfzn
+  REAL :: sns0(ntra),snsz(ntra),snsm
+  REAL :: s1v(llm),slatv(llm)
+  REAL :: qy1(iim,llm,ntra),qylat(iim,llm,ntra)
+  REAL :: cx1(llm,ntra), cxLAT(llm,ntra)
+  REAL :: cy1(llm,ntra), cyLAT(llm,ntra)
+  REAL :: z1(iim), zcos(iim), zsin(iim)
+  real :: smpn,smps,s0pn,s0ps
+  REAL :: SSUM
+  EXTERNAL SSUM
+  !
+  REAL :: sqi,sqf
+  LOGICAL :: LIMIT
+
+  lon = iim         ! rem : Il est possible qu'un pbl. arrive ici
+  lat = jjp1        ! a cause des dim. differentes entre les
+  niv=llm
+
+  !
+  !  the moments Fi are used as temporary storage for
+  !  portions of the grid boxes in transit at the current level
+  !
+  !  work arrays
+  !
+
+  DO l = 1,llm
+     DO j = 1,jjm
+        DO i = 1,iip1
+        vgri (i,j,llm+1-l)=-1.*pbarv(i,j,l)
+        enddo
+     enddo
+     do i=1,iip1
+         vgri(i,0,l) = 0.
+         vgri(i,jjp1,l) = 0.
+     enddo
+  enddo
+
+  DO L=1,NIV
+  !
+  !  place limits on appropriate moments before transport
+  !  (if flux-limiting is to be applied)
+  !
+  IF(.NOT.LIMIT) GO TO 11
+  !
+  DO JV=1,NTRA
+  DO K=1,LAT
+  DO I=1,LON
+     sy(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.), &
+           ABS(sy(I,K,L,JV))),sy(I,K,L,JV))
+  END DO
+  END DO
+  END DO
+  !
  11   CONTINUE
-C
-C  le flux a travers le pole Nord est traite separement
-C
-      SM0=0.
-      DO 20 JV=1,NTRA
-         S00(JV)=0.
- 20   CONTINUE
-C
-      DO 21 I=1,LON
-C
-         IF(VGRI(I,0,L).LE.0.) THEN
-           FM(I,0)=-VGRI(I,0,L)*DTY
-           ALF(I,0)=FM(I,0)/SM(I,1,L)
-           SM(I,1,L)=SM(I,1,L)-FM(I,0)
-           SM0=SM0+FM(I,0)
-         ENDIF
-C
-         ALFQ(I,0)=ALF(I,0)*ALF(I,0)
-         ALF1(I,0)=1.-ALF(I,0)
-         ALF1Q(I,0)=ALF1(I,0)*ALF1(I,0)
-C
- 21   CONTINUE
-C
-      DO 22 JV=1,NTRA
-      DO 220 I=1,LON
-C
-         IF(VGRI(I,0,L).LE.0.) THEN
-C
-           F0(I,0,JV)=ALF(I,0)*
-     +               ( S0(I,1,L,JV)-ALF1(I,0)*sy(I,1,L,JV) )
-C
-           S00(JV)=S00(JV)+F0(I,0,JV)
-           S0(I,1,L,JV)=S0(I,1,L,JV)-F0(I,0,JV)
-           sy(I,1,L,JV)=ALF1Q(I,0)*sy(I,1,L,JV)
-           sx(I,1,L,JV)=ALF1 (I,0)*sx(I,1,L,JV)
-           sz(I,1,L,JV)=ALF1 (I,0)*sz(I,1,L,JV)
-C
-         ENDIF
-C
- 220  CONTINUE
- 22   CONTINUE
-C
-      DO 23 I=1,LON
-         IF(VGRI(I,0,L).GT.0.) THEN
-           FM(I,0)=VGRI(I,0,L)*DTY
-           ALF(I,0)=FM(I,0)/SM0
-         ENDIF
- 23   CONTINUE
-C
-      DO 24 JV=1,NTRA
-      DO 240 I=1,LON
-         IF(VGRI(I,0,L).GT.0.) THEN
-           F0(I,0,JV)=ALF(I,0)*S00(JV)
-         ENDIF
- 240  CONTINUE
- 24   CONTINUE
-C
-C  puts the temporary moments Fi into appropriate neighboring boxes
-C
-      DO 25 I=1,LON
-C
-         IF(VGRI(I,0,L).GT.0.) THEN
-           SM(I,1,L)=SM(I,1,L)+FM(I,0)
-           ALF(I,0)=FM(I,0)/SM(I,1,L)
-         ENDIF
-C
-         ALF1(I,0)=1.-ALF(I,0)
-C
- 25   CONTINUE
-C
-      DO 26 JV=1,NTRA
-      DO 260 I=1,LON
-C
-         IF(VGRI(I,0,L).GT.0.) THEN
-C
-         TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV)
-         S0(I,1,L,JV)=S0(I,1,L,JV)+F0(I,0,JV)
-         sy(I,1,L,JV)=ALF1(I,0)*sy(I,1,L,JV)+3.*TEMPTM
-C
-         ENDIF
-C
- 260  CONTINUE
- 26   CONTINUE
-C
-C  calculate flux and moments between adjacent boxes
-C  1- create temporary moments/masses for partial boxes in transit
-C  2- reajusts moments remaining in the box
-C
-C  flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
-C
-      DO 30 K=1,LAT-1
-      KP=K+1
-      DO 300 I=1,LON
-C
-         IF(VGRI(I,K,L).LT.0.) THEN
-           FM(I,K)=-VGRI(I,K,L)*DTY
-           ALF(I,K)=FM(I,K)/SM(I,KP,L)
-           SM(I,KP,L)=SM(I,KP,L)-FM(I,K)
-         ELSE
-           FM(I,K)=VGRI(I,K,L)*DTY
-           ALF(I,K)=FM(I,K)/SM(I,K,L)
-           SM(I,K,L)=SM(I,K,L)-FM(I,K)
-         ENDIF
-C
-         ALFQ(I,K)=ALF(I,K)*ALF(I,K)
-         ALF1(I,K)=1.-ALF(I,K)
-         ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
-C
- 300  CONTINUE
- 30   CONTINUE
-C
-      DO 31 JV=1,NTRA
-      DO 31 K=1,LAT-1
-      KP=K+1
-      DO 310 I=1,LON
-C
-         IF(VGRI(I,K,L).LT.0.) THEN
-C
-           F0(I,K,JV)=ALF (I,K)*
-     +                ( S0(I,KP,L,JV)-ALF1(I,K)*sy(I,KP,L,JV) )
-           FY(I,K,JV)=ALFQ(I,K)*sy(I,KP,L,JV)
-           FX(I,K,JV)=ALF (I,K)*sx(I,KP,L,JV)
-           FZ(I,K,JV)=ALF (I,K)*sz(I,KP,L,JV)
-C
-           S0(I,KP,L,JV)=S0(I,KP,L,JV)-F0(I,K,JV)
-           sy(I,KP,L,JV)=ALF1Q(I,K)*sy(I,KP,L,JV)
-           sx(I,KP,L,JV)=sx(I,KP,L,JV)-FX(I,K,JV)
-           sz(I,KP,L,JV)=sz(I,KP,L,JV)-FZ(I,K,JV)
-C
-         ELSE
-C
-           F0(I,K,JV)=ALF (I,K)*
-     +               ( S0(I,K,L,JV)+ALF1(I,K)*sy(I,K,L,JV) )
-           FY(I,K,JV)=ALFQ(I,K)*sy(I,K,L,JV)
-           FX(I,K,JV)=ALF(I,K)*sx(I,K,L,JV)
-           FZ(I,K,JV)=ALF(I,K)*sz(I,K,L,JV)
-C
-           S0(I,K,L,JV)=S0(I,K,L,JV)-F0(I,K,JV)
-           sy(I,K,L,JV)=ALF1Q(I,K)*sy(I,K,L,JV)
-           sx(I,K,L,JV)=sx(I,K,L,JV)-FX(I,K,JV)
-           sz(I,K,L,JV)=sz(I,K,L,JV)-FZ(I,K,JV)
-C
-         ENDIF
-C
- 310  CONTINUE
- 31   CONTINUE
-C
-C  puts the temporary moments Fi into appropriate neighboring boxes
-C
-      DO 32 K=1,LAT-1
-      KP=K+1
-      DO 320 I=1,LON
-C
-         IF(VGRI(I,K,L).LT.0.) THEN
-           SM(I,K,L)=SM(I,K,L)+FM(I,K)
-           ALF(I,K)=FM(I,K)/SM(I,K,L)
-         ELSE
-           SM(I,KP,L)=SM(I,KP,L)+FM(I,K)
-           ALF(I,K)=FM(I,K)/SM(I,KP,L)
-         ENDIF
-C
-         ALF1(I,K)=1.-ALF(I,K)
-C
- 320  CONTINUE
- 32   CONTINUE
-C
-      DO 33 JV=1,NTRA
-      DO 33 K=1,LAT-1
-      KP=K+1
-      DO 330 I=1,LON
-C
-         IF(VGRI(I,K,L).LT.0.) THEN
-C
-         TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
-         S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,K,JV)
-         sy(I,K,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*sy(I,K,L,JV)
-     +               +3.*TEMPTM
-         sx(I,K,L,JV)=sx(I,K,L,JV)+FX(I,K,JV)
-         sz(I,K,L,JV)=sz(I,K,L,JV)+FZ(I,K,JV)
-C
-         ELSE
-C
-         TEMPTM=ALF(I,K)*S0(I,KP,L,JV)-ALF1(I,K)*F0(I,K,JV)
-         S0(I,KP,L,JV)=S0(I,KP,L,JV)+F0(I,K,JV)
-         sy(I,KP,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*sy(I,KP,L,JV)
-     +                +3.*TEMPTM
-         sx(I,KP,L,JV)=sx(I,KP,L,JV)+FX(I,K,JV)
-         sz(I,KP,L,JV)=sz(I,KP,L,JV)+FZ(I,K,JV)
-C
-         ENDIF
-C
- 330  CONTINUE
- 33   CONTINUE
-C
-C  traitement special pour le pole Sud (idem pole Nord)
-C
-      K=LAT
-C
-      SM0=0.
-      DO 40 JV=1,NTRA
-         S00(JV)=0.
- 40   CONTINUE
-C
-      DO 41 I=1,LON
-C
-         IF(VGRI(I,K,L).GE.0.) THEN
-           FM(I,K)=VGRI(I,K,L)*DTY
-           ALF(I,K)=FM(I,K)/SM(I,K,L)
-           SM(I,K,L)=SM(I,K,L)-FM(I,K)
-           SM0=SM0+FM(I,K)
-         ENDIF
-C
-         ALFQ(I,K)=ALF(I,K)*ALF(I,K)
-         ALF1(I,K)=1.-ALF(I,K)
-         ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
-C
- 41   CONTINUE
-C
-      DO 42 JV=1,NTRA
-      DO 420 I=1,LON
-C
-         IF(VGRI(I,K,L).GE.0.) THEN
-           F0 (I,K,JV)=ALF(I,K)*
-     +                ( S0(I,K,L,JV)+ALF1(I,K)*sy(I,K,L,JV) )
-           S00(JV)=S00(JV)+F0(I,K,JV)
-C
-           S0(I,K,L,JV)=S0 (I,K,L,JV)-F0 (I,K,JV)
-           sy(I,K,L,JV)=ALF1Q(I,K)*sy(I,K,L,JV)
-           sx(I,K,L,JV)=ALF1(I,K)*sx(I,K,L,JV)
-           sz(I,K,L,JV)=ALF1(I,K)*sz(I,K,L,JV)
-         ENDIF
-C
- 420  CONTINUE
- 42   CONTINUE
-C
-      DO 43 I=1,LON
-         IF(VGRI(I,K,L).LT.0.) THEN
-           FM(I,K)=-VGRI(I,K,L)*DTY
-           ALF(I,K)=FM(I,K)/SM0
-         ENDIF
- 43   CONTINUE
-C
-      DO 44 JV=1,NTRA
-      DO 440 I=1,LON
-         IF(VGRI(I,K,L).LT.0.) THEN
-           F0(I,K,JV)=ALF(I,K)*S00(JV)
-         ENDIF
- 440  CONTINUE
- 44   CONTINUE
-C
-C  puts the temporary moments Fi into appropriate neighboring boxes
-C
-      DO 45 I=1,LON
-C
-         IF(VGRI(I,K,L).LT.0.) THEN
-           SM(I,K,L)=SM(I,K,L)+FM(I,K)
-           ALF(I,K)=FM(I,K)/SM(I,K,L)
-         ENDIF
-C
-         ALF1(I,K)=1.-ALF(I,K)
-C
- 45   CONTINUE
-C
-      DO 46 JV=1,NTRA
-      DO 460 I=1,LON
-C
-         IF(VGRI(I,K,L).LT.0.) THEN
-C
-         TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
-         S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,K,JV)
-         sy(I,K,L,JV)=ALF1(I,K)*sy(I,K,L,JV)+3.*TEMPTM
-C
-         ENDIF
-C
- 460  CONTINUE
- 46   CONTINUE
-C
- 1    CONTINUE
-C
-      RETURN
-      END
-
+  !
+  !  le flux a travers le pole Nord est traite separement
+  !
+  SM0=0.
+  DO JV=1,NTRA
+     S00(JV)=0.
+  END DO
+  !
+  DO I=1,LON
+  !
+     IF(VGRI(I,0,L).LE.0.) THEN
+       FM(I,0)=-VGRI(I,0,L)*DTY
+       ALF(I,0)=FM(I,0)/SM(I,1,L)
+       SM(I,1,L)=SM(I,1,L)-FM(I,0)
+       SM0=SM0+FM(I,0)
+     ENDIF
+  !
+     ALFQ(I,0)=ALF(I,0)*ALF(I,0)
+     ALF1(I,0)=1.-ALF(I,0)
+     ALF1Q(I,0)=ALF1(I,0)*ALF1(I,0)
+  !
+  END DO
+  !
+  DO JV=1,NTRA
+  DO I=1,LON
+  !
+     IF(VGRI(I,0,L).LE.0.) THEN
+  !
+       F0(I,0,JV)=ALF(I,0)* &
+             ( S0(I,1,L,JV)-ALF1(I,0)*sy(I,1,L,JV) )
+  !
+       S00(JV)=S00(JV)+F0(I,0,JV)
+       S0(I,1,L,JV)=S0(I,1,L,JV)-F0(I,0,JV)
+       sy(I,1,L,JV)=ALF1Q(I,0)*sy(I,1,L,JV)
+       sx(I,1,L,JV)=ALF1 (I,0)*sx(I,1,L,JV)
+       sz(I,1,L,JV)=ALF1 (I,0)*sz(I,1,L,JV)
+  !
+     ENDIF
+  !
+  END DO
+  END DO
+  !
+  DO I=1,LON
+     IF(VGRI(I,0,L).GT.0.) THEN
+       FM(I,0)=VGRI(I,0,L)*DTY
+       ALF(I,0)=FM(I,0)/SM0
+     ENDIF
+  END DO
+  !
+  DO JV=1,NTRA
+  DO I=1,LON
+     IF(VGRI(I,0,L).GT.0.) THEN
+       F0(I,0,JV)=ALF(I,0)*S00(JV)
+     ENDIF
+  END DO
+  END DO
+  !
+  !  puts the temporary moments Fi into appropriate neighboring boxes
+  !
+  DO I=1,LON
+  !
+     IF(VGRI(I,0,L).GT.0.) THEN
+       SM(I,1,L)=SM(I,1,L)+FM(I,0)
+       ALF(I,0)=FM(I,0)/SM(I,1,L)
+     ENDIF
+  !
+     ALF1(I,0)=1.-ALF(I,0)
+  !
+  END DO
+  !
+  DO JV=1,NTRA
+  DO I=1,LON
+  !
+     IF(VGRI(I,0,L).GT.0.) THEN
+  !
+     TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV)
+     S0(I,1,L,JV)=S0(I,1,L,JV)+F0(I,0,JV)
+     sy(I,1,L,JV)=ALF1(I,0)*sy(I,1,L,JV)+3.*TEMPTM
+  !
+     ENDIF
+  !
+  END DO
+  END DO
+  !
+  !  calculate flux and moments between adjacent boxes
+  !  1- create temporary moments/masses for partial boxes in transit
+  !  2- reajusts moments remaining in the box
+  !
+  !  flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
+  !
+  DO K=1,LAT-1
+  KP=K+1
+  DO I=1,LON
+  !
+     IF(VGRI(I,K,L).LT.0.) THEN
+       FM(I,K)=-VGRI(I,K,L)*DTY
+       ALF(I,K)=FM(I,K)/SM(I,KP,L)
+       SM(I,KP,L)=SM(I,KP,L)-FM(I,K)
+     ELSE
+       FM(I,K)=VGRI(I,K,L)*DTY
+       ALF(I,K)=FM(I,K)/SM(I,K,L)
+       SM(I,K,L)=SM(I,K,L)-FM(I,K)
+     ENDIF
+  !
+     ALFQ(I,K)=ALF(I,K)*ALF(I,K)
+     ALF1(I,K)=1.-ALF(I,K)
+     ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
+  !
+  END DO
+  END DO
+  !
+  DO JV=1,NTRA
+  DO K=1,LAT-1
+  KP=K+1
+  DO I=1,LON
+  !
+     IF(VGRI(I,K,L).LT.0.) THEN
+  !
+       F0(I,K,JV)=ALF (I,K)* &
+             ( S0(I,KP,L,JV)-ALF1(I,K)*sy(I,KP,L,JV) )
+       FY(I,K,JV)=ALFQ(I,K)*sy(I,KP,L,JV)
+       FX(I,K,JV)=ALF (I,K)*sx(I,KP,L,JV)
+       FZ(I,K,JV)=ALF (I,K)*sz(I,KP,L,JV)
+  !
+       S0(I,KP,L,JV)=S0(I,KP,L,JV)-F0(I,K,JV)
+       sy(I,KP,L,JV)=ALF1Q(I,K)*sy(I,KP,L,JV)
+       sx(I,KP,L,JV)=sx(I,KP,L,JV)-FX(I,K,JV)
+       sz(I,KP,L,JV)=sz(I,KP,L,JV)-FZ(I,K,JV)
+  !
+     ELSE
+  !
+       F0(I,K,JV)=ALF (I,K)* &
+             ( S0(I,K,L,JV)+ALF1(I,K)*sy(I,K,L,JV) )
+       FY(I,K,JV)=ALFQ(I,K)*sy(I,K,L,JV)
+       FX(I,K,JV)=ALF(I,K)*sx(I,K,L,JV)
+       FZ(I,K,JV)=ALF(I,K)*sz(I,K,L,JV)
+  !
+       S0(I,K,L,JV)=S0(I,K,L,JV)-F0(I,K,JV)
+       sy(I,K,L,JV)=ALF1Q(I,K)*sy(I,K,L,JV)
+       sx(I,K,L,JV)=sx(I,K,L,JV)-FX(I,K,JV)
+       sz(I,K,L,JV)=sz(I,K,L,JV)-FZ(I,K,JV)
+  !
+     ENDIF
+  !
+  END DO
+  END DO
+  END DO
+  !
+  !  puts the temporary moments Fi into appropriate neighboring boxes
+  !
+  DO K=1,LAT-1
+  KP=K+1
+  DO I=1,LON
+  !
+     IF(VGRI(I,K,L).LT.0.) THEN
+       SM(I,K,L)=SM(I,K,L)+FM(I,K)
+       ALF(I,K)=FM(I,K)/SM(I,K,L)
+     ELSE
+       SM(I,KP,L)=SM(I,KP,L)+FM(I,K)
+       ALF(I,K)=FM(I,K)/SM(I,KP,L)
+     ENDIF
+  !
+     ALF1(I,K)=1.-ALF(I,K)
+  !
+  END DO
+  END DO
+  !
+  DO JV=1,NTRA
+  DO K=1,LAT-1
+  KP=K+1
+  DO I=1,LON
+  !
+     IF(VGRI(I,K,L).LT.0.) THEN
+  !
+     TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
+     S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,K,JV)
+     sy(I,K,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*sy(I,K,L,JV) &
+           +3.*TEMPTM
+     sx(I,K,L,JV)=sx(I,K,L,JV)+FX(I,K,JV)
+     sz(I,K,L,JV)=sz(I,K,L,JV)+FZ(I,K,JV)
+  !
+     ELSE
+  !
+     TEMPTM=ALF(I,K)*S0(I,KP,L,JV)-ALF1(I,K)*F0(I,K,JV)
+     S0(I,KP,L,JV)=S0(I,KP,L,JV)+F0(I,K,JV)
+     sy(I,KP,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*sy(I,KP,L,JV) &
+           +3.*TEMPTM
+     sx(I,KP,L,JV)=sx(I,KP,L,JV)+FX(I,K,JV)
+     sz(I,KP,L,JV)=sz(I,KP,L,JV)+FZ(I,K,JV)
+  !
+     ENDIF
+  !
+  END DO
+  END DO
+  END DO
+  !
+  !  traitement special pour le pole Sud (idem pole Nord)
+  !
+  K=LAT
+  !
+  SM0=0.
+  DO JV=1,NTRA
+     S00(JV)=0.
+  END DO
+  !
+  DO I=1,LON
+  !
+     IF(VGRI(I,K,L).GE.0.) THEN
+       FM(I,K)=VGRI(I,K,L)*DTY
+       ALF(I,K)=FM(I,K)/SM(I,K,L)
+       SM(I,K,L)=SM(I,K,L)-FM(I,K)
+       SM0=SM0+FM(I,K)
+     ENDIF
+  !
+     ALFQ(I,K)=ALF(I,K)*ALF(I,K)
+     ALF1(I,K)=1.-ALF(I,K)
+     ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
+  !
+  END DO
+  !
+  DO JV=1,NTRA
+  DO I=1,LON
+  !
+     IF(VGRI(I,K,L).GE.0.) THEN
+       F0 (I,K,JV)=ALF(I,K)* &
+             ( S0(I,K,L,JV)+ALF1(I,K)*sy(I,K,L,JV) )
+       S00(JV)=S00(JV)+F0(I,K,JV)
+  !
+       S0(I,K,L,JV)=S0 (I,K,L,JV)-F0 (I,K,JV)
+       sy(I,K,L,JV)=ALF1Q(I,K)*sy(I,K,L,JV)
+       sx(I,K,L,JV)=ALF1(I,K)*sx(I,K,L,JV)
+       sz(I,K,L,JV)=ALF1(I,K)*sz(I,K,L,JV)
+     ENDIF
+  !
+  END DO
+  END DO
+  !
+  DO I=1,LON
+     IF(VGRI(I,K,L).LT.0.) THEN
+       FM(I,K)=-VGRI(I,K,L)*DTY
+       ALF(I,K)=FM(I,K)/SM0
+     ENDIF
+  END DO
+  !
+  DO JV=1,NTRA
+  DO I=1,LON
+     IF(VGRI(I,K,L).LT.0.) THEN
+       F0(I,K,JV)=ALF(I,K)*S00(JV)
+     ENDIF
+  END DO
+  END DO
+  !
+  !  puts the temporary moments Fi into appropriate neighboring boxes
+  !
+  DO I=1,LON
+  !
+     IF(VGRI(I,K,L).LT.0.) THEN
+       SM(I,K,L)=SM(I,K,L)+FM(I,K)
+       ALF(I,K)=FM(I,K)/SM(I,K,L)
+     ENDIF
+  !
+     ALF1(I,K)=1.-ALF(I,K)
+  !
+  END DO
+  !
+  DO JV=1,NTRA
+  DO I=1,LON
+  !
+     IF(VGRI(I,K,L).LT.0.) THEN
+  !
+     TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
+     S0(I,K,L,JV)=S0(I,K,L,JV)+F0(I,K,JV)
+     sy(I,K,L,JV)=ALF1(I,K)*sy(I,K,L,JV)+3.*TEMPTM
+  !
+     ENDIF
+  !
+  END DO
+  END DO
+  !
+  END DO
+  !
+  RETURN
+END SUBROUTINE advy
+
