Index: /NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_oce_interp.F90
===================================================================
--- /NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_oce_interp.F90	(revision 13229)
+++ /NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_oce_interp.F90	(revision 13230)
@@ -34,4 +34,5 @@
    USE lib_mpp
    USE vremap
+   USE lbclnk
  
    IMPLICIT NONE
@@ -44,7 +45,7 @@
    PUBLIC   interpunb, interpvnb , interpub2b, interpvb2b
    PUBLIC   interpe3t, interpglamt, interpgphit
-#if defined key_vertical
    PUBLIC   interpht0, interpmbkt
-# endif
+   PUBLIC   agrif_initts, agrif_initssh
+
    INTEGER ::   bdy_tinterp = 0
 
@@ -89,6 +90,9 @@
       Agrif_UseSpecialValue = ln_spc_dyn
       !
+      use_sign_north = .TRUE.
+      sign_north = -1.
       CALL Agrif_Bc_variable( un_interp_id, procname=interpun )
       CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn )
+      use_sign_north = .FALSE.
       !
       Agrif_UseSpecialValue = .FALSE.
@@ -99,5 +103,5 @@
          ibdy2 = nn_hls + 1 + nbghostcells   ! halo + land + nbghostcells
          !
-         IF( .NOT.ln_dynspg_ts ) THEN  ! Store tangential transport
+         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
             DO ji = mi0(ibdy1), mi1(ibdy2)
                uu_b(ji,:,Krhs_a) = 0._wp
@@ -116,5 +120,5 @@
          !
          DO ji = mi0(ibdy1), mi1(ibdy2)
-            zub(ji,:) = 0._wp    ! Correct tangential transport
+            zub(ji,:) = 0._wp    ! Correct transport
             DO jk = 1, jpkm1
                DO jj = 1, jpj
@@ -126,5 +130,5 @@
                zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
             END DO
-            
+               
             DO jk = 1, jpkm1
                DO jj = 1, jpj
@@ -133,5 +137,5 @@
             END DO
          END DO
-            
+               
          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
             DO ji = mi0(ibdy1), mi1(ibdy2)
@@ -155,5 +159,5 @@
 
       ! --- East --- !
-      IF( lk_east ) THEN
+      IF( lk_east) THEN
          ibdy1 = jpiglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells
          ibdy2 = jpiglo - ( nn_hls + 2 )                 ! halo + land + 1
@@ -185,5 +189,5 @@
                zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
             END DO
-            
+               
             DO jk = 1, jpkm1
                DO jj = 1, jpj
@@ -193,8 +197,8 @@
             END DO
          END DO
-            
+               
          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
-            ibdy1 = jpiglo - ( nn_hls + nbghostcells )   ! halo + land + nbghostcells - 1
-            ibdy2 = jpiglo - ( nn_hls + 1 )              ! halo + land + 1            - 1
+            ibdy1 = jpiglo-nbghostcells
+            ibdy2 = jpiglo-1 
             DO ji = mi0(ibdy1), mi1(ibdy2)
                zvb(ji,:) = 0._wp
@@ -257,5 +261,5 @@
             END DO
          END DO
-            
+               
          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
             DO jj = mj0(jbdy1), mj1(jbdy2)
@@ -270,9 +274,9 @@
                   zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
                END DO
-               
+                  
                DO jk = 1, jpkm1
                   DO ji = 1, jpi
                      uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 
-                        & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk)
+                       & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk)
                   END DO
                END DO
@@ -316,12 +320,12 @@
                DO ji = 1, jpi
                   vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 
-                     & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-            
+                    & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+               
          IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
-            jbdy1 = jpjglo - ( nn_hls + nbghostcells )   ! halo + land + nbghostcells - 1
-            jbdy2 = jpjglo - ( nn_hls + 1 )              ! halo + land + 1            - 1
+            jbdy1 = jpjglo-nbghostcells
+            jbdy2 = jpjglo-1
             DO jj = mj0(jbdy1), mj1(jbdy2)
                zub(:,jj) = 0._wp
@@ -335,5 +339,5 @@
                   zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a)
                END DO
-               
+                  
                DO jk = 1, jpkm1
                   DO ji = 1, jpi
@@ -378,4 +382,5 @@
          iend   = jpiglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1
          DO ji = mi0(istart), mi1(iend)
+
             DO jj=1,jpj
                va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
@@ -389,5 +394,5 @@
             END DO
          END DO
-      ENDIF
+      ENDIF 
       !
       !--- South ---!
@@ -395,5 +400,8 @@
          jstart = nn_hls + 2                              ! halo + land + 1
          jend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
+         jstart = 2
+         jend   = nbghostcells+1
          DO jj = mj0(jstart), mj1(jend)
+
             DO ji=1,jpi
                ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
@@ -401,5 +409,5 @@
             END DO
          END DO
-      ENDIF
+      ENDIF       
       !
       !--- North ---!
@@ -419,8 +427,9 @@
             END DO
          END DO
-      ENDIF
+      ENDIF 
       !
    END SUBROUTINE Agrif_dyn_ts
 
+   
    SUBROUTINE Agrif_dyn_ts_flux( jn, zu, zv )
       !!----------------------------------------------------------------------
@@ -498,4 +507,5 @@
    END SUBROUTINE Agrif_dyn_ts_flux
 
+   
    SUBROUTINE Agrif_dta_ts( kt )
       !!----------------------------------------------------------------------
@@ -513,9 +523,13 @@
       !
       ! Enforce volume conservation if no time refinement:  
-      IF    ( Agrif_rhot()==1 ) ll_int_cons=.TRUE.  
+      IF ( Agrif_rhot()==1 ) ll_int_cons=.TRUE.  
       !
       ! Interpolate barotropic fluxes
       Agrif_SpecialValue = 0._wp
       Agrif_UseSpecialValue = ln_spc_dyn
