From 78481193dba1a9aeeb25b2a7b0a8f154b97b6ba6 Mon Sep 17 00:00:00 2001
From: Sebastien Masson <sebastien.masson@locean.ipsl.fr>
Date: Tue, 28 Jun 2022 09:06:17 +0000
Subject: [PATCH] dev_2021_rk3_halo into 68-summer-body-2022

---
 cfgs/SHARED/field_def_nemo-oce.xml | 112 +++++++-------
 cfgs/SHARED/grid_def_nemo.xml      |  33 ++--
 src/OCE/DIA/diadetide.F90          |  21 ++-
 src/OCE/DIA/diahsb.F90             | 136 ++++++++++++-----
 src/OCE/DIA/diahth.F90             |  78 ++++++----
 src/OCE/DIA/diamlr.F90             |   9 +-
 src/OCE/DIA/diaptr.F90             | 237 +++++++++++++++--------------
 src/OCE/DIA/diawri.F90             |   4 +-
 src/OCE/IOM/iom.F90                |  36 +++--
 src/OCE/IOM/iom_nf90.F90           |  37 ++---
 src/OCE/IOM/prtctl.F90             |  44 ++++--
 src/OCE/TRA/eosbn2.F90             |  14 +-
 src/OCE/TRA/traadv.F90             |   2 +-
 src/OCE/TRD/trddyn.F90             |  17 ++-
 src/OCE/TRD/trdglo.F90             | 111 +++++++-------
 src/OCE/TRD/trdken.F90             | 101 ++++++------
 src/OCE/TRD/trdmxl.F90             | 115 +++++---------
 src/OCE/TRD/trdmxl_oce.F90         |  48 +++---
 src/OCE/TRD/trdpen.F90             |  36 +++--
 src/OCE/TRD/trdtra.F90             | 122 +++++++--------
 src/OCE/TRD/trdvor.F90             |  10 +-
 src/OCE/lib_fortran_generic.h90    |   4 +-
 22 files changed, 705 insertions(+), 622 deletions(-)

diff --git a/cfgs/SHARED/field_def_nemo-oce.xml b/cfgs/SHARED/field_def_nemo-oce.xml
index e5bff5cc7..838cd2834 100644
--- a/cfgs/SHARED/field_def_nemo-oce.xml
+++ b/cfgs/SHARED/field_def_nemo-oce.xml
@@ -25,7 +25,7 @@ that are available in the tidal-forcing implementation (see
     -->
 
     <!-- Time -->
-    <field id="diamlr_time" grid_ref="diamlr_grid_T_2D" prec="8" />
+    <field id="diamlr_time" grid_ref="diamlr_grid_T_2D_inner" prec="8" />
 
     <!-- Regressors for tidal harmonic analysis -->
     <field id="diamlr_r001"  field_ref="diamlr_time" expr="sin( __TDE_M2_omega__   * diamlr_time )" enabled=".TRUE."  comment="harmonic:sin:M2"   />
@@ -199,22 +199,22 @@ that are available in the tidal-forcing implementation (see
     <field id="wi_cff"       long_name="Fraction of implicit vertical velocity"                                         unit="#"   grid_ref="grid_T_3D" />
 
     <!-- next variables available with key_diahth -->
-    <field id="mlddzt"       long_name="Thermocline Depth (depth of max dT/dz)"         standard_name="depth_at_maximum_upward_derivative_of_sea_water_potential_temperature"             unit="m"                         />
-    <field id="mldr10_3"     long_name="Mixed Layer Depth (dsigma = 0.03 wrt 10m)"      standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"                         />
-    <field id="mldr0_1"      long_name="Mixed Layer Depth (dsigma = 0.01 wrt sfc)"      standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"                         />
-    <field id="mldr0_3"      long_name="Mixed Layer Depth (dsigma = 0.03 wrt sfc)"      standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"                         />
-    <field id="mld_dt02"     long_name="Mixed Layer Depth (|dT| = 0.2 wrt 10m)"         standard_name="ocean_mixed_layer_thickness_defined_by_temperature"                                unit="m"                         />
-    <field id="topthdep"     long_name="Top of Thermocline Depth (dT = -0.2 wrt 10m)"   standard_name="ocean_mixed_layer_thickness_defined_by_temperature"                                unit="m"                         />
-    <field id="pycndep"      long_name="Pycnocline Depth (dsigma[dT=-0.2] wrt 10m)"     standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"                         />
-    <field id="BLT"          long_name="Barrier Layer Thickness"                                                                                                                          unit="m"   > topthdep - pycndep </field>
-    <field id="tinv"         long_name="Max of vertical invertion of temperature"                                                                                                         unit="degC"                      />
-    <field id="depti"        long_name="Depth of max. vert. inv. of temperature"                                                                                                          unit="m"                         />
-    <field id="20d"          long_name="Depth of 20C isotherm"                          standard_name="depth_of_isosurface_of_sea_water_potential_temperature"                            unit="m"   axis_ref="iax_20C" />
-    <field id="26d"          long_name="Depth of 26C isotherm"                          standard_name="depth_of_isosurface_of_sea_water_potential_temperature"                            unit="m"   axis_ref="iax_26C"   />
-    <field id="28d"          long_name="Depth of 28C isotherm"                          standard_name="depth_of_isosurface_of_sea_water_potential_temperature"                            unit="m"   axis_ref="iax_28C" />
-    <field id="hc300"        long_name="Heat content 0-300m"                            standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"                      />
-    <field id="hc700"        long_name="Heat content 0-700m"                            standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"                      />
-    <field id="hc2000"       long_name="Heat content 0-2000m"                           standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"                      />
+    <field id="mlddzt"       long_name="Thermocline Depth (depth of max dT/dz)"         standard_name="depth_at_maximum_upward_derivative_of_sea_water_potential_temperature"             unit="m"      grid_ref="grid_T_2D_inner" />
+    <field id="mldr10_3"     long_name="Mixed Layer Depth (dsigma = 0.03 wrt 10m)"      standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"      grid_ref="grid_T_2D_inner" />
+    <field id="mldr0_1"      long_name="Mixed Layer Depth (dsigma = 0.01 wrt sfc)"      standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"      grid_ref="grid_T_2D_inner" />
+    <field id="mldr0_3"      long_name="Mixed Layer Depth (dsigma = 0.03 wrt sfc)"      standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"      grid_ref="grid_T_2D_inner" />
+    <field id="mld_dt02"     long_name="Mixed Layer Depth (|dT| = 0.2 wrt 10m)"         standard_name="ocean_mixed_layer_thickness_defined_by_temperature"                                unit="m"      grid_ref="grid_T_2D_inner" />
+    <field id="topthdep"     long_name="Top of Thermocline Depth (dT = -0.2 wrt 10m)"   standard_name="ocean_mixed_layer_thickness_defined_by_temperature"                                unit="m"      grid_ref="grid_T_2D_inner" />
+    <field id="pycndep"      long_name="Pycnocline Depth (dsigma[dT=-0.2] wrt 10m)"     standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"      grid_ref="grid_T_2D_inner" />
+    <field id="BLT"          long_name="Barrier Layer Thickness"                                                                                                                          unit="m"      grid_ref="grid_T_2D_inner" > topthdep - pycndep </field>
+    <field id="tinv"         long_name="Max of vertical invertion of temperature"                                                                                                         unit="degC"   grid_ref="grid_T_2D_inner" />
+    <field id="depti"        long_name="Depth of max. vert. inv. of temperature"                                                                                                          unit="m"      grid_ref="grid_T_2D_inner" />
+    <field id="20d"          long_name="Depth of 20C isotherm"                          standard_name="depth_of_isosurface_of_sea_water_potential_temperature"                            unit="m"      grid_ref="grid_T_2D_inner"  axis_ref="iax_20C" />
+    <field id="26d"          long_name="Depth of 26C isotherm"                          standard_name="depth_of_isosurface_of_sea_water_potential_temperature"                            unit="m"      grid_ref="grid_T_2D_inner"  axis_ref="iax_26C"   />
+    <field id="28d"          long_name="Depth of 28C isotherm"                          standard_name="depth_of_isosurface_of_sea_water_potential_temperature"                            unit="m"      grid_ref="grid_T_2D_inner"  axis_ref="iax_28C" />
+    <field id="hc300"        long_name="Heat content 0-300m"                            standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"   grid_ref="grid_T_2D_inner"  />
+    <field id="hc700"        long_name="Heat content 0-700m"                            standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"   grid_ref="grid_T_2D_inner"  />
+    <field id="hc2000"       long_name="Heat content 0-2000m"                           standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"   grid_ref="grid_T_2D_inner"  />
 
     <!-- variables available with diaar5 -->
     <field id="botpres"      long_name="Sea Water Pressure at Sea Floor"          standard_name="sea_water_pressure_at_sea_floor"                    unit="dbar" />
@@ -1018,7 +1018,7 @@ that are available in the tidal-forcing implementation (see
   </field_group>
 
   <!--  Total trends calculated every time step-->
-  <field_group id="trendT" grid_ref="grid_T_3D">
+  <field_group id="trendT" grid_ref="grid_T_3D_inner">
     <field id="ttrd_tot"      long_name="temperature-trend: total model trend"         unit="degC/s" />
     <field id="strd_tot"      long_name="salinity   -trend: total model trend"         unit="1e-3/s" />
     <!-- Thickness weighted versions: -->
@@ -1034,10 +1034,10 @@ that are available in the tidal-forcing implementation (see
     <field id="ketrd_spg"     long_name="ke-trend: surface     pressure gradient"          unit="W/s^3"                        />
     <field id="ketrd_spgexp"  long_name="ke-trend: surface pressure gradient (explicit)"   unit="W/s^3"                        />
     <field id="ketrd_spgflt"  long_name="ke-trend: surface pressure gradient (filter)"     unit="W/s^3"                        />
-    <field id="ssh_flt"       long_name="filtered contribution to ssh (dynspg_flt)"        unit="m"       grid_ref="grid_T_2D" />
-    <field id="w0"            long_name="surface vertical velocity"                        unit="m/s"     grid_ref="grid_T_2D" />
-    <field id="pw0_exp"       long_name="surface pressure flux due to ssh"                 unit="W/s^2"   grid_ref="grid_T_2D" />
-    <field id="pw0_flt"       long_name="surface pressure flux due to filtered ssh"        unit="W/s^2"   grid_ref="grid_T_2D" />
+    <field id="ssh_flt"       long_name="filtered contribution to ssh (dynspg_flt)"        unit="m"       grid_ref="grid_T_2D_inner" />
+    <field id="w0"            long_name="surface vertical velocity"                        unit="m/s"     grid_ref="grid_T_2D_inner" />
+    <field id="pw0_exp"       long_name="surface pressure flux due to ssh"                 unit="W/s^2"   grid_ref="grid_T_2D_inner" />
+    <field id="pw0_flt"       long_name="surface pressure flux due to filtered ssh"        unit="W/s^2"   grid_ref="grid_T_2D_inner" />
     <field id="ketrd_keg"     long_name="ke-trend: KE gradient         or hor. adv."       unit="W/s^3"                        />
     <field id="ketrd_rvo"     long_name="ke-trend: relative  vorticity or metric term"     unit="W/s^3"                        />
     <field id="ketrd_pvo"     long_name="ke-trend: planetary vorticity"                    unit="W/s^3"                        />
@@ -1045,7 +1045,7 @@ that are available in the tidal-forcing implementation (see
     <field id="ketrd_udx"     long_name="ke-trend: U.dx[U]"                                unit="W/s^3"                        />
     <field id="ketrd_ldf"     long_name="ke-trend: lateral   diffusion"                    unit="W/s^3"                        />
     <field id="ketrd_zdf"     long_name="ke-trend: vertical  diffusion"                    unit="W/s^3"                        />
-    <field id="ketrd_tau"     long_name="ke-trend: wind stress "                           unit="W/s^3"   grid_ref="grid_T_2D" />
+    <field id="ketrd_tau"     long_name="ke-trend: wind stress "                           unit="W/s^3"   grid_ref="grid_T_2D_inner" />
     <field id="ketrd_bfr"     long_name="ke-trend: bottom friction (explicit)"             unit="W/s^3"                        />
     <field id="ketrd_bfri"    long_name="ke-trend: bottom friction (implicit)"             unit="W/s^3"                        />
     <field id="ketrd_atf"     long_name="ke-trend: asselin time filter trend"              unit="W/s^3"                        />
@@ -1062,7 +1062,7 @@ that are available in the tidal-forcing implementation (see
     <field id="petrd_xad"     long_name="pe-trend: i-advection"                unit="W/m^3"                        />
     <field id="petrd_yad"     long_name="pe-trend: j-advection"                unit="W/m^3"                        />
     <field id="petrd_zad"     long_name="pe-trend: k-advection"                unit="W/m^3"                        />
-    <field id="petrd_sad"     long_name="pe-trend: surface adv. (linssh true)" unit="W/m^3"   grid_ref="grid_T_2D" />
+    <field id="petrd_sad"     long_name="pe-trend: surface adv. (linssh true)" unit="W/m^3"   grid_ref="grid_T_2D_inner" />
     <field id="petrd_ldf"     long_name="pe-trend: lateral  diffusion"         unit="W/m^3"                        />
     <field id="petrd_zdf"     long_name="pe-trend: vertical diffusion"         unit="W/m^3"                        />
     <field id="petrd_zdfp"    long_name="pe-trend: pure vert. diffusion"       unit="W/m^3"                        />
@@ -1080,42 +1080,42 @@ that are available in the tidal-forcing implementation (see
 
   <field_group id="trendU" grid_ref="grid_U_3D">
     <!-- variables available with ln_dyn_trd -->
-    <field id="utrd_hpg"       long_name="i-trend: hydrostatic pressure gradient"          unit="m/s^2"                        />
-    <field id="utrd_spg"       long_name="i-trend: surface     pressure gradient"          unit="m/s^2"                        />
-    <field id="utrd_spgexp"    long_name="i-trend: surface pressure gradient (explicit)"   unit="m/s^2"                        />
-    <field id="utrd_spgflt"    long_name="i-trend: surface pressure gradient (filtered)"   unit="m/s^2"                        />
-    <field id="utrd_keg"       long_name="i-trend: KE gradient         or hor. adv."       unit="m/s^2"                        />
-    <field id="utrd_rvo"       long_name="i-trend: relative  vorticity or metric term"     unit="m/s^2"                        />
-    <field id="utrd_pvo"       long_name="i-trend: planetary vorticity"                    unit="m/s^2"                        />
-    <field id="utrd_zad"       long_name="i-trend: vertical  advection"                    unit="m/s^2"                        />
-    <field id="utrd_udx"       long_name="i-trend: U.dx[U]"                                unit="m/s^2"                        />
-    <field id="utrd_ldf"       long_name="i-trend: lateral   diffusion"                    unit="m/s^2"                        />
-    <field id="utrd_zdf"       long_name="i-trend: vertical  diffusion"                    unit="m/s^2"                        />
-    <field id="utrd_tau"       long_name="i-trend: wind stress "                           unit="m/s^2"   grid_ref="grid_U_2D" />
-    <field id="utrd_bfr"       long_name="i-trend: bottom friction (explicit)"             unit="m/s^2"                        />
-    <field id="utrd_bfri"      long_name="i-trend: bottom friction (implicit)"             unit="m/s^2"                        />
-    <field id="utrd_tot"       long_name="i-trend: total momentum trend before atf"        unit="m/s^2"                        />
-    <field id="utrd_atf"       long_name="i-trend: asselin time filter trend"              unit="m/s^2"                        />
+    <field id="utrd_hpg"       long_name="i-trend: hydrostatic pressure gradient"          unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_spg"       long_name="i-trend: surface     pressure gradient"          unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_spgexp"    long_name="i-trend: surface pressure gradient (explicit)"   unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_spgflt"    long_name="i-trend: surface pressure gradient (filtered)"   unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_keg"       long_name="i-trend: KE gradient         or hor. adv."       unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_rvo"       long_name="i-trend: relative  vorticity or metric term"     unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_pvo"       long_name="i-trend: planetary vorticity"                    unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_zad"       long_name="i-trend: vertical  advection"                    unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_udx"       long_name="i-trend: U.dx[U]"                                unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_ldf"       long_name="i-trend: lateral   diffusion"                    unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_zdf"       long_name="i-trend: vertical  diffusion"                    unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_tau"       long_name="i-trend: wind stress "                           unit="m/s^2"   grid_ref="grid_U_2D_inner"  />
+    <field id="utrd_bfr"       long_name="i-trend: bottom friction (explicit)"             unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_bfri"      long_name="i-trend: bottom friction (implicit)"             unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_tot"       long_name="i-trend: total momentum trend before atf"        unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="utrd_atf"       long_name="i-trend: asselin time filter trend"              unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
   </field_group>
 
   <field_group id="trendV" grid_ref="grid_V_3D">
     <!-- variables available with ln_dyn_trd -->
-    <field id="vtrd_hpg"       long_name="j-trend: hydrostatic pressure gradient"          unit="m/s^2"                        />
-    <field id="vtrd_spg"       long_name="j-trend: surface     pressure gradient"          unit="m/s^2"                        />
-    <field id="vtrd_spgexp"    long_name="j-trend: surface pressure gradient (explicit)"   unit="m/s^2"                        />
-    <field id="vtrd_spgflt"    long_name="j-trend: surface pressure gradient (filtered)"   unit="m/s^2"                        />
-    <field id="vtrd_keg"       long_name="j-trend: KE gradient         or hor. adv."       unit="m/s^2"                        />
-    <field id="vtrd_rvo"       long_name="j-trend: relative  vorticity or metric term"     unit="m/s^2"                        />
-    <field id="vtrd_pvo"       long_name="j-trend: planetary vorticity"                    unit="m/s^2"                        />
-    <field id="vtrd_zad"       long_name="j-trend: vertical  advection"                    unit="m/s^2"                        />
-    <field id="vtrd_vdy"       long_name="i-trend: V.dx[V]"                                unit="m/s^2"                        />
-    <field id="vtrd_ldf"       long_name="j-trend: lateral   diffusion"                    unit="m/s^2"                        />
-    <field id="vtrd_zdf"       long_name="j-trend: vertical  diffusion"                    unit="m/s^2"                        />
-    <field id="vtrd_tau"       long_name="j-trend: wind stress "                           unit="m/s^2"   grid_ref="grid_V_2D" />
-    <field id="vtrd_bfr"       long_name="j-trend: bottom friction (explicit)"             unit="m/s^2"                        />
-    <field id="vtrd_bfri"      long_name="j-trend: bottom friction (implicit)"             unit="m/s^2"                        />
-    <field id="vtrd_tot"       long_name="j-trend: total momentum trend before atf"        unit="m/s^2"                        />
-    <field id="vtrd_atf"       long_name="j-trend: asselin time filter trend"              unit="m/s^2"                        />
+    <field id="vtrd_hpg"       long_name="j-trend: hydrostatic pressure gradient"          unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_spg"       long_name="j-trend: surface     pressure gradient"          unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_spgexp"    long_name="j-trend: surface pressure gradient (explicit)"   unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_spgflt"    long_name="j-trend: surface pressure gradient (filtered)"   unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_keg"       long_name="j-trend: KE gradient         or hor. adv."       unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_rvo"       long_name="j-trend: relative  vorticity or metric term"     unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_pvo"       long_name="j-trend: planetary vorticity"                    unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_zad"       long_name="j-trend: vertical  advection"                    unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_vdy"       long_name="i-trend: V.dx[V]"                                unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_ldf"       long_name="j-trend: lateral   diffusion"                    unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_zdf"       long_name="j-trend: vertical  diffusion"                    unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_tau"       long_name="j-trend: wind stress "                           unit="m/s^2"   grid_ref="grid_V_2D_inner"  />
+    <field id="vtrd_bfr"       long_name="j-trend: bottom friction (explicit)"             unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_bfri"      long_name="j-trend: bottom friction (implicit)"             unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_tot"       long_name="j-trend: total momentum trend before atf"        unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
+    <field id="vtrd_atf"       long_name="j-trend: asselin time filter trend"              unit="m/s^2"   grid_ref="grid_T_3D_inner"  />
   </field_group>
 
   <!-- shared variables available with TOP interface -->
diff --git a/cfgs/SHARED/grid_def_nemo.xml b/cfgs/SHARED/grid_def_nemo.xml
index d66b314f3..374f67c65 100644
--- a/cfgs/SHARED/grid_def_nemo.xml
+++ b/cfgs/SHARED/grid_def_nemo.xml
@@ -336,41 +336,44 @@
     <scalar />
   </grid>
   <grid id="diamlr_grid_T_2D" >
-    <domain domain_ref="grid_T" />
+    <domain domain_ref="grid_T_inner" />
     <scalar />
   </grid>
+  <grid id="diamlr_grid_T_2D_inner" >
+    <domain domain_ref="grid_T_inner" />
+  </grid>
   <grid id="diamlr_grid_U_2D" >
-    <domain domain_ref="grid_U" />
+    <domain domain_ref="grid_U_inner" />
     <scalar />
   </grid>
   <grid id="diamlr_grid_V_2D" >
-    <domain domain_ref="grid_V" />
+    <domain domain_ref="grid_V_inner" />
     <scalar />
   </grid>
   <grid id="diamlr_grid_W_2D" >
-    <domain domain_ref="grid_W" />
+    <domain domain_ref="grid_W_inner" />
     <scalar />
   </grid>
   <grid id="diamlr_grid_2D_to_grid_T_3D" >
-    <domain domain_ref="grid_T" />
+    <domain domain_ref="grid_T_inner" />
     <axis axis_ref="deptht">
       <duplicate_scalar />
     </axis>
   </grid>
   <grid id="diamlr_grid_2D_to_grid_U_3D" >
-    <domain domain_ref="grid_U" />
+    <domain domain_ref="grid_U_inner" />
     <axis axis_ref="depthu">
       <duplicate_scalar />
     </axis>
   </grid>
   <grid id="diamlr_grid_2D_to_grid_V_3D" >
-    <domain domain_ref="grid_V" />
+    <domain domain_ref="grid_V_inner" />
     <axis axis_ref="depthv">
       <duplicate_scalar />
     </axis>
   </grid>
   <grid id="diamlr_grid_2D_to_grid_W_3D" >
-    <domain domain_ref="grid_W" />
+    <domain domain_ref="grid_W_inner" />
     <axis axis_ref="depthw">
       <duplicate_scalar />
     </axis>
@@ -383,37 +386,37 @@
   </grid>
   <!-- grid definitions for the computation of daily detided model diagnostics (diadetide) -->
   <grid id="diadetide_grid_T_2D" >
-    <domain domain_ref="grid_T" />
+    <domain domain_ref="grid_T_inner" />
     <scalar />
   </grid>
   <grid id="diadetide_grid_U_2D" >
-    <domain domain_ref="grid_U" />
+    <domain domain_ref="grid_U_inner" />
     <scalar />
   </grid>
   <grid id="diadetide_grid_V_2D" >
-    <domain domain_ref="grid_V" />
+    <domain domain_ref="grid_V_inner" />
     <scalar />
   </grid>
   <grid id="diadetide_grid_2D_to_grid_T_3D" >
-    <domain domain_ref="grid_T" />
+    <domain domain_ref="grid_T_inner" />
     <axis axis_ref="deptht">
       <duplicate_scalar />
     </axis>
   </grid>
   <grid id="diadetide_grid_2D_to_grid_U_3D" >
-    <domain domain_ref="grid_U" />
+    <domain domain_ref="grid_U_inner" />
     <axis axis_ref="depthu">
       <duplicate_scalar />
     </axis>
   </grid>
   <grid id="diadetide_grid_2D_to_grid_V_3D" >
-    <domain domain_ref="grid_V" />
+    <domain domain_ref="grid_V_inner" />
     <axis axis_ref="depthv">
       <duplicate_scalar />
     </axis>
   </grid>
   <grid id="diadetide_grid_2D_to_grid_W_3D" >
-    <domain domain_ref="grid_W" />
+    <domain domain_ref="grid_W_inner" />
     <axis axis_ref="depthw">
       <duplicate_scalar />
     </axis>
diff --git a/src/OCE/DIA/diadetide.F90 b/src/OCE/DIA/diadetide.F90
index 9d6760039..43dcc23de 100644
--- a/src/OCE/DIA/diadetide.F90
+++ b/src/OCE/DIA/diadetide.F90
@@ -5,11 +5,11 @@ MODULE diadetide
    !!======================================================================
    !! History :       !  2019  (S. Mueller)
    !!----------------------------------------------------------------------
-   USE par_oce        , ONLY :   wp, jpi, jpj
-   USE in_out_manager , ONLY :   lwp, numout
-   USE iom            , ONLY :   iom_put
-   USE dom_oce        , ONLY :   rn_Dt, nsec_day
-   USE phycst         , ONLY :   rpi
+   USE par_oce        
+   USE in_out_manager 
+   USE iom            
+   USE dom_oce        
+   USE phycst        
    USE tide_mod
 #if defined key_xios
    USE xios
@@ -24,6 +24,8 @@ MODULE diadetide
 
    PUBLIC ::   dia_detide_init, dia_detide
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2019)
    !! $Id$
@@ -90,9 +92,9 @@ CONTAINS
       !!----------------------------------------------------------------------
 
       INTEGER, INTENT(in)          ::   kt
-      REAL(wp), DIMENSION(jpi,jpj) ::   zwght_2D
+      REAL(wp), DIMENSION(A2D(0))  ::   zwght_2D
       REAL(wp)                     ::   zwght, ztmp
-      INTEGER                      ::   jn
+      INTEGER                      ::   ji, jj, jn
 
       ! Compute detiding weight at the current time-step; the daily total weight
       ! is one, and the daily summation of a diagnosed field multiplied by this
@@ -104,7 +106,10 @@ CONTAINS
             zwght = zwght + 1.0_wp / REAL( ndiadetide, KIND=wp )
          END IF
       END DO
-      zwght_2D(:,:) = zwght
+ 
+      DO_2D( 0, 0, 0, 0 )
+         zwght_2D(ji,jj) = zwght
+      END_2D
       CALL iom_put( "diadetide_weight", zwght_2D)
 
    END SUBROUTINE dia_detide
diff --git a/src/OCE/DIA/diahsb.F90 b/src/OCE/DIA/diahsb.F90
index f9307c0af..1039d5c3c 100644
--- a/src/OCE/DIA/diahsb.F90
+++ b/src/OCE/DIA/diahsb.F90
@@ -50,6 +50,7 @@ MODULE diahsb
    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   tmask_ini
 
    !! * Substitutions
+#  include "do_loop_substitute.h90"
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -82,30 +83,61 @@ CONTAINS
       REAL(wp)   ::   z_frc_trd_v                 !    -     -
       REAL(wp)   ::   z_wn_trd_t , z_wn_trd_s     !    -     -
       REAL(wp)   ::   z_ssh_hc , z_ssh_sc         !    -     -
-      REAL(wp), DIMENSION(jpi,jpj,13)      ::   ztmp
-      REAL(wp), DIMENSION(jpi,jpj,jpkm1,4) ::   ztmpk
+      REAL(wp), DIMENSION(A2D(0),13)      ::   ztmp
+      REAL(wp), DIMENSION(A2D(0),jpkm1,4) ::   ztmpk
       REAL(wp), DIMENSION(17)              ::   zbg          
       !!---------------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('dia_hsb')
       !
-      ztmp (:,:,:)   = 0._wp ! should be better coded
-      ztmpk(:,:,:,:) = 0._wp ! should be better coded
-      !
-      ts(:,:,:,1,Kmm) = ts(:,:,:,1,Kmm) * tmask(:,:,:) ; ts(:,:,:,1,Kbb) = ts(:,:,:,1,Kbb) * tmask(:,:,:) ;
-      ts(:,:,:,2,Kmm) = ts(:,:,:,2,Kmm) * tmask(:,:,:) ; ts(:,:,:,2,Kbb) = ts(:,:,:,2,Kbb) * tmask(:,:,:) ;
+      DO_2D( 0, 0, 0, 0 )
+         ztmp (ji,jj,:)   = 0._wp ! should be better coded
+         ztmpk(ji,jj,:,:) = 0._wp ! should be better coded
+         !
+         ts(ji,jj,:,1,Kmm) = ts(ji,jj,:,1,Kmm) * tmask(ji,jj,:)
+         ts(ji,jj,:,1,Kbb) = ts(ji,jj,:,1,Kbb) * tmask(ji,jj,:)
+         !
+         ts(ji,jj,:,2,Kmm) = ts(ji,jj,:,2,Kmm) * tmask(ji,jj,:)
+         ts(ji,jj,:,2,Kbb) = ts(ji,jj,:,2,Kbb) * tmask(ji,jj,:)
+      END_2D
       !
       ! ------------------------- !
       ! 1 - Trends due to forcing !
       ! ------------------------- !
       ! prepare trends
-      ztmp(:,:,1)  = - r1_rho0 * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) ) * surf(:,:)    ! volume
-      ztmp(:,:,2)  =   sbc_tsc(:,:,jp_tem) * surf(:,:)                                                      ! heat
-      ztmp(:,:,3)  =   sbc_tsc(:,:,jp_sal) * surf(:,:)                                                      ! salt
-      IF( ln_rnf     )    ztmp(:,:,4) =   rnf_tsc(:,:,jp_tem) * surf(:,:)                                   ! runoff temp
-      IF( ln_rnf_sal )    ztmp(:,:,5) =   rnf_tsc(:,:,jp_sal) * surf(:,:)                                   ! runoff salt
-      IF( ln_isf     )    ztmp(:,:,6) = ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) ! isf temp
-      IF( ln_traqsr  )    ztmp(:,:,7) =   r1_rho0_rcp * qsr(:,:) * surf(:,:)                                ! penetrative solar radiation
-      IF( ln_trabbc  )    ztmp(:,:,8) =   qgh_trd0(:,:) * surf(:,:)                                         ! geothermal heat
+      DO_2D( 0, 0, 0, 0 )
+         ztmp(ji,jj,1)  = - r1_rho0 * (   emp(ji,jj)        &                         ! volume
+             &                          - rnf(ji,jj)        &
+             &                          - fwfisf_cav(ji,jj) &
+             &                          - fwfisf_par(ji,jj) ) * surf(ji,jj)
+         ztmp(ji,jj,2)  =   sbc_tsc(ji,jj,jp_tem) * surf(ji,jj)                       ! heat
+         ztmp(ji,jj,3)  =   sbc_tsc(ji,jj,jp_sal) * surf(ji,jj)                       ! salt
+      END_2D
+      IF( ln_rnf     ) THEN
+         DO_2D( 0, 0, 0, 0 )
+            ztmp(ji,jj,4) = rnf_tsc(ji,jj,jp_tem) * surf(ji,jj)    ! runoff temp
+         END_2D
+      END IF
+      IF( ln_rnf_sal ) THEN
+         DO_2D( 0, 0, 0, 0 )
+            ztmp(ji,jj,5) = rnf_tsc(ji,jj,jp_sal) * surf(ji,jj)    ! runoff salt
+         END_2D
+      END IF
+      IF( ln_isf     ) THEN
+         DO_2D( 0, 0, 0, 0 )
+            ztmp(ji,jj,6) = (   risf_cav_tsc(ji,jj,jp_tem) &
+                &             + risf_par_tsc(ji,jj,jp_tem) ) * surf(ji,jj) ! isf temp
+         END_2D
+      END IF
+      IF( ln_traqsr  ) THEN
+         DO_2D( 0, 0, 0, 0 )
+            ztmp(ji,jj,7) = r1_rho0_rcp * qsr(ji,jj) * surf(ji,jj) ! penetrative solar radiation
+         END_2D
+      END IF
+      IF( ln_trabbc  ) THEN
+         DO_2D( 0, 0, 0, 0 )
+            ztmp(ji,jj,8) = qgh_trd0(ji,jj) * surf(ji,jj)             ! geothermal heat
+         END_2D
+      END IF
       !
       IF( ln_linssh ) THEN   ! Advection flux through fixed surface (z=0)
          IF( ln_isfcav ) THEN
@@ -116,8 +148,10 @@ CONTAINS
                END DO
             END DO
          ELSE
-            ztmp(:,:,9 ) = - surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_tem,Kbb)
-            ztmp(:,:,10) = - surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_sal,Kbb)
+            DO_2D( 0, 0, 0, 0 )
+               ztmp(ji,jj,9 ) = - surf(ji,jj) * ww(ji,jj,1) * ts(ji,jj,1,jp_tem,Kbb)
+               ztmp(ji,jj,10) = - surf(ji,jj) * ww(ji,jj,1) * ts(ji,jj,1,jp_sal,Kbb)
+            END_2D
          END IF
       ENDIF
       
@@ -152,7 +186,9 @@ CONTAINS
       ! glob_sum is needed because you keep only the interior domain to compute the sum (iscpl)
       !
       !                    ! volume variation (calculated with ssh)
-      ztmp(:,:,11) = surf(:,:)*ssh(:,:,Kmm) - surf_ini(:,:)*ssh_ini(:,:)
+      DO_2D( 0, 0, 0, 0 )
+         ztmp(ji,jj,11) = surf(ji,jj)*ssh(ji,jj,Kmm) - surf_ini(ji,jj)*ssh_ini(ji,jj)
+      END_2D
 
       !                    ! heat & salt content variation (associated with ssh)
       IF( ln_linssh ) THEN       ! linear free surface case
@@ -164,8 +200,10 @@ CONTAINS
                END DO
             END DO
          ELSE                          ! no under ice-shelf seas
-            ztmp(:,:,12) = surf(:,:) * ( ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) - ssh_hc_loc_ini(:,:) )
-            ztmp(:,:,13) = surf(:,:) * ( ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) - ssh_sc_loc_ini(:,:) )
+            DO_2D( 0, 0, 0, 0 )
+               ztmp(ji,jj,12) = surf(ji,jj) * ( ts(ji,jj,1,jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) )
+               ztmp(ji,jj,13) = surf(ji,jj) * ( ts(ji,jj,1,jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) )
+            END_2D
          END IF
       ENDIF
 
