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

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE tssub_advdifx_ultimatequickestmacho  ***
   !&E
   !&E ** Purpose : Estimate of zonal fluxes of temperature and salinity
   !&E              (advection and diffusion)
   !&E              Advection uses QUICKEST scheme with ULTIMATE limiter and
   !&E              its MACHO multidimensional extension.
   !&E
   !&E ** Description : No horizontal diffusion at lower depths ( < hminkxky)
   !&E                  - case of cyclicity in 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,hty,hey,hy,fxeau,fyeau,sal,temp,c(cv_wat),hzex,dx_u,dy,
   !&E                     sal2,temp2,c2,fsalx,ftempx,fx 
   !&E
   !&E ** Modified variables : fsalx,ftempx,fx
   !&E
   !&E ** Reference : Leonard B. P., Lock A. P. and M. K. MacVean (1996) :
   !&E                Conservative Explicit Unrestricted-Time-Step
   !&E                Multidimensional COnstancy-Preserving Advections Schemes.
   !&E                Monthly Weather Review, 124, 2588:2605
   !&E
   !&E ** History :
   !&E           !  2008-09 (L. Debreu, T. Duhaut, M. Honnorat) original code
   !&E           !  2008-09 (V. Garnier) generalized sigma, correction for western boundary
   !&E           !  2008-10 (V. Garnier) optimization
   !&E           !  2008-12 (F. Dumas) correction test on ig and id
   !&E           !  2009-04 (A. Mano) MPI2D
   !&E           !  2010-05 (V. Garnier) OMP
   !&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,h0fond,hminkxky,dx_u,dy,hty,hey,    &
                            jb,jh,hy,halfdt,l_obc_cycl_x,dy_u
#if defined key_siggen || defined key_gencoord
   USE comsiggen,    ONLY : hzex,hzey
#endif
   USE comfluxdo,    ONLY : fsalx,fxeau,fyeau
   USE cometemp,     ONLY : ftempx
   USE comvarp2d,    ONLY : sal2,temp2
#ifdef key_substance
   USE comsubstance, ONLY : nv_adv,nv_grav,nv_tot,fx,c2
#ifdef key_sedim_MUSTANG
   USE comMUSTANG,    ONLY : corflux
#endif
#endif
   USE_MPI toolmpi,      ONLY : ex_i_rsh,ex_i_rlg,ex_peri_x_rsh

   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
   INTEGER                   :: i,j,k
   REAL(KIND=rsh), PARAMETER :: dsmall = 1.0d-20
   REAL(KIND=rsh)            :: uu,dyds,dygs,dydt,d0,d1,uuu,courant,vv,  &
                                dygt,kxb,yas,yat,ybs,ybt,ycs,yct,thetap, &
                                psip,deltat,deltas
   REAL(KIND=rlg), DIMENSION(kmax) :: fvb,fvt
   REAL(KIND=rsh), DIMENSION(kmax) :: salb,salt,tempb,tempt
#ifdef key_substance
   INTEGER                              :: iv
   REAL(KIND=rsh), DIMENSION(nv_adv)    :: dyd,dyg,ya,yb,yc,delta
   REAL(KIND=rsh), DIMENSION(1:nv_adv,kmax) :: cb,ct
#endif

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

   ! estimate of intermediate tracers (from upwind scheme)
   ! -----------------------------------------------------

   ! at open boundaries (east and west)

!$OMP SINGLE

   IF_MPI (limin==imin) THEN
#if ! defined key_offline
     temp2(:,imin:imin+1,ljmin:ljmax) = temp(:,imin:imin+1,ljmin:ljmax)
     sal2(:,imin:imin+1,ljmin:ljmax) = sal(:,imin:imin+1,ljmin:ljmax)
#endif
#ifdef key_substance
     !refonte subtances advestion avril 2016
     c2(1:nv_adv,:,imin:imin+1,ljmin:ljmax) = c(1:nv_adv,:,imin:imin+1,ljmin:ljmax)
#endif
   ENDIF_MPI
   IF_MPI (limax==imax) THEN
#if ! defined key_offline
     temp2(:,imax,ljmin:ljmax) = temp(:,imax,ljmin:ljmax)
     sal2(:,imax,ljmin:ljmax) = sal(:,imax,ljmin:ljmax)
