From 64c002477b70685072ce1e6a631d8052a72e324c Mon Sep 17 00:00:00 2001
From: clem <clement.rousset@locean.ipsl.fr>
Date: Tue, 6 Sep 2022 22:07:23 +0200
Subject: [PATCH] debug diagnostics and some ice options. And adapt
 field_def_oce.xml to the inner domains when necessary

---
 cfgs/SHARED/field_def_nemo-oce.xml | 329 ++++++++++++++---------------
 src/ICE/icedyn_rdgrft.F90          |  71 ++++---
 src/ICE/icedyn_rhg_eap.F90         | 262 ++++++++++++-----------
 src/OCE/DIA/diaar5.F90             | 179 ++++++++--------
 src/OCE/DIA/diacfl.F90             |  18 +-
 src/OCE/DIA/diahsb.F90             |  10 +-
 src/OCE/DIA/diaptr.F90             |  10 +-
 src/OCE/DIA/diawri.F90             | 277 ++++++++++++------------
 src/OCE/DYN/sshwzv.F90             |   3 +-
 src/OCE/LDF/ldftra.F90             | 143 +++++++------
 src/OCE/SBC/sbccpl.F90             |   2 +-
 src/OCE/TRA/traadv_cen.F90         |   2 +-
 src/OCE/TRA/traadv_cen_lf.F90      |   2 +-
 src/OCE/TRA/traadv_fct.F90         |  12 +-
 src/OCE/TRA/traadv_mus.F90         |   2 +-
 src/OCE/TRA/traadv_qck.F90         |   2 +-
 src/OCE/TRA/traadv_qck_lf.F90      |   2 +-
 src/OCE/TRA/traadv_ubs.F90         |   2 +-
 src/OCE/TRA/traadv_ubs_lf.F90      |   2 +-
 src/OCE/TRA/tradmp.F90             |   4 +-
 src/OCE/TRA/traldf_iso.F90         |   2 +-
 src/OCE/TRA/traldf_lap_blp.F90     |   2 +-
 src/OCE/TRA/traldf_triad.F90       |   2 +-
 src/OCE/ZDF/zdfphy.F90             |   2 +-
 24 files changed, 677 insertions(+), 665 deletions(-)

