Index: /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp
===================================================================
--- /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp	(revision 757)
+++ /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp	(revision 758)
@@ -35,4 +35,5 @@
 	double* vy_obs=NULL;
 	double* control_parameter=NULL;
+	double* param_g=NULL; 
 
 	double* u_g_obs=NULL;
@@ -124,10 +125,10 @@
 
 		u_g_obs=(double*)xcalloc(numberofnodes*2,sizeof(double));
-		p_g=(double*)xcalloc(numberofnodes*2,sizeof(double));
+		param_g=(double*)xcalloc(numberofnodes*2,sizeof(double));
 
 		for(i=0;i<numberofnodes;i++){
 			u_g_obs[(int)(2*partition[i]+0)]=vx_obs[i];  
 			u_g_obs[(int)(2*partition[i]+1)]=vy_obs[i];  
-			p_g[(int)(2*partition[i]+0)]=control_parameter[i];
+			param_g[(int)(2*partition[i]+0)]=control_parameter[i];
 		}
 
@@ -139,6 +140,6 @@
 
 		count++;
-		param= new Param(count,"p_g",DOUBLEVEC);
-		param->SetDoubleVec(p_g,2*numberofnodes);
+		param= new Param(count,"param_g",DOUBLEVEC);
+		param->SetDoubleVec(param_g,2*numberofnodes);
 		parameters->AddObject(param);
 
@@ -272,4 +273,5 @@
 	xfree((void**)&vx_obs);
 	xfree((void**)&vy_obs);
+	xfree((void**)&param_g);
 	xfree((void**)&pressure);
 	xfree((void**)&control_parameter);
Index: /issm/trunk/src/c/objects/OptArgs.h
===================================================================
--- /issm/trunk/src/c/objects/OptArgs.h	(revision 757)
+++ /issm/trunk/src/c/objects/OptArgs.h	(revision 758)
@@ -15,5 +15,5 @@
 	mxArray* m;
 	mxArray* inputs;
-	mxArray* p_g;
+	mxArray* param_g;
 	mxArray* u_g_obs;
 	mxArray* grad_g;
@@ -31,5 +31,5 @@
 struct OptArgs{
 	FemModel* femmodel;
-	double* p_g;
+	double* param_g;
 	double* u_g_obs;
 	double* grad_g;
Index: /issm/trunk/src/c/objects/Result.cpp
===================================================================
--- /issm/trunk/src/c/objects/Result.cpp	(revision 757)
+++ /issm/trunk/src/c/objects/Result.cpp	(revision 758)
@@ -20,4 +20,27 @@
 Result::Result(){
 	return;
+}
+
+Result::Result(const Result& result){
+
+	id=result.id;
+	time=result.time;
+	step=result.step;
+	size=result.size;
+	
+	/*copy constructor: copy dynamically allocated fields: */
+	if(result.fieldname){
+		fieldname=(char*)xmalloc((strlen(result.fieldname)+1)*sizeof(char));
+		strcpy(fieldname,result.fieldname);
+	}
+	if(result.field){
+		VecDuplicatePatch(&field,result.field);
+		dfield=NULL;
+	}
+	if(result.dfield){
+		dfield=(double*)xmalloc(result.size*sizeof(double));
+		memcpy(dfield,result.dfield,result.size*sizeof(double));
+		field=NULL;
+	}
 }
 
@@ -144,5 +167,4 @@
 	/*First write field name :*/
 	length=(strlen(fieldname)+1)*sizeof(char);
-
 	fwrite(&length,sizeof(int),1,fid);
 	fwrite(fieldname,length,1,fid);
@@ -169,2 +191,4 @@
 	return new Result(*this); 
 }
+		
+
Index: /issm/trunk/src/c/objects/Result.h
===================================================================
--- /issm/trunk/src/c/objects/Result.h	(revision 757)
+++ /issm/trunk/src/c/objects/Result.h	(revision 758)
@@ -24,4 +24,5 @@
 
 		Result();
