Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18841)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18842)
@@ -1315,6 +1315,48 @@
 	*pxyz_zero= xyz_zero;
 }
 /*}}}*/
+void       Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
+	
+	/* Intermediaries */
+	const int dim=3;
+	int i, dir,nrfrontnodes;
+	IssmDouble  levelset[NUMVERTICES];
+
+	/*Recover parameters and values*/
+	GetInputListOnVertices(&levelset[0],levelsetenum);
+
+	int* indicesfront = xNew<int>(NUMVERTICES);
+	/* Get basal nodes where there is no ice */
+	nrfrontnodes=0;
+	for(i=0;i<NUMVERTICES2D;i++){
+		if(levelset[i]>=0.){
+			indicesfront[nrfrontnodes]=i;
+			nrfrontnodes++;
+		}
+	}
+	_assert_(nrfrontnodes==2);
+
+	/* arrange order of basal frontnodes such that they are oriented counterclockwise */
+	if((NUMVERTICES2D+indicesfront[0]-indicesfront[1])%NUMVERTICES2D!=NUMVERTICES2D-1){
+		int index=indicesfront[0];
+		indicesfront[0]=indicesfront[1];
+		indicesfront[1]=index;
+	}	
+
+	IssmDouble* xyz_front = xNew<IssmDouble>(2*dim*nrfrontnodes);
+	/* Return basal and top front nodes */
+	for(i=0;i<nrfrontnodes;i++){
+		for(dir=0;dir<dim;dir++){
+			int ind1=i*dim+dir, ind2=(2*nrfrontnodes-1-i)*dim+dir; // vertex structure front segment: base0, base1, top1, top0
+			xyz_front[ind1]=xyz_list[dim*indicesfront[i]+dir];
+			xyz_front[ind2]=xyz_list[dim*(indicesfront[i]+NUMVERTICES2D)+dir];
+		}
+	}
+
+	*pxyz_front=xyz_front;
+
+	xDelete<int>(indicesfront);
+}/*}}}*/
 void       Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){/*{{{*/
 
 	int  step,i;
Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18841)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18842)
@@ -95,7 +95,7 @@
 		bool   IsIcefront(void);
 		bool   IsFaceOnBoundary(void){_error_("not implemented yet");};
 		void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
-		void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
+		void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
 		void   GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented yet");};
 		void   PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
 		void   ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
