Index: ../trunk-jpl/m4/analyses.m4
===================================================================
--- ../trunk-jpl/m4/analyses.m4	(revision 18491)
+++ ../trunk-jpl/m4/analyses.m4	(revision 18492)
@@ -10,466 +10,888 @@
 [
 
 dnl with-AdjointBalancethickness{{{
+
 AC_ARG_WITH([AdjointBalancethickness],
+
 	AS_HELP_STRING([--with-AdjointBalancethickness = YES], [compile with AdjointBalancethickness capabilities (default is yes)]),
+
 	[ADJOINTBALANCETHICKNESS=$withval],[ADJOINTBALANCETHICKNESS=yes])
+
 AC_MSG_CHECKING(for AdjointBalancethickness capability compilation)
 
+
 HAVE_ADJOINTBALANCETHICKNESS=no 
+
 if test "x$ADJOINTBALANCETHICKNESS" = "xyes"; then
+
 	HAVE_ADJOINTBALANCETHICKNESS=yes
+
 	AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS_],[1],[with AdjointBalancethicknesscapability])
+
 fi
+
 AM_CONDITIONAL([ADJOINTBALANCETHICKNESS], [test x$HAVE_ADJOINTBALANCETHICKNESS = xyes])
+
 AC_MSG_RESULT($HAVE_ADJOINTBALANCETHICKNESS)
+
 dnl }}}
 dnl with-AdjointBalancethickness2{{{
+
 AC_ARG_WITH([AdjointBalancethickness2],
+
 	AS_HELP_STRING([--with-AdjointBalancethickness2 = YES], [compile with AdjointBalancethickness2 capabilities (default is yes)]),
+
 	[ADJOINTBALANCETHICKNESS2=$withval],[ADJOINTBALANCETHICKNESS2=yes])
+
 AC_MSG_CHECKING(for AdjointBalancethickness2 capability compilation)
 
+
 HAVE_ADJOINTBALANCETHICKNESS2=no 
+
 if test "x$ADJOINTBALANCETHICKNESS2" = "xyes"; then
+
 	HAVE_ADJOINTBALANCETHICKNESS2=yes
+
 	AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS2_],[1],[with AdjointBalancethickness2capability])
+
 fi
+
 AM_CONDITIONAL([ADJOINTBALANCETHICKNESS2], [test x$HAVE_ADJOINTBALANCETHICKNESS2 = xyes])
+
 AC_MSG_RESULT($HAVE_ADJOINTBALANCETHICKNESS2)
+
 dnl }}}
 dnl with-AdjointHoriz{{{
+
 AC_ARG_WITH([AdjointHoriz],
+
 	AS_HELP_STRING([--with-AdjointHoriz = YES], [compile with AdjointHoriz capabilities (default is yes)]),
+
 	[ADJOINTHORIZ=$withval],[ADJOINTHORIZ=yes])
+
 AC_MSG_CHECKING(for AdjointHoriz capability compilation)
 
+
 HAVE_ADJOINTHORIZ=no 
+
 if test "x$ADJOINTHORIZ" = "xyes"; then
+
 	HAVE_ADJOINTHORIZ=yes
+
 	AC_DEFINE([_HAVE_ADJOINTHORIZ_],[1],[with AdjointHorizcapability])
+
 fi
+
 AM_CONDITIONAL([ADJOINTHORIZ], [test x$HAVE_ADJOINTHORIZ = xyes])
+
 AC_MSG_RESULT($HAVE_ADJOINTHORIZ)
+
 dnl }}}
 dnl with-Balancethickness{{{
+
 AC_ARG_WITH([Balancethickness],
+
 	AS_HELP_STRING([--with-Balancethickness = YES], [compile with Balancethickness capabilities (default is yes)]),
+
 	[BALANCETHICKNESS=$withval],[BALANCETHICKNESS=yes])
+
 AC_MSG_CHECKING(for Balancethickness capability compilation)
 
+
 HAVE_BALANCETHICKNESS=no 
+
 if test "x$BALANCETHICKNESS" = "xyes"; then
+
 	HAVE_BALANCETHICKNESS=yes
+
 	AC_DEFINE([_HAVE_BALANCETHICKNESS_],[1],[with Balancethicknesscapability])
+
 fi
+
 AM_CONDITIONAL([BALANCETHICKNESS], [test x$HAVE_BALANCETHICKNESS = xyes])
+
 AC_MSG_RESULT($HAVE_BALANCETHICKNESS)
+
 dnl }}}
 dnl with-Balancethickness2{{{
+
 AC_ARG_WITH([Balancethickness2],
+
 	AS_HELP_STRING([--with-Balancethickness2 = YES], [compile with Balancethickness2 capabilities (default is yes)]),
+
 	[BALANCETHICKNESS2=$withval],[BALANCETHICKNESS2=yes])
+
 AC_MSG_CHECKING(for Balancethickness2 capability compilation)
 
+
 HAVE_BALANCETHICKNESS2=no 
+
 if test "x$BALANCETHICKNESS2" = "xyes"; then
+
 	HAVE_BALANCETHICKNESS2=yes
+
 	AC_DEFINE([_HAVE_BALANCETHICKNESS2_],[1],[with Balancethickness2capability])
+
 fi
+
 AM_CONDITIONAL([BALANCETHICKNESS2], [test x$HAVE_BALANCETHICKNESS2 = xyes])
+
 AC_MSG_RESULT($HAVE_BALANCETHICKNESS2)
+
 dnl }}}
 dnl with-BalancethicknessSoft{{{
+
 AC_ARG_WITH([BalancethicknessSoft],
+
 	AS_HELP_STRING([--with-BalancethicknessSoft = YES], [compile with BalancethicknessSoft capabilities (default is yes)]),
+
 	[BALANCETHICKNESSSOFT=$withval],[BALANCETHICKNESSSOFT=yes])
+
 AC_MSG_CHECKING(for BalancethicknessSoft capability compilation)
 
+
 HAVE_BALANCETHICKNESSSOFT=no 
+
 if test "x$BALANCETHICKNESSSOFT" = "xyes"; then
+
 	HAVE_BALANCETHICKNESSSOFT=yes
+
 	AC_DEFINE([_HAVE_BALANCETHICKNESSSOFT_],[1],[with BalancethicknessSoftcapability])
+
 fi
+
 AM_CONDITIONAL([BALANCETHICKNESSSOFT], [test x$HAVE_BALANCETHICKNESSSOFT = xyes])
+
 AC_MSG_RESULT($HAVE_BALANCETHICKNESSSOFT)
+
 dnl }}}
 dnl with-Balancevelocity{{{
+
 AC_ARG_WITH([Balancevelocity],
+
 	AS_HELP_STRING([--with-Balancevelocity = YES], [compile with Balancevelocity capabilities (default is yes)]),
+
 	[BALANCEVELOCITY=$withval],[BALANCEVELOCITY=yes])
+
 AC_MSG_CHECKING(for Balancevelocity capability compilation)
 
+
 HAVE_BALANCEVELOCITY=no 
+
 if test "x$BALANCEVELOCITY" = "xyes"; then
+
 	HAVE_BALANCEVELOCITY=yes
+
 	AC_DEFINE([_HAVE_BALANCEVELOCITY_],[1],[with Balancevelocitycapability])
+
 fi
+
 AM_CONDITIONAL([BALANCEVELOCITY], [test x$HAVE_BALANCEVELOCITY = xyes])
+
 AC_MSG_RESULT($HAVE_BALANCEVELOCITY)
+
 dnl }}}
 dnl with-L2ProjectionEPL{{{
+
 AC_ARG_WITH([L2ProjectionEPL],
+
 	AS_HELP_STRING([--with-L2ProjectionEPL = YES], [compile with L2ProjectionEPL capabilities (default is yes)]),
+
 	[L2PROJECTIONEPL=$withval],[L2PROJECTIONEPL=yes])
+
 AC_MSG_CHECKING(for L2ProjectionEPL capability compilation)
 
+
 HAVE_L2PROJECTIONEPL=no 
+
 if test "x$L2PROJECTIONEPL" = "xyes"; then
+
 	HAVE_L2PROJECTIONEPL=yes
+
 	AC_DEFINE([_HAVE_L2PROJECTIONEPL_],[1],[with L2ProjectionEPLcapability])
+
 fi
+
 AM_CONDITIONAL([L2PROJECTIONEPL], [test x$HAVE_L2PROJECTIONEPL = xyes])
+
 AC_MSG_RESULT($HAVE_L2PROJECTIONEPL)
+
 dnl }}}
 dnl with-L2ProjectionBase{{{
+
 AC_ARG_WITH([L2ProjectionBase],
+
 	AS_HELP_STRING([--with-L2ProjectionBase = YES], [compile with L2ProjectionBase capabilities (default is yes)]),
+
 	[L2PROJECTIONBASE=$withval],[L2PROJECTIONBASE=yes])
+
 AC_MSG_CHECKING(for L2ProjectionBase capability compilation)
 
+
 HAVE_L2PROJECTIONBASE=no 
+
 if test "x$L2PROJECTIONBASE" = "xyes"; then
+
 	HAVE_L2PROJECTIONBASE=yes
+
 	AC_DEFINE([_HAVE_L2PROJECTIONBASE_],[1],[with L2ProjectionBasecapability])
+
 fi
+
 AM_CONDITIONAL([L2PROJECTIONBASE], [test x$HAVE_L2PROJECTIONBASE = xyes])
+
 AC_MSG_RESULT($HAVE_L2PROJECTIONBASE)
+
 dnl }}}
 dnl with-DamageEvolution{{{
+
 AC_ARG_WITH([DamageEvolution],
+
 	AS_HELP_STRING([--with-DamageEvolution = YES], [compile with DamageEvolution capabilities (default is yes)]),
+
 	[DAMAGEEVOLUTION=$withval],[DAMAGEEVOLUTION=yes])
+
 AC_MSG_CHECKING(for DamageEvolution capability compilation)
 
+
 HAVE_DAMAGEEVOLUTION=no 
+
 if test "x$DAMAGEEVOLUTION" = "xyes"; then
+
 	HAVE_DAMAGEEVOLUTION=yes
+
 	AC_DEFINE([_HAVE_DAMAGEEVOLUTION_],[1],[with DamageEvolutioncapability])
+
 fi
+
 AM_CONDITIONAL([DAMAGEEVOLUTION], [test x$HAVE_DAMAGEEVOLUTION = xyes])
+
 AC_MSG_RESULT($HAVE_DAMAGEEVOLUTION)
+
 dnl }}}
 dnl with-Stressbalance{{{
+
 AC_ARG_WITH([Stressbalance],
+
 	AS_HELP_STRING([--with-Stressbalance = YES], [compile with Stressbalance capabilities (default is yes)]),
+
 	[STRESSBALANCE=$withval],[STRESSBALANCE=yes])
+
 AC_MSG_CHECKING(for Stressbalance capability compilation)
 
+
 HAVE_STRESSBALANCE=no 
+
 if test "x$STRESSBALANCE" = "xyes"; then
+
 	HAVE_STRESSBALANCE=yes
+
 	AC_DEFINE([_HAVE_STRESSBALANCE_],[1],[with Stressbalancecapability])
+
 fi
+
 AM_CONDITIONAL([STRESSBALANCE], [test x$HAVE_STRESSBALANCE = xyes])
+
 AC_MSG_RESULT($HAVE_STRESSBALANCE)
+
 dnl }}}
 dnl with-StressbalanceSIA{{{
+
 AC_ARG_WITH([StressbalanceSIA],
+
 	AS_HELP_STRING([--with-StressbalanceSIA = YES], [compile with StressbalanceSIA capabilities (default is yes)]),
+
 	[STRESSBALANCESIA=$withval],[STRESSBALANCESIA=yes])
+
 AC_MSG_CHECKING(for StressbalanceSIA capability compilation)
 
+
 HAVE_STRESSBALANCESIA=no 
+
 if test "x$STRESSBALANCESIA" = "xyes"; then
+
 	HAVE_STRESSBALANCESIA=yes
+
 	AC_DEFINE([_HAVE_STRESSBALANCESIA_],[1],[with StressbalanceSIAcapability])
+
 fi
+
 AM_CONDITIONAL([STRESSBALANCESIA], [test x$HAVE_STRESSBALANCESIA = xyes])
+
 AC_MSG_RESULT($HAVE_STRESSBALANCESIA)
+
 dnl }}}
 dnl with-StressbalanceVertical{{{
+
 AC_ARG_WITH([StressbalanceVertical],
+
 	AS_HELP_STRING([--with-StressbalanceVertical = YES], [compile with StressbalanceVertical capabilities (default is yes)]),
+
 	[STRESSBALANCEVERTICAL=$withval],[STRESSBALANCEVERTICAL=yes])
+
 AC_MSG_CHECKING(for StressbalanceVertical capability compilation)
 
+
 HAVE_STRESSBALANCEVERTICAL=no 
+
 if test "x$STRESSBALANCEVERTICAL" = "xyes"; then
+
 	HAVE_STRESSBALANCEVERTICAL=yes
+
 	AC_DEFINE([_HAVE_STRESSBALANCEVERTICAL_],[1],[with StressbalanceVerticalcapability])
+
 fi
+
 AM_CONDITIONAL([STRESSBALANCEVERTICAL], [test x$HAVE_STRESSBALANCEVERTICAL = xyes])
+
 AC_MSG_RESULT($HAVE_STRESSBALANCEVERTICAL)
+
 dnl }}}
 dnl with-Enthalpy{{{
+
 AC_ARG_WITH([Enthalpy],
+
 	AS_HELP_STRING([--with-Enthalpy = YES], [compile with Enthalpy capabilities (default is yes)]),
+
 	[ENTHALPY=$withval],[ENTHALPY=yes])
+
 AC_MSG_CHECKING(for Enthalpy capability compilation)
 
+
 HAVE_ENTHALPY=no 
+
 if test "x$ENTHALPY" = "xyes"; then
+
 	HAVE_ENTHALPY=yes
+
 	AC_DEFINE([_HAVE_ENTHALPY_],[1],[with Enthalpycapability])
+
 fi
+
 AM_CONDITIONAL([ENTHALPY], [test x$HAVE_ENTHALPY = xyes])
+
 AC_MSG_RESULT($HAVE_ENTHALPY)
+
 dnl }}}
 dnl with-HydrologyShreve{{{
+
 AC_ARG_WITH([HydrologyShreve],
+
 	AS_HELP_STRING([--with-HydrologyShreve = YES], [compile with HydrologyShreve capabilities (default is yes)]),
+
 	[HYDROLOGYSHREVE=$withval],[HYDROLOGYSHREVE=yes])
+
 AC_MSG_CHECKING(for HydrologyShreve capability compilation)
 
+
 HAVE_HYDROLOGYSHREVE=no 
+
 if test "x$HYDROLOGYSHREVE" = "xyes"; then
+
 	HAVE_HYDROLOGYSHREVE=yes
+
 	AC_DEFINE([_HAVE_HYDROLOGYSHREVE_],[1],[with HydrologyShrevecapability])
+
 fi
+
 AM_CONDITIONAL([HYDROLOGYSHREVE], [test x$HAVE_HYDROLOGYSHREVE = xyes])
+
 AC_MSG_RESULT($HAVE_HYDROLOGYSHREVE)
+
 dnl }}}
 dnl with-HydrologyDCInefficient{{{
+
 AC_ARG_WITH([HydrologyDCInefficient],
+
 	AS_HELP_STRING([--with-HydrologyDCInefficient = YES], [compile with HydrologyDCInefficient capabilities (default is yes)]),
+
 	[HYDROLOGYDCINEFFICIENT=$withval],[HYDROLOGYDCINEFFICIENT=yes])
+
 AC_MSG_CHECKING(for HydrologyDCInefficient capability compilation)
 
+
 HAVE_HYDROLOGYDCINEFFICIENT=no 
+
 if test "x$HYDROLOGYDCINEFFICIENT" = "xyes"; then
+
 	HAVE_HYDROLOGYDCINEFFICIENT=yes
+
 	AC_DEFINE([_HAVE_HYDROLOGYDCINEFFICIENT_],[1],[with HydrologyDCInefficientcapability])
+
 fi
+
 AM_CONDITIONAL([HYDROLOGYDCINEFFICIENT], [test x$HAVE_HYDROLOGYDCINEFFICIENT = xyes])
+
 AC_MSG_RESULT($HAVE_HYDROLOGYDCINEFFICIENT)
+
 dnl }}}
 dnl with-HydrologyDCEfficient{{{
+
 AC_ARG_WITH([HydrologyDCEfficient],
+
 	AS_HELP_STRING([--with-HydrologyDCEfficient = YES], [compile with HydrologyDCEfficient capabilities (default is yes)]),
+
 	[HYDROLOGYDCEFFICIENT=$withval],[HYDROLOGYDCEFFICIENT=yes])
+
 AC_MSG_CHECKING(for HydrologyDCEfficient capability compilation)
 
+
 HAVE_HYDROLOGYDCEFFICIENT=no 
+
 if test "x$HYDROLOGYDCEFFICIENT" = "xyes"; then
+
 	HAVE_HYDROLOGYDCEFFICIENT=yes
+
 	AC_DEFINE([_HAVE_HYDROLOGYDCEFFICIENT_],[1],[with HydrologyDCEfficientcapability])
+
 fi
+
 AM_CONDITIONAL([HYDROLOGYDCEFFICIENT], [test x$HAVE_HYDROLOGYDCEFFICIENT = xyes])
+
 AC_MSG_RESULT($HAVE_HYDROLOGYDCEFFICIENT)
+
 dnl }}}
 dnl with-Melting{{{
+
 AC_ARG_WITH([Melting],
+
 	AS_HELP_STRING([--with-Melting = YES], [compile with Melting capabilities (default is yes)]),
+
 	[MELTING=$withval],[MELTING=yes])
+
 AC_MSG_CHECKING(for Melting capability compilation)
 
+
 HAVE_MELTING=no 
+
 if test "x$MELTING" = "xyes"; then
+
 	HAVE_MELTING=yes
+
 	AC_DEFINE([_HAVE_MELTING_],[1],[with Meltingcapability])
+
 fi
+
 AM_CONDITIONAL([MELTING], [test x$HAVE_MELTING = xyes])
+
 AC_MSG_RESULT($HAVE_MELTING)
+
 dnl }}}
 dnl with-Masstransport{{{
+
 AC_ARG_WITH([Masstransport],
+
 	AS_HELP_STRING([--with-Masstransport = YES], [compile with Masstransport capabilities (default is yes)]),
+
 	[MASSTRANSPORT=$withval],[MASSTRANSPORT=yes])
+
 AC_MSG_CHECKING(for Masstransport capability compilation)
 
+
 HAVE_MASSTRANSPORT=no 
+
 if test "x$MASSTRANSPORT" = "xyes"; then
+
 	HAVE_MASSTRANSPORT=yes
+
 	AC_DEFINE([_HAVE_MASSTRANSPORT_],[1],[with Masstransportcapability])
+
 fi
+
 AM_CONDITIONAL([MASSTRANSPORT], [test x$HAVE_MASSTRANSPORT = xyes])
+
 AC_MSG_RESULT($HAVE_MASSTRANSPORT)
+
 dnl }}}
 dnl with-FreeSurfaceBase{{{
+
 AC_ARG_WITH([FreeSurfaceBase],
+
 	AS_HELP_STRING([--with-FreeSurfaceBase = YES], [compile with FreeSurfaceBase capabilities (default is yes)]),
+
 	[FREESURFACEBASE=$withval],[FREESURFACEBASE=yes])
+
 AC_MSG_CHECKING(for FreeSurfaceBase capability compilation)
 
+
 HAVE_FREESURFACEBASE=no 
+
 if test "x$FREESURFACEBASE" = "xyes"; then
+
 	HAVE_FREESURFACEBASE=yes
+
 	AC_DEFINE([_HAVE_FREESURFACEBASE_],[1],[with FreeSurfaceBasecapability])
+
 fi
+
 AM_CONDITIONAL([FREESURFACEBASE], [test x$HAVE_FREESURFACEBASE = xyes])
+
 AC_MSG_RESULT($HAVE_FREESURFACEBASE)
+
 dnl }}}
 dnl with-FreeSurfaceTop{{{
+
 AC_ARG_WITH([FreeSurfaceTop],
+
 	AS_HELP_STRING([--with-FreeSurfaceTop = YES], [compile with FreeSurfaceTop capabilities (default is yes)]),
+
 	[FREESURFACETOP=$withval],[FREESURFACETOP=yes])
+
 AC_MSG_CHECKING(for FreeSurfaceTop capability compilation)
 
+
 HAVE_FREESURFACETOP=no 
+
 if test "x$FREESURFACETOP" = "xyes"; then
+
 	HAVE_FREESURFACETOP=yes
+
 	AC_DEFINE([_HAVE_FREESURFACETOP_],[1],[with FreeSurfaceTopcapability])
+
 fi
+
 AM_CONDITIONAL([FREESURFACETOP], [test x$HAVE_FREESURFACETOP = xyes])
+
 AC_MSG_RESULT($HAVE_FREESURFACETOP)
+
 dnl }}}
 dnl with-ExtrudeFromBase{{{
+
 AC_ARG_WITH([ExtrudeFromBase],
+
 	AS_HELP_STRING([--with-ExtrudeFromBase = YES], [compile with ExtrudeFromBase capabilities (default is yes)]),
+
 	[EXTRUDEFROMBASE=$withval],[EXTRUDEFROMBASE=yes])
+
 AC_MSG_CHECKING(for ExtrudeFromBase capability compilation)
 
+
 HAVE_EXTRUDEFROMBASE=no 
+
 if test "x$EXTRUDEFROMBASE" = "xyes"; then
+
 	HAVE_EXTRUDEFROMBASE=yes
+
 	AC_DEFINE([_HAVE_EXTRUDEFROMBASE_],[1],[with ExtrudeFromBasecapability])
+
 fi
+
 AM_CONDITIONAL([EXTRUDEFROMBASE], [test x$HAVE_EXTRUDEFROMBASE = xyes])
+
 AC_MSG_RESULT($HAVE_EXTRUDEFROMBASE)
+
 dnl }}}
 dnl with-ExtrudeFromTop{{{
+
 AC_ARG_WITH([ExtrudeFromTop],
+
 	AS_HELP_STRING([--with-ExtrudeFromTop = YES], [compile with ExtrudeFromTop capabilities (default is yes)]),
+
 	[EXTRUDEFROMTOP=$withval],[EXTRUDEFROMTOP=yes])
+
 AC_MSG_CHECKING(for ExtrudeFromTop capability compilation)
 
+
 HAVE_EXTRUDEFROMTOP=no 
+
 if test "x$EXTRUDEFROMTOP" = "xyes"; then
+
 	HAVE_EXTRUDEFROMTOP=yes
+
 	AC_DEFINE([_HAVE_EXTRUDEFROMTOP_],[1],[with ExtrudeFromTopcapability])
+
 fi
+
 AM_CONDITIONAL([EXTRUDEFROMTOP], [test x$HAVE_EXTRUDEFROMTOP = xyes])
+
 AC_MSG_RESULT($HAVE_EXTRUDEFROMTOP)
+
 dnl }}}
 dnl with-DepthAverage{{{
+
 AC_ARG_WITH([DepthAverage],
+
 	AS_HELP_STRING([--with-DepthAverage = YES], [compile with DepthAverage capabilities (default is yes)]),
+
 	[DEPTHAVERAGE=$withval],[DEPTHAVERAGE=yes])
+
 AC_MSG_CHECKING(for DepthAverage capability compilation)
 
+
 HAVE_DEPTHAVERAGE=no 
+
 if test "x$DEPTHAVERAGE" = "xyes"; then
+
 	HAVE_DEPTHAVERAGE=yes
+
 	AC_DEFINE([_HAVE_DEPTHAVERAGE_],[1],[with DepthAveragecapability])
+
 fi
+
 AM_CONDITIONAL([DEPTHAVERAGE], [test x$HAVE_DEPTHAVERAGE = xyes])
+
 AC_MSG_RESULT($HAVE_DEPTHAVERAGE)
+
 dnl }}}
 dnl with-SmoothedSurfaceSlopeX{{{
+
 AC_ARG_WITH([SmoothedSurfaceSlopeX],
+
 	AS_HELP_STRING([--with-SmoothedSurfaceSlopeX = YES], [compile with SmoothedSurfaceSlopeX capabilities (default is yes)]),
+
 	[SMOOTHEDSURFACESLOPEX=$withval],[SMOOTHEDSURFACESLOPEX=yes])
+
 AC_MSG_CHECKING(for SmoothedSurfaceSlopeX capability compilation)
 
+
 HAVE_SMOOTHEDSURFACESLOPEX=no 
+
 if test "x$SMOOTHEDSURFACESLOPEX" = "xyes"; then
+
 	HAVE_SMOOTHEDSURFACESLOPEX=yes
+
 	AC_DEFINE([_HAVE_SMOOTHEDSURFACESLOPEX_],[1],[with SmoothedSurfaceSlopeXcapability])
+
 fi
+
 AM_CONDITIONAL([SMOOTHEDSURFACESLOPEX], [test x$HAVE_SMOOTHEDSURFACESLOPEX = xyes])
+
 AC_MSG_RESULT($HAVE_SMOOTHEDSURFACESLOPEX)
+
 dnl }}}
 dnl with-SmoothedSurfaceSlopeY{{{
+
 AC_ARG_WITH([SmoothedSurfaceSlopeY],
+
 	AS_HELP_STRING([--with-SmoothedSurfaceSlopeY = YES], [compile with SmoothedSurfaceSlopeY capabilities (default is yes)]),
+
 	[SMOOTHEDSURFACESLOPEY=$withval],[SMOOTHEDSURFACESLOPEY=yes])
+
 AC_MSG_CHECKING(for SmoothedSurfaceSlopeY capability compilation)
 
+
 HAVE_SMOOTHEDSURFACESLOPEY=no 
+
 if test "x$SMOOTHEDSURFACESLOPEY" = "xyes"; then
+
 	HAVE_SMOOTHEDSURFACESLOPEY=yes
+
 	AC_DEFINE([_HAVE_SMOOTHEDSURFACESLOPEY_],[1],[with SmoothedSurfaceSlopeYcapability])
+
 fi
+
 AM_CONDITIONAL([SMOOTHEDSURFACESLOPEY], [test x$HAVE_SMOOTHEDSURFACESLOPEY = xyes])
+
 AC_MSG_RESULT($HAVE_SMOOTHEDSURFACESLOPEY)
-dnl }}}
-dnl with-UzawaPressure{{{
-AC_ARG_WITH([UzawaPressure],
-	AS_HELP_STRING([--with-UzawaPressure = YES], [compile with UzawaPressure capabilities (default is yes)]),
-	[UZAWAPRESSURE=$withval],[UZAWAPRESSURE=yes])
-AC_MSG_CHECKING(for UzawaPressure capability compilation)
 
-HAVE_UZAWAPRESSURE=no 
-if test "x$UZAWAPRESSURE" = "xyes"; then
-	HAVE_UZAWAPRESSURE=yes
-	AC_DEFINE([_HAVE_UZAWAPRESSURE_],[1],[with UzawaPressurecapability])
-fi
-AM_CONDITIONAL([UZAWAPRESSURE], [test x$HAVE_UZAWAPRESSURE = xyes])
-AC_MSG_RESULT($HAVE_UZAWAPRESSURE)
 dnl }}}
 dnl with-Thermal{{{
+
 AC_ARG_WITH([Thermal],
+
 	AS_HELP_STRING([--with-Thermal = YES], [compile with Thermal capabilities (default is yes)]),
+
 	[THERMAL=$withval],[THERMAL=yes])
+
 AC_MSG_CHECKING(for Thermal capability compilation)
 
+
 HAVE_THERMAL=no 
+
 if test "x$THERMAL" = "xyes"; then
+
 	HAVE_THERMAL=yes
+
 	AC_DEFINE([_HAVE_THERMAL_],[1],[with Thermalcapability])
+
 fi
+
 AM_CONDITIONAL([THERMAL], [test x$HAVE_THERMAL = xyes])
+
 AC_MSG_RESULT($HAVE_THERMAL)
+
 dnl }}}
