Index: /NEMO/trunk/src/OCE/LDF/ldfdyn.F90
===================================================================
--- /NEMO/trunk/src/OCE/LDF/ldfdyn.F90	(revision 14200)
+++ /NEMO/trunk/src/OCE/LDF/ldfdyn.F90	(revision 14201)
@@ -99,7 +99,7 @@
       !!    nn_ahm_ijk_t  =  0 => = constant
       !!                  = 10 => = F(z) :     = constant with a reduction of 1/4 with depth 
-      !!                  =-20 => = F(i,j)     = shape read in 'eddy_viscosity.nc' file
+      !!                  =-20 => = F(i,j)     = shape read in 'eddy_viscosity_2D.nc' file
       !!                  = 20    = F(i,j)     = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
-      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_viscosity.nc'  file
+      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_viscosity_3D.nc'  file
       !!                  = 30    = F(i,j,k)   = 2D (case 20) + decrease with depth (case 10)
       !!                  = 31    = F(i,j,k,t) = F(local velocity) (  |u|e  /12   laplacian operator
Index: /NEMO/trunk/src/OCE/LDF/ldftra.F90
===================================================================
--- /NEMO/trunk/src/OCE/LDF/ldftra.F90	(revision 14200)
+++ /NEMO/trunk/src/OCE/LDF/ldftra.F90	(revision 14201)
@@ -118,9 +118,9 @@
       !!                  = 10 => = F(z) : constant with a reduction of 1/4 with depth
       !!                  !
-      !!                  =-20 => = F(i,j)   = shape read in 'eddy_diffusivity.nc' file
+      !!                  =-20 => = F(i,j)   = shape read in 'eddy_diffusivity_2D.nc' file
       !!                  = 20    = F(i,j)   = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
       !!                  = 21    = F(i,j,t) = F(growth rate of baroclinic instability)
       !!                  !
-      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_diffusivity.nc' file
+      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_diffusivity_3D.nc' file
       !!                  = 30    = F(i,j,k)   = 2D (case 20) + decrease with depth (case 10)
       !!                  = 31    = F(i,j,k,t) = F(local velocity) (  1/2  |u|e     laplacian operator
@@ -481,9 +481,9 @@
       !!                  = 10 => = F(z) : constant with a reduction of 1/4 with depth
       !!                  !
-      !!                  =-20 => = F(i,j)   = shape read in 'eddy_diffusivity.nc' file
+      !!                  =-20 => = F(i,j)   = shape read in 'eddy_induced_velocity_2D.nc' file
       !!                  = 20    = F(i,j)   = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
       !!                  = 21    = F(i,j,t) = F(growth rate of baroclinic instability)
       !!                  !
-      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_diffusivity.nc' file
+      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_induced_velocity_3D.nc' file
       !!                  = 30    = F(i,j,k)   = 2D (case 20) + decrease with depth (case 10)
       !!
Index: /NEMO/trunk/src/OCE/step.F90
===================================================================
--- /NEMO/trunk/src/OCE/step.F90	(revision 14200)
+++ /NEMO/trunk/src/OCE/step.F90	(revision 14201)
@@ -50,8 +50,6 @@
    PUBLIC   stp   ! called by nemogcm.F90
 
-   !!----------------------------------------------------------------------
-   !! time level indices
-   !!----------------------------------------------------------------------
-   INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs          !! used by nemo_init
+   !                                          !**  time level indices  **!
+   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !: used by nemo_init
 
    !! * Substitutions
@@ -88,6 +86,4 @@
       !!----------------------------------------------------------------------
       INTEGER ::   ji, jj, jk, jtile   ! dummy loop indice
-!!gm kcall can be removed, I guess
-      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
       !! ---------------------------------------------------------------------
 #if defined key_agrif
@@ -111,7 +107,6 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      IF( l_1st_euler ) THEN
-         ! start or restart with Euler 1st time-step
-         rDt =  rn_Dt
+      IF( l_1st_euler ) THEN     ! start or restart with Euler 1st time-step
+         rDt   = rn_Dt   
          r1_Dt = 1._wp / rDt
       ENDIF
@@ -120,4 +115,5 @@
       ! update I/O and calendar
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
       IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
                              CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including possible AGRIF zoom)
@@ -126,19 +122,19 @@
          IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
       ENDIF
-      IF((kstp == nitrst) .AND. lwxios) THEN
-         CALL iom_swap(      cw_ocerst_cxt          )
-         CALL iom_init_closedef(cw_ocerst_cxt)
-         CALL iom_setkt( kstp - nit000 + 1,      cw_ocerst_cxt          )
+      IF( kstp == nitrst .AND. lwxios ) THEN
+                             CALL iom_swap(                     cw_ocerst_cxt )
+                             CALL iom_init_closedef(            cw_ocerst_cxt )
+                             CALL iom_setkt( kstp - nit000 + 1, cw_ocerst_cxt )
 #if defined key_top
-         CALL iom_swap(      cw_toprst_cxt          )
-         CALL iom_init_closedef(cw_toprst_cxt)
-         CALL iom_setkt( kstp - nit000 + 1,      cw_toprst_cxt          )
+                             CALL iom_swap(                     cw_toprst_cxt )
+                             CALL iom_init_closedef(            cw_toprst_cxt )
+                             CALL iom_setkt( kstp - nit000 + 1, cw_toprst_cxt )
 #endif
       ENDIF
 #if defined key_si3
-      IF(((kstp + nn_fsbc - 1) == nitrst) .AND. lwxios) THEN
-         CALL iom_swap(      cw_icerst_cxt          )
-         CALL iom_init_closedef(cw_icerst_cxt)
-         CALL iom_setkt( kstp - nit000 + 1,      cw_icerst_cxt          )
+      IF( kstp + nn_fsbc - 1 == nitrst .AND. lwxios ) THEN
+                             CALL iom_swap(                     cw_icerst_cxt )
+                             CALL iom_init_closedef(            cw_icerst_cxt )
+                             CALL iom_setkt( kstp - nit000 + 1, cw_icerst_cxt )
       ENDIF
 #endif
@@ -159,6 +155,6 @@
       ! Update stochastic parameters and random T/S fluctuations
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters
-      IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn)  )          ! Random T/S fluctuations
+      IF( ln_sto_eos )   CALL sto_par( kstp )                         ! Stochastic parameters
+      IF( ln_sto_eos )   CALL sto_pts( ts(:,:,:,:,Nnn)  )             ! Random T/S fluctuations
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -200,9 +196,10 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
-      IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors
-                            CALL wzv           ( kstp, Nbb, Nnn, Naa, ww  )    ! now cross-level velocity
-      IF( ln_zad_Aimp )     CALL wAimp         ( kstp,      Nnn           )  ! Adaptive-implicit vertical advection partitioning
-                            CALL eos    ( ts(:,:,:,:,Nnn), rhd, rhop, gdept(:,:,:,Nnn) )  ! now in situ density for hpg computation
+                         CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )   ! after ssh (includes call to div_hor)
+      IF( .NOT.ln_linssh )   &
+                       & CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )   ! after vertical scale factors
+                         CALL wzv           ( kstp, Nbb, Nnn, Naa, ww  )   ! now cross-level velocity
+      IF( ln_zad_Aimp )  CALL wAimp         ( kstp,      Nnn           )   ! Adaptive-implicit vertical advection partitioning
+                         CALL eos    ( ts(:,:,:,:,Nnn), rhd, rhop, gdept(:,:,:,Nnn) )  ! now in situ density for hpg computation
 
 
@@ -231,5 +228,5 @@
                             CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion
       IF( ln_dynspg_ts ) THEN                                                       ! vertical scale factors and vertical velocity need to be updated
