Index: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90	(revision 5950)
+++ branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90	(revision 5951)
@@ -1,5 +1,5 @@
 #define SPONGE && define SPONGE_TOP
 
-Module agrif_opa_sponge
+MODULE agrif_opa_sponge
 #if defined key_agrif  && ! defined key_offline
    USE par_oce
@@ -9,11 +9,13 @@
    USE agrif_oce
    USE wrk_nemo  
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn, interptsn, interpun, interpvn
-
-  !! * Substitutions
+   PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
+   PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
+
+   !! * Substitutions
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
@@ -23,5 +25,5 @@
    !!----------------------------------------------------------------------
 
-   CONTAINS
+CONTAINS
 
    SUBROUTINE Agrif_Sponge_Tra
@@ -30,55 +32,17 @@
       !!---------------------------------------------
       !!
-      INTEGER :: ji,jj,jk,jn
       REAL(wp) :: timecoeff
-      REAL(wp) :: ztsa, zabe1, zabe2, zbtr
-      REAL(wp), POINTER, DIMENSION(:,:    ) :: ztu, ztv
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztab
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: tsbdiff
 
 #if defined SPONGE
-      CALL wrk_alloc( jpi, jpj, ztu, ztv )
-      CALL wrk_alloc( jpi, jpj, jpk, jpts, ztab, tsbdiff  )
-
       timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
 
+      CALL Agrif_Sponge
       Agrif_SpecialValue=0.
       Agrif_UseSpecialValue = .TRUE.
-      ztab = 0.e0
-      CALL Agrif_Bc_Variable(ztab, tsa_id,calledweight=timecoeff,procname=interptsn)
+      tabspongedone_tsn = .FALSE.
+
+      CALL Agrif_Bc_Variable(tsn_sponge_id,calledweight=timecoeff,procname=interptsn_sponge)
+
       Agrif_UseSpecialValue = .FALSE.
-
-      tsbdiff(:,:,:,:) = tsb(:,:,:,:) - ztab(:,:,:,:)
-
-      CALL Agrif_Sponge
-
-      DO jn = 1, jpts
-         DO jk = 1, jpkm1
-            !
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  zabe1 = umask(ji,jj,jk) * spe1ur(ji,jj) * fse3u(ji,jj,jk)
-                  zabe2 = vmask(ji,jj,jk) * spe2vr(ji,jj) * fse3v(ji,jj,jk)
-                  ztu(ji,jj) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) )
-                  ztv(ji,jj) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
-               ENDDO
-            ENDDO
-
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk)
-                  ! horizontal diffusive trends
-                  ztsa = zbtr * (  ztu(ji,jj) - ztu(ji-1,jj  )   &
-                  &              + ztv(ji,jj) - ztv(ji  ,jj-1)  )
-                  ! add it to the general tracer trends
-                  tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
-               END DO
-            END DO
-            !
-         ENDDO
-      ENDDO
-
-      CALL wrk_dealloc( jpi, jpj, ztu, ztv )
-      CALL wrk_dealloc( jpi, jpj, jpk, jpts, ztab, tsbdiff  )
 #endif
 
@@ -90,89 +54,21 @@
       !!---------------------------------------------
       !!
-      INTEGER :: ji,jj,jk
       REAL(wp) :: timecoeff
-      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ubdiff, vbdiff
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: rotdiff, hdivdiff
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztab
 
 #if defined SPONGE
-      CALL wrk_alloc( jpi, jpj, jpk, ztab, ubdiff, vbdiff, rotdiff, hdivdiff )
-
       timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
 
       Agrif_SpecialValue=0.
       Agrif_UseSpecialValue = ln_spc_dyn
-      ztab = 0.e0
-      CALL Agrif_Bc_Variable(ztab, ua_id,calledweight=timecoeff,procname=interpun)
+
+      tabspongedone_u = .FALSE.
+      tabspongedone_v = .FALSE.         
+      CALL Agrif_Bc_Variable(un_sponge_id,calledweight=timecoeff,procname=interpun_sponge)
+
+      tabspongedone_u = .FALSE.
+      tabspongedone_v = .FALSE.
+      CALL Agrif_Bc_Variable(vn_sponge_id,calledweight=timecoeff,procname=interpvn_sponge)
+
       Agrif_UseSpecialValue = .FALSE.