diff --git a/cfgs/SHARED/field_def_nemo-oce.xml b/cfgs/SHARED/field_def_nemo-oce.xml
index f1c958b87..bc6b20728 100644
--- a/cfgs/SHARED/field_def_nemo-oce.xml
+++ b/cfgs/SHARED/field_def_nemo-oce.xml
@@ -109,62 +109,60 @@ that are available in the tidal-forcing implementation (see
   <!-- T grid -->
 
   <field_group id="grid_T" grid_ref="grid_T_2D" >
-    <field id="e3t"          long_name="T-cell thickness"                    standard_name="cell_thickness"        unit="m"   grid_ref="grid_T_3D" />
-    <field id="e3ts"         long_name="T-cell thickness"   field_ref="e3t"  standard_name="cell_thickness"        unit="m"   grid_ref="grid_T_SFC"/>
-    <field id="e3t_0"        long_name="Initial T-cell thickness"            standard_name="ref_cell_thickness"    unit="m"   grid_ref="grid_T_3D" />
-    <field id="e3tb"         long_name="bottom T-cell thickness"             standard_name="bottom_cell_thickness" unit="m"   grid_ref="grid_T_2D"/>
-    <field id="e3t_300"      field_ref="e3t"                grid_ref="grid_T_zoom_300"       detect_missing_value="true" />
-    <field id="e3t_vsum300"  field_ref="e3t_300"            grid_ref="grid_T_vsum"   detect_missing_value="true" />
-    <field id="masscello"    long_name="Sea Water Mass per unit area"   standard_name="sea_water_mass_per_unit_area"   unit="kg/m2"   grid_ref="grid_T_3D"/>
-    <field id="volcello"     long_name="Ocean Volume"                   standard_name="ocean_volume"   unit="m3"       grid_ref="grid_T_3D"/>
-    <field id="toce"         long_name="temperature"                         standard_name="sea_water_potential_temperature"   unit="degC"     grid_ref="grid_T_3D"/>
-    <field id="toce_e3t"     long_name="temperature (thickness weighted)"                                                      unit="degC"     grid_ref="grid_T_3D" > toce * e3t </field >
-    <field id="soce"         long_name="salinity"                            standard_name="sea_water_practical_salinity"      unit="1e-3"     grid_ref="grid_T_3D"/>
-    <field id="soce_e3t"     long_name="salinity    (thickness weighted)"                                                      unit="1e-3"     grid_ref="grid_T_3D" > soce * e3t </field >
-
-    <field id="toce_e3t_300"      field_ref="toce_e3t"          unit="degree_C"     grid_ref="grid_T_zoom_300"      detect_missing_value="true" />
-    <field id="toce_e3t_vsum300"  field_ref="toce_e3t_300"      unit="degress_C*m"  grid_ref="grid_T_vsum"  detect_missing_value="true" />
-    <field id="toce_vmean300"     field_ref="toce_e3t_vsum300"  unit="degree_C"     grid_ref="grid_T_vsum"  detect_missing_value="true" > toce_e3t_vsum300/e3t_vsum300 </field>
+    <field id="e3t"          long_name="T-cell thickness"                    standard_name="cell_thickness"        unit="m"   grid_ref="grid_T_3D"       />
+    <field id="e3ts"         long_name="T-cell thickness"   field_ref="e3t"  standard_name="cell_thickness"        unit="m"   grid_ref="grid_T_SFC"      />
+    <field id="e3t_0"        long_name="Initial T-cell thickness"            standard_name="ref_cell_thickness"    unit="m"   grid_ref="grid_T_3D"       />
+    <field id="e3tb"         long_name="bottom T-cell thickness"             standard_name="bottom_cell_thickness" unit="m"   grid_ref="grid_T_2D_inner" />
+    <field id="e3t_300"      field_ref="e3t"                                                                                  grid_ref="grid_T_zoom_300" detect_missing_value="true" />
+    <field id="e3t_vsum300"  field_ref="e3t_300"                                                                              grid_ref="grid_T_vsum"     detect_missing_value="true" />
+
+    <field id="masscello"    long_name="Sea Water Mass per unit area"        standard_name="sea_water_mass_per_unit_area"    unit="kg/m2" grid_ref="grid_T_3D_inner"/>
+    <field id="volcello"     long_name="Ocean Volume"                        standard_name="ocean_volume"                    unit="m3"    grid_ref="grid_T_3D_inner"/>
+    <field id="toce"         long_name="temperature"                         standard_name="sea_water_potential_temperature" unit="degC"  grid_ref="grid_T_3D"/>
+    <field id="toce_e3t"     long_name="temperature (thickness weighted)"                                                    unit="degC"  grid_ref="grid_T_3D" > toce * e3t </field >
+    <field id="soce"         long_name="salinity"                            standard_name="sea_water_practical_salinity"    unit="1e-3"  grid_ref="grid_T_3D"/>
+    <field id="soce_e3t"     long_name="salinity    (thickness weighted)"                                                    unit="1e-3"  grid_ref="grid_T_3D" > soce * e3t </field >
+
+    <field id="toce_e3t_300"      field_ref="toce_e3t"          unit="degree_C"     grid_ref="grid_T_zoom_300" detect_missing_value="true" />
+    <field id="toce_e3t_vsum300"  field_ref="toce_e3t_300"      unit="degress_C*m"  grid_ref="grid_T_vsum"     detect_missing_value="true" />
+    <field id="toce_vmean300"     field_ref="toce_e3t_vsum300"  unit="degree_C"     grid_ref="grid_T_vsum"     detect_missing_value="true"  > toce_e3t_vsum300/e3t_vsum300 </field>
 
     <!-- AGRIF sponge -->
     <field id="agrif_spt"         long_name=" AGRIF t-sponge coefficient"   unit=" " />
 
     <!-- additions to diawri.F90 -->
-    <field id="sssgrad"     long_name="module of surface salinity gradient"              unit="1e-3/m"   grid_ref="grid_T_2D_inner"/>
-    <field id="sssgrad2"    long_name="square of module of surface salinity gradient"    unit="1e-6/m2"  grid_ref="grid_T_2D_inner"/>
-    <field id="ke"          long_name="kinetic energy"          standard_name="specific_kinetic_energy_of_sea_water"   unit="m2/s2"  grid_ref="grid_T_3D" />
-    <field id="ke_int"      long_name="vertical integration of kinetic energy"   unit="m3/s2"  grid_ref="grid_T_2D_inner" />
+    <field id="sssgrad"     long_name="module of surface salinity gradient"                                    unit="1e-3/m"   grid_ref="grid_T_2D_inner" />
+    <field id="sssgrad2"    long_name="square of module of surface salinity gradient"                          unit="1e-6/m2"  grid_ref="grid_T_2D_inner" />
+    <field id="ke"          long_name="kinetic energy"  standard_name="specific_kinetic_energy_of_sea_water"   unit="m2/s2"    grid_ref="grid_T_3D_inner" />
+    <field id="ke_int"      long_name="vertical integration of kinetic energy"                                 unit="m3/s2"    grid_ref="grid_T_2D_inner" />
+    <field id="taubot"      long_name="bottom stress module"                                                   unit="N/m2"     grid_ref="grid_T_2D_inner" />
 
     <!-- t-eddy viscosity coefficients (ldfdyn) -->
     <field id="ahmt_2d"      long_name=" surface t-eddy viscosity coefficient"   unit="m2/s or m4/s"                      />
     <field id="ahmt_3d"      long_name=" 3D      t-eddy viscosity coefficient"   unit="m2/s or m4/s"  grid_ref="grid_T_3D"/>
 
     <field id="sst"          long_name="Bulk sea surface temperature"                       standard_name="bulk_sea_surface_temperature"        unit="degC"     />
-    <field id="t_skin"       long_name="Skin temperature aka SSST"                          standard_name="skin_temperature"                    unit="degC"     />
+    <field id="sss"          long_name="sea surface salinity"                               standard_name="sea_surface_salinity"                unit="1e-3"     />
     <field id="sst2"         long_name="square of sea surface temperature"                  standard_name="square_of_sea_surface_temperature"   unit="degC2"     > sst * sst </field >
+    <field id="sss2"         long_name="square of sea surface salinity"                                                                         unit="1e-6"      > sss * sss </field >
     <field id="sstmax"       long_name="max of sea surface temperature"   field_ref="sst"   operation="maximum"                                                 />
+    <field id="sssmax"       long_name="max of sea surface salinity"      field_ref="sss"   operation="maximum"                                                 />
     <field id="sstmin"       long_name="min of sea surface temperature"   field_ref="sst"   operation="minimum"                                                 />
+    <field id="sssmin"       long_name="min of sea surface salinity"      field_ref="sss"   operation="minimum"                                                 />
     <field id="sstgrad"      long_name="module of sst gradient"                                                                                 unit="degC/m"   grid_ref="grid_T_2D_inner" />
     <field id="sstgrad2"     long_name="square of module of sst gradient"                                                                       unit="degC2/m2" grid_ref="grid_T_2D_inner" />
     <field id="sbt"          long_name="sea bottom temperature"                                                                                 unit="degC"     grid_ref="grid_T_2D_inner" />
-    <field id="tosmint"      long_name="vertical integral of temperature times density"     standard_name="integral_wrt_depth_of_product_of_density_and_potential_temperature"  unit="(kg m2) degree_C" grid_ref="grid_T_2D_inner" />
-    <field id="sst_wl"       long_name="Delta SST of warm layer"                                                                                unit="degC"     />
-    <field id="sst_cs"       long_name="Delta SST of cool skin"                                                                                 unit="degC"     />
-    <field id="temp_3m"      long_name="temperature at 3m"                                                                                      unit="degC"     />
-
-    <field id="sss"          long_name="sea surface salinity"                               standard_name="sea_surface_salinity"                unit="1e-3"     />
-    <field id="sss2"         long_name="square of sea surface salinity"                                                                         unit="1e-6"      > sss * sss </field >
-    <field id="sssmax"       long_name="max of sea surface salinity"      field_ref="sss"   operation="maximum"                                                 />
-    <field id="sssmin"       long_name="min of sea surface salinity"      field_ref="sss"   operation="minimum"                                                 />
-    <field id="sbs"          long_name="sea bottom salinity"                                                                                    unit="0.001"   grid_ref="grid_T_2D_inner" />
-    <field id="somint"       long_name="vertical integral of salinity times density"        standard_name="integral_wrt_depth_of_product_of_density_and_salinity"  unit="(kg m2) x (1e-3)" grid_ref="grid_T_2D_inner" />
+    <field id="sbs"          long_name="sea bottom salinity"                                                                                    unit="0.001"    grid_ref="grid_T_2D_inner" />
+    <field id="sst_wl"       long_name="Delta SST of warm layer"                                                                                unit="degC"     grid_ref="grid_T_2D_inner" />
+    <field id="sst_cs"       long_name="Delta SST of cool skin"                                                                                 unit="degC"     grid_ref="grid_T_2D_inner" />
 
-    <field id="taubot"       long_name="bottom stress module"   unit="N/m2"   grid_ref="grid_T_2D_inner" />
+    <field id="tosmint"      long_name="vertical integral of temperature times density" standard_name="integral_wrt_depth_of_product_of_density_and_potential_temperature" unit="(kg m2) degree_C" grid_ref="grid_T_2D_inner" />
+    <field id="somint"       long_name="vertical integral of salinity times density"    standard_name="integral_wrt_depth_of_product_of_density_and_salinity"              unit="(kg m2) x (1e-3)" grid_ref="grid_T_2D_inner" />
 
     <!-- Case EOS = TEOS-10 : output potential temperature -->
-    <field id="toce_pot"     long_name="Sea Water Potential Temperature"              standard_name="sea_water_potential_temperature"   unit="degC"     grid_ref="grid_T_3D"/>
-    <field id="sst_pot"      long_name="potential sea surface temperature"             standard_name="sea_surface_temperature"             unit="degC"     />
-    <field id="tosmint_pot"  long_name="vertical integral of potential temperature times density"   standard_name="integral_wrt_depth_of_product_of_density_and_potential_temperature"  unit="(kg m2) degree_C" />
+    <field id="toce_pot"     long_name="Sea Water Potential Temperature"                          standard_name="sea_water_potential_temperature"                                    unit="degC"             grid_ref="grid_T_3D_inner"/>
+    <field id="sst_pot"      long_name="potential sea surface temperature"                        standard_name="sea_surface_temperature"                                            unit="degC"             grid_ref="grid_T_2D_inner"/>
+    <field id="tosmint_pot"  long_name="vertical integral of potential temperature times density" standard_name="integral_wrt_depth_of_product_of_density_and_potential_temperature" unit="(kg m2) degree_C" grid_ref="grid_T_2D_inner"/>
 
     <field id="ht"           long_name="water column height at T point"                     standard_name="water_column_height_T"                      unit="m" />
     <field id="ssh"          long_name="sea surface height"                                 standard_name="sea_surface_height_above_geoid"             unit="m" />
@@ -172,13 +170,13 @@ that are available in the tidal-forcing implementation (see
     <field id="wetdep"       long_name="wet depth"                                          standard_name="wet_depth"                                  unit="m" />
     <field id="sshmax"       long_name="max of sea surface height"        field_ref="ssh"   operation="maximum"                                                 />
 
-    <field id="mldkz5"       long_name="Turbocline depth (Kz = 5e-4)"                       standard_name="ocean_mixed_layer_thickness_defined_by_vertical_tracer_diffusivity"                unit="m"  grid_ref="grid_T_2D_inner" />
-    <field id="mldr10_1"     long_name="Mixed Layer Depth (dsigma = 0.01 wrt 10m)"          standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                                unit="m"  grid_ref="grid_T_2D_inner" />
-    <field id="mldr10_1max"  long_name="Max of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="maximum"                                                                  grid_ref="grid_T_2D_inner" />
-    <field id="mldr10_1min"  long_name="Min of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="minimum"                                                                  grid_ref="grid_T_2D_inner" />
-    <field id="heatc"        long_name="Heat content vertically integrated"                 standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"      grid_ref="grid_T_2D_inner" />
-    <field id="saltc"        long_name="Salt content vertically integrated"                                                                                                                   unit="PSU*kg/m2"    grid_ref="grid_T_2D_inner" />
-    <field id="salt2c"       long_name="square of Salt content vertically integrated"                                                                                                         unit="PSU2*kg/m2"   grid_ref="grid_T_2D_inner" />
+    <field id="mldkz5"       long_name="Turbocline depth (Kz = 5e-4)"                     standard_name="ocean_mixed_layer_thickness_defined_by_vertical_tracer_diffusivity"              unit="m"          grid_ref="grid_T_2D_inner" />
+    <field id="mldr10_1"     long_name="Mixed Layer Depth (dsigma = 0.01 wrt 10m)"        standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta"                              unit="m"          grid_ref="grid_T_2D_inner" />
+    <field id="mldr10_1max"  long_name="Max of Mixed Layer Depth (dsigma = 0.01 wrt 10m)" field_ref="mldr10_1"   operation="maximum"                                                                        grid_ref="grid_T_2D_inner" />
+    <field id="mldr10_1min"  long_name="Min of Mixed Layer Depth (dsigma = 0.01 wrt 10m)" field_ref="mldr10_1"   operation="minimum"                                                                        grid_ref="grid_T_2D_inner" />
+    <field id="heatc"        long_name="Heat content vertically integrated"               standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content" unit="J/m2"       grid_ref="grid_T_2D_inner" />
+    <field id="saltc"        long_name="Salt content vertically integrated"                                                                                                               unit="PSU*kg/m2"  grid_ref="grid_T_2D_inner" />
+    <field id="salt2c"       long_name="square of Salt content vertically integrated"                                                                                                     unit="PSU2*kg/m2" grid_ref="grid_T_2D_inner" />
 
     <!-- EOS -->
     <field id="alpha"        long_name="thermal expansion"                                                         unit="degC-1" grid_ref="grid_T_3D" />
@@ -217,28 +215,28 @@ that are available in the tidal-forcing implementation (see
     <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" />
+    <field id="botpres"      long_name="Sea Water Pressure at Sea Floor"          standard_name="sea_water_pressure_at_sea_floor"                    unit="dbar" grid_ref="grid_T_2D_inner" />
     <field id="sshdyn"       long_name="dynamic sea surface height"               standard_name="dynamic_sea_surface_height_above_geoid"             unit="m"    />
     <field id="sshdyn2"      long_name="square of dynamic sea surface height"     standard_name="dynamic_sea_surface_height_above_geoid_squared"     unit="m2"    > sshdyn * sshdyn </field>
-    <field id="tnpeo"      long_name="Tendency of ocean potential energy content"                                                                    unit="W/m2" />
+    <field id="tnpeo"      long_name="Tendency of ocean potential energy content"                                                                    unit="W/m2" grid_ref="grid_T_2D_inner" />
 
     <!-- variables available ln_linssh=.FALSE. -->
-    <field id="tpt_dep"      long_name="T-point depth"                  standard_name="depth_below_geoid"   unit="m"   grid_ref="grid_T_3D" />
+    <field id="tpt_dep"      long_name="T-point depth"                  standard_name="depth_below_geoid"   unit="m"   grid_ref="grid_T_3D_inner" />
     <field id="e3tdef"       long_name="T-cell thickness deformation"                                       unit="%"   grid_ref="grid_T_3D" />
 
     <!-- variables available with ln_diacfl=.true. -->
-    <field id="cfl_cu"       long_name="u-courant number"   unit="#" />
-    <field id="cfl_cv"       long_name="v-courant number"   unit="#" />
-    <field id="cfl_cw"       long_name="w-courant number"   unit="#" />
+    <field id="cfl_cu"       long_name="u-courant number"   unit="#" grid_ref="grid_T_2D_inner" />
+    <field id="cfl_cv"       long_name="v-courant number"   unit="#" grid_ref="grid_T_2D_inner" />
+    <field id="cfl_cw"       long_name="w-courant number"   unit="#" grid_ref="grid_T_2D_inner" />
 
     <!-- variables available with ln_zdfmfc=.true. -->
     <field id="mf_Tp"       long_name="plume_temperature"      standard_name="plume_temperature"     unit="degC"   grid_ref="grid_T_3D" />
     <field id="mf_Sp"       long_name="plume_salinity"         standard_name="plume_salinity"        unit="1e-3"   grid_ref="grid_T_3D" />
-    <field id="mf_mf"       long_name="mass flux"              standard_name="mf_mass_flux"          unit="m"      grid_ref="grid_T_3D" />
+    <field id="mf_mf"       long_name="mass flux"              standard_name="mf_mass_flux"          unit="m"      grid_ref="grid_T_3D_inner" />
 
     <!-- fluxes from damping -->
-    <field id="sflx_dmp_cea"  long_name="salt flux due to damping"  standard_name="salt_flux_due_to_damping"          unit="g/m2/s"   />
-    <field id="hflx_dmp_cea"  long_name="heat flux due to damping"  standard_name="heat_flux_due_to_damping"          unit="W/m2"     />
+    <field id="sflx_dmp_cea"  long_name="salt flux due to damping"  standard_name="salt_flux_due_to_damping"          unit="g/m2/s"  grid_ref="grid_T_2D_inner" />
+    <field id="hflx_dmp_cea"  long_name="heat flux due to damping"  standard_name="heat_flux_due_to_damping"          unit="W/m2"    grid_ref="grid_T_2D_inner" />
 
     <!-- * variable related to ice shelf forcing * -->
 
@@ -439,43 +437,44 @@ that are available in the tidal-forcing implementation (see
       <field id="ssh_ib"       long_name="Inverse barometer sea surface height"  standard_name="sea_surface_height_correction_due_to_air_pressure_at_low_frequency"   unit="m" grid_ref="grid_T_2D" />
 
       <!-- *_oce variables available with ln_blk_clio or ln_blk_core -->
-      <field id="rho_air"      long_name="Air density at 10m above sea surface"         standard_name="rho_air_10m"                                        unit="kg/m3" />
-      <field id="dt_skin"      long_name="SSST-SST temperature difference"              standard_name="SSST-SST"                                             unit="K"   />
-      <field id="qlw_oce"      long_name="Longwave Downward Heat Flux over open ocean"  standard_name="surface_net_downward_longwave_flux"                 unit="W/m2"  />
-      <field id="qsb_oce"      long_name="Sensible Downward Heat Flux over open ocean"  standard_name="surface_downward_sensible_heat_flux"                unit="W/m2"  />
-      <field id="qla_oce"      long_name="Latent Downward Heat Flux over open ocean"    standard_name="surface_downward_latent_heat_flux"                  unit="W/m2"  />
-      <field id="evap_oce"     long_name="Evaporation over open ocean"                  standard_name="evaporation"                                        unit="kg/m2/s" />
-      <field id="qt_oce"       long_name="total flux at ocean surface"                  standard_name="surface_downward_heat_flux_in_sea_water"            unit="W/m2"  />
-      <field id="qsr_oce"      long_name="solar heat flux at ocean surface"             standard_name="net_downward_shortwave_flux_at_sea_water_surface"   unit="W/m2"  />
-      <field id="qns_oce"      long_name="non-solar heat flux at ocean surface (including E-P)"                                                            unit="W/m2"  />
-      <field id="qemp_oce"     long_name="Downward Heat Flux from E-P over open ocean"                                                                     unit="W/m2"  />
-      <field id="taum_oce"     long_name="wind stress module over open ocean"           standard_name="magnitude_of_surface_downward_stress"               unit="N/m2"  />
+      <field id="rho_air"      long_name="Air density at 10m above sea surface"         standard_name="rho_air_10m"                                        unit="kg/m3"  />
+      <field id="t_skin"       long_name="Skin temperature aka SSST"                    standard_name="skin_temperature"                                   unit="degC"   />
+      <field id="dt_skin"      long_name="SSST-SST temperature difference"              standard_name="SSST-SST"                                           unit="K"      />
+      <field id="qlw_oce"      long_name="Longwave Downward Heat Flux over open ocean"  standard_name="surface_net_downward_longwave_flux"                 unit="W/m2"   />
+      <field id="qsb_oce"      long_name="Sensible Downward Heat Flux over open ocean"  standard_name="surface_downward_sensible_heat_flux"                unit="W/m2"   />
+      <field id="qla_oce"      long_name="Latent Downward Heat Flux over open ocean"    standard_name="surface_downward_latent_heat_flux"                  unit="W/m2"   />
+      <field id="evap_oce"     long_name="Evaporation over open ocean"                  standard_name="evaporation"                                        unit="kg/m2/s"/>
+      <field id="qt_oce"       long_name="total flux at ocean surface"                  standard_name="surface_downward_heat_flux_in_sea_water"            unit="W/m2"   />
+      <field id="qsr_oce"      long_name="solar heat flux at ocean surface"             standard_name="net_downward_shortwave_flux_at_sea_water_surface"   unit="W/m2"   />
+      <field id="qns_oce"      long_name="non-solar heat flux at ocean surface (including E-P)"                                                            unit="W/m2"   />
+      <field id="qemp_oce"     long_name="Downward Heat Flux from E-P over open ocean"                                                                     unit="W/m2"   />
+      <field id="taum_oce"     long_name="wind stress module over open ocean"           standard_name="magnitude_of_surface_downward_stress"               unit="N/m2"   />
       <field id="utau_oce"     long_name="Wind Stress along i-axis over open ocean (T-points)"  standard_name="surf_down_x_stress_open_oce_Tpoints"        unit="N/m2"   />
-      <field id="vtau_oce"     long_name="Wind Stress along j-axis over open ocean (T-points)"  standard_name="surf_down_y_stress_open_oce_Tpoints"        unit="N/m2"    />
+      <field id="vtau_oce"     long_name="Wind Stress along j-axis over open ocean (T-points)"  standard_name="surf_down_y_stress_open_oce_Tpoints"        unit="N/m2"   />
 
       <!-- variables computed by the bulk parameterization algorithms (ln_blk) -->
-      <field id="Cd_oce"      long_name="Drag coefficient over open ocean"              standard_name="drag_coefficient_water"                unit=""  />
-      <field id="Ce_oce"      long_name="Evaporaion coefficient over open ocean"        standard_name="evap_coefficient_water"                unit=""  />
-      <field id="Ch_oce"      long_name="Sensible heat coefficient over open ocean"     standard_name="sensible_heat_coefficient_water"       unit=""  />
+      <field id="Cd_oce"      long_name="Drag coefficient over open ocean"              standard_name="drag_coefficient_water"                unit=""     />
+      <field id="Ce_oce"      long_name="Evaporaion coefficient over open ocean"        standard_name="evap_coefficient_water"                unit=""     />
+      <field id="Ch_oce"      long_name="Sensible heat coefficient over open ocean"     standard_name="sensible_heat_coefficient_water"       unit=""     />
       <field id="theta_zt"    long_name="Potential air temperature at z=zt"             standard_name="potential_air_temperature_at_zt"       unit="degC" />
-      <field id="q_zt"        long_name="Specific air humidity at z=zt"                 standard_name="specific_air_humidity_at_zt"           unit="kg/kg" />
+      <field id="q_zt"        long_name="Specific air humidity at z=zt"                 standard_name="specific_air_humidity_at_zt"           unit="kg/kg"/>
       <field id="theta_zu"    long_name="Potential air temperature at z=zu"             standard_name="potential_air_temperature_at_zu"       unit="degC" />
-      <field id="q_zu"        long_name="Specific air humidity at z=zu"                 standard_name="specific_air_humidity_at_zu"           unit="kg/kg" />
-      <field id="ssq"         long_name="Saturation specific humidity of air at z=0"    standard_name="surface_air_saturation_spec_humidity"  unit="kg/kg" />
-      <field id="wspd_blk"    long_name="Bulk wind speed at z=zu"                       standard_name="bulk_wind_speed_at_zu"                 unit="m/s"   />
+      <field id="q_zu"        long_name="Specific air humidity at z=zu"                 standard_name="specific_air_humidity_at_zu"           unit="kg/kg"/>
+      <field id="ssq"         long_name="Saturation specific humidity of air at z=0"    standard_name="surface_air_saturation_spec_humidity"  unit="kg/kg"/>
+      <field id="wspd_blk"    long_name="Bulk wind speed at z=zu"                       standard_name="bulk_wind_speed_at_zu"                 unit="m/s"  />
       <!-- ln_blk + key_si3 -->
-      <field id="Cd_ice"      long_name="Drag coefficient over ice"                     standard_name="drag_coefficient_ice"                 unit=""  />
-      <field id="Ce_ice"      long_name="Evaporaion coefficient over ice"               standard_name="evap_coefficient_ice"                 unit=""  />
-      <field id="Ch_ice"      long_name="Sensible heat coefficient over ice"            standard_name="sensible_heat_coefficient_ice"        unit=""  />
+      <field id="Cd_ice"      long_name="Drag coefficient over ice"                     standard_name="drag_coefficient_ice"                  unit=""  />
+      <field id="Ce_ice"      long_name="Evaporaion coefficient over ice"               standard_name="evap_coefficient_ice"                  unit=""  />
+      <field id="Ch_ice"      long_name="Sensible heat coefficient over ice"            standard_name="sensible_heat_coefficient_ice"         unit=""  />
 
       <!-- available key_oasis3 -->
       <field id="snow_ao_cea"  long_name="Snow over ice-free ocean (cell average)"   standard_name="snowfall_flux"                             unit="kg/m2/s"  />
       <field id="snow_ai_cea"  long_name="Snow over sea-ice (cell average)"          standard_name="snowfall_flux"                             unit="kg/m2/s"  />
       <field id="subl_ai_cea"  long_name="Sublimation over sea-ice (cell average)"   standard_name="surface_snow_and_ice_sublimation_flux"     unit="kg/m2/s"  />
       <field id="icealb_cea"   long_name="Ice albedo (cell average)"                 standard_name="sea_ice_albedo"                            unit="1"        />
-      <field id="calving_cea"  long_name="Calving"                                   standard_name="water_flux_into_sea_water_from_icebergs"   unit="kg/m2/s" grid_ref="grid_T_2D" />
-      <field id="iceberg_cea"  long_name="Iceberg"                                   standard_name="water_flux_into_sea_water_from_icebergs"   unit="kg/m2/s" grid_ref="grid_T_2D" />
-      <field id="iceshelf_cea" long_name="Iceshelf"                                  standard_name="water_flux_into_sea_water_from_iceshelf"   unit="kg/m2/s" grid_ref="grid_T_2D" />
+      <field id="calving_cea"  long_name="Calving"                                   standard_name="water_flux_into_sea_water_from_icebergs"   unit="kg/m2/s"  />
+      <field id="iceberg_cea"  long_name="Iceberg"                                   standard_name="water_flux_into_sea_water_from_icebergs"   unit="kg/m2/s"  />
+      <field id="iceshelf_cea" long_name="Iceshelf"                                  standard_name="water_flux_into_sea_water_from_iceshelf"   unit="kg/m2/s"  />
 
 
       <!-- available if key_oasis3 + conservative method -->
@@ -515,10 +514,10 @@ that are available in the tidal-forcing implementation (see
       <field id="ice_cover"    long_name="Ice fraction"  standard_name="sea_ice_area_fraction"  unit="1"  grid_ref="grid_T_2D" />
 
       <!-- dilution -->
-      <field id="emp_x_sst"    long_name="Concentration/Dilution term on SST"                                                                                              unit="kg*degC/m2/s" />
-      <field id="emp_x_sss"    long_name="Concentration/Dilution term on SSS"                                                                                              unit="kg*1e-3/m2/s" />
-      <field id="rnf_x_sst"    long_name="Runoff term on SST"                                                                                                              unit="kg*degC/m2/s" />
-      <field id="rnf_x_sss"    long_name="Runoff term on SSS"                                                                                                              unit="kg*1e-3/m2/s" />
+      <field id="emp_x_sst"    long_name="Concentration/Dilution term on SST"   unit="kg*degC/m2/s" grid_ref="grid_T_2D" />
+      <field id="emp_x_sss"    long_name="Concentration/Dilution term on SSS"   unit="kg*1e-3/m2/s" grid_ref="grid_T_2D" />
+      <field id="rnf_x_sst"    long_name="Runoff term on SST"                   unit="kg*degC/m2/s" grid_ref="grid_T_2D" />
+      <field id="rnf_x_sss"    long_name="Runoff term on SSS"                   unit="kg*1e-3/m2/s" grid_ref="grid_T_2D" />
 
       <!-- sbcssm variables -->
       <field id="sst_m"    unit="degC" grid_ref="grid_T_2D" />
@@ -582,17 +581,17 @@ that are available in the tidal-forcing implementation (see
 
   <field_group id="grid_U"   grid_ref="grid_U_2D">
     <field id="hu"            long_name="water column height at U point"                         standard_name="water_column_height_U"       unit="m" />
-    <field id="e2u"           long_name="U-cell width in meridional direction"                   standard_name="cell_width"                  unit="m"                               />
-    <field id="e3u"           long_name="U-cell thickness"                                       standard_name="cell_thickness"              unit="m"          grid_ref="grid_U_3D" />
-    <field id="e3u_0"         long_name="Initial U-cell thickness"                               standard_name="ref_cell_thickness"          unit="m"          grid_ref="grid_U_3D"/>
-    <field id="uoce"          long_name="ocean current along i-axis"                             standard_name="sea_water_x_velocity"        unit="m/s"        grid_ref="grid_U_3D" />
-    <field id="uoce_e3u"      long_name="ocean current along i-axis  (thickness weighted)"                                                   unit="m/s"        grid_ref="grid_U_3D"  > uoce * e3u </field>
-    <field id="uoce_e3u_vsum" long_name="ocean current along i-axis * e3u summed on the vertical"  field_ref="uoce_e3u"    unit="m3/s"       grid_ref="grid_U_vsum"/>
-    <field id="uocetr_vsum"   long_name="ocean transport along i-axis  summed on the vertical"         field_ref="e2u"       unit="m3/s"> this * uoce_e3u_vsum  </field>
-
-    <field id="uocetr_vsum_op"    long_name="ocean current along i-axis * e3u * e2u summed on the vertical"  read_access="true"  freq_op="1mo"    field_ref="e2u"       unit="m3/s"> @uocetr_vsum </field>
+    <field id="e2u"           long_name="U-cell width in meridional direction"                   standard_name="cell_width"                  unit="m"                                 />
+    <field id="e3u"           long_name="U-cell thickness"                                       standard_name="cell_thickness"              unit="m"          grid_ref="grid_U_3D"   />
+    <field id="e3u_0"         long_name="Initial U-cell thickness"                               standard_name="ref_cell_thickness"          unit="m"          grid_ref="grid_U_3D"   />
+    <field id="uoce"          long_name="ocean current along i-axis"                             standard_name="sea_water_x_velocity"        unit="m/s"        grid_ref="grid_U_3D"   />
+    <field id="uoce_e3u"      long_name="ocean current along i-axis  (thickness weighted)"                                                   unit="m/s"        grid_ref="grid_U_3D"    > uoce * e3u </field>
+    <field id="uoce_e3u_vsum" long_name="ocean current along i-axis * e3u summed on the vertical"  field_ref="uoce_e3u"                      unit="m3/s"       grid_ref="grid_U_vsum" />
+    <field id="uocetr_vsum"   long_name="ocean transport along i-axis  summed on the vertical"     field_ref="e2u"                           unit="m3/s"                               > this * uoce_e3u_vsum  </field>
+
+    <field id="uocetr_vsum_op"    long_name="ocean current along i-axis * e3u * e2u summed on the vertical"  read_access="true"  freq_op="1mo" field_ref="e2u" unit="m3/s"             > @uocetr_vsum </field>
     <field id="uocetr_vsum_cumul" long_name="ocean current along i-axis * e3u * e2u cumulated from southwest point" freq_offset="_reset_" operation="instant" freq_op="1mo"  unit="m3/s" />
-    <field id="msftbarot"         long_name="ocean_barotropic_mass_streamfunction"   unit="kg s-1" > uocetr_vsum_cumul * $rho0 </field>
+    <field id="msftbarot"         long_name="ocean_barotropic_mass_streamfunction"                                                                             unit="kg s-1"           > uocetr_vsum_cumul * $rho0 </field>
 
 
     <field id="ssu"          long_name="ocean surface current along i-axis"                                                                 unit="m/s"                             />
@@ -606,21 +605,21 @@ that are available in the tidal-forcing implementation (see
     <field id="agrif_spu"    long_name=" AGRIF u-sponge coefficient"   unit=" " />
     <!-- u-eddy diffusivity coefficients (available if ln_traldf_OFF=F) -->
     <field id="ahtu_2d"      long_name=" surface u-eddy diffusivity coefficient"   unit="m2/s or m4/s" />
-    <field id="ahtu_3d"      long_name=" 3D u-EIV coefficient"                     unit="m2/s or m4/s"      grid_ref="grid_U_3D"/>
+    <field id="ahtu_3d"      long_name=" 3D u-EIV coefficient"                     unit="m2/s or m4/s" grid_ref="grid_U_3D" />
     <!-- u-eiv diffusivity coefficients (available if ln_ldfeiv=F) -->
-    <field id="aeiu_2d"      long_name=" surface u-EIV coefficient"                unit="m2/s" />
-    <field id="aeiu_3d"      long_name=" 3D u-EIV coefficient"                     unit="m2/s"              grid_ref="grid_U_3D"/>
+    <field id="aeiu_2d"      long_name=" surface u-EIV coefficient"                unit="m2/s"         />
+    <field id="aeiu_3d"      long_name=" 3D u-EIV coefficient"                     unit="m2/s"         grid_ref="grid_U_3D" />
 
     <!-- variables available with MLE (ln_mle=T) -->
-    <field id="psiu_mle"     long_name="MLE streamfunction along i-axis"   unit="m3/s"   grid_ref="grid_U_3D" />
+    <field id="psiu_mle"     long_name="MLE streamfunction along i-axis"           unit="m3/s"         grid_ref="grid_U_3D" />
 
     <!-- uoce_eiv: available EIV (ln_ldfeiv=T and ln_ldfeiv_dia=T) -->
-    <field id="uoce_eiv"      long_name="EIV ocean current along i-axis"                                  standard_name="bolus_sea_water_x_velocity"                     unit="m/s"   grid_ref="grid_U_3D" />
-    <field id="ueiv_masstr"   long_name="EIV Ocean Mass X Transport"                                      standard_name="bolus_ocean_mass_x_transport"                   unit="kg/s"  grid_ref="grid_U_3D" />
-    <field id="ueiv_heattr"   long_name="ocean bolus heat transport along i-axis"                         standard_name="ocean_heat_x_transport_due_to_bolus_advection"  unit="W"                         />
-    <field id="ueiv_salttr"   long_name="ocean bolus salt transport along i-axis"                         standard_name="ocean_salt_x_transport_due_to_bolus_advection"  unit="Kg"                        />
-    <field id="ueiv_heattr3d" long_name="ocean bolus heat transport along i-axis"                         standard_name="ocean_heat_x_transport_due_to_bolus_advection"  unit="W"     grid_ref="grid_U_3D" />
-    <field id="ueiv_salttr3d" long_name="ocean bolus salt transport along i-axis"                         standard_name="ocean_salt_x_transport_due_to_bolus_advection"  unit="kg"    grid_ref="grid_U_3D" />
+    <field id="uoce_eiv"      long_name="EIV ocean current along i-axis"            standard_name="bolus_sea_water_x_velocity"                     unit="m/s"   grid_ref="grid_U_3D_inner" />
+    <field id="ueiv_masstr"   long_name="EIV Ocean Mass X Transport"                standard_name="bolus_ocean_mass_x_transport"                   unit="kg/s"  grid_ref="grid_U_3D_inner" />
+    <field id="ueiv_heattr"   long_name="ocean bolus heat transport along i-axis"   standard_name="ocean_heat_x_transport_due_to_bolus_advection"  unit="W"     grid_ref="grid_U_2D_inner" />
+    <field id="ueiv_salttr"   long_name="ocean bolus salt transport along i-axis"   standard_name="ocean_salt_x_transport_due_to_bolus_advection"  unit="Kg"    grid_ref="grid_U_2D_inner" />
+    <field id="ueiv_heattr3d" long_name="ocean bolus heat transport along i-axis"   standard_name="ocean_heat_x_transport_due_to_bolus_advection"  unit="W"     grid_ref="grid_U_3D_inner" />
+    <field id="ueiv_salttr3d" long_name="ocean bolus salt transport along i-axis"   standard_name="ocean_salt_x_transport_due_to_bolus_advection"  unit="kg"    grid_ref="grid_U_3D_inner" />
 
     <!-- uoce_bbl: available with ln_trabbl=T and nn_bbl_adv=1 -->
     <field id="uoce_bbl"     long_name="BBL ocean current along i-axis"    unit="m/s"  />
@@ -635,20 +634,20 @@ that are available in the tidal-forcing implementation (see
     <field id="utbl"         long_name="zonal current in the Losh tbl"     unit="m/s" />
 
     <!-- variables available with diaar5 -->
-    <field id="u_masstr"      long_name="Ocean Mass X Transport"                                          standard_name="ocean_mass_x_transport"                         unit="kg/s"   grid_ref="grid_U_3D" />
-    <field id="u_masstr_vint" long_name="vertical integral of ocean eulerian mass transport along i-axis" standard_name="vertical_integral_of_ocean_mass_x_transport"    unit="kg/s"   grid_ref="grid_U_2D_inner" />
-    <field id="u_heattr"      long_name="ocean eulerian heat transport along i-axis"                      standard_name="ocean_heat_x_transport"                         unit="W"      grid_ref="grid_U_2D_inner" />
+    <field id="u_masstr"      long_name="Ocean Mass X Transport"                                          standard_name="ocean_mass_x_transport"                         unit="kg/s"      grid_ref="grid_U_3D_inner" />
+    <field id="u_masstr_vint" long_name="vertical integral of ocean eulerian mass transport along i-axis" standard_name="vertical_integral_of_ocean_mass_x_transport"    unit="kg/s"      grid_ref="grid_U_2D_inner" />
+    <field id="u_heattr"      long_name="ocean eulerian heat transport along i-axis"                      standard_name="ocean_heat_x_transport"                         unit="W"         grid_ref="grid_U_2D_inner" />
     <field id="u_salttr"      long_name="ocean eulerian salt transport along i-axis"                      standard_name="ocean_salt_x_transport"                         unit="1e-3*kg/s" grid_ref="grid_U_2D_inner" />
-    <field id="uadv_heattr"   long_name="ocean advective heat transport along i-axis"                     standard_name="advectice_ocean_heat_x_transport"               unit="W"                         />
-    <field id="uadv_salttr"   long_name="ocean advective salt transport along i-axis"                     standard_name="advectice_ocean_salt_x_transport"               unit="1e-3*kg/s"                 />
-    <field id="udiff_heattr"  long_name="ocean diffusion heat transport along i-axis"                     standard_name="ocean_heat_x_transport_due_to_diffusion"        unit="W"                         />
-    <field id="udiff_salttr"  long_name="ocean diffusion salt transport along i-axis"                     standard_name="ocean_salt_x_transport_due_to_diffusion"        unit="1e-3*kg/s"                 />
+    <field id="uadv_heattr"   long_name="ocean advective heat transport along i-axis"                     standard_name="advectice_ocean_heat_x_transport"               unit="W"         grid_ref="grid_U_2D_inner" />
+    <field id="uadv_salttr"   long_name="ocean advective salt transport along i-axis"                     standard_name="advectice_ocean_salt_x_transport"               unit="1e-3*kg/s" grid_ref="grid_U_2D_inner" />
+    <field id="udiff_heattr"  long_name="ocean diffusion heat transport along i-axis"                     standard_name="ocean_heat_x_transport_due_to_diffusion"        unit="W"         grid_ref="grid_U_2D_inner" />
+    <field id="udiff_salttr"  long_name="ocean diffusion salt transport along i-axis"                     standard_name="ocean_salt_x_transport_due_to_diffusion"        unit="1e-3*kg/s" grid_ref="grid_U_2D_inner" />
   </field_group>
 
   <!-- V grid -->
 
   <field_group id="grid_V"   grid_ref="grid_V_2D">
-    <field id="e1v"          long_name="V-cell width in longitudinal direction"                 standard_name="cell_width"                  unit="m"                              />
+    <field id="e1v"          long_name="V-cell width in longitudinal direction"                 standard_name="cell_width"                  unit="m"                               />
     <field id="e3v"          long_name="V-cell thickness"                                       standard_name="cell_thickness"              unit="m"          grid_ref="grid_V_3D" />
     <field id="e3v_0"        long_name="Initial V-cell thickness"                               standard_name="ref_cell_thickness"          unit="m"          grid_ref="grid_V_3D" />
     <field id="hv"            long_name="water column height at V point"                        standard_name="water_column_height_V"       unit="m" />
@@ -665,21 +664,21 @@ that are available in the tidal-forcing implementation (see
     <field id="agrif_spv"    long_name=" AGRIF v-sponge coefficient"   unit=" " />
     <!-- v-eddy diffusivity coefficients (available if ln_traldf_OFF=F) -->
     <field id="ahtv_2d"      long_name=" surface v-eddy diffusivity coefficient"     unit="m2/s or (m4/s)^1/2" />
-    <field id="ahtv_3d"      long_name=" 3D v-eddy diffusivity coefficient"          unit="m2/s or (m4/s)^1/2"           grid_ref="grid_V_3D"/>
+    <field id="ahtv_3d"      long_name=" 3D v-eddy diffusivity coefficient"          unit="m2/s or (m4/s)^1/2" grid_ref="grid_V_3D" />
     <!-- v-eiv diffusivity coefficients (available if ln_ldfeiv=F) -->
-    <field id="aeiv_2d"      long_name=" surface v-EIV coefficient"                  unit="m2/s" />
-    <field id="aeiv_3d"      long_name=" 3D v-EIV coefficient"                       unit="m2/s"                         grid_ref="grid_V_3D" />
+    <field id="aeiv_2d"      long_name=" surface v-EIV coefficient"                  unit="m2/s"               />
+    <field id="aeiv_3d"      long_name=" 3D v-EIV coefficient"                       unit="m2/s"               grid_ref="grid_V_3D" />
 
     <!-- variables available with MLE (ln_mle=T) -->
-    <field id="psiv_mle"     long_name="MLE streamfunction along j-axis"   unit="m3/s"   grid_ref="grid_V_3D" />
+    <field id="psiv_mle"     long_name="MLE streamfunction along j-axis"             unit="m3/s"               grid_ref="grid_V_3D" />
 
     <!-- voce_eiv: available EIV (ln_ldfeiv=T and ln_ldfeiv_dia=T)  -->
-    <field id="voce_eiv"     long_name="EIV ocean current along j-axis"  standard_name="bolus_sea_water_y_velocity"     unit="m/s"   grid_ref="grid_V_3D" />
-    <field id="veiv_masstr"  long_name="EIV Ocean Mass Y Transport"      standard_name="bolus_ocean_mass_y_transport"   unit="kg/s"  grid_ref="grid_V_3D" />
-    <field id="veiv_heattr"   long_name="ocean bolus heat transport along j-axis"       standard_name="ocean_heat_y_transport_due_to_bolus_advection"   unit="W"                         />
-    <field id="veiv_salttr"   long_name="ocean bolus salt transport along j-axis"       standard_name="ocean_salt_x_transport_due_to_bolus_advection"   unit="Kg"                        />
-    <field id="veiv_heattr3d" long_name="ocean bolus heat transport along j-axis"       standard_name="ocean_heat_y_transport_due_to_bolus_advection"   unit="W"    grid_ref="grid_V_3D" />
-    <field id="veiv_salttr3d" long_name="ocean bolus salt transport along j-axis"       standard_name="ocean_salt_y_transport_due_to_bolus_advection"   unit="kg"   grid_ref="grid_V_3D" />
+    <field id="voce_eiv"      long_name="EIV ocean current along j-axis"           standard_name="bolus_sea_water_y_velocity"                     unit="m/s"  grid_ref="grid_V_3D_inner" />
+    <field id="veiv_masstr"   long_name="EIV Ocean Mass Y Transport"               standard_name="bolus_ocean_mass_y_transport"                   unit="kg/s" grid_ref="grid_V_3D_inner" />
+    <field id="veiv_heattr"   long_name="ocean bolus heat transport along j-axis"  standard_name="ocean_heat_y_transport_due_to_bolus_advection"  unit="W"    grid_ref="grid_V_2D_inner" />
+    <field id="veiv_salttr"   long_name="ocean bolus salt transport along j-axis"  standard_name="ocean_salt_x_transport_due_to_bolus_advection"  unit="Kg"   grid_ref="grid_V_2D_inner" />
+    <field id="veiv_heattr3d" long_name="ocean bolus heat transport along j-axis"  standard_name="ocean_heat_y_transport_due_to_bolus_advection"  unit="W"    grid_ref="grid_V_3D_inner" />
+    <field id="veiv_salttr3d" long_name="ocean bolus salt transport along j-axis"  standard_name="ocean_salt_y_transport_due_to_bolus_advection"  unit="kg"   grid_ref="grid_V_3D_inner" />
 
 
     <!-- voce_bbl: available with ln_trabbl=T and nn_bbl_adv=1 -->
@@ -695,13 +694,13 @@ that are available in the tidal-forcing implementation (see
     <field id="vtbl"         long_name="meridional current in the Losh tbl"   unit="m/s" />
 
     <!-- variables available with diaar5 -->
-    <field id="v_masstr"      long_name="ocean eulerian mass transport along j-axis"    standard_name="ocean_mass_y_transport"                          unit="kg/s" grid_ref="grid_V_3D" />
+    <field id="v_masstr"      long_name="ocean eulerian mass transport along j-axis"    standard_name="ocean_mass_y_transport"                          unit="kg/s"       grid_ref="grid_V_3D_inner" />
     <field id="v_heattr"      long_name="ocean eulerian heat transport along j-axis"    standard_name="ocean_heat_y_transport"                          unit="W"          grid_ref="grid_V_2D_inner" />
     <field id="v_salttr"      long_name="ocean eulerian salt transport along i-axis"    standard_name="ocean_salt_y_transport"                          unit="1e-3*kg/s"  grid_ref="grid_V_2D_inner" />
-    <field id="vadv_heattr"   long_name="ocean advective heat transport along j-axis"   standard_name="advectice_ocean_heat_y_transport"                unit="W"                         />
-    <field id="vadv_salttr"   long_name="ocean advective salt transport along j-axis"   standard_name="advectice_ocean_salt_y_transport"                unit="1e-3*kg/s"                 />
-    <field id="vdiff_heattr"  long_name="ocean diffusion heat transport along j-axis"   standard_name="ocean_heat_y_transport_due_to_diffusion"         unit="W"                         />
-    <field id="vdiff_salttr"  long_name="ocean diffusion salt transport along j-axis"   standard_name="ocean_salt_y_transport_due_to_diffusion"         unit="1e-3*kg/s"                 />
+    <field id="vadv_heattr"   long_name="ocean advective heat transport along j-axis"   standard_name="advectice_ocean_heat_y_transport"                unit="W"          grid_ref="grid_V_2D_inner" />
+    <field id="vadv_salttr"   long_name="ocean advective salt transport along j-axis"   standard_name="advectice_ocean_salt_y_transport"                unit="1e-3*kg/s"  grid_ref="grid_V_2D_inner" />
+    <field id="vdiff_heattr"  long_name="ocean diffusion heat transport along j-axis"   standard_name="ocean_heat_y_transport_due_to_diffusion"         unit="W"          grid_ref="grid_V_2D_inner" />
+    <field id="vdiff_salttr"  long_name="ocean diffusion salt transport along j-axis"   standard_name="ocean_salt_y_transport_due_to_diffusion"         unit="1e-3*kg/s"  grid_ref="grid_V_2D_inner" />
   </field_group>
 
   <!-- W grid -->
@@ -712,32 +711,35 @@ that are available in the tidal-forcing implementation (see
     <field id="woce_e3w"     long_name="ocean vertical velocity * e3w"                                                                        unit="m2/s"  > woce * e3w </field>
     <field id="wocetr_eff"   long_name="effective ocean vertical transport"                                                                   unit="m3/s" />
 
-    <!-- woce_eiv: available with EIV  (ln_ldfeiv=T and ln_ldfeiv_dia=T)  -->
-    <field id="woce_eiv"     long_name="EIV ocean vertical velocity"                    standard_name="bolus_upward_sea_water_velocity"       unit="m/s"  />
-    <field id="weiv_masstr"  long_name="EIV Upward Ocean Mass Transport"  standard_name="bolus_upward_ocean_mass_transport"             unit="kg/s"   />
-    <field id="weiv_heattr3d" long_name="ocean bolus heat transport"    standard_name="ocean_heat_z_transport_due_to_bolus_advection"   unit="W"    />
-    <field id="weiv_salttr3d" long_name="ocean bolus salt transport"    standard_name="ocean_salt_z_transport_due_to_bolus_advection"   unit="kg"   />
+    <!-- variables available with WAVE (ln_wave=T) -->
+    <field id="wstokes"      long_name="Stokes Drift vertical velocity"                 standard_name="upward_StokesDrift_velocity"           unit="m/s"  />
 
-    <field id="avt"          long_name="vertical eddy diffusivity"                      standard_name="ocean_vertical_heat_diffusivity"       unit="m2/s" grid_ref="grid_W_3D_inner" />
-    <field id="avt_e3w"      long_name="vertical heat diffusivity * e3w"                unit="m3/s" > avt * e3w </field>
-    <field id="logavt"       long_name="logarithm of vertical eddy diffusivity"         standard_name="ocean_vertical_heat_diffusivity"       unit="m2/s" grid_ref="grid_W_3D_inner" />
-    <field id="avm"          long_name="vertical eddy viscosity"                        standard_name="ocean_vertical_momentum_diffusivity"   unit="m2/s" />
-    <field id="avm_e3w"      long_name="vertical eddy viscosity * e3w"   unit="m3/s" > avm * e3w </field>
+    <!-- woce_eiv: available with EIV  (ln_ldfeiv=T and ln_ldfeiv_dia=T)  -->
+    <field id="woce_eiv"     long_name="EIV ocean vertical velocity"                    standard_name="bolus_upward_sea_water_velocity"               unit="m/s"  grid_ref="grid_W_3D_inner" />
+    <field id="weiv_masstr"  long_name="EIV Upward Ocean Mass Transport"                standard_name="bolus_upward_ocean_mass_transport"             unit="kg/s" grid_ref="grid_W_3D_inner" />
+    <field id="weiv_heattr3d" long_name="ocean bolus heat transport"                    standard_name="ocean_heat_z_transport_due_to_bolus_advection" unit="W"    grid_ref="grid_W_3D_inner" />
+    <field id="weiv_salttr3d" long_name="ocean bolus salt transport"                    standard_name="ocean_salt_z_transport_due_to_bolus_advection" unit="kg"   grid_ref="grid_W_3D_inner" />
+
+    <!-- avt, avm -->
+    <field id="avt"          long_name="vertical eddy diffusivity"                      standard_name="ocean_vertical_heat_diffusivity"     unit="m2/s" grid_ref="grid_W_3D_inner" />
+    <field id="avt_e3w"      long_name="vertical heat diffusivity * e3w"                                                                    unit="m3/s"                             > avt * e3w </field>
+    <field id="logavt"       long_name="logarithm of vertical eddy diffusivity"         standard_name="ocean_vertical_heat_diffusivity"     unit="m2/s" grid_ref="grid_W_3D_inner" />
+    <field id="avm"          long_name="vertical eddy viscosity"                        standard_name="ocean_vertical_momentum_diffusivity" unit="m2/s"                            />
+    <field id="avm_e3w"      long_name="vertical eddy viscosity * e3w"                                                                      unit="m3/s"                             > avm * e3w </field>
 
     <!-- avs: /= avt with ln_zdfddm=T -->
-    <field id="avs"          long_name="salt vertical eddy diffusivity"                 standard_name="ocean_vertical_salt_diffusivity"       unit="m2/s" grid_ref="grid_W_3D_inner" />
-    <field id="avs_e3w"      long_name="vertical salt diffusivity * e3w"   unit="m3/s" > avs * e3w </field>
-    <field id="logavs"       long_name="logarithm of salt vertical eddy diffusivity"    standard_name="ocean_vertical_heat_diffusivity"       unit="m2/s" grid_ref="grid_W_3D_inner" />
+    <field id="avs"          long_name="salt vertical eddy diffusivity"                 standard_name="ocean_vertical_salt_diffusivity"     unit="m2/s" grid_ref="grid_W_3D_inner" />
+    <field id="avs_e3w"      long_name="vertical salt diffusivity * e3w"                                                                    unit="m3/s"                             > avs * e3w </field>
+    <field id="logavs"       long_name="logarithm of salt vertical eddy diffusivity"    standard_name="ocean_vertical_heat_diffusivity"     unit="m2/s" grid_ref="grid_W_3D_inner" />
 
     <!-- avt_evd and avm_evd: available with ln_zdfevd -->
-    <field id="avt_evd"      long_name="convective enhancement of vertical diffusivity" standard_name="ocean_vertical_tracer_diffusivity_due_to_convection"     unit="m2/s" />
-    <field id="avt_evd_e3w"  long_name="convective enhancement to vertical diffusivity * e3w "    unit="m3/s" > avt_evd * e3w </field>
-    <field id="avm_evd"      long_name="convective enhancement of vertical viscosity"   standard_name="ocean_vertical_momentum_diffusivity_due_to_convection"   unit="m2/s" />
+    <field id="avt_evd"      long_name="convective enhancement of vertical diffusivity" standard_name="ocean_vertical_tracer_diffusivity_due_to_convection"   unit="m2/s" grid_ref="grid_W_3D_inner" />
+    <field id="avt_evd_e3w"  long_name="convective enhancement to vertical diffusivity * e3w "                                                                unit="m3/s"                             > avt_evd * e3w </field>
+    <field id="avm_evd"      long_name="convective enhancement of vertical viscosity"   standard_name="ocean_vertical_momentum_diffusivity_due_to_convection" unit="m2/s" grid_ref="grid_W_3D_inner" />
 
     <!-- mf_app and mf_wp: available with ln_zdfmfc -->
-    <field id="mf_app"      long_name="convective area"        standard_name="mf_convective_area"    unit="%"      grid_ref="grid_W_3D" />
-    <field id="mf_wp"       long_name="convective velocity"    standard_name="mf_convective_velo"    unit="m/s"    grid_ref="grid_W_3D" />
-
+    <field id="mf_app"      long_name="convective area"        standard_name="mf_convective_area"    unit="%"      grid_ref="grid_W_3D_inner" />
+    <field id="mf_wp"       long_name="convective velocity"    standard_name="mf_convective_velo"    unit="m/s"    grid_ref="grid_W_3D_inner" />
 
     <!-- avt_tide: available with ln_zdfiwm=T -->
     <field id="av_ratio"     long_name="S over T diffusivity ratio"                     standard_name="salinity_over_temperature_diffusivity_ratio"                     unit="1"    grid_ref="grid_W_3D_inner" />
@@ -746,12 +748,9 @@ that are available in the tidal-forcing implementation (see
     <field id="pcmap_iwm"    long_name="power consumed by wave-driven mixing"           standard_name="vertically_integrated_power_consumption_by_wave_driven_mixing"   unit="W/m2" grid_ref="grid_W_2D_inner" />
     <field id="emix_iwm"     long_name="power density available for mixing"             standard_name="power_available_for_mixing_from_breaking_internal_waves"         unit="W/kg" grid_ref="grid_W_3D_inner" />
 
-    <!-- variables available with WAVE (ln_wave=T) -->
-    <field id="wstokes"      long_name="Stokes Drift vertical velocity"                 standard_name="upward_StokesDrift_velocity"   unit="m/s" />
-
     <!-- variables available with diaar5 -->
-    <field id="w_masstr"     long_name="vertical mass transport"                        standard_name="upward_ocean_mass_transport"             unit="kg/s"   />
-    <field id="w_masstr2"    long_name="square of vertical mass transport"              standard_name="square_of_upward_ocean_mass_transport"   unit="kg2/s2" />
+    <field id="w_masstr"     long_name="vertical mass transport"                        standard_name="upward_ocean_mass_transport"             unit="kg/s"   grid_ref="grid_W_3D_inner" />
+    <field id="w_masstr2"    long_name="square of vertical mass transport"              standard_name="square_of_upward_ocean_mass_transport"   unit="kg2/s2" grid_ref="grid_W_3D_inner" />
 
     <!-- EOS -->
     <field id="bn2"          long_name="squared Brunt-Vaisala frequency"                unit="s-2" />
@@ -767,15 +766,15 @@ that are available in the tidal-forcing implementation (see
 
   <!-- F grid -->
   <field_group id="grid_F" grid_ref="grid_F_2D">
-    <field id="e3f"          long_name="F-cell thickness"                      standard_name="cell_thickness"         unit="m"   grid_ref="grid_F_3D" />
-    <field id="e3f_0"        long_name="F-cell thickness"                      standard_name="cell_thickness"         unit="m"   grid_ref="grid_F_3D" />
-    <field id="hf"           long_name="water column height at F point"        standard_name="water_column_height_F"  unit="m"                     />
-    <field id="ssKEf"        long_name="surface kinetic energy at F point"     standard_name="specific_kinetic_energy_of_sea_water"   unit="m2/s2" grid_ref="grid_F_2D_inner" />
-    <field id="ssrelvor"     long_name="surface relative vorticity"            standard_name="relative_vorticity"     unit="1/s"       grid_ref="grid_F_2D_inner" />
-    <field id="ssplavor"     long_name="surface planetary vorticity"           standard_name="planetary_vorticity"    unit="1/s"       />
-    <field id="ssrelpotvor"  long_name="surface relative potential vorticity"  standard_name="relpot_vorticity"       unit="1/m.s"     grid_ref="grid_F_2D_inner" />
-    <field id="ssabspotvor"  long_name="surface absolute potential vorticity"  standard_name="abspot_vorticity"       unit="1/m.s"     grid_ref="grid_F_2D_inner" />
-    <field id="ssEns"        long_name="surface enstrophy"                     standard_name="enstrophy"              unit="1/m2.s2"   grid_ref="grid_F_2D_inner" />
+    <field id="e3f"          long_name="F-cell thickness"                      standard_name="cell_thickness"                       unit="m"       grid_ref="grid_F_3D"       />
+    <field id="e3f_0"        long_name="F-cell thickness"                      standard_name="cell_thickness"                       unit="m"       grid_ref="grid_F_3D"       />
+    <field id="hf"           long_name="water column height at F point"        standard_name="water_column_height_F"                unit="m"                                  />
+    <field id="ssKEf"        long_name="surface kinetic energy at F point"     standard_name="specific_kinetic_energy_of_sea_water" unit="m2/s2"   grid_ref="grid_F_2D_inner" />
+    <field id="ssrelvor"     long_name="surface relative vorticity"            standard_name="relative_vorticity"                   unit="1/s"     grid_ref="grid_F_2D_inner" />
+    <field id="ssplavor"     long_name="surface planetary vorticity"           standard_name="planetary_vorticity"                  unit="1/s"     grid_ref="grid_F_2D_inner" />
+    <field id="ssrelpotvor"  long_name="surface relative potential vorticity"  standard_name="relpot_vorticity"                     unit="1/m.s"   grid_ref="grid_F_2D_inner" />
+    <field id="ssabspotvor"  long_name="surface absolute potential vorticity"  standard_name="abspot_vorticity"                     unit="1/m.s"   grid_ref="grid_F_2D_inner" />
+    <field id="ssEns"        long_name="surface enstrophy"                     standard_name="enstrophy"                            unit="1/m2.s2" grid_ref="grid_F_2D_inner" />
   </field_group>
 
   <!-- AGRIF sponge -->
@@ -816,16 +815,16 @@ that are available in the tidal-forcing implementation (see
 
   <!-- transects -->
   <field_group id="oce_straits">
-    <field id="uoce_e3u_ave"         long_name="Monthly average of u*e3u"                        field_ref="uoce_e3u"                    freq_op="1mo"   freq_offset="_reset_" > @uoce_e3u </field>
-    <field id="uoce_e3u_ave_vsum"    long_name="Vertical sum of u*e3u"                           field_ref="uoce_e3u_ave"         grid_ref="grid_U_vsum"     />
+    <field id="uoce_e3u_ave"         long_name="Monthly average of u*e3u"                        field_ref="uoce_e3u"             freq_op="1mo"   freq_offset="_reset_" > @uoce_e3u </field>
+    <field id="uoce_e3u_ave_vsum"    long_name="Vertical sum of u*e3u"                           field_ref="uoce_e3u_ave"         grid_ref="grid_U_vsum"    />
     <field id="uocetr_vsum_section"  long_name="Total 2D transport in i-direction"               field_ref="uoce_e3u_ave_vsum"    grid_ref="grid_U_scalar"  detect_missing_value="true"> this * e2u </field>
     <field id="uocetr_strait"        long_name="Total transport across lines in i-direction"     field_ref="uocetr_vsum_section"  grid_ref="grid_U_4strait" />
     <field id="u_masstr_strait"      long_name="Sea water transport across line in i-direction"  field_ref="uocetr_strait"        grid_ref="grid_U_4strait_hsum" unit="kg/s"> this * maskMFO_u * $rho0 </field>
 
-    <field id="voce_e3v_ave"         long_name="Monthly average of v*e3v"                        field_ref="voce_e3v"                    freq_op="1mo"   freq_offset="_reset_" > @voce_e3v </field>
-    <field id="voce_e3v_ave_vsum"    long_name="Vertical sum of v*e3v"                           field_ref="voce_e3v_ave"         grid_ref="grid_V_vsum"      />
+    <field id="voce_e3v_ave"         long_name="Monthly average of v*e3v"                        field_ref="voce_e3v"             freq_op="1mo"   freq_offset="_reset_" > @voce_e3v </field>
+    <field id="voce_e3v_ave_vsum"    long_name="Vertical sum of v*e3v"                           field_ref="voce_e3v_ave"         grid_ref="grid_V_vsum"    />
     <field id="vocetr_vsum_section"  long_name="Total 2D transport of in j-direction"            field_ref="voce_e3v_ave_vsum"    grid_ref="grid_V_scalar"  detect_missing_value="true"> this * e1v </field>
-    <field id="vocetr_strait"        long_name="Total transport across lines in j-direction"     field_ref="vocetr_vsum_section"  grid_ref="grid_V_4strait"  />
+    <field id="vocetr_strait"        long_name="Total transport across lines in j-direction"     field_ref="vocetr_vsum_section"  grid_ref="grid_V_4strait" />
     <field id="v_masstr_strait"      long_name="Sea water transport across line in j-direction"  field_ref="vocetr_strait"        grid_ref="grid_V_4strait_hsum" unit="kg/s"> this * maskMFO_v * $rho0 </field>
 
     <field id="masstr_strait"        long_name="Sea water transport across line"                                                  grid_ref="grid_4strait"  > u_masstr_strait + v_masstr_strait </field>
@@ -1052,12 +1051,6 @@ that are available in the tidal-forcing implementation (see
     <field id="ketrd_convP2K" long_name="ke-trend: conversion (potential to kinetic)"      unit="W/s^3"                        />
     <field id="KE"            long_name="kinetic energy: u(n)*u(n+1)/2"                    unit="W/s^2"                        />
 
-    <!-- variables available when explicit lateral mixing is used (ln_dynldf_OFF=F) -->
-    <field id="dispkexyfo"    long_name="KE-trend: lateral  mixing induced dissipation"   standard_name="ocean_kinetic_energy_dissipation_per_unit_area_due_to_xy_friction"                   unit="W/m^2" grid_ref="grid_T_2D" />
-    <field id="dispkevfo"     long_name="KE-trend: vertical mixing induced dissipation"   standard_name="ocean_kinetic_energy_dissipation_per_unit_area_due_to_vertical_friction"             unit="W/m^2" grid_ref="grid_T_2D" />
-    <!-- variables available with ln_traadv_eiv=T and ln_diaeiv=T -->
-    <field id="eketrd_eiv"    long_name="EKE-trend due to parameterized eddy advection"   standard_name="tendency_of_ocean_eddy_kinetic_energy_content_due_to_parameterized_eddy_advection"   unit="W/m^2" grid_ref="grid_T_2D" />
-
     <!-- variables available with ln_PE_trd -->
     <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"                        />
diff --git a/src/ICE/icedyn_rdgrft.F90 b/src/ICE/icedyn_rdgrft.F90
index edf8f8215..0f6afe0fb 100644
--- a/src/ICE/icedyn_rdgrft.F90
+++ b/src/ICE/icedyn_rdgrft.F90
@@ -309,8 +309,9 @@ CONTAINS
       !! ** Method  :   Compute the thickness distribution of the ice and open water
       !!                participating in ridging and of the resulting ridges.
       !!-------------------------------------------------------------------
-      REAL(wp), DIMENSION(:)  , INTENT(in) ::   pato_i, pclosing_net
-      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pa_i, pv_i
+      REAL(wp), DIMENSION(:,:), INTENT(in)           ::   pa_i, pv_i
+      REAL(wp), DIMENSION(:)  , INTENT(in)           ::   pato_i
+      REAL(wp), DIMENSION(:)  , INTENT(in), OPTIONAL ::   pclosing_net
       !!
       INTEGER  ::   ji, jl                     ! dummy loop indices
       REAL(wp) ::   z1_gstar, z1_astar, zhmean, zfac   ! local scalar
@@ -511,39 +512,43 @@ CONTAINS
          END DO
       END DO
       !
-      ! 3) closing_gross
-      !-----------------
-      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.
-      ! NOTE: 0 < aksum <= 1
-      WHERE( zaksum(1:npti) > epsi10 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
-      ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp
-      END WHERE
-
-      ! correction to closing rate if excessive ice removal
-      !----------------------------------------------------
-      ! Reduce the closing rate if more than 100% of any ice category would be removed
-      ! Reduce the opening rate in proportion
-      DO jl = 1, jpl
+      IF( PRESENT( pclosing_net ) ) THEN
+         !
+         ! 3) closing_gross
+         !-----------------
+         ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.
+         ! NOTE: 0 < aksum <= 1
+         WHERE( zaksum(1:npti) > epsi10 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
+         ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp
+         END WHERE
+         
+         ! correction to closing rate if excessive ice removal
+         !----------------------------------------------------
+         ! Reduce the closing rate if more than 100% of any ice category would be removed
+         ! Reduce the opening rate in proportion
+         DO jl = 1, jpl
+            DO ji = 1, npti
+               zfac = apartf(ji,jl) * closing_gross(ji) * rDt_ice
+               IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
+                  closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_Dt_ice
+               ENDIF
+            END DO
+         END DO
+         
+         ! 4) correction to opening if excessive open water removal
+         !---------------------------------------------------------
+         ! Reduce the closing rate if more than 100% of the open water would be removed
+         ! Reduce the opening rate in proportion
          DO ji = 1, npti
-            zfac = apartf(ji,jl) * closing_gross(ji) * rDt_ice
-            IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
-               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_Dt_ice
+            zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice
+            IF( zfac < 0._wp ) THEN           ! would lead to negative ato_i
+               opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice
+            ELSEIF( zfac > zasum(ji) ) THEN   ! would lead to ato_i > asum
+               opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice
             ENDIF
          END DO
-      END DO
-
-      ! 4) correction to opening if excessive open water removal
-      !---------------------------------------------------------
-      ! Reduce the closing rate if more than 100% of the open water would be removed
-      ! Reduce the opening rate in proportion
-      DO ji = 1, npti
-         zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice
-         IF( zfac < 0._wp ) THEN           ! would lead to negative ato_i
-            opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice
-         ELSEIF( zfac > zasum(ji) ) THEN   ! would lead to ato_i > asum
-            opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice
-         ENDIF
-      END DO
+         !
+      ENDIF
       !
    END SUBROUTINE rdgrft_prep
 
@@ -907,7 +912,7 @@ CONTAINS
             CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i )
             CALL tab_2d_1d( npti, nptidx(1:npti), zstrength(1:npti)     , strength )
 
-            CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
+            CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d )
             !
             zaksum(1:npti) = apartf(1:npti,0) !clem: aksum should be defined in the header => local to module
             DO jl = 1, jpl
diff --git a/src/ICE/icedyn_rhg_eap.F90 b/src/ICE/icedyn_rhg_eap.F90
index 7f3ac9e9a..59fb26ffd 100644
--- a/src/ICE/icedyn_rhg_eap.F90
+++ b/src/ICE/icedyn_rhg_eap.F90
@@ -125,12 +125,12 @@ CONTAINS
       !!              Bouillon et al., Ocean Modelling 2013
       !!              Kimmritz et al., Ocean Modelling 2016 & 2017
       !!-------------------------------------------------------------------
-      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
-      INTEGER                 , INTENT(in   ) ::   Kmm                                   ! ocean time level index
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
-      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   paniso_11 , paniso_12                 ! structure tensor components
-      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   prdg_conv                             ! for ridging
+      INTEGER                    , INTENT(in   ) ::   kt                                    ! time step
+      INTEGER                    , INTENT(in   ) ::   Kmm                                   ! ocean time level index
+      REAL(wp), DIMENSION(:,:)   , INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
+      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
+      REAL(wp), DIMENSION(:,:)   , INTENT(inout) ::   paniso_11 , paniso_12                 ! structure tensor components
+      REAL(wp), DIMENSION(:,:)   , INTENT(inout) ::   prdg_conv                             ! for ridging
       !!
       INTEGER ::   ji, jj       ! dummy loop indices
       INTEGER ::   jter         ! local integers
@@ -148,15 +148,15 @@ CONTAINS
       !
       REAL(wp) ::   zintb, zintn                                        ! dummy argument
       REAL(wp) ::   zfac_x, zfac_y
-      REAL(wp) ::   zshear, zdum1, zdum2
+      REAL(wp) ::   zdum1, zdum2
       REAL(wp) ::   zstressptmp, zstressmtmp, zstress12tmpF             ! anisotropic stress tensor components
       REAL(wp) ::   zalphar, zalphas                                    ! for mechanical redistribution
       REAL(wp) ::   zmresult11, zmresult12, z1dtevpkth, zp5kth, z1_dtevp_A  ! for structure tensor evolution
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zstress12tmp                    ! anisotropic stress tensor component for regridding
-      REAL(wp), DIMENSION(jpi,jpj) ::   zyield11, zyield22, zyield12    ! yield surface tensor for history
+      REAL(wp), DIMENSION(A2D(0))  ::   zyield11, zyield22, zyield12    ! yield surface tensor for history
       REAL(wp), DIMENSION(jpi,jpj) ::   zdelta, zp_delt                 ! delta and P/delta at T points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! tension
+      REAL(wp), DIMENSION(A2D(0))  ::   zten_i, zshear                  ! tension, shear
       REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
@@ -178,7 +178,6 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast)
       REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
       !
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00, zmsk15
       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence
 
@@ -186,7 +185,8 @@ CONTAINS
       REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small
       REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
       !! --- check convergence
-      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice
+      REAL(wp), DIMENSION(A2D(0))  ::   zmsk00, zmsk15
+      REAL(wp), DIMENSION(A2D(0))  ::   zu_ice, zv_ice
       !! --- diags
       REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p
@@ -202,11 +202,11 @@ CONTAINS
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology'
       !
       ! for diagnostics and convergence tests
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( 0, 0, 0, 0 )
          zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
       END_2D
       IF( nn_rhg_chkcvg > 0 ) THEN
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 0, 0, 0, 0 )
             zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
          END_2D
       ENDIF
@@ -283,7 +283,14 @@ CONTAINS
       !    non-embedded sea ice: use ocean surface for slope calculation
       zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
 
-      DO_2D( 0, 0, 0, 0 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         zm1          = ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) )  ! Ice/snow mass at U-V points
+!!$         zm1          = ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) + rhow * (vt_ip(ji,jj) + vt_il(ji,jj)) ) ! clem: this should replace the above
+         zmf  (ji,jj) = zm1 * ff_t(ji,jj)                            ! Coriolis at T points (m*f)
+         zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin )                   ! dt/m at T points (for alpha and beta coefficients)
+      END_2D
+
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
 
          ! ice fraction at U-V points
          zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
@@ -291,8 +298,11 @@ CONTAINS
 
          ! Ice/snow mass at U-V points
          zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
+!!$         zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) + rhow * (vt_ip(ji  ,jj  ) + vt_il(ji  ,jj  )) ) ! clem: this should replace the above
          zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
+!!$         zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) + rhow * (vt_ip(ji+1,jj  ) + vt_il(ji+1,jj  )) ) ! clem: this should replace the above
          zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
+!!$         zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) + rhow * (vt_ip(ji  ,jj+1) + vt_il(ji  ,jj+1)) ) ! clem: this should replace the above
          zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
          zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
 
