Index: /issm/trunk/src/m/solutions/cielo/processresults.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/processresults.m	(revision 967)
+++ /issm/trunk/src/m/solutions/cielo/processresults.m	(revision 968)
@@ -12,5 +12,5 @@
 
 %recover models first
-if strcmpi(analysis_type,'diagnostic'),
+if (strcmpi(analysis_type,'diagnostic') | strcmpi(analysis_type,'transient')),
 	m_dh=models.dh;
 	m_ds=models.ds;
@@ -22,6 +22,6 @@
 	ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
 	isstokes=m_ds.parameters.isstokes;
-
-elseif strcmpi(analysis_type,'thermal'),
+end
+if (strcmpi(analysis_type,'thermal') | strcmpi(analysis_type,'transient')),
 	m_m=models.m;
 end
Index: /issm/trunk/src/m/solutions/cielo/prognostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/prognostic.m	(revision 967)
+++ /issm/trunk/src/m/solutions/cielo/prognostic.m	(revision 968)
@@ -25,11 +25,9 @@
 
 	displaystring(md.debug,'\n%s',['call computational core:']);
-	rawresults=prognostic_core(models,inputs,'prognostic','');
+	results=prognostic_core(models,inputs,'prognostic','');
 
 	displaystring(md.debug,'\n%s',['load results...']);
 	if ~isstruct(md.results), md.results=struct(); end
-	md.results.prognostic.step=1;
-	md.results.prognostic.time=0;
-	md.results.prognostic.thickness=rawresults.h_g;
+	md.results.prognostic=processresults(results,models,'prognostic');
 
 	%stop timing
Index: /issm/trunk/src/m/solutions/cielo/prognostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/prognostic_core.m	(revision 967)
+++ /issm/trunk/src/m/solutions/cielo/prognostic_core.m	(revision 968)
@@ -7,4 +7,6 @@
 	%get FE model
 	m=models.p;
+	results.time=0;
+	results.step=1;
 
 	displaystring(m.parameters.debug,'\n%s',['depth averaging velocity...']);
Index: /issm/trunk/src/m/solutions/cielo/thermal.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 967)
+++ /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 968)
@@ -26,19 +26,9 @@
 	
 	%call core
-	rawresults=thermal_core(models,inputs);
-
-	%process results
-	results=struct('step',{},'time',{},'temperature',{},'melting',{});
-	yts=models.t.parameters.yts;
-	for n=1:length(rawresults),
-		results(n).temperature=rawresults(n).t_g;
-		results(n).melting=rawresults(n).m_g*yts; %in m/year
-		results(n).time=rawresults(n).time/yts;      %in years
-		results(n).step=n;
-	end
+	results=thermal_core(models,inputs);
 
 	%plug onto the model
 	if ~isstruct(md.results), md.results=struct(); end
-	md.results.thermal=results;
+	md.results.thermal=processresults(results,models,'thermal');
 
 	%stop timing
Index: /issm/trunk/src/m/solutions/cielo/thermal_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal_core.m	(revision 967)
+++ /issm/trunk/src/m/solutions/cielo/thermal_core.m	(revision 968)
@@ -12,4 +12,5 @@
 
 	results.time=0;
+	results.step=1;
 
 	displaystring(m_t.parameters.debug,'\n%s',['computing temperatures...']);
@@ -29,11 +30,13 @@
 
 	%initialize temperature and melting
+	results.step=1;
+	results.time=0;
 	results.t_g=t_g;
 	results.m_g=m_g;
-	results.time=0;
 
 	for n=1:nsteps, 
 
 		displaystring(m_t.parameters.debug,'\n%s%i/%i\n','time step: ',n,nsteps);
+		results(n+1).step=n+1;
 		results(n+1).time=n*m_t.parameters.dt;
 
Index: /issm/trunk/src/m/solutions/cielo/transient3d.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/transient3d.m	(revision 967)
+++ /issm/trunk/src/m/solutions/cielo/transient3d.m	(revision 968)
@@ -35,23 +35,23 @@
 
 %initialize solution
