Index: /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce.F90
===================================================================
--- /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce.F90	(revision 11740)
+++ /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce.F90	(revision 11741)
@@ -46,4 +46,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ubdy, vbdy, hbdy
 
+# if defined key_vertical
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht0_parent, hu0_parent, hv0_parent
+   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt_parent, mbku_parent, mbkv_parent
+# endif
 
    INTEGER, PUBLIC :: tsn_id                                                  ! AGRIF profile for tracers interpolation and update
@@ -60,4 +64,5 @@
    INTEGER, PUBLIC :: avt_id, avm_id, en_id                ! TKE related identificators
    INTEGER, PUBLIC :: umsk_id, vmsk_id
+   INTEGER, PUBLIC :: mbkt_id, ht0_id
    INTEGER, PUBLIC :: kindic_agr
    
@@ -83,5 +88,10 @@
 # if defined key_top         
          &      tabspongedone_trn(jpi,jpj),           &
-# endif         
+# endif   
+# if defined key_vertical
+         &      ht0_parent(jpi,jpj), mbkt_parent(jpi,jpj),  &
+         &      hu0_parent(jpi,jpj), mbku_parent(jpi,jpj),  &
+         &      hv0_parent(jpi,jpj), mbkv_parent(jpi,jpj),  &
+# endif      
          &      tabspongedone_u  (jpi,jpj),           &
          &      tabspongedone_v  (jpi,jpj), STAT = ierr(1) )
Index: /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_interp.F90
===================================================================
--- /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_interp.F90	(revision 11740)
+++ /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_interp.F90	(revision 11741)
@@ -44,5 +44,7 @@
    PUBLIC   interpunb, interpvnb , interpub2b, interpvb2b
    PUBLIC   interpe3t, interpumsk, interpvmsk
-
+#if defined key_vertical
+   PUBLIC   interpht0, interpmbkt
+# endif
    INTEGER ::   bdy_tinterp = 0
 
@@ -662,5 +664,5 @@
       INTEGER  ::   N_in, N_out
       ! vertical interpolation:
-      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jpts) :: ptab_child
+      REAL(wp) :: zhtot
       REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin
       REAL(wp), DIMENSION(k1:k2) :: h_in
@@ -680,24 +682,50 @@
 
 # if defined key_vertical
-        DO jk=k1,k2
+        ! Interpolate thicknesses
+        ! Warning: these are masked, hence extrapolated prior interpolation.
+        DO jk=k1,k2-1
            DO jj=j1,j2
               DO ji=i1,i2
-                 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
+                  ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)
               END DO
            END DO
         END DO
+
+        ! Extrapolate thicknesses in partial bottom cells:
+        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
+        IF (ln_zps) THEN
+           DO jj=j1,j2
+              DO ji=i1,i2
+                  jk = 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) = sshn(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 
+        ELSE
+           ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp
+        END IF      
 # endif
       ELSE 
 
-# if defined key_vertical              
+# if defined key_vertical 
+         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 
+            
          DO jj=j1,j2
             DO ji=i1,i2
-               ptab_child(ji,jj,:,:) = 0._wp
-               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
+               tsa(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) + 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)
-                  h_in(N_in) = ptab(ji,jj,jk,n2)
                END DO
                N_out = 0
@@ -708,15 +736,14 @@
                ENDDO
                IF (N_in*N_out > 0) THEN
-                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ptab_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
+                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tsa(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
                ENDIF
             ENDDO
          ENDDO
 # else
-         ptab_child(i1:i2,j1:j2,1:jpk,1:jpts) = ptab(i1:i2,j1:j2,1:jpk,1:jpts)
-# endif
          !
          DO jn=1, jpts
-            tsa(i1:i2,j1:j2,1:jpk,jn)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
-         END DO
+            tsa(i1:i2,j1:j2,1:jpk,jn)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
+         END DO
+# endif
 
       ENDIF
@@ -752,5 +779,5 @@
       !!
       INTEGER :: ji,jj,jk
-      REAL(wp) :: zrhoy
+      REAL(wp) :: zrhoy, zhtot
       ! vertical interpolation:
       REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
@@ -766,9 +793,32 @@
                   ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) 
 # if defined key_vertical
