Index: /issm/trunk/src/c/objects/Icefront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Icefront.cpp	(revision 441)
+++ /issm/trunk/src/c/objects/Icefront.cpp	(revision 442)
@@ -686,5 +686,4 @@
 	QuadPressureLoadStokes(pe_g,matpar->GetRhoWater(),matpar->GetRhoIce(),matpar->GetG(),thickness_list_quad,bed_list_quad,normal1,normal2,normal3,normal4,&xyz_list_quad[0][0],fill);
 
-
 	#ifdef _DEBUGELEMENTS_
 	if(my_rank==RANK && eid==ELID){ 
@@ -702,5 +701,5 @@
 		printf("fill %i\n",fill);
 		printf("pe_g->terms\n");
-		for(i=0;i<numgridsload*NDOF4;i++){
+		for(i=0;i<numdofloads*NDOF4;i++){
 			printf("%g ",*(pe_g->terms+i));
 		}
@@ -1180,9 +1179,8 @@
 	double J[4];
 	double z_g[4];
-	double ice_pressure_tria[4];
 	double air_pressure_tria;
 	double water_level_above_g_tria;
 	double water_pressure_tria;
-	double pressure_tria[4];
+	double pressure_tria;
 
 	/*To use tria functionalities: */
@@ -1331,7 +1329,4 @@
 			//Loop on the grids of the quad
 			//Calculate the ice pressure
-			for(j=0;j<4;j++){
-				ice_pressure_tria[j]=gravity*rho_ice*(surface_list[j]-z_g[i]);
-			}
 			air_pressure_tria=0;
 
@@ -1348,29 +1343,23 @@
 			}
 
-			//Add pressure from triangle i
-			//Loop on the grids of the quad
-			for(j=0;j<4;j++){
-				pressure_tria[j] = ice_pressure_tria[j] + water_pressure_tria + air_pressure_tria;
-			}
-
-
-			pe_g[0]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*nx[i];
-			pe_g[1]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*ny[i];
-			pe_g[2]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*nz[i];
+			//Add pressure 
+			pressure_tria = water_pressure_tria + air_pressure_tria;
+
+			pe_g[0]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][0]*nx[i];
+			pe_g[1]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][0]*ny[i];
+			pe_g[2]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][0]*nz[i];
 			pe_g[3]+= 0;
-			pe_g[4]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*nx[i];
-			pe_g[5]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*ny[i];
-			pe_g[6]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*nz[i];
+			pe_g[4]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][1]*nx[i];
+			pe_g[5]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][1]*ny[i];
+			pe_g[6]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][1]*nz[i];
 			pe_g[7]+= 0;
-			pe_g[8]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*nx[i];
-			pe_g[9]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*ny[i];
-			pe_g[10]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*nz[i];
+			pe_g[8]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][2]*nx[i];
+			pe_g[9]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][2]*ny[i];
+			pe_g[10]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][2]*nz[i];
 			pe_g[11]+= 0;
-			pe_g[12]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*nx[i];
-			pe_g[13]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*ny[i];
-			pe_g[14]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*nz[i];
+			pe_g[12]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][3]*nx[i];
+			pe_g[13]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][3]*ny[i];
+			pe_g[14]+= J[i]*gauss_weight* pressure_tria*l1l4_tria[i][3]*nz[i];
 			pe_g[15]+= 0;
-
-			
 
 		} //for(i=0;i<4;i++)
Index: /issm/trunk/src/c/objects/Matice.cpp
===================================================================
--- /issm/trunk/src/c/objects/Matice.cpp	(revision 441)
+++ /issm/trunk/src/c/objects/Matice.cpp	(revision 442)
@@ -333,5 +333,5 @@
 	double eps0;
 
