Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/field_def.xml
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/field_def.xml	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/field_def.xml	(revision 4649)
@@ -152,5 +152,4 @@
 	
          <!-- *_oce variables available with ln_blk_clio or ln_blk_core -->
-         <field id="qns_oce"      long_name="Non solar Downward Heat Flux over open ocean"                 unit="W/m2"     />
          <field id="qlw_oce"      long_name="Longwave Downward Heat Flux over open ocean"                  unit="W/m2"     />
          <field id="qsb_oce"      long_name="Sensible Downward Heat Flux over open ocean"                  unit="W/m2"     />
@@ -222,8 +221,10 @@
          <field id="qsr_oce"      long_name="solar heat flux at ocean surface"                             unit="W/m2"     />
          <field id="qns_oce"      long_name="non-solar heat flux at ocean surface"                         unit="W/m2"     />
+         <field id="qt_ice"       long_name="total flux at ice surface"                                    unit="W/m2"     />
+         <field id="qsr_ice"      long_name="solar heat flux at ice surface"                               unit="W/m2"     />
+         <field id="qns_ice"      long_name="non-solar heat flux at ice surface"                           unit="W/m2"     />
+         <field id="qtr_ice"      long_name="solar heat flux transmitted thru the ice"                     unit="W/m2"     />
          <field id="utau_ice"     long_name="Wind stress along i-axis over the ice at i-point"             unit="N/m2"     />
          <field id="vtau_ice"     long_name="Wind stress along j-axis over the ice at i-point"             unit="N/m2"     />
-	 <field id="qsr_io"       long_name="Ice-Oce downward solar heat flux"                             unit="W/m2"     />
-	 <field id="qns_io"       long_name="Ice-Oce downward non-solar heat flux"                         unit="W/m2"     />
          <field id="micesalt"     long_name="Mean ice salinity"                                            unit="psu"      />
          <field id="miceage"      long_name="Mean ice age"                                                 unit="years"    />
@@ -237,5 +238,5 @@
 
          <field id="micet"        long_name="Mean ice temperature"                                         unit="degC"     />
-         <field id="icehc"        long_name="ice total heat content"                                       unit="10^9 J"   /> 
+         <field id="icehc"        long_name="ice total heat content"                                       unit="10^9J"   /> 
          <field id="isnowhc"      long_name="snow total heat content"                                      unit="10^9J"    />
          <field id="icest"        long_name="ice surface temperature"                                      unit="degC"     />
@@ -251,6 +252,6 @@
          <field id="icetrp"       long_name="ice volume transport"                                         unit="m/day"   />
          <field id="snwtrp"       long_name="snw volume transport"                                         unit="m/day"   />
-         <field id="deitrp"       long_name="advected ice enhalpy"                                         unit="W/2"   />
-         <field id="destrp"       long_name="advected snw enhalpy"                                         unit="W/2"   />
+         <field id="deitrp"       long_name="advected ice enhalpy"                                         unit="W/m2"   />
+         <field id="destrp"       long_name="advected snw enhalpy"                                         unit="W/m2"   />
 
          <field id="sfxbri"       long_name="brine salt flux"                                              unit="psu*kg/m2/day" />
@@ -274,20 +275,27 @@
          <field id="vfxsnw"       long_name="snw melt/growth"                                              unit="m/day"   />
          <field id="vfxsub"       long_name="snw sublimation"                                              unit="m/day"   />
-
-         <field id="hfxdhc1"   long_name="Heat content variation in snow and ice"   unit="W/m2" />
-         <field id="hfxspr"    long_name="Heat content of snow precip"              unit="W/m2" />
-         <field id="hfxqsr"    long_name="solar fluxes given to ocean"             unit="W/m2"  />
-         <field id="hfxqns"    long_name="non solar fluxes given to ocean"         unit="W/m2"  />
-
-         <field id="hfxthd"   long_name="heat fluxes from ice-ocean exchange during thermo"              unit="W/m2"  />
-         <field id="hfxdyn"   long_name="heat fluxes from ice-ocean exchange during dynamic"             unit="W/m2"  />
-         <field id="hfxres"   long_name="heat fluxes from ice-ocean exchange during resultant"           unit="W/m2"  />
-         <field id="hfxsnw"   long_name="heat fluxes from snow-ocean exchange"                           unit="W/m2"  />
-         <field id="hfxsub"   long_name="heat fluxes from sublimation"                                   unit="W/m2"  />
-         <field id="hfxerr"   long_name="heat fluxes error after heat diffusion"                         unit="W/m2"  />
-         <field id="hfxerr_rem" long_name="heat fluxes error after remapping"                         unit="W/m2"  />
-         <field id="hfxtot"   long_name="heat fluxes total used by ice"                                  unit="W/m2"  />
-         <field id="hfxout"   long_name="non solar heat fluxes received by the ocean"             unit="W/m2"  />
-         <field id="hfxin"    long_name="total     heat fluxes at the ice surface"             unit="W/m2"  />
+         <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="m/day"   />
+
+         <field id="hfxsum"   long_name="heat fluxes causing surface ice melt"            unit="W/m2"  />
+         <field id="hfxbom"   long_name="heat fluxes causing bottom ice melt"             unit="W/m2"  />
+         <field id="hfxbog"   long_name="heat fluxes causing bottom ice growth"           unit="W/m2"  />
+         <field id="hfxdif"   long_name="heat fluxes causing ice temperature change"      unit="W/m2"  />
+         <field id="hfxopw"   long_name="heat fluxes causing open water ice formation"    unit="W/m2"  />
+         <field id="hfxsnw"   long_name="heat fluxes causing snow melt"                   unit="W/m2"  />
+         <field id="hfxerr"   long_name="heat fluxes error after heat diffusion"          unit="W/m2"  />
+         <field id="hfxerr_rem" long_name="heat fluxes error after remapping"             unit="W/m2"  />
+         <field id="hfxout"   long_name="non solar heat fluxes received by the ocean"     unit="W/m2"  />
+         <field id="hfxin"    long_name="total     heat fluxes at the ice surface"        unit="W/m2"  />
+
+	 <!-- heat flux associated with mass exchange -->
+         <field id="hfxthd"   long_name="heat fluxes from ice-ocean mass exchange during thermo"              unit="W/m2"  />
+         <field id="hfxdyn"   long_name="heat fluxes from ice-ocean mass exchange during dynamic"             unit="W/m2"  />
+         <field id="hfxres"   long_name="heat fluxes from ice-ocean mass exchange during resultant"           unit="W/m2"  />
+         <field id="hfxsub"   long_name="heat fluxes from ice-atm. mass exchange during sublimation"          unit="W/m2"  />
+         <field id="hfxspr"   long_name="heat fluxes from ice-atm. mass exchange during snow precip"          unit="W/m2" />
+
+	 <!-- diags -->
+         <field id="hfxdhc"    long_name="Heat content variation in snow and ice"   unit="W/m2" />
+         <field id="hfxtur"    long_name="turbulent heat flux at the ice base"      unit="W/m2"  />
 
       </field_group>
@@ -403,4 +411,7 @@
 	 <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="psu*m/day"   />
@@ -415,15 +426,17 @@
 
 
-        <field id="ibghfxdhc1"   long_name="Heat content variation in snow and ice"   unit="W" />
+        <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="ibghfxqsr"    long_name="solar fluxes given to ocean"             unit="W"  />
-        <field id="ibghfxqns"    long_name="non solar fluxes given to ocean"         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="ibghfxtot"   long_name="heat fluxes total used by ice"                                  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"  />
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml	(revision 4649)
@@ -53,6 +53,11 @@
 	<file id="file2" name_suffix="_SBC" description="surface fluxes variables" > <!-- time step automaticaly defined based on nn_fsbc -->
 	  <field field_ref="empmr"        name="wfo"      long_name="water_flux_into_sea_water"                     />
-	  <field field_ref="qsr"          name="hfxqsr"   long_name="surface_net_downward_shortwave_flux"           />
-	  <field field_ref="qt"           name="hfxtot"   long_name="surface_net_downward_total_heat_flux"          />
+	  <field field_ref="qsr_oce"      name="qsr_oce"  long_name="downward shortwave flux at ocean surface"           />
+	  <field field_ref="qns_oce"      name="qns_oce"  long_name="downward non solar flux at ocean surface"           />
+	  <field field_ref="qt_oce"       name="qt_oce"   long_name="downward total flux at ocean surface"           />
+	  <field field_ref="qsr_ice"      name="qsr_ice"  long_name="downward shortwave flux at ice surface"           />
+	  <field field_ref="qns_ice"      name="qns_ice"  long_name="downward non solar flux at ice surface"           />
+	  <field field_ref="qtr_ice"      name="qtr_ice"  long_name="shortwave flux transmitted thru the ice"           />
+	  <field field_ref="qt_ice"       name="qt_ice"   long_name="downward total flux at ice surface"           />
 	  <field field_ref="saltflx"      name="sfx"  />
 	  <field field_ref="taum"         name="taum" />
@@ -63,8 +68,4 @@
           <field field_ref="utau_ice"         name="utau_ice" />
           <field field_ref="vtau_ice"         name="vtau_ice" />
-	  <!-- clem
-          <field field_ref="qsr_io"           name="iicesflx" />
-          <field field_ref="qns_io"           name="iicenflx" />
-	  -->
 
 	</file>
@@ -97,5 +98,5 @@
 
 	<file id="file6" name_suffix="_icemod" description="ice variables" enabled=".true." >
-	  <field field_ref="snowthic_cea"     name="sndept"     long_name="surface_snow_thickness"   />
+	  <field field_ref="snowthic_cea"     name="snthic"     long_name="surface_snow_thickness"   />
 	  <field field_ref="icethic_cea"      name="sithic"     long_name="sea_ice_thickness"        />
           <field field_ref="icevolu"          name="sivolu" />
@@ -113,4 +114,5 @@
           <field field_ref="vfxsnw"          name="vfxsnw" />
           <field field_ref="vfxsub"          name="vfxsub" />
+          <field field_ref="vfxspr"          name="vfxspr" />
 
           <field field_ref="icetrp"          name="sivtrp" />
@@ -129,19 +131,26 @@
           <field field_ref="sfx"              name="sfx" />
 
-          <field field_ref="hfxdhc1"          name="hfxdhc1"    />
-          <field field_ref="hfxspr"           name="hfxspr"    />
-          <field field_ref="hfxqsr"           name="hfxqsr"    />
-          <field field_ref="hfxqns"           name="hfxqns"    />
-
-          <field field_ref="hfxthd"          name="hfxthd"    />
+          <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="hfxout"          name="hfxout"    />
-          <field field_ref="hfxin"           name="hfxin"    />
-          <field field_ref="hfxtot"          name="hfxtot"    />
-          <field field_ref="hfxsnw"          name="hfxsnw"    />
+          <field field_ref="hfxthd"          name="hfxthd"    />
+	  <!-- ice-atm. heat flux from mass exchange -->
           <field field_ref="hfxsub"          name="hfxsub"    />
-          <field field_ref="hfxerr"          name="hfxerr"    />
-          <field field_ref="hfxerr_rem"      name="hfxerr_rem"    />
+          <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" />
@@ -167,7 +176,7 @@
           <field field_ref="iceconc_cat"      name="siconcat"/>
           <field field_ref="icethic_cat"      name="sithicat"/>
-          <field field_ref="snowthic_cat"     name="sndeptcat"/>
+          <field field_ref="snowthic_cat"     name="snthicat"/>
           <field field_ref="salinity_cat"     name="salincat"/>
-          <field field_ref="brinevol_cat"     name="sibrincat"/>
+          <field field_ref="brinevol_cat"     name="sibricat"/>
 
 	</file>
@@ -191,8 +200,5 @@
           <field field_ref="bgfrctem"     name="bgfrctem"    />
           <field field_ref="bgfrcsal"     name="bgfrcsal"    />
-	  <!--
-          <field field_ref="bgmistem"     name="bgmistem"    />
-          <field field_ref="bgmissal"     name="bgmissal"    />
-	  -->
+
         </file>
 
@@ -215,4 +221,7 @@
           <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"     />
@@ -226,8 +235,6 @@
           <field field_ref="ibgsfxsum"    name="ibgsfxsum"      />
 
-          <field field_ref="ibghfxdhc1"    name="ibghfxdhc1"    />
+          <field field_ref="ibghfxdhc"    name="ibghfxdhc"    />
           <field field_ref="ibghfxspr"    name="ibghfxspr"    />
-          <field field_ref="ibghfxqsr"    name="ibghfxqsr"    />
-          <field field_ref="ibghfxqns"    name="ibghfxqns"    />
 
           <field field_ref="ibghfxres"    name="ibghfxres"    />
@@ -235,7 +242,11 @@
           <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="ibghfxtot"    name="ibghfxtot"    />
           <field field_ref="ibghfxsnw"    name="ibghfxsnw"    />
 
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef_default.xml
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef_default.xml	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef_default.xml	(revision 4649)
@@ -53,6 +53,11 @@
 	<file id="file2" name_suffix="_SBC" description="surface fluxes variables" > <!-- time step automaticaly defined based on nn_fsbc -->
 	  <field field_ref="empmr"        name="wfo"      long_name="water_flux_into_sea_water"                     />
-	  <field field_ref="qsr"          name="hfxqsr"   long_name="surface_net_downward_shortwave_flux"           />
-	  <field field_ref="qt"           name="hfxtot"   long_name="surface_net_downward_total_heat_flux"          />
+	  <field field_ref="qsr_oce"      name="qsr_oce"  long_name="downward shortwave flux at ocean surface"           />
+	  <field field_ref="qns_oce"      name="qns_oce"  long_name="downward non solar flux at ocean surface"           />
+	  <field field_ref="qt_oce"       name="qt_oce"   long_name="downward total flux at ocean surface"           />
+	  <field field_ref="qsr_ice"      name="qsr_ice"  long_name="downward shortwave flux at ice surface"           />
+	  <field field_ref="qns_ice"      name="qns_ice"  long_name="downward non solar flux at ice surface"           />
+	  <field field_ref="qtr_ice"      name="qtr_ice"  long_name="shortwave flux transmitted thru the ice"           />
+	  <field field_ref="qt_ice"       name="qt_ice"   long_name="downward total flux at ice surface"           />
 	  <field field_ref="saltflx"      name="sfx"  />
 	  <field field_ref="taum"         name="taum" />
@@ -63,8 +68,4 @@
           <field field_ref="utau_ice"         name="utau_ice" />
           <field field_ref="vtau_ice"         name="vtau_ice" />
-	  <!-- clem
-          <field field_ref="qsr_io"           name="iicesflx" />
-          <field field_ref="qns_io"           name="iicenflx" />
-	  -->
 
 	</file>
@@ -97,5 +98,5 @@
 
 	<file id="file6" name_suffix="_icemod" description="ice variables" enabled=".true." >
-	  <field field_ref="snowthic_cea"     name="sndept"     long_name="surface_snow_thickness"   />
+	  <field field_ref="snowthic_cea"     name="snthic"     long_name="surface_snow_thickness"   />
 	  <field field_ref="icethic_cea"      name="sithic"     long_name="sea_ice_thickness"        />
           <field field_ref="icevolu"          name="sivolu" />
@@ -113,4 +114,5 @@
           <field field_ref="vfxsnw"          name="vfxsnw" />
           <field field_ref="vfxsub"          name="vfxsub" />
+          <field field_ref="vfxspr"          name="vfxspr" />
 
           <field field_ref="icetrp"          name="sivtrp" />
@@ -129,19 +131,26 @@
           <field field_ref="sfx"              name="sfx" />
 
-          <field field_ref="hfxdhc1"          name="hfxdhc1"    />
-          <field field_ref="hfxspr"           name="hfxspr"    />
-          <field field_ref="hfxqsr"           name="hfxqsr"    />
-          <field field_ref="hfxqns"           name="hfxqns"    />
-
-          <field field_ref="hfxthd"          name="hfxthd"    />
+          <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="hfxout"          name="hfxout"    />
-          <field field_ref="hfxin"           name="hfxin"    />
-          <field field_ref="hfxtot"          name="hfxtot"    />
-          <field field_ref="hfxsnw"          name="hfxsnw"    />
+          <field field_ref="hfxthd"          name="hfxthd"    />
+	  <!-- ice-atm. heat flux from mass exchange -->
           <field field_ref="hfxsub"          name="hfxsub"    />
-          <field field_ref="hfxerr"          name="hfxerr"    />
-          <field field_ref="hfxerr_rem"      name="hfxerr_rem"    />
+          <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" />
@@ -167,7 +176,7 @@
           <field field_ref="iceconc_cat"      name="siconcat"/>
           <field field_ref="icethic_cat"      name="sithicat"/>
-          <field field_ref="snowthic_cat"     name="sndeptcat"/>
+          <field field_ref="snowthic_cat"     name="snthicat"/>
           <field field_ref="salinity_cat"     name="salincat"/>
-          <field field_ref="brinevol_cat"     name="sibrincat"/>
+          <field field_ref="brinevol_cat"     name="sibricat"/>
 
 	</file>
@@ -191,8 +200,5 @@
           <field field_ref="bgfrctem"     name="bgfrctem"    />
           <field field_ref="bgfrcsal"     name="bgfrcsal"    />
-	  <!--
-          <field field_ref="bgmistem"     name="bgmistem"    />
-          <field field_ref="bgmissal"     name="bgmissal"    />
-	  -->
+
         </file>
 