-                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk))
+                  ! Interpolate thicknesses (masked for subsequent extrapolation)
+                  ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(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)
+                  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_n(i1:i2,j1:j2,jk) * 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
+         !
       ELSE
          zrhoy = Agrif_rhoy()
@@ -776,19 +826,21 @@
 ! VERTICAL REFINEMENT BEGIN
 
+         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
+
          DO ji=i1,i2
             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) 
+               ua(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) + 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
-         
-              IF (N_in == 0) THEN
-                 ua(ji,jj,:) = 0._wp
-                 CYCLE
-              ENDIF
-         
+                  
               N_out = 0
               DO jk=1,jpk
@@ -797,11 +849,7 @@
                  h_out(N_out) = e3u_a(ji,jj,jk)
               ENDDO
-         
-              IF (N_out == 0) THEN
-                 ua(ji,jj,:) = 0._wp
-                 CYCLE
+              IF (N_in*N_out > 0) THEN
+                 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,1)
               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,1)
             ENDDO
          ENDDO
@@ -834,5 +882,5 @@
       REAL(wp), DIMENSION(1:jpk) :: h_out
       INTEGER  :: N_in, N_out
-      REAL(wp) :: h_diff
+      REAL(wp) :: h_diff, zhtot
       !!---------------------------------------------    
       !      
@@ -843,4 +891,5 @@
                   ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*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_n(ji,jj,jk)
 # endif
@@ -848,21 +897,45 @@
             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 = 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_n(i1:i2,j1:j2,jk) * 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
-               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
+               va(ji,jj,:) = 0._wp
+               N_in = mbkv_parent(ji,jj)
+               zhtot = 0._wp
+               DO jk=1,N_in
+                  IF (jk==N_in) THEN
+                     h_in(jk) = hv0_parent(ji,jj) + 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)
+                  tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk))
+              ENDDO
          
                N_out = 0
@@ -872,9 +945,8 @@
                   h_out(N_out) = e3v_a(ji,jj,jk)
                END DO
-               IF (N_out == 0) THEN
-                 va(ji,jj,:) = 0._wp
-                 CYCLE
+
+               IF (N_in*N_out > 0) THEN
+                  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,1)
                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,1)
             END DO
          END DO
@@ -1257,4 +1329,42 @@
    END SUBROUTINE interpavm
 
+# if defined key_vertical
+   SUBROUTINE interpmbkt( 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) = REAL(mbkt(i1:i2,j1:j2),wp)
+      ELSE
+         mbkt_parent(i1:i2,j1:j2) = NINT(ptab(i1:i2,j1:j2))
+      ENDIF
+      !
+   END SUBROUTINE interpmbkt
+
+   SUBROUTINE interpht0( 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) = ht_0(i1:i2,j1:j2)
+      ELSE
+         ht0_parent(i1:i2,j1:j2) = ptab(i1:i2,j1:j2)
+      ENDIF
+      !
+   END SUBROUTINE interpht0
+#endif
+
 #else
    !!----------------------------------------------------------------------
Index: /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90
===================================================================
--- /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90	(revision 11740)
+++ /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90	(revision 11741)
@@ -258,5 +258,5 @@
          CALL lbc_lnk( 'agrif_Sponge', fsaht_spu, 'U', 1. )   ! Lateral boundary conditions
          CALL lbc_lnk( 'agrif_Sponge', fsaht_spv, 'V', 1. )
-         
+
          spongedoneT = .TRUE.
       ENDIF
@@ -293,5 +293,5 @@
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
       INTEGER  ::   iku, ikv
