Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_oce_update.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_oce_update.F90	(revision 11947)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_oce_update.F90	(revision 11949)
@@ -230,10 +230,10 @@
       ! -----------------------
       !
-      e3u_a(:,:,:) = e3u_n(:,:,:)
-      e3v_a(:,:,:) = e3v_n(:,:,:)
-!      ua(:,:,:) = e3u_b(:,:,:)
-!      va(:,:,:) = e3v_b(:,:,:)
-      hu_a(:,:) = hu_n(:,:)
-      hv_a(:,:) = hv_n(:,:)
+      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(:,:,Krhs_a) = hu(:,:,Kmm_a)
+      hv(:,:,Krhs_a) = hv(:,:,Kmm_a)
 
       ! 1) NOW fields
@@ -242,22 +242,22 @@
          ! Vertical scale factor interpolations
          ! ------------------------------------
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) ,  'U' )
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) ,  'V' )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) ,  'F' )
-
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), '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:
          ! --------------------
-      hu_n(:,:) = 0._wp                        ! Ocean depth at U-points
-      hv_n(:,:) = 0._wp                        ! Ocean depth at V-points
+      hu(:,:,Kmm_a) = 0._wp                        ! Ocean depth at U-points
+      hv(:,:,Kmm_a) = 0._wp                        ! Ocean depth at V-points
       DO jk = 1, jpkm1
-         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
+         hu(:,:,Kmm_a) = hu(:,:,Kmm_a) + e3u(:,:,jk,Kmm_a) * umask(:,:,jk)
+         hv(:,:,Kmm_a) = hv(:,:,Kmm_a) + e3v(:,:,jk,Kmm_a) * vmask(:,:,jk)
       END DO
       !                                        ! Inverse of the local depth
-      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
-      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
+      r1_hu(:,:,Kmm_a) = ssumask(:,:) / ( hu(:,:,Kmm_a) + 1._wp - ssumask(:,:) )
+      r1_hv(:,:,Kmm_a) = ssvmask(:,:) / ( hv(:,:,Kmm_a) + 1._wp - ssvmask(:,:) )
 
 
@@ -268,21 +268,21 @@
          ! Vertical scale factor interpolations
          ! ------------------------------------
-         CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:),  'U'  )
-         CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:),  'V'  )
-
-         CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-         CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), '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:
          ! --------------------
-         hu_b(:,:) = 0._wp                     ! Ocean depth at U-points
-         hv_b(:,:) = 0._wp                     ! Ocean depth at V-points
+         hu(:,:,Kbb_a) = 0._wp                     ! Ocean depth at U-points
+         hv(:,:,Kbb_a) = 0._wp                     ! Ocean depth at V-points
          DO jk = 1, jpkm1
-            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
-            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
+            hu(:,:,Kbb_a) = hu(:,:,Kbb_a) + e3u(:,:,jk,Kbb_a) * umask(:,:,jk)
+            hv(:,:,Kbb_a) = hv(:,:,Kbb_a) + e3v(:,:,jk,Kbb_a) * vmask(:,:,jk)
          END DO
          !                                     ! Inverse of the local depth
-         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
-         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
+         r1_hu(:,:,Kbb_a) = ssumask(:,:) / ( hu(:,:,Kbb_a) + 1._wp - ssumask(:,:) )
+         r1_hv(:,:,Kbb_a) = ssvmask(:,:) / ( hv(:,:,Kbb_a) + 1._wp - ssvmask(:,:) )
       ENDIF
       !
@@ -315,6 +315,6 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) &
-                                           * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp
+                     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._wp)*999._wp
                   END DO
                END DO
@@ -324,11 +324,11 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) &
-                                           + (tmask(ji,jj,jk)-1)*999._wp
+                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) &
+                    &                      + (tmask(ji,jj,jk) - 1._wp)*999._wp
                END DO
             END DO
          END DO
       ELSE
-         tabres_child(:,:,:,:) = 0.
+         tabres_child(:,:,:,:) = 0._wp
          AGRIF_SpecialValue = 0._wp
          DO jj=j1,j2
