  SUBROUTINE tssub_advz_macho(xe,sal,temp,sal2,temp2 &

                                      ,c,c2       &

                                           )

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE tssub_advz_macho  ***
!&E
!&E ** Purpose : Estimate of intermediate tracers
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : hey,hy,wz,sal,temp,c(cv_wat),dsigu
!&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
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d,    ONLY : h0,hm,h0fond,ig,id,halfdt
   USE comvars3d,    ONLY : wz,dsigu

   USE comsubstance, ONLY : nv_adv,nvpi,nv_tot

   USE toolmpi,  ONLY : ex_i_rsh,ex_j_rsh

   IMPLICIT NONE

!! * Arguments
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),         INTENT(in)  :: xe
   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=rsh)            :: h3
   REAL(KIND=rlg), DIMENSION(kmax) :: wadv

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

   !PRINT*, 'enter TSSUB_ADVZ_MACHO.F90 routine'

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

!cval toseelater HUM que se passe t til pour les partie non initialiesee lors des echanges ?
!cval que fait on au fond

! 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(nvpi+1:nv_adv,:,MAX0(imin,imin+ishift_obc_west-1):imin+1+ishift_obc_west,ljmin:ljmax) = &
                        c(nvpi+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(nvpi+1:nv_adv,:,imax-ishift_obc_east:MIN0(imax,imax-ishift_obc_east+2),ljmin:ljmax) = &
                        c(nvpi+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(nvpi+1:nv_adv,:,limin:limax,MAX0(jmin,jmin+jshift_obc_south-1):jmin+1+jshift_obc_south) = &
                        c(nvpi+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(nvpi+1:nv_adv,:,limin:limax,jmax-jshift_obc_north:MIN0(jmax,jmax-jshift_obc_north+2)) = &
                        c(nvpi+1:nv_adv,:,limin:limax,jmax-jshift_obc_north:MIN0(jmax,jmax-jshift_obc_north+2))

   ENDIF

! at the bottom
   temp2(1,limin:limax,ljmin:ljmax) = temp(1,limin:limax,ljmin:ljmax)
   sal2(1,limin:limax,ljmin:ljmax) = sal(1,limin:limax,ljmin:ljmax)

! at the sea surface
   temp2(kmax,limin:limax,ljmin:ljmax) = temp(kmax,limin:limax,ljmin:ljmax)
   sal2(kmax,limin:limax,ljmin:ljmax) = sal(kmax,limin:limax,ljmin:ljmax)

!$OMP END SINGLE

! inside the domain

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

   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)

       h3=h0(i,j)+xe(i,j)
       IF(h3 > hm) THEN

!         k=1
!         wadv(k)=0.5_rsh*(wz(k,i,j))

         DO k=2,kmax-1

! estimate of the advective velocity
! a revoir car pas sur que mpyenne soit correct en siggen
!wadv(k)=0.5_rsh*(wz(k-1,i,j)+wz(k,i,j))
           IF (wz(k,i,j)*wz(k-1,i,j) > 0.0_rsh) THEN
             IF (wz(k,i,j) > 0.0_rsh) THEN
               wadv(k)=wz(k-1,i,j)
             ELSE
               wadv(k)=wz(k,i,j)
             END IF
           ELSE
             wadv(k)=0.0_rsh
           END IF

! estimate of intermediate tracers

           temp2(k,i,j) = temp(k,i,j) - (halfdt/dsigu(k))*wadv(k)*0.5_rsh*(temp(k+1,i,j)-temp(k-1,i,j))
           sal2(k,i,j)  = sal(k,i,j)  - (halfdt/dsigu(k))*wadv(k)*0.5_rsh*(sal(k+1,i,j)-sal(k-1,i,j))


           c2(nvpi+1:nv_adv,k,i,j)  = c(nvpi+1:nv_adv,k,i,j)  - (halfdt/dsigu(k))*wadv(k) &
                                         *0.5_rsh*(c(nvpi+1:nv_adv,k+1,i,j)-c(nvpi+1:nv_adv,k-1,i,j))


         END DO   ! ends loop on k

       ELSE  ! to initialize correctly (for nbouc_macho/=0 and coherent values of tracers above small depths)
         temp2(:,i,j) = temp(:,i,j)
         sal2(:,i,j) = sal(:,i,j)

         c2(nvpi+1:nv_adv,:,i,j)  = c(nvpi+1:nv_adv,:,i,j)

       END IF   ! ends test for selecting right j loop

     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-nvpi)*kmax,liminm1,limaxp2,ljmin,ljmax,c2(nvpi+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_ADVZ_MACHO.F90 routine'

  END SUBROUTINE tssub_advz_macho
