Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/step_c1d.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/step_c1d.F90	(revision 10977)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/step_c1d.F90	(revision 10978)
@@ -79,8 +79,8 @@
                          CALL zdf_phy( kstp )         ! vertical physics update (bfr, avt, avs, avm + MLD)
 
-      IF(.NOT.ln_linssh )   CALL ssh_nxt       ( kstp, Nbb, Nnn )  ! after ssh (includes call to div_hor)
-      IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
+      IF(.NOT.ln_linssh )   CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )  ! after ssh (includes call to div_hor)
+      IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )  ! after vertical scale factors 
 
-      IF(.NOT.ln_linssh )   CALL wzv           ( kstp )  ! now cross-level velocity 
+      IF(.NOT.ln_linssh )   CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )  ! now cross-level velocity 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! diagnostics and outputs             (ua, va, ta, sa used as workspace)
@@ -125,7 +125,7 @@
                         CALL dyn_zdf    ( kstp )   ! vertical diffusion
                         CALL dyn_nxt    ( kstp )   ! lateral velocity at next time step
-      IF(.NOT.ln_linssh)CALL ssh_swp    ( kstp )   ! swap of sea surface height
+      IF(.NOT.ln_linssh)CALL ssh_swp    ( kstp, Nbb, Nnn, Naa )   ! swap of sea surface height
 
-      IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp )! swap of vertical scale factors
+      IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa )! swap of vertical scale factors
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/domain.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/domain.F90	(revision 10977)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/domain.F90	(revision 10978)
@@ -58,5 +58,5 @@
 CONTAINS
 
-   SUBROUTINE dom_init(cdstr)
+   SUBROUTINE dom_init( Kbb, Kmm, Kaa, cdstr )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dom_init  ***
@@ -73,8 +73,10 @@
       !!              - 1D configuration, move Coriolis, u and v at T-point
       !!----------------------------------------------------------------------
+      INTEGER          , INTENT(in) :: Kbb, Kmm, Kaa          ! ocean time level indices
+      CHARACTER (len=*), INTENT(in) :: cdstr                  ! model: NEMO or SAS. Determines core restart variables
+      !
       INTEGER ::   ji, jj, jk, ik   ! dummy loop indices
       INTEGER ::   iconf = 0    ! local integers
       CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7))" 
-      CHARACTER (len=*), INTENT(IN) :: cdstr                  ! model: NEMO or SAS. Determines core restart variables
       INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level
       REAL(wp), DIMENSION(jpi,jpj) ::   z1_hu_0, z1_hv_0
@@ -161,15 +163,15 @@
       !
          !       before        !          now          !       after         !
-            gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points
-            gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          !
-                                   gde3w_n = gde3w_0   !        ---          !
+            gdept(:,:,:,Kbb) = gdept_0  ;   gdept(:,:,:,Kmm) = gdept_0   !        ---          ! depth of grid-points
+            gdepw(:,:,:,Kbb) = gdepw_0  ;   gdepw(:,:,:,Kmm) = gdepw_0   !        ---          !
+                                   gde3w = gde3w_0   !        ---          !
          !                                                                  
-              e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors
-              e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    !
-              e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
-                                     e3f_n =   e3f_0   !        ---          !
-              e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          !
-             e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          !
-             e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          !
+              e3t(:,:,:,Kbb) =   e3t_0  ;     e3t(:,:,:,Kmm) =   e3t_0   ;   e3t(:,:,:,Kaa) =  e3t_0    ! scale factors
+              e3u(:,:,:,Kbb) =   e3u_0  ;     e3u(:,:,:,Kmm) =   e3u_0   ;   e3u(:,:,:,Kaa) =  e3u_0    !
+              e3v(:,:,:,Kbb) =   e3v_0  ;     e3v(:,:,:,Kmm) =   e3v_0   ;   e3v(:,:,:,Kaa) =  e3v_0    !
+                                     e3f =   e3f_0   !        ---          !
+              e3w(:,:,:,Kbb) =   e3w_0  ;     e3w(:,:,:,Kmm) =   e3w_0   !        ---          !
+             e3uw(:,:,:,Kbb) =  e3uw_0  ;    e3uw(:,:,:,Kmm) =  e3uw_0   !        ---          !
+             e3vw(:,:,:,Kbb) =  e3vw_0  ;    e3vw(:,:,:,Kmm) =  e3vw_0   !        ---          !
          !
          z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF
@@ -186,5 +188,5 @@
       ELSE                       != time varying : initialize before/now/after variables
          !
-         IF( .NOT.l_offline )  CALL dom_vvl_init 
+         IF( .NOT.l_offline )  CALL dom_vvl_init( Kbb, Kmm, Kaa )
          !
       ENDIF
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/domvvl.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/domvvl.F90	(revision 10977)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/domvvl.F90	(revision 10978)
@@ -93,5 +93,5 @@
 
 
-   SUBROUTINE dom_vvl_init
+   SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dom_vvl_init  ***
@@ -104,10 +104,10 @@
       !!
       !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
-      !!              - Regrid: e3(u/v)_n
-      !!                        e3(u/v)_b       
-      !!                        e3w_n           
-      !!                        e3(u/v)w_b      
-      !!                        e3(u/v)w_n      
-      !!                        gdept_n, gdepw_n and gde3w_n
+      !!              - Regrid: e3[u/v](:,:,:,Kmm)
+      !!                        e3[u/v](:,:,:,Kmm)       
+      !!                        e3w(:,:,:,Kmm)           
+      !!                        e3[u/v]w_b
+      !!                        e3[u/v]w_n      
+      !!                        gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w
       !!              - h(t/u/v)_0
       !!              - frq_rst_e3t and frq_rst_hdv
@@ -115,4 +115,6 @@
       !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
+      !
       INTEGER ::   ji, jj, jk
       INTEGER ::   ii0, ii1, ij0, ij1
@@ -130,33 +132,33 @@
       !
       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
-      CALL dom_vvl_rst( nit000, 'READ' )
-      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
+      CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' )
+      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
       !
       !                    !== Set of all other vertical scale factors  ==!  (now and before)
       !                                ! Horizontal interpolation of e3t
-      CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )    ! from T to U
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-      CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )    ! from T to V 
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F
+      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V 
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F
       !                                ! Vertical interpolation of e3t,u,v 
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )  ! from T to W
-      CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )  ! from U to UW
-      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )  ! from V to UW
-      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W
+      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW
+      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
+      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW
+      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
 
       ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