-                            CALL wzv        ( kstp, Nbb, Nnn, Naa, ww )             ! now cross-level velocity
+                            CALL wzv        ( kstp, Nbb, Nnn, Naa, ww )             ! Nnn cross-level velocity
          IF( ln_zad_Aimp )  CALL wAimp      ( kstp,      Nnn )                      ! Adaptive-implicit vertical advection partitioning
       ENDIF
@@ -259,5 +256,5 @@
       ! Passive Tracer Model
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL trc_stp       ( kstp, Nbb, Nnn, Nrhs, Naa )  ! time-stepping
+                         CALL trc_stp    ( kstp, Nbb, Nnn, Nrhs, Naa )           ! time-stepping
 #endif
 
@@ -267,8 +264,8 @@
       ! Loop over tile domains
       DO jtile = 1, nijtile
-         IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
+         IF( ln_tile    )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
 
          DO_3D( 0, 0, 0, 0, 1, jpk )
-            ts(ji,jj,jk,:,Nrhs) = 0._wp                                         ! set tracer trends to zero
+            ts(ji,jj,jk,:,Nrhs) = 0._wp                                   ! set tracer trends to zero
          END_3D
 
@@ -285,19 +282,19 @@
 
 #if defined key_agrif
-         IF(.NOT. Agrif_Root()) THEN
-            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )
-            CALL Agrif_Sponge_tra        ! tracers sponge
-         ENDIF
+      IF(.NOT. Agrif_Root()) THEN
+         IF( ln_tile    )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )
+                            CALL Agrif_Sponge_tra        ! tracers sponge
+      ENDIF
 #endif
 
       ! TEMP: [tiling] Separate loop over tile domains (due to tra_adv workarounds for tiling)
       DO jtile = 1, nijtile
