  SUBROUTINE ts_advz_quickupwind(xe,sal,temp)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE ts_advz_quickupwind  ***
   !&E
   !&E ** Purpose : compute vertical fluxes of 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 in "adv" only if the key
   !&E              key_tssub_wquickupwind is activated (see tssub_eq_wquick.h).
   !&E
   !&E ** Called by : adv
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Used ij-arrays : xe,h0,fzeau,sal,temp
   !&E
   !&E ** Modified variables : ftempz,fsalz
   !&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
   !&E---------------------------------------------------------------------
   !! * Modules used

   USE parameters
   USE comvars2d,  ONLY : id,ig,h0,hm
   USE comvars3d,  ONLY : sig,dsigw,qmaxz
   USE comfluxdo,  ONLY : fzeau
   USE comfluvert, ONLY : fsalz,ftempz
#if defined key_siggen || defined key_gencoord
   USE toolgeom,  ONLY   : f_zu
#endif
   IMPLICIT NONE

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

   !! * Local declarations
   LOGICAL                                  :: decentt,decents
   INTEGER                                  :: i,j,k
   REAL(kind=rsh)                           :: unsqmaxz
   REAL(kind=rsh)                           :: qot,dyds,dygs,dydt,dygt,h3,yr
   REAL(kind=rsh)                           :: yas,yat,ybs,ybt,ycs,yct,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) &
!$OMP& private(i,j,k,h3,ybs,ycs,ybt,yct,yas,yat,dyds,dygs,dydt &
!$OMP&, dygt,decents,qot,decentt,yr,x1,x2,x3,xx,ap,bp,cp,am,bm,cm)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)),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
#if defined key_siggen || defined key_gencoord
           x1=f_zu(h0(i,j),xe(i,j),k-1,i,j)
           x2=f_zu(h0(i,j),xe(i,j),k  ,i,j)
           x3=f_zu(h0(i,j),xe(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),xe(i,j),k  ,i,j)
           x2=f_zu(h0(i,j),xe(i,j),k+1,i,j)
           x3=f_zu(h0(i,j),xe(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 (fond et surface)
!---------------------------------------------------------------------

         ftempz(0,i,j)=0.0_rsh
         ftempz(kmax,i,j)=0.0_rsh
         fsalz(0,i,j)=0.0_rsh
         fsalz(kmax,i,j)=0.0_rsh
         k=1
         IF (fzeau(k,i,j) > 0.0_rlg) THEN
           fsalz(k,i,j)=fzeau(k,i,j)*sal(k,i,j)
           ftempz(k,i,j)=fzeau(k,i,j)*temp(k,i,j)
         ELSE
           fsalz(k,i,j)=fzeau(k,i,j)*sal(k+1,i,j)
           ftempz(k,i,j)=fzeau(k,i,j)*temp(k+1,i,j)
         END IF

         k=kmax-1
         IF (fzeau(k,i,j) > 0.0_rlg) THEN
           fsalz(k,i,j)=fzeau(k,i,j)*sal(k,i,j)
           ftempz(k,i,j)=fzeau(k,i,j)*temp(k,i,j)
         ELSE
           fsalz(k,i,j)=fzeau(k,i,j)*sal(k+1,i,j)
           ftempz(k,i,j)=fzeau(k,i,j)*temp(k+1,i,j)
         END IF

! au milieu de la colonne d eau

         DO k=2,kmax-2
           IF(fzeau(k,i,j) <= 0.0_rlg) THEN
             ybs=sal(k+1,i,j)
             ycs=sal(k,i,j)
             ybt=temp(k+1,i,j)
             yct=temp(k,i,j)
             IF(k <= kmax-2) THEN
               yas=sal(k+2,i,j)
               yat=temp(k+2,i,j)
             ELSE
               yas=ybs
               yat=ybt
             END IF

             dyds=(ybs-ycs)/dsigw(k)
             dygs=(ybs-yas)/dsigw(k+1)
             dydt=(ybt-yct)/dsigw(k)
             dygt=(ybt-yat)/dsigw(k+1)

             decents=(dyds*dygs) >= 0.0_rsh
             IF (.NOT.decents) THEN
               qot=ABS(dyds/dygs)
               decents=.TRUE.
               IF (qot <= qmaxz .AND. qot >= unsqmaxz) decents=.FALSE.
             END IF

             decentt=(dydt*dygt) >= 0.0_rsh
             IF (.NOT.decentt) THEN
               qot=ABS(dydt/dygt)
               decentt=.TRUE.
               IF (qot <= qmaxz .AND. qot >= unsqmaxz) decentt=.FALSE.
             END IF

!  salinite
!-------------

             IF (decents) THEN
               fsalz(k,i,j)=fzeau(k,i,j)*ybs
             ELSE
               yr=cm(k)*ycs+bm(k)*ybs+am(k)*yas
               fsalz(k,i,j)=fzeau(k,i,j)*yr
             END IF

! temperature
!-------------

             IF (decentt) THEN
               ftempz(k,i,j)=fzeau(k,i,j)*ybt
             ELSE
               yr=cm(k)*yct+bm(k)*ybt+am(k)*yat
               ftempz(k,i,j)=fzeau(k,i,j)*yr
             END IF   

           ELSE ! deuxieme option du test sur la vitesse verticale (cas w >0)

             ybs=sal(k,i,j)
             ycs=sal(k+1,i,j)
             ybt=temp(k,i,j)
             yct=temp(k+1,i,j)
             IF (k >= 2) THEN
               yas=sal(k-1,i,j)
               yat=temp(k-1,i,j)
             ELSE
               yas=ybs
               yat=ybt
             END IF

             dyds=(ycs-ybs)/dsigw(k)
             dygs=(ybs-yas)/dsigw(k-1)
             dydt=(yct-ybt)/dsigw(k)
             dygt=(ybt-yat)/dsigw(k-1)

             decents=(dyds*dygs) <= 0.0_rsh
             IF (.NOT.decents) THEN
               qot=ABS(dyds/dygs)
               decents=.TRUE.
               IF (qot <= qmaxz .AND. qot >= unsqmaxz) decents=.FALSE.
             END IF

             decentt=(dydt*dygt) <= 0.0_rsh
             IF (.NOT.decentt) THEN
               qot=ABS(dydt/dygt)
               decentt=.TRUE.
               IF (qot <= qmaxz .AND. qot >= unsqmaxz) decentt=.FALSE.
             END IF

!  salinite
!-------------
!
             IF (decents) THEN
               fsalz(k,i,j)=fzeau(k,i,j)*ybs
             ELSE
               yr=cp(k)*ycs+bp(k)*ybs+ap(k)*yas
               fsalz(k,i,j)=fzeau(k,i,j)*yr
             END IF

! temperature
!-------------

             IF (decentt) THEN
               ftempz(k,i,j)=fzeau(k,i,j)*ybt
             ELSE
               yr=cp(k)*yct+bp(k)*ybt+ap(k)*yat
               ftempz(k,i,j)=fzeau(k,i,j)*yr
             END IF
        
           END IF     ! fin du test sur les vitesses

         END DO    ! fin de la boucle en k

       ENDIF     ! fin du 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 ts_advz_quickupwind
