Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce.F90	(revision 11027)
@@ -66,4 +66,5 @@
    INTEGER, PUBLIC :: umsk_id, vmsk_id
    INTEGER, PUBLIC :: kindic_agr
+   INTEGER, PUBLIC :: Kbb_a, Kmm_a, Krhs_a
    
    !!----------------------------------------------------------------------
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_interp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_interp.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_interp.F90	(revision 11027)
@@ -107,13 +107,13 @@
          !
          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            uu_b(ibdy1:ibdy2,:,Krhs) = 0._wp
+            uu_b(ibdy1:ibdy2,:,Krhs_a) = 0._wp
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  uu_b(ibdy1:ibdy2,jj,Krhs) = uu_b(ibdy1:ibdy2,jj,Krhs) & 
-                      & + e3u(ibdy1:ibdy2,jj,jk,Krhs) * uu(ibdy1:ibdy2,jj,jk,Krhs) * umask(ibdy1:ibdy2,jj,jk)
+                  uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) & 
+                      & + e3u(ibdy1:ibdy2,jj,jk,Krhs_a) * uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
             DO jj = 1, jpj
-               uu_b(ibdy1:ibdy2,jj,Krhs) = uu_b(ibdy1:ibdy2,jj,Krhs) * r1_hu_a(ibdy1:ibdy2,jj)
+               uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu_a(ibdy1:ibdy2,jj)
             END DO
          ENDIF
@@ -122,5 +122,5 @@
             DO jk=1,jpkm1              ! Smooth
                DO jj=j1,j2
-                  uu(ibdy2,jj,jk,Krhs) = 0.25_wp*(uu(ibdy2-1,jj,jk,Krhs)+2._wp*uu(ibdy2,jj,jk,Krhs)+uu(ibdy2+1,jj,jk,Krhs))
+                  uu(ibdy2,jj,jk,Krhs_a) = 0.25_wp*(uu(ibdy2-1,jj,jk,Krhs_a)+2._wp*uu(ibdy2,jj,jk,Krhs_a)+uu(ibdy2+1,jj,jk,Krhs_a))
                END DO
             END DO
@@ -131,5 +131,5 @@
             DO jj = 1, jpj
                zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) & 
-                  & + e3u(ibdy1:ibdy2,jj,jk,Krhs)  * uu(ibdy1:ibdy2,jj,jk,Krhs)*umask(ibdy1:ibdy2,jj,jk)
+                  & + e3u(ibdy1:ibdy2,jj,jk,Krhs_a)  * uu(ibdy1:ibdy2,jj,jk,Krhs_a)*umask(ibdy1:ibdy2,jj,jk)
             END DO
          END DO
@@ -140,6 +140,6 @@
          DO jk = 1, jpkm1
             DO jj = 1, jpj
-               uu(ibdy1:ibdy2,jj,jk,Krhs) = ( uu(ibdy1:ibdy2,jj,jk,Krhs) &
-                 & + uu_b(ibdy1:ibdy2,jj,Krhs)-zub(ibdy1:ibdy2,jj)) * umask(ibdy1:ibdy2,jj,jk)
+               uu(ibdy1:ibdy2,jj,jk,Krhs_a) = ( uu(ibdy1:ibdy2,jj,jk,Krhs_a) &
+                 & + uu_b(ibdy1:ibdy2,jj,Krhs_a)-zub(ibdy1:ibdy2,jj)) * umask(ibdy1:ibdy2,jj,jk)
             END DO
          END DO
@@ -150,5 +150,5 @@
                DO jj = 1, jpj
                   zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) & 
-                     & + e3v(ibdy1:ibdy2,jj,jk,Krhs) * vv(ibdy1:ibdy2,jj,jk,Krhs) * vmask(ibdy1:ibdy2,jj,jk)
+                     & + e3v(ibdy1:ibdy2,jj,jk,Krhs_a) * vv(ibdy1:ibdy2,jj,jk,Krhs_a) * vmask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -158,6 +158,6 @@
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  vv(ibdy1:ibdy2,jj,jk,Krhs) = ( vv(ibdy1:ibdy2,jj,jk,Krhs) & 
-                    & + vv_b(ibdy1:ibdy2,jj,Krhs)-zvb(ibdy1:ibdy2,jj))*vmask(ibdy1:ibdy2,jj,jk)
+                  vv(ibdy1:ibdy2,jj,jk,Krhs_a) = ( vv(ibdy1:ibdy2,jj,jk,Krhs_a) & 
+                    & + vv_b(ibdy1:ibdy2,jj,Krhs_a)-zvb(ibdy1:ibdy2,jj))*vmask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -166,6 +166,6 @@
          DO jk = 1, jpkm1              ! Mask domain edges
             DO jj = 1, jpj
-               uu(1,jj,jk,Krhs) = 0._wp
-               vv(1,jj,jk,Krhs) = 0._wp
+               uu(1,jj,jk,Krhs_a) = 0._wp
+               vv(1,jj,jk,Krhs_a) = 0._wp
             END DO
          END DO 
@@ -178,13 +178,13 @@
          !
          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            uu_b(ibdy1:ibdy2,:,Krhs) = 0._wp
+            uu_b(ibdy1:ibdy2,:,Krhs_a) = 0._wp
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  uu_b(ibdy1:ibdy2,jj,Krhs) = uu_b(ibdy1:ibdy2,jj,Krhs) & 
-                      & + e3u(ibdy1:ibdy2,jj,jk,Krhs) * uu(ibdy1:ibdy2,jj,jk,Krhs) * umask(ibdy1:ibdy2,jj,jk)
+                  uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) & 
+                      & + e3u(ibdy1:ibdy2,jj,jk,Krhs_a) * uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
             DO jj = 1, jpj
-               uu_b(ibdy1:ibdy2,jj,Krhs) = uu_b(ibdy1:ibdy2,jj,Krhs) * r1_hu_a(ibdy1:ibdy2,jj)
+               uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu_a(ibdy1:ibdy2,jj)
             END DO
          ENDIF
@@ -193,5 +193,5 @@
             DO jk=1,jpkm1              ! Smooth
                DO jj=j1,j2
-                  uu(ibdy1,jj,jk,Krhs) = 0.25_wp*(uu(ibdy1-1,jj,jk,Krhs)+2._wp*uu(ibdy1,jj,jk,Krhs)+uu(ibdy1+1,jj,jk,Krhs))
+                  uu(ibdy1,jj,jk,Krhs_a) = 0.25_wp*(uu(ibdy1-1,jj,jk,Krhs_a)+2._wp*uu(ibdy1,jj,jk,Krhs_a)+uu(ibdy1+1,jj,jk,Krhs_a))
                END DO
             END DO
@@ -202,5 +202,5 @@
             DO jj = 1, jpj
                zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) & 
-                  & + e3u(ibdy1:ibdy2,jj,jk,Krhs)  * uu(ibdy1:ibdy2,jj,jk,Krhs) * umask(ibdy1:ibdy2,jj,jk)
+                  & + e3u(ibdy1:ibdy2,jj,jk,Krhs_a)  * uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk)
             END DO
          END DO
@@ -211,6 +211,6 @@
          DO jk = 1, jpkm1
             DO jj = 1, jpj
-               uu(ibdy1:ibdy2,jj,jk,Krhs) = ( uu(ibdy1:ibdy2,jj,jk,Krhs) & 
-                 & + uu_b(ibdy1:ibdy2,jj,Krhs)-zub(ibdy1:ibdy2,jj))*umask(ibdy1:ibdy2,jj,jk)
+               uu(ibdy1:ibdy2,jj,jk,Krhs_a) = ( uu(ibdy1:ibdy2,jj,jk,Krhs_a) & 
+                 & + uu_b(ibdy1:ibdy2,jj,Krhs_a)-zub(ibdy1:ibdy2,jj))*umask(ibdy1:ibdy2,jj,jk)
             END DO
          END DO
@@ -223,5 +223,5 @@
                DO jj = 1, jpj
                   zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) &
-                     & + e3v(ibdy1:ibdy2,jj,jk,Krhs) * vv(ibdy1:ibdy2,jj,jk,Krhs) * vmask(ibdy1:ibdy2,jj,jk)
+                     & + e3v(ibdy1:ibdy2,jj,jk,Krhs_a) * vv(ibdy1:ibdy2,jj,jk,Krhs_a) * vmask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -231,6 +231,6 @@
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  vv(ibdy1:ibdy2,jj,jk,Krhs) = ( vv(ibdy1:ibdy2,jj,jk,Krhs) & 
-                      & + vv_b(ibdy1:ibdy2,jj,Krhs)-zvb(ibdy1:ibdy2,jj)) * vmask(ibdy1:ibdy2,jj,jk)
+                  vv(ibdy1:ibdy2,jj,jk,Krhs_a) = ( vv(ibdy1:ibdy2,jj,jk,Krhs_a) & 
+                      & + vv_b(ibdy1:ibdy2,jj,Krhs_a)-zvb(ibdy1:ibdy2,jj)) * vmask(ibdy1:ibdy2,jj,jk)
                END DO
             END DO
@@ -239,6 +239,6 @@
          DO jk = 1, jpkm1              ! Mask domain edges
             DO jj = 1, jpj
-               uu(nlci-1,jj,jk,Krhs) = 0._wp
-               vv(nlci  ,jj,jk,Krhs) = 0._wp
+               uu(nlci-1,jj,jk,Krhs_a) = 0._wp
+               vv(nlci  ,jj,jk,Krhs_a) = 0._wp
             END DO
          END DO 
@@ -251,13 +251,13 @@
          !
          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            vv_b(:,jbdy1:jbdy2,Krhs) = 0._wp
+            vv_b(:,jbdy1:jbdy2,Krhs_a) = 0._wp
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  vv_b(ji,jbdy1:jbdy2,Krhs) = vv_b(ji,jbdy1:jbdy2,Krhs) & 
-                      & + e3v(ji,jbdy1:jbdy2,jk,Krhs) * vv(ji,jbdy1:jbdy2,jk,Krhs) * vmask(ji,jbdy1:jbdy2,jk)
+                  vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) & 
+                      & + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) * vv(ji,jbdy1:jbdy2,jk,Krhs_a) * vmask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
             DO ji=1,jpi
-               vv_b(ji,jbdy1:jbdy2,Krhs) = vv_b(ji,jbdy1:jbdy2,Krhs) * r1_hv_a(ji,jbdy1:jbdy2)
+               vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv_a(ji,jbdy1:jbdy2)
             END DO
          ENDIF
@@ -266,5 +266,5 @@
             DO jk = 1, jpkm1           ! Smooth
                DO ji = i1, i2
-                  vv(ji,jbdy2,jk,Krhs) = 0.25_wp*(vv(ji,jbdy2-1,jk,Krhs)+2._wp*vv(ji,jbdy2,jk,Krhs)+vv(ji,jbdy2+1,jk,Krhs))
+                  vv(ji,jbdy2,jk,Krhs_a) = 0.25_wp*(vv(ji,jbdy2-1,jk,Krhs_a)+2._wp*vv(ji,jbdy2,jk,Krhs_a)+vv(ji,jbdy2+1,jk,Krhs_a))
                END DO
             END DO
@@ -275,5 +275,5 @@
             DO ji=1,jpi
                zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) & 
-                  & + e3v(ji,jbdy1:jbdy2,jk,Krhs) * vv(ji,jbdy1:jbdy2,jk,Krhs) * vmask(ji,jbdy1:jbdy2,jk)
+                  & + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) * vv(ji,jbdy1:jbdy2,jk,Krhs_a) * vmask(ji,jbdy1:jbdy2,jk)
             END DO
          END DO
@@ -284,6 +284,6 @@
          DO jk = 1, jpkm1
             DO ji = 1, jpi
-               vv(ji,jbdy1:jbdy2,jk,Krhs) = ( vv(ji,jbdy1:jbdy2,jk,Krhs) & 
-                 & + vv_b(ji,jbdy1:jbdy2,Krhs) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
+               vv(ji,jbdy1:jbdy2,jk,Krhs_a) = ( vv(ji,jbdy1:jbdy2,jk,Krhs_a) & 
+                 & + vv_b(ji,jbdy1:jbdy2,Krhs_a) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
             END DO
          END DO
@@ -294,5 +294,5 @@
                DO ji = 1, jpi
                   zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) & 
-                     & + e3u(ji,jbdy1:jbdy2,jk,Krhs) * uu(ji,jbdy1:jbdy2,jk,Krhs) * umask(ji,jbdy1:jbdy2,jk)
+                     & + e3u(ji,jbdy1:jbdy2,jk,Krhs_a) * uu(ji,jbdy1:jbdy2,jk,Krhs_a) * umask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -303,6 +303,6 @@
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  uu(ji,jbdy1:jbdy2,jk,Krhs) = ( uu(ji,jbdy1:jbdy2,jk,Krhs) & 
-                    & + uu_b(ji,jbdy1:jbdy2,Krhs) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
+                  uu(ji,jbdy1:jbdy2,jk,Krhs_a) = ( uu(ji,jbdy1:jbdy2,jk,Krhs_a) & 
+                    & + uu_b(ji,jbdy1:jbdy2,Krhs_a) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -311,6 +311,6 @@
          DO jk = 1, jpkm1              ! Mask domain edges
             DO ji = 1, jpi
-               uu(ji,1,jk,Krhs) = 0._wp
-               vv(ji,1,jk,Krhs) = 0._wp
+               uu(ji,1,jk,Krhs_a) = 0._wp
+               vv(ji,1,jk,Krhs_a) = 0._wp
             END DO
          END DO 
@@ -323,13 +323,13 @@
          !
          IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            vv_b(:,jbdy1:jbdy2,Krhs) = 0._wp
+            vv_b(:,jbdy1:jbdy2,Krhs_a) = 0._wp
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  vv_b(ji,jbdy1:jbdy2,Krhs) = vv_b(ji,jbdy1:jbdy2,Krhs) & 
-                      & + e3v(ji,jbdy1:jbdy2,jk,Krhs) * vv(ji,jbdy1:jbdy2,jk,Krhs) * vmask(ji,jbdy1:jbdy2,jk)
+                  vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) & 
+                      & + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) * vv(ji,jbdy1:jbdy2,jk,Krhs_a) * vmask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
             DO ji=1,jpi
-               vv_b(ji,jbdy1:jbdy2,Krhs) = vv_b(ji,jbdy1:jbdy2,Krhs) * r1_hv_a(ji,jbdy1:jbdy2)
+               vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv_a(ji,jbdy1:jbdy2)
             END DO
          ENDIF
@@ -338,5 +338,5 @@
             DO jk = 1, jpkm1           ! Smooth
                DO ji = i1, i2
-                  vv(ji,jbdy1,jk,Krhs) = 0.25_wp*(vv(ji,jbdy1-1,jk,Krhs)+2._wp*vv(ji,jbdy1,jk,Krhs)+vv(ji,jbdy1+1,jk,Krhs))
+                  vv(ji,jbdy1,jk,Krhs_a) = 0.25_wp*(vv(ji,jbdy1-1,jk,Krhs_a)+2._wp*vv(ji,jbdy1,jk,Krhs_a)+vv(ji,jbdy1+1,jk,Krhs_a))
                END DO
             END DO
@@ -347,5 +347,5 @@
             DO ji=1,jpi
                zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) & 
-                  & + e3v(ji,jbdy1:jbdy2,jk,Krhs) * vv(ji,jbdy1:jbdy2,jk,Krhs) * vmask(ji,jbdy1:jbdy2,jk)
+                  & + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) * vv(ji,jbdy1:jbdy2,jk,Krhs_a) * vmask(ji,jbdy1:jbdy2,jk)
             END DO
          END DO
@@ -356,6 +356,6 @@
          DO jk = 1, jpkm1
             DO ji = 1, jpi
