Index: branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90
===================================================================
--- branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90	(revision 8997)
+++ branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90	(revision 8998)
@@ -17,5 +17,7 @@
 
    PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90
-
+#if defined key_vertical
+   PUBLIC reconstructandremap ! remapping routine
+#endif   
    !                                              !!* Namelist namagrif: AGRIF parameters
    LOGICAL , PUBLIC ::   ln_spc_dyn    = .FALSE.   !:
@@ -24,6 +26,5 @@
    REAL(wp), PUBLIC ::   rn_sponge_dyn = 2800.     !: sponge coeff. for dynamics
    LOGICAL , PUBLIC ::   ln_chk_bathy  = .FALSE.   !: check of parent bathymetry 
-   LOGICAL , PUBLIC ::   lk_agrif_clp  = .FALSE.   !: Flag to retrieve clamped open boundaries
-
+   LOGICAL , PUBLIC ::   lk_agrif_clp  = .TRUE.    !: Force clamped bcs
    !                                              !!! OLD namelist names
    REAL(wp), PUBLIC ::   visc_tra                  !: sponge coeff. for tracers
@@ -101,4 +102,222 @@
    END FUNCTION agrif_oce_alloc
 
+#if defined key_vertical
+   SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)      
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION reconstructandremap  ***
+      !!----------------------------------------------------------------------
+      IMPLICIT NONE
+      INTEGER N, Nout
+      REAL(wp) tabin(N), tabout(Nout)
+      REAL(wp) hin(N), hout(Nout)
+      REAL(wp) coeffremap(N,3),zwork(N,3)
+      REAL(wp) zwork2(N+1,3)
+      INTEGER jk
+      DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8  
+      REAL(wp) q,q01,q02,q001,q002,q0
+      REAL(wp) z_win(1:N+1), z_wout(1:Nout+1)
+      REAL(wp),PARAMETER :: dpthin = 1.D-3
+      INTEGER :: k1, kbox, ktop, ka, kbot
+      REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
+
+      z_win(1)=0.; z_wout(1)= 0.
+      DO jk=1,N
+         z_win(jk+1)=z_win(jk)+hin(jk)
+      ENDDO 
+      
+      DO jk=1,Nout
+         z_wout(jk+1)=z_wout(jk)+hout(jk)       
+      ENDDO       
+
+      DO jk=2,N
+         zwork(jk,1)=1./(hin(jk-1)+hin(jk))
+      ENDDO
+        
+      DO jk=2,N-1
+         q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1))
+         zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1)
+         zwork(jk,3)=q0
+      ENDDO       
+     
+      DO jk= 2,N
+         zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1))
+      ENDDO
+        
+      coeffremap(:,1) = tabin(:)
+ 
+      DO jk=2,N-1
+         q001 = hin(jk)*zwork2(jk+1,1)
+         q002 = hin(jk)*zwork2(jk,1)        
+         IF (q001*q002 < 0) then
+            q001 = 0.
+            q002 = 0.
+         ENDIF
+         q=zwork(jk,2)
+         q01=q*zwork2(jk+1,1)
+         q02=q*zwork2(jk,1)
+         IF (abs(q001) > abs(q02)) q001 = q02
+         IF (abs(q002) > abs(q01)) q002 = q01
+        
+         q=(q001-q002)*zwork(jk,3)
+         q001=q001-q*hin(jk+1)
+         q002=q002+q*hin(jk-1)
+        
+         coeffremap(jk,3)=coeffremap(jk,1)+q001
+         coeffremap(jk,2)=coeffremap(jk,1)-q002
+        
+         zwork2(jk,1)=(2.*q001-q002)**2
+         zwork2(jk,2)=(2.*q002-q001)**2
+      ENDDO
+        
+      DO jk=1,N
+         IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN
+            coeffremap(jk,3) = coeffremap(jk,1)
+            coeffremap(jk,2) = coeffremap(jk,1)
+            zwork2(jk,1) = 0.
+            zwork2(jk,2) = 0.
+         ENDIF
+      ENDDO
+        
+      DO jk=2,N
+         q002=max(zwork2(jk-1,2),dsmll)
+         q001=max(zwork2(jk,1),dsmll)
+         zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002)
+      ENDDO
+        
+      zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
+      zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)
+ 
+      DO jk=1,N
+         q01=zwork2(jk+1,3)-coeffremap(jk,1)
+         q02=coeffremap(jk,1)-zwork2(jk,3)
+         q001=2.*q01
+         q002=2.*q02
+         IF (q01*q02<0) then
+            q01=0.
+            q02=0.
+         ELSEIF (abs(q01)>abs(q002)) then
+            q01=q002
+         ELSEIF (abs(q02)>abs(q001)) then
+            q02=q001
+         ENDIF
+         coeffremap(jk,2)=coeffremap(jk,1)-q02
+         coeffremap(jk,3)=coeffremap(jk,1)+q01
+      ENDDO
+
+      zbot=0.0
+      kbot=1
+      DO jk=1,Nout
+         ztop=zbot  !top is bottom of previous layer
+         ktop=kbot
+         IF     (ztop.GE.z_win(ktop+1)) then
+            ktop=ktop+1
+         ENDIF
+        
+         zbot=z_wout(jk+1)
+         zthk=zbot-ztop
+
+         IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN
+
+            kbot=ktop
+            DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)
+               kbot=kbot+1
+            ENDDO
+            zbox=zbot
+            DO k1= jk+1,Nout
+               IF     (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN
+                  exit !thick layer
+               ELSE
+                  zbox=z_wout(k1+1)  !include thin adjacent layers
+                  IF(zbox.EQ.z_wout(Nout+1)) THEN
+                     exit !at bottom
+                  ENDIF
+               ENDIF
+            ENDDO
+            zthk=zbox-ztop
+
+            kbox=ktop
+            DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)
+               kbox=kbox+1
+            ENDDO
+          
+            IF(ktop.EQ.kbox) THEN
+               IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN
+                  IF(hin(kbox).GT.dpthin) THEN
+                     q001 = (zbox-z_win(kbox))/hin(kbox)
+                     q002 = (ztop-z_win(kbox))/hin(kbox)
+                     q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)
+                     q02=q01-1.+(q001+q002)
+                     q0=1.-q01-q02
+                  ELSE
+                     q0 = 1.0
+                     q01 = 0.
+                     q02 = 0.
+                  ENDIF
+                  tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
+               ELSE
+                  tabout(jk) = tabin(kbox)
+               ENDIF 
+            ELSE
+               IF(ktop.LE.jk .AND. kbox.GE.jk) THEN
+                  ka = jk
+               ELSEIF (kbox-ktop.GE.3) THEN
+                  ka = (kbox+ktop)/2
+               ELSEIF (hin(ktop).GE.hin(kbox)) THEN
+                  ka = ktop
+               ELSE
+                  ka = kbox
+               ENDIF !choose ka
+    
+               offset=coeffremap(ka,1)
+    
+               qtop = z_win(ktop+1)-ztop !partial layer thickness
+               IF(hin(ktop).GT.dpthin) THEN
+                  q=(ztop-z_win(ktop))/hin(ktop)
+                  q01=q*(q-1.)
+                  q02=q01+q
+                  q0=1-q01-q02            
+               ELSE
+                  q0 = 1.
+                  q01 = 0.
+                  q02 = 0.
+               ENDIF
+               
+               tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop
+    
+               DO k1= ktop+1,kbox-1
+                  tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
+               ENDDO !k1
+               
+               qbot = zbox-z_win(kbox) !partial layer thickness
+               IF(hin(kbox).GT.dpthin) THEN
+                  q=qbot/hin(kbox)
+                  q01=(q-1.)**2
+                  q02=q01-1.+q
+                  q0=1-q01-q02                            
+               ELSE
+                  q0 = 1.0
+                  q01 = 0.
+                  q02 = 0.
+               ENDIF
+              
+               tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot
+               
+               rpsum=1.0d0/zthk
+               tabout(jk)=offset+tsum*rpsum
+                 
+            ENDIF !single or multiple layers
+         ELSE
+            IF (jk==1) THEN
+               write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1)
+            ENDIF
+            tabout(jk) = tabout(jk-1)
+             
+         ENDIF !normal:thin layer
+      ENDDO !jk
+            
+      return
+      end subroutine reconstructandremap
+#endif
+
 #endif
    !!======================================================================
