Index: /issm/trunk/src/m/solutions/cielo/modelsize.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/modelsize.m	(revision 935)
+++ /issm/trunk/src/m/solutions/cielo/modelsize.m	(revision 936)
@@ -3,25 +3,17 @@
 %
 %   Usage:
-%      dof=modelsize(m_dh,m_dv,m_ds,m_dhu,m_sl)
+%      dof=modelsize(models)
 
-%recover models
-m_dh=models.dh;
-m_dv=models.dv;
-m_ds=models.ds;
-m_dhu=models.dhu;
-m_sl=models.sl;
+%initialize dof
+dof=0;
 
-dof=0;
-if ~isempty(m_ds.nodesets),
-	dof=m_ds.nodesets.fsize; %biggest dof number
-end
-if ~isempty(m_dh.nodesets),
-	if dof<m_dh.nodesets.fsize,
-		dof=m_dh.nodesets.fsize; %biggest dof number
+%get all models
+modelsname=fieldnames(models);
+
+%get max dof
+for i=1:length(modelsname);
+	modelsname{i}
+	if ~isempty(models.(modelsname{i}).nodesets),
+		dof=max(dof,models.(modelsname{i}).nodesets.fsize)
 	end
 end
-if ~isempty(m_dhu.nodesets),
-	if dof<m_dhu.nodesets.fsize,
-		dof=m_dhu.nodesets.fsize;
-	end
-end
Index: /issm/trunk/src/m/solutions/cielo/prognostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/prognostic.m	(revision 935)
+++ /issm/trunk/src/m/solutions/cielo/prognostic.m	(revision 936)
@@ -10,20 +10,20 @@
 	%Build all models requested for diagnostic simulation
 	displaystring(md.debug,'%s',['reading prognostic model data']);
-	md.analysis_type='prognostic'; m_p=CreateFemModel(md);
+	md.analysis_type='prognostic'; models.p=CreateFemModel(md);
 
 	% figure out number of dof: just for information purposes.
-	md.dof=m_p.nodesets.fsize;
+	md.dof=models.p.nodesets.fsize;
 
 	%initialize inputs
 	displaystring(md.debug,'\n%s',['setup inputs...']);
 	inputs=inputlist;
-	inputs=add(inputs,'velocity',m_p.parameters.u_g,'doublevec',3,m_p.parameters.numberofnodes);
-	inputs=add(inputs,'thickness',m_p.parameters.h_g,'doublevec',1,m_p.parameters.numberofnodes);
-	inputs=add(inputs,'melting',m_p.parameters.m_g,'doublevec',1,m_p.parameters.numberofnodes);
-	inputs=add(inputs,'accumulation',m_p.parameters.a_g,'doublevec',1,m_p.parameters.numberofnodes);
-	inputs=add(inputs,'dt',m_p.parameters.dt,'double');
+	inputs=add(inputs,'velocity',models.p.parameters.u_g,'doublevec',3,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'thickness',models.p.parameters.h_g,'doublevec',1,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'dt',models.p.parameters.dt,'double');
 
 	displaystring(md.debug,'\n%s',['call computational core:']);
-	h_g=prognostic_core(m_p,inputs,'prognostic','');
+	rawresults=prognostic_core(models,inputs,'prognostic','');
 
 	displaystring(md.debug,'\n%s',['load results...']);
@@ -31,5 +31,5 @@
 	md.results.prognostic.step=1;
 	md.results.prognostic.time=0;
-	md.results.prognostic.thickness=h_g;
+	md.results.prognostic.thickness=rawresults.h_g;
 
 	%stop timing
Index: /issm/trunk/src/m/solutions/cielo/prognostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/prognostic_core.m	(revision 935)
+++ /issm/trunk/src/m/solutions/cielo/prognostic_core.m	(revision 936)
@@ -1,7 +1,10 @@
-function h_g=prognostic_core(m,inputs,analysis_type,sub_analysis_type)
+function results=prognostic_core(models,inputs,analysis_type,sub_analysis_type)
 %PROGNOSTIC_CORE - linear solution sequence
 %
 %   Usage:
 %      h_g=prognostic_core(m,inputs,analysis_type,sub_analysis_type)
+
+	%get FE model
+	m=models.p;
 
 	displaystring(m.parameters.debug,'\n%s',['depth averaging velocity...']);
@@ -12,8 +15,8 @@
 
 	displaystring(m.parameters.debug,'\n%s',['call computational core:']);