-      e3t_a(:,:,:) = e3t_n(:,:,:)
-      e3u_a(:,:,:) = e3u_n(:,:,:)
-      e3v_a(:,:,:) = e3v_n(:,:,:)
+      e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
+      e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
+      e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
       !
       !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
-      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration)
-      gdepw_n(:,:,1) = 0.0_wp
-      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)  ! reference to a common level z=0 for hpg
-      gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1)
-      gdepw_b(:,:,1) = 0.0_wp
+      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)       ! reference to the ocean surface (used for MLD and light penetration)
+      gdepw(:,:,1,Kmm) = 0.0_wp
+      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
+      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb)
+      gdepw(:,:,1,Kbb) = 0.0_wp
       DO jk = 2, jpk                               ! vertical sum
          DO jj = 1,jpj
@@ -165,13 +167,13 @@
                !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
                !                             ! 0.5 where jk = mikt     
-!!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
+!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ??
                zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
-               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
-                  &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
-               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
-               gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
-               gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
-                  &                + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
+               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
+               gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
+                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm)) 
+               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
+               gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb)
+               gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  &
+                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb)) 
             END DO
          END DO
@@ -179,15 +181,15 @@
       !
       !                    !==  thickness of the water column  !!   (ocean portion only)
-      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
-      hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
-      hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)
-      hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
-      hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)
+      ht_n(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
+      hu_b(:,:) = e3u(:,:,1,Kbb) * umask(:,:,1)
+      hu_n(:,:) = e3u(:,:,1,Kmm) * umask(:,:,1)
+      hv_b(:,:) = e3v(:,:,1,Kbb) * vmask(:,:,1)
+      hv_n(:,:) = e3v(:,:,1,Kmm) * vmask(:,:,1)
       DO jk = 2, jpkm1
-         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
-         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-         hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
-         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
+         ht_n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+         hu_b(:,:) = hu_b(:,:) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
+         hu_n(:,:) = hu_n(:,:) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
+         hv_b(:,:) = hv_b(:,:) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
+         hv_n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
       END DO
       !
@@ -266,5 +268,5 @@
 
 
-   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
+   SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dom_vvl_sf_nxt  ***
@@ -288,6 +290,7 @@
       !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in )           ::   kt      ! time step
-      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
+      INTEGER, INTENT( in )           ::   kt             ! time step
+      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
+      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
       !
       INTEGER                ::   ji, jj, jk            ! dummy loop indices
@@ -321,8 +324,8 @@
       !                                                ! --------------------------------------------- !
       !
-      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
+      z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
       DO jk = 1, jpkm1
-         ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
-         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+         ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0)
+         e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
       END DO
       !
@@ -337,6 +340,6 @@
          zht(:,:)   = 0._wp
          DO jk = 1, jpkm1
-            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
-            zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
+            zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
+            zht  (:,:) = zht  (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
          END DO
          zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
@@ -348,5 +351,5 @@
                DO jk = 1, jpkm1
                   hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
-                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
+                     &          * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) )
                END DO
             ENDIF
@@ -361,9 +364,9 @@
          IF( ln_vvl_ztilde ) THEN     ! z_tilde case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
             END DO
          ELSE                         ! layer case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
             END DO
          ENDIF
@@ -476,7 +479,7 @@
             zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
-         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
+         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
          DO jk = 1, jpkm1
-            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
          END DO
 
@@ -486,5 +489,5 @@
       !                                           ! ---baroclinic part--------- !
          DO jk = 1, jpkm1
-            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
+            e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
       ENDIF
@@ -501,37 +504,37 @@
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
+            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+         END DO
+         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax
          !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
+            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk)
+         END DO
+         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax
          !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
+            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk)
+         END DO
+         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kbb) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kmm) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kaa) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax
       END IF
 
@@ -540,6 +543,6 @@
       ! *********************************** !
 
-      CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
-      CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
 
       ! *********************************** !
@@ -547,9 +550,9 @@
       ! *********************************** !
 
-      hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
-      hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
+      hu_a(:,:) = e3u(:,:,1,Kaa) * umask(:,:,1)
+      hv_a(:,:) = e3v(:,:,1,Kaa) * vmask(:,:,1)
       DO jk = 2, jpkm1
-         hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
-         hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
+         hu_a(:,:) = hu_a(:,:) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
+         hv_a(:,:) = hv_a(:,:) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
       END DO
       !                                        ! Inverse of the local depth
@@ -563,5 +566,5 @@
 
 
-   SUBROUTINE dom_vvl_sf_swp( kt )
+   SUBROUTINE dom_vvl_sf_swp( kt, Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dom_vvl_sf_swp  ***
@@ -578,8 +581,8 @@
       !!               - Recompute:
       !!                    e3(u/v)_b       
-      !!                    e3w_n           
+      !!                    e3w(:,:,:,Kmm)           
       !!                    e3(u/v)w_b      
       !!                    e3(u/v)w_n      
-      !!                    gdept_n, gdepw_n  and gde3w_n
+      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
       !!                    h(u/v) and h(u/v)r
       !!
@@ -587,5 +590,6 @@
       !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! time step
+      INTEGER, INTENT( in ) ::   kt              ! time step
+      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -615,10 +619,10 @@
          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
       ENDIF
-      gdept_b(:,:,:) = gdept_n(:,:,:)
-      gdepw_b(:,:,:) = gdepw_n(:,:,:)
-
-      e3t_n(:,:,:) = e3t_a(:,:,:)
-      e3u_n(:,:,:) = e3u_a(:,:,:)
-      e3v_n(:,:,:) = e3v_a(:,:,:)
+      gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)
+      gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
+
+      e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa)
+      e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa)
+      e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa)
 
       ! Compute all missing vertical scale factor and depths
@@ -626,21 +630,21 @@
       ! Horizontal scale factor interpolations
       ! --------------------------------------
-      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
+      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are allready computed in dynnxt
       ! - JC - hu_b, hv_b, hur_b, hvr_b also
       
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  )
       
       ! Vertical scale factor interpolations
-      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
-      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
+      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
+      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
 
       ! t- and w- points depth (set the isf depth as it is in the initial step)
-      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-      gdepw_n(:,:,1) = 0.0_wp
-      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
+      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
+      gdepw(:,:,1,Kmm) = 0.0_wp
+      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
       DO jk = 2, jpk
          DO jj = 1,jpj
@@ -649,8 +653,8 @@
                                                                  ! 1 for jk = mikt
                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  &
-                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) ) 
-               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
+               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
+               gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  &
+                   &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) ) 
+               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
             END DO
          END DO