Index: branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
===================================================================
--- branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 8997)
+++ branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 8998)
@@ -118,5 +118,5 @@
          !
          IF (.NOT.lk_agrif_clp) THEN
-            DO jk=1,jpkm1              ! Smooth
+            DO jk=1,jpkm1                 ! Smooth
                DO jj=j1,j2
                   ua(2,jj,jk) = 0.25_wp*(ua(1,jj,jk)+2._wp*ua(2,jj,jk)+ua(3,jj,jk))
@@ -124,5 +124,5 @@
                END DO
             END DO
-         END IF
+         ENDIF
          !
          zub(2,:) = 0._wp              ! Correct transport
@@ -189,5 +189,5 @@
 
          IF (.NOT.lk_agrif_clp) THEN
-            DO jk = 1, jpkm1           ! Smooth
+            DO jk = 1, jpkm1              ! Smooth
                DO jj = j1, j2
                   ua(nlci-2,jj,jk) = 0.25_wp * umask(nlci-2,jj,jk)      &
@@ -196,5 +196,4 @@
             END DO
          ENDIF
-
          zub(nlci-2,:) = 0._wp        ! Correct transport
          DO jk = 1, jpkm1
@@ -331,5 +330,5 @@
          !
          IF (.NOT.lk_agrif_clp) THEN
-            DO jk = 1, jpkm1           ! Smooth
+            DO jk = 1, jpkm1              ! Smooth
                DO ji = i1, i2
                   va(ji,nlcj-2,jk) = 0.25_wp * vmask(ji,nlcj-2,jk)   &
@@ -337,5 +336,5 @@
                END DO
             END DO
-         ENDIF
+         END IF
          !
          zvb(:,nlcj-2) = 0._wp         ! Correct transport
@@ -614,14 +613,78 @@
       INTEGER                                     , INTENT(in   ) ::   nb , ndir
       !
-      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      INTEGER  ::   imin, imax, jmin, jmax
+      INTEGER  ::   ji, jj, jk, jn, iref, jref   ! dummy loop indices
+      INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out
       REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
       REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
-      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
-      !!----------------------------------------------------------------------
-      !
+      LOGICAL :: western_side, eastern_side,northern_side,southern_side
+      ! vertical interpolation:
+      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
+      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
+      REAL(wp), DIMENSION(k1:k2) :: h_in
+      REAL(wp), DIMENSION(1:jpk) :: h_out(1:jpk)
+      REAL(wp) :: h_diff, zrhoxy
+
+      zrhoxy = Agrif_rhox()*Agrif_rhoy()
       IF (before) THEN         