@@ -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_n(ji,jj,jk) 
+                  h_out(N_out) = e3t(ji,jj,jk,Kmm_a) 
                ENDDO
                IF (N_in > 0) THEN !Remove this?
@@ -356,5 +356,6 @@
                   ENDIF
                   DO jn=n1,n2-1
-                     CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out)
+                     CALL reconstructandremap( tabin(1:N_in,jn), h_in(1:N_in), tabres_child(ji,jj,1:N_out,jn),  &
+                       &                       h_out(1:N_out)  , N_in        , N_out                           )
                   ENDDO
                ENDIF
@@ -369,7 +370,8 @@
                      DO ji=i1,i2
                         IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
-                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
-                                 & + atfp * ( tabres_child(ji,jj,jk,jn) &
-                                 &          - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
+                           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_a) &
+                                 &                         ) * tmask(ji,jj,jk)
                         ENDIF
                      ENDDO
@@ -383,5 +385,5 @@
                   DO ji=i1,i2
                      IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
-                        tsn(ji,jj,jk,jn) = 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 +415,6 @@
                   DO ji=i1,i2
 !> jc tmp
-                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
-!                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk)
+                     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 +436,9 @@
                      DO ji = i1, i2
                         IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
-                           ztb  = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! 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 = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
-                           tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
-                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk)
+                           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 +452,5 @@
                   DO ji=i1,i2
                      IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 
-                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)
+                        ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a)
                      END IF
                   END DO
@@ -458,5 +460,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            tsb(i1:i2,j1:j2,k1:k2,1:jpts)  = tsn(i1:i2,j1:j2,k1:k2,1:jpts)
+            ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb_a)  = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm_a)
          ENDIF
          !
@@ -495,8 +497,8 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  &
-                                       + (umask(ji,jj,jk)-1)*999._wp
-                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  &
-                                       + (umask(ji,jj,jk)-1)*999._wp
+                  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._wp)*999._wp
+                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)  &
+                    &                  + (umask(ji,jj,jk)-1._wp)*999._wp
                END DO
             END DO
@@ -513,5 +515,5 @@
                   IF( tabres(ji,jj,jk,2) < -900) EXIT
                   N_in = N_in + 1
-                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
+                  tabin(jk)  = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
                   h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
                ENDDO
@@ -520,8 +522,8 @@
                   IF (umask(ji,jj,jk) == 0) EXIT
                   N_out = N_out + 1
-                  h_out(N_out) = e3u_n(ji,jj,jk)
+                  h_out(N_out) = e3u(ji,jj,jk,Kmm_a)
                ENDDO
                IF (N_in * N_out > 0) THEN
-                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                  h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) )
                   IF (h_diff < -1.e-4) THEN
 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 
@@ -540,6 +542,7 @@
                      ENDDO
                   ENDIF
-                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
-                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
+                  CALL reconstructandremap( tabin(1:N_in) , h_in(1:N_in), tabres_child(ji,jj,1:N_out),       &
+                    &                       h_out(1:N_out), N_in        , N_out                            )
+                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/( e2u(ji,jj)*h_out(N_out) )
                ENDIF
             ENDDO
@@ -550,9 +553,9 @@
                DO ji=i1,i2
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     ub(ji,jj,jk) = ub(ji,jj,jk) & 
-                           & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * 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
                   !
-                  un(ji,jj,jk) = 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 +582,6 @@
          zrhoy = Agrif_Rhoy()
          DO jk = k1, k2
-            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
+            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 +592,12 @@
                   !
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)  ! fse3t_b prior update should be used
-                     zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk)
+                     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)
-                     ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &      
-                                    & * umask(ji,jj,jk) / e3u_b(ji,jj,jk)
+                     uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno ) )       &      
+                       &                  * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)
                   ENDIF
                   !
-                  un(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)
+                  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 +605,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2)
+            uu(i1:i2,j1:j2,k1:k2,Kbb_a)  = uu(i1:i2,j1:j2,k1:k2,Kmm_a)
          ENDIF
          !
@@ -632,8 +636,8 @@
          IF (western_side) THEN
             DO jj=j1,j2