@@ -662,12 +666,12 @@
       hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
       !
-      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
+      ht_n(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
       DO jk = 2, jpkm1
-         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
+         ht_n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
       END DO
 
       ! write restart file
       ! ==================
-      IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' )
+      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
       !
       IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp')
@@ -783,5 +787,5 @@
 
 
-   SUBROUTINE dom_vvl_rst( kt, cdrw )
+   SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw )
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE dom_vvl_rst  ***
@@ -795,6 +799,7 @@
       !!                they are set to 0.
       !!----------------------------------------------------------------------
-      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
-      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      INTEGER         , INTENT(in) ::   kt        ! ocean time-step
+      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
+      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag
       !
       INTEGER ::   ji, jj, jk
@@ -806,5 +811,5 @@
          IF( ln_rstart ) THEN                   !* Read the restart file
             CALL rst_read_open                  !  open the restart file if necessary
-            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    )
+            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
             !
             id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
@@ -817,39 +822,39 @@
             !                             ! --------- !
             IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
                ! needed to restart if land processor not computed 
-               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
+               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
                WHERE ( tmask(:,:,:) == 0.0_wp ) 
-                  e3t_n(:,:,:) = e3t_0(:,:,:)
-                  e3t_b(:,:,:) = e3t_0(:,:,:)
+                  e3t(:,:,:,Kmm) = e3t_0(:,:,:)
+                  e3t(:,:,:,Kbb) = e3t_0(:,:,:)
                END WHERE
                IF( neuler == 0 ) THEN
-                  e3t_b(:,:,:) = e3t_n(:,:,:)
+                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
                ENDIF
             ELSE IF( id1 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
                IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
-               e3t_n(:,:,:) = e3t_b(:,:,:)
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
+               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
                neuler = 0
             ELSE IF( id2 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
                IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
-               e3t_b(:,:,:) = e3t_n(:,:,:)
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
+               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
                neuler = 0
             ELSE
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
                IF(lwp) write(numout,*) 'Compute scale factor from sshn'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
                DO jk = 1, jpk
-                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
+                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
                       &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
                       &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
                END DO
-               e3t_b(:,:,:) = e3t_n(:,:,:)
+               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
                neuler = 0
             ENDIF
@@ -888,21 +893,19 @@
                IF( cn_cfg == 'wad' ) THEN
                   ! Wetting and drying test case
-                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )
-                  tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
-                  sshn (:,:)     = sshb(:,:)
-                  un   (:,:,:)   = ub  (:,:,:)
-                  vn   (:,:,:)   = vb  (:,:,:)
+                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
+                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
+                  ssh (:,:,Kmm)     = ssh(:,:,Kbb)
+                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
+                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
                ELSE
                   ! if not test case
-                  sshn(:,:) = -ssh_ref
-                  sshb(:,:) = -ssh_ref
+                  ssh(:,:,Kmm) = -ssh_ref
+                  ssh(:,:,Kbb) = -ssh_ref
 
                   DO jj = 1, jpj
                      DO ji = 1, jpi
                         IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
-
-                           sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
-                           sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
-                           ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
+                           ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
+                           ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
                         ENDIF
                      ENDDO
@@ -912,9 +915,9 @@
                ! Adjust vertical metrics for all wad
                DO jk = 1, jpk
-                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) &
+                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) &
                     &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
                     &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
                END DO
-               e3t_b(:,:,:) = e3t_n(:,:,:)
+               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
 
                DO ji = 1, jpi
@@ -930,14 +933,14 @@
                ! Just to read set ssh in fact, called latter once vertical grid
                ! is set up:
-!               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )
+!               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
 !               !
 !               DO jk=1,jpk
-!                  e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &
+!                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
 !                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
 !               END DO
-!               e3t_n(:,:,:) = e3t_b(:,:,:)
-                sshn(:,:)=0._wp
-                e3t_n(:,:,:)=e3t_0(:,:,:)
-                e3t_b(:,:,:)=e3t_0(:,:,:)
+!               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
+                ssh(:,:,Kmm)=0._wp
+                e3t(:,:,:,Kmm)=e3t_0(:,:,:)
+                e3t(:,:,:,Kbb)=e3t_0(:,:,:)
                !
             END IF           ! end of ll_wd edits
@@ -957,6 +960,6 @@
          !                                           ! all cases !
          !                                           ! --------- !
-         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios )
          !                                           ! ----------------------- !
          IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/iscplhsb.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/iscplhsb.F90	(revision 10977)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/iscplhsb.F90	(revision 10978)
@@ -40,5 +40,5 @@
 CONTAINS
 
-   SUBROUTINE iscpl_cons(ptmask_b, psmask_b, pe3t_b, pts_flx, pvol_flx, prdt_iscpl)
+   SUBROUTINE iscpl_cons( Kbb, Kmm, ptmask_b, psmask_b, pe3t_b, pts_flx, pvol_flx, prdt_iscpl )
       !!---------------------------------------------------------------------- 
       !!                   ***  ROUTINE iscpl_cons  ***
@@ -48,6 +48,7 @@
       !!                compute where the correction have to be applied
       !! 
-      !! ** Method  :   compute tsn*e3t-tsb*e3tb and e3t-e3t_b
-      !!----------------------------------------------------------------------
+      !! ** Method  :   compute tsn*e3tn-tsb*e3tb and e3tn-e3tb
+      !!----------------------------------------------------------------------
+      INTEGER                     , INTENT(in ) :: Kbb, Kmm    !! time level indices
       REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b    !! mask before
       REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b      !! scale factor before
@@ -78,9 +79,9 @@
       z1_rdtiscpl = 1._wp / prdt_iscpl 
 
-      ! mask tsn and tsb 
-      tsb(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) * ptmask_b(:,:,:)
-      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) *  tmask  (:,:,:)
-      tsb(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) * ptmask_b(:,:,:)
-      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) *  tmask  (:,:,:)
+      ! mask ts(:,:,:,:,Kmm) and ts(:,:,:,:,Kbb) 
+      ts(:,:,:,jp_tem,Kbb) = ts(:,:,:,jp_tem,Kbb) * ptmask_b(:,:,:)
+      ts(:,:,:,jp_tem,Kmm) = ts(:,:,:,jp_tem,Kmm) *  tmask  (:,:,:)
+      ts(:,:,:,jp_sal,Kbb) = ts(:,:,:,jp_sal,Kbb) * ptmask_b(:,:,:)
+      ts(:,:,:,jp_sal,Kmm) = ts(:,:,:,jp_sal,Kmm) *  tmask  (:,:,:)
 
       !==============================================================================
