  SUBROUTINE parsub_advz_quickupwind(xe3,setvel,c,ll,niter,num_ws)
 
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE parsub_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_parsub_wquick
   !&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_parsub_wquick CPP key
   !&E                  is activated (see parsub_eq_wquick.h).
   !&E
   !&E ** Important remark about the settling velocity with the use of key_subs_part_eqsubdt cpp key
   !&E              To save computation time, we assume the substances have the same settling velocity
   !&E              in a group (since the group and the corresponding sub time step depend on
   !&E              the diameter and the density of the particule).
   !&E              This assumption is not necessary valid for the biology. If so, modify slightly
   !&E              the maximum settling velocity in variable.dat file, to create a new group.
   !&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       !  20004-10-04
   !&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
   !&E       !  2009-04  (A. Mano, V. Garnier) MPI2D
   !&E       !  2010-06  (V.M. Faure, N. Bousserez, M. Baklouti) coupling with eco3m
   !&E       !  2011-02 (B. Thouvenin) correction for key_subs_part_eqsubdt key + nb of iterations
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2016-04  (B.Thouvenin) only vertical transport due to settling velocoty
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
#if defined key_substance
   USE parameters
   USE comvars2d,   ONLY : ig,id,h0,hm,halfdt
   USE comvars3d,   ONLY : sig,dsigw,qmaxz
   USE comsubstance,ONLY : iv_ws_max,nv_ws_max,fz
   USE comfluxdo,   ONLY : fzeau_ws
#if defined key_siggen || defined key_gencoord
   USE toolgeom,  ONLY   : f_hzw,f_zu
#endif
#if defined key_sedim_MUSTANG || defined key_biolo
   USE comsubstance,ONLY : ws3
#endif

   IMPLICIT NONE

   !! * Arguments
   INTEGER,                                                              INTENT(in) :: ll,niter,num_ws
   REAL(KIND=rsh),                                                       INTENT(in) :: setvel
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),            INTENT(in) :: xe3
   REAL(KIND=rsh),DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: c

   !! * Local declarations
   LOGICAL                               :: decent
   INTEGER                               :: i,j,k,l1,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
   REAL(KIND=rsh),DIMENSION(1:kmax)      :: ws
#if defined key_siggen || defined key_gencoord
   REAL(KIND=rsh)                        :: hz3_w
#endif

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

   unsqmaxz=1.0_rsh/qmaxz


!$OMP DO SCHEDULE(RUNTIME),PRIVATE(i,j,k,ws,l1,iv,h3,        &
#if defined key_siggen || defined key_gencoord
!$OMP& hz3_w,                                                &
#endif
!$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)

     ! get the settling velocity
     ! -------------------------

   ! we assume all the substances of a group have the same settling velocity (iv=1 for each group)

#if defined key_sedim_MUSTANG || defined key_biolo
       iv=iv_ws_max(num_ws,1)
       ws(:)=-ws3(:,iv,i,j)
#else
       ! constant settling velocity
       ws(:)=-setvel
#endif


   !  sea surface height at t+
   ! --------------------------------------------
       h3=h0(i,j)+xe3(i,j)
       IF (h3 > hm) THEN

   ! estimate flux of particulate substance
   ! --------------------------------------
         DO k=1,kmax-1
#if defined key_siggen || defined key_gencoord
           hz3_w=f_hzw(h0(i,j),xe3(i,j),k,i,j)
           fzeau_ws(k,i,j)=(0.5_rsh*(ws(k)+ws(k+1)))*REAL(halfdt,rsh)/REAL(niter,rsh)
#else
           fzeau_ws(k,i,j)=(0.5_rsh*(ws(k)+ws(k+1)))*REAL(halfdt,rsh)/REAL(niter,rsh)
#endif
         END DO

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

         DO k=2,kmax-1
#if defined key_siggen || defined key_gencoord
           x1=f_zu(h0(i,j),xe3(i,j),k-1,i,j)
           x2=f_zu(h0(i,j),xe3(i,j),k  ,i,j)
           x3=f_zu(h0(i,j),xe3(i,j),k+1,i,j)
#else
           x1=sig(k-1)
           x2=sig(k)
           x3=sig(k+1)
#endif
           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
#if defined key_siggen || defined key_gencoord
           x1=f_zu(h0(i,j),xe3(i,j),k  ,i,j)
           x2=f_zu(h0(i,j),xe3(i,j),k+1,i,j)
           x3=f_zu(h0(i,j),xe3(i,j),k+2,i,j)
#else
           x1=sig(k)
           x2=sig(k+1)
           x3=sig(k+2)
#endif
           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
   ! ---------------------------------------------------

         DO l1=1,nv_ws_max(num_ws)
           iv=iv_ws_max(num_ws,l1)
           fz(iv,0:kmax,i,j)=0.0_rsh

           k=1
           IF (fzeau_ws(k,i,j) > 0.0_rsh) THEN
             fz(iv,k,i,j)=fzeau_ws(k,i,j)*c(iv,k,i,j)
           ELSE
             fz(iv,k,i,j)=fzeau_ws(k,i,j)*c(iv,k+1,i,j)
           END IF

           k=kmax-1
           IF (fzeau_ws(k,i,j) > 0.0_rsh) THEN
             fz(iv,k,i,j)=fzeau_ws(k,i,j)*c(iv,k,i,j)
           ELSE
             fz(iv,k,i,j)=fzeau_ws(k,i,j)*c(iv,k+1,i,j)
           END IF

           DO k=2,kmax-2
             IF (fzeau_ws(k,i,j) <= 0.0_rsh) THEN
               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_ws(k,i,j)*yb
               ELSE
                 yr=cm(k)*yc+bm(k)*yb+am(k)*ya
                 fz(iv,k,i,j)=fzeau_ws(k,i,j)*yr
               END IF

             ELSE ! test sur la vitesse cas w >0

               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_ws(k,i,j)*yb
               ELSE
                 yr=cp(k)*yc+bp(k)*yb+ap(k)*ya
                 fz(iv,k,i,j)=fzeau_ws(k,i,j)*yr
               END IF

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

         END DO   ! fin de boucle sur les traceurs

       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

#else

   USE parameters, ONLY : rsh,kmax,limin,limax,ljmin,ljmax,liminm1,liminm2,limaxp2,ljminm1,ljminm2,ljmaxp2
   INTEGER, INTENT(in) :: ll,niter,num_ws
   REAL(KIND=rsh), INTENT(in) :: setvel
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in) :: xe3
   REAL(KIND=rsh),DIMENSION(kmax,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: c

#endif

  END SUBROUTINE parsub_advz_quickupwind
