Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/field_def_nemo-lim.xml
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/field_def_nemo-lim.xml	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/field_def_nemo-lim.xml	(revision 8239)
@@ -26,16 +26,107 @@
          <field id="v_imasstr"    long_name="Sea-ice mass transport along j-axis"                          standard_name="sea_ice_y_transport"                                unit="kg/s"         />
 
+         <!-- SIMIP LIM fields -->
+         <field id="icethic"      long_name="Sea-ice thickness"                                    standard_name="sea_ice_thickness"                          unit="m"            />
+         <field id="icevel"       long_name="Sea-ice speed"                                        standard_name="sea_ice_speed"                              unit="m/s"          />
+
+         <field id="icepres"      long_name="Fraction of time steps with sea ice"                  standard_name="sea_ice_time_fraction"                      unit=""             />
+         <field id="iceconc"      long_name="Sea-ice area fraction"                                standard_name="sea_ice_area_fraction"                      unit=""             />
+         <field id="iceconc_pct"  long_name="Sea-ice area fraction in percent"                     standard_name="sea_ice_area_fraction_in_percent"           unit="%"             > iceconc * 100. </field >
+         <field id="icemass"      long_name="Sea-ice mass per area"                                standard_name="sea_ice_amount"                             unit="kg/m2"        />
+         <field id="icevolu"      long_name="Sea-ice volume per area"                              standard_name="sea_ice_volume"                             unit="m"            />
+         <field id="snomass"      long_name="Snow mass per area"                                   standard_name="liquid_water_content_of_surface_snow"       unit="kg/m2"        />
+         <field id="snowthic"     long_name="Snow thickness"                                       standard_name="surface_snow_thickness"                     unit="m"            />
+
+         <field id="iceconc_cat"  long_name="Sea-ice area fractions in thickness categories"       standard_name="sea_ice_area_fraction_over_categories"      unit=""   grid_ref="grid_T_3D_ncatice" />
+         <field id="icethic_cat"  long_name="Sea-ice thickness in thickness categories"            standard_name="sea_ice_thickness_over_categories"          unit="m"  grid_ref="grid_T_3D_ncatice" />
+         <field id="snowthic_cat" long_name="Snow thickness in thickness categories"               standard_name="snow_thickness_over_categories"             unit="m"  grid_ref="grid_T_3D_ncatice" />
+
+         <field id="icestK"       long_name="ice surface temperature"                              standard_name="sea_ice_surface_temperature        "        unit="C"     > icest + 273.15  </field >
+         <field id="icesntK"      long_name="Temperature at snow-ice interface"                    standard_name="sea_ice_snow_interface_temperature"         unit="degC"  > icesnt + 273.15 </field >
+         <field id="icebotK"      long_name="Temperature at ice-ocean interface"                   standard_name="sea_ice_bottom_temperature"                 unit="degC"  > icebot + 273.15 </field >
+         <field id="icebotK"      long_name="Temperature at ice-ocean interface"                   standard_name="sea_ice_bottom_temperature"                 unit="degC"  > icebot + 273.15 </field >
+
+         <field id="miceage"      long_name="Age of sea ice"                                       standard_name="age_of_sea_ice"                             unit="s"            />
+         <field id="alb_ice"      long_name="Sea-ice or snow albedo"                               standard_name="sea_ice_albedo"                             unit=""             />
+         <field id="icesmass"     long_name="Mass of salt in sea ice per area"                     standard_name="sea_ice_salt_mass"                          unit="kg/m2"        />
+
+         <field id="micesalt"     long_name="Sea ice salinity"                                     standard_name="sea_ice_salinity"                           unit="g/kg"         />
+         <field id="icefb"        long_name="Sea-ice freeboard"                                    standard_name="sea_ice_freeboard"                          unit="m"            />
+
+         <field id="icehcneg"     long_name="Sea-ice heat content per unit area"  standard_name="integral_of_sea_ice_temperature_wrt_depth_expressed_as_heat_content" unit="J/m2" > (-1.0)*icehc </field>
+         <field id="isnhcneg"     long_name="Snow-heat content per unit area"     standard_name="thermal_energy_content_of_surface_snow"                              unit="J/m2" > (-1.0)*isnowhc      </field>
+
+         <field id="wfxsum"       long_name="Freshwater flux from sea-ice surface"               standard_name="freshwater_flux_from_ice_surface"                   unit="kg/m2/s"      />
+
+         <field id="afxthd"       long_name="sea-ice area fraction change from thermodynamics"   standard_name="tendency_of_sea_ice_area_fraction_due_to_dynamics"  unit="s-1"   />
+         <field id="afxdyn"       long_name="sea-ice area fraction change from dynamics"         standard_name="tendency_of_sea_ice_area_fraction_due_to_dynamics"  unit="s-1"   />
+
+         <field id="dmithd"       long_name="sea-ice mass change from thermodynamics"            standard_name="tendency_of_sea_ice_amount_due_to_thermodynamics"   unit="kg/m2/s" />
+         <field id="dmidyn"       long_name="sea-ice mass change from dynamics"                  standard_name="tendency_of_sea_ice_amount_due_to_dynamics"         unit="kg/m2/s" />
+         <field id="dmiopw"       long_name="sea-ice mass change through growth in supercooled open water (aka frazil)"   standard_name="tendency_of_sea_ice_amount_due_to_freezing_in_open_water"   unit="kg/m2/s" />
+         <field id="dmibog"       long_name="sea-ice mass change through basal growth"           standard_name="tendency_of_sea_ice_amount_due_to_congelation_ice_accumulation"  unit="kg/m2/s" />
+         <field id="dmisni"       long_name="sea-ice mass change through snow-to-ice conversion" standard_name="tendency_of_sea_ice_amount_due_to_snow_conversion"               unit="kg/m2/s" />
+         <field id="dmtsub"       long_name="snow and sea ice mass change through sublimation"   standard_name="tendency_of_snow_and_ice_amount_due_to_sublimation"              unit="kg/m2/s" />
+         <field id="dmssub"       long_name="sea-ice mass change through evaporation and sublimation" standard_name="water_evaporation_flux"                                     unit="kg/m2/s" />
+         <field id="dmisub"       long_name="snow mass change through evaporation or sublimation"     standard_name="surface_snow_sublimation_flux"           unit="kg/m2/s" />
+         <field id="dmisum"       long_name="sea-ice mass change through surface melting"        standard_name="tendency_of_sea_ice_amount_due_to_surface_melting"               unit="kg/m2/s" />
+         <field id="dmibom"       long_name="sea-ice mass change through bottom melting"         standard_name="tendency_of_sea_ice_amount_due_to_basal_melting"                 unit="kg/m2/s" />
+
+         <field id="dmsspr"       long_name="snow mass change through snow fall"                 standard_name="snowfall_flux"                                                   unit="kg/m2/s" />
+         <field id="dmsmel"       long_name="snow mass change through melt"                      standard_name="surface_snow_melt_flux"                                          unit="kg/m2/s" />
+         <field id="dmsdyn"       long_name="snow mass change through advection by sea-ice dynamics" standard_name="tendency_of_snow_mass_due_to_sea_ice_dynamics"               unit="kg/m2/s" />
+         <field id="dmsssi"       long_name="snow mass change through snow-to-ice conversion"    standard_name="tendency_of_snow_mass_due_to_snow_to_ice_conversion"             unit="kg/m2/s" />
+
+         <field id="hfxturneg"    long_name="Net sensible heat flux under sea ice"               standard_name="ice_ocean_heat_flux"                          unit="W/m2" > (-1.0)*hfxtur </field>
+         <field id="hfxconsu"     long_name="Net conductive heat flux in ice at the surface"     standard_name="conductive_heat_flux_at_sea_ice_surface"      unit="W/m2" /> 
+         <field id="hfxconbo"     long_name="Net conductive heat fluxes in ice at the bottom"    standard_name="conductive_heat_flux_at_sea_ice_bottom"       unit="W/m2" />
+
+         <field id="msfx"         long_name="Salt flux from sea ice"                             standard_name="salt_flux_from_ice"                           unit="kg/m2/s" > sfx/86400000. </field>
+         <field id="wfxtot"       long_name="Freshwater flux from sea ice"                       standard_name="freshwater_flux_from_ice"                     unit="kg/m2/s" /> 
+
+         <field id="utau_ice"     long_name="X-component of atmospheric stress on sea ice"       standard_name="surface_downward_x_stress"                    unit="N/m2"    />
+         <field id="vtau_ice"     long_name="Y-component of atmospheric stress on sea ice"       standard_name="surface_downward_y_stress"                    unit="N/m2"    />
+
+         <field id="utau_oi"      long_name="X-component of ocean stress on sea ice"             standard_name="sea_ice_base_upward_x_stress"                 unit="N/m2"    />
+         <field id="vtau_oi"      long_name="Y-component of ocean stress on sea ice"             standard_name="sea_ice_base_upward_y_stress"                 unit="N/m2"    />
+
+         <field id="icestr"       long_name="Compressive sea ice strength"                       standard_name="compressive_strength_of_sea_ice"              unit="N/m"     />
+
+         <field id="dssh_dx"      long_name="Sea-surface tilt term in force balance (x-component)"  standard_name="sea_surface_tilt_force_on_sea_ice_x"       unit="N/m2"    />
+         <field id="dssh_dy"      long_name="Sea-surface tilt term in force balance (y-component)"  standard_name="sea_surface_tilt_force_on_sea_ice_y"       unit="N/m2"    />
+
+         <field id="corstrx"      long_name="Coriolis force term in force balance (x-component)"    standard_name="coriolis_force_on_sea_ice_x"               unit="N/m2"    />
+         <field id="corstry"      long_name="Coriolis force term in force balance (y-component)"    standard_name="coriolis_force_on_sea_ice_y"               unit="N/m2"    />
+
+         <field id="intstrx"      long_name="Internal stress term in force balance (x-component)"   standard_name="internal_stress_in_sea_ice_x"              unit="N/m2"    />
+         <field id="intstry"      long_name="Internal stress term in force balance (y-component)"   standard_name="internal_stress_in_sea_ice_y"              unit="N/m2"    />
+
+         <field id="xmtrpice"     long_name="X-component of ice mass transport"                  standard_name="ice_x_transport"                              unit="kg/s" />
+         <field id="ymtrpice"     long_name="Y-component of ice mass transport"                  standard_name="ice_y_transport"                              unit="kg/s" />
+
+         <field id="xmtrpsnw"     long_name="X-component of snw mass transport"                  standard_name="snw_x_transport"                              unit="kg/s" />
+         <field id="ymtrpsnw"     long_name="Y-component of snw mass transport"                  standard_name="snw_y_transport"                              unit="kg/s" />
+
+         <field id="xatrp"        long_name="X-component of ice area transport"                  standard_name="area_x_transport"                             unit="m2/s" />
+         <field id="yatrp"        long_name="Y-component of ice area transport"                  standard_name="area_y_transport"                             unit="m2/s" />
+
+         <field id="xmtrptot"     long_name="X-component of sea-ice mass transport"              standard_name="sea_ice_x_transport"                          unit="kg/s" > xmtrpice + xmtrpsnw </field>
+         <field id="ymtrptot"     long_name="Y-component of sea-ice mass transport"              standard_name="sea_ice_y_transport"                          unit="kg/s" > ymtrpice + ymtrpsnw </field>
+
+         <field id="normstr"      long_name="Average normal stress in sea ice"                   standard_name="average_normal_stress"                        unit="N/m"     />
+         <field id="sheastr"      long_name="Maximum shear stress in sea ice"                    standard_name="maximum_shear_stress"                         unit="N/m"     />
+
+         <field id="idive"        long_name="Divergence of the sea-ice velocity field"           standard_name="divergence_of_sea_ice_velocity"               unit="s-1"     />
+         <field id="ishear"       long_name="Maximum shear of sea-ice velocity field"            standard_name="maximum_shear_of_sea_ice_velocity"            unit="s-1"     />
+
+
+
          <!-- LIM3 fields -->  
