diff --git a/cfgs/AGRIF_DEMO/EXPREF/file_def_nemo-oce.xml b/cfgs/AGRIF_DEMO/EXPREF/file_def_nemo-oce.xml
index 0d0368d72066347fcc66532477e50d40941ecf95..6cb04ff08edd17dbda063ce517c65ec3582c7805 100644
--- a/cfgs/AGRIF_DEMO/EXPREF/file_def_nemo-oce.xml
+++ b/cfgs/AGRIF_DEMO/EXPREF/file_def_nemo-oce.xml
@@ -34,6 +34,8 @@
 	  <field field_ref="qt_oce"       name="qt_oce"   />
 	  <field field_ref="saltflx"      name="sfx"      />
 	  <field field_ref="taum"         name="taum"     />
+	  <field field_ref="utau"         name="tauuo"    />
+	  <field field_ref="vtau"         name="tauvo"    />
 	  <field field_ref="wspd"         name="windsp"   />
 	  <field field_ref="precip"       name="precip"   />
 	  <!-- ice and snow -->
@@ -44,7 +46,6 @@
 	  <field field_ref="e3u" />
 	  <field field_ref="ssu"          name="uos"      />
 	  <field field_ref="uoce"         name="uo"       operation="instant" freq_op="5d" > @uoce_e3u / @e3u </field>
-	  <field field_ref="utau"         name="tauuo"    />
           <field field_ref="uocetr_eff"   name="uocetr_eff"  />
           <!-- available with diaar5 -->
 	  <field field_ref="u_masstr"     name="vozomatr" />
@@ -56,7 +57,6 @@
 	  <field field_ref="e3v" />
 	  <field field_ref="ssv"          name="vos"      />
 	  <field field_ref="voce"         name="vo"       operation="instant" freq_op="5d" > @voce_e3v / @e3v </field>
-	  <field field_ref="vtau"         name="tauvo"    />
           <field field_ref="vocetr_eff"   name="vocetr_eff"  />
           <!-- available with diaar5 -->
 	  <field field_ref="v_masstr"     name="vomematr" />
diff --git a/cfgs/AMM12/EXPREF/file_def_nemo-oce.xml b/cfgs/AMM12/EXPREF/file_def_nemo-oce.xml
index 776def802a18a6cc9f7ee4e9dd83f95132c4f5c7..d8c7291419b3cf1e34be1995cfffc2b50512c4eb 100644
--- a/cfgs/AMM12/EXPREF/file_def_nemo-oce.xml
+++ b/cfgs/AMM12/EXPREF/file_def_nemo-oce.xml
@@ -100,6 +100,8 @@
 	  <field field_ref="qsr"          name="rsntds" />
 	  <field field_ref="qt"           name="tohfls" />
 	  <field field_ref="taum"     />
+	  <field field_ref="utau"         name="tauuo" />
+	  <field field_ref="vtau"         name="tauvo" />
 	  <field field_ref="mldkz5"   />
 	  <field field_ref="mldr10_1" />
 	</file>
@@ -107,13 +109,11 @@
 	<file id="file5" name_suffix="_grid_U" description="ocean U grid variables" >
 	  <field field_ref="uoce"         name="uo"    />
 	  <field field_ref="ssu"          name="uos"   />
-	  <field field_ref="utau"         name="tauuo" />
 	</file>
 	
 	<file id="file6" name_suffix="_grid_V" description="ocean V grid variables" >
 	  <field field_ref="voce"         name="vo"    />
 	  <field field_ref="ssv"          name="vos"   />
-	  <field field_ref="vtau"         name="tauvo" />
 	</file>
 	
 	<file id="file7" name_suffix="_grid_W" description="ocean W grid variables" >
diff --git a/cfgs/GYRE_PISCES/EXPREF/file_def_nemo.xml b/cfgs/GYRE_PISCES/EXPREF/file_def_nemo.xml
index e24e2271a290446e9abcd4371329e210b83e0ad5..57144055d7774e652660e39313bf7036a4573a88 100644
--- a/cfgs/GYRE_PISCES/EXPREF/file_def_nemo.xml
+++ b/cfgs/GYRE_PISCES/EXPREF/file_def_nemo.xml
@@ -32,16 +32,16 @@
 	  <field field_ref="qt"           name="sohefldo"  />
 	  <field field_ref="mldr10_1"     name="somxl010"  />
 	  <field field_ref="mldkz5"       name="somixhgt"  />
+          <field field_ref="utau"         name="sozotaux"  />
+          <field field_ref="vtau"         name="sometauy"  /> 
         </file>
 	
         <file id="file2" name_suffix="_grid_U" description="ocean U grid variables" >
           <field field_ref="uoce"         name="vozocrtx"  />
-          <field field_ref="utau"         name="sozotaux"  />
         </file>
 	
         <file id="file3" name_suffix="_grid_V" description="ocean V grid variables" >
           <field field_ref="voce"         name="vomecrty"  /> 
-          <field field_ref="vtau"         name="sometauy"  /> 
         </file>
 	
         <file id="file4" name_suffix="_grid_W" description="ocean W grid variables" >
diff --git a/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-oce.xml b/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-oce.xml
index 6d6df82d2902518fb8e6e7e45a9e1d182726462d..9ba0c304fe37810dc0178058ac87ca15dff327df 100644
--- a/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-oce.xml
+++ b/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-oce.xml
@@ -30,6 +30,8 @@
 	  <field field_ref="qt_oce"       name="qt_oce"   />
 	  <field field_ref="saltflx"      name="sfx"      />
 	  <field field_ref="taum"         name="taum"     />
+	  <field field_ref="utau"         name="tauuo"    />
+	  <field field_ref="vtau"         name="tauvo"    />
 	  <field field_ref="wspd"         name="windsp"   />
 	  <field field_ref="precip"       name="precip"   />
 	  <!-- ice and snow -->
@@ -40,14 +42,12 @@
 	  <field field_ref="e3u" />
 	  <field field_ref="ssu"          name="uos"      />
 	  <field field_ref="uoce"         name="uo"       />
-	  <field field_ref="utau"         name="tauuo"    />
 	</file>
 	
 	<file id="file13" name_suffix="_grid_V" description="ocean V grid variables" >
 	  <field field_ref="e3v" />
 	  <field field_ref="ssv"          name="vos"      />
 	  <field field_ref="voce"         name="vo"       />
-	  <field field_ref="vtau"         name="tauvo"    />
 	</file>
 	
 	<file id="file14" name_suffix="_grid_ABL" description="ABL grid variables" >
diff --git a/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-oce.xml b/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-oce.xml
index c043f56fe5e9ae6fc252ae34d78b6c2660381ec9..04590b55f29789842576a1dff286ad9ec572d172 100644
--- a/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-oce.xml
+++ b/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-oce.xml
@@ -34,6 +34,8 @@
 	  <field field_ref="qt_oce"       name="qt_oce"   />
 	  <field field_ref="saltflx"      name="sfx"      />
 	  <field field_ref="taum"         name="taum"     />
+	  <field field_ref="utau"         name="tauuo"    />
+	  <field field_ref="vtau"         name="tauvo"    />
 	  <field field_ref="wspd"         name="windsp"   />
 	  <field field_ref="precip"       name="precip"   />
 	  <!-- ice and snow -->
@@ -44,7 +46,6 @@
 	  <field field_ref="e3u" />
 	  <field field_ref="ssu"          name="uos"      />
 	  <field field_ref="uoce"         name="uo"       operation="instant" freq_op="5d" > @uoce_e3u / @e3u </field>
-	  <field field_ref="utau"         name="tauuo"    />
           <field field_ref="uocetr_eff"   name="uocetr_eff"  />
           <!-- available with diaar5 -->
 	  <field field_ref="u_masstr"     name="vozomatr" />
@@ -56,7 +57,6 @@
 	  <field field_ref="e3v" />
 	  <field field_ref="ssv"          name="vos"      />
 	  <field field_ref="voce"         name="vo"       operation="instant" freq_op="5d" > @voce_e3v / @e3v </field>
-	  <field field_ref="vtau"         name="tauvo"    />
           <field field_ref="vocetr_eff"   name="vocetr_eff"  />
           <!-- available with diaar5 -->
 	  <field field_ref="v_masstr"     name="vomematr" />
diff --git a/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-oce.xml b/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-oce.xml
index 5c75b89204e6b01528a9f6d0310d71949f58b0b4..2cdea4e2397740c0b0d8933093314a70e98891ef 100644
--- a/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-oce.xml
+++ b/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-oce.xml
@@ -34,6 +34,8 @@
 	  <field field_ref="qt_oce"       name="qt_oce"   />
 	  <field field_ref="saltflx"      name="sfx"      />
 	  <field field_ref="taum"         name="taum"     />
+	  <field field_ref="utau"         name="tauuo"    />
+	  <field field_ref="vtau"         name="tauvo"    />
 	  <field field_ref="wspd"         name="windsp"   />
 	  <field field_ref="precip"       name="precip"   />
 	  <!-- ice and snow -->
@@ -44,7 +46,6 @@
 	  <field field_ref="e3u" />
 	  <field field_ref="ssu"          name="uos"      />
 	  <field field_ref="uoce"         name="uo"       operation="instant" freq_op="5d" > @uoce_e3u / @e3u </field>
-	  <field field_ref="utau"         name="tauuo"    />
           <field field_ref="uocetr_eff"   name="uocetr_eff"  />
           <!-- available with diaar5 -->
 	  <field field_ref="u_masstr"     name="vozomatr" />
@@ -56,7 +57,6 @@
 	  <field field_ref="e3v" />
 	  <field field_ref="ssv"          name="vos"      />
 	  <field field_ref="voce"         name="vo"       operation="instant" freq_op="5d" > @voce_e3v / @e3v </field>
-	  <field field_ref="vtau"         name="tauvo"    />
           <field field_ref="vocetr_eff"   name="vocetr_eff"  />
           <!-- available with diaar5 -->
 	  <field field_ref="v_masstr"     name="vomematr" />
diff --git a/cfgs/SHARED/field_def_nemo-oce.xml b/cfgs/SHARED/field_def_nemo-oce.xml
index ea75b3ec77a6eceb29995558d315213dcc6df9de..e5bff5cc7052509381db14633c700780af2adc61 100644
--- a/cfgs/SHARED/field_def_nemo-oce.xml
+++ b/cfgs/SHARED/field_def_nemo-oce.xml
@@ -432,6 +432,8 @@ that are available in the tidal-forcing implementation (see
       <field id="erp"          long_name="Surface Water Flux: Damping"           standard_name="water_flux_out_of_sea_water_due_to_newtonian_relaxation"              unit="kg/m2/s"                        />
       <field id="taum"         long_name="wind stress module"                    standard_name="magnitude_of_surface_downward_stress"                                 unit="N/m2"                           />
       <field id="wspd"         long_name="wind speed module"                     standard_name="wind_speed"                                                           unit="m/s"                            />
+      <field id="utau"         long_name="Wind Stress along i-axis"              standard_name="surface_downward_x_stress"                                            unit="N/m2"                           />
+      <field id="vtau"         long_name="Wind Stress along j-axis"              standard_name="surface_downward_y_stress"                                            unit="N/m2"                           />
 
       <!-- * variable relative to atmospheric pressure forcing : available with ln_apr_dyn -->
       <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"        />
@@ -583,7 +585,6 @@ that are available in the tidal-forcing implementation (see
     <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="utau"          long_name="Wind Stress along i-axis"                               standard_name="surface_downward_x_stress"   unit="N/m2"                            />
     <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"/>
@@ -651,7 +652,6 @@ that are available in the tidal-forcing implementation (see
     <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" />
-    <field id="vtau"         long_name="Wind Stress along j-axis"                               standard_name="surface_downward_y_stress"   unit="N/m2"                            />
     <field id="voce"         long_name="ocean current along j-axis"                             standard_name="sea_water_y_velocity"        unit="m/s"        grid_ref="grid_V_3D" />
     <field id="voce_e3v"     long_name="ocean current along j-axis  (thickness weighted)"                                                   unit="m/s"        grid_ref="grid_V_3D"  > voce * e3v </field>
     <field id="ssv"          long_name="ocean surface current along j-axis"                                                                 unit="m/s"                             />
@@ -1189,6 +1189,8 @@ that are available in the tidal-forcing implementation (see
     <field field_ref="qsr"          name="rsntds"   long_name="surface_net_downward_shortwave_flux"           />
     <field field_ref="qt"           name="tohfls"   long_name="surface_net_downward_total_heat_flux"          />
     <field field_ref="taum"                                                                                   />
+    <field field_ref="utau"         name="tauuo"    long_name="surface_downward_x_stress"                     />
+    <field field_ref="vtau"         name="tauvo"    long_name="surface_downward_y_stress"                     />
     <field field_ref="20d"                                                                                    />
     <field field_ref="mldkz5"                                                                                 />
     <field field_ref="mldr10_1"                                                                               />
@@ -1203,12 +1205,10 @@ that are available in the tidal-forcing implementation (see
 
   <field_group id="groupU" >
     <field field_ref="uoce"         name="uo"      long_name="sea_water_x_velocity"      />
-    <field field_ref="utau"         name="tauuo"   long_name="surface_downward_x_stress" />
   </field_group>
 
   <field_group id="groupV" >
     <field field_ref="voce"         name="vo"      long_name="sea_water_y_velocity"      />
-    <field field_ref="vtau"         name="tauvo"   long_name="surface_downward_y_stress" />
   </field_group>
 
   <field_group id="groupW" >
diff --git a/cfgs/SHARED/namelist_ref b/cfgs/SHARED/namelist_ref
index ae9d42ee79e125df420a6142699d2a3f4de78754..dc77ca9a56b679a9a8b7c1dd80edbac11191ea24 100644
--- a/cfgs/SHARED/namelist_ref
+++ b/cfgs/SHARED/namelist_ref
@@ -369,7 +369,7 @@
 !***  receive  ***
    sn_rcv_w10m   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
    sn_rcv_taumod =   'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
-   sn_rcv_tau    =   'oce only'             ,    'no'    , 'cartesian' , 'eastward-northward' ,  'U,V'
+   sn_rcv_tau    =   'oce only'             ,    'no'    , 'cartesian' , 'eastward-northward' ,   ''
    sn_rcv_dqnsdt =   'coupled'              ,    'no'    ,     ''      ,         ''           ,   ''
    sn_rcv_qsr    =   'oce and ice'          ,    'no'    ,     ''      ,         ''           ,   ''
    sn_rcv_qns    =   'oce and ice'          ,    'no'    ,     ''      ,         ''           ,   ''
@@ -380,21 +380,22 @@
    sn_rcv_iceflx =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
    sn_rcv_mslp   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
    sn_rcv_ts_ice =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_qtrice =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
    sn_rcv_isf    =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
    sn_rcv_icb    =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
-   sn_rcv_hsig   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'T'
-   sn_rcv_phioc  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'T'
-   sn_rcv_sdrfx  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'T'
-   sn_rcv_sdrfy  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'T'
-   sn_rcv_wper   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'T'
-   sn_rcv_wnum   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'T'
-   sn_rcv_wstrf  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'T'
-   sn_rcv_wdrag  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'T'
-   sn_rcv_charn  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'T'
-   sn_rcv_taw    =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'U,V'
-   sn_rcv_bhd    =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'T'
-   sn_rcv_tusd   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'T'
-   sn_rcv_tvsd   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'T'
+   sn_rcv_hsig   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_phioc  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_sdrfx  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_sdrfy  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_wper   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_wnum   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_wstrf  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_wdrag  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_charn  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_taw    =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_bhd    =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_tusd   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_tvsd   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
 /
 !-----------------------------------------------------------------------
 &namsbc_sas    !   Stand-Alone Surface module: ocean data               (SAS_SRC  only)
diff --git a/cfgs/SPITZ12/EXPREF/file_def_nemo-oce.xml b/cfgs/SPITZ12/EXPREF/file_def_nemo-oce.xml
index aa6deff58516c5d377db0c0f9193c1ad4260bb4f..ee0982291a1ef1a38d24837c6a912f81346fdf84 100644
--- a/cfgs/SPITZ12/EXPREF/file_def_nemo-oce.xml
+++ b/cfgs/SPITZ12/EXPREF/file_def_nemo-oce.xml
@@ -34,6 +34,8 @@
 	  <field field_ref="qt_oce"       name="qt_oce"   />
 	  <field field_ref="saltflx"      name="sfx"      />
 	  <field field_ref="taum"         name="taum"     />
+	  <field field_ref="utau"         name="tauuo"    />
+	  <field field_ref="vtau"         name="tauvo"    />
 	  <field field_ref="wspd"         name="windsp"   />
 	  <field field_ref="precip"       name="precip"   />
 	  <!-- ice and snow -->
@@ -44,7 +46,6 @@
 	  <field field_ref="e3u" />
 	  <field field_ref="ssu"          name="uos"      />
 	  <field field_ref="uoce"         name="uo"       operation="instant" freq_op="5d" > @uoce_e3u / @e3u </field>
-	  <field field_ref="utau"         name="tauuo"    />
           <field field_ref="uocetr_eff"   name="uocetr_eff"  />
           <!-- available with diaar5 -->
 	  <field field_ref="u_masstr"     name="vozomatr" />
@@ -56,7 +57,6 @@
 	  <field field_ref="e3v" />
 	  <field field_ref="ssv"          name="vos"      />
 	  <field field_ref="voce"         name="vo"       operation="instant" freq_op="5d" > @voce_e3v / @e3v </field>
-	  <field field_ref="vtau"         name="tauvo"    />
           <field field_ref="vocetr_eff"   name="vocetr_eff"  />
           <!-- available with diaar5 -->
 	  <field field_ref="v_masstr"     name="vomematr" />
diff --git a/cfgs/WED025/EXPREF/file_def_nemo-oce.xml b/cfgs/WED025/EXPREF/file_def_nemo-oce.xml
index 8d9d89e5fb2ddc8a404054b9b15ac9bf4c19cd3e..d08f13e1f8bae0dfaf39b75951205adbd5c0c194 100644
--- a/cfgs/WED025/EXPREF/file_def_nemo-oce.xml
+++ b/cfgs/WED025/EXPREF/file_def_nemo-oce.xml
@@ -31,6 +31,8 @@
 	  <field field_ref="qt_oce"       name="qt_oce"   />
 	  <field field_ref="saltflx"      name="sfx"      />
 	  <field field_ref="taum"         name="taum"     />
+	  <field field_ref="utau"         name="tauuo"    />
+	  <field field_ref="vtau"         name="tauvo"    />
 	  <field field_ref="wspd"         name="windsp"   />
 	  <field field_ref="precip"       name="precip"   />
 	  <field field_ref="snowpre"                      />
@@ -67,14 +69,12 @@
 	  <field field_ref="e3u" />
 	  <field field_ref="ssu"          name="uos"      />
 	  <field field_ref="uoce"         name="uo"       operation="instant" freq_op="5d" > @uoce_e3u / @e3u </field>
-	  <field field_ref="utau"         name="tauuo"    />
 	</file>
 	
 	<file id="file13" name_suffix="_grid_V" description="ocean V grid variables" >
 	  <field field_ref="e3v" />
 	  <field field_ref="ssv"          name="vos"      />
 	  <field field_ref="voce"         name="vo"       operation="instant" freq_op="5d" > @voce_e3v / @e3v </field>
-	  <field field_ref="vtau"         name="tauvo"    />
 	</file>
 	
 	<file id="file14" name_suffix="_grid_W" description="ocean W grid variables" >
diff --git a/sette/BATCH_TEMPLATE/batch-X64_JEANZAY b/sette/BATCH_TEMPLATE/batch-X64_JEANZAY
index e7af8cea409ad356de3c20c0fc777b3fcb767e2d..cafaff218d723068ced9983a689be4ccef41bd4e 100644
--- a/sette/BATCH_TEMPLATE/batch-X64_JEANZAY
+++ b/sette/BATCH_TEMPLATE/batch-X64_JEANZAY
@@ -2,6 +2,7 @@
 #SBATCH -A GROUP_IDRIS@cpu
 #SBATCH --job-name=SETTE_JOB      # nom du job
 #SBATCH --partition=cpu_p1        # Nom de la partition d'exécution
+#SBATCH --qos=qos_cpu-dev           # quality of service
 #SBATCH --ntasks=NPROCS                # Nombre total de processus MPI
 #SBATCH --ntasks-per-node=40       # Nombre de processus MPI par noeud
 # /!\ Attention, la ligne suivante est trompeuse mais dans le vocabulaire
diff --git a/src/ABL/ablmod.F90 b/src/ABL/ablmod.F90
index b5b66b8a6d129be51693c33859efd507f543476d..ffaddbbcd36d7635404b1ef1d294790d4846ea31 100644
--- a/src/ABL/ablmod.F90
+++ b/src/ABL/ablmod.F90
@@ -98,8 +98,8 @@ CONTAINS
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pevp_ice     ! Ce x Du over ice (T-point)
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pwndm_ice    ! ||uwnd - uice||
       REAL(wp) , INTENT(in   ), DIMENSION(:,:  ) ::   pfrac_oce    ! ocean fraction
-      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui_ice    ! ice-surface taux stress (U-point)
-      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj_ice    ! ice-surface tauy stress (V-point)
+      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptaui_ice    ! ice-surface taux stress (T-point)
+      REAL(wp) , INTENT(  out), DIMENSION(:,:  ) ::   ptauj_ice    ! ice-surface tauy stress (T-point)
 #endif
       !
       REAL(wp), DIMENSION(1:jpi,1:jpj       )    ::   zwnd_i, zwnd_j
@@ -617,25 +617,17 @@ CONTAINS
          zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
          zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
       END_2D
-      ! ... utau, vtau at U- and V_points, resp.
-      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
-      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
-      DO_2D( 0, 0, 0, 0 )
-         zcff  = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji+1,jj) )
-         zztmp = MAX(msk_abl(ji,jj),msk_abl(ji+1,jj))
-         ptaui(ji,jj) = zcff * zztmp * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) )
-         zcff  = 0.5_wp * ( 2._wp - msk_abl(ji,jj)*msk_abl(ji,jj+1) )
-         zztmp = MAX(msk_abl(ji,jj),msk_abl(ji,jj+1))
-         ptauj(ji,jj) = zcff * zztmp * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) )
+      ! ... utau, vtau at T-points
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         ptaui(ji,jj) = zwnd_i(ji,jj) * msk_abl(ji,jj) !!clem tau: check
+         ptauj(ji,jj) = zwnd_j(ji,jj) * msk_abl(ji,jj) !!clem tau: check
       END_2D
