Index: /issm/trunk/src/c/HoleFillerx/HoleFillerx.cpp
===================================================================
--- /issm/trunk/src/c/HoleFillerx/HoleFillerx.cpp	(revision 1906)
+++ /issm/trunk/src/c/HoleFillerx/HoleFillerx.cpp	(revision 1907)
@@ -61,6 +61,6 @@
  			if ( *(image+i*samps+j) == 0 ){ 
   				if ( (j > 3) && (j < samps-4) && (i > 3) && (i < lines-4)){
-					for ( k = 0; k < 4; k++ ){
-						for ( l = 0; l < 4; l++ ){
+					for ( k = 0; k < 30; k++ ){
+						for ( l = 0; l < 30; l++ ){
 							*(image2+samps*(i+k)+j+l)=0;
 							*(image2+samps*(i-k)+j+l)=0;
Index: /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 1906)
+++ /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 1907)
@@ -42,122 +42,130 @@
 	count=parameters->Size();
 	
-	/*control_type: */
+	//control analysis?
 	count++;
-	param= new Param(count,"control_type",STRING);
-	param->SetString(iomodel->control_type);
-	parameters->AddObject(param);
-
-	/*extrude_param: */
-	count++;
-	param= new Param(count,"extrude_param",DOUBLE);
-	if (strcmp(iomodel->control_type,"drag")==0)   param->SetDouble(0);
-	else if (strcmp(iomodel->control_type,"B")==0) param->SetDouble(1);
-	else throw ErrorException(__FUNCT__,exprintf("control_type %s not supported yet!",iomodel->control_type));
-	parameters->AddObject(param);
-
-	/*nsteps: */
-	count++;
-	param= new Param(count,"nsteps",INTEGER);
-	param->SetInteger(iomodel->nsteps);
-	parameters->AddObject(param);
-
-	/*tolx: */
-	count++;
-	param= new Param(count,"tolx",DOUBLE);
-	param->SetDouble(iomodel->tolx);
-	parameters->AddObject(param);
-
-	/*mincontrolconstraint: */
-	count++;
-	param= new Param(count,"mincontrolconstraint",DOUBLE);
-	param->SetDouble(iomodel->mincontrolconstraint);
-	parameters->AddObject(param);
-
-	/*maxcontrolconstraint: */
-	count++;
-	param= new Param(count,"maxcontrolconstraint",DOUBLE);
-	param->SetDouble(iomodel->maxcontrolconstraint);
+	param= new Param(count,"control_analysis",INTEGER);
+	param->SetInteger(iomodel->control_analysis);
 	parameters->AddObject(param);
 	
-	/*epsvel: */
-	count++;
-	param= new Param(count,"epsvel",DOUBLE);
-	param->SetDouble(iomodel->epsvel);
-	parameters->AddObject(param);
-	
-	/*meanvel: */
-	count++;
-	param= new Param(count,"meanvel",DOUBLE);
-	param->SetDouble(iomodel->meanvel);
-	parameters->AddObject(param);
+	if(iomodel->control_analysis){
+		/*control_type: */
+		count++;
+		param= new Param(count,"control_type",STRING);
+		param->SetString(iomodel->control_type);
+		parameters->AddObject(param);
 
-	/*Now, recover fit, optscal and maxiter as vectors: */
-	IoModelFetchData((void**)&iomodel->fit,NULL,NULL,iomodel_handle,"fit","Matrix","Mat");
-	IoModelFetchData((void**)&iomodel->optscal,NULL,NULL,iomodel_handle,"optscal","Matrix","Mat");
-	IoModelFetchData((void**)&iomodel->maxiter,NULL,NULL,iomodel_handle,"maxiter","Matrix","Mat");
+		/*extrude_param: */
+		count++;
+		param= new Param(count,"extrude_param",DOUBLE);
+		if (strcmp(iomodel->control_type,"drag")==0)   param->SetDouble(0);
+		else if (strcmp(iomodel->control_type,"B")==0) param->SetDouble(1);
+		else throw ErrorException(__FUNCT__,exprintf("control_type %s not supported yet!",iomodel->control_type));
+		parameters->AddObject(param);
 
-	count++;
-	param= new Param(count,"fit",DOUBLEVEC);
-	param->SetDoubleVec(iomodel->fit,iomodel->nsteps);
-	parameters->AddObject(param);
+		/*nsteps: */
+		count++;
+		param= new Param(count,"nsteps",INTEGER);
+		param->SetInteger(iomodel->nsteps);
+		parameters->AddObject(param);
 
-	count++;
-	param= new Param(count,"optscal",DOUBLEVEC);
-	param->SetDoubleVec(iomodel->optscal,iomodel->nsteps);
-	parameters->AddObject(param);
+		/*tolx: */
+		count++;
+		param= new Param(count,"tolx",DOUBLE);
+		param->SetDouble(iomodel->tolx);
+		parameters->AddObject(param);
 
-	count++;
-	param= new Param(count,"maxiter",DOUBLEVEC);
-	param->SetDoubleVec(iomodel->maxiter,iomodel->nsteps);
-	parameters->AddObject(param);
+		/*mincontrolconstraint: */
+		count++;
+		param= new Param(count,"mincontrolconstraint",DOUBLE);
+		param->SetDouble(iomodel->mincontrolconstraint);
+		parameters->AddObject(param);
 
-	xfree((void**)&iomodel->fit);
-	xfree((void**)&iomodel->optscal);
-	xfree((void**)&iomodel->maxiter);
+		/*maxcontrolconstraint: */
+		count++;
+		param= new Param(count,"maxcontrolconstraint",DOUBLE);
+		param->SetDouble(iomodel->maxcontrolconstraint);
+		parameters->AddObject(param);
+		
+		/*epsvel: */
+		count++;
+		param= new Param(count,"epsvel",DOUBLE);
+		param->SetDouble(iomodel->epsvel);
+		parameters->AddObject(param);
+		
+		/*meanvel: */
+		count++;
+		param= new Param(count,"meanvel",DOUBLE);
+		param->SetDouble(iomodel->meanvel);
+		parameters->AddObject(param);
 
-	/*Get vx, vx_obs, vy, vy_obs, and the parameter value: */
-	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");
-	IoModelFetchData((void**)&vx_obs,NULL,NULL,iomodel_handle,"vx_obs","Matrix","Mat");
-	IoModelFetchData((void**)&vy_obs,NULL,NULL,iomodel_handle,"vy_obs","Matrix","Mat");
-	IoModelFetchData((void**)&control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type,"Matrix","Mat");
+		/*Now, recover fit, optscal and maxiter as vectors: */
+		IoModelFetchData((void**)&iomodel->fit,NULL,NULL,iomodel_handle,"fit","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->optscal,NULL,NULL,iomodel_handle,"optscal","Matrix","Mat");
+		IoModelFetchData((void**)&iomodel->maxiter,NULL,NULL,iomodel_handle,"maxiter","Matrix","Mat");
 
-	u_g=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double));
-	if(vx)for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+0]=vx[i]/iomodel->yts;
-	if(vy)for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+1]=vy[i]/iomodel->yts;
-	if(vz)for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+2]=vz[i]/iomodel->yts;
+		count++;
+		param= new Param(count,"fit",DOUBLEVEC);
+		param->SetDoubleVec(iomodel->fit,iomodel->nsteps);
+		parameters->AddObject(param);
 
-	count++;
-	param= new Param(count,"u_g",DOUBLEVEC);
-	param->SetDoubleVec(u_g,3*iomodel->numberofnodes,3);
-	parameters->AddObject(param);
+		count++;
+		param= new Param(count,"optscal",DOUBLEVEC);
+		param->SetDoubleVec(iomodel->optscal,iomodel->nsteps);
+		parameters->AddObject(param);
 
-	u_g_obs=(double*)xcalloc(iomodel->numberofnodes*2,sizeof(double));
-	if(vx_obs)for(i=0;i<iomodel->numberofnodes;i++)u_g_obs[2*i+0]=vx_obs[i]/iomodel->yts;
-	if(vy_obs)for(i=0;i<iomodel->numberofnodes;i++)u_g_obs[2*i+1]=vy_obs[i]/iomodel->yts;
+		count++;
+		param= new Param(count,"maxiter",DOUBLEVEC);
+		param->SetDoubleVec(iomodel->maxiter,iomodel->nsteps);
+		parameters->AddObject(param);
 
-	count++;
-	param= new Param(count,"u_g_obs",DOUBLEVEC);
-	param->SetDoubleVec(u_g_obs,2*iomodel->numberofnodes,2);
-	parameters->AddObject(param);
-	
-	param_g=(double*)xcalloc(iomodel->numberofnodes,sizeof(double));
-	for(i=0;i<iomodel->numberofnodes;i++)param_g[i]=control_parameter[i];
+		xfree((void**)&iomodel->fit);
+		xfree((void**)&iomodel->optscal);
+		xfree((void**)&iomodel->maxiter);
 
-	count++;
-	param= new Param(count,"param_g",DOUBLEVEC);
-	param->SetDoubleVec(param_g,iomodel->numberofnodes,1);
-	parameters->AddObject(param);
+		/*Get vx, vx_obs, vy, vy_obs, and the parameter value: */
+		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");
+		IoModelFetchData((void**)&vx_obs,NULL,NULL,iomodel_handle,"vx_obs","Matrix","Mat");
+		IoModelFetchData((void**)&vy_obs,NULL,NULL,iomodel_handle,"vy_obs","Matrix","Mat");
+		IoModelFetchData((void**)&control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type,"Matrix","Mat");
 
-	xfree((void**)&vx);
-	xfree((void**)&vy);
-	xfree((void**)&vz);
-	xfree((void**)&u_g);
-	xfree((void**)&vx_obs);
-	xfree((void**)&vy_obs);
-	xfree((void**)&u_g_obs);
-	xfree((void**)&param_g);
-	xfree((void**)&control_parameter);
+		u_g=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double));
+		if(vx)for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+0]=vx[i]/iomodel->yts;
+		if(vy)for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+1]=vy[i]/iomodel->yts;
+		if(vz)for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+2]=vz[i]/iomodel->yts;
+
+		count++;
+		param= new Param(count,"u_g",DOUBLEVEC);
+		param->SetDoubleVec(u_g,3*iomodel->numberofnodes,3);
+		parameters->AddObject(param);
+
+		u_g_obs=(double*)xcalloc(iomodel->numberofnodes*2,sizeof(double));
+		if(vx_obs)for(i=0;i<iomodel->numberofnodes;i++)u_g_obs[2*i+0]=vx_obs[i]/iomodel->yts;
+		if(vy_obs)for(i=0;i<iomodel->numberofnodes;i++)u_g_obs[2*i+1]=vy_obs[i]/iomodel->yts;
+
+		count++;
+		param= new Param(count,"u_g_obs",DOUBLEVEC);
+		param->SetDoubleVec(u_g_obs,2*iomodel->numberofnodes,2);
+		parameters->AddObject(param);
+		
+		param_g=(double*)xcalloc(iomodel->numberofnodes,sizeof(double));
+		for(i=0;i<iomodel->numberofnodes;i++)param_g[i]=control_parameter[i];
+
+		count++;
+		param= new Param(count,"param_g",DOUBLEVEC);
+		param->SetDoubleVec(param_g,iomodel->numberofnodes,1);
+		parameters->AddObject(param);
+
+		xfree((void**)&vx);
+		xfree((void**)&vy);
+		xfree((void**)&vz);
+		xfree((void**)&u_g);
+		xfree((void**)&vx_obs);
+		xfree((void**)&vy_obs);
+		xfree((void**)&u_g_obs);
+		xfree((void**)&param_g);
+		xfree((void**)&control_parameter);
+	}
 
 	/*Assign output pointer: */
Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp	(revision 1906)
+++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp	(revision 1907)
@@ -23,42 +23,8 @@
 	CreateParameters(pparameters,iomodel,iomodel_handle);
 	CreateParametersQmu(pparameters,iomodel,iomodel_handle);
