  SUBROUTINE turbadv_wcompact(xe2,xe3,ect,psi)

#ifdef key_turbadv
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE turbadv_wcompact  ***
   !&E
   !&E ** Purpose : other possibility of subroutine body for the computation 
   !&E              of transport of ect, psi. 
   !&E
   !&E ** Description : An Euler-Compact scheme is used for depths higher than hm.
   !&E                  For water heights lower than hm, a special computation
   !&E                  is performed in the subroutine tssub_eq_hm.
   !&E                  This subroutine body replaces the default body
   !&E                  when the key key_turbadv_wcompact is actived.
   !&E
   !&E ** Called by : stepu,stepv
   !&E
   !&E ** External calls : turbadv_difx,turbadv_dify,turbadvzcompact
   !&E
   !&E ** Used ij-arrays : ect,psi,cv_wat,h0,xe2,xe3,kz,flx_w2s,
   !&E                     fxect,fyect,fdifect,phiect,
   !&E                     fxpsi,fypsi,fdifpsi,phipsi,
   !&E                     fx,fy,phicon,dy_u,dx_v,surf
   !&E                     fzpsi,fzect,fz
   !&E  
   !&E ** Modified variables : ect,psi,cv_wat
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2015-08  (F. Ganthy) adaptation from subroutine tssub_eq_wcompact
   !&E  
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d,   ONLY : ig,id,h0,hm,dy_u,dx_v,surf,l_tssubadv_over_smalldepth,halfdt,fwet,fwetp
   USE comvars3d

   IMPLICIT NONE

   !! * Arguments
   REAL(kind=rsh),DIMENSION(limin:limax,ljmin:ljmax),                    INTENT(in)    :: xe2
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),            INTENT(in)    :: xe3
!cval xe2 not ready for tsdiffh
   REAL(kind=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),       INTENT(inout) :: ect,psi

   !! * Local declarations
   INTEGER                                   :: n,i,j,k,l
   REAL(kind=rsh)                            :: dtm,tfleuve,d3,d2,unsd3
   REAL(kind=rlg)                            :: xk,yk
   REAL(kind=rlg),DIMENSION(kmax)            :: a,b,c,yect,u0,v0,p0,q0,ypsi

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

   PRINT_DBG*, 'enter TURB_EQ_WCOMPACT.H routine'

!$OMP parallel default(none)                               &
!$OMP& PRIVATE(d2,d3,unsd3,i,j,k,l,xk,yk,                  &
!$OMP&         dtm,p0,q0,u0,v0,a,b,c,                      &
!$OMP&                                 yect,ypsi)         &
!$OMP& shared(ect,psi,xe3,phiect,phipsi,hm,ig,id,h0      &
!$OMP& ,kz,dsigw,dsigu,dy_u                                &
!$OMP& ,limin,limax,ljmin,ljmax,ljminp2,ljmaxm1            &
!$OMP& ,liminm1,limaxp2,ljminm1,ljmaxp2,fwet,fwetp           &
!$OMP& ,fxect,dx_v,fyect,fzect,surf,fxpsi,fypsi,fzpsi   &
!$OMP& ,l_tssubadv_over_smalldepth,xe2,halfdt)

   dtm=halfdt

   ! calcul des flux horizontaux et verticaux (advectifs et diffusifs)
   !------------------------------------------------------------------

   ! advection (+ horizontal diffusion if not key_ts_diffh_geo or key_ts_diffh_rho)
   ! ------------------------------------------------------------------------------

#if ! defined key_1DV
   CALL turbadv_difx(ect,psi)
   CALL turbadv_dify(ect,psi)
   CALL turb_advz_compact(xe3,ect,psi)
#endif
!$OMP DO schedule(runtime)

   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)

       d3=h0(i,j)+xe3(i,j)
       d2=h0(i,j)+xe2(i,j)
       IF(d3 > hm) THEN
         unsd3=1.0_rsh/d3

! matrice tridiagonale du modele vertical
!----------------------------------------

!-------------------------
         k=1
!-------------------------
         c(k)=-dtm*kz(k,i,j)*unsd3/dsigw(k)/dsigu(k)
         b(k)=d3*fwet(i,j)+dtm*kz(k,i,j)*unsd3/dsigw(k)/dsigu(k)
         a(k)=0.0_rlg

         yect(k)=d2*fwetp(i,j)*ect(k,i,j)-(dy_u(i,j)*fxect(k,i,j)-dy_u(i-1,j)*fxect(k,i-1,j) &
                +dx_v(i,j)*fyect(k,i,j)-dx_v(i,j-1)*fyect(k,i,j-1))/surf(i,j)     

         yect(k)=yect(k)-(fzect(k,i,j)-fzect(k-1,i,j))/dsigu(k)

         ypsi(k)=d2*fwetp(i,j)*psi(k,i,j)-(dy_u(i,j)*fxpsi(k,i,j)-dy_u(i-1,j)*fxpsi(k,i-1,j) &
                 +dx_v(i,j)*fypsi(k,i,j)-dx_v(i,j-1)*fypsi(k,i,j-1))/surf(i,j) 
                 
         ypsi(k)=ypsi(k)-(fzpsi(k,i,j)-fzpsi(k-1,i,j))/dsigu(k)


         q0(k)=-c(k)/b(k)        ! factorisation LU

         u0(k)=yect(k)/b(k)      ! ectinite
         v0(k)=ypsi(k)/b(k)     ! psierature