-      !
-      CALL lbc_lnk( 'ablmod', ptaui(:,:), 'U', -1.0_wp, ptauj(:,:), 'V', -1.0_wp )
 
       CALL iom_put( "taum_oce", ptaum )
 
       IF(sn_cfctl%l_prtctl) THEN
-         CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau   : ', mask1=umask,   &
-            &          tab2d_2=ptauj , clinfo2='          vtau   : ', mask2=vmask )
+         CALL prt_ctl( tab2d_1=ptaui , clinfo1=' abl_stp: utau   : ', mask1=tmask,   &
+            &          tab2d_2=ptauj , clinfo2='          vtau   : ', mask2=tmask )
          CALL prt_ctl( tab2d_1=pwndm , clinfo1=' abl_stp: wndm   : ' )
       ENDIF
 
@@ -657,23 +649,14 @@ CONTAINS
       ! ------------------------------------------------------------ !
       !    Wind stress relative to the moving ice ( U10m - U_ice )   !
       ! ------------------------------------------------------------ !
-      DO_2D( 0, 0, 0, 0 )
-
-         zztmp1 = 0.5_wp * ( u_abl(ji+1,jj  ,2,nt_a) + u_abl(ji,jj,2,nt_a) )
-         zztmp2 = 0.5_wp * ( v_abl(ji  ,jj+1,2,nt_a) + v_abl(ji,jj,2,nt_a) )
-
-         ptaui_ice(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * pCd_du_ice(ji+1,jj)             &
-            &                      +    rhoa(ji  ,jj) * pCd_du_ice(ji  ,jj)  )          &
-            &         * ( zztmp1 - pssu_ice(ji,jj) * rn_vfac )
-         ptauj_ice(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * pCd_du_ice(ji,jj+1)             &
-            &                      +    rhoa(ji,jj  ) * pCd_du_ice(ji,jj  )  )          &
-            &         * ( zztmp2 - pssv_ice(ji,jj) * rn_vfac )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         ptaui_ice(ji,jj) = rhoa(ji,jj) * pCd_du_ice(ji,jj) * ( u_abl(ji,jj,2,nt_a) - pssu_ice(ji,jj) * rn_vfac )
+         ptauj_ice(ji,jj) = rhoa(ji,jj) * pCd_du_ice(ji,jj) * ( v_abl(ji,jj,2,nt_a) - pssv_ice(ji,jj) * rn_vfac )
       END_2D
-      CALL lbc_lnk( 'ablmod', ptaui_ice, 'U', -1.0_wp, ptauj_ice,'V', -1.0_wp )
       !
       IF(sn_cfctl%l_prtctl) THEN
-         CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: utau_ice : ', mask1=umask,   &
-            &          tab2d_2=ptauj_ice , clinfo2='          vtau_ice : ', mask2=vmask )
+         CALL prt_ctl( tab2d_1=ptaui_ice , clinfo1=' abl_stp: utau_ice : ', mask1=tmask,   &
+            &          tab2d_2=ptauj_ice , clinfo2='          vtau_ice : ', mask2=tmask )
       END IF
 #endif
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
diff --git a/src/ABL/sbcabl.F90 b/src/ABL/sbcabl.F90
index f087614e39166824aa5084b19c44237ccc957963..c2c1d13c12941ac034f32e782f553b30bba719c2 100644
--- a/src/ABL/sbcabl.F90
+++ b/src/ABL/sbcabl.F90
@@ -308,8 +308,8 @@ CONTAINS
       !!              - Perform 1 time-step of the ABL model
       !!              - Finalize flux computation in blk_oce_2
       !!
-      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
-      !!              - vtau    : j-component of the stress at V-point  (N/m2)
+      !! ** Outputs : - utau    : i-component of the stress at T-point  (N/m2)
+      !!              - vtau    : j-component of the stress at T-point  (N/m2)
       !!              - taum    : Wind stress module at T-point         (N/m2)
       !!              - wndm    : Wind speed module at T-point          (m/s)
       !!              - qsr     : Solar heat flux over the ocean        (W/m2)
diff --git a/src/ICE/icectl.F90 b/src/ICE/icectl.F90
index ad20e6734e4f96454372ef4f8ea548145b473d4a..9593be4f8f0fd7c38b2cd4f486cf676435fff189 100644
--- a/src/ICE/icectl.F90
+++ b/src/ICE/icectl.F90
@@ -733,10 +733,10 @@ CONTAINS
       CALL prt_ctl_info(' ')
       CALL prt_ctl_info(' - Stresses : ')
       CALL prt_ctl_info('   ~~~~~~~~~~ ')
-      CALL prt_ctl(tab2d_1=utau       , clinfo1= ' utau      : ', mask1 = umask,  &
-         &         tab2d_2=vtau       , clinfo2= ' vtau      : ', mask2 = vmask)
-      CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' utau_ice  : ', mask1 = umask,  &
-         &         tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ', mask2 = vmask)
+      CALL prt_ctl(tab2d_1=utau       , clinfo1= ' utau      : ', mask1 = tmask,  &
+         &         tab2d_2=vtau       , clinfo2= ' vtau      : ', mask2 = tmask)
+      CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' utau_ice  : ', mask1 = tmask,  &
+         &         tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ', mask2 = tmask)
 
    END SUBROUTINE ice_prt3D
 
diff --git a/src/ICE/icedyn_rhg_eap.F90 b/src/ICE/icedyn_rhg_eap.F90
index 2dcfc57595b72ea417470e89f16baeeae274f7f1..797ceba713de4241a41a013adc9add3387d9e756 100644
--- a/src/ICE/icedyn_rhg_eap.F90
+++ b/src/ICE/icedyn_rhg_eap.F90
@@ -311,8 +311,12 @@ CONTAINS
          zmV_t(ji,jj)    = zmassV * z1_dtevp
 
          ! Drag ice-atm.
-         ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
-         ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
+         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
+         !      and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
+         ztaux_ai(ji,jj) = zaU(ji,jj) * 0.5_wp * ( utau_ice(ji,jj) + utau_ice(ji+1,jj) ) * &
+            &                                    ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji+1,jj,1) )
+         ztauy_ai(ji,jj) = zaV(ji,jj) * 0.5_wp * ( vtau_ice(ji,jj) + vtau_ice(ji,jj+1) ) * &
+            &                                    ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji,jj+1,1) )
 
          ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
          zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
diff --git a/src/ICE/icedyn_rhg_evp.F90 b/src/ICE/icedyn_rhg_evp.F90
index 52f15aae8612760424ab298201a0a89508e67d59..a260797a2a630a096e4e389072ab020f676c6a32 100644
--- a/src/ICE/icedyn_rhg_evp.F90
+++ b/src/ICE/icedyn_rhg_evp.F90
@@ -292,8 +292,12 @@ CONTAINS
          zmV_t(ji,jj)    = zmassV * z1_dtevp
 
          ! Drag ice-atm.
-         ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
-         ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
+         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
+         !      and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
+         ztaux_ai(ji,jj) = zaU(ji,jj) * 0.5_wp * ( utau_ice(ji,jj) + utau_ice(ji+1,jj) ) * &
+            &                                    ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji+1,jj,1) )
+         ztauy_ai(ji,jj) = zaV(ji,jj) * 0.5_wp * ( vtau_ice(ji,jj) + vtau_ice(ji,jj+1) ) * &
+            &                                    ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji,jj+1,1) )
 
          ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
          zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
diff --git a/src/ICE/icedyn_rhg_vp.F90 b/src/ICE/icedyn_rhg_vp.F90
index c087f97a8c2cad4cd03c234b646906e34eeaf515..300336ccb889921aa7898f852a37201fce2cf34b 100644
--- a/src/ICE/icedyn_rhg_vp.F90
+++ b/src/ICE/icedyn_rhg_vp.F90
@@ -332,7 +332,10 @@ 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)
          
          ! Wind stress
-         ztaux_ai(ji,jj) = za_iU(ji,jj) * utau_ice(ji,jj)
+         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
+         !      and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
+         ztaux_ai(ji,jj) = za_iU(ji,jj) * 0.5_wp * ( utau_ice(ji,jj) + utau_ice(ji+1,jj) ) * &
+            &                                      ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji+1,jj,1) )
          
          ! Force due to sea surface tilt(- m*g*GRAD(ssh))
          zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
@@ -369,7 +372,10 @@ CONTAINS
          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)
          
          ! Wind stress
-         ztauy_ai(ji,jj) = za_iV(ji,jj) * vtau_ice(ji,jj)
+         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
+         !      and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
+         ztauy_ai(ji,jj) = za_iV(ji,jj) * 0.5_wp * ( vtau_ice(ji,jj) + vtau_ice(ji,jj+1) ) * &
+            &                                      ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji,jj+1,1) )
          
          ! Force due to sea surface tilt(- m*g*GRAD(ssh))
          zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
diff --git a/src/ICE/icesbc.F90 b/src/ICE/icesbc.F90
index 8191540529a5316c0de590c7ab89e60b916d4af5..553e449977ecb918dceb80d381e687c85a001c72 100644
--- a/src/ICE/icesbc.F90
+++ b/src/ICE/icesbc.F90
@@ -272,15 +272,13 @@ CONTAINS
             zvel (ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj  ) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) ) + &
                &                          ( v_ice(ji  ,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) ) )
          END_2D
+         CALL lbc_lnk( 'icesbc', zfric, 'T',  1.0_wp, zvel, 'T', 1.0_wp )
       ELSE      !  if no ice dynamics => transfer directly the atmospheric stress to the ocean
-         DO_2D( 0, 0, 0, 0 )
-            zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp *  &
-               &                         (  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
-               &                          + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1)
-            zvel(ji,jj) = 0._wp
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            zfric(ji,jj) = r1_rho0 * SQRT( utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj) ) * tmask(ji,jj,1)
+            zvel (ji,jj) = 0._wp
          END_2D
       ENDIF
-      CALL lbc_lnk( 'icesbc', zfric, 'T',  1.0_wp, zvel, 'T', 1.0_wp )
       !
       !--------------------------------------------------------------------!
       ! Partial computation of forcing for the thermodynamic sea ice model
diff --git a/src/ICE/icethd_do.F90 b/src/ICE/icethd_do.F90
index 3cd1bcd7ba3875bf6704a267ac63b4fd84f15981..ace1000aa497d30ca52321d59e4a898dcff48e90 100644
--- a/src/ICE/icethd_do.F90
+++ b/src/ICE/icethd_do.F90
@@ -372,8 +372,8 @@ CONTAINS
          DO_2D( 0, 0, 0, 0 )
             IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling
                ! -- Wind stress -- !
-               ztaux = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1) + utau_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
-               ztauy = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1) + vtau_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
+               ztaux = utau_ice(ji,jj) * tmask(ji,jj,1)
+               ztauy = vtau_ice(ji,jj) * tmask(ji,jj,1)
                ! Square root of wind stress
                ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
 
diff --git a/src/ICE/iceupdate.F90 b/src/ICE/iceupdate.F90
index a973312f04364781465403618f175cf52c6170ed..f6d047adf7c6abfbef6897311e8161cce3caa5f5 100644
--- a/src/ICE/iceupdate.F90
+++ b/src/ICE/iceupdate.F90
@@ -314,7 +314,7 @@ CONTAINS
       !!
       !! ** Action  : * at each ice time step (every nn_fsbc time step):
       !!                - compute the modulus of ice-ocean relative velocity
-      !!                  (*rho*Cd) at T-point (C-grid) or I-point (B-grid)
+      !!                  (*rho*Cd) at T-point (C-grid)
       !!                      tmod_io = rhoco * | U_ice-U_oce |
       !!                - update the modulus of stress at ocean surface
       !!                      taum = (1-a) * taum + a * tmod_io * | U_ice-U_oce |
@@ -325,19 +325,19 @@ CONTAINS
       !!
       !!    NB: - ice-ocean rotation angle no more allowed
       !!        - here we make an approximation: taum is only computed every ice time step
-      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
-      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
+      !!          This avoids mutiple average to pass from U,V grids to T grids
+      !!          taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
       !!
-      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
+      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (T-pts) updated with ice-ocean fluxes
       !!              - taum         : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
       !!---------------------------------------------------------------------
       INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
       !
       INTEGER  ::   ji, jj   ! dummy loop indices
-      REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar
-      REAL(wp) ::   zat_v, zvtau_ice, zv_t, zrhoco  !   -      -
-      REAL(wp) ::   zflagi                          !   -      -
+      REAL(wp) ::   zutau_ice, zu_t, zmodt   ! local scalar
+      REAL(wp) ::   zvtau_ice, zv_t, zrhoco  !   -      -
+      REAL(wp) ::   zflagi                   !   -      -
       !!---------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('iceupdate')
 
@@ -352,8 +352,8 @@ CONTAINS
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
          DO_2D( 0, 0, 0, 0 )                          !* update the modulus of stress at ocean surface (T-point)
             !                                               ! 2*(U_ice-U_oce) at T-point
-            zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)
-            zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)
+            zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj) ! u_oce = ssu_m
+            zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) ! v_oce = ssv_m
             !                                              ! |U_ice-U_oce|^2
             zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )
             !                                               ! update the ocean stress modulus
@@ -377,19 +377,14 @@ CONTAINS
       ENDIF
       !
       DO_2D( 0, 0, 0, 0 )                             !* update the stress WITHOUT an ice-ocean rotation angle
-         ! ice area at u and v-points
-         zat_u  = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji+1,jj    ) * tmask(ji+1,jj  ,1) )  &
-            &     / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji+1,jj  ,1) )
-         zat_v  = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji  ,jj+1  ) * tmask(ji  ,jj+1,1) )  &
-            &     / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji  ,jj+1,1) )
          !                                                   ! linearized quadratic drag formulation
-         zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
-         zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
+         zutau_ice   = 0.5_wp * tmod_io(ji,jj) * ( u_ice(ji,jj) + u_ice(ji-1,jj) - pu_oce(ji,jj) - pu_oce(ji-1,jj) )
+         zvtau_ice   = 0.5_wp * tmod_io(ji,jj) * ( v_ice(ji,jj) + v_ice(ji,jj-1) - pv_oce(ji,jj) - pv_oce(ji,jj-1) )
          !                                                   ! stresses at the ocean surface
-         utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
-         vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
+         utau(ji,jj) = ( 1._wp - at_i(ji,jj) ) * utau_oce(ji,jj) + at_i(ji,jj) * zutau_ice
+         vtau(ji,jj) = ( 1._wp - at_i(ji,jj) ) * vtau_oce(ji,jj) + at_i(ji,jj) * zvtau_ice
       END_2D
-      CALL lbc_lnk( 'iceupdate', utau, 'U', -1.0_wp, vtau, 'V', -1.0_wp )   ! lateral boundary condition
+      CALL lbc_lnk( 'iceupdate', utau, 'T', -1.0_wp, vtau, 'T', -1.0_wp )   ! lateral boundary condition
       !
       IF( ln_timing )   CALL timing_stop('iceupdate')
       !
diff --git a/src/OCE/CRS/crsfld.F90 b/src/OCE/CRS/crsfld.F90
index b8115048d68e1faff79a13d65c186ca77fb2a250..efa11cf50690eaa4458ab7d57cad81ec47fab453 100644
--- a/src/OCE/CRS/crsfld.F90
+++ b/src/OCE/CRS/crsfld.F90
@@ -211,8 +211,8 @@ CONTAINS
       
       !  sbc fields  
       CALL crs_dom_ope( ssh(:,:,Kmm) , 'VOL', 'T', tmask, sshn_crs , p_e12=e1e2t, p_e3=ze3t           , psgn=1.0_wp )  
-      CALL crs_dom_ope( utau , 'SUM', 'U', umask, utau_crs , p_e12=e2u  , p_surf_crs=e2u_crs  , psgn=1.0_wp )
-      CALL crs_dom_ope( vtau , 'SUM', 'V', vmask, vtau_crs , p_e12=e1v  , p_surf_crs=e1v_crs  , psgn=1.0_wp )
+      CALL crs_dom_ope( utau , 'SUM', 'T', tmask, utau_crs , p_e12=e2t  , p_surf_crs=e2t_crs  , psgn=1.0_wp ) !clem tau: check psgn ??
+      CALL crs_dom_ope( vtau , 'SUM', 'T', tmask, vtau_crs , p_e12=e1t  , p_surf_crs=e1t_crs  , psgn=1.0_wp ) !
       CALL crs_dom_ope( wndm , 'SUM', 'T', tmask, wndm_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0_wp )
       CALL crs_dom_ope( rnf  , 'MAX', 'T', tmask, rnf_crs                                     , psgn=1.0_wp )
       CALL crs_dom_ope( qsr  , 'SUM', 'T', tmask, qsr_crs  , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0_wp )
diff --git a/src/OCE/DIA/diawri.F90 b/src/OCE/DIA/diawri.F90
index db5da93a7cc7f9a1d61498a55c689bc90aa7c231..0deeb3370630962fda61e471ed735bda60f26807 100644
--- a/src/OCE/DIA/diawri.F90
+++ b/src/OCE/DIA/diawri.F90
@@ -868,7 +868,11 @@ CONTAINS
          CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "J/m2"   ,   & ! htc3
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
 #endif
-
+         CALL histdef( nid_T, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   & ! utau
+            &          jpi, jpj, nh_T, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
+         CALL histdef( nid_T, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   & ! vtau
+            &          jpi, jpj, nh_T, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
+         !
          CALL histend( nid_T, snc4chunks=snc4set )
 
          !                                                                                      !!! nid_U : 3D
@@ -878,10 +882,7 @@ CONTAINS
             CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd
                &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
          ENDIF
-         !                                                                                      !!! nid_U : 2D
-         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
-            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
-
+         !
          CALL histend( nid_U, snc4chunks=snc4set )
 
          !                                                                                      !!! nid_V : 3D
@@ -891,10 +892,7 @@ CONTAINS
             CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd
                &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
          ENDIF
-         !                                                                                      !!! nid_V : 2D
-         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
-            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
-
+         !
          CALL histend( nid_V, snc4chunks=snc4set )
 
          !                                                                                      !!! nid_W : 3D
@@ -1066,12 +1064,12 @@ CONTAINS
       CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
       CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
 #endif
+      CALL histwrite( nid_T, "sozotaux", it, utau          , ndim_hT, ndex_hT )   ! i-wind stress
+      CALL histwrite( nid_T, "sometauy", it, vtau          , ndim_hT, ndex_hT )   ! j-wind stress
 
       CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm) , ndim_U , ndex_U )    ! i-current