+dnl with-UzawaPressure{{{
+
+AC_ARG_WITH([UzawaPressure],
+
+	AS_HELP_STRING([--with-UzawaPressure = YES], [compile with UzawaPressure capabilities (default is yes)]),
+
+	[UZAWAPRESSURE=$withval],[UZAWAPRESSURE=yes])
+
+AC_MSG_CHECKING(for UzawaPressure capability compilation)
+
+
+HAVE_UZAWAPRESSURE=no 
+
+if test "x$UZAWAPRESSURE" = "xyes"; then
+
+	HAVE_UZAWAPRESSURE=yes
+
+	AC_DEFINE([_HAVE_UZAWAPRESSURE_],[1],[with UzawaPressurecapability])
+
+fi
+
+AM_CONDITIONAL([UZAWAPRESSURE], [test x$HAVE_UZAWAPRESSURE = xyes])
+
+AC_MSG_RESULT($HAVE_UZAWAPRESSURE)
+
+dnl }}}
 dnl with-Gia{{{
+
 AC_ARG_WITH([Gia],
+
 	AS_HELP_STRING([--with-Gia = YES], [compile with Gia capabilities (default is yes)]),
+
 	[GIA=$withval],[GIA=yes])
+
 AC_MSG_CHECKING(for Gia capability compilation)
 
+
 HAVE_GIA=no 
+
 if test "x$GIA" = "xyes"; then
+
 	HAVE_GIA=yes
+
 	AC_DEFINE([_HAVE_GIA_],[1],[with Giacapability])
+
 fi
+
 AM_CONDITIONAL([GIA], [test x$HAVE_GIA = xyes])
+
 AC_MSG_RESULT($HAVE_GIA)
+
 dnl }}}
+dnl with-Seaice{{{
+
+AC_ARG_WITH([Seaice],
+
+	AS_HELP_STRING([--with-Seaice = YES], [compile with Seaice capabilities (default is yes)]),
+
+	[SEAICE=$withval],[SEAICE=yes])
+
+AC_MSG_CHECKING(for Seaice capability compilation)
+
+
+HAVE_SEAICE=no 
+
+if test "x$SEAICE" = "xyes"; then
+
+	HAVE_SEAICE=yes
+
+	AC_DEFINE([_HAVE_SEAICE_],[1],[with Seaicecapability])
+
+fi
+
+AM_CONDITIONAL([SEAICE], [test x$HAVE_SEAICE = xyes])
+
+AC_MSG_RESULT($HAVE_SEAICE)
+
+dnl }}}
 dnl with-Meshdeformation{{{
+
 AC_ARG_WITH([Meshdeformation],
+
 	AS_HELP_STRING([--with-Meshdeformation = YES], [compile with Meshdeformation capabilities (default is yes)]),
+
 	[MESHDEFORMATION=$withval],[MESHDEFORMATION=yes])
+
 AC_MSG_CHECKING(for Meshdeformation capability compilation)
 
+
 HAVE_MESHDEFORMATION=no 
+
 if test "x$MESHDEFORMATION" = "xyes"; then
+
 	HAVE_MESHDEFORMATION=yes
+
 	AC_DEFINE([_HAVE_MESHDEFORMATION_],[1],[with Meshdeformationcapability])
+
 fi
+
 AM_CONDITIONAL([MESHDEFORMATION], [test x$HAVE_MESHDEFORMATION = xyes])
+
 AC_MSG_RESULT($HAVE_MESHDEFORMATION)
+
 dnl }}}
 dnl with-Levelset{{{
+
 AC_ARG_WITH([Levelset],
+
 	AS_HELP_STRING([--with-Levelset = YES], [compile with Levelset capabilities (default is yes)]),
+
 	[LEVELSET=$withval],[LEVELSET=yes])
+
 AC_MSG_CHECKING(for Levelset capability compilation)
 
+
 HAVE_LEVELSET=no 
+
 if test "x$LEVELSET" = "xyes"; then
+
 	HAVE_LEVELSET=yes
+
 	AC_DEFINE([_HAVE_LEVELSET_],[1],[with Levelsetcapability])
+
 fi
+
 AM_CONDITIONAL([LEVELSET], [test x$HAVE_LEVELSET = xyes])
+
 AC_MSG_RESULT($HAVE_LEVELSET)
+
 dnl }}}
 dnl with-Extrapolation{{{
+
 AC_ARG_WITH([Extrapolation],
+
 	AS_HELP_STRING([--with-Extrapolation = YES], [compile with Extrapolation capabilities (default is yes)]),
+
 	[EXTRAPOLATION=$withval],[EXTRAPOLATION=yes])
+
 AC_MSG_CHECKING(for Extrapolation capability compilation)
 
+
 HAVE_EXTRAPOLATION=no 
+
 if test "x$EXTRAPOLATION" = "xyes"; then
+
 	HAVE_EXTRAPOLATION=yes
+
 	AC_DEFINE([_HAVE_EXTRAPOLATION_],[1],[with Extrapolationcapability])
+
 fi
+
 AM_CONDITIONAL([EXTRAPOLATION], [test x$HAVE_EXTRAPOLATION = xyes])
+
 AC_MSG_RESULT($HAVE_EXTRAPOLATION)
+
 dnl }}}
 dnl with-LsfReinitialization{{{
+
 AC_ARG_WITH([LsfReinitialization],
+
 	AS_HELP_STRING([--with-LsfReinitialization = YES], [compile with LsfReinitialization capabilities (default is yes)]),
+
 	[LSFREINITIALIZATION=$withval],[LSFREINITIALIZATION=yes])
+
 AC_MSG_CHECKING(for LsfReinitialization capability compilation)
 
+
 HAVE_LSFREINITIALIZATION=no 
+
 if test "x$LSFREINITIALIZATION" = "xyes"; then
+
 	HAVE_LSFREINITIALIZATION=yes
+
 	AC_DEFINE([_HAVE_LSFREINITIALIZATION_],[1],[with LsfReinitializationcapability])
+
 fi
+
 AM_CONDITIONAL([LSFREINITIALIZATION], [test x$HAVE_LSFREINITIALIZATION = xyes])
+
 AC_MSG_RESULT($HAVE_LSFREINITIALIZATION)
+
 dnl }}}
 
 ])