+		Result(const Result& result);
 		Result(int result_id,double result_time,int result_step,char* result_fieldname,Vec result_field);
 		Result(int result_id,double result_time,int result_step,char* result_fieldname,double* result_field,int result_size);
Index: /issm/trunk/src/c/parallel/OutputControl.cpp
===================================================================
--- /issm/trunk/src/c/parallel/OutputControl.cpp	(revision 757)
+++ /issm/trunk/src/c/parallel/OutputControl.cpp	(revision 758)
@@ -16,5 +16,5 @@
 #define __FUNCT__ "OutputControl"
 
-void OutputControl(Vec u_g,double* p_g,double* J,int nsteps,FemModel* fem,char* outputfilename){
+void OutputControl(Vec u_g,double* param_g,double* J,int nsteps,FemModel* fem,char* outputfilename){
 
 	/*intermediary: */
@@ -26,4 +26,5 @@
 	fem->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
 
+	/*Initialize results*/
 	results=new DataSet(ResultsEnum());
 
@@ -32,5 +33,5 @@
 	results->AddObject(result);
 
-	result=new Result(results->Size()+1,0,1,"param_g",p_g,numberofnodes);
+	result=new Result(results->Size()+1,0,1,"param_g",param_g,numberofnodes);
 	results->AddObject(result);
 
@@ -40,5 +41,5 @@
 	//process results
 	ProcessResults(&results,fem,ControlAnalysisEnum());
-	
+
 	//write results to disk:
 	OutputResults(results,outputfilename);
@@ -46,5 +47,3 @@
 	/*Free ressources:*/
 	delete results;
-
 }
-
Index: /issm/trunk/src/c/parallel/ProcessResults.cpp
===================================================================
--- /issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 757)
+++ /issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 758)
@@ -45,4 +45,7 @@
 	/*fem prognostic models: */
 	FemModel* fem_p=NULL;
+
+	/*fem control models: */
+	FemModel* fem_c=NULL;
 
 	/*some parameters*/
@@ -62,7 +65,9 @@
 	double* p_g_serial=NULL;
 	double* pressure=NULL;
-	double* parameter=NULL;
 	double* partition=NULL;
 	double  yts;
+
+	double* param_g=NULL;
+	double* parameter=NULL;
 
 	Vec     h_g=NULL;
@@ -108,4 +113,9 @@
 	if(analysis_type==ThermalAnalysisEnum()){
 		fem_t=fems+0;
+	}
+
+	if(analysis_type==ControlAnalysisEnum()){
+		fem_dh=fems+0; //load u_g
+		fem_c=fems+0;  //load param_g
 	}
 
@@ -209,5 +219,4 @@
 			xfree((void**)&partition);
 		}
-
 		else if(strcmp(result->GetFieldName(),"p_g")==0){
 			/*easy, p_g is of size numberofnodes, on 1 dof, just repartition: */
@@ -244,5 +253,4 @@
 			xfree((void**)&partition);
 		}
-
 		else if(strcmp(result->GetFieldName(),"t_g")==0){
 			/*easy, t_g is of size numberofnodes, on 1 dof, just repartition: */
@@ -266,5 +274,4 @@
 			xfree((void**)&partition);
 		}
-
 		else if(strcmp(result->GetFieldName(),"m_g")==0){
 			/*easy, m_g is of size numberofnodes, on 1 dof, just repartition: */
@@ -289,5 +296,4 @@
 			xfree((void**)&partition);
 		}
-
 		else if(strcmp(result->GetFieldName(),"h_g")==0){
 			/*easy, h_g is of size numberofnodes, on 1 dof, just repartition: */
@@ -311,16 +317,14 @@
 			xfree((void**)&partition);
 		}