-               zcor = un_b(i1-1,jj) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - un_b(i1-1,jj)
-               un_b(i1-1,jj) = un_b(i1-1,jj) + zcor
+               zcor = uu_b(i1-1,jj,Kmm_a) * hu(i1-1,jj,Krhs_a) * r1_hu(i1-1,jj,Kmm_a) - 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
-                  un(i1-1,jj,jk) = un(i1-1,jj,jk) + 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 +646,8 @@
          IF (eastern_side) THEN
             DO jj=j1,j2
-               zcor = un_b(i2+1,jj) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - un_b(i2+1,jj)
-               un_b(i2+1,jj) = un_b(i2+1,jj) + zcor
+               zcor = uu_b(i2+1,jj,Kmm_a) * hu(i2+1,jj,Krhs_a) * r1_hu(i2+1,jj,Kmm_a) - 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
-                  un(i2+1,jj,jk) = un(i2+1,jj,jk) + 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 +686,8 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) &
-                                       + (vmask(ji,jj,jk)-1)*999._wp
-                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) &
+                  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_a) &
                                        + (vmask(ji,jj,jk)-1)*999._wp
                END DO
@@ -705,8 +710,8 @@
                   IF (vmask(ji,jj,jk) == 0) EXIT
                   N_out = N_out + 1
-                  h_out(N_out) = e3v_n(ji,jj,jk)
+                  h_out(N_out) = e3v(ji,jj,jk,Kmm_a)
                ENDDO
                IF (N_in * N_out > 0) THEN
-                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                  h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) )
                   IF (h_diff < -1.e-4) then
 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 
@@ -714,7 +719,7 @@
                      excess = 0._wp
                      DO jk=N_in,1,-1
-                        thick = MIN(-1*h_diff, h_in(jk))
+                        thick = MIN( -1._wp*h_diff, h_in(jk) )
                         excess = excess + tabin(jk)*thick*e2u(ji,jj)
-                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
+                        tabin(jk) = tabin(jk)*(1._wp - thick/h_in(jk))
                         h_diff = h_diff + thick
                         IF ( h_diff == 0) THEN
@@ -725,5 +730,6 @@
                      ENDDO
                   ENDIF
-                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
+                  CALL reconstructandremap( tabin(1:N_in) , h_in(1:N_in), tabres_child(ji,jj,1:N_out),     &
+                    &                       h_out(1:N_out), N_in        , N_out                          )
                   tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
                ENDIF
@@ -736,9 +742,10 @@
                   !
                   IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part
-                     vb(ji,jj,jk) = vb(ji,jj,jk) & 
-                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * 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
                   !
-                  vn(ji,jj,jk) = 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 +774,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
+                  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 +785,12 @@
                   !
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
-                     zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used
-                     zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk)
+                     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)
-                     vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &      
-                                    & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk)
+                     vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) )        &      
+                       &                  * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)
                   ENDIF
                   !
-                  vn(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)
+                  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 +798,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2)
+            vv(i1:i2,j1:j2,k1:k2,Kbb_a)  = vv(i1:i2,j1:j2,k1:k2,Kmm_a)
          ENDIF
          !
@@ -822,8 +829,8 @@
          IF (southern_side) THEN
             DO ji=i1,i2
-               zcor = vn_b(ji,j1-1) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vn_b(ji,j1-1)
-               vn_b(ji,j1-1) = vn_b(ji,j1-1) + zcor
+               zcor = vv_b(ji,j1-1,Kmm_a) * hv(ji,j1-1,Krhs_a) * r1_hv(ji,j1-1,Kmm_a) - 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
-                  vn(ji,j1-1,jk) = vn(ji,j1-1,jk) + 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 +839,8 @@
          IF (northern_side) THEN
             DO ji=i1,i2
-               zcor = vn_b(ji,j2+1) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vn_b(ji,j2+1)
-               vn_b(ji,j2+1) = vn_b(ji,j2+1) + zcor
+               zcor = vv_b(ji,j2+1,Kmm_a) * hv(ji,j2+1,Krhs_a) * r1_hv(ji,j2+1,Kmm_a) - 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
-                  vn(ji,j2+1,jk) = vn(ji,j2+1,jk) + 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 +869,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)
+               tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm_a) * hu(ji,jj,Kmm_a) * e2u(ji,jj)
             END DO
          END DO