-      REAL(wp) :: ztsa, zabe1, zabe2, zbtr
+      REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot
       REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
       REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
@@ -302,5 +302,4 @@
       REAL(wp), DIMENSION(1:jpk) :: h_out
       INTEGER :: N_in, N_out
-      REAL(wp) :: h_diff
       !!----------------------------------------------------------------------
       !
@@ -317,11 +316,31 @@
 
 # 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_b(ji,jj,jk) 
-               END DO
-            END DO
-         END DO
+        ! Interpolate thicknesses
+        ! Warning: these are masked, hence extrapolated prior interpolation.
+        DO jk=k1,k2-1
+           DO jj=j1,j2
+              DO ji=i1,i2
+                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_b(ji,jj,jk)
+              END DO
+           END DO
+        END DO
+
+        ! Extrapolate thicknesses in partial bottom cells:
+        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
+        IF (ln_zps) THEN
+           DO jj=j1,j2
+              DO ji=i1,i2
+                  jk = 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) = sshb(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 
+        ELSE
+           tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp
+        END IF      
 # endif
 
@@ -329,13 +348,20 @@
          !
 # 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 = 0
-               DO jk=k1,k2 !k2 = jpk of parent grid
-                  IF (tabres(ji,jj,jk,n2) == 0) EXIT
-                  N_in = N_in + 1
+               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)
-                  h_in(N_in) = tabres(ji,jj,jk,n2)
                END DO
                N_out = 0
@@ -345,4 +371,11 @@
                   h_out(jk) = e3t_b(ji,jj,jk) !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)
@@ -356,7 +389,7 @@
                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)
+                  tsbdiff(ji,jj,jk,1:jpts) = (tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk)
 # else
-                  tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts)
+                  tsbdiff(ji,jj,jk,1:jpts) = (tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk)
 # endif
                ENDDO
@@ -427,5 +460,5 @@
 
       ! sponge parameters 
-      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff
+      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
@@ -434,5 +467,5 @@
       REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
       REAL(wp), DIMENSION(1:jpk) :: h_out
-      INTEGER ::N_in,N_out
+      INTEGER ::N_in, N_out
       !!---------------------------------------------    
       !
@@ -449,43 +482,63 @@
          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_b(i1:i2,j1:j2,jk) * 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
-         tabres_child(:,:,:) = 0._wp
+         IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 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
+               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)
-                  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_b(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,1)
-         
+               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_b(ji,jj,jk)
+               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
@@ -580,5 +633,5 @@
       !
       INTEGER  ::   ji, jj, jk, imax
-      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, h_diff
+      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, zhtot
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
@@ -601,37 +654,63 @@
             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 = mbkv(ji,jj)
+                  tabres(ji,jj,jk,m2) = 0._wp
+               END DO
+            END DO           
+         END IF
+        ! Save ssh at last level:
+        tabres(i1:i2,j1:j2,k2,m2) = 0._wp
+        IF (.NOT.ln_linssh) THEN
+           ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
+           DO jk=1,jpk
+              tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v_b(i1:i2,j1:j2,jk) * vmask(i1:i2,j1:j2,jk)
+           END DO
+           tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2)
+        END IF 
+# endif
+
       ELSE
 
 # if defined key_vertical
-         tabres_child(:,:,:) = 0._wp
+         IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 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
+               tabres_child(ji,jj,:) = 0._wp
+               N_in = mbkv_parent(ji,jj)
+               zhtot = 0._wp
+               DO jk=1,N_in
+                  IF (jk==N_in) THEN
+                     h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
+                  ELSE
+                     h_in(jk) = tabres(ji,jj,jk,m2)
+                  ENDIF
+                  zhtot = zhtot + h_in(jk)
                   tabin(jk) = tabres(ji,jj,jk,m1)