-               vv(ji,jbdy1:jbdy2,jk,Krhs) = ( vv(ji,jbdy1:jbdy2,jk,Krhs) & 
-                 & + vv_b(ji,jbdy1:jbdy2,Krhs) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
+               vv(ji,jbdy1:jbdy2,jk,Krhs_a) = ( vv(ji,jbdy1:jbdy2,jk,Krhs_a) & 
+                 & + vv_b(ji,jbdy1:jbdy2,Krhs_a) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
             END DO
          END DO
@@ -368,5 +368,5 @@
                DO ji = 1, jpi
                   zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) & 
-                     & + e3u(ji,jbdy1:jbdy2,jk,Krhs) * uu(ji,jbdy1:jbdy2,jk,Krhs) * umask(ji,jbdy1:jbdy2,jk)
+                     & + e3u(ji,jbdy1:jbdy2,jk,Krhs_a) * uu(ji,jbdy1:jbdy2,jk,Krhs_a) * umask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -377,6 +377,6 @@
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  uu(ji,jbdy1:jbdy2,jk,Krhs) = ( uu(ji,jbdy1:jbdy2,jk,Krhs) & 
-                    & + uu_b(ji,jbdy1:jbdy2,Krhs) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
+                  uu(ji,jbdy1:jbdy2,jk,Krhs_a) = ( uu(ji,jbdy1:jbdy2,jk,Krhs_a) & 
+                    & + uu_b(ji,jbdy1:jbdy2,Krhs_a) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
                END DO
             END DO
@@ -385,6 +385,6 @@
          DO jk = 1, jpkm1              ! Mask domain edges
             DO ji = 1, jpi
-               uu(ji,nlcj  ,jk,Krhs) = 0._wp
-               vv(ji,nlcj-1,jk,Krhs) = 0._wp
+               uu(ji,nlcj  ,jk,Krhs_a) = 0._wp
+               vv(ji,nlcj-1,jk,Krhs_a) = 0._wp
             END DO
          END DO 
@@ -520,5 +520,5 @@
          DO jj = 1, jpj
             DO ji = 2, indx
-               ssh(ji,jj,Krhs) = hbdy_w(ji-1,jj)
+               ssh(ji,jj,Krhs_a) = hbdy_w(ji-1,jj)
             ENDDO
          ENDDO
@@ -530,5 +530,5 @@
          DO jj = 1, jpj
             DO ji = indx, nlci-1
-               ssh(ji,jj,Krhs) = hbdy_e(ji-indx+1,jj)
+               ssh(ji,jj,Krhs_a) = hbdy_e(ji-indx+1,jj)
             ENDDO
          ENDDO
@@ -540,5 +540,5 @@
          DO jj = 2, indy
             DO ji = 1, jpi
-               ssh(ji,jj,Krhs) = hbdy_s(ji,jj-1)
+               ssh(ji,jj,Krhs_a) = hbdy_s(ji,jj-1)
             ENDDO
          ENDDO
@@ -550,5 +550,5 @@
          DO jj = indy, nlcj-1
             DO ji = 1, jpi
-               ssh(ji,jj,Krhs) = hbdy_n(ji,jj-indy+1)
+               ssh(ji,jj,Krhs_a) = hbdy_n(ji,jj-indy+1)
             ENDDO
          ENDDO
@@ -659,5 +659,5 @@
                DO jj=j1,j2
                  DO ji=i1,i2
-                       ptab(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm)
+                       ptab(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)
                  END DO
               END DO
@@ -669,5 +669,5 @@
            DO jj=j1,j2
               DO ji=i1,i2
-                 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) 
+                 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
               END DO
            END DO
@@ -699,5 +699,5 @@
                   IF (tmask(iref,jref,jk) == 0) EXIT 
                   N_out = N_out + 1
-                  h_out(jk) = e3t(iref,jref,jk,Kmm)
+                  h_out(jk) = e3t(iref,jref,jk,Kmm_a)
                ENDDO
                IF (N_in > 0) THEN
@@ -713,5 +713,5 @@
          !
          DO jn=1, jpts
-            ts(i1:i2,j1:j2,1:jpk,jn,Krhs)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
+            ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
          END DO
 
@@ -737,14 +737,14 @@
                ibdy = nlci-nbghostcells
                DO jn = 1, jpts
-                  ts(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
+                  ts(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
                   DO jk = 1, jpkm1
                      DO jj = jmin,jmax
                         IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN
-                           ts(ibdy,jj,jk,jn,Krhs) = ts(ibdy+1,jj,jk,jn,Krhs) * tmask(ibdy,jj,jk)
+                           ts(ibdy,jj,jk,jn,Krhs_a) = ts(ibdy+1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)
                         ELSE
-                           ts(ibdy,jj,jk,jn,Krhs)=(z4*ts(ibdy+1,jj,jk,jn,Krhs)+z3*ts(ibdy-1,jj,jk,jn,Krhs))*tmask(ibdy,jj,jk)
-                           IF( uu(ibdy-1,jj,jk,Kmm) > 0._wp ) THEN
-                              ts(ibdy,jj,jk,jn,Krhs)=( z6*ts(ibdy-1,jj,jk,jn,Krhs)+z5*ts(ibdy+1,jj,jk,jn,Krhs) & 
-                                                 + z7*ts(ibdy-2,jj,jk,jn,Krhs) ) * tmask(ibdy,jj,jk)
+                           ts(ibdy,jj,jk,jn,Krhs_a)=(z4*ts(ibdy+1,jj,jk,jn,Krhs_a)+z3*ts(ibdy-1,jj,jk,jn,Krhs_a))*tmask(ibdy,jj,jk)
+                           IF( uu(ibdy-1,jj,jk,Kmm_a) > 0._wp ) THEN
+                              ts(ibdy,jj,jk,jn,Krhs_a)=( z6*ts(ibdy-1,jj,jk,jn,Krhs_a)+z5*ts(ibdy+1,jj,jk,jn,Krhs_a) & 
+                                                 + z7*ts(ibdy-2,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)
                            ENDIF
                         ENDIF
@@ -752,5 +752,5 @@
                   END DO
                   ! Restore ghost points:
-                  ts(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy+1,jmin:jmax,1:jpkm1)
+                  ts(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy+1,jmin:jmax,1:jpkm1)
                END DO
             ENDIF
@@ -766,14 +766,14 @@
                jbdy = nlcj-nbghostcells         
                DO jn = 1, jpts
-                  ts(imin:imax,jbdy+1,1:jpkm1,jn,Krhs) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
+                  ts(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
                   DO jk = 1, jpkm1
                      DO ji = imin,imax
                         IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN
-                           ts(ji,jbdy,jk,jn,Krhs) = ts(ji,jbdy+1,jk,jn,Krhs) * tmask(ji,jbdy,jk)
+                           ts(ji,jbdy,jk,jn,Krhs_a) = ts(ji,jbdy+1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
                         ELSE
-                           ts(ji,jbdy,jk,jn,Krhs)=(z4*ts(ji,jbdy+1,jk,jn,Krhs)+z3*ts(ji,jbdy-1,jk,jn,Krhs))*tmask(ji,jbdy,jk)        
-                           IF (vv(ji,jbdy-1,jk,Kmm) > 0._wp ) THEN
-                              ts(ji,jbdy,jk,jn,Krhs)=( z6*ts(ji,jbdy-1,jk,jn,Krhs)+z5*ts(ji,jbdy+1,jk,jn,Krhs)  &
-                                                 + z7*ts(ji,jbdy-2,jk,jn,Krhs) ) * tmask(ji,jbdy,jk)
+                           ts(ji,jbdy,jk,jn,Krhs_a)=(z4*ts(ji,jbdy+1,jk,jn,Krhs_a)+z3*ts(ji,jbdy-1,jk,jn,Krhs_a))*tmask(ji,jbdy,jk)        
+                           IF (vv(ji,jbdy-1,jk,Kmm_a) > 0._wp ) THEN
+                              ts(ji,jbdy,jk,jn,Krhs_a)=( z6*ts(ji,jbdy-1,jk,jn,Krhs_a)+z5*ts(ji,jbdy+1,jk,jn,Krhs_a)  &
+                                                 + z7*ts(ji,jbdy-2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
                            ENDIF
                         ENDIF
@@ -781,5 +781,5 @@
                   END DO
                   ! Restore ghost points:
-                  ts(imin:imax,jbdy+1,1:jpkm1,jn,Krhs) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) * tmask(imin:imax,jbdy+1,1:jpkm1)
+                  ts(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) * tmask(imin:imax,jbdy+1,1:jpkm1)
                END DO
             ENDIF
@@ -795,14 +795,14 @@
                ibdy = 1+nbghostcells       
                DO jn = 1, jpts
-                  ts(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
+                  ts(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
                   DO jk = 1, jpkm1
                      DO jj = jmin,jmax
                         IF( umask(ibdy,jj,jk) == 0._wp ) THEN
-                           ts(ibdy,jj,jk,jn,Krhs) = ts(ibdy-1,jj,jk,jn,Krhs) * tmask(ibdy,jj,jk)
+                           ts(ibdy,jj,jk,jn,Krhs_a) = ts(ibdy-1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)
                         ELSE
-                           ts(ibdy,jj,jk,jn,Krhs)=(z4*ts(ibdy-1,jj,jk,jn,Krhs)+z3*ts(ibdy+1,jj,jk,jn,Krhs))*tmask(ibdy,jj,jk)        
-                           IF( uu(ibdy,jj,jk,Kmm) < 0._wp ) THEN
-                              ts(ibdy,jj,jk,jn,Krhs)=( z6*ts(ibdy+1,jj,jk,jn,Krhs)+z5*ts(ibdy-1,jj,jk,jn,Krhs) &
-                                                 + z7*ts(ibdy+2,jj,jk,jn,Krhs) ) * tmask(ibdy,jj,jk)
+                           ts(ibdy,jj,jk,jn,Krhs_a)=(z4*ts(ibdy-1,jj,jk,jn,Krhs_a)+z3*ts(ibdy+1,jj,jk,jn,Krhs_a))*tmask(ibdy,jj,jk)        
+                           IF( uu(ibdy,jj,jk,Kmm_a) < 0._wp ) THEN
+                              ts(ibdy,jj,jk,jn,Krhs_a)=( z6*ts(ibdy+1,jj,jk,jn,Krhs_a)+z5*ts(ibdy-1,jj,jk,jn,Krhs_a) &
+                                                 + z7*ts(ibdy+2,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)
                            ENDIF
                         ENDIF
@@ -810,5 +810,5 @@
                   END DO
                   ! Restore ghost points:
-                  ts(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy-1,jmin:jmax,1:jpkm1)
+                  ts(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy-1,jmin:jmax,1:jpkm1)
                END DO
             ENDIF
@@ -824,14 +824,14 @@
                jbdy=1+nbghostcells        
                DO jn = 1, jpts
-                  ts(imin:imax,jbdy-1,1:jpkm1,jn,Krhs) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
+                  ts(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
                   DO jk = 1, jpkm1      
                      DO ji = imin,imax
                         IF( vmask(ji,jbdy,jk) == 0._wp ) THEN
-                           ts(ji,jbdy,jk,jn,Krhs)=ts(ji,jbdy-1,jk,jn,Krhs) * tmask(ji,jbdy,jk)
+                           ts(ji,jbdy,jk,jn,Krhs_a)=ts(ji,jbdy-1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
                         ELSE
-                           ts(ji,jbdy,jk,jn,Krhs)=(z4*ts(ji,jbdy-1,jk,jn,Krhs)+z3*ts(ji,jbdy+1,jk,jn,Krhs))*tmask(ji,jbdy,jk)
-                           IF( vv(ji,jbdy,jk,Kmm) < 0._wp ) THEN
-                              ts(ji,jbdy,jk,jn,Krhs)=( z6*ts(ji,jbdy+1,jk,jn,Krhs)+z5*ts(ji,jbdy-1,jk,jn,Krhs) & 
-                                                 + z7*ts(ji,jbdy+2,jk,jn,Krhs) ) * tmask(ji,jbdy,jk)
+                           ts(ji,jbdy,jk,jn,Krhs_a)=(z4*ts(ji,jbdy-1,jk,jn,Krhs_a)+z3*ts(ji,jbdy+1,jk,jn,Krhs_a))*tmask(ji,jbdy,jk)
+                           IF( vv(ji,jbdy,jk,Kmm_a) < 0._wp ) THEN
+                              ts(ji,jbdy,jk,jn,Krhs_a)=( z6*ts(ji,jbdy+1,jk,jn,Krhs_a)+z5*ts(ji,jbdy-1,jk,jn,Krhs_a) & 
+                                                 + z7*ts(ji,jbdy+2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
                            ENDIF
                         ENDIF
@@ -839,5 +839,5 @@
                   END DO
                   ! Restore ghost points:
-                  ts(imin:imax,jbdy-1,1:jpkm1,jn,Krhs) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) * tmask(imin:imax,jbdy-1,1:jpkm1)
+                  ts(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) * tmask(imin:imax,jbdy-1,1:jpkm1)
                END DO
             ENDIF
@@ -861,5 +861,5 @@
       !
       IF( before) THEN
-         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm)
+         ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a)
       ELSE
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -900,7 +900,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm)*umask(ji,jj,jk)) 
+                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)*umask(ji,jj,jk)) 
 # if defined key_vertical
-                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm))
+                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a))
 # endif
                END DO
@@ -928,5 +928,5 @@
          
               IF (N_in == 0) THEN
-                 uu(ji,jj,:,Krhs) = 0._wp
+                 uu(ji,jj,:,Krhs_a) = 0._wp
                  CYCLE
               ENDIF
@@ -936,9 +936,9 @@
                  if (umask(iref,jj,jk) == 0) EXIT
                  N_out = N_out + 1
-                 h_out(N_out) = e3u(iref,jj,jk,Krhs)
+                 h_out(N_out) = e3u(iref,jj,jk,Krhs_a)
               ENDDO
          
               IF (N_out == 0) THEN
-                 uu(ji,jj,:,Krhs) = 0._wp
+                 uu(ji,jj,:,Krhs_a) = 0._wp
                  CYCLE
               ENDIF
@@ -952,5 +952,5 @@
                  endif
               ENDIF
-              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs),h_out(1:N_out),N_in,N_out)
+              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out)
             ENDDO
          ENDDO
@@ -959,5 +959,5 @@
          DO jk = 1, jpkm1
             DO jj=j1,j2
-               uu(i1:i2,jj,jk,Krhs) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs) )
+               uu(i1:i2,jj,jk,Krhs_a) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs_a) )
             END DO
          END DO
@@ -992,7 +992,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm)*vmask(ji,jj,jk))
+                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)*vmask(ji,jj,jk))
 # if defined key_vertical
-                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
+                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
 # endif
                END DO
@@ -1019,5 +1019,5 @@
                END DO
                IF (N_in == 0) THEN
-                  vv(ji,jj,:,Krhs) = 0._wp
+                  vv(ji,jj,:,Krhs_a) = 0._wp
                   CYCLE
                ENDIF
@@ -1027,16 +1027,16 @@
                   if (vmask(ji,jref,jk) == 0) EXIT
                   N_out = N_out + 1
-                  h_out(N_out) = e3v(ji,jref,jk,Krhs)
+                  h_out(N_out) = e3v(ji,jref,jk,Krhs_a)
                END DO
                IF (N_out == 0) THEN
-                 vv(ji,jj,:,Krhs) = 0._wp
+                 vv(ji,jj,:,Krhs_a) = 0._wp
                  CYCLE
                ENDIF
-               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs),h_out(1:N_out),N_in,N_out)
+               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out)
             END DO
          END DO
 # else
          DO jk = 1, jpkm1
-            vv(i1:i2,j1:j2,jk,Krhs) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs) )
+            vv(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs_a) )
          END DO
 # endif
@@ -1060,5 +1060,5 @@
       !
       IF( before ) THEN 
-         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * uu_b(i1:i2,j1:j2,Kmm)
+         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * uu_b(i1:i2,j1:j2,Kmm_a)
       ELSE
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -1113,5 +1113,5 @@
       ! 
       IF( before ) THEN 
-         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vv_b(i1:i2,j1:j2,Kmm)
+         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vv_b(i1:i2,j1:j2,Kmm_a)
       ELSE
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -1394,5 +1394,5 @@
            DO jj=j1,j2
               DO ji=i1,i2
-                 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w(ji,jj,jk,Kmm) 
+                 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w(ji,jj,jk,Kmm_a) 
               END DO
            END DO
@@ -1415,5 +1415,5 @@
                   IF (wmask(ji,jj,jk) == 0) EXIT 
                   N_out = N_out + 1
-                  h_out(jk) = e3t(ji,jj,jk,Kmm)
+                  h_out(jk) = e3t(ji,jj,jk,Kmm_a)
                ENDDO
                IF (N_in > 0) THEN
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_sponge.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_sponge.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_sponge.F90	(revision 11027)
@@ -191,5 +191,5 @@
    END SUBROUTINE Agrif_Sponge
 
-   SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
+   SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)
       !!----------------------------------------------------------------------
       !!                 *** ROUTINE interptsn_sponge ***
@@ -218,5 +218,5 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb)
+                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a)
                   END DO
                END DO
@@ -228,5 +228,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) 
+                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
                END DO
             END DO
@@ -251,5 +251,5 @@
                   IF (tmask(ji,jj,jk) == 0) EXIT 
                   N_out = N_out + 1
-                  h_out(jk) = e3t(ji,jj,jk,Kmm) !Child grid scale factors. Could multiply by e1e2t here instead of division above
+                  h_out(jk) = e3t(ji,jj,jk,Kmm_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above
                ENDDO
                IF (N_in > 0) THEN
@@ -268,7 +268,7 @@
                DO jk=1,jpkm1
 # if defined key_vertical
-                  tsbdiff(ji,jj,jk,1:jpts) = ts(ji,jj,jk,1:jpts,Kbb) - tabres_child(ji,jj,jk,1:jpts)
+                  tsbdiff(ji,jj,jk,1:jpts) = ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)
 # else
-                  tsbdiff(ji,jj,jk,1:jpts) = ts(ji,jj,jk,1:jpts,Kbb) - tabres(ji,jj,jk,1:jpts)
+                  tsbdiff(ji,jj,jk,1:jpts) = ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts)
 # endif
                ENDDO
@@ -281,5 +281,5 @@
                DO jj = j1,j2
                   DO ji = i1,i2-1
-                     zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)
+                     zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
                      ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
                   END DO
@@ -288,5 +288,5 @@
                DO ji = i1,i2
                   DO jj = j1,j2-1
-                     zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
+                     zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
                      ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
                   END DO
@@ -310,9 +310,9 @@
                   DO ji = i1+1,i2-1
                      IF (.NOT. tabspongedone_tsn(ji,jj)) THEN 
-                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
                         ! horizontal diffusive trends
                         ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
                         ! add it to the general tracer trends
-                        ts(ji,jj,jk,jn,Krhs) = ts(ji,jj,jk,jn,Krhs) + ztsa
+                        ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa
                      ENDIF
                   END DO
@@ -353,7 +353,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb)
+                  tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a)
 # if defined key_vertical
-                  tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kmm)*umask(ji,jj,jk)
+                  tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kmm_a)*umask(ji,jj,jk)
 # endif
                END DO
@@ -384,5 +384,5 @@
                  if (umask(ji,jj,jk) == 0) EXIT
                  N_out = N_out + 1
-                 h_out(N_out) = e3u(ji,jj,jk,Kmm)
+                 h_out(N_out) = e3u(ji,jj,jk,Kmm_a)
               ENDDO
          
@@ -403,7 +403,7 @@
          ENDDO
 
-         ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
+         ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
 #else