-
 		else if(strcmp(result->GetFieldName(),"param_g")==0){
 			/*easy, param_g is of size numberofnodes, on 1 dof, just repartition: */
-			result->GetField(&p_g);
-			VecToMPISerial(&p_g_serial,p_g);
-			fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
-			VecToMPISerial(&partition,fem_p->partition);
+			result->GetField(&param_g);
+			fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+			VecToMPISerial(&partition,fem_dh->partition);
 
 			parameter=(double*)xmalloc(numberofnodes*sizeof(double));
 
 			for(i=0;i<numberofnodes;i++){
-				parameter[i]=p_g_serial[(int)partition[i]];
+				parameter[i]=param_g[2*(int)partition[i]];
 			}
 			
@@ -330,11 +334,9 @@
 
 			/*do some cleanup: */
-			xfree((void**)&p_g_serial);
-			xfree((void**)&partition);
-		}
-
+			xfree((void**)&partition);
+		}
 		else{
 			/*Just copy the result into the new results dataset: */
-			newresults->AddObject(result);
+			newresults->AddObject(result->copy());
 		}
 	}
Index: /issm/trunk/src/c/parallel/control.cpp
===================================================================
--- /issm/trunk/src/c/parallel/control.cpp	(revision 757)
+++ /issm/trunk/src/c/parallel/control.cpp	(revision 758)
@@ -43,5 +43,5 @@
 	double*    maxiter=NULL;
 	double  tolx;
-	double*   p_g=NULL;
+	double*   param_g=NULL;
 	Vec       grad_g=NULL;
 	Vec       new_grad_g=NULL;
@@ -90,5 +90,5 @@
 	femmodel.parameters->FindParam((void*)&mincontrolconstraint,"mincontrolconstraint");
 	femmodel.parameters->FindParam((void*)&maxcontrolconstraint,"maxcontrolconstraint");
-	femmodel.parameters->FindParam((void*)&p_g,"p_g");
+	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");
@@ -107,5 +107,5 @@
 	
 	/*erase useless parameters: */
-	param=(Param*)femmodel.parameters->FindParamObject("p_g");
+	param=(Param*)femmodel.parameters->FindParamObject("param_g");
 	femmodel.parameters->DeleteObject((Object*)param);
 
@@ -121,5 +121,5 @@
 			
 		_printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
-		inputs->Add(control_type,p_g,2,numberofnodes);
+		inputs->Add(control_type,param_g,2,numberofnodes);
 		inputs->Add("fit",fit[n]);
 
@@ -139,5 +139,5 @@
 	
 		_printf_("%s\n","      optimizing along gradient direction...");
-		optargs.femmodel=&femmodel; optargs.p_g=p_g; optargs.u_g_obs=u_g_obs; optargs.grad_g=grad_g_double; optargs.inputs=inputs;optargs.n=n;
+		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);
@@ -145,9 +145,12 @@
 	
 		_printf_("%s\n","      updating parameter using optimized search scalar...");
