Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv.F90	(revision 10788)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv.F90	(revision 10789)
@@ -53,5 +53,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_adv( kt )
+   SUBROUTINE dyn_adv( kt, ktlev1, ktlev2, pu_rhs, pv_rhs )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_adv  ***
@@ -67,4 +67,6 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   ktlev1, ktlev2   ! time level indices for source terms
+      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends
       !!----------------------------------------------------------------------
       !
@@ -73,10 +75,10 @@
       SELECT CASE( n_dynadv )    !==  compute advection trend and add it to general trend  ==!
       CASE( np_VEC_c2  )     
-         CALL dyn_keg     ( kt, nn_dynkeg )    ! vector form : horizontal gradient of kinetic energy
-         CALL dyn_zad     ( kt )               ! vector form : vertical advection
+         CALL dyn_keg     ( kt, ktlev2, nn_dynkeg, pu_rhs, pv_rhs )    ! vector form : horizontal gradient of kinetic energy
+         CALL dyn_zad     ( kt, ktlev2,            pu_rhs, pv_rhs )    ! vector form : vertical advection
       CASE( np_FLX_c2  ) 
-         CALL dyn_adv_cen2( kt )               ! 2nd order centered scheme
+         CALL dyn_adv_cen2( kt, ktlev2,            pu_rhs, pv_rhs )    ! 2nd order centered scheme
       CASE( np_FLX_ubs )   
-         CALL dyn_adv_ubs ( kt )               ! 3rd order UBS      scheme (UP3)
+         CALL dyn_adv_ubs ( kt, ktlev1, ktlev2,    pu_rhs, pv_rhs  )               ! 3rd order UBS      scheme (UP3)
       END SELECT
       !
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_cen2.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_cen2.F90	(revision 10788)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_cen2.F90	(revision 10789)
@@ -35,5 +35,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_adv_cen2( kt )
+   SUBROUTINE dyn_adv_cen2( kt, ktlev, pu_rhs, pv_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_adv_cen2  ***
@@ -44,7 +44,9 @@
       !! ** Method  :   Trend evaluated using now fields (centered in time) 
       !!
-      !! ** Action  :   (ua,va) updated with the now vorticity term trend
+      !! ** Action  :   (pu_rhs,pv_rhs) updated with the now vorticity term trend
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      INTEGER, INTENT( in ) ::   ktlev   ! time level index for source terms
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -60,6 +62,6 @@
       !
       IF( l_trddyn ) THEN           ! trends: store the input trends
-         zfu_uw(:,:,:) = ua(:,:,:)
-         zfv_vw(:,:,:) = va(:,:,:)
+         zfu_uw(:,:,:) = pu_rhs(:,:,:)
+         zfv_vw(:,:,:) = pv_rhs(:,:,:)
       ENDIF
       !
@@ -67,20 +69,20 @@
       !
       DO jk = 1, jpkm1                    ! horizontal transport
-         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
-         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
+         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,ktlev) * uu(:,:,jk,ktlev)
+         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,ktlev) * vv(:,:,jk,ktlev)
          DO jj = 1, jpjm1                 ! horizontal momentum fluxes (at T- and F-point)
             DO ji = 1, fs_jpim1   ! vector opt.
-               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
-               zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) )
-               zfu_f(ji  ,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) )
-               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
+               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( uu(ji,jj,jk,ktlev) + uu(ji+1,jj  ,jk,ktlev) )
+               zfv_f(ji  ,jj  ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( uu(ji,jj,jk,ktlev) + uu(ji  ,jj+1,jk,ktlev) )
+               zfu_f(ji  ,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( vv(ji,jj,jk,ktlev) + vv(ji+1,jj  ,jk,ktlev) )
+               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( vv(ji,jj,jk,ktlev) + vv(ji  ,jj+1,jk,ktlev) )
             END DO
          END DO
          DO jj = 2, jpjm1                 ! divergence of horizontal momentum fluxes
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
-                  &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
-               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
-                  &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
+                  &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev)
+               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
+                  &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev)
             END DO
          END DO
@@ -88,9 +90,9 @@
       !
       IF( l_trddyn ) THEN           ! trends: send trend to trddyn for diagnostic
-         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
-         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
+         zfu_uw(:,:,:) = pu_rhs(:,:,:) - zfu_uw(:,:,:)
+         zfv_vw(:,:,:) = pv_rhs(:,:,:) - zfv_vw(:,:,:)
          CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
-         zfu_t(:,:,:) = ua(:,:,:)
-         zfv_t(:,:,:) = va(:,:,:)
+         zfu_t(:,:,:) = pu_rhs(:,:,:)
+         zfv_t(:,:,:) = pv_rhs(:,:,:)
       ENDIF
       !
@@ -106,6 +108,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1)
-               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1)
+               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * uu(ji,jj,1,ktlev)
+               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * vv(ji,jj,1,ktlev)
             END DO
          END DO
@@ -114,11 +116,11 @@
          DO jj = 2, jpj                       ! 1/4 * Vertical transport
             DO ji = 2, jpi
-               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
+               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
             END DO
          END DO
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
-               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
+               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj  ,jk) ) * ( uu(ji,jj,jk,ktlev) + uu(ji,jj,jk-1,ktlev) )
+               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji  ,jj+1,jk) ) * ( vv(ji,jj,jk,ktlev) + vv(ji,jj,jk-1,ktlev) )
             END DO
          END DO