@@ -215,4 +221,7 @@
           <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"     />
@@ -226,8 +235,6 @@
           <field field_ref="ibgsfxsum"    name="ibgsfxsum"      />
 
-          <field field_ref="ibghfxdhc1"    name="ibghfxdhc1"    />
+          <field field_ref="ibghfxdhc"    name="ibghfxdhc"    />
           <field field_ref="ibghfxspr"    name="ibghfxspr"    />
-          <field field_ref="ibghfxqsr"    name="ibghfxqsr"    />
-          <field field_ref="ibghfxqns"    name="ibghfxqns"    />
 
           <field field_ref="ibghfxres"    name="ibghfxres"    />
@@ -235,7 +242,11 @@
           <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="ibghfxtot"    name="ibghfxtot"    />
           <field field_ref="ibghfxsnw"    name="ibghfxsnw"    />
 
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice	(revision 4649)
@@ -15,12 +15,12 @@
 !-----------------------------------------------------------------------
    cn_icerst_in  = "restart_ice"   !  suffix of ice restart name (input)
-   cn_icerst_out = "restart_ice"      !  suffix of ice restart name (output)
-   ln_limdyn   = .true.    !  ice dynamics (T) or thermodynamics only (F)
-   amax        = 0.999      !  maximum ice concentration
-   cai         =  1.40e-3  !  atmospheric drag over sea ice
-   cao         =  1.00e-3  !  atmospheric drag over ocean
-   ln_nicep    = .false.   !  Ice points output for debug (yes or no)
-   ln_limdiahsb  = .false.    !  check the heat and salt budgets (T) or not (F)
-   ln_limdiaout  = .true.    !  output the heat and salt budgets (T) or not (F)
+   cn_icerst_out = "restart_ice"   !  suffix of ice restart name (output)
+   ln_limdyn     = .true.          !  ice dynamics (T) or thermodynamics only (F)
+   amax          = 0.999           !  maximum ice concentration
+   cai           = 1.40e-3         !  atmospheric drag over sea ice (clio)
+   cao           = 1.00e-3         !  atmospheric drag over ocean   (clio)
+   ln_nicep      = .false.         !  Ice points output for debug (yes or no)
+   ln_limdiahsb  = .false.          !  check the heat and salt budgets (T) or not (F)
+   ln_limdiaout  = .true.          !  output the heat and salt budgets (T) or not (F)
 /
 !-----------------------------------------------------------------------
@@ -147,15 +147,4 @@
 !&namicehsb       !  Heat and salt budgets 
 !!-----------------------------------------------------------------------
-!
 !/
-!-----------------------------------------------------------------------
-&namiceout     !   parameters for outputs
-!-----------------------------------------------------------------------
-!   noumef      =   43      !  number of fields
-!SF   add_diag_swi=    1      !  1 -> diagnose distribution in thickness space
-!SF                           !  0 -> only simple diagnostics
-!
-!           !         title of the field           !  name     !   units   !  save  ! multipl. ! additive !
-!           !                                      !           !           ! or not !  factor  !  factor  !
-/ 
 
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice_lim3
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice_lim3	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice_lim3	(revision 4649)
@@ -15,12 +15,12 @@
 !-----------------------------------------------------------------------
    cn_icerst_in  = "restart_ice"   !  suffix of ice restart name (input)
-   cn_icerst_out = "restart_ice"      !  suffix of ice restart name (output)
-   ln_limdyn   = .true.    !  ice dynamics (T) or thermodynamics only (F)
-   amax        = 0.999      !  maximum ice concentration
-   cai         =  1.40e-3  !  atmospheric drag over sea ice
-   cao         =  1.00e-3  !  atmospheric drag over ocean
-   ln_nicep    = .false.   !  Ice points output for debug (yes or no)
-   ln_limdiahsb  = .false.    !  check the heat and salt budgets (T) or not (F)
-   ln_limdiaout  = .true.    !  output the heat and salt budgets (T) or not (F)
+   cn_icerst_out = "restart_ice"   !  suffix of ice restart name (output)
+   ln_limdyn     = .true.          !  ice dynamics (T) or thermodynamics only (F)
+   amax          = 0.999           !  maximum ice concentration
+   cai           = 1.40e-3         !  atmospheric drag over sea ice (clio)
+   cao           = 1.00e-3         !  atmospheric drag over ocean   (clio)
+   ln_nicep      = .false.         !  Ice points output for debug (yes or no)
+   ln_limdiahsb  = .false.          !  check the heat and salt budgets (T) or not (F)
+   ln_limdiaout  = .true.          !  output the heat and salt budgets (T) or not (F)
 /
 !-----------------------------------------------------------------------
@@ -147,15 +147,4 @@
 !&namicehsb       !  Heat and salt budgets 
 !!-----------------------------------------------------------------------
-!
 !/
-!-----------------------------------------------------------------------
-&namiceout     !   parameters for outputs
-!-----------------------------------------------------------------------
-!   noumef      =   43      !  number of fields
-!SF   add_diag_swi=    1      !  1 -> diagnose distribution in thickness space
-!SF                           !  0 -> only simple diagnostics
-!
-!           !         title of the field           !  name     !   units   !  save  ! multipl. ! additive !
-!           !                                      !           !           ! or not !  factor  !  factor  !
-/ 
 
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4649)
@@ -270,4 +270,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sum   ! vertical surface melt
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res   ! production (growth+melt) due to limupdate
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr   ! snow precipitation on ice
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice growth/melt                      [PSU/m2/s]
@@ -280,15 +281,23 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_res     !: residual salt flux due to correction of ice thickness [PSU/m2/s]
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_thd     !: ice-ocean heat flux from thermo processes (limthd_dh) 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dyn     !: ice-ocean heat flux from mecanical processes (limitd_me) 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_tot     !: total heat flux lost/gained by ice 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_spr     !: heat flux of the snow precipitation 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: residual heat flux due to correction of ice thickness
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bom     !: total heat flux causing bottom ice melt 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sum     !: total heat flux causing surface ice melt 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_opw     !: total heat flux causing open water ice formation
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dif     !: total heat flux causing Temp change in the ice 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err     !: heat flux error after heat diffusion 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_out     !: heat flux remaining at the end of thermo transformations 
+
+   ! heat flux associated with ice-atmosphere mass exchange
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_spr     !: heat flux of the snow precipitation 
+
+   ! heat flux associated with ice-ocean mass exchange
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_thd     !: ice-ocean heat flux from thermo processes (limthd_dh) 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dyn     !: ice-ocean heat flux from mecanical processes (limitd_me) 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: residual heat flux due to correction of ice thickness
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice     !: transmitted solar radiation under ice
@@ -416,5 +425,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_trp_es   ! transport of snw enthalpy (W/m2)
    !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_heat_dhc1    ! snw/ice heat content variation   [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_heat_dhc    ! snw/ice heat content variation   [W/m2] 
    !
    INTEGER , PUBLIC ::   jiindx, jjindx        !: indexes of the debugging point
@@ -450,15 +459,16 @@
 
       ii = ii + 1
-      ALLOCATE( sist     (jpi,jpj) , icethi (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) ,    &
-         &      wfx_bog(jpi,jpj)  , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     &
-         &      wfx_res(jpi,jpj)  , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) ,  qlead  (jpi,jpj) ,     &
-         &      fhtur    (jpi,jpj) , ftr_ice(jpi,jpj,jpl) ,      &
-         &      sfx_res  (jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,      &
-         &      sfx_bog  (jpi,jpj) , sfx_bom  (jpi,jpj) , sfx_sum  (jpi,jpj) ,  sfx_sni  (jpi,jpj) ,  sfx_opw  (jpi,jpj) ,   &
-         &      hfx_res  (jpi,jpj) , hfx_snw  (jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj), hfx_err_rem(jpi,jpj), &
-         &      hfx_in   (jpi,jpj) , hfx_out(jpi,jpj) , fhld(jpi,jpj) ,  &
-         &      hfx_tot  (jpi,jpj) , hfx_thd  (jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj),  STAT=ierr(ii) )
+      ALLOCATE( sist   (jpi,jpj) , icethi (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) ,    &
+         &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     &
+         &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,  qlead  (jpi,jpj) ,     &
+         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl) ,      &
+         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,      &
+         &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,   &
+         &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) , hfx_err_rem(jpi,jpj), &
+         &      hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld(jpi,jpj) ,  &
+         &      hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , hfx_opw(jpi,jpj) , &
+         &      hfx_thd(jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) ,  STAT=ierr(ii) )
 
       ii = ii + 1
@@ -526,5 +536,5 @@
       ALLOCATE( dv_dt_thd(jpi,jpj,jpl) ,     &
          &      izero    (jpi,jpj,jpl)  , diag_trp_vi(jpi,jpj) , diag_trp_vs(jpi,jpj), diag_trp_ei(jpi,jpj), diag_trp_es(jpi,jpj),     & 
-         &      diag_heat_dhc1(jpi,jpj) ,  STAT=ierr(ii) )
+         &      diag_heat_dhc(jpi,jpj) ,  STAT=ierr(ii) )
 
       ice_alloc = MAXVAL( ierr(:) )
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 4649)
@@ -7,4 +7,5 @@
    !!            3.0  ! 2004-06  (M. Vancoppenolle)   Energy Conservation 
    !!            4.0  ! 2011-02  (G. Madec)  add mpp considerations
+   !!             -   ! 2014-05  (C. Rousset) add lim_cons_hsm
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -14,4 +15,5 @@
    !!    lim_cons     :   checks whether energy, mass and salt are conserved 
    !!----------------------------------------------------------------------
+   USE phycst         ! physical constants
    USE par_ice        ! LIM-3 parameter
    USE ice            ! LIM-3 variables
@@ -28,4 +30,5 @@
    PUBLIC   lim_column_sum_energy
    PUBLIC   lim_cons_check
+   PUBLIC   lim_cons_hsm
 
    !!----------------------------------------------------------------------
@@ -151,4 +154,63 @@
    END SUBROUTINE lim_cons_check
 
