Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 1831)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 1832)
@@ -11,8 +11,8 @@
 #include "../../objects/objects.h"
 #include "../../shared/shared.h"
-#include "../Model.h"
+#include "../IoModel.h"
 
 
-void	CreateConstraintsDiagnosticHoriz(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){
+void	CreateConstraintsDiagnosticHoriz(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
 
@@ -44,17 +44,17 @@
 
 	/*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
-	if (!model->ismacayealpattyn)goto cleanup_and_return;
+	if (!iomodel->ismacayealpattyn)goto cleanup_and_return;
 
 	/*Fetch data: */
-	ModelFetchData((void**)&spcvelocity,NULL,NULL,model_handle,"spcvelocity","Matrix","Mat");
-	ModelFetchData((void**)&gridonhutter,NULL,NULL,model_handle,"gridonhutter","Matrix","Mat");
+	IoModelFetchData((void**)&spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity","Matrix","Mat");
+	IoModelFetchData((void**)&gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter","Matrix","Mat");
 
 	count=0;
 
 	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
-	for (i=0;i<model->numberofnodes;i++){
+	for (i=0;i<iomodel->numberofnodes;i++){
 	#ifdef _PARALLEL_
 	/*keep only this partition's nodes:*/
-	if((model->my_grids[i]==1)){
+	if((iomodel->my_grids[i]==1)){
 	#endif
 
@@ -66,5 +66,5 @@
 			spc_node=i+1;
 			spc_dof=1; //we enforce first x translation degree of freedom
-			spc_value=*(spcvelocity+6*i+3)/model->yts;
+			spc_value=*(spcvelocity+6*i+3)/iomodel->yts;
 
 			spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
@@ -80,5 +80,5 @@
 			spc_node=i+1;
 			spc_dof=2; //we enforce first y translation degree of freedom
-			spc_value=*(spcvelocity+6*i+4)/model->yts;
+			spc_value=*(spcvelocity+6*i+4)/iomodel->yts;
 			
 			spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
@@ -95,19 +95,19 @@
 
 	/*First fetch penalties: */
-	if (strcmp(model->meshtype,"3d")==0){
-		ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
+	if (strcmp(iomodel->meshtype,"3d")==0){
+		IoModelFetchData((void**)&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties","Matrix","Mat");
 	}
 
 	count=0;
-	if(model->numpenalties){
+	if(iomodel->numpenalties){
 
 		/*First deal with internal grids: */
-		for (i=0;i<model->numpenalties;i++){
-			if (model->penaltypartitioning[i]>=0){ //this penalty belongs to at least this CPU
+		for (i=0;i<iomodel->numpenalties;i++){
+			if (iomodel->penaltypartitioning[i]>=0){ //this penalty belongs to at least this CPU
 
-				for(j=1;j<model->numlayers;j++){
+				for(j=1;j<iomodel->numlayers;j++){
 					/*We are pairing grids along a vertical profile.*/
-					grid1=(int)*(model->penalties+model->numlayers*i);
-					grid2=(int)*(model->penalties+model->numlayers*i+j);
+					grid1=(int)*(iomodel->penalties+iomodel->numlayers*i);
+					grid2=(int)*(iomodel->penalties+iomodel->numlayers*i+j);
 
 					rgb_id=count+1; //matlab indexing
@@ -132,5 +132,5 @@
 
 	/*Free data: */
-	xfree((void**)&model->penalties);
+	xfree((void**)&iomodel->penalties);
 
 	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 1831)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 1832)
@@ -12,8 +12,8 @@
 #include "../../shared/shared.h"
 #include "../../MeshPartitionx/MeshPartitionx.h"
-#include "../Model.h"
-
-
-void	CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){
+#include "../IoModel.h"
+
+
+void	CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
 
@@ -149,8 +149,8 @@
 
 	/*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
-	if (!model->ismacayealpattyn)goto cleanup_and_return;
+	if (!iomodel->ismacayealpattyn)goto cleanup_and_return;
 
 	/*Width of elements: */
-	if(strcmp(model->meshtype,"2d")==0){
+	if(strcmp(iomodel->meshtype,"2d")==0){
 		elements_width=3; //tria elements
 	}
@@ -163,33 +163,33 @@
 	#ifdef _PARALLEL_
 	/*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
-	if(strcmp(model->meshtype,"2d")==0){
+	if(strcmp(iomodel->meshtype,"2d")==0){
 		/*load elements: */
-		ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->elements,NULL,NULL,iomodel_handle,"elements","Matrix","Mat");
 	}
 	else{
 		/*load elements2d: */
-		ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
-	}
-
-	MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs);
+		IoModelFetchData((void**)&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d","Matrix","Mat");
+	}
+
+	MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs);
 
 	/*Free elements and elements2d: */
-	xfree((void**)&model->elements);
-	xfree((void**)&model->elements2d);
+	xfree((void**)&iomodel->elements);
+	xfree((void**)&iomodel->elements2d);
 
 
 	/*Deal with rifts, they have to be included into one partition only, not several: */
-	ModelFetchData((void**)&model->riftinfo,&model->numrifts,NULL,model_handle,"riftinfo","Matrix","Mat");
-
-	for(i=0;i<model->numrifts;i++){
-		el1=(int)*(model->riftinfo+RIFTINFOSIZE*i+2)-1; //matlab indexing to c indexing
-		el2=(int)*(model->riftinfo+RIFTINFOSIZE*i+3)-1; //matlab indexing to c indexing
+	IoModelFetchData((void**)&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo","Matrix","Mat");
+
+	for(i=0;i<iomodel->numrifts;i++){
+		el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2)-1; //matlab indexing to c indexing
+		el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+3)-1; //matlab indexing to c indexing
 		epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids;
 	}
 	/*Free rifts: */
-	xfree((void**)&model->riftinfo); 
+	xfree((void**)&iomodel->riftinfo); 
 
 	/*Used later on: */
-	my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int));
+	my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int));
 	#endif
 
@@ -199,23 +199,23 @@
 
 	/*2d mesh: */
-	if (strcmp(model->meshtype,"2d")==0){
+	if (strcmp(iomodel->meshtype,"2d")==0){
 
 		/*Fetch data needed: */
-		ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
-		ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
-		ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
-		ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
-		ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat");
-		ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat");
-		ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat");
-		ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
-		ModelFetchData((void**)&model->elementonwater,NULL,NULL,model_handle,"elementonwater","Matrix","Mat");
-		ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
-		ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
-		ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
-		ModelFetchData((void**)&model->accumulation,NULL,NULL,model_handle,"accumulation","Matrix","Mat");
-		ModelFetchData((void**)&model->melting,NULL,NULL,model_handle,"melting","Matrix","Mat");
-		
-		for (i=0;i<model->numberofelements;i++){
+		IoModelFetchData((void**)&iomodel->elements,NULL,NULL,iomodel_handle,"elements","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->surface,NULL,NULL,iomodel_handle,"surface","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->bed,NULL,NULL,iomodel_handle,"bed","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->drag,NULL,NULL,iomodel_handle,"drag","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->p,NULL,NULL,iomodel_handle,"p","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->q,NULL,NULL,iomodel_handle,"q","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->B,NULL,NULL,iomodel_handle,"B","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->n,NULL,NULL,iomodel_handle,"n","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->accumulation,NULL,NULL,iomodel_handle,"accumulation","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->melting,NULL,NULL,iomodel_handle,"melting","Matrix","Mat");
+		
+		for (i=0;i<iomodel->numberofelements;i++){
 
 		#ifdef _PARALLEL_
@@ -224,57 +224,57 @@
 		#endif
 			
-			if (*(model->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
+			if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
 				
 				/*ids: */
 				tria_id=i+1; //matlab indexing.
 				tria_mid=i+1; //refers to the corresponding material property card
-				tria_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
+				tria_mparid=iomodel->numberofelements+1;//refers to the corresponding parmat property card
 
 				/*vertices offsets: */
-				tria_g[0]=(int)*(model->elements+elements_width*i+0);
-				tria_g[1]=(int)*(model->elements+elements_width*i+1);
-				tria_g[2]=(int)*(model->elements+elements_width*i+2);
+				tria_g[0]=(int)*(iomodel->elements+elements_width*i+0);
+				tria_g[1]=(int)*(iomodel->elements+elements_width*i+1);
+				tria_g[2]=(int)*(iomodel->elements+elements_width*i+2);
 
 				/*thickness,surface and bed:*/
-				tria_h[0]= *(model->thickness+ ((int)*(model->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
-				tria_h[1]=*(model->thickness+  ((int)*(model->elements+elements_width*i+1)-1)); 
-				tria_h[2]=*(model->thickness+  ((int)*(model->elements+elements_width*i+2)-1)) ;
-
-				tria_s[0]=*(model->surface+    ((int)*(model->elements+elements_width*i+0)-1)); 
-				tria_s[1]=*(model->surface+    ((int)*(model->elements+elements_width*i+1)-1)); 
-				tria_s[2]=*(model->surface+    ((int)*(model->elements+elements_width*i+2)-1)); 
-
-				tria_b[0]=*(model->bed+        ((int)*(model->elements+elements_width*i+0)-1)); 
-				tria_b[1]=*(model->bed+        ((int)*(model->elements+elements_width*i+1)-1)); 
-				tria_b[2]=*(model->bed+        ((int)*(model->elements+elements_width*i+2)-1)); 
+				tria_h[0]= *(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
+				tria_h[1]=*(iomodel->thickness+  ((int)*(iomodel->elements+elements_width*i+1)-1)); 
+				tria_h[2]=*(iomodel->thickness+  ((int)*(iomodel->elements+elements_width*i+2)-1)) ;
+
+				tria_s[0]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+0)-1)); 
+				tria_s[1]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+1)-1)); 
+				tria_s[2]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+2)-1)); 
+
+				tria_b[0]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+0)-1)); 
+				tria_b[1]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+1)-1)); 
+				tria_b[2]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+2)-1)); 
 
 				/*basal drag:*/
-				tria_friction_type=(int)model->drag_type;
-
-				tria_k[0]=*(model->drag+        ((int)*(model->elements+elements_width*i+0)-1)); 
-				tria_k[1]=*(model->drag+        ((int)*(model->elements+elements_width*i+1)-1)); 
-				tria_k[2]=*(model->drag+        ((int)*(model->elements+elements_width*i+2)-1)); 
+				tria_friction_type=(int)iomodel->drag_type;
+
+				tria_k[0]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+0)-1)); 
+				tria_k[1]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+1)-1)); 
+				tria_k[2]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+2)-1)); 
 				