@@ -127,6 +129,6 @@
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
-               va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev)
+               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev)
             END DO
          END DO
@@ -134,6 +136,6 @@
       !
       IF( l_trddyn ) THEN                 ! trends: send trend to trddyn for diagnostic
-         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
-         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
+         zfu_t(:,:,:) = pu_rhs(:,:,:) - zfu_t(:,:,:)
+         zfv_t(:,:,:) = pv_rhs(:,:,:) - zfv_t(:,:,:)
          CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
       ENDIF
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90	(revision 10788)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90	(revision 10789)
@@ -41,5 +41,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_adv_ubs( kt )
+   SUBROUTINE dyn_adv_ubs( kt, ktlev1, ktlev2, pu_rhs, pv_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_adv_ubs  ***
@@ -64,9 +64,11 @@
       !!      gamma1=1/3 and gamma2=1/32.
       !!
-      !! ** Action : - (ua,va) updated with the 3D advective momentum trends
+      !! ** Action : - (pu_rhs,pv_rhs) updated with the 3D advective momentum trends
       !!
       !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling. 
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt               ! ocean time-step index
+      INTEGER, INTENT(in) ::   ktlev1, ktlev2   ! time level indices for source terms
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -95,6 +97,6 @@
       !
       IF( l_trddyn ) THEN           ! trends: store the input trends
-         zfu_uw(:,:,:) = ua(:,:,:)
-         zfv_vw(:,:,:) = va(:,:,:)
+         zfu_uw(:,:,:) = pu_rhs(:,:,:)
+         zfv_vw(:,:,:) = pv_rhs(:,:,:)
       ENDIF
       !                                      ! =========================== !
@@ -102,15 +104,15 @@
          !                                   ! =========================== !
          !                                         ! horizontal volume fluxes
-         zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
-         zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
+         zfu(:,:,jk) = e2u(:,:) * e3u(:,:,jk,ktlev2) * uu(:,:,jk,ktlev2)
+         zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,ktlev2) * vv(:,:,jk,ktlev2)
          !            
          DO jj = 2, jpjm1                          ! laplacian
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj  ,jk) - 2.*ub (ji,jj,jk) + ub (ji-1,jj  ,jk) ) * umask(ji,jj,jk)
-               zlv_vv(ji,jj,jk,1) = ( vb (ji  ,jj+1,jk) - 2.*vb (ji,jj,jk) + vb (ji  ,jj-1,jk) ) * vmask(ji,jj,jk)
-               zlu_uv(ji,jj,jk,1) = ( ub (ji  ,jj+1,jk) - ub (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
-                  &               - ( ub (ji  ,jj  ,jk) - ub (ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
-               zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj  ,jk) - vb (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
-                  &               - ( vb (ji  ,jj  ,jk) - vb (ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
+               zlu_uu(ji,jj,jk,1) = ( uu (ji+1,jj  ,jk,ktlev1) - 2.*uu (ji,jj,jk,ktlev1) + uu (ji-1,jj  ,jk,ktlev1) ) * umask(ji,jj,jk)
+               zlv_vv(ji,jj,jk,1) = ( vv (ji  ,jj+1,jk,ktlev1) - 2.*vv (ji,jj,jk,ktlev1) + vv (ji  ,jj-1,jk,ktlev1) ) * vmask(ji,jj,jk)
+               zlu_uv(ji,jj,jk,1) = ( uu (ji  ,jj+1,jk,ktlev1) - uu (ji  ,jj  ,jk,ktlev1) ) * fmask(ji  ,jj  ,jk)   &
+                  &               - ( uu (ji  ,jj  ,jk,ktlev1) - uu (ji  ,jj-1,jk,ktlev1) ) * fmask(ji  ,jj-1,jk)
+               zlv_vu(ji,jj,jk,1) = ( vv (ji+1,jj  ,jk,ktlev1) - vv (ji  ,jj  ,jk,ktlev1) ) * fmask(ji  ,jj  ,jk)   &
+                  &               - ( vv (ji  ,jj  ,jk,ktlev1) - vv (ji-1,jj  ,jk,ktlev1) ) * fmask(ji-1,jj  ,jk)
                !
                zlu_uu(ji,jj,jk,2) = ( zfu(ji+1,jj  ,jk) - 2.*zfu(ji,jj,jk) + zfu(ji-1,jj  ,jk) ) * umask(ji,jj,jk)
@@ -132,11 +134,11 @@
       DO jk = 1, jpkm1                       ! ====================== !
          !                                         ! horizontal volume fluxes
-         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
-         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
+         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,ktlev2) * uu(:,:,jk,ktlev2)
+         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,ktlev2) * vv(:,:,jk,ktlev2)
          !
          DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
             DO ji = 1, fs_jpim1   ! vector opt.
-               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
-               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
+               zui = ( uu(ji,jj,jk,ktlev2) + uu(ji+1,jj  ,jk,ktlev2) )
+               zvj = ( vv(ji,jj,jk,ktlev2) + vv(ji  ,jj+1,jk,ktlev2) )
                !
                IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
@@ -164,24 +166,24 @@
                !
                zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
-                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
+                  &                * ( uu(ji,jj,jk,ktlev2) + uu(ji  ,jj+1,jk,ktlev2) - gamma1 * zl_u )
                zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
-                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
+                  &                * ( vv(ji,jj,jk,ktlev2) + vv(ji+1,jj  ,jk,ktlev2) - gamma1 * zl_v )
             END DO
          END DO
          DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
-                  &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
-               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
-                  &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
+                  &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev2)
+               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
+                  &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev2)
             END DO
          END DO
       END DO
       IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic
-         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
-         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
+         zfu_uw(:,:,:) = pu_rhs(:,:,:) - zfu_uw(:,:,:)
+         zfv_vw(:,:,:) = pv_rhs(:,:,:) - zfv_vw(:,:,:)
          CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
-         zfu_t(:,:,:) = ua(:,:,:)
-         zfv_t(:,:,:) = va(:,:,:)
+         zfu_t(:,:,:) = pu_rhs(:,:,:)
+         zfv_t(:,:,:) = pv_rhs(:,:,:)
       ENDIF
       !                                      ! ==================== !
@@ -199,6 +201,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1)
-               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1)
+               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * uu(ji,jj,1,ktlev2)
+               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * vv(ji,jj,1,ktlev2)
             END DO
          END DO
