Index: /issm/trunk/src/m/solutions/cielo/prognostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/prognostic.m	(revision 726)
+++ /issm/trunk/src/m/solutions/cielo/prognostic.m	(revision 727)
@@ -12,9 +12,4 @@
 	md.analysis_type='prognostic'; m_p=CreateFemModel(md);
 
-	%Take only the first two dofs of m_p.parameters.u_g
-	u_g=m_p.parameters.u_g(dofsetgen([1,2],3,m_p.parameters.numberofnodes*3));
-	displaystring(md.debug,'\n%s',['depth averaging velocity...']);
-	u_g=VelocityDepthAverage(m_p.elements,m_p.nodes,m_p.loads,m_p.materials,u_g);
-
 	% figure out number of dof: just for information purposes.
 	md.dof=m_p.nodesets.fsize;
@@ -23,5 +18,5 @@
 	displaystring(md.debug,'\n%s',['setup inputs...']);
 	inputs=inputlist;
-	inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m_p.parameters.numberofnodes);
+	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);
@@ -29,10 +24,7 @@
 	inputs=add(inputs,'dt',m_p.parameters.dt,'double');
 
-	displaystring(md.debug,'\n%s',['computing new thickness...']);
+	displaystring(md.debug,'\n%s',['call computational core:']);
 	h_g=prognostic_core(m_p,inputs,'prognostic','');
-	
-	displaystring(md.debug,'\n%s',['extrude new thickness...']);
-	h_g=ThicknessExtrude(m_p.elements,m_p.nodes,m_p.loads,m_p.materials,h_g);
-	
+
 	displaystring(md.debug,'\n%s',['load results...']);
 	md.results.prognostic.step=1;
Index: /issm/trunk/src/m/solutions/cielo/prognostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/prognostic_core.m	(revision 726)
+++ /issm/trunk/src/m/solutions/cielo/prognostic_core.m	(revision 727)
@@ -1,29 +1,19 @@
-function u_g=prognostic_core(m,inputs,analysis_type,sub_analysis_type)
+function h_g=prognostic_core(m,inputs,analysis_type,sub_analysis_type)
 %PROGNOSTIC_CORE - linear solution sequence
 %
 %   Usage:
-%      u_g=prognostic_core(m,inputs,analysis_type,sub_analysis_type)
+%      h_g=prognostic_core(m,inputs,analysis_type,sub_analysis_type)
 
-	%stiffness and load generation only:
-	m.parameters.kflag=1; m.parameters.pflag=1;
+	displaystring(m.parameters.debug,'\n%s',['depth averaging velocity...']);
+	%Take only the first two dofs of m.parameters.u_g
+	u_g=get(inputs,'velocity',[1 1 0 0]);
+	u_g=VelocityDepthAverage(m.elements,m.nodes,m.loads,m.materials,u_g);
+	inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.numberofnodes);
 
-	%Update inputs in datasets
-	[m.elements,m.nodes, m.loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, m.loads,m.materials,inputs);
-	
-	%system matrices
-	[K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-	[K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-	
-	%Reduce tangent matrix from g size to f size
-	[K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); 
-	
-	%Reduce load from g size to f size
-	[p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets);
-	
-	%Solve	
-	[u_f]=Solver(K_ff,p_f,[],m.parameters);
-	
-	%Merge back to g set
-	[u_g]= Mergesolutionfromftog( u_f, m.Gmn, m.ys, m.nodesets ); 
+	displaystring(m.parameters.debug,'\n%s',['call computational core:']);
+	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=ThicknessExtrude(m.elements,m.nodes,m.loads,m.materials,h_g);
 
 end %end function
Index: /issm/trunk/src/m/solutions/cielo/thermal.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 726)
+++ /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 727)
@@ -36,8 +36,10 @@
 		
 		displaystring(md.debug,'\n%s',['load results...']);
-		md.results.thermal.step=1;
-		md.results.thermal.time=0;
-		md.results.thermal.temperature=t_g;
-		md.results.thermal.melting=m_g*md.yts; %from m/s to m/a
+		solution=struct('step',[],'time',[],'temperature',[],'melting',[]);
+		solution.step=1;
+		solution.time=0;
+		solution.temperature=t_g;
+		solution.melting=m_g*md.yts; %from m/s to m/a;
+		md.results.thermal=solution;
 
 	else
@@ -70,5 +72,5 @@
 		
 		%Wrap up
-		solution=struct('time',{},'temperature',{},'melting',{});
+		solution=struct('step',{},'time',{},'temperature',{},'melting',{});
 		for n=1:nsteps+1,
 			solution(n).temperature=soln(n).t_g;