Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 18735)
+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 18736)
@@ -104,8 +104,8 @@
 	IssmDouble kappa;
 	IssmDouble Jdet, dt, D_scalar;
 	IssmDouble h,hx,hy,hz;
-	IssmDouble vel, calvingrate;
-	IssmDouble norm_dlsf;
+	IssmDouble vel;
+//	IssmDouble norm_dlsf;
 	IssmDouble* xyz_list = NULL;
 
 	/*Get problem dimension*/
@@ -129,30 +129,37 @@
 	IssmDouble*    v        = xNew<IssmDouble>(dim);
 	IssmDouble*    w        = xNew<IssmDouble>(dim);
 	IssmDouble*    c        = xNew<IssmDouble>(dim);
-	IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
+	//IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
 
 	/*Retrieve all inputs and parameters*/
 	basalelement->GetVerticesCoordinates(&xyz_list);
 	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
 	Input* vx_input=NULL;
 	Input* vy_input=NULL;
+	Input* calvingratex_input=NULL;
+	Input* calvingratey_input=NULL;
 	if(domaintype==Domain2DhorizontalEnum){
 		vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
 		vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
+		calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+		calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
 	}
 	else{
 		if(dim==1){
 			vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
+			calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
 		}
 		if(dim==2){
 			vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
 			vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
+			calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
+			calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
 		}
 	}
 
-	Input* lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
-	Input* lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
-	Input* calvingrate_input  = basalelement->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
+//	Input* lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
+//	Input* lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
+	//Input* calvingrate_input  = basalelement->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
 	
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
@@ -177,18 +184,20 @@
 		GetBprime(Bprime,basalelement,xyz_list,gauss); 
 		vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity 
 		vy_input->GetInputValue(&v[1],gauss); 
-		lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
-		lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
-		calvingrate_input->GetInputValue(&calvingrate,gauss);
+		calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity 
+		calvingratey_input->GetInputValue(&c[1],gauss); 
+		//lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
+		//lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
+		//calvingrate_input->GetInputValue(&calvingrate,gauss);
 
-		norm_dlsf=0.;
-		for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
-		norm_dlsf=sqrt(norm_dlsf);
+		//norm_dlsf=0.;
+		//for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
+		//norm_dlsf=sqrt(norm_dlsf);
 
-		if(norm_dlsf>1.e-10)
-			for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
-		else
-			for(i=0;i<dim;i++) c[i]=0.;
+		//if(norm_dlsf>1.e-10)
+		//	for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
+		//else
+		//for(i=0;i<dim;i++) c[i]=0.;
 		
 		for(i=0;i<dim;i++) w[i]=v[i]-c[i];
 
@@ -205,9 +214,9 @@
 					&Ke->values[0],1);
 
 		/* Stabilization */
-		int stabilization=2;
+		int stabilization=1;
 		vel=0.;
-		for(i=0;i<dim;i++) vel+=v[i]*v[i];
+		for(i=0;i<dim;i++) vel+=w[i]*w[i];
 		vel=sqrt(vel)+1.e-14;
 		switch(stabilization){
 			case 0:
@@ -216,7 +225,7 @@
 			case 1:
 				/* Artificial Diffusion */
 				basalelement->ElementSizes(&hx,&hy,&hz);
-				h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) ); 
+				h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) ); 
 				kappa=h*vel/2.;
 				for(row=0;row<dim;row++)
 					for(col=0;col<dim;col++)
@@ -233,10 +242,10 @@
 			case 2:
 				/* Streamline Upwinding */
 				basalelement->ElementSizes(&hx,&hy,&hz);
-				h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) );
+				h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
 				for(row=0;row<dim;row++) 
 					for(col=0;col<dim;col++) 
-						D[row*dim+col] = D_scalar*h/(2.*vel)*v[row]*v[col];
+						D[row*dim+col] = D_scalar*h/(2.*vel)*w[row]*w[col];
 
 				TripleMultiply(Bprime,dim,numnodes,1,
 							D,dim,dim,0,
@@ -257,7 +266,7 @@
 	xDelete<IssmDouble>(v);
 	xDelete<IssmDouble>(w);
 	xDelete<IssmDouble>(c);
-	xDelete<IssmDouble>(dlsf);
+	//xDelete<IssmDouble>(dlsf);
 	delete gauss;
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return Ke;
Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 18735)
+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 18736)
@@ -47,8 +47,10 @@
 void DamageEvolutionAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
 	int finiteelement;
