Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15813)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15814)
@@ -1024,6 +1024,6 @@
 	const IssmPDouble  epsilon= 1.e-15;
 	IssmDouble         phi,s1,s2,area_init,area_grounded;
-	IssmDouble         gl[3];
-	IssmDouble         xyz_bis[3][3];
+	IssmDouble         gl[NUMVERTICES];
+	IssmDouble         xyz_bis[NUMVERTICES2D][3];
 
 	/*Recover parameters and values*/
@@ -1031,7 +1031,7 @@
 
 	/*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;
+	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 grounded or floating*/
@@ -1048,7 +1048,7 @@
 		if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
 			/*Coordinates of point 2: same as initial point 2*/
-			xyz_bis[2][0]=*(xyz_list+3*2+0);
-			xyz_bis[2][1]=*(xyz_list+3*2+1);
-			xyz_bis[2][2]=*(xyz_list+3*2+2);
+			xyz_bis[2][0]=xyz_list[3*2+0];
+			xyz_bis[2][1]=xyz_list[3*2+1];
+			xyz_bis[2][2]=xyz_list[3*2+2];
 
 			/*Portion of the segments*/
@@ -1057,12 +1057,12 @@
 
 			/*New point 1*/
-			xyz_bis[1][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*2+0));
-			xyz_bis[1][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*2+1));
-			xyz_bis[1][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*2+2));
+			xyz_bis[1][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
+			xyz_bis[1][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
+			xyz_bis[1][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
 
 			/*New point 0*/
-			xyz_bis[0][0]=*(xyz_list+3*2+0)+s2*(*(xyz_list+3*0+0)-*(xyz_list+3*2+0));
-			xyz_bis[0][1]=*(xyz_list+3*2+1)+s2*(*(xyz_list+3*0+1)-*(xyz_list+3*2+1));
-			xyz_bis[0][2]=*(xyz_list+3*2+2)+s2*(*(xyz_list+3*0+2)-*(xyz_list+3*2+2));
+			xyz_bis[0][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
+			xyz_bis[0][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
+			xyz_bis[0][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
 		}
 		else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
@@ -1108,11 +1108,11 @@
 
 		/*Compute fraction of grounded element*/
-		GetJacobianDeterminant(&area_init, xyz_list,NULL);
-		GetJacobianDeterminant(&area_grounded, &xyz_bis[0][0],NULL);
+		GetTriaJacobianDeterminant(&area_init, xyz_list,NULL);
+		GetTriaJacobianDeterminant(&area_grounded, &xyz_bis[0][0],NULL);
 		if(mainlyfloating==true) area_grounded=area_init-area_grounded;
 		phi=area_grounded/area_init;
 	}
 
-	if(phi>1 || phi<0) _error_("Error. Problem with portion of grounded element: value should be between 0 and 1");
+	if(phi>1. || phi<0.) _error_("Error. Problem with portion of grounded element: value should be between 0 and 1");
 
 	return phi;
@@ -7660,5 +7660,5 @@
 	int         analysis_type,migration_style;
 	IssmDouble  xyz_list[NUMVERTICES][3];
-	IssmDouble  xyz_list_tria[NUMVERTICES2D][3]={0.};
+	IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
 	IssmDouble  alpha2,Jdet;
 	IssmDouble  phi=1.0;
@@ -7693,5 +7693,5 @@
 
 	/*Recover portion of element that is grounded*/
-	if(migration_style==SubelementMigrationEnum) phi=this->GetGroundedPortion(&xyz_list[0][0]);
+	if(migration_style==SubelementMigrationEnum) phi=this->GetGroundedPortion(&xyz_list_tria[0][0]);
 
 	/* Start  looping on the number of gaussian points: */
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15813)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15814)
@@ -818,5 +818,5 @@
 	int                point;
 	const IssmPDouble  epsilon= 1.e-15;
-	IssmDouble         gl[3];
+	IssmDouble         gl[NUMVERTICES];
 	IssmDouble         f1,f2;
 
@@ -825,7 +825,7 @@
 
 	/*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;
+	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 grounded or floating*/
@@ -869,8 +869,8 @@
 	/*Computeportion of the element that is grounded*/ 
 
-	bool               mainlyfloating = true;
-	const IssmPDouble  epsilon= 1.e-15;
+	bool              mainlyfloating = true;
+	const IssmPDouble epsilon        = 1.e-15;
 	IssmDouble         phi,s1,s2,area_init,area_grounded;
-	IssmDouble         gl[3];
+	IssmDouble         gl[NUMVERTICES];
 	IssmDouble         xyz_bis[3][3];
 
@@ -879,7 +879,7 @@
 
 	/*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;
+	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 grounded or floating*/