  SUBROUTINE turb_advz_compact(xe,ect,psi)

#ifdef key_turbadv
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE turb_advz_compact  ***
   !&E
   !&E ** Purpose : compute vertical fluxes of psi and ect
   !&E              using an Upwind Compact and Conservative scheme
   !&E              adapted to an euler time scheme
   !&E
   !&E ** Description : For very low water heights ( < hm), no vertical fluxes are computed
   !&E              because the scheme is degenerated vertically and vertical fluxes are useless.
   !&E              This routine is called in "adv" only if the key
   !&E              key_turbadv_wcompact is activated.
   !&E
   !&E ** Called by : adv
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Used ij-arrays : xe,h0,fzeau,wz,ect,psi
   !&E
   !&E ** Modified variables : fzpsi,fzect
   !&E
   !&E ** Reference : Debreu, 2008 (note : UCC : Upwind Compact and Conservative)
   !&E
   !&E ** History :
   !&E       !  2015-08 (F. Ganthy) Original code: adaptation from tsub_advz_compact
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used

   USE parameters
   USE comvars2d,  ONLY : id,ig,h0,hm,halfdt
   USE comfluxdo,  ONLY : fzeau
   USE comvars3d,  ONLY : dsigu,wz,fzect,fzpsi

   IMPLICIT NONE

   !! * Arguments
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in)      :: xe
   REAL(kind=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: ect,psi

   !! * Local declarations
   INTEGER                                  :: i,j,k,n
   REAL(kind=rsh)                           :: h3,cfl,alpha0,alpha1,alpha2,gamma1,gamma2,mu1,mu2
   REAL(kind=rsh),DIMENSION(kmax)           :: ats,bts,cts,y2dt,y2ds,ecthat,psihat

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

!$OMP DO SCHEDULE(RUNTIME) &
!$OMP& PRIVATE(i,j,k,n,h3,cfl,alpha0,alpha1,alpha2,gamma1,gamma2,mu1,mu2 &
!$OMP&, ats,bts,cts,y2dt,y2ds,ecthat,psihat)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)),MIN0(limax,id(j)-1)

   ! coefficients of Upwind Compact and Conservative scheme
   !-------------------------------------------------------

       h3=h0(i,j)+xe(i,j)
       IF (h3 > hm) THEN
       DO k=1,kmax-1
         !cfl=(wz(k,i,j)+wz(k-1,i,j))*0.5_rsh*REAL(halfdt,rsh)/dsigu(k)
         cfl=wz(k,i,j)*REAL(halfdt,rsh)/dsigu(k)
         IF (cfl > 0.0_rsh) THEN
           alpha1=3.0_rsh*(cfl-3.0_rsh)*(cfl-2.0_rsh)
           alpha2=3.0_rsh*(cfl+1.0_rsh)*(cfl+2.0_rsh)
           alpha0=-6.0_rsh*(cfl-3.0_rsh)*(cfl+2.0_rsh)
           gamma1=38.0_rsh-2.0_rsh*cfl**2
           gamma2=20.0_rsh+(cfl-3.0_rsh)*cfl
           mu1=(cfl+1.0_rsh)*(cfl+2.0_rsh)
           mu2=0.0_rsh
         ELSE
           alpha1=3.0_rsh*(cfl-1.0_rsh)*(cfl-2.0_rsh)
           alpha2=3.0_rsh*(cfl+3.0_rsh)*(cfl+2.0_rsh)
           alpha0=-6.0_rsh*(cfl-2.0_rsh)*(cfl+3.0_rsh)
           gamma1=20.0_rsh+(cfl+3.0_rsh)*cfl
           gamma2=38.0_rsh-2.0_rsh*cfl**2
           mu1=0.0_rsh
           mu2=(cfl-1.0_rsh)*(cfl-2.0_rsh)
         END IF

         ats(k)=alpha1
         bts(k)=alpha0
         cts(k)=alpha2

         IF (k >= 2 .AND. k <= kmax-2) THEN
           y2ds(k)=mu1*ect(k-1,i,j) +gamma1*ect(k,i,j) +gamma2*ect(k+1,i,j) +mu2*ect(k+2,i,j)
           y2dt(k)=mu1*psi(k-1,i,j)+gamma1*psi(k,i,j)+gamma2*psi(k+1,i,j)+mu2*psi(k+2,i,j)
         END IF
         IF (k==1) THEN
           y2ds(k)=(mu1+gamma1-alpha1)*ect(k,i,j) +gamma2*ect(k+1,i,j) +mu2*ect(k+2,i,j)
           y2dt(k)=(mu1+gamma1-alpha1)*psi(k,i,j)+gamma2*psi(k+1,i,j)+mu2*psi(k+2,i,j)
         END IF
         IF (k==kmax-1) THEN
           y2ds(k)= mu1*ect(k-1,i,j)+gamma1*ect(k,i,j) +(mu2+gamma2-alpha2)*ect(k+1,i,j)
           y2dt(k)=mu1*psi(k-1,i,j)+gamma1*psi(k,i,j)+(mu2+gamma2-alpha2)*psi(k+1,i,j)
         END IF