+
+      use_sign_north = .TRUE.
+      sign_north = -1.
+
       !
       ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners)
@@ -542,4 +556,5 @@
       ENDIF
       Agrif_UseSpecialValue = .FALSE.
+      use_sign_north = .FALSE.
       ! 
    END SUBROUTINE Agrif_dta_ts
@@ -566,16 +581,16 @@
       !
       ! --- West --- !
-      IF( lk_west ) THEN
-        istart = nn_hls + 2                              ! halo + land + 1
+      IF(lk_west) THEN
+         istart = nn_hls + 2                              ! halo + land + 1
          iend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
          DO ji = mi0(istart), mi1(iend)
             DO jj = 1, jpj
                ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
       ENDIF
       !
       ! --- East --- !
-      IF( lk_east ) THEN
+      IF(lk_east) THEN
          istart = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
          iend   = jpiglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1
@@ -583,10 +598,10 @@
             DO jj = 1, jpj
                ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
       ENDIF
       !
       ! --- South --- !
-      IF( lk_south ) THEN
+      IF(lk_south) THEN
          jstart = nn_hls + 2                              ! halo + land + 1
          jend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
@@ -594,10 +609,10 @@
             DO ji = 1, jpi
                ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
       ENDIF
       !
       ! --- North --- !
-      IF( lk_north ) THEN
+      IF(lk_north) THEN
          jstart = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
          jend   = jpjglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1
@@ -605,6 +620,6 @@
             DO ji = 1, jpi
                ssh(ji,jj,Krhs_a) = hbdy(ji,jj)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
       ENDIF
       !
@@ -625,5 +640,5 @@
       !
       ! --- West --- !
-      IF( lk_west ) THEN
+      IF(lk_west) THEN
          istart = nn_hls + 2                              ! halo + land + 1
          iend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
@@ -631,10 +646,10 @@
             DO jj = 1, jpj
                ssha_e(ji,jj) = hbdy(ji,jj)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
       ENDIF
       !
       ! --- East --- !
-      IF( lk_east ) THEN
+      IF(lk_east) THEN
          istart = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
          iend   = jpiglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1
@@ -642,10 +657,10 @@
             DO jj = 1, jpj
                ssha_e(ji,jj) = hbdy(ji,jj)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
       ENDIF
       !
       ! --- South --- !
-      IF( lk_south ) THEN
+      IF(lk_south) THEN
          jstart = nn_hls + 2                              ! halo + land + 1
          jend   = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
@@ -653,10 +668,10 @@
             DO ji = 1, jpi
                ssha_e(ji,jj) = hbdy(ji,jj)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
       ENDIF
       !
       ! --- North --- !
-      IF( lk_north ) THEN
+      IF(lk_north) THEN
          jstart = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
          jend   = jpjglo - ( nn_hls + 1 )                 ! halo + land + 1            - 1
@@ -664,10 +679,11 @@
             DO ji = 1, jpi
                ssha_e(ji,jj) = hbdy(ji,jj)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
       ENDIF
       !
    END SUBROUTINE Agrif_ssh_ts
 
+   
    SUBROUTINE Agrif_avm
       !!----------------------------------------------------------------------
@@ -690,5 +706,5 @@
       !
    END SUBROUTINE Agrif_avm
-   
+
 
    SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
@@ -702,12 +718,17 @@
       INTEGER  ::   ji, jj, jk, jn  ! dummy loop indices
       INTEGER  ::   N_in, N_out
+      INTEGER  :: item
       ! vertical interpolation:
       REAL(wp) :: zhtot
       REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin
-      REAL(wp), DIMENSION(k1:k2) :: h_in
-      REAL(wp), DIMENSION(1:jpk) :: h_out
-      !!----------------------------------------------------------------------
-
-      IF( before ) THEN         
+      REAL(wp), DIMENSION(k1:k2) :: h_in, z_in
+      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
+      !!----------------------------------------------------------------------
+
+      IF( before ) THEN
+
+         item = Kmm_a
+         IF( l_ini_child )   Kmm_a = Kbb_a  
+
          DO jn = 1,jpts
             DO jk=k1,k2
@@ -718,71 +739,96 @@
               END DO
            END DO
-        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
-                  ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * 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,jpts+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,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
-        ELSE
-           ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp
-        END IF      
-# endif
+         END DO
+
+         IF( l_vremap .OR. l_ini_child) THEN
+            ! Interpolate thicknesses
+            ! Warning: these are masked, hence extrapolated prior interpolation.
+            DO jk=k1,k2
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                      ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * 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,jpts+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,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
+            ELSE
+               ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp
+            END IF      
+         ENDIF
+         Kmm_a = item
+
       ELSE 
-
-# if defined key_vertical 
-         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 
-            
-         DO jj=j1,j2
-            DO ji=i1,i2
-               ts(ji,jj,:,:,Krhs_a) = 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) + ptab(ji,jj,k2,n2) - zhtot
-                  ELSE
-                     h_in(jk) = ptab(ji,jj,jk,n2)
+         item = Krhs_a
+         IF( l_ini_child )   Krhs_a = Kbb_a  
+
+         IF( l_vremap .OR. l_ini_child ) THEN
+            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 
+               
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  ts(ji,jj,:,:,Krhs_a) = 0.                  
+               !   IF( l_ini_child) ts(ji,jj,:,:,Krhs_a) = ptab(ji,jj,:,1:jpts)
+                  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)
+                  DO jk=2,N_in
+                     z_in(jk) = z_in(jk-1) + 0.5_wp * 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
+                     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)
+                  DO jk=2,N_out
+                     z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk)
+                  END DO
+
+                  IF (N_in*N_out > 0) THEN
+                     IF( l_ini_child ) THEN
+                        CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),          &
+                                      &   z_out(1:N_out),N_in,N_out,jpts)  
+                     ELSE 
+                        CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),   &
+                                      &   h_out(1:N_out),N_in,N_out,jpts)  
+                     ENDIF
                   ENDIF