-         IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
+         IF( ln_tile    )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
 
                             CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection	==> RHS
          IF( ln_zdfmfc  )   CALL tra_mfc    ( kstp, Nbb,      ts, Nrhs )  ! Mass Flux Convection
          IF( ln_zdfosm  )   CALL tra_osm    ( kstp,      Nnn, ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
-         IF( lrst_oce .AND. ln_zdfosm ) &
-              &             CALL osm_rst    ( kstp,      Nnn, 'WRITE'  )  ! write OSMOSIS outputs + ww (so must do here) to restarts
+         IF( lrst_oce .AND. ln_zdfosm )   &
+            &               CALL osm_rst    ( kstp,      Nnn, 'WRITE'  )  ! write OSMOSIS outputs + ww (so must do here) to restarts
                             CALL tra_ldf    ( kstp, Nbb, Nnn, ts, Nrhs )  ! lateral mixing
 
@@ -363,10 +360,10 @@
       ! AGRIF update
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( Agrif_NbStepint() == 0 .AND. nstop == 0 ) THEN
-                         CALL Agrif_update_all( )                  ! Update all components
-      ENDIF
-
-#endif
-      IF( ln_diaobs .AND. nstop == 0 )  CALL dia_obs( kstp, Nnn )  ! obs-minus-model (assimilation) diags (after dynamics update)
+      IF( Agrif_NbStepint() == 0 .AND. nstop == 0 )   &
+         &               CALL Agrif_update_all( )                  ! Update all components
+
+#endif
+      IF( ln_diaobs .AND. nstop == 0 )   &
+         &               CALL dia_obs( kstp, Nnn )  ! obs-minus-model (assimilation) diags (after dynamics update)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -375,5 +372,5 @@
       IF( kstp == nit000 ) THEN                          ! 1st time step only
                                         CALL iom_close( numror )   ! close input  ocean restart file
-         IF( lrxios )                   CALL iom_context_finalize(      cr_ocerst_cxt         )
+         IF( lrxios )                   CALL iom_context_finalize( cr_ocerst_cxt )
          IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
          IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
@@ -396,5 +393,5 @@
       !
       IF( l_1st_euler ) THEN         ! recover Leap-frog timestep
-         rDt = 2._wp * rn_Dt
+         rDt   = 2._wp * rn_Dt
          r1_Dt = 1._wp / rDt
          l_1st_euler = .FALSE.
Index: /NEMO/trunk/src/OCE/stpmlf.F90
===================================================================
--- /NEMO/trunk/src/OCE/stpmlf.F90	(revision 14200)
+++ /NEMO/trunk/src/OCE/stpmlf.F90	(revision 14201)
@@ -35,5 +35,4 @@
    !!            4.x  !  2020-08  (S. Techene, G. Madec)  quasi eulerian coordinate time stepping
    !!----------------------------------------------------------------------
-
 #if defined key_qco   ||   defined key_linssh
    !!----------------------------------------------------------------------
@@ -42,5 +41,5 @@
    !!   'key_linssh                       Fixed in time vertical coordinate
    !!----------------------------------------------------------------------
-
+   !!
    !!----------------------------------------------------------------------
    !!   stp_MLF       : NEMO modified Leap Frog time-stepping with qco or linssh
@@ -49,6 +48,6 @@
    !
    USE domqco         ! quasi-eulerian coordinate
-   USE traatf_qco     ! time filtering                   (tra_atf_qco routine)
-   USE dynatf_qco     ! time filtering                   (dyn_atf_qco routine)
+   USE traatf_qco     ! time filtering                 (tra_atf_qco routine)
+   USE dynatf_qco     ! time filtering                 (dyn_atf_qco routine)
    
    USE bdydyn         ! ocean open boundary conditions (define bdy_dyn)
@@ -60,5 +59,5 @@
    IMPLICIT NONE
    PRIVATE
-   
+
    PUBLIC   stp_MLF   ! called by nemogcm.F90
 
@@ -67,4 +66,5 @@
 
    !! * Substitutions
+#  include "do_loop_substitute.h90"
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
@@ -98,10 +98,10 @@
       !!              -8- Outputs and diagnostics
       !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj, jk   ! dummy loop indice
-      INTEGER ::   indic        ! error indicator if < 0
+      INTEGER ::   ji, jj, jk, jtile   ! dummy loop indice
       REAL(wp),              DIMENSION(jpi,jpj,jpk) ::   zgdept
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zssh_f
       !! ---------------------------------------------------------------------
 #if defined key_agrif
+      IF( nstop > 0 ) RETURN   ! avoid to go further if an error was detected during previous time step (child grid)
       kstp = nit000 + Agrif_Nb_Step()
       Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
@@ -130,12 +130,28 @@
       ! update I/O and calendar
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                             indic = 0                ! reset to no error condition
-
+      !
       IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
-                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including passible AGRIF zoom)
+                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including possible AGRIF zoom)
          IF( lk_diamlr   )   CALL dia_mlr_iom_init    ! with additional setup for multiple-linear-regression analysis
                              CALL iom_init_closedef
          IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
       ENDIF
