Index: issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16837)
+++ issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16838)
@@ -1110,4 +1110,63 @@
 ElementVector* StressbalanceAnalysis::CreatePVectorSSAFront(Element* element){/*{{{*/
 
+	/*If no front, return NULL*/
+	if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
+
+	/*Intermediaries*/
+	IssmDouble  rho_ice,rho_water,gravity;
+	IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure;
+	IssmDouble  surface_under_water,base_under_water,pressure;
+	IssmDouble *xyz_list = NULL;
+	IssmDouble *xyz_list_front = NULL;
+	IssmDouble  normal[2];
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = element->NewElementVector(SSAApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* bed_input      =element->GetInput(BedEnum);       _assert_(bed_input);
+	rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
+	rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	gravity   = element->GetMaterialParameter(ConstantsGEnum);
+	element->GetVerticesCoordinates(&xyz_list);
+	element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
+	element->NormalBase(&normal[0],xyz_list_front);
+
+	/*Start looping on Gaussian points*/
+	Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		thickness_input->GetInputValue(&thickness,gauss);
+		bed_input->GetInputValue(&bed,gauss);
+		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		surface_under_water = min(0.,thickness+bed); // 0 if the top of the glacier is above water level
+		base_under_water    = min(0.,bed);           // 0 if the bottom of the glacier is above water level
+		water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
+		ice_pressure   = 1.0/2.0*gravity*rho_ice*thickness*thickness;
+		pressure = ice_pressure + water_pressure;
+
+		for (int i=0;i<numnodes;i++){
+			pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
+			pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,XYEnum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(xyz_list_front);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
 	return NULL;
 }/*}}}*/
Index: issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16837)
+++ issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16838)
@@ -82,4 +82,5 @@
 		virtual void	GetDofListPressure(int** pdoflist,int setenum)=0;
 		virtual void   JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss)=0;
+		virtual void   JacobianDeterminantSurface(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss)=0;
 		virtual void   JacobianDeterminantBase(IssmDouble* Jdet,IssmDouble* xyz_list_base,Gauss* gauss)=0;
 		virtual void   JacobianDeterminantTop(IssmDouble* Jdet,IssmDouble* xyz_list_base,Gauss* gauss)=0;
@@ -130,4 +131,5 @@
 		virtual Gauss* NewGauss(void)=0;
 		virtual Gauss* NewGauss(int order)=0;
+      virtual Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order)=0;
 		virtual Gauss* NewGaussBase(int order)=0;
 		virtual Gauss* NewGaussTop(int order)=0;
@@ -152,4 +154,5 @@
 		virtual int    PressureInterpolation()=0;
 		virtual bool   IsZeroLevelset(int levelset_enum)=0;
+		virtual void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0;
 
 		#ifdef _HAVE_RESPONSES_
Index: issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16837)
+++ issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16838)
@@ -884,5 +884,5 @@
 /*}}}*/
 /*FUNCTION Penta::GetAreaCoordinates{{{*/
-void Penta::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble xyz_zero[4][3],IssmDouble xyz_list[6][3],int numpoints){
+void Penta::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){
 	/*Computeportion of the element that is grounded*/ 
 
@@ -891,10 +891,10 @@
 	IssmDouble  xyz_bis[3][3];
 
-	area_init=fabs(xyz_list[1][0]*xyz_list[2][1] - xyz_list[1][1]*xyz_list[2][0] + xyz_list[0][0]*xyz_list[1][1] - xyz_list[0][1]*xyz_list[1][0] + xyz_list[2][0]*xyz_list[0][1] - xyz_list[2][1]*xyz_list[0][0])/2.;
+	area_init=fabs(xyz_list[1*3+0]*xyz_list[2*3+1] - xyz_list[1*3+1]*xyz_list[2*3+0] + xyz_list[0*3+0]*xyz_list[1*3+1] - xyz_list[0*3+1]*xyz_list[1*3+0] + xyz_list[2*3+0]*xyz_list[0*3+1] - xyz_list[2*3+1]*xyz_list[0*3+0])/2.;
 
 	/*Initialize xyz_list with original xyz_list of triangle coordinates*/
 	for(j=0;j<3;j++){ 
 		for(k=0;k<3;k++){
-			xyz_bis[j][k]=xyz_list[j][k];
+			xyz_bis[j][k]=xyz_list[j*3+k];
 		}
 	}
@@ -903,5 +903,5 @@
 			for(k=0;k<3;k++){
 				/*Change appropriate line*/
-				xyz_bis[j][k]=xyz_zero[i][k];
+				xyz_bis[j][k]=xyz_zero[i*3+k];
 			}
 
@@ -913,5 +913,5 @@
 			for(k=0;k<3;k++){
 				/*Reinitialize xyz_list with original coordinates*/
-				xyz_bis[j][k]=xyz_list[j][k];
+				xyz_bis[j][k]=xyz_list[j*3+k];
 			}
 		}