+	CreateParametersControl(pparameters,iomodel,iomodel_handle);
 
 	/*This is just a high level driver: */
-	if (iomodel->analysis_type==ControlAnalysisEnum()){
-
-		if (iomodel->sub_analysis_type==HorizAnalysisEnum()){
-
-			CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
-			CreateConstraintsDiagnosticHoriz(pconstraints,iomodel,iomodel_handle);
-			CreateLoadsDiagnosticHoriz(ploads,iomodel,iomodel_handle);
-			CreateParametersControl(pparameters,iomodel,iomodel_handle);
-
-		}
-		else if (iomodel->sub_analysis_type==VertAnalysisEnum()){
-
-			CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
-			CreateConstraintsDiagnosticVert(pconstraints,iomodel,iomodel_handle);
-			CreateLoadsDiagnosticVert(ploads,iomodel,iomodel_handle);
-
-		}
-		else if (iomodel->sub_analysis_type==StokesAnalysisEnum()){
-
-			CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
-			CreateConstraintsDiagnosticStokes(pconstraints,iomodel,iomodel_handle);
-			CreateLoadsDiagnosticStokes(ploads,iomodel,iomodel_handle);
-			CreateParametersControl(pparameters,iomodel,iomodel_handle);
-
-		}
-		else if (iomodel->sub_analysis_type==HutterAnalysisEnum()){
-
-			CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
-			CreateConstraintsDiagnosticHutter(pconstraints,iomodel,iomodel_handle);
-			CreateLoadsDiagnosticHutter(ploads,iomodel,iomodel_handle);
-			CreateParametersControl(pparameters,iomodel,iomodel_handle);
-
-		}
-
-	}
-	else if (iomodel->analysis_type==DiagnosticAnalysisEnum()){
+	if (iomodel->analysis_type==DiagnosticAnalysisEnum()){
 
 		if (iomodel->sub_analysis_type==HorizAnalysisEnum()){
Index: /issm/trunk/src/c/ModelProcessorx/IoModel.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/IoModel.cpp	(revision 1906)
+++ /issm/trunk/src/c/ModelProcessorx/IoModel.cpp	(revision 1907)
@@ -38,4 +38,5 @@
 	iomodel->sub_analysis_type=0;
 	iomodel->qmu_analysis=0;
+	iomodel->control_analysis=0;
 	iomodel->solverstring=NULL;
 	iomodel->numberofresponses=0;
@@ -286,4 +287,5 @@
 	IoModelFetchData((void**)&iomodel->sub_analysis_type,NULL,NULL,iomodel_handle,"sub_analysis_type","Integer",NULL); 
 	IoModelFetchData((void**)&iomodel->qmu_analysis,NULL,NULL,iomodel_handle,"qmu_analysis","Integer",NULL); 
+	IoModelFetchData((void**)&iomodel->control_analysis,NULL,NULL,iomodel_handle,"control_analysis","Integer",NULL); 
 	IoModelFetchData((void**)&iomodel->meshtype,NULL,NULL,iomodel_handle,"type","String",NULL);
 	/*!Get numberofelements and numberofnodes: */
Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 1906)
+++ /issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 1907)
@@ -21,4 +21,5 @@
 	int     sub_analysis_type;
 	int     qmu_analysis;
+	int     control_analysis;
 	char*   solverstring;
 
Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 1906)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 1907)
@@ -1119,6 +1119,8 @@
 	if(inputs->Recover("temperature",&temperature_list[0],1,dofs,6,(void**)nodes)){
 		if(matice && !collapse){
-			B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2])
-						+Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0;
+			//B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2])
+			//			+Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0;
+			temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2]+temperature_list[3]+temperature_list[4]+temperature_list[5])/6.0;
+			B_average=Paterson(temperature_average);
 			matice->SetB(B_average);
 		}
