Index: /issm/trunk/src/m/solutions/cielo/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 162)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 163)
@@ -9,5 +9,4 @@
 
 	%Build all models requested for diagnostic simulation
-	%remark, partitions are all identical.
 	md.analysis_type='diagnostic_horiz'; m_dh=CreateFemModel(md);
 
Index: /issm/trunk/src/m/solutions/cielo/thermal.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 163)
+++ /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 163)
@@ -0,0 +1,82 @@
+function md=thermal(md)
+%THERMAL - thermal solution sequence: steady state and transient
+%
+%   Usage:
+%      md=thermal(md)
+
+	%timing
+	t1=clock;
+	
+	%md.analysis_type is set to thermalsteady or thermaltransient. save it.
+	solutiontype=md.analysis_type;
+
+	%Build all models requested for thermal and melting  simulation
+	md.analysis_type=solutiontype; m_t=CreateFemModel(md); 
+	md.analysis_type='melting' m_m=CreateFEMModel(md);
+
+	% figure out number of dof: just for information purposes.
+	md.dof=m_t.nodesets.fsize; %biggest dof number
+
+
+	if strcmpi(solutiontype,'thermalsteady'),
+
+		%initialize velocity and pressure 
+		u_g=m_t.parameters.u_g;
+		p_g=m_t.parameters.p_g;
+
+		inputs=struct('pressure',p_g,'velocity',u_g);
+
+		disp(sprintf('\n%s\n','computing temperature...'));
+		[t_g,m_t]=thermal_core(m_t,inputs);
+
+		disp(sprintf('\n%s\n','computing melting...'));
+		inputs=struct('pressure',p_g,'temperature',t_g);
+		
+		melting_g=melting_core(m_m,inputs);
+
+		%load results onto model
+		md.temperature=t_g;
+		md.melting=melting_g*md.yts; %from m/s to m/a
+
+	else
+
+		%initialize velocity, pressure  and temperature
+		u_g=m.parameters.u_g;
+		p_g=m.parameters.p_g;
+		t_g=m.parameters.t_g;
+		melting_g=m.parameters.melting_g;
+
+		nsteps=m_t.parameters.ndt/m_t.parameters.dt;
+
+		%initialize temperature and melting
+		soln.t_g=t_g;
+		soln.melting_g=melting_g;
+
+		for n=1:nsteps, 
+
+			disp(sprintf('\n%s%i/%i\n','time step: ',n,md.ndt/md.dt));
+			disp('   computing temperature...');
+		
+			inputs=struct('pressure',p_g,'temperature',soln(n).t_g,'velocity',u_g,'dt',m_t.parameters.dt);
+			[soln(n+1).t_g,m_t]=thermal_core(m_t,inputs);
+			
+			disp('   computing melting...');
+			inputs=struct('pressure',p_g,'temperature',soln(n).t_g,'dt',m_t.parameters.dt);
+			soln(n+1).melting_g=melting_core(m_m,inputs);
+			
+		end
+		
+		%Wrap up
+		solution_temperature=struct('temperature',{});
+		for n=1:nsteps+1,
+			solution_temperature(n).temperature=soln(n).t_g;
+		end
+		md.temperature=solution_temperature;
+		
+		solution_melting=struct('melting',{});
+		for n=1:nsteps+1,
+			solution_melting(n).melting=soln(n).melting_g*md.yts; %in m/a
+		end
+		md.melting=solution_melting;
+	end
+	
Index: /issm/trunk/src/m/solutions/cielo/thermal_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal_core.m	(revision 163)
+++ /issm/trunk/src/m/solutions/cielo/thermal_core.m	(revision 163)
@@ -0,0 +1,59 @@
+function [t_g ,m]=thermal_core(m,inputs)
+%THERMAL_CORE - core of thermal solution sequence.
+%   model is return together with temperature
+%
+%
+%   Usage:
+%      [t_g,m]=thermal_core(m,inputs)
+%    
+%      
+
+	count=1;
+	converged=0;
+
+	%stiffness and load generation only:
+	m.parameters.kflag=1; m.parameters.pflag=1;
+
+	disp(sprintf('%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 
+		[K_gg_ , p_g_]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs);
+
+		%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	
+		if m.parameters.debug,
+			disp(sprintf('%s%g','      condition number of stiffness matrix: ',condest(K_ff)));
+		end
+		[t_f]=Solver(K_ff,p_f,[],m.parameters);
+
+		%Merge back to g set
+		[t_g]= Mergesolutionfromftog( t_f, m.Gmn, m.ys, m.nodesets ); 
+
+		%Update inputs in datasets
+		inputs.temperature=t_g;
+		[m.elements,m.nodes, loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,inputs);
+	
+		%penalty constraints
+		[m.loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, m.loads, m.materials,m.parameters,inputs);
+	
+		if ~converged,
+			if m.parameters.debug, disp(sprintf('   %s %i','#unstable constraints ',num_unstable_constraints));end;
+			
+			if num_unstable_constraints<thermalparams.min_thermal_constraints,
+				converged=1;
+			end
+		end
+
+		count=count+1;
+	end
+
+end