
! $Header$

      SUBROUTINE ADVZP(LIMIT,DTZ,W,SM,S0,SSX,SY,SZ
     .                 ,SSXX,SSXY,SSXZ,SYY,SYZ,SZZ,ntra )

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                 C
C  second-order moments (SOM) advection of tracer in Z direction  C
C                                                                 C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                 C
C  Source : Pascal Simon ( Meteo, CNRM )                          C
C  Adaptation : A.A. (LGGE)                                       C
C  Derniere Modif : 19/11/95 LAST                                 C
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
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Rem : Probleme aux poles il faut reecrire ce cas specifique
C        Attention au sens de l'indexation
C

C
C  parametres principaux du modele
C
      include "dimensions.h"
      include "paramet.h"
      include "comgeom.h"
C
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
c
        INTEGER lon,lat,niv
        INTEGER i,j,jv,k,kp,l,lp
        INTEGER ntra
c        PARAMETER (ntra = 1)
c
        REAL dtz
        REAL w ( iip1,jjp1,llm )
c
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 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)
C
C  Local :
C  -------
C
C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
C  mass fluxes in kg
C  declaration :
C
      REAL WGRI(iip1,jjp1,0:llm)

C Rem : UGRI et VGRI ne sont pas utilises dans
C  cette subroutine ( advection en z uniquement )
C  Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
C         attention a celui de WGRI
C
C  the moments F are similarly defined and used as temporary
C  storage for portions of the grid boxes in transit
C
C  the moments Fij are used as temporary storage for
C  portions of the grid boxes in transit at the current level
C
C  work arrays
C
C
      REAL F0(iim,llm,ntra),FM(iim,llm)
      REAL FX(iim,llm,ntra),FY(iim,llm,ntra)
      REAL FZ(iim,llm,ntra)
      REAL FXX(iim,llm,ntra),FXY(iim,llm,ntra)
      REAL FXZ(iim,llm,ntra),FYY(iim,llm,ntra)
      REAL FYZ(iim,llm,ntra),FZZ(iim,llm,ntra)
      REAL S00(ntra)
      REAL SM0             ! Just temporal variable
C
C  work arrays
C
      REAL ALF(iim),ALF1(iim)
      REAL ALFQ(iim),ALF1Q(iim)
      REAL ALF2(iim),ALF3(iim)
      REAL ALF4(iim)
      REAL TEMPTM          ! Just temporal variable
      REAL SLPMAX,S1MAX,S1NEW,S2NEW
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         !       tab. S et VGRI 
                    
c-----------------------------------------------------------------
C *** Test : diag de la qtite totale de traceur dans
C            l'atmosphere avant l'advection en Y
c 
      sqi = 0.
      sqf = 0.
c
      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 ADVZP - ENTREE --------'
      PRINT*,'sqi=',sqi

c-----------------------------------------------------------------
C  Interface : adaptation nouveau modele
C  -------------------------------------
C
C  Conversion des flux de masses en kg

      DO l = 1,llm
         DO j = 1,jjp1
            DO i = 1,iip1  
            wgri (i,j,llm+1-l) = w (i,j,l)  
      END DO
      END DO
      END DO
      do j=1,jjp1
         do i=1,iip1
            wgri(i,j,0)=0.
         enddo
      enddo
c
cAA rem : Je ne suis pas sur du signe  
cAA       Je ne suis pas sur pour le 0:llm
c
c-----------------------------------------------------------------
C---------------------- START HERE -------------------------------
C
C  boucle sur les latitudes
C
      DO K=1,LAT
C
C  place limits on appropriate moments before transport
C      (if flux-limiting is to be applied)
C
      IF(.NOT.LIMIT) GO TO 101
