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 10874)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90	(revision 10877)
@@ -96,15 +96,17 @@
 CONTAINS
 
-   SUBROUTINE dyn_vor( kt )
+   SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
       !!----------------------------------------------------------------------
       !!
       !! ** Purpose :   compute the lateral ocean tracer physics.
       !!
-      !! ** Action : - Update (ua,va) with the now vorticity term trend
+      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) 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  ) ::   Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! ocean velocity field and RHS of momentum equation
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv
@@ -117,34 +119,34 @@
          ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
          !
-         ztrdu(:,:,:) = ua(:,:,:)            !* planetary vorticity trend (including Stokes-Coriolis force)
-         ztrdv(:,:,:) = va(:,:,:)
+         ztrdu(:,:,:) = puu(:,:,:,Krhs)            !* planetary vorticity trend (including Stokes-Coriolis force)
+         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
          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, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! enstrophy conserving scheme
+            IF( ln_stcor )            CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+         CASE( np_ENE, np_MIX )   ;   CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme
+            IF( ln_stcor )            CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+         CASE( np_ENT )           ;   CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (T-pts)
+            IF( ln_stcor )            CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+         CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy conserving scheme (een with e3t)
+            IF( ln_stcor )            CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
+         CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! energy & enstrophy scheme
+            IF( ln_stcor )            CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend
          END SELECT
-         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          !
          IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case)
-            ztrdu(:,:,:) = ua(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:)
+            ztrdu(:,:,:) = puu(:,:,:,Krhs)
+            ztrdv(:,:,:) = pvv(:,:,:,Krhs)
             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, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (T-pts)
+            CASE( np_EET )           ;   CALL vor_eeT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme (een with e3t)
+            CASE( np_ENE )           ;   CALL vor_ene( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy conserving scheme
+            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! enstrophy conserving scheme
+            CASE( np_EEN )           ;   CALL vor_een( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! energy & enstrophy scheme
             END SELECT
-            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+            ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
+            ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
             CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
          ENDIF
@@ -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, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
+            IF( ln_stcor )   CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! 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, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
+            IF( ln_stcor )   CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! 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, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
+            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! 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, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend
+            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! 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, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! relative vorticity or metric trend (ens)
+                             CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! planetary vorticity trend (ene)
+            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! 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, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend
+            IF( ln_stcor )   CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! 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, Kmm, 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 (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index
       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,Kmm)
          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,Kmm)
                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,Kmm)
                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,Kmm)
                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,Kmm)
                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,Kmm)                    &
                   &                                * (  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,Kmm)                    &
                   &                                * (  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, Kmm, kvor, pu_mm, pv_mm, pu_rhs, pv_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 pvv(:,:,:,Kmm)) ]
+      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f  mj(e2u*e3u puu(:,:,:,Kmm)) ]
       !!       where rvor is the relative vorticity
       !!
-      !! ** Action : - Update (ua,va) with the now vorticity term trend
+      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) 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   ) ::   Kmm              ! ocean time level index
       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_mm, pv_mm    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_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_mm(ji+1,jj  ,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)    &
+                     &          - e1u(ji  ,jj+1) * pu_mm(ji  ,jj+1,jk) + e1u(ji,jj) * pu_mm(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_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &       - ( pu_mm(ji  ,jj+1,jk) + pu_mm(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_mm(ji+1,jj,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)      &
+                     &                        - e1u(ji,jj+1) * pu_mm(ji,jj+1,jk) + e1u(ji,jj) * pu_mm(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_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &                     - ( pu_mm(ji  ,jj+1,jk) + pu_mm(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
                END DO
             END DO
@@ -406,10 +410,10 @@
 
          IF( ln_sco ) THEN
-            zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk)
-            zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk)
-            zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk)
+            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
+            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu_mm(:,:,jk)
+            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv_mm(:,:,jk)
          ELSE
-            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
-            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
+            zwx(:,:) = e2u(:,:) * pu_mm(:,:,jk)
+            zwy(:,:) = e1v(:,:) * pv_mm(:,:,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 )
+               pv_rhs(ji,jj,jk) = pv_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, Kmm, kvor, pu_mm, pv_mm, pu_rhs, pv_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 pvv(:,:,:,Kmm)) ]
+      !!          vorv = 1/e2v  mi-1[ (rvor+f)/e3f ]  mi-1[ mj(e2u*e3u puu(:,:,:,Kmm)) ]
+      !!      Add this trend to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)):
+      !!          (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) = (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) + ( voru , vorv )
+      !!
+      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) 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   ) ::   Kmm              ! ocean time level index
       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_mm, pv_mm    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_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_mm(ji+1,jj  ,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)    &
