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 10884)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv.F90	(revision 10893)
@@ -75,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, Kmm, puu, pvv, Krhs )    ! vector form : horizontal gradient of kinetic energy
-         CALL dyn_zad     ( kt,            Kmm, puu, pvv, Krhs )    ! vector form : vertical advection
+         CALL dyn_keg     ( kt, nn_dynkeg,      Kmm, puu, pvv, Krhs )    ! vector form : horizontal gradient of kinetic energy
+         CALL dyn_zad     ( kt,                 Kmm, puu, pvv, Krhs )    ! vector form : vertical advection
       CASE( np_FLX_c2  ) 
-         CALL dyn_adv_cen2( kt,      Kmm, puu, pvv, Krhs )          ! 2nd order centered scheme
+         CALL dyn_adv_cen2( kt,                 Kmm, puu, pvv, Krhs )    ! 2nd order centered scheme
       CASE( np_FLX_ubs )   
-         CALL dyn_adv_ubs ( kt, Kbb, Kmm, puu, pvv, Krhs )          ! 3rd order UBS      scheme (UP3)
+         CALL dyn_adv_ubs ( kt,            Kbb, Kmm, puu, pvv, Krhs )    ! 3rd order UBS      scheme (UP3)
       END SELECT
       !
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 10884)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynkeg.F90	(revision 10893)
@@ -57,11 +57,11 @@
       !!              * kscheme = nkeg_HW : Hollingsworth correction following
       !!      Arakawa (2001). The now horizontal kinetic energy is given by:
-      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((puu(j+1,Kmm)+puu(j-1,Kmm))/2)^2  )
-      !!                    + mj-1(  2 * vn^2 + ((pvv(i+1,Kmm)+pvv(i-1,Kmm))/2)^2  ) ]
+      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((u(j+1)+u(j-1))/2)^2  )
+      !!                    + mj-1(  2 * vn^2 + ((v(i+1)+v(i-1))/2)^2  ) ]
       !!      
       !!      Take its horizontal gradient and add it to the general momentum
-      !!      trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
-      !!         puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u di[ zhke ]
-      !!         pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v dj[ zhke ]
+      !!      trend.
+      !!         u(rhs) = u(rhs) - 1/e1u di[ zhke ]
+      !!         v(rhs) = v(rhs) - 1/e2v dj[ zhke ]
       !!
       !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend
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 10884)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90	(revision 10893)
@@ -205,5 +205,5 @@
       !!       where rvor is the relative vorticity at f-point
       !!
-      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
+      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
@@ -325,5 +325,5 @@
 
 
-   SUBROUTINE vor_ene( kt, Kmm, kvor, pu_mm, pv_mm, pu_rhs, pv_rhs )
+   SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE vor_ene  ***
@@ -341,5 +341,5 @@
       !!       where rvor is the relative vorticity
       !!
-      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
+      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
       !!
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
@@ -348,5 +348,5 @@
       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_mm, pv_mm    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
       !
@@ -372,6 +372,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  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)
+                  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
@@ -379,6 +379,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  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)
+                  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
@@ -386,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) * 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)
+                  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
@@ -393,6 +393,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  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)
+                  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
@@ -411,9 +411,9 @@
          IF( ln_sco ) THEN
             zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
-            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu_mm(:,:,jk)
-            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv_mm(:,:,jk)
+            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
+            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
          ELSE
-            zwx(:,:) = e2u(:,:) * pu_mm(:,:,jk)
-            zwy(:,:) = e1v(:,:) * pv_mm(:,:,jk)
+            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
+            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
          ENDIF
          !                                   !==  compute and add the vorticity term trend  =!
