Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 510)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 511)
@@ -4,5 +4,5 @@
 
 #ifdef HAVE_CONFIG_H
-	#include "config.h"
+#include "config.h"
 #else
 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
@@ -20,8 +20,8 @@
 }
 Penta::Penta( int penta_id, int penta_mid, int penta_mparid, int penta_node_ids[6], double penta_h[6], double penta_s[6], double penta_b[6], double penta_k[6], int penta_friction_type, 
-				double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, double penta_meanvel,double penta_epsvel, 
-				int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6], 
-				int penta_artdiff, int penta_thermal_steadystate,double penta_viscosity_overshoot,double penta_stokesreconditioning){
-	
+			double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, double penta_meanvel,double penta_epsvel, 
+			int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6], 
+			int penta_artdiff, int penta_thermal_steadystate,double penta_viscosity_overshoot,double penta_stokesreconditioning){
+
 	int i;
 
@@ -83,5 +83,5 @@
 	printf("   b=[%i,%i,%i,%i,%i,%i]\n",b[0],b[1],b[2],b[3],b[4],b[5]);
 	printf("   k=[%i,%i,%i,%i,%i,%i]\n",k[0],k[1],k[2],k[3],k[4],k[5]);
-	
+
 	printf("   friction_type: %i\n",friction_type);
 	printf("   p: %g\n",p);
@@ -93,5 +93,5 @@
 	printf("   epsvel: %g\n",epsvel);
 	printf("   collapse: %i\n",collapse);
-	
+
 	printf("   melting=[%i,%i,%i,%i,%i,%i]\n",melting[0],melting[1],melting[2],melting[3],melting[4],melting[5]);
 	printf("   accumulation=[%i,%i,%i,%i,%i,%i]\n",accumulation[0],accumulation[1],accumulation[2],accumulation[3],accumulation[4],accumulation[5]);
@@ -114,8 +114,8 @@
 	/*get enum type of Penta: */
 	enum_type=PentaEnum();
-	
+
 	/*marshall enum: */
 	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
-	
+
 	/*marshall Penta data: */
 	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
@@ -149,42 +149,42 @@
 	memcpy(marshalled_dataset,&viscosity_overshoot,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
 	memcpy(marshalled_dataset,&stokesreconditioning,sizeof(stokesreconditioning));marshalled_dataset+=sizeof(stokesreconditioning);
-	
+
 	*pmarshalled_dataset=marshalled_dataset;
 	return;
 }
-		
+
 int   Penta::MarshallSize(){
 
 	return sizeof(id)+
-		sizeof(mid)+
-		sizeof(mparid)+
-		sizeof(node_ids)+
-		sizeof(nodes)+
-		sizeof(node_offsets)+
-		sizeof(matice)+
-		sizeof(matice_offset)+
-		sizeof(matpar)+
-		sizeof(matpar_offset)+
-		sizeof(h)+
-		sizeof(s)+
-		sizeof(b)+
-		sizeof(k)+
-		sizeof(friction_type)+
-		sizeof(p)+
-		sizeof(q)+
-		sizeof(shelf)+
-		sizeof(onbed)+
-		sizeof(onsurface)+
-		sizeof(meanvel)+
-		sizeof(epsvel)+
-		sizeof(collapse)+
-		sizeof(melting)+
-		sizeof(accumulation)+
-		sizeof(geothermalflux)+
-		sizeof(artdiff)+
-		sizeof(thermal_steadystate) +
-		sizeof(viscosity_overshoot) +
-		sizeof(stokesreconditioning) +
-		sizeof(int); //sizeof(int) for enum type
+	  sizeof(mid)+
+	  sizeof(mparid)+
+	  sizeof(node_ids)+
+	  sizeof(nodes)+
+	  sizeof(node_offsets)+
+	  sizeof(matice)+
+	  sizeof(matice_offset)+
+	  sizeof(matpar)+
+	  sizeof(matpar_offset)+
+	  sizeof(h)+
+	  sizeof(s)+
+	  sizeof(b)+
+	  sizeof(k)+
+	  sizeof(friction_type)+
+	  sizeof(p)+
+	  sizeof(q)+
+	  sizeof(shelf)+
+	  sizeof(onbed)+
+	  sizeof(onsurface)+
+	  sizeof(meanvel)+
+	  sizeof(epsvel)+
+	  sizeof(collapse)+
+	  sizeof(melting)+
+	  sizeof(accumulation)+
+	  sizeof(geothermalflux)+
+	  sizeof(artdiff)+
+	  sizeof(thermal_steadystate) +
+	  sizeof(viscosity_overshoot) +
+	  sizeof(stokesreconditioning) +
+	  sizeof(int); //sizeof(int) for enum type
 }
 
@@ -262,5 +262,5 @@
 
 	int i;
-	
+
 	DataSet* loadsin=NULL;
 	DataSet* nodesin=NULL;
@@ -274,5 +274,5 @@
 	/*Link this element with its nodes, ie find pointers to the nodes in the nodes dataset.: */
 	ResolvePointers((Object**)nodes,node_ids,node_offsets,6,nodesin);
-	
+
 	/*Same for materials: */
 	ResolvePointers((Object**)&matice,&mid,&matice_offset,1,materialsin);
@@ -293,30 +293,30 @@
 	}
 	else if (analysis_type==DiagnosticAnalysisEnum()){
-	
+
 		if (sub_analysis_type==HorizAnalysisEnum()){
-		
+
 			CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
 		}
 		else if (sub_analysis_type==VertAnalysisEnum()){
-		
+
 			CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type,sub_analysis_type);
 		}
 		else if (sub_analysis_type==StokesAnalysisEnum()){
-		
+
 			CreateKMatrixDiagnosticStokes( Kgg,inputs,analysis_type,sub_analysis_type);
-		
+
 		}
 		else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
 	}
 	else if (analysis_type==SlopeComputeAnalysisEnum()){
-		
+
 		CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type);
 	}
 	else if (analysis_type==ThermalAnalysisEnum()){
-		
+
 		CreateKMatrixThermal( Kgg,inputs,analysis_type,sub_analysis_type);
 	}
 	else if (analysis_type==MeltingAnalysisEnum()){
-		
+
 		CreateKMatrixMelting( Kgg,inputs,analysis_type,sub_analysis_type);
 	}
@@ -342,6 +342,6 @@
 	int          doflist[numdof];
 	int          numberofdofspernode;
-	
-	
+
+
 	/* 3d gaussian points: */
 	int     num_gauss,ig;
@@ -359,5 +359,5 @@
 	int     num_area_gauss;
 	double  gauss_weight;
-	
+
 	/* 2d gaussian point: */
 	int     num_gauss2d;
@@ -391,5 +391,5 @@
 	double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
 	double Jdet;
-	
+
 	/*slope: */
 	double  slope[2]={0.0,0.0};
