Index: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_sponge.F90
===================================================================
--- NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_sponge.F90	(revision 13313)
+++ NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_sponge.F90	(revision 13337)
@@ -55,4 +55,5 @@
       Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = .TRUE.
+      l_vremap              = ln_vremap
       tabspongedone_tsn     = .FALSE.
       !
@@ -60,4 +61,5 @@
       !
       Agrif_UseSpecialValue = .FALSE.
+      l_vremap              = .FALSE.
 #endif
       !
@@ -80,4 +82,5 @@
       Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = ln_spc_dyn
+      l_vremap              = ln_vremap
       use_sign_north        = .TRUE.
       sign_north            = -1._wp
@@ -93,4 +96,5 @@
       Agrif_UseSpecialValue = .FALSE.
       use_sign_north        = .FALSE.
+      l_vremap              = .FALSE.
 #endif
       !
@@ -109,10 +113,6 @@
       REAL(wp) ::   z1_ispongearea, z1_jspongearea
       REAL(wp), DIMENSION(jpi,jpj) :: ztabramp
-#if defined key_vertical
       REAL(wp), DIMENSION(jpi,jpj) :: ztabrampu
       REAL(wp), DIMENSION(jpi,jpj) :: ztabrampv
-#endif
-      REAL(wp), DIMENSION(jpjmax)  :: zmskwest,  zmskeast
-      REAL(wp), DIMENSION(jpimax)  :: zmsknorth, zmsksouth
       !!----------------------------------------------------------------------
       !
@@ -130,56 +130,4 @@
 #if defined SPONGE || defined SPONGE_TOP
       IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
-         !
-         ! Retrieve masks at open boundaries:
-
-         IF( lk_west ) THEN                             ! --- West --- !
-            ztabramp(:,:) = 0._wp
-            ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
-            DO ji = mi0(ind1), mi1(ind1)                
-               ztabramp(ji,:) = ssumask(ji,:)
-            END DO
-            zmskwest(    1:jpj   ) = MAXVAL(ztabramp(:,:), dim=1)
-            zmskwest(jpj+1:jpjmax) = 0._wp
-         ENDIF
-         IF( lk_east ) THEN                             ! --- East --- !
-            ztabramp(:,:) = 0._wp
-            ind1 = jpiglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells
-            DO ji = mi0(ind1), mi1(ind1)                 
-               ztabramp(ji,:) = ssumask(ji,:)
-            END DO
-            zmskeast(    1:jpj   ) = MAXVAL(ztabramp(:,:), dim=1)
-            zmskeast(jpj+1:jpjmax) = 0._wp
-         ENDIF
-         IF( lk_south ) THEN                            ! --- South --- !
-            ztabramp(:,:) = 0._wp
-            ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
-            DO jj = mj0(ind1), mj1(ind1)                 
-               ztabramp(:,jj) = ssvmask(:,jj)
-            END DO
-            zmsksouth(    1:jpi   ) = MAXVAL(ztabramp(:,:), dim=2)
-            zmsksouth(jpi+1:jpimax) = 0._wp
-         ENDIF
-         IF( lk_north ) THEN                            ! --- North --- !
-            ztabramp(:,:) = 0._wp
-            ind1 = jpjglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells
-            DO jj = mj0(ind1), mj1(ind1)                 
-               ztabramp(:,jj) = ssvmask(:,jj)
-            END DO
-            zmsknorth(    1:jpi   ) = MAXVAL(ztabramp(:,:), dim=2)
-            zmsknorth(jpi+1:jpimax) = 0._wp
-         ENDIF
-
-         ! JC: SPONGE MASKING TO BE SORTED OUT:
-         zmskwest(:)  = 1._wp
-         zmskeast(:)  = 1._wp
-         zmsksouth(:) = 1._wp
-         zmsknorth(:) = 1._wp
-#if defined key_mpp_mpi
-!         CALL mpp_max( 'AGRIF_sponge', zmskwest(:) , jpjmax )
-!         CALL mpp_max( 'AGRIF_Sponge', zmskeast(:) , jpjmax )
-!         CALL mpp_max( 'AGRIF_Sponge', zmsksouth(:), jpimax )
-!         CALL mpp_max( 'AGRIF_Sponge', zmsknorth(:), jpimax )
-#endif
-
          ! Define ramp from boundaries towards domain interior at T-points
          ! Store it in ztabramp
