Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 18952)
+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 18953)
@@ -63,6 +63,10 @@
 				iomodel->FetchDataToInput(elements,CalvinglevermannCoeffEnum);
 				iomodel->FetchDataToInput(elements,CalvinglevermannMeltingrateEnum);
 				break;
+			case CalvingPiEnum:
+				iomodel->FetchDataToInput(elements,CalvingpiCoeffEnum);
+				iomodel->FetchDataToInput(elements,CalvingpiMeltingrateEnum);
+				break;
 			default:
 				_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
 		}
@@ -207,6 +211,22 @@
 				}
 				meltingrate_input = basalelement->GetInput(CalvinglevermannMeltingrateEnum);     _assert_(meltingrate_input);
 				break;
+			case CalvingPiEnum:
+				if(domaintype==Domain2DhorizontalEnum){
+					calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+					calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
+				}
+				else{
+					if(dim==1){
+						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+					}
+					if(dim==2){
+						calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
+						calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
+					}
+				}
+				meltingrate_input = basalelement->GetInput(CalvingpiMeltingrateEnum);     _assert_(meltingrate_input);
+				break;
 			default:
 				_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
 		}
@@ -261,6 +281,7 @@
 					 }
 
 					break;
+
 				case CalvingLevermannEnum:
 					calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity 
 					if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
@@ -273,6 +294,20 @@
 					for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
 					
 					break;
+
+				case CalvingPiEnum:
+					calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity
+					if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
+					meltingrate_input->GetInputValue(&meltingrate,gauss);
+
+					norm_calving=0.;
+					for(i=0;i<dim;i++) norm_calving+=pow(c[i],2);
+					norm_calving=sqrt(norm_calving)+1.e-14;
+
+					for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
+
+					break;
+
 				default:
 					_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
 			}
@@ -348,7 +383,6 @@
 	xDelete<IssmDouble>(v);
 	xDelete<IssmDouble>(w);
 	xDelete<IssmDouble>(c);
-	xDelete<IssmDouble>(m);
 	xDelete<IssmDouble>(dlsf);
 	delete gauss;
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18952)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18953)
@@ -214,10 +214,13 @@
 	CalvingCalvingrateEnum,
 	CalvingMeltingrateEnum,
 	CalvingLevermannEnum,
+	CalvingPiEnum,
 	DefaultCalvingEnum,
 	CalvingRequestedOutputsEnum,
 	CalvinglevermannCoeffEnum,
 	CalvinglevermannMeltingrateEnum,
+	CalvingpiCoeffEnum,
+	CalvingpiMeltingrateEnum,
 	CalvingratexEnum,
 	CalvingrateyEnum,
 	CalvingratexAverageEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18952)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18953)
@@ -222,10 +222,13 @@
 		case CalvingCalvingrateEnum : return "CalvingCalvingrate";
 		case CalvingMeltingrateEnum : return "CalvingMeltingrate";
 		case CalvingLevermannEnum : return "CalvingLevermann";
+		case CalvingPiEnum : return "CalvingPi";
 		case DefaultCalvingEnum : return "DefaultCalving";
 		case CalvingRequestedOutputsEnum : return "CalvingRequestedOutputs";
 		case CalvinglevermannCoeffEnum : return "CalvinglevermannCoeff";
 		case CalvinglevermannMeltingrateEnum : return "CalvinglevermannMeltingrate";
+		case CalvingpiCoeffEnum : return "CalvingpiCoeff";
+		case CalvingpiMeltingrateEnum : return "CalvingpiMeltingrate";
 		case CalvingratexEnum : return "Calvingratex";
 		case CalvingrateyEnum : return "Calvingratey";
 		case CalvingratexAverageEnum : return "CalvingratexAverage";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18952)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18953)
@@ -225,10 +225,13 @@
 	      else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
 	      else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
 	      else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
+	      else if (strcmp(name,"CalvingPi")==0) return CalvingPiEnum;
 	      else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
 	      else if (strcmp(name,"CalvingRequestedOutputs")==0) return CalvingRequestedOutputsEnum;
 	      else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
 	      else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum;