@@ -207,11 +209,11 @@
          DO jj = 2, jpj
             DO ji = 2, jpi
-               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
+               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
             END DO
          END DO
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
-               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
+               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( uu(ji,jj,jk,ktlev2) + uu(ji,jj,jk-1,ktlev2) )
+               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vv(ji,jj,jk,ktlev2) + vv(ji,jj,jk-1,ktlev2) )
             END DO
          END DO
@@ -220,6 +222,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
-               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
+               pu_rhs(ji,jj,jk) =  pu_rhs(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev2)
+               pv_rhs(ji,jj,jk) =  pv_rhs(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev2)
             END DO
          END DO
@@ -227,10 +229,10 @@
       !
       IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic
-         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
-         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
+         zfu_t(:,:,:) = pu_rhs(:,:,:) - zfu_t(:,:,:)
+         zfv_t(:,:,:) = pv_rhs(:,:,:) - zfv_t(:,:,:)
          CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
       ENDIF
       !                                         ! Control print
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' uu(:,:,:,ktlev1)s2 adv - Ua: ', mask1=umask,   &
          &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynkeg.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynkeg.F90	(revision 10788)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynkeg.F90	(revision 10789)
@@ -44,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_keg( kt, kscheme )
+   SUBROUTINE dyn_keg( kt, ktlev, kscheme, pu_rhs, pv_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_keg  ***
@@ -54,16 +54,16 @@
       !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that 
       !!      conserve kinetic energy. Compute the now horizontal kinetic energy 
-      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
+      !!         zhke = 1/2 [ mi-1( uu(:,:,:,ktlev)^2 ) + mj-1( vv(:,:,:,ktlev)^2 ) ]
       !!              * kscheme = nkeg_HW : Hollingsworth correction following
       !!      Arakawa (2001). The now horizontal kinetic energy is given by:
-      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((un(j+1)+un(j-1))/2)^2  )
-      !!                    + mj-1(  2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2  ) ]
+      !!         zhke = 1/6 [ mi-1(  2 * uu(:,:,:,ktlev)^2 + ((uu(j+1,ktlev)+uu(j-1,ktlev))/2)^2  )
+      !!                    + mj-1(  2 * vv(:,:,:,ktlev)^2 + ((vv(i+1,ktlev)+vv(i-1,ktlev))/2)^2  ) ]
       !!      
       !!      Take its horizontal gradient and add it to the general momentum
-      !!      trend (ua,va).
-      !!         ua = ua - 1/e1u di[ zhke ]
-      !!         va = va - 1/e2v dj[ zhke ]
+      !!      trend (pu_rhs,pv_rhs).
+      !!         pu_rhs = pu_rhs - 1/e1u di[ zhke ]
+      !!         pv_rhs = pv_rhs - 1/e2v dj[ zhke ]
       !!
-      !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
+      !! ** Action : - Update the (pu_rhs, pv_rhs) with the hor. ke gradient trend
       !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
       !!
@@ -72,5 +72,7 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
+      INTEGER, INTENT( in ) ::   ktlev     ! time level index for source terms
       INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme 
+      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends
       !
       INTEGER  ::   ji, jj, jk, jb    ! dummy loop indices
@@ -92,6 +94,6 @@
       IF( l_trddyn ) THEN           ! Save the input trends
          ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
-         ztrdu(:,:,:) = ua(:,:,:) 
-         ztrdv(:,:,:) = va(:,:,:) 
+         ztrdu(:,:,:) = pu_rhs(:,:,:) 
+         ztrdv(:,:,:) = pv_rhs(:,:,:) 
       ENDIF
       
@@ -109,5 +111,5 @@
                      ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
                      ifu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )
-                     un(ii-ifu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk)
+                     uu(ii-ifu,ij,jk,ktlev) = uu(ii,ij,jk,ktlev) * umask(ii,ij,jk)
                   END DO
                END DO
@@ -119,5 +121,5 @@
                      ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
                      ifv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )
-                     vn(ii,ij-ifv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk)
+                     vv(ii,ij-ifv,jk,ktlev) = vv(ii,ij,jk,ktlev) * vmask(ii,ij,jk)
                   END DO
                END DO
@@ -132,8 +134,8 @@
             DO jj = 2, jpj
                DO ji = fs_2, jpi   ! vector opt.