@@ -89,5 +90,5 @@
 
       ! 
-      zdssh(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
+      zdssh(:,:) = ssh(:,:,Kmm) * ssmask(:,:) - ssh(:,:,Kbb) * psmask_b(:,:)
       IF (.NOT. ln_linssh ) zdssh = 0.0_wp ! already included in the levels by definition
       
@@ -98,12 +99,12 @@
 
                   ! volume differences
-                  zde3t = e3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
+                  zde3t = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
 
                   ! heat diff
-                  zdtem = tsn(ji,jj,jk,jp_tem) * e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
-                        - tsb(ji,jj,jk,jp_tem) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
+                  zdtem = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm) *  tmask  (ji,jj,jk)   &
+                        - ts(ji,jj,jk,jp_tem,Kbb) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
                   ! salt diff
-                  zdsal = tsn(ji,jj,jk,jp_sal) * e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
-                        - tsb(ji,jj,jk,jp_sal) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
+                  zdsal = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm) *  tmask  (ji,jj,jk)   &
+                        - ts(ji,jj,jk,jp_sal,Kbb) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
                
                   ! shh changes
@@ -296,5 +297,5 @@
 
 
-   SUBROUTINE iscpl_div( phdivn )
+   SUBROUTINE iscpl_div( Kmm, phdivn )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE iscpl_div  ***
@@ -307,4 +308,5 @@
       !! ** Action  :   phdivn   increase by the iscpl correction term
       !!----------------------------------------------------------------------
+      INTEGER                   , INTENT(in   ) ::   Kmm      ! time level index
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
       !!
@@ -315,5 +317,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / e3t_n(ji,jj,jk)
+               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / e3t(ji,jj,jk,Kmm)
             END DO
          END DO
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/iscplrst.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/iscplrst.F90	(revision 10977)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/iscplrst.F90	(revision 10978)
@@ -40,13 +40,15 @@
 CONTAINS
 
-   SUBROUTINE iscpl_stp
+   SUBROUTINE iscpl_stp( Kbb, Kmm )
       !!---------------------------------------------------------------------- 
       !!                   ***  ROUTINE iscpl_stp  ***
       !! 
       !! ** Purpose : compute initialisation
-      !!              compute extrapolation of restart variable un, vn, tsn, sshn (wetting/drying)   
+      !!              compute extrapolation of restart variable uu(Kmm), vv(Kmm), ts(Kmm), ssh(Kmm) (wetting/drying)   
       !!              compute correction term if needed
       !! 
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kbb, Kmm   ! time level indices
+      !
       INTEGER  ::   inum0
       REAL(wp), DIMENSION(jpi,jpj)     ::   zsmask_b
@@ -69,9 +71,9 @@
       CALL iscpl_init()       ! read namelist
       !                       ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
-      CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
+      CALL iscpl_rst_interpol( Kbb, Kmm, ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
       !
       IF ( ln_hsb ) THEN      ! compute correction if conservation needed
          IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
-         CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl)
+         CALL iscpl_cons( Kbb, Kmm, ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl )
       END IF
          
@@ -92,18 +94,18 @@
       !
       !                       ! set _b and _n variables equal
-      tsb (:,:,:,:) = tsn (:,:,:,:)
-      ub  (:,:,:)   = un  (:,:,:)
-      vb  (:,:,:)   = vn  (:,:,:)
-      sshb(:,:)     = sshn(:,:)
+      ts (:,:,:,:,Kbb) = ts (:,:,:,:,Kmm)
+      uu  (:,:,:,Kbb)   = uu  (:,:,:,Kmm)
+      vv  (:,:,:,Kbb)   = vv  (:,:,:,Kmm)
+      ssh(:,:,Kbb)     = ssh(:,:,Kmm)
       !
       !                       ! set _b and _n vertical scale factor equal
-      e3t_b (:,:,:) = e3t_n (:,:,:)
-      e3u_b (:,:,:) = e3u_n (:,:,:)
-      e3v_b (:,:,:) = e3v_n (:,:,:)
-      !
-      e3uw_b (:,:,:) = e3uw_n (:,:,:)
-      e3vw_b (:,:,:) = e3vw_n (:,:,:)
-      gdept_b(:,:,:) = gdept_n(:,:,:)
-      gdepw_b(:,:,:) = gdepw_n(:,:,:)
+      e3t (:,:,:,Kbb) = e3t (:,:,:,Kmm)
+      e3u (:,:,:,Kbb) = e3u (:,:,:,Kmm)
+      e3v (:,:,:,Kbb) = e3v (:,:,:,Kmm)
+      !
+      e3uw (:,:,:,Kbb) = e3uw (:,:,:,Kmm)
+      e3vw (:,:,:,Kbb) = e3vw (:,:,:,Kmm)
+      gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)
+      gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
       hu_b   (:,:)   = hu_n   (:,:)
       hv_b   (:,:)   = hv_n   (:,:)
@@ -114,14 +116,15 @@
 
 
-   SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b)
+   SUBROUTINE iscpl_rst_interpol ( Kbb, Kmm, ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b )
       !!---------------------------------------------------------------------- 
       !!                   ***  ROUTINE iscpl_rst_interpol  ***
       !! 
-      !! ** Purpose :   compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves 
+      !! ** Purpose :   compute new ts(Kmm), uu(Kmm), vv(Kmm) and ssh(Kmm) in case of evolving geometry of ice shelves 
       !!                compute 2d fields of heat, salt and volume correction
       !! 
-      !! ** Method  :   tn, sn : extrapolation from neigbourg cells
-      !!                un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
+      !! ** Method  :   ts(Kmm) : extrapolation from neigbourg cells
+      !!                uu(Kmm), vv(Kmm) : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
       !!----------------------------------------------------------------------
+      INTEGER                     , INTENT(in ) :: Kbb, Kmm                        !! time level indices
       REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before
       REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
@@ -143,6 +146,6 @@
       !
       !                 ! mask value to be sure
-      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
-      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
+      ts(:,:,:,jp_tem,Kmm) = ts(:,:,:,jp_tem,Kmm) * ptmask_b(:,:,:)
+      ts(:,:,:,jp_sal,Kmm) = ts(:,:,:,jp_sal,Kmm) * ptmask_b(:,:,:)
       !
       !                 ! compute wmask
@@ -155,6 +158,6 @@
       !    
       !                 ! compute new ssh if we open a full water column (average of the closest neigbourgs)  