-
-      ubdiff(:,:,:) = ( ub(:,:,:) - ztab(:,:,:) ) * umask(:,:,:)
-
-      ztab = 0.e0
-      Agrif_SpecialValue=0.
-      Agrif_UseSpecialValue = ln_spc_dyn
-      CALL Agrif_Bc_Variable(ztab, va_id,calledweight=timecoeff,procname=interpvn)
-      Agrif_UseSpecialValue = .FALSE.
-
-      vbdiff(:,:,:) = ( vb(:,:,:) - ztab(:,:,:) ) * vmask(:,:,:)
-
-      CALL Agrif_Sponge
-
-      DO jk = 1,jpkm1
-         ubdiff(:,:,jk) = ubdiff(:,:,jk) * spe1ur2(:,:)
-         vbdiff(:,:,jk) = vbdiff(:,:,jk) * spe2vr2(:,:)
-      ENDDO
-      
-      hdivdiff = 0.
-      rotdiff = 0.
-
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-
-         !                                             ! --------
-         ! Horizontal divergence                       !   div
-         !                                             ! --------
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1   ! vector opt.
-               zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk)
-               hdivdiff(ji,jj,jk) =  (  e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * ubdiff(ji  ,jj  ,jk)     &
-                  &                   - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * ubdiff(ji-1,jj  ,jk)     &
-                  &                   + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * vbdiff(ji  ,jj  ,jk)     &
-                  &                   - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * vbdiff(ji  ,jj-1,jk)  ) * zbtr
-            END DO
-         END DO
-
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1   ! vector opt.
-               zbtr = spbtr3(ji,jj) * fse3f(ji,jj,jk)
-               rotdiff(ji,jj,jk) = (  e2v(ji+1,jj  ) * vbdiff(ji+1,jj  ,jk) - e2v(ji,jj) * vbdiff(ji,jj,jk)    &
-                  &                 - e1u(ji  ,jj+1) * ubdiff(ji  ,jj+1,jk) + e1u(ji,jj) * ubdiff(ji,jj,jk)  ) &
-                  &               * fmask(ji,jj,jk) * zbtr
-            END DO
-         END DO
-
-      ENDDO
-
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1   ! vector opt.
-               ! horizontal diffusive trends
-               zua = - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
-                     + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk) ) / e1u(ji,jj)
-
-               zva = + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
-                     + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) / e2v(ji,jj)
-               ! add it to the general momentum trends
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
-               va(ji,jj,jk) = va(ji,jj,jk) + zva
-            END DO
-         END DO
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-      CALL wrk_dealloc( jpi, jpj, jpk, ztab, ubdiff, vbdiff, rotdiff, hdivdiff )
 #endif
 
@@ -199,11 +95,10 @@
          CALL wrk_alloc( jpi, jpj, ztabramp )
 
-         ispongearea  = 2 + 2 * Agrif_irhox()
+         ispongearea  = 2 + nn_sponge_len * Agrif_irhox()
          ilci = nlci - ispongearea
          ilcj = nlcj - ispongearea 
          z1spongearea = 1._wp / REAL( ispongearea - 2 )
-         spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
-
-         ztabramp(:,:) = 0.
+
+         ztabramp(:,:) = 0._wp
 
          IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
@@ -254,54 +149,15 @@
       ! Tracers
       IF( .NOT. spongedoneT ) THEN
