Index: NEMO/branches/2021/ticket2680_C1D_PAPA/cfgs/SHARED/domain_def_nemo.xml
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/cfgs/SHARED/domain_def_nemo.xml	(revision 15020)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/cfgs/SHARED/domain_def_nemo.xml	(revision 15022)
@@ -4,4 +4,5 @@
   <!--   T grid   -->
   <domain id="grid_T" long_name="grid T"/>
+  <domain id="grid_T_inner" long_name="grid T inner"/>
 
   <!--   My zoom: example of hand defined zoom   -->
@@ -167,4 +168,5 @@
   <!--   U grid   -->
   <domain id="grid_U" long_name="grid U"/>
+  <domain id="grid_U_inner" long_name="grid U inner"/>
   <!--   Eq section   -->
   <domain id="EqU" domain_ref="grid_U" > <zoom_domain id="EqU"/> </domain>
@@ -173,4 +175,5 @@
   <!--   V grid   -->
   <domain id="grid_V" long_name="grid V"/>
+  <domain id="grid_V_inner" long_name="grid V inner"/>
   <!--   Eq section : no V point on the Equator...   -->
 
@@ -178,4 +181,5 @@
   <!--   W grid   -->
   <domain id="grid_W" long_name="grid W"/>
+  <domain id="grid_W_inner" long_name="grid W inner"/>
   <!--   Eq section   -->
   <domain id="EqW" domain_ref="grid_W" > <zoom_domain id="EqW"/> </domain>
@@ -183,4 +187,5 @@
   <!--   F grid   -->
   <domain id="grid_F" long_name="grid F"/>
+  <domain id="grid_F_inner" long_name="grid F inner"/>
 
   <!--   zonal mean grid   -->
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/cfgs/SHARED/field_def_nemo-oce.xml	(revision 15020)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/cfgs/SHARED/field_def_nemo-oce.xml	(revision 15022)
@@ -131,11 +131,8 @@
 
     <!-- additions to diawri.F90 -->
-    <field id="socegrad"    long_name="module of salinity gradient"              unit="psu/m"   grid_ref="grid_T_3D"/>
-    <field id="socegrad2"   long_name="square of module of salinity gradient"    unit="psu2/m2" grid_ref="grid_T_3D"/>
+    <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"   />
-    <field id="relvor"      long_name="relative vorticity"                       unit="s-1"    grid_ref="grid_T_3D"/>
-    <field id="absvor"      long_name="absolute vorticity"                       unit="s-1"    grid_ref="grid_T_3D"/>
-    <field id="potvor"      long_name="potential vorticity"                      unit="s-1"    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="salt2c"      long_name="Salt content vertically integrated"       unit="1e-3*kg/m2" />
 
@@ -149,8 +146,8 @@
     <field id="sstmax"       long_name="max of sea surface temperature"   field_ref="sst"   operation="maximum"                                                 />
     <field id="sstmin"       long_name="min of sea surface temperature"   field_ref="sst"   operation="minimum"                                                 />
-    <field id="sstgrad"      long_name="module of sst gradient"                                                                                 unit="degC/m"   />
-    <field id="sstgrad2"     long_name="square of module of sst gradient"                                                                       unit="degC2/m2" />
-    <field id="sbt"          long_name="sea bottom temperature"                                                                                 unit="degC"     />
-    <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" />
+    <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"     />
@@ -161,8 +158,8 @@
     <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"    />
-    <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)" />
-
-    <field id="taubot"       long_name="bottom stress module"                                                                                   unit="N/m2"     />
+    <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="taubot"       long_name="bottom stress module"   unit="N/m2"   grid_ref="grid_T_2D_inner" />
 
     <!-- Case EOS = TEOS-10 : output potential temperature -->
@@ -181,7 +178,7 @@
     <field id="mldr10_1max"  long_name="Max of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="maximum"                                                                          />
     <field id="mldr10_1min"  long_name="Min of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="minimum"                                                                          />
-    <field id="heatc"        long_name="Heat content vertically integrated"                 standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       />
-    <field id="saltc"        long_name="Salt content vertically integrated"                                                                                                                   unit="1e-3*kg/m2" />
-    <field id="saltc2"       long_name="Salt content vertically integrated"                                                                                                                   unit="1e-3*kg/m2" > saltc * saltc </field >
+    <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="1e-3.m"    grid_ref="grid_T_2D_inner" />
+    <field id="saltc2"       long_name="square of Salt content vertically integrated"                                                                                                         unit="1e-6.m2"   grid_ref="grid_T_2D_inner" />
 
     <!-- EOS -->
@@ -191,5 +188,5 @@
 
     <!-- Energy - horizontal divergence -->