@@ -434,5 +434,5 @@
 
 
-   SUBROUTINE vor_ens( kt, Kmm, kvor, pu_mm, pv_mm, pu_rhs, pv_rhs )
+   SUBROUTINE vor_ens( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE vor_ens  ***
@@ -447,8 +447,8 @@
       !!          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
+      !!      Add this trend to the general momentum trend:
+      !!          (u(rhs),v(Krhs)) = (u(rhs),v(Krhs)) + ( voru , vorv )
+      !!
+      !! ** Action : - Update (pu_rhs,pv_rhs)) arrays with the now vorticity term trend
       !!
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
@@ -457,5 +457,5 @@
       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_mm, pv_mm    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
       !
@@ -480,6 +480,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  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)
+                  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
@@ -487,6 +487,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  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)
+                  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
@@ -494,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  ) * 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)
+                  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
@@ -501,6 +501,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  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)
+                  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
@@ -519,9 +519,9 @@
          IF( ln_sco ) THEN                   !==  horizontal fluxes  ==!
             zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
-            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu_mm(:,:,jk)
-            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv_mm(:,:,jk)
+            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
+            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
          ELSE
-            zwx(:,:) = e2u(:,:) * pu_mm(:,:,jk)
-            zwy(:,:) = e1v(:,:) * pv_mm(:,:,jk)
+            zwx(:,:) = e2u(:,:) * pu(:,:,jk)
+            zwy(:,:) = e1v(:,:) * pv(:,:,jk)
          ENDIF
          !                                   !==  compute and add the vorticity term trend  =!
@@ -542,5 +542,5 @@
 
 
-   SUBROUTINE vor_een( kt, Kmm, kvor, pu_mm, pv_mm, pu_rhs, pv_rhs )
+   SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE vor_een  ***
@@ -553,7 +553,7 @@
       !!      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 (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
-      !!
-      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
+      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
+      !!
+      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
       !!
       !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
@@ -562,5 +562,5 @@
       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_mm, pv_mm    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
       !
@@ -619,6 +619,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  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)
+                  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
@@ -626,6 +626,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  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)
+                  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
@@ -633,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  ) * 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)  )   &
+                  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
@@ -641,6 +641,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  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)
+                  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
@@ -663,6 +663,6 @@
          !
          !                                   !==  horizontal fluxes  ==!
-         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu_mm(:,:,jk)
-         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv_mm(:,:,jk)
+         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
+         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
 
          !                                   !==  compute and add the vorticity term trend  =!
@@ -700,5 +700,5 @@
 
 
-   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu_mm, pv_mm, pu_rhs, pv_rhs )
+   SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE vor_eeT  ***
@@ -711,7 +711,7 @@
       !!      a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
       !!      The change consists in 
-      !!      Add this trend to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
-      !!
-      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
+      !!      Add this trend to the general momentum trend (pu_rhs,pv_rhs).
+      !!
+      !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
       !!
       !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
@@ -720,5 +720,5 @@
       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_mm, pv_mm    ! now velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend
       !
@@ -753,6 +753,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  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)  ) &
+                  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
@@ -761,6 +761,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  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)
+                  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
@@ -768,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  ) * 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)  ) &
+                  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
@@ -776,6 +776,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  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)
+                  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
@@ -798,6 +798,6 @@
 
       !                                   !==  horizontal fluxes  ==!
-         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu_mm(:,:,jk)
-         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv_mm(:,:,jk)
+         zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
+         zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
 
          !                                   !==  compute and add the vorticity term trend  =!
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 10884)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzad.F90	(revision 10893)
@@ -44,8 +44,8 @@
       !!
       !! ** Method  :   The now vertical advection of momentum is given by:
-      !!         w dz(u) = puu(:,:,:,Krhs) + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*ww) dk(puu(:,:,:,Kmm)) ]
-      !!         w dz(v) = pvv(:,:,:,Krhs) + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*ww) dk(pvv(:,:,:,Kmm)) ]
+      !!         w dz(u) = u(rhs) + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*ww) dk(u) ]
+      !!         w dz(v) = v(rhs) + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*ww) dk(v) ]
       !!      Add this trend to the general trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)):