@@ -1127,6 +1129,8 @@
 	if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,6,(void**)nodes)){
 		if(matice && collapse){
-			B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2])
-						+Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0;
+			temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2]+temperature_list[3]+temperature_list[4]+temperature_list[5])/6.0;
+			B_average=Paterson(temperature_average);
+			//B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2])
+			//			+Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0;
 			matice->SetB(B_average);
 		}
Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m
===================================================================
--- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 1906)
+++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 1907)
@@ -335,5 +335,5 @@
 
 %CONTROL
-if md.analysis_type==ControlAnalysisEnum,
+if md.control_analysis,
 
 	%CONTROL TYPE
Index: /issm/trunk/src/m/classes/public/solve.m
===================================================================
--- /issm/trunk/src/m/classes/public/solve.m	(revision 1906)
+++ /issm/trunk/src/m/classes/public/solve.m	(revision 1907)
@@ -11,5 +11,4 @@
 %   Examples:
 %      md=solve(md,'analysis_type','diagnostic','package','cielo');
-%      md=solve(md,'analysis_type','control','package','ice');
 %      md=solve(md,'analysis_type','thermal','sub_analysis_type','transient','package','ice');
 %      md=solve(md,'analysis_type','thermal','sub_analysis_type','steady','package','cielo');
@@ -67,7 +66,4 @@
 	md=prognostic(md);;
 
