Index: /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp	(revision 15078)
+++ /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp	(revision 15079)
@@ -16,5 +16,5 @@
 
 /*Local prototypes*/
-bool steadystateconvergence(FemModel* femmodel);
+bool steadystateconvergence(Vector<IssmDouble>* tg,Vector<IssmDouble>* tg_old,Vector<IssmDouble>* ug,Vector<IssmDouble>* ug_old,IssmDouble reltol);
 
 void steadystate_core(FemModel* femmodel){
@@ -22,10 +22,15 @@
 	/*intermediary: */
 	int step; 
+	Vector<IssmDouble>* ug     = NULL;
+	Vector<IssmDouble>* ug_old = NULL;
+	Vector<IssmDouble>* tg     = NULL;
+	Vector<IssmDouble>* tg_old = NULL;
 
 	/*parameters: */
-	bool save_results,isenthalpy;
-	int  maxiter;
-	int  numoutputs         = 0;
-	int  *requested_outputs = NULL;
+	bool        save_results,isenthalpy;
+	int         maxiter;
+	IssmDouble  reltol;
+	int         numoutputs        = 0;
+	int        *requested_outputs = NULL;
 
 	/* recover parameters:*/
@@ -34,4 +39,5 @@
 	femmodel->parameters->FindParam(&numoutputs,SteadystateNumRequestedOutputsEnum);
 	femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
+	femmodel->parameters->FindParam(&reltol,SteadystateReltolEnum);
 	femmodel->parameters->SetParam(false,SaveResultsEnum);
 	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SteadystateRequestedOutputsEnum);
@@ -46,7 +52,10 @@
 		if(isenthalpy==0){
 			thermal_core(femmodel);
+			femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
+			GetSolutionFromInputsx(&tg,femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
 		}
 		else{
 			enthalpy_core(femmodel);
+			GetSolutionFromInputsx(&tg,femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
 		}
 		#else
@@ -56,8 +65,9 @@
 		if(VerboseSolution()) _pprintLine_("   computing new velocity");
 		diagnostic_core(femmodel);
+		GetSolutionFromInputsx(&ug,femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
 
-		if (step>1){
-			if(VerboseSolution()) _pprintLine_("   checking velocity, temperature and pressure convergence");
-			if(steadystateconvergence(femmodel)) break;
+		if(step>1){
+			if(VerboseSolution()) _pprintLine_("   checking steadystate convergence");
+			if(steadystateconvergence(tg,tg_old,ug,ug_old,reltol)) break;
 		}
 		if(step>maxiter){
@@ -66,12 +76,7 @@
 		}
 
-		if(VerboseSolution()) _pprintLine_("   saving velocity, temperature and pressure to check for convergence at next step");
-		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,VxPicardEnum);
-		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,VyPicardEnum);
-		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,VzPicardEnum);
-		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressurePicardEnum);
-		InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureOldEnum);
-
-		//increase counter
+		/*update results and increase counter*/
+		delete tg_old;tg_old=tg;
+		delete ug_old;ug_old=ug;
 		step++;
 	}
@@ -92,29 +97,46 @@
 
 	/*Free ressources:*/
+	delete tg;
+	delete ug;
 	xDelete<int>(requested_outputs);
 }
-bool steadystateconvergence(FemModel* femmodel){
+bool steadystateconvergence(Vector<IssmDouble>* tg,Vector<IssmDouble>* tg_old,Vector<IssmDouble>* ug,Vector<IssmDouble>* ug_old,IssmDouble reltol){
 
-	/*output: */
-	bool converged=false;
-	bool velocity_converged=false;
-	bool temperature_converged=false;
+	/*Output*/
+	bool converged = true;
 
-	/*intermediary: */
-	int velocityenums[8]={VxEnum,VxPicardEnum,VyEnum,VyPicardEnum,VzEnum,VzPicardEnum,PressureEnum,PressurePicardEnum}; //pairs of enums (new and old) on which to carry out the converence tests
-	int temperatureenums[2]={TemperatureEnum,TemperatureOldEnum};
-	int convergencecriterion[1]={RelativeEnum}; //criterions for convergence, RelativeEnum or AbsoluteEnum
-	IssmDouble convergencecriterionvalue[1]; //value of criterion to be respected
+	/*Intermediary*/
+	Vector<IssmDouble>* dug    = NULL;
+	Vector<IssmDouble>* dtg    = NULL;
+	IssmDouble          ndt,nt;
+	IssmDouble          ndu,nu;
 
-	/*retrieve parameters: */
-	femmodel->parameters->FindParam(&convergencecriterionvalue[0],SteadystateReltolEnum);
+	/*compute norm(du)/norm(u)*/
+	dug=ug_old->Duplicate(); ug_old->Copy(dug); dug->AYPX(ug,-1.0);
+	ndu=dug->Norm(NORM_TWO); nu=ug_old->Norm(NORM_TWO);
+	if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
+	if((ndu/nu)<reltol){
+		if(VerboseConvergence()) _pprintLine_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<ndu/nu*100<<" < "<<reltol*100<<" %");
+	}
+	else{ 
+		if(VerboseConvergence()) _pprintLine_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<ndu/nu*100<<" > "<<reltol*100<<" %");
+		converged=false;
+	}
 
-	/*figure out convergence at the input level, because we don't have the solution vectors!: */
-	velocity_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,&velocityenums[0],8,&convergencecriterion[0],&convergencecriterionvalue[0],1);
-	temperature_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,&temperatureenums[0],2,&convergencecriterion[0],&convergencecriterionvalue[0],1);
+	/*compute norm(dt)/norm(t)*/
+	dtg=tg_old->Duplicate(); tg_old->Copy(dtg); dtg->AYPX(tg,-1.0);
+	ndt=dtg->Norm(NORM_TWO); nt=tg_old->Norm(NORM_TWO);
+	if (xIsNan<IssmDouble>(ndt) || xIsNan<IssmDouble>(nt)) _error_("convergence criterion is NaN!");
+	if((ndt/nt)<reltol){
+		if(VerboseConvergence()) _pprintLine_(setw(50)<<left<<"   Temperature convergence: norm(dt)/norm(t)"<<ndt/nt*100<<" < "<<reltol*100<<" %");
+	}
+	else{ 
+		if(VerboseConvergence()) _pprintLine_(setw(50)<<left<<"   Temperature convergence: norm(dt)/norm(t)"<<ndt/nt*100<<" > "<<reltol*100<<" %");
+		converged=false;
+	}
 
-	if(velocity_converged && temperature_converged) converged=true;
-
-	/*return: */
+	/*clean up and return*/
+	delete dtg;
+	delete dug;
 	return converged;
 }