-      !!         (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) = (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) + w dz(u,v)
+      !!         (u(rhs),v(rhs)) = (u(rhs),v(rhs)) + w dz(u,v)
       !!
       !! ** Action  : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the vert. momentum adv. trends
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90	(revision 10884)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90	(revision 10893)
@@ -54,9 +54,9 @@
       !!
       !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
-      !!         puu(:,:,:,Kaa) =         puu(:,:,:,Kbb) + 2*dt *       puu(:,:,:,Krhs)             vector form or linear free surf.
-      !!         puu(:,:,:,Kaa) = ( e3u_b*puu(:,:,:,Kbb) + 2*dt * e3u_n*puu(:,:,:,Krhs) ) / e3u(:,:,:,Kaa)   otherwise
+      !!         u(after) =         u(before) + 2*dt *       u(rhs)                vector form or linear free surf.
+      !!         u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u(after)   otherwise
       !!               - update the after velocity with the implicit vertical mixing.
       !!      This requires to solver the following system: 
-      !!         puu(:,:,:,Kaa) = puu(:,:,:,Kaa) + 1/e3u(:,:,:,Kaa) dk+1[ mi(avm) / e3uw(:,:,:,Kaa) dk[ua] ]
+      !!         u(after) = u(after) + 1/e3u(after) dk+1[ mi(avm) / e3uw(after) dk[ua] ]
       !!      with the following surface/top/bottom boundary condition:
       !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_cen.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_cen.F90	(revision 10884)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_cen.F90	(revision 10893)
@@ -44,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pu_mm, pv_mm, pww,     &
+   SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pU, pV, pW,     &
       &                    Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v ) 
       !!----------------------------------------------------------------------
@@ -70,6 +70,6 @@
       INTEGER                                  , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme)
       INTEGER                                  , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme)
-      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm, pv_mm, pww   ! 3 ocean velocity components
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -106,6 +106,6 @@
                DO jj = 1, jpjm1
                   DO ji = 1, fs_jpim1   ! vector opt.
-                     zwx(ji,jj,jk) = 0.5_wp * pu_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
-                     zwy(ji,jj,jk) = 0.5_wp * pv_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
+                     zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
+                     zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
                   END DO
                END DO
@@ -134,6 +134,6 @@
                      zC4t_v =  zC2t_v + r1_6 * ( ztv(ji,jj-1,jk) - ztv(ji,jj+1,jk) )
                      !                                                  ! C4 fluxes
-                     zwx(ji,jj,jk) =  0.5_wp * pu_mm(ji,jj,jk) * zC4t_u
-                     zwy(ji,jj,jk) =  0.5_wp * pv_mm(ji,jj,jk) * zC4t_v
+                     zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u
+                     zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v
                   END DO
                END DO
@@ -150,5 +150,5 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     zwz(ji,jj,jk) = 0.5 * pww(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) * wmask(ji,jj,jk)
+                     zwz(ji,jj,jk) = 0.5 * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) * wmask(ji,jj,jk)
                   END DO
                END DO
@@ -160,5 +160,5 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     zwz(ji,jj,jk) = pww(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
+                     zwz(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
                   END DO
                END DO
@@ -171,9 +171,9 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zwz(ji,jj, mikt(ji,jj) ) = pww(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) 
+                     zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) 
                   END DO
                END DO   
             ELSE                                   ! no ice-shelf cavities (only ocean surface)
-               zwz(:,:,1) = pww(:,:,1) * pt(:,:,1,jn,Kmm)
+               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kmm)
             ENDIF
          ENDIF
@@ -191,7 +191,7 @@
          !                             ! trend diagnostics
          IF( l_trd ) THEN
-            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pu_mm, pt(:,:,:,jn,Kmm) )
-            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv_mm, pt(:,:,:,jn,Kmm) )
-            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pww, pt(:,:,:,jn,Kmm) )
+            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
+            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
+            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
          END IF
          !                                 ! "Poleward" heat and salt transports 
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_fct.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_fct.F90	(revision 10884)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_fct.F90	(revision 10893)
@@ -52,5 +52,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pu_mm, pv_mm, pww,       &
+   SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pU, pV, pW,       &
       &                    Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v )
       !!----------------------------------------------------------------------
