Index: /issm/trunk/src/m/solutions/cielo/control.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/control.m	(revision 1032)
+++ /issm/trunk/src/m/solutions/cielo/control.m	(revision 1033)
@@ -1,3 +1,8 @@
 function md=control(md)
+%CONTROL - launch control solution sequence
+%
+%   Usage:
+%      md=control(md)
+%
 
 	%timing
@@ -5,94 +10,30 @@
 	
 	%Build all models requested for control simulation
-	models.analysis_type='control'; 
-	md.analysis_type='control';
-	md.sub_analysis_type='';
-	models.dh=CreateFemModel(md); 
+	models.analysis_type='control'; %needed for processresults
+
+	displaystring(md.debug,'%s',['reading diagnostic horiz model data']);
+	md.analysis_type='control'; md.sub_analysis_type=''; models.dh=CreateFemModel(md); 
 
 	% figure out number of dof: just for information purposes.
 	md.dof=modelsize(models);
 
-	%initialize control parameters, gradients and observations
-	u_g_obs=models.dh.parameters.u_g_obs;
-	u_g=models.dh.parameters.u_g;
-	param_g=models.dh.parameters.param_g;
-	grad_g=zeros(models.dh.nodesets.gsize,1);
-		
-	%set optimization options.
-	options=ControlOptions(models.dh.parameters);
+	%initialize inputs
+	inputs=inputlist;
+	inputs=add(inputs,'velocity',models.dh.parameters.u_g,'doublevec',3,models.dh.parameters.numberofnodes);
 
-	%initialize inputs, ie models.dh.nparameters on which we invert.
-	inputs=inputlist;
-	inputs=add(inputs,models.dh.parameters.control_type,param_g,'doublevec',2,models.dh.parameters.numberofnodes);
-	inputs=add(inputs,'velocity',u_g,'doublevec',3,models.dh.parameters.numberofnodes);
+	%compute solution
+	if ~models.dh.parameters.qmu_analysis,
+		%launch core of control solution.
+		results=control_core(models,inputs);
 
-	for n=1:models.dh.parameters.nsteps,
-		
-		%set options
-		options=optimset(options,'MaxFunEvals',models.dh.parameters.maxiter(n));
-
-		disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(models.dh.parameters.nsteps)]));
-
-		%update inputs with new fit
-		inputs=add(inputs,'fit',models.dh.parameters.fit(n),'double');
-		inputs=add(inputs,models.dh.parameters.control_type,param_g,'doublevec',2,models.dh.parameters.numberofnodes);
-
-		%Update inputs in datasets
-		[models.dh.elements,models.dh.nodes,models.dh.loads,models.dh.materials]=UpdateFromInputs(models.dh.elements,models.dh.nodes,models.dh.loads,models.dh.materials,inputs);
-
-		disp('      computing gradJ...');
-		[u_g c(n).grad_g]=GradJCompute(models.dh,inputs,u_g_obs,md.analysis_type,md.sub_analysis_type);
-		inputs=add(inputs,'velocity',u_g,'doublevec',2,models.dh.parameters.numberofnodes);
-		disp('      done.');
-
-		disp('      normalizing directions...');
-		if n>=2,
-			c(n).grad_g=Orth(c(n).grad_g,c(n-1).grad_g);
-		else
-			c(n).grad_g=Orth(c(n).grad_g,{});
-		end
-		disp('      done.');
-		
-		%visualize direction.
-		if models.dh.parameters.plot
-			plot_direction;
-		end
-		
-		disp('      optimizing along gradient direction...'); 
-		[search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,models.dh,inputs,param_g,u_g_obs,c(n).grad_g,n,md.analysis_type,md.sub_analysis_type);
-		disp('      done.');
-
-		disp('      updating parameter using optimized search scalar...');
-		param_g=param_g+search_scalar*models.dh.parameters.optscal(n)*c(n).grad_g;
-		disp('      done.');
-
-		disp('      constraining the new distribution...');    
-		param_g=ControlConstrain(param_g,models.dh.parameters);
-		disp('      done.');
-
-		%visualize direction.
-		if models.dh.parameters.plot,
-			plot_newdistribution;
-		end
-
-		%some temporary saving 
-		if(mod(n,5)==0),
-			solution=controlfinalsol(c,models.dh,param_g,inputs,md.analysis_type,md.sub_analysis_type);
-			save temporary_control_results solution
-		end
-		disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);
-
+		%process results
+		if ~isstruct(md.results), md.results=struct(); end
+		md.results.control=processresults(models,results);
+	else
+		%launch dakota driver for diagnostic core solution
+		Qmu(models,inputs,models.dh.parameters);
 	end
 
