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

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE tssub_advdifx_fifthorder  ***
   !&E
   !&E ** Purpose : Estimate of zonal fluxes of temperature and salinity
   !&E              and other substances (advection and diffusion)
   !&E              Advection uses FIFTH order 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 : 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, Duhaut (2011) contrat PREVIMER 09/2 211 097 
   !&E
   !&E ** History :
   !&E           !  2012-07 (V. Garnier) Original code (from tssubadvdifxultimatequickest.F90)
   !&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---------------------------------------------------------------------

#include "toolcpp.h"

#ifdef key_tssub_adv_fifthorder

   !! * Modules used
   USE parameters
   USE comvars2d, ONLY : ig,id,htx,hex,hx,kx,h0fond,hminkxky,dx_u,halfdt,l_obc_cycl_x,dy_u
#if defined key_siggen || defined key_gencoord
   USE comsiggen, ONLY : hzex
#endif
   USE comfluxdo, ONLY : fsalx,fxeau
   USE cometemp,  ONLY : ftempx
#ifdef key_substance
   USE comsubstance, ONLY : nv_adv,nv_grav,fx
   USE comsubstance,ONLY : c5
#ifdef key_sedim_MUSTANG
   USE comMUSTANG,    ONLY : corflux
#endif
#endif
   USE_MPI toolmpi, ONLY : ex_i_rsh,ex_i_rlg
   USE_MPI toolmpi, ONLY : ex_j_rsh
#ifdef key_tssub_adv_fifthorder
   USE comvarp2d, ONLY : sal5,temp5  !VALERIE PROVISOIRE
#endif

   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)            :: uu,da,db,dc,dd,de,uuu,cfl,     &
                                kxb,yas,yat,ybs,ybt,ycs,yct,yds,ydt, &
                                yes,yet,deltat,deltas
#ifdef key_substance
   INTEGER                              :: iv
   REAL(KIND=rsh), DIMENSION(nv_adv)    :: ya,yb,yc,yd,ye,delta
#endif

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

   PRINT_DBG*, 'enter TSSUB_ADVDIFX_FIFTHORDER.F90 routine'

!$OMP SINGLE
   sal5(:,liminm1:limaxp2,ljminm1:ljmaxp2) = sal(:,liminm1:limaxp2,ljminm1:ljmaxp2)
   temp5(:,liminm1:limaxp2,ljminm1:ljmaxp2) = temp(:,liminm1:limaxp2,ljminm1:ljmaxp2)
#ifdef key_substance
   c5(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2) = c(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2)
#endif
!$OMP END SINGLE
OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-2,3,kmax,liminm2,limaxp3,ljminm2,ljmaxp3,sal5(:,liminm2:limaxp3,ljminm2:ljmaxp3))
   CALL_MPI ex_j_rsh(-2,3,kmax,liminm2,limaxp3,ljminm2,ljmaxp3,sal5(:,liminm2:limaxp3,ljminm2:ljmaxp3))
   CALL_MPI ex_i_rsh(-2,3,kmax,liminm2,limaxp3,ljminm2,ljmaxp3,temp5(:,liminm2:limaxp3,ljminm2:ljmaxp3))
   CALL_MPI ex_j_rsh(-2,3,kmax,liminm2,limaxp3,ljminm2,ljmaxp3,temp5(:,liminm2:limaxp3,ljminm2:ljmaxp3))
#ifdef key_substance
   CALL_MPI ex_i_rsh(-2,3,nv_adv*kmax,liminm2,limaxp3,ljminm2,ljmaxp3,c5(:,:,liminm2:limaxp3,ljminm2:ljmaxp3))
   CALL_MPI ex_j_rsh(-2,3,nv_adv*kmax,liminm2,limaxp3,ljminm2,ljmaxp3,c5(:,:,liminm2:limaxp3,ljminm2:ljmaxp3))
