Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/cfgs/SHARED/field_def_nemo-oce.xml	(revision 11089)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/cfgs/SHARED/field_def_nemo-oce.xml	(revision 11099)
@@ -62,4 +62,14 @@
         <field id="mldr10_1max"  long_name="Max of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="maximum"                                                                          />
         <field id="mldr10_1min"  long_name="Min of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="minimum"                                                                          />
+         <field id="mldzint_1"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          />
+         <field id="mldzint_2"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          />
+         <field id="mldzint_3"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          />
+         <field id="mldzint_4"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          />
+         <field id="mldzint_5"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          />
+         <field id="mldhtc_1"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
+         <field id="mldhtc_2"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
+         <field id="mldhtc_3"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
+         <field id="mldhtc_4"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
+         <field id="mldhtc_5"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
         <field id="heatc"        long_name="Heat content vertically integrated"                 standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
         <field id="saltc"        long_name="Salt content vertically integrated"                                                                                                                   unit="1e-3*kg/m2" />
@@ -77,10 +87,4 @@
         <!-- variables available with MLE -->
         <field id="Lf_NHpf"      long_name="MLE: Lf = N H / f"   unit="m" />
-
-        <!-- next variables available with ln_zad_Aimp=.true. -->
-        <field id="Courant"      long_name="Courant number"                                                                 unit="#"   grid_ref="grid_T_3D" />
-        <field id="wimp"         long_name="Implicit vertical velocity"                                                     unit="m/s" grid_ref="grid_T_3D" />
-        <field id="wexp"         long_name="Explicit vertical velocity"                                                     unit="m/s" grid_ref="grid_T_3D" />
-        <field id="wi_cff"       long_name="Fraction of implicit vertical velocity"                                         unit="#"   grid_ref="grid_T_3D" />
 
         <!-- next variables available with key_diahth -->
@@ -351,4 +355,5 @@
         <field id="uoce"         long_name="ocean current along i-axis"                             standard_name="sea_water_x_velocity"        unit="m/s"        grid_ref="grid_U_3D" />
         <field id="uoce_e3u"     long_name="ocean current along i-axis  (thickness weighted)"                                                   unit="m/s"        grid_ref="grid_U_3D"  > uoce * e3u </field>
+        <field id="uoce2_e3u"    long_name="ocean current along i-axis squared (thickness weighted)"                                            unit="m3/s2"      grid_ref="grid_U_3D"  > uoce * uoce * e3u </field>
         <field id="ssu"          long_name="ocean surface current along i-axis"                                                                 unit="m/s"                             />
         <field id="sbu"          long_name="ocean bottom current along i-axis"                                                                  unit="m/s"                             />
@@ -405,4 +410,5 @@
         <field id="voce"         long_name="ocean current along j-axis"                             standard_name="sea_water_y_velocity"        unit="m/s"        grid_ref="grid_V_3D" />
         <field id="voce_e3v"     long_name="ocean current along j-axis  (thickness weighted)"                                                   unit="m/s"        grid_ref="grid_V_3D"  > voce * e3v </field>
+        <field id="voce2_e3v"    long_name="ocean current along j-axis squared (thickness weighted)"                                            unit="m3/s2"      grid_ref="grid_V_3D"  > voce * voce * e3v </field>
         <field id="ssv"          long_name="ocean surface current along j-axis"                                                                 unit="m/s"                             />
         <field id="sbv"          long_name="ocean bottom current along j-axis"                                                                  unit="m/s"                             />
@@ -975,5 +981,5 @@
 
     <field_group id="25h_grid_W" grid_ref="grid_W_3D" operation="instant">
-      <field id="vovecrtz25h"         name="k current 25h mean"                 unit="m/s"      />
+      <field id="vomecrtz25h"         name="k current 25h mean"                 unit="m/s"      />
       <field id="avt25h"              name="vertical diffusivity25h mean"       unit="m2/s" />
       <field id="avm25h"              name="vertical viscosity 25h mean"        unit="m2/s" />
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_oce_interp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_oce_interp.F90	(revision 11089)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_oce_interp.F90	(revision 11099)
@@ -115,5 +115,5 @@
             END DO
             DO jj = 1, jpj