+      IF( kstp == nitrst .AND. lwxios ) THEN
+                             CALL iom_swap(                     cw_ocerst_cxt )
+                             CALL iom_init_closedef(            cw_ocerst_cxt )
+                             CALL iom_setkt( kstp - nit000 + 1, cw_ocerst_cxt )
+#if defined key_top
+                             CALL iom_swap(                     cw_toprst_cxt )
+                             CALL iom_init_closedef(            cw_toprst_cxt )
+                             CALL iom_setkt( kstp - nit000 + 1, cw_toprst_cxt )
+#endif
+      ENDIF
+#if defined key_si3
+      IF( kstp + nn_fsbc - 1 == nitrst .AND. lwxios ) THEN
+                             CALL iom_swap(                     cw_icerst_cxt )
+                             CALL iom_init_closedef(            cw_icerst_cxt )
+                             CALL iom_setkt( kstp - nit000 + 1, cw_icerst_cxt )
+      ENDIF
+#endif
       IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
                              CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
@@ -154,6 +170,6 @@
       ! Update stochastic parameters and random T/S fluctuations
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters
-      IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn)  )          ! Random T/S fluctuations
+      IF( ln_sto_eos )   CALL sto_par( kstp )                         ! Stochastic parameters
+      IF( ln_sto_eos )   CALL sto_pts( ts(:,:,:,:,Nnn)  )             ! Random T/S fluctuations
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -197,12 +213,13 @@
          zgdept(:,:,jk) = gdept(:,:,jk,Nnn)
       END DO
-                            CALL ssh_nxt    ( kstp, Nbb, Nnn, ssh,  Naa )   ! after ssh (includes call to div_hor)
-      IF( .NOT.lk_linssh )  THEN
-                             CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa)           )   ! "after" ssh/h_0 ratio at t,u,v pts
-         IF( ln_dynspg_exp ) CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn), r3f(:,:) )   ! spg_exp : needed only for "now" ssh/h_0 ratio at f point
-      ENDIF
-                            CALL wzv        ( kstp, Nbb, Nnn, Naa, ww  )    ! Nnn cross-level velocity
-      IF( ln_zad_Aimp )     CALL wAimp      ( kstp,      Nnn           )    ! Adaptive-implicit vertical advection partitioning
-                            CALL eos        ( ts(:,:,:,:,Nnn), rhd, rhop, zgdept ) ! now in situ density for hpg computation
+                         CALL ssh_nxt    ( kstp, Nbb, Nnn, ssh,  Naa )   ! after ssh (includes call to div_hor)
+      IF( .NOT.lk_linssh ) THEN
+                         CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa)           )   ! "after" ssh/h_0 ratio at t,u,v pts
+         IF( ln_dynspg_exp )   &
+            &            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn), r3f(:,:) )   ! spg_exp : needed only for "now" ssh/h_0 ratio at f point
+      ENDIF
+                         CALL wzv        ( kstp, Nbb, Nnn, Naa, ww  )    ! Nnn cross-level velocity
+      IF( ln_zad_Aimp )  CALL wAimp      ( kstp,      Nnn           )    ! Adaptive-implicit vertical advection partitioning
+                         CALL eos        ( ts(:,:,:,:,Nnn), rhd, rhop, zgdept ) ! now in situ density for hpg computation
 
 
@@ -223,6 +240,6 @@
                          CALL dyn_hpg( kstp,      Nnn      , uu, vv, Nrhs )  ! horizontal gradient of Hydrostatic pressure
                          CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient
+
                                                       ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well
-
       IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated
                             CALL div_hor    ( kstp, Nbb, Nnn )                ! Horizontal divergence  (2nd call in time-split case)
@@ -270,28 +287,46 @@
       ! Active tracers
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         ts(:,:,:,:,Nrhs) = 0._wp         ! set tracer trends to zero
-
-      IF(  lk_asminc .AND. ln_asmiau .AND. &
-         & ln_trainc )   CALL tra_asm_inc( kstp, Nbb, Nnn, ts, Nrhs )  ! apply tracer assimilation increment
-                         CALL tra_sbc    ( kstp,      Nnn, ts, Nrhs )  ! surface boundary condition
-      IF( ln_traqsr  )   CALL tra_qsr    ( kstp,      Nnn, ts, Nrhs )  ! penetrative solar radiation qsr
-      IF( ln_isf     )   CALL tra_isf    ( kstp,      Nnn, ts, Nrhs )  ! ice shelf heat flux
-      IF( ln_trabbc  )   CALL tra_bbc    ( kstp,      Nnn, ts, Nrhs )  ! bottom heat flux
-      IF( ln_trabbl  )   CALL tra_bbl    ( kstp, Nbb, Nnn, ts, Nrhs )  ! advective (and/or diffusive) bottom boundary layer scheme
-      IF( ln_tradmp  )   CALL tra_dmp    ( kstp, Nbb, Nnn, ts, Nrhs )  ! internal damping trends
-      IF( ln_bdy     )   CALL bdy_tra_dmp( kstp, Nbb,      ts, Nrhs )  ! bdy damping trends
-#if defined key_agrif
-      IF(.NOT. Agrif_Root())  &
-               &         CALL Agrif_Sponge_tra        ! tracers sponge
-#endif
-                         CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection	==> RHS
-      IF( ln_zdfosm  )   CALL tra_osm    ( kstp,      Nnn, ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
-      IF( lrst_oce .AND. ln_zdfosm ) &
-           &             CALL osm_rst    ( kstp,      Nnn, 'WRITE'  )  ! write OSMOSIS outputs + ww (so must do here) to restarts
-                         CALL tra_ldf    ( kstp, Nbb, Nnn, ts, Nrhs )  ! lateral mixing
-
-                         CALL tra_zdf    ( kstp, Nbb, Nnn, Nrhs, ts, Naa  )  ! vertical mixing and after tracer fields
-      IF( ln_zdfnpc  )   CALL tra_npc    ( kstp,      Nnn, Nrhs, ts, Naa  )  ! update after fields by non-penetrative convection
-
+      ! Loop over tile domains
+      DO jtile = 1, nijtile
+         IF( ln_tile    )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
+
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            ts(ji,jj,jk,:,Nrhs) = 0._wp                                   ! set tracer trends to zero
+         END_3D
+
+         IF(  lk_asminc .AND. ln_asmiau .AND. &
+            & ln_trainc )   CALL tra_asm_inc( kstp, Nbb, Nnn, ts, Nrhs )  ! apply tracer assimilation increment
+                            CALL tra_sbc    ( kstp,      Nnn, ts, Nrhs )  ! surface boundary condition
+         IF( ln_traqsr  )   CALL tra_qsr    ( kstp,      Nnn, ts, Nrhs )  ! penetrative solar radiation qsr
+         IF( ln_isf     )   CALL tra_isf    ( kstp,      Nnn, ts, Nrhs )  ! ice shelf heat flux
+         IF( ln_trabbc  )   CALL tra_bbc    ( kstp,      Nnn, ts, Nrhs )  ! bottom heat flux
+         IF( ln_trabbl  )   CALL tra_bbl    ( kstp, Nbb, Nnn, ts, Nrhs )  ! advective (and/or diffusive) bottom boundary layer scheme
+         IF( ln_tradmp  )   CALL tra_dmp    ( kstp, Nbb, Nnn, ts, Nrhs )  ! internal damping trends
+         IF( ln_bdy     )   CALL bdy_tra_dmp( kstp, Nbb,      ts, Nrhs )  ! bdy damping trends
+      END DO
+
+#if defined key_agrif
+      IF(.NOT. Agrif_Root()) THEN
+         IF( ln_tile    )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )
+                            CALL Agrif_Sponge_tra        ! tracers sponge
+      ENDIF
+#endif
+
+      ! TEMP: [tiling] Separate loop over tile domains (due to tra_adv workarounds for tiling)
+      DO jtile = 1, nijtile
+         IF( ln_tile    )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
+
+                            CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection	==> RHS
+         IF( ln_zdfmfc  )   CALL tra_mfc    ( kstp, Nbb,      ts, Nrhs )  ! Mass Flux Convection
+         IF( ln_zdfosm  )   CALL tra_osm    ( kstp,      Nnn, ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
+         IF( lrst_oce .AND. ln_zdfosm )   &
+            &               CALL osm_rst    ( kstp,      Nnn, 'WRITE'  )  ! write OSMOSIS outputs + ww (so must do here) to restarts
+                            CALL tra_ldf    ( kstp, Nbb, Nnn, ts, Nrhs )  ! lateral mixing
+
+                            CALL tra_zdf    ( kstp, Nbb, Nnn, Nrhs, ts, Naa  )  ! vertical mixing and after tracer fields
+         IF( ln_zdfnpc  )   CALL tra_npc    ( kstp,      Nnn, Nrhs, ts, Naa  )  ! update after fields by non-penetrative convection
+      END DO
+
+      IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) ! Revert to tile over full domain
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Set boundary conditions, time filter and swap time levels
@@ -319,5 +354,4 @@
                          r3v(:,:,Nnn) = r3v_f(:,:)
       ENDIF