-elseif md.analysis_type==ControlAnalysisEnum,
-	md=control(md);
-
 elseif md.analysis_type==ThermalAnalysisEnum,
 	md=thermal(md);
@@ -91,7 +87,7 @@
 %Check result is consistent
 displaystring(md.debug,'%s\n','checking result consistency');
-if ~isresultconsistent(md,options.analysis_type),
-	disp('!! results not consistent correct the model !!') %it would be very cruel to put an error, it would kill the computed results (even if not consistent...)
-end
+%if ~isresultconsistent(md,options.analysis_type),
+%	disp('!! results not consistent correct the model !!') %it would be very cruel to put an error, it would kill the computed results (even if not consistent...)
+%end
 
 %convert analysis type to string finally
Index: /issm/trunk/src/m/solutions/cielo/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 1906)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 1907)
@@ -33,13 +33,25 @@
 	inputs=inputlist;
 	inputs=add(inputs,'velocity',models.dh.parameters.u_g,'doublevec',3,models.dh.parameters.numberofnodes);
+	if md.control_analysis,
+		inputs=add(inputs,'velocity_obs',models.dh.parameters.u_g_obs,'doublevec',2,models.dh.parameters.numberofnodes);
+	end
 
 	%compute solution
 	if ~models.dh.parameters.qmu_analysis,