-               uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu_a(ibdy1:ibdy2,jj)
+               uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu(ibdy1:ibdy2,jj,Krhs_a)
             END DO
          ENDIF
@@ -135,5 +135,5 @@
          END DO
          DO jj=1,jpj
-            zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
+            zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu(ibdy1:ibdy2,jj,Krhs_a)
          END DO
             
@@ -154,5 +154,5 @@
             END DO
             DO jj = 1, jpj
-               zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj)
+               zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv(ibdy1:ibdy2,jj,Krhs_a)
             END DO
             DO jk = 1, jpkm1
@@ -186,5 +186,5 @@
             END DO
             DO jj = 1, jpj
-               uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu_a(ibdy1:ibdy2,jj)
+               uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu(ibdy1:ibdy2,jj,Krhs_a)
             END DO
          ENDIF
@@ -206,5 +206,5 @@
          END DO
          DO jj=1,jpj
-            zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
+            zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu(ibdy1:ibdy2,jj,Krhs_a)
          END DO
             
@@ -227,5 +227,5 @@
             END DO
             DO jj = 1, jpj
-               zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj)
+               zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv(ibdy1:ibdy2,jj,Krhs_a)
             END DO
             DO jk = 1, jpkm1
@@ -259,5 +259,5 @@
             END DO
             DO ji=1,jpi
-               vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv_a(ji,jbdy1:jbdy2)
+               vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv(ji,jbdy1:jbdy2,Krhs_a)
             END DO
          ENDIF
@@ -279,5 +279,5 @@
          END DO
          DO ji = 1, jpi
-            zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
+            zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv(ji,jbdy1:jbdy2,Krhs_a)
          END DO
 
@@ -298,5 +298,5 @@
             END DO
             DO ji = 1, jpi
-               zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2)
+               zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu(ji,jbdy1:jbdy2,Krhs_a)
             END DO
                
@@ -331,5 +331,5 @@
             END DO
             DO ji=1,jpi
-               vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv_a(ji,jbdy1:jbdy2)
+               vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv(ji,jbdy1:jbdy2,Krhs_a)
             END DO
          ENDIF
@@ -351,5 +351,5 @@
          END DO
          DO ji = 1, jpi
-            zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
+            zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv(ji,jbdy1:jbdy2,Krhs_a)
          END DO
 
@@ -372,5 +372,5 @@
             END DO
             DO ji = 1, jpi
-               zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2)
+               zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu(ji,jbdy1:jbdy2,Krhs_a)
             END DO
                
@@ -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_a)
+         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu(i1:i2,j1:j2,Kmm_a) * 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_a)
+         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv(i1:i2,j1:j2,Kmm_a) * vv_b(i1:i2,j1:j2,Kmm_a)
       ELSE
          western_side  = (nb == 1).AND.(ndir == 1)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_oce_update.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_oce_update.F90	(revision 11089)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_oce_update.F90	(revision 11099)
@@ -234,6 +234,6 @@
 !      uu(:,:,:,Krhs_a) = e3u(:,:,:,Kbb_a)
 !      vv(:,:,:,Krhs_a) = e3v(:,:,:,Kbb_a)
-      hu_a(:,:) = hu_n(:,:)
-      hv_a(:,:) = hv_n(:,:)
+      hu(:,:,Krhs_a) = hu(:,:,Kmm_a)
+      hv(:,:,Krhs_a) = hv(:,:,Kmm_a)
 
       ! 1) NOW fields
@@ -251,13 +251,13 @@
          ! 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(:,:,jk,Kmm_a) * umask(:,:,jk)
-         hv_n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm_a) * 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(:,:) )
 
 
@@ -276,13 +276,13 @@
          ! 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(:,:,jk,Kbb_a) * umask(:,:,jk)
-            hv_b(:,:) = hv_b(:,:) + e3v(:,:,jk,Kbb_a) * 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
       !
@@ -632,5 +632,5 @@
          IF (western_side) THEN
             DO jj=j1,j2
-               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)
+               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
@@ -642,5 +642,5 @@
          IF (eastern_side) THEN
             DO jj=j1,j2
-               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)
+               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
@@ -822,5 +822,5 @@
          IF (southern_side) THEN
             DO ji=i1,i2