Index: ../trunk-jpl/src/c/analyses/SeaiceAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/SeaiceAnalysis.h	(revision 0)
+++ ../trunk-jpl/src/c/analyses/SeaiceAnalysis.h	(revision 18492)
@@ -0,0 +1,38 @@
+/*! \file SeaiceAnalysis.h 
+ *  \brief: header file for generic external result object
+ */
+
+#ifndef _SeaiceAnalysis_
+#define _SeaiceAnalysis_
+
+/*Headers*/
+#include "./Analysis.h"
+
+class SeaiceAnalysis: public Analysis{
+
+	public:
+		/*Model processing*/
+		int  DofsPerNode(int** doflist,int domaintype,int approximation);
+		void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
+		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
+		void CreateNodes(Nodes* nodes,IoModel* iomodel);
+		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
+		void CreateLoads(Loads* loads, IoModel* iomodel);
+
+		/*Finite element Analysis*/
+		void           Core(FemModel* femmodel);
+		ElementVector* CreateDVector(Element* element);
+		ElementMatrix* CreateJacobianMatrix(Element* element);
+		ElementMatrix* CreateKMatrix(Element* element);
+		ElementVector* CreatePVector(Element* element);
+		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
+		void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
+		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
+		void UpdateConstraints(FemModel* femmodel);
+
+		/*Sea ice specifics*/
+		void CreateCTensor(IssmDouble* C,Element* element,Gauss* gauss);
+		void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void GetM(IssmDouble* M,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+};
+#endif
Index: ../trunk-jpl/src/c/analyses/EnumToAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/EnumToAnalysis.cpp	(revision 18491)
+++ ../trunk-jpl/src/c/analyses/EnumToAnalysis.cpp	(revision 18492)
@@ -100,6 +100,9 @@
 		#ifdef _HAVE_GIA_
 		case GiaAnalysisEnum : return new GiaAnalysis();
 		#endif
+		#ifdef _HAVE_SEAICE_
+		case SeaiceAnalysisEnum : return new SeaiceAnalysis();
+		#endif
 		#ifdef _HAVE_MESHDEFORMATION_
 		case MeshdeformationAnalysisEnum : return new MeshdeformationAnalysis();
 		#endif
Index: ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 0)
+++ ../trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp	(revision 18492)
@@ -0,0 +1,343 @@
+#include "./SeaiceAnalysis.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+
+/*Model processing*/
+int  SeaiceAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
+	return 2;
+}/*}}}*/
+void SeaiceAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
+}/*}}}*/
+void SeaiceAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
+
+	/*Make sure we are in 2D horizontal*/
+	if(iomodel->domaintype!=Domain2DhorizontalEnum) _error_("Only 2D model are supported");
+
+	/*Update the elements so that they are aware of this analysis*/
+	int counter=0;
+	for(int i=0;i<iomodel->numberofelements;i++){
+		if(iomodel->my_elements[i]){
+			Element* element=(Element*)elements->GetObjectByOffset(counter);
+			element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
+			counter++;
+		}
+	}
+
+	/*Push to the elements all the inputs that we need*/
+	iomodel->FetchDataToInput(elements,SeaiceThicknessEnum);
+	iomodel->FetchDataToInput(elements,SeaiceConcentrationEnum);
+	iomodel->FetchDataToInput(elements,BasalforcingsOceanSshEnum);
+	iomodel->FetchDataToInput(elements,BasalforcingsOceanVxEnum);
+	iomodel->FetchDataToInput(elements,BasalforcingsOceanVyEnum);
+	iomodel->FetchDataToInput(elements,SurfaceforcingsWindVxEnum);
+	iomodel->FetchDataToInput(elements,SurfaceforcingsWindVyEnum);
+	iomodel->FetchDataToInput(elements,MaterialsDamageEnum);
+	iomodel->FetchDataToInput(elements,VxStarEnum);
+	iomodel->FetchDataToInput(elements,VyStarEnum);
+	iomodel->FetchDataToInput(elements,VxEnum);
+	iomodel->FetchDataToInput(elements,VyEnum);
+}/*}}}*/
+void SeaiceAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+
+	::CreateNodes(nodes,iomodel,SeaiceAnalysisEnum,P1Enum);
+
+}/*}}}*/
+void SeaiceAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+
+	const int finiteelement = P1Enum;
+	IoModelToConstraintsx(constraints,iomodel,SeaiceSpcvxEnum,SeaiceAnalysisEnum,finiteelement,0);
+	IoModelToConstraintsx(constraints,iomodel,SeaiceSpcvyEnum,SeaiceAnalysisEnum,finiteelement,1);
+}/*}}}*/
+void SeaiceAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+	return;
+}/*}}}*/
+
+/*Finite Element Analysis*/
+void           SeaiceAnalysis::Core(FemModel* femmodel){/*{{{*/
+	_error_("not implemented");
+}/*}}}*/
+ElementVector* SeaiceAnalysis::CreateDVector(Element* element){/*{{{*/
+	/*Default, return NULL*/
+	return NULL;
+}/*}}}*/
+ElementMatrix* SeaiceAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
+_error_("Not implemented");
+}/*}}}*/
+ElementMatrix* SeaiceAnalysis::CreateKMatrix(Element* element){/*{{{*/
+
+	/* Check if there is ice in this element */
+	Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
+	if(concentration_input->Max()==0.) return NULL;
+
+	/*Intermediaries */
+	IssmDouble  Jdet,D_scalar,dt;
+	IssmDouble  C[3][3];
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+	int numdof   = numnodes*2;
+
+	/*Initialize Element vector and other vectors*/
+	ElementMatrix* Ke      = element->NewElementMatrix();
+	IssmDouble*    B       = xNew<IssmDouble>(3*numdof);
+	IssmDouble*    M       = xNew<IssmDouble>(numdof);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+
+		/*Create first part of the stiffness matrix*/
+		this->GetM(M,element,xyz_list,gauss);
+		D_scalar=gauss->weight*Jdet;
+		TripleMultiply(M,1,numdof,1,
+					&D_scalar,1,1,0,
+					M,1,numdof,0,
+					&Ke->values[0],1);
+
+		/*Create Elasitc part*/
+		this->GetB(B,element,xyz_list,gauss);
+		this->CreateCTensor(&C[0][0],element,gauss);
+		for(int i=0;i<3;i++){
+			for(int j=0;j<3;j++){
+				C[i][j] = gauss->weight*Jdet*dt*C[i][j];
+			}
+		}
+		TripleMultiply(B,3,numdof,1,
+					&C[0][0],3,3,0,
+					B,3,numdof,0,
+					&Ke->values[0],1);
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(B);
+	delete gauss;
+	return Ke;
+}
+/*}}}*/
+ElementVector* SeaiceAnalysis::CreatePVector(Element* element){/*{{{*/
+
+	/* Check if there is ice in this element */
+	Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
+	if(concentration_input->Max()==0.) return NULL;
+
+	/*Intermediaries */
+	IssmDouble  concentration,thickness,dt,Jdet,vx,vy;
+	IssmDouble  rho_air,air_coef,air_lin_drag_coef,air_quad_drag_coef;
+	IssmDouble  windvx,windvy,wind_vnorm;
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+	int numdof   = 2*numnodes;
+
+	/*Initialize Element vector and vectors*/
+	ElementVector* pe    = element->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	element->FindParam(&rho_air,SurfaceforcingsRhoAirEnum);
+	element->FindParam(&air_coef,SurfaceforcingsAirCoefEnum);
+	element->FindParam(&air_lin_drag_coef,SurfaceforcingsAirLinDragCoefEnum);
+	element->FindParam(&air_quad_drag_coef,SurfaceforcingsAirQuadDragCoefEnum);
+	Input* vx_input            = element->GetInput(VxEnum);                    _assert_(vx_input);
+	Input* vy_input            = element->GetInput(VyEnum);                    _assert_(vy_input);
+	Input* windvx_input        = element->GetInput(SurfaceforcingsWindVxEnum); _assert_(windvx_input);
+	Input* windvy_input        = element->GetInput(SurfaceforcingsWindVyEnum); _assert_(windvy_input);
+	Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);       _assert_(thickness_input);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis, gauss);
+
+		thickness_input->GetInputValue(&thickness,gauss);
+		concentration_input->GetInputValue(&concentration,gauss);
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		windvx_input->GetInputValue(&windvx,gauss);
+		windvy_input->GetInputValue(&windvy,gauss);
+
+		/*Atmespheric forcing*/
+		wind_vnorm = sqrt(windvx*windvx + windvy*windvy);
+
+		for(int i=0;i<numnodes;i++){
+			pe->values[i*2+0] += dt*concentration*air_coef*rho_air*(air_lin_drag_coef+air_quad_drag_coef*wind_vnorm)*windvx*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1] += dt*concentration*air_coef*rho_air*(air_lin_drag_coef+air_quad_drag_coef*wind_vnorm)*windvy*Jdet*gauss->weight*basis[i];
+		}
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
+}/*}}}*/
+void SeaiceAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+	   _error_("not implemented yet");
+}/*}}}*/
+void SeaiceAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+	_error_("Not implemented yet");
+}/*}}}*/
+void SeaiceAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	int*        doflist=NULL;
+	IssmDouble* xyz_list=NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+	int numdof   = numnodes*2;
+
+	/*Fetch dof list and allocate solution vectors*/
+	element->GetDofList(&doflist,SSAApproximationEnum,GsetEnum);
+	IssmDouble* values    = xNew<IssmDouble>(numdof);
+	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
+
+	/*Use the dof list to index into the solution vector: */
+	for(int i=0;i<numdof;i++) values[i]=solution[doflist[i]];
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+	for(int i=0;i<numnodes;i++){
+		vx[i]=values[i*2+0];
+		vy[i]=values[i*2+1];
+		if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
+	}
+
+	/*Compute vel*/
+	for(int i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
+
+	/*Add vx and vy as inputs to the tria element: */
+	element->AddInput(VxEnum,vx,element->GetElementType());
+	element->AddInput(VyEnum,vy,element->GetElementType());
+	element->AddInput(VelEnum,vel,element->GetElementType());
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(vel);
+	xDelete<IssmDouble>(vy);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(values);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<int>(doflist);
+}/*}}}*/
+void SeaiceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+	/*Default, do nothing*/
+	return;
+}/*}}}*/
+
+/*Intermediary functions*/
+void SeaiceAnalysis::CreateCTensor(IssmDouble* C,Element* element,Gauss* gauss){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble D[3][3];
+	IssmDouble damage,thickness,concentration;
+
+	/*Get Material parameters*/
+	IssmDouble poisson          = element->GetMaterialParameter(MaterialsPoissonEnum);
+	IssmDouble young_modulus    = element->GetMaterialParameter(MaterialsYoungModulusEnum);
+	IssmDouble ridging_exponent = element->GetMaterialParameter(MaterialsRidgingExponentEnum);
+
+	/*Get damage input at this location*/
+	Input* damage_input        = element->GetInput(MaterialsDamageEnum);     _assert_(damage_input);
+	Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);     _assert_(thickness_input);
+	Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
+	damage_input->GetInputValue(&damage,gauss);
+	thickness_input->GetInputValue(&thickness,gauss);
+	concentration_input->GetInputValue(&concentration,gauss);
+
+	/*Build D tensor once for all*/
+	D[0][0] = 1.;      D[0][1] = poisson; D[0][2] = 0.;
+	D[1][0] = poisson; D[1][1] = 1.;      D[1][2] = 0.;
+	D[2][0] = 0.;      D[2][1] = 0;       D[2][2] = (1.-poisson)/2.;
+	for(int i=0;i<3;i++){
+		for(int j=0;j<3;j++){
+			D[i][j] = D[i][j]/(1-poisson*poisson);
+		}
+	}
+
+	/*Get f(A), E*/
+	IssmDouble f = exp(ridging_exponent*(1.-concentration));
+	IssmDouble E = young_modulus*(1.-damage)*f;
+
+	/*Let's build C = hE*D */
+	for(int i=0;i<3;i++){
+		for(int j=0;j<3;j++){
+			C[i*3+j] = thickness*E*D[i][j];
+		}
+	}
+}/*}}}*/
+void SeaiceAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+	 * For node i, Bi can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi=[ dN/dx           0    ] 
+	 *          [   0           dN/dy  ]      
+	 *          [ 1/2*dN/dy  1/2*dN/dx ]     
+	 * where N is the basis function for node i.
+	 *
+	 * We assume B has been allocated already           */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Build B: */
+	for(int i=0;i<numnodes;i++){
+		B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
+		B[2*numnodes*0+2*i+1] = 0.;
+		B[2*numnodes*1+2*i+0] = 0.;
+		B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
+		B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
+		B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
+}/*}}}*/
+void SeaiceAnalysis::GetM(IssmDouble* M,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute M  matrix. M=[M1 M2 M3] where Mi is of size 3*NDOF2. 
+	 * For node i, Mi can be expressed in the actual coordinate system
+	 * by: 
+	 *       Mi=[ N  N ] 
+	 * where N is the basis function for node i.
+	 *
+	 * We assume M has been allocated already           */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	element->NodalFunctions(basis,gauss);
+
+	/*Build M: */
+	for(int i=0;i<numnodes;i++){
+		M[2*i+0] = basis[i];
+		M[2*i+1] = basis[i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
+}/*}}}*/
Index: ../trunk-jpl/src/c/analyses/analyses.h
===================================================================
--- ../trunk-jpl/src/c/analyses/analyses.h	(revision 18491)
+++ ../trunk-jpl/src/c/analyses/analyses.h	(revision 18492)
@@ -33,6 +33,7 @@
 #include "./MeshdeformationAnalysis.h"
 #include "./SmoothedSurfaceSlopeXAnalysis.h"
 #include "./SmoothedSurfaceSlopeYAnalysis.h"
+#include "./SeaiceAnalysis.h"
 #include "./StressbalanceAnalysis.h"
 #include "./StressbalanceSIAAnalysis.h"
 #include "./StressbalanceVerticalAnalysis.h"
Index: ../trunk-jpl/src/c/bamg/Mesh.cpp
===================================================================
--- ../trunk-jpl/src/c/bamg/Mesh.cpp	(revision 18491)
+++ ../trunk-jpl/src/c/bamg/Mesh.cpp	(revision 18492)
@@ -241,7 +241,7 @@
 			if (Gh.NbRef>0) Gh.NbRef--;
 			else if (Gh.NbRef==0) delete &Gh;
 		}
-		if(&BTh && (&BTh != this)){
+		if(&BTh && (&BTh != this) && false){
 			if (BTh.NbRef>0) BTh.NbRef--;
 			else if (BTh.NbRef==0) delete &BTh;
 		}
@@ -4504,308 +4504,295 @@
 		/*Get options*/
 		int verbose=bamgopts->verbose;
 
-		Gh.NbRef++;// add a ref to Gh
+		/*Intermediaries*/
+		int                i,k;
+		int                nbcurves    = 0;
+		int                NbNewPoints,NbEdgeCurve;
+		double             lcurve,lstep,s;
+		const int          MaxSubEdge  = 10;
 
-		/************************************************************************* 
-		 * method in 2 steps
-		 * 1 - compute the number of new edges to allocate
-		 * 2 - construct the edges
-		 * remark:
-		 * in this part we suppose to have a background mesh with the same geometry 
-		 * 
-		 * To construct the discretization of the new mesh we have to 
-		 * rediscretize the boundary of background Mesh 
-		 * because we have only the pointeur from the background mesh to the geometry.
-		 * We need the abcisse of the background mesh vertices on geometry
-		 * so a vertex is 
-		 * 0 on GeomVertex ;
-		 * 1 on GeomEdge + abcisse
-		 * 2 internal 
-		 *************************************************************************/
+		R2          AB;
+		GeomVertex *a, *b;
+		BamgVertex *va,*vb;
+		GeomEdge   *e;
 
-		//Check that background mesh and current mesh do have the same geometry
-		_assert_(&BTh.Gh==&Gh);
-		BTh.NbRef++; // add a ref to BackGround Mesh
+		// add a ref to GH to make sure that it is not destroyed by mistake
+		Gh.NbRef++;
 
-		//Initialize new mesh
-		BTh.SetVertexFieldOn();
-		int* bcurve = new int[Gh.nbcurves]; // 
+		//build background mesh flag (1 if background, else 0)
+		bool background=(&BTh != this);
 
-		/* There are 2 ways to make the loop 
-		 * 1) on the geometry 
-		 * 2) on the background mesh
-		 *  if you do the loop on geometry, we don't have the pointeur on background,
-		 *  and if you do the loop in background we have the pointeur on geometry
-		 * so do the walk on  background */
+		/*Build VerticesOnGeomVertex*/
 
-		NbVerticesOnGeomVertex=0;
-		NbVerticesOnGeomEdge=0;
+		//Compute the number of geometrical vertices that we are going to use to mesh
+		for (i=0;i<Gh.nbv;i++){
+			if (Gh[i].Required()) NbVerticesOnGeomVertex++;
+		}
+		//allocate
+		VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];  
+		if(NbVerticesOnGeomVertex >= maxnbv) _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
+		_assert_(nbv==0);
+		//Build VerticesOnGeomVertex
+		for (i=0;i<Gh.nbv;i++){
+			/* Add vertex only if required*/
+			if (Gh[i].Required()) {//Gh  vertices Required
 
-		/*STEP 1 copy of Required vertices*/
+				//Add the vertex
+				_assert_(nbv<maxnbv);
+				vertices[nbv]=Gh[i];
 
-		int i; 
-		for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++;
-		printf("\n");
-		if(NbVerticesOnGeomVertex >= maxnbv){
-			_error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
-		}
+				//Add pointer from geometry (Gh) to vertex from mesh (Th)
+				Gh[i].MeshVertexHook=vertices+nbv;
 
-		VerticesOnGeomVertex = new VertexOnGeom[  NbVerticesOnGeomVertex];
-		VertexOnBThVertex    = new VertexOnVertex[NbVerticesOnGeomVertex];
+				//Build VerticesOnGeomVertex for current point
+				VerticesOnGeomVertex[nbv]=VertexOnGeom(*Gh[i].MeshVertexHook,Gh[i]);
 
-		//At this point there is NO vertex but vertices should have been allocated by Init
-		_assert_(vertices);
-		for (i=0;i<Gh.nbv;i++){
-			if (Gh[i].Required()) {//Gh vertices Required
-				vertices[nbv]  =Gh[i];
-				vertices[nbv].i=I2(0,0);
-				Gh[i].MeshVertexHook = vertices + nbv;// save Geom -> Th
-				VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]);
+				//nbv increment
 				nbv++;
 			}
-			else Gh[i].MeshVertexHook=0;
-		} 
-		for (i=0;i<BTh.NbVerticesOnGeomVertex;i++){ 
-			VertexOnGeom &vog=BTh.VerticesOnGeomVertex[i];
-			if (vog.IsRequiredVertex()){
-				GeomVertex* gv=vog;
-				BamgVertex *bv = vog;
-				_assert_(gv->MeshVertexHook); // use of Geom -> Th
-				VertexOnBThVertex[NbVertexOnBThVertex++]=VertexOnVertex(gv->MeshVertexHook,bv);
-				gv->MeshVertexHook->m = bv->m; // for taking the metric of the background mesh
-			}
 		}
-		_assert_(NbVertexOnBThVertex==NbVerticesOnGeomVertex); /*This might be due to MaxCornerAngle too small*/
 
-		/*STEP 2: reseed boundary edges*/
+		/*Build VerticesOnGeomEdge*/
 
-		//  find the begining of the curve in BTh
-		Gh.UnMarkEdges();	
-		int bfind=0;
-		for (int i=0;i<Gh.nbcurves;i++) bcurve[i]=-1; 
+		//check that edges is still empty (Init)
+		_assert_(!edges);
 
-		/*Loop over the backgrounf mesh BTh edges*/
-		for (int iedge=0;iedge<BTh.nbe;iedge++){      
-			Edge &ei = BTh.edges[iedge];
+		/* Now we are going to create the first edges corresponding
+		 * to the one present in the geometry provided.
+		 * We proceed in 2 steps
+		 *  -step 0: we count all the edges
+		 *           we allocate the number of edges at the end of step 0
+		 *  -step 1: the edges are created */
+		for (int step=0;step<2;step++){
 
-			/*Loop over the 2 vertices of the current edge*/
-			for(int je=0;je<2;je++){
+			//initialize number of edges and number of edges max
+			long nbex=0;
+			nbe=0;
+			long NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge;
+			Gh.UnMarkEdges();	
+			nbcurves=0;
 
-				/* If one of the vertex is required we are in a new curve*/
-				if (ei[je].GeomEdgeHook->IsRequiredVertex()){ 
+			//go through the edges of the geometry
+			for (i=0;i<Gh.nbe;i++){
 
-					/*Get curve number*/
-					int nc=ei.GeomEdgeHook->CurveNumber;
+				//ei = current Geometrical edge
+				GeomEdge &ei=Gh.edges[i];   
 
-					//_printf_("Dealing with curve number " << nc << "\n");
-					//_printf_("edge on geometry is same as GhCurve? " << (ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge)?"yes":"no\n");
-					//if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge){
-					//	_printf_("Do we have the right extremity? curve first vertex -> " << ((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex])?"yes":"no\n");
-					//	_printf_("Do we have the right extremity? curve last  vertex -> " << ((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex])?"yes":"no\n");
-					//}
-					//BUG FIX from original bamg
-					/*Check that we are on the same edge and right vertex (0 or 1) */
-					if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge  && (GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex]){
-						bcurve[nc]=iedge*2+je;
-						bfind++;	
-					}
-					else if ((ei.GeomEdgeHook==Gh.curves[nc].LastEdge  && (GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex]) && bcurve[nc]==-1){
-						bcurve[nc]=iedge*2+je;
-						bfind++;	
-					}
-				}
-			}
-		} 
-		if (bfind!=Gh.nbcurves){
-			delete [] bcurve;
-			_error_("problem generating number of curves (" << Gh.nbcurves << " found in the geometry but " << bfind << " curve found in the mesh)");
-		}
+				//loop over the two vertices of the edge ei
+				for(int j=0;j<2;j++) {
 
-		// method in 2 + 1 step 
-		//  0.0) compute the length and the number of vertex to do allocation
-		//  1.0) recompute the length
-		//  1.1) compute the  vertex 
+					/*Take only required vertices (corner->beginning of a new curve)*/
+					if (!ei.Mark() && ei[j].Required()){ 
 
-		long nbex=0,NbVerticesOnGeomEdgex=0;
-		for (int step=0; step <2;step++){
+						long  nbvend=0;
+						Edge* PreviousNewEdge=NULL;
+						lstep = -1;
 
-			long NbOfNewPoints=0;
-			long NbOfNewEdge=0;
-			long iedge;
-			Gh.UnMarkEdges();	
-			double L=0;
+						/*If Edge is required (do that only once for the 2 vertices)*/
+						if(ei.Required()){
+							if (j==0){
+								//do not create internal points if required (take it as is)
+								if(step==0) nbe++;
+								else{ 
+									e=&ei;
+									a=ei(0);
+									b=ei(1);
 
-			/*Go through all geometrical curve*/
-			for (int icurve=0;icurve<Gh.nbcurves;icurve++){ 
+									//check that edges has been allocated
+									_assert_(edges);
+									edges[nbe].v[0]=a->MeshVertexHook;
+									edges[nbe].v[1]=b->MeshVertexHook;;
+									edges[nbe].ReferenceNumber = e->ReferenceNumber;
+									edges[nbe].GeomEdgeHook = e;
+									edges[nbe].adj[0] = 0;
+									edges[nbe].adj[1] = 0;
+									nbe++;
+								}
+							}
+						}
 
-				/*Get edge and vertex (index) of background mesh on this curve*/
-				iedge=bcurve[icurve]/2;
-				int jedge=bcurve[icurve]%2;
+						/*If Edge is not required: we are on a curve*/
+						else {
+							for (int kstep=0;kstep<=step;kstep++){
+								//kstep=0, compute number of edges (discretize curve)
+								//kstep=1  create the points and edge
+								PreviousNewEdge=0;
+								NbNewPoints=0;
+								NbEdgeCurve=0;
+								if (nbvend>=maxnbv) _error_("maximum number of vertices too low! Check the domain outline or increase maxnbv");
+								lcurve =0;
+								s = lstep; //-1 initially, then length of each sub edge
 
-				/*Get edge of Bth with index iedge*/
-				Edge &ei = BTh.edges[iedge];
+								/*reminder: i = edge number, j=[0;1] vertex index in edge*/
+								k=j;            // k = vertex index in edge (0 or 1)
+								e=&ei;          // e = reference of current edge
+								a=ei(k);        // a = pointer toward the kth vertex of the current edge
+								va = a->MeshVertexHook; // va = pointer toward mesh vertex associated
+								e->SetMark();   // Mark edge
 
-				/*Initialize variables*/
-				double Lstep=0;             // step between two points   (phase==1) 
-				long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1) 
+								/*Loop until we reach the end of the curve*/
+								for(;;){ 
+									k = 1-k;            // other vertx index of the curve
+									b = (*e)(k);        // b = pointer toward the other vertex of the current edge
+									AB= b->r - a->r;   // AB = vector of the current edge
+									Metric MA = background ? BTh.MetricAt(a->r) :a->m ;  //Get metric associated to A
+									Metric MB = background ? BTh.MetricAt(b->r) :b->m ;  //Get metric associated to B
+									double ledge = (MA(AB) + MB(AB))/2;                  //Get edge length in metric
 
-				/*Do phase 0 to step*/
-				for(int phase=0;phase<=step;phase++){
+									/* We are now creating the mesh edges from the geometrical edge selected above.
+									 * The edge will be divided according to the metric previously computed and cannot
+									 * be divided more than 10 times (MaxSubEdge). */
 
-					/*Current curve pointer*/
-					Curve *curve= Gh.curves+icurve;
+									//By default, there is only one subedge that is the geometrical edge itself
+									int NbSubEdge = 1;
 
-					/*Get index of current curve*/
-					int icurveequi= Gh.GetId(curve);
+									//initialize lSubEdge, holding the length of each subedge (cannot be higher than 10)
+									double lSubEdge[MaxSubEdge];
 
-					/*For phase 0, check that we are at the begining of the curve only*/
-					if(phase==0 &&  icurveequi!=icurve)  continue;
+									//Build Subedges according to the edge length
+									if (ledge < 1.5){
+										//if ledge < 1.5 (between one and 2), take the edge as is
+										lSubEdge[0] = ledge;
+									}
+									//else, divide the edge
+									else {
+										//compute number of subedges (division of the edge), Maximum is 10
+										NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5));
+										/*Now, we are going to divide the edge according to the metric.
+										 * Get segment by sement along the edge.
+										 * Build lSubEdge, which holds the distance between the first vertex
+										 * of the edge and the next point on the edge according to the 
+										 * discretization (each SubEdge is AB)*/
+										R2 A,B;
+										A=a->r;
+										Metric MAs=MA,MBs;
+										ledge=0; 
+										double x =0, xstep= 1./NbSubEdge;
+										for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) {
+											x += xstep;
+											B =  e->F(k ? x : 1-x);
+											MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB);
+											AB = A-B;
+											lSubEdge[kk]=(ledge+=(MAs(AB)+MBs(AB))/2);
+										}
+									}
 
-					int   k0=jedge,k1;
-					Edge* pe=  BTh.edges+iedge;
-					int   iedgeequi=bcurve[icurveequi]/2;
-					int   jedgeequi=bcurve[icurveequi]%2;
+									double lcurveb = lcurve+ledge;
 
-					int k0equi=jedgeequi,k1equi;		  
-					Edge * peequi= BTh.edges+iedgeequi;
-					GeomEdge *ongequi = peequi->GeomEdgeHook;
+									/*Now, create corresponding points*/
+									while(s>=lcurve && s<=lcurveb && nbv<nbvend){
 
-					double sNew=Lstep;// abscisse of the new points (phase==1) 
-					L=0;// length of the curve
-					long i=0;// index of new points on the curve
-					GeomVertex * GA0 = *(*peequi)[k0equi].GeomEdgeHook;
-					BamgVertex *A0;
-					A0 = GA0->MeshVertexHook;  // the vertex in new mesh
-					BamgVertex *A1;
-					VertexOnGeom *GA1;
-					Edge* PreviousNewEdge = 0;
+										/*Schematic of current curve
+										 *
+										 *  a                   vb                  b          // vertex
+										 *  0              ll0     ll1              ledge      // length from a
+										 *  + --- + - ... - + --S-- + --- + - ... - +          // where is S
+										 *  0              kk0     kk1              NbSubEdge  // Sub edge index
+										 *
+										 */
 
-					// New Curve phase 
-					_assert_(A0-vertices>=0 && A0-vertices<nbv);
-					if(ongequi->Required()){
-						GeomVertex *GA1 = *(*peequi)[1-k0equi].GeomEdgeHook;
-						A1 = GA1->MeshVertexHook;  //
-					}       
-					else {
-						for(;;){
-							Edge &ee=*pe; 
-							Edge &eeequi=*peequi; 
-							k1 = 1-k0; // next vertex of the edge 
-							k1equi= 1 - k0equi;
-							_assert_(pe && ee.GeomEdgeHook);
-							ee.GeomEdgeHook->SetMark();
-							BamgVertex & v0=ee[0], & v1=ee[1];
-							R2 AB=(R2)v1-(R2)v0;
-							double L0=L,LAB;
-							LAB=LengthInterpole(v0.m,v1.m,AB);
-							L+= LAB;
+										double ss = s-lcurve;
 
-							if (phase){
-								// computation of the new points for the given curve
-								while ((i!=NbCreatePointOnCurve) && sNew<=L) { 
+										/*Find the SubEdge containing ss using Dichotomy*/
+										int kk0=-1,kk1=NbSubEdge-1,kkk;
+										double ll0=0,ll1=ledge,llk;
+										while (kk1-kk0>1){
+											if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2]))
+											 kk1=kkk,ll1=llk;
+											else
+											 kk0=kkk,ll0=llk;
+										}
+										_assert_(kk1!=kk0);
 
-									//some checks
-									_assert_(sNew>=L0);
-									_assert_(LAB);
-									_assert_(vertices && nbv<maxnbv);
-									_assert_(edges && nbe<nbex);
-									_assert_(VerticesOnGeomEdge && NbVerticesOnGeomEdge<NbVerticesOnGeomEdgex);
+										/*Curvilinear coordinate in [0 1] of ss in current edge*/
+										// WARNING: This is what we would do
+										// ssa = (ss-ll0)/(ll1-ll0);
+										// aa = (kk0+ssa)/NbSubEdge
+										// This is what Bamg does:
+										double sbb = (ss-ll0)/(ll1-ll0);
+										/*Curvilinear coordinate in [0 1] of ss in current curve*/
+										double bb = (kk1+sbb)/NbSubEdge;
+										double aa = 1-bb;
 
-									// new vertex on edge
-									A1=vertices+nbv++;
-									GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge;
-									Edge* e = edges + nbe++;
-									double se= (sNew-L0)/LAB;
-									if (se<0 || se>=1.000000001){
-										_error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]");
+										// new vertex on edge
+										vb = &vertices[nbv++];
+										vb->m = Metric(aa,a->m,bb,b->m);
+										vb->ReferenceNumber = e->ReferenceNumber;
+										vb->DirOfSearch =NoDirOfSearch;
+										double abcisse = k ? bb : aa;
+										vb->r =  e->F(abcisse);
+										VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse);        
+
+										// to take into account the direction of the edge
+										s += lstep;
+										edges[nbe].v[0]=va;
+										edges[nbe].v[1]=vb;
+										edges[nbe].ReferenceNumber =e->ReferenceNumber;
+										edges[nbe].GeomEdgeHook = e;
+										edges[nbe].adj[0] = PreviousNewEdge;
+										if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe];
+										PreviousNewEdge=edges+nbe;
+										nbe++;
+										va = vb;
 									}
-									se = abscisseInterpole(v0.m,v1.m,AB,se,1);
-									if (se<0 || se>1){
-										_error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]");
-									}
-									se = k1         ? se : 1. - se;
-									se = k1==k1equi ? se : 1. - se;
-									VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save 
-									ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1); 
-									A1->ReferenceNumber = eeequi.ReferenceNumber;
-									A1->DirOfSearch =NoDirOfSearch;
-									e->GeomEdgeHook = ongequi;
-									e->v[0]=A0;
-									e->v[1]=A1;
-									e->ReferenceNumber = eeequi.ReferenceNumber;
-									e->adj[0]=PreviousNewEdge;
 
-									if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
-									PreviousNewEdge=e;
-									A0=A1;
-									sNew += Lstep;
-									if (++i== NbCreatePointOnCurve) break;
+									/*We just added one edge to the curve: Go to the next one*/
+									lcurve = lcurveb;
+									e->SetMark();
+									a=b;
+
+									/*If b is required, we are on a new curve->break*/
+									if (b->Required()) break;
+									int kprev=k;
+									k = e->AdjVertexIndex[kprev];// next vertices
+									e = e->Adj[kprev];
+									_assert_(e);
+								}// for(;;)
+								vb = b->MeshVertexHook;
+
+								/*Number of edges in the last disretized curve*/
+								NbEdgeCurve = Max((long) (lcurve +0.5), (long) 1);
+								/*Number of internal vertices in the last disretized curve*/
+								NbNewPoints = NbEdgeCurve-1;
+								if(!kstep){
+									NbVerticesOnGeomEdge0 += NbNewPoints;
+									nbcurves++;
 								}
+								nbvend=nbv+NbNewPoints; 
+								lstep = lcurve / NbEdgeCurve; //approximately one
+							}// end of curve --
+							if (edges) { // last edges of the curves 
+								edges[nbe].v[0]=va;
+								edges[nbe].v[1]=vb;
+								edges[nbe].ReferenceNumber = e->ReferenceNumber;
+								edges[nbe].GeomEdgeHook = e;
+								edges[nbe].adj[0] = PreviousNewEdge;
+								edges[nbe].adj[1] = 0;
+								if(PreviousNewEdge) PreviousNewEdge->adj[1] = & edges[nbe];
+								nbe++;
 							}
