Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 24047)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 24048)
@@ -3292,4 +3292,5 @@
 			}
 			break;
+		case CalvingFluxLevelsetEnum: this->CalvingFluxLevelset(); break;
 		case StrainRateparallelEnum: this->StrainRateparallel(); break;
 		case StrainRateperpendicularEnum: this->StrainRateperpendicular(); break;
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 24047)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 24048)
@@ -208,4 +208,5 @@
 		virtual void       CalvingCrevasseDepth(void){_error_("not implemented yet");};
 		virtual void	    CalvingRateLevermann(void)=0;
+		virtual void       CalvingFluxLevelset(void){_error_("not implemented yet");};
 		virtual IssmDouble CharacteristicLength(void)=0;
 		virtual void       ComputeBasalStress(void){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24047)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24048)
@@ -482,4 +482,138 @@
 	delete gauss;
 
+}
+/*}}}*/
+void       Tria::CalvingFluxLevelset(){/*{{{*/
+
+	/*Make sure there is an ice front here*/
+	if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)){
+		IssmDouble flux_per_area=0;
+		this->inputs->AddInput(new TriaInput(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum));
+	}
+	else{
+	int               domaintype,index1,index2;
+	const IssmPDouble epsilon = 1.e-15;
+	IssmDouble        s1,s2;
+	IssmDouble        gl[NUMVERTICES];
+	IssmDouble        xyz_front[2][3];
+
+
+	IssmDouble *xyz_list = NULL;
+	this->GetVerticesCoordinates(&xyz_list);
+
+	/*Recover parameters and values*/
+	parameters->FindParam(&domaintype,DomainTypeEnum);
+	GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+
+	/*Be sure that values are not zero*/
+	if(gl[0]==0.) gl[0]=gl[0]+epsilon;
+	if(gl[1]==0.) gl[1]=gl[1]+epsilon;
+	if(gl[2]==0.) gl[2]=gl[2]+epsilon;
+
+	if(domaintype==Domain2DverticalEnum){
+		_error_("not implemented");
+	}
+	else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
+		int pt1 = 0;
+		int pt2 = 1;
+		if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+
+			/*Portion of the segments*/
+			s1=gl[2]/(gl[2]-gl[1]);
+			s2=gl[2]/(gl[2]-gl[0]);
+			if(gl[2]<0.){
+				pt1 = 1; pt2 = 0;
+			}
+			xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
+			xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
+			xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
+			xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
+			xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
+			xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
+		}
+		else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+
+			/*Portion of the segments*/
+			s1=gl[0]/(gl[0]-gl[1]);
+			s2=gl[0]/(gl[0]-gl[2]);
+			if(gl[0]<0.){
+				pt1 = 1; pt2 = 0;
+			}
+
+			xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
+			xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
+			xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
+			xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
+			xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
+			xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
+		}
+		else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+
+			/*Portion of the segments*/
+			s1=gl[1]/(gl[1]-gl[0]);
+			s2=gl[1]/(gl[1]-gl[2]);
+			if(gl[1]<0.){
+				pt1 = 1; pt2 = 0;
+			}
+
+			xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
+			xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
+			xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
+			xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
+			xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
+			xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+		}
+		else{
+			_error_("case not possible");
+		}
+
+	}
+	else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
+
+	/*Some checks in debugging mode*/
+	_assert_(s1>=0 && s1<=1.); 
+	_assert_(s2>=0 && s2<=1.); 
+
+	/*Get normal vector*/
+	IssmDouble normal[3];
+	this->NormalSection(&normal[0],&xyz_front[0][0]);
+	normal[0] = -normal[0];
+	normal[1] = -normal[1];
+
+	/*Get inputs*/
+	IssmDouble flux = 0.;
+	IssmDouble area = 0.;
+	IssmDouble calvingratex,calvingratey,thickness,Jdet,flux_per_area;
+	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
+	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* calvingratex_input=NULL;
+	Input* calvingratey_input=NULL;
+	if(domaintype==Domain2DhorizontalEnum){
+		calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
+		calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
+	}
+	else{
+		calvingratex_input=inputs->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
+		calvingratey_input=inputs->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
+	}
+
+	/*Start looping on Gaussian points*/
+	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		thickness_input->GetInputValue(&thickness,gauss);
+		calvingratex_input->GetInputValue(&calvingratex,gauss);
+		calvingratey_input->GetInputValue(&calvingratey,gauss);
+		this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
+
+		flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
+		area += Jdet*gauss->weight*thickness; 
+
+		flux_per_area=flux/area;
+	}
+	
+	this->inputs->AddInput(new TriaInput(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum));
+	}
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 24047)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 24048)
@@ -54,4 +54,5 @@
 		void        CalvingCrevasseDepth();
 		void			CalvingRateLevermann();