-    <field id="sKE"          long_name="surface kinetic energy"  standard_name="specific_kinetic_energy_of_sea_water"   unit="m2/s2"  grid_ref="grid_T_2D" />
+    <field id="sKE"          long_name="surface kinetic energy"  standard_name="specific_kinetic_energy_of_sea_water"   unit="m2/s2"  grid_ref="grid_T_2D_inner" />
     <field id="hdiv"         long_name="horizontal divergence"                                                          unit="s-1"    grid_ref="grid_T_3D" />
 
@@ -600,5 +597,5 @@
 
     <field id="ssu"          long_name="ocean surface current along i-axis"                                                                 unit="m/s"                             />
-    <field id="sbu"          long_name="ocean bottom current along i-axis"                                                                  unit="m/s"                             />
+    <field id="sbu"          long_name="ocean bottom current along i-axis"                                                                  unit="m/s"  grid_ref="grid_U_2D_inner" />
     <field id="ubar"         long_name="ocean barotropic current along i-axis"                                                              unit="m/s"                             />
     <field id="uocetr_eff"   long_name="Effective ocean transport along i-axis"                 standard_name="ocean_volume_x_transport"    unit="m3/s"       grid_ref="grid_U_3D" />
@@ -639,8 +636,8 @@
 
     <!-- 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"                      />
-    <field id="u_heattr"      long_name="ocean eulerian heat transport along i-axis"                      standard_name="ocean_heat_x_transport"                         unit="W"                         />
-    <field id="u_salttr"      long_name="ocean eulerian salt transport along i-axis"                      standard_name="ocean_salt_x_transport"                         unit="1e-3*kg/s"                 />
+    <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_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"                 />
@@ -660,5 +657,5 @@
     <field id="voce_e3v"     long_name="ocean current along j-axis  (thickness weighted)"                                                   unit="m/s"        grid_ref="grid_V_3D"  > voce * e3v </field>
     <field id="ssv"          long_name="ocean surface current along j-axis"                                                                 unit="m/s"                             />
-    <field id="sbv"          long_name="ocean bottom current along j-axis"                                                                  unit="m/s"                             />
+    <field id="sbv"          long_name="ocean bottom current along j-axis"                                                                  unit="m/s"  grid_ref="grid_V_2D_inner" />
     <field id="vbar"         long_name="ocean barotropic current along j-axis"                                                              unit="m/s"                             />
     <field id="vocetr_eff"   long_name="Effective ocean transport along j-axis"                 standard_name="ocean_volume_y_transport"    unit="m3/s"       grid_ref="grid_V_3D" />
@@ -701,6 +698,6 @@
     <!-- 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_heattr"      long_name="ocean eulerian heat transport along j-axis"    standard_name="ocean_heat_y_transport"                          unit="W"                         />
-    <field id="v_salttr"      long_name="ocean eulerian salt transport along i-axis"    standard_name="ocean_salt_y_transport"                          unit="1e-3*kg/s"                 />
+    <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"                 />
@@ -775,10 +772,10 @@
     <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" />
-    <field id="ssrelvor"     long_name="surface relative vorticity"            standard_name="relative_vorticity"     unit="1/s"                   />
-    <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"                 />
-    <field id="ssabspotvor"  long_name="surface absolute potential vorticity"  standard_name="abspot_vorticity"       unit="1/m.s"                 />
-    <field id="ssEns"        long_name="surface enstrophy"                     standard_name="enstrophy"              unit="1/m2.s2"               />
+    <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_group>
 
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/cfgs/SHARED/grid_def_nemo.xml
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/cfgs/SHARED/grid_def_nemo.xml	(revision 15020)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/cfgs/SHARED/grid_def_nemo.xml	(revision 15022)
@@ -12,4 +12,7 @@
     <domain domain_ref="grid_T" />
   </grid>
+  <grid id="grid_T_2D_inner" >
+    <domain domain_ref="grid_T_inner" />
+  </grid>
   <!--  -->
   <grid id="grid_T_ncatice" >
@@ -17,4 +20,8 @@
     <axis axis_ref="ncatice" />
   </grid>
+  <grid id="grid_T_ncatice_inner" >
+    <domain domain_ref="grid_T_inner" />
+    <axis axis_ref="ncatice" />
+  </grid>
   <!--  -->
   <grid id="grid_T_3D" >
@@ -22,4 +29,8 @@
     <axis axis_ref="deptht" />
   </grid>
+  <grid id="grid_T_3D_inner" >
+    <domain domain_ref="grid_T_inner" />
+    <axis axis_ref="deptht" />
+  </grid>
   <!--  -->
   <grid id="grid_T_3DS" >
