Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 26625)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 26626)
@@ -1236,4 +1236,89 @@
 	/*return PentaRef field*/
 	return this->element_type;
+}
+/*}}}*/
+void       Penta::GetFractionGeometry2D(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){/*{{{*/
+
+  /*Compute portion of element that is grounded based on levelset at the 3 lower vertices of the Penta element*/
+   bool               trapezeisnegative=true;
+   int                point;
+   const IssmPDouble  epsilon= 1.e-15;
+   IssmDouble         f1,f2,phi;
+
+   /*Weights*/
+   Gauss* gauss = NULL;
+   IssmDouble loadweights_g[NUMVERTICES2D];
+   IssmDouble total_weight = 0;
+
+   _assert_(!xIsNan<IssmDouble>(gl[0]));
+   _assert_(!xIsNan<IssmDouble>(gl[1]));
+   _assert_(!xIsNan<IssmDouble>(gl[2]));
+
+   /*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;
+
+   /*Check that not all nodes are positive or negative: */
+   if(gl[0]>0 && gl[1]>0 && gl[2]>0){
+      point = 0;
+      f1    = 1.;
+      f2    = 1.;
+   }
+   else if(gl[0]<0 && gl[1]<0 && gl[2]<0){
+      point = 0;
+      f1    = 0.;
+      f2    = 0.;
+   }
+	else{
+		if(gl[0]*gl[1]*gl[2]<0) trapezeisnegative = false;
+
+		/*Find the similar nodes*/
+		if(gl[0]*gl[1]>0){ 
+			point = 2;
+			f1    = gl[2]/(gl[2]-gl[0]);
+			f2    = gl[2]/(gl[2]-gl[1]);
+		}
+		else if(gl[1]*gl[2]>0){ 
+			point = 0;
+			f1    = gl[0]/(gl[0]-gl[1]);
+			f2    = gl[0]/(gl[0]-gl[2]);
+		}
+		else if(gl[0]*gl[2]>0){ 
+			point = 1;
+			f1    = gl[1]/(gl[1]-gl[2]);
+			f2    = gl[1]/(gl[1]-gl[0]);
+		}
+		else _error_("case not possible");
+	}
+	if(trapezeisnegative) phi = 1-f1*f2;
+	else                  phi = f1*f2;
+	
+	/*Compute weights*/
+	gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2);
+
+	total_weight = 0.0;
+	for(int i=0;i<NUMVERTICES2D;i++)weights[i] = 0;
+	while(gauss->next()){
+		GetNodalFunctions(&loadweights_g[0],gauss,P1Enum);
+		for(int i=0;i<NUMVERTICES2D;i++)weights[i] += loadweights_g[i]*gauss->weight;
+		total_weight += gauss->weight;
+	}
+
+	/*Normalize to phi*/
+	if(total_weight>0.){
+		for(int i=0;i<NUMVERTICES2D;i++) weights[i] = weights[i]*phi/total_weight;
+	}
+	else for(int i=0;i<NUMVERTICES2D;i++) weights[i] = 0.0;
+
+	/*Cleanup*/
+	delete gauss;
+	
+	/*Assign output pointers*/
+	*pphi               = phi;
+	*ppoint1            = point;
+	*pfraction1         = f1;
+	*pfraction2         = f2;
+	*ptrapezeisnegative = trapezeisnegative;
 }
 /*}}}*/
@@ -2099,4 +2184,5 @@
 	IssmDouble base,height,scalefactor;
 	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble lsf[NUMVERTICES];
 
 	if(!IsIceInElement())return 0;
@@ -2104,17 +2190,53 @@
 	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