-               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)
+               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
@@ -832,5 +832,5 @@
          IF (northern_side) THEN
             DO ji=i1,i2
-               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)
+               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
@@ -862,5 +862,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm_a) * 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
@@ -876,5 +876,5 @@
                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              
                   uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk)           
@@ -884,9 +884,9 @@
                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_a) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
+                     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    
-               uu_b(ji,jj,Kmm_a) = 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:
@@ -896,5 +896,5 @@
                END DO
                !
-               zcorr = uu_b(ji,jj,Kbb_a) - 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              
                   uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk)           
@@ -928,5 +928,5 @@
          DO jj=j1,j2
             DO ji=i1,i2
-               tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm_a) * 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
@@ -942,5 +942,5 @@
                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              
                   vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk)           
@@ -950,9 +950,9 @@
                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_a(ji,jj)) * r1_hv_b(ji,jj)
+                     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              
-               vv_b(ji,jj,Kmm_a) = 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:
@@ -962,5 +962,5 @@
                END DO
                !
-               zcorr = vv_b(ji,jj,Kbb_a) - 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              
                   vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk)           
@@ -1373,7 +1373,7 @@
          !
          ! 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(i1:i2,j1:j2,jk,Kmm_a) * 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
          !
@@ -1382,5 +1382,5 @@
          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
+         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
@@ -1393,5 +1393,5 @@
                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
+               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
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/ASM/asminc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/ASM/asminc.F90	(revision 11089)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/ASM/asminc.F90	(revision 11099)
@@ -805,5 +805,5 @@
       ELSE 
          ALLOCATE( ztim(jpi,jpj) )
-         ztim(:,:) = ssh_iau(:,:) / ( ht_n(:,:) + 1.0 - ssmask(:,:) )
+         ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) )
          DO jk = 1, jpkm1                                 
             phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DOM/domain.F90	(revision 11089)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DOM/domain.F90	(revision 11099)
@@ -163,6 +163,6 @@
       !
          !       before        !          now          !       after         !
-            gdept(:,:,:,Kbb) = gdept_0  ;   gdept(:,:,:,Kmm) = gdept_0   !        ---          ! depth of grid-points
-            gdepw(:,:,:,Kbb) = gdepw_0  ;   gdepw(:,:,:,Kmm) = gdepw_0   !        ---          !
+            gdept(:,:,:,Kbb) = gdept_0  ;   gdept(:,:,:,Kmm) = gdept_0   ;   gdept(:,:,:,Kaa) = gdept_0   ! depth of grid-points
+            gdepw(:,:,:,Kbb) = gdepw_0  ;   gdepw(:,:,:,Kmm) = gdepw_0   ;   gdepw(:,:,:,Kaa) = gdepw_0   !
                                    gde3w = gde3w_0   !        ---          !
          !                                                                  
@@ -171,7 +171,7 @@
               e3v(:,:,:,Kbb) =   e3v_0  ;     e3v(:,:,:,Kmm) =   e3v_0   ;   e3v(:,:,:,Kaa) =  e3v_0    !
                                      e3f =   e3f_0   !        ---          !
-              e3w(:,:,:,Kbb) =   e3w_0  ;     e3w(:,:,:,Kmm) =   e3w_0   !        ---          !
-             e3uw(:,:,:,Kbb) =  e3uw_0  ;    e3uw(:,:,:,Kmm) =  e3uw_0   !        ---          !
-             e3vw(:,:,:,Kbb) =  e3vw_0  ;    e3vw(:,:,:,Kmm) =  e3vw_0   !        ---          !
+              e3w(:,:,:,Kbb) =   e3w_0  ;     e3w(:,:,:,Kmm) =   e3w_0   ;    e3w(:,:,:,Kaa) =   e3w_0   ! 
+             e3uw(:,:,:,Kbb) =  e3uw_0  ;    e3uw(:,:,:,Kmm) =  e3uw_0   ;   e3uw(:,:,:,Kaa) =  e3uw_0   !  
+             e3vw(:,:,:,Kbb) =  e3vw_0  ;    e3vw(:,:,:,Kmm) =  e3vw_0   ;   e3vw(:,:,:,Kaa) =  e3vw_0   !
          !
          z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/traatf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/traatf.F90	(revision 11089)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/traatf.F90	(revision 11099)