@@ -27,8 +38,15 @@
     <axis axis_ref="profsed" />
   </grid>
+  <grid id="grid_T_3DS_inner" >
+    <domain domain_ref="grid_T_inner" />
+    <axis axis_ref="profsed" />
+  </grid>
   <!--  -->
   <grid id="grid_U_2D" >
     <domain domain_ref="grid_U" />
   </grid>
+  <grid id="grid_U_2D_inner" >
+    <domain domain_ref="grid_U_inner" />
+  </grid>
   <!--  -->
   <grid id="grid_U_3D" >
@@ -36,8 +54,15 @@
     <axis axis_ref="depthu" />
   </grid>
+  <grid id="grid_U_3D_inner" >
+    <domain domain_ref="grid_U_inner" />
+    <axis axis_ref="depthu" />
+  </grid>
   <!--  -->
   <grid id="grid_V_2D" >
     <domain domain_ref="grid_V" />
   </grid>
+  <grid id="grid_V_2D_inner" >
+    <domain domain_ref="grid_V_inner" />
+  </grid>
   <!--  -->
   <grid id="grid_V_3D" >
@@ -45,16 +70,30 @@
     <axis axis_ref="depthv" />
   </grid>
+  <grid id="grid_V_3D_inner" >
+    <domain domain_ref="grid_V_inner" />
+    <axis axis_ref="depthv" />
+  </grid>
   <!--  -->
   <grid id="grid_W_2D" >
     <domain domain_ref="grid_W" />
   </grid>
+  <grid id="grid_W_2D_inner" >
+    <domain domain_ref="grid_W_inner" />
+  </grid>
   <!--  -->
   <grid id="grid_W_3D" >
     <domain domain_ref="grid_W" />
     <axis axis_ref="depthw" />
+  </grid>
+  <grid id="grid_W_3D_inner" >
+    <domain domain_ref="grid_W" />
+    <axis axis_ref="depthw_inner" />
   </grid>
   <!--  -->
   <grid id="grid_F_2D" >
     <domain domain_ref="grid_F" />
+  </grid>
+  <grid id="grid_F_2D_inner" >
+    <domain domain_ref="grid_F_inner" />
   </grid>
   <!--  -->
@@ -63,4 +102,8 @@
     <axis axis_ref="depthf" />
   </grid>
+  <grid id="grid_F_3D_inner" >
+    <domain domain_ref="grid_F_inner" />
+    <axis axis_ref="depthf" />
+  </grid>
   <!--  -->
   <grid id="grid_1point" >
@@ -77,5 +120,4 @@
   </grid>
   <!--  -->
-
 
   <grid id="grid_znl_T_2D">
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DIA/diawri.F90	(revision 15020)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DIA/diawri.F90	(revision 15022)
@@ -123,6 +123,6 @@
       REAL(wp)::   zztmp2, zztmpy   !   -      -
       REAL(wp)::   ze3
-      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
+      REAL(wp), DIMENSION(A2D(     0))     ::   z2d   ! 2D workspace
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   z3d   ! 3D workspace
       !!----------------------------------------------------------------------
       ! 
@@ -145,5 +145,5 @@
             z3d(:,:,jk) = gdept(:,:,jk,Kmm)
          END DO
-         CALL iom_put( "tpt_dep",     z3d(:,:,:) )
+         CALL iom_put( "tpt_dep", z3d )
       ENDIF
 
@@ -152,6 +152,9 @@
             z3d(:,:,jk) =  e3t(:,:,jk,Kmm)
          END DO