@@ -300,12 +310,6 @@ CONTAINS
          v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1)
          u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
 
-         ! Coriolis at T points (m*f)
-         zmf(ji,jj)      = zm1 * ff_t(ji,jj)
-
-         ! dt/m at T points (for alpha and beta coefficients)
-         zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
-
          ! m/dt
          zmU_t(ji,jj)    = zmassU * z1_dtevp
          zmV_t(ji,jj)    = zmassV * z1_dtevp
@@ -333,12 +337,11 @@ CONTAINS
          ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF
 
       END_2D
-      CALL lbc_lnk( 'icedyn_rhg_eap', zmf, 'T', 1.0_wp, zdt_m, 'T', 1.0_wp )
       !
       !                                  !== Landfast ice parameterization ==!
       !
       IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ! ice thickness at U-V points
             zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
             zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
@@ -354,7 +357,7 @@ CONTAINS
          END_2D
          !
       ELSE                               !-- no landfast
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ztaux_base(ji,jj) = 0._wp
             ztauy_base(ji,jj) = 0._wp
          END_2D
@@ -371,14 +374,14 @@ CONTAINS
          !
          ! convergence test
          IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2  ) THEN
-            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            DO_2D( 0, 0, 0, 0 )
                zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step
                zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1)
             END_2D
          ENDIF
 
          ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
 
             ! shear at F points
             zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