-
-							//some checks
-							_assert_(ee.GeomEdgeHook->CurveNumber==ei.GeomEdgeHook->CurveNumber);
-							if (ee[k1].GeomEdgeHook->IsRequiredVertex()) {
-								_assert_(eeequi[k1equi].GeomEdgeHook->IsRequiredVertex());
-								GeomVertex * GA1 = *eeequi[k1equi].GeomEdgeHook;
-								A1=GA1->MeshVertexHook;// the vertex in new mesh
-								_assert_(A1-vertices>=0 && A1-vertices<nbv);
-								break;
-							}
-							if (!ee.adj[k1]) {
-								_error_("adj edge " << BTh.GetId(ee) << ", nbe=" << nbe << ", Gh.vertices=" << Gh.vertices);
-							}
-							pe = ee.adj[k1]; // next edge
-							k0 = pe->Intersection(ee); 
-							peequi= eeequi.adj[k1equi];  // next edge
-							k0equi=peequi->Intersection(eeequi);            
-						}// for(;;) end of the curve
+							else nbe += NbEdgeCurve;
+						} // end on  curve ---
 					}
+				}
+			} // for (i=0;i<nbe;i++)
+			if(!step) {
+				_assert_(!edges);
+				_assert_(!VerticesOnGeomEdge);
 
-					if (phase){ // construction of the last edge
-						Edge* e=edges + nbe++;
-						e->GeomEdgeHook  = ongequi;
-						e->v[0]=A0;
-						e->v[1]=A1;
-						e->ReferenceNumber = peequi->ReferenceNumber;
-						e->adj[0]=PreviousNewEdge;
-						e->adj[1]=0;
-						if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
-						PreviousNewEdge = e;
+				edges = new Edge[nbex=nbe];
+				if(NbVerticesOnGeomEdge0) VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0];
 
-						_assert_(i==NbCreatePointOnCurve);
-					}
-
-					if (!phase)  { // 
-						long NbSegOnCurve = Max((long)(L+0.5),(long) 1);// nb of seg
-						Lstep = L/NbSegOnCurve; 
-						NbCreatePointOnCurve = NbSegOnCurve-1;
-						NbOfNewEdge += NbSegOnCurve;
-						NbOfNewPoints += NbCreatePointOnCurve;
-					}
-				}
-			}//  end of curve loop 
-
-			//Allocate memory
-			if(step==0){
-				if(nbv+NbOfNewPoints > maxnbv) {
-					_error_("too many vertices on geometry: " << nbv+NbOfNewPoints << " >= " << maxnbv);
-				}
-				edges = new Edge[NbOfNewEdge];
-				nbex = NbOfNewEdge;
-				if(NbOfNewPoints) {
-					VerticesOnGeomEdge    = new VertexOnGeom[NbOfNewPoints];
-					NbVertexOnBThEdge     = NbOfNewPoints;
-					VertexOnBThEdge       = new  VertexOnEdge[NbOfNewPoints];
-					NbVerticesOnGeomEdgex = NbOfNewPoints;
-				}
-				NbOfNewPoints =0;
-				NbOfNewEdge = 0;
+				// do the vertex on a geometrical vertex
+				_assert_(VerticesOnGeomEdge || NbVerticesOnGeomEdge0==0);
+				NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge;       
 			}
+			else{
+				_assert_(NbVerticesOnGeomEdge==NbVerticesOnGeomEdge0);
+			}
 		}
-		_assert_(nbe!=0);
-		delete [] bcurve;
 
+
 		//Insert points inside existing triangles
 		if (verbose>4) _printf_("      -- current number of vertices = " << nbv << "\n");
 		if (verbose>3) _printf_("      Creating initial Constrained Delaunay Triangulation...\n");
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18491)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18492)
@@ -386,6 +386,8 @@
 	UzawaPressureAnalysisEnum,
 	GiaSolutionEnum,
 	GiaAnalysisEnum,
+	SeaiceSolutionEnum,
+	SeaiceAnalysisEnum,
 	MeshdeformationSolutionEnum,
 	MeshdeformationAnalysisEnum,
 	LevelsetAnalysisEnum,
@@ -741,6 +743,39 @@
 	LevelsetfunctionSlopeYEnum,
 	LevelsetfunctionPicardEnum,
 	/*}}}*/
+	/*Sea Ice {{{*/
+	SeaiceatmEnum,
+	SeaiceoceanEnum,
+	SeaiceThicknessEnum,
+	SeaiceConcentrationEnum,
+	SeaiceSpcvxEnum,
+	SeaiceSpcvyEnum,
+	BasalforcingsRhoOceanEnum,
+	BasalforcingsOceanCoefEnum,
+	BasalforcingsOceanLinDragCoefEnum,
+	BasalforcingsOceanQuadDragCoefEnum,
+	BasalforcingsOceanTurningAngleEnum,
+	BasalforcingsOceanSshEnum,
+	BasalforcingsOceanVxEnum,
+	BasalforcingsOceanVyEnum,
+	SurfaceforcingsRhoAirEnum,
+	SurfaceforcingsAirCoefEnum,
+	SurfaceforcingsAirLinDragCoefEnum,
+	SurfaceforcingsAirQuadDragCoefEnum,
+	SurfaceforcingsWindVxEnum,
+	SurfaceforcingsWindVyEnum,
+	MatseaiceEnum,
+	MaterialsPoissonEnum,
+	MaterialsYoungModulusEnum,
+	MaterialsDamageEnum,
+	MaterialsRidgingExponentEnum,
+	MaterialsCohesionEnum,
+	MaterialsInternalFrictionCoefEnum,
+	MaterialsCompressionCoefEnum,
+	MaterialsTractionCoefEnum,
+	VxStarEnum,
+	VyStarEnum,
+	/*}}}*/
 	MaximumNumberOfDefinitionsEnum
 };
 
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18491)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18492)
@@ -389,6 +389,8 @@
 		case UzawaPressureAnalysisEnum : return "UzawaPressureAnalysis";
 		case GiaSolutionEnum : return "GiaSolution";
 		case GiaAnalysisEnum : return "GiaAnalysis";
+		case SeaiceSolutionEnum : return "SeaiceSolution";
+		case SeaiceAnalysisEnum : return "SeaiceAnalysis";
 		case MeshdeformationSolutionEnum : return "MeshdeformationSolution";
 		case MeshdeformationAnalysisEnum : return "MeshdeformationAnalysis";
 		case LevelsetAnalysisEnum : return "LevelsetAnalysis";
@@ -701,6 +703,37 @@
 		case LevelsetfunctionSlopeXEnum : return "LevelsetfunctionSlopeX";
 		case LevelsetfunctionSlopeYEnum : return "LevelsetfunctionSlopeY";
 		case LevelsetfunctionPicardEnum : return "LevelsetfunctionPicard";