-                  zhtot = zhtot + h_in(jk)
-                  tabin(jk,:) = ptab(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._wp) EXIT 
-                  N_out = N_out + 1
-                  h_out(jk) = e3t(ji,jj,jk,Krhs_a)
-               ENDDO
-               IF (N_in*N_out > 0) THEN
-                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),h_out(1:N_out),N_in,N_out,jpts)
-               ENDIF
-            ENDDO
-         ENDDO
-# else
-         !
-         DO jn=1, jpts
-            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
+               END DO
+            END DO
+            Krhs_a = item
+ 
+         ELSE
+         
+            DO jn=1, jpts
+                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
 
       ENDIF
@@ -790,4 +836,5 @@
    END SUBROUTINE interptsn
 
+   
    SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before )
       !!----------------------------------------------------------------------
@@ -808,4 +855,5 @@
    END SUBROUTINE interpsshn
 
+   
    SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
       !!----------------------------------------------------------------------
@@ -820,85 +868,111 @@
       REAL(wp) :: zrhoy, zhtot
       ! vertical interpolation:
-      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
-      REAL(wp), DIMENSION(1:jpk) :: h_out
-      INTEGER  :: N_in, N_out
+      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in
+      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
+      INTEGER  :: N_in, N_out,item
       REAL(wp) :: h_diff
       !!---------------------------------------------    
       !
       IF (before) THEN 
+
+         item = Kmm_a
+         IF( l_ini_child )   Kmm_a = Kbb_a     
+
          DO jk=1,jpk
             DO jj=j1,j2
                DO ji=i1,i2
                   ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)*umask(ji,jj,jk)) 
-# if defined key_vertical
-                  ! Interpolate thicknesses (masked for subsequent extrapolation)
-                  ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
-# endif
-               END DO
-            END DO
-         END DO
-# if defined key_vertical
+                  IF( l_vremap .OR. l_ini_child) THEN
+                     ! Interpolate thicknesses (masked for subsequent extrapolation)
+                     ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
+                  ENDIF
+               END DO
+            END DO
+         END DO
+
+        IF( l_vremap .OR. l_ini_child) THEN
          ! 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)
-                  ptab(ji,jj,jk,2) = 0._wp
-               END DO
-            END DO           
-         END IF
-        ! Save ssh at last level:
-        ptab(i1:i2,j1:j2,k2,2) = 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
-              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u(i1:i2,j1:j2,jk,Kmm_a) * umask(i1:i2,j1:j2,jk)
-           END DO
-           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2)
-        END IF 
-# endif
+            IF (ln_zps) THEN
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     jk = mbku(ji,jj)
+                     ptab(ji,jj,jk,2) = 0._wp
+                  END DO
+               END DO           
+            END IF
+
+           ! Save ssh at last level:
+           ptab(i1:i2,j1:j2,k2,2) = 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
+                 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u(i1:i2,j1:j2,jk,Kmm_a) * umask(i1:i2,j1:j2,jk)
+              END DO
+              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2)
+           END IF
+        ENDIF
+
+         Kmm_a = item
          !
       ELSE
          zrhoy = Agrif_rhoy()
-# if defined key_vertical
+
+        IF( l_vremap .OR. l_ini_child) THEN
 ! VERTICAL REFINEMENT BEGIN
 
-         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
-
-         DO ji=i1,i2
-            DO jj=j1,j2
-               uu(ji,jj,:,Krhs_a) = 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) + ptab(ji,jj,k2,2) - zhtot
-                  ELSE
-                     h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
-                  ENDIF
-                  zhtot = zhtot + h_in(jk)
-                  tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk))
-              ENDDO
-                  
-              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,Krhs_a)
-              ENDDO
-              IF (N_in*N_out > 0) THEN
-                 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1)
-              ENDIF
-            ENDDO
-         ENDDO
-
-# else
-         DO jk = 1, jpkm1
-            DO jj=j1,j2
-               uu(i1:i2,jj,jk,Krhs_a) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs_a) )
-            END DO
-         END DO
-# endif
+            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
+
+            DO ji=i1,i2
+               DO jj=j1,j2
+                  uu(ji,jj,:,Krhs_a) = 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) + ptab(ji,jj,k2,2) - zhtot
+                     ELSE
+                        h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
+                     ENDIF
+                     zhtot = zhtot + h_in(jk)
+                     IF( h_in(jk) .GT. 0. ) THEN
+                     tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk))
+                     ELSE
+                     tabin(jk) = 0.
+                     ENDIF
+                 END DO
+                 z_in(1) = 0.5_wp * h_in(1) - zhtot + hu0_parent(ji,jj) 
+                 DO jk=2,N_in
+                    z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk)
+                 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,Krhs_a)
+                 END DO
+
+                 z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hu_0(ji,jj)
+                 DO jk=2,N_out
+                    z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk) 
+                 END DO  
+
+                 IF (N_in*N_out > 0) THEN
+                     IF( l_ini_child ) THEN
+                        CALL remap_linear       (tabin(1:N_in),z_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),z_out(1:N_out),N_in,N_out,1)
+                     ELSE
+                        CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1)
+                     ENDIF   
+                 ENDIF
+               END DO
+            END DO
+         ELSE
+            DO jk = 1, jpkm1
+               DO jj=j1,j2
+                  uu(i1:i2,jj,jk,Krhs_a) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs_a) )
+               END DO
+            END DO
+         ENDIF
 
       ENDIF
