Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml	(revision 6962)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml	(revision 6963)
@@ -88,5 +88,5 @@
 	  <field field_ref="u_masstr"     name="vozomatr" />
 	  <field field_ref="u_heattr"     name="sozohetr" />
-      <field field_ref="u_salttr"     name="sozosatr" />
+	  <field field_ref="u_salttr"     name="sozosatr" />
 	</file>
 	
@@ -99,5 +99,5 @@
 	  <field field_ref="v_masstr"     name="vomematr" />
 	  <field field_ref="v_heattr"     name="somehetr" />
-      <field field_ref="v_salttr"     name="somesatr" />
+	  <field field_ref="v_salttr"     name="somesatr" />
 	</file>
 	
@@ -116,14 +116,8 @@
           <field field_ref="iceconc"         name="siconc" />
 
-          <field field_ref="vfxbog"          name="vfxbog" />
-          <field field_ref="vfxdyn"          name="vfxdyn" />
-          <field field_ref="vfxopw"          name="vfxopw" />
-          <field field_ref="vfxsni"          name="vfxsni" />
-          <field field_ref="vfxsum"          name="vfxsum" />
-          <field field_ref="vfxbom"          name="vfxbom" />
-          <field field_ref="vfxres"          name="vfxres" />
           <field field_ref="vfxice"          name="vfxice" />
           <field field_ref="vfxsnw"          name="vfxsnw" />
           <field field_ref="vfxsub"          name="vfxsub" />
+          <field field_ref="vfxsub_err"      name="vfxsub_err" />
           <field field_ref="vfxspr"          name="vfxspr" />
 
@@ -134,37 +128,11 @@
           <field field_ref="destrp"          name="destrp" />
 
-          <field field_ref="sfxbri"          name="sfxbri" />
-          <field field_ref="sfxdyn"          name="sfxdyn" />
-          <field field_ref="sfxres"          name="sfxres" />
-          <field field_ref="sfxbog"          name="sfxbog" />
-          <field field_ref="sfxbom"          name="sfxbom" />
-          <field field_ref="sfxsum"          name="sfxsum" />
-          <field field_ref="sfxsni"          name="sfxsni" />
-          <field field_ref="sfxopw"          name="sfxopw" />
           <field field_ref="sfx"             name="sfx"    />
 
-          <field field_ref="hfxsum"          name="hfxsum"     />
-          <field field_ref="hfxbom"          name="hfxbom"     />
-          <field field_ref="hfxbog"          name="hfxbog"     />
-          <field field_ref="hfxdif"          name="hfxdif"     />
-          <field field_ref="hfxopw"          name="hfxopw"     />
           <field field_ref="hfxout"          name="hfxout"     />
           <field field_ref="hfxin"           name="hfxin"      />
-          <field field_ref="hfxsnw"          name="hfxsnw"     />
-          <field field_ref="hfxerr"          name="hfxerr"     />
-          <field field_ref="hfxerr_rem"      name="hfxerr_rem" />
-
-	  <!-- ice-ocean heat flux from mass exchange -->
-          <field field_ref="hfxdyn"          name="hfxdyn" />
-          <field field_ref="hfxres"          name="hfxres" />
-          <field field_ref="hfxthd"          name="hfxthd" />
-	  <!-- ice-atm. heat flux from mass exchange -->
-          <field field_ref="hfxsub"          name="hfxsub" />
-          <field field_ref="hfxspr"          name="hfxspr" />
+
 
 	  <!-- diags -->
-          <field field_ref="hfxdhc"          name="hfxdhc" />
-          <field field_ref="hfxtur"          name="hfxtur" />
-
           <field field_ref="isst"            name="sst"    />
           <field field_ref="isss"            name="sss"    />
@@ -209,4 +177,5 @@
           <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"   /> 
@@ -214,56 +183,26 @@
           <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="ibgvoltot"    name="ibgvoltot"  />
-          <field field_ref="sbgvoltot"    name="sbgvoltot"  />
-          <field field_ref="ibgarea"      name="ibgarea"    />
-          <field field_ref="ibgsaline"    name="ibgsaline"  />
-          <field field_ref="ibgtemper"    name="ibgtemper"  />
-          <field field_ref="ibgheatco"    name="ibgheatco"  />
-          <field field_ref="sbgheatco"    name="sbgheatco"  />
-          <field field_ref="ibgsaltco"    name="ibgsaltco"  />
-
-          <field field_ref="ibgvfx"       name="ibgvfx"     />
-          <field field_ref="ibgvfxbog"    name="ibgvfxbog"  />
-          <field field_ref="ibgvfxopw"    name="ibgvfxopw"  />
-          <field field_ref="ibgvfxsni"    name="ibgvfxsni"  />
-          <field field_ref="ibgvfxdyn"    name="ibgvfxdyn"  />
-          <field field_ref="ibgvfxbom"    name="ibgvfxbom"  />
-          <field field_ref="ibgvfxsum"    name="ibgvfxsum"  />
-          <field field_ref="ibgvfxres"    name="ibgvfxres"  />
-          <field field_ref="ibgvfxspr"    name="ibgvfxspr"  />
-          <field field_ref="ibgvfxsnw"    name="ibgvfxsnw"  />
-          <field field_ref="ibgvfxsub"    name="ibgvfxsub"  />
-
-          <field field_ref="ibgsfx"       name="ibgsfx"     />
-          <field field_ref="ibgsfxbri"    name="ibgsfxbri"  />
-          <field field_ref="ibgsfxdyn"    name="ibgsfxdyn"  />
-          <field field_ref="ibgsfxres"    name="ibgsfxres"  />
-          <field field_ref="ibgsfxbog"    name="ibgsfxbog"  />
-          <field field_ref="ibgsfxopw"    name="ibgsfxopw"  />
-          <field field_ref="ibgsfxsni"    name="ibgsfxsni"  />
-          <field field_ref="ibgsfxbom"    name="ibgsfxbom"  />
-          <field field_ref="ibgsfxsum"    name="ibgsfxsum"  />
-
-          <field field_ref="ibghfxdhc"    name="ibghfxdhc"  />
-          <field field_ref="ibghfxspr"    name="ibghfxspr"  />
-
-          <field field_ref="ibghfxres"    name="ibghfxres"  />
-          <field field_ref="ibghfxsub"    name="ibghfxsub"  />
-          <field field_ref="ibghfxdyn"    name="ibghfxdyn"  />
-          <field field_ref="ibghfxthd"    name="ibghfxthd"  />
-          <field field_ref="ibghfxsum"    name="ibghfxsum"  />
-          <field field_ref="ibghfxbom"    name="ibghfxbom"  />
-          <field field_ref="ibghfxbog"    name="ibghfxbog"  />
-          <field field_ref="ibghfxdif"    name="ibghfxdif"  />
-          <field field_ref="ibghfxopw"    name="ibghfxopw"  />
-          <field field_ref="ibghfxout"    name="ibghfxout"  />
-          <field field_ref="ibghfxin"     name="ibghfxin"   />
-          <field field_ref="ibghfxsnw"    name="ibghfxsnw"  />
-
-          <field field_ref="ibgfrcvol"    name="ibgfrcvol"  />
-          <field field_ref="ibgfrcsfx"    name="ibgfrcsfx"  />
-          <field field_ref="ibgvolgrm"    name="ibgvolgrm"  />
+	  <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" />
+	  
+	  <field field_ref="ibgvolume"    name="ibgvolume"   />
+	  <field field_ref="ibgsaltco"    name="ibgsaltco"   />
+	  <field field_ref="ibgheatco"    name="ibgheatco"   />
+          <field field_ref="ibgheatfx"    name="ibgheatfx"   />
+	  
+	  <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>
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6962)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6963)
@@ -225,5 +225,7 @@
          <field id="hflx_rain_cea" long_name="heat flux due to rainfall"                                standard_name="temperature_flux_due_to_rainfall_expressed_as_heat_flux_into_sea_water"        unit="W/m2"     />
          <field id="hflx_evap_cea" long_name="heat flux due to evaporation"                             standard_name="temperature_flux_due_to_evaporation_expressed_as_heat_flux_out_of_sea_water"   unit="W/m2"     />
-         <field id="hflx_snow_cea" long_name="heat flux due to snow falling over ice-free ocean"        standard_name="heat_flux_into_sea_water_due_to_snow_thermodynamics"                           unit="W/m2"     />
+         <field id="hflx_snow_cea" long_name="heat flux due to snow falling"                            standard_name="heat_flux_onto_ocean_and_ice_due_to_snow_thermodynamics"                       unit="W/m2"     />
+         <field id="hflx_snow_ai_cea" long_name="heat flux due to snow falling over ice"                standard_name="heat_flux_onto_ice_due_to_snow_thermodynamics"                                 unit="W/m2"     />
+         <field id="hflx_snow_ao_cea" long_name="heat flux due to snow falling over ice-free ocean"     standard_name="heat_flux_onto_sea_water_due_to_snow_thermodynamics"                           unit="W/m2"     />
          <field id="hflx_ice_cea"  long_name="heat flux due to ice thermodynamics"                      standard_name="heat_flux_into_sea_water_due_to_sea_ice_thermodynamics"                        unit="W/m2"     />
          <field id="hflx_rnf_cea"  long_name="heat flux due to runoffs"                                 standard_name="temperature_flux_due_to_runoff_expressed_as_heat_flux_into_sea_water"          unit="W/m2"     />
@@ -331,4 +333,5 @@
          <field id="vfxsnw"       long_name="snw melt/growth"                                              unit="m/day"   />
          <field id="vfxsub"       long_name="snw sublimation"                                              unit="m/day"   />
+         <field id="vfxsub_err"   long_name="excess of snw sublimation sent to ocean"                      unit="m/day"   />
          <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="m/day"   />
          <field id="vfxthin"      long_name="daily thermo ice prod. for thin ice(<20cm) + open water"      unit="m/day"   />
@@ -499,16 +502,15 @@
     	 <field id="bgtemper"     long_name="drift in global mean temperature wrt timestep 1"                 standard_name="change_over_time_in_sea_water_potential_temperature"   unit="degC"     />
     	 <field id="bgsaline"     long_name="drift in global mean salinity wrt timestep 1"                    standard_name="change_over_time_in_sea_water_practical_salinity"      unit="1e-3"     />