-         ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
-      ELSE
+         DO jn = 1,jpts
+            DO jk=k1,k2
+               DO jj=j1,j2
+                 DO ji=i1,i2
+                       ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)
+                 END DO
+              END DO
+           END DO
+        END DO
+
+# if defined key_vertical
+        DO jk=k1,k2
+           DO jj=j1,j2
+              DO ji=i1,i2
+                 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
+              END DO
+           END DO
+        END DO
+# endif
+      ELSE 
+
+         western_side  = (nb == 1).AND.(ndir == 1)
+         eastern_side  = (nb == 1).AND.(ndir == 2)
+         southern_side = (nb == 2).AND.(ndir == 1)
+         northern_side = (nb == 2).AND.(ndir == 2)
+
+# if defined key_vertical              
+         DO jj=j1,j2
+            DO ji=i1,i2
+               iref = ji
+               jref = jj
+               if(western_side) iref=MAX(2,ji)
+               if(eastern_side) iref=MIN(nlci-1,ji)
+               if(southern_side) jref=MAX(2,jj)
+               if(northern_side) jref=MIN(nlcj-1,jj)
+               N_in = 0
+               DO jk=k1,k2 !k2 = jpk of parent grid
+                  IF (ptab(ji,jj,jk,n2) == 0) EXIT
+                  N_in = N_in + 1
+                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
+                  h_in(N_in) = ptab(ji,jj,jk,n2)
+               END DO
+               N_out = 0
+               DO jk=1,jpk ! jpk of child grid
+                  IF (tmask(iref,jref,jk) == 0) EXIT 
+                  N_out = N_out + 1
+                  h_out(jk) = e3t_n(iref,jref,jk)
+               ENDDO
+               IF (N_in > 0) THEN
+                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                  DO jn=1,jpts
+                     call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
+                  ENDDO
+               ENDIF
+            ENDDO
+         ENDDO
+# else
+         ptab_child(i1:i2,j1:j2,1:jpk,1:jpts) = ptab(i1:i2,j1:j2,1:jpk,1:jpts)
+# endif
+         !
          IF (lk_agrif_clp) THEN
             DO jn = 1, jpts
@@ -629,5 +692,5 @@
                   DO ji = i1,i2
                      DO jj = j1,j2
-                        tsa(ji,jj,jk,jn) = ptab(ji,jj,jk,jn)
+                        tsa(ji,jj,jk,jn) = ptab_child(ji,jj,jk,jn)
                      END DO
                   END DO
@@ -636,9 +699,4 @@
             return
          ENDIF
-         !
-         western_side  = (nb == 1).AND.(ndir == 1)
-         eastern_side  = (nb == 1).AND.(ndir == 2)
-         southern_side = (nb == 2).AND.(ndir == 1)
-         northern_side = (nb == 2).AND.(ndir == 2)
          !
          zrhox = Agrif_Rhox()
@@ -667,5 +725,5 @@
          IF( eastern_side ) THEN
             DO jn = 1, jpts
-               tsa(nlci,j1:j2,k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn)
+               tsa(nlci,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(nlci,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(nlci-1,j1:j2,1:jpk,jn)
                DO jk = 1, jpkm1
                   DO jj = jmin,jmax
@@ -687,5 +745,5 @@
          IF( northern_side ) THEN            
             DO jn = 1, jpts
-               tsa(i1:i2,nlcj,k1:k2,jn) = zalpha1 * ptab(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn)
+               tsa(i1:i2,nlcj,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,nlcj,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,nlcj-1,1:jpk,jn)
                DO jk = 1, jpkm1
                   DO ji = imin,imax
@@ -707,5 +765,5 @@
          IF( western_side ) THEN            
             DO jn = 1, jpts
-               tsa(1,j1:j2,k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn)
+               tsa(1,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(1,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(2,j1:j2,1:jpk,jn)
                DO jk = 1, jpkm1
                   DO jj = jmin,jmax
@@ -726,5 +784,5 @@
          IF( southern_side ) THEN           
             DO jn = 1, jpts
-               tsa(i1:i2,1,k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn)
+               tsa(i1:i2,1,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,1,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,2,1:jpk,jn)
                DO jk = 1, jpk      
                   DO ji=imin,imax
@@ -747,17 +805,17 @@
          ! East south
          IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
-            tsa(nlci-1,2,:,:) = ptab(nlci-1,2,:,:)
+            tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,1:jpts)
          ENDIF
          ! East north
          IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
-            tsa(nlci-1,nlcj-1,:,:) = ptab(nlci-1,nlcj-1,:,:)
+            tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,1:jpts)
          ENDIF
          ! West south
          IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
-            tsa(2,2,:,:) = ptab(2,2,:,:)
+            tsa(2,2,:,:) = ptab_child(2,2,:,1:jpts)
          ENDIF
          ! West north
          IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
-            tsa(2,nlcj-1,:,:) = ptab(2,nlcj-1,:,:)
+            tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,1:jpts)
          ENDIF
          !
@@ -765,5 +823,4 @@
       !
    END SUBROUTINE interptsn
-
 
    SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir )
@@ -794,58 +851,192 @@
    END SUBROUTINE interpsshn
 
-
-   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, before )
+   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
       !!----------------------------------------------------------------------
       !!   *** ROUTINE interpun ***