-	h_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type);
+	results.h_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type);
 
 	displaystring(m.parameters.debug,'\n%s',['extrude computed thickness on all layers:']);
-	h_g=FieldExtrude(m.elements,m.nodes,m.loads,m.materials,h_g,'thickness',0);
+	results.h_g=FieldExtrude(m.elements,m.nodes,m.loads,m.materials,results.h_g,'thickness',0);
 
 end %end function
Index: /issm/trunk/src/m/solutions/cielo/thermal.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 935)
+++ /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 936)
@@ -26,5 +26,5 @@
 	
 	%call core
-	results=thermal_core(models.t,models.m,inputs,md.sub_analysis_type);
+	results=thermal_core(models,inputs,md.sub_analysis_type);
 
 	%plug onto the model
Index: /issm/trunk/src/m/solutions/cielo/thermal_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal_core.m	(revision 935)
+++ /issm/trunk/src/m/solutions/cielo/thermal_core.m	(revision 936)
@@ -1,8 +1,11 @@
-function solution=thermal_core(m_t,m_m,inputs,sub_analysis_type)
+function results=thermal_core(models,inputs,sub_analysis_type)
 %THERMAL_CORE - core of thermal solution
 %
 %   Usage:
-%      solution=thermal_core(m_t,m_m,inputs,sub_analysis_type)
+%      solution=thermal_core(models,inputs,sub_analysis_type)
 
+%get FE model
+m_t=models.t;
+m_m=models.m;
 
 if strcmpi(sub_analysis_type,'steady'),
@@ -17,9 +20,9 @@
 
 	displaystring(m_t.parameters.debug,'\n%s',['load results...']);
-	solution=struct('step',[],'time',[],'temperature',[],'melting',[]);
-	solution.step=1;
-	solution.time=0;
-	solution.temperature=t_g;
-	solution.melting=m_g*m_t.parameters.yts; %from m/s to m/a;
+	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;
 
 else
@@ -52,10 +55,10 @@
 
 	%Wrap up
-	solution=struct('step',{},'time',{},'temperature',{},'melting',{});
+	results=struct('step',{},'time',{},'temperature',{},'melting',{});
 	for n=1:nsteps+1,
-		solution(n).temperature=soln(n).t_g;
-		solution(n).melting=soln(n).m_g*m_t.parameters.yts; %in m/year
-		solution(n).time=soln(n).time/m_t.parameters.yts;         %in years
-		solution(n).step=n;
+		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
Index: /issm/trunk/src/m/solutions/cielo/transient2d.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/transient2d.m	(revision 935)
+++ /issm/trunk/src/m/solutions/cielo/transient2d.m	(revision 936)
@@ -69,11 +69,12 @@
 
 	%Get horizontal solution. 
-	resultsdiag=diagnostic_core(models,inputs);
-	solution(n+1).u_g=resultsdiag.u_g; solution(n+1).p_g=resultsdiag.p_g;
+	rawresults=diagnostic_core(models,inputs);
+	solution(n+1).u_g=rawresults.u_g; solution(n+1).p_g=rawresults.p_g;
 
 	%compute new thickness
 	disp(sprintf('%s','   computing new thickness...'));
 	inputs=add(inputs,'velocity',solution(n+1).u_g,'doublevec',2,models.p.parameters.numberofnodes);
-	new_thickness=prognostic_core(models.p,inputs,'prognostic','');
+	rawresults=prognostic_core(models,inputs,'prognostic','');
+	new_thickness=rawresults.h_g;
 
 	%update surface and bed using the new thickness
Index: /issm/trunk/src/m/solutions/cielo/transient3d.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/transient3d.m	(revision 935)
+++ /issm/trunk/src/m/solutions/cielo/transient3d.m	(revision 936)
@@ -83,5 +83,5 @@
 	%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');
+	[solution(n+1).t_g models.t.loads melting_offset]=thermal_core_nonlinear(models,inputs,'thermal','transient');
 	inputs=add(inputs,'temperature',solution(n+1).t_g,'doublevec',1,models.t.parameters.numberofnodes);
 	
@@ -95,6 +95,6 @@
 
 	%Deal with velocities.
-	resultsdiag=diagnostic_core(models,inputs);
-	solution(n+1).u_g=resultsdiag.u_g; solution(n+1).p_g=resultsdiag.p_g;
+	rawresults=diagnostic_core(models,inputs);
+	solution(n+1).u_g=rawresults.u_g; solution(n+1).p_g=rawresults.p_g;
 
 	%compute new thickness