+		case SeaiceatmEnum : return "Seaiceatm";
+		case SeaiceoceanEnum : return "Seaiceocean";
+		case SeaiceThicknessEnum : return "SeaiceThickness";
+		case SeaiceConcentrationEnum : return "SeaiceConcentration";
+		case SeaiceSpcvxEnum : return "SeaiceSpcvx";
+		case SeaiceSpcvyEnum : return "SeaiceSpcvy";
+		case BasalforcingsRhoOceanEnum : return "BasalforcingsRhoOcean";
+		case BasalforcingsOceanCoefEnum : return "BasalforcingsOceanCoef";
+		case BasalforcingsOceanLinDragCoefEnum : return "BasalforcingsOceanLinDragCoef";
+		case BasalforcingsOceanQuadDragCoefEnum : return "BasalforcingsOceanQuadDragCoef";
+		case BasalforcingsOceanTurningAngleEnum : return "BasalforcingsOceanTurningAngle";
+		case BasalforcingsOceanSshEnum : return "BasalforcingsOceanSsh";
+		case BasalforcingsOceanVxEnum : return "BasalforcingsOceanVx";
+		case BasalforcingsOceanVyEnum : return "BasalforcingsOceanVy";
+		case SurfaceforcingsRhoAirEnum : return "SurfaceforcingsRhoAir";
+		case SurfaceforcingsAirCoefEnum : return "SurfaceforcingsAirCoef";
+		case SurfaceforcingsAirLinDragCoefEnum : return "SurfaceforcingsAirLinDragCoef";
+		case SurfaceforcingsAirQuadDragCoefEnum : return "SurfaceforcingsAirQuadDragCoef";
+		case SurfaceforcingsWindVxEnum : return "SurfaceforcingsWindVx";
+		case SurfaceforcingsWindVyEnum : return "SurfaceforcingsWindVy";
+		case MatseaiceEnum : return "Matseaice";
+		case MaterialsPoissonEnum : return "MaterialsPoisson";
+		case MaterialsYoungModulusEnum : return "MaterialsYoungModulus";
+		case MaterialsDamageEnum : return "MaterialsDamage";
+		case MaterialsRidgingExponentEnum : return "MaterialsRidgingExponent";
+		case MaterialsCohesionEnum : return "MaterialsCohesion";
+		case MaterialsInternalFrictionCoefEnum : return "MaterialsInternalFrictionCoef";
+		case MaterialsCompressionCoefEnum : return "MaterialsCompressionCoef";
+		case MaterialsTractionCoefEnum : return "MaterialsTractionCoef";
+		case VxStarEnum : return "VxStar";
+		case VyStarEnum : return "VyStar";
 		case MaximumNumberOfDefinitionsEnum : return "MaximumNumberOfDefinitions";
 		default : return "unknown";
 
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18491)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18492)
@@ -398,6 +398,8 @@
 	      else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum;
 	      else if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum;
 	      else if (strcmp(name,"GiaAnalysis")==0) return GiaAnalysisEnum;
+	      else if (strcmp(name,"SeaiceSolution")==0) return SeaiceSolutionEnum;
+	      else if (strcmp(name,"SeaiceAnalysis")==0) return SeaiceAnalysisEnum;
 	      else if (strcmp(name,"MeshdeformationSolution")==0) return MeshdeformationSolutionEnum;
 	      else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum;
 	      else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
@@ -503,12 +505,12 @@
 	      else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
 	      else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
 	      else if (strcmp(name,"Misfit")==0) return MisfitEnum;
-	      else if (strcmp(name,"Pressure")==0) return PressureEnum;
-	      else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
+	      if (strcmp(name,"Pressure")==0) return PressureEnum;
+	      else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
+	      else if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
 	      else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum;
 	      else if (strcmp(name,"SegmentOnIceShelf")==0) return SegmentOnIceShelfEnum;
 	      else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
@@ -626,12 +628,12 @@
 	      else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
 	      else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
 	      else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
-	      else if (strcmp(name,"MinVel")==0) return MinVelEnum;
-	      else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"MinVx")==0) return MinVxEnum;
+	      if (strcmp(name,"MinVel")==0) return MinVelEnum;
+	      else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
+	      else if (strcmp(name,"MinVx")==0) return MinVxEnum;
 	      else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
 	      else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
 	      else if (strcmp(name,"MinVy")==0) return MinVyEnum;
@@ -716,6 +718,37 @@
 	      else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
 	      else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
 	      else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
+	      else if (strcmp(name,"Seaiceatm")==0) return SeaiceatmEnum;
+	      else if (strcmp(name,"Seaiceocean")==0) return SeaiceoceanEnum;
+	      else if (strcmp(name,"SeaiceThickness")==0) return SeaiceThicknessEnum;
+	      else if (strcmp(name,"SeaiceConcentration")==0) return SeaiceConcentrationEnum;
+	      else if (strcmp(name,"SeaiceSpcvx")==0) return SeaiceSpcvxEnum;
+	      else if (strcmp(name,"SeaiceSpcvy")==0) return SeaiceSpcvyEnum;
+	      else if (strcmp(name,"BasalforcingsRhoOcean")==0) return BasalforcingsRhoOceanEnum;
+	      else if (strcmp(name,"BasalforcingsOceanCoef")==0) return BasalforcingsOceanCoefEnum;
+	      else if (strcmp(name,"BasalforcingsOceanLinDragCoef")==0) return BasalforcingsOceanLinDragCoefEnum;
+	      else if (strcmp(name,"BasalforcingsOceanQuadDragCoef")==0) return BasalforcingsOceanQuadDragCoefEnum;
+	      else if (strcmp(name,"BasalforcingsOceanTurningAngle")==0) return BasalforcingsOceanTurningAngleEnum;
+	      else if (strcmp(name,"BasalforcingsOceanSsh")==0) return BasalforcingsOceanSshEnum;
+	      else if (strcmp(name,"BasalforcingsOceanVx")==0) return BasalforcingsOceanVxEnum;
+	      else if (strcmp(name,"BasalforcingsOceanVy")==0) return BasalforcingsOceanVyEnum;
+	      else if (strcmp(name,"SurfaceforcingsRhoAir")==0) return SurfaceforcingsRhoAirEnum;
+	      else if (strcmp(name,"SurfaceforcingsAirCoef")==0) return SurfaceforcingsAirCoefEnum;
+	      else if (strcmp(name,"SurfaceforcingsAirLinDragCoef")==0) return SurfaceforcingsAirLinDragCoefEnum;
+	      else if (strcmp(name,"SurfaceforcingsAirQuadDragCoef")==0) return SurfaceforcingsAirQuadDragCoefEnum;
+	      else if (strcmp(name,"SurfaceforcingsWindVx")==0) return SurfaceforcingsWindVxEnum;
+	      else if (strcmp(name,"SurfaceforcingsWindVy")==0) return SurfaceforcingsWindVyEnum;
+	      else if (strcmp(name,"Matseaice")==0) return MatseaiceEnum;
+	      else if (strcmp(name,"MaterialsPoisson")==0) return MaterialsPoissonEnum;
+	      else if (strcmp(name,"MaterialsYoungModulus")==0) return MaterialsYoungModulusEnum;
+	      else if (strcmp(name,"MaterialsDamage")==0) return MaterialsDamageEnum;
+	      else if (strcmp(name,"MaterialsRidgingExponent")==0) return MaterialsRidgingExponentEnum;
+	      else if (strcmp(name,"MaterialsCohesion")==0) return MaterialsCohesionEnum;
+	      else if (strcmp(name,"MaterialsInternalFrictionCoef")==0) return MaterialsInternalFrictionCoefEnum;
+	      else if (strcmp(name,"MaterialsCompressionCoef")==0) return MaterialsCompressionCoefEnum;
+	      else if (strcmp(name,"MaterialsTractionCoef")==0) return MaterialsTractionCoefEnum;
+	      else if (strcmp(name,"VxStar")==0) return VxStarEnum;
+	      else if (strcmp(name,"VyStar")==0) return VyStarEnum;
 	      else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
          else stage=7;
    }
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 18491)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 18492)
@@ -42,10 +42,7 @@
 	/*Copy some constants from iomodel */
 	parameters->AddObject(iomodel->CopyConstantObject(DomainTypeEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(DomainDimensionEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(MeshElementtypeEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(SettingsOutputFrequencyEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(SteadystateReltolEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(SteadystateMaxiterEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(ConstantsYtsEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(TimesteppingStartTimeEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(TimesteppingFinalTimeEnum));
@@ -56,37 +53,44 @@
 	parameters->AddObject(iomodel->CopyConstantObject(SettingsLowmemEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(DebugProfilingEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(MeshAverageVertexConnectivityEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(ConstantsReferencetemperatureEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(SettingsWaitonlockEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(MeshNumberofelementsEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(MeshNumberofverticesEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(SettingsResultsOnNodesEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(SettingsIoGatherEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(GroundinglineMigrationEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(TransientIsstressbalanceEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(TransientIsmasstransportEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(TransientIsthermalEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(TransientIsgroundinglineEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(TransientIsgiaEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(TransientIslevelsetEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(TransientIsdamageevolutionEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(TransientIshydrologyEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(AutodiffIsautodiffEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(QmuIsdakotaEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(InversionIscontrolEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(InversionTypeEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum));
+	if(solution_type==SeaiceSolutionEnum){
+	}
+	else{
+		parameters->AddObject(iomodel->CopyConstantObject(MeshElementtypeEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(SteadystateReltolEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(SteadystateMaxiterEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(ConstantsReferencetemperatureEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(GroundinglineMigrationEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(TransientIsstressbalanceEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(TransientIsmasstransportEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(TransientIsthermalEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(TransientIsgroundinglineEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(TransientIsgiaEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(TransientIslevelsetEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(TransientIsdamageevolutionEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(TransientIshydrologyEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum));
 
-	/*For stress balance only*/
-	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum));
-	if(iomodel->domaintype==Domain3DEnum)
-	 parameters->AddObject(iomodel->CopyConstantObject(MeshNumberoflayersEnum));
+		/*For stress balance only*/
+		parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum));
+		if(iomodel->domaintype==Domain3DEnum)
+		 parameters->AddObject(iomodel->CopyConstantObject(MeshNumberoflayersEnum));
+	}
 
 	/*Surface mass balance parameters*/
 	parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsEnum));
@@ -128,6 +132,12 @@
 		case SMBmeltcomponentsEnum:
 			/*Nothing to add to parameters*/
 			break;
+		case SeaiceatmEnum:
+			parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsRhoAirEnum));
+			parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsAirCoefEnum));
+			parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsAirLinDragCoefEnum));
+			parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsAirQuadDragCoefEnum));
+			break;
 		default:
 			_error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
 	}
@@ -144,6 +154,13 @@
 			parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsDeepwaterElevationEnum));
 			parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsUpperwaterElevationEnum));
 			break;
+		case SeaiceoceanEnum:
+			parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsRhoOceanEnum));
+			parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsOceanCoefEnum));
+			parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsOceanLinDragCoefEnum));
+			parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsOceanQuadDragCoefEnum));
+			parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsOceanTurningAngleEnum));
+			break;
 		default:
 			_error_("Basal forcing model "<<EnumToStringx(smb_model)<<" not supported yet");
 	}
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 18491)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 18492)
@@ -78,6 +78,10 @@
 					_error_("Mesh not supported yet");
 			}
 			break;
+		case MatseaiceEnum:
+			iomodel->FetchDataToInput(elements,MaterialsDamageEnum);
+			for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matseaice(i+1,i,iomodel));
+			break;
 		default:
 			_error_("Materials "<<EnumToStringx(materials_type)<<" not supported");
 	}
Index: ../trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp	(revision 18491)
+++ ../trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp	(revision 18492)
@@ -199,13 +199,16 @@
 		BTh.WriteMetric(bamgopts);
 
 		/*clean up*/
+		printf("-------------- file: Bamgx.cpp line: %i\n",__LINE__); 
 		delete &Th;
+		printf("-------------- file: Bamgx.cpp line: %i\n",__LINE__); 
 		//delete &BTh;
 		/*}}}*/
 	}
 
 	/*No error return*/
 	if (verbosity>1) _printf_("   Exiting Bamg.\n");
+	printf("-------------- file: Bamgx.cpp line: %i\n",__LINE__); 
 	return noerr;
 
 }
Index: ../trunk-jpl/src/c/Makefile.am
===================================================================
--- ../trunk-jpl/src/c/Makefile.am	(revision 18491)
+++ ../trunk-jpl/src/c/Makefile.am	(revision 18492)
@@ -107,6 +107,8 @@
 					./classes/Materials/Materials.cpp\
 					./classes/Materials/Matice.h\
 					./classes/Materials/Matice.cpp\
+					./classes/Materials/Matseaice.h\
+					./classes/Materials/Matseaice.cpp\
 					./classes/Materials/Matpar.h\
 					./classes/Materials/Matpar.cpp\
 					./classes/Constraints/Constraints.h\
@@ -648,6 +650,10 @@
 if DEPTHAVERAGE
 issm_sources += ./analyses/DepthAverageAnalysis.cpp
 endif
+if SEAICE
+issm_sources += ./cores/seaice_core.cpp
+issm_sources += ./analyses/SeaiceAnalysis.cpp
+endif
 if THERMAL
 issm_sources += ./analyses/ThermalAnalysis.cpp
 endif
Index: ../trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp	(revision 18491)
+++ ../trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp	(revision 18492)
@@ -68,6 +68,9 @@
 		case DamageEvolutionSolutionEnum:
 			solutioncore=&damage_core;
 			break;
+		case SeaiceSolutionEnum:
+			solutioncore=&seaice_core;
+			break;
 		default:
 			_error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!");
 			break;
Index: ../trunk-jpl/src/c/cores/seaice_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/seaice_core.cpp	(revision 0)
+++ ../trunk-jpl/src/c/cores/seaice_core.cpp	(revision 18492)
@@ -0,0 +1,30 @@
+/*!\file: seaice_core.cpp
+ * \brief: core of the seaice solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void seaice_core(FemModel* femmodel){
+
+	/*parameters: */
+	bool save_results;
+
+	/*recover parameters: */
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+
+	if(VerboseSolution()) _printf0_("call computational core:\n");
+	femmodel->SetCurrentConfiguration(SeaiceAnalysisEnum);
+	solutionsequence_linear(femmodel);
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		int outputs[6] = {VxEnum,VyEnum,VelEnum,SeaiceConcentrationEnum,SeaiceThicknessEnum,MaterialsDamageEnum};
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],6);
+	}
+
+}
Index: ../trunk-jpl/src/c/cores/AnalysisConfiguration.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/AnalysisConfiguration.cpp	(revision 18491)
+++ ../trunk-jpl/src/c/cores/AnalysisConfiguration.cpp	(revision 18492)
@@ -147,6 +147,12 @@
 			analyses[19]=L2ProjectionEPLAnalysisEnum;
 			break;
 
+		case SeaiceSolutionEnum:
+			numanalyses=1;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=SeaiceAnalysisEnum;
+			break;
+
 		default:
 			_error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!");
 			break;
Index: ../trunk-jpl/src/c/cores/cores.h
===================================================================
--- ../trunk-jpl/src/c/cores/cores.h	(revision 18491)
+++ ../trunk-jpl/src/c/cores/cores.h	(revision 18492)
@@ -44,6 +44,7 @@
 void dummy_core(FemModel* femmodel);
 void gia_core(FemModel* femmodel);
 void damage_core(FemModel* femmodel);
+void seaice_core(FemModel* femmodel);
 IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel);
 
 //optimization
