Index: /NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_interp.F90
===================================================================
--- /NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_interp.F90	(revision 14012)
+++ /NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_interp.F90	(revision 14013)
@@ -940,34 +940,31 @@
          ELSE
          
+            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells 
+                                             ! linear vertical interpolation
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     !
+                     N_in  = mbkt(ji,jj)
+                     N_out = mbkt(ji,jj)
+                     z_in(1) = ptab(ji,jj,1,n2)
+                     tabin(1,1:jpts) = ptab(ji,jj,1,1:jpts)
+                     DO jk=2, N_in
+                        z_in(jk) = ptab(ji,jj,jk,n2)
+                        tabin(jk,1:jpts) = ptab(ji,jj,jk,1:jpts)
+                     END DO
+                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - ptab(ji,jj,k2,n2)
+                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Krhs_a)
+                     DO jk=2, N_out
+                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Krhs_a) + e3t(ji,jj,jk,Krhs_a))
+                     END DO
+                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Krhs_a)
+                     CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ptab(ji,jj,1:N_out,1:jpts), &
+                                   &   z_out(1:N_out),N_in,N_out,jpts)  
+                  END DO
+               END DO
+            ENDIF
+
             DO jn =1, jpts
-
-               IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells 
-                                                ! linear vertical interpolation
-                  DO jj=j1,j2
-                     DO ji=i1,i2
-                        !
-                        N_in  = mbkt(ji,jj)
-                        N_out = mbkt(ji,jj)
-                        z_in(1) = ptab(ji,jj,1,n2)
-                        tabin(1,1:jpts) = ptab(ji,jj,1,1:jpts)
-                        DO jk=2, N_in
-                           z_in(jk) = ptab(ji,jj,jk,n2)
-                           tabin(jk,1:jpts) = ptab(ji,jj,jk,1:jpts)
-                        END DO
-                        IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - ptab(ji,jj,k2,n2)
-                        z_out(1) = 0.5_wp * e3t(ji,jj,1,Krhs_a)
-                        DO jk=2, N_out
-                           z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Krhs_a) + e3t(ji,jj,jk,Krhs_a))
-                        END DO
-                        IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Krhs_a)
-                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ptab(ji,jj,1:N_out,1:jpts), &
-                                      &   z_out(1:N_out),N_in,N_out,jpts)  
-                     END DO
-                  END DO
-
-               ENDIF
-
                ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a) = ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk)
-
             END DO
          ENDIF
Index: /NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_top_interp.F90
===================================================================
--- /NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_top_interp.F90	(revision 14012)
+++ /NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_top_interp.F90	(revision 14013)
@@ -43,5 +43,5 @@
       Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = .TRUE.
-      l_vremap = ln_vert_remap
+      l_vremap              = ln_vert_remap
       !
       CALL Agrif_Bc_variable( trn_id, procname=interptrn )
@@ -64,7 +64,7 @@
       INTEGER  :: item
       ! vertical interpolation:
-      REAL(wp) :: zhtot
-      REAL(wp), DIMENSION(k1:k2,1:jptra) :: tabin
-      REAL(wp), DIMENSION(k1:k2) :: h_in, z_in
+      REAL(wp) :: zhtot, zwgt
+      REAL(wp), DIMENSION(k1:k2,1:jptra) :: tabin, tabin_i
+      REAL(wp), DIMENSION(k1:k2) :: z_in, h_in_i, z_in_i
       REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
       !!----------------------------------------------------------------------
@@ -85,32 +85,20 @@
          END DO
 
-         IF( l_vremap .OR. l_ini_child) THEN
-            ! Interpolate thicknesses
+         IF( l_vremap .OR. l_ini_child .OR. ln_zps ) THEN
+            ! Fill cell depths (i.e. gdept) to be interpolated
             ! Warning: these are masked, hence extrapolated prior interpolation.
-            DO jk=k1,k2
-               DO jj=j1,j2
-                  DO ji=i1,i2
-                      ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
-
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  ptab(ji,jj,k1,jptra+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kmm_a)
+                  DO jk=k1+1,k2
+                     ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * &
+                        & ( ptab(ji,jj,jk-1,jptra+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kmm_a)+e3t(ji,jj,jk,Kmm_a)) )
                   END DO
                END DO
             END DO
 