-         CALL iom_put( "e3t"     ,     z3d(:,:,:) )
-         CALL iom_put( "e3tdef"  , ( ( z3d(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
+         CALL iom_put( "e3t", z3d )
+         IF ( iom_use("e3tdef") ) THEN
+            z3d(:,:,:) = ( ( z3d(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100._wp * tmask(:,:,:) ) ** 2 
+            CALL iom_put( "e3tdef", z3d ) 
+         ENDIF
       ENDIF 
       IF ( iom_use("e3u") ) THEN                         ! time-varying e3u
@@ -159,5 +162,5 @@
             z3d(:,:,jk) =  e3u(:,:,jk,Kmm)
          END DO 
-         CALL iom_put( "e3u" , z3d(:,:,:) )
+         CALL iom_put( "e3u" , z3d )
       ENDIF
       IF ( iom_use("e3v") ) THEN                         ! time-varying e3v
@@ -165,5 +168,5 @@
             z3d(:,:,jk) =  e3v(:,:,jk,Kmm)
          END DO 
-         CALL iom_put( "e3v" , z3d(:,:,:) )
+         CALL iom_put( "e3v" , z3d )
       ENDIF
       IF ( iom_use("e3w") ) THEN                         ! time-varying e3w
@@ -171,5 +174,5 @@
             z3d(:,:,jk) =  e3w(:,:,jk,Kmm)
          END DO 
-         CALL iom_put( "e3w" , z3d(:,:,:) )
+         CALL iom_put( "e3w" , z3d )
       ENDIF
       IF ( iom_use("e3f") ) THEN                         ! time-varying e3f caution here at Kaa
@@ -177,11 +180,13 @@
             z3d(:,:,jk) =  e3f(:,:,jk)
          END DO
-         CALL iom_put( "e3f" , z3d(:,:,:) )
-      ENDIF
-
-      IF( ll_wd ) THEN                                   ! sea surface height (brought back to the reference used for wetting and drying)
-         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*ssmask(:,:) )
-      ELSE
-         CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height
+         CALL iom_put( "e3f" , z3d )
+      ENDIF
+
+      IF ( iom_use("ssh") ) THEN
+         IF( ll_wd ) THEN                                ! sea surface height (brought back to the reference used for wetting and drying)
+            CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*ssmask(:,:) )
+         ELSE
+            CALL iom_put( "ssh" ,  ssh(:,:,Kmm) )        ! sea surface height
+         ENDIF
       ENDIF
 
@@ -218,5 +223,5 @@
 
       IF ( iom_use("taubot") ) THEN                ! bottom stress
-         zztmp = rho0 * 0.25
+         zztmp = rho0 * 0.25_wp
          z2d(:,:) = 0._wp
          DO_2D( 0, 0, 0, 0 )
@@ -252,6 +257,6 @@
 
       !                                            ! vertical velocity
-      IF( ln_zad_Aimp ) THEN   ;   CALL iom_put( "woce", ww + wi )   ! explicit plus implicit parts
-      ELSE                     ;   CALL iom_put( "woce", ww )
+      IF( ln_zad_Aimp ) THEN   ;   IF( iom_use('woce') )   CALL iom_put( "woce", ww + wi )   ! explicit plus implicit parts
+      ELSE                     ;                           CALL iom_put( "woce", ww )
       ENDIF
 
@@ -276,33 +281,36 @@
       IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
 
-      IF ( iom_use("socegrad") .OR. iom_use("socegrad2") ) THEN
-         z3d(:,:,jpk) = 0.
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            zztmp  = ts(ji,jj,jk,jp_sal,Kmm)
-            zztmpx = (ts(ji+1,jj,jk,jp_sal,Kmm) - zztmp) * r1_e1u(ji,jj) + (zztmp - ts(ji-1,jj  ,jk,jp_sal,Kmm)) * r1_e1u(ji-1,jj)
-            zztmpy = (ts(ji,jj+1,jk,jp_sal,Kmm) - zztmp) * r1_e2v(ji,jj) + (zztmp - ts(ji  ,jj-1,jk,jp_sal,Kmm)) * r1_e2v(ji,jj-1)
-            z3d(ji,jj,jk) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
-               &                 * umask(ji,jj,jk) * umask(ji-1,jj,jk) * vmask(ji,jj,jk) * umask(ji,jj-1,jk)
-         END_3D
-         CALL iom_put( "socegrad2",  z3d )          ! square of module of sal gradient
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            z3d(ji,jj,jk) = SQRT( z3d(ji,jj,jk) )
-         END_3D
-         CALL iom_put( "socegrad" ,  z3d )          ! module of sal gradient
+      IF ( iom_use("sssgrad") .OR. iom_use("sssgrad2") ) THEN
+         DO_2D( 0, 0, 0, 0 )                       ! sss gradient
+            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 )   &
+               &                 * 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
+         IF ( iom_use("sssgrad") ) THEN
+            DO_2D( 0, 0, 0, 0 )
+               z2d(ji,jj) = SQRT( z2d(ji,jj) )
+            END_2D
+            CALL iom_put( "sssgrad",  z2d )        ! module of sss gradient
+         ENDIF
       ENDIF
          
       IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
-         DO_2D( 0, 0, 0, 0 )                                 ! sst gradient
+         DO_2D( 0, 0, 0, 0 )                       ! sst gradient
             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 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
-               &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
+            z2d(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
-         DO_2D( 0, 0, 0, 0 )
-            z2d(ji,jj) = SQRT( z2d(ji,jj) )
-         END_2D
-         CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient
+         IF ( iom_use("sstgrad") ) THEN
+            DO_2D( 0, 0, 0, 0 )
+               z2d(ji,jj) = SQRT( z2d(ji,jj) )
+            END_2D
+            CALL iom_put( "sstgrad",  z2d )        ! module of sst gradient
+         ENDIF
       ENDIF
          
@@ -321,5 +329,5 @@
             z2d(ji,jj) = z2d(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 * z2d )       ! vertically integrated salt content (PSU*kg/m2)
       ENDIF
       !
@@ -329,5 +337,5 @@
             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)
          END_3D
-         CALL iom_put( "salt2c", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
+         CALL iom_put( "salt2c", rho0 * z2d )      ! vertically integrated salt content (PSU*kg/m2)
       ENDIF
       !
@@ -335,7 +343,7 @@
          z3d(:,:,jpk) = 0._wp 
          DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            zztmpx = 0.5 * ( uu(ji-1,jj  ,jk,Kmm) + uu(ji,jj,jk,Kmm) )
-            zztmpy = 0.5 * ( vv(ji  ,jj-1,jk,Kmm) + vv(ji,jj,jk,Kmm) )
-            z3d(ji,jj,jk) = 0.5 * ( zztmpx*zztmpx + zztmpy*zztmpy )
+            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 )
          END_3D
          CALL iom_put( "ke", z3d )                 ! kinetic energy
@@ -345,8 +353,8 @@
             z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
          END_3D
-         CALL iom_put( "ke_int", z2d )   ! vertically integrated kinetic energy
-      ENDIF
-      !
-      IF ( iom_use("sKE") ) THEN                        ! surface kinetic energy at T point
+         CALL iom_put( "ke_int", z2d )             ! vertically integrated kinetic energy
+      ENDIF
+      !
+      IF ( iom_use("sKE") ) THEN                   ! surface kinetic energy at T point
          z2d(:,:) = 0._wp
          DO_2D( 0, 0, 0, 0 )
@@ -357,10 +365,9 @@
                &                 * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * ssmask(ji,jj)
          END_2D
-         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
          IF ( iom_use("sKE" ) )  CALL iom_put( "sKE" , z2d )   
       ENDIF
       !    
-      IF ( iom_use("ssKEf") ) THEN                        ! surface kinetic energy at F point
-         z2d(:,:) = 0._wp                                ! CAUTION : only valid in SWE, not with bathymetry
+      IF ( iom_use("ssKEf") ) THEN                 ! surface kinetic energy at F point
+         z2d(:,:) = 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)  &
@@ -370,86 +377,66 @@
                &                 * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj)
          END_2D
-         CALL lbc_lnk( 'diawri', z2d, 'F', 1. )
          CALL iom_put( "ssKEf", z2d )                     
       ENDIF
       !
-      CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence
-
-      IF ( iom_use("relvor") .OR. iom_use("absvor") .OR. iom_use("potvor") ) THEN
+      CALL iom_put( "hdiv", hdiv )                 ! Horizontal divergence
+      !
+      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
          
-         z3d(:,:,jpk) = 0._wp 
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            z3d(ji,jj,jk) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,jk,Kmm) - e2v(ji,jj) * vv(ji,jj,jk,Kmm)    &
-               &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,jk,Kmm) + e1u(ji,jj) * uu(ji,jj,jk,Kmm)  ) * r1_e1e2f(ji,jj)
-         END_3D
-         CALL iom_put( "relvor", z3d )                  ! relative vorticity
-
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            z3d(ji,jj,jk) = ff_f(ji,jj) + z3d(ji,jj,jk) 
-         END_3D
-         CALL iom_put( "absvor", z3d )                  ! absolute vorticity
-
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            ze3  = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
-               &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
-            IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3
-            ELSE                      ;   ze3 = 0._wp
-            ENDIF
-            z3d(ji,jj,jk) = ze3 * z3d(ji,jj,jk) 
-         END_3D
-         CALL iom_put( "potvor", z3d )                  ! potential vorticity
-
-      ENDIF
-      !
-      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
-         z3d(:,:,jpk) = 0.e0
-         z2d(:,:) = 0.e0
+         z3d(:,:,jpk) = 0._wp
          DO jk = 1, jpkm1
             z3d(:,:,jk) = rho0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk)