@@ -1512,5 +1512,5 @@
 /*}}}*/
 /*FUNCTION Penta::GetQuadNormal {{{*/
-void Penta:: GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){
+void Penta:: GetQuadNormal(IssmDouble* normal,IssmDouble* xyz_list){
 
 	/*Build unit outward pointing vector*/
@@ -1519,10 +1519,10 @@
 	IssmDouble norm;
 
-	AB[0]=xyz_list[1][0] - xyz_list[0][0];
-	AB[1]=xyz_list[1][1] - xyz_list[0][1];
-	AB[2]=xyz_list[1][2] - xyz_list[0][2];
-	AC[0]=xyz_list[2][0] - xyz_list[0][0];
-	AC[1]=xyz_list[2][1] - xyz_list[0][1];
-	AC[2]=xyz_list[2][2] - xyz_list[0][2];
+	AB[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
+	AB[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
+	AB[2]=xyz_list[1*3+2] - xyz_list[0*3+2];
+	AC[0]=xyz_list[2*3+0] - xyz_list[0*3+0];
+	AC[1]=xyz_list[2*3+1] - xyz_list[0*3+1];
+	AC[2]=xyz_list[2*3+2] - xyz_list[0*3+2];
 
 	cross(normal,AB,AC);
@@ -1648,6 +1648,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetZeroLevelsetCoordinates{{{*/
-void Penta::GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[6][3],int levelsetenum){
+/*FUNCTION Penta::ZeroLevelsetCoordinates{{{*/
+void Penta::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){
 	/*Computeportion of the element that is grounded*/ 
 
@@ -1655,4 +1655,5 @@
 	IssmDouble  s1,s2;
 	IssmDouble  levelset[NUMVERTICES];
+	IssmDouble* xyz_zero = xNew<IssmDouble>(4*3);
 
 	/*Recover parameters and values*/
@@ -1666,22 +1667,22 @@
 		if(levelset[2]>0.) normal_orientation=1;
 		/*New point 1*/
-		xyz_zero[3*normal_orientation+0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
-		xyz_zero[3*normal_orientation+1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]);
-		xyz_zero[3*normal_orientation+2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]);
+		xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
+		xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
+		xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
 
 		/*New point 0*/
-		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]);
-		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]);
-		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]);
+		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
+		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
+		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]);
 
 		/*New point 3*/
-		xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[5][0]+s1*(xyz_list[4][0]-xyz_list[5][0]);
-		xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[5][1]+s1*(xyz_list[4][1]-xyz_list[5][1]);
-		xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[5][2]+s1*(xyz_list[4][2]-xyz_list[5][2]);
+		xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[5*3+0]+s1*(xyz_list[4*3+0]-xyz_list[5*3+0]);
+		xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[5*3+1]+s1*(xyz_list[4*3+1]-xyz_list[5*3+1]);
+		xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[5*3+2]+s1*(xyz_list[4*3+2]-xyz_list[5*3+2]);
 
 		/*New point 4*/
-		xyz_zero[3*(2+normal_orientation)+0]=xyz_list[5][0]+s2*(xyz_list[3][0]-xyz_list[5][0]);
-		xyz_zero[3*(2+normal_orientation)+1]=xyz_list[5][1]+s2*(xyz_list[3][1]-xyz_list[5][1]);
-		xyz_zero[3*(2+normal_orientation)+2]=xyz_list[5][2]+s2*(xyz_list[3][2]-xyz_list[5][2]);
+		xyz_zero[3*(2+normal_orientation)+0]=xyz_list[5*3+0]+s2*(xyz_list[3*3+0]-xyz_list[5*3+0]);
+		xyz_zero[3*(2+normal_orientation)+1]=xyz_list[5*3+1]+s2*(xyz_list[3*3+1]-xyz_list[5*3+1]);
+		xyz_zero[3*(2+normal_orientation)+2]=xyz_list[5*3+2]+s2*(xyz_list[3*3+2]-xyz_list[5*3+2]);
 	}
 	else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
@@ -1692,22 +1693,22 @@
 		if(levelset[0]>0.) normal_orientation=1;
 		/*New point 1*/
