Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15640)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15641)
@@ -8083,18 +8083,23 @@
 	if(!isfront) return NULL;
 
-	/*Fetch number of nodes and dof for this finite element*/
-	int         numnodes = this->NumberofNodes();
-	int         numdof   = numnodes*NDOF2;
+	/*Intermediaries*/
 	IssmDouble  rho_ice,rho_water,gravity;
 	IssmDouble  Jdet,surface,z_g,water_pressure,ice_pressure,air_pressure;
 	IssmDouble  surface_under_water,base_under_water,pressure;
-	GaussPenta*  gauss;
-	IssmDouble* basis = xNew<IssmDouble>(numnodes);
 	IssmDouble  xyz_list_front[4][3];
 	IssmDouble  area_coordinates[4][3];
 	IssmDouble  normal[3];
-
+	GaussPenta* gauss;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+	int numdof   = numnodes*NDOF2;
+
+	/*Initialize Element vector and other vectors*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,HOApproximationEnum);
+	ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters,HOApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
 	Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
 	rho_water=matpar->GetRhoWater();
@@ -8105,5 +8110,5 @@
 	GetQuadNormal(&normal[0],xyz_list_front);
 
-	/*Start looping on Gaussian points*/
+	/*Initialize gauss points*/
 	IssmDouble zmax=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
 	IssmDouble zmin=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
@@ -8111,4 +8116,5 @@
 	else                 gauss=new GaussPenta(area_coordinates,3,3);
 
+	/* Start  looping on the number of gaussian points: */
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
@@ -8119,8 +8125,8 @@
 		GetQuadJacobianDeterminant(&Jdet,xyz_list_front,gauss);
 
-		water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
-		ice_pressure=rho_ice*gravity*(surface-z_g);
-		air_pressure=0;
-		pressure = ice_pressure + water_pressure + air_pressure;
+		water_pressure = rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
+		ice_pressure   = rho_ice*gravity*(surface-z_g);
+		air_pressure   = 0;
+		pressure       = ice_pressure + water_pressure + air_pressure;
 
 		for (int i=0;i<numnodes;i++){
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15640)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15641)
@@ -3158,23 +3158,27 @@
 	if(!isfront) return NULL;
 
-	/*Fetch number of nodes and dof for this finite element*/
-	int         numnodes = this->NumberofNodes();
-	int         numdof   = numnodes*NDOF2;
+	/*Intermediaries*/
 	IssmDouble  rho_ice,rho_water,gravity;
 	IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure,air_pressure;
 	IssmDouble  surface_under_water,base_under_water,pressure;
-	GaussTria*  gauss;
-	IssmDouble* basis = xNew<IssmDouble>(numnodes);
 	IssmDouble  xyz_list_front[2][3];
 	IssmDouble  area_coordinates[2][3];
 	IssmDouble  normal[2];
 
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+	int numdof   = numnodes*NDOF2;
+
+	/*Initialize Element vector and other vectors*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,SSAApproximationEnum);
+	ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters,HOApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
 	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
 	Input* bed_input      =inputs->GetInput(BedEnum);       _assert_(bed_input);
-	rho_water=matpar->GetRhoWater();
-	rho_ice  =matpar->GetRhoIce();
-	gravity  =matpar->GetG();
+	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,2);
@@ -3182,6 +3186,5 @@
 
 	/*Start looping on Gaussian points*/
-	gauss=new GaussTria(area_coordinates,3);
-
+	GaussTria* gauss=new GaussTria(area_coordinates,3);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
@@ -3192,9 +3195,9 @@
 		GetNodalFunctions(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*(pow(surface_under_water,2) - pow(base_under_water,2));
-		ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
-		air_pressure=0;
+		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*(pow(surface_under_water,2) - pow(base_under_water,2));
+		ice_pressure   = 1.0/2.0*gravity*rho_ice*pow(thickness,2);
+		air_pressure   = 0;
 		pressure = ice_pressure + water_pressure + air_pressure;