Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 23887)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 23888)
@@ -241,4 +241,5 @@
 		virtual IssmDouble IceVolumeAboveFloatation(bool scaled)=0;
 		virtual IssmDouble IcefrontMassFlux(bool scaled){_error_("not implemented");};
+		virtual IssmDouble GroundinglineMassFlux(bool scaled){_error_("not implemented");};
 		virtual void       InputDepthAverageAtBase(int enum_type,int average_enum_type)=0;
 		virtual void       InputExtrude(int input_enum,int start)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 23887)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 23888)
@@ -1080,4 +1080,7 @@
 	int nrfrontbed,numiceverts;
 
+	if(!this->IsOnBase()) return 0;
+	if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
+
 	/*Retrieve all inputs and parameters*/
 	GetInputListOnVertices(&bed[0],BedEnum);
@@ -1085,7 +1088,4 @@
 	GetInputListOnVertices(&bases[0],BaseEnum);
 	GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
-
-	if(!this->IsOnBase()) return 0;
-	if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
 
 	nrfrontbed=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 23887)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 23888)
@@ -1409,4 +1409,6 @@
 	int nrfrontbed,numiceverts;
 
+	if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
+
 	/*Retrieve all inputs and parameters*/
 	GetInputListOnVertices(&bed[0],BedEnum);
@@ -1414,6 +1416,4 @@
 	GetInputListOnVertices(&bases[0],BaseEnum);
 	GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
-
-	if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
 
 	nrfrontbed=0;
@@ -1945,4 +1945,260 @@
 		return false;
 	}
+}
+/*}}}*/
+IssmDouble Tria::IcefrontMassFlux(bool scaled){/*{{{*/
+
+	/*Make sure there is an ice front here*/
+	if(!IsIceInElement() || !IsIcefront()) 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*/
+	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 vx,vy,thickness,Jdet;
+	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
+	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	if(domaintype==Domain2DhorizontalEnum){
+		vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
+		vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
+	}
+	else{
+		vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
+		vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_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);
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
+
+		flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
+	}
+
+	return flux;
+}
+/*}}}*/
+IssmDouble Tria::GroundinglineMassFlux(bool scaled){/*{{{*/
+
+	/*Make sure there is an ice front here*/
+	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*/
+	parameters->FindParam(&domaintype,DomainTypeEnum);
+	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;
+
+	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]);
+
+	/*Get inputs*/
+	IssmDouble flux = 0.;
+	IssmDouble vx,vy,thickness,Jdet;
+	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
+	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	if(domaintype==Domain2DhorizontalEnum){
+		vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
+		vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
+	}
+	else{
+		vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
+		vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_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);
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
+
+		flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
+	}
+
+	return flux;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 23887)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 23888)
@@ -95,4 +95,6 @@
 		IssmDouble  IceVolume(bool scaled);
 		IssmDouble  IceVolumeAboveFloatation(bool scaled);
+		IssmDouble  IcefrontMassFlux(bool scaled);
+		IssmDouble  GroundinglineMassFlux(bool scaled);
 		void        Ismip6FloatingiceMeltingRate(void);
 		void        InputDepthAverageAtBase(int enum_type,int average_enum_type);
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23887)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23888)
@@ -1546,4 +1546,20 @@
 
 }/*}}}*/
+void FemModel::GroundinglineMassFluxx(IssmDouble* pM, bool scaled){/*{{{*/
+
+	IssmDouble local_mass_flux = 0;
+	IssmDouble total_mass_flux;
+
+	for(int i=0;i<this->elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+		local_mass_flux+=element->GroundinglineMassFlux(scaled);
+	}
+	ISSM_MPI_Reduce(&local_mass_flux,&total_mass_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
+	ISSM_MPI_Bcast(&total_mass_flux,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+	/*Assign output pointers: */
+	*pM=total_mass_flux;
+
+}/*}}}*/
 void FemModel::IceMassx(IssmDouble* pM, bool scaled){/*{{{*/
 
@@ -2206,4 +2222,5 @@
 					case GroundedAreaEnum:                   this->GroundedAreax(&double_result,false);             break;
 					case GroundedAreaScaledEnum:             this->GroundedAreax(&double_result,true);              break;
+					case GroundinglineMassFluxEnum:          this->GroundinglineMassFluxx(&double_result,false);    break;
 					case FloatingAreaEnum:                   this->FloatingAreax(&double_result,false);             break;
 					case FloatingAreaScaledEnum:             this->FloatingAreax(&double_result,true);              break;
@@ -2403,5 +2420,5 @@
 		case IceVolumeAboveFloatationEnum:       this->IceVolumeAboveFloatationx(responses, false); break;
 		case IceVolumeAboveFloatationScaledEnum: this->IceVolumeAboveFloatationx(responses, true); break;
-		case IcefrontMassFluxEnum:               this->IcefrontMassFluxx(responses, true); break;
+		case IcefrontMassFluxEnum:               this->IcefrontMassFluxx(responses, false); break;
 		case GroundedAreaEnum:                   this->GroundedAreax(responses, false); break;
 		case GroundedAreaScaledEnum:             this->GroundedAreax(responses, true); break;
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 23887)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 23888)
@@ -102,4 +102,5 @@
 		void IcefrontAreax();
 		void IcefrontMassFluxx(IssmDouble* presponse, bool scaled);
+		void GroundinglineMassFluxx(IssmDouble* presponse, bool scaled);
 		void IceMassx(IssmDouble* pV, bool scaled);
 		void IceVolumex(IssmDouble* pV, bool scaled);