Index: /issm/trunk/src/m/solutions/steadystate_core.m
===================================================================
--- /issm/trunk/src/m/solutions/steadystate_core.m	(revision 4515)
+++ /issm/trunk/src/m/solutions/steadystate_core.m	(revision 4516)
@@ -25,7 +25,8 @@
 		femmodel=diagnostic_core(femmodel); 
 
-		displaystring(verbose,'\n%s',['checking temperature, velocity and pressure convergence...']);
 		if step>1,
-			if steadystateconvergence(femmodel),
+			displaystring(verbose,'\n%s',['checking temperature, velocity and pressure convergence...']);
+			converged=steadystateconvergence(femmodel);
+			if converged,
 				break;
 			end
@@ -36,10 +37,9 @@
 		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyOldEnum);
 		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum,VzOldEnum);
-		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,PressureEnum,PressureOlddnum);
+		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,PressureEnum,PressureOldEnum);
 		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,TemperatureEnum,TemperatureOldEnum);
 
 		%Increase counter
 		step=step+1;
-
 	end
 
@@ -49,6 +49,8 @@
 		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum);
 		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum);
+		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VelEnum);
 		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,PressureEnum);
 		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,TemperatureEnum);
+		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,MeltingRateEnum);
 	end
 
Index: /issm/trunk/src/m/solutions/steadystateconvergence.m
===================================================================
--- /issm/trunk/src/m/solutions/steadystateconvergence.m	(revision 4516)
+++ /issm/trunk/src/m/solutions/steadystateconvergence.m	(revision 4516)
@@ -0,0 +1,23 @@
+function converged=steadystateconvergence(femmodel),
+%STEADYSTATECONVERGENCE - determine convergence of steady state solution
+% 
+%   Usage:
+%      steadystateconvergence(femmodel)
+	
+	converged=0;
+	velocity_convergence=0;
+	temperature_convergence=0;
+
+	%Retrieve parameters
+	convergencecriterionvalue=femmodel.parameters.EpsRel;
+	velocityenums=[VxEnum;VxOldEnum;VyEnum;VyOldEnum;VzEnum;VzOldEnum;PressureEnum;PressureOldEnum]; %pairs of enums (new and old) on which to carry out the converence tests
+	temperatureenums=[TemperatureEnum;TemperatureOldEnum];
+	convergencecriterion=[RelativeEnum]; %criterions for convergence, RelativeEnum or AbsoluteEnum 
+
+	%Figure out convergence at the input level, because we don't have the solution vectors
+	temperature_convergence=InputConvergence(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,temperatureenums,convergencecriterion,convergencecriterionvalue);
+	velocity_convergence=InputConvergence(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,velocityenums,convergencecriterion,convergencecriterionvalue);
+
+	if (velocity_convergence & temperature_convergence),
+		converged=1;
+	end