@@ -409,15 +412,13 @@ CONTAINS
             ! delta at T points
             zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
 
-         END_2D
-         CALL lbc_lnk( 'icedyn_rhg_eap', zdelta, 'T', 1.0_wp )
-
-         ! P/delta at T points
-         DO_2D( 1, 1, 1, 1 )
+            ! P/delta at T points
             zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
+            
          END_2D
+         CALL lbc_lnk( 'icedyn_rhg_eap', zdelta, 'T', 1.0_wp, zp_delt, 'T', 1.0_wp )
 
-         DO_2D( 0, 1, 0, 1 )   ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
+         DO_2D( 0, 0, 0, 0 )
 
              ! shear at T points
             zdsT = ( zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
@@ -470,16 +471,17 @@ CONTAINS
             zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zstressptmp ) * z1_alph1
             zs2(ji,jj) = ( zs2(ji,jj) * zalph1 + zstressmtmp ) * z1_alph1
          END_2D
-         CALL lbc_lnk( 'icedyn_rhg_eap', zstress12tmp, 'T', 1.0_wp , paniso_11, 'T', 1.0_wp , paniso_12, 'T', 1.0_wp)
+         CALL lbc_lnk( 'icedyn_rhg_eap', zstress12tmp, 'T', 1.0_wp , paniso_11, 'T', 1.0_wp , paniso_12, 'T', 1.0_wp, &
+            &                            zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp )
 
         ! Save beta at T-points for further computations
          IF( ln_aEVP ) THEN
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
                zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
             END_2D
          ENDIF
 
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
             ! stress12tmp at F points
             zstress12tmpF = ( zstress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + zstress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  &
                &            + zstress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + zstress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  &
@@ -498,10 +500,9 @@ CONTAINS
             zs12(ji,jj) = ( zs12(ji,jj) * zalph1 + zstress12tmpF ) * z1_alph1
 
          END_2D
-         CALL lbc_lnk( 'icedyn_rhg_eap', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )
 
          ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             !                   !--- U points
             zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
                &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
@@ -529,7 +530,7 @@ CONTAINS
          !  Bouillon et al. 2009 (eq 34-35) => stable
          IF( MOD(jter,2) == 0 ) THEN ! even iterations
             !
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                !                 !--- tau_io/(v_oce - v_ice)
                zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
                   &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
@@ -573,13 +574,7 @@ CONTAINS
                      &            )   * zmsk00y(ji,jj)
                ENDIF
             END_2D
-            CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1.0_wp )
-            !
-#if defined key_agrif
-!!          CALL agrif_interp_ice( 'V', jter, nn_nevp )
-            CALL agrif_interp_ice( 'V' )
-#endif
-            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
+            IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1.0_wp )
             !
             DO_2D( 0, 0, 0, 0 )
                !                 !--- tau_io/(u_oce - u_ice)
@@ -625,17 +620,13 @@ CONTAINS
                      &            )   * zmsk00x(ji,jj)
                ENDIF
             END_2D
-            CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp )
-            !
-#if defined key_agrif
-!!          CALL agrif_interp_ice( 'U', jter, nn_nevp )
-            CALL agrif_interp_ice( 'U' )
-#endif
-            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
+            IF( nn_hls == 1 ) THEN   ;   CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp )
+            ELSE                     ;   CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp )
+            ENDIF
             !
          ELSE ! odd iterations
             !
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                !                 !--- tau_io/(u_oce - u_ice)
                zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
                   &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
@@ -679,13 +670,7 @@ CONTAINS
                      &            )   * zmsk00x(ji,jj)
                ENDIF
             END_2D
-            CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp )
-            !
-#if defined key_agrif
-!!          CALL agrif_interp_ice( 'U', jter, nn_nevp )
-            CALL agrif_interp_ice( 'U' )
-#endif
-            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
+            IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp )
             !
             DO_2D( 0, 0, 0, 0 )
                !                 !--- tau_io/(v_oce - v_ice)
@@ -731,15 +716,19 @@ CONTAINS
                      &            )   * zmsk00y(ji,jj)
                ENDIF
             END_2D
-            CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1.0_wp )
+            IF( nn_hls == 1 ) THEN   ;   CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1.0_wp )
+            ELSE                     ;   CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp )
+            ENDIF
             !
+         ENDIF
 #if defined key_agrif
-!!          CALL agrif_interp_ice( 'V', jter, nn_nevp )
-            CALL agrif_interp_ice( 'V' )
+!!       CALL agrif_interp_ice( 'U', jter, nn_nevp )
+!!       CALL agrif_interp_ice( 'V', jter, nn_nevp )
+         CALL agrif_interp_ice( 'U' )
+         CALL agrif_interp_ice( 'V' )
 #endif
-            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
-            !
-         ENDIF
+         IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
+         IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
 
          ! convergence test
          IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg_eap( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 )
@@ -754,7 +743,7 @@ CONTAINS
       !------------------------------------------------------------------------------!
       ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
       !------------------------------------------------------------------------------!
-      DO_2D( 1, 0, 1, 0 )
+      DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
 
          ! shear at F points
          zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
@@ -778,22 +767,30 @@ CONTAINS
             &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
             &   ) * 0.25_wp * r1_e1e2t(ji,jj)
 
-         ! shear at T points
+         ! maximum shear rate at T points (includes tension, output only)
          pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
 
+        ! shear at T-points
+         zshear(ji,jj)   = SQRT( zds2 )
+
          ! divergence at T points
          pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
             &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
             &             ) * r1_e1e2t(ji,jj)
 
          ! delta at T points
-         zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta
-         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0
-         pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl
+         zdelta(ji,jj)   = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta
+
+         ! delta at T points
+         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0
+         pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * rswitch  
+                           ! it seems that deformation used for advection and mech redistribution is delta*
+                           ! MV in principle adding creep limit is a regularization for viscosity not for delta
+                           ! delta_star should not (in my view) be used in a replacement for delta
 
       END_2D
       CALL lbc_lnk( 'icedyn_rhg_eap', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp, &
-         &                              zten_i, 'T', 1.0_wp, zs1    , 'T', 1.0_wp, zs2     , 'T', 1.0_wp, &
+         &                                 zs1, 'T', 1.0_wp, zs2    , 'T', 1.0_wp, &
          &                                zs12, 'F', 1.0_wp )
 
       ! --- Store the stress tensor for the next time step --- !
@@ -806,42 +803,38 @@ CONTAINS
       ! 5) diagnostics
       !------------------------------------------------------------------------------!
       ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
-      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
-         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
-         !
-         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
-         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
-         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
-         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
-         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
-         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
-      ENDIF
+      IF( iom_use('utau_oi') )   CALL iom_put( 'utau_oi' , ztaux_oi(A2D(0)) * zmsk00 )
+      IF( iom_use('vtau_oi') )   CALL iom_put( 'vtau_oi' , ztauy_oi(A2D(0)) * zmsk00 )
+      IF( iom_use('utau_ai') )   CALL iom_put( 'utau_ai' , ztaux_ai(A2D(0)) * zmsk00 )
+      IF( iom_use('vtau_ai') )   CALL iom_put( 'vtau_ai' , ztauy_ai(A2D(0)) * zmsk00 )
+      IF( iom_use('utau_bi') )   CALL iom_put( 'utau_bi' , ztaux_bi(A2D(0)) * zmsk00 )
+      IF( iom_use('vtau_bi') )   CALL iom_put( 'vtau_bi' , ztauy_bi(A2D(0)) * zmsk00 )
 
       ! --- divergence, shear and strength --- !
-      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
-      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
-      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
-      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
+      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i (A2D(0)) * zmsk00 )   ! divergence
+      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i(A2D(0)) * zmsk00 )   ! shear
+      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength(A2D(0)) * zmsk00 )   ! strength
+      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , zdelta  (A2D(0)) * zmsk00 )   ! delta
 
       ! --- Stress tensor invariants (SIMIP diags) --- !
       IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
          !
-         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
+         ALLOCATE( zsig_I(A2D(0)) , zsig_II(A2D(0)) )
          !
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 0, 0, 0, 0 )
 
             ! Ice stresses
             ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
             ! These are NOT stress tensor components, neither stress invariants, neither stress principal components
             ! I know, this can be confusing...
-            zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl )
-            zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
+            zfac             =   strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )          ! viscosity
+            zsig1            =   zfac * ( pdivu_i(ji,jj) - zdelta(ji,jj) )
             zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
-            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
+            zsig12           =   zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp
 
             ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
-            zsig_I (ji,jj)   =   zsig1 * 0.5_wp                                      ! 1st stress invariant, aka average normal stress, aka negative pressure
-            zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )  ! 2nd  ''       ''    , aka maximum shear stress
+            zsig_I (ji,jj)   =   0.5_wp * zsig1 
+            zsig_II(ji,jj)   =   0.5_wp * SQRT ( zsig2 * zsig2 + 4._wp * zsig12 * zsig12 )
 
          END_2D
          !