#endif
OMPMPI end master
OMPMPI barrier
#if ! defined key_offline
OMPMPI flush(sal5,temp5)
#endif
#ifdef key_substance
OMPMPI flush(c5)
#endif


!$OMP DO schedule(runtime),private(i,j,k,kxb               &
!$OMP& ,yas,yat,ybs,ybt,ycs,yct,yds,ydt,yes,yet        &
#ifdef key_substance
!$OMP& ,iv,ya,yb,yc,yd,ye,delta                          &
#endif
!$OMP& ,deltat,deltas,da,db,dc,dd,de,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
#if ! defined key_offline
           fsalx(k,i,j)=0.0_rlg
           ftempx(k,i,j)=0.0_rlg
#endif
#ifdef key_substance
           fx(nv_grav+1:nv_adv,k,i,j)=0.0_rsh
#endif
         ELSE
           kxb=kx
           IF (htx(i,j) < hminkxky) kxb=0.0_rsh ! no diffusion for small water depth
           uu=fxeau(k,i,j)
#if defined key_siggen || defined key_gencoord
           uuu = uu/(halfdt*hzex(k,i,j))
#else
           uuu = uu/(halfdt*hex(i,j))
#endif
           cfl = uuu*halfdt/dx_u(i,j)

           IF (uu <= 0.0_rsh) THEN  ! westward advection

             de=(1.0_rsh/120.0_rsh)*(4.0_rsh-5.0_rsh*cfl**2+cfl**4)
             dd=(1.0_rsh/120.0_rsh)*(-26.0_rsh+5.0_rsh*cfl+30.0_rsh*cfl**2-5.0_rsh*cfl**3-4.0_rsh*cfl**4)
             dc=(1.0_rsh/120.0_rsh)*(94.0_rsh-75.0_rsh*cfl-40.0_rsh*cfl**2+15.0_rsh*cfl**3+6.0_rsh*cfl**4)
             db=(1.0_rsh/120.0_rsh)*(54.0_rsh+75.0_rsh*cfl+10.0_rsh*cfl**2-15.0_rsh*cfl**3-4.0_rsh*cfl**4)
             da=(1.0_rsh/120.0_rsh)*(-6.0_rsh-5.0_rsh*cfl+5.0_rsh*cfl**2+5.0_rsh*cfl**3+cfl**4)

#if ! defined key_offline
             ycs=sal5(k,i+1,j)
             ybs=sal5(k,i,j)
             yct=temp5(k,i+1,j)
             ybt=temp5(k,i,j)
#endif
#ifdef key_substance
             yc(nv_grav+1:nv_adv)=c5(nv_grav+1:nv_adv,k,i+1,j)
             yb(nv_grav+1:nv_adv)=c5(nv_grav+1:nv_adv,k,i,j)
#endif
             IF (htx(i+1,j) > 0.0_rsh .AND. (i+1 < id(j) .OR. l_obc_cycl_x)) THEN
#if ! defined key_offline
               yds=sal5(k,i+2,j)
               ydt=temp5(k,i+2,j)
#endif
#ifdef key_substance
               yd(nv_grav+1:nv_adv)=c5(nv_grav+1:nv_adv,k,i+2,j)
#endif
             ELSE
#if ! defined key_offline
               yds=ycs
               ydt=yct
#endif
#ifdef key_substance
               yd(nv_grav+1:nv_adv)=yc(nv_grav+1:nv_adv)
#endif
             END IF
             IF (i+2<=imax .AND. htx(MIN0(i+2,imax),j) > 0.0_rsh .AND. (i+2 < id(j) .OR. l_obc_cycl_x)) THEN
#if ! defined key_offline
               yes=sal5(k,i+3,j)
               yet=temp5(k,i+3,j)
#endif
#ifdef key_substance
               ye(nv_grav+1:nv_adv)=c5(nv_grav+1:nv_adv,k,i+3,j)
#endif
             ELSE
#if ! defined key_offline
               yes=yds
               yet=ydt
#endif
#ifdef key_substance
               ye(nv_grav+1:nv_adv)=yd(nv_grav+1:nv_adv)
#endif
             END IF
             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 .AND. i-1 > ig(j)) THEN