-                  h_in(N_in) = tabres(ji,jj,jk,m2)
-              ENDDO
+               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_b(ji,jj,jk)
+               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 == 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_b(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,1)
+               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
Index: /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_update.F90
===================================================================
--- /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_update.F90	(revision 11740)
+++ /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_update.F90	(revision 11741)
@@ -1,4 +1,5 @@
-#undef DECAL_FEEDBACK  /* SEPARATION of INTERFACES*/
-#undef VOL_REFLUX      /* VOLUME REFLUXING*/
+#undef DECAL_FEEDBACK     /* SEPARATION of INTERFACES */
+#undef DECAL_FEEDBACK_2D  /* SEPARATION of INTERFACES (Barotropic mode) */
+#undef VOL_REFLUX         /* VOLUME REFLUXING*/
  
 MODULE agrif_oce_update
@@ -48,4 +49,5 @@
       IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed()
 
+! jc_alt      Agrif_UseSpecialValueInUpdate = .FALSE.
       Agrif_UseSpecialValueInUpdate = .TRUE.
       Agrif_SpecialValueFineGrid    = 0._wp
@@ -92,5 +94,5 @@
 # endif
 
-# if ! defined DECAL_FEEDBACK
+# if ! defined DECAL_FEEDBACK_2D
       CALL Agrif_Update_Variable(e1u_id,procname = updateU2d)
       CALL Agrif_Update_Variable(e2v_id,procname = updateV2d)  
@@ -100,5 +102,5 @@
 # endif
       !
-# if ! defined DECAL_FEEDBACK
+# if ! defined DECAL_FEEDBACK_2D
       ! Account for updated thicknesses at boundary edges
       IF (.NOT.ln_linssh) THEN
@@ -110,5 +112,5 @@
       IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN
          ! Update time integrated transports
-#  if ! defined DECAL_FEEDBACK
+#  if ! defined DECAL_FEEDBACK_2D
          CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b)
          CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b)
@@ -130,5 +132,5 @@
       Agrif_UseSpecialValueInUpdate = .TRUE.
       Agrif_SpecialValueFineGrid = 0.
-# if ! defined DECAL_FEEDBACK
+# if ! defined DECAL_FEEDBACK_2D
       CALL Agrif_Update_Variable(sshn_id,procname = updateSSH)
 # else
@@ -141,5 +143,5 @@
       IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
          ! Refluxing on ssh:
-#  if defined DECAL_FEEDBACK
+#  if defined DECAL_FEEDBACK_2D
          CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu)
          CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv)
@@ -299,5 +301,6 @@
                   DO ji=i1,i2
                      tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) &
-                                           * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp)*999._wp
+                                         &  * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp
+!jc_alt                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk)
                   END DO
                END DO
@@ -307,6 +310,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) &
-                                           + (tmask(ji,jj,jk)-1)*999._wp
+                  tabres(ji,jj,jk,n2) =      tmask(ji,jj,jk) * e3t_n(ji,jj,jk) &
+                                      &   + (tmask(ji,jj,jk) - 1._wp) * 999._wp
+! jc_alt                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)
                END DO
             END DO
@@ -320,4 +324,5 @@
                DO jk=k1,k2 !k2 = jpk of child grid
                   IF (tabres(ji,jj,jk,n2) < -900._wp  ) EXIT
+! jc_alt                  IF (tabres(ji,jj,jk,n2) == 0._wp  ) EXIT
                   N_in = N_in + 1
                   tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
@@ -476,7 +481,9 @@
                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
+                                     &  + (umask(ji,jj,jk)-1._wp)*999._wp
+! jc_alt                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  
                   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
+                                     &  + (umask(ji,jj,jk)-1._wp)*999._wp
+! jc_alt                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)
                END DO
             END DO
@@ -491,5 +498,6 @@
                tabin(:) = 0._wp
                DO jk=k1,k2 !k2=jpk of child grid
-                  IF( tabres(ji,jj,jk,2) < -900) EXIT
+                  IF( tabres(ji,jj,jk,2) < -900._wp) EXIT
+! jc_alt                  IF( tabres(ji,jj,jk,2) == 0.) EXIT
                   N_in = N_in + 1
                   tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