@@ -201,5 +149,5 @@
             DO ji = mi0(ind1), mi1(ind2)   
                DO jj = 1, jpj               
-                  ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea   * zmskwest(jj)
+                  ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea
                END DO
             END DO
@@ -209,5 +157,5 @@
             DO ji = mi0(ind1), mi1(ind2)   
                DO jj = 1, jpj               
-                  ztabramp(ji,jj) = zmskwest(jj)
+                  ztabramp(ji,jj) = 1._wp
                END DO
             END DO
@@ -218,5 +166,5 @@
             DO ji = mi0(ind1), mi1(ind2)
                DO jj = 1, jpj
-                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea ) * zmskeast(jj)
+                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea ) 
                END DO
             END DO
@@ -226,5 +174,5 @@
             DO ji = mi0(ind1), mi1(ind2)
                DO jj = 1, jpj
-                  ztabramp(ji,jj) = zmskeast(jj)
+                  ztabramp(ji,jj) = 1._wp
                END DO
             END DO
@@ -235,5 +183,5 @@
             DO jj = mj0(ind1), mj1(ind2) 
                DO ji = 1, jpi
-                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea ) * zmsksouth(ji)
+                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea )
                END DO
             END DO
@@ -243,5 +191,5 @@
             DO jj = mj0(ind1), mj1(ind2) 
                DO ji = 1, jpi
-                  ztabramp(ji,jj) = zmsksouth(ji)
+                  ztabramp(ji,jj) = 1._wp
                END DO
             END DO
@@ -252,5 +200,5 @@
             DO jj = mj0(ind1), mj1(ind2)
                DO ji = 1, jpi
-                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea ) * zmsknorth(ji)
+                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea ) 
                END DO
             END DO
@@ -260,5 +208,5 @@
             DO jj = mj0(ind1), mj1(ind2)
                DO ji = 1, jpi
-                  ztabramp(ji,jj) = zmsknorth(ji)
+                  ztabramp(ji,jj) = 1._wp
                END DO
             END DO
@@ -303,5 +251,4 @@
       ENDIF
 
-#if defined key_vertical
       ! Remove vertical interpolation where not needed:
       DO_2D( 0, 0, 0, 0 )
@@ -327,5 +274,4 @@
       mbku_parent(:,:) = NINT( ztabrampu(:,:) )
       mbkv_parent(:,:) = NINT( ztabrampv(:,:) )
-#endif
       !
 #endif
@@ -366,85 +312,92 @@
          END DO
 
-# if defined key_vertical
-        ! Interpolate thicknesses
-        ! Warning: these are masked, hence extrapolated prior interpolation.
-        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)
+        IF ( l_vremap ) THEN
+
+           ! Interpolate thicknesses
+           ! Warning: these are masked, hence extrapolated prior interpolation.
+           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
-        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)
-                  tabres(ji,jj,jk,jpts+1) = 0._wp
-              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)
+                    tabres(ji,jj,jk,jpts+1) = 0._wp
+                 END DO
+              END DO           
+           END IF
+     
+           ! Save ssh at last level:
+           IF (.NOT.ln_linssh) THEN
+              tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1) 
+           ELSE
+              tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp
+           END IF      
         END IF
-     
-        ! Save ssh at last level:
-        IF (.NOT.ln_linssh) THEN
-           tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1) 
-        ELSE
-           tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp
-        END IF      
-# endif
 
       ELSE   
          !
-# if defined key_vertical
-
-         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 
-               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) + tabres(ji,jj,k2,n2) - zhtot
+         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 
+                  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) + tabres(ji,jj,k2,n2) - zhtot
+                     ELSE
+                        h_in(jk) = tabres(ji,jj,jk,n2)
+                     END IF
+                     zhtot = zhtot + h_in(jk)
+                     tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)
+                  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
+                  END DO
+
+                  ! Account for small differences in the free-surface
+                  IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
+                     h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
                   ELSE
-                     h_in(jk) = tabres(ji,jj,jk,n2)
+                     h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
+                  END IF
+                  IF (N_in*N_out > 0) THEN
+                     CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
                   ENDIF
