Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 18588)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 18589)
@@ -15,4 +15,6 @@
 
 	parameters->AddObject(iomodel->CopyConstantObject(ThermalStabilizationEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(ThermalMaxiterEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(ThermalReltolEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(ThermalIsdynamicbasalspcEnum));
Index: /issm/trunk-jpl/src/c/cores/thermal_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/thermal_core.cpp	(revision 18588)
+++ /issm/trunk-jpl/src/c/cores/thermal_core.cpp	(revision 18589)
@@ -31,5 +31,5 @@
 		if(VerboseSolution()) _printf0_("   computing enthalpy\n");
 		femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
-		solutionsequence_nonlinear(femmodel,true);
+		solutionsequence_thermal_nonlinear(femmodel);
 
 		/*transfer enthalpy to enthalpy picard for the next step: */
Index: /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp	(revision 18588)
+++ /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp	(revision 18589)
@@ -23,4 +23,5 @@
 	IssmDouble nF;
 	IssmDouble solver_residue,res;
+	int analysis_type;
 
 	/*convergence options*/
@@ -41,5 +42,10 @@
 	/*get convergence options*/
 	parameters->FindParam(&eps_res,StressbalanceRestolEnum);
-	parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
+	// get analysis type
+	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+	if(analysis_type==StressbalanceAnalysisEnum)
+		parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
+	else if(analysis_type==EnthalpyAnalysisEnum)
+		parameters->FindParam(&eps_rel,ThermalReltolEnum);
 	parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
 	parameters->FindParam(&yts,ConstantsYtsEnum);
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp	(revision 18588)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp	(revision 18589)
@@ -3,4 +3,5 @@
  */ 
 
+#include "./solutionsequences.h"
 #include "../toolkits/toolkits.h"
 #include "../classes/classes.h"
@@ -24,4 +25,5 @@
 
 	bool converged;
+	bool isenthalpy, isdynamicbasalspc;
 	int constraints_converged;
 	int num_unstable_constraints;
@@ -29,4 +31,5 @@
 	int thermal_penalty_threshold;
 	int thermal_maxiter;
+	IssmDouble thermal_reltol;
 
 	/*parameters:*/
@@ -34,45 +37,92 @@
 
 	/*Recover parameters: */
-	femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum);
+	femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
 	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
 	femmodel->parameters->FindParam(&thermal_maxiter,ThermalMaxiterEnum);
 
+	converged=false;
+	InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
+
+	if(isenthalpy){
+		femmodel->parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
+		femmodel->parameters->FindParam(&thermal_reltol,ThermalReltolEnum);
+		femmodel->UpdateConstraintsx();
+
+		//Update the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
+		GetSolutionFromInputsx(&tg,femmodel);
+		Reducevectorgtofx(&tf, tg, femmodel->nodes,femmodel->parameters);
+		InputUpdateFromSolutionx(femmodel,tg);
+	}
+	else{
+		femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum);
+		InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum);
+		femmodel->UpdateConstraintsx();
+	}
+
 	count=1;
-	converged=false;
+	
+	for(;;){
+		delete tf_old;tf_old=tf;
+		delete tg;
 
-	InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum);
-	InputUpdateFromConstantx(femmodel,false,ConvergedEnum);
-	femmodel->UpdateConstraintsx();
-
-	for(;;){
-
-		delete tf_old; tf_old=tf;
-		SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel);
+		if(isenthalpy) SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
+		else SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf, Kfs, ys); delete Kfs;
-		Solverx(&tf, Kff, pf,tf_old, df, femmodel->parameters);
-		delete Kff;delete pf;delete tg; delete df;
+		Solverx(&tf, Kff, pf, tf_old, df, femmodel->parameters);
 		Mergesolutionfromftogx(&tg, tf,ys,femmodel->nodes,femmodel->parameters); delete ys;
+		if(isenthalpy){ 
+			convergence(&converged,Kff,pf,tf,tf_old,femmodel->parameters); delete Kff; delete pf; delete df;
+			InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
+		}
 		InputUpdateFromSolutionx(femmodel,tg);
+		ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
+		if(VerboseConvergence()) _printf0_("   number of unstable constraints: " << num_unstable_constraints << "\n");
 
-		ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
-
-		if (!converged){
-			if(VerboseConvergence()) _printf0_("   #unstable constraints = " << num_unstable_constraints << "\n");
-			if (num_unstable_constraints <= thermal_penalty_threshold)converged=true;
-			if (count>=thermal_maxiter){
+		if(isenthalpy){
+			/*Increase count: */
+			count++;
+			if(converged==true){
+				bool max_iteration_state=false;
+				int step; IssmDouble time;
+				femmodel->parameters->FindParam(&time,TimeEnum);
+				femmodel->parameters->FindParam(&step,StepEnum);
+				femmodel->results->AddObject(new GenericExternalResult<bool>(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, step, time));
+				break;
+			}
+			if(count>=thermal_maxiter){
+				_printf0_("   maximum number of nonlinear iterations (" << thermal_maxiter << ") exceeded\n"); 
 				converged=true;
-				_printf0_("   maximum number of iterations (" << thermal_maxiter << ") exceeded\n"); 
+				InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
+				InputUpdateFromSolutionx(femmodel,tg);		
+				bool max_iteration_state=true;
+				int step; IssmDouble time;
+				femmodel->parameters->FindParam(&time,TimeEnum);
+				femmodel->parameters->FindParam(&step,StepEnum);
+				femmodel->results->AddObject(new GenericExternalResult<bool>(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, step, time));
+				break;
 			}
 		}
-		count++;
-
-		InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
-
-		if(converged)break;
+		else{
+			if(!converged){
+				if(num_unstable_constraints<=thermal_penalty_threshold) converged=true;
+				if(count>=thermal_maxiter){
+					converged=true;
+					_printf0_("   maximum number of iterations (" << thermal_maxiter << ") exceeded\n"); 
+				}
+			}
+			count++;
+			InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
+			if(converged)break;
+		}
 	}
 
-	InputUpdateFromSolutionx(femmodel,tg);
-	femmodel->parameters->SetParam(melting_offset,MeltingOffsetEnum);
+	if(isenthalpy){
+		if(VerboseConvergence()) _printf0_("\n   total number of iterations: " << count-1 << "\n");
+	}
+	else{
+		InputUpdateFromSolutionx(femmodel,tg);
+		femmodel->parameters->SetParam(melting_offset,MeltingOffsetEnum);
+	}
 
 	/*Free ressources: */