C
      DO JV=1,NTRA
      DO L=1,NIV
         DO I=1,LON
            IF(S0(I,K,L,JV)>0.) THEN
              SLPMAX=S0(I,K,L,JV)
              S1MAX =1.5*SLPMAX
              S1NEW =AMIN1(S1MAX,AMAX1(-S1MAX,SZ(I,K,L,JV)))
              S2NEW =AMIN1( 2.*SLPMAX-ABS(S1NEW)/3. ,
     +                     AMAX1(ABS(S1NEW)-SLPMAX,SZZ(I,K,L,JV)) )
              SZ (I,K,L,JV)=S1NEW
              SZZ(I,K,L,JV)=S2NEW
              SSXZ(I,K,L,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,SSXZ(I,K,L,JV)))
              SYZ(I,K,L,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,SYZ(I,K,L,JV)))
            ELSE
              SZ (I,K,L,JV)=0.
              SZZ(I,K,L,JV)=0.
              SSXZ(I,K,L,JV)=0.
              SYZ(I,K,L,JV)=0.
            ENDIF
      END DO
      END DO
      END DO
C
 101  CONTINUE
C
C  boucle sur les niveaux intercouches de 1 a NIV-1
C   (flux nul au sommet L=0 et a la base L=NIV)
C
C  calculate flux and moments between adjacent boxes
C     (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
C  1- create temporary moments/masses for partial boxes in transit
C  2- reajusts moments remaining in the box
C
      DO L=1,NIV-1
      LP=L+1
C
      DO I=1,LON
C
         IF(WGRI(I,K,L)<0.) THEN
           FM(I,L)=-WGRI(I,K,L)*DTZ
           ALF(I)=FM(I,L)/SM(I,K,LP)
           SM(I,K,LP)=SM(I,K,LP)-FM(I,L)
         ELSE
           FM(I,L)=WGRI(I,K,L)*DTZ
           ALF(I)=FM(I,L)/SM(I,K,L)
           SM(I,K,L)=SM(I,K,L)-FM(I,L)
         ENDIF
C
         ALFQ (I)=ALF(I)*ALF(I)
         ALF1 (I)=1.-ALF(I)
         ALF1Q(I)=ALF1(I)*ALF1(I)
         ALF2 (I)=ALF1(I)-ALF(I)
         ALF3 (I)=ALF(I)*ALFQ(I)
         ALF4 (I)=ALF1(I)*ALF1Q(I)
C
      END DO
C
      DO JV=1,NTRA
      DO I=1,LON
C
         IF(WGRI(I,K,L)<0.) THEN
C
           F0 (I,L,JV)=ALF (I)* ( S0(I,K,LP,JV)-ALF1(I)*
     +          ( SZ(I,K,LP,JV)-ALF2(I)*SZZ(I,K,LP,JV) ) )
           FZ (I,L,JV)=ALFQ(I)*(SZ(I,K,LP,JV)-3.*ALF1(I)*SZZ(I,K,LP,JV))
           FZZ(I,L,JV)=ALF3(I)*SZZ(I,K,LP,JV)
           FXZ(I,L,JV)=ALFQ(I)*SSXZ(I,K,LP,JV)
           FYZ(I,L,JV)=ALFQ(I)*SYZ(I,K,LP,JV)
           FX (I,L,JV)=ALF (I)*(SSX(I,K,LP,JV)-ALF1(I)*SSXZ(I,K,LP,JV))
           FY (I,L,JV)=ALF (I)*(SY(I,K,LP,JV)-ALF1(I)*SYZ(I,K,LP,JV))
           FXX(I,L,JV)=ALF (I)*SSXX(I,K,LP,JV)
           FXY(I,L,JV)=ALF (I)*SSXY(I,K,LP,JV)
           FYY(I,L,JV)=ALF (I)*SYY(I,K,LP,JV)
C
           S0 (I,K,LP,JV)=S0 (I,K,LP,JV)-F0 (I,L,JV)
           SZ (I,K,LP,JV)=ALF1Q(I)
     +                   *(SZ(I,K,LP,JV)+3.*ALF(I)*SZZ(I,K,LP,JV))
           SZZ(I,K,LP,JV)=ALF4 (I)*SZZ(I,K,LP,JV)
           SSXZ(I,K,LP,JV)=ALF1Q(I)*SSXZ(I,K,LP,JV)
           SYZ(I,K,LP,JV)=ALF1Q(I)*SYZ(I,K,LP,JV)
           SSX (I,K,LP,JV)=SSX (I,K,LP,JV)-FX (I,L,JV)
           SY (I,K,LP,JV)=SY (I,K,LP,JV)-FY (I,L,JV)
           SSXX(I,K,LP,JV)=SSXX(I,K,LP,JV)-FXX(I,L,JV)
           SSXY(I,K,LP,JV)=SSXY(I,K,LP,JV)-FXY(I,L,JV)
           SYY(I,K,LP,JV)=SYY(I,K,LP,JV)-FYY(I,L,JV)
C
         ELSE
C
           F0 (I,L,JV)=ALF (I)*(S0(I,K,L,JV)
     +           +ALF1(I) * (SZ(I,K,L,JV)+ALF2(I)*SZZ(I,K,L,JV)) )
           FZ (I,L,JV)=ALFQ(I)*(SZ(I,K,L,JV)+3.*ALF1(I)*SZZ(I,K,L,JV))
           FZZ(I,L,JV)=ALF3(I)*SZZ(I,K,L,JV)
           FXZ(I,L,JV)=ALFQ(I)*SSXZ(I,K,L,JV)
           FYZ(I,L,JV)=ALFQ(I)*SYZ(I,K,L,JV)
           FX (I,L,JV)=ALF (I)*(SSX(I,K,L,JV)+ALF1(I)*SSXZ(I,K,L,JV))
           FY (I,L,JV)=ALF (I)*(SY(I,K,L,JV)+ALF1(I)*SYZ(I,K,L,JV))
           FXX(I,L,JV)=ALF (I)*SSXX(I,K,L,JV)
           FXY(I,L,JV)=ALF (I)*SSXY(I,K,L,JV)
           FYY(I,L,JV)=ALF (I)*SYY(I,K,L,JV)
C
           S0 (I,K,L,JV)=S0 (I,K,L,JV)-F0(I,L,JV)
           SZ (I,K,L,JV)=ALF1Q(I)*(SZ(I,K,L,JV)-3.*ALF(I)*SZZ(I,K,L,JV))
           SZZ(I,K,L,JV)=ALF4 (I)*SZZ(I,K,L,JV)
           SSXZ(I,K,L,JV)=ALF1Q(I)*SSXZ(I,K,L,JV)
           SYZ(I,K,L,JV)=ALF1Q(I)*SYZ(I,K,L,JV)
           SSX (I,K,L,JV)=SSX (I,K,L,JV)-FX (I,L,JV)
           SY (I,K,L,JV)=SY (I,K,L,JV)-FY (I,L,JV)
           SSXX(I,K,L,JV)=SSXX(I,K,L,JV)-FXX(I,L,JV)
           SSXY(I,K,L,JV)=SSXY(I,K,L,JV)-FXY(I,L,JV)
           SYY(I,K,L,JV)=SYY(I,K,L,JV)-FYY(I,L,JV)
C
         ENDIF
C
      END DO
      END DO
C
      END DO
C
C  puts the temporary moments Fi into appropriate neighboring boxes
C
      DO L=1,NIV-1
      LP=L+1
C
      DO I=1,LON
C
         IF(WGRI(I,K,L)<0.) THEN
           SM(I,K,L)=SM(I,K,L)+FM(I,L)
           ALF(I)=FM(I,L)/SM(I,K,L)
         ELSE
           SM(I,K,LP)=SM(I,K,LP)+FM(I,L)
           ALF(I)=FM(I,L)/SM(I,K,LP)
         ENDIF
C
         ALF1(I)=1.-ALF(I)
         ALFQ(I)=ALF(I)*ALF(I)
         ALF1Q(I)=ALF1(I)*ALF1(I)
         ALF2(I)=ALF(I)*ALF1(I)
         ALF3(I)=ALF1(I)-ALF(I)
C
      END DO
C
      DO JV=1,NTRA
      DO I=1,LON
C
         IF(WGRI(I,K,L)<0.) THEN
C
           TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)
           S0 (I,K,L,JV)=S0(I,K,L,JV)+F0(I,L,JV)
           SZZ(I,K,L,JV)=ALFQ(I)*FZZ(I,L,JV)+ALF1Q(I)*SZZ(I,K,L,JV)
     +        +5.*( ALF2(I)*(FZ(I,L,JV)-SZ(I,K,L,JV))+ALF3(I)*TEMPTM )
           SZ (I,K,L,JV)=ALF (I)*FZ (I,L,JV)+ALF1 (I)*SZ (I,K,L,JV)
     +                  +3.*TEMPTM
           SSXZ(I,K,L,JV)=ALF (I)*FXZ(I,L,JV)+ALF1 (I)*SSXZ(I,K,L,JV)
     +              +3.*(ALF1(I)*FX (I,L,JV)-ALF  (I)*SSX (I,K,L,JV))
           SYZ(I,K,L,JV)=ALF (I)*FYZ(I,L,JV)+ALF1 (I)*SYZ(I,K,L,JV)
     +              +3.*(ALF1(I)*FY (I,L,JV)-ALF  (I)*SY (I,K,L,JV))
           SSX (I,K,L,JV)=SSX (I,K,L,JV)+FX (I,L,JV)
           SY (I,K,L,JV)=SY (I,K,L,JV)+FY (I,L,JV)
           SSXX(I,K,L,JV)=SSXX(I,K,L,JV)+FXX(I,L,JV)
           SSXY(I,K,L,JV)=SSXY(I,K,L,JV)+FXY(I,L,JV)
           SYY(I,K,L,JV)=SYY(I,K,L,JV)+FYY(I,L,JV)