-      sshb (:,:)=sshn(:,:)
-      zssh0(:,:)=sshn(:,:)
+      ssh (:,:,Kbb)=ssh(:,:,Kmm)
+      zssh0(:,:)=ssh(:,:,Kmm)
       zsmask0(:,:) = psmask_b(:,:)
       zsmask1(:,:) = psmask_b(:,:)
@@ -167,5 +170,5 @@
                summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
                IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
-                  sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
+                  ssh(ji,jj,Kmm)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
                   &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     &
                   &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
@@ -175,9 +178,9 @@
             END DO
          END DO
-         CALL lbc_lnk_multi( 'iscplrst', sshn, 'T', 1., zsmask1, 'T', 1. )
-         zssh0   = sshn
+         CALL lbc_lnk_multi( 'iscplrst', ssh(:,:,Kmm), 'T', 1., zsmask1, 'T', 1. )
+         zssh0   = ssh(:,:,Kmm)
          zsmask0 = zsmask1
       END DO
-      sshn(:,:) = sshn(:,:) * ssmask(:,:)
+      ssh(:,:,Kmm) = ssh(:,:,Kmm) * ssmask(:,:)
 
 !=============================================================================
@@ -192,5 +195,5 @@
                DO ji=1,jpi
                   IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
-                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) /   &
+                     e3t(ji,jj,jk,Kmm) = e3t_0(ji,jj,jk) * ( 1._wp + ssh(ji,jj,Kmm) /   &
                      &   ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) )
                   ENDIF
@@ -199,36 +202,36 @@
          END DO
          !
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-         CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
+         CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
+         CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
+         CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
 
          ! Vertical scale factor interpolations
          ! ------------------------------------
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
-         CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-         CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
+         CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )
+         CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
+         CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
          
          ! t- and w- points depth
          ! ----------------------
-         gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-         gdepw_n(:,:,1) = 0.0_wp
-         gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
+         gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
+         gdepw(:,:,1,Kmm) = 0.0_wp
+         gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
          DO jj = 1,jpj
             DO ji = 1,jpi
                DO jk = 2,mikt(ji,jj)-1
-                  gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
-                  gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
+                  gdept(ji,jj,jk,Kmm) = gdept_0(ji,jj,jk)
+                  gdepw(ji,jj,jk,Kmm) = gdepw_0(ji,jj,jk)
+                  gde3w(ji,jj,jk) = gdept_0(ji,jj,jk) - ssh(ji,jj,Kmm)
                END DO
                IF (mikt(ji,jj) > 1) THEN
                   jk = mikt(ji,jj)
-                  gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
-                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn   (ji,jj)
+                  gdept(ji,jj,jk,Kmm) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w(ji,jj,jk,Kmm)
+                  gdepw(ji,jj,jk,Kmm) = gdepw_0(ji,jj,jk)
+                  gde3w(ji,jj,jk) = gdept(ji,jj,jk  ,Kmm) - ssh   (ji,jj,Kmm)
                END IF
                DO jk = mikt(ji,jj)+1, jpk
-                  gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn (ji,jj)
+                  gdept(ji,jj,jk,Kmm) = gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)
+                  gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
+                  gde3w(ji,jj,jk) = gdept(ji,jj,jk  ,Kmm) - ssh (ji,jj,Kmm)
                END DO
             END DO
@@ -239,7 +242,7 @@
          ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp
          DO jk = 1, jpkm1
-            hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-            hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
-            ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
+            hu_n(:,:) = hu_n(:,:) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
+            hv_n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
+            ht_n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
          END DO
          !                                        ! Inverse of the local depth
@@ -252,11 +255,11 @@
 ! compute velocity
 ! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
-      ub(:,:,:)=un(:,:,:)
-      vb(:,:,:)=vn(:,:,:)
+      uu(:,:,:,Kbb)=uu(:,:,:,Kmm)
+      vv(:,:,:,Kbb)=vv(:,:,:,Kmm)
       DO jk = 1,jpk
          DO jj = 1,jpj
             DO ji = 1,jpi
-               un(ji,jj,jk) = ub(ji,jj,jk)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/e3u_n(ji,jj,jk)*umask(ji,jj,jk);
-               vn(ji,jj,jk) = vb(ji,jj,jk)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/e3v_n(ji,jj,jk)*vmask(ji,jj,jk);
+               uu(ji,jj,jk,Kmm) = uu(ji,jj,jk,Kbb)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/e3u(ji,jj,jk,Kmm)*umask(ji,jj,jk);
+               vv(ji,jj,jk,Kmm) = vv(ji,jj,jk,Kbb)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/e3v(ji,jj,jk,Kmm)*vmask(ji,jj,jk);
             END DO
          END DO
@@ -265,11 +268,11 @@
 ! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
 ! compute barotropic velocity now and after 
