Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15813)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15814)
@@ -817,16 +817,16 @@
 	bool               floating=true;
 	int                point;
 	const IssmPDouble  epsilon= 1.e-15;
-	IssmDouble         gl[3];
+	IssmDouble         gl[NUMVERTICES];
 	IssmDouble         f1,f2;
 
 	/*Recover parameters and values*/
 	GetInputListOnVertices(&gl[0],GLlevelsetEnum);
 
 	/*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*/
 	if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
@@ -868,19 +868,19 @@
 IssmDouble Tria::GetGroundedPortion(IssmDouble* xyz_list){
 	/*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];
 
 	/*Recover parameters and values*/
 	GetInputListOnVertices(&gl[0],GLlevelsetEnum);
 
 	/*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*/
 	if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15813)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15814)
@@ -1023,16 +1023,16 @@
 	bool               mainlyfloating = true;
 	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*/
 	GetInputListOnVertices(&gl[0],GLlevelsetEnum);
 
 	/*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*/
 	if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
@@ -1047,23 +1047,23 @@
 
 		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*/
 			s1=gl[2]/(gl[2]-gl[1]);
 			s2=gl[2]/(gl[2]-gl[0]);
 
 			/*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
 			/*Coordinates of point 0: same as initial point 2*/
@@ -1107,13 +1107,13 @@
 		}
 
 		/*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;
 }
@@ -7659,7 +7659,7 @@
 	int         i,j;
 	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;
 	IssmDouble  DL_scalar;
@@ -7692,7 +7692,7 @@
 	friction=new Friction("2d",inputs,matpar,analysis_type);
 
 	/*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: */
 	gauss=new GaussPenta(0,1,2,2);