Index: /issm/trunk/src/m/solutions/cielo/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 921)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 922)
@@ -11,30 +11,30 @@
 	%Build all models requested for diagnostic simulation
 	displaystring(md.debug,'%s',['reading diagnostic horiz model data']);
-	md.analysis_type='diagnostic'; md.sub_analysis_type='horiz'; m_dh=CreateFemModel(md);
+	md.analysis_type='diagnostic'; md.sub_analysis_type='horiz'; models.dh=CreateFemModel(md);
 	
 	displaystring(md.debug,'\n%s',['reading diagnostic vert model data']);
-	md.analysis_type='diagnostic'; md.sub_analysis_type='vert'; m_dv=CreateFemModel(md);
+	md.analysis_type='diagnostic'; md.sub_analysis_type='vert'; models.dv=CreateFemModel(md);
 	
 	displaystring(md.debug,'\n%s',['reading diagnostic stokes model data']);
-	md.analysis_type='diagnostic'; md.sub_analysis_type='stokes'; m_ds=CreateFemModel(md);
+	md.analysis_type='diagnostic'; md.sub_analysis_type='stokes'; models.ds=CreateFemModel(md);
 
 	displaystring(md.debug,'\n%s',['reading diagnostic hutter model data']);
-	md.analysis_type='diagnostic'; md.sub_analysis_type='hutter'; m_dhu=CreateFemModel(md);
+	md.analysis_type='diagnostic'; md.sub_analysis_type='hutter'; models.dhu=CreateFemModel(md);
 	
 	displaystring(md.debug,'\n%s',['reading surface and bed slope computation model data']);
-	md.analysis_type='slope_compute'; md.sub_analysis_type=''; m_sl=CreateFemModel(md);
+	md.analysis_type='slope_compute'; md.sub_analysis_type=''; models.sl=CreateFemModel(md);
 	
 	% figure out number of dof: just for information purposes.
-	md.dof=modelsize(m_dh,m_dv,m_ds,m_dhu,m_sl);
+	md.dof=modelsize(models);
 
 	%initialize inputs
 	inputs=inputlist;
-	inputs=add(inputs,'velocity',m_dh.parameters.u_g,'doublevec',3,m_dh.parameters.numberofnodes);
+	inputs=add(inputs,'velocity',models.dh.parameters.u_g,'doublevec',3,models.dh.parameters.numberofnodes);
 
 	%compute solution
-	[u_g,p_g]=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_sl,inputs);
+	results=diagnostic_core(models,inputs);
 
-	%Load results onto model
-	md=loadresults(md,u_g,p_g,m_dh,m_ds,m_dhu);
+	%load results onto model
+	md=loadresults(md,models,results);
 
 	%stop timing
Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 921)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 922)
@@ -1,8 +1,15 @@
-function [u_g p_g]=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_sl,inputs);
+function results=diagnostic_core(models,inputs);
 %DIAGNOSTIC_CORE - compute the core velocity field 
 %
 %   Usage:
-%      u_g=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_sl,inputs);
+%      results=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_sl,inputs);
 %
+
+%recover models
+m_dh=models.dh;
+m_dv=models.dv;
+m_ds=models.ds;
+m_dhu=models.dhu;
+m_sl=models.sl;
 
 %recover parameters common to all solutions
@@ -102,2 +109,6 @@
 	end
 end
+
+%load onto results
+results.u_g=u_g;
+results.p_g=p_g;
Index: /issm/trunk/src/m/solutions/cielo/loadcontrolfinalsol.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/loadcontrolfinalsol.m	(revision 921)
+++ /issm/trunk/src/m/solutions/cielo/loadcontrolfinalsol.m	(revision 922)
@@ -1,5 +1,9 @@
 function md=loadfinalcontrolsol(md,solution);
 
-if ~isstruct(md.results), md.results=struct(); end
+if ~isstruct(md.results),
+	if isnan(md.results), 
+		md.results=struct(); 
+	end
+end
 md.results.control.vx=solution.vx*md.yts;
 md.results.control.vy=solution.vy*md.yts;
Index: /issm/trunk/src/m/solutions/cielo/loadresults.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/loadresults.m	(revision 921)
+++ /issm/trunk/src/m/solutions/cielo/loadresults.m	(revision 922)
@@ -1,3 +1,3 @@
-function md=loadresults(md,u_g,p_g,m_dh,m_ds,m_dhu);
+function md=loadresults(md,models,results);
 %LOADRESULTS - load results onto model
 %
@@ -6,6 +6,12 @@
 %
 %   Usage:
-%      md=Loadresults(md,u_g,p_g,m_dh,m_ds,m_dhu)
+%      md=Loadresults(md,results,m_dh,m_ds,m_dhu)
 
+%recover models
+m_dh=models.dh;
+m_ds=models.ds;
+m_dhu=models.dhu;
+
+%recover parameters
 debug=m_dhu.parameters.debug;
 dim=m_dhu.parameters.dim;
@@ -13,4 +19,8 @@
 ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
 isstokes=m_ds.parameters.isstokes;
+
+%recover results
+u_g=results.u_g;
+p_g=results.p_g;
 
 if ~isstruct(md.results), md.results=struct(); end 
Index: /issm/trunk/src/m/solutions/cielo/modelsize.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/modelsize.m	(revision 921)
+++ /issm/trunk/src/m/solutions/cielo/modelsize.m	(revision 922)
@@ -1,7 +1,14 @@
-function dof=modelsize(m_dh,m_dv,m_ds,m_dhu,m_sl);
+function dof=modelsize(models)
 %DOF - return the maximum number of degrees of freedom of the model
 %
 %   Usage:
 %      dof=modelsize(m_dh,m_dv,m_ds,m_dhu,m_sl)
+
+%recover models
+m_dh=models.dh;
+m_dv=models.dv;
+m_ds=models.ds;
+m_dhu=models.dhu;
+m_sl=models.sl;
 
 dof=0;