-            ! Extrapolate thicknesses in partial bottom cells:
-            ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
-            IF (ln_zps) THEN
-               DO jj=j1,j2
-                  DO ji=i1,i2
-                      jk = mbkt(ji,jj)
-                      ptab(ji,jj,jk,jptra+1) = 0._wp
-                  END DO
-               END DO           
-            END IF
-        
             ! Save ssh at last level:
             IF (.NOT.ln_linssh) THEN
                ptab(i1:i2,j1:j2,k2,jptra+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
-            ELSE
-               ptab(i1:i2,j1:j2,k2,jptra+1) = 0._wp
             END IF      
          ENDIF
@@ -126,32 +114,35 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tr(ji,jj,:,:,Krhs_a) = 0.                  
+                  tr(ji,jj,:,:,Krhs_a) = 0.  
+                  !
+                  ! Build vertical grids:
                   N_in = mbkt_parent(ji,jj)
-                  zhtot = 0._wp
-                  DO jk=1,N_in !k2 = jpk of parent grid
-                     IF (jk==N_in) THEN
-                        h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot
-                     ELSE
-                        h_in(jk) = ptab(ji,jj,jk,n2)
-                     ENDIF
-                     zhtot = zhtot + h_in(jk)
-                     tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
-                  END DO
-                  z_in(1) = 0.5_wp * h_in(1) - zhtot + ht0_parent(ji,jj)
+                  ! Input grid (account for partial cells if any):
+                  DO jk=1,N_in
+                     z_in(jk) = ptab(ji,jj,jk,n2) - ptab(ji,jj,k2,n2)
+                     tabin(jk,1:jptra) = ptab(ji,jj,jk,1:jptra)
+                  END DO
+                  
+                  ! Intermediate grid:
+                  DO jk = 1, N_in
+                     h_in_i(jk) = e3t0_parent(ji,jj,jk) * & 
+                       &       (1._wp + ptab(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj)))
+                  END DO
+                  z_in_i(1) = 0.5_wp * h_in_i(1)
                   DO jk=2,N_in
-                     z_in(jk) = z_in(jk-1) + 0.5_wp * (h_in(jk-1)+h_in(jk))
-                  END DO
-
-                  N_out = 0
-                  DO jk=1,jpk ! jpk of child grid
-                     IF (tmask(ji,jj,jk) == 0._wp) EXIT 
-                     N_out = N_out + 1
+                     z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) )
+                  END DO
+                  z_in_i(1:N_in) = z_in_i(1:N_in)  - ptab(ji,jj,k2,n2)
+
+                  ! Output (Child) grid:
+                  N_out = mbkt(ji,jj)
+                  DO jk=1,N_out
                      h_out(jk) = e3t(ji,jj,jk,Krhs_a)
                   END DO
-
-                  z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + ht_0(ji,jj)
+                  z_out(1) = 0.5_wp * h_out(1)
                   DO jk=2,N_out
-                     z_out(jk) = z_out(jk-1) + 0.5_wp * (h_out(jk-1)+h_out(jk))
-                  END DO
+                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) )
+                  END DO
+                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Krhs_a)               
 
                   IF (N_in*N_out > 0) THEN
@@ -159,7 +150,9 @@
                         CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tr(ji,jj,1:N_out,1:jptra,Krhs_a),          &
                                       &   z_out(1:N_out),N_in,N_out,jptra)  
-                     ELSE 
-                        CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in(1:N_in),tr(ji,jj,1:N_out,1:jptra,Krhs_a),   &
-                                      &   h_out(1:N_out),N_in,N_out,jptra)  
+                     ELSE  
+                        CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tabin_i(1:N_in,1:jptra),                     &
+                                     &   z_in_i(1:N_in),N_in,N_in,jptra)
+                        CALL reconstructandremap(tabin_i(1:N_in,1:jptra),h_in_i(1:N_in),tr(ji,jj,1:N_out,1:jptra,Krhs_a), &
+                                      &   h_out(1:N_out),N_in,N_out,jptra)   
                      ENDIF
                   ENDIF
@@ -169,5 +162,31 @@
  
          ELSE