-                  zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
-                     &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
-                  zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
-                     &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
+                  zu =    uu(ji-1,jj  ,jk,ktlev) * uu(ji-1,jj  ,jk,ktlev)   &
+                     &  + uu(ji  ,jj  ,jk,ktlev) * uu(ji  ,jj  ,jk,ktlev)
+                  zv =    vv(ji  ,jj-1,jk,ktlev) * vv(ji  ,jj-1,jk,ktlev)   &
+                     &  + vv(ji  ,jj  ,jk,ktlev) * vv(ji  ,jj  ,jk,ktlev)
                   zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
                END DO  
@@ -145,13 +147,13 @@
             DO jj = 2, jpjm1       
                DO ji = fs_2, jpim1   ! vector opt.
-                  zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
-                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
-                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
-                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
+                  zu = 8._wp * ( uu(ji-1,jj  ,jk,ktlev) * uu(ji-1,jj  ,jk,ktlev)    &
+                     &         + uu(ji  ,jj  ,jk,ktlev) * uu(ji  ,jj  ,jk,ktlev) )  &
+                     &   +     ( uu(ji-1,jj-1,jk,ktlev) + uu(ji-1,jj+1,jk,ktlev) ) * ( uu(ji-1,jj-1,jk,ktlev) + uu(ji-1,jj+1,jk,ktlev) )   &
+                     &   +     ( uu(ji  ,jj-1,jk,ktlev) + uu(ji  ,jj+1,jk,ktlev) ) * ( uu(ji  ,jj-1,jk,ktlev) + uu(ji  ,jj+1,jk,ktlev) )
                      !
-                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
-                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
-                     &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
-                     &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
+                  zv = 8._wp * ( vv(ji  ,jj-1,jk,ktlev) * vv(ji  ,jj-1,jk,ktlev)    &
+                     &         + vv(ji  ,jj  ,jk,ktlev) * vv(ji  ,jj  ,jk,ktlev) )  &
+                     &  +      ( vv(ji-1,jj-1,jk,ktlev) + vv(ji+1,jj-1,jk,ktlev) ) * ( vv(ji-1,jj-1,jk,ktlev) + vv(ji+1,jj-1,jk,ktlev) )   &
+                     &  +      ( vv(ji-1,jj  ,jk,ktlev) + vv(ji+1,jj  ,jk,ktlev) ) * ( vv(ji-1,jj  ,jk,ktlev) + vv(ji+1,jj  ,jk,ktlev) )
                   zhke(ji,jj,jk) = r1_48 * ( zv + zu )
                END DO  
@@ -164,6 +166,6 @@
       IF (ln_bdy) THEN
          ! restore velocity masks at points outside boundary
-         un(:,:,:) = un(:,:,:) * umask(:,:,:)
-         vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)
+         uu(:,:,:,ktlev) = uu(:,:,:,ktlev) * umask(:,:,:)
+         vv(:,:,:,ktlev) = vv(:,:,:,ktlev) * vmask(:,:,:)
       ENDIF      
 
@@ -172,6 +174,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
-               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
+               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
             END DO 
          END DO
@@ -179,6 +181,6 @@
       !
       IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
-         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
          DEALLOCATE( ztrdu , ztrdv )
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90	(revision 10788)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90	(revision 10789)
@@ -96,15 +96,17 @@
 CONTAINS
 
-   SUBROUTINE dyn_vor( kt )
+   SUBROUTINE dyn_vor( kt, ktlev, pu_rhs, pv_rhs )
       !!----------------------------------------------------------------------
       !!
       !! ** Purpose :   compute the lateral ocean tracer physics.
       !!
-      !! ** Action : - Update (ua,va) with the now vorticity term trend
+      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
       !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
       !!               and planetary vorticity trends) and send them to trd_dyn 
       !!               for futher diagnostics (l_trddyn=T)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      INTEGER, INTENT( in ) ::   ktlev   ! time level index for source terms
+      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv
@@ -120,14 +122,14 @@
          ztrdv(:,:,:) = va(:,:,:)
          SELECT CASE( nvor_scheme )
