Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 5714)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 5715)
@@ -115,5 +115,6 @@
 	RiftfrontEnum,
 	SegmentRiftfrontEnum,
-	MacAyealIceFrontEnum,
+	MacAyeal2dIceFrontEnum,
+	MacAyeal3dIceFrontEnum,
 	PattynIceFrontEnum,
 	StokesIceFrontEnum,
Index: /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp	(revision 5714)
+++ /issm/trunk/src/c/EnumDefinitions/EnumToString.cpp	(revision 5715)
@@ -100,5 +100,6 @@
 		case RiftfrontEnum : return "Riftfront";
 		case SegmentRiftfrontEnum : return "SegmentRiftfront";
-		case MacAyealIceFrontEnum : return "MacAyealIceFront";
+		case MacAyeal2dIceFrontEnum : return "MacAyeal2dIceFront";
+		case MacAyeal3dIceFrontEnum : return "MacAyeal3dIceFront";
 		case PattynIceFrontEnum : return "PattynIceFront";
 		case StokesIceFrontEnum : return "StokesIceFront";
Index: /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp	(revision 5714)
+++ /issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp	(revision 5715)
@@ -98,5 +98,6 @@
 	else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
 	else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
-	else if (strcmp(name,"MacAyealIceFront")==0) return MacAyealIceFrontEnum;
+	else if (strcmp(name,"MacAyeal2dIceFront")==0) return MacAyeal2dIceFrontEnum;
+	else if (strcmp(name,"MacAyeal3dIceFront")==0) return MacAyeal3dIceFrontEnum;
 	else if (strcmp(name,"PattynIceFront")==0) return PattynIceFrontEnum;
 	else if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum;
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp	(revision 5714)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp	(revision 5715)
@@ -60,6 +60,10 @@
 
 		/*Create and  add load: */
-		if ((int)*(iomodel->elements_type+element)==(MacAyealApproximationEnum)){
-			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyealIceFrontEnum,DiagnosticHorizAnalysisEnum));
+		if ((int)*(iomodel->elements_type+element)==(MacAyealApproximationEnum) && iomodel->dim==2){
+			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal2dIceFrontEnum,DiagnosticHorizAnalysisEnum));
+			count++;
+		}
+		else if ((int)*(iomodel->elements_type+element)==(MacAyealApproximationEnum) && iomodel->dim==3){
+			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
 			count++;
 		}
@@ -73,5 +77,5 @@
 		}
 		else if ((int)*(iomodel->elements_type+element)==(MacAyealPattynApproximationEnum)){
-			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyealIceFrontEnum,DiagnosticHorizAnalysisEnum));
+			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
 			count++;
 			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,PattynIceFrontEnum,DiagnosticHorizAnalysisEnum));
Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5714)
+++ /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5715)
@@ -68,5 +68,5 @@
 	icefront_mparid=iomodel->numberofelements+1; //matlab indexing
 
