Index: /issm/trunk-jpl/src/c/bamg/Geometry.cpp
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Geometry.cpp	(revision 16165)
+++ /issm/trunk-jpl/src/c/bamg/Geometry.cpp	(revision 16166)
@@ -499,5 +499,4 @@
 		double        eps      = 1e-20;
 		BamgQuadtree  quadtree;                            // build quadtree to find duplicates
-		BamgVertex   *v0       = vertices;
 
 		k=0;
Index: /issm/trunk-jpl/src/c/bamg/Mesh.cpp
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 16165)
+++ /issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 16166)
@@ -522,7 +522,9 @@
 				head=(int)bamgmesh->SubDomains[i*3+1]-1;//C indexing
 				direction=(int)bamgmesh->SubDomains[i*3+2];
-				if (i3!=23) _error_("Bad Subdomain definition: first number should be 3");
+				if (i3!=3) _error_("Bad Subdomain definition: first number should be 3");
 				if (head<0 || head>=nbt) _error_("Bad Subdomain definition: head should in [1 " << nbt << "] (triangle number)");
 				subdomains[i].head = triangles+head;
+				subdomains[i].direction = direction;
+				subdomains[i].ReferenceNumber = i3;
 			}
 		}
@@ -1051,7 +1053,7 @@
 
 		/*Get options*/
-		int    verbose=bamgopts->verbose;
-		double anisomax =bamgopts->anisomax;
-		double errg     =bamgopts->errg;
+		int    verbose  = bamgopts->verbose;
+		double anisomax = bamgopts->anisomax;
+		double errg     = bamgopts->errg;
 
 		double ss[2]={0.00001,0.99999};
@@ -1061,10 +1063,8 @@
 
 		//check that hmax is positive
-		if (hmax<=0){
-			_error_("hmax<=0");
-		}
+		_assert_(hmax>0); 
 
 		//errC cannot be higher than 1
-		if (errC>1) errC=1;
+		if(errC>1) errC=1;
 
 		//Set all vertices to "on"
@@ -1072,5 +1072,5 @@
 
 		//loop over all the vertices on edges
-		for (int  i=0;i<nbe;i++){
+		for(int  i=0;i<nbe;i++){
 			for (int j=0;j<2;j++){
 
@@ -2621,6 +2621,5 @@
 
 			  }
-			else
-			  { // find the head for all sub domaine
+			else{ // find the head for all subdomains
 				if (Gh.nbsubdomains != nbsubdomains && subdomains)
 				 delete [] subdomains, subdomains=0;
@@ -2628,5 +2627,4 @@
 				 subdomains = new SubDomain[ Gh.nbsubdomains];
 				nbsubdomains =Gh.nbsubdomains;
-				long err=0;
 				CreateSingleVertexToTriangleConnectivity();
 				long * mark = new long[nbt];
@@ -2639,5 +2637,4 @@
 					GeomEdge &eg = *Gh.subdomains[i].edge;
 					subdomains[i].ReferenceNumber = Gh.subdomains[i].ReferenceNumber;
-					int ssdlab = subdomains[i].ReferenceNumber;
 					// by carefull is not easy to find a edge create from a GeomEdge 
 					// see routine MakeGeomEdgeToEdge
@@ -3109,8 +3106,5 @@
 	void Mesh::MakeBamgQuadtree() {  
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeBamgQuadtree)*/
-
-		long int verbose=0;
-		if (  !quadtree )  quadtree = new BamgQuadtree(this);
-
+		if(!quadtree) quadtree = new BamgQuadtree(this);
 	}
 	/*}}}*/
@@ -3314,16 +3308,11 @@
 			}
 
-		} while (nbv!=nbvold);
-
-		delete []  first_np_or_next_t;
-
-		long NbSwapf =0,NbSwp;
-
-		NbSwp = NbSwapf;
-		for (i=0;i<nbv;i++)
-		 NbSwapf += vertices[i].Optim(0);
-		NbTSwap +=  NbSwapf ;
-	}
-	/*}}}*/
+		}while(nbv!=nbvold);
+		delete [] first_np_or_next_t;
+
+		long NbSwapf =0;
+
+		for(i=0;i<nbv;i++) NbSwapf += vertices[i].Optim(0);
+	}/*}}}*/
 	/*FUNCTION Mesh::ProjectOnCurve{{{*/
 	GeomEdge*   Mesh::ProjectOnCurve( Edge & BhAB, BamgVertex &  vA, BamgVertex & vB,
@@ -3766,5 +3755,4 @@
 		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/
 
-		long int verbose=0;
 		long *renu= new long[nbt];
 		register Triangle *t0,*t,*te=triangles+nbt;
@@ -4889,5 +4877,4 @@
 	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ConsRefTriangle)*/
 
