  SUBROUTINE tssub_advx_macho(sal,temp,sal2,temp2 &
#ifdef key_substance
                                      ,c,c2       &
#endif
                                           )

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE tssub_advx_macho  ***
   !&E
   !&E ** Purpose : Intermediate estimate of the tracers
   !&E
   !&E ** Description : 
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays : htx,hx,uz,sal,temp,c(cv_wat),dx
   !&E
   !&E ** Modified variables : temp2,sal2,c2
   !&E
   !&E ** Reference : Leonard B. P., Lock A. P. and M. K. MacVean (1996) :
   !&E                Conservative Explicit Unrestricted-Time-Step
   !&E                Multidimensional COnstancy-Preserving Advections Schemes.
   !&E                Monthly Weather Review, 124, 2588:2605
   !&E
   !&E ** History :
   !&E       !  2010-05 (V. Garnier) Original code from tssubadvdifxultimatequickestmacho.F90
   !&E       !  2016-04 (B.Thouvenin) applied to all substances (dissolved and particulate)
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d,    ONLY : h0fond,dx,htx,ig,id,hx,halfdt
   USE comvarp3d,    ONLY : uz
#ifdef key_substance
   USE comsubstance, ONLY : nv_adv,nv_grav,nv_tot
#endif
   USE_MPI toolmpi,  ONLY : ex_i_rsh,ex_j_rsh

   IMPLICIT NONE

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in)  :: sal,temp
   REAL(KIND=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(inout) :: sal2,temp2
#ifdef key_substance
   REAL(KIND=rsh),DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in)  :: c
   REAL(KIND=rsh),DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(inout) :: c2
#endif

   !! * Local declarations
   INTEGER                   :: i,j,k,iv
   REAL(KIND=rlg), DIMENSION(kmax) :: uadv

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

   PRINT_DBG*, 'enter TSSUB_ADVX_MACHO.F90 routine'

   ! estimate of intermediate tracers (from upwind scheme)
   ! -----------------------------------------------------

   ! at open boundaries (north and south)

!$OMP SINGLE

   IF_MPI (limin==imin) THEN
#if ! defined key_offline
     temp2(:,MAX0(imin,imin+ishift_obc_west-1):imin+1+ishift_obc_west,ljmin:ljmax) = &
                        temp(:,MAX0(imin,imin+ishift_obc_west-1):imin+1+ishift_obc_west,ljmin:ljmax)
     sal2(:,MAX0(imin,imin+ishift_obc_west-1):imin+1+ishift_obc_west,ljmin:ljmax) =  &
                        sal(:,MAX0(imin,imin+ishift_obc_west-1):imin+1+ishift_obc_west,ljmin:ljmax)
#endif
#ifdef key_substance
     c2(nv_grav+1:nv_adv,:,MAX0(imin,imin+ishift_obc_west-1):imin+1+ishift_obc_west,ljmin:ljmax) = &
                        c(nv_grav+1:nv_adv,:,MAX0(imin,imin+ishift_obc_west-1):imin+1+ishift_obc_west,ljmin:ljmax)
#endif
   ENDIF_MPI
   IF_MPI (limax==imax) THEN
#if ! defined key_offline
     temp2(:,imax-ishift_obc_east:MIN0(imax,imax-ishift_obc_east+2),ljmin:ljmax) = &
                        temp(:,imax-ishift_obc_east:MIN0(imax,imax-ishift_obc_east+2),ljmin:ljmax)
     sal2(:,imax-ishift_obc_east:MIN0(imax,imax-ishift_obc_east+2),ljmin:ljmax) =  &
                        sal(:,imax-ishift_obc_east:MIN0(imax,imax-ishift_obc_east+2),ljmin:ljmax)
#endif
#ifdef key_substance
     c2(nv_grav+1:nv_adv,:,imax-ishift_obc_east:MIN0(imax,imax-ishift_obc_east+2),ljmin:ljmax) = &
                        c(nv_grav+1:nv_adv,:,imax-ishift_obc_east:MIN0(imax,imax-ishift_obc_east+2),ljmin:ljmax)
#endif
   ENDIF_MPI
   IF_MPI (ljmin==jmin) THEN
#if ! defined key_offline
     temp2(:,limin:limax,MAX0(jmin,jmin+jshift_obc_south-1):jmin+1+jshift_obc_south) = &
                        temp(:,limin:limax,MAX0(jmin,jmin+jshift_obc_south-1):jmin+1+jshift_obc_south)
     sal2(:,limin:limax,MAX0(jmin,jmin+jshift_obc_south-1):jmin+1+jshift_obc_south) = &
                        sal(:,limin:limax,MAX0(jmin,jmin+jshift_obc_south-1):jmin+1+jshift_obc_south)