+	      else if (strcmp(name,"CalvingpiCoeff")==0) return CalvingpiCoeffEnum;
+	      else if (strcmp(name,"CalvingpiMeltingrate")==0) return CalvingpiMeltingrateEnum;
 	      else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
 	      else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
 	      else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
@@ -256,13 +259,13 @@
 	      else if (strcmp(name,"MeshNumberofvertices2d")==0) return MeshNumberofvertices2dEnum;
 	      else if (strcmp(name,"MeshNumberofvertices")==0) return MeshNumberofverticesEnum;
 	      else if (strcmp(name,"MeshUpperelements")==0) return MeshUpperelementsEnum;
-	      else if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum;
-	      else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
-	      else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum;
          else stage=3;
    }
    if(stage==3){
-	      if (strcmp(name,"MeshX")==0) return MeshXEnum;
+	      if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum;
+	      else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
+	      else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum;
+	      else if (strcmp(name,"MeshX")==0) return MeshXEnum;
 	      else if (strcmp(name,"MeshY")==0) return MeshYEnum;
 	      else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
 	      else if (strcmp(name,"MeshElementtype")==0) return MeshElementtypeEnum;
@@ -379,13 +382,13 @@
 	      else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;
 	      else if (strcmp(name,"SurfaceforcingsMelt")==0) return SurfaceforcingsMeltEnum;
 	      else if (strcmp(name,"SurfaceforcingsRefreeze")==0) return SurfaceforcingsRefreezeEnum;
-	      else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
-	      else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum;
-	      else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
+	      if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
+	      else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum;
+	      else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
+	      else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
 	      else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
 	      else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
 	      else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
@@ -502,13 +505,13 @@
 	      else if (strcmp(name,"Masscon")==0) return MassconEnum;
 	      else if (strcmp(name,"MassconName")==0) return MassconNameEnum;
 	      else if (strcmp(name,"MassconDefinitionenum")==0) return MassconDefinitionenumEnum;
-	      else if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum;
-	      else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
-	      else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"MassconaxpbyDefinitionenum")==0) return MassconaxpbyDefinitionenumEnum;
+	      if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum;
+	      else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
+	      else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum;
+	      else if (strcmp(name,"MassconaxpbyDefinitionenum")==0) return MassconaxpbyDefinitionenumEnum;
 	      else if (strcmp(name,"MassconaxpbyNamex")==0) return MassconaxpbyNamexEnum;
 	      else if (strcmp(name,"MassconaxpbyNamey")==0) return MassconaxpbyNameyEnum;
 	      else if (strcmp(name,"MassconaxpbyAlpha")==0) return MassconaxpbyAlphaEnum;
@@ -625,13 +628,13 @@
 	      else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
 	      else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
 	      else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
-	      else if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
-	      else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum;
-	      else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum;
+	      if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
+	      else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum;
+	      else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum;
+	      else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum;
 	      else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum;
 	      else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum;
 	      else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum;
@@ -748,13 +751,13 @@
 	      else if (strcmp(name,"Gset")==0) return GsetEnum;
 	      else if (strcmp(name,"Index")==0) return IndexEnum;
 	      else if (strcmp(name,"Indexed")==0) return IndexedEnum;
-	      else if (strcmp(name,"Intersect")==0) return IntersectEnum;
-	      else if (strcmp(name,"Nodal")==0) return NodalEnum;
-	      else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
+	      if (strcmp(name,"Intersect")==0) return IntersectEnum;
+	      else if (strcmp(name,"Nodal")==0) return NodalEnum;
+	      else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
+	      else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
 	      else if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum;
 	      else if (strcmp(name,"RootPath")==0) return RootPathEnum;
 	      else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum;
Index: ../trunk-jpl/src/c/cores/transient_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/transient_core.cpp	(revision 18952)
+++ ../trunk-jpl/src/c/cores/transient_core.cpp	(revision 18953)
@@ -108,6 +108,11 @@
 				femmodel->StrainRateperpendicularx();
 				femmodel->CalvingRateLevermannx();
 			}
+			if(iscalving && calvinglaw==CalvingPiEnum){
+				if(VerboseSolution()) _printf0_("   computing calving rate\n");
+				femmodel->StrainRateparallelx();
+				femmodel->CalvingRatePix();
+			}
 			if(VerboseSolution()) _printf0_("   computing movement of ice boundaries\n");
 			/* smoothen slope of lsf for computation of normal on ice domain*/
 			levelsetfunctionslope_core(femmodel);
Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18952)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18953)
@@ -1734,6 +1734,14 @@
 							   }
 }
 /*}}}*/
