   !!======================================================================

  SUBROUTINE ts_obc3d_apply_agrif

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE ts_obc3d_apply_agrif  ***
   !&E
   !&E ** Purpose : apply open boundary conditions for temperature
   !&E              and salinity field for AGRIF
   !&E
   !&E ** Description : apply external value when water mass enters the domain
   !&E
   !&E ** Called by : stepu, stepv
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !                       Original code
   !&E       !  2007-05 (V. Garnier) Mise en forme - introduction de obc_* variables
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
#ifdef key_agrif

   USE parameters
   USE comvars2d, ONLY : c_agrif_interp
   USE comvarp2d
   USE agrifcommonid
   USE obccombine, ONLY : l_obc_char,l_obc_south,l_obc_north,l_obc_east,l_obc_west

   USE_MPI toolmpi, ONLY : ex_i_rsh,ex_j_rsh

   !! * Arguments

   !! * Local declarations
   External :: interp3dtraceur
   INTEGER :: i,j

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

   PRINT_DBG*, 'enter TS_OBC3D_APPLY_AGRIF routine'

   Agrif_UseSpecialValue = .TRUE.
   Agrif_SpecialValue = valmanq

   c_agrif_interp='traceur_bc'
   CALL Agrif_Bc_Variable(traceur_id,procname=interp3dtraceur)
   Agrif_UseSpecialValue = .FALSE.

   IF (Agrif_GrandMother()) THEN
     IF ((l_obc_west .OR. l_obc_east)) THEN

       DO j=ljminp2,ljmaxm1

   ! Western boundary
         i=ig(j)
         IF_MPI (limin==i-1) THEN
           IF (hx(i,j) /= -valmanq-h0fond) THEN
              sal(:,i,j) = sal(:,i+1,j)
              temp(:,i,j) = temp(:,i+1,j)
           ELSE
              sal(:,i,j) = valmanq
              temp(:,i,j) = valmanq
           END IF
         ENDIF_MPI

   ! Eastern boundary
         i=id(j)
         IF_MPI (limax==i) THEN
           IF (hx(i-1,j) /= -valmanq-h0fond) THEN
              sal(:,i,j) = sal(:,i-1,j)
              temp(:,i,j) = temp(:,i-1,j)
           ELSE
              sal(:,i,j) = valmanq
              temp(:,i,j) = valmanq
           END IF
         ENDIF_MPI

       END DO

     END IF  ! ends test on .not.l_obc_cycl_x

     IF ((l_obc_north .OR. l_obc_south)) THEN

       DO j=ljmin,ljmax
       DO i=liminp2,limaxm1

   ! Southern boundary
         IF (j == jb(i))THEN
           IF (hy(i,j) /= -valmanq-h0fond) THEN
             sal(:,i,j) = sal(:,i,j+1)
             temp(:,i,j) = temp(:,i,j+1)
           ELSE
             sal(:,i,j) = valmanq
             temp(:,i,j) = valmanq
           END IF

   ! Northern boundary
         ELSE IF (j == jh(i))THEN
           IF (hy(i,j-1) /= -valmanq-h0fond) THEN
             sal(:,i,j) = sal(:,i,j-1)
             temp(:,i,j) = temp(:,i,j-1)
           ELSE
             sal(:,i,j) = valmanq
             temp(:,i,j) = valmanq
           END IF
         END IF

       END DO
       END DO

     END IF  ! ends test on .not.l_obc_cycl_y

   END IF

! keep these exchanges whatever the use of sponge_ts or not
OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(sal,temp)

#endif

   PRINT_DBG*, 'exit TS_OBC3D_APPLY_AGRIF routine'

  END SUBROUTINE ts_obc3d_apply_agrif

   !!======================================================================

  SUBROUTINE ts_obc2d_apply_agrif

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE ts_obc2d_apply_agrif  ***
   !&E
   !&E ** Purpose : apply2d open boundary conditions for temperature
   !&E              and salinity field for AGRIF
   !&E
   !&E ** Description : apply2d external value when water mass enters the domain
   !&E
   !&E ** Called by : stepu, stepv
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !                        Original code
   !&E       !  2007-05  (V. Garnier) Mise en forme - introduction de obc_* variables
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
#ifdef key_agrif

   USE parameters
   USE comvars2d, ONLY : c_agrif_interp
   USE comvarp2d
   USE agrifcommonid
   USE_MPI toolmpi, ONLY : ex_i_rsh,ex_j_rsh


   !! * Arguments

   !! * Local declarations
   External :: interp3dtraceur

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

   PRINT_DBG*, 'enter TS_OBC2D_APPLY_AGRIF routine'

   Agrif_UseSpecialValue = .TRUE.
   Agrif_SpecialValue = valmanq

   c_agrif_interp='traceur_bc'
   CALL Agrif_Bc_Variable(traceur_id,procname=interp3dtraceur)
   Agrif_UseSpecialValue = .FALSE.