-    	 <field id="bgheatco"     long_name="drift in global mean heat content wrt timestep 1"                                                                                      unit="10^9J"    />
-    	 <field id="bgsaltco"     long_name="drift in global mean salt content wrt timestep 1"                                                                                      unit="1e-3*m3"  />
+    	 <field id="bgheatco"     long_name="drift in global mean heat content wrt timestep 1"                                                                                      unit="1.e20J"   />
+     	 <field id="bgheatfx"     long_name="drift in global mean heat flux    wrt timestep 1"                                                                                      unit="W/m2"     />
+   	 <field id="bgsaltco"     long_name="drift in global mean salt content wrt timestep 1"                                                                                      unit="1e-3*km3" />
     	 <field id="bgvolssh"     long_name="drift in global mean ssh volume wrt timestep 1"                                                                                        unit="km3"      />
          <field id="bgvole3t"     long_name="drift in global mean volume variation (e3t) wrt timestep 1"                                                                            unit="km3"      />
-    	 <field id="bgvoltot"     long_name="drift in global mean volume wrt timestep 1"                                                                                            unit="km3"      />
-         <!-- NOTE: No matching iom_put call -->
-    	 <field id="bgsshtot"     long_name="drift in global mean ssh wrt timestep 1"                         standard_name="global_average_sea_level_change"                       unit="m"        />
-    	 <field id="bgfrcvol"     long_name="drift in global mean volume from forcing wrt timestep 1"                                                                               unit="km3"      />
-    	 <field id="bgfrctem"     long_name="drift in global mean heat content from forcing wrt timestep 1"                                                                         unit="10^9J"    />
-    	 <field id="bgfrcsal"     long_name="drift in global mean salt content from forcing wrt timestep 1"                                                                         unit="1e-3*km3" />
-    	 <field id="bgmistem"     long_name="global mean temperature error due to free surface"                                                                                     unit="degC"     />
-    	 <field id="bgmissal"     long_name="global mean salinity error due to free surface"                                                                                        unit="1e-3"     />
+    	 <field id="bgfrcvol"     long_name="global mean volume from forcing"                                                                                                       unit="km3"      />
+    	 <field id="bgfrctem"     long_name="global mean heat content from forcing"                                                                                                 unit="1.e20J"   />
+         <field id="bgfrchfx"     long_name="global mean heat flux from forcing"                                                                                                    unit="W/m2"     />
+    	 <field id="bgfrcsal"     long_name="global mean salt content from forcing"                                                                                                 unit="1e-3*km3" />
+    	 <field id="bgmistem"     long_name="global mean temperature error due to free surface (no vvl)"                                                                            unit="degC"     />
+    	 <field id="bgmissal"     long_name="global mean salinity error due to free surface (no vvl)"                                                                               unit="1e-3"     />
       </field_group>
 
@@ -517,55 +519,23 @@
       <field_group id="SBC_scalar"  domain_ref="1point" >
          <!-- available with ln_limdiaout -->
-         <field id="ibgvoltot"    long_name="global mean ice volume"                                 unit="km3"        />
-         <field id="sbgvoltot"    long_name="global mean snow volume"                                unit="km3"        />
-         <field id="ibgarea"      long_name="global mean ice area"                                   unit="km2"        />
-         <field id="ibgsaline"    long_name="global mean ice salinity"                               unit="1e-3"       />
-         <field id="ibgtemper"    long_name="global mean ice temperature"                            unit="degC"       />
-         <field id="ibgheatco"    long_name="global mean ice heat content"                           unit="10^20J"     />
-         <field id="sbgheatco"    long_name="global mean snow heat content"                          unit="10^20J"     />
-         <field id="ibgsaltco"    long_name="global mean ice salt content"                           unit="1e-3*km3"   />
-
-         <field id="ibgvfx"       long_name="global mean volume flux (emp)"                          unit="m/day"      />
-         <field id="ibgvfxbog"    long_name="global mean volume flux (bottom growth)"                unit="m/day"      />
-         <field id="ibgvfxopw"    long_name="global mean volume flux (open water growth)"            unit="m/day"      />
-         <field id="ibgvfxsni"    long_name="global mean volume flux (snow-ice growth)"              unit="m/day"      />
-         <field id="ibgvfxdyn"    long_name="global mean volume flux (dynamic growth)"               unit="m/day"      />
-         <field id="ibgvfxbom"    long_name="global mean volume flux (bottom melt)"                  unit="m/day"      />
-         <field id="ibgvfxsum"    long_name="global mean volume flux (surface melt)"                 unit="m/day"      />
-         <field id="ibgvfxres"    long_name="global mean volume flux (resultant)"                    unit="m/day"      />
-         <field id="ibgvfxspr"    long_name="global mean volume flux (snow precip)"                  unit="m/day"      />
-         <field id="ibgvfxsnw"    long_name="global mean volume flux (snow melt)"                    unit="m/day"      />
-         <field id="ibgvfxsub"    long_name="global mean volume flux (snow sublimation)"             unit="m/day"      />
-
-         <field id="ibgsfx"       long_name="global mean salt flux (total)"                          unit="1e-3*m/day" />
-         <field id="ibgsfxbri"    long_name="global mean salt flux (brines)"                         unit="1e-3*m/day" />
-         <field id="ibgsfxdyn"    long_name="global mean salt flux (dynamic)"                        unit="1e-3*m/day" />
-         <field id="ibgsfxres"    long_name="global mean salt flux (resultant)"                      unit="1e-3*m/day" />
-         <field id="ibgsfxbog"    long_name="global mean salt flux (thermo)"                         unit="1e-3*m/day" />
-         <field id="ibgsfxopw"    long_name="global mean salt flux (thermo)"                         unit="1e-3*m/day" />
-         <field id="ibgsfxsni"    long_name="global mean salt flux (thermo)"                         unit="1e-3*m/day" />
-         <field id="ibgsfxbom"    long_name="global mean salt flux (thermo)"                         unit="1e-3*m/day" />
-         <field id="ibgsfxsum"    long_name="global mean salt flux (thermo)"                         unit="1e-3*m/day" />
-         <field id="ibgsfxsub"    long_name="global mean salt flux (thermo)"                         unit="1e-3*m/day" />
-
-         <field id="ibghfxdhc"    long_name="Heat content variation in snow and ice"                 unit="W"          />
-         <field id="ibghfxspr"    long_name="Heat content of snow precip"                            unit="W"          />
-
-         <field id="ibghfxthd"    long_name="heat fluxes from ice-ocean exchange during thermo"      unit="W"          />
-         <field id="ibghfxsum"    long_name="heat fluxes causing surface ice melt"                   unit="W"          />
-         <field id="ibghfxbom"    long_name="heat fluxes causing bottom ice melt"                    unit="W"          />
-         <field id="ibghfxbog"    long_name="heat fluxes causing bottom ice growth"                  unit="W"          />
-         <field id="ibghfxdif"    long_name="heat fluxes causing ice temperature change"             unit="W"          />
-         <field id="ibghfxopw"    long_name="heat fluxes causing open water ice formation"           unit="W"          />
-         <field id="ibghfxdyn"    long_name="heat fluxes from ice-ocean exchange during dynamic"     unit="W"          />
-         <field id="ibghfxres"    long_name="heat fluxes from ice-ocean exchange during resultant"   unit="W"          />
-         <field id="ibghfxsub"    long_name="heat fluxes from sublimation"                           unit="W"          />
-         <field id="ibghfxsnw"    long_name="heat fluxes from snow-ocean exchange"                   unit="W"          />
-         <field id="ibghfxout"    long_name="non solar heat fluxes received by the ocean"            unit="W"          />
-         <field id="ibghfxin"     long_name="total heat fluxes at the ice surface"                   unit="W"          />
-
-         <field id="ibgfrcvol"    long_name="global mean forcing volume (emp)"                       unit="km3"        />
-         <field id="ibgfrcsfx"    long_name="global mean forcing salt   (sfx)"                       unit="1e-3*km3"   />
-         <field id="ibgvolgrm"    long_name="global mean ice growth+melt volume"                     unit="km3"        />
+         <field id="ibgfrcvoltop"    long_name="global mean ice/snow forcing at interface ice/snow-atm (volume equivalent ocean volume)"   unit="km3"       />
+         <field id="ibgfrcvolbot"    long_name="global mean ice/snow forcing at interface ice/snow-ocean (volume equivalent ocean volume)" unit="km3"       />
+         <field id="ibgfrctemtop"    long_name="global mean heat on top of ice/snw/ocean-atm "                                             unit="1e20J"     />
+         <field id="ibgfrctembot"    long_name="global mean heat below ice (on top of ocean) "                                             unit="1e20J"     />
+         <field id="ibgfrcsal"       long_name="global mean ice/snow forcing (salt equivalent ocean volume)"                               unit="pss*km3"   />
+         <field id="ibgfrchfxtop"    long_name="global mean heat flux on top of ice/snw/ocean-atm "                                        unit="W/m2"      />
+         <field id="ibgfrchfxbot"    long_name="global mean heat flux below ice (on top of ocean) "                                        unit="W/m2"      />
+ 
+         <field id="ibgvolume"    long_name="drift in ice/snow volume (equivalent ocean volume)"            unit="km3"        />
+         <field id="ibgsaltco"    long_name="drift in ice salt content (equivalent ocean volume)"           unit="pss*km3"    />
+         <field id="ibgheatco"    long_name="drift in ice/snow heat content"                                unit="1e20J"      />
+         <field id="ibgheatfx"    long_name="drift in ice/snow heat flux"                                   unit="W/m2"       />
+
+         <field id="ibgvol_tot"    long_name="global mean ice volume"                                       unit="km3"        />
+         <field id="sbgvol_tot"    long_name="global mean snow volume"                                      unit="km3"        />
+         <field id="ibgarea_tot"   long_name="global mean ice area"                                         unit="km2"        />
+         <field id="ibgsalt_tot"   long_name="global mean ice salt content"                                 unit="1e-3*km3"   />
+         <field id="ibgheat_tot"   long_name="global mean ice heat content"                                 unit="1e20J"      />
+         <field id="sbgheat_tot"   long_name="global mean snow heat content"                                unit="1e20J"      />
       </field_group>
   
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6962)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6963)
@@ -243,5 +243,4 @@
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sist        !: Average Sea-Ice Surface Temperature [Kelvin]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   icethi      !: total ice thickness (for all categories) (diag only)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_bo        !: Sea-Ice bottom temperature [Kelvin]     
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frld        !: Leads fraction = 1 - ice fraction
@@ -320,18 +319,20 @@
    !                                                                  !  this is an extensive variable that has to be transported
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i     !: Sea-Ice Age (days)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ov_i    !: Sea-Ice Age times volume per area (days.m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i    !: Sea-Ice Age times ice area (days)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bv_i    !: brine volume
 
    !! Variables summed over all categories, or associated to all the ice in a single grid cell
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_ice, v_ice   !: components of the ice velocity (m/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tio_u, tio_v   !: components of the ice-ocean stress (N/m2)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vt_i , vt_s    !: ice and snow total volume per unit area (m)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   at_i           !: ice total fractional area (ice concentration)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ato_i          !: =1-at_i ; total open water fractional area
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   et_i , et_s    !: ice and snow total heat content
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ot_i           !: mean age over all categories
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_i           !: mean ice temperature over all categories
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bv_i           !: brine volume averaged over all categories
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   smt_i          !: mean sea ice salinity averaged over all categories [PSU]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_i         !: mean ice temperature over all categories
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bvm_i        !: brine volume averaged over all categories
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   smt_i        !: mean sea ice salinity averaged over all categories [PSU]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_su        !: mean surface temperature over all categories
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htm_i        !: mean ice  thickness over all categories
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htm_s        !: mean snow thickness over all categories
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   om_i         !: mean ice age over all categories
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_s        !: Snow temperatures [K]
@@ -435,5 +436,5 @@
 
       ii = ii + 1
-      ALLOCATE( sist   (jpi,jpj) , icethi (jpi,jpj) , t_bo   (jpi,jpj) ,                        &
+      ALLOCATE( sist   (jpi,jpj) , t_bo   (jpi,jpj) ,                        &
          &      frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,                        &
          &      wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) ,                        &
@@ -456,10 +457,11 @@
          &      v_s  (jpi,jpj,jpl) , ht_s (jpi,jpj,jpl) , t_su (jpi,jpj,jpl) ,     &
          &      sm_i (jpi,jpj,jpl) , smv_i(jpi,jpj,jpl) , o_i  (jpi,jpj,jpl) ,     &
-         &      ov_i (jpi,jpj,jpl) , oa_i (jpi,jpj,jpl)                      , STAT=ierr(ii) )
-      ii = ii + 1
-      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , tio_u(jpi,jpj) , tio_v(jpi,jpj) ,     &
+         &      oa_i (jpi,jpj,jpl) , bv_i (jpi,jpj,jpl) , STAT=ierr(ii) )
+      ii = ii + 1
+      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,      &
          &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i (jpi,jpj) , ato_i(jpi,jpj) ,     &
-         &      et_i (jpi,jpj) , et_s (jpi,jpj) , ot_i (jpi,jpj) , tm_i (jpi,jpj) ,     &
-         &      bv_i (jpi,jpj) , smt_i(jpi,jpj)                                   , STAT=ierr(ii) )
+         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i (jpi,jpj) , bvm_i(jpi,jpj) ,     &
+         &      smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) ,     &
+         &      om_i (jpi,jpj)                              , STAT=ierr(ii) )
       ii = ii + 1
       ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) )
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 6962)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 6963)
@@ -288,6 +288,7 @@
 #if ! defined key_bdy
       ! heat flux