@@ -185,19 +223,27 @@ CONTAINS
       ! glob_sum is needed because you keep only the interior domain to compute the sum (iscpl)
       !
       DO jk = 1, jpkm1           ! volume
-         ztmpk(:,:,jk,1) =   surf    (:,:) * e3t(:,:,jk,Kmm)*tmask(:,:,jk)   &
-            &              - surf_ini(:,:) * e3t_ini(:,:,jk    )*tmask_ini(:,:,jk)
+         DO_2D( 0, 0, 0, 0 )
+            ztmpk(ji,jj,jk,1) =   surf    (ji,jj) * e3t(ji,jj,jk,Kmm)*tmask(ji,jj,jk)   &
+               &                - surf_ini(ji,jj) * e3t_ini(ji,jj,jk    )*tmask_ini(ji,jj,jk)
+         END_2D
       END DO
       DO jk = 1, jpkm1           ! heat
-         ztmpk(:,:,jk,2) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm)   &
-            &              - surf_ini(:,:) *         hc_loc_ini(:,:,jk) )
+         DO_2D( 0, 0, 0, 0 )
+            ztmpk(ji,jj,jk,2) = ( surf    (ji,jj) * e3t(ji,jj,jk,Kmm)*ts(ji,jj,jk,jp_tem,Kmm)   &
+               &                - surf_ini(ji,jj) *         hc_loc_ini(ji,jj,jk) )
+         END_2D
       END DO
       DO jk = 1, jpkm1           ! salt
-         ztmpk(:,:,jk,3) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm)   &
-            &              - surf_ini(:,:) *         sc_loc_ini(:,:,jk) )
+         DO_2D( 0, 0, 0, 0 )
+            ztmpk(ji,jj,jk,3) = ( surf    (ji,jj) * e3t(ji,jj,jk,Kmm)*ts(ji,jj,jk,jp_sal,Kmm)   &
+               &                - surf_ini(ji,jj) *         sc_loc_ini(ji,jj,jk) )
+         END_2D
       END DO
       DO jk = 1, jpkm1           ! total ocean volume
-         ztmpk(:,:,jk,4) = surf(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+         DO_2D( 0, 0, 0, 0 )
+            ztmpk(ji,jj,jk,4) = surf(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
+         END_2D
       END DO
       
       ! global sum
@@ -315,14 +361,18 @@ CONTAINS
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '   dia_hsb_rst : initialise hsb at initial state '
             IF(lwp) WRITE(numout,*)
-            surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:)         ! initial ocean surface
-            ssh_ini(:,:) = ssh(:,:,Kmm)                          ! initial ssh
+            DO_2D( 0, 0, 0, 0 )
+               surf_ini(ji,jj) = e1e2t(ji,jj) * tmask_i(ji,jj)         ! initial ocean surface
+               ssh_ini(ji,jj) = ssh(ji,jj,Kmm)                          ! initial ssh
+            END_2D
             DO jk = 1, jpk
               ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
-               e3t_ini   (:,:,jk) = e3t(:,:,jk,Kmm)                      * tmask(:,:,jk)  ! initial vertical scale factors
-               tmask_ini (:,:,jk) = tmask(:,:,jk)                                       ! initial mask
-               hc_loc_ini(:,:,jk) = ts(:,:,jk,jp_tem,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial heat content
-               sc_loc_ini(:,:,jk) = ts(:,:,jk,jp_sal,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial salt content
+               DO_2D( 0, 0, 0, 0 )
+                  e3t_ini   (ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)  ! initial vertical scale factors
+                  tmask_ini (ji,jj,jk) = tmask(ji,jj,jk)                      ! initial mask
+                  hc_loc_ini(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)  ! initial heat content
+                  sc_loc_ini(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)  ! initial salt content
+               END_2D
             END DO
             frc_v = 0._wp                                           ! volume       trend due to forcing
             frc_t = 0._wp                                           ! heat content   -    -   -    -
@@ -334,13 +384,15 @@ CONTAINS
                         ssh_hc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm)   ! initial heat content in ssh
                         ssh_sc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm)   ! initial salt content in ssh
                      END DO
-                   END DO
-                ELSE
-                  ssh_hc_loc_ini(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm)   ! initial heat content in ssh
-                  ssh_sc_loc_ini(:,:) = ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm)   ! initial salt content in ssh
+                  END DO
+               ELSE
+                  DO_2D( 0, 0, 0, 0 )
+                     ssh_hc_loc_ini(ji,jj) = ts(ji,jj,1,jp_tem,Kmm) * ssh(ji,jj,Kmm)   ! initial heat content in ssh
+                     ssh_sc_loc_ini(ji,jj) = ts(ji,jj,1,jp_sal,Kmm) * ssh(ji,jj,Kmm)   ! initial salt content in ssh
+                  END_2D
                END IF
-               frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface
-               frc_wn_s = 0._wp                                       ! initial salt content misfit due to free surface
+               frc_wn_t = 0._wp    ! initial heat content misfit due to free surface
+               frc_wn_s = 0._wp    ! initial salt content misfit due to free surface
             ENDIF
          ENDIF
          !
@@ -388,6 +440,7 @@ CONTAINS
       INTEGER, INTENT(in) :: Kmm ! time level index
       !
       INTEGER ::   ierror, ios   ! local integer
+      INTEGER ::   ji, jj        ! loop index
       !!
       NAMELIST/namhsb/ ln_diahsb
       !!----------------------------------------------------------------------