-      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
+      ztrp(:,:,:) = uu(:,:,:,Kbb)*pe3u_b(:,:,:); 
       zbub(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
+      ztrp(:,:,:) = vv(:,:,:,Kbb)*pe3v_b(:,:,:); 
       zbvb(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = un(:,:,:)*e3u_n(:,:,:); 
+      ztrp(:,:,:) = uu(:,:,:,Kmm)*e3u(:,:,:,Kmm); 
       zbun(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:); 
+      ztrp(:,:,:) = vv(:,:,:,Kmm)*e3v(:,:,:,Kmm); 
       zbvn(:,:)   = SUM(ztrp,DIM=3)
 
@@ -278,6 +281,6 @@
       zhv1=0.0_wp ;
       DO jk  = 1,jpk
-        zhu1(:,:) = zhu1(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-        zhv1(:,:) = zhv1(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
+        zhu1(:,:) = zhu1(:,:) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
+        zhv1(:,:) = zhv1(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
       END DO
       
@@ -298,6 +301,6 @@
       ! update velocity
       DO jk  = 1,jpk
-         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
-         vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
+         uu(:,:,jk,Kmm)=(uu(:,:,jk,Kmm) - zucorr(:,:))*umask(:,:,jk)
+         vv(:,:,jk,Kmm)=(vv(:,:,jk,Kmm) - zvcorr(:,:))*vmask(:,:,jk)
       END DO
 
@@ -305,6 +308,6 @@
       ! compute temp and salt
       ! compute new tn and sn if we open a new cell
-      tsb (:,:,:,:) = tsn(:,:,:,:)
-      zts0(:,:,:,:) = tsn(:,:,:,:)
+      ts (:,:,:,:,Kbb) = ts(:,:,:,:,Kmm)
+      zts0(:,:,:,:) = ts(:,:,:,:,Kmm)
       ztmask1(:,:,:) = ptmask_b(:,:,:)
       ztmask0(:,:,:) = ptmask_b(:,:,:)
@@ -319,9 +322,9 @@
                       IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
                       !! horizontal basic extrapolation
-                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
+                         ts(ji,jj,jk,1,Kmm)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
                          &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
                          &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
                          &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
-                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
+                         ts(ji,jj,jk,2,Kmm)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
                          &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
                          &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
@@ -333,7 +336,7 @@
                          summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
                          IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
-                            tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
+                            ts(ji,jj,jk,1,Kmm)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
                             &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
-                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
+                            ts(ji,jj,jk,2,Kmm)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
                             &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
                             ztmask1(ji,jj,jk)=1._wp
@@ -344,15 +347,15 @@
           END DO
           
-          CALL lbc_lnk_multi( 'iscplrst', tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.)
+          CALL lbc_lnk_multi( 'iscplrst', ts(:,:,:,jp_tem,Kmm), 'T', 1., ts(:,:,:,jp_sal,Kmm), 'T', 1., ztmask1, 'T', 1.)
 
           ! update
-          zts0(:,:,:,:) = tsn(:,:,:,:)
+          zts0(:,:,:,:) = ts(:,:,:,:,Kmm)
           ztmask0 = ztmask1
 
       END DO
 
-      ! mask new tsn field
-      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
-      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
+      ! mask new ts(:,:,:,:,Kmm) field
+      ts(:,:,:,jp_tem,Kmm) = ts(:,:,:,jp_tem,Kmm) * tmask(:,:,:)
+      ts(:,:,:,jp_sal,Kmm) = ts(:,:,:,jp_sal,Kmm) * tmask(:,:,:)
 
       ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
@@ -365,16 +368,16 @@
                   &      (tmask(ji,jj,1)==0._wp .OR. ptmask_b(ji,jj,1)==0._wp) ) THEN
                      !compute weight
-                     zdzp1 = MAX(0._wp,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
-                     zdz   =           gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
-                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - gdepw_n(ji,jj,jk  ))
+                     zdzp1 = MAX(0._wp,gdepw(ji,jj,jk+1,Kmm) - pdepw_b(ji,jj,jk+1))
+                     zdz   =           gdepw(ji,jj,jk+1,Kmm) - pdepw_b(ji,jj,jk  ) 
+                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - gdepw(ji,jj,jk  ,Kmm))
                      IF (zdz .LT. 0._wp) THEN 
                         CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
                      END IF
-                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
-                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
-                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/e3t_n(ji,jj,jk)
-                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
-                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
-                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/e3t_n(ji,jj,jk)
+                     ts(ji,jj,jk,jp_tem,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_tem,Kbb)      &
+                        &                   + zdz  *ts(ji,jj,jk  ,jp_tem,Kbb)      &
+                        &                   + zdzm1*ts(ji,jj,jk-1,jp_tem,Kbb)      )/e3t(ji,jj,jk,Kmm)
+                     ts(ji,jj,jk,jp_sal,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_sal,Kbb)      &
+                        &                   + zdz  *ts(ji,jj,jk  ,jp_sal,Kbb)      &
+                        &                   + zdzm1*ts(ji,jj,jk-1,jp_sal,Kbb)      )/e3t(ji,jj,jk,Kmm)
                   END IF
                END DO
@@ -386,6 +389,6 @@
       ! -----------------------------------------------------------------------------------------
       ! case we open a cell but no neigbour cells available to get an estimate of T and S
-      WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp) 
-         tsn(:,:,:,2) = -99._wp  ! Special value for closed pool (checking purpose in output.init)
+      WHERE (tmask(:,:,:) == 1._wp .AND. ts(:,:,:,2,Kmm) == 0._wp) 
+         ts(:,:,:,2,Kmm) = -99._wp  ! Special value for closed pool (checking purpose in output.init)
          tmask(:,:,:) = 0._wp    ! set mask to 0 to run
          umask(:,:,:) = 0._wp
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90	(revision 10977)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90	(revision 10978)
@@ -88,5 +88,5 @@
          !                                    ! -------------------
          CALL rst_read( Kbb, Kmm )            ! Read the restart file
-         IF (ln_iscpl)       CALL iscpl_stp   ! extrapolate restart to wet and dry
+         IF (ln_iscpl)       CALL iscpl_stp( Kbb, Kmm )   ! extrapolate restart to wet and dry
          CALL day_init                        ! model calendar (using both namelist and restart infos)
          !
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/divhor.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/divhor.F90	(revision 10977)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/divhor.F90	(revision 10978)
@@ -103,5 +103,5 @@
       IF( ln_isf )   CALL sbc_isf_div( hdiv, Kmm )                     !==  ice shelf  ==!   (update hdiv field)
       !
-      IF( ln_iscpl .AND. ln_hsb )   CALL iscpl_div( hdiv ) !==  ice sheet  ==!   (update hdiv field)
+      IF( ln_iscpl .AND. ln_hsb )   CALL iscpl_div( Kmm, hdiv ) !==  ice sheet  ==!   (update hdiv field)
       !
       CALL lbc_lnk( 'divhor', hdiv, 'T', 1. )   !   (no sign change)
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90	(revision 10977)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90	(revision 10978)
@@ -54,9 +54,9 @@
 CONTAINS
 
-   SUBROUTINE ssh_nxt( kt, Kbb, Kmm )
+   SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE ssh_nxt  ***
       !!                   
-      !! ** Purpose :   compute the after ssh (ssha)
+      !! ** Purpose :   compute the after ssh (ssh(Kaa))
       !!
       !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
@@ -64,10 +64,11 @@
       !!      by the time step.
       !!
-      !! ** action  :   ssha, after sea surface height
+      !! ** action  :   ssh(:,:,Kaa), after sea surface height
       !!
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt        ! time step
-      INTEGER, INTENT(in) ::   Kbb, Kmm  ! time level index
+      INTEGER                         , INTENT(in   ) ::   kt             ! time step
+      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level index
+      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height
       ! 
       INTEGER  ::   jk            ! dummy loop indice
@@ -92,5 +93,5 @@
       !                                           !------------------------------!
       IF(ln_wd_il) THEN
-         CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)
+         CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), z2dt)
       ENDIF
 
@@ -99,5 +100,5 @@
       zhdiv(:,:) = 0._wp
       DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
-        zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
+        zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
       END DO
       !                                                ! Sea surface elevation time stepping
