Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 26652)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 26653)
@@ -1295,7 +1295,7 @@
 	else                  phi = f1*f2;
 	
 	/*Compute weights*/
-	gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2);
+	gauss = this->NewGauss(point,f1,f2,1-trapezeisnegative,2);
 
 	total_weight = 0.0;
 	for(int i=0;i<NUMVERTICES2D;i++)weights[i] = 0;
@@ -1305,10 +1305,8 @@
 		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;
-	}
+	/*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;
 
 	/*Cleanup*/
@@ -2206,7 +2204,7 @@
 		
 		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;	
 	
@@ -2216,7 +2214,7 @@
 			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];
+			for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
 			base = base*scalefactor;
 			xDelete<IssmDouble>(scalefactor_vertices);
 		}
@@ -4388,6 +4386,7 @@
 	/*The smb[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
 	IssmDouble base,smb,rho_ice,scalefactor;
 	IssmDouble Total_Smb=0;
+	IssmDouble lsf[NUMVERTICES];
 	IssmDouble xyz_list[NUMVERTICES][3];
 
 	/*Get material parameters :*/
@@ -4403,16 +4402,47 @@
 	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]));
 
 	/*Now get the average SMB over the element*/
-	Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
+	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];
 
-	smb_input->GetInputAverage(&smb);
-	if(scaled==true){
-		Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
-		scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
+		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
 
 	/*Return: */
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26652)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26653)
@@ -1816,7 +1816,7 @@
 	else phi=f1*f2;
 
 	/*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;
 	for(int i=0;i<NUMVERTICES;i++)weights[i]=0;
@@ -1825,7 +1825,7 @@
 		for(int i=0;i<NUMVERTICES;i++)weights[i]+=loadweights_g[i]*gauss->weight;
 		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;
 
@@ -3545,8 +3545,8 @@
 		Haverage/=IssmDouble(numthk);
 	}*/
 
-   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)){
 		bool istrapneg;
@@ -3562,7 +3562,8 @@
       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];
+		/*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();
 		area_base    = phi*area_basetot;
@@ -3572,7 +3573,7 @@
 			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];
+			for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
 			area_base = area_base*scalefactor;
 			xDelete<IssmDouble>(scalefactor_vertices);
 		}
@@ -3605,7 +3606,6 @@
 	xDelete<int>(indices);
 	xDelete<IssmDouble>(H);
 	xDelete<IssmDouble>(SF);
-	xDelete<IssmDouble>(lsf);
 
 	if(domaintype==Domain2DverticalEnum){
 	  return area_base;
@@ -5478,6 +5478,7 @@
 	/*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/
 	IssmDouble base,smb,rho_ice,scalefactor;
 	IssmDouble Total_Smb=0;
+	IssmDouble lsf[NUMVERTICES];
 	IssmDouble xyz_list[NUMVERTICES][3];
 
 	/*Get material parameters :*/
@@ -5491,17 +5492,47 @@
 	 * http://en.wikipedia.org/wiki/Triangle
 	 * 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*/
+	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;
 
-	/*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
+		/*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
 
 	/*Return: */