diff --git a/models/GreyscaleColorModel.cpp b/models/GreyscaleColorModel.cpp
index 6525c09c..d97c844d 100644
--- a/models/GreyscaleColorModel.cpp
+++ b/models/GreyscaleColorModel.cpp
@@ -596,15 +596,15 @@ void ScaLBL_GreyscaleColorModel::AssignGreyPoroPermLabels()
 
 	if (rank==0){
         printf("Image resolution: %.5g [um/voxel]\n",Dm->voxel_length);
-		printf("Number of component labels: %lu \n",NLABELS);
+		printf("Number of Grey-fluid labels: %lu \n",NLABELS);
 		for (unsigned int idx=0; idx<NLABELS; idx++){
 			VALUE=LabelList[idx];
 			POROSITY=PorosityList[idx];
 			PERMEABILITY=PermeabilityList[idx];
 			double volume_fraction  = double(label_count_global[idx])/double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
-			printf("   label=%d: porosity=%.3g, permeability=%.3g [um^2] (=%.3g [voxel^2]), volume fraction=%.3g\n",
+			printf("   grey-fluid label=%d, porosity=%.3g, permeability=%.3g [um^2] (=%.3g [voxel^2]), volume fraction=%.3g\n",
                     VALUE,POROSITY,PERMEABILITY,PERMEABILITY/Dm->voxel_length/Dm->voxel_length,volume_fraction); 
-            printf("             effective porosity=%.3g\n",volume_fraction*POROSITY);
+            printf("                        effective porosity=%.3g\n",volume_fraction*POROSITY);
 		}
         printf("The weighted porosity, considering both open and grey voxels, is %.3g\n",GreyPorosity);
 	}
diff --git a/models/GreyscaleModel.cpp b/models/GreyscaleModel.cpp
index f2f94f9d..85832a4b 100644
--- a/models/GreyscaleModel.cpp
+++ b/models/GreyscaleModel.cpp
@@ -274,7 +274,7 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
 
 	if (rank==0){
         printf("Image resolution: %.5g [um/voxel]\n",Dm->voxel_length);
-		printf("Component labels: %lu \n",NLABELS);
+		printf("Number of component labels: %lu \n",NLABELS);
 		for (unsigned int idx=0; idx<NLABELS; idx++){
 			VALUE=LabelList[idx];
 			POROSITY=PorosityList[idx];
diff --git a/tests/lbpm_greyscaleColor_simulator.cpp b/tests/lbpm_greyscaleColor_simulator.cpp
index 3f9fdd19..fec85c0e 100644
--- a/tests/lbpm_greyscaleColor_simulator.cpp
+++ b/tests/lbpm_greyscaleColor_simulator.cpp
@@ -6,53 +6,66 @@
 #include <stdexcept>
 #include <fstream>
 
-#include "common/ScaLBL.h"
-#include "common/Communication.h"
-#include "common/MPI.h"
 #include "models/GreyscaleColorModel.h"
+#include "common/Utilities.h"
 //#define WRITE_SURFACES
 
+//*************************************************************************
+// Implementation of Greyscale Two-Fluid Color LBM using CUDA
+//*************************************************************************
+
 using namespace std;
 
 
 int main(int argc, char **argv)
 {
-	//*****************************************
-	// ***** MPI STUFF ****************
-	//*****************************************
-	// Initialize MPI
-	int rank,nprocs;
-	MPI_Init(&argc,&argv);
-	MPI_Comm comm = MPI_COMM_WORLD;
-	MPI_Comm_rank(comm,&rank);
-	MPI_Comm_size(comm,&nprocs);
-	{
-		// parallel domain size (# of sub-domains)
-		int nprocx,nprocy,nprocz;
-		int iproc,jproc,kproc;
 
-		if (rank == 0){
-			printf("****************************************\n");
-			printf("Running Greyscale Two-Phase Calculation \n");
-			printf("****************************************\n");
-		}
-		// Initialize compute device
-		int device=ScaLBL_SetDevice(rank);
-		ScaLBL_DeviceBarrier();
-		MPI_Barrier(comm);
-		
-		ScaLBL_GreyscaleColorModel GreyscaleColor(rank,nprocs,comm);
-		auto filename = argv[1];
-		GreyscaleColor.ReadParams(filename);
-		GreyscaleColor.SetDomain();    // this reads in the domain 
-		GreyscaleColor.ReadInput();
-		GreyscaleColor.Create();       // creating the model will create data structure to match the pore structure and allocate variables
-		GreyscaleColor.Initialize();   // initializing the model will set initial conditions for variables
-		GreyscaleColor.Run();	 
-		GreyscaleColor.WriteDebug();
-	}
-	// ****************************************************
-	MPI_Barrier(comm);
-	MPI_Finalize();
-	// ****************************************************
+    // Initialize MPI and error handlers
+    Utilities::startup( argc, argv );
+
+    { // Limit scope so variables that contain communicators will free before MPI_Finialize
+  
+      MPI_Comm comm;
+      MPI_Comm_dup(MPI_COMM_WORLD,&comm);
+      int rank = comm_rank(comm);
+      int nprocs = comm_size(comm);
+  
+      if (rank == 0){
+	      printf("****************************************\n");
+	      printf("Running Greyscale Two-Phase Calculation \n");
+	      printf("****************************************\n");
+      }
+      // Initialize compute device
+      ScaLBL_SetDevice(rank);
+      ScaLBL_DeviceBarrier();
+      MPI_Barrier(comm);
+  
+      PROFILE_ENABLE(1);
+      //PROFILE_ENABLE_TRACE();
+      //PROFILE_ENABLE_MEMORY();
+      PROFILE_SYNCHRONIZE();
+      PROFILE_START("Main");
+      Utilities::setErrorHandlers();
+  
+      auto filename = argv[1];
+      ScaLBL_GreyscaleColorModel GreyscaleColor(rank,nprocs,comm);
+      GreyscaleColor.ReadParams(filename);
+      GreyscaleColor.SetDomain();    
+      GreyscaleColor.ReadInput();    
+      GreyscaleColor.Create();       // creating the model will create data structure to match the pore structure and allocate variables
+      GreyscaleColor.Initialize();   // initializing the model will set initial conditions for variables
+      GreyscaleColor.Run();	       
+      GreyscaleColor.WriteDebug();
+  
+      PROFILE_STOP("Main");
+      PROFILE_SAVE("lbpm_greyscaleColor_simulator",1);
+      // ****************************************************
+  
+      MPI_Barrier(comm);
+      MPI_Comm_free(&comm);
+  
+    } // Limit scope so variables that contain communicators will free before MPI_Finialize
+
+    Utilities::shutdown();
+
 }
diff --git a/tests/lbpm_greyscale_simulator.cpp b/tests/lbpm_greyscale_simulator.cpp
index 782ed5b5..df8cb3cb 100644
--- a/tests/lbpm_greyscale_simulator.cpp
+++ b/tests/lbpm_greyscale_simulator.cpp
@@ -6,58 +6,67 @@
 #include <stdexcept>
 #include <fstream>
 
-#include "common/ScaLBL.h"
-#include "common/Communication.h"
-#include "common/MPI_Helpers.h"
 #include "models/GreyscaleModel.h"
+#include "common/Utilities.h"
 //#define WRITE_SURFACES
 
-/*
- * Simulator for two-phase flow in porous media
- * James E. McClure 2013-2014
- */
+//****************************************************************
+// Implementation of Greyscale Single-Fluid LBM using CUDA
+//****************************************************************
 
 using namespace std;
 
 
 int main(int argc, char **argv)
 {
-	//*****************************************
-	// ***** MPI STUFF ****************
-	//*****************************************
-	// Initialize MPI
-	int rank,nprocs;
-	MPI_Init(&argc,&argv);
-	MPI_Comm comm = MPI_COMM_WORLD;
-	MPI_Comm_rank(comm,&rank);
-	MPI_Comm_size(comm,&nprocs);
-	{
-		// parallel domain size (# of sub-domains)
 
-		if (rank == 0){
-			printf("********************************************************\n");
-			printf("Running Greyscale Single Phase Permeability Calculation \n");
-			printf("********************************************************\n");
-		}
-		// Initialize compute device
-		int device=ScaLBL_SetDevice(rank);
-		NULL_USE(device);
-		ScaLBL_DeviceBarrier();
-		MPI_Barrier(comm);
-		
-		ScaLBL_GreyscaleModel Greyscale(rank,nprocs,comm);
-		auto filename = argv[1];
-		Greyscale.ReadParams(filename);
-		Greyscale.SetDomain();    // this reads in the domain 
-		Greyscale.ReadInput();
-		Greyscale.Create();       // creating the model will create data structure to match the pore structure and allocate variables
-		Greyscale.Initialize();   // initializing the model will set initial conditions for variables
-		Greyscale.Run();	 
-		Greyscale.VelocityField();
-		//Greyscale.WriteDebug();
-	}
-	// ****************************************************
-	MPI_Barrier(comm);
-	MPI_Finalize();
-	// ****************************************************
+    // Initialize MPI and error handlers
+    Utilities::startup( argc, argv );
+
+    { // Limit scope so variables that contain communicators will free before MPI_Finialize
+  
+      MPI_Comm comm;
+      MPI_Comm_dup(MPI_COMM_WORLD,&comm);
+      int rank = comm_rank(comm);
+      int nprocs = comm_size(comm);
+  
+      if (rank == 0){
+	      printf("********************************************************\n");
+	 	  printf("Running Greyscale Single Phase Permeability Calculation \n");
+	 	  printf("********************************************************\n");
+      }
+      // Initialize compute device
+      ScaLBL_SetDevice(rank);
+      ScaLBL_DeviceBarrier();
+      MPI_Barrier(comm);
+  
+      PROFILE_ENABLE(1);
+      //PROFILE_ENABLE_TRACE();
+      //PROFILE_ENABLE_MEMORY();
+      PROFILE_SYNCHRONIZE();
+      PROFILE_START("Main");
+      Utilities::setErrorHandlers();
+  
+      auto filename = argv[1];
+      ScaLBL_GreyscaleModel Greyscale(rank,nprocs,comm);
+      Greyscale.ReadParams(filename);
+      Greyscale.SetDomain();    
+      Greyscale.ReadInput();    
+      Greyscale.Create();       // creating the model will create data structure to match the pore structure and allocate variables
+      Greyscale.Initialize();   // initializing the model will set initial conditions for variables
+      Greyscale.Run();	       
+	  Greyscale.VelocityField();
+      //Greyscale.WriteDebug();
+  
+      PROFILE_STOP("Main");
+      PROFILE_SAVE("lbpm_greyscale_simulator",1);
+      // ****************************************************
+  
+      MPI_Barrier(comm);
+      MPI_Comm_free(&comm);
+  
+    } // Limit scope so variables that contain communicators will free before MPI_Finialize
+
+    Utilities::shutdown();
+
 }