-         CASE( np_ENS )           ;   CALL vor_ens( kt, ncor, un , vn , ua, va )   ! enstrophy conserving scheme
-            IF( ln_stcor )            CALL vor_ens( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, ncor, un , vn , ua, va )   ! energy conserving scheme
-            IF( ln_stcor )            CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         CASE( np_ENT )           ;   CALL vor_enT( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (T-pts)
-            IF( ln_stcor )            CALL vor_enT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         CASE( np_EET )           ;   CALL vor_eeT( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (een with e3t)
-            IF( ln_stcor )            CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
-         CASE( np_EEN )           ;   CALL vor_een( kt, ncor, un , vn , ua, va )   ! energy & enstrophy scheme
-            IF( ln_stcor )            CALL vor_een( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
+         CASE( np_ENS )           ;   CALL vor_ens( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! enstrophy conserving scheme
+            IF( ln_stcor )            CALL vor_ens( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend
+         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! energy conserving scheme
+            IF( ln_stcor )            CALL vor_ene( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend
+         CASE( np_ENT )           ;   CALL vor_enT( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! energy conserving scheme (T-pts)
+            IF( ln_stcor )            CALL vor_enT( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend
+         CASE( np_EET )           ;   CALL vor_eeT( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! energy conserving scheme (een with e3t)
+            IF( ln_stcor )            CALL vor_eeT( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend
+         CASE( np_EEN )           ;   CALL vor_een( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! energy & enstrophy scheme
+            IF( ln_stcor )            CALL vor_een( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend
          END SELECT
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
@@ -139,9 +141,9 @@
             ztrdv(:,:,:) = va(:,:,:)
             SELECT CASE( nvor_scheme )
-            CASE( np_ENT )           ;   CALL vor_enT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (T-pts)
-            CASE( np_EET )           ;   CALL vor_eeT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (een with e3t)
-            CASE( np_ENE )           ;   CALL vor_ene( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme
-            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, nrvm, un , vn , ua, va )  ! enstrophy conserving scheme
-            CASE( np_EEN )           ;   CALL vor_een( kt, nrvm, un , vn , ua, va )  ! energy & enstrophy scheme
+            CASE( np_ENT )           ;   CALL vor_enT( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! energy conserving scheme (T-pts)
+            CASE( np_EET )           ;   CALL vor_eeT( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! energy conserving scheme (een with e3t)
+            CASE( np_ENE )           ;   CALL vor_ene( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! energy conserving scheme
+            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! enstrophy conserving scheme
+            CASE( np_EEN )           ;   CALL vor_een( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! energy & enstrophy scheme
             END SELECT
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
@@ -156,22 +158,22 @@
          SELECT CASE ( nvor_scheme )      !==  vorticity trend added to the general trend  ==!
          CASE( np_ENT )                        !* energy conserving scheme  (T-pts)
-                             CALL vor_enT( kt, ntot, un , vn , ua, va )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_enT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
+                             CALL vor_enT( kt, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! total vorticity trend
+            IF( ln_stcor )   CALL vor_enT( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend
          CASE( np_EET )                        !* energy conserving scheme (een scheme using e3t)
-                             CALL vor_eeT( kt, ntot, un , vn , ua, va )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
+                             CALL vor_eeT( kt, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! total vorticity trend
+            IF( ln_stcor )   CALL vor_eeT( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend
          CASE( np_ENE )                        !* energy conserving scheme
-                             CALL vor_ene( kt, ntot, un , vn , ua, va )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
+                             CALL vor_ene( kt, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! total vorticity trend
+            IF( ln_stcor )   CALL vor_ene( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend
          CASE( np_ENS )                        !* enstrophy conserving scheme
-                             CALL vor_ens( kt, ntot, un , vn , ua, va )  ! total vorticity trend
-            IF( ln_stcor )   CALL vor_ens( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
+                             CALL vor_ens( kt, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )  ! total vorticity trend
+            IF( ln_stcor )   CALL vor_ens( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )  ! add the Stokes-Coriolis trend
          CASE( np_MIX )                        !* mixed ene-ens scheme
-                             CALL vor_ens( kt, nrvm, un , vn , ua, va )   ! relative vorticity or metric trend (ens)
-                             CALL vor_ene( kt, ncor, un , vn , ua, va )   ! planetary vorticity trend (ene)
-            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
+                             CALL vor_ens( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! relative vorticity or metric trend (ens)
+                             CALL vor_ene( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! planetary vorticity trend (ene)
+            IF( ln_stcor )   CALL vor_ene( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend
          CASE( np_EEN )                        !* energy and enstrophy conserving scheme
-                             CALL vor_een( kt, ntot, un , vn , ua, va )   ! total vorticity trend
-            IF( ln_stcor )   CALL vor_een( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend
+                             CALL vor_een( kt, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs )   ! total vorticity trend
+            IF( ln_stcor )   CALL vor_een( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs )   ! add the Stokes-Coriolis trend
          END SELECT
          !
@@ -187,5 +189,5 @@
 
 
-   SUBROUTINE vor_enT( kt, kvor, pu, pv, pu_rhs, pv_rhs )
+   SUBROUTINE vor_enT( kt, ktlev, kvor, pu, pv, pu_rhs, pv_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE vor_enT  ***
@@ -203,7 +205,8 @@
       !!       where rvor is the relative vorticity at f-point
       !!
-      !! ** Action : - Update (ua,va) with the now vorticity term trend
+      !! ** Action : - Update (u_rhs,v_rhs) with the now vorticity term trend
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
+      INTEGER                         , INTENT( in )  ::   ktlev            ! time level index for source terms
       INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
@@ -270,10 +273,10 @@
          SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
          CASE ( np_COR )                           !* Coriolis (planetary vorticity)
-            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t_n(:,:,jk)
+            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,ktlev)
          CASE ( np_RVO )                           !* relative vorticity
             DO jj = 2, jpj
                DO ji = 2, jpi   ! vector opt.
                   zwt(ji,jj) = r1_4 * (   zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)   &
-                     &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk)
+                     &                  + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,ktlev)
                END DO
             END DO
@@ -282,5 +285,5 @@
                DO ji = 2, jpi
                   zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   &
-                     &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t_n(ji,jj,jk)
+                     &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t(ji,jj,jk,ktlev)
                END DO
             END DO
@@ -289,5 +292,5 @@
                DO ji = 2, jpi   ! vector opt.
                   zwt(ji,jj) = (  ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj  ,jk) + zwz(ji,jj  ,jk)    &
-                     &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk)
+                     &                                 + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) )  ) * e1e2t(ji,jj)*e3t(ji,jj,jk,ktlev)
                END DO
             END DO
@@ -297,5 +300,5 @@
                   zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           &
                        &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  &
-                       &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t_n(ji,jj,jk)
+                       &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t(ji,jj,jk,ktlev)
                END DO
             END DO
@@ -307,9 +310,9 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1   ! vector opt.
-               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)                    &
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev)                    &
                   &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
                   &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   )
                   !
-               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)                    &
+               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev)                    &
                   &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
                   &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
@@ -322,5 +325,5 @@
 
 
-   SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva )
+   SUBROUTINE vor_ene( kt, ktlev, kvor, pu, pv, pu_rhs, pva_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE vor_ene  ***
@@ -334,16 +337,17 @@
       !!         The general trend of momentum is increased due to the vorticity 
       !!       term which is given by:
-      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v vn) ]
-      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u un) ]
+      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f  mi(e1v*e3v vv(:,:,:,ktlev)) ]
+      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u uu(:,:,:,ktlev)) ]
       !!       where rvor is the relative vorticity
       !!
-      !! ** Action : - Update (ua,va) with the now vorticity term trend
+      !! ** Action : - Update (u_rhs,v_rhs) with the now vorticity term trend
       !!
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
+      INTEGER                         , INTENT( in )  ::   ktlev            ! time level index for source terms
       INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_rhs    ! total v-trend
       !
       INTEGER  ::   ji, jj, jk           ! dummy loop indices
@@ -368,6 +372,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
-                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
+                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
+                     &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
                END DO
             END DO
@@ -375,6 +379,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
+                  zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
                END DO
             END DO
@@ -382,6 +386,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      &
-                     &                        - e1u(ji,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
+                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
+                     &                        - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
                END DO
             END DO
@@ -389,6 +393,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
+                  zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
                END DO
             END DO
@@ -407,9 +411,9 @@
          IF( ln_sco ) THEN
             zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
-            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
-            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
+            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,ktlev) * pu(:,:,jk)
+            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,ktlev) * pv(:,:,jk)
          ELSE
-            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
-            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
+            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
+            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
          ENDIF
          !                                   !==  compute and add the vorticity term trend  =!
@@ -420,6 +424,6 @@
                zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
                zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1)
-               pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-               pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
+               pva_rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
             END DO  
          END DO  
@@ -430,5 +434,5 @@
 
 
-   SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva )
+   SUBROUTINE vor_ens( kt, ktlev, kvor, pu, pv, pu_rhs, pva_rhs )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE vor_ens  ***
@@ -441,17 +445,18 @@
       !!      potential enstrophy of a horizontally non-divergent flow. the
       !!      trend of the vorticity term is given by:
-      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v vn) ]
-      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u un) ]
-      !!      Add this trend to the general momentum trend (ua,va):
-      !!          (ua,va) = (ua,va) + ( voru , vorv )
-      !!
-      !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
+      !!          voru = 1/e1u  mj-1[ (rvor+f)/e3f ]  mj-1[ mi(e1v*e3v vv(:,:,:,ktlev)) ]
+      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u uu(:,:,:,ktlev)) ]
+      !!      Add this trend to the general momentum trend (u_rhs,v_rhs):
+      !!          (u_rhs,v_rhs) = (u_rhs,v_rhs) + ( voru , vorv )
+      !!
+      !! ** Action : - Update (u_rhs,v_rhs) arrays with the now vorticity term trend
       !!
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
+      INTEGER                         , INTENT( in )  ::   ktlev            ! time level index for source terms
       INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_rhs    ! total v-trend
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -475,6 +480,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
-                     &          - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
+                  zwz(ji,jj) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
+                     &          - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
                END DO
             END DO