-         spe1ur(:,:) = 0.
-         spe2vr(:,:) = 0.
-
-         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
-            spe1ur(2:ispongearea-1,:       ) = visc_tra                                        &
-               &                             *    0.5 * (  ztabramp(2:ispongearea-1,:      )   &
-               &                                         + ztabramp(3:ispongearea  ,:      ) ) &
-               &                             * e2u(2:ispongearea-1,:) / e1u(2:ispongearea-1,:)
-
-            spe2vr(2:ispongearea  ,1:jpjm1 ) = visc_tra                                        &
-               &                             *    0.5 * (  ztabramp(2:ispongearea  ,1:jpjm1)   &
-               &                                         + ztabramp(2:ispongearea,2  :jpj  ) ) &
-               &                             * e1v(2:ispongearea,1:jpjm1) / e2v(2:ispongearea,1:jpjm1)
-         ENDIF
-
-         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
-            spe1ur(ilci+1:nlci-2,:        ) = visc_tra                                   &
-               &                            * 0.5 * (  ztabramp(ilci+1:nlci-2,:      )   & 
-               &                                     + ztabramp(ilci+2:nlci-1,:      ) ) &
-               &                            * e2u(ilci+1:nlci-2,:) / e1u(ilci+1:nlci-2,:)
-
-            spe2vr(ilci+1:nlci-1,1:jpjm1  )  = visc_tra                                  &
-               &                            * 0.5 * (  ztabramp(ilci+1:nlci-1,1:jpjm1)   & 
-               &                                     + ztabramp(ilci+1:nlci-1,2:jpj  ) ) & 
-               &                            * e1v(ilci+1:nlci-1,1:jpjm1) / e2v(ilci+1:nlci-1,1:jpjm1)
-         ENDIF
-
-         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
-            spe1ur(1:jpim1,2:ispongearea  ) = visc_tra                                     &
-               &                            * 0.5 * (  ztabramp(1:jpim1,2:ispongearea  )   & 
-               &                                     + ztabramp(2:jpi  ,2:ispongearea  ) ) &
-               &                            * e2u(1:jpim1,2:ispongearea) / e1u(1:jpim1,2:ispongearea)
-   
-            spe2vr(:      ,2:ispongearea-1) = visc_tra                                     &
-               &                            * 0.5 * (  ztabramp(:      ,2:ispongearea-1)   &
-               &                                     + ztabramp(:      ,3:ispongearea  ) ) &
-               &                            * e1v(:,2:ispongearea-1) / e2v(:,2:ispongearea-1)
-         ENDIF
-
-         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
-            spe1ur(1:jpim1,ilcj+1:nlcj-1) = visc_tra                                   &
-               &                          * 0.5 * (  ztabramp(1:jpim1,ilcj+1:nlcj-1)   &
-               &                                   + ztabramp(2:jpi  ,ilcj+1:nlcj-1) ) &
-               &                                * e2u(1:jpim1,ilcj+1:nlcj-1) / e1u(1:jpim1,ilcj+1:nlcj-1)
-
-            spe2vr(:      ,ilcj+1:nlcj-2) = visc_tra                                   &
-               &                          * 0.5 * (  ztabramp(:      ,ilcj+1:nlcj-2)   &
-               &                                   + ztabramp(:      ,ilcj+2:nlcj-1) ) &
-               &                                * e1v(:,ilcj+1:nlcj-2) / e2v(:,ilcj+1:nlcj-2)
-         ENDIF
+         fsaht_spu(:,:) = 0._wp
+         fsaht_spv(:,:) = 0._wp
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1   ! vector opt.
+               fsaht_spu(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji+1,jj  ))
+               fsaht_spv(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji  ,jj+1))
+            END DO
+         END DO
+
+         CALL lbc_lnk( fsaht_spu, 'U', 1. )   ! Lateral boundary conditions
+         CALL lbc_lnk( fsaht_spv, 'V', 1. )
          spongedoneT = .TRUE.
       ENDIF
@@ -309,44 +165,17 @@
       ! Dynamics
       IF( .NOT. spongedoneU ) THEN