@@ -427,7 +480,10 @@ CONTAINS
       ! ----------------------------------------------- !
       ! 2 - Time independant variables and file opening !
       ! ----------------------------------------------- !
-      surf(:,:) = e1e2t(:,:) * tmask_i(:,:)               ! masked surface grid cell area
+ 
+      DO_2D( 0, 0, 0, 0 )
+         surf(ji,jj) = e1e2t(ji,jj) * tmask_i(ji,jj)               ! masked surface grid cell area
+      END_2D
       surf_tot  = glob_sum( 'diahsb', surf(:,:) )         ! total ocean surface area
 
       IF( ln_bdy ) CALL ctl_warn( 'dia_hsb_init: heat/salt budget does not consider open boundary fluxes' )
diff --git a/src/OCE/DIA/diahth.F90 b/src/OCE/DIA/diahth.F90
index 0b8650acf..dcec1291f 100644
--- a/src/OCE/DIA/diahth.F90
+++ b/src/OCE/DIA/diahth.F90
@@ -86,22 +86,22 @@ CONTAINS
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
       INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index
       !!
-      INTEGER                      ::   ji, jj, jk            ! dummy loop arguments
-      REAL(wp)                     ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth
-      REAL(wp)                     ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth
-      REAL(wp)                     ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth
-      REAL(wp)                     ::   zztmp, zzdep          ! temporary scalars inside do loop
-      REAL(wp)                     ::   zu, zv, zw, zut, zvt  ! temporary workspace
-      REAL(wp), DIMENSION(jpi,jpj) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2 
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2     
-      REAL(wp), DIMENSION(jpi,jpj) ::   zrho10_3   ! MLD: rho = rho10m + zrho3      
-      REAL(wp), DIMENSION(jpi,jpj) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztinv      ! max of temperature inversion
-      REAL(wp), DIMENSION(jpi,jpj) ::   zdepinv    ! depth of temperature inversion
-      REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03
-      REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmaxdzT    ! max of dT/dz
-      REAL(wp), DIMENSION(jpi,jpj) ::   zdelr      ! delta rho equivalent to deltaT = 0.2
+      INTEGER                     ::   ji, jj, jk            ! dummy loop arguments
+      REAL(wp)                    ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth
+      REAL(wp)                    ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth
+      REAL(wp)                    ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth
+      REAL(wp)                    ::   zztmp, zzdep          ! temporary scalars inside do loop
+      REAL(wp)                    ::   zu, zv, zw, zut, zvt  ! temporary workspace
+      REAL(wp), DIMENSION(A2D(0)) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2 
+      REAL(wp), DIMENSION(A2D(0)) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2     
+      REAL(wp), DIMENSION(A2D(0)) ::   zrho10_3   ! MLD: rho = rho10m + zrho3      
+      REAL(wp), DIMENSION(A2D(0)) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
+      REAL(wp), DIMENSION(A2D(0)) ::   ztinv      ! max of temperature inversion
+      REAL(wp), DIMENSION(A2D(0)) ::   zdepinv    ! depth of temperature inversion
+      REAL(wp), DIMENSION(A2D(0)) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03
+      REAL(wp), DIMENSION(A2D(0)) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01
+      REAL(wp), DIMENSION(A2D(0)) ::   zmaxdzT    ! max of dT/dz
+      REAL(wp), DIMENSION(A2D(0)) ::   zdelr      ! delta rho equivalent to deltaT = 0.2
       !!----------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('dia_hth')
 
@@ -131,7 +131,7 @@ CONTAINS
          IF( iom_use( 'mlddzt' ) )   zmaxdzT(:,:) = 0._wp  
          IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' )   &
             &                    .OR. iom_use( 'mldr10_3' ) .OR. iom_use( 'pycndep'  ) ) THEN
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( 0, 0, 0, 0 )
                zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 
                hth     (ji,jj) = zztmp
                zabs2   (ji,jj) = zztmp
@@ -142,7 +142,7 @@ CONTAINS
          ENDIF
          IF( iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) ) THEN
             IF( nla10 > 1 ) THEN 
-               DO_2D( 1, 1, 1, 1 )
+               DO_2D( 0, 0, 0, 0 )
                   zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 
                   zrho0_3(ji,jj) = zztmp
                   zrho0_1(ji,jj) = zztmp
@@ -157,7 +157,7 @@ CONTAINS
             ! MLD: rho = rho(1) + zrho3                                     !
             ! MLD: rho = rho(1) + zrho1                                     !
             ! ------------------------------------------------------------- !
-            DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )   ! loop from bottom to 2
+            DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! loop from bottom to 2
                !
                zzdep = gdepw(ji,jj,jk,Kmm)
                zztmp = ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ) &
@@ -189,7 +189,7 @@ CONTAINS
             !
             ! Preliminary computation
             ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( 0, 0, 0, 0 )
                IF( tmask(ji,jj,nla10) == 1. ) THEN
                   zu  =  1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80   * ts(ji,jj,nla10,jp_sal,Kmm)  &
                      &           - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)   &
@@ -213,7 +213,7 @@ CONTAINS
             ! temperature inversion: max( 0, max of tn - tn(10m) )          !
             ! depth of temperature inversion                                !
             ! ------------------------------------------------------------- !
-            DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )   ! loop from bottom to nlb10
+            DO_3DS( 0, 0, 0, 0, jpkm1, nlb10, -1 )   ! loop from bottom to nlb10
                !
                zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1)
                !
@@ -305,13 +305,16 @@ CONTAINS
       !
       INTEGER  :: ji, jj, jk, iid
       REAL(wp) :: zztmp, zzdep
-      INTEGER, DIMENSION(jpi,jpj) :: iktem
+      INTEGER, DIMENSION(A2D(0)) :: iktem
       
       ! --------------------------------------- !
       ! search deepest level above ptem         !
       ! --------------------------------------- !
-      iktem(:,:) = 1
-      DO_3D( 1, 1, 1, 1, 1, jpkm1 )   ! beware temperature is not always decreasing with depth => loop from top to bottom
+      DO_2D( 0, 0, 0, 0 )
+         iktem(ji,jj) = 1
+      END_2D
+
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! beware temperature is not always decreasing with depth => loop from top to bottom
          zztmp = ts(ji,jj,jk,jp_tem,Kmm)
          IF( zztmp >= ptem )   iktem(ji,jj) = jk
       END_3D
@@ -319,7 +322,7 @@ CONTAINS
       ! ------------------------------- !
       !  Depth of ptem isotherm         !
       ! ------------------------------- !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( 0, 0, 0, 0 )
          !
          zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)       ! depth of the ocean bottom
          !
@@ -346,18 +349,29 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj),     INTENT(inout) ::   phtc  
       !
       INTEGER  ::   ji, jj, jk, ik
-      REAL(wp), DIMENSION(jpi,jpj) ::   zthick
-      INTEGER , DIMENSION(jpi,jpj) ::   ilevel
+      REAL(wp), DIMENSION(A2D(0)) ::   zthick
+      INTEGER , DIMENSION(A2D(0)) ::   ilevel
 
 
       ! surface boundary condition
       
-      IF( .NOT. ln_linssh ) THEN   ;   zthick(:,:) = 0._wp          ;   phtc(:,:) = 0._wp                                   
-      ELSE                         ;   zthick(:,:) = ssh(:,:,Kmm)   ;   phtc(:,:) = pt(:,:,1) * ssh(:,:,Kmm) * tmask(:,:,1)   
+      IF( .NOT. ln_linssh ) THEN 
+         DO_2D( 0, 0, 0, 0 )
+            zthick(ji,jj) = 0._wp 
+            phtc  (ji,jj) = 0._wp                                   
+         END_2D
+      ELSE                
+         DO_2D( 0, 0, 0, 0 )
+            zthick(ji,jj) = ssh(ji,jj,Kmm)   
+            phtc  (ji,jj) = pt(ji,jj,1) * ssh(ji,jj,Kmm) * tmask(ji,jj,1)   
+         END_2D
       ENDIF
       !
-      ilevel(:,:) = 1
-      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+      DO_2D( 0, 0, 0, 0 )
+         ilevel(ji,jj) = 1
+      END_2D
+      !
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          IF( ( gdepw(ji,jj,jk+1,Kmm) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN
              ilevel(ji,jj) = jk+1
              zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm)
@@ -365,7 +379,7 @@ CONTAINS
          ENDIF
       END_3D
       !
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( 0, 0, 0, 0 )
          ik = ilevel(ji,jj)
          IF( tmask(ji,jj,ik) == 1 ) THEN
             zthick(ji,jj) = MIN ( gdepw(ji,jj,ik+1,Kmm), pdep ) - zthick(ji,jj)   ! remaining thickness to reach dephw pdep
diff --git a/src/OCE/DIA/diamlr.F90 b/src/OCE/DIA/diamlr.F90
index edd8e19f9..dbcf5bbea 100644
--- a/src/OCE/DIA/diamlr.F90
+++ b/src/OCE/DIA/diamlr.F90
@@ -6,7 +6,7 @@ MODULE diamlr
    !! History :  4.0  !  2019  (S. Mueller)   Original code
    !!----------------------------------------------------------------------
 
-   USE par_oce        , ONLY :   wp, jpi, jpj
+   USE par_oce        , ONLY :   wp, jpi, jpj, ntsi, ntei, ntsj, ntej 
    USE phycst         , ONLY :   rpi
    USE dom_oce        , ONLY :   adatrj
    USE tide_mod
@@ -407,8 +407,9 @@ CONTAINS
       !! ** Purpose : update time used in multiple-linear-regression analysis
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj) ::   zadatrj2d
+      REAL(wp), DIMENSION(A2D(0)) ::   zadatrj2d
       !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj
 
       IF( ln_timing )   CALL timing_start('dia_mlr')
 
@@ -417,7 +418,9 @@ CONTAINS
       !
       ! A 2-dimensional field of constant value is sent, and subsequently used directly 
       ! or transformed to a scalar or a constant 3-dimensional field as required.
-      zadatrj2d(:,:) = adatrj*86400.0_wp
+      DO_2D( 0, 0, 0, 0 )
+         zadatrj2d(ji,jj) = adatrj*86400.0_wp
+      END_2D
       IF ( iom_use('diamlr_time') ) CALL iom_put('diamlr_time', zadatrj2d)
       !
       IF( ln_timing )   CALL timing_stop('dia_mlr')
diff --git a/src/OCE/DIA/diaptr.F90 b/src/OCE/DIA/diaptr.F90
index 8962bd371..6c4f48a3d 100644
--- a/src/OCE/DIA/diaptr.F90
+++ b/src/OCE/DIA/diaptr.F90
@@ -78,7 +78,7 @@ CONTAINS
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in)           ::   kt     ! ocean time-step index
       INTEGER                         , INTENT(in)           ::   Kmm    ! time level index
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
+      REAL(wp), DIMENSION(A2D(0),jpk) , INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
       !!----------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('dia_ptr')
@@ -110,13 +110,13 @@ CONTAINS
       !!----------------------------------------------------------------------
       !! ** Purpose : Calculate diagnostics and send to XIOS
       !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in)           ::   kt     ! ocean time-step index
-      INTEGER                         , INTENT(in)           ::   Kmm    ! time level index
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
+      INTEGER                   , INTENT(in)           ::   kt     ! ocean time-step index
+      INTEGER                   , INTENT(in)           ::   Kmm    ! time level index
+      REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport (used only by PRESENT)
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace
-      REAL(wp), DIMENSION(jpj)      ::  zvsum, ztsum, zssum   ! 1D workspace
+      REAL(wp), DIMENSION(Nis0:Nie0,Njs0:Nje0)   ::  z2d                   ! 2D workspace
+      REAL(wp), DIMENSION(          Njs0:Nje0)   ::  zvsum, ztsum, zssum   ! 1D workspace
       !
       !overturning calculation
       REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   sjk, r1_sjk, v_msf  ! i-mean i-k-surface and its inverse
@@ -126,19 +126,19 @@ CONTAINS
       REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   z3dtr
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( z3dtr(jpi,jpj,nbasin) )
+      ALLOCATE( z3dtr(Nis0:Nie0,Njs0:Nje0,nbasin) )
 
       IF( PRESENT( pvtr ) ) THEN
          IF( iom_use( 'zomsf' ) ) THEN    ! effective MSF
-            ALLOCATE( z4d1(jpi,jpj,jpk,nbasin) )
+            ALLOCATE( z4d1(Nis0:Nie0,Njs0:Nje0,jpk,nbasin) )
             !
             DO jn = 1, nbasin                                    ! by sub-basins
-               z4d1(1,:,:,jn) =  pvtr_int(:,:,jp_vtr,jn)                  ! zonal cumulative effective transport excluding closed seas
+               z4d1(Nis0,:,:,jn) =  pvtr_int(:,:,jp_vtr,jn)                  ! zonal cumulative effective transport excluding closed seas
                DO jk = jpkm1, 1, -1
-                  z4d1(1,:,jk,jn) = z4d1(1,:,jk+1,jn) - z4d1(1,:,jk,jn)    ! effective j-Stream-Function (MSF)
+                  z4d1(Nis0,:,jk,jn) = z4d1(Nis0,:,jk+1,jn) - z4d1(Nis0,:,jk,jn)    ! effective j-Stream-Function (MSF)
                END DO
-               DO ji = 2, jpi
-                  z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
+               DO ji = Nis0+1, Nie0
+                  z4d1(ji,:,:,jn) = z4d1(Nis0,:,:,jn)
                ENDDO
             END DO
             CALL iom_put( 'zomsf', z4d1 * rc_sv )
@@ -146,8 +146,8 @@ CONTAINS
             DEALLOCATE( z4d1 )
          ENDIF
          IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
-            ALLOCATE( sjk(jpj,jpk,nbasin), r1_sjk(jpj,jpk,nbasin), v_msf(jpj,jpk,nbasin),   &
-               &      zt_jk(jpj,jpk,nbasin), zs_jk(jpj,jpk,nbasin) )
+            ALLOCATE( sjk(  Njs0:Nje0,jpk,nbasin), r1_sjk(Njs0:Nje0,jpk,nbasin), v_msf(Njs0:Nje0,jpk,nbasin),   &
+               &      zt_jk(Njs0:Nje0,jpk,nbasin), zs_jk( Njs0:Nje0,jpk,nbasin) )
             !
             DO jn = 1, nbasin
                sjk(:,:,jn) = pvtr_int(:,:,jp_msk,jn)
@@ -162,16 +162,16 @@ CONTAINS
                !
             ENDDO
             DO jn = 1, nbasin