@@ -105,5 +106,5 @@
       ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
       ! 
-      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
+      pssh(:,:,Kaa) = (  pssh(:,:,Kbb) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
       !
 #if defined key_agrif
@@ -113,6 +114,6 @@
       IF ( .NOT.ln_dynspg_ts ) THEN
          IF( ln_bdy ) THEN
-            CALL lbc_lnk( 'sshwzv', ssha, 'T', 1. )    ! Not sure that's necessary
-            CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries
+            CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1. )    ! Not sure that's necessary
+            CALL bdy_ssh( pssh(:,:,Kaa) )             ! Duplicate sea level across open boundaries
          ENDIF
       ENDIF
@@ -121,5 +122,5 @@
       !                                           !------------------------------!
       !
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask )
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa)  - : ', mask1=tmask )
       !
       IF( ln_timing )   CALL timing_stop('ssh_nxt')
@@ -128,5 +129,5 @@
 
    
-   SUBROUTINE wzv( kt )
+   SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE wzv  ***
@@ -139,9 +140,11 @@
       !!        The boundary conditions are w=0 at the bottom (no flux) and.
       !!
-      !! ** action  :   wn      : now vertical velocity
+      !! ** action  :   pww      : now vertical velocity
       !!
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! time step
+      INTEGER                         , INTENT(in)    ::   kt             ! time step
+      INTEGER                         , INTENT(in)    ::   Kbb, Kmm, Kaa  ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pww            ! now vertical velocity
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -157,5 +160,5 @@
          IF(lwp) WRITE(numout,*) '~~~~~ '
          !
-         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
+         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
       ENDIF
       !                                           !------------------------------!
@@ -179,17 +182,17 @@
          CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
          !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
-         !                             ! Same question holds for hdivn. Perhaps just for security
+         !                             ! Same question holds for hdiv. Perhaps just for security
          DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
             ! computation of w
-            wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)    &
-               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )     ) * tmask(:,:,jk)
-         END DO
-         !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
+            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk)    &
+               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )     ) * tmask(:,:,jk)
+         END DO
+         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
          DEALLOCATE( zhdiv ) 
       ELSE   ! z_star and linear free surface cases
          DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
             ! computation of w
-            wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk)                 &
-               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )  ) * tmask(:,:,jk)
+            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                 &
+               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk)
          END DO
       ENDIF
@@ -197,5 +200,5 @@
       IF( ln_bdy ) THEN
          DO jk = 1, jpkm1
-            wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
+            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
          END DO
       ENDIF
@@ -203,8 +206,8 @@
 #if defined key_agrif 
       IF( .NOT. AGRIF_Root() ) THEN 
-         IF ((nbondi ==  1).OR.(nbondi == 2)) wn(nlci-1 , :     ,:) = 0.e0      ! east 
-         IF ((nbondi == -1).OR.(nbondi == 2)) wn(2      , :     ,:) = 0.e0      ! west 
-         IF ((nbondj ==  1).OR.(nbondj == 2)) wn(:      ,nlcj-1 ,:) = 0.e0      ! north 
-         IF ((nbondj == -1).OR.(nbondj == 2)) wn(:      ,2      ,:) = 0.e0      ! south 
+         IF ((nbondi ==  1).OR.(nbondi == 2)) pww(nlci-1 , :     ,:) = 0.e0      ! east 
+         IF ((nbondi == -1).OR.(nbondi == 2)) pww(2      , :     ,:) = 0.e0      ! west 
+         IF ((nbondj ==  1).OR.(nbondj == 2)) pww(:      ,nlcj-1 ,:) = 0.e0      ! north 
+         IF ((nbondj == -1).OR.(nbondj == 2)) pww(:      ,2      ,:) = 0.e0      ! south 
       ENDIF 
 #endif 
@@ -215,5 +218,5 @@
 
 
-   SUBROUTINE ssh_swp( kt )
+   SUBROUTINE ssh_swp( kt, Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ssh_nxt  ***
@@ -221,18 +224,19 @@
       !! ** Purpose :   achieve the sea surface  height time stepping by 
       !!              applying Asselin time filter and swapping the arrays
-      !!              ssha  already computed in ssh_nxt  
+      !!              ssh(:,:,Kaa)  already computed in ssh_nxt  
       !!
       !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
       !!              from the filter, see Leclair and Madec 2010) and swap :
-      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
+      !!                ssh(:,:,Kmm) = ssh(:,:,Kaa) + atfp * ( ssh(:,:,Kbb) -2 ssh(:,:,Kmm) + ssh(:,:,Kaa) )
       !!                            - atfp * rdt * ( emp_b - emp ) / rau0
-      !!                sshn = ssha
-      !!
-      !! ** action  : - sshb, sshn   : before & now sea surface height
+      !!                ssh(:,:,Kmm) = ssh(:,:,Kaa)
+      !!
+      !! ** action  : - ssh(:,:,Kbb), ssh(:,:,Kmm)   : before & now sea surface height
       !!                               ready for the next time step
       !!
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt              ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! ocean time-step index
       !
       REAL(wp) ::   zcoef   ! local scalar
@@ -248,19 +252,19 @@
       !              !==  Euler time-stepping: no filter, just swap  ==!
       IF ( neuler == 0 .AND. kt == nit000 ) THEN
-         sshn(:,:) = ssha(:,:)                              ! now    <-- after  (before already = now)
+         ssh(:,:,Kmm) = ssh(:,:,Kaa)                              ! now    <-- after  (before already = now)
          !
       ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==!
          !                                                  ! before <-- now filtered
-         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )
+         ssh(:,:,Kbb) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa) )
          IF( .NOT.ln_linssh ) THEN                          ! before <-- with forcing removed
             zcoef = atfp * rdt * r1_rau0
-            sshb(:,:) = sshb(:,:) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
+            ssh(:,:,Kbb) = ssh(:,:,Kbb) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
                &                             -    rnf_b(:,:) + rnf   (:,:)   &
                &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:)
          ENDIF
-         sshn(:,:) = ssha(:,:)                              ! now <-- after
-      ENDIF
-      !
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask )
+         ssh(:,:,Kmm) = ssh(:,:,Kaa)                              ! now <-- after
+      ENDIF
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssh(:,:,Kbb), clinfo1=' ssh(:,:,Kbb)  - : ', mask1=tmask )
       !
       IF( ln_timing )   CALL timing_stop('ssh_swp')
@@ -268,5 +272,5 @@
    END SUBROUTINE ssh_swp
 