@@ -859,21 +852,20 @@ CONTAINS
       ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
       IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
          !
-         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
+         ALLOCATE( zsig1_p(A2D(0)) , zsig2_p(A2D(0)) , zsig_I(A2D(0)) , zsig_II(A2D(0)) )
          !
-         DO_2D( 1, 1, 1, 1 )
+        DO_2D( 0, 0, 0, 0 )
 
-            ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
-            !                        and **deformations** at current iterates
+            ! For EVP solvers, ice stresses at current iterates can be used
             !                        following Lemieux & Dupont (2020)
-            zfac             =   zp_delt(ji,jj)
-            zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) )
+            zfac             =   strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
+            zsig1            =   zfac * ( pdivu_i(ji,jj) - zdelta(ji,jj) )
             zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
-            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
+            zsig12           =   zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp
 
             ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
-            zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                      ! 1st stress invariant, aka average normal stress, aka negative pressure
-            zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )  ! 2nd  ''       ''    , aka maximum shear stress
+            zsig_I(ji,jj)    =   0.5_wp * zsig1                                         ! normal stress
+            zsig_II(ji,jj)   =   0.5_wp * SQRT ( zsig2 * zsig2 + 4._wp * zsig12 * zsig12 ) ! max shear stress
 
             ! Normalized  principal stresses (used to display the ellipse)
             z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) )
@@ -881,8 +873,8 @@ CONTAINS
             zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
          END_2D
          !
-         CALL iom_put( 'sig1_pnorm' , zsig1_p )
-         CALL iom_put( 'sig2_pnorm' , zsig2_p )
+         CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00 )
+         CALL iom_put( 'sig2_pnorm' , zsig2_p(:,:) * zmsk00 )
 
          DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
 
@@ -901,21 +893,18 @@ CONTAINS
       ENDIF
 
       ! --- SIMIP --- !
-      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
-         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
-         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
-         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
-         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
-         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
-         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
-         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
-      ENDIF
+      IF( iom_use('dssh_dx') )   CALL iom_put( 'dssh_dx' , zspgU(A2D(0)) * zmsk00 )   ! Sea-surface tilt term in force balance (x)
+      IF( iom_use('dssh_dy') )   CALL iom_put( 'dssh_dy' , zspgV(A2D(0)) * zmsk00 )   ! Sea-surface tilt term in force balance (y)
+      IF( iom_use('corstrx') )   CALL iom_put( 'corstrx' , zCorU(A2D(0)) * zmsk00 )   ! Coriolis force term in force balance (x)
+      IF( iom_use('corstry') )   CALL iom_put( 'corstry' , zCorV(A2D(0)) * zmsk00 )   ! Coriolis force term in force balance (y)
+      IF( iom_use('intstrx') )   CALL iom_put( 'intstrx' , zfU  (A2D(0)) * zmsk00 )   ! Internal force term in force balance (x)
+      IF( iom_use('intstry') )   CALL iom_put( 'intstry' , zfV  (A2D(0)) * zmsk00 )   ! Internal force term in force balance (y)
 
       IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
          & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
          !
-         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
-            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
+         ALLOCATE( zdiag_xmtrp_ice(A2D(0)) , zdiag_ymtrp_ice(A2D(0)) , &
+            &      zdiag_xmtrp_snw(A2D(0)) , zdiag_ymtrp_snw(A2D(0)) , zdiag_xatrp(A2D(0)) , zdiag_yatrp(A2D(0)) )
          !
          DO_2D( 0, 0, 0, 0 )
             ! 2D ice mass, snow mass, area transport arrays (X, Y)
@@ -949,11 +938,11 @@ CONTAINS
       IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN
          IF( iom_use('uice_cvg') ) THEN
             IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
-               CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , &
-                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) )
+               CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(A2D(0)) - zu_ice(:,:) ) * zbeta(A2D(0)) * umask(A2D(0),1) , &
+                  &                           ABS( v_ice(A2D(0)) - zv_ice(:,:) ) * zbeta(A2D(0)) * vmask(A2D(0),1) ) * zmsk15(:,:) )
             ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
-               CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &
-                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
+               CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(A2D(0)) - zu_ice(:,:) ) * umask(A2D(0),1) , &
+                  &                                             ABS( v_ice(A2D(0)) - zv_ice(:,:) ) * vmask(A2D(0),1) ) * zmsk15(:,:) )
             ENDIF
          ENDIF
       ENDIF
@@ -974,22 +963,27 @@ CONTAINS
       !!
       !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)
       !!----------------------------------------------------------------------
-      INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index
-      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb          ! now and before velocities
-      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pmsk15
+      INTEGER ,                    INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index
+      REAL(wp), DIMENSION(:,:)   , INTENT(in) ::   pu, pv                    ! now velocities
+      REAL(wp), DIMENSION(A2D(0)), INTENT(in) ::   pub, pvb                  ! before velocities
+      REAL(wp), DIMENSION(A2D(0)), INTENT(in) ::   pmsk15
       !!
       INTEGER           ::   it, idtime, istatus
       INTEGER           ::   ji, jj          ! dummy loop indices
       REAL(wp)          ::   zresm           ! local real
       CHARACTER(len=20) ::   clname
+      LOGICAL           ::   ll_maxcvg
+      REAL(wp), DIMENSION(A2D(0),2) ::   zres
+      REAL(wp), DIMENSION(2)        ::   ztmp
       !!----------------------------------------------------------------------
-
+      ll_maxcvg = .FALSE.
+      !
       ! create file
       IF( kt == nit000 .AND. kiter == 1 ) THEN
          !
          IF( lwp ) THEN
             WRITE(numout,*)
-            WRITE(numout,*) 'rhg_cvg_eap : ice rheology convergence control'
+            WRITE(numout,*) 'rhg_cvg : ice rheology convergence control'
             WRITE(numout,*) '~~~~~~~'
          ENDIF
          !
@@ -998,7 +992,7 @@ CONTAINS
             IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
             istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
             istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime )
-            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid )
+            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid )
             istatus = NF90_ENDDEF(ncvgid)
          ENDIF
          !
@@ -1012,11 +1006,21 @@ CONTAINS
          zresm = 0._wp
       ELSE
          zresm = 0._wp
-         DO_2D( 0, 0, 0, 0 )
-            zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
-               &                     ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) )
-         END_2D
-         CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
+         IF( ll_maxcvg ) THEN   ! error max over the domain
+            DO_2D( 0, 0, 0, 0 )
+               zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
+                  &                     ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) )
+            END_2D
+            CALL mpp_max( 'icedyn_rhg_eap', zresm )
+         ELSE                   ! error averaged over the domain
+            DO_2D( 0, 0, 0, 0 )
+               zres(ji,jj,1) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
+                  &                 ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj)
+               zres(ji,jj,2) = pmsk15(ji,jj)
+            END_2D
+            ztmp(:) = glob_sum_vec( 'icedyn_rhg_eap', zres )
+            IF( ztmp(2) /= 0._wp )   zresm = ztmp(1) / ztmp(2)
+         ENDIF
       ENDIF
 
       IF( lwm ) THEN
diff --git a/src/OCE/DIA/diaar5.F90 b/src/OCE/DIA/diaar5.F90
index da501c5c1..a22ad1640 100644
--- a/src/OCE/DIA/diaar5.F90
+++ b/src/OCE/DIA/diaar5.F90
@@ -53,7 +53,7 @@ CONTAINS
       INTEGER :: dia_ar5_alloc
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( thick0(jpi,jpj) , sn0(jpi,jpj,jpk), STAT=dia_ar5_alloc )
+      ALLOCATE( thick0(A2D(0)) , sn0(A2D(0),jpk), STAT=dia_ar5_alloc )
       !
       CALL mpp_sum ( 'diaar5', dia_ar5_alloc )
       IF( dia_ar5_alloc /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_alloc: failed to allocate arrays' )
@@ -75,9 +75,10 @@ CONTAINS
       REAL(wp) ::   zvolssh, zvol, zssh_steric, zztmp, zarho, ztemp, zsal, zmass, zsst
       REAL(wp) ::   zaw, zbw, zrw
       !
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z2d, zpe                   ! 2D workspace
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: z3d, zrhd, ztpot, zgdept   ! 3D workspace (zgdept: needed to use the substitute)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zarea_ssh             ! 2D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z2d                   ! 2D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: z3d
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: zrhd, zgdept          ! 3D workspace (zgdept: needed to use the substitute)
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace
       !!--------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('dia_ar5')
@@ -85,10 +86,12 @@ CONTAINS
       IF( kt == nit000 )     CALL dia_ar5_init
 
       IF( l_ar5 ) THEN
-         ALLOCATE( zarea_ssh(jpi,jpj), zbotpres(jpi,jpj), z2d(jpi,jpj) )
+         ALLOCATE( zarea_ssh(A2D(0)), z2d(A2D(0)), z3d(A2D(0),jpk) )
          ALLOCATE( zrhd(jpi,jpj,jpk) )
          ALLOCATE( ztsn(jpi,jpj,jpk,jpts) )
-         zarea_ssh(:,:) = e1e2t(:,:) * ssh(:,:,Kmm)
+         zarea_ssh(:,:) = e1e2t(A2D(0)) * ssh(A2D(0),Kmm)
+         ztsn(:,:,:,:) = 0._wp
+         zrhd(:,:,:) = 0._wp
       ENDIF
       !
       CALL iom_put( 'e2u'      , e2u  (:,:) )
@@ -96,19 +99,19 @@ CONTAINS
       CALL iom_put( 'areacello', e1e2t(:,:) )
       !
       IF( iom_use( 'volcello' ) .OR. iom_use( 'masscello' )  ) THEN
-         zrhd(:,:,jpk) = 0._wp        ! ocean volume ; rhd is used as workspace
-         DO jk = 1, jpkm1
-            zrhd(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
-         END DO
-         DO jk = 1, jpk
-            z3d(:,:,jk) =  rho0 * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
-         END DO
-         CALL iom_put( 'volcello'  , zrhd(:,:,:)  )  ! WARNING not consistent with CMIP DR where volcello is at ca. 2000
+         z3d(:,:,jpk) = 0._wp        ! ocean volume ; rhd is used as workspace
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            z3d(ji,jj,jk) = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
+         END_3D
+         CALL iom_put( 'volcello'  , z3d(:,:,:)  )  ! WARNING not consistent with CMIP DR where volcello is at ca. 2000
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            z3d(ji,jj,jk) =  rho0 * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
+         END_3D
          CALL iom_put( 'masscello' , z3d (:,:,:) )   ! ocean mass
       ENDIF
       !
       IF( iom_use( 'e3tb' ) )  THEN    ! bottom layer thickness
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 0, 0, 0, 0 )
             ikb = mbkt(ji,jj)
             z2d(ji,jj) = e3t(ji,jj,ikb,Kmm)
          END_2D
@@ -128,61 +131,63 @@ CONTAINS
 
       IF( iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) THEN
          !
-         ztsn(:,:,:,jp_tem) = ts(:,:,:,jp_tem,Kmm)                    ! thermosteric ssh
-         ztsn(:,:,:,jp_sal) = sn0(:,:,:)
+         ztsn(A2D(0),:,jp_tem) = ts(A2D(0),:,jp_tem,Kmm)                    ! thermosteric ssh
+         ztsn(A2D(0),:,jp_sal) = sn0(:,:,:)
          ALLOCATE( zgdept(jpi,jpj,jpk) )
-         DO jk = 1, jpk
-            zgdept(:,:,jk) = gdept(:,:,jk,Kmm)
-         END DO
-         CALL eos( ztsn, zrhd, zgdept)                       ! now in situ density using initial salinity
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+            zgdept(ji,jj,jk) = gdept(ji,jj,jk,Kmm)
+         END_3D
+         CALL eos( ztsn, zrhd, zgdept )                       ! now in situ density using initial salinity
          !
-         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
-         DO jk = 1, jpkm1
-            zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * zrhd(:,:,jk)
-         END DO
+         z2d(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * zrhd(ji,jj,jk)
+         END_3D
          IF( ln_linssh ) THEN
             IF( ln_isfcav ) THEN
-               DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+               DO_2D( 0, 0, 0, 0 )
                   iks = mikt(ji,jj)
-                  zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,iks) + riceload(ji,jj)
+                  z2d(ji,jj) = z2d(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,iks) + riceload(ji,jj)
                END_2D
             ELSE
-               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * zrhd(:,:,1)
+               DO_2D( 0, 0, 0, 0 )                    
+                  z2d(ji,jj) = z2d(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,1)
+               END_2D
             END IF
 !!gm
 !!gm   riceload should be added in both ln_linssh=T or F, no?
 !!gm
          END IF
          !
-         zarho = glob_sum( 'diaar5', e1e2t(:,:) * zbotpres(:,:) )
+         zarho = glob_sum( 'diaar5', e1e2t(A2D(0)) * z2d(:,:) )
          zssh_steric = - zarho / area_tot
          CALL iom_put( 'sshthster', zssh_steric )
 
          !                                         ! steric sea surface height
-         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
-         DO jk = 1, jpkm1
-            zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * rhd(:,:,jk)
-         END DO
+         z2d(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * rhd(ji,jj,jk)
+         END_3D
          IF( ln_linssh ) THEN
             IF ( ln_isfcav ) THEN
-               DO ji = 1,jpi
-                  DO jj = 1,jpj
-                     iks = mikt(ji,jj)
-                     zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * rhd(ji,jj,iks) + riceload(ji,jj)
-                  END DO
-               END DO
+               DO_2D( 0, 0, 0, 0 )
+                  iks = mikt(ji,jj)
+                  z2d(ji,jj) = z2d(ji,jj) + ssh(ji,jj,Kmm) * rhd(ji,jj,iks) + riceload(ji,jj)
+               END_2D
             ELSE
-               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * rhd(:,:,1)
+               DO_2D( 0, 0, 0, 0 )
+                  z2d(ji,jj) = z2d(ji,jj) + ssh(ji,jj,Kmm) * rhd(ji,jj,1)
+               END_2D
             END IF
          END IF
          !
-         zarho = glob_sum( 'diaar5', e1e2t(:,:) * zbotpres(:,:) )
+         zarho = glob_sum( 'diaar5', e1e2t(A2D(0)) * z2d(:,:) )
          zssh_steric = - zarho / area_tot
          CALL iom_put( 'sshsteric', zssh_steric )
          !                                         ! ocean bottom pressure
          zztmp = rho0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
-         zbotpres(:,:) = zztmp * ( zbotpres(:,:) + ssh(:,:,Kmm) + thick0(:,:) )
-         CALL iom_put( 'botpres', zbotpres )
+         z2d(:,:) = zztmp * ( z2d(:,:) + ssh(A2D(0),Kmm) + thick0(:,:) )
+         CALL iom_put( 'botpres', z2d )
          !
          DEALLOCATE( zgdept )
          !
@@ -191,7 +196,7 @@ CONTAINS
       IF( iom_use( 'masstot' ) .OR. iom_use( 'temptot' )  .OR. iom_use( 'saltot' )  ) THEN
           !                                         ! Mean density anomalie, temperature and salinity
           ztsn(:,:,:,:) = 0._wp                    ! ztsn(:,:,1,jp_tem/sal) is used here as 2D Workspace for temperature & salinity
-          DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+          DO_3D( 0, 0, 0, 0, 1, jpkm1 )
              zztmp = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)
              ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zztmp * ts(ji,jj,jk,jp_tem,Kmm)
              ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zztmp * ts(ji,jj,jk,jp_sal,Kmm)
@@ -199,16 +204,16 @@ CONTAINS
 
           IF( ln_linssh ) THEN
             IF( ln_isfcav ) THEN
-               DO ji = 1, jpi
-                  DO jj = 1, jpj
-                     iks = mikt(ji,jj)
-                     ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zarea_ssh(ji,jj) * ts(ji,jj,iks,jp_tem,Kmm)
-                     ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zarea_ssh(ji,jj) * ts(ji,jj,iks,jp_sal,Kmm)
-                  END DO
-               END DO
+               DO_2D( 0, 0, 0, 0 )
+                  iks = mikt(ji,jj)
+                  ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zarea_ssh(ji,jj) * ts(ji,jj,iks,jp_tem,Kmm)
+                  ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zarea_ssh(ji,jj) * ts(ji,jj,iks,jp_sal,Kmm)
+               END_2D
             ELSE
-               ztsn(:,:,1,jp_tem) = ztsn(:,:,1,jp_tem) + zarea_ssh(:,:) * ts(:,:,1,jp_tem,Kmm)
-               ztsn(:,:,1,jp_sal) = ztsn(:,:,1,jp_sal) + zarea_ssh(:,:) * ts(:,:,1,jp_sal,Kmm)
+               DO_2D( 0, 0, 0, 0 )
+                  ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zarea_ssh(ji,jj) * ts(ji,jj,1,jp_tem,Kmm)
+                  ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zarea_ssh(ji,jj) * ts(ji,jj,1,jp_sal,Kmm)
+               END_2D
             END IF
          ENDIF
          !
@@ -226,37 +231,35 @@ CONTAINS
          IF( iom_use( 'toce_pot') .OR. iom_use( 'temptot_pot' ) .OR. iom_use( 'sst_pot' )  &
                                   .OR. iom_use( 'ssttot' ) .OR.  iom_use( 'tosmint_pot' ) ) THEN
             !
-            ALLOCATE( ztpot(jpi,jpj,jpk) )
-            ztpot(:,:,jpk) = 0._wp
+            z3d(:,:,jpk) = 0._wp
             DO jk = 1, jpkm1
-               ztpot(:,:,jk) = eos_pt_from_ct( ts(:,:,jk,jp_tem,Kmm), ts(:,:,jk,jp_sal,Kmm) )
+               z3d(:,:,jk) = eos_pt_from_ct( ts(A2D(0),jk,jp_tem,Kmm), ts(A2D(0),jk,jp_sal,Kmm) )
             END DO
             !
-            CALL iom_put( 'toce_pot', ztpot(:,:,:) )  ! potential temperature (TEOS-10 case)
-            CALL iom_put( 'sst_pot' , ztpot(:,:,1) )  ! surface temperature
+            CALL iom_put( 'toce_pot', z3d(:,:,:) )  ! potential temperature (TEOS-10 case)
+            CALL iom_put( 'sst_pot' , z3d(:,:,1) )  ! surface temperature
             !
             IF( iom_use( 'temptot_pot' ) ) THEN   ! Output potential temperature in case we use TEOS-10
                z2d(:,:) = 0._wp
-               DO jk = 1, jpkm1
-                 z2d(:,:) = z2d(:,:) + e1e2t(:,:) * e3t(:,:,jk,Kmm) * ztpot(:,:,jk)
-               END DO
+               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+                 z2d(ji,jj) = z2d(ji,jj) + e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk)
+               END_3D
                ztemp = glob_sum( 'diaar5', z2d(:,:)  )
                CALL iom_put( 'temptot_pot', ztemp / zvol )
              ENDIF
              !
              IF( iom_use( 'ssttot' ) ) THEN   ! Output potential temperature in case we use TEOS-10
-               zsst = glob_sum( 'diaar5',  e1e2t(:,:) * ztpot(:,:,1)  )
+               zsst = glob_sum( 'diaar5',  e1e2t(A2D(0)) * z3d(:,:,1)  )
                CALL iom_put( 'ssttot', zsst / area_tot )
              ENDIF
              ! Vertical integral of temperature
              IF( iom_use( 'tosmint_pot') ) THEN
                z2d(:,:) = 0._wp
-               DO_3D( 1, 1, 1, 1, 1, jpkm1 )
-                  z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) *  ztpot(ji,jj,jk)
+               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+                  z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) *  z3d(ji,jj,jk)
                END_3D
                CALL iom_put( 'tosmint_pot', z2d )
             ENDIF
-            DEALLOCATE( ztpot )
         ENDIF
       ELSE
          IF( iom_use('ssttot') ) THEN   ! Output sst in case we use EOS-80
@@ -269,8 +272,7 @@ CONTAINS
         ! Work done against stratification by vertical mixing
         ! Exclude points where rn2 is negative as convection kicks in here and
         ! work is not being done against stratification
-         ALLOCATE( zpe(A2D(0)) )
-         zpe(:,:) = 0._wp
+         z2d(:,:) = 0._wp
          IF( ln_zdfddm ) THEN
             DO_3D( 0, 0, 0, 0, 2, jpk )
                IF( rn2(ji,jj,jk) > 0._wp ) THEN
@@ -279,23 +281,22 @@ CONTAINS
                   zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw
                   zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
                   !
-                  zpe(ji, jj) = zpe(ji,jj)   &
+                  z2d(ji, jj) = z2d(ji,jj)   &
                      &        -  grav * (  avt(ji,jj,jk) * zaw * (ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) )  &
                      &                   - avs(ji,jj,jk) * zbw * (ts(ji,jj,jk-1,jp_sal,Kmm) - ts(ji,jj,jk,jp_sal,Kmm) ) )
                ENDIF
             END_3D
           ELSE
             DO_3D( 0, 0, 0, 0, 1, jpk )
-               zpe(ji,jj) = zpe(ji,jj) + avt(ji,jj,jk) * MIN(0._wp,rn2(ji,jj,jk)) * rho0 * e3w(ji,jj,jk,Kmm)
+               z2d(ji,jj) = z2d(ji,jj) + avt(ji,jj,jk) * MIN(0._wp,rn2(ji,jj,jk)) * rho0 * e3w(ji,jj,jk,Kmm)
             END_3D
          ENDIF
-          CALL iom_put( 'tnpeo', zpe )
-          DEALLOCATE( zpe )
+          CALL iom_put( 'tnpeo', z2d )
       ENDIF
 
       IF( l_ar5 ) THEN