-         ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
+         ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
 #endif
          !
@@ -416,7 +416,7 @@
             DO jj = j1,j2
                DO ji = i1+1,i2   ! vector opt.
-                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * fsahm_spt(ji,jj)
-                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm) * ubdiff(ji  ,jj,jk) &
-                                     &   -e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm) * ubdiff(ji-1,jj,jk) ) * zbtr
+                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a) * fsahm_spt(ji,jj)
+                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm_a) * ubdiff(ji  ,jj,jk) &
+                                     &   -e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm_a) * ubdiff(ji-1,jj,jk) ) * zbtr
                END DO
             END DO
@@ -439,9 +439,9 @@
                      ze1v = hdivdiff(ji,jj,jk)
                      ! horizontal diffusive trends
-                     zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm) )   &
+                     zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )   &
                            + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj)
 
                      ! add it to the general momentum trends
-                     uu(ji,jj,jk,Krhs) = uu(ji,jj,jk,Krhs) + zua
+                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua
 
                   END DO
@@ -465,9 +465,9 @@
 
                      ! horizontal diffusive trends
-                     zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm) )   &
+                     zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) )   &
                            + ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
 
                      ! add it to the general momentum trends
-                     vv(ji,jj,jk,Krhs) = vv(ji,jj,jk,Krhs) + zva
+                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) + zva
                   END DO
                ENDIF
@@ -506,7 +506,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb)
+                  tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a)
 # if defined key_vertical
-                  tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kmm)
+                  tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kmm_a)
 # endif
                END DO
@@ -536,5 +536,5 @@
                  if (vmask(ji,jj,jk) == 0) EXIT
                  N_out = N_out + 1
-                 h_out(N_out) = e3v(ji,jj,jk,Kmm)
+                 h_out(N_out) = e3v(ji,jj,jk,Kmm_a)
               ENDDO
          
@@ -549,7 +549,7 @@
          ENDDO
 
-         vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)  
+         vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)  
 # else
-         vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
+         vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
 # endif
          !
@@ -562,7 +562,7 @@
             DO jj = j1+1,j2
                DO ji = i1,i2   ! vector opt.
-                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * fsahm_spt(ji,jj)
-                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vbdiff(ji,jj  ,jk)  &
-                                     &  -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vbdiff(ji,jj-1,jk)  ) * zbtr
+                  zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a) * fsahm_spt(ji,jj)
+                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm_a) * vbdiff(ji,jj  ,jk)  &
+                                     &  -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm_a) * vbdiff(ji,jj-1,jk)  ) * zbtr
                END DO
             END DO
@@ -586,6 +586,6 @@
                IF( .NOT. tabspongedone_u(ji,jj) ) THEN
                   DO jk = 1, jpkm1
-                     uu(ji,jj,jk,Krhs) = uu(ji,jj,jk,Krhs)                                                               &
-                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm) )  &
+                     uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a)                                                               &
+                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) )  &
                         & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj)
                   END DO
@@ -600,6 +600,6 @@
                IF( .NOT. tabspongedone_v(ji,jj) ) THEN
                   DO jk = 1, jpkm1
-                     vv(ji,jj,jk,Krhs) = vv(ji,jj,jk,Krhs)                                                                  &
-                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm) )   &
+                     vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a)                                                                  &
+                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) )   &
                         &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)
                   END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_update.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_update.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_update.F90	(revision 11027)
@@ -230,8 +230,8 @@
       ! -----------------------
       !
-      e3u(:,:,:,Krhs) = e3u(:,:,:,Kmm)
-      e3v(:,:,:,Krhs) = e3v(:,:,:,Kmm)
-!      uu(:,:,:,Krhs) = e3u(:,:,:,Kbb)
-!      vv(:,:,:,Krhs) = e3v(:,:,:,Kbb)
+      e3u(:,:,:,Krhs_a) = e3u(:,:,:,Kmm_a)
+      e3v(:,:,:,Krhs_a) = e3v(:,:,:,Kmm_a)
+!      uu(:,:,:,Krhs_a) = e3u(:,:,:,Kbb_a)
+!      vv(:,:,:,Krhs_a) = e3v(:,:,:,Kbb_a)
       hu_a(:,:) = hu_n(:,:)
       hv_a(:,:) = hv_n(:,:)
@@ -242,10 +242,10 @@
          ! Vertical scale factor interpolations
          ! ------------------------------------
-      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm) ,  'U' )
-      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm) ,  'V' )
-      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:) ,  'F' )
-
-      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
-      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3u(:,:,:,Kmm_a) ,  'U' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3v(:,:,:,Kmm_a) ,  'V' )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3f(:,:,:) ,  'F' )
+
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3uw(:,:,:,Kmm_a), 'UW' )
+      CALL dom_vvl_interpol( e3v(:,:,:,Kmm_a), e3vw(:,:,:,Kmm_a), 'VW' )
 
          ! Update total depths:
@@ -254,6 +254,6 @@
       hv_n(:,:) = 0._wp                        ! Ocean depth at V-points
       DO jk = 1, jpkm1
-         hu_n(:,:) = hu_n(:,:) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
-         hv_n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
+         hu_n(:,:) = hu_n(:,:) + e3u(:,:,jk,Kmm_a) * umask(:,:,jk)
+         hv_n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm_a) * vmask(:,:,jk)
       END DO
       !                                        ! Inverse of the local depth
@@ -268,9 +268,9 @@
          ! Vertical scale factor interpolations
          ! ------------------------------------
-         CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb),  'U'  )
-         CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb),  'V'  )
-
-         CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
-         CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
+         CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3u(:,:,:,Kbb_a),  'U'  )
+         CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3v(:,:,:,Kbb_a),  'V'  )
+
+         CALL dom_vvl_interpol( e3u(:,:,:,Kbb_a), e3uw(:,:,:,Kbb_a), 'UW' )
+         CALL dom_vvl_interpol( e3v(:,:,:,Kbb_a), e3vw(:,:,:,Kbb_a), 'VW' )
 
          ! Update total depths:
@@ -279,6 +279,6 @@
          hv_b(:,:) = 0._wp                     ! Ocean depth at V-points
          DO jk = 1, jpkm1
-            hu_b(:,:) = hu_b(:,:) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
-            hv_b(:,:) = hv_b(:,:) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
+            hu_b(:,:) = hu_b(:,:) + e3u(:,:,jk,Kbb_a) * umask(:,:,jk)
+            hv_b(:,:) = hv_b(:,:) + e3v(:,:,jk,Kbb_a) * vmask(:,:,jk)
          END DO
          !                                     ! Inverse of the local depth
@@ -315,5 +315,5 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm) * e3t(ji,jj,jk,Kmm) ) &
+                     tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) &
                                            * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp
                   END DO
@@ -324,5 +324,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) &
+                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) &
                                            + (tmask(ji,jj,jk)-1)*999._wp
                END DO
@@ -345,5 +345,5 @@
                   IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF
                   N_out = N_out + 1
-                  h_out(N_out) = e3t(ji,jj,jk,Kmm) 
+                  h_out(N_out) = e3t(ji,jj,jk,Kmm_a) 
                ENDDO
                IF (N_in > 0) THEN !Remove this?
@@ -369,7 +369,7 @@
                      DO ji=i1,i2
                         IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
-                           ts(ji,jj,jk,jn,Kbb) = ts(ji,jj,jk,jn,Kbb) & 
+                           ts(ji,jj,jk,jn,Kbb_a) = ts(ji,jj,jk,jn,Kbb_a) & 
                                  & + atfp * ( tabres_child(ji,jj,jk,jn) &
-                                 &          - ts(ji,jj,jk,jn,Kmm) ) * tmask(ji,jj,jk)
+                                 &          - ts(ji,jj,jk,jn,Kmm_a) ) * tmask(ji,jj,jk)
                         ENDIF
                      ENDDO
@@ -383,5 +383,5 @@
                   DO ji=i1,i2
                      IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
-                        ts(ji,jj,jk,jn,Kmm) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
+                        ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
                      END IF
                   END DO
@@ -413,6 +413,6 @@
                   DO ji=i1,i2
 !> jc tmp
-                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm)  * e3t(ji,jj,jk,Kmm) / e3t_0(ji,jj,jk)
-!                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm)  * e3t(ji,jj,jk,Kmm)
+                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk)
+!                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a)
 !< jc tmp
                   END DO
@@ -434,9 +434,9 @@
                      DO ji = i1, i2
                         IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
-                           ztb  = ts(ji,jj,jk,jn,Kbb) * e3t(ji,jj,jk,Kbb) ! fse3t_b prior update should be used
+                           ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
                            ztnu = tabres(ji,jj,jk,jn)
-                           ztno = ts(ji,jj,jk,jn,Kmm) * e3t(ji,jj,jk,Krhs)
-                           ts(ji,jj,jk,jn,Kbb) = ( ztb + atfp * ( ztnu - ztno) )  & 
-                                     &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb)
+                           ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
+                           ts(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) )  & 
+                                     &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
                         ENDIF
                      END DO
@@ -450,5 +450,5 @@
                   DO ji=i1,i2
                      IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 
-                        ts(ji,jj,jk,jn,Kmm) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm)
+                        ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a)
                      END IF
                   END DO
@@ -458,5 +458,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb)  = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm)
+            ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb_a)  = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm_a)
          ENDIF
          !
@@ -495,7 +495,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm)  &
+                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a)  &
                                        + (umask(ji,jj,jk)-1)*999._wp
-                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)  &
+                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)  &
                                        + (umask(ji,jj,jk)-1)*999._wp
                END DO
@@ -520,5 +520,5 @@
                   IF (umask(ji,jj,jk) == 0) EXIT
                   N_out = N_out + 1
-                  h_out(N_out) = e3u(ji,jj,jk,Kmm)
+                  h_out(N_out) = e3u(ji,jj,jk,Kmm_a)
                ENDDO
                IF (N_in * N_out > 0) THEN
@@ -550,9 +550,9 @@
                DO ji=i1,i2
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     uu(ji,jj,jk,Kbb) = uu(ji,jj,jk,Kbb) & 
-                           & + atfp * ( tabres_child(ji,jj,jk) - uu(ji,jj,jk,Kmm) ) * umask(ji,jj,jk)
+                     uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) & 
+                           & + atfp * ( tabres_child(ji,jj,jk) - uu(ji,jj,jk,Kmm_a) ) * umask(ji,jj,jk)
                   ENDIF
                   !
-                  uu(ji,jj,jk,Kmm) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
+                  uu(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
                END DO
             END DO
@@ -579,5 +579,5 @@
          zrhoy = Agrif_Rhoy()
          DO jk = k1, k2
-            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm) * uu(i1:i2,j1:j2,jk,Kmm)
+            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) * uu(i1:i2,j1:j2,jk,Kmm_a)
          END DO
       ELSE
@@ -588,12 +588,12 @@
                   !
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     zub  = uu(ji,jj,jk,Kbb) * e3u(ji,jj,jk,Kbb)  ! fse3t_b prior update should be used
-                     zuno = uu(ji,jj,jk,Kmm) * e3u(ji,jj,jk,Krhs)
+                     zub  = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a)  ! fse3t_b prior update should be used
+                     zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a)
                      zunu = tabres(ji,jj,jk,1)
-                     uu(ji,jj,jk,Kbb) = ( zub + atfp * ( zunu - zuno) ) &      
-                                    & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb)
+                     uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno) ) &      
+                                    & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)
                   ENDIF
                   !
-                  uu(ji,jj,jk,Kmm) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm)
+                  uu(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm_a)
                END DO
             END DO
@@ -601,5 +601,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            uu(i1:i2,j1:j2,k1:k2,Kbb)  = uu(i1:i2,j1:j2,k1:k2,Kmm)
+            uu(i1:i2,j1:j2,k1:k2,Kbb_a)  = uu(i1:i2,j1:j2,k1:k2,Kmm_a)
          ENDIF
          !
@@ -632,8 +632,8 @@
          IF (western_side) THEN
             DO jj=j1,j2
-               zcor = uu_b(i1-1,jj,Kmm) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - uu_b(i1-1,jj,Kmm)
-               uu_b(i1-1,jj,Kmm) = uu_b(i1-1,jj,Kmm) + zcor
+               zcor = uu_b(i1-1,jj,Kmm_a) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - uu_b(i1-1,jj,Kmm_a)
+               uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor
                DO jk=1,jpkm1
-                  uu(i1-1,jj,jk,Kmm) = uu(i1-1,jj,jk,Kmm) + zcor * umask(i1-1,jj,jk)
+                  uu(i1-1,jj,jk,Kmm_a) = uu(i1-1,jj,jk,Kmm_a) + zcor * umask(i1-1,jj,jk)
                END DO 
             END DO
@@ -642,8 +642,8 @@
          IF (eastern_side) THEN
             DO jj=j1,j2
-               zcor = uu_b(i2+1,jj,Kmm) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - uu_b(i2+1,jj,Kmm)
-               uu_b(i2+1,jj,Kmm) = uu_b(i2+1,jj,Kmm) + zcor
+               zcor = uu_b(i2+1,jj,Kmm_a) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - uu_b(i2+1,jj,Kmm_a)
+               uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor
                DO jk=1,jpkm1
-                  uu(i2+1,jj,jk,Kmm) = uu(i2+1,jj,jk,Kmm) + zcor * umask(i2+1,jj,jk)
+                  uu(i2+1,jj,jk,Kmm_a) = uu(i2+1,jj,jk,Kmm_a) + zcor * umask(i2+1,jj,jk)
                END DO 
             END DO
@@ -682,7 +682,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm) &
+                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) &
                                        + (vmask(ji,jj,jk)-1)*999._wp
-                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) &
+                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) &
                                        + (vmask(ji,jj,jk)-1)*999._wp
                END DO
@@ -705,5 +705,5 @@
                   IF (vmask(ji,jj,jk) == 0) EXIT
                   N_out = N_out + 1
-                  h_out(N_out) = e3v(ji,jj,jk,Kmm)
+                  h_out(N_out) = e3v(ji,jj,jk,Kmm_a)
                ENDDO
                IF (N_in * N_out > 0) THEN
@@ -736,9 +736,9 @@
                   !
                   IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part
-                     vv(ji,jj,jk,Kbb) = vv(ji,jj,jk,Kbb) & 
-                           & + atfp * ( tabres_child(ji,jj,jk) - vv(ji,jj,jk,Kmm) ) * vmask(ji,jj,jk)
+                     vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) & 
+                           & + atfp * ( tabres_child(ji,jj,jk) - vv(ji,jj,jk,Kmm_a) ) * vmask(ji,jj,jk)
                   ENDIF
                   !
-                  vv(ji,jj,jk,Kmm) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
+                  vv(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
                END DO
             END DO
@@ -767,5 +767,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm)
+                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)
                END DO
             END DO
@@ -778,12 +778,12 @@
                   !
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     zvb  = vv(ji,jj,jk,Kbb) * e3v(ji,jj,jk,Kbb) ! fse3t_b prior update should be used
-                     zvno = vv(ji,jj,jk,Kmm) * e3v(ji,jj,jk,Krhs)
+                     zvb  = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
+                     zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a)
                      zvnu = tabres(ji,jj,jk,1)
-                     vv(ji,jj,jk,Kbb) = ( zvb + atfp * ( zvnu - zvno) ) &      
-                                    & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb)
+                     vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) ) &      
+                                    & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)
                   ENDIF
                   !
-                  vv(ji,jj,jk,Kmm) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm)
+                  vv(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm_a)
                END DO
             END DO
@@ -791,5 +791,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            vv(i1:i2,j1:j2,k1:k2,Kbb)  = vv(i1:i2,j1:j2,k1:k2,Kmm)
+            vv(i1:i2,j1:j2,k1:k2,Kbb_a)  = vv(i1:i2,j1:j2,k1:k2,Kmm_a)
          ENDIF
          !
@@ -822,8 +822,8 @@
          IF (southern_side) THEN
             DO ji=i1,i2
-               zcor = vv_b(ji,j1-1,Kmm) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vv_b(ji,j1-1,Kmm)
-               vv_b(ji,j1-1,Kmm) = vv_b(ji,j1-1,Kmm) + zcor
+               zcor = vv_b(ji,j1-1,Kmm_a) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vv_b(ji,j1-1,Kmm_a)
+               vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor
                DO jk=1,jpkm1
-                  vv(ji,j1-1,jk,Kmm) = vv(ji,j1-1,jk,Kmm) + zcor * vmask(ji,j1-1,jk)
+                  vv(ji,j1-1,jk,Kmm_a) = vv(ji,j1-1,jk,Kmm_a) + zcor * vmask(ji,j1-1,jk)
                END DO 
             END DO
@@ -832,8 +832,8 @@
          IF (northern_side) THEN
             DO ji=i1,i2
-               zcor = vv_b(ji,j2+1,Kmm) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vv_b(ji,j2+1,Kmm)
-               vv_b(ji,j2+1,Kmm) = vv_b(ji,j2+1,Kmm) + zcor
+               zcor = vv_b(ji,j2+1,Kmm_a) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vv_b(ji,j2+1,Kmm_a)
+               vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor
                DO jk=1,jpkm1
-                  vv(ji,j2+1,jk,Kmm) = vv(ji,j2+1,jk,Kmm) + zcor * vmask(ji,j2+1,jk)
+                  vv(ji,j2+1,jk,Kmm_a) = vv(ji,j2+1,jk,Kmm_a) + zcor * vmask(ji,j2+1,jk)
                END DO 
             END DO
@@ -862,5 +862,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm) * hu_n(ji,jj) * e2u(ji,jj)
+               tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm_a) * hu_n(ji,jj) * e2u(ji,jj)
             END DO
          END DO
@@ -873,10 +873,10 @@
                spgu(ji,jj) = 0._wp
                DO jk=1,jpkm1