@@ -146,10 +146,6 @@
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
-         DO jn = 1, jpts
-            DO jk = 1, jpkm1
-               pts(:,:,jk,jn,Kmm) = pts(:,:,jk,jn,Kaa)    
-            END DO
-         END DO
+      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping 
+         !
          IF (l_trdtra .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
             !                                        ! Asselin filter is output by tra_atf_vvl that is not called on this time step
@@ -175,6 +171,6 @@
          zfact = 1._wp / r2dt             
          DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kbb) - ztrdt(:,:,jk) ) * zfact
-            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kbb) - ztrds(:,:,jk) ) * zfact
+            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * zfact
+            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kmm) - ztrds(:,:,jk) ) * zfact
          END DO
          CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt )
@@ -192,5 +188,5 @@
 
 
-   SUBROUTINE tra_atf_fix( kt, kit000, Kbb, Kmm, Kaa, cdtype, pt, kjpt )
+   SUBROUTINE tra_atf_fix( kt, Kbb, Kmm, Kaa, kit000, cdtype, pt, kjpt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_atf_fix  ***
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/trasbc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/trasbc.F90	(revision 11089)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/trasbc.F90	(revision 11099)
@@ -233,5 +233,5 @@
             DO jj = 2, jpj 
                DO ji = fs_2, fs_jpim1
-                  ztim = ssh_iau(ji,jj) / ( ht_n(ji,jj) + 1. - ssmask(ji, jj) )
+                  ztim = ssh_iau(ji,jj) / ( ht(ji,jj) + 1. - ssmask(ji, jj) )
                   pts(ji,jj,:,jp_tem,Krhs) = pts(ji,jj,:,jp_tem,Krhs) + pts(ji,jj,:,jp_tem,Kmm) * ztim
                   pts(ji,jj,:,jp_sal,Krhs) = pts(ji,jj,:,jp_sal,Krhs) + pts(ji,jj,:,jp_sal,Kmm) * ztim
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/PISCES/SED/trcdmp_sed.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/PISCES/SED/trcdmp_sed.F90	(revision 11089)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/PISCES/SED/trcdmp_sed.F90	(revision 11099)
@@ -149,6 +149,7 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trc_dmp_sed( kt )        ! Empty routine
+   SUBROUTINE trc_dmp_sed( kt, Kbb, Kmm, Krhs )   ! Empty routine
       INTEGER, INTENT(in) :: kt
+      INTEGER, INTENT(in) :: Kbb, Kmm, Krhs
       WRITE(*,*) 'trc_dmp_sed: You should not have seen this print! error?', kt
    END SUBROUTINE trc_dmp_sed
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcatf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcatf.F90	(revision 11089)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcatf.F90	(revision 11099)
@@ -30,5 +30,5 @@
    USE trd_oce
    USE trdtra
-   USE tranxt
+   USE traatf
    USE bdy_oce   , ONLY: ln_bdy
    USE trcbdy          ! BDY open boundaries
@@ -143,11 +143,21 @@
       ENDIF
       !                                ! Leap-Frog + Asselin filter time stepping
-      IF( .NOT.( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) ) THEN    ! Only time filter if not an Euler timestep
+      IF( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) THEN    ! Euler time-stepping 
+         !
+         IF (l_trdtrc .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
+            !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
+            ztrdt(:,:,:,:) = 0._wp            
+            DO jn = 1, jptra
+               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
+            ENDDO
+         END IF
+         !
+      ELSE     
          IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
-            IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, Kbb, Kmm, Kaa, nittrc000,         'TRC', ptr, jptra )                     !     linear ssh
-            ELSE                   ;   CALL tra_nxt_vvl( kt, Kbb, Kmm, Kaa, nittrc000, rdttrc, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
+            IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nittrc000,         'TRC', ptr, jptra )                     !     linear ssh
+            ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nittrc000, rdttrc, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
             ENDIF
          ELSE
-                                       CALL trc_atf_off( kt, Kbb, Kmm, Kaa )       ! offline 
+                                       CALL trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )       ! offline 
          ENDIF
          !
@@ -260,9 +270,11 @@
    !!   Default option                                         Empty module
    !!----------------------------------------------------------------------
