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

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE tssub_advdify_ultimatequickestmacho  ***
   !&E
   !&E ** Purpose : Estimate of meridional 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 southward 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),hzey,dx,dy_v
   !&E                     fsaly,ftempy,fy
   !&E
   !&E ** Modified variables : fsaly,ftempy,fy
   !&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 southern boundary
   !&E           !  2008-10 (V. Garnier) optimization
   !&E           !  2009-04 (A. Mano) MPI2D
   !&E           !  2010-05 (V. Garnier) OMP
   !&E           !  2010-10 (T. Duhaut) introduction of diffusion along isopycnes
   !&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,dx,htx,hex,    &
                            ig,id,hx,halfdt,l_obc_cycl_y,dx_v
#if defined key_siggen || defined key_gencoord
   USE comsiggen,    ONLY : hzex,hzey
#endif
   USE comfluxdo,    ONLY : fsaly,fyeau,fxeau
   USE cometemp,     ONLY : ftempy
   USE comvarp2d,    ONLY : sal2,temp2
#ifdef key_substance
   USE comsubstance, ONLY : nv_adv,nv_grav,nv_tot,fy,c2
#ifdef key_sedim_MUSTANG
   USE comMUSTANG,    ONLY : corfluy
#endif
#endif
   USE_MPI toolmpi,      ONLY : ex_j_rlg,ex_j_rsh,ex_peri_y_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)            :: vv,dyds,dygs,dydt,d0,d1,vvv,courant,uu,  &
                                dygt,kyb,yas,yat,ybs,ybt,ycs,yct,thetap, &
                                psip,deltat,deltas
   REAL(KIND=rlg), DIMENSION(kmax) :: ful,fur
   REAL(KIND=rsh), DIMENSION(kmax) :: sall,salr,templ,tempr
#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) :: cl,cr
#endif

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

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

   ! at open boundaries (north and south)

!$OMP SINGLE

   IF_MPI (ljmin==jmin) THEN
#if ! defined key_offline
     temp2(:,limin:limax,jmin:jmin+1) = temp(:,limin:limax,jmin:jmin+1)
     sal2(:,limin:limax,jmin:jmin+1) = sal(:,limin:limax,jmin:jmin+1)
#endif
#ifdef key_substance
     c2(nv_grav+1:nv_adv,:,limin:limax,jmin:jmin+1) = c(nv_grav+1:nv_adv,:,limin:limax,jmin:jmin+1)
#endif
   ENDIF_MPI
   IF_MPI (ljmax==jmax) THEN
#if ! defined key_offline
     temp2(:,limin:limax,jmax) = temp(:,limin:limax,jmax)
     sal2(:,limin:limax,jmax) = sal(:,limin:limax,jmax)
#endif
#ifdef key_substance
     c2(nv_grav+1:nv_adv,:,limin:limax,jmax) = c(nv_grav+1:nv_adv,:,limin:limax,jmax)
#endif
   ENDIF_MPI

!$OMP END SINGLE

   ! inside the domain

!$OMP DO SCHEDULE(RUNTIME),PRIVATE(i,j,k,     &
!$OMP& salr,tempr,fur,sall,templ,ful,         &
#ifdef key_substance
!$OMP& iv,cr,cl,                              &
#endif
!$OMP& uu)

   DO j=ljmin,ljmax
#ifdef key_MPI_2D
     DO i=liminp1m1,limax
     IF (i==ig(j) .OR. i==liminp1m1) THEN
#else
     i=ig(j)
#endif

       DO k=1,kmax
         IF (htx(i,j) <= 0.0_rsh .OR. hx(i,j) == -valmanq-h0fond) THEN
#if ! defined key_offline
           salr(k)=0.0_rsh       ! no flux advection from this cell
           tempr(k)=0.0_rsh
#endif
#ifdef key_substance
           cr(nv_grav+1:nv_adv,k)=0.0_rsh
