Index: /issm/trunk/src/m/solutions/cielo/modelsize.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/modelsize.m	(revision 939)
+++ /issm/trunk/src/m/solutions/cielo/modelsize.m	(revision 940)
@@ -14,5 +14,5 @@
 for i=1:length(modelsname);
 	if ~isempty(models.(modelsname{i}).nodesets),
-		dof=max(dof,models.(modelsname{i}).nodesets.fsize)
+		dof=max(dof,models.(modelsname{i}).nodesets.fsize);
 	end
 end
Index: /issm/trunk/src/m/solutions/cielo/thermal.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 939)
+++ /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 940)
@@ -24,8 +24,17 @@
 	inputs=add(inputs,'pressure',models.t.parameters.p_g,'doublevec',1,models.t.parameters.numberofnodes);
 	inputs=add(inputs,'dt',models.t.parameters.dt,'double');
-	models.t.parameters
 	
 	%call core
-	results=thermal_core(models,inputs,md.sub_analysis_type);
+	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
 
 	%plug onto the model
Index: /issm/trunk/src/m/solutions/cielo/thermal_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal_core.m	(revision 939)
+++ /issm/trunk/src/m/solutions/cielo/thermal_core.m	(revision 940)
@@ -1,7 +1,7 @@
-function results=thermal_core(models,inputs,sub_analysis_type)
+function results=thermal_core(models,inputs)
 %THERMAL_CORE - core of thermal solution
 %
 %   Usage:
-%      solution=thermal_core(models,inputs,sub_analysis_type)
+%      solution=thermal_core(models,inputs)
 
 %get FE model
@@ -9,20 +9,15 @@
 m_m=models.m;
 
-if strcmpi(sub_analysis_type,'steady'),
+if m_t.parameters.sub_analysis_type==steadyenum,
+
+	results.time=0;
 
 	displaystring(m_t.parameters.debug,'\n%s',['computing temperatures...']);
-	[t_g m_t.loads melting_offset]=thermal_core_nonlinear(m_t,inputs,'thermal','steady');
+	[results.t_g m_t.loads melting_offset]=thermal_core_nonlinear(m_t,inputs,'thermal','steady');
 
 	displaystring(m_t.parameters.debug,'\n%s',['computing melting...']);
 	inputs=add(inputs,'melting_offset',melting_offset,'double');
-	inputs=add(inputs,'temperature',t_g,'doublevec',1,m_t.parameters.numberofnodes);
-	m_g=diagnostic_core_linear(m_m,inputs,'melting','steady');
-
-	displaystring(m_t.parameters.debug,'\n%s',['load results...']);
-	results=struct('step',[],'time',[],'temperature',[],'melting',[]);
-	results.step=1;
-	results.time=0;
-	results.temperature=t_g;
-	results.melting=m_g*m_t.parameters.yts; %from m/s to m/a;
+	inputs=add(inputs,'temperature',results.t_g,'doublevec',1,m_t.parameters.numberofnodes);
+	results.m_g=diagnostic_core_linear(m_m,inputs,'melting','steady');
 
 else
@@ -34,31 +29,22 @@
 
 	%initialize temperature and melting
-	soln.t_g=t_g;
-	soln.m_g=m_g;
-	soln.time=0;
+	results.t_g=t_g;
+	results.m_g=m_g;
 
 	for n=1:nsteps, 
 
 		displaystring(m_t.parameters.debug,'\n%s%i/%i\n','time step: ',n,nsteps);
-		soln(n+1).time=n*m_t.parameters.dt;
+		results(n+1).time=n*m_t.parameters.dt;
 
 		displaystring(m_t.parameters.debug,'\n%s',['    computing temperatures...']);
-		inputs=add(inputs,'temperature',soln(n).t_g,'doublevec',1,m_t.parameters.numberofnodes);
-		[soln(n+1).t_g m_t.loads melting_offset]=thermal_core_nonlinear(m_t,inputs,'thermal','transient');
+		inputs=add(inputs,'temperature',results(n).t_g,'doublevec',1,m_t.parameters.numberofnodes);
+		[results(n+1).t_g m_t.loads melting_offset]=thermal_core_nonlinear(m_t,inputs,'thermal','transient');
 
 		disp('   computing melting...');
-		inputs=add(inputs,'temperature',soln(n+1).t_g,'doublevec',1,m_t.parameters.numberofnodes);
+		inputs=add(inputs,'temperature',results(n+1).t_g,'doublevec',1,m_t.parameters.numberofnodes);
 		inputs=add(inputs,'melting_offset',melting_offset,'double');
-		soln(n+1).m_g=diagnostic_core_linear(m_m,inputs,'melting','transient');
+		results(n+1).m_g=diagnostic_core_linear(m_m,inputs,'melting','transient');
 
 	end
 
-	%Wrap up
-	results=struct('step',{},'time',{},'temperature',{},'melting',{});
-	for n=1:nsteps+1,
-		results(n).temperature=soln(n).t_g;
-		results(n).melting=soln(n).m_g*m_t.parameters.yts; %in m/year
-		results(n).time=soln(n).time/m_t.parameters.yts;         %in years
-		results(n).step=n;
-	end
 end