  SUBROUTINE tssub_advdifx_uqm(sal,temp   &

                                       ,c &

                                         )

!&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---------------------------------------------------------------------





























































!! * 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

   USE comfluxdo,    ONLY : fsalx,fxeau,fyeau
   USE cometemp,     ONLY : ftempx
   USE comvarp2d,    ONLY : sal2,temp2

   USE comsubstance, ONLY : nv_adv,nv_grav,nv_tot,fx,c2

   USE comMUSTANG,    ONLY : corflux


   USE 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

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


!! * 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

   INTEGER                              :: iv
   REAL(KIND=rsh), DIMENSION(nv_adv)    :: dyd,dyg,ya,yb,yc,delta
   REAL(KIND=rsh), DIMENSION(1:nv_adv,kmax) :: cb,ct


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

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

! at open boundaries (east and west)

!$OMP SINGLE

   IF (limin==imin) THEN

     temp2(:,imin:imin+1,ljmin:ljmax) = temp(:,imin:imin+1,ljmin:ljmax)
     sal2(:,imin:imin+1,ljmin:ljmax) = sal(:,imin:imin+1,ljmin:ljmax)


!refonte subtances advestion avril 2016
     c2(1:nv_adv,:,imin:imin+1,ljmin:ljmax) = c(1:nv_adv,:,imin:imin+1,ljmin:ljmax)

   ENDIF
   IF (limax==imax) THEN

     temp2(:,imax,ljmin:ljmax) = temp(:,imax,ljmin:ljmax)
     sal2(:,imax,ljmin:ljmax) = sal(:,imax,ljmin:ljmax)


     c2(nv_grav+1:nv_adv,:,imax,ljmin:ljmax) = c(nv_grav+1:nv_adv,:,imax,ljmin:ljmax)

   ENDIF

!$OMP END SINGLE

! inside the domain

!$OMP DO SCHEDULE(RUNTIME),PRIVATE(i,j,k,     &
!$OMP& salt,tempt,fvt,salb,tempb,fvb,         &

!$OMP& iv,ct,cb,                              &

!$OMP& vv)

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


         IF (j == jb(i) .OR. j==ljminp1m1) THEN


           IF (hty(i,j) <= 0.0_rsh .OR. hy(i,j)==-valmanq-h0fond) THEN

             salt(k)=0.0_rsh       ! no flux advection from this cell
             tempt(k)=0.0_rsh


             ct(nv_grav+1:nv_adv,k)=0.0_rsh

             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

               salt(k)=sal(k,i,j+1)
               tempt(k)=temp(k,i,j+1)


               ct(nv_grav+1:nv_adv,k)=c(nv_grav+1:nv_adv,k,i,j+1)

             ELSE

               salt(k)=sal(k,i,j)
               tempt(k)=temp(k,i,j)


               ct(nv_grav+1:nv_adv,k)=c(nv_grav+1:nv_adv,k,i,j)

             END IF   ! ends test on the direction of propagation

           END IF   ! ends test on the water depth


         ELSE IF (j >= jb(i)+1 .AND. j <= jh(i)-1 .AND. j>ljminp1m1) THEN


! bottom of the stencil

           salb(k)=salt(k)
           tempb(k)=tempt(k)


           cb(nv_grav+1:nv_adv,k)=ct(nv_grav+1:nv_adv,k)

           fvb(k)=fvt(k)

! right of the stencil
           IF (hty(i,j) <= 0.0_rsh .OR. hy(i,j)==-valmanq-h0fond) THEN

             salt(k)=0.0_rsh       ! no flux advection from this cell
             tempt(k)=0.0_rsh


             ct(nv_grav+1:nv_adv,k)=0.0_rsh

             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

               salt(k)=sal(k,i,j+1)
               tempt(k)=temp(k,i,j+1)


               ct(nv_grav+1:nv_adv,k)=c(nv_grav+1:nv_adv,k,i,j+1)

             ELSE

               salt(k)=sal(k,i,j)
               tempt(k)=temp(k,i,j)


               ct(nv_grav+1:nv_adv,k)=c(nv_grav+1:nv_adv,k,i,j)

             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 ?

             vv = fvb(k)/(halfdt*hey(i,j-1))

           ELSE IF ( (fvb(k)<0.0_rlg) .AND. (fvt(k)<0.0_rlg) ) THEN

             vv = fvt(k)/(halfdt*hey(i,j))

           ELSE