-               z3dtr(1,:,jn) = hstr_ove(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
-               DO ji = 2, jpi
-                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               z3dtr(Nis0,:,jn) = hstr_ove(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
+               DO ji = Nis0+1, Nie0
+                  z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'sophtove', z3dtr )
             DO jn = 1, nbasin
-               z3dtr(1,:,jn) = hstr_ove(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
-               DO ji = 2, jpi
-                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               z3dtr(Nis0,:,jn) = hstr_ove(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
+               DO ji = Nis0+1, Nie0
+                  z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'sopstove', z3dtr )
@@ -181,7 +181,7 @@ CONTAINS
 
          IF( iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
             ! Calculate barotropic heat and salt transport here
-            ALLOCATE( sjk(jpj,1,nbasin), r1_sjk(jpj,1,nbasin) )
+            ALLOCATE( sjk(A1Dj(0),1,nbasin), r1_sjk(A1Dj(0),1,nbasin) )
             !
             DO jn = 1, nbasin
                sjk(:,1,jn) = SUM( pvtr_int(:,:,jp_msk,jn), 2 )
@@ -196,16 +196,16 @@ CONTAINS
                !
             ENDDO
             DO jn = 1, nbasin
-               z3dtr(1,:,jn) = hstr_btr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
-               DO ji = 2, jpi
-                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               z3dtr(Nis0,:,jn) = hstr_btr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
+               DO ji = Nis0+1, Nie0
+                  z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'sophtbtr', z3dtr )
             DO jn = 1, nbasin
-               z3dtr(1,:,jn) = hstr_btr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
-               DO ji = 2, jpi
-                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               z3dtr(Nis0,:,jn) = hstr_btr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
+               DO ji = Nis0+1, Nie0
+                  z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'sopstbtr', z3dtr )
@@ -218,28 +218,28 @@ CONTAINS
          pvtr_int(:,:,:,:) = 0._wp
       ELSE
          IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN    ! i-mean i-k-surface
-            ALLOCATE( z4d1(jpi,jpj,jpk,nbasin), z4d2(jpi,jpj,jpk,nbasin) )
+            ALLOCATE( z4d1(Nis0:Nie0,Njs0:Nje0,jpk,nbasin), z4d2(Nis0:Nie0,Njs0:Nje0,jpk,nbasin) )
             !
             DO jn = 1, nbasin
-               z4d1(1,:,:,jn) = pzon_int(:,:,jp_msk,jn)
-               DO ji = 2, jpi
-                  z4d1(ji,:,:,jn) = z4d1(1,:,:,jn)
+               z4d1(Nis0,:,:,jn) = pzon_int(:,:,jp_msk,jn)
+               DO ji = Nis0+1, Nie0
+                  z4d1(ji,:,:,jn) = z4d1(Nis0,:,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'zosrf', z4d1 )
             !
             DO jn = 1, nbasin
-               z4d2(1,:,:,jn) = pzon_int(:,:,jp_tem,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 )
-               DO ji = 2, jpi
-                  z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
+               z4d2(Nis0,:,:,jn) = pzon_int(:,:,jp_tem,jn) / MAX( z4d1(Nis0,:,:,jn), 10.e-15 )
+               DO ji = Nis0+1, Nie0
+                  z4d2(ji,:,:,jn) = z4d2(Nis0,:,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'zotem', z4d2 )
             !
             DO jn = 1, nbasin
-               z4d2(1,:,:,jn) = pzon_int(:,:,jp_sal,jn) / MAX( z4d1(1,:,:,jn), 10.e-15 )
-               DO ji = 2, jpi
-                  z4d2(ji,:,:,jn) = z4d2(1,:,:,jn)
+               z4d2(Nis0,:,:,jn) = pzon_int(:,:,jp_sal,jn) / MAX( z4d1(Nis0,:,:,jn), 10.e-15 )
+               DO ji = Nis0+1, Nie0
+                  z4d2(ji,:,:,jn) = z4d2(Nis0,:,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'zosal', z4d2 )
@@ -251,16 +251,16 @@ CONTAINS
          IF( iom_use( 'sophtadv' ) .OR. iom_use( 'sopstadv' ) ) THEN
             !
             DO jn = 1, nbasin
-               z3dtr(1,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
-               DO ji = 2, jpi
-                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               z3dtr(Nis0,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
+               DO ji = Nis0+1, Nie0
+                  z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'sophtadv', z3dtr )
             DO jn = 1, nbasin
-               z3dtr(1,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
-               DO ji = 2, jpi
-                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               z3dtr(Nis0,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
+               DO ji = Nis0+1, Nie0
+                  z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'sopstadv', z3dtr )
@@ -269,16 +269,16 @@ CONTAINS
          IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) THEN
             !
             DO jn = 1, nbasin
-               z3dtr(1,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
-               DO ji = 2, jpi
-                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               z3dtr(Nis0,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
+               DO ji = Nis0+1, Nie0
+                  z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'sophtldf', z3dtr )
             DO jn = 1, nbasin
-               z3dtr(1,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
-               DO ji = 2, jpi
-                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               z3dtr(Nis0,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
+               DO ji = Nis0+1, Nie0
+                  z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'sopstldf', z3dtr )
@@ -287,16 +287,16 @@ CONTAINS
          IF( iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) ) THEN
             !
             DO jn = 1, nbasin
-               z3dtr(1,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
-               DO ji = 2, jpi
-                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               z3dtr(Nis0,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
+               DO ji = Nis0+1, Nie0
+                  z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'sophteiv', z3dtr )
             DO jn = 1, nbasin
-               z3dtr(1,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
-               DO ji = 2, jpi
-                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               z3dtr(Nis0,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
+               DO ji = Nis0+1, Nie0
+                  z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'sopsteiv', z3dtr )
@@ -304,16 +304,16 @@ CONTAINS
          !
          IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
              DO jn = 1, nbasin
-                z3dtr(1,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
-                DO ji = 2, jpi
-                   z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+                z3dtr(Nis0,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW)
+                DO ji = Nis0+1, Nie0
+                   z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
                 ENDDO
              ENDDO
              CALL iom_put( 'sophtvtr', z3dtr )
              DO jn = 1, nbasin
-               z3dtr(1,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
-               DO ji = 2, jpi
-                  z3dtr(ji,:,jn) = z3dtr(1,:,jn)
+               z3dtr(Nis0,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg)
+               DO ji = Nis0+1, Nie0
+                  z3dtr(ji,:,jn) = z3dtr(Nis0,:,jn)
                ENDDO
             ENDDO
             CALL iom_put( 'sopstvtr', z3dtr )
@@ -349,8 +349,8 @@ CONTAINS
       !! ** Action  : pvtr_int - terms for volume streamfunction, heat/salt transport barotropic/overturning terms
       !!              pzon_int - terms for i mean temperature/salinity
       !!----------------------------------------------------------------------
-      INTEGER                     , INTENT(in)           :: Kmm          ! time level index
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in), OPTIONAL :: pvtr         ! j-effective transport
+      INTEGER                        , INTENT(in)           :: Kmm          ! time level index
+      REAL(wp), DIMENSION(A2D(0),jpk), INTENT(in), OPTIONAL :: pvtr         ! j-effective transport
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE            :: zmask        ! 3D workspace
       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE          :: zts          ! 4D workspace
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE            :: sjk, v_msf   ! Zonal sum: i-k surface area, j-effective transport
@@ -362,7 +362,7 @@ CONTAINS
       IF( PRESENT( pvtr ) ) THEN
          ! i sum of effective j transport excluding closed seas
          IF( iom_use( 'zomsf' ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
-            ALLOCATE( v_msf(A1Dj(nn_hls),jpk,nbasin) )
+            ALLOCATE( v_msf(A1Dj(0),jpk,nbasin) )
 
             DO jn = 1, nbasin
                v_msf(:,:,jn) = ptr_sjk( pvtr(:,:,:), btmsk34(:,:,jn) )
@@ -374,16 +374,16 @@ CONTAINS
          ENDIF
 
          ! i sum of j surface area, j surface area - temperature/salinity product on V grid
-         IF(  iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.   &
+         IF(  iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.          &
             & iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN
-            ALLOCATE( zmask(A2D(nn_hls),jpk), zts(A2D(nn_hls),jpk,jpts), &
-               &      sjk(A1Dj(nn_hls),jpk,nbasin), &
-               &      zt_jk(A1Dj(nn_hls),jpk,nbasin), zs_jk(A1Dj(nn_hls),jpk,nbasin) )
+            ALLOCATE( zmask( A2D(0),jpk       ), zts(   A2D(0),jpk,jpts  ),  &
+               &      sjk(  A1Dj(0),jpk,nbasin),                             &
+               &      zt_jk(A1Dj(0),jpk,nbasin), zs_jk(A1Dj(0),jpk,nbasin) )
 
             zmask(:,:,:) = 0._wp
             zts(:,:,:,:) = 0._wp
 
-            DO_3D( 1, 1, 1, 0, 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
                zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc
                zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc !Tracers averaged onto V grid
@@ -405,14 +405,14 @@ CONTAINS
       ELSE
          ! i sum of j surface area - temperature/salinity product on T grid
          IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN
-            ALLOCATE( zmask(A2D(nn_hls),jpk), zts(A2D(nn_hls),jpk,jpts), &
-               &      sjk(A1Dj(nn_hls),jpk,nbasin), &
-               &      zt_jk(A1Dj(nn_hls),jpk,nbasin), zs_jk(A1Dj(nn_hls),jpk,nbasin) )
+            ALLOCATE( zmask( A2D(0),jpk       ), zts(   A2D(0),jpk,jpts  ),   &
+               &      sjk(  A1Dj(0),jpk,nbasin),                              &
+               &      zt_jk(A1Dj(0),jpk,nbasin), zs_jk(A1Dj(0),jpk,nbasin) )
 
             zmask(:,:,:) = 0._wp
             zts(:,:,:,:) = 0._wp
 
-            DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                zsfc = e1t(ji,jj) * e3t(ji,jj,jk,Kmm)
                zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
                zts(ji,jj,jk,jp_tem) = ts(ji,jj,jk,jp_tem,Kmm) * zsfc
@@ -434,11 +434,12 @@ CONTAINS
 
          ! i-k sum of j surface area - temperature/salinity product on V grid
          IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
+            ! TODO: Can be A2D(0) once all dia_ptr_hst calls have arguments with consistent declarations
             ALLOCATE( zts(A2D(nn_hls),jpk,jpts) )
 
             zts(:,:,:,:) = 0._wp
 
-            DO_3D( 1, 1, 1, 0, 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
                zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc  !Tracers averaged onto V grid
                zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
@@ -538,11 +539,12 @@ CONTAINS
       !! Wrapper for heat and salt transport calculations to calculate them for each basin
       !! Called from all advection and/or diffusion routines
       !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in )  :: ktra  ! tracer index
-      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'/'eiv'
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in)   :: pvflx ! 3D input array of advection/diffusion
-      REAL(wp), DIMENSION(A1Dj(nn_hls),nbasin)                 :: zsj   !
-      INTEGER                                        :: jn    !
+      INTEGER                             , INTENT(in)  ::  ktra  ! tracer index
+      CHARACTER(len=3)                    , INTENT(in)  ::  cptr  ! transport type  'adv'/'ldf'/'eiv'
+      ! TODO: Can be A2D(0) once all dia_ptr_hst calls have arguments with consistent declarations
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in)  ::  pvflx ! 3D input array of advection/diffusion
+      REAL(wp), DIMENSION(A1Dj(0),nbasin)               ::  zsj   !
+      INTEGER                                           ::  jn    !
 
       DO jn = 1, nbasin
          zsj(:,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
@@ -576,13 +578,13 @@ CONTAINS
       !!
       !! ** Action  : phstr
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpj,nbasin) , INTENT(inout)         ::  phstr  !
-      REAL(wp), DIMENSION(A1Dj(nn_hls),nbasin), INTENT(in)            ::  pva    !
+      REAL(wp), DIMENSION(Njs0:Nje0,nbasin), INTENT(inout)  ::  phstr  !
+      REAL(wp), DIMENSION(A1Dj(0)  ,nbasin), INTENT(in   )  ::  pva    !
       INTEGER                                               ::  jj
 #if ! defined key_mpi_off
-      INTEGER, DIMENSION(1)           ::  ish1d
-      INTEGER, DIMENSION(2)           ::  ish2d
-      REAL(wp), DIMENSION(jpj*nbasin) ::  zwork
+      INTEGER,  DIMENSION(1)               ::  ish1d
+      INTEGER,  DIMENSION(2)               ::  ish2d
+      REAL(wp), DIMENSION(:), ALLOCATABLE  ::  zwork
 #endif
 
       DO jj = ntsj, ntej
@@ -591,11 +593,13 @@ CONTAINS
 
 #if ! defined key_mpi_off
       IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
-         ish1d(1) = jpj*nbasin
-         ish2d(1) = jpj ; ish2d(2) = nbasin
+         ALLOCATE( zwork(Nj_0*nbasin) )
+         ish1d(1) = Nj_0*nbasin
+         ish2d(1) = Nj_0 ; ish2d(2) = nbasin
          zwork(:) = RESHAPE( phstr(:,:), ish1d )
          CALL mpp_sum( 'diaptr', zwork, ish1d(1), ncomm_znl )
          phstr(:,:) = RESHAPE( zwork, ish2d )
+         DEALLOCATE( zwork )
       ENDIF
 #endif
    END SUBROUTINE ptr_sum_2d
@@ -612,13 +616,13 @@ CONTAINS
       !!
       !! ** Action  : phstr
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpj,jpk,nbasin) , INTENT(inout)     ::  phstr  !
-      REAL(wp), DIMENSION(A1Dj(nn_hls),jpk,nbasin), INTENT(in)        ::  pva    !
-      INTEGER                                               ::  jj, jk
+      REAL(wp), DIMENSION(Njs0:Nje0,jpk,nbasin), INTENT(inout)   ::  phstr  !
+      REAL(wp), DIMENSION(A1Dj(0)  ,jpk,nbasin), INTENT(in   )   ::  pva    !
+      INTEGER                                                    ::  jj, jk
 #if ! defined key_mpi_off
-      INTEGER, DIMENSION(1)              ::  ish1d
-      INTEGER, DIMENSION(3)              ::  ish3d
-      REAL(wp), DIMENSION(jpj*jpk*nbasin)  ::  zwork
+      INTEGER,  DIMENSION(1)               ::  ish1d
+      INTEGER,  DIMENSION(3)               ::  ish3d
+      REAL(wp), DIMENSION(:), ALLOCATABLE  ::  zwork
 #endif
 
       DO jk = 1, jpk
@@ -629,11 +633,13 @@ CONTAINS
 
 #if ! defined key_mpi_off
       IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN
-         ish1d(1) = jpj*jpk*nbasin
-         ish3d(1) = jpj ; ish3d(2) = jpk ; ish3d(3) = nbasin
+         ALLOCATE( zwork(Nj_0*jpk*nbasin) )
+         ish1d(1) = Nj_0*jpk*nbasin
+         ish3d(1) = Nj_0 ; ish3d(2) = jpk ; ish3d(3) = nbasin
          zwork(:) = RESHAPE( phstr(:,:,:), ish1d )
          CALL mpp_sum( 'diaptr', zwork, ish1d(1), ncomm_znl )
          phstr(:,:,:) = RESHAPE( zwork, ish3d )
+         DEALLOCATE( zwork )
       ENDIF
 #endif
    END SUBROUTINE ptr_sum_3d
@@ -651,13 +657,13 @@ CONTAINS
       ! nbasin has been initialized in iom_init to define the axis "basin"
       !
       IF( .NOT. ALLOCATED( btmsk ) ) THEN
-         ALLOCATE( btmsk(jpi,jpj,nbasin)    , btmsk34(jpi,jpj,nbasin),   &
-            &      hstr_adv(jpj,jpts,nbasin), hstr_eiv(jpj,jpts,nbasin), &
-            &      hstr_ove(jpj,jpts,nbasin), hstr_btr(jpj,jpts,nbasin), &
-            &      hstr_ldf(jpj,jpts,nbasin), hstr_vtr(jpj,jpts,nbasin), STAT=ierr(1)  )
+         ALLOCATE( btmsk(jpi,jpj,nbasin)          , btmsk34(jpi,jpj,nbasin),         &
+            &      hstr_adv(Njs0:Nje0,jpts,nbasin), hstr_eiv(Njs0:Nje0,jpts,nbasin), &
+            &      hstr_ove(Njs0:Nje0,jpts,nbasin), hstr_btr(Njs0:Nje0,jpts,nbasin), &
+            &      hstr_ldf(Njs0:Nje0,jpts,nbasin), hstr_vtr(Njs0:Nje0,jpts,nbasin), STAT=ierr(1)  )
             !
-         ALLOCATE( pvtr_int(jpj,jpk,jpts+2,nbasin), &
-            &      pzon_int(jpj,jpk,jpts+1,nbasin), STAT=ierr(2) )
+         ALLOCATE( pvtr_int(Njs0:Nje0,jpk,jpts+2,nbasin), &
+            &      pzon_int(Njs0:Nje0,jpk,jpts+1,nbasin), STAT=ierr(2) )
          !
          dia_ptr_alloc = MAXVAL( ierr )
          CALL mpp_sum( 'diaptr', dia_ptr_alloc )
@@ -677,11 +683,12 @@ CONTAINS
       !!
       !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
       !!----------------------------------------------------------------------
+      ! TODO: Can be A2D(0) once all dia_ptr_hst calls have arguments with consistent declarations
       REAL(wp), INTENT(in), DIMENSION(A2D(nn_hls),jpk)  ::   pvflx  ! mask flux array at V-point
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj)  ::   pmsk   ! Optional 2D basin mask
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj        )  ::   pmsk   ! Optional 2D basin mask
       !
-      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments
-      REAL(wp), DIMENSION(A1Dj(nn_hls)) :: p_fval  ! function value
+      INTEGER                      :: ji, jj, jk  ! dummy loop arguments
+      REAL(wp), DIMENSION(A1Dj(0)) :: p_fval      ! function value
       !!--------------------------------------------------------------------
       !
       p_fval(:) = 0._wp
@@ -702,11 +709,12 @@ CONTAINS
       !!
       !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
       !!----------------------------------------------------------------------
-      REAL(wp) , INTENT(in), DIMENSION(A2D(nn_hls))     ::   pvflx  ! mask flux array at V-point
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask
+      ! TODO: Can be A2D(0) once all dia_ptr_hst calls have arguments with consistent declarations
+      REAL(wp) , INTENT(in), DIMENSION(A2D(nn_hls))  ::   pvflx  ! mask flux array at V-point
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj    )  ::   pmsk   ! Optional 2D basin mask
       !
-      INTEGER                  ::   ji,jj       ! dummy loop arguments
-      REAL(wp), DIMENSION(A1Dj(nn_hls)) :: p_fval ! function value
+      INTEGER                      :: ji, jj  ! dummy loop arguments
+      REAL(wp), DIMENSION(A1Dj(0)) :: p_fval  ! function value
       !!--------------------------------------------------------------------
       !
       p_fval(:) = 0._wp
@@ -725,14 +733,13 @@ CONTAINS
       !!
       !! ** Action  : - p_fval: j-cumulated sum of pva
       !!----------------------------------------------------------------------
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)  ::   pva   ! mask flux array at V-point
+      REAL(wp) , INTENT(in), DIMENSION(A2D(0))  ::   pva   ! mask flux array at V-point
       !
-      INTEGER                  ::   ji,jj,jc       ! dummy loop arguments
-      INTEGER                  ::   ijpj        ! ???
-      REAL(wp), DIMENSION(jpi,jpj) :: p_fval ! function value
+      INTEGER                     ::   ji,jj,jc  ! dummy loop arguments
+      INTEGER                     ::   ijpj      ! ???
+      REAL(wp), DIMENSION(A2D(0)) :: p_fval      ! function value
       !!--------------------------------------------------------------------
       !
-      ijpj = jpj  ! ???
       p_fval(:,:) = 0._wp
       DO jc = 1, jpnj ! looping over all processors in j axis
          DO_2D( 0, 0, 0, 0 )
@@ -756,11 +763,11 @@ CONTAINS
       !!----------------------------------------------------------------------
       !!
       IMPLICIT none
-      REAL(wp) , INTENT(in), DIMENSION(A2D(nn_hls),jpk) ::   pta    ! mask flux array at V-point
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask
+      REAL(wp) , INTENT(in), DIMENSION(A2D(0) ,jpk) ::   pta    ! mask flux array at V-point
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj    ) ::   pmsk   ! Optional 2D basin mask
       !!
       INTEGER                           :: ji, jj, jk ! dummy loop arguments
-      REAL(wp), DIMENSION(A1Dj(nn_hls),jpk) :: p_fval     ! return function value
+      REAL(wp), DIMENSION(A1Dj(0),jpk)  :: p_fval     ! return function value
       !!--------------------------------------------------------------------
       !
       p_fval(:,:) = 0._wp
diff --git a/src/OCE/DIA/diawri.F90 b/src/OCE/DIA/diawri.F90
index 0deeb3370..b6dd96cda 100644
--- a/src/OCE/DIA/diawri.F90
+++ b/src/OCE/DIA/diawri.F90
@@ -135,8 +135,8 @@ CONTAINS
       ENDIF
 
       ! initialize arrays
-      z2d(:,:)   = 0._wp
-      z3d(:,:,:) = 0._wp
+      z2d(A2D(0))   = 0._wp
+      z3d(A2D(0),:) = 0._wp
       
       ! Output of initial vertical scale factor
       CALL iom_put("e3t_0", e3t_0(:,:,:) )
diff --git a/src/OCE/IOM/iom.F90 b/src/OCE/IOM/iom.F90
index d66b13d19..728c75ebc 100644
--- a/src/OCE/IOM/iom.F90
+++ b/src/OCE/IOM/iom.F90
@@ -1327,14 +1327,21 @@ CONTAINS
                IF( irankpv == 1 )        ishape(1:1) = SHAPE(pv_r1d)
                IF( irankpv == 2 )        ishape(1:2) = SHAPE(pv_r2d)
                IF( irankpv == 3 )        ishape(1:3) = SHAPE(pv_r3d)
+               ix1 = 1      ;   ix2 = icnt(1)   ;   iy1 = 1   ;   iy2 = icnt(2)         ! index of the array to be read
                ctmp1 = 'd'
             ELSE
-               IF( irankpv == 2 ) THEN
-                  ishape(1:2) = SHAPE(pv_r2d(Nis0:Nie0,Njs0:Nje0  ))   ;   ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0)'
-               ENDIF
-               IF( irankpv == 3 ) THEN
-                  ishape(1:3) = SHAPE(pv_r3d(Nis0:Nie0,Njs0:Nje0,:))   ;   ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0,:)'
+               IF( irankpv == 2 )   ishape(1:2) = SHAPE(pv_r2d)
+               IF( irankpv == 3 )   ishape(1:3) = SHAPE(pv_r3d)
+               IF( ishape(1) == jpi .AND. ishape(2) == jpj ) THEN
+                  ishape(1:2) = (/ Ni_0, Nj_0 /)
+                  ix1 = Nis0   ;   ix2 = Nie0      ;   iy1 = Njs0   ;   iy2 = Nje0      ! index of the array to be read
+                  ctmp1 = 'd(Nis0:Nie0,Njs0:Nje0'
+               ELSE
+                  ix1 = 1      ;   ix2 = icnt(1)   ;   iy1 = 1      ;   iy2 = icnt(2)   ! index of the array to be read
+                  ctmp1 = 'd(:,:'
                ENDIF
+               IF( irankpv == 3 )   ctmp1 = TRIM(ctmp1)//',:'
+               ctmp1 = TRIM(ctmp1)//')'
             ENDIF
             DO jl = 1, irankpv
                WRITE( ctmp2, FMT="(', ', i1,'): ', i5,' /= icnt(', i1,'):', i5)" ) jl, ishape(jl), jl, icnt(jl)
@@ -1347,11 +1354,6 @@ CONTAINS
          !-
          IF( idvar > 0 .AND. istop == nstop ) THEN   ! no additional errors until this point...
             !
-            ! find the right index of the array to be read
-            IF( idom /= jpdom_unknown ) THEN   ;   ix1 = Nis0   ;   ix2 = Nie0      ;   iy1 = Njs0   ;   iy2 = Nje0
-            ELSE                               ;   ix1 = 1      ;   ix2 = icnt(1)   ;   iy1 = 1      ;   iy2 = icnt(2)
-            ENDIF
-
             CALL iom_nf90_get( kiomid, idvar, inbdim, istart, icnt, ix1, ix2, iy1, iy2, pv_r1d, pv_r2d, pv_r3d )
 
             IF( istop == nstop ) THEN   ! no additional errors until this point...
@@ -1362,9 +1364,11 @@ CONTAINS
                zsgn = 1._wp
                IF( PRESENT(psgn   ) )   zsgn    = psgn
                !--- overlap areas and extra hallows (mpp)
-               IF(     PRESENT(pv_r2d) .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+               llok = idom /= jpdom_unknown .AND. cl_type /= 'Z'  &
+                  &   .AND.   ix1 == Nis0   .AND.   ix2 == Nie0   .AND.   iy1 == Njs0   .AND.   iy2 == Nje0
+               IF(     PRESENT(pv_r2d) .AND. llok ) THEN
                   CALL lbc_lnk( 'iom', pv_r2d, cl_type, zsgn, kfillmode = kfill )
-               ELSEIF( PRESENT(pv_r3d) .AND. idom /= jpdom_unknown .AND. cl_type /= 'Z' ) THEN
+               ELSEIF( PRESENT(pv_r3d) .AND. llok ) THEN
                   CALL lbc_lnk( 'iom', pv_r3d, cl_type, zsgn, kfillmode = kfill )
                ENDIF
                !
@@ -2336,14 +2340,14 @@ CONTAINS
             idb(jn) = -nn_hls                         ! Tile data offset (halo size)
          END DO
 
-         ! Tile_[ij]begin are defined with respect to the processor data domain, so data_[ij]begin is added
          CALL iom_set_domain_attr("grid_"//cdgrd, ntiles=nijtile,                                     &
-            & tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
+            & tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 1, &
             & tile_ni=ini(:), tile_nj=inj(:),                                                         &
             & tile_data_ibegin=idb(:), tile_data_jbegin=idb(:),                                       &
             & tile_data_ni=ini(:) - 2 * idb(:), tile_data_nj=inj(:) - 2 * idb(:))
+         idb(:) = 0
          CALL iom_set_domain_attr("grid_"//cdgrd//"_inner", ntiles=nijtile,                           &
-            & tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
+            & tile_ibegin=ntsi_a(1:nijtile) - nn_hls - 1, tile_jbegin=ntsj_a(1:nijtile) - nn_hls - 1, &
             & tile_ni=ini(:), tile_nj=inj(:),                                                         &
             & tile_data_ibegin=idb(:), tile_data_jbegin=idb(:),                                       &
             & tile_data_ni=ini(:) - 2 * idb(:), tile_data_nj=inj(:) - 2 * idb(:))
@@ -2453,7 +2457,7 @@ CONTAINS
 !      CALL dom_ngb( -168.53_wp, 65.03_wp, ix, iy, 'T' ) !  i-line that passes through Bering Strait: Reference latitude (used in plots)
       CALL dom_ngb( 180.0_wp, 90.0_wp, ix, iy, 'T' ) !  i-line that passes near the North Pole : Reference latitude (used in plots)
       CALL iom_set_domain_attr("gznl", ni_glo=Ni0glo, nj_glo=Nj0glo, ibegin=mig0(Nis0)-1, jbegin=mjg0(Njs0)-1, ni=Ni_0, nj=Nj_0)
-      CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj)
+      CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin=0, data_ni=Ni_0, data_jbegin=0, data_nj=Nj_0)
       CALL iom_set_domain_attr("gznl", lonvalue = real(zlon, dp),   &
          &                             latvalue = real(RESHAPE(plat(Nis0:Nie0, Njs0:Nje0),(/ Ni_0*Nj_0 /)),dp))
       CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=Nj0glo)