+	bool   islevelset;
 
 	iomodel->Constant(&finiteelement,DamageElementinterpEnum);
+	iomodel->Constant(&islevelset,TransientIslevelsetEnum);
 
 	/*Update elements: */
 	int counter=0;
@@ -68,6 +70,10 @@
 	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
 	iomodel->FetchDataToInput(elements,PressureEnum);
 
+	if(islevelset){
+		iomodel->FetchDataToInput(elements,IceMaskNodeActivationEnum);
+	}
+
 }/*}}}*/
 void DamageEvolutionAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
@@ -440,6 +446,11 @@
 }/*}}}*/
 void DamageEvolutionAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
+	bool islevelset;
+	femmodel->parameters->FindParam(&islevelset,TransientIslevelsetEnum);
+	if(islevelset){
+		SetActiveNodesLSMx(femmodel);
+	}
 	return;
 }/*}}}*/
 
Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18735)
+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 18736)
@@ -207,6 +207,10 @@
 	DamageEnum,
 	NewDamageEnum,
 	StressIntensityFactorEnum,
+	CalvingratexEnum,
+	CalvingrateyEnum,
+	CalvingratexAverageEnum,
+	CalvingrateyAverageEnum,
 	StrainRateparallelEnum,
 	StrainRateperpendicularEnum,
 	MaterialsRhoIceEnum,
@@ -244,6 +248,8 @@
 	Domain3DEnum,
 	MiscellaneousNameEnum, //FIXME: only used by qmu, should not be marshalled (already in queueing script)
 	MasstransportHydrostaticAdjustmentEnum,
+	MasstransportIscalvingrateEnum,
+	MasstransportLevermannCalvingCoeffEnum,
 	MasstransportIsfreesurfaceEnum,
 	MasstransportMinThicknessEnum,
 	MasstransportPenaltyFactorEnum,
Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18735)
+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 18736)
@@ -215,6 +215,10 @@
 		case DamageEnum : return "Damage";
 		case NewDamageEnum : return "NewDamage";
 		case StressIntensityFactorEnum : return "StressIntensityFactor";
+		case CalvingratexEnum : return "Calvingratex";
+		case CalvingrateyEnum : return "Calvingratey";
+		case CalvingratexAverageEnum : return "CalvingratexAverage";
+		case CalvingrateyAverageEnum : return "CalvingrateyAverage";
 		case StrainRateparallelEnum : return "StrainRateparallel";
 		case StrainRateperpendicularEnum : return "StrainRateperpendicular";
 		case MaterialsRhoIceEnum : return "MaterialsRhoIce";
@@ -252,6 +256,8 @@
 		case Domain3DEnum : return "Domain3D";
 		case MiscellaneousNameEnum : return "MiscellaneousName";
 		case MasstransportHydrostaticAdjustmentEnum : return "MasstransportHydrostaticAdjustment";
+		case MasstransportIscalvingrateEnum : return "MasstransportIscalvingrate";
+		case MasstransportLevermannCalvingCoeffEnum : return "MasstransportLevermannCalvingCoeff";
 		case MasstransportIsfreesurfaceEnum : return "MasstransportIsfreesurface";
 		case MasstransportMinThicknessEnum : return "MasstransportMinThickness";
 		case MasstransportPenaltyFactorEnum : return "MasstransportPenaltyFactor";
Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18735)
+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 18736)
@@ -218,6 +218,10 @@
 	      else if (strcmp(name,"Damage")==0) return DamageEnum;
 	      else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
 	      else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
+	      else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
+	      else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
+	      else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
+	      else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum;
 	      else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
 	      else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum;
 	      else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
@@ -255,14 +259,16 @@
 	      else if (strcmp(name,"Domain3D")==0) return Domain3DEnum;
 	      else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum;
 	      else if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum;
