Index: /issm/trunk/src/c/BuildNodeSetsx/PartitionSets.cpp
===================================================================
--- /issm/trunk/src/c/BuildNodeSetsx/PartitionSets.cpp	(revision 3594)
+++ /issm/trunk/src/c/BuildNodeSetsx/PartitionSets.cpp	(revision 3595)
@@ -13,6 +13,8 @@
 	/*output: */
 	Vec partitionb=NULL;
+	Vec partitionb_vec=NULL;
 	double* partitionb_local=NULL;
 	Vec partitionc=NULL;
+	Vec partitionc_vec=NULL;
 	double* partitionc_local=NULL;
 
@@ -84,6 +86,23 @@
 
 	/*Now, using the local partitions for b and c, create parallel vectors: */
-	VecCreateMPIWithArray(MPI_COMM_WORLD,bsize,PETSC_DECIDE,partitionb_local,&partitionb);
-	VecCreateMPIWithArray(MPI_COMM_WORLD,csize,PETSC_DECIDE,partitionc_local,&partitionc);
+	VecCreateMPIWithArray(MPI_COMM_WORLD,bsize,PETSC_DECIDE,partitionb_local,&partitionb_vec);
+	VecCreateMPIWithArray(MPI_COMM_WORLD,csize,PETSC_DECIDE,partitionc_local,&partitionc_vec);
+
+	VecAssemblyBegin(partitionb_vec); 
+	VecAssemblyEnd(partitionb_vec);
+
+	VecAssemblyBegin(partitionc_vec); 
+	VecAssemblyEnd(partitionc_vec);
+
+	/*Now we must duplicate these vectors because:
+	 *
+	 * Petsc specifies on http://www.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Vec/VecCreateMPIWithArray.html
+	 * "PETSc does NOT free the array when the vector is destroyed via VecDestroy().
+	 * The user should not free the array until the vector is destroyed"
+	 *
+	 * The only way we can keep the vector partitionb_vec while destroying the array partitionb_local
+	 * is to duplicate the vector and then destroy the initial vector as well as the array*/
+	VecDuplicate(partitionb_vec,&partitionb); VecCopy(partitionb_vec,partitionb); VecFree(&partitionb_vec); xfree((void**)&partitionb_local);
+	VecDuplicate(partitionc_vec,&partitionc); VecCopy(partitionc_vec,partitionc); VecFree(&partitionc_vec); xfree((void**)&partitionc_local);
 
 	/*Free ressources:*/
@@ -93,10 +112,4 @@
 	xfree((void**)&flags_c_local);
 	
-	VecAssemblyBegin(partitionb); 
-	VecAssemblyEnd(partitionb);
-
-	VecAssemblyBegin(partitionc); 
-	VecAssemblyEnd(partitionc);
-
 	/*Assign output pointers*/
 	*ppartitionb=partitionb;
Index: /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp
===================================================================
--- /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp	(revision 3594)
+++ /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp	(revision 3595)
@@ -18,11 +18,11 @@
 	extern int my_rank;
 	
-	_printf_("      Configuring elements...\n");
+	//_printf_("      Configuring elements...\n");
 	elements->Configure(elements,loads,nodes,vertices,materials,parameters);
-	_printf_("      Configuring loads...\n");
+	//_printf_("      Configuring loads...\n");
 	loads->Configure(elements,loads,nodes,vertices,materials,parameters);
-	_printf_("      Configuring nodes...\n");
+	//_printf_("      Configuring nodes...\n");
 	nodes->Configure(elements,loads,nodes,vertices,materials,parameters);
-	_printf_("      Configuring parameters...\n");
+	//_printf_("      Configuring parameters...\n");
 	parameters->Configure(elements,loads, nodes,vertices, materials,parameters);
 
Index: /issm/trunk/src/c/ContourToMeshx/ContourToMeshxt.cpp
===================================================================
--- /issm/trunk/src/c/ContourToMeshx/ContourToMeshxt.cpp	(revision 3594)
+++ /issm/trunk/src/c/ContourToMeshx/ContourToMeshxt.cpp	(revision 3595)
@@ -62,7 +62,4 @@
 	/*Loop through all contours: */
 	for (i=0;i<numcontours;i++){
-		#ifdef _ISSM_DEBUG_
-			printf("\nHandling contour %i/%i\n",i,numcontours);
-		#endif
 		contouri=*(contours+i);
 		numgrids=contouri->nods;
Index: /issm/trunk/src/c/ContourToNodesx/ContourToNodesx.cpp
===================================================================
--- /issm/trunk/src/c/ContourToNodesx/ContourToNodesx.cpp	(revision 3594)
+++ /issm/trunk/src/c/ContourToNodesx/ContourToNodesx.cpp	(revision 3595)
@@ -23,7 +23,4 @@
 	/*Loop through all contours: */
 	for (i=0;i<numcontours;i++){
-		#ifdef _ISSM_DEBUG_
-			printf("\nHandling contour %i/%i\n",i,numcontours);
-		#endif
 		contouri=*(contours+i);
 		numgrids=contouri->nods;
Index: /issm/trunk/src/c/DataSet/DataSet.cpp
===================================================================
--- /issm/trunk/src/c/DataSet/DataSet.cpp	(revision 3594)
+++ /issm/trunk/src/c/DataSet/DataSet.cpp	(revision 3595)
@@ -179,8 +179,4 @@
 	}
 
-#ifdef _ISSM_DEBUG_
-	_printf_("Number of objects in dataset being demarshalled: %i\n",numobjects);
-#endif
-
 	for(i=0;i<numobjects;i++){
 
@@ -1152,10 +1148,4 @@
 	Ranks(ranks);
 
-	#ifdef _ISSM_DEBUG_
-	for(i=0;i<numberofobjects;i++){
-		_printf_("%i\n",ranks[i]);
-	}
-	#endif
-
 	/*We need to take the minimum rank for each vertex, and every cpu needs to get that result. That way, 
 	 * when we start building the dof list for all vertexs, a cpu can check whether its vertex already has been 
@@ -1163,10 +1153,4 @@
 	 * order of cpu rank. This is also why we initialized this array to num_procs.*/
 	MPI_Allreduce ( (void*)ranks,(void*)minranks,numberofobjects,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
-
-	#ifdef _ISSM_DEBUG_
-		for(i=0;i<numberofobjects;i++){
-			_printf_("%i\n",minranks[i]);
-		}
-	#endif
 
 	/*Now go through all objects, and use minranks to flag which objects are cloned: */
Index: /issm/trunk/src/c/HoleFillerx/HoleFillerx.cpp
===================================================================
--- /issm/trunk/src/c/HoleFillerx/HoleFillerx.cpp	(revision 3594)
+++ /issm/trunk/src/c/HoleFillerx/HoleFillerx.cpp	(revision 3595)
@@ -43,13 +43,4 @@
 	int             imageoutsize;
 
-	#ifdef _ISSM_DEBUG_
-		if ( smooth == 1 ){
-			printf("Data patches will be SMOOTHED. \n");
-		}
-		if ( smooth == 0 ){ 
-			printf("Data patches will NOT be smoothed. \n");
-		}
-	#endif
-	
 	 /*^^^^^^^^^^^^^  Remove pixels close to the holes ^^^^^^^^^^^^^*/
 	image2 = (double*) xmalloc( lines*samps*sizeof(double));
@@ -92,9 +83,4 @@
 	#endif
 
-	#ifdef _ISSM_DEBUG_
-		printf( "\n" );
-		printf( "Iterations %5ld   and %5ld.", count-1, count );
-		fflush( stdout );
-	#endif
 	goto afterfirst2;
 
@@ -115,11 +101,4 @@
 		printf( "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b" );
 		printf("Number of zeroes remaining: %10ld",lines*samps-counter);
-		fflush( stdout );
-	#endif
-
-	#ifdef _ISSM_DEBUG_
-		printf( "\b\b\b\b\b\b\b\b\b\b\b\b\b\b" );
-		printf( "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b" );
-		printf( "Iterations %5ld   and %5ld.", count-1, count );
 		fflush( stdout );
 	#endif
@@ -366,17 +345,6 @@
 		}
 
-		#ifdef _ISSM_DEBUG_
-			if ( !(i%100) ){
-			fprintf(stdout, "\b\b\b\b\b\b%5ld ", i );   /* Count lines on output screen */
-			fflush(stdout);	
-			}
-		#endif
-	}
-
-	#ifdef _ISSM_DEBUG_
-		fprintf(stdout, "\b\b\b\b\b\b%5ld ", lines );
-		fflush(stdout);	
-	#endif
-	
+	}
+
 		
 there2:		