-                  spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm)
+                  spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)
                END DO
                !
                zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj)
                DO jk=1,jpkm1              
-                  uu(ji,jj,jk,Kmm) = uu(ji,jj,jk,Kmm) + zcorr * umask(ji,jj,jk)           
+                  uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk)           
                END DO
                !
@@ -884,19 +884,19 @@
                IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
-                     uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + atfp * zcorr * umask(ji,jj,1)
+                     zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
+                     uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + atfp * zcorr * umask(ji,jj,1)
                   END IF
                ENDIF    
-               uu_b(ji,jj,Kmm) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
+               uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
                !       
                ! Correct "before" velocities to hold correct bt component:
                spgu(ji,jj) = 0.e0
                DO jk=1,jpkm1
-                  spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb)
+                  spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a)
                END DO
                !
-               zcorr = uu_b(ji,jj,Kbb) - spgu(ji,jj) * r1_hu_b(ji,jj)
+               zcorr = uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu_b(ji,jj)
                DO jk=1,jpkm1              
-                  uu(ji,jj,jk,Kbb) = uu(ji,jj,jk,Kbb) + zcorr * umask(ji,jj,jk)           
+                  uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk)           
                END DO
                !
@@ -905,5 +905,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            uu_b(i1:i2,j1:j2,Kbb)  = uu_b(i1:i2,j1:j2,Kmm)
+            uu_b(i1:i2,j1:j2,Kbb_a)  = uu_b(i1:i2,j1:j2,Kmm_a)
          ENDIF
       ENDIF
@@ -928,5 +928,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm) * hv_n(ji,jj) * e1v(ji,jj) 
+               tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm_a) * hv_n(ji,jj) * e1v(ji,jj) 
             END DO
          END DO
@@ -939,10 +939,10 @@
                spgv(ji,jj) = 0.e0
                DO jk=1,jpkm1
-                  spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm)
+                  spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)
                END DO
                !
                zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj)
                DO jk=1,jpkm1              
-                  vv(ji,jj,jk,Kmm) = vv(ji,jj,jk,Kmm) + zcorr * vmask(ji,jj,jk)           
+                  vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk)           
                END DO
                !
@@ -950,19 +950,19 @@
                IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
-                     vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + atfp * zcorr * vmask(ji,jj,1)
+                     zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
+                     vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + atfp * zcorr * vmask(ji,jj,1)
                   END IF
                ENDIF              
-               vv_b(ji,jj,Kmm) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
+               vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
                !       
                ! Correct "before" velocities to hold correct bt component:
                spgv(ji,jj) = 0.e0
                DO jk=1,jpkm1
-                  spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb)
+                  spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a)
                END DO
                !
-               zcorr = vv_b(ji,jj,Kbb) - spgv(ji,jj) * r1_hv_b(ji,jj)
+               zcorr = vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv_b(ji,jj)
                DO jk=1,jpkm1              
-                  vv(ji,jj,jk,Kbb) = vv(ji,jj,jk,Kbb) + zcorr * vmask(ji,jj,jk)           
+                  vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk)           
                END DO
                !
@@ -971,5 +971,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            vv_b(i1:i2,j1:j2,Kbb)  = vv_b(i1:i2,j1:j2,Kmm)
+            vv_b(i1:i2,j1:j2,Kbb_a)  = vv_b(i1:i2,j1:j2,Kmm_a)
          ENDIF
          !
@@ -993,5 +993,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) = ssh(ji,jj,Kmm)
+               tabres(ji,jj) = ssh(ji,jj,Kmm_a)
             END DO
          END DO
@@ -1000,6 +1000,6 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ssh(ji,jj,Kbb) =   ssh(ji,jj,Kbb) &
-                        & + atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm) ) * tmask(ji,jj,1)
+                  ssh(ji,jj,Kbb_a) =   ssh(ji,jj,Kbb_a) &
+                        & + atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) * tmask(ji,jj,1)
                END DO
             END DO
@@ -1008,10 +1008,10 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               ssh(ji,jj,Kmm) = tabres(ji,jj) * tmask(ji,jj,1)
+               ssh(ji,jj,Kmm_a) = tabres(ji,jj) * tmask(ji,jj,1)
             END DO
          END DO
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            ssh(i1:i2,j1:j2,Kbb)  = ssh(i1:i2,j1:j2,Kmm)
+            ssh(i1:i2,j1:j2,Kbb_a)  = ssh(i1:i2,j1:j2,Kmm_a)
          ENDIF
          !
@@ -1094,6 +1094,6 @@
             DO jj=j1,j2
                zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 
-               ssh(i1  ,jj,Kmm) = ssh(i1  ,jj,Kmm) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i1  ,jj,Kbb) = ssh(i1  ,jj,Kbb) + atfp * zcor
+               ssh(i1  ,jj,Kmm_a) = ssh(i1  ,jj,Kmm_a) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i1  ,jj,Kbb_a) = ssh(i1  ,jj,Kbb_a) + atfp * zcor
             END DO
          ENDIF
@@ -1101,6 +1101,6 @@
             DO jj=j1,j2
                zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
-               ssh(i2+1,jj,Kmm) = ssh(i2+1,jj,Kmm) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i2+1,jj,Kbb) = ssh(i2+1,jj,Kbb) + atfp * zcor
+               ssh(i2+1,jj,Kmm_a) = ssh(i2+1,jj,Kmm_a) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + atfp * zcor
             END DO
          ENDIF
@@ -1182,6 +1182,6 @@
             DO ji=i1,i2
                zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
-               ssh(ji,j1  ,Kmm) = ssh(ji,j1  ,Kmm) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j1  ,Kbb) = ssh(ji,j1,Kbb) + atfp * zcor
+               ssh(ji,j1  ,Kmm_a) = ssh(ji,j1  ,Kmm_a) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j1  ,Kbb_a) = ssh(ji,j1,Kbb_a) + atfp * zcor
             END DO
          ENDIF
@@ -1189,6 +1189,6 @@
             DO ji=i1,i2
                zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
-               ssh(ji,j2+1,Kmm) = ssh(ji,j2+1,Kmm) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j2+1,Kbb) = ssh(ji,j2+1,Kbb) + atfp * zcor
+               ssh(ji,j2+1,Kmm_a) = ssh(ji,j2+1,Kmm_a) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + atfp * zcor
             END DO
          ENDIF
@@ -1319,5 +1319,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Kmm) &
+                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Kmm_a) &
                                      & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
                END DO
@@ -1330,8 +1330,8 @@
          ! Save "old" scale factor (prior update) for subsequent asselin correction
          ! of prognostic variables
-         e3t(i1:i2,j1:j2,1:jpkm1,Krhs) = e3t(i1:i2,j1:j2,1:jpkm1,Kmm)
-
-         ! One should also save e3t(:,:,:,Kbb), but lacking of workspace...
-!         hdiv(i1:i2,j1:j2,1:jpkm1)   = e3t(i1:i2,j1:j2,1:jpkm1,Kbb)
+         e3t(i1:i2,j1:j2,1:jpkm1,Krhs_a) = e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a)
+
+         ! One should also save e3t(:,:,:,Kbb_a), but lacking of workspace...
+!         hdiv(i1:i2,j1:j2,1:jpkm1)   = e3t(i1:i2,j1:j2,1:jpkm1,Kbb_a)
 
          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
@@ -1339,13 +1339,13 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     e3t(ji,jj,jk,Kbb) =  e3t(ji,jj,jk,Kbb) &
-                           & + atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm) )
+                     e3t(ji,jj,jk,Kbb_a) =  e3t(ji,jj,jk,Kbb_a) &
+                           & + atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) )
                   END DO
                END DO
             END DO
             !
-            e3w  (i1:i2,j1:j2,1,Kbb) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kbb) - e3t_0(i1:i2,j1:j2,1)
-            gdepw(i1:i2,j1:j2,1,Kbb) = 0.0_wp
-            gdept(i1:i2,j1:j2,1,Kbb) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kbb)
+            e3w  (i1:i2,j1:j2,1,Kbb_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kbb_a) - e3t_0(i1:i2,j1:j2,1)
+            gdepw(i1:i2,j1:j2,1,Kbb_a) = 0.0_wp
+            gdept(i1:i2,j1:j2,1,Kbb_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kbb_a)
             !
             DO jk = 2, jpk
@@ -1353,11 +1353,11 @@
                   DO ji = i1,i2            
                      zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-                     e3w(ji,jj,jk,Kbb)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
-                     &                                        ( e3t(ji,jj,jk-1,Kbb) - e3t_0(ji,jj,jk-1) )  &
+                     e3w(ji,jj,jk,Kbb_a)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
+                     &                                        ( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) )  &
                      &                                  +            0.5_wp * tmask(ji,jj,jk)   *        &
-                     &                                        ( e3t(ji,jj,jk  ,Kbb) - e3t_0(ji,jj,jk  ) )
-                     gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb)
-                     gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  &
-                         &               + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb)) 
+                     &                                        ( e3t(ji,jj,jk  ,Kbb_a) - e3t_0(ji,jj,jk  ) )
+                     gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk-1,Kbb_a)
+                     gdept(ji,jj,jk,Kbb_a) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb_a) + 0.5 * e3w(ji,jj,jk,Kbb_a))  &
+                         &               + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) +       e3w(ji,jj,jk,Kbb_a)) 
                   END DO
                END DO
@@ -1370,17 +1370,17 @@
          !
          ! Update vertical scale factor at T-points:
-         e3t(i1:i2,j1:j2,1:jpkm1,Kmm) = ptab(i1:i2,j1:j2,1:jpkm1)
+         e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a) = ptab(i1:i2,j1:j2,1:jpkm1)
          !
          ! Update total depth:
          ht_n(i1:i2,j1:j2) = 0._wp
          DO jk = 1, jpkm1
-            ht_n(i1:i2,j1:j2) = ht_n(i1:i2,j1:j2) + e3t(i1:i2,j1:j2,jk,Kmm) * tmask(i1:i2,j1:j2,jk)
+            ht_n(i1:i2,j1:j2) = ht_n(i1:i2,j1:j2) + e3t(i1:i2,j1:j2,jk,Kmm_a) * tmask(i1:i2,j1:j2,jk)
          END DO
          !
          ! Update vertical scale factor at W-points and depths:
-         e3w (i1:i2,j1:j2,1,Kmm) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm) - e3t_0(i1:i2,j1:j2,1)
-         gdept(i1:i2,j1:j2,1,Kmm) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm)
-         gdepw(i1:i2,j1:j2,1,Kmm) = 0.0_wp
-         gde3w(i1:i2,j1:j2,1) = gdept(i1:i2,j1:j2,1,Kmm) - (ht_n(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh
+         e3w (i1:i2,j1:j2,1,Kmm_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm_a) - e3t_0(i1:i2,j1:j2,1)
+         gdept(i1:i2,j1:j2,1,Kmm_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm_a)
+         gdepw(i1:i2,j1:j2,1,Kmm_a) = 0.0_wp
+         gde3w(i1:i2,j1:j2,1) = gdept(i1:i2,j1:j2,1,Kmm_a) - (ht_n(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh
          !
          DO jk = 2, jpk
@@ -1388,10 +1388,10 @@
                DO ji = i1,i2            
                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               e3w(ji,jj,jk,Kmm)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t(ji,jj,jk-1,Kmm) - e3t_0(ji,jj,jk-1) )   &
-               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t(ji,jj,jk  ,Kmm) - e3t_0(ji,jj,jk  ) )
-               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
-               gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
-                   &               + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm)) 
-               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - (ht_n(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
+               e3w(ji,jj,jk,Kmm_a)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t(ji,jj,jk-1,Kmm_a) - e3t_0(ji,jj,jk-1) )   &
+               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t(ji,jj,jk  ,Kmm_a) - e3t_0(ji,jj,jk  ) )
+               gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a)
+               gdept(ji,jj,jk,Kmm_a) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm_a) + 0.5 * e3w(ji,jj,jk,Kmm_a))  &
+                   &               + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm_a) +       e3w(ji,jj,jk,Kmm_a)) 
+               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm_a) - (ht_n(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
                END DO
             END DO
@@ -1399,8 +1399,8 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            e3t (i1:i2,j1:j2,1:jpk,Kbb)  = e3t (i1:i2,j1:j2,1:jpk,Kmm)
-            e3w (i1:i2,j1:j2,1:jpk,Kbb)  = e3w (i1:i2,j1:j2,1:jpk,Kmm)
-            gdepw(i1:i2,j1:j2,1:jpk,Kbb) = gdepw(i1:i2,j1:j2,1:jpk,Kmm)
-            gdept(i1:i2,j1:j2,1:jpk,Kbb) = gdept(i1:i2,j1:j2,1:jpk,Kmm)
+            e3t (i1:i2,j1:j2,1:jpk,Kbb_a)  = e3t (i1:i2,j1:j2,1:jpk,Kmm_a)
+            e3w (i1:i2,j1:j2,1:jpk,Kbb_a)  = e3w (i1:i2,j1:j2,1:jpk,Kmm_a)
+            gdepw(i1:i2,j1:j2,1:jpk,Kbb_a) = gdepw(i1:i2,j1:j2,1:jpk,Kmm_a)
+            gdept(i1:i2,j1:j2,1:jpk,Kbb_a) = gdept(i1:i2,j1:j2,1:jpk,Kmm_a)
          ENDIF
          !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_interp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_interp.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_interp.F90	(revision 11027)
@@ -73,5 +73,5 @@
                DO jj=j1,j2
                  DO ji=i1,i2
-                       ptab(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm)
+                       ptab(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a)
                  END DO
               END DO
@@ -83,5 +83,5 @@
            DO jj=j1,j2
               DO ji=i1,i2
-                 ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) 
+                 ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
               END DO
            END DO
@@ -113,5 +113,5 @@
                   IF (tmask(iref,jref,jk) == 0) EXIT 
                   N_out = N_out + 1
-                  h_out(jk) = e3t(iref,jref,jk,Kmm)
+                  h_out(jk) = e3t(iref,jref,jk,Kmm_a)
                ENDDO
                IF (N_in > 0) THEN
@@ -127,5 +127,5 @@
          !
          DO jn=1, jptra
-            tr(i1:i2,j1:j2,1:jpk,jn,Krhs)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
+            tr(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
          END DO
 
@@ -151,20 +151,20 @@
                ibdy = nlci-nbghostcells
                DO jn = 1, jptra
-                  tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
+                  tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
                   DO jk = 1, jpkm1
                      DO jj = jmin,jmax
                         IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN
-                           tr(ibdy,jj,jk,jn,Krhs) = tr(ibdy+1,jj,jk,jn,Krhs) * tmask(ibdy,jj,jk)
-                        ELSE
-                           tr(ibdy,jj,jk,jn,Krhs)=(z4*tr(ibdy+1,jj,jk,jn,Krhs)+z3*tr(ibdy-1,jj,jk,jn,Krhs))*tmask(ibdy,jj,jk)
-                           IF( uu(ibdy-1,jj,jk,Kmm) > 0._wp ) THEN
-                              tr(ibdy,jj,jk,jn,Krhs)=( z6*tr(ibdy-1,jj,jk,jn,Krhs)+z5*tr(ibdy+1,jj,jk,jn,Krhs) & 
-                                                 + z7*tr(ibdy-2,jj,jk,jn,Krhs) ) * tmask(ibdy,jj,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy+1,jmin:jmax,1:jpkm1)
+                           tr(ibdy,jj,jk,jn,Krhs_a) = tr(ibdy+1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)
+                        ELSE
+                           tr(ibdy,jj,jk,jn,Krhs_a)=(z4*tr(ibdy+1,jj,jk,jn,Krhs_a)+z3*tr(ibdy-1,jj,jk,jn,Krhs_a))*tmask(ibdy,jj,jk)
+                           IF( uu(ibdy-1,jj,jk,Kmm_a) > 0._wp ) THEN
+                              tr(ibdy,jj,jk,jn,Krhs_a)=( z6*tr(ibdy-1,jj,jk,jn,Krhs_a)+z5*tr(ibdy+1,jj,jk,jn,Krhs_a) & 
+                                                 + z7*tr(ibdy-2,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)
+                           ENDIF
+                        ENDIF
+                     END DO
+                  END DO
+                  ! Restore ghost points:
+                  tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy+1,jmin:jmax,1:jpkm1)
                END DO
             ENDIF
@@ -180,20 +180,20 @@
                jbdy = nlcj-nbghostcells         
                DO jn = 1, jptra
-                  tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
+                  tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
                   DO jk = 1, jpkm1
                      DO ji = imin,imax
                         IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN
-                           tr(ji,jbdy,jk,jn,Krhs) = tr(ji,jbdy+1,jk,jn,Krhs) * tmask(ji,jbdy,jk)
-                        ELSE
-                           tr(ji,jbdy,jk,jn,Krhs)=(z4*tr(ji,jbdy+1,jk,jn,Krhs)+z3*tr(ji,jbdy-1,jk,jn,Krhs))*tmask(ji,jbdy,jk)        
-                           IF (vv(ji,jbdy-1,jk,Kmm) > 0._wp ) THEN
-                              tr(ji,jbdy,jk,jn,Krhs)=( z6*tr(ji,jbdy-1,jk,jn,Krhs)+z5*tr(ji,jbdy+1,jk,jn,Krhs)  &
-                                                 + z7*tr(ji,jbdy-2,jk,jn,Krhs) ) * tmask(ji,jbdy,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) * tmask(imin:imax,jbdy+1,1:jpkm1)
+                           tr(ji,jbdy,jk,jn,Krhs_a) = tr(ji,jbdy+1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
+                        ELSE
+                           tr(ji,jbdy,jk,jn,Krhs_a)=(z4*tr(ji,jbdy+1,jk,jn,Krhs_a)+z3*tr(ji,jbdy-1,jk,jn,Krhs_a))*tmask(ji,jbdy,jk)        
+                           IF (vv(ji,jbdy-1,jk,Kmm_a) > 0._wp ) THEN
+                              tr(ji,jbdy,jk,jn,Krhs_a)=( z6*tr(ji,jbdy-1,jk,jn,Krhs_a)+z5*tr(ji,jbdy+1,jk,jn,Krhs_a)  &
+                                                 + z7*tr(ji,jbdy-2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
+                           ENDIF
+                        ENDIF
+                     END DO
+                  END DO
+                  ! Restore ghost points:
+                  tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) * tmask(imin:imax,jbdy+1,1:jpkm1)
                END DO
             ENDIF
@@ -209,20 +209,20 @@
                ibdy = 1+nbghostcells       
                DO jn = 1, jptra
-                  tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
+                  tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
                   DO jk = 1, jpkm1
                      DO jj = jmin,jmax
                         IF( umask(ibdy,jj,jk) == 0._wp ) THEN
-                           tr(ibdy,jj,jk,jn,Krhs) = tr(ibdy-1,jj,jk,jn,Krhs) * tmask(ibdy,jj,jk)
-                        ELSE
-                           tr(ibdy,jj,jk,jn,Krhs)=(z4*tr(ibdy-1,jj,jk,jn,Krhs)+z3*tr(ibdy+1,jj,jk,jn,Krhs))*tmask(ibdy,jj,jk)        
-                           IF( uu(ibdy,jj,jk,Kmm) < 0._wp ) THEN
-                              tr(ibdy,jj,jk,jn,Krhs)=( z6*tr(ibdy+1,jj,jk,jn,Krhs)+z5*tr(ibdy-1,jj,jk,jn,Krhs) &
-                                                 + z7*tr(ibdy+2,jj,jk,jn,Krhs) ) * tmask(ibdy,jj,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy-1,jmin:jmax,1:jpkm1)
+                           tr(ibdy,jj,jk,jn,Krhs_a) = tr(ibdy-1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)
+                        ELSE
+                           tr(ibdy,jj,jk,jn,Krhs_a)=(z4*tr(ibdy-1,jj,jk,jn,Krhs_a)+z3*tr(ibdy+1,jj,jk,jn,Krhs_a))*tmask(ibdy,jj,jk)        
+                           IF( uu(ibdy,jj,jk,Kmm_a) < 0._wp ) THEN
+                              tr(ibdy,jj,jk,jn,Krhs_a)=( z6*tr(ibdy+1,jj,jk,jn,Krhs_a)+z5*tr(ibdy-1,jj,jk,jn,Krhs_a) &
+                                                 + z7*tr(ibdy+2,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)
+                           ENDIF
+                        ENDIF
+                     END DO
+                  END DO
+                  ! Restore ghost points:
+                  tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy-1,jmin:jmax,1:jpkm1)
                END DO
             ENDIF
@@ -238,20 +238,20 @@
                jbdy=1+nbghostcells        
                DO jn = 1, jptra
-                  tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
+                  tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
                   DO jk = 1, jpkm1      
                      DO ji = imin,imax
                         IF( vmask(ji,jbdy,jk) == 0._wp ) THEN
-                           tr(ji,jbdy,jk,jn,Krhs)=tr(ji,jbdy-1,jk,jn,Krhs) * tmask(ji,jbdy,jk)
-                        ELSE
-                           tr(ji,jbdy,jk,jn,Krhs)=(z4*tr(ji,jbdy-1,jk,jn,Krhs)+z3*tr(ji,jbdy+1,jk,jn,Krhs))*tmask(ji,jbdy,jk)
-                           IF( vv(ji,jbdy,jk,Kmm) < 0._wp ) THEN
-                              tr(ji,jbdy,jk,jn,Krhs)=( z6*tr(ji,jbdy+1,jk,jn,Krhs)+z5*tr(ji,jbdy-1,jk,jn,Krhs) & 
-                                                 + z7*tr(ji,jbdy+2,jk,jn,Krhs) ) * tmask(ji,jbdy,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) * tmask(imin:imax,jbdy-1,1:jpkm1)
+                           tr(ji,jbdy,jk,jn,Krhs_a)=tr(ji,jbdy-1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
+                        ELSE
+                           tr(ji,jbdy,jk,jn,Krhs_a)=(z4*tr(ji,jbdy-1,jk,jn,Krhs_a)+z3*tr(ji,jbdy+1,jk,jn,Krhs_a))*tmask(ji,jbdy,jk)
+                           IF( vv(ji,jbdy,jk,Kmm_a) < 0._wp ) THEN
+                              tr(ji,jbdy,jk,jn,Krhs_a)=( z6*tr(ji,jbdy+1,jk,jn,Krhs_a)+z5*tr(ji,jbdy-1,jk,jn,Krhs_a) & 
+                                                 + z7*tr(ji,jbdy+2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
+                           ENDIF
+                        ENDIF
+                     END DO
+                  END DO
+                  ! Restore ghost points:
+                  tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) * tmask(imin:imax,jbdy-1,1:jpkm1)
                END DO
             ENDIF
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_user.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_user.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_user.F90	(revision 11027)
@@ -175,6 +175,6 @@
    tabspongedone_tsn = .FALSE.
    CALL Agrif_Bc_variable(tsn_sponge_id,calledweight=1.,procname=interptsn_sponge)
-   ! reset ts(:,:,:,:,Krhs) to zero
-   ts(:,:,:,:,Krhs) = 0.
+   ! reset ts(:,:,:,:,Krhs_a) to zero
+   ts(:,:,:,:,Krhs_a) = 0.
 
    Agrif_UseSpecialValue = ln_spc_dyn
@@ -191,5 +191,5 @@
    CALL Agrif_Bc_variable(sshn_id,calledweight=1., procname=interpsshn )
    hbdy_w(:,:) = 0.e0 ; hbdy_e(:,:) = 0.e0 ; hbdy_n(:,:) = 0.e0 ; hbdy_s(:,:) = 0.e0
-   ssh(:,:,Krhs) = 0.e0
+   ssh(:,:,Krhs_a) = 0.e0
 
    IF ( ln_dynspg_ts ) THEN
@@ -207,6 +207,6 @@
    Agrif_UseSpecialValue = .FALSE. 
    ! reset velocities to zero
-   uu(:,:,:,Krhs) = 0.
-   vv(:,:,:,Krhs) = 0.
+   uu(:,:,:,Krhs_a) = 0.
+   vv(:,:,:,Krhs_a) = 0.
 
    ! 3. Some controls
@@ -591,6 +591,6 @@
    tabspongedone_trn = .FALSE.
    CALL Agrif_Bc_variable(trn_sponge_id,calledweight=1.,procname=interptrn_sponge)
-   ! reset ts(:,:,:,:,Krhs) to zero
-   tr(:,:,:,:,Krhs) = 0.
+   ! reset ts(:,:,:,:,Krhs_a) to zero
+   tr(:,:,:,:,Krhs_a) = 0.
 
 
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diawri.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diawri.F90	(revision 11027)
@@ -138,5 +138,5 @@
 
       IF( ll_wd ) THEN
-         CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
+         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
       ELSE
          CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/dom_oce.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/dom_oce.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/dom_oce.F90	(revision 11027)
@@ -155,17 +155,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
-
-   !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::     e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::     e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::     e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::               e3f_n            !: f- vert. scale factor [m]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::     e3w_b ,   e3w_n            !: w- vert. scale factor [m]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::    e3uw_b ,  e3uw_n            !: uw-vert. scale factor [m]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::    e3vw_b ,  e3vw_n            !: vw-vert. scale factor [m]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::    gdept_b , gdept_n           !: t- depth              [m]    
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::    gdepw_b , gdepw_n           !: w- depth              [m]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:) ::    gde3w_n                     !: w- depth (sum of e3w) [m]
-   !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY
 
    INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90	(revision 11027)
@@ -64,5 +64,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_nxt ( kt, Kbb, Kmm, Kaa )
+   SUBROUTINE dyn_nxt ( kt, Kbb, Kmm, Krhs, puu, pvv, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_nxt  ***
@@ -83,15 +83,16 @@
       !!              * Apply the time filter applied and swap of the dynamics
       !!             arrays to start the next time step:
-      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
-      !!                (un,vn) = (ua,va).
+      !!                (puu(Kbb),pvv(Kbb)) = (puu(Kmm),pvv(Kmm)) + atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Krhs),pvv(Krhs)) - 2 (puu(Kmm),pvv(Kmm)) ]
+      !!                (puu(Kmm),pvv(Kmm)) = (puu(Krhs),pvv(Krhs)).
       !!             Note that with flux form advection and non linear free surface,
       !!             the time filter is applied on thickness weighted velocity.
       !!             As a result, dyn_nxt MUST be called after tra_nxt.
       !!
-      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
-      !!               un,vn   now horizontal velocity of next time-step
+      !! ** Action :   puu(Kbb),pvv(Kbb)   filtered before horizontal velocity of next time-step
+      !!               puu(Kmm),pvv(Kmm)   now horizontal velocity of next time-step
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt             ! ocean time-step index
-      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa  ! time level indices
+      INTEGER, INTENT( in ) ::   kt                   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs, Kaa  ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv        !  velocity arrays
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -116,13 +117,13 @@
          ! Ensure below that barotropic velocities match time splitting estimate
          ! Compute actual transport and replace it with ts estimate at "after" time step
-         zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
-         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
+         zue(:,:) = e3u(:,:,1,Krhs) * puu(:,:,1,Krhs) * umask(:,:,1)
+         zve(:,:) = e3v(:,:,1,Krhs) * pvv(:,:,1,Krhs) * vmask(:,:,1)
          DO jk = 2, jpkm1
-            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
-            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
+            zue(:,:) = zue(:,:) + e3u(:,:,jk,Krhs) * puu(:,:,jk,Krhs) * umask(:,:,jk)
+            zve(:,:) = zve(:,:) + e3v(:,:,jk,Krhs) * pvv(:,:,jk,Krhs) * vmask(:,:,jk)
          END DO
          DO jk = 1, jpkm1
-            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
-            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
+            puu(:,:,jk,Krhs) = ( puu(:,:,jk,Krhs) - zue(:,:) * r1_hu_a(:,:) + uu_b(:,:,Krhs) ) * umask(:,:,jk)
+            pvv(:,:,jk,Krhs) = ( pvv(:,:,jk,Krhs) - zve(:,:) * r1_hv_a(:,:) + vv_b(:,:,Krhs) ) * vmask(:,:,jk)
          END DO
          !
@@ -133,6 +134,6 @@
             ! so that asselin contribution is removed at the same time 
             DO jk = 1, jpkm1
-               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
-               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk)
+               puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu_n(:,:) + uu_b(:,:,Kmm) )*umask(:,:,jk)
+               pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv_n(:,:) + vv_b(:,:,Kmm) )*vmask(:,:,jk)
             END DO  
          ENDIF
@@ -142,8 +143,8 @@
       ! --------------------------------------------------      
 # if defined key_agrif
-      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
+      Krhs_a = Krhs ; CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
 # endif
       !
-      CALL lbc_lnk_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries
+      CALL lbc_lnk_multi( 'dynnxt', puu(:,:,:,Krhs), 'U', -1., pvv(:,:,:,Krhs), 'V', -1. )     !* local domain boundaries
       !
       !                                !* BDY open boundaries
@@ -162,15 +163,15 @@
          !
          !                                  ! Kinetic energy and Conversion
-         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt, Kmm )
+         IF( ln_KE_trd  )   CALL trd_dyn( puu(:,:,:,Krhs), pvv(:,:,:,Krhs), jpdyn_ken, kt, Kmm )
          !
          IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
-            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
-            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
+            zua(:,:,:) = ( puu(:,:,:,Krhs) - puu(:,:,:,Kbb) ) * z1_2dt
+            zva(:,:,:) = ( pvv(:,:,:,Krhs) - pvv(:,:,:,Kbb) ) * z1_2dt
             CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
             CALL iom_put( "vtrd_tot", zva )
          ENDIF
          !
-         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
-         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
+         zua(:,:,:) = puu(:,:,:,Kmm)             ! save the now velocity before the asselin filter
+         zva(:,:,:) = pvv(:,:,:,Kmm)             ! (caution: there will be a shift by 1 timestep in the
          !                                  !  computation of the asselin filter trends)
       ENDIF
@@ -180,17 +181,17 @@
       IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
          DO jk = 1, jpkm1
-            un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua
-            vn(:,:,jk) = va(:,:,jk)
+            puu(:,:,jk,Kmm) = puu(:,:,jk,Krhs)                         ! puu(:,:,:,Kmm) <-- puu(:,:,:,Krhs)
+            pvv(:,:,jk,Kmm) = pvv(:,:,jk,Krhs)
          END DO
          IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
 !!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a  
             DO jk = 1, jpkm1
-!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
-!               e3u_b(:,:,jk) = e3u_n(:,:,jk)
-!               e3v_b(:,:,jk) = e3v_n(:,:,jk)
+!               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
+!               e3u(:,:,jk,Kbb) = e3u(:,:,jk,Kmm)
+!               e3v(:,:,jk,Kbb) = e3v(:,:,jk,Kmm)
                !
-               e3t_n(:,:,jk) = e3t_a(:,:,jk)
-               e3u_n(:,:,jk) = e3u_a(:,:,jk)
-               e3v_n(:,:,jk) = e3v_a(:,:,jk)
+               e3t(:,:,jk,Kmm) = e3t(:,:,jk,Krhs)
+               e3u(:,:,jk,Kmm) = e3u(:,:,jk,Krhs)
+               e3v(:,:,jk,Kmm) = e3v(:,:,jk,Krhs)
             END DO
 !!gm BUG end
@@ -205,11 +206,11 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi    
-                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
-                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
+                     zuf = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Krhs) )
+                     zvf = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Krhs) )
                      !
-                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
-                     vb(ji,jj,jk) = zvf
-                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
-                     vn(ji,jj,jk) = va(ji,jj,jk)
+                     puu(ji,jj,jk,Kbb) = zuf                      ! puu(:,:,:,Kbb) <-- filtered velocity
+                     pvv(ji,jj,jk,Kbb) = zvf
+                     puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Krhs)             ! puu(:,:,:,Kmm) <-- puu(:,:,:,Krhs)
+                     pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Krhs)
                   END DO
                END DO
@@ -222,5 +223,5 @@
             ! ----------------------------------------------------
             DO jk = 1, jpkm1
-               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
+               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm) + atfp * ( e3t(:,:,jk,Kbb) - 2._wp * e3t(:,:,jk,Kmm) + e3t(:,:,jk,Krhs) )
             END DO
             ! Add volume filter correction: compatibility with tracer advection scheme
@@ -228,5 +229,5 @@
             zcoef = atfp * rdt * r1_rau0
 
-            e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+            e3t(:,:,1,Kbb) = e3t(:,:,1,Kbb) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
 
             IF ( ln_rnf ) THEN
@@ -236,6 +237,6 @@
                         DO ji = 1, jpi
                            IF( jk <=  nk_rnf(ji,jj)  ) THEN
-                               e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
-                                      &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
+                               e3t(ji,jj,jk,Kbb) =   e3t(ji,jj,jk,Kbb) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
+                                      &          * ( e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
                            ENDIF
                         ENDDO
@@ -243,5 +244,5 @@
                   ENDDO
                ELSE
-                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
+                  e3t(:,:,1,Kbb) = e3t(:,:,1,Kbb) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
                ENDIF
             END IF
@@ -252,9 +253,9 @@
                      DO ji = 1, jpi
                         IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj)  ) THEN
-                           e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
-                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)
+                           e3t(ji,jj,jk,Kbb) = e3t(ji,jj,jk,Kbb) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
+                                &          * ( e3t(ji,jj,jk,Kmm) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)
                         ELSEIF ( jk==misfkb(ji,jj) ) THEN
-                           e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
-                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk)
+                           e3t(ji,jj,jk,Kbb) = e3t(ji,jj,jk,Kbb) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
+                                &          * ( e3t(ji,jj,jk,Kmm) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk)
                         ENDIF
                      END DO
@@ -265,16 +266,16 @@
             IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
                ! Before filtered scale factor at (u/v)-points
-               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
-               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
+               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )
+               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )
                DO jk = 1, jpkm1
                   DO jj = 1, jpj
                      DO ji = 1, jpi
-                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
-                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
+                        zuf = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Krhs) )
+                        zvf = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Krhs) )
                         !
-                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
-                        vb(ji,jj,jk) = zvf
-                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
-                        vn(ji,jj,jk) = va(ji,jj,jk)
+                        puu(ji,jj,jk,Kbb) = zuf                      ! puu(:,:,:,Kbb) <-- filtered velocity
+                        pvv(ji,jj,jk,Kbb) = zvf
+                        puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Krhs)             ! puu(:,:,:,Kmm) <-- puu(:,:,:,Krhs)
+                        pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Krhs)
                      END DO
                   END DO
@@ -285,28 +286,28 @@
                ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
                ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
-               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
-               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
+               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), ze3u_f, 'U' )
+               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), ze3v_f, 'V' )
                DO jk = 1, jpkm1
                   DO jj = 1, jpj
                      DO ji = 1, jpi                  
-                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
-                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
-                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
-                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
-                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
-                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
+                        zue3a = e3u(ji,jj,jk,Krhs) * puu(ji,jj,jk,Krhs)
+                        zve3a = e3v(ji,jj,jk,Krhs) * pvv(ji,jj,jk,Krhs)
+                        zue3n = e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
+                        zve3n = e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
+                        zue3b = e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb)
+                        zve3b = e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb)
                         !
                         zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
                         zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
                         !