-      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
 
       CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm) , ndim_V , ndex_V  )   ! j-current
-      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
 
       IF( ln_zad_Aimp ) THEN
          DO_3D( 0, 0, 0, 0, 1, jpk )
diff --git a/src/OCE/DYN/dynatf.F90 b/src/OCE/DYN/dynatf.F90
index 997692694ddf1bf7dd1fa04a5c50e96759b4e3c2..1a5df99676dabf0493e33335dc2564cc3ce6a8fb 100644
--- a/src/OCE/DYN/dynatf.F90
+++ b/src/OCE/DYN/dynatf.F90
@@ -331,7 +331,7 @@ CONTAINS
             ALLOCATE(zutau(jpi,jpj))
             DO_2D( 0, 0, 0, 0 )
                jk = miku(ji,jj)
-               zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa)
+               zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( puu(ji-1,jj,jk,Kaa) + puu(ji,jj,jk,Kaa) )
             END_2D
             CALL iom_put(  "utau", zutau(:,:) )
             DEALLOCATE(zutau)
@@ -345,7 +345,7 @@ CONTAINS
             ALLOCATE(zvtau(jpi,jpj))
             DO_2D( 0, 0, 0, 0 )
                jk = mikv(ji,jj)
-               zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa)
+               zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( pvv(ji,jj-1,jk,Kaa) + pvv(ji,jj,jk,Kaa) )
             END_2D
             CALL iom_put(  "vtau", zvtau(:,:) )
             DEALLOCATE(zvtau)
diff --git a/src/OCE/DYN/dynatf_qco.F90 b/src/OCE/DYN/dynatf_qco.F90
index f9ed9f4644d525d233b884f53ad0e7f37760253c..3dd46032e952fc59270a5047155ac46c1570516a 100644
--- a/src/OCE/DYN/dynatf_qco.F90
+++ b/src/OCE/DYN/dynatf_qco.F90
@@ -248,7 +248,7 @@ CONTAINS
             ALLOCATE(zutau(jpi,jpj))
             DO_2D( 0, 0, 0, 0 )
                jk = miku(ji,jj)
-               zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa)
+               zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( puu(ji-1,jj,jk,Kaa) + puu(ji,jj,jk,Kaa) )
             END_2D
             CALL iom_put(  "utau", zutau(:,:) )
             DEALLOCATE(zutau)
@@ -262,7 +262,7 @@ CONTAINS
             ALLOCATE(zvtau(jpi,jpj))
             DO_2D( 0, 0, 0, 0 )
                jk = mikv(ji,jj)
-               zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa)
+               zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( pvv(ji,jj-1,jk,Kaa) + pvv(ji,jj,jk,Kaa) )
             END_2D
             CALL iom_put(  "vtau", zvtau(:,:) )
             DEALLOCATE(zvtau)
diff --git a/src/OCE/DYN/dynspg_ts.F90 b/src/OCE/DYN/dynspg_ts.F90
index 83b2d99cac15b0fec655681371eb81ca929b9064..16c2fcc1fa8e3244dc5b5c4a5741340ce0acf0fc 100644
--- a/src/OCE/DYN/dynspg_ts.F90
+++ b/src/OCE/DYN/dynspg_ts.F90
@@ -334,14 +334,14 @@ CONTAINS
       !                                   !  ------------------  !
       IF( ln_bt_fw ) THEN
          DO_2D( 0, 0, 0, 0 )
-            zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)
-            zv_frc(ji,jj) =  zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm)
+            zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rho0 * utauU(ji,jj) * r1_hu(ji,jj,Kmm)
+            zv_frc(ji,jj) =  zv_frc(ji,jj) + r1_rho0 * vtauV(ji,jj) * r1_hv(ji,jj,Kmm)
          END_2D
       ELSE
          zztmp = r1_rho0 * r1_2
          DO_2D( 0, 0, 0, 0 )
-            zu_frc(ji,jj) =  zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm)
-            zv_frc(ji,jj) =  zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm)
+            zu_frc(ji,jj) =  zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utauU(ji,jj) ) * r1_hu(ji,jj,Kmm)
+            zv_frc(ji,jj) =  zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtauV(ji,jj) ) * r1_hv(ji,jj,Kmm)
          END_2D
       ENDIF  
       !
diff --git a/src/OCE/DYN/dynzdf.F90 b/src/OCE/DYN/dynzdf.F90
index 1e026021941c8a5891459385a52d899c8c88123e..98c65e4604e951f0119fbd4767381487ec3596a0 100644
--- a/src/OCE/DYN/dynzdf.F90
+++ b/src/OCE/DYN/dynzdf.F90
@@ -267,10 +267,10 @@ CONTAINS
       DO_2D( 0, 0, 0, 0 )             !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
 #if defined key_RK3
          !                                  ! RK3: use only utau (not utau_b)
-         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utau(ji,jj)   &
+         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utauU(ji,jj)   &
               &                                 / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
 #else
-         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utau(ji,jj) )   &
+         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utauU(ji,jj) )   &
               &                                 / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
 #endif
       END_2D
@@ -397,10 +397,10 @@ CONTAINS
       DO_2D( 0, 0, 0, 0 )             !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
 #if defined key_RK3
          !                                  ! RK3: use only vtau (not vtau_b)
-         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtau(ji,jj)   &
+         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtauV(ji,jj)   &
             &                                   / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
 #else
-         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2*( vtau_b(ji,jj) + vtau(ji,jj) )   &
+         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2 * ( vtau_b(ji,jj) + vtauV(ji,jj) )   &
               &                                 / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
 #endif
       END_2D
diff --git a/src/OCE/SBC/cpl_oasis3.F90 b/src/OCE/SBC/cpl_oasis3.F90
index 0147ffdac802e5234499df6ae87b268afefcab8b..3c7e15ba498c8dfea73eaf4365356bd7b82089ce 100644
--- a/src/OCE/SBC/cpl_oasis3.F90
+++ b/src/OCE/SBC/cpl_oasis3.F90
@@ -419,6 +419,7 @@ CONTAINS
          IF( .NOT. ll_1st ) THEN
             CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,jc), srcv(kid)%clgrid, srcv(kid)%nsgn )
          ENDIF
+         !!clem: mettre T instead of clgrid
 
       ENDDO
       !