+void FemModel::CalvingRatePix(){/*{{{*/
+
+	   for(int i=0;i<elements->Size();i++){
+			      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+					      element->CalvingRatePi();
+							   }
+}
+/*}}}*/
 void FemModel::StrainRateparallelx(){/*{{{*/
 
 	   for(int i=0;i<elements->Size();i++){
Index: ../trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 18952)
+++ ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 18953)
@@ -169,6 +169,7 @@
 		virtual void       AddInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
 		virtual void       AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0;
 		virtual void	    CalvingRateLevermann(void)=0;
+		virtual void       CalvingRatePi(void)=0;
 		virtual IssmDouble CharacteristicLength(void)=0;
 		virtual void       ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;
 		virtual void       ComputeDeviatoricStressTensor(void)=0;
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18952)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18953)
@@ -233,8 +233,8 @@
 		if(calvingrate[iv]<0){
 			calvingrate[iv]=0;
 		}
-		calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6);
-		calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6);
+		calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14);
+		calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14);
 	}
 
 	/*Add input*/
@@ -247,6 +247,84 @@
 
 }
 /*}}}*/
+void       Tria::CalvingRatePi(){/*{{{*/
+
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	GaussTria* gauss=NULL;
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainparallel;
+	IssmDouble  sxx;
+	IssmDouble  sxy;
+	IssmDouble  syy;
+	IssmDouble  sigVM;
+	IssmDouble  thickness;
+	IssmDouble  base;
+	IssmDouble  hAB;
+	IssmDouble  propcoeff;
+	IssmDouble  calvingratex[NUMVERTICES];
+	IssmDouble  calvingratey[NUMVERTICES];
+	IssmDouble  calvingrate[NUMVERTICES];
+
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*Retrieve all inputs and parameters we will need*/
+	Input* vx_input=inputs->GetInput(VxEnum);                                                    _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);                                                    _assert_(vy_input);
+	Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);                        _assert_(strainparallel_input);
+	Input* sxx_input=inputs->GetInput(DeviatoricStressxxEnum);                                   _assert_(sxx_input);
+	Input* sxy_input=inputs->GetInput(DeviatoricStressxyEnum);                                   _assert_(sxy_input);
+	Input* syy_input=inputs->GetInput(DeviatoricStressyyEnum);                                   _assert_(syy_input);
+	Input* thickness_input=inputs->GetInput(ThicknessEnum);                                      _assert_(thickness_input);
+	Input* base_input=inputs->GetInput(BaseEnum);                                                _assert_(base_input);
+	Input* picoeff_input=inputs->GetInput(CalvingpiCoeffEnum);                                   _assert_(picoeff_input);
+
+
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussTria();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/* Get the value we need*/
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vel=vx*vx+vy*vy;
+		strainparallel_input->GetInputValue(&strainparallel,gauss);
+		sxx_input->GetInputValue(&sxx,gauss);
+		sxy_input->GetInputValue(&sxy,gauss);
+		syy_input->GetInputValue(&syy,gauss);
+		thickness_input->GetInputValue(&thickness,gauss);
+		base_input->GetInputValue(&base,gauss);
+		picoeff_input->GetInputValue(&propcoeff,gauss);
+
+		/* Computing sigma Von Mises*/
+		sigVM=sqrt(sxx*sxx+syy*syy+3*sxy*sxy-sxx*syy);
+
+		/* Computing heigth above buoyancy*/
+		hAB=thickness+1028/920*base;
+
+		/*Calving rate for Pi criterion proportionnal to the product of the strain rate along the ice flow direction and the Von Mises stress and the square of the glacier width (hardcoded) divided by the height above buoyancy and the max of the ice velocity power 3 and the ice density (ignored here)*/
+
+		calvingrate[iv]=propcoeff*strainparallel*sigVM*25.e6/hAB/1e9;
+		if(calvingrate[iv]<0){
+			calvingrate[iv]=0;
+		}
+		calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14);
+		calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14);
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+
+}
+/*}}}*/
 IssmDouble Tria::CharacteristicLength(void){/*{{{*/
 
 	return sqrt(2*this->GetArea());
@@ -2583,7 +2661,7 @@
 		strainxy=epsilon[2];
 
 		/*strainparallel= Strain rate along the ice flow direction */
-		strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
+		strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-14);
 	}
 
 	/*Add input*/