diff --git a/src/OCE/IOM/iom_nf90.F90 b/src/OCE/IOM/iom_nf90.F90
index c49d9538e..287d943a9 100644
--- a/src/OCE/IOM/iom_nf90.F90
+++ b/src/OCE/IOM/iom_nf90.F90
@@ -40,6 +40,8 @@ MODULE iom_nf90
       MODULE PROCEDURE iom_nf90_rp0123d_dp
    END INTERFACE
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
    !! $Id: iom_nf90.F90 14433 2021-02-11 08:06:49Z smasson $
@@ -544,7 +546,7 @@ CONTAINS
       INTEGER               :: idvar                ! variable id
       INTEGER               :: jd                   ! dimension loop counter
       INTEGER               :: ix1, ix2, iy1, iy2   ! subdomain indexes
-      INTEGER, DIMENSION(4) :: idimsz               ! dimensions size
+      INTEGER, DIMENSION(3) :: ishape               ! dimensions size
       INTEGER, DIMENSION(4) :: idimid               ! dimensions id
       CHARACTER(LEN=256)    :: clinfo               ! info character
       INTEGER               :: if90id               ! nf90 file identifier
@@ -627,11 +629,9 @@ CONTAINS
             itype = NF90_DOUBLE
          ENDIF
          IF( PRESENT(pv_r0d) ) THEN
-            CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype,                    &
-               &                              iom_file(kiomid)%nvid(idvar) ), clinfo )
+            CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype,                  iom_file(kiomid)%nvid(idvar) ), clinfo )
          ELSE
-            CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, idimid(1:idims),   &
-               &                              iom_file(kiomid)%nvid(idvar) ), clinfo )
+            CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cdvar), itype, idimid(1:idims), iom_file(kiomid)%nvid(idvar) ), clinfo )
          ENDIF
          lchunk = .false.
          IF( snc4set%luse .AND. idims == 4 )   lchunk = .true.
@@ -673,23 +673,13 @@ CONTAINS
          ENDIF
          ! on what kind of domain must the data be written?
          IF( PRESENT(pv_r2d) .OR. PRESENT(pv_r3d) ) THEN
-            idimsz(1:2) = iom_file(kiomid)%dimsz(1:2,idvar)
-            IF(     idimsz(1) == Ni_0 .AND. idimsz(2) == Nj_0 ) THEN
-               ix1 = Nis0   ;   ix2 = Nie0   ;   iy1 = Njs0   ;   iy2 = Nje0
-            ELSEIF( idimsz(1) == jpi  .AND. idimsz(2) == jpj  ) THEN
-               ix1 = 1      ;   ix2 = jpi    ;   iy1 = 1      ;   iy2 = jpj
-            ELSEIF( idimsz(1) == jpi  .AND. idimsz(2) == jpj  ) THEN
-               ix1 = 1      ;   ix2 = jpi    ;   iy1 = 1      ;   iy2 = jpj
-            ELSE
-               CALL ctl_stop( 'iom_nf90_rp0123d: should have been an impossible case...' )
-            ENDIF
 
             ! write dimension variables if it is not already done
             ! =============
             ! trick: is defined to 0 => dimension variable are defined but not yet written
             IF( iom_file(kiomid)%dimsz(1, 4) == 0 ) THEN   ! time_counter = 0
-               CALL iom_nf90_check(    NF90_PUT_VAR( if90id, 1,                            glamt(ix1:ix2, iy1:iy2) ), clinfo )
-               CALL iom_nf90_check(    NF90_PUT_VAR( if90id, 2,                            gphit(ix1:ix2, iy1:iy2) ), clinfo )
+               CALL iom_nf90_check(    NF90_PUT_VAR( if90id, 1,                                      glamt(A2D(0)) ), clinfo )
+               CALL iom_nf90_check(    NF90_PUT_VAR( if90id, 2,                                      gphit(A2D(0)) ), clinfo )
                SELECT CASE (iom_file(kiomid)%comp)
                CASE ('OCE')
                   CALL iom_nf90_check( NF90_PUT_VAR( if90id, 3,                                           gdept_1d ), clinfo )
@@ -704,6 +694,17 @@ CONTAINS
                iom_file(kiomid)%dimsz(1, 4) = 1   ! so we don't enter this IF case any more...
                IF(lwp) WRITE(numout,*) TRIM(clinfo)//' write dimension variables done'
             ENDIF
+
+            IF( PRESENT(pv_r2d) )  ishape(1:2) = SHAPE(pv_r2d)
+            IF( PRESENT(pv_r3d) )  ishape(1:3) = SHAPE(pv_r3d)
+            IF(     ishape(1) == Ni_0 .AND. ishape(2) == Nj_0 ) THEN
+               ix1 = 1      ;   ix2 = Ni_0   ;   iy1 = 1      ;   iy2 = Nj_0
+            ELSEIF( ishape(1) == jpi  .AND. ishape(2) == jpj  ) THEN
+               ix1 = Nis0   ;   ix2 = Nie0   ;   iy1 = Njs0   ;   iy2 = Nje0
+            ELSE
+               CALL ctl_stop( 'iom_nf90_rp0123d: should have been an impossible case...' )
+            ENDIF
+
          ENDIF
 
          ! write the data
@@ -712,7 +713,7 @@ CONTAINS
             CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r0d                    ), clinfo )
          ELSEIF( PRESENT(pv_r1d) ) THEN
             CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r1d(:)                 ), clinfo )
-         ELSEIF( PRESENT(pv_r2d) ) THEN
+         ELSEIF( PRESENT(pv_r2d) ) THEN     
             CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r2d(ix1:ix2,iy1:iy2)   ), clinfo )
          ELSEIF( PRESENT(pv_r3d) ) THEN
             CALL iom_nf90_check( NF90_PUT_VAR( if90id, idvar, pv_r3d(ix1:ix2,iy1:iy2,:) ), clinfo )
diff --git a/src/OCE/IOM/prtctl.F90 b/src/OCE/IOM/prtctl.F90
index fae7c0d80..6fcca35f9 100644
--- a/src/OCE/IOM/prtctl.F90
+++ b/src/OCE/IOM/prtctl.F90
@@ -137,9 +137,14 @@ CONTAINS
       INTEGER ::  jn, jl, kdir
       INTEGER ::  iis, iie, jjs, jje
       INTEGER ::  itra, inum
+      INTEGER, DIMENSION(4) ::  ishape
       REAL(2*wp) :: zsum1, zsum2, zvctl1, zvctl2
       !!----------------------------------------------------------------------
       !
+      IF( ( ktab2d_1 * ktab3d_1 * ktab4d_1 * ktab2d_2 * ktab3d_2 ) /= 0 ) THEN
+         CALL ctl_stop( 'prt_ctl is not working with tiles' )
+      ENDIF
+      
       ! Arrays, scalars initialization
       cl1  = ''
       cl2  = ''
@@ -157,12 +162,19 @@ CONTAINS
       
       ! Loop over each sub-domain, i.e. the total number of processors ijsplt
       DO jl = 1, SIZE(nall_ictls)
-
-         ! define shoter names...
-         iis = MAX( nall_ictls(jl), ntsi )
-         iie = MIN( nall_ictle(jl), ntei )
-         jjs = MAX( nall_jctls(jl), ntsj )
-         jje = MIN( nall_jctle(jl), ntej )
+         
+         IF( PRESENT(tab2d_1) )   ishape(1:2) = SHAPE(tab2d_1)
+         IF( PRESENT(tab3d_1) )   ishape(1:3) = SHAPE(tab3d_1)
+         IF( PRESENT(tab4d_1) )   ishape(1:4) = SHAPE(tab4d_1)
+         IF( ishape(1) == jpi .AND. ishape(2) == jpj ) THEN
+            iis = Nis0   ;   iie = Nie0        ;   jjs = Njs0   ;   jje = Nje0
+         ELSE
+            iis = 1      ;   iie = ishape(1)   ;   jjs = 1      ;   jje = ishape(2)
+         ENDIF
+         iis = MAX( nall_ictls(jl), iis )
+         iie = MIN( nall_ictle(jl), iie )
+         jjs = MAX( nall_jctls(jl), jjs )
+         jje = MIN( nall_jctle(jl), jje )
 
          IF( PRESENT(clinfo) ) THEN   ;   inum = numprt_top(jl)
          ELSE                         ;   inum = numprt_oce(jl)
@@ -188,32 +200,32 @@ CONTAINS
 
                ! 2D arrays
                IF( PRESENT(tab2d_1) ) THEN
-                  IF( PRESENT(mask1) ) THEN   ;   zsum1 = SUM( tab2d_1(iis:iie,jjs:jje) * mask1(iis:iie,jjs:jje,1) )
-                  ELSE                        ;   zsum1 = SUM( tab2d_1(iis:iie,jjs:jje)                            )
+                  IF( PRESENT(mask1) ) THEN   ;   zsum1 = SUM( tab2d_1(iis:iie,jjs:jje) * mask1(A2D(0),1) )
+                  ELSE                        ;   zsum1 = SUM( tab2d_1(iis:iie,jjs:jje)                   )
                   ENDIF
                ENDIF
                IF( PRESENT(tab2d_2) ) THEN
-                  IF( PRESENT(mask2) ) THEN   ;   zsum2 = SUM( tab2d_2(iis:iie,jjs:jje) * mask2(iis:iie,jjs:jje,1) )
-                  ELSE                        ;   zsum2 = SUM( tab2d_2(iis:iie,jjs:jje)                            )
+                  IF( PRESENT(mask2) ) THEN   ;   zsum2 = SUM( tab2d_2(iis:iie,jjs:jje) * mask2(A2D(0),1) )
+                  ELSE                        ;   zsum2 = SUM( tab2d_2(iis:iie,jjs:jje)                   )
                   ENDIF
                ENDIF
 
                ! 3D arrays
                IF( PRESENT(tab3d_1) ) THEN
-                  IF( PRESENT(mask1) ) THEN   ;   zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir) * mask1(iis:iie,jjs:jje,1:kdir) )
-                  ELSE                        ;   zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir)                                 )
+                  IF( PRESENT(mask1) ) THEN   ;   zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir) * mask1(A2D(0),1:kdir) )
+                  ELSE                        ;   zsum1 = SUM( tab3d_1(iis:iie,jjs:jje,1:kdir)                        )
                   ENDIF
                ENDIF
                IF( PRESENT(tab3d_2) ) THEN
-                  IF( PRESENT(mask2) ) THEN   ;   zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir) * mask2(iis:iie,jjs:jje,1:kdir) )
-                  ELSE                        ;   zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir)                                 )
+                  IF( PRESENT(mask2) ) THEN   ;   zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir) * mask2(A2D(0),1:kdir) )
+                  ELSE                        ;   zsum2 = SUM( tab3d_2(iis:iie,jjs:jje,1:kdir)                        )
                   ENDIF
                ENDIF
 
                ! 4D arrays
                IF( PRESENT(tab4d_1) ) THEN
-                  IF( PRESENT(mask1) ) THEN   ;   zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn) * mask1(iis:iie,jjs:jje,1:kdir) )
-                  ELSE                        ;   zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn)                                 )
+                  IF( PRESENT(mask1) ) THEN   ;   zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn) * mask1(A2D(0),1:kdir) )
+                  ELSE                        ;   zsum1 = SUM( tab4d_1(iis:iie,jjs:jje,1:kdir,jn)                        )
                   ENDIF
                ENDIF
 
diff --git a/src/OCE/TRA/eosbn2.F90 b/src/OCE/TRA/eosbn2.F90
index 81fda25e0..1f808dfb9 100644
--- a/src/OCE/TRA/eosbn2.F90
+++ b/src/OCE/TRA/eosbn2.F90
@@ -1336,10 +1336,10 @@ CONTAINS
       !!                    pab_pe(:,:,:,jp_tem) is alpha_pe
       !!                    pab_pe(:,:,:,jp_sal) is beta_pe
       !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts     ! pot. temperature & salinity
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab_pe  ! alpha_pe and beta_pe
-      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   ppen     ! potential energy anomaly
+      INTEGER                     , INTENT(in   ) ::   Kmm   ! time level index
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts     ! pot. temperature & salinity
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(  out) ::   pab_pe  ! alpha_pe and beta_pe
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(  out) ::   ppen     ! potential energy anomaly
       !
       INTEGER  ::   ji, jj, jk                ! dummy loop indices
       REAL(wp) ::   zt , zh , zs , ztm        ! local scalars
@@ -1352,7 +1352,7 @@ CONTAINS
       !
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             !
             zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
             zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
@@ -1411,9 +1411,9 @@ CONTAINS
          !
       CASE( np_seos )                !==  Vallis (2006) simplified EOS  ==!
          !
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
             zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0)
-            zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0)
+            zs  = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0)
             zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point
             ztm = tmask(ji,jj,jk)                ! tmask
             zn  = 0.5_wp * zh * r1_rho0 * ztm
diff --git a/src/OCE/TRA/traadv.F90 b/src/OCE/TRA/traadv.F90
index a55ffff58..190a6a7ad 100644
--- a/src/OCE/TRA/traadv.F90
+++ b/src/OCE/TRA/traadv.F90
@@ -178,7 +178,7 @@ CONTAINS
          !
 !!gm ???
          ! TEMP: [tiling] This copy-in not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
-         IF( l_diaptr ) CALL dia_ptr( kt, Kmm, zvv(A2D(nn_hls),:) )                                    ! diagnose the effective MSF
+         IF( l_diaptr ) CALL dia_ptr( kt, Kmm, zvv(A2D(0),:) )                                    ! diagnose the effective MSF
 !!gm ???
          !
 
diff --git a/src/OCE/TRD/trddyn.F90 b/src/OCE/TRD/trddyn.F90
index 32e95a865..dd2daf7ef 100644
--- a/src/OCE/TRD/trddyn.F90
+++ b/src/OCE/TRD/trddyn.F90
@@ -56,10 +56,13 @@ CONTAINS
       INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
       INTEGER                   , INTENT(in   ) ::   kt             ! time step
       INTEGER                   , INTENT(in   ) ::   Kmm            ! time level index
+      INTEGER                                   ::   ji, jj, jk     ! lopp indices
       !!----------------------------------------------------------------------
       !
-      putrd(:,:,:) = putrd(:,:,:) * umask(:,:,:)                       ! mask the trends
-      pvtrd(:,:,:) = pvtrd(:,:,:) * vmask(:,:,:)
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         putrd(ji,jj,jk) = putrd(ji,jj,jk) * umask(ji,jj,jk)                       ! mask the trends
+         pvtrd(ji,jj,jk) = pvtrd(ji,jj,jk) * vmask(ji,jj,jk)
+      END_3D
       !
 
 !!gm NB : here a lbc_lnk should probably be added
@@ -120,10 +123,10 @@ CONTAINS
                               CALL iom_put( "vtrd_rvo", pvtrd )
       CASE( jpdyn_keg )   ;   CALL iom_put( "utrd_keg", putrd )    ! Kinetic Energy gradient (or had)
                               CALL iom_put( "vtrd_keg", pvtrd )
-                              ALLOCATE( z3dx(jpi,jpj,jpk) , z3dy(jpi,jpj,jpk) )
-                              z3dx(:,:,:) = 0._wp                  ! U.dxU & V.dyV (approximation)
-                              z3dy(:,:,:) = 0._wp
-                              DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! no mask as un,vn are masked
+                              ALLOCATE( z3dx(A2D(0),jpk) , z3dy(A2D(0),jpk) )   ! U.dxU & V.dyV (approximation)
+                              z3dx(A2D(0),jpk) = 0._wp
+                              z3dy(A2D(0),jpk) = 0._wp
+                              DO_3D( 0, 0, 0, 0, 1, jpkm1 )                          ! no mask as un,vn are masked
                                  z3dx(ji,jj,jk) = uu(ji,jj,jk,Kmm) * ( uu(ji+1,jj,jk,Kmm) - uu(ji-1,jj,jk,Kmm) ) / ( 2._wp * e1u(ji,jj) )
                                  z3dy(ji,jj,jk) = vv(ji,jj,jk,Kmm) * ( vv(ji,jj+1,jk,Kmm) - vv(ji,jj-1,jk,Kmm) ) / ( 2._wp * e2v(ji,jj) )
                               END_3D
@@ -139,7 +142,7 @@ CONTAINS
                               CALL iom_put( "vtrd_zdf", pvtrd )
                               !
                               !                                    ! wind stress trends