-        DEALLOCATE( zarea_ssh , zbotpres, z2d )
-        DEALLOCATE( ztsn                 )
+        DEALLOCATE( zarea_ssh , z2d, z3d )
+        DEALLOCATE( ztsn )
       ENDIF
       !
       IF( ln_timing )   CALL timing_stop('dia_ar5')
@@ -316,33 +317,33 @@ CONTAINS
       REAL(wp), DIMENSION(A2D(nn_hls),jpk)    , INTENT(in)   :: pvflx  ! v-flux of advection/diffusion
       !
       INTEGER    ::  ji, jj, jk
-      REAL(wp), DIMENSION(A2D(nn_hls))  :: z2d
+      REAL(wp), DIMENSION(A2D(0))  :: z2d
       !!----------------------------------------------------------------------
 
-      z2d(:,:) = puflx(:,:,1)
+      z2d(:,:) = 0._wp
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          z2d(ji,jj) = z2d(ji,jj) + puflx(ji,jj,jk)
       END_3D
 
       IF( cptr == 'adv' ) THEN
-         IF( ktra == jp_tem ) CALL iom_put( 'uadv_heattr' , rho0_rcp * z2d(:,:) )  ! advective heat transport in i-direction
-         IF( ktra == jp_sal ) CALL iom_put( 'uadv_salttr' , rho0     * z2d(:,:) )  ! advective salt transport in i-direction
+         IF( ktra == jp_tem )   CALL iom_put( 'uadv_heattr'  , rho0_rcp * z2d(:,:) ) ! advective heat transport in i-direction
+         IF( ktra == jp_sal )   CALL iom_put( 'uadv_salttr'  , rho0     * z2d(:,:) ) ! advective salt transport in i-direction
       ELSE IF( cptr == 'ldf' ) THEN
-         IF( ktra == jp_tem ) CALL iom_put( 'udiff_heattr' , rho0_rcp * z2d(:,:) ) ! diffusive heat transport in i-direction
-         IF( ktra == jp_sal ) CALL iom_put( 'udiff_salttr' , rho0     * z2d(:,:) ) ! diffusive salt transport in i-direction
+         IF( ktra == jp_tem )   CALL iom_put( 'udiff_heattr' , rho0_rcp * z2d(:,:) ) ! diffusive heat transport in i-direction
+         IF( ktra == jp_sal )   CALL iom_put( 'udiff_salttr' , rho0     * z2d(:,:) ) ! diffusive salt transport in i-direction
       ENDIF
       !
-      z2d(:,:) = pvflx(:,:,1)
+      z2d(:,:) = 0._wp
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          z2d(ji,jj) = z2d(ji,jj) + pvflx(ji,jj,jk)
       END_3D
 
       IF( cptr == 'adv' ) THEN
-         IF( ktra == jp_tem ) CALL iom_put( 'vadv_heattr' , rho0_rcp * z2d(:,:) )  ! advective heat transport in j-direction
-         IF( ktra == jp_sal ) CALL iom_put( 'vadv_salttr' , rho0     * z2d(:,:) )  ! advective salt transport in j-direction
+         IF( ktra == jp_tem )   CALL iom_put( 'vadv_heattr'  , rho0_rcp * z2d(:,:) ) ! advective heat transport in j-direction
+         IF( ktra == jp_sal )   CALL iom_put( 'vadv_salttr'  , rho0     * z2d(:,:) ) ! advective salt transport in j-direction
       ELSE IF( cptr == 'ldf' ) THEN
-         IF( ktra == jp_tem ) CALL iom_put( 'vdiff_heattr' , rho0_rcp * z2d(:,:) ) ! diffusive heat transport in j-direction
-         IF( ktra == jp_sal ) CALL iom_put( 'vdiff_salttr' , rho0     * z2d(:,:) ) ! diffusive salt transport in j-direction
+         IF( ktra == jp_tem )   CALL iom_put( 'vdiff_heattr' , rho0_rcp * z2d(:,:) ) ! diffusive heat transport in j-direction
+         IF( ktra == jp_sal )   CALL iom_put( 'vdiff_salttr' , rho0     * z2d(:,:) ) ! diffusive salt transport in j-direction
       ENDIF
 
    END SUBROUTINE dia_ar5_hst
@@ -380,10 +381,10 @@ CONTAINS
 
          area_tot  = glob_sum( 'diaar5', e1e2t(:,:) )
 
-         ALLOCATE( zvol0(jpi,jpj) )
+         ALLOCATE( zvol0(A2D(0)) )
          zvol0 (:,:) = 0._wp
          thick0(:,:) = 0._wp
-         DO_3D( 1, 1, 1, 1, 1, jpkm1 )   ! interpolation of salinity at the last ocean level (i.e. the partial step)
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   ! interpolation of salinity at the last ocean level (i.e. the partial step)
             zztmp = tmask(ji,jj,jk) * e3t_0(ji,jj,jk)
             zvol0 (ji,jj) = zvol0 (ji,jj) + zztmp * e1e2t(ji,jj)
             thick0(ji,jj) = thick0(ji,jj) + zztmp
@@ -392,16 +393,16 @@ CONTAINS
          DEALLOCATE( zvol0 )
 
          IF( iom_use( 'sshthster' ) ) THEN
-            ALLOCATE( zsaldta(jpi,jpj,jpk,jpts) )
+            ALLOCATE( zsaldta(A2D(0),jpk,jpts) )
             CALL iom_open ( 'sali_ref_clim_monthly', inum )
             CALL iom_get  ( inum, jpdom_global, 'vosaline' , zsaldta(:,:,:,1), 1  )
             CALL iom_get  ( inum, jpdom_global, 'vosaline' , zsaldta(:,:,:,2), 12 )
             CALL iom_close( inum )
 
             sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )
-            sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
+            sn0(:,:,:) = sn0(:,:,:) * tmask(A2D(0),:)
             IF( ln_zps ) THEN               ! z-coord. partial steps
-               DO_2D( 1, 1, 1, 1 )          ! interpolation of salinity at the last ocean level (i.e. the partial step)
+               DO_2D( 0, 0, 0, 0 )          ! interpolation of salinity at the last ocean level (i.e. the partial step)
                   ik = mbkt(ji,jj)
                   IF( ik > 1 ) THEN
                      zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
diff --git a/src/OCE/DIA/diacfl.F90 b/src/OCE/DIA/diacfl.F90
index e35764b3b..7d0ee952d 100644
--- a/src/OCE/DIA/diacfl.F90
+++ b/src/OCE/DIA/diacfl.F90
@@ -51,19 +51,19 @@ CONTAINS
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
       !
-      INTEGER                          ::   ji, jj, jk                       ! dummy loop indices
-      REAL(wp)                         ::   zCu_max, zCv_max, zCw_max        ! local scalars
-      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc  ! workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace
-      LOGICAL , DIMENSION(jpi,jpj,jpk) ::   llmsk
+      INTEGER                         ::   ji, jj, jk                       ! dummy loop indices
+      REAL(wp)                        ::   zCu_max, zCv_max, zCw_max        ! local scalars
+      INTEGER , DIMENSION(3)          ::   iloc_u , iloc_v , iloc_w , iloc  ! workspace
+      REAL(wp), DIMENSION(A2D(0),jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace
+      LOGICAL , DIMENSION(A2D(0),jpk) ::   llmsk
       !!----------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('dia_cfl')
       !
-      llmsk(     1:nn_hls,:,:) = .FALSE.   ! exclude halos from the checked region
-      llmsk(Nie0+1:   jpi,:,:) = .FALSE.
-      llmsk(:,     1:nn_hls,:) = .FALSE.
-      llmsk(:,Nje0+1:   jpj,:) = .FALSE.
+      !llmsk(     1:nn_hls,:,:) = .FALSE.   ! exclude halos from the checked region
+      !llmsk(Nie0+1:   jpi,:,:) = .FALSE.
+      !llmsk(:,     1:nn_hls,:) = .FALSE.
+      !llmsk(:,Nje0+1:   jpj,:) = .FALSE.
       !
       DO_3D( 0, 0, 0, 0, 1, jpk )      ! calculate Courant numbers
          zCu_cfl(ji,jj,jk) = ABS( uu(ji,jj,jk,Kmm) ) * rDt / e1u  (ji,jj)      ! for i-direction
diff --git a/src/OCE/DIA/diahsb.F90 b/src/OCE/DIA/diahsb.F90
index 71c049260..29e80ec7a 100644
--- a/src/OCE/DIA/diahsb.F90
+++ b/src/OCE/DIA/diahsb.F90
@@ -141,12 +141,10 @@ CONTAINS
       !
       IF( ln_linssh ) THEN   ! Advection flux through fixed surface (z=0)
          IF( ln_isfcav ) THEN
-            DO ji=1,jpi
-               DO jj=1,jpj
-                  ztmp(ji,jj,9 ) = - surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb)
-                  ztmp(ji,jj,10) = - surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb)
-               END DO
-            END DO
+            DO_2D( 0, 0, 0, 0 )
+               ztmp(ji,jj,9 ) = - surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb)
+               ztmp(ji,jj,10) = - surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb)
+            END_2D
          ELSE
             DO_2D( 0, 0, 0, 0 )
                ztmp(ji,jj,9 ) = - surf(ji,jj) * ww(ji,jj,1) * ts(ji,jj,1,jp_tem,Kbb)
diff --git a/src/OCE/DIA/diaptr.F90 b/src/OCE/DIA/diaptr.F90
index 6c4f48a3d..68fe7f3a6 100644
--- a/src/OCE/DIA/diaptr.F90
+++ b/src/OCE/DIA/diaptr.F90
@@ -542,7 +542,7 @@ CONTAINS
       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(A2D(0),jpk)     , INTENT(in)  ::  pvflx ! 3D input array of advection/diffusion
       REAL(wp), DIMENSION(A1Dj(0),nbasin)               ::  zsj   !
       INTEGER                                           ::  jn    !
 
@@ -684,8 +684,8 @@ 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(A2D(0),jpk)  ::   pvflx  ! 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(0)) :: p_fval      ! function value
@@ -710,8 +710,8 @@ 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))  ::   pvflx  ! 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) )  ::   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(0)) :: p_fval  ! function value
diff --git a/src/OCE/DIA/diawri.F90 b/src/OCE/DIA/diawri.F90
index b6dd96cda..cf3c014ad 100644
--- a/src/OCE/DIA/diawri.F90
+++ b/src/OCE/DIA/diawri.F90
@@ -122,8 +122,9 @@ CONTAINS
       REAL(wp)::   zztmp , zztmpx   ! local scalar
       REAL(wp)::   zztmp2, zztmpy   !   -      -
       REAL(wp)::   ze3
-      REAL(wp), DIMENSION(A2D(     0))     ::   z2d   ! 2D workspace
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   z3d   ! 3D workspace
+      REAL(wp), DIMENSION(A2D(0))          ::   z2d0   ! 2D workspace
+      REAL(wp), DIMENSION(A2D(0)     ,jpk) ::   z3d0   ! 3D workspace
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   z3d    ! 3D workspace
       !!----------------------------------------------------------------------
       ! 
       IF( ln_timing )   CALL timing_start('dia_wri')
@@ -135,8 +136,9 @@ CONTAINS
       ENDIF
 
       ! initialize arrays
-      z2d(A2D(0))   = 0._wp
-      z3d(A2D(0),:) = 0._wp
+      z2d0(:,:)   = 0._wp
+      z3d0(:,:,:) = 0._wp
+      z3d (:,:,:) = 0._wp
       
       ! Output of initial vertical scale factor
       CALL iom_put("e3t_0", e3t_0(:,:,:) )
@@ -146,44 +148,44 @@ CONTAINS
       !
       IF ( iom_use("tpt_dep") ) THEN
          DO_3D( 0, 0, 0, 0, 1, jpk )
-            z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm)
+            z3d0(ji,jj,jk) = gdept(ji,jj,jk,Kmm)
          END_3D
-         CALL iom_put( "tpt_dep", z3d )
+         CALL iom_put( "tpt_dep", z3d0 )
       ENDIF
 
       ! --- vertical scale factors --- !
       IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN  ! time-varying e3t
-         DO_3D( 0, 0, 0, 0, 1, jpk )
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
             z3d(ji,jj,jk) =  e3t(ji,jj,jk,Kmm)
          END_3D
          CALL iom_put( "e3t", z3d )
          IF ( iom_use("e3tdef") ) THEN
-            DO_3D( 0, 0, 0, 0, 1, jpk )
+            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
                z3d(ji,jj,jk) = ( ( z3d(ji,jj,jk) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
             END_3D
             CALL iom_put( "e3tdef", z3d ) 
          ENDIF
       ENDIF 
       IF ( iom_use("e3u") ) THEN                         ! time-varying e3u
-         DO_3D( 0, 0, 0, 0, 1, jpk )
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
             z3d(ji,jj,jk) =  e3u(ji,jj,jk,Kmm)
          END_3D 
          CALL iom_put( "e3u" , z3d )
       ENDIF
       IF ( iom_use("e3v") ) THEN                         ! time-varying e3v
-         DO_3D( 0, 0, 0, 0, 1, jpk )
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
             z3d(ji,jj,jk) =  e3v(ji,jj,jk,Kmm)
          END_3D
          CALL iom_put( "e3v" , z3d )
       ENDIF
       IF ( iom_use("e3w") ) THEN                         ! time-varying e3w
-         DO_3D( 0, 0, 0, 0, 1, jpk )
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
             z3d(ji,jj,jk) =  e3w(ji,jj,jk,Kmm)
          END_3D
          CALL iom_put( "e3w" , z3d )
       ENDIF
       IF ( iom_use("e3f") ) THEN                         ! time-varying e3f caution here at Kaa
-         DO_3D( 0, 0, 0, 0, 1, jpk )
+         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
             z3d(ji,jj,jk) =  e3f(ji,jj,jk)
          END_3D
          CALL iom_put( "e3f" , z3d )
@@ -213,9 +215,9 @@ CONTAINS
       IF ( iom_use("sbt") ) THEN
          DO_2D( 0, 0, 0, 0 )
             ikbot = mbkt(ji,jj)
-            z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
+            z2d0(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
          END_2D
-         CALL iom_put( "sbt", z2d )                ! bottom temperature
+         CALL iom_put( "sbt", z2d0 )                ! bottom temperature
       ENDIF
       
       CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) )    ! 3D salinity
@@ -223,9 +225,9 @@ CONTAINS
       IF ( iom_use("sbs") ) THEN
          DO_2D( 0, 0, 0, 0 )
             ikbot = mbkt(ji,jj)
-            z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
+            z2d0(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
          END_2D
-         CALL iom_put( "sbs", z2d )                ! bottom salinity
+         CALL iom_put( "sbs", z2d0 )                ! bottom salinity
       ENDIF
 
       IF( .NOT.lk_SWE )   CALL iom_put( "rhop", rhop(:,:,:) )          ! 3D potential density (sigma0)
@@ -233,16 +235,16 @@ CONTAINS
       ! --- momentum --- !
       IF ( iom_use("taubot") ) THEN                ! bottom stress
          zztmp = rho0 * 0.25_wp
-         z2d(:,:) = 0._wp
+         z2d0(:,:) = 0._wp
          DO_2D( 0, 0, 0, 0 )
             zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Kmm)  )**2   &
                &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm)  )**2   &
                &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vv(ji,jj  ,mbkv(ji,jj  ),Kmm)  )**2   &
                &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm)  )**2