-            z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
          END DO
-         CALL iom_put( "u_masstr"     , z3d )         ! mass transport in i-direction
-         CALL iom_put( "u_masstr_vint", z2d )         ! mass transport in i-direction vertical sum
-      ENDIF
-      
-      IF( iom_use("u_heattr") ) THEN
-         z2d(:,:) = 0._wp 
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            z2d(ji,jj) = z2d(ji,jj) + z3d(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", 0.5*rcp * z2d )    ! heat transport in i-direction
-      ENDIF
-
-      IF( iom_use("u_salttr") ) THEN
-         z2d(:,:) = 0.e0 
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            z2d(ji,jj) = z2d(ji,jj) + z3d(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", 0.5 * z2d )        ! heat transport in i-direction
-      ENDIF
-
+         CALL iom_put( "u_masstr"     , z3d )      ! mass transport in i-direction
+         
+         IF( iom_use("u_masstr_vint") ) THEN
+            z2d(:,:) = 0._wp 
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+               z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
+            END_3D
+            CALL iom_put( "u_masstr_vint", z2d )   ! mass transport in i-direction vertical sum
+         ENDIF
+         IF( iom_use("u_heattr") ) THEN
+            z2d(:,:) = 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) )
+            END_3D
+            CALL iom_put( "u_heattr", z2d )        ! heat transport in i-direction
+         ENDIF
+         IF( iom_use("u_salttr") ) THEN
+            z2d(:,:) = 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) )
+            END_3D
+            CALL iom_put( "u_salttr", z2d )        ! heat transport in i-direction
+         ENDIF
+         
+      ENDIF
       
       IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