! keep these exchanges whatever the use of sponge_ts or not
OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(sal,temp)

#endif

   PRINT_DBG*, 'exit TS_OBC2D_APPLY_AGRIF routine'

   END SUBROUTINE ts_obc2d_apply_agrif

   !!======================================================================

#ifdef key_agrif

  SUBROUTINE interp3dtraceur(tabres,l1,l2,k1,k2,i1,i2,j1,j2,before)

   USE parameters
   USE comvarp3d
   USE comvars2d
   USE comvarp2d
#if defined key_siggen || defined key_gencoord
   USE toolgeom, ONLY : f_hzu
   USE comsiggen,ONLY : hzex,hzey,hzdcusds
#endif
   IMPLICIT NONE

   INTEGER i1,i2,j1,j2,k1,k2,l1,l2
   REAL tabres(l1:l2,k1:k2,i1:i2,j1:j2)
   REAL tra_diff(l1:l2,k1:k2,i1:i2,j1:j2)
   REAL(KIND=rsh),DIMENSION(l1:l2,k1:k2,i1:i2,j1:j2) :: Fx, Fy
   INTEGER i,j,k
   LOGICAL :: before
   REAL(KIND=rsh),DIMENSION(i1:i2,j1:j2) :: coeff
   REAL(KIND=rsh) :: visc_sponge
   INTEGER :: dist
   REAL :: f_hzuinv

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

   PRINT_DBG*, 'enter INTERP3DTRACEUR routine'

   IF (before) THEN
     DO j=j1,j2
       DO i=i1,i2
         IF (h0(i,j)==-valmanq) THEN
           DO k=k1,k2
             tabres(1:2,k,i,j) = Agrif_SpecialValue
!            tabres(k,i,j) = 10.0
           END DO
         ELSE
           DO k=k1,k2
             tabres(1,k,i,j) = temp(k,i,j)
             tabres(2,k,i,j) = sal(k,i,j)
           END DO
         END IF

       END DO
     END DO
   ELSE
     SELECT CASE (c_agrif_interp)
       CASE('traceur_init')

         DO j=j1,j2
           DO i=i1,i2
             IF (h0(i,j)==-valmanq) THEN
               DO k=k1,k2
                 temp(k,i,j) = valmanq
                 sal(k,i,j)  = valmanq
               END DO
             ELSE
               DO k=k1,k2
                 temp(k,i,j) = tabres(1,k,i,j)
                 sal(k,i,j)  = tabres(2,k,i,j)
               END DO
             END IF
           END DO
         END DO

       CASE('traceur_bc')

         ! Initialize tra_diff : T/S (fine grid) - T/S (interpolated from mother to child grid)
         DO j=j1,j2
           DO i=i1,i2
!         IF (h0(i,j)+ssh(i,j)>=h0fond) THEN
             ! mask
             IF (h0(i,j)==-valmanq) THEN
               ! when the child is located along a boundary of the mother grid
               IF (i==imin+1.or.i==imax.or.j==jmax.or.j==jmin+1) THEN
                 temp(:,i,j) = valmanq
                 sal(:,i,j)  = valmanq
               END IF
               tra_diff(:,:,i,j) = 0.
             ELSE
               ! when the child is located along a boundary of the mother grid
               ! set interpolated field from mother to child grid
               IF (i==imin+1.or.i==imax.or.j==jmax.or.j==jmin+1) THEN
                 DO k=k1,k2
                   temp(k,i,j) = tabres(1,k,i,j)
                   sal(k,i,j)  = tabres(2,k,i,j)
                 END DO
               END IF
#if defined key_agrif_obc_sponge_ts
               ! T/S (fine grid) - T/S (interpolated from mother to child grid)
               tra_diff(1,:,i,j) = temp(:,i,j) - tabres(1,:,i,j)
               tra_diff(2,:,i,j) = sal(:,i,j) - tabres(2,:,i,j)
#endif
             END IF
           END DO
         END DO


#if ! defined key_agrif_obc_sponge_ts
         RETURN