#endif
           fur(k)=0.0_rlg
         ELSE

           fur(k)=fxeau(k,i,j)   ! h*uz*dt/2

           ! upwind scheme
           IF (fur(k)<= 0.0_rlg) THEN
#if ! defined key_offline
             salr(k)=sal(k,i+1,j)
             tempr(k)=temp(k,i+1,j)
#endif
#ifdef key_substance
             cr(nv_grav+1:nv_adv,k)=c(nv_grav+1:nv_adv,k,i+1,j)
#endif
           ELSE
#if ! defined key_offline
             salr(k)=sal(k,i,j)
             tempr(k)=temp(k,i,j)
#endif
#ifdef key_substance
             cr(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
       END DO   ! ends loop on k

#ifdef key_MPI_2D
     ELSE IF(i >= ig(j)+1 .AND. i <=id(j)-1 .AND. i>liminp1m1) THEN
#else
     DO i=ig(j)+1,id(j)-1
#endif

       DO k=1,kmax

         ! left of the stencil
#if ! defined key_offline
         sall(k)=salr(k)
         templ(k)=tempr(k)
#endif
#ifdef key_substance
         cl(nv_grav+1:nv_adv,k)=cr(nv_grav+1:nv_adv,k)
#endif
         ful(k)=fur(k)

         ! right of the stencil
         IF (htx(i,j) <= 0.0_rsh .OR. hx(i,j) == -valmanq-h0fond) THEN
#if ! defined key_offline
           salr(k)=0.0_rsh       ! no flux advection from this cell
           tempr(k)=0.0_rsh
#endif
#ifdef key_substance
           cr(nv_grav+1:nv_adv,k)=0.0_rsh
#endif
           fur(k)=0.0_rlg
         ELSE

           fur(k)=fxeau(k,i,j)   ! h*uz*dt/2

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

         END IF   ! ends test on the water depth

     ! estimate of intermediate tracers (from upwind scheme)

         IF ( (ful(k)>0.0_rlg) .AND. (fur(k)>0.0_rlg) ) THEN
#if defined key_siggen || defined key_gencoord
           uu = ful(k)/(halfdt*hzex(k,i-1,j))
#else
           uu = ful(k)/(halfdt*hex(i-1,j))
#endif
         ELSE IF ( (ful(k)<0.0_rlg) .AND. (fur(k)<0.0_rlg) ) THEN
#if defined key_siggen || defined key_gencoord
           uu = fur(k)/(halfdt*hzex(k,i,j))
#else
           uu = fur(k)/(halfdt*hex(i,j))
#endif
         ELSE
           ! no time evolution of the intermediate temperature
           ! along the boundaries
           uu = 0.0_rsh
         END IF

         ! apply euler (time) - upwind (x) scheme
#if ! defined key_offline
         temp2(k,i,j) = temp(k,i,j) - (halfdt/dx(i,j))*uu*(tempr(k)-templ(k))
         sal2(k,i,j)  = sal(k,i,j)  - (halfdt/dx(i,j))*uu*(salr(k)-sall(k))  
#endif
#ifdef key_substance
         c2(nv_grav+1:nv_adv,k,i,j)  = c(nv_grav+1:nv_adv,k,i,j)  &
           - (halfdt/dx(i,j))*uu*(cr(nv_grav+1:nv_adv,k)-cl(nv_grav+1:nv_adv,k))  
#endif

       END DO   ! ends loop on k
#ifdef key_MPI_2D
     END IF
#endif
     END DO   ! ends loop on i
   END DO   ! ends loop on j
!$OMP END DO

!$OMP SINGLE
   ! periodic conditions along y-axis
   IF (l_obc_cycl_y) THEN

#ifdef key_MPI_2D
#ifdef key_shiftobc
     CALL_MPI ex_peri_y_rsh(3,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
     CALL_MPI ex_peri_y_rsh(3,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#ifdef key_substance
     CALL_MPI ex_peri_y_rsh(3,(nv_adv)*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,c2(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#else
     CALL_MPI ex_peri_y_rsh(1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
     CALL_MPI ex_peri_y_rsh(1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#ifdef key_substance
     CALL_MPI ex_peri_y_rsh(1,(nv_adv)*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,c2(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#endif
#else
     DO i=limin,limax
   ! southern boundary
       sal2(:,i,jb(i)) = sal2(:,i,jh(i)-1)
       sal2(:,i,jb(i)-1) = sal2(:,i,jh(i)-2)
       temp2(:,i,jb(i)) = temp2(:,i,jh(i)-1)
       temp2(:,i,jb(i)-1) = temp2(:,i,jh(i)-2)
#ifdef key_substance
       c2(nv_grav+1:nv_adv,:,i,jb(i))=c2(nv_grav+1:nv_adv,:,i,jh(i)-1)
       c2(nv_grav+1:nv_adv,:,i,jb(i)-1)=c2(nv_grav+1:nv_adv,:,i,jh(i)-2)
#endif
#ifdef key_shiftobc
       sal2(:,i,jb(i)-2) = sal2(:,i,jh(i)-3)
       temp2(:,i,jb(i)-2) = temp2(:,i,jh(i)-3)
#ifdef key_substance
       c2(nv_grav+1:nv_adv,:,i,jb(i)-2)=c2(nv_grav+1:nv_adv,:,i,jh(i)-3)
#endif
#endif
   ! northern boundary
       sal2(:,i,jh(i)) = sal2(:,i,jb(i)+1)
       temp2(:,i,jh(i)) = temp2(:,i,jb(i)+1)
#ifdef key_substance
       c2(nv_grav+1:nv_adv,:,i,jh(i))=c2(nv_grav+1:nv_adv,:,i,jb(i)+1)
#endif
#ifdef key_shiftobc
       sal2(:,i,jh(i)+1) = sal2(:,i,jb(i)+2)
       sal2(:,i,jh(i)+2) = sal2(:,i,jb(i)+3)
       temp2(:,i,jh(i)+1) = temp2(:,i,jb(i)+2)
       temp2(:,i,jh(i)+2) = temp2(:,i,jb(i)+3)
#ifdef key_substance
       c2(nv_grav+1:nv_adv,:,i,jh(i)+1)=c2(nv_grav+1:nv_adv,:,i,jb(i)+2)
       c2(nv_grav+1:nv_adv,:,i,jh(i)+2)=c2(nv_grav+1:nv_adv,:,i,jb(i)+3)
#endif
#endif
       END DO
#endif

     END IF  ! ends periodic conditions along y-axis
!$OMP END SINGLE

OMPMPI barrier
OMPMPI master
#if ! defined key_offline
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#ifdef key_substance
!to avoid temporary copy 1st dimension is not contiguous   CALL_MPI ex_j_rsh(-1,2,(nv_adv)*kmax,limin,limax,ljminm1,ljmaxp2,c2(1:nv_adv,:,limin:limax,ljminm1:ljmaxp2))
   CALL_MPI ex_j_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,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=sal2(k,i,j+1)
               ycs=sal2(k,i,j)
               ybt=temp2(k,i,j+1)
               yct=temp2(k,i,j)
#endif
#ifdef key_substance
               yb(nv_grav+1:nv_adv)=c2(nv_grav+1:nv_adv,k,i,j+1)
               yc(nv_grav+1:nv_adv)=c2(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=sal2(k,i,j+2)
                 yat=temp2(k,i,j+2)
#endif
#ifdef key_substance
                 ya(nv_grav+1:nv_adv)=c2(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=sal2(k,i,j)
               ycs=sal2(k,i,j+1)
               ybt=temp2(k,i,j)
               yct=temp2(k,i,j+1)
#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,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=sal2(k,i,j-1)
                 yat=temp2(k,i,j-1)
#endif
#ifdef key_substance
                 ya(nv_grav+1:nv_adv)=c2(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 (iso-sigma) 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
               fy(iv,k,i,j) = fy(iv,k,i,j)*dx_v(i,j)
#else
#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_uqm