-		xyz_zero[3*normal_orientation+0]=xyz_list[0][0]+s1*(xyz_list[2][0]-xyz_list[0][0]);
-		xyz_zero[3*normal_orientation+1]=xyz_list[0][1]+s1*(xyz_list[2][1]-xyz_list[0][1]);
-		xyz_zero[3*normal_orientation+2]=xyz_list[0][2]+s1*(xyz_list[2][2]-xyz_list[0][2]);
+		xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
+		xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
+		xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
 
 		/*New point 2*/
-		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0][0]+s2*(xyz_list[1][0]-xyz_list[0][0]);
-		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0][1]+s2*(xyz_list[1][1]-xyz_list[0][1]);
-		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0][2]+s2*(xyz_list[1][2]-xyz_list[0][2]);
+		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
+		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
+		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]);
 
 		/*New point 3*/
-		xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[3][0]+s1*(xyz_list[5][0]-xyz_list[3][0]);
-		xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[3][1]+s1*(xyz_list[5][1]-xyz_list[3][1]);
-		xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[3][2]+s1*(xyz_list[5][2]-xyz_list[3][2]);
+		xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[3*3+0]+s1*(xyz_list[5*3+0]-xyz_list[3*3+0]);
+		xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[3*3+1]+s1*(xyz_list[5*3+1]-xyz_list[3*3+1]);
+		xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[3*3+2]+s1*(xyz_list[5*3+2]-xyz_list[3*3+2]);
 
 		/*New point 4*/
-		xyz_zero[3*(2+normal_orientation)+0]=xyz_list[3][0]+s2*(xyz_list[4][0]-xyz_list[3][0]);
-		xyz_zero[3*(2+normal_orientation)+1]=xyz_list[3][1]+s2*(xyz_list[4][1]-xyz_list[3][1]);
-		xyz_zero[3*(2+normal_orientation)+2]=xyz_list[3][2]+s2*(xyz_list[4][2]-xyz_list[3][2]);
+		xyz_zero[3*(2+normal_orientation)+0]=xyz_list[3*3+0]+s2*(xyz_list[4*3+0]-xyz_list[3*3+0]);
+		xyz_zero[3*(2+normal_orientation)+1]=xyz_list[3*3+1]+s2*(xyz_list[4*3+1]-xyz_list[3*3+1]);
+		xyz_zero[3*(2+normal_orientation)+2]=xyz_list[3*3+2]+s2*(xyz_list[4*3+2]-xyz_list[3*3+2]);
 	}
 	else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
@@ -1718,84 +1719,87 @@
 		if(levelset[1]>0.) normal_orientation=1;
 		/*New point 0*/
-		xyz_zero[3*normal_orientation+0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
-		xyz_zero[3*normal_orientation+1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]);
-		xyz_zero[3*normal_orientation+2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]);
+		xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
+		xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
+		xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
 
 		/*New point 2*/
-		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]);
-		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]);
-		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]);
+		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
+		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
+		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
 
 		/*New point 3*/
-		xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[4][0]+s1*(xyz_list[3][0]-xyz_list[4][0]);
-		xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[4][1]+s1*(xyz_list[3][1]-xyz_list[4][1]);
-		xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[4][2]+s1*(xyz_list[3][2]-xyz_list[4][2]);
+		xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[4*3+0]+s1*(xyz_list[3*3+0]-xyz_list[4*3+0]);
+		xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[4*3+1]+s1*(xyz_list[3*3+1]-xyz_list[4*3+1]);
+		xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[4*3+2]+s1*(xyz_list[3*3+2]-xyz_list[4*3+2]);
 
 		/*New point 4*/
-		xyz_zero[3*(2+normal_orientation)+0]=xyz_list[4][0]+s2*(xyz_list[5][0]-xyz_list[4][0]);
-		xyz_zero[3*(2+normal_orientation)+1]=xyz_list[4][1]+s2*(xyz_list[5][1]-xyz_list[4][1]);
-		xyz_zero[3*(2+normal_orientation)+2]=xyz_list[4][2]+s2*(xyz_list[5][2]-xyz_list[4][2]);
+		xyz_zero[3*(2+normal_orientation)+0]=xyz_list[4*3+0]+s2*(xyz_list[5*3+0]-xyz_list[4*3+0]);
+		xyz_zero[3*(2+normal_orientation)+1]=xyz_list[4*3+1]+s2*(xyz_list[5*3+1]-xyz_list[4*3+1]);
+		xyz_zero[3*(2+normal_orientation)+2]=xyz_list[4*3+2]+s2*(xyz_list[5*3+2]-xyz_list[4*3+2]);
 	}
 	else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