-	eps0=10^-27;
+	eps0=pow(10,-27);
 	
 	if (n==1){
@@ -362,5 +362,4 @@
 			else{
 				e=(n-1)/2/n;
-			
 				viscosity3d=2*B/(2*pow(A,e));
 			}
Index: /issm/trunk/src/c/objects/Pengrid.cpp
===================================================================
--- /issm/trunk/src/c/objects/Pengrid.cpp	(revision 441)
+++ /issm/trunk/src/c/objects/Pengrid.cpp	(revision 442)
@@ -102,5 +102,4 @@
 void  Pengrid::Demarshall(char** pmarshalled_dataset){
 
-	int i;
 	char* marshalled_dataset=NULL;
 
@@ -207,8 +206,7 @@
 	double slope[2];
 	int found=0;
-	double Ke[4][4]={{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}};
+	double Ke[4][4]={0.0};
 	
 	ParameterInputs* inputs=NULL;
-	
 
 	/*recover pointers: */
@@ -229,5 +227,4 @@
 	Ke[2][2]=kmax*pow(10,penalty_offset);
 	
-	
 	/*Add Ke to global matrix Kgg: */
 	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
@@ -237,5 +234,8 @@
 #define __FUNCT__ "Pengrid::PenaltyCreatePVector"
 void  Pengrid::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){
-	throw ErrorException(__FUNCT__," not implemented yet!");
+
+	/*No penalty applied, do nothing: */
+	return;
+
 }
 
Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 441)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 442)
@@ -679,5 +679,4 @@
 			/* Get Jacobian determinant: */
 			GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
-
 			DL_scalar=gauss_weight*Jdet;
 
@@ -746,5 +745,5 @@
 
 	/*matrices: */
-	double     Ke_temp[27][27]; //for the six nodes and the bubble 
+	double     Ke_temp[27][27]={0.0}; //for the six nodes and the bubble 
 	double     Ke_reduced[numdof][numdof]; //for the six nodes only
 	double     Ke_gaussian[27][27];
@@ -756,8 +755,8 @@
 	double     Jdet;
 	double     Jdet2d;
-	double     D[8][8]={{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 }};          
+	double     D[8][8]={0.0};
 	double     D_scalar;
 	double     tBD[numdof][8];
-	double     DLStokes[14][14];
+	double     DLStokes[14][14]={0.0};
 	double     tLDStokes[numdof2d][14];
 	
@@ -801,5 +800,4 @@
 	gravity=matpar->GetG();
 
-
 	/* Get node coordinates and dof list: */
 	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
@@ -808,11 +806,4 @@
 	/*recover extra inputs from users, at current convergence iteration: */
 	inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
-
-	/*Initialize Ke_temp to zero befor adding terms */
-	for (i=0;i<27;i++){
-		for (j=0;j<27;j++){
-			Ke_temp[i][j]=0;
-		}
-	}
 
 	/* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
@@ -821,5 +812,4 @@
 	   which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
 	   points, same deal, which yields 3 gaussian points.*/
-
 
 	area_order=5;
@@ -841,5 +831,4 @@
 			gauss_coord[3]=*(vert_gauss_coord+igvert);
 
-
 			/*Compute thickness: */
 			GetParameterValue(&thickness, &h[0],gauss_coord);
@@ -868,5 +857,4 @@
 			}
 
-
 			/*  Do the triple product tB*D*Bprime: */
 			MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0);
@@ -883,5 +871,4 @@
 
 	if((onbed==1) && (shelf==0)){
-
 
 		/*Build alpha2_list used by drag stiffness matrix*/
@@ -915,5 +902,4 @@
 		}
 		
-		
 		GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
 
@@ -974,5 +960,5 @@
 			for(i=0;i<numdof2d;i++){
 				for(j=0;j<numdof2d;j++){
-					K_terms[i][j]+=Ke_drag_gaussian[i][j];
+					Ke_temp[i][j]+=Ke_drag_gaussian[i][j];
 				}
 			}
@@ -993,5 +979,4 @@
 	
 	cleanup_and_return: 
-
 
 	return;