@@ -873,30 +880,30 @@
                spgu(ji,jj) = 0._wp
                DO jk=1,jpkm1
-                  spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
+                  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)
+               zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu(ji,jj,Kmm_a)
                DO jk=1,jpkm1              
-                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
+                  uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk)           
                END DO
                !
                ! Update barotropic velocities:
                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) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
-                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
+                  IF ( .NOT.( lk_agrif_fstep .AND. (neuler==0) ) ) THEN                          ! Add asselin part
+                     zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu(ji,jj,Krhs_a)) * r1_hu(ji,jj,Kbb_a)
+                     uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + atfp * zcorr * umask(ji,jj,1)
                   END IF
                ENDIF    
-               un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
+               uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * 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_b(ji,jj,jk) * ub(ji,jj,jk)
+                  spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a)
                END DO
                !
-               zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj)
+               zcorr = uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu(ji,jj,Kbb_a)
                DO jk=1,jpkm1              
-                  ub(ji,jj,jk) = ub(ji,jj,jk) + 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 +912,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2)
+            uu_b(i1:i2,j1:j2,Kbb_a)  = uu_b(i1:i2,j1:j2,Kmm_a)
          ENDIF
       ENDIF
@@ -928,5 +935,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj) 
+               tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm_a) * hv(ji,jj,Kmm_a) * e1v(ji,jj) 
             END DO
          END DO
@@ -939,30 +946,30 @@
                spgv(ji,jj) = 0.e0
                DO jk=1,jpkm1
-                  spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
+                  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)
+               zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv(ji,jj,Kmm_a)
                DO jk=1,jpkm1              
-                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
+                  vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk)           
                END DO
                !
                ! Update barotropic velocities:
-               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) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
-                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
+               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_a) * hv(ji,jj,Krhs_a)) * r1_hv(ji,jj,Kbb_a)
+                     vv_b(ji,jj,Kbb_a) =       vv_b(ji,jj,Kbb_a) + atfp * zcorr  *   vmask(ji,jj,1)
                   END IF
                ENDIF              
-               vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
+               vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * 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_b(ji,jj,jk) * vb(ji,jj,jk)
+                  spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a)
                END DO
                !
-               zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj)
+               zcorr = vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv(ji,jj,Kbb_a)
                DO jk=1,jpkm1              
-                  vb(ji,jj,jk) = vb(ji,jj,jk) + 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 +978,5 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2)
+            vv_b(i1:i2,j1:j2,Kbb_a)  = vv_b(i1:i2,j1:j2,Kmm_a)
          ENDIF
          !
@@ -993,5 +1000,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) = sshn(ji,jj)
+               tabres(ji,jj) = ssh(ji,jj,Kmm_a)
             END DO
          END DO
@@ -1000,6 +1007,7 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  sshb(ji,jj) =   sshb(ji,jj) &
-                        & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * 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 +1016,10 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               sshn(ji,jj) = 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
-            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2)
+            ssh(i1:i2,j1:j2,Kbb_a)  = ssh(i1:i2,j1:j2,Kmm_a)
          ENDIF
          !
@@ -1094,6 +1102,7 @@
             DO jj=j1,j2
                zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 
-               sshn(i1  ,jj) = sshn(i1  ,jj) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + 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 +1110,7 @@
             DO jj=j1,j2
                zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
-               sshn(i2+1,jj) = sshn(i2+1,jj) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + 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
@@ -1181,14 +1191,16 @@
          IF (southern_side) THEN
             DO ji=i1,i2
-               zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
-               sshn(ji,j1  ) = sshn(ji,j1  ) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor
+               zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * ( vb2_b(ji,j1)-tabres(ji,j1) )
+               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
          IF (northern_side) THEN               
             DO ji=i1,i2
-               zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
-               sshn(ji,j2+1) = sshn(ji,j2+1) + zcor
-               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor
+               zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * ( vb2_b(ji,j2)-tabres(ji,j2) )
+               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
@@ -1223,7 +1235,7 @@
             END DO
          END DO