-            z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
+            z2d0(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
             !
          END_2D
-         CALL iom_put( "taubot", z2d )           
+         CALL iom_put( "taubot", z2d0 )           
       ENDIF
          
       CALL iom_put( "uoce", uu(:,:,:,Kmm) )            ! 3D i-current
@@ -250,9 +252,9 @@ CONTAINS
       IF ( iom_use("sbu") ) THEN
          DO_2D( 0, 0, 0, 0 )
             ikbot = mbku(ji,jj)
-            z2d(ji,jj) = uu(ji,jj,ikbot,Kmm)
+            z2d0(ji,jj) = uu(ji,jj,ikbot,Kmm)
          END_2D
-         CALL iom_put( "sbu", z2d )                ! bottom i-current
+         CALL iom_put( "sbu", z2d0 )                ! bottom i-current
       ENDIF
       
       CALL iom_put( "voce", vv(:,:,:,Kmm) )            ! 3D j-current
@@ -260,18 +262,15 @@ CONTAINS
       IF ( iom_use("sbv") ) THEN
          DO_2D( 0, 0, 0, 0 )
             ikbot = mbkv(ji,jj)
-            z2d(ji,jj) = vv(ji,jj,ikbot,Kmm)
+            z2d0(ji,jj) = vv(ji,jj,ikbot,Kmm)
          END_2D
-         CALL iom_put( "sbv", z2d )                ! bottom j-current
+         CALL iom_put( "sbv", z2d0 )                ! bottom j-current
       ENDIF
 
       !                                            ! vertical velocity
       IF( ln_zad_Aimp ) THEN
          IF( iom_use('woce') ) THEN
-            DO_3D( 0, 0, 0, 0, 1, jpk )
-               z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
-            END_3D
-            CALL iom_put( "woce", z3d )   ! explicit plus implicit parts
+            CALL iom_put( "woce", ww+wi )   ! explicit plus implicit parts
          ENDIF
       ELSE
          CALL iom_put( "woce", ww )
@@ -281,15 +280,15 @@ CONTAINS
          !                     ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
          IF( ln_zad_Aimp ) THEN
             DO_3D( 0, 0, 0, 0, 1, jpk )
-               z3d(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wi(ji,jj,jk) )
+               z3d0(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wi(ji,jj,jk) )
             END_3D
          ELSE
             DO_3D( 0, 0, 0, 0, 1, jpk )
-               z3d(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ww(ji,jj,jk)
+               z3d0(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ww(ji,jj,jk)
             END_3D
          ENDIF
-         CALL iom_put( "w_masstr" , z3d )  
-         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d * z3d )
+         CALL iom_put( "w_masstr" , z3d0 )  
+         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d0 * z3d0 )
       ENDIF
 
       CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
@@ -304,15 +303,15 @@ CONTAINS
             zztmp  = ts(ji,jj,1,jp_sal,Kmm)
             zztmpx = (ts(ji+1,jj,1,jp_sal,Kmm) - zztmp) * r1_e1u(ji,jj) + (zztmp - ts(ji-1,jj  ,1,jp_sal,Kmm)) * r1_e1u(ji-1,jj)
             zztmpy = (ts(ji,jj+1,1,jp_sal,Kmm) - zztmp) * r1_e2v(ji,jj) + (zztmp - ts(ji  ,jj-1,1,jp_sal,Kmm)) * r1_e2v(ji,jj-1)
-            z2d(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
+            z2d0(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
                &                 * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
          END_2D
-         CALL iom_put( "sssgrad2",  z2d )          ! square of module of sss gradient
+         CALL iom_put( "sssgrad2",  z2d0 )          ! square of module of sss gradient
          IF ( iom_use("sssgrad") ) THEN
             DO_2D( 0, 0, 0, 0 )
-               z2d(ji,jj) = SQRT( z2d(ji,jj) )
+               z2d0(ji,jj) = SQRT( z2d0(ji,jj) )
             END_2D
-            CALL iom_put( "sssgrad",  z2d )        ! module of sss gradient
+            CALL iom_put( "sssgrad",  z2d0 )        ! module of sss gradient
          ENDIF
       ENDIF
          
@@ -321,80 +320,80 @@ CONTAINS
             zztmp  = ts(ji,jj,1,jp_tem,Kmm)
             zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj)
             zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1)
-            z2d(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
+            z2d0(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
                &                 * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
          END_2D
-         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient
+         CALL iom_put( "sstgrad2",  z2d0 )          ! square of module of sst gradient
          IF ( iom_use("sstgrad") ) THEN
             DO_2D( 0, 0, 0, 0 )
-               z2d(ji,jj) = SQRT( z2d(ji,jj) )
+               z2d0(ji,jj) = SQRT( z2d0(ji,jj) )
             END_2D
-            CALL iom_put( "sstgrad",  z2d )        ! module of sst gradient
+            CALL iom_put( "sstgrad",  z2d0 )        ! module of sst gradient
          ENDIF
       ENDIF
          
       ! heat and salt contents
       IF( iom_use("heatc") ) THEN
-         z2d(:,:)  = 0._wp 
+         z2d0(:,:)  = 0._wp 
          DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
+            z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
          END_3D
-         CALL iom_put( "heatc", rho0_rcp * z2d )   ! vertically integrated heat content (J/m2)
+         CALL iom_put( "heatc", rho0_rcp * z2d0 )   ! vertically integrated heat content (J/m2)
       ENDIF
 
       IF( iom_use("saltc") ) THEN
-         z2d(:,:)  = 0._wp 
+         z2d0(:,:)  = 0._wp 
          DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
+            z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
          END_3D
-         CALL iom_put( "saltc", rho0 * z2d )       ! vertically integrated salt content (PSU*kg/m2)
+         CALL iom_put( "saltc", rho0 * z2d0 )       ! vertically integrated salt content (PSU*kg/m2)
       ENDIF
       !
       IF( iom_use("salt2c") ) THEN
-         z2d(:,:)  = 0._wp 
+         z2d0(:,:)  = 0._wp 
          DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
+            z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
          END_3D
-         CALL iom_put( "salt2c", rho0 * z2d )      ! vertically integrated square of salt content (PSU2*kg/m2)
+         CALL iom_put( "salt2c", rho0 * z2d0 )      ! vertically integrated square of salt content (PSU2*kg/m2)
       ENDIF
       !
       IF ( iom_use("ke") .OR. iom_use("ke_int") ) THEN
          DO_3D( 0, 0, 0, 0, 1, jpk )
             zztmpx = uu(ji-1,jj  ,jk,Kmm) + uu(ji,jj,jk,Kmm)
             zztmpy = vv(ji  ,jj-1,jk,Kmm) + vv(ji,jj,jk,Kmm)
-            z3d(ji,jj,jk) = 0.25_wp * ( zztmpx*zztmpx + zztmpy*zztmpy )
+            z3d0(ji,jj,jk) = 0.25_wp * ( zztmpx*zztmpx + zztmpy*zztmpy )
          END_3D
-         CALL iom_put( "ke", z3d )                 ! kinetic energy
+         CALL iom_put( "ke", z3d0 )                 ! kinetic energy
 
-         z2d(:,:)  = 0._wp 
+         z2d0(:,:)  = 0._wp 
          DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
+            z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d0(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
          END_3D
-         CALL iom_put( "ke_int", z2d )             ! vertically integrated kinetic energy
+         CALL iom_put( "ke_int", z2d0 )             ! vertically integrated kinetic energy
       ENDIF
       !
       IF ( iom_use("sKE") ) THEN                   ! surface kinetic energy at T point
-         z2d(:,:) = 0._wp
+         z2d0(:,:) = 0._wp
          DO_2D( 0, 0, 0, 0 )
-            z2d(ji,jj) = 0.25_wp * ( uu(ji  ,jj,1,Kmm) * uu(ji  ,jj,1,Kmm) * e1e2u(ji  ,jj) * e3u(ji  ,jj,1,Kmm)  &
+            z2d0(ji,jj) = 0.25_wp * ( uu(ji  ,jj,1,Kmm) * uu(ji  ,jj,1,Kmm) * e1e2u(ji  ,jj) * e3u(ji  ,jj,1,Kmm)  &
                &                   + uu(ji-1,jj,1,Kmm) * uu(ji-1,jj,1,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,1,Kmm)  &
                &                   + vv(ji,jj  ,1,Kmm) * vv(ji,jj  ,1,Kmm) * e1e2v(ji,jj  ) * e3v(ji,jj  ,1,Kmm)  & 
                &                   + vv(ji,jj-1,1,Kmm) * vv(ji,jj-1,1,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,1,Kmm)  )  &
                &                 * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * ssmask(ji,jj)
          END_2D
-         IF ( iom_use("sKE" ) )  CALL iom_put( "sKE" , z2d )   
+         IF ( iom_use("sKE" ) )  CALL iom_put( "sKE" , z2d0 )   
       ENDIF
       !    
       IF ( iom_use("ssKEf") ) THEN                 ! surface kinetic energy at F point
-         z2d(:,:) = 0._wp                          ! CAUTION : only valid in SWE, not with bathymetry
+         z2d0(:,:) = 0._wp                          ! CAUTION : only valid in SWE, not with bathymetry
          DO_2D( 0, 0, 0, 0 )
-            z2d(ji,jj) = 0.25_wp * ( uu(ji,jj  ,1,Kmm) * uu(ji,jj  ,1,Kmm) * e1e2u(ji,jj  ) * e3u(ji,jj  ,1,Kmm)  &
+            z2d0(ji,jj) = 0.25_wp * ( uu(ji,jj  ,1,Kmm) * uu(ji,jj  ,1,Kmm) * e1e2u(ji,jj  ) * e3u(ji,jj  ,1,Kmm)  &
                &                   + uu(ji,jj+1,1,Kmm) * uu(ji,jj+1,1,Kmm) * e1e2u(ji,jj+1) * e3u(ji,jj+1,1,Kmm)  &
                &                   + vv(ji  ,jj,1,Kmm) * vv(ji,jj  ,1,Kmm) * e1e2v(ji  ,jj) * e3v(ji  ,jj,1,Kmm)  & 
                &                   + vv(ji+1,jj,1,Kmm) * vv(ji+1,jj,1,Kmm) * e1e2v(ji+1,jj) * e3v(ji+1,jj,1,Kmm)  )  &
                &                 * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj)
          END_2D
-         CALL iom_put( "ssKEf", z2d )                     
+         CALL iom_put( "ssKEf", z2d0 )                     
       ENDIF
       !
       CALL iom_put( "hdiv", hdiv )                 ! Horizontal divergence
@@ -402,31 +401,31 @@ CONTAINS
       IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
          
          DO_3D( 0, 0, 0, 0, 1, jpk )
-            z3d(ji,jj,jk) = rho0 * uu(ji,jj,jk,Kmm) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
+            z3d0(ji,jj,jk) = rho0 * uu(ji,jj,jk,Kmm) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
          END_3D
-         CALL iom_put( "u_masstr"     , z3d )      ! mass transport in i-direction
+         CALL iom_put( "u_masstr"     , z3d0 )      ! mass transport in i-direction
          
          IF( iom_use("u_masstr_vint") ) THEN
-            z2d(:,:) = 0._wp 
+            z2d0(:,:) = 0._wp 
             DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-               z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk)
+               z2d0(ji,jj) = z2d0(ji,jj) + z3d0(ji,jj,jk)
             END_3D
-            CALL iom_put( "u_masstr_vint", z2d )   ! mass transport in i-direction vertical sum
+            CALL iom_put( "u_masstr_vint", z2d0 )   ! mass transport in i-direction vertical sum
          ENDIF
          IF( iom_use("u_heattr") ) THEN
-            z2d(:,:) = 0._wp 
+            z2d0(:,:) = 0._wp 
             zztmp = 0.5_wp * rcp
             DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-               z2d(ji,jj) = z2d(ji,jj) + zztmp * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
+               z2d0(ji,jj) = z2d0(ji,jj) + zztmp * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
             END_3D
-            CALL iom_put( "u_heattr", z2d )        ! heat transport in i-direction
+            CALL iom_put( "u_heattr", z2d0 )        ! heat transport in i-direction
          ENDIF
          IF( iom_use("u_salttr") ) THEN
-            z2d(:,:) = 0._wp 
+            z2d0(:,:) = 0._wp 
             DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-               z2d(ji,jj) = z2d(ji,jj) +   0.5 * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
+               z2d0(ji,jj) = z2d0(ji,jj) +   0.5 * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
             END_3D
-            CALL iom_put( "u_salttr", z2d )        ! heat transport in i-direction
+            CALL iom_put( "u_salttr", z2d0 )        ! heat transport in i-direction
          ENDIF
          
       ENDIF
@@ -434,41 +433,41 @@ CONTAINS
       IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
          
          DO_3D( 0, 0, 0, 0, 1, jpk )
-            z3d(ji,jj,jk) = rho0 * vv(ji,jj,jk,Kmm) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
+            z3d0(ji,jj,jk) = rho0 * vv(ji,jj,jk,Kmm) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
          END_3D
-         CALL iom_put( "v_masstr", z3d )           ! mass transport in j-direction
+         CALL iom_put( "v_masstr", z3d0 )           ! mass transport in j-direction
          
          IF( iom_use("v_heattr") ) THEN
-            z2d(:,:) = 0._wp
+            z2d0(:,:) = 0._wp
             zztmp = 0.5_wp * rcp
             DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-               z2d(ji,jj) = z2d(ji,jj) + zztmp * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
+               z2d0(ji,jj) = z2d0(ji,jj) + zztmp * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
             END_3D
-            CALL iom_put( "v_heattr", z2d )        !  heat transport in j-direction
+            CALL iom_put( "v_heattr", z2d0 )        !  heat transport in j-direction
          ENDIF
          IF( iom_use("v_salttr") ) THEN
-            z2d(:,:) = 0._wp 
+            z2d0(:,:) = 0._wp 
             DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-               z2d(ji,jj) = z2d(ji,jj) +   0.5 * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
+               z2d0(ji,jj) = z2d0(ji,jj) +   0.5 * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
             END_3D
-            CALL iom_put( "v_salttr", z2d )        !  heat transport in j-direction
+            CALL iom_put( "v_salttr", z2d0 )        !  heat transport in j-direction
          ENDIF
 
       ENDIF
 
       IF( iom_use("tosmint") ) THEN
-         z2d(:,:) = 0._wp
+         z2d0(:,:) = 0._wp
          DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
+            z2d0(ji,jj) = z2d0(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
          END_3D
-         CALL iom_put( "tosmint", z2d )            ! Vertical integral of temperature
+         CALL iom_put( "tosmint", z2d0 )            ! Vertical integral of temperature
       ENDIF
       IF( iom_use("somint") ) THEN
-         z2d(:,:) = 0._wp
+         z2d0(:,:) = 0._wp
          DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
+            z2d0(ji,jj) = z2d0(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
          END_3D
-         CALL iom_put( "somint", z2d )             ! Vertical integral of salinity
+         CALL iom_put( "somint", z2d0 )             ! Vertical integral of salinity
       ENDIF
 
       CALL iom_put( "bn2", rn2 )                   ! Brunt-Vaisala buoyancy frequency (N^2)
@@ -477,45 +476,47 @@ CONTAINS
       
       ! Output of surface vorticity terms
       !
-      CALL iom_put( "ssplavor", ff_f )             ! planetary vorticity ( f )
-      !
-      IF ( iom_use("ssrelvor")    .OR. iom_use("ssEns")    .OR.   &
+      IF ( iom_use("ssplavor")    .OR. iom_use("ssrelvor")    .OR. iom_use("ssEns")    .OR.   &
          & iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
          !
-         z2d(:,:) = 0._wp 
          DO_2D( 0, 0, 0, 0 )
-            z2d(ji,jj) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm)    &
+            z2d0(ji,jj) = ff_f(ji,jj)
+         END_2D
+         CALL iom_put( "ssplavor", z2d0 )           ! planetary vorticity ( f )
+         
+         DO_2D( 0, 0, 0, 0 )
+            z2d0(ji,jj) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm)    &
             &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,1,Kmm) + e1u(ji,jj) * uu(ji,jj,1,Kmm)  ) * r1_e1e2f(ji,jj)
          END_2D
-         CALL iom_put( "ssrelvor", z2d )           ! relative vorticity ( zeta ) 
+         CALL iom_put( "ssrelvor", z2d0 )           ! relative vorticity ( zeta ) 
          !
          IF ( iom_use("ssEns") .OR. iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
             DO_2D( 0, 0, 0, 0 )  
                ze3 = (  e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1)    &
-                  &    + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  ) + e3t(ji+1,jj  ,1,Kmm) * e1e2t(ji+1,jj  )  ) * r1_e1e2f(ji,jj)
+                  &   + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  ) + e3t(ji+1,jj  ,1,Kmm) * e1e2t(ji+1,jj  )  ) * r1_e1e2f(ji,jj)
                IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3
                ELSE                      ;   ze3 = 0._wp
                ENDIF
-               z2d(ji,jj) = ze3 * z2d(ji,jj) 
+               z2d0(ji,jj) = ze3 * z2d0(ji,jj) 
             END_2D
-            CALL iom_put( "ssrelpotvor", z2d )     ! relative potential vorticity (zeta/h)
+            CALL iom_put( "ssrelpotvor", z2d0 )     ! relative potential vorticity (zeta/h)
             !
             IF ( iom_use("ssEns") .OR. iom_use("ssabspotvor") ) THEN
                DO_2D( 0, 0, 0, 0 )
                   ze3 = (  e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1)    &
-                     &    + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  ) + e3t(ji+1,jj  ,1,Kmm) * e1e2t(ji+1,jj  )  ) * r1_e1e2f(ji,jj)
+                     &   + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  ) + e3t(ji+1,jj  ,1,Kmm) * e1e2t(ji+1,jj  )  ) * r1_e1e2f(ji,jj)
                   IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3
                   ELSE                      ;   ze3 = 0._wp
                   ENDIF
-                  z2d(ji,jj) = ze3 * ff_f(ji,jj) + z2d(ji,jj) 
+                  z2d0(ji,jj) = ze3 * ff_f(ji,jj) + z2d0(ji,jj) 
                END_2D
-               CALL iom_put( "ssabspotvor", z2d )  ! absolute potential vorticity ( q )
+               CALL iom_put( "ssabspotvor", z2d0 )  ! absolute potential vorticity ( q )
                !
                IF ( iom_use("ssEns") ) THEN
                   DO_2D( 0, 0, 0, 0 )  
-                     z2d(ji,jj) = 0.5_wp * z2d(ji,jj) * z2d(ji,jj) 
+                     z2d0(ji,jj) = 0.5_wp * z2d0(ji,jj) * z2d0(ji,jj) 
                   END_2D
-                  CALL iom_put( "ssEns", z2d )     ! potential enstrophy ( 1/2*q2 )
+                  CALL iom_put( "ssEns", z2d0 )     ! potential enstrophy ( 1/2*q2 )
                ENDIF
             ENDIF
          ENDIF
@@ -582,8 +583,8 @@ CONTAINS
       INTEGER  ::   jn, ierror                               ! local integers
       REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
       !
-      REAL(wp), DIMENSION(jpi,jpj    ) :: z2d     ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d     ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj    ) :: z2d0     ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d0     ! 3D workspace
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! ABL 3D workspace
       !!----------------------------------------------------------------------
       !
@@ -935,21 +936,21 @@ CONTAINS
 
       IF( .NOT.ln_linssh ) THEN
          DO_3D( 0, 0, 0, 0, 1, jpk )
-            z3d(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm)
+            z3d0(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm)
          END_3D
-         CALL histwrite( nid_T, "votemper", it, z3d, ndim_T , ndex_T  )   ! heat content
+         CALL histwrite( nid_T, "votemper", it, z3d0, ndim_T , ndex_T  )   ! heat content
          DO_3D( 0, 0, 0, 0, 1, jpk )
-            z3d(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm)
+            z3d0(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm)
          END_3D
-         CALL histwrite( nid_T, "vosaline", it, z3d, ndim_T , ndex_T  )   ! salt content
+         CALL histwrite( nid_T, "vosaline", it, z3d0, ndim_T , ndex_T  )   ! salt content
          DO_2D( 0, 0, 0, 0 )
-            z2d(ji,jj   ) = ts(ji,jj, 1,jp_tem,Kmm) * e3t(ji,jj, 1,Kmm)
+            z2d0(ji,jj   ) = ts(ji,jj, 1,jp_tem,Kmm) * e3t(ji,jj, 1,Kmm)
          END_2D
-         CALL histwrite( nid_T, "sosstsst", it, z2d, ndim_hT, ndex_hT )   ! sea surface heat content
+         CALL histwrite( nid_T, "sosstsst", it, z2d0, ndim_hT, ndex_hT )   ! sea surface heat content
          DO_2D( 0, 0, 0, 0 )
-            z2d(ji,jj   ) = ts(ji,jj, 1,jp_sal,Kmm) * e3t(ji,jj, 1,Kmm)
+            z2d0(ji,jj   ) = ts(ji,jj, 1,jp_sal,Kmm) * e3t(ji,jj, 1,Kmm)
          END_2D
-         CALL histwrite( nid_T, "sosaline", it, z2d, ndim_hT, ndex_hT )   ! sea surface salinity content
+         CALL histwrite( nid_T, "sosaline", it, z2d0, ndim_hT, ndex_hT )   ! sea surface salinity content
       ELSE
          CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature
          CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T  )   ! salinity
@@ -958,41 +959,41 @@ CONTAINS
       ENDIF
       IF( .NOT.ln_linssh ) THEN
          DO_3D( 0, 0, 0, 0, 1, jpk )
-           z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm)     ! 3D workspace for qco substitution
+           z3d0(ji,jj,jk) = e3t(ji,jj,jk,Kmm)     ! 3D workspace for qco substitution
          END_3D
-         CALL histwrite( nid_T, "vovvle3t", it, z3d        , ndim_T , ndex_T  )   ! level thickness
+         CALL histwrite( nid_T, "vovvle3t", it, z3d0        , ndim_T , ndex_T  )   ! level thickness
          DO_3D( 0, 0, 0, 0, 1, jpk )
-           z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm)   ! 3D workspace for qco substitution
+           z3d0(ji,jj,jk) = gdept(ji,jj,jk,Kmm)   ! 3D workspace for qco substitution
          END_3D
-         CALL histwrite( nid_T, "vovvldep", it, z3d        , ndim_T , ndex_T  )   ! t-point depth 
+         CALL histwrite( nid_T, "vovvldep", it, z3d0        , ndim_T , ndex_T  )   ! t-point depth 
          DO_3D( 0, 0, 0, 0, 1, jpk )
-            z3d(ji,jj,jk) = ( ( e3t(ji,jj,jk,Kmm) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
+            z3d0(ji,jj,jk) = ( ( e3t(ji,jj,jk,Kmm) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
          END_3D         
-         CALL histwrite( nid_T, "vovvldef", it, z3d        , ndim_T , ndex_T  )   ! level thickness deformation
+         CALL histwrite( nid_T, "vovvldef", it, z3d0        , ndim_T , ndex_T  )   ! level thickness deformation
       ENDIF
       CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm)  , ndim_hT, ndex_hT )   ! sea surface height
       DO_2D( 0, 0, 0, 0 )
-         z2d(ji,jj) = emp(ji,jj) - rnf(ji,jj)
+         z2d0(ji,jj) = emp(ji,jj) - rnf(ji,jj)
       END_2D
-      CALL histwrite( nid_T, "sowaflup", it, z2d           , ndim_hT, ndex_hT )   ! upward water flux
+      CALL histwrite( nid_T, "sowaflup", it, z2d0           , ndim_hT, ndex_hT )   ! upward water flux
       CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
       CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux 
                                                                                   ! (includes virtual salt flux beneath ice 
                                                                                   ! in linear free surface case)
       IF( ln_linssh ) THEN
          DO_2D( 0, 0, 0, 0 )
-            z2d(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_tem,Kmm)
+            z2d0(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_tem,Kmm)
          END_2D
-         CALL histwrite( nid_T, "sosst_cd", it, z2d, ndim_hT, ndex_hT )          ! c/d term on sst
+         CALL histwrite( nid_T, "sosst_cd", it, z2d0, ndim_hT, ndex_hT )          ! c/d term on sst
          DO_2D( 0, 0, 0, 0 )
-            z2d(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_sal,Kmm)
+            z2d0(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_sal,Kmm)
          END_2D
-         CALL histwrite( nid_T, "sosss_cd", it, z2d, ndim_hT, ndex_hT )          ! c/d term on sss
+         CALL histwrite( nid_T, "sosss_cd", it, z2d0, ndim_hT, ndex_hT )          ! c/d term on sss
       ENDIF
       DO_2D( 0, 0, 0, 0 )
-         z2d(ji,jj) = qsr(ji,jj) + qns(ji,jj)
+         z2d0(ji,jj) = qsr(ji,jj) + qns(ji,jj)
       END_2D
-      CALL histwrite( nid_T, "sohefldo", it, z2d           , ndim_hT, ndex_hT )   ! total heat flux
+      CALL histwrite( nid_T, "sohefldo", it, z2d0           , ndim_hT, ndex_hT )   ! total heat flux
       CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
       IF( ALLOCATED(hmld) ) THEN   ! zdf_mxl not called by SWE
          CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
@@ -1051,9 +1052,9 @@ CONTAINS
          CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
          CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
          DO_2D( 0, 0, 0, 0 )
-            z2d(ji,jj) = erp(ji,jj) * ts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1)
+            z2d0(ji,jj) = erp(ji,jj) * ts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1)
          END_2D
-         CALL histwrite( nid_T, "sosafldp", it, z2d           , ndim_hT, ndex_hT )   ! salt flux damping
+         CALL histwrite( nid_T, "sosafldp", it, z2d0           , ndim_hT, ndex_hT )   ! salt flux damping
       ENDIF
 !      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
 !      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
@@ -1073,9 +1074,9 @@ CONTAINS
 
       IF( ln_zad_Aimp ) THEN
          DO_3D( 0, 0, 0, 0, 1, jpk )
-           z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
+           z3d0(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
          END_3D         
-         CALL histwrite( nid_W, "vovecrtz", it, z3d         , ndim_T, ndex_T )    ! vert. current
+         CALL histwrite( nid_W, "vovecrtz", it, z3d0         , ndim_T, ndex_T )    ! vert. current
       ELSE
          CALL histwrite( nid_W, "vovecrtz", it, ww          , ndim_T, ndex_T )    ! vert. current
       ENDIF
@@ -1124,8 +1125,8 @@ CONTAINS
       !!
       INTEGER ::   ji, jj, jk       ! dummy loop indices
       INTEGER ::   inum
-      REAL(wp), DIMENSION(jpi,jpj)     :: z2d      
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d      
+      REAL(wp), DIMENSION(A2D(0))     :: z2d0      
+      REAL(wp), DIMENSION(A2D(0),jpk) :: z3d0      
       !!----------------------------------------------------------------------
       ! 
       IF(lwp) THEN
@@ -1144,9 +1145,9 @@ CONTAINS
       CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)        )    ! now j-velocity
       IF( ln_zad_Aimp ) THEN
          DO_3D( 0, 0, 0, 0, 1, jpk )
-           z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
+           z3d0(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
          END_3D
-         CALL iom_rstput( 0, 0, inum, 'vovecrtz', z3d            )    ! now k-velocity
+         CALL iom_rstput( 0, 0, inum, 'vovecrtz', z3d0           )    ! now k-velocity
       ELSE
          CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity
       ENDIF
@@ -1182,26 +1183,26 @@ CONTAINS
          CALL iom_rstput( 0, 0, inum,  'ahmf', ahmf              )    ! ahmf at v-point
       ENDIF
       DO_2D( 0, 0, 0, 0 )
-         z2d(ji,jj) = emp(ji,jj) - rnf(ji,jj)
+         z2d0(ji,jj) = emp(ji,jj) - rnf(ji,jj)
       END_2D
-      CALL iom_rstput( 0, 0, inum, 'sowaflup', z2d               )    ! freshwater budget
+      CALL iom_rstput( 0, 0, inum, 'sowaflup', z2d0              )    ! freshwater budget
       DO_2D( 0, 0, 0, 0 )
-         z2d(ji,jj) = qsr(ji,jj) + qns(ji,jj)
+         z2d0(ji,jj) = qsr(ji,jj) + qns(ji,jj)
       END_2D
-      CALL iom_rstput( 0, 0, inum, 'sohefldo', z2d               )    ! total heat flux
+      CALL iom_rstput( 0, 0, inum, 'sohefldo', z2d0              )    ! total heat flux
       CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
       CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
       CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
       CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
       IF(  .NOT.ln_linssh  ) THEN
          DO_3D( 0, 0, 0, 0, 1, jpk )
-           z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm)   ! 3D workspace for qco substitution
+           z3d0(ji,jj,jk) = gdept(ji,jj,jk,Kmm)   ! 3D workspace for qco substitution
          END_3D