@@ -482,6 +487,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &       - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
+                  zwz(ji,jj) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &       - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
                END DO
             END DO
@@ -489,6 +494,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  &
-                     &                        - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
+                  zwz(ji,jj) = ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
+                     &                        - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
                END DO
             END DO
@@ -496,6 +501,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &                     - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
+                  zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &                     - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
                END DO
             END DO
@@ -514,9 +519,9 @@
          IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
             zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
-            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
-            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
+            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,ktlev) * pu(:,:,jk)
+            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,ktlev) * pv(:,:,jk)
          ELSE
-            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
-            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
+            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
+            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
          ENDIF
          !                                   !==  compute and add the vorticity term trend  =!
@@ -527,6 +532,6 @@
                zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  &
                   &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  )
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
-               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
+               pva_rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
             END DO  
          END DO  
@@ -537,5 +542,5 @@
 
 
-   SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva )
+   SUBROUTINE vor_een( kt, ktlev, kvor, pu, pv, pu_rhs, pva_rhs )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE vor_een  ***
@@ -548,14 +553,15 @@
       !!      both the horizontal kinetic energy and the potential enstrophy
       !!      when horizontal divergence is zero (see the NEMO documentation)
-      !!      Add this trend to the general momentum trend (ua,va).
-      !!
-      !! ** Action : - Update (ua,va) with the now vorticity term trend
+      !!      Add this trend to the general momentum trend (u_rhs,v_rhs).
+      !!
+      !! ** Action : - Update (u_rhs,v_rhs) with the now vorticity term trend
       !!
       !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
+      INTEGER                         , INTENT( in )  ::   ktlev            ! time level index for source terms
       INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_rhs    ! total v-trend
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -582,6 +588,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  ze3f = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
-                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
+                  ze3f = (  e3t(ji,jj+1,jk,ktlev)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,ktlev)*tmask(ji+1,jj+1,jk)   &
+                     &    + e3t(ji,jj  ,jk,ktlev)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,ktlev)*tmask(ji+1,jj  ,jk)  )
                   IF( ze3f /= 0._wp ) THEN   ;   z1_e3f(ji,jj) = 4._wp / ze3f
                   ELSE                       ;   z1_e3f(ji,jj) = 0._wp