+
          <field id="ice_cover"    long_name="Ice fraction"                                                 standard_name="sea_ice_area_fraction"                              unit="1"            />
- 
-         <field id="snowthic_cea" long_name="Snow thickness (cell average)"                                standard_name="surface_snow_thickness"                             unit="m"            />
-         <field id="icethic_cea"  long_name="Ice thickness (cell average)"                                 standard_name="sea_ice_thickness"                                  unit="m"            />
-
-         <field id="uice_ipa"     long_name="Ice velocity along i-axis at I-point (ice presence average)"  standard_name="sea_ice_x_velocity"                                 unit="m/s"          />      
-         <field id="vice_ipa"     long_name="Ice velocity along j-axis at I-point (ice presence average)"  standard_name="sea_ice_y_velocity"                                 unit="m/s"          />               
-         <field id="utau_ice"     long_name="Wind stress along i-axis over the ice at i-point"             standard_name="surface_downward_x_stress"                          unit="N/m2"         />
-         <field id="vtau_ice"     long_name="Wind stress along j-axis over the ice at i-point"             standard_name="surface_downward_y_stress"                          unit="N/m2"         />
-         
-         <field id="iceconc"      long_name="ice concentration"                                            standard_name="sea_ice_area_fraction"                              unit="%"            />
+         <field id="icepres"      long_name="Ice presence"                                                                                                                    unit=""             />
+
+         <field id="uice_ipa"     long_name="X-component of sea ice velocity"                              standard_name="sea_ice_x_velocity"                                 unit="m/s"          /> 
+         <field id="vice_ipa"     long_name="Y-component of sea ice velocity"                              standard_name="sea_ice_y_velocity"                                 unit="m/s"          />      
        	 <field id="isst"         long_name="sea surface temperature"                                      standard_name="sea_surface_temperature"                            unit="degC"         />
          <field id="isss"         long_name="sea surface salinity"                                         standard_name="sea_surface_salinity"                               unit="1e-3"         /> 
@@ -48,7 +139,4 @@
          <field id="qtr_ice"      long_name="solar heat flux transmitted through ice: sum over categories"                                                                    unit="W/m2"         />
          <field id="qemp_ice"     long_name="Downward Heat Flux from E-P over ice"                                                                                            unit="W/m2"         />
-         <field id="micesalt"     long_name="Mean ice salinity"                                                                                                               unit="1e-3"         />
-         <field id="miceage"      long_name="Mean ice age"                                                                                                                    unit="years"        />
-         <field id="alb_ice"      long_name="Mean albedo over sea ice"                                                                                                        unit=""             />
          <field id="albedo"       long_name="Mean albedo over sea ice and ocean"                                                                                              unit=""             />
 
@@ -57,7 +145,4 @@
 
          <field id="iceage_cat"   long_name="Ice age for categories"                                       unit="days"   grid_ref="grid_T_3D_ncatice" />
-         <field id="iceconc_cat"  long_name="Ice concentration for categories"                             unit="%"      grid_ref="grid_T_3D_ncatice" />
-         <field id="icethic_cat"  long_name="Ice thickness for categories"                                 unit="m"      grid_ref="grid_T_3D_ncatice" />
-         <field id="snowthic_cat" long_name="Snow thicknessi for categories"                               unit="m"      grid_ref="grid_T_3D_ncatice" />
          <field id="salinity_cat" long_name="Sea-Ice Bulk salinity for categories"                         unit="g/kg"   grid_ref="grid_T_3D_ncatice" />
          <field id="brinevol_cat" long_name="Brine volume for categories"                                  unit="%"      grid_ref="grid_T_3D_ncatice" />
@@ -68,14 +153,11 @@
 
          <field id="micet"        long_name="Mean ice temperature"                                         unit="degC"     />
-         <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="icehc"        long_name="ice total heat content"                                       unit="J/m2"     /> 
+         <field id="isnowhc"      long_name="snow total heat content"                                      unit="J/m2"     />
          <field id="icest"        long_name="ice surface temperature"                                      unit="degC"     />
+         <field id="icesnt"       long_name="Temperature at snow-ice interface"                            unit="degC"     />
+         <field id="icebot"       long_name="Temperature at ice-ocean interface"                           unit="degC"     />
          <field id="ibrinv"       long_name="brine volume"                                                 unit="%"        />
          <field id="icecolf"      long_name="frazil ice collection thickness"                              unit="m"        />
-         <field id="icestr"       long_name="ice strength"                                                 unit="N/m"      />
-         <field id="icevel"       long_name="ice velocity"                                                 unit="m/s"      />
-         <field id="idive"        long_name="divergence"                                                   unit="1e-8s-1"  />
-         <field id="ishear"       long_name="shear"                                                        unit="1e-8s-1"  />
-         <field id="icevolu"      long_name="ice volume"                                                   unit="m"        />
          <field id="snowvol"      long_name="snow volume"                                                  unit="m"        />
          <field id="tau_icebfr"   long_name="ice friction on ocean bottom for landfast ice"                unit="N/2"      />
@@ -97,5 +179,5 @@
          <field id="sfxopw"       long_name="salt flux from open water ice formation"                      unit="1e-3*kg/m2/day" />
          <field id="sfxsub"       long_name="salt flux from sublimation"                                   unit="1e-3*kg/m2/day" />
-         <field id="sfx"          long_name="salt flux total"                                              unit="1e-3*kg/m2/day" />
+         <field id="sfx"          long_name="Salt flux from sea ice"                                       unit="1e-3*kg/m2/day" />
 
          <field id="vfxbog"       long_name="daily bottom thermo ice prod."                                unit="m/day"   />
@@ -112,9 +194,7 @@
          <field id="vfxsub_err"   long_name="excess of snw sublimation sent to ocean"                      unit="m/day"   />
          <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="m/day"   />
-         <field id="vfxthin"      long_name="daily thermo ice prod. for thin ice(20cm) + open water"      unit="m/day"   />
-
-         <field id="afxtot"       long_name="area tendency (total)"                                        unit="day-1"   />
-         <field id="afxdyn"       long_name="area tendency (dynamics)"                                     unit="day-1"   />
-         <field id="afxthd"       long_name="area tendency (thermo)"                                       unit="day-1"   />
+         <field id="vfxthin"      long_name="daily thermo ice prod. for thin ice(20cm) + open water"       unit="m/day"   />
+
+         <field id="afxtot"       long_name="area tendency (total)"                                        unit="s-1"   />
 
          <field id="hfxsum"       long_name="heat fluxes causing surface ice melt"                         unit="W/m2"  />
@@ -140,5 +220,5 @@
          <field id="hfxtur"       long_name="turbulent heat flux at the ice base"                          unit="W/m2" />
 
-	 <!-- sbcssm variables -->
+	     <!-- sbcssm variables -->
          <field id="sst_m"    unit="degC" />
          <field id="sss_m"    unit="psu"  />
@@ -154,23 +234,35 @@
       <field_group id="SBC_scalar"  grid_ref="grid_T_2D" >
          <!-- available with ln_limdiaout -->