-      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - SUM( qevap_ice * a_i_b, dim=3 ) )  &
-         &              * e12t * tmask(:,:,1) * zconv ) 
+      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es  &
+      !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check)
+         &              ) * e12t * tmask(:,:,1) * zconv ) 
       ! salt flux
       zsfx  = glob_sum( ( sfx + diag_smvi ) * e12t * tmask(:,:,1) * zconv ) * rday
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 6962)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 6963)
@@ -31,8 +31,9 @@
 
    PUBLIC   lim_diahsb        ! routine called by ice_step.F90
-
-   real(wp) ::   frc_sal, frc_vol   ! global forcing trends
-   real(wp) ::   bg_grme            ! global ice growth+melt trends
-
+   PUBLIC   lim_diahsb_init   ! routine called in sbcice_lim.F90
+
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   vol_loc_ini, sal_loc_ini, tem_loc_ini ! initial volume, salt and heat contents
+   REAL(wp)                              ::   frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot  ! global forcing trends
+   
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -46,5 +47,5 @@
 CONTAINS
 
-   SUBROUTINE lim_diahsb
+   SUBROUTINE lim_diahsb( kt )
       !!---------------------------------------------------------------------------
       !!                  ***  ROUTINE lim_diahsb  ***
@@ -53,169 +54,86 @@
       !!	
       !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt    ! number of iteration
       !!