+
+   SUBROUTINE lim_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
+      !!-------------------------------------------------------------------
+      !!               ***  ROUTINE lim_cons_hsm ***
+      !!
+      !! ** Purpose : Test the conservation of heat, salt and mass for each routine
+      !!
+      !! ** Method  :
+      !!---------------------------------------------------------------------
+      INTEGER         , INTENT(in)    :: icount      ! determine wether this is the beggining of the routine (0) or the end (1)
+      CHARACTER(len=*), INTENT(in)    :: cd_routine  ! name of the routine
+      REAL(wp)        , INTENT(inout) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
+      REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,   zft
+      REAL(wp)                        :: zvmin, zamin, zamax 
+
+      IF( icount == 0 ) THEN
+
+         zvi_b  = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) )
+         zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
+         zei_b  = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) )
+         zfw_b  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
+            &                   wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) ) * area(:,:) * tms(:,:) )
+         zfs_b  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
+            &                   sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) )
+         zft_b  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
+            &                 - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) ) * area(:,:) / unit_fac * tms(:,:) )
+
+      ELSEIF( icount == 1 ) THEN
+
+         zfs  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
+            &                 sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zfs_b
+         zfw  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
+            &                 wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) ) * area(:,:) * tms(:,:) ) - zfw_b
+         zft  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
+            &               - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zft_b
+ 
+         zvi  = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zvi_b ) * r1_rdtice - zfw 
+         zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic )
+         zei  =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zei_b * r1_rdtice + zft
+
+         zvmin = glob_min(v_i)
+         zamax = glob_max(SUM(a_i,dim=3))
+         zamin = glob_min(a_i)
+       
+         IF(lwp) THEN
+            IF ( ABS( zvi    ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (',cd_routine,') = ',(zvi * rday)
+            IF ( ABS( zsmv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday)
+            IF ( ABS( zei    ) >  1.    ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (',cd_routine,') = ',(zei)
+            IF ( zvmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin)
+            IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > amax+1.e-10 ) THEN
+                                          WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
+            ENDIF
+            IF ( zamin <  0.            ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
+         ENDIF
+
+      ENDIF
+
+   END SUBROUTINE lim_cons_hsm
+
 #else
    !!----------------------------------------------------------------------
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 4649)
@@ -35,6 +35,6 @@
    !!PUBLIC   lim_diahsb_rst   ! routine called by ice_init.F90
 
-   REAL(wp) ::   frc_sal, frc_vol   ! global forcing trends
-   REAL(wp) ::   bg_grme            ! global ice growth+melt trends
+   REAL(dp) ::   frc_sal, frc_vol   ! global forcing trends
+   REAL(dp) ::   bg_grme            ! global ice growth+melt trends
    REAL(wp) ::   epsi06 = 1.e-6_wp  ! small number
 
@@ -47,4 +47,5 @@
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
+
 CONTAINS
 
@@ -57,12 +58,14 @@
       !!---------------------------------------------------------------------------
       !!
-      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 
-      REAL(wp)   ::   zbg_vfx, zbg_vfx_bog, zbg_vfx_opw, zbg_vfx_sni, zbg_vfx_dyn, zbg_vfx_bom, zbg_vfx_sum, zbg_vfx_res 
-      REAL(wp)   ::   zbg_hfx_dhc1, zbg_hfx_spr, zbg_hfx_qsr, zbg_hfx_qns
-      REAL(wp)   ::   zbg_hfx_res, zbg_hfx_sub, zbg_hfx_dyn, zbg_hfx_thd, zbg_hfx_snw, zbg_hfx_tot, zbg_hfx_out, zbg_hfx_in   
-      REAL(wp)   ::   z_frc_vol, z_frc_sal, z_bg_grme 
-      REAL(wp)   ::   z1_area, zcoef                    
-      REAL(wp)   ::   zinda, zindb
+      REAL(dp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc
+      REAL(dp)   ::   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 
+      REAL(dp)   ::   zbg_vfx, zbg_vfx_bog, zbg_vfx_opw, zbg_vfx_sni, zbg_vfx_dyn
+      REAL(dp)   ::   zbg_vfx_bom, zbg_vfx_sum, zbg_vfx_res, zbg_vfx_spr, zbg_vfx_snw, zbg_vfx_sub  
+      REAL(dp)   ::   zbg_hfx_dhc, zbg_hfx_spr
+      REAL(dp)   ::   zbg_hfx_res, zbg_hfx_sub, zbg_hfx_dyn, zbg_hfx_thd, zbg_hfx_snw, zbg_hfx_out, zbg_hfx_in   
+      REAL(dp)   ::   zbg_hfx_sum, zbg_hfx_bom, zbg_hfx_bog, zbg_hfx_dif, zbg_hfx_opw
+      REAL(dp)   ::   z_frc_vol, z_frc_sal, z_bg_grme 
+      REAL(dp)   ::   z1_area                     !    -     -
+      REAL(dp)   ::   zinda, zindb
       !!---------------------------------------------------------------------------
       IF( nn_timing == 1 )   CALL timing_start('lim_diahsb')
@@ -70,7 +73,8 @@
       IF( numit == nstart ) CALL lim_diahsb_init 
 
-      z1_area = 1._wp / MAX( glob_sum( area(:,:) * tms(:,:) ), epsi06 )
-
-      zinda = MAX( 0._wp , SIGN( 1._wp , glob_sum( area(:,:) * tms(:,:) ) - epsi06 ) )
+      ! 1/area
+      z1_area = 1.d0 / MAX( glob_sum( area(:,:) * tms(:,:) ), epsi06 )
+
+      zinda = MAX( 0.d0 , SIGN( 1.d0 , glob_sum( area(:,:) * tms(:,:) ) - epsi06 ) )
       ! ----------------------- !
       ! 1 -  Content variations !
@@ -79,37 +83,40 @@
       zbg_svo = glob_sum( vt_s(:,:) * area(:,:) * tms(:,:) ) ! volume snow
       zbg_are = glob_sum( at_i(:,:) * area(:,:) * tms(:,:) ) ! area
-      zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 )      * area(:,:) * tms(:,:) )  ! mean salt content
+      zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )       ! mean salt content
       zbg_tem = glob_sum( ( tm_i(:,:) - rtt ) * vt_i(:,:) * area(:,:) * tms(:,:) )  ! mean temp content
 
-      zcoef = zinda * z1_area * r1_rau0 * rday
+      !zbg_ihc = glob_sum( et_i(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content
+      !zbg_shc = glob_sum( et_s(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_svo,epsi06 ) ! snow heat content
+
       ! Volume
-      zbg_vfx     = glob_sum(     emp(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_vfx_bog = glob_sum( wfx_bog(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_vfx_opw = glob_sum( wfx_opw(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_vfx_sni = glob_sum( wfx_sni(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_vfx_dyn = glob_sum( wfx_dyn(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_vfx_bom = glob_sum( wfx_bom(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_vfx_sum = glob_sum( wfx_sum(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_vfx_res = glob_sum( wfx_res(:,:) * area(:,:) * tms(:,:) ) * zcoef
+      zbg_vfx     = zinda * glob_sum(      emp(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_vfx_bog = zinda * glob_sum( wfx_bog(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_vfx_opw = zinda * glob_sum( wfx_opw(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_vfx_sni = zinda * glob_sum( wfx_sni(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_vfx_dyn = zinda * glob_sum( wfx_dyn(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_vfx_bom = zinda * glob_sum( wfx_bom(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_vfx_sum = zinda * glob_sum( wfx_sum(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_vfx_res = zinda * glob_sum( wfx_res(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_vfx_spr = zinda * glob_sum( wfx_spr(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_vfx_snw = zinda * glob_sum( wfx_snw(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_vfx_sub = zinda * glob_sum( wfx_sub(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
 
       ! Salt
-      zbg_sfx     = glob_sum(     sfx(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_sfx_bri = glob_sum( sfx_bri(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_sfx_res = glob_sum( sfx_res(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_sfx_dyn = glob_sum( sfx_dyn(:,:) * area(:,:) * tms(:,:) ) * zcoef
-
-      zbg_sfx_bog = glob_sum( sfx_bog(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_sfx_opw = glob_sum( sfx_opw(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_sfx_sni = glob_sum( sfx_sni(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_sfx_bom = glob_sum( sfx_bom(:,:) * area(:,:) * tms(:,:) ) * zcoef
-      zbg_sfx_sum = glob_sum( sfx_sum(:,:) * area(:,:) * tms(:,:) ) * zcoef
+      zbg_sfx     = zinda * glob_sum(     sfx(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_sfx_bri = zinda * glob_sum( sfx_bri(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_sfx_res = zinda * glob_sum( sfx_res(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_sfx_dyn = zinda * glob_sum( sfx_dyn(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+
+      zbg_sfx_bog = zinda * glob_sum( sfx_bog(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_sfx_opw = zinda * glob_sum( sfx_opw(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_sfx_sni = zinda * glob_sum( sfx_sni(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_sfx_bom = zinda * glob_sum( sfx_bom(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
+      zbg_sfx_sum = zinda * glob_sum( sfx_sum(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday
 
       ! Heat budget
       zbg_ihc      = glob_sum( et_i(:,:) * 1.e-20 ) ! ice heat content  [1.e-20 J]
       zbg_shc      = glob_sum( et_s(:,:) * 1.e-20 ) ! snow heat content [1.e-20 J]
-      zbg_hfx_dhc1 = glob_sum( diag_heat_dhc1(:,:) * area(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_dhc  = glob_sum( diag_heat_dhc(:,:) * area(:,:) * tms(:,:) ) ! [in W]
       zbg_hfx_spr  = glob_sum( hfx_spr(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_qsr  = glob_sum( qsr(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_qns  = glob_sum( qns(:,:) * area(:,:) * tms(:,:) ) ! [in W]
 
       zbg_hfx_thd  = glob_sum( hfx_thd(:,:) * area(:,:) * tms(:,:) ) ! [in W]
@@ -118,5 +125,9 @@
       zbg_hfx_sub  = glob_sum( hfx_sub(:,:) * area(:,:) * tms(:,:) ) ! [in W]
       zbg_hfx_snw  = glob_sum( hfx_snw(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_tot  = glob_sum( hfx_tot(:,:) * area(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_sum  = glob_sum( hfx_sum(:,:) * area(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_bom  = glob_sum( hfx_bom(:,:) * area(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_bog  = glob_sum( hfx_bog(:,:) * area(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_dif  = glob_sum( hfx_dif(:,:) * area(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_opw  = glob_sum( hfx_opw(:,:) * area(:,:) * tms(:,:) ) ! [in W]
       zbg_hfx_out  = glob_sum( hfx_out(:,:) * area(:,:) * tms(:,:) ) ! [in W]
       zbg_hfx_in   = glob_sum(  hfx_in(:,:) * area(:,:) * tms(:,:) ) ! [in W]
@@ -127,15 +138,19 @@
       z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * area(:,:) * tms(:,:) ) ! volume fluxes
       z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * area(:,:) * tms(:,:) ) ! salt fluxes
-      z_bg_grme = glob_sum( ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + &
-                          &   wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) ) / rhoic * area(:,:) * tms(:,:) ) ! 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
-            
+      z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + &
+                          &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * area(:,:) * tms(:,:) ) ! 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
+      
+      ! difference
+      !frc_vol = zbg_ivo - frc_vol
+      !frc_sal = zbg_sal - frc_sal
+      
       ! ----------------------- !
       ! 3 - Diagnostics writing !
       ! ----------------------- !
-      zindb = MAX( 0._wp , SIGN( 1._wp , zbg_ivo - epsi06 ) )
+      zindb = MAX( 0.d0 , SIGN( 1.d0 , zbg_ivo - epsi06 ) )
       !
       CALL iom_put( 'ibgvoltot' , zbg_ivo * rhoic * r1_rau0 * 1.e-9        )   ! ice volume (km3 equivalent liquid)         
@@ -156,4 +171,7 @@
       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)       
@@ -167,8 +185,6 @@
       CALL iom_put( 'ibgsfxsum' , zbg_sfx_sum                              )   ! salt flux surface melt      -
 
-      CALL iom_put( 'ibghfxdhc1', zbg_hfx_dhc1                            )   ! Heat content variation in snow and ice [W]
+      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( 'ibghfxqsr' , zbg_hfx_qsr                              )   !     solar fluxes used by snw/ice [W]
-      CALL iom_put( 'ibghfxqns' , zbg_hfx_qns                              )   ! non solar fluxes used by snw/ice [W]
 
       CALL iom_put( 'ibghfxres' , zbg_hfx_res                              )   ! 
@@ -177,11 +193,15 @@
       CALL iom_put( 'ibghfxthd' , zbg_hfx_thd                              )   ! 
       CALL iom_put( 'ibghfxsnw' , zbg_hfx_snw                              )   ! 
-      CALL iom_put( 'ibghfxtot' , zbg_hfx_tot                              )   ! 
+      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( '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)   
-      CALL iom_put( 'ibgvolgrm' , bg_grme * rhoic * r1_rau0 * 1.e-9        )   ! vol growth + melt (km3 equivalent liquid)         
+      CALL iom_put( 'ibgvolgrm' , bg_grme * r1_rau0 * 1.e-9                )   ! vol growth + melt (km3 equivalent liquid)         
 
       !
@@ -189,5 +209,5 @@
       !
       IF( nn_timing == 1 )   CALL timing_stop('lim_diahsb')
-      !
+!
    END SUBROUTINE lim_diahsb
 
@@ -223,4 +243,7 @@
       ! 2 - initial conservation variables !
       ! ---------------------------------- !
+      !frc_vol = 0.d0                                           ! volume       trend due to forcing
+      !frc_sal = 0.d0                                           ! salt content   -    -   -    -         
+      !bg_grme = 0.d0                                           ! ice growth + melt volume trend
       !
       CALL lim_diahsb_rst( nstart, 'READ' )  !* read or initialize all required files
@@ -256,7 +279,7 @@
            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                                        
+           frc_vol  = 0.d0                                           
+           frc_sal  = 0.d0                                                  
+           bg_grme  = 0.d0                                        
        ENDIF	
 
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 4649)
@@ -30,4 +30,5 @@
    USE lib_fortran      ! glob_sum
    USE timing          ! Timing
+   USE limcons        ! conservation tests
 
    IMPLICIT NONE
@@ -66,6 +67,6 @@
       REAL(wp), POINTER, DIMENSION(:)   ::   zmsk           ! i-averaged of tmask
       REAL(wp), POINTER, DIMENSION(:,:) ::   zu_io, zv_io   ! ice-ocean velocity
-      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
-      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
+      !
+      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
      !!---------------------------------------------------------------------
 
@@ -75,19 +76,11 @@
       CALL wrk_alloc( jpj, zind, zmsk )
 
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF (ln_limdiahsb) THEN
-         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) )
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
-
       IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
 
       IF( ln_limdyn ) THEN
          !
+         ! conservation test
+         IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
          old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)
          old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)
@@ -171,4 +164,7 @@
             END DO
          END DO
+         !
+         ! conservation test
+         IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
          !
       ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
@@ -224,28 +220,4 @@
       ENDIF
       !
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF (ln_limdiahsb) THEN
-         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
-         zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b
- 
-         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - ( zchk_fw / rhoic )
-         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
-
-         zchk_vmin = glob_min(v_i)
-         zchk_amax = glob_max(SUM(a_i,dim=3))
-         zchk_amin = glob_min(a_i)
-
-         IF(lwp) THEN
-            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limdyn) = ',(zchk_v_i * rday)
-            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limdyn) = ',(zchk_smv * rday)
-            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limdyn) = ',(zchk_vmin * 1.e-3)
-            !IF ( zchk_amax >  amax+1.e-10   ) WRITE(numout,*) 'violation a_i>amax            (limdyn) = ',zchk_amax
-            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limdyn) = ',zchk_amin
-         ENDIF
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
-
       CALL wrk_dealloc( jpi, jpj, zu_io, zv_io )
       CALL wrk_dealloc( jpj, zind, zmsk )
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4649)
@@ -26,5 +26,4 @@
    USE dom_ice          ! sea-ice domain
    USE in_out_manager   ! I/O manager
-   USE lbclnk           ! lateral boundary condition - MPP exchanges
    USE lib_mpp          ! MPP library
    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
@@ -390,33 +389,4 @@
       tn_ice (:,:,:) = t_su (:,:,:)
 
-      !--------------------------------------------------------------------
-      ! 4) Lateral boundary conditions                                    | 
-      !--------------------------------------------------------------------
-      DO jl = 1, jpl
-
-         CALL lbc_lnk( a_i(:,:,jl)  , 'T', 1. )
-         CALL lbc_lnk( v_i(:,:,jl)  , 'T', 1. )
-         CALL lbc_lnk( v_s(:,:,jl)  , 'T', 1. )
-         CALL lbc_lnk( smv_i(:,:,jl), 'T', 1. )
-         CALL lbc_lnk( oa_i(:,:,jl) , 'T', 1. )
-
-         CALL lbc_lnk( ht_i(:,:,jl) , 'T', 1. )
-         CALL lbc_lnk( ht_s(:,:,jl) , 'T', 1. )
-         CALL lbc_lnk( sm_i(:,:,jl) , 'T', 1. )
-         CALL lbc_lnk( o_i(:,:,jl)  , 'T', 1. )
-         CALL lbc_lnk( t_su(:,:,jl) , 'T', 1. )
-         CALL lbc_lnk( tn_ice(:,:,jl) , 'T', 1. )
-         DO jk = 1, nlay_s
-            CALL lbc_lnk(t_s(:,:,jk,jl), 'T', 1. )
-            CALL lbc_lnk(e_s(:,:,jk,jl), 'T', 1. )
-         END DO
-         DO jk = 1, nlay_i
-            CALL lbc_lnk(t_i(:,:,jk,jl), 'T', 1. )
-            CALL lbc_lnk(e_i(:,:,jk,jl), 'T', 1. )
-         END DO
-         !
-!         a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl)
-      END DO
-
       ELSE 
          ! if ln_limini=false
@@ -449,9 +419,7 @@
          at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
       END DO
-      CALL lbc_lnk( at_i , 'T', 1. )
-!      at_i(:,:) = tms(:,:) * at_i(:,:)
       !
       !--------------------------------------------------------------------
-      ! 5) Global ice variables for output diagnostics                    | 
+      ! 4) Global ice variables for output diagnostics                    | 
       !--------------------------------------------------------------------
       u_ice (:,:)     = 0._wp
@@ -462,5 +430,5 @@
 
       !--------------------------------------------------------------------
-      ! 6) Moments for advection
+      ! 5) Moments for advection
       !--------------------------------------------------------------------
 
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4649)
@@ -22,5 +22,4 @@
    USE limthd_lac       ! LIM
    USE limvar           ! LIM
-   USE limcons          ! LIM
    USE in_out_manager   ! I/O manager
    USE lbclnk           ! lateral boundary condition - MPP exchanges
@@ -33,4 +32,5 @@
    USE limdiahsb
    USE timing          ! Timing
+   USE limcons        ! conservation tests
 
    IMPLICIT NONE
@@ -143,6 +143,6 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2)
       REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final  !  ice volume summed over categories
-      REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b 
-      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
+      !
+      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!-----------------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limitd_me')
@@ -158,16 +158,7 @@
 
       IF( ln_limdyn ) THEN          !   Start ridging and rafting   !
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF (ln_limdiahsb) THEN
-         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) )
-         zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) )
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
+
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       !-----------------------------------------------------------------------------!
@@ -355,5 +346,5 @@
             ! 5) Heat, salt and freshwater fluxes
             !-----------------------------------------------------------------------------!
-            wfx_snw(ji,jj) = wfx_snw(ji,jj) - msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean
+            wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean
             hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * unit_fac / area(ji,jj) * r1_rdtice  ! heat sink for ocean (<0, W.m-2)
 
@@ -425,30 +416,6 @@
       ENDIF
 
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF (ln_limdiahsb) THEN
-         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
-         zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b
-         zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b
- 
-         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 
-         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
-         zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft
-
-         zchk_vmin = glob_min(v_i)
-         zchk_amax = glob_max(SUM(a_i,dim=3))
-         zchk_amin = glob_min(a_i)
-       
-         IF(lwp) THEN
-            IF ( ABS( zchk_v_i   ) >  1.e-2 ) WRITE(numout,*) 'violation volume [kg/day]     (limitd_me) = ',(zchk_v_i * rday)
-            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday)
-            IF ( ABS( zchk_e_i   ) >  1.e-4 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limitd_me) = ',(zchk_e_i)
-            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(zchk_vmin * 1.e-3)
-            IF ( zchk_amax >  kamax+epsi10  ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',zchk_amax
-            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_me) = ',zchk_amin
-         ENDIF
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       ENDIF  ! ln_limdyn=.true.
@@ -1108,6 +1075,5 @@
             
             sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice
-            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) + vsw (ji,jj) * rhoic * r1_rdtice   ! gurvan: increase in ice volume du to seawater frozen in voids             
-            ! MV HC 2014 this previous line seems ok, i'm not sure at this moment of the sign convention
+            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice   ! gurvan: increase in ice volume du to seawater frozen in voids             
 
             !------------------------------------            
@@ -1440,9 +1406,10 @@
       CALL wrk_alloc( jpi, jpj, zmask )
 
+      ! to be sure that at_i is the sum of a_i(jl)
+      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
+
       DO jl = 1, jpl
-
          !-----------------------------------------------------------------
          ! Count categories to be zapped.
-         ! Abort model in case of negative area.
          !-----------------------------------------------------------------
          icells = 0
@@ -1450,10 +1417,7 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-!               IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <  0._wp   ) .OR.   &
-!                  & ( a_i(ji,jj,jl) >  0._wp   .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   &
-!                  & ( v_i(ji,jj,jl) == 0._wp   .AND. a_i(ji,jj,jl) >  0._wp   ) .OR.   &
-!                  & ( v_i(ji,jj,jl) >  0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp
-               IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   &
-                  & ( v_i(ji,jj,jl) >= 0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp
+               IF( a_i(ji,jj,jl) <= epsi10 .OR. v_i(ji,jj,jl) <= epsi10 .OR. at_i(ji,jj) <= epsi10 ) THEN
+                  zmask(ji,jj) = 1._wp
+               ENDIF
             END DO
          END DO
@@ -1470,4 +1434,5 @@
                   zei  = e_i(ji,jj,jk,jl)
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) )
+                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) + rtt * zmask(ji,jj)
                   ! update exchanges with ocean
                   hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
@@ -1497,13 +1462,5 @@
                ! zap ice and snow volume, add water and salt to ocean
                !-----------------------------------------------------------------
-
-               !           xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt
-               !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj)                  )   &
-               !                                            * rhosn * v_s(ji,jj,jl) * r1_rdtice
-               !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) )   & 
-               !                                            * rhoic * v_i(ji,jj,jl) * r1_rdtice
-               !           sfx (i,j)      = sfx (i,j)      + xtmp
-
-               ato_i(ji,jj)    = a_i  (ji,jj,jl) *       zmask(ji,jj)   + ato_i(ji,jj)
+               ato_i(ji,jj)    = a_i  (ji,jj,jl) *           zmask(ji,jj)   + ato_i(ji,jj)
                a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
                v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )
@@ -1520,11 +1477,13 @@
                ! update exchanges with ocean
                sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
-               wfx_res(ji,jj)  = wfx_res(ji,jj) + ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
-               wfx_snw(ji,jj)  = wfx_snw(ji,jj) + ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
+               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
+               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
                hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
             END DO
          END DO
-         !
-      END DO                 ! jl 
+      END DO ! jl 
+
+      ! to be sure that at_i is the sum of a_i(jl)
+      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
       !
       CALL wrk_dealloc( jpi, jpj, zmask )
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 4649)
@@ -35,4 +35,5 @@
    USE lib_fortran      ! to use key_nosignedzero
    USE timing          ! Timing
+   USE limcons        ! conservation tests
 
    IMPLICIT NONE
@@ -66,21 +67,11 @@
       !
       INTEGER ::   ji,jj, jk, jl, ja, jm, jbnd1, jbnd2   ! ice types    dummy loop index         
-      REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset)
-      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
+      !
+      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limitd_th')
 
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF (ln_limdiahsb) THEN
-         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) )
-         zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) )
-       ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       IF( kt == nit000 .AND. lwp ) THEN
@@ -107,5 +98,4 @@
       !  3) Add frazil ice growing in leads.
       !------------------------------------------------------------------------------|
-
       CALL lim_thd_lac
       CALL lim_var_glo2eqv    ! only for info
@@ -143,30 +133,6 @@
       ENDIF
       !
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF( ln_limdiahsb ) THEN
-         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
-         zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b
-         zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b
- 
-         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 
-         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
-         zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft
-
-         zchk_vmin = glob_min(v_i)
-         zchk_amax = glob_max(SUM(a_i,dim=3))
-         zchk_amin = glob_min(a_i)
-
-         IF(lwp) THEN
-            IF ( ABS( zchk_v_i   ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (limitd_th) = ',(zchk_v_i * rday)
-            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * rday)
-            IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limitd_th) = ',(zchk_e_i)
-            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_th) = ',(zchk_vmin * 1.e-3)
-            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limitd_th) = ',zchk_amax
-            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_th) = ',zchk_amin
-         ENDIF
-       ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !
      IF( nn_timing == 1 )  CALL timing_stop('limitd_th')
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4649)
@@ -42,4 +42,5 @@
    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE traqsr           ! clem: add penetration of solar flux into the calculation of heat budget
+   USE iom
 
    IMPLICIT NONE
@@ -104,12 +105,19 @@
       INTEGER  ::   ji, jj, jl, jk           ! dummy loop indices
       REAL(wp) ::   zinda, zemp      ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace
-      REAL(wp) ::   ztmelts         ! clem 2014: for HC diags
-
       REAL(wp) ::   zf_mass         ! Heat flux associated with mass exchange ice->ocean (W.m-2)
       REAL(wp) ::   zfcm1           ! New solar flux received by the ocean
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace
       !!---------------------------------------------------------------------
       
       IF( lk_cpl )   CALL wrk_alloc( jpi, jpj, jpl, zalb, zalbp )
+
+      ! make calls for heat fluxes before it is modified
+      CALL iom_put( "qsr_oce" , qsr(:,:) * pfrld(:,:) )   !     solar flux at ocean surface
+      CALL iom_put( "qns_oce" , qns(:,:) * pfrld(:,:) )   ! non-solar flux at ocean surface
+      CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  !     solar flux at ice surface
+      CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  ! non-solar flux at ice surface
+      CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  !     solar flux transmitted thru ice
+      CALL iom_put( "qt_oce"  , ( qsr(:,:) + qns(:,:) ) * pfrld(:,:) )  
+      CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * old_a_i(:,:,:), dim=3 ) )
 
       ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
@@ -165,5 +173,5 @@
             IF( lk_cpl ) THEN 
                zemp = - emp_tot(ji,jj) + emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !
-                  &   - wfx_snw(ji,jj)
+                  &   + wfx_snw(ji,jj)
             ELSE
                zemp =   emp(ji,jj)     *           pfrld(ji,jj)            &   ! evaporation over oceanic fraction
@@ -176,6 +184,6 @@
 
             ! mass flux at the ocean/ice interface
-            fmmflx(ji,jj) = wfx_ice(ji,jj) * rdt_ice                   ! F/M mass flux save at least for biogeochemical model
-            emp(ji,jj)    = zemp + wfx_ice(ji,jj) + wfx_snw(ji,jj)     ! mass flux + F/M mass flux (always ice/ocean mass exchange)
+            fmmflx(ji,jj) = - wfx_ice(ji,jj) * rdt_ice                   ! F/M mass flux save at least for biogeochemical model
+            emp(ji,jj)    = zemp - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_sub(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange)
             
          END DO
@@ -213,16 +221,4 @@
       ENDIF
 
-      ! -------------------------------------------------
-      ! C. Rousset Begin Diagnostics for heat in W/m2
-      ! -------------------------------------------------
-      DO jj = 1, jpj
-         DO ji = 1, jpi            
-            diag_heat_dhc1(ji,jj) = ( SUM( d_e_i_trp(ji,jj,1:nlay_i,:) + d_e_i_thd(ji,jj,1:nlay_i,:) ) +  & 
-               &                      SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) ) ) * unit_fac * r1_rdtice / area(ji,jj)   
-         END DO
-      END DO
-      ! -------------------------------------------------
-      ! C. Rousset End Diagnostics
-      ! -------------------------------------------------
 
       IF(ln_ctl) THEN
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4649)
@@ -44,4 +44,5 @@
    USE timing         ! Timing
    USE cpl_oasis3, ONLY : lk_cpl
+   USE limcons        ! conservation tests
 
    IMPLICIT NONE
@@ -86,13 +87,12 @@
       INTEGER  :: nbpb             ! nb of icy pts for thermo. cal.
       INTEGER  :: ii, ij           ! temporary dummy loop index
-      REAL(wp) :: zfric_umin = 5e-03_wp    ! lower bound for the friction velocity
-      REAL(wp) :: zfric_umax = 2e-02_wp    ! upper bound for the friction velocity
-      REAL(wp) :: zinda, zindb, zfric_u     ! local scalar
-      REAL(wp) :: zareamin  !    -         -
-      REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b 
-      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
-      REAL(wp) :: zqld, zqfr
+      REAL(wp) :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04)
+      REAL(wp) :: zch        = 0.0057_wp    ! heat transfer coefficient
+      REAL(wp) :: zinda, zindb, zareamin 
+      REAL(wp) :: zfric_u, zqld, zqfr
       REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx, zqfx
       REAL(wp)                        :: zhfx_err, ztest
+      !
+      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!-------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limthd')
@@ -103,16 +103,6 @@
       zdq(:) = 0._wp ; zq_ini(:) = 0._wp ; zhfx(:) = 0._wp ; zqfx(:) = 0._wp      
 
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF (ln_limdiahsb) THEN
-         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) )
-         zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) )
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       !------------------------------------------------------------------------------!
@@ -158,5 +148,5 @@
 !CDIR NOVERRCHK
          DO ji = 1, jpi
-            zinda          = tms(ji,jj) * ( 1.0 - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice
+            zinda          = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice
             !
             !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
@@ -165,24 +155,14 @@
             !           !  net downward heat flux from the ice to the ocean, expressed as a function of ocean 
             !           !  temperature and turbulent mixing (McPhee, 1992)
-
-            ! friction velocity
-            zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) 
-
-            !-- Energy from the turbulent oceanic heat flux. here the drag will depend on ice thickness and type (0.006)
-            fhtur(ji,jj)  = zinda * rau0 * rcp * 0.006  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 
-            ! clem: why not the following? 
-            !fhtur(ji,jj)  = zinda * rau0 * rcp * 0.006  * SQRT( ust2s(ji,jj) ) * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) )
-
-            !-- Energy received in the lead, zqld is defined everywhere (J.m-2)
-            !   It includes turbulent ocean heat flux (only in the leads, the rest is used for bottom melting)
-            zqld =  tms(ji,jj) * rdt_ice *                               &
-               &  ( pfrld(ji,jj)        * ( qsr(ji,jj) * oatte(ji,jj)           &   ! solar heat + clem modif
-               &                          + qns(ji,jj)                          &   ! non solar heat
-               &                          + fhtur(ji,jj) )                      &   ! turbulent ice-ocean heat (0 if no ice)
+            !
+            ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- !
+            zqld =  tms(ji,jj) * rdt_ice *                                       &
+               &  ( pfrld(ji,jj)         * ( qsr(ji,jj) * oatte(ji,jj)           &   ! solar heat + clem modif
+               &                           + qns(ji,jj) )                        &   ! non solar heat
                ! latent heat of precip (note that precip is included in qns but not in qns_ice)
                &    + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )  &
                &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) )
 
-            !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2)
+            !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
             zqfr = tms(ji,jj) * rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
 
@@ -192,10 +172,17 @@
             ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting 
             IF( at_i(ji,jj) > epsi10 .AND. zqld > 0._wp ) THEN
-               fhld (ji,jj) = zqld * r1_rdtice / at_i(ji,jj) ! divided by a_i since this is (re)multiplied by a_i in limthd_dh.F90
+               fhld (ji,jj) = zqld * r1_rdtice / at_i(ji,jj) ! divided by at_i since this is (re)multiplied by a_i in limthd_dh.F90
                qlead(ji,jj) = 0._wp
             ENDIF
             !
-            IF( qlead(ji,jj) == 0._wp )  zqld = 0._wp ; zqfr = 0._wp
-            !
+            !-- Energy from the turbulent oceanic heat flux --- !
+            !clem zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )
+            zfric_u      = MAX( SQRT( ust2s(ji,jj) ), zfric_umin ) 
+            fhtur(ji,jj) = MAX( 0._wp, zinda * rau0 * rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2 
+            ! upper bound for fhtur: we do not want SST to drop below Tfreeze. 
+            ! So we say that the heat retrieved from the ocean (fhtur+fhld) must be < to the heat necessary to reach Tfreeze (zqfr)   
+            ! This is not a clean budget, so that should be corrected at some point
+            fhtur(ji,jj) = zinda * MIN( fhtur(ji,jj), - fhld(ji,jj) - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
+
             ! -----------------------------------------
             ! Net heat flux on top of ice-ocean [W.m-2]
@@ -317,4 +304,6 @@
             CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum         , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni         , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res         , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr         , jpi, jpj, npb(1:nbpb) )
 
             CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog         , jpi, jpj, npb(1:nbpb) )
@@ -323,4 +312,5 @@
             CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni         , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri         , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res         , jpi, jpj, npb(1:nbpb) )
 
             CALL tab_2d_1d( nbpb, iatte_1d   (1:nbpb), iatte           , jpi, jpj, npb(1:nbpb) ) 
@@ -329,5 +319,9 @@
             CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd         , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr         , jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, hfx_tot_1d (1:nbpb), hfx_tot         , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum         , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom         , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog         , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif         , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw         , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw         , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub         , jpi, jpj, npb(1:nbpb) )
@@ -365,5 +359,4 @@
                ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
                hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )
-
             END DO
 
@@ -379,6 +372,6 @@
             CALL lim_thd_dh( 1, nbpb, jl )    
 
-            ! --- Ice/Snow enthalpy remapping --- !
-            CALL lim_thd_ent( 1, nbpb, jl ) 
+            ! --- Ice enthalpy remapping --- !
+            CALL lim_thd_ent( 1, nbpb, jl, q_i_b(1:nbpb,:) ) 
             !                                
             ! --- diag error on heat remapping - PART 2 --- !
@@ -391,9 +384,13 @@
 
             !---------------------------------!
-            ! Ice salinity                    !
+            ! --- Ice salinity --- !
             !---------------------------------!
             CALL lim_thd_sal( 1, nbpb )    
 
-            !           CALL lim_thd_enmelt(1,nbpb)   ! computes sea ice energy of melting
+            !---------------------------------!
+            ! --- temperature update --- !
+            !---------------------------------!
+            CALL lim_thd_temp( 1, nbpb )
+
             !--------------------------------
             ! 4.4) Move 1D to 2D vectors
@@ -424,4 +421,6 @@
                CALL tab_1d_2d( nbpb, wfx_sum       , npb, wfx_sum_1d(1:nbpb)   , jpi, jpj )
                CALL tab_1d_2d( nbpb, wfx_sni       , npb, wfx_sni_1d(1:nbpb)   , jpi, jpj )
+               CALL tab_1d_2d( nbpb, wfx_res       , npb, wfx_res_1d(1:nbpb)   , jpi, jpj )
+               CALL tab_1d_2d( nbpb, wfx_spr       , npb, wfx_spr_1d(1:nbpb)   , jpi, jpj )
 
                CALL tab_1d_2d( nbpb, sfx_bog       , npb, sfx_bog_1d(1:nbpb)   , jpi, jpj )
@@ -429,4 +428,5 @@
                CALL tab_1d_2d( nbpb, sfx_sum       , npb, sfx_sum_1d(1:nbpb)   , jpi, jpj )
                CALL tab_1d_2d( nbpb, sfx_sni       , npb, sfx_sni_1d(1:nbpb)   , jpi, jpj )
+               CALL tab_1d_2d( nbpb, sfx_res       , npb, sfx_res_1d(1:nbpb)   , jpi, jpj )
             !
             IF( num_sal == 2 ) THEN
@@ -436,5 +436,9 @@
               CALL tab_1d_2d( nbpb, hfx_thd       , npb, hfx_thd_1d(1:nbpb)   , jpi, jpj )
               CALL tab_1d_2d( nbpb, hfx_spr       , npb, hfx_spr_1d(1:nbpb)   , jpi, jpj )
-              CALL tab_1d_2d( nbpb, hfx_tot       , npb, hfx_tot_1d(1:nbpb)   , jpi, jpj )
+              CALL tab_1d_2d( nbpb, hfx_sum       , npb, hfx_sum_1d(1:nbpb)   , jpi, jpj )
+              CALL tab_1d_2d( nbpb, hfx_bom       , npb, hfx_bom_1d(1:nbpb)   , jpi, jpj )
+              CALL tab_1d_2d( nbpb, hfx_bog       , npb, hfx_bog_1d(1:nbpb)   , jpi, jpj )
+              CALL tab_1d_2d( nbpb, hfx_dif       , npb, hfx_dif_1d(1:nbpb)   , jpi, jpj )
+              CALL tab_1d_2d( nbpb, hfx_opw       , npb, hfx_opw_1d(1:nbpb)   , jpi, jpj )
               CALL tab_1d_2d( nbpb, hfx_snw       , npb, hfx_snw_1d(1:nbpb)   , jpi, jpj )
               CALL tab_1d_2d( nbpb, hfx_sub       , npb, hfx_sub_1d(1:nbpb)   , jpi, jpj )
@@ -521,30 +525,6 @@
       ENDIF
       !
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF (ln_limdiahsb) THEN
-         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
-         zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b
-         zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b
- 
-         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 
-         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
-         zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft
-
-         zchk_vmin = glob_min(v_i)
-         zchk_amax = glob_max(SUM(a_i,dim=3))
-         zchk_amin = glob_min(a_i)
-       
-         IF(lwp) THEN
-            IF ( ABS( zchk_v_i   ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (limthd) = ',(zchk_v_i * rday)
-            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * rday)
-            IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limthd) = ',(zchk_e_i)
-            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limthd) = ',(zchk_vmin * 1.e-3)
-            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limthd) = ',zchk_amax
-            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limthd) = ',zchk_amin
-         ENDIF
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !
       CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx, zqfx )
@@ -558,5 +538,5 @@
       !!                   ***  ROUTINE lim_thd_enmelt *** 
       !!                 
-      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3)
+      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
       !!
       !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
@@ -565,13 +545,14 @@
       !!
       INTEGER  ::   ji, jk   ! dummy loop indices
-      REAL(wp) ::   ztmelts  ! local scalar 
+      REAL(wp) ::   ztmelts, zindb  ! local scalar 
       !!-------------------------------------------------------------------
       !
       DO jk = 1, nlay_i             ! Sea ice energy of melting
          DO ji = kideb, kiut
-            ztmelts      =  - tmut  * s_i_b(ji,jk) + rtt 
-            q_i_b(ji,jk) =    rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) )                                 &
-               &                      + lfus * ( 1.0 - (ztmelts-rtt) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) )   &
-               &                      - rcp  * ( ztmelts-rtt  )  ) 
+            ztmelts      = - tmut  * s_i_b(ji,jk) + rtt 
+            zindb        = MAX( 0._wp , SIGN( 1._wp , -(t_i_b(ji,jk) - rtt) - epsi10 ) )
+            q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) )                                             &
+               &                   + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) )   &
+               &                   - rcp  *                 ( ztmelts-rtt )  ) 
          END DO
       END DO
@@ -584,4 +565,35 @@
    END SUBROUTINE lim_thd_enmelt
 
+   SUBROUTINE lim_thd_temp( kideb, kiut )
+      !!-----------------------------------------------------------------------
+      !!                   ***  ROUTINE lim_thd_temp *** 
+      !!                 
+      !! ** Purpose :   Computes sea ice temperature (Kelvin) from enthalpy
+      !!
+      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
+      !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
+      !!
+      INTEGER  ::   ji, jk   ! dummy loop indices
+      REAL(wp) ::   ztmelts, zswitch, zaaa, zbbb, zccc, zdiscrim  ! local scalar 
+      !!-------------------------------------------------------------------
+      ! Recover ice temperature
+      DO jk = 1, nlay_i
+         DO ji = kideb, kiut
+            ztmelts       =  -tmut * s_i_b(ji,jk) + rtt
+            ! Conversion q(S,T) -> T (second order equation)
+            zaaa          =  cpic
+            zbbb          =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus
+            zccc          =  lfus * ( ztmelts - rtt )
+            zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) )
+            t_i_b(ji,jk)  =  rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )
+            
+            ! mask temperature
+            zswitch      =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) ) 
+            t_i_b(ji,jk) =  zswitch * t_i_b(ji,jk) + ( 1._wp - zswitch ) * rtt
+         END DO 
+      END DO 
+
+   END SUBROUTINE lim_thd_temp
 
    SUBROUTINE lim_thd_init
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4649)
@@ -72,5 +72,4 @@
       INTEGER  ::   ji , jk        ! dummy loop indices
       INTEGER  ::   ii, ij         ! 2D corresponding indices to ji
-      INTEGER  ::   i_ice_switch   ! ice thickness above a certain treshold or not
       INTEGER  ::   iter
 
@@ -114,5 +113,5 @@
 
       ! mass and salt flux (clem)
-      REAL(wp) :: zdvres, zswitch_sal
+      REAL(wp) :: zdvres, zswitch_sal, zswitch
 
       ! Heat conservation 
@@ -180,5 +179,5 @@
          zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_b(ji) - ztmelts ) )
          zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
