Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 1036)
+++ /issm/trunk/src/c/Makefile.am	(revision 1037)
@@ -559,4 +559,5 @@
 					./parallel/WriteLockFile.cpp\
 					./parallel/control.cpp\
+					./parallel/control_core.cpp\
 					./parallel/objectivefunctionC.cpp\
 					./parallel/GradJCompute.cpp\
Index: /issm/trunk/src/c/parallel/ProcessResults.cpp
===================================================================
--- /issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 1036)
+++ /issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 1037)
@@ -146,4 +146,10 @@
 		fem_dh=fems+0; //load u_g
 		fem_c=fems+0;  //load param_g
+
+		/*some flags needed: */
+		fem_dh->parameters->FindParam((void*)&dim,"dim");
+		fem_dh->parameters->FindParam((void*)&ishutter,"ishutter");
+		fem_dh->parameters->FindParam((void*)&isstokes,"isstokes");
+		fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn");
 	}
 
Index: /issm/trunk/src/c/parallel/control.cpp
===================================================================
--- /issm/trunk/src/c/parallel/control.cpp	(revision 1036)
+++ /issm/trunk/src/c/parallel/control.cpp	(revision 1037)
@@ -17,7 +17,4 @@
 int main(int argc,char* *argv){
 
-	int n;
-	int i;
-
 	/*I/O: */
 	FILE* fid=NULL;
@@ -25,34 +22,19 @@
 	char* outputfilename=NULL;
 	char* lockname=NULL;
-	int   analysis_type;
-	int   sub_analysis_type;
+	int   numberofnodes;
+	int   qmu_analysis=0;
+
+	/*Fem models : */
+	FemModel femmodels[1];
+
+	/*Results: */
+	DataSet* results=NULL;
+
+	ParameterInputs* inputs=NULL;
+	int waitonlock=0;
 
 	/*Intermediary: */
-	FemModel femmodel;
-	Vec u_g=NULL;
-	ParameterInputs* inputs=NULL;
-	int waitonlock=0;
-	double search_scalar;
-	char* control_type=NULL;
-	double* fit=NULL;
-	double* optscal=NULL;
-	double* u_g_obs=NULL;
 	double* u_g_initial=NULL;
-	int  gsize;
-	int  numberofnodes;
-	double*    maxiter=NULL;
-	double  tolx;
-	double*   param_g=NULL;
-	Vec       grad_g=NULL;
-	Vec       new_grad_g=NULL;
-	Vec       grad_g_old=NULL;
-	double*   grad_g_double=NULL;
-	double    mincontrolconstraint;
-	double    maxcontrolconstraint;
-	int       nsteps;
-	double*   J=NULL;
-	OptArgs   optargs;
-	OptPars   optpars;
-	Param*    param=NULL;
+	Param*  param=NULL;
 
 
@@ -76,116 +58,55 @@
 	/*Open handle to data on disk: */
 	fid=pfopen(inputfilename,"rb");
-	
+
 	_printf_("read and create finite element model:\n");
-	CreateFemModel(&femmodel,fid,"control","");
+	_printf_("\n   reading diagnostic horiz model data:\n");
+	CreateFemModel(&femmodels[0],fid,"control","");
 
-	/*Recover parameters used throughout the solution:*/
-	femmodel.parameters->FindParam((void*)&nsteps,"nsteps");
-	femmodel.parameters->FindParam((void*)&control_type,"control_type");
-	femmodel.parameters->FindParam((void*)&fit,"fit");
-	femmodel.parameters->FindParam((void*)&optscal,"optscal");
-	femmodel.parameters->FindParam((void*)&maxiter,"maxiter");
-	femmodel.parameters->FindParam((void*)&tolx,"tolx");
-	femmodel.parameters->FindParam((void*)&waitonlock,"waitonlock");
-	femmodel.parameters->FindParam((void*)&mincontrolconstraint,"mincontrolconstraint");
-	femmodel.parameters->FindParam((void*)&maxcontrolconstraint,"maxcontrolconstraint");
-	femmodel.parameters->FindParam((void*)&param_g,"param_g");
-	femmodel.parameters->FindParam((void*)&u_g_obs,"u_g_obs");
-	femmodel.parameters->FindParam((void*)&analysis_type,"analysis_type");
-	femmodel.parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
-	gsize=femmodel.nodes->NumberOfDofs();
-
-	/*Initialize misfit: */
-	J=(double*)xmalloc(nsteps*sizeof(double));
-
-	/*Initialize inputs:*/
-	femmodel.parameters->FindParam((void*)&u_g_initial,"u_g");
-	femmodel.parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+	_printf_("initialize inputs:\n");
+	femmodels[0].parameters->FindParam((void*)&u_g_initial,"u_g");
+	femmodels[0].parameters->FindParam((void*)&numberofnodes,"numberofnodes");
 
 	inputs=new ParameterInputs;
 	inputs->Add("velocity",u_g_initial,3,numberofnodes);
-	
-	/*erase useless parameters: */
-	param=(Param*)femmodel.parameters->FindParamObject("param_g");
-	femmodel.parameters->DeleteObject((Object*)param);
 
-	param=(Param*)femmodel.parameters->FindParamObject("u_g_obs");
-	femmodel.parameters->DeleteObject((Object*)param);
+	/*erase velocities: */
+	param=(Param*)femmodels[0].parameters->FindParamObject("u_g");
+	femmodels[0].parameters->DeleteObject((Object*)param);
 
-	param=(Param*)femmodel.parameters->FindParamObject("u_g");
-	femmodel.parameters->DeleteObject((Object*)param);
+	_printf_("initialize results:\n");
+	results=new DataSet(ResultsEnum());
 
+	/*are we running the solution sequence, or a qmu wrapper around it? : */
+	femmodels[0].parameters->FindParam((void*)&qmu_analysis,"qmu_analysis");
+	if(!qmu_analysis){
 
-	/*Start looping: */
-	for(n=0;n<nsteps;n++){
-			
-		_printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
-		inputs->Add(control_type,param_g,2,numberofnodes);
-		inputs->Add("fit",fit[n]);
+		/*run control analysis: */
+		_printf_("call computational core:\n");
+		control_core(results,femmodels,inputs);
 
-		/*Update parameters: */
-		UpdateFromInputsx(femmodel.elements,femmodel.nodes,femmodel.loads, femmodel.materials,inputs);
+	}
+	else{
+		/*run qmu analysis: */
+		_printf_("calling qmu analysis on control core:\n");
 
-		_printf_("%s\n","      computing gradJ...");
-		grad_g=GradJCompute(inputs,&femmodel,u_g_obs);
-		_printf_("%s\n","      done.");
-			
-		_printf_("%s\n","      normalizing directions...");
-		Orthx(&new_grad_g,grad_g,grad_g_old);
-		VecFree(&grad_g); grad_g=new_grad_g;
-		VecFree(&grad_g_old); grad_g_old=grad_g;
-		VecToMPISerial(&grad_g_double,grad_g);
-		_printf_("%s\n","      done.");
-	
-		_printf_("%s\n","      optimizing along gradient direction...");
-		optargs.femmodel=&femmodel; optargs.param_g=param_g; optargs.u_g_obs=u_g_obs; optargs.grad_g=grad_g_double; optargs.inputs=inputs;optargs.n=n;
-		optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx; optpars.maxiter=(int)maxiter[n];
-		BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
-		_printf_("%s\n","      done.");
-	
-		_printf_("%s\n","      updating parameter using optimized search scalar...");
-		for(i=0;i<gsize;i++)param_g[i]=param_g[i]+search_scalar*optscal[n]*grad_g_double[i];
-		_printf_("%s\n","      done.");
-
-		_printf_("%s\n","      constraining the new distribution...");    
-		ControlConstrainx( param_g,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
-		_printf_("%s\n","      done.");
-
-		//some temporary saving 
-		/*if (((n+1)%5)==0){
-			_printf_("%s\n","      saving temporary results...");
-			inputs->Add(control_type,param_g,2,numberofnodes);
-			inputs->Add("fit",fit[n]);
-			UpdateFromInputsx(femmodel.elements,femmodel.nodes,femmodel.loads, femmodel.materials,inputs);
-			diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type);
-			OutputControl(u_g,param_g,J,nsteps,&femmodel,outputfilename);
-			_printf_("%s\n","      done.");
-		}*/
-
-		_printf_("%s%i%s%g\n","      value of misfit J after optimization #",n,": ",J[n]);
-		
-		/*some freeing:*/
-		xfree((void**)&grad_g_double);
+#ifdef _HAVE_DAKOTA_ 
+		Qmux(&femmodels[0],inputs,DiagnosticAnalysisEnum(),NoneAnalysisEnum());
+#else
+		throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!");
+#endif
 	}
 
-	/*Write results to disk: */
-	_printf_("%s\n","      preparing final velocity solution...");
-	inputs->Add(control_type,param_g,2,numberofnodes);
-	inputs->Add("fit",fit[n]);
+	_printf_("process results:\n");
+	ProcessResults(&results,&femmodels[0],ControlAnalysisEnum());
 
-	/*Update parameters: */
-	UpdateFromInputsx(femmodel.elements,femmodel.nodes,femmodel.loads, femmodel.materials,inputs);
-	
-	diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type);
-	
-	_printf_("%s\n","      saving final results...");
-	OutputControl(u_g,param_g,J,nsteps,&femmodel,outputfilename);
-	_printf_("%s\n","      done.");
+	_printf_("write results to disk:\n");
+	OutputResults(results,outputfilename);
 
-	/*Write lock file if requested: */
+	_printf_("write lock file:\n");
+	femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock");
 	if (waitonlock){
 		WriteLockFile(lockname);
 	}
-			
+
 	_printf_("closing MPI and Petsc\n");
 	PetscFinalize(); 
Index: /issm/trunk/src/c/parallel/control_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/control_core.cpp	(revision 1037)
+++ /issm/trunk/src/c/parallel/control_core.cpp	(revision 1037)
@@ -0,0 +1,155 @@
+/*!\file: control_core.cpp
+ * \brief: core of the control solution 
+ */ 
+
+#undef __FUNCT__ 
+#define __FUNCT__ "control_core"
+
+#include "./parallel.h"
+#include "../issm.h"
+
+void control_core(DataSet* results,FemModel* fems, ParameterInputs* inputs){
+
+	extern int my_rank;
+
+	/*fem models: */
+	FemModel* fem_dh=NULL;
+
+	/*output: */
+	Result* result=NULL;
+
+	/*Intermediary: */
+	Vec     u_g=NULL;
+	double  search_scalar;
+	char*   control_type=NULL;
+	double* fit=NULL;
+	double* optscal=NULL;
+	double* u_g_obs=NULL;
+	int     gsize;
+	double* maxiter=NULL;
+	double  tolx;
+	double* param_g=NULL;
+	Vec     grad_g=NULL;
+	Vec     new_grad_g=NULL;
+	Vec     grad_g_old=NULL;
+	double* grad_g_double=NULL;
+	double  mincontrolconstraint;
+	double  maxcontrolconstraint;
+	int     nsteps,n,i;
+	double* J=NULL;
+	OptArgs optargs;
+	OptPars optpars;
+	Param*  param=NULL;
+
+	/*flags: */
+	int   analysis_type;
+	int   sub_analysis_type;
+	int debug=0;
+	int dim=-1;
+	int ishutter=0;
+	int ismacayealpattyn=0;
+	int isstokes=0;
+	int numberofdofspernode_sl;
+	int numberofdofspernode_dh;
+	int numberofdofspernode_ds;
+	int numberofnodes;
+
+	/*recover fem models: */
+	fem_dh=fems+0;
+
+	//first recover parameters common to all solutions
+	fem_dh->parameters->FindParam((void*)&debug,"debug");
+	fem_dh->parameters->FindParam((void*)&dim,"dim");
+	fem_dh->parameters->FindParam((void*)&ishutter,"ishutter");
+	fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn");
+	fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+	fem_dh->parameters->FindParam((void*)&isstokes,"isstokes");
+
+	/*Recover parameters used throughout the solution:*/
+	fem_dh->parameters->FindParam((void*)&nsteps,"nsteps");
+	fem_dh->parameters->FindParam((void*)&control_type,"control_type");
+	fem_dh->parameters->FindParam((void*)&fit,"fit");
+	fem_dh->parameters->FindParam((void*)&optscal,"optscal");
+	fem_dh->parameters->FindParam((void*)&maxiter,"maxiter");
+	fem_dh->parameters->FindParam((void*)&tolx,"tolx");
+	fem_dh->parameters->FindParam((void*)&mincontrolconstraint,"mincontrolconstraint");
+	fem_dh->parameters->FindParam((void*)&maxcontrolconstraint,"maxcontrolconstraint");
+	fem_dh->parameters->FindParam((void*)&param_g,"param_g");
+	fem_dh->parameters->FindParam((void*)&u_g_obs,"u_g_obs");
+	fem_dh->parameters->FindParam((void*)&analysis_type,"analysis_type");
+	fem_dh->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
+	fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+	gsize=fem_dh->nodes->NumberOfDofs();
+
+	/*Initialize misfit: */
+	J=(double*)xmalloc(nsteps*sizeof(double));
+
+	/*erase useless parameters: */
+	param=(Param*)fem_dh->parameters->FindParamObject("param_g");
+	fem_dh->parameters->DeleteObject((Object*)param);
+
+	param=(Param*)fem_dh->parameters->FindParamObject("u_g_obs");
+	fem_dh->parameters->DeleteObject((Object*)param);
+
+	param=(Param*)fem_dh->parameters->FindParamObject("u_g");
+	fem_dh->parameters->DeleteObject((Object*)param);
+
+	/*Start looping: */
+	for(n=0;n<nsteps;n++){
+
+		_printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
+		inputs->Add(control_type,param_g,2,numberofnodes);
+		inputs->Add("fit",fit[n]);
+
+		/*Update parameters: */
+		UpdateFromInputsx(fem_dh->elements,fem_dh->nodes,fem_dh->loads, fem_dh->materials,inputs);
+
+		_printf_("%s\n","      computing gradJ...");
+		grad_g=GradJCompute(inputs,fem_dh,u_g_obs);
+		_printf_("%s\n","      done.");
+
+		_printf_("%s\n","      normalizing directions...");
+		Orthx(&new_grad_g,grad_g,grad_g_old);
+		VecFree(&grad_g); grad_g=new_grad_g;
+		VecFree(&grad_g_old); grad_g_old=grad_g;
+		VecToMPISerial(&grad_g_double,grad_g);
+		_printf_("%s\n","      done.");
+
+		_printf_("%s\n","      optimizing along gradient direction...");
+		optargs.femmodel=fem_dh; optargs.param_g=param_g; optargs.u_g_obs=u_g_obs; optargs.grad_g=grad_g_double; optargs.inputs=inputs;optargs.n=n;
+		optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx; optpars.maxiter=(int)maxiter[n];
+		BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
+		_printf_("%s\n","      done.");
+
+		_printf_("%s\n","      updating parameter using optimized search scalar...");
+		for(i=0;i<gsize;i++)param_g[i]=param_g[i]+search_scalar*optscal[n]*grad_g_double[i];
+		_printf_("%s\n","      done.");
+
+		_printf_("%s\n","      constraining the new distribution...");    
+		ControlConstrainx( param_g,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
+		_printf_("%s\n","      done.");
+
+		_printf_("%s%i%s%g\n","      value of misfit J after optimization #",n,": ",J[n]);
+
+		/*some freeing:*/
+		xfree((void**)&grad_g_double);
+	}
+
+	/*Write results to disk: */
+	_printf_("%s\n","      preparing final velocity solution...");
+	inputs->Add(control_type,param_g,2,numberofnodes);
+	inputs->Add("fit",fit[n]);
+
+	/*Update parameters: */
+	UpdateFromInputsx(fem_dh->elements,fem_dh->nodes,fem_dh->loads, fem_dh->materials,inputs);
+
+	diagnostic_core_nonlinear(&u_g,NULL,NULL,fem_dh,inputs,analysis_type,sub_analysis_type);
+
+	/*Plug results into output dataset: */
+	result=new Result(results->Size()+1,0,1,"u_g",u_g);
+	results->AddObject(result);
+	result=new Result(results->Size()+1,0,1,"param_g",param_g,gsize);
+	results->AddObject(result);
+	result=new Result(results->Size()+1,0,1,"J",J,nsteps);
+	results->AddObject(result);
+}
Index: /issm/trunk/src/c/parallel/parallel.h
===================================================================
--- /issm/trunk/src/c/parallel/parallel.h	(revision 1036)
+++ /issm/trunk/src/c/parallel/parallel.h	(revision 1037)
@@ -13,4 +13,5 @@
 void diagnostic_core(DataSet* results,FemModel* fems, ParameterInputs* inputs);
 void prognostic_core(DataSet* results,FemModel* fems, ParameterInputs* inputs);
+void control_core(DataSet* results,FemModel* fems, ParameterInputs* inputs);
 
 void thermal_core(DataSet* results,FemModel* fems, ParameterInputs* inputs);
@@ -37,5 +38,4 @@
 //int ParameterUpdate(double* search_vector,int step, WorkspaceParams* workspaceparams,BatchParams* batchparams);
 void OutputResults(DataSet* results,char* filename);
-void OutputControl(Vec u_g,double* p_g,double* J,int nsteps,FemModel* fem,char* outputfilename);
 void WriteLockFile(char* filename);