@@ -906,4 +980,5 @@
    END SUBROUTINE interpun
 
+   
    SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
       !!----------------------------------------------------------------------
@@ -918,80 +993,107 @@
       REAL(wp) :: zrhox
       ! vertical interpolation:
-      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
-      REAL(wp), DIMENSION(1:jpk) :: h_out
-      INTEGER  :: N_in, N_out
+      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in
+      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
+      INTEGER  :: N_in, N_out, item
       REAL(wp) :: h_diff, zhtot
       !!---------------------------------------------    
       !      
-      IF (before) THEN          
+      IF (before) THEN   
+
+         item = Kmm_a
+         IF( l_ini_child )   Kmm_a = Kbb_a     
+       
          DO jk=k1,k2
             DO jj=j1,j2
                DO ji=i1,i2
                   ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)*vmask(ji,jj,jk))
-# if defined key_vertical
-                  ! Interpolate thicknesses (masked for subsequent extrapolation)
-                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
-# endif
-               END DO
-            END DO
-         END DO
-# if defined key_vertical
+                  IF( l_vremap .OR. l_ini_child) THEN
+                     ! Interpolate thicknesses (masked for subsequent extrapolation)
+                     ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
+                  ENDIF
+               END DO
+            END DO
+         END DO
+
+         IF( l_vremap .OR. l_ini_child) THEN
          ! Extrapolate thicknesses in partial bottom cells:
          ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
-         IF (ln_zps) THEN
+            IF (ln_zps) THEN
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     jk = mbkv(ji,jj)
+                     ptab(ji,jj,jk,2) = 0._wp
+                  END DO
+               END DO           
+            END IF
+            ! Save ssh at last level:
+            ptab(i1:i2,j1:j2,k2,2) = 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
+                  ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk)
+               END DO
+               ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2)
+            END IF 
+         ENDIF
+         item = Kmm_a
+
+      ELSE       
+         zrhox = Agrif_rhox()
+
+         IF( l_vremap .OR. l_ini_child ) THEN
+
+            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
+
             DO jj=j1,j2
                DO ji=i1,i2
-                  jk = mbkv(ji,jj)
-                  ptab(ji,jj,jk,2) = 0._wp
-               END DO
-            END DO           
-         END IF
-        ! Save ssh at last level:
-        ptab(i1:i2,j1:j2,k2,2) = 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
-              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk)
-           END DO
-           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2)
-        END IF 
-# endif
-      ELSE       
-         zrhox = Agrif_rhox()
-# if defined key_vertical
-
-         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
-
-         DO jj=j1,j2
-            DO ji=i1,i2
-               vv(ji,jj,:,Krhs_a) = 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) + ptab(ji,jj,k2,2) - zhtot
-                  ELSE
-                     h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 
+                  vv(ji,jj,:,Krhs_a) = 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) + ptab(ji,jj,k2,2) - zhtot
+                     ELSE
+                        h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 
+                     ENDIF
+                     zhtot = zhtot + h_in(jk)
+                     IF( h_in(jk) .GT. 0. ) THEN
+                       tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk))
+                     ELSE
+                       tabin(jk)  = 0.
+                     ENDIF 
+                  END DO
+
+                  z_in(1) = 0.5_wp * h_in(1) - zhtot + hv0_parent(ji,jj)
+                  DO jk=2,N_in
+                     z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk)
+                  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,Krhs_a)
+                  END DO
+
+                  z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hv_0(ji,jj)
+                  DO jk=2,N_out
+                     z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk)
+                  END DO
+ 
+                  IF (N_in*N_out > 0) THEN
+                     IF( l_ini_child ) THEN
+                        CALL remap_linear       (tabin(1:N_in),z_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),z_out(1:N_out),N_in,N_out,1)
+                     ELSE
+                        CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1)
+                     ENDIF   
                   ENDIF
-                  zhtot = zhtot + h_in(jk)
-                  tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk))
-              ENDDO
-         
-               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,Krhs_a)
-               END DO
-               IF (N_in*N_out > 0) THEN
-                  call reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1)
-               ENDIF
-            END DO
-         END DO
-# else
-         DO jk = 1, jpkm1
-            vv(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs_a) )
-         END DO
-# endif
+               END DO
+            END DO
+         ELSE
+            DO jk = 1, jpkm1
+               vv(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs_a) )
+            END DO
+         ENDIF
       ENDIF
       !        
@@ -1203,5 +1305,4 @@
    END SUBROUTINE interpe3t
 
-
    SUBROUTINE interpglamt( ptab, i1, i2, j1, j2, before )
       !!----------------------------------------------------------------------
@@ -1283,66 +1384,68 @@
               END DO
            END DO
-        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
-                  ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * 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,2) = 0._wp
-              END DO
-           END DO           
-        END IF
-     
-        ! Save ssh at last level:
-        IF (.NOT.ln_linssh) THEN
-           ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
-        ELSE
-           ptab(i1:i2,j1:j2,k2,2) = 0._wp
-        END IF      
-# endif
+         END DO
+
+         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
+                      ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * 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,2) = 0._wp
+                  END DO
+               END DO           
+            END IF
+        
+           ! Save ssh at last level:
+            IF (.NOT.ln_linssh) THEN
+               ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
+            ELSE
+               ptab(i1:i2,j1:j2,k2,2) = 0._wp
+            END IF      
+          ENDIF
+
       ELSE 