-                  zhtot = zhtot + h_in(jk)
-                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)
-               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
-               END DO
-
-               ! Account for small differences in free-surface
-               IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
-                  h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
-               ELSE
-                  h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
-               ENDIF
-               IF (N_in*N_out > 0) THEN
-                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
-               ENDIF
-            END DO
-         END DO
-# endif
-
-         DO jj=j1,j2
-            DO ji=i1,i2
-               DO jk=1,jpkm1
-# if defined key_vertical
-                  tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk)
-# else
-                  tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk)
-# endif
-               END DO
-            END DO
-         END DO
+               END DO
+            END DO
+
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  DO jk=1,jpkm1
+                     tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+
+         ELSE
+
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  DO jk=1,jpkm1
+                     tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         END IF
 
          DO jn = 1, jpts            
@@ -528,77 +481,87 @@
                DO ji=i1,i2
                   tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a)
-# if defined key_vertical
-                  tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk)
-# endif
-               END DO
-            END DO
-         END DO
-
-# if defined key_vertical
-         ! Extrapolate thicknesses in partial bottom cells:
-         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
-         IF (ln_zps) THEN
+               END DO
+            END DO
+         END DO
+
+         IF ( l_vremap ) THEN
+
+            DO jk=k1,k2
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk)
+                  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 = mbku(ji,jj)
+                     tabres(ji,jj,jk,m2) = 0._wp
+                  END DO
+               END DO           
+            END IF
+            ! Save ssh at last level:
+            tabres(i1:i2,j1:j2,k2,m2) = 0._wp
+            IF (.NOT.ln_linssh) THEN
+               ! This vertical sum below should be replaced by the sea-level at U-points (optimization):
+               DO jk=1,jpk
+                  tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk)
+               END DO
+               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2)
+            END IF 
+         END IF
+
+      ELSE
+
+         IF ( l_vremap ) THEN
+
+            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
+
             DO jj=j1,j2
                DO ji=i1,i2
-                  jk = mbku(ji,jj)
-                  tabres(ji,jj,jk,m2) = 0._wp
-               END DO
-            END DO           
-         END IF
-        ! Save ssh at last level:
-        tabres(i1:i2,j1:j2,k2,m2) = 0._wp
-        IF (.NOT.ln_linssh) THEN
-           ! This vertical sum below should be replaced by the sea-level at U-points (optimization):
-           DO jk=1,jpk
-              tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk)
-           END DO
-           tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2)
-        END IF 
-# endif
-
-      ELSE
-
-# if defined key_vertical
-         IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
-
-         DO jj=j1,j2
-            DO ji=i1,i2
-               tabres_child(ji,jj,:) = 0._wp
-               N_in = mbku_parent(ji,jj)
-               zhtot = 0._wp
-               DO jk=1,N_in
-                  IF (jk==N_in) THEN
-                     h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
+                  tabres_child(ji,jj,:) = 0._wp
+                  N_in = mbku_parent(ji,jj)
+                  zhtot = 0._wp
+                  DO jk=1,N_in
+                     IF (jk==N_in) THEN
+                        h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
+                     ELSE
+                        h_in(jk) = tabres(ji,jj,jk,m2)
+                     ENDIF
+                     zhtot = zhtot + h_in(jk)
+                     tabin(jk) = tabres(ji,jj,jk,m1)
+                  END DO
+                  !         
+                  N_out = 0
+                  DO jk=1,jpk
+                     IF (umask(ji,jj,jk) == 0) EXIT
+                     N_out = N_out + 1
+                     h_out(N_out) = e3u(ji,jj,jk,Kbb_a)
+                  END DO
+
+                  ! Account for small differences in free-surface
+                  IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
+                     h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
                   ELSE
-                     h_in(jk) = tabres(ji,jj,jk,m2)
+                     h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
                   ENDIF
-                  zhtot = zhtot + h_in(jk)
-                  tabin(jk) = tabres(ji,jj,jk,m1)
-               END DO
-               !         
-               N_out = 0
-               DO jk=1,jpk
-                  IF (umask(ji,jj,jk) == 0) EXIT
-                  N_out = N_out + 1
-                  h_out(N_out) = e3u(ji,jj,jk,Kbb_a)
-               END DO
-
-               ! Account for small differences in free-surface
-               IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
-                  h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
-               ELSE
-                  h_in(1)   = h_in(1)   - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
-               ENDIF
                   
