Index: /issm/trunk/src/m/solutions/cielo/thermal.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 730)
+++ /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 731)
@@ -25,60 +25,9 @@
 	inputs=add(inputs,'dt',m_t.parameters.dt,'double');
 	
-	if strcmpi(md.sub_analysis_type,'steady'),
-	
-		displaystring(md.debug,'\n%s',['computing temperatures...']);
-		[t_g m_t.loads melting_offset]=thermal_core(m_t,inputs,'thermal','steady');
-		
-		displaystring(md.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(md.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*md.yts; %from m/s to m/a;
-		md.results.thermal=solution;
+	%call core
+	results=thermal_core(m_t,m_m,inputs,md.sub_analysis_type);
 
-	else
-
-		%initialize temperature and melting
-		t_g=m_t.parameters.t_g;
-		m_g=m_m.parameters.m_g;
-		nsteps=m_t.parameters.ndt/m_t.parameters.dt;
-
-		%initialize temperature and melting
-		soln.t_g=t_g;
-		soln.m_g=m_g;
-		soln.time=0;
-
-		for n=1:nsteps, 
-
-			displaystring(md.debug,'\n%s%i/%i\n','time step: ',n,nsteps);
-			soln(n+1).time=n*m_t.parameters.dt;
-			
-			displaystring(md.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(m_t,inputs,'thermal','transient');
-			
-			disp('   computing melting...');
-			inputs=add(inputs,'temperature',soln(n).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');
-			
-		end
-		
-		%Wrap up
-		solution=struct('step',{},'time',{},'temperature',{},'melting',{});
-		for n=1:nsteps+1,
-			solution(n).temperature=soln(n).t_g;
-			solution(n).melting=soln(n).m_g*md.yts; %in m/year
-			solution(n).time=soln(n).time/md.yts;         %in years
-			solution(n).step=n;
-		end
-		md.results.thermal=solution;
-	end
+	%plug onto the model
+	md.results.thermal=results;
 
 	%stop timing
Index: /issm/trunk/src/m/solutions/cielo/thermal_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal_core.m	(revision 730)
+++ /issm/trunk/src/m/solutions/cielo/thermal_core.m	(revision 731)
@@ -1,76 +1,61 @@
-function [t_g ,loads, melting_offset]=thermal_core(m,inputs,analysis_type,sub_analysis_type)
-%THERMAL_CORE - core of thermal solution sequence.
-%   model is return together with temperature
-%
+function solution=thermal_core(m_t,m_m,inputs,sub_analysis_type)
+%THERMAL_CORE - core of thermal solution
 %
 %   Usage:
-%      [t_g ,loads, melting_offset]=thermal_core(m,inputs,analysis_type,sub_analysis_type);
-%    
-%      
+%      solution=thermal_core(m_t,m_m,inputs,sub_analysis_type)
 
-	count=1;
-	converged=0;
 
-%   we are going to return the loads, make them a variable of this routine
-	loads=m.loads;
+if strcmpi(sub_analysis_type,'steady'),
 
-	%stiffness and load generation only:
-	m.parameters.kflag=1; m.parameters.pflag=1;
+	displaystring(m_t.parameters.debug,'\n%s',['computing temperatures...']);
+	[t_g m_t.loads melting_offset]=thermal_core_nonlinear(m_t,inputs,'thermal','steady');
 
-	displaystring(m.parameters.debug,'\n%s',['   starting direct shooting method']);
-	
-	while(~converged),
+	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');
 
-		%Update inputs in datasets
-		[m.elements,m.nodes, loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,inputs);
+	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;
 
-		%system matrices 
-		if ~m.parameters.lowmem 
-			if count==1
-				displaystring(m.parameters.debug,'%s',['   system matrices']);
-				[K_gg_nopenalty, p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-			end
-			displaystring(m.parameters.debug,'%s',['   penalty system matrices']);
-			[K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-		else
-			displaystring(m.parameters.debug,'%s',['   system matrices']);
-			[K_gg , p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-			displaystring(m.parameters.debug,'%s',['   penalty system matrices']);
-			[K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-		end
+else
 
-		%Reduce tangent matrix from g size to f size
-		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); 
+	%initialize temperature and melting
+	t_g=m_t.parameters.t_g;
+	m_g=m_m.parameters.m_g;
+	nsteps=m_t.parameters.ndt/m_t.parameters.dt;
 
-		%Reduce load from g size to f size
-		[p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets);
+	%initialize temperature and melting
+	soln.t_g=t_g;
+	soln.m_g=m_g;
+	soln.time=0;
 
-		%Solve	
-		displaystring(m.parameters.debug,'%s%g','   condition number of stiffness matrix: ',condest(K_ff));
-		[t_f]=Solver(K_ff,p_f,[],m.parameters);
+	for n=1:nsteps, 
 
-		%Merge back to g set
-		displaystring(m.parameters.debug,'%s',['   merging solution back to g set']);
-		[t_g]= Mergesolutionfromftog( t_f, m.Gmn, m.ys, m.nodesets ); 
+		displaystring(m_t.parameters.debug,'\n%s%i/%i\n','time step: ',n,nsteps);
+		soln(n+1).time=n*m_t.parameters.dt;
 
-		%Update inputs in datasets
-		inputs=add(inputs,'temperature',t_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
-		displaystring(m.parameters.debug,'%s',['   update inputs']);
-		[m.elements,m.nodes, loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,inputs);
-	
-		%penalty constraints
-		displaystring(m.parameters.debug,'%s',['   penalty constraints']);
-		[loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
-	
-		if ~converged,
-			displaystring(m.parameters.debug,'%s%i','   #unstable constraints ',num_unstable_constraints);
-			
-			if num_unstable_constraints<=m.parameters.min_thermal_constraints,
-				converged=1;
-			end
-		end
+		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');
 
-		count=count+1;
+		disp('   computing melting...');
+		inputs=add(inputs,'temperature',soln(n).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');
+
 	end
 
+	%Wrap up
+	solution=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;
+	end
 end
Index: /issm/trunk/src/m/solutions/cielo/thermal_core_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal_core_nonlinear.m	(revision 731)
+++ /issm/trunk/src/m/solutions/cielo/thermal_core_nonlinear.m	(revision 731)
@@ -0,0 +1,76 @@
+function [t_g ,loads, melting_offset]=thermal_core_nonlinear(m,inputs,analysis_type,sub_analysis_type)
+%THERMAL_CORE_NONLINEAR - core of thermal solution sequence.
+%   model is return together with temperature
+%
+%
+%   Usage:
+%      [t_g ,loads, melting_offset]=thermal_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
+%    
+%      
+
+	count=1;
+	converged=0;
+
+%   we are going to return the loads, make them a variable of this routine
+	loads=m.loads;
+
+	%stiffness and load generation only:
+	m.parameters.kflag=1; m.parameters.pflag=1;
+
+	displaystring(m.parameters.debug,'\n%s',['   starting direct shooting method']);
+	
+	while(~converged),
+
+		%Update inputs in datasets
+		[m.elements,m.nodes, loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,inputs);
+
+		%system matrices 
+		if ~m.parameters.lowmem 
+			if count==1
+				displaystring(m.parameters.debug,'%s',['   system matrices']);
+				[K_gg_nopenalty, p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+			end
+			displaystring(m.parameters.debug,'%s',['   penalty system matrices']);
+			[K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+		else
+			displaystring(m.parameters.debug,'%s',['   system matrices']);
+			[K_gg , p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+			displaystring(m.parameters.debug,'%s',['   penalty system matrices']);
+			[K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+		end
+
+		%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	
+		displaystring(m.parameters.debug,'%s%g','   condition number of stiffness matrix: ',condest(K_ff));
+		[t_f]=Solver(K_ff,p_f,[],m.parameters);
+
+		%Merge back to g set
+		displaystring(m.parameters.debug,'%s',['   merging solution back to g set']);
+		[t_g]= Mergesolutionfromftog( t_f, m.Gmn, m.ys, m.nodesets ); 
+
+		%Update inputs in datasets
+		inputs=add(inputs,'temperature',t_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
+		displaystring(m.parameters.debug,'%s',['   update inputs']);
+		[m.elements,m.nodes, loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,inputs);
+	
+		%penalty constraints
+		displaystring(m.parameters.debug,'%s',['   penalty constraints']);
+		[loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+	
+		if ~converged,
+			displaystring(m.parameters.debug,'%s%i','   #unstable constraints ',num_unstable_constraints);
+			
+			if num_unstable_constraints<=m.parameters.min_thermal_constraints,
+				converged=1;
+			end
+		end
+
+		count=count+1;
+	end
+
+end