  SUBROUTINE sub_advz_quickupwind(xe,c)
 
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sub_advz_quickupwind  ***
!&E
!&E ** Description : A limiter switch to the QUICK or upwind scheme.
!&E                  No estimate above small depths ( < hm)
!&E                  Used if key_tssub_wquickupwind
!&E ** Purpose : compute vertical fluxes of tracers other than temperature and salinity
!&E              using TVD scheme (quick and upwind decentred) with all-or-nothing limitator
!&E
!&E ** Description : For very low water heights ( < hm), no vertical fluxes are computed
!&E                  because the scheme is degenerated vertically and vertical fluxes are useless.
!&E                  This routine is called only if the key_tssub_wquickupwind CPP key
!&E                  is activated (see tssub_eq_wquick.h).
!&E
!&E
!&E ** Called by : adv
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xe,h0,fzeau,c(cv_wat)
!&E
!&E ** Modified variables : fz
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2004-10-04
!&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
!&E       !  2009-04  (A. Mano, V. Garnier) MPI2D
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2016-04  (B.Thouvenin) applied to all substances (dissolved and particulate, include sand)
!&E
!&E---------------------------------------------------------------------
!! * Modules used

   USE parameters
   USE comvars2d,   ONLY : ig,id,h0,hm
   USE comvars3d,   ONLY : sig,dsigw,qmaxz
   USE comsubstance,ONLY : nv_adv,nv_grav,fz
   USE comfluxdo,   ONLY : fzeau

   IMPLICIT NONE

!! * Arguments
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in)             :: xe
   REAL(kind=rsh),DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: c

!! * Local declarations
   LOGICAL                               :: decent
   INTEGER                               :: i,j,k,iv
   REAL(kind=rsh)                        :: unsqmaxz
   REAL(kind=rsh)                        :: yr,h3,qot,dyd,dyg,ya,yb,yc,x1,x2,x3,xx
   REAL(kind=rsh),DIMENSION(kmax)        :: ap,bp,cp,am,bm,cm

!!----------------------------------------------------------------------
!! * Executable part

   unsqmaxz=1.0_rsh/qmaxz

!$OMP DO schedule(runtime),private(i,j,k,iv,h3 &
!$OMP& ,ya,yb,yc,yr,qot,decent,dyd,dyg,x1,x2,x3,xx,ap,bp,cp,am,bm,cm)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)

       h3=h0(i,j)+xe(i,j)
       IF (h3 > hm) THEN

! calcul des coef du schema QUICK pour w positif
!-----------------------------------------------

         DO k=2,kmax-1

           x1=sig(k-1)
           x2=sig(k)
           x3=sig(k+1)

           xx=(x1-(x2+x3)/2.0_rsh)*2.0_rsh/(x3-x2)
           ap(k)=1.0_rsh/(1.0_rsh-xx**2)
           bp(k)=(xx-xx**2)/(2.0_rsh*(1.0_rsh-xx**2))
           cp(k)=-(xx+xx**2)/(2.0_rsh*(1.0_rsh-xx**2))
         END DO


! cas w negatif
!-------------------

         DO k=1,kmax-2

           x1=sig(k)
           x2=sig(k+1)
           x3=sig(k+2)

           xx=(x1-(x2+x3)/2.0_rsh)*2.0_rsh/(x2-x3)
           am(k)=1.0_rsh/(1.0_rsh-xx**2)
           bm(k)=-(xx+xx**2)/(2.0_rsh*(1.0_rsh-xx**2))
           cm(k)=(xx-xx**2)/(2.0_rsh*(1.0_rsh-xx**2))
         END DO

! pour faire simple : decentre amont pres des limites
! ---------------------------------------------------

         fz(1:nv_adv,0:kmax,i,j)=0.0_rsh

         k=1
         IF (fzeau(k,i,j) > 0.0_rlg) THEN
           fz(nv_grav+1:nv_adv,k,i,j)=fzeau(k,i,j)*c(nv_grav+1:nv_adv,k,i,j)
         ELSE
           fz(nv_grav+1:nv_adv,k,i,j)=fzeau(k,i,j)*c(nv_grav+1:nv_adv,k+1,i,j)
         END IF

         k=kmax-1
         IF (fzeau(k,i,j) > 0.0_rlg) THEN
           fz(nv_grav+1:nv_adv,k,i,j)=fzeau(k,i,j)*c(nv_grav+1:nv_adv,k,i,j)
         ELSE
           fz(nv_grav+1:nv_adv,k,i,j)=fzeau(k,i,j)*c(nv_grav+1:nv_adv,k+1,i,j)
         END IF

         DO k=2,kmax-2
           IF (fzeau(k,i,j) <= 0.0_rlg) THEN
             DO iv=nv_grav+1,nv_adv
               yb=c(iv,k+1,i,j)
               yc=c(iv,k,i,j)
               IF (k <= kmax-2) THEN
                 ya=c(iv,k+2,i,j)
               ELSE
                 ya=yb
               END IF
               dyd=(yb-yc)/dsigw(k)
               dyg=(yb-ya)/dsigw(k+1)

               decent=(dyd*dyg) >= 0.0_rsh
               IF (.NOT.decent) THEN
                 qot=ABS(dyd/dyg)
                 decent=.TRUE.
                 IF (qot <= qmaxz .AND. qot >= unsqmaxz) decent=.FALSE.
               END IF

               IF (decent) THEN
                 fz(iv,k,i,j)=fzeau(k,i,j)*yb
               ELSE
                 yr=cm(k)*yc+bm(k)*yb+am(k)*ya
                 fz(iv,k,i,j)=fzeau(k,i,j)*yr
               END IF
             END DO

           ELSE ! test sur la vitesse cas w >0

             DO iv=nv_grav+1,nv_adv
               yb=c(iv,k,i,j)
               yc=c(iv,k+1,i,j)
               IF (k >= 2) THEN
                 ya=c(iv,k-1,i,j)
               ELSE
                 ya=yb
               ENDIF

               dyd=(yc-yb)/dsigw(k)
               dyg=(yb-ya)/dsigw(k-1)

               decent=(dyd*dyg) <= 0.0_rsh
               IF (.NOT.decent) THEN
                 qot=ABS(dyd/dyg)
                 decent=.TRUE.
                 IF (qot <= qmaxz .AND. qot >= unsqmaxz) decent=.FALSE.
               ENDIF

               IF (decent) THEN
                 fz(iv,k,i,j)=fzeau(k,i,j)*yb
               ELSE
                 yr=cp(k)*yc+bp(k)*yb+ap(k)*ya
                 fz(iv,k,i,j)=fzeau(k,i,j)*yr
               END IF
             END DO   ! fin de boucle sur les traceurs

           ENDIF    ! fin de test sur la vitesse
         END DO   ! fin de boucle en k

       ENDIF     ! fin de test sur la hauteur d eau
     END DO     ! fin de la boucle en i
   END DO      ! fin de la boucle en j
!$OMP end do



  END SUBROUTINE sub_advz_quickupwind