@@ -454,5 +454,5 @@
 		}
 
-		#ifdef _DEBUGELEMENTS_
+#ifdef _DEBUGELEMENTS_
 		if(my_rank==RANK && id==ELID){ 
 			printf("El id %i Rank %i PentaElement input list before gaussian loop: \n",ELID,RANK); 
@@ -471,5 +471,5 @@
 			printf("   temperature_average [%g %g %g %g %g %g]\n",temperature_average_list[0],temperature_average_list[1],temperature_average_list[2],temperature_average_list[3],temperature_average_list[4],temperature_average_list[5]);
 		}
-		#endif
+#endif
 
 		/*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
@@ -484,5 +484,5 @@
 		GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
 
-		#ifdef _DEBUGGAUSS_
+#ifdef _DEBUGGAUSS_
 		if(my_rank==RANK && id==ELID){ 
 			printf("El id %i Rank %i PentaElement gauss points\n",ELID,RANK); 
@@ -494,15 +494,15 @@
 			}	
 		}
-		#endif
+#endif
 		/* Start  looping on the number of gaussian points: */
 		for (ig1=0; ig1<num_area_gauss; ig1++){
 			for (ig2=0; ig2<num_vert_gauss; ig2++){
-			
+
 				/*Pick up the gaussian point: */
 				gauss_weight1=*(area_gauss_weights+ig1);
 				gauss_weight2=*(vert_gauss_weights+ig2);
 				gauss_weight=gauss_weight1*gauss_weight2;
-				
-				
+
+
 				gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
 				gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
@@ -514,5 +514,5 @@
 				GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
 				GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
-			
+
 				/*Get viscosity: */
 				matice->GetViscosity3d(&viscosity, &epsilon[0]);
@@ -525,8 +525,8 @@
 				/* Get Jacobian determinant: */
 				GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
-	
+
 				/*Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant 
 				  onto this scalar matrix, so that we win some computational time: */
-				
+
 				newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
 				D_scalar=newviscosity*gauss_weight*Jdet;
@@ -534,10 +534,10 @@
 					D[i][i]=D_scalar;
 				}
-		
+
 				/*  Do the triple product tB*D*Bprime: */
 				TripleMultiply( &B[0][0],5,numdof,1,
-						&D[0][0],5,5,0,
-						&Bprime[0][0],5,numdof,0,
-						&Ke_gg_gaussian[0][0],0);
+							&D[0][0],5,5,0,
+							&Bprime[0][0],5,numdof,0,
+							&Ke_gg_gaussian[0][0],0);
 
 				/* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
@@ -549,9 +549,9 @@
 			} //for (ig2=0; ig2<num_vert_gauss; ig2++)
 		} //for (ig1=0; ig1<num_area_gauss; ig1++)
-		
+
 
 		/*Add Ke_gg to global matrix Kgg: */
 		MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
-	
+
 		//Deal with 2d friction at the bedrock interface
 		if((onbed && !shelf)){
@@ -567,6 +567,6 @@
 
 	} 
-		
-	cleanup_and_return: 
+
+cleanup_and_return: 
 	xfree((void**)&first_gauss_area_coord);
 	xfree((void**)&second_gauss_area_coord);
@@ -627,5 +627,5 @@
 	/*recover pointers: */
 	inputs=(ParameterInputs*)vinputs;
-	
+
 
 	/*If this element  is on the surface, we have a dynamic boundary condition that applies, as a stiffness 
@@ -660,5 +660,5 @@
 
 	GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
-	#ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
 	for (i=0;i<num_area_gauss;i++){
 		_printf_("Area Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(area_gauss_weights+i));
@@ -667,15 +667,15 @@
 		_printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
 	}
-	#endif
-	
+#endif
+
 	/* Start  looping on the number of gaussian points: */
 	for (ig1=0; ig1<num_area_gauss; ig1++){
 		for (ig2=0; ig2<num_vert_gauss; ig2++){
-		
+
 			/*Pick up the gaussian point: */
 			gauss_weight1=*(area_gauss_weights+ig1);
 			gauss_weight2=*(vert_gauss_weights+ig2);
 			gauss_weight=gauss_weight1*gauss_weight2;
-			
+
 			gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
 			gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
@@ -693,7 +693,7 @@
 			/*  Do the triple product tB*D*Bprime: */
 			TripleMultiply( &B[0][0],1,numgrids,1,
-					&DL_scalar,1,1,0,
-					&Bprime[0][0],1,numgrids,0,
-					&Ke_gg_gaussian[0][0],0);
+						&DL_scalar,1,1,0,
+						&Bprime[0][0],1,numgrids,0,
+						&Ke_gg_gaussian[0][0],0);
 
 			/* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
@@ -708,6 +708,6 @@
 	/*Add Ke_gg to global matrix Kgg: */
 	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
-	
-	cleanup_and_return: 
+
+cleanup_and_return: 
 	xfree((void**)&first_gauss_area_coord);
 	xfree((void**)&second_gauss_area_coord);
@@ -740,5 +740,5 @@
 	double         gravity,rho_ice,rho_water;
 
-	
+
 	/*Collapsed formulation: */
 	Tria*  tria=NULL;
@@ -746,5 +746,5 @@
 	/*Grid data: */
 	double        xyz_list[numgrids][3];
-	
+
 	/*parameters: */
 	double		   xyz_list_tria[numgrids2d][3];
@@ -770,5 +770,5 @@
 	double     DLStokes[14][14]={0.0};
 	double     tLDStokes[numdof2d][14];
-	
+
 	/* gaussian points: */
 	int     num_area_gauss;
@@ -792,5 +792,5 @@
 	double  alpha2_list[numgrids2d];
 	double  alpha2_gauss;
-	
+
 	ParameterInputs* inputs=NULL;
 
@@ -804,5 +804,5 @@
 		}
 	}
-	
+
 	/*recovre material parameters: */
 	rho_water=matpar->GetRhoWater();
@@ -818,8 +818,8 @@
 
 	/* 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 
-	   points, order of integration is 2, because we need to integrate the product tB*D*B' 
-	   which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
-	   points, same deal, which yields 3 gaussian points.*/
+		get tria gaussian points as well as segment gaussian points. For tria gaussian 
+		points, order of integration is 2, because we need to integrate the product tB*D*B' 
+		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;
@@ -845,5 +845,5 @@
 			/*Compute strain rate: */
 			GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
-		
+
 			/*Get viscosity: */
 			matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
@@ -883,9 +883,9 @@
 		/*Build alpha2_list used by drag stiffness matrix*/
 		Friction* friction=NewFriction();
-		
+
 		/*Initialize all fields: */
 		friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
 		strcpy(friction->element_type,"2d");
-		
+
 		friction->gravity=matpar->GetG();
 		friction->rho_ice=matpar->GetRhoIce();