-         z3d(:,:,jpk) = 0.e0
+         
+         z3d(:,:,jpk) = 0._wp
          DO jk = 1, jpkm1
             z3d(:,:,jk) = rho0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk)
          END DO
-         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
-      ENDIF
-      
-      IF( iom_use("v_heattr") ) THEN
-         z2d(:,:) = 0.e0 
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            z2d(ji,jj) = z2d(ji,jj) + z3d(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", 0.5*rcp * z2d )    !  heat transport in j-direction
-      ENDIF
-
-      IF( iom_use("v_salttr") ) THEN
-         z2d(:,:) = 0._wp 
-         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            z2d(ji,jj) = z2d(ji,jj) + z3d(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", 0.5 * z2d )        !  heat transport in j-direction
+         CALL iom_put( "v_masstr", z3d )           ! mass transport in j-direction
+         
+         IF( iom_use("v_heattr") ) THEN
+            z2d(:,:) = 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) )
+            END_3D
+            CALL iom_put( "v_heattr", z2d )        !  heat transport in j-direction
+         ENDIF
+         IF( iom_use("v_salttr") ) THEN
+            z2d(:,:) = 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) )
+            END_3D
+            CALL iom_put( "v_salttr", z2d )        !  heat transport in j-direction
+         ENDIF
+
       ENDIF
 
@@ -457,65 +444,64 @@
          z2d(:,:) = 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)
+            z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
          END_3D
-         CALL iom_put( "tosmint", rho0 * z2d )        ! Vertical integral of temperature
+         CALL iom_put( "tosmint", z2d )            ! Vertical integral of temperature
       ENDIF
       IF( iom_use("somint") ) THEN
-         z2d(:,:)=0._wp
+         z2d(:,:) = 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)
+            z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
          END_3D
-         CALL iom_put( "somint", rho0 * z2d )         ! Vertical integral of salinity
-      ENDIF
-
-      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2)
-      !
+         CALL iom_put( "somint", z2d )             ! Vertical integral of salinity
+      ENDIF
+
+      CALL iom_put( "bn2", rn2 )                   ! Brunt-Vaisala buoyancy frequency (N^2)
       
-      IF (ln_dia25h)   CALL dia_25h( kt, Kmm )        ! 25h averaging
+      IF (ln_dia25h)   CALL dia_25h( kt, Kmm )     ! 25h averaging
       
       ! Output of surface vorticity terms
-      IF ( iom_use("ssrelvor")    .OR. iom_use("ssplavor")    .OR.   &
-         & iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") .OR.   &
-         & iom_use("ssEns")                                        ) THEN
+      !
+      CALL iom_put( "ssplavor", ff_f )             ! planetary vorticity ( f )
+      !
+      IF ( iom_use("ssrelvor")    .OR. iom_use("ssEns")    .OR.   &
+         & iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
          !
          z2d(:,:) = 0._wp 
-         ze3 = 0._wp 
-         DO_2D( 1, 0, 1, 0 )
+         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)    &
             &              - 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 lbc_lnk( 'diawri', z2d, 'F', 1. )
-         CALL iom_put( "ssrelvor", z2d )                  ! relative vorticity ( zeta ) 
+         CALL iom_put( "ssrelvor", z2d )           ! relative vorticity ( zeta ) 
          !