-	long int verbose=0;
 	register Triangle *t0,*t;
 	register long k=0, num;   
@@ -5416,5 +5403,5 @@
 
 				/*Initialize variables*/
-				double Lstep=0,Lcurve=0;    // step between two points   (phase==1) 
+				double Lstep=0;             // step between two points   (phase==1) 
 				long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1) 
 
@@ -5549,5 +5536,4 @@
 						long NbSegOnCurve = Max((long)(L+0.5),(long) 1);// nb of seg
 						Lstep = L/NbSegOnCurve; 
-						Lcurve = L;
 						NbCreatePointOnCurve = NbSegOnCurve-1;
 						NbOfNewEdge += NbSegOnCurve;
@@ -5680,6 +5666,6 @@
 				_error_("!v1 || !v2");
 			}
-			Icoor2 detss = 0,l=0,ks;
-			while ((ks=SwapForForcingEdge(  va,  vb, tc, detss, det1,det2,NbSwap)))
+			Icoor2 detss = 0,l=0;
+			while ((SwapForForcingEdge(  va,  vb, tc, detss, det1,det2,NbSwap)))
 			 if(l++ > 10000000) {
 				 _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big");
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16165)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16166)
@@ -157,13 +157,10 @@
 void Penta::BasalFrictionCreateInput(void){
 
-	/*Constants*/
-	const int  numdof=NUMVERTICES*NDOF1;
-
 	/*Intermediaries */
-	int    count;
-	IssmDouble basalfriction[NUMVERTICES]={0,0,0,0,0,0};
-	IssmDouble alpha2,vx,vy;
-	Friction*  friction=NULL;
-	GaussPenta* gauss=NULL;
+	int         count;
+	IssmDouble  basalfriction[NUMVERTICES];
+	IssmDouble  alpha2                       ,vx,vy;
+	Friction   *friction                   = NULL;
+	GaussPenta *gauss                      = NULL;
 
 	/* Basal friction can only be found at the base of an ice sheet: */
@@ -567,12 +564,10 @@
 
 	/*retrieve parameters: */
-	ElementMatrix* Ke=NULL;
 	ElementVector* De=NULL;
 	int analysis_type;
 	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 
-	/*Checks in debugging {{{*/
+	/*Checks in debugging*/
 	_assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
-	/*}}}*/
 
 	switch(analysis_type){
@@ -2144,11 +2139,13 @@
 void  Penta::InputToResult(int enum_type,int step,IssmDouble time){
 
-	bool   found = false;
 	Input *input = NULL;
 
 	/*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
-	if (enum_type==MaterialsRheologyBbarEnum) input=this->material->inputs->GetInput(MaterialsRheologyBEnum);
-	else if (enum_type==MaterialsRheologyZbarEnum) input=this->material->inputs->GetInput(MaterialsRheologyZEnum);
-	else input=this->inputs->GetInput(enum_type);
+	if      (enum_type==MaterialsRheologyBbarEnum)
+	 input=this->material->inputs->GetInput(MaterialsRheologyBEnum);
+	else if (enum_type==MaterialsRheologyZbarEnum)
+	 input=this->material->inputs->GetInput(MaterialsRheologyZEnum);
+	else
+	 input=this->inputs->GetInput(enum_type);
 	//if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found in penta->inputs"); why error out? if the requested input does not exist, we should still 
 	//try and output whatever we can instead of just failing.
@@ -2852,10 +2849,8 @@
 
 	/*Intermediaries*/
-	int     i;
-	int     numberofresults = 0;
-	int     *resultsenums   = NULL;
-	int     *resultssizes   = NULL;
-	IssmDouble  *resultstimes   = NULL;
-	int     *resultssteps   = NULL;
+	int        *resultsenums    = NULL;
+	int        *resultssizes    = NULL;
+	IssmDouble *resultstimes    = NULL;
+	int        *resultssteps    = NULL;
 
 	/*Checks*/
@@ -2863,11 +2858,7 @@
 
 	/*Count number of results*/
-	for(i=0;i<this->results->Size();i++){
-		ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
-		numberofresults++;
-	}
+	int numberofresults = this->results->Size();
 
 	if(numberofresults){
-
 		/*Allocate output*/
 		resultsenums=xNew<int>(numberofresults);
@@ -2877,5 +2868,5 @@
 
 		/*populate enums*/
-		for(i=0;i<this->results->Size();i++){
+		for(int i=0;i<this->results->Size();i++){
 			ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
 			resultsenums[i]=elementresult->InstanceEnum();
@@ -3781,7 +3772,4 @@
 void Penta::ViscousHeatingCreateInput(void){
 
-	/*Constants*/
-	const int  numdof=NUMVERTICES*NDOF1;
-
 	/*Intermediaries*/
 	IssmDouble phi;
@@ -3789,10 +3777,6 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble epsilon[6];
-	IssmDouble viscousheating[NUMVERTICES]={0,0,0,0,0,0};
+	IssmDouble viscousheating[NUMVERTICES];
 	IssmDouble thickness;
-	GaussPenta *gauss=NULL;
-
-	/*Initialize Element vector*/
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
 
 	/*Retrieve all inputs and parameters*/
@@ -3804,5 +3788,5 @@
 
 	/*loop over vertices: */
-	gauss=new GaussPenta();
+	GaussPenta* gauss=new GaussPenta();
 	for (int iv=0;iv<NUMVERTICES;iv++){
 		gauss->GaussVertex(iv);
@@ -4101,5 +4085,4 @@
 	IssmDouble D[3][3];
 	IssmDouble K[3][3]={0.0};
-	Tria*      tria=NULL;
 	GaussPenta *gauss=NULL;
 
@@ -4334,5 +4317,4 @@
 	IssmDouble D[3][3];
 	IssmDouble K[3][3]={0.0};
-	Tria*      tria=NULL;
 	GaussPenta *gauss=NULL;
 
Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 16165)
+++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 16166)
@@ -441,5 +441,4 @@
 
 	//   The penalty is stable if it doesn't change during to successive iterations.   
-	const int  numnodes=1;
 	IssmDouble pressure;
 	IssmDouble temperature;
@@ -623,5 +622,4 @@
 
 	//   The penalty is stable if it doesn't change during two consecutive iterations.   
-	const int  numnodes        = 1;
 	int        unstable        = 0;
 	int        new_active;
Index: /issm/trunk-jpl/src/c/classes/kriging/Quadtree.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/kriging/Quadtree.cpp	(revision 16165)
+++ /issm/trunk-jpl/src/c/classes/kriging/Quadtree.cpp	(revision 16166)
@@ -205,5 +205,5 @@
 	QuadtreeBox  *box  = NULL;
 	int           xi,yi;
-	int           level,levelbin;
+	int           levelbin;
 	int           index;
 	double        length,length2;
@@ -212,6 +212,5 @@
 	this->IntergerCoordinates(&xi,&yi,x,y);
 
-	/*Initialize levels*/
-	level    = 0;
+	/*Initialize level*/
 	levelbin = (1L<<this->MaxDepth);// = 2^30
 
@@ -221,7 +220,5 @@
 	/*Find the smallest box where this point is located*/
 	while((box=*pbox) && (box->nbitems<0)){ 
-
-		levelbin>>=1; level+=1; 
-
+		levelbin>>=1;
 		pbox = &box->box[IJ(xi,yi,levelbin)];
 	}
Index: /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.cpp	(revision 16165)
+++ /issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshx.cpp	(revision 16166)
@@ -13,12 +13,5 @@
 
 	/*Contour:*/
-	double   value;
-
-	/*threading: */
-	ContourToMeshxThreadStruct gate;
-	int num=1;
-	#ifdef _MULTITHREADING_
-	num=_NUMTHREADS_;
-	#endif
+	double value;
 
 	/*output: */
@@ -29,13 +22,14 @@
 
 	/*initialize thread parameters: */
-	gate.contours=contours;
-	gate.nods=nods;
-	gate.edgevalue=edgevalue;
-	gate.in_nod=in_nod;
-	gate.x=x;
-	gate.y=y;
+	ContourToMeshxThreadStruct gate;
+	gate.contours  = contours;
+	gate.nods      = nods;
+	gate.edgevalue = edgevalue;
+	gate.in_nod    = in_nod;
+	gate.x         = x;
+	gate.y         = y;
 
 	/*launch the thread manager with ContourToMeshxt as a core: */
-	LaunchThread(ContourToMeshxt,(void*)&gate,num);
+	LaunchThread(ContourToMeshxt,(void*)&gate,_NUMTHREADS_);
 
 	/*Take care of the case where an element interpolation has been requested: */
@@ -43,5 +37,5 @@
 		for(int n=0;n<nel;n++){
 			if ( (in_nod[ (int)*(index+3*n+0) -1] == 1) && (in_nod[ (int)*(index+3*n+1) -1] == 1) && (in_nod[ (int)*(index+3*n+2) -1] == 1) ){
-				value=1; in_elem[n]=value;
+				value=1.; in_elem[n]=value;
 			}
 		}
Index: /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp	(revision 16165)
+++ /issm/trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp	(revision 16166)
@@ -26,11 +26,4 @@
 	int     i;
 
-	/*threading: */
-	InterpFromGridToMeshxThreadStruct gate;
-	int num=1;
-	#ifdef _MULTITHREADING_
-	num=_NUMTHREADS_;
-	#endif
-
 	/*Some checks on arguments: */
 	if ((M<2) || (N<2) || (nods<=0)){
@@ -53,6 +46,6 @@
 		x=xNew<double>(N);
 		y=xNew<double>(M);
-		for (i=0;i<N;i++) x[i]=(x_in[i]+x_in[i+1])/2;
-		for (i=0;i<M;i++) y[i]=(y_in[i]+y_in[i+1])/2;
+		for(i=0;i<N;i++) x[i]=(x_in[i]+x_in[i+1])/2.;
+		for(i=0;i<M;i++) y[i]=(y_in[i]+y_in[i+1])/2.;
 		x_rows=x_rows-1;
 		y_rows=y_rows-1;
@@ -63,6 +56,6 @@
 		x=xNew<double>(N);
 		y=xNew<double>(M);
-		for (i=0;i<N;i++) x[i]=x_in[i];
-		for (i=0;i<M;i++) y[i]=y_in[i];
+		for(i=0;i<N;i++) x[i]=x_in[i];
+		for(i=0;i<M;i++) y[i]=y_in[i];
 	}
 	else{
@@ -71,20 +64,21 @@
 
 	/*initialize thread parameters: */
-	gate.x_mesh=x_mesh;
-	gate.y_mesh=y_mesh;
-	gate.x_rows=x_rows;
-	gate.y_rows=y_rows;
-	gate.x=x;
-	gate.y=y;
-	gate.nods=nods;
-	gate.data_mesh=data_mesh;
-	gate.data=data;
-	gate.default_value=default_value;
-	gate.interp=interpenum;
-	gate.M=M;
-	gate.N=N;
+	InterpFromGridToMeshxThreadStruct gate;
+	gate.x_mesh        = x_mesh;
+	gate.y_mesh        = y_mesh;
+	gate.x_rows        = x_rows;
+	gate.y_rows        = y_rows;
+	gate.x             = x;
+	gate.y             = y;
+	gate.nods          = nods;
+	gate.data_mesh     = data_mesh;
+	gate.data          = data;
+	gate.default_value = default_value;
+	gate.interp        = interpenum;
+	gate.M             = M;
+	gate.N             = N;
 
 	/*launch the thread manager with InterpFromGridToMeshxt as a core: */
-	LaunchThread(InterpFromGridToMeshxt,(void*)&gate,num);
+	LaunchThread(InterpFromGridToMeshxt,(void*)&gate,_NUMTHREADS_);
 	_printf_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%  \n");
 
Index: /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp	(revision 16165)
+++ /issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp	(revision 16166)
@@ -25,13 +25,5 @@
 
 	/*contours: */
-	double              *incontour     = NULL;
-
-	/*threading: */
-	InterpFromMesh2dxThreadStruct gate;
-	int num=1;
-
-	#ifdef _MULTITHREADING_
-	num=_NUMTHREADS_;
-	#endif
+	double *incontour     = NULL;
 
 	/*some checks*/
@@ -84,25 +76,26 @@
 
 	/*initialize thread parameters: */
-	gate.interpolation_type=interpolation_type;
-	gate.debug=debug;
-	gate.nels_data=nels_data;
-	gate.index_data=index_data;
-	gate.x_data=x_data;
-	gate.y_data=y_data;
-	gate.data=data;
-	gate.xmin=xmin;
-	gate.xmax=xmax;
-	gate.ymin=ymin;
-	gate.ymax=ymax;
-	gate.nods_prime=nods_prime;
-	gate.data_prime=data_prime;
-	gate.x_prime=x_prime;
-	gate.y_prime=y_prime;
-	gate.default_values=default_values;
-	gate.num_default_values=num_default_values;
-	gate.incontour=incontour;
+	InterpFromMesh2dxThreadStruct gate;
+	gate.interpolation_type = interpolation_type;
+	gate.debug              = debug;
+	gate.nels_data          = nels_data;
+	gate.index_data         = index_data;
+	gate.x_data             = x_data;
+	gate.y_data             = y_data;
+	gate.data               = data;
+	gate.xmin               = xmin;
+	gate.xmax               = xmax;
+	gate.ymin               = ymin;
+	gate.ymax               = ymax;
+	gate.nods_prime         = nods_prime;
+	gate.data_prime         = data_prime;
+	gate.x_prime            = x_prime;
+	gate.y_prime            = y_prime;
+	gate.default_values     = default_values;
+	gate.num_default_values = num_default_values;
+	gate.incontour          = incontour;
 
 	/*launch the thread manager with InterpFromGridToMeshxt as a core: */
-	LaunchThread(InterpFromMesh2dxt,(void*)&gate,num);
+	LaunchThread(InterpFromMesh2dxt,(void*)&gate,_NUMTHREADS_);
 
 	/*Assign output pointers:*/
Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 16165)
+++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 16166)
@@ -28,8 +28,5 @@
 	/*threading: */
 	KrigingxThreadStruct gate;
-	int num=1;
-#ifdef _MULTITHREADING_
-	num=_NUMTHREADS_;
-#endif
+	int num = _NUMTHREADS_;
 
 	/*Get Variogram from Options*/
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Stressbalance/CreateConstraintsStressbalance.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Stressbalance/CreateConstraintsStressbalance.cpp	(revision 16165)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Stressbalance/CreateConstraintsStressbalance.cpp	(revision 16166)
@@ -40,5 +40,4 @@
 	/*Output*/
 	Constraints *constraints      = NULL;
-	SpcStatic   *spcstatic        = NULL;
 
 	/*Fetch parameters: */
Index: /issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp	(revision 16165)
+++ /issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp	(revision 16166)
@@ -11,20 +11,14 @@
 
 	/*threading: */
-	PointCloudFindNeighborsThreadStruct gate;
-
-	#ifdef _MULTITHREADING_
 	int num=_NUMTHREADS_;
-	#else
-	int num=1;
-	#endif
-
 	if(!multithread)num=1;
 
 	/*initialize thread parameters: */
-	gate.x=x;
-	gate.y=y;
-	gate.nods=nods;
-	gate.mindistance=mindistance;
-	gate.flags=flags;
+	PointCloudFindNeighborsThreadStruct gate;
+	gate.x           = x;
+	gate.y           = y;
+	gate.nods        = nods;
+	gate.mindistance = mindistance;
+	gate.flags       = flags;
 
 	/*launch the thread manager with InterpFromGridToMeshxt as a core: */