-         <field id="ibgfrcvoltop"    long_name="global mean ice/snow forcing at interface ice/snow-atm (volume equivalent ocean volume)"   unit="km3"       />
-         <field id="ibgfrcvolbot"    long_name="global mean ice/snow forcing at interface ice/snow-ocean (volume equivalent ocean volume)" unit="km3"       />
-         <field id="ibgfrctemtop"    long_name="global mean heat on top of ice/snw/ocean-atm "                                             unit="1e20J"     />
-         <field id="ibgfrctembot"    long_name="global mean heat below ice (on top of ocean) "                                             unit="1e20J"     />
-         <field id="ibgfrcsal"       long_name="global mean ice/snow forcing (salt equivalent ocean volume)"                               unit="pss*km3"   />
-         <field id="ibgfrchfxtop"    long_name="global mean heat flux on top of ice/snw/ocean-atm "                                        unit="W/m2"      />
-         <field id="ibgfrchfxbot"    long_name="global mean heat flux below ice (on top of ocean) "                                        unit="W/m2"      />
+         <field id="ibgfrcvoltop"    long_name="global mean ice/snow forcing at interface ice/snow-atm (volume equivalent ocean volume)"   unit="km3"        />
+         <field id="ibgfrcvolbot"    long_name="global mean ice/snow forcing at interface ice/snow-ocean (volume equivalent ocean volume)" unit="km3"        />
+         <field id="ibgfrctemtop"    long_name="global mean heat on top of ice/snw/ocean-atm "                                             unit="1e20J"      />
+         <field id="ibgfrctembot"    long_name="global mean heat below ice (on top of ocean) "                                             unit="1e20J"      />
+         <field id="ibgfrcsal"       long_name="global mean ice/snow forcing (salt equivalent ocean volume)"                               unit="pss*km3"    />
+         <field id="ibgfrchfxtop"    long_name="global mean heat flux on top of ice/snw/ocean-atm "                                        unit="W/m2"       />
+         <field id="ibgfrchfxbot"    long_name="global mean heat flux below ice (on top of ocean) "                                        unit="W/m2"       />
  
-         <field id="ibgvolume"       long_name="drift in ice/snow volume (equivalent ocean volume)"            unit="km3"        />
-         <field id="ibgsaltco"       long_name="drift in ice salt content (equivalent ocean volume)"           unit="pss*km3"    />
-         <field id="ibgheatco"       long_name="drift in ice/snow heat content"                                unit="1e20J"      />
-         <field id="ibgheatfx"       long_name="drift in ice/snow heat flux"                                   unit="W/m2"       />
-
-         <field id="ibgvol_tot"      long_name="global mean ice volume"                                        unit="km3"        />
-         <field id="sbgvol_tot"      long_name="global mean snow volume"                                       unit="km3"        />
-         <field id="ibgarea_tot"     long_name="global mean ice area"                                          unit="km2"        />
-         <field id="ibgsalt_tot"     long_name="global mean ice salt content"                                  unit="1e-3*km3"   />
-         <field id="ibgheat_tot"     long_name="global mean ice heat content"                                  unit="1e20J"      />
-         <field id="sbgheat_tot"     long_name="global mean snow heat content"                                 unit="1e20J"      />
+         <field id="ibgvolume"       long_name="drift in ice/snow volume (equivalent ocean volume)"                                        unit="km3"        />
+         <field id="ibgsaltco"       long_name="drift in ice salt content (equivalent ocean volume)"                                       unit="pss*km3"    />
+         <field id="ibgheatco"       long_name="drift in ice/snow heat content"                                                            unit="1e20J"      />
+         <field id="ibgheatfx"       long_name="drift in ice/snow heat flux"                                                               unit="W/m2"       />
+
+         <field id="ibgvol_tot"      long_name="global mean ice volume"                                                                    unit="km3"        />
+         <field id="sbgvol_tot"      long_name="global mean snow volume"                                                                   unit="km3"        />
+         <field id="ibgarea_tot"     long_name="global mean ice area"                                                                      unit="km2"        />
+         <field id="ibgsalt_tot"     long_name="global mean ice salt content"                                                              unit="1e-3*km3"   />
+         <field id="ibgheat_tot"     long_name="global mean ice heat content"                                                              unit="1e20J"      />
+         <field id="sbgheat_tot"     long_name="global mean snow heat content"                                                             unit="1e20J"      />
+
+         <field id="NH_iceextt"      long_name="Sea ice extent North"                            standard_name="sea_ice_extent_n"          unit="1e6_km2"    />
+         <field id="SH_iceextt"      long_name="Sea ice extent South"                            standard_name="sea_ice_extent_s"          unit="1e6_km2"    />
+         <field id="NH_icevolu"      long_name="Sea ice volume North"                            standard_name="sea_ice_volume_n"          unit="1e3_km3"    />
+         <field id="SH_icevolu"      long_name="Sea ice volume South"                            standard_name="sea_ice_volume_s"          unit="1e3_km3"    />
+         <field id="NH_icearea"      long_name="Sea ice area North"                              standard_name="sea_ice_area_n"            unit="1e6_km2"    />
+         <field id="SH_icearea"      long_name="Sea ice area South"                              standard_name="sea_ice_area_s"            unit="1e6_km2"    />
+
+         <field id="strait_mifl"     long_name="Sea ice mass flux through straits"      standard_name="sea_ice_mass_transport_across_line" unit="kg/s"  grid_ref="grid_4strait"   />
+         <field id="strait_arfl"     long_name="Sea ice area flux through straits"      standard_name="sea_ice_area_transport_across_line" unit="m2/s"  grid_ref="grid_4strait"   />  
+         <field id="strait_msfl"     long_name="Sea ice snow flux through straits"      standard_name="snow_mass_transport_across_line"    unit="kg/s"  grid_ref="grid_4strait"   />
+
       </field_group>
   
@@ -181,8 +273,162 @@
       
       <field_group id="myvarICE" >
-         <field field_ref="icethic_cea"      name="sithic"     long_name="sea_ice_thickness"        />
+         <field field_ref="icethic"          name="sithic"     long_name="sea_ice_thickness"        />
          <field field_ref="icevolu"          name="sivolu"  />
          <field field_ref="iceconc"          name="siconc"  />
       </field_group>    