-         CALL iom_put( "ssplavor", ff_f )                 ! planetary vorticity ( f )
-         !
-         DO_2D( 1, 0, 1, 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)
-            IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3
-            ELSE                      ;   ze3 = 0._wp
+         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)
+               IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3
+               ELSE                      ;   ze3 = 0._wp
+               ENDIF
+               z2d(ji,jj) = ze3 * z2d(ji,jj) 
+            END_2D
+            CALL iom_put( "ssrelpotvor", z2d )     ! 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)
+                  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) 
+               END_2D
+               CALL iom_put( "ssabspotvor", z2d )  ! 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) 
+                  END_2D
+                  CALL iom_put( "ssEns", z2d )     ! potential enstrophy ( 1/2*q2 )
+               ENDIF
             ENDIF
-            z2d(ji,jj) = ze3 * z2d(ji,jj) 
-         END_2D
-         CALL lbc_lnk( 'diawri', z2d, 'F', 1. )
-         CALL iom_put( "ssrelpotvor", z2d )                  ! relative potential vorticity (zeta/h)
-         !
-         DO_2D( 1, 0, 1, 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)
-            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) 
-         END_2D
-         CALL lbc_lnk( 'diawri', z2d, 'F', 1. )
-         CALL iom_put( "ssabspotvor", z2d )                  ! absolute potential vorticity ( q )
-         !
-         DO_2D( 1, 0, 1, 0 )  
-            z2d(ji,jj) = 0.5_wp * z2d(ji,jj)  * z2d(ji,jj) 
-         END_2D
-         CALL lbc_lnk( 'diawri', z2d, 'F', 1. )
-         CALL iom_put( "ssEns", z2d )                        ! potential enstrophy ( 1/2*q2 )
-         !
+         ENDIF
       ENDIF
 
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DYN/dynspg_ts.F90	(revision 15020)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DYN/dynspg_ts.F90	(revision 15022)
@@ -499,6 +499,10 @@
          !      
          !                             ! resulting flux at mid-step (not over the full domain)
-         zhU(1:jpim1,1:jpj  ) = e2u(1:jpim1,1:jpj  ) * ua_e(1:jpim1,1:jpj  ) * zhup2_e(1:jpim1,1:jpj  )   ! not jpi-column
-         zhV(1:jpi  ,1:jpjm1) = e1v(1:jpi  ,1:jpjm1) * va_e(1:jpi  ,1:jpjm1) * zhvp2_e(1:jpi  ,1:jpjm1)   ! not jpj-row
+         DO_2D( 1, 0, 1, 1 )   ! not jpi-column
+            zhU(ji,jj) = e2u(ji,jj) * ua_e(ji,jj) * zhup2_e(ji,jj)
+         END_2D
+         DO_2D( 1, 1, 1, 0 )   ! not jpj-row
+            zhV(ji,jj) = e1v(ji,jj) * va_e(ji,jj) * zhvp2_e(ji,jj)
+         END_2D
          !
 #if defined key_agrif
@@ -538,6 +542,10 @@
          ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 
          IF ( ln_wd_dl_bc ) THEN
-            zuwdav2(1:jpim1,1:jpj  ) = zuwdav2(1:jpim1,1:jpj  ) + za2 * zuwdmask(1:jpim1,1:jpj  )   ! not jpi-column
-            zvwdav2(1:jpi  ,1:jpjm1) = zvwdav2(1:jpi  ,1:jpjm1) + za2 * zvwdmask(1:jpi  ,1:jpjm1)   ! not jpj-row
+            DO_2D( 1, 0, 1, 1 )   ! not jpi-column
+               zuwdav2(ji,jj) = zuwdav2(ji,jj) + za2 * zuwdmask(ji,jj)
+            END_2D
+            DO_2D( 1, 1, 1, 0 )   ! not jpj-row
+               zvwdav2(ji,jj) = zvwdav2(ji,jj) + za2 * zvwdmask(ji,jj)
+            END_2D
          END IF
          !
@@ -579,6 +587,8 @@
          IF( ln_wd_il ) THEN        ! W/D : gravity filters applied on pressure gradient
             CALL wad_spg( zsshp2_e, zcpx, zcpy )   ! Calculating W/D gravity filters
-            zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1)
-            zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1)
+            DO_2D( 0, 0, 0, 0 )
+               zu_spg(ji,jj) = zu_spg(ji,jj) * zcpx(ji,jj)
+               zv_spg(ji,jj) = zv_spg(ji,jj) * zcpy(ji,jj)
+            END_2D
          ENDIF
          !
@@ -670,8 +680,10 @@
        
          IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only)
