Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 2211)
+++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 2212)
@@ -166,4 +166,10 @@
 	parameters->AddObject(param);
 
+	/*stabilize_constraints: */
+	count++;
+	param= new Param(count,"stabilize_constraints",INTEGER);
+	param->SetInteger(iomodel->stabilize_constraints);
+	parameters->AddObject(param);
+
 	/*stokesreconditioning: */
 	count++;
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp	(revision 2211)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp	(revision 2212)
@@ -49,4 +49,5 @@
 	int pengrid_active=0;
 	int pengrid_thermal_steadystate=1;
+	int pengrid_stabilize_constraints=0;
 
 	int numberofpressureloads_stokes;
@@ -143,6 +144,6 @@
 			pengrid_penalty_offset=iomodel->penalty_offset;
 			pengrid_mparid=iomodel->numberofelements+1;//refers to the corresponding parmat property card
-
-			pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_mparid,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate);
+			pengrid_stabilize_constraints=iomodel->stabilize_constraints;
+			pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_mparid,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate,pengrid_stabilize_constraints);
 			
 			loads->AddObject(pengrid);
Index: /issm/trunk/src/c/ModelProcessorx/IoModel.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/IoModel.cpp	(revision 2211)
+++ /issm/trunk/src/c/ModelProcessorx/IoModel.cpp	(revision 2212)
@@ -144,4 +144,5 @@
 	iomodel->thermalconductivity=0;
 	iomodel->min_thermal_constraints=0;
+	iomodel->stabilize_constraints=0;
 	iomodel->mixed_layer_capacity=0;
 	iomodel->thermal_exchange_velocity=0;
@@ -355,7 +356,8 @@
 	IoModelFetchData((void**)&iomodel->thermalconductivity,NULL,NULL,iomodel_handle,"thermalconductivity","Scalar",NULL);
 	IoModelFetchData((void**)&iomodel->min_thermal_constraints,NULL,NULL,iomodel_handle,"min_thermal_constraints","Integer",NULL);
+	IoModelFetchData((void**)&iomodel->stabilize_constraints,NULL,NULL,iomodel_handle,"stabilize_constraints","Integer",NULL);
 	IoModelFetchData((void**)&iomodel->mixed_layer_capacity,NULL,NULL,iomodel_handle,"mixed_layer_capacity","Scalar",NULL);
 	IoModelFetchData((void**)&iomodel->thermal_exchange_velocity,NULL,NULL,iomodel_handle,"thermal_exchange_velocity","Scalar",NULL);
-
+	
 	/*qmu: */
 	if(iomodel->qmu_analysis){
Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 2211)
+++ /issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 2212)
@@ -140,4 +140,5 @@
 	double  heatcapacity,thermalconductivity;
 	int    min_thermal_constraints;
+	int    stabilize_constraints;
 	double mixed_layer_capacity;
 	double thermal_exchange_velocity;
Index: /issm/trunk/src/c/ModelProcessorx/Melting/CreateLoadsMelting.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Melting/CreateLoadsMelting.cpp	(revision 2211)
+++ /issm/trunk/src/c/ModelProcessorx/Melting/CreateLoadsMelting.cpp	(revision 2212)
@@ -38,4 +38,5 @@
 	int pengrid_active=0;
 	int pengrid_thermal_steadystate=1;
+	int pengrid_stabilize_constraints=0;
 
 	int count;
@@ -65,6 +66,7 @@
 			pengrid_active=0;
 			pengrid_mparid=iomodel->numberofelements+1;//refers to the corresponding parmat property card
+			pengrid_stabilize_constraints=iomodel->stabilize_constraints;
 			
-			pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_mparid,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate);
+			pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_mparid,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate,pengrid_stabilize_constraints);
 			
 			loads->AddObject(pengrid);
Index: /issm/trunk/src/c/ModelProcessorx/Thermal/CreateLoadsThermal.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Thermal/CreateLoadsThermal.cpp	(revision 2211)
+++ /issm/trunk/src/c/ModelProcessorx/Thermal/CreateLoadsThermal.cpp	(revision 2212)
@@ -38,4 +38,5 @@
 	int pengrid_active=0;
 	int pengrid_thermal_steadystate=1;
+	int pengrid_stabilize_constraints=0;
 
 	int numberofsegs_diag_stokes;
@@ -66,6 +67,7 @@
 			pengrid_penalty_offset=iomodel->penalty_offset;
 			pengrid_active=0;
+			pengrid_stabilize_constraints=iomodel->stabilize_constraints;
 			
-			pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_mparid,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate);
+			pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_mparid,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate,pengrid_stabilize_constraints);
 			
 			loads->AddObject(pengrid);
