  SUBROUTINE tssub_advdifx_quickupwind(sal,temp &
#ifdef key_substance
                                       ,c       &
#endif
                                      )

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE tssub_advdifx_quickupwind  ***
   !&E
   !&E ** Purpose : Estimate of zonal fluxes of temperature and salinity
   !&E              and other substances (advection and diffusion)
   !&E              Advection uses TVD scheme (QUICK and upwind) owing to limiter.
   !&E
   !&E ** Description : No horizontal diffusion at lower depths ( < hminkxky)
   !&E
   !&E ** Cas cyclique en y : correct for eastward propagation
   !&E                        uncorrect for westward propagation (stencil not long enough)
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays : htx,hex,hx,fxeau,sal,temp,c(cv_wat),hzex,dx_u 
   !&E                     fsalx,ftempx,fx 
   !&E  
   !&E ** Modified variables : fsalx,ftempx,fx
   !&E
   !&E ** Reference : Kowalik Z. and T. S. Murty : Numerical modeling of ocean
   !&E                dynamics in Advanced series on ocean engineering, vol 5.
   !&E
   !&E ** History :
   !&E           !  2007-09 (V. Garnier) homogeneneity of dx dy...
   !&E           !  2008-09 (V. Garnier) correction for western boundary
   !&E           !  2009-04 (A. Mano) MPI2D
   !&E           !  2009-10  (F. Cayocca) correction for sedimentology
   !&E           !  2010-10 (T. Duhaut) introduction of diffusion along isopycnes
   !&E           !  2016-04 (B.Thouvenin) applied to all substances (dissolved and particulate)
   !&E           !  2016-12  (B.Thouvenin) version sedim_MUSTANG : only mixsed (mudsand not operationnal)
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d,    ONLY : ig,id,htx,hex,hx,kx,qmax,h0fond,hminkxky,dx_u,halfdt,l_obc_cycl_x,dy_u
#if defined key_siggen || defined key_gencoord
   USE comsiggen,    ONLY : hzex
#endif
   USE comfluxdo,    ONLY : fsalx,fxeau
   USE cometemp,     ONLY : ftempx
#ifdef key_substance
   USE comsubstance, ONLY : nv_adv,nv_grav,fx
#ifdef key_sedim_MUSTANG
   USE comMUSTANG,    ONLY : corflux
#endif
#endif
   USE_MPI toolmpi,      ONLY : ex_i_rsh,ex_i_rlg

   IMPLICIT NONE

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

#ifdef key_substance
   REAL(KIND=rsh), DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT(in) :: c
#endif

   !! * Local declarations
   LOGICAL                 :: decentt,decents
   INTEGER                 :: i,j,k
   REAL(KIND=rsh)          :: unsqmax,uu,yr,qot,dyds,dygs,dydt
   REAL(KIND=rsh)          :: dygt,kxb,yas,yat,ybs,ybt,ycs,yct,epsi
#ifdef key_substance
   LOGICAL, DIMENSION(nv_adv)           :: decent
   INTEGER                              :: iv
   REAL(KIND=rsh)                       :: dyg
   REAL(KIND=rsh), DIMENSION(nv_adv)    :: ya,yb,yc,dyd
#endif

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

   unsqmax=1.0_rsh/qmax
   epsi=1.0e-10

!$OMP DO schedule(runtime),private(i,j,k,kxb               &
!$OMP& ,uu,yr,qot,decents,decentt,dyds,dygs,dydt,dygt      &
#ifdef key_substance
!$OMP& ,iv,ya,yb,yc,decent,dyd,dyg                         &
#endif
!$OMP& ,yas,yat,ybs,ybt,ycs,yct)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)),MIN0(limax,id(j)-1)
       DO k=1,kmax
         IF (htx(i,j) <= 0.0_rsh .OR. hx(i,j) == -valmanq-h0fond) THEN
#if ! defined key_offline
           fsalx(k,i,j)=0.0_rlg
           ftempx(k,i,j)=0.0_rlg
#endif
#ifdef key_substance
           !refonte subtances advection avril 2016
           fx(nv_grav+1:nv_adv,k,i,j)=0.0_rsh
#endif
         ELSE
           kxb=kx
           IF (htx(i,j) < hminkxky) kxb=0.0_rsh ! no diffusion for small water depth
           uu=fxeau(k,i,j)
           IF (uu <= 0.0_rsh) THEN  ! westward advection
#if ! defined key_offline
             ybs=sal(k,i+1,j)
             ycs=sal(k,i,j)
             ybt=temp(k,i+1,j)
             yct=temp(k,i,j)
#endif
#ifdef key_substance
             yb(nv_grav+1:nv_adv)=c(nv_grav+1:nv_adv,k,i+1,j)
             yc(nv_grav+1:nv_adv)=c(nv_grav+1:nv_adv,k,i,j)
#endif
             IF (htx(i+1,j) > 0.0_rsh .AND. (i+1 < id(j) .OR. l_obc_cycl_x)) THEN
#if ! defined key_offline
               yas=sal(k,i+2,j)
               yat=temp(k,i+2,j)
