!
! $Header$
!
SUBROUTINE ADVYP(LIMIT,DTY,PBARV,SM,S0,SSX,SY,SZ &
        ,SSXX,SSXY,SSXZ,SYY,SYZ,SZZ,ntra )
  IMPLICIT NONE
  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  !                                                                 C
  !  second-order moments (SOM) advection of tracer in Y direction  C
  !                                                                 C
  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                                                             ! C
  !  Source : Pascal Simon ( Meteo, CNRM )                         C
  !  Adaptation : A.A. (LGGE)                                      C
  !  Derniere Modif : 19/10/95 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 "comgeom.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 :: SSX(iip1,jjp1,llm,ntra) &
        ,SY(iip1,jjp1,llm,ntra) &
        ,SZ(iip1,jjp1,llm,ntra) &
        ,SSXX(iip1,jjp1,llm,ntra) &
        ,SSXY(iip1,jjp1,llm,ntra) &
        ,SSXZ(iip1,jjp1,llm,ntra) &
        ,SYY(iip1,jjp1,llm,ntra) &
        ,SYZ(iip1,jjp1,llm,ntra) &
        ,SZZ(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
  !
  !  the moments Fij are used as temporary storage for
  !  portions of the grid boxes in transit at the current level
  !
  !  work arrays
  !
  !
  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 :: FXX(iim,jjm,ntra),FXY(iim,jjm,ntra)
  REAL :: FXZ(iim,jjm,ntra),FYY(iim,jjm,ntra)
  REAL :: FYZ(iim,jjm,ntra),FZZ(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 :: ALF2(iim,0:jjp1),ALF3(iim,0:jjp1)
  REAL :: ALF4(iim,0:jjp1)
  REAL :: TEMPTM          ! Just temporal variable
  REAL :: SLPMAX,S1MAX,S1NEW,S2NEW
  !
  !  Special pour poles
  !
  REAL :: sbms,sfms,sfzs,sbmn,sfmn,sfzn
  REAL :: sns0(ntra),snsz(ntra),snsm
  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 :: 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         !       tab. S et VGRI

  !-----------------------------------------------------------------
  ! initialisations

  sbms = 0.
  sfms = 0.
  sfzs = 0.
  sbmn = 0.
  sfmn = 0.
  sfzn = 0.

  !-----------------------------------------------------------------
  ! *** Test : diag de la qtite totale de traceur dans
         ! l'atmosphere avant l'advection en Y
  !
  sqi = 0.
  sqf = 0.

  DO l = 1,llm
     DO j = 1,jjp1
       DO i = 1,iim
          sqi = sqi + S0(i,j,l,ntra)
       END DO
     END DO
  END DO
  PRINT*,'---------- DIAG DANS ADVY - ENTREE --------'
  PRINT*,'sqi=',sqi

  !-----------------------------------------------------------------
  !  Interface : adaptation nouveau modele
  !  -------------------------------------
  !
  !  Conversion des flux de masses en kg
  !-AA 20/10/94  le signe -1 est necessaire car indexation opposee

  DO l = 1,llm
     DO j = 1,jjm
        DO i = 1,iip1
        vgri (i,j,llm+1-l)=-1.*pbarv (i,j,l)
        END DO
     END DO
  END DO

  !AA Initialisation de flux fictifs aux bords sup. des boites pol.

  DO l = 1,llm
     DO i = 1,iip1
         vgri(i,0,l) = 0.
         vgri(i,jjp1,l) = 0.
     ENDDO
  ENDDO
  !
  !----------------- START HERE -----------------------
  !  boucle sur les niveaux
  !
  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
     IF(S0(I,K,L,JV).GT.0.) THEN
       SLPMAX=AMAX1(S0(I,K,L,JV),0.)
       S1MAX=1.5*SLPMAX
       S1NEW=AMIN1(S1MAX,AMAX1(-S1MAX,SY(I,K,L,JV)))
       S2NEW=AMIN1( 2.*SLPMAX-ABS(S1NEW)/3. , &
             AMAX1(ABS(S1NEW)-SLPMAX,SYY(I,K,L,JV)) )
       SY (I,K,L,JV)=S1NEW
       SYY(I,K,L,JV)=S2NEW
   SSXY(I,K,L,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,SSXY(I,K,L,JV)))
   SYZ(I,K,L,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,SYZ(I,K,L,JV)))
     ELSE
       SY (I,K,L,JV)=0.
       SYY(I,K,L,JV)=0.
       SSXY(I,K,L,JV)=0.
       SYZ(I,K,L,JV)=0.
     ENDIF
  END DO
  END DO
  END DO
  !
 11   CONTINUE
  !
  !  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)
     ALF2(I,0)=ALF1(I,0)-ALF(I,0)
     ALF3(I,0)=ALF(I,0)*ALFQ(I,0)
     ALF4(I,0)=ALF1(I,0)*ALF1Q(I,0)
  !
  END DO
  ! print*,'ADVYP 21'
  !
  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)-ALF2(I,0)*SYY(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)+3.*ALF(I,0)*SYY(I,1,L,JV))
       SYY(I,1,L,JV)=ALF4 (I,0)*SYY(I,1,L,JV)
       SSX (I,1,L,JV)=ALF1 (I,0)* &
             (SSX(I,1,L,JV)+ALF(I,0)*SSXY(I,1,L,JV) )
       SZ (I,1,L,JV)=ALF1 (I,0)* &
             (SZ(I,1,L,JV)+ALF(I,0)*SSXZ(I,1,L,JV) )
       SSXX(I,1,L,JV)=ALF1 (I,0)*SSXX(I,1,L,JV)
       SSXZ(I,1,L,JV)=ALF1 (I,0)*SSXZ(I,1,L,JV)
       SZZ(I,1,L,JV)=ALF1 (I,0)*SZZ(I,1,L,JV)
       SSXY(I,1,L,JV)=ALF1Q(I,0)*SSXY(I,1,L,JV)
       SYZ(I,1,L,JV)=ALF1Q(I,0)*SYZ(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
  !
  ! print*,'av ADVYP 25'
  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
  !
     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)
     ALF2(I,0)=ALF1(I,0)-ALF(I,0)
     ALF3(I,0)=ALF1(I,0)*ALF(I,0)
  !
  END DO
  ! print*,'av ADVYP 25'
  !
  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)
     SYY(I,1,L,JV)=ALF1Q(I,0)*SYY(I,1,L,JV) &
           +5.*( ALF3 (I,0)*SY (I,1,L,JV)-ALF2(I,0)*TEMPTM )
     SY (I,1,L,JV)=ALF1 (I,0)*SY (I,1,L,JV)+3.*TEMPTM
  SSXY(I,1,L,JV)=ALF1 (I,0)*SSXY(I,1,L,JV)+3.*ALF(I,0)*SSX(I,1,L,JV)
  SYZ(I,1,L,JV)=ALF1 (I,0)*SYZ(I,1,L,JV)+3.*ALF(I,0)*SZ(I,1,L,JV)
  !
     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
  !
  ! print*,'av ADVYP 30'
  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)
     ALF2(I,K)=ALF1(I,K)-ALF(I,K)
     ALF3(I,K)=ALF(I,K)*ALFQ(I,K)
     ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K)
  !
  END DO
  END DO
  ! print*,'ap ADVYP 30'
  !
  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)-ALF2(I,K)*SYY(I,KP,L,JV) ) )
       FY (I,K,JV)=ALFQ(I,K)* &
             (SY(I,KP,L,JV)-3.*ALF1(I,K)*SYY(I,KP,L,JV))
       FYY(I,K,JV)=ALF3(I,K)*SYY(I,KP,L,JV)
       FX (I,K,JV)=ALF (I,K)* &
             (SSX(I,KP,L,JV)-ALF1(I,K)*SSXY(I,KP,L,JV))
       FZ (I,K,JV)=ALF (I,K)* &
             (SZ(I,KP,L,JV)-ALF1(I,K)*SYZ(I,KP,L,JV))
       FXY(I,K,JV)=ALFQ(I,K)*SSXY(I,KP,L,JV)
       FYZ(I,K,JV)=ALFQ(I,K)*SYZ(I,KP,L,JV)
       FXX(I,K,JV)=ALF (I,K)*SSXX(I,KP,L,JV)
       FXZ(I,K,JV)=ALF (I,K)*SSXZ(I,KP,L,JV)
       FZZ(I,K,JV)=ALF (I,K)*SZZ(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)+3.*ALF(I,K)*SYY(I,KP,L,JV))
       SYY(I,KP,L,JV)=ALF4(I,K)*SYY(I,KP,L,JV)
       SSX (I,KP,L,JV)=SSX (I,KP,L,JV)-FX (I,K,JV)
       SZ (I,KP,L,JV)=SZ (I,KP,L,JV)-FZ (I,K,JV)
       SSXX(I,KP,L,JV)=SSXX(I,KP,L,JV)-FXX(I,K,JV)
       SSXZ(I,KP,L,JV)=SSXZ(I,KP,L,JV)-FXZ(I,K,JV)
       SZZ(I,KP,L,JV)=SZZ(I,KP,L,JV)-FZZ(I,K,JV)
       SSXY(I,KP,L,JV)=ALF1Q(I,K)*SSXY(I,KP,L,JV)
       SYZ(I,KP,L,JV)=ALF1Q(I,K)*SYZ(I,KP,L,JV)
  !
     ELSE
  !
       F0 (I,K,JV)=ALF (I,K)* ( S0(I,K,L,JV)+ALF1(I,K)* &
             ( SY(I,K,L,JV)+ALF2(I,K)*SYY(I,K,L,JV) ) )
       FY (I,K,JV)=ALFQ(I,K)* &
             (SY(I,K,L,JV)+3.*ALF1(I,K)*SYY(I,K,L,JV))
       FYY(I,K,JV)=ALF3(I,K)*SYY(I,K,L,JV)
  FX (I,K,JV)=ALF (I,K)*(SSX(I,K,L,JV)+ALF1(I,K)*SSXY(I,K,L,JV))
  FZ (I,K,JV)=ALF (I,K)*(SZ(I,K,L,JV)+ALF1(I,K)*SYZ(I,K,L,JV))
       FXY(I,K,JV)=ALFQ(I,K)*SSXY(I,K,L,JV)
       FYZ(I,K,JV)=ALFQ(I,K)*SYZ(I,K,L,JV)
       FXX(I,K,JV)=ALF (I,K)*SSXX(I,K,L,JV)
       FXZ(I,K,JV)=ALF (I,K)*SSXZ(I,K,L,JV)
       FZZ(I,K,JV)=ALF (I,K)*SZZ(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)-3.*ALF(I,K)*SYY(I,K,L,JV))
       SYY(I,K,L,JV)=ALF4(I,K)*SYY(I,K,L,JV)
       SSX (I,K,L,JV)=SSX (I,K,L,JV)-FX (I,K,JV)
       SZ (I,K,L,JV)=SZ (I,K,L,JV)-FZ (I,K,JV)
       SSXX(I,K,L,JV)=SSXX(I,K,L,JV)-FXX(I,K,JV)
       SSXZ(I,K,L,JV)=SSXZ(I,K,L,JV)-FXZ(I,K,JV)
       SZZ(I,K,L,JV)=SZZ(I,K,L,JV)-FZZ(I,K,JV)
       SSXY(I,K,L,JV)=ALF1Q(I,K)*SSXY(I,K,L,JV)
       SYZ(I,K,L,JV)=ALF1Q(I,K)*SYZ(I,K,L,JV)
  !
     ENDIF
  !
  END DO
  END DO
  END DO
  ! print*,'ap ADVYP 31'
  !
  !  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
  !
     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)
     ALF2(I,K)=ALF1(I,K)-ALF(I,K)
     ALF3(I,K)=ALF1(I,K)*ALF(I,K)
  !
  END DO
  END DO
  ! print*,'ap ADVYP 32'
  !
  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)
   SYY(I,K,L,JV)=ALFQ(I,K)*FYY(I,K,JV)+ALF1Q(I,K)*SYY(I,K,L,JV) &
         +5.*( ALF3(I,K)*(FY(I,K,JV)-SY(I,K,L,JV))+ALF2(I,K)*TEMPTM )
     SY (I,K,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*SY(I,K,L,JV) &
           +3.*TEMPTM
   SSXY(I,K,L,JV)=ALF (I,K)*FXY(I,K,JV)+ALF1(I,K)*SSXY(I,K,L,JV) &
         +3.*(ALF1(I,K)*FX (I,K,JV)-ALF (I,K)*SSX (I,K,L,JV))
   SYZ(I,K,L,JV)=ALF (I,K)*FYZ(I,K,JV)+ALF1(I,K)*SYZ(I,K,L,JV) &
         +3.*(ALF1(I,K)*FZ (I,K,JV)-ALF (I,K)*SZ (I,K,L,JV))
     SSX (I,K,L,JV)=SSX (I,K,L,JV)+FX (I,K,JV)
     SZ (I,K,L,JV)=SZ (I,K,L,JV)+FZ (I,K,JV)
     SSXX(I,K,L,JV)=SSXX(I,K,L,JV)+FXX(I,K,JV)
     SSXZ(I,K,L,JV)=SSXZ(I,K,L,JV)+FXZ(I,K,JV)
     SZZ(I,K,L,JV)=SZZ(I,K,L,JV)+FZZ(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)
   SYY(I,KP,L,JV)=ALFQ(I,K)*FYY(I,K,JV)+ALF1Q(I,K)*SYY(I,KP,L,JV) &
         +5.*( ALF3(I,K)*(SY(I,KP,L,JV)-FY(I,K,JV))-ALF2(I,K)*TEMPTM )
     SY (I,KP,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*SY(I,KP,L,JV) &
           +3.*TEMPTM
   SSXY(I,KP,L,JV)=ALF(I,K)*FXY(I,K,JV)+ALF1(I,K)*SSXY(I,KP,L,JV) &
         +3.*(ALF(I,K)*SSX(I,KP,L,JV)-ALF1(I,K)*FX(I,K,JV))
     SYZ(I,KP,L,JV)=ALF(I,K)*FYZ(I,K,JV)+ALF1(I,K)*SYZ(I,KP,L,JV) &
           +3.*(ALF(I,K)*SZ(I,KP,L,JV)-ALF1(I,K)*FZ(I,K,JV))
     SSX (I,KP,L,JV)=SSX (I,KP,L,JV)+FX (I,K,JV)
     SZ (I,KP,L,JV)=SZ (I,KP,L,JV)+FZ (I,K,JV)
     SSXX(I,KP,L,JV)=SSXX(I,KP,L,JV)+FXX(I,K,JV)
     SSXZ(I,KP,L,JV)=SSXZ(I,KP,L,JV)+FXZ(I,K,JV)
     SZZ(I,KP,L,JV)=SZZ(I,KP,L,JV)+FZZ(I,K,JV)
  !
     ENDIF
  !
  END DO
  END DO
  END DO
  ! print*,'ap ADVYP 33'
  !
  !  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)
     ALF2(I,K)=ALF1(I,K)-ALF(I,K)
     ALF3(I,K)=ALF(I,K)*ALFQ(I,K)
     ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K)
  !
  END DO
  ! print*,'ap ADVYP 41'
  !
  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)+ALF2(I,K)*SYY(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)-3.*ALF(I,K)*SYY(I,K,L,JV))
       SYY(I,K,L,JV)=ALF4 (I,K)*SYY(I,K,L,JV)
  SSX (I,K,L,JV)=ALF1(I,K)*(SSX(I,K,L,JV)-ALF(I,K)*SSXY(I,K,L,JV))
  SZ (I,K,L,JV)=ALF1(I,K)*(SZ(I,K,L,JV)-ALF(I,K)*SYZ(I,K,L,JV))
       SSXX(I,K,L,JV)=ALF1 (I,K)*SSXX(I,K,L,JV)
       SSXZ(I,K,L,JV)=ALF1 (I,K)*SSXZ(I,K,L,JV)
       SZZ(I,K,L,JV)=ALF1 (I,K)*SZZ(I,K,L,JV)
       SSXY(I,K,L,JV)=ALF1Q(I,K)*SSXY(I,K,L,JV)
       SYZ(I,K,L,JV)=ALF1Q(I,K)*SYZ(I,K,L,JV)
     ENDIF
  !
  END DO
  END DO
  ! print*,'ap ADVYP 42'
  !
  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
  ! print*,'ap ADVYP 43'
  !
  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
  !
     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)
     ALF2(I,K)=ALF1(I,K)-ALF(I,K)
     ALF3(I,K)=ALF1(I,K)*ALF(I,K)
  !
  END DO
  ! print*,'ap ADVYP 45'
  !
  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)
     SYY(I,K,L,JV)=ALF1Q(I,K)*SYY(I,K,L,JV) &
           +5.*(-ALF3 (I,K)*SY (I,K,L,JV)+ALF2(I,K)*TEMPTM )
     SY (I,K,L,JV)=ALF1(I,K)*SY (I,K,L,JV)+3.*TEMPTM
  SSXY(I,K,L,JV)=ALF1(I,K)*SSXY(I,K,L,JV)-3.*ALF(I,K)*SSX(I,K,L,JV)
  SYZ(I,K,L,JV)=ALF1(I,K)*SYZ(I,K,L,JV)-3.*ALF(I,K)*SZ(I,K,L,JV)
  !
     ENDIF
  !
  END DO
  END DO
  ! print*,'ap ADVYP 46'
  !
  END DO

  !--------------------------------------------------
  ! bouclage cyclique horizontal .

  DO l = 1,llm
     DO jv = 1,ntra
        DO j = 1,jjp1
           SM(iip1,j,l) = SM(1,j,l)
           S0(iip1,j,l,jv) = S0(1,j,l,jv)
           SSX(iip1,j,l,jv) = SSX(1,j,l,jv)
           SY(iip1,j,l,jv) = SY(1,j,l,jv)
           SZ(iip1,j,l,jv) = SZ(1,j,l,jv)
        END DO
     END DO
  END DO

  ! -------------------------------------------------------------------
  ! *** Test  negativite:

   ! DO jv = 1,ntra
   !  DO l = 1,llm
   !    DO j = 1,jjp1
   !      DO i = 1,iip1
   !         IF (s0( i,j,l,jv ).lt.0.) THEN
   !            PRINT*, '------ S0 < 0 en FIN ADVYP ---'
   !            PRINT*, 'S0(',i,j,l,jv,')=', S0(i,j,l,jv)
  !c                 STOP
   !         ENDIF
   !      ENDDO
   !    ENDDO
   !  ENDDO
   ! ENDDO


  ! -------------------------------------------------------------------
  ! *** Test : diag de la qtite totale de traceur dans
   !       l'atmosphere avant l'advection en Y

   DO l = 1,llm
     DO j = 1,jjp1
       DO i = 1,iim
          sqf = sqf + S0(i,j,l,ntra)
       END DO
     END DO
   END DO
  PRINT*,'---------- DIAG DANS ADVY - SORTIE --------'
  PRINT*,'sqf=',sqf
  ! print*,'ap ADVYP fin'

  !-----------------------------------------------------------------
  !
  RETURN
END SUBROUTINE ADVYP