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

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE tssub_advdify_ultimatequickest  ***
   !&E
   !&E ** Purpose : Estimate of meridional fluxes of temperature and salinity
   !&E              and other substances (advection and diffusion)
   !&E              Advection uses QUICKEST scheme with ULTIMATE limiter.
   !&E
   !&E ** Description : No horizontal diffusion at lower depths ( < hminkxky)
   !&E                 - case of cyclicity in y : correct for northward propagation,
   !&E                    uncorrect for southward propagation (stencil not long enough)
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays : hty,hey,hy,fyeau,sal,temp,c(cv_wat),hzey,dy_v
   !&E                     fsaly,ftempy,fy
   !&E
   !&E ** Modified variables : fsaly,ftempy,fy
   !&E
   !&E ** Reference : Leonard B. P. (1991) : The ULTIMATE conservative
   !&E                difference scheme applied to unsteady one-dimensional
   !&E                advection. Computer Methods in Applied Mechanics and
   !&E                Engineering, 88, 17-74
   !&E
   !&E ** History :
   !&E           !  2008-09 (L. Debreu, T. Duhaut, M. Honnorat) original code
   !&E           !  2008-09 (V. Garnier) generalized sigma, correction for southern boundary
   !&E           !  2009-04 (A. Mano) MPI2D
   !&E           !  2010-10 (T. Duhaut) introduction of diffusion along isopycnes
   !&E           !  2016-04  (B.Thouvenin) applied to all substances (dissolved and particulate, include sand)
   !&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 : jb,jh,hty,hey,hy,h0fond,ky,hminkxky,dy_v,halfdt,l_obc_cycl_y,dx_v
#if defined key_siggen || defined key_gencoord
   USE comsiggen,   ONLY : hzey
#endif
   USE comfluxdo,   ONLY : fsaly,fyeau
   USE cometemp,    ONLY : ftempy
#ifdef key_substance
   USE comsubstance,ONLY : nv_adv,nv_grav,fy
#ifdef key_sedim_MUSTANG
   USE comMUSTANG,    ONLY : corfluy
#endif
#endif
   USE_MPI toolmpi, ONLY : ex_j_rsh,ex_j_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
   INTEGER                   :: i,j,k
   REAL(KIND=rsh), PARAMETER :: dsmall = 1.0d-20
   REAL(KIND=rsh)            :: vv,dyds,dygs,dydt,d0,d1,vvv,courant,     &
                                dygt,kyb,yas,yat,ybs,ybt,ycs,yct,thetap, &
                                psip,deltat,deltas
#ifdef key_substance
   INTEGER                              :: iv
   REAL(KIND=rsh), DIMENSION(nv_adv)    :: dyd,dyg,ya,yb,yc,delta
#endif

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

!$OMP DO schedule(runtime),private(i,j,k,dyds,dygs,dydt    &
!$OMP&,dygt,kyb,yas,yat,ybs,ybt,ycs,yct,deltat,deltas      &
#ifdef key_substance
!$OMP& ,iv,dyd,dyg,ya,yb,yc,delta                          &
#endif
!$OMP&,d0,d1,vv,vvv,courant,thetap,psip)

   DO j=ljmin,ljmax
     DO i=liminp2,limaxm1

       IF (j >= jb(i) .AND. j <= jh(i)-1) THEN
         DO k=1,kmax
           IF (hty(i,j) <= 0.0_rsh .OR. hy(i,j)==-valmanq-h0fond) THEN
#if ! defined key_offline
             fsaly(k,i,j)=0.0_rlg
             ftempy(k,i,j)=0.0_rlg
#endif
#ifdef key_substance
             fy(nv_grav+1:nv_adv,k,i,j)=0.0_rsh
#endif
           ELSE
             kyb=ky
             IF (hty(i,j) < hminkxky) kyb=0.0_rsh ! no diffusion for small water depth
             vv=fyeau(k,i,j)
#if defined key_siggen || defined key_gencoord
             vvv = vv/(halfdt*hzey(k,i,j))
#else
             vvv = vv/(halfdt*hey(i,j))