@@ -671,7 +679,9 @@
                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
+                                     & + (vmask(ji,jj,jk)-1._wp) * 999._wp
+! jc_alt                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) 
+                  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._wp) * 999._wp
+! jc_alt                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk)
                END DO
             END DO
@@ -684,5 +694,6 @@
                N_in = 0
                DO jk=k1,k2
-                  IF (tabres(ji,jj,jk,2) < -900) EXIT
+                  IF (tabres(ji,jj,jk,2) < -900._wp) EXIT
+! jc_alt                  IF (tabres(ji,jj,jk,2) == 0) EXIT
                   N_in = N_in + 1
                   tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
Index: /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_user.F90
===================================================================
--- /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_user.F90	(revision 11740)
+++ /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_user.F90	(revision 11741)
@@ -97,9 +97,14 @@
       USE oce
       USE lib_mpp
-      !
-      IMPLICIT NONE
-      !
+      USE lbclnk
+      !
+      IMPLICIT NONE
+      !
+      INTEGER :: ji, jj
       LOGICAL :: check_namelist
       CHARACTER(len=15) :: cl_check1, cl_check2, cl_check3, cl_check4 
+#if defined key_vertical
+      REAL(wp), DIMENSION(jpi,jpj) ::   zk   ! workspace
+#endif
       !!----------------------------------------------------------------------
 
@@ -110,4 +115,32 @@
       ! 2. First interpolations of potentially non zero fields
       !-------------------------------------------------------
+
+#if defined key_vertical
+      ! Build consistent parent bathymetry and number of levels
+      ! on the child grid 
+      Agrif_UseSpecialValue = .FALSE.
+      ht0_parent(:,:) = 0._wp
+      mbkt_parent(:,:) = 0
+      !
+      CALL Agrif_Bc_variable(ht0_id ,calledweight=1.,procname=interpht0 )
+      CALL Agrif_Bc_variable(mbkt_id,calledweight=1.,procname=interpmbkt)
+      !
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+            hu0_parent(ji,jj)  = MIN( ht0_parent(ji,jj), ht0_parent(ji+1,jj) ) 
+            hv0_parent(ji,jj)  = MIN( ht0_parent(ji,jj), ht0_parent(ji,jj+1) ) 
+            mbku_parent(ji,jj) = MIN(  mbkt_parent(ji+1,jj  ) , mbkt_parent(ji,jj)  )
+            mbkv_parent(ji,jj) = MIN(  mbkt_parent(ji  ,jj+1) , mbkt_parent(ji,jj)  )
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( 'Agrif_InitValues_cont', hu0_parent, 'U', 1. )
+      CALL lbc_lnk( 'Agrif_InitValues_cont', hv0_parent, 'V', 1. )
+      zk(:,:) = REAL( mbku_parent(:,:), wp )   ;   CALL lbc_lnk( 'Agrif_InitValues_cont', zk, 'U', 1. )
+      mbku_parent(:,:) = MAX( NINT( zk(:,:) ), 1 )
+      zk(:,:) = REAL( mbkv_parent(:,:), wp )   ;   CALL lbc_lnk( 'Agrif_InitValues_cont', zk, 'V', 1. )
+      mbkv_parent(:,:) = MAX( NINT( zk(:,:) ), 1 )      
+#endif
+
       Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = .TRUE.
@@ -225,4 +258,10 @@
             END IF
          ENDIF
+
+#if defined key_vertical
+	 IF ( Agrif_Parent(jpk).GT.jpk ) THEN
+	    CALL ctl_stop( ' With key_vertical, child grids must have jpk greater or equal to the parent value' )
+         ENDIF
+#endif
          !
       ENDIF