diff --git a/src/OCE/SBC/sbc_ice.F90 b/src/OCE/SBC/sbc_ice.F90
index acc6adf1fd10939f9a8e0b742f951b1330fbe844..cac087f79246701a2fb0197b1b7e3bf934d87a30 100644
--- a/src/OCE/SBC/sbc_ice.F90
+++ b/src/OCE/SBC/sbc_ice.F90
@@ -50,8 +50,8 @@ MODULE sbc_ice
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice        !: heat conduction flux in the layer below surface   [W/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qtr_ice_top    !: solar flux transmitted below the ice surface      [W/m2]
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   utau_ice       !: atmos-ice u-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vtau_ice       !: atmos-ice v-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   utau_ice       !: atmos-ice u-stress. T-pts                  [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vtau_ice       !: atmos-ice v-stress. T-pts                  [N/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice        !: sublimation - precip over sea ice          [kg/m2/s]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   topmelt            !: category topmelt
diff --git a/src/OCE/SBC/sbc_oce.F90 b/src/OCE/SBC/sbc_oce.F90
index 3098b5865b6719eb5250180257f7d9e1ab143c75..f4cd6c97ee11a04d1aa9af9db8dba2b9242fd774 100644
--- a/src/OCE/SBC/sbc_oce.F90
+++ b/src/OCE/SBC/sbc_oce.F90
@@ -103,34 +103,36 @@ MODULE sbc_oce
    INTEGER , PUBLIC ::  ncpl_qsr_freq = 0        !: qsr coupling frequency per days from atmosphere (used by top)
    !
    !!                                   !!   now    ! before   !!
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau   , utau_b   !: sea surface i-stress (ocean referential)     [N/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vtau   , vtau_b   !: sea surface j-stress (ocean referential)     [N/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_icb, vtau_icb !: sea surface (i,j)-stress used by icebergs   [N/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   taum              !: module of sea surface stress (at T-point)    [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau               !: sea surface i-stress (ocean referential) T-pt [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vtau               !: sea surface j-stress (ocean referential) T-pt [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utauU  , utau_b    !: sea surface i-stress (ocean referential) U-pt [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vtauV  , vtau_b    !: sea surface j-stress (ocean referential) V-pt [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_icb, vtau_icb !: sea surface (i,j)-stress used by icebergs     [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   taum               !: module of sea surface stress (at T-point)     [N/m2]
    !! wndm is used compute surface gases exchanges in ice-free ocean or leads
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rhoa              !: air density at "rn_zu" m above the sea       [kg/m3]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qsr               !: sea heat flux:     solar                     [W/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qns    , qns_b    !: sea heat flux: non solar                     [W/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qsr_tot           !: total     solar heat flux (over sea and ice) [W/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qns_tot           !: total non solar heat flux (over sea and ice) [W/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emp    , emp_b    !: freshwater budget: volume flux               [Kg/m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx    , sfx_b    !: salt flux                                    [PSS.kg/m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emp_tot           !: total E-P over ocean and ice                 [Kg/m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmflx            !: freshwater budget: freezing/melting          [Kg/m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff                                 [Kg/m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwficb , fwficb_b !: iceberg melting                              [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wndm               !: wind speed module at T-point (=|U10m-Uoce|)   [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rhoa               !: air density at "rn_zu" m above the sea        [kg/m3]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qsr                !: sea heat flux:     solar                      [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qns    , qns_b     !: sea heat flux: non solar                      [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qsr_tot            !: total     solar heat flux (over sea and ice)  [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qns_tot            !: total non solar heat flux (over sea and ice)  [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emp    , emp_b     !: freshwater budget: volume flux                [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx    , sfx_b     !: salt flux                                     [PSS.kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emp_tot            !: total E-P over ocean and ice                  [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmflx             !: freshwater budget: freezing/melting           [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b     !: river runoff                                  [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwficb , fwficb_b  !: iceberg melting                               [Kg/m2/s]
    !!
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  sbc_tsc, sbc_tsc_b  !: sbc content trend                      [K.m/s] jpi,jpj,jpts
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  qsr_hc , qsr_hc_b   !: heat content trend due to qsr flux     [K.m/s] jpi,jpj,jpk
    !!
    !!
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tprecip           !: total precipitation                          [Kg/m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sprecip           !: solid precipitation                          [Kg/m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fr_i              !: ice fraction = 1 - lead fraction      (between 0 to 1)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cloud_fra         !: cloud cover (fraction of cloud in a gridcell) [-]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tprecip            !: total precipitation                           [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sprecip            !: solid precipitation                           [Kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fr_i               !: ice fraction = 1 - lead fraction       (between 0 to 1)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2            !: atmospheric pCO2                              [ppm]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask           !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cloud_fra          !: cloud cover (fraction of cloud in a gridcell) [-]
 
    !!---------------------------------------------------------------------
    !! ABL Vertical Domain size
@@ -177,8 +179,8 @@ CONTAINS
       !!---------------------------------------------------------------------
       ierr(:) = 0
       !
-      ALLOCATE( utau(jpi,jpj) , utau_b(jpi,jpj) , taum(jpi,jpj) ,     &
-         &      vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , rhoa(jpi,jpj) , STAT=ierr(1) )
+      ALLOCATE( utau(jpi,jpj) , utau_b(jpi,jpj) , utauU(jpi,jpj) , taum(jpi,jpj) ,     &
+         &      vtau(jpi,jpj) , vtau_b(jpi,jpj) , vtauV(jpi,jpj) , wndm(jpi,jpj) , rhoa(jpi,jpj) , STAT=ierr(1) )
       !
       ALLOCATE( qns_tot(jpi,jpj) , qns  (jpi,jpj) , qns_b(jpi,jpj),        &
          &      qsr_tot(jpi,jpj) , qsr  (jpi,jpj) ,                        &
@@ -205,9 +207,10 @@ CONTAINS
    END FUNCTION sbc_oce_alloc
 
 
+   !!clem => this subroutine is never used in nemo
    SUBROUTINE sbc_tau2wnd
       !!---------------------------------------------------------------------
-      !!                    ***  ROUTINE sbc_tau2wnd  ***
+      !!                    ***  ROUTINE   ***
       !!
       !! ** Purpose : Estimation of wind speed as a function of wind stress
       !!
@@ -217,17 +220,14 @@ CONTAINS
       USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
       REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
       REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
-      REAL(wp) ::   ztx, zty, ztau, zcoef ! temporary variables
+      REAL(wp) ::   ztau, zcoef           ! temporary variables
       INTEGER  ::   ji, jj                ! dummy indices
       !!---------------------------------------------------------------------
       zcoef = 0.5 / ( zrhoa * zcdrag )
-      DO_2D( 0, 0, 0, 0 )
-         ztx = utau(ji-1,jj  ) + utau(ji,jj)
-         zty = vtau(ji  ,jj-1) + vtau(ji,jj)
-         ztau = SQRT( ztx * ztx + zty * zty )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         ztau = SQRT( utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj) )
          wndm(ji,jj) = SQRT ( ztau * zcoef ) * tmask(ji,jj,1)
       END_2D
-      CALL lbc_lnk( 'sbc_oce', wndm(:,:) , 'T', 1.0_wp )
       !
    END SUBROUTINE sbc_tau2wnd
 
diff --git a/src/OCE/SBC/sbcblk.F90 b/src/OCE/SBC/sbcblk.F90
index aa3668c3300bb8d3a46998a76f2f56014332d540..397155a58091501fd54d1da3fedcb6bf39d23797 100644
--- a/src/OCE/SBC/sbcblk.F90
+++ b/src/OCE/SBC/sbcblk.F90
@@ -493,7 +493,7 @@ CONTAINS
       !!                      the stress is assumed to be in the (i,j) mesh referential
       !!
       !! ** Action  :   defined at each time-step at the air-sea interface
-      !!              - utau, vtau  i- and j-component of the wind stress
+      !!              - utau, vtau  i- and j-component of the wind stress at T-point
       !!              - taum        wind stress module at T-point
       !!              - wndm        wind speed  module at T-point over free ocean or leads in presence of sea-ice
       !!              - qns, qsr    non-solar and solar heat fluxes
@@ -611,10 +611,10 @@ CONTAINS
    END SUBROUTINE sbc_blk
 
 
-   SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, pqair,         &  ! inp
-      &                      pslp , pst  , pu   , pv,            &  ! inp
-      &                      puatm, pvatm, pdqsr , pdqlw ,       &  ! inp
-      &                      ptsk , pssq , pcd_du, psen, plat, pevp ) ! out
+   SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, pqair,         &    ! <<= in
+      &                      pslp , pst  , pu   , pv,            &    ! <<= in
+      &                      puatm, pvatm, pdqsr , pdqlw ,       &    ! <<= in
+      &                      ptsk , pssq , pcd_du, psen, plat, pevp ) ! =>> out
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE blk_oce_1  ***
       !!
@@ -657,7 +657,6 @@ CONTAINS
 #if defined key_cyclone
       REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
 #endif
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztau_i, ztau_j    ! wind stress components at T-point
       REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s]
       REAL(wp), DIMENSION(jpi,jpj) ::   zcd_oce           ! momentum transfert coefficient over ocean
       REAL(wp), DIMENSION(jpi,jpj) ::   zch_oce           ! sensible heat transfert coefficient over ocean
@@ -695,22 +694,20 @@ CONTAINS
 #else
       ! ... scalar wind module at T-point (not masked)
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-         wndm(ji,jj) = SQRT(  pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj)  )
+         wndm(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
       END_2D
 #endif
       ! ----------------------------------------------------------------------------- !
       !      I   Solar FLUX                                                           !
       ! ----------------------------------------------------------------------------- !
 
-      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
-      zztmp = 1. - albo
+      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle
       IF( ln_dm2dc ) THEN
-         qsr(:,:) = zztmp * sbc_dcy( pdqsr(:,:) ) * tmask(:,:,1)
+         qsr(:,:) = ( 1._wp - albo ) * sbc_dcy( pdqsr(:,:) ) * tmask(:,:,1)
       ELSE
-         qsr(:,:) = zztmp *          pdqsr(:,:)   * tmask(:,:,1)
+         qsr(:,:) = ( 1._wp - albo ) *          pdqsr(:,:)   * tmask(:,:,1)
       ENDIF
 
-
       ! ----------------------------------------------------------------------------- !
       !     II   Turbulent FLUXES                                                     !
       ! ----------------------------------------------------------------------------- !
@@ -718,69 +715,62 @@ CONTAINS
       ! specific humidity at SST
       pssq(:,:) = rdct_qsat_salt * q_sat( ptsk(:,:), pslp(:,:) )
 
+      ! Backup "bulk SST" and associated spec. hum.
       IF( ln_skin_cs .OR. ln_skin_wl ) THEN
-         !! Backup "bulk SST" and associated spec. hum.
          zztmp1(:,:) = zsspt(:,:)
-         zztmp2(:,:) = pssq(:,:)
+         zztmp2(:,:) = pssq (:,:)
       ENDIF
 
-      !! Time to call the user-selected bulk parameterization for
-      !!  ==  transfer coefficients  ==!   Cd, Ch, Ce at T-point, and more...
-      SELECT CASE( nblk )
-
+      ! transfer coefficients (Cd, Ch, Ce at T-point, and more)
+      SELECT CASE( nblk )   ! user-selected bulk parameterization
+         !
       CASE( np_NCAR      )
          CALL turb_ncar    (     rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
-            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
+            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu ,  &
             &                nb_iter=nn_iter_algo )
-         !
       CASE( np_COARE_3p0 )
          CALL turb_coare3p0( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
-            &                ln_skin_cs, ln_skin_wl,                            &
-            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu,  &
-            &                nb_iter=nn_iter_algo,                              &
+            &                ln_skin_cs, ln_skin_wl,                             &
+            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu,   &
+            &                nb_iter=nn_iter_algo,                               &
             &                Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
-         !
       CASE( np_COARE_3p6 )
          CALL turb_coare3p6( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
-            &                ln_skin_cs, ln_skin_wl,                            &
-            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu,  &
-            &                nb_iter=nn_iter_algo,                              &
+            &                ln_skin_cs, ln_skin_wl,                             &
+            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu,   &
+            &                nb_iter=nn_iter_algo,                               &
             &                Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
-         !
       CASE( np_ECMWF     )
          CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
-            &                ln_skin_cs, ln_skin_wl,                            &
-            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu,  &
-            &                nb_iter=nn_iter_algo,                              &
+            &                ln_skin_cs, ln_skin_wl,                             &
+            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu,   &
+            &                nb_iter=nn_iter_algo,                               &
             &                Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) )
-         !
       CASE( np_ANDREAS   )
          CALL turb_andreas (     rn_zqt, rn_zu, zsspt, ptair, pssq, pqair, wndm, &
-            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , &
+            &                zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu,   &
             &                nb_iter=nn_iter_algo   )
-         !
       CASE DEFAULT
          CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk parameterizaton selected' )
-         !
       END SELECT
 
-      IF( iom_use('Cd_oce') )   CALL iom_put("Cd_oce",   zcd_oce * tmask(:,:,1))
-      IF( iom_use('Ce_oce') )   CALL iom_put("Ce_oce",   zce_oce * tmask(:,:,1))
-      IF( iom_use('Ch_oce') )   CALL iom_put("Ch_oce",   zch_oce * tmask(:,:,1))
+      ! outputs
+      IF( iom_use('Cd_oce') )   CALL iom_put( "Cd_oce",   zcd_oce * tmask(:,:,1) )
+      IF( iom_use('Ce_oce') )   CALL iom_put( "Ce_oce",   zce_oce * tmask(:,:,1) )
+      IF( iom_use('Ch_oce') )   CALL iom_put( "Ch_oce",   zch_oce * tmask(:,:,1) )
       !! LB: mainly here for debugging purpose:
-      IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ptair-rt0) * tmask(:,:,1)) ! potential temperature at z=zt
-      IF( iom_use('q_zt') )     CALL iom_put("q_zt",     pqair       * tmask(:,:,1)) ! specific humidity       "
-      IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (theta_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu
-      IF( iom_use('q_zu') )     CALL iom_put("q_zu",     q_zu        * tmask(:,:,1)) ! specific humidity       "
-      IF( iom_use('ssq') )      CALL iom_put("ssq",      pssq        * tmask(:,:,1)) ! saturation specific humidity at z=0
-      IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu       * tmask(:,:,1)) ! bulk wind speed at z=zu
-
+      IF( iom_use('theta_zt') ) CALL iom_put( "theta_zt", (ptair-rt0) * tmask(:,:,1) ) ! potential temperature at z=zt
+      IF( iom_use('q_zt') )     CALL iom_put( "q_zt",     pqair       * tmask(:,:,1) ) ! specific humidity       "
+      IF( iom_use('theta_zu') ) CALL iom_put( "theta_zu", (theta_zu -rt0) * tmask(:,:,1) ) ! potential temperature at z=zu
+      IF( iom_use('q_zu') )     CALL iom_put( "q_zu",     q_zu        * tmask(:,:,1) ) ! specific humidity       "
+      IF( iom_use('ssq') )      CALL iom_put( "ssq",      pssq        * tmask(:,:,1) ) ! saturation specific humidity at z=0
+      IF( iom_use('wspd_blk') ) CALL iom_put( "wspd_blk", zU_zu       * tmask(:,:,1) ) ! bulk wind speed at z=zu
+
+      ! In the presence of sea-ice we do not use the cool-skin/warm-layer update of zsspt, pssq & ptsk from turb_*()
       IF( ln_skin_cs .OR. ln_skin_wl ) THEN
-         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zsspt, pssq & ptsk:
          WHERE ( fr_i(:,:) > 0.001_wp )
-            ! sea-ice present, we forget about the update, using what we backed up before call to turb_*()
             zsspt(:,:) = zztmp1(:,:)
-            pssq(:,:)  = zztmp2(:,:)
+            pssq (:,:) = zztmp2(:,:)
          END WHERE
          ! apply potential temperature increment to abolute SST
          ptsk(:,:) = ptsk(:,:) + ( zsspt(:,:) - zztmp1(:,:) )
@@ -809,10 +799,10 @@ CONTAINS
          END_2D
 
          CALL BULK_FORMULA( rn_zu, zsspt(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), &
-            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),          &
-            &               wndm(:,:), zU_zu(:,:), pslp(:,:), rhoa(:,:),       &
-            &               taum(:,:), psen(:,:), plat(:,:),                   &
-            &               pEvap=pevp(:,:), pfact_evap=rn_efac )
+            &                      zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),        &
+            &                      wndm(:,:), zU_zu(:,:), pslp(:,:), rhoa(:,:),     &
+            &                      taum(:,:), psen(:,:), plat(:,:),                 &
+            &                      pEvap=pevp(:,:), pfact_evap=rn_efac )
 
          psen(:,:) = psen(:,:) * tmask(:,:,1)
          plat(:,:) = plat(:,:) * tmask(:,:,1)
@@ -821,57 +811,42 @@ CONTAINS
 
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             IF( wndm(ji,jj) > 0._wp ) THEN
-              zztmp = taum(ji,jj) / wndm(ji,jj)
+               zztmp = taum(ji,jj) / wndm(ji,jj)
 #if defined key_cyclone
-               ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
-               ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
+               utau(ji,jj) = zztmp * zwnd_i(ji,jj)
+               vtau(ji,jj) = zztmp * zwnd_j(ji,jj)
 #else
-               ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
-               ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
+               utau(ji,jj) = zztmp * pwndi(ji,jj)
+               vtau(ji,jj) = zztmp * pwndj(ji,jj)
 #endif
             ELSE
-               ztau_i(ji,jj) = 0._wp
-               ztau_j(ji,jj) = 0._wp
+               utau(ji,jj) = 0._wp
+               vtau(ji,jj) = 0._wp
             ENDIF
          END_2D
 
          IF( ln_crt_fbk ) THEN   ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715)
             zstmax = MIN( rn_stau_a * 3._wp + rn_stau_b, 0._wp )   ! set the max value of Stau corresponding to a wind of 3 m/s (<0)
-            DO_2D( 0, 1, 0, 1 )   ! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop
-               zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax )   ! stau (<0) must be smaller than zstmax
-               ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj  ) + pu(ji,jj) ) - puatm(ji,jj) )
-               ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji  ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) )
-               taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) )
+            DO_2D( 0, 0, 0, 0 )
+               zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax ) * tmask(ji,jj,1)   ! stau (<0) must be smaller than zstmax
+               utau(ji,jj) = utau(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj  ) + pu(ji,jj) ) - puatm(ji,jj) )
+               vtau(ji,jj) = vtau(ji,jj) + zstau * ( 0.5_wp * ( pv(ji  ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) )
+               taum(ji,jj) = SQRT( utau(ji,jj) * utau(ji,jj) + vtau(ji,jj) * vtau(ji,jj) )
             END_2D
+            CALL lbc_lnk( 'sbcblk', utau, 'T', -1._wp, vtau, 'T', -1._wp, taum, 'T', 1._wp )
          ENDIF
 
-         ! ... utau, vtau at U- and V_points, resp.
-         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
-         !     Note that coastal wind stress is not used in the code... so this extra care has no effect
-         DO_2D( 0, 0, 0, 0 )              ! start loop at 2, in case ln_crt_fbk = T
-            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj  ) ) &
-               &              * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
-            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji  ,jj+1) ) &
-               &              * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
-         END_2D
-
-         IF( ln_crt_fbk ) THEN
-            CALL lbc_lnk( 'sbcblk', utau, 'U', -1._wp, vtau, 'V', -1._wp, taum, 'T', 1._wp )
-         ELSE
-            CALL lbc_lnk( 'sbcblk', utau, 'U', -1._wp, vtau, 'V', -1._wp )
-         ENDIF
-
-         ! Saving open-ocean wind-stress (module and components) on T-points:
-         CALL iom_put( "taum_oce",   taum(:,:)*tmask(:,:,1) )   ! output wind stress module
-         !#LB: These 2 lines below mostly here for 'STATION_ASF' test-case, otherwize "utau" (U-grid) and vtau" (V-grid) does the job in: [DYN/dynatf.F90])
-         CALL iom_put( "utau_oce", ztau_i(:,:)*tmask(:,:,1) )  ! utau at T-points!
-         CALL iom_put( "vtau_oce", ztau_j(:,:)*tmask(:,:,1) )  ! vtau at T-points!
+         ! Saving open-ocean wind-stress (module and components)
+         CALL iom_put( "taum_oce", taum(:,:) )   ! wind stress module
+         !                                       ! LB: These 2 lines below mostly here for 'STATION_ASF' test-case
+         CALL iom_put( "utau_oce", utau(:,:) )   ! utau
+         CALL iom_put( "vtau_oce", vtau(:,:) )   ! vtau
 
          IF(sn_cfctl%l_prtctl) THEN
             CALL prt_ctl( tab2d_1=pssq   , clinfo1=' blk_oce_1: pssq   : ', mask1=tmask )
             CALL prt_ctl( tab2d_1=wndm   , clinfo1=' blk_oce_1: wndm   : ', mask1=tmask )
-            CALL prt_ctl( tab2d_1=utau   , clinfo1=' blk_oce_1: utau   : ', mask1=umask,   &
-               &          tab2d_2=vtau   , clinfo2='            vtau   : ', mask2=vmask )
+            CALL prt_ctl( tab2d_1=utau   , clinfo1=' blk_oce_1: utau   : ', mask1=tmask,   &
+               &          tab2d_2=vtau   , clinfo2='            vtau   : ', mask2=tmask )
             CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd     : ', mask1=tmask )
          ENDIF
          !
@@ -896,8 +871,8 @@ CONTAINS
       !!                at the ocean surface at each time step knowing Cd, Ch, Ce and
       !!                atmospheric variables (from ABL or external data)
       !!
-      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
-      !!              - vtau    : j-component of the stress at V-point  (N/m2)
+      !! ** Outputs : - utau    : i-component of the stress at T-point  (N/m2)
+      !!              - vtau    : j-component of the stress at T-point  (N/m2)
       !!              - taum    : Wind stress module at T-point         (N/m2)
       !!              - wndm    : Wind speed module at T-point          (m/s)
       !!              - qsr     : Solar heat flux over the ocean        (W/m2)
@@ -1015,11 +990,16 @@ CONTAINS
       REAL(wp) , INTENT(  out), DIMENSION(:,:  ), OPTIONAL ::   pcd_dui ! if ln_abl
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) ::   zootm_su                      ! sea-ice surface mean temperature
-      REAL(wp) ::   zztmp1, zztmp2                ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zsipt ! temporary array
+      REAL(wp) ::   zztmp                           ! temporary scalars
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp, zsipt ! temporary array
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00      ! O% concentration ice mask
       !!---------------------------------------------------------------------
       !
+      ! treshold for outputs
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , fr_i(ji,jj) - 1.e-6_wp  ) ) ! 1 if ice, 0 if no ice
+      END_2D
+
       ! ------------------------------------------------------------ !
       !    Wind module relative to the moving ice ( U10m - U_ice )   !
       ! ------------------------------------------------------------ !
@@ -1032,9 +1012,9 @@ CONTAINS
       zsipt(:,:) = theta_exner( ptsui(:,:), pslp(:,:) )
 
       ! sea-ice <-> atmosphere bulk transfer coefficients
-      SELECT CASE( nblk_ice )
-
-      CASE( np_ice_cst      )
+      SELECT CASE( nblk_ice )   ! user-selected bulk parameterization
+         !
+      CASE( np_ice_cst  )
          ! Constant bulk transfer coefficients over sea-ice:
          Cd_ice(:,:) = rn_Cd_i
          Ch_ice(:,:) = rn_Ch_i
@@ -1042,62 +1022,45 @@ CONTAINS
          ! no height adjustment, keeping zt values:
          theta_zu_i(:,:) = ptair(:,:)
          q_zu_i(:,:)     = pqair(:,:)
-
-      CASE( np_ice_an05 )  ! calculate new drag from Lupkes(2015) equations
+         !
+      CASE( np_ice_an05 )  ! from Andreas(2005) equations
          ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
          CALL turb_ice_an05( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice,       &
             &                      Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
-         !!
-      CASE( np_ice_lu12 )
+         !
+      CASE( np_ice_lu12 )  ! from Lupkes(2012) equations
          ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
          CALL turb_ice_lu12( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, fr_i, &
             &                      Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
-         !!
-      CASE( np_ice_lg15 )  ! calculate new drag from Lupkes(2015) equations
+         !
+      CASE( np_ice_lg15 )  ! from Lupkes and Gryanik (2015) equations
          ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ
          CALL turb_ice_lg15( rn_zqt, rn_zu, zsipt, ptair, ztmp, pqair, wndm_ice, fr_i, &
             &                      Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i )
-         !!
+         !
       END SELECT
 
-      IF( iom_use('Cd_ice').OR.iom_use('Ce_ice').OR.iom_use('Ch_ice').OR.iom_use('taum_ice').OR.iom_use('utau_ice').OR.iom_use('vtau_ice') ) &
-         & ztmp(:,:) = ( 1._wp - MAX(0._wp, SIGN( 1._wp, 1.E-6_wp - fr_i )) )*tmask(:,:,1) ! mask for presence of ice !
-
-      IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice*ztmp)
-      IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice*ztmp)
-      IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice*ztmp)
-
 
       IF( ln_blk ) THEN
          ! ---------------------------------------------------- !
          !    Wind stress relative to nonmoving ice ( U10m )    !
          ! ---------------------------------------------------- !
          ! supress moving ice in wind stress computation as we don't know how to do it properly...
-         DO_2D( 0, 1, 0, 1 )    ! at T point
-            zztmp1        = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj)
-            putaui(ji,jj) = zztmp1 * pwndi(ji,jj)
-            pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj)
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            zztmp         = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj)
+            putaui(ji,jj) = zztmp * pwndi(ji,jj)
+            pvtaui(ji,jj) = zztmp * pwndj(ji,jj)
          END_2D
 
-         !#LB: saving the module, and x-y components, of the ai wind-stress at T-points: NOT weighted by the ice concentration !!!
-         IF(iom_use('taum_ice')) CALL iom_put('taum_ice', SQRT( putaui*putaui + pvtaui*pvtaui )*ztmp )
-         !#LB: These 2 lines below mostly here for 'STATION_ASF' test-case, otherwize "utau_oi" (U-grid) and vtau_oi" (V-grid) does the job in: [ICE/icedyn_rhg_evp.F90])
-         IF(iom_use('utau_ice')) CALL iom_put("utau_ice", putaui*ztmp)  ! utau at T-points!
-         IF(iom_use('vtau_ice')) CALL iom_put("vtau_ice", pvtaui*ztmp)  ! vtau at T-points!
-
-         !
-         DO_2D( 0, 0, 0, 0 )    ! U & V-points (same as ocean).
-            !#LB: QUESTION?? so SI3 expects wind stress vector to be provided at U & V points? Not at T-points ?
-            ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology
-            zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) )
-            zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji  ,jj+1,1) )
-            putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj  ) )
-            pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji  ,jj+1) )
-         END_2D
-         CALL lbc_lnk( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp )
+         ! outputs
+         !        LB: not weighted by the ice concentration
+         IF( iom_use('taum_ice') ) CALL iom_put( 'taum_ice', SQRT( putaui*putaui + pvtaui*pvtaui ) * zmsk00 )
+         !        LB: These 2 lines below mostly here for 'STATION_ASF' test-case
+         IF( iom_use('utau_ice') ) CALL iom_put( "utau_ice", putaui * zmsk00 )
+         IF( iom_use('vtau_ice') ) CALL iom_put( "vtau_ice", pvtaui * zmsk00 )
          !
-         IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : ', mask1=umask   &
-            &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ', mask2=vmask )
+         IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : ', mask1=tmask   &
+            &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ', mask2=tmask )
       ELSE ! ln_abl
 
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
@@ -1105,10 +1068,15 @@ CONTAINS
             pseni  (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj)
             pevpi  (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj)
          END_2D
-         pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ; ! more accurate way to obtain ssq !
+         pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! more accurate way to obtain ssq
 
       ENDIF ! ln_blk  / ln_abl
       !
+      ! outputs
+      IF( iom_use('Cd_ice') ) CALL iom_put( "Cd_ice", Cd_ice * zmsk00 )
+      IF( iom_use('Ce_ice') ) CALL iom_put( "Ce_ice", Ce_ice * zmsk00 )
+      IF( iom_use('Ch_ice') ) CALL iom_put( "Ch_ice", Ch_ice * zmsk00 )
+      !
       IF(sn_cfctl%l_prtctl)  CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ', mask1=tmask )
       !
    END SUBROUTINE blk_ice_1
diff --git a/src/OCE/SBC/sbccpl.F90 b/src/OCE/SBC/sbccpl.F90
index eddf626a522f3479f787e155a6dc76ed02900634..c3feee0dc7ad84ed850fcc57d988dc13772ef686 100644
--- a/src/OCE/SBC/sbccpl.F90
+++ b/src/OCE/SBC/sbccpl.F90
@@ -68,15 +68,15 @@ MODULE sbccpl
    INTEGER, PARAMETER ::   jpr_otx1   =  1   ! 3 atmosphere-ocean stress components on grid 1
    INTEGER, PARAMETER ::   jpr_oty1   =  2   !
    INTEGER, PARAMETER ::   jpr_otz1   =  3   !
-   INTEGER, PARAMETER ::   jpr_otx2   =  4   ! 3 atmosphere-ocean stress components on grid 2
-   INTEGER, PARAMETER ::   jpr_oty2   =  5   !
-   INTEGER, PARAMETER ::   jpr_otz2   =  6   !
+!!$   INTEGER, PARAMETER ::   jpr_otx2   =  4   ! 3 atmosphere-ocean stress components on grid 2
+!!$   INTEGER, PARAMETER ::   jpr_oty2   =  5   !
+!!$   INTEGER, PARAMETER ::   jpr_otz2   =  6   !
    INTEGER, PARAMETER ::   jpr_itx1   =  7   ! 3 atmosphere-ice   stress components on grid 1
    INTEGER, PARAMETER ::   jpr_ity1   =  8   !
    INTEGER, PARAMETER ::   jpr_itz1   =  9   !
-   INTEGER, PARAMETER ::   jpr_itx2   = 10   ! 3 atmosphere-ice   stress components on grid 2
-   INTEGER, PARAMETER ::   jpr_ity2   = 11   !
-   INTEGER, PARAMETER ::   jpr_itz2   = 12   !
+!!$   INTEGER, PARAMETER ::   jpr_itx2   = 10   ! 3 atmosphere-ice   stress components on grid 2
+!!$   INTEGER, PARAMETER ::   jpr_ity2   = 11   !
+!!$   INTEGER, PARAMETER ::   jpr_itz2   = 12   !
    INTEGER, PARAMETER ::   jpr_qsroce = 13   ! Qsr above the ocean
    INTEGER, PARAMETER ::   jpr_qsrice = 14   ! Qsr above the ice
    INTEGER, PARAMETER ::   jpr_qsrmix = 15
@@ -128,9 +128,9 @@ MODULE sbccpl
    INTEGER, PARAMETER ::   jpr_isf    = 60
    INTEGER, PARAMETER ::   jpr_icb    = 61
    INTEGER, PARAMETER ::   jpr_ts_ice = 62   ! Sea ice surface temp
-   !!INTEGER, PARAMETER ::   jpr_qtrice = 63   ! Transmitted solar thru sea-ice
+   INTEGER, PARAMETER ::   jpr_qtrice = 63   ! Transmitted solar thru sea-ice
 
-   INTEGER, PARAMETER ::   jprcv      = 62   ! total number of fields received
+   INTEGER, PARAMETER ::   jprcv      = 63   ! total number of fields received
 
    INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
    INTEGER, PARAMETER ::   jps_toce   =  2   ! ocean temperature
@@ -194,7 +194,7 @@ MODULE sbccpl
       &             sn_snd_thick1, sn_snd_cond, sn_snd_mpnd , sn_snd_sstfrz, sn_snd_ttilyr
    !                                   ! Received from the atmosphere
    TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr,  &
-      &             sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf, sn_rcv_ts_ice
+      &             sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf, sn_rcv_ts_ice, sn_rcv_qtrice
    TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf
    !                                   ! Send to waves
    TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev
@@ -202,7 +202,6 @@ MODULE sbccpl
    TYPE(FLD_C) ::   sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, &
       &             sn_rcv_wstrf, sn_rcv_wdrag, sn_rcv_charn, sn_rcv_taw, sn_rcv_bhd, sn_rcv_tusd, sn_rcv_tvsd
    !                                   ! Other namelist parameters
-!!   TYPE(FLD_C) ::   sn_rcv_qtrice
    INTEGER     ::   nn_cplmodel           ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
    LOGICAL     ::   ln_usecplmask         !  use a coupling mask file to merge data received from several models
                                           !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
@@ -281,7 +280,7 @@ CONTAINS
          &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_wstrf ,  &
          &                  sn_rcv_charn , sn_rcv_taw   , sn_rcv_bhd  , sn_rcv_tusd  , sn_rcv_tvsd,    &
          &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal  ,   &
-         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_icb  , sn_rcv_isf   , sn_rcv_ts_ice,  & !!, sn_rcv_qtrice
+         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_icb  , sn_rcv_isf   , sn_rcv_ts_ice, sn_rcv_qtrice, &
          &                  sn_rcv_mslp
 
       !!---------------------------------------------------------------------