! no time evolution of the intermediate temperature
! along the boundaries
             vv = 0.0_rsh
           END IF

! apply euler (time) - upwind (x) scheme

           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))


           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))


         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


     CALL ex_peri_x_rsh(1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
     CALL ex_peri_x_rsh(1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))

     CALL ex_peri_x_rsh(1,nv_adv*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,c2(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2))



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

!$OMP barrier
!$OMP master

   CALL ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))


! to avoid temporary copy   CALL ex_i_rsh(-1,2,(nv_adv)*kmax,liminm1,limaxp2,ljmin,ljmax,c2(1:nv_adv,:,liminm1:limaxp2,ljmin:ljmax))
   CALL ex_i_rsh(-1,2,nv_tot*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,c2(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))

!$OMP end master
!$OMP barrier

!$OMP flush(sal2,temp2)


!$OMP flush(c2)


! 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        &

!$OMP& ,iv,dyd,dyg,ya,yb,yc,delta                          &

!$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

           fsalx(k,i,j)=0.0_rlg
           ftempx(k,i,j)=0.0_rlg


!refonte subtances advection avril 2016
           fx(nv_grav+1:nv_adv,k,i,j)=0.0_rsh

         ELSE
           kxb=kx
           IF (htx(i,j) < hminkxky) kxb=0.0_rsh ! no diffusion for small water depth
           uu=fxeau(k,i,j)

           uuu = uu/(halfdt*hex(i,j))

           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

             ybs=sal2(k,i+1,j)
             ycs=sal2(k,i,j)
             ybt=temp2(k,i+1,j)
             yct=temp2(k,i,j)


             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)

             IF (htx(i+1,j) > 0.0_rsh .AND. (i+1 < id(j) .OR. l_obc_cycl_x)) THEN

               yas=sal2(k,i+2,j)
               yat=temp2(k,i+2,j)


               ya(nv_grav+1:nv_adv)=c2(nv_grav+1:nv_adv,k,i+2,j)

             ELSE

               yas=ybs
               yat=ybt


               ya(nv_grav+1:nv_adv)=yb(nv_grav+1:nv_adv)

             END IF

             dyds=ycs-ybs
             dygs=ybs-yas
             dydt=yct-ybt
             dygt=ybt-yat
             deltas=ybs-ycs
             deltat=ybt-yct


             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)


           ELSE  ! eastward advection


             ybs=sal2(k,i,j)
             ycs=sal2(k,i+1,j)
             ybt=temp2(k,i,j)
             yct=temp2(k,i+1,j)


             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)

             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

               yas=sal2(k,i-1,j)
               yat=temp2(k,i-1,j)


               ya(nv_grav+1:nv_adv)=c2(nv_grav+1:nv_adv,k,i-1,j)

             ELSE

               yas=ybs
               yat=ybt


               ya(nv_grav+1:nv_adv)=yb(nv_grav+1:nv_adv)

             END IF

             dyds=ycs-ybs
             dygs=ybs-yas
             dydt=yct-ybt
             dygt=ybt-yat
             deltas=ycs-ybs
             deltat=yct-ybt


             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)


           END IF    ! ends test on velocity sign



! 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




! 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





! Advection/diffusion term (iso-sigma)
! ------------------------------------

           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            





           DO iv=nv_grav+1,nv_adv

! Advection/diffusion term (iso-sigma)
! ------------------------------------

             fx(iv,k,i,j) = fx(iv,k,i,j) - hex(i,j)*kxb*delta(iv)/dx_u(i,j)*halfdt



             IF (k == 1) fx(iv,k,i,j)=fx(iv,k,i,j)*corflux(iv,i,j)

           END DO


         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

!$OMP barrier
!$OMP master

   CALL ex_i_rlg(-1,0,kmax,liminm1,limax,ljmin,ljmax,fsalx(:,liminm1:limax,ljmin:ljmax))
   CALL ex_i_rlg(-1,0,kmax,liminm1,limax,ljmin,ljmax,ftempx(:,liminm1:limax,ljmin:ljmax))


   CALL ex_i_rsh(-1,0,nv_adv*kmax,liminm1,limax,ljmin,ljmax,fx(:,:,liminm1:limax,ljmin:ljmax))

!$OMP end master
!$OMP barrier

!$OMP flush(fsalx,ftempx)


!$OMP flush(fx)


  END SUBROUTINE tssub_advdifx_uqm