-      END DO ! ji
+      END DO
 
       !
@@ -192,5 +191,5 @@
             hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_b(ji,1) * ht_s_b(ji) * a_i_b(ji) * r1_rdtice
             ! Contribution to mass flux
-            wfx_snw_1d(ji) =  wfx_snw_1d(ji) - rhosn * ht_s_b(ji) * a_i_b(ji) * r1_rdtice
+            wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_b(ji) * a_i_b(ji) * r1_rdtice
             ! updates
             ht_s_b(ji)   = 0._wp
@@ -247,8 +246,10 @@
          zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
          ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K)
-         zqprec    (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )   
+         zqprec   (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )   
          IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
          ! heat flux from snow precip (>0, W.m-2)
          hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
+         ! mass flux, <0
+         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_b(ji) * zdh_s_pre(ji) * r1_rdtice
          ! update thickness
          ht_s_b    (ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_pre(ji) )
@@ -258,13 +259,12 @@
          !---------------------
          ! thickness change
+         IF( zdh_s_pre(ji) > 0._wp ) THEN
          zindq          = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) )
          zdh_s_mel (ji) = - zindq * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
          zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting 
-         ! Heat flux associated with snow melt of falling snow (W.m-2, <0)
-         hfx_snw_1d(ji) = hfx_snw_1d(ji) + zdh_s_mel(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice 
          ! heat used to melt snow (W.m-2, >0)
-         hfx_tot_1d(ji) = hfx_tot_1d(ji) - zdh_s_mel(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
-         ! snow melting only = water into the ocean (then without snow precip)
-         wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdh_s_mel(ji) * r1_rdtice
+         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
+         ! snow melting only = water into the ocean (then without snow precip), >0
+         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdh_s_mel(ji) * r1_rdtice
          
          ! updates available heat + thickness
@@ -275,5 +275,5 @@
          ! clem debug: variation of enthalpy (J.m-2)
          dq_s(ji) = dq_s(ji) + ( zdh_s_pre(ji) + zdh_s_mel(ji) ) * zqprec(ji)  
-
+         ENDIF
       END DO
 
@@ -288,10 +288,8 @@
             zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting
             zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)    
-            ! heat flux associated with snow melt(W.m-2, <0)
-            hfx_snw_1d(ji)   = hfx_snw_1d(ji) + zdeltah(ji,jk) * a_i_b(ji) * q_s_b(ji,jk) * r1_rdtice
             ! heat used to melt snow(W.m-2, >0)
-            hfx_tot_1d(ji)   = hfx_tot_1d(ji) - zdeltah(ji,jk) * a_i_b(ji) * q_s_b(ji,jk) * r1_rdtice 
+            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_b(ji) * q_s_b(ji,jk) * r1_rdtice 
             ! snow melting only = water into the ocean (then without snow precip)
-            wfx_snw_1d(ji)   = wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
+            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
 
             ! updates available heat + thickness
@@ -308,4 +306,6 @@
       !----------------------
       ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
+      ! clem comment: not counted in mass exchange in limsbc since this is an exchange with atm. (not ocean)
+      ! clem comment: ice should also sublimate
       IF( lk_cpl ) THEN
          ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice)