#endif
             courant = ABS(vvv)*halfdt/dy_v(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 (vv <= 0.0_rsh) THEN  ! southward advection
#if ! defined key_offline
               ybs=sal(k,i,j+1)
               ycs=sal(k,i,j)
               ybt=temp(k,i,j+1)
               yct=temp(k,i,j)
#endif
#ifdef key_substance
               yb(nv_grav+1:nv_adv)=c(nv_grav+1:nv_adv,k,i,j+1)
               yc(nv_grav+1:nv_adv)=c(nv_grav+1:nv_adv,k,i,j)
#endif

               IF (hty(i,j+1) > 0.0_rsh .AND. (j+1 < jh(i) .OR. l_obc_cycl_y)) THEN
#if ! defined key_offline
                 yas=sal(k,i,j+2)
                 yat=temp(k,i,j+2)
#endif
#ifdef key_substance
                 ya(nv_grav+1:nv_adv)=c(nv_grav+1:nv_adv,k,i,j+2)
#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
               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)=yb(nv_grav+1:nv_adv)-yc(nv_grav+1:nv_adv)
#endif

             ELSE      ! northward advection

#if ! defined key_offline
               ybs=sal(k,i,j)
               ycs=sal(k,i,j+1)
               ybt=temp(k,i,j)
               yct=temp(k,i,j+1)
#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,j+1)
#endif
               IF (hty(i,j-1) > 0.0_rsh .AND. (j > jb(i) .OR. l_obc_cycl_y)) THEN
#if ! defined key_offline
                 yas=sal(k,i,j-1)
                 yat=temp(k,i,j-1)
#endif
#ifdef key_substance
                 ya(nv_grav+1:nv_adv)=c(nv_grav+1:nv_adv,k,i,j-1)
#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

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

             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
             fsaly(k,i,j) = vv*(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
             ftempy(k,i,j) = vv*(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
               fy(iv,k,i,j) = vv*(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 alone
   ! --------------------
             ftempy(k,i,j)=ftempy(k,i,j)*dx_v(i,j)
             fsaly(k,i,j)=fsaly(k,i,j)*dx_v(i,j)
#else
   ! Advection/diffusion term (iso-sigma)
   ! ------------------------------------
#if defined key_siggen || defined key_gencoord
             ftempy(k,i,j)=ftempy(k,i,j)-hzey(k,i,j)*kyb*deltat/dy_v(i,j)*halfdt
             fsaly(k,i,j)=fsaly(k,i,j)-hzey(k,i,j)*kyb*deltas/dy_v(i,j)*halfdt              
#else
             ftempy(k,i,j)=ftempy(k,i,j)-hey(i,j)*kyb*deltat/dy_v(i,j)*halfdt
             fsaly(k,i,j)=fsaly(k,i,j)-hey(i,j)*kyb*deltas/dy_v(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 alone
   ! --------------------
               fy(iv,k,i,j) = fy(iv,k,i,j)*dx_v(i,j)
#else
   ! Advection/diffusion term (iso-sigma)
   ! ------------------------------------
#if defined key_siggen || defined key_gencoord
               fy(iv,k,i,j) = fy(iv,k,i,j) - hzey(k,i,j)*kyb*delta(iv)/dy_v(i,j)*halfdt
#else
               fy(iv,k,i,j) = fy(iv,k,i,j) - hey(i,j)*kyb*delta(iv)/dy_v(i,j)*halfdt
#endif
#endif
#ifdef key_sedim_MUSTANG
               IF (k.eq.1) fy(iv,k,i,j)=fy(iv,k,i,j)*corfluy(iv,i,j)
#endif
             END DO
#endif

           END IF ! ends test on water depth
         END DO ! ends loop on k
       END IF ! for selecting right j loop
     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_j_rlg(-1,0,kmax,limin,limax,ljminm1,ljmax,fsaly(:,limin:limax,ljminm1:ljmax))
   CALL_MPI ex_j_rlg(-1,0,kmax,limin,limax,ljminm1,ljmax,ftempy(:,limin:limax,ljminm1:ljmax))
#endif
#ifdef key_substance
   CALL_MPI ex_j_rsh(-1,0,nv_adv*kmax,limin,limax,ljminm1,ljmax,fy(:,:,limin:limax,ljminm1:ljmax))
#endif
OMPMPI end master
OMPMPI barrier
#if ! defined key_offline
OMPMPI flush(fsaly,ftempy)
#endif
#ifdef key_substance
OMPMPI flush(fy)
#endif

  END SUBROUTINE tssub_advdify_uq