-#ifdef key_vertical         
-         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
-         avm_k(i1:i2,j1:j2,k1:k2) = 0._wp
-            
-         DO jj = j1, j2
-            DO ji =i1, i2
-               N_in = mbkt_parent(ji,jj)
-               IF ( tmask(ji,jj,1) == 0._wp) N_in = 0
-               z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2)
-               DO jk = N_in, 1, -1  ! Parent vertical grid               
-                     z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2)
-                    tabin(jk) = ptab(ji,jj,jk,1)
-               END DO
-               N_out = mbkt(ji,jj) 
-               DO jk = 1, N_out        ! Child vertical grid
-                  z_out(jk) = gdepw(ji,jj,jk,Kmm_a)
-               ENDDO
-               IF (N_in*N_out > 0) THEN
-                  CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1)
-               ENDIF
-            ENDDO
-         ENDDO
-#else
-         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
-#endif
+
+         IF( l_vremap ) THEN
+            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
+            avm_k(i1:i2,j1:j2,k1:k2) = 0._wp
+               
+            DO jj = j1, j2
+               DO ji =i1, i2
+                  N_in = mbkt_parent(ji,jj)
+                  IF ( tmask(ji,jj,1) == 0._wp) N_in = 0
+                  z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2)
+                  DO jk = N_in, 1, -1  ! Parent vertical grid               
+                        z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2)
+                       tabin(jk) = ptab(ji,jj,jk,1)
+                  END DO
+                  N_out = mbkt(ji,jj) 
+                  DO jk = 1, N_out        ! Child vertical grid
+                     z_out(jk) = gdepw(ji,jj,jk,Kmm_a)
+                  END DO
+                  IF (N_in*N_out > 0) THEN
+                     CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1)
+                  ENDIF
+               END DO
+            END DO
+         ELSE
+            avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
+         ENDIF
       ENDIF
       !
    END SUBROUTINE interpavm
 
-# if defined key_vertical
+   
    SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before )
       !!----------------------------------------------------------------------
@@ -1363,4 +1466,5 @@
    END SUBROUTINE interpmbkt
 
+   
    SUBROUTINE interpht0( ptab, i1, i2, j1, j2, before )
       !!----------------------------------------------------------------------
@@ -1380,6 +1484,43 @@
       !
    END SUBROUTINE interpht0
-#endif
-
+
+   
+   SUBROUTINE agrif_initts(tabres,i1,i2,j1,j2,k1,k2,m1,m2,before)
+       INTEGER :: i1, i2, j1, j2, k1, k2, m1, m2
+       REAL(wp):: tabres(i1:i2,j1:j2,k1:k2,m1:m2)
+       LOGICAL :: before
+
+       INTEGER :: jm
+
+       IF (before) THEN
+         DO jm=1,jpts
+             tabres(i1:i2,j1:j2,k1:k2,jm) = ts(i1:i2,j1:j2,k1:k2,jm,Kbb_a)
+         END DO
+       ELSE
+         DO jm=1,jpts
+             ts(i1:i2,j1:j2,k1:k2,jm,Kbb_a)=tabres(i1:i2,j1:j2,k1:k2,jm)
+         END DO
+       ENDIF
+   END SUBROUTINE agrif_initts 
+
+   
+   SUBROUTINE agrif_initssh( ptab, i1, i2, j1, j2, before )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE interpsshn  ***
+      !!----------------------------------------------------------------------  
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
+      LOGICAL                         , INTENT(in   ) ::   before
+      !
+      !!----------------------------------------------------------------------  
+      !
+      IF( before) THEN
+         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kbb_a)
+      ELSE
+         ssh(i1:i2,j1:j2,Kbb_a) = ptab(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1)
+      ENDIF
+      !
+   END SUBROUTINE agrif_initssh
+   
 #else
    !!----------------------------------------------------------------------
Index: /NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_oce_sponge.F90
===================================================================
--- /NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_oce_sponge.F90	(revision 13229)
+++ /NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_oce_sponge.F90	(revision 13230)
@@ -32,702 +32,4 @@
 
    !! * Substitutions