@@ -388,9 +356,4 @@
 
 	time(&t2);	
-	#ifdef _ISSM_DEBUG_
-		printf ( "\n\nEnd                  "); printf( ctime(&t2) );
-		printf( "\n                           Total time:        ");
-		printf( "%4.2f minutes. \n", difftime(t2,t1)/60 );		
-	#endif
 	
 	#ifdef _DEBUG2_
@@ -402,8 +365,4 @@
 	end:
 
-	#ifdef _ISSM_DEBUG_
-		printf( "Done.\a\n" );
-	#endif
-
 	/*Assign output pointers: */
 	*pimageout=imageout;
Index: /issm/trunk/src/c/MeshPartitionx/MeshPartitionx.cpp
===================================================================
--- /issm/trunk/src/c/MeshPartitionx/MeshPartitionx.cpp	(revision 3594)
+++ /issm/trunk/src/c/MeshPartitionx/MeshPartitionx.cpp	(revision 3595)
@@ -81,12 +81,4 @@
 		}
 	}
-
-	#ifdef _ISSM_DEBUG_
-	if(my_rank==0){
-		for (i=0;i<numberofelements;i++){
-			printf("El %i My rank %i\n",i+1,epart[i]);
-		}
-	}
-	#endif
 	
 	/*Assign output pointer:*/
Index: /issm/trunk/src/c/ModelProcessorx/Partitioning.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Partitioning.cpp	(revision 3594)
+++ /issm/trunk/src/c/ModelProcessorx/Partitioning.cpp	(revision 3595)
@@ -150,4 +150,5 @@
 	xfree((void**)&npart);
 	xfree((void**)&epart);
+	xfree((void**)&serial_bordervertices);
 	VecFree(&bordervertices);
 
Index: /issm/trunk/src/c/ModelProcessorx/Qmu/CreateParametersQmu.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Qmu/CreateParametersQmu.cpp	(revision 3594)
+++ /issm/trunk/src/c/ModelProcessorx/Qmu/CreateParametersQmu.cpp	(revision 3595)
@@ -155,15 +155,4 @@
 		param->SetStringArray(responsedescriptors,iomodel->numberofresponses);
 		parameters->AddObject(param);
-
-		#ifdef _ISSM_DEBUG_
-			for(i=0;i<iomodel->numberofvariables;i++){
-				_printf_("variable descriptor %s\n",variabledescriptors[i]);
-			}
-			
-			for(i=0;i<iomodel->numberofresponses;i++){
-				_printf_("response descriptor %s\n",responsedescriptors[i]);
-			}
-		#endif
-
 
 		/*partition vertices in iomodel->qmu_npart parts, unless a partition is already present: */
Index: /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp
===================================================================
--- /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp	(revision 3594)
+++ /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp	(revision 3595)
@@ -307,8 +307,4 @@
 	#endif
 
-	#ifdef _ISSM_DEBUG_
-		printf("Maximum rift penetration2: %g\n",max_penetration);
-	#endif
-
 	/*feed max_penetration to inputs: */
 	inputs->Add("max_penetration",max_penetration);
Index: /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp
===================================================================
--- /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp	(revision 3594)
+++ /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp	(revision 3595)
@@ -27,7 +27,4 @@
 	/*Now, figure out maximum value of K_gg, so that we can penalize it correctly: */
 	MatNorm(Kgg,NORM_INFINITY,&kmax);
-	#ifdef _ISSM_DEBUG_
-	 _printf_("   K_gg infinity norm: %g\n",kmax);
-	#endif
 
 	/*Add penalties to stiffnesses, from loads: */
@@ -46,9 +43,4 @@
 	}
 
-	#ifdef _ISSM_DEBUG_
-	MatNorm(Kgg,NORM_INFINITY,&kmax2);
-	 _printf_("   K_gg infinity norm after penalties: %g\n",kmax2);
-	#endif
-
 	/*Assign output pointers:*/
 	if(pkmax)*pkmax=kmax;
Index: /issm/trunk/src/c/Qmux/SpawnCoreParallel.cpp
===================================================================
--- /issm/trunk/src/c/Qmux/SpawnCoreParallel.cpp	(revision 3594)
+++ /issm/trunk/src/c/Qmux/SpawnCoreParallel.cpp	(revision 3595)
@@ -68,10 +68,4 @@
 	model->FindParam(&qmu_npart,"qmu_npart");
 	model->FindParam(&qmu_part,NULL,NULL,"qmu_part");