-            hu_e (2:jpim1,2:jpjm1) =    hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)
-            hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / (  hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1)  )
-            hv_e (2:jpim1,2:jpjm1) =    hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)
-            hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / (  hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1)  )
+            DO_2D( 0, 0, 0, 0 )
+               hu_e (ji,jj) =    hu_0(ji,jj) + zsshu_a(ji,jj)
+               hur_e(ji,jj) = ssumask(ji,jj) / (  hu_e(ji,jj) + 1._wp - ssumask(ji,jj)  )
+               hv_e (ji,jj) =    hv_0(ji,jj) + zsshv_a(ji,jj)
+               hvr_e(ji,jj) = ssvmask(ji,jj) / (  hv_e(ji,jj) + 1._wp - ssvmask(ji,jj)  )
+            END_2D
          ENDIF
          !
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/IOM/iom.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/IOM/iom.F90	(revision 15020)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/IOM/iom.F90	(revision 15022)
@@ -2325,4 +2325,8 @@
       CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = -nn_hls, data_ni=jpi, data_jbegin = -nn_hls, data_nj=jpj)
 
+      CALL iom_set_domain_attr("grid_"//cdgrd//"_inner", ni_glo = Ni0glo, nj_glo = Nj0glo,   &
+         &                     ibegin = mig0(Nis0) - 1, jbegin = mjg0(Njs0) - 1, ni = Ni_0, nj = Nj_0)
+      CALL iom_set_domain_attr("grid_"//cdgrd//"_inner", data_dim=2, data_ibegin = 0, data_ni=Ni_0, data_jbegin = 0, data_nj=Nj_0)
+
       IF( ln_tile ) THEN
          DO jn = 1, nijtile
@@ -2334,4 +2338,9 @@
          ! Tile_[ij]begin are defined with respect to the processor data domain, so data_[ij]begin is added
          CALL iom_set_domain_attr("grid_"//cdgrd, ntiles=nijtile,                                     &
+            & tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
+            & tile_ni=ini(:), tile_nj=inj(:),                                                         &
+            & tile_data_ibegin=idb(:), tile_data_jbegin=idb(:),                                       &
+            & tile_data_ni=ini(:) - 2 * idb(:), tile_data_nj=inj(:) - 2 * idb(:))
+         CALL iom_set_domain_attr("grid_"//cdgrd//"_inner", ntiles=nijtile,                           &
             & tile_ibegin=ntsi_a(1:nijtile) + idb(:) - 1, tile_jbegin=ntsj_a(1:nijtile) + idb(:) - 1, &
             & tile_ni=ini(:), tile_nj=inj(:),                                                         &
@@ -2356,6 +2365,8 @@
          END SELECT
          !
-         CALL iom_set_domain_attr( "grid_"//cdgrd       , mask = RESHAPE(zmask(Nis0:Nie0,Njs0:Nje0,1),(/Ni_0*Nj_0    /)) /= 0. )
-         CALL iom_set_grid_attr  ( "grid_"//cdgrd//"_3D", mask = RESHAPE(zmask(Nis0:Nie0,Njs0:Nje0,:),(/Ni_0,Nj_0,jpk/)) /= 0. )
+         CALL iom_set_domain_attr( "grid_"//cdgrd             , mask=RESHAPE(zmask(Nis0:Nie0,Njs0:Nje0,1),(/Ni_0*Nj_0    /)) /= 0. )
+         CALL iom_set_grid_attr  ( "grid_"//cdgrd//"_3D"      , mask=RESHAPE(zmask(Nis0:Nie0,Njs0:Nje0,:),(/Ni_0,Nj_0,jpk/)) /= 0. )
+         CALL iom_set_domain_attr( "grid_"//cdgrd//"_inner"   , mask=RESHAPE(zmask(Nis0:Nie0,Njs0:Nje0,1),(/Ni_0*Nj_0    /)) /= 0. )
+         CALL iom_set_grid_attr  ( "grid_"//cdgrd//"_3D_inner", mask=RESHAPE(zmask(Nis0:Nie0,Njs0:Nje0,:),(/Ni_0,Nj_0,jpk/)) /= 0. )
       ENDIF
       !
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/tests/CANAL/EXPREF/file_def_nemo-oce.xml
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/tests/CANAL/EXPREF/file_def_nemo-oce.xml	(revision 15020)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/tests/CANAL/EXPREF/file_def_nemo-oce.xml	(revision 15022)
@@ -15,8 +15,8 @@
 	  <field field_ref="soce" /> 
 	  <field field_ref="ssh"  />
-	  <field field_ref="socegrad"  />
+	  <field field_ref="sssgrad"  />
 	  <field field_ref="ke_int"  />
-	  <field field_ref="relvor"  />
-	  <field field_ref="potvor"  />
+	  <field field_ref="ssrelvor"  />
+	  <field field_ref="ssrelpotvor"  />
 	  <field field_ref="saltc"  />
 	  <field field_ref="salt2c"  />
