Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp	(revision 2958)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp	(revision 2959)
@@ -108,6 +108,6 @@
 		#endif
 	
-		/*Do not create ice front if Hutter elements*/
-		if ((int)*(iomodel->elements_type+2*element+0)==HutterFormulationEnum()){ 
+		/*Do not create ice front if Hutter or Stokes elements*/
+		if ((int)*(iomodel->elements_type+2*element+0)==(HutterFormulationEnum() || StokesFormulationEnum())){ 
 			continue;
 		}
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp	(revision 2958)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp	(revision 2959)
@@ -51,5 +51,5 @@
 	int pengrid_stabilize_constraints=0;
 
-	int numberofpressureloads_stokes;
+	int numberofpressureloads;
 	int count;
 
@@ -63,5 +63,5 @@
 	/*Create pressure loads as boundary conditions. Pay attention to the partitioning if we are running in parallel (the grids
 	 * referenced by a certain load must belong to the cluster node): */
-	IoModelFetchData(&iomodel->pressureload_stokes,&numberofpressureloads_stokes,NULL,iomodel_handle,"pressureload_stokes");
+	IoModelFetchData(&iomodel->pressureload,&numberofpressureloads,NULL,iomodel_handle,"pressureload");
 	IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type");
 	IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
@@ -71,9 +71,9 @@
 
 	/*First load data:*/
-	for (i=0;i<numberofpressureloads_stokes;i++){
+	for (i=0;i<numberofpressureloads;i++){
 		
 		segment_width=5;
 
-		element=(int)(*(iomodel->pressureload_stokes+segment_width*i+segment_width-1)-1); //element is in the last column
+		element=(int)(*(iomodel->pressureload+segment_width*i+segment_width-1)-1); //element is in the last column
 
 		#ifdef _PARALLEL_
@@ -85,13 +85,18 @@
 		#endif
 	
+		/*Do not create ice front if it is not a Stokes element*/
+		if (!(int)*(iomodel->elements_type+2*element+1)==StokesFormulationEnum()){ 
+			continue;
+		}
+
 		icefront_mparid=iomodel->numberofelements+1; //matlab indexing
 		icefront_sid=count+1; //matlab indexing
-		icefront_eid=(int)*(iomodel->pressureload_stokes+segment_width*i+segment_width-1); //matlab indexing
+		icefront_eid=(int)*(iomodel->pressureload+segment_width*i+segment_width-1); //matlab indexing
 		icefront_element_type=PentaEnum();
 
-		i1=(int)*(iomodel->pressureload_stokes+segment_width*i+0);
-		i2=(int)*(iomodel->pressureload_stokes+segment_width*i+1);
-		i3=(int)*(iomodel->pressureload_stokes+segment_width*i+2);
-		i4=(int)*(iomodel->pressureload_stokes+segment_width*i+3);
+		i1=(int)*(iomodel->pressureload+segment_width*i+0);
+		i2=(int)*(iomodel->pressureload+segment_width*i+1);
+		i3=(int)*(iomodel->pressureload+segment_width*i+2);
+		i4=(int)*(iomodel->pressureload+segment_width*i+3);
 	
 		icefront_node_ids[0]=i1;
@@ -117,5 +122,5 @@
 	}
 	/*Free data: */
-	xfree((void**)&iomodel->pressureload_stokes);
+	xfree((void**)&iomodel->pressureload);
 	xfree((void**)&iomodel->elements_type);
 	xfree((void**)&iomodel->thickness);
Index: /issm/trunk/src/c/ModelProcessorx/IoModel.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/IoModel.cpp	(revision 2958)
+++ /issm/trunk/src/c/ModelProcessorx/IoModel.cpp	(revision 2959)
@@ -91,7 +91,5 @@
 	
 	iomodel->numberofpressureloads=0;
-	iomodel->numberofpressureloads_stokes=0;
 	iomodel->pressureload=NULL;
-	iomodel->pressureload_stokes=NULL;
 	iomodel-> spcvelocity=NULL;
 	iomodel-> spctemperature=NULL;
@@ -243,5 +241,4 @@
 	xfree((void**)&iomodel->gridoniceshelf);
 	xfree((void**)&iomodel->pressureload);
-	xfree((void**)&iomodel->pressureload_stokes);
 	xfree((void**)&iomodel->spcvelocity);
 	xfree((void**)&iomodel->spcthickness);
Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 2958)
+++ /issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 2959)
@@ -90,7 +90,5 @@
 	/*boundary conditions: */
 	int     numberofpressureloads;