@@ -320,10 +320,7 @@
                &  ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_b(ji,1) )  &
                &  * a_i_b(ji) * r1_rdtice
-            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff ! diag only (to close heat budget)
-            ! heat used for sublimation (>0, W.m-2)
-            !!? hfx_tot_1d(ji) = hfx_tot_1d(ji) - zcoeff
+            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff
             ! Mass flux by sublimation
-            wfx_sub_1d(ji) =  wfx_sub_1d(ji) + rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice ! diag only
-            wfx_snw_1d(ji) =  wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice
+            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice
             ! new snow thickness
             ht_s_b(ji)     =  MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) )
@@ -388,9 +385,9 @@
             hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
 
-            ! Total heat flux used in this process [W.m-2], < 0  
-            hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
+            ! Total heat flux used in this process [W.m-2], > 0  
+            hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
 
             ! Contribution to mass flux
-            wfx_sum_1d(ji) =  wfx_sum_1d(ji) + rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
+            wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
            
             ! record which layers have disappeared (for bottom melting) 
@@ -424,5 +421,5 @@
       ! Basal growth is driven by heat imbalance at the ice-ocean interface,
       ! between the inner conductive flux  (fc_bo_i), from the open water heat flux 
-      ! (fhldb) and the turbulent ocean flux (fhtur). 
+      ! (fhld) and the turbulent ocean flux (fhtur). 
       ! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice 
 
@@ -490,16 +487,29 @@
             ! New ice growth
                                     
-            zfmdt          = - rhoic * dh_i_bott(ji)                       ! Mass flux x time step (kg/m2, < 0)
+            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0)
+
+            ztmelts        = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K)
+            
+            zt_i_new       = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)
+            
+            zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)      
+               &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   &
+               &               + rcp  * ( ztmelts-rtt )          
+            
+            zEw            = rcp  * ( t_bo_b(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
+            
+            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
             
             ! Contribution to heat flux to the ocean [W.m-2], >0  
             hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
-            ! Total heat flux used in this process [W.m-2]  
-            hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
+
+            ! Total heat flux used in this process [W.m-2], <0  
+            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
             
-            ! Contribution to salt flux  ()
+            ! Contribution to salt flux, <0
             sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt * r1_rdtice
 
-            ! Contribution to mass flux
-            wfx_bog_1d(ji) =  wfx_bog_1d(ji) + rhoic * a_i_b(ji) * dh_i_bott(ji) * r1_rdtice
+            ! Contribution to mass flux, <0
+            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_b(ji) * dh_i_bott(ji) * r1_rdtice
 
             ! clem debug: variation of enthalpy (J.m-2)
@@ -542,4 +552,10 @@
                   hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_b(ji) * zEi * r1_rdtice
 
+                  ! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
+                  sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
+                                    
+                  ! Contribution to mass flux
+                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
+
                   ! clem debug: variation of enthalpy (J.m-2)
                   dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)  
@@ -573,4 +589,13 @@
                   ! Contribution to heat flux to the ocean [W.m-2], <0  
                   hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
+
+                  ! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
+                  sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
+                  
+                  ! Total heat flux used in this process [W.m-2], >0  
+                  hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
+                  
+                  ! Contribution to mass flux
+                  wfx_bom_1d(ji) =  wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
 
                   ! clem debug: variation of enthalpy (J.m-2)
@@ -581,13 +606,4 @@
                   h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
                ENDIF
-
-               ! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
-               sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
-
-               ! Total heat flux used in this process [W.m-2]  
-               hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
-
-               ! Contribution to mass flux
-               wfx_bom_1d(ji) =  wfx_bom_1d(ji) + rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
            
             ENDIF
@@ -602,21 +618,21 @@
       ! clem bug: I think this should be included above, so we would not have to 
       !           track heat/salt/mass fluxes backwards
-      IF( jpl == 1 ) THEN
-         DO ji = kideb, kiut
-            IF(  zf_tt(ji)  >=  0._wp  ) THEN
-               zdh            = MAX( hmelt , dh_i_bott(ji) )
-               zdvres         = zdh - dh_i_bott(ji) ! >=0
-               dh_i_bott(ji)  = zdh
-
-               ! excessive energy is sent to lateral ablation
-               zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_b(ji) - epsi20 ) )
-               zq_1cat(ji) =  zinda * rhoic * lfus * at_i_b(ji) / MAX( 1._wp - at_i_b(ji) , epsi20 ) * zdvres ! J.m-2 >=0
-
-               ! correct salt and mass fluxes
-               sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation
-               wfx_bom_1d(ji) = wfx_bom_1d(ji) + rhoic * a_i_b(ji) * zdvres * r1_rdtice
-            ENDIF
-         END DO
-      ENDIF
+!      IF( jpl == 1 ) THEN
+!         DO ji = kideb, kiut
+!            IF(  zf_tt(ji)  >=  0._wp  ) THEN
+!               zdh            = MAX( hmelt , dh_i_bott(ji) )
+!               zdvres         = zdh - dh_i_bott(ji) ! >=0
+!               dh_i_bott(ji)  = zdh
+!
+!               ! excessive energy is sent to lateral ablation
+!               zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_b(ji) - epsi20 ) )
+!               zq_1cat(ji) =  zinda * rhoic * lfus * at_i_b(ji) / MAX( 1._wp - at_i_b(ji) , epsi20 ) * zdvres ! J.m-2 >=0
+!
+!               ! correct salt and mass fluxes
+!               sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation
+!               wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdvres * r1_rdtice
+!            ENDIF
+!         END DO
+!      ENDIF
 
       !-------------------------------------------
@@ -644,10 +660,8 @@
 !        
 !         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji)                ! update available heat (J.m-2)
-!         ! Heat flux associated with snow melt
-!         hfx_snw_1d(ji)  = hfx_snw_1d(ji) + zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (<0)
 !         ! heat used to melt snow
-!         hfx_tot_1d(ji)  = hfx_tot_1d(ji) - zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)
+!         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)
 !         ! Contribution to mass flux
-!         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdeltah(ji,1) * r1_rdtice
+!         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,1) * r1_rdtice
 !         ! clem debug: variation of enthalpy (J.m-2)
 !         dq_s(ji) = dq_s(ji) + zdeltah(ji,1) * q_s_b(ji,1)  
@@ -680,10 +694,10 @@
          ! new salinity difference stored (to be used in limthd_ent.F90)
          IF (  num_sal == 2  ) THEN
-            i_ice_switch = MAX( 0._wp , SIGN( 1._wp , ht_i_b(ji) - epsi10 ) )
+            zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_b(ji) - epsi10 ) )
             ! salinity dif due to snow-ice formation
-            dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( ht_i_b(ji), epsi10 ) * i_ice_switch     
+            dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch     
             ! salinity dif due to bottom growth 
             IF (  zf_tt(ji)  < 0._wp ) THEN
-               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( ht_i_b(ji), epsi10 ) * i_ice_switch
+               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch
             ENDIF
          ENDIF
@@ -704,6 +718,6 @@
          ! Contribution to mass flux
          ! All snow is thrown in the ocean, and seawater is taken to replace the volume
-         wfx_sni_1d(ji) = wfx_sni_1d(ji) + a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice
-         wfx_snw_1d(ji) = wfx_snw_1d(ji) - a_i_b(ji) * dh_snowice(ji) * rhosn * r1_rdtice
+         wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice
+         wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_b(ji) * dh_snowice(ji) * rhosn * r1_rdtice
 
          ! clem debug: variation of enthalpy (J.m-2)
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4649)
@@ -32,5 +32,5 @@
    PUBLIC   lim_thd_dif   ! called by lim_thd
 
-   REAL(wp) ::   epsi10      =  1.e-10_wp    !
+   REAL(wp) ::   epsi10 = 1.e-10_wp    !
    !!----------------------------------------------------------------------
    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
@@ -92,7 +92,6 @@
       !!           (04-2007) Energy conservation tested by M. Vancoppenolle
       !!------------------------------------------------------------------
-      INTEGER , INTENT (in) ::   kideb   ! Start point on which the  the computation is applied
-      INTEGER , INTENT (in) ::   kiut    ! End point on which the  the computation is applied
-      INTEGER , INTENT (in) ::   jl      ! Category number
+      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
+      INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number
 
       !! * Local variables
@@ -103,7 +102,7 @@
       INTEGER ::   nconv       ! number of iterations in iterative procedure
       INTEGER ::   minnumeqmin, maxnumeqmax
-      INTEGER, DIMENSION(kiut) ::   numeqmin   ! reference number of top equation
-      INTEGER, DIMENSION(kiut) ::   numeqmax   ! reference number of bottom equation
-      INTEGER, DIMENSION(kiut) ::   isnow      ! switch for presence (1) or absence (0) of snow
+      INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation
+      INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation
+      INTEGER, POINTER, DIMENSION(:) ::   isnow      ! switch for presence (1) or absence (0) of snow
       REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
       REAL(wp) ::   zg1       =  2._wp        !
@@ -115,39 +114,46 @@
       REAL(wp) ::   ztmelt_i    ! ice melting temperature
       REAL(wp) ::   zerritmax   ! current maximal error on temperature 
-      REAL(wp), DIMENSION(kiut) ::   ztfs        ! ice melting point
-      REAL(wp), DIMENSION(kiut) ::   ztsuold     ! old surface temperature (before the iterative procedure )
-      REAL(wp), DIMENSION(kiut) ::   ztsuoldit   ! surface temperature at previous iteration
-      REAL(wp), DIMENSION(kiut) ::   zh_i        ! ice layer thickness
-      REAL(wp), DIMENSION(kiut) ::   zh_s        ! snow layer thickness
-      REAL(wp), DIMENSION(kiut) ::   zfsw        ! solar radiation absorbed at the surface
-      REAL(wp), DIMENSION(kiut) ::   zf          ! surface flux function
-      REAL(wp), DIMENSION(kiut) ::   dzf         ! derivative of the surface flux function
-      REAL(wp), DIMENSION(kiut) ::   zerrit      ! current error on temperature
-      REAL(wp), DIMENSION(kiut) ::   zdifcase    ! case of the equation resolution (1->4)
-      REAL(wp), DIMENSION(kiut) ::   zftrice     ! solar radiation transmitted through the ice
-      REAL(wp), DIMENSION(kiut) ::   zihic, zhsu
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztcond_i    ! Ice thermal conductivity
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradtr_i    ! Radiation transmitted through the ice
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradab_i    ! Radiation absorbed in the ice
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zkappa_i    ! Kappa factor in the ice
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztiold      ! Old temperature in the ice
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zeta_i      ! Eta factor in the ice
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztitemp     ! Temporary temperature in the ice to check the convergence
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zspeche_i   ! Ice specific heat
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   z_i         ! Vertical cotes of the layers in the ice
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradtr_s    ! Radiation transmited through the snow
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradab_s    ! Radiation absorbed in the snow
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zkappa_s    ! Kappa factor in the snow
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zeta_s       ! Eta factor in the snow
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztstemp      ! Temporary temperature in the snow to check the convergence
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztsold       ! Temporary temperature in the snow
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   z_s          ! Vertical cotes of the layers in the snow
-      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindterm   ! Independent term
-      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindtbis   ! temporary independent term
-      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zdiagbis
-      REAL(wp), DIMENSION(kiut,jkmax+2,3) ::   ztrid   ! tridiagonal system terms
-      REAL(wp) ::   ztemp   ! local scalar 
+      REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! ice melting point
+      REAL(wp), POINTER, DIMENSION(:) ::   ztsuold     ! old surface temperature (before the iterative procedure )
+      REAL(wp), POINTER, DIMENSION(:) ::   ztsuoldit   ! surface temperature at previous iteration
+      REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness
+      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness
+      REAL(wp), POINTER, DIMENSION(:) ::   zfsw        ! solar radiation absorbed at the surface
+      REAL(wp), POINTER, DIMENSION(:) ::   zf          ! surface flux function
+      REAL(wp), POINTER, DIMENSION(:) ::   dzf         ! derivative of the surface flux function
+      REAL(wp), POINTER, DIMENSION(:) ::   zerrit      ! current error on temperature
+      REAL(wp), POINTER, DIMENSION(:) ::   zdifcase    ! case of the equation resolution (1->4)
+      REAL(wp), POINTER, DIMENSION(:) ::   zftrice     ! solar radiation transmitted through the ice
+      REAL(wp), POINTER, DIMENSION(:) ::   zihic, zhsu
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztcond_i    ! Ice thermal conductivity
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_i    ! Radiation transmitted through the ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_i    ! Radiation absorbed in the ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_i    ! Kappa factor in the ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztiold      ! Old temperature in the ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_i      ! Eta factor in the ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztitemp     ! Temporary temperature in the ice to check the convergence
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zspeche_i   ! Ice specific heat
+      REAL(wp), POINTER, DIMENSION(:,:) ::   z_i         ! Vertical cotes of the layers in the ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_s    ! Radiation transmited through the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_s    ! Radiation absorbed in the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_s    ! Kappa factor in the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_s       ! Eta factor in the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztstemp      ! Temporary temperature in the snow to check the convergence
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztsold       ! Temporary temperature in the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   z_s          ! Vertical cotes of the layers in the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zindterm   ! Independent term
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zindtbis   ! temporary independent term
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiagbis
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid   ! tridiagonal system terms
       !!------------------------------------------------------------------     
       ! 
+      CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow )
+      CALL wrk_alloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw )
+      CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu )
+      CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0)
+      CALL wrk_alloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart=0)
+      CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis  )
+      CALL wrk_alloc( jpij, jkmax+2, 3, ztrid )
+
       !------------------------------------------------------------------------------!
       ! 1) Initialization                                                            !
@@ -318,15 +324,13 @@
             DO layer = 1, nlay_i-1
                DO ji = kideb , kiut
-                  ztemp = t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2._wp * rtt
-                  ztcond_i(ji,layer) = rcdic + 0.0900_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   &
-                     &                                   / MIN( -2.0_wp * epsi10, ztemp )   &
-                     &                       - 0.0055_wp * ztemp
+                  ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   &
+                     &                                  / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   &
+                     &                       - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )  
                   ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin )
                END DO
             END DO
             DO ji = kideb , kiut