#endif
#ifdef key_substance
     c2(nv_grav+1:nv_adv,:,imax,ljmin:ljmax) = c(nv_grav+1:nv_adv,:,imax,ljmin:ljmax)
#endif
   ENDIF_MPI

!$OMP END SINGLE

   ! inside the domain

!$OMP DO SCHEDULE(RUNTIME),PRIVATE(i,j,k,     &
!$OMP& salt,tempt,fvt,salb,tempb,fvb,         &
#ifdef key_substance
!$OMP& iv,ct,cb,                              &
#endif
!$OMP& vv)

   DO i=liminp2,limaxm1
     DO j=ljminp1m1,ljmax
       DO k=1,kmax

#ifdef key_MPI_2D
         IF (j == jb(i) .OR. j==ljminp1m1) THEN
#else
         IF (j == jb(i) .OR. j==ljminp1-1) THEN
#endif

           IF (hty(i,j) <= 0.0_rsh .OR. hy(i,j)==-valmanq-h0fond) THEN
#if ! defined key_offline
             salt(k)=0.0_rsh       ! no flux advection from this cell
             tempt(k)=0.0_rsh
#endif
#ifdef key_substance
             ct(nv_grav+1:nv_adv,k)=0.0_rsh
#endif
             fvt(k)=0.0_rlg
           ELSE

             fvt(k)=fyeau(k,i,j)   ! h*vz*dt/2

             ! upwind scheme 
             IF (fvt(k) <= 0.0_rlg) THEN
#if ! defined key_offline
               salt(k)=sal(k,i,j+1)
               tempt(k)=temp(k,i,j+1)
#endif
#ifdef key_substance
               ct(nv_grav+1:nv_adv,k)=c(nv_grav+1:nv_adv,k,i,j+1)
#endif
             ELSE
#if ! defined key_offline
               salt(k)=sal(k,i,j)
               tempt(k)=temp(k,i,j)
#endif
#ifdef key_substance
               ct(nv_grav+1:nv_adv,k)=c(nv_grav+1:nv_adv,k,i,j)
#endif
             END IF   ! ends test on the direction of propagation

           END IF   ! ends test on the water depth

#ifdef key_MPI_2D
         ELSE IF (j >= jb(i)+1 .AND. j <= jh(i)-1 .AND. j>ljminp1m1) THEN
#else
         ELSE IF (j >= jb(i)+1 .AND. j <= jh(i)-1 .AND. j>ljminp1-1) THEN
#endif

           ! bottom of the stencil
#if ! defined key_offline
           salb(k)=salt(k)
           tempb(k)=tempt(k)
#endif
#ifdef key_substance
           cb(nv_grav+1:nv_adv,k)=ct(nv_grav+1:nv_adv,k)
#endif
           fvb(k)=fvt(k)

           ! right of the stencil
           IF (hty(i,j) <= 0.0_rsh .OR. hy(i,j)==-valmanq-h0fond) THEN
#if ! defined key_offline
             salt(k)=0.0_rsh       ! no flux advection from this cell
             tempt(k)=0.0_rsh
#endif
#ifdef key_substance
             ct(nv_grav+1:nv_adv,k)=0.0_rsh
#endif
             fvt(k)=0.0_rlg
           ELSE

             fvt(k)=fyeau(k,i,j)   ! h*vz*dt/2

             ! upwind scheme 
             IF (fvt(k) <= 0.0_rlg) THEN
#if ! defined key_offline
               salt(k)=sal(k,i,j+1)
               tempt(k)=temp(k,i,j+1)
#endif
#ifdef key_substance
               ct(nv_grav+1:nv_adv,k)=c(nv_grav+1:nv_adv,k,i,j+1)
#endif
             ELSE
#if ! defined key_offline
               salt(k)=sal(k,i,j)
               tempt(k)=temp(k,i,j)
#endif
#ifdef key_substance
               ct(nv_grav+1:nv_adv,k)=c(nv_grav+1:nv_adv,k,i,j)
#endif
             END IF   ! ends test on the direction of propagation

           END IF   ! ends test on the water depth

   ! estimate of intermediate tracers (from upwind scheme)
           IF ( (fvb(k)>0.0_rlg) .AND. (fvt(k)>0.0_rlg) ) THEN