@@ -910,5 +910,5 @@
 			}
 		}
-		
+
 		GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
 
@@ -920,25 +920,25 @@
 			gauss_coord[2]=*(third_gauss_area_coord+igarea);
 			gauss_coord[3]=-1;
-			
+
 			gauss_coord_tria[0]=*(first_gauss_area_coord+igarea); 
 			gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
 			gauss_coord_tria[2]=*(third_gauss_area_coord+igarea);
-	
+
 			/*Get the Jacobian determinant */
 			tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria);
-			
+
 			/*Get L matrix if viscous basal drag present: */
 			GetLStokes(&LStokes[0][0],  gauss_coord_tria);
 			GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss_coord_tria, gauss_coord);
-				
+
 			/*Compute strain rate: */
 			GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
-		
+
 			/*Get viscosity at last iteration: */
 			matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
-	
+
 			/*Get normal vecyor to the bed */
 			SurfaceNormal(&surface_normal[0],xyz_list_tria);
-			
+
 			bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
 			bed_normal[1]=-surface_normal[1];
@@ -974,5 +974,5 @@
 		}
 	} //if ( (onbed==1) && (shelf==0))
-	
+
 	/*Reduce the matrix */
 	ReduceMatrixStokes(&Ke_reduced[0][0], &Ke_temp[0][0]);
@@ -986,6 +986,6 @@
 	/*Add Ke_gg to global matrix Kgg: */
 	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
-	
-	cleanup_and_return: 
+
+cleanup_and_return: 
 
 	return;
@@ -998,9 +998,9 @@
 	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
 	if (analysis_type==ControlAnalysisEnum()){
-		
+
 		CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
 	}
 	else if (analysis_type==DiagnosticAnalysisEnum()){
-	
+
 		if (sub_analysis_type==HorizAnalysisEnum()){
 
@@ -1008,9 +1008,9 @@
 		}
 		else if (sub_analysis_type==VertAnalysisEnum()){
-		
+
 			CreatePVectorDiagnosticVert( pg,inputs,analysis_type,sub_analysis_type);
 		}
 		else if (sub_analysis_type==StokesAnalysisEnum()){
-		
+
 			CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type);
 		}
@@ -1018,13 +1018,13 @@
 	}
 	else if (analysis_type==SlopeComputeAnalysisEnum()){
-		
+
 		CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
 	}
 	else if (analysis_type==ThermalAnalysisEnum()){
-		
+
 		CreatePVectorThermal( pg,inputs,analysis_type,sub_analysis_type);
 	}
 	else if (analysis_type==MeltingAnalysisEnum()){
-		
+
 		CreatePVectorMelting( pg,inputs,analysis_type,sub_analysis_type);
 	}
@@ -1056,5 +1056,5 @@
 	inputs->Recover("melting",&melting[0],1,dofs,6,(void**)nodes);
 	inputs->Recover("accumulation_param",&accumulation[0],1,dofs,6,(void**)nodes);
-	
+
 	//Update material if necessary
 	if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,6,(void**)nodes)){
@@ -1063,5 +1063,5 @@
 		matice->SetB(B_average);
 	}
-	
+
 	if(inputs->Recover("B",&B_list[0],1,dofs,6,(void**)nodes)){
 		B_average=(B_list[0]+B_list[1]+B_list[2]+B_list[3]+B_list[4]+B_list[5])/6.0;
@@ -1088,5 +1088,5 @@
 	}
 }
-		
+
 int Penta::GetOnBed(){
 	return onbed;
@@ -1099,5 +1099,5 @@
 }
 void          Penta::GetBedList(double* bed_list){
-	
+
 	int i;
 	for(i=0;i<6;i++)bed_list[i]=b[i];
@@ -1115,5 +1115,5 @@
 	int i;
 	Tria* tria=NULL;
-	
+
 	/*Bail out if this element if:
 	 * -> Non collapsed and not on the surface
@@ -1132,5 +1132,5 @@
 	}
 	else{
-		
+
 		tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
 		tria->Du(du_g,u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
@@ -1156,8 +1156,8 @@
 #define __FUNCT__ "Penta::GradjDrag"
 void  Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type){
-	
-	
+
+
 	Tria* tria=NULL;
-	
+
 	/*Bail out if this element does not touch the bedrock: */
 	if (!onbed){
@@ -1165,5 +1165,5 @@
 	}
 	else{
-		
+
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 		tria->GradjDrag( grad_g,u_g,lambda_g,inputs,analysis_type,sub_analysis_type);
@@ -1178,12 +1178,12 @@
 	throw ErrorException(__FUNCT__," not supported yet!");
 }
-        
+
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::Misfit"
 double Penta::Misfit(double* velocity,double* obs_velocity,void* inputs,int analysis_type,int sub_analysis_type){
-	
+
 	double J;
 	Tria* tria=NULL;
-	
+
 	/*Bail out if this element if:
 	 * -> Non collapsed and not on the surface
@@ -1209,5 +1209,5 @@
 	}
 }
-		
+
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::SpawnTria"
@@ -1224,5 +1224,5 @@
 	double   tria_accumulation[3];
 	double   tria_geothermalflux[3];
-	
+
 	/*configuration: */
 	int tria_node_ids[3];
@@ -1249,5 +1249,5 @@
 	tria_melting[1]=melting[g1];
 	tria_melting[2]=melting[g2];
-	
+
 	tria_accumulation[0]=accumulation[g0];
 	tria_accumulation[1]=accumulation[g1];
@@ -1286,5 +1286,5 @@
 	int doflist_per_node[MAXDOFSPERNODE];
 	int numberofdofspernode;
-	
+
 	for(i=0;i<6;i++){
 		nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
@@ -1320,5 +1320,5 @@
 	GetB(&B[0][0], xyz_list, gauss_l1l2l3l4);
 
-	#ifdef _DEBUG_
+#ifdef _DEBUG_
 	_printf_("B for grid1 : [ %lf   %lf  \n",B[0][0],B[0][1]);
 	_printf_("              [ %lf   %lf  \n",B[1][0],B[1][1]);
@@ -1326,5 +1326,5 @@
 	_printf_("              [ %lf   %lf  ]\n",B[3][0],B[3][1]);
 	_printf_("              [ %lf   %lf  ]\n",B[4][0],B[4][1]);
-	
+
 	_printf_("B for grid2 : [ %lf   %lf  \n",B[0][2],B[0][3]);
 	_printf_("              [ %lf   %lf  \n",B[1][2],B[1][3]);
@@ -1332,5 +1332,5 @@
 	_printf_("              [ %lf   %lf  ]\n",B[3][2],B[3][3]);
 	_printf_("              [ %lf   %lf  ]\n",B[4][2],B[4][3]);
-	
+
 	_printf_("B for grid3 : [ %lf   %lf  \n", B[0][4],B[0][5]);
 	_printf_("              [ %lf   %lf  \n", B[1][4],B[1][5]);
@@ -1338,5 +1338,5 @@
 	_printf_("              [ %lf   %lf  ]\n",B[3][4],B[3][5]);
 	_printf_("              [ %lf   %lf  ]\n",B[4][4],B[4][5]);
-	
+
 	_printf_("B for grid4 : [ %lf   %lf  \n", B[0][6],B[0][7]);
 	_printf_("              [ %lf   %lf  \n", B[1][6],B[1][7]);
@@ -1344,5 +1344,5 @@
 	_printf_("              [ %lf   %lf  ]\n",B[3][6],B[3][7]);
 	_printf_("              [ %lf   %lf  ]\n",B[4][6],B[4][7]);
-				
+
 	_printf_("B for grid5 : [ %lf   %lf  \n", B[0][8],B[0][9]);
 	_printf_("              [ %lf   %lf  \n", B[1][8],B[1][9]);
@@ -1360,10 +1360,10 @@
 		_printf_("Velocity for grid %i %lf %lf\n",i,*(vxvy_list+2*i+0),*(vxvy_list+2*i+1));
 	}
-	#endif
+#endif
 
 	/*Multiply B by velocity, to get strain rate: */
 	MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
-			              velocity,NDOF2*numgrids,1,0,
-						  epsilon,0);
+				velocity,NDOF2*numgrids,1,0,
+				epsilon,0);
 
 }