-               ztemp = t_bo_b(ji) - rtt
-               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN( -epsi10, ztemp )   &
-                  &                        - 0.011_wp * ztemp  
+               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt)   &
+                  &                        - 0.011_wp * ( t_bo_b(ji) - rtt )  
                ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
             END DO
@@ -395,5 +399,4 @@
          !
          DO ji = kideb , kiut
-
             ! update of the non solar flux according to the update in T_su
             qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) )
@@ -403,7 +406,4 @@
                + qns_ice_1d(ji)                  ! non solar total flux 
             ! (LWup, LWdw, SH, LH)
-
-            ! heat flux used to change surface temperature
-            !hfx_tot_1d(ji) = hfx_tot_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) ) * a_i_b(ji)
          END DO
 
@@ -721,11 +721,19 @@
       DO ji = kideb, kiut
          IF( t_su_b(ji) < rtt ) THEN  ! case T_su < 0degC
-            hfx_tot_1d(ji) = hfx_tot_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
+            hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
          ELSE                                    ! case T_su = 0degC
-            hfx_tot_1d(ji) = hfx_tot_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
+            hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
          ENDIF
       END DO
 
       !
+      CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow )
+      CALL wrk_dealloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw )
+      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu )
+      CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )
+      CALL wrk_dealloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart = 0 )
+      CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis )
+      CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid )
+
    END SUBROUTINE lim_thd_dif
 
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 4649)
@@ -48,5 +48,5 @@
 CONTAINS
  
-   SUBROUTINE lim_thd_ent( kideb, kiut, jl )
+   SUBROUTINE lim_thd_ent( kideb, kiut, jl, qnew )
       !!-------------------------------------------------------------------
       !!               ***   ROUTINE lim_thd_ent  ***
@@ -62,5 +62,12 @@
       !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses
       !!            2) linear remapping on the new layers
-      !!            3) Ice salinity update + recover temperature from enthalpies
+      !!
+      !! ------------ cum0(0)                        ------------- cum1(0)
+      !!                                    NEW      -------------
+      !! ------------ cum0(1)               ==>      -------------
+      !!     ...                                     -------------
+      !! ------------                                -------------
+      !! ------------ cum0(nlay_i+2)                 ------------- cum1(nlay_i)
+      !!
       !!
       !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
@@ -68,4 +75,6 @@
       INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
       INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number
+
+      REAL(wp), INTENT(inout), DIMENSION(:,:) :: qnew          ! new enthlapies (remapped)
 
       INTEGER  :: ji,ii,ij   !  dummy loop indices
@@ -76,8 +85,10 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: zqh_cum0, zh_cum0   ! old cumulative enthlapies and layers interfaces
       REAL(wp), POINTER, DIMENSION(:,:) :: zqh_cum1, zh_cum1   ! new cumulative enthlapies and layers interfaces
+      REAL(wp), POINTER, DIMENSION(:)   :: zhnew               ! new layers thicknesses
       !!-------------------------------------------------------------------
 
       CALL wrk_alloc( jpij, nlay_i+3, zqh_cum0, zh_cum0, kjstart = 0 )
       CALL wrk_alloc( jpij, nlay_i+1, zqh_cum1, zh_cum1, kjstart = 0 )
+      CALL wrk_alloc( jpij, zhnew )
 
       !--------------------------------------------------------------------------
@@ -96,9 +107,14 @@
       !  2) Interpolation on the new layers
       !------------------------------------
+      ! new layer thickesses
+      DO ji = kideb, kiut
+         zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) / REAL( nlay_i )  
+      ENDDO
+
       ! new layers interfaces
       zh_cum1(:,0:nlay_i) = 0._wp
       DO jk1 = 1, nlay_i
          DO ji = kideb, kiut
-            zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + ht_i_b(ji) / REAL( nlay_i )
+            zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji)
          ENDDO
       ENDDO
@@ -123,37 +139,12 @@
       DO jk1 = 1, nlay_i
          DO ji = kideb, kiut
-            zswitch       =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) + epsi20 ) ) 
-            q_i_b(ji,jk1) = zswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) * REAL( nlay_i ) / MAX( ht_i_b(ji), epsi20 )
+            zswitch      =  1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi10 ) ) 
+            qnew(ji,jk1) = zswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi10 )
          ENDDO
       ENDDO
-
-      !---------------------------------------------------------
-      !  3) Update ice salinity and recover ice temperature
-      !---------------------------------------------------------
-      ! Update salinity (basal entrapment, snow ice formation)
-      DO ji = kideb, kiut
-         sm_i_b(ji) = sm_i_b(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji)
-      END DO !ji
-
-      ! Recover ice temperature
-      DO jk1 = 1, nlay_i
-         DO ji = kideb, kiut
-            ztmelts       =  -tmut * s_i_b(ji,jk1) + rtt
-            ! Conversion q(S,T) -> T (second order equation)
-            zaaa          =  cpic
-            zbbb          =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk1) / rhoic - lfus
-            zccc          =  lfus * ( ztmelts - rtt )
-            zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) )
-            t_i_b(ji,jk1) =  rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )
-            
-            ! mask temperature
-            zswitch       =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) ) 
-            t_i_b(ji,jk1) =  zswitch * t_i_b(ji,jk1) + ( 1._wp - zswitch ) * rtt
-         END DO 
-      END DO 
-
       !
       CALL wrk_dealloc( jpij, nlay_i+3, zqh_cum0, zh_cum0, kjstart = 0 )
       CALL wrk_dealloc( jpij, nlay_i+1, zqh_cum1, zh_cum1, kjstart = 0 )
+      CALL wrk_dealloc( jpij, zhnew )
       !
    END SUBROUTINE lim_thd_ent
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4649)
@@ -30,4 +30,5 @@
    USE wrk_nemo       ! work arrays
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE limthd_ent
 
    IMPLICIT NONE
@@ -85,6 +86,8 @@
       REAL(wp) ::   zEw          ! seawater specific enthalpy (J/kg)
       REAL(wp) ::   zfmdt        ! mass flux x time step (kg/m2, >0 towards ocean)
-       
-      INTEGER , POINTER, DIMENSION(:) ::   zcatac      ! indexes of categories where new ice grows
+     
+      REAL(wp) ::   zv_newfra
+  
+      INTEGER , POINTER, DIMENSION(:) ::   jcat      ! indexes of categories where new ice grows
       REAL(wp), POINTER, DIMENSION(:) ::   zswinew     ! switch for new ice or not
 
@@ -97,50 +100,27 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zdv_res     ! residual volume in case of excessive heat budget
       REAL(wp), POINTER, DIMENSION(:) ::   zda_res     ! residual area in case of excessive heat budget
-      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_ac    ! total ice fraction    
+      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_1d    ! total ice fraction    
       REAL(wp), POINTER, DIMENSION(:) ::   zat_i_lev   ! total ice fraction for level ice only (type 1)   
       REAL(wp), POINTER, DIMENSION(:) ::   zv_frazb   ! accretion of frazil ice at the ice bottom
-      REAL(wp), POINTER, DIMENSION(:) ::   zvrel_ac    ! relative ice / frazil velocity (1D vector)
+      REAL(wp), POINTER, DIMENSION(:) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector)
 
       REAL(wp), POINTER, DIMENSION(:,:) ::   zv_old      ! old volume of ice in category jl
       REAL(wp), POINTER, DIMENSION(:,:) ::   za_old      ! old area of ice in category jl
-      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_ac     ! 1-D version of a_i
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_ac     ! 1-D version of v_i
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_ac    ! 1-D version of oa_i
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_ac   ! 1-D version of smv_i
-
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_ac   !: 1-D version of e_i
-
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqm0      ! old layer-system heat content
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zthick0   ! old ice thickness
-
-      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final   ! ice volume summed over categories
-      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   !  snow volume summed over categories
-      REAL(wp), POINTER, DIMENSION(:,:) ::   et_i_init, et_i_final   !  ice energy summed over categories
-      REAL(wp), POINTER, DIMENSION(:,:) ::   et_s_init               !  snow energy summed over categories
+      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_1d     ! 1-D version of a_i
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_1d     ! 1-D version of v_i
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_1d    ! 1-D version of oa_i
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d   ! 1-D version of smv_i
+
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_1d   !: 1-D version of e_i
+
       REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel                   ! relative ice / frazil velocity
       !!-----------------------------------------------------------------------!
 
-      CALL wrk_alloc( jpij, zcatac )   ! integer
+      CALL wrk_alloc( jpij, jcat )   ! integer
       CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
-      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zv_frazb, zvrel_ac )
-      CALL wrk_alloc( jpij,jpl, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )
-      CALL wrk_alloc( jpij,jkmax,jpl, ze_i_ac )
-      CALL wrk_alloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
-      CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel )
-
-      et_i_init(:,:) = 0._wp
-      et_s_init(:,:) = 0._wp
-      vt_i_init(:,:) = 0._wp
-      vt_s_init(:,:) = 0._wp
-
-      !------------------------------------------------------------------------------!
-      ! 1) Conservation check and changes in each ice category
-      !------------------------------------------------------------------------------!
-      IF( con_i ) THEN
-         CALL lim_column_sum        ( jpl, v_i          , vt_i_init)
-         CALL lim_column_sum        ( jpl, v_s          , vt_s_init)
-         CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init)
-         CALL lim_column_sum        ( jpl, e_s(:,:,1,:) , et_s_init)
-      ENDIF
+      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d )
+      CALL wrk_alloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )
+      CALL wrk_alloc( jpij,jkmax,jpl, ze_i_1d )
+      CALL wrk_alloc( jpi,jpj, zvrel )
 
       !------------------------------------------------------------------------------|
@@ -204,11 +184,12 @@
                      &          +   vtau_ice(ji  ,jj  ) * tmv(ji  ,jj  ) ) * 0.5_wp
                   ! Square root of wind stress
-                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
+                  ztenagm       =  SQRT( SQRT( ztaux**2 + ztauy**2 ) )
 
                   !---------------------
                   ! Frazil ice velocity
                   !---------------------
-                  zvfrx         = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10)
-                  zvfry         = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10)
+                  zindb = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
+                  zvfrx = zindb * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
+                  zvfry = zindb * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
 
                   !-------------------
@@ -305,12 +286,12 @@
       IF ( nbpac > 0 ) THEN
 
-         CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, zat_i_1d  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) )
          DO jl = 1, jpl
-            CALL tab_2d_1d( nbpac, za_i_ac  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
-            CALL tab_2d_1d( nbpac, zv_i_ac  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
-            CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
-            CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
+            CALL tab_2d_1d( nbpac, za_i_1d  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
+            CALL tab_2d_1d( nbpac, zv_i_1d  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
+            CALL tab_2d_1d( nbpac, zoa_i_1d (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
+            CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
             DO jk = 1, nlay_i
-               CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
+               CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
             END DO ! jk
          END DO ! jl
@@ -322,8 +303,8 @@
          CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac)     , wfx_opw, jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, zvrel_1d  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) )
 
          CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac)     , hfx_thd, jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, hfx_tot_1d(1:nbpac)     , hfx_tot, jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac)     , hfx_opw, jpi, jpj, npac(1:nbpac) )
 
          !------------------------------------------------------------------------------!
@@ -334,6 +315,6 @@
          ! Keep old ice areas and volume in memory
          !-----------------------------------------
-         zv_old(:,:) = zv_i_ac(:,:) 
-         za_old(:,:) = za_i_ac(:,:)
+         zv_old(:,:) = zv_i_1d(:,:) 
+         za_old(:,:) = za_i_1d(:,:)
 
          !----------------------
@@ -343,5 +324,5 @@
             zh_newice(ji) = hiccrit(1)
          END DO
-         IF( fraz_swi == 1.0 )   zh_newice(:) = hicol_b(:)
+         IF( fraz_swi == 1.0 ) zh_newice(:) = hicol_b(:)
 
          !----------------------
@@ -370,10 +351,6 @@
             ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K)
             ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             &
-               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   &
+               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_b(ji) - rtt, -epsi10 ) )   &
                &                       - rcp  *         ( ztmelts - rtt )  )
-            ! MV HC 2014 comment I dont see why this line below is here... ?
-            ! This implies that ze_newice gets to rhoic*Lfus if it was negative, but this should never happen
-            ze_newice(ji) =   MAX( ze_newice(ji) , 0._wp )    &
-               &          +   MAX(  0.0 , SIGN( 1.0 , - ze_newice(ji) )  ) * rhoic * lfus
          END DO ! ji
          !----------------
@@ -405,31 +382,17 @@
             hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice
             ! Total heat flux used in this process [W.m-2]  
-            hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * zdE * r1_rdtice
+            hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice
             ! mass flux
-            wfx_opw_1d(ji) = wfx_opw_1d(ji) + zv_newice(ji) * rhoic * r1_rdtice
+            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice
             ! salt flux
             sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice
 
             ! A fraction zfrazb of frazil ice is accreted at the ice bottom
-            zfrazb        = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
+            zinda         = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
+            zfrazb        = zinda * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
             zv_frazb(ji)  =         zfrazb   * zv_newice(ji)
             zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
          END DO
 
-         !------------------------------------
-         ! Diags for energy conservation test
-         !------------------------------------
-         DO ji = 1, nbpac
-            ii = MOD( npac(ji) - 1 , jpi ) + 1
-            ij =    ( npac(ji) - 1 ) / jpi + 1
-            !
-            zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji)
-            !zde = ze_newice(ji) * area(ii,ij) * zv_newice(ji)
-            !
-            ! clem: change that?
-            vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji)             ! volume
-            et_i_init(ii,ij) = et_i_init(ii,ij) + zde                       ! Energy
-
-         END DO
 
          !-----------------
@@ -438,5 +401,5 @@
          DO ji = 1, nbpac
             za_newice(ji) = zv_newice(ji) / zh_newice(ji)
-         END DO !ji
+         END DO
 
          !------------------------------------------------------------------------------!
@@ -450,6 +413,6 @@
          ! we keep the excessive volume in memory and attribute it later to bottom accretion
          DO ji = 1, nbpac
-            IF ( za_newice(ji) >  ( amax - zat_i_ac(ji) ) ) THEN
-               zda_res(ji)   = za_newice(ji) - ( amax - zat_i_ac(ji) )
+            IF ( za_newice(ji) >  ( amax - zat_i_1d(ji) ) ) THEN
+               zda_res(ji)   = za_newice(ji) - ( amax - zat_i_1d(ji) )
                zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
                za_newice(ji) = za_newice(ji) - zda_res  (ji)
@@ -464,14 +427,13 @@
          ! Laterally redistribute new ice volume and area
          !------------------------------------------------
-         zat_i_ac(:) = 0._wp
+         zat_i_1d(:) = 0._wp
          DO jl = 1, jpl
             DO ji = 1, nbpac
-               IF(  hi_max   (jl-1)  <   zh_newice(ji)   .AND.   &
-                  & zh_newice(ji)    <=  hi_max   (jl)         ) THEN
-                  za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji)
-                  zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji)
-                  zcatac  (ji)    = jl
+               IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
+                  za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji)
+                  zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji)
+                  jcat  (ji)    = jl
                ENDIF
-               zat_i_ac(ji)    = zat_i_ac(ji)    + za_i_ac  (ji,jl)
+               zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d  (ji,jl)
             END DO
          END DO
@@ -481,13 +443,15 @@
          !----------------------------------
          DO ji = 1, nbpac
-            jl = zcatac(ji)                                                           ! categroy in which new ice is put
-            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) )   ! 0 if old ice
+            jl = jcat(ji)                                                  ! categroy in which new ice is put
+            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) ) )   ! 0 if old ice
          END DO
 
          DO jk = 1, nlay_i
             DO ji = 1, nbpac
-               jl = zcatac(ji)
-               ze_i_ac(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji)                &
-                  &      + ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_ac(ji,jk,jl) * zv_old(ji,jl) ) / zv_i_ac(ji,jl)
+               jl = jcat(ji)
+               zinda = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )
+               ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                      &
+                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_old(ji,jl) )  &
+                  &        * zinda / MAX( zv_i_1d(ji,jl), epsi20 )
             END DO
          END DO
@@ -496,67 +460,30 @@
          ! Add excessive volume of new ice at the bottom
          !-----------------------------------------------
-         jm = 1
-
-         ! ---  Redistributing energy on the new grid (energy is equally distributed in every layer) --- !
-!         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
-!            DO jk = 1, nlay_i
-!               DO ji = 1, nbpac
-!                  ze_i_ac(ji,jk,jl) = ( ze_i_ac(ji,jk,jl) * zv_i_ac(ji,jl) + ze_newice(ji) * ( zdv_res(ji) + zv_frazb(ji) ) ) / &
-!                     &                ( zv_i_ac(ji,jl) + ( zdv_res(ji) + zv_frazb(ji) ) ) 
-!               END DO
-!            END DO
-!         END DO
-
-         ! --- Redistributing energy on the new grid (energy is sent to the bottom) PART 1 --- !
-         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
+         DO jl = 1, jpl
+            h_i_old (1:nbpac,0:nlay_i+1) = 0._wp
+            qh_i_old(1:nbpac,0:nlay_i+1) = 0._wp
+
             DO jk = 1, nlay_i
                DO ji = 1, nbpac