+         else stage=3;
+   }
+   if(stage==3){
+	      if (strcmp(name,"MasstransportIscalvingrate")==0) return MasstransportIscalvingrateEnum;
+	      else if (strcmp(name,"MasstransportLevermannCalvingCoeff")==0) return MasstransportLevermannCalvingCoeffEnum;
 	      else if (strcmp(name,"MasstransportIsfreesurface")==0) return MasstransportIsfreesurfaceEnum;
 	      else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum;
 	      else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
 	      else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum;
-         else stage=3;
-   }
-   if(stage==3){
-	      if (strcmp(name,"MasstransportCalvingrate")==0) return MasstransportCalvingrateEnum;
+	      else if (strcmp(name,"MasstransportCalvingrate")==0) return MasstransportCalvingrateEnum;
 	      else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
 	      else if (strcmp(name,"MasstransportVertexPairing")==0) return MasstransportVertexPairingEnum;
 	      else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum;
@@ -376,16 +382,16 @@
 	      else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
 	      else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
 	      else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
-	      else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
+         else stage=4;
+   }
+   if(stage==4){
+	      if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
 	      else if (strcmp(name,"BalancethicknessSoftAnalysis")==0) return BalancethicknessSoftAnalysisEnum;
 	      else if (strcmp(name,"BalancethicknessSoftSolution")==0) return BalancethicknessSoftSolutionEnum;
 	      else if (strcmp(name,"BalancevelocityAnalysis")==0) return BalancevelocityAnalysisEnum;
 	      else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
 	      else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
-         else stage=4;
-   }
-   if(stage==4){
-	      if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
+	      else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
 	      else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
 	      else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
 	      else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
@@ -499,16 +505,16 @@
 	      else if (strcmp(name,"Tria")==0) return TriaEnum;
 	      else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
 	      else if (strcmp(name,"Tetra")==0) return TetraEnum;
-	      else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
+         else stage=5;
+   }
+   if(stage==5){
+	      if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
 	      else if (strcmp(name,"Penta")==0) return PentaEnum;
 	      else if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
 	      else if (strcmp(name,"Vertex")==0) return VertexEnum;
 	      else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
 	      else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum;
-         else stage=5;
-   }
-   if(stage==5){
-	      if (strcmp(name,"Air")==0) return AirEnum;
+	      else if (strcmp(name,"Air")==0) return AirEnum;
 	      else if (strcmp(name,"Ice")==0) return IceEnum;
 	      else if (strcmp(name,"Melange")==0) return MelangeEnum;
 	      else if (strcmp(name,"Water")==0) return WaterEnum;
@@ -622,16 +628,16 @@
 	      else if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
 	      else if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
 	      else if (strcmp(name,"P1P1")==0) return P1P1Enum;
-	      else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
+         else stage=6;
+   }
+   if(stage==6){
+	      if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
 	      else if (strcmp(name,"MINI")==0) return MINIEnum;
 	      else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
 	      else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
 	      else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
 	      else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
-         else stage=6;
-   }
-   if(stage==6){
-	      if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
+	      else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
 	      else if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum;
 	      else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
 	      else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
@@ -745,16 +751,16 @@
 	      else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
 	      else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
 	      else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
-	      else if (strcmp(name,"Seaiceatm")==0) return SeaiceatmEnum;
+         else stage=7;
+   }
+   if(stage==7){
+	      if (strcmp(name,"Seaiceatm")==0) return SeaiceatmEnum;
 	      else if (strcmp(name,"Seaiceocean")==0) return SeaiceoceanEnum;
 	      else if (strcmp(name,"SeaiceThickness")==0) return SeaiceThicknessEnum;
 	      else if (strcmp(name,"SeaiceConcentration")==0) return SeaiceConcentrationEnum;
 	      else if (strcmp(name,"SeaiceMinConcentration")==0) return SeaiceMinConcentrationEnum;
 	      else if (strcmp(name,"SeaiceMinThickness")==0) return SeaiceMinThicknessEnum;
-         else stage=7;
-   }
-   if(stage==7){
-	      if (strcmp(name,"SeaiceMaxThickness")==0) return SeaiceMaxThicknessEnum;
+	      else if (strcmp(name,"SeaiceMaxThickness")==0) return SeaiceMaxThicknessEnum;
 	      else if (strcmp(name,"SeaiceSpcvx")==0) return SeaiceSpcvxEnum;
 	      else if (strcmp(name,"SeaiceSpcvy")==0) return SeaiceSpcvyEnum;
 	      else if (strcmp(name,"SeaiceCoriolisFactor")==0) return SeaiceCoriolisFactorEnum;
Index: ../trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp	(revision 18735)
+++ ../trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp	(revision 18736)
@@ -58,7 +58,7 @@
 					analysis_type==L2ProjectionBaseAnalysisEnum || 
 					analysis_type==BalancethicknessAnalysisEnum ||
 					analysis_type==HydrologyDCInefficientAnalysisEnum ||
-					analysis_type==DamageEvolutionAnalysisEnum || 
+					//analysis_type==DamageEvolutionAnalysisEnum || 
 					analysis_type==HydrologyDCEfficientAnalysisEnum ||
 					analysis_type==LevelsetAnalysisEnum ||
 					analysis_type==ExtrapolationAnalysisEnum
Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 18735)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 18736)
@@ -82,6 +82,8 @@
 		parameters->AddObject(iomodel->CopyConstantObject(TransientIslevelsetEnum));
 		parameters->AddObject(iomodel->CopyConstantObject(TransientIsdamageevolutionEnum));
 		parameters->AddObject(iomodel->CopyConstantObject(TransientIshydrologyEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(MasstransportIscalvingrateEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(MasstransportLevermannCalvingCoeffEnum));
 		parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
 		parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum));
 
Index: ../trunk-jpl/src/c/cores/transient_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/transient_core.cpp	(revision 18735)
+++ ../trunk-jpl/src/c/cores/transient_core.cpp	(revision 18736)
@@ -20,7 +20,7 @@
 
 	/*parameters: */
 	IssmDouble starttime,finaltime,dt,yts;
-	bool       isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology;
+	bool       isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology,iscalvingrate;
 	bool       save_results,dakota_analysis;
 	bool       time_adapt=false;
 	int        output_frequency;
@@ -52,6 +52,7 @@
 	femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
 	femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum);
 	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
+	femmodel->parameters->FindParam(&iscalvingrate,MasstransportIscalvingrateEnum);
 	if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
 	femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
 	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum);
@@ -100,6 +101,12 @@
 		}
 
 		if(islevelset){
+			if(iscalvingrate){
+				if(VerboseSolution()) _printf0_("   computing calving rate\n");
+				femmodel->StrainRateparallelx();
+				femmodel->StrainRateperpendicularx();
+				femmodel->CalvingRateLevermannx();
+			}
 			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 18735)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18736)
@@ -1694,6 +1694,30 @@
 	}
 }
 	/*}}}*/
+void FemModel::CalvingRateLevermannx(){/*{{{*/
+
+	   for(int i=0;i<elements->Size();i++){
+			      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+					      element->CalvingRateLevermann();
+							   }
+}
+/*}}}*/
+void FemModel::StrainRateparallelx(){/*{{{*/
+
+	   for(int i=0;i<elements->Size();i++){
+			      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+					      element->StrainRateparallel();
+							   }
+}
+/*}}}*/
+void FemModel::StrainRateperpendicularx(){/*{{{*/
+
+	   for(int i=0;i<elements->Size();i++){
+			      Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+					      element->StrainRateperpendicular();
+							   }
+}
+/*}}}*/
 #ifdef  _HAVE_DAKOTA_
 void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/
 
Index: ../trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 18735)
+++ ../trunk-jpl/src/c/classes/Elements/Element.h	(revision 18736)
@@ -220,7 +220,8 @@
 		virtual void   ComputeDeviatoricStressTensor(void)=0;
 		virtual void	StressIntensityFactor(void)=0;
 		virtual void	StrainRateparallel(void)=0;
-		virtual void   StrainRateperpendicular(void)=0;
+		virtual void	StrainRateperpendicular(void)=0;
+		virtual void	CalvingRateLevermann(void)=0;
 
 		virtual void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0;
 		virtual void   InputExtrude(int input_enum,int start)=0;
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18735)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18736)
@@ -275,29 +275,76 @@
 	delete gauss;
 }
 /*}}}*/
-void       Tria::StrainRateperpendicular(){/*{{{*/
+void       Tria::ComputeDeviatoricStressTensor(){/*{{{*/
 
 	IssmDouble  xyz_list[NUMVERTICES][3];
-	GaussPenta* gauss=NULL;
+	IssmDouble  viscosity;
+	IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+	IssmDouble  tau_xx[NUMVERTICES];
+	IssmDouble	tau_yy[NUMVERTICES];
+	IssmDouble	tau_zz[NUMVERTICES]={0,0,0};
+	IssmDouble  tau_xy[NUMVERTICES];
+	IssmDouble	tau_xz[NUMVERTICES]={0,0,0};
+	IssmDouble	tau_yz[NUMVERTICES]={0,0,0};
+	GaussTria*  gauss=NULL;
+	int domaintype,dim=2;
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*Retrieve all inputs we will be needing: */
+	this->FindParam(&domaintype,DomainTypeEnum);
+	if(domaintype!=Domain2DhorizontalEnum) _error_("deviatoric stress tensor calculation not implemented for mesh of type " <<EnumToStringx(domaintype));
+	Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussTria();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/*Compute strain rate and viscosity: */
+		this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+		this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
+
+		/*Compute Stress*/
+		tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps
+		tau_yy[iv]=2*viscosity*epsilon[1];
+		tau_xy[iv]=2*viscosity*epsilon[2];
+	}
+
+	/*Add Stress tensor components into inputs*/
+	this->inputs->AddInput(new TriaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+}
+/*}}}*/
+void       Tria::StrainRateparallel(){/*{{{*/
+
+	IssmDouble *xyz_list = NULL;
+	IssmDouble  epsilon[3];
+	GaussTria* gauss=NULL;
 	IssmDouble  vx,vy,vel;
 	IssmDouble  strainxx;
 	IssmDouble  strainxy;
 	IssmDouble  strainyy;
-	IssmDouble  strainperpendicular[NUMVERTICES];
+	IssmDouble  strainparallel[NUMVERTICES];
 
 	/* Get node coordinates and dof list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	this->GetVerticesCoordinates(&xyz_list);
 
 	/*Retrieve all inputs we will need*/
 	Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
-	Input* strainxx_input=inputs->GetInput(StrainRatexxEnum);             _assert_(strainxx_input);
-	Input* strainxy_input=inputs->GetInput(StrainRatexyEnum);             _assert_(strainxy_input);
-	Input* strainyy_input=inputs->GetInput(StrainRateyyEnum);             _assert_(strainyy_input);
 
-
 	/* Start looping on the number of vertices: */
-	gauss=new GaussPenta();
+	gauss=new GaussTria();
 	for (int iv=0;iv<NUMVERTICES;iv++){
 		gauss->GaussVertex(iv);
 
@@ -305,43 +352,45 @@
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
 		vel=vx*vx+vy*vy;
-		strainxx_input->GetInputValue(&strainxx,gauss);
-		strainxy_input->GetInputValue(&strainxy,gauss);
-		strainyy_input->GetInputValue(&strainyy,gauss);
 
-		/*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);
+		/*Compute strain rate viscosity and pressure: */
+		this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+		strainxx=epsilon[0];
+		strainyy=epsilon[1];
+		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);
 	}
 
 	/*Add input*/
-	this->inputs->AddInput(new PentaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
 
 	/*Clean up and return*/
 	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
 }
 /*}}}*/
-void       Tria::StrainRateparallel(){/*{{{*/
+void       Tria::StrainRateperpendicular(){/*{{{*/
 
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	GaussPenta* gauss=NULL;
+	IssmDouble *xyz_list = NULL;
+	GaussTria* gauss=NULL;
+	IssmDouble  epsilon[3];
 	IssmDouble  vx,vy,vel;
 	IssmDouble  strainxx;
 	IssmDouble  strainxy;
 	IssmDouble  strainyy;
-	IssmDouble  strainparallel[NUMVERTICES];
+	IssmDouble  strainperpendicular[NUMVERTICES];
 
 	/* Get node coordinates and dof list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	
+	this->GetVerticesCoordinates(&xyz_list);
+
 	/*Retrieve all inputs we will need*/
 	Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
-	Input* strainxx_input=inputs->GetInput(StrainRatexxEnum);             _assert_(strainxx_input);
-	Input* strainxy_input=inputs->GetInput(StrainRatexyEnum);             _assert_(strainxy_input);
-	Input* strainyy_input=inputs->GetInput(StrainRateyyEnum);             _assert_(strainyy_input);
 
 	/* Start looping on the number of vertices: */
-	gauss=new GaussPenta();
+	gauss=new GaussTria();
 	for (int iv=0;iv<NUMVERTICES;iv++){
 		gauss->GaussVertex(iv);
 
@@ -349,69 +398,77 @@
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
 		vel=vx*vx+vy*vy;
-		strainxx_input->GetInputValue(&strainxx,gauss);
-		strainxy_input->GetInputValue(&strainxy,gauss);
-		strainyy_input->GetInputValue(&strainyy,gauss);
 
-		/*strainparallel= Strain rate along the ice flow direction */
-		strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
+		/*Compute strain rate viscosity and pressure: */
+		this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+		strainxx=epsilon[0];
+		strainyy=epsilon[1];
+		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);
 	}
 
 	/*Add input*/
-	this->inputs->AddInput(new PentaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
 
 	/*Clean up and return*/
 	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
 }
 /*}}}*/
-void       Tria::ComputeDeviatoricStressTensor(){/*{{{*/
+void       Tria::CalvingRateLevermann(){/*{{{*/
 
 	IssmDouble  xyz_list[NUMVERTICES][3];
-	IssmDouble  viscosity;
-	IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
-	IssmDouble  tau_xx[NUMVERTICES];
-	IssmDouble	tau_yy[NUMVERTICES];
-	IssmDouble	tau_zz[NUMVERTICES]={0,0,0};
-	IssmDouble  tau_xy[NUMVERTICES];
-	IssmDouble	tau_xz[NUMVERTICES]={0,0,0};
-	IssmDouble	tau_yz[NUMVERTICES]={0,0,0};
-	GaussTria*  gauss=NULL;
-	int domaintype,dim=2;
+	GaussTria* gauss=NULL;
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainparallel;
+	IssmDouble  propcoeff;
+	IssmDouble  strainperpendicular;
+	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 we will be needing: */
-	this->FindParam(&domaintype,DomainTypeEnum);
-	if(domaintype!=Domain2DhorizontalEnum) _error_("deviatoric stress tensor calculation not implemented for mesh of type " <<EnumToStringx(domaintype));
-	Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
-	Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
+	/*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* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);					_assert_(strainperpendicular_input);
+	this->parameters->FindParam(&propcoeff,MasstransportLevermannCalvingCoeffEnum);
 
 	/* Start looping on the number of vertices: */
 	gauss=new GaussTria();
 	for (int iv=0;iv<NUMVERTICES;iv++){
 		gauss->GaussVertex(iv);
 
-		/*Compute strain rate and viscosity: */
-		this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
-		this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
+		/* 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);
+		strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
 
-		/*Compute Stress*/
-		tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps
-		tau_yy[iv]=2*viscosity*epsilon[1];
-		tau_xy[iv]=2*viscosity*epsilon[2];
+		/*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
+		calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
+		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);
 	}
 
-	/*Add Stress tensor components into inputs*/
-	this->inputs->AddInput(new TriaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
+	/*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(MasstransportCalvingrateEnum,&calvingrate[0],P1Enum));
 
 	/*Clean up and return*/
 	delete gauss;
+
 }
 /*}}}*/
 void       Tria::Configure(Elements* elementsin, Loads* loadsin,Nodes* nodesin,Vertices *verticesin,Materials* materialsin, Parameters* parametersin){/*{{{*/
Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18735)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18736)
@@ -54,10 +54,11 @@
 		void        ComputeSigmaNN();
 		void        ComputeStressTensor();
 		void        ComputeDeviatoricStressTensor();
+		void			StrainRateparallel();
+		void			StrainRateperpendicular();
 		void        ComputeSurfaceNormalVelocity();
 		void        StressIntensityFactor(void){_error_("not implemented yet");};
-		void			StrainRateparallel();
-		void			StrainRateperpendicular();
+		void			CalvingRateLevermann();
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
 		void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
 		void        ResetHooks();
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18735)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18736)
@@ -346,6 +346,60 @@
 	delete gauss;
 }
 /*}}}*/
+void       Penta::CalvingRateLevermann(){/*{{{*/
+
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	GaussPenta* gauss=NULL;
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainparallel;
+	IssmDouble  propcoeff;
+	IssmDouble  strainperpendicular;
+	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* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);             _assert_(strainperpendicular_input);
+	this->parameters->FindParam(&propcoeff,MasstransportLevermannCalvingCoeffEnum);
+
+	/* 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);
+		strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
+
+		/*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
+		calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;	
+		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);
+	}
+
+	/*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(MasstransportCalvingrateEnum,&calvingrate[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+
+}
+/*}}}*/
 void       Penta::StressIntensityFactor(){/*{{{*/
 
 	/* Check if we are on the base */
@@ -425,7 +479,8 @@
 /*}}}*/
 void       Penta::StrainRateparallel(){/*{{{*/
 
-	IssmDouble  xyz_list[NUMVERTICES][3];
+	IssmDouble *xyz_list = NULL;
+	IssmDouble  epsilon[6];
 	GaussPenta* gauss=NULL;
 	IssmDouble  vx,vy,vel;
 	IssmDouble  strainxx;
@@ -434,14 +489,12 @@
 	IssmDouble  strainparallel[NUMVERTICES];
 
 	/* Get node coordinates and dof list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	
+	this->GetVerticesCoordinates(&xyz_list);
+
 	/*Retrieve all inputs we will need*/
 	Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
-	Input* strainxx_input=inputs->GetInput(StrainRatexxEnum);             _assert_(strainxx_input);
-	Input* strainxy_input=inputs->GetInput(StrainRatexyEnum);             _assert_(strainxy_input);
-	Input* strainyy_input=inputs->GetInput(StrainRateyyEnum);             _assert_(strainyy_input);
+	Input* vz_input=inputs->GetInput(VzEnum);												_assert_(vz_input);
 
 	/* Start looping on the number of vertices: */
 	gauss=new GaussPenta();
@@ -452,10 +505,13 @@
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
 		vel=vx*vx+vy*vy;
-		strainxx_input->GetInputValue(&strainxx,gauss);
-		strainxy_input->GetInputValue(&strainxy,gauss);
-		strainyy_input->GetInputValue(&strainyy,gauss);
 
+		/*Compute strain rate and viscosity: */
+		this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
+		strainxx=epsilon[0];
+		strainyy=epsilon[1];
+		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);
 	}
@@ -465,11 +521,13 @@
 
 	/*Clean up and return*/
 	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
 }
 /*}}}*/
 void       Penta::StrainRateperpendicular(){/*{{{*/
 
-	IssmDouble  xyz_list[NUMVERTICES][3];
+	IssmDouble *xyz_list = NULL;
+	IssmDouble  epsilon[6];
 	GaussPenta* gauss=NULL;
 	IssmDouble  vx,vy,vel;
 	IssmDouble  strainxx;
@@ -478,16 +536,13 @@
 	IssmDouble  strainperpendicular[NUMVERTICES];
 
 	/* Get node coordinates and dof list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	this->GetVerticesCoordinates(&xyz_list);
 
 	/*Retrieve all inputs we will need*/
 	Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
-	Input* strainxx_input=inputs->GetInput(StrainRatexxEnum);             _assert_(strainxx_input);
-	Input* strainxy_input=inputs->GetInput(StrainRatexyEnum);             _assert_(strainxy_input);
-	Input* strainyy_input=inputs->GetInput(StrainRateyyEnum);             _assert_(strainyy_input);
+	Input* vz_input=inputs->GetInput(VzEnum);												_assert_(vz_input);
 
-
 	/* Start looping on the number of vertices: */
 	gauss=new GaussPenta();
 	for (int iv=0;iv<NUMVERTICES;iv++){
@@ -497,10 +552,13 @@
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
 		vel=vx*vx+vy*vy;
-		strainxx_input->GetInputValue(&strainxx,gauss);
-		strainxy_input->GetInputValue(&strainxy,gauss);
-		strainyy_input->GetInputValue(&strainyy,gauss);
 
+		/*Compute strain rate and viscosity: */
+		this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
+		strainxx=epsilon[0];
+		strainyy=epsilon[1];
+		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);
 	}
@@ -510,6 +568,7 @@
 
 	/*Clean up and return*/
 	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
 }
 /*}}}*/
 void       Penta::Configure(Elements* elementsin, Loads* loadsin, Nodes* nodesin,Vertices* verticesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
@@ -2175,11 +2234,15 @@
 	if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
 	if(this->inputs->GetInput(VxEnum)) this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
 	if(this->inputs->GetInput(VyEnum)) this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+	if(this->inputs->GetInput(CalvingratexEnum)) this->InputDepthAverageAtBase(CalvingratexEnum,CalvingratexAverageEnum);
+	if(this->inputs->GetInput(CalvingrateyEnum)) this->InputDepthAverageAtBase(CalvingrateyEnum,CalvingrateyAverageEnum);
 	Tria* tria=(Tria*)SpawnTria(0,1,2);
 	this->inputs->DeleteInput(MaterialsRheologyBbarEnum);
 	this->inputs->DeleteInput(DamageDbarEnum);
 	this->inputs->DeleteInput(VxAverageEnum);
 	this->inputs->DeleteInput(VyAverageEnum);
+	this->inputs->DeleteInput(CalvingratexAverageEnum);
+	this->inputs->DeleteInput(CalvingrateyAverageEnum);
 
 	return tria;
 }
Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18735)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18736)
@@ -59,6 +59,7 @@
 		void   StressIntensityFactor();
 		void   StrainRateparallel();
 		void   StrainRateperpendicular();
