Index: /issm/trunk/src/m/solutions/cielo/control.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/control.m	(revision 979)
+++ /issm/trunk/src/m/solutions/cielo/control.m	(revision 980)
@@ -5,42 +5,43 @@
 	
 	%Build all models requested for control simulation
+	models.analysis_type='control'; 
 	md.analysis_type='control';
 	md.sub_analysis_type='';
-	m=CreateFemModel(md); 
+	models.dh=CreateFemModel(md); 
 
 	% figure out number of dof: just for information purposes.
-	md.dof=m.nodesets.fsize; %biggest dof number
+	md.dof=modelsize(models);
 
 	%initialize control parameters, gradients and observations
-	u_g_obs=m.parameters.u_g_obs;
-	u_g=m.parameters.u_g;
-	param_g=m.parameters.param_g;
-	grad_g=zeros(m.nodesets.gsize,1);
+	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(m.parameters);
+	options=ControlOptions(models.dh.parameters);
 
-	%initialize inputs, ie m.nparameters on which we invert.
+	%initialize inputs, ie models.dh.nparameters on which we invert.
 	inputs=inputlist;
-	inputs=add(inputs,m.parameters.control_type,param_g,'doublevec',2,m.parameters.numberofnodes);
-	inputs=add(inputs,'velocity',u_g,'doublevec',3,m.parameters.numberofnodes);
+	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);
 
-	for n=1:m.parameters.nsteps,
+	for n=1:models.dh.parameters.nsteps,
 		
 		%set options
-		options=optimset(options,'MaxFunEvals',m.parameters.maxiter(n));
+		options=optimset(options,'MaxFunEvals',models.dh.parameters.maxiter(n));
 
-		disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(m.parameters.nsteps)]));
+		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',m.parameters.fit(n),'double');
-		inputs=add(inputs,m.parameters.control_type,param_g,'doublevec',2,m.parameters.numberofnodes);
+		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
-		[m.elements,m.nodes,m.loads,m.materials]=UpdateFromInputs(m.elements,m.nodes,m.loads,m.materials,inputs);
+		[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(m,inputs,u_g_obs,md.analysis_type,md.sub_analysis_type);
-		inputs=add(inputs,'velocity',u_g,'doublevec',2,m.parameters.numberofnodes);
+		[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.');
 
@@ -54,22 +55,22 @@
 		
 		%visualize direction.
-		if m.parameters.plot
+		if models.dh.parameters.plot
 			plot_direction;
 		end
 		
 		disp('      optimizing along gradient direction...'); 
-		[search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m,inputs,param_g,u_g_obs,c(n).grad_g,n,md.analysis_type,md.sub_analysis_type);
+		[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*m.parameters.optscal(n)*c(n).grad_g;
+		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,m.parameters);
+		param_g=ControlConstrain(param_g,models.dh.parameters);
 		disp('      done.');
 
 		%visualize direction.
-		if m.parameters.plot,
+		if models.dh.parameters.plot,
 			plot_newdistribution;
 		end
@@ -77,5 +78,5 @@
 		%some temporary saving 
 		if(mod(n,5)==0),
-			solution=controlfinalsol(c,m,param_g,inputs,md.analysis_type,md.sub_analysis_type);
+			solution=controlfinalsol(c,models.dh,param_g,inputs,md.analysis_type,md.sub_analysis_type);
 			save temporary_control_results solution
 		end
@@ -86,9 +87,10 @@
 	%Create final solution
 	disp('      preparing final velocity solution...');
-	solution=controlfinalsol(c,m,param_g,inputs,md.analysis_type,md.sub_analysis_type);
+	results=controlfinalsol(c,models.dh,param_g,inputs,md.analysis_type,md.sub_analysis_type);
 	disp('      done.');
 	
 	disp('      load results ontol model...');
-	md=loadcontrolfinalsol(md,solution);
+	if ~isstruct(md.results), md.results=struct(); end
+	md.results.control=processresults(models,results);
 	disp('      done.');
 	
Index: /issm/trunk/src/m/solutions/cielo/controlfinalsol.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/controlfinalsol.m	(revision 979)
+++ /issm/trunk/src/m/solutions/cielo/controlfinalsol.m	(revision 980)
@@ -1,17 +1,11 @@
-function solution=controlfinalsol(c,m,p_g,inputs,analysis_type,sub_analysis_type);
+function results=controlfinalsol(c,m,param_g,inputs,analysis_type,sub_analysis_type);
+
+%initialize results
+results.time=0;
+results.step=1;
 
 %From parameters, build inputs for icediagnostic_core, using the final parameters
-inputs=add(inputs,m.parameters.control_type,p_g,'doublevec',2,m.parameters.numberofnodes);
-u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
-
-%Build partitioning vectors to recover solution
-indx=[1:2:m.nodesets.gsize];
-indy=[2:2:m.nodesets.gsize];
-
-%Recover velocity, and parameters, in the correct partitioning.
-vx=u_g(indx);
-vy=u_g(indy);
-vel=sqrt(vx.^2+vy.^2);
-parameter=p_g(indx);
+inputs=add(inputs,m.parameters.control_type,param_g,'doublevec',2,m.parameters.numberofnodes);
+results.u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
 
 %Recover misfit at each iteration of the control method 
@@ -20,9 +14,6 @@
 	J(i)=c(i).J;
 end
+results.J=J;
 
-%Store in solution
-solution.vx=vx;
-solution.vy=vy;
-solution.vel=vel;
-solution.J=J;
-solution.parameter=parameter;
+%store param_g
+results.param_g=param_g;
Index: /issm/trunk/src/m/solutions/cielo/processresults.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/processresults.m	(revision 979)
+++ /issm/trunk/src/m/solutions/cielo/processresults.m	(revision 980)
@@ -26,4 +26,13 @@
 	isstokes=m_ds.parameters.isstokes;
 end
+if strcmpi(analysis_type,'control'),
+	m_dh=models.dh;
+
+	%some flags needed
+	dim=m_dh.parameters.dim;
+	ishutter=m_dh.parameters.ishutter;
+	ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
+	isstokes=m_dh.parameters.isstokes;
+end
 if (strcmpi(analysis_type,'thermal') | strcmpi(analysis_type,'transient')),
 	m_m=models.m;
@@ -43,5 +52,5 @@
 
 			u_g=results(i).u_g;
-			yts=m_ds.parameters.yts;
+			yts=m_dh.parameters.yts;
 
 			if dim==2,