-         CALL iom_rstput( 0, 0, inum, 'vovvldep', z3d            )    !  T-cell depth 
+         CALL iom_rstput( 0, 0, inum, 'vovvldep', z3d0           )    !  T-cell depth 
          DO_3D( 0, 0, 0, 0, 1, jpk )
-           z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm)     ! 3D workspace for qco substitution
+           z3d0(ji,jj,jk) = e3t(ji,jj,jk,Kmm)     ! 3D workspace for qco substitution
          END_3D
-         CALL iom_rstput( 0, 0, inum, 'vovvle3t', z3d            )    !  T-cell thickness  
+         CALL iom_rstput( 0, 0, inum, 'vovvle3t', z3d0           )    !  T-cell thickness  
       END IF
       IF( ln_wave .AND. ln_sdw ) THEN
          CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd            )    ! now StokesDrift i-velocity
diff --git a/src/OCE/DYN/sshwzv.F90 b/src/OCE/DYN/sshwzv.F90
index 4d2ca527f..aed60c986 100644
--- a/src/OCE/DYN/sshwzv.F90
+++ b/src/OCE/DYN/sshwzv.F90
@@ -175,7 +175,8 @@ CONTAINS
          IF(lwp) WRITE(numout,*) 'wzv_MLF : now vertical velocity '
          IF(lwp) WRITE(numout,*) '~~~~~~~'
          !
-         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
+         pww(:,:,:) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
+         !                                   ! needed over the halos for the output (ww+wi) in diawri.F90
       ENDIF
       !                                           !------------------------------!
       !                                           !     Now Vertical Velocity    !
diff --git a/src/OCE/LDF/ldftra.F90 b/src/OCE/LDF/ldftra.F90
index 87e1d7465..5536f96ff 100644
--- a/src/OCE/LDF/ldftra.F90
+++ b/src/OCE/LDF/ldftra.F90
@@ -794,8 +794,8 @@ CONTAINS
       !
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
       REAL(wp) ::   zztmp   ! local scalar
-      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zw2d   ! 2D workspace
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zw3d   ! 3D workspace
+      REAL(wp), DIMENSION(A2D(0))     ::   zw2d   ! 2D workspace
+      REAL(wp), DIMENSION(A2D(0),jpk) ::   zw3d   ! 3D workspace
       !!----------------------------------------------------------------------
       !
 !!gm I don't like this routine....   Crazy  way of doing things, not optimal at all...
@@ -806,42 +806,53 @@ CONTAINS
       !
       !                                                  !==  eiv velocities: calculate and output  ==!
       !
-      zw3d(:,:,jpk) = 0._wp                                    ! bottom value always 0
+      zw3d(:,:,jpk) = 0._wp                                   ! bottom value always 0
       !
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                                  ! e2u e3u u_eiv = -dk[psi_uw]
-         zw3d(ji,jj,jk) = ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm) )
-      END_3D
-      CALL iom_put( "uoce_eiv", zw3d )
+      IF( iom_use('uoce_eiv') ) THEN
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )                                  ! e2u e3u u_eiv = -dk[psi_uw]
+            zw3d(ji,jj,jk) = ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm) )
+         END_3D
+         CALL iom_put( "uoce_eiv", zw3d )
+      ENDIF
       !
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                                  ! e1v e3v v_eiv = -dk[psi_vw]
-         zw3d(ji,jj,jk) = ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm) )
-      END_3D
-      CALL iom_put( "voce_eiv", zw3d )
+      IF( iom_use('ueiv_masstr') ) THEN
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zw3d(ji,jj,jk) = rho0 * ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) )
+         END_3D
+         CALL iom_put( "ueiv_masstr", zw3d )                  ! mass transport in i-direction
+      ENDIF
       !
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                            ! e1 e2 w_eiv = dk[psix] + dk[psix]
-         zw3d(ji,jj,jk) = (  psi_vw(ji,jj,jk) - psi_vw(ji  ,jj-1,jk)    &
-            &              + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj  ,jk)  ) / e1e2t(ji,jj)
-      END_3D
-      CALL iom_put( "woce_eiv", zw3d )
+      IF( iom_use('uoce_eiv') ) THEN
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )                                  ! e1v e3v v_eiv = -dk[psi_vw]
+            zw3d(ji,jj,jk) = ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm) )
+         END_3D
+         CALL iom_put( "voce_eiv", zw3d )
+      ENDIF
       !
-      IF( iom_use('weiv_masstr') ) THEN   ! vertical mass transport & its square value
-         DO_2D( 0, 0, 0, 0 )
-            zw2d(ji,jj) = rho0 * e1e2t(ji,jj)
-         END_2D
-         DO jk = 1, jpk
-            zw3d(:,:,jk) = zw3d(:,:,jk) * zw2d(:,:)
-         END DO
-         CALL iom_put( "weiv_masstr" , zw3d )
+      IF( iom_use('veiv_masstr') ) THEN
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zw3d(ji,jj,jk) = rho0 * ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) )
+         END_3D
+         CALL iom_put( "veiv_masstr", zw3d )                  ! mass transport in j-direction
+      ENDIF
+       !
+      IF( iom_use('woce_eiv') ) THEN
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )                                  ! e1 e2 w_eiv = dk[psix] + dk[psix]
+            zw3d(ji,jj,jk) = (  psi_vw(ji,jj,jk) - psi_vw(ji  ,jj-1,jk)    &
+               &              + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj  ,jk)  ) / e1e2t(ji,jj)
+         END_3D
+         CALL iom_put( "woce_eiv", zw3d )
       ENDIF
       !
-      IF( iom_use('ueiv_masstr') ) THEN
-         zw3d(:,:,:) = 0.e0
-         DO jk = 1, jpkm1
-            zw3d(:,:,jk) = rho0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) )
-         END DO
-         CALL iom_put( "ueiv_masstr", zw3d )                  ! mass transport in i-direction
+      IF( iom_use('weiv_masstr') ) THEN
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zw3d(ji,jj,jk) = rho0 * (  psi_vw(ji,jj,jk) - psi_vw(ji  ,jj-1,jk)  &
+               &                     + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj  ,jk)  )
+         END_3D
+         CALL iom_put( "weiv_masstr" , zw3d )                 ! mass transport in z-direction
       ENDIF
       !
+      !
       zztmp = 0.5_wp * rho0 * rcp
       IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN
         zw2d(:,:)   = 0._wp
@@ -855,49 +866,47 @@ CONTAINS
         CALL iom_put( "ueiv_heattr3d", zztmp * zw3d )                  ! heat transport in i-direction
       ENDIF
       !
-      IF( iom_use('veiv_masstr') ) THEN
-         zw3d(:,:,:) = 0.e0
-         DO jk = 1, jpkm1
-            zw3d(:,:,jk) = rho0 * ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) )
-         END DO
-         CALL iom_put( "veiv_masstr", zw3d )                  ! mass transport in i-direction
+      IF( iom_use('veiv_heattr') .OR. iom_use('veiv_heattr3d') .OR. iom_use('sophteiv') ) THEN
+         zw2d(:,:)   = 0._wp
+         zw3d(:,:,:) = 0._wp
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
+               &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji,jj+1,jk,jp_tem,Kmm) )
+            zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
+         END_3D
+         CALL iom_put( "veiv_heattr"  , zztmp * zw2d )                  !  heat transport in j-direction
+         CALL iom_put( "veiv_heattr3d", zztmp * zw3d )                  !  heat transport in j-direction
+         !
+         IF( iom_use( 'sophteiv' ) .AND. l_diaptr )   CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )
       ENDIF
       !
-      zw2d(:,:)   = 0._wp
-      zw3d(:,:,:) = 0._wp
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-         zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
-            &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji,jj+1,jk,jp_tem,Kmm) )
-         zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
-      END_3D
-      CALL iom_put( "veiv_heattr"  , zztmp * zw2d )                  !  heat transport in j-direction
-      CALL iom_put( "veiv_heattr3d", zztmp * zw3d )                  !  heat transport in j-direction
-      !
-      IF( iom_use( 'sophteiv' ) .AND. l_diaptr )   CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )
       !
       zztmp = 0.5_wp * 0.5
-      IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN
-        zw2d(:,:) = 0._wp
-        zw3d(:,:,:) = 0._wp
-        DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-           zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   &
-              &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji+1,jj,jk,jp_sal,Kmm) )
-           zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
-        END_3D
-        CALL iom_put( "ueiv_salttr", zztmp * zw2d )                  ! salt transport in i-direction
-        CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                ! salt transport in i-direction
+      IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d') ) THEN
+         zw2d(:,:) = 0._wp
+         zw3d(:,:,:) = 0._wp
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   &
+               &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji+1,jj,jk,jp_sal,Kmm) )
+            zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
+         END_3D
+         CALL iom_put( "ueiv_salttr", zztmp * zw2d )                    ! salt transport in i-direction
+         CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                  ! salt transport in i-direction
       ENDIF
-      zw2d(:,:) = 0._wp
-      zw3d(:,:,:) = 0._wp
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-         zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
-            &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji,jj+1,jk,jp_sal,Kmm) )
-         zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
-      END_3D
-      CALL iom_put( "veiv_salttr"  , zztmp * zw2d )                  !  salt transport in j-direction
-      CALL iom_put( "veiv_salttr3d", zztmp * zw3d )                  !  salt transport in j-direction
       !
-      IF( iom_use( 'sopsteiv' ) .AND. l_diaptr ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )
+      IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d') .OR. iom_use('sopsteiv') ) THEN
+         zw2d(:,:) = 0._wp
+         zw3d(:,:,:) = 0._wp
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   &
+               &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji,jj+1,jk,jp_sal,Kmm) )
+            zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk)
+         END_3D
+         CALL iom_put( "veiv_salttr"  , zztmp * zw2d )                  !  salt transport in j-direction
+         CALL iom_put( "veiv_salttr3d", zztmp * zw3d )                  !  salt transport in j-direction
+         !
+         IF( iom_use( 'sopsteiv' ) .AND. l_diaptr )   CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )
+      ENDIF
       !
       !
    END SUBROUTINE ldf_eiv_dia
diff --git a/src/OCE/SBC/sbccpl.F90 b/src/OCE/SBC/sbccpl.F90
index dcb56415d..f681b98be 100644
--- a/src/OCE/SBC/sbccpl.F90
+++ b/src/OCE/SBC/sbccpl.F90
@@ -1820,7 +1820,7 @@ CONTAINS
          &                                                         - zevap_ice_total(:,:) * picefr(:,:) ) * smask0(:,:)  )  ! ice-free oce evap (cell average)
       ! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf
 !!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff
-!!      IF( srcv(jpr_isf)%laction )   CALL iom_put( 'iceshelf_cea', fwfisf(:,:) * tmask(:,:,1)                          )  ! iceshelf
+      IF( srcv(jpr_isf)%laction )    CALL iom_put( 'iceshelf_cea', frcv(jpr_isf)%z3(:,:,1) * smask0(:,:)                 )  ! iceshelf
       !
       !                                                      ! ========================= !
       SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !  ice topmelt and botmelt  !
diff --git a/src/OCE/TRA/traadv_cen.F90 b/src/OCE/TRA/traadv_cen.F90
index 1666eee41..72761ac41 100644
--- a/src/OCE/TRA/traadv_cen.F90
+++ b/src/OCE/TRA/traadv_cen.F90
@@ -511,7 +511,7 @@ CONTAINS
             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
          ENDIF
          !                                 ! "Poleward" heat and salt transports
-         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(A2D(0),:) )
          !                                 !  heat and salt transport
          IF( l_hst )   CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
          !
diff --git a/src/OCE/TRA/traadv_cen_lf.F90 b/src/OCE/TRA/traadv_cen_lf.F90
index 6d1f08fb9..0c8e2bf76 100644
--- a/src/OCE/TRA/traadv_cen_lf.F90
+++ b/src/OCE/TRA/traadv_cen_lf.F90
@@ -175,7 +175,7 @@ CONTAINS
             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
          ENDIF
          !                                 ! "Poleward" heat and salt transports
-         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(A2D(0),:) )
          !                                 !  heat and salt transport
          IF( l_hst )   CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
          !
diff --git a/src/OCE/TRA/traadv_fct.F90 b/src/OCE/TRA/traadv_fct.F90
index e66936ade..7ab4ae15e 100644
--- a/src/OCE/TRA/traadv_fct.F90
+++ b/src/OCE/TRA/traadv_fct.F90
@@ -130,7 +130,7 @@ CONTAINS
       ENDIF
       !
       IF( l_ptr ) THEN
-         ALLOCATE( zptry(A2D(nn_hls),jpk) )
+         ALLOCATE( zptry(A2D(0),jpk) )
          zptry(:,:,:) = 0._wp
       ENDIF
       !
@@ -213,7 +213,7 @@ CONTAINS
             ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:)
          END IF
          !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( l_ptr )   zptry(:,:,:) = zwy(:,:,:)
+         IF( l_ptr )   zptry(:,:,:) = zwy(A2D(0),:)
          !
          !        !==  anti-diffusive flux : high order minus low order  ==!
          !
@@ -364,7 +364,7 @@ CONTAINS
             !
          ENDIF
          IF( l_ptr ) THEN              ! "Poleward" transports
-            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< add anti-diffusive fluxes
+            zptry(:,:,:) = zptry(:,:,:) + zwy(A2D(0),:)  ! <<< add anti-diffusive fluxes
             CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
          ENDIF
          !
@@ -753,7 +753,7 @@ CONTAINS
       ENDIF
       !
       IF( l_ptr ) THEN
-         ALLOCATE( zptry(jpi,jpj,jpk) )
+         ALLOCATE( zptry(A2D(0),jpk) )
          zptry(:,:,:) = 0._wp
       ENDIF
       !
@@ -838,7 +838,7 @@ CONTAINS
             ztrdx(:,:,:) = zwx_3d(:,:,:)   ;   ztrdy(:,:,:) = zwy_3d(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:)
          END IF
          !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( l_ptr )   zptry(:,:,:) = zwy_3d(:,:,:)
+         IF( l_ptr )   zptry(:,:,:) = zwy_3d(A2D(0),:)
          !
          !        !==  anti-diffusive flux : high order minus low order  ==!
          !
@@ -986,7 +986,7 @@ CONTAINS
          ENDIF
          ! NOT TESTED - NEED l_ptr TRUE
          IF( l_ptr ) THEN              ! "Poleward" transports
-            zptry(:,:,:) = zptry(:,:,:) + zwy_3d(:,:,:)  ! <<< add anti-diffusive fluxes
+            zptry(:,:,:) = zptry(:,:,:) + zwy_3d(A2D(0),:)  ! <<< add anti-diffusive fluxes
             CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
          ENDIF
          !
diff --git a/src/OCE/TRA/traadv_mus.F90 b/src/OCE/TRA/traadv_mus.F90
index c22fcfd01..314f4c992 100644
--- a/src/OCE/TRA/traadv_mus.F90
+++ b/src/OCE/TRA/traadv_mus.F90
@@ -401,7 +401,7 @@ CONTAINS
             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) )
          END IF
          !                                 ! "Poleward" heat and salt transports
-         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(A2D(0),:) )
          !                                 !  heat transport
          IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
          !
diff --git a/src/OCE/TRA/traadv_qck.F90 b/src/OCE/TRA/traadv_qck.F90
index cdb96902b..f3976bbd0 100644
--- a/src/OCE/TRA/traadv_qck.F90
+++ b/src/OCE/TRA/traadv_qck.F90
@@ -301,7 +301,7 @@ CONTAINS
          !                                 ! trend diagnostics
          IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(A2D(0),:) )
          !
       END DO
       !
diff --git a/src/OCE/TRA/traadv_qck_lf.F90 b/src/OCE/TRA/traadv_qck_lf.F90
index e866bd809..e3746d601 100644
--- a/src/OCE/TRA/traadv_qck_lf.F90
+++ b/src/OCE/TRA/traadv_qck_lf.F90
@@ -268,7 +268,7 @@ CONTAINS
          !                                 ! trend diagnostics
          IF( l_trd )   CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(A2D(0),:) )
          !
       END DO
       !
diff --git a/src/OCE/TRA/traadv_ubs.F90 b/src/OCE/TRA/traadv_ubs.F90
index 6d65af219..94d783330 100644
--- a/src/OCE/TRA/traadv_ubs.F90
+++ b/src/OCE/TRA/traadv_ubs.F90
@@ -184,7 +184,7 @@ CONTAINS
          END IF
          !
          !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) )
+         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(A2D(0),:) )
          !                                !  heati/salt transport
          IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) )
          !
diff --git a/src/OCE/TRA/traadv_ubs_lf.F90 b/src/OCE/TRA/traadv_ubs_lf.F90
index ee8acfe29..07b3d5e34 100644
--- a/src/OCE/TRA/traadv_ubs_lf.F90
+++ b/src/OCE/TRA/traadv_ubs_lf.F90
@@ -190,7 +190,7 @@ CONTAINS
          END IF
          !
          !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) )
+         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(A2D(0),:) )
          !                                !  heati/salt transport
          IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) )
          !
diff --git a/src/OCE/TRA/tradmp.F90 b/src/OCE/TRA/tradmp.F90
index 16cab127e..164203528 100644
--- a/src/OCE/TRA/tradmp.F90
+++ b/src/OCE/TRA/tradmp.F90
@@ -94,7 +94,7 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts             ! active tracers and RHS of tracer equation
       !
       INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts)     ::  zts_dta
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts) ::  zts_dta
       REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE ::  zwrk
       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  ztrdts
       !!----------------------------------------------------------------------
@@ -146,7 +146,7 @@ CONTAINS
       !
       ! outputs (clem trunk)
       IF( iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN
-         ALLOCATE( zwrk(A2D(nn_hls),jpk) )          ! Needed to handle expressions containing e3t when using key_qco or key_linssh
+         ALLOCATE( zwrk(A2D(0),jpk) )          ! Needed to handle expressions containing e3t when using key_qco or key_linssh
          zwrk(:,:,:) = 0._wp
 
          IF( iom_use('hflx_dmp_cea') ) THEN
diff --git a/src/OCE/TRA/traldf_iso.F90 b/src/OCE/TRA/traldf_iso.F90
index 8d84e7689..d814e0584 100644
--- a/src/OCE/TRA/traldf_iso.F90
+++ b/src/OCE/TRA/traldf_iso.F90
@@ -388,7 +388,7 @@ CONTAINS
             !
             !                             ! "Poleward" diffusive heat or salt transports (T-S case only)
                ! note sign is reversed to give down-gradient diffusive transports )
-            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -zftv(:,:,:)  )
+            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -zftv(A2D(0),:)  )
             !                          ! Diffusive heat transports
             IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', -zftu(:,:,:), -zftv(:,:,:) )
             !
diff --git a/src/OCE/TRA/traldf_lap_blp.F90 b/src/OCE/TRA/traldf_lap_blp.F90
index 16e5df16c..1bb51569d 100644
--- a/src/OCE/TRA/traldf_lap_blp.F90
+++ b/src/OCE/TRA/traldf_lap_blp.F90
@@ -171,7 +171,7 @@ CONTAINS
          IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR.  &     !==  first pass only (  laplacian)  ==!
              ( kpass == 2 .AND.      ln_traldf_blp ) ) THEN      !==  2nd   pass only (bilaplacian)  ==!
 
-            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -ztv(:,:,:)  )
+            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -ztv(A2D(0),:)  )
             IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', -ztu(:,:,:), -ztv(:,:,:) )
          ENDIF
          !                          ! ==================
diff --git a/src/OCE/TRA/traldf_triad.F90 b/src/OCE/TRA/traldf_triad.F90
index 19039b88c..0b341db58 100644
--- a/src/OCE/TRA/traldf_triad.F90
+++ b/src/OCE/TRA/traldf_triad.F90
@@ -485,7 +485,7 @@ CONTAINS
              ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==!
             !
             !                          ! "Poleward" diffusive heat or salt transports (T-S case only)
-            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:)  )
+            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', zftv(A2D(0),:)  )
             !                          ! Diffusive heat transports
             IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:), zftv(:,:,:) )
             !
diff --git a/src/OCE/ZDF/zdfphy.F90 b/src/OCE/ZDF/zdfphy.F90
index 746614a3b..ed3c93b47 100644
--- a/src/OCE/ZDF/zdfphy.F90
+++ b/src/OCE/ZDF/zdfphy.F90
@@ -387,7 +387,7 @@ CONTAINS
       IF( iom_use('avt_k') .OR. iom_use('avm_k') .OR. iom_use('eshear_k') .OR. iom_use('estrat_k') ) THEN
          IF( l_zdfsh2 ) THEN
             CALL iom_put( 'avt_k'   ,   avt_k       * wmask(A2D(0),:) )
-            CALL iom_put( 'avm_k'   ,   avm_k       * wmask(A2D(0),:) )
+            CALL iom_put( 'avm_k'   ,   avm_k       * wmask(:,:,:)    )
             CALL iom_put( 'eshear_k',   zsh2        * wmask(A2D(0),:) )
             CALL iom_put( 'estrat_k', - avt_k * rn2 * wmask(A2D(0),:) )
          ENDIF
-- 
GitLab