! $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 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 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
C
 11   CONTINUE
C
C  le flux a travers le pole Nord est traite separement
C
      SM0=0.
      DO JV=1,NTRA
         S00(JV)=0.
      END DO
C
      DO I=1,LON
C
         IF(VGRI(I,0,L)<=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
      END DO
C
      DO JV=1,NTRA
      DO I=1,LON
C
         IF(VGRI(I,0,L)<=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
      END DO
      END DO
C
      DO I=1,LON
         IF(VGRI(I,0,L)>0.) THEN
           FM(I,0)=VGRI(I,0,L)*DTY
           ALF(I,0)=FM(I,0)/SM0
         ENDIF
      END DO
C
      DO JV=1,NTRA
      DO I=1,LON
         IF(VGRI(I,0,L)>0.) THEN
           F0(I,0,JV)=ALF(I,0)*S00(JV)
         ENDIF
      END DO
      END DO
C
C  puts the temporary moments Fi into appropriate neighboring boxes
C
      DO I=1,LON
C
         IF(VGRI(I,0,L)>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
      END DO
C
      DO JV=1,NTRA
      DO I=1,LON
C
         IF(VGRI(I,0,L)>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
      END DO
      END DO
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 K=1,LAT-1
      KP=K+1
      DO I=1,LON
C
         IF(VGRI(I,K,L)<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
      END DO
      END DO
C
      DO JV=1,NTRA
      DO K=1,LAT-1
      KP=K+1
      DO I=1,LON
C
         IF(VGRI(I,K,L)<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
      END DO
      END DO
      END DO
C
C  puts the temporary moments Fi into appropriate neighboring boxes
C
      DO K=1,LAT-1
      KP=K+1
      DO I=1,LON
C
         IF(VGRI(I,K,L)<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
      END DO
      END DO
C
      DO JV=1,NTRA
      DO K=1,LAT-1
      KP=K+1
      DO I=1,LON
C
         IF(VGRI(I,K,L)<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
      END DO
      END DO
      END DO
C
C  traitement special pour le pole Sud (idem pole Nord)
C
      K=LAT
C
      SM0=0.
      DO JV=1,NTRA
         S00(JV)=0.
      END DO
C
      DO I=1,LON
C
         IF(VGRI(I,K,L)>=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
      END DO
C
      DO JV=1,NTRA
      DO I=1,LON
C
         IF(VGRI(I,K,L)>=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
      END DO
      END DO
C
      DO I=1,LON
         IF(VGRI(I,K,L)<0.) THEN
           FM(I,K)=-VGRI(I,K,L)*DTY
           ALF(I,K)=FM(I,K)/SM0
         ENDIF
      END DO
C
      DO JV=1,NTRA
      DO I=1,LON
         IF(VGRI(I,K,L)<0.) THEN
           F0(I,K,JV)=ALF(I,K)*S00(JV)
         ENDIF
      END DO
      END DO
C
C  puts the temporary moments Fi into appropriate neighboring boxes
C
      DO I=1,LON
C
         IF(VGRI(I,K,L)<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
      END DO
C
      DO JV=1,NTRA
      DO I=1,LON
C
         IF(VGRI(I,K,L)<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
      END DO
      END DO
C
      END DO
C
      RETURN
      END