Index: ../trunk-jpl/src/c/classes/Materials/Matseaice.h
===================================================================
--- ../trunk-jpl/src/c/classes/Materials/Matseaice.h	(revision 0)
+++ ../trunk-jpl/src/c/classes/Materials/Matseaice.h	(revision 18492)
@@ -0,0 +1,75 @@
+/*!\file Matseaice.h
+ * \brief: header file for matice object
+ */
+
+#ifndef MATSEAICE_H_
+#define MATSEAICE_H_
+
+/*Headers:*/
+/*{{{*/
+#include "./Material.h"
+#include "../Hook.h"
+class IoModel;
+class Elements;
+class Element;
+class Loads;
+class Nodes;
+class Vertices;
+class Materials;
+class Parameters;
+/*}}}*/
+
+class Matseaice: public Material{
+
+	private: 
+		int      mid;
+		Hook    *helement;
+		Element *element;
+
+	public:
+		/*Matseaice constructors, destructors: {{{*/
+		Matseaice();
+		Matseaice(int mid,int i, IoModel* iomodel);
+		~Matseaice();
+		/*}}}*/
+		/*Object virtual functions definitions:{{{ */
+		void  Echo();
+		void  DeepEcho();
+		int   Id(); 
+		int   ObjectEnum();
+		Object* copy();
+		/*}}}*/
+		/*Update virtual funictions definitions: {{{*/
+		void  InputUpdateFromVector(IssmDouble* vector, int name, int type){_error_("not implemented");};
+		void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrow, int ncols, int name, int type){_error_("not implemented");};
+		void  InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("not implemented");};
+		void  InputUpdateFromConstant(IssmDouble constant, int name){_error_("not implemented");};
+		void  InputUpdateFromConstant(int constant, int name){_error_("not implemented");};
+		void  InputUpdateFromConstant(bool constant, int name){_error_("not implemented");};
+		void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented");};
+		/*}}}*/
+		/*Material virtual functions resolution: {{{*/
+		void   Configure(Elements* elements);
+		Material*  copy2(Element* element);
+		void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
+		void       GetViscosity(IssmDouble* pviscosity, IssmDouble eps_eff){_error_("not supported for this material");};
+		void       GetViscosity_B(IssmDouble* pviscosity, IssmDouble eps_eff){_error_("not supported for this material");};
+		void       GetViscosity_D(IssmDouble* pviscosity, IssmDouble eps_eff){_error_("not supported for this material");};
+		void       GetViscosityBar(IssmDouble* pviscosity, IssmDouble eps_eff){_error_("not supported for this material");};
+		void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported for this material");};
+		void       GetViscosityDComplement(IssmDouble*, IssmDouble*){_error_("not supported for this material");};
+		void       GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported for this material");};
+		void       GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported for this material");};
+		IssmDouble GetA(){_error_("not supported for this material");};
+		IssmDouble GetAbar(){_error_("not supported for this material");};
+		IssmDouble GetB(){_error_("not supported for this material");};
+		IssmDouble GetBbar(){_error_("not supported for this material");};
+		IssmDouble GetD(){_error_("not supported for this material");};
+		IssmDouble GetDbar(){_error_("not supported for this material");};
+		IssmDouble GetN(){_error_("not supported for this material");};
+		bool       IsDamage(){_error_("not supported for this material");};
+		void       ResetHooks();
+		/*}}}*/
+};
+
+#endif  /* _MATICE_H_ */
Index: ../trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 18491)
+++ ../trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 18492)
@@ -19,76 +19,94 @@
 Matpar::Matpar(int matpar_mid, IoModel* iomodel){/*{{{*/
 
 	bool isefficientlayer;
-	int  hydrology_model,smb_model;
+	int  hydrology_model,smb_model,materials_type;
 	iomodel->Constant(&hydrology_model,HydrologyModelEnum);
 	iomodel->Constant(&smb_model,SurfaceforcingsEnum);
+	iomodel->Constant(&materials_type,MaterialsEnum);
 
 	this->mid = matpar_mid;
-	iomodel->Constant(&this->rho_ice,MaterialsRhoIceEnum);
-	iomodel->Constant(&this->rho_water,MaterialsRhoSeawaterEnum);
-	iomodel->Constant(&this->rho_freshwater,MaterialsRhoFreshwaterEnum);
-	iomodel->Constant(&this->mu_water,MaterialsMuWaterEnum);
-	iomodel->Constant(&this->heatcapacity,MaterialsHeatcapacityEnum);
-	iomodel->Constant(&this->thermalconductivity,MaterialsThermalconductivityEnum);
-	iomodel->Constant(&this->temperateiceconductivity,MaterialsTemperateiceconductivityEnum);
-	iomodel->Constant(&this->latentheat,MaterialsLatentheatEnum);
-	iomodel->Constant(&this->beta,MaterialsBetaEnum);
-	iomodel->Constant(&this->meltingpoint,MaterialsMeltingpointEnum);
-	iomodel->Constant(&this->referencetemperature,ConstantsReferencetemperatureEnum);
-	iomodel->Constant(&this->mixed_layer_capacity,MaterialsMixedLayerCapacityEnum);
-	iomodel->Constant(&this->thermal_exchange_velocity,MaterialsThermalExchangeVelocityEnum);
-	iomodel->Constant(&this->g,ConstantsGEnum);
 
-	switch(smb_model){
-		case SMBEnum:
-			/*Nothing to add*/
+	switch(materials_type){
+		case MaticeEnum:
+			iomodel->Constant(&this->rho_ice,MaterialsRhoIceEnum);
+			iomodel->Constant(&this->rho_water,MaterialsRhoSeawaterEnum);
+			iomodel->Constant(&this->rho_freshwater,MaterialsRhoFreshwaterEnum);
+			iomodel->Constant(&this->mu_water,MaterialsMuWaterEnum);
+			iomodel->Constant(&this->heatcapacity,MaterialsHeatcapacityEnum);
+			iomodel->Constant(&this->thermalconductivity,MaterialsThermalconductivityEnum);
+			iomodel->Constant(&this->temperateiceconductivity,MaterialsTemperateiceconductivityEnum);
+			iomodel->Constant(&this->latentheat,MaterialsLatentheatEnum);
+			iomodel->Constant(&this->beta,MaterialsBetaEnum);
+			iomodel->Constant(&this->meltingpoint,MaterialsMeltingpointEnum);
+			iomodel->Constant(&this->referencetemperature,ConstantsReferencetemperatureEnum);
+			iomodel->Constant(&this->mixed_layer_capacity,MaterialsMixedLayerCapacityEnum);
+			iomodel->Constant(&this->thermal_exchange_velocity,MaterialsThermalExchangeVelocityEnum);
+			iomodel->Constant(&this->g,ConstantsGEnum);
+
+			switch(smb_model){
+				case SMBEnum:
+					/*Nothing to add*/
+					break;
+				case SMBpddEnum:
+					iomodel->Constant(&this->desfac,SurfaceforcingsDesfacEnum);
+					iomodel->Constant(&this->s0p,SurfaceforcingsS0pEnum);
+					break;
+				case SMBgradientsEnum:
+					/*Nothing to add*/
+					break;
+				case SMBhenningEnum:
+					/*Nothing to add*/
+					break;
+				case SMBcomponentsEnum:
+					/*Nothing to add*/
+					break;
+				case SMBmeltcomponentsEnum:
+					/*Nothing to add*/
+					break;
+				default:
+					_error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
+			}
+			if(hydrology_model==HydrologydcEnum){
+				iomodel->Constant(&this->sediment_compressibility,HydrologydcSedimentCompressibilityEnum);
+				iomodel->Constant(&this->sediment_porosity,HydrologydcSedimentPorosityEnum);
+				iomodel->Constant(&this->sediment_thickness,HydrologydcSedimentThicknessEnum);
+				iomodel->Constant(&this->water_compressibility,HydrologydcWaterCompressibilityEnum);
+				iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
+
+				if(isefficientlayer){
+					iomodel->Constant(&this->epl_compressibility,HydrologydcEplCompressibilityEnum);
+					iomodel->Constant(&this->epl_porosity,HydrologydcEplPorosityEnum);
+					iomodel->Constant(&this->epl_init_thickness,HydrologydcEplInitialThicknessEnum);
+					iomodel->Constant(&this->epl_max_thickness,HydrologydcEplMaxThicknessEnum);
+					iomodel->Constant(&this->epl_conductivity,HydrologydcEplConductivityEnum);
+				}
+			}
+			else if(hydrology_model==HydrologyshreveEnum){
+				/*Nothing to add*/
+			}
+			else{
+				_error_("Hydrology model "<<EnumToStringx(hydrology_model)<<" not supported yet");
+			}
+
+			/*gia: */
+			iomodel->Constant(&this->lithosphere_shear_modulus,MaterialsLithosphereShearModulusEnum);
+			iomodel->Constant(&this->lithosphere_density,MaterialsLithosphereDensityEnum);
+			iomodel->Constant(&this->mantle_shear_modulus,MaterialsMantleShearModulusEnum);
+			iomodel->Constant(&this->mantle_density,MaterialsMantleDensityEnum);
 			break;
-		case SMBpddEnum:
-			iomodel->Constant(&this->desfac,SurfaceforcingsDesfacEnum);
-			iomodel->Constant(&this->s0p,SurfaceforcingsS0pEnum);
+		case MatseaiceEnum:
+			iomodel->Constant(&this->poisson,MaterialsPoissonEnum);
+			iomodel->Constant(&this->young_modulus,MaterialsYoungModulusEnum);
+			iomodel->Constant(&this->ridging_exponent,MaterialsRidgingExponentEnum);
+			iomodel->Constant(&this->cohesion,MaterialsCohesionEnum);
+			iomodel->Constant(&this->internal_friction_coef,MaterialsInternalFrictionCoefEnum);
+			iomodel->Constant(&this->compression_coef,MaterialsCompressionCoefEnum);
+			iomodel->Constant(&this->traction_coef,MaterialsTractionCoefEnum);
+			iomodel->Constant(&this->g,ConstantsGEnum);
 			break;
-		case SMBgradientsEnum:
-			/*Nothing to add*/
-			break;
-		case SMBhenningEnum:
-			/*Nothing to add*/
-			break;
-		case SMBcomponentsEnum:
-			/*Nothing to add*/
-			break;
-		case SMBmeltcomponentsEnum:
-			/*Nothing to add*/
-			break;
 		default:
-			_error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
+			_error_("Material not supported yet");
 	}
-	if(hydrology_model==HydrologydcEnum){
-		iomodel->Constant(&this->sediment_compressibility,HydrologydcSedimentCompressibilityEnum);
-		iomodel->Constant(&this->sediment_porosity,HydrologydcSedimentPorosityEnum);
-		iomodel->Constant(&this->sediment_thickness,HydrologydcSedimentThicknessEnum);
-		iomodel->Constant(&this->water_compressibility,HydrologydcWaterCompressibilityEnum);
-		iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
-
-		if(isefficientlayer){
-				iomodel->Constant(&this->epl_compressibility,HydrologydcEplCompressibilityEnum);
-				iomodel->Constant(&this->epl_porosity,HydrologydcEplPorosityEnum);
-				iomodel->Constant(&this->epl_init_thickness,HydrologydcEplInitialThicknessEnum);
-				iomodel->Constant(&this->epl_max_thickness,HydrologydcEplMaxThicknessEnum);
-				iomodel->Constant(&this->epl_conductivity,HydrologydcEplConductivityEnum);
-		}
-	}
-	else if(hydrology_model==HydrologyshreveEnum){
-		/*Nothing to add*/
-	}
-	else{
-		_error_("Hydrology model "<<EnumToStringx(hydrology_model)<<" not supported yet");
-	}
-
-	/*gia: */
-	iomodel->Constant(&this->lithosphere_shear_modulus,MaterialsLithosphereShearModulusEnum);
-	iomodel->Constant(&this->lithosphere_density,MaterialsLithosphereDensityEnum);
-	iomodel->Constant(&this->mantle_shear_modulus,MaterialsMantleShearModulusEnum);
-	iomodel->Constant(&this->mantle_density,MaterialsMantleDensityEnum);
 }
 /*}}}*/
 Matpar::~Matpar(){/*{{{*/
@@ -292,6 +310,13 @@
 		case HydrologydcEplMaxThicknessEnum:         return this->epl_max_thickness;
 		case HydrologydcWaterCompressibilityEnum:    return this->water_compressibility;
 		case ConstantsGEnum:                         return this->g;
+		case MaterialsPoissonEnum:                   return this->poisson;
+		case MaterialsYoungModulusEnum:              return this->young_modulus;
+		case MaterialsRidgingExponentEnum:				return this->ridging_exponent;
+		case MaterialsCohesionEnum:		            return this->cohesion;
+		case MaterialsInternalFrictionCoefEnum:		return this->internal_friction_coef;
+		case MaterialsCompressionCoefEnum:				return this->compression_coef;
+		case MaterialsTractionCoefEnum:					return this->traction_coef;
 		default: _error_("Enum "<<EnumToStringx(enum_in)<<" not supported yet");
 	}
 
Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h
===================================================================
--- ../trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 18491)
+++ ../trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 18492)
@@ -50,6 +50,15 @@
 		IssmDouble mantle_shear_modulus;
 		IssmDouble mantle_density;
 
+		/*Sea ice*/
+		IssmDouble poisson;
+		IssmDouble young_modulus;
+		IssmDouble ridging_exponent;
+		IssmDouble cohesion;
+		IssmDouble internal_friction_coef;
+		IssmDouble compression_coef;
+		IssmDouble traction_coef;
+
 	public:
 		Matpar();
 		Matpar(int matpar_id, IoModel* iomodel);
Index: ../trunk-jpl/src/c/classes/Materials/Matseaice.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Materials/Matseaice.cpp	(revision 0)
+++ ../trunk-jpl/src/c/classes/Materials/Matseaice.cpp	(revision 18492)
@@ -0,0 +1,136 @@
+/*!\file Matseaice.c
+ * \brief: implementation of the Matseaice object
+ */
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./Matseaice.h"
+#include "./Materials.h"
+#include "../Inputs/Input.h"
+#include "../Inputs/Inputs.h"
+#include "../Inputs/TriaInput.h"
+#include "../Inputs/PentaInput.h"
+#include "../Inputs/ControlInput.h"
+#include "../Elements/Element.h"
+#include "../Elements/Tria.h"
+#include "../Elements/Penta.h"
+#include "../Params/Parameters.h"
+#include "../Vertex.h"
+#include "../Hook.h"
+#include "../Node.h"
+#include "../IoModel.h"
+#include "../../shared/shared.h"
+
+/*Matseaice constructors and destructor*/
+Matseaice::Matseaice(){/*{{{*/
+	this->helement=NULL;
+	this->element=NULL;
+	return;
+}
+/*}}}*/
+Matseaice::Matseaice(int matice_mid,int index, IoModel* iomodel){/*{{{*/
+
+	/*Intermediaries:*/
+	int    matice_eid;
+
+	/*Initialize id*/
+	this->mid=matice_mid;
+
+	/*Hooks: */
+	matice_eid=index+1;
+	this->helement=new Hook(&matice_eid,1);
+	this->element=NULL;
+
+	return;
+}
+/*}}}*/
+Matseaice::~Matseaice(){/*{{{*/
+	delete helement;
+	return;
+}
+/*}}}*/
+
+/*Object virtual functions definitions:*/
+void Matseaice::Echo(void){/*{{{*/
+
+	_printf_("Matseaice:\n");
+	_printf_("   mid: " << mid << "\n");
+	_printf_("   element:\n");
+	helement->Echo();
+}
+/*}}}*/
+void Matseaice::DeepEcho(void){/*{{{*/
+
+	_printf_("Matseaice:\n");
+	_printf_("   mid: " << mid << "\n");
+	_printf_("   element:\n");
+	helement->Echo();
+}		
+/*}}}*/
+int  Matseaice::Id(void){ return mid; }/*{{{*/
+/*}}}*/
+int Matseaice::ObjectEnum(void){/*{{{*/
+
+	return MatseaiceEnum;
+
+}
+/*}}}*/
+Object* Matseaice::copy() {/*{{{*/
+
+	/*Output*/
+	Matseaice* matice=NULL;
+
+	/*Initialize output*/
+	matice=new Matseaice();
+
+	/*copy fields: */
+	matice->mid=this->mid;
+	matice->helement=(Hook*)this->helement->copy();
+	matice->element =(Element*)this->helement->delivers();
+
+	return matice;
+}
+/*}}}*/
+Material* Matseaice::copy2(Element* element_in) {/*{{{*/
+
+	/*Output*/
+	Matseaice* matice=NULL;
+
+	/*Initialize output*/
+	matice=new Matseaice();
+
+	/*copy fields: */
+	matice->mid=this->mid;
+	matice->helement=(Hook*)this->helement->copy();
+	matice->element =element_in;
+
+	return matice;
+}
+/*}}}*/
+
+/*Matseaice management*/
+void  Matseaice::Configure(Elements* elementsin){/*{{{*/
+
+	/*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective 
+	 * datasets, using internal ids and offsets hidden in hooks: */
+	helement->configure((DataSet*)elementsin);
+	this->element  = (Element*)helement->delivers();
+}
+/*}}}*/
+void  Matseaice::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
+
+}
+/*}}}*/
+void  Matseaice::ResetHooks(){/*{{{*/
+
+	this->element=NULL;
+
+	/*Get Element type*/
+	this->helement->reset();
+
+}
+/*}}}*/
Index: ../trunk-jpl/src/c/classes/classes.h
===================================================================
--- ../trunk-jpl/src/c/classes/classes.h	(revision 18491)
+++ ../trunk-jpl/src/c/classes/classes.h	(revision 18492)
@@ -77,6 +77,7 @@
 #include "./Materials/Materials.h"
 #include "./Materials/Material.h"
 #include "./Materials/Matice.h"
+#include "./Materials/Matseaice.h"
 #include "./Materials/Matpar.h"
 
 /*Params: */
Index: ../trunk-jpl/src/c/classes/Vertex.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Vertex.cpp	(revision 18491)
+++ ../trunk-jpl/src/c/classes/Vertex.cpp	(revision 18492)
@@ -31,17 +31,21 @@
 	this->z            = iomodel->Data(MeshZEnum)[i];
 	this->domaintype     = iomodel->domaintype;
 
-	_assert_(iomodel->Data(BaseEnum) && iomodel->Data(ThicknessEnum) && iomodel->numbernodetoelementconnectivity);
 	switch(iomodel->domaintype){
 		case Domain3DEnum:
-		case Domain2DhorizontalEnum:
+			_assert_(iomodel->Data(BaseEnum) && iomodel->Data(ThicknessEnum));
 			this->sigma = (iomodel->Data(MeshZEnum)[i]-iomodel->Data(BaseEnum)[i])/(iomodel->Data(ThicknessEnum)[i]);
 			break;
+		case Domain2DhorizontalEnum:
+			this->sigma = 0.;
+			break;
 		case Domain2DverticalEnum:
+			_assert_(iomodel->Data(BaseEnum) && iomodel->Data(ThicknessEnum));
 			this->sigma = (iomodel->Data(MeshYEnum)[i]-iomodel->Data(BaseEnum)[i])/(iomodel->Data(ThicknessEnum)[i]);
 			break;
 	}
 
+	_assert_(iomodel->numbernodetoelementconnectivity);
 	this->connectivity = iomodel->numbernodetoelementconnectivity[i];
 
 }
Index: ../trunk-jpl/src/m/consistency/ismodelselfconsistent.m
===================================================================
--- ../trunk-jpl/src/m/consistency/ismodelselfconsistent.m	(revision 18491)
+++ ../trunk-jpl/src/m/consistency/ismodelselfconsistent.m	(revision 18492)
@@ -90,6 +90,9 @@
 	case DamageEvolutionSolutionEnum(),
 		analyses=[DamageEvolutionAnalysisEnum()];
 
+	case SeaiceSolutionEnum(),
+		analyses=[SeaiceAnalysisEnum()];
+
 	otherwise
 		error('%s%s%s',' solution type: ',EnumToString(solutiontype),' not supported yet!');
 
