  SUBROUTINE tssub_eq_hm(xe2,xe3,sal,temp        &
#ifdef key_substance
                                         ,cv_wat &
#endif
                                                )

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE tssub_eq_hm  ***
   !&E
   !&E ** Purpose : compute the transport of temperature, salinity and
   !&E              other tracers without 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,xe2,sal,fsalx,fsaly,phisal,dy_u,dx_v,
   !&E                     surf,flx_s2w,flx_w2s,temp,ftempx,ftempy,phitemp,
   !&E                     cv_wat,fx,fy,phicon,dcdt
   !&E
   !&E ** Modified variables : sal,temp,cv_wat
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2004-12-27
   !&E       !  2008-05-07 (M. Sourisseau, V. Garnier) correction vertical interpolation of dcdt
   !&E       !  2009-04    (A. Mano, V. Garnier) MPI2D
   !&E       !  2010-06    (V.M. Faure, N. Bousserez, M. Baklouti) coupling with eco3m
   !&E       !  2010-06  (F. Cayocca) coupling with the module of sedimentology
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2016-04  (B.Thouvenin) applied to all substances (dissolved and particulate, include sand)
   !&E       !  2016-12  (B.Thouvenin) version sedim_MUSTANG : only mixsed (mudsand not operationnal)
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d,    ONLY : ig,id,h0,h0fond,hm,dy_u,dx_v,surf,dt,halfdt,fwet,fwetp
   USE comvars3d,    ONLY : dsigu
#if defined key_siggen || defined key_gencoord
   USE toolgeom,     ONLY : f_dzu
#endif
   USE comfluxdo,    ONLY : fsalx,fsaly
   USE cometemp,     ONLY : ftempx,ftempy,phitemp,phisal
#ifdef key_substance
   USE comsubstance, ONLY : fx,fy,nvp,nv_adv,phicon,f90,typdiss,l_subs2D
#if defined key_sedim_MUSTANG
!!#if defined key_sedim_mudsand
!!   USE comsubstance, ONLY : flx_s2w,flx_w2s
!!#else
   USE comMUSTANG, ONLY : flx_s2w
!!#endif
#endif
#if defined key_contaminant || defined key_biolo || defined key_eco3m
   USE comsubstance, ONLY : dcdt
#endif
#endif

   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
   REAL(KIND=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),       INTENT(inout) :: sal,temp
#ifdef key_substance
   REAL(KIND=rsh),DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(inout) :: cv_wat
#endif

   !! * Local declarations
   INTEGER                                         :: i,iv,j,k
   REAL(KIND=rsh)                                  :: cnouv,d3c3,d3,d2
   REAL(KIND=rsh)                                  :: fxmg,fxmd,fymb,fymh,cm,phi,dcdthm
   REAL(KIND=rlg)                                  :: dtm
   REAL(KIND=rsh),DIMENSION(kmax)                  :: dz


   !!----------------------------------------------------------------------
   !! * Executable part
   
   dtm=0.5_rlg*dt

#ifdef key_MPI_2D
!$OMP DO schedule(runtime) private(i,j,k,cnouv,d3c3,d3,d2 &
!$OMP& ,fxmg,fxmd,fymb,fymh,cm,phi,dcdthm,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)
       d2=h0(i,j)+xe2(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),xe2(i,j),k,i,j)
         END DO
#else
         dz(:)=dsigu(:)   !  no unit
#endif

!--------------------------------------
! integration des flux de salinite
!--------------------------------------

         !if(i==101 .and.j==203) then
!        !  write(*,*)'tssubeqhm',d2,d3
         !endif
         cm=SUM(sal(:,i,j)*dz(:))
         fxmg=SUM(fsalx(:,i-1,j)*dsigu(:))
         fxmd=SUM(fsalx(:,i,j)*dsigu(:))
         fymb=SUM(fsaly(:,i,j-1)*dsigu(:))
         fymh=SUM(fsaly(:,i,j)*dsigu(:))
         phi =SUM(phisal(:,i,j))
#if defined key_siggen || defined key_gencoord
         d3c3=(cm*fwetp(i,j)-(fxmd*dy_u(i,j)-fxmg*dy_u(i-1,j)+fymh*dx_v(i,j)-fymb*dx_v(i,j-1))/surf(i,j)+phi)/fwet(i,j)
#else
         d3c3=(d2*cm*fwetp(i,j)-(fxmd*dy_u(i,j)-fxmg*dy_u(i-1,j)+fymh*dx_v(i,j)-fymb*dx_v(i,j-1))/surf(i,j)+phi)/fwet(i,j)
#endif
#ifdef key_sedim_MUSTANG
         d3c3=d3c3+dtm*flx_s2w(0,i,j)
! modificatin des gestions des variables dissoutes dans le sediment
! flux w2s traites explicitement et deja rajouter dans flx_s2w
#ifdef key_sedim_mudsand
         cnouv=d3c3/(d3+dtm*flx_w2s(0,i,j))
#else
         cnouv=d3c3/d3
#endif
#else
         cnouv=d3c3/d3
#endif
         IF(cnouv < 1.e-10) cnouv=0.0_rsh
         sal(:,i,j)=cnouv