-#  include "do_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/NST 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE Agrif_Sponge_Tra
-      !!----------------------------------------------------------------------
-      !!                 *** ROUTINE Agrif_Sponge_Tra ***
-      !!----------------------------------------------------------------------
-      REAL(wp) ::   zcoef   ! local scalar
-      
-      !!----------------------------------------------------------------------
-      !
-#if defined SPONGE
-      !! Assume persistence:
-      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
-
-      CALL Agrif_Sponge
-      Agrif_SpecialValue    = 0._wp
-      Agrif_UseSpecialValue = .TRUE.
-      tabspongedone_tsn     = .FALSE.
-      !
-      CALL Agrif_Bc_Variable( tsn_sponge_id, calledweight=zcoef, procname=interptsn_sponge )
-      !
-      Agrif_UseSpecialValue = .FALSE.
-#endif
-      !
-      CALL iom_put( 'agrif_spu', fspu(:,:))
-      CALL iom_put( 'agrif_spv', fspv(:,:))
-      !
-   END SUBROUTINE Agrif_Sponge_Tra
-
-
-   SUBROUTINE Agrif_Sponge_dyn
-      !!----------------------------------------------------------------------
-      !!                 *** ROUTINE Agrif_Sponge_dyn ***
-      !!----------------------------------------------------------------------
-      REAL(wp) ::   zcoef   ! local scalar
-      !!----------------------------------------------------------------------
-      !
-#if defined SPONGE
-      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
-
-      Agrif_SpecialValue=0.
-      Agrif_UseSpecialValue = ln_spc_dyn
-      !
-      tabspongedone_u = .FALSE.
-      tabspongedone_v = .FALSE.         
-      CALL Agrif_Bc_Variable( un_sponge_id, calledweight=zcoef, procname=interpun_sponge )
-      !
-      tabspongedone_u = .FALSE.
-      tabspongedone_v = .FALSE.
-      CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge )
-      !
-      Agrif_UseSpecialValue = .FALSE.
-#endif
-      !
-      CALL iom_put( 'agrif_spt', fspt(:,:))
-      CALL iom_put( 'agrif_spf', fspf(:,:))
-      !
-   END SUBROUTINE Agrif_Sponge_dyn
-
-
-   SUBROUTINE Agrif_Sponge
-      !!----------------------------------------------------------------------
-      !!                 *** ROUTINE  Agrif_Sponge ***
-      !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, ind1, ind2
-      INTEGER  ::   ispongearea, jspongearea
-      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
-      !!----------------------------------------------------------------------
-      !
-      ! Sponge 1d example with:
-      !      iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2
-      !                        
-      !coarse :     U     T     U     T     U     T     U
-      !|            |           |           |           |
-      !fine :     t u t u t u t u t u t u t u t u t u t u t
-      !sponge val:0   0   0   1  5/6 4/6 3/6 2/6 1/6  0   0
-      !           |   ghost     | <-- sponge area  -- > |
-      !           |   points    |                       |
-      !                         |--> dynamical interface
-
-#if defined SPONGE || defined SPONGE_TOP
-      IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
-         !
-         ! Retrieve masks at open boundaries:
-
-         ! --- West --- !
-         IF( lk_west ) THEN
-            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
-
-         ! --- East --- !
-         IF( lk_east ) THEN
-            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
-
-         ! --- South --- !
-         IF( lk_south ) THEN
-            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
-
-         ! --- North --- !
-         IF( lk_north ) THEN
-            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
-
-         ispongearea  = nn_sponge_len * Agrif_irhox()
-         z1_ispongearea = 1._wp / REAL( ispongearea )
-         jspongearea  = nn_sponge_len * Agrif_irhoy()
-         z1_jspongearea = 1._wp / REAL( jspongearea )
-         
-         ztabramp(:,:) = 0._wp
-
-         ! Trick to remove sponge in 2DV domains:
-         IF ( nbcellsx <= 3 ) ispongearea = -1
-         IF ( nbcellsy <= 3 ) jspongearea = -1
-
-         ! --- West --- !
-         IF( lk_west ) THEN
-            ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
-            ind2 = nn_hls + 1 + nbghostcells + ispongearea 
-            DO ji = mi0(ind1), mi1(ind2)   
-               DO jj = 1, jpj               
-                  ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj)
-               END DO
-            END DO
-
-            ! ghost cells:
-            ind1 = 1
-            ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
-            DO ji = mi0(ind1), mi1(ind2)   
-               DO jj = 1, jpj               
-                  ztabramp(ji,jj) = zmskwest(jj)
-               END DO
-            END DO
-         ENDIF
-
-         ! --- East --- !
-         IF( lk_east ) THEN
-            ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea
-            ind2 = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
-            DO ji = mi0(ind1), mi1(ind2)
-               DO jj = 1, jpj
-                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj)
-               ENDDO
-            END DO
-
-            ! ghost cells:
-            ind1 = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
-            ind2 = jpiglo
-            DO ji = mi0(ind1), mi1(ind2)
-               DO jj = 1, jpj
-                  ztabramp(ji,jj) = zmskeast(jj)
-               ENDDO
-            END DO
-         ENDIF
-
-         ! --- South --- !
-         IF( lk_south ) THEN
-            ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
-            ind2 = nn_hls + 1 + nbghostcells + jspongearea 
-            DO jj = mj0(ind1), mj1(ind2) 
-               DO ji = 1, jpi
-                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji)
-               END DO
-            END DO
-
-            ! ghost cells:
-            ind1 = 1
-            ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
-            DO jj = mj0(ind1), mj1(ind2) 
-               DO ji = 1, jpi
-                  ztabramp(ji,jj) = zmsksouth(ji)
-               END DO
-            END DO
-         ENDIF
-
-         ! --- North --- !
-         IF( lk_north ) THEN
-            ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea
-            ind2 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
-            DO jj = mj0(ind1), mj1(ind2)
-               DO ji = 1, jpi
-                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mjg(jj) - ind1 ) * z1_jspongearea) * zmsknorth(ji)
-               END DO
-            END DO
-
-            ! ghost cells:
-            ind1 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
-            ind2 = jpjglo
-            DO jj = mj0(ind1), mj1(ind2)
-               DO ji = 1, jpi
-                  ztabramp(ji,jj) = zmsknorth(ji)
-               END DO
-            END DO
-         ENDIF
-
-      ENDIF
-
-      ! Tracers
-      IF( .NOT. spongedoneT ) THEN
-         fspu(:,:) = 0._wp
-         fspv(:,:) = 0._wp
-         DO_2D_00_00
-            fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji+1,jj  ) ) * ssumask(ji,jj)
-            fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji  ,jj+1) ) * ssvmask(ji,jj)
-         END_2D
-      ENDIF
-
-      ! Dynamics
-      IF( .NOT. spongedoneU ) THEN
-         fspt(:,:) = 0._wp
-         fspf(:,:) = 0._wp
-         DO_2D_00_00
-            fspt(ji,jj) = ztabramp(ji,jj) * ssmask(ji,jj)
-            fspf(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji  ,jj+1)   &
-                                  &  +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj  ) ) &
-                                  &  * ssvmask(ji,jj) * ssvmask(ji,jj+1)
-         END_2D
-      ENDIF
-      
-      IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN
-         CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1., fspv, 'V', 1., fspt, 'T', 1., fspf, 'F', 1. )
-         spongedoneT = .TRUE.
-         spongedoneU = .TRUE.
-      ENDIF
-      IF( .NOT. spongedoneT ) THEN
-         CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1., fspv, 'V', 1. )
-         spongedoneT = .TRUE.
-      ENDIF
-      IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN
-         CALL lbc_lnk_multi( 'agrif_Sponge', fspt, 'T', 1., fspf, 'F', 1. )
-         spongedoneU = .TRUE.
-      ENDIF
-
-#if defined key_vertical
-      ! Remove vertical interpolation where not needed:
-      DO_2D_00_00
-         IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp).AND. &
-         &   (fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0
-!
-         IF ((fspt(ji+1,jj)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. &
-         &   (fspf(ji,jj-1)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0
-!
-         IF ((fspt(ji,jj+1)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. &
-         &   (fspf(ji-1,jj)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0
-!
-         IF ( ssmask(ji,jj) == 0._wp) mbkt_parent(ji,jj) = 0
-         IF (ssumask(ji,jj) == 0._wp) mbku_parent(ji,jj) = 0
-         IF (ssvmask(ji,jj) == 0._wp) mbkv_parent(ji,jj) = 0
-      END_2D
-      !
-      ztabramp (:,:) = REAL( mbkt_parent (:,:), wp )
-      ztabrampu(:,:) = REAL( mbku_parentu(:,:), wp )
-      ztabrampv(:,:) = REAL( mbkv_parentv(:,:), wp )
-      CALL lbc_lnk_multi( 'Agrif_Sponge', ztabramp, 'T', 1., ztabrampu, 'U', 1., ztabrampv, 'V', 1. )
-      mbkt_parent(:,:) = NINT( ztabramp (:,:) )
-      mbku_parent(:,:) = NINT( ztabrampu(:,:) )
-      mbkv_parent(:,:) = NINT( ztabrampv(:,:) )
-#endif
-      !
-#endif
-      !
-   END SUBROUTINE Agrif_Sponge
-
-   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
-      LOGICAL                                     , INTENT(in   ) ::   before
-      !
-      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      INTEGER  ::   iku, ikv
-      REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot, ztrelax
-      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
-      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
-      ! vertical interpolation:
-      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
-      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
-      REAL(wp), DIMENSION(k1:k2) :: h_in
-      REAL(wp), DIMENSION(1:jpk) :: h_out
-      INTEGER :: N_in, N_out
-      !!----------------------------------------------------------------------
-      !
-      IF( before ) THEN
-         DO jn = 1, jpts
-            DO jk=k1,k2
-               DO jj=j1,j2
-                  DO ji=i1,i2
-                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a)
-                  END DO
-               END DO
-            END DO
-         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)
-              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           
-        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
-                  ELSE
-                     h_in(jk) = tabres(ji,jj,jk,n2)
-                  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
-               ENDDO
-
-               ! 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
-            ENDDO
-         ENDDO
-# 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
-               ENDDO
-            ENDDO
-         ENDDO
-
-         !* set relaxation time scale
-         IF( l_1st_euler .AND. lk_agrif_fstep ) THEN   ;   ztrelax =   rn_trelax_tra  / (        rn_Dt )
-         ELSE                                          ;   ztrelax =   rn_trelax_tra  / (2._wp * rn_Dt )
-         ENDIF
-
-         DO jn = 1, jpts            
-            DO jk = 1, jpkm1
-               ztu(i1:i2,j1:j2,jk) = 0._wp
-               DO jj = j1,j2
-                  DO ji = i1,i2-1
-                     zabe1 = rn_sponge_tra * fspu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
-                     ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
-                  END DO
-               END DO
-               ztv(i1:i2,j1:j2,jk) = 0._wp
-               DO ji = i1,i2
-                  DO jj = j1,j2-1
-                     zabe2 = rn_sponge_tra * fspv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
-                     ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
-                  END DO
-               END DO
-               !
-               IF( ln_zps ) THEN      ! set gradient at partial step level
-                  DO jj = j1,j2
-                     DO ji = i1,i2
-                        ! last level
-                        iku = mbku(ji,jj)
-                        ikv = mbkv(ji,jj)
-                        IF( iku == jk )   ztu(ji,jj,jk) = 0._wp
-                        IF( ikv == jk )   ztv(ji,jj,jk) = 0._wp
-                     END DO
-                  END DO
-               ENDIF
-            END DO
-            !
-            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) / e3t(ji,jj,jk,Kmm_a)
-                        ! horizontal diffusive trends
-                        ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) &
-                             &  - ztrelax * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn) 
-                        ! add it to the general tracer trends
-                        ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa
-                     ENDIF
-                  END DO
-               END DO
-            END DO
-            !
-         END DO
-         !
-         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,m1,m2, before)
-      !!---------------------------------------------
-      !!   *** ROUTINE interpun_sponge ***
-      !!---------------------------------------------    
-      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
-      LOGICAL, INTENT(in) :: before
-
-      INTEGER :: ji,jj,jk,jmax
-      INTEGER :: ind1
-      ! sponge parameters 
-      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax
-      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
-      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
-      ! vertical interpolation:
-      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
-      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
-      REAL(wp), DIMENSION(1:jpk) :: h_out
-      INTEGER ::N_in, N_out
-      !!---------------------------------------------    
-      !
-      IF( before ) THEN
-         DO jk=k1,k2
-            DO jj=j1,j2
-               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
-            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
-                  ELSE
-                     h_in(jk) = tabres(ji,jj,jk,m2)
-                  ENDIF
-                  zhtot = zhtot + h_in(jk)
-                  tabin(jk) = tabres(ji,jj,jk,m1)
-               ENDDO
-               !         
-               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)
-               ENDDO
-
-               ! 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 
-            ENDDO
-         ENDDO
-
-         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
-         !* set relaxation time scale
-         IF( l_1st_euler .AND. lk_agrif_fstep ) THEN   ;   ztrelax =   rn_trelax_dyn  / (        rn_Dt )
-         ELSE                                          ;   ztrelax =   rn_trelax_dyn  / (2._wp * rn_Dt )
-         ENDIF
-         !
-         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) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj)
-                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kbb_a) * ubdiff(ji  ,jj,jk) &
-                                     &   -e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb_a) * 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) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(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
-         END DO
-         !
-         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) * e3u(ji,jj,jk,Kmm_a) )   &
-                         & + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) & 
-                         & - ztrelax  * fspu(ji,jj) * ubdiff(ji,jj,jk)
-
-                     ! add it to the general momentum trends
-                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua                                 
-                  END DO
-               ENDIF
-
-            END DO
-         END DO
-
-         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
-
-         jmax = j2-1
-         ind1 = jpjglo - ( nn_hls + nbghostcells + 2 )   ! North
-         DO jj = mj0(ind1), mj1(ind1)                 
-            jmax = MIN(jmax,jj)
-         END DO
-
-         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) * e3v(ji,jj,jk,Kmm_a) )   &
-                           + ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
-
-                     ! add it to the general momentum trends
-                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) + 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,m1,m2, before,nb,ndir)
-      !!---------------------------------------------
-      !!   *** ROUTINE interpvn_sponge ***
-      !!--------------------------------------------- 
-      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
-      LOGICAL, INTENT(in) :: before
-      INTEGER, INTENT(in) :: nb , ndir
-      !
-      INTEGER  ::   ji, jj, jk, imax
-      INTEGER  ::   ind1
-      ! sponge parameters 
-      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax
-      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
-      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
-      ! vertical interpolation:
-      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
-      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
-      REAL(wp), DIMENSION(1:jpk) :: h_out
-      INTEGER :: N_in, N_out
-      !!--------------------------------------------- 
-      
-      IF( before ) THEN 
-         DO jk=k1,k2
-            DO jj=j1,j2
                DO ji=i1,i2
                   tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a)