@@ -2629,7 +2707,7 @@
 		strainxy=epsilon[2];
 
 		/*strainperpendicular= Strain rate perpendicular to the ice flow direction */
-		strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
+		strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-14);
 	}
 
 	/*Add input*/
Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18952)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18953)
@@ -57,6 +57,7 @@
 		void        ComputeSigmaNN();
 		void        ComputeStressTensor();
 		void        ComputeSurfaceNormalVelocity();
+		void			CalvingRatePi();
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
 		void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
 		void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18952)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18953)
@@ -200,8 +200,8 @@
 		if(calvingrate[iv]<0){
 			calvingrate[iv]=0;
 		}
-		calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6);
-		calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6);
+		calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14);
+		calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14);
 	}
 
 	/*Add input*/
@@ -214,6 +214,81 @@
 
 }
 /*}}}*/
+void       Penta::CalvingRatePi(){/*{{{*/
+
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	GaussPenta* gauss=NULL;
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainparallel;
+	IssmDouble  sxx;
+	IssmDouble  sxy;
+	IssmDouble  syy;
+	IssmDouble  sigVM;
+	IssmDouble  thickness;
+	IssmDouble  base;
+	IssmDouble  hAB;
+	IssmDouble  propcoeff;
+	IssmDouble  calvingratex[NUMVERTICES];
+	IssmDouble  calvingratey[NUMVERTICES];
+	IssmDouble  calvingrate[NUMVERTICES];
+
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*Retrieve all inputs and parameters we will need*/
+	Input* vx_input=inputs->GetInput(VxEnum);                                                    _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);                                                    _assert_(vy_input);
+	Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);                        _assert_(strainparallel_input);
+	Input* sxx_input=inputs->GetInput(DeviatoricStressxxEnum);                                   _assert_(sxx_input);
+	Input* sxy_input=inputs->GetInput(DeviatoricStressxyEnum);                                   _assert_(sxy_input);
+	Input* syy_input=inputs->GetInput(DeviatoricStressyyEnum);                                   _assert_(syy_input);
+	Input* thickness_input=inputs->GetInput(ThicknessEnum);                                      _assert_(thickness_input);
+	Input* base_input=inputs->GetInput(BaseEnum);                                                _assert_(base_input);
+	Input* picoeff_input=inputs->GetInput(CalvingpiCoeffEnum);                                   _assert_(picoeff_input);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussPenta();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/* Get the value we need*/
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vel=vx*vx+vy*vy;
+		strainparallel_input->GetInputValue(&strainparallel,gauss);
+		sxx_input->GetInputValue(&sxx,gauss);
+		sxy_input->GetInputValue(&sxy,gauss);
+		syy_input->GetInputValue(&syy,gauss);
+		thickness_input->GetInputValue(&thickness,gauss);
+		base_input->GetInputValue(&base,gauss);
+		picoeff_input->GetInputValue(&propcoeff,gauss);
+
+		/* Computing sigma Von Mises*/
+		sigVM=sqrt(sxx*sxx+syy*syy+3*sxy*sxy-sxx*syy);
+
+		/* Computing heigth above buoyancy*/
+		hAB=thickness+1028/920*base;
+
+		/*Calving rate for Pi criterion proportionnal to the product of the strain rate along the ice flow direction and the Von Mises stress and the square of the glacier width (hardcoded) divided by the height above buoyancy and the max of the ice velocity power 3 and the ice density (ignored here)*/
+		calvingrate[iv]=propcoeff*strainparallel*sigVM*25.e6/hAB/1e9;
+		if(calvingrate[iv]<0){
+			calvingrate[iv]=0;
+		}
+		calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-6);
+		calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-6);
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new PentaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+
+}
+/*}}}*/
 void       Penta::ComputeBasalStress(Vector<IssmDouble>* sigma_b){/*{{{*/
 
 	int         i,j;
@@ -2432,7 +2507,7 @@
 		strainxy=epsilon[3];
 
 		/*strainparallel= Strain rate along the ice flow direction */
-		strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
+		strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-14);
 	}
 
 	/*Add input*/