-         spe1ur2(:,:) = 0.
-         spe2vr2(:,:) = 0.
-
-         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
-            spe1ur2(2:ispongearea-1,:      ) = visc_dyn                                   &
-               &                             * 0.5 * (  ztabramp(2:ispongearea-1,:      ) &
-               &                                      + ztabramp(3:ispongearea  ,:      ) )
-            spe2vr2(2:ispongearea  ,1:jpjm1) = visc_dyn                                   &
-               &                             * 0.5 * (  ztabramp(2:ispongearea  ,1:jpjm1) &
-               &                                      + ztabramp(2:ispongearea  ,2:jpj  ) ) 
-         ENDIF
-
-         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
-            spe1ur2(ilci+1:nlci-2  ,:      ) = visc_dyn                                   &
-               &                             * 0.5 * (  ztabramp(ilci+1:nlci-2, :       ) &
-               &                                      + ztabramp(ilci+2:nlci-1, :       ) )                      
-            spe2vr2(ilci+1:nlci-1  ,1:jpjm1) = visc_dyn                                   &
-               &                             * 0.5 * (  ztabramp(ilci+1:nlci-1,1:jpjm1  ) &
-               &                                      + ztabramp(ilci+1:nlci-1,2:jpj    ) ) 
-         ENDIF
-
-         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
-            spe1ur2(1:jpim1,2:ispongearea  ) = visc_dyn                                   &  
-               &                             * 0.5 * (  ztabramp(1:jpim1,2:ispongearea  ) &
-               &                                      + ztabramp(2:jpi  ,2:ispongearea  ) ) 
-            spe2vr2(:      ,2:ispongearea-1) = visc_dyn                                   &
-               &                             * 0.5 * (  ztabramp(:      ,2:ispongearea-1) &
-               &                                      + ztabramp(:      ,3:ispongearea  ) )
-         ENDIF
-
-         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
-            spe1ur2(1:jpim1,ilcj+1:nlcj-1  ) = visc_dyn                                   &
-               &                             * 0.5 * (  ztabramp(1:jpim1,ilcj+1:nlcj-1  ) &
-               &                                      + ztabramp(2:jpi  ,ilcj+1:nlcj-1  ) ) 
-            spe2vr2(:      ,ilcj+1:nlcj-2  ) = visc_dyn                                   &
-               &                             * 0.5 * (  ztabramp(:      ,ilcj+1:nlcj-2  ) &
-               &                                      + ztabramp(:      ,ilcj+2:nlcj-1  ) )
-         ENDIF
+         fsahm_spt(:,:) = 0._wp
+         fsahm_spf(:,:) = 0._wp
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1   ! vector opt.
+               fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj)
+               fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji,jj) + ztabramp(ji  ,jj+1) &
+                                                     &  +ztabramp(ji,jj) + ztabramp(ji+1,jj  ) )
+            END DO
+         END DO
+
+         CALL lbc_lnk( fsahm_spt, 'T', 1. )   ! Lateral boundary conditions
+         CALL lbc_lnk( fsahm_spf, 'F', 1. )
          spongedoneU = .TRUE.
-         spbtr3(:,:) = 1. / ( e1f(:,:) * e2f(:,:) )
       ENDIF
       !
@@ -357,36 +186,277 @@
    END SUBROUTINE Agrif_Sponge
 
-   SUBROUTINE interptsn(tabres,i1,i2,j1,j2,k1,k2,n1,n2)
-      !!---------------------------------------------
-      !!   *** ROUTINE interptsn ***
+   SUBROUTINE interptsn_sponge(tabres,i1,i2,j1,j2,k1,k2,n1,n2,before)
+      !!---------------------------------------------
+      !!   *** ROUTINE interptsn_sponge ***
       !!---------------------------------------------
       INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