-		%launch core of diagnostic solution.
-		results=diagnostic_core(models,inputs);
-	
-		%process results
-		if ~isstruct(md.results), md.results=struct(); end
-		md.results.diagnostic=processresults(models,results);
+		if md.control_analysis,
+			%launch core of control solution.
+			results=control_core(models,inputs);
+
+			%process results
+			if ~isstruct(md.results), md.results=struct(); end
+			md.results.control=processresults(models,results);
+		else,
+			%launch core of diagnostic solution.
+			results=diagnostic_core(models,inputs);
+
+			%process results
+			if ~isstruct(md.results), md.results=struct(); end
+			md.results.diagnostic=processresults(models,results);
+		end
 	else
 		%launch dakota driver for diagnostic core solution
Index: /issm/trunk/src/m/solutions/cielo/gradjcompute_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/gradjcompute_core.m	(revision 1906)
+++ /issm/trunk/src/m/solutions/cielo/gradjcompute_core.m	(revision 1907)
@@ -44,6 +44,6 @@
 %compute initial velocity from diagnostic_core (horiz+vertical)
 displaystring(debug,'%s','          compute initial velocity...');
-results_diag=diagnostic_core(models,inputs);
-u_g=results_diag.u_g;
+%results_diag=diagnostic_core(models,inputs);
+%u_g=results_diag.u_g;
 
 %Assign output