-                  zthick0(ji,jk,jl) =  zv_i_ac(ji,jl) / REAL( nlay_i )
-                  zqm0   (ji,jk,jl) =  ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)
+                  h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i )
+                  qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk)
                END DO
             END DO
-         END DO
-         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
+
             DO ji = 1, nbpac
-               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_ac(ji) - epsi10 ) )
-               zthick0(ji,nlay_i+1,jl) =  zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_ac(ji,jl) / MAX( zat_i_ac(ji) , epsi10 )
-               zqm0   (ji,nlay_i+1,jl) =  ze_newice(ji) * zthick0(ji,nlay_i+1,jl)
-            END DO ! ji
-         END DO ! jl
-
-         ze_i_ac(:,:,:) = 0._wp
-         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
-            DO jk = 1, nlay_i
-               DO layer = 1, nlay_i + 1
-                  DO ji = 1, nbpac
-                      zindb   =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - zthick0(ji,layer,jl) + epsi10 ) ) 
-                      zweight = zindb * MAX( 0._wp,                                                                                              &
-                         &      MIN( zv_i_ac(ji,jl) * REAL( layer     ), ( zv_i_ac(ji,jl) +  zthick0(ji,nlay_i+1,jl) ) * REAL( jk     ) )     &
-                         &    - MAX( zv_i_ac(ji,jl) * REAL( layer - 1 ), ( zv_i_ac(ji,jl) +  zthick0(ji,nlay_i+1,jl) ) * REAL( jk - 1 ) ) )   &
-                         &  / ( REAL( nlay_i ) * MAX( zthick0(ji,layer,jl), epsi10 ) )
-                      ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl)  
-                  END DO 
-               END DO 
-            END DO 
-         END DO
-
-         ! --- new volumes and layer thickness ---
-         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
-            DO ji = 1, nbpac
-               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_ac(ji) - epsi10 ) )
-               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_ac(ji,jl) / MAX( zat_i_ac(ji) , epsi10 )
-            END DO
-         END DO
-
-         ! --- Redistributing energy on the new grid (energy is sent to the bottom) PART 2 --- !
-         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
-            DO jk = 1, nlay_i
-               DO ji = 1, nbpac
-                  zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 
-                  ze_i_ac(ji,jk,jl) = zindb * ze_i_ac(ji,jk,jl) / MAX( zv_i_ac(ji,jl), epsi10 ) * REAL( nlay_i )
-               END DO
-            END DO
-         END DO
-
+               zinda          = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )
+               zv_newfra      = zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )
+               za_i_1d(ji,jl) = zinda * za_i_1d(ji,jl)               
+
+               h_i_old (ji,nlay_i+1) = zv_newfra
+               qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
+
+               zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra
+            ENDDO
+
+            ! --- Ice enthalpy remapping --- !
+            CALL lim_thd_ent( 1, nbpac, jl, ze_i_1d(1:nbpac,:,jl) ) 
+
+         ENDDO
 
          !------------
@@ -565,6 +492,6 @@
          DO jl = 1, jpl
             DO ji = 1, nbpac
-               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
-               zoa_i_ac(ji,jl)  = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb   
+               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) )  ! 0 if no ice and 1 if yes
+               zoa_i_1d(ji,jl)  = za_old(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb   
             END DO 
          END DO   
@@ -576,6 +503,6 @@
             DO jl = 1, jpl
                DO ji = 1, nbpac
-                  zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
-                  zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji)
+                  zdv   = zv_i_1d(ji,jl) - zv_old(ji,jl)
+                  zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji)
                END DO
             END DO   
@@ -586,10 +513,10 @@
          !------------------------------------------------------------------------------!
          DO jl = 1, jpl
-            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ac (1:nbpac,jl), jpi, jpj )
-            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj )
-            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj )
-            CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj )
+            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj )
+            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj )
+            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_1d(1:nbpac,jl), jpi, jpj )
+            CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj )
             DO jk = 1, nlay_i
-               CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj )
+               CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_1d(1:nbpac,jk,jl), jpi, jpj )
             END DO
          END DO
@@ -599,5 +526,5 @@
 
          CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj )
-         CALL tab_1d_2d( nbpac, hfx_tot, npac(1:nbpac), hfx_tot_1d(1:nbpac), jpi, jpj )
+         CALL tab_1d_2d( nbpac, hfx_opw, npac(1:nbpac), hfx_opw_1d(1:nbpac), jpi, jpj )
          !
       ENDIF ! nbpac > 0
@@ -612,5 +539,4 @@
                   ! heat content in Joules
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ) * unit_fac ) 
-                  !e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ) ) 
                END DO
             END DO
@@ -618,43 +544,11 @@
       END DO
 
-      !------------------------------------------------------------------------------|
-      ! 10) Conservation check and changes in each ice category
-      !------------------------------------------------------------------------------|
-      IF( con_i ) THEN 
-         CALL lim_column_sum (jpl,   v_i, vt_i_final)
-         fieldid = 'v_i, limthd_lac'
-         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
-         !
-         CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)
-         fieldid = 'e_i, limthd_lac'
-         CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid) 
-         !
-         CALL lim_column_sum (jpl,   v_s, vt_s_final)
-         fieldid = 'v_s, limthd_lac'
-         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
-         !
-         !     CALL lim_column_sum (jpl,   e_s(:,:,1,:) , et_s_init)
-         !     fieldid = 'e_s, limthd_lac'
-         !     CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid) 
-         IF( ln_nicep ) THEN
-            DO ji = mi0(jiindx), mi1(jiindx)
-               DO jj = mj0(jjindx), mj1(jjindx)
-                  WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj)
-                  WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj)
-                  WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj)
-                  WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj)
-               END DO
-            END DO
-         ENDIF
-         !
-      ENDIF
       !
-      CALL wrk_dealloc( jpij, zcatac )   ! integer
+      CALL wrk_dealloc( jpij, jcat )   ! integer
       CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
-      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zv_frazb, zvrel_ac )
-      CALL wrk_dealloc( jpij,jpl, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )
-      CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac )
-      CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
-      CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel )
+      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d )
+      CALL wrk_dealloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )
+      CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_1d )
+      CALL wrk_dealloc( jpi,jpj, zvrel )
       !
    END SUBROUTINE lim_thd_lac
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 4649)
@@ -53,8 +53,14 @@
       !
       INTEGER  ::   ji, jk     ! dummy loop indices 
-      REAL(wp) ::   iflush, igravdr, ztmelts   ! local scalars
-      REAL(wp) ::   zaaa, zbbb, zccc, zdiscrim   ! local scalars
+      REAL(wp) ::   iflush, igravdr   ! local scalars
       !!---------------------------------------------------------------------
 
+      !---------------------------------------------------------
+      !  0) Update ice salinity from snow-ice and bottom growth
+      !---------------------------------------------------------
+      DO ji = kideb, kiut
+         sm_i_b(ji) = sm_i_b(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji)
+      END DO
+ 
       !------------------------------------------------------------------------------|
       ! 1) Constant salinity, constant in time                                       |
@@ -71,5 +77,4 @@
       !  Module 2 : Constant salinity varying in time                                |
       !------------------------------------------------------------------------------|
-
       IF(  num_sal == 2  ) THEN
 
@@ -78,14 +83,15 @@
             ! Switches 
             !----------
-            iflush       =         MAX( 0._wp , SIGN( 1._wp , t_su_b(ji) - rtt )        )    ! =1 if summer 
-            igravdr      =         MAX( 0._wp , SIGN( 1._wp , t_bo_b(ji) - t_su_b(ji) ) )    ! =1 if t_su < t_bo
+            iflush  = MAX( 0._wp , SIGN( 1._wp , t_su_b(ji) - rtt )        )    ! =1 if summer 
+            igravdr = MAX( 0._wp , SIGN( 1._wp , t_bo_b(ji) - t_su_b(ji) ) )    ! =1 if t_su < t_bo
 
             !---------------------
             ! Salinity tendencies
             !---------------------
-            !                               
-            dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice  ! drainage by gravity
-            dsm_i_fl_1d(ji) = - iflush  * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice  ! drainage by flushing 
-            !
+            ! drainage by gravity drainage
+            dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice 
+            ! drainage by flushing  
+            dsm_i_fl_1d(ji) = - iflush  * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice
+
             !-----------------
             ! Update salinity   
@@ -104,21 +110,4 @@
          ! Salinity profile
          CALL lim_var_salprof1d( kideb, kiut )
-
-
-         ! Only necessary for conservation check since salinity is modified
-         !--------------------
-         ! Temperature update
-         !--------------------
-         DO jk = 1, nlay_i
-            DO ji = kideb, kiut
-               ztmelts    =  -tmut*s_i_b(ji,jk) + rtt
-               !Conversion q(S,T) -> T (second order equation)
-               zaaa         =  cpic
-               zbbb         =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus
-               zccc         =  lfus * ( ztmelts - rtt )
-               zdiscrim     =  SQRT(  MAX( zbbb*zbbb - 4.0*zaaa*zccc, 0._wp )  )
-               t_i_b(ji,jk) =  rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa )
-            END DO
-         END DO
          !
       ENDIF 
@@ -127,5 +116,4 @@
       !  Module 3 : Profile of salinity, constant in time                            |
       !------------------------------------------------------------------------------|
-
       IF(  num_sal == 3  )   CALL lim_var_salprof1d( kideb, kiut )
 
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 4649)
@@ -30,4 +30,5 @@
    USE limvar          ! clem for ice thickness correction
    USE timing          ! Timing
+   USE limcons        ! conservation tests
 
    IMPLICIT NONE
@@ -72,6 +73,4 @@
       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi
       REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e
-      REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset)
-      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax, zchk_umax ! Check errors (C Rousset)
       ! mass and salt flux (clem)
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold   ! old ice volume...
@@ -79,4 +78,6 @@
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zeiold, zesold   ! old enthalpies
       REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei
+      !
+      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!---------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limtrp')
@@ -87,17 +88,4 @@
 
       CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold )   ! clem
-
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF( ln_limdiahsb ) THEN
-         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) )
-         zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) )
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
 
       IF( numit == nstart .AND. lwp ) THEN
@@ -114,4 +102,8 @@
       IF( ln_limdyn ) THEN          !   Advection of sea ice properties   !
          !                          !-------------------------------------!
+
+         ! conservation test
+         IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
          ! mass and salt flux init (clem)
          zviold(:,:,:)  = v_i(:,:,:)
@@ -368,6 +360,6 @@
 
                  ! Update fluxes
-                  wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 
-                  wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
+                  wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 
+                  wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
                   sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
                   hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
@@ -425,6 +417,6 @@
 
                      ! Update mass fluxes
-                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
-                     wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
+                     wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
+                     wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
                      sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
                      hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
@@ -474,4 +466,7 @@
             END DO
          END DO      
+
+         ! conservation test
+         IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       ENDIF
@@ -508,34 +503,4 @@
          END DO
       ENDIF
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF( ln_limdiahsb ) THEN
-         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
-         zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b
-         zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b
- 
-         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 
-         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
-         zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft
-
-         zchk_vmin = glob_min(v_i)
-         zchk_amax = glob_max(SUM(a_i,dim=3))
-         zchk_amin = glob_min(a_i)
-         zchk_umax = glob_max(SQRT(u_ice**2 + v_ice**2))
-
-         IF(lwp) THEN
-            IF ( ABS( zchk_v_i   ) >  1.e-4 ) THEN
-               WRITE(numout,*) 'violation volume [kg/day]     (limtrp) = ',(zchk_v_i * rday)
-               WRITE(numout,*) 'u_ice max [m/s]               (limtrp) = ',zchk_umax
-               WRITE(numout,*) 'number of time steps          (limtrp) =',kt
-            ENDIF
-            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * rday)
-            IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limtrp) = ',(zchk_e_i)
-            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limtrp) = ',(zchk_vmin * 1.e-3)
-            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limtrp) = ',zchk_amin
-         ENDIF
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
       !
       CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold )
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 4649)
@@ -38,9 +38,9 @@
    USE wrk_nemo         ! work arrays
    USE lib_fortran     ! glob_sum
-   ! Check budget (Rousset)
    USE in_out_manager   ! I/O manager
    USE iom              ! I/O manager
    USE lib_mpp          ! MPP library
    USE timing          ! Timing
+   USE limcons        ! conservation tests
 
    IMPLICIT NONE
@@ -82,31 +82,15 @@
       REAL(wp) ::   zh, zdvres, zsal
       REAL(wp) ::   zat_i_old, ztmelts
-
-      REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset)
-      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
+      !
+      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!-------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limupdate1')
 
-      !------------------------------------------------------------------------------
-      ! 1. Update of Global variables                                               |
-      !------------------------------------------------------------------------------
-
-      !-----------------
-      !  Trend terms
-      !-----------------
-
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF (ln_limdiahsb) THEN
-         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) )
-         zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) )
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
-
+      IF( ln_limdyn ) THEN 
+
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
+      !-----------------
       ! zap small values
       !-----------------
@@ -115,11 +99,6 @@
       CALL lim_var_glo2eqv
      
-      at_i(:,:) = 0._wp
-      DO jl = 1, jpl
-         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
-      END DO
-
       !----------------------------------------------------
-      ! 2.2) Rebin categories with thickness out of bounds
+      ! Rebin categories with thickness out of bounds
       !----------------------------------------------------
       DO jm = 1, jpm
@@ -134,5 +113,6 @@
       END DO
 
-      !--- 2.13 ice concentration should not exceed amax 
+      !----------------------------------------------------
+      ! ice concentration should not exceed amax 
       !-----------------------------------------------------
       DO jl  = 1, jpl
@@ -147,10 +127,10 @@
       END DO
 
-      at_i(:,:) = 0.0
+      at_i(:,:) = 0._wp
       DO jl = 1, jpl
          at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
       END DO
-  
-  
+    
+      ! --------------------------------------
       ! Final thickness distribution rebinning
       ! --------------------------------------
@@ -163,5 +143,5 @@
       END DO
 
-
+      !-----------------
       ! zap small values
       !-----------------
@@ -169,5 +149,5 @@
 
       !---------------------
-      ! 2.11) Ice salinity
+      ! Ice salinity bounds
       !---------------------
       IF (  num_sal == 2  ) THEN 
@@ -179,18 +159,11 @@
                   ! salinity stays in bounds
                   i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) )
-                  smv_i(ji,jj,jl) = i_ice_switch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl)
+                  smv_i(ji,jj,jl) = i_ice_switch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) )
                   ! associated salt flux
                   sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
-               END DO ! ji
-            END DO ! jj
-         END DO !jl
+               END DO
+            END DO
+         END DO
       ENDIF
-
-      ! -------------------
-      at_i(:,:) = a_i(:,:,1)
-      DO jl = 2, jpl
-         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
-      END DO
-  
 
       ! -------------------------------------------------
@@ -206,32 +179,8 @@
       d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:)
       d_smv_i_trp(:,:,:)   = 0._wp
-      IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:)
-
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF( ln_limdiahsb ) THEN
-         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
-         zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b
-         zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b
- 
-         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 
-         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
-         zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft
-
-         zchk_vmin = glob_min(v_i)
-         zchk_amax = glob_max(SUM(a_i,dim=3))
-         zchk_amin = glob_min(a_i)
-
-         IF(lwp) THEN
-            IF ( ABS( zchk_v_i   ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (limupdate1) = ',(zchk_v_i * rday)
-            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * rday)
-            IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limupdate1) = ',(zchk_e_i)
-            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate1) = ',(zchk_vmin * 1.e-3)
-            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate1) = ',zchk_amax
-            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate1) = ',zchk_amin
-         ENDIF
-       ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
+      IF(  num_sal == 2  ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)
+
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       IF(ln_ctl) THEN   ! Control print
@@ -302,4 +251,6 @@
       ENDIF
    
+      ENDIF ! ln_limdyn
+
       IF( nn_timing == 1 )  CALL timing_stop('limupdate1')
    END SUBROUTINE lim_update1
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 4649)
@@ -39,4 +39,5 @@
    USE lib_fortran     ! glob_sum
    USE timing          ! Timing
+   USE limcons        ! conservation tests
 
    IMPLICIT NONE
@@ -84,27 +85,13 @@
       REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
       REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
-
-      REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset)
-      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
+      !
+      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!-------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limupdate2')
 
-      !----------------------------------------------------------------------------------------
-      ! 1. Computation of trend terms      
-      !----------------------------------------------------------------------------------------
-
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF (ln_limdiahsb) THEN
-         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) )
-         zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) )
-         zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) )
-      ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
-
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
+      !-----------------
       ! zap small values
       !-----------------
@@ -113,14 +100,6 @@
       CALL lim_var_glo2eqv
 
-      !--------------------------------------
-      ! 2. Review of all pathological cases
-      !--------------------------------------
-      at_i(:,:) = 0._wp
-      DO jl = 1, jpl
-         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
-      END DO
-
       !----------------------------------------------------
-      ! 2.2) Rebin categories with thickness out of bounds
+      ! Rebin categories with thickness out of bounds
       !----------------------------------------------------
       DO jm = 1, jpm
@@ -130,35 +109,6 @@
       END DO
 