-      real(wp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc
-      real(wp)   ::   zbg_sfx, zbg_sfx_bri, zbg_sfx_bog, zbg_sfx_bom, zbg_sfx_sum, zbg_sfx_sni,   &
-      &               zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn, zbg_sfx_sub 
-      real(wp)   ::   zbg_vfx, zbg_vfx_bog, zbg_vfx_opw, zbg_vfx_sni, zbg_vfx_dyn
-      real(wp)   ::   zbg_vfx_bom, zbg_vfx_sum, zbg_vfx_res, zbg_vfx_spr, zbg_vfx_snw, zbg_vfx_sub  
-      real(wp)   ::   zbg_hfx_dhc, zbg_hfx_spr
-      real(wp)   ::   zbg_hfx_res, zbg_hfx_sub, zbg_hfx_dyn, zbg_hfx_thd, zbg_hfx_snw, zbg_hfx_out, zbg_hfx_in   
-      real(wp)   ::   zbg_hfx_sum, zbg_hfx_bom, zbg_hfx_bog, zbg_hfx_dif, zbg_hfx_opw
-      real(wp)   ::   z_frc_vol, z_frc_sal, z_bg_grme 
-      real(wp)   ::   z1_area                     !    -     -
-      REAL(wp)   ::   ztmp
+      real(wp)   ::   zbg_ivol, zbg_svol, zbg_area, zbg_isal, zbg_item ,zbg_stem
+      REAL(wp)   ::   z_frc_voltop, z_frc_volbot, z_frc_sal, z_frc_temtop, z_frc_tembot  
+      REAL(wp)   ::   zdiff_vol, zdiff_sal, zdiff_tem  
       !!---------------------------------------------------------------------------
       IF( nn_timing == 1 )   CALL timing_start('lim_diahsb')
 
-      IF( numit == nstart ) CALL lim_diahsb_init 
-
-      ! 1/area
-      z1_area = 1._wp / MAX( glob_sum( e12t(:,:) * tmask(:,:,1) ), epsi06 )
-
-      rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( e12t(:,:) * tmask(:,:,1) ) - epsi06 ) )
-      ! ----------------------- !
-      ! 1 -  Content variations !
-      ! ----------------------- !
-      zbg_ivo = glob_sum( vt_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume ice 
-      zbg_svo = glob_sum( vt_s(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume snow
-      zbg_are = glob_sum( at_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! area
-      zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) )       ! mean salt content
-      zbg_tem = glob_sum( ( tm_i(:,:) - rt0 ) * vt_i(:,:) * e12t(:,:) * tmask(:,:,1) )  ! mean temp content
-
-      !zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content
-      !zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_svo,epsi06 ) ! snow heat content
-
-      ! Volume
-      ztmp = rswitch * z1_area * r1_rau0 * rday
-      zbg_vfx     = ztmp * glob_sum(     emp(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) )
-
-      ! Salt
-      zbg_sfx     = ztmp * glob_sum(     sfx(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) )
-
-      zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) )
-      zbg_sfx_sub = ztmp * glob_sum( sfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) )
-
-      ! Heat budget
-      zbg_ihc      = glob_sum( et_i(:,:) * e12t(:,:) * 1.e-20 ) ! ice heat content  [1.e20 J]
-      zbg_shc      = glob_sum( et_s(:,:) * e12t(:,:) * 1.e-20 ) ! snow heat content [1.e20 J]
-      zbg_hfx_dhc  = glob_sum( diag_heat(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_spr  = glob_sum( hfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-
-      zbg_hfx_thd  = glob_sum( hfx_thd(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_dyn  = glob_sum( hfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_res  = glob_sum( hfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_sub  = glob_sum( hfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_snw  = glob_sum( hfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_sum  = glob_sum( hfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_bom  = glob_sum( hfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_bog  = glob_sum( hfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_dif  = glob_sum( hfx_dif(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_opw  = glob_sum( hfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_out  = glob_sum( hfx_out(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-      zbg_hfx_in   = glob_sum(  hfx_in(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
-    
-      ! --------------------------------------------- !
-      ! 2 - Trends due to forcing and ice growth/melt !
-      ! --------------------------------------------- !
-      z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes
-      z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) ! salt fluxes
-      z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + &
-                          &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + &
-                          &     wfx_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes
-      !
-      frc_vol  = frc_vol  + z_frc_vol  * rdt_ice
-      frc_sal  = frc_sal  + z_frc_sal  * rdt_ice
-      bg_grme  = bg_grme  + z_bg_grme  * rdt_ice
+      ! ----------------------- !
+      ! 1 -  Contents !
+      ! ----------------------- !
+      zbg_ivol = glob_sum( vt_i(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! ice volume (km3)
+      zbg_svol = glob_sum( vt_s(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! snow volume (km3)
+      zbg_area = glob_sum( at_i(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-6 )                  ! area (km2)
+      zbg_isal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) ! salt content (pss*km3)
+      zbg_item = glob_sum( et_i * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J)
+      zbg_stem = glob_sum( et_s * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J)
       
-      ! difference
-      !frc_vol = zbg_ivo - frc_vol
-      !frc_sal = zbg_sal - frc_sal
-      
-      ! ----------------------- !
-      ! 3 - Diagnostics writing !
-      ! ----------------------- !
-      rswitch = MAX( 0._wp , SIGN( 1._wp , zbg_ivo - epsi06 ) )
-      !
-      IF( iom_use('ibgvoltot') )   &
-      CALL iom_put( 'ibgvoltot' , zbg_ivo * rhoic * r1_rau0 * 1.e-9        )   ! ice volume (km3 equivalent liquid)         
-      IF( iom_use('sbgvoltot') )   &
-      CALL iom_put( 'sbgvoltot' , zbg_svo * rhosn * r1_rau0 * 1.e-9        )   ! snw volume (km3 equivalent liquid)       
-      IF( iom_use('ibgarea') )   &
-      CALL iom_put( 'ibgarea'   , zbg_are * 1.e-6                          )   ! ice area   (km2)
-      IF( iom_use('ibgsaline') )   &
-      CALL iom_put( 'ibgsaline' , rswitch * zbg_sal / MAX( zbg_ivo, epsi06 ) )   ! ice saline (psu)
-      IF( iom_use('ibgtemper') )   &
-      CALL iom_put( 'ibgtemper' , rswitch * zbg_tem / MAX( zbg_ivo, epsi06 ) )   ! ice temper (C)
-      CALL iom_put( 'ibgheatco' , zbg_ihc                                  )   ! ice heat content (1.e20 J)        
-      CALL iom_put( 'sbgheatco' , zbg_shc                                  )   ! snw heat content (1.e20 J)
-      IF( iom_use('ibgsaltco') )   &
-      CALL iom_put( 'ibgsaltco' , zbg_sal * rhoic * r1_rau0 * 1.e-9        )   ! ice salt content (psu*km3 equivalent liquid)        
-
-      CALL iom_put( 'ibgvfx'    , zbg_vfx                                  )   ! volume flux emp (m/day liquid)
-      CALL iom_put( 'ibgvfxbog' , zbg_vfx_bog                              )   ! volume flux bottom growth     -(m/day equivalent liquid)
-      CALL iom_put( 'ibgvfxopw' , zbg_vfx_opw                              )   ! volume flux open water growth -
-      CALL iom_put( 'ibgvfxsni' , zbg_vfx_sni                              )   ! volume flux snow ice growth   -
-      CALL iom_put( 'ibgvfxdyn' , zbg_vfx_dyn                              )   ! volume flux dynamic growth    -
-      CALL iom_put( 'ibgvfxbom' , zbg_vfx_bom                              )   ! volume flux bottom melt       -
-      CALL iom_put( 'ibgvfxsum' , zbg_vfx_sum                              )   ! volume flux surface melt      -
-      CALL iom_put( 'ibgvfxres' , zbg_vfx_res                              )   ! volume flux resultant         -
-      CALL iom_put( 'ibgvfxspr' , zbg_vfx_spr                              )   ! volume flux from snow precip         -
-      CALL iom_put( 'ibgvfxsnw' , zbg_vfx_snw                              )   ! volume flux from snow melt         -
-      CALL iom_put( 'ibgvfxsub' , zbg_vfx_sub                              )   ! volume flux from sublimation         -
-          
-      CALL iom_put( 'ibgsfx'    , zbg_sfx                                  )   ! salt flux         -(psu*m/day equivalent liquid)       
-      CALL iom_put( 'ibgsfxbri' , zbg_sfx_bri                              )   ! salt flux brines  -      
-      CALL iom_put( 'ibgsfxdyn' , zbg_sfx_dyn                              )   ! salt flux dynamic -    
-      CALL iom_put( 'ibgsfxres' , zbg_sfx_res                              )   ! salt flux result  -    
-      CALL iom_put( 'ibgsfxbog' , zbg_sfx_bog                              )   ! salt flux bottom growth   
-      CALL iom_put( 'ibgsfxopw' , zbg_sfx_opw                              )   ! salt flux open water growth -
-      CALL iom_put( 'ibgsfxsni' , zbg_sfx_sni                              )   ! salt flux snow ice growth   -
-      CALL iom_put( 'ibgsfxbom' , zbg_sfx_bom                              )   ! salt flux bottom melt       -
-      CALL iom_put( 'ibgsfxsum' , zbg_sfx_sum                              )   ! salt flux surface melt      -
-      CALL iom_put( 'ibgsfxsub' , zbg_sfx_sub                              )   ! salt flux sublimation      -
-
-      CALL iom_put( 'ibghfxdhc' , zbg_hfx_dhc                              )   ! Heat content variation in snow and ice [W]
-      CALL iom_put( 'ibghfxspr' , zbg_hfx_spr                              )   ! Heat content of snow precip [W]
-
-      CALL iom_put( 'ibghfxres' , zbg_hfx_res                              )   ! 
-      CALL iom_put( 'ibghfxsub' , zbg_hfx_sub                              )   ! 
-      CALL iom_put( 'ibghfxdyn' , zbg_hfx_dyn                              )   ! 
-      CALL iom_put( 'ibghfxthd' , zbg_hfx_thd                              )   ! 
-      CALL iom_put( 'ibghfxsnw' , zbg_hfx_snw                              )   ! 
-      CALL iom_put( 'ibghfxsum' , zbg_hfx_sum                              )   ! 
-      CALL iom_put( 'ibghfxbom' , zbg_hfx_bom                              )   ! 
-      CALL iom_put( 'ibghfxbog' , zbg_hfx_bog                              )   ! 
-      CALL iom_put( 'ibghfxdif' , zbg_hfx_dif                              )   ! 
-      CALL iom_put( 'ibghfxopw' , zbg_hfx_opw                              )   ! 
-      CALL iom_put( 'ibghfxout' , zbg_hfx_out                              )   ! 
-      CALL iom_put( 'ibghfxin'  , zbg_hfx_in                               )   ! 
-
-      CALL iom_put( 'ibgfrcvol' , frc_vol * 1.e-9                          )   ! vol - forcing     (km3 equivalent liquid) 
-      CALL iom_put( 'ibgfrcsfx' , frc_sal * 1.e-9                          )   ! sal - forcing     (psu*km3 equivalent liquid)   
-      IF( iom_use('ibgvolgrm') )   &
-      CALL iom_put( 'ibgvolgrm' , bg_grme * r1_rau0 * 1.e-9                )   ! vol growth + melt (km3 equivalent liquid)         
-
+      ! ---------------------------!
+      ! 2 - Trends due to forcing  !
+      ! ---------------------------!
+      z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )  ! freshwater flux ice/snow-ocean 
+      z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                     ! freshwater flux ice/snow-atm
+      z_frc_sal    = r1_rau0 * glob_sum( - sfx(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                                            ! salt fluxes ice/snow-ocean
+      z_frc_tembot =           glob_sum( hfx_out(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ocean (and below ice)
+      z_frc_temtop =           glob_sum( hfx_in (:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ice-coean
+      !
+      frc_voltop  = frc_voltop  + z_frc_voltop  * rdt_ice ! km3
+      frc_volbot  = frc_volbot  + z_frc_volbot  * rdt_ice ! km3
+      frc_sal     = frc_sal     + z_frc_sal     * rdt_ice ! km3*pss
+      frc_temtop  = frc_temtop  + z_frc_temtop  * rdt_ice ! 1.e20 J
+      frc_tembot  = frc_tembot  + z_frc_tembot  * rdt_ice ! 1.e20 J
+            
+      ! ----------------------- !
+      ! 3 -  Content variations !
+      ! ----------------------- !
+      zdiff_vol = r1_rau0 * glob_sum( ( rhoic * vt_i(:,:) + rhosn * vt_s(:,:) - vol_loc_ini(:,:)  &  ! freshwater trend (km3) 
+         &                            ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) 
+      zdiff_sal = r1_rau0 * glob_sum( ( rhoic * SUM( smv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:)     &  ! salt content trend (km3*pss)
+         &                            ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )
+      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:)                  &  ! heat content trend (1.e20 J)
+      !  &                            + SUM( qevap_ice * a_i_b, dim=3 ) &     !! clem: I think this line should be commented (but needs a check)
+         &                            ) * e12t(:,:) * tmask(:,:,1) * 1.e-20 )
+
+      ! ----------------------- !
+      ! 4 -  Drifts             !
+      ! ----------------------- !
+      zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot )
+      zdiff_sal = zdiff_sal - frc_sal
+      zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop )
+
+      ! ----------------------- !
+      ! 5 - Diagnostics writing !
+      ! ----------------------- !
+      !
+      IF( iom_use('ibgvolume') )  CALL iom_put( 'ibgvolume' , zdiff_vol        )   ! ice/snow volume  drift            (km3 equivalent ocean water)         
+      IF( iom_use('ibgsaltco') )  CALL iom_put( 'ibgsaltco' , zdiff_sal        )   ! ice salt content drift            (psu*km3 equivalent ocean water)
+      IF( iom_use('ibgheatco') )  CALL iom_put( 'ibgheatco' , zdiff_tem        )   ! ice/snow heat content drift       (1.e20 J)
+      IF( iom_use('ibgheatfx') )  CALL iom_put( 'ibgheatfx' , zdiff_tem /      &   ! ice/snow heat flux drift          (W/m2)
+         &                                                    glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) )
+
+      IF( iom_use('ibgfrcvoltop') )  CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water) 
+      IF( iom_use('ibgfrcvolbot') )  CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water) 
+      IF( iom_use('ibgfrcsal') )     CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)   
+      IF( iom_use('ibgfrctemtop') )  CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)   
+      IF( iom_use('ibgfrctembot') )  CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)   
+      IF( iom_use('ibgfrchfxtop') )  CALL iom_put( 'ibgfrchfxtop' , frc_temtop / & ! heat on top of ice/snw/ocean      (W/m2) 
+         &                                                    glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) )
+      IF( iom_use('ibgfrchfxbot') )  CALL iom_put( 'ibgfrchfxbot' , frc_tembot / & ! heat on top of ocean(below ice)   (W/m2) 
+         &                                                    glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) )
+
+      IF( iom_use('ibgvol_tot' ) )  CALL iom_put( 'ibgvol_tot'  , zbg_ivol     )   ! ice volume                        (km3)
+      IF( iom_use('sbgvol_tot' ) )  CALL iom_put( 'sbgvol_tot'  , zbg_svol     )   ! snow volume                       (km3)
+      IF( iom_use('ibgarea_tot') )  CALL iom_put( 'ibgarea_tot' , zbg_area     )   ! ice area                          (km2)
+      IF( iom_use('ibgsalt_tot') )  CALL iom_put( 'ibgsalt_tot' , zbg_isal     )   ! ice salinity content              (pss*km3)
+      IF( iom_use('ibgheat_tot') )  CALL iom_put( 'ibgheat_tot' , zbg_item     )   ! ice heat content                  (1.e20 J)
+      IF( iom_use('sbgheat_tot') )  CALL iom_put( 'sbgheat_tot' , zbg_stem     )   ! snow heat content                 (1.e20 J)
       !
       IF( lrst_ice )   CALL lim_diahsb_rst( numit, 'WRITE' )
       !
       IF( nn_timing == 1 )   CALL timing_stop('lim_diahsb')
-!
+      !
    END SUBROUTINE lim_diahsb
 
@@ -233,5 +151,4 @@
       !!             - Compute coefficients for conversion
       !!---------------------------------------------------------------------------
-      INTEGER            ::   jk       ! dummy loop indice
       INTEGER            ::   ierror   ! local integer
       !!
@@ -247,5 +164,11 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
       ENDIF
-      !
+      !      
+      ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ierror )
+      IF( ierror > 0 )  THEN
+         CALL ctl_stop( 'lim_diahsb: unable to allocate vol_loc_ini' )
+         RETURN
+      ENDIF
+
       CALL lim_diahsb_rst( nstart, 'READ' )  !* read or initialize all required files
       !
@@ -263,24 +186,35 @@
      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
      !
-     INTEGER ::   id1, id2, id3   ! local integers
      !!----------------------------------------------------------------------
      !
      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
         IF( ln_rstart ) THEN                   !* Read the restart file
-           !id1 = iom_varid( numrir, 'frc_vol'  , ldstop = .TRUE. )
            !
            IF(lwp) WRITE(numout,*) '~~~~~~~'
-           IF(lwp) WRITE(numout,*) ' lim_diahsb_rst at it= ', kt,' date= ', ndastp
-           IF(lwp) WRITE(numout,*) '~~~~~~~'
-           CALL iom_get( numrir, 'frc_vol', frc_vol )
-           CALL iom_get( numrir, 'frc_sal', frc_sal )
-           CALL iom_get( numrir, 'bg_grme', bg_grme )
+           IF(lwp) WRITE(numout,*) ' lim_diahsb_rst read at it= ', kt,' date= ', ndastp
+           IF(lwp) WRITE(numout,*) '~~~~~~~'
+           CALL iom_get( numrir, 'frc_voltop' , frc_voltop  )
+           CALL iom_get( numrir, 'frc_volbot' , frc_volbot  )
+           CALL iom_get( numrir, 'frc_temtop' , frc_temtop  )
+           CALL iom_get( numrir, 'frc_tembot' , frc_tembot  )
+           CALL iom_get( numrir, 'frc_sal'    , frc_sal     )
+           CALL iom_get( numrir, jpdom_autoglo, 'vol_loc_ini', vol_loc_ini )
+           CALL iom_get( numrir, jpdom_autoglo, 'tem_loc_ini', tem_loc_ini )
+           CALL iom_get( numrir, jpdom_autoglo, 'sal_loc_ini', sal_loc_ini )
         ELSE
            IF(lwp) WRITE(numout,*) '~~~~~~~'
            IF(lwp) WRITE(numout,*) ' lim_diahsb at initial state '
            IF(lwp) WRITE(numout,*) '~~~~~~~'
-           frc_vol  = 0._wp                                          
-           frc_sal  = 0._wp                                                 
-           bg_grme  = 0._wp                                       
+           ! set trends to 0
+           frc_voltop  = 0._wp                                          
+           frc_volbot  = 0._wp                                          
+           frc_temtop  = 0._wp                                                 
+           frc_tembot  = 0._wp                                                 
+           frc_sal     = 0._wp                                                 
+           ! record initial ice volume, salt and temp
+           vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:)  ! ice/snow volume (kg/m2)
+           tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                  ! ice/snow heat content (J)
+           sal_loc_ini(:,:) = rhoic * SUM( smv_i(:,:,:), dim=3 )     ! ice salt content (pss*kg/m2)
+           
        ENDIF
 
@@ -288,9 +222,14 @@
         !                                   ! -------------------
         IF(lwp) WRITE(numout,*) '~~~~~~~'
-        IF(lwp) WRITE(numout,*) ' lim_diahsb_rst at it= ', kt,' date= ', ndastp
+        IF(lwp) WRITE(numout,*) ' lim_diahsb_rst write at it= ', kt,' date= ', ndastp
         IF(lwp) WRITE(numout,*) '~~~~~~~'
-        CALL iom_rstput( kt, nitrst, numriw, 'frc_vol'   , frc_vol     )
-        CALL iom_rstput( kt, nitrst, numriw, 'frc_sal'   , frc_sal     )
-        CALL iom_rstput( kt, nitrst, numriw, 'bg_grme'   , bg_grme     )
+        CALL iom_rstput( kt, nitrst, numriw, 'frc_voltop' , frc_voltop  )
+        CALL iom_rstput( kt, nitrst, numriw, 'frc_volbot' , frc_volbot  )
+        CALL iom_rstput( kt, nitrst, numriw, 'frc_temtop' , frc_temtop  )
+        CALL iom_rstput( kt, nitrst, numriw, 'frc_tembot' , frc_tembot  )
+        CALL iom_rstput( kt, nitrst, numriw, 'frc_sal'    , frc_sal     )
+        CALL iom_rstput( kt, nitrst, numriw, 'vol_loc_ini', vol_loc_ini )
+        CALL iom_rstput( kt, nitrst, numriw, 'tem_loc_ini', tem_loc_ini )
+        CALL iom_rstput( kt, nitrst, numriw, 'sal_loc_ini', sal_loc_ini )
         !
      ENDIF
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 6962)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 6963)
@@ -866,5 +866,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
+               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0)))
             END DO
          END DO
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 6962)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 6963)
@@ -110,20 +110,5 @@
       !!---------------------------------------------------------------------
 
-      ! make calls for heat fluxes before it is modified
-      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
-      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface
-      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface
-      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface
-      IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface
-      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice
-      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) )  
-      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   &
-         &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) )
-      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )  
-      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )  
-      IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce(:,:) )   ! emp over ocean (taking into account the snow blown away from the ice)
-      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice(:,:) )   ! emp over ice   (taking into account the snow blown away from the ice)
-
-      ! albedo output
+      ! make call for albedo output before it is modified
       CALL wrk_alloc( jpi,jpj, zalb )    
 
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 6962)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 6963)
@@ -54,5 +54,4 @@
    PUBLIC   lim_var_eqv2glo      
    PUBLIC   lim_var_salprof      