-		for(i=0;i<gsize;i++)p_g[i]=p_g[i]+search_scalar*optscal[n]*grad_g_double[i];
+		PetscSynchronizedPrintf(MPI_COMM_WORLD,"gsize=%i\n",gsize);
+		PetscSynchronizedFlush(MPI_COMM_WORLD);
+
+		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( p_g,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
+		ControlConstrainx( param_g,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
 		_printf_("%s\n","      done.");
 	
@@ -155,8 +158,8 @@
 		if (((n+1)%5)==0){
 			_printf_("%s\n","      saving temporary results...");
-			inputs->Add(control_type,p_g,2,numberofnodes);
+			inputs->Add(control_type,param_g,2,numberofnodes);
 			inputs->Add("fit",fit[n]);
 			diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type);
-			OutputControl(u_g,p_g,J,nsteps,&femmodel,outputfilename);
+			OutputControl(u_g,param_g,J,nsteps,&femmodel,outputfilename);
 			_printf_("%s\n","      done.");
 		}
@@ -170,5 +173,5 @@
 	/*Write results to disk: */
 	_printf_("%s\n","      preparing final velocity solution...");
-	inputs->Add(control_type,p_g,2,numberofnodes);
+	inputs->Add(control_type,param_g,2,numberofnodes);
 	inputs->Add("fit",fit[n]);
 
@@ -179,5 +182,5 @@
 	
 	_printf_("%s\n","      saving final results...");
-	OutputControl(u_g,p_g,J,nsteps,&femmodel,outputfilename);
+	OutputControl(u_g,param_g,J,nsteps,&femmodel,outputfilename);
 	_printf_("%s\n","      done.");
 
Index: /issm/trunk/src/c/parallel/objectivefunctionC.cpp
===================================================================
--- /issm/trunk/src/c/parallel/objectivefunctionC.cpp	(revision 757)
+++ /issm/trunk/src/c/parallel/objectivefunctionC.cpp	(revision 758)
@@ -24,5 +24,5 @@
 	/*parameters: */
 	FemModel* femmodel=NULL;
-	double* p_g=NULL;
+	double* param_g=NULL;
 	double* u_g_obs=NULL;
 	double* grad_g=NULL;
@@ -37,5 +37,5 @@
 	double  maxcontrolconstraint;
 	char*   control_type=NULL;
-	double* p_g_copy=NULL;
+	double* param_g_copy=NULL;
 	int     analysis_type;
 	int     sub_analysis_type;
@@ -47,5 +47,5 @@
 	/*Recover parameters: */
 	femmodel=optargs->femmodel;
-	p_g=optargs->p_g;
+	param_g=optargs->param_g;
 	u_g_obs=optargs->u_g_obs;
 	grad_g=optargs->grad_g;
@@ -63,16 +63,16 @@
 	femmodel->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
 
-	/*First copy p_g so we don't modify it: */
-	p_g_copy=(double*)xmalloc(gsize*sizeof(double));
-	memcpy(p_g_copy,p_g,gsize*sizeof(double));
+	/*First copy param_g so we don't modify it: */
+	param_g_copy=(double*)xmalloc(gsize*sizeof(double));
+	memcpy(param_g_copy,param_g,gsize*sizeof(double));
 
-	/*First, update p_g using search_scalar: */
-	for(i=0;i<gsize;i++)p_g_copy[i]=p_g_copy[i]+search_scalar*optscal[n]*grad_g[i];
+	/*First, update param_g using search_scalar: */
+	for(i=0;i<gsize;i++)param_g_copy[i]=param_g_copy[i]+search_scalar*optscal[n]*grad_g[i];
 
 	/*Constrain:*/
-	ControlConstrainx( p_g_copy,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
+	ControlConstrainx( param_g_copy,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
 
 	/*Add new parameter to inputs: */
-	inputs->Add(control_type,p_g_copy,2,numberofnodes);
+	inputs->Add(control_type,param_g_copy,2,numberofnodes);
 
 	//Run diagnostic with updated parameters.
@@ -89,5 +89,5 @@
 	xfree((void**)&optscal);
 	xfree((void**)&control_type);
-	xfree((void**)&p_g_copy);
+	xfree((void**)&param_g_copy);
 	xfree((void**)&u_g_double);
 
Index: /issm/trunk/src/c/shared/Numerics/OptFunc.cpp
===================================================================
--- /issm/trunk/src/c/shared/Numerics/OptFunc.cpp	(revision 757)
+++ /issm/trunk/src/c/shared/Numerics/OptFunc.cpp	(revision 758)
@@ -32,5 +32,5 @@
 	inputs[1]=optargs->m;
 	inputs[2]=optargs->inputs;
-	inputs[3]=optargs->p_g;
+	inputs[3]=optargs->param_g;
 	inputs[4]=optargs->u_g_obs;
 	inputs[5]=optargs->grad_g;