#ifdef key_plustard
! comprendre le schema de ces OBC
         IF (k==1) THEN
           IF (cfl > 0.0_rsh) THEN
             alpha1=0.0_rsh
             alpha2=(cfl+1.0_rsh)*(1.0_rsh+5.0_rsh*cfl)
             alpha0=(cfl-1.0_rsh)*(1.0_rsh-5.0_rsh*cfl)
             y2ds(k)=(cfl**2-1.0_rsh)*ect(k,i,j) +(1.0_rsh-(cfl-12.0_rsh)*cfl)*ect(k+1,i,j) 
             y2dt(k)=(cfl**2-1.0_rsh)*psi(k,i,j)+(1.0_rsh-(cfl-12.0_rsh)*cfl)*psi(k+1,i,j)
           ELSE
             alpha1=0.0_rsh
             alpha2=2.0_rsh*(1.0_rsh+cfl)
             alpha0=4.0_rsh-2.0_rsh*cfl
             y2ds(k)=(1.0_rsh+cfl)*ect(k,i,j) +(5.0_rsh-cfl)*ect(k+1,i,j)
             y2dt(k)=(1.0_rsh+cfl)*psi(k,i,j)+(5.0_rsh-cfl)*psi(k+1,i,j)
           END IF
         END IF
! revoir ce schema comme ci-dessus
         IF (k==kmax-1) THEN
           alpha1=0.0_rsh
           alpha2=0.0_rsh
           alpha0=1.0_rsh
           IF (cfl > 0.0_rsh) THEN
             y2ds(k)=ect(k,i,j)
             y2dt(k)=psi(k,i,j)
           ELSE
             y2ds(k)=ect(k+1,i,j)
             y2dt(k)=psi(k+1,i,j)
           END IF
         END IF
#endif

       END DO

       n=kmax-1
       CALL  inversion_tridiag(ats(1:kmax-1),bts(1:kmax-1),cts(1:kmax-1), ecthat(1:kmax-1),y2ds(1:kmax-1),n)
       CALL  inversion_tridiag(ats(1:kmax-1),bts(1:kmax-1),cts(1:kmax-1),psihat(1:kmax-1),y2dt(1:kmax-1),n)

   ! pour faire simple : decentre amont pres des limites (fond et surface)
   !----------------------------------------------------------------------

       fzpsi(0,i,j)=0.0_rsh
       fzpsi(kmax,i,j)=0.0_rsh
       fzect(0,i,j)=0.0_rsh
       fzect(kmax,i,j)=0.0_rsh

       DO k=1,kmax-1
         fzect(k,i,j) =fzeau(k,i,j)*ecthat(k)
         fzpsi(k,i,j)=fzeau(k,i,j)*psihat(k)
       END DO

       ENDIF     ! fin du test sur la hauteur d eau
     END DO     ! fin de la boucle en i
   END DO     ! fin de la boucle en j
!$OMP end do
#endif

  END SUBROUTINE turb_advz_compact