-   PUBLIC   lim_var_icetm        
    PUBLIC   lim_var_bv           
    PUBLIC   lim_var_salprof1d    
@@ -89,23 +88,11 @@
       ! Compute variables
       !--------------------
-      vt_i (:,:) = 0._wp
-      vt_s (:,:) = 0._wp
-      at_i (:,:) = 0._wp
-      ato_i(:,:) = 1._wp
-      !
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume
-               vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume
-               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
-               !
-               rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
-               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch  ! ice thickness
-            END DO
-         END DO
-      END DO
-
+      ! integrated values
+      vt_i (:,:) = SUM( v_i, dim=3 )
+      vt_s (:,:) = SUM( v_s, dim=3 )
+      at_i (:,:) = SUM( a_i, dim=3 )
+      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
+      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
+      !
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -115,26 +102,43 @@
 
       IF( kn > 1 ) THEN
-         et_s (:,:) = 0._wp
-         ot_i (:,:) = 0._wp
+         !
+         ! mean ice/snow thickness
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
+               htm_i(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
+               htm_s(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
+            ENDDO
+         ENDDO
+
+         ! mean temperature (K), salinity and age
          smt_i(:,:) = 0._wp
-         et_i (:,:) = 0._wp
-         !
+         tm_i(:,:)  = 0._wp
+         tm_su(:,:) = 0._wp
+         om_i (:,:) = 0._wp
          DO jl = 1, jpl
+            
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                           ! snow heat content
-                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi20 ) ) 
-                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi20 ) * rswitch   ! ice salinity
-                  rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi20 ) ) 
-                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi20 ) * rswitch   ! ice age
-               END DO
-            END DO
-         END DO
-         !
-         DO jl = 1, jpl
+                  rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
+                  tm_su(ji,jj) = tm_su(ji,jj) + rswitch * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 )
+                  om_i (ji,jj) = om_i (ji,jj) + rswitch *   oa_i(ji,jj,jl)                         / MAX( at_i(ji,jj) , epsi10 )
+               END DO
+            END DO
+            
             DO jk = 1, nlay_i
-               et_i(:,:) = et_i(:,:) + e_i(:,:,jk,jl)       ! ice heat content
-            END DO
-         END DO
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
+                     tm_i(ji,jj)  = tm_i(ji,jj)  + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  &
+                        &            / MAX( vt_i(ji,jj) , epsi10 )
+                     smt_i(ji,jj) = smt_i(ji,jj) + r1_nlay_i * rswitch * s_i(ji,jj,jk,jl) * v_i(ji,jj,jl)  &
+                        &            / MAX( vt_i(ji,jj) , epsi10 )
+                  END DO
+               END DO
+            END DO
+         END DO
+         tm_i  = tm_i  + rt0
+         tm_su = tm_su + rt0
          !
       ENDIF
@@ -246,8 +250,8 @@
       ! Mean temperature
       !-------------------
-      vt_i (:,:) = 0._wp
-      DO jl = 1, jpl
-         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
-      END DO
+      ! integrated values
+      vt_i (:,:) = SUM( v_i, dim=3 )
+      vt_s (:,:) = SUM( v_s, dim=3 )
+      at_i (:,:) = SUM( a_i, dim=3 )
 
       tm_i(:,:) = 0._wp
@@ -397,67 +401,34 @@
    END SUBROUTINE lim_var_salprof
 
-
-   SUBROUTINE lim_var_icetm
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_var_icetm ***
-      !!
-      !! ** Purpose :   computes mean sea ice temperature
+   SUBROUTINE lim_var_bv
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE lim_var_bv ***
+      !!
+      !! ** Purpose :   computes mean brine volume (%) in sea ice
+      !!
+      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
+      !!
+      !! References : Vancoppenolle et al., JGR, 2007
       !!------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
       !!------------------------------------------------------------------
