#ifdef key_substance

  SUBROUTINE parsub_eq_hm(xe3,cv_wat,ll,niter,num_ws)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE tssub_eq_hm  ***
   !&E
   !&E ** Purpose : compute the transport of tracers with fall velocity, 
   !&E              for meshes with low water heights (< hm).
   !&E
   !&E ** Description : WARNING : The resulting vertical profiles are homogenous.
   !&E
   !&E ** Called by : adv
   !&E
   !&E ** External calls :                  :
   !&E
   !&E ** Used ij-arrays : h0,xe3
   !&E                     flx_s2w,flx_w2s,flx_w2s_sum
   !&E                     cv_wat
   !&E
   !&E ** Modified variables : cv_wat
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2010-03  (P. Douillet) substances with settling velocities
   !&E       !                         issued from tssub_eq_hm
   !&E       !  2010-06  (V.M. Faure, N. Bousserez, M. Baklouti) coupling with eco3m
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2016-04  (B.Thouvenin) only vertical transport due to settling velocoty
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d,    ONLY : ig,id,h0,h0fond,hm,dt,halfdt,fwet
   USE comvars3d,    ONLY : dsigu
#if defined key_siggen || defined key_gencoord
   USE toolgeom,     ONLY : f_dzu
#endif
   USE comsubstance, ONLY : nvp,nv_adv,nv_ws_max,iv_ws_max,l_subs2D
#if defined key_sedim_MUSTANG
   USE comsubstance, ONLY : isand1,isand2
   USE comMUSTANG,   ONLY : flx_s2w,flx_w2s,flx_w2s_sum
#endif

   IMPLICIT NONE

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

   !! * Local declarations
   INTEGER                                         :: i,iv,j,l1,k
   REAL(KIND=rsh)                                  :: cnouv,d3c3,d3
   REAL(KIND=rsh)                                  :: cm
   REAL(KIND=rlg)                                  :: dtm
   REAL(KIND=rsh),DIMENSION(kmax)                  :: dz

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

   dtm=0.5_rlg*dt/REAL(niter,rlg)

#ifdef key_MPI_2D
!$OMP DO schedule(runtime) private(i,l1,j,k,cnouv,d3c3,d3 &
!$OMP& ,cm,dz)
#else
!$OMP DO schedule(dynamic,3)
#endif
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
       d3=h0(i,j)+xe3(i,j)
       IF(d3 > h0fond.AND.d3 <= hm) THEN

#if defined key_siggen || defined key_gencoord
         DO k=1,kmax
           dz(k)=f_dzu(h0(i,j),xe3(i,j),k,i,j)   ! thickness of the layer [m]
         END DO
#else
         dz(:)=dsigu(:)   !  no unit
#endif

   !-----------------------------------------------
   ! integration des flux de traceurs particulaires
   !-----------------------------------------------

         DO l1=1,nv_ws_max(num_ws) ! number of variables in the group of settling velocity
           iv=iv_ws_max(num_ws,l1)  ! initialize index of the sustance
                                   ! according to the number of settling velocity group
                                   ! and the number of variables in this group
           IF( l_subs2D(iv)) THEN     ! for variables "sable" or "sand"

#if defined key_siggen || defined key_gencoord
              cm=cv_wat(iv,1,i,j)*d3
#else
              cm=cv_wat(iv,1,i,j)
#endif
           ELSE
              cm=SUM(cv_wat(iv,:,i,j)*dz(:))
           ENDIF

#if defined key_siggen || defined key_gencoord
           d3c3=cm
#else
           d3c3=d3*cm
#endif
#if defined key_sedim_MUSTANG
           d3c3=d3c3+dtm*flx_s2w(iv,i,j)
           cnouv=d3c3/(d3+dtm*flx_w2s(iv,i,j))
           IF(cnouv < 1.e-10) cnouv=0.0_rsh
           flx_w2s_sum(iv,i,j)=flx_w2s_sum(iv,i,j)+flx_w2s(iv,i,j)*cnouv*dtm
#else
           cnouv=d3c3/d3
           IF(cnouv < 1.e-10) cnouv=0.0_rsh
#endif
           IF( l_subs2D(iv)) THEN     ! for variables "sable" or "sand"
              cv_wat(iv,1,i,j)=cnouv
              cv_wat(iv,2:kmax,i,j)=0.0_rsh
           ELSE
              cv_wat(iv,:,i,j)=cnouv
           ENDIF

         ENDDO

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

!$OMP end DO

  END SUBROUTINE parsub_eq_hm

#else

  SUBROUTINE parsub_eq_hm_empty
    WRITE(*,*) 'parsub_eq_hm_empty routine :'
    WRITE(*,*) 'You should not have seen this print!  error?'
    return
  END SUBROUTINE parsub_eq_hm_empty

#endif