+   USE par_oce
+   USE par_trc
 CONTAINS
    SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr )  
-      INTEGER                                   , INTENT(in) :: kt
-      INTEGER                                   , INTENT(in) :: Kbb, Kmm, Kaa  ! time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr            ! passive tracers
+      INTEGER                                   , INTENT(in)    :: kt
+      INTEGER,                                    INTENT(in   ) :: Kbb, Kmm, Kaa ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr           ! passive tracers and RHS of tracer equation
       WRITE(*,*) 'trc_atf: You should not have seen this print! error?', kt
    END SUBROUTINE trc_atf
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcsbc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcsbc.F90	(revision 11089)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/TRP/trcsbc.F90	(revision 11099)
@@ -197,4 +197,6 @@
    !!   Dummy module :                      NO passive tracer
    !!----------------------------------------------------------------------
+   USE par_oce
+   USE par_trc
 CONTAINS
    SUBROUTINE trc_sbc ( kt, Kmm, ptr, Krhs )      ! Empty routine
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/tests/VORTEX/MY_SRC/domvvl.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/tests/VORTEX/MY_SRC/domvvl.F90	(revision 11089)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/tests/VORTEX/MY_SRC/domvvl.F90	(revision 11099)
@@ -106,8 +106,8 @@
       !!              - Regrid: e3(u/v)_n
       !!                        e3(u/v)_b       
-      !!                        e3w_n           
+      !!                        e3w(:,:,:,Kmm)           
       !!                        e3(u/v)w_b      
       !!                        e3(u/v)w_n      
-      !!                        gdept_n, gdepw_n and gde3w_n
+      !!                        gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w
       !!              - h(t/u/v)_0
       !!              - frq_rst_e3t and frq_rst_hdv
@@ -131,32 +131,32 @@
       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
       CALL dom_vvl_rst( nit000, 'READ' )
-      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
+      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
       !
       !                    !== Set of all other vertical scale factors  ==!  (now and before)
       !                                ! Horizontal interpolation of e3t
-      CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )    ! from T to U
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-      CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )    ! from T to V 
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F
+      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V 
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F
       !                                ! Vertical interpolation of e3t,u,v 
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )  ! from T to W
-      CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )  ! from U to UW
-      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )  ! from V to UW
-      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W
+      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW
+      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
+      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW
+      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
 
       ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
-      e3t_a(:,:,:) = e3t_n(:,:,:)
-      e3u_a(:,:,:) = e3u_n(:,:,:)
-      e3v_a(:,:,:) = e3v_n(:,:,:)
+      e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
+      e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
+      e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
       !
       !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
-      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration)
-      gdepw_n(:,:,1) = 0.0_wp
-      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)  ! reference to a common level z=0 for hpg
-      gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1)
-      gdepw_b(:,:,1) = 0.0_wp
+      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)       ! reference to the ocean surface (used for MLD and light penetration)
+      gdepw(:,:,1,Kmm) = 0.0_wp
+      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
+      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb)
+      gdepw(:,:,1,Kbb) = 0.0_wp
       DO jk = 2, jpk                               ! vertical sum
          DO jj = 1,jpj
@@ -165,13 +165,13 @@
                !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
                !                             ! 0.5 where jk = mikt     
-!!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
+!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ??
                zcoef = ( tmask(ji,jj,jk) - wmask(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) - sshn(ji,jj)
-               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)) 
+               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) - ssh(ji,jj,Kmm)
+               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)) 
             END DO
          END DO
@@ -179,22 +179,22 @@
       !
       !                    !==  thickness of the water column  !!   (ocean portion only)
-      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
-      hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
-      hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)
-      hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
-      hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)
+      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
+      hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1)
+      hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1)
+      hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1)
+      hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1)
       DO jk = 2, jpkm1
-         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
-         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-         hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
-         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
+         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+         hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
+         hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
+         hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
+         hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
       END DO
       !
       !                    !==  inverse of water column thickness   ==!   (u- and v- points)
-      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
-      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
-      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
-      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
+      r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
+      r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
+      r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
+      r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
 
       !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies)
@@ -321,8 +321,8 @@
       !                                                ! --------------------------------------------- !
       !
-      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
+      z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
       DO jk = 1, jpkm1