-	int     numberofpressureloads_stokes;
 	double* pressureload;
-	double* pressureload_stokes;
 	double* spcvelocity;
 	double* spctemperature;
Index: /issm/trunk/src/m/classes/@model/model.m
===================================================================
--- /issm/trunk/src/m/classes/@model/model.m	(revision 2958)
+++ /issm/trunk/src/m/classes/@model/model.m	(revision 2959)
@@ -141,5 +141,4 @@
 	md.gridonboundary=NaN;
 	md.pressureload=NaN;
-	md.pressureload_stokes=NaN;
 	md.spcvelocity=NaN;
 	md.spctemperature=NaN;
Index: /issm/trunk/src/m/classes/public/marshall.m
===================================================================
--- /issm/trunk/src/m/classes/public/marshall.m	(revision 2958)
+++ /issm/trunk/src/m/classes/public/marshall.m	(revision 2959)
@@ -82,5 +82,4 @@
 
 WriteData(fid,md.pressureload,'Mat','pressureload');
-WriteData(fid,md.pressureload_stokes,'Mat','pressureload_stokes');
 
 WriteData(fid,md.geothermalflux,'Mat','geothermalflux');
Index: /issm/trunk/src/m/classes/public/modelextract.m
===================================================================
--- /issm/trunk/src/m/classes/public/modelextract.m	(revision 2958)
+++ /issm/trunk/src/m/classes/public/modelextract.m	(revision 2959)
@@ -234,14 +234,4 @@
 		end
 		md2.pressureload=md2.pressureload(find(md2.pressureload(:,1) & md2.pressureload(:,2) & md2.pressureload(:,end)),:);
-	end
-	if ~isnan(md2.pressureload_stokes)
-		md2.pressureload_stokes(:,1)=Pgrid(md1.pressureload_stokes(:,1)); 
-		md2.pressureload_stokes(:,2)=Pgrid(md1.pressureload_stokes(:,2)); 
-		md2.pressureload_stokes(:,end)=Pelem(md1.pressureload_stokes(:,end)); 
-		if strcmpi(md1.type,'3d')
-			md2.pressureload_stokes(:,3)=Pgrid(md1.pressureload_stokes(:,3)); 
-			md2.pressureload_stokes(:,4)=Pgrid(md1.pressureload_stokes(:,4)); 
-		end
-		md2.pressureload_stokes=md2.pressureload_stokes(find(md2.pressureload_stokes(:,1) & md2.pressureload_stokes(:,2) &  md2.pressureload_stokes(:,end)),:);
 	end
 
Index: /issm/trunk/src/m/classes/public/setelementstype.m
===================================================================
--- /issm/trunk/src/m/classes/public/setelementstype.m	(revision 2958)
+++ /issm/trunk/src/m/classes/public/setelementstype.m	(revision 2959)
@@ -126,9 +126,4 @@
 md.borderstokes=borderflags-md.gridonstokes; %remove stokes grids from this list
 
-%icefront for stokes
-if ~isnan(md.pressureload),
-	md.pressureload_stokes=md.pressureload(find(md.elements_type(md.pressureload(:,end),2)==StokesFormulationEnum()),:);
-end
-
 %figure out solution types
 md.ishutter=double(any(md.elements_type(:,1)==HutterFormulationEnum));