-	if (in_icefront_type==MacAyealIceFrontEnum){
+	if (in_icefront_type==MacAyeal2dIceFrontEnum || in_icefront_type==MacAyeal3dIceFrontEnum){
 		icefront_node_ids[0]=iomodel->nodecounter+(int)*(iomodel->pressureload+segment_width*i+0);
 		icefront_node_ids[1]=iomodel->nodecounter+(int)*(iomodel->pressureload+segment_width*i+1);
@@ -330,6 +330,9 @@
 		int type;
 		inputs->GetParameterValue(&type,TypeEnum);
-		if (type==MacAyealIceFrontEnum){
-			CreatePVectorDiagnosticMacAyeal( pg);
+		if (type==MacAyeal2dIceFrontEnum){
+			CreatePVectorDiagnosticMacAyeal2d( pg);
+		}
+		else if (type==MacAyeal3dIceFrontEnum){
+			CreatePVectorDiagnosticMacAyeal3d( pg);
 		}
 		else if (type==PattynIceFrontEnum){
@@ -419,6 +422,6 @@
 
 /*Icefront numerics: */
-/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal{{{1*/
-void Icefront::CreatePVectorDiagnosticMacAyeal( Vec pg){
+/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal2d{{{1*/
+void Icefront::CreatePVectorDiagnosticMacAyeal2d( Vec pg){
 
 	int       i;
@@ -450,20 +453,10 @@
 	Input*  bed_input=NULL;
 	Tria*   tria=NULL;
-	Penta*  penta=NULL;
 
 	BeamRef* beam=NULL;
 
-	//check that the element is onbed (collapsed formulation) otherwise:pe=0
-	if(element->Enum()==PentaEnum){
-		if(!element->GetOnBed()) return;
-		thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum);
-		bed_input      =((Penta*)element)->inputs->GetInput(BedEnum);
-	}
-	else{
-		thickness_input=((Tria*)element)->inputs->GetInput(ThicknessEnum);
-		bed_input      =((Tria*)element)->inputs->GetInput(BedEnum);
-	}
-
 	/*Recover inputs: */
+	thickness_input=((Tria*)element)->inputs->GetInput(ThicknessEnum);
+	bed_input      =((Tria*)element)->inputs->GetInput(BedEnum);
 	inputs->GetParameterValue(&fill,FillEnum);
 
@@ -514,5 +507,5 @@
 		}
 		else{
-			ISSMERROR("fill type %i not supported yet",fill);
+			ISSMERROR("fill type %s not supported yet",EnumToString(fill));
 		}
 		pressure = ice_pressure + water_pressure + air_pressure;
@@ -528,4 +521,116 @@
 
 			
+	/*Plug pe_g into vector: */
+	VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
+
+	/*Free ressources:*/
+	xfree((void**)&segment_gauss_coord);
+	xfree((void**)&gauss_weights);
+	xfree((void**)&doflist);
+
+}
+/*}}}*/
+/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal3d{{{1*/
+void Icefront::CreatePVectorDiagnosticMacAyeal3d( Vec pg){
+
+	int       i;
+	const int numnodes= NUMVERTICESSEG;
+	const int numdofs = numnodes *NDOF2;
+	int*      doflist=NULL;
+	double    xyz_list[numnodes][3];
+
+	/*Objects: */
+	double    pe_g[numdofs] = {0.0};
+
+	/*Segment*/
+	double   normal[2];
+	int      num_gauss;
+	double*  segment_gauss_coord=NULL;
+	double*  gauss_weights=NULL;
+	double   gauss_weight;
+	double   gauss_coord;
+	int      ig;
+	double   Jdet;
+	double   L[2];
+	double   thickness,bed,pressure;
+	double   ice_pressure,water_pressure,air_pressure,rho_water,rho_ice,gravity;
+	double   surface_under_water,base_under_water;
+	int      fill;
+
+	/*Inputs*/
+	Input*  thickness_input=NULL;
+	Input*  bed_input=NULL;
+	Penta*  penta=NULL;
+
+	BeamRef* beam=NULL;
+
+	//check that the element is onbed (collapsed formulation) otherwise:pe=0
+	if(!element->GetOnBed()) return;
+
+	/*Recover inputs: */
+	thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum);
+	bed_input      =((Penta*)element)->inputs->GetInput(BedEnum);
+	inputs->GetParameterValue(&fill,FillEnum);
+
+	/* Get dof list and node coordinates: */
+	GetDofList(&doflist,MacAyealApproximationEnum);
+	GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes);
+
+	/*Get Matpar parameters*/
+	rho_water=matpar->GetRhoWater();
+	rho_ice  =matpar->GetRhoIce();
+	gravity  =matpar->GetG();
+
+	/*Get normal*/
+	GetNormal(&normal[0],xyz_list);
+
+	//Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
+	num_gauss=3;
+	gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
+
+	for(ig=0;ig<num_gauss;ig++){
+
+		/*Pick up the gaussian point: */
+		gauss_weight=gauss_weights[ig];
+		gauss_coord=segment_gauss_coord[ig];
+
+		//compute Jacobian of segment
+		beam->GetJacobianDeterminant2d(&Jdet,&xyz_list[0][0],gauss_coord);
+
+		/*Get thickness and bed*/
+		this->GetParameterValue(&thickness,gauss_coord,thickness_input);
+		this->GetParameterValue(&bed,      gauss_coord,bed_input);
+
+		/*Compute total pressure applied to ice front*/
+		if (fill==WaterEnum){
+			//icefront ends in water: 
+			ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
+			air_pressure=0;
+
+			//Now deal with water pressure
+			surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level
+			base_under_water=min(0,bed);              // 0 if the bottom of the glacier is above water level
+			water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
+		}
+		else if (fill==AirEnum){
+			ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
+			air_pressure=0;
+			water_pressure=0;
+		}
+		else{
+			ISSMERROR("fill type %i not supported yet",fill);
+		}
+		pressure = ice_pressure + water_pressure + air_pressure;
+
+		/*Get L*/
+		beam->GetNodalFunctions(&L[0],gauss_coord);
+
+		for (i=0;i<numnodes;i++){
+			pe_g[2*i+0]+= pressure*Jdet*gauss_weights[ig]*normal[0]*L[i];
+			pe_g[2*i+1]+= pressure*Jdet*gauss_weights[ig]*normal[1]*L[i];
+		}
+	} //for(ig=0;ig<num_gauss;ig++)
+
+
 	/*Plug pe_g into vector: */
 	VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
@@ -818,6 +923,8 @@
 		
 	/*How many nodes? :*/
-	if(type==MacAyealIceFrontEnum)numberofnodes=2;
-	else numberofnodes=4;
+	if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum)
+	 numberofnodes=2;
+	else 
+	 numberofnodes=4;
 	
 	/*Figure out size of doflist: */
Index: /issm/trunk/src/c/objects/Loads/Icefront.h
===================================================================
--- /issm/trunk/src/c/objects/Loads/Icefront.h	(revision 5714)
+++ /issm/trunk/src/c/objects/Loads/Icefront.h	(revision 5715)
@@ -76,7 +76,8 @@
 		/*}}}*/
 		/*Load management: {{{1*/
-		void  CreatePVectorDiagnosticMacAyeal( Vec pg);
-		void  CreatePVectorDiagnosticPattyn( Vec pg);
-		void  CreatePVectorDiagnosticStokes( Vec pg);
+		void  CreatePVectorDiagnosticMacAyeal2d(Vec pg);
+		void  CreatePVectorDiagnosticMacAyeal3d(Vec pg);
+		void  CreatePVectorDiagnosticPattyn(Vec pg);
+		void  CreatePVectorDiagnosticStokes(Vec pg);
 		void  GetDofList(int** pdoflist,int approximation_enum=0);
 		void  QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,