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

  USE comgeom_mod_h
  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
USE paramet_mod_h
IMPLICIT NONE

  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  !                                                                 C
  !  second-order moments (SOM) advection of tracer in Z direction  C
  !                                                                 C
  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  !                                                                 C
  !  Source : Pascal Simon ( Meteo, CNRM )                          C
  !  Adaptation : A.A. (LGGE)                                       C
  !  Derniere Modif : 19/11/95 LAST                                 C
  !                                                                 C
  !  sont les arguments d'entree pour le s-pg                       C
  !                                                                 C
  !  argument de sortie du s-pg                                     C
  !                                                                 C
  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  !
  ! Rem : Probleme aux poles il faut reecrire ce cas specifique
  !    Attention au sens de l'indexation
  !

  !
  !  parametres principaux du modele
  !


  !
  !  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,lp
    INTEGER :: ntra
     ! PARAMETER (ntra = 1)
  !
    REAL :: dtz
    REAL :: w ( iip1,jjp1,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 :: WGRI(iip1,jjp1,0:llm)

  ! Rem : UGRI et VGRI ne sont pas utilises dans
  !  cette subroutine ( advection en z uniquement )
  !  Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
      ! attention a celui de WGRI
  !
  !  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,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
  !
  !  work arrays
  !
  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
  !
  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

  !-----------------------------------------------------------------
  ! *** 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 ADVZP - ENTREE --------'
  PRINT*,'sqi=',sqi

  !-----------------------------------------------------------------
  !  Interface : adaptation nouveau modele
  !  -------------------------------------
  !
  !  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
  !
  !AA rem : Je ne suis pas sur du signe
  !AA       Je ne suis pas sur pour le 0:llm
  !
  !-----------------------------------------------------------------
  !---------------------- START HERE -------------------------------
  !
  !  boucle sur les latitudes
  !
  DO K=1,LAT
  !
  !  place limits on appropriate moments before transport
  !  (if flux-limiting is to be applied)
  !
  IF(.NOT.LIMIT) GO TO 101
  !
  DO JV=1,NTRA
  DO L=1,NIV
     DO I=1,LON
        IF(S0(I,K,L,JV).GT.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
  !
 101   CONTINUE
  !
  !  boucle sur les niveaux intercouches de 1 a NIV-1
  !   (flux nul au sommet L=0 et a la base L=NIV)
  !
  !  calculate flux and moments between adjacent boxes
  ! (flux from LP to L if WGRI(L).lt.0, from L to LP if WGRI(L).gt.0)
  !  1- create temporary moments/masses for partial boxes in transit
  !  2- reajusts moments remaining in the box
  !
  DO L=1,NIV-1
  LP=L+1
  !
  DO I=1,LON
  !
     IF(WGRI(I,K,L).LT.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
  !
     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)
  !
  END DO
  !
  DO JV=1,NTRA
  DO I=1,LON
  !
     IF(WGRI(I,K,L).LT.0.) THEN
  !
       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)
  !
       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)
  !
     ELSE
  !
       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)
  !
       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)
  !
     ENDIF
  !
  END DO
  END DO
  !
  END DO
  !
  !  puts the temporary moments Fi into appropriate neighboring boxes
  !
  DO L=1,NIV-1
  LP=L+1
  !
  DO I=1,LON
  !
     IF(WGRI(I,K,L).LT.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
  !
     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)
  !
  END DO
  !
  DO JV=1,NTRA
  DO I=1,LON
  !
     IF(WGRI(I,K,L).LT.0.) THEN
  !
       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)
  !
     ELSE
  !
       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)
  !
     ENDIF
  !
  END DO
  END DO
  !
  END DO
  !
  !  fin de la boucle principale sur les latitudes
  !
  END DO
  !
  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-------------------------------------------------------------
  ! *** Test : diag de la qqtite totale de tarceur
         ! 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 SUBROUTINE ADVZP