@@ -2479,7 +2554,7 @@
 		strainxy=epsilon[3];
 
 		/*strainperpendicular= Strain rate perpendicular to the ice flow direction */
-		strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
+		strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-14);
 	}
 
 	/*Add input*/
Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18952)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18953)
@@ -48,6 +48,7 @@
 		void           AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
 		IssmDouble     CharacteristicLength(void){_error_("not implemented yet");};
 		void           CalvingRateLevermann();
+		void           CalvingRatePi();
 		void           ComputeBasalStress(Vector<IssmDouble>* sigma_b);
 		void           ComputeDeviatoricStressTensor();
 		void           ComputeSigmaNN(){_error_("not implemented yet");};
Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Seg.h	(revision 18952)
+++ ../trunk-jpl/src/c/classes/Elements/Seg.h	(revision 18953)
@@ -48,6 +48,7 @@
 		void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
 		void        ComputeSigmaNN(){_error_("not implemented yet");};
 		void        ComputeStressTensor(){_error_("not implemented yet");};
+		void        CalvingRatePi(void){_error_("not implemented yet");};
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
 		void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
 		void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 18952)
+++ ../trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 18953)
@@ -48,6 +48,7 @@
 		void        ComputeSigmaNN(){_error_("not implemented yet");};
 		void        ComputeStressTensor(){_error_("not implemented yet");};
 		void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
+		void        CalvingRatePi(void){_error_("not implemented yet");};
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
 		void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
 		void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18952)
+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18953)
@@ -1368,7 +1368,18 @@
 			case CalvingCalvingrateEnum:
 				this->StrainRateparallel();
 				this->StrainRateperpendicular();
-				this->CalvingRateLevermann();
+				int calvinglaw;
+				this->FindParam(&calvinglaw,CalvingLawEnum);
+				switch(calvinglaw){
+					         case CalvingLevermannEnum:
+									this->CalvingRateLevermann();
+									break;
+								case CalvingPiEnum:
+									this->CalvingRatePi();
+									break;
+								default:
+									_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
+				}
 				input=this->inputs->GetInput(output_enum);
 				break;
 			case StrainRateparallelEnum:
Index: ../trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.h	(revision 18952)
+++ ../trunk-jpl/src/c/classes/FemModel.h	(revision 18953)
@@ -85,6 +85,7 @@
 		void StrainRateparallelx();
 		void StrainRateperpendicularx();
 		void CalvingRateLevermannx();
+		void CalvingRatePix();
 		#ifdef  _HAVE_DAKOTA_
 		void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);
 		#endif
Index: ../trunk-jpl/src/m/classes/model.m
===================================================================
--- ../trunk-jpl/src/m/classes/model.m	(revision 18952)
+++ ../trunk-jpl/src/m/classes/model.m	(revision 18953)
@@ -835,6 +835,7 @@
 			% Calving variables
 			if isa(md.calving,'calving'),md.calving.calvingrate=project3d(md,'vector',md.calving.calvingrate,'type','node');end;
 			if isa(md.calving,'calvinglevermann'),md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node');end;
+			if isa(md.calving,'calvingpi'),md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node');end;
 
 			% Hydrologydc variables
 			if isa(md.hydrology,'hydrologydc');