-         ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
-         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+         ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0)
+         e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
       END DO
       !
@@ -337,6 +337,6 @@
          zht(:,:)   = 0._wp
          DO jk = 1, jpkm1
-            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
-            zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
+            zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
+            zht  (:,:) = zht  (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
          END DO
          zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
@@ -348,5 +348,5 @@
                DO jk = 1, jpkm1
                   hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
-                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
+                     &          * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) )
                END DO
             ENDIF
@@ -361,9 +361,9 @@
          IF( ln_vvl_ztilde ) THEN     ! z_tilde case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
             END DO
          ELSE                         ! layer case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
             END DO
          ENDIF
@@ -476,7 +476,7 @@
             zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
-         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
+         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
          DO jk = 1, jpkm1
-            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
          END DO
 
@@ -486,5 +486,5 @@
       !                                           ! ---baroclinic part--------- !
          DO jk = 1, jpkm1
-            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
+            e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
       ENDIF
@@ -501,37 +501,37 @@
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
+            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+         END DO
+         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax
          !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
+            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk)
+         END DO
+         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax
          !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
+            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk)
+         END DO
+         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kbb) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kmm) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kaa) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax
       END IF
 
@@ -540,6 +540,6 @@
       ! *********************************** !
 
-      CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
-      CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
 
       ! *********************************** !
@@ -547,14 +547,14 @@
       ! *********************************** !
 
-      hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
-      hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
+      hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
+      hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
       DO jk = 2, jpkm1
-         hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
-         hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
+         hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
+         hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
       END DO
       !                                        ! Inverse of the local depth
 !!gm BUG ?  don't understand the use of umask_i here .....
-      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
-      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
+      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
+      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
       !
       IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
@@ -578,8 +578,8 @@
       !!               - Recompute:
       !!                    e3(u/v)_b       
-      !!                    e3w_n           
+      !!                    e3w(:,:,:,Kmm)           
       !!                    e3(u/v)w_b      
       !!                    e3(u/v)w_n      
-      !!                    gdept_n, gdepw_n  and gde3w_n
+      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
       !!                    h(u/v) and h(u/v)r
       !!
@@ -615,10 +615,10 @@
          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
       ENDIF
-      gdept_b(:,:,:) = gdept_n(:,:,:)
-      gdepw_b(:,:,:) = gdepw_n(:,:,:)
-
-      e3t_n(:,:,:) = e3t_a(:,:,:)
-      e3u_n(:,:,:) = e3u_a(:,:,:)
-      e3v_n(:,:,:) = e3v_a(:,:,:)
+      gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)
+      gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
+
+      e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa)
+      e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa)
+      e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa)
 
       ! Compute all missing vertical scale factor and depths
@@ -626,21 +626,21 @@
       ! Horizontal scale factor interpolations
       ! --------------------------------------
-      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
-      ! - JC - hu_b, hv_b, hur_b, hvr_b also
+      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are allready computed in dynnxt
+      ! - JC - hu(:,:,Kbb), hv(:,:,Kbb), hur_b, hvr_b also
       
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  )
       
       ! Vertical scale factor interpolations
-      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
-      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
+      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
+      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
 
       ! t- and w- points depth (set the isf depth as it is in the initial step)
-      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-      gdepw_n(:,:,1) = 0.0_wp
-      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
+      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
+      gdepw(:,:,1,Kmm) = 0.0_wp
+      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
       DO jk = 2, jpk
          DO jj = 1,jpj
@@ -649,8 +649,8 @@
                                                                  ! 1 for jk = mikt
                zcoef = (tmask(ji,jj,jk) - wmask(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) - sshn(ji,jj)
+               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) - ssh(ji,jj,Kmm)
             END DO
          END DO
@@ -659,10 +659,10 @@
       ! Local depth and Inverse of the local depth of the water
       ! -------------------------------------------------------
-      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
-      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
-      !
-      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
+      hu(:,:,Kmm) = hu(:,:,Kaa)   ;   r1_hu(:,:,Kmm) = r1_hu(:,:,Kaa)
+      hv(:,:,Kmm) = hv(:,:,Kaa)   ;   r1_hv(:,:,Kmm) = r1_hv(:,:,Kaa)
+      !
+      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
       DO jk = 2, jpkm1