-	/*First calculate the area of the base (cross section triangle)
-	 * http://en.wikipedia.org/wiki/Pentangle
-	 * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
-	base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
-
-	if(scaled==true){ //scale for area projection correction
-		Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
-		scalefactor_input->GetInputAverage(&scalefactor);
-		base=base*scalefactor;
-	}
-
-	/*Now get the average height*/
-	height = 1./3.*((xyz_list[3][2]-xyz_list[0][2])+(xyz_list[4][2]-xyz_list[1][2])+(xyz_list[5][2]-xyz_list[2][2]));
+	Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+	/*Deal with partially ice-covered elements*/
+	if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
+		bool        istrapneg;
+		int         point;
+		IssmDouble  f1,f2,phi;
+		IssmDouble* heights = xNew<IssmDouble>(NUMVERTICES2D);
+		IssmDouble  weights[NUMVERTICES2D];
+		IssmDouble  lsf2d[NUMVERTICES2D];
+		for(int i=0;i<NUMVERTICES2D;i++){
+			heights[i] = xyz_list[i+NUMVERTICES2D][2]-xyz_list[i][2];
+			lsf2d[i]   = lsf[i];
+		}
+		GetFractionGeometry2D(weights,&phi,&point,&f1,&f2,&istrapneg,lsf2d);
+		
+		IssmDouble basetot;
+		height = 0.0;
+		for(int i=0;i<NUMVERTICES2D;i++) height += weights[i]*heights[i];
+		basetot = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
+		base    = basetot*phi;	
+	
+		/*Account for scaling factor averaged over subelement 2D area*/
+		if(scaled==true){
+			IssmDouble* scalefactor_vertices   = xNew<IssmDouble>(NUMVERTICES);
+			Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
+			/*Compute loop only over lower vertices: i<NUMVERTICES2D*/
+			scalefactor = 0.0;
+			for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]*scalefactor_vertices[i];
+			base = base*scalefactor;
+			xDelete<IssmDouble>(scalefactor_vertices);
+		}
+		xDelete<IssmDouble>(heights); 
+	}
+
+	else{ 
+		/*First calculate the area of the base (cross section triangle)
+		 * http://en.wikipedia.org/wiki/Pentangle
+		 * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
+		base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
+	
+		if(scaled==true){ //scale for area projection correction
+			Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
+			scalefactor_input->GetInputAverage(&scalefactor);
+			base=base*scalefactor;
+		}
+	
+		/*Now get the average height*/
+		height = 1./3.*((xyz_list[3][2]-xyz_list[0][2])+(xyz_list[4][2]-xyz_list[1][2])+(xyz_list[5][2]-xyz_list[2][2]));
+	}
 
 	/*Return: */
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 26625)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 26626)
@@ -86,4 +86,5 @@
 		IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
 		void           GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){_error_("not implemented yet");};
+		void           GetFractionGeometry2D(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl);
 		void       GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area,  int levelsetenum){_error_("not implemented yet");};
 		void       GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelset1enum, int levelset2enum){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26625)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26626)
@@ -3489,4 +3489,5 @@
 	parameters->FindParam(&domaintype,DomainTypeEnum);
 
+	/*Relict code
 	if(false && IsIcefront()){
 		//Assumption: linear ice thickness profile on element.
@@ -3543,5 +3544,44 @@
 		for(i=0;i<numthk;i++)	Haverage+=H[i];
 		Haverage/=IssmDouble(numthk);
-	}
+	}*/
+
+   IssmDouble* lsf = xNew<IssmDouble>(NUMVERTICES);
+   Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+   /*Deal with partially ice-covered elements*/
+	if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
+		bool istrapneg;
+      int point;
+      IssmDouble* weights  = xNew<IssmDouble>(NUMVERTICES);
+      IssmDouble* surfaces = xNew<IssmDouble>(NUMVERTICES);
+      IssmDouble* bases    = xNew<IssmDouble>(NUMVERTICES);
+      IssmDouble* Hice     = xNew<IssmDouble>(NUMVERTICES);
+      IssmDouble area_basetot,f1,f2,phi;
+      /*Average thickness over subelement*/
+		Element::GetInputListOnVertices(&surfaces[0],SurfaceEnum);
+      Element::GetInputListOnVertices(&bases[0],BaseEnum);
+      GetFractionGeometry(weights,&phi,&point,&f1,&f2,&istrapneg,lsf);
+      for(int i=0;i<NUMVERTICES;i++) Hice[i] = surfaces[i]-bases[i];
+      Haverage = 0.0;
+      for(int i=0;i<NUMVERTICES;i++) Haverage += weights[i]*Hice[i];
+		/*Get back area of ice-covered base*/
+		area_basetot = this->GetArea();
+		area_base    = phi*area_basetot;
+
+		/*Account for scaling factor averaged over subelement*/
+		if(scaled==true){
+			IssmDouble* scalefactor_vertices = xNew<IssmDouble>(NUMVERTICES);
+			Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
+			scalefactor = 0.0;
+			for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]*scalefactor_vertices[i];
+			area_base = area_base*scalefactor;
+			xDelete<IssmDouble>(scalefactor_vertices);
+		}
+
+		/*Cleanup*/
+		xDelete<IssmDouble>(weights);
+		xDelete<IssmDouble>(surfaces);
+		xDelete<IssmDouble>(bases);
+		xDelete<IssmDouble>(Hice);
+   }
 	else{
 		/*First get back the area of the base*/
@@ -3561,8 +3601,10 @@
 	}
 
+	
 	/*Cleanup & return: */
 	xDelete<int>(indices);
 	xDelete<IssmDouble>(H);
 	xDelete<IssmDouble>(SF);
+	xDelete<IssmDouble>(lsf);
 
 	if(domaintype==Domain2DverticalEnum){