-		xyz_zero[3*0+0]=xyz_list[0][0];
-		xyz_zero[3*0+1]=xyz_list[0][1];
-		xyz_zero[3*0+2]=xyz_list[0][2];
+		xyz_zero[3*0+0]=xyz_list[0*3+0];
+		xyz_zero[3*0+1]=xyz_list[0*3+1];
+		xyz_zero[3*0+2]=xyz_list[0*3+2];
 
 		/*New point 2*/
-		xyz_zero[3*1+0]=xyz_list[1][0];
-		xyz_zero[3*1+1]=xyz_list[1][1];
-		xyz_zero[3*1+2]=xyz_list[1][2];
+		xyz_zero[3*1+0]=xyz_list[1*3+0];
+		xyz_zero[3*1+1]=xyz_list[1*3+1];
+		xyz_zero[3*1+2]=xyz_list[1*3+2];
 
 		/*New point 3*/
-		xyz_zero[3*2+0]=xyz_list[4][0];
-		xyz_zero[3*2+1]=xyz_list[4][1];
-		xyz_zero[3*2+2]=xyz_list[4][2];
+		xyz_zero[3*2+0]=xyz_list[4*3+0];
+		xyz_zero[3*2+1]=xyz_list[4*3+1];
+		xyz_zero[3*2+2]=xyz_list[4*3+2];
 
 		/*New point 4*/
-		xyz_zero[3*3+0]=xyz_list[3][0];
-		xyz_zero[3*3+1]=xyz_list[3][1];
-		xyz_zero[3*3+2]=xyz_list[3][2];
+		xyz_zero[3*3+0]=xyz_list[3*3+0];
+		xyz_zero[3*3+1]=xyz_list[3*3+1];
+		xyz_zero[3*3+2]=xyz_list[3*3+2];
 	}
 	else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
-		xyz_zero[3*0+0]=xyz_list[2][0];
-		xyz_zero[3*0+1]=xyz_list[2][1];
-		xyz_zero[3*0+2]=xyz_list[2][2];
+		xyz_zero[3*0+0]=xyz_list[2*3+0];
+		xyz_zero[3*0+1]=xyz_list[2*3+1];
+		xyz_zero[3*0+2]=xyz_list[2*3+2];
 
 		/*New point 2*/
-		xyz_zero[3*1+0]=xyz_list[0][0];
-		xyz_zero[3*1+1]=xyz_list[0][1];
-		xyz_zero[3*1+2]=xyz_list[0][2];
+		xyz_zero[3*1+0]=xyz_list[0*3+0];
+		xyz_zero[3*1+1]=xyz_list[0*3+1];
+		xyz_zero[3*1+2]=xyz_list[0*3+2];
 
 		/*New point 3*/
-		xyz_zero[3*2+0]=xyz_list[3][0];
-		xyz_zero[3*2+1]=xyz_list[3][1];
-		xyz_zero[3*2+2]=xyz_list[3][2];
+		xyz_zero[3*2+0]=xyz_list[3*3+0];
+		xyz_zero[3*2+1]=xyz_list[3*3+1];
+		xyz_zero[3*2+2]=xyz_list[3*3+2];
 
 		/*New point 4*/
-		xyz_zero[3*3+0]=xyz_list[5][0];
-		xyz_zero[3*3+1]=xyz_list[5][1];
-		xyz_zero[3*3+2]=xyz_list[5][2];
+		xyz_zero[3*3+0]=xyz_list[5*3+0];
+		xyz_zero[3*3+1]=xyz_list[5*3+1];
+		xyz_zero[3*3+2]=xyz_list[5*3+2];
 	}
 	else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
-		xyz_zero[3*0+0]=xyz_list[1][0];
-		xyz_zero[3*0+1]=xyz_list[1][1];
-		xyz_zero[3*0+2]=xyz_list[1][2];
+		xyz_zero[3*0+0]=xyz_list[1*3+0];
+		xyz_zero[3*0+1]=xyz_list[1*3+1];
+		xyz_zero[3*0+2]=xyz_list[1*3+2];
 
 		/*New point 2*/
-		xyz_zero[3*1+0]=xyz_list[2][0];
-		xyz_zero[3*1+1]=xyz_list[2][1];
-		xyz_zero[3*1+2]=xyz_list[2][2];
+		xyz_zero[3*1+0]=xyz_list[2*3+0];
+		xyz_zero[3*1+1]=xyz_list[2*3+1];
+		xyz_zero[3*1+2]=xyz_list[2*3+2];
 
 		/*New point 3*/
-		xyz_zero[3*2+0]=xyz_list[5][0];
-		xyz_zero[3*2+1]=xyz_list[5][1];
-		xyz_zero[3*2+2]=xyz_list[5][2];
+		xyz_zero[3*2+0]=xyz_list[5*3+0];
+		xyz_zero[3*2+1]=xyz_list[5*3+1];
+		xyz_zero[3*2+2]=xyz_list[5*3+2];
 
 		/*New point 4*/