+		void			CalvingFluxLevelset();
 		IssmDouble  CharacteristicLength(void);
 		void        ComputeBasalStress(void);
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 24047)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 24048)
@@ -1029,4 +1029,12 @@
 		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
 		element->CalvingRateLevermann();
+	}
+}
+/*}}}*/
+void FemModel::CalvingFluxLevelsetx(){/*{{{*/
+
+	for(int i=0;i<elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		element->CalvingFluxLevelset();
 	}
 }
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 24047)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 24048)
@@ -90,4 +90,5 @@
 		void CalvingRateVonmisesx();
 		void CalvingRateLevermannx();
+		void CalvingFluxLevelsetx();
 		void DeviatoricStressx();
 		void Divergencex(IssmDouble* pdiv);
Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 24047)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 24048)
@@ -482,4 +482,5 @@
 syn keyword cConstant CalvingrateyAverageEnum
 syn keyword cConstant CalvingrateyEnum
+syn keyword cConstant CalvingFluxLevelsetEnum
 syn keyword cConstant ConvergedEnum
 syn keyword cConstant CrevasseDepthEnum
@@ -1271,4 +1272,5 @@
 syn keyword cType Cfsurfacesquare
 syn keyword cType Channel
+syn keyword cType classes
 syn keyword cType Constraint
 syn keyword cType Constraints
@@ -1277,6 +1279,6 @@
 syn keyword cType ControlInput
 syn keyword cType Covertree
+syn keyword cType DatasetInput
 syn keyword cType DataSetParam
-syn keyword cType DatasetInput
 syn keyword cType Definition
 syn keyword cType DependentObject
@@ -1291,6 +1293,6 @@
 syn keyword cType ElementHook
 syn keyword cType ElementMatrix
+syn keyword cType Elements
 syn keyword cType ElementVector
-syn keyword cType Elements
 syn keyword cType ExponentialVariogram
 syn keyword cType ExternalResult
@@ -1299,9 +1301,10 @@
 syn keyword cType Friction
 syn keyword cType Gauss
+syn keyword cType GaussianVariogram
+syn keyword cType gaussobjects
 syn keyword cType GaussPenta
 syn keyword cType GaussSeg
 syn keyword cType GaussTetra
 syn keyword cType GaussTria
-syn keyword cType GaussianVariogram
 syn keyword cType GenericExternalResult
 syn keyword cType GenericOption
@@ -1318,4 +1321,5 @@
 syn keyword cType IssmDirectApplicInterface
 syn keyword cType IssmParallelDirectApplicInterface
+syn keyword cType krigingobjects
 syn keyword cType Load
 syn keyword cType Loads
@@ -1328,4 +1332,5 @@
 syn keyword cType Matice
 syn keyword cType Matlitho
+syn keyword cType matrixobjects
 syn keyword cType MatrixParam
 syn keyword cType Misfit
@@ -1340,6 +1345,6 @@
 syn keyword cType Observations
 syn keyword cType Option
+syn keyword cType Options
 syn keyword cType OptionUtilities
-syn keyword cType Options
 syn keyword cType Param
 syn keyword cType Parameters
@@ -1354,10 +1359,10 @@
 syn keyword cType Regionaloutput
 syn keyword cType Results
+syn keyword cType Riftfront
 syn keyword cType RiftStruct
-syn keyword cType Riftfront
 syn keyword cType Seg
 syn keyword cType SegInput
+syn keyword cType Segment
 syn keyword cType SegRef
-syn keyword cType Segment
 syn keyword cType SpcDynamic
 syn keyword cType SpcStatic
@@ -1379,8 +1384,4 @@
 syn keyword cType Vertex
 syn keyword cType Vertices