-	#ifdef _ISSM_DEBUG_
-	for(i=0;i<numresponses;i++){
-		PetscSynchronizedPrintf(MPI_COMM_WORLD,"response descriptor %i: %s\n",i,responses_descriptors[i]);
-		PetscSynchronizedFlush(MPI_COMM_WORLD);
-	}
-	#endif
 
 	/*broadcast variables: only cpu 0 has correct values*/
@@ -79,11 +73,4 @@
 	if(my_rank!=0)variables=(double*)xmalloc(numvariables*sizeof(double));
 	MPI_Bcast(variables,numvariables,MPI_DOUBLE,0,MPI_COMM_WORLD); 
-
-	#ifdef _ISSM_DEBUG_
-	for(i=0;i<numvariables;i++){
-		PetscSynchronizedPrintf(MPI_COMM_WORLD,"variable %i: %g\n",i,variables[i]);
-		PetscSynchronizedFlush(MPI_COMM_WORLD);
-	}
-	#endif
 
 	/*broadcast variables_descriptors: */
@@ -102,21 +89,6 @@
 	}
 
-	#ifdef _ISSM_DEBUG_
-	for(i=0;i<numvariables;i++){
-		PetscSynchronizedPrintf(MPI_COMM_WORLD,"variable descriptor %i: %s value: %g\n",i,variables_descriptors[i],variables[i]);
-		PetscSynchronizedFlush(MPI_COMM_WORLD);
-	}
-	#endif
-
 	/*broadcast numresponses: */
 	MPI_Bcast(&numresponses,1,MPI_INT,0,MPI_COMM_WORLD); 
-
-	#ifdef _ISSM_DEBUG_
-	for(i=0;i<numresponses;i++){
-		PetscSynchronizedPrintf(MPI_COMM_WORLD,"variable descriptor %i: %s value: %g\n",i,responses_descriptors[i],responses[i]);
-		PetscSynchronizedFlush(MPI_COMM_WORLD);
-	}
-	#endif
-
 
 	_printf_("qmu iteration: %i\n",counter);
Index: /issm/trunk/src/c/Solverx/Solverx.cpp
===================================================================
--- /issm/trunk/src/c/Solverx/Solverx.cpp	(revision 3594)
+++ /issm/trunk/src/c/Solverx/Solverx.cpp	(revision 3595)
@@ -89,11 +89,4 @@
 		}
 	}
-		
-
-	#ifdef _ISSM_DEBUG_ 
-		KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
-		KSPGetInitialGuessNonzero(ksp,&flag); 
-		_printf_("User provided initial guess solution: %i\n",flag);
-	#endif
 
 	KSPSolve(ksp,pf,uf);
@@ -103,9 +96,4 @@
 	if (iteration_number<0){
 		ISSMERROR("%s%i"," Solver diverged at iteration number: ",-iteration_number);
-	}
-	else{
-		#ifdef _ISSM_DEBUG_ 
-		_printf_("%s%i%s\n","Solver converged after ",iteration_number," iterations");
-		#endif
 	}
 
Index: /issm/trunk/src/c/include/macros.h
===================================================================
--- /issm/trunk/src/c/include/macros.h	(revision 3594)
+++ /issm/trunk/src/c/include/macros.h	(revision 3595)
@@ -18,7 +18,12 @@
 #endif
 
-/*Assertion macro*/
+/*Assertion macro: do nothing if macro _ISSM_DEBUG_ undefined*/
+#ifdef _ISSM_DEBUG_ 
 #define ISSMASSERT(statement)\
   if (!(statement)) ISSMERROR("Assertion \"%s\" failed, please report bug to an ISSM developer",#statement)
+#else
+#define ISSMASSERT(ignore)\
+  ((void) 0)
+#endif
 
 /*The following macros hide the error exception handling in a matlab module. Just put 
Index: /issm/trunk/src/c/objects/Matice.cpp
===================================================================
--- /issm/trunk/src/c/objects/Matice.cpp	(revision 3594)
+++ /issm/trunk/src/c/objects/Matice.cpp	(revision 3595)
@@ -250,7 +250,4 @@
 		}
 	}
-	#ifdef _ISSM_DEBUG_
-	printf("Viscosity %lf\n",viscosity);
-	#endif
 
 	/*Return: */
@@ -314,7 +311,4 @@
 		}
 	}
-	#ifdef _ISSM_DEBUG_
-	printf("Viscosity %lf\n",viscosity3d);
-	#endif
 
 	/*Assign output pointers:*/
@@ -382,8 +376,4 @@
 	}
 
-	#ifdef _ISSM_DEBUG_
-	printf("Viscosity %lf\n",viscosity3d);
-	#endif
-
 	/*Assign output pointers:*/
 	*pviscosity3d=viscosity3d;
@@ -436,8 +426,4 @@
 	}
 		
-	#ifdef _ISSM_DEBUG_
-	printf("viscosity_complement %lf\n",viscosity_complement);
-	#endif
-
 	/*Return: */
 	*pviscosity_complement=viscosity_complement;
Index: /issm/trunk/src/c/objects/Numericalflux.cpp
===================================================================
--- /issm/trunk/src/c/objects/Numericalflux.cpp	(revision 3594)
+++ /issm/trunk/src/c/objects/Numericalflux.cpp	(revision 3595)
@@ -256,7 +256,5 @@
 		GetParameterValue(&vy, &vy_list[0],gauss_coord);
 		UdotN=vx*normal[0]+vy*normal[1];
-		if (fabs(UdotN)<1.0e-9) printf("Edge number %i has a flux very small (u.n = %g ), which could lead to unaccurate results\n",id,UdotN);
-		//if (fabs(UdotN)>0 && fabs(UdotN)< 1.0e-7) UdotN= 1.0e-7;
-		//if (fabs(UdotN)<0 && fabs(UdotN)>-1.0e-7) UdotN=-1.0e-7;
+		if (fabs(UdotN)<1.0e-9 && analysis_type==Balancedthickness2AnalysisEnum) printf("Edge number %i has a flux very small (u.n = %g ), which could lead to unaccurate results\n",id,UdotN);
 
 		/*Get L and B: */
Index: /issm/trunk/src/c/objects/Numpar.cpp
===================================================================
--- /issm/trunk/src/c/objects/Numpar.cpp	(revision 3594)
+++ /issm/trunk/src/c/objects/Numpar.cpp	(revision 3595)
@@ -45,4 +45,8 @@
 /*FUNCTION Numpar::destructor {{{1*/
 Numpar::~Numpar(){
+
+	/*Free the only pointer*/
+	xfree((void**)&control_type);
+
 	return;
 }