@@ -1385,5 +1385,5 @@
 	 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
 	 */
-	
+
 	int i;
 	const int numgrids=6;
@@ -1396,9 +1396,9 @@
 	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
 
-	#ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
 	for (i=0;i<numgrids;i++){
 		_printf_("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
 	}
-	#endif
+#endif
 
 	/*Build B: */
@@ -1415,5 +1415,5 @@
 		*(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i]; 
 		*(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
-		
+
 		*(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
 		*(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i]; 
@@ -1439,5 +1439,5 @@
 	 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
 	 */
-	
+
 	int i;
 	const int NDOF3=3;
@@ -1450,9 +1450,9 @@
 	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
 
-	#ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
 	for (i=0;i<numgrids;i++){
 		_printf_("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
 	}
-	#endif
+#endif
 
 	/*Build BPrime: */
@@ -1469,5 +1469,5 @@
 		*(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[2][i]; 
 		*(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
-		
+
 		*(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
 		*(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[2][i]; 
@@ -1482,5 +1482,5 @@
 	 * the determinant of it: */
 	const int NDOF3=3;
-	
+
 	double J[NDOF3][NDOF3];
 
@@ -1496,8 +1496,8 @@
 #define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasic" 
 void Penta::GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3dh4dh5dh6_basic,double* xyz_list, double* gauss_l1l2l3l4){
-	
+
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the basic coordinate system: */
 
-	
+
 	int i;
 	const int NDOF3=3;
@@ -1534,5 +1534,5 @@
 	const int NDOF3=3;
 	int i,j;
-	
+
 	/*The Jacobian is constant over the element, discard the gaussian points. 
 	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
@@ -1544,7 +1544,7 @@
 	double y1,y2,y3,y4,y5,y6;
 	double z1,z2,z3,z4,z5,z6;
-	
+
 	double sqrt3=sqrt(3.0);
-	
+
 	/*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
 	A1=gauss_l1l2l3l4[0];
@@ -1562,5 +1562,5 @@
 	x5=*(xyz_list+3*4+0);
 	x6=*(xyz_list+3*5+0);
-	
+
 	y1=*(xyz_list+3*0+1);
 	y2=*(xyz_list+3*1+1);
@@ -1590,5 +1590,5 @@
 	*(J+NDOF3*2+2)=sqrt3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+1.0/4.0*(z1-z2-z4+z5)*xi+1.0/4.0*(-z1+z5-z2+z4);
 
-	#ifdef _DEBUG_
+#ifdef _DEBUG_
 	for(i=0;i<3;i++){
 		for (j=0;j<3;j++){
@@ -1597,5 +1597,5 @@
 		printf("\n");
 	}
-	#endif
+#endif
 }
 
@@ -1603,5 +1603,5 @@
 #define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParams" 
 void Penta::GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4){
-	
+
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
 	 * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */
@@ -1610,5 +1610,5 @@
 	double A1,A2,A3,z;
 	double sqrt3=sqrt(3.0);
-	
+
 	A1=gauss_l1l2l3l4[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*sqrt(3));
 	A2=gauss_l1l2l3l4[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*sqrt(3));
@@ -1677,5 +1677,5 @@
 	int          doflist[numdof];
 	int          numberofdofspernode;
-	
+
 	/* parameters: */
 	double  slope[NDOF2];
@@ -1752,5 +1752,5 @@
 
 		GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
-		#ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
 		for (i=0;i<num_area_gauss;i++){
 			_printf_("Area Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(area_gauss_weights+i));
@@ -1759,20 +1759,20 @@
 			_printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
 		}
-		#endif
-	
+#endif
+
 		/* Start  looping on the number of gaussian points: */
 		for (ig1=0; ig1<num_area_gauss; ig1++){
 			for (ig2=0; ig2<num_vert_gauss; ig2++){
-			
+
 				/*Pick up the gaussian point: */
 				gauss_weight1=*(area_gauss_weights+ig1);
 				gauss_weight2=*(vert_gauss_weights+ig2);
 				gauss_weight=gauss_weight1*gauss_weight2;
-				
+
 				gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
 				gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
 				gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
 				gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
-		
+
 				/*Compute thickness at gaussian point: */
 				GetParameterValue(&thickness, &h[0],gauss_l1l2l3l4);
@@ -1783,6 +1783,6 @@
 				/* Get Jacobian determinant: */
 				GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
-		
-				 /*Get nodal functions: */
+
+				/*Get nodal functions: */
 				GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
 
@@ -1808,5 +1808,5 @@
 	VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
 
-	cleanup_and_return: 
+cleanup_and_return: 
 	xfree((void**)&first_gauss_area_coord);
 	xfree((void**)&second_gauss_area_coord);
@@ -1822,5 +1822,5 @@
 #define __FUNCT__ "Penta::CreateKMatrix"
 void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4){
-	
+
 	const int numgrids=6;
 	double l1l2l3l4l5l6[numgrids];
@@ -1834,5 +1834,5 @@
 #define __FUNCT__ "Penta::GetParameterDerivativeValue"
 void Penta::GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4){
-				
+
 	/*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_l1l2l3l4:
 	 *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
@@ -1849,7 +1849,7 @@
 	/*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
 	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
-	
+
 	*(p+0)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[0][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[0][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[0][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[0][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[0][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[0][5];
-;
+	;
 	*(p+1)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[1][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[1][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[1][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[1][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[1][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[1][5];
 
@@ -1861,5 +1861,5 @@
 #define __FUNCT__ "Penta::GetNodalFunctions"
 void Penta::GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4){
-	
+
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
 
@@ -1867,5 +1867,5 @@
 
 	l1l2l3l4l5l6[1]=gauss_l1l2l3l4[1]*(1-gauss_l1l2l3l4[3])/2.0;
-	
+
 	l1l2l3l4l5l6[2]=gauss_l1l2l3l4[2]*(1-gauss_l1l2l3l4[3])/2.0;
 
@@ -1873,5 +1873,5 @@
 
 	l1l2l3l4l5l6[4]=gauss_l1l2l3l4[1]*(1+gauss_l1l2l3l4[3])/2.0;
-	
+
 	l1l2l3l4l5l6[5]=gauss_l1l2l3l4[2]*(1+gauss_l1l2l3l4[3])/2.0;
 
@@ -1888,12 +1888,12 @@
 	int          nodedofs[2];
 	int          numberofdofspernode;
-	
+
 	Node* node=NULL;
 	int   i;
 	double velocity[2];
-		
+
 
 	if((collapse==1) && (onbed==1)){
-			
+
 		GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1902,7 +1902,7 @@
 		 * inserting the same velocity value into ug, until we reach the surface: */
 		for(i=0;i<3;i++){
-	
+
 			node=nodes[i]; //base nodes
-	
+
 			/*get velocity for this base node: */
 			velocity[0]=ug_serial[doflist[numberofdofspernode*i+0]];
@@ -1938,12 +1938,12 @@
 	int          nodedof;
 	int          numberofdofspernode;
-	
+
 	Node* node=NULL;
 	int   i;
 	double slope;
-		
+
 
 	if(onbed==1){
-			
+
 		GetDofList(&doflist[0],&numberofdofspernode);
 
@@ -1952,9 +1952,9 @@
 		/*For each node on the base of this penta,  we grab the slope. Once we know the slope, we follow the upper nodes, 
 		 * inserting the same slope value into sg, until we reach the surface: */
-		
+
 		for(i=0;i<3;i++){
-	
+
 			node=nodes[i]; //base nodes
-	
+
 			/*get velocity for this base node: */
 			slope=sg_serial[doflist[i]];
@@ -1984,5 +1984,5 @@
 
 	/*	Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
-	where hi is the interpolation function for grid i.*/
+		where hi is the interpolation function for grid i.*/
 
 	int i;
@@ -1994,9 +1994,9 @@
 	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
 
-	#ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
 	for (i=0;i<numgrids;i++){
 		_printf_("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
 	}
-	#endif
+#endif
 
 	/*Build B: */
@@ -2004,5 +2004,5 @@
 		B[i]=dh1dh2dh3dh4dh5dh6_basic[2][i];  
 	}
-	
+
 }
 
@@ -2015,5 +2015,5 @@
 
 	int i;
-				
+
 	GetNodalFunctions(B, gauss_l1l2l3l4);
 
@@ -2034,5 +2034,5 @@
 	int          doflist[numdof];
 	int          numberofdofspernode;
-	
+
 	/* gaussian points: */
 	int     num_gauss,ig;
@@ -2094,5 +2094,5 @@
 	/* recover input parameters: */
 	if(!inputs->Recover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity!");
-	    inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);
+	inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);
 
 	/*Get gaussian points and weights :*/
@@ -2101,5 +2101,5 @@
 
 	GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
-	#ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
 	for (i=0;i<num_area_gauss;i++){
 		_printf_("Area Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(area_gauss_weights+i));
@@ -2108,20 +2108,20 @@
 		_printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
 	}
-	#endif
+#endif
 
 	/* Start  looping on the number of gaussian points: */
 	for (ig1=0; ig1<num_area_gauss; ig1++){
 		for (ig2=0; ig2<num_vert_gauss; ig2++){
-		
+
 			/*Pick up the gaussian point: */
 			gauss_weight1=*(area_gauss_weights+ig1);
 			gauss_weight2=*(vert_gauss_weights+ig2);
 			gauss_weight=gauss_weight1*gauss_weight2;
-			
+
 			gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1); 
 			gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
 			gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
 			gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
-	
+
 			/*Get velocity derivative, with respect to x and y: */
 			GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3l4);
@@ -2129,13 +2129,13 @@
 			dudx=du[0];
 			dvdy=dv[1];
-			
+
 
 			/* Get Jacobian determinant: */
 			GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
-			#ifdef _DEBUG_ 
+#ifdef _DEBUG_ 
 			_printf_("Element id %i Jacobian determinant: %lf\n",PentaElementGetID(this),Jdet);
-			#endif
-		
-			 /*Get nodal functions: */
+#endif
+
+			/*Get nodal functions: */
 			GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
 
@@ -2154,5 +2154,5 @@
 	VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
 
-	cleanup_and_return: 
+cleanup_and_return: 
 	xfree((void**)&first_gauss_area_coord);
 	xfree((void**)&second_gauss_area_coord);
@@ -2175,8 +2175,8 @@
 	double rho_ice,g;
 	double       xyz_list[numgrids][3];
-		
+
 	/*Get node data: */
 	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
-        
+
 	/*pressure is lithostatic: */
 	//md.pressure=md.rho_ice*md.g*(md.surface-md.z); a la matlab
@@ -2191,5 +2191,5 @@
 		pressure[i]=rho_ice*g*(s[i]-xyz_list[i][2]);
 	}
-	
+
 	/*plug local pressure values into global pressure vector: */
 	VecSetValues(pg,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
@@ -2205,5 +2205,5 @@
 	/*Collapsed formulation: */
 	Tria*  tria=NULL;
-	
+
 	/*Is this element on the bed? :*/
 	if(!onbed)return;
@@ -2221,8 +2221,8 @@
 
 void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
-	
+
 	/*Collapsed formulation: */
 	Tria*  tria=NULL;
-	
+
 	/*Is this element on the bed? :*/
 	if(!onbed)return;
@@ -2238,5 +2238,5 @@
 #define __FUNCT__ "ReduceMatrixStokes" 
 void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){
-		
+
 	int i,j;
 
@@ -2292,5 +2292,5 @@
 	double det;
 	int calculationdof=3;
-	
+
 	/*Take the matrix components: */
 	a=*(Ke+calculationdof*0+0);
@@ -2305,5 +2305,5 @@
 
 	det=a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g);
-	
+
 	*(Ke_invert+calculationdof*0+0)=(e*i-f*h)/det;
 	*(Ke_invert+calculationdof*0+1)=(c*h-b*i)/det;
@@ -2392,15 +2392,15 @@
 	 * For grid i, Bi can be expressed in the basic coordinate system
 	 * by: 				Bi_basic=[ dh/dx          0             0       0  ]
-		*					[   0           dh/dy           0       0  ]
-		*					[   0             0           dh/dy     0  ]
-		*					[ 1/2*dh/dy    1/2*dh/dx        0       0  ]
-		*					[ 1/2*dh/dz       0         1/2*dh/dx   0  ]
-		*					[   0          1/2*dh/dz    1/2*dh/dy   0  ]
-		*					[   0             0             0       h  ]
-		*					[ dh/dx         dh/dy         dh/dz     0  ]
-		*	where h is the interpolation function for grid i.
-		*	Same thing for Bb except the last column that does not exist.
-		*/
-	
+	 *					[   0           dh/dy           0       0  ]
+	 *					[   0             0           dh/dy     0  ]
+	 *					[ 1/2*dh/dy    1/2*dh/dx        0       0  ]
+	 *					[ 1/2*dh/dz       0         1/2*dh/dx   0  ]
+	 *					[   0          1/2*dh/dz    1/2*dh/dy   0  ]
+	 *					[   0             0             0       h  ]
+	 *					[ dh/dx         dh/dy         dh/dz     0  ]
+	 *	where h is the interpolation function for grid i.
+	 *	Same thing for Bb except the last column that does not exist.
+	 */
+
 	int i;
 	int calculationdof=3;
@@ -2416,11 +2416,11 @@
 
 	GetNodalFunctions(l1l6, gauss_coord);
-	
-	#ifdef _DEBUG_ 
+
+#ifdef _DEBUG_ 
 	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]);
 	}
 
-	#endif
+#endif
 
 	/*Build B: */
@@ -2470,18 +2470,18 @@
 
 	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 
-	*	For grid i, Bi' can be expressed in the basic coordinate system
-	*	by: 
-	*				Bi_basic'=[ dh/dx   0          0       0]
-	*					 [   0      dh/dy      0       0]
-	*					 [   0      0         dh/dz    0]
-	*					 [  dh/dy   dh/dx      0       0]
-	*					 [  dh/dz   0        dh/dx     0]
-	*					 [   0      dh/dz    dh/dy     0]
-	*					 [  dh/dx   dh/dy    dh/dz     0]
-	*					 [   0      0          0       h]
-	*	where h is the interpolation function for grid i.
-	*
-	* 	Same thing for the bubble fonction except that there is no fourth column
-	*/
+	 *	For grid i, Bi' can be expressed in the basic coordinate system
+	 *	by: 
+	 *				Bi_basic'=[ dh/dx   0          0       0]
+	 *					 [   0      dh/dy      0       0]
+	 *					 [   0      0         dh/dz    0]
+	 *					 [  dh/dy   dh/dx      0       0]
+	 *					 [  dh/dz   0        dh/dx     0]
+	 *					 [   0      dh/dz    dh/dy     0]
+	 *					 [  dh/dx   dh/dy    dh/dz     0]
+	 *					 [   0      0          0       h]
+	 *	where h is the interpolation function for grid i.
+	 *
+	 * 	Same thing for the bubble fonction except that there is no fourth column
+	 */
 
 	int i;
@@ -2498,11 +2498,11 @@
 
 	GetNodalFunctions(l1l6, gauss_coord);
-	
-	#ifdef _DEBUG_ 
+
+#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]);
 	}
 
-	#endif
+#endif
 
 	/*B_primeuild B_prime: */
@@ -2551,23 +2551,23 @@
 
 	/*
-	* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
-	* For grid i, Li can be expressed in the basic coordinate system
-	* by: 
-	*       Li_basic=[ h    0    0   0]
-	*	 	 [ 0    h    0   0]
-	*		 [ 0    0    h   0]
-	*		 [ 0    0    h   0]
-	*	 	 [ h    0    0   0]
-	*	 	 [ 0    h    0   0]
-	*	 	 [ h    0    0   0]
-	*	 	 [ 0    h    0   0]
-	*		 [ 0    0    h   0]
-	*		 [ 0    0    h   0]
-	*		 [ 0    0    h   0]
-	*	 	 [ h    0    0   0]
-	*	 	 [ 0    h    0   0]
-	*		 [ 0    0    h   0]
-	* where h is the interpolation function for grid i.
-	*/
+	 * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+	 * For grid i, Li can be expressed in the basic coordinate system
+	 * by: 
+	 *       Li_basic=[ h    0    0   0]
+	 *	 	 [ 0    h    0   0]
+	 *		 [ 0    0    h   0]
+	 *		 [ 0    0    h   0]
+	 *	 	 [ h    0    0   0]
+	 *	 	 [ 0    h    0   0]
+	 *	 	 [ h    0    0   0]
+	 *	 	 [ 0    h    0   0]
+	 *		 [ 0    0    h   0]
+	 *		 [ 0    0    h   0]
+	 *		 [ 0    0    h   0]
+	 *	 	 [ h    0    0   0]
+	 *	 	 [ 0    h    0   0]
+	 *		 [ 0    0    h   0]
+	 * where h is the interpolation function for grid i.
+	 */
 
 	int i;
@@ -2582,11 +2582,11 @@
 	l1l2l3[1]=gauss_coord_tria[1];
 	l1l2l3[2]=gauss_coord_tria[2];
-	
-
-	#ifdef _DELUG_ 
+
+
+#ifdef _DELUG_ 
 	for (i=0;i<3;i++){
 		printf("Node %i  h=%lf \n",i,l1l2l3[i]);
 	}
-	#endif
+#endif
 
 	/*Build LStokes: */
@@ -2648,5 +2648,5 @@
 		*(LStokes+num_dof*numgrids2d*13+num_dof*i+2)=l1l2l3[i];
 		*(LStokes+num_dof*numgrids2d*13+num_dof*i+3)=0;
-	
+
 	}
 }
@@ -2658,23 +2658,23 @@
 
 	/*
-	* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
-	* For grid i, Lpi can be expressed in the basic coordinate system
-	* by: 
-	*       Lpi_basic=[ h    0    0   0]
-	*		 [ 0    h    0   0]
-	*		 [ h    0    0   0]
-	*		 [ 0    h    0   0]
-	*		 [ 0    0    h   0]
-	*		 [ 0    0    h   0]
-	*		 [ 0    0  dh/dz 0]
-	*		 [ 0    0  dh/dz 0]
-	*		 [ 0    0  dh/dz 0]
-	*		 [dh/dz 0  dh/dx 0]
-	*		 [ 0 dh/dz dh/dy 0]
-	* 		 [ 0    0    0   h]
-	* 		 [ 0    0    0   h]
-	* 		 [ 0    0    0   h]
-	* where h is the interpolation function for grid i.
-	*/
+	 * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
+	 * For grid i, Lpi can be expressed in the basic coordinate system
+	 * by: 
+	 *       Lpi_basic=[ h    0    0   0]
+	 *		 [ 0    h    0   0]
+	 *		 [ h    0    0   0]
+	 *		 [ 0    h    0   0]
+	 *		 [ 0    0    h   0]
+	 *		 [ 0    0    h   0]
+	 *		 [ 0    0  dh/dz 0]
+	 *		 [ 0    0  dh/dz 0]
+	 *		 [ 0    0  dh/dz 0]
+	 *		 [dh/dz 0  dh/dx 0]
+	 *		 [ 0 dh/dz dh/dy 0]
+	 * 		 [ 0    0    0   h]
+	 * 		 [ 0    0    0   h]
+	 * 		 [ 0    0    0   h]
+	 * where h is the interpolation function for grid i.
+	 */
 
 	int i;
@@ -2690,12 +2690,12 @@
 	l1l2l3[1]=gauss_coord_tria[1];
 	l1l2l3[2]=gauss_coord_tria[2];
-	
+
 	GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
 
-	#ifdef _DELUG_ 
+#ifdef _DELUG_ 
 	for (i=0;i<3;i++){
 		printf("Node %i  h=%lf \n",i,l1l2l3[i]);
 	}
-	#endif
+#endif
 
 	/*Build LprimeStokes: */
@@ -2757,7 +2757,7 @@
 		*(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0;
 		*(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i];
-	
-	}
-	
+
+	}
+
 }
 
@@ -2765,5 +2765,5 @@
 #define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasicStokes"
 void Penta::GetNodalFunctionsDerivativesBasicStokes(double* dh1dh7_basic,double* xyz_list, double* gauss_coord){
-	
+
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
 	 * basic coordinate system: */
@@ -2801,5 +2801,5 @@
 #define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParamsStokes"
 void Penta::GetNodalFunctionsDerivativesParamsStokes(double* dl1dl7,double* gauss_coord){
-	
+
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
 	 * natural coordinate system) at the gaussian point. */
@@ -2893,5 +2893,5 @@
 	int     area_order=5;
 	int	  num_vert_gauss=5;
-	
+
 	double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
 	double  viscosity;
@@ -2916,5 +2916,5 @@
 	double     tBD[27][8];
 	double     P_terms[numdof];
-	
+
 	ParameterInputs* inputs=NULL;
 	Tria*            tria=NULL;
@@ -2941,8 +2941,8 @@
 
 	/* 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 
-	   points, order of integration is 2, because we need to integrate the product tB*D*B' 
-	   which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
-	   points, same deal, which yields 3 gaussian points.*/
+		get tria gaussian points as well as segment gaussian points. For tria gaussian 
+		points, order of integration is 2, because we need to integrate the product tB*D*B' 
+		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): */
@@ -2964,5 +2964,5 @@
 			GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
 			matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
-		
+
 			/* Get Jacobian determinant: */
 			GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
@@ -2984,5 +2984,5 @@
 			GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord); 
 			GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord); 
-			
+
 			/*Get bubble part of Bprime */
 			for(i=0;i<8;i++){
@@ -3033,5 +3033,5 @@
 			gauss_coord[2]=*(third_gauss_area_coord+igarea);
 			gauss_coord[3]=-1;
-			
+
 			gauss_coord_tria[0]=*(first_gauss_area_coord+igarea); 
 			gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
@@ -3040,5 +3040,5 @@
 			/*Get the Jacobian determinant */
 			tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria);
-			
+
 			/* Get bed at gaussian point */
 			GetParameterValue(&bed,&b[0],gauss_coord);
@@ -3046,5 +3046,5 @@
 			/*Get L matrix : */
 			tria->GetL(&L[0], &xyz_list[0][0], gauss_coord_tria,1);
-				
+
 			/*Get water_pressure at gaussian point */
 			water_pressure=gravity*rho_water*bed;
@@ -3052,5 +3052,5 @@
 			/*Get normal vecyor to the bed */
 			SurfaceNormal(&surface_normal[0],xyz_list_tria);
-			
+
 			bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
 			bed_normal[1]=-surface_normal[1];
@@ -3064,5 +3064,5 @@
 		}
 	} //if ( (onbed==1) && (shelf==1))