-
-      ! Mean sea ice temperature
-      vt_i (:,:) = 0._wp
-      DO jl = 1, jpl
-         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
-      END DO
-
-      tm_i(:,:) = 0._wp
+      !
+      bvm_i(:,:)   = 0._wp
+      bv_i (:,:,:) = 0._wp
       DO jl = 1, jpl
          DO jk = 1, nlay_i
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
-                  tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  &
-                     &            / MAX( vt_i(ji,jj) , epsi10 )
-               END DO
-            END DO
-         END DO
-      END DO
-      tm_i = tm_i + rt0
-
-   END SUBROUTINE lim_var_icetm
-
-
-   SUBROUTINE lim_var_bv
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_var_bv ***
-      !!
-      !! ** Purpose :   computes mean brine volume (%) in sea ice
-      !!
-      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
-      !!
-      !! References : Vancoppenolle et al., JGR, 2007
-      !!------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) ::   zbvi             ! local scalars
-      !!------------------------------------------------------------------
-      !
-      vt_i (:,:) = 0._wp
-      DO jl = 1, jpl
-         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
-      END DO
-
-      bv_i(:,:) = 0._wp
-      DO jl = 1, jpl
-         DO jk = 1, nlay_i
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  )
-                  zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )   &
-                     &                   * v_i(ji,jj,jl) * r1_nlay_i
-                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi20 ) )  )
-                  bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi20 )
-               END DO
+                  rswitch        = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  )
+                  bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i  &
+                     &                            / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )
+               END DO
+            END DO
+         END DO
+         
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               rswitch      = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
+               bvm_i(ji,jj) = bvm_i(ji,jj) + rswitch * bv_i(ji,jj,jl) * v_i(ji,jj,jl) / MAX( vt_i(ji,jj), epsi10 )
             END DO
          END DO
@@ -715,5 +686,6 @@
             zht_i(ji,1:jpl) = 0._wp
             za_i (ji,1:jpl) = 0._wp
-            
+            itest(:)        = 0      
+      
             ! *** case very thin ice: fill only category 1
             IF ( i_fill == 1 ) THEN
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 6962)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 6963)
@@ -17,5 +17,4 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
-   USE dom_ice
    USE ice
    USE limvar
@@ -40,8 +39,4 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
-#if defined key_dimgout
-# include "limwri_dimg.h90"
-#else
 
    SUBROUTINE lim_wri( kindic )
@@ -59,13 +54,13 @@
       INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices
       REAL(wp) ::  z1_365
-      REAL(wp) ::  ztmp
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zoi, zei, zt_i, zt_s
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, z2da, z2db, zswi    ! 2D workspace
+      REAL(wp) ::  z2da, z2db, ztmp
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zswi2
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, zswi    ! 2D workspace
       !!-------------------------------------------------------------------
 
       IF( nn_timing == 1 )  CALL timing_start('limwri')
 
-      CALL wrk_alloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s )
-      CALL wrk_alloc( jpi, jpj     , z2d, z2da, z2db, zswi )
+      CALL wrk_alloc( jpi, jpj, jpl, zswi2 )
+      CALL wrk_alloc( jpi, jpj     , z2d, zswi )
 
       !-----------------------------
@@ -74,105 +69,59 @@
       z1_365 = 1._wp / 365._wp
 
-      CALL lim_var_icetm      ! mean sea ice temperature
-
-      CALL lim_var_bv         ! brine volume
-
-      DO jj = 1, jpj          ! presence indicator of ice
+      ! brine volume
+      CALL lim_var_bv 
+
+      ! tresholds for outputs
+      DO jj = 1, jpj
          DO ji = 1, jpi
             zswi(ji,jj)  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) )
          END DO
       END DO
-      !
-      !
-      !                                             
-      IF ( iom_use( "icethic_cea" ) ) THEN                       ! mean ice thickness
-         DO jj = 1, jpj 
+      DO jl = 1, jpl
+         DO jj = 1, jpj
             DO ji = 1, jpi
-               z2d(ji,jj)  = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj)
+               zswi2(ji,jj,jl)  = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) )
             END DO
          END DO
-         CALL iom_put( "icethic_cea"  , z2d              )
-      ENDIF
-
-      IF ( iom_use( "snowthic_cea" ) ) THEN                      ! snow thickness = mean snow thickness over the cell 
-         DO jj = 1, jpj                                            
-            DO ji = 1, jpi
-               z2d(ji,jj)  = vt_s(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj)
-            END DO
-         END DO
-         CALL iom_put( "snowthic_cea" , z2d              )       
-      ENDIF
+      END DO
       !
+      ! velocity
       IF ( iom_use( "uice_ipa" ) .OR. iom_use( "vice_ipa" ) .OR. iom_use( "icevel" ) ) THEN 
          DO jj = 2 , jpjm1
             DO ji = 2 , jpim1
-               z2da(ji,jj)  = (  u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp
-               z2db(ji,jj)  = (  v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp
+               z2da  = ( u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp
+               z2db  = ( v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp
+               z2d(ji,jj) = SQRT( z2da * z2da + z2db * z2db )
            END DO
          END DO
-         CALL lbc_lnk( z2da, 'T', -1. )
-         CALL lbc_lnk( z2db, 'T', -1. )
-         CALL iom_put( "uice_ipa"     , z2da             )       ! ice velocity u component
-         CALL iom_put( "vice_ipa"     , z2db             )       ! ice velocity v component
-         DO jj = 1, jpj                                 
-            DO ji = 1, jpi
-               z2d(ji,jj)  = SQRT( z2da(ji,jj) * z2da(ji,jj) + z2db(ji,jj) * z2db(ji,jj) ) 
-            END DO
-         END DO
-         CALL iom_put( "icevel"       , z2d              )       ! ice velocity module
+         CALL lbc_lnk( z2d, 'T', 1. )
+         CALL iom_put( "uice_ipa"     , u_ice      )       ! ice velocity u component
+         CALL iom_put( "vice_ipa"     , v_ice      )       ! ice velocity v component
+         CALL iom_put( "icevel"       , z2d        )       ! ice velocity module
       ENDIF
       !
-      IF ( iom_use( "miceage" ) ) THEN 
-         z2d(:,:) = 0.e0
-         DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch    = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) )
-                  z2d(ji,jj) = z2d(ji,jj) + rswitch * oa_i(ji,jj,jl) / MAX( at_i(ji,jj), 0.1 )
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "miceage"     , z2d * z1_365      )        ! mean ice age
-      ENDIF
-
-      IF ( iom_use( "micet" ) ) THEN 
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               z2d(ji,jj) = ( tm_i(ji,jj) - rt0 ) * zswi(ji,jj)
-            END DO
-         END DO
-         CALL iom_put( "micet"       , z2d               )        ! mean ice temperature
-      ENDIF
+      IF ( iom_use( "miceage" ) )       CALL iom_put( "miceage"     , om_i * zswi * z1_365   )  ! mean ice age
+      IF ( iom_use( "icethic_cea" ) )   CALL iom_put( "icethic_cea" , htm_i * zswi           )  ! ice thickness mean
+      IF ( iom_use( "snowthic_cea" ) )  CALL iom_put( "snowthic_cea", htm_s * zswi           )  ! snow thickness mean
+      IF ( iom_use( "micet" ) )         CALL iom_put( "micet"       , ( tm_i  - rt0 ) * zswi )  ! ice mean    temperature
+      IF ( iom_use( "icest" ) )         CALL iom_put( "icest"       , ( tm_su - rt0 ) * zswi )  ! ice surface temperature
+      IF ( iom_use( "icecolf" ) )       CALL iom_put( "icecolf"     , hicol                  )  ! frazil ice collection thickness
       !
-      IF ( iom_use( "icest" ) ) THEN 
-         z2d(:,:) = 0.e0
-         DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  z2d(ji,jj) = z2d(ji,jj) + zswi(ji,jj) * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi06 )
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "icest"       , z2d              )        ! ice surface temperature
-      ENDIF
-
-      IF ( iom_use( "icecolf" ) )   CALL iom_put( "icecolf", hicol )  ! frazil ice collection thickness
-
       CALL iom_put( "isst"        , sst_m               )        ! sea surface temperature
       CALL iom_put( "isss"        , sss_m               )        ! sea surface salinity
-      CALL iom_put( "iceconc"     , at_i                )        ! ice concentration
-      CALL iom_put( "icevolu"     , vt_i                )        ! ice volume = mean ice thickness over the cell
-      CALL iom_put( "icehc"       , et_i                )        ! ice total heat content
-      CALL iom_put( "isnowhc"     , et_s                )        ! snow total heat content
-      CALL iom_put( "ibrinv"      , bv_i * 100._wp      )        ! brine volume
+      CALL iom_put( "iceconc"     , at_i  * zswi        )        ! ice concentration
+      CALL iom_put( "icevolu"     , vt_i  * zswi        )        ! ice volume = mean ice thickness over the cell
+      CALL iom_put( "icehc"       , et_i  * zswi        )        ! ice total heat content
+      CALL iom_put( "isnowhc"     , et_s  * zswi        )        ! snow total heat content
+      CALL iom_put( "ibrinv"      , bvm_i * zswi * 100. )        ! brine volume
       CALL iom_put( "utau_ice"    , utau_ice            )        ! wind stress over ice along i-axis at I-point
       CALL iom_put( "vtau_ice"    , vtau_ice            )        ! wind stress over ice along j-axis at I-point
       CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation 
-      CALL iom_put( "micesalt"    , smt_i               )        ! mean ice salinity
-
-      CALL iom_put( "icestr"      , strength * 0.001    )        ! ice strength
-      CALL iom_put( "idive"       , divu_i * 1.0e8      )        ! divergence
-      CALL iom_put( "ishear"      , shear_i * 1.0e8     )        ! shear
-      CALL iom_put( "snowvol"     , vt_s                )        ! snow volume
+      CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity
+
+      CALL iom_put( "icestr"      , strength * zswi )    ! ice strength
+      CALL iom_put( "idive"       , divu_i * 1.0e8      )    ! divergence
+      CALL iom_put( "ishear"      , shear_i * 1.0e8     )    ! shear
+      CALL iom_put( "snowvol"     , vt_s   * zswi       )        ! snow volume
       
       CALL iom_put( "icetrp"      , diag_trp_vi * rday  )        ! ice volume transport
@@ -183,10 +132,10 @@
 
       CALL iom_put( "sfxbog"      , sfx_bog * rday      )        ! salt flux from bottom growth
-      CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melt
-      CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from surface melt
+      CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melting
+      CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from surface melting
       CALL iom_put( "sfxsni"      , sfx_sni * rday      )        ! salt flux from snow ice formation
       CALL iom_put( "sfxopw"      , sfx_opw * rday      )        ! salt flux from open water formation
       CALL iom_put( "sfxdyn"      , sfx_dyn * rday      )        ! salt flux from ridging rafting
-      CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from residual
+      CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from limupdate (resultant)
       CALL iom_put( "sfxbri"      , sfx_bri * rday      )        ! salt flux from brines
       CALL iom_put( "sfxsub"      , sfx_sub * rday      )        ! salt flux from sublimation
@@ -202,8 +151,18 @@
       CALL iom_put( "vfxbom"     , wfx_bom * ztmp       )        ! bottom melt 
       CALL iom_put( "vfxice"     , wfx_ice * ztmp       )        ! total ice growth/melt 
