  SUBROUTINE tssub2d_eq(                       &
#ifdef key_substance
                        cv_wat,                &
#endif
                               sal,temp,xem,xep)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE tssub2d_eq  ***
   !&E
   !&E ** Purpose : estimate the advection-diffusion and decay of substances
   !&E              for the 2d model 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : step2dxy, step2dyx
   !&E
   !&E ** External calls : tssub_advdifx, tssub_advdify
   !&E
   !&E ** Used ij-arrays : xep h0 xem cv_wat dy_u fx dx_v fy surf sal fsalx
   !&E                     fsaly temp ftempx ftempy phitemp phicon
   !&E
   !&E ** Modified variables : cv_wat, sal, temp
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2004-06-13
   !&E       !  2007-09  (V. Garnier) homogeneneity of dx dy...
   !&E       !  2009-04  (A. Mano, V. Garnier) MPI2D
   !&E       !  2010-06  (P. Lazure) adding of tests on fluxes to update the tracers
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d,   ONLY : ig,id,h0,h0fond,dy_u,dx_v,surf,fwet,fwetp
   USE comfluxdo,   ONLY : fsalx,fsaly,nbfl
   USE cometemp,    ONLY : ftempx,ftempy,phitemp
#ifdef key_substance
   USE comsubstance,ONLY : fx,fy,nv_adv,f90,nbrej,xrej,yrej,nucon,phicon,nv_adv
   USE comvars2d,   ONLY : l_champroche
#endif
#ifdef key_MPI_2D
   USE toolmpi,     ONLY : ex_i_rsh,ex_j_rsh
#endif

   IMPLICIT NONE

   !! * Arguments
   REAL(kind=rsh),DIMENSION(limin:limax,ljmin:ljmax)                    ,INTENT(in)    :: xem
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2)            ,INTENT(in)    :: xep
   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,j,l
   REAL(kind=rsh)                                           :: hmo,hpl

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

!$OMP PARALLEL

#ifdef key_substance
   CALL tssub_advdifx(sal,temp,cv_wat)
   CALL tssub_advdify(sal,temp,cv_wat)
#else
   CALL tssub_advdifx(sal,temp)
   CALL tssub_advdify(sal,temp)
#endif

#ifdef key_substance
!$OMP DO SCHEDULE(RUNTIME) PRIVATE(hpl,hmo)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
       hpl=xep(i,j)+h0(i,j)
       hmo=xem(i,j)+h0(i,j)

       ! attention tel que c est fait ce qui est sur les bancs
       ! decouvrants ne decroit pas pendant la phase de decouvrement
       !------------------------------------------------------------

       IF (hpl >= h0fond .AND. hpl > 0.0_rsh) THEN
         ! .or. hmo >= h0fond) then
         DO l=1,nv_adv
           IF (fx(l,1,i,j)/=0.0_rsh .OR. fx(l,1,i-1,j)/=0.0_rsh .OR. &
               fy(l,1,i,j)/=0.0_rsh .OR. fy(l,1,i,j-1)/=0.0_rsh) THEN
             cv_wat(l,1,i,j)=f90(l)*(hmo*fwetp(i,j)*cv_wat(l,1,i,j) &
                                  -(dy_u(i,j)*fx(l,1,i,j)-dy_u(i-1,j)*fx(l,1,i-1,j)  &
                                   +dx_v(i,j)*fy(l,1,i,j)-dx_v(i,j-1)*fy(l,1,i,j-1)) &
                                  /surf(i,j))/(hpl*fwet(i,j))
           END IF
         END DO
       END IF
     END DO
   END DO
!$OMP END DO
#endif

   IF ( nbfl >= 1 .OR. num_testcase==18 .OR. num_testcase==19 .OR. num_testcase==20 ) THEN
!$OMP DO SCHEDULE(RUNTIME) PRIVATE(hpl,hmo)
     DO j=ljminp2,ljmaxm1
       DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
         hpl=xep(i,j)+h0(i,j)
         hmo=xem(i,j)+h0(i,j)
         !if(hpl >= h0fond .or. hmo >= h0fond) then
         IF(hpl >= h0fond .AND. hpl > 0.0_rsh) THEN
           IF (fsalx(1,i,j)/=0.0_rsh .OR. fsalx(1,i-1,j)/=0.0_rsh .OR. &
               fsaly(1,i,j)/=0.0_rsh .OR. fsaly(1,i,j-1)/=0.0_rsh ) THEN
             sal(1,i,j)=(hmo*fwetp(i,j)*sal(1,i,j)-(dy_u(i,j)*fsalx(1,i,j)-dy_u(i-1,j)*fsalx(1,i-1,j)  &
                                      +dx_v(i,j)*fsaly(1,i,j)-dx_v(i,j-1)*fsaly(1,i,j-1)) &
                                      /surf(i,j))/(hpl*fwet(i,j))
           END IF
           IF (ftempx(1,i,j)/=0.0_rsh .OR. ftempx(1,i-1,j)/=0.0_rsh .OR. &
               ftempy(1,i,j)/=0.0_rsh .OR. ftempy(1,i,j-1)/=0.0_rsh ) THEN
             temp(1,i,j)=(hmo*fwetp(i,j)*temp(1,i,j)-(dy_u(i,j)*ftempx(1,i,j)-dy_u(i-1,j)*ftempx(1,i-1,j)  &
                                        +dx_v(i,j)*ftempy(1,i,j)-dx_v(i,j-1)*ftempy(1,i,j-1)) &
                                        /surf(i,j)+SUM(phitemp(:,i,j)))/(hpl*fwet(i,j))
           END IF
         ENDIF
       ENDDO
     ENDDO
!$OMP END DO
   ENDIF

!$OMP END PARALLEL

#ifdef key_substance
   IF(.NOT.l_champroche) THEN
     DO l=1,nbrej
       i=NINT(xrej(l))
       j=NINT(yrej(l))
       IF_MPI((j.GE.ljmin).AND.(j.LE.ljmax)) THEN
       IF_MPI((i.GE.limin).AND.(i.LE.limax)) THEN
         hpl=xep(i,j)+h0(i,j)
         ! ajout traceur meme sur estran
         IF(hpl > 0.0_rsh) THEN
            cv_wat(nucon(l),1,i,j) = cv_wat(nucon(l),1,i,j)+phicon(nucon(l),1,i,j)/(surf(i,j)*hpl*fwet(i,j))
         ELSE
            PRINT*,'The outflow',l,' is located at a wetting cell : ',i,j
         ENDIF
       ENDIF_MPI
       ENDIF_MPI
     ENDDO
   ENDIF
#endif

  END SUBROUTINE tssub2d_eq
