  SUBROUTINE ts_advh(sal,temp)

   !&E---------------------------------------------------------------------
   !&E ** Purpose : compute horizontal advection of temperature and salinity
   !&E              using TVD scheme (quick and upwind decentred) 
   !&E              with all-or-nothing limiter.
   !&E
   !&E ** Description : For very low water heights ( < hm), 
   !&E                  no horizontal diffusion is applied.
   !&E
   !&E ** Called by : adv
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Used ij-arrays : temp, sal, hex, ig, id, qmax, kx, hm, halfdt, fxeau
   !&E
   !&E ** Modified variables : fadsal, fadtemp
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2006-06  (C. Mazoyer, V. Garnier) original code
   !&E       !  2007-09  (V. Garnier) homogeneneity of dx dy...
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters
   USE comvars2d, ONLY : ig,id,hex,hx,kx,qmax,h0fond,hm,dy_u,hey,hy,ky,jb,jh,dx_v,halfdt
   USE comfluxdo, ONLY : fsalx,fsaly,fxeau,fyeau
   USE cometemp,  ONLY : ftempx,ftempy
   IMPLICIT NONE

   !! * Arguments
   REAL(kind=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(inout) :: sal,temp

   !! * Local declarations
   LOGICAL            :: decentt,decents
   INTEGER            :: i,j,k
   REAL(kind=rsh)     :: unsqmax,uu,yr,qot,dyds,dygs,dydt
   REAL(kind=rsh)     :: dygt,kxb,yas,yat,ybs,ybt,ycs,yct,epsi
   REAL(kind=rsh)     :: kyb,vv

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

   unsqmax=1.0_rsh/qmax
   epsi=1.0e-10

!$OMP DO schedule(runtime)
   DO j=ljminp2,ljmaxm1
     DO i=ig(j),id(j)-1

       !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       !~~~~~~~~~~~~~~~~~~~~ADVECTION EN X~~~~~~~~~~~~~~~~~~~~
       !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       IF(hex(i,j) <= 0.0_rsh .OR. hx(i,j) == -valmanq-h0fond) THEN
         fsalx(:,i,j)=0.0_rlg
         ftempx(:,i,j)=0.0_rlg
       ELSE
         DO k=1,kmax
           kxb=kx
           uu=fxeau(k,i,j)
           IF(uu <= 0.0_rsh) THEN
             ybs=sal(k,i+1,j)
             ycs=sal(k,i,j)
             ybt=temp(k,i+1,j)
             yct=temp(k,i,j)
             IF(hex(i+1,j) > 0.0_rsh .AND. i+1 < id(j)) THEN
               yas=sal(k,i+2,j)
               yat=temp(k,i+2,j)
             ELSE
               yas=ybs
               yat=ybt
             ENDIF
             dyds=ybs-ycs
             dygs=ybs-yas
             dydt=ybt-yct
             dygt=ybt-yat

             decents=(dyds*dygs) >= 0.0_rsh
             IF(.NOT.decents) THEN
               qot=ABS(dyds/MAX(ABS(dygs),epsi))
               IF(qot <= qmax .AND. qot >= unsqmax) THEN
                 decents=.FALSE.
               ELSE
                 decents=.TRUE.
               ENDIF
             ENDIF

             decentt=(dydt*dygt) >= 0.0_rsh
             IF(.NOT.decentt) THEN
               qot=ABS(dydt/MAX(ABS(dygt),epsi))
               IF(qot <= qmax .AND. qot >= unsqmax) THEN
                 decentt=.FALSE.
               ELSE
                 decentt=.TRUE.
               ENDIF
             ENDIF

           ELSE ! option du test sur la vitesse

             ybs=sal(k,i,j)
             ycs=sal(k,i+1,j)
             ybt=temp(k,i,j)
             yct=temp(k,i+1,j)
             IF(hex(i-1,j) > 0.0_rsh) THEN
               yas=sal(k,i-1,j)
               yat=temp(k,i-1,j)
             ELSE
               yas=ybs
               yat=ybt
             ENDIF
             dyds=ycs-ybs
             dygs=ybs-yas
             dydt=yct-ybt
             dygt=ybt-yat

             decents=(dyds*dygs) <= 0.0_rsh
             IF(.NOT.decents) THEN
               qot=ABS(dyds/MAX(ABS(dygs),epsi))
               IF(qot <= qmax .AND. qot >= unsqmax) THEN
                 decents=.FALSE.
               ELSE
                 decents=.TRUE.
               ENDIF
             ENDIF

             decentt=(dydt*dygt) <= 0.0_rsh
             IF(.NOT.decentt) THEN
               qot=ABS(dydt/MAX(ABS(dygt),epsi))
               IF(qot <= qmax .AND. qot >= unsqmax) THEN
                 decentt=.FALSE.
               ELSE
                 decentt=.TRUE.
               ENDIF
             ENDIF
           ENDIF ! fin de test sur la vitesse.

           !  salinite
           !-------------
           IF(decents) THEN
             fsalx(k,i,j)=uu*ybs*dy_u(i,j)
           ELSE
             yr=0.375_rsh*ycs+0.75_rsh*ybs-0.125_rsh*yas
             fsalx(k,i,j)=uu*yr*dy_u(i,j)
           ENDIF

           ! temperature
           !-------------
           IF(decentt) THEN
             ftempx(k,i,j)=uu*ybt*dy_u(i,j)
           ELSE
             yr=0.375_rsh*yct+0.75_rsh*ybt-0.125_rsh*yat
             ftempx(k,i,j)=uu*yr*dy_u(i,j)
           ENDIF

         ENDDO ! boucle en k
       ENDIF ! test sur la hauteur d eau
     ENDDO ! boucle en i
   ENDDO ! boucle en j
!$OMP end do

   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   !~~~~~~~~~~~~~~~~~~~~ADVECTION EN Y~~~~~~~~~~~~~~~~~~~~
   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!$OMP DO schedule(runtime)
   DO j=ljminp1-1,ljmax
     DO i=imin+2,imax-1
       IF(j >= jb(i) .AND. j <= jh(i)-1) THEN
         IF(hey(i,j) <= 0.0_rsh .OR. hy(i,j)==-valmanq-h0fond) THEN
           fsaly(:,i,j)=0.0_rlg
           ftempy(:,i,j)=0.0_rlg
         ELSE
           DO k=1,kmax
             kyb=ky
             vv=fyeau(k,i,j)
             IF(vv <= 0.0_rsh) THEN
               ybs=sal(k,i,j+1)
               ycs=sal(k,i,j)
               ybt=temp(k,i,j+1)
               yct=temp(k,i,j)

               IF(hey(i,j+1) > 0.0_rsh .AND. j+1 < jh(i)) THEN
                 yas=sal(k,i,j+2)
                 yat=temp(k,i,j+2)
               ELSE
                 yas=ybs
                 yat=ybt
               ENDIF
               dyds=ybs-ycs
               dygs=ybs-yas
               dydt=ybt-yct
               dygt=ybt-yat
               decents=(dyds*dygs) >= 0.0_rsh
               IF(.NOT.decents) THEN
                 qot=ABS(dyds/MAX(ABS(dygs),epsi))
                 IF(qot <= qmax .AND. qot >= unsqmax) THEN
                   decents=.FALSE.
                 ELSE
                   decents=.TRUE.
                 ENDIF
               ENDIF

               decentt=(dydt*dygt) >= 0.0_rsh
               IF(.NOT.decentt) THEN
                 qot=ABS(dydt/MAX(ABS(dygt),epsi))
                 IF(qot <= qmax .AND. qot >= unsqmax) THEN
                   decentt=.FALSE.
                 ELSE
                   decentt=.TRUE.
                 ENDIF
               ENDIF

             ELSE ! deuxieme option du test sur les vitesses.

               ybs=sal(k,i,j)
               ycs=sal(k,i,j+1)
               ybt=temp(k,i,j)
               yct=temp(k,i,j+1)
               IF(hey(i,j-1) > 0.0_rsh .AND. j > jb(i)) THEN
                 yas=sal(k,i,j-1)
                 yat=temp(k,i,j-1)
               ELSE
                 yas=ybs
                 yat=ybt
               ENDIF
               dyds=ycs-ybs
               dygs=ybs-yas
               dydt=yct-ybt
               dygt=ybt-yat

               decents=(dyds*dygs) <= 0.0_rsh
               IF(.NOT.decents) THEN
                 qot=ABS(dyds/MAX(ABS(dygs),epsi))
                 IF(qot <= qmax .AND. qot >= unsqmax) THEN
                   decents=.FALSE.
                 ELSE
                   decents=.TRUE.
                 ENDIF
               ENDIF

               decentt=(dydt*dygt) <= 0.0_rsh
               IF(.NOT.decentt) THEN
                 qot=ABS(dydt/MAX(ABS(dygt),epsi))
                 IF(qot <= qmax .AND. qot >= unsqmax) THEN
                   decentt=.FALSE.
                 ELSE
                   decentt=.TRUE.
                 ENDIF
               ENDIF

             ENDIF ! fin du test sur la vitesse

             !  salinite
             !-----------
             IF(decents) THEN
               fsaly(k,i,j)=vv*ybs*dx_v(i,j)
             ELSE
               yr=0.375_rsh*ycs+0.75_rsh*ybs-0.125_rsh*yas
               fsaly(k,i,j)=vv*yr*dx_v(i,j)
             ENDIF

             ! temperature
             !-------------
             IF(decentt) THEN
               ftempy(k,i,j)=vv*ybt*dx_v(i,j)
             ELSE
               yr=0.375_rsh*yct+0.75_rsh*ybt-0.125_rsh*yat
               ftempy(k,i,j)=vv*yr*dx_v(i,j)
             ENDIF

           ENDDO ! fin de la boucle en k
         ENDIF ! fin du test sur la hauteur d eau
       ENDIF ! for selecting right j loop
     ENDDO ! fin de la boucle en i
   ENDDO ! fin de la boucle en j
!$OMP end do

END SUBROUTINE ts_advh