@@ -778,5 +80,5 @@
                   zhtot = zhtot + h_in(jk)
                   tabin(jk) = tabres(ji,jj,jk,m1)
-               ENDDO
+               END DO
                !          
                N_out = 0
@@ -785,5 +87,5 @@
                   N_out = N_out + 1
                   h_out(N_out) = e3v(ji,jj,jk,Kbb_a)
-               ENDDO
+               END DO
 
                ! Account for small differences in free-surface
@@ -797,6 +99,6 @@
                   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
-            ENDDO
-         ENDDO
+            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,:)  
@@ -804,8 +106,4 @@
          vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
 # endif
-         !* set relaxation time scale
-         IF( l_1st_euler .AND. lk_agrif_fstep ) THEN   ;   ztrelax =   rn_trelax_dyn  / (        rn_Dt )
-         ELSE                                          ;   ztrelax =   rn_trelax_dyn  / (2._wp * rn_Dt )
-         ENDIF
          !
          DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -817,5 +115,5 @@
             DO jj = j1+1,j2
                DO ji = i1,i2   ! vector opt.
-                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj)
+                  zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a)
                   hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kbb_a) * vbdiff(ji,jj  ,jk)  &
                                      &  -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kbb_a) * vbdiff(ji,jj-1,jk)  ) * zbtr