@@ -2408,6 +2393,6 @@
 	
 	#ifdef _DEBUG_ 
-	for (i=0;i<6;i++){
-		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i]);
+	for (i=0;i<7;i++){
+		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i],dh1dh7_basic[2][i]);
 	}
 
@@ -2780,7 +2765,7 @@
 
 	for (i=0;i<numgrids;i++){
-		*(dh1dh7_basic+numgrids*0+i)=Jinv[0][0]*dh1dh7_param[0][i]+Jinv[1][0]*dh1dh7_param[1][i]+Jinv[2][0]*dh1dh7_param[2][i];
-		*(dh1dh7_basic+numgrids*1+i)=Jinv[0][1]*dh1dh7_param[0][i]+Jinv[1][1]*dh1dh7_param[1][i]+Jinv[2][1]*dh1dh7_param[2][i];
-		*(dh1dh7_basic+numgrids*2+i)=Jinv[0][2]*dh1dh7_param[0][i]+Jinv[1][2]*dh1dh7_param[1][i]+Jinv[2][2]*dh1dh7_param[2][i];
+		*(dh1dh7_basic+numgrids*0+i)=Jinv[0][0]*dh1dh7_param[0][i]+Jinv[0][1]*dh1dh7_param[1][i]+Jinv[0][2]*dh1dh7_param[2][i];
+		*(dh1dh7_basic+numgrids*1+i)=Jinv[1][0]*dh1dh7_param[0][i]+Jinv[1][1]*dh1dh7_param[1][i]+Jinv[1][2]*dh1dh7_param[2][i];
+		*(dh1dh7_basic+numgrids*2+i)=Jinv[2][0]*dh1dh7_param[0][i]+Jinv[2][1]*dh1dh7_param[1][i]+Jinv[2][2]*dh1dh7_param[2][i];
 	}
 
@@ -2890,7 +2875,7 @@
 
 	/*matrices: */
-	double     Pe_temp[27]; //for the six nodes and the bubble 
-	double     Pe_gaussian[27]; //for the six nodes and the bubble 
-	double     Ke_temp[27][3]; //for the six nodes and the bubble 
+	double     Pe_temp[27]={0.0}; //for the six nodes and the bubble 
+	double     Pe_gaussian[27]={0.0}; //for the six nodes and the bubble 
+	double     Ke_temp[27][3]={0.0}; //for the six nodes and the bubble 
 	double     Pe_reduced[numdof]; //for the six nodes only
 	double     Ke_gaussian[27][3];
@@ -2902,5 +2887,5 @@
 	double     Jdet;
 	double     Jdet2d;
-	double     D[8][8]={{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 },{ 0,0,0,0,0,0,0,0 }};          
+	double     D[8][8]={0.0};
 	double     D_scalar;
 	double     tBD[27][8];
@@ -2930,14 +2915,4 @@
 	GetDofList(&doflist[0],&numberofdofspernode);
 
-	/*Initialize Ke_temp to zero befor adding terms */
-	for (i=0;i<27;i++){
-		Pe_temp[i]=0;
-		for (j=0;j<3;j++){
-			Ke_temp[i][j]=0;
-		}
-	}
-
-	
-	
 	/* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
 	   get tria gaussian points as well as segment gaussian points. For tria gaussian 
@@ -2945,6 +2920,4 @@
 	   which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
 	   points, same deal, which yields 3 gaussian points.*/
-
-	
 
 	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
@@ -3017,5 +2990,4 @@
 	}
 
-
 	/*Deal with 2d friction at the bedrock interface: */
 	if ( (onbed==1) && (shelf==1)){
@@ -3040,5 +3012,4 @@
 			gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
 			gauss_coord_tria[2]=*(third_gauss_area_coord+igarea);
-			
 
 			/*Get the Jacobian determinant */
@@ -3076,5 +3047,4 @@
 	}
 
-
 	/*Add P_terms to global vector pg: */
 	VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);