+		void   CalvingRateLevermann();
 		void   Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
 		void   ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
 		int    FiniteElement(void);
Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Seg.h	(revision 18735)
+++ ../trunk-jpl/src/c/classes/Elements/Seg.h	(revision 18736)
@@ -57,8 +57,9 @@
 		void        ComputeStressTensor(){_error_("not implemented yet");};
 		void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
 		void        StressIntensityFactor(void){_error_("not implemented yet");};
-		void			StrainRateparallel(void){_error_("not implemented yet");};
-		void			StrainRateperpendicular(void){_error_("not implemented yet");};
+		void        StrainRateparallel(void){_error_("not implemented yet");};
+		void        StrainRateperpendicular(void){_error_("not implemented yet");};
+		void        CalvingRateLevermann(void){_error_("not implemented yet");};
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
 		void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
 		void        ResetHooks(){_error_("not implemented yet");};
Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 18735)
+++ ../trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 18736)
@@ -59,6 +59,7 @@
 		void        StressIntensityFactor(void){_error_("not implemented yet");};
 		void        StrainRateparallel(void){_error_("not implemented yet");};
 		void        StrainRateperpendicular(void){_error_("not implemented yet");};
+		void        CalvingRateLevermann(void){_error_("not implemented yet");};
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
 		void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
 		void        ResetHooks();
Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18735)
+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18736)
@@ -1209,6 +1209,14 @@
 				this->StressIntensityFactor();
 				input=this->inputs->GetInput(output_enum);
 				break;
+			case CalvingratexEnum:
+			case CalvingrateyEnum:
+			case MasstransportCalvingrateEnum:
+				this->StrainRateparallel();
+				this->StrainRateperpendicular();
+				this->CalvingRateLevermann();
+				input=this->inputs->GetInput(output_enum);
+				break;
 			case StrainRateparallelEnum:
 				this->StrainRateparallel();
 				input=this->inputs->GetInput(output_enum);
Index: ../trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.h	(revision 18735)
+++ ../trunk-jpl/src/c/classes/FemModel.h	(revision 18736)
@@ -81,6 +81,9 @@
 		void ElementResponsex(IssmDouble* presponse,int response_enum);
 		void BalancethicknessMisfitx(IssmDouble* pV);
 		void StressIntensityFactorx();
+		void StrainRateparallelx();
+		void StrainRateperpendicularx();
+		void CalvingRateLevermannx();
 		#ifdef  _HAVE_DAKOTA_
 		void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);
 		#endif