!cval toseelater Normal qu il n y ait pas de sous pas de temps ?
#if defined key_siggen || defined key_gencoord
             vv = fvb(k)/(halfdt*hzey(k,i,j-1))
#else
             vv = fvb(k)/(halfdt*hey(i,j-1))
#endif
           ELSE IF ( (fvb(k)<0.0_rlg) .AND. (fvt(k)<0.0_rlg) ) THEN
#if defined key_siggen || defined key_gencoord
             vv = fvt(k)/(halfdt*hzey(k,i,j))
#else
             vv = fvt(k)/(halfdt*hey(i,j))
#endif
           ELSE
             ! no time evolution of the intermediate temperature
             ! along the boundaries
             vv = 0.0_rsh
           END IF

           ! apply euler (time) - upwind (x) scheme
#if ! defined key_offline
           temp2(k,i,j) = temp(k,i,j) - (halfdt/dy(i,j))*vv*(tempt(k)-tempb(k))
           sal2(k,i,j)  = sal(k,i,j)  - (halfdt/dy(i,j))*vv*(salt(k)-salb(k))
#endif
#ifdef key_substance
           c2(nv_grav+1:nv_adv,k,i,j)  = c(nv_grav+1:nv_adv,k,i,j)  &
             - (halfdt/dy(i,j))*vv*(ct(nv_grav+1:nv_adv,k)-cb(nv_grav+1:nv_adv,k))
#endif

         END IF   ! ends test for selecting right j loop

       END DO   ! ends loop on k
     END DO   ! ends loop on i
   END DO   ! ends loop on j
!$OMP END DO

!$OMP SINGLE
   ! periodic conditions along x-axis
   IF (l_obc_cycl_x) THEN