Index: /issm/trunk/src/c/objects/ParameterInputs.cpp
===================================================================
--- /issm/trunk/src/c/objects/ParameterInputs.cpp	(revision 3594)
+++ /issm/trunk/src/c/objects/ParameterInputs.cpp	(revision 3595)
@@ -361,14 +361,4 @@
 			}
 
-			#ifdef _ISSM_DEBUG_
-				PetscSynchronizedPrintf(MPI_COMM_WORLD,"Parameter vetor:");
-				PetscSynchronizedFlush(MPI_COMM_WORLD);
-				for(k=0;k<numberofnodes;k++){
-					PetscSynchronizedPrintf(MPI_COMM_WORLD," node %i value %g\n",k+1,parameter[k]);
-					PetscSynchronizedFlush(MPI_COMM_WORLD);
-				}
-			#endif
-			  
-
 			/*Add parameter to inputs: */
 			this->Add(root,parameter,1,numberofnodes);
Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 3594)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 3595)
@@ -1247,12 +1247,4 @@
 
 	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 _ISSM_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));
-	}
-	for (i=0;i<num_vert_gauss;i++){
-		printf("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
-	}
-#endif
 
 	/* Start  looping on the number of gaussian points: */
@@ -1773,12 +1765,4 @@
 
 		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 _ISSM_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));
-		}
-		for (i=0;i<num_vert_gauss;i++){
-			printf("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
-		}
-#endif
 
 		/* Start  looping on the number of gaussian points: */
@@ -2177,12 +2161,4 @@
 
 	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 _ISSM_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));
-	}
-	for (i=0;i<num_vert_gauss;i++){
-		printf("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
-	}
-#endif
 
 	/* Start  looping on the number of gaussian points: */
@@ -2209,7 +2185,4 @@
 			/* Get Jacobian determinant: */
 			GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
-#ifdef _ISSM_DEBUG_ 
-			printf("Element id %i Jacobian determinant: %lf\n",GetId(),Jdet);
-#endif
 
 			/*Get nodal functions: */
@@ -2697,10 +2670,4 @@
 	GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
 
-#ifdef _ISSM_DEBUG_ 
-	for (i=0;i<numgrids;i++){
-		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh6[0][i],dh1dh6[1][i],dh1dh6[2][i]);
-	}
-#endif
-
 	/*Build B: */
 	for (i=0;i<numgrids;i++){
@@ -2834,10 +2801,4 @@
 	/*Get dh1dh6 in actual coordinate system: */
 	GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
-
-#ifdef _ISSM_DEBUG_ 
-	for (i=0;i<numgrids;i++){
-		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh6[0][i],dh1dh6[1][i],dh1dh6[2][i]);
-	}
-#endif
 
 	/*Build B: */
@@ -2881,10 +2842,4 @@
 	/*Get dh1dh6 in actual coordinate system: */
 	GetNodalFunctionsDerivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
-
-#ifdef _ISSM_DEBUG_ 
-	for (i=0;i<numgrids;i++){
-		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh6[0][i],dh1dh6[1][i],dh1dh6[2][i]);
-	}
-#endif
 
 	/*Build BPrime: */
@@ -2983,11 +2938,4 @@
 
 	GetNodalFunctions(l1l6, gauss_coord);
-
-#ifdef _ISSM_DEBUG_ 
-	for (i=0;i<6;i++){
-		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7[0][i],dh1dh7[1][i]);
-	}
-
-#endif
 
 	/*B_primeuild B_prime: */
@@ -3063,11 +3011,4 @@
 	GetNodalFunctions(l1l6, gauss_coord);
 
-#ifdef _ISSM_DEBUG_ 
-	for (i=0;i<7;i++){
-		printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7[0][i],dh1dh7[1][i],dh1dh7[2][i]);
-	}
-
-#endif
-
 	/*Build B: */
 	for (i=0;i<numgrids+1;i++){
@@ -3216,12 +3157,4 @@
 	*(J+NDOF3*2+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
 
-#ifdef _ISSM_DEBUG_
-	for(i=0;i<3;i++){
-		for (j=0;j<3;j++){
-			printf("%lf ",*(J+NDOF3*i+j));
-		}
-		printf("\n");
-	}
-#endif
 }
 /*}}}*/
@@ -3841,43 +3774,4 @@
 	GetB(&B[0][0], xyz_list, gauss_coord);
 
-#ifdef _ISSM_DEBUG_
-	printf("B for grid1 : [ %lf   %lf  \n",B[0][0],B[0][1]);
-	printf("              [ %lf   %lf  \n",B[1][0],B[1][1]);
-	printf("              [ %lf   %lf  ]\n",B[2][0],B[2][1]);
-	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]);
-	printf("              [ %lf   %lf  ]\n",B[2][2],B[2][3]);
-	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]);
-	printf("              [ %lf   %lf  ]\n",B[2][4],B[2][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]);
-	printf("              [ %lf   %lf  ]\n",B[2][6],B[2][7]);
-	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]);
-	printf("              [ %lf   %lf  ]\n",B[2][8],B[2][9]);
-	printf("              [ %lf   %lf  ]\n",B[3][8],B[3][9]);
-	printf("              [ %lf   %lf  ]\n",B[4][8],B[4][9]);
-
-	printf("B for grid6 : [ %lf   %lf  \n", B[0][10],B[0][11]);
-	printf("              [ %lf   %lf  \n", B[1][10],B[1][11]);
-	printf("              [ %lf   %lf  ]\n",B[2][10],B[2][11]);
-	printf("              [ %lf   %lf  ]\n",B[3][10],B[3][11]);
-	printf("              [ %lf   %lf  ]\n",B[4][10],B[4][11]);
-
-#endif
-
 	/*Multiply B by velocity, to get strain rate: */
 	MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
Index: /issm/trunk/src/c/objects/Riftfront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Riftfront.cpp	(revision 3594)
+++ /issm/trunk/src/c/objects/Riftfront.cpp	(revision 3595)
@@ -604,16 +604,8 @@
 		 *contact slip friction. */
 		  
-		#ifdef _ISSM_DEBUG_
-		printf("Dealing with grid pair (%i,%i)\n",nodes[0]->GetId(),nodes[1]->GetId());
-		#endif
-
 		/*Recover input parameters: */
 		inputs->Recover("thickness",&h[0],1,dofs,MAX_RIFTFRONT_GRIDS,(void**)nodes);
 		if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts");
 		thickness=h[0];