!--------------------------------------
! integration des flux de temperature
!--------------------------------------

         cm=SUM(temp(:,i,j)*dz(:))
         fxmg=SUM(ftempx(:,i-1,j)*dsigu(:))
         fxmd=SUM(ftempx(:,i,j)*dsigu(:))
         fymb=SUM(ftempy(:,i,j-1)*dsigu(:))
         fymh=SUM(ftempy(:,i,j)*dsigu(:))
         phi =SUM(phitemp(:,i,j))
#if defined key_siggen || defined key_gencoord
         d3c3=(cm*fwetp(i,j)-(fxmd*dy_u(i,j)-fxmg*dy_u(i-1,j)+fymh*dx_v(i,j)-fymb*dx_v(i,j-1))/surf(i,j)+phi)/fwet(i,j)
#else
         d3c3=(d2*cm*fwetp(i,j)-(fxmd*dy_u(i,j)-fxmg*dy_u(i-1,j)+fymh*dx_v(i,j)-fymb*dx_v(i,j-1))/surf(i,j)+phi)/fwet(i,j)
#endif
#ifdef key_sedim_MUSTANG
         d3c3=d3c3+dtm*flx_s2w(-1,i,j)
! modificatin des gestions des variables dissoutes dans le sediment
! flux w2s traites explicitement et deja rajouter dans flx_s2w
!!#ifdef key_sedim_mudsand
!!         cnouv=d3c3/(d3+dtm*flx_w2s(-1,i,j))
!!#else
         cnouv=d3c3/d3
!!#endif
#else
         cnouv=d3c3/d3
#endif
         IF(cnouv < 1.e-10) cnouv=0.0_rsh
         temp(:,i,j)=cnouv

#ifdef key_substance
!-----------------------------------------
! integration des flux de traceurs autres
!----------------------------------------- 

         DO iv=1,nv_adv
           IF(l_subs2D(iv)) THEN
#if defined key_siggen || defined key_gencoord
             cm=cv_wat(iv,1,i,j)*d2
#else
             cm=cv_wat(iv,1,i,j)
#endif
             fxmg=fx(iv,1,i-1,j)
             fxmd=fx(iv,1,i,j)
             fymb=fy(iv,1,i,j-1)
             fymh=fy(iv,1,i,j)
           ELSE
             cm=SUM(cv_wat(iv,:,i,j)*dz(:))
             fxmg=SUM(fx(iv,:,i-1,j)*dsigu(:))
             fxmd=SUM(fx(iv,:,i,j)*dsigu(:))
             fymb=SUM(fy(iv,:,i,j-1)*dsigu(:))
             fymh=SUM(fy(iv,:,i,j)*dsigu(:))
           ENDIF
           phi =SUM(phicon(iv,:,i,j))
#if defined key_biolo || defined key_contaminant || defined key_eco3m
#ifdef key_biolo
           !dcdthm=SUM(dcdt(iv,:,i,j)*dz(:))
           ! correction BT juillet 2015 : dcdt calcule dans biolosinksource et biolodyncellwat
           !         dans la couche de fond mais correspondant a une moyenne sur la profondeur
           dcdthm=dcdt(iv,1,i,j)*d2
#else
#if defined key_siggen || defined key_gencoord
           dcdthm=SUM(dcdt(iv,:,i,j)*dz(:))
#else
           dcdthm=d2*SUM(dcdt(iv,:,i,j)*dsigu(:))
#endif
#endif
#endif
#if defined key_siggen || defined key_gencoord
           d3c3=(cm*fwetp(i,j)-(fxmd*dy_u(i,j)-fxmg*dy_u(i-1,j)+fymh*dx_v(i,j)-fymb*dx_v(i,j-1)-phi)/surf(i,j))/fwet(i,j)
#else
           d3c3=(d2*cm*fwetp(i,j)-(fxmd*dy_u(i,j)-fxmg*dy_u(i-1,j)+fymh*dx_v(i,j)-fymb*dx_v(i,j-1)-phi)/surf(i,j))/fwet(i,j)
#endif
#if defined key_biolo || defined key_contaminant || defined key_eco3m
           d3c3=d3c3+dtm*dcdthm
#endif
#ifdef key_sedim_MUSTANG
! pas de flux a l interface eau sediment pour les particulaires (traites avec eq vert.)
           d3c3=d3c3+dtm*flx_s2w(iv,i,j)*typdiss(iv)
! modification des gestions des variables dissoutes dans le sediment
! flux w2s traites explicitement et deja rajouter dans flx_s2w pour mixsed
!!#ifdef key_sedim_mudsand
!!           cnouv=d3c3/(d3+dtm*flx_w2s(iv,i,j)*typdiss(iv))
!!#else
           cnouv=d3c3/d3
!!#endif
#else
           cnouv=d3c3/d3
#endif
           IF(cnouv < 1.e-10) cnouv=0.0_rsh
           IF(l_subs2D(iv)) THEN     ! for variables "sable" or "sand"
              cv_wat(iv,1,i,j)=f90(iv)*(cnouv)
              cv_wat(iv,2:kmax,i,j)=0.0_rsh
           ELSE
              cv_wat(iv,:,i,j)=f90(iv)*(cnouv)
           ENDIF
         ENDDO
#endif

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

!$OMP end DO

  END SUBROUTINE tssub_eq_hm
