Index: NEMO/trunk/tests/DIA_GPU/EXPREF/axis_def_nemo.xml
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/axis_def_nemo.xml	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/axis_def_nemo.xml	(revision 14854)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/axis_def_nemo.xml
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/context_nemo.xml
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/context_nemo.xml	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/context_nemo.xml	(revision 14854)
@@ -0,0 +1,42 @@
+<!--
+ ============================================================================================== 
+    NEMO context
+============================================================================================== 
+-->
+<context id="nemo">
+    <!-- $id$ -->
+    <variable_definition>
+       <!-- Year/Month/Day of time origin for NetCDF files; defaults to 1800-01-01 -->
+       <variable id="ref_year"  type="int"> 1900 </variable>
+       <variable id="ref_month" type="int"> 01 </variable>
+       <variable id="ref_day"   type="int"> 01 </variable>
+       <variable id="rho0"      type="float" > 1026.0 </variable>
+       <variable id="cpocean"   type="float" > 3991.86795711963 </variable>
+       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable>
+       <variable id="rhoic"     type="float" > 917.0 </variable>
+       <variable id="rhosn"     type="float" > 330.0 </variable>
+       <variable id="missval"   type="float" > 1.e20 </variable>
+    </variable_definition>
+
+<!-- Fields definition -->
+    <field_definition src="./field_def_nemo-oce.xml"/>    <!--  NEMO ocean dynamics     -->
+    <field_definition src="./field_def_nemo-ice.xml"/>    <!--  NEMO sea-ice model      -->
+    <field_definition src="./field_def_nemo-pisces.xml"/> <!--  NEMO ocean biology      -->
+
+<!-- Files definition -->
+    <file_definition src="./file_def_nemo-oce.xml"/>     <!--  NEMO ocean dynamics      -->
+    <file_definition src="./file_def_nemo-ice.xml"/>     <!--  NEMO sea-ice model       -->
+    <file_definition src="./file_def_nemo-pisces.xml"/>  <!--  NEMO ocean biology       -->
+
+
+<!-- Axis definition -->
+    <axis_definition src="./axis_def_nemo.xml"/>
+ 
+<!-- Domain definition -->
+    <domain_definition src="./domain_def_nemo.xml"/>
+
+<!-- Grids definition -->
+    <grid_definition   src="./grid_def_nemo.xml"/>
+  
+
+</context>
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/domain_def_nemo.xml
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/domain_def_nemo.xml	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/domain_def_nemo.xml	(revision 14854)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/domain_def_nemo.xml
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/field_def_nemo-ice.xml
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/field_def_nemo-ice.xml	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/field_def_nemo-ice.xml	(revision 14854)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/field_def_nemo-ice.xml
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/field_def_nemo-innerttrc.xml
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/field_def_nemo-innerttrc.xml	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/field_def_nemo-innerttrc.xml	(revision 14854)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/field_def_nemo-innerttrc.xml
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/field_def_nemo-oce.xml
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/field_def_nemo-oce.xml	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/field_def_nemo-oce.xml	(revision 14854)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/field_def_nemo-oce.xml
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/field_def_nemo-pisces.xml
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/field_def_nemo-pisces.xml	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/field_def_nemo-pisces.xml	(revision 14854)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/field_def_nemo-pisces.xml
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/file_def_nemo-ice.xml
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/file_def_nemo-ice.xml	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/file_def_nemo-ice.xml	(revision 14854)
@@ -0,0 +1,155 @@
+<?xml version="1.0"?>
+    <!-- $id$ -->
+    
+ <!--
+============================================================================================================
+=                                           output files definition                                        =
+=                                      Define your own files for sea ice                                   =
+=                                         put the variables you want...                                    =
+============================================================================================================
+ -->
+    
+ <file_definition type="one_file" name="@expname@_@freq@_@startdate@_@enddate@" sync_freq="1mo" min_digits="4">
+   
+   <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE.">  <!-- 5d files -->   	
+     <file id="file21" name_suffix="_icemod" description="ice variables" enabled=".true." >
+
+       <!-- ice mask -->
+       <field field_ref="icemask"          name="simsk"   />
+       <field field_ref="icemask05"        name="simsk05" />
+       <field field_ref="icemask15"        name="simsk15" />
+       
+       <!-- general -->
+       <field field_ref="snwvolu"          name="snvolu" />
+       <field field_ref="icethic"          name="sithic" />
+       <field field_ref="icevolu"          name="sivolu" />
+       <field field_ref="iceconc"          name="siconc" />
+       <field field_ref="icesalt"          name="sisali" />
+       <field field_ref="iceapnd"          name="siapnd" />
+       <field field_ref="icevpnd"          name="sivpnd" />
+       <field field_ref="sst_m"            name="sst_m"  />
+       <field field_ref="sss_m"            name="sss_m"  />
+       
+       <!-- heat -->
+       <field field_ref="icetemp"          name="sitemp" />
+       <field field_ref="snwtemp"          name="sntemp" />
+       <field field_ref="icettop"          name="sittop" />
+       <field field_ref="icetbot"          name="sitbot" />
+       <field field_ref="icetsni"          name="sitsni" />
+       
+       <!-- momentum -->
+       <field field_ref="uice"             name="sivelu"  />
+       <field field_ref="vice"             name="sivelv"  />
+       <field field_ref="icevel"           name="sivelo"  />
+       <field field_ref="utau_ai"          name="utau_ai" />
+       <field field_ref="vtau_ai"          name="vtau_ai" />
+       <field field_ref="utau_oi"          name="utau_oi" />
+       <field field_ref="vtau_oi"          name="vtau_oi" />
+       
+       <!-- rheology -->
+       <field field_ref="icediv"           name="sidive"  />
+       <field field_ref="iceshe"           name="sishea"  />
+       <field field_ref="icestr"           name="sistre"  />
+       <field field_ref="normstr"          name="normstr" />
+       <field field_ref="sheastr"          name="sheastr" />
+       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
+       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
+       
+       <!-- heat fluxes -->
+       <field field_ref="qt_oce_ai"        name="qt_oce_ai"  />
+       <field field_ref="qt_atm_oi"        name="qt_atm_oi"  />
+       <field field_ref="qtr_ice_top"      name="qtr_ice_top"/>
+       <field field_ref="qtr_ice_bot"      name="qtr_ice_bot"/>
+       <field field_ref="qt_ice"           name="qt_ice"     />
+       <field field_ref="qsr_ice"          name="qsr_ice"    />
+       <field field_ref="qns_ice"          name="qns_ice"    />
+       <field field_ref="qemp_ice"         name="qemp_ice"   />
+       <field field_ref="albedo"           name="albedo"     />
+       
+       <field field_ref="hfxcndtop"        name="hfxcndtop"  />
+       <field field_ref="hfxcndbot"        name="hfxcndbot"  />
+       <field field_ref="hfxsensib"        name="hfxsensib"  />
+       
+       <!-- salt fluxes -->
+       <field field_ref="sfxice"           name="sfxice" />
+       
+       <!-- mass fluxes -->
+       <field field_ref="vfxice"           name="vfxice" />
+       <field field_ref="vfxsnw"           name="vfxsnw" />
+
+       <!-- categories -->
+       <field field_ref="icemask_cat"      name="simskcat"/>
+       <field field_ref="snwthic_cat"      name="snthicat"/>
+       <field field_ref="iceconc_cat"      name="siconcat"/>
+       <field field_ref="icethic_cat"      name="sithicat"/>
+       <field field_ref="icesalt_cat"      name="sisalcat"/>
+       <field field_ref="icetemp_cat"      name="sitemcat"/>
+       <field field_ref="snwtemp_cat"      name="sntemcat"/>
+       <field field_ref="icettop_cat"      name="sitopcat"/>
+
+     </file>
+     
+     <file id="file22" name_suffix="_SBC_scalar" description="scalar variables" enabled=".true." >
+       <!-- global contents -->
+       <field field_ref="ibgvol_tot"     name="ibgvol_tot"   />
+       <field field_ref="sbgvol_tot"     name="sbgvol_tot"   />
+       <field field_ref="ibgarea_tot"    name="ibgarea_tot"  />
+       <field field_ref="ibgsalt_tot"    name="ibgsalt_tot"  />
+       <field field_ref="ibgheat_tot"    name="ibgheat_tot"  />
+       <field field_ref="sbgheat_tot"    name="sbgheat_tot"  />
+       
+       <!-- global drifts (conservation checks) -->
+       <field field_ref="ibgvolume"      name="ibgvolume"    />
+       <field field_ref="ibgsaltco"      name="ibgsaltco"    />
+       <field field_ref="ibgheatco"      name="ibgheatco"    />
+       <field field_ref="ibgheatfx"      name="ibgheatfx"    />
+       
+       <!-- global forcings  -->
+       <field field_ref="ibgfrcvoltop"   name="ibgfrcvoltop" />
+       <field field_ref="ibgfrcvolbot"   name="ibgfrcvolbot" />
+       <field field_ref="ibgfrctemtop"   name="ibgfrctemtop" />
+       <field field_ref="ibgfrctembot"   name="ibgfrctembot" />
+       <field field_ref="ibgfrcsal"      name="ibgfrcsal"    />
+       <field field_ref="ibgfrchfxtop"   name="ibgfrchfxtop" />
+       <field field_ref="ibgfrchfxbot"   name="ibgfrchfxbot" />
+     </file>
+     
+   </file_group>
+   
+   <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_group id="2h"  output_freq="2h"  output_level="10" enabled=".TRUE."/> <!-- 2h files -->
+   <file_group id="3h"  output_freq="3h"  output_level="10" enabled=".TRUE."/> <!-- 3h files -->     
+   <file_group id="4h"  output_freq="4h"  output_level="10" enabled=".TRUE."/> <!-- 4h files -->
+   <file_group id="6h"  output_freq="6h"  output_level="10" enabled=".TRUE."/> <!-- 6h files -->      
+
+   <file_group id="1m"  output_freq="1mo" output_level="10" enabled=".TRUE."> <!-- real monthly files -->
+
+     <!-- To compute transport through straits : need to read ice mask at ice iteration at freq_offset = 1mo - nn_fsbc
+      <file id="file23" name_suffix="_strait_ice" description="transport variables through straits" >
+          <field field_ref="strait_mifl"      name="simassacrossline"  freq_offset="1mo-4ts" />
+          <field field_ref="strait_msfl"      name="snmassacrossline"  freq_offset="1mo-4ts" />
+          <field field_ref="strait_arfl"      name="siareaacrossline"  freq_offset="1mo-4ts" />
+        </file>
+     -->
+
+    </file_group>
+
+   <file_group id="2m"  output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files -->
+   <file_group id="3m"  output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files -->
+   <file_group id="4m"  output_freq="4mo" output_level="10" enabled=".TRUE."/> <!-- real 4m files -->
+   <file_group id="6m"  output_freq="6mo" output_level="10" enabled=".TRUE."/> <!-- real 6m files -->
+   <file_group id="1y"  output_freq="1y"  output_level="10" enabled=".TRUE."/> <!-- real yearly files -->
+   <file_group id="2y"  output_freq="2y"  output_level="10" enabled=".TRUE."/> <!-- real 2y files -->
+   <file_group id="5y"  output_freq="5y"  output_level="10" enabled=".TRUE."/> <!-- real 5y files -->
+   <file_group id="10y" output_freq="10y" output_level="10" enabled=".TRUE."/> <!-- real 10y files -->
+ 
+   <!-- To compute transport through straits : need to read ice mask at ice iteration at freq_offset = - nn_fsbc + 1
+   <file id="maskMFO" name="maskMFO" enabled="true" mode="read" output_freq="1mo" cyclic="true"  >
+      <field id="maskMFO_u_ice" operation="instant"  freq_offset="-3ts"  grid_ref="grid_U_4strait_ice" />
+      <field id="maskMFO_v_ice" operation="instant"  freq_offset="-3ts"  grid_ref="grid_V_4strait_ice" />
+   </file>
+
+     -->
+ </file_definition>
+ 
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/file_def_nemo-oce.xml
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/file_def_nemo-oce.xml	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/file_def_nemo-oce.xml	(revision 14854)
@@ -0,0 +1,173 @@
+<?xml version="1.0"?>
+
+<!--
+============================================================================================================
+=                                           output files definition                                        =
+=                                            Define your own filesfor ocean dynamics context               =
+=                                         put the variables you want...                                    =
+============================================================================================================
+    -->
+   
+    <file_definition type="one_file" name="@expname@_@freq@_@startdate@_@enddate@" sync_freq="1mo" min_digits="4">
+    
+      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE.">  <!-- 5d files -->   
+        <file id="file11" name_suffix="_grid_T" description="ocean T grid variables" >
+          <field field_ref="e3t"      />
+          <field field_ref="toce"         name="thetao"                                                                      operation="instant" freq_op="5d" > @toce_e3t / @e3t </field>
+          <field field_ref="soce"         name="so"                                                                          operation="instant" freq_op="5d" > @soce_e3t / @e3t </field>
+          <field field_ref="sst"          name="tos"   />
+          <field field_ref="sss"          name="sos"   />
+          <field field_ref="ssh"          name="zos"   />
+          <field field_ref="sst"          name="tosstd"       long_name="sea surface temperature standard deviation"         operation="instant" freq_op="5d" > sqrt( @sst2 - @sst * @sst ) </field>
+          <field field_ref="ssh"          name="zosstd"       long_name="sea surface height above geoid standard deviation"  operation="instant" freq_op="5d" > sqrt( @ssh2 - @ssh * @ssh ) </field>
+          <field field_ref="sst"          name="sstdcy"       long_name="amplitude of sst diurnal cycle"                     operation="average" freq_op="1d" > @sstmax - @sstmin </field>
+          <field field_ref="mldkz5"   />
+          <field field_ref="mldr10_1" />
+          <field field_ref="mldr10_1"     name="mldr10_1dcy"  long_name="amplitude of mldr10_1 diurnal cycle"                operation="average" freq_op="1d" > @mldr10_1max - @mldr10_1min </field>
+          <field field_ref="sbt"      />
+          <field field_ref="heatc"        name="heatc" />
+          <field field_ref="saltc"        name="saltc" />
+
+	  <field field_ref="empmr"        name="wfo"      />
+	  <field field_ref="qsr_oce"      name="qsr_oce"  />
+	  <field field_ref="qns_oce"      name="qns_oce"  />
+	  <field field_ref="qt_oce"       name="qt_oce"   />
+	  <field field_ref="saltflx"      name="sfx"      />
+	  <field field_ref="taum"         name="taum"     />
+	  <field field_ref="wspd"         name="windsp"   />
+	  <field field_ref="precip"       name="precip"   />
+	  <!-- ice and snow -->
+	  <field field_ref="snowpre" />
+	</file>
+
+	<file id="file12" name_suffix="_grid_U" description="ocean U grid variables" >
+	  <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" />
+	  <field field_ref="u_heattr"     name="sozohetr" />
+	  <field field_ref="u_salttr"     name="sozosatr" />
+	</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"    />
+          <field field_ref="vocetr_eff"   name="vocetr_eff"  />
+          <!-- available with diaar5 -->
+	  <field field_ref="v_masstr"     name="vomematr" />
+	  <field field_ref="v_heattr"     name="somehetr" />
+          <field field_ref="v_salttr"     name="somesatr" />
+	</file>
+	
+	<file id="file14" name_suffix="_grid_W" description="ocean W grid variables" >
+	  <field field_ref="e3w" />
+	  <field field_ref="woce"         name="wo"       />
+	  <field field_ref="avt"          name="difvho"   />
+	  <field field_ref="w_masstr"     name="vovematr" />
+         <!-- variables available with ln_zdfiwm =T -->
+          <field field_ref="av_wave"      name="av_wave"    />
+          <field field_ref="bn2"          name="bn2"        />
+          <field field_ref="bflx_iwm"     name="bflx_iwm"   />
+          <field field_ref="pcmap_iwm"    name="pcmap_iwm"  />
+          <field field_ref="emix_iwm"     name="emix_iwm"   />
+          <field field_ref="av_ratio"     name="av_ratio"   />
+	</file>
+
+        <file id="file15" name_suffix="_scalar" description="scalar variables" >
+          <!-- global drifts (conservation checks) -->
+          <field field_ref="bgtemper"     name="bgtemper"    />
+          <field field_ref="bgsaline"     name="bgsaline"    />
+          <field field_ref="bgheatco"     name="bgheatco"    />
+          <field field_ref="bgheatfx"     name="bgheatfx"    />
+          <field field_ref="bgsaltco"     name="bgsaltco"    />
+          <field field_ref="bgvolssh"     name="bgvolssh"    />
+          <field field_ref="bgvole3t"     name="bgvole3t"    />
+
+          <!-- global surface forcings  -->
+          <field field_ref="bgfrcvol"     name="bgfrcvol"    />
+          <field field_ref="bgfrctem"     name="bgfrctem"    />
+          <field field_ref="bgfrchfx"     name="bgfrchfx"    />
+	  <field field_ref="bgfrcsal"     name="bgfrcsal"    />
+
+	  <field field_ref="masstot"      name="masso"  />
+          <field field_ref="voltot"       name="volo"  />
+          <field field_ref="sshthster"    name="zostoga"  />
+          <field field_ref="temptot"      name="bigthetaoga"  />
+          <field field_ref="saltot"       name="soga"  />
+          <field field_ref="ssttot"       name="tosga"  />
+        </file>
+
+      </file_group> 
+
+
+      <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."> <!-- real monthly files -->
+
+        <file id="file16" name_suffix="_diaptr2D" description="zonal mean variables" >
+          <field field_ref="sophtove"    name="htovovrt"      grid_ref="grid_ptr_T_2D"  />
+          <field field_ref="sopstove"    name="sltovovrt"     grid_ref="grid_ptr_T_2D"  />
+          <field field_ref="sophtgyre"   name="htovgyre"      grid_ref="grid_ptr_T_2D"  />
+          <field field_ref="sopstgyre"   name="sltogyre"      grid_ref="grid_ptr_T_2D"  />
+          <field field_ref="sophtbtr"    name="htbtr"         grid_ref="grid_ptr_T_2D"  />
+          <field field_ref="sopstbtr"    name="sltbtr"        grid_ref="grid_ptr_T_2D"  />
+          <field field_ref="sophtadv"    name="htadv"         grid_ref="grid_ptr_T_2D"  />
+          <field field_ref="sopstadv"    name="sltadv"        grid_ref="grid_ptr_T_2D"  />
+          <field field_ref="sophtldf"    name="htldf"         grid_ref="grid_ptr_T_2D"  />
+          <field field_ref="sopstldf"    name="sltldf"        grid_ref="grid_ptr_T_2D"  />
+          <field field_ref="sophtvtr"    name="hfbasin"       grid_ref="grid_ptr_T_2D"  />
+          <field field_ref="sopstvtr"    name="sltbasin"      grid_ref="grid_ptr_T_2D"  />
+          <field field_ref="sophteiv"    name="hfbasinpmadv"  grid_ref="grid_ptr_T_2D"  />
+          <field field_ref="sopsteiv"    name="sltbasinpmadv" grid_ref="grid_ptr_T_2D"  />
+        </file>
+
+
+         <file id="file17" name_suffix="_diaptr3D" description="zonal mean variables" >
+          <field field_ref="zomsf"       name="msftyz"        grid_ref="grid_ptr_W_3D"  />
+          <field field_ref="zotem"       name="znltem"        grid_ref="grid_ptr_T_3D"  />
+          <field field_ref="zosal"       name="znlsal"        grid_ref="grid_ptr_T_3D"  />
+          <field field_ref="zosrf"       name="znlsrf"        grid_ref="grid_ptr_T_3D"  />
+        </file>
+
+      <!--  
+         <file id="file18" name_suffix="_strait_oce" description="transport variables through straits" >
+           <field field_ref="masstr_strait"        name="mfo"               />
+         </file>	   
+      -->
+
+      </file_group>
+
+      <file_group id="1y"  output_freq="1y" output_level="10" enabled=".TRUE."> <!-- real yearly files -->
+      </file_group>
+
+      <!--  To compute transport through straits : need to read mask file ( every month is the best otherwise costly )
+
+      <file id="maskMFO"  name="maskMFO" enabled="true" mode="read" output_freq="1mo" cyclic="true"  >
+        <field id="maskMFO_u" operation="instant" freq_offset="1mo" grid_ref="grid_U_4strait" />
+        <field id="maskMFO_v" operation="instant" freq_offset="1mo" grid_ref="grid_V_4strait" />
+      </file>
+
+    -->
+
+      
+      <file_group id="1ts" output_freq="1ts" output_level="10" enabled=".TRUE."/> <!-- 1 time step files -->
+      <file_group id="2h"  output_freq="2h"  output_level="10" enabled=".TRUE."/> <!-- 2h files -->
+      <file_group id="3h"  output_freq="3h"  output_level="10" enabled=".TRUE."/> <!-- 3h files -->     
+      <file_group id="4h"  output_freq="4h"  output_level="10" enabled=".TRUE."/> <!-- 4h files -->
+      <file_group id="6h"  output_freq="6h"  output_level="10" enabled=".TRUE."/> <!-- 6h files -->     
+      <file_group id="2m"  output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files -->
+      <file_group id="3m"  output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files -->
+      <file_group id="4m"  output_freq="4mo" output_level="10" enabled=".TRUE."/> <!-- real 4m files -->
+      <file_group id="6m"  output_freq="6mo" output_level="10" enabled=".TRUE."/> <!-- real 6m files -->
+      <file_group id="2y"  output_freq="2y"  output_level="10" enabled=".TRUE."/> <!-- real 2y files -->
+      <file_group id="5y"  output_freq="5y"  output_level="10" enabled=".TRUE."/> <!-- real 5y files -->
+      <file_group id="10y" output_freq="10y" output_level="10" enabled=".TRUE."/> <!-- real 10y files -->
+
+   </file_definition>
+    
+
+
+
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/file_def_nemo-pisces.xml
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/file_def_nemo-pisces.xml	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/file_def_nemo-pisces.xml	(revision 14854)
@@ -0,0 +1,128 @@
+<?xml version="1.0"?>
+
+<!--
+============================================================================================================
+=                                           output files definition                                        =
+=                                            Define your own filesfor ocean dynamics context               =
+=                                         put the variables you want...                                    =
+============================================================================================================
+    -->
+   
+    <file_definition type="one_file" name="@expname@_@freq@_@startdate@_@enddate@" sync_freq="1mo" min_digits="4">
+    
+      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE.">  <!-- 5d files -->   
+
+        <file id="file31" name_suffix="_bioscalar" description="pisces sms variables" >
+          <field field_ref="tdenit"   name="tdenit"    unit="TgN/yr" operation="instant" > tdenit * 14. * 86400. * 365. / 1e12 </field>
+          <field field_ref="tnfix"    name="tnfix"     unit="TgN/yr" operation="instant" > tnfix * 14. * 86400. * 365. / 1e12 </field>
+          <field field_ref="tcflx"    name="tcflx"     unit="PgC/yr" operation="instant" > tcflx * -1. * 12. * 86400. * 365. / 1e15 </field>
+          <field field_ref="tcflxcum" name="tcflxcum"  unit="PgC"    operation="instant" > tcflxcum * -1. * 12. / 1e15 </field>
+          <field field_ref="tcexp"    name="tcexp"     unit="PgC/yr" operation="instant" > tcexp * 12. * 86400. * 365. / 1e15 </field>
+          <field field_ref="tintpp"   name="tintpp"    unit="PgC/yr" operation="instant" > tintpp * 12. * 86400. * 365. / 1e15 </field>
+          <field field_ref="pno3tot"  name="pno3tot"   unit="umolN"  > pno3tot * 16. / 122. * 1e6 </field>
+          <field field_ref="ppo4tot"  name="ppo4tot"   unit="umolP"  > ppo4tot * 1. / 122. * 1e6 </field>
+          <field field_ref="psiltot"  name="psiltot"   unit="umolC"  > psiltot * 1e6  </field>
+          <field field_ref="palktot"  name="palktot"   unit="umolC"  > palktot * 1e6  </field>
+          <field field_ref="pfertot"  name="pfertot"   unit="nmolFe" > pfertot * 1e9  </field>
+        </file>
+
+      </file_group> 
+
+
+      <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."> <!-- real monthly files -->
+
+        <file id="file32" name_suffix="_ptrc_T" description="pisces sms variables" >
+          <field field_ref="DIC"       name="DIC"      operation="average" freq_op="1mo" > @DIC_e3t / @e3t </field>
+          <field field_ref="Alkalini"  name="Alkalini" operation="average" freq_op="1mo" > @Alkalini_e3t / @e3t </field>
+          <field field_ref="O2"        name="O2"       operation="average" freq_op="1mo" > @O2_e3t / @e3t </field>
+          <field field_ref="PO4"       name="PO4"      operation="average" freq_op="1mo" > @PO4_e3t / @e3t </field>
+          <field field_ref="Si"        name="Si"       operation="average" freq_op="1mo" > @Si_e3t / @e3t </field>
+          <field field_ref="Fer"       name="Fer"      operation="average" freq_op="1mo" > @Fer_e3t / @e3t </field>
+          <field field_ref="NCHL"      name="NCHL"     operation="average" freq_op="1mo" > @NCHL_e3t / @e3t </field>
+          <field field_ref="DCHL"      name="DCHL"     operation="average" freq_op="1mo" > @DCHL_e3t / @e3t </field>
+          <field field_ref="NO3"       name="NO3"      operation="average" freq_op="1mo" > @NO3_e3t / @e3t </field>
+        </file>
+
+        <file id="file33" name_suffix="_diad_T" description="additional pisces diagnostics" >
+          <field field_ref="Cflx"     />
+          <field field_ref="Dpco2"    />
+        </file>
+
+      </file_group>
+
+      <file_group id="1y"  output_freq="1y" output_level="10" enabled=".TRUE."> <!-- real yearly files -->
+
+	<file id="file34" name_suffix="_ptrc_T" description="pisces sms variables" >
+          <field field_ref="DIC"       name="DIC"      operation="average" freq_op="1y" > @DIC_e3t / @e3t </field>
+          <field field_ref="Alkalini"  name="Alkalini" operation="average" freq_op="1y" > @Alkalini_e3t / @e3t </field>
+          <field field_ref="O2"        name="O2"       operation="average" freq_op="1y" > @O2_e3t / @e3t </field>
+          <field field_ref="CaCO3"     name="CaCO3"    operation="average" freq_op="1y" > @CaCO3_e3t / @e3t </field>
+          <field field_ref="PO4"       name="PO4"      operation="average" freq_op="1y" > @PO4_e3t / @e3t </field>
+          <field field_ref="POC"       name="POC"      operation="average" freq_op="1y" > @POC_e3t / @e3t </field>
+          <field field_ref="Si"        name="Si"       operation="average" freq_op="1y" > @Si_e3t / @e3t </field>
+          <field field_ref="PHY"       name="PHY"      operation="average" freq_op="1y" > @PHY_e3t / @e3t </field>
+          <field field_ref="ZOO"       name="ZOO"      operation="average" freq_op="1y" > @ZOO_e3t / @e3t </field>
+          <field field_ref="DOC"       name="DOC"      operation="average" freq_op="1y" > @DOC_e3t / @e3t </field>
+          <field field_ref="PHY2"      name="PHY2"     operation="average" freq_op="1y" > @PHY2_e3t / @e3t </field>
+          <field field_ref="ZOO2"      name="ZOO2"     operation="average" freq_op="1y" > @ZOO2_e3t / @e3t </field>
+          <field field_ref="DSi"       name="DSi"      operation="average" freq_op="1y" > @DSi_e3t / @e3t </field>
+          <field field_ref="Fer"       name="Fer"      operation="average" freq_op="1y" > @Fer_e3t / @e3t </field>
+          <field field_ref="BFe"       name="BFe"      operation="average" freq_op="1y" > @BFe_e3t / @e3t </field>
+          <field field_ref="GOC"       name="GOC"      operation="average" freq_op="1y" > @GOC_e3t / @e3t </field>
+          <field field_ref="SFe"       name="SFe"      operation="average" freq_op="1y" > @SFe_e3t / @e3t </field>
+          <field field_ref="DFe"       name="DFe"      operation="average" freq_op="1y" > @DFe_e3t / @e3t </field>
+          <field field_ref="GSi"       name="GSi"      operation="average" freq_op="1y" > @GSi_e3t / @e3t </field>
+          <field field_ref="NFe"       name="NFe"      operation="average" freq_op="1y" > @NFe_e3t / @e3t </field>
+          <field field_ref="NCHL"      name="NCHL"     operation="average" freq_op="1y" > @NCHL_e3t / @e3t </field>
+          <field field_ref="DCHL"      name="DCHL"     operation="average" freq_op="1y" > @DCHL_e3t / @e3t </field>
+          <field field_ref="NO3"       name="NO3"      operation="average" freq_op="1y" > @NO3_e3t / @e3t </field>
+          <field field_ref="NH4"       name="NH4"      operation="average" freq_op="1y" > @NH4_e3t / @e3t </field>
+        </file>
+
+        <file id="file35" name_suffix="_diad_T" description="additional pisces diagnostics" >
+          <field field_ref="PH"       />
+          <field field_ref="PAR"      />
+          <field field_ref="PPPHYN"    />
+          <field field_ref="PPPHYD"   />
+          <field field_ref="PPPHYP"   />
+          <field field_ref="PPNEWN"   />
+          <field field_ref="PPNEWD"   />
+          <field field_ref="PPNEWP"   />
+          <field field_ref="PBSi"     />
+          <field field_ref="PFeN"     />
+          <field field_ref="PFeD"     />
+          <field field_ref="PCAL"     />
+          <field field_ref="DCAL"     />
+          <field field_ref="GRAZ1"    />
+          <field field_ref="GRAZ2"    />
+          <field field_ref="EPC100"   />
+          <field field_ref="EPFE100"  />
+          <field field_ref="EPSI100"  />
+          <field field_ref="EPCAL100" />
+          <field field_ref="Cflx"     />
+          <field field_ref="Dpco2"    />
+          <field field_ref="Heup"     />
+          <field field_ref="Nfix"     />
+        </file>
+     </file_group>
+	
+	
+      
+      <file_group id="1ts" output_freq="1ts" output_level="10" enabled=".TRUE."/> <!-- 1 time step files -->
+      <file_group id="2h"  output_freq="2h"  output_level="10" enabled=".TRUE."/> <!-- 2h files -->
+      <file_group id="3h"  output_freq="3h"  output_level="10" enabled=".TRUE."/> <!-- 3h files -->     
+      <file_group id="4h"  output_freq="4h"  output_level="10" enabled=".TRUE."/> <!-- 4h files -->
+      <file_group id="6h"  output_freq="6h"  output_level="10" enabled=".TRUE."/> <!-- 6h files -->     
+      <file_group id="2m"  output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files -->
+      <file_group id="3m"  output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files -->
+      <file_group id="4m"  output_freq="4mo" output_level="10" enabled=".TRUE."/> <!-- real 4m files -->
+      <file_group id="6m"  output_freq="6mo" output_level="10" enabled=".TRUE."/> <!-- real 6m files -->
+      <file_group id="2y"  output_freq="2y"  output_level="10" enabled=".TRUE."/> <!-- real 2y files -->
+      <file_group id="5y"  output_freq="5y"  output_level="10" enabled=".TRUE."/> <!-- real 5y files -->
+      <file_group id="10y" output_freq="10y" output_level="10" enabled=".TRUE."/> <!-- real 10y files -->
+
+   </file_definition>
+    
+
+
+
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/grid_def_nemo.xml
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/grid_def_nemo.xml	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/grid_def_nemo.xml	(revision 14854)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/grid_def_nemo.xml
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/iodef.xml
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/iodef.xml	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/iodef.xml	(revision 14854)
@@ -0,0 +1,26 @@
+<?xml version="1.0"?>
+<simulation>
+
+<!-- ============================================================================================ -->
+<!-- XIOS context                                                                                 -->
+<!-- ============================================================================================ -->
+
+  <context id="xios" >
+
+      <variable_definition>
+
+          <variable id="info_level"                type="int">10</variable>
+          <variable id="using_server"              type="bool">false</variable>
+          <variable id="using_oasis"               type="bool">false</variable>
+          <variable id="oasis_codes_id"            type="string" >oceanx</variable>
+
+      </variable_definition>
+  </context>
+
+<!-- ============================================================================================ -->
+<!-- NEMO  CONTEXT add and suppress the components you need                                       -->
+<!-- ============================================================================================ -->
+
+  <context id="nemo" src="./context_nemo.xml"/>       <!--  NEMO       -->
+
+</simulation>
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_cfg
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_cfg	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_cfg	(revision 14854)
@@ -0,0 +1,464 @@
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! NEMO/OCE  Configuration namelist : overwrite default values defined in SHARED/namelist_ref
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!!             ORCA2 - ICE - PISCES configuration                     !!
+!!======================================================================
+!!              ***  Domain & Run management namelists  ***           !!
+!!                                                                    !!
+!!   namrun       parameters of the run
+!!   namdom       space and time domain
+!!   namcfg       parameters of the configuration                       (default: user defined GYRE)
+!!   namwad       Wetting and drying                                    (default: OFF)
+!!   namtsd       data: temperature & salinity                          (default: OFF)
+!!   namcrs       coarsened grid (for outputs and/or TOP)               (ln_crs =T)
+!!   namc1d       1D configuration options                              ("key_c1d")
+!!   namc1d_dyndmp 1D newtonian damping applied on currents             ("key_c1d")
+!!   namc1d_uvd   1D data (currents)                                    ("key_c1d")
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namrun        !   parameters of the run
+!-----------------------------------------------------------------------
+   cn_exp      =  "ORCA2"  !  experience name
+   nn_it000    =       1   !  first time step
+   nn_itend    =    5840   !  last  time step (std 5475)
+   nn_istate   =       0   !  output the initial state (1) or not (0)
+/
+!-----------------------------------------------------------------------
+&namdom        !   time and space domain
+!-----------------------------------------------------------------------
+   rn_Dt       = 5400.     !  time step for the dynamics and tracer
+/
+!-----------------------------------------------------------------------
+&namcfg        !   parameters of the configuration                      (default: use namusr_def in namelist_cfg)
+!-----------------------------------------------------------------------
+   ln_read_cfg = .true.    !  (=T) read the domain configuration file
+      cn_domcfg = "ORCA_R2_zps_domcfg"    ! domain configuration filename
+      !
+      ln_closea    = .false.    !  F => suppress closed seas (defined by closea_mask field) 
+      !                         !       from the bathymetry at runtime.
+/
+!-----------------------------------------------------------------------
+&namtile        !   parameters of the tiling
+!-----------------------------------------------------------------------
+   ln_tile = .false.     !  Use tiling (T) or not (F)
+   nn_ltile_i = 10       !  Length of tiles in i
+   nn_ltile_j = 10       !  Length of tiles in j
+/
+!-----------------------------------------------------------------------
+&namtsd        !    Temperature & Salinity Data  (init/dmp)             (default: OFF)
+!-----------------------------------------------------------------------
+   !                       ! =T  read T-S fields for:
+   ln_tsd_init = .true.          !  ocean initialisation
+   ln_tsd_dmp  = .true.          !  T-S restoring   (see namtra_dmp)
+
+   cn_dir = './'     !  root directory for the T-S data location
+   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
+   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
+   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   !
+   sn_tem = 'data_1m_potential_temperature_nomask',     -1.     ,'votemper',   .true.    , .true. , 'yearly'  ,    ''    ,    ''    ,    ''
+   sn_sal = 'data_1m_salinity_nomask'             ,     -1.     ,'vosaline',   .true.    , .true. , 'yearly'  ,    ''    ,    ''    ,    ''
+/
+!!======================================================================
+!!            ***  Surface Boundary Condition namelists  ***          !!
+!!                                                                    !!
+!!   namsbc          surface boundary condition manager                 (default: NO selection)
+!!   namsbc_flx      flux               formulation                     (ln_flx     =T)
+!!   namsbc_blk      Bulk formulae formulation                          (ln_blk     =T)
+!!   namsbc_cpl      CouPLed            formulation                     ("key_oasis3" )
+!!   namsbc_sas      Stand-Alone Surface module                         (SAS_SRC  only)
+!!   namsbc_iif      Ice-IF: use observed ice cover                     (nn_ice = 1   )
+!!   namtra_qsr      penetrative solar radiation                        (ln_traqsr  =T)
+!!   namsbc_ssr      sea surface restoring term (for T and/or S)        (ln_ssr     =T)
+!!   namsbc_rnf      river runoffs                                      (ln_rnf     =T)
+!!   namsbc_apr      Atmospheric Pressure                               (ln_apr_dyn =T)
+!!   namsbc_isf      ice shelf melting/freezing                         (ln_isfcav  =T : read (ln_read_cfg=T) or set or usr_def_zgr )
+!!   namsbc_iscpl    coupling option between land ice model and ocean   (ln_isfcav  =T)
+!!   namsbc_wave     external fields from wave model                    (ln_wave    =T)
+!!   namberg         iceberg floats                                     (ln_icebergs=T)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namsbc        !   Surface Boundary Condition manager                   (default: NO selection)
+!-----------------------------------------------------------------------
+   nn_fsbc     = 4         !  frequency of SBC module call
+                           !     (also = the frequency of sea-ice & iceberg model call)
+                     ! Type of air-sea fluxes 
+   ln_blk      = .true.    !  Bulk formulation                          (T => fill namsbc_blk )
+                     ! Sea-ice :
+   nn_ice      = 2         !  =0 no ice boundary condition
+      !                    !  =1 use observed ice-cover                 (  => fill namsbc_iif )
+      !                    !  =2 or 3 for SI3 and CICE, respectively
+                     ! Misc. options of sbc : 
+   ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr)
+   ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr)
+   ln_rnf      = .true.    !  runoffs                                   (T => fill namsbc_rnf)
+   nn_fwb      = 2         !  FreshWater Budget: 
+   !                       !    =2 annual global mean of e-p-r set to zero
+   ln_wave     = .false.   !  Activate coupling with wave  (T => fill namsbc_wave)
+/
+!-----------------------------------------------------------------------
+&namsbc_blk    !   namsbc_blk  generic Bulk formula                     (ln_blk =T)
+!-----------------------------------------------------------------------
+   !                    !  bulk algorithm :
+   ln_NCAR    = .true.     ! "NCAR"      algorithm   (Large and Yeager 2008)
+
+   cn_dir = './'  !  root directory for the bulk data location
+   !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________!
+   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ !          weights filename            ! rotation ! land/sea mask !
+   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                                      ! pairing  !    filename   !
+   sn_wndi     = 'u_10.15JUNE2009_fill'       ,    6.        , 'U_10_MOD',   .false.   , .true. , 'yearly'  , 'weights_core2_orca2_bicub'  , 'Uwnd'   , ''
+   sn_wndj     = 'v_10.15JUNE2009_fill'       ,    6.        , 'V_10_MOD',   .false.   , .true. , 'yearly'  , 'weights_core2_orca2_bicub'  , 'Vwnd'   , ''
+   sn_qsr      = 'ncar_rad.15JUNE2009_fill'   ,   24.        , 'SWDN_MOD',   .false.   , .true. , 'yearly'  , 'weights_core2_orca2_bilin' , ''       , ''
+   sn_qlw      = 'ncar_rad.15JUNE2009_fill'   ,   24.        , 'LWDN_MOD',   .false.   , .true. , 'yearly'  , 'weights_core2_orca2_bilin' , ''       , ''
+   sn_tair     = 't_10.15JUNE2009_fill'       ,    6.        , 'T_10_MOD',   .false.   , .true. , 'yearly'  , 'weights_core2_orca2_bilin' , ''       , ''
+   sn_humi     = 'q_10.15JUNE2009_fill'       ,    6.        , 'Q_10_MOD',   .false.   , .true. , 'yearly'  , 'weights_core2_orca2_bilin' , ''       , ''
+   sn_prec     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'PRC_MOD1',   .false.   , .true. , 'yearly'  , 'weights_core2_orca2_bilin' , ''       , ''
+   sn_snow     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'SNOW'    ,   .false.   , .true. , 'yearly'  , 'weights_core2_orca2_bilin' , ''       , ''
+   sn_slp      = 'slp.15JUNE2009_fill'        ,    6.        , 'SLP'     ,   .false.   , .true. , 'yearly'  , 'weights_core2_orca2_bilin' , ''       , ''
+/
+!-----------------------------------------------------------------------
+&namsbc_abl    !   Atmospheric Boundary Layer formulation           (ln_abl = T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtra_qsr    !   penetrative solar radiation                          (ln_traqsr =T)
+!-----------------------------------------------------------------------
+   !                       ! type of penetration                        (default: NO selection)
+   ln_qsr_rgb  = .true.       !  RGB light penetration (Red-Green-Blue)
+   !
+   nn_chldta   =      1       !  RGB : Chl data (=1) or cst value (=0)
+
+   cn_dir = './'  !  root directory for the chlorophyl data location
+   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
+   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
+   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   !
+   sn_chl      ='chlorophyll'            ,        -1.        , 'CHLA'    ,   .true.    , .true. , 'yearly'  , ''               , ''       , ''
+/
+!-----------------------------------------------------------------------
+&namsbc_ssr    !   surface boundary condition : sea surface restoring   (ln_ssr =T)
+!-----------------------------------------------------------------------
+   nn_sssr     =     2     !  add a damping term to the surface freshwater flux
+      rn_deds     =  -166.67  !  magnitude of the damping on salinity   [mm/day]
+      ln_sssr_bnd =  .true.   !  flag to bound erp term (associated with nn_sssr=2)
+      rn_sssr_bnd =   4.e0    !  ABS(Max/Min) value of the damping erp term [mm/day]
+/
+!-----------------------------------------------------------------------
+&namsbc_rnf    !   runoffs                                              (ln_rnf =T)
+!-----------------------------------------------------------------------
+   ln_rnf_mouth = .true.    !  specific treatment at rivers mouths
+      rn_hrnf     =  15.e0    !  depth over which enhanced vertical mixing is used    (ln_rnf_mouth=T)
+      rn_avt_rnf  =   1.e-3   !  value of the additional vertical mixing coef. [m2/s] (ln_rnf_mouth=T)
+   rn_rfact    =   1.e0    !  multiplicative factor for runoff
+
+   cn_dir = './'  !  root directory for the location of the runoff files
+   !___________!_____________!___________________!___________!_____________!_________!___________!__________!__________!_______________!
+   !           !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+   !           !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
+   sn_rnf      = 'runoff_core_monthly',    -1.   , 'sorunoff',   .true.     , .true. , 'yearly'  , ''       , ''       , ''
+   sn_cnf      = 'runoff_core_monthly',     0.   , 'socoefr0',   .false.    , .true. , 'yearly'  , ''       , ''       , ''
+   sn_s_rnf    = 'runoffs'            ,    24.   , 'rosaline',   .true.     , .true. , 'yearly'  , ''       , ''       , ''
+   sn_t_rnf    = 'runoffs'            ,    24.   , 'rotemper',   .true.     , .true. , 'yearly'  , ''       , ''       , ''
+   sn_dep_rnf  = 'runoffs'            ,     0.   , 'rodepth' ,   .false.    , .true. , 'yearly'  , ''       , ''       , ''
+/
+!-----------------------------------------------------------------------
+&namsbc_wave   ! External fields from wave model                        (ln_wave=T)
+!-----------------------------------------------------------------------
+   ln_sdw      = .false.    !  get the 2D Surf Stokes Drift & Compute the 3D stokes drift
+   ln_stcor    = .false.    !  add Stokes Coriolis and tracer advection terms
+   ln_cdgw     = .false.    !  Neutral drag coefficient read from wave model
+   ln_tauoc    = .false.    !  ocean stress is modified by wave induced stress
+   ln_wave_test= .false.    !  Test case with constant wave fields
+!
+   ln_charn    = .false.     !  Charnock coefficient read from wave model (IFS only)
+   ln_taw      = .false.     !  ocean stress is modified by wave induced stress (coupled mode)
+   ln_phioc    = .false.     !  TKE flux from wave model
+   ln_bern_srfc= .false.     ! wave induced pressure. Bernoulli head J term
+   ln_breivikFV_2016 = .false. ! breivik 2016 vertical stokes profile
+   ln_vortex_force = .false.
+!
+   cn_dir      = './'      !  root directory for the waves data location
+   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
+   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
+   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   !
+   sn_cdg      =  'sdw_ecwaves_orca2'    ,        6.         , 'drag_coeff' ,  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
+   sn_usd      =  'sdw_ecwaves_orca2'    ,        6.         , 'u_sd2d'     ,  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
+   sn_vsd      =  'sdw_ecwaves_orca2'    ,        6.         , 'v_sd2d'     ,  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
+   sn_hsw      =  'sdw_ecwaves_orca2'    ,        6.         , 'hs'         ,  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
+   sn_wmp      =  'sdw_ecwaves_orca2'    ,        6.         , 'wmp'        ,  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
+   sn_wnum     =  'sdw_ecwaves_orca2'    ,        6.         , 'wave_num'   ,  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
+/
+!-----------------------------------------------------------------------
+&namberg       !   iceberg parameters                                   (default: OFF)
+!-----------------------------------------------------------------------
+   ln_icebergs = .true.    ! activate iceberg floats (force =F with "key_agrif")
+
+   cn_dir = './'  !  root directory for the location of drag coefficient files
+   !______!___________!___________________!______________!______________!_________!___________!__________!__________!_______________!
+   !      ! file name ! frequency (hours) !   variable   ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+   !      !           !  (if <0  months)  !     name     !   (logical)  !  (T/F ) ! 'monthly' ! filename ! pairing  ! filename      !
+   sn_icb =  'calving',       -1.         , 'calving'    ,   .true.     , .true.  , 'yearly'  , ''       , ''       , ''
+/
+!!======================================================================
+!!               ***  Lateral boundary condition  ***                 !!
+!!                                                                    !!
+!!   namlbc        lateral momentum boundary condition                  (default: NO selection)
+!!   namagrif      agrif nested grid   (read by child model only)       ("key_agrif")
+!!   nam_tide      Tidal forcing                                        (default: OFF)
+!!   nambdy        Unstructured open boundaries                         (default: OFF)
+!!   nambdy_dta    Unstructured open boundaries - external data         (see  nambdy)
+!!   nambdy_tide   tidal forcing at open boundaries                     (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namlbc        !   lateral momentum boundary condition                  (default: NO selection)
+!-----------------------------------------------------------------------
+   rn_shlat    =    2.     !  no slip
+/
+!-----------------------------------------------------------------------
+&namagrif      !  AGRIF zoom                                            ("key_agrif")
+!-----------------------------------------------------------------------
+   ln_spc_dyn    = .true.  !  use 0 as special value for dynamics
+   rn_sponge_tra = 2880.   !  coefficient for tracer   sponge layer [m2/s]
+   rn_sponge_dyn = 2880.   !  coefficient for dynamics sponge layer [m2/s]
+   ln_chk_bathy  = .false. !  =T  check the parent bathymetry
+/
+!!======================================================================
+!!                ***  Top/Bottom boundary condition  ***             !!
+!!                                                                    !!
+!!   namdrg        top/bottom drag coefficient                          (default: NO selection)
+!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T)
+!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F)
+!!   nambbc        bottom temperature boundary condition                (default: OFF)
+!!   nambbl        bottom boundary layer scheme                         (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namdrg        !   top/bottom drag coefficient                          (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_lin = .true.         !     linear  drag: Cd = Cd0 Uc0
+/
+!-----------------------------------------------------------------------
+&nambbc        !   bottom temperature boundary condition                (default: OFF)
+!-----------------------------------------------------------------------
+   ln_trabbc   = .true.    !  Apply a geothermal heating at the ocean bottom
+      nn_geoflx   =    2      !  geothermal heat flux: = 2 read variable flux [mW/m2]
+
+   cn_dir = './'  !  root directory for the geothermal data location
+   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
+   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
+   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   !
+   sn_qgh      ='geothermal_heating.nc'  ,       -12.        , 'heatflow',   .false.   , .true. , 'yearly'  ,   ''             ,   ''     ,   ''
+/
+!-----------------------------------------------------------------------
+&nambbl        !   bottom boundary layer scheme                         (default: OFF)
+!-----------------------------------------------------------------------
+   ln_trabbl   = .true.    !  Bottom Boundary Layer parameterisation flag
+      nn_bbl_ldf  =  1        !  diffusive bbl (=1)   or not (=0)
+      nn_bbl_adv  =  0        !  advective bbl (=1/2) or not (=0)
+      rn_ahtbbl   =  1000.    !  lateral mixing coefficient in the bbl  [m2/s]
+      rn_gambbl   =  10.      !  advective bbl coefficient                 [s]
+/
+!!======================================================================
+!!                        Tracer (T-S) namelists                      !!
+!!                                                                    !!
+!!   nameos        equation of state                                    (default: NO selection)
+!!   namtra_adv    advection scheme                                     (default: NO selection)
+!!   namtra_ldf    lateral diffusion scheme                             (default: NO selection)
+!!   namtra_mle    mixed layer eddy param. (Fox-Kemper param.)          (default: OFF)
+!!   namtra_eiv    eddy induced velocity param.                         (default: OFF)
+!!   namtra_dmp    T & S newtonian damping                              (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nameos        !   ocean Equation Of Seawater                           (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_eos80    = .true.          !  = Use EOS80
+/
+!-----------------------------------------------------------------------
+&namtra_adv    !   advection scheme for tracer                          (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_traadv_fct = .true.     !  FCT scheme
+      nn_fct_h   =  2               !  =2/4, horizontal 2nd / 4th order 
+      nn_fct_v   =  2               !  =2/4, vertical   2nd / COMPACT 4th order 
+/
+!-----------------------------------------------------------------------
+&namtra_ldf    !   lateral diffusion scheme for tracers                 (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_traldf_lap   = .true.    !    laplacian operator
+   ln_traldf_iso   = .true.    !  iso-neutral (Standard operator)
+   ln_traldf_msc   = .true.    !  Method of Stabilizing Correction      (both operators)
+   !                       !  Coefficients:
+   nn_aht_ijk_t    = 20        !  space/time variation of eddy coef
+      !                             !   = 20     aht = 1/2  Ud. max(e1,e2)
+      rn_Ud        = 0.018          !  lateral diffusive velocity [m/s] (nn_aht_ijk_t= 0, 10, 20, 30)
+      rn_Ld        = 200.e+3        !  lateral diffusive length   [m]   (nn_aht_ijk_t= 0, 10)
+/
+!-----------------------------------------------------------------------
+&namtra_mle    !   mixed layer eddy parametrisation (Fox-Kemper)        (default: OFF)
+!-----------------------------------------------------------------------
+   ln_mle      = .true.   ! (T) use the Mixed Layer Eddy (MLE) parameterisation
+/
+!-----------------------------------------------------------------------
+&namtra_eiv    !   eddy induced velocity param.                         (default: OFF)
+!-----------------------------------------------------------------------
+   ln_ldfeiv   = .true.    ! use eddy induced velocity parameterization
+      !                        !  Coefficients:
+      nn_aei_ijk_t  = 21          ! space/time variation of the eiv coeficient
+      !                                !   = 21 F(i,j,t)  =Treguier et al. JPO 1997 formulation
+      !                           !  time invariant coefficients:  aei0 = 1/2  Ue*Le 
+      rn_Ue        = 0.03              !  lateral diffusive velocity [m/s] (nn_aht_ijk_t= 0, 10, 20, 30)
+      rn_Le        = 200.e+3           !  lateral diffusive length   [m]   (nn_aht_ijk_t= 0, 10)
+      !
+      ln_ldfeiv_dia =.true.   ! diagnose eiv stream function and velocities
+/
+!-----------------------------------------------------------------------
+&namtra_dmp    !   tracer: T & S newtonian damping                      (default: OFF)
+!-----------------------------------------------------------------------
+   ln_tradmp   =  .true.   !  add a damping term (using resto.nc coef.)
+      nn_zdmp  =    0         !  vertical shape =0    damping throughout the water column
+/
+!!======================================================================
+!!                      ***  Dynamics namelists  ***                  !!
+!!                                                                    !!
+!!   nam_vvl       vertical coordinate options                          (default: z-star)
+!!   namdyn_adv    formulation of the momentum advection                (default: NO selection)
+!!   namdyn_vor    advection scheme                                     (default: NO selection)
+!!   namdyn_hpg    hydrostatic pressure gradient                        (default: NO selection)
+!!   namdyn_spg    surface pressure gradient                            (default: NO selection)
+!!   namdyn_ldf    lateral diffusion scheme                             (default: NO selection)
+!!   namdta_dyn    offline TOP: dynamics read in files                  (OFF_SRC only)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namdyn_adv    !   formulation of the momentum advection                (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_dynadv_vec = .true.  !  vector form - 2nd centered scheme
+     nn_dynkeg     = 0        ! grad(KE) scheme: =0   C2  ;  =1   Hollingsworth correction
+/
+!-----------------------------------------------------------------------
+&namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_dynvor_een = .true.  !  energy & enstrophy scheme
+/
+!-----------------------------------------------------------------------
+&namdyn_hpg    !   Hydrostatic pressure gradient option                 (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_hpg_sco  = .true.   !  s-coordinate (standard jacobian formulation)
+/
+!-----------------------------------------------------------------------
+&namdyn_spg    !   surface pressure gradient                            (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_dynspg_ts  = .true.  !  split-explicit free surface
+/
+!-----------------------------------------------------------------------
+&namdyn_ldf    !   lateral diffusion on momentum                        (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_dynldf_lap =  .true.     !    laplacian operator
+   ln_dynldf_lev =  .true.     !  iso-level
+   nn_ahm_ijk_t  = -30         !  =-30  3D coeff. read in eddy_diffusivity_3D.nc
+/
+!!======================================================================
+!!                     vertical physics namelists                     !!
+!!                                                                    !!
+!!    namzdf        vertical physics manager                            (default: NO selection)
+!!    namzdf_ric    richardson number vertical mixing                   (ln_zdfric=T)
+!!    namzdf_tke    TKE vertical mixing                                 (ln_zdftke=T)
+!!    namzdf_gls    GLS vertical mixing                                 (ln_zdfgls=T)
+!!    namzdf_osm    OSM vertical diffusion                              (ln_zdfosm=T)
+!!    namzdf_iwm    tidal mixing parameterization                       (ln_zdfiwm=T)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&namzdf        !   vertical physics manager                             (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_zdftke   = .true.       !  Turbulent Kinetic Energy closure       (T =>   fill namzdf_tke)
+   ln_zdfevd   = .true.       !  Enhanced Vertical Diffusion scheme
+      nn_evdm  =    0            !  evd apply on tracer (=0) or on tracer and momentum (=1)
+      rn_evd   =  100.           !  evd mixing coefficient [m2/s]
+   ln_zdfddm   = .true.    ! double diffusive mixing
+      rn_avts  =    1.e-4     !  maximum avs (vertical mixing on salinity)
+      rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio
+   ln_zdfiwm   = .true.       ! internal wave-induced mixing            (T =>   fill namzdf_iwm)
+   !                       !  Coefficients
+   rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
+   rn_avt0     =   1.2e-5     !  vertical eddy diffusivity [m2/s]       (background Kz if ln_zdfcst=F)
+   nn_avb      =    0         !  profile for background avt & avm (=1) or not (=0)
+   nn_havtb    =    1         !  horizontal shape for avtb (=1) or not (=0)
+/
+!-----------------------------------------------------------------------
+&namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  (ln_zdftke =T)
+!-----------------------------------------------------------------------
+      nn_mxl      =   3       !  mixing length: = 0 bounded by the distance to surface and bottom
+      !                       !                 = 1 bounded by the local vertical scale factor
+      !                       !                 = 2 first vertical derivative of mixing length bounded by 1
+      !                       !                 = 3 as =2 with distinct dissipative an mixing length scale
+      nn_etau     =   1       !  penetration of tke below the mixed layer (ML) due to NIWs
+      !                       !        = 0 none ; = 1 add a tke source below the ML
+      !                       !        = 2 add a tke source just at the base of the ML
+      !                       !        = 3 as = 1 applied on HF part of the stress           (ln_cpl=T)
+      ln_mxhsw    = .false.   !  surface mixing length scale = F(wave height)
+/
+!-----------------------------------------------------------------------
+&namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T)
+!-----------------------------------------------------------------------
+   nn_zpyc     = 2         !  pycnocline-intensified dissipation scales as N (=1) or N^2 (=2)
+   ln_mevar    = .true.    !  variable (T) or constant (F) mixing efficiency
+   ln_tsdiff   = .true.    !  account for differential T/S mixing (T) or not (F)
+
+   cn_dir      = './'      !  root directory for the iwm data location                                                                           
+   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
+   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
+   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   !
+   sn_mpb      = 'int_wave_mix'      , -12.       , 'mixing_power_bot'   , .false.  , .true. , 'yearly' , '' , ''  , ''
+   sn_mpp      = 'int_wave_mix'      , -12.       , 'mixing_power_pyc'   , .false.  , .true. , 'yearly' , '' , ''  , ''
+   sn_mpc      = 'int_wave_mix'      , -12.       , 'mixing_power_cri'   , .false.  , .true. , 'yearly' , '' , ''  , ''
+   sn_dsb      = 'int_wave_mix'      , -12.       , 'decay_scale_bot'    , .false.  , .true. , 'yearly' , '' , ''  , ''
+   sn_dsc      = 'int_wave_mix'      , -12.       , 'decay_scale_cri'    , .false.  , .true. , 'yearly' , '' , ''  , ''
+/
+!!======================================================================
+!!                  ***  Diagnostics namelists  ***                   !!
+!!                                                                    !!
+!!   namtrd       dynamics and/or tracer trends                         (default: OFF)
+!!   namhsb       Heat and salt budgets                                 (default: OFF)
+!!   namdiu       Cool skin and warm layer models                       (default: OFF)
+!!   namdiu       Cool skin and warm layer models                       (default: OFF)
+!!   namflo       float parameters                                      (default: OFF)
+!!   nam_diadct   transports through some sections                      (default: OFF)
+!!   nam_dia25h   25h Mean Output                                       (default: OFF)
+!!   namnc4       netcdf4 chunking and compression settings             ("key_netcdf4")
+!!======================================================================
+!
+
+!-----------------------------------------------------------------------
+&namhsb        !  Heat and salt budgets                                 (default: OFF)
+!-----------------------------------------------------------------------
+   ln_diahsb   = .true.   !  output the heat and salt budgets (T) or not (F)
+/
+!!======================================================================
+!!               ***  Observation & Assimilation  ***                 !!
+!!                                                                    !!
+!!   namobs       observation and model comparison                      (default: OFF)
+!!   nam_asminc   assimilation increments                               ('key_asminc')
+!!======================================================================
+!
+!!======================================================================
+!!                  ***  Miscellaneous namelists  ***                 !!
+!!                                                                    !!
+!!   nammpp            Massively Parallel Processing
+!!   namctl            Control prints                                   (default: OFF)
+!!   namsto            Stochastic parametrization of EOS                (default: OFF)
+!!======================================================================
+!
+!-----------------------------------------------------------------------
+&nammpp        !   Massively Parallel Processing
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namctl        !   Control prints                                       (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsto        ! Stochastic parametrization of EOS                      (default: OFF)
+!-----------------------------------------------------------------------
+/
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_ice_cfg
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_ice_cfg	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_ice_cfg	(revision 14854)
@@ -0,0 +1,84 @@
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! SI3 configuration namelist: Overwrites SHARED/namelist_ice_ref
+!!              1 - Generic parameters                 (nampar)
+!!              2 - Ice thickness discretization       (namitd)
+!!              3 - Ice dynamics                       (namdyn)
+!!              4 - Ice ridging/rafting                (namdyn_rdgrft)
+!!              5 - Ice rheology                       (namdyn_rhg)
+!!              6 - Ice advection                      (namdyn_adv)
+!!              7 - Ice surface boundary conditions    (namsbc)
+!!              8 - Ice thermodynamics                 (namthd)
+!!              9 - Ice heat diffusion                 (namthd_zdf)
+!!             10 - Ice lateral melting                (namthd_da)
+!!             11 - Ice growth in open water           (namthd_do)
+!!             12 - Ice salinity                       (namthd_sal)
+!!             13 - Ice melt ponds                     (namthd_pnd)
+!!             14 - Ice initialization                 (namini)
+!!             15 - Ice/snow albedos                   (namalb)
+!!             16 - Ice diagnostics                    (namdia)
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!
+!------------------------------------------------------------------------------
+&nampar         !   Generic parameters
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namitd         !   Ice discretization
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdyn         !   Ice dynamics
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdyn_rdgrft  !   Ice ridging/rafting
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdyn_rhg     !   Ice rheology
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdyn_adv     !   Ice advection
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namsbc         !   Ice surface boundary conditions
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd         !   Ice thermodynamics
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_zdf     !   Ice heat diffusion
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_da      !   Ice lateral melting
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_do      !   Ice growth in open water
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_sal     !   Ice salinity
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namthd_pnd     !   Melt ponds
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namini         !   Ice initialization
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namalb         !   albedo parameters
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdia         !   Diagnostics
+!------------------------------------------------------------------------------
+/
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_ice_ref
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_ice_ref	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_ice_ref	(revision 14854)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_ice_ref
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_pisces_cfg
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_pisces_cfg	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_pisces_cfg	(revision 14854)
@@ -0,0 +1,141 @@
+!-----------------------------------------------------------------------
+&nampismod     !  Model used 
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nampisext     !   air-sea exchange
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nampisatm     !  Atmospheric prrssure 
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nampisbio     !   biological parameters
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namp4zlim     !   parameters for nutrient limitations for PISCES std  - ln_p4z
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namp5zlim     !   parameters for nutrient limitations PISCES QUOTA    - ln_p5z
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namp5zquota   !   parameters for nutrient limitations PISCES quota    - ln_p5z
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nampisopt     !   parameters for optics
+!-----------------------------------------------------------------------
+/ 
+!-----------------------------------------------------------------------
+&namp4zprod    !   parameters for phytoplankton growth for PISCES std  - ln_p4z
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namp5zprod    !   parameters for phytoplankton growth for PISCES quota- ln_p5z
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namp4zmort    !   parameters for phytoplankton sinks for PISCES std   - ln_p4z
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namp5zmort    !   parameters for phytoplankton sinks for PISCES quota - ln_p5z
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namp4zmes     !   parameters for mesozooplankton for PISCES std       - ln_p4z
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namp5zmes     !   parameters for mesozooplankton
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namp4zzoo     !   parameters for microzooplankton for PISCES std      - ln_p4z
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namp5zzoo     !   parameters for microzooplankton
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nampisfer     !   parameters for iron chemistry
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------  
+&nampisrem     !   parameters for remineralization
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nampispoc     !   parameters for organic particles
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nampiscal     !   parameters for Calcite chemistry
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nampisbc      !   parameters for inputs deposition
+!-----------------------------------------------------------------------
+   sn_dust     = 'dust.orca.new'       ,     -1            , 'dust'        ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
+   ln_ironsed  =  .true.   ! boolean for Fe input from sediments
+   ln_ironice  =  .true.   ! boolean for Fe input from sea ice
+   ln_hydrofe  =  .true.   ! boolean for from hydrothermal vents
+/
+!-----------------------------------------------------------------------
+&nampissed     !   parameters for sediments mobilization
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nampislig     !   Namelist parameters for ligands, nampislig
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nampisice     !   Prescribed sea ice tracers
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&nampisdmp     !   Damping 
+!-----------------------------------------------------------------------
+   nn_pisdmp    =  5840       !  Frequency of Relaxation
+/
+!-----------------------------------------------------------------------
+&nampismass    !   Mass conservation
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namlobphy     !   biological parameters for phytoplankton
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namlobnut     !   biological parameters for nutrients
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namlobzoo     !   biological parameters for zooplankton
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namlobdet     !   biological parameters for detritus
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namlobdom     !   biological parameters for DOM
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namlobsed     !   parameters from aphotic layers to sediment
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namlobrat     !   general coefficients
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namlobopt     !   optical parameters
+!-----------------------------------------------------------------------
+/
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_pisces_ref
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_pisces_ref	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_pisces_ref	(revision 14854)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_pisces_ref
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_ref
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_ref	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_ref	(revision 14854)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_ref
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_top_cfg
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_top_cfg	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_top_cfg	(revision 14854)
@@ -0,0 +1,145 @@
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! NEMO/TOP1 :  Configuration namelist : used to overwrite defaults values defined in SHARED/namelist_top_ref
+!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!-----------------------------------------------------------------------
+&namtrc_run      !   run information
+!-----------------------------------------------------------------------
+   ln_top_euler  = .true.
+/
+!-----------------------------------------------------------------------
+&namtrc          !   tracers definition
+!-----------------------------------------------------------------------
+   jp_bgc        =  24
+!
+   ln_pisces     =  .true.
+   ln_my_trc     =  .false.
+   ln_age        =  .false.
+   ln_cfc11      =  .false.
+   ln_cfc12      =  .false.
+   ln_c14        =  .false.
+!
+   ln_trcdta     =  .true.  !  Initialisation from data input file (T) or not (F)
+   ln_trcbc      =  .true.  !  Enables Boundary conditions
+!                !           !                                           !             !         !
+!                !    name   !           title of the field              !   units     ! init    ! sbc    ! cbc    !  obc    !  ais
+   sn_tracer(1)   = 'DIC     ' , 'Dissolved inorganic Concentration      ',  'mol-C/L' , .true.  , .false., .true. , .false. , .false.
+   sn_tracer(2)   = 'Alkalini' , 'Total Alkalinity Concentration         ',  'eq/L '   , .true.  , .false., .true. , .false. , .false.
+   sn_tracer(3)   = 'O2      ' , 'Dissolved Oxygen Concentration         ',  'mol-C/L' , .true.  , .false., .false., .false. , .false.
+   sn_tracer(4)   = 'CaCO3   ' , 'Calcite Concentration                  ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(5)   = 'PO4     ' , 'Phosphate Concentration                ',  'mol-C/L' , .true.  , .true. , .true. , .false. , .false.
+   sn_tracer(6)   = 'POC     ' , 'Small organic carbon Concentration     ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(7)   = 'Si      ' , 'Silicate Concentration                 ',  'mol-C/L' , .true.  , .true. , .true. , .false. , .false.
+   sn_tracer(8)   = 'PHY     ' , 'Nanophytoplankton Concentration        ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(9)   = 'ZOO     ' , 'Microzooplankton Concentration         ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(10)  = 'DOC     ' , 'Dissolved organic Concentration        ',  'mol-C/L' , .true.  , .false., .true. , .false. , .false.
+   sn_tracer(11)  = 'PHY2    ' , 'Diatoms Concentration                  ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(12)  = 'ZOO2    ' , 'Mesozooplankton Concentration          ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(13)  = 'DSi     ' , 'Diatoms Silicate Concentration         ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(14)  = 'Fer     ' , 'Dissolved Iron Concentration           ',  'mol-C/L' , .true.  , .true. , .true. , .false. , .true.
+   sn_tracer(15)  = 'BFe     ' , 'Big iron particles Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(16)  = 'GOC     ' , 'Big organic carbon Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(17)  = 'SFe     ' , 'Small iron particles Concentration     ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(18)  = 'DFe     ' , 'Diatoms iron  Concentration            ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(19)  = 'GSi     ' , 'Sinking biogenic Silicate Concentration',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(20)  = 'NFe     ' , 'Nano iron Concentration                ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(21)  = 'NCHL    ' , 'Nano chlorophyl Concentration          ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(22)  = 'DCHL    ' , 'Diatoms chlorophyl Concentration       ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+   sn_tracer(23)  = 'NO3     ' , 'Nitrates Concentration                 ',  'mol-C/L' , .true.  , .true. , .true. , .false. , .false.
+   sn_tracer(24)  = 'NH4     ' , 'Ammonium Concentration                 ',  'mol-C/L' , .false. , .false., .false., .false. , .false.
+/
+!-----------------------------------------------------------------------
+&namage          !   AGE 
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtrc_dta      !    Initialisation from data input file
+!-----------------------------------------------------------------------
+!          !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!          !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
+   sn_trcdta(1)  = 'data_DIC_nomask.nc',        -12        ,  'PiDIC'  ,    .false.   , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , ''
+   sn_trcdta(2)  = 'data_ALK_nomask.nc',        -12        ,  'TALK'   ,    .false.   , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , ''
+   sn_trcdta(3)  = 'data_OXY_nomask.nc',        -1         ,  'O2'     ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , ''
+   sn_trcdta(5)  = 'data_PO4_nomask.nc',        -1         ,  'PO4'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , ''
+   sn_trcdta(7)  = 'data_SIL_nomask.nc',        -1         ,  'Si'     ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , ''
+   sn_trcdta(10) = 'data_DOC_nomask.nc',        -1         ,  'DOC'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , ''
+   sn_trcdta(14) = 'data_FER_nomask.nc',        -1         ,  'Fer'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , ''
+   sn_trcdta(23) = 'data_NO3_nomask.nc',        -1         ,  'NO3'    ,    .true.    , .true. , 'yearly'  , 'weights_3D_r360x180_bilin.nc'       , ''   , ''
+   rn_trfac(1)   =   1.028e-06  !  multiplicative factor
+   rn_trfac(2)   =   1.028e-06  !  -      -      -     -
+   rn_trfac(3)   =  44.6e-06  !  -      -      -     -
+   rn_trfac(5)   = 122.0e-06  !  -      -      -     -
+   rn_trfac(7)   =   1.0e-06  !  -      -      -     -
+   rn_trfac(10)  =   1.0e-06  !  -      -      -     -
+   rn_trfac(14)  =   1.0e-06  !  -      -      -     -
+   rn_trfac(23)  =   7.6e-06  !  -      -      -     -
+/
+!-----------------------------------------------------------------------
+&namtrc_adv      !   advection scheme for passive tracer                (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_trcadv_mus =  .true.   !  MUSCL scheme
+      ln_mus_ups =  .false.         !  use upstream scheme near river mouths
+/
+!-----------------------------------------------------------------------
+&namtrc_ldf      !   lateral diffusion scheme for passive tracer        (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_trcldf_tra   =  .true.     !  use active tracer setting
+/
+!-----------------------------------------------------------------------
+&namtrc_rad      !  treatment of negative concentrations 
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtrc_snk      !  sedimentation of particles
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtrc_dmp      !   passive tracer newtonian damping   
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtrc_ice      !    Representation of sea ice growth & melt effects
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtrc_trd      !   diagnostics on tracer trends                       ('key_trdtrc')
+!----------------------------------------------------------------------
+/
+!----------------------------------------------------------------------
+&namtrc_bc       !   data for boundary conditions
+!-----------------------------------------------------------------------
+!                !  file name        ! frequency (hours) ! variable      ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!                !                   !  (if <0  months)  !   name        !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
+   sn_trcsbc(5)  = 'dust.orca.new'   ,       -1          , 'dustpo4'     ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
+   sn_trcsbc(7)  = 'dust.orca.new'   ,       -1          , 'dustsi'      ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
+   sn_trcsbc(14) = 'dust.orca.new'   ,       -1          , 'dustfer'     ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
+   sn_trcsbc(23) = 'ndeposition.orca',      -12          , 'ndep'        ,  .false.     , .true. , 'yearly'  , ''       , ''    , ''
+   rn_trsfac(5)  = 8.264e-02   !  (  0.021 / 31. * 122 )
+   rn_trsfac(7)  = 3.313e-01     !  ( 8.8   / 28.1 )
+   rn_trsfac(14) = 6.266e-04   !  (  0.035 / 55.85 )
+   rn_trsfac(23) =  5.4464e-01  !  ( From kgN m-2 s-1 to molC l-1 ====> zfact = 7.625/14 )
+   rn_sbc_time   =  1.          !  Time scaling factor for SBC and CBC data (seconds in a day)
+   !
+   sn_trccbc(1)  = 'river.orca'      ,    -12            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
+   sn_trccbc(2)  = 'river.orca'      ,    -12            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
+   sn_trccbc(5)  = 'river.orca'      ,    -12            , 'riverdip'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
+   sn_trccbc(7)  = 'river.orca'      ,    -12            , 'riverdsi'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
+   sn_trccbc(10) = 'river.orca'      ,    -12            , 'riverdoc'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
+   sn_trccbc(14) = 'river.orca'      ,    -12            , 'riverdic'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
+   sn_trccbc(23) = 'river.orca'      ,    -12            , 'riverdin'    ,  .true.      , .true. , 'yearly'  , ''       , ''    , ''
+   rn_trcfac(1)  = 8.333e+01   !  ( data in Mg/m2/yr : 1e3/12/ryyss)
+   rn_trcfac(2)  = 8.333e+01   !  ( 1e3 /12 )
+   rn_trcfac(5)  = 3.935e+04   !  ( 1e3 / 31. * 122 )
+   rn_trcfac(7)  = 3.588e+01   !  ( 1e3 / 28.1 )
+   rn_trcfac(10) = 8.333e+01   !  ( 1e3 / 12
+   rn_trcfac(14) = 4.166e-03   !  ( 1e3 / 12 * 5e-5 )
+   rn_trcfac(23) = 5.446e+02   !  (  1e3 / 14 * 7.625 )
+   rn_cbc_time   = 3.1536e+7   !  Time scaling factor for CBC data (seconds in a year)
+/
+!----------------------------------------------------------------------
+&namtrc_bdy      !   Setup of tracer boundary conditions
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namtrc_ais      !  Representation of Antarctic Ice Sheet tracers supply
+!-----------------------------------------------------------------------
+/
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_top_ref
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_top_ref	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_top_ref	(revision 14854)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_top_ref
Index: NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_trc_ref
===================================================================
--- NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_trc_ref	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/EXPREF/namelist_trc_ref	(revision 14854)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_trc_ref
Index: NEMO/trunk/tests/DIA_GPU/MY_SRC/cuda_fortran.F90
===================================================================
--- NEMO/trunk/tests/DIA_GPU/MY_SRC/cuda_fortran.F90	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/MY_SRC/cuda_fortran.F90	(revision 14854)
@@ -0,0 +1,118 @@
+#ifdef key_gpu
+MODULE cuda_fortran
+    USE cudafor
+    USE dom_oce
+    USE lib_mpp
+    CONTAINS
+      ATTRIBUTES(global) &
+      SUBROUTINE dia_hsb_kernel(surf, e3t, surf_ini, e3t_ini, ts, hc_loc_ini, sc_loc_ini,     &
+                            & tmask, tmask_ini, zwrkv, zwrkh, zwrks, zwrk, jpi, jpj, jpk, jpt, Kmm)
+            !
+            INTEGER, VALUE :: jpi, jpj, jpk, jpt                                                                         ! MPI process sub-domain and stream tile
+            INTEGER, VALUE :: Kmm                                                                                   !ocean time level indice
+            REAL(8)        :: surf(jpi, jpj), e3t(jpi, jpj, jpk, jpt), surf_ini(jpi, jpj), e3t_ini(jpi, jpj, jpk), ts(jpi, jpj, jpk, 2, jpt)
+            REAL(8)        :: hc_loc_ini(jpi, jpj, jpk), sc_loc_ini(jpi, jpj, jpk), tmask(jpi, jpj, jpk), tmask_ini(jpi, jpj, jpk)
+            REAL(8)        :: zwrkv(jpi, jpj, jpk), zwrkh(jpi, jpj, jpk), zwrks(jpi, jpj, jpk), zwrk(jpi, jpj, jpk)
+            !
+
+            INTEGER        :: i, j, k                       ! dummy indexes
+
+            !
+            i = blockDim%x * (blockIdx%x -1) + threadIdx%x
+            j = blockDim%y * (blockIdx%y -1) + threadIdx%y
+            k = blockDim%z * (blockIdx%z -1) + threadIdx%z
+            !
+            IF ( (i .le. jpi ) .AND. (j .le. jpj) .AND. (k .le. jpk-1) ) THEN
+                zwrkv(i, j, k) =   surf    (i,j) * e3t    (i,j,k,Kmm) * tmask    (i,j,k)      &
+                &                - surf_ini(i,j) * e3t_ini(i,j,k)     * tmask_ini(i,j,k)
+                
+                zwrkh(i, j, k) = ( surf(i, j) * e3t(i, j, k, Kmm) * ts(i, j, k, 1, Kmm) - surf_ini(i, j) * hc_loc_ini(i, j, k) )
+                
+                zwrks(i, j, k) = ( surf(i, j) * e3t(i, j, k, Kmm) * ts(i, j, k, 2, Kmm) - surf_ini(i, j) * sc_loc_ini(i, j, k) )
+                
+                zwrk(i, j, k)  =   surf(i, j) * e3t(i, j, k, Kmm) * tmask(i, j, k)
+            END IF
+      END SUBROUTINE dia_hsb_kernel
+
+!     ATTRIBUTES(global) &
+!     SUBROUTINE dia_hsb_kernel1d(surf, e3t_n, surf_ini, e3t_ini, &
+!        &  tsn, hc_loc_ini, sc_loc_ini, tmask, zwrkv, zwrkh, zwrks, zwrk, jpi, jpj, jpk)
+!        IMPLICIT NONE
+!        REAL(kind=8) :: surf(:, :), e3t_n(:, :, :), surf_ini(:, :), e3t_ini(:, :, :), tsn(:, :, :, :), &
+!                            & hc_loc_ini(:, :, :), sc_loc_ini(:, :, :), tmask(:, :, :), zwrkv(:, :, :),      &
+!                            & zwrkh(:, :, :), zwrks(:, :, :), zwrk(:, :, :)
+!        REAL(kind=8), SHARED       :: sdata(*)
+!        INTEGER, VALUE             :: jpi, jpj, jpk
+!        INTEGER                    :: i, si, ti, globsize
+!
+!        globsize = jpi*jpj*jpk
+!
+!        i  = blockDim%x * (blockIdx%x-1) + threadIdx%x
+!        ti = threadIdx%x
+!        si = MOD(i, jpk)+1 !Shared memory indexing
+!
+!        !Compute volume (zwrk), volume variation (zwrkv), heat deviation (zwrkh)
+!        !and salinity deviation (zwrks)
+!        IF ( i .le. globsize ) THEN
+!            sdata(ti) = surf(si) * e3t_n(i)
+!            zwrkv(i) = ( sdata(ti) - surf_ini(si) * e3t_ini(i) ) * tmask(i) * surf(si)
+!            zwrkh(i) = ( sdata(ti) * tsn(i) - surf_ini(si) * hc_loc_ini(i) ) * tmask(i) * surf(si)
+!            zwrks(i) = ( sdata(ti) * tsn(globsize + i) - surf_ini(si) * sc_loc_ini(i) ) &
+!                    * tmask(i) * surf(si)
+!            zwrk(i) = sdata(ti) * tmask(i) * surf(si)
+!        END IF
+!     END SUBROUTINE dia_hsb_kernel1d
+
+
+       ATTRIBUTES(global) &
+       SUBROUTINE filter_cuda(ptab, mask, jpi, jpj, jpk)
+           !
+           REAL(kind=8)         :: ptab(jpi,jpj,jpk) , mask(jpi,jpj)
+           !
+           INTEGER, VALUE  :: jpi, jpj, jpk                 ! MPI process sub-domain
+           INTEGER         :: i, j, k                 ! dummy indexes
+           !
+           i = blockDim%x * (blockIdx%x -1) + threadIdx%x
+           j = blockDim%y * (blockIdx%y -1) + threadIdx%y
+           k = blockDim%z * (blockIdx%z -1) + threadIdx%z
+           !
+           tile = 1
+           IF ( (i .le. jpi ) .AND. (j .le. jpj) .AND. (k .le. jpk-1) ) THEN
+               ptab(i, j, k)  = ptab(i, j, k) * mask(i, j)
+           END IF
+       END SUBROUTINE filter_cuda
+
+       ATTRIBUTES(global) SUBROUTINE array3dto1d(d_inp,d_out,ipi, ipj, ipk)
+           !
+           REAL(kind=8)                 :: d_inp(:,:,:)
+           REAL(kind=8), intent(out)    :: d_out(:)
+           !
+           INTEGER, VALUE  :: ipi, ipj, ipk                 ! MPI process sub-domain
+           INTEGER         :: i, j, k                       ! dummy indexes
+           !
+           i = blockDim%x * (blockIdx%x -1) + threadIdx%x
+           j = blockDim%y * (blockIdx%y -1) + threadIdx%y
+           k = blockDim%z * (blockIdx%z -1) + threadIdx%z
+           !
+           IF ( (i .le. ipi ) .AND. (j .le. ipj) .AND. (k .le. ipk) ) THEN
+               d_out(i + ipi*(j-1) + ipi*ipj*(k-1) )  = d_inp(i, j, k)
+           END IF
+       END SUBROUTINE array3dto1d
+
+       !Knuth's trick
+       ATTRIBUTES(device) &
+       SUBROUTINE DDPDD_d(ydda, yddb)
+           IMPLICIT NONE
+           COMPLEX(kind=8), INTENT(in   ) :: ydda              !Scalar to add
+           COMPLEX(kind=8), INTENT(inout) :: yddb              !Total sum
+           REAL(kind=8)                   :: zerr, zt1, zt2
+       
+           zt1  = REAL(ydda) + REAL(yddb)
+           zerr = zt1 - REAL(ydda)
+           zt2  = ( (REAL(yddb) - zerr) + (REAL(ydda) - (zt1 - zerr)) ) &
+                 & + AIMAG(ydda) + AIMAG(yddb)
+           !The result is t1 + t2 after normalization
+           yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), 8 )
+       END SUBROUTINE DDPDD_d       
+END MODULE
+#endif
Index: NEMO/trunk/tests/DIA_GPU/MY_SRC/diahsb.F90
===================================================================
--- NEMO/trunk/tests/DIA_GPU/MY_SRC/diahsb.F90	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/MY_SRC/diahsb.F90	(revision 14854)
@@ -0,0 +1,419 @@
+#ifdef key_gpu
+#include "diahsb_gpu.h90"
+#else
+MODULE diahsb
+   !!======================================================================
+   !!                       ***  MODULE  diahsb  ***
+   !! Ocean diagnostics: Heat, salt and volume budgets
+   !!======================================================================
+   !! History :  3.3  ! 2010-09  (M. Leclair)  Original code
+   !!                 ! 2012-10  (C. Rousset)  add iom_put
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dia_hsb       : Diagnose the conservation of ocean heat and salt contents, and volume
+   !!   dia_hsb_rst   : Read or write DIA file in restart file
+   !!   dia_hsb_init  : Initialization of the conservation diagnostic
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE sbc_oce        ! surface thermohaline fluxes
+   USE isf_oce        ! ice shelf fluxes
+   USE sbcrnf         ! river runoff
+   USE domvvl         ! vertical scale factors
+   USE traqsr         ! penetrative solar radiation
+   USE trabbc         ! bottom boundary condition
+   USE trabbc         ! bottom boundary condition
+   USE restart        ! ocean restart
+   USE bdy_oce , ONLY : ln_bdy
+   !
+   USE iom            ! I/O manager
+   USE in_out_manager ! I/O manager
+   USE lib_fortran    ! glob_sum
+   USE lib_mpp        ! distributed memory computing library
+   USE timing         ! preformance summary
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   dia_hsb        ! routine called by step.F90
+   PUBLIC   dia_hsb_init   ! routine called by nemogcm.F90
+
+   LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets
+
+   REAL(wp) ::   surf_tot              ! ocean surface
+   REAL(wp) ::   frc_t, frc_s, frc_v   ! global forcing trends
+   REAL(wp) ::   frc_wn_t, frc_wn_s    ! global forcing trends
+   !
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf_ini      , ssh_ini          !
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   !
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   tmask_ini
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_hsb( kt, Kbb, Kmm )
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_hsb  ***
+      !!
+      !! ** Purpose: Compute the ocean global heat content, salt content and volume conservation
+      !!
+      !! ** Method : - Compute the deviation of heat content, salt content and volume
+      !!	            at the current time step from their values at nit000
+      !!	            - Compute the contribution of forcing and remove it from these deviations
+      !!
+      !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt         ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
+      !
+      INTEGER    ::   ji, jj, jk                  ! dummy loop indice
+      REAL(wp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
+      REAL(wp)   ::   zdiff_hc1   , zdiff_sc1     !  -         -     -        -
+      REAL(wp)   ::   zdiff_v1    , zdiff_v2      ! volume variation
+      REAL(wp)   ::   zerr_hc1    , zerr_sc1      ! heat and salt content misfit
+      REAL(wp)   ::   zvol_tot                    ! volume
+      REAL(wp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
+      REAL(wp)   ::   z_frc_trd_v                 !    -     -
+      REAL(wp)   ::   z_wn_trd_t , z_wn_trd_s     !    -     -
+      REAL(wp)   ::   z_ssh_hc , z_ssh_sc         !    -     -
+      REAL(wp), DIMENSION(jpi,jpj)       ::   z2d0, z2d1   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwrk         ! 3D workspace
+      !!---------------------------------------------------------------------------
+      IF( ln_timing )   CALL timing_start('dia_hsb')
+      !
+      ts(:,:,:,1,Kmm) = ts(:,:,:,1,Kmm) * tmask(:,:,:) ; ts(:,:,:,1,Kbb) = ts(:,:,:,1,Kbb) * tmask(:,:,:) ;
+      ts(:,:,:,2,Kmm) = ts(:,:,:,2,Kmm) * tmask(:,:,:) ; ts(:,:,:,2,Kbb) = ts(:,:,:,2,Kbb) * tmask(:,:,:) ;
+      ! ------------------------- !
+      ! 1 - Trends due to forcing !
+      ! ------------------------- !
+      z_frc_trd_v = r1_rho0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * surf(:,:) )   ! volume fluxes
+      z_frc_trd_t =           glob_sum( 'diahsb', sbc_tsc(:,:,jp_tem) * surf(:,:) )                       ! heat fluxes
+      z_frc_trd_s =           glob_sum( 'diahsb', sbc_tsc(:,:,jp_sal) * surf(:,:) )                       ! salt fluxes
+      !                    !  Add runoff    heat & salt input
+      IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( 'diahsb', rnf_tsc(:,:,jp_tem) * surf(:,:) )
+      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( 'diahsb', rnf_tsc(:,:,jp_sal) * surf(:,:) )
+      !                    ! Add ice shelf heat & salt input
+      IF( ln_isf    )   z_frc_trd_t = z_frc_trd_t &
+         &                          + glob_sum( 'diahsb', ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) )
+      !                    ! Add penetrative solar radiation
+      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rho0_rcp * glob_sum( 'diahsb', qsr     (:,:) * surf(:,:) )
+      !                    ! Add geothermal heat flux
+      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( 'diahsb', qgh_trd0(:,:) * surf(:,:) )
+      !
+      IF( ln_linssh ) THEN
+         IF( ln_isfcav ) THEN
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  z2d0(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb)
+                  z2d1(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb)
+               END DO
+            END DO
+         ELSE
+            z2d0(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_tem,Kbb)
+            z2d1(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_sal,Kbb)
+         END IF
+         z_wn_trd_t = - glob_sum( 'diahsb', z2d0 )
+         z_wn_trd_s = - glob_sum( 'diahsb', z2d1 )
+      ENDIF
+
+      frc_v = frc_v + z_frc_trd_v * rn_Dt
+      frc_t = frc_t + z_frc_trd_t * rn_Dt
+      frc_s = frc_s + z_frc_trd_s * rn_Dt
+      !                                          ! Advection flux through fixed surface (z=0)
+      IF( ln_linssh ) THEN
+         frc_wn_t = frc_wn_t + z_wn_trd_t * rn_Dt
+         frc_wn_s = frc_wn_s + z_wn_trd_s * rn_Dt
+      ENDIF
+
+      ! ------------------------ !
+      ! 2 -  Content variations  !
+      ! ------------------------ !
+      ! glob_sum_full is needed because you keep the full interior domain to compute the sum (iscpl)
+
+      !                    ! volume variation (calculated with ssh)
+      zdiff_v1 = glob_sum_full( 'diahsb', surf(:,:)*ssh(:,:,Kmm) - surf_ini(:,:)*ssh_ini(:,:) )
+
+      !                    ! heat & salt content variation (associated with ssh)
+      IF( ln_linssh ) THEN       ! linear free surface case
+         IF( ln_isfcav ) THEN          ! ISF case
+            DO ji = 1, jpi
+               DO jj = 1, jpj
+                  z2d0(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) )
+                  z2d1(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) )
+               END DO
+            END DO
+         ELSE                          ! no under ice-shelf seas
+            z2d0(:,:) = surf(:,:) * ( ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) - ssh_hc_loc_ini(:,:) )
+            z2d1(:,:) = surf(:,:) * ( ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) - ssh_sc_loc_ini(:,:) )
+         END IF
+         z_ssh_hc = glob_sum_full( 'diahsb', z2d0 )
+         z_ssh_sc = glob_sum_full( 'diahsb', z2d1 )
+      ENDIF
+      !
+      DO jk = 1, jpkm1           ! volume variation (calculated with scale factors)
+         zwrk(:,:,jk) =   surf    (:,:) * e3t    (:,:,jk,Kmm)*tmask    (:,:,jk)   &
+            &           - surf_ini(:,:) * e3t_ini(:,:,jk    )*tmask_ini(:,:,jk)
+      END DO
+      zdiff_v2 = glob_sum_full( 'diahsb', zwrk(:,:,:) )     ! glob_sum_full needed as tmask and tmask_ini could be different
+      DO jk = 1, jpkm1           ! heat content variation
+         zwrk(:,:,jk) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm)   &
+            &           - surf_ini(:,:) *         hc_loc_ini(:,:,jk) )
+      END DO
+      zdiff_hc = glob_sum_full( 'diahsb', zwrk(:,:,:) )
+      DO jk = 1, jpkm1           ! salt content variation
+         zwrk(:,:,jk) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm)   &
+            &           - surf_ini(:,:) *         sc_loc_ini(:,:,jk) )
+      END DO
+      zdiff_sc = glob_sum_full( 'diahsb', zwrk(:,:,:) )
+
+      ! ------------------------ !
+      ! 3 -  Drifts              !
+      ! ------------------------ !
+      zdiff_v1 = zdiff_v1 - frc_v
+      IF( .NOT.ln_linssh )   zdiff_v2 = zdiff_v2 - frc_v
+      zdiff_hc = zdiff_hc - frc_t
+      zdiff_sc = zdiff_sc - frc_s
+      IF( ln_linssh ) THEN
+         zdiff_hc1 = zdiff_hc + z_ssh_hc
+         zdiff_sc1 = zdiff_sc + z_ssh_sc
+         zerr_hc1  = z_ssh_hc - frc_wn_t
+         zerr_sc1  = z_ssh_sc - frc_wn_s
+      ENDIF
+
+      ! ----------------------- !
+      ! 4 - Diagnostics writing !
+      ! ----------------------- !
+      DO jk = 1, jpkm1           ! total ocean volume (calculated with scale factors)
+         zwrk(:,:,jk) = surf(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+      END DO
+      zvol_tot = glob_sum( 'diahsb', zwrk(:,:,:) )
+
+!!gm to be added ?
+!      IF( ln_linssh ) THEN            ! fixed volume, add the ssh contribution
+!        zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * ssh(:,:,Kmm) )
+!      ENDIF
+!!gm end
+
+      CALL iom_put(   'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3)
+      CALL iom_put(   'bgfrctem' , frc_t    * rho0 * rcp * 1.e-20 )   ! hc  - surface forcing (1.e20 J)
+      CALL iom_put(   'bgfrchfx' , frc_t    * rho0 * rcp /  &         ! hc  - surface forcing (W/m2)
+         &                       ( surf_tot * kt * rn_Dt )        )
+      CALL iom_put(   'bgfrcsal' , frc_s    * 1.e-9    )              ! sc  - surface forcing (psu*km3)
+
+      IF( .NOT. ln_linssh ) THEN
+         CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature drift     (C)
+         CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot )              ! Salinity    drift     (PSU)
+         CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rho0 * rcp )   ! Heat content drift    (1.e20 J)
+         CALL iom_put( 'bgheatfx' , zdiff_hc * rho0 * rcp /  &         ! Heat flux drift       (W/m2)
+            &                       ( surf_tot * kt * rn_Dt )        )
+         CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9    )              ! Salt content drift    (psu*km3)
+         CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh drift      (km3)
+         CALL iom_put( 'bgvole3t' , zdiff_v2 * 1.e-9    )              ! volume e3t drift      (km3)
+         !
+         IF( kt == nitend .AND. lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'dia_hsb : last time step hsb diagnostics: at it= ', kt,' date= ', ndastp
+            WRITE(numout,*) '~~~~~~~'
+            WRITE(numout,*) '   Temperature drift = ', zdiff_hc / zvol_tot, ' C'
+            WRITE(numout,*) '   Salinity    drift = ', zdiff_sc / zvol_tot, ' PSU'
+            WRITE(numout,*) '   volume ssh  drift = ', zdiff_v1 * 1.e-9   , ' km^3'
+            WRITE(numout,*) '   volume e3t  drift = ', zdiff_v2 * 1.e-9   , ' km^3'
+         ENDIF
+         !
+      ELSE
+         CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content drift    (C)
+         CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content drift    (PSU)
+         CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rho0 * rcp )  ! Heat content drift    (1.e20 J)
+         CALL iom_put( 'bgheatfx' , zdiff_hc1 * rho0 * rcp /  &        ! Heat flux drift       (W/m2)
+            &                       ( surf_tot * kt * rn_Dt )         )
+         CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content drift    (psu*km3)
+         CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh drift      (km3)
+         CALL iom_put( 'bgmistem' , zerr_hc1 / zvol_tot )              ! hc  - error due to free surface (C)
+         CALL iom_put( 'bgmissal' , zerr_sc1 / zvol_tot )              ! sc  - error due to free surface (psu)
+      ENDIF
+      !
+      IF( lrst_oce )   CALL dia_hsb_rst( kt, Kmm, 'WRITE' )
+      !
+      IF( ln_timing )   CALL timing_stop('dia_hsb')
+      !
+   END SUBROUTINE dia_hsb
+
+
+   SUBROUTINE dia_hsb_rst( kt, Kmm, cdrw )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE dia_hsb_rst  ***
+      !!
+      !! ** Purpose :   Read or write DIA file in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      INTEGER         , INTENT(in) ::   Kmm    ! ocean time level index
+      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      !
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            !
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   dia_hsb_rst : read hsb restart at it= ', kt,' date= ', ndastp
+            IF(lwp) WRITE(numout,*)
+            CALL iom_get( numror, 'frc_v', frc_v )
+            CALL iom_get( numror, 'frc_t', frc_t )
+            CALL iom_get( numror, 'frc_s', frc_s )
+            IF( ln_linssh ) THEN
+               CALL iom_get( numror, 'frc_wn_t', frc_wn_t )
+               CALL iom_get( numror, 'frc_wn_s', frc_wn_s )
+            ENDIF
+            CALL iom_get( numror, jpdom_auto, 'surf_ini'  , surf_ini   ) ! ice sheet coupling
+            CALL iom_get( numror, jpdom_auto, 'ssh_ini'   , ssh_ini    )
+            CALL iom_get( numror, jpdom_auto, 'e3t_ini'   , e3t_ini    )
+            CALL iom_get( numror, jpdom_auto, 'tmask_ini' , tmask_ini  )
+            CALL iom_get( numror, jpdom_auto, 'hc_loc_ini', hc_loc_ini )
+            CALL iom_get( numror, jpdom_auto, 'sc_loc_ini', sc_loc_ini )
+            IF( ln_linssh ) THEN
+               CALL iom_get( numror, jpdom_auto, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
+               CALL iom_get( numror, jpdom_auto, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
+            ENDIF
+         ELSE
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   dia_hsb_rst : initialise hsb at initial state '
+            IF(lwp) WRITE(numout,*)
+            surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:)         ! initial ocean surface
+            ssh_ini(:,:) = ssh(:,:,Kmm)                          ! initial ssh
+            DO jk = 1, jpk
+              ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
+               e3t_ini   (:,:,jk) = e3t(:,:,jk,Kmm)                      * tmask(:,:,jk)  ! initial vertical scale factors
+               tmask_ini (:,:,jk) = tmask(:,:,jk)                                       ! initial mask
+               hc_loc_ini(:,:,jk) = ts(:,:,jk,jp_tem,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial heat content
+               sc_loc_ini(:,:,jk) = ts(:,:,jk,jp_sal,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial salt content
+            END DO
+            frc_v = 0._wp                                           ! volume       trend due to forcing
+            frc_t = 0._wp                                           ! heat content   -    -   -    -
+            frc_s = 0._wp                                           ! salt content   -    -   -    -
+            IF( ln_linssh ) THEN
+               IF( ln_isfcav ) THEN
+                  DO ji = 1, jpi
+                     DO jj = 1, jpj
+                        ssh_hc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm)   ! initial heat content in ssh
+                        ssh_sc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm)   ! initial salt content in ssh
+                     END DO
+                   END DO
+                ELSE
+                  ssh_hc_loc_ini(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm)   ! initial heat content in ssh
+                  ssh_sc_loc_ini(:,:) = ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm)   ! initial salt content in ssh
+               END IF
+               frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface
+               frc_wn_s = 0._wp                                       ! initial salt content misfit due to free surface
+            ENDIF
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   dia_hsb_rst : write restart at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*)
+         !
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_v', frc_v )
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_t', frc_t )
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_s', frc_s )
+         IF( ln_linssh ) THEN
+            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_t', frc_wn_t )
+            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s )
+         ENDIF
+         CALL iom_rstput( kt, nitrst, numrow, 'surf_ini'  , surf_ini   ) ! ice sheet coupling
+         CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini'   , ssh_ini    )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini'   , e3t_ini    )
+         CALL iom_rstput( kt, nitrst, numrow, 'tmask_ini' , tmask_ini  )
+         CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini )
+         CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini )
+         IF( ln_linssh ) THEN
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE dia_hsb_rst
+
+
+   SUBROUTINE dia_hsb_init( Kmm )
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_hsb  ***
+      !!
+      !! ** Purpose: Initialization for the heat salt volume budgets
+      !!
+      !! ** Method : Compute initial heat content, salt content and volume
+      !!
+      !! ** Action : - Compute initial heat content, salt content and volume
+      !!             - Initialize forcing trends
+      !!             - Compute coefficients for conversion
+      !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kmm ! time level index
+      !
+      INTEGER ::   ierror, ios   ! local integer
+      !!
+      NAMELIST/namhsb/ ln_diahsb
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_hsb_init : heat and salt budgets diagnostics'
+         WRITE(numout,*) '~~~~~~~~~~~~ '
+      ENDIF
+      READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namhsb in reference namelist' )
+      READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namhsb in configuration namelist' )
+      IF(lwm) WRITE( numond, namhsb )
+
+      IF(lwp) THEN
+         WRITE(numout,*) '   Namelist  namhsb :'
+         WRITE(numout,*) '      check the heat and salt budgets (T) or not (F)       ln_diahsb = ', ln_diahsb
+      ENDIF
+      !
+      IF( .NOT. ln_diahsb )   RETURN
+
+      ! ------------------- !
+      ! 1 - Allocate memory !
+      ! ------------------- !
+      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), &
+         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), tmask_ini(jpi,jpj,jpk),STAT=ierror  )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb_init: unable to allocate hc_loc_ini' )   ;   RETURN
+      ENDIF
+
+      IF( ln_linssh )   ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to allocate ssh_hc_loc_ini' )   ;   RETURN
+      ENDIF
+
+      ! ----------------------------------------------- !
+      ! 2 - Time independant variables and file opening !
+      ! ----------------------------------------------- !
+      surf(:,:) = e1e2t(:,:) * tmask_i(:,:)               ! masked surface grid cell area
+      surf_tot  = glob_sum( 'diahsb', surf(:,:) )         ! total ocean surface area
+
+      IF( ln_bdy ) CALL ctl_warn( 'dia_hsb_init: heat/salt budget does not consider open boundary fluxes' )
+      !
+      ! ---------------------------------- !
+      ! 4 - initial conservation variables !
+      ! ---------------------------------- !
+      CALL dia_hsb_rst( nit000, Kmm, 'READ' )  !* read or initialize all required files
+      !
+   END SUBROUTINE dia_hsb_init
+
+   !!======================================================================
+END MODULE diahsb
+#endif
Index: NEMO/trunk/tests/DIA_GPU/MY_SRC/diahsb_gpu.h90
===================================================================
--- NEMO/trunk/tests/DIA_GPU/MY_SRC/diahsb_gpu.h90	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/MY_SRC/diahsb_gpu.h90	(revision 14854)
@@ -0,0 +1,666 @@
+MODULE diahsb
+   !!======================================================================
+   !!                       ***  MODULE  diahsb  ***
+   !! Ocean diagnostics: Heat, salt and volume budgets
+   !!======================================================================
+   !! History :  3.3  ! 2010-09  (M. Leclair)  Original code
+   !!                 ! 2012-10  (C. Rousset)  add iom_put
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dia_hsb       : Diagnose the conservation of ocean heat and salt contents, and volume
+   !!   dia_hsb_rst   : Read or write DIA file in restart file
+   !!   dia_hsb_init  : Initialization of the conservation diagnostic
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE sbc_oce        ! surface thermohaline fluxes
+   USE isf_oce        ! ice shelf fluxes
+   USE sbcrnf         ! river runoff
+   USE domvvl         ! vertical scale factors
+   USE traqsr         ! penetrative solar radiation
+   USE trabbc         ! bottom boundary condition
+   USE trabbc         ! bottom boundary condition
+   USE restart        ! ocean restart
+   USE bdy_oce , ONLY : ln_bdy
+   !
+   USE iom            ! I/O manager
+   USE in_out_manager ! I/O manager
+   USE gpu_manager    ! GPU manager
+   USE cudafor        ! CUDA toolkit libs
+   USE cuda_fortran   ! CUDA routines
+   !USE nvtx          ! CUDA profiling/DEGUG tools
+   USE lib_fortran    ! glob_sum
+   USE lib_mpp        ! distributed memory computing library
+   USE timing         ! preformance summary
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   dia_hsb        ! routine called by step.F90
+   PUBLIC   dia_hsb_init   ! routine called by nemogcm.F90
+
+   LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets
+
+   REAL(wp)                      :: surf_tot ! ocean surface
+   REAL(wp) , DIMENSION(2), SAVE :: frc_t, frc_s, frc_v ! global forcing trends
+   REAL(wp) , DIMENSION(2), SAVE :: frc_wn_t, frc_wn_s ! global forcing trends
+   !
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf_ini      , ssh_ini          !
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   !
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, PINNED :: hc_loc_ini, sc_loc_ini !
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: e3t_ini !
+   REAL(wp), DIMENSION(:) , ALLOCATABLE, PINNED, SAVE :: h_ztmpv, h_ztmph, h_ztmps, h_ztmp !
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmask_ini
+
+
+   !Device data associate to PUBLIC arrays
+   REAL(8), DIMENSION(:,:,:,:) , ALLOCATABLE, DEVICE :: d_e3t !
+   REAL(8), DIMENSION(:,:,:) , ALLOCATABLE, DEVICE :: d_tmask !
+   REAL(8), DIMENSION(:,:) , ALLOCATABLE, DEVICE :: d_tmask_h !
+   REAL(8), DIMENSION(:,:,:) , ALLOCATABLE, DEVICE :: d_tmask_ini !
+   REAL(8), DIMENSION(:,:,:,:,:), ALLOCATABLE, DEVICE :: d_ts !
+   !Device data associate to LOCAL/DEVICE arrays
+   REAL(8), DEVICE , DIMENSION(:,:) , ALLOCATABLE :: d_surf !
+   REAL(8), DEVICE , DIMENSION(:,:) , ALLOCATABLE :: d_surf_ini !
+   REAL(8), DEVICE , DIMENSION(:,:,:) , ALLOCATABLE :: d_hc_loc_ini !
+   REAL(8), DEVICE , DIMENSION(:,:,:) , ALLOCATABLE :: d_sc_loc_ini !
+   REAL(8), DEVICE , DIMENSION(:,:,:) , ALLOCATABLE :: d_e3t_ini !
+   REAL(8), DEVICE , DIMENSION(:,:,:) , ALLOCATABLE :: d_zwrkv, d_zwrkh, d_zwrks, d_zwrk ! 3D GPU workspace
+   REAL(8), DEVICE :: ztmpv, ztmph, ztmps, ztmp ! Device Reduction
+   !
+   INTEGER :: globsize ! 3D workspace size
+   type(dim3) :: dimGrid, dimBlock ! cuda parameters
+   INTEGER, parameter :: nstreams = 3 ! Streams Number
+   INTEGER(kind=cuda_stream_kind) :: stream(nstreams), str ! Stream ID
+   !DEBUG
+   !REAL(8) , save , DIMENSION(:,:,:) , ALLOCATABLE :: prev_3d
+   !REAL(8) :: accum
+
+
+
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_hsb( kt, Kbb, Kmm )
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_hsb  ***
+      !!
+      !! ** Purpose: Compute the ocean global heat content, salt content and volume conservation
+      !!
+      !! ** Method : - Compute the deviation of heat content, salt content and volume
+      !!	            at the current time step from their values at nit000
+      !!	            - Compute the contribution of forcing and remove it from these deviations
+      !!
+      !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt         ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
+      !
+      INTEGER, VALUE                 :: ji, jj, jk, kts ! dummy loop indice
+      INTEGER, VALUE                 :: localsize ! jpi * jpj * jpk
+      INTEGER                        :: istat ! CUDA error check
+      COMPLEX                        :: ctmp ! dummy complex number
+      INTEGER(kind=cuda_stream_kind) :: str ! dummy kernel stream
+      INTEGER                        :: tile_n, tile_b ! tile indexe. _n now, _b before
+      REAL(wp) , DIMENSION(2), SAVE  :: zdiff_hc1, zdiff_sc1 ! heat and salt content variations
+      REAL(wp) , DIMENSION(2), SAVE  :: zdiff_hc, zdiff_sc ! - - - -
+      REAL(wp) , DIMENSION(2), SAVE  :: zdiff_v2 ! volume variation
+      REAL(wp) , DIMENSION(2), SAVE  :: zdiff_v1 ! volume variation
+      REAL(wp) , DIMENSION(2), SAVE  :: zerr_hc1, zerr_sc1 ! heat and salt content misfit
+      REAL(wp) , DIMENSION(2), SAVE  :: zvol_tot ! volume
+      REAL(wp) , DIMENSION(2), SAVE  :: z_frc_trd_t, z_frc_trd_s ! - -
+      REAL(wp) , DIMENSION(2), SAVE  :: z_frc_trd_v ! - -
+      REAL(wp) , DIMENSION(2), SAVE  :: z_wn_trd_t, z_wn_trd_s ! - -
+      REAL(wp) , DIMENSION(2), SAVE  :: z_ssh_hc, z_ssh_sc ! - -
+# 147 "diahsb_new.F90"
+      REAL(wp), DIMENSION(jpi,jpj) :: z2d0, z2d1 ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwrk ! 3D workspace
+      !!---------------------------------------------------------------------------
+      IF( ln_timing )   CALL timing_start('dia_hsb')
+
+      localsize = jpi * jpj * jpk
+      kts = kt
+      IF (kts == 1) THEN
+          tile_n = 1
+          tile_b = 1
+      ELSE
+          IF( MOD(kts,2) == 0) THEN
+              tile_n = 2
+              tile_b = 1
+          ELSE IF( MOD(kts,2) == 1 ) THEN
+              tile_n = 1
+              tile_b = 2
+          END IF
+      END IF
+
+      !
+      ts(:,:,:,1,Kmm) = ts(:,:,:,1,Kmm) * tmask(:,:,:) ; ts(:,:,:,1,Kbb) = ts(:,:,:,1,Kbb) * tmask(:,:,:) ;
+      ts(:,:,:,2,Kmm) = ts(:,:,:,2,Kmm) * tmask(:,:,:) ; ts(:,:,:,2,Kbb) = ts(:,:,:,2,Kbb) * tmask(:,:,:) ;
+      ! ------------------------- !
+      ! 1 - Trends due to forcing !
+      ! ------------------------- !
+      z_frc_trd_v(tile_n) = r1_rho0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * surf(:,:) )! volume fluxes
+      z_frc_trd_t(tile_n) = glob_sum( 'diahsb', sbc_tsc(:,:,jp_tem) * surf(:,:) ) ! heat fluxes
+      z_frc_trd_s(tile_n) = glob_sum( 'diahsb', sbc_tsc(:,:,jp_sal) * surf(:,:) ) ! salt fluxes
+      !                    !  Add runoff    heat & salt input
+      IF( ln_rnf    )   z_frc_trd_t(tile_n) = z_frc_trd_t(tile_n) + glob_sum( 'diahsb', rnf_tsc(:,:,jp_tem) * surf(:,:) )
+      IF( ln_rnf_sal)   z_frc_trd_s(tile_n) = z_frc_trd_s(tile_n) + glob_sum( 'diahsb', rnf_tsc(:,:,jp_sal) * surf(:,:) ) ! Add ice shelf heat & salt input
+      !                    ! Add ice shelf heat & salt input
+      IF( ln_isf    )   z_frc_trd_t(tile_n) = z_frc_trd_t(tile_n) &
+         &                          + glob_sum( 'diahsb', ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) ) ! Add penetrative solar radiation
+      !                    ! Add penetrative solar radiation
+      IF( ln_traqsr )   z_frc_trd_t(tile_n) = z_frc_trd_t(tile_n) + r1_rho0_rcp * glob_sum( 'diahsb', qsr (:,:) * surf(:,:) ) ! Add geothermal heat flux
+      !                    ! Add geothermal heat flux
+      IF( ln_trabbc )   z_frc_trd_t(tile_n) = z_frc_trd_t(tile_n) + glob_sum( 'diahsb', qgh_trd0(:,:) * surf(:,:) )
+      !
+      IF( ln_linssh ) THEN
+         IF( ln_isfcav ) THEN
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  z2d0(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb)
+                  z2d1(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb)
+               END DO
+            END DO
+         ELSE
+            z2d0(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_tem,Kbb)
+            z2d1(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_sal,Kbb)
+         END IF
+         z_wn_trd_t(tile_n) = - glob_sum( 'diahsb', z2d0 )
+         z_wn_trd_s(tile_n) = - glob_sum( 'diahsb', z2d1 )
+      ENDIF
+
+      IF (kts>1) THEN
+         frc_v(tile_n) = frc_v(tile_b)
+         frc_t(tile_n) = frc_t(tile_b)
+         frc_s(tile_n) = frc_s(tile_b)
+         frc_wn_t(tile_n) = frc_wn_t(tile_b)
+         frc_wn_s(tile_n) = frc_wn_s(tile_b)
+      END IF
+      frc_v(tile_n) = frc_v(tile_n) + z_frc_trd_v(tile_n) * rn_Dt
+      frc_t(tile_n) = frc_t(tile_n) + z_frc_trd_t(tile_n) * rn_Dt
+      frc_s(tile_n) = frc_s(tile_n) + z_frc_trd_s(tile_n) * rn_Dt
+      ! ! Advection flux through fixed surface (z=0)
+      IF( ln_linssh ) THEN
+         frc_wn_t(tile_n) = frc_wn_t(tile_n) + z_wn_trd_t(tile_n) * rn_Dt
+         frc_wn_s(tile_n) = frc_wn_s(tile_n) + z_wn_trd_s(tile_n) * rn_Dt
+      ENDIF
+
+      ! ------------------------ !
+      ! 2 -  Content variations  !
+      ! ------------------------ !
+      ! glob_sum_full is needed because you keep the full interior domain to compute the sum (iscpl)
+
+      !                    ! volume variation (calculated with ssh)
+
+      zdiff_v1(tile_n) = glob_sum_full( 'diahsb', surf(:,:)*ssh(:,:,Kmm) - surf_ini(:,:)*ssh_ini(:,:) )
+
+      !                    ! heat & salt content variation (associated with ssh)
+      IF( ln_linssh ) THEN       ! linear free surface case
+         IF( ln_isfcav ) THEN          ! ISF case
+            DO ji = 1, jpi
+               DO jj = 1, jpj
+                  z2d0(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) )
+                  z2d1(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) )
+               END DO
+            END DO
+         ELSE                          ! no under ice-shelf seas
+            z2d0(:,:) = surf(:,:) * ( ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) - ssh_hc_loc_ini(:,:) )
+            z2d1(:,:) = surf(:,:) * ( ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) - ssh_sc_loc_ini(:,:) )
+         END IF
+         z_ssh_hc(tile_n) = glob_sum_full( 'diahsb', z2d0 )
+         z_ssh_sc(tile_n) = glob_sum_full( 'diahsb', z2d1 )
+      ENDIF
+
+      str = stream(tile_n)
+      istat = 0
+      istat = cudaMemcpyAsync( d_e3t, e3t , jpi*jpj*jpk*jpt , str ) + istat
+      istat = cudaMemcpyAsync(d_hc_loc_ini, hc_loc_ini, jpi*jpj*jpk , str ) + istat
+      istat = cudaMemcpyAsync(d_sc_loc_ini, sc_loc_ini, jpi*jpj*jpk , str ) + istat
+      istat = cudaMemcpyAsync( d_ts, ts , jpi*jpj*jpk*2*jpt, str ) + istat
+      IF( istat /= 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to async GPU copy H2D' ) ; RETURN
+      ENDIF
+      dimBlock = dim3(4,4,4)
+      dimGrid = dim3( ceiling( real( jpi ) / dimBlock%x ) , ceiling( real( jpj ) / dimBlock%y ) , &
+                                                             ceiling( real( jpkm1 ) / dimBlock%z ) )
+      !
+      CALL dia_hsb_kernel<<<dimGrid, dimBlock, 0, str>>> (d_surf , d_e3t, d_surf_ini, d_e3t_ini, &
+      & d_ts, d_hc_loc_ini, d_sc_loc_ini, d_tmask, d_tmask_ini, d_zwrkv, d_zwrkh, d_zwrks, d_zwrk, jpi, jpj, jpk, jpt, Kmm)
+
+
+      CALL filter_cuda<<<dimGrid, dimBlock, 0, str>>>(d_zwrkv , d_tmask_h , jpi, jpj, jpk)
+      CALL filter_cuda<<<dimGrid, dimBlock, 0, str>>>(d_zwrkh , d_tmask_h , jpi, jpj, jpk)
+      CALL filter_cuda<<<dimGrid, dimBlock, 0, str>>>(d_zwrks , d_tmask_h , jpi, jpj, jpk)
+      CALL filter_cuda<<<dimGrid, dimBlock, 0, str>>>(d_zwrk , d_tmask_h , jpi, jpj, jpk)
+
+      ztmpv = 0.e0
+      ztmph = 0.e0
+      ztmps = 0.e0
+      ztmp = 0.e0
+
+      istat = 0
+      !$cuf kernel do <<< *, *, stream=str >>>
+      do ji = 1, localsize
+        ztmpv = ztmpv + d_zwrkv(ji)
+      end do
+      istat = cudaMemcpyAsync( h_ztmpv(tile_n) , ztmpv , 1 , str ) + istat
+      !$cuf kernel do <<< *, *, stream=str >>>
+      do ji = 1, localsize
+        ztmph = ztmph + d_zwrkh(ji)
+      end do
+      istat = cudaMemcpyAsync( h_ztmph(tile_n) , ztmph , 1 , str ) + istat
+      !$cuf kernel do <<< *, *, stream=str >>>
+      do ji = 1, localsize
+        ztmps = ztmps + d_zwrks(ji)
+      end do
+      istat = cudaMemcpyAsync( h_ztmps(tile_n) , ztmps , 1 , str ) + istat
+      !$cuf kernel do <<< *, *, stream=str >>>
+      do ji = 1, localsize
+        ztmp = ztmp + d_zwrk(ji)
+      end do
+      istat = cudaMemcpyAsync( h_ztmp (tile_n) , ztmp , 1 , str ) + istat
+      !
+      IF( istat /= 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to async GPU copy D2H' ) ; RETURN
+      ENDIF
+      !
+      istat = cudaStreamSynchronize(stream(tile_b))
+      !
+      IF( istat /= 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to stream synchronize' ) ; RETURN
+      ENDIF
+      !
+      ctmp = CMPLX( h_ztmpv(tile_b) , 0.e0, 8 )
+      CALL mpp_sum('diahsb', ctmp )
+      zdiff_v2(tile_b) = REAL( ctmp, 8 )
+
+      ctmp = CMPLX( h_ztmph(tile_b) , 0.e0, 8 )
+      CALL mpp_sum('diahsb', ctmp )
+      zdiff_hc(tile_b) = REAL( ctmp, 8 )
+
+      ctmp = CMPLX( h_ztmps(tile_b) , 0.e0, 8 )
+      CALL mpp_sum('diahsb', ctmp )
+      zdiff_sc(tile_b) = REAL( ctmp, 8 )
+
+      ctmp = CMPLX( h_ztmp(tile_b) , 0.e0, 8 )
+      CALL mpp_sum('diahsb', ctmp )
+      zvol_tot(tile_b) = REAL( ctmp, 8 )
+
+      IF ( kt == nitend ) THEN
+         !
+         istat = cudaStreamSynchronize(stream(tile_n))
+         IF( istat /= 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to stream synchronize' ) ; RETURN
+         ENDIF
+         !
+         ctmp = CMPLX( h_ztmpv(tile_n) , 0.e0, 8 )
+         CALL mpp_sum('diahsb', ctmp )
+         zdiff_v2(tile_n) = REAL( ctmp, 8 )
+
+         ctmp = CMPLX( h_ztmph(tile_n) , 0.e0, 8 )
+         CALL mpp_sum('diahsb', ctmp )
+         zdiff_hc(tile_n) = REAL( ctmp, 8 )
+
+         ctmp = CMPLX( h_ztmps(tile_n) , 0.e0, 8 )
+         CALL mpp_sum('diahsb', ctmp )
+         zdiff_sc(tile_n) = REAL( ctmp, 8 )
+
+         ctmp = CMPLX( h_ztmp(tile_n) , 0.e0, 8 )
+         CALL mpp_sum('diahsb', ctmp )
+         zvol_tot(tile_n) = REAL( ctmp, 8 )
+      ENDIF
+      ! ------------------------ !
+      ! 3 - Drifts !
+      ! ------------------------ !
+
+      IF ( kt > 1 ) THEN
+            kts = kts - 1
+            zdiff_v1(tile_b) = zdiff_v1(tile_b) - frc_v(tile_b)
+            IF( .NOT.ln_linssh ) zdiff_v2(tile_b) = zdiff_v2(tile_b) - frc_v(tile_b)
+            zdiff_hc(tile_b) = zdiff_hc(tile_b) - frc_t(tile_b)
+            zdiff_sc(tile_b) = zdiff_sc(tile_b) - frc_s(tile_b)
+            IF( ln_linssh ) THEN
+                zdiff_hc1(tile_b) = zdiff_hc (tile_b) + z_ssh_hc(tile_b)
+                zdiff_sc1(tile_b) = zdiff_sc (tile_b) + z_ssh_sc(tile_b)
+                zerr_hc1 (tile_b) = z_ssh_hc(tile_b) - frc_wn_t(tile_b)
+                zerr_sc1 (tile_b) = z_ssh_sc(tile_b) - frc_wn_s(tile_b)
+            ENDIF
+            !!gm to be added ?
+            ! IF( ln_linssh ) THEN ! fixed volume, add the ssh contribution
+            ! zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * sshn(:,:) )
+            ! ENDIF
+            !!gm end
+
+            CALL iom_put( 'bgfrcvol' , frc_v(tile_b) * 1.e-9 ) ! vol - surface forcing (km3)
+            CALL iom_put( 'bgfrctem' , frc_t(tile_b) * rho0 * rcp * 1.e-20 ) ! hc - surface forcing (1.e20 J)
+            CALL iom_put( 'bgfrchfx' , frc_t(tile_b) * rho0 * rcp / & ! hc - surface forcing (W/m2)
+                & ( surf_tot * kts * rn_Dt ) )
+            CALL iom_put( 'bgfrcsal' , frc_s(tile_b) * 1.e-9 ) ! sc - surface forcing (psu*km3)
+            IF( .NOT. ln_linssh ) THEN
+                CALL iom_put( 'bgtemper' , zdiff_hc(tile_b) / zvol_tot(tile_b) ) ! Temperature drift (C)
+                CALL iom_put( 'bgsaline' , zdiff_sc(tile_b) / zvol_tot(tile_b) ) ! Salinity drift (PSU)
+                CALL iom_put( 'bgheatco' , zdiff_hc(tile_b) * 1.e-20 * rho0 * rcp ) ! Heat content drift (1.e20 J)
+                CALL iom_put( 'bgheatfx' , zdiff_hc(tile_b) * rho0 * rcp / & ! Heat flux drift (W/m2)
+                    & ( surf_tot * kts * rn_Dt ) )
+                CALL iom_put( 'bgsaltco' , zdiff_sc(tile_b) * 1.e-9 ) ! Salt content drift (psu*km3)
+                CALL iom_put( 'bgvolssh' , zdiff_v1(tile_b) * 1.e-9 ) ! volume ssh drift (km3)
+                CALL iom_put( 'bgvole3t' , zdiff_v2(tile_b) * 1.e-9 ) ! volume e3t drift (km3)
+                !
+! IF( lwp ) THEN
+! WRITE(numout,*)
+! WRITE(numout,*) 'dia_hsb : last time step hsb diagnostics: at it= ', kt,' date= ', ndastp
+! WRITE(numout,*) '~~~~~~~'
+! WRITE(numout,*) '   Temperature drift = ', zdiff_hc(tile_b) / zvol_tot(tile_b), ' C'
+! WRITE(numout,*) '   Salinity12  drift = ', zdiff_sc(tile_b) / zvol_tot(tile_b), ' PSU'
+! WRITE(numout,*) '   volume ssh  drift = ', zdiff_v1(tile_b) * 1.e-9 , ' km^3'
+! WRITE(numout,*) '   volume e3t  drift = ', zdiff_v2(tile_b) * 1.e-9 , ' km^3'
+! ENDIF
+            ELSE
+                CALL iom_put( 'bgtemper' , zdiff_hc1(tile_b) / zvol_tot(tile_b)) ! Heat content drift (C)
+                CALL iom_put( 'bgsaline' , zdiff_sc1(tile_b) / zvol_tot(tile_b)) ! Salt content drift (PSU)
+                CALL iom_put( 'bgheatco' , zdiff_hc1(tile_b) * 1.e-20 * rho0 * rcp ) ! Heat content drift (1.e20 J)
+                CALL iom_put( 'bgheatfx' , zdiff_hc1(tile_b) * rho0 * rcp / & ! Heat flux drift (W/m2)
+                    & ( surf_tot * kts * rn_Dt ) )
+                CALL iom_put( 'bgsaltco' , zdiff_sc1(tile_b) * 1.e-9 ) ! Salt content drift (psu*km3)
+                CALL iom_put( 'bgvolssh' , zdiff_v1(tile_b) * 1.e-9 ) ! volume ssh drift (km3)
+                CALL iom_put( 'bgmistem' , zerr_hc1(tile_b) / zvol_tot(tile_b) ) ! hc - error due to free surface (C)
+                CALL iom_put( 'bgmissal' , zerr_sc1(tile_b) / zvol_tot(tile_b) ) ! sc - error due to free surface (psu)
+            ENDIF
+            !
+            IF( lrst_oce ) CALL dia_hsb_rst( kts, Kmm, tile_n, 'WRITE' )
+        !
+        END IF
+        IF ( kt == nitend ) THEN
+
+            zdiff_v1(tile_n) = zdiff_v1(tile_n) - frc_v(tile_n)
+            IF( .NOT.ln_linssh ) zdiff_v2(tile_n) = zdiff_v2(tile_n) - frc_v(tile_n)
+            zdiff_hc(tile_n) = zdiff_hc(tile_n) - frc_t(tile_n)
+            zdiff_sc(tile_n) = zdiff_sc(tile_n) - frc_s(tile_n)
+            IF( ln_linssh ) THEN
+                zdiff_hc1(tile_n) = zdiff_hc (tile_n) + z_ssh_hc(tile_n)
+                zdiff_sc1(tile_n) = zdiff_sc (tile_n) + z_ssh_sc(tile_n)
+                zerr_hc1 (tile_n) = z_ssh_hc(tile_n) - frc_wn_t(tile_n)
+                zerr_sc1 (tile_n) = z_ssh_sc(tile_n) - frc_wn_s(tile_n)
+            ENDIF
+            !!gm to be added ?
+            ! IF( ln_linssh ) THEN ! fixed volume, add the ssh contribution
+            ! zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * sshn(:,:) )
+            ! ENDIF
+            !!gm end
+
+            CALL iom_put( 'bgfrcvol' , frc_v(tile_n) * 1.e-9 ) ! vol - surface forcing (km3)
+            CALL iom_put( 'bgfrctem' , frc_t(tile_n) * rho0 * rcp * 1.e-20 ) ! hc - surface forcing (1.e20 J)
+            CALL iom_put( 'bgfrchfx' , frc_t(tile_n) * rho0 * rcp / & ! hc - surface forcing (W/m2)
+                & ( surf_tot * kt * rn_Dt ) )
+            CALL iom_put( 'bgfrcsal' , frc_s(tile_n) * 1.e-9 ) ! sc - surface forcing (psu*km3)
+            IF( .NOT. ln_linssh ) THEN
+                CALL iom_put( 'bgtemper' , zdiff_hc(tile_n) / zvol_tot(tile_n) ) ! Temperature drift (C)
+                CALL iom_put( 'bgsaline' , zdiff_sc(tile_n) / zvol_tot(tile_n) ) ! Salinity drift (PSU)
+                CALL iom_put( 'bgheatco' , zdiff_hc(tile_n) * 1.e-20 * rho0 * rcp ) ! Heat content drift (1.e20 J)
+                CALL iom_put( 'bgheatfx' , zdiff_hc(tile_n) * rho0 * rcp / & ! Heat flux drift (W/m2)
+                    & ( surf_tot * kt * rn_Dt ) )
+                CALL iom_put( 'bgsaltco' , zdiff_sc(tile_n) * 1.e-9 ) ! Salt content drift (psu*km3)
+                CALL iom_put( 'bgvolssh' , zdiff_v1(tile_n) * 1.e-9 ) ! volume ssh drift (km3)
+                CALL iom_put( 'bgvole3t' , zdiff_v2(tile_n) * 1.e-9 ) ! volume e3t drift (km3)
+                !
+                IF( kt == nitend .AND. lwp ) THEN
+                    WRITE(numout,*)
+                    WRITE(numout,*) 'dia_hsb : last time step hsb diagnostics: at it= ', kt,' date= ', ndastp
+                    WRITE(numout,*) '~~~~~~~'
+                    WRITE(numout,*) '   Temperature drift = ', zdiff_hc(tile_n) / zvol_tot(tile_n), ' C'
+                    WRITE(numout,*) '   Salinity  drift   = ', zdiff_sc(tile_n) / zvol_tot(tile_n), ' PSU'
+                    WRITE(numout,*) '   volume ssh  drift = ', zdiff_v1(tile_n) * 1.e-9 , ' km^3'
+                    WRITE(numout,*) '   volume e3t  drift = ', zdiff_v2(tile_n) * 1.e-9 , ' km^3'
+                !
+                ENDIF
+               !
+            ELSE
+                CALL iom_put( 'bgtemper' , zdiff_hc1(tile_n) / zvol_tot(tile_n)) ! Heat content drift (C)
+                CALL iom_put( 'bgsaline' , zdiff_sc1(tile_n) / zvol_tot(tile_n)) ! Salt content drift (PSU)
+                CALL iom_put( 'bgheatco' , zdiff_hc1(tile_n) * 1.e-20 * rho0 * rcp ) ! Heat content drift (1.e20 J)
+                CALL iom_put( 'bgheatfx' , zdiff_hc1(tile_n) * rho0 * rcp / & ! Heat flux drift (W/m2)
+                    & ( surf_tot * kt * rn_Dt ) )
+                CALL iom_put( 'bgsaltco' , zdiff_sc1(tile_n) * 1.e-9 ) ! Salt content drift (psu*km3)
+                CALL iom_put( 'bgvolssh' , zdiff_v1(tile_n) * 1.e-9 ) ! volume ssh drift (km3)
+                CALL iom_put( 'bgmistem' , zerr_hc1(tile_n) / zvol_tot(tile_n) ) ! hc - error due to free surface (C)
+                CALL iom_put( 'bgmissal' , zerr_sc1(tile_n) / zvol_tot(tile_n) ) ! sc - error due to free surface (psu)
+            ENDIF
+            !
+            !Last step, don't need restart
+            !IF( lrst_oce ) CALL dia_hsb_rst( kts, Kmm, tile_n, 'WRITE' )
+        !
+        END IF
+      IF( ln_timing ) CALL timing_stop('dia_hsb')
+      !
+   END SUBROUTINE dia_hsb
+
+
+   SUBROUTINE dia_hsb_rst( kt, Kmm, tile, cdrw )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE dia_hsb_rst  ***
+      !!
+      !! ** Purpose : Read or write DIA file in restart file
+      !!
+      !! ** Method  : use of IOM library
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      INTEGER         , INTENT(in) ::   Kmm    ! ocean time level index
+      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      INTEGER , INTENT(in)         :: tile     ! host tile
+      !
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            !
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   dia_hsb_rst : read hsb restart at it= ', kt,' date= ', ndastp
+            IF(lwp) WRITE(numout,*)
+            CALL iom_get( numror, 'frc_v', frc_v(tile) )
+            CALL iom_get( numror, 'frc_t', frc_t(tile) )
+            CALL iom_get( numror, 'frc_s', frc_s(tile) )
+            IF( ln_linssh ) THEN
+               CALL iom_get( numror, 'frc_wn_t', frc_wn_t(tile) )
+               CALL iom_get( numror, 'frc_wn_s', frc_wn_s(tile) )
+            ENDIF
+            CALL iom_get( numror, jpdom_auto, 'surf_ini'  , surf_ini   ) ! ice sheet coupling
+            CALL iom_get( numror, jpdom_auto, 'ssh_ini'   , ssh_ini    )
+            CALL iom_get( numror, jpdom_auto, 'e3t_ini'   , e3t_ini    )
+            CALL iom_get( numror, jpdom_auto, 'tmask_ini' , tmask_ini  )
+            CALL iom_get( numror, jpdom_auto, 'hc_loc_ini', hc_loc_ini )
+            CALL iom_get( numror, jpdom_auto, 'sc_loc_ini', sc_loc_ini )
+            IF( ln_linssh ) THEN
+               CALL iom_get( numror, jpdom_auto, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
+               CALL iom_get( numror, jpdom_auto, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
+            ENDIF
+         ELSE
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   dia_hsb_rst : initialise hsb at initial state '
+            IF(lwp) WRITE(numout,*)
+            surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:)         ! initial ocean surface
+            ssh_ini(:,:) = ssh(:,:,Kmm)                          ! initial ssh
+            DO jk = 1, jpk
+              ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
+               e3t_ini   (:,:,jk) = e3t(:,:,jk,Kmm)                      * tmask(:,:,jk)  ! initial vertical scale factors
+               tmask_ini (:,:,jk) = tmask(:,:,jk)                                       ! initial mask
+               hc_loc_ini(:,:,jk) = ts(:,:,jk,jp_tem,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial heat content
+               sc_loc_ini(:,:,jk) = ts(:,:,jk,jp_sal,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial salt content
+            END DO
+            d_surf_ini = surf_ini
+            d_e3t_ini = e3t_ini
+            d_tmask_ini = tmask_ini
+            d_hc_loc_ini = hc_loc_ini
+            d_sc_loc_ini = sc_loc_ini
+            frc_v(tile) = 0._wp                                           ! volume       trend due to forcing
+            frc_t(tile) = 0._wp                                           ! heat content   -    -   - -
+            frc_s(tile) = 0._wp                                           ! salt content   -    -   - -
+            IF( ln_linssh ) THEN
+               IF( ln_isfcav ) THEN
+                  DO ji = 1, jpi
+                     DO jj = 1, jpj
+                        ssh_hc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) ! initial heat content in ssh
+                        ssh_sc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) ! initial salt content in ssh
+                     END DO
+                   END DO
+                ELSE
+                  ssh_hc_loc_ini(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) ! initial heat content in ssh
+                  ssh_sc_loc_ini(:,:) = ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) ! initial salt content in ssh
+               END IF
+               frc_wn_t(tile) = 0._wp ! initial heat content misfit due to free surface
+               frc_wn_s(tile) = 0._wp ! initial salt content misfit due to free surface
+            ENDIF
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   dia_hsb_rst : write restart at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*)
+         !
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_v', frc_v(tile) )
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_t', frc_t(tile) )
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_s', frc_s(tile) )
+         IF( ln_linssh ) THEN
+            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_t', frc_wn_t(tile) )
+            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s(tile) )
+         ENDIF
+         CALL iom_rstput( kt, nitrst, numrow, 'surf_ini'  , surf_ini   ) ! ice sheet coupling
+         CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini'   , ssh_ini    )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini'   , e3t_ini    )
+         CALL iom_rstput( kt, nitrst, numrow, 'tmask_ini' , tmask_ini  )
+         CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini )
+         CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini )
+         IF( ln_linssh ) THEN
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE dia_hsb_rst
+
+
+   SUBROUTINE dia_hsb_init( Kmm )
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_hsb  ***
+      !!
+      !! ** Purpose: Initialization for the heat salt volume budgets
+      !!
+      !! ** Method : Compute initial heat content, salt content and volume
+      !!
+      !! ** Action : - Compute initial heat content, salt content and volume
+      !!             - Initialize forcing trends
+      !!             - Compute coefficients for conversion
+      !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kmm ! time level index
+      !
+      INTEGER ::   ierror, ios   ! local integer
+      INTEGER ::   i, istat      ! local integer
+      !!
+      NAMELIST/namhsb/ ln_diahsb
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_hsb_init : heat and salt budgets diagnostics'
+         WRITE(numout,*) '~~~~~~~~~~~~ '
+      ENDIF
+      READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namhsb in reference namelist' )
+      READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namhsb in configuration namelist' )
+      IF(lwm) WRITE( numond, namhsb )
+
+      IF(lwp) THEN
+         WRITE(numout,*) '   Namelist  namhsb :'
+         WRITE(numout,*) '      check the heat and salt budgets (T) or not (F)       ln_diahsb = ', ln_diahsb
+      ENDIF
+      !
+      IF( .NOT. ln_diahsb )   RETURN
+
+      ! ------------------- !
+      ! 1 - Allocate memory !
+      ! ------------------- !
+
+      CALL setdevice()
+      !Device data associate to PUBLIC arrays
+      ALLOCATE(d_e3t (jpi,jpj,jpk,jpt) ) !
+      ALLOCATE(d_tmask (jpi,jpj,jpk) ) !
+      ALLOCATE(d_tmask_ini (jpi,jpj,jpk) ) !
+      ALLOCATE(d_tmask_h (jpi,jpj) ) !
+      ALLOCATE(d_ts (jpi,jpj,jpk,2,jpj) ) !
+      !Device data associate to LOCAL/DEVICE arrays !
+      ALLOCATE(d_surf (jpi,jpj) ) !
+      ALLOCATE(d_surf_ini (jpi,jpj) ) !
+      ALLOCATE(d_hc_loc_ini (jpi,jpj,jpk) ) !
+      ALLOCATE(d_sc_loc_ini (jpi,jpj,jpk) ) !
+      ALLOCATE(d_e3t_ini (jpi,jpj,jpk) ) !
+      ALLOCATE(d_zwrkv (jpi,jpj,jpkm1) ) !
+      ALLOCATE(d_zwrkh (jpi,jpj,jpkm1) ) !
+      ALLOCATE(d_zwrks (jpi,jpj,jpkm1) ) !
+      ALLOCATE(d_zwrk (jpi,jpj,jpkm1) ) !
+      ALLOCATE(h_ztmpv(2),h_ztmph(2),h_ztmps(2),h_ztmp(2)) !
+
+      DO i = 1, nstreams !Create Streams
+          istat = cudaStreamCreate(stream(i))
+         IF( istat /= 0 ) THEN
+         CALL ctl_stop( 'dia_hsb_init: error in Stream creation' ) ; RETURN
+         ENDIF
+      END DO
+            !
+      !Pinned reallocation step non constant
+      istat = cudaHostRegister(C_LOC(ts ), sizeof(ts ), cudaHostRegisterMapped)
+      istat = cudaHostRegister(C_LOC(e3t), sizeof(e3t), cudaHostRegisterMapped)
+      IF( istat /= 0 ) THEN
+         CALL ctl_stop( 'dia_hsb_init: unable to pin host memory to GPU' ) ; RETURN
+      ENDIF
+
+      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), &
+         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), tmask_ini(jpi,jpj,jpk),STAT=ierror  )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb_init: unable to allocate hc_loc_ini' )   ;   RETURN
+      ENDIF
+
+      IF( ln_linssh )   ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to allocate ssh_hc_loc_ini' )   ;   RETURN
+      ENDIF
+
+      ! ----------------------------------------------- !
+      ! 2 - Time independant variables and file opening !
+      ! ----------------------------------------------- !
+      surf(:,:) = e1e2t(:,:) * tmask_i(:,:)               ! masked surface grid cell area
+      surf_tot  = glob_sum( 'diahsb', surf(:,:) )         ! total ocean surface area
+
+       d_surf = surf
+       d_surf_ini = surf_ini
+       d_e3t_ini = e3t_ini
+       d_tmask = tmask
+       d_tmask_ini = tmask_ini
+       d_tmask_h = tmask_h
+       h_ztmp = 0.0
+
+      IF( ln_bdy ) CALL ctl_warn( 'dia_hsb_init: heat/salt budget does not consider open boundary fluxes' )
+      !
+      ! ---------------------------------- !
+      ! 4 - initial conservation variables !
+      ! ---------------------------------- !
+      CALL dia_hsb_rst( nit000, Kmm, 1, 'READ' ) !* read or initialize all required files
+
+
+
+      !
+   END SUBROUTINE dia_hsb_init
+
+   !!======================================================================
+END MODULE diahsb
Index: NEMO/trunk/tests/DIA_GPU/MY_SRC/gpu_manager.F90
===================================================================
--- NEMO/trunk/tests/DIA_GPU/MY_SRC/gpu_manager.F90	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/MY_SRC/gpu_manager.F90	(revision 14854)
@@ -0,0 +1,34 @@
+#ifdef key_gpu
+MODULE gpu_manager
+    USE lib_mpp        ! distributed memory computing library
+    USE cudafor
+    USE par_oce
+    !!======================================================================
+    !!                       ***  MODULE  gpu_manager  ***
+    !! GPU CUDA utilities : Defines run parameters for GPGPU routines
+    !!=====================================================================
+    !! History :   1.0  !  2019-10  (M. Faria)   original code
+    !!----------------------------------------------------------------------
+
+    !!----------------------------------------------------------------------
+    IMPLICIT NONE
+    PUBLIC
+    !!----------------------------------------------------------------------
+    !!                   gpu configuration parameters
+    !!----------------------------------------------------------------------
+    INTEGER          :: cudaistat, usecudev
+    INTEGER          :: nodesize
+    INTEGER          :: cudadev
+    INTEGER          :: cusocket
+CONTAINS
+    SUBROUTINE setdevice()
+        nodesize  = 1                                                    !Node number of MPI processes, CTE-Power default 40
+        cusocket  = 1                                                    !Number of NVlink or PCI-E nodes, CTE-Power default 2
+        cudadev   = 1                                                    !Number of GPUs on a node, CTE-Power default 4
+        usecudev  = MOD( mpprank /  ( nodesize / cusocket )  , cudadev ) !Bind MPI RANK to GPU, consider balanced NVlink and constant number of devices per node 
+        cudaistat = cudaSetDevice( usecudev  )                           !Set which GPU use. Ex for 4 GPUs. cpus_per_socket / sockets_per_node for "bind-to core"
+        print *, 'CUDA istat ', cudaistat, 'rank', mpprank, 'cuda device', cudadev 
+    END SUBROUTINE setdevice
+
+END MODULE gpu_manager
+#endif
Index: NEMO/trunk/tests/DIA_GPU/MY_SRC/nvtx.F90
===================================================================
--- NEMO/trunk/tests/DIA_GPU/MY_SRC/nvtx.F90	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/MY_SRC/nvtx.F90	(revision 14854)
@@ -0,0 +1,66 @@
+#ifdef key_gpu
+module nvtx
+
+use iso_c_binding
+implicit none
+
+integer,private :: col(7) = [ Z'0000ff00', Z'000000ff', Z'00ffff00', Z'00ff00ff', Z'0000ffff', Z'00ff0000', Z'00ffffff']
+character(len=256),private :: tempName
+
+type, bind(C):: nvtxEventAttributes
+  integer(C_INT16_T):: version=1
+  integer(C_INT16_T):: size=48 !
+  integer(C_INT):: category=0
+  integer(C_INT):: colorType=1 ! NVTX_COLOR_ARGB = 1
+  integer(C_INT):: color
+  integer(C_INT):: payloadType=0 ! NVTX_PAYLOAD_UNKNOWN = 0
+  integer(C_INT):: reserved0
+  integer(C_INT64_T):: payload   ! union uint,int,double
+  integer(C_INT):: messageType=1  ! NVTX_MESSAGE_TYPE_ASCII     = 1 
+  type(C_PTR):: message  ! ascii char
+end type
+
+interface nvtxRangePush
+  ! push range with custom label and standard color
+  subroutine nvtxRangePushA(name) bind(C, name='nvtxRangePushA')
+  use iso_c_binding
+  character(kind=C_CHAR,len=*) :: name
+  end subroutine
+
+  ! push range with custom label and custom color
+  subroutine nvtxRangePushEx(event) bind(C, name='nvtxRangePushEx')
+  use iso_c_binding
+  import:: nvtxEventAttributes
+  type(nvtxEventAttributes):: event
+  end subroutine
+end interface
+
+interface nvtxRangePop
+  subroutine nvtxRangePop() bind(C, name='nvtxRangePop')
+  end subroutine
+end interface
+
+contains
+
+subroutine nvtxStartRange(name,id)
+  character(kind=c_char,len=*) :: name
+  integer, optional:: id
+  type(nvtxEventAttributes):: event
+
+  tempName=trim(name)//c_null_char
+
+  if ( .not. present(id)) then
+    call nvtxRangePush(tempName)
+  else
+    event%color=col(mod(id,7)+1)
+    event%message=c_loc(tempName)
+    call nvtxRangePushEx(event)
+  end if
+end subroutine
+
+subroutine nvtxEndRange
+  call nvtxRangePop
+end subroutine
+
+end module nvtx
+#endif
Index: NEMO/trunk/tests/DIA_GPU/MY_SRC/stpctl.F90
===================================================================
--- NEMO/trunk/tests/DIA_GPU/MY_SRC/stpctl.F90	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/MY_SRC/stpctl.F90	(revision 14854)
@@ -0,0 +1,330 @@
+MODULE stpctl
+   !!======================================================================
+   !!                       ***  MODULE  stpctl  ***
+   !! Ocean run control :  gross check of the ocean time stepping
+   !!======================================================================
+   !! History :  OPA  ! 1991-03  (G. Madec) Original code
+   !!            6.0  ! 1992-06  (M. Imbard)
+   !!            8.0  ! 1997-06  (A.M. Treguier)
+   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
+   !!            2.0  ! 2009-07  (G. Madec)  Add statistic for time-spliting
+   !!            3.7  ! 2016-09  (G. Madec)  Remove solver
+   !!            4.0  ! 2017-04  (G. Madec)  regroup global communications
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   stp_ctl      : Control the run
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables 
+   USE c1d             ! 1D vertical configuration
+   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables
+   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy
+   !  
+   USE diawri          ! Standard run outputs       (dia_wri_state routine)
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp         ! distributed memory computing
+   USE eosbn2, ONLY: ln_SEOS, rn_b0
+   !
+   USE netcdf          ! NetCDF library
+   USE, INTRINSIC :: IEEE_ARITHMETIC, ONLY : IEEE_IS_NAN
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC stp_ctl           ! routine called by step.F90
+
+   INTEGER, PARAMETER         ::   jpvar = 8
+   INTEGER                    ::   nrunid   ! netcdf file id
+   INTEGER, DIMENSION(jpvar)  ::   nvarid   ! netcdf variable id
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: stpctl.F90 13616 2020-10-15 15:58:33Z smasson $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE stp_ctl( kt, Kmm )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE stp_ctl  ***
+      !!
+      !! ** Purpose :   Control the run
+      !!
+      !! ** Method  : - Save the time step in numstp
+      !!              - Stop the run IF problem encountered by setting nstop > 0
+      !!                Problems checked: |ssh| maximum larger than 10 m
+      !!                                  |U|   maximum larger than 10 m/s 
+      !!                                  negative sea surface salinity
+      !!
+      !! ** Actions :   "time.step" file = last ocean time-step
+      !!                "run.stat"  file = run statistics
+      !!                 nstop indicator sheared among all local domain
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
+      !!
+      INTEGER, PARAMETER              ::   jptst = 4
+      INTEGER                         ::   ji                                    ! dummy loop indices
+      INTEGER                         ::   idtime, istatus
+      INTEGER , DIMENSION(jptst)      ::   iareasum, iareamin, iareamax
+      INTEGER , DIMENSION(3,jptst)    ::   iloc                                  ! min/max loc indices
+      REAL(wp)                        ::   zzz, zminsal, zmaxsal                 ! local real 
+      REAL(wp), DIMENSION(jpvar+1)    ::   zmax
+      REAL(wp), DIMENSION(jptst)      ::   zmaxlocal
+      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce
+      LOGICAL, DIMENSION(jpi,jpj,jpk) ::   llmsk
+      CHARACTER(len=20)               ::   clname
+      !!----------------------------------------------------------------------
+      IF( nstop > 0 .AND. ngrdstop > -1 )   RETURN   !   stpctl was already called by a child grid
+      !
+      ll_wrtstp  = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
+      ll_colruns = ll_wrtstp .AND. sn_cfctl%l_runstat .AND. jpnij > 1
+      ll_wrtruns = ( ll_colruns .OR. jpnij == 1 ) .AND. lwm
+      !
+      IF( kt == nit000 ) THEN
+         !
+         IF( lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'stp_ctl : time-stepping control'
+            WRITE(numout,*) '~~~~~~~'
+         ENDIF
+         !                                ! open time.step    ascii file, done only by 1st subdomain
+         IF( lwm )   CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+         !
+         IF( ll_wrtruns ) THEN
+            !                             ! open run.stat     ascii file, done only by 1st subdomain
+            CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+            !                             ! open run.stat.nc netcdf file, done only by 1st subdomain
+            clname = 'run.stat.nc'
+            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
+            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, nrunid )
+            istatus = NF90_DEF_DIM( nrunid, 'time', NF90_UNLIMITED, idtime )
+            istatus = NF90_DEF_VAR( nrunid, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), nvarid(1) )
+            istatus = NF90_DEF_VAR( nrunid,   'abs_u_max', NF90_DOUBLE, (/ idtime /), nvarid(2) )
+            istatus = NF90_DEF_VAR( nrunid,       's_min', NF90_DOUBLE, (/ idtime /), nvarid(3) )
+            istatus = NF90_DEF_VAR( nrunid,       's_max', NF90_DOUBLE, (/ idtime /), nvarid(4) )
+            istatus = NF90_DEF_VAR( nrunid,       't_min', NF90_DOUBLE, (/ idtime /), nvarid(5) )
+            istatus = NF90_DEF_VAR( nrunid,       't_max', NF90_DOUBLE, (/ idtime /), nvarid(6) )
+            IF( ln_zad_Aimp ) THEN
+               istatus = NF90_DEF_VAR( nrunid,   'Cf_max', NF90_DOUBLE, (/ idtime /), nvarid(7) )
+               istatus = NF90_DEF_VAR( nrunid,'abs_wi_max',NF90_DOUBLE, (/ idtime /), nvarid(8) )
+            ENDIF
+            istatus = NF90_ENDDEF(nrunid)
+         ENDIF
+         !
+      ENDIF
+      !
+      !                                   !==              write current time step              ==!
+      !                                   !==  done only by 1st subdomain at writting timestep  ==!
+      IF( lwm .AND. ll_wrtstp ) THEN
+         WRITE ( numstp, '(1x, i8)' )   kt
+         REWIND( numstp )
+      ENDIF
+      !                                   !==            test of local extrema           ==!
+      !                                   !==  done by all processes at every time step  ==!
+      !
+      llmsk(     1:nn_hls,:,:) = .FALSE.                                          ! exclude halos from the checked region
+      llmsk(Nie0+1:   jpi,:,:) = .FALSE.
+      llmsk(:,     1:nn_hls,:) = .FALSE.
+      llmsk(:,Nje0+1:   jpj,:) = .FALSE.
+      !
+      llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0) == 1._wp         ! define only the inner domain
+      !
+      ll_0oce = .NOT. ANY( llmsk(:,:,1) )                                         ! no ocean point in the inner domain?
+      !
+      IF( ll_wd ) THEN
+         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max
+      ELSE
+         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm)           ), mask = llmsk(:,:,1) )   ! ssh max
+      ENDIF
+      llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) ), mask = llmsk )                     ! velocity max (zonal only)
+      llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+      zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     ! minus salinity max
+      zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     !       salinity max
+      IF( ll_colruns .OR. jpnij == 1 ) THEN     ! following variables are used only in the netcdf file
+         zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  ! minus temperature max
+         zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  !       temperature max
+         IF( ln_zad_Aimp ) THEN
+            zmax(7) = MAXVAL(   Cu_adv(:,:,:)   , mask = llmsk )                  ! partitioning coeff. max
+            llmsk(:,:,:) = wmask(:,:,:) == 1._wp
+            zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = llmsk )                  ! implicit vertical vel. max
+         ELSE
+            zmax(7:8) = 0._wp
+         ENDIF
+      ELSE
+         zmax(5:8) = 0._wp
+      ENDIF
+      zmax(jpvar+1) = REAL( nstop, wp )                                           ! stop indicator
+      !
+      !                                   !==               get global extrema             ==!
+      !                                   !==  done by all processes if writting run.stat  ==!
+      IF( ll_colruns ) THEN
+         zmaxlocal(:) = zmax(1:jptst)
+         CALL mpp_max( "stpctl", zmax )          ! max over the global domain: ok even of ll_0oce = .true. 
+         nstop = NINT( zmax(jpvar+1) )           ! update nstop indicator (now sheared among all local domains)
+      ELSE
+         ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow.
+         IF( ll_0oce )   zmax(1:jptst) = (/ 0._wp, 0._wp, -1._wp, 1._wp /)   ! default "valid" values...
+      ENDIF
+      !
+      zmax(3) = -zmax(3)                              ! move back from max(-zz) to min(zz) : easier to manage! 
+      zmax(5) = -zmax(5)                              ! move back from max(-zz) to min(zz) : easier to manage!
+      IF( ll_colruns ) zmaxlocal(3) = -zmaxlocal(3)   ! move back from max(-zz) to min(zz) : easier to manage!
+      !
+      !                                   !==              write "run.stat" files              ==!
+      !                                   !==  done only by 1st subdomain at writting timestep  ==!
+      IF( ll_wrtruns ) THEN
+         WRITE(numrun,9500) kt, zmax(1:jptst)
+         DO ji = 1, jpvar - 2 * COUNT( .NOT. (/ln_zad_Aimp/) )
+            istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) )
+         END DO
+         IF( kt == nitend )   istatus = NF90_CLOSE(nrunid)
+      END IF
+      !                                   !==               error handling               ==!
+      !                                   !==  done by all processes at every time step  ==!
+      !
+      IF ( ln_SEOS.AND.(rn_b0==0._wp) ) THEN             ! Discard checks on salinity
+         zmaxsal =  HUGE(1._wp)                               ! if not used in eos
+         zminsal = -HUGE(1._wp)
+      ELSE
+         zmaxsal = 100._wp
+         zminsal =   0._wp
+      ENDIF 
+      ! 
+      IF(  zmax(1) >   20._wp .OR.   &                        ! too large sea surface height ( > 20 m )
+         & zmax(2) >   10._wp .OR.   &                        ! too large velocity ( > 10 m/s)
+         & zmax(3) <= zminsal .OR.   &                        ! negative or zero sea surface salinity
+         & zmax(4) >= zmaxsal .OR.   &                        ! too large sea surface salinity ( > 100 )
+         & zmax(4) <  zminsal .OR.   &                        ! too large sea surface salinity (keep this line for sea-ice)
+         & IEEE_IS_NAN( SUM(zmax(1:jptst)) ) .OR.   &         ! NaN encounter in the tests
+         & ABS(   SUM(zmax(1:jptst)) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests
+         !
+         iloc(:,:) = 0
+         IF( ll_colruns ) THEN   ! zmax is global, so it is the same on all subdomains -> no dead lock with mpp_maxloc
+            ! first: close the netcdf file, so we can read it
+            IF( lwm .AND. kt /= nitend )   istatus = NF90_CLOSE(nrunid)
+            ! get global loc on the min/max
+            llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp         ! define only the inner domain
+            CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,         Kmm)), llmsk(:,:,1), zzz, iloc(1:2,1) )   ! mpp_maxloc ok if mask = F 
+            llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+            CALL mpp_maxloc( 'stpctl', ABS( uu(:,:,:,       Kmm)), llmsk(:,:,:), zzz, iloc(1:3,2) )
+            llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+            CALL mpp_minloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) , llmsk(:,:,:), zzz, iloc(1:3,3) )
+            CALL mpp_maxloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) , llmsk(:,:,:), zzz, iloc(1:3,4) )
+            ! find which subdomain has the max.
+            iareamin(:) = jpnij+1   ;   iareamax(:) = 0   ;   iareasum(:) = 0
+            DO ji = 1, jptst
+               IF( zmaxlocal(ji) == zmax(ji) ) THEN
+                  iareamin(ji) = narea   ;   iareamax(ji) = narea   ;   iareasum(ji) = 1
+               ENDIF
+            END DO
+            CALL mpp_min( "stpctl", iareamin )         ! min over the global domain
+            CALL mpp_max( "stpctl", iareamax )         ! max over the global domain
+            CALL mpp_sum( "stpctl", iareasum )         ! sum over the global domain
+         ELSE                    ! find local min and max locations:
+            ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc
+            llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp        ! define only the inner domain
+            iloc(1:2,1) = MAXLOC( ABS( ssh(:,:,         Kmm)), mask = llmsk(:,:,1) )
+            llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+            iloc(1:3,2) = MAXLOC( ABS(  uu(:,:,:,       Kmm)), mask = llmsk(:,:,:) )
+            llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+            iloc(1:3,3) = MINLOC(       ts(:,:,:,jp_sal,Kmm) , mask = llmsk(:,:,:) )
+            iloc(1:3,4) = MAXLOC(       ts(:,:,:,jp_sal,Kmm) , mask = llmsk(:,:,:) )
+            DO ji = 1, jptst   ! local domain indices ==> global domain indices, excluding halos
+               iloc(1:2,ji) = (/ mig0(iloc(1,ji)), mjg0(iloc(2,ji)) /)
+            END DO
+            iareamin(:) = narea   ;   iareamax(:) = narea   ;   iareasum(:) = 1         ! this is local information
+         ENDIF
+         !
+         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
+         CALL wrt_line( ctmp2, kt, '|ssh| max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) )
+         CALL wrt_line( ctmp3, kt, '|U|   max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) )
+         CALL wrt_line( ctmp4, kt, 'Sal   min', zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )
+         CALL wrt_line( ctmp5, kt, 'Sal   max', zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) )
+         IF( Agrif_Root() ) THEN
+            WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files'
+         ELSE
+            WRITE(ctmp6,*) '      ===> output of last computed fields in '//TRIM(Agrif_CFixed())//'_output.abort* files'
+         ENDIF
+         !
+         CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file
+         !
+         IF( ll_colruns .OR. jpnij == 1 ) THEN   ! all processes synchronized -> use lwp to print in opened ocean.output files
+            IF(lwp) THEN   ;   CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 )
+            ELSE           ;   nstop = MAX(1, nstop)   ! make sure nstop > 0 (automatically done when calling ctl_stop)
+            ENDIF
+         ELSE                                    ! only mpi subdomains with errors are here -> STOP now
+            CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 )
+         ENDIF
+         !
+      ENDIF
+      !
+      IF( nstop > 0 ) THEN                                                  ! an error was detected and we did not abort yet...
+         ngrdstop = Agrif_Fixed()                                           ! store which grid got this error
+         IF( .NOT. ll_colruns .AND. jpnij > 1 )   CALL ctl_stop( 'STOP' )   ! we must abort here to avoid MPI deadlock
+      ENDIF
+      !
+9500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16)
+      !
+   END SUBROUTINE stp_ctl
+
+
+   SUBROUTINE wrt_line( cdline, kt, cdprefix, pval, kloc, ksum, kmin, kmax )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE wrt_line  ***
+      !!
+      !! ** Purpose :   write information line
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*),      INTENT(  out) ::   cdline
+      CHARACTER(len=*),      INTENT(in   ) ::   cdprefix
+      REAL(wp),              INTENT(in   ) ::   pval
+      INTEGER, DIMENSION(3), INTENT(in   ) ::   kloc
+      INTEGER,               INTENT(in   ) ::   kt, ksum, kmin, kmax
+      !
+      CHARACTER(len=80) ::   clsuff
+      CHARACTER(len=9 ) ::   clkt, clsum, clmin, clmax
+      CHARACTER(len=9 ) ::   cli, clj, clk
+      CHARACTER(len=1 ) ::   clfmt
+      CHARACTER(len=4 ) ::   cl4   ! needed to be able to compile with Agrif, I don't know why
+      INTEGER           ::   ifmtk
+      !!----------------------------------------------------------------------
+      WRITE(clkt , '(i9)') kt
+      
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpnij  ,wp))) + 1     ! how many digits to we need to write ? (we decide max = 9)
+      !!! WRITE(clsum, '(i'//clfmt//')') ksum                   ! this is creating a compilation error with AGRIF
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clsum, cl4) ksum
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1    ! how many digits to we need to write ? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clmin, cl4) kmin-1
+                                   WRITE(clmax, cl4) kmax-1
+      !
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpiglo,wp))) + 1      ! how many digits to we need to write jpiglo? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(cli, cl4) kloc(1)      ! this is ok with AGRIF
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpjglo,wp))) + 1      ! how many digits to we need to write jpjglo? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clj, cl4) kloc(2)      ! this is ok with AGRIF
+      !
+      IF( ksum == 1 ) THEN   ;   WRITE(clsuff,9100) TRIM(clmin)
+      ELSE                   ;   WRITE(clsuff,9200) TRIM(clsum), TRIM(clmin), TRIM(clmax)
+      ENDIF
+      IF(kloc(3) == 0) THEN
+         ifmtk = INT(LOG10(REAL(jpk,wp))) + 1                   ! how many digits to we need to write jpk? (we decide max = 9)
+         clk = REPEAT(' ', ifmtk)                               ! create the equivalent in blank string
+         WRITE(cdline,9300) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), clk(1:ifmtk), TRIM(clsuff)
+      ELSE
+         WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpk,wp))) + 1      ! how many digits to we need to write jpk? (we decide max = 9)
+         !!! WRITE(clk, '(i'//clfmt//')') kloc(3)               ! this is creating a compilation error with AGRIF
+         cl4 = '(i'//clfmt//')'   ;   WRITE(clk, cl4) kloc(3)   ! this is ok with AGRIF
+         WRITE(cdline,9400) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj),    TRIM(clk), TRIM(clsuff)
+      ENDIF
+      !
+9100  FORMAT('MPI rank ', a)
+9200  FORMAT('found in ', a, ' MPI tasks, spread out among ranks ', a, ' to ', a)
+9300  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j   ', a, ' ', a, ' ', a, ' ', a)
+9400  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j k ', a, ' ', a, ' ', a, ' ', a)
+      !
+   END SUBROUTINE wrt_line
+
+
+   !!======================================================================
+END MODULE stpctl
Index: NEMO/trunk/tests/DIA_GPU/cpp_DIA_GPU.fcm
===================================================================
--- NEMO/trunk/tests/DIA_GPU/cpp_DIA_GPU.fcm	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/cpp_DIA_GPU.fcm	(revision 14854)
@@ -0,0 +1,1 @@
+bld::tool::fppkeys   key_si3 key_xios key_gpu key_qco
Index: NEMO/trunk/tests/DIA_GPU/internal_nl_workaround.pl
===================================================================
--- NEMO/trunk/tests/DIA_GPU/internal_nl_workaround.pl	(revision 14854)
+++ NEMO/trunk/tests/DIA_GPU/internal_nl_workaround.pl	(revision 14854)
@@ -0,0 +1,23 @@
+INUNITS=( numnam_ref numnam_cfg numnat_ref numnat_cfg numtrc_ref numtrc_cfg numnam_ice_ref numnam_ice_cfg numnamsed_ref numnamsed_cfg numnatp_cfg numnatp_ref )
+#
+# build a list of files that need to be changed
+#
+listfile=tmplistfile$$.txt
+for iunit in ${INUNITS[@]}
+do
+  grep -l $iunit `find ./ -name '*.[fFh]90'` >>  $listfile
+done
+allfiles=`cat $listfile | sort -u`
+echo $allfiles
+#
+if [ -f $listfile ] ; then rm $listfile; fi
+for f in  $allfiles
+do
+ echo "Working on " $f
+ n=0
+ for n in `seq 0 1 $(( ${#INUNITS[*]} - 1 ))`
+ do
+   numnam=${INUNITS[$n]}
+   perl -ni -e 'unless ( m@.*\s*READ\s*\(\s*'${INUNITS[$n]}'\s*,\s*[a-z,0-9]*.*@) { print  } else { $line= $_ ; $line=~s@(.*\s*READ\s*\()(\s*)('${INUNITS[$n]}')(\s*,\s*)([a-z0-9_]*)(.*)@\1\2\3(INDEX(\3,"\5 ")-1:)\4\5\6@i ; print $line }' $f
+ done
+done
Index: NEMO/trunk/tests/README.rst
===================================================================
--- NEMO/trunk/tests/README.rst	(revision 14843)
+++ NEMO/trunk/tests/README.rst	(revision 14854)
@@ -224,4 +224,10 @@
   a coupled configuration through OASIS. See CPL_OASIS/README.md for more information.
 
+DIA_GPU
+---------
+| This is a demonstrator of diagnostic DIAHSB ported to GPU using CUDA Fortran. 
+  Memory communications between host and device are asynchronous given the device has that capability. 
+  This experiment is target for ORCA2_ICE_PISCES
+
 TSUNAMI
 ---------
Index: NEMO/trunk/tests/demo_cfgs.txt
===================================================================
--- NEMO/trunk/tests/demo_cfgs.txt	(revision 14843)
+++ NEMO/trunk/tests/demo_cfgs.txt	(revision 14854)
@@ -12,4 +12,5 @@
 STATION_ASF OCE ICE
 CPL_OASIS  OCE TOP ICE NST
+DIA_GPU OCE ICE
 SWG OCE SWE
 C1D_ASICS OCE
