Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11298)
+++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11299)
@@ -2050,4 +2050,6 @@
 				name==InversionVzObsEnum ||
 				name==TemperatureEnum ||
+				name==EnthalpyEnum ||
+				name==EnthalpyPicardEnum ||
 				name==WaterfractionEnum||
 				name==FrictionCoefficientEnum ||
@@ -3241,5 +3243,5 @@
 	this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
 	Input* pressure_input=inputs->GetInput(PressureEnum);      _assert_(pressure_input);
-	Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);      _assert_(enthalpy_input);
+	Input* enthalpy_input=inputs->GetInput(EnthalpyPicardEnum);_assert_(enthalpy_input); //for this iteration of the step
 	Input* vx_input=inputs->GetInput(VxEnum);                  _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);                  _assert_(vy_input);
@@ -4239,5 +4241,5 @@
 	const int    numdof=NDOF1*NUMVERTICES;
 
-	bool   converged;
+	bool   converged=false;
 	int    i,rheology_law;
 	double xyz_list[NUMVERTICES][3];
@@ -4266,7 +4268,6 @@
 	Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
 	
-
-//	this->inputs->GetInputValue(&converged,ConvergedEnum);
-//	if(converged){
+	this->inputs->GetInputValue(&converged,ConvergedEnum);
+	if(converged){
 		/*Convert enthalpy into temperature and water fraction*/
 		for(i=0;i<numdof;i++) matpar->EnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]);
@@ -4300,8 +4301,8 @@
 
 		}
-//	}
-//	else{
-//		this->inputs->AddInput(new PentaP1Input(EnthalpyPicardEnum,values));
-//	}
+	}
+	else{
+		this->inputs->AddInput(new PentaP1Input(EnthalpyPicardEnum,values));
+	}
 
 	/*Free ressources:*/
Index: /issm/trunk-jpl/src/c/objects/Materials/Matpar.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Materials/Matpar.cpp	(revision 11298)
+++ /issm/trunk-jpl/src/c/objects/Materials/Matpar.cpp	(revision 11299)
@@ -420,4 +420,5 @@
 	}
 	else{
+		printf("-------------- file: Matpar.cpp line: %i\n",__LINE__); 
 		return 0.1*thermalconductivity/(rho_ice*heatcapacity);
 	}
Index: /issm/trunk-jpl/src/c/solutions/enthalpy_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/enthalpy_core.cpp	(revision 11298)
+++ /issm/trunk-jpl/src/c/solutions/enthalpy_core.cpp	(revision 11299)
@@ -46,6 +46,5 @@
 
 		/*call enthalpy_core_step: */
-		femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
-		solver_linear(femmodel);
+		enthalpy_core_step(femmodel,i,time);
 
 		if(solution_type==EnthalpySolutionEnum && !control_analysis){
Index: /issm/trunk-jpl/src/c/solutions/enthalpy_core_step.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/enthalpy_core_step.cpp	(revision 11298)
+++ /issm/trunk-jpl/src/c/solutions/enthalpy_core_step.cpp	(revision 11299)
@@ -21,6 +21,6 @@
 	solver_nonlinear(femmodel,modify_loads);
 
-	/*transfer water column thickness to old water column thickness: */
-	//InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum,EnthalpyOldEnum);
+	/*transfer enthalpy to enthalpy picard for the next step: */
+	InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum,EnthalpyPicardEnum);
 
 }
Index: /issm/trunk-jpl/src/m/solutions/enthalpy_core.m
===================================================================
--- /issm/trunk-jpl/src/m/solutions/enthalpy_core.m	(revision 11298)
+++ /issm/trunk-jpl/src/m/solutions/enthalpy_core.m	(revision 11299)
@@ -4,5 +4,4 @@
 %   Usage:
 %      femmodel=enthalpy_core(femmodel)
-
 
 	%recover parameters common to all solutions
@@ -25,7 +24,5 @@
 		time=(i+1)*dt;
 
-		issmprintf(VerboseSolution,'\n%s',['   computing enthalpy']);
-		femmodel=SetCurrentConfiguration(femmodel,EnthalpyAnalysisEnum);
-		femmodel=solver_linear(femmodel);
+		femmodel=enthalpy_core_step(femmodel,i,time); 
 
 		if (solution_type==EnthalpySolutionEnum & ~control_analysis),
Index: /issm/trunk-jpl/src/m/solutions/enthalpy_core_step.m
===================================================================
--- /issm/trunk-jpl/src/m/solutions/enthalpy_core_step.m	(revision 11299)
+++ /issm/trunk-jpl/src/m/solutions/enthalpy_core_step.m	(revision 11299)
@@ -0,0 +1,14 @@
+function femmodel=enthalpy_core_step(femmodel,step,time)
+%ENTHALPY_CORE_STEP - core of the enthalpy solution for one step 
+%
+%   Usage:
+%      femmodel=enthalpy_core_step(femmodel,step,time)
+
+	%some parameters
+	modify_loads=true;
+
+	issmprintf(VerboseSolution,'\n%s',['   computing enthalpy']);
+	femmodel=SetCurrentConfiguration(femmodel,EnthalpyAnalysisEnum);
+	femmodel=solver_nonlinear(femmodel,modify_loads);
+
+end %end function
Index: /issm/trunk-jpl/src/m/solvers/solver_nonlinear.m
===================================================================
--- /issm/trunk-jpl/src/m/solvers/solver_nonlinear.m	(revision 11298)
+++ /issm/trunk-jpl/src/m/solvers/solver_nonlinear.m	(revision 11299)
@@ -22,4 +22,5 @@
 
 	%Update the solution to make sure that vx and vxold are similar
+	[femmodel.elements loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,double(converged),ConvergedEnum);
 	[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug);
 
@@ -38,9 +39,4 @@
 		ug= Mergesolutionfromftog( uf, ys, femmodel.nodes,femmodel.parameters); 
 
-		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug);
-		[loads,constraints_converged,num_unstable_constraints] =ConstraintsState( femmodel.elements,femmodel.nodes,femmodel.vertices,loads, femmodel.materials,femmodel.parameters);
-
-		issmprintf(VerboseConvergence(),'%s%i','      number of unstable constraints: ',num_unstable_constraints);
-
 		%Figure out if convergence have been reached
 		converged=convergence(K_ff,p_f,uf,old_uf,femmodel.parameters);
@@ -48,4 +44,9 @@
 		%add convergence status into  status
 		[femmodel.elements loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,double(converged),ConvergedEnum);
+
+		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug);
+		[loads,constraints_converged,num_unstable_constraints] =ConstraintsState( femmodel.elements,femmodel.nodes,femmodel.vertices,loads, femmodel.materials,femmodel.parameters);
+
+		issmprintf(VerboseConvergence(),'%s%i','      number of unstable constraints: ',num_unstable_constraints);
 
 		%rift convergence