#endif
#ifdef key_substance
               ya(nv_grav+1:nv_adv)=c(nv_grav+1:nv_adv,k,i+2,j)
#endif
             ELSE
#if ! defined key_offline
               yas=ybs
               yat=ybt
#endif
#ifdef key_substance
               ya(nv_grav+1:nv_adv)=yb(nv_grav+1:nv_adv)
#endif
             END IF

#if ! defined key_offline
             dyds=ybs-ycs
             dygs=ybs-yas
             dydt=ybt-yct
             dygt=ybt-yat

             decents=(dyds*dygs) >= 0.0_rsh
             IF (.NOT.decents) THEN
               qot=ABS(dyds/MAX(ABS(dygs),epsi))
               IF (qot <= qmax .AND. qot >= unsqmax) THEN
                 decents=.FALSE.
               ELSE
                 decents=.TRUE.
               END IF
             END IF

             decentt=(dydt*dygt) >= 0.0_rsh
             IF (.NOT.decentt) THEN
               qot=ABS(dydt/MAX(ABS(dygt),epsi))
               IF (qot <= qmax .AND. qot >= unsqmax) THEN
                 decentt=.FALSE.
               ELSE
                 decentt=.TRUE.
               END IF
             END IF
#endif

#ifdef key_substance
             DO iv=nv_grav+1,nv_adv
               dyd(iv)=yb(iv)-yc(iv)
               dyg=yb(iv)-ya(iv)
               decent(iv)=(dyd(iv)*dyg) >= 0.0_rsh
               IF (.NOT.decent(iv)) THEN
                 qot=ABS(dyd(iv)/MAX(ABS(dyg),epsi))
                 IF (qot <= qmax .AND. qot >= unsqmax) THEN
                   decent(iv)=.FALSE.
                 ELSE
                   decent(iv)=.TRUE.
                 END IF
               END IF
             END DO
#endif

           ELSE  ! eastward advection

#if ! defined key_offline
             ybs=sal(k,i,j)
             ycs=sal(k,i+1,j)
             ybt=temp(k,i,j)
             yct=temp(k,i+1,j)
#endif
#ifdef key_substance
             yb(nv_grav+1:nv_adv)=c(nv_grav+1:nv_adv,k,i,j)
             yc(nv_grav+1:nv_adv)=c(nv_grav+1:nv_adv,k,i+1,j)
#endif
             IF (htx(i-1,j) > 0.0_rsh .AND. (i > ig(j) .OR. l_obc_cycl_x)) THEN
!              IF (hex(i-1,j) > 0.0_rsh) THEN
#if ! defined key_offline
               yas=sal(k,i-1,j)
               yat=temp(k,i-1,j)
#endif
#ifdef key_substance
               ya(nv_grav+1:nv_adv)=c(nv_grav+1:nv_adv,k,i-1,j)
#endif
             ELSE
#if ! defined key_offline
               yas=ybs
               yat=ybt
#endif
#ifdef key_substance
               ya(nv_grav+1:nv_adv)=yb(nv_grav+1:nv_adv)
#endif
             END IF

#if ! defined key_offline
             dyds=ycs-ybs
             dygs=ybs-yas
             dydt=yct-ybt
             dygt=ybt-yat

             decents=(dyds*dygs) <= 0.0_rsh
             IF (.NOT.decents) THEN
               qot=ABS(dyds/MAX(ABS(dygs),epsi))
               IF (qot <= qmax .AND. qot >= unsqmax) THEN
                 decents=.FALSE.
               ELSE
                 decents=.TRUE.
               END IF
             END IF

             decentt=(dydt*dygt) <= 0.0_rsh
             IF (.NOT.decentt) THEN
               qot=ABS(dydt/MAX(ABS(dygt),epsi))
               IF (qot <= qmax .AND. qot >= unsqmax) THEN
                 decentt=.FALSE.
               ELSE
                 decentt=.TRUE.
               END IF
             END IF
#endif

#ifdef key_substance
             DO iv=nv_grav+1,nv_adv
               dyd(iv)=yc(iv)-yb(iv)
               dyg=yb(iv)-ya(iv)
               decent(iv)=(dyd(iv)*dyg) <= 0.0_rsh
               IF (.NOT.decent(iv)) THEN
                 qot=ABS(dyd(iv)/MAX(ABS(dyg),epsi))
                 IF (qot <= qmax .AND. qot >= unsqmax) THEN
                   decent(iv)=.FALSE.
                 ELSE
                   decent(iv)=.TRUE.
                 END IF
               END IF
             END DO
#endif
           END IF    ! ends test on velocity sign


#if ! defined key_offline

#if defined key_ts_diffh_geo || defined key_ts_diffh_rho
   ! advection flux of salinity alone
   ! --------------------------------

           IF (decents) THEN
             fsalx(k,i,j)=uu*ybs*dy_u(i,j)
           ELSE
             yr=0.375_rsh*ycs+0.75_rsh*ybs-0.125_rsh*yas
             fsalx(k,i,j)=uu*yr*dy_u(i,j)
           END IF

   ! advection flux of temperature alone
   ! -----------------------------------

           IF (decentt) THEN
             ftempx(k,i,j)=uu*ybt*dy_u(i,j)
           ELSE
             yr=0.375_rsh*yct+0.75_rsh*ybt-0.125_rsh*yat
             ftempx(k,i,j)=uu*yr*dy_u(i,j)
           END IF