-      
+
+      <field_group id="ICE_globalbudget">
+
+    	<field field_ref="ibgvol_tot"     grid_ref="grid_1point"  name="ibgvol_tot"   />
+	    <field field_ref="sbgvol_tot"     grid_ref="grid_1point"  name="sbgvol_tot"   />
+    	<field field_ref="ibgarea_tot"    grid_ref="grid_1point"  name="ibgarea_tot"  />
+    	<field field_ref="ibgsalt_tot"    grid_ref="grid_1point"  name="ibgsalt_tot"  />
+    	<field field_ref="ibgheat_tot"    grid_ref="grid_1point"  name="ibgheat_tot"  />
+    	<field field_ref="sbgheat_tot"    grid_ref="grid_1point"  name="sbgheat_tot"  />
+
+     	<field field_ref="ibgvolume"      grid_ref="grid_1point"  name="ibgvolume"    />
+    	<field field_ref="ibgsaltco"      grid_ref="grid_1point"  name="ibgsaltco"    />
+    	<field field_ref="ibgheatco"      grid_ref="grid_1point"  name="ibgheatco"    />
+        <field field_ref="ibgheatfx"      grid_ref="grid_1point"  name="ibgheatfx"    />
+
+    	<field field_ref="ibgfrcvoltop"   grid_ref="grid_1point"  name="ibgfrcvoltop" />
+    	<field field_ref="ibgfrcvolbot"   grid_ref="grid_1point"  name="ibgfrcvolbot" />
+    	<field field_ref="ibgfrctemtop"   grid_ref="grid_1point"  name="ibgfrctemtop" />
+    	<field field_ref="ibgfrctembot"   grid_ref="grid_1point"  name="ibgfrctembot" />
+    	<field field_ref="ibgfrcsal"      grid_ref="grid_1point"  name="ibgfrcsal"    />
+        <field field_ref="ibgfrchfxtop"   grid_ref="grid_1point"  name="ibgfrchfxtop" />
+        <field field_ref="ibgfrchfxbot"   grid_ref="grid_1point"  name="ibgfrchfxbot" />
+
+      </field_group>
+
+
+      <!--============================-->
+      <!-- SIMIP sea ice field groups -->
+      <!-- Notz et al 2016            -->
+      <!--============================-->
+
+      <!-- SIMIP monthly scalar variables -->
+      <field_group id="SImon_scalars">
+        <!-- Integrated quantities -->
+        <field field_ref="NH_iceextt"     grid_ref="grid_1point"  name="siextentn"    />
+        <field field_ref="SH_iceextt"     grid_ref="grid_1point"  name="siextents"    />
+        <field field_ref="NH_icevolu"     grid_ref="grid_1point"  name="sivoln"       />
+        <field field_ref="SH_icevolu"     grid_ref="grid_1point"  name="sivols"       />
+        <field field_ref="NH_icearea"     grid_ref="grid_1point"  name="siarean"      />
+        <field field_ref="SH_icearea"     grid_ref="grid_1point"  name="siareas"      />
+      </field_group>
+
+      <!-- SIMIP daily fields -->
+      <field_group id="SIday_fields">
+          <field field_ref="iceconc_pct"      name="siconc"       />
+          <field field_ref="snowthic"         name="sisnthick"    />
+          <field field_ref="icevel"           name="sispeed" />
+          <field field_ref="icestK"           name="sitemptop"    />
+          <field field_ref="icethic"          name="sithick"      />
+          <field field_ref="icepres"          name="sitimefrac"   />
+          <field field_ref="uice_ipa"         name="siu"              />
+          <field field_ref="vice_ipa"         name="siv"              />
+      </field_group>
+
+      <!-- SIMIP monthly fields -->
+      <field_group id="SImon_fields">
+          <!-- Sea-ice state variables -->
+          <field field_ref="icepres"          name="sitimefrac"   />
+          <field field_ref="iceconc_pct"      name="siconc"       />
+          <field field_ref="icemass"          name="simass"       />
+          <field field_ref="icevolu"          name="sivol"        />
+          <field field_ref="snomass"          name="sisnmass"     />
+          <field field_ref="snowthic"         name="sisnthick"    />
+          <field field_ref="icethic"          name="sithick"      />
+          <field field_ref="icevel"           name="sispeed" />
+
+          <!-- category-dependent fields -->
+          <field field_ref="iceconc_cat"      name="siitdconc"    />
+          <field field_ref="icethic_cat"      name="siitdthick"   />
+          <field field_ref="snowthic_cat"     name="siitdsnthick" />
+
+          <!-- additional state variables-->
+          <field field_ref="icestK"           name="sitemptop"    />
+          <field field_ref="icesntK"          name="sitempsnic"   />
+          <field field_ref="icebotK"          name="sitempbot"    />
+          <field field_ref="miceage"          name="siage"        />
+          <field field_ref="alb_ice"          name="sialb"        />
+          <field field_ref="icesmass"         name="sisaltmass"   />
+          <field field_ref="micesalt"         name="sisali"       />
+          <field field_ref="icefb"            name="sifb"         />
+          <field field_ref="icehcneg"         name="sihc"         />
+          <field field_ref="isnhcneg"         name="sisnhc"       />
+
+          <!-- freshwater fluxes -->
+          <field field_ref="wfxtot"           name="siflfwbot"    />
+          <field field_ref="wfxsum"           name="siflfwdrain"  />
+
+          <!-- area balance -->
+          <field field_ref="afxthd"           name="sidconcth"    />
+          <field field_ref="afxdyn"           name="sidconcdyn"   />
+
+          <!-- mass balance -->
+          <field field_ref="dmithd"           name="sidmassth"    />
+          <field field_ref="dmidyn"           name="sidmassdyn"   />
+          <field field_ref="dmiopw"           name="sidmassgrowthwat" />
+          <field field_ref="dmibog"           name="sidmassgrowthbot" />
+          <field field_ref="dmisni"           name="sidmasssi"        />
+          <field field_ref="dmisub"           name="sidmassevapsubl"  />
+          <field field_ref="dmisum"           name="sidmassmelttop"   />
+          <field field_ref="dmibom"           name="sidmassmeltbot"   />
+          <field field_ref="dmsspr"           name="sndmasssnf"       />
+          <field field_ref="dmsmel"           name="sndmassmelt"      />
+          <field field_ref="dmssub"           name="sndmasssubl"      />
+          <field field_ref="dmsdyn"           name="sndmassdyn"       />
+          <field field_ref="dmsssi"           name="sndmasssi"        />
+
+          <!-- heat balance -->
+          <field field_ref="hfxturneg"        name="siflsensupbot"    />
+          <field field_ref="hfxconsu"         name="siflcondtop"      />
+          <field field_ref="hfxconbo"         name="siflcondbot"      />
+
+          <!-- salt balance -->
+          <field field_ref="msfx"             name="siflsaltbot"      />
+
+          <!-- sea-ice dynamics -->
+          <field field_ref="uice_ipa"         name="siu"              />
+          <field field_ref="vice_ipa"         name="siv"              />
+
+          <field field_ref="utau_ice"         name="sistrxdtop"       />
+          <field field_ref="vtau_ice"         name="sistrydtop"       />
+
+          <field field_ref="utau_oi"          name="sistrxubot"       />
+          <field field_ref="vtau_oi"          name="sistryubot"       />
+
+          <field field_ref="icestr"           name="sicompstren"      />
+
+          <field field_ref="dssh_dx"          name="siforcetiltx"     />
+          <field field_ref="dssh_dy"          name="siforcetilty"     />
+
+          <field field_ref="corstrx"          name="siforcecoriolx"   />
+          <field field_ref="corstry"          name="siforcecorioly"   />
+
+          <field field_ref="intstrx"          name="siforceintstrx"   />
+          <field field_ref="intstry"          name="siforceintstry"   />
+
+          <field field_ref="xmtrptot"         name="sidmasstranx"     />
+          <field field_ref="ymtrptot"         name="sidmasstrany"     />
+
+          <field field_ref="normstr"          name="sistresave"       operation="instant" />
+          <field field_ref="sheastr"          name="sistremax"        operation="instant" />
+
+          <field field_ref="idive"            name="sidivvel"         operation="instant" />
+          <field field_ref="ishear"           name="sishevel"         operation="instant" />
+
+          <!-- sea ice fluxes across straits: 2D arrays to be post-processed  -->
+          <field field_ref="xmtrpice"         name="simassacrossline_x" />
+          <field field_ref="ymtrpice"         name="simassacrossline_y" />
+          <field field_ref="xmtrpsnw"         name="snmassacrossline_x" />
+          <field field_ref="ymtrpsnw"         name="snmassacrossline_y" />
+          <field field_ref="xatrp"            name="siareaacrossline_x" />
+          <field field_ref="yatrp"            name="siareaacrossline_y" />
+
+      </field_group>
+
+     
    </field_definition>
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8239)
@@ -108,5 +108,5 @@
    !! v_s         |      -      |    Snow volume per unit area    | m     |
    !! smv_i       |      -      |    Sea ice salt content         | ppt.m |
-   !! oa_i        !      -      !    Sea ice areal age content    | day   |
+   !! oa_i        !      -      !    Sea ice areal age content    | s     |
    !! e_i         !      -      !    Ice enthalpy                 | J/m2  | 
    !!      -      ! q_i_1d      !    Ice enthalpy per unit vol.   | J/m3  | 
@@ -123,5 +123,5 @@
    !! sm_i        ! sm_i_1d     |    Sea ice bulk salinity        ! ppt   |
    !! s_i         ! s_i_1d      |    Sea ice salinity profile     ! ppt   |
-   !! o_i         !      -      |    Sea ice Age                  ! days  |
+   !! o_i         !      -      |    Sea ice Age                  ! s     |
    !! t_i         ! t_i_1d      |    Sea ice temperature          ! K     |
    !! t_s         ! t_s_1d      |    Snow temperature             ! K     |