-syn keyword cType classes
-syn keyword cType gaussobjects
-syn keyword cType krigingobjects
-syn keyword cType matrixobjects
 syn keyword cType AdjointBalancethickness2Analysis
 syn keyword cType AdjointBalancethicknessAnalysis
@@ -1391,4 +1392,5 @@
 syn keyword cType BalancethicknessSoftAnalysis
 syn keyword cType BalancevelocityAnalysis
+syn keyword cType DamageCalvingAnalysis
 syn keyword cType DamageEvolutionAnalysis
 syn keyword cType DepthAverageAnalysis
@@ -1401,6 +1403,6 @@
 syn keyword cType FreeSurfaceBaseAnalysis
 syn keyword cType FreeSurfaceTopAnalysis
+syn keyword cType GiaIvinsAnalysis
 syn keyword cType GLheightadvectionAnalysis
-syn keyword cType GiaIvinsAnalysis
 syn keyword cType HydrologyDCEfficientAnalysis
 syn keyword cType HydrologyDCInefficientAnalysis
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 24047)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 24048)
@@ -478,4 +478,5 @@
 	CalvingrateyAverageEnum,
 	CalvingrateyEnum,
+	CalvingFluxLevelsetEnum,
 	ConvergedEnum,
 	CrevasseDepthEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 24047)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 24048)
@@ -484,4 +484,5 @@
 		case CalvingrateyAverageEnum : return "CalvingrateyAverage";
 		case CalvingrateyEnum : return "Calvingratey";
+		case CalvingFluxLevelsetEnum : return "CalvingFluxLevelset";
 		case ConvergedEnum : return "Converged";
 		case CrevasseDepthEnum : return "CrevasseDepth";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 24047)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 24048)
@@ -493,4 +493,5 @@
 	      else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum;
 	      else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
+	      else if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum;
 	      else if (strcmp(name,"Converged")==0) return ConvergedEnum;
 	      else if (strcmp(name,"CrevasseDepth")==0) return CrevasseDepthEnum;
@@ -505,9 +506,9 @@
 	      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 stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"DeviatoricStress1")==0) return DeviatoricStress1Enum;
+	      if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
+	      else if (strcmp(name,"DeviatoricStress1")==0) return DeviatoricStress1Enum;
 	      else if (strcmp(name,"DeviatoricStress2")==0) return DeviatoricStress2Enum;
 	      else if (strcmp(name,"DistanceToCalvingfront")==0) return DistanceToCalvingfrontEnum;
@@ -628,9 +629,9 @@
 	      else if (strcmp(name,"Misfit")==0) return MisfitEnum;
 	      else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
-	      else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"Node")==0) return NodeEnum;
+	      if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
+	      else if (strcmp(name,"Node")==0) return NodeEnum;
 	      else if (strcmp(name,"OmegaAbsGradient")==0) return OmegaAbsGradientEnum;
 	      else if (strcmp(name,"P0")==0) return P0Enum;
@@ -751,9 +752,9 @@
 	      else if (strcmp(name,"SmbVz")==0) return SmbVzEnum;
 	      else if (strcmp(name,"SmbW")==0) return SmbWEnum;
-	      else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum;
+	      if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
+	      else if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum;
 	      else if (strcmp(name,"SmbZMin")==0) return SmbZMinEnum;
 	      else if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum;
@@ -874,9 +875,9 @@
 	      else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum;
 	      else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum;
-	      else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
+	      if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
+	      else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
 	      else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
 	      else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum;
@@ -997,9 +998,9 @@
 	      else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
 	      else if (strcmp(name,"Element")==0) return ElementEnum;
-	      else if (strcmp(name,"ElementHook")==0) return ElementHookEnum;
          else stage=9;
    }
    if(stage==9){
-	      if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
+	      if (strcmp(name,"ElementHook")==0) return ElementHookEnum;
+	      else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
 	      else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
 	      else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum;
@@ -1120,9 +1121,9 @@
 	      else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
 	      else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;
-	      else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum;
          else stage=10;
    }
    if(stage==10){
-	      if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
+	      if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum;
+	      else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
 	      else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
 	      else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
@@ -1243,9 +1244,9 @@
 	      else if (strcmp(name,"Tetra")==0) return TetraEnum;
 	      else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
-	      else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
          else stage=11;
    }
    if(stage==11){
-	      if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
+	      if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
+	      else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
 	      else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum;
 	      else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum;