Index: issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15400)
+++ issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15401)
@@ -2178,4 +2178,5 @@
 				name==FrictionCoefficientEnum ||
 				name==GLlevelsetEnum ||
+				name==IcelevelsetEnum ||
 				name==GradientEnum ||
 				name==OldGradientEnum  ||
Index: issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15400)
+++ issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15401)
@@ -1733,4 +1733,5 @@
 				name==BedEnum ||
 				name==GLlevelsetEnum ||
+				name==IcelevelsetEnum ||
 				name==SurfaceSlopeXEnum ||
 				name==SurfaceSlopeYEnum ||
@@ -2942,4 +2943,18 @@
 ElementVector* Tria::CreatePVectorDiagnosticMacAyeal(){
 
+	/*compute all load vectors for this element*/
+	ElementVector* pe1=CreatePVectorDiagnosticMacAyealDrivingStress();
+	ElementVector* pe2=CreatePVectorDiagnosticMacAyealFront();
+	ElementVector* pe =new ElementVector(pe1,pe2);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	return pe;
+}
+/*}}}*/
+/*FUNCTION Tria::CreatePVectorDiagnosticMacAyealDrivingStress {{{*/
+ElementVector* Tria::CreatePVectorDiagnosticMacAyealDrivingStress(){
+
 	/*Intermediaries */
 	int            i,j;
@@ -2948,4 +2963,5 @@
 	IssmDouble     xyz_list[NUMVERTICES][3];
 	IssmDouble     slope[2];
+	IssmDouble     icefrontlevel[3];
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -2963,4 +2979,5 @@
 	Input* surface_input=inputs->GetInput(SurfaceEnum);     _assert_(surface_input);
 	Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);_assert_(drag_input);
+	GetInputListOnVertices(&icefrontlevel[0],BedEnum);
 
 	/* Start  looping on the number of gaussian points: */
@@ -2981,4 +2998,41 @@
 				pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
 			}
+		}
+	}
+
+	/*Transform coordinate system*/
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYEnum);
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(basis);
+	return pe;
+}
+/*}}}*/
+/*FUNCTION Tria::CreatePVectorDiagnosticMacAyealFront {{{*/
+ElementVector* Tria::CreatePVectorDiagnosticMacAyealFront(){
+
+	/*Intermediaries */
+	int            i,j;
+	IssmDouble     ls[3];
+	IssmDouble     xyz_list[NUMVERTICES][3];
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+	int numdof   = numnodes*NDOF2;
+	Icefront *icefront=NULL;
+
+	return NULL;
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	GetInputListOnVertices(&ls[0],IcelevelsetEnum);
+	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
+	GaussTria*     gauss  = new GaussTria(2);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Create Ice Front if necessary*/
+	if(ls[0]>0. || ls[1]>0. || ls[2]>0.){
+		if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]==0.)){
+			//icefront=new Icefront("2d",inputs,matpar,MacAyealApproximationEnum,analysis_type);
 		}
 	}
Index: issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 15400)
+++ issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 15401)
@@ -224,4 +224,6 @@
 		ElementMatrix* CreateKMatrixDiagnosticHutter(void);
 		ElementVector* CreatePVectorDiagnosticMacAyeal(void);
+		ElementVector* CreatePVectorDiagnosticMacAyealDrivingStress(void);
+		ElementVector* CreatePVectorDiagnosticMacAyealFront(void);
 		ElementVector* CreatePVectorDiagnosticHutter(void);
 		ElementMatrix* CreateJacobianDiagnosticMacayeal(void);
Index: issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp	(revision 15400)
+++ issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp	(revision 15401)
@@ -122,4 +122,87 @@
 	this->element    = NULL;
 	this->matpar     = NULL;
+}
+
+/*}}}*/
+/*FUNCTION Icefront::Icefront(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in, int icefront_type, int in_analysis_type) {{{*/
+Icefront::Icefront(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in,int in_icefront_type,  int in_analysis_type){
+
+	int segment_width;
+	int element;
+	int numnodes; 
+	int numvertices; 
+	int dim;
+	int numberofelements;
+
+	/*icefront constructor data: */
+	int  icefront_eid;
+	int  icefront_mparid;
+	int  icefront_node_ids[NUMVERTICESQUA]; //initialize with largest size
+	int  icefront_vertex_ids[NUMVERTICESQUA]; //initialize with largest size
+
+//	/*find parameters: */
+//	iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
+//
+	/*First, retrieve element index and element type: */
+	if(strcmp(element_type_in,"2d")==0){
+		segment_width=4;
+	}
+	else{
+		segment_width=6;
+	}
+//	element=reCast<int,IssmDouble>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)-1); //element is in the penultimate column (node1 node2 ... elem fill)
+//
+//	/*Build ids for hook constructors: */
+//	icefront_eid=reCast<int,IssmDouble>( *(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)); //matlab indexing
+//	icefront_mparid=numberofelements+1; //matlab indexing
+//
+	if (in_icefront_type==MacAyeal2dIceFrontEnum || in_icefront_type==MacAyeal3dIceFrontEnum){
+//		icefront_node_ids[0]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
+//		icefront_node_ids[1]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
+//		icefront_vertex_ids[0]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
+//		icefront_vertex_ids[1]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
+	}
+	else if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
+//		icefront_node_ids[0]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
+//		icefront_node_ids[1]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
+//		icefront_node_ids[2]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2));
+//		icefront_node_ids[3]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3));
+//		icefront_vertex_ids[0]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
+//		icefront_vertex_ids[1]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
+//		icefront_vertex_ids[2]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2));
+//		icefront_vertex_ids[3]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3));
+	}
+	else _error_("in_icefront_type " << EnumToStringx(in_icefront_type) << " not supported yet!");
+
+	if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
+		numnodes=4;
+		numvertices=4;
+	}
+	else{
+		numnodes=2;
+		numvertices=2;
+	}
+
+	/*Ok, we have everything to build the object: */
+	this->id=1;
+	this->analysis_type=in_analysis_type;
+
+	/*Hooks: */
+	this->hnodes=new Hook(icefront_node_ids,numnodes);
+	this->hvertices=new Hook(icefront_vertex_ids,numvertices);
+	this->helement=new Hook(&icefront_eid,1);
+	this->hmatpar=new Hook(&icefront_mparid,1);
+
+	//intialize  and add as many inputs per element as requested: 
+	this->inputs=inputs_in;
+	this->inputs->AddInput(new IntInput(FillEnum,1)); //We always consider we have water, if above sea level, only air will be applied
+	this->inputs->AddInput(new IntInput(IceFrontTypeEnum,in_icefront_type)); 
+
+	//parameters and hooked fields: we still can't point to them, they may not even exist. Configure will handle this.
+	this->parameters = NULL;
+	this->nodes      = NULL;
+	this->vertices   = NULL;
+	this->element    = NULL;
+	this->matpar     = matpar_in;
 }
 
Index: issm/trunk-jpl/src/c/classes/Loads/Icefront.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Loads/Icefront.h	(revision 15400)
+++ issm/trunk-jpl/src/c/classes/Loads/Icefront.h	(revision 15401)
@@ -45,4 +45,5 @@
 		Icefront();
 		Icefront(int icefront_id,int i, IoModel* iomodel,int in_icefront_type, int analysis_type);
+		Icefront(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in, int icefront_type, int in_analysis_type);
 		~Icefront();
 		/*}}}*/
Index: issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
===================================================================
--- issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 15400)
+++ issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 15401)
@@ -92,5 +92,5 @@
 
 	/*Fetch data:*/
-	iomodel->FetchData(6,MeshElementsEnum,MeshXEnum,MeshYEnum,MeshZEnum,BedEnum,ThicknessEnum);
+	iomodel->FetchData(7,MeshElementsEnum,MeshXEnum,MeshYEnum,MeshZEnum,BedEnum,ThicknessEnum,IcelevelsetEnum);
 	CreateNumberNodeToElementConnectivity(iomodel);
 
Index: issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 15400)
+++ issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 15401)
@@ -144,4 +144,5 @@
 	MaskVertexongroundediceEnum,
 	MaskVertexonwaterEnum,
+	IcelevelsetEnum,
 	MaterialsBetaEnum,
 	MaterialsHeatcapacityEnum,