#else

         ! Sponge layer for the child grid (over nb_connect*Agrif_Irhox)

         !cval WARNING AGRIF en dur !!!!! (depend on advection scheme ?)
         visc_sponge = 1./32.

         ! Estimate of the coefficient of diffusion (m2/s)
         coeff(:,:) = 0.0_rsh
         DO j=j1,j2
           DO i=i1,i2
             IF (h0(i,j)+ssh(i,j) > h0fond) THEN
               dist = nb_connect*Agrif_Irhox()
               dist = min(dist,i-(imin+1))
               dist = min(dist,imax-i)
               dist = min(dist,j-(jmin+1))
               dist = min(dist,jmax-j)
               coeff(i,j) = (1.-REAL(dist)/real(nb_connect*Agrif_Irhox()))*visc_sponge* &
                            (dx(i,j)*dy(i,j))**2/(halfdt*(dx(i,j)**2+dy(i,j)**2))
             END IF
           END DO
         END DO

         ! nu d TSdiff /dx
         Fx(:,:,:,:)=0.0_rsh
         DO j=MAX(j1,ljminp2),MIN(j2,ljmaxm1)
           DO i=MAX(i1,liminp1m1),MIN(i2-1,limaxm1)
             IF (h0(i,j)+ssh(i,j) > h0fond) THEN
               DO k=k1,k2
#if defined key_siggen || defined key_gencoord
                 Fx(:,k,i,j) = hzex(k,i,j)*0.5_rsh*(coeff(i,j)+coeff(i+1,j))* &
                               (tra_diff(:,k,i+1,j)-tra_diff(:,k,i,j))/dx_u(i,j)
#else
                 Fx(:,k,i,j) = hex(i,j)*0.5_rsh*(coeff(i,j)+coeff(i+1,j))* &
                               (tra_diff(:,k,i+1,j)-tra_diff(:,k,i,j))/dx_u(i,j)
#endif
               END DO
             END IF
           END DO
         END DO

         ! nu d TSdiff /dy
         Fy(:,:,:,:)=0.0_rsh
         DO j=MAX(j1,ljminp1m1),MIN(j2-1,ljmaxm1)
           DO i=MAX(i1,liminp2),MIN(i2,limaxm1)
             IF (h0(i,j)+ssh(i,j) > h0fond) THEN
               DO k=k1,k2
#if defined key_siggen || defined key_gencoord
                 Fy(:,k,i,j) = hzey(k,i,j)*0.5_rsh*(coeff(i,j)+coeff(i,j+1))* &
                              (tra_diff(:,k,i,j+1)-tra_diff(:,k,i,j))/dy_v(i,j)
#else
                 Fy(:,k,i,j) = hey(i,j)*0.5_rsh*(coeff(i,j)+coeff(i,j+1))* &
                              (tra_diff(:,k,i,j+1)-tra_diff(:,k,i,j))/dy_v(i,j)
#endif
               END DO
             END IF
           END DO
         END DO

         ! Add d / dx ( nu d TSdiff/dx ) and  d / dx ( nu d TSdiff/dy )

         DO j=MAX(j1+1,ljminp2),MIN(j2-1,ljmaxm1)
           DO i=MAX(i1+1,liminp2),MIN(i2-1,limaxm1)
             IF (h0(i,j)+ssh(i,j) > hm) THEN
               DO k=k1,k2
#if defined key_siggen || defined key_gencoord
                 f_hzuinv = halfdt/(ssh(i,j)+hzdcusds(k,i,j))
!                 f_hzuinv=halfdt/f_hzu(h0(i,j),ssh(i,j),k,i,j)
#else
                 f_hzuinv=halfdt/(ssh(i,j)+h0(i,j))
#endif
                 temp(k,i,j) = temp(k,i,j) &
                               + ((1./dx(i,j))*(Fx(1,k,i,j)-Fx(1,k,i-1,j))  &
                               +  (1./dy(i,j))*(Fy(1,k,i,j)-Fy(1,k,i,j-1))) &
                                *f_hzuinv

                 sal(k,i,j) = sal(k,i,j) &
                              + ((1./dx(i,j))*(Fx(2,k,i,j)-Fx(2,k,i-1,j))  &
                              +  (1./dy(i,j))*(Fy(2,k,i,j)-Fy(2,k,i,j-1))) &
                                *f_hzuinv
               END DO
             END IF
           END DO
         END DO

#endif /* key_agrif_obc_sponge_ts */

     END SELECT

   END IF

   PRINT_DBG*, 'exit INTERP3DTRACEUR routine'

   RETURN

  END SUBROUTINE interp3dtraceur
#endif