@@ -77,6 +77,6 @@
       INTEGER                                  , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4)
       REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm, pv_mm, pww   ! 3 ocean velocity components
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices  
@@ -125,8 +125,8 @@
                DO ji = 1, fs_jpim1   ! vector opt.
                   ! upstream scheme
-                  zfp_ui = pu_mm(ji,jj,jk) + ABS( pu_mm(ji,jj,jk) )
-                  zfm_ui = pu_mm(ji,jj,jk) - ABS( pu_mm(ji,jj,jk) )
-                  zfp_vj = pv_mm(ji,jj,jk) + ABS( pv_mm(ji,jj,jk) )
-                  zfm_vj = pv_mm(ji,jj,jk) - ABS( pv_mm(ji,jj,jk) )
+                  zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )
+                  zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
+                  zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
+                  zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) )
                   zwx(ji,jj,jk) = 0.5 * ( zfp_ui * pt(ji,jj,jk,jn,Kbb) + zfm_ui * pt(ji+1,jj  ,jk,jn,Kbb) )
                   zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji  ,jj+1,jk,jn,Kbb) )
@@ -138,6 +138,6 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zfp_wk = pww(ji,jj,jk) + ABS( pww(ji,jj,jk) )
-                  zfm_wk = pww(ji,jj,jk) - ABS( pww(ji,jj,jk) )
+                  zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
+                  zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
                   zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb) ) * wmask(ji,jj,jk)
                END DO