Index: issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 15400)
+++ issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 15401)
@@ -152,4 +152,5 @@
 		case MaskVertexongroundediceEnum : return "MaskVertexongroundedice";
 		case MaskVertexonwaterEnum : return "MaskVertexonwater";
+		case IcelevelsetEnum : return "Icelevelset";
 		case MaterialsBetaEnum : return "MaterialsBeta";
 		case MaterialsHeatcapacityEnum : return "MaterialsHeatcapacity";
Index: issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 15400)
+++ issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 15401)
@@ -155,4 +155,5 @@
 	      else if (strcmp(name,"MaskVertexongroundedice")==0) return MaskVertexongroundediceEnum;
 	      else if (strcmp(name,"MaskVertexonwater")==0) return MaskVertexonwaterEnum;
+	      else if (strcmp(name,"Icelevelset")==0) return IcelevelsetEnum;
 	      else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
 	      else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
@@ -259,9 +260,9 @@
 	      else if (strcmp(name,"GiaMantleViscosity")==0) return GiaMantleViscosityEnum;
 	      else if (strcmp(name,"GiaLithosphereThickness")==0) return GiaLithosphereThicknessEnum;
-	      else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
          else stage=3;
    }
    if(stage==3){
-	      if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum;
+	      if (strcmp(name,"Thickness")==0) return ThicknessEnum;
+	      else if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum;
 	      else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum;
 	      else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;
@@ -382,9 +383,9 @@
 	      else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
 	      else if (strcmp(name,"RiftfrontType")==0) return RiftfrontTypeEnum;
-	      else if (strcmp(name,"Segment")==0) return SegmentEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
+	      if (strcmp(name,"Segment")==0) return SegmentEnum;
+	      else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
 	      else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
 	      else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
@@ -505,9 +506,9 @@
 	      else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
 	      else if (strcmp(name,"DoubleElementResult")==0) return DoubleElementResultEnum;
-	      else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
+	      if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
+	      else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
 	      else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
 	      else if (strcmp(name,"J")==0) return JEnum;
Index: issm/trunk-jpl/src/m/classes/mask.m
===================================================================
--- issm/trunk-jpl/src/m/classes/mask.m	(revision 15400)
+++ issm/trunk-jpl/src/m/classes/mask.m	(revision 15401)
@@ -54,4 +54,8 @@
 			WriteData(fid,'object',obj,'fieldname','vertexongroundedice','format','DoubleMat','mattype',1);
 			WriteData(fid,'object',obj,'fieldname','vertexonwater','format','DoubleMat','mattype',1);
+			icelevelset=ones(md.mesh.numberofvertices,1);
+			pos=md.diagnostic.icefront(:,1:end-1);
+			icelevelset(pos(:))=0;
+			WriteData(fid,'data',icelevelset,'format','DoubleMat','mattype',1,'enum',IcelevelsetEnum());
 		end % }}}
 	end
Index: issm/trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- issm/trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 15400)
+++ issm/trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 15401)
@@ -1899,4 +1899,18 @@
 	return StringToEnum('MaskVertexonwater')[0]
 
+def IcelevelsetEnum():
+	"""
+	ICELEVELSETENUM - Enum of Icelevelset
+
+	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=IcelevelsetEnum()
+	"""
+
+	return StringToEnum('Icelevelset')[0]
+
 def MaterialsBetaEnum():
 	"""
@@ -7861,4 +7875,4 @@
 	"""
 
-	return 560
-
+	return 561
+
Index: issm/trunk-jpl/src/m/enum/IcelevelsetEnum.m
===================================================================
--- issm/trunk-jpl/src/m/enum/IcelevelsetEnum.m	(revision 15401)
+++ issm/trunk-jpl/src/m/enum/IcelevelsetEnum.m	(revision 15401)
@@ -0,0 +1,11 @@
+function macro=IcelevelsetEnum()
+%ICELEVELSETENUM - Enum of Icelevelset
+%
+%   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=IcelevelsetEnum()
+
+macro=StringToEnum('Icelevelset');
Index: issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m
===================================================================
--- issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m	(revision 15400)
+++ issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m	(revision 15401)
@@ -9,3 +9,3 @@
 %      macro=MaximumNumberOfEnums()
 
-macro=560;
+macro=561;