-         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
+         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
       END DO
 
@@ -806,5 +806,5 @@
          IF( ln_rstart ) THEN                   !* Read the restart file
             CALL rst_read_open                  !  open the restart file if necessary
-            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    )
+            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
             !
             id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
@@ -817,39 +817,39 @@
             !                             ! --------- !
             IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
                ! needed to restart if land processor not computed 
-               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
+               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
                WHERE ( tmask(:,:,:) == 0.0_wp ) 
-                  e3t_n(:,:,:) = e3t_0(:,:,:)
-                  e3t_b(:,:,:) = e3t_0(:,:,:)
+                  e3t(:,:,:,Kmm) = e3t_0(:,:,:)
+                  e3t(:,:,:,Kbb) = e3t_0(:,:,:)
                END WHERE
                IF( neuler == 0 ) THEN
-                  e3t_b(:,:,:) = e3t_n(:,:,:)
+                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
                ENDIF
             ELSE IF( id1 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
                IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
-               e3t_n(:,:,:) = e3t_b(:,:,:)
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
+               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
                neuler = 0
             ELSE IF( id2 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
                IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
-               e3t_b(:,:,:) = e3t_n(:,:,:)
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
+               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
                neuler = 0
             ELSE
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
                IF(lwp) write(numout,*) 'Compute scale factor from sshn'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
                DO jk = 1, jpk
-                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
+                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
                       &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
                       &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
                END DO
-               e3t_b(:,:,:) = e3t_n(:,:,:)
+               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
                neuler = 0
             ENDIF
@@ -888,13 +888,13 @@
                IF( cn_cfg == 'wad' ) THEN
                   ! Wetting and drying test case
-                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )
-                  tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
-                  sshn (:,:)     = sshb(:,:)
-                  un   (:,:,:)   = ub  (:,:,:)
-                  vn   (:,:,:)   = vb  (:,:,:)
+                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
+                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
+                  ssh (:,:,Kmm)     = ssh(:,:,Kbb)
+                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
+                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
                ELSE
                   ! if not test case
-                  sshn(:,:) = -ssh_ref
-                  sshb(:,:) = -ssh_ref
+                  ssh(:,:,Kmm) = -ssh_ref
+                  ssh(:,:,Kbb) = -ssh_ref
 
                   DO jj = 1, jpj
@@ -902,7 +902,7 @@
                         IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
 
-                           sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
-                           sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
-                           ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
+                           ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
+                           ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
+                           ssh(ji,jj,Kaa) = rn_wdmin1 - (ht_0(ji,jj) )
                         ENDIF
                      ENDDO
@@ -912,9 +912,9 @@
                ! Adjust vertical metrics for all wad
                DO jk = 1, jpk
-                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) &
+                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) &
                     &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
                     &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
                END DO
-               e3t_b(:,:,:) = e3t_n(:,:,:)
+               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
 
                DO ji = 1, jpi
@@ -928,17 +928,17 @@
             ELSE
                !
-               ! usr_def_istate called here only to get sshb, that is needed to initialize e3t_b and e3t_n
-               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )  
+               ! usr_def_istate called here only to get ssh(:,:,Kbb), that is needed to initialize e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm)
+               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )  
                ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn)
                !
                DO jk=1,jpk
-                  e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &
+                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
                     &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
-                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t_b != 0 on land points
+                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t(:,:,:,Kbb) != 0 on land points
                END DO
-               e3t_n(:,:,:) = e3t_b(:,:,:)
-               sshn(:,:) = sshb(:,:)   ! needed later for gde3w
-!!$                e3t_n(:,:,:)=e3t_0(:,:,:)
-!!$                e3t_b(:,:,:)=e3t_0(:,:,:)
+               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
+               ssh(:,:,Kmm) = ssh(:,:,Kbb)   ! needed later for gde3w
+!!$                e3t(:,:,:,Kmm)=e3t_0(:,:,:)
+!!$                e3t(:,:,:,Kbb)=e3t_0(:,:,:)
                !
             END IF           ! end of ll_wd edits
@@ -958,6 +958,6 @@
          !                                           ! all cases !
          !                                           ! --------- !
-         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios )
          !                                           ! ----------------------- !
          IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