@@ -328,5 +328,10 @@
    ! END MV MP 2016
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr     !: snow precipitation on ice  [kg.m-2.s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub     !: snow/ice sublimation       [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub     !: sublimation of snow/ice    [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sub !: snow sublimation           [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice_sub !: ice sublimation            [kg.m-2.s-1]
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_dyn !: dynamical component of wfx_snw    [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sum !: surface melt component of wfx_snw [kg.m-2.s-1]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice     !: ice-ocean mass exchange                   [kg.m-2.s-1]
@@ -395,6 +400,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   smv_i     !: Sea-Ice Bulk salinity times volume per area (ppt.m)
    !                                                                    !  this is an extensive variable that has to be transported
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i       !: Sea-Ice Age (days)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i      !: Sea-Ice Age times ice area (days)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i       !: Sea-Ice Age (s)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i      !: Sea-Ice Age times ice area (s)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bv_i      !: brine volume
 
@@ -479,4 +484,34 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vice     !: ice volume variation   [m/s] 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vsnw     !: snw volume variation   [m/s] 
+
+   !
+   !!--------------------------------------------------------------------------
+   !! * SIMIP extra diagnostics
+   !!--------------------------------------------------------------------------
+   ! Extra sea ice diagnostics to address the data request
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   t_si          !: Temperature at Snow-ice interface (K) 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   tm_si         !: mean temperature at the snow-ice interface (K) 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dmi_dyn  !: Change in ice mass due to ice dynamics (kg/m2/s)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dms_dyn  !: Change in snow mass due to ice dynamics (kg/m2/s)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_xmtrp_ice !: X-component of ice mass transport (kg/s)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_ymtrp_ice !: Y-component of ice mass transport (kg/s)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_xmtrp_snw !: X-component of snow mass transport (kg/s)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_ymtrp_snw !: Y-component of snow mass transport (kg/s)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_xatrp    !: X-component of area transport (m2/s)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_yatrp    !: Y-component of area transport (m2/s)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_bo    !: Bottom conduction flux (W/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_su    !: Surface conduction flux (W/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_utau_oi  !: X-direction ocean-ice stress
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vtau_oi  !: Y-direction ocean-ice stress  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dssh_dx  !: X-direction sea-surface tilt term (N/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dssh_dy  !: X-direction sea-surface tilt term (N/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_corstrx  !: X-direction coriolis stress (N/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_corstry  !: Y-direction coriolis stress (N/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_intstrx  !: X-direction internal stress (N/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_intstry  !: Y-direction internal stress (N/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_sig1     !: Average normal stress in sea ice   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_sig2     !: Maximum shear stress in sea ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_shear    !: Maximum shear of sea-ice velocity field
+
    !
    !!----------------------------------------------------------------------
@@ -493,5 +528,5 @@
       INTEGER :: ice_alloc
       !
-      INTEGER :: ierr(15), ii
+      INTEGER :: ierr(18), ii
       !!-----------------------------------------------------------------
 
@@ -507,6 +542,7 @@
 
       ii = ii + 1
-      ALLOCATE( 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_lam(jpi,jpj) ,     &
+      ALLOCATE( t_bo   (jpi,jpj) , frld       (jpi,jpj) , pfrld      (jpi,jpj) , phicif     (jpi,jpj) ,  &
+         &      wfx_snw(jpi,jpj) , wfx_snw_dyn(jpi,jpj) , wfx_snw_sum(jpi,jpj) , wfx_snw_sub(jpi,jpj) ,  &
+         &      wfx_ice(jpi,jpj) , wfx_sub    (jpi,jpj) , wfx_ice_sub(jpi,jpj) , wfx_lam    (jpi,jpj) ,  &
          ! MV MP 2016
          &      wfx_pnd(jpi,jpj) , wfx_snw_sum(jpi,jpj) ,                                       &
@@ -568,4 +604,5 @@
 
       ! MV MP 2016
+      ii = ii + 1
       ALLOCATE( sxap(jpi,jpj,jpl) , syap(jpi,jpj,jpl) , sxxap(jpi,jpj,jpl) , syyap(jpi,jpj,jpl) , sxyap(jpi,jpj,jpl) ,   &
          &      sxvp(jpi,jpj,jpl) , syvp(jpi,jpj,jpl) , sxxvp(jpi,jpj,jpl) , syyvp(jpi,jpj,jpl) , sxyvp(jpi,jpj,jpl) ,   &
@@ -590,4 +627,19 @@
          &      diag_trp_es(jpi,jpj) , diag_trp_smv(jpi,jpj) , diag_heat  (jpi,jpj),   &
          &      diag_smvi  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), STAT=ierr(ii) )
+
+      ! * SIMIP diagnostics
+      ii = ii + 1
+      ALLOCATE( t_si (jpi,jpj,jpl)    , tm_si(jpi,jpj)        ,    & 
+                diag_dmi_dyn(jpi,jpj) , diag_dms_dyn(jpi,jpj) ,    &
+                diag_xmtrp_ice(jpi,jpj), diag_ymtrp_ice(jpi,jpj),  &
+                diag_xmtrp_snw(jpi,jpj), diag_ymtrp_snw(jpi,jpj),  &
+                diag_xatrp(jpi,jpj)    , diag_yatrp(jpi,jpj)    ,  &
+                diag_fc_bo(jpi,jpj)   , diag_fc_su(jpi,jpj)   ,    &
+                diag_utau_oi(jpi,jpj) , diag_vtau_oi(jpi,jpj) ,    &
+                diag_dssh_dx(jpi,jpj) , diag_dssh_dy(jpi,jpj) ,    &
+                diag_corstrx(jpi,jpj) , diag_corstry(jpi,jpj) ,    &
+                diag_intstrx(jpi,jpj) , diag_intstry(jpi,jpj) ,    &
+                diag_sig1(jpi,jpj)    , diag_sig2(jpi,jpj)    ,    &
+                STAT = ierr(ii) )
 
       ice_alloc = MAXVAL( ierr(:) )
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 8239)
@@ -302,5 +302,5 @@
                   ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness
                   sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity
-                  o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                                    ! age (1 day)
+                  o_i(ji,jj,jl)   = 0._wp                                                     ! age (0 day)
                   t_su(ji,jj,jl)  = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 8239)
@@ -500,4 +500,5 @@
       INTEGER ::   icells            ! number of cells with a_i > puny
       REAL(wp) ::   hL, hR, farea    ! left and right limits of integration
+      REAL(wp) ::   zwfx_snw         ! snow mass flux increment
 
       INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices
@@ -654,5 +655,5 @@
                   &                            - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice     ! and get  sm_i  from the ocean 
             ENDIF
-               
+
             !------------------------------------------            
             ! 3.7 Put the snow somewhere in the ocean
@@ -663,6 +664,9 @@
             !  During the next time step, thermo_rates will determine whether
             !  the ocean cools or new ice grows.
-            wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
-               &                              + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
+            zwfx_snw           = ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
+               &                 + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
+
+            wfx_snw_dyn(ji,jj)  =   wfx_snw_dyn(ji,jj) + zwfx_snw
+            wfx_snw(ji,jj)      =   wfx_snw(ji,jj)     + zwfx_snw
 
             hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
@@ -789,4 +793,8 @@
       END DO ! jl1 (deforming categories)
 
+      ! SIMIP diagnostics
+      diag_dmi_dyn(:,:) = - wfx_dyn(:,:)     + rhoic * diag_trp_vi(:,:)
+      diag_dms_dyn(:,:) = - wfx_snw_dyn(:,:) + rhosn * diag_trp_vs(:,:)
+      
       !
       CALL wrk_dealloc( jpij,        indxi, indxj )
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 8239)
@@ -117,9 +117,10 @@
       REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass
       REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars
-      REAL(wp) ::   zTauO, zTauB, zTauE, zCor, zvel                          ! temporary scalars
+      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                                ! temporary scalars
 
       REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress
       REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity
       REAL(wp) ::   zintb, zintn                                             ! dummy argument
+      REAL(wp) ::   zfac_x, zfac_y
       
       REAL(wp), POINTER, DIMENSION(:,:) ::   z1_e1t0, z1_e2t0                ! scale factors
@@ -140,4 +141,7 @@
                                                                              !   ocean surface (ssh_m) if ice is not embedded
                                                                              !   ice top surface if ice is embedded   
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zCorx, zCory                    ! Coriolis stress array
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
+
       REAL(wp), POINTER, DIMENSION(:,:) ::   zswitchU, zswitchV              ! dummy arrays
       REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence
@@ -154,4 +158,6 @@
       CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
       CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
+      CALL wrk_alloc( jpi,jpj, zCorx, zCory)
+      CALL wrk_alloc( jpi,jpj, ztaux_oi, ztauy_oi)
 
 #if defined key_agrif 
@@ -427,4 +433,7 @@
                      &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
 
+                  ! Ocean-to-Ice stress
+                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+
                   ! tau_bottom/v_ice
                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
@@ -432,10 +441,10 @@
 
                   ! Coriolis at V-points (energy conserving formulation)
-                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
 
                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
 
                   ! landfast switch => 0 = static friction ; 1 = sliding friction
@@ -451,5 +460,5 @@
                   ! Bouillon 2013
                   !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
-                  !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
+                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
                   !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
 
@@ -471,4 +480,7 @@
                      &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
 
+                  ! Ocean-to-Ice stress
+                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+
                   ! tau_bottom/u_ice
                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
@@ -476,10 +488,10 @@
 
                   ! Coriolis at U-points (energy conserving formulation)
-                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  &
+                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
                   
                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
 
                   ! landfast switch => 0 = static friction ; 1 = sliding friction
@@ -495,5 +507,5 @@
                   ! Bouillon 2013
                   !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
-                  !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
+                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
                   !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
                END DO
@@ -516,4 +528,7 @@
                      &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
 
+                  ! Ocean-to-Ice stress
+                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+
                   ! tau_bottom/u_ice
                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
@@ -521,10 +536,10 @@
 
                   ! Coriolis at U-points (energy conserving formulation)
-                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  &
+                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
 
                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
 
                   ! landfast switch => 0 = static friction ; 1 = sliding friction
@@ -540,5 +555,5 @@
                   ! Bouillon 2013
                   !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
-                  !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
+                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
                   !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
                END DO
@@ -559,4 +574,7 @@
                      &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
 
+                  ! Ocean-to-Ice stress
+                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+
                   ! tau_bottom/v_ice
                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
@@ -564,10 +582,10 @@
                   
                   ! Coriolis at V-points (energy conserving formulation)
-                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
 
                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
 
                   ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
@@ -583,5 +601,5 @@
                   ! Bouillon 2013
                   !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
-                  !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
+                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
                   !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
                END DO
@@ -662,5 +680,54 @@
 
       !------------------------------------------------------------------------------!
-      ! 5) Control prints of residual and charge ellipse
+      ! 5) SIMIP diagnostics
+      !------------------------------------------------------------------------------!
+                           
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+             rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
+
+             ! Stress tensor invariants (normal and shear stress N/m)
+             diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
+             diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
+
+             ! Stress terms of the momentum equation (N/m2)
+             diag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term
+             diag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
+
+             diag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term
+             diag_corstry(ji,jj) = zCory(ji,jj) * rswitch
+
+             diag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term
+             diag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
+           
+             diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
+             diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
+
+             ! 2D ice mass, snow mass, area transport arrays (X, Y)
+             zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
+             zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
+
+             diag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
+             diag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
+
+             diag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
+             diag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
+
+             diag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component
+             diag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   ''
+
+         END DO
+      END DO
+
+      CALL lbc_lnk_multi(   diag_sig1   , 'T',  1., diag_sig2   , 'T',  1.,   &
+                 &          diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1.,   &
+                 &          diag_corstrx, 'U', -1., diag_corstry, 'V', -1.,   & 
+                 &          diag_intstrx, 'U', -1., diag_intstry, 'V', -1.    )
+
+      CALL lbc_lnk_multi(   diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1.    )
+
+      !
+      !------------------------------------------------------------------------------!
+      ! 6) Control prints of residual and charge ellipse
       !------------------------------------------------------------------------------!
       !
@@ -703,4 +770,6 @@
       CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
       CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
+      CALL wrk_dealloc( jpi,jpj, zCorx, zCory )
+      CALL wrk_dealloc( jpi,jpj, ztaux_oi, ztauy_oi )
 
    END SUBROUTINE lim_rhg
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8239)
@@ -83,5 +83,5 @@
       INTEGER  :: ji, jj, jk, jl   ! dummy loop indices
       INTEGER  :: nbpb             ! nb of icy pts for vertical thermo calculations
-      REAL(wp) :: zfric_u, zqld, zqfr
+      REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg
       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       REAL(wp), PARAMETER :: zfric_umin = 0._wp           ! lower bound for the friction velocity (cice value=5.e-04)
@@ -175,15 +175,20 @@
 
             ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
+            ! includes supercooling potential energy (>0) or "above-freezing" energy (<0)
             zqfr = tmask(ji,jj,1) * rau0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) )
 
-            ! --- Energy from the turbulent oceanic heat flux (W/m2) --- !
+            ! --- Above-freezing sensible heat content (J/m2 grid)
+            zqfr_neg = tmask(ji,jj,1) * rau0 * rcp * e3t_m(ji,jj) * MIN( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ), 0._wp )
+
+            ! --- Sensible ocean-to-ice heat flux (W/m2)
             zfric_u      = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 
-            fhtur(ji,jj) = MAX( 0._wp, rswitch * rau0 * rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2
-            fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
+            fhtur(ji,jj) = rswitch * rau0 * rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2
+
+            fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
             ! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach 
             !                        the freezing point, so that we do not have SST < T_freeze
             !                        This implies: - ( fhtur(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0
 
-            !-- Energy Budget of the leads (J.m-2). Must be < 0 to form ice
+            !-- Energy Budget of the leads (J.m-2), source of lateral accretion. Must be < 0 to form ice
             qlead(ji,jj) = MIN( 0._wp , zqld - ( fhtur(ji,jj) * at_i(ji,jj) * rdt_ice ) - zqfr )
 
@@ -470,4 +475,6 @@
          CALL tab_2d_1d( nbpb, wfx_snw_sum_1d(1:nbpb), wfx_snw_sum  , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub         , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, wfx_snw_sub_1d(1:nbpb), wfx_snw_sub  , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, wfx_ice_sub_1d(1:nbpb), wfx_ice_sub  , jpi, jpj, npb(1:nbpb) )
          !
          CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog         , jpi, jpj, npb(1:nbpb) )
@@ -500,4 +507,8 @@
          CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) )
          !
+         ! SIMIP diagnostics
+         CALL tab_2d_1d( nbpb, diag_fc_bo_1d   (1:nbpb), diag_fc_bo  , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, diag_fc_su_1d   (1:nbpb), diag_fc_su  , jpi, jpj, npb(1:nbpb) )
+         !
       CASE( 2 )            ! from 1D to 2D
          !
@@ -522,4 +533,6 @@
          CALL tab_1d_2d( nbpb, wfx_snw_sum   , npb, wfx_snw_sum_1d(1:nbpb),jpi, jpj )
          CALL tab_1d_2d( nbpb, wfx_sub       , npb, wfx_sub_1d(1:nbpb)   , jpi, jpj )
+         CALL tab_1d_2d( nbpb, wfx_snw_sub   , npb, wfx_snw_sub_1d(1:nbpb), jpi, jpj )
+         CALL tab_1d_2d( nbpb, wfx_ice_sub   , npb, wfx_ice_sub_1d(1:nbpb), jpi, jpj )
          !
          CALL tab_1d_2d( nbpb, wfx_bog       , npb, wfx_bog_1d(1:nbpb)   , jpi, jpj )
@@ -554,5 +567,9 @@
          CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
          CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj )
-         !         
+         !
+         ! SIMIP diagnostics         
+         CALL tab_1d_2d( nbpb, t_si(:,:,jl)   , npb, t_si_1d    (1:nbpb)     , jpi, jpj )
+         CALL tab_1d_2d( nbpb, diag_fc_bo     , npb, diag_fc_bo_1d(1:nbpb)   , jpi, jpj )
+         CALL tab_1d_2d( nbpb, diag_fc_su     , npb, diag_fc_su_1d(1:nbpb)   , jpi, jpj )
       END SELECT
       !
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8239)
@@ -87,5 +87,5 @@
       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) ::   zsstK        ! SST in Kelvin
+      REAL(wp) ::   zsstK        ! SST (K)
 
       REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow                       (J.m-3)
@@ -277,5 +277,6 @@
             &                              ) * a_i_1d(ji) * r1_rdtice
          ! Mass flux by sublimation