-   SUBROUTINE wAimp( kt )
+   SUBROUTINE wAimp( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE wAimp  ***
@@ -277,5 +281,5 @@
       !! ** Method  : - 
       !!
-      !! ** action  :   wn      : now vertical velocity (to be handled explicitly)
+      !! ** action  :   ww      : now vertical velocity (to be handled explicitly)
       !!            :   wi      : now vertical velocity (for implicit treatment)
       !!
@@ -283,4 +287,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! time step
+      INTEGER, INTENT(in) ::   Kmm  ! time level index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -305,11 +310,11 @@
          DO jj = 2, jpjm1
             DO ji = 2, fs_jpim1   ! vector opt.
-               z1_e3w = 1._wp / e3w_n(ji,jj,jk)
-               Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    &
-                  &                      + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
-                  &                          MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
+               z1_e3w = 1._wp / e3w(ji,jj,jk,Kmm)
+               Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )    &
+                  &                      + ( MAX( e2u(ji  ,jj)*e3uw(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   &
+                  &                          MIN( e2u(ji-1,jj)*e3uw(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   &
                   &                        * r1_e1e2t(ji,jj)                                                  &
-                  &                      + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
-                  &                          MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
+                  &                      + ( MAX( e1v(ji,jj  )*e3vw(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm), 0._wp ) -   &
+                  &                          MIN( e1v(ji,jj-1)*e3vw(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) )   &
                   &                        * r1_e1e2t(ji,jj)                                                  &
                   &                      ) * z1_e3w
@@ -338,6 +343,6 @@
                   zcff = MIN(1._wp, zcff)
                   !
-                  wi(ji,jj,jk) =           zcff   * wn(ji,jj,jk)
-                  wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk)
+                  wi(ji,jj,jk) =           zcff   * ww(ji,jj,jk)
+                  ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
                   !
                   Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient
@@ -351,5 +356,5 @@
       CALL iom_put("wimp",wi) 
       CALL iom_put("wi_cff",Cu_adv)
-      CALL iom_put("wexp",wn)
+      CALL iom_put("wexp",ww)
       !
       IF( ln_timing )   CALL timing_stop('wAimp')
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90	(revision 10977)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90	(revision 10978)
@@ -420,6 +420,6 @@
       IF( lk_c1d       )   CALL     c1d_init        ! 1D column configuration
                            CALL     wad_init        ! Wetting and drying options
-                           CALL     dom_init("OPA") ! Domain
-      IF( ln_crs       )   CALL     crs_init( Nnn ) ! coarsened grid: domain initialization 
+                           CALL     dom_init( Nbb, Nnn, Naa, "OPA") ! Domain
+      IF( ln_crs       )   CALL     crs_init(      Nnn )       ! coarsened grid: domain initialization 
       IF( ln_ctl       )   CALL prt_ctl_init        ! Print control
       
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10977)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10978)
@@ -165,8 +165,8 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-                            CALL ssh_nxt       ( kstp, Nbb, Nnn )  ! after ssh (includes call to div_hor)
-      IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
-                            CALL wzv           ( kstp )  ! now cross-level velocity 
-      IF( ln_zad_Aimp )     CALL wAimp         ( kstp )  ! Adaptive-implicit vertical advection partitioning
+                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
+      IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors 
+                            CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )    ! now cross-level velocity 
+      IF( ln_zad_Aimp )     CALL wAimp         ( kstp,      Nnn           )  ! Adaptive-implicit vertical advection partitioning
                             CALL eos    ( tsn, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation
                             
@@ -202,8 +202,8 @@
                                                       ! With split-explicit free surface, since now transports have been updated and ssha as well
       IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated
-                            CALL div_hor    ( kstp, Nbb, Nnn )    ! Horizontal divergence  (2nd call in time-split case)
-         IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component)
-                            CALL wzv        ( kstp )              ! now cross-level velocity 
-         IF( ln_zad_Aimp )  CALL wAimp      ( kstp )  ! Adaptive-implicit vertical advection partitioning
+                            CALL div_hor       ( kstp, Nbb, Nnn )    ! Horizontal divergence  (2nd call in time-split case)
+         IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn, Naa, kcall=2 )  ! after vertical scale factors (update depth average component)
+                            CALL wzv        ( kstp, Nbb, Nnn, ww, Naa )  ! now cross-level velocity 
+         IF( ln_zad_Aimp )  CALL wAimp      ( kstp,      Nnn )           ! Adaptive-implicit vertical advection partitioning
       ENDIF
       
@@ -283,6 +283,6 @@
                          CALL tra_nxt       ( kstp, Nbb, Nnn, Nrhs, Naa )  ! finalize (bcs) tracer fields at next time step and swap
                          CALL dyn_nxt       ( kstp, Nbb, Nnn, Naa  )  ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
-                         CALL ssh_swp       ( kstp )  ! swap of sea surface height
-      IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
+                         CALL ssh_swp       ( kstp, Nbb, Nnn, Naa )  ! swap of sea surface height
+      IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa )  ! swap of vertical scale factors
       !
       IF( ln_diahsb  )   CALL dia_hsb       ( kstp, Nbb, Nnn )  ! - ML - global conservation diagnostics
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90	(revision 10977)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90	(revision 10978)
@@ -312,5 +312,5 @@
                            CALL     eos_init        ! Equation of state
       IF( lk_c1d       )   CALL     c1d_init        ! 1D column configuration
-                           CALL     dom_init("OPA") ! Domain
+                           CALL     dom_init( Nbb, Nnn, Naa, "OPA") ! Domain
       IF( ln_ctl       )   CALL prt_ctl_init        ! Print control
 
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAO/nemogcm.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAO/nemogcm.F90	(revision 10977)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAO/nemogcm.F90	(revision 10978)
@@ -246,5 +246,5 @@
                            CALL phy_cst            ! Physical constants
                            CALL eos_init           ! Equation of state
-                           CALL dom_init('SAO')    ! Domain
+                           CALL dom_init( Nbb, Nnn, Naa, 'SAO')    ! Domain
 
 
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/nemogcm.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/nemogcm.F90	(revision 10977)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/nemogcm.F90	(revision 10978)
@@ -369,5 +369,5 @@
                            CALL phy_cst         ! Physical constants
                            CALL eos_init        ! Equation of seawater
-                           CALL dom_init('SAS') ! Domain
+                           CALL dom_init( Nbb, Nnn, Naa, 'SAS') ! Domain
       IF( ln_ctl      )    CALL prt_ctl_init    ! Print control
       