#endif
#ifdef key_substance
     c2(nv_grav+1:nv_adv,:,limin:limax,MAX0(jmin,jmin+jshift_obc_south-1):jmin+1+jshift_obc_south) = &
                        c(nv_grav+1:nv_adv,:,limin:limax,MAX0(jmin,jmin+jshift_obc_south-1):jmin+1+jshift_obc_south)
#endif
   ENDIF_MPI
   IF_MPI (ljmax==jmax) THEN
#if ! defined key_offline
     temp2(:,limin:limax,jmax-jshift_obc_north:MIN0(jmax,jmax-jshift_obc_north+2)) = &
                        temp(:,limin:limax,jmax-jshift_obc_north:MIN0(jmax,jmax-jshift_obc_north+2))
     sal2(:,limin:limax,jmax-jshift_obc_north:MIN0(jmax,jmax-jshift_obc_north+2)) = &
                        sal(:,limin:limax,jmax-jshift_obc_north:MIN0(jmax,jmax-jshift_obc_north+2))
#endif
#ifdef key_substance
     c2(nv_grav+1:nv_adv,:,limin:limax,jmax-jshift_obc_north:MIN0(jmax,jmax-jshift_obc_north+2)) = &
                        c(nv_grav+1:nv_adv,:,limin:limax,jmax-jshift_obc_north:MIN0(jmax,jmax-jshift_obc_north+2))
#endif
   ENDIF_MPI

!$OMP END SINGLE

   ! inside the domain

!$OMP DO SCHEDULE(RUNTIME),PRIVATE(i,j,k,uadv)

   DO j=ljmin,ljmax
#ifdef key_MPI_2D
     DO i=limin,limax
     IF (i >= ig(j)+1 .AND. i <=id(j)-1) THEN
#else
     DO i=ig(j)+1,id(j)-1
#endif

!       DO k=1,kmax

         ! estimate of the velocities
         IF (htx(i  ,j) <= 0.0_rsh .OR. hx(i  ,j) == -valmanq-h0fond .OR. &
             htx(i-1,j) <= 0.0_rsh .OR. hx(i-1,j) == -valmanq-h0fond) THEN
           uadv(:)=0.0_rlg
         ELSE
           uadv(:)=0.5_rsh*(uz(:,i-1,j)+uz(:,i,j))
#ifdef key_test
           DO k=1,kmax
             IF (uz(k,i,j)*uz(k,i-1,j) > 0.0_rsh) THEN
               IF (uz(k,i,j) > 0.0_rsh) THEN
                 uadv(k)=uz(k,i-1,j)
               ELSE
                 uadv(k)=uz(k,i,j)
               END IF
             ELSE
               uadv(k)=0.0_rsh
             END IF
           END DO
#endif
         END IF

         ! estimate of the intermediate tracers
#if ! defined key_offline
         temp2(:,i,j) = temp(:,i,j) - (halfdt/dx(i,j))*uadv(:)*0.5_rsh*(temp(:,i+1,j)-temp(:,i-1,j))
         sal2(:,i,j)  = sal(:,i,j)  - (halfdt/dx(i,j))*uadv(:)*0.5_rsh*(sal(:,i+1,j)-sal(:,i-1,j))
#endif
#ifdef key_substance
         DO iv=nv_grav+1,nv_adv
           c2(iv,:,i,j)  = c(iv,:,i,j)-(halfdt/dx(i,j))*uadv(:)*0.5_rsh*(c(iv,:,i+1,j)-c(iv,:,i-1,j))
         END DO
#endif

!       END DO   ! ends loop on k
#ifdef key_MPI_2D
     END IF
#endif
     END DO   ! ends loop on i
   END DO   ! ends loop on j
!$OMP END DO

OMPMPI barrier
OMPMPI master
#if ! defined key_offline
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#ifdef key_substance
!to avoid temporary copy 1st dimension is not contiguous   CALL_MPI ex_j_rsh(-1,2,(nv_adv)*kmax,limin,limax,ljminm1,ljmaxp2,c2(1:nv_adv,:,limin:limax,ljminm1:ljmaxp2))
   CALL_MPI ex_i_rsh(-1,2,nv_tot*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,c2(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,nv_tot*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,c2(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
OMPMPI end master
OMPMPI barrier
#if ! defined key_offline
OMPMPI flush(sal2,temp2)
#endif
#ifdef key_substance
OMPMPI flush(c2)
#endif

   PRINT_DBG*, 'exit TSSUB_ADVX_MACHO.F90 routine'

  END SUBROUTINE tssub_advx_macho