@@ -148,9 +148,9 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zwz(ji,jj, mikt(ji,jj) ) = pww(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
+                     zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
                   END DO
                END DO   
             ELSE                             ! no cavities: only at the ocean surface
-               zwz(:,:,1) = pww(:,:,1) * pt(:,:,1,jn,Kbb)
+               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
             ENDIF
          ENDIF
@@ -184,6 +184,6 @@
                DO jj = 1, jpjm1
                   DO ji = 1, fs_jpim1   ! vector opt.
-                     zwx(ji,jj,jk) = 0.5_wp * pu_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx(ji,jj,jk)
-                     zwy(ji,jj,jk) = 0.5_wp * pv_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy(ji,jj,jk)
+                     zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx(ji,jj,jk)
+                     zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy(ji,jj,jk)
                   END DO
                END DO
@@ -215,6 +215,6 @@
                      zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
                      !                                                  ! C4 minus upstream advective fluxes 
-                     zwx(ji,jj,jk) =  0.5_wp * pu_mm(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
-                     zwy(ji,jj,jk) =  0.5_wp * pv_mm(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
+                     zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
+                     zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
                   END DO
                END DO
@@ -243,6 +243,6 @@
                      zC4t_v =  zC2t_v + r1_6 * ( ztv(ji  ,jj-1,jk) - ztv(ji  ,jj+1,jk) )
                      !                                                  ! C4 minus upstream advective fluxes 
-                     zwx(ji,jj,jk) =  0.5_wp * pu_mm(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
-                     zwy(ji,jj,jk) =  0.5_wp * pv_mm(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
+                     zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
+                     zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
                   END DO
                END DO
@@ -257,5 +257,5 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     zwz(ji,jj,jk) =  (  pww(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
+                     zwz(ji,jj,jk) =  (  pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
                         &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk)
                   END DO
@@ -268,5 +268,5 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     zwz(ji,jj,jk) = ( pww(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
+                     zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
                   END DO
                END DO
@@ -303,7 +303,7 @@
             !
             IF( l_trd ) THEN              ! trend diagnostics
-               CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pu_mm, pt(:,:,:,jn,Kmm) )
-               CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pv_mm, pt(:,:,:,jn,Kmm) )
-               CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pww, pt(:,:,:,jn,Kmm) )
+               CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pU, pt(:,:,:,jn,Kmm) )
+               CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pV, pt(:,:,:,jn,Kmm) )
+               CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pW, pt(:,:,:,jn,Kmm) )
             ENDIF
             !                             ! heat/salt transport
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_mus.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_mus.F90	(revision 10884)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_mus.F90	(revision 10893)
@@ -54,5 +54,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pu_mm, pv_mm, pww,             &
+   SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pU, pV, pW,             &
       &                    Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
       !!----------------------------------------------------------------------
@@ -80,6 +80,6 @@
       LOGICAL                                  , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl
       REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm, pv_mm, pww   ! 3 ocean velocity components
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -172,17 +172,17 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! MUSCL fluxes
-                  z0u = SIGN( 0.5, pu_mm(ji,jj,jk) )
+                  z0u = SIGN( 0.5, pU(ji,jj,jk) )
                   zalpha = 0.5 - z0u
-                  zu  = z0u - 0.5 * pu_mm(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+                  zu  = z0u - 0.5 * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
                   zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
                   zzwy = pt(ji  ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk)
-                  zwx(ji,jj,jk) = pu_mm(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
+                  zwx(ji,jj,jk) = pU(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
                   !
-                  z0v = SIGN( 0.5, pv_mm(ji,jj,jk) )
+                  z0v = SIGN( 0.5, pV(ji,jj,jk) )
                   zalpha = 0.5 - z0v
-                  zv  = z0v - 0.5 * pv_mm(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
+                  zv  = z0v - 0.5 * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
                   zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
                   zzwy = pt(ji,jj  ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk)
-                  zwy(ji,jj,jk) = pv_mm(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
+                  zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
                END DO
             END DO
@@ -201,6 +201,6 @@
          !                                ! trend diagnostics
          IF( l_trd )  THEN
-            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pu_mm, pt(:,:,:,jn,Kbb) )
-            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv_mm, pt(:,:,:,jn,Kbb) )
+            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kbb) )
+            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) )
          END IF
          !                                 ! "Poleward" heat and salt transports 
@@ -239,10 +239,10 @@
             DO jj = 2, jpjm1      
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z0w = SIGN( 0.5, pww(ji,jj,jk+1) )
+                  z0w = SIGN( 0.5, pW(ji,jj,jk+1) )
                   zalpha = 0.5 + z0w
-                  zw  = z0w - 0.5 * pww(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm)
+                  zw  = z0w - 0.5 * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm)
                   zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
                   zzwy = pt(ji,jj,jk  ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  )
-                  zwx(ji,jj,jk+1) = pww(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
+                  zwx(ji,jj,jk+1) = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
                END DO 
             END DO
@@ -252,9 +252,9 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zwx(ji,jj, mikt(ji,jj) ) = pww(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)
+                     zwx(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)
                   END DO
                END DO   
             ELSE                                      ! no cavities: only at the ocean surface
-               zwx(:,:,1) = pww(:,:,1) * pt(:,:,1,jn,Kbb)
+               zwx(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
             ENDIF
          ENDIF
@@ -268,5 +268,5 @@
          END DO
          !                                ! send trends for diagnostic
-         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pww, pt(:,:,:,jn,Kbb) )
+         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pW, pt(:,:,:,jn,Kbb) )
          !
       END DO                     ! end of tracer loop
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90	(revision 10884)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90	(revision 10893)
@@ -47,5 +47,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pu_mm, pv_mm, pww, Kbb, Kmm, pt, kjpt, Krhs )
+   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_adv_qck  ***
@@ -89,6 +89,6 @@
       INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
       REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm, pv_mm, pww   ! 3 ocean velocity components
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
       !!----------------------------------------------------------------------
       !
@@ -107,14 +107,14 @@
       !
       !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