+                     &          - e1u(ji  ,jj+1) * pu_mm(ji  ,jj+1,jk) + e1u(ji,jj) * pu_mm(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_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &       - ( pu_mm(ji  ,jj+1,jk) + pu_mm(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_mm(ji+1,jj  ,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)  &
+                     &                        - e1u(ji  ,jj+1) * pu_mm(ji  ,jj+1,jk) + e1u(ji,jj) * pu_mm(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_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &                     - ( pu_mm(ji  ,jj+1,jk) + pu_mm(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
                END DO
             END DO
@@ -513,10 +518,10 @@
          !
          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)
+            zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
+            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu_mm(:,:,jk)
+            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv_mm(:,:,jk)
          ELSE
-            zwx(:,:) = e2u(:,:) * pun(:,:,jk)
-            zwy(:,:) = e1v(:,:) * pvn(:,:,jk)
+            zwx(:,:) = e2u(:,:) * pu_mm(:,:,jk)
+            zwy(:,:) = e1v(:,:) * pv_mm(:,:,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) )
+               pv_rhs(ji,jj,jk) = pv_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, Kmm, kvor, pu_mm, pv_mm, pu_rhs, pv_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 (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
+      !!
+      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) 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   ) ::   Kmm              ! ocean time level index
       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_mm, pv_mm    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_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,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
+                     &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*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,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
+                     &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*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_mm(ji+1,jj,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)  &
+                     &            - e1u(ji  ,jj+1) * pu_mm(ji,jj+1,jk) + e1u(ji,jj) * pu_mm(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_mm(ji+1,jj,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &              - ( pu_mm(ji,jj+1,jk) + pu_mm(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_mm(ji+1,jj,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)      &
+                     &                              - e1u(ji  ,jj+1) * pu_mm(ji,jj+1,jk) + e1u(ji,jj) * pu_mm(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_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &                            - ( pu_mm(ji  ,jj+1,jk) + pu_mm(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,Kmm) * pu_mm(:,:,jk)
+         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv_mm(:,:,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
+               pv_rhs(ji,jj,jk) = pv_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, Kmm, kvor, pu_mm, pv_mm, pu_rhs, pv_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 (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
+      !!
+      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) 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   ) ::   Kmm              ! ocean time level index
       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_mm, pv_mm    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_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_mm(ji+1,jj  ,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)    &
+                     &             - e1u(ji  ,jj+1) * pu_mm(ji  ,jj+1,jk) + e1u(ji,jj) * pu_mm(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_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &          - ( pu_mm(ji  ,jj+1,jk) + pu_mm(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_mm(ji+1,jj  ,jk) - e2v(ji,jj) * pv_mm(ji,jj,jk)    &
+                     &                              - e1u(ji  ,jj+1) * pu_mm(ji  ,jj+1,jk) + e1u(ji,jj) * pu_mm(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_mm(ji+1,jj  ,jk) + pv_mm(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   &
+                     &                        - ( pu_mm(ji  ,jj+1,jk) + pu_mm(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,Kmm) * pu_mm(:,:,jk)
+         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv_mm(:,:,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,Kmm)
                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,Kmm)
                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
+               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
             END DO  
          END DO  