@@ -313,7 +312,6 @@ CONTAINS
          WRITE(numout,*)'      surface stress                  = ', TRIM(sn_rcv_tau%cldes   ), ' (', TRIM(sn_rcv_tau%clcat   ), ')'
          WRITE(numout,*)'                     - referential    = ', sn_rcv_tau%clvref
          WRITE(numout,*)'                     - orientation    = ', sn_rcv_tau%clvor
-         WRITE(numout,*)'                     - mesh           = ', sn_rcv_tau%clvgrd
          WRITE(numout,*)'      non-solar heat flux sensitivity = ', TRIM(sn_rcv_dqnsdt%cldes), ' (', TRIM(sn_rcv_dqnsdt%clcat), ')'
          WRITE(numout,*)'      solar heat flux                 = ', TRIM(sn_rcv_qsr%cldes   ), ' (', TRIM(sn_rcv_qsr%clcat   ), ')'
          WRITE(numout,*)'      non-solar heat flux             = ', TRIM(sn_rcv_qns%cldes   ), ' (', TRIM(sn_rcv_qns%clcat   ), ')'
@@ -323,7 +321,7 @@ CONTAINS
          WRITE(numout,*)'      iceberg                         = ', TRIM(sn_rcv_icb%cldes   ), ' (', TRIM(sn_rcv_icb%clcat   ), ')'
          WRITE(numout,*)'      ice shelf                       = ', TRIM(sn_rcv_isf%cldes   ), ' (', TRIM(sn_rcv_isf%clcat   ), ')'
          WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
-!!       WRITE(numout,*)'      transmitted solar thru sea-ice  = ', TRIM(sn_rcv_qtrice%cldes), ' (', TRIM(sn_rcv_qtrice%clcat), ')'
+         WRITE(numout,*)'      transmitted solar thru sea-ice  = ', TRIM(sn_rcv_qtrice%cldes), ' (', TRIM(sn_rcv_qtrice%clcat), ')'
          WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
          WRITE(numout,*)'      Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')'
          WRITE(numout,*)'      surface waves:'
@@ -358,7 +356,7 @@ CONTAINS
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd
       ENDIF
       IF( lwp .AND. ln_wave) THEN                        ! control print
-      WRITE(numout,*)'      surface waves:'
+         WRITE(numout,*)'      surface waves:'
          WRITE(numout,*)'      Significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')'
          WRITE(numout,*)'      Wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')'
          WRITE(numout,*)'      Surface Stokes drift grid u     = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')'
@@ -368,8 +366,8 @@ CONTAINS
          WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')'
          WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')'
          WRITE(numout,*)'      Charnock coefficient            = ', TRIM(sn_rcv_charn%cldes ), ' (', TRIM(sn_rcv_charn%clcat ), ')'
-         WRITE(numout,*)' Transport associated to Stokes drift grid u = ', TRIM(sn_rcv_tusd%cldes ), ' (', TRIM(sn_rcv_tusd%clcat ), ')'
-         WRITE(numout,*)' Transport associated to Stokes drift grid v = ', TRIM(sn_rcv_tvsd%cldes ), ' (', TRIM(sn_rcv_tvsd%clcat ), ')'
+         WRITE(numout,*)' Transport associated to Stokes drift u = ', TRIM(sn_rcv_tusd%cldes ), ' (', TRIM(sn_rcv_tusd%clcat ), ')'
+         WRITE(numout,*)' Transport associated to Stokes drift v = ', TRIM(sn_rcv_tvsd%cldes ), ' (', TRIM(sn_rcv_tvsd%clcat ), ')'
          WRITE(numout,*)'      Bernouilli pressure head        = ', TRIM(sn_rcv_bhd%cldes   ), ' (', TRIM(sn_rcv_bhd%clcat  ), ')'
          WRITE(numout,*)'Wave to ocean momentum flux and Net wave-supported stress = ', TRIM(sn_rcv_taw%cldes ), ' (', TRIM(sn_rcv_taw%clcat ), ')'
          WRITE(numout,*)'      Surface current to waves        = ', TRIM(sn_snd_crtw%cldes  ), ' (', TRIM(sn_snd_crtw%clcat  ), ')'
@@ -399,87 +397,26 @@ CONTAINS
       srcv(jpr_otx1)%clname = 'O_OTaux1'      ! 1st ocean component on grid ONE (T or U)
       srcv(jpr_oty1)%clname = 'O_OTauy1'      ! 2nd   -      -         -     -
       srcv(jpr_otz1)%clname = 'O_OTauz1'      ! 3rd   -      -         -     -
-      srcv(jpr_otx2)%clname = 'O_OTaux2'      ! 1st ocean component on grid TWO (V)
-      srcv(jpr_oty2)%clname = 'O_OTauy2'      ! 2nd   -      -         -     -
-      srcv(jpr_otz2)%clname = 'O_OTauz2'      ! 3rd   -      -         -     -
       !
       srcv(jpr_itx1)%clname = 'O_ITaux1'      ! 1st  ice  component on grid ONE (T, F, I or U)
       srcv(jpr_ity1)%clname = 'O_ITauy1'      ! 2nd   -      -         -     -
       srcv(jpr_itz1)%clname = 'O_ITauz1'      ! 3rd   -      -         -     -
-      srcv(jpr_itx2)%clname = 'O_ITaux2'      ! 1st  ice  component on grid TWO (V)
-      srcv(jpr_ity2)%clname = 'O_ITauy2'      ! 2nd   -      -         -     -
-      srcv(jpr_itz2)%clname = 'O_ITauz2'      ! 3rd   -      -         -     -
       !
       ! Vectors: change of sign at north fold ONLY if on the local grid
       IF(       TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM( sn_rcv_tau%cldes ) == 'oce and ice'  &
            .OR. TRIM( sn_rcv_tau%cldes ) == 'mixed oce-ice' ) THEN ! avoid working with the atmospheric fields if they are not coupled
       !
-      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
-
-      !                                                           ! Set grid and action
-      SELECT CASE( TRIM( sn_rcv_tau%clvgrd ) )      !  'T', 'U,V', 'U,V,I', 'U,V,F', 'T,I', 'T,F', or 'T,U,V'
-      CASE( 'T' )
-         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
-         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
-         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
-      CASE( 'U,V' )
-         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
-         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
-         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
-         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
-         srcv(jpr_otx1:jpr_itz2)%laction = .TRUE.     ! receive oce and ice components on both grid 1 & 2
-      CASE( 'U,V,T' )
-         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
-         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
-         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'T'        ! ice components given at T-point
-         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
-         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
-      CASE( 'U,V,I' )
-         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
-         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
-         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
-         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
-         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
-      CASE( 'U,V,F' )
-         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
-         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
-         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
-         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
-         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
-      CASE( 'T,I' )
-         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
-         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
-         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
-         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
-      CASE( 'T,F' )
-         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
-         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
-         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
-         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
-      CASE( 'T,U,V' )
-         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'T'        ! oce components given at T-point
-         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
-         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
-         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1 only
-         srcv(jpr_itx1:jpr_itz2)%laction = .TRUE.     ! receive ice components on grid 1 & 2
-      CASE default
-         CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_tau%clvgrd' )
-      END SELECT
+      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz1)%nsgn = -1.
+
+      !                                                            ! Set grid and action
+      srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
+      srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
       !
       IF( TRIM( sn_rcv_tau%clvref ) == 'spherical' )   &           ! spherical: 3rd component not received
-         &     srcv( (/jpr_otz1, jpr_otz2, jpr_itz1, jpr_itz2/) )%laction = .FALSE.
-      !
-      IF( TRIM( sn_rcv_tau%clvor  ) == 'local grid' ) THEN        ! already on local grid -> no need of the second grid
-            srcv(jpr_otx2:jpr_otz2)%laction = .FALSE.
-            srcv(jpr_itx2:jpr_itz2)%laction = .FALSE.
-            srcv(jpr_oty1)%clgrid = srcv(jpr_oty2)%clgrid   ! not needed but cleaner...
-            srcv(jpr_ity1)%clgrid = srcv(jpr_ity2)%clgrid   ! not needed but cleaner...
-      ENDIF
+         &     srcv( (/jpr_otz1, jpr_itz1/) )%laction = .FALSE.
       !
-      IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN        ! 'oce and ice' case ocean stress on ocean mesh used
-         srcv(jpr_itx1:jpr_itz2)%laction = .FALSE.    ! ice components not received
-         srcv(jpr_itx1)%clgrid = 'U'                  ! ocean stress used after its transformation
-         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
+      IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN         ! 'oce and ice' case ocean stress on ocean mesh used
+         srcv(jpr_itx1:jpr_itz1)%laction = .FALSE.    ! ice components not received
       ENDIF
       ENDIF
 
@@ -612,18 +549,18 @@ CONTAINS
          ENDIF
          srcv(jpr_topm:jpr_botm)%laction = .TRUE.
       ENDIF
-!!      !                                                      ! --------------------------- !
-!!      !                                                      ! transmitted solar thru ice  !   
-!!      !                                                      ! --------------------------- !
-!!      srcv(jpr_qtrice)%clname = 'OQtr'
-!!      IF( TRIM(sn_rcv_qtrice%cldes) == 'coupled' ) THEN
-!!         IF ( TRIM( sn_rcv_qtrice%clcat ) == 'yes' ) THEN
-!!            srcv(jpr_qtrice)%nct = nn_cats_cpl
-!!         ELSE
-!!           CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qtrice%clcat should always be set to yes currently' )
-!!         ENDIF
-!!         srcv(jpr_qtrice)%laction = .TRUE.
-!!      ENDIF
+      !                                                      ! --------------------------- !
+      !                                                      ! transmitted solar thru ice  !   
+      !                                                      ! --------------------------- !
+      srcv(jpr_qtrice)%clname = 'OQtr'
+      IF( TRIM(sn_rcv_qtrice%cldes) == 'coupled' ) THEN
+         IF ( TRIM( sn_rcv_qtrice%clcat ) == 'yes' ) THEN
+            srcv(jpr_qtrice)%nct = nn_cats_cpl
+         ELSE
+           CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qtrice%clcat should always be set to yes currently' )
+         ENDIF
+         srcv(jpr_qtrice)%laction = .TRUE.
+      ENDIF
       !                                                      ! ------------------------- !
       !                                                      !    ice skin temperature   !
       !                                                      ! ------------------------- !
@@ -725,11 +662,10 @@ CONTAINS
          srcv(:)%clgrid  = 'T'       ! force default definition in case of opa <-> sas coupling
          srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
          srcv( (/jpr_qsroce, jpr_qnsoce, jpr_oemp, jpr_sflx, jpr_fice, jpr_otx1, jpr_oty1, jpr_taum/) )%laction = .TRUE.
-         srcv(jpr_otx1)%clgrid = 'U'        ! oce components given at U-point
-         srcv(jpr_oty1)%clgrid = 'V'        !           and           V-point
+         srcv(jpr_otx1)%clgrid = 'T'        ! oce components given at T-point
+         srcv(jpr_oty1)%clgrid = 'T'        
          ! Vectors: change of sign at north fold ONLY if on the local grid
          srcv( (/jpr_otx1,jpr_oty1/) )%nsgn = -1.
-         sn_rcv_tau%clvgrd = 'U,V'
          sn_rcv_tau%clvor = 'local grid'
          sn_rcv_tau%clvref = 'spherical'
          sn_rcv_emp%cldes = 'oce only'
@@ -1162,7 +1098,7 @@ CONTAINS
       !! ** Method  :   receive all fields from the atmosphere and transform
       !!              them into ocean surface boundary condition fields
       !!
-      !! ** Action  :   update  utau, vtau   ocean stress at U,V grid
+      !! ** Action  :   update  utau, vtau   ocean stress at T-point
       !!                        taum         wind stress module at T-point
       !!                        wndm         wind speed  module at T-point over free ocean or leads in presence of sea-ice
       !!                        qns          non solar heat fluxes including emp heat content    (ocean only case)
@@ -1221,39 +1157,20 @@ CONTAINS
             IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
                !                                                       ! (cartesian to spherical -> 3 to 2 components)
                !
-               CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1),   &
-                  &          srcv(jpr_otx1)%clgrid, ztx, zty )
+               CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1), 'T', ztx, zty )
                frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
                frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
                !
-               IF( srcv(jpr_otx2)%laction ) THEN
-                  CALL geo2oce( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), frcv(jpr_otz2)%z3(:,:,1),   &
-                     &          srcv(jpr_otx2)%clgrid, ztx, zty )
-                  frcv(jpr_otx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
-                  frcv(jpr_oty2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
-               ENDIF
-               !
             ENDIF
             !
             IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
                !                                                       ! (geographical to local grid -> rotate the components)
-               CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )
-               IF( srcv(jpr_otx2)%laction ) THEN
-                  CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )
-               ELSE
-                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )
-               ENDIF
+               CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), 'T', 'en->i', ztx )
+               CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), 'T', 'en->j', zty )
                frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
                frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid
             ENDIF
             !
-            IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
-               DO_2D( 0, 0, 0, 0 )                                        ! T ==> (U,V)
-                  frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
-                  frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
-               END_2D
-               CALL lbc_lnk( 'sbccpl', frcv(jpr_otx1)%z3(:,:,1), 'U',  -1.0_wp, frcv(jpr_oty1)%z3(:,:,1), 'V',  -1.0_wp )
-            ENDIF
             llnewtx = .TRUE.
          ELSE
             llnewtx = .FALSE.
@@ -1272,12 +1189,9 @@ CONTAINS
       IF( .NOT. srcv(jpr_taum)%laction ) THEN                    ! compute wind stress module from its components if not received
          ! => need to be done only when otx1 was changed
          IF( llnewtx ) THEN
-            DO_2D( 0, 0, 0, 0 )
-               zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
-               zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
-               frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
-            END_2D
-            CALL lbc_lnk( 'sbccpl', frcv(jpr_taum)%z3(:,:,1), 'T', 1.0_wp )
+            zzx = frcv(jpr_otx1)%z3(ji,jj,1)
+            zzy = frcv(jpr_oty1)%z3(ji,jj,1) 
+            frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
             llnewtau = .TRUE.
          ELSE
             llnewtau = .FALSE.
@@ -1594,7 +1508,7 @@ CONTAINS
       !! ** Action  :   return ptau_i, ptau_j, the stress over the ice
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(inout), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
-      REAL(wp), INTENT(inout), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
+      REAL(wp), INTENT(inout), DIMENSION(:,:) ::   p_tauj   !                   at T-point
       !!
       INTEGER ::   ji, jj   ! dummy loop indices
       INTEGER ::   itx      ! index of taux over ice
@@ -1616,28 +1530,16 @@ CONTAINS
             !
             IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
                !                                                       ! (cartesian to spherical -> 3 to 2 components)
-               CALL geo2oce(  frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1),   &
-                  &          srcv(jpr_itx1)%clgrid, ztx, zty )
+               CALL geo2oce(  frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1), 'T', ztx, zty )
                frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
                frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
                !
-               IF( srcv(jpr_itx2)%laction ) THEN
-                  CALL geo2oce( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), frcv(jpr_itz2)%z3(:,:,1),   &
-                     &          srcv(jpr_itx2)%clgrid, ztx, zty )
-                  frcv(jpr_itx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
-                  frcv(jpr_ity2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
-               ENDIF
-               !
             ENDIF
             !
             IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
                !                                                       ! (geographical to local grid -> rotate the components)
-               CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->i', ztx )
-               IF( srcv(jpr_itx2)%laction ) THEN
-                  CALL rot_rep( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), srcv(jpr_itx2)%clgrid, 'en->j', zty )
-               ELSE
-                  CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->j', zty )
-               ENDIF
+               CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), 'T', 'en->i', ztx )
+               CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), 'T', 'en->j', zty )
                frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
                frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 1st grid
             ENDIF
@@ -1651,29 +1553,8 @@ CONTAINS
          !                                                      ! ======================= !
          !                                                      !     put on ice grid     !
          !                                                      ! ======================= !
-         !
-         !                                                  j+1   j     -----V---F
-         ! ice stress on ice velocity point                              !       |
-         ! (C-grid ==>(U,V))                                      j      |   T   U
-         !                                                               |       |
-         !                                                   j    j-1   -I-------|
-         !                                               (for I)         |       |
-         !                                                              i-1  i   i
-         !                                                               i      i+1 (for I)
-         SELECT CASE ( srcv(jpr_itx1)%clgrid )
-         CASE( 'U' )
-            p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V)
-            p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
-         CASE( 'T' )
-            DO_2D( 0, 0, 0, 0 )                    ! T ==> (U,V)
-               ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology
-               zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) )
-               zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji  ,jj+1,1) )
-               p_taui(ji,jj) = zztmp1 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
-               p_tauj(ji,jj) = zztmp2 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
-            END_2D
-            CALL lbc_lnk( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. )
-         END SELECT
+         p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)
+         p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
 
       ENDIF
       !
@@ -1794,8 +1675,8 @@ CONTAINS
       END SELECT
 
       ! --- evaporation over ice (kg/m2/s) --- !
-      IF (ln_scale_ice_flux) THEN ! typically met-office requirements
-         IF (sn_rcv_emp%clcat == 'yes') THEN
+      IF( ln_scale_ice_flux ) THEN ! typically met-office requirements
+         IF( sn_rcv_emp%clcat == 'yes' ) THEN
             WHERE( a_i(:,:,:) > 1.e-10 )  ; zevap_ice(:,:,:) = frcv(jpr_ievp)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
             ELSEWHERE                     ; zevap_ice(:,:,:) = 0._wp
             END WHERE
@@ -1812,7 +1693,7 @@ CONTAINS
             ENDDO
          ENDIF
       ELSE
-         IF (sn_rcv_emp%clcat == 'yes') THEN
+         IF( sn_rcv_emp%clcat == 'yes' ) THEN
             zevap_ice(:,:,1:jpl) = frcv(jpr_ievp)%z3(:,:,1:jpl)
             WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:)
             ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp
@@ -1826,7 +1707,7 @@ CONTAINS
          ENDIF
       ENDIF
 
-      IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN
+      IF( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN
          ! For conservative case zemp_ice has not been defined yet. Do it now.
          zemp_ice(:,:) = zevap_ice_total(:,:) * picefr(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:)
       ENDIF
@@ -1926,13 +1807,13 @@ CONTAINS
          &                                                         - zevap_ice_total(:,:) * picefr(:,:) ) * tmask(:,:,1) )  ! 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', fwfisf(:,:) * tmask(:,:,1)                          )  ! iceshelf
       !
       !                                                      ! ========================= !
       SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !  ice topmelt and botmelt  !
       !                                                      ! ========================= !
-      CASE ('coupled')
-         IF (ln_scale_ice_flux) THEN
+      CASE( 'coupled' )
+         IF( ln_scale_ice_flux ) THEN
             WHERE( a_i(:,:,:) > 1.e-10_wp )
                qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
                qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
@@ -2213,28 +2094,29 @@ CONTAINS
          !
       ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==!
          !
-!!         SELECT CASE( TRIM( sn_rcv_qtrice%cldes ) )
-!!            !
-!!            !      ! ===> here we receive the qtr_ice_top array from the coupler
-!!         CASE ('coupled')
-!!            IF (ln_scale_ice_flux) THEN
-!!               WHERE( a_i(:,:,:) > 1.e-10_wp )
-!!                  zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
-!!               ELSEWHERE
-!!                  zqtr_ice_top(:,:,:) = 0.0_wp
-!!               ENDWHERE
-!!            ELSE
-!!               zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:)
-!!            ENDIF
-!!           
-!!            ! Add retrieved transmitted solar radiation onto the ice and total solar radiation
-!!            zqsr_ice(:,:,:) = zqsr_ice(:,:,:) + zqtr_ice_top(:,:,:)
-!!            zqsr_tot(:,:)   = zqsr_tot(:,:) + SUM( zqtr_ice_top(:,:,:) * a_i(:,:,:), dim=3 )
-!!            
-!!            !      if we are not getting this data from the coupler then assume zero (fully opaque ice)
-!!         CASE ('none')
-         zqtr_ice_top(:,:,:) = 0._wp
-!!         END SELECT
+         !
+         SELECT CASE( TRIM( sn_rcv_qtrice%cldes ) )
+            !
+            !      ! ===> here we receive the qtr_ice_top array from the coupler
+         CASE ('coupled')
+            IF (ln_scale_ice_flux) THEN
+               WHERE( a_i(:,:,:) > 1.e-10_wp )
+                  zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+               ELSEWHERE
+                  zqtr_ice_top(:,:,:) = 0.0_wp
+               ENDWHERE
+            ELSE
+               zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:)
+            ENDIF
+           
+            ! Add retrieved transmitted solar radiation onto the ice and total solar radiation
+            zqsr_ice(:,:,:) = zqsr_ice(:,:,:) + zqtr_ice_top(:,:,:)
+            zqsr_tot(:,:)   = zqsr_tot(:,:) + SUM( zqtr_ice_top(:,:,:) * a_i(:,:,:), dim=3 )
+            
+            !      if we are not getting this data from the coupler then assume zero (fully opaque ice)
+         CASE ('none')
+            zqtr_ice_top(:,:,:) = 0._wp
+         END SELECT
             !
       ENDIF
 