#if ! defined key_offline
               yas=sal5(k,i-1,j)
               yat=temp5(k,i-1,j)
#endif
#ifdef key_substance
               ya(nv_grav+1:nv_adv)=c5(nv_grav+1:nv_adv,k,i-1,j)
#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
             deltas=ycs-ybs
             deltat=yct-ybt
#endif
#ifdef key_substance
             delta(nv_grav+1:nv_adv)=yc(nv_grav+1:nv_adv)-yb(nv_grav+1:nv_adv)
#endif


           ELSE  ! eastward advection

             da=(1.0_rsh/120.0_rsh)*(4.0_rsh-5.0_rsh*cfl**2+cfl**4)
             db=(1.0_rsh/120.0_rsh)*(-26.0_rsh-5.0_rsh*cfl+30.0_rsh*cfl**2+5.0_rsh*cfl**3-4.0_rsh*cfl**4)
             dc=(1.0_rsh/120.0_rsh)*(94.0_rsh+75.0_rsh*cfl-40.0_rsh*cfl**2-15.0_rsh*cfl**3+6.0_rsh*cfl**4)
             dd=(1.0_rsh/120.0_rsh)*(54.0_rsh-75.0_rsh*cfl+10.0_rsh*cfl**2+15.0_rsh*cfl**3-4.0_rsh*cfl**4)
             de=(1.0_rsh/120.0_rsh)*(-6.0_rsh+5.0_rsh*cfl+5.0_rsh*cfl**2-5.0_rsh*cfl**3+cfl**4)

#if ! defined key_offline
             ycs=sal5(k,i,j)
             yds=sal5(k,i+1,j)
             yct=temp5(k,i,j)
             ydt=temp5(k,i+1,j)
#endif
#ifdef key_substance
             yc(nv_grav+1:nv_adv)=c5(nv_grav+1:nv_adv,k,i,j)
             yd(nv_grav+1:nv_adv)=c5(nv_grav+1:nv_adv,k,i+1,j)
#endif
             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
#if ! defined key_offline
               ybs=sal5(k,i-1,j)
               ybt=temp5(k,i-1,j)
#endif
#ifdef key_substance
               yb(nv_grav+1:nv_adv)=c5(nv_grav+1:nv_adv,k,i-1,j)
#endif
             ELSE
#if ! defined key_offline
               ybs=ycs
               ybt=yct
#endif
#ifdef key_substance
               yb(nv_grav+1:nv_adv)=yc(nv_grav+1:nv_adv)
#endif
             END IF
             IF (i-2>=0 .AND. htx(MAX0(i-2,imin),j) > 0.0_rsh .AND. (i-1 > ig(j) .OR. l_obc_cycl_x)) THEN
#if ! defined key_offline
               yas=sal5(k,i-2,j)
               yat=temp5(k,i-2,j)
#endif
#ifdef key_substance
               ya(nv_grav+1:nv_adv)=c5(nv_grav+1:nv_adv,k,i-2,j)
#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 (htx(i+1,j) > 0.0_rsh .AND. (i+1 < id(j) .OR. l_obc_cycl_x)) THEN
#if ! defined key_offline
               yes=sal5(k,i+2,j)
               yet=temp5(k,i+2,j)
#endif
#ifdef key_substance
               ye(nv_grav+1:nv_adv)=c5(nv_grav+1:nv_adv,k,i+2,j)
#endif
             ELSE
#if ! defined key_offline
               yes=yds
               yet=ydt
#endif
#ifdef key_substance
               ye(nv_grav+1:nv_adv)=yd(nv_grav+1:nv_adv)
#endif
             END IF

#if ! defined key_offline
             deltas=yds-ycs
             deltat=ydt-yct