-				tria_p=model->p[i];
-				tria_q=model->q[i];
+				tria_p=iomodel->p[i];
+				tria_q=iomodel->q[i];
 
 				/*meling and accumulation*/
-				tria_melting[0]=*(model->melting+        ((int)*(model->elements+elements_width*i+0)-1));
-				tria_melting[1]=*(model->melting+        ((int)*(model->elements+elements_width*i+1)-1));
-				tria_melting[2]=*(model->melting+        ((int)*(model->elements+elements_width*i+2)-1));
-
-				tria_accumulation[0]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+0)-1));
-				tria_accumulation[1]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+1)-1));
-				tria_accumulation[2]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+2)-1));
+				tria_melting[0]=*(iomodel->melting+        ((int)*(iomodel->elements+elements_width*i+0)-1));
+				tria_melting[1]=*(iomodel->melting+        ((int)*(iomodel->elements+elements_width*i+1)-1));
+				tria_melting[2]=*(iomodel->melting+        ((int)*(iomodel->elements+elements_width*i+2)-1));
+
+				tria_accumulation[0]=*(iomodel->accumulation+        ((int)*(iomodel->elements+elements_width*i+0)-1));
+				tria_accumulation[1]=*(iomodel->accumulation+        ((int)*(iomodel->elements+elements_width*i+1)-1));
+				tria_accumulation[2]=*(iomodel->accumulation+        ((int)*(iomodel->elements+elements_width*i+2)-1));
 
 				/*element on iceshelf, water?:*/