@@ -2403,10 +2285,10 @@ CONTAINS
           CASE( 'weighted ice' )   ;
              SELECT CASE( sn_snd_alb%clcat )
              CASE( 'yes' )
-                ztmp3(:,:,1:jpl) =  alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
+                ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
              CASE( 'no' )
                 WHERE( fr_i (:,:) > 0. )
-                   ztmp1(:,:) = SUM (  alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 )
+                   ztmp1(:,:) = SUM ( alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 )
                 ELSEWHERE
                    ztmp1(:,:) = 0.
                 END WHERE
@@ -2572,17 +2454,18 @@ CONTAINS
       IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
          !                                                   ! ------------------------- !
          !
-         !                                                  j+1   j     -----V---F
-         ! surface velocity always sent from T point                     !       |
-         !                                                        j      |   T   U
-         !                                                               |       |
-         !                                                   j    j-1   -I-------|
-         !                                               (for I)         |       |
-         !                                                              i-1  i   i
-         !                                                               i      i+1 (for I)
+         !                                                  j     -----V---F
+         ! surface velocity always sent from T point               !       |
+         !                                                  j      |   T   U
+         !                                                         |       |
+         !                                                  j-1   -I-------|
+         !                                                         |       |
+         !                                                        i-1  i   i
+         !!clem: make a new variable at T-point to replace uu and vv => uuT and vvT for instance
          IF( nn_components == jp_iam_oce ) THEN
             zotx1(:,:) = uu(:,:,1,Kmm)
             zoty1(:,:) = vv(:,:,1,Kmm)
+            !!clem : should be demi sum, no? Or uuT and vvT
          ELSE
             SELECT CASE( TRIM( sn_snd_crt%cldes ) )
             CASE( 'oce only'             )      ! C-grid ==> T
@@ -2652,42 +2535,42 @@ CONTAINS
       !                                                      !  Surface current to waves !
       !                                                      ! ------------------------- !
       IF( ssnd(jps_ocxw)%laction .OR. ssnd(jps_ocyw)%laction ) THEN
-          !
-          !                                                  j+1  j     -----V---F
-          ! surface velocity always sent from T point                    !       |
-          !                                                       j      |   T   U
-          !                                                              |       |
-          !                                                   j   j-1   -I-------|
-          !                                               (for I)        |       |
-          !                                                             i-1  i   i
-          !                                                              i      i+1 (for I)
-          SELECT CASE( TRIM( sn_snd_crtw%cldes ) )
-          CASE( 'oce only'             )      ! C-grid ==> T
-             DO_2D( 0, 0, 0, 0 )
-                zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj  ,1,Kmm) )
-                zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji , jj-1,1,Kmm) )
-             END_2D
-          CASE( 'weighted oce and ice' )      ! Ocean and Ice on C-grid ==> T
-             DO_2D( 0, 0, 0, 0 )
-                zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)
-                zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Kmm) + vv   (ji  ,jj-1,1,Kmm) ) * zfr_l(ji,jj)
-                zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
-                zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
-             END_2D
-             CALL lbc_lnk( 'sbccpl', zitx1, 'T', -1.0_wp,  zity1, 'T', -1.0_wp )
-          CASE( 'mixed oce-ice'        )      ! Ocean and Ice on C-grid ==> T
-             DO_2D( 0, 0, 0, 0 )
-                zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)   &
-                   &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
-                zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Kmm) + vv   (ji  ,jj-1,1,Kmm) ) * zfr_l(ji,jj)   &
-                   &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
-             END_2D
-          END SELECT
+         !
+         !                                                  j     -----V---F
+         ! surface velocity always sent from T point               !       |
+         !                                                  j      |   T   U
+         !                                                         |       |
+         !                                                  j-1   -I-------|
+         !                                                         |       |
+         !                                                        i-1  i   i
+         !!clem: make a new variable at T-point to replace uu and vv => uuT and vvT for instance
+         SELECT CASE( TRIM( sn_snd_crtw%cldes ) )
+         CASE( 'oce only'             )      ! C-grid ==> T
+            DO_2D( 0, 0, 0, 0 )
+               zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj  ,1,Kmm) )
+               zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji , jj-1,1,Kmm) )
+            END_2D
+         CASE( 'weighted oce and ice' )      ! Ocean and Ice on C-grid ==> T
+            DO_2D( 0, 0, 0, 0 )
+               zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)
+               zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Kmm) + vv   (ji  ,jj-1,1,Kmm) ) * zfr_l(ji,jj)
+               zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
+               zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
+            END_2D
+            CALL lbc_lnk( 'sbccpl', zitx1, 'T', -1.0_wp,  zity1, 'T', -1.0_wp )
+         CASE( 'mixed oce-ice'        )      ! Ocean and Ice on C-grid ==> T
+            DO_2D( 0, 0, 0, 0 )
+               zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)   &
+                  &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
+               zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Kmm) + vv   (ji  ,jj-1,1,Kmm) ) * zfr_l(ji,jj)   &
+                  &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
+            END_2D
+         END SELECT
          CALL lbc_lnk( 'sbccpl', zotx1, ssnd(jps_ocxw)%clgrid, -1.0_wp, zoty1, ssnd(jps_ocyw)%clgrid, -1.0_wp )
          !
          !
          IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
-         !                                                                        ! Ocean component
+            !                                                                        ! Ocean component
             CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->e', ztmp1 )       ! 1st component
             CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
             zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
@@ -2700,18 +2583,18 @@ CONTAINS
             ENDIF
          ENDIF
          !
-!         ! spherical coordinates to cartesian -> 2 components to 3 components
-!         IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN
-!            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
-!            ztmp2(:,:) = zoty1(:,:)
-!            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
-!            !
-!            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
-!               ztmp1(:,:) = zitx1(:,:)
-!               ztmp1(:,:) = zity1(:,:)
-!               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
-!            ENDIF
-!         ENDIF
+         !         ! spherical coordinates to cartesian -> 2 components to 3 components
+         !         IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN
+         !            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
+         !            ztmp2(:,:) = zoty1(:,:)
+         !            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
+         !            !
+         !            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
+         !               ztmp1(:,:) = zitx1(:,:)
+         !               ztmp1(:,:) = zity1(:,:)
+         !               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
+         !            ENDIF
+         !         ENDIF
          !
          IF( ssnd(jps_ocxw)%laction )   CALL cpl_snd( jps_ocxw, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
          IF( ssnd(jps_ocyw)%laction )   CALL cpl_snd( jps_ocyw, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
diff --git a/src/OCE/SBC/sbcflx.F90 b/src/OCE/SBC/sbcflx.F90
index c9ad4fcb5df00ee1384f884615a47ecd87470b2d..7607f8e24a2e6fd59993d8765a4e53b1dd1c4fe8 100644
--- a/src/OCE/SBC/sbcflx.F90
+++ b/src/OCE/SBC/sbcflx.F90
@@ -119,8 +119,8 @@ CONTAINS
          END DO
          !                                         ! fill sf with slf_i and control print
          CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' )
-         sf(jp_utau)%cltype = 'U'   ;   sf(jp_utau)%zsgn = -1._wp   ! vector field at U point: overwrite default definition of cltype and zsgn
-         sf(jp_vtau)%cltype = 'V'   ;   sf(jp_vtau)%zsgn = -1._wp   ! vector field at V point: overwrite default definition of cltype and zsgn
+         sf(jp_utau)%cltype = 'T'   ;   sf(jp_utau)%zsgn = -1._wp   ! vector field at T point: overwrite default definition of cltype and zsgn
+         sf(jp_vtau)%cltype = 'T'   ;   sf(jp_vtau)%zsgn = -1._wp   ! vector field at T point: overwrite default definition of cltype and zsgn
          !
       ENDIF
 
@@ -147,8 +147,8 @@ CONTAINS
       ENDIF
 #endif
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )                  ! set the ocean fluxes from read fields
-            utau(ji,jj) =   sf(jp_utau)%fnow(ji,jj,1)                              * umask(ji,jj,1)
-            vtau(ji,jj) =   sf(jp_vtau)%fnow(ji,jj,1)                              * vmask(ji,jj,1)
+            utau(ji,jj) =   sf(jp_utau)%fnow(ji,jj,1)                              * tmask(ji,jj,1)
+            vtau(ji,jj) =   sf(jp_vtau)%fnow(ji,jj,1)                              * tmask(ji,jj,1)
             qns (ji,jj) = ( sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
             emp (ji,jj) =   sf(jp_emp )%fnow(ji,jj,1)                              * tmask(ji,jj,1)
             !!sfx (ji,jj) = sf(jp_sfx )%fnow(ji,jj,1)                              * tmask(ji,jj,1)
@@ -170,19 +170,15 @@ CONTAINS
          ENDIF
          !
       ENDIF
-      !                                                           ! module of wind stress and wind speed at T-point
-      ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
+      !
+      ! module of wind stress and wind speed at T-point
       zcoef = 1. / ( zrhoa * zcdrag )
-      DO_2D( 0, 0, 0, 0 )
-         ztx = ( utau(ji-1,jj  ) + utau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( umask(ji-1,jj  ,1), umask(ji,jj,1) ) )
-         zty = ( vtau(ji  ,jj-1) + vtau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( vmask(ji  ,jj-1,1), vmask(ji,jj,1) ) )
-         zmod = SQRT( ztx * ztx + zty * zty ) * tmask(ji,jj,1)
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         zmod = SQRT( utau(ji,jj) * utau(ji,jj) + vtau(ji,jj) * vtau(ji,jj) ) * tmask(ji,jj,1)
          taum(ji,jj) = zmod
          wndm(ji,jj) = SQRT( zmod * zcoef )  !!clem: not used?
       END_2D
       !
-      CALL lbc_lnk( 'sbcflx', taum, 'T', 1._wp, wndm, 'T', 1._wp )
-      !
    END SUBROUTINE sbc_flx
 
    !!======================================================================
diff --git a/src/OCE/SBC/sbcmod.F90 b/src/OCE/SBC/sbcmod.F90
index 8881e2f4aff9d366301bb136bc3d40ea2e2c8d5f..81f018a49b41a78742f19cd35af95216d31e43a4 100644
--- a/src/OCE/SBC/sbcmod.F90
+++ b/src/OCE/SBC/sbcmod.F90
@@ -158,8 +158,8 @@ CONTAINS
       !
       IF( .NOT.ln_usr ) THEN     ! the model calendar needs some specificities (except in user defined case)
          IF( MOD( rday , rn_Dt ) /= 0. )   CALL ctl_stop( 'the time step must devide the number of second of in a day' )
-         IF( MOD( rday , 2.  ) /= 0. )   CALL ctl_stop( 'the number of second of in a day must be an even number'    )
-         IF( MOD( rn_Dt  , 2.  ) /= 0. )   CALL ctl_stop( 'the time step (in second) must be an even number'           )
+         IF( MOD( rday , 2.    ) /= 0. )   CALL ctl_stop( 'the number of second of in a day must be an even number'    )
+         IF( MOD( rn_Dt, 2.    ) /= 0. )   CALL ctl_stop( 'the time step (in second) must be an even number'           )
       ENDIF
       !                       !**  check option consistency
       !
@@ -395,16 +395,16 @@ CONTAINS
       !                                            ! ---------------------------------------- !
       IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
          !                                         ! ---------------------------------------- !
-         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields
-         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
-         qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine)
-         emp_b (:,:) = emp (:,:)
-         sfx_b (:,:) = sfx (:,:)
+         utau_b(:,:) = utauU(:,:)                        ! Swap the ocean forcing fields
+         vtau_b(:,:) = vtauV(:,:)                        ! (except at nit000 where before fields
+         qns_b (:,:) = qns  (:,:)                        !  are set at the end of the routine)
+         emp_b (:,:) = emp  (:,:)
+         sfx_b (:,:) = sfx  (:,:)
          IF( ln_rnf ) THEN
             rnf_b    (:,:  ) = rnf    (:,:  )
             rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
          ENDIF
-        !
+         !
       ENDIF
       !                                            ! ---------------------------------------- !
       !                                            !        forcing field computation         !
@@ -420,60 +420,62 @@ CONTAINS
       !
       SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
       !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx)
-      CASE( jp_usr   )     ;   CALL usrdef_sbc_oce( kt, Kbb )                        ! user defined formulation
-      CASE( jp_flx     )   ;   CALL sbc_flx       ( kt )                             ! flux formulation
+      CASE( jp_usr   )     ;   CALL usrdef_sbc_oce( kt, Kbb )                         ! user defined formulation
+      CASE( jp_flx     )   ;   CALL sbc_flx       ( kt )                              ! flux formulation
       CASE( jp_blk     )
          IF( ll_sas    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )   ! OCE-SAS coupling: SAS receiving fields from OCE
 !!!!!!!!!!! ATTENTION:ln_wave is not only used for oasis coupling !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          IF( ln_wave )   THEN
-             IF ( lk_oasis )  CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice, Kbb, Kmm )   ! OCE-wave coupling
+             IF ( lk_oasis )   CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )   ! OCE-wave coupling
              CALL sbc_wave ( kt, Kmm )
          ENDIF
-                               CALL sbc_blk       ( kt )                    ! bulk formulation for the ocean
+                               CALL sbc_blk       ( kt )                              ! bulk formulation for the ocean
                                !
       CASE( jp_abl     )
          IF( ll_sas    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )   ! OCE-SAS coupling: SAS receiving fields from OCE
-                               CALL sbc_abl       ( kt )                    ! ABL  formulation for the ocean
+                               CALL sbc_abl       ( kt )                              ! ABL  formulation for the ocean
                                !
       CASE( jp_purecpl )   ;   CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )   ! pure coupled formulation
       CASE( jp_none    )
-         IF( ll_opa    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! OCE-SAS coupling: OCE receiving fields from SAS
+         IF( ll_opa    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )   ! OCE-SAS coupling: OCE receiving fields from SAS
       END SELECT
       !
-      IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! forced-coupled mixed formulation after forcing
+      IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )   ! forced-coupled mixed formulation after forcing
       !
       IF( ln_wave .AND. ln_tauoc )  THEN            ! Wave stress reduction
-         DO_2D( 0, 0, 0, 0)
-            utau(ji,jj) = utau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji+1,jj) ) * 0.5_wp
-            vtau(ji,jj) = vtau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji,jj+1) ) * 0.5_wp
-         END_2D
          !
-         CALL lbc_lnk( 'sbcwave', utau, 'U', -1. )
-         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. )
-         !
-         taum(:,:) = taum(:,:)*tauoc_wave(:,:)
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            utau(ji,jj) = utau(ji,jj) * tauoc_wave(ji,jj)
+            vtau(ji,jj) = vtau(ji,jj) * tauoc_wave(ji,jj)
+            taum(ji,jj) = taum(ji,jj) * tauoc_wave(ji,jj)
+         END_2D
          !
          IF( kt == nit000 )   CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.',   &
             &                                'If not requested select ln_tauoc=.false.' )
          !
       ELSEIF( ln_wave .AND. ln_taw ) THEN                  ! Wave stress reduction
-         utau(:,:) = utau(:,:) - tawx(:,:) + twox(:,:)
-         vtau(:,:) = vtau(:,:) - tawy(:,:) + twoy(:,:)
-         CALL lbc_lnk( 'sbcwave', utau, 'U', -1. )
-         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. )
-         !
-         DO_2D( 0, 0, 0, 0)
-             taum(ji,jj) = sqrt((.5*(utau(ji-1,jj)+utau(ji,jj)))**2 + (.5*(vtau(ji,jj-1)+vtau(ji,jj)))**2)
+         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+            utau(ji,jj) = utau(ji,jj) - tawx(ji,jj) + twox(ji,jj)
+            vtau(ji,jj) = vtau(ji,jj) - tawy(ji,jj) + twoy(ji,jj)
+            taum(ji,jj) = SQRT( utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj) )
          END_2D
          !
          IF( kt == nit000 )   CALL ctl_warn( 'sbc: You are subtracting the wave stress to the ocean.',   &
             &                                'If not requested select ln_taw=.false.' )
          !
       ENDIF
-      CALL lbc_lnk( 'sbcmod', taum(:,:), 'T', 1. )
+      !
       !
       IF( ln_icebergs ) THEN  ! save pure stresses (with no ice-ocean stress) for use by icebergs
-         utau_icb(:,:) = utau(:,:) ; vtau_icb(:,:) = vtau(:,:) 
+         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
+         !      and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
+         DO_2D( 0, 0, 0, 0 )
+            utau_icb(ji,jj) = 0.5_wp * ( utau(ji,jj) + utau(ji+1,jj) ) * &
+               &                       ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji+1,jj,1) )
+            vtau_icb(ji,jj) = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj+1) ) * &
+               &                       ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji,jj+1,1) )
+         END_2D
+         CALL lbc_lnk( 'sbcmod', utau_icb, 'U', -1.0_wp, vtau_icb, 'V', -1.0_wp )
       ENDIF
       !
       !                                            !==  Misc. Options  ==!
@@ -507,9 +509,6 @@ CONTAINS
       ! Should not be run if ln_diurnal_only
       IF( l_sbc_clo      )   CALL sbc_clo( kt )
 
-!!$!RBbug do not understand why see ticket 667
-!!$!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why.
-!!$      CALL lbc_lnk( 'sbcmod', emp, 'T', 1.0_wp )
       IF( ll_wd ) THEN     ! If near WAD point limit the flux for now
          zthscl = atanh(rn_wd_sbcfra)                     ! taper frac default is .999
          zwdht(:,:) = ssh(:,:,Kmm) + ht_0(:,:) - rn_wdmin1   ! do this calc of water
@@ -534,6 +533,16 @@ CONTAINS
             emp  (:,:) =  emp(:,:)  * zwght(:,:)
          END WHERE
       ENDIF
+      ! --- calculate utau and vtau on U,V-points --- !
+      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
+      !      and the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
+      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         utauU  (ji,jj) = 0.5_wp * ( utau(ji,jj) + utau(ji+1,jj) ) * &
+            &                      ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji+1,jj,1) )
+         vtauV  (ji,jj) = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj+1) ) * &
+            &                      ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1), tmask(ji,jj+1,1) )
+      END_2D
+      IF( nn_hls == 1 )   CALL lbc_lnk( 'sbcmod', utauU, 'U', -1.0_wp, vtauV, 'V', -1.0_wp )
       !
       IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
          !                                             ! ---------------------------------------- !
@@ -543,27 +552,28 @@ CONTAINS
          IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN            !* MLF: Restart: read in restart file
 #endif
             IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields read in the restart file'
-            CALL iom_get( numror, jpdom_auto, 'utau_b', utau_b )   ! i-stress
-            CALL iom_get( numror, jpdom_auto, 'vtau_b', vtau_b )   ! j-stress
-            CALL iom_get( numror, jpdom_auto,  'qns_b',  qns_b )   ! non solar heat flux
-            CALL iom_get( numror, jpdom_auto,  'emp_b',  emp_b )   ! freshwater flux
+            CALL iom_get( numror, jpdom_auto, 'utau_b', utau_b, cd_type = 'U', psgn = -1._wp )   ! i-stress
+            CALL iom_get( numror, jpdom_auto, 'vtau_b', vtau_b, cd_type = 'V', psgn = -1._wp )   ! j-stress
+            CALL iom_get( numror, jpdom_auto,  'qns_b',  qns_b, cd_type = 'T', psgn =  1._wp )   ! non solar heat flux
+            CALL iom_get( numror, jpdom_auto,  'emp_b',  emp_b, cd_type = 'T', psgn =  1._wp )   ! freshwater flux
             ! NB: The 3D heat content due to qsr forcing (qsr_hc_b) is treated in traqsr
             ! To ensure restart capability with 3.3x/3.4 restart files    !! to be removed in v3.6
             IF( iom_varid( numror, 'sfx_b', ldstop = .FALSE. ) > 0 ) THEN