-	
+
 	/*Reduce the matrix */
 	ReduceVectorStokes(&Pe_reduced[0], &Ke_temp[0][0], &Pe_temp[0]);
@@ -3080,5 +3080,5 @@
 #define __FUNCT__ "Penta::ReduceVectorStokes" 
 void Penta::ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp){
-		
+
 	int i,j;
 
@@ -3123,5 +3123,5 @@
 #define __FUNCT__ "Penta::GetNodalFunctionsStokes" 
 void Penta::GetNodalFunctionsStokes(double* l1l7, double* gauss_coord){
-	
+
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
 
@@ -3161,7 +3161,7 @@
 	const int    NDOF1=1;
 	const int    numdof=NDOF1*numgrids;
-	double       xyz_list[numgrids][3];
-	int          doflist[numdof];
-	int          numberofdofspernode;
+	double xyz_list[numgrids][3];
+	int    doflist[numdof];
+	int    numberofdofspernode;
 
 	/* gaussian points: */
@@ -3177,52 +3177,47 @@
 
 	int area_order=5;
-	int	num_vert_gauss=5;
-	
-	int            dofs[3]={0,1,2};
-	double         dt;
-	int            dt_exists;
-	
-	double         vxvyvz_list[numgrids][3];
-	double         vx_list[numgrids];
-	int            vx_list_exists;
-	double         u;
-	double         vy_list[numgrids];
-	int            vy_list_exists;
-	double         v;
-	double         vz_list[numgrids];
-	int            vz_list_exists;
-	double         w;
-
-	/* element parameters: */
-	int	    onbed;
-	int	    shelf;
-	int     steady_state;
+	int num_vert_gauss=5;
+
+	int     dofs[3]={0,1,2};
+	double  dt;
+	int     dt_exists;
+
+	double  vxvyvz_list[numgrids][3];
+	double  vx_list[numgrids];
+	int     vx_list_exists;
+	double  u;
+	double  vy_list[numgrids];
+	int     vy_list_exists;
+	double  v;
+	double  vz_list[numgrids];
+	int     vz_list_exists;
+	double  w;
 
 	/*matrices: */
-	double     K_terms[numdof][numdof]={0.0};
-	double     Ke_gaussian_conduct[numdof][numdof];
-	double     Ke_gaussian_advec[numdof][numdof];
-	double     Ke_gaussian_transient[numdof][numdof];
-	double     B[3][numdof];
-	double     Bprime[3][numdof];
-	double     B_conduct[3][numdof];
-	double     B_advec[3][numdof];
-	double     Bprime_advec[3][numdof];
-	double     L[numdof];
-	double     D_scalar;
-	double     D[3][3];
-	double     l1l2l3[3];
-	double     tl1l2l3D[3];
-	double     tBD[3][numdof];
-	double     tBD_conduct[3][numdof];
-	double     tBD_advec[3][numdof];
-	double     tLD[numdof];
-
-	double     Jdet;
+	double  K_terms[numdof][numdof]={0.0};
+	double  Ke_gaussian_conduct[numdof][numdof];
+	double  Ke_gaussian_advec[numdof][numdof];
+	double  Ke_gaussian_transient[numdof][numdof];
+	double  B[3][numdof];
+	double  Bprime[3][numdof];
+	double  B_conduct[3][numdof];
+	double  B_advec[3][numdof];
+	double  Bprime_advec[3][numdof];
+	double  L[numdof];
+	double  D_scalar;
+	double  D[3][3];
+	double  l1l2l3[3];
+	double  tl1l2l3D[3];
+	double  tBD[3][numdof];
+	double  tBD_conduct[3][numdof];
+	double  tBD_advec[3][numdof];
+	double  tLD[numdof];
+
+	double  Jdet;
 
 	/*Material properties: */
-	double         gravity,rho_ice,rho_water;
-	double         heatcapacity,thermalconductivity;
-	double         mixed_layer_capacity,thermal_exchange_velocity;
+	double  gravity,rho_ice,rho_water;
+	double  heatcapacity,thermalconductivity;
+	double  mixed_layer_capacity,thermal_exchange_velocity;
 
 	/*Collapsed formulation: */
@@ -3237,5 +3232,5 @@
 	GetDofList(&doflist[0],&numberofdofspernode);
 	
-	// /*recovre material parameters: */
+	/*recovre material parameters: */
 	rho_water=matpar->GetRhoWater();
 	rho_ice=matpar->GetRhoIce();
@@ -3243,5 +3238,4 @@
 	heatcapacity=matpar->GetHeatCapacity();
 	thermalconductivity=matpar->GetThermalConductivity();
-
 		
 	/*recover extra inputs from users, dt and velocity: */
@@ -3300,5 +3294,4 @@
 			D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar;
 
-
 			/*  Do the triple product B'*D*B: */
 			MatrixMultiply(&B_conduct[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_conduct[0][0],0);
@@ -3310,5 +3303,4 @@
 			GetB_advec(&B_advec[0][0],&xyz_list[0][0],gauss_coord); 
 			GetBprime_advec(&Bprime_advec[0][0],&xyz_list[0][0],gauss_coord); 
-
 
 			//Build the D matrix
@@ -3330,5 +3322,4 @@
 			MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0);
 			MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0);
-
 
 			/*Transient: */
@@ -3359,14 +3350,4 @@
 		}
 	}