-               IF (N_in * N_out > 0) THEN
-                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
-               ENDIF 
-            END DO
-         END DO
-
-         ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
-#else
-         ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
-#endif
+                  IF (N_in * N_out > 0) THEN
+                     CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
+                  ENDIF 
+               END DO
+            END DO
+
+            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
+
+         ELSE
+
+            ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
+
+         ENDIF
          !
          DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -705,76 +668,85 @@
                DO ji=i1,i2
                   tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a)
-# if defined key_vertical
-                  tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a)
-# endif
-               END DO
-            END DO
-         END DO
-
-# if defined key_vertical
-         ! Extrapolate thicknesses in partial bottom cells:
-         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
-         IF (ln_zps) THEN
+               END DO
+            END DO
+         END DO
+
+         IF ( l_vremap ) THEN
+
+            DO jk=k1,k2
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_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 = mbkv(ji,jj)
+                     tabres(ji,jj,jk,m2) = 0._wp
+                  END DO
+               END DO           
+            END IF
+            ! Save ssh at last level:
+            tabres(i1:i2,j1:j2,k2,m2) = 0._wp
+            IF (.NOT.ln_linssh) THEN
+               ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
+               DO jk=1,jpk
+                  tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v(i1:i2,j1:j2,jk,Kbb_a) * vmask(i1:i2,j1:j2,jk)
+               END DO
+               tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2)
+            END IF
+
+         END IF 
+
+      ELSE
+
+         IF ( l_vremap ) THEN
+            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
             DO jj=j1,j2
                DO ji=i1,i2
-                  jk = mbkv(ji,jj)
-                  tabres(ji,jj,jk,m2) = 0._wp
-               END DO
-            END DO           
-         END IF
-        ! Save ssh at last level:
-        tabres(i1:i2,j1:j2,k2,m2) = 0._wp
-        IF (.NOT.ln_linssh) THEN
-           ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
-           DO jk=1,jpk
-              tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v(i1:i2,j1:j2,jk,Kbb_a) * vmask(i1:i2,j1:j2,jk)
-           END DO
-           tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2)
-        END IF 
-# endif
-
-      ELSE
-
-# if defined key_vertical
-         IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
-         DO jj=j1,j2
-            DO ji=i1,i2
-               tabres_child(ji,jj,:) = 0._wp
-               N_in = mbkv_parent(ji,jj)
-               zhtot = 0._wp
-               DO jk=1,N_in
-                  IF (jk==N_in) THEN
-                     h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
+                  tabres_child(ji,jj,:) = 0._wp
+                  N_in = mbkv_parent(ji,jj)
+                  zhtot = 0._wp
+                  DO jk=1,N_in
+                     IF (jk==N_in) THEN
+                        h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
+                     ELSE
+                        h_in(jk) = tabres(ji,jj,jk,m2)
+                     ENDIF
+                     zhtot = zhtot + h_in(jk)
+                     tabin(jk) = tabres(ji,jj,jk,m1)
+                  END DO
+                  !          
+                  N_out = 0
+                  DO jk=1,jpk
+                     IF (vmask(ji,jj,jk) == 0) EXIT
+                     N_out = N_out + 1
+                     h_out(N_out) = e3v(ji,jj,jk,Kbb_a)
+                  END DO
+
+                  ! Account for small differences in free-surface
+                  IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
+                     h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
                   ELSE
-                     h_in(jk) = tabres(ji,jj,jk,m2)
+                     h_in(1)   = h_in(1) - (  sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
                   ENDIF
-                  zhtot = zhtot + h_in(jk)
-                  tabin(jk) = tabres(ji,jj,jk,m1)
-               END DO
-               !          
-               N_out = 0
-               DO jk=1,jpk
-                  IF (vmask(ji,jj,jk) == 0) EXIT
-                  N_out = N_out + 1
-                  h_out(N_out) = e3v(ji,jj,jk,Kbb_a)
-               END DO
-
-               ! Account for small differences in free-surface
-               IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
-                  h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
-               ELSE
-                  h_in(1)   = h_in(1) - (  sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
-               ENDIF
          
-               IF (N_in * N_out > 0) THEN
-                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
-               ENDIF
-            END DO
-         END DO
-
-         vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)  
-# else
-         vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
-# endif
+                  IF (N_in * N_out > 0) THEN
+                     CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
+                  ENDIF
+               END DO
+            END DO
+
+            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)  
+
+         ELSE
+
+            vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
+
+         ENDIF
          !
          DO jk = 1, jpkm1                                 ! Horizontal slab