-
-		#ifdef _ISSM_DEBUG_
-			printf("Thickness at grid (%i,%i): %lg\n",nodes[0]->GetId(),nodes[1]->GetID(),thickness);
-		#endif
 
 		/*From Peter Wriggers book (Computational Contact Mechanics, p191): */
@@ -708,8 +700,4 @@
 		 * and we want to avoid zigzagging of the loads, we want lump the loads onto grids, not onto surfaces between grids.:*/
 	
-		#ifdef _ISSM_DEBUG_
-		_printf_("Grids  (%i,%i) are free of constraints\n",nodes[0]->GetId(),nodes[1]->GetID());
-		#endif
-
 		/*Ok, to compute the pressure, we are going to need material properties, thickness and bed for the two grids. We assume those properties to 
 		 * be the same across the rift.: */
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 3594)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 3595)
@@ -511,10 +511,4 @@
 	  /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
 	  GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
-
-#ifdef _ISSM_DEBUG_ 
-	for (i=0;i<num_gauss;i++){
-		printf("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),*(gauss_weights+i));
-	}
-#endif
 
 	/* Start  looping on the number of gaussian points: */
@@ -3804,7 +3798,4 @@
 		/* Get Jacobian determinant: */
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-#ifdef _ISSM_DEBUG_ 
-		printf("Element id %i Jacobian determinant: %g\n",GetId(),Jdet);
-#endif
 
 		/* Get nodal functions value at gaussian point:*/
@@ -3932,10 +3923,4 @@
 	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss_l1l2l3);
 
-	#ifdef _ISSM_DEBUG_ 
-	for (i=0;i<3;i++){
-		printf("Node %i  dh/dx=%lf dh/dy=%lf \n",i,dh1dh3[0][i],dh1dh3[1][i]);
-	}
-	#endif
-
 	/*Build B: */
 	for (i=0;i<numgrids;i++){
@@ -3972,10 +3957,4 @@
 	/*Get dh1dh2dh3 in actual coordinate system: */
 	GetNodalFunctions(&l1l2l3[0],gauss_l1l2l3);
-
-#ifdef _ISSM_DEBUG_ 
-	for (i=0;i<3;i++){
-		printf("Node %i  h=%lf \n",i,l1l2l3[i]);
-	}
-#endif
 
 	/*Build B_prog: */
@@ -4556,9 +4535,4 @@
 	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
 	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
-#ifdef _ISSM_DEBUG_ 
-	for (i=0;i<num_gauss;i++){
-		printf("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),*(gauss_weights+i));
-	}
-#endif
 
 	/* Start  looping on the number of gaussian points: */
@@ -4590,7 +4564,4 @@
 		/* Get nodal functions value at gaussian point:*/
 		GetNodalFunctions(l1l2l3, gauss_l1l2l3);
-#ifdef _ISSM_DEBUG_
-		printf("viscositycomp %g thickness %g dvx [%g %g] dvy [%g %g]  dadjx [%g %g] dadjy[%g %g]\n",viscosity_complement,thickness,dvx[0],dvx[1],dvy[0],dvy[1],dadjx[0],dadjx[1],dadjy[0],dadjy[1]);
-#endif
 
 		/*Get nodal functions derivatives*/
@@ -4797,28 +4768,16 @@
 		GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
 		GetParameterValue(&drag, &this->properties.k[0],gauss_l1l2l3);
-		#ifdef _ISSM_DEBUG_ 
-			printf("Drag complement: %20.20lf Drag: %20.20lf\n",alpha_complement,drag);
-		#endif
 
 		/*recover lambda and mu: */
 		GetParameterValue(&lambda, &adjx_list[0],gauss_l1l2l3);
 		GetParameterValue(&mu, &adjy_list[0],gauss_l1l2l3);
-		#ifdef _ISSM_DEBUG_ 
-			printf("Adjoint vector %20.20lf %20.20lf\n",lambda,mu);
-		#endif
 			
 		/*recover vx and vy: */
 		GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
 		GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
-		#ifdef _ISSM_DEBUG_ 
-			printf("Velocity vector %20.20lf %20.20lf\n",vx,vy);
-		#endif
 
 		/* Get Jacobian determinant: */
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-		#ifdef _ISSM_DEBUG_ 
-		printf("Element id %i Jacobian determinant: %lf\n",GetId(),Jdet);
-		#endif
-		
+	
 		/* Get nodal functions value at gaussian point:*/
 		GetNodalFunctions(l1l2l3, gauss_l1l2l3);
@@ -5031,7 +4990,4 @@
 		GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
 		GetParameterValue(&drag, &this->properties.k[0],gauss_l1l2l3);
-#ifdef _ISSM_DEBUG_ 
-		printf("Drag complement: %20.20lf Drag: %20.20lf\n",alpha_complement,drag);
-#endif
 
 		/*recover lambda mu and xi: */
@@ -5039,7 +4995,4 @@
 		GetParameterValue(&mu, &adjy_list[0],gauss_l1l2l3);
 		GetParameterValue(&xi, &adjz_list[0],gauss_l1l2l3);
-#ifdef _ISSM_DEBUG_ 
-		printf("Adjoint vector %20.20lf %20.20lf\n",lambda,mu);
-#endif
 
 		/*recover vx vy and vz: */
@@ -5047,20 +5000,7 @@
 		GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
 		GetParameterValue(&vz, &vz_list[0],gauss_l1l2l3);
-#ifdef _ISSM_DEBUG_ 
-		printf("Velocity vector %20.20lf %20.20lf\n",vx,vy);
-
-		/*Get normal vecyor to the bed */
-		SurfaceNormal(&surface_normal[0],xyz_list);
-
-		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];
-		bed_normal[2]=-surface_normal[2];
-#endif
 
 		/* Get Jacobian determinant: */
 		GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-#ifdef _ISSM_DEBUG_ 
-		printf("Element id %i Jacobian determinant: %lf\n",GetId(),Jdet);
-#endif
 
 		/* Get nodal functions value at gaussian point:*/
Index: /issm/trunk/src/c/parallel/balancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/parallel/balancedthickness.cpp	(revision 3594)
+++ /issm/trunk/src/c/parallel/balancedthickness.cpp	(revision 3595)
@@ -140,4 +140,5 @@
 	delete processedresults;
 	delete results;
+	delete model;
 	delete inputs;
 
Index: /issm/trunk/src/c/parallel/diagnostic.cpp
===================================================================
--- /issm/trunk/src/c/parallel/diagnostic.cpp	(revision 3594)
+++ /issm/trunk/src/c/parallel/diagnostic.cpp	(revision 3595)
@@ -173,4 +173,5 @@
 	xfree((void**)&u_g_obs);
 	xfree((void**)&weights);
+	xfree((void**)&control_type);
 	delete model;
 	delete inputs;
Index: /issm/trunk/src/c/parallel/prognostic.cpp
===================================================================
--- /issm/trunk/src/c/parallel/prognostic.cpp	(revision 3594)
+++ /issm/trunk/src/c/parallel/prognostic.cpp	(revision 3595)
@@ -25,5 +25,4 @@
 	Model* model=NULL;
 
-	Vec     u_g=NULL;
 	double* u_g_serial=NULL;
 	double* h_g_initial=NULL;
@@ -141,4 +140,8 @@
 
 	/*Free ressources:*/
+	xfree((void**)&u_g_serial);
+	xfree((void**)&h_g_initial);
+	xfree((void**)&melting_g);
+	xfree((void**)&accumulation_g);
 	delete processedresults;
 	delete results;
Index: /issm/trunk/src/c/shared/Exp/DomainOutlineRead.cpp
===================================================================
--- /issm/trunk/src/c/shared/Exp/DomainOutlineRead.cpp	(revision 3594)
+++ /issm/trunk/src/c/shared/Exp/DomainOutlineRead.cpp	(revision 3595)
@@ -59,7 +59,4 @@
 		nprof++;
 	}
