Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 16698)
+++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 16699)
@@ -116,4 +116,81 @@
 }/*}}}*/
 void ThermalAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-	_error_("not implemented yet");
+
+	bool        converged;
+	int         i,rheology_law;
+	IssmDouble  B_average,s_average,T_average=0.;
+	int        *doflist   = NULL;
+	IssmDouble *xyz_list  = NULL;
+	bool        hack      = false;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Fetch dof list and allocate solution vector*/
+	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	IssmDouble* values    = xNew<IssmDouble>(numnodes);
+
+	/*Use the dof list to index into the solution vector: */
+	for(i=0;i<numnodes;i++){
+		values[i]=solution[doflist[i]];
+
+		/*Check solution*/
+		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
+		//if(values[i]<0)      _printf_("temperature < 0°K found in solution vector\n");
+		//if(values[i]>275)    _printf_("temperature > 275°K found in solution vector (Paterson's rheology associated is negative)\n");
+	}
+
+	/*Force temperature between [Tpmp-50 Tpmp] to disable penalties*/
+	if(hack){
+		IssmDouble* pressure = xNew<IssmDouble>(numnodes);
+		element->GetInputListOnVertices(pressure,PressureEnum);
+		for(i=0;i<numnodes;i++){
+			if(values[i]>element->TMeltingPoint(pressure[i]))     values[i]=element->TMeltingPoint(pressure[i]);
+			if(values[i]<element->TMeltingPoint(pressure[i])-50.) values[i]=element->TMeltingPoint(pressure[i])-50.;
+		}
+		xDelete<IssmDouble>(pressure);
+	}
+
+	/*Get all inputs and parameters*/
+
+	element->GetInputValue(&converged,ConvergedEnum);
+	if(converged){
+		element->AddInput(TemperatureEnum,values,P1Enum);
+
+		/*Update Rheology only if converged (we must make sure that the temperature is below melting point
+		 * otherwise the rheology could be negative*/
+		element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
+		switch(rheology_law){
+			case NoneEnum:
+				/*Do nothing: B is not temperature dependent*/
+				break;
+			case PatersonEnum:
+				for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
+				B_average=Paterson(T_average);
+				element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
+				break;
+			case ArrheniusEnum:{
+				Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input);
+				surface_input->GetInputAverage(&s_average);
+				element->GetVerticesCoordinates(&xyz_list);
+				for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
+				//B_average=Arrhenius(T_average,
+							//s_average-((xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2]+xyz_list[3][2]+xyz_list[4][2]+xyz_list[5][2])/6.0),
+							//element->GetMaticeParameter(MaterialsRheologyNEnum));
+				element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
+				break;
+				}
+			default:
+				_error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
+
+		}
+	}
+	else{
+		element->AddInput(TemperaturePicardEnum,values,P1Enum);
+	}
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(values);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<int>(doflist);
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16699)
@@ -40,4 +40,5 @@
 		virtual void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
 		virtual void   AddInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
+		virtual void   AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
 		virtual void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>*  Kfs)=0;
 		virtual void   CreateDVector(Vector<IssmDouble>* df)=0;
@@ -95,4 +96,5 @@
 		virtual void   Delta18oParameterization(void)=0;
 		virtual void   SmbGradients()=0;
+		virtual IssmDouble TMeltingPoint(IssmDouble pressure)=0;
 		virtual void   ResetCoordinateSystem()=0;
 		virtual int    VelocityInterpolation()=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16699)
@@ -124,8 +124,16 @@
 void  Penta::AddInput(int input_enum,IssmDouble* values, int interpolation_enum){
 
-	/*Call inputs method*/
+	_assert_(this->inputs);
 	this->inputs->AddInput(new PentaInput(input_enum,values,interpolation_enum));
 }
 /*}}}*/
+/*FUNCTION Penta::AddMaterialInput{{{*/
+void  Penta::AddMaterialInput(int input_enum,IssmDouble* values, int interpolation_enum){
+
+	_assert_(this->material);
+	this->material->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
+}
+/*}}}*/
+
 /*FUNCTION Penta::BedNormal {{{*/
 void Penta::BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]){
@@ -3151,4 +3159,12 @@
 	  /*Update inputs*/
 	  this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&smb[0],P1Enum));
+}
+/*}}}*/
+/*FUNCTION Penta::TMeltingPoint{{{*/
+IssmDouble Penta::TMeltingPoint(IssmDouble pressure){
+
+	_assert_(matpar);
+	return this->matpar->TMeltingPoint(pressure);
+
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16699)
@@ -112,4 +112,5 @@
 		void   ResetCoordinateSystem(void);
 		void   SmbGradients();
+		IssmDouble  TMeltingPoint(IssmDouble pressure);
 		IssmDouble SurfaceArea(void);
 		void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
@@ -182,4 +183,5 @@
 		/*Penta specific routines:{{{*/
 		void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
+		void           AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum);
 		void	         BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]);
 		ElementMatrix* CreateBasalMassMatrix(void);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16699)
@@ -67,4 +67,5 @@
 		/*Element virtual functions definitions: {{{*/
 		void        AddInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
+		void        AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
 		void        ComputeBasalStress(Vector<IssmDouble>* sigma_b){_error_("not implemented yet");};
 		void        ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};
@@ -137,4 +138,5 @@
 		void        ResetCoordinateSystem(void){_error_("not implemented yet");};
 		void	      SmbGradients(){_error_("not implemented yet");};
+		IssmDouble  TMeltingPoint(IssmDouble pressure){_error_("not implemented yet");};
 		IssmDouble  SurfaceArea(void){_error_("not implemented yet");};
 		void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16699)
@@ -170,5 +170,13 @@
 
 	/*Call inputs method*/
+	_assert_(this->inputs);
 	this->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
+}
+/*}}}*/
+/*FUNCTION Tria::AddMaterialInput{{{*/
+void  Tria::AddMaterialInput(int input_enum,IssmDouble* values, int interpolation_enum){
+
+	_assert_(this->material);
+	this->material->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
 }
 /*}}}*/
@@ -2641,4 +2649,12 @@
 	  /*Update inputs*/
 	  this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&smb[0],P1Enum));
+}
+/*}}}*/
+/*FUNCTION Tria::TMeltingPoint{{{*/
+IssmDouble Tria::TMeltingPoint(IssmDouble pressure){
+
+	_assert_(matpar);
+	return this->matpar->TMeltingPoint(pressure);
+
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16699)
@@ -118,4 +118,5 @@
 		void        ResetCoordinateSystem(void);
 		void	      SmbGradients();
+		IssmDouble  TMeltingPoint(IssmDouble pressure);
 		int         VelocityInterpolation();
 		int         PressureInterpolation();
@@ -193,4 +194,5 @@
 		/*Tria specific routines:{{{*/
 		void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
+		void           AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum);
 		ElementMatrix* CreateKMatrix(void);
 		ElementMatrix* CreateKMatrixBalancethickness(void);