Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18384)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18385)
@@ -562,6 +562,7 @@
 	StressTensoryyEnum,
 	StressTensoryzEnum,
 	StressTensorzzEnum,
+	StressMaxPrincipalEnum,
 	DeviatoricStressEnum,
 	DeviatoricStressxxEnum,
 	DeviatoricStressxyEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18384)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18385)
@@ -553,6 +553,7 @@
 		case StressTensoryyEnum : return "StressTensoryy";
 		case StressTensoryzEnum : return "StressTensoryz";
 		case StressTensorzzEnum : return "StressTensorzz";
+		case StressMaxPrincipalEnum : return "StressMaxPrincipal";
 		case DeviatoricStressEnum : return "DeviatoricStress";
 		case DeviatoricStressxxEnum : return "DeviatoricStressxx";
 		case DeviatoricStressxyEnum : return "DeviatoricStressxy";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18384)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18385)
@@ -565,6 +565,7 @@
 	      else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
 	      else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
 	      else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
+	      else if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum;
 	      else if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum;
 	      else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum;
 	      else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
@@ -627,11 +628,11 @@
 	      else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
 	      else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
 	      else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
-	      else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"MinVel")==0) return MinVelEnum;
+	      if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
+	      else if (strcmp(name,"MinVel")==0) return MinVelEnum;
 	      else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
 	      else if (strcmp(name,"MinVx")==0) return MinVxEnum;
 	      else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
Index: ../trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 18384)
+++ ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 18385)
@@ -148,6 +148,7 @@
 		void       TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum);
 		void       TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array);
 		void       ViscousHeatingCreateInput(void);
+		void       StressMaxPrincipalCreateInput(void);
 		void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
 		void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
 		void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18384)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18385)
@@ -244,15 +244,15 @@
 /*}}}*/
 void       Penta::ComputeStressTensor(){/*{{{*/
 
-	IssmDouble      xyz_list[NUMVERTICES][3];
-	IssmDouble      pressure,viscosity;
-	IssmDouble      epsilon[6]; /* epsilon=[exx,eyy,exy];*/
-	IssmDouble      sigma_xx[NUMVERTICES];
-	IssmDouble		sigma_yy[NUMVERTICES];
-	IssmDouble		sigma_zz[NUMVERTICES];
-	IssmDouble      sigma_xy[NUMVERTICES];
-	IssmDouble		sigma_xz[NUMVERTICES];
-	IssmDouble		sigma_yz[NUMVERTICES];
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	IssmDouble  pressure,viscosity;
+	IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,exy];*/
+	IssmDouble  sigma_xx[NUMVERTICES];
+	IssmDouble	sigma_yy[NUMVERTICES];
+	IssmDouble	sigma_zz[NUMVERTICES];
+	IssmDouble  sigma_xy[NUMVERTICES];
+	IssmDouble	sigma_xz[NUMVERTICES];
+	IssmDouble	sigma_yz[NUMVERTICES];
 	GaussPenta* gauss=NULL;
 
 	/* Get node coordinates and dof list: */
Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18384)
+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18385)
@@ -1157,6 +1157,10 @@
 				this->ViscousHeatingCreateInput();
 				input=this->inputs->GetInput(output_enum);
 				break;
+			case StressMaxPrincipalEnum:
+				this->StressMaxPrincipalCreateInput();
+				input=this->inputs->GetInput(output_enum);
+				break;
 			case StressTensorxxEnum: 
 			case StressTensorxyEnum: 
 			case StressTensorxzEnum: 
@@ -1337,6 +1341,92 @@
 	_assert_(matpar);
 	return this->matpar->TMeltingPoint(pressure);
 }/*}}}*/