-
-      tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
-
-   END SUBROUTINE interptsn
-
-   SUBROUTINE interpun(tabres,i1,i2,j1,j2,k1,k2)
-      !!---------------------------------------------
-      !!   *** ROUTINE interpun ***
-      !!---------------------------------------------
+      LOGICAL, INTENT(in) :: before
+
+
+      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
+      INTEGER  ::   iku, ikv
+      REAL(wp) :: ztsa, zabe1, zabe2, zbtr
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ztu, ztv
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::tsbdiff
+      !
+      IF (before) THEN
+         tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
+      ELSE   
+   
+         tsbdiff(:,:,:,:) = tsb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:)    
+         DO jn = 1, jpts            
+            DO jk = 1, jpkm1
+               DO jj = j1,j2-1
+                  DO ji = i1,i2-1
+                     zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk)
+                     zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk)
+                     ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
+                     ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
+                  ENDDO
+               ENDDO
+
+               IF( ln_zps ) THEN      ! set gradient at partial step level
+                  DO jj = j1,j2-1
+                     DO ji = i1,i2-1
+                        ! last level
+                        iku = mbku(ji,jj)
+                        ikv = mbkv(ji,jj)
+                        IF( iku == jk ) THEN
+                           ztu(ji,jj,jk) = 0._wp
+                        ENDIF
+                        IF( ikv == jk ) THEN
+                           ztv(ji,jj,jk) = 0._wp
+                        ENDIF
+                     END DO
+                  END DO
+               ENDIF
+            ENDDO
+
+            DO jk = 1, jpkm1
+               DO jj = j1+1,j2-1
+                  DO ji = i1+1,i2-1
+
+                     IF (.NOT. tabspongedone_tsn(ji,jj)) THEN 
+                        zbtr = r1_e1e2t(ji,jj) / fse3t_n(ji,jj,jk)
+                        ! horizontal diffusive trends
+                        ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  )
+                        ! add it to the general tracer trends
+                        tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
+                     ENDIF
+
+                  ENDDO
+               ENDDO
+
+            ENDDO
+         ENDDO
+
+         tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
+
+      ENDIF
+
+   END SUBROUTINE interptsn_sponge
+
+   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2, before)
+      !!---------------------------------------------
+      !!   *** ROUTINE interpun_sponge ***
+      !!---------------------------------------------    
       INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
-
-      tabres(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2)
-
-   END SUBROUTINE interpun
-
-   SUBROUTINE interpvn(tabres,i1,i2,j1,j2,k1,k2)
-      !!---------------------------------------------
-      !!   *** ROUTINE interpvn ***
-      !!---------------------------------------------
+      LOGICAL, INTENT(in) :: before
+
+      INTEGER :: ji,jj,jk
+
+      ! sponge parameters 
+      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ubdiff
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
+      INTEGER :: jmax
+      !
+
+
+      IF (before) THEN
+         tabres = un(i1:i2,j1:j2,:)
+      ELSE
+
+         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:))*umask(i1:i2,j1:j2,:)
+
+         DO jk = 1, jpkm1                                 ! Horizontal slab
+            !                                             ! ===============
+
+            !                                             ! --------
+            ! Horizontal divergence                       !   div
+            !                                             ! --------
+            DO jj = j1,j2
+               DO ji = i1+1,i2   ! vector opt.
+                  zbtr = r1_e1e2t(ji,jj) / fse3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
+                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*fse3u_n(ji  ,jj,jk) * ubdiff(ji  ,jj,jk) &
+                                     &   -e2u(ji-1,jj)*fse3u_n(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr
+               END DO
+            END DO
+
+            DO jj = j1,j2-1
+               DO ji = i1,i2   ! vector opt.
+                  zbtr = r1_e1e2f(ji,jj) * fse3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
+                  rotdiff(ji,jj,jk) = (-e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &
+                                       +e1u(ji,jj  ) * ubdiff(ji,jj  ,jk) & 
+                                    & ) * fmask(ji,jj,jk) * zbtr 
+               END DO
+            END DO
+         ENDDO
+
+         !
+
+
+
+         DO jj = j1+1, j2-1
+            DO ji = i1+1, i2-1   ! vector opt.
+
+               IF (.NOT. tabspongedone_u(ji,jj)) THEN
+                  DO jk = 1, jpkm1                                 ! Horizontal slab
+                     ze2u = rotdiff (ji,jj,jk)
+                     ze1v = hdivdiff(ji,jj,jk)
+                     ! horizontal diffusive trends
+                     zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * fse3u_n(ji,jj,jk) )   &
+                           + ( hdivdiff(ji+1,jj,jk) - ze1v  ) / e1u(ji,jj)
+
+                     ! add it to the general momentum trends
+                     ua(ji,jj,jk) = ua(ji,jj,jk) + zua
+
+                  END DO
+               ENDIF
+
+            END DO
+         END DO
+
+         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
+
+         jmax = j2-1
+         IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-3)
+
+         DO jj = j1+1, jmax
+            DO ji = i1+1, i2   ! vector opt.
+
+               IF (.NOT. tabspongedone_v(ji,jj)) THEN
+                  DO jk = 1, jpkm1                                 ! Horizontal slab
+                     ze2u = rotdiff (ji,jj,jk)
+                     ze1v = hdivdiff(ji,jj,jk)
+
+                     ! horizontal diffusive trends
+                     zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * fse3v_n(ji,jj,jk) )   &
+                           + ( hdivdiff(ji,jj+1,jk) - ze1v  ) / e2v(ji,jj)
+
+                     ! add it to the general momentum trends
+                     va(ji,jj,jk) = va(ji,jj,jk) + zva
+                  END DO
+               ENDIF
+
+            END DO
+         END DO
+
+
+         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
+
+      ENDIF
+
+
+   END SUBROUTINE interpun_sponge
+
+
+   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2, before,nb,ndir)
+      !!---------------------------------------------
+      !!   *** ROUTINE interpvn_sponge ***
+      !!--------------------------------------------- 
       INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