-solution=struct('step',[],'time',[],'u_g',[],'p_g',[],'h_g',[],'s_g',[],'b_g',[]);
-solution.step=1;
-solution.time=0;
-solution.u_g=models.dh.parameters.u_g;
-%solution.u_g=models.dh.parameters.u_g(dofsetgen([1,2],3,length(models.dh.parameters.u_g)));
-solution.p_g=models.p.parameters.p_g;
-solution.h_g=models.p.parameters.h_g;
-solution.s_g=models.p.parameters.s_g;
-solution.b_g=models.p.parameters.b_g;
-solution.t_g=models.p.parameters.t_g;
-solution.m_g=models.p.parameters.m_g;
-solution.a_g=models.p.parameters.a_g;
+results=struct('step',[],'time',[],'u_g',[],'p_g',[],'h_g',[],'s_g',[],'b_g',[]);
+results.step=1;
+results.time=0;
+results.u_g=models.dh.parameters.u_g;
+%results.u_g=models.dh.parameters.u_g(dofsetgen([1,2],3,length(models.dh.parameters.u_g)));
+results.p_g=models.p.parameters.p_g;
+results.h_g=models.p.parameters.h_g;
+results.s_g=models.p.parameters.s_g;
+results.b_g=models.p.parameters.b_g;
+results.t_g=models.p.parameters.t_g;
+results.m_g=models.p.parameters.m_g;
+results.a_g=models.p.parameters.a_g;
 
 %initialize inputs
 displaystring(md.debug,'\n%s',['setup inputs...']);
 inputs=inputlist;
-inputs=add(inputs,'velocity',solution.u_g,'doublevec',3,models.p.parameters.numberofnodes);
-inputs=add(inputs,'melting',solution.m_g,'doublevec',1,models.p.parameters.numberofnodes);
-inputs=add(inputs,'accumulation',solution.a_g,'doublevec',1,models.p.parameters.numberofnodes);
+inputs=add(inputs,'velocity',results.u_g,'doublevec',3,models.p.parameters.numberofnodes);
+inputs=add(inputs,'melting',results.m_g,'doublevec',1,models.p.parameters.numberofnodes);
+inputs=add(inputs,'accumulation',results.a_g,'doublevec',1,models.p.parameters.numberofnodes);
 inputs=add(inputs,'dt',models.p.parameters.dt,'double');
 
@@ -70,35 +70,35 @@
 	displaystring(md.debug,'\n%s%g%s%g%s%g\n','time [yr]: ',time/yts,'    iteration number: ',n,'/',floor(finaltime/dt));
 
-	solution(n+1).step=n+1;
-	solution(n+1).time=time;
+	results(n+1).step=n+1;
+	results(n+1).time=time;
 
 	%update inputs
-	inputs=add(inputs,'thickness',solution(n).h_g,'doublevec',1,models.p.parameters.numberofnodes);
-	inputs=add(inputs,'surface',solution(n).s_g,'doublevec',1,models.p.parameters.numberofnodes);
-	inputs=add(inputs,'bed',solution(n).b_g,'doublevec',1,models.p.parameters.numberofnodes);
-	inputs=add(inputs,'velocity',solution(n).u_g,'doublevec',3,models.p.parameters.numberofnodes);
-	inputs=add(inputs,'pressure',solution(n).p_g,'doublevec',1,models.p.parameters.numberofnodes);
-	inputs=add(inputs,'temperature',solution(n).t_g,'doublevec',1,models.t.parameters.numberofnodes);
+	inputs=add(inputs,'thickness',results(n).h_g,'doublevec',1,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'surface',results(n).s_g,'doublevec',1,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'bed',results(n).b_g,'doublevec',1,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'velocity',results(n).u_g,'doublevec',3,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'pressure',results(n).p_g,'doublevec',1,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'temperature',results(n).t_g,'doublevec',1,models.t.parameters.numberofnodes);
 
 	%Deal with temperature first 
 	displaystring(md.debug,'\n%s',['    computing temperatures...']);
-	[solution(n+1).t_g models.t.loads melting_offset]=thermal_core_nonlinear(models.t,inputs,'thermal','transient');
-	inputs=add(inputs,'temperature',solution(n+1).t_g,'doublevec',1,models.t.parameters.numberofnodes);
+	[results(n+1).t_g models.t.loads melting_offset]=thermal_core_nonlinear(models.t,inputs,'thermal','transient');
+	inputs=add(inputs,'temperature',results(n+1).t_g,'doublevec',1,models.t.parameters.numberofnodes);
 	
 	displaystring(md.debug,'\n%s',['    computing melting...']);
 	inputs=add(inputs,'melting_offset',melting_offset,'double');