-                              ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) )
+                              ALLOCATE( z2dx(A2D(0)) , z2dy(A2D(0)) )
                               DO_2D( 0, 0, 0, 0 )
                                  z2dx(ji,jj) = ( utau_b(ji,jj) + utauU(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rho0 )
                                  z2dy(ji,jj) = ( vtau_b(ji,jj) + vtauV(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rho0 )
diff --git a/src/OCE/TRD/trdglo.F90 b/src/OCE/TRD/trdglo.F90
index a19770158..40c7b05b9 100644
--- a/src/OCE/TRD/trdglo.F90
+++ b/src/OCE/TRD/trdglo.F90
@@ -42,6 +42,7 @@ MODULE trdglo
    REAL(wp) ::   tvolv    ! volume of the whole ocean computed at v-points
    REAL(wp) ::   rpktrd   ! potential to kinetic energy conversion
    REAL(wp) ::   peke     ! conversion potential energy - kinetic energy trend
+   REAL(wp) ::   xcof
 
    !                     !!! domain averaged trends
    REAL(wp), DIMENSION(jptot_tra) ::   tmo, smo   ! temperature and salinity trends 
@@ -77,44 +78,47 @@ CONTAINS
       INTEGER ::   ji, jj, jk      ! dummy loop indices
       INTEGER ::   ikbu, ikbv      ! local integers
       REAL(wp)::   zvm, zvt, zvs, z1_2rho0   ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace 
+      REAL(wp), DIMENSION(A2D(0))  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace 
       !!----------------------------------------------------------------------
       !
       IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
          !
-      	SELECT CASE( ctype )
-      	!
-      	CASE( 'TRA' )          !==  Tracers (T & S)  ==!
-      	   DO_3D( 1, 1, 1, 1, 1, jpkm1 )   ! global sum of mask volume trend and trend*T (including interior mask)
-	            zvm = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
-	            zvt = ptrdx(ji,jj,jk) * zvm
-	            zvs = ptrdy(ji,jj,jk) * zvm
-   	         tmo(ktrd) = tmo(ktrd) + zvt   
-   	         smo(ktrd) = smo(ktrd) + zvs
-   	         t2 (ktrd) = t2(ktrd)  + zvt * ts(ji,jj,jk,jp_tem,Kmm)
-   	         s2 (ktrd) = s2(ktrd)  + zvs * ts(ji,jj,jk,jp_sal,Kmm)
-      	   END_3D
+     SELECT CASE( ctype )
+     !
+     CASE( 'TRA' )          !==  Tracers (T & S)  ==!
+        DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! global sum of mask volume trend and trend*T (including interior mask)
+           zvm = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
+           zvt = ptrdx(ji,jj,jk) * zvm
+           zvs = ptrdy(ji,jj,jk) * zvm
+           tmo(ktrd) = tmo(ktrd) + zvt   
+           smo(ktrd) = smo(ktrd) + zvs
+           t2 (ktrd) = t2(ktrd)  + zvt * ts(ji,jj,jk,jp_tem,Kmm)
+           s2 (ktrd) = s2(ktrd)  + zvs * ts(ji,jj,jk,jp_sal,Kmm)
+        END_3D
             !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
-            IF( ln_linssh .AND. ktrd == jptra_zad ) THEN  
-               tmo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
-               smo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
-               t2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
-               s2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
-            ENDIF
+        IF( ln_linssh .AND. ktrd == jptra_zad ) THEN  
+            DO_2D( 0, 0, 0, 0 )   ! global sum of mask volume trend and trend*T (including interior mask)
+               zvm = ww(ji,jj,1) * e1e2t(ji,jj) * tmask_i(ji,jj)
+               tmo(jptra_sad) = tmo(jptra_sad) + ts(ji,jj,1,jp_tem,Kmm) * zvm
+               smo(jptra_sad) = smo(jptra_sad) + ts(ji,jj,1,jp_sal,Kmm) * zvm
+               t2 (jptra_sad) = t2 (jptra_sad) + ts(ji,jj,1,jp_tem,Kmm) * ts(ji,jj,1,jp_tem,Kmm) * zvm
+               s2 (jptra_sad) = s2 (jptra_sad) + ts(ji,jj,1,jp_sal,Kmm) * ts(ji,jj,1,jp_sal,Kmm) * zvm
+            END_2D
+        ENDIF
             !
-            IF( ktrd == jptra_atf ) THEN     ! last trend (asselin time filter)
-               ! 
-               CALL glo_tra_wri( kt )             ! print the results in ocean.output
-               !                
-	            tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero)
-   	         smo(:) = 0._wp
-               t2 (:) = 0._wp
-               s2 (:) = 0._wp
-               !
-            ENDIF
+        IF( ktrd == jptra_atf ) THEN     ! last trend (asselin time filter)
+            ! 
+            CALL glo_tra_wri( kt )             ! print the results in ocean.output
+            !                
+            tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero)
+            smo(:) = 0._wp
+            t2 (:) = 0._wp
+            s2 (:) = 0._wp
+            !
+         ENDIF
             !
          CASE( 'DYN' )          !==  Momentum and KE  ==!        
-            DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
                   &                                     * e1e2u  (ji,jj) * e3u(ji,jj,jk,Kmm)
                zvs = ptrdy(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
@@ -126,7 +130,7 @@ CONTAINS
             !                 
             IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
                z1_2rho0 = 0.5_wp / rho0
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( 0, 0, 0, 0 )
                   zvt = ( utau_b(ji,jj) + utauU(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
                      &                                                      * z1_2rho0       * e1e2u(ji,jj)
                   zvs = ( vtau_b(ji,jj) + vtauV(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
@@ -184,8 +188,7 @@ CONTAINS
       INTEGER, INTENT(in) ::   Kmm  ! time level index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zcof         ! local scalar
-      REAL(wp), DIMENSION(jpi,jpj,jpk)  ::  zkx, zky, zkz, zkepe  
+      REAL(wp), DIMENSION(A2D(0),jpk)  ::  zkx, zky, zkz, zkepe  
       !!----------------------------------------------------------------------
 
       ! I. Momentum trends
@@ -196,24 +199,24 @@ CONTAINS
          ! I.1 Conversion potential energy - kinetic energy
          ! --------------------------------------------------
          ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
-         zkx  (:,:,:) = 0._wp
-         zky  (:,:,:) = 0._wp
-         zkz  (:,:,:) = 0._wp
-         zkepe(:,:,:) = 0._wp
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! loop from top to bottom
+            zkx  (ji,jj,jk) = 0._wp
+            zky  (ji,jj,jk) = 0._wp
+            zkz  (ji,jj,jk) = 0._wp
+            zkepe(ji,jj,jk) = 0._wp
+         END_3D
    
          CALL eos( ts(:,:,:,:,Kmm), rhd, rhop )       ! now potential density
 
-         zcof = 0.5_wp / rho0             ! Density flux at w-point
-         zkz(:,:,1) = 0._wp
-         DO jk = 2, jpk
-            zkz(:,:,jk) = zcof * e1e2t(:,:) * ww(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
-         END DO
+         zkz(A2D(0),1) = 0._wp
+         DO_3D( 0, 0, 0, 0, 1, jpk ) ! loop from top to bottom
+            zkz(ji,jj,jk) = xcof * e1e2t(ji,jj) * ww(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji,jj,jk-1) ) * tmask_i(ji,jj)
+         END_3D
          
-         zcof   = 0.5_wp / rho0           ! Density flux at u and v-points
-         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
-            zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)   &
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zkx(ji,jj,jk) = xcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)   &
                &                              *  uu(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
-            zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)   &
+            zky(ji,jj,jk) = xcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)   &
                &                              *  vv(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
          END_3D
          
@@ -227,10 +230,9 @@ CONTAINS
          ! I.2 Basin averaged kinetic energy trend
          ! ----------------------------------------
          peke = 0._wp
-         DO jk = 1, jpkm1
-            peke = peke + SUM( zkepe(:,:,jk) * gdept(:,:,jk,Kmm) * e1e2t(:,:)   &
-               &                               * e3t(:,:,jk,Kmm) )
-         END DO
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! loop from top to bottom
+            peke = peke + zkepe(ji,jj,jk) * gdept(ji,jj,jk,Kmm) * e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)
+         END_3D
          peke = grav * peke
 
          ! I.3 Sums over the global domain
@@ -509,11 +511,14 @@ CONTAINS
          WRITE(numout,*) '~~~~~~~~~~~~~'
       ENDIF
 
+      xcof = 0.5_wp / rho0
+
       ! Total volume at t-points:
       tvolt = 0._wp
-      DO jk = 1, jpkm1
-         tvolt = tvolt + SUM( e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) * tmask_i(:,:) )
-      END DO
+
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )    ! Density flux divergence at t-point
+         tvolt = tvolt + e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
+      END_3D
       CALL mpp_sum( 'trdglo', tvolt )   ! sum over the global domain
 
       IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
diff --git a/src/OCE/TRD/trdken.F90 b/src/OCE/TRD/trdken.F90
index 1920e3e72..d397a4d15 100644
--- a/src/OCE/TRD/trdken.F90
+++ b/src/OCE/TRD/trdken.F90
@@ -30,6 +30,10 @@ MODULE trdken
    IMPLICIT NONE
    PRIVATE
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+
    PUBLIC   trd_ken       ! called by trddyn module
    PUBLIC   trd_ken_init  ! called by trdini module
 
@@ -38,9 +42,6 @@ MODULE trdken
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   bu, bv   ! volume of u- and v-boxes
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   r1_bt    ! inverse of t-box volume
 
-   !! * Substitutions
-#  include "do_loop_substitute.h90"
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
    !! $Id: trdken.F90 15104 2021-07-07 14:36:00Z clem $
@@ -52,7 +53,7 @@ CONTAINS
       !!---------------------------------------------------------------------
       !!                  ***  FUNCTION trd_ken_alloc  ***
       !!---------------------------------------------------------------------
-      ALLOCATE( bu(jpi,jpj,jpk) , bv(jpi,jpj,jpk) , r1_bt(jpi,jpj,jpk) , STAT= trd_ken_alloc )
+      ALLOCATE( bu(A2D(0),jpk) , bv(A2D(0),jpk) , r1_bt(A2D(0),jpk) , STAT= trd_ken_alloc )
       !
       CALL mpp_sum ( 'trdken', trd_ken_alloc )
       IF( trd_ken_alloc /= 0 )   CALL ctl_stop( 'STOP', 'trd_ken_alloc: failed to allocate arrays' )
@@ -77,31 +78,32 @@ CONTAINS
       !
       !
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V masked trends
-      INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
-      INTEGER                   , INTENT(in   ) ::   kt             ! time step
-      INTEGER                   , INTENT(in   ) ::   Kmm            ! time level index
+      REAL(wp), DIMENSION(:,:,:)   , INTENT(inout) ::   putrd, pvtrd   ! U and V masked trends
+      INTEGER                      , INTENT(in   ) ::   ktrd           ! trend index
+      INTEGER                      , INTENT(in   ) ::   kt             ! time step
+      INTEGER                      , INTENT(in   ) ::   Kmm            ! time level index
       !
       INTEGER ::   ji, jj, jk       ! dummy loop indices
       INTEGER ::   ikbu  , ikbv     ! local integers
       INTEGER ::   ikbum1, ikbvm1   !   -       -
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   z2dx, z2dy, zke2d   ! 2D workspace 
-      REAL(wp), DIMENSION(jpi,jpj,jpk)      ::   zke                 ! 3D workspace 
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE    ::   zke2d   ! 2D workspace 
+      REAL(wp)                                 ::   z2dx, z2dy, z2dxm1, z2dym1   ! 2D workspace 
+      REAL(wp), DIMENSION(A2D(0),jpk)          ::   zke                 ! 3D workspace 
       !!----------------------------------------------------------------------
       !
       CALL lbc_lnk( 'trdken', putrd, 'U', -1.0_wp , pvtrd, 'V', -1.0_wp )      ! lateral boundary conditions
       !
       nkstp = kt
       DO jk = 1, jpkm1
-         bu   (:,:,jk) =    e1e2u(:,:) * e3u(:,:,jk,Kmm)
-         bv   (:,:,jk) =    e1e2v(:,:) * e3v(:,:,jk,Kmm)
-         r1_bt(:,:,jk) = r1_e1e2t(:,:) / e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+         DO_2D( 0, 0, 0, 0 )
+            bu   (ji,jj,jk) =    e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
+            bv   (ji,jj,jk) =    e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
+            r1_bt(ji,jj,jk) = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
+         END_2D
       END DO
       !
-      zke(:,:,jpk) = 0._wp
-      zke(1:nn_hls,:, : ) = 0._wp
-      zke(:,1:nn_hls, : ) = 0._wp
-      DO_3D( 0, nn_hls, 0, nn_hls, 1, jpkm1 )
+      zke(A2D(0),:) = 0._wp
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          zke(ji,jj,jk) = 0.5_wp * rho0 *( uu(ji  ,jj,jk,Kmm) * putrd(ji  ,jj,jk) * bu(ji  ,jj,jk)  &
             &                           + uu(ji-1,jj,jk,Kmm) * putrd(ji-1,jj,jk) * bu(ji-1,jj,jk)  &
             &                           + vv(ji,jj  ,jk,Kmm) * pvtrd(ji,jj  ,jk) * bv(ji,jj  ,jk)  &
@@ -118,37 +120,31 @@ CONTAINS
          CASE( jpdyn_ldf )   ;   CALL iom_put( "ketrd_ldf"   , zke )    ! lateral diffusion
          CASE( jpdyn_zdf )   ;   CALL iom_put( "ketrd_zdf"   , zke )    ! vertical diffusion 
          !                   !                                          ! wind stress trends
-            ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) , zke2d(jpi,jpj) )
-            DO_2D( 1, 0, 1, 0 )
-               z2dx(ji,jj) = uu(ji,jj,1,Kmm) * ( utau_b(ji,jj) + utauU(ji,jj) ) * e1e2u(ji,jj) * umask(ji,jj,1)
-               z2dy(ji,jj) = vv(ji,jj,1,Kmm) * ( vtau_b(ji,jj) + vtauV(ji,jj) ) * e1e2v(ji,jj) * vmask(ji,jj,1)
-            END_2D
-            DO_2D( 0, 0, 0, 0 )
-               zke2d(ji,jj) = r1_rho0 * 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
-                  &                                + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,1)
-            END_2D
-            CALL iom_put( "ketrd_tau"   , zke2d )  ! 
-            DEALLOCATE( z2dx , z2dy , zke2d )
-            !
+             ALLOCATE( zke2d(A2D(0)) ) ; zke2d(A2D(0)) = 0._wp 
+             DO_2D( 0, 0, 0, 0 )
+                z2dx = uu(ji,jj,1,Kmm) * ( utau_b(ji,jj) + utauU(ji,jj) ) * e1e2u(ji,jj) * umask(ji,jj,1)
+                z2dy = vv(ji,jj,1,Kmm) * ( vtau_b(ji,jj) + vtauV(ji,jj) ) * e1e2v(ji,jj) * vmask(ji,jj,1)
+                z2dxm1 = uu(ji-1,jj,1,Kmm) * ( utau_b(ji-1,jj) + utauU(ji-1,jj) ) * e1e2u(ji-1,jj) * umask(ji-1,jj,1)
+                z2dym1 = vv(ji,jj-1,1,Kmm) * ( vtau_b(ji,jj-1) + vtauV(ji,jj-1) ) * e1e2v(ji,jj-1) * vmask(ji,jj-1,1)
+                zke2d(ji,jj) = r1_rho0 * 0.5_wp * ( z2dx + z2dxm1 + z2dy + z2dym1 ) * r1_bt(ji,jj,1)
+             END_2D
+                                 CALL iom_put( "ketrd_tau"   , zke2d )  ! 
+                                 DEALLOCATE( zke2d )
          CASE( jpdyn_bfr )   ;   CALL iom_put( "ketrd_bfr"   , zke )    ! bottom friction (explicit case) 
 !!gm TO BE DONE properly
 !!gm only valid if ln_drgimp=F otherwise the bottom stress as to be recomputed at the end of the computation....
 !         IF(.NOT. ln_drgimp) THEN
-!            DO jj = 1, jpj    !   
-!               DO ji = 1, jpi
-!                  ikbu = mbku(ji,jj)         ! deepest ocean u- & v-levels
-!                  ikbv = mbkv(ji,jj)   
-!                  z2dx(ji,jj) = uu(ji,jj,ikbu,Kmm) * bfrua(ji,jj) * uu(ji,jj,ikbu,Kmm)
-!                  z2dy(ji,jj) = vv(ji,jj,ikbu,Kmm) * bfrva(ji,jj) * vv(ji,jj,ikbv,Kmm)
-!               END DO
-!            END DO
+!            zke2d(A2D(0)) = 0._wp
+!            DO_2D( 0, 0, 0, 0 )   
+!               ikbu   = mbku(ji,jj)         ! deepest ocean u- & v-levels
+!               ikbv   = mbkv(ji,jj)   
+!               z2dx   = uu(ji,jj,ikbu,Kmm) * bfrua(ji,jj) * uu(ji,jj,ikbu,Kmm)
+!               z2dy   = vv(ji,jj,ikbu,Kmm) * bfrva(ji,jj) * vv(ji,jj,ikbv,Kmm)
+!               z2dxm1 = uu(ji-1,jj,ikbu,Kmm) * bfrua(ji-1,jj) * uu(ji-1,jj,ikbu,Kmm)
+!               z2dym1 = vv(ji,jj-1,ikbu,Kmm) * bfrva(ji,jj-1) * vv(ji,jj-1,ikbv,Kmm)
+!               zke2d(ji,jj) = 0.5_wp * ( z2dx + z2dxm1 + z2dy + z2dym1 ) * r1_bt(ji,jj,1)
+!            END_2D
 !            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
-!            DO jj = 2, jpj
-!               DO ji = 2, jpi
-!                  zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
-!                     &                      + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,  BEURK!!!
-!               END DO
-!            END DO
 !                                    CALL iom_put( "ketrd_bfr"  , zke2d )   ! bottom friction (explicit case)
 !         ENDIF
 !!gm end
@@ -177,11 +173,11 @@ CONTAINS
         CASE( jpdyn_ken )   ;   ! kinetic energy
                     ! called in dynnxt.F90 before asselin time filter
                     ! with putrd=uu(Krhs) and pvtrd=vv(Krhs)
-                    zke(:,:,:) = 0.5_wp * zke(:,:,:)
+                    zke(A2D(0),:) = 0.5_wp * zke(A2D(0),:)
                     CALL iom_put( "KE", zke )
                     !
                     CALL ken_p2k( kt , zke, Kmm )
-                      CALL iom_put( "ketrd_convP2K", zke )     ! conversion -rau*g*w
+                    CALL iom_put( "ketrd_convP2K", zke )     ! conversion -rau*g*w
          !
       END SELECT
       !
@@ -205,22 +201,23 @@ CONTAINS
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   iku, ikv     ! local integers
       REAL(wp) ::   zcoef        ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zconv  ! 3D workspace
+      REAL(wp), DIMENSION(A2D(0),jpk) ::  zconv  ! 3D workspace
       !!----------------------------------------------------------------------
       !
       ! Local constant initialization 
       zcoef = - rho0 * grav * 0.5_wp      
       
       !  Surface value (also valid in partial step case)
-      zconv(:,:,1) = zcoef * ( 2._wp * rhd(:,:,1) ) * ww(:,:,1) * e3w(:,:,1,Kmm)
-
+      DO_2D( 0, 0, 0, 0 )
+         zconv(ji,jj,1) = zcoef * ( 2._wp * rhd(ji,jj,1) ) * ww(ji,jj,1) * e3w(ji,jj,1,Kmm)
+      END_2D
       ! interior value (2=<jk=<jpkm1)