-               CALL iom_get( numror, jpdom_auto, 'sfx_b', sfx_b )  ! before salt flux (T-point)
+               CALL iom_get( numror, jpdom_auto, 'sfx_b', sfx_b, cd_type = 'T', psgn = 1._wp )   ! before salt flux (T-point)
             ELSE
                sfx_b (:,:) = sfx(:,:)
             ENDIF
          ELSE                                                   !* no restart: set from nit000 values
             IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
-            utau_b(:,:) = utau(:,:)
-            vtau_b(:,:) = vtau(:,:)
+            utau_b(:,:) = utauU(:,:)
+            vtau_b(:,:) = vtauV(:,:)
             qns_b (:,:) = qns (:,:)
             emp_b (:,:) = emp (:,:)
             sfx_b (:,:) = sfx (:,:)
          ENDIF
       ENDIF
       !
+      !
 #if defined key_RK3
       !                                                ! ---------------------------------------- !
       IF( lrst_oce .AND. lk_SWE ) THEN                 !   RK3: Write in the ocean restart file   !
@@ -578,13 +588,13 @@ CONTAINS
          IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   &
             &                    'at it= ', kt,' date= ', ndastp
          IF(lwp) WRITE(numout,*) '~~~~'
-         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau )
-         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau )
-         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  )
+         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utauU )
+         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtauV )
+         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns   )
          ! The 3D heat content due to qsr forcing is treated in traqsr
          ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
-         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  )
-         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx  )
+         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp   )
+         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx   )
       ENDIF
       !                                                ! ---------------------------------------- !
       !                                                !        Outputs and control print         !
@@ -628,8 +638,8 @@ CONTAINS
          CALL prt_ctl(tab3d_1=tmask               , clinfo1=' tmask    - : ', mask1=tmask, kdim=jpk )
          CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' sst      - : ', mask1=tmask, kdim=1   )
          CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sss      - : ', mask1=tmask, kdim=1   )
-         CALL prt_ctl(tab2d_1=utau                , clinfo1=' utau     - : ', mask1=umask,                      &
-            &         tab2d_2=vtau                , clinfo2=' vtau     - : ', mask2=vmask )
+         CALL prt_ctl(tab2d_1=utau                , clinfo1=' utau     - : ', mask1=tmask,                      &
+            &         tab2d_2=vtau                , clinfo2=' vtau     - : ', mask2=tmask )
       ENDIF
 
       IF( kt == nitend )   CALL sbc_final         ! Close down surface module if necessary
diff --git a/src/OCE/TRD/trddyn.F90 b/src/OCE/TRD/trddyn.F90
index 8816056e52b877dbd4118eaba52ceb326bc1c8b4..32e95a8653d448e0dd2afdc350c7a05efe02811a 100644
--- a/src/OCE/TRD/trddyn.F90
+++ b/src/OCE/TRD/trddyn.F90
@@ -140,8 +140,10 @@ CONTAINS
                               !
                               !                                    ! wind stress trends
                               ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) )
-                              z2dx(:,:) = ( utau_b(:,:) + utau(:,:) ) / ( e3u(:,:,1,Kmm) * rho0 )
-                              z2dy(:,:) = ( vtau_b(:,:) + vtau(:,:) ) / ( e3v(:,:,1,Kmm) * rho0 )
+                              DO_2D( 0, 0, 0, 0 )
+                                 z2dx(ji,jj) = ( utau_b(ji,jj) + utauU(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rho0 )
+                                 z2dy(ji,jj) = ( vtau_b(ji,jj) + vtauV(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rho0 )
+                              END_2D
                               CALL iom_put( "utrd_tau", z2dx )
                               CALL iom_put( "vtrd_tau", z2dy )
                               DEALLOCATE( z2dx , z2dy )
diff --git a/src/OCE/TRD/trdglo.F90 b/src/OCE/TRD/trdglo.F90
index 0e4a2aea0f1d16b4f7299c6d414d292d8b45b027..a19770158b0ef874869a857f323bfcdec3818151 100644
--- a/src/OCE/TRD/trdglo.F90
+++ b/src/OCE/TRD/trdglo.F90
@@ -127,10 +127,10 @@ CONTAINS
             IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
                z1_2rho0 = 0.5_wp / rho0
                DO_2D( 1, 0, 1, 0 )
-                  zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
-                     &                                                     * z1_2rho0       * e1e2u(ji,jj)
-                  zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
-                     &                                                     * z1_2rho0       * e1e2v(ji,jj)
+                  zvt = ( utau_b(ji,jj) + utauU(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
+                     &                                                      * z1_2rho0       * e1e2u(ji,jj)
+                  zvs = ( vtau_b(ji,jj) + vtauV(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
+                     &                                                      * z1_2rho0       * e1e2v(ji,jj)
                   umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
                   vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
                   hke(jpdyn_tau) = hke(jpdyn_tau) + uu(ji,jj,1,Kmm) * zvt + vv(ji,jj,1,Kmm) * zvs
diff --git a/src/OCE/TRD/trdken.F90 b/src/OCE/TRD/trdken.F90
index d97b3eb80a37e9db31884014b55435e2ba3547cc..1920e3e72199202116ed3103f0649a1c241e3886 100644
--- a/src/OCE/TRD/trdken.F90
+++ b/src/OCE/TRD/trdken.F90
@@ -118,16 +118,18 @@ CONTAINS
          CASE( jpdyn_ldf )   ;   CALL iom_put( "ketrd_ldf"   , zke )    ! lateral diffusion
          CASE( jpdyn_zdf )   ;   CALL iom_put( "ketrd_zdf"   , zke )    ! vertical diffusion 
          !                   !                                          ! wind stress trends
-                                 ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) , zke2d(jpi,jpj) )
-                           z2dx(:,:) = uu(:,:,1,Kmm) * ( utau_b(:,:) + utau(:,:) ) * e1e2u(:,:) * umask(:,:,1)
-                           z2dy(:,:) = vv(:,:,1,Kmm) * ( vtau_b(:,:) + vtau(:,:) ) * e1e2v(:,:) * vmask(:,:,1)
-                           zke2d(1:nn_hls,:) = 0._wp   ;   zke2d(:,1:nn_hls) = 0._wp
-                           DO_2D( 0, nn_hls, 0, nn_hls )
-                              zke2d(ji,jj) = r1_rho0 * 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
-                              &                                   + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,1)
-                           END_2D
-                                 CALL iom_put( "ketrd_tau"   , zke2d )  ! 
-                                 DEALLOCATE( z2dx , z2dy , zke2d )
+            ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) , zke2d(jpi,jpj) )
+            DO_2D( 1, 0, 1, 0 )
+               z2dx(ji,jj) = uu(ji,jj,1,Kmm) * ( utau_b(ji,jj) + utauU(ji,jj) ) * e1e2u(ji,jj) * umask(ji,jj,1)
+               z2dy(ji,jj) = vv(ji,jj,1,Kmm) * ( vtau_b(ji,jj) + vtauV(ji,jj) ) * e1e2v(ji,jj) * vmask(ji,jj,1)
+            END_2D
+            DO_2D( 0, 0, 0, 0 )
+               zke2d(ji,jj) = r1_rho0 * 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
+                  &                                + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,1)
+            END_2D
+            CALL iom_put( "ketrd_tau"   , zke2d )  ! 
+            DEALLOCATE( z2dx , z2dy , zke2d )
+            !
          CASE( jpdyn_bfr )   ;   CALL iom_put( "ketrd_bfr"   , zke )    ! bottom friction (explicit case) 
 !!gm TO BE DONE properly
 !!gm only valid if ln_drgimp=F otherwise the bottom stress as to be recomputed at the end of the computation....
diff --git a/src/OCE/TRD/trdvor.F90 b/src/OCE/TRD/trdvor.F90
index 974f5b1a9b06f1a16fcdbee8d7457b71f9a5ea6c..8515d666d398a18fc170669dde8a9c13160a5cbd 100644
--- a/src/OCE/TRD/trdvor.F90
+++ b/src/OCE/TRD/trdvor.F90
@@ -105,9 +105,9 @@ CONTAINS
       CASE( jpdyn_zad )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_zad, Kmm )   ! Vertical Advection 
       CASE( jpdyn_spg )   ;   CALL trd_vor_zint( putrd, pvtrd, jpvor_spg, Kmm )   ! Surface Pressure Grad. 
       CASE( jpdyn_zdf )                                                           ! Vertical Diffusion 
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )                                                               ! wind stress trends
-            ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rho0 )
-            ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rho0 )
+         DO_2D( 0, 1, 0, 1 )                                                      ! wind stress trends
+            ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utauU(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rho0 )
+            ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtauV(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rho0 )
          END_2D
          CALL trd_vor_zint( putrd, pvtrd, jpvor_zdf, Kmm )                             ! zdf trend including surf./bot. stresses 
          CALL trd_vor_zint( ztswu, ztswv, jpvor_swf, Kmm )                             ! surface wind stress 
@@ -165,7 +165,7 @@ CONTAINS
       SELECT CASE( ktrd ) 
       !
       CASE( jpvor_bfr )        ! bottom friction
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 1, 0, 1 )
             ikbu = mbkv(ji,jj)
             ikbv = mbkv(ji,jj)            
             zudpvor(ji,jj) = putrdvor(ji,jj) * e3u(ji,jj,ikbu,Kmm) * e1u(ji,jj) * umask(ji,jj,ikbu)
@@ -173,14 +173,18 @@ CONTAINS
          END_2D
          !
       CASE( jpvor_swf )        ! wind stress
-         zudpvor(:,:) = putrdvor(:,:) * e3u(:,:,1,Kmm) * e1u(:,:) * umask(:,:,1)
-         zvdpvor(:,:) = pvtrdvor(:,:) * e3v(:,:,1,Kmm) * e2v(:,:) * vmask(:,:,1)
+         DO_2D( 0, 1, 0, 1 )
+            zudpvor(ji,jj) = putrdvor(ji,jj) * e3u(ji,jj,1,Kmm) * e1u(ji,jj) * umask(ji,jj,1)
+            zvdpvor(ji,jj) = pvtrdvor(ji,jj) * e3v(ji,jj,1,Kmm) * e2v(ji,jj) * vmask(ji,jj,1)
+         END_2D
          !
       END SELECT
 
       ! Average except for Beta.V
-      zudpvor(:,:) = zudpvor(:,:) * r1_hu(:,:,Kmm)
-      zvdpvor(:,:) = zvdpvor(:,:) * r1_hv(:,:,Kmm)
+      DO_2D( 0, 1, 0, 1 )
+         zudpvor(ji,jj) = zudpvor(ji,jj) * r1_hu(ji,jj,Kmm)
+         zvdpvor(ji,jj) = zvdpvor(ji,jj) * r1_hv(ji,jj,Kmm)
+      END_2D
    
       ! Curl
       DO_2D( 0, 0, 0, 0 )
@@ -238,10 +242,11 @@ CONTAINS
       !  I vertical integration of 3D trends
       !  =====================================
       ! putrdvor and pvtrdvor terms
-      DO jk = 1,jpk
-        zudpvor(:,:) = zudpvor(:,:) + putrdvor(:,:,jk) * e3u(:,:,jk,Kmm) * e1u(:,:) * umask(:,:,jk)
-        zvdpvor(:,:) = zvdpvor(:,:) + pvtrdvor(:,:,jk) * e3v(:,:,jk,Kmm) * e2v(:,:) * vmask(:,:,jk)
-      END DO
+      zudpvor(:,:) = 0._wp ; zvdpvor(:,:) = 0._wp
+      DO_3D( 0, 1, 0, 1, 1, jpk )
+        zudpvor(ji,jj) = zudpvor(ji,jj) + putrdvor(ji,jj,jk) * e3u(ji,jj,jk,Kmm) * e1u(ji,jj) * umask(ji,jj,jk)
+        zvdpvor(ji,jj) = zvdpvor(ji,jj) + pvtrdvor(ji,jj,jk) * e3v(ji,jj,jk,Kmm) * e2v(ji,jj) * vmask(ji,jj,jk)
+      END_3D
 
       ! Planetary vorticity: 2nd computation (Beta.V term) store the vertical sum
       ! as Beta.V term need intergration, not average
@@ -254,8 +259,10 @@ CONTAINS
       ENDIF
       !
       ! Average 
-      zudpvor(:,:) = zudpvor(:,:) * r1_hu(:,:,Kmm)
-      zvdpvor(:,:) = zvdpvor(:,:) * r1_hv(:,:,Kmm)
+      DO_2D( 0, 1, 0, 1 )
+         zudpvor(ji,jj) = zudpvor(ji,jj) * r1_hu(ji,jj,Kmm)
+         zvdpvor(ji,jj) = zvdpvor(ji,jj) * r1_hv(ji,jj,Kmm)
+      END_2D
       !
       ! Curl
       DO_2D( 0, 0, 0, 0 )
diff --git a/src/OCE/USR/usrdef_sbc.F90 b/src/OCE/USR/usrdef_sbc.F90
index ce5b785a973082a9dd81b952667d965cd21a2bfc..2af3f30d2789fe3cfafadda55d4141c9c38c95b8 100644
--- a/src/OCE/USR/usrdef_sbc.F90
+++ b/src/OCE/USR/usrdef_sbc.F90
@@ -166,19 +166,17 @@ CONTAINS
       ! seasonal oscillation intensity
       ztau_sais = 0.015
       ztaun = ztau - ztau_sais * COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi )
-      DO_2D( 1, 1, 1, 1 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
         ! domain from 15deg to 50deg and 1/2 period along 14deg
         ! so 5/4 of half period with seasonal cycle
-        utau(ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) )
-        vtau(ji,jj) =   ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) )
+        utau(ji,jj) = - ztaun * SIN( rpi * (gphit(ji,jj) - 15.) / (29.-15.) )
+        vtau(ji,jj) =   ztaun * SIN( rpi * (gphit(ji,jj) - 15.) / (29.-15.) )
       END_2D
 
       ! module of wind stress and wind speed at T-point
       zcoef = 1. / ( zrhoa * zcdrag ) 
-      DO_2D( 0, 0, 0, 0 )
-         ztx = utau(ji-1,jj  ) + utau(ji,jj) 
-         zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
-         zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         zmod = SQRT( utau(ji,jj) * utau(ji,jj) + vtau(ji,jj) * vtau(ji,jj) )
          taum(ji,jj) = zmod
          wndm(ji,jj) = SQRT( zmod * zcoef )
       END_2D
diff --git a/src/OCE/ZDF/zdfosm.F90 b/src/OCE/ZDF/zdfosm.F90
index df648a8c08f42edbfb5ec326af5ce3e9f96b1b08..490a5d56db2dbadc90b52939ed5b68ebe1ae63c1 100644
--- a/src/OCE/ZDF/zdfosm.F90
+++ b/src/OCE/ZDF/zdfosm.F90
@@ -486,8 +486,8 @@ CONTAINS
             &             grav  * zbeta * swsav(ji,jj)                      ! OBSBL
       END_2D
       DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
-         suw0(ji,jj)    = -0.5_wp * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1)   ! Surface upward velocity fluxes
-         zvw0           = -0.5_wp * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1)
+         suw0(ji,jj)    = - utau(ji,jj) * r1_rho0 * tmask(ji,jj,1)          ! Surface upward velocity fluxes
+         zvw0           = - vtau(ji,jj) * r1_rho0 * tmask(ji,jj,1)
          sustar(ji,jj)  = MAX( SQRT( SQRT( suw0(ji,jj) * suw0(ji,jj) + zvw0 * zvw0 ) ),   &   ! Friction velocity (sustar), at
             &                  1e-8_wp )                                                      !    T-point : LMD94 eq. 2
          scos_wind(ji,jj) = -1.0_wp * suw0(ji,jj) / ( sustar(ji,jj) * sustar(ji,jj) )
diff --git a/src/OCE/ZDF/zdftke.F90 b/src/OCE/ZDF/zdftke.F90
index 68dd9f6ed89537ea357a57cb963ec863cfb94e1d..4d3adbd3e83c02c0f3c65f59097220966c61052a 100644
--- a/src/OCE/ZDF/zdftke.F90
+++ b/src/OCE/ZDF/zdftke.F90
@@ -210,7 +210,7 @@ CONTAINS
       REAL(wp) ::   zcdrag = 1.5e-3            ! drag coefficient
       REAL(wp) ::   zbbrau, zbbirau, zri       ! local scalars
       REAL(wp) ::   zfact1, zfact2, zfact3     !   -      -
-      REAL(wp) ::   ztx2  , zty2  , zcof       !   -      -
+      REAL(wp) ::   zcof                       !   -      -
       REAL(wp) ::   ztau  , zdif               !   -      -
       REAL(wp) ::   zus   , zwlc  , zind       !   -      -
       REAL(wp) ::   zzd_up, zzd_lw             !   -      -
@@ -302,8 +302,8 @@ CONTAINS
 !  Projection of Stokes drift in the wind stress direction
 !
             DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
-                  ztaui   = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) )
-                  ztauj   = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) )
+                  ztaui   = utau(ji,jj)
+                  ztauj   = vtau(ji,jj)
                   z1_norm = 1._wp / MAX( SQRT(ztaui*ztaui+ztauj*ztauj), 1.e-12 ) * tmask(ji,jj,1)
                   zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2
             END_2D
@@ -483,9 +483,7 @@ CONTAINS
          END_2D
       ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
          DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
-            ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
-            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
-            ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress
+            ztau = SQRT( utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj) ) * tmask(ji,jj,1)    ! module of the mean stress
             zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
             zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
             en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
diff --git a/src/OCE/stp2d.F90 b/src/OCE/stp2d.F90
index 96170dc207405395a7248c429d57e8516a227d51..06c08b84b7ce3051c0bbb3a5a99f7aa44b99c548 100644
--- a/src/OCE/stp2d.F90
+++ b/src/OCE/stp2d.F90
@@ -151,14 +151,14 @@ CONTAINS
       !                             !* wind forcing *!
       IF( ln_bt_fw ) THEN
          DO_2D( 0, 0, 0, 0 )
-            Ue_rhs(ji,jj) =  Ue_rhs(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kbb)
-            Ve_rhs(ji,jj) =  Ve_rhs(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kbb)
+            Ue_rhs(ji,jj) =  Ue_rhs(ji,jj) + r1_rho0 * utauU(ji,jj) * r1_hu(ji,jj,Kbb)
+            Ve_rhs(ji,jj) =  Ve_rhs(ji,jj) + r1_rho0 * vtauV(ji,jj) * r1_hv(ji,jj,Kbb)
          END_2D
       ELSE
          zztmp = r1_rho0 * r1_2
          DO_2D( 0, 0, 0, 0 )
-            Ue_rhs(ji,jj) =  Ue_rhs(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kbb)
-            Ve_rhs(ji,jj) =  Ve_rhs(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kbb)
+            Ue_rhs(ji,jj) =  Ue_rhs(ji,jj) + zztmp * ( utau_b(ji,jj) + utauU(ji,jj) ) * r1_hu(ji,jj,Kbb)
+            Ve_rhs(ji,jj) =  Ve_rhs(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtauV(ji,jj) ) * r1_hv(ji,jj,Kbb)
          END_2D
       ENDIF
       !
diff --git a/src/SAS/diawri.F90 b/src/SAS/diawri.F90
index 68c13c0b4bf688444072644edc48f3936fe620c1..8d0c2ece3c67635cd2aa4e3e619c95e181949c68 100644
--- a/src/SAS/diawri.F90
+++ b/src/SAS/diawri.F90
@@ -317,22 +317,22 @@ CONTAINS
 		 !
 		 ENDIF
 !