-
-!clem debug: it is done in limthd_dh now
-!      ! Melt of snow
-!      !--------------
-!      DO jl = 1, jpl
-!         DO jj = 1, jpj 
-!            DO ji = 1, jpi
-!               IF( v_s(ji,jj,jl) >= epsi20 ) THEN
-!                  ! If snow energy of melting smaller then Lf
-!                  ! Then all snow melts and heat go to the ocean
-!                  !IF ( zEs <= lfus ) THEN 
-!                  IF( t_s(ji,jj,1,jl) >= rtt ) THEN
-!                     zdvres = - v_s(ji,jj,jl)
-!                     zEs    = - e_s(ji,jj,1,jl) * unit_fac / ( area(ji,jj) * MAX( v_s(ji,jj,jl), epsi20 ) )  ! snow energy of melting (J.m-3)
-!                     ! Contribution to heat flux to the ocean [W.m-2], < 0  
-!                     hfx_res(ji,jj) = hfx_res(ji,jj) - zEs * zdvres * r1_rdtice
-!                     ! Contribution to mass flux
-!                     wfx_snw(ji,jj) =  wfx_snw(ji,jj) + rhosn * zdvres * r1_rdtice
-!                     ! updates
-!                     v_s (ji,jj,jl)   = 0._wp
-!                     ht_s(ji,jj,jl)   = 0._wp
-!                     e_s (ji,jj,1,jl) = 0._wp
-!                     t_s (ji,jj,1,jl) = rtt
-!                  ENDIF
-!               ENDIF
-!            END DO
-!         END DO
-!      END DO
-!clem debug
-
-      !--- 2.12 Constrain the thickness of the smallest category above 10 cm
+      !----------------------------------------------------------------------
+      ! Constrain the thickness of the smallest category above hiclim
       !----------------------------------------------------------------------
       DO jm = 1, jpm
@@ -176,7 +126,8 @@
       END DO !jm
       
-      !--- 2.13 ice concentration should not exceed amax 
       !-----------------------------------------------------
-      at_i(:,:) = 0.0
+      ! ice concentration should not exceed amax 
+      !-----------------------------------------------------
+      at_i(:,:) = 0._wp
       DO jl = 1, jpl
          at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
@@ -199,4 +150,5 @@
       END DO
 
+      ! --------------------------------------
       ! Final thickness distribution rebinning
       ! --------------------------------------
@@ -209,4 +161,5 @@
       END DO
 
+      !-----------------
       ! zap small values
       !-----------------
@@ -232,11 +185,4 @@
       ENDIF
 
-
-      ! -------------------
-      at_i(:,:) = a_i(:,:,1)
-      DO jl = 2, jpl
-         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
-      END DO
-
       !------------------------------------------------------------------------------
       ! 2) Corrections to avoid wrong values                                        |
@@ -246,9 +192,9 @@
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
-            IF ( at_i(ji,jj) .EQ. 0.0 ) THEN ! what to do if there is no ice
-               IF ( at_i(ji+1,jj) .EQ. 0.0 ) u_ice(ji,jj)   = 0.0 ! right side
-               IF ( at_i(ji-1,jj) .EQ. 0.0 ) u_ice(ji-1,jj) = 0.0 ! left side
-               IF ( at_i(ji,jj+1) .EQ. 0.0 ) v_ice(ji,jj)   = 0.0 ! upper side
-               IF ( at_i(ji,jj-1) .EQ. 0.0 ) v_ice(ji,jj-1) = 0.0 ! bottom side
+            IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice
+               IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj)   = 0._wp ! right side
+               IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side
+               IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj)   = 0._wp ! upper side
+               IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side
             ENDIF
          END DO
@@ -261,5 +207,4 @@
       v_ice(:,:) = v_ice(:,:) * tmv(:,:)
  
-
       ! -------------------------------------------------
       ! Diagnostics
@@ -276,30 +221,14 @@
       dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday
 
-      ! -------------------------------
-      !- check conservation (C Rousset)
-      IF( ln_limdiahsb ) THEN
-         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
-         zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b
-         zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b
- 
-         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 
-         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
-         zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft
-
-         zchk_vmin = glob_min(v_i)
-         zchk_amax = glob_max(SUM(a_i,dim=3))
-         zchk_amin = glob_min(a_i)
-
-         IF(lwp) THEN
-            IF ( ABS( zchk_v_i   ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (limupdate2) = ',(zchk_v_i * rday)
-            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate2) = ',(zchk_smv * rday)
-            IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limupdate2) = ',(zchk_e_i)
-            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate2) = ',(zchk_vmin * 1.e-3)
-            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate2) = ',zchk_amax
-            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate2) = ',zchk_amin
-         ENDIF
-       ENDIF
-      !- check conservation (C Rousset)
-      ! -------------------------------
+      ! heat content variation (W.m-2)
+      DO jj = 1, jpj
+         DO ji = 1, jpi            
+            diag_heat_dhc(ji,jj) = ( SUM( d_e_i_trp(ji,jj,1:nlay_i,:) + d_e_i_thd(ji,jj,1:nlay_i,:) ) +  & 
+               &                     SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) ) ) * unit_fac * r1_rdtice / area(ji,jj)   
+         END DO
+      END DO
+
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       IF(ln_ctl) THEN   ! Control print
@@ -370,4 +299,5 @@
    
       IF( nn_timing == 1 )  CALL timing_stop('limupdate2')
+
    END SUBROUTINE lim_update2
 #else
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 4649)
@@ -102,5 +102,4 @@
       CALL iom_put( "isst"         , sst_m               )       ! sea surface temperature
       CALL iom_put( "isss"         , sss_m               )       ! sea surface salinity
-      CALL iom_put( "qt_oce"       , qns + qsr           )       ! total flux at ocean surface
       !
       DO jj = 2 , jpjm1
@@ -120,25 +119,6 @@
       CALL iom_put( "vice_ipa"     , z2db                )       ! ice velocity v component
       CALL iom_put( "icevel"       , z2d                 )       ! ice velocity module
-!!SF BE CAREFUL : qsr_oce qnd qns_oce are after penetration over ice
-      CALL iom_put( "qsr_oce"      , qsr                 )       ! solar flux at ocean surface
-      CALL iom_put( "qns_oce"      , qns                 )       ! non-solar flux at ocean surface
-!!SF end be careful
       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
-!!SF commented because this computation is not ok
- !SF because qsr is not qsr_ocean but it contains already qsr_ice
-      !SF
-      !SF DO jj = 1 , jpj
-      !SF    DO ji = 1 , jpi
-      !SF       z2d(ji,jj)  = ( 1._wp - at_i(ji,jj) ) * qsr(ji,jj)
-      !SF   END DO
-      !SF END DO
-      !SF CALL iom_put( "qsr_io"       , z2d                 )        ! solar flux at ice/ocean surface
-      !SF DO jj = 1 , jpj
-      !SF    DO ji = 1 , jpi
-      !SF       z2d(ji,jj)  = ( 1._wp - at_i(ji,jj) ) * qns(ji,jj)
-      !SF   END DO
-      !SF END DO
-      !SF CALL iom_put( "qns_io"       , z2d                 )        ! non-solar flux at ice/ocean surface
       CALL iom_put( "snowpre"      , sprecip             )        ! snow precipitation 
       CALL iom_put( "micesalt"     , smt_i               )        ! mean ice salinity
@@ -211,9 +191,5 @@
       CALL iom_put( "vfxsnw"     , wfx_snw * rday / rhoic  )        ! total snw growth/melt 
       CALL iom_put( "vfxsub"     , wfx_sub * rday / rhoic  )        ! sublimation (snow) 
-
-       CALL iom_put ('hfxdhc1', diag_heat_dhc1(:,:) )          ! Heat content variation in snow and ice 
-       CALL iom_put ('hfxspr', hfx_spr(:,:) )          ! Heat content of snow precip 
-       CALL iom_put ('hfxqsr', qsr(:,:) )          !  solar fluxes used by snw/ice
-       CALL iom_put ('hfxqns', qns(:,:) )          !  non solar fluxes used by snw/ice
+      CALL iom_put( "vfxspr"     , wfx_spr * rday / rhoic  )        ! precip (snow) 
 
        CALL iom_put ('hfxthd', hfx_thd(:,:) )   !  
@@ -222,9 +198,17 @@
        CALL iom_put ('hfxout', hfx_out(:,:) )   !  
        CALL iom_put ('hfxin' , hfx_in(:,:) )   !  
-       CALL iom_put ('hfxtot', hfx_tot(:,:) )   !  
        CALL iom_put ('hfxsnw', hfx_snw(:,:) )   !  
        CALL iom_put ('hfxsub', hfx_sub(:,:) )   !  
        CALL iom_put ('hfxerr', hfx_err(:,:) )   !  
        CALL iom_put ('hfxerr_rem', hfx_err_rem(:,:) )   !  
+
+       CALL iom_put ('hfxsum', hfx_sum(:,:) )   !  
+       CALL iom_put ('hfxbom', hfx_bom(:,:) )   !  
+       CALL iom_put ('hfxbog', hfx_bog(:,:) )   !  
+       CALL iom_put ('hfxdif', hfx_dif(:,:) )   !  
+       CALL iom_put ('hfxopw', hfx_opw(:,:) )   !  
+       CALL iom_put ('hfxtur', fhtur(:,:) * at_i(:,:) )   ! turbulent heat flux at ice base 
+       CALL iom_put ('hfxdhc', diag_heat_dhc(:,:) )          ! Heat content variation in snow and ice 
+       CALL iom_put ('hfxspr', hfx_spr(:,:) )          ! Heat content of snow precip 
 
       !--------------------------------
@@ -327,5 +311,4 @@
       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, "iisfxthd", "Salt flux from thermo"   , ""      , 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 )
@@ -354,5 +337,4 @@
       CALL histwrite( kid, "iicebome", kt, wfx_bom        , jpi*jpj, (/1/) )
       CALL histwrite( kid, "iicesume", kt, wfx_sum        , jpi*jpj, (/1/) )
-      !CALL histwrite( kid, "iisfxthd", kt, sfx_thd        , jpi*jpj, (/1/) )
       CALL histwrite( kid, "iisfxdyn", kt, sfx_dyn        , jpi*jpj, (/1/) )
       CALL histwrite( kid, "iisfxres", kt, sfx_res        , jpi*jpj, (/1/) )
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 4649)
@@ -57,12 +57,20 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_bo_b        !: <==> the 2D  t_bo
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_tot_1d
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_thd_1d
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_spr_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_sum_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_bom_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_bog_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_dif_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_opw_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_snw_1d
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_sub_1d
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_res_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_rem_1d
+
+   ! heat flux associated with ice-atmosphere mass exchange
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_sub_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_spr_1d
+
+   ! heat flux associated with ice-ocean mass exchange
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_thd_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_res_1d
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_ice_1d    !: <==> the 2D  wfx_ice
@@ -75,4 +83,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sni_1d    !: <==> the 2D  wfx_ice
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_opw_1d    !: <==> the 2D  wfx_ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_res_1d    !: <==> the 2D  wfx_ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_spr_1d    !: <==> the 2D  wfx_ice
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_bri_1d    !: <==> the 2D sfx_bri
@@ -82,8 +92,9 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_sni_1d    !: 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_opw_1d    !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_res_1d    !:
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d    !: <==> the 2D  sprecip
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   frld_1d       !: <==> the 2D  frld
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_b        !: <==> the 2D  frld
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_b        !: <==> the 2D  at_i
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhtur_1d       !: <==> the 2D  fhtur
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhld_1d       !: <==> the 2D  fhld
@@ -150,15 +161,16 @@
          &      qsr_ice_1d (jpij) ,     &
          &      fr1_i0_1d(jpij) , fr2_i0_1d(jpij) , qns_ice_1d(jpij) ,     &
-         &      t_bo_b   (jpij) , iatte_1d   (jpij) ,     &
-         &      oatte_1d (jpij) , hfx_tot_1d(jpij), hfx_thd_1d(jpij) , hfx_spr_1d(jpij) , &
-         &      hfx_snw_1d(jpij), hfx_sub_1d(jpij), hfx_err_1d(jpij) , hfx_res_1d(jpij) , hfx_err_rem_1d(jpij),       STAT=ierr(1) )
+         &      t_bo_b   (jpij) , iatte_1d  (jpij) , oatte_1d (jpij) ,     &
+         &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,hfx_dif_1d(jpij) ,hfx_opw_1d(jpij) , &
+         &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) , &
+         &      hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) , hfx_res_1d(jpij) , hfx_err_rem_1d(jpij),       STAT=ierr(1) )
       !
       ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_b     (jpij) ,     &
-         &      fhtur_1d   (jpij) , wfx_ice_1d (jpij) , wfx_snw_1d (jpij) ,     &
-         &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d(jpij) , wfx_bom_1d(jpij) , wfx_sum_1d(jpij) , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) ,  &
+         &      fhtur_1d   (jpij) , wfx_ice_1d (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) ,     &
+         &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d(jpij) , wfx_bom_1d(jpij) , wfx_sum_1d(jpij) , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) ,  wfx_res_1d (jpij) ,  &
          &      dqns_ice_1d(jpij) , qla_ice_1d (jpij) , dqla_ice_1d(jpij) ,     &
          &      tatm_ice_1d(jpij) ,      &   
          &      i0         (jpij) ,     &  
-         &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) ,sfx_sum_1d (jpij) ,sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , &
+         &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) ,sfx_sum_1d (jpij) ,sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , &
          &      dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) ,     &     
          &      dsm_i_si_1d(jpij) , hicol_b    (jpij)                     , STAT=ierr(2) )
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4649)
@@ -69,4 +69,5 @@
 
    PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90
+   PUBLIC lim_prt_state
    
    !! * Substitutions
@@ -170,5 +171,5 @@
          !
          IF( ln_nicep ) THEN      ! control print at a given point
-            jiindx = 3    ;   jjindx =  49
+            jiindx = 15    ;   jjindx =  44
             IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx
          ENDIF
@@ -183,7 +184,8 @@
          u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point
          v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean)
-         !
-         t_bo(:,:) = tfreez( sss_m ) +  rt0          ! masked sea surface freezing temperature [Kelvin]
-         !                                           ! (set to rt0 over land)
+
+         ! masked sea surface freezing temperature [Kelvin]
+         t_bo(:,:) = ( tfreez( sss_m ) +  rt0 ) * tmask(:,:,1) + rt0 * ( 1. - tmask(:,:,1) )
+
          CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os )  ! ... ice albedo
 
@@ -310,27 +312,35 @@
          ! salt, heat and mass fluxes
          sfx    (:,:) = 0._wp   ;
-         sfx_bri(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp   ;   sfx_res(:,:) = 0._wp
+         sfx_bri(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
          sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp
          sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
          sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
-
-         hfx_thd(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp   ;   hfx_snw(:,:) = 0._wp
-         hfx_tot(:,:) = 0._wp   ;   hfx_spr(:,:) = 0._wp   ;   hfx_res(:,:) = 0._wp
-         hfx_sub(:,:) = 0._wp   ;   hfx_err(:,:) = 0._wp   ;   hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp
-         hfx_err_rem(:,:) = 0._wp
-
-         wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
+         sfx_res(:,:) = 0._wp
+
+         wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
          wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp
          wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp
          wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp
-         wfx_res(:,:) = 0._wp   ;   
-         !
-         fhld  (:,:) = 0._wp 
-         fmmflx(:,:) = 0._wp     
-         ftr_ice(:,:,:) = 0._wp   ! part of solar radiation transmitted through the ice
+         wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
+         wfx_spr(:,:) = 0._wp   ;   
+
+         hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp
+         hfx_thd(:,:) = 0._wp   ;   
+         hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
+         hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp
+         hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp
+         hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp
+         hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp 
+         hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp
+
+         !
+         fhld  (:,:)    = 0._wp 
+         fmmflx(:,:)    = 0._wp     
+         ! part of solar radiation transmitted through the ice
+         ftr_ice(:,:,:) = 0._wp
 
          ! diags
-         diag_trp_vi(:,:) = 0._wp  ; diag_trp_vs(:,:) = 0._wp  ;  diag_trp_ei(:,:) = 0._wp  ;  diag_trp_es(:,:) = 0._wp  ; 
-         diag_heat_dhc1(:,:) = 0._wp   ;   
+         diag_trp_vi  (:,:) = 0._wp  ; diag_trp_vs(:,:) = 0._wp  ;  diag_trp_ei(:,:) = 0._wp  ;  diag_trp_es(:,:) = 0._wp
+         diag_heat_dhc(:,:) = 0._wp  
 
          ! dynamical invariants
@@ -820,6 +830,5 @@
                WRITE(numout,*) ' hfx_in       : ', hfx_in(ji,jj)
                WRITE(numout,*) ' hfx_out      : ', hfx_out(ji,jj)
-               WRITE(numout,*) ' hfx_tot      : ', hfx_tot(ji,jj)
-               WRITE(numout,*) ' dhc          : ', diag_heat_dhc1(ji,jj)              
+               WRITE(numout,*) ' dhc          : ', diag_heat_dhc(ji,jj)              
                WRITE(numout,*)
                WRITE(numout,*) ' hfx_dyn      : ', hfx_dyn(ji,jj)
Index: /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4648)
+++ /branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4649)
@@ -391,5 +391,5 @@
          CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux 
          CALL iom_put( "qns"   , qns        )                   ! solar heat flux
-         CALL iom_put( "qsr"   ,       qsr  )                   ! solar heat flux
+         CALL iom_put( "qsr"   ,        qsr )                   ! solar heat flux
          IF( nn_ice > 0 )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction 
       ENDIF
