Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16276)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16277)
@@ -4808,14 +4808,14 @@
 
 	/*Intermediaries */
-	int        i,j;
-	int        analysis_type;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0};
-	IssmDouble Jdet2d,dt;
-	IssmDouble rho_ice,heatcapacity,geothermalflux_value;
-	IssmDouble basalfriction,alpha2,vx,vy;
-	IssmDouble basis[NUMVERTICES];
-	IssmDouble scalar;
-	Friction*  friction=NULL;
+	int         i,j;
+	int         analysis_type;
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	IssmDouble  xyz_list_tria[NUMVERTICES2D][3]={0.0};
+	IssmDouble  Jdet2d,dt;
+	IssmDouble  rho_ice,heatcapacity,geothermalflux_value;
+	IssmDouble  basalfriction,alpha2,vx,vy;
+	IssmDouble  basis[NUMVERTICES];
+	IssmDouble  scalar;
+	Friction*   friction=NULL;
 	GaussPenta* gauss=NULL;
 
@@ -4880,5 +4880,5 @@
 	IssmDouble  B_average,s_average;
 	int        *doflist = NULL;
-	IssmDouble  pressure[numdof];
+	bool        hack    = false;
 
 	/*Get dof list: */
@@ -4888,7 +4888,4 @@
 	for(i=0;i<numdof;i++){
 		values[i]=solution[doflist[i]];
-		GetInputListOnVertices(&pressure[0],PressureEnum);
-		if(values[i]>matpar->TMeltingPoint(pressure[i])) values[i]=matpar->TMeltingPoint(pressure[i]);
-		if(values[i]<matpar->TMeltingPoint(pressure[i])-50) values[i]=matpar->TMeltingPoint(pressure[i])-50;
 
 		/*Check solution*/
@@ -4897,4 +4894,13 @@
 		//if(values[i]>275)    _printf_("temperature > 275°K found in solution vector (Paterson's rheology associated is negative)\n");
 	}
+
+	if(hack){
+		/*Force temperature between [Tpmp-50 Tpmp] to disable penalties*/
+		IssmDouble pressure[numdof];
+		GetInputListOnVertices(&pressure[0],PressureEnum);
+		if(values[i]>matpar->TMeltingPoint(pressure[i]))     values[i]=matpar->TMeltingPoint(pressure[i]);
+		if(values[i]<matpar->TMeltingPoint(pressure[i])-50.) values[i]=matpar->TMeltingPoint(pressure[i])-50.;
+	}
+
 
 	/*Get all inputs and parameters*/