-
-      tabres(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2)
-
-   END SUBROUTINE interpvn
+      LOGICAL, INTENT(in) :: before
+      INTEGER, INTENT(in) :: nb , ndir
+
+      INTEGER :: ji,jj,jk
+
+      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
+
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: vbdiff
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
+      INTEGER :: imax
+      !
+
+      IF (before) THEN 
+         tabres = vn(i1:i2,j1:j2,:)
+      ELSE
+
+         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:)
+
+         DO jk = 1, jpkm1                                 ! Horizontal slab
+            !                                             ! ===============
+
+            !                                             ! --------
+            ! Horizontal divergence                       !   div
+            !                                             ! --------
+            DO jj = j1+1,j2
+               DO ji = i1,i2   ! vector opt.
+                  zbtr = r1_e1e2t(ji,jj) / fse3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
+                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  &
+                                     &  -e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr
+               END DO
+            END DO
+            DO jj = j1,j2
+               DO ji = i1,i2-1   ! vector opt.
+                  zbtr = r1_e1e2f(ji,jj) * fse3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
+                  rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 
+                                    &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk) &
+                                    & ) * fmask(ji,jj,jk) * zbtr
+               END DO
+            END DO
+         ENDDO
+
+         !                                                ! ===============
+         !                                                
+
+         imax = i2-1
+         IF ((nbondi == 1).OR.(nbondi == 2)) imax = MIN(imax,nlci-3)
+
+         DO jj = j1+1, j2
+            DO ji = i1+1, imax   ! vector opt.
+               IF (.NOT. tabspongedone_u(ji,jj)) THEN
+                  DO jk = 1, jpkm1                                 ! Horizontal slab
+                     ze2u = rotdiff (ji,jj,jk)
+                     ze1v = hdivdiff(ji,jj,jk)
+                     ! horizontal diffusive trends
+                     zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * fse3u_n(ji,jj,jk) ) + ( hdivdiff(ji+1,jj,jk) - ze1v) &
+                           / e1u(ji,jj)
+
+
+                     ! add it to the general momentum trends
+                     ua(ji,jj,jk) = ua(ji,jj,jk) + zua
+                  END DO
+
+               ENDIF
+            END DO
+         END DO
+
+         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
+
+         DO jj = j1+1, j2-1
+            DO ji = i1+1, i2-1   ! vector opt.
+               IF (.NOT. tabspongedone_v(ji,jj)) THEN
+                  DO jk = 1, jpkm1                                 ! Horizontal slab
+                     ze2u = rotdiff (ji,jj,jk)
+                     ze1v = hdivdiff(ji,jj,jk)
+                     ! horizontal diffusive trends
+
+                     zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * fse3v_n(ji,jj,jk) ) + ( hdivdiff(ji,jj+1,jk) - ze1v) &
+                           / e2v(ji,jj)
+
+                     ! add it to the general momentum trends
+                     va(ji,jj,jk) = va(ji,jj,jk) + zva
+                  END DO
+               ENDIF
+            END DO
+         END DO
+         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
+      ENDIF
+
+   END SUBROUTINE interpvn_sponge
 
 #else