-         
+
+            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells 
+                                             ! linear vertical interpolation
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     !
+                     N_in  = mbkt(ji,jj)
+                     N_out = mbkt(ji,jj)
+                     z_in(1) = ptab(ji,jj,1,n2)
+                     tabin(1,1:jptra) = ptab(ji,jj,1,1:jptra)
+                     DO jk=2, N_in
+                        z_in(jk) = ptab(ji,jj,jk,n2)
+                        tabin(jk,1:jptra) = ptab(ji,jj,jk,1:jptra)
+                     END DO
+                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - ptab(ji,jj,k2,n2)
+                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Krhs_a)
+                     DO jk=2, N_out
+                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Krhs_a) + e3t(ji,jj,jk,Krhs_a))
+                     END DO
+                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Krhs_a)
+                     CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),ptab(ji,jj,1:N_out,1:jptra), &
+                                   &   z_out(1:N_out),N_in,N_out,jptra)  
+                  END DO
+               END DO
+
+            ENDIF
+
             DO jn=1, jptra
                 tr(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
Index: /NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_top_sponge.F90
===================================================================
--- /NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_top_sponge.F90	(revision 14012)
+++ /NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_top_sponge.F90	(revision 14013)
@@ -45,12 +45,16 @@
       !
 #if defined SPONGE_TOP
-!! Assume persistence 
+      !! Assume persistence:
       zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
-      CALL Agrif_sponge
+
       Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = .TRUE.
+      l_vremap              = ln_vert_remap
       tabspongedone_trn     = .FALSE.
+      !
       CALL Agrif_Bc_Variable( trn_sponge_id, calledweight=zcoef, procname=interptrn_sponge )
+      !
       Agrif_UseSpecialValue = .FALSE.
+      l_vremap              = .FALSE.
 #endif
       !
@@ -58,7 +62,7 @@
 
 
-   SUBROUTINE interptrn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
-      !!----------------------------------------------------------------------
-      !!                   *** ROUTINE interptrn_sponge ***
+   SUBROUTINE interptrn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before) 
+      !!----------------------------------------------------------------------
+      !!                 *** ROUTINE interptrn_sponge ***
       !!----------------------------------------------------------------------
       INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
@@ -68,12 +72,12 @@
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
       INTEGER  ::   iku, ikv
-      REAL(wp) ::   ztra, zabe1, zabe2, zbtr
-      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk)  ::   ztu, ztv
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::   trbdiff
+      REAL(wp) :: ztra, zabe1, zabe2, zbtr, zhtot
+      REAL(wp), DIMENSION(i1-1:i2,j1-1:j2,jpk) :: ztu, ztv
+      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::trbdiff
       ! vertical interpolation:
       REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
-      REAL(wp), DIMENSION(k1:k2,1:jptra) :: tabin
-      REAL(wp), DIMENSION(k1:k2) :: h_in
-      REAL(wp), DIMENSION(1:jpk) :: h_out
+      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin, tabin_i
+      REAL(wp), DIMENSION(k1:k2) :: z_in, z_in_i, h_in_i
+      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
       INTEGER :: N_in, N_out
       !!----------------------------------------------------------------------
@@ -90,64 +94,140 @@
          END DO
 
