   SUBROUTINE turbadv_difx(ect,psi)

#ifdef key_turbadv
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE turbadv_difx  ***
   !&E
   !&E ** Purpose : Estimate of zonal fluxes of turbulent kinetic energy and
   !&E              turbulence dissipation (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,fxsw,fysw,ect,psi,hzex,dx_u,dy,
   !&E                     ect2,psi2,c2,fxect,fxpsi,fx 
   !&E
   !&E ** Modified variables : fxect,fxpsi,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           !  2015-08 (F. Ganthy) original code adapted from tssubadvdifxultimatequickestmacho.F90
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * 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 comvars3d
   USE comvarp3d, ONLY : ect2,psi2

   USE_MPI toolmpi,      ONLY : ex_i_rsh,ex_i_rlg

   IMPLICIT NONE

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(kmax,limin:limax,ljmin:ljmax),INTENT(inout) :: ect,psi

   !! * Local declarations
   INTEGER                   :: i,j,k
   REAL(KIND=rsh), PARAMETER :: dsmall = 1.0d-20
   REAL(KIND=rsh)            :: uu,dyde,dyge,dydp,d0,d1,uuu,courant,vv,  &
                                dygp,kxb,yae,yap,ybe,ybp,yce,ycp,thetap, &
                                psip,deltae,deltap
   REAL(KIND=rlg), DIMENSION(kmax) :: fvb,fvt
   REAL(KIND=rsh), DIMENSION(kmax) :: ectt,ectb,psit,psib

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

   ! estimate of intermediate TKE and PSI (from upwind scheme)
   !----------------------------------------------------------

!$OMP SINGLE

   ! * At open boundaries
   IF_MPI (limin.EQ.imin) THEN
     ect2(:,imin:imin+1,ljmin:ljmax) = ect(:,imin:imin+1,ljmin:ljmax)
     psi2(:,imin:imin+1,ljmin:ljmax) = psi(:,imin:imin+1,ljmin:ljmax)
   ENDIF_MPI
   IF_MPI (limax.EQ.imax) THEN
     ect2(:,imax,ljmin:ljmax) = ect(:,imax,ljmin:ljmax)
     psi2(:,imax,ljmin:ljmax) = psi(:,imax,ljmin:ljmax)
   ENDIF_MPI

!$OMP END SINGLE

   ! * Inside the domaine

!$OMP DO SCHEDULE(RUNTIME),PRIVATE(i,j,k,ectt,ectb,psit,psib,fvt,fvb,vv)

   DO i = liminp2,limaxm1
     DO j = ljminp1m1,ljmax
       DO k = 1,kmax-1
#ifdef key_MPI_2D
         IF (j.EQ.jb(i) .OR. j.EQ.ljminp1m1) THEN
#else
         IF (j.EQ.jb(i) .OR. j.EQ.ljminp1-1) THEN
#endif
           IF (hty(i,j).LE.0.0_rsh .OR. hy(i,j).EQ.-valmanq-h0fond) THEN
             ! No fluxe advection from this cell
             ectt(k) = 0.0_rsh
             psit(k) = 0.0_rsh
             fvt(k) = 0.0_rlg
           ELSE
             fvt(k) = fysw(k,i,j) ! h*vz*dt/2
             ! upwind scheme
             IF (fvt(k).LE.0.0_rlg) THEN
               ectt(k) = ect(k,i,j+1)
               psit(k) = psi(k,i,j+1)
             ELSE
               ectt(k) = ect(k,i,j)
               psit(k) = psi(k,i,j)
             ENDIF ! End test on fvt<=0.0
           ENDIF ! End test on hty or hy

#ifdef key_MPI_2D
         ELSEIF (j.GE.jb(i)+1 .AND. j.LE.jh(i)-1 .AND. j.GT.ljminp1m1) THEN
#else
         ELSEIF (j.GE.jb(i)+1 .AND. j.LE.jh(i)-1 .AND. j.GT.ljminp1-1) THEN
#endif
           ! Bottom of the stencil
           ectb(k) = ectt(k)
           psib(k) = psit(k)
           fvb(k)  = fvt(k)

           ! Right of the stencil
           IF (hty(i,j).LE.0.0_rsh .OR. hy(i,j).EQ.-valmanq-h0fond) THEN
             ! No fluxe advection from this cell
             ectt(k) = 0.0_rsh
             psit(k) = 0.0_rsh
             fvt(k) = 0.0_rlg
           ELSE
             fvt(k) = fysw(k,i,j) ! h*vz*dt/2
             ! upwind scheme
             IF (fvt(k).LE.0.0_rlg) THEN
               ectt(k) = ect(k,i,j+1)
               psit(k) = psi(k,i,j+1)
             ELSE
               ectt(k) = ect(k,i,j)
               psit(k) = psi(k,i,j)
             ENDIF ! End test on fvt<=0.0
           ENDIF ! End test on hty or hy

           ! Estimate of intermediate tke and psi (from upwind scheme)
           IF (fvb(k).GT.0.0_rlg .AND. fvt(k).GT.0.0_rlg) THEN
             vv = fvb(k)/(halfdt*hey(i,j-1))
           ELSEIF (fvb(k).LT.0.0_rlg .AND. fvt(k).LT.0.0_rlg) THEN
             vv = fvt(k)/(halfdt*hey(i,j))
           ELSE
             vv = 0.0_rsh ! No time evolution along the boundaries
           ENDIF

           ! Apply euler (time) - upwind (x) scheme
           ect2(k,i,j) = ect(k,i,j) - (halfdt/dy(i,j))*vv*(ectt(k)-ectb(k))
           psi2(k,i,j) = psi(k,i,j) - (halfdt/dy(i,j))*vv*(psit(k)-psib(k))

         ENDIF ! End test on j=jb or j=ljminp1m1

       ENDDO ! End loop on k
     ENDDO ! End loop on j
   ENDDO ! End loop on i

!$OMP END DO

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-1,2,kmax-1,liminm1,limaxp2,ljminm1,ljmaxp2,ect2(:,liminm1:limaxp2,ljminm1:ljmaxp2)
   CALL_MPI ex_i_rsh(-1,2,kmax-1,liminm1,limaxp2,ljminm1,ljmaxp2,psi2(:,liminm1:limaxp2,ljminm1:ljmaxp2)
OMPMPI end master
OMPMPI barrier
OMPMPI flush(ect2,psi2)

   ! Estimate fluxes (from ultimate quickest scheme)
   !------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME),PRIVATE(i,j,k,kxb,dyde,dydp,duge,dygp, &
!$OMP& yae,yap,ybe,ybp,yce,ycp,deltae,deltap,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-1
         IF(htx(i,j).LE.0.0_rsh .OR. hx(i,j).EQ.-valmanq-h0fond) THEN
           fxect(k,i,j) = 0.0_rlg
           fxpsi(k,i,j) = 0.0_rlg
         ELSE
           kxb=0.0_rsh!kx
           IF (htx(i,j).LT.hminkxky) kxb = 0.0_rsh ! No diffusion for small water depth
           uu  = fxsw(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.LE.0.0_rsh) THEN ! Westward advection
             ybe = ect2(k,i+1,j)
             yce = ect2(k,i,j)
             ybp = psi2(k,i+1,j)
             ycp = psi2(k,i,j)
             IF (htx(i+1,j).GT.0.0_rsh .AND. (i+1.LT.id(j) .OR. l_obc_cycl_x)) THEN
               yae = ect2(k,i+2,j)
               yap = psi2(k,i+2,j)
             ELSE
               yae = ybe
               yap = ybp
             ENDIF
             dyde = yce-ybe
             dyge = ybe-yae
             dydp = ycp-ybp
             dygp = ybp-yap
             deltae = ybe-yce
             deltap = ybp-ycp
           ELSE ! Eastward advection
             ybe = ect2(k,i,j)
             yce = ect2(k,i+1,j)
             ybp = psi2(k,i,j)
             ycp = psi2(k,i+1,j)
             IF (htx(i-1,j).GT.0.0_rsh .AND. (i.GT.ig(j) .OR. l_obc_cycl_x)) THEN
               yae = ect2(k,i-1,j)
               yap = psi2(k,i-1,j)
             ELSE
               yae = ybe
               yap = ybp
             ENDIF
             dyde = yce-ybe
             dyge = ybe-yae
             dydp = ycp-ybp
             dygp = ybp-yap
             deltae = ybe-yce
             deltap = ybp-ycp
           ENDIF ! End test on velocity sign


           ! Advection fluxe of TKE
           !------------------------
           IF (ABS(dyde).LE.dsmall .OR. ABS(dyge).LE.dsmall) THEN
             ! No horizontal gradient
             thetap = 0.0_rsh
             psip   = 0.0_rsh
           ELSE
             ! Ultimate limiter
             thetap = (dyge+dsmall)/(dyde+dsmall)
             psip   = d0+d1*thetap
             psip=MAX(0.0_rsh,MIN( MIN(1.0_rsh,psip),((1.0_rsh-courant)/(courant+dsmall))*thetap) )
           ENDIF
           fxect(k,i,j) = uu*(ybe+psip*dyde)

           ! Advection fluxe of PSI
           !------------------------
           IF (ABS(dydp).LE.dsmall .OR. ABS(dygp).LE.dsmall) THEN
             ! No horizontal gradient
             thetap = 0.0_rsh
             psip   = 0.0_rsh
           ELSE
             ! Ultimate limiter
             thetap = (dygp+dsmall)/(dydp+dsmall)
             psip   = d0+d1*thetap
             psip=MAX(0.0_rsh,MIN( MIN(1.0_rsh,psip),((1.0_rsh-courant)/(courant+dsmall))*thetap) )
           ENDIF
           fxpsi(k,i,j) = uu*(ybp+psip*dydp)


           ! Advection/diffusion term
           !--------------------------
           fxect(k,i,j) = fxect(k,i,j)-hex(i,j)*kxb*deltae/dx_u(i,j)*halfdt
           fxpsi(k,i,j) = fxpsi(k,i,j)-hex(i,j)*kxb*deltap/dx_u(i,j)*halfdt

         ENDIF ! End test on water depth
       ENDDO ! End loop on k
     ENDDO ! End loop on i
   ENDDO ! End loop on j
!$OMP END DO

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rlg(-1,0,kmax-1,liminm1,limax,ljmin,ljmax,fxect(:,liminm1:limax,ljmin:ljmax))
   CALL_MPI ex_i_rlg(-1,0,kmax-1,liminm1,limax,ljmin,ljmax,fxpsi(:,liminm1:limax,ljmin:ljmax))
OMPMPI end master
OMPMPI barrier

OMPMPI flush(fxect,fxpsi)
#endif

   END SUBROUTINE turbadv_difx