-                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
-                        vb(ji,jj,jk) = zvf
-                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
-                        vn(ji,jj,jk) = va(ji,jj,jk)
+                        puu(ji,jj,jk,Kbb) = zuf                     ! puu(:,:,:,Kbb) <-- filtered velocity
+                        pvv(ji,jj,jk,Kbb) = zvf
+                        puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Krhs)            ! puu(:,:,:,Kmm) <-- puu(:,:,:,Krhs)
+                        pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Krhs)
                      END DO
                   END DO
                END DO
-               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
-               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
+               e3u(:,:,1:jpkm1,Kbb) = ze3u_f(:,:,1:jpkm1)        ! e3u(:,:,:,Kbb) <-- filtered scale factor
+               e3v(:,:,1:jpkm1,Kbb) = ze3v_f(:,:,1:jpkm1)
                !
                DEALLOCATE( ze3u_f , ze3v_f )
@@ -318,13 +319,13 @@
             ! Revert "before" velocities to time split estimate
             ! Doing it here also means that asselin filter contribution is removed  
-            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
-            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)    
+            zue(:,:) = e3u(:,:,1,Kbb) * puu(:,:,1,Kbb) * umask(:,:,1)
+            zve(:,:) = e3v(:,:,1,Kbb) * pvv(:,:,1,Kbb) * vmask(:,:,1)    
             DO jk = 2, jpkm1
-               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
-               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)    
+               zue(:,:) = zue(:,:) + e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb) * umask(:,:,jk)
+               zve(:,:) = zve(:,:) + e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb) * vmask(:,:,jk)    
             END DO
             DO jk = 1, jpkm1
-               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
-               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
+               puu(:,:,jk,Kbb) = puu(:,:,jk,Kbb) - (zue(:,:) * r1_hu_n(:,:) - uu_b(:,:,Kmm)) * umask(:,:,jk)
+               pvv(:,:,jk,Kbb) = pvv(:,:,jk,Kbb) - (zve(:,:) * r1_hv_n(:,:) - vv_b(:,:,Kmm)) * vmask(:,:,jk)
             END DO
          ENDIF
@@ -338,9 +339,9 @@
       !
       IF(.NOT.ln_linssh ) THEN
-         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
-         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
+         hu_b(:,:) = e3u(:,:,1,Kbb) * umask(:,:,1)
+         hv_b(:,:) = e3v(:,:,1,Kbb) * vmask(:,:,1)
          DO jk = 2, jpkm1
-            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
-            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
+            hu_b(:,:) = hu_b(:,:) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
+            hv_b(:,:) = hv_b(:,:) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
          END DO
          r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
@@ -348,31 +349,31 @@
       ENDIF
       !
-      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
-      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
-      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
-      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
+      uu_b(:,:,Kmm) = e3u(:,:,1,Krhs) * puu(:,:,1,Kmm) * umask(:,:,1)
+      uu_b(:,:,Kbb) = e3u(:,:,1,Kbb) * puu(:,:,1,Kbb) * umask(:,:,1)
+      vv_b(:,:,Kmm) = e3v(:,:,1,Krhs) * pvv(:,:,1,Kmm) * vmask(:,:,1)
+      vv_b(:,:,Kbb) = e3v(:,:,1,Kbb) * pvv(:,:,1,Kbb) * vmask(:,:,1)
       DO jk = 2, jpkm1
-         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
-         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
-         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
-         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
+         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + e3u(:,:,jk,Krhs) * puu(:,:,jk,Kmm) * umask(:,:,jk)
+         uu_b(:,:,Kbb) = uu_b(:,:,Kbb) + e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb) * umask(:,:,jk)
+         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + e3v(:,:,jk,Krhs) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
+         vv_b(:,:,Kbb) = vv_b(:,:,Kbb) + e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb) * vmask(:,:,jk)
       END DO
-      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
-      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
-      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
-      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
+      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu_a(:,:)
+      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv_a(:,:)
+      uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu_b(:,:)
+      vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv_b(:,:)
       !
       IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
-         CALL iom_put(  "ubar", un_b(:,:) )
-         CALL iom_put(  "vbar", vn_b(:,:) )
+         CALL iom_put(  "ubar", uu_b(:,:,Kmm) )
+         CALL iom_put(  "vbar", vv_b(:,:,Kmm) )
       ENDIF
       IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
-         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
-         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
+         zua(:,:,:) = ( puu(:,:,:,Kbb) - zua(:,:,:) ) * z1_2dt
+         zva(:,:,:) = ( pvv(:,:,:,Kbb) - zva(:,:,:) ) * z1_2dt
          CALL trd_dyn( zua, zva, jpdyn_atf, kt, Kmm )
       ENDIF
       !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
-         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kmm), clinfo1=' nxt  - Un: ', mask1=umask,   &
+         &                       tab3d_2=pvv(:,:,:,Kmm), clinfo2=' Vn: '       , mask2=vmask )
       ! 
       IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg_ts.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg_ts.F90	(revision 11027)
@@ -1307,6 +1307,6 @@
       !
       IF( ln_diatmb ) THEN
-         CALL iom_put( "baro_u" , un_b*ssumask(:,:)+zmdi*(1.-ssumask(:,:) ) )  ! Barotropic  U Velocity
-         CALL iom_put( "baro_v" , vn_b*ssvmask(:,:)+zmdi*(1.-ssvmask(:,:) ) )  ! Barotropic  V Velocity
+         CALL iom_put( "baro_u" , uu_b(:,:,Kmm)*ssumask(:,:)+zmdi*(1.-ssumask(:,:) ) )  ! Barotropic  U Velocity
+         CALL iom_put( "baro_v" , vv_b(:,:,Kmm)*ssvmask(:,:)+zmdi*(1.-ssvmask(:,:) ) )  ! Barotropic  V Velocity
       ENDIF
       !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90	(revision 11027)
@@ -93,5 +93,5 @@
       !                                           !------------------------------!
       IF(ln_wd_il) THEN
-         CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), z2dt)
+         CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), z2dt, Kmm, uu, vv )
       ENDIF
 
@@ -109,5 +109,5 @@
       !
 #if defined key_agrif
-      CALL agrif_ssh( kt )
+      Krhs_a = Kaa ; CALL agrif_ssh( kt )
 #endif
       !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/wet_dry.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/wet_dry.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/wet_dry.F90	(revision 11027)
@@ -122,5 +122,5 @@
 
 
-   SUBROUTINE wad_lmt( sshb1, sshemp, z2dt )
+   SUBROUTINE wad_lmt( psshb1, psshemp, z2dt, Kmm, puu, pvv )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE wad_lmt  ***
@@ -132,7 +132,9 @@
       !! ** Action  : - calculate flux limiter and W/D flag
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   sshb1        !!gm DOCTOR names: should start with p !
-      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   sshemp
-      REAL(wp)                , INTENT(in   ) ::   z2dt
+      REAL(wp), DIMENSION(:,:)            , INTENT(inout) ::   psshb1
+      REAL(wp), DIMENSION(:,:)            , INTENT(in   ) ::   psshemp
+      REAL(wp)                            , INTENT(in   ) ::   z2dt
+      INTEGER                             , INTENT(in   ) ::   Kmm       ! time level index
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv  ! velocity arrays
       !
       INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
@@ -150,6 +152,6 @@
       !
       DO jk = 1, jpkm1
-         un(:,:,jk) = un(:,:,jk) * zwdlmtu(:,:) 
-         vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:,:) 
+         puu(:,:,jk,Kmm) = puu(:,:,jk,Kmm) * zwdlmtu(:,:) 
+         pvv(:,:,jk,Kmm) = pvv(:,:,jk,Kmm) * zwdlmtv(:,:) 
       END DO
       jflag  = 0
@@ -165,6 +167,6 @@
       !
       DO jk = 1, jpkm1     ! Horizontal Flux in u and v direction
-         zflxu(:,:) = zflxu(:,:) + e3u_n(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
-         zflxv(:,:) = zflxv(:,:) + e3v_n(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
+         zflxu(:,:) = zflxu(:,:) + e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
+         zflxv(:,:) = zflxv(:,:) + e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
       END DO
       zflxu(:,:) = zflxu(:,:) * e2u(:,:)
@@ -183,7 +185,7 @@
                &         + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji,  jj-1) , 0._wp ) 
             !
-            zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1
+            zdep2 = ht_0(ji,jj) + psshb1(ji,jj) - rn_wdmin1
             IF( zdep2 <= 0._wp ) THEN     ! add more safty, but not necessary
-               sshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj)
+               psshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj)
                IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp
                IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp
@@ -196,5 +198,5 @@
       !
       !           ! HPG limiter from jholt
-      wdramp(:,:) = min((ht_0(:,:) + sshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp)
+      wdramp(:,:) = min((ht_0(:,:) + psshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp)
       !jth assume don't need a lbc_lnk here
       DO jj = 1, jpjm1
@@ -226,5 +228,5 @@
                !
                zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
-               zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)
+               zdep2 = ht_0(ji,jj) + psshb1(ji,jj) - rn_wdmin1 - z2dt * psshemp(ji,jj)
                !
                IF( zdep1 > zdep2 ) THEN
@@ -255,18 +257,18 @@
       !
       DO jk = 1, jpkm1
-         un(:,:,jk) = un(:,:,jk) * zwdlmtu(:,:) 
-         vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:,:) 
-      END DO
-      un_b(:,:) = un_b(:,:) * zwdlmtu(:, :)
-      vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :)
+         puu(:,:,jk,Kmm) = puu(:,:,jk,Kmm) * zwdlmtu(:,:) 
+         pvv(:,:,jk,Kmm) = pvv(:,:,jk,Kmm) * zwdlmtv(:,:) 
+      END DO
+      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * zwdlmtu(:, :)
+      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * zwdlmtv(:, :)
       !
 !!gm TO BE SUPPRESSED ?  these lbc_lnk are useless since zwdlmtu and zwdlmtv are defined everywhere !
-      CALL lbc_lnk_multi( 'wet_dry', un  , 'U', -1., vn  , 'V', -1. )
-      CALL lbc_lnk_multi( 'wet_dry', un_b, 'U', -1., vn_b, 'V', -1. )
+      CALL lbc_lnk_multi( 'wet_dry', puu(:,:,:,Kmm)  , 'U', -1., pvv(:,:,:,Kmm)  , 'V', -1. )
+      CALL lbc_lnk_multi( 'wet_dry', uu_b(:,:,Kmm), 'U', -1., vv_b(:,:,Kmm), 'V', -1. )
 !!gm
       !
       IF(jflag == 1 .AND. lwp)   WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
       !
-      !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
+      !IF( ln_rnf      )   CALL sbc_rnf_div( hdiv )          ! runoffs (update hdiv field)
       !
       IF( ln_timing )   CALL timing_stop('wad_lmt')      !
@@ -392,5 +394,5 @@
       IF( jflag == 1 .AND. lwp )   WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
       !
-      !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
+      !IF( ln_rnf      )   CALL sbc_rnf_div( hdiv )          ! runoffs (update hdiv field)
       !
       IF( ln_timing )   CALL timing_stop('wad_lmt_bt')      !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/IOM/restart.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/IOM/restart.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/IOM/restart.F90	(revision 11027)
@@ -148,15 +148,15 @@
 
       IF ( .NOT. ln_diurnal_only ) THEN
-                     CALL iom_rstput( kt, nitrst, numrow, 'ub'     , uu(:,:,:,Kbb), ldxios = lwxios        )     ! before fields
-                     CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vv(:,:,:,Kbb), ldxios = lwxios        )
+                     CALL iom_rstput( kt, nitrst, numrow, 'ub'     , uu(:,:,:       ,Kbb), ldxios = lwxios        )     ! before fields
+                     CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vv(:,:,:       ,Kbb), ldxios = lwxios        )
                      CALL iom_rstput( kt, nitrst, numrow, 'tb'     , ts(:,:,:,jp_tem,Kbb), ldxios = lwxios )
                      CALL iom_rstput( kt, nitrst, numrow, 'sb'     , ts(:,:,:,jp_sal,Kbb), ldxios = lwxios )
-                     CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb, ldxios = lwxios      )
+                     CALL iom_rstput( kt, nitrst, numrow, 'sshb'   ,ssh(:,:         ,Kbb), ldxios = lwxios      )
                      !
-                     CALL iom_rstput( kt, nitrst, numrow, 'un'     , uu(:,:,:,Kmm), ldxios = lwxios        )     ! now fields
-                     CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vv(:,:,:,Kmm), ldxios = lwxios        )
+                     CALL iom_rstput( kt, nitrst, numrow, 'un'     , uu(:,:,:       ,Kmm), ldxios = lwxios        )     ! now fields
+                     CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vv(:,:,:       ,Kmm), ldxios = lwxios        )
                      CALL iom_rstput( kt, nitrst, numrow, 'tn'     , ts(:,:,:,jp_tem,Kmm), ldxios = lwxios )
                      CALL iom_rstput( kt, nitrst, numrow, 'sn'     , ts(:,:,:,jp_sal,Kmm), ldxios = lwxios )
-                     CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , sshn, ldxios = lwxios      )
+                     CALL iom_rstput( kt, nitrst, numrow, 'sshn'   ,ssh(:,:         ,Kmm), ldxios = lwxios      )
                      CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop, ldxios = lwxios      )
                   ! extra variable needed for the ice sheet coupling
@@ -275,18 +275,18 @@
       
       IF( iom_varid( numror, 'ub', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'ub'     , uu(:,:,:,Kbb), ldxios = lrxios        )   ! before fields
-         CALL iom_get( numror, jpdom_autoglo, 'vb'     , vv(:,:,:,Kbb), ldxios = lrxios        )
+         CALL iom_get( numror, jpdom_autoglo, 'ub'     , uu(:,:,:       ,Kbb), ldxios = lrxios )   ! before fields
+         CALL iom_get( numror, jpdom_autoglo, 'vb'     , vv(:,:,:       ,Kbb), ldxios = lrxios )
          CALL iom_get( numror, jpdom_autoglo, 'tb'     , ts(:,:,:,jp_tem,Kbb), ldxios = lrxios )
          CALL iom_get( numror, jpdom_autoglo, 'sb'     , ts(:,:,:,jp_sal,Kbb), ldxios = lrxios )
-         CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb, ldxios = lrxios                 )
+         CALL iom_get( numror, jpdom_autoglo, 'sshb'   ,ssh(:,:         ,Kbb), ldxios = lrxios )
       ELSE
          neuler = 0
       ENDIF
       !
-      CALL iom_get( numror, jpdom_autoglo, 'un'     , uu(:,:,:,Kmm), ldxios = lrxios        )   ! now    fields
-      CALL iom_get( numror, jpdom_autoglo, 'vn'     , vv(:,:,:,Kmm), ldxios = lrxios        )
+      CALL iom_get( numror, jpdom_autoglo, 'un'     , uu(:,:,:       ,Kmm), ldxios = lrxios )       ! now    fields
+      CALL iom_get( numror, jpdom_autoglo, 'vn'     , vv(:,:,:       ,Kmm), ldxios = lrxios )
       CALL iom_get( numror, jpdom_autoglo, 'tn'     , ts(:,:,:,jp_tem,Kmm), ldxios = lrxios )
       CALL iom_get( numror, jpdom_autoglo, 'sn'     , ts(:,:,:,jp_sal,Kmm), ldxios = lrxios )
-      CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios )
+      CALL iom_get( numror, jpdom_autoglo, 'sshn'   ,ssh(:,:         ,Kmm), ldxios = lrxios )
       IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN
          CALL iom_get( numror, jpdom_autoglo, 'rhop'   , rhop, ldxios = lrxios )   ! now    potential density
@@ -297,7 +297,7 @@
       IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0)
          ts   (:,:,:,:,Kbb) = ts   (:,:,:,:,Kmm)              ! all before fields set to now values
-         uu   (:,:,:,Kbb)   = uu   (:,:,:,Kmm)
-         vv   (:,:,:,Kbb)   = vv   (:,:,:,Kmm)
-         sshb (:,:)         = sshn (:,:)
+         uu   (:,:,:  ,Kbb) = uu   (:,:,:  ,Kmm)
+         vv   (:,:,:  ,Kbb) = vv   (:,:,:  ,Kmm)
+         ssh  (:,:    ,Kbb) = ssh  (:,:    ,Kmm)
          !
          IF( .NOT.ln_linssh ) THEN
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/diaobs.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/diaobs.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/diaobs.F90	(revision 11027)
@@ -460,5 +460,5 @@
             !
             IF( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
-               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) )
+               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype), Kmm )
                IF( ln_altbias )   &
                   & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile )
@@ -499,5 +499,5 @@
       USE dom_oce, ONLY : gdept, gdept_1d     ! Ocean space domain variables (Kmm time-level only)
       USE phycst , ONLY : rday                ! Physical constants
-      USE oce    , ONLY : ts, uu, vv, sshn    ! Ocean dynamics and tracers variables (Kmm time-level only)
+      USE oce    , ONLY : ts, uu, vv, ssh     ! Ocean dynamics and tracers variables (Kmm time-level only)
       USE phycst , ONLY : rday                ! Physical constants
 #if defined  key_si3
@@ -598,5 +598,5 @@
                zsurfvar(:,:) = ts(:,:,1,jp_tem,Kmm)
             CASE('sla')
-               zsurfvar(:,:) = sshn(:,:)
+               zsurfvar(:,:) = ssh(:,:,Kmm)
             CASE('sss')
                zsurfvar(:,:) = ts(:,:,1,jp_sal,Kmm)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/obs_read_altbias.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/obs_read_altbias.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/obs_read_altbias.F90	(revision 11027)