Index: /issm/trunk/src/c/objects/Pengrid.cpp
===================================================================
--- /issm/trunk/src/c/objects/Pengrid.cpp	(revision 2211)
+++ /issm/trunk/src/c/objects/Pengrid.cpp	(revision 2212)
@@ -22,5 +22,5 @@
 }
 
-Pengrid::Pengrid(int	pengrid_id, int pengrid_node_id,int pengrid_mparid, int pengrid_dof, int pengrid_active, double pengrid_penalty_offset,int pengrid_thermal_steadystate){
+Pengrid::Pengrid(int	pengrid_id, int pengrid_node_id,int pengrid_mparid, int pengrid_dof, int pengrid_active, double pengrid_penalty_offset,int pengrid_thermal_steadystate,int pengrid_stabilize_constraints){
 	
 	id=pengrid_id;
@@ -31,4 +31,5 @@
 	penalty_offset =pengrid_penalty_offset;
 	thermal_steadystate=pengrid_thermal_steadystate;
+	stabilize_constraints=pengrid_stabilize_constraints;
 
 	node_offset=UNDEF;
@@ -36,4 +37,6 @@
 	matpar=NULL;
 	matpar_offset=UNDEF;
+
+	zigzag_counter=0;
 
 	return;
@@ -600,4 +603,13 @@
 	else{
 		unstable=1;
+		if(stabilize_constraints)zigzag_counter++;
+	}
+
+	/*If penalty keeps zigzagging more than 5 times: */
+	if(stabilize_constraints){
+		if(zigzag_counter>5){
+			unstable=0;
+			active=1;
+		}
 	}
 
Index: /issm/trunk/src/c/objects/Pengrid.h
===================================================================
--- /issm/trunk/src/c/objects/Pengrid.h	(revision 2211)
+++ /issm/trunk/src/c/objects/Pengrid.h	(revision 2212)
@@ -30,8 +30,11 @@
 		int   matpar_offset;
 
+		int stabilize_constraints;
+		int zigzag_counter;
+
 	public:
 
 		Pengrid();
-		Pengrid(int	id, int node_id,int mparid,int dof, int active, double penalty_offset,int thermal_steadystate);
+		Pengrid(int	id, int node_id,int mparid,int dof, int active, double penalty_offset,int thermal_steadystate,int stabilize_constraints);
 		~Pengrid();
 
Index: /issm/trunk/src/m/classes/@model/model.m
===================================================================
--- /issm/trunk/src/m/classes/@model/model.m	(revision 2211)
+++ /issm/trunk/src/m/classes/@model/model.m	(revision 2212)
@@ -110,4 +110,6 @@
 	md.mixed_layer_capacity=0;
 	md.thermal_exchange_velocity=0;
+	md.min_thermal_constraints=0;
+	md.stabilize_constraints=0;
 	
 	%Physical parameters
@@ -139,5 +141,4 @@
 	md.spctemperature=NaN;
 	md.spcthickness=NaN;
-	md.min_thermal_constraints=0;
 
 	%Observations
Index: /issm/trunk/src/m/classes/public/display/displaythermal.m
===================================================================
--- /issm/trunk/src/m/classes/public/display/displaythermal.m	(revision 2211)
+++ /issm/trunk/src/m/classes/public/display/displaythermal.m	(revision 2212)
@@ -11,4 +11,8 @@
 disp(sprintf('   Thermal solution parameters:'));
 
+disp(sprintf('\n      parameters:'));
+fielddisplay(md,'stabilize_constraints','stabilize unstable thermal constraints that keep zigzagging (default is 0)');
+fielddisplay(md,'min_thermal_constraints','threshold to declare convergence of thermal solution (default is 0)');
+
 disp(sprintf('\n      boundary conditions:'));
 fielddisplay(md,'spctemperature','constraints flag list (first column) and values (second column)');
Index: /issm/trunk/src/m/classes/public/marshall.m
===================================================================
--- /issm/trunk/src/m/classes/public/marshall.m	(revision 2211)
+++ /issm/trunk/src/m/classes/public/marshall.m	(revision 2212)
@@ -140,4 +140,5 @@
 WriteData(fid,md.mixed_layer_capacity,'Scalar','mixed_layer_capacity');
 WriteData(fid,md.thermal_exchange_velocity,'Scalar','thermal_exchange_velocity');
+WriteData(fid,md.stabilize_constraints,'Integer','stabilize_constraints');
 
 %elements type
Index: /issm/trunk/src/m/utils/OS/oshostname.m
===================================================================
--- /issm/trunk/src/m/utils/OS/oshostname.m	(revision 2211)
+++ /issm/trunk/src/m/utils/OS/oshostname.m	(revision 2212)
@@ -19,3 +19,4 @@
 	end
 	hostname=hostname(1:end-1);
+	[hostname '|']
 end