-      !!----------------------------------------------------------------------
-      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
-      LOGICAL                               , INTENT(in   ) ::   before
-      !
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zrhoy	
-      !!----------------------------------------------------------------------
-      !
-      IF( before ) THEN 
-         DO jk = 1, jpkm1
-            ptab(i1:i2,j1:j2,jk) = e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
+      !!---------------------------------------------    
+      !!
+      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
+      LOGICAL, INTENT(in) :: before
+      INTEGER, INTENT(in) :: nb , ndir
+      !!
+      INTEGER :: ji,jj,jk
+      REAL(wp) :: zrhoy
+      ! vertical interpolation:
+      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
+      REAL(wp), DIMENSION(1:jpk) :: h_out
+      INTEGER  :: N_in, N_out, iref
+      REAL(wp) :: h_diff
+      LOGICAL  :: western_side, eastern_side
+      !!---------------------------------------------    
+      !
+      zrhoy = Agrif_rhoy()
+      IF (before) THEN 
+         !We can't use zero as the special value because we need to include zeros
+         !when interpolating the scale factors
+         IF(Agrif_UseSpecialValue) THEN 
+!             Agrif_SpecialValue = -999._wp
+             Agrif_SpecialValue = 0._wp
+         ELSE
+             Agrif_SpecialValue = 0._wp
+         ENDIF
+         DO jk=1,jpk
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) - &
+                                   & ((umask(ji,jj,jk)-1) * Agrif_SpecialValue)
+# if defined key_vertical
+                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk))
+# endif
+               END DO
+            END DO
          END DO
       ELSE
-         zrhoy = Agrif_Rhoy()
+         zrhoy = Agrif_rhoy()
+# if defined key_vertical
+! VERTICAL REFINEMENT BEGIN
+         western_side  = (nb == 1).AND.(ndir == 1)
+         eastern_side  = (nb == 1).AND.(ndir == 2)
+
+         Agrif_SpecialValue = 0._wp ! reset specialvalue to zero now interpolation completed
+
+         DO ji=i1,i2
+            iref = ji
+            IF (western_side) iref = MAX(2,ji)
+            IF (eastern_side) iref = MIN(nlci-2,ji)
+            DO jj=j1,j2
+               N_in = 0
+               DO jk=k1,k2
+                  IF (ptab(ji,jj,jk,2) == 0) EXIT
+                  N_in = N_in + 1
+                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
+                  h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
+              ENDDO
+         
+              IF (N_in == 0) THEN
+                 ua(ji,jj,:) = 0._wp
+                 CYCLE
+              ENDIF
+         
+              N_out = 0
+              DO jk=1,jpk
+                 if (umask(iref,jj,jk) == 0) EXIT
+                 N_out = N_out + 1
+                 h_out(N_out) = e3u_a(iref,jj,jk)
+              ENDDO
+         
+              IF (N_out == 0) THEN
+                 ua(ji,jj,:) = 0._wp
+                 CYCLE
+              ENDIF
+         
+              IF (N_in * N_out > 0) THEN
+                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
+                 if (h_diff < -1.e4) then
+                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
+!                    stop
+                 endif
+              ENDIF
+              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
+            ENDDO
+         ENDDO
+
+# else
          DO jk = 1, jpkm1
             DO jj=j1,j2
-               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) )
-            END DO
-         END DO
+               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) )
+            END DO
+         END DO
+# endif
+
       ENDIF
       ! 
    END SUBROUTINE interpun
 
-
-   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, before )
+   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
       !!----------------------------------------------------------------------
       !!   *** ROUTINE interpvn ***
       !!----------------------------------------------------------------------
-      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
-      LOGICAL                               , INTENT(in   ) ::   before
-      !
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zrhox	
-      !!----------------------------------------------------------------------
+      !
+      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
+      LOGICAL, INTENT(in) :: before
+      INTEGER, INTENT(in) :: nb , ndir
+      !
+      INTEGER :: ji,jj,jk
+      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, jref
+      REAL(wp) :: h_diff
+      LOGICAL  :: northern_side,southern_side
+      !!---------------------------------------------    
       !      
-      IF( before ) THEN  
+      IF (before) THEN          
+         IF(Agrif_UseSpecialValue) THEN 
+!             Agrif_SpecialValue = -999._wp
+             Agrif_SpecialValue = 0._wp
+         ELSE
+             Agrif_SpecialValue = 0._wp
+         ENDIF
+         DO jk=k1,k2
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk)) - &
+                                   & ((vmask(ji,jj,jk)-1) * Agrif_SpecialValue)
+# if defined key_vertical
+                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)
+# endif
+               END DO
+            END DO
+         END DO
+      ELSE       
+         zrhox = Agrif_rhox()
+# if defined key_vertical
+         Agrif_SpecialValue = 0._wp !Reset special value to zero now interpolation is done
+
+         southern_side = (nb == 2).AND.(ndir == 1)
+         northern_side = (nb == 2).AND.(ndir == 2)
+
+         DO jj=j1,j2
+            jref = jj
+            IF (southern_side) jref = MAX(2,jj)
+            IF (northern_side) jref = MIN(nlcj-2,jj)
+            DO ji=i1,i2
+               N_in = 0
+               DO jk=k1,k2
+                  if (ptab(ji,jj,jk,2) == 0) EXIT
+                  N_in = N_in + 1
+                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
+                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)
+               END DO
+               IF (N_in == 0) THEN
+                  va(ji,jj,:) = 0._wp
+                  CYCLE
+               ENDIF
+         
+               N_out = 0
+               DO jk=1,jpk
+                  if (vmask(ji,jref,jk) == 0) EXIT
+                  N_out = N_out + 1
+                  h_out(N_out) = e3v_a(ji,jref,jk)
+               END DO
+               IF (N_out == 0) THEN
+                 va(ji,jj,:) = 0._wp
+                 CYCLE
+               ENDIF
+               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
+            END DO
+         END DO
+# else
          DO jk = 1, jpkm1