-		xyz_zero[3*3+0]=xyz_list[4][0];
-		xyz_zero[3*3+1]=xyz_list[4][1];
-		xyz_zero[3*3+2]=xyz_list[4][2];
+		xyz_zero[3*3+0]=xyz_list[4*3+0];
+		xyz_zero[3*3+1]=xyz_list[4*3+1];
+		xyz_zero[3*3+2]=xyz_list[4*3+2];
 	}
 	else _error_("Case not covered");
+
+	/*Assign output pointer*/
+	*pxyz_zero= xyz_zero;
 }
 /*}}}*/
@@ -8764,5 +8768,5 @@
 	IssmDouble  surface_under_water,base_under_water,pressure;
 	IssmDouble  xyz_list[NUMVERTICES][3];
-	IssmDouble  xyz_list_front[4][3];
+	IssmDouble* xyz_list_front = NULL;
 	IssmDouble  area_coordinates[4][3];
 	IssmDouble  normal[3];
@@ -8782,6 +8786,6 @@
 	rho_ice  =matpar->GetRhoIce();
 	gravity  =matpar->GetG();
-	GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum);
-	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
+	ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
+	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],4);
 	GetQuadNormal(&normal[0],xyz_list_front);
 
@@ -8854,5 +8858,5 @@
 	IssmDouble  surface_under_water,base_under_water,pressure;
 	IssmDouble  xyz_list[NUMVERTICES][3];
-	IssmDouble  xyz_list_front[4][3];
+	IssmDouble* xyz_list_front = NULL;
 	IssmDouble  area_coordinates[4][3];
 	IssmDouble  normal[3];
@@ -8869,14 +8873,14 @@
 
 	/*Initialize Element matrix and vectors*/
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
 	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
 
 	/*Retrieve all inputs and parameters*/
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	rho_water=matpar->GetRhoWater();
 	rho_ice  =matpar->GetRhoIce();
 	gravity  =matpar->GetG();
-	GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum);
-	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
+	ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
+	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],4);
 	GetQuadNormal(&normal[0],xyz_list_front);
 
@@ -8910,4 +8914,5 @@
 	/*Clean up and return*/
 	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(xyz_list_front);
 	xDelete<IssmDouble>(vbasis);
 	delete gauss;
Index: issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16837)
+++ issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16838)
@@ -114,4 +114,5 @@
 		int    PressureInterpolation();
 		bool   IsZeroLevelset(int levelset_enum);
+		void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
 
 		void   ResultInterpolation(int* pinterpolation,int output_enum);
@@ -211,5 +212,5 @@
 		IssmDouble     EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
 		IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
-		void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[6][3],int numpoints);
+		void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
 
 		void	         GetVertexPidList(int* doflist);
@@ -231,9 +232,8 @@
 		void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		void	         GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
-		void           GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
+		void           GetQuadNormal(IssmDouble* normal,IssmDouble* xyz_list);
 		void           GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input);
 		void           GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
 		Penta*         GetUpperElement(void);
-		void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[6][3],int levelsetenum);
 		Penta*         GetLowerElement(void);
 		void           InputChangeName(int input_enum, int enum_type_old);
@@ -247,4 +247,5 @@
 		bool           IsNodeOnShelfFromFlags(IssmDouble* flags);
 		void           JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
+		void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
 		void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
 		void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
@@ -252,4 +253,5 @@
 		Gauss*         NewGauss(void);
 		Gauss*         NewGauss(int order);
+      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
 		Gauss*         NewGaussBase(int order);
 		Gauss*         NewGaussTop(int order);
Index: issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 16837)
+++ issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 16838)
@@ -1945,21 +1945,21 @@
 /*}}}*/
 /*FUNCTION PentaRef::GetQuadJacobianDeterminant{{{*/