-         wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
+         wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
+
          ! new snow thickness
          ht_s_1d(ji)    =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) )
@@ -378,5 +379,6 @@
             hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * q_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice
             ! Mass flux > 0
-            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice
+            wfx_ice_sub_1d(ji) =  wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice
+
             ! update remaining mass flux
             zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic
@@ -615,5 +617,10 @@
          IF( ln_limctl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
       END DO
-      
+
+      ! Water fluxes
+      DO ji = kideb, kiut
+         wfx_sub_1d(ji) = wfx_snw_sub_1d(ji) + wfx_ice_sub_1d(ji) ! sum ice and snow sublimation contributions
+      END DO
+
       !
       !------------------------------------------------------------------------------|
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 8239)
@@ -734,4 +734,14 @@
       END DO  ! End of the do while iterative procedure
 
+      ! MV SIMIP 2016
+      !--- Snow-ice interfacial temperature (diagnostic SIMIP)
+      DO ji = kideb, kiut
+         zfac        = 1. / MAX( epsi10 , rn_cdsn * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) )
+         t_si_1d(ji) = ( rn_cdsn        * zh_i(ji) * t_s_1d(ji,1) + &
+            &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) * zfac
+      END DO
+      WHERE( ( zh_s < 1.0e-3 ) ) ; t_si_1d(:) = t_su_1d(:) ; END WHERE
+      ! END MV SIMIP 2016
+
       IF( ln_limctl .AND. lwp ) THEN
          WRITE(numout,*) ' zerritmax : ', zerritmax
@@ -751,4 +761,10 @@
          fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
       END DO
+
+      ! MV SIMIP 2016
+      !--- Conduction fluxes (positive downwards)
+      diag_fc_bo_1d(:)   = diag_fc_bo_1d(:) + fc_bo_i(:) * a_i_1d(:) / at_i_1d(:)
+      diag_fc_su_1d(:)   = diag_fc_su_1d(:) + fc_su(:)   * a_i_1d(:) / at_i_1d(:)
+      ! END MV SIMIP 2016
 
       ! --- computes sea ice energy of melting compulsory for limthd_dh --- !
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 8239)
@@ -80,5 +80,4 @@
                IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN
                   a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
-                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
                ENDIF
             END DO
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 8239)
@@ -77,5 +77,4 @@
             IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN
                a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin
-               oa_i(ji,jj,1) = oa_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin
             ENDIF
          END DO
@@ -95,5 +94,4 @@
                IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN
                   a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
-                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
                ENDIF
             END DO
@@ -155,5 +153,5 @@
       ! -------------------------------------------------
       DO jl  = 1, jpl
-         oa_i(:,:,jl) = oa_i(:,:,jl) + a_i(:,:,jl) * rdt_ice / rday   ! ice natural aging
+         oa_i(:,:,jl) = oa_i(:,:,jl) + a_i(:,:,jl) * rdt_ice          ! ice natural aging
          afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
       END DO
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 8239)
@@ -120,4 +120,5 @@
          tm_i(:,:)  = 0._wp
          tm_su(:,:) = 0._wp
+         tm_si(:,:) = 0._wp
          om_i (:,:) = 0._wp
          DO jl = 1, jpl
@@ -127,5 +128,6 @@
                   rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
                   tm_su(ji,jj) = tm_su(ji,jj) + rswitch * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 )
-                  om_i (ji,jj) = om_i (ji,jj) + rswitch *   oa_i(ji,jj,jl)                         / MAX( at_i(ji,jj) , epsi10 )
+                  tm_si(ji,jj) = tm_si(ji,jj) + rswitch * ( t_si(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 )
+                  om_i (ji,jj) = om_i (ji,jj) + rswitch *   oa_i(ji,jj,jl)                         / MAX( at_i(ji,jj) , epsi06 )
                END DO
             END DO
@@ -145,4 +147,5 @@
          tm_i  = tm_i  + rt0
          tm_su = tm_su + rt0
+         tm_si = tm_si + rt0
          !
       ENDIF
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 8239)
@@ -54,8 +54,18 @@
       !
       INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices
-      REAL(wp) ::  z1_365
-      REAL(wp) ::  z2da, z2db, ztmp
+      REAL(wp) ::  z2da, z2db, ztmp, zrho1, zrho2
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zswi2
       REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, zswi    ! 2D workspace
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zfb          ! ice freeboard
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zamask, zamask15 ! 15% concentration mask
+
+      ! Global ice diagnostics (SIMIP)
+      REAL(wp) ::  zdiag_area_nh, &   ! area, extent, volume
+         &         zdiag_extt_nh, &
+         &         zdiag_area_sh, & 
+         &         zdiag_extt_sh, & 
+         &         zdiag_volu_nh, & 
+         &         zdiag_volu_sh 
+
       !!-------------------------------------------------------------------
 
@@ -64,9 +74,8 @@
       CALL wrk_alloc( jpi,jpj,jpl, zswi2 )
       CALL wrk_alloc( jpi,jpj    , z2d, zswi )
-
+      CALL wrk_alloc( jpi,jpj    , zfb, zamask, zamask15 )
       !-----------------------------
       ! Mean category values
       !-----------------------------
-      z1_365 = 1._wp / 365._wp
 
       ! brine volume
@@ -76,5 +85,7 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zswi(ji,jj)  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) )
+            zswi(ji,jj)      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
+            zamask(ji,jj)    = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.05   ) ) ! 1 if 5% ice, 0 if less - required to mask thickness and snow depth
+            zamask15(ji,jj)  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15   ) ) ! 1 if 15% ice, 0 if less
          END DO
       END DO
@@ -119,7 +130,5 @@
       IF ( iom_use( "tau_icebfr" ) )    CALL iom_put( "tau_icebfr"  , tau_icebfr             )  ! ice friction with ocean bottom (landfast ice)  
       !
-      IF ( iom_use( "miceage" ) )       CALL iom_put( "miceage"     , om_i * zswi * z1_365   )  ! mean ice age
-      IF ( iom_use( "icethic_cea" ) )   CALL iom_put( "icethic_cea" , htm_i * zswi           )  ! ice thickness mean
-      IF ( iom_use( "snowthic_cea" ) )  CALL iom_put( "snowthic_cea", htm_s * zswi           )  ! snow thickness mean
+      IF ( iom_use( "miceage" ) )       CALL iom_put( "miceage"     , om_i * zswi * zamask15 )  ! mean ice age
       IF ( iom_use( "micet" ) )         CALL iom_put( "micet"       , ( tm_i  - rt0 ) * zswi )  ! ice mean    temperature
       IF ( iom_use( "icest" ) )         CALL iom_put( "icest"       , ( tm_su - rt0 ) * zswi )  ! ice surface temperature