@@ -592,6 +598,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  ze3f = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
-                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
+                  ze3f = (  e3t(ji,jj+1,jk,ktlev)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,ktlev)*tmask(ji+1,jj+1,jk)   &
+                     &    + e3t(ji,jj  ,jk,ktlev)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,ktlev)*tmask(ji+1,jj  ,jk)  )
                   zmsk = (                    tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
                      &                      + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk)  )
@@ -613,6 +619,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  &
-                     &            - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
+                  zwz(ji,jj,jk) = ( e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
+                     &            - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
                END DO
             END DO
@@ -620,6 +626,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj,jk) = (   ( pvn(ji+1,jj,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &              - ( pun(ji,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
+                  zwz(ji,jj,jk) = (   ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &              - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
                END DO
             END DO
@@ -627,6 +633,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      &
-                     &                              - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  )   &
+                  zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)      &
+                     &                              - e1u(ji  ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  )   &
                      &                           * r1_e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
                END DO
@@ -635,6 +641,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &                            - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
+                  zwz(ji,jj,jk) = (   ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &                            - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) * z1_e3f(ji,jj)
                END DO
             END DO
@@ -657,6 +663,6 @@
          !
          !                                   !==  horizontal fluxes  ==!
-         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
-         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
+         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,ktlev) * pu(:,:,jk)
+         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,ktlev) * pv(:,:,jk)
 
          !                                   !==  compute and add the vorticity term trend  =!
@@ -683,6 +689,6 @@
                zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
                   &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
-               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
+               pva_rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) + zva
             END DO  
          END DO  
@@ -694,5 +700,5 @@
 
 
-   SUBROUTINE vor_eeT( kt, kvor, pun, pvn, pua, pva )
+   SUBROUTINE vor_eeT( kt, ktlev, kvor, pu, pv, pu_rhs, pva_rhs )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE vor_eeT  ***
@@ -705,14 +711,15 @@
       !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
       !!      The change consists in 
-      !!      Add this trend to the general momentum trend (ua,va).
-      !!
-      !! ** Action : - Update (ua,va) with the now vorticity term trend
+      !!      Add this trend to the general momentum trend (u_rhs,v_rhs).
+      !!
+      !! ** Action : - Update (u_rhs,v_rhs) with the now vorticity term trend
       !!
       !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index
+      INTEGER                         , INTENT( in )  ::   ktlev            ! time level index for source terms
       INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_rhs    ! total v-trend
       !
       INTEGER  ::   ji, jj, jk     ! dummy loop indices
@@ -746,6 +753,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
-                     &             - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
+                  zwz(ji,jj,jk) = (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
+                     &             - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
                      &          * r1_e1e2f(ji,jj)
                END DO
@@ -754,6 +761,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj,jk) = ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &          - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
+                  zwz(ji,jj,jk) = ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &          - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
                END DO
             END DO
@@ -761,6 +768,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj  ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)    &
-                     &                              - e1u(ji  ,jj+1) * pun(ji  ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) &
+                  zwz(ji,jj,jk) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pv(ji+1,jj  ,jk) - e2v(ji,jj) * pv(ji,jj,jk)    &
+                     &                              - e1u(ji  ,jj+1) * pu(ji  ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk)  ) &
                      &                         * r1_e1e2f(ji,jj)    )
                END DO
@@ -769,6 +776,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zwz(ji,jj,jk) = ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
-                     &                        - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
+                  zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj  ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &                        - ( pu(ji  ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
                END DO
             END DO
@@ -791,6 +798,6 @@
 
       !                                   !==  horizontal fluxes  ==!
-         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
-         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
+         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,ktlev) * pu(:,:,jk)
+         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,ktlev) * pv(:,:,jk)
 
          !                                   !==  compute and add the vorticity term trend  =!
@@ -798,5 +805,5 @@
          ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0
          DO ji = 2, jpi          ! split in 2 parts due to vector opt.
-               z1_e3t = 1._wp / e3t_n(ji,jj,jk)
+               z1_e3t = 1._wp / e3t(ji,jj,jk,ktlev)
                ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
                ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
@@ -806,5 +813,5 @@
          DO jj = 3, jpj
             DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3
-               z1_e3t = 1._wp / e3t_n(ji,jj,jk)
+               z1_e3t = 1._wp / e3t(ji,jj,jk,ktlev)
                ztne(ji,jj) = ( zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) + zwz(ji  ,jj-1,jk) ) * z1_e3t
                ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj  ,jk) + zwz(ji  ,jj  ,jk) ) * z1_e3t
@@ -819,6 +826,6 @@
                zva = - r1_12 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   &
                   &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zua
-               pva(ji,jj,jk) = pva(ji,jj,jk) + zva
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
+               pva_rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) + zva
             END DO  
          END DO  
@@ -866,10 +873,10 @@
          WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
          WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
-         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
+         WRITE(numout,*) '      masked (=T) or uu(:,:,:,ktlev)masked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
       ENDIF
 
       IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available')
 