C
         ELSE
C
           TEMPTM=ALF(I)*S0(I,K,LP,JV)-ALF1(I)*F0(I,L,JV)
           S0 (I,K,LP,JV)=S0(I,K,LP,JV)+F0(I,L,JV)
           SZZ(I,K,LP,JV)=ALFQ(I)*FZZ(I,L,JV)+ALF1Q(I)*SZZ(I,K,LP,JV)
     +        +5.*( ALF2(I)*(SZ(I,K,LP,JV)-FZ(I,L,JV))-ALF3(I)*TEMPTM )
           SZ (I,K,LP,JV)=ALF (I)*FZ(I,L,JV)+ALF1(I)*SZ(I,K,LP,JV)
     +                   +3.*TEMPTM
           SSXZ(I,K,LP,JV)=ALF(I)*FXZ(I,L,JV)+ALF1(I)*SSXZ(I,K,LP,JV)
     +                   +3.*(ALF(I)*SSX(I,K,LP,JV)-ALF1(I)*FX(I,L,JV))
           SYZ(I,K,LP,JV)=ALF(I)*FYZ(I,L,JV)+ALF1(I)*SYZ(I,K,LP,JV)
     +                   +3.*(ALF(I)*SY(I,K,LP,JV)-ALF1(I)*FY(I,L,JV))
           SSX (I,K,LP,JV)=SSX (I,K,LP,JV)+FX (I,L,JV)
           SY (I,K,LP,JV)=SY (I,K,LP,JV)+FY (I,L,JV)
           SSXX(I,K,LP,JV)=SSXX(I,K,LP,JV)+FXX(I,L,JV)
           SSXY(I,K,LP,JV)=SSXY(I,K,LP,JV)+FXY(I,L,JV)
           SYY(I,K,LP,JV)=SYY(I,K,LP,JV)+FYY(I,L,JV)
C
         ENDIF
C
      END DO
      END DO
C
      END DO
C
C  fin de la boucle principale sur les latitudes
C
      END DO
C
      DO l = 1,llm
      DO j = 1,jjp1
          SM(iip1,j,l) = SM(1,j,l)
          S0(iip1,j,l,ntra) = S0(1,j,l,ntra)
          SSX(iip1,j,l,ntra) = SSX(1,j,l,ntra)
          SY(iip1,j,l,ntra) = SY(1,j,l,ntra)
          SZ(iip1,j,l,ntra) = SZ(1,j,l,ntra)
      ENDDO
      ENDDO
c                                                                                C-------------------------------------------------------------
C *** Test : diag de la qqtite totale de tarceur
C            dans l'atmosphere avant l'advection en z
       DO l = 1,llm
       DO j = 1,jjp1
       DO i = 1,iim
          sqf = sqf + S0(i,j,l,ntra)
       ENDDO
       ENDDO
       ENDDO
       PRINT*,'-------- DIAG DANS ADVZ - SORTIE ---------'
       PRINT*,'sqf=', sqf

      RETURN
      END