-         tabres(:,:,:,1)=tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy()
-         tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
-         tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
+         tabres(:,:,:,1) = tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy()
+         tabres(:,:,:,2) = tabres(:,:,:,2)*Agrif_Rhox()
+         tabres(:,:,:,3) = tabres(:,:,:,3)*Agrif_Rhoy()
       ELSE
          DO jk=k1,k2
@@ -1234,10 +1246,10 @@
                      print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk)
                      print *,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)*tmask(ji,jj,jk)
-                     ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
+                     ztemp = SQRT( tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)) )
                      print *,'CORR = ',ztemp-1.
                      print *,'NEW VALS = ',tabres(ji,jj,jk,2)*ztemp,tabres(ji,jj,jk,3)*ztemp, &
                            tabres(ji,jj,jk,2)*ztemp*tabres(ji,jj,jk,3)*ztemp
-                     e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
-                     e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
+                     e1t(ji,jj) = tabres(ji,jj,jk,2) * ztemp
+                     e2t(ji,jj) = tabres(ji,jj,jk,3) * ztemp
                   END IF
                END DO
@@ -1319,6 +1331,8 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + sshn(ji,jj) &
-                                     & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
+                  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
             END DO
@@ -1330,8 +1344,8 @@
          ! Save "old" scale factor (prior update) for subsequent asselin correction
          ! of prognostic variables
-         e3t_a(i1:i2,j1:j2,1:jpkm1) = e3t_n(i1:i2,j1:j2,1:jpkm1)
-
-         ! One should also save e3t_b, but lacking of workspace...
-!         hdivn(i1:i2,j1:j2,1:jpkm1)   = e3t_b(i1:i2,j1:j2,1:jpkm1)
+         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 +1353,13 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     e3t_b(ji,jj,jk) =  e3t_b(ji,jj,jk) &
-                           & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) )
+                     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_b  (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_b(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
-            gdepw_b(i1:i2,j1:j2,1) = 0.0_wp
-            gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1)
+            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 +1367,11 @@
                   DO ji = i1,i2            
                      zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-                     e3w_b(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
-                     &                                        ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )  &
-                     &                                  +            0.5_wp * tmask(ji,jj,jk)   *        &
-                     &                                        ( e3t_b(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
-                     gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
-                     gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
-                         &               + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
+                     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_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._wp - zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) +       e3w(ji,jj,jk  ,Kbb_a) ) 
                   END DO
                END DO
@@ -1370,17 +1384,17 @@
          !
          ! Update vertical scale factor at T-points:
-         e3t_n(i1:i2,j1:j2,1:jpkm1) = 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
+         ht(i1:i2,j1:j2) = 0._wp
          DO jk = 1, jpkm1
-            ht_n(i1:i2,j1:j2) = ht_n(i1:i2,j1:j2) + e3t_n(i1:i2,j1:j2,jk) * tmask(i1:i2,j1:j2,jk)
+            ht(i1:i2,j1:j2) = ht(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_n (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_n(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
-         gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1)
-         gdepw_n(i1:i2,j1:j2,1) = 0.0_wp
-         gde3w_n(i1:i2,j1:j2,1) = gdept_n(i1:i2,j1:j2,1) - (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(i1:i2,j1:j2) - ht_0(i1:i2,j1:j2) ) ! Last term in the rhs is ssh
          !
          DO jk = 2, jpk
@@ -1388,10 +1402,12 @@
                DO ji = i1,i2            
                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               e3w_n(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t_n(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )   &
-               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t_n(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
-               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
-                   &               + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
-               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - (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._wp - 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(ji,jj)-ht_0(ji,jj) )    ! Last term in the rhs is ssh
                END DO
             END DO
@@ -1399,8 +1415,8 @@
          !
          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
-            e3t_b (i1:i2,j1:j2,1:jpk)  = e3t_n (i1:i2,j1:j2,1:jpk)
-            e3w_b (i1:i2,j1:j2,1:jpk)  = e3w_n (i1:i2,j1:j2,1:jpk)
-            gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk)
-            gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk)
+              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
          !
