Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 24367)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 24368)
@@ -1928,4 +1928,125 @@
 	/*Clean up and return*/
 	return groundedarea;
+}
+/*}}}*/
+IssmDouble Penta::GroundinglineMassFlux(bool scaled){/*{{{*/
+
+	/*Make sure there is a grounding line here*/
+	if(!IsOnBase()) return 0;
+	if(!IsIceInElement()) return 0;
+	if(!IsZeroLevelset(MaskGroundediceLevelsetEnum)) return 0;
+
+	/*Scaled not implemented yet...*/
+	_assert_(!scaled);
+
+	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*/
+	GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
+
+	/*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;
+
+	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");
+	}
+
+
+	/*Some checks in debugging mode*/
+	_assert_(s1>=0 && s1<=1.);
+	_assert_(s2>=0 && s2<=1.);
+
+	/*Get normal vector*/
+	IssmDouble normal[3];
+	this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
+	normal[0] = -normal[0];
+	normal[1] = -normal[1];
+
+	this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+	this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+
+	/*Get inputs*/
+	IssmDouble flux = 0.;
+	IssmDouble vx,vy,thickness,Jdet;
+	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
+	Input2 *thickness_input = this->GetInput2(ThicknessEnum); _assert_(thickness_input);
+	Input2 *vx_input        = this->GetInput2(VxAverageEnum); _assert_(vx_input);
+	Input2 *vy_input        = this->GetInput2(VyAverageEnum); _assert_(vy_input);
+
+	/*Start looping on Gaussian points*/
+	Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		thickness_input->GetInputValue(&thickness,gauss);
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
+
+		flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	return flux;
+
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 24367)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 24368)
@@ -104,4 +104,5 @@
 		void           GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
 		IssmDouble     GroundedArea(bool scaled);
+		IssmDouble     GroundinglineMassFlux(bool scaled);
 		IssmDouble		IcefrontMassFluxLevelset(bool scaled);
 		IssmDouble     IceVolume(bool scaled);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24367)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24368)
@@ -2872,5 +2872,5 @@
 IssmDouble Tria::GroundinglineMassFlux(bool scaled){/*{{{*/
 
-	/*Make sure there is an ice front here*/
+	/*Make sure there is a grounding line here*/
 	if(!IsIceInElement()) return 0;
 	if(!IsZeroLevelset(MaskGroundediceLevelsetEnum)) return 0;
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 24367)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 24368)
@@ -1620,4 +1620,16 @@
 }/*}}}*/
 void FemModel::GroundinglineMassFluxx(IssmDouble* pM, bool scaled){/*{{{*/
+
+	/*First we need to depth average the velocities*/
+	int domaintype;
+	this->parameters->FindParam(&domaintype,DomainTypeEnum);
+	this->parameters->SetParam(VxEnum,InputToDepthaverageInEnum);
+	this->parameters->SetParam(VxAverageEnum,InputToDepthaverageOutEnum);
+	depthaverage_core(this);
+	if(domaintype!=Domain2DverticalEnum){
+		this->parameters->SetParam(VyEnum,InputToDepthaverageInEnum);
+		this->parameters->SetParam(VyAverageEnum,InputToDepthaverageOutEnum);
+		depthaverage_core(this);
+	}
 
 	IssmDouble local_mass_flux = 0;