Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18955)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18956)
@@ -879,4 +879,37 @@
 }
 /*}}}*/
+IssmDouble Tria::GetAreaIce(void){/*{{{*/
+
+	/*return area of element covered by ice*/
+	/*Intermediaries*/
+	int numiceverts;
+	IssmDouble area_fraction;
+	IssmDouble s[2]; // s:fraction of intersected triangle edges that lie inside ice
+	int* indices=NULL;
+
+	this->GetLevelsetIntersection(&indices, &numiceverts, s, MaskIceLevelsetEnum, 0.);
+
+	switch (numiceverts){
+		case 0: // no vertex has ice: element is ice free
+			area_fraction=0.;
+			break;
+		case 1: // one vertex has ice: get area of triangle
+			area_fraction=s[0]*s[1];
+			break;
+		case 2: // two vertices have ice: get area of quadrangle
+			area_fraction=s[0]+s[1]-s[0]*s[1];
+			break;
+		case NUMVERTICES: // all vertices have ice: return triangle area
+			area_fraction=1.;
+			break;
+		default:
+			_error_("Wrong number of ice vertices in Tria::GetAreaIce!");
+			break;
+	}
+	_assert_((area_fraction>=0.) && (area_fraction<=1.));
+
+	xDelete<int>(indices);
+	return area_fraction*this->GetArea();
+}/*}}}*/
 void       Tria::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){/*{{{*/
 	/*Computeportion of the element that is grounded*/ 
@@ -1184,4 +1217,80 @@
 	xDelete<int>(indicesfront);
 }/*}}}*/
+void			Tria::GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level){/*{{{*/
+	
+	/* GetLevelsetIntersection computes: 
+	 * 1. indices of element, sorted in [iceverts, noiceverts] in counterclockwise fashion,
+	 * 2. fraction of intersected triangle edges intersected by levelset, lying below level*/
+
+	/*Intermediaries*/
+	int i, numiceverts, numnoiceverts;
+	int ind0, ind1, lastindex;
+	int indices_ice[NUMVERTICES],indices_noice[NUMVERTICES];
+	IssmDouble lsf[NUMVERTICES];
+	int* indices = xNew<int>(NUMVERTICES);
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&lsf[0],levelset_enum);
+
+	/* Determine distribution of ice over element.
+	 * Exploit: ice/no-ice parts are connected, so find starting vertex of segment*/
+	lastindex=0;
+	for(i=0;i<NUMVERTICES;i++){ // go backwards along vertices, and check for sign change
+		ind0=(NUMVERTICES-i)%NUMVERTICES;
+		ind1=(ind0-1+NUMVERTICES)%NUMVERTICES;
+		if((lsf[ind0]-level)*(lsf[ind1]-level)<=0.){ // levelset has been crossed, find last index belonging to segment
+			if(lsf[ind1]==level) //if levelset intersects 2nd vertex, choose this vertex as last
+				lastindex=ind1;
+			else
+				lastindex=ind0;
+			break;
+		}
+	}
+
+	numiceverts=0;
+	numnoiceverts=0;
+	for(i=0;i<NUMVERTICES;i++){
+		ind0=(lastindex+i)%NUMVERTICES;
+		if(lsf[i]<=level){
+			indices_ice[numiceverts]=i;
+			numiceverts++;
+		}
+		else{
+			indices_noice[numnoiceverts]=i;
+			numnoiceverts++;
+		}
+	}
+	//merge indices 
+	for(i=0;i<numiceverts;i++){indices[i]=indices_ice[i];}
+	for(i=0;i<numnoiceverts;i++){indices[numiceverts+i]=indices_noice[i];}
+
+	switch (numiceverts){
+		case 0: // no vertex has ice: element is ice free, no intersection
+			for(i=0;i<2;i++)
+				fraction[i]=0.;
+			break;
+		case 1: // one vertex has ice:
+			for(i=0;i<2;i++){
+				fraction[i]=(level-lsf[indices[0]])/(lsf[indices[numiceverts+i]]-lsf[indices[0]]);
+			}
+			break;
+		case 2: // two vertices have ice: get area of quadrangle
+			for(i=0;i<2;i++){
+				fraction[i]=(level-lsf[indices[i]])/(lsf[indices[numiceverts]]-lsf[indices[i]]);
+			}
+			break;
+		case NUMVERTICES: // all vertices have ice: return triangle area
+			for(i=0;i<2;i++)
+				fraction[i]=1.;
+			break;
+		default:
+			_error_("Wrong number of ice vertices in Tria::GetLevelsetIntersection!");
+			break;
+	}
+
+	*pindices=indices;
+	*pnumiceverts=numiceverts;
+}
+/*}}}*/
 void       Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/
 	
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18955)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18956)
@@ -143,5 +143,7 @@
 		void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
 		IssmDouble     GetArea(void);
+		IssmDouble 	GetAreaIce(void);
 		void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
+		void		GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level);
 		int            GetElementType(void);
 		void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);