@@ -824,5 +122,5 @@
             DO jj = j1,j2
                DO ji = i1,i2-1   ! vector opt.
-                  zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj)
+                  zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk) 
                   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
@@ -844,5 +142,5 @@
                IF( .NOT. tabspongedone_u(ji,jj) ) THEN
                   DO jk = 1, jpkm1
-                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a)                                                               &
+                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a)                                                     &
                         & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )  &
                         & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj)
@@ -858,8 +156,8 @@
                IF( .NOT. tabspongedone_v(ji,jj) ) THEN
                   DO jk = 1, jpkm1
-                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a)                                                                  &
+                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a)                                                        &
                         &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) )   &
-                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)                      &
-                        &  - ztrelax * fspv(ji,jj) * vbdiff(ji,jj,jk)
+                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)                          &
+                        &  - rn_trelax_dyn * r1_Dt * fspv(ji,jj) * vbdiff(ji,jj,jk)
                   END DO
                ENDIF
Index: /NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN/sshwzv.F90	(revision 13229)
+++ /NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN/sshwzv.F90	(revision 13230)
@@ -200,5 +200,4 @@
       ENDIF
       !
-#if defined key_agrif 
       IF( .NOT. AGRIF_Root() ) THEN
          !
@@ -207,5 +206,5 @@
          DO jk = 1, jpkm1
             ! --- West --- !
-            IF( lk_west ) THEN
+            IF( lk_west) THEN
                DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
                   DO jj = 1, jpj
@@ -216,5 +215,5 @@
             !
             ! --- East --- !
-            IF( lk_east ) THEN
+            IF( lk_east) THEN
                DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
                   DO jj = 1, jpj
@@ -225,5 +224,5 @@
             !
             ! --- South --- !
-            IF( lk_south ) THEN
+            IF( lk_south) THEN
                DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
                   DO ji = 1, jpi
@@ -234,5 +233,5 @@
             !
             ! --- North --- !
-            IF( lk_north ) THEN
+            IF( lk_north) THEN
                DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
                   DO ji = 1, jpi
@@ -241,8 +240,8 @@
                END DO
             ENDIF
+            !
          END DO
          !
       ENDIF 
-#endif
       !
       IF( ln_timing )   CALL timing_stop('wzv')