-# if defined key_vertical
-         DO jk=k1,k2
-            DO jj=j1,j2
-               DO ji=i1,i2
-                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a)
-               END DO
-            END DO
-         END DO
-# endif
-      ELSE      
-# if defined key_vertical
-         tabres_child(:,:,:,:) = 0.
-         DO jj=j1,j2
-            DO ji=i1,i2
-               N_in = 0
-               DO jk=k1,k2 !k2 = jpk of parent grid
-                  IF (tabres(ji,jj,jk,n2) == 0) EXIT
-                  N_in = N_in + 1
-                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)
-                  h_in(N_in) = tabres(ji,jj,jk,n2)
-               END DO
-               N_out = 0
-               DO jk=1,jpk ! jpk of child grid
-                  IF (tmask(ji,jj,jk) == 0) EXIT 
-                  N_out = N_out + 1
-                  h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above
-               ENDDO
-               IF (N_in > 0) THEN
-                  CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in,tabres_child(ji,jj,1:N_out,1:jptra),h_out,N_in,N_out,jptra)
-               ENDIF
-            ENDDO
-         ENDDO
-# endif
-
-         DO jj=j1,j2
-            DO ji=i1,i2
-               DO jk=1,jpkm1
-# if defined key_vertical
-                  trbdiff(ji,jj,jk,1:jptra) = ( tr(ji,jj,jk,1:jptra,Kbb_a) - tabres_child(ji,jj,jk,1:jptra) ) * tmask(ji,jj,jk)
-# else
-                  trbdiff(ji,jj,jk,1:jptra) = ( tr(ji,jj,jk,1:jptra,Kbb_a) - tabres(ji,jj,jk,1:jptra)) * tmask(ji,jj,jk)
-# endif
-               ENDDO
-            ENDDO
-         ENDDO
-
-         DO jn = 1, jptra
+         IF ( l_vremap.OR.ln_zps ) THEN
+
+            ! Fill cell depths (i.e. gdept) to be interpolated
+            ! Warning: these are masked, hence extrapolated prior interpolation.
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  tabres(ji,jj,k1,jptra+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kbb_a)
+                  DO jk=k1+1,k2
+                     tabres(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * &
+                        & ( tabres(ji,jj,jk-1,jptra+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a)+e3t(ji,jj,jk,Kbb_a)) )
+                  END DO
+               END DO
+            END DO
+
+            ! Save ssh at last level:
+            IF ( .NOT.ln_linssh ) THEN
+               tabres(i1:i2,j1:j2,k2,jptra+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1) 
+            END IF  
+    
+         END IF
+
+      ELSE 
+         !
+         IF ( l_vremap ) THEN
+
+            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp
+
+            DO jj=j1,j2
+               DO ji=i1,i2
+
+                  tabres_child(ji,jj,:,:) = 0._wp 
+                  ! Build vertical grids:
+                  N_in = mbkt_parent(ji,jj)
+                  ! Input grid (account for partial cells if any):
+                  DO jk=1,N_in
+                     z_in(jk) = tabres(ji,jj,jk,n2) - tabres(ji,jj,k2,n2)
+                     tabin(jk,1:jptra) = tabres(ji,jj,jk,1:jptra)
+                  END DO
+                  
+                  ! Intermediate grid:
+                  DO jk = 1, N_in
+                     h_in_i(jk) = e3t0_parent(ji,jj,jk) * & 
+                       &       (1._wp + tabres(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj)))
+                  END DO
+                  z_in_i(1) = 0.5_wp * h_in_i(1)
+                  DO jk=2,N_in
+                     z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) )
+                  END DO
+                  z_in_i(1:N_in) = z_in_i(1:N_in)  - tabres(ji,jj,k2,n2)
+
+                  ! Output (Child) grid:
+                  N_out = mbkt(ji,jj)
+                  DO jk=1,N_out
+                     h_out(jk) = e3t(ji,jj,jk,Kbb_a)
+                  END DO
+                  z_out(1) = 0.5_wp * h_out(1)
+                  DO jk=2,N_out
+                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) )
+                  END DO
+                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Kbb_a)
+
+                  ! Account for small differences in the free-surface
+                  IF ( sum(h_out(1:N_out)) > sum(h_in_i(1:N_in) )) THEN
+                     h_out(1) = h_out(1)  - ( sum(h_out(1:N_out))-sum(h_in_i(1:N_in)) )
+                  ELSE
+                     h_in_i(1)= h_in_i(1) - ( sum(h_in_i(1:N_in))-sum(h_out(1:N_out)) )
+                  END IF
+                  IF (N_in*N_out > 0) THEN
+                     CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tabin_i(1:N_in,1:jptra),z_in_i(1:N_in),N_in,N_in,jptra)
+                     CALL reconstructandremap(tabin_i(1:N_in,1:jptra),h_in_i(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),h_out(1:N_out),N_in,N_out,jptra)
+!                     CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),z_out(1:N_in),N_in,N_out,jptra)  
+                  ENDIF
+               END DO
+            END DO
+
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  DO jk=1,jpkm1
+                     trbdiff(ji,jj,jk,1:jptra) = (tr(ji,jj,jk,1:jptra,Kbb_a) - tabres_child(ji,jj,jk,1:jptra)) * tmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+
+         ELSE
+
+            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells
+
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     !
+                     N_in  = mbkt(ji,jj) 
+                     N_out = mbkt(ji,jj) 
+                     z_in(1) = tabres(ji,jj,1,n2)
+                     tabin(1,1:jptra) = tabres(ji,jj,1,1:jptra)
+                     DO jk=2, N_in
+                        z_in(jk) = tabres(ji,jj,jk,n2)
+                        tabin(jk,1:jptra) = tabres(ji,jj,jk,1:jptra)
+                     END DO 
+                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - tabres(ji,jj,k2,n2)
+
+                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Kbb_a)
+                     DO jk=2, N_out
+                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk,Kbb_a)) 
+                     END DO 
+                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Kbb_a)
+
+                     CALL remap_linear(tabin(1:N_in,1:jptra), z_in(1:N_in), tabres(ji,jj,1:N_out,1:jptra), &
+                                         &   z_out(1:N_out), N_in, N_out, jptra)
+                  END DO
+               END DO
+            ENDIF
+
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  DO jk=1,jpkm1
+                     trbdiff(ji,jj,jk,1:jptra) = (tr(ji,jj,jk,1:jptra,Kbb_a) - tabres(ji,jj,jk,1:jptra))*tmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+
+         END IF
+
+         DO jn = 1, jptra            
             DO jk = 1, jpkm1
