  SUBROUTINE tssub_advy_macho(sal,temp,sal2,temp2 &

                                      ,c,c2       &

                                           )

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE tssub_advy_macho  ***
!&E
!&E ** Purpose : Estimate of intermediate tracers
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : hty,hy,vz,sal,temp,c(cv_wat),dy
!&E
!&E ** Modified variables : sal2,temp2,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, include sand)
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d,    ONLY : h0fond,dy,hty,jb,jh,hy,halfdt
   USE comvarp3d,    ONLY : vz

   USE comsubstance, ONLY : nv_adv,nv_grav,nv_tot

   USE 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

   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


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

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

   !PRINT*, 'enter TSSUB_ADVY_MACHO.F90 routine'

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

! at open boundaries (east and west)

!$OMP SINGLE

   IF (limin==imin) THEN

     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)


     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
   IF (limax==imax) THEN

     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)


     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
   IF (ljmin==jmin) THEN

     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)


     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
   IF (ljmax==jmax) THEN

     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))


     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

!$OMP END SINGLE

! inside the domain

!$OMP DO SCHEDULE(RUNTIME),PRIVATE(i,j,k,vadv)
   DO i=liminp2,limaxm1
     DO j=ljminp1,ljmax
!       DO k=1,kmax

         IF (j >= jb(i)+1 .AND. j <= jh(i)-1) THEN

! estimate of the advective velocity
           IF (hty(i,j  ) <= 0.0_rsh .OR. hy(i,j  )==-valmanq-h0fond .OR. &
               hty(i,j-1) <= 0.0_rsh .OR. hy(i,j-1)==-valmanq-h0fond) THEN
             vadv(:)=0.0_rlg
           ELSE
             vadv(:)=0.5_rsh*(vz(:,i,j-1)+vz(:,i,j))
           END IF

! estimate of intermediate tracers

           temp2(:,i,j) = temp(:,i,j) - (halfdt/dy(i,j))*vadv(:)*0.5_rsh*(temp(:,i,j+1)-temp(:,i,j-1))
           sal2(:,i,j)  = sal(:,i,j)  - (halfdt/dy(i,j))*vadv(:)*0.5_rsh*(sal(:,i,j+1)-sal(:,i,j-1))


           DO iv=nv_grav+1,nv_adv
             c2(iv,:,i,j)  = c(iv,:,i,j)  - (halfdt/dy(i,j))*vadv(:)*0.5_rsh*(c(iv,:,i,j+1)-c(iv,:,i,j-1))
           END DO


         END IF   ! ends test for selecting right j loop

!       END DO   ! ends loop on k
     END DO   ! ends loop on i
   END DO   ! ends loop on j
!$OMP END DO

!$OMP barrier
!$OMP master

   CALL ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))


! to avoid temporary copy   CALL ex_i_rsh(-1,2,(nv_adv)*kmax,liminm1,limaxp2,ljmin,ljmax,c2(1:nv_adv,:,liminm1:limaxp2,ljmin:ljmax))
   CALL ex_i_rsh(-1,2,nv_tot*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,c2(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL ex_j_rsh(-1,2,nv_tot*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,c2(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))

!$OMP end master
!$OMP barrier

!$OMP flush(sal2,temp2)


!$OMP flush(c2)


   !PRINT*, 'exit TSSUB_ADVY_MACHO.F90 routine'

  END SUBROUTINE tssub_advy_macho
