   SUBROUTINE turbadv_dify(ect,psi)

#ifdef key_turbadv
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE turbadv_dify  ***
   !&E
   !&E ** Purpose : Estimate of meridional 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 : fyect,fypsi,fy
   !&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,ky,h0fond,hminkxky,dy_u,dx,hty,hey,    &
                            jb,jh,hy,halfdt,l_obc_cycl_y,dy_v
   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,dygp,d0,d1,vvv,courant,vv,  &
                                kyb,yae,yap,ybe,ybp,yce,ycp,thetap,psip,      &
                                deltae,deltap
   REAL(KIND=rlg), DIMENSION(kmax) :: ful,fur
   REAL(KIND=rlg), DIMENSION(kmax) :: ectl,ectr,psil,psir

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

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

!$OMP SINGLE

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

!$OMP END SINGLE

   ! * Inside the domaine

!$OMP DO SCHEDULE(RUNTIME),PRIVATE(i,j,k,ectl,ectr,psil,psir,ful,fur,uu)

   DO j = ljmin,ljmax
#ifdef key_MPI_2D
     DO i = liminp1m1,limax
     IF (i.EQ.ig(j) .OR. i.EQ.liminp1m1) THEN
#else
     i=ig(j)
#endif
       DO k = 1,kmax-1
         IF (htx(i,j).LE.0.0_rsh .OR. hx(i,j).EQ.-valmanq-h0fond) THEN
           ! No adevection from this cell
           ectr(k) = 0.0_rsh
           psir(k) = 0.0_rsh
         ELSE
           fur(k)=fxsw(k,i,j)   ! h*uz*dt/2
           ! Upwind scheme
           IF (fur(k).LE.0.0_rlg) THEN
             ectr(k) = ect(k,i+1,j)
             psir(k) = psi(k,i+1,j)
           ELSE
             ectr(k) = ect(k,i,j)
             psir(k) = psi(k,i,j)
           ENDIF ! End test on fur<=0.0
         ENDIF ! End test on htx of hx
       ENDDO ! End loop on k

#ifdef key_MPI_2D
     ELSE IF(i.GE.ig(j)+1 .AND. i.LE.id(j)-1 .AND. i.GT.liminp1m1) THEN
#else
     DO i = ig(j)+1,id(j)-1
#endif

       DO k = 1,kmax-1
         ! Left of the stencil
         ectl(k) = ectr(k)
         psil(k) = psir(k)
         ful(k) = fur(k)
         ! Right of the stencil
         IF (htx(i,j).LE.0.0_rsh .OR. hx(i,j).EQ.valmanq-h0fond) THEN
           ! No adevection from this cell
           ectr(k) = 0.0_rsh
           psir(k) = 0.0_rsh
           fur(k)  = 0.0_rlg
         ELSE
           fur(k) = fxsw(k,i,j) ! h*uz*dt/2
           ! Upwind scheme
           IF (fur(k).LE.0.0_rlg) THEN
             ectr(k) = ect(k,i+1,j)
             psir(k) = psi(k,i+1,j)
           ELSE
             ectr(k) = ect(k,i,j)
             psir(k) = psi(k,i,j)
           ENDIF ! End test on fur<=0.0
         ENDIF ! End test on htx or hx

         ! Estimate of intermediate tke and psi (from upwind scheme)
         IF (ful(k).GT.0.0_rlg .AND. fur(k).GT.0.0_rlg) THEN
           uu = ful(k)/(halfdt*hex(i,j))
         ELSEIF (ful(k).LT.0.0_rlg .AND. fur(k).LT.0.0_rlg) THEN
           uu = fur(k)/(halfdt*hex(i,j))
         ELSE
           uu = 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/dx(i,j))*uu*(ectr(k)-ectl(k))
         psi2(k,i,j) = psi(k,i,j) - (halfdt/dx(i,j))*uu*(psir(k)-psil(k))

       ENDDO ! End loop on k
#ifdef key_MPI_2D
     END IF ! End test on ig, id
#endif
     END DO   ! ends loop on i
   END DO   ! ends loop on j
!$OMP END DO

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_j_rsh(-1,2,kmax-1,liminm1,limaxp2,ljminm1,ljmaxp2,ect2(:,liminm1:limaxp2,ljminm1:ljmaxp2)
   CALL_MPI ex_j_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,dyde,dyge,dydp,dyge,kyb,yae,yap,ybe,ybp, &
!$OMP& yce,ycp,deltae,deltap,d0,d1,vv,vvv,courant,thetap,psip)

   DO j = ljmin,ljmax
     DO i = liminp2,limaxm1
       IF (j.GE.jb(i) .AND. j.LE.jh(i)-1) THEN
         DO k = 1,kmax-1
           IF (hty(i,j).LE.0.0_rsh .OR. hy(i,j).EQ.-valmanq-h0fond) THEN
             fyect(k,i,j) = 0.0_rlg
             fypsi(k,i,j) = 0.0_rlg
           ELSE
             kyb = 0.0_rsh!ky
             IF (hty(i,j).LT.hminkxky) kyb = 0.0_rsh ! No diffusion for small water depth
             vv = fysw(k,i,j)
             vvv = vv/(halfdt*hey(i,j))
             courant = ABS(vvv)*halfdt/dy_v(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 (vv.LE.0.0_rsh) THEN ! Southward advection
               ybe = ect2(k,i,j+1)
               yce = ect2(k,i,j)
               ybp = psi2(k,i,j+1)
               ycp = psi2(k,i,j)
               IF (hty(i,j+1).GT.0.0_rsh .AND. (j+1.LT.jh(i) .OR. l_obc_cycl_y)) THEN
                 yae = ect2(k,i,j+2)
                 yap = psi2(k,i,j+2)
               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 ! Northward advection
               ybe = ect2(k,i,j)
               yce = ect2(k,i,j+1)
               ybp = ect2(k,i,j)
               ycp = ect2(k,i,j+1)
               IF (hty(i,j-1).GT.0.0_rsh .AND. (j.GT.jb(i) .OR. l_obc_cycl_y)) THEN
                 yae = ect2(k,i,j-1)
                 yap = ect2(k,i,j-1)
               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
             fyect(k,i,j) = vv*(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
             fypsi(k,i,j) = vv*(ybp+psip*dydp) 

             ! Advection/diffusion term
             !--------------------------
             fyect(k,i,j) = fyect(k,i,j)-hey(i,j)*kyb*deltae/dy_v(i,j)*halfdt
             fypsi(k,i,j) = fypsi(k,i,j)-hey(i,j)*kyb*deltap/dy_v(i,j)*halfdt

           ENDIF ! End test on water depth
         ENDDO ! End loop on k
       ENDIF ! End test for selecting right j loop
     ENDDO ! End loop on i
   ENDDO  ! End loop on j

! $OMP END DO

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_j_rlg(-1,0,kmax-1,limin,limax,ljminm1,ljmax,fyect(:,limin:limax,ljminm1:ljmax))
   CALL_MPI ex_j_rlg(-1,0,kmax-1,limin,limax,ljminm1,ljmax,fypsi(:,limin:limax,ljminm1:ljmax))
OMPMPI end master
OMPMPI barrier

OMPMPI flush(fyect,fypsi)
#endif

   END SUBROUTINE turbadv_dify