-!!gm  this should be removed when choosing a unique strategy for fmask at the coast
+!!gm  this should be removed when choosing a uu(:,:,:,ktlev)ique strategy for fmask at the coast
       ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
       ! at angles with three ocean points and one land point
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzad.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzad.F90	(revision 10788)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzad.F90	(revision 10789)
@@ -36,5 +36,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_zad ( kt )
+   SUBROUTINE dyn_zad ( kt, ktlev, pu_rhs, pv_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dynzad  ***
@@ -44,13 +44,15 @@
       !!
       !! ** Method  :   The now vertical advection of momentum is given by:
-      !!         w dz(u) = ua + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*wn) dk(un) ]
-      !!         w dz(v) = va + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*wn) dk(vn) ]
-      !!      Add this trend to the general trend (ua,va):
-      !!         (ua,va) = (ua,va) + w dz(u,v)
+      !!         w dz(u) = pu_rhs + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*ww) dk(uu(:,:,:,ktlev)) ]
+      !!         w dz(v) = pv_rhs + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*ww) dk(vv(:,:,:,ktlev)) ]
+      !!      Add this trend to the general trend (pu_rhs,pv_rhs):
+      !!         (pu_rhs,pv_rhs) = (pu_rhs,pv_rhs) + w dz(u,v)
       !!
-      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
+      !! ** Action  : - Update (pu_rhs,pv_rhs) with the vert. momentum adv. trends
       !!              - Send the trends to trddyn for diagnostics (l_trddyn=T)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
+      INTEGER, INTENT(in) ::   kt             ! ocean time-step index
+      INTEGER, INTENT(in) ::   ktlev          ! time level index for source terms
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -68,8 +70,8 @@
       ENDIF
 
-      IF( l_trddyn )   THEN         ! Save ua and va trends
+      IF( l_trddyn )   THEN         ! Save pu_rhs and pv_rhs trends
          ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
-         ztrdu(:,:,:) = ua(:,:,:) 
-         ztrdv(:,:,:) = va(:,:,:) 
+         ztrdu(:,:,:) = pu_rhs(:,:,:) 
+         ztrdv(:,:,:) = pv_rhs(:,:,:) 
       ENDIF
       
@@ -77,11 +79,11 @@
          DO jj = 2, jpj                   ! vertical fluxes 
             DO ji = fs_2, jpi             ! vector opt.
-               zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
+               zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
             END DO
          END DO
          DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
             DO ji = fs_2, fs_jpim1        ! vector opt.
-               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )
-               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )
+               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( uu(ji,jj,jk-1,ktlev) - uu(ji,jj,jk,ktlev) )
+               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vv(ji,jj,jk-1,ktlev) - vv(ji,jj,jk,ktlev) )
             END DO  
          END DO   
@@ -101,6 +103,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1       ! vector opt.
-               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
-               va(ji,jj,jk) = va(ji,jj,jk) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev)
+               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev)
             END DO  
          END DO  
@@ -108,6 +110,6 @@
 
       IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
-         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
          DEALLOCATE( ztrdu, ztrdv ) 
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 10788)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/oce.F90	(revision 10789)
@@ -20,5 +20,5 @@
    !! --------------------------                            
    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(:,:,:),   TARGET ::   ww             !: vertical velocity            [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           wi             !: vertical vel. (adaptive-implicit) [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::           hdivn          !: horizontal divergence        [s-1]
@@ -71,4 +71,5 @@
    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]
+   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::           wn             !: k-vertical   velocity        [m/s]
 
    !!----------------------------------------------------------------------
@@ -88,5 +89,5 @@
       ierr(:) = 0 
       ALLOCATE( uu   (jpi,jpj,jpk, jpt) , vv   (jpi,jpj,jpk,jpt)  ,                             &
-         &      wn   (jpi,jpj,jpk)      , hdivn(jpi,jpj,jpk)      ,                             &
+         &      ww   (jpi,jpj,jpk)      , hdivn(jpi,jpj,jpk)      ,                             &
          &      tsb  (jpi,jpj,jpk,jpts) , tsn  (jpi,jpj,jpk,jpts) , tsa(jpi,jpj,jpk,jpts) ,     &
          &      rab_b(jpi,jpj,jpk,jpts) , rab_n(jpi,jpj,jpk,jpts) ,                             &
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 10788)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10789)
@@ -188,6 +188,6 @@
                &         CALL Agrif_Sponge_dyn        ! momentum sponge
 #endif
-                         CALL dyn_adv       ( kstp )  ! advection (vector or flux form)
-                         CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis
+                         CALL dyn_adv       ( kstp, Nm1, Nnn, uu(:,:,:,Nrhs), vv(:,:,:,Nrhs) )  ! advection (vector or flux form)
+                         CALL dyn_vor       ( kstp,      Nnn, uu(:,:,:,Nrhs), vv(:,:,:,Nrhs) )  ! vorticity term including Coriolis
                          CALL dyn_ldf       ( kstp )  ! lateral mixing
       IF( ln_zdfosm  )   CALL dyn_osm       ( kstp )  ! OSMOSIS non-local velocity fluxes
@@ -355,4 +355,5 @@
       ub => uu(:,:,:,Nm1); un => uu(:,:,:,Nnn); ua => uu(:,:,:,Np1)
       vb => vv(:,:,:,Nm1); vn => vv(:,:,:,Nnn); va => vv(:,:,:,Np1)
+      wn => ww(:,:,:)
 
       e3t_b => e3t(:,:,:,Nm1); e3t_n => e3t(:,:,:,Nnn); e3t_a => e3t(:,:,:,Np1)