-	solution(n+1).m_g=diagnostic_core_linear(models.m,inputs,'melting','transient');
+	results(n+1).m_g=diagnostic_core_linear(models.m,inputs,'melting','transient');
 
 	%Compute depth averaged temperature
-	temperature_average=FieldDepthAverage(models.t.elements,models.t.nodes,models.t.loads,models.t.materials,solution(n+1).t_g,'temperature');
+	temperature_average=FieldDepthAverage(models.t.elements,models.t.nodes,models.t.loads,models.t.materials,results(n+1).t_g,'temperature');
 	inputs=add(inputs,'temperature_average',temperature_average,'doublevec',1,models.t.parameters.numberofnodes);
 
 	%Deal with velocities.
 	rawresults=diagnostic_core(models,inputs);
-	solution(n+1).u_g=rawresults.u_g; solution(n+1).p_g=rawresults.p_g;
+	results(n+1).u_g=rawresults.u_g; results(n+1).p_g=rawresults.p_g;
 
 	%compute new thickness
 	displaystring(md.debug,'\n%s',['    computing new thickness...']);
-	inputs=add(inputs,'velocity',solution(n+1).u_g,'doublevec',3,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'velocity',results(n+1).u_g,'doublevec',3,models.p.parameters.numberofnodes);
 	rawresults=prognostic_core(models,inputs,'prognostic','');
 	new_thickness=rawresults.h_g;
@@ -106,14 +106,14 @@
 	%update surface and bed using the new thickness
 	displaystring(md.debug,'\n%s',['    updating geometry...']);
-	[new_thickness,new_bed,new_surface]=UpdateGeometry(models.p.elements,models.p.nodes,models.p.loads,models.p.materials,models.p.parameters,new_thickness,solution(n).b_g,solution(n).s_g);
+	[new_thickness,new_bed,new_surface]=UpdateGeometry(models.p.elements,models.p.nodes,models.p.loads,models.p.materials,models.p.parameters,new_thickness,results(n).b_g,results(n).s_g);
 
-	%Record bed surface and thickness in the solution
-	solution(n+1).h_g=new_thickness;
-	solution(n+1).b_g=new_bed;
-	solution(n+1).s_g=new_surface;
+	%Record bed surface and thickness in the results
+	results(n+1).h_g=new_thickness;
+	results(n+1).b_g=new_bed;
+	results(n+1).s_g=new_surface;
 
 	%figure out if time stepping is good
 	%displaystring(md.debug,'   checking time stepping...'));
-	%[back,dt,time]=TimeStepping(md,solution,dt,time);
+	%[back,dt,time]=TimeStepping(md,results,dt,time);
 	%if back,
 	%	continue;
@@ -135,20 +135,5 @@
 end
 
-%Load results onto model
-results=struct('step',[],'time',[],'vx',[],'vy',[],'vel',[],'pressure',[],'thickness',[],'surface',[],'bed',[]);
-for i=1:length(solution),
-	results(i).step=solution(i).step;
-	results(i).time=solution(i).time/yts;
-	results(i).vx=solution(i).u_g(1:3:end)*yts;
-	results(i).vy=solution(i).u_g(2:3:end)*yts;
-	results(i).vz=solution(i).u_g(3:3:end)*yts;
-	results(i).vel=sqrt(solution(i).u_g(1:3:end).^2+solution(i).u_g(2:3:end).^2+solution(i).u_g(3:3:end).^2)*yts;
-	results(i).pressure=solution(i).p_g;
-	results(i).bed=solution(i).b_g;
-	results(i).surface=solution(i).s_g;
-	results(i).thickness=solution(i).h_g;
-	results(i).temperature=solution(i).t_g;
-	results(i).melting=solution(i).m_g;
-end
+%process results
 if ~isstruct(md.results), md.results=struct(); end
-md.results.transient=results;
+md.results.transient=processresults(results,models,'transient');