Index: ../trunk-jpl/src/m/classes/calvingpi.m
===================================================================
--- ../trunk-jpl/src/m/classes/calvingpi.m	(revision 0)
+++ ../trunk-jpl/src/m/classes/calvingpi.m	(revision 18953)
@@ -0,0 +1,53 @@
+%CALVINGPI class definition
+%
+%   Usage:
+%      calvingpi=calvingpi();
+
+classdef calvingpi
+	properties (SetAccess=public) 
+		 coeff       = NaN;
+		 meltingrate = NaN;
+	end
+	methods
+		function obj = calvingpi(varargin) % {{{
+			switch nargin
+				case 0
+					obj=setdefaultparameters(obj);
+				case 1
+					inputstruct=varargin{1};
+					list1 = properties('calvingpi');
+					list2 = fieldnames(inputstruct);
+					for i=1:length(list1)
+						fieldname = list1{i};
+						if ismember(fieldname,list2),
+							obj.(fieldname) = inputstruct.(fieldname);
+						end
+					end
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function obj = setdefaultparameters(obj) % {{{
+			%Proportionality coefficient in Pi model
+			obj.coeff=2e13;
+		end % }}}
+		function md = checkconsistency(obj,md,solution,analyses) % {{{
+			%Early return
+			if (solution~=TransientSolutionEnum() | md.transient.iscalving==0), return; end
+			md = checkfield(md,'fieldname','calving.coeff','>',0,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'size',[md.mesh.numberofvertices 1],'>=',0);
+		end % }}}
+		function disp(obj) % {{{
+			disp(sprintf('   Calving Pi parameters:'));
+			fielddisplay(obj,'coeff','proportionality coefficient in Pi model');
+			fielddisplay(obj,'meltingrate','melting rate at given location [m/a]');
+
+		end % }}}
+		function marshall(obj,md,fid) % {{{
+			yts=365.0*24.0*3600.0;
+			WriteData(fid,'enum',CalvingLawEnum(),'data',CalvingPiEnum(),'format','Integer');
+			WriteData(fid,'enum',CalvingpiCoeffEnum(),'data',obj.coeff,'format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts);
+		end % }}}
+	end
+end
Index: ../trunk-jpl/src/m/enum/CalvingpiCoeffEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/CalvingpiCoeffEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/CalvingpiCoeffEnum.m	(revision 18953)
@@ -0,0 +1,11 @@
+function macro=CalvingpiCoeffEnum()
+%CALVINGPICOEFFENUM - Enum of CalvingpiCoeff
+%
+%   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=CalvingpiCoeffEnum()
+
+macro=StringToEnum('CalvingpiCoeff');
Index: ../trunk-jpl/src/m/enum/CalvingPiEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/CalvingPiEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/CalvingPiEnum.m	(revision 18953)
@@ -0,0 +1,11 @@
+function macro=CalvingPiEnum()
+%CALVINGPIENUM - Enum of CalvingPi
+%
+%   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=CalvingPiEnum()
+
+macro=StringToEnum('CalvingPi');
Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18952)
+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 18953)
@@ -214,10 +214,13 @@
 def CalvingCalvingrateEnum(): return StringToEnum("CalvingCalvingrate")[0]
 def CalvingMeltingrateEnum(): return StringToEnum("CalvingMeltingrate")[0]
 def CalvingLevermannEnum(): return StringToEnum("CalvingLevermann")[0]
+def CalvingPiEnum(): return StringToEnum("CalvingPi")[0]
 def DefaultCalvingEnum(): return StringToEnum("DefaultCalving")[0]
 def CalvingRequestedOutputsEnum(): return StringToEnum("CalvingRequestedOutputs")[0]
 def CalvinglevermannCoeffEnum(): return StringToEnum("CalvinglevermannCoeff")[0]
 def CalvinglevermannMeltingrateEnum(): return StringToEnum("CalvinglevermannMeltingrate")[0]
+def CalvingpiCoeffEnum(): return StringToEnum("CalvingpiCoeff")[0]
+def CalvingpiMeltingrateEnum(): return StringToEnum("CalvingpiMeltingrate")[0]
 def CalvingratexEnum(): return StringToEnum("Calvingratex")[0]
 def CalvingrateyEnum(): return StringToEnum("Calvingratey")[0]
 def CalvingratexAverageEnum(): return StringToEnum("CalvingratexAverage")[0]
Index: ../trunk-jpl/src/m/enum/CalvingpiMeltingrateEnum.m
===================================================================
--- ../trunk-jpl/src/m/enum/CalvingpiMeltingrateEnum.m	(revision 0)
+++ ../trunk-jpl/src/m/enum/CalvingpiMeltingrateEnum.m	(revision 18953)
@@ -0,0 +1,11 @@
+function macro=CalvingpiMeltingrateEnum()
+%CALVINGPIMELTINGRATEENUM - Enum of CalvingpiMeltingrate
+%
+%   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=CalvingpiMeltingrateEnum()
+
+macro=StringToEnum('CalvingpiMeltingrate');
