  SUBROUTINE tssub_advdifx_q(sal,temp   &

                                      ,c &

                                        )

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE tssub_advdifx_quickest  ***
!&E
!&E ** Purpose : Estimate of zonal fluxes of temperature and salinity
!&E              and other substances (advection and diffusion)
!&E              Advection uses QUICKEST scheme.
!&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,fxeau,sal,temp,c(cv_wat),hzex,dx_u
!&E                     fsalx,ftempx,fx
!&E
!&E ** Modified variables : fsalx,ftempx,fx
!&E
!&E ** Reference : Debreu : Caspasthomas.....
!&E
!&E ** History :
!&E       !  2010-12 (V. Garnier) Original code
!&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,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
   INTEGER                   :: i,j,k
   REAL(KIND=rsh), PARAMETER :: dsmall = 1.0d-20
   REAL(KIND=rsh)            :: uu,da,db,dc,uuu,cfl,     &
                                kxb,yas,yat,ybs,ybt,ycs,yct,        &
                                deltat,deltas

   INTEGER                              :: iv
   REAL(KIND=rsh), DIMENSION(nv_adv)    :: ya,yb,yc,delta


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

!$OMP DO schedule(runtime),private(i,j,k,kxb,  &
!$OMP& yas,yat,ybs,ybt,ycs,yct,deltat,deltas,  &

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

!$OMP& da,db,dc,uu,uuu,cfl)

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

           cfl = uuu*halfdt/dx_u(i,j)

           IF (uu <= 0.0_rsh) THEN  ! westward advection
             da=(1.0_rsh/6.0_rsh)*(-1.0_rsh+cfl**2)
             db=(1.0_rsh/6.0_rsh)*(5.0_rsh-3.0_rsh*cfl-2.0_rsh*cfl**2)
             dc=(1.0_rsh/6.0_rsh)*(2.0_rsh+3.0_rsh*cfl+cfl**2)

             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

             deltas=ybs-ycs
             deltat=ybt-yct


             delta(nv_grav+1:nv_adv)=yb(nv_grav+1:nv_adv)-yc(nv_grav+1:nv_adv)


           ELSE  ! eastward advection

             da=(1.0_rsh/6.0_rsh)*(-1.0_rsh+cfl**2)
             db=(1.0_rsh/6.0_rsh)*(5.0_rsh+3.0_rsh*cfl-2.0_rsh*cfl**2)
             dc=(1.0_rsh/6.0_rsh)*(2.0_rsh-3.0_rsh*cfl+cfl**2)

             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

             deltas=ycs-ybs
             deltat=yct-ybt


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

           fsalx(k,i,j)  = uu*(da*yas+db*ybs+dc*ycs)

! flux of temperature advection
! -----------------------------

           ftempx(k,i,j)  = uu*(da*yat+db*ybt+dc*yct)

! key_offline




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

           DO iv=nv_grav+1,nv_adv
             fx(iv,k,i,j) = uu*(da*ya(iv)+db*yb(iv)+dc*yc(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_q
