diff --git a/Exec/LyA/inputs b/Exec/LyA/inputs
index 66d4c8f3..51144b63 100644
--- a/Exec/LyA/inputs
+++ b/Exec/LyA/inputs
@@ -73,7 +73,6 @@ nyx.comoving_h   = 0.71d0
 
 # PARTICLES
 nyx.do_dm_particles = 1
-nyx.do_tiling = 1
 
 # >>>>>>>>>>>>>  PARTICLE INIT OPTIONS <<<<<<<<<<<<<<<<
 #  "AsciiFile"        "Random"	    "Cosmological"
diff --git a/Exec/LyA/inputs.rt b/Exec/LyA/inputs.rt
index 18e64339..4ecdc100 100644
--- a/Exec/LyA/inputs.rt
+++ b/Exec/LyA/inputs.rt
@@ -74,7 +74,6 @@ nyx.comoving_h   = 0.71d0
 
 # PARTICLES
 nyx.do_dm_particles = 1
-nyx.do_tiling = 1
 
 # >>>>>>>>>>>>>  PARTICLE INIT OPTIONS <<<<<<<<<<<<<<<<
 #  "AsciiFile"        "Random"	    "Cosmological"
diff --git a/Exec/Make.Nyx b/Exec/Make.Nyx
index 9aaf8436..4308625e 100644
--- a/Exec/Make.Nyx
+++ b/Exec/Make.Nyx
@@ -102,7 +102,7 @@ VPATH_LOCATIONS   += $(Blocs)
 
 #These are the directories in AMReX
 
-Pdirs   := Base AmrCore Amr Boundary Particle Extern/amrdata
+Pdirs   := Base AmrCore Amr Boundary Particle Extern/amrdata 
 Ppack	+= $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir)/Make.package)
 Plocs	+= $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir))
 
diff --git a/Exec/Scaling/inputs b/Exec/Scaling/inputs
index ea7bd1d2..2c6b1d5c 100644
--- a/Exec/Scaling/inputs
+++ b/Exec/Scaling/inputs
@@ -83,8 +83,6 @@ nyx.comoving_h   = 0.702d0
 
 # PARTICLES
 nyx.do_dm_particles = 1
-particles.do_tiling = 1
-particles.tile_size = 2048 4 4
 
 # >>>>>>>>>>>>>  PARTICLE INIT OPTIONS <<<<<<<<<<<<<<<<
 #  "AsciiFile"        "Random"	    "Cosmological"
diff --git a/Source/Gravity/Gravity.cpp b/Source/Gravity/Gravity.cpp
index de5988bb..2965563b 100644
--- a/Source/Gravity/Gravity.cpp
+++ b/Source/Gravity/Gravity.cpp
@@ -143,11 +143,11 @@ Gravity::read_params ()
         Real Gconst;
         fort_get_grav_const(&Gconst);
         Ggravity = -4.0 * M_PI * Gconst;
-        if (verbose > 0 && ParallelDescriptor::IOProcessor())
+        if (verbose > 0)
         {
-            std::cout << "Getting Gconst from nyx_constants: " << Gconst
+            amrex::Print() << "Getting Gconst from nyx_constants: " << Gconst
                       << '\n';
-            std::cout << "Using " << Ggravity << " for 4 pi G in Gravity.cpp "
+            amrex::Print() << "Using " << Ggravity << " for 4 pi G in Gravity.cpp "
                       << '\n';
         }
         done = true;
