  SUBROUTINE sub_advz_compact(xe,c)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE sub_advz_compact  ***
   !&E
   !&E ** Purpose : compute vertical fluxes of dissolved substances AND PARTICULATE (WITHOUT SELLTING VELOCITY)
   !&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_tssub_wcompact is activated (see tssub_eq_wcompact.h).
   !&E
   !&E ** Called by : adv
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Used ij-arrays : xe,h0,fzeau,wz,c(cv_wat)
   !&E
   !&E ** Modified variables : fz
   !&E
   !&E ** Reference : Debreu, 2008 (note : UCC : Upwind Compact and Conservative)
   !&E
   !&E ** History :
   !&E       !  2010-05 (L. Debreu, F. Dumas, V. Garnier) Original code
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2012-08  (V. Garnier) Use of w(k,i,j) to estimate cfl (Courant–Friedrichs–Lewy) variable
   !&E                                rewritting of the coefficients
   !&E                                review of the OBC for the estimate of interpolated concentrations
   !&E       !  2016-04  (B.Thouvenin) applied to all substances (dissolved and particulate, include sand)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
#ifdef key_substance
   USE parameters
   USE comvars2d,   ONLY : id,ig,h0,hm,halfdt
   USE comvars3d,   ONLY : dsigu,wz
   USE comsubstance,ONLY : nv_adv,nv_grav,fz
   USE comfluxdo,   ONLY : fzeau
   IMPLICIT NONE

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in)              :: xe
   REAL(KIND=rsh),DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: c

   !! * Local declarations
   INTEGER                                  :: i,j,k,n,iv
   REAL(KIND=rsh)                           :: h3,cfl,alpha0,alpha1,alpha2,gamma1,gamma2,mu1,mu2
   REAL(KIND=rsh),DIMENSION(kmax)           :: atc,btc,ctc,chat
   REAL(KIND=rsh),DIMENSION(kmax,nb_var)    :: y2dc

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

!$OMP DO SCHEDULE(RUNTIME) &
!$OMP& PRIVATE(i,j,k,n,iv,h3,cfl,alpha0,alpha1,alpha2,gamma1,gamma2,mu1,mu2 &
!$OMP&, atc,btc,ctc,y2dc,chat)
   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

         atc(k)=alpha1
         btc(k)=alpha0
         ctc(k)=alpha2

         DO iv=nv_grav+1,nv_adv
           IF (k >= 2 .and. k <= kmax-2) THEN 
             y2dc(k,iv)=mu1*c(iv,k-1,i,j) +gamma1*c(iv,k,i,j) +gamma2*c(iv,k+1,i,j) +mu2*c(iv,k+2,i,j)
           END IF
           IF (k==1) THEN
             y2dc(k,iv)=(mu1+gamma1-alpha1)*c(iv,k,i,j) +gamma2*c(iv,k+1,i,j) +mu2*c(iv,k+2,i,j)
           END IF
           IF (k==kmax-1) THEN
             y2dc(k,iv)= mu1*c(iv,k-1,i,j)+gamma1*c(iv,k,i,j) +(mu2+gamma2-alpha2)*c(iv,k+1,i,j)
           END IF
         END DO

       END DO

       n=kmax-1


       DO iv=nv_grav+1,nv_adv
         CALL inversion_tridiag(atc(1:kmax-1),btc(1:kmax-1),ctc(1:kmax-1),chat(1:kmax-1),y2dc(1:kmax-1,iv),n)

         ! pour faire simple : decentre amont pres des limites (fond et surface)
         !----------------------------------------------------------------------
         fz(iv,0:kmax,i,j)=0.0_rsh
         DO k=1,kmax-1
           fz(iv,k,i,j) =fzeau(k,i,j)*chat(k)
         END DO

       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

#else

   USE parameters, ONLY : rsh,kmax,liminm1,ljminm1,liminm2,limaxp2,ljminm2,ljmaxp2
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in)             :: xe
   REAL(kind=rsh),DIMENSION(kmax,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: c

#endif

  END SUBROUTINE sub_advz_compact