-      CALL tra_adv_qck_i( kt, cdtype, p2dt, pu_mm, Kbb, Kmm, pt, kjpt, Krhs ) 
-      CALL tra_adv_qck_j( kt, cdtype, p2dt, pv_mm, Kbb, Kmm, pt, kjpt, Krhs ) 
+      CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs ) 
+      CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs ) 
 
       !        ! vertical fluxes are computed with the 2nd order centered scheme
-      CALL tra_adv_cen2_k( kt, cdtype, pww, Kmm, pt, kjpt, Krhs )
+      CALL tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
       !
    END SUBROUTINE tra_adv_qck
 
 
-   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pu_mm, Kbb, Kmm, pt, kjpt, Krhs )
+   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs )
       !!----------------------------------------------------------------------
       !!
@@ -125,5 +125,5 @@
       INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
       REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm        ! i-velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU        ! i-velocity components
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
       !!
@@ -156,5 +156,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.         
-                  zdir = 0.5 + SIGN( 0.5, pu_mm(ji,jj,jk) )   ! if pu_mm > 0 : zdir = 1 otherwise zdir = 0 
+                  zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
                   zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T 
                END DO
@@ -165,7 +165,7 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.   
-                  zdir = 0.5 + SIGN( 0.5, pu_mm(ji,jj,jk) )   ! if pu_mm > 0 : zdir = 1 otherwise zdir = 0 
+                  zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
                   zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