-      DO jk = 2, jpk
-         zconv(:,:,jk) = zcoef * ( rhd(:,:,jk) + rhd(:,:,jk-1) ) * ww(:,:,jk) * e3w(:,:,jk,Kmm)
-      END DO
+      DO_3D( 0, 0, 0, 0 , 2, jpk )
+            zconv(ji,jj,jk) = zcoef * ( rhd(ji,jj,jk) + rhd(ji,jj,jk-1) ) * ww(ji,jj,jk) * e3w(ji,jj,jk,Kmm)
+      END_3D
 
       ! conv value on T-point
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          zcoef = 0.5_wp / e3t(ji,jj,jk,Kmm)
          pconv(ji,jj,jk) = zcoef * ( zconv(ji,jj,jk) + zconv(ji,jj,jk+1) ) * tmask(ji,jj,jk)
       END_3D
diff --git a/src/OCE/TRD/trdmxl.F90 b/src/OCE/TRD/trdmxl.F90
index bce0dd536..8edb0c54b 100644
--- a/src/OCE/TRD/trdmxl.F90
+++ b/src/OCE/TRD/trdmxl.F90
@@ -45,20 +45,7 @@ MODULE trdmxl
    PUBLIC   trd_mxl_init   ! routine called by opa.F90
    PUBLIC   trd_mxl_zint   ! routine called by tracers routines
 
-   INTEGER ::   nkstp       ! current time step 
-
-!!gm  to be moved from trdmxl_oce
-!   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   hml                ! ML depth (sum of e3t over nmln-1 levels) [m] 
-!   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   tml    , sml       ! now ML averaged T & S 
-!   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   tmlb_nf, smlb_nf   ! not filtered before ML averaged T & S
-!
-!
-!   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   hmlb, hmln         ! before, now, and after Mixed Layer depths [m]
-!   
-!   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   tb_mlb, tb_mln     ! before (not filtered) tracer averaged over before and now ML 
-!
-!   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   tn_mln             ! now tracer averaged over now ML
-!!gm end   
+   INTEGER ::   nkstp      ! current time step 
    
    CHARACTER (LEN=40) ::  clhstnam         ! name of the trends NetCDF file
    INTEGER ::   nh_t, nmoymltrd
@@ -111,44 +98,41 @@ CONTAINS
       IF ( kt /= nkstp ) THEN   !=  1st call at kt time step  =!
          !                      !==============================!
          nkstp = kt
-         
-         
          !                          !==  reset trend arrays to zero  ==!
-         tmltrd(:,:,:) = 0._wp    ;    smltrd(:,:,:) = 0._wp    
-         
-         
+         DO_3D( 0, 0, 0, 0, 1, jpktrd )
+            tmltrd(ji,jj,jk) = 0._wp
+            smltrd(ji,jj,jk) = 0._wp
+         END_3D
          !
-         wkx(:,:,:) = 0._wp         !==  now ML weights for vertical averaging  ==!
-         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpktrd )  ! initialize wkx with vertical scale factor in mixed-layer
+         !                          !==  now ML weights for vertical averaging  ==!
+         DO_3D( 0, 0, 0, 0, 1, jpktrd )  ! initialize wkx with vertical scale factor in mixed-layer
             IF( jk - kmxln(ji,jj) < 0 )   THEN
                wkx(ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
+            ELSE
+               wkx(ji,jj,jk) = 0._wp
             ENDIF
          END_3D
+         !
          hmxl(:,:) = 0._wp               ! NOW mixed-layer depth
-         DO jk = 1, jpktrd
-            hmxl(:,:) = hmxl(:,:) + wkx(:,:,jk)
-         END DO
-         DO jk = 1, jpktrd               ! integration weights
-            wkx(:,:,jk) = wkx(:,:,jk) / MAX( 1.e-20_wp, hmxl(:,:) ) * tmask(:,:,1)
-         END DO
-         
-         
+         DO_3D( 0, 0, 0, 0, 1, jpktrd )
+            hmxl(ji,jj)    = hmxl(ji,jj) + wkx(ji,jj,jk)
+            wkx (ji,jj,jk) = wkx (ji,jj,jk) / MAX( 1.e-20_wp, hmxl(ji,jj) ) * tmask(ji,jj,1)
+         END_3D
          !
          !                          !==  Vertically averaged T and S  ==!
          tml(:,:) = 0._wp   ;   sml(:,:) = 0._wp
-         DO jk = 1, jpktrd
-            tml(:,:) = tml(:,:) + wkx(:,:,jk) * ts(:,:,jk,jp_tem,Kmm)
-            sml(:,:) = sml(:,:) + wkx(:,:,jk) * ts(:,:,jk,jp_sal,Kmm)
-         END DO
+         DO_3D( 0, 0, 0, 0, 1, jpktrd )
+            tml(ji,jj) = tml(ji,jj) + wkx(ji,jj,jk) * ts(ji,jj,jk,jp_tem,Kmm)
+            sml(ji,jj) = sml(ji,jj) + wkx(ji,jj,jk) * ts(ji,jj,jk,jp_sal,Kmm)
+         END_3D
          !
       ENDIF
 
-
-
       ! mean now trends over the now ML 
-      tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + ptrdx(:,:,jk) * wkx(:,:,jk)   ! temperature
-      smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + ptrdy(:,:,jk) * wkx(:,:,jk)   ! salinity
- 
+      DO_2D( 0, 0, 0, 0 )
+         tmltrd(ji,jj,ktrd) = tmltrd(ji,jj,ktrd) + ptrdx(ji,jj,jk) * wkx(ji,jj,jk)   ! temperature
+         smltrd(ji,jj,ktrd) = smltrd(ji,jj,ktrd) + ptrdy(ji,jj,jk) * wkx(ji,jj,jk)   ! salinity
+      END_2D
 
  
 !!gm to be put juste before the output !
@@ -249,11 +233,11 @@ CONTAINS
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT( in ) ::   ktrd       ! ocean trend index
       CHARACTER(len=2)                , INTENT( in ) ::   ctype      ! 2D surface/bottom or 3D interior physics
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   pttrdmxl   ! temperature trend 
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   pstrdmxl   ! salinity trend 
+      REAL(wp), DIMENSION(A2D(0),jpk), INTENT( in ) ::   pttrdmxl   ! temperature trend 
+      REAL(wp), DIMENSION(A2D(0),jpk), INTENT( in ) ::   pstrdmxl   ! salinity trend 
       !
       INTEGER ::   ji, jj, jk, isum
-      REAL(wp), DIMENSION(jpi,jpj)  :: zvlmsk 
+      REAL(wp), DIMENSION(A2D(0))  :: zvlmsk 
       !!----------------------------------------------------------------------
 
       ! I. Definition of control surface and associated fields
@@ -268,11 +252,13 @@ CONTAINS
 
          
          ! ... Set nmxl(ji,jj) = index of first T point below control surf. or outside mixed-layer
-         IF( nn_ctls == 0 ) THEN       ! * control surface = mixed-layer with density criterion 
-            nmxl(:,:) = nmln(:,:)    ! array nmln computed in zdfmxl.F90
-         ELSEIF( nn_ctls == 1 ) THEN   ! * control surface = read index from file 
+         IF( nn_ctls == 0 ) THEN         ! * control surface = mixed-layer with density criterion 
+            DO_2D( 0, 0, 0, 0 )
+               nmxl(ji,jj) = nmln(ji,jj) !   array nmln computed in zdfmxl.F90
+            END_2D
+         ELSEIF( nn_ctls == 1 ) THEN     ! * control surface = read index from file 
             nmxl(:,:) = nbol(:,:)
-         ELSEIF( nn_ctls >= 2 ) THEN   ! * control surface = model level
+         ELSEIF( nn_ctls >= 2 ) THEN     ! * control surface = model level
             nn_ctls = MIN( nn_ctls, jpktrd - 1 )
             nmxl(:,:) = nn_ctls + 1
          ENDIF
@@ -335,9 +321,9 @@ CONTAINS
       REAL(wp) :: zavt, zfn, zfn2
       !                                              ! z(ts)mltot : dT/dt over the anlysis window (including Asselin)
       !                                              ! z(ts)mlres : residual = dh/dt entrainment term
-      REAL(wp), DIMENSION(jpi,jpj  )   ::  ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf
-      REAL(wp), DIMENSION(jpi,jpj  )   ::  ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2  
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  ztmltrd2, zsmltrd2   ! only needed for mean diagnostics
+      REAL(wp), DIMENSION(A2D(0)  )   ::  ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf
+      REAL(wp), DIMENSION(A2D(0)  )   ::  ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2  
+      REAL(wp), DIMENSION(A2D(0),jpk) ::  ztmltrd2, zsmltrd2   ! only needed for mean diagnostics
       !!----------------------------------------------------------------------
   
       ! ======================================================================
@@ -446,12 +432,7 @@ CONTAINS
       itmod = kt - nit000 + 1
 
       MODULO_NTRD : IF( MOD( itmod, nn_trd ) == 0 ) THEN        ! nitend MUST be multiple of nn_trd
-         !
-         ztmltot (:,:) = 0.e0   ;   zsmltot (:,:) = 0.e0   ! reset arrays to zero
-         ztmlres (:,:) = 0.e0   ;   zsmlres (:,:) = 0.e0
-         ztmltot2(:,:) = 0.e0   ;   zsmltot2(:,:) = 0.e0
-         ztmlres2(:,:) = 0.e0   ;   zsmlres2(:,:) = 0.e0
-      
+         !      
          zfn  = REAL( nmoymltrd, wp )   ;   zfn2 = zfn * zfn
          
          ! III.1 Prepare fields for output ("instantaneous" diagnostics) 
@@ -469,25 +450,18 @@ CONTAINS
          ztmlatf(:,:) = tmlatfm(:,:) - tmlatfn(:,:) + tmlatfb(:,:)
          zsmlatf(:,:) = smlatfm(:,:) - smlatfn(:,:) + smlatfb(:,:)
          
-         !-- Lateral boundary conditions
-         !         ... temperature ...                    ... salinity ...
-         CALL lbc_lnk( 'trdmxl', ztmltot , 'T', 1.0_wp, zsmltot , 'T', 1.0_wp, &
-            &                    ztmlres , 'T', 1.0_wp, zsmlres , 'T', 1.0_wp, &
-            &                    ztmlatf , 'T', 1.0_wp, zsmlatf , 'T', 1.0_wp )
-
-
          ! III.2 Prepare fields for output ("mean" diagnostics) 
          ! ----------------------------------------------------
          
          !-- Update the ML depth time sum (to build the Leap-Frog time mean)
-         hmxl_sum(:,:) = hmxlbn(:,:) + 2 * ( hmxl_sum(:,:) - hmxl(:,:) ) + hmxl(:,:)
+         hmxl_sum(:,:) = hmxlbn(:,:) + 2._wp * ( hmxl_sum(:,:) - hmxl(:,:) ) + hmxl(:,:)
 
          !-- Compute temperature total trends
-         tml_sum (:,:) = tmlbn(:,:) + 2 * ( tml_sum(:,:) - tml(:,:) ) + tml(:,:)
+         tml_sum (:,:) = tmlbn(:,:) + 2._wp * ( tml_sum(:,:) - tml(:,:) ) + tml(:,:)
          ztmltot2(:,:) = ( tml_sum(:,:) - tml_sumb(:,:) ) / p2dt    ! now in degC/s
          
          !-- Compute salinity total trends
-         sml_sum (:,:) = smlbn(:,:) + 2 * ( sml_sum(:,:) - sml(:,:) ) + sml(:,:)
+         sml_sum (:,:) = smlbn(:,:) + 2._wp * ( sml_sum(:,:) - sml(:,:) ) + sml(:,:)
          zsmltot2(:,:) = ( sml_sum(:,:) - sml_sumb(:,:) ) / p2dt    ! now in psu/s
          
          !-- Compute temperature residuals
@@ -495,7 +469,7 @@ CONTAINS
             ztmltrd2(:,:,jl) = tmltrd_csum_ub(:,:,jl) + tmltrd_csum_ln(:,:,jl)
          END DO
 
-         ztmltrdm2(:,:) = 0.e0
+         ztmltrdm2(:,:) = 0._wp
          DO jl = 1, jpltrd
             ztmltrdm2(:,:) = ztmltrdm2(:,:) + ztmltrd2(:,:,jl)
          END DO
@@ -508,7 +482,7 @@ CONTAINS
             zsmltrd2(:,:,jl) = smltrd_csum_ub(:,:,jl) + smltrd_csum_ln(:,:,jl)
          END DO
 
-         zsmltrdm2(:,:) = 0.
+         zsmltrdm2(:,:) = 0._wp
          DO jl = 1, jpltrd
             zsmltrdm2(:,:) = zsmltrdm2(:,:) + zsmltrd2(:,:,jl)
          END DO
@@ -519,13 +493,6 @@ CONTAINS
          !-- Diagnose Asselin trend over the analysis window
          ztmlatf2(:,:) = ztmltrd2(:,:,jpmxl_atf) - tmltrd_sum(:,:,jpmxl_atf) + tmltrd_atf_sumb(:,:)
          zsmlatf2(:,:) = zsmltrd2(:,:,jpmxl_atf) - smltrd_sum(:,:,jpmxl_atf) + smltrd_atf_sumb(:,:)
-
-         !-- Lateral boundary conditions
-         !         ... temperature ...                    ... salinity ...
-         CALL lbc_lnk( 'trdmxl', ztmltot2, 'T', 1.0_wp, zsmltot2, 'T', 1.0_wp, &
-            &                    ztmlres2, 'T', 1.0_wp, zsmlres2, 'T', 1.0_wp )
-         !
-         CALL lbc_lnk( 'trdmxl', ztmltrd2(:,:,:), 'T', 1.0_wp, zsmltrd2(:,:,:), 'T', 1.0_wp ) ! /  in the NetCDF trends file
          
          ! III.3 Time evolution array swap
          ! -------------------------------
@@ -622,6 +589,8 @@ CONTAINS
       ! ======================================================================
 
       !-- Write the trends for T/S instantaneous diagnostics 
+
+      ! clem => these fields do not exist in field_def
       
       IF( ln_trdmxl_instant ) THEN           
 
diff --git a/src/OCE/TRD/trdmxl_oce.F90 b/src/OCE/TRD/trdmxl_oce.F90
index c2d4d4be9..b1266b5ce 100644
--- a/src/OCE/TRD/trdmxl_oce.F90
+++ b/src/OCE/TRD/trdmxl_oce.F90
@@ -80,6 +80,8 @@ MODULE trdmxl_oce
       smltrd_csum_ln,               & !:    ( idem for salinity )
       smltrd_csum_ub                  !: 
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
    !! $Id: trdmxl_oce.F90 10425 2018-12-19 21:54:16Z smasson $ 
@@ -101,29 +103,29 @@ CONTAINS
 
      ierr(:) = 0
 
-     ALLOCATE( nmxl (jpi,jpj)    , nbol (jpi,jpj),    &
-        &      wkx  (jpi,jpj,jpk), hmxl (jpi,jpj),    & 
-        &      tml  (jpi,jpj)    , sml  (jpi,jpj),    & 
-        &      tmlb (jpi,jpj)    , smlb (jpi,jpj),    &
-        &      tmlbb(jpi,jpj)    , smlbb(jpi,jpj), STAT = ierr(1) )
-
-     ALLOCATE( tmlbn(jpi,jpj)  , smlbn(jpi,jpj),   &
-        &      tmltrdm(jpi,jpj), smltrdm(jpi,jpj), &
-        &      tml_sum(jpi,jpj), tml_sumb(jpi,jpj),&
-        &      tmltrd_atf_sumb(jpi,jpj)           , STAT=ierr(2) )
-
-     ALLOCATE( sml_sum(jpi,jpj), sml_sumb(jpi,jpj), &
-        &      smltrd_atf_sumb(jpi,jpj),            &
-        &      hmxl_sum(jpi,jpj), hmxlbn(jpi,jpj),  &
-        &      tmlatfb(jpi,jpj), tmlatfn(jpi,jpj), STAT = ierr(3) )
-
-     ALLOCATE( smlatfb(jpi,jpj), smlatfn(jpi,jpj), & 
-        &      tmlatfm(jpi,jpj), smlatfm(jpi,jpj), &
-        &      tmltrd(jpi,jpj,jpltrd),   smltrd(jpi,jpj,jpltrd), STAT=ierr(4))
-
-     ALLOCATE( tmltrd_sum(jpi,jpj,jpltrd),tmltrd_csum_ln(jpi,jpj,jpltrd),      &
-        &      tmltrd_csum_ub(jpi,jpj,jpltrd), smltrd_sum(jpi,jpj,jpltrd),     &
-        &      smltrd_csum_ln(jpi,jpj,jpltrd), smltrd_csum_ub(jpi,jpj,jpltrd), STAT=ierr(5) )
+     ALLOCATE( nmxl (A2D(0))    , nbol (A2D(0)),    &
+        &      wkx  (A2D(0),jpk), hmxl (A2D(0)),    & 
+        &      tml  (A2D(0))    , sml  (A2D(0)),    & 
+        &      tmlb (A2D(0))    , smlb (A2D(0)),    &
+        &      tmlbb(A2D(0))    , smlbb(A2D(0)), STAT = ierr(1) )
+
+     ALLOCATE( tmlbn(A2D(0))  , smlbn(A2D(0)),   &
+        &      tmltrdm(A2D(0)), smltrdm(A2D(0)), &
+        &      tml_sum(A2D(0)), tml_sumb(A2D(0)),&
+        &      tmltrd_atf_sumb(A2D(0))           , STAT=ierr(2) )
+
+     ALLOCATE( sml_sum(A2D(0)), sml_sumb(A2D(0)), &
+        &      smltrd_atf_sumb(A2D(0)),           &
+        &      hmxl_sum(A2D(0)), hmxlbn(A2D(0)),  &
+        &      tmlatfb(A2D(0)), tmlatfn(A2D(0)), STAT = ierr(3) )
+
+     ALLOCATE( smlatfb(A2D(0)), smlatfn(A2D(0)), & 
+        &      tmlatfm(A2D(0)), smlatfm(A2D(0)), &
+        &      tmltrd(A2D(0),jpltrd),   smltrd(A2D(0),jpltrd), STAT=ierr(4))
+
+     ALLOCATE( tmltrd_sum(A2D(0),jpltrd),tmltrd_csum_ln(A2D(0),jpltrd),      &
+        &      tmltrd_csum_ub(A2D(0),jpltrd), smltrd_sum(A2D(0),jpltrd),     &
+        &      smltrd_csum_ln(A2D(0),jpltrd), smltrd_csum_ub(A2D(0),jpltrd), STAT=ierr(5) )
       !
       trdmxl_oce_alloc = MAXVAL( ierr )
       CALL mpp_sum ( 'trdmxl_oce', trdmxl_oce_alloc )
diff --git a/src/OCE/TRD/trdpen.F90 b/src/OCE/TRD/trdpen.F90
index 0c356f841..1a1fdec15 100644
--- a/src/OCE/TRD/trdpen.F90
+++ b/src/OCE/TRD/trdpen.F90
@@ -35,6 +35,7 @@ MODULE trdpen
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   rab_pe   ! partial derivatives of PE anomaly with respect to T and S
 
    !! * Substitutions
+#  include "do_loop_substitute.h90"
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -48,7 +49,7 @@ CONTAINS
       !!---------------------------------------------------------------------
       !!                  ***  FUNCTION trd_tra_alloc  ***
       !!---------------------------------------------------------------------
-      ALLOCATE( rab_pe(jpi,jpj,jpk,jpts) , STAT= trd_pen_alloc )
+      ALLOCATE( rab_pe(A2D(0),jpk,jpts) , STAT= trd_pen_alloc )
       !
       CALL mpp_sum ( 'trdpen', trd_pen_alloc )
       IF( trd_pen_alloc /= 0 )   CALL ctl_stop( 'STOP',  'trd_pen_alloc: failed to allocate arrays'  )
@@ -69,37 +70,40 @@ CONTAINS
       INTEGER                   , INTENT(in) ::   Kmm            ! time level index
       REAL(wp)                  , INTENT(in) ::   pdt            ! time step [s]
       !
-      INTEGER ::   jk                                            ! dummy loop indices
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)      ::   z2d            ! 2D workspace 
-      REAL(wp), DIMENSION(jpi,jpj,jpk)           ::   zpe            ! 3D workspace 
+      INTEGER ::   ji, jj, jk                                            ! dummy loop indices
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)  ::   z2d            ! 2D workspace 
+      REAL(wp), DIMENSION(A2D(0),jpk)        ::   zpe            ! 3D workspace 
       !!----------------------------------------------------------------------
       !