+
+      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations  
+         WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog
+         ELSEWHERE                                       ; z2d = 0._wp
+         END WHERE
+         CALL iom_put( "vfxthin", ( wfx_opw + z2d ) * ztmp )
+      ENDIF
+
+      ztmp = rday / rhosn
+      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow)
       CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt 
-      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow) 
-      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow)
-      
+      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow/ice) 
+      CALL iom_put( "vfxsub_err" , wfx_err_sub * ztmp   )        ! "excess" of sublimation sent to ocean      
+ 
       CALL iom_put( "afxtot"     , afx_tot * rday       )        ! concentration tendency (total)
       CALL iom_put( "afxdyn"     , afx_dyn * rday       )        ! concentration tendency (dynamics)
@@ -225,72 +184,24 @@
       CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   !  
       CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   !  
-      CALL iom_put ('hfxtur'     , fhtur(:,:) * SUM(a_i_b(:,:,:), dim=3) ) ! turbulent heat flux at ice base 
+      CALL iom_put ('hfxtur'     , fhtur(:,:) * SUM( a_i_b(:,:,:), dim=3 ) ) ! turbulent heat flux at ice base 
       CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice 
       CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip 
 
-
-      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations  
-         DO jj = 1, jpj 
-            DO ji = 1, jpi
-               z2d(ji,jj)  = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) ! mean ice thickness
-            END DO
-         END DO
-         WHERE( z2d(:,:) < 0.2 .AND. z2d(:,:) > 0. ) ; z2da = wfx_bog
-         ELSEWHERE                                   ; z2da = 0._wp
-         END WHERE
-         CALL iom_put( "vfxthin", ( wfx_opw + z2da ) * ztmp )
-      ENDIF
       
       !--------------------------------
       ! Output values for each category
       !--------------------------------
-      CALL iom_put( "iceconc_cat"      , a_i         )        ! area for categories
-      CALL iom_put( "icethic_cat"      , ht_i        )        ! thickness for categories
-      CALL iom_put( "snowthic_cat"     , ht_s        )        ! snow depth for categories
-      CALL iom_put( "salinity_cat"     , sm_i        )        ! salinity for categories
-
+      IF ( iom_use( "iceconc_cat"  ) )  CALL iom_put( "iceconc_cat"      , a_i   * zswi2   )        ! area for categories
+      IF ( iom_use( "icethic_cat"  ) )  CALL iom_put( "icethic_cat"      , ht_i  * zswi2   )        ! thickness for categories
+      IF ( iom_use( "snowthic_cat" ) )  CALL iom_put( "snowthic_cat"     , ht_s  * zswi2   )        ! snow depth for categories
+      IF ( iom_use( "salinity_cat" ) )  CALL iom_put( "salinity_cat"     , sm_i  * zswi2   )        ! salinity for categories
       ! ice temperature
-      IF ( iom_use( "icetemp_cat" ) ) THEN 
-         zt_i(:,:,:) = SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i
-         CALL iom_put( "icetemp_cat"   , zt_i - rt0  )
-      ENDIF
-      
+      IF ( iom_use( "icetemp_cat"  ) )  CALL iom_put( "icetemp_cat", ( SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i - rt0 ) * zswi2 )
       ! snow temperature
-      IF ( iom_use( "snwtemp_cat" ) ) THEN 
-         zt_s(:,:,:) = SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s
-         CALL iom_put( "snwtemp_cat"   , zt_s - rt0  )
-      ENDIF
-
-      ! Compute ice age
-      IF ( iom_use( "iceage_cat" ) ) THEN 
-         DO jl = 1, jpl 
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) )
-                  rswitch = rswitch * MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - 0.1 ) )
-                  zoi(ji,jj,jl) = oa_i(ji,jj,jl)  / MAX( a_i(ji,jj,jl) , 0.1 ) * rswitch
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "iceage_cat"   , zoi * z1_365 )        ! ice age for categories
-      ENDIF
-
-      ! Compute brine volume
-      IF ( iom_use( "brinevol_cat" ) ) THEN 
-         zei(:,:,:) = 0._wp
-         DO jl = 1, jpl 
-            DO jk = 1, nlay_i
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     rswitch = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) )
-                     zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0 *  &
-                        ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rt0 ), - epsi06 ) ) * &
-                        rswitch * r1_nlay_i
-                  END DO
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "brinevol_cat"     , zei      )        ! brine volume for categories
-      ENDIF
+      IF ( iom_use( "snwtemp_cat"  ) )  CALL iom_put( "snwtemp_cat", ( SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s - rt0 ) * zswi2 )
+      ! ice age
+      IF ( iom_use( "iceage_cat"   ) )  CALL iom_put( "iceage_cat" , o_i * zswi2 * z1_365 )
+      ! brine volume
+      IF ( iom_use( "brinevol_cat" ) )  CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 )
 
       !     !  Create an output files (output.lim.abort.nc) if S < 0 or u > 20 m/s
@@ -298,11 +209,10 @@
       !     not yet implemented
       
-      CALL wrk_dealloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s )
-      CALL wrk_dealloc( jpi, jpj     , z2d, zswi, z2da, z2db )
+      CALL wrk_dealloc( jpi, jpj, jpl, zswi2 )
+      CALL wrk_dealloc( jpi, jpj     , z2d, zswi )
 
       IF( nn_timing == 1 )  CALL timing_stop('limwri')
       
    END SUBROUTINE lim_wri
-#endif
 
  
@@ -319,78 +229,100 @@
       !!   4.0  !  2013-06  (C. Rousset)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt               ! ocean time-step index)
-      INTEGER, INTENT( in ) ::   kid , kh_i       
+      INTEGER, INTENT( in )   ::   kt               ! ocean time-step index)
+      INTEGER, INTENT( in )   ::   kid , kh_i
+      INTEGER                 ::   nz_i, jl
+      REAL(wp), DIMENSION(jpl) :: jcat
       !!----------------------------------------------------------------------
-
-      CALL histdef( kid, "iicethic", "Ice thickness"           , "m"      ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iiceconc", "Ice concentration"       , "%"      ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicetemp", "Ice temperature"         , "C"      ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)"   , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)"   , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa",   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa",   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicesflx", "Solar flux over ocean"     , "w/m2" ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2" ,   &
+      DO jl = 1, jpl
+         jcat(jl) = REAL(jl)
+      ENDDO
+      
+      CALL histvert( kid, "ncatice", "Ice Categories","", jpl, jcat, nz_i, "up")
+
+      CALL histdef( kid, "sithic", "Ice thickness"           , "m"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "siconc", "Ice concentration"       , "%"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sitemp", "Ice temperature"         , "C"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sivelu", "i-Ice speed "            , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sivelv", "j-Ice speed "            , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sistru", "i-Wind stress over ice " , "Pa"     ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sistrv", "j-Wind stress over ice " , "Pa"     ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sisflx", "Solar flux over ocean"     , "w/m2" ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sinflx", "Non-solar flux over ocean" , "w/m2" ,   &
       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
       CALL histdef( kid, "isnowpre", "Snow precipitation"      , "kg/m2/s",   &
       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicesali", "Ice salinity"            , "PSU"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicevolu", "Ice volume"              , "m"      ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicedive", "Ice divergence"          , "10-8s-1",   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicebopr", "Ice bottom production"   , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicedypr", "Ice dynamic production"  , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicelapr", "Ice open water prod"     , "m/s"    ,   &
+      CALL histdef( kid, "sisali", "Ice salinity"            , "PSU"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sivolu", "Ice volume"              , "m"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "sidive", "Ice divergence"          , "10-8s-1",   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+
+      CALL histdef( kid, "vfxbog", "Ice bottom production"   , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxdyn", "Ice dynamic production"  , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxopw", "Ice open water prod"     , "m/s"    ,   &
       &       jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicesipr", "Snow ice production "    , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicebome", "Ice bottom melt"         , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicesume", "Ice surface melt"        , "m/s"    ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iisfxdyn", "Salt flux from dynmics"  , ""       ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iisfxres", "Salt flux from limupdate", ""       ,   &
-      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxsni", "Snow ice production "    , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxres", "Ice prod from limupdate" , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxbom", "Ice bottom melt"         , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "vfxsum", "Ice surface melt"        , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+
+      CALL histdef( kid, "sithicat", "Ice thickness"         , "m"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "siconcat", "Ice concentration"     , "%"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sisalcat", "Ice salinity"           , ""      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sitemcat", "Ice temperature"       , "C"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "snthicat", "Snw thickness"         , "m"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "sntemcat", "Snw temperature"       , "C"      ,   &
+      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
 
       CALL histend( kid, snc4set )   ! end of the file definition
 
-      CALL histwrite( kid, "iicethic", kt, icethi        , jpi*jpj, (/1/) )    
-      CALL histwrite( kid, "iiceconc", kt, at_i          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicetemp", kt, tm_i - rt0    , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevelu", kt, u_ice          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevelv", kt, v_ice          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicestru", kt, utau_ice       , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicestrv", kt, vtau_ice       , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesflx", kt, qsr , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicenflx", kt, qns , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sithic", kt, htm_i         , jpi*jpj, (/1/) )    
+      CALL histwrite( kid, "siconc", kt, at_i          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sitemp", kt, tm_i - rt0    , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sivelu", kt, u_ice          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sivelv", kt, v_ice          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sistru", kt, utau_ice       , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sistrv", kt, vtau_ice       , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sisflx", kt, qsr , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sinflx", kt, qns , jpi*jpj, (/1/) )
       CALL histwrite( kid, "isnowpre", kt, sprecip        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesali", kt, smt_i          , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicevolu", kt, vt_i           , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) )
-
-      CALL histwrite( kid, "iicebopr", kt, wfx_bog        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicedypr", kt, wfx_dyn        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicelapr", kt, wfx_opw        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesipr", kt, wfx_sni        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicerepr", kt, wfx_res        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicebome", kt, wfx_bom        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iicesume", kt, wfx_sum        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iisfxdyn", kt, sfx_dyn        , jpi*jpj, (/1/) )
-      CALL histwrite( kid, "iisfxres", kt, sfx_res        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sisali", kt, smt_i          , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sivolu", kt, vt_i           , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "sidive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) )
+
+      CALL histwrite( kid, "vfxbog", kt, wfx_bog        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxdyn", kt, wfx_dyn        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxopw", kt, wfx_opw        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxsni", kt, wfx_sni        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxres", kt, wfx_res        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxbom", kt, wfx_bom        , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxsum", kt, wfx_sum        , jpi*jpj, (/1/) )
+
+      CALL histwrite( kid, "sithicat", kt, ht_i        , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "siconcat", kt, a_i         , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "sisalcat", kt, sm_i        , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "sitemcat", kt, tm_i - rt0  , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "snthicat", kt, ht_s        , jpi*jpj*jpl, (/1/) )    
+      CALL histwrite( kid, "sntemcat", kt, tm_su - rt0 , jpi*jpj*jpl, (/1/) )    
 
       ! Close the file
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 6962)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 6963)
@@ -38,5 +38,4 @@
    PUBLIC   dia_hsb        ! routine called by step.F90
    PUBLIC   dia_hsb_init   ! routine called by nemogcm.F90
-   PUBLIC   dia_hsb_rst    ! routine called by step.F90
 
    LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets
@@ -86,4 +85,5 @@
       !!---------------------------------------------------------------------------
       IF( nn_timing == 1 )   CALL timing_start('dia_hsb')      
+      !
       CALL wrk_alloc( jpi,jpj,   z2d0, z2d1 )
       !
@@ -174,5 +174,7 @@
       ENDDO
 
-      ! Substract forcing from heat content, salt content and volume variations
+      ! ------------------------ !
+      ! 3 -  Drifts              !
+      ! ------------------------ !
       zdiff_v1 = zdiff_v1 - frc_v
       IF( lk_vvl )   zdiff_v2 = zdiff_v2 - frc_v
@@ -187,5 +189,5 @@
 
       ! ----------------------- !
-      ! 3 - Diagnostics writing !
+      ! 4 - Diagnostics writing !
       ! ----------------------- !
       zvol_tot = 0._wp                    ! total ocean volume (calculated with scale factors)
@@ -200,23 +202,27 @@
 !!gm end
 
+      CALL iom_put(   'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3) 
+      CALL iom_put(   'bgfrctem' , frc_t    * rau0 * rcp * 1.e-20 )   ! hc  - surface forcing (1.e20 J) 
+      CALL iom_put(   'bgfrchfx' , frc_t    * rau0 * rcp /  &         ! hc  - surface forcing (W/m2) 
+         &                       ( surf_tot * kt * rdt )        )
+      CALL iom_put(   'bgfrcsal' , frc_s    * 1.e-9    )              ! sc  - surface forcing (psu*km3) 
+
       IF( lk_vvl ) THEN
-        CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature variation (C) 
-        CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot )              ! Salinity    variation (psu)
-        CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rau0 * rcp )   ! Heat content variation (1.e20 J) 
-        CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9    )              ! Salt content variation (psu*km3)
-        CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh variation (km3)  
-        CALL iom_put( 'bgvole3t' , zdiff_v2 * 1.e-9    )              ! volume e3t variation (km3)  
-        CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3) 
-        CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C) 
-        CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot    )              ! sc  - surface forcing (psu) 
+        CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature drift     (C) 
+        CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot )              ! Salinity    drift     (pss)
+        CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rau0 * rcp )   ! Heat content drift    (1.e20 J) 
+        CALL iom_put( 'bgheatfx' , zdiff_hc * rau0 * rcp /  &         ! Heat flux drift       (W/m2) 
+           &                       ( surf_tot * kt * rdt )        )
+        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)  
       ELSE