@@ -158,8 +158,8 @@ void
 Gravity::install_level (int       level,
                         AmrLevel* level_data_to_install)
 {
-    if (verbose > 1 && ParallelDescriptor::IOProcessor())
-        std::cout << "Installing Gravity level " << level << '\n';
+    if (verbose > 1)
+        amrex::Print() << "Installing Gravity level " << level << '\n';
 
     LevelData[level] = level_data_to_install;
 
@@ -280,8 +280,8 @@ Gravity::solve_for_old_phi (int               level,
         return;
 #endif
 
-    if (verbose && ParallelDescriptor::IOProcessor())
-        std::cout << "Gravity ... single level solve for old phi at level "
+    if (verbose)
+        amrex::Print() << "Gravity ... single level solve for old phi at level "
                   << level << std::endl;
     MultiFab Rhs(grids[level], dmap[level], 1, 0);
     Rhs.setVal(0.0);
@@ -315,8 +315,8 @@ Gravity::solve_for_new_phi (int               level,
         return;
 #endif
 
-    if (verbose && ParallelDescriptor::IOProcessor())
-        std::cout << "Gravity ... single level solve for new phi at level "
+    if (verbose)
+        amrex::Print() << "Gravity ... single level solve for new phi at level "
                   << level << std::endl;
 
     MultiFab Rhs(grids[level], dmap[level], 1, 0);
@@ -348,8 +348,8 @@ Gravity::solve_for_phi (int               level,
 
 {
     BL_PROFILE("Gravity::solve_for_phi()");
-    if (verbose && ParallelDescriptor::IOProcessor())
-        std::cout << " ... solve for phi at level " << level << '\n';
+    if (verbose)
+        amrex::Print() << " ... solve for phi at level " << level << '\n';
 
     // This is a correction for fully periodic domains only
     if (Geometry::isAllPeriodic())
@@ -449,11 +449,9 @@ Gravity::solve_for_phi (int               level,
 
             Rhs.plus(-sum, 0, 1, 0);
 
-            if (verbose && ParallelDescriptor::IOProcessor())
-                std::cout << " ... subtracting "
-                          << sum
-                          << " to ensure solvability "
-                          << '\n';
+            if (verbose)
+                amrex::Print() << " ... subtracting " << sum
+                               << " to ensure solvability " << '\n';
         }
     }
 
@@ -503,12 +501,10 @@ Gravity::solve_for_delta_phi (int                        crse_level,
     BL_ASSERT(grad_delta_phi.size() == num_levels);
     BL_ASSERT(delta_phi.size() == num_levels);
 
-    if (verbose && ParallelDescriptor::IOProcessor())
+    if (verbose)
     {
-        std::cout << "... solving for delta_phi at crse_level = " << crse_level
-                  << '\n';
-        std::cout << "...                    up to fine_level = " << fine_level
-                  << '\n';
+        amrex::Print() << "... solving for delta_phi at crse_level = " << crse_level << '\n';
+        amrex::Print() << "...                    up to fine_level = " << fine_level << '\n';
     }
 
     const Geometry& geom = parent->Geom(crse_level);
@@ -595,8 +591,8 @@ Gravity::solve_for_delta_phi (int                        crse_level,
 
        local_correction = local_correction / grids[crse_level].numPts();
 
-       if (verbose && ParallelDescriptor::IOProcessor())
-          std::cout << "WARNING: Adjusting RHS in solve_for_delta_phi by " << local_correction << std::endl;
+       if (verbose)
+          amrex::Print() << "WARNING: Adjusting RHS in solve_for_delta_phi by " << local_correction << std::endl;
 
        for (int lev = crse_level; lev <= fine_level; lev++) {
          Rhs_p[lev-crse_level]->plus(-local_correction,0,1,0);
@@ -631,10 +627,10 @@ Gravity::gravity_sync (int crse_level, int fine_level, int iteration, int ncycle
     BL_PROFILE("Gravity::gravity_sync()");
     BL_ASSERT(parent->finestLevel()>crse_level);
 
-    if (verbose && ParallelDescriptor::IOProcessor())
+    if (verbose)
     {
-        std::cout << " ... gravity_sync at crse_level " << crse_level << '\n';
-        std::cout << " ...     up to finest_level     " << fine_level << '\n';
+        amrex::Print() << " ... gravity_sync at crse_level " << crse_level << '\n';
+        amrex::Print() << " ...     up to finest_level     " << fine_level << '\n';
     }
 
     // Build Rhs for solve for delta_phi
@@ -659,8 +655,8 @@ Gravity::gravity_sync (int crse_level, int fine_level, int iteration, int ncycle
 
         local_correction /= grids[crse_level].numPts();
 
-        if (verbose && ParallelDescriptor::IOProcessor())
-            std::cout << "WARNING: Adjusting RHS in gravity_sync solve by " << local_correction << '\n';
+        if (verbose)
+            amrex::Print() << "WARNING: Adjusting RHS in gravity_sync solve by " << local_correction << '\n';
 
         crse_rhs.plus(-local_correction,0,1,0);
     }
@@ -851,9 +847,9 @@ Gravity::multilevel_solve_for_new_phi (int level,
                                        int use_previous_phi_as_guess)
 {
     BL_PROFILE("Gravity::multilevel_solve_for_new_phi()");
-    if (verbose && ParallelDescriptor::IOProcessor())
-        std::cout << "Gravity ... multilevel solve for new phi at base level " << level
-                  << " to finest level " << finest_level << '\n';
+    if (verbose)
+        amrex::Print() << "Gravity ... multilevel solve for new phi at base level " << level
+                       << " to finest level " << finest_level << '\n';
 
     for (int lev = level; lev <= finest_level; lev++)
     {
@@ -877,9 +873,9 @@ Gravity::multilevel_solve_for_old_phi (int level,
                                        int use_previous_phi_as_guess)
 {
     BL_PROFILE("Gravity::multilevel_solve_for_old_phi()");
-    if (verbose && ParallelDescriptor::IOProcessor())
-        std::cout << "Gravity ... multilevel solve for old phi at base level " << level
-                  << " to finest level " << finest_level << '\n';
+    if (verbose)
+        amrex::Print() << "Gravity ... multilevel solve for old phi at base level " << level
+                       << " to finest level " << finest_level << '\n';
 
     for (int lev = level; lev <= finest_level; lev++)
     {
@@ -912,7 +908,6 @@ Gravity::actual_multilevel_solve (int                       level,
                                   int                       use_previous_phi_as_guess)
 {
     BL_PROFILE("Gravity::actual_multilevel_solve()");
-    const int IOProc = ParallelDescriptor::IOProcessorNumber();
 
     const int num_levels = finest_level - level + 1;
 
@@ -1047,16 +1042,15 @@ Gravity::actual_multilevel_solve (int                       level,
     // This correction is for fully periodic domains only.
     if (Geometry::isAllPeriodic())
     {
-        if (verbose && ParallelDescriptor::IOProcessor())
-            std::cout << " ... subtracting average density " << mass_offset
-                      << " from RHS at each level " << '\n';
+        if (verbose)
+            amrex::Print() << " ... subtracting average density " << mass_offset
+                           << " from RHS at each level " << '\n';
 
         for (int lev = 0; lev < num_levels; lev++)
             for (MFIter mfi(*Rhs_p[lev]); mfi.isValid(); ++mfi)
                 (*Rhs_p[lev])[mfi].plus(-mass_offset);
 
        // This is used to enforce solvability if appropriate.
-//     if ( grids[level].contains(parent->Geom(level).Domain()) )
        if ( parent->Geom(level).Domain().numPts() == grids[level].numPts() )
        {
            Real sum = 0;
@@ -1071,17 +1065,12 @@ Gravity::actual_multilevel_solve (int                       level,
             const Real eps = 1.e-10 * std::abs(mass_offset);
             if (std::abs(sum) > eps)
             {
-                 if (ParallelDescriptor::IOProcessor())
-                 {
-                     std::cout << " ... current avg differs from mass_offset by " << sum << " " << '\n';
-                     std::cout << " ... Gravity::actual_multilevel_solve -- total mass has changed!"
-                               << '\n';;
-                 }
+               amrex::Print() << " ... current avg differs from mass_offset by " << sum << " " << '\n';
+               amrex::Print() << " ... Gravity::actual_multilevel_solve -- total mass has changed!" << '\n';;
             }
 
-            if (verbose && ParallelDescriptor::IOProcessor())
-                std::cout << " ... subtracting " << sum << " to ensure solvability "
-                          << '\n';
+            if (verbose)
+                amrex::Print() << " ... subtracting " << sum << " to ensure solvability " << '\n';
 
             for (int lev = 0; lev < num_levels; lev++)
                 (*Rhs_p[lev]).plus(-sum, 0, 1, 0);
@@ -1671,23 +1660,19 @@ Gravity::set_mass_offset (Real time)
 
         mass_offset /= geom.ProbSize();
 
-        if (verbose && ParallelDescriptor::IOProcessor())
-            std::cout << "Gravity ... defining average density in Gravity::set_mass_offset to be " << mass_offset
-                      << '\n';
+        if (verbose)
+            amrex::Print() << "Gravity ... defining average density in Gravity::set_mass_offset to be " 
+                           << mass_offset << '\n';
     }
 
     const Real diff = std::abs(mass_offset - old_mass_offset);
     const Real eps  = 1.e-10 * std::abs(old_mass_offset);
     if (diff > eps && old_mass_offset > 0)
     {
-        if (ParallelDescriptor::IOProcessor())
-        {
-            std::cout << " ... new vs old mass_offset " << mass_offset << " "
-                      << old_mass_offset
-                      << " ... diff is " << diff <<  '\n';
-            std::cout << " ... Gravity::set_mass_offset -- total mass has changed!"
-                      << '\n';;
-        }
+        amrex::Print() << " ... new vs old mass_offset " << mass_offset << " "
+                       << old_mass_offset << " ... diff is " << diff <<  '\n';
+        amrex::Print() << " ... Gravity::set_mass_offset -- total mass has changed!"
+                       << '\n';;
     }
 }
 
@@ -1896,9 +1881,9 @@ void
 Gravity::CorrectRhsUsingOffset(int level, MultiFab& Rhs)
 {
     BL_PROFILE("Gravity::CorrectRhsUsingOffset()");
-    if (verbose && ParallelDescriptor::IOProcessor())
-        std::cout << " ... subtracting average density from RHS in solve ... "
-                  << mass_offset << '\n';
+    if (verbose)
+        amrex::Print() << " ... subtracting average density from RHS in solve ... "
+                       << mass_offset << '\n';
 
     for (MFIter mfi(Rhs); mfi.isValid(); ++mfi)
         Rhs[mfi].plus(-mass_offset);
@@ -1912,17 +1897,13 @@ Gravity::CorrectRhsUsingOffset(int level, MultiFab& Rhs)
         const Real eps = 1.e-10 * std::abs(mass_offset);
         if (std::abs(sum) > eps)
         {
-            if (ParallelDescriptor::IOProcessor())
-            {
-                std::cout << " ... current avg differs from mass_offset by " << sum << " " << '\n';
-                std::cout << " ... Gravity : single_level solve for phi -- total mass has changed!"
+            amrex::Print() << " ... current avg differs from mass_offset by " << sum << " " << '\n';
+            amrex::Print() << " ... Gravity : single_level solve for phi -- total mass has changed!"
                           << '\n';;
-            }
         }
 
-        if (verbose && ParallelDescriptor::IOProcessor())
-            std::cout << " ... subtracting " << sum << " to ensure solvability "
-                      << '\n';
+        if (verbose)
+            amrex::Print() << " ... subtracting " << sum << " to ensure solvability " << '\n';
 
         Rhs.plus(-sum, 0, 1, 0);
     }
@@ -2021,7 +2002,8 @@ Gravity::solve_with_HPGMG(int level,
     {
       if (ParallelDescriptor::IOProcessor())
       {
-        std::cerr << "WARNING: Periodic boundary conditions, but f does not sum to zero... mean(f)=" << average_value_of_f << std::endl;
+        std::cerr << "WARNING: Periodic boundary conditions, but f does not sum to zero... mean(f)=" 
+                << average_value_of_f << std::endl;
         std::cerr << "Subtracting mean(f) from RHS ..." << std::endl;
       }
       shift_vector(&level_h,VECTOR_F,VECTOR_F,-average_value_of_f);
@@ -2031,8 +2013,7 @@ Gravity::solve_with_HPGMG(int level,
   rebuild_operator(&level_h,NULL,a,b);    // i.e. calculate Dinv and lambda_max
   MGBuild(&MG_h,&level_h,a,b,minCoarseDim,ParallelDescriptor::Communicator()); // build the Multigrid Hierarchy
   //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-  if (ParallelDescriptor::IOProcessor())
-    std::cout << std::endl << std::endl << "===== STARTING SOLVE =====" << std::endl << std::flush;
+  amrex::Print() << std::endl << std::endl << "===== STARTING SOLVE =====" << std::endl << std::flush;
 
   MGResetTimers (&MG_h);
   zero_vector (MG_h.levels[0], VECTOR_U);
diff --git a/Source/Nyx.cpp b/Source/Nyx.cpp
index 0ef310d2..ccf5057c 100644
--- a/Source/Nyx.cpp
+++ b/Source/Nyx.cpp
@@ -804,12 +804,11 @@ Nyx::est_time_step (Real dt_old)
         Real a = get_comoving_a(cur_time);
         const Real* dx = geom.CellSize();
 
-	Real dt = est_dt;
-
 #ifdef _OPENMP
-#pragma omp parallel firstprivate(dt)
+#pragma omp parallel reduction(min:est_dt)
 #endif
 	{
+          Real dt = 1.e200;
 	  for (MFIter mfi(stateMF,true); mfi.isValid(); ++mfi)
 	    {
 	      const Box& box = mfi.tilebox();
@@ -818,12 +817,7 @@ Nyx::est_time_step (Real dt_old)
                 (BL_TO_FORTRAN(stateMF[mfi]), box.loVect(), box.hiVect(), dx,
                  &dt, &a);
 	    }
-#ifdef _OPENMP
-#pragma omp critical (nyx_estdt)
-#endif
-	  {
-	    est_dt = std::min(est_dt, dt);
-	  }
+          est_dt = std::min(est_dt, dt);
 	}
 
         // If in comoving coordinates, then scale dt (based on u and c) by a
diff --git a/Source/NyxParticles.cpp b/Source/NyxParticles.cpp
index 7ea25208..3d73594d 100644
--- a/Source/NyxParticles.cpp
+++ b/Source/NyxParticles.cpp
@@ -385,12 +385,12 @@ Nyx::init_particles ()
                 amrex::Abort("Nyx::init_particles(): particle_initrandom_iseed must be > 0");
             }
 
-            if (verbose && ParallelDescriptor::IOProcessor())
+            if (verbose)
             {
-                std::cout << "\nInitializing DM with cloud of "
-                          << particle_initrandom_count
-                          << " random particles with initial seed: "
-                          << particle_initrandom_iseed << "\n\n";
+                amrex::Print() << "\nInitializing DM with cloud of "
+                               << particle_initrandom_count
+                               << " random particles with initial seed: "
+                               << particle_initrandom_iseed << "\n\n";
             }
 
             DMPC->InitRandom(particle_initrandom_count,
@@ -408,12 +408,10 @@ Nyx::init_particles ()
                 amrex::Abort("Nyx::init_particles(): particle_initrandom_iseed must be > 0");
             }
 
-            if (verbose && ParallelDescriptor::IOProcessor())
-            {
-                std::cout << "\nInitializing DM with of " << particle_initrandom_count_per_box
-                          << " random particles per box with initial seed: "
-                          << particle_initrandom_iseed << "\n\n";
-            }
+            if (verbose)
+                amrex::Print() << "\nInitializing DM with of " << particle_initrandom_count_per_box
+                               << " random particles per box with initial seed: "
+                               << particle_initrandom_iseed << "\n\n";
 
             DMPC->InitRandomPerBox(particle_initrandom_count_per_box,
                                    particle_initrandom_iseed, pdata);
@@ -421,21 +419,21 @@ Nyx::init_particles ()
         }
         else if (particle_init_type == "RandomPerCell")
         {
-            if (verbose && ParallelDescriptor::IOProcessor())
-                std::cout << "\nInitializing DM with 1 random particle per cell " << "\n";
+            if (verbose)
+                amrex::Print() << "\nInitializing DM with 1 random particle per cell " << "\n";
 
             int n_per_cell = 1;
             DMPC->InitNRandomPerCell(n_per_cell, pdata);
         }
         else if (particle_init_type == "AsciiFile")
         {
-            if (verbose && ParallelDescriptor::IOProcessor())
+            if (verbose)
             {
-                std::cout << "\nInitializing DM particles from \""
-                          << ascii_particle_file << "\" ...\n\n";
+                amrex::Print() << "\nInitializing DM particles from \""
+                               << ascii_particle_file << "\" ...\n\n";
                 if (init_with_sph_particles == 1)
-                    std::cout << "\nInitializing SPH particles from ascii \""
-                              << sph_particle_file << "\" ...\n\n";
+                    amrex::Print() << "\nInitializing SPH particles from ascii \""
+                                   << sph_particle_file << "\" ...\n\n";
             }
             //
             // The second argument is how many Reals we read into `m_data[]`
@@ -448,13 +446,13 @@ Nyx::init_particles ()
         }
         else if (particle_init_type == "BinaryFile")
         {
-            if (verbose && ParallelDescriptor::IOProcessor())
+            if (verbose)
             {
-                std::cout << "\nInitializing DM particles from \""
-                          << binary_particle_file << "\" ...\n\n";
+                amrex::Print() << "\nInitializing DM particles from \""
+                               << binary_particle_file << "\" ...\n\n";
                 if (init_with_sph_particles == 1)
-                    std::cout << "\nInitializing SPH particles from binary \""
-                              << sph_particle_file << "\" ...\n\n";
+                    amrex::Print() << "\nInitializing SPH particles from binary \""
+                                   << sph_particle_file << "\" ...\n\n";
             }
             //
             // The second argument is how many Reals we read into `m_data[]`
@@ -467,13 +465,13 @@ Nyx::init_particles ()
         }
         else if (particle_init_type == "BinaryMetaFile")
         {
-            if (verbose && ParallelDescriptor::IOProcessor())
+            if (verbose)
             {
-                std::cout << "\nInitializing DM particles from meta file\""
-                          << binary_particle_file << "\" ...\n\n";
+                amrex::Print() << "\nInitializing DM particles from meta file\""
+                               << binary_particle_file << "\" ...\n\n";
                 if (init_with_sph_particles == 1)
-                    std::cout << "\nInitializing SPH particles from meta file\""
-                              << sph_particle_file << "\" ...\n\n";
+                    amrex::Print() << "\nInitializing SPH particles from meta file\""
+                                   << sph_particle_file << "\" ...\n\n";
             }
             //
             // The second argument is how many Reals we read into `m_data[]`
@@ -567,9 +565,9 @@ Nyx::init_particles ()
         NPC->SetVerbose(particle_verbose);
         if (particle_init_type == "AsciiFile")
         {
-            if (verbose && ParallelDescriptor::IOProcessor())
-                std::cout << "\nInitializing Neutrino particles from \""
-                          << neutrino_particle_file << "\" ...\n\n";
+            if (verbose)
+                amrex::Print() << "\nInitializing Neutrino particles from \""
+                               << neutrino_particle_file << "\" ...\n\n";
             //
             // The second argument is how many Reals we read into `m_data[]`
             // after reading in `m_pos[]`. Here we're reading in the particle
@@ -579,10 +577,10 @@ Nyx::init_particles ()
         }
         else if (particle_init_type == "BinaryFile")
         {
-            if (verbose && ParallelDescriptor::IOProcessor())
+            if (verbose)
             {
-                std::cout << "\nInitializing Neutrino particles from \""
-                          << neutrino_particle_file << "\" ...\n\n";
+                amrex::Print() << "\nInitializing Neutrino particles from \""
+                               << neutrino_particle_file << "\" ...\n\n";
             }
             //
             // The second argument is how many Reals we read into `m_data[]`
@@ -609,10 +607,8 @@ Nyx::init_santa_barbara (int init_sb_vels)
     Real cur_time = state[State_Type].curTime();
     Real a = old_a;
 
-    if (ParallelDescriptor::IOProcessor()) {
-        std::cout << "... time and comoving a when data is initialized at level " 
-                  << level << " " << cur_time << " " << a << '\n';
-    }
+    amrex::Print() << "... time and comoving a when data is initialized at level " 
+                   << level << " " << cur_time << " " << a << '\n';
 
     if (level == 0)
     {
@@ -864,22 +860,22 @@ Nyx::particle_est_time_step (Real& est_dt)
 	}
 #endif
 
-        if (verbose && ParallelDescriptor::IOProcessor())
+        if (verbose)
         {
             if (est_dt_particle > 0)
             {
-                std::cout << "...estdt from particles at level "
+                amrex::Print() << "...estdt from particles at level "
                           << level << ": " << est_dt_particle << '\n';
             }
             else
             {
-                std::cout << "...there are no particles at level "
+                amrex::Print() << "...there are no particles at level "
                           << level << '\n';
             }
 #ifdef NEUTRINO_PARTICLES
             if (est_dt_neutrino > 0)
             {
-                std::cout << "...estdt from neutrinos at level "
+                amrex::Print() << "...estdt from neutrinos at level "
                           << level << ": " << est_dt_neutrino << '\n';
             }
 #endif
@@ -910,7 +906,10 @@ Nyx::particle_redistribute (int lbase, bool init)
         static Array<BoxArray>            ba;
         static Array<DistributionMapping> dm;
 
-        bool changed = false;
+        bool    changed      = false;
+        bool dm_changed      = false;
+        bool ba_changed      = false;
+        bool ba_size_changed = false;
 
         int flev = parent->finestLevel();
 	
@@ -920,27 +919,36 @@ Nyx::particle_redistribute (int lbase, bool init)
  
         if (ba.size() != flev+1)
         {
+            amrex::Print() << "BA SIZE " << ba.size() << std::endl;
+            amrex::Print() << "FLEV    " << flev << std::endl;
             ba.resize(flev+1);
             dm.resize(flev+1);
             changed = true;
+            ba_size_changed = true;
         }
         else
         {
             for (int i = 0; i <= flev && !changed; i++)
             {
                 if (ba[i] != parent->boxArray(i))
+                {
                     //
                     // The BoxArrays have changed in the regridding.
                     //
                     changed = true;
+                    ba_changed = true;
+                }
 
                 if ( ! changed)
                 {
                     if (dm[i] != parent->getLevel(i).get_new_data(0).DistributionMap())
+                    {
                         //
                         // The DistributionMaps have changed in the regridding.
                         //
                         changed = true;
+                        dm_changed = true;
+                    }
                 }
             }
         }
@@ -953,8 +961,15 @@ Nyx::particle_redistribute (int lbase, bool init)
 	    // because of we called redistribute during a subcycle, there may be particles not in
 	    // the proper position on coarser levels.
             //
-            if (verbose && ParallelDescriptor::IOProcessor())
-                std::cout << "Calling redistribute because changed " << '\n';
+            if (verbose)
+            {
+                if (ba_size_changed)
+                   amrex::Print() << "Calling redistribute because the size of BoxArray changed " << '\n';
+                else if (ba_changed)
+                   amrex::Print() << "Calling redistribute because BoxArray changed " << '\n';
+                else if (dm_changed)
+                   amrex::Print() << "Calling redistribute because DistMap changed " << '\n';
+            }
 
             DMPC->Redistribute(lbase);
             //
@@ -968,8 +983,8 @@ Nyx::particle_redistribute (int lbase, bool init)
         }
         else
         {
-            if (verbose && ParallelDescriptor::IOProcessor())
-                std::cout << "NOT calling redistribute because NOT changed " << '\n';
+            if (verbose)
+                amrex::Print() << "NOT calling redistribute because NOT changed " << '\n';
         }
     }
 }
@@ -1064,12 +1079,6 @@ Nyx::NyxParticlesAddProcsToComp(Amr *aptr, int nSidecarProcs, int prevSidecarPro
                     int ioProcNumSCS, int ioProcNumAll, int scsMyId,
 		                        MPI_Comm scsCommn)
 {
-#if 0
-// What is this doing here???
-if(ParallelDescriptor::IOProcessor()) {
-  std::cout << "PPPPPPPP:  DMPC SPC APC NPC = " << DMPC << "  " << SPC << "  " << APC << "  " << NPC << std::endl;
-}
-#endif
 }
 
 //NyxParticleContainerBase::~NyxParticleContainerBase() {}
diff --git a/Source/Nyx_hydro.cpp b/Source/Nyx_hydro.cpp
index b792fc75..604ac4a4 100644
--- a/Source/Nyx_hydro.cpp
+++ b/Source/Nyx_hydro.cpp
@@ -129,7 +129,7 @@ Nyx::just_the_hydro (Real time,
         strang_first_step(time,dt,S_old_tmp,D_old_tmp);
 
 #ifdef _OPENMP
-#pragma omp parallel
+#pragma omp parallel reduction(max:courno)
 #endif
        {
        FArrayBox flux[BL_SPACEDIM], u_gdnv[BL_SPACEDIM];
@@ -182,12 +182,7 @@ Nyx::just_the_hydro (Real time,
         ke_added += ske;
        } // end of MFIter loop
 
-#ifdef _OPENMP
-#pragma omp critical (hydro_courno)
-#endif
-       {
         courno = std::max(courno, cflLoc);
-       }
 
        } // end of omp parallel region
 
diff --git a/Source/particle_mod.f90 b/Source/particle_mod.f90
index 09dfa2a5..90ed3238 100644
--- a/Source/particle_mod.f90
+++ b/Source/particle_mod.f90
@@ -1,6 +1,6 @@
 module particle_mod
 
-  use amrex_fort_module, only: c_real => amrex_real
+  use amrex_fort_module, only: c_real => amrex_particle_real
   use iso_c_binding ,    only: c_int
 
   implicit none