-				tria_shelf=(int)*(model->elementoniceshelf+i);
-				tria_onwater=(bool)*(model->elementonwater+i);
-
-				tria_meanvel=model->meanvel;
-				tria_epsvel=model->epsvel;
+				tria_shelf=(int)*(iomodel->elementoniceshelf+i);
+				tria_onwater=(bool)*(iomodel->elementonwater+i);
+
+				tria_meanvel=iomodel->meanvel;
+				tria_epsvel=iomodel->epsvel;
 
 				/*viscosity_overshoot*/
-				tria_viscosity_overshoot=model->viscosity_overshoot;
+				tria_viscosity_overshoot=iomodel->viscosity_overshoot;
 
 				/*Create tria element using its constructor:*/
@@ -290,9 +290,9 @@
 				B_avg=0;
 				for(j=0;j<3;j++){
-					B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
+					B_avg+=*(iomodel->B+((int)*(iomodel->elements+elements_width*i+j)-1));
 				}
 				B_avg=B_avg/3;
 				matice_B=B_avg;
-				matice_n=(double)*(model->n+i);
+				matice_n=(double)*(iomodel->n+i);
 				
 				/*Create matice using its constructor:*/
@@ -309,7 +309,7 @@
 			 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
 			 will hold which grids belong to this partition*/
-			my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
-			my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
-			my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
 			#endif
 
@@ -322,16 +322,16 @@
 	
 		/*Free data : */
-		xfree((void**)&model->elements);
-		xfree((void**)&model->thickness);
-		xfree((void**)&model->surface);
-		xfree((void**)&model->bed);
-		xfree((void**)&model->drag);
-		xfree((void**)&model->p);
-		xfree((void**)&model->q);
-		xfree((void**)&model->elementoniceshelf);
-		xfree((void**)&model->elementonwater);
-		xfree((void**)&model->B);
-		xfree((void**)&model->n);
-		xfree((void**)&model->elements_type);
+		xfree((void**)&iomodel->elements);
+		xfree((void**)&iomodel->thickness);
+		xfree((void**)&iomodel->surface);
+		xfree((void**)&iomodel->bed);
+		xfree((void**)&iomodel->drag);
+		xfree((void**)&iomodel->p);
+		xfree((void**)&iomodel->q);
+		xfree((void**)&iomodel->elementoniceshelf);
+		xfree((void**)&iomodel->elementonwater);
+		xfree((void**)&iomodel->B);
+		xfree((void**)&iomodel->n);
+		xfree((void**)&iomodel->elements_type);
 
 	}
@@ -339,22 +339,22 @@
 
 		/*Fetch data needed: */