-void PentaRef::GetQuadJacobianDeterminant(IssmDouble* Jdet,IssmDouble xyz_list[4][3],Gauss* gauss){
+void PentaRef::GetQuadJacobianDeterminant(IssmDouble* Jdet,IssmDouble* xyz_list,Gauss* gauss){
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
 
 	IssmDouble x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4;
 
-	x1=xyz_list[0][0];
-	y1=xyz_list[0][1];
-	z1=xyz_list[0][2];
-	x2=xyz_list[1][0];
-	y2=xyz_list[1][1];
-	z2=xyz_list[1][2];
-	x3=xyz_list[2][0];
-	y3=xyz_list[2][1];
-	z3=xyz_list[2][2];
-	x4=xyz_list[3][0];
-	y4=xyz_list[3][1];
-	z4=xyz_list[3][2];
+	x1=xyz_list[0*3+0];
+	y1=xyz_list[0*3+1];
+	z1=xyz_list[0*3+2];
+	x2=xyz_list[1*3+0];
+	y2=xyz_list[1*3+1];
+	z2=xyz_list[1*3+2];
+	x3=xyz_list[2*3+0];
+	y3=xyz_list[2*3+1];
+	z3=xyz_list[2*3+2];
+	x4=xyz_list[3*3+0];
+	y4=xyz_list[3*3+1];
+	z4=xyz_list[3*3+2];
 
 	/*Jdet = (Area of the trapezoid)/(Area trapezoid ref) with AreaRef = 4*/
Index: issm/trunk-jpl/src/c/classes/Elements/PentaRef.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 16837)
+++ issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 16838)
@@ -35,5 +35,5 @@
 		void GetNodalFunctionsP1DerivativesReference(IssmDouble* dl1dl6,Gauss* gauss);
 		void GetNodalFunctionsMINIDerivativesReference(IssmDouble* dl1dl7,Gauss* gauss);
-		void GetQuadJacobianDeterminant(IssmDouble*  Jdet, IssmDouble xyz_list[4][3],Gauss* gauss);
+		void GetQuadJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
 		void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss);
 		void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
Index: issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16837)
+++ issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16838)
@@ -114,4 +114,5 @@
 		bool        IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
 		void        JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
+		void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
 		void        JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
 		void        JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
@@ -142,4 +143,5 @@
 		Gauss*      NewGauss(void){_error_("not implemented yet");};
 		Gauss*      NewGauss(int order){_error_("not implemented yet");};
+      Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
 		Gauss*      NewGaussBase(int order){_error_("not implemented yet");};
 		Gauss*      NewGaussTop(int order){_error_("not implemented yet");};
@@ -151,4 +153,5 @@
 		void        ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
 		bool        IsZeroLevelset(int levelset_enum){_error_("not implemented");};
+		void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");};
 
 		#ifdef _HAVE_THERMAL_
Index: issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16837)
+++ issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16838)
@@ -879,5 +879,5 @@
 /*}}}*/
 /*FUNCTION Tria::GetAreaCoordinates{{{*/
-void Tria::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints){
+void Tria::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){
 	/*Computeportion of the element that is grounded*/ 
 
@@ -891,5 +891,5 @@
 	for(j=0;j<3;j++){ 
 		for(k=0;k<3;k++){
-			xyz_bis[j][k]=xyz_list[j][k];
+			xyz_bis[j][k]=xyz_list[j*3+k];
 		}
 	}
@@ -898,5 +898,5 @@
 			for(k=0;k<3;k++){
 				/*Change appropriate line*/
-				xyz_bis[j][k]=xyz_zero[i][k];
+				xyz_bis[j][k]=xyz_zero[i*3+k];
 			}
 
@@ -908,5 +908,5 @@
 			for(k=0;k<3;k++){
 				/*Reinitialize xyz_list with original coordinates*/
-				xyz_bis[j][k]=xyz_list[j][k];
+				xyz_bis[j][k]=xyz_list[j*3+k];
 			}
 		}
@@ -1184,7 +1184,6 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetZeroLevelsetCoordinates{{{*/
-void Tria::GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum){
-	/*Computeportion of the element that is grounded*/ 
+/*FUNCTION Tria::ZeroLevelsetCoordinates{{{*/
+void Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){
 
 	int         normal_orientation=0;
@@ -1193,4 +1192,5 @@
 
 	/*Recover parameters and values*/
+	IssmDouble* xyz_zero = xNew<IssmDouble>(2*3);
 	GetInputListOnVertices(&levelset[0],levelsetenum);
 
@@ -1202,12 +1202,12 @@
 		if(levelset[2]>0.) normal_orientation=1;
 		/*New point 1*/
-		xyz_zero[3*normal_orientation+0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
-		xyz_zero[3*normal_orientation+1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]);
-		xyz_zero[3*normal_orientation+2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]);
+		xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
+		xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
+		xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
 
 		/*New point 0*/
-		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]);
-		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]);
-		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]);
+		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
+		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
+		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]);
 	}
 	else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
@@ -1218,12 +1218,12 @@
 		if(levelset[0]>0.) normal_orientation=1;
 		/*New point 1*/