!-------------------------
         DO k=2,kmax-1
!-------------------------
           c(k)=-dtm*kz(k,i,j)*unsd3/dsigw(k)/dsigu(k)
           a(k)=-dtm*kz(k-1,i,j)*unsd3/dsigw(k-1)/dsigu(k)
           b(k)=d3*fwet(i,j)+dtm*(kz(k,i,j)*unsd3/dsigw(k)+kz(k-1,i,j)*unsd3/dsigw(k-1))/dsigu(k)

           yect(k)=d2*fwetp(i,j)*ect(k,i,j)-(dy_u(i,j)*fxect(k,i,j)-dy_u(i-1,j)*fxect(k,i-1,j) &
                  +dx_v(i,j)*fyect(k,i,j)-dx_v(i,j-1)*fyect(k,i,j-1))/surf(i,j)

           yect(k)=yect(k)-(fzect(k,i,j)-fzect(k-1,i,j))/dsigu(k)

           ypsi(k)=d2*fwetp(i,j)*psi(k,i,j)-(dy_u(i,j)*fxpsi(k,i,j)-dy_u(i-1,j)*fxpsi(k,i-1,j) &
                   +dx_v(i,j)*fypsi(k,i,j)-dx_v(i,j-1)*fypsi(k,i,j-1))/surf(i,j)

           ypsi(k)=ypsi(k)-(fzpsi(k,i,j)-fzpsi(k-1,i,j))/dsigu(k)

           p0(k)=a(k)*q0(k-1)+b(k)                ! factorisation LU
           q0(k)=-c(k)/p0(k)

           u0(k)=(yect(k)-a(k)*u0(k-1))/p0(k)     ! ectinite
           v0(k)=(ypsi(k)-a(k)*v0(k-1))/p0(k)    ! psierature

!-------------------------
         ENDDO
!-------------------------

!-------------------------
         k=kmax
!-------------------------
         a(k)=-dtm*kz(k-1,i,j)*unsd3/dsigw(k-1)/dsigu(k)
         b(k)=d3*fwet(i,j)+dtm*kz(k-1,i,j)*unsd3/dsigw(k-1)/dsigu(k)
         c(k)=0.0_rlg

         yect(k)=d2*fwetp(i,j)*ect(k,i,j)-(dy_u(i,j)*fxect(k,i,j)-dy_u(i-1,j)*fxect(k,i-1,j) &
                +dx_v(i,j)*fyect(k,i,j)-dx_v(i,j-1)*fyect(k,i,j-1))/surf(i,j)

         yect(k)=yect(k)-(fzect(k,i,j)-fzect(k-1,i,j))/dsigu(k)

         ypsi(k)=d2*fwetp(i,j)*psi(k,i,j)-(dy_u(i,j)*fxpsi(k,i,j)-dy_u(i-1,j)*fxpsi(k,i-1,j) &
                 +dx_v(i,j)*fypsi(k,i,j)-dx_v(i,j-1)*fypsi(k,i,j-1))/surf(i,j)

         ypsi(k)=ypsi(k)-(fzpsi(k,i,j)-fzpsi(k-1,i,j))/dsigu(k)


         p0(k)=a(k)*q0(k-1)+b(k)                ! factorisation LU

         u0(k)=(yect(k)-a(k)*u0(k-1))/p0(k)     ! ectinite
         v0(k)=(ypsi(k)-a(k)*v0(k-1))/p0(k)    ! psierature

!-------------------------------------------------
!   resolution des systemes lineaires tridiagonaux
!-------------------------------------------------

         xk=u0(kmax)
         ect(kmax,i,j)=xk
         yk=v0(kmax)
         psi(kmax,i,j)=yk
         DO k=kmax-1,1,-1
           xk=q0(k)*xk+u0(k)
           ect(k,i,j)=xk
           yk=q0(k)*yk+v0(k)
           psi(k,i,j)=yk
         ENDDO

       ENDIF ! fin test sur la hauteur
     ENDDO ! fin boucle en i
   ENDDO ! fin boucle en j

!$OMP END DO
!$OMP END PARALLEL

   PRINT_DBG*, 'exit TURB_ADV_WCOMPACT.F90 routine'
#endif

  END SUBROUTINE turbadv_wcompact