@@ -29,5 +29,5 @@
       & gphit
    USE oce, ONLY : &           ! Model variables
-      & sshn
+      & ssh
    USE obs_inter_h2d
    USE obs_utils               ! Various observation tools
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/obs_readmdt.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/obs_readmdt.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/obs_readmdt.F90	(revision 11027)
@@ -25,5 +25,5 @@
       &                    tmask, tmask_i, e1e2t, gphit, glamt
    USE obs_const, ONLY :   obfillflt      ! Fillvalue
-   USE oce      , ONLY :   sshn           ! Model variables
+   USE oce      , ONLY :   ssh            ! Model variables
 
    IMPLICIT NONE
@@ -44,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE obs_rea_mdt( sladata, k2dint )
+   SUBROUTINE obs_rea_mdt( sladata, k2dint, Kmm )
       !!---------------------------------------------------------------------
       !!
@@ -59,4 +59,5 @@
       TYPE(obs_surf), INTENT(inout) ::   sladata   ! SLA data
       INTEGER       , INTENT(in)    ::   k2dint    ! ?
+      INTEGER       , INTENT(in)    ::   Kmm       ! ?
       !
       CHARACTER(LEN=12), PARAMETER ::   cpname  = 'obs_rea_mdt'
@@ -106,5 +107,5 @@
       ! Remove the offset between the MDT used with the sla and the model MDT
       IF( nn_msshc == 1 .OR. nn_msshc == 2 ) &
-         & CALL obs_offset_mdt( jpi, jpj, z_mdt, zfill )
+         & CALL obs_offset_mdt( jpi, jpj, z_mdt, zfill, Kmm )
 
       ! Intepolate the MDT already on the model grid at the observation point
@@ -169,5 +170,5 @@
 
 
-   SUBROUTINE obs_offset_mdt( kpi, kpj, mdt, zfill )
+   SUBROUTINE obs_offset_mdt( kpi, kpj, mdt, zfill, Kmm )
       !!---------------------------------------------------------------------
       !!
@@ -183,4 +184,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(IN) ::  kpi, kpj
+      INTEGER, INTENT(IN) ::  Kmm
       REAL(wp), DIMENSION(kpi,kpj), INTENT(INOUT) ::   mdt     ! MDT used on the model grid
       REAL(wp)                    , INTENT(IN   ) ::   zfill 
@@ -216,5 +218,5 @@
           zarea = zarea + zdxdy
           zeta1 = zeta1 + mdt(ji,jj) * zdxdy
-          zeta2 = zeta2 + sshn (ji,jj) * zdxdy
+          zeta2 = zeta2 + ssh(ji,jj,Kmm) * zdxdy
         END DO      
       END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/obs_sstbias.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/obs_sstbias.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/OBS/obs_sstbias.F90	(revision 11027)
@@ -28,5 +28,5 @@
       & glamt
    USE oce, ONLY : &           ! Model variables
-      & sshn
+      & ssh
    USE obs_inter_h2d
    USE obs_utils               ! Various observation tools
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbccpl.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbccpl.F90	(revision 11027)
@@ -32,5 +32,5 @@
    USE cpl_oasis3     ! OASIS3 coupling
    USE geo2ocean      ! 
-   USE oce     , ONLY : ts, uu, vv, sshn, sshb, fraqsr_1lev
+   USE oce     , ONLY : ts, uu, vv, ssh, fraqsr_1lev
    USE ocealb         ! 
    USE eosbn2         ! 
@@ -2037,5 +2037,5 @@
    
    
-   SUBROUTINE sbc_cpl_snd( kt, Kmm )
+   SUBROUTINE sbc_cpl_snd( kt, Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!             ***  ROUTINE sbc_cpl_snd  ***
@@ -2047,5 +2047,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt
-      INTEGER, INTENT(in) ::   Kmm    ! ocean model time level index
+      INTEGER, INTENT(in) ::   Kbb, Kmm    ! ocean model time level index
       !
       INTEGER ::   ji, jj, jl   ! dummy loop indices
@@ -2476,10 +2476,10 @@
          IF( ln_apr_dyn ) THEN  
             IF( kt /= nit000 ) THEN  
-               ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )  
+               ztmp1(:,:) = ssh(:,:,Kbb) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )  
             ELSE  
-               ztmp1(:,:) = sshb(:,:)  
+               ztmp1(:,:) = ssh(:,:,Kbb)  
             ENDIF  
          ELSE  
-            ztmp1(:,:) = sshn(:,:)  
+            ztmp1(:,:) = ssh(:,:,Kmm)  
          ENDIF  
          CALL cpl_snd( jps_wlev  , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
@@ -2491,6 +2491,6 @@
          !                          ! removed inverse barometer ssh when Patm
          !                          forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
-         ELSE                    ;   ztmp1(:,:) = sshn(:,:)
+         IF( ln_apr_dyn ) THEN   ;   ztmp1(:,:) = ssh(:,:,Kbb) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
+         ELSE                    ;   ztmp1(:,:) = ssh(:,:,Kmm)
          ENDIF
          CALL cpl_snd( jps_ssh   , isec, RESHAPE ( ztmp1            , (/jpi,jpj,1/) ), info )
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcfwb.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcfwb.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcfwb.F90	(revision 11027)
@@ -48,5 +48,5 @@
 CONTAINS
 
-   SUBROUTINE sbc_fwb( kt, kn_fwb, kn_fsbc )
+   SUBROUTINE sbc_fwb( kt, kn_fwb, kn_fsbc, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_fwb  ***
@@ -65,4 +65,5 @@
       INTEGER, INTENT( in ) ::   kn_fsbc  ! 
       INTEGER, INTENT( in ) ::   kn_fwb   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm      ! ocean time level index
       !
       INTEGER  ::   inum, ikty, iyear     ! local integers
@@ -131,5 +132,5 @@
             a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
                                                       ! sum over the global domain
-            a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) )
+            a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rau0 ) )
             a_fwb   = a_fwb * 1.e+3 / ( area * rday * 365. )     ! convert in Kg/m3/s = mm/s
 !!gm        !                                                      !!bug 365d year 
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcice_cice.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcice_cice.F90	(revision 11027)
@@ -147,5 +147,5 @@
 
 
-   SUBROUTINE cice_sbc_init( ksbc )
+   SUBROUTINE cice_sbc_init( ksbc, Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE cice_sbc_init  ***
@@ -154,4 +154,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT( in  ) ::   ksbc                ! surface forcing type
+      INTEGER, INTENT( in  ) ::   Kbb, Kmm            ! time level indices
       REAL(wp), DIMENSION(jpi,jpj) :: ztmp1, ztmp2
       REAL(wp) ::   zcoefu, zcoefv, zcoeff            ! local scalar
@@ -227,6 +228,6 @@
       IF( .NOT.ln_rstart ) THEN
          IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
-            sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
-            sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
+            ssh(:,:,Kmm) = ssh(:,:,Kmm) - snwice_mass(:,:) * r1_rau0
+            ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rau0
 
 !!gm This should be put elsewhere....   (same remark for limsbc)
@@ -235,6 +236,6 @@
                !
                DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
-                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-                  e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+                  e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
                ENDDO
                e3t(:,:,:,Krhs) = e3t(:,:,:,Kbb)
@@ -259,5 +260,5 @@
                gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
                gdepw(:,:,1,Kmm) = 0.0_wp
-               gde3w(:,:,1)     = gdept(:,:,1,Kmm) - sshn(:,:)
+               gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
                DO jk = 2, jpk
                   gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk,Kmm)
@@ -1056,7 +1057,8 @@
    END SUBROUTINE sbc_ice_cice
 
-   SUBROUTINE cice_sbc_init (ksbc)    ! Dummy routine
+   SUBROUTINE cice_sbc_init (ksbc, Kbb, Kmm)    ! Dummy routine
       IMPLICIT NONE
       INTEGER, INTENT( in ) :: ksbc
+      INTEGER, INTENT( in ) :: Kbb, Kmm
       WRITE(*,*) 'cice_sbc_init: You should not have seen this print! error?', ksbc
    END SUBROUTINE cice_sbc_init
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcisf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcisf.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcisf.F90	(revision 11027)
@@ -149,5 +149,5 @@
          DO jj = 1,jpj
             DO ji = 1,jpi
-               zdep(ji,jj)=gdepw_n(ji,jj,misfkt(ji,jj))
+               zdep(ji,jj)=gdepw(ji,jj,misfkt(ji,jj),Kmm)
             END DO
          END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcmod.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcmod.F90	(revision 11027)
@@ -341,10 +341,10 @@
                           IF( sbc_ice_alloc() /= 0 )   CALL ctl_stop('STOP', 'sbc_ice_alloc : unable to allocate arrays' )
       ELSEIF( nn_ice == 2 ) THEN
-                          CALL ice_init( Kbb, Kmm, Kaa )  ! ICE initialization
+                          CALL ice_init( Kbb, Kmm, Kaa )         ! ICE initialization
       ENDIF
 #endif
-      IF( nn_ice == 3 )   CALL cice_sbc_init( nsbc )   ! CICE initialization
-      !
-      IF( ln_wave     )   CALL sbc_wave_init           ! surface wave initialisation
+      IF( nn_ice == 3 )   CALL cice_sbc_init( nsbc, Kbb, Kmm )   ! CICE initialization
+      !
+      IF( ln_wave     )   CALL sbc_wave_init                     ! surface wave initialisation
       !
       IF( lwxios ) THEN
@@ -454,11 +454,11 @@
       ENDIF
 
-      IF( ln_isf         )   CALL sbc_isf( kt, Kmm )              ! compute iceshelves
-
-      IF( ln_rnf         )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes
-
-      IF( ln_ssr         )   CALL sbc_ssr( kt )                   ! add SST/SSS damping term
-
-      IF( nn_fwb    /= 0 )   CALL sbc_fwb( kt, nn_fwb, nn_fsbc )  ! control the freshwater budget
+      IF( ln_isf         )   CALL sbc_isf( kt, Kmm )                   ! compute iceshelves
+
+      IF( ln_rnf         )   CALL sbc_rnf( kt )                        ! add runoffs to fresh water fluxes
+
+      IF( ln_ssr         )   CALL sbc_ssr( kt )                        ! add SST/SSS damping term
+
+      IF( nn_fwb    /= 0 )   CALL sbc_fwb( kt, nn_fwb, nn_fsbc, Kmm )  ! control the freshwater budget
 
       ! Special treatment of freshwater fluxes over closed seas in the model domain
@@ -471,5 +471,5 @@
       IF ( ll_wd ) THEN     ! If near WAD point limit the flux for now
          zthscl = atanh(rn_wd_sbcfra)                     ! taper frac default is .999 
-         zwdht(:,:) = sshn(:,:) + ht_0(:,:) - rn_wdmin1   ! do this calc of water
+         zwdht(:,:) = ssh(:,:,Kmm) + ht_0(:,:) - rn_wdmin1   ! do this calc of water
                                                      ! depth above wd limit once
          WHERE( zwdht(:,:) <= 0.0 )
@@ -557,14 +557,14 @@
       !
       IF(ln_ctl) THEN         ! print mean trends (used for debugging)
-         CALL prt_ctl(tab2d_1=fr_i              , clinfo1=' fr_i    - : ', mask1=tmask )
-         CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf), clinfo1=' emp-rnf - : ', mask1=tmask )
-         CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf), clinfo1=' sfx-rnf - : ', mask1=tmask )
-         CALL prt_ctl(tab2d_1=qns              , clinfo1=' qns      - : ', mask1=tmask )
-         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr      - : ', mask1=tmask )
-         CALL prt_ctl(tab3d_1=tmask            , clinfo1=' tmask    - : ', mask1=tmask, kdim=jpk )
-         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' sst      - : ', mask1=tmask, kdim=1   )
-         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sss      - : ', mask1=tmask, kdim=1   )
-         CALL prt_ctl(tab2d_1=utau             , clinfo1=' utau     - : ', mask1=umask,                      &
-            &         tab2d_2=vtau             , clinfo2=' vtau     - : ', mask2=vmask )
+         CALL prt_ctl(tab2d_1=fr_i                , clinfo1=' fr_i    - : ' , mask1=tmask )
+         CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf)  , clinfo1=' emp-rnf - : ' , mask1=tmask )
+         CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf)  , clinfo1=' sfx-rnf - : ' , mask1=tmask )
+         CALL prt_ctl(tab2d_1=qns                 , clinfo1=' qns      - : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=qsr                 , clinfo1=' qsr      - : ', mask1=tmask )
+         CALL prt_ctl(tab3d_1=tmask               , clinfo1=' tmask    - : ', mask1=tmask, kdim=jpk )
+         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' sst      - : ', mask1=tmask, kdim=1   )
+         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sss      - : ', mask1=tmask, kdim=1   )
+         CALL prt_ctl(tab2d_1=utau                , clinfo1=' utau     - : ', mask1=umask,                      &
+            &         tab2d_2=vtau                , clinfo2=' vtau     - : ', mask2=vmask )
       ENDIF
 
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcssm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcssm.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcssm.F90	(revision 11027)
@@ -77,6 +77,6 @@
          sss_m(:,:) = zts(:,:,jp_sal)
          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
-         ELSE                    ;   ssh_m(:,:) = sshn(:,:)
+         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
+         ELSE                    ;   ssh_m(:,:) = ssh(:,:,Kmm)
          ENDIF
          !
@@ -100,6 +100,6 @@
             sss_m(:,:) = zcoef * zts(:,:,jp_sal)
             !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) )
-            ELSE                    ;   ssh_m(:,:) = zcoef * sshn(:,:)
+            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) )
+            ELSE                    ;   ssh_m(:,:) = zcoef *   ssh(:,:,Kmm)
             ENDIF
             !
@@ -128,6 +128,6 @@
          sss_m(:,:) = sss_m(:,:) + zts(:,:,jp_sal)
          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
-         ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:)
+         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
+         ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + ssh(:,:,Kmm)
          ENDIF
          !
@@ -250,5 +250,5 @@
          ENDIF
          sss_m(:,:) = ts  (:,:,1,jp_sal,Kmm)
-         ssh_m(:,:) = sshn(:,:)
+         ssh_m(:,:) = ssh(:,:,Kmm)
          e3t_m(:,:) = e3t (:,:,1,Kmm)
          frq_m(:,:) = 1._wp
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90	(revision 11027)
@@ -137,4 +137,5 @@
       !                            !-----------------------!
 #if defined key_agrif
+      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
       CALL Agrif_Declare_Var_dom   ! AGRIF: set the meshes for DOM
       CALL Agrif_Declare_Var       !  "      "   "   "      "  DYN/TRA 
@@ -169,4 +170,5 @@
       !
       ! Recursive update from highest nested level to lowest:
+      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nnn    ! agrif_oce module copies of time level indices
       CALL Agrif_step_child_adj(Agrif_Update_All)
       !
@@ -404,6 +406,4 @@
       Nbb = 1; Nnn = 2; Naa = 3; Nrhs = Naa
 
-      ! Initialisation of temporary pointers (to be deleted after development finished)
-      CALL update_pointers( Nbb, Nnn, Naa )
       !                             !-------------------------------!
       !                             !  NEMO general initialization  !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/oce.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/oce.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/oce.F90	(revision 11027)
@@ -66,15 +66,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fraqsr_1lev        !: fraction of solar net radiation absorbed in the first ocean level [-]
 
-   !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::   ub   ,  un    , ua       !: i-horizontal velocity        [m/s]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::   vb   ,  vn    , va       !: j-horizontal velocity        [m/s]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::           wn               !: k-vertical   velocity        [m/s]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:)   ::           hdivn            !: horizontal divergence        [s-1]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:,:) ::   tsb  ,  tsn   , tsa      !: 4D T-S fields                [Celsius,psu]         
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)     ::   ub_b   ,  un_b  ,  ua_b  !: Barotropic velocities at u-point [m/s]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)     ::   vb_b   ,  vn_b  ,  va_b  !: Barotropic velocities at v-point [m/s]
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:)     ::   sshb   ,  sshn  ,  ssha  !: sea surface height at t-point [m]
-   !! TEMPORARY POINTERS FOR DEVELOPMENT ONLY
-
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 11027)
@@ -44,5 +44,4 @@
 
    PUBLIC   stp   ! called by nemogcm.F90
-   PUBLIC   update_pointers ! called by nemo_init
 
    !!----------------------------------------------------------------------
@@ -253,11 +252,11 @@
 #endif
                          CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection	==> RHS