-		xyz_zero[3*normal_orientation+0]=xyz_list[0][0]+s1*(xyz_list[2][0]-xyz_list[0][0]);
-		xyz_zero[3*normal_orientation+1]=xyz_list[0][1]+s1*(xyz_list[2][1]-xyz_list[0][1]);
-		xyz_zero[3*normal_orientation+2]=xyz_list[0][2]+s1*(xyz_list[2][2]-xyz_list[0][2]);
+		xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
+		xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
+		xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
 
 		/*New point 2*/
-		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0][0]+s2*(xyz_list[1][0]-xyz_list[0][0]);
-		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0][1]+s2*(xyz_list[1][1]-xyz_list[0][1]);
-		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0][2]+s2*(xyz_list[1][2]-xyz_list[0][2]);
+		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
+		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
+		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]);
 	}
 	else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
@@ -1234,44 +1234,47 @@
 		if(levelset[1]>0.) normal_orientation=1;
 		/*New point 0*/
-		xyz_zero[3*normal_orientation+0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
-		xyz_zero[3*normal_orientation+1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]);
-		xyz_zero[3*normal_orientation+2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]);
+		xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
+		xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
+		xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
 
 		/*New point 2*/
-		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]);
-		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]);
-		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]);
+		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
+		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
+		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
 	}
 	else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
-		xyz_zero[3*0+0]=xyz_list[0][0];
-		xyz_zero[3*0+1]=xyz_list[0][1];
-		xyz_zero[3*0+2]=xyz_list[0][2];
+		xyz_zero[3*0+0]=xyz_list[0*3+0];
+		xyz_zero[3*0+1]=xyz_list[0*3+1];
+		xyz_zero[3*0+2]=xyz_list[0*3+2];
 
 		/*New point 2*/
-		xyz_zero[3*1+0]=xyz_list[1][0];
-		xyz_zero[3*1+1]=xyz_list[1][1];
-		xyz_zero[3*1+2]=xyz_list[1][2];
+		xyz_zero[3*1+0]=xyz_list[1*3+0];
+		xyz_zero[3*1+1]=xyz_list[1*3+1];
+		xyz_zero[3*1+2]=xyz_list[1*3+2];
 	}
 	else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
-		xyz_zero[3*0+0]=xyz_list[2][0];
-		xyz_zero[3*0+1]=xyz_list[2][1];
-		xyz_zero[3*0+2]=xyz_list[2][2];
+		xyz_zero[3*0+0]=xyz_list[2*3+0];
+		xyz_zero[3*0+1]=xyz_list[2*3+1];
+		xyz_zero[3*0+2]=xyz_list[2*3+2];
 
 		/*New point 2*/
-		xyz_zero[3*1+0]=xyz_list[0][0];
-		xyz_zero[3*1+1]=xyz_list[0][1];
-		xyz_zero[3*1+2]=xyz_list[0][2];
+		xyz_zero[3*1+0]=xyz_list[0*3+0];
+		xyz_zero[3*1+1]=xyz_list[0*3+1];
+		xyz_zero[3*1+2]=xyz_list[0*3+2];
 	}
 	else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
-		xyz_zero[3*0+0]=xyz_list[1][0];
-		xyz_zero[3*0+1]=xyz_list[1][1];
-		xyz_zero[3*0+2]=xyz_list[1][2];
+		xyz_zero[3*0+0]=xyz_list[1*3+0];
+		xyz_zero[3*0+1]=xyz_list[1*3+1];
+		xyz_zero[3*0+2]=xyz_list[1*3+2];
 
 		/*New point 2*/
-		xyz_zero[3*1+0]=xyz_list[2][0];
-		xyz_zero[3*1+1]=xyz_list[2][1];
-		xyz_zero[3*1+2]=xyz_list[2][2];
+		xyz_zero[3*1+0]=xyz_list[2*3+0];
+		xyz_zero[3*1+1]=xyz_list[2*3+1];
+		xyz_zero[3*1+2]=xyz_list[2*3+2];
 	}
 	else _error_("Case not covered");
+
+	/*Assign output pointer*/
+	*pxyz_zero= xyz_zero;
 }
 /*}}}*/
@@ -2014,4 +2017,12 @@
 }
 /*}}}*/
+/*FUNCTION Tria::JacobianDeterminantSurface{{{*/
+void Tria::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){
+
+	_assert_(gauss->Enum()==GaussTriaEnum);
+	this->GetSegmentJacobianDeterminant(pJdet,xyz_list,(GaussTria*)gauss);
+
+}
+/*}}}*/
 /*FUNCTION Tria::HasEdgeOnBed {{{*/
 bool Tria::HasEdgeOnBed(){
@@ -2207,4 +2218,14 @@
 Gauss* Tria::NewGauss(int order){
 	return new GaussTria(order);
+}
+/*}}}*/
+/*FUNCTION Tria::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){{{*/
+Gauss* Tria::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){
+
+	IssmDouble  area_coordinates[2][3];
+
+	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
+
+	return new GaussTria(area_coordinates,order);
 }
 /*}}}*/