-      zpe(:,:,:) = 0._wp
+      zpe(A2D(0),:) = 0._wp
       !
       IF( kt /= nkstp ) THEN     ! full eos: set partial derivatives at the 1st call of kt time step
          nkstp = kt
-         CALL eos_pen( ts(:,:,:,:,Kmm), rab_PE, zpe, Kmm )
+         CALL eos_pen( ts(A2D(0),:,:,Kmm), rab_pe, zpe, Kmm )
          CALL iom_put( "alphaPE", rab_pe(:,:,:,jp_tem) )
          CALL iom_put( "betaPE" , rab_pe(:,:,:,jp_sal) )
          CALL iom_put( "PEanom" , zpe )
       ENDIF
       !
-      zpe(:,:,jpk) = 0._wp
-      DO jk = 1, jpkm1
-         zpe(:,:,jk) = ( - ( rab_n(:,:,jk,jp_tem) + rab_pe(:,:,jk,jp_tem) ) * ptrdx(:,:,jk)   &
-            &            + ( rab_n(:,:,jk,jp_sal) + rab_pe(:,:,jk,jp_sal) ) * ptrdy(:,:,jk)  )
-      END DO
+      zpe(A2D(0),jpk) = 0._wp
+      !
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         zpe(ji,jj,jk) = ( - ( rab_n(ji,jj,jk,jp_tem) + rab_pe(ji,jj,jk,jp_tem) ) * ptrdx(ji,jj,jk)   &
+            &              + ( rab_n(ji,jj,jk,jp_sal) + rab_pe(ji,jj,jk,jp_sal) ) * ptrdy(ji,jj,jk)  )
+      END_3D
 
       SELECT CASE ( ktrd )
       CASE ( jptra_xad  )   ;   CALL iom_put( "petrd_xad", zpe )   ! zonal    advection
       CASE ( jptra_yad  )   ;   CALL iom_put( "petrd_yad", zpe )   ! merid.   advection
       CASE ( jptra_zad  )   ;   CALL iom_put( "petrd_zad", zpe )   ! vertical advection
                                 IF( ln_linssh ) THEN                   ! cst volume : adv flux through z=0 surface
-                                   ALLOCATE( z2d(jpi,jpj) )
-                                   z2d(:,:) = ww(:,:,1) * ( &
-                                     &   - ( rab_n(:,:,1,jp_tem) + rab_pe(:,:,1,jp_tem) ) * ts(:,:,1,jp_tem,Kmm)    &
-                                     &   + ( rab_n(:,:,1,jp_sal) + rab_pe(:,:,1,jp_sal) ) * ts(:,:,1,jp_sal,Kmm)    &
-                                     & ) / e3t(:,:,1,Kmm)
+                                   ALLOCATE( z2d(A2D(0)) )
+                                   DO_2D( 0, 0, 0, 0 )
+                                      z2d(ji,jj) = ww(ji,jj,1) * ( &
+                                        &   - ( rab_n(ji,jj,1,jp_tem) + rab_pe(ji,jj,1,jp_tem) ) * ts(ji,jj,1,jp_tem,Kmm)    &
+                                        &   + ( rab_n(ji,jj,1,jp_sal) + rab_pe(ji,jj,1,jp_sal) ) * ts(ji,jj,1,jp_sal,Kmm)    &
+                                        & ) / e3t(ji,jj,1,Kmm)
+                                   END_2D
                                    CALL iom_put( "petrd_sad" , z2d )
                                    DEALLOCATE( z2d )
                                 ENDIF
diff --git a/src/OCE/TRD/trdtra.F90 b/src/OCE/TRD/trdtra.F90
index a33778632..6c978ddd1 100644
--- a/src/OCE/TRD/trdtra.F90
+++ b/src/OCE/TRD/trdtra.F90
@@ -53,7 +53,7 @@ CONTAINS
       !!---------------------------------------------------------------------
       !!                  ***  FUNCTION trd_tra_alloc  ***
       !!---------------------------------------------------------------------
-      ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , avt_evd(jpi,jpj,jpk), STAT= trd_tra_alloc )
+      ALLOCATE( trdtx(A2D(0),jpk), trdty(A2D(0),jpk), trdt(A2D(0),jpk), avt_evd(A2D(0),jpk), STAT= trd_tra_alloc )
       !
       CALL mpp_sum ( 'trdtra', trd_tra_alloc )
       IF( trd_tra_alloc /= 0 )   CALL ctl_stop( 'STOP', 'trd_tra_alloc: failed to allocate arrays' )
@@ -73,24 +73,25 @@ CONTAINS
       !!              - 'TRA' case : regroup T & S trends
       !!              - send the trends to trd_tra_mng (trdtrc) for further processing
       !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in)           ::   kt      ! time step
-      CHARACTER(len=3)                , INTENT(in)           ::   ctype   ! tracers trends type 'TRA'/'TRC'
-      INTEGER                         , INTENT(in)           ::   ktra    ! tracer index
-      INTEGER                         , INTENT(in)           ::   ktrd    ! tracer trend index
-      INTEGER                         , INTENT(in)           ::   Kmm, Krhs ! time level indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)           ::   ptrd    ! tracer trend  or flux
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pu      ! now velocity 
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   ptra    ! now tracer variable
+      INTEGER                   , INTENT(in)           ::   kt      ! time step
+      CHARACTER(len=3)          , INTENT(in)           ::   ctype   ! tracers trends type 'TRA'/'TRC'
+      INTEGER                   , INTENT(in)           ::   ktra    ! tracer index
+      INTEGER                   , INTENT(in)           ::   ktrd    ! tracer trend index
+      INTEGER                   , INTENT(in)           ::   Kmm, Krhs ! time level indices
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)           ::   ptrd    ! tracer trend  or flux
+      REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL ::   pu      ! now velocity 
+      REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL ::   ptra    ! now tracer variable
       !
-      INTEGER ::   jk    ! loop indices
-      INTEGER ::   i01   ! 0 or 1
-      REAL(wp),        DIMENSION(jpi,jpj,jpk) ::   ztrds             ! 3D workspace
+      INTEGER  ::   ji,jj,jk    ! loop indices
+      INTEGER  ::   i01   ! 0 or 1
+      REAL(wp) ::  z1e3w
+      REAL(wp),              DIMENSION(A2D(0),jpk) ::   ztrds             ! 3D workspace
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zwt, zws, ztrdt   ! 3D workspace
       !!----------------------------------------------------------------------
       !      
       IF( .NOT. ALLOCATED( trdtx ) ) THEN      ! allocate trdtra arrays
          IF( trd_tra_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trd_tra : unable to allocate arrays' )
-         avt_evd(:,:,:) = 0._wp
+         avt_evd(A2D(0),:) = 0._wp
       ENDIF
       !
       i01 = COUNT( (/ PRESENT(pu) .OR. ( ktrd /= jptra_xad .AND. ktrd /= jptra_yad .AND. ktrd /= jptra_zad ) /) )
@@ -103,13 +104,13 @@ CONTAINS
          CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd, pu, ptra, 'Y', trdty, Kmm ) 
          CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd, pu, ptra, 'Z', trdt, Kmm )
          CASE( jptra_bbc,    &        ! qsr, bbc: on temperature only, send to trd_tra_mng
-            &  jptra_qsr )   ;   trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
-                                 ztrds(:,:,:) = 0._wp
+            &  jptra_qsr )   ;   trdt(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:)
+                                 ztrds(A2D(0),:) = 0._wp
                                  CALL trd_tra_mng( trdt, ztrds, ktrd, kt, Kmm )
  !!gm Gurvan, verify the jptra_evd trend please !
-         CASE( jptra_evd )   ;   avt_evd(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
+         CASE( jptra_evd )   ;   avt_evd(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:)
          CASE DEFAULT                 ! other trends: masked trends
-            trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)              ! mask & store
+            trdt(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:)              ! mask & store
          END SELECT
          !
       ENDIF
@@ -127,44 +128,44 @@ CONTAINS
                                   CALL trd_tra_mng( trdt , ztrds, ktrd, kt, Kmm   )
          CASE( jptra_zdfp )           ! diagnose the "PURE" Kz trend (here: just before the swap)
             !                         ! iso-neutral diffusion case otherwise jptra_zdf is "PURE"
-            ALLOCATE( zwt(jpi,jpj,jpk), zws(jpi,jpj,jpk), ztrdt(jpi,jpj,jpk) )
+            ALLOCATE( zwt(A2D(0),jpk), zws(A2D(0),jpk), ztrdt(A2D(0),jpk) )
             !
-            zwt(:,:, 1 ) = 0._wp   ;   zws(:,:, 1 ) = 0._wp            ! vertical diffusive fluxes
-            zwt(:,:,jpk) = 0._wp   ;   zws(:,:,jpk) = 0._wp
-            DO jk = 2, jpk
-               zwt(:,:,jk) = avt(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) )   &
-                  &        / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
-               zws(:,:,jk) = avs(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) )   &
-                  &        / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
-            END DO
+            zwt(A2D(0), 1 ) = 0._wp   ;   zws(A2D(0), 1 ) = 0._wp            ! vertical diffusive fluxes
+            zwt(A2D(0),jpk) = 0._wp   ;   zws(A2D(0),jpk) = 0._wp
+            DO_3D( 0, 0, 0, 0, 2, jpk )
+               z1e3w = 1._wp / ( e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) )
+               zwt(ji,jj,jk) = avt(ji,jj,jk) * ( ts(ji,jj,jk-1,jp_tem,Krhs) - ts(ji,jj,jk,jp_tem,Krhs) ) * z1e3w
+               zws(ji,jj,jk) = avs(ji,jj,jk) * ( ts(ji,jj,jk-1,jp_sal,Krhs) - ts(ji,jj,jk,jp_sal,Krhs) ) * z1e3w
+            END_3D
             !
-            ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp
-            DO jk = 1, jpkm1
-               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
-               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t(:,:,jk,Kmm) 
-            END DO
+            ztrdt(A2D(0),jpk) = 0._wp   ;   ztrds(A2D(0),jpk) = 0._wp
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+               z1e3w = 1._wp / e3w(ji,jj,jk,Kmm) 
+               ztrdt(ji,jj,jk) = ( zwt(ji,jj,jk) - zwt(ji,jj,jk+1) ) * z1e3w
+               ztrds(ji,jj,jk) = ( zws(ji,jj,jk) - zws(ji,jj,jk+1) ) * z1e3w
+            END_3D
             CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt, Kmm )  
             !
             !                         ! Also calculate EVD trend at this point. 
-            zwt(:,:,:) = 0._wp   ;   zws(:,:,:) = 0._wp            ! vertical diffusive fluxes
-            DO jk = 2, jpk
-               zwt(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) )   &
-                  &            / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
-               zws(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) )   &
-                  &            / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
-            END DO
+            zwt(A2D(0), :) = 0._wp   ;   zws(A2D(0), :) = 0._wp            ! vertical diffusive fluxes
+            DO_3D( 0, 0, 0, 0, 2, jpk )
+               z1e3w = 1._wp / ( e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) )
+               zwt(ji,jj,jk) = avt_evd(ji,jj,jk) * ( ts(ji,jj,jk-1,jp_tem,Krhs) - ts(ji,jj,jk,jp_tem,Krhs) ) * z1e3w
+               zws(ji,jj,jk) = avt_evd(ji,jj,jk) * ( ts(ji,jj,jk-1,jp_sal,Krhs) - ts(ji,jj,jk,jp_sal,Krhs) ) * z1e3w
+            END_3D
             !
-            ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp
-            DO jk = 1, jpkm1
-               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
-               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t(:,:,jk,Kmm) 
-            END DO
+            ztrdt(A2D(0),jpk) = 0._wp   ;   ztrds(A2D(0),jpk) = 0._wp
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+               z1e3w = 1._wp / e3w(ji,jj,jk,Kmm) 
+               ztrdt(ji,jj,jk) = ( zwt(ji,jj,jk) - zwt(ji,jj,jk+1) ) * z1e3w
+               ztrds(ji,jj,jk) = ( zws(ji,jj,jk) - zws(ji,jj,jk+1) ) * z1e3w
+            END_3D
             CALL trd_tra_mng( ztrdt, ztrds, jptra_evd, kt, Kmm )  
             !
             DEALLOCATE( zwt, zws, ztrdt )
             !
          CASE DEFAULT                 ! other trends: mask and send T & S trends to trd_tra_mng
-            ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
+            ztrds(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:)
             CALL trd_tra_mng( trdt, ztrds, ktrd, kt, Kmm )  
          END SELECT
       ENDIF
@@ -177,7 +178,7 @@ CONTAINS
          CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd , pu , ptra, 'Y', ztrds, Kmm ) 
          CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd , pu , ptra, 'Z', ztrds, Kmm ) 
          CASE DEFAULT                 ! other trends: just masked 
-                                 ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
+                                 ztrds(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:)
          END SELECT
          !                            ! send trend to trd_trc
          CALL trd_trc( ztrds, ktra, ktrd, kt, Kmm ) 
@@ -199,12 +200,12 @@ CONTAINS
       !!       k-advective trends = -un. di+1[T] = -( dk+1[fi] - tn dk+1[un] )
       !!                where fi is the incoming advective flux.
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pf      ! advective flux in one direction
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu      ! now velocity   in one direction
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt      ! now or before tracer 
-      CHARACTER(len=1)                , INTENT(in   ) ::   cdir    ! X/Y/Z direction
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   ptrd    ! advective trend in one direction
-      INTEGER,  INTENT(in)                            ::   Kmm     ! time level index
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pf      ! advective flux in one direction
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pu      ! now velocity   in one direction
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pt      ! now or before tracer 
+      CHARACTER(len=1)          , INTENT(in   ) ::   cdir    ! X/Y/Z direction
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   ptrd    ! advective trend in one direction
+      INTEGER,  INTENT(in)                      ::   Kmm     ! time level index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   ii, ij, ik   ! index shift as function of the direction
@@ -217,9 +218,7 @@ CONTAINS
       END SELECT
       !
       !                               ! set to zero uncomputed values
-      ptrd(jpi,:,:) = 0._wp   ;   ptrd(1,:,:) = 0._wp
-      ptrd(:,jpj,:) = 0._wp   ;   ptrd(:,1,:) = 0._wp
-      ptrd(:,:,jpk) = 0._wp
+      ptrd(:,:,:) = 0._wp 
       !
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! advective trend
          ptrd(ji,jj,jk) = - (     pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)                        &
@@ -248,7 +247,7 @@ CONTAINS
       !                   ! 3D output of tracers trends using IOM interface
       IF( ln_tra_trd )   CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt, Kmm )
 
-      !                   ! Integral Constraints Properties for tracers trends                                       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                   ! Integral Constraints Properties for tracers trends           
       IF( ln_glo_trd )   CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt, Kmm )
 
       !                   ! Potential ENergy trends
@@ -305,8 +304,8 @@ CONTAINS
       INTEGER                   , INTENT(in   ) ::   kt      ! time step
       INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
       !!
-      INTEGER ::   ji, jj, jk   ! dummy loop indices
-      INTEGER ::   ikbu, ikbv   ! local integers
+      INTEGER  :: ji, jj
+      REAL(wp) :: z1e3
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace 
       !!----------------------------------------------------------------------
       !
@@ -329,9 +328,12 @@ CONTAINS
          CASE( jptra_zad  )   ;   CALL iom_put( "ttrd_zad"  , ptrdx )        ! z- vertical   advection
                                   CALL iom_put( "strd_zad"  , ptrdy )
                                   IF( ln_linssh ) THEN                   ! cst volume : adv flux through z=0 surface
-                                     ALLOCATE( z2dx(jpi,jpj), z2dy(jpi,jpj) )
-                                     z2dx(:,:) = ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) / e3t(:,:,1,Kmm)
-                                     z2dy(:,:) = ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) / e3t(:,:,1,Kmm)
+                                     ALLOCATE( z2dx(A2D(0)), z2dy(A2D(0)) )
+                                     DO_2D( 0, 0, 0, 0 )
+                                        z1e3 = ww(ji,jj,1) / e3t(ji,jj,1,Kmm)
+                                        z2dx(ji,jj) = ts(ji,jj,1,jp_tem,Kmm) * z1e3
+                                        z2dy(ji,jj) = ts(ji,jj,1,jp_sal,Kmm) * z1e3
+                                     END_2D
                                      CALL iom_put( "ttrd_sad", z2dx )
                                      CALL iom_put( "strd_sad", z2dy )
                                      DEALLOCATE( z2dx, z2dy )
diff --git a/src/OCE/TRD/trdvor.F90 b/src/OCE/TRD/trdvor.F90
index 8515d666d..b01062697 100644
--- a/src/OCE/TRD/trdvor.F90
+++ b/src/OCE/TRD/trdvor.F90
@@ -68,9 +68,9 @@ CONTAINS
       !!----------------------------------------------------------------------------
       !!                  ***  ROUTINE trd_vor_alloc  ***
       !!----------------------------------------------------------------------------
-      ALLOCATE( vor_avr   (jpi,jpj) , vor_avrb(jpi,jpj) , vor_avrbb (jpi,jpj) ,   &
-         &      vor_avrbn (jpi,jpj) , rotot   (jpi,jpj) , vor_avrtot(jpi,jpj) ,   &
-         &      vor_avrres(jpi,jpj) , vortrd  (jpi,jpj,jpltot_vor) ,              &
+      ALLOCATE( vor_avr   (A2D(0)) , vor_avrb(A2D(0)) , vor_avrbb (A2D(0)) ,   &
+         &      vor_avrbn (A2D(0)) , rotot   (A2D(0)) , vor_avrtot(A2D(0)) ,   &
+         &      vor_avrres(A2D(0)) , vortrd  (A2D(0),jpltot_vor) ,              &
          &      ndexvor1  (jpi*jpj)                                ,   STAT= trd_vor_alloc )
          !
       CALL mpp_sum ( 'trdvor', trd_vor_alloc )
@@ -375,10 +375,6 @@ CONTAINS
          zmean = 1._wp / REAL( nmoydpvor, wp )
          vor_avrres(:,:) = vor_avrtot(:,:) - rotot(:,:) / zmean
 
-         ! Boundary conditions
-         CALL lbc_lnk( 'trdvor', vor_avrtot, 'F', 1.0_wp , vor_avrres, 'F', 1.0_wp )
-
-
          ! III.3 time evolution array swap
          ! ------------------------------
          vor_avrbb(:,:) = vor_avrb(:,:)
diff --git a/src/OCE/lib_fortran_generic.h90 b/src/OCE/lib_fortran_generic.h90
index 5201ad6d0..d09c2cd02 100644
--- a/src/OCE/lib_fortran_generic.h90
+++ b/src/OCE/lib_fortran_generic.h90
@@ -51,8 +51,8 @@
          iis = Nis0   ;   iie = Nie0
          ijs = Njs0   ;   ije = Nje0
       ELSE
-         iis = 1   ;   iie = jpi
-         ijs = 1   ;   ije = jpj
+         iis = 1   ;   iie = ipi
+         ijs = 1   ;   ije = ipj
       ENDIF
       !
       ctmp = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
-- 
GitLab