+         CALL histdef( nid_T, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
+            &          jpi, jpj, nh_T, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
+         CALL histdef( nid_T, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
+            &          jpi, jpj, nh_T, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
 
          CALL histend( nid_T, snc4chunks=snc4set )
 
          !                                                                                      !!! nid_U : 3D
          CALL histdef( nid_U, "ssu_m", "Velocity component in x-direction", "m/s"   ,         &  ! ssu
             &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
-         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
-            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
 
          CALL histend( nid_U, snc4chunks=snc4set )
 
          !                                                                                      !!! nid_V : 3D
          CALL histdef( nid_V, "ssv_m", "Velocity component in y-direction", "m/s",            &  ! ssv_m
             &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
-         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
-            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
 
          CALL histend( nid_V, snc4chunks=snc4set )
 
@@ -366,6 +366,8 @@ CONTAINS
       CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
       CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
       CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
+      CALL histwrite( nid_T, "sozotaux", it, utau          , ndim_hT, ndex_hT )   ! i-wind stress
+      CALL histwrite( nid_T, "sometauy", it, vtau          , ndim_hT, ndex_hT )   ! j-wind stress
 !
       IF( ln_abl ) THEN 
 	     ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
@@ -393,11 +395,9 @@ CONTAINS
 
          ! Write fields on U grid
       CALL histwrite( nid_U, "ssu_m"   , it, ssu_m         , ndim_hU, ndex_hU )   ! i-current speed
-      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
 
          ! Write fields on V grid
       CALL histwrite( nid_V, "ssv_m"   , it, ssv_m         , ndim_hV, ndex_hV )   ! j-current speed
-      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
 
       ! 3. Close all files
       ! ---------------------------------------
diff --git a/src/SWE/stpmlf.F90 b/src/SWE/stpmlf.F90
index 5494248c2da3d7d583c2377f3228ce1d73cbae4c..4c49e55aaaab7aa745a49bcbfee63991a4c92743 100644
--- a/src/SWE/stpmlf.F90
+++ b/src/SWE/stpmlf.F90
@@ -114,9 +114,9 @@ CONTAINS
          zrhs_u =        - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
          zrhs_v =        - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
          !                                                              ! wind stress and layer friction
-         zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
+         zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utauU(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
             &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
-         zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
+         zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtauV(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
             &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
          !                                                              ! ==> RHS
          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
diff --git a/src/SWE/stprk3.F90 b/src/SWE/stprk3.F90
index 0e1a1910fc971a67837aafea9430c53efa272d2e..2fded6aea76692509a6c7d870a9c91b6497f4928 100644
--- a/src/SWE/stprk3.F90
+++ b/src/SWE/stprk3.F90
@@ -135,9 +135,9 @@ CONTAINS
          zrhs_v =        - grav    * ( ssh(ji,jj+1,Nbb) - ssh(ji,jj,Nbb) ) * r1_e2v(ji,jj)
 #if defined key_RK3all
          !                                          ! wind stress and layer friction
-         zrhs_u = zrhs_u + r1_rho0 * ( z5_6*utau_b(ji,jj) + (1._wp - z5_6)*utau(ji,jj) ) / e3u(ji,jj,jk,Nbb)   &
+         zrhs_u = zrhs_u + r1_rho0 * ( z5_6*utau_b(ji,jj) + (1._wp - z5_6)*utauU(ji,jj) ) / e3u(ji,jj,jk,Nbb)   &
             &            - rn_rfr  * uu(ji,jj,jk,Nbb)
-         zrhs_v = zrhs_v + r1_rho0 * ( z5_6*vtau_b(ji,jj) + (1._wp - z5_6)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb)   &
+         zrhs_v = zrhs_v + r1_rho0 * ( z5_6*vtau_b(ji,jj) + (1._wp - z5_6)*vtauV(ji,jj) ) / e3v(ji,jj,jk,Nbb)   &
             &            - rn_rfr  * vv(ji,jj,jk,Nbb)
 #endif
          !                                          ! ==> RHS
@@ -201,9 +201,9 @@ CONTAINS
          zrhs_v =        - grav    * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
 #if defined key_RK3all
          !                                          ! wind stress and layer friction
-         zrhs_u = zrhs_u + r1_rho0 * ( z3_4*utau_b(ji,jj) + (1._wp - z3_4)*utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
+         zrhs_u = zrhs_u + r1_rho0 * ( z3_4*utau_b(ji,jj) + (1._wp - z3_4)*utauU(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
             &            - rn_rfr  * uu(ji,jj,jk,Nbb)
-         zrhs_v = zrhs_v + r1_rho0 * ( z3_4*vtau_b(ji,jj) + (1._wp - z3_4)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
+         zrhs_v = zrhs_v + r1_rho0 * ( z3_4*vtau_b(ji,jj) + (1._wp - z3_4)*vtauV(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
             &            - rn_rfr  * vv(ji,jj,jk,Nbb)
 #endif
          !                                          ! ==> RHS
@@ -265,9 +265,9 @@ CONTAINS
          zrhs_u =        - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
          zrhs_v =        - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
          !                                          ! wind stress and layer friction
-         zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
+         zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utauU(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
             &            - rn_rfr   * uu(ji,jj,jk,Nbb)
-         zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
+         zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtauV(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
             &            - rn_rfr   * vv(ji,jj,jk,Nbb)
          !                                          ! ==> RHS
          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
diff --git a/tests/BENCH/MY_SRC/usrdef_sbc.F90 b/tests/BENCH/MY_SRC/usrdef_sbc.F90
index e9ac757b669a92d58a27241281a7abf2db635ba8..4bd1292280807e173740a890e4ed19048225b51f 100644
--- a/tests/BENCH/MY_SRC/usrdef_sbc.F90
+++ b/tests/BENCH/MY_SRC/usrdef_sbc.F90
@@ -109,7 +109,7 @@ CONTAINS
          vtau_ice(ji,jj) = 0.1_wp + zztmp
       END_2D
 
-      CALL lbc_lnk( 'usrdef_sbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
+      CALL lbc_lnk( 'usrdef_sbc', utau_ice, 'T', -1., vtau_ice, 'T', -1. )
 #endif
       !
    END SUBROUTINE usrdef_sbc_ice_tau
diff --git a/tests/CANAL/EXPREF/file_def_nemo-oce.xml b/tests/CANAL/EXPREF/file_def_nemo-oce.xml
index 689b95dd6d9b1e960650cd83922b4ed7448d8ae5..45712ad6d87398fee622f6ab7574d39787508da4 100644
--- a/tests/CANAL/EXPREF/file_def_nemo-oce.xml
+++ b/tests/CANAL/EXPREF/file_def_nemo-oce.xml
@@ -20,15 +20,15 @@
 	  <field field_ref="ssrelpotvor"  />
 	  <field field_ref="saltc"  />
 	  <field field_ref="salt2c"  />
+	  <field field_ref="utau"  />
+	  <field field_ref="vtau"  />
 	</file>
 
 	<file id="file3" name_suffix="_grid_U" description="ocean U grid variables" >
-	  <field field_ref="utau"  />
 	  <field field_ref="uoce" />
 	</file>
 	
 	<file id="file4" name_suffix="_grid_V" description="ocean V grid variables" >
-	  <field field_ref="vtau"  />
 	  <field field_ref="voce" />
 	</file>
 	
diff --git a/tests/DIA_GPU/EXPREF/file_def_nemo-oce.xml b/tests/DIA_GPU/EXPREF/file_def_nemo-oce.xml
index c043f56fe5e9ae6fc252ae34d78b6c2660381ec9..04590b55f29789842576a1dff286ad9ec572d172 100644
--- a/tests/DIA_GPU/EXPREF/file_def_nemo-oce.xml
+++ b/tests/DIA_GPU/EXPREF/file_def_nemo-oce.xml
@@ -34,6 +34,8 @@
 	  <field field_ref="qt_oce"       name="qt_oce"   />
 	  <field field_ref="saltflx"      name="sfx"      />
 	  <field field_ref="taum"         name="taum"     />
+	  <field field_ref="utau"         name="tauuo"    />
+	  <field field_ref="vtau"         name="tauvo"    />
 	  <field field_ref="wspd"         name="windsp"   />
 	  <field field_ref="precip"       name="precip"   />
 	  <!-- ice and snow -->
@@ -44,7 +46,6 @@
 	  <field field_ref="e3u" />
 	  <field field_ref="ssu"          name="uos"      />
 	  <field field_ref="uoce"         name="uo"       operation="instant" freq_op="5d" > @uoce_e3u / @e3u </field>
-	  <field field_ref="utau"         name="tauuo"    />
           <field field_ref="uocetr_eff"   name="uocetr_eff"  />
           <!-- available with diaar5 -->
 	  <field field_ref="u_masstr"     name="vozomatr" />
@@ -56,7 +57,6 @@
 	  <field field_ref="e3v" />
 	  <field field_ref="ssv"          name="vos"      />
 	  <field field_ref="voce"         name="vo"       operation="instant" freq_op="5d" > @voce_e3v / @e3v </field>
-	  <field field_ref="vtau"         name="tauvo"    />
           <field field_ref="vocetr_eff"   name="vocetr_eff"  />
           <!-- available with diaar5 -->
 	  <field field_ref="v_masstr"     name="vomematr" />
diff --git a/tests/DONUT/EXPREF/file_def_nemo-oce.xml b/tests/DONUT/EXPREF/file_def_nemo-oce.xml
index 40a640ffbc616de642a88e5a7976cb0633b3a40d..0f5dc0b4e0ca47c227f1ed46d96122ad66f4481e 100644
--- a/tests/DONUT/EXPREF/file_def_nemo-oce.xml
+++ b/tests/DONUT/EXPREF/file_def_nemo-oce.xml
@@ -28,6 +28,8 @@
 	  <field field_ref="qt_oce"       name="qt_oce"   />
 	  <field field_ref="saltflx"      name="sfx"      />
 	  <field field_ref="taum"         name="taum"     />
+	  <field field_ref="utau"         name="tauuo"    />
+	  <field field_ref="vtau"         name="tauvo"    />
 	  <field field_ref="wspd"         name="windsp"   />
 	  <field field_ref="precip"       name="precip"   />
 	</file>
@@ -36,14 +38,12 @@
 	  <field field_ref="e3u" />
 	  <field field_ref="ssu"          name="uos"      />
 	  <field field_ref="uoce"         name="uo"       operation="instant" freq_op="1d" > @uoce_e3u / @e3u </field>
-	  <field field_ref="utau"         name="tauuo"    />
 	</file>
 	
 	<file id="file13" name_suffix="_grid_V" description="ocean V grid variables" >
 	  <field field_ref="e3v" />
 	  <field field_ref="ssv"          name="vos"      />
 	  <field field_ref="voce"         name="vo"       operation="instant" freq_op="1d" > @voce_e3v / @e3v </field>
-	  <field field_ref="vtau"         name="tauvo"    />
 	</file>
 	
       </file_group> 
diff --git a/tests/ICE_RHEO/EXPREF/file_def_nemo-oce.xml b/tests/ICE_RHEO/EXPREF/file_def_nemo-oce.xml
index 0dbbd2144297dcf224fc26e69fcf0c10ee7978bc..b307458f5d70ec49c7634a845887c617d487c101 100644
--- a/tests/ICE_RHEO/EXPREF/file_def_nemo-oce.xml
+++ b/tests/ICE_RHEO/EXPREF/file_def_nemo-oce.xml
@@ -13,14 +13,17 @@
 
       <file_group id="1ts" output_freq="1ts"  output_level="10" enabled=".TRUE."/> <!-- 1 time step files -->
       <file_group id="1h" output_freq="1h"  output_level="10" enabled=".TRUE."> <!-- 1h files -->
+        <file id="file1" name_suffix="_grid_T" description="ocean T grid variables" >
+          <field field_ref="utau"         name="sozotaux"  />
+          <field field_ref="vtau"         name="sometauy"  />
+        </file>
+
         <file id="file2" name_suffix="_grid_U" description="ocean U grid variables" >
           <field field_ref="uoce"         name="vozocrtx"  />
-          <field field_ref="utau"         name="sozotaux"  />
         </file>
 
         <file id="file3" name_suffix="_grid_V" description="ocean V grid variables" >
           <field field_ref="voce"         name="vomecrty"  />
-          <field field_ref="vtau"         name="sometauy"  />
         </file>
       </file_group>
 
diff --git a/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90 b/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90
index 8d5fe4eb615de8df2477cec3d0c65b87c3c64b02..2a22cfc80c11be0f67c616b8aa3e31cb72242066 100644
--- a/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90
+++ b/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90
@@ -125,37 +125,21 @@ CONTAINS
          windu(ji,jj) = Umax/sqrt(d*1000)*(d-2*mig(ji)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*min(kt*30./21600,1.)
          windv(ji,jj) = Umax/sqrt(d*1000)*(d-2*mjg(jj)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*Rwind*min(kt*30./21600,1.)
       END_2D
-      CALL lbc_lnk( 'usrdef_sbc', windu, 'U', -1., windv, 'V', -1. )
-
-      wndm_ice(:,:) = 0._wp      !!gm brutal....
 
       ! ------------------------------------------------------------ !
-      !    Wind module relative to the moving ice ( U10m - U_ice )   !
+      !    Wind module and stress relative to the moving ice ( U10m - U_ice )   !
       ! ------------------------------------------------------------ !
-      ! C-grid ice dynamics :   U & V-points (same as ocean)
       DO_2D( 0, 0, 0, 0 )
-         zwndi_t = (  windu(ji,jj) - r_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
+         zwndi_t = (  windu(ji,jj) - r_vfac * 0.5 * ( u_ice(ji-1,jj) + u_ice(ji,jj) )  )
          zwndj_t = (  windv(ji,jj) - r_vfac * 0.5 * ( v_ice(ji,jj-1) + v_ice(ji,jj) )  )
+         !
          wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
+         !
+         utau_ice(ji,jj) = zrhoa * Cd_atm * wndm_ice(ji,jj) * zwndi_t
+         vtau_ice(ji,jj) = zrhoa * Cd_atm * wndm_ice(ji,jj) * zwndj_t
       END_2D
-      CALL lbc_lnk( 'usrdef_sbc', wndm_ice, 'T',  1. )
-
-      !!gm brutal....
-      utau_ice  (:,:) = 0._wp
-      vtau_ice  (:,:) = 0._wp
-      !!gm end
+      CALL lbc_lnk( 'usrdef_sbc', wndm_ice, 'T',  1., utau_ice, 'T', -1., vtau_ice, 'T', -1. )
 
-      ! ------------------------------------------------------------ !
-      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
-      ! ------------------------------------------------------------ !
-      ! C-grid ice dynamics :   U & V-points (same as ocean)
-      DO_2D( 0, 0, 0, 0 )
-         utau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
-            &          * ( 0.5 * (windu(ji+1,jj) + windu(ji,jj) ) - r_vfac * u_ice(ji,jj) )
-         vtau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
-            &          * ( 0.5 * (windv(ji,jj+1) + windv(ji,jj) ) - r_vfac * v_ice(ji,jj) )
-      END_2D
-      CALL lbc_lnk( 'usrdef_sbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
       !
    END SUBROUTINE usrdef_sbc_ice_tau
 
diff --git a/tests/SWG/EXPREF/file_def_nemo-oce.xml b/tests/SWG/EXPREF/file_def_nemo-oce.xml
index fb16bbbd5abcd9873c2bbb7356a437604f2564bd..3b1c4d50b1da28e119c57f856dc623c5651651f5 100644
--- a/tests/SWG/EXPREF/file_def_nemo-oce.xml
+++ b/tests/SWG/EXPREF/file_def_nemo-oce.xml
@@ -27,6 +27,8 @@
           <field field_ref="sKE"          name="sKE"       operation = "instant" />
           <field field_ref="ssh"          name="sossheig"  operation = "instant" />
           <field field_ref="wetdep"       name="hswe_wd"   operation = "instant" />
+          <field field_ref="utau"         name="sozotaux"   operation = "instant"  />
+          <field field_ref="vtau"         name="sometauy"   operation = "instant"  />
         </file>
 
         <file id="file2" name_suffix="_grid_U" description="ocean U grid variables" >
@@ -34,7 +36,6 @@
           <field field_ref="e3u_0"        name="e3u_0"      operation = "instant" />
           <field field_ref="hu"           name="hu"         operation = "instant" />
           <field field_ref="ssu"          name="ssu"        operation = "instant"  />
-          <field field_ref="utau"         name="sozotaux"   operation = "instant"  />
         </file>
 
         <file id="file3" name_suffix="_grid_V" description="ocean V grid variables" >
@@ -42,7 +43,6 @@
           <field field_ref="e3v_0"        name="e3v_0"      operation = "instant" />
           <field field_ref="hv"           name="hv"         operation = "instant" />
           <field field_ref="ssv"          name="ssv"        operation = "instant"  />
-          <field field_ref="vtau"         name="sometauy"   operation = "instant"  />
         </file>
         
         <file id="file5" name_suffix="_grid_F" description="ocean F grid variables" >
diff --git a/tests/SWG/MY_SRC/usrdef_sbc.F90 b/tests/SWG/MY_SRC/usrdef_sbc.F90
index a3069954721d86aa464b7bb297e5d233f808c8b3..66ea77ad2ea3d1a8fc6f939af2f1336c644b7032 100644
--- a/tests/SWG/MY_SRC/usrdef_sbc.F90
+++ b/tests/SWG/MY_SRC/usrdef_sbc.F90
@@ -61,7 +61,7 @@ CONTAINS
       REAL(wp) ::   ztauu, ztauv          ! wind intensity projeted
       REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
       REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
-      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
+      REAL(wp) ::   zmod, zcoef           ! temporary variables
       !!---------------------------------------------------------------------
 
       ! ---------------------------- !
@@ -88,20 +88,17 @@ CONTAINS
       
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          ! length of the domain : 2000km x 2000km 
-         utau(ji,jj) = - ztauu * COS( rpi * gphiu(ji,jj) / 2000000._wp)
-         vtau(ji,jj) = - ztauv * COS( rpi * gphiv(ji,jj) / 2000000._wp)
+         utau(ji,jj) = - ztauu * COS( rpi * gphit(ji,jj) / 2000000._wp)
+         vtau(ji,jj) = - ztauv * COS( rpi * gphit(ji,jj) / 2000000._wp)
       END_2D
       
       ! module of wind stress and wind speed at T-point
       zcoef = 1. / ( zrhoa * zcdrag ) 
-      DO_2D( 0, 0, 0, 0 )
-         ztx = utau(ji-1,jj  ) + utau(ji,jj) 
-         zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
-         zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         zmod = SQRT( utau(ji,jj) * utau(ji,jj) + vtau(ji,jj) * vtau(ji,jj) )
          taum(ji,jj) = zmod
          wndm(ji,jj) = SQRT( zmod * zcoef )
       END_2D
-      CALL lbc_lnk( 'usrdef_sbc', taum(:,:), 'T', 1. , wndm(:,:), 'T', 1. )
       !
    END SUBROUTINE usrdef_sbc_oce
 
diff --git a/tests/WAD/EXPREF/file_def_nemo-oce.xml b/tests/WAD/EXPREF/file_def_nemo-oce.xml
index a9c10b7e0ec45084087a8614c64a3d13fc5392a0..90124a63601afba0a6d011cf701d31ec48cfa39d 100644
--- a/tests/WAD/EXPREF/file_def_nemo-oce.xml
+++ b/tests/WAD/EXPREF/file_def_nemo-oce.xml
@@ -69,16 +69,16 @@
 	  <field field_ref="qt"           name="sohefldo"  />
 	  <field field_ref="mldr10_1"     name="somxl010"  />
 	  <field field_ref="mldkz5"       name="somixhgt"  />
+          <field field_ref="utau"         name="sozotaux"  />
+          <field field_ref="vtau"         name="sometauy"  /> 
     </file>
 	
 	<file id="file2" name_suffix="_grid_U" description="ocean U grid variables" >
           <field field_ref="uoce"         name="vozocrtx"  />
-          <field field_ref="utau"         name="sozotaux"  />
     </file>
 	
 	<file id="file3" name_suffix="_grid_V" description="ocean V grid variables" >
           <field field_ref="voce"         name="vomecrty"  /> 
-          <field field_ref="vtau"         name="sometauy"  /> 
     </file>
 	
 	<file id="file4" name_suffix="_grid_W" description="ocean W grid variables" >