-               ztu(i1:i2,j1:j2,jk) = 0._wp
+               ztu(i1-1:i2,j1-1:j2,jk) = 0._wp
                DO jj = j1,j2
                   DO ji = i1,i2-1
-                     zabe1 = rn_sponge_tra * r1_Dt * fspu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
-                     ztu(ji,jj,jk) = zabe1 * ( trbdiff(ji+1,jj,jk,jn) - trbdiff(ji,jj,jk,jn) )
-                  END DO
-               END DO
-               ztv(i1:i2,j1:j2,jk) = 0._wp
+                     zabe1 = rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
+                     ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( trbdiff(ji+1,jj  ,jk,jn) - trbdiff(ji,jj,jk,jn) ) 
+                  END DO
+               END DO
+               ztv(i1-1:i2,j1-1:j2,jk) = 0._wp
                DO ji = i1,i2
                   DO jj = j1,j2-1
-                     zabe2 = rn_sponge_tra * r1_Dt * fspv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
-                     ztv(ji,jj,jk) = zabe2 * ( trbdiff(ji,jj+1,jk,jn) - trbdiff(ji,jj,jk,jn) )
+                     zabe2 = rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
+                     ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( trbdiff(ji  ,jj+1,jk,jn) - trbdiff(ji,jj,jk,jn) )
                   END DO
                END DO
@@ -166,12 +246,14 @@
             END DO
             !
+! JC: there is something wrong with the Laplacian in corners
             DO jk = 1, jpkm1
-               DO jj = j1+1,j2-1
-                  DO ji = i1+1,i2-1
-                     IF( .NOT. tabspongedone_trn(ji,jj) ) THEN
+               DO jj = j1,j2
+                  DO ji = i1,i2
+                     IF (.NOT. tabspongedone_trn(ji,jj)) THEN 
                         zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
-                                                ! horizontal diffusive trends
-                        ztra = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) &
-                             &  - rn_trelax_tra * r1_Dt * fspt(ji,jj) * trbdiff(ji,jj,jk,jn)
+                        ! horizontal diffusive trends
+                        ztra = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   & 
+                          &           + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) &
+                          &   - rn_trelax_tra * r1_Dt * fspt(ji,jj) * trbdiff(ji,jj,jk,jn)
                         ! add it to the general tracer trends
                         tr(ji,jj,jk,jn,Krhs_a) = tr(ji,jj,jk,jn,Krhs_a) + ztra
@@ -179,10 +261,13 @@
                   END DO
                END DO
-            END DO
+
+            END DO
+            !
          END DO
          !
-         tabspongedone_trn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
+         tabspongedone_trn(i1:i2,j1:j2) = .TRUE.
+         !
       ENDIF
-      !                 
+      !
    END SUBROUTINE interptrn_sponge
 
