  SUBROUTINE tssub_advdify_quickupwind(sal,temp &

                                       ,c      &

                                       )

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE tssub_advdify_quickupwind  ***
!&E
!&E ** Purpose : Estimate of meridional 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                  - 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 : 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 southern 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, include sand)
!&E           !  2016-12  (B.Thouvenin) version sedim_MUSTANG : only mixsed (mudsand not operationnal)
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d,    ONLY : jb,jh,hty,hey,hy,h0fond,qmax,ky,hminkxky,dy_v,halfdt,l_obc_cycl_y,dx_v

   USE comfluxdo,    ONLY : fsaly,fyeau
   USE cometemp,     ONLY  : ftempy

   USE comsubstance, ONLY : nv_adv,nv_grav,fy

   USE comMUSTANG,    ONLY : corfluy


   USE toolmpi,      ONLY : ex_j_rsh,ex_j_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,vv,yr,qot,dyds,dygs,dydt
   REAL(KIND=rsh)          :: dygt,kyb,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,kyb               &
!$OMP& ,vv,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=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

             fsaly(k,i,j)=0.0_rlg
             ftempy(k,i,j)=0.0_rlg


             fy(nv_grav+1:nv_adv,k,i,j)=0.0_rsh

           ELSE
             kyb=ky
             IF (hty(i,j) < hminkxky) kyb=0.0_rsh ! no diffusion for small water depth
             vv=fyeau(k,i,j)
             IF (vv <= 0.0_rsh) THEN  ! southward advection

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


               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)


               IF (hty(i,j+1) > 0.0_rsh .AND. (j+1 < jh(i) .OR. l_obc_cycl_y)) THEN

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


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

               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      ! northward advection


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


               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)

               IF (hty(i,j-1) > 0.0_rsh .AND. (j > jb(i) .OR. l_obc_cycl_y)) THEN

                 yas=sal(k,i,j-1)
                 yat=temp(k,i,j-1)


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

               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

               fsaly(k,i,j)=vv*ybs-hey(i,j)*kyb*dyds/dy_v(i,j)*halfdt

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

               fsaly(k,i,j)=vv*yr-hey(i,j)*kyb*dyds/dy_v(i,j)*halfdt

             END IF

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

             IF (decentt) THEN

               ftempy(k,i,j)=vv*ybt-hey(i,j)*kyb*dydt/dy_v(i,j)*halfdt

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

               ftempy(k,i,j)=vv*yr -hey(i,j)*kyb*dydt/dy_v(i,j)*halfdt

             END IF



! key_offline




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

             DO iv=nv_grav+1,nv_adv


               IF (decent(iv)) THEN

                 fy(iv,k,i,j)=fyeau(k,i,j)*yb(iv)-hey(i,j)*kyb*dyd(iv)/dy_v(i,j)*halfdt

               ELSE
                 yr=0.375_rsh*yc(iv)+0.75_rsh*yb(iv)-0.125_rsh*ya(iv)

                 fy(iv,k,i,j)=fyeau(k,i,j)*yr-hey(i,j)*kyb*dyd(iv)/dy_v(i,j)*halfdt

               END IF



               IF (k.eq.1) fy(iv,k,i,j)=fy(iv,k,i,j)*corfluy(iv,i,j)

             END DO


           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

!$OMP barrier
!$OMP master

   CALL ex_j_rlg(-1,0,kmax,limin,limax,ljminm1,ljmax,fsaly(:,limin:limax,ljminm1:ljmax))
   CALL ex_j_rlg(-1,0,kmax,limin,limax,ljminm1,ljmax,ftempy(:,limin:limax,ljminm1:ljmax))


   CALL ex_j_rsh(-1,0,nv_adv*kmax,limin,limax,ljminm1,ljmax,fy(:,:,limin:limax,ljminm1:ljmax))

!$OMP end master
!$OMP barrier

!$OMP flush(fsaly,ftempy)


!$OMP flush(fy)


  END SUBROUTINE tssub_advdify_quickupwind