#else

   ! advection/diffusion flux of salinity
   ! ------------------------------------

           IF (decents) THEN
#if defined key_siggen || defined key_gencoord
             fsalx(k,i,j)=uu*ybs-hzex(k,i,j)*kxb*dyds/dx_u(i,j)*halfdt
#else
             fsalx(k,i,j)=uu*ybs-hex(i,j)*kxb*dyds/dx_u(i,j)*halfdt
#endif
           ELSE
             yr=0.375_rsh*ycs+0.75_rsh*ybs-0.125_rsh*yas
#if defined key_siggen || defined key_gencoord
             fsalx(k,i,j)=uu*yr-hzex(k,i,j)*kxb*dyds/dx_u(i,j)*halfdt
#else
             fsalx(k,i,j)=uu*yr-hex(i,j)*kxb*dyds/dx_u(i,j)*halfdt
#endif
           END IF

   ! advection/diffusion flux of temperature
   ! ---------------------------------------

           IF (decentt) THEN
#if defined key_siggen || defined key_gencoord
             ftempx(k,i,j)=uu*ybt-hzex(k,i,j)*kxb*dydt/dx_u(i,j)*halfdt
#else
             ftempx(k,i,j)=uu*ybt-hex(i,j)*kxb*dydt/dx_u(i,j)*halfdt
#endif
           ELSE
             yr=0.375_rsh*yct+0.75_rsh*ybt-0.125_rsh*yat
#if defined key_siggen || defined key_gencoord
             ftempx(k,i,j)=uu*yr -hzex(k,i,j)*kxb*dydt/dx_u(i,j)*halfdt
#else
             ftempx(k,i,j)=uu*yr -hex(i,j)*kxb*dydt/dx_u(i,j)*halfdt
#endif
           END IF

#endif

   ! key_offline
#endif

#ifdef key_substance

   ! advection flux of other substances
   ! ----------------------------------

           DO iv=nv_grav+1,nv_adv

#if defined key_ts_diffh_geo || defined key_ts_diffh_rho
             IF (decent(iv)) THEN
               fx(iv,k,i,j)=fxeau(k,i,j)*yb(iv)*dy_u(i,j)
             ELSE
               yr=0.375_rsh*yc(iv)+0.75_rsh*yb(iv)-0.125_rsh*ya(iv)
               fx(iv,k,i,j)=fxeau(k,i,j)*yr*dy_u(i,j)
             END IF

#else
             IF (decent(iv)) THEN
#if defined key_siggen || defined key_gencoord
               fx(iv,k,i,j)=fxeau(k,i,j)*yb(iv)-hzex(k,i,j)*kxb*dyd(iv)/dx_u(i,j)*halfdt
#else
               fx(iv,k,i,j)=fxeau(k,i,j)*yb(iv)-hex(i,j)*kxb*dyd(iv)/dx_u(i,j)*halfdt
#endif
             ELSE
               yr=0.375_rsh*yc(iv)+0.75_rsh*yb(iv)-0.125_rsh*ya(iv)
#if defined key_siggen || defined key_gencoord
               fx(iv,k,i,j)=fxeau(k,i,j)*yr - hzex(k,i,j)*kxb*dyd(iv)/dx_u(i,j)*halfdt
#else
               fx(iv,k,i,j)=fxeau(k,i,j)*yr - hex(i,j)*kxb*dyd(iv)/dx_u(i,j)*halfdt
#endif
             END IF
#endif
#ifdef key_sedim_MUSTANG
               IF (k.eq.1) fx(iv,k,i,j)=fx(iv,k,i,j)*corflux(iv,i,j)
#endif

           END DO
#endif

          END IF ! ends test on water depth
        END DO ! ends loop on k
      END DO ! ends loop on i
   END DO ! ends loop on j
!$OMP END DO

OMPMPI barrier
OMPMPI master
# if ! defined key_offline
   CALL_MPI ex_i_rlg(-1,0,kmax,liminm1,limax,ljmin,ljmax,fsalx(:,liminm1:limax,ljmin:ljmax))
   CALL_MPI ex_i_rlg(-1,0,kmax,liminm1,limax,ljmin,ljmax,ftempx(:,liminm1:limax,ljmin:ljmax))
#endif
#ifdef key_substance
   CALL_MPI ex_i_rsh(-1,0,nv_adv*kmax,liminm1,limax,ljmin,ljmax,fx(:,:,liminm1:limax,ljmin:ljmax))
#endif
OMPMPI end master
OMPMPI barrier
# if ! defined key_offline
OMPMPI flush(fsalx,ftempx)
#endif
#ifdef key_substance
OMPMPI flush(fx)
#endif

  END SUBROUTINE tssub_advdifx_quickupwind