-      IF( ln_zdfosm  )   CALL tra_osm    ( kstp, Nnn     , ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
+      IF( ln_zdfosm  )   CALL tra_osm    ( kstp,      Nnn, ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
       IF( lrst_oce .AND. ln_zdfosm ) &
-           &             CALL osm_rst    ( kstp, Nnn, 'WRITE' )        ! write OSMOSIS outputs + ww (so must do here) to restarts
+           &             CALL osm_rst    ( kstp,      Nnn, 'WRITE'  )  ! write OSMOSIS outputs + ww (so must do here) to restarts
                          CALL tra_ldf    ( kstp, Nbb, Nnn, ts, Nrhs )  ! lateral mixing
 
 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
-      IF( ln_diaptr  )   CALL dia_ptr( Nnn )                 ! Poleward adv/ldf TRansports diagnostics
+      IF( ln_diaptr  )   CALL dia_ptr( Nnn )                           ! Poleward adv/ldf TRansports diagnostics
 !!gm
                          CALL tra_zdf    ( kstp, Nbb, Nnn, Nrhs, ts, Naa  )  ! vert. mixing & after tracer	==> after
@@ -282,5 +281,5 @@
 !!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine
                          CALL tra_nxt       ( kstp, Nbb, Nnn, Nrhs, Naa )  ! finalize (bcs) tracer fields at next time step and swap
-                         CALL dyn_nxt       ( kstp, Nbb, Nnn, Naa  )  ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
+                         CALL dyn_nxt       ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
                          CALL ssh_swp       ( kstp, Nbb, Nnn, Naa )  ! swap of sea surface height
       IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa )  ! swap of vertical scale factors
@@ -300,7 +299,10 @@
       ! AGRIF
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<      
-                         CALL Agrif_Integrate_ChildGrids( stp )  ! allows to finish all the Child Grids before updating
-
-                         IF( Agrif_NbStepint() == 0 ) CALL Agrif_update_all( ) ! Update all components
+                         CALL Agrif_Integrate_ChildGrids( stp )       ! allows to finish all the Child Grids before updating
+
+                         IF( Agrif_NbStepint() == 0 ) THEN
+                            Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+                            CALL Agrif_update_all( )                  ! Update all components
+                         ENDIF
 #endif
       IF( ln_diaobs  )   CALL dia_obs      ( kstp, Nnn )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
@@ -321,5 +323,5 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 !!gm why lk_oasis and not lk_cpl ????
-      IF( lk_oasis   )   CALL sbc_cpl_snd( kstp, Nnn )     ! coupled mode : field exchanges
+      IF( lk_oasis   )   CALL sbc_cpl_snd( kstp, Nbb, Nnn )        ! coupled mode : field exchanges
       !
 #if defined key_iomput
@@ -335,42 +337,4 @@
    END SUBROUTINE stp
    
-   SUBROUTINE update_pointers( Kbb, Kmm, Kaa )
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE update_pointers  ***
-      !!
-      !! ** Purpose :   Associate temporary pointer arrays.
-      !!                For IMMERSE development phase only - to be deleted
-      !!
-      !! ** Method  :
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices
-
-      ub => uu(:,:,:,Kbb); un => uu(:,:,:,Kmm); ua => uu(:,:,:,Kaa)
-      vb => vv(:,:,:,Kbb); vn => vv(:,:,:,Kmm); va => vv(:,:,:,Kaa)
-      wn => ww(:,:,:)
-      hdivn => hdiv(:,:,:)
-
-      sshb =>  ssh(:,:,Kbb); sshn =>  ssh(:,:,Kmm); ssha =>  ssh(:,:,Kaa)
-      ub_b => uu_b(:,:,Kbb); un_b => uu_b(:,:,Kmm); ua_b => uu_b(:,:,Kaa)
-      vb_b => vv_b(:,:,Kbb); vn_b => vv_b(:,:,Kmm); va_b => vv_b(:,:,Kaa)
-
-      tsb => ts(:,:,:,:,Kbb); tsn => ts(:,:,:,:,Kmm); tsa => ts(:,:,:,:,Kaa)
-
-      e3t_b => e3t(:,:,:,Kbb); e3t_n => e3t(:,:,:,Kmm); e3t_a => e3t(:,:,:,Kaa)
-      e3u_b => e3u(:,:,:,Kbb); e3u_n => e3u(:,:,:,Kmm); e3u_a => e3u(:,:,:,Kaa)
-      e3v_b => e3v(:,:,:,Kbb); e3v_n => e3v(:,:,:,Kmm); e3v_a => e3v(:,:,:,Kaa)
-
-      e3f_n => e3f(:,:,:)
-
-      e3w_b  => e3w (:,:,:,Kbb); e3w_n  => e3w (:,:,:,Kmm)
-      e3uw_b => e3uw(:,:,:,Kbb); e3uw_n => e3uw(:,:,:,Kmm)
-      e3vw_b => e3vw(:,:,:,Kbb); e3vw_n => e3vw(:,:,:,Kmm)
-
-      gdept_b => gdept(:,:,:,Kbb); gdept_n => gdept(:,:,:,Kmm) 
-      gdepw_b => gdepw(:,:,:,Kbb); gdepw_n => gdepw(:,:,:,Kmm) 
-      gde3w_n => gde3w(:,:,:)
-
-   END SUBROUTINE update_pointers
-
    !!======================================================================
 END MODULE step
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90	(revision 11027)
@@ -182,5 +182,5 @@
          CALL prt_ctl(tab3d_1=uu(:,:,:,Kmm)               , clinfo1=' uu(:,:,:,Kmm)      - : ', mask1=umask,  kdim=jpk   )
          CALL prt_ctl(tab3d_1=vv(:,:,:,Kmm)               , clinfo1=' vv(:,:,:,Kmm)      - : ', mask1=vmask,  kdim=jpk   )
-         CALL prt_ctl(tab3d_1=wn               , clinfo1=' ww      - : ', mask1=tmask,  kdim=jpk   )
+         CALL prt_ctl(tab3d_1=ww               , clinfo1=' ww      - : ', mask1=tmask,  kdim=jpk   )
          CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask,  kdim=jpk   )
          CALL prt_ctl(tab3d_1=uslp             , clinfo1=' slp  - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90	(revision 11027)
@@ -59,5 +59,4 @@
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
    USE lbcnfd  , ONLY : isendto, nsndto, nfsloop, nfeloop   ! Setup of north fold exchanges
-   USE step, ONLY : update_pointers 
 
    IMPLICIT NONE
@@ -296,6 +295,4 @@
       Nbb = 1; Nnn = 2; Naa = 3; Nrhs = Naa
 	
-      ! Initialisation of temporary pointers (to be deleted after development finished)
-      CALL update_pointers( Nbb, Nnn, Naa )
 
       !                             !-------------------------------!
@@ -536,5 +533,5 @@
       vv   (:,:,:,Kmm)   = 0._wp   ;   vv(:,:,:,Kaa) = 0._wp   !
       ww   (:,:,:)   = 0._wp   !                       !
-      hdivn(:,:,:)   = 0._wp   !                       !
+      hdiv (:,:,:)   = 0._wp   !                       !
       ts  (:,:,:,:,Kmm) = 0._wp   !                       !
       !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/diawri.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/diawri.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/diawri.F90	(revision 11027)
@@ -78,5 +78,5 @@
 
    
-   SUBROUTINE dia_wri( kt )
+   SUBROUTINE dia_wri( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE dia_wri  ***
@@ -90,9 +90,10 @@
       !!
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm     ! ocean time levelindex
       !!----------------------------------------------------------------------
       ! 
       ! Output the initial state and forcings
       IF( ninist == 1 ) THEN
-         CALL dia_wri_state( 'output.init' )
+         CALL dia_wri_state( 'output.init', Kmm )
          ninist = 0
       ENDIF
@@ -330,5 +331,5 @@
 #endif
 
-   SUBROUTINE dia_wri_state( cdfile_name )
+   SUBROUTINE dia_wri_state( cdfile_name, Kmm )
       !!---------------------------------------------------------------------
       !!                 ***  ROUTINE dia_wri_state  ***
@@ -344,4 +345,5 @@
       !!----------------------------------------------------------------------
       CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
+      INTEGER           , INTENT( in ) ::   Kmm              ! ocean time levelindex
       !!
       INTEGER :: inum
@@ -359,16 +361,16 @@
 #endif
 
-      CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) )    ! now temperature
-      CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) )    ! now salinity
-      CALL iom_rstput( 0, 0, inum, 'sossheig', sshn              )    ! sea surface height
-      CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity
-      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity
-      CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! now k-velocity
-      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
-      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
-      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
-      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
-      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
-      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
+      CALL iom_rstput( 0, 0, inum, 'votemper', ts (:,:,:,jp_tem,Kmm) )    ! now temperature
+      CALL iom_rstput( 0, 0, inum, 'vosaline', ts (:,:,:,jp_sal,Kmm) )    ! now salinity
+      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,         Kmm) )    ! sea surface height
+      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu (:,:,:,       Kmm) )    ! now i-velocity
+      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv (:,:,:,       Kmm) )    ! now j-velocity
+      CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww                    )    ! now k-velocity
+      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf             )    ! freshwater budget
+      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns             )    ! total heat flux
+      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr                   )    ! solar heat flux
+      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i                  )    ! ice fraction
+      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau                  )    ! i-wind stress
+      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau                  )    ! j-wind stress
  
 #if defined key_si3
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/nemogcm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/nemogcm.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/nemogcm.F90	(revision 11027)
@@ -355,6 +355,4 @@
       Nbb = 1; Nnn = 2; Naa = 3; Nrhs = Naa
 
-      ! Initialisation of temporary pointers (to be deleted after development finished)
-      CALL update_pointers( Nbb, Nnn, Naa )
       !                             !-------------------------------!
       !                             !  NEMO general initialization  !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/sbcssm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/sbcssm.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/sbcssm.F90	(revision 11027)
@@ -121,15 +121,15 @@
          IF( .NOT. ln_linssh ) e3t_m(:,:) = e3t_0(:,:,1) !clem: necessary at least for sas2D
          frq_m(:,:) = 1._wp                              !              - -
-         sshn (:,:) = 0._wp                              !              - -
+         ssh  (:,:,Kmm) = 0._wp                              !              - -
       ENDIF
       
       IF ( nn_ice == 1 ) THEN
-         tsn(:,:,1,jp_tem) = sst_m(:,:)
-         tsn(:,:,1,jp_sal) = sss_m(:,:)
-         tsb(:,:,1,jp_tem) = sst_m(:,:)
-         tsb(:,:,1,jp_sal) = sss_m(:,:)
-      ENDIF
-      ub (:,:,1) = ssu_m(:,:)
-      vb (:,:,1) = ssv_m(:,:)
+         ts(:,:,1,jp_tem,Kmm) = sst_m(:,:)
+         ts(:,:,1,jp_sal,Kmm) = sss_m(:,:)
+         ts(:,:,1,jp_tem,Kbb) = sst_m(:,:)
+         ts(:,:,1,jp_sal,Kbb) = sss_m(:,:)
+      ENDIF
+      uu (:,:,1,Kbb) = ssu_m(:,:)
+      vv (:,:,1,Kbb) = ssv_m(:,:)
  
       IF(ln_ctl) THEN                  ! print control
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/step.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/step.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/step.F90	(revision 11027)
@@ -47,5 +47,4 @@
 
    PUBLIC   stp   ! called by nemogcm.F90
-   PUBLIC   update_pointers ! called by nemo_init
 
    !!----------------------------------------------------------------------
@@ -105,5 +104,5 @@
                              CALL sbc    ( kstp, Nbb, Nnn )                   ! Sea Boundary Condition (including sea-ice)
 
-                             CALL dia_wri( kstp )         ! ocean model: outputs
+                             CALL dia_wri( kstp,      Nnn )                   ! ocean model: outputs
 
 #if defined key_agrif
@@ -126,12 +125,12 @@
       IF( indic < 0  )  THEN
                              CALL ctl_stop( 'step: indic < 0' )
-                             CALL dia_wri_state( 'output.abort' )
+                             CALL dia_wri_state( 'output.abort', Nnn )
       ENDIF
-      IF( kstp == nit000   ) CALL iom_close( numror )     ! close input  ocean restart file
+      IF( kstp == nit000   ) CALL iom_close( numror )           ! close input  ocean restart file
       
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Coupled mode
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_oasis    )  CALL sbc_cpl_snd( kstp, Nnn )     ! coupled mode : field exchanges if OASIS-coupled ice
+      IF( lk_oasis    )  CALL sbc_cpl_snd( kstp, Nbb, Nnn )     ! coupled mode : field exchanges if OASIS-coupled ice
 
 #if defined key_iomput
@@ -153,42 +152,4 @@
    END SUBROUTINE stp
 
-   SUBROUTINE update_pointers( Kbb, Kmm, Kaa )
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE update_pointers  ***
-      !!
-      !! ** Purpose :   Associate temporary pointer arrays.
-      !!                For IMMERSE development phase only - to be deleted
-      !!
-      !! ** Method  :
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices
-
-      ub => uu(:,:,:,Kbb); un => uu(:,:,:,Kmm); ua => uu(:,:,:,Kaa)
-      vb => vv(:,:,:,Kbb); vn => vv(:,:,:,Kmm); va => vv(:,:,:,Kaa)
-      wn => ww(:,:,:)
-      hdivn => hdiv(:,:,:)
-
-      sshb =>  ssh(:,:,Kbb); sshn =>  ssh(:,:,Kmm); ssha =>  ssh(:,:,Kaa)
-      ub_b => uu_b(:,:,Kbb); un_b => uu_b(:,:,Kmm); ua_b => uu_b(:,:,Kaa)
-      vb_b => vv_b(:,:,Kbb); vn_b => vv_b(:,:,Kmm); va_b => vv_b(:,:,Kaa)
-
-      tsb => ts(:,:,:,:,Kbb); tsn => ts(:,:,:,:,Kmm); tsa => ts(:,:,:,:,Kaa)
-
-      e3t_b => e3t(:,:,:,Kbb); e3t_n => e3t(:,:,:,Kmm); e3t_a => e3t(:,:,:,Kaa)
-      e3u_b => e3u(:,:,:,Kbb); e3u_n => e3u(:,:,:,Kmm); e3u_a => e3u(:,:,:,Kaa)
-      e3v_b => e3v(:,:,:,Kbb); e3v_n => e3v(:,:,:,Kmm); e3v_a => e3v(:,:,:,Kaa)
-
-      e3f_n => e3f(:,:,:)
-
-      e3w_b  => e3w (:,:,:,Kbb); e3w_n  => e3w (:,:,:,Kmm)
-      e3uw_b => e3uw(:,:,:,Kbb); e3uw_n => e3uw(:,:,:,Kmm)
-      e3vw_b => e3vw(:,:,:,Kbb); e3vw_n => e3vw(:,:,:,Kmm)
-
-      gdept_b => gdept(:,:,:,Kbb); gdept_n => gdept(:,:,:,Kmm) 
-      gdepw_b => gdepw(:,:,:,Kbb); gdepw_n => gdepw(:,:,:,Kmm) 
-      gde3w_n => gde3w(:,:,:)
-
-   END SUBROUTINE update_pointers
-
    !!======================================================================
 END MODULE step
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/oce_trc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/oce_trc.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/oce_trc.F90	(revision 11027)
@@ -34,13 +34,4 @@
 
    !* ocean fields: here now and after fields *
-   USE oce , ONLY :   tsn     =>    tsn     !: 4D array contaning ( tn, sn ) !TEMPORARY
-   USE oce , ONLY :   tsb     =>    tsb     !: 4D array contaning ( tb, sb ) !TEMPORARY
-   USE oce , ONLY :   tsa     =>    tsa     !: 4D array contaning ( ta, sa ) !TEMPORARY
-   USE oce , ONLY :   sshn    =>    sshn    !: sea surface height at t-point [m]    !TEMPORARY
-   USE oce , ONLY :   sshb    =>    sshb    !: sea surface height at t-point [m]    !TEMPORARY
-   USE oce , ONLY :   ssha    =>    ssha    !: sea surface height at t-point [m]    !TEMPORARY
-   USE oce , ONLY :    un     =>     un     !: 4D array  !TEMPORARY
-   USE oce , ONLY :    vn     =>     vn     !: 4D array  !TEMPORARY
-   USE oce , ONLY :    wn     =>     wn     !: 4D array  !TEMPORARY
    USE oce , ONLY :   uu     =>    uu     !: i-horizontal velocity (m s-1) 
    USE oce , ONLY :   vv     =>    vv     !: j-horizontal velocity (m s-1)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trc.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trc.F90	(revision 11027)
@@ -40,10 +40,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) ::  trc_o          !: prescribed tracer concentration in ocean for SBC
    INTEGER             , PUBLIC                            ::  nn_ice_tr      !: handling of sea ice tracers
-
-   !! TEMPORARY POINTERS - TO BE DELETED AFTER IMMERSE DEVELOPMENT COMPLETE
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:,:)   ::  trn            !: tracer concentration for now time step
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:,:)   ::  tra            !: tracer concentration for next time step
-   REAL(wp), PUBLIC, POINTER, SAVE, DIMENSION(:,:,:,:)   ::  trb            !: tracer concentration for before time step
-   !! TEMPORARY POINTERS - TO BE DELETED AFTER IMMERSE DEVELOPMENT COMPLETE
 
    !! interpolated gradient
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcini.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcini.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcini.F90	(revision 11027)
@@ -52,6 +52,4 @@
       !!                or read data or analytical formulation
       !!---------------------------------------------------------------------
-      !! Time level indices only required for call to update_pointers_trc
-      !! To be removed after IMMERSE development finished. 
       INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices
       !
@@ -66,6 +64,4 @@
       CALL top_alloc()   ! allocate TOP arrays
 
-      ! Initialisation of temporary pointers (to be deleted after development finished)
-      CALL update_pointers_trc( Kbb, Kmm, Kaa )
       !
       IF(.NOT.ln_trcdta )   ln_trc_ini(:) = .FALSE.
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90	(revision 11024)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90	(revision 11027)
@@ -30,5 +30,4 @@
 
    PUBLIC   trc_stp    ! called by step
-   PUBLIC   update_pointers_trc ! called in initialisation
 
    LOGICAL  ::   llnew                   ! ???
@@ -126,19 +125,4 @@
       !
    END SUBROUTINE trc_stp
-
-   SUBROUTINE update_pointers_trc( Kbb, Kmm, Kaa )
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE update_pointers_trc  ***
-      !!
-      !! ** Purpose :   Associate temporary pointer arrays.
-      !!                For IMMERSE development phase only - to be deleted
-      !!
-      !! ** Method  :
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices
-
-      trb => tr(:,:,:,:,Kbb); trn => tr(:,:,:,:,Kmm); tra => tr(:,:,:,:,Kaa)
-
-   END SUBROUTINE update_pointers_trc
 
    SUBROUTINE trc_mean_qsr( kt )