#ifdef key_MPI_2D
#ifdef key_shiftobc
     CALL_MPI ex_peri_x_rsh(3,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
     CALL_MPI ex_peri_x_rsh(3,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#ifdef key_substance
     CALL_MPI ex_peri_x_rsh(3,nv_adv*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,c2(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#else
     CALL_MPI ex_peri_x_rsh(1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
     CALL_MPI ex_peri_x_rsh(1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#ifdef key_substance
     CALL_MPI ex_peri_x_rsh(1,nv_adv*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,c2(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#endif
#else
     DO j=ljmin,ljmax
   ! western boundary
       sal2(:,ig(j),j) = sal2(:,id(j)-1,j)
       sal2(:,ig(j)-1,j) = sal2(:,id(j)-2,j)
       temp2(:,ig(j),j) = temp2(:,id(j)-1,j)
       temp2(:,ig(j)-1,j) = temp2(:,id(j)-2,j)
#ifdef key_substance
       c2(nv_grav+1:nv_adv,:,ig(j),j)=c2(nv_grav+1:nv_adv,:,id(j)-1,j)
       c2(nv_grav+1:nv_adv,:,ig(j)-1,j)=c2(nv_grav+1:nv_adv,:,id(j)-2,j)
#endif
#ifdef key_shiftobc
       sal2(:,ig(j)-2,j) = sal2(:,id(j)-3,j)
       temp2(:,ig(j)-2,j) = temp2(:,id(j)-3,j)
#ifdef key_substance
       c2(nv_grav+1:nv_adv,:,ig(j)-2,j)=c2(nv_grav+1:nv_adv,:,id(j)-3,j)
#endif
#endif
   ! eastern boundary
       sal2(:,id(j),j) = sal2(:,ig(j)+1,j)
       temp2(:,id(j),j) = temp2(:,ig(j)+1,j)
#ifdef key_substance
       c2(nv_grav+1:nv_adv,:,id(j),j)=c2(nv_grav+1:nv_adv,:,ig(j)+1,j)
#endif
#ifdef key_shiftobc
       sal2(:,id(j)+1,j) = sal2(:,ig(j)+2,j)
       sal2(:,id(j)+2,j) = sal2(:,ig(j)+3,j)
       temp2(:,id(j)+1,j) = temp2(:,ig(j)+2,j)
       temp2(:,id(j)+2,j) = temp2(:,ig(j)+3,j)
#ifdef key_substance
       c2(nv_grav+1:nv_adv,:,id(j)+1,j)=c2(nv_grav+1:nv_adv,:,ig(j)+2,j)
       c2(nv_grav+1:nv_adv,:,id(j)+2,j)=c2(nv_grav+1:nv_adv,:,ig(j)+3,j)
#endif
#endif
     END DO
#endif
   END IF  ! ends periodic conditions along x-axis
!$OMP END SINGLE

OMPMPI barrier
OMPMPI master
#if ! defined key_offline
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#ifdef key_substance
! to avoid temporary copy   CALL_MPI ex_i_rsh(-1,2,(nv_adv)*kmax,liminm1,limaxp2,ljmin,ljmax,c2(1:nv_adv,:,liminm1:limaxp2,ljmin:ljmax))
   CALL_MPI ex_i_rsh(-1,2,nv_tot*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,c2(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
OMPMPI end master
OMPMPI barrier
#if ! defined key_offline
OMPMPI flush(sal2,temp2)
#endif
#ifdef key_substance
OMPMPI flush(c2)
#endif

   ! estimate of fluxes (from ultimate quickest scheme)
   ! --------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME),PRIVATE(i,j,k,kxb               &
!$OMP& ,dyds,dygs,dydt,dygt,yas,yat,ybs,ybt,ycs,yct        &
#ifdef key_substance
!$OMP& ,iv,dyd,dyg,ya,yb,yc,delta                          &
#endif
!$OMP& ,deltat,deltas,d0,d1,uu,uuu,courant,thetap,psip)

   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 defined key_siggen || defined key_gencoord
           uuu = uu/(halfdt*hzex(k,i,j))
#else
           uuu = uu/(halfdt*hex(i,j))
#endif
           courant = ABS(uuu)*halfdt/dx_u(i,j)
           d0=(1.0_rsh/6.0_rsh)*(2.0_rsh-courant)*(1.0_rsh-courant)
           d1=(1.0_rsh/6.0_rsh)*(1.0_rsh-courant)*(1.0_rsh+courant)
             
           IF (uu <= 0.0_rsh) THEN  ! westward advection
#if ! defined key_offline
             ybs=sal2(k,i+1,j)
             ycs=sal2(k,i,j)
             ybt=temp2(k,i+1,j)
             yct=temp2(k,i,j)
#endif
#ifdef key_substance
             yb(nv_grav+1:nv_adv)=c2(nv_grav+1:nv_adv,k,i+1,j)
             yc(nv_grav+1:nv_adv)=c2(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=sal2(k,i+2,j)
               yat=temp2(k,i+2,j)
#endif
#ifdef key_substance
               ya(nv_grav+1:nv_adv)=c2(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=ycs-ybs
             dygs=ybs-yas
             dydt=yct-ybt
             dygt=ybt-yat
             deltas=ybs-ycs
             deltat=ybt-yct
#endif
#ifdef key_substance
             dyd(nv_grav+1:nv_adv)=yc(nv_grav+1:nv_adv)-yb(nv_grav+1:nv_adv)
             dyg(nv_grav+1:nv_adv)=yb(nv_grav+1:nv_adv)-ya(nv_grav+1:nv_adv)
             delta(nv_grav+1:nv_adv)=yb(nv_grav+1:nv_adv)-yc(nv_grav+1:nv_adv)
#endif

           ELSE  ! eastward advection

#if ! defined key_offline
             ybs=sal2(k,i,j)
             ycs=sal2(k,i+1,j)
             ybt=temp2(k,i,j)
             yct=temp2(k,i+1,j)
#endif
#ifdef key_substance
             yb(nv_grav+1:nv_adv)=c2(nv_grav+1:nv_adv,k,i,j)
             yc(nv_grav+1:nv_adv)=c2(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=sal2(k,i-1,j)
               yat=temp2(k,i-1,j)
#endif
#ifdef key_substance
               ya(nv_grav+1:nv_adv)=c2(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
             deltas=ycs-ybs
             deltat=yct-ybt
#endif
#ifdef key_substance
             dyg(nv_grav+1:nv_adv)=yb(nv_grav+1:nv_adv)-ya(nv_grav+1:nv_adv)
             dyd(nv_grav+1:nv_adv)=yc(nv_grav+1:nv_adv)-yb(nv_grav+1:nv_adv)
             delta(nv_grav+1:nv_adv)=yc(nv_grav+1:nv_adv)-yb(nv_grav+1:nv_adv)
#endif

           END IF    ! ends test on velocity sign

#if ! defined key_offline

   ! advection flux of salinity
   ! --------------------------

           IF ((ABS(dyds).LE.dsmall).OR.(ABS(dygs).LE.dsmall)) THEN
             ! no horizontal gradient
             thetap = 0.0_rsh
             psip = 0.0_rsh
           ELSE
             ! ultimate limiter
             thetap=(dygs+dsmall)/(dyds+dsmall)
             psip = d0+d1*thetap
             psip=MAX(0.0_rsh,MIN( MIN(1.0_rsh,psip),((1.0_rsh-courant)/(courant+dsmall))*thetap) )
           END IF
           fsalx(k,i,j)  = uu*(ybs+psip*dyds)

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

           IF ((ABS(dydt).LE.dsmall).OR.(ABS(dygt).LE.dsmall)) THEN
             ! no horizontal gradient
             thetap = 0.0_rsh
             psip = 0.0_rsh
           ELSE
             ! ultimate limiter
             thetap=(dygt+dsmall)/(dydt+dsmall)
             psip = d0+d1*thetap
             psip=MAX(0.0_rsh,MIN(MIN(1.0_rsh,psip),((1.0_rsh-courant)/(courant+dsmall))*thetap))
           END IF
           ftempx(k,i,j)  = uu*(ybt+psip*dydt)

   ! key_offline
#endif

#ifdef key_substance

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

           DO iv=nv_grav+1,nv_adv
             IF ((ABS(dyd(iv)).LE.dsmall).OR.(ABS(dyg(iv)).LE.dsmall)) THEN
               ! no horizontal gradient
               thetap = 0.0_rsh
               psip = 0.0_rsh
             ELSE
               ! ultimate limiter
               thetap=(dyg(iv)+dsmall)/(dyd(iv)+dsmall)
               psip = d0+d1*thetap
               psip=MAX(0.0_rsh,MIN(MIN(1.0_rsh,psip),((1.0_rsh-courant)/(courant+dsmall))*thetap))
             END IF
             fx(iv,k,i,j) = uu*(yb(iv)+psip*dyd(iv))
           END DO
#endif


#if ! defined key_offline
#if defined key_ts_diffh_geo || defined key_ts_diffh_rho
   ! Advection term (iso-sigma) alone
   ! --------------------------------
           ftempx(k,i,j) = ftempx(k,i,j)*dy_u(i,j)
           fsalx(k,i,j) = fsalx(k,i,j)*dy_u(i,j)
#else
   ! Advection/diffusion term (iso-sigma)
   ! ------------------------------------
#if defined key_siggen || defined key_gencoord
           ftempx(k,i,j) = ftempx(k,i,j)-hzex(k,i,j)*kxb*deltat/dx_u(i,j)*halfdt
           fsalx(k,i,j) = fsalx(k,i,j)-hzex(k,i,j)*kxb*deltas/dx_u(i,j)*halfdt            
#else
           ftempx(k,i,j) = ftempx(k,i,j)-hex(i,j)*kxb*deltat/dx_u(i,j)*halfdt
           fsalx(k,i,j) = fsalx(k,i,j)-hex(i,j)*kxb*deltas/dx_u(i,j)*halfdt            
#endif
#endif
#endif

#ifdef key_substance
           DO iv=nv_grav+1,nv_adv
#if defined key_ts_diffh_geo || defined key_ts_diffh_rho
   ! Advection term (iso-sigma) alone
   ! --------------------------------
             fx(iv,k,i,j) = fx(iv,k,i,j)*dy_u(i,j)
#else
   ! Advection/diffusion term (iso-sigma)
   ! ------------------------------------
#if defined key_siggen || defined key_gencoord
             fx(iv,k,i,j) = fx(iv,k,i,j) - hzex(k,i,j)*kxb*delta(iv)/dx_u(i,j)*halfdt
#else
             fx(iv,k,i,j) = fx(iv,k,i,j) - hex(i,j)*kxb*delta(iv)/dx_u(i,j)*halfdt
#endif
#endif
#ifdef key_sedim_MUSTANG
             IF (k == 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_uqm