-	%Create final solution
-	disp('      preparing final velocity solution...');
-	results=controlfinalsol(c,models.dh,param_g,inputs,md.analysis_type,md.sub_analysis_type);
-	disp('      done.');
-	
-	disp('      load results ontol model...');
-	if ~isstruct(md.results), md.results=struct(); end
-	md.results.control=processresults(models,results);
-	disp('      done.');
-	
+	%stop timing
 	t2=clock;
-	disp(['      overall time spent on control code ' num2str(etime(t2,t1)) ' seconds'])
+	displaystring(md.debug,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);
Index: /issm/trunk/src/m/solutions/cielo/control_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/control_core.m	(revision 1033)
+++ /issm/trunk/src/m/solutions/cielo/control_core.m	(revision 1033)
@@ -0,0 +1,75 @@
+function results=control_core(models,inputs)
+%CONTROL_CORE - compute the core inversion
+%
+%   Usage:
+%      results=control_core(models,inputs);
+%
+
+%recover models
+m_dh=models.dh;
+
+%recover parameters common to all solutions
+debug=m_dh.parameters.debug;
+dim=m_dh.parameters.dim;
+ishutter=m_dh.parameters.ishutter;
+ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
+isstokes=m_dh.parameters.isstokes;
+
+%initialize control parameters, gradients and observations
+u_g_obs=m_dh.parameters.u_g_obs;
+grad_g=zeros(m_dh.nodesets.gsize,1);
+param_g=models.dh.parameters.param_g;
+
+%set optimization options.
+options=ControlOptions(m_dh.parameters);
+
+for n=1:m_dh.parameters.nsteps,
+
+	%set options
+	options=optimset(options,'MaxFunEvals',m_dh.parameters.maxiter(n));
+
+	disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(m_dh.parameters.nsteps)]));
+
+	%update inputs with new fit
+	inputs=add(inputs,'fit',m_dh.parameters.fit(n),'double');
+	inputs=add(inputs,m_dh.parameters.control_type,param_g,'doublevec',2,m_dh.parameters.numberofnodes);
+
+	%Update inputs in datasets
+	[m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials]=UpdateFromInputs(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,inputs);
+
+	displaystring(debug,'\n%s',['      computing gradJ...']);
+	[u_g c(n).grad_g]=GradJCompute(m_dh,inputs,u_g_obs,'diagnostic','horiz');
+	inputs=add(inputs,'velocity',u_g,'doublevec',2,m_dh.parameters.numberofnodes);
+
+	displaystring(debug,'\n%s',['      normalizing directions...']);
+	if n>=2,
+		c(n).grad_g=Orth(c(n).grad_g,c(n-1).grad_g);
+	else
+		c(n).grad_g=Orth(c(n).grad_g,{});
+	end
+
+	%visualize direction.
+	if m_dh.parameters.plot
+		plot_direction;
+	end
+
+	displaystring(debug,'\n%s',['      optimizing along gradient direction...']);
+	[search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m_dh,inputs,param_g,u_g_obs,c(n).grad_g,n,'diagnostic','horiz');
+
+	displaystring(debug,'\n%s',['      updating parameter using optimized search scalar...']);
+	param_g=param_g+search_scalar*m_dh.parameters.optscal(n)*c(n).grad_g;
+
+	displaystring(debug,'\n%s',['      constraining the new distribution...']);
+	param_g=ControlConstrain(param_g,m_dh.parameters);
+
+	%visualize direction.
+	if m_dh.parameters.plot,
+		plot_newdistribution;
+	end
+
+	disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);
+end
+
+%generate output
+displaystring(debug,'\n%s',['      preparing final velocity solution...']);
+results=controlfinalsol(c,m_dh,param_g,inputs,'diagnostic','horiz');
Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 1032)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 1033)
@@ -3,5 +3,5 @@
 %
 %   Usage:
-%      results=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_sl,inputs);
+%      results=diagnostic_core(models,inputs);
 %