-		ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
-		ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
-		ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
-		ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
-		ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat");
-		ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat");
-		ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat");
-		ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
-		ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat");
-		ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat");
-		ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
-		ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
-		ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
-		ModelFetchData((void**)&model->accumulation,NULL,NULL,model_handle,"accumulation","Matrix","Mat");
-		ModelFetchData((void**)&model->melting,NULL,NULL,model_handle,"melting","Matrix","Mat");
-		ModelFetchData((void**)&model->elementonwater,NULL,NULL,model_handle,"elementonwater","Matrix","Mat");
-		
-		for (i=0;i<model->numberofelements;i++){
+		IoModelFetchData((void**)&iomodel->elements,NULL,NULL,iomodel_handle,"elements","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->surface,NULL,NULL,iomodel_handle,"surface","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->bed,NULL,NULL,iomodel_handle,"bed","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->drag,NULL,NULL,iomodel_handle,"drag","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->p,NULL,NULL,iomodel_handle,"p","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->q,NULL,NULL,iomodel_handle,"q","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->B,NULL,NULL,iomodel_handle,"B","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->n,NULL,NULL,iomodel_handle,"n","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->accumulation,NULL,NULL,iomodel_handle,"accumulation","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->melting,NULL,NULL,iomodel_handle,"melting","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater","Matrix","Mat");
+		
+		for (i=0;i<iomodel->numberofelements;i++){
 		#ifdef _PARALLEL_
 		/*We are using our element partition to decide which elements will be created on this node: */
@@ -362,40 +362,40 @@
 		#endif
 
-			if (*(model->elements_type+2*i+0)==MacAyealFormulationEnum() | *(model->elements_type+2*i+0)==PattynFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
+			if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum() | *(iomodel->elements_type+2*i+0)==PattynFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements.
 			
 				/*name and id: */
 				penta_id=i+1; //matlab indexing.
 				penta_mid=i+1; //refers to the corresponding material property card
-				penta_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
+				penta_mparid=iomodel->numberofelements+1;//refers to the corresponding parmat property card
 
 				/*vertices,thickness,surface,bed and drag: */
 				for(j=0;j<6;j++){
-					penta_g[j]=(int)*(model->elements+elements_width*i+j);
-					penta_h[j]=*(model->thickness+    ((int)*(model->elements+elements_width*i+j)-1)); 
-					penta_s[j]=*(model->surface+    ((int)*(model->elements+elements_width*i+j)-1)); 
-					penta_b[j]=*(model->bed+    ((int)*(model->elements+elements_width*i+j)-1)); 
-					penta_k[j]=*(model->drag+        ((int)*(model->elements+elements_width*i+j)-1)); 
-					penta_melting[j]=*(model->melting+        ((int)*(model->elements+elements_width*i+j)-1));
-					penta_accumulation[j]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+j)-1));
+					penta_g[j]=(int)*(iomodel->elements+elements_width*i+j);
+					penta_h[j]=*(iomodel->thickness+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
+					penta_s[j]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
+					penta_b[j]=*(iomodel->bed+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
+					penta_k[j]=*(iomodel->drag+        ((int)*(iomodel->elements+elements_width*i+j)-1)); 
+					penta_melting[j]=*(iomodel->melting+        ((int)*(iomodel->elements+elements_width*i+j)-1));
+					penta_accumulation[j]=*(iomodel->accumulation+        ((int)*(iomodel->elements+elements_width*i+j)-1));
 				}
 
 				/*basal drag:*/
-				penta_friction_type=(int)model->drag_type;
-		
-				penta_p=model->p[i];
-				penta_q=model->q[i];
+				penta_friction_type=(int)iomodel->drag_type;
+		
+				penta_p=iomodel->p[i];
+				penta_q=iomodel->q[i];
 
 				/*diverse: */
-				penta_shelf=(int)*(model->elementoniceshelf+i);
-				penta_onbed=(int)*(model->elementonbed+i);
-				penta_onsurface=(int)*(model->elementonsurface+i);
-				penta_meanvel=model->meanvel;
-				penta_epsvel=model->epsvel;
-				penta_onwater=(bool)*(model->elementonwater+i);
+				penta_shelf=(int)*(iomodel->elementoniceshelf+i);
+				penta_onbed=(int)*(iomodel->elementonbed+i);
+				penta_onsurface=(int)*(iomodel->elementonsurface+i);
+				penta_meanvel=iomodel->meanvel;
+				penta_epsvel=iomodel->epsvel;
+				penta_onwater=(bool)*(iomodel->elementonwater+i);
 				
 				/*viscosity_overshoot*/
-				penta_viscosity_overshoot=model->viscosity_overshoot;
-
-				if (*(model->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base.
+				penta_viscosity_overshoot=iomodel->viscosity_overshoot;
+
+				if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base.
 					penta_collapse=1;
 				}
@@ -420,9 +420,9 @@
 				B_avg=0;
 				for(j=0;j<6;j++){
-					B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
+					B_avg+=*(iomodel->B+((int)*(iomodel->elements+elements_width*i+j)-1));
 				}
 				B_avg=B_avg/6;
 				matice_B= B_avg;
-				matice_n=(double)*(model->n+i);
+				matice_n=(double)*(iomodel->n+i);
 		
 				/*Create matice using its constructor:*/
@@ -439,10 +439,10 @@
 			 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
 			 will hold which grids belong to this partition*/
-			my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
-			my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
-			my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
-			my_grids[(int)*(model->elements+elements_width*i+3)-1]=1;
-			my_grids[(int)*(model->elements+elements_width*i+4)-1]=1;
-			my_grids[(int)*(model->elements+elements_width*i+5)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;
+			my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;
 			#endif
 
@@ -454,35 +454,35 @@
 
 		/*Free data: */
-		xfree((void**)&model->elements);
-		xfree((void**)&model->thickness);
-		xfree((void**)&model->surface);
-		xfree((void**)&model->bed);
-		xfree((void**)&model->drag);
-		xfree((void**)&model->p);
-		xfree((void**)&model->q);
-		xfree((void**)&model->elementoniceshelf);
-		xfree((void**)&model->elementonbed);
-		xfree((void**)&model->elementonsurface);
-		xfree((void**)&model->elements_type);
-		xfree((void**)&model->n);
-		xfree((void**)&model->B);
-		xfree((void**)&model->accumulation);
-		xfree((void**)&model->melting);
-		xfree((void**)&model->elementonwater);
+		xfree((void**)&iomodel->elements);
+		xfree((void**)&iomodel->thickness);
+		xfree((void**)&iomodel->surface);
+		xfree((void**)&iomodel->bed);
+		xfree((void**)&iomodel->drag);
+		xfree((void**)&iomodel->p);
+		xfree((void**)&iomodel->q);
+		xfree((void**)&iomodel->elementoniceshelf);
+		xfree((void**)&iomodel->elementonbed);
+		xfree((void**)&iomodel->elementonsurface);
+		xfree((void**)&iomodel->elements_type);
+		xfree((void**)&iomodel->n);
+		xfree((void**)&iomodel->B);
+		xfree((void**)&iomodel->accumulation);
+		xfree((void**)&iomodel->melting);
+		xfree((void**)&iomodel->elementonwater);
 
 	} //if (strcmp(meshtype,"2d")==0)
 
 	/*Add one constant material property to materials: */
-	matpar_mid=model->numberofelements+1; //put it at the end of the materials
-	matpar_g=model->g; 
-	matpar_rho_ice=model->rho_ice; 
-	matpar_rho_water=model->rho_water; 
-	matpar_thermalconductivity=model->thermalconductivity; 
-	matpar_heatcapacity=model->heatcapacity; 
-	matpar_latentheat=model->latentheat; 
-	matpar_beta=model->beta; 
-	matpar_meltingpoint=model->meltingpoint; 
-	matpar_mixed_layer_capacity=model->mixed_layer_capacity; 
-	matpar_thermal_exchange_velocity=model->thermal_exchange_velocity; 
+	matpar_mid=iomodel->numberofelements+1; //put it at the end of the materials
+	matpar_g=iomodel->g; 
+	matpar_rho_ice=iomodel->rho_ice; 
+	matpar_rho_water=iomodel->rho_water; 
+	matpar_thermalconductivity=iomodel->thermalconductivity; 
+	matpar_heatcapacity=iomodel->heatcapacity; 
+	matpar_latentheat=iomodel->latentheat; 
+	matpar_beta=iomodel->beta; 
+	matpar_meltingpoint=iomodel->meltingpoint; 
+	matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; 
+	matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; 
 
 	/*Create matpar object using its constructor: */
@@ -497,7 +497,7 @@
 		/*From the element partitioning, we can determine which grids are on the inside of this cpu's 
 		 *element partition, and which are on its border with other nodes:*/
-		gridborder=NewVec(model->numberofnodes);
-
-		for (i=0;i<model->numberofnodes;i++){
+		gridborder=NewVec(iomodel->numberofnodes);
+
+		for (i=0;i<iomodel->numberofnodes;i++){
 			if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
 		}
@@ -513,5 +513,5 @@
 		#ifdef _DEBUG_
 		if(my_rank==0){
-			for (i=0;i<model->numberofnodes;i++){
+			for (i=0;i<iomodel->numberofnodes;i++){
 				printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
 			}
@@ -521,25 +521,25 @@
 
 	/*Partition penalties in 3d: */
-	if(strcmp(model->meshtype,"3d")==0){
+	if(strcmp(iomodel->meshtype,"3d")==0){
 	
 		/*Get penalties: */
-		ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
-
-		if(model->numpenalties){
-
-			model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int));
+		IoModelFetchData((void**)&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties","Matrix","Mat");
+
+		if(iomodel->numpenalties){
+
+			iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int));
 			#ifdef _SERIAL_
-			for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=1;
+			for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1;
 			#else
-			for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=-1;
-
-			for(i=0;i<model->numpenalties;i++){
-				first_grid_index=(int)(*(model->penalties+i*model->numlayers+0)-1);
+			for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1;
+
+			for(i=0;i<iomodel->numpenalties;i++){
+				first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1);
 				if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
 					/*All grids that are being penalised belong to this node's internal grid partition.:*/
-					model->penaltypartitioning[i]=1;
+					iomodel->penaltypartitioning[i]=1;
 				}
 				if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
-					model->penaltypartitioning[i]=0;
+					iomodel->penaltypartitioning[i]=0;
 				}
 			}
@@ -548,5 +548,5 @@
 
 		/*Free penalties: */
-		xfree((void**)&model->penalties);
+		xfree((void**)&iomodel->penalties);
 	}
 
@@ -559,23 +559,23 @@
 		
 	/*First fetch data: */
-	if (strcmp(model->meshtype,"3d")==0){
-		ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat");
-		ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat");
-	}
-	ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat");
-	ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat");
-	ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat");
-	ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
-	ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
-	ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
-	ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
-	ModelFetchData((void**)&model->gridonhutter,NULL,NULL,model_handle,"gridonhutter","Matrix","Mat");
-	ModelFetchData((void**)&model->gridonicesheet,NULL,NULL,model_handle,"gridonicesheet","Matrix","Mat");
-	ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
+	if (strcmp(iomodel->meshtype,"3d")==0){
+		IoModelFetchData((void**)&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids","Matrix","Mat");
+	}
+	IoModelFetchData((void**)&iomodel->x,NULL,NULL,iomodel_handle,"x","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->y,NULL,NULL,iomodel_handle,"y","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->z,NULL,NULL,iomodel_handle,"z","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->bed,NULL,NULL,iomodel_handle,"bed","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf","Matrix","Mat");
 
 	/*Get number of dofs per node: */
-	DistributeNumDofs(&node_numdofs,model->analysis_type,model->sub_analysis_type);
-
-	for (i=0;i<model->numberofnodes;i++){
+	DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
+
+	for (i=0;i<iomodel->numberofnodes;i++){
 	#ifdef _PARALLEL_
 
@@ -598,20 +598,20 @@
 		#endif
 
-		node_x[0]=model->x[i];
-		node_x[1]=model->y[i];
-		node_x[2]=model->z[i];
-		node_sigma=(model->z[i]-model->bed[i])/(model->thickness[i]);
-		
-		node_onbed=(int)model->gridonbed[i];
-		node_onsurface=(int)model->gridonsurface[i];	
-		node_onshelf=(int)model->gridoniceshelf[i];	
-		node_onsheet=(int)model->gridonicesheet[i];	
-
-		if (strcmp(model->meshtype,"3d")==0){
-			if (isnan(model->uppernodes[i])){
+		node_x[0]=iomodel->x[i];
+		node_x[1]=iomodel->y[i];
+		node_x[2]=iomodel->z[i];
+		node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
+		
+		node_onbed=(int)iomodel->gridonbed[i];
+		node_onsurface=(int)iomodel->gridonsurface[i];	
+		node_onshelf=(int)iomodel->gridoniceshelf[i];	
+		node_onsheet=(int)iomodel->gridonicesheet[i];	
+
+		if (strcmp(iomodel->meshtype,"3d")==0){
+			if (isnan(iomodel->uppernodes[i])){
 				node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
 			}
 			else{
-				node_upper_node_id=(int)model->uppernodes[i];
+				node_upper_node_id=(int)iomodel->uppernodes[i];
 			}
 		}
@@ -625,7 +625,7 @@
 
 		/*set single point constraints.: */
-		if (strcmp(model->meshtype,"3d")==0){
+		if (strcmp(iomodel->meshtype,"3d")==0){
 			/*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
-			if (model->deadgrids[i]){
+			if (iomodel->deadgrids[i]){
 				for(k=1;k<=node_numdofs;k++){
 					node->FreezeDof(k);
@@ -633,5 +633,5 @@
 			}
 		}
-		if (model->gridonhutter[i]){
+		if (iomodel->gridonhutter[i]){
 			for(k=1;k<=node_numdofs;k++){
 				node->FreezeDof(k);
@@ -653,22 +653,22 @@
 
 	/*Clean fetched data: */
-	xfree((void**)&model->deadgrids);
-	xfree((void**)&model->x);
-	xfree((void**)&model->y);
-	xfree((void**)&model->z);
-	xfree((void**)&model->thickness);
-	xfree((void**)&model->bed);
-	xfree((void**)&model->gridonbed);
-	xfree((void**)&model->gridonsurface);
-	xfree((void**)&model->gridonhutter);
-	xfree((void**)&model->uppernodes);
-	xfree((void**)&model->gridonicesheet);
-	xfree((void**)&model->gridoniceshelf);
+	xfree((void**)&iomodel->deadgrids);
+	xfree((void**)&iomodel->x);
+	xfree((void**)&iomodel->y);
+	xfree((void**)&iomodel->z);
+	xfree((void**)&iomodel->thickness);
+	xfree((void**)&iomodel->bed);
+	xfree((void**)&iomodel->gridonbed);
+	xfree((void**)&iomodel->gridonsurface);
+	xfree((void**)&iomodel->gridonhutter);
+	xfree((void**)&iomodel->uppernodes);
+	xfree((void**)&iomodel->gridonicesheet);
+	xfree((void**)&iomodel->gridoniceshelf);
 	
 
-	/*Keep partitioning information into model*/
-	model->epart=epart;
-	model->my_grids=my_grids;
-	model->my_bordergrids=my_bordergrids;
+	/*Keep partitioning information into iomodel*/
+	iomodel->epart=epart;
+	iomodel->my_grids=my_grids;
+	iomodel->my_bordergrids=my_bordergrids;
 
 	/*Free ressources:*/
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp	(revision 1831)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp	(revision 1832)
@@ -11,8 +11,8 @@
 #include "../../shared/shared.h"
 #include "../../include/macros.h"
-#include "../Model.h"
-
-
-void	CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model,ConstDataHandle model_handle){
+#include "../IoModel.h"
+
+
+void	CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
 	int i;
@@ -55,5 +55,4 @@
 	double riftfront_friction;
 	double riftfront_fraction;
-	double riftfront_fractionincrement;
 	bool   riftfront_shelf;
 	double riftfront_penalty_offset;
@@ -69,5 +68,4 @@
 	double friction;
 	double fraction;
-	double fractionincrement;
 	
 	/*Create loads: */
@@ -75,17 +73,17 @@
 
 	/*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
-	if (!model->ismacayealpattyn)goto cleanup_and_return;
+	if (!iomodel->ismacayealpattyn)goto cleanup_and_return;
 	
 	/*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): */
-	ModelFetchData((void**)&model->pressureload,&model->numberofpressureloads,NULL,model_handle,"pressureload","Matrix","Mat");
-	ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
-	ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
-	ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->pressureload,&iomodel->numberofpressureloads,NULL,iomodel_handle,"pressureload","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->bed,NULL,NULL,iomodel_handle,"bed","Matrix","Mat");
 
 	/*First load data:*/
-	for (i=0;i<model->numberofpressureloads;i++){
-		
-		if (strcmp(model->meshtype,"2d")==0){
+	for (i=0;i<iomodel->numberofpressureloads;i++){
+		
+		if (strcmp(iomodel->meshtype,"2d")==0){
 			segment_width=3;
 			element_type=TriaEnum();
@@ -97,8 +95,8 @@
 
 
-		element=(int)(*(model->pressureload+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_
-		if (model->epart[element]!=my_rank){
+		if (iomodel->epart[element]!=my_rank){
 			/*This load does not belong to this cluster node, as it references an element 
 			 *that does not belong to this node's partition. Drop this 'i':*/
@@ -107,39 +105,39 @@
 		#endif
 	
-		icefront_mparid=model->numberofelements+1; //matlab indexing
+		icefront_mparid=iomodel->numberofelements+1; //matlab indexing
 		icefront_sid=i+1; //matlab indexing
-		icefront_eid=(int)*(model->pressureload+segment_width*i+segment_width-1); //matlab indexing
+		icefront_eid=(int)*(iomodel->pressureload+segment_width*i+segment_width-1); //matlab indexing
 		icefront_element_type=element_type;
 
-		i1=(int)*(model->pressureload+segment_width*i+0);
-		i2=(int)*(model->pressureload+segment_width*i+1);
+		i1=(int)*(iomodel->pressureload+segment_width*i+0);
+		i2=(int)*(iomodel->pressureload+segment_width*i+1);
 			
 		icefront_node_ids[0]=i1;
 		icefront_node_ids[1]=i2;
 		
-		icefront_h[0]=model->thickness[i1-1];
-		icefront_h[1]=model->thickness[i2-1];
-
-		icefront_b[0]=model->bed[i1-1];
-		icefront_b[1]=model->bed[i2-1];
-		
-		if ((int)*(model->elements_type+2*element+0)==MacAyealFormulationEnum()){ //this is a collapsed 3d element, icefront will be 2d
+		icefront_h[0]=iomodel->thickness[i1-1];
+		icefront_h[1]=iomodel->thickness[i2-1];
+
+		icefront_b[0]=iomodel->bed[i1-1];
+		icefront_b[1]=iomodel->bed[i2-1];
+		
+		if ((int)*(iomodel->elements_type+2*element+0)==MacAyealFormulationEnum()){ //this is a collapsed 3d element, icefront will be 2d
 			strcpy(icefront_type,"segment");
 		}
-		else if ((int)*(model->elements_type+2*element+0)==PattynFormulationEnum()){ //this is a real 3d element, icefront will be 3d.
+		else if ((int)*(iomodel->elements_type+2*element+0)==PattynFormulationEnum()){ //this is a real 3d element, icefront will be 3d.
 			strcpy(icefront_type,"quad");
-			i3=(int)*(model->pressureload+segment_width*i+2);
-			i4=(int)*(model->pressureload+segment_width*i+3);
+			i3=(int)*(iomodel->pressureload+segment_width*i+2);
+			i4=(int)*(iomodel->pressureload+segment_width*i+3);
 			icefront_node_ids[2]=i3;
 			icefront_node_ids[3]=i4;
 			
-			icefront_h[2]=model->thickness[i3-1];
-			icefront_h[3]=model->thickness[i4-1];
-
-			icefront_b[2]=model->bed[i3-1];
-			icefront_b[3]=model->bed[i4-1];
+			icefront_h[2]=iomodel->thickness[i3-1];
+			icefront_h[3]=iomodel->thickness[i4-1];
+
+			icefront_b[2]=iomodel->bed[i3-1];
+			icefront_b[3]=iomodel->bed[i4-1];
 		}
 		else{
-			throw ErrorException(__FUNCT__,exprintf(" element type %i not supported yet",(int)*(model->elements_type+2*element+0)));
+			throw ErrorException(__FUNCT__,exprintf(" element type %i not supported yet",(int)*(iomodel->elements_type+2*element+0)));
 		}
 
@@ -150,23 +148,23 @@
 	}
 	/*Free data: */
-	xfree((void**)&model->pressureload);
-	xfree((void**)&model->elements_type);
-	xfree((void**)&model->thickness);
-	xfree((void**)&model->bed);
+	xfree((void**)&iomodel->pressureload);
+	xfree((void**)&iomodel->elements_type);
+	xfree((void**)&iomodel->thickness);
+	xfree((void**)&iomodel->bed);
 
 	/*Create penpair loads also for rift grids, so that they don't penetrate one another, if needed: */
 	/*First fetch rifts: */
-	ModelFetchData((void**)&model->riftinfo,&model->numrifts,NULL,model_handle,"riftinfo","Matrix","Mat");
-	ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
-	ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
-	ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
-	ModelFetchData((void**)&model->gridoniceshelf,NULL,NULL,model_handle,"gridoniceshelf","Matrix","Mat");
-	
-	if(model->numrifts){
-		for(i=0;i<model->numrifts;i++){
+	IoModelFetchData((void**)&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->bed,NULL,NULL,iomodel_handle,"bed","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->surface,NULL,NULL,iomodel_handle,"surface","Matrix","Mat");
+	IoModelFetchData((void**)&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf","Matrix","Mat");
+	
+	if(iomodel->numrifts){
+		for(i=0;i<iomodel->numrifts;i++){
 				
-			el1=(int)*(model->riftinfo+RIFTINFOSIZE*i+2);
+			el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2);
 			#ifdef _PARALLEL_
-			if (model->epart[el1-1]!=my_rank){
+			if (iomodel->epart[el1-1]!=my_rank){
 				/*This load does not belong to this cluster node, as it references an element 
 				 *that does not belong to this node's partition. Drop this 'j':*/
@@ -176,16 +174,16 @@
 
 			/*Ok, retrieve all the data needed to add a penalty between the two grids: */
-			el2=(int)*(model->riftinfo+RIFTINFOSIZE*i+3); 
-
-			grid1=(int)*(model->riftinfo+RIFTINFOSIZE*i+0); 
-			grid2=(int)*(model->riftinfo+RIFTINFOSIZE*i+1);
-			
-			normal[0]=*(model->riftinfo+RIFTINFOSIZE*i+4);
-			normal[1]=*(model->riftinfo+RIFTINFOSIZE*i+5);
-			length=*(model->riftinfo+RIFTINFOSIZE*i+6);
-			
-			fill = (int)*(model->riftinfo+RIFTINFOSIZE*i+7);
-			friction=*(model->riftinfo+RIFTINFOSIZE*i+8);
-			fraction=*(model->riftinfo+RIFTINFOSIZE*i+9);
+			el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+3); 
+
+			grid1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+0); 
+			grid2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+1);
+			
+			normal[0]=*(iomodel->riftinfo+RIFTINFOSIZE*i+4);
+			normal[1]=*(iomodel->riftinfo+RIFTINFOSIZE*i+5);
+			length=*(iomodel->riftinfo+RIFTINFOSIZE*i+6);
+			
+			fill = (int)*(iomodel->riftinfo+RIFTINFOSIZE*i+7);
+			friction=*(iomodel->riftinfo+RIFTINFOSIZE*i+8);
+			fraction=*(iomodel->riftinfo+RIFTINFOSIZE*i+9);
 			fractionincrement=*(model->riftinfo+RIFTINFOSIZE*i+10);
 	
@@ -194,14 +192,14 @@
 			riftfront_node_ids[0]=grid1;
 			riftfront_node_ids[1]=grid2;
-			riftfront_mparid=model->numberofelements+1; //matlab indexing
-
-			riftfront_h[0]=model->thickness[grid1-1];
-			riftfront_h[1]=model->thickness[grid2-1];
-
-			riftfront_b[0]=model->bed[grid1-1];
-			riftfront_b[1]=model->bed[grid2-1];
-
-			riftfront_s[0]=model->surface[grid1-1];
-			riftfront_s[1]=model->surface[grid2-1];
+			riftfront_mparid=iomodel->numberofelements+1; //matlab indexing
+
+			riftfront_h[0]=iomodel->thickness[grid1-1];
+			riftfront_h[1]=iomodel->thickness[grid2-1];
+
+			riftfront_b[0]=iomodel->bed[grid1-1];
+			riftfront_b[1]=iomodel->bed[grid2-1];
+
+			riftfront_s[0]=iomodel->surface[grid1-1];
+			riftfront_s[1]=iomodel->surface[grid2-1];
 
 			riftfront_normal[0]=normal[0];
@@ -215,6 +213,6 @@
 			riftfront_shelf=(bool)model->gridoniceshelf[grid1-1];
 
-			riftfront_penalty_offset=model->penalty_offset;
-			riftfront_penalty_lock=model->penalty_lock;
+			riftfront_penalty_offset=iomodel->penalty_offset;
+			riftfront_penalty_lock=iomodel->penalty_lock;
 
 			riftfront_active=0;
@@ -222,5 +220,5 @@
 			riftfront_prestable=0;
 			
-			riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction, riftfront_fractionincrement, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_counter,riftfront_prestable,riftfront_shelf);
+			riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_counter,riftfront_prestable,riftfront_shelf);
 
 			loads->AddObject(riftfront);
@@ -230,9 +228,9 @@
 
 	/*free ressources: */
-	xfree((void**)&model->riftinfo);
-	xfree((void**)&model->thickness);
-	xfree((void**)&model->bed);
-	xfree((void**)&model->surface);
-	xfree((void**)&model->gridoniceshelf);
+	xfree((void**)&iomodel->riftinfo);
+	xfree((void**)&iomodel->thickness);
+	xfree((void**)&iomodel->bed);
+	xfree((void**)&iomodel->surface);
+	xfree((void**)&iomodel->gridoniceshelf);
 
 
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp	(revision 1831)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp	(revision 1832)
@@ -11,7 +11,7 @@
 #include "../../objects/objects.h"
 #include "../../shared/shared.h"
-#include "../Model.h"
+#include "../IoModel.h"
 
-void CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle){
+void CreateParametersDiagnosticHoriz(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
 	
 	Param*   param = NULL;
@@ -31,20 +31,20 @@
 
 	/*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
-//	if (!model->ismacayealpattyn)return;
+//	if (!iomodel->ismacayealpattyn)return;
 
 	/*Get vx and vy: */
-	ModelFetchData((void**)&vx,NULL,NULL,model_handle,"vx","Matrix","Mat");
-	ModelFetchData((void**)&vy,NULL,NULL,model_handle,"vy","Matrix","Mat");
-	ModelFetchData((void**)&vz,NULL,NULL,model_handle,"vz","Matrix","Mat");
+	IoModelFetchData((void**)&vx,NULL,NULL,iomodel_handle,"vx","Matrix","Mat");
+	IoModelFetchData((void**)&vy,NULL,NULL,iomodel_handle,"vy","Matrix","Mat");
+	IoModelFetchData((void**)&vz,NULL,NULL,iomodel_handle,"vz","Matrix","Mat");
 
-	ug=(double*)xcalloc(model->numberofnodes*3,sizeof(double));
+	ug=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double));
 
-	if(vx) for(i=0;i<model->numberofnodes;i++)ug[3*i+0]=vx[i]/model->yts;
-	if(vy) for(i=0;i<model->numberofnodes;i++)ug[3*i+1]=vy[i]/model->yts;
-	if(vz) for(i=0;i<model->numberofnodes;i++)ug[3*i+2]=vz[i]/model->yts;
+	if(vx) for(i=0;i<iomodel->numberofnodes;i++)ug[3*i+0]=vx[i]/iomodel->yts;
+	if(vy) for(i=0;i<iomodel->numberofnodes;i++)ug[3*i+1]=vy[i]/iomodel->yts;
+	if(vz) for(i=0;i<iomodel->numberofnodes;i++)ug[3*i+2]=vz[i]/iomodel->yts;
 
 	count++;
 	param= new Param(count,"u_g",DOUBLEVEC);
-	param->SetDoubleVec(ug,3*model->numberofnodes,3);
+	param->SetDoubleVec(ug,3*iomodel->numberofnodes,3);
 	parameters->AddObject(param);