-	#ifdef _ISSM_DEBUG_
-		printf("Number of profiles in domain outline file: %i\n",nprof);
-	#endif
 	
 	/*Allocate and initialize all the profiles: */
@@ -113,17 +110,4 @@
 	fclose(fid);
 
-	
-	#ifdef _ISSM_DEBUG_
-	for (i=0;i<nprof;i++){
-		printf("Profile #%i\n",i);
-		x=pprofx[i];
-		y=pprofy[i];
-		for (j=0;j<profngrids[i];j++){
-			printf("   %lf %lf\n",x[j],y[j]);
-		}
-	}
-	#endif
-
-
 	cleanupandreturn: 
 	/*Free ressources: */
Index: /issm/trunk/src/c/shared/Sorting/binary_search.cpp
===================================================================
--- /issm/trunk/src/c/shared/Sorting/binary_search.cpp	(revision 3594)
+++ /issm/trunk/src/c/shared/Sorting/binary_search.cpp	(revision 3595)
@@ -37,7 +37,4 @@
 	else{
 		while((beg <= end) && (*mid != target)){
-			#ifdef _ISSM_DEBUG_
-				printf("1: %i %i %i\n",*beg,*mid,*(end-1));
-			#endif
 			// is the target in lower or upper half?
 			if (target < *mid) {
@@ -49,7 +46,4 @@
 				mid = beg + (end-beg)/2;  //new middle
 			}
-			#ifdef _ISSM_DEBUG_
-				printf("2: %i %i %i\n",*beg,*mid,*(end-1));
-			#endif
 		}
 			  
Index: /issm/trunk/src/c/shared/TriMesh/TriMeshUtils.cpp
===================================================================
--- /issm/trunk/src/c/shared/TriMesh/TriMeshUtils.cpp	(revision 3594)
+++ /issm/trunk/src/c/shared/TriMesh/TriMeshUtils.cpp	(revision 3595)
@@ -161,7 +161,4 @@
 		if (el2!=-1){
 			/*el and el2 are on a segment rift, facing one another, plug them into riftsegments_uncompressed: */
-			#ifdef _ISSM_DEBUG_
-				printf("Elements %i and %i are on a rift\n",el+1,el2+1);
-			#endif
 		    *(riftsegments_uncompressed+5*i+0)=1;
 		    *(riftsegments_uncompressed+5*i+1)=el;
@@ -186,9 +183,4 @@
 	}
 
- 	#ifdef _ISSM_DEBUG_
-	for (i=0;i<nriftsegs;i++){
-		printf("Elements %i and %i are on a rift, joined by grids %i and %i\n",*(riftsegments+4*i+0)+1,*(riftsegments+4*i+1)+1,*(riftsegments+4*i+2),*(riftsegments+4*i+3));
-	}
-	#endif
 	xfree((void**)&riftsegments_uncompressed);
 	
@@ -217,11 +209,5 @@
 	/*Build a list of all the elements connected to this grid: */
 	GridElementsList(&GridElements,&NumGridElements,grid,index,nel);
-	#ifdef _ISSM_DEBUG_
-		printf("Connected elements for grid %i\n",grid);
-		for (k=0;k<NumGridElements;k++){
-			printf("El %i\n",GridElements[k]+1);
-		}
-	#endif
-			
+
 	/*Figure out the list of elements  that are on the same side of the rift. To do so, we start from one 
 	 * side of the rift and keep rotating in the same direction:*/
@@ -232,7 +218,4 @@
 															   done rotating*/
 	GridElementListOnOneSideOfRift[1]=*(riftsegments+4*segmentnumber+1);
-	#ifdef _ISSM_DEBUG_
-	printf("Starting with elements %i and %i for grid %i\n",GridElementListOnOneSideOfRift[0]+1,GridElementListOnOneSideOfRift[1]+1,grid);
-	#endif
 	counter=1;
 	for (;;){
@@ -240,8 +223,4 @@
 		 * equal to GridElementListOnOneSideOfRift[counter-1]*/
 		for (k=0;k<NumGridElements;k++){
-			#ifdef _ISSM_DEBUG_
-			printf("k: %i GridElements[k]: %i GridElementListOnOneSideOfRift[counter-1]: %i GridElementListOnOneSideOfRift[counter] %i Neighboor %i\n",k,GridElements[k]+1,
-				GridElementListOnOneSideOfRift[counter-1]+1,GridElementListOnOneSideOfRift[counter]+1,IsNeighbor(GridElements[k],GridElementListOnOneSideOfRift[counter],index));
-			#endif
 			if(IsNeighbor(GridElements[k],GridElementListOnOneSideOfRift[counter],index)){
 				/*Verify this element is not already in our list of element on the same side of the rift: */
@@ -268,11 +247,4 @@
 			GridElementListOnOneSideOfRift[l]=GridElementListOnOneSideOfRift[l+1];
 		}
-		#ifdef _ISSM_DEBUG_
-		printf("Grid %i is owned by  the following elements  on the same side of the rift: \n",grid);
-		for (l=0;l<NumGridElementListOnOneSideOfRift;l++){
-			printf("%i ",GridElementListOnOneSideOfRift[l]+1);
-		}
-		printf("\n");
-		#endif
 		break;
 	}// for (;;)
@@ -315,7 +287,4 @@
 		for (j=0;j<nsegs;j++){
 			if (*(segments+3*j+2)==(el1+1)){
-				#ifdef _ISSM_DEBUG_
-					printf("Segment %i is the same as rift segment %i\n",j,i);
-				#endif
 				/*segment j is the same as rift segment i.Let's update segments[j][:] using  element el1 and the corresponding rift segment.
 				 *Because riftsegments does not represent a list of rift segments anymore (it got heavily modified in SplitElementsForRifts, 
@@ -350,7 +319,4 @@
 			}
 			if (*(segments+3*j+2)==(el2+1)){
-				#ifdef _ISSM_DEBUG_
-					printf("Segment %i is the same as rift segment %i\n",j,i);
-				#endif
 				/*segment j is the same as rift segment i.*/
 				/*Let's update segments[j][:] using  element el2 and the corresponding rift segment: */
@@ -440,13 +406,7 @@
 	int n;
 	int el=-1;
-	#ifdef _ISSM_DEBUG_
-		printf("Looking for %lf %lf\n",A,B);
-	#endif
 	for (n=0;n<nel;n++){
 		if (((*(index+3*n+0)==A)  || (*(index+3*n+1)==A) || (*(index+3*n+2)==A) ) && ((*(index+3*n+0)==B)  || (*(index+3*n+1)==B) || (*(index+3*n+2)==B))){
 			el=n;
-			#ifdef _ISSM_DEBUG_
-				printf("Found them: %lf %lf %lf\n",*(index+3*n+0),*(index+3*n+1),*(index+3*n+2));
-			#endif
 			break;
 		}
@@ -921,7 +881,4 @@
 		}
 
-		#ifdef _ISSM_DEBUG_
-			printf("Tips for rift#%i   (%i-%i)\n",i,tip1,tip2);
-		#endif
 		/*Record tips in riftstips: */
 		*(riftstips+2*i+0)=(double)tip1;
@@ -960,10 +917,4 @@
 		}
 
-		#ifdef _ISSM_DEBUG_
-		for (j=0;j<numsegs;j++){
-			printf("%i\n",order[j]);
-		}
-		#endif
-
 		/*Using the order vector, and the riftsegments_copy and riftspairs_copy, reorder the segments and the pairs: */
 		for (j=0;j<numsegs;j++){
Index: /issm/trunk/src/c/toolkits/petsc/patches/MatPartition.cpp
===================================================================
--- /issm/trunk/src/c/toolkits/petsc/patches/MatPartition.cpp	(revision 3594)
+++ /issm/trunk/src/c/toolkits/petsc/patches/MatPartition.cpp	(revision 3595)
@@ -62,9 +62,5 @@
 		upper_row--;
 		range=upper_row-lower_row+1;
-		#ifdef _ISSM_DEBUG_
-		 	PetscSynchronizedPrintf(MPI_COMM_WORLD,"My rank %i Range %i",my_rank,range);
-			PetscSynchronizedFlush(MPI_COMM_WORLD);
-		#endif
-			
+
 		count=0;
 		if (range){
@@ -84,7 +80,4 @@
 		/*Now each node has a node_rows vectors holding which rows they should extract from matrixA. Create an Index Set from node_rows.*/
 		ISCreateGeneral(MPI_COMM_WORLD,count,node_rows,&row_index);
-		#ifdef _ISSM_DEBUG_
-			ISView(row_index,PETSC_VIEWER_STDOUT_WORLD);
-		#endif
 		
 		/*Same deal for columns*/
@@ -94,7 +87,4 @@
 		}
 		ISCreateGeneral(MPI_COMM_WORLD,col_partition_vector_size,node_cols,&col_index);
-		#ifdef _ISSM_DEBUG_
-			ISView(col_index,PETSC_VIEWER_STDOUT_WORLD);
-		#endif
 
 		/*Call MatGetSubMatrix*/
Index: /issm/trunk/src/c/toolkits/petsc/patches/MatlabMatrixToDoubleMatrix.cpp
===================================================================
--- /issm/trunk/src/c/toolkits/petsc/patches/MatlabMatrixToDoubleMatrix.cpp	(revision 3594)
+++ /issm/trunk/src/c/toolkits/petsc/patches/MatlabMatrixToDoubleMatrix.cpp	(revision 3595)
@@ -86,13 +86,4 @@
 		}
 		
-		#ifdef _ISSM_DEBUG_
-		for(i=0;i<rows;i++){
-			for(j=0;j<cols;j++){
-				printf("%g ",*(matrix+cols*i+j));
-			}
-			printf("\n");
-		}
-		#endif
-
 	}
 
Index: /issm/trunk/src/c/toolkits/petsc/patches/MatlabMatrixToPetscMatrix.cpp
===================================================================
--- /issm/trunk/src/c/toolkits/petsc/patches/MatlabMatrixToPetscMatrix.cpp	(revision 3594)
+++ /issm/trunk/src/c/toolkits/petsc/patches/MatlabMatrixToPetscMatrix.cpp	(revision 3595)
@@ -75,18 +75,8 @@
 	else{
 
-
 		/*Dealing with dense matrix: recover pointer and size: */
 		mxmatrix_ptr=(double*)mxGetPr(mxmatrix);
 		rows=mxGetM(mxmatrix);
 		cols=mxGetN(mxmatrix);
-
-		#ifdef _ISSM_DEBUG_
-		for(i=0;i<rows;i++){
-			for(j=0;j<cols;j++){
-				printf("%g ",*(mxmatrix_ptr+cols*j+i));
-			}
-			printf("\n");
-		}
-		#endif
 
 		/*transpose, as Petsc now does not allows MAT_COLUMN_ORIENTED matrices in MatSetValues: */
Index: /issm/trunk/src/c/toolkits/petsc/patches/PetscMatrixToMatlabMatrix.cpp
===================================================================
--- /issm/trunk/src/c/toolkits/petsc/patches/PetscMatrixToMatlabMatrix.cpp	(revision 3594)
+++ /issm/trunk/src/c/toolkits/petsc/patches/PetscMatrixToMatlabMatrix.cpp	(revision 3595)
@@ -85,10 +85,4 @@
 			MatGetRow(matrix,i,&ncols,&columns,&column_values);
 
-			#ifdef _ISSM_DEBUG_
-			for(j=0;j<ncols;j++){
-				printf("%i %i: %g\n",i,columns[j],column_values[j]);
-			}
-			#endif
-
 			/*copy values: */
 			if(ncols)memcpy( val+j, column_values,ncols*sizeof(double));
Index: /issm/trunk/src/c/toolkits/petsc/patches/VecMerge.cpp
===================================================================
--- /issm/trunk/src/c/toolkits/petsc/patches/VecMerge.cpp	(revision 3594)
+++ /issm/trunk/src/c/toolkits/petsc/patches/VecMerge.cpp	(revision 3595)
@@ -33,14 +33,4 @@
 	VecGetSize(B,&MB);
 
-	#ifdef _ISSM_DEBUG_
-		OutputServerMessages(0,0,"\n%s%s\n",STROFFSET,"Row partition vector:");
-		if (my_rank==0){
-			int i;
-			for (i=0;i<row_partition_vector_size;i++){
-				OutputServerMessages(0,0,"%s%lf\n",STROFFSET,*(row_partition_vector+i));
-			}
-		}
-	#endif
-
 	/*If the dimension of the partitioning vector is not the same as that of vector B, we have a problem: */
 	if ( (row_partition_size !=MB) ){
@@ -72,13 +62,4 @@
 	VecAssemblyEnd(A);
 
-	#ifdef _ISSM_DEBUG_
-		_printf_("Vector B:\n");
-		fflush(stdout);
-		VecView(B,PETSC_VIEWER_STDOUT_WORLD);
-		_printf_("Merge into  vector A:\n");
-		fflush(stdout);
-		VecView(A,PETSC_VIEWER_STDOUT_WORLD);
-	#endif
-
 	/*Free ressources:*/
 	xfree((void**)&idxm);
Index: /issm/trunk/src/c/toolkits/petsc/patches/VecPartition.cpp
===================================================================
--- /issm/trunk/src/c/toolkits/petsc/patches/VecPartition.cpp	(revision 3594)
+++ /issm/trunk/src/c/toolkits/petsc/patches/VecPartition.cpp	(revision 3595)
@@ -54,8 +54,4 @@
 		upper_row--;
 		range=upper_row-lower_row+1;
-		#ifdef _ISSM_DEBUG_
-		 	PetscSynchronizedPrintf(MPI_COMM_WORLD,"My rank %i Range %i\n",my_rank,range);
-			PetscSynchronizedFlush(MPI_COMM_WORLD);
-		#endif
 
 		if (range){
@@ -81,9 +77,4 @@
 			values=NULL;
 		}
-
-		#ifdef _ISSM_DEBUG_
-			PetscSynchronizedPrintf(MPI_COMM_WORLD,"My rank: %i My count: %i node_rows: %p values: %p\n",my_rank,count,node_rows,values);
-			PetscSynchronizedFlush(MPI_COMM_WORLD);
-		#endif 
 
 		if (count){
@@ -115,16 +106,6 @@
 		VecAssemblyBegin(outvector);
 		VecAssemblyEnd(outvector);
-
 		
 	}
-	
-	#ifdef _ISSM_DEBUG_
-		_printf_("Vector A:\n");
-		fflush(stdout);
-		VecView(vectorA,PETSC_VIEWER_STDOUT_WORLD);
-		_printf_("Partition of vector A:\n");
-		fflush(stdout);
-		VecView(outvector,PETSC_VIEWER_STDOUT_WORLD);
-	#endif
 	
 	/*Assign output pointers:*/
Index: /issm/trunk/src/c/toolkits/petsc/patches/VecToMPISerial.cpp
===================================================================
--- /issm/trunk/src/c/toolkits/petsc/patches/VecToMPISerial.cpp	(revision 3594)
+++ /issm/trunk/src/c/toolkits/petsc/patches/VecToMPISerial.cpp	(revision 3595)
@@ -64,7 +64,4 @@
 		if (my_rank==0){
 			MPI_Recv(buffer,3,MPI_INT,i,1,MPI_COMM_WORLD,&status); 
-			#ifdef _ISSM_DEBUG_
-				_printf_("Received from node %i: %i-%i\n",buffer[0],buffer[1],buffer[2]); 
-			#endif
 			if (buffer[2])MPI_Recv(gathered_vector+buffer[1],buffer[2],MPI_DOUBLE,i,1,MPI_COMM_WORLD,&status);
 		}
Index: /issm/trunk/src/c/toolkits/plapack/patches/CyclicalFactorization.cpp
===================================================================
--- /issm/trunk/src/c/toolkits/plapack/patches/CyclicalFactorization.cpp	(revision 3594)
+++ /issm/trunk/src/c/toolkits/plapack/patches/CyclicalFactorization.cpp	(revision 3595)
@@ -35,7 +35,4 @@
 	*pnprows=nprows;
 	*pnpcols=npcols;
-	#ifdef _ISSM_DEBUG_
-		_printf_("Decomposition: %i-%i\n",nprows,npcols);
-	#endif
 
 }
@@ -51,7 +48,4 @@
 	for (i=0;i<input;i++){
 		SmallestPrimeFactor(&prime_factor,*(decomp+i));
-		#ifdef _ISSM_DEBUG_
-			_printf_("Smallest prime factor for term %i : %i\n",i,prime_factor);
-		#endif
 		if (prime_factor==*(decomp+i)){
 			*pdecomp_size=i;
@@ -63,10 +57,4 @@
 		}
 	}
-	#ifdef _ISSM_DEBUG_
-		_printf_("Prime factor decomposition for integer %i: \n",input);
-		for(i=0;i<*pdecomp_size;i++){
-			_printf_("%i ",*(decomp+i));
-		}
-	#endif
 
 	*pdecomp=decomp;
Index: /issm/trunk/src/c/toolkits/plapack/patches/PlapackInvertMatrix.cpp
===================================================================
--- /issm/trunk/src/c/toolkits/plapack/patches/PlapackInvertMatrix.cpp	(revision 3594)
+++ /issm/trunk/src/c/toolkits/plapack/patches/PlapackInvertMatrix.cpp	(revision 3595)
@@ -118,9 +118,4 @@
 	PLA_API_end(); 
 
-	#ifdef _ISSM_DEBUG_
-		PLA_Global_show("Matrix A",a," %lf","Done with A");
-		MatView(*A,PETSC_VIEWER_STDOUT_WORLD);
-	#endif
-
 	/*Call the plapack invert routine*/
 	PLA_General_invert(PLA_METHOD_INV,a);
@@ -129,9 +124,4 @@
 	MatGetType(*A,&type);
 	PlapackToPetsc(inv_A,local_mA,local_nA,mA,nA,type,a,templ,nprows,npcols,nb);
-
-	#ifdef _ISSM_DEBUG_
-		PLA_Global_show("Inverse of A",a," %lf","Done...");
-		MatView(*inv_A,PETSC_VIEWER_STDOUT_WORLD);
-	#endif
 
 	/*Free ressources:*/