-	
-	//Ice/ocean heat exchange flux on ice shelf base 
-
-	if(onbed && shelf){
-
-		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-		tria->CreateKMatrixThermal(Kgg,inputs, analysis_type,sub_analysis_type);
-		delete tria;
-		return;
-	}
 
 	/*Add Ke_gg to global matrix Kgg: */
@@ -3381,4 +3362,11 @@
 	xfree((void**)&vert_gauss_coord);
 
+	//Ice/ocean heat exchange flux on ice shelf base 
+	if(onbed && shelf){
+
+		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+		tria->CreateKMatrixThermal(Kgg,inputs, analysis_type,sub_analysis_type);
+		delete tria;
+	}
 }
 
@@ -3397,5 +3385,5 @@
 	 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
 	 */
-	
+
 	int i;
 	int calculationdof=3;
@@ -3431,5 +3419,5 @@
 	 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
 	 */
-	
+
 	int i;
 	int calculationdof=3;
@@ -3466,5 +3454,5 @@
 	 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
 	 */
-	
+
 	int i;
 	int calculationdof=3;
@@ -3489,5 +3477,5 @@
 #define __FUNCT__ "Penta::CreateKMatrixMelting"
 void  Penta::CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
-	
+
 	Tria* tria=NULL;
 	if (!onbed){
@@ -3495,5 +3483,5 @@
 	}
 	else{
-		
+
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 		tria->CreateKMatrixMelting(Kgg,inputs, analysis_type,sub_analysis_type);
@@ -3520,5 +3508,5 @@
 	/*Grid data: */
 	double        xyz_list[numgrids][3];
-	
+
 	/* gaussian points: */
 	int     num_area_gauss,igarea,igvert;
@@ -3533,5 +3521,5 @@
 	int area_order=2;
 	int	num_vert_gauss=3;
-	
+
 	double         dt;
 	double         vx_list[numgrids];
@@ -3549,5 +3537,5 @@
 	/* element parameters: */
 	int     friction_type;
-	
+
 	int            dofs[3]={0,1,2};
 	int            dofs1[1]={0};
@@ -3584,5 +3572,5 @@
 	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
 	GetDofList(&doflist[0],&numberofdofspernode);
-	
+
 	// /*recovre material parameters: */
 	rho_water=matpar->GetRhoWater();
@@ -3604,5 +3592,5 @@
 		if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
 	}