@@ -3930,5 +3951,5 @@
 	IssmDouble  Jdet,water_pressure,air_pressure,pressure;
 	IssmDouble  xyz_list[NUMVERTICES][3];
-	IssmDouble  xyz_list_front[2][3];
+	IssmDouble* xyz_list_front = NULL;
 	IssmDouble  area_coordinates[2][3];
 	IssmDouble  normal[2];
@@ -3953,11 +3974,11 @@
 	rho_ice = matpar->GetRhoIce();
 	gravity   = matpar->GetG();
-	GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum);
-	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
-	NormalBase(&normal[0],&xyz_list_front[0][0]);
+	ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
+	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2);
+	NormalBase(&normal[0],xyz_list_front);
 
 	/*Start looping on Gaussian points*/
-	IssmDouble ymax=max(xyz_list_front[0][1],xyz_list_front[1][1]);
-	IssmDouble ymin=min(xyz_list_front[0][1],xyz_list_front[1][1]);
+	IssmDouble ymax=max(xyz_list_front[0*3+1],xyz_list_front[1*3+1]);
+	IssmDouble ymin=min(xyz_list_front[0*3+1],xyz_list_front[1*3+1]);
 	if(ymax>0. && ymin<0.) gauss=new GaussTria(area_coordinates,30); //refined in vertical because of the sea level discontinuity
 	else                   gauss=new GaussTria(area_coordinates,3);
@@ -3968,5 +3989,5 @@
 		y_g=GetYcoord(gauss);
 		GetNodalFunctionsVelocity(vbasis,gauss);
-		GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss);
+		GetSegmentJacobianDeterminant(&Jdet,xyz_list_front,gauss);
 
 		water_pressure = rho_water*gravity*min(0.,y_g); //0 if the gaussian point is above water level
@@ -3980,5 +4001,4 @@
 	}
 
-
 	/*Transform coordinate system*/
 	::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
@@ -3986,4 +4006,5 @@
 	/*Clean up and return*/
 	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(xyz_list_front);
 	xDelete<IssmDouble>(vbasis);
 	delete gauss;
@@ -4201,5 +4222,5 @@
 	IssmDouble  surface_under_water,base_under_water,pressure;
 	IssmDouble  xyz_list[NUMVERTICES][3];
-	IssmDouble  xyz_list_front[2][3];
+	IssmDouble  *xyz_list_front = NULL;
 	IssmDouble  area_coordinates[2][3];
 	IssmDouble  normal[2];
@@ -4219,7 +4240,7 @@
 	rho_ice   = matpar->GetRhoIce();
 	gravity   = matpar->GetG();
-	GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum);
-	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
-	NormalBase(&normal[0],&xyz_list_front[0][0]);
+	ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
+	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2);
+	NormalBase(&normal[0],xyz_list_front);
 
 	/*Start looping on Gaussian points*/
@@ -4230,5 +4251,5 @@
 		thickness_input->GetInputValue(&thickness,gauss);
 		bed_input->GetInputValue(&bed,gauss);
-		GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss);
+		GetSegmentJacobianDeterminant(&Jdet,xyz_list_front,gauss);
 		GetNodalFunctions(basis,gauss);
 
@@ -4251,4 +4272,5 @@
 	/*Clean up and return*/
 	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list_front);
 	delete gauss;
 	return pe;
Index: issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16837)
+++ issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16838)
@@ -135,4 +135,5 @@
 		void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
 		IssmDouble  TimeAdapt();
+		void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
 		bool        IsZeroLevelset(int levelset_enum);
 
@@ -248,5 +249,5 @@
 		IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
 		IssmDouble     GetArea(void);
-		void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints);
+		void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
 		int            GetElementType(void);
 
@@ -259,5 +260,4 @@
 		void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
 		IssmDouble     GetMaterialParameter(int enum_in);
-		void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum);
 		Input*         GetInput(int inputenum);
 		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
@@ -277,4 +277,5 @@
 		bool	         IsInput(int name);
 		void           JacobianDeterminant(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
+		void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
 		void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
 		void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
@@ -282,4 +283,5 @@
 		Gauss*         NewGauss(void);
 		Gauss*         NewGauss(int order);
+      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order);
 		Gauss*         NewGaussBase(int order){_error_("not implemented yet");};
 		Gauss*         NewGaussTop(int order){_error_("not implemented yet");};