-                  zwx(ji,jj,jk)  = ABS( pu_mm(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
+                  zwx(ji,jj,jk)  = ABS( pU(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
                   zfc(ji,jj,jk)  = zdir * pt(ji  ,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji+1,jj,jk,jn,Kbb)  ! FC in the x-direction for T
                   zfd(ji,jj,jk)  = zdir * pt(ji+1,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji  ,jj,jk,jn,Kbb)  ! FD in the x-direction for T
@@ -195,10 +195,10 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.               
-                  zdir = 0.5 + SIGN( 0.5, pu_mm(ji,jj,jk) )   ! if pu_mm > 0 : zdir = 1 otherwise zdir = 0 
+                  zdir = 0.5 + SIGN( 0.5, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
                   !--- If the second ustream point is a land point
                   !--- the flux is computed by the 1st order UPWIND scheme
                   zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
                   zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
-                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pu_mm(ji,jj,jk)
+                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk)
                END DO
             END DO
@@ -220,5 +220,5 @@
          END DO
          !                                 ! trend diagnostics
-         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pu_mm, pt(:,:,:,jn,Kmm) )
+         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
          !
       END DO
@@ -227,5 +227,5 @@
 
 
-   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pv_mm, Kbb, Kmm, pt, kjpt, Krhs )
+   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )
       !!----------------------------------------------------------------------
       !!
@@ -236,5 +236,5 @@
       INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
       REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pv_mm        ! j-velocity components
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pV        ! j-velocity components
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
       !!
@@ -272,5 +272,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.         
-                  zdir = 0.5 + SIGN( 0.5, pv_mm(ji,jj,jk) )   ! if pu_mm > 0 : zdir = 1 otherwise zdir = 0 
+                  zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
                   zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T 
                END DO
@@ -281,7 +281,7 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.   
-                  zdir = 0.5 + SIGN( 0.5, pv_mm(ji,jj,jk) )   ! if pu_mm > 0 : zdir = 1 otherwise zdir = 0 
+                  zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
                   zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
-                  zwy(ji,jj,jk)  = ABS( pv_mm(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
+                  zwy(ji,jj,jk)  = ABS( pV(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
                   zfc(ji,jj,jk)  = zdir * pt(ji,jj  ,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji,jj+1,jk,jn,Kbb)  ! FC in the x-direction for T
                   zfd(ji,jj,jk)  = zdir * pt(ji,jj+1,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji,jj  ,jk,jn,Kbb)  ! FD in the x-direction for T
@@ -311,10 +311,10 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.               
-                  zdir = 0.5 + SIGN( 0.5, pv_mm(ji,jj,jk) )   ! if pu_mm > 0 : zdir = 1 otherwise zdir = 0 
+                  zdir = 0.5 + SIGN( 0.5, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0 
                   !--- If the second ustream point is a land point
                   !--- the flux is computed by the 1st order UPWIND scheme
                   zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
                   zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
-                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pv_mm(ji,jj,jk)
+                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk)
                END DO
             END DO
@@ -336,5 +336,5 @@
          END DO
          !                                 ! trend diagnostics
-         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv_mm, pt(:,:,:,jn,Kmm) )
+         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
          IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
@@ -345,5 +345,5 @@
 
 
-   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pww, Kmm, pt, kjpt, Krhs )
+   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
       !!----------------------------------------------------------------------
       !!
@@ -353,5 +353,5 @@
       CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
       INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pww      ! vertical velocity 
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity 
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
       !
@@ -370,5 +370,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zwz(ji,jj,jk) = 0.5 * pww(ji,jj,jk) * ( pt(ji,jj,jk-1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk)
+                  zwz(ji,jj,jk) = 0.5 * pW(ji,jj,jk) * ( pt(ji,jj,jk-1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk)
                END DO
             END DO
@@ -378,9 +378,9 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zwz(ji,jj, mikt(ji,jj) ) = pww(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface 
+                     zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface 
                   END DO
                END DO   
             ELSE                                   ! no ocean cavities (only ocean surface)
-               zwz(:,:,1) = pww(:,:,1) * pt(:,:,1,jn,Kmm)
+               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kmm)
             ENDIF
          ENDIF
@@ -395,5 +395,5 @@
          END DO
          !                                 ! Send trends for diagnostic
-         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pww, pt(:,:,:,jn,Kmm) )
+         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
          !
       END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90	(revision 10884)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90	(revision 10893)
@@ -46,5 +46,5 @@
 CONTAINS
 
-   SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pu_mm, pv_mm, pww,          &
+   SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pU, pV, pW,          &
       &                    Kbb, Kmm, pt, kjpt, Krhs, kn_ubs_v )
       !!----------------------------------------------------------------------
@@ -58,7 +58,7 @@
       !!      It is only used in the horizontal direction.
       !!      For example the i-component of the advective fluxes are given by :
-      !!                !  e2u e3u uu ( mi(Tn) - zltu(i  ) ,Kmm)   if uu(i,Kmm) >= 0
+      !!                !  e2u e3u un ( mi(Tn) - zltu(i  ) )   if un(i) >= 0
       !!          ztu = !  or 
-      !!                !  e2u e3u uu ( mi(Tn) - zltu(i+1) ,Kmm)   if uu(i,Kmm) < 0
+      !!                !  e2u e3u un ( mi(Tn) - zltu(i+1) )   if un(i) < 0
       !!      where zltu is the second derivative of the before temperature field:
       !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
@@ -91,5 +91,5 @@
       INTEGER                                  , INTENT(in   ) ::   kn_ubs_v        ! number of tracers
       REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pu_mm, pv_mm, pww   ! 3 ocean transport components
+      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
       !
@@ -146,11 +146,11 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zfp_ui = pu_mm(ji,jj,jk) + ABS( pu_mm(ji,jj,jk) )      ! upstream transport (x2)
-                  zfm_ui = pu_mm(ji,jj,jk) - ABS( pu_mm(ji,jj,jk) )
-                  zfp_vj = pv_mm(ji,jj,jk) + ABS( pv_mm(ji,jj,jk) )
-                  zfm_vj = pv_mm(ji,jj,jk) - ABS( pv_mm(ji,jj,jk) )
+                  zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )      ! upstream transport (x2)
+                  zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
+                  zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
+                  zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) )
                   !                                                  ! 2nd order centered advective fluxes (x2)
-                  zcenut = pu_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
-                  zcenvt = pv_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
+                  zcenut = pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
+                  zcenvt = pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
                   !                                                  ! UBS advective fluxes
                   ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) )