-            ptab(i1:i2,j1:j2,jk) = e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) * vn(i1:i2,j1:j2,jk)
-         END DO
-      ELSE          
-         zrhox= Agrif_Rhox()
-         DO jk = 1, jpkm1
-            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) )
-         END DO
+            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) )
+         END DO
+# endif
       ENDIF
       !        
    END SUBROUTINE interpvn
-   
 
    SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir )
Index: branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90
===================================================================
--- branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90	(revision 8997)
+++ branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90	(revision 8998)
@@ -39,5 +39,5 @@
 #if defined SPONGE
 !!      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
-!! Assume persistence:
+      !! Assume persistence:
       timecoeff = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
 
@@ -195,5 +195,4 @@
    END SUBROUTINE Agrif_Sponge
 
-
    SUBROUTINE interptsn_sponge(tabres,i1,i2,j1,j2,k1,k2,n1,n2,before)
       !!---------------------------------------------
@@ -207,12 +206,77 @@
       INTEGER  ::   iku, ikv
       REAL(wp) :: ztsa, zabe1, zabe2, zbtr
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ztu, ztv
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::tsbdiff
+      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
+      REAL(wp) :: h_diff
       !
       IF( before ) THEN
-         tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsb(i1:i2,j1:j2,k1:k2,n1:n2)
+         DO jn = 1, jpts
+            DO jk=k1,k2
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     tabres(ji,jj,jk,jn) = tsb(ji,jj,jk,jn)
+                  END DO
+               END DO
+            END DO
+         END DO
+
+# if defined key_vertical
+         DO jk=k1,k2
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
+               END DO
+            END DO
+         END DO
+# endif
+
       ELSE   
          !
-         tsbdiff(:,:,:,:) = tsb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:)    
+# if defined key_vertical
+         tabres_child(:,:,:,:) = 0.
+         DO jj=j1,j2
+            DO ji=i1,i2
+               N_in = 0
+               DO jk=k1,k2 !k2 = jpk of parent grid
+                  IF (tabres(ji,jj,jk,n2) == 0) EXIT
+                  N_in = N_in + 1
+                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)
+                  h_in(N_in) = tabres(ji,jj,jk,n2)
+               END DO
+               N_out = 0
+               DO jk=1,jpk ! jpk of child grid
+                  IF (tmask(ji,jj,jk) == 0) EXIT 
+                  N_out = N_out + 1
+                  h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above
+               ENDDO
+               IF (N_in > 0) THEN
+                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                  tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for?????
+                  DO jn=1,jpts
+                     call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
+                  ENDDO
+               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) = tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)
+# else
+                  tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts)
+# endif
+               ENDDO
+            ENDDO
+         ENDDO
+
          DO jn = 1, jpts            
             DO jk = 1, jpkm1
@@ -261,26 +325,83 @@
    END SUBROUTINE interptsn_sponge
 
-
-   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2, before)
+   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
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
+      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
+      INTEGER :: ji,jj,jk,jmax
 
       ! sponge parameters 
-      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ubdiff
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
-      INTEGER :: jmax
+      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff
+      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
-         tabres = ub(i1:i2,j1:j2,:)
+         DO jk=1,jpkm1
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  tabres(ji,jj,jk,m1) = ub(ji,jj,jk)
+# if defined key_vertical
+                  tabres(ji,jj,jk,m2) = e3u_n(ji,jj,jk)*umask(ji,jj,jk)
+# endif
+               END DO
+            END DO
+         END DO
+
       ELSE
-         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:))*umask(i1:i2,j1:j2,:)
+
+# if defined key_vertical
+         tabres_child(:,:,:) = 0._wp
+         DO jj=j1,j2
+            DO ji=i1,i2
+               N_in = 0
+               DO jk=k1,k2
+                  IF (tabres(ji,jj,jk,m2) == 0) EXIT
+                  N_in = N_in + 1
+                  tabin(jk) = tabres(ji,jj,jk,m1)
+                  h_in(N_in) = tabres(ji,jj,jk,m2)
+              ENDDO
+              !
+              IF (N_in == 0) THEN
+                 tabres_child(ji,jj,:) = 0.
+                 CYCLE
+              ENDIF
+         
+              N_out = 0
+              DO jk=1,jpk
+                 if (umask(ji,jj,jk) == 0) EXIT
+                 N_out = N_out + 1
+                 h_out(N_out) = e3u_n(ji,jj,jk)
+              ENDDO
+         
+              IF (N_out == 0) THEN
+                 tabres_child(ji,jj,:) = 0.
+                 CYCLE
+              ENDIF
+         
+              IF (N_in * N_out > 0) THEN
+                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                 if (h_diff < -1.e4) then
+                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
+                 endif
+              ENDIF
+              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)
+         
+            ENDDO
+         ENDDO
+
+         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
+#else
+         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
+#endif
          !
          DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -359,26 +480,75 @@
    END SUBROUTINE interpun_sponge
 
