Index: /issm/trunk/src/c/parallel/control.cpp
===================================================================
--- /issm/trunk/src/c/parallel/control.cpp	(revision 1823)
+++ /issm/trunk/src/c/parallel/control.cpp	(revision 1824)
@@ -27,5 +27,5 @@
 
 	/*Fem models : */
-	FemModel femmodels[5];
+	Model* model=NULL;
 
 	/*Results: */
@@ -62,20 +62,23 @@
 	fid=pfopen(inputfilename,"rb");
 
+	/*Initialize model structure: */
+	model=new Model();
+
 	_printf_("read and create finite element model:\n");
 	_printf_("\n   reading control horiz model data:\n");
-	CreateFemModel(&femmodels[0],fid,ControlAnalysisEnum(),HorizAnalysisEnum());
+	CreateFemModel(model->DiagnosticHorizontal(),fid,ControlAnalysisEnum(),HorizAnalysisEnum());
 	_printf_("\n   reading control vert model data:\n");
-	CreateFemModel(&femmodels[1],fid,ControlAnalysisEnum(),VertAnalysisEnum());
+	CreateFemModel(model->DiagnosticVertical(),fid,ControlAnalysisEnum(),VertAnalysisEnum());
 	_printf_("\n   reading control stokes model data:\n");
-	CreateFemModel(&femmodels[2],fid,ControlAnalysisEnum(),StokesAnalysisEnum());
+	CreateFemModel(model->DiagnosticStokes(),fid,ControlAnalysisEnum(),StokesAnalysisEnum());
 	_printf_("\n   reading control hutter model data:\n");
-	CreateFemModel(&femmodels[3],fid,ControlAnalysisEnum(),HutterAnalysisEnum());
+	CreateFemModel(model->DiagnosticHutter(),fid,ControlAnalysisEnum(),HutterAnalysisEnum());
 	_printf_("\n   reading surface and bed slope computation model data:\n");
-	CreateFemModel(&femmodels[4],fid,SlopeComputeAnalysisEnum(),NoneAnalysisEnum());
+	CreateFemModel(model->Slope(),fid,SlopeComputeAnalysisEnum(),NoneAnalysisEnum());
 
 	_printf_("initialize inputs:\n");
-	femmodels[0].parameters->FindParam((void*)&u_g_initial,"u_g");
-	femmodels[0].parameters->FindParam((void*)&u_g_obs,"u_g_obs");
-	femmodels[0].parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+	model->DiagnosticHorizontal()->parameters->FindParam((void*)&u_g_initial,"u_g");
+	model->DiagnosticHorizontal()->parameters->FindParam((void*)&u_g_obs,"u_g_obs");
+	model->DiagnosticHorizontal()->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
 
 	inputs=new ParameterInputs;
@@ -85,9 +88,9 @@
 
 	/*erase velocities: */
-	param=(Param*)femmodels[0].parameters->FindParamObject("u_g");
-	femmodels[0].parameters->DeleteObject((Object*)param);
+	param=(Param*)model->DiagnosticHorizontal()->parameters->FindParamObject("u_g");
+	model->DiagnosticHorizontal()->parameters->DeleteObject((Object*)param);
 
-	param=(Param*)femmodels[0].parameters->FindParamObject("u_g_obs");
-	femmodels[0].parameters->DeleteObject((Object*)param);
+	param=(Param*)model->DiagnosticHorizontal()->parameters->FindParamObject("u_g_obs");
+	model->DiagnosticHorizontal()->parameters->DeleteObject((Object*)param);
 
 	_printf_("initialize results:\n");
@@ -95,16 +98,17 @@
 
 	//Add output file name to parameters of femmodels[0]
-	count=femmodels[0].parameters->Size()+1;
+	count=model->DiagnosticHorizontal()->parameters->Size()+1;
 	param= new Param(count,"outputfilename",STRING);
 	param->SetString(outputfilename);
-	femmodels[0].parameters->AddObject(param);
+	model->DiagnosticHorizontal()->parameters->AddObject(param);
 
 	/*are we running the solution sequence, or a qmu wrapper around it? : */
-	femmodels[0].parameters->FindParam((void*)&qmu_analysis,"qmu_analysis");
+	model->DiagnosticHorizontal()->parameters->FindParam((void*)&qmu_analysis,"qmu_analysis");
+	
 	if(!qmu_analysis){
 
 		/*run control analysis: */
 		_printf_("call computational core:\n");
-		control_core(results,femmodels,inputs);
+		control_core(results,model,inputs);
 
 	}
@@ -114,5 +118,5 @@
 
 #ifdef _HAVE_DAKOTA_ 
-		Qmux(&femmodels[0],inputs,DiagnosticAnalysisEnum(),NoneAnalysisEnum());
+		Qmux(model,inputs,DiagnosticAnalysisEnum(),NoneAnalysisEnum());
 #else
 		throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!");
@@ -124,10 +128,10 @@
 	results->AddObject(result);
 	
-	femmodels[0].parameters->FindParam((void*)&control_type,"control_type");
+	model->DiagnosticHorizontal()->parameters->FindParam((void*)&control_type,"control_type");
 	result=new Result(results->Size()+1,0,1,"control_type",control_type);
 	results->AddObject(result);
 
 	_printf_("process results:\n");
-	ProcessResults(&results,&femmodels[0],ControlAnalysisEnum());
+	ProcessResults(&results,model,ControlAnalysisEnum());
 
 	_printf_("write results to disk:\n");
@@ -135,5 +139,5 @@
 
 	_printf_("write lock file:\n");
-	femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock");
+	model->DiagnosticHorizontal()->parameters->FindParam((void*)&waitonlock,"waitonlock");
 	if (waitonlock){
 		WriteLockFile(lockname);