@@ -177,6 +177,6 @@
          !                
          IF( l_trd ) THEN                  ! trend diagnostics
-             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztu, pu_mm, pt(:,:,:,jn,Kmm) )
-             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pv_mm, pt(:,:,:,jn,Kmm) )
+             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztu, pU, pt(:,:,:,jn,Kmm) )
+             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pV, pt(:,:,:,jn,Kmm) )
          END IF
          !     
@@ -199,6 +199,6 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zfp_wk = pww(ji,jj,jk) + ABS( pww(ji,jj,jk) )
-                     zfm_wk = pww(ji,jj,jk) - ABS( pww(ji,jj,jk) )
+                     zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
+                     zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
                      ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb)  ) * wmask(ji,jj,jk)
                   END DO
@@ -209,9 +209,9 @@
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        ztw(ji,jj, mikt(ji,jj) ) = pww(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
+                        ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
                      END DO
                   END DO   
                ELSE                                ! no cavities: only at the ocean surface
-                  ztw(:,:,1) = pww(:,:,1) * pt(:,:,1,jn,Kbb)
+                  ztw(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
                ENDIF
             ENDIF
@@ -232,5 +232,5 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     ztw(ji,jj,jk) = (   0.5_wp * pww(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
+                     ztw(ji,jj,jk) = (   0.5_wp * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
                         &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk)
                   END DO
@@ -247,9 +247,9 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1
-                     ztw(ji,jj,jk) = pww(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-            IF( ln_linssh )   ztw(:,:, 1 ) = pww(:,:,1) * pt(:,:,1,jn,Kmm)     !!gm ISF & 4th COMPACT doesn't work
+                     ztw(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            IF( ln_linssh )   ztw(:,:, 1 ) = pW(:,:,1) * pt(:,:,1,jn,Kmm)     !!gm ISF & 4th COMPACT doesn't work
             !
          END SELECT
@@ -268,5 +268,5 @@
                   DO ji = fs_2, fs_jpim1   ! vector opt.
                      zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltv(ji,jj,jk)                          &
-                        &           + pt(ji,jj,jk,jn,Kmm) * (  pww(ji,jj,jk) - pww(ji,jj,jk+1)  )   &
+                        &           + pt(ji,jj,jk,jn,Kmm) * (  pW(ji,jj,jk) - pW(ji,jj,jk+1)  )   &
                         &                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                   END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90	(revision 10884)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90	(revision 10893)
@@ -77,5 +77,5 @@
       !
       !                                      !* compute lateral mixing trend and add it to the general trend
-      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, pts, Kaa, jpts ) 
+      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 
 
 !!gm WHY here !   and I don't like that !
@@ -109,5 +109,5 @@
 
  
-   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, pt, Kaa, kjpt ) 
+   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt ) 
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_zdf_imp  ***
@@ -130,5 +130,5 @@
       !!---------------------------------------------------------------------
       INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index
-      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! ocean time level indices
+      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs, Kaa  ! ocean time level indices
       INTEGER                                  , INTENT(in   ) ::   kit000   ! first time step index
       CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
@@ -240,5 +240,5 @@
          DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
             DO ji = fs_2, fs_jpim1
-               pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Kaa)
+               pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
             END DO
          END DO
@@ -246,5 +246,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Kaa)   ! zrhs=right hand side
+                  zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side
                   pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
                END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trczdf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trczdf.F90	(revision 10884)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trczdf.F90	(revision 10893)
@@ -56,5 +56,5 @@
       IF( l_trdtrc )   ztrtrd(:,:,:,:)  = ptr(:,:,:,:,Krhs)
       !
-      CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dttrc, Kbb, Kmm, ptr, Kaa, jptra )    !   implicit scheme          
+      CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dttrc, Kbb, Kmm, Krhs, ptr, Kaa, jptra )    !   implicit scheme          
       !
       IF( l_trdtrc )   THEN                      ! save the vertical diffusive trends for further diagnostics