@@ -133,12 +142,12 @@
       CALL iom_put( "isnowhc"     , et_s  * zswi        )        ! snow total heat content
       CALL iom_put( "ibrinv"      , bvm_i * zswi * 100. )        ! brine volume
-      CALL iom_put( "utau_ice"    , utau_ice            )        ! wind stress over ice along i-axis at I-point
-      CALL iom_put( "vtau_ice"    , vtau_ice            )        ! wind stress over ice along j-axis at I-point
+      CALL iom_put( "utau_ice"    , utau_ice*zswi       )        ! wind stress over ice along i-axis at I-point
+      CALL iom_put( "vtau_ice"    , vtau_ice*zswi       )        ! wind stress over ice along j-axis at I-point
       CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation 
       CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity
 
-      CALL iom_put( "icestr"      , strength * zswi )    ! ice strength
-      CALL iom_put( "idive"       , divu_i * 1.0e8      )    ! divergence
-      CALL iom_put( "ishear"      , shear_i * 1.0e8     )    ! shear
+      CALL iom_put( "icestr"      , strength * zswi )            ! ice strength
+      CALL iom_put( "idive"       , divu_i              )        ! divergence
+      CALL iom_put( "ishear"      , shear_i             )        ! shear
       CALL iom_put( "snowvol"     , vt_s   * zswi       )        ! snow volume
       
@@ -186,9 +195,9 @@
       CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt 
       CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow/ice) 
-      CALL iom_put( "vfxsub_err" , wfx_err_sub * ztmp   )        ! "excess" of sublimation sent to ocean 
-      
-      CALL iom_put( "afxtot"     , afx_tot * rday       )        ! concentration tendency (total)
-      CALL iom_put( "afxdyn"     , afx_dyn * rday       )        ! concentration tendency (dynamics)
-      CALL iom_put( "afxthd"     , afx_thd * rday       )        ! concentration tendency (thermo)
+      CALL iom_put( "vfxsub_err" , wfx_err_sub * ztmp   )        ! "excess" of sublimation sent to ocean      
+ 
+      CALL iom_put( "afxtot"     , afx_tot              )        ! concentration tendency (total)
+      CALL iom_put( "afxdyn"     , afx_dyn              )        ! concentration tendency (dynamics)
+      CALL iom_put( "afxthd"     , afx_thd              )        ! concentration tendency (thermo)
 
       CALL iom_put ('hfxthd'     , hfx_thd(:,:)         )   !  
@@ -218,8 +227,7 @@
       ! END MV MP 2016
 
-      
-      !--------------------------------
-      ! Output values for each category
-      !--------------------------------
+      !----------------------------------
+      ! Output category-dependent fields
+      !----------------------------------
       IF ( iom_use( "iceconc_cat"  ) )  CALL iom_put( "iceconc_cat"      , a_i   * zswi2   )        ! area for categories
       IF ( iom_use( "icethic_cat"  ) )  CALL iom_put( "icethic_cat"      , ht_i  * zswi2   )        ! thickness for categories
@@ -231,5 +239,5 @@
       IF ( iom_use( "snwtemp_cat"  ) )  CALL iom_put( "snwtemp_cat", ( SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s - rt0 ) * zswi2 )
       ! ice age
-      IF ( iom_use( "iceage_cat"   ) )  CALL iom_put( "iceage_cat" , o_i * zswi2 * z1_365 )
+      IF ( iom_use( "iceage_cat"   ) )  CALL iom_put( "iceage_cat" , o_i * zswi2 ) 
       ! brine volume
       IF ( iom_use( "brinevol_cat" ) )  CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 )
@@ -244,4 +252,117 @@
       ! END MV MP 2016
 
+      !--------------------------------
+      ! Add-ons for SIMIP
+      !--------------------------------
+      zrho1 = ( rau0 - rhoic ) / rau0; zrho2 = rhosn / rau0
+
+      IF  ( iom_use( "icethic"  ) ) CALL iom_put( "icethic"     , htm_i * zamask             )          ! Ice thickness 
+      IF  ( iom_use( "icepres"  ) ) CALL iom_put( "icepres"     , zswi                       )          ! Ice presence (1 or 0) 
+      IF  ( iom_use( "snowthic" ) ) CALL iom_put( "snowthic"    , htm_s * zamask             )          ! Snow thickness       
+      IF  ( iom_use( "icemass"  ) ) CALL iom_put( "icemass"     , rhoic * vt_i(:,:) * zswi   )          ! Ice mass per cell area 
+      IF  ( iom_use( "snomass"  ) ) CALL iom_put( "snomass"     , rhosn * vt_s(:,:) * zswi   )          ! Snow mass per cell area
+      IF  ( iom_use( "icesnt"   ) ) CALL iom_put( "icesnt"      , ( tm_si - rt0 ) * zswi     )          ! Snow-ice interface temperature
+      IF  ( iom_use( "icebot"   ) ) CALL iom_put( "icebot"      , ( t_bo  - rt0 ) * zswi     )          ! Ice bottom temperature
+      IF  ( iom_use( "icesmass" ) ) CALL iom_put( "icesmass"    , SUM( smv_i, DIM = 3 ) * rhoic * 1.0e-3 * zswi )   ! Mass of salt in sea ice per cell area
+      IF  ( iom_use( "icefb"    ) ) THEN
+         zfb(:,:) = ( zrho1 * htm_i(:,:) - zrho2 * htm_s(:,:) ) * zswi(:,:)                             
+         WHERE( zfb < 0._wp ) ;   zfb = 0._wp ;   END WHERE
+                                    CALL iom_put( "icefb"       , zfb                        )          ! Ice freeboard
+      ENDIF
+      IF  ( iom_use( "dmithd"   ) ) CALL iom_put( "dmithd"      , - wfx_bog - wfx_bom - wfx_sum   &     ! Sea-ice mass change from thermodynamics
+              &                     - wfx_sni - wfx_opw - wfx_res )
+      IF  ( iom_use( "dmidyn"   ) ) CALL iom_put( "dmidyn"      ,   diag_dmi_dyn             )          ! Sea-ice mass change from dynamics
+      IF  ( iom_use( "dmiopw"   ) ) CALL iom_put( "dmiopw"      , - wfx_opw                  )          ! Sea-ice mass change through growth in open water
+      IF  ( iom_use( "dmibog"   ) ) CALL iom_put( "dmibog"      , - wfx_bog                  )          ! Sea-ice mass change through basal growth
+      IF  ( iom_use( "dmisni"   ) ) CALL iom_put( "dmisni"      , - wfx_sni                  )          ! Sea-ice mass change through snow-to-ice conversion
+      IF  ( iom_use( "dmisum"   ) ) CALL iom_put( "dmisum"      , - wfx_sum                  )          ! Sea-ice mass change through surface melting
+      IF  ( iom_use( "dmibom"   ) ) CALL iom_put( "dmibom"      , - wfx_bom                  )          ! Sea-ice mass change through bottom melting
+
+      IF  ( iom_use( "dmtsub"   ) ) CALL iom_put( "dmtsub"      , - wfx_sub                  )          ! Sea-ice mass change through evaporation and sublimation
+      IF  ( iom_use( "dmssub"   ) ) CALL iom_put( "dmssub"      , - wfx_snw_sub              )          ! Snow mass change through sublimation
+      IF  ( iom_use( "dmisub"   ) ) CALL iom_put( "dmisub"      , - wfx_ice_sub              )          ! Sea-ice mass change through sublimation
+
+      IF  ( iom_use( "dmsspr"   ) ) CALL iom_put( "dmsspr"      , - wfx_spr                  )          ! Snow mass change through snow fall
+      IF  ( iom_use( "dmsssi"   ) ) CALL iom_put( "dmsssi"      ,   wfx_sni*rhosn/rhoic      )          ! Snow mass change through snow-to-ice conversion
+
+      IF  ( iom_use( "dmsmel"   ) ) CALL iom_put( "dmsmel"      , - wfx_snw_sum              )          ! Snow mass change through melt
+      IF  ( iom_use( "dmsdyn"   ) ) CALL iom_put( "dmsdyn"      ,   diag_dms_dyn             )          ! Snow mass change through dynamics
+
+      IF  ( iom_use( "hfxconbo" ) ) CALL iom_put( "hfxconbo"    ,   diag_fc_bo               )          ! Bottom conduction flux
+      IF  ( iom_use( "hfxconsu" ) ) CALL iom_put( "hfxconsu"    ,   diag_fc_su               )          ! Surface conduction flux
+
+      IF  ( iom_use( "wfxtot"   ) ) CALL iom_put( "wfxtot"      ,   wfx_ice                  )          ! Total freshwater flux from sea ice
+      IF  ( iom_use( "wfxsum"   ) ) CALL iom_put( "wfxsum"      ,   wfx_sum                  )          ! Freshwater flux from sea-ice surface
+
+      IF  ( iom_use( "utau_oi"  ) ) CALL iom_put( "utau_oi"     ,   diag_utau_oi*zswi        )          ! X-component of ocean stress on sea ice
+      IF  ( iom_use( "vtau_oi"  ) ) CALL iom_put( "vtau_oi"     ,   diag_vtau_oi*zswi        )          ! Y-component of ocean stress on sea ice
+
+      IF  ( iom_use( "dssh_dx"  ) ) CALL iom_put( "dssh_dx"     ,   diag_dssh_dx*zswi        )          ! Sea-surface tilt term in force balance (x-component)
+      IF  ( iom_use( "dssh_dy"  ) ) CALL iom_put( "dssh_dy"     ,   diag_dssh_dy*zswi        )          ! Sea-surface tilt term in force balance (y-component)
+
+      IF  ( iom_use( "corstrx"  ) ) CALL iom_put( "corstrx"     ,   diag_corstrx*zswi        )          ! Coriolis force term in force balance (x-component)
+      IF  ( iom_use( "corstry"  ) ) CALL iom_put( "corstry"     ,   diag_corstry*zswi        )          ! Coriolis force term in force balance (y-component)
+
+      IF  ( iom_use( "intstrx"  ) ) CALL iom_put( "intstrx"     ,   diag_intstrx*zswi        )          ! Internal force term in force balance (x-component)
+      IF  ( iom_use( "intstry"  ) ) CALL iom_put( "intstry"     ,   diag_intstry*zswi        )          ! Internal force term in force balance (y-component)
+
+      IF  ( iom_use( "normstr"  ) ) CALL iom_put( "normstr"     ,   diag_sig1   *zswi        )          ! Normal stress
+      IF  ( iom_use( "sheastr"  ) ) CALL iom_put( "sheastr"     ,   diag_sig2   *zswi        )          ! Shear stress
+
+      IF  ( iom_use( "xmtrpice" ) ) CALL iom_put( "xmtrpice"     ,  diag_xmtrp_ice           )          ! X-component of sea-ice mass transport
+      IF  ( iom_use( "ymtrpice" ) ) CALL iom_put( "ymtrpice"     ,  diag_ymtrp_ice           )          ! Y-component of sea-ice mass transport
+
+      IF  ( iom_use( "xmtrpsnw" ) ) CALL iom_put( "xmtrpsnw"     ,  diag_xmtrp_snw           )          ! X-component of snow mass transport
+      IF  ( iom_use( "ymtrpsnw" ) ) CALL iom_put( "ymtrpsnw"     ,  diag_ymtrp_snw           )          ! Y-component of snow mass transport
+
+      IF  ( iom_use( "xatrp"    ) ) CALL iom_put( "xatrp"        ,  diag_xatrp               )          ! X-component of ice area transport
+      IF  ( iom_use( "yatrp"    ) ) CALL iom_put( "yatrp"        ,  diag_yatrp               )          ! Y-component of ice area transport
+
+      !--------------------------------
+      ! Global ice diagnostics (SIMIP)
+      !--------------------------------
+
+      IF ( iom_use( "NH_icearea" ) .OR. iom_use( "NH_icevolu" ) .OR. iom_use( "NH_iceextt" ) )   THEN   ! NH integrated diagnostics
+ 
+         WHERE( fcor > 0._wp ); zswi(:,:) = 1.0e-12
+         ELSEWHERE            ; zswi(:,:) = 0.
+         END WHERE 
+
+         zdiag_area_nh = glob_sum( at_i(:,:) * zswi(:,:) * e12t(:,:) )
+         zdiag_volu_nh = glob_sum( vt_i(:,:) * zswi(:,:) * e12t(:,:) )
+
+         WHERE( fcor > 0._wp .AND. at_i > 0.15 ); zswi(:,:) = 1.0e-12
+         ELSEWHERE                              ; zswi(:,:) = 0.
+         END WHERE 
+
+         zdiag_extt_nh = glob_sum( zswi(:,:) * e12t(:,:) )
+
+         IF ( iom_use( "NH_icearea" ) ) CALL iom_put( "NH_icearea" ,  zdiag_area_nh  )
+         IF ( iom_use( "NH_icevolu" ) ) CALL iom_put( "NH_icevolu" ,  zdiag_volu_nh  )
+         IF ( iom_use( "NH_iceextt" ) ) CALL iom_put( "NH_iceextt" ,  zdiag_extt_nh  )
+
+      ENDIF
+
+      IF ( iom_use( "SH_icearea" ) .OR. iom_use( "SH_icevolu" ) .OR. iom_use( "SH_iceextt" ) )   THEN   ! SH integrated diagnostics
+
+         WHERE( fcor < 0._wp ); zswi(:,:) = 1.0e-12; 
+         ELSEWHERE            ; zswi(:,:) = 0.
+         END WHERE 
+
+         zdiag_area_sh = glob_sum( at_i(:,:) * zswi(:,:) * e12t(:,:) ) 
+         zdiag_volu_sh = glob_sum( vt_i(:,:) * zswi(:,:) * e12t(:,:) )
+
+         WHERE( fcor < 0._wp .AND. at_i > 0.15 ); zswi(:,:) = 1.0e-12
+         ELSEWHERE                              ; zswi(:,:) = 0.
+         END WHERE 
+
+         zdiag_extt_sh = glob_sum( zswi(:,:) * e12t(:,:) )
+
+         IF ( iom_use( "SH_icearea" ) ) CALL iom_put( "SH_icearea", zdiag_area_sh )
+         IF ( iom_use( "SH_icevolu" ) ) CALL iom_put( "SH_icevolu", zdiag_volu_sh )
+         IF ( iom_use( "SH_iceextt" ) ) CALL iom_put( "SH_iceextt", zdiag_extt_sh )
+
+      ENDIF 
+
       !     !  Create an output files (output.lim.abort.nc) if S < 0 or u > 20 m/s
       !     IF( kindic < 0 )   CALL lim_wri_state( 'output.abort' )