@@ -275,4 +314,8 @@
       CALL agrif_declare_variable((/1,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),umsk_id)
       CALL agrif_declare_variable((/2,1,0/),(/ind3,ind2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vmsk_id)
+# if defined key_vertical
+      CALL agrif_declare_variable((/2,2/),(/ind3,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),mbkt_id)
+      CALL agrif_declare_variable((/2,2/),(/ind3,ind3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),ht0_id)
+# endif
 
       CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,3/),scales_t_id)
@@ -311,5 +354,12 @@
       CALL Agrif_Set_bcinterp(ub2b_interp_id,interp1=Agrif_linear,interp2=AGRIF_ppm)
       CALL Agrif_Set_bcinterp(vb2b_interp_id,interp1=AGRIF_ppm,interp2=Agrif_linear)
-
+!
+! > Divergence conserving alternative:
+!      CALL Agrif_Set_bcinterp(sshn_id,interp=AGRIF_constant)
+!      CALL Agrif_Set_bcinterp(unb_id,interp1=Agrif_linear,interp2=AGRIF_constant)
+!      CALL Agrif_Set_bcinterp(vnb_id,interp1=AGRIF_constant,interp2=Agrif_linear)
+!      CALL Agrif_Set_bcinterp(ub2b_interp_id,interp1=Agrif_linear,interp2=AGRIF_constant)
+!      CALL Agrif_Set_bcinterp(vb2b_interp_id,interp1=AGRIF_constant,interp2=Agrif_linear)
+!<
 
       CALL Agrif_Set_bcinterp(un_sponge_id,interp1=Agrif_linear,interp2=AGRIF_ppm)
@@ -319,4 +369,9 @@
       CALL Agrif_Set_bcinterp(umsk_id,interp=AGRIF_constant)
       CALL Agrif_Set_bcinterp(vmsk_id,interp=AGRIF_constant)
+
+# if defined key_vertical
+      CALL Agrif_Set_bcinterp(mbkt_id,interp=AGRIF_constant)
+      CALL Agrif_Set_bcinterp(ht0_id ,interp=AGRIF_constant)
+# endif
 
       IF( ln_zdftke.OR.ln_zdfgls )   CALL Agrif_Set_bcinterp( avm_id, interp=AGRIF_linear )
@@ -341,5 +396,9 @@
       CALL Agrif_Set_bc( umsk_id, (/0,0/) )
       CALL Agrif_Set_bc( vmsk_id, (/0,0/) )
-
+# if defined key_vertical
+      ! extend the interpolation zone by 1 more point than necessary:
+      CALL Agrif_Set_bc(  mbkt_id, (/-nn_sponge_len*Agrif_irhox()-2,ind1/) )
+      CALL Agrif_Set_bc(  ht0_id,  (/-nn_sponge_len*Agrif_irhox()-2,ind1/) )
+# endif
 
       IF( ln_zdftke.OR.ln_zdfgls )   CALL Agrif_Set_bc( avm_id, (/0,ind1/) )
Index: /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/vremap.F90
===================================================================
--- /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/vremap.F90	(revision 11740)
+++ /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/vremap.F90	(revision 11741)
@@ -317,10 +317,14 @@
 
       ! specify methods
+!      opts%edge_meth = p1e_method     ! 1st-order edge interp.
+!      opts%cell_meth = plm_method     ! PLM method in cells
       opts%edge_meth = p3e_method     ! 3rd-order edge interp.
-      opts%cell_meth = ppm_method     ! PPM method in cells
-      opts%cell_lims = null_limit     ! no lim.     
+      opts%cell_meth = ppm_method     ! PPM method in cells    
 !      opts%edge_meth = p5e_method     ! 5th-order edge interp.
 !      opts%cell_meth = pqm_method     ! PQM method in cells
-!      opts%cell_lims = mono_limit     ! monotone limiter   
+
+      ! limiter
+!      opts%cell_lims = null_limit     ! no lim.
+      opts%cell_lims = mono_limit     ! monotone limiter   
  
       ! set boundary conditions
