Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18956)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18957)
@@ -1504,27 +1504,70 @@
 IssmDouble Tria::IceVolume(void){/*{{{*/
 
-	/*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
-	IssmDouble base,surface,bed;
-	IssmDouble xyz_list[NUMVERTICES][3];
-
-	if(!IsIceInElement())return 0;
-
-	/*First get back the area of the base*/
-	base=this->GetArea();
-
-	/*Now get the average height*/
-	Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
-	Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
-	surface_input->GetInputAverage(&surface);
-	base_input->GetInputAverage(&bed);
-
-	/*Return: */
+	/*The volume of a truncated prism is area_base * 1/numedges sum(length of edges)*/
+
+	/*Intermediaries*/
+	int i, numiceverts;
+	IssmDouble area_base,surface,base,Haverage;
+	IssmDouble Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
+	IssmDouble s[2]; // s:fraction of intersected triangle edges, that lies inside ice
+	int* indices=NULL;
+	IssmDouble* H=NULL;
+
+	if(!IsIceInElement())return 0.;
+
 	int domaintype;
 	parameters->FindParam(&domaintype,DomainTypeEnum);
+
+	if(IsIcefront()){
+		area_base=this->GetAreaIce();
+		//Assumption: linear ice thickness profile on element. 
+		//Hence ice thickness at intersection of levelset function with triangle edge is linear interpolation of ice thickness at vertices.
+		this->GetLevelsetIntersection(&indices, &numiceverts, s, MaskIceLevelsetEnum, 0.);
+		GetInputListOnVertices(&surfaces[0],SurfaceEnum);
+		GetInputListOnVertices(&bases[0],BaseEnum);
+		for(i=0;i<NUMVERTICES;i++) Haux[i]= surfaces[indices[i]]-bases[indices[i]]; //sort thicknesses in ice/noice
+		int numthk=numiceverts+2;
+		H=xNew<IssmDouble>(numthk);
+		switch(numiceverts){
+			case 1: // average over triangle 
+				H[0]=Haux[0];
+				H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
+				H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
+				break;
+			case 2: // average over quadrangle
+				H[0]=Haux[0];
+				H[1]=Haux[1];
+				H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
+				H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
+				break;
+			default:
+				_error_("Number of ice covered vertices wrong in Tria::IceVolume()");
+				break;
+		}
+		Haverage=0.;
+		for(i=0;i<numthk;i++)	Haverage+=H[i];
+		Haverage/=IssmDouble(numthk);
+	}
+	else{
+		/*First get back the area of the base*/
+		area_base=this->GetArea();
+
+		/*Now get the average height*/
+		Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
+		Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
+		surface_input->GetInputAverage(&surface);
+		base_input->GetInputAverage(&base);
+		Haverage=surface-base;
+	}
+
+	/*Cleanup & return: */
+	xDelete<int>(indices);
+	xDelete<IssmDouble>(H);
+
 	if(domaintype==Domain2DverticalEnum){
-	  return base;
+	  return area_base;
 	}
 	else{
-	  return base*(surface-bed);
+	  return area_base*Haverage;
 	}
 }