Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/dom_oce.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/dom_oce.F90	(revision 10740)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/dom_oce.F90	(revision 10743)
@@ -62,4 +62,9 @@
    INTEGER , PUBLIC ::   neuler         !: restart euler forward option (0=Euler)
    REAL(wp), PUBLIC ::   r2dt           !: = 2*rdt except at nit000 (=rdt) if neuler=0
+
+   !!----------------------------------------------------------------------
+   !! time level indices
+   !!----------------------------------------------------------------------
+   INTEGER, PUBLIC :: Nm1, Nnn, Np1, Nrhs, Nm1_2lev, Nnn_2lev
 
    !!----------------------------------------------------------------------
@@ -129,12 +134,15 @@
    LOGICAL, PUBLIC ::   ln_sco       !: s-coordinate or hybrid z-s coordinate
    LOGICAL, PUBLIC ::   ln_isfcav    !: presence of ISF 
-   !                                                        !  ref.   ! before  !   now   ! after  !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3w_0 ,   e3w_b ,   e3w_n            !: w- vert. scale factor [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3uw_0 ,  e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0 ,  e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+   !                                                        !  ref. scale factors
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3t_0  !: t- vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3u_0  !: u- vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3v_0  !: v- vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3f_0  !: f- vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::     e3w_0  !: w- vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3uw_0  !: uw-vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::    e3vw_0  !: vw-vert. scale factor [m]
+   !                                                        !  time-dependent scale factors
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:,:) ::   e3t, e3u, e3v, e3w, e3uw, e3vw  !: vert. scale factor [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:,:) ::     e3f                             !: F-point vert. scale factor [m]
 
    !                                                        !  ref.   ! before  !   now   !
@@ -149,4 +157,14 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
+
+   !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::     e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::     e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::     e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::               e3f_n            !: f- vert. scale factor [m]
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::     e3w_b ,   e3w_n            !: w- vert. scale factor [m]
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::    e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::    e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
+   !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY
 
    INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1)
@@ -260,12 +278,9 @@
          &      gdept_n(jpi,jpj,jpk) , gdepw_n(jpi,jpj,jpk) , gde3w_n(jpi,jpj,jpk) , STAT=ierr(4) )
          !
-      ALLOCATE( e3t_0(jpi,jpj,jpk) , e3u_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) ,   &
-         &      e3t_b(jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b(jpi,jpj,jpk) ,                      e3w_b(jpi,jpj,jpk) ,   & 
-         &      e3t_n(jpi,jpj,jpk) , e3u_n(jpi,jpj,jpk) , e3v_n(jpi,jpj,jpk) , e3f_n(jpi,jpj,jpk) , e3w_n(jpi,jpj,jpk) ,   & 
-         &      e3t_a(jpi,jpj,jpk) , e3u_a(jpi,jpj,jpk) , e3v_a(jpi,jpj,jpk) ,                                             &
-         !                                                          !
-         &      e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) ,         &
-         &      e3uw_b(jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) ,         &               
-         &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,     STAT=ierr(5) )                       
+      ALLOCATE( e3t_0(jpi,jpj,jpk)       , e3u_0(jpi,jpj,jpk)       , e3v_0(jpi,jpj,jpk)       , e3f_0(jpi,jpj,jpk) ,    &
+         &      e3w_0(jpi,jpj,jpk)       , e3uw_0(jpi,jpj,jpk)      , e3vw_0(jpi,jpj,jpk)      ,                         &
+         &      e3t  (jpi,jpj,jpk,jpt)   , e3u  (jpi,jpj,jpk,jpt)   , e3v  (jpi,jpj,jpk,jpt)   ,                         & 
+         &      e3w  (jpi,jpj,jpk,jpt-1) , e3uw (jpi,jpj,jpk,jpt-1) , e3vw (jpi,jpj,jpk,jpt-1) ,                         & 
+         &      e3f  (jpi,jpj,jpk),     STAT=ierr(5) )                       
          !
       ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,                                           &
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 10740)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/domvvl.F90	(revision 10743)
@@ -563,5 +563,5 @@
 
 
-   SUBROUTINE dom_vvl_sf_swp( kt )
+   SUBROUTINE dom_vvl_sf_swp( kt, kNm1, kNnn, kNm1_2lev, kNnn_2lev )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dom_vvl_sf_swp  ***
@@ -588,4 +588,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! time step
+      INTEGER, INTENT( in ) ::   kNm1, kNnn, kNm1_2lev, kNnn_2lev ! time level indices
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -605,5 +606,5 @@
       ! Time filter and swap of scale factors
       ! =====================================
-      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
+      ! - ML - e3(t/u/v)_b are already computed in dynnxt.
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
          IF( neuler == 0 .AND. kt == nit000 ) THEN
@@ -618,7 +619,9 @@
       gdepw_b(:,:,:) = gdepw_n(:,:,:)
 
-      e3t_n(:,:,:) = e3t_a(:,:,:)
-      e3u_n(:,:,:) = e3u_a(:,:,:)
-      e3v_n(:,:,:) = e3v_a(:,:,:)
+!! TO BE DELETED
+!!$      e3t_n(:,:,:) = e3t_a(:,:,:)
+!!$      e3u_n(:,:,:) = e3u_a(:,:,:)
+!!$      e3v_n(:,:,:) = e3v_a(:,:,:)
+!! TO BE DELETED
 
       ! Compute all missing vertical scale factor and depths
@@ -629,16 +632,16 @@
       ! - JC - hu_b, hv_b, hur_b, hvr_b also
       
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,kNnn), 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(:,:,:,kNnn),  e3w(:,:,:,kNnn_2lev), 'W'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,kNnn), e3uw(:,:,:,kNnn_2lev), 'UW' )
+      CALL dom_vvl_interpol( e3v(:,:,:,kNnn), e3vw(:,:,:,kNnn_2lev), 'VW' )
+      CALL dom_vvl_interpol( e3t(:,:,:,kNm1),  e3w(:,:,:,kNm1_2lev), 'W'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,kNm1), e3uw(:,:,:,kNm1_2lev), 'UW' )
+      CALL dom_vvl_interpol( e3v(:,:,:,kNm1), e3vw(:,:,:,kNm1_2lev), '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)
+      gdept_n(:,:,1) = 0.5_wp * e3w(:,:,1,kNnn_2lev)
       gdepw_n(:,:,1) = 0.0_wp
       gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
@@ -649,7 +652,7 @@
                                                                  ! 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) ) 
+               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t(ji,jj,jk-1,kNnn)
+               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w(ji,jj,jk,kNnn_2lev) )  &
+                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w(ji,jj,jk,kNnn_2lev) ) 
                gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
             END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90	(revision 10740)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90	(revision 10743)
@@ -64,5 +64,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_nxt ( kt )
+   SUBROUTINE dyn_nxt ( kt, kNm1, kNp1, puu, pvv, pe3t, pe3u, pe3v )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_nxt  ***
@@ -92,12 +92,15 @@
       !!               un,vn   now horizontal velocity of next time-step
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
+      INTEGER, INTENT( in ) ::   kNm1, kNp1 ! before and after time level indices
+      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: puu, pvv ! now velocities to be time filtered
+      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: pe3t, pe3u, pe3v ! now scale factors to be time filtered
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   ikt          ! local integers
-      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
-      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
+      REAL(wp) ::   zue3a, zue3n, zue3b, zcoef    ! local scalars
+      REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      -
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva 
       !!----------------------------------------------------------------------
       !
@@ -115,13 +118,13 @@
          ! Ensure below that barotropic velocities match time splitting estimate
          ! Compute actual transport and replace it with ts estimate at "after" time step
-         zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
-         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
+         zue(:,:) = e3u(:,:,1,kNp1) * uu(:,:,1,kNp1) * umask(:,:,1)
+         zve(:,:) = e3v(:,:,1,kNp1) * vv(:,:,1,kNp1) * vmask(:,:,1)
          DO jk = 2, jpkm1
-            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
-            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
+            zue(:,:) = zue(:,:) + e3u(:,:,jk,kNp1) * uu(:,:,jk,kNp1) * umask(:,:,jk)
+            zve(:,:) = zve(:,:) + e3v(:,:,jk,kNp1) * vv(:,:,jk,kNp1) * vmask(:,:,jk)
          END DO
          DO jk = 1, jpkm1
-            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
-            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
+            uu(:,:,jk,kNp1) = ( uu(:,:,jk,kNp1) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
+            vv(:,:,jk,kNp1) = ( vv(:,:,jk,kNp1) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
          END DO
          !
@@ -132,6 +135,6 @@
             ! so that asselin contribution is removed at the same time 
             DO jk = 1, jpkm1
-               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
-               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk)
+               puu(:,:,jk) = ( puu(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
+               pvv(:,:,jk) = ( pvv(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk)
             END DO  
          ENDIF
@@ -144,5 +147,5 @@
 # endif
       !
-      CALL lbc_lnk_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries
+      CALL lbc_lnk_multi( 'dynnxt', uu(:,:,:,kNp1), 'U', -1., vv(:,:,:,kNp1), 'V', -1. )     !* local domain boundaries
       !
       !                                !* BDY open boundaries
@@ -157,15 +160,15 @@
          !
          !                                  ! Kinetic energy and Conversion
-         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
+         IF( ln_KE_trd  )   CALL trd_dyn( uu(:,:,:,kNp1), vv(:,:,:,kNp1), jpdyn_ken, kt )
          !
          IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
-            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
-            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
+            zua(:,:,:) = ( uu(:,:,:,kNp1) - uu(:,:,:,kNm1) ) * z1_2dt
+            zva(:,:,:) = ( vv(:,:,:,kNp1) - vv(:,:,:,kNm1) ) * z1_2dt
             CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
             CALL iom_put( "vtrd_tot", zva )
          ENDIF
          !
-         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
-         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
+         zua(:,:,:) = puu(:,:,:)             ! save the now velocity before the asselin filter
+         zva(:,:,:) = pvv(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
          !                                  !  computation of the asselin filter trends)
       ENDIF
@@ -173,25 +176,27 @@
       ! Time filter and swap of dynamics arrays
       ! ------------------------------------------
-      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
-         DO jk = 1, jpkm1
-            un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua
-            vn(:,:,jk) = va(:,:,jk)
-         END DO
-         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
+!! TO BE DELETED
+!!$      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
+!!$         DO jk = 1, jpkm1
+!!$            un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua
+!!$            vn(:,:,jk) = va(:,:,jk)
+!!$         END DO 
+!!$         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
 !!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a  
-            DO jk = 1, jpkm1
+!!$            DO jk = 1, jpkm1
 !               e3t_b(:,:,jk) = e3t_n(:,:,jk)
 !               e3u_b(:,:,jk) = e3u_n(:,:,jk)
 !               e3v_b(:,:,jk) = e3v_n(:,:,jk)
                !
-               e3t_n(:,:,jk) = e3t_a(:,:,jk)
-               e3u_n(:,:,jk) = e3u_a(:,:,jk)
-               e3v_n(:,:,jk) = e3v_a(:,:,jk)
-            END DO
+!!$               e3t_n(:,:,jk) = e3t_a(:,:,jk)
+!!$               e3u_n(:,:,jk) = e3u_a(:,:,jk)
+!!$               e3v_n(:,:,jk) = e3v_a(:,:,jk)
+!!$            END DO
 !!gm BUG end
-         ENDIF
+!!$         ENDIF
+!! TO BE DELETED
                                                             ! 
          
-      ELSE                                             !* Leap-Frog : Asselin filter and swap
+      IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin filter and swap
          !                                ! =============!
          IF( ln_linssh ) THEN             ! Fixed volume !
@@ -200,11 +205,13 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi    
-                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
-                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
+                     puu(ji,jj,jk) = puu(ji,jj,jk) + atfp * ( uu(ji,jj,jk,kNm1) - 2._wp * puu(ji,jj,jk) + uu(ji,jj,jk,kNp1) )
+                     pvv(ji,jj,jk) = pvv(ji,jj,jk) + atfp * ( vv(ji,jj,jk,kNm1) - 2._wp * pvv(ji,jj,jk) + vv(ji,jj,jk,kNp1) )
                      !
-                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
-                     vb(ji,jj,jk) = zvf
-                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
-                     vn(ji,jj,jk) = va(ji,jj,jk)
+!! TO BE DELETED
+!!$                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
+!!$                     vb(ji,jj,jk) = zvf
+!!$                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
+!!$                     vn(ji,jj,jk) = va(ji,jj,jk)
+!! TO BE DELETED
                   END DO
                END DO
@@ -213,9 +220,9 @@
          ELSE                             ! Variable volume !
             !                             ! ================!
-            ! Before scale factor at t-points
-            ! (used as a now filtered scale factor until the swap)
+            ! Time-filtered scale factor at t-points
             ! ----------------------------------------------------
+            ALLOCATE( ze3t_f(jpi,jpj,jpk) )
             DO jk = 1, jpkm1
-               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
+               ze3t_f(:,:,jk) = pe3t(:,:,jk) + atfp * ( e3t(:,:,jk,kNm1) - 2._wp * pe3t(:,:,jk) + e3t(:,:,jk,kNp1) )
             END DO
             ! Add volume filter correction: compatibility with tracer advection scheme
@@ -223,5 +230,5 @@
             zcoef = atfp * rdt * r1_rau0
 
-            e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+            ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
 
             IF ( ln_rnf ) THEN
@@ -231,6 +238,6 @@
                         DO ji = 1, jpi
                            IF( jk <=  nk_rnf(ji,jj)  ) THEN
-                               e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
-                                      &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
+                               ze3t_f(ji,jj,jk) =   ze3t_f(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
+                                      &          * ( pe3t(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
                            ENDIF
                         ENDDO
@@ -238,5 +245,5 @@
                   ENDDO
                ELSE
-                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
+                  ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
                ENDIF
             END IF
@@ -247,9 +254,9 @@
                      DO ji = 1, jpi
                         IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj)  ) THEN
-                           e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
-                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)
+                           ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
+                                &          * ( pe3t(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)
                         ELSEIF ( jk==misfkb(ji,jj) ) THEN
-                           e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
-                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk)
+                           ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
+                                &          * ( pe3t(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk)
                         ENDIF
                      END DO
@@ -258,18 +265,22 @@
             END IF
             !
+            pe3t(:,:,1:jpkm1) = ze3t_f(:,:,1:jpkm1)        ! filtered scale factor at T-points
+            !
             IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
                ! Before filtered scale factor at (u/v)-points
-               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
-               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
+               CALL dom_vvl_interpol( pe3t(:,:,:), pe3u(:,:,:), 'U' )
+               CALL dom_vvl_interpol( pe3t(:,:,:), pe3v(:,:,:), 'V' )
                DO jk = 1, jpkm1
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
-                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
+                        puu(ji,jj,jk) = puu(ji,jj,jk) + atfp * ( uu(ji,jj,jk,kNm1) - 2._wp * puu(ji,jj,jk) + uu(ji,jj,jk,kNp1) )
+                        pvv(ji,jj,jk) = pvv(ji,jj,jk) + atfp * ( vv(ji,jj,jk,kNm1) - 2._wp * pvv(ji,jj,jk) + vv(ji,jj,jk,kNp1) )
                         !
-                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
-                        vb(ji,jj,jk) = zvf
-                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
-                        vn(ji,jj,jk) = va(ji,jj,jk)
+!! TO BE DELETED
+!!$                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
+!!$                        vb(ji,jj,jk) = zvf
+!!$                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
+!!$                        vn(ji,jj,jk) = va(ji,jj,jk)
+!! TO BE DELETED
                      END DO
                   END DO
@@ -279,47 +290,50 @@
                !
                ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
-               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
-               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
-               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
+               ! Now filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
+               CALL dom_vvl_interpol( pe3t(:,:,:), ze3u_f, 'U' )
+               CALL dom_vvl_interpol( pe3t(:,:,:), ze3v_f, 'V' )
                DO jk = 1, jpkm1
                   DO jj = 1, jpj
                      DO ji = 1, jpi                  
-                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
-                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
-                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
-                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
-                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
-                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
+                        zue3a = e3u(ji,jj,jk,kNp1) * uu(ji,jj,jk,kNp1)
+                        zve3a = e3v(ji,jj,jk,kNp1) * vv(ji,jj,jk,kNp1)
+                        zue3n = pe3u(ji,jj,jk)     * puu(ji,jj,jk)
+                        zve3n = pe3v(ji,jj,jk)     * pvv(ji,jj,jk)
+                        zue3b = e3u(ji,jj,jk,kNm1) * uu(ji,jj,jk,kNm1)
+                        zve3b = e3v(ji,jj,jk,kNm1) * vv(ji,jj,jk,kNm1)
                         !
-                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
-                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
+                        puu(ji,jj,jk) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
+                        pvv(ji,jj,jk) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
                         !
-                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
-                        vb(ji,jj,jk) = zvf
-                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
-                        vn(ji,jj,jk) = va(ji,jj,jk)
+!! TO BE DELETED
+!!$                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
+!!$                        vb(ji,jj,jk) = zvf
+!!$                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
+!!$                        vn(ji,jj,jk) = va(ji,jj,jk)
+!! TO BE DELETED
                      END DO
                   END DO
                END DO
-               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
-               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
+               pe3u(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)  
+               pe3v(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
                !
                DEALLOCATE( ze3u_f , ze3v_f )
             ENDIF
             !
+            DEALLOCATE( ze3t_f )
          ENDIF
          !
          IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
-            ! Revert "before" velocities to time split estimate
+            ! Revert filtered "now" velocities to time split estimate
             ! Doing it here also means that asselin filter contribution is removed  
-            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
-            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)    
+            zue(:,:) = pe3u(:,:,1) * puu(:,:,1) * umask(:,:,1)
+            zve(:,:) = pe3v(:,:,1) * pvv(:,:,1) * vmask(:,:,1)    
             DO jk = 2, jpkm1
-               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
-               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)    
+               zue(:,:) = zue(:,:) + pe3u(:,:,jk) * puu(:,:,jk) * umask(:,:,jk)
+               zve(:,:) = zve(:,:) + pe3v(:,:,jk) * pvv(:,:,jk) * vmask(:,:,jk)    
             END DO
             DO jk = 1, jpkm1
-               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
-               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
+               puu(:,:,jk) = puu(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
+               pvv(:,:,jk) = pvv(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
             END DO
          ENDIF
@@ -331,11 +345,14 @@
       ! integration
       !
+      ! DS IMMERSE :: This is very confusing as it stands. But should 
+      ! hopefully be simpler when we do the time-level swapping for the
+      ! external mode solution. 
       !
       IF(.NOT.ln_linssh ) THEN
-         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
-         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
+         hu_b(:,:) = pe3u(:,:,1) * umask(:,:,1)
+         hv_b(:,:) = pe3v(:,:,1) * vmask(:,:,1)
          DO jk = 2, jpkm1
-            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
-            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
+            hu_b(:,:) = hu_b(:,:) + pe3u(:,:,jk) * umask(:,:,jk)
+            hv_b(:,:) = hv_b(:,:) + pe3v(:,:,jk) * vmask(:,:,jk)
          END DO
          r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
@@ -343,13 +360,13 @@
       ENDIF
       !
-      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
-      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
-      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
-      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
+      un_b(:,:) = e3u(:,:,1,kNp1) * uu(:,:,1,kNp1) * umask(:,:,1)
+      ub_b(:,:) = pe3u(:,:,1)     * puu(:,:,1) * umask(:,:,1)
+      vn_b(:,:) = e3v(:,:,1,kNp1) * vv(:,:,1,kNp1) * vmask(:,:,1)
+      vb_b(:,:) = pe3v(:,:,1)     * pvv(:,:,1) * vmask(:,:,1)
       DO jk = 2, jpkm1
-         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
-         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
-         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
-         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
+         un_b(:,:) = un_b(:,:) + e3u(:,:,jk,kNp1) * uu(:,:,jk,kNp1) * umask(:,:,jk)
+         ub_b(:,:) = ub_b(:,:) + pe3u(:,:,jk)     * puu(:,:,jk) * umask(:,:,jk)
+         vn_b(:,:) = vn_b(:,:) + e3v(:,:,jk,kNp1) * vv(:,:,jk,kNp1) * vmask(:,:,jk)
+         vb_b(:,:) = vb_b(:,:) + pe3v(:,:,jk)     * pvv(:,:,jk) * vmask(:,:,jk)
       END DO
       un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
@@ -363,11 +380,11 @@
       ENDIF
       IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
-         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
-         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
+         zua(:,:,:) = ( puu(:,:,:) - zua(:,:,:) ) * z1_2dt
+         zva(:,:,:) = ( pvv(:,:,:) - zva(:,:,:) ) * z1_2dt
          CALL trd_dyn( zua, zva, jpdyn_atf, kt )
       ENDIF
       !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
-         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=uu(:,:,:,kNp1), clinfo1=' nxt  - uu(:,:,:,kNp1): ', mask1=umask,   &
+         &                       tab3d_2=vv(:,:,:,kNp1), clinfo2=' vv(:,:,:,kNp1): '       , mask2=vmask )
       ! 
       IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
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 10740)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90	(revision 10743)
@@ -401,4 +401,11 @@
       CALL nemo_alloc()
 
+      ! Initialise time level indices
+      Nm1 = 1; Nnn = 2; Np1 = 3; Nrhs = Np1
+      Nm1_2lev = 1; Nnn_2lev = 2
+
+      ! Initialisation of temporary pointers (to be deleted after development finished)
+      CALL update_pointers()
+
       !                             !-------------------------------!
       !                             !  NEMO general initialization  !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/oce.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/oce.F90	(revision 10740)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/oce.F90	(revision 10743)
@@ -17,8 +17,7 @@
    PUBLIC oce_alloc ! routine called by nemo_init in nemogcm.F90
 
-   !! dynamics and tracer fields                            ! before ! now    ! after  ! the after trends becomes the fields
-   !! --------------------------                            ! fields ! fields ! trends ! only after tra_zdf and dyn_spg
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   ub   ,  un    , ua     !: i-horizontal velocity        [m/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity        [m/s]
+   !! dynamics and tracer fields  NOTE THAT "TARGET" ATTRIBUTE CAN BE REMOVED AFTER IMMERSE DEVELOPMENT FINISHED                            
+   !! --------------------------                            
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:), TARGET ::   uu   ,  vv     !: horizontal velocities        [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wn             !: vertical velocity            [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wi             !: vertical vel. (adaptive-implicit) [m/s]
@@ -69,4 +68,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fraqsr_1lev        !: fraction of solar net radiation absorbed in the first ocean level [-]
 
+   !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::   ub   ,  un    , ua     !: i-horizontal velocity        [m/s]
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va     !: j-horizontal velocity        [m/s]
+
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -84,6 +87,5 @@
       !
       ierr(:) = 0 
-      ALLOCATE( ub   (jpi,jpj,jpk)      , un   (jpi,jpj,jpk)      , ua(jpi,jpj,jpk)       ,     &
-         &      vb   (jpi,jpj,jpk)      , vn   (jpi,jpj,jpk)      , va(jpi,jpj,jpk)       ,     &          
+      ALLOCATE( uu   (jpi,jpj,jpk, jpt) , vv   (jpi,jpj,jpk,jpt)  ,                             &
          &      wn   (jpi,jpj,jpk)      , hdivn(jpi,jpj,jpk)      ,                             &
          &      tsb  (jpi,jpj,jpk,jpts) , tsn  (jpi,jpj,jpk,jpts) , tsa(jpi,jpj,jpk,jpts) ,     &
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/par_oce.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/par_oce.F90	(revision 10740)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/par_oce.F90	(revision 10743)
@@ -51,4 +51,5 @@
    INTEGER, PUBLIC ::   jpj   !                                                    !: second dimension
    INTEGER, PUBLIC ::   jpk   ! = jpkglo                                           !: third  dimension
+   INTEGER, PUBLIC, PARAMETER :: jpt = 3                                           !: fourth dimension (time level)
    INTEGER, PUBLIC ::   jpim1 ! = jpi-1                                            !: inner domain indices
    INTEGER, PUBLIC ::   jpjm1 ! = jpj-1                                            !:   -     -      -
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 10740)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10743)
@@ -44,4 +44,5 @@
 
    PUBLIC   stp   ! called by nemogcm.F90
+   PUBLIC   update_pointers ! called by nemo_init
 
    !!----------------------------------------------------------------------
@@ -276,7 +277,18 @@
 !!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine
                          CALL tra_nxt       ( kstp )  ! finalize (bcs) tracer fields at next time step and swap
-                         CALL dyn_nxt       ( kstp )  ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
+                         CALL dyn_nxt       ( kstp, Nm1, Np1, uu(:,:,:,Nnn), vv(:,:,:,Nnn), e3t(:,:,:,Nnn), e3u(:,:,:,Nnn), e3v(:,:,:,Nnn) )  ! 
                          CALL ssh_swp       ( kstp )  ! swap of sea surface height
-      IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
+      !
+      ! Swap time levels
+      Nrhs = Nm1
+      Nm1 = Nnn
+      Nnn = Np1
+      Np1 = Nrhs
+      !
+      ! Update temporary pointers
+      CALL update_pointers()
+      !
+      ! Note that 2-time-level indices don't need to be swapped because both "before" and "now" fields are derived in dom_vvl_sf_swp
+      IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp, Nm1, Nnn, Nm1_2lev, Nnn_2lev )  ! interpolate vertical scale factors for Nnn time level
       !
       IF( ln_diahsb  )   CALL dia_hsb       ( kstp )  ! - ML - global conservation diagnostics
@@ -329,4 +341,29 @@
    END SUBROUTINE stp
    
+   SUBROUTINE update_pointers
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE update_pointers  ***
+      !!
+      !! ** Purpose :   Associate temporary pointer arrays.
+      !!                For IMMERSE development phase only - to be deleted
+      !!
+      !! ** Method  :
+      !!----------------------------------------------------------------------
+
+      ub => uu(:,:,:,Nm1); un => uu(:,:,:,Nnn); ua => uu(:,:,:,Np1)
+      vb => vv(:,:,:,Nm1); vn => vv(:,:,:,Nnn); va => vv(:,:,:,Np1)
+
+      e3t_b => e3t(:,:,:,Nm1); e3t_n => e3t(:,:,:,Nnn); e3t_a => e3t(:,:,:,Np1)
+      e3u_b => e3u(:,:,:,Nm1); e3u_n => e3u(:,:,:,Nnn); e3u_a => e3u(:,:,:,Np1)
+      e3v_b => e3v(:,:,:,Nm1); e3v_n => e3v(:,:,:,Nnn); e3v_a => e3v(:,:,:,Np1)
+
+      e3f_n => e3f(:,:,:)
+
+      e3w_b  => e3w (:,:,:,Nm1_2lev); e3w_n  => e3w (:,:,:,Nnn_2lev)
+      e3uw_b => e3uw(:,:,:,Nm1_2lev); e3uw_n => e3uw(:,:,:,Nnn_2lev)
+      e3vw_b => e3vw(:,:,:,Nm1_2lev); e3vw_n => e3vw(:,:,:,Nnn_2lev)
+
+   END SUBROUTINE update_pointers
+
    !!======================================================================
 END MODULE step
