  SUBROUTINE parsub_advz_compact(xe3,setvel,c,ll,niter,num_ws)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE parsub_advz_compact  ***
   !&E
   !&E ** Purpose : compute vertical fluxes of substances with settling velocities
   !&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_parsub_wcompact is activated (see parsub_eq_wcompact.h).
   !&E              niter : total number of iteration or sub time step
   !&E              ll : numerous of the iteration
   !&E              num_ws : index of the settling velocity
   !&E                       (index of variable for 4D settling velocities)
   !&E
   !&E ** Important remark about the settling velocity with the use of key_subs_part_eqsubdt cpp key
   !&E              To save computation time, we assume the substances have the same settling velocity
   !&E              in a group (since the group and the corresponding sub time step depend on
   !&E              the diameter and the density of the particule).
   !&E              This assumption is not necessary valid for the biology. If so, modify slightly
   !&E              the maximum settling velocity in variable.dat file, to create a new group.
   !&E
   !&E ** Called by : adv
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Used ij-arrays : xe,h0,fzeau_ws,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       !  2010-06 (V.M. Faure, N. Bousserez, M. Baklouti) coupling with eco3m
   !&E       !  2011-02 (B. Thouvenin) correction for key_subs_part_eqsubdt key + sub time step
   !&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) only vertical transport due to settling velocoty
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
#if defined key_substance && defined key_parsub_wcompact
   USE parameters
   USE comvars2d,   ONLY : id,ig,h0,hm,halfdt
   USE comvars3d,   ONLY : dsigu
   USE comsubstance,ONLY : iv_ws_max,nv_ws_max,fz
   USE comfluxdo,   ONLY : fzeau_ws
#if defined key_siggen || defined key_gencoord
   USE toolgeom,  ONLY   : f_hzw
#endif
#if defined key_sedim_MUSTANG || defined key_biolo
   USE comsubstance,ONLY : ws3
#endif

   IMPLICIT NONE

   !! * Arguments
   INTEGER,                                          INTENT(in) :: ll,niter,num_ws
   REAL(KIND=rsh),                                   INTENT(in) :: setvel
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),            INTENT(in) :: xe3
   REAL(KIND=rsh),DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: c

   !! * Local declarations
   INTEGER                                  :: i,j,k,n,iv,l1
   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
   REAL(KIND=rsh),DIMENSION(1:kmax)         :: ws
#if defined key_siggen || defined key_gencoord
   REAL(KIND=rsh)                        :: hz3_w
#endif

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

!$OMP DO SCHEDULE(RUNTIME) &
!$OMP& PRIVATE(i,j,k,n,l1,iv,h3,cfl,alpha0,alpha1,alpha2,gamma1,gamma2,mu1,mu2,   &
#if defined key_siggen || defined key_gencoord
!$OMP&         hz3_w,                                                     &
#endif
!$OMP&         atc,btc,ctc,y2dc,chat)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)

   ! get the settling velocity
   ! -------------------------

#if defined key_sedim_MUSTANG || defined key_biolo
       iv=iv_ws_max(num_ws,1)
       ws(:)=-ws3(:,iv,i,j)
!!!#elif defined key_subs_part_eqsubdt || defined key_eco3m
!!!       ws(:)=-setvel
#else
   ! constant settling velocity
       ws(:)=-setvel
#endif


   ! Time interpolation of the sea surface height
   ! --------------------------------------------

       h3=h0(i,j)+xe3(i,j)
       IF (h3 > hm) THEN
       DO k=1,kmax-1

   ! estimate flux of particulate substance
   ! --------------------------------------
#if defined key_siggen || defined key_gencoord
         hz3_w=f_hzw(h0(i,j),xe3(i,j),k,i,j)
         fzeau_ws(k,i,j)=(0.5_rsh*(ws(k)+ws(k+1)))*REAL(halfdt,rsh)/REAL(niter,rsh)
#else
         fzeau_ws(k,i,j)=(0.5_rsh*(ws(k)+ws(k+1)))*REAL(halfdt,rsh)/REAL(niter,rsh)
#endif

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

#if defined key_siggen || defined key_gencoord
         cfl=(ws(k)/hz3_w)*REAL(halfdt,rsh)/dsigu(k)/REAL(niter,rsh)
#else
         cfl=(ws(k)/h3)*REAL(halfdt,rsh)/dsigu(k)/REAL(niter,rsh)
#endif
         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 l1=1,nv_ws_max(num_ws)
           iv=iv_ws_max(num_ws,l1)
           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 l1=1,nv_ws_max(num_ws)
         iv=iv_ws_max(num_ws,l1)
         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_ws(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,limin,limax,ljmin,ljmax,liminm1,liminm2,limaxp2,ljminm1,ljminm2,ljmaxp2
   INTEGER,                                          INTENT(in) :: ll,niter,num_ws
   REAL(KIND=rsh),                                   INTENT(in) :: setvel
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in) :: xe3
   REAL(KIND=rsh),DIMENSION(kmax,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: c


#endif

  END SUBROUTINE parsub_advz_compact