-
-   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2, before,nb,ndir)
+   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
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
+      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
-      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: vbdiff
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
-      INTEGER :: imax
+      INTEGER  ::   ji, jj, jk, imax
+      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, h_diff
+      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 
-         tabres = vb(i1:i2,j1:j2,:)
+         DO jk=1,jpkm1
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  tabres(ji,jj,jk,m1) = vb(ji,jj,jk)
+# if defined key_vertical
+                  tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_n(ji,jj,jk)
+# endif
+               END DO
+            END DO
+         END DO
       ELSE
-         !
-         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:)
+
+# if defined key_vertical
+         tabres_child(:,:,:) = 0._wp
+         DO jj=j1,j2
+            DO ji=i1,i2
+               N_in = 0
+               DO jk=k1,k2
+                  IF (tabres(ji,jj,jk,m2) == 0) EXIT
+                  N_in = N_in + 1
+                  tabin(jk) = tabres(ji,jj,jk,m1)
+                  h_in(N_in) = tabres(ji,jj,jk,m2)
+              ENDDO
+         
+              IF (N_in == 0) THEN
+                 tabres_child(ji,jj,:) = 0.
+                 CYCLE
+              ENDIF
+         
+              N_out = 0
+              DO jk=1,jpk
+                 if (vmask(ji,jj,jk) == 0) EXIT
+                 N_out = N_out + 1
+                 h_out(N_out) = e3v_n(ji,jj,jk)
+              ENDDO
+         
+              IF (N_in * N_out > 0) THEN
+                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                 if (h_diff < -1.e4) then
+                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
+                 endif
+              ENDIF
+              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)
+            ENDDO
+         ENDDO
+
+         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)  
+# else
+         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
+# endif
          !
          DO jk = 1, jpkm1                                 ! Horizontal slab
Index: branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
===================================================================
--- branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 8997)
+++ branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 8998)
@@ -1,4 +1,4 @@
 #define TWO_WAY        /* TWO WAY NESTING */
-#undef DECAL_FEEDBACK  /* SEPARATION of INTERFACES*/
+#define DECAL_FEEDBACK /* SEPARATION of INTERFACES*/
 #undef VOL_REFLUX      /* VOLUME REFLUXING*/
  
@@ -97,5 +97,4 @@
       ! Account for updated thicknesses at boundary edges
       IF (.NOT.ln_linssh) THEN
-! For the time being calls below do not ensure reproducible results 
 !         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy)
 !         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy)
@@ -283,4 +282,111 @@
    END SUBROUTINE dom_vvl_update_UVF
 
+#if defined key_vertical
+
+   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
+      !!---------------------------------------------
+      !!           *** ROUTINE updateT ***
+      !!---------------------------------------------
+      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
+      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child
+      REAL(wp) :: h_in(k1:k2)
+      REAL(wp) :: h_out(1:jpk)
+      INTEGER  :: N_in, N_out
+      REAL(wp) :: h_diff
+      REAL(wp) :: zrho_xy
+      REAL(wp) :: tabin(k1:k2,n1:n2)
+      !!---------------------------------------------
+      !
+      IF (before) THEN
+         AGRIF_SpecialValue = -999._wp
+         zrho_xy = Agrif_rhox() * Agrif_rhoy() 
+         DO jn = n1,n2-1
+            DO jk=k1,k2
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) &
+                                           * tmask(ji,jj,jk) * zrho_xy + (tmask(ji,jj,jk)-1)*999._wp
+                  END DO
+               END DO
+            END DO
+         END DO
+         DO jk=k1,k2
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) * zrho_xy  &
+                                           + (tmask(ji,jj,jk)-1)*999._wp
+               END DO
+            END DO
+         END DO
+      ELSE
+         tabres_child(:,:,:,:) = 0.
+         AGRIF_SpecialValue = 0._wp
+         DO jj=j1,j2
+            DO ji=i1,i2
+               N_in = 0
+               DO jk=k1,k2 !k2 = jpk of child grid
+                  IF (tabres(ji,jj,jk,n2) == 0  ) EXIT
+                  N_in = N_in + 1
+                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
+                  h_in(N_in) = tabres(ji,jj,jk,n2)/e1e2t(ji,jj)
+               ENDDO
+               N_out = 0
+               DO jk=1,jpk ! jpk of parent grid
+                  IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF
+                  N_out = N_out + 1
+                  h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above
+               ENDDO
+               IF (N_in > 0) THEN !Remove this?
+                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                  IF (h_diff < -1.e-4) THEN
+                     print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out))
+                     print *,h_in(1:N_in)
+                     print *,h_out(1:N_out)
+                     STOP
+                  ENDIF
+                  DO jn=n1,n2-1
+                     CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out)
+                  ENDDO
+               ENDIF
+            ENDDO
+         ENDDO
+
+         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
+            ! Add asselin part
+            DO jn = n1,n2-1
+               DO jk=1,jpk
+                  DO jj=j1,j2
+                     DO ji=i1,i2
+                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
+                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
+                                 & + atfp * ( tabres_child(ji,jj,jk,jn) &
+                                 &          - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
+                        ENDIF
+                     ENDDO
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+         DO jn = n1,n2-1
+            DO jk=1,jpk
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
+                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
+                     END IF
+                  END DO
+               END DO
+            END DO
+         END DO
+      ENDIF
+      ! 
+   END SUBROUTINE updateTS
+
+# else
+
    SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
       !!---------------------------------------------
@@ -296,5 +402,5 @@
       !
       IF (before) THEN
-         DO jn = n1,n2
+         DO jn = 1,jpts
             DO jk=k1,k2
                DO jj=j1,j2
@@ -310,5 +416,5 @@
       ELSE
 !> jc tmp
-         DO jn = n1,n2
+         DO jn = 1,jpts
             tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
                                          & * tmask(i1:i2,j1:j2,k1:k2)
@@ -317,5 +423,5 @@
          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
             ! Add asselin part
-            DO jn = n1,n2
+            DO jn = 1,jpts
                DO jk=k1,k2
                   DO jj=j1,j2
@@ -333,5 +439,5 @@
             ENDDO
          ENDIF
-         DO jn = n1,n2
+         DO jn = 1,jpts
             DO jk=k1,k2
                DO jj=j1,j2
@@ -346,5 +452,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            tsb(i1:i2,j1:j2,k1:k2,n1:n2)  = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
+            tsb(i1:i2,j1:j2,k1:k2,1:jpts)  = tsn(i1:i2,j1:j2,k1:k2,1:jpts)
          ENDIF
          !
@@ -353,12 +459,110 @@
    END SUBROUTINE updateTS
 
-
-   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, before )
+# endif
+
+# if defined key_vertical
+
+   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
       !!---------------------------------------------
       !!           *** ROUTINE updateu ***
       !!---------------------------------------------