#endif
#ifdef key_substance
             delta(nv_grav+1:nv_adv)=yd(nv_grav+1:nv_adv)-yc(nv_grav+1:nv_adv)
#endif

           END IF    ! ends test on velocity sign


   ! advection flux of salinity/temperature/substances
   ! -------------------------------------------------

#if ! defined key_offline
           fsalx(k,i,j)  = uu*(da*yas+db*ybs+dc*ycs+dd*yds+de*yes)
           ftempx(k,i,j)  = uu*(da*yat+db*ybt+dc*yct+dd*ydt+de*yet)
#endif

#ifdef key_substance
           DO iv=nv_grav+1,nv_adv
             fx(iv,k,i,j) = uu*(da*ya(iv)+db*yb(iv)+dc*yc(iv)+dd*yd(iv)+de*ye(iv))
           END DO
#endif

#if ! defined key_offline
#if defined key_ts_diffh_geo || defined key_ts_diffh_rho
   ! Advection term (iso-sigma) alone
   ! --------------------------------
           ftempx(k,i,j) = ftempx(k,i,j)*dy_u(i,j)
           fsalx(k,i,j)  = fsalx(k,i,j)*dy_u(i,j)
#else
   ! Advection/diffusion term (iso-sigma)
   ! ----------------------------------
#if defined key_siggen || defined key_gencoord
           ftempx(k,i,j) = ftempx(k,i,j)-hzex(k,i,j)*kxb*deltat/dx_u(i,j)*halfdt
           fsalx(k,i,j) = fsalx(k,i,j)-hzex(k,i,j)*kxb*deltas/dx_u(i,j)*halfdt            
#else
           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            
#endif
#endif
#endif

#ifdef key_substance
           DO iv=nv_grav+1,nv_adv
#if defined key_ts_diffh_geo || defined key_ts_diffh_rho
   ! Advection term (iso-sigma) alone
   ! --------------------------------
             fx(iv,k,i,j) = fx(iv,k,i,j)*dy_u(i,j)
#else
   ! Advection/diffusion term (iso-sigma)
   ! ----------------------------------
#if defined key_siggen || defined key_gencoord
             fx(iv,k,i,j) = fx(iv,k,i,j) - hzex(k,i,j)*kxb*delta(iv)/dx_u(i,j)*halfdt
#else
             fx(iv,k,i,j) = fx(iv,k,i,j) - hex(i,j)*kxb*delta(iv)/dx_u(i,j)*halfdt
#endif
#endif
#ifdef key_sedim_MUSTANG
             IF (k == 1) fx(iv,k,i,j)=fx(iv,k,i,j)*corflux(iv,i,j)
#endif

           END DO
#endif

         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

OMPMPI barrier
OMPMPI master
#if ! defined key_offline
   CALL_MPI ex_i_rlg(-1,0,kmax,liminm1,limax,ljmin,ljmax,fsalx(:,liminm1:limax,ljmin:ljmax))
   CALL_MPI ex_i_rlg(-1,0,kmax,liminm1,limax,ljmin,ljmax,ftempx(:,liminm1:limax,ljmin:ljmax))
#endif
#ifdef key_substance
   CALL_MPI ex_i_rsh(-1,0,nv_adv*kmax,liminm1,limax,ljmin,ljmax,fx(:,:,liminm1:limax,ljmin:ljmax))
#endif
OMPMPI end master
OMPMPI barrier
#if ! defined key_offline
OMPMPI flush(fsalx,ftempx)
#endif
#ifdef key_substance
OMPMPI flush(fx)
#endif

#else

   USE parameters, ONLY : rsh,kmax,liminm1,limaxp2,ljminm1,ljmaxp2
   REAL(KIND=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: sal,temp
#ifdef key_substance
   REAL(KIND=rsh),DIMENSION(kmax,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: c
#endif

#endif

   PRINT_DBG*, 'exit TSSUB_ADVDIFX_FIFTHORDER.F90 routine'

  END SUBROUTINE tssub_advdifx_fifthorder
