Index: /issm/trunk/src/m/solutions/cielo/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 966)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 967)
@@ -38,5 +38,6 @@
 	
 		%process results
-		md.results.diagnostic=processresults(models,results);
+		if ~isstruct(md.results), md.results=struct(); end
+		md.results.diagnostic=processresults(results,models,'diagnostic');
 	else
 		%launch dakota driver for diagnostic core solution
Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 966)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 967)
@@ -111,4 +111,6 @@
 
 %load onto results
+results.step=1;
+results.time=0;
 results.u_g=u_g;
 results.p_g=p_g;
Index: /issm/trunk/src/m/solutions/cielo/processresults.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/processresults.m	(revision 966)
+++ /issm/trunk/src/m/solutions/cielo/processresults.m	(revision 967)
@@ -1,3 +1,3 @@
-function newresults=processresults(models,results);
+function newresults=processresults(results,models,analysis_type);
 %PROCESSRESULTS - process results from core solution
 %
@@ -6,47 +6,119 @@
 %
 %   Usage:
-%      results=processresults(md,results,m_dh,m_ds,m_dhu)
+%      results=processresults(results,models,analysis_type)
 
-%recover models
-m_dh=models.dh;
-m_ds=models.ds;
-m_dhu=models.dhu;
+%initialize output
+newresults=struct();
 
-%recover parameters
-debug=m_dhu.parameters.debug;
-dim=m_dhu.parameters.dim;
-ishutter=m_dhu.parameters.ishutter;
-ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
-isstokes=m_ds.parameters.isstokes;
-yts=m_ds.parameters.yts;
+%recover models first
+if strcmpi(analysis_type,'diagnostic'),
+	m_dh=models.dh;
+	m_ds=models.ds;
+	m_dhu=models.dhu;
 
-%recover results
-u_g=results.u_g;
-p_g=results.p_g;
+	%some flags needed
+	dim=m_dhu.parameters.dim;
+	ishutter=m_dhu.parameters.ishutter;
+	ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
+	isstokes=m_ds.parameters.isstokes;
 
-newresults=struct();
-if dim==2,
-	newresults.step=1;
-	newresults.time=1;
-	newresults.vx=u_g(1:2:end)*yts;
-	newresults.vy=u_g(2:2:end)*yts;
-	newresults.vel=sqrt(newresults.vx.^2+newresults.vy.^2);
-	newresults.pressure=p_g;
+elseif strcmpi(analysis_type,'thermal'),
+	m_m=models.m;
+end
 
-else
-	
-	newresults.step=1;
-	newresults.time=1;
-	if isstokes,
-		newresults.vx=u_g(1:4:end)*yts;
-		newresults.vy=u_g(2:4:end)*yts;
-		newresults.vz=u_g(3:4:end)*yts;
-		newresults.pressure=p_g;
-	else
-		newresults.vx=u_g(1:3:end)*yts;
-		newresults.vy=u_g(2:3:end)*yts;
-		newresults.vz=u_g(3:3:end)*yts;
-		newresults.pressure=p_g;
+
+%go through results
+for i=1:length(results),
+
+	%get field names of results
+	resultsname=fieldnames(results(i));
+
+	%go through all field
+	for j=1:length(resultsname),
+
+		if strcmpi(resultsname{j},'u_g'),
+
+			u_g=results(i).u_g;
+			yts=m_ds.parameters.yts;
+
+			if dim==2,
+				newresults(i).step=results(i).step;
+				newresults(i).time=results(i).time;
+				newresults(i).vx=u_g(1:2:end)*yts;
+				newresults(i).vy=u_g(2:2:end)*yts;
+				newresults(i).vel=sqrt(newresults(i).vx.^2+newresults(i).vy.^2);
+
+			else
+				newresults(i).step=results(i).step;
+				newresults(i).time=results(i).time;
+				if isstokes,
+					newresults(i).vx=u_g(1:4:end)*yts;
+					newresults(i).vy=u_g(2:4:end)*yts;
+					newresults(i).vz=u_g(3:4:end)*yts;
+				else
+					newresults(i).vx=u_g(1:3:end)*yts;
+					newresults(i).vy=u_g(2:3:end)*yts;
+					newresults(i).vz=u_g(3:3:end)*yts;
+				end
+				newresults(i).vel=sqrt(newresults(i).vx.^2+newresults(i).vy.^2+newresults(i).vz.^2);
+			end
+
+		elseif strcmpi(resultsname{j},'p_g'),
+
+			p_g=results(i).p_g;
+			newresults(i).step=results(i).step;
+			newresults(i).time=results(i).time;
+			newresults(i).pressure=p_g;
+
+		elseif strcmpi(resultsname{j},'t_g'),
+
+			t_g=results(i).t_g;
+			newresults(i).step=results(i).step;
+			newresults(i).time=results(i).time;
+			newresults(i).temperature=t_g;
+
+		elseif strcmpi(resultsname{j},'m_g'),
+
+			m_g=results(i).m_g;
+			yts=m_m.parameters.yts;
+			newresults(i).step=results(i).step;
+			newresults(i).time=results(i).time;
+			newresults(i).melting=m_g*yts;
+
+		elseif strcmpi(resultsname{j},'h_g'),
+
+			h_g=results(i).h_g;
+			newresults(i).step=results(i).step;
+			newresults(i).time=results(i).time;
+			newresults(i).thickness=h_g;
+
+		elseif strcmpi(resultsname{j},'s_g'),
+
+			s_g=results(i).s_g;
+			newresults(i).step=results(i).step;
+			newresults(i).time=results(i).time;
+			newresults(i).surface=s_g;
+
+		elseif strcmpi(resultsname{j},'b_g'),
+
+			b_g=results(i).b_g;
+			newresults(i).step=results(i).step;
+			newresults(i).time=results(i).time;
+			newresults(i).bed=b_g;
+
+		elseif strcmpi(resultsname{j},'param_g'),
+
+			param_g=results(i).param_g;
+			newresults(i).step=results(i).step;
+			newresults(i).time=results(i).time;
+			newresults(i).parameter=param_g(1:2:end);
+
+		else
+
+			newresults(i).step=results(i).step;
+			newresults(i).time=results(i).time;
+			newresults(i).(resultsname{j})=results(i).(resultsname{j});
+
+		end
 	end
-	newresults.vel=sqrt(newresults.vx.^2+newresults.vy.^2+newresults.vz.^2);
 end