  SUBROUTINE tssub_advdifx_quickupwind(sal,temp &

                                       ,c       &

                                      )

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





























































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

   USE comfluxdo,    ONLY : fsalx,fxeau
   USE cometemp,     ONLY : ftempx

   USE comsubstance, ONLY : nv_adv,nv_grav,fx

   USE comMUSTANG,    ONLY : corflux


   USE toolmpi,      ONLY : ex_i_rsh,ex_i_rlg

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

   LOGICAL, DIMENSION(nv_adv)           :: decent
   INTEGER                              :: iv
   REAL(KIND=rsh)                       :: dyg
   REAL(KIND=rsh), DIMENSION(nv_adv)    :: ya,yb,yc,dyd


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

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

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

           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)
           IF (uu <= 0.0_rsh) THEN  ! westward advection

             ybs=sal(k,i+1,j)
             ycs=sal(k,i,j)
             ybt=temp(k,i+1,j)
             yct=temp(k,i,j)


             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)

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

               yas=sal(k,i+2,j)
               yat=temp(k,i+2,j)


               ya(nv_grav+1:nv_adv)=c(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=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



             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


           ELSE  ! eastward advection


             ybs=sal(k,i,j)
             ycs=sal(k,i+1,j)
             ybt=temp(k,i,j)
             yct=temp(k,i+1,j)


             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)

             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=sal(k,i-1,j)
               yat=temp(k,i-1,j)


               ya(nv_grav+1:nv_adv)=c(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

             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



             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

           END IF    ! ends test on velocity sign






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

           IF (decents) THEN

             fsalx(k,i,j)=uu*ybs-hex(i,j)*kxb*dyds/dx_u(i,j)*halfdt

           ELSE
             yr=0.375_rsh*ycs+0.75_rsh*ybs-0.125_rsh*yas

             fsalx(k,i,j)=uu*yr-hex(i,j)*kxb*dyds/dx_u(i,j)*halfdt

           END IF

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

           IF (decentt) THEN

             ftempx(k,i,j)=uu*ybt-hex(i,j)*kxb*dydt/dx_u(i,j)*halfdt

           ELSE
             yr=0.375_rsh*yct+0.75_rsh*ybt-0.125_rsh*yat

             ftempx(k,i,j)=uu*yr -hex(i,j)*kxb*dydt/dx_u(i,j)*halfdt

           END IF



! key_offline




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

           DO iv=nv_grav+1,nv_adv


             IF (decent(iv)) THEN

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

             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 - hex(i,j)*kxb*dyd(iv)/dx_u(i,j)*halfdt

             END IF


               IF (k.eq.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_quickupwind