@@ -250,4 +371,5 @@
       CALL wrk_dealloc( jpi, jpj, jpl, zswi2 )
       CALL wrk_dealloc( jpi, jpj     , z2d, zswi )
+      CALL wrk_dealloc( jpi, jpj     , zfb, zamask, zamask15 )
 
       IF( nn_timing == 1 )  CALL timing_stop('limwri')
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8239)
@@ -60,4 +60,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_sum_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sub_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_sub_1d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_ice_sub_1d 
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_bog_1d    
@@ -93,4 +95,5 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_su_1d       !: <==> the 2D  t_su
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_si_1d       !: <==> the 2D  t_si
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_i_1d        !: <==> the 2D  a_i
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_i_1d       !: <==> the 2D  ht_s
@@ -115,4 +118,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_old  !: ice thickness layer (m)
 
+   ! Conduction flux diagnostics (SIMIP)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   diag_fc_bo_1d      !: <==> the 2D  diag_fc_bo
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   diag_fc_su_1d      !: <==> the 2D  diag_fc_su
+
    INTEGER , PUBLIC ::   jiindex_1d   ! 1D index of debugging point
 
@@ -129,5 +136,5 @@
       !!---------------------------------------------------------------------!
       INTEGER ::   thd_ice_alloc   ! return value
-      INTEGER ::   ierr(4), ii
+      INTEGER ::   ierr(5), ii
       !!---------------------------------------------------------------------!
       ierr(:) = 0
@@ -150,12 +157,14 @@
          &      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) ,  &
+         &      wfx_snw_sub_1d(jpij), wfx_ice_sub_1d(jpij) ,                                    &
          &      dqns_ice_1d(jpij) , evap_ice_1d (jpij),                                         &
          &      qprec_ice_1d(jpij), qevap_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_res_1d (jpij) , sfx_sub_1d (jpij),  &
-         &      hicol_1d   (jpij) , STAT=ierr(ii) )
+         &      dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) ,                     &     
+         &      dsm_i_si_1d(jpij) , hicol_1d    (jpij)                     , STAT=ierr(ii) )
       !
       ii = ii + 1
-      ALLOCATE( t_su_1d   (jpij) , a_i_1d    (jpij) , ht_i_1d  (jpij) ,                      &
+      ALLOCATE( t_su_1d   (jpij) , t_si_1d  (jpij)  , a_i_1d   (jpij) , ht_i_1d  (jpij) ,    &
          &      ht_s_1d   (jpij) , fc_su     (jpij) , fc_bo_i  (jpij) ,                      &    
          &      dh_s_tot  (jpij) , dh_i_surf (jpij) , dh_i_sub (jpij) ,                      &    
@@ -168,4 +177,7 @@
          &      qh_i_old(jpij,0:nlay_i+1) , h_i_old(jpij,0:nlay_i+1) , STAT=ierr(ii) )
       !
+      ii = ii + 1
+      ALLOCATE( diag_fc_bo_1d(jpij)      , diag_fc_su_1d(jpij)      , STAT=ierr(ii) )
+
       thd_ice_alloc = MAXVAL( ierr(:) )
       IF( thd_ice_alloc /= 0 )   CALL ctl_warn( 'thd_ice_alloc: failed to allocate arrays.' )
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 8239)
@@ -195,4 +195,6 @@
 #if defined key_lim3 || defined key_lim2
       CALL iom_set_axis_attr( "ncatice", (/ (REAL(ji,wp), ji=1,jpl) /) )
+      ! SIMIP diagnostics (4 main arctic straits)
+      CALL iom_set_axis_attr( "nstrait", (/ (REAL(ji,wp), ji=1,4) /) )
 #endif
       CALL iom_set_axis_attr( "icbcla", class_num )
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8237)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8239)
@@ -634,4 +634,6 @@
       wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
       wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp  
+      wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp
+      wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp
 
       ! MV MP 2016
@@ -654,4 +656,8 @@
       diag_heat(:,:) = 0._wp ;   diag_smvi(:,:) = 0._wp
       diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp
+
+      ! SIMIP diagnostics
+      diag_dms_dyn(:,:)  = 0._wp ; diag_dmi_dyn(:,:)  = 0._wp
+      diag_fc_bo(:,:)    = 0._wp ; diag_fc_su(:,:)    = 0._wp
 
       tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