-        CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content variation (C) 
-        CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content variation (psu)
-        CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content variation (1.e20 J) 
-        CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content variation (psu*km3)
-        CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh variation (km3)  
-        CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3) 
-        CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C) 
-        CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot    )              ! sc  - surface forcing (psu) 
+        CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content drift    (C) 
+        CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content drift    (pss)
+        CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content drift    (1.e20 J) 
+        CALL iom_put( 'bgheatfx' , zdiff_hc1 * rau0 * rcp /  &        ! Heat flux drift       (W/m2) 
+           &                       ( surf_tot * kt * rdt )         )
+        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)
@@ -244,10 +250,8 @@
      !
      INTEGER ::   ji, jj, jk   ! dummy loop indices
-     INTEGER ::   id1          ! local integers
      !!----------------------------------------------------------------------
      !
      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
         IF( ln_rstart ) THEN                   !* Read the restart file
-           !id1 = iom_varid( numror, 'frc_vol'  , ldstop = .FALSE. )
            !
            IF(lwp) WRITE(numout,*) '~~~~~~~'
@@ -261,11 +265,11 @@
               CALL iom_get( numror, 'frc_wn_s', frc_wn_s )
            ENDIF
-           CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini )
-           CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini )
-           CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini )
-           CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini )
+           CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini(:,:) )
+           CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini(:,:,:) )
+           CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini(:,:,:) )
+           CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini(:,:,:) )
            IF( .NOT. lk_vvl ) THEN
-              CALL iom_get( numror, jpdom_autoglo, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
-              CALL iom_get( numror, jpdom_autoglo, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
+              CALL iom_get( numror, jpdom_autoglo, 'ssh_hc_loc_ini', ssh_hc_loc_ini(:,:) )
+              CALL iom_get( numror, jpdom_autoglo, 'ssh_sc_loc_ini', ssh_sc_loc_ini(:,:) )
            ENDIF
        ELSE
@@ -312,12 +316,13 @@
            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s )
         ENDIF
-        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, 'hc_loc_ini', hc_loc_ini )
-        CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini )
+        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, 'hc_loc_ini', hc_loc_ini(:,:,:) )
+        CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini(:,:,:) )
         IF( .NOT. lk_vvl ) 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 )
+           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
@@ -338,5 +343,4 @@
       !!             - Compute coefficients for conversion
       !!---------------------------------------------------------------------------
-      INTEGER ::   jk       ! dummy loop indice
       INTEGER ::   ierror   ! local integer
       INTEGER ::   ios
@@ -344,10 +348,4 @@
       NAMELIST/namhsb/ ln_diahsb
       !!----------------------------------------------------------------------
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
-         WRITE(numout,*) '~~~~~~~~ '
-      ENDIF
 
       REWIND( numnam_ref )              ! Namelist namhsb in reference namelist
@@ -360,14 +358,11 @@
       IF(lwm) WRITE ( numond, namhsb )
 
-      !
-      IF(lwp) THEN                   ! Control print
+      IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namhsb : set hsb parameters'
-         WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb
-         WRITE(numout,*)
-      ENDIF
-
+         WRITE(numout,*) 'dia_hsb_init'
+         WRITE(numout,*) '~~~~~~~~ '
+         WRITE(numout,*) '  check the heat and salt budgets (T) or not (F)       ln_diahsb = ', ln_diahsb
+      ENDIF
+      !
       IF( .NOT. ln_diahsb )   RETURN
          !      IF( .NOT. lk_mpp_rep ) &
@@ -382,10 +377,14 @@
          &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror )
       IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
-      ENDIF
-
-      IF(.NOT. lk_vvl ) 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 hc_loc_ini' )   ;   RETURN
+         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )
+         RETURN
+      ENDIF
+
+      IF( .NOT. lk_vvl ) THEN
+         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 hc_loc_ini' )
+            RETURN
+         ENDIF
       ENDIF
 
@@ -393,8 +392,6 @@
       ! 2 - Time independant variables and file opening !
       ! ----------------------------------------------- !
-      IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"
-      IF(lwp) WRITE(numout,*) '~~~~~~~'
       surf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)      ! masked surface grid cell area
-      surf_tot  = glob_sum( surf(:,:) )                                       ! total ocean surface area
+      surf_tot  = glob_sum( surf(:,:) )                   ! total ocean surface area
 
       IF( lk_bdy ) CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )         
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6962)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 6963)
@@ -1607,10 +1607,9 @@
       ENDIF
 
-      !! clem: we should output qemp_oce and qemp_ice (at least)
-      IF( iom_use('hflx_snow_cea') )   CALL iom_put( 'hflx_snow_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) )   ! heat flux from snow (cell average)
-      !! these diags are not outputed yet
-!!      IF( iom_use('hflx_rain_cea') )   CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptn(:,:) )   ! heat flux from rain (cell average)
-!!      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put( 'hflx_snow_ao_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * (1._wp - zsnw(:,:)) ) ! heat flux from snow (cell average)
-!!      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put( 'hflx_snow_ai_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * zsnw(:,:) ) ! heat flux from snow (cell average)
+      ! some more outputs
+      IF( iom_use('hflx_snow_cea') )    CALL iom_put('hflx_snow_cea',   sprecip(:,:) * ( zcptn(:,:) - Lfus ) )                       ! heat flux from snow (cell average)
+      IF( iom_use('hflx_rain_cea') )    CALL iom_put('hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptn(:,:) )                 ! heat flux from rain (cell average)
+      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea',sprecip(:,:) * ( zcptn(:,:) - Lfus ) * (1._wp - zsnw(:,:)) ) ! heat flux from snow (cell average)
+      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea',sprecip(:,:) * ( zcptn(:,:) - Lfus ) * zsnw(:,:) )           ! heat flux from snow (cell average)
 
 #else
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 6962)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 6963)
@@ -229,5 +229,5 @@
          CALL lim_sbc_flx( kt )                     ! Update surface ocean mass, heat and salt fluxes
          !
-         IF(ln_limdiaout) CALL lim_diahsb           ! Diagnostics and outputs 
+         IF(ln_limdiaout) CALL lim_diahsb( kt )     ! Diagnostics and outputs 
          !
          CALL lim_wri( 1 )                          ! Ice outputs 
@@ -310,8 +310,10 @@
          numit = nit000 - 1
       ENDIF
-      CALL lim_var_agg(1)
+      CALL lim_var_agg(2)
       CALL lim_var_glo2eqv
       !
       CALL lim_sbc_init                 ! ice surface boundary condition   
+      !
+      IF( ln_limdiaout) CALL lim_diahsb_init  ! initialization for diags
       !
       fr_i(:,:)     = at_i(:,:)         ! initialisation of sea-ice fraction
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6962)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6963)
@@ -337,4 +337,6 @@
       IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
       !
+      IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics
+
       IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file
       IF( ln_sto_eos       )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters
@@ -351,6 +353,5 @@
       ENDIF
 #endif
-      IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics
-      IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
+      IF( lk_diaobs        )   CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