-	
+
 	for(i=0;i<numgrids;i++){
 		vx_list[i]=vxvyvz_list[i][0];
@@ -3612,9 +3600,9 @@
 
 	/* 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 
-	   points, order of integration is 2, because we need to integrate the product tB*D*B' 
-	   which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
-	   points, same deal, which yields 3 gaussian points.: */
-	
+		get tria gaussian points as well as segment gaussian points. For tria gaussian 
+		points, order of integration is 2, because we need to integrate the product tB*D*B' 
+		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: */
 	GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
@@ -3631,5 +3619,5 @@
 			gauss_coord[2]=*(third_gauss_area_coord+igarea);
 			gauss_coord[3]=*(vert_gauss_coord+igvert);
-	
+
 			/*Compute strain rate and viscosity: */
 			GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
@@ -3673,5 +3661,4 @@
 	/* Ice/ocean heat exchange flux on ice shelf base */
 	if(onbed && shelf){
-
 
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 510)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 511)
@@ -16,5 +16,4 @@
 #include "../DataSet/DataSet.h"
 #include "../include/typedefs.h"
-
 
 /*For debugging purposes: */
@@ -51,5 +50,6 @@
 		melting[i]=tria_melting[i];
 		accumulation[i]=tria_accumulation[i];
-		geothermalflux[i]=tria_geothermalflux[i]; }
+		geothermalflux[i]=tria_geothermalflux[i];
+	}
 	matice=NULL;
 	matice_offset=UNDEF;
@@ -242,5 +242,4 @@
 }
 
-
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::Configure"
@@ -294,7 +293,5 @@
 		throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
 	}
-
-}
-
+}
 
 #undef __FUNCT__ 
@@ -302,5 +299,4 @@
 
 void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
-
 
 	/* local declarations */