-      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
-      LOGICAL                               , INTENT(in   ) :: before
+      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
+      REAL(wp) ::   zrhoy
+! VERTICAL REFINEMENT BEGIN
+      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
+      REAL(wp) :: h_in(k1:k2)
+      REAL(wp) :: h_out(1:jpk)
+      INTEGER  :: N_in, N_out
+      REAL(wp) :: h_diff, excess, thick
+      REAL(wp) :: tabin(k1:k2)
+! VERTICAL REFINEMENT END
+      !!---------------------------------------------
+      ! 
+      IF( before ) THEN
+         zrhoy = Agrif_Rhoy()
+         AGRIF_SpecialValue = -999._wp
+         DO jk=k1,k2
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  &
+                                       + (umask(ji,jj,jk)-1)*999._wp
+                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  &
+                                       + (umask(ji,jj,jk)-1)*999._wp
+               END DO
+            END DO
+         END DO
+      ELSE
+         tabres_child(:,:,:) = 0.
+         AGRIF_SpecialValue = 0._wp
+         DO jj=j1,j2
+            DO ji=i1,i2
+               N_in = 0
+               h_in(:) = 0._wp
+               tabin(:) = 0._wp
+               DO jk=k1,k2 !k2=jpk of child grid
+                  IF( tabres(ji,jj,jk,2) < -900) EXIT
+                  N_in = N_in + 1
+                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
+                  h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
+               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_n(ji,jj,jk)
+               ENDDO
+               IF (N_in * N_out > 0) THEN
+                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                  IF (h_diff < -1.e-4) THEN
+!Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 
+!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
+                     excess = 0._wp
+                     DO jk=N_in,1,-1
+                        thick = MIN(-1*h_diff, h_in(jk))
+                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
+                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
+                        h_diff = h_diff + thick
+                        IF ( h_diff == 0) THEN
+                           N_in = jk
+                           h_in(jk) = h_in(jk) - thick
+                           EXIT
+                        ENDIF
+                     ENDDO
+                  ENDIF
+                  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)
+                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
+               ENDIF
+            ENDDO
+         ENDDO
+
+         DO jk=1,jpk
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
+                     ub(ji,jj,jk) = ub(ji,jj,jk) & 
+                           & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
+                  ENDIF
+                  !
+                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
+      ! 
+   END SUBROUTINE updateu
+
+#else
+
+   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
+      !!---------------------------------------------
+      !!           *** ROUTINE updateu ***
+      !!---------------------------------------------
+      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
@@ -369,5 +573,5 @@
          zrhoy = Agrif_Rhoy()
          DO jk = k1, k2
-            tabres(i1:i2,j1:j2,jk) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
+            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
          END DO
       ELSE
@@ -375,15 +579,15 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) 
+                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) 
                   !
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
                      zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)  ! fse3t_b prior update should be used
                      zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk)
-                     zunu = tabres(ji,jj,jk)
+                     zunu = tabres(ji,jj,jk,1)
                      ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &      
                                     & * umask(ji,jj,jk) / e3u_b(ji,jj,jk)
                   ENDIF
                   !
-                  un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)
+                  un(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)
                END DO
             END DO
@@ -398,12 +602,14 @@
    END SUBROUTINE updateu
 
-   SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, before, nb, ndir )
+# endif
+
+   SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
       !!---------------------------------------------
       !!           *** ROUTINE correct_u_bdy ***
       !!---------------------------------------------
-      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
-      LOGICAL                               , INTENT(in   ) :: before
-      INTEGER                               , INTENT(in)    :: nb, ndir
+      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                                     , INTENT(in)    :: nb, ndir
       !!
       LOGICAL :: western_side, eastern_side 
@@ -442,12 +648,107 @@
    END SUBROUTINE correct_u_bdy
 
-
-   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before)
+# if  defined key_vertical
+
+   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
       !!---------------------------------------------
       !!           *** ROUTINE updatev ***
       !!---------------------------------------------
