Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 26652)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 26653)
@@ -1296,5 +1296,5 @@
 	
 	/*Compute weights*/
-	gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2);
+	gauss = this->NewGauss(point,f1,f2,1-trapezeisnegative,2);
 
 	total_weight = 0.0;
@@ -1306,8 +1306,6 @@
 	}
 
-	/*Normalize to phi*/
-	if(total_weight>0.){
-		for(int i=0;i<NUMVERTICES2D;i++) weights[i] = weights[i]*phi/total_weight;
-	}
+	/*Normalizing to phi such that weights provide coefficients for integration over subelement (for averaging:phi*weights)*/
+   if(total_weight>0.) for(int i=0;i<NUMVERTICES;i++) weights[i] = weights[i]*phi/total_weight;
 	else for(int i=0;i<NUMVERTICES2D;i++) weights[i] = 0.0;
 
@@ -2207,5 +2205,5 @@
 		IssmDouble basetot;
 		height = 0.0;
-		for(int i=0;i<NUMVERTICES2D;i++) height += weights[i]*heights[i];
+		for(int i=0;i<NUMVERTICES2D;i++) height += weights[i]/phi*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;	
@@ -2217,5 +2215,5 @@
 			/*Compute loop only over lower vertices: i<NUMVERTICES2D*/
 			scalefactor = 0.0;
-			for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]*scalefactor_vertices[i];
+			for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
 			base = base*scalefactor;
 			xDelete<IssmDouble>(scalefactor_vertices);
@@ -4389,4 +4387,5 @@
 	IssmDouble base,smb,rho_ice,scalefactor;
 	IssmDouble Total_Smb=0;
+	IssmDouble lsf[NUMVERTICES];
 	IssmDouble xyz_list[NUMVERTICES][3];
 
@@ -4404,14 +4403,45 @@
 
 	/*Now get the average SMB over the element*/
-	Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
-
-	smb_input->GetInputAverage(&smb);
-	if(scaled==true){
-		Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
-		scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
-	}
+	Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+   if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
+		/*Partially ice-covered element*/
+		bool mainlyice;
+      int point;
+      IssmDouble* smb_vertices = xNew<IssmDouble>(NUMVERTICES);
+		IssmDouble  weights[NUMVERTICES2D];
+		IssmDouble  lsf2d[NUMVERTICES2D];
+      IssmDouble f1,f2,phi;
+      Element::GetInputListOnVertices(&smb_vertices[0],SmbMassBalanceEnum);
+		for(int i=0;i<NUMVERTICES2D;i++) lsf2d[i] = lsf[i];
+		GetFractionGeometry2D(weights,&phi,&point,&f1,&f2,&mainlyice,lsf2d);
+		smb = 0.0;
+		for(int i=0;i<NUMVERTICES2D;i++) smb += weights[i]*smb_vertices[i];
+
+		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]/phi*scalefactor_vertices[i];
+         xDelete<IssmDouble>(scalefactor_vertices);
+		}
+		else scalefactor = 1.0;
+
+		/*Cleanup*/
+      xDelete<IssmDouble>(smb_vertices);
+	}
+
 	else{
-		scalefactor=1.;
-	}
+		/*Fully ice-covered element*/
+		Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
+		smb_input->GetInputAverage(&smb);
+
+		if(scaled==true){
+			Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
+			scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
+		}
+		else scalefactor=1.0;
+	}
+
 	Total_Smb=rho_ice*base*smb*scalefactor;// smb on element in kg s-1
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26652)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26653)
@@ -1817,5 +1817,5 @@
 
 	/*Compute weights:*/
-	gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2); 
+	gauss = this->NewGauss(point,f1,f2,1-trapezeisnegative,2); //VV correction (16Nov2021)
 
 	total_weight=0;
@@ -1826,5 +1826,5 @@
 		total_weight+=gauss->weight;
 	}
-	//normalize to phi. 
+	/*Normalizing to phi such that weights provide coefficients for integration over subelement (for averaging:phi*weights)*/  
 	if(total_weight>0.) for(int i=0;i<NUMVERTICES;i++)weights[i]/=total_weight/phi; 
 	else for(int i=0;i<NUMVERTICES;i++)weights[i]=0;
@@ -3546,6 +3546,6 @@
 	}*/
 
-   IssmDouble* lsf = xNew<IssmDouble>(NUMVERTICES);
-   Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+   IssmDouble lsf[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)){
@@ -3563,5 +3563,6 @@
       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];
+		/*Use weights[i]/phi to get average thickness over subelement*/
+      for(int i=0;i<NUMVERTICES;i++) Haverage += weights[i]/phi*Hice[i];
 		/*Get back area of ice-covered base*/
 		area_basetot = this->GetArea();
@@ -3573,5 +3574,5 @@
 			Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
 			scalefactor = 0.0;
-			for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]*scalefactor_vertices[i];
+			for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
 			area_base = area_base*scalefactor;
 			xDelete<IssmDouble>(scalefactor_vertices);
@@ -3606,5 +3607,4 @@
 	xDelete<IssmDouble>(H);
 	xDelete<IssmDouble>(SF);
-	xDelete<IssmDouble>(lsf);
 
 	if(domaintype==Domain2DverticalEnum){
@@ -5479,4 +5479,5 @@
 	IssmDouble base,smb,rho_ice,scalefactor;
 	IssmDouble Total_Smb=0;
+	IssmDouble lsf[NUMVERTICES];
 	IssmDouble xyz_list[NUMVERTICES][3];
 
@@ -5492,15 +5493,45 @@
 	 * 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]));	// area of element in m2
-
+	
 	/*Now get the average SMB over the element*/
-	Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
-	smb_input->GetInputAverage(&smb);	// average smb on element in m ice s-1
-	if(scaled==true){
-		Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
-		scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
+	Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+	if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
+		/*Partially ice-covered element*/
+		bool mainlyice;
+      int point;
+      IssmDouble* weights       = xNew<IssmDouble>(NUMVERTICES);
+      IssmDouble* smb_vertices  = xNew<IssmDouble>(NUMVERTICES);
+      IssmDouble f1,f2,phi;
+		
+		Element::GetInputListOnVertices(&smb_vertices[0],SmbMassBalanceEnum);
+		GetFractionGeometry(weights,&phi,&point,&f1,&f2,&mainlyice,lsf);
+		smb = 0.0;
+		for(int i=0;i<NUMVERTICES;i++) smb += weights[i]*smb_vertices[i];
+	
+		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]/phi*scalefactor_vertices[i];
+         xDelete<IssmDouble>(scalefactor_vertices);
+      }
+		else scalefactor = 1.0;
+
+		/*Cleanup*/
+      xDelete<IssmDouble>(weights);
+      xDelete<IssmDouble>(smb_vertices);
 	}
 	else{
-		scalefactor=1.;
-	}
+		/*Fully ice-covered element*/
+		Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
+		smb_input->GetInputAverage(&smb);   // average smb on element in m ice s-1
+
+		if(scaled==true){
+			Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
+			scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
+		}
+		else scalefactor=1.0;
+	}
+	
    Total_Smb=rho_ice*base*smb*scalefactor;	// smb on element in kg s-1