+void       Element::StressMaxPrincipalCreateInput(void){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble *xyz_list = NULL;
+	IssmDouble  sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz;
+	IssmDouble  a,b,c,d,x[3],max;
+	int         dim,numroots;
+
+	/*First: get stress tensor*/
+	this->ComputeStressTensor();
+
+	/*Get domain dimension*/
+	this->FindParam(&dim,DomainDimensionEnum);
+
+	/*Fetch number vertices and allocate memory*/
+	int         numvertices  = this->GetNumberOfVertices();
+	IssmDouble* maxprincipal = xNew<IssmDouble>(numvertices);
+
+	/*Retrieve all inputs and parameters*/
+	this->GetVerticesCoordinatesBase(&xyz_list);
+	Input* sigma_xx_input  = this->GetInput(StressTensorxxEnum); _assert_(sigma_xx_input);
+	Input* sigma_yy_input  = this->GetInput(StressTensoryyEnum); _assert_(sigma_yy_input);
+	Input* sigma_xy_input  = this->GetInput(StressTensorxyEnum); _assert_(sigma_xy_input);
+	Input* sigma_xz_input  = NULL;
+	Input* sigma_yz_input  = NULL;
+	Input* sigma_zz_input  = NULL;
+	if(dim==3){
+		sigma_xz_input  = this->GetInput(StressTensorxzEnum); _assert_(sigma_xz_input);
+		sigma_yz_input  = this->GetInput(StressTensoryzEnum); _assert_(sigma_yz_input);
+		sigma_zz_input  = this->GetInput(StressTensorzzEnum); _assert_(sigma_zz_input);
+	}
+
+	/*loop over vertices: */
+	Gauss* gauss=this->NewGauss();
+	for (int iv=0;iv<numvertices;iv++){
+		gauss->GaussVertex(iv);
+
+		sigma_xx_input->GetInputValue(&sigma_xx,gauss);
+		sigma_yy_input->GetInputValue(&sigma_yy,gauss);
+		sigma_xy_input->GetInputValue(&sigma_xy,gauss);
+		if(dim==3){
+			sigma_xz_input->GetInputValue(&sigma_xz,gauss);
+			sigma_yz_input->GetInputValue(&sigma_yz,gauss);
+			sigma_zz_input->GetInputValue(&sigma_zz,gauss);
+		}
+
+		if(dim==2){
+			a = 0.;
+			b = 1.;
+			c = -sigma_yy -sigma_xx;
+			d = sigma_xx*sigma_yy - sigma_xy*sigma_xy;
+		}
+		else{
+			a = -1.;
+			b = sigma_xx+sigma_yy+sigma_zz;
+			c = -sigma_xx*sigma_yy -sigma_xx*sigma_zz -sigma_yy*sigma_zz + sigma_xy*sigma_xy +sigma_xz*sigma_xz +sigma_yz*sigma_yz;
+			d = sigma_xx*sigma_yy*sigma_zz - sigma_xx*sigma_yz*sigma_yz -sigma_yy*sigma_xz*sigma_xz - sigma_zz*sigma_xy*sigma_xy + 2.*sigma_xy*sigma_xz*sigma_yz;
+		}
+
+		/*Get roots of polynomials*/
+		cubic(a,b,c,d,x,&numroots);
+
+		/*Initialize maximum eigne value*/
+		if(numroots>0){
+			max = x[0];
+		}
+		else{
+			_error_("No eigen value found");
+		}
+
+		/*Get max*/
+		for(int i=1;i<numroots;i++){
+			if(x[i]>max) max = x[i];
+		}
+
+		maxprincipal[iv]=max;
+	}
+
+	/*Create input*/
+	this->AddInput(StressMaxPrincipalEnum,maxprincipal,P1Enum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(maxprincipal);
+	delete gauss;
+}
+/*}}}*/
 void       Element::ViscousHeatingCreateInput(void){/*{{{*/
 
 	/*Intermediaries*/
Index: ../trunk-jpl/src/m/enum/StressMaxPrincipalEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/StressMaxPrincipalEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/StressMaxPrincipalEnum.m	(revision 18385)
@@ -0,0 +1,11 @@
+function macro=StressMaxPrincipalEnum()
+%STRESSMAXPRINCIPALENUM - Enum of StressMaxPrincipal
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=StressMaxPrincipalEnum()
+
+macro=StringToEnum('StressMaxPrincipal');
Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18384)
+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18385)
@@ -545,6 +545,7 @@
 def StressTensoryyEnum(): return StringToEnum("StressTensoryy")[0]
 def StressTensoryzEnum(): return StringToEnum("StressTensoryz")[0]
 def StressTensorzzEnum(): return StringToEnum("StressTensorzz")[0]
+def StressMaxPrincipalEnum(): return StringToEnum("StressMaxPrincipal")[0]
 def DeviatoricStressEnum(): return StringToEnum("DeviatoricStress")[0]
 def DeviatoricStressxxEnum(): return StringToEnum("DeviatoricStressxx")[0]
 def DeviatoricStressxyEnum(): return StringToEnum("DeviatoricStressxy")[0]