-      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
-      LOGICAL                               , INTENT(in   ) :: before
+      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
+      REAL(wp) ::   zrhox
+! VERTICAL REFINEMENT BEGIN
+      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
+      REAL(wp) :: h_in(k1:k2)
+      REAL(wp) :: h_out(1:jpk)
+      INTEGER :: N_in, N_out
+      REAL(wp) :: h_diff, excess, thick
+      REAL(wp) :: tabin(k1:k2)
+! VERTICAL REFINEMENT END
+      !!---------------------------------------------      
+      !
+      IF (before) THEN
+         zrhox = Agrif_Rhox()
+         AGRIF_SpecialValue = -999._wp
+         DO jk=k1,k2
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) &
+                                       + (vmask(ji,jj,jk)-1)*999._wp
+                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) &
+                                       + (vmask(ji,jj,jk)-1)*999._wp
+               END DO
+            END DO
+         END DO
+      ELSE
+         tabres_child(:,:,:) = 0.
+         AGRIF_SpecialValue = 0._wp
+         DO jj=j1,j2
+            DO ji=i1,i2
+               N_in = 0
+               DO jk=k1,k2
+                  IF (tabres(ji,jj,jk,2) < -900) EXIT
+                  N_in = N_in + 1
+                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
+                  h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)
+               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_n(ji,jj,jk)
+               ENDDO
+               IF (N_in * N_out > 0) THEN
+                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                  IF (h_diff < -1.e-4) then
+!Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 
+!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
+                     excess = 0._wp
+                     DO jk=N_in,1,-1
+                        thick = MIN(-1*h_diff, h_in(jk))
+                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
+                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
+                        h_diff = h_diff + thick
+                        IF ( h_diff == 0) THEN
+                           N_in = jk
+                           h_in(jk) = h_in(jk) - thick
+                           EXIT
+                        ENDIF
+                     ENDDO
+                  ENDIF
+                  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)
+                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
+               ENDIF
+            ENDDO
+         ENDDO
+
+         DO jk=1,jpk
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  !
+                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
+                     vb(ji,jj,jk) = vb(ji,jj,jk) & 
+                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
+                  ENDIF
+                  !
+                  vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
+      ! 
+   END SUBROUTINE updatev
+
+# else
+
+   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)
+      !!---------------------------------------------
+      !!           *** ROUTINE updatev ***
+      !!---------------------------------------------
+      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
@@ -460,5 +761,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
+                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
                END DO
             END DO
@@ -468,15 +769,15 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj)
+                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj)
                   !
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
                      zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used
                      zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk)
-                     zvnu = tabres(ji,jj,jk)
+                     zvnu = tabres(ji,jj,jk,1)
                      vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &      
                                     & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk)
                   ENDIF
                   !
-                  vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)
+                  vn(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)
                END DO
             END DO
@@ -491,12 +792,14 @@
    END SUBROUTINE updatev
 
-   SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, before, nb, ndir )
+# endif
+
+   SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
       !!---------------------------------------------
       !!           *** ROUTINE correct_u_bdy ***
       !!---------------------------------------------
-      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
-      LOGICAL                               , INTENT(in   ) :: before
-      INTEGER                               , INTENT(in)    :: nb, ndir
+      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                                     , INTENT(in)    :: nb, ndir
       !!
       LOGICAL :: southern_side, northern_side 
Index: branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
===================================================================
--- branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 8997)
+++ branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 8998)
@@ -1,3 +1,3 @@
-#undef UPD_HIGH   /* MIX HIGH UPDATE */
+#define UPD_HIGH   /* MIX HIGH UPDATE */
 #if defined key_agrif
 !!----------------------------------------------------------------------
@@ -330,14 +330,13 @@
    IF (Agrif_Root()) RETURN
    !
-                       CALL Agrif_Update_ssh()
    IF (.NOT.ln_linssh) CALL Agrif_Update_vvl()
-                       CALL Agrif_Update_tra()
+   CALL Agrif_Update_tra()
 #if defined key_top
-                       CALL Agrif_Update_Trc()
+   CALL Agrif_Update_Trc()
 #endif
-                       CALL Agrif_Update_dyn()
+   CALL Agrif_Update_dyn()
 # if defined key_zdftke
 ! JC remove update because this precludes from perfect restartability
-!!                       CALL Agrif_Update_tke(0)
+!!   CALL Agrif_Update_tke(0)
 # endif
 
@@ -364,13 +363,25 @@
    ! 1. Declaration of the type of variable which have to be interpolated
    !---------------------------------------------------------------------
+# if defined key_vertical
+   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_id)
+   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_sponge_id)
+
+   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_interp_id)
+   CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_interp_id)
+   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_update_id)
+   CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_update_id)
+   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_sponge_id)
+   CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_sponge_id)
+# else
    CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_id)
    CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_sponge_id)
 
-   CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_interp_id)
-   CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_interp_id)
-   CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_update_id)
-   CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_update_id)
-   CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_sponge_id)
-   CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_sponge_id)
+   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_interp_id)
+   CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_interp_id)
+   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_update_id)
+   CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_update_id)
+   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_sponge_id)
+   CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_sponge_id)
+# endif
 
    CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),e3t_id)
@@ -392,5 +403,9 @@
    CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/), en_id)
    CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),avt_id)
+# if defined key_vertical
+   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/jpi,jpj,jpk,2/),avm_id)
+# else
    CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),avm_id)
+# endif
 # endif
 
@@ -777,6 +792,11 @@
    ! 1. Declaration of the type of variable which have to be interpolated
    !---------------------------------------------------------------------
+# if defined key_vertical
+   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_id)
+   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_sponge_id)
+# else
    CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_id)
    CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_sponge_id)
+# endif
 
    ! 2. Type of interpolation