-
       !
       ! Swap time levels
@@ -340,15 +374,10 @@
 #if defined key_agrif
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! AGRIF
+      ! AGRIF recursive integration
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                          Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs      ! agrif_oce module copies of time level indices
                          CALL Agrif_Integrate_ChildGrids( stp_MLF )       ! allows to finish all the Child Grids before updating
 
-                         IF( Agrif_NbStepint() == 0 ) THEN
-                            CALL Agrif_update_all( )                  ! Update all components
-                         ENDIF
-#endif
-      IF( ln_diaobs  )   CALL dia_obs      ( kstp, Nnn )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
-
+#endif
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Control
@@ -356,6 +385,22 @@
                          CALL stp_ctl      ( kstp, Nnn )
 
+#if defined key_agrif
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! AGRIF update
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( Agrif_NbStepint() == 0 .AND. nstop == 0 )   &
+         &               CALL Agrif_update_all( )                  ! Update all components
+      ENDIF
+
+#endif
+      IF( ln_diaobs .AND. nstop == 0 )   &
+         &               CALL dia_obs( kstp, Nnn )  ! obs-minus-model (assimilation) diags (after dynamics update)
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! File manipulation at the end of the first time step
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( kstp == nit000 ) THEN                          ! 1st time step only
                                         CALL iom_close( numror )   ! close input  ocean restart file
+         IF( lrxios )                   CALL iom_context_finalize( cr_ocerst_cxt )
          IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
          IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
@@ -365,9 +410,11 @@
       ! Coupled mode
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-!!gm why lk_oasis and not lk_cpl ????
-      IF( lk_oasis   )   CALL sbc_cpl_snd( kstp, Nbb, Nnn )        ! coupled mode : field exchanges
+      IF( lk_oasis .AND. nstop == 0 )   CALL sbc_cpl_snd( kstp, Nbb, Nnn )     ! coupled mode : field exchanges
       !
 #if defined key_iomput
-      IF( kstp == nitend .OR. indic < 0 ) THEN
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Finalize contextes if end of simulation or error detected
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( kstp == nitend .OR. nstop > 0 ) THEN
                       CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
          IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
@@ -376,7 +423,7 @@
       !
       IF( l_1st_euler ) THEN         ! recover Leap-frog timestep
-         rDt   = 2._wp * rn_Dt   
+         rDt   = 2._wp * rn_Dt
          r1_Dt = 1._wp / rDt
-         l_1st_euler = .FALSE.      
+         l_1st_euler = .FALSE.
       ENDIF
       !