Index: ../trunk-jpl/src/m/classes/mask.m
===================================================================
--- ../trunk-jpl/src/m/classes/mask.m	(revision 18491)
+++ ../trunk-jpl/src/m/classes/mask.m	(revision 18492)
@@ -39,6 +39,9 @@
 
 		end % }}}
 		function md = checkconsistency(obj,md,solution,analyses) % {{{
+			if (solution==SeaiceSolutionEnum()),
+				return;
+			end
 
 			md = checkfield(md,'fieldname','mask.groundedice_levelset','size',[md.mesh.numberofvertices 1]);
 			md = checkfield(md,'fieldname','mask.ice_levelset'        ,'size',[md.mesh.numberofvertices 1]);
@@ -60,6 +63,10 @@
 			WriteData(fid,'object',obj,'fieldname','groundedice_levelset','format','DoubleMat','mattype',1);
 			WriteData(fid,'object',obj,'fieldname','ice_levelset','format','DoubleMat','mattype',1);
 
+			if(md.private.solution==SeaiceSolutionEnum()),
+				return;
+			end
+
 			% get mask of vertices of elements with ice
 			isice=md.mask.ice_levelset<=0.;
 			vlist = zeros(md.mesh.numberofvertices,1);
Index: ../trunk-jpl/src/m/classes/seaice.m
===================================================================
--- ../trunk-jpl/src/m/classes/seaice.m	(revision 0)
+++ ../trunk-jpl/src/m/classes/seaice.m	(revision 18492)
@@ -0,0 +1,41 @@
+%SEAICE class definition
+%
+%   Usage:
+%      seaice=seaice();
+
+classdef seaice
+	properties (SetAccess=public) 
+		thickness     = NaN;
+		concentration = NaN;
+		spcvx         = NaN;
+		spcvy         = NaN;
+	end
+	methods
+		function obj = seaice(varargin) % {{{
+			switch nargin
+				case 0
+					return;
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function md = checkconsistency(obj,md,solution,analyses) % {{{
+			md = checkfield(md,'fieldname','seaice.thickness','size',[md.mesh.numberofelements 1],'NaN',1,'>=',0);
+			md = checkfield(md,'fieldname','seaice.concentration','size',[md.mesh.numberofelements 1],'NaN',1,'>=',0);
+			md = checkfield(md,'fieldname','seaice.spcvx','size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'fieldname','seaice.spcvy','size',[md.mesh.numberofvertices 1]);
+		end % }}}
+		function disp(obj) % {{{
+			disp(sprintf('   seaice parameters:'));
+			fielddisplay(obj,'thickness','sea ice thickness [m]');
+			fielddisplay(obj,'spcvx','x-axis velocity constraint (NaN means no constraint) [m/s]');
+			fielddisplay(obj,'spcvy','y-axis velocity constraint (NaN means no constraint) [m/s]');
+		end % }}}
+		function marshall(obj,md,fid) % {{{
+			WriteData(fid,'object',obj,'fieldname','thickness','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',obj,'fieldname','concentration','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',obj,'fieldname','spcvx','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','spcvy','format','DoubleMat','mattype',1);
+		end % }}}
+	end
+end
Index: ../trunk-jpl/src/m/classes/geometry.m
===================================================================
--- ../trunk-jpl/src/m/classes/geometry.m	(revision 18491)
+++ ../trunk-jpl/src/m/classes/geometry.m	(revision 18492)
@@ -57,6 +57,10 @@
 		end % }}}
 		function md = checkconsistency(obj,md,solution,analyses) % {{{
 
+			if (solution==SeaiceSolutionEnum()),
+				return;
+			end
+
 			if (solution==TransientSolutionEnum() & md.transient.isgia) | (solution==GiaSolutionEnum()),
 				md = checkfield(md,'fieldname','geometry.thickness','forcing',1,'NaN',1,'>=',0);
 			else
Index: ../trunk-jpl/src/m/classes/matseaice.m
===================================================================
--- ../trunk-jpl/src/m/classes/matseaice.m	(revision 0)
+++ ../trunk-jpl/src/m/classes/matseaice.m	(revision 18492)
@@ -0,0 +1,94 @@
+%MATSEAICE class definition
+%
+%   Usage:
+%      matseaice=matseaice();
+
+classdef matseaice
+	properties (SetAccess=public) 
+		poisson                = 0.;
+		young_modulus          = 0.;
+		damage                 = NaN;
+		ridging_exponent       = 0.;
+		cohesion               = 0.;
+		internal_friction_coef = 0.;
+		compression_coef       = 0.;
+		traction_coef          = 0.;
+	end
+	methods
+		function obj = matseaice(varargin) % {{{
+			switch nargin
+				case 0
+					obj=setdefaultparameters(obj);
+				case 1
+					inputstruct=varargin{1};
+					list1 = properties('matseaice');
+					list2 = fieldnames(inputstruct);
+					for i=1:length(list1)
+						fieldname = list1{i};
+						if ismember(fieldname,list2),
+							obj.(fieldname) = inputstruct.(fieldname);
+						end
+					end
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function obj = setdefaultparameters(obj) % {{{
+
+			%poisson modulus
+			obj.poisson=0.3;
+
+			%Young modulus
+			obj.young_modulus=9.e+9;
+
+			%ridging exponent
+			obj.ridging_exponent=-20.;
+
+			%Cohesion (Weiss 2007)
+			obj.cohesion = 40.e+3; %40 kPa
+
+			%internal_friction_coef 
+			obj.internal_friction_coef=0.7;
+
+			%compression_coef (Bouillon and Rampal 2014)
+			obj.compression_coef=5./2;
+
+			%traction_coef (Bouillon and Rampal 2014)
+			obj.traction_coef=5./6.;
+
+		end % }}}
+		function md = checkconsistency(obj,md,solution,analyses) % {{{
+			md = checkfield(md,'fieldname','materials.poisson','NaN',1,'>',0,'<',.5,'numel',1);
+			md = checkfield(md,'fieldname','materials.young_modulus','NaN',1,'>',0,'numel',1);
+			md = checkfield(md,'fieldname','materials.damage','NaN',1,'>=',0,'<',1,'size',[md.mesh.numberofelements 1]);
+			md = checkfield(md,'fieldname','materials.ridging_exponent','NaN',1,'<',0,'numel',1);
+			md = checkfield(md,'fieldname','materials.cohesion','NaN',1,'>',0,'numel',1);
+			md = checkfield(md,'fieldname','materials.internal_friction_coef','NaN',1,'>',0,'numel',1);
+			md = checkfield(md,'fieldname','materials.compression_coef','NaN',1,'>',1,'numel',1);
+			md = checkfield(md,'fieldname','materials.traction_coef','NaN',1,'>',0,'<',1,'numel',1);
+		end % }}}
+		function disp(obj) % {{{
+			disp(sprintf('   Sea Ice Material:'));
+
+			fielddisplay(obj,'poisson','poisson ratio for undamaged ice [no unit]');
+			fielddisplay(obj,'young_modulus','Young modulus for undamaged ice [Pa]');
+			fielddisplay(obj,'damage','damage, between 0 (no damage) and 1 (fully damaged)');
+			fielddisplay(obj,'ridging_exponent','Riging exponent (c, Hibler parameter) [no unit]');
+			fielddisplay(obj,'cohesion','cohesion (C) [Pa]');
+			fielddisplay(obj,'internal_friction_coef','Mohr-Coulomb internal friction coefficient (mu=tan(phi)) [no unit]');
+			fielddisplay(obj,'compression_coef','Ratio between cutoff compressive strength and the cohesion [no unit]');
+			fielddisplay(obj,'traction_coef','Ratio between cutoff tensile strength and Mohr-Coulomb tensile strength [no unit]');
+		end % }}}
+		function marshall(obj,md,fid) % {{{
+			WriteData(fid,'enum',MaterialsEnum(),'data',MatseaiceEnum(),'format','Integer');
+			WriteData(fid,'object',obj,'class','materials','fieldname','poisson','format','Double');
+			WriteData(fid,'object',obj,'class','materials','fieldname','young_modulus','format','Double');
+			WriteData(fid,'object',obj,'class','materials','fieldname','damage','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',obj,'class','materials','fieldname','ridging_exponent','format','Double');
+			WriteData(fid,'object',obj,'class','materials','fieldname','cohesion','format','Double');
+			WriteData(fid,'object',obj,'class','materials','fieldname','internal_friction_coef','format','Double');
+			WriteData(fid,'object',obj,'class','materials','fieldname','compression_coef','format','Double');
+			WriteData(fid,'object',obj,'class','materials','fieldname','traction_coef','format','Double');
+		end % }}}
+	end
+end
Index: ../trunk-jpl/src/m/classes/seaiceinitialization.m
===================================================================
--- ../trunk-jpl/src/m/classes/seaiceinitialization.m	(revision 0)
+++ ../trunk-jpl/src/m/classes/seaiceinitialization.m	(revision 18492)
@@ -0,0 +1,42 @@
+%SEAICEINITIALIZATION class definition
+%
+%   Usage:
+%      seaiceinitialization=seaiceinitialization();
+
+classdef seaiceinitialization
+	properties (SetAccess=public) 
+		vx          = NaN;
+		vy          = NaN;
+		vx_coriolis = NaN;
+		vy_coriolis = NaN;
+	end
+	methods
+		function obj = seaiceinitialization(varargin) % {{{
+			switch nargin
+				case 0
+					return;
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function md = checkconsistency(obj,md,solution,analyses) % {{{
+			md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'fieldname','initialization.vx_coriolis','NaN',1,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'fieldname','initialization.vy_coriolis','NaN',1,'size',[md.mesh.numberofvertices 1]);
+		end % }}}
+		function disp(obj) % {{{
+			disp(sprintf('   initial field values:'));
+			fielddisplay(obj,'vx','x component of the ice velocity [m/s]');
+			fielddisplay(obj,'vy','y component of the ice velocity [m/s]');
+			fielddisplay(obj,'vx_coriolis','x component of the ice velocity used to calculate coriolis forces [m/s]');
+			fielddisplay(obj,'vy_coriolis','y component of the ice velocity used to calculate coriolis forces [m/s]');
+		end % }}}
+		function marshall(obj,md,fid) % {{{
+			WriteData(fid,'data',obj.vx,'format','DoubleMat','mattype',1,'enum',VxEnum);
+			WriteData(fid,'data',obj.vy,'format','DoubleMat','mattype',1,'enum',VyEnum);
+			WriteData(fid,'data',obj.vx_coriolis,'format','DoubleMat','mattype',1,'enum',VxStarEnum);
+			WriteData(fid,'data',obj.vy_coriolis,'format','DoubleMat','mattype',1,'enum',VyStarEnum);
+		end % }}}
+	end
+end
Index: ../trunk-jpl/src/m/classes/stressbalance.m
===================================================================
--- ../trunk-jpl/src/m/classes/stressbalance.m	(revision 18491)
+++ ../trunk-jpl/src/m/classes/stressbalance.m	(revision 18492)
@@ -209,6 +209,12 @@
 		end % }}}
 		function marshall(obj,md,fid) % {{{
 
+			WriteData(fid,'object',obj,'class','stressbalance','fieldname','vertex_pairing','format','DoubleMat','mattype',3);
+
+			if md.private.solution==SeaiceSolutionEnum,
+				return;
+			end
+
 			yts=365.0*24.0*3600.0;
 
 			WriteData(fid,'object',obj,'class','stressbalance','fieldname','spcvx','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
@@ -222,7 +228,6 @@
 			WriteData(fid,'object',obj,'class','stressbalance','fieldname','viscosity_overshoot','format','Double');
 			WriteData(fid,'object',obj,'class','stressbalance','fieldname','maxiter','format','Integer');
 			WriteData(fid,'object',obj,'class','stressbalance','fieldname','shelf_dampening','format','Integer');
-			WriteData(fid,'object',obj,'class','stressbalance','fieldname','vertex_pairing','format','DoubleMat','mattype',3);
 			WriteData(fid,'object',obj,'class','stressbalance','fieldname','penalty_factor','format','Double');
 			WriteData(fid,'object',obj,'class','stressbalance','fieldname','rift_penalty_lock','format','Integer');
 			WriteData(fid,'object',obj,'class','stressbalance','fieldname','rift_penalty_threshold','format','Integer');
Index: ../trunk-jpl/src/m/classes/seaiceatm.m
===================================================================
--- ../trunk-jpl/src/m/classes/seaiceatm.m	(revision 0)
+++ ../trunk-jpl/src/m/classes/seaiceatm.m	(revision 18492)
@@ -0,0 +1,65 @@
+%SEAICEATM Class definition
+%
+%   Usage:
+%      seaiceatm=seaiceatm();
+
+classdef seaiceatm
+	properties (SetAccess=public) 
+		rho_air            = 0.;
+		air_coef           = 0.;
+		air_lin_drag_coef  = 0.;
+		air_quad_drag_coef = 0.;
+		wind_vx            = NaN;
+		wind_vy            = NaN;
+	end
+	methods
+		function obj = seaiceatm(varargin) % {{{
+			switch nargin
+				case 0
+					obj=setdefaultparameters(obj);
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function obj = setdefaultparameters(obj) % {{{
+
+			%air density (kg/m^3)
+			obj.rho_air=1.3;
+
+			%By default, we apply the full stress
+			obj.air_coef = 1.;
+
+			%drag coefficients (different from Hibler 0.0012)
+			obj.air_lin_drag_coef  = 0.;
+			obj.air_quad_drag_coef = 0.003;
+
+		end % }}}
+		function md = checkconsistency(obj,md,solution,analyses) % {{{
+			md = checkfield(md,'fieldname','surfaceforcings.rho_air','NaN',1,'>',0.,'numel',1);
+			md = checkfield(md,'fieldname','surfaceforcings.air_coef','NaN',1,'>=',0.,'<=',1,'numel',1);
+			md = checkfield(md,'fieldname','surfaceforcings.air_lin_drag_coef','NaN',1,'>=',0.,'numel',1);
+			md = checkfield(md,'fieldname','surfaceforcings.air_quad_drag_coef','NaN',1,'>=',0.,'numel',1);
+			md = checkfield(md,'fieldname','surfaceforcings.wind_vx','NaN',1,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'fieldname','surfaceforcings.wind_vy','NaN',1,'size',[md.mesh.numberofvertices 1]);
+		end % }}}
+		function disp(obj) % {{{
+			disp(sprintf('   Atmospheric forcings for sea ice:'));
+			disp(sprintf('\n     tau_a = air_coef*rho_air*(air_lin_drag_coef+air_quad_drag_coef*|v_wind|)*v_wind\n'));
+			fielddisplay(obj,'rho_air','air density [kg/m^3]');
+			fielddisplay(obj,'air_coef','atmospheric stress spin up coefficient (between 0 and 1)');
+			fielddisplay(obj,'air_lin_drag_coef','air linear drag coefficient (no unit)');
+			fielddisplay(obj,'air_quad_drag_coef','air quadratic drag coefficient [s/m]');
+			fielddisplay(obj,'wind_vx','wind speed x-component at 10 m elevation [m/s]');
+			fielddisplay(obj,'wind_vy','wind speed y-component at 10 m elevation [m/s]');
+		end % }}}
+		function marshall(obj,md,fid) % {{{
+			WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SeaiceatmEnum(),'format','Integer');
+			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','rho_air','format','Double');
+			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','air_coef','format','Double');
+			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','air_lin_drag_coef','format','Double');
+			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','air_quad_drag_coef','format','Double');
+			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','wind_vx','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','wind_vy','format','DoubleMat','mattype',1);
+		end % }}}
+	end
+end
Index: ../trunk-jpl/src/m/classes/seaiceocean.m
===================================================================
--- ../trunk-jpl/src/m/classes/seaiceocean.m	(revision 0)
+++ ../trunk-jpl/src/m/classes/seaiceocean.m	(revision 18492)
@@ -0,0 +1,76 @@
+%SEAICEOCEAN Class definition
+%
+%   Usage:
+%      seaiceocean=seaiceocean();
+
+classdef seaiceocean
+	properties (SetAccess=public) 
+		rho_ocean            = 0.;
+		ocean_coef           = 0.;
+		ocean_lin_drag_coef  = 0.;
+		ocean_quad_drag_coef = 0.;
+		ocean_turning_angle  = 0.;
+		ocean_ssh            = NaN;
+		ocean_vx             = NaN;
+		ocean_vy             = NaN;
+	end
+	methods
+		function obj = seaiceocean(varargin) % {{{
+			switch nargin
+				case 0
+					obj=setdefaultparameters(obj);
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function obj = setdefaultparameters(obj) % {{{
+
+			%ocean density (kg/m^3)
+			obj.rho_ocean=1025.;
+
+			%By default, we apply the full stress
+			obj.ocean_coef = 1.;
+
+			%drag coefficients (different from Hibler 0.0055)
+			obj.ocean_lin_drag_coef  = 0.;
+			obj.ocean_quad_drag_coef = 0.004;
+
+			%Turning angle in degrees (McPhee 1998)
+			obj.ocean_turning_angle = 25.;
+
+		end % }}}
+		function md = checkconsistency(obj,md,solution,analyses) % {{{
+			md = checkfield(md,'fieldname','basalforcings.rho_ocean','NaN',1,'>',0.,'numel',1);
+			md = checkfield(md,'fieldname','basalforcings.ocean_coef','NaN',1,'>=',0.,'<=',1.,'numel',1);
+			md = checkfield(md,'fieldname','basalforcings.ocean_lin_drag_coef','NaN',1,'>=',0.,'numel',1);
+			md = checkfield(md,'fieldname','basalforcings.ocean_quad_drag_coef','NaN',1,'>=',0.,'numel',1);
+			md = checkfield(md,'fieldname','basalforcings.ocean_turning_angle','NaN',1,'>=',0.,'<',90,'numel',1);
+			md = checkfield(md,'fieldname','basalforcings.ocean_ssh','NaN',1,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'fieldname','basalforcings.ocean_vx','NaN',1,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'fieldname','basalforcings.ocean_vy','NaN',1,'size',[md.mesh.numberofvertices 1]);
+		end % }}}
+		function disp(obj) % {{{
+			disp(sprintf('   Ocean forcings for sea ice:'));
+			%disp(sprintf('\n     tau_b = ocean_coef*rho_ocean*...'));
+			fielddisplay(obj,'rho_ocean','ocean density [kg/m^3]');
+			fielddisplay(obj,'ocean_coef','ocean stress spin up coefficient (between 0 and 1)');
+			fielddisplay(obj,'ocean_lin_drag_coef','ocean linear drag coefficient [Pa/(m/s)]');
+			fielddisplay(obj,'ocean_quad_drag_coef','ocean quadratic drag coefficient [Pa/(m/s)^2]');
+			fielddisplay(obj,'ocean_turning_angle','ocean turning angle [degree]');
+			fielddisplay(obj,'ocean_ssh','ocean sea surface height [m]');
+			fielddisplay(obj,'ocean_vx','ocean speed x-component [m/s]');
+			fielddisplay(obj,'ocean_vy','ocean speed y-component [m/s]');
+		end % }}}
+		function marshall(obj,md,fid) % {{{
+			WriteData(fid,'enum',BasalforcingsEnum(),'data',SeaiceoceanEnum(),'format','Integer');
+			WriteData(fid,'object',obj,'class','basalforcings','fieldname','rho_ocean','format','Double');
+			WriteData(fid,'object',obj,'class','basalforcings','fieldname','ocean_coef','format','Double');
+			WriteData(fid,'object',obj,'class','basalforcings','fieldname','ocean_lin_drag_coef','format','Double');
+			WriteData(fid,'object',obj,'class','basalforcings','fieldname','ocean_quad_drag_coef','format','Double');
+			WriteData(fid,'object',obj,'class','basalforcings','fieldname','ocean_turning_angle','format','Double');
+			WriteData(fid,'object',obj,'class','basalforcings','fieldname','ocean_ssh','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'class','basalforcings','fieldname','ocean_vx','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'class','basalforcings','fieldname','ocean_vy','format','DoubleMat','mattype',1);
+		end % }}}
+	end
+end
Index: ../trunk-jpl/src/m/classes/model.m
===================================================================
--- ../trunk-jpl/src/m/classes/model.m	(revision 18491)
+++ ../trunk-jpl/src/m/classes/model.m	(revision 18492)
@@ -38,6 +38,8 @@
 		transient        = 0;
 		gia              = 0;
 
+		seaice           = 0;
+
 		autodiff         = 0;
 		flaim            = 0;
 		inversion        = 0;
@@ -131,6 +133,17 @@
 			switch nargin
 				case 0
 					md=setdefaultparameters(md);
+				case 1
+					if strcmpi(varargin{1},'seaice'),
+						md=setdefaultparameters(md);
+						md.materials = matseaice();
+						md.surfaceforcings = seaiceatm();
+						md.basalforcings   = seaiceocean();
+						md.initialization  = seaiceinitialization();
+					else
+						error('model constructor not supported yet');
+					end
+
 				otherwise
 					error('model constructor error message: 0 of 1 argument only in input.');
 				end
@@ -1117,6 +1130,7 @@
 			md.steadystate      = steadystate();
 			md.transient        = transient();
 			md.gia              = gia();
+			md.seaice           = seaice();
 			md.autodiff         = autodiff();
 			md.flaim            = flaim();
 			md.inversion        = inversion();
@@ -1294,6 +1308,7 @@
 			disp(sprintf('%19s: %-22s -- %s','results'         ,['[1x1 ' class(obj.results) ']'],'model results'));
 			disp(sprintf('%19s: %-22s -- %s','radaroverlay'    ,['[1x1 ' class(obj.radaroverlay) ']'],'radar image for plot overlay'));
 			disp(sprintf('%19s: %-22s -- %s','miscellaneous'   ,['[1x1 ' class(obj.miscellaneous) ']'],'miscellaneous fields'));
+			disp(sprintf('%19s: %-22s -- %s','seaice'       ,['[1x1 ' class(obj.seaice) ']'],'parameters for Sea Ice solution'));
 		end % }}}
 		function memory(obj) % {{{
 
Index: ../trunk-jpl/src/m/enum/MaterialsYoungModulusEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/MaterialsYoungModulusEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/MaterialsYoungModulusEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=MaterialsYoungModulusEnum()
+%MATERIALSYOUNGMODULUSENUM - Enum of MaterialsYoungModulus
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=MaterialsYoungModulusEnum()
+
+macro=StringToEnum('MaterialsYoungModulus');
Index: ../trunk-jpl/src/m/enum/BasalforcingsOceanLinDragCoefEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BasalforcingsOceanLinDragCoefEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/BasalforcingsOceanLinDragCoefEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=BasalforcingsOceanLinDragCoefEnum()
+%BASALFORCINGSOCEANLINDRAGCOEFENUM - Enum of BasalforcingsOceanLinDragCoef
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=BasalforcingsOceanLinDragCoefEnum()
+
+macro=StringToEnum('BasalforcingsOceanLinDragCoef');
Index: ../trunk-jpl/src/m/enum/MaterialsCohesionEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/MaterialsCohesionEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/MaterialsCohesionEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=MaterialsCohesionEnum()
+%MATERIALSCOHESIONENUM - Enum of MaterialsCohesion
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=MaterialsCohesionEnum()
+
+macro=StringToEnum('MaterialsCohesion');
Index: ../trunk-jpl/src/m/enum/MaterialsDamageEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/MaterialsDamageEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/MaterialsDamageEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=MaterialsDamageEnum()
+%MATERIALSDAMAGEENUM - Enum of MaterialsDamage
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=MaterialsDamageEnum()
+
+macro=StringToEnum('MaterialsDamage');
Index: ../trunk-jpl/src/m/enum/BasalforcingsRhoOceanEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BasalforcingsRhoOceanEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/BasalforcingsRhoOceanEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=BasalforcingsRhoOceanEnum()
+%BASALFORCINGSRHOOCEANENUM - Enum of BasalforcingsRhoOcean
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=BasalforcingsRhoOceanEnum()
+
+macro=StringToEnum('BasalforcingsRhoOcean');
Index: ../trunk-jpl/src/m/enum/MaterialsInternalFrictionCoefEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/MaterialsInternalFrictionCoefEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/MaterialsInternalFrictionCoefEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=MaterialsInternalFrictionCoefEnum()
+%MATERIALSINTERNALFRICTIONCOEFENUM - Enum of MaterialsInternalFrictionCoef
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=MaterialsInternalFrictionCoefEnum()
+
+macro=StringToEnum('MaterialsInternalFrictionCoef');
Index: ../trunk-jpl/src/m/enum/SeaiceConcentrationEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SeaiceConcentrationEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SeaiceConcentrationEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SeaiceConcentrationEnum()
+%SEAICECONCENTRATIONENUM - Enum of SeaiceConcentration
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SeaiceConcentrationEnum()
+
+macro=StringToEnum('SeaiceConcentration');
Index: ../trunk-jpl/src/m/enum/SurfaceforcingsWindVxEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SurfaceforcingsWindVxEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SurfaceforcingsWindVxEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SurfaceforcingsWindVxEnum()
+%SURFACEFORCINGSWINDVXENUM - Enum of SurfaceforcingsWindVx
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SurfaceforcingsWindVxEnum()
+
+macro=StringToEnum('SurfaceforcingsWindVx');
Index: ../trunk-jpl/src/m/enum/VyStarEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/VyStarEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/VyStarEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=VyStarEnum()
+%VYSTARENUM - Enum of VyStar
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=VyStarEnum()
+
+macro=StringToEnum('VyStar');
Index: ../trunk-jpl/src/m/enum/SurfaceforcingsAirLinDragCoefEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SurfaceforcingsAirLinDragCoefEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SurfaceforcingsAirLinDragCoefEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SurfaceforcingsAirLinDragCoefEnum()
+%SURFACEFORCINGSAIRLINDRAGCOEFENUM - Enum of SurfaceforcingsAirLinDragCoef
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SurfaceforcingsAirLinDragCoefEnum()
+
+macro=StringToEnum('SurfaceforcingsAirLinDragCoef');
Index: ../trunk-jpl/src/m/enum/SeaiceSolutionEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SeaiceSolutionEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SeaiceSolutionEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SeaiceSolutionEnum()
+%SEAICESOLUTIONENUM - Enum of SeaiceSolution
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SeaiceSolutionEnum()
+
+macro=StringToEnum('SeaiceSolution');
Index: ../trunk-jpl/src/m/enum/MaterialsPoissonEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/MaterialsPoissonEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/MaterialsPoissonEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=MaterialsPoissonEnum()
+%MATERIALSPOISSONENUM - Enum of MaterialsPoisson
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=MaterialsPoissonEnum()
+
+macro=StringToEnum('MaterialsPoisson');
Index: ../trunk-jpl/src/m/enum/MaterialsCompressionCoefEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/MaterialsCompressionCoefEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/MaterialsCompressionCoefEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=MaterialsCompressionCoefEnum()
+%MATERIALSCOMPRESSIONCOEFENUM - Enum of MaterialsCompressionCoef
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=MaterialsCompressionCoefEnum()
+
+macro=StringToEnum('MaterialsCompressionCoef');
Index: ../trunk-jpl/src/m/enum/BasalforcingsOceanCoefEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BasalforcingsOceanCoefEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/BasalforcingsOceanCoefEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=BasalforcingsOceanCoefEnum()
+%BASALFORCINGSOCEANCOEFENUM - Enum of BasalforcingsOceanCoef
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=BasalforcingsOceanCoefEnum()
+
+macro=StringToEnum('BasalforcingsOceanCoef');
Index: ../trunk-jpl/src/m/enum/MaterialsTractionCoefEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/MaterialsTractionCoefEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/MaterialsTractionCoefEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=MaterialsTractionCoefEnum()
+%MATERIALSTRACTIONCOEFENUM - Enum of MaterialsTractionCoef
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=MaterialsTractionCoefEnum()
+
+macro=StringToEnum('MaterialsTractionCoef');
Index: ../trunk-jpl/src/m/enum/SurfaceforcingsAirCoefEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SurfaceforcingsAirCoefEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SurfaceforcingsAirCoefEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SurfaceforcingsAirCoefEnum()
+%SURFACEFORCINGSAIRCOEFENUM - Enum of SurfaceforcingsAirCoef
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SurfaceforcingsAirCoefEnum()
+
+macro=StringToEnum('SurfaceforcingsAirCoef');
Index: ../trunk-jpl/src/m/enum/BasalforcingsOceanVxEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BasalforcingsOceanVxEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/BasalforcingsOceanVxEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=BasalforcingsOceanVxEnum()
+%BASALFORCINGSOCEANVXENUM - Enum of BasalforcingsOceanVx
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=BasalforcingsOceanVxEnum()
+
+macro=StringToEnum('BasalforcingsOceanVx');
Index: ../trunk-jpl/src/m/enum/SurfaceforcingsWindVyEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SurfaceforcingsWindVyEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SurfaceforcingsWindVyEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SurfaceforcingsWindVyEnum()
+%SURFACEFORCINGSWINDVYENUM - Enum of SurfaceforcingsWindVy
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SurfaceforcingsWindVyEnum()
+
+macro=StringToEnum('SurfaceforcingsWindVy');
Index: ../trunk-jpl/src/m/enum/SeaiceatmEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SeaiceatmEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SeaiceatmEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SeaiceatmEnum()
+%SEAICEATMENUM - Enum of Seaiceatm
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SeaiceatmEnum()
+
+macro=StringToEnum('Seaiceatm');
Index: ../trunk-jpl/src/m/enum/SeaiceAnalysisEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SeaiceAnalysisEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SeaiceAnalysisEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SeaiceAnalysisEnum()
+%SEAICEANALYSISENUM - Enum of SeaiceAnalysis
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SeaiceAnalysisEnum()
+
+macro=StringToEnum('SeaiceAnalysis');
Index: ../trunk-jpl/src/m/enum/MaterialsRidgingExponentEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/MaterialsRidgingExponentEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/MaterialsRidgingExponentEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=MaterialsRidgingExponentEnum()
+%MATERIALSRIDGINGEXPONENTENUM - Enum of MaterialsRidgingExponent
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=MaterialsRidgingExponentEnum()
+
+macro=StringToEnum('MaterialsRidgingExponent');
Index: ../trunk-jpl/src/m/enum/BasalforcingsOceanVyEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BasalforcingsOceanVyEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/BasalforcingsOceanVyEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=BasalforcingsOceanVyEnum()
+%BASALFORCINGSOCEANVYENUM - Enum of BasalforcingsOceanVy
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=BasalforcingsOceanVyEnum()
+
+macro=StringToEnum('BasalforcingsOceanVy');
Index: ../trunk-jpl/src/m/enum/SeaiceThicknessEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SeaiceThicknessEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SeaiceThicknessEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SeaiceThicknessEnum()
+%SEAICETHICKNESSENUM - Enum of SeaiceThickness
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SeaiceThicknessEnum()
+
+macro=StringToEnum('SeaiceThickness');
Index: ../trunk-jpl/src/m/enum/BasalforcingsOceanSshEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BasalforcingsOceanSshEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/BasalforcingsOceanSshEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=BasalforcingsOceanSshEnum()
+%BASALFORCINGSOCEANSSHENUM - Enum of BasalforcingsOceanSsh
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=BasalforcingsOceanSshEnum()
+
+macro=StringToEnum('BasalforcingsOceanSsh');
Index: ../trunk-jpl/src/m/enum/SeaiceSpcvxEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SeaiceSpcvxEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SeaiceSpcvxEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SeaiceSpcvxEnum()
+%SEAICESPCVXENUM - Enum of SeaiceSpcvx
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SeaiceSpcvxEnum()
+
+macro=StringToEnum('SeaiceSpcvx');
Index: ../trunk-jpl/src/m/enum/BasalforcingsOceanQuadDragCoefEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BasalforcingsOceanQuadDragCoefEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/BasalforcingsOceanQuadDragCoefEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=BasalforcingsOceanQuadDragCoefEnum()
+%BASALFORCINGSOCEANQUADDRAGCOEFENUM - Enum of BasalforcingsOceanQuadDragCoef
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=BasalforcingsOceanQuadDragCoefEnum()
+
+macro=StringToEnum('BasalforcingsOceanQuadDragCoef');
Index: ../trunk-jpl/src/m/enum/SurfaceforcingsRhoAirEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SurfaceforcingsRhoAirEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SurfaceforcingsRhoAirEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SurfaceforcingsRhoAirEnum()
+%SURFACEFORCINGSRHOAIRENUM - Enum of SurfaceforcingsRhoAir
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SurfaceforcingsRhoAirEnum()
+
+macro=StringToEnum('SurfaceforcingsRhoAir');
Index: ../trunk-jpl/src/m/enum/SurfaceforcingsAirQuadDragCoefEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SurfaceforcingsAirQuadDragCoefEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SurfaceforcingsAirQuadDragCoefEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SurfaceforcingsAirQuadDragCoefEnum()
+%SURFACEFORCINGSAIRQUADDRAGCOEFENUM - Enum of SurfaceforcingsAirQuadDragCoef
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SurfaceforcingsAirQuadDragCoefEnum()
+
+macro=StringToEnum('SurfaceforcingsAirQuadDragCoef');
Index: ../trunk-jpl/src/m/enum/SeaiceoceanEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SeaiceoceanEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SeaiceoceanEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SeaiceoceanEnum()
+%SEAICEOCEANENUM - Enum of Seaiceocean
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SeaiceoceanEnum()
+
+macro=StringToEnum('Seaiceocean');
Index: ../trunk-jpl/src/m/enum/BasalforcingsOceanTurningAngleEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/BasalforcingsOceanTurningAngleEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/BasalforcingsOceanTurningAngleEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=BasalforcingsOceanTurningAngleEnum()
+%BASALFORCINGSOCEANTURNINGANGLEENUM - Enum of BasalforcingsOceanTurningAngle
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=BasalforcingsOceanTurningAngleEnum()
+
+macro=StringToEnum('BasalforcingsOceanTurningAngle');
Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18491)
+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18492)
@@ -381,6 +381,8 @@
 def UzawaPressureAnalysisEnum(): return StringToEnum("UzawaPressureAnalysis")[0]
 def GiaSolutionEnum(): return StringToEnum("GiaSolution")[0]
 def GiaAnalysisEnum(): return StringToEnum("GiaAnalysis")[0]
+def SeaiceSolutionEnum(): return StringToEnum("SeaiceSolution")[0]
+def SeaiceAnalysisEnum(): return StringToEnum("SeaiceAnalysis")[0]
 def MeshdeformationSolutionEnum(): return StringToEnum("MeshdeformationSolution")[0]
 def MeshdeformationAnalysisEnum(): return StringToEnum("MeshdeformationAnalysis")[0]
 def LevelsetAnalysisEnum(): return StringToEnum("LevelsetAnalysis")[0]
@@ -693,4 +695,35 @@
 def LevelsetfunctionSlopeXEnum(): return StringToEnum("LevelsetfunctionSlopeX")[0]
 def LevelsetfunctionSlopeYEnum(): return StringToEnum("LevelsetfunctionSlopeY")[0]
 def LevelsetfunctionPicardEnum(): return StringToEnum("LevelsetfunctionPicard")[0]
+def SeaiceatmEnum(): return StringToEnum("Seaiceatm")[0]
+def SeaiceoceanEnum(): return StringToEnum("Seaiceocean")[0]
+def SeaiceThicknessEnum(): return StringToEnum("SeaiceThickness")[0]
+def SeaiceConcentrationEnum(): return StringToEnum("SeaiceConcentration")[0]
+def SeaiceSpcvxEnum(): return StringToEnum("SeaiceSpcvx")[0]
+def SeaiceSpcvyEnum(): return StringToEnum("SeaiceSpcvy")[0]
+def BasalforcingsRhoOceanEnum(): return StringToEnum("BasalforcingsRhoOcean")[0]
+def BasalforcingsOceanCoefEnum(): return StringToEnum("BasalforcingsOceanCoef")[0]
+def BasalforcingsOceanLinDragCoefEnum(): return StringToEnum("BasalforcingsOceanLinDragCoef")[0]
+def BasalforcingsOceanQuadDragCoefEnum(): return StringToEnum("BasalforcingsOceanQuadDragCoef")[0]
+def BasalforcingsOceanTurningAngleEnum(): return StringToEnum("BasalforcingsOceanTurningAngle")[0]
+def BasalforcingsOceanSshEnum(): return StringToEnum("BasalforcingsOceanSsh")[0]
+def BasalforcingsOceanVxEnum(): return StringToEnum("BasalforcingsOceanVx")[0]
+def BasalforcingsOceanVyEnum(): return StringToEnum("BasalforcingsOceanVy")[0]
+def SurfaceforcingsRhoAirEnum(): return StringToEnum("SurfaceforcingsRhoAir")[0]
+def SurfaceforcingsAirCoefEnum(): return StringToEnum("SurfaceforcingsAirCoef")[0]
+def SurfaceforcingsAirLinDragCoefEnum(): return StringToEnum("SurfaceforcingsAirLinDragCoef")[0]
+def SurfaceforcingsAirQuadDragCoefEnum(): return StringToEnum("SurfaceforcingsAirQuadDragCoef")[0]
+def SurfaceforcingsWindVxEnum(): return StringToEnum("SurfaceforcingsWindVx")[0]
+def SurfaceforcingsWindVyEnum(): return StringToEnum("SurfaceforcingsWindVy")[0]
+def MatseaiceEnum(): return StringToEnum("Matseaice")[0]
+def MaterialsPoissonEnum(): return StringToEnum("MaterialsPoisson")[0]
+def MaterialsYoungModulusEnum(): return StringToEnum("MaterialsYoungModulus")[0]
+def MaterialsDamageEnum(): return StringToEnum("MaterialsDamage")[0]
+def MaterialsRidgingExponentEnum(): return StringToEnum("MaterialsRidgingExponent")[0]
+def MaterialsCohesionEnum(): return StringToEnum("MaterialsCohesion")[0]
+def MaterialsInternalFrictionCoefEnum(): return StringToEnum("MaterialsInternalFrictionCoef")[0]
+def MaterialsCompressionCoefEnum(): return StringToEnum("MaterialsCompressionCoef")[0]
+def MaterialsTractionCoefEnum(): return StringToEnum("MaterialsTractionCoef")[0]
+def VxStarEnum(): return StringToEnum("VxStar")[0]
+def VyStarEnum(): return StringToEnum("VyStar")[0]
 def MaximumNumberOfDefinitionsEnum(): return StringToEnum("MaximumNumberOfDefinitions")[0]
Index: ../trunk-jpl/src/m/enum/SeaiceSpcvyEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/SeaiceSpcvyEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/SeaiceSpcvyEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=SeaiceSpcvyEnum()
+%SEAICESPCVYENUM - Enum of SeaiceSpcvy
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SeaiceSpcvyEnum()
+
+macro=StringToEnum('SeaiceSpcvy');
Index: ../trunk-jpl/src/m/enum/VxStarEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/VxStarEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/VxStarEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=VxStarEnum()
+%VXSTARENUM - Enum of VxStar
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=VxStarEnum()
+
+macro=StringToEnum('VxStar');
Index: ../trunk-jpl/src/m/enum/MatseaiceEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/MatseaiceEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/MatseaiceEnum.m	(revision 18492)
@@ -0,0 +1,11 @@
+function macro=MatseaiceEnum()
+%MATSEAICEENUM - Enum of Matseaice
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=MatseaiceEnum()
+
+macro=StringToEnum('Matseaice');
