diff --git a/src/boundary_conditions/incflo_set_bcs.cpp b/src/boundary_conditions/incflo_set_bcs.cpp
index 993822d95..6d930080c 100644
--- a/src/boundary_conditions/incflo_set_bcs.cpp
+++ b/src/boundary_conditions/incflo_set_bcs.cpp
@@ -89,7 +89,7 @@ incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector<amrex::BCRec> const& bcs,
             if (m_bc_type[olo] == BC::mixed) {
                 prob_set_BC_MF(olo, blo, bc_arr, lev, inflow, outflow, field);
             } else {
-                amrex::ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     for (int n = 0; n < ncomp; n++) {
                         bc_arr(i,j,k,n) = bc_ptr[n].lo(dir);
@@ -99,7 +99,7 @@ incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector<amrex::BCRec> const& bcs,
             if (m_bc_type[ohi] == BC::mixed) {
                 prob_set_BC_MF(ohi, bhi, bc_arr, lev, inflow, outflow, field);
             } else {
-                amrex::ParallelFor(bhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     for (int n = 0; n < ncomp; n++) {
                         bc_arr(i,j,k,n) = bc_ptr[n].hi(dir);
@@ -189,7 +189,7 @@ incflo::make_robinBC_MFs(int lev, MultiFab* state)
 
 #ifdef AMREX_USE_EB
 void
-incflo::set_eb_velocity (int lev, amrex::Real /*time*/, MultiFab& eb_vel, int nghost)
+incflo::set_eb_velocity (int lev, Real /*time*/, MultiFab& eb_vel, int nghost)
 {
     Geometry const& gm = Geom(lev);
     eb_vel.setVal(0.);
@@ -224,7 +224,7 @@ incflo::set_eb_velocity (int lev, amrex::Real /*time*/, MultiFab& eb_vel, int ng
 
           bool has_comps = false;
           Real eb_vel_mag(0.0);
-          GpuArray<amrex::Real,3> eb_vel_comps{0.};
+          GpuArray<Real,3> eb_vel_comps{0.};
           // Flow is specified as a velocity magnitude
           if ( m_eb_flow.is_mag ) {
              eb_vel_mag = m_eb_flow.vel_mag;
@@ -281,7 +281,7 @@ incflo::set_eb_velocity (int lev, amrex::Real /*time*/, MultiFab& eb_vel, int ng
 }
 
 void
-incflo::set_eb_density (int lev, amrex::Real /*time*/, MultiFab& eb_density, int nghost)
+incflo::set_eb_density (int lev, Real /*time*/, MultiFab& eb_density, int nghost)
 {
     Geometry const& gm = Geom(lev);
     eb_density.setVal(0.);
@@ -350,7 +350,7 @@ incflo::set_eb_density (int lev, amrex::Real /*time*/, MultiFab& eb_density, int
 }
 
 void
-incflo::set_eb_tracer (int lev, amrex::Real /*time*/, MultiFab& eb_tracer, int nghost)
+incflo::set_eb_tracer (int lev, Real /*time*/, MultiFab& eb_tracer, int nghost)
 {
     Geometry const& gm = Geom(lev);
     eb_tracer.setVal(0.);
diff --git a/src/convection/incflo_compute_MAC_projected_velocities.cpp b/src/convection/incflo_compute_MAC_projected_velocities.cpp
index 31680fde9..1fd291920 100644
--- a/src/convection/incflo_compute_MAC_projected_velocities.cpp
+++ b/src/convection/incflo_compute_MAC_projected_velocities.cpp
@@ -36,7 +36,7 @@ incflo::compute_MAC_projected_velocities (
                     Array4<Real const> const& rho      = density[lev]->array(mfi);
                     Array4<Real const> const& divtau   = ld.divtau_o.const_array(mfi);
                     if (m_advect_momentum) {
-                        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                         {
                             AMREX_D_TERM(vel_f(i,j,k,0) += divtau(i,j,k,0)/rho(i,j,k);,
                                          vel_f(i,j,k,1) += divtau(i,j,k,1)/rho(i,j,k);,
@@ -44,7 +44,7 @@ incflo::compute_MAC_projected_velocities (
                         });
                     }
                     else {
-                        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                         {
                             AMREX_D_TERM(vel_f(i,j,k,0) += divtau(i,j,k,0);,
                                          vel_f(i,j,k,1) += divtau(i,j,k,1);,
@@ -80,7 +80,7 @@ incflo::compute_MAC_projected_velocities (
         EB_interp_CellCentroid_to_FaceCentroid (*density[lev], GetArrOfPtrs(inv_rho[lev]),
                                                 0, 0, 1, geom[lev], get_density_bcrec());
 #else
-        amrex::average_cellcenter_to_face(GetArrOfPtrs(inv_rho[lev]), *density[lev], geom[lev]);
+        average_cellcenter_to_face(GetArrOfPtrs(inv_rho[lev]), *density[lev], geom[lev]);
 #endif
 
         for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
diff --git a/src/convection/incflo_compute_advection_term.cpp b/src/convection/incflo_compute_advection_term.cpp
index 0d64f979c..30a4c8f8d 100644
--- a/src/convection/incflo_compute_advection_term.cpp
+++ b/src/convection/incflo_compute_advection_term.cpp
@@ -154,7 +154,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
                     Array4<Real const> const& rho      = density[lev]->array(mfi);
                     Array4<Real const> const& divtau   = ld.divtau_o.const_array(mfi);
                     if (m_advect_momentum) {
-                        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                         {
                             AMREX_D_TERM(vel_f(i,j,k,0) += divtau(i,j,k,0)/rho(i,j,k);,
                                          vel_f(i,j,k,1) += divtau(i,j,k,1)/rho(i,j,k);,
@@ -162,7 +162,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
                         });
                     }
                     else {
-                        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                         {
                             AMREX_D_TERM(vel_f(i,j,k,0) += divtau(i,j,k,0);,
                                          vel_f(i,j,k,1) += divtau(i,j,k,1);,
@@ -345,7 +345,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
                 Array4<Real const> rho     = density[lev]->const_array(mfi);
                 Array4<Real      > rho_vel =  rhovel[lev].array(mfi);
 
-                amrex::ParallelFor(bxg, AMREX_SPACEDIM,
+                ParallelFor(bxg, AMREX_SPACEDIM,
                 [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                 {
                     rho_vel(i,j,k,n) = rho(i,j,k) * U(i,j,k,n);
@@ -358,7 +358,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
                     Array4<Real const> vf        = vel_forces[lev]->const_array(mfi);
                     Array4<Real      > rho_vel_f =  rhovel_f.array();
 
-                    amrex::ParallelFor(fbx, AMREX_SPACEDIM,
+                    ParallelFor(fbx, AMREX_SPACEDIM,
                     [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                     {
                         rho_vel_f(i,j,k,n) = rho(i,j,k) * vf(i,j,k,n);
@@ -456,7 +456,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
                 if ( any_conserv_trac ) {
                     trac_tmp = rhotrac[lev].array(mfi);
 
-                    amrex::ParallelFor(bxg, m_ntrac,
+                    ParallelFor(bxg, m_ntrac,
                     [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                     {
                         if ( iconserv[n] ){
@@ -560,7 +560,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
                 Array4<Real const> rho     =  get_density_eb()[lev]->const_array(mfi);
                 Array4<Real      > rho_vel =  rhovel_fab.array();
 
-                amrex::ParallelFor(bxg, AMREX_SPACEDIM,
+                ParallelFor(bxg, AMREX_SPACEDIM,
                 [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                 {
                     rho_vel(i,j,k,n) = rho(i,j,k) * U(i,j,k,n);
@@ -741,7 +741,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
                                   *density[lev], bc_den, lev);
             } else {
                 auto const& drdt = conv_r[lev]->array(mfi);
-                amrex::ParallelFor(bx,
+                ParallelFor(bx,
                 [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     drdt(i,j,k) = 0.;
diff --git a/src/derive/incflo_derive.cpp b/src/derive/incflo_derive.cpp
index f7e6d940d..66ca35171 100644
--- a/src/derive/incflo_derive.cpp
+++ b/src/derive/incflo_derive.cpp
@@ -73,7 +73,7 @@ void incflo::compute_strainrate_at_level (int /*lev*/,
                 auto typ = flag_fab.getType(bx);
                 if (typ == FabType::covered)
                 {
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         sr_arr(i,j,k) = Real(0.0);
                     });
@@ -81,7 +81,7 @@ void incflo::compute_strainrate_at_level (int /*lev*/,
                 else if (typ == FabType::singlevalued)
                 {
                     auto const& flag_arr = flag_fab.const_array();
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         sr_arr(i,j,k) = incflo_strainrate_eb(i,j,k,AMREX_D_DECL(idx,idy,idz),vel_arr,flag_arr(i,j,k));
                     });
@@ -89,7 +89,7 @@ void incflo::compute_strainrate_at_level (int /*lev*/,
                 else
 #endif
                 {
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         sr_arr(i,j,k) = incflo_strainrate(i,j,k,AMREX_D_DECL(idx,idy,idz),vel_arr);
                     });
@@ -171,7 +171,7 @@ void incflo::ComputeMagVel (int /*lev*/,
         auto typ = flags.getType(bx);
         if (typ == FabType::covered)
         {
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 magvel_fab(i,j,k) = Real(0.0);
             });
@@ -240,7 +240,7 @@ void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab
         auto typ = flags.getType(bx);
         if (typ == FabType::covered)
         {
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 vort_fab(i,j,k) = Real(0.0);
             });
@@ -248,7 +248,7 @@ void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab
         else if (typ == FabType::singlevalued)
         {
             const auto& flag_fab = flags.const_array();
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 constexpr Real c0 = Real(-1.5);
                 constexpr Real c1 = Real( 2.0);
@@ -306,7 +306,7 @@ void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab
         else
 #endif
         {
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 Real vx = Real(0.5) * (ccvel_fab(i+1,j,k,1) - ccvel_fab(i-1,j,k,1)) * idx;
                 Real uy = Real(0.5) * (ccvel_fab(i,j+1,k,0) - ccvel_fab(i,j-1,k,0)) * idy;
@@ -343,7 +343,7 @@ void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab
         auto typ = flags.getType(bx);
         if (typ == FabType::covered)
         {
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 vort_fab(i,j,k) = Real(0.0);
             });
@@ -351,7 +351,7 @@ void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab
         else if (typ == FabType::singlevalued)
         {
             const auto& flag_fab = flags.const_array();
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 constexpr Real c0 = Real(-1.5);
                 constexpr Real c1 = Real( 2.0);
@@ -447,7 +447,7 @@ void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab
         else
 #endif
         {
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 Real vx = Real(0.5) * (ccvel_fab(i+1,j,k,1) - ccvel_fab(i-1,j,k,1)) * idx;
                 Real wx = Real(0.5) * (ccvel_fab(i+1,j,k,2) - ccvel_fab(i-1,j,k,2)) * idx;
diff --git a/src/derive/incflo_error.cpp b/src/derive/incflo_error.cpp
index c41f8dd4d..1d31a9468 100644
--- a/src/derive/incflo_error.cpp
+++ b/src/derive/incflo_error.cpp
@@ -20,7 +20,7 @@ void incflo::DiffFromExact (int /*lev*/, Geometry& lev_geom, Real time, Real dt,
             // When we enter this routine, this holds the computed solution
             Array4<Real> const& err = error.array(mfi);
 
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 constexpr Real  twopi = Real(2.0)*Real(3.1415926535897932);
                 constexpr Real fourpi = Real(4.0)*Real(3.1415926535897932);
@@ -75,7 +75,7 @@ void incflo::DiffFromExact (int /*lev*/, Geometry& lev_geom, Real time, Real dt,
 
             Real omega = pi * pi * visc_coef;
 
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
 
                 Real x = Real(i+0.5)*dx[0];
@@ -117,7 +117,7 @@ void incflo::DiffFromExact (int /*lev*/, Geometry& lev_geom, Real time, Real dt,
             // When we enter this routine, this holds the computed solution
             Array4<Real> const& err = error.array(mfi);
 
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 Real x = Real(i+0.5)*dx[0];
                 Real y = Real(j+0.5)*dx[1];
diff --git a/src/diffusion/DiffusionScalarOp.cpp b/src/diffusion/DiffusionScalarOp.cpp
index 032da0c11..1c5c8bf2d 100644
--- a/src/diffusion/DiffusionScalarOp.cpp
+++ b/src/diffusion/DiffusionScalarOp.cpp
@@ -260,8 +260,7 @@ DiffusionScalarOp::diffuse_scalar (Vector<MultiFab*> const& tracer,
                     Array4<Real> const& rhs_a = rhs[lev].array(mfi);
                     Array4<Real const> const& tra_a = tracer[lev]->const_array(mfi,comp);
                     Array4<Real const> const& rho_a = density[lev]->const_array(mfi);
-                    amrex::ParallelFor(bx,
-                    [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         rhs_a(i,j,k) = rho_a(i,j,k) * tra_a(i,j,k);
                     });
@@ -412,8 +411,7 @@ DiffusionScalarOp::diffuse_vel_components (Vector<MultiFab*> const& vel,
                 Array4<Real> const& rhs_a = rhs[lev].array(mfi);
                 Array4<Real const> const& vel_a = vel[lev]->const_array(mfi,comp);
                 Array4<Real const> const& rho_a = density[lev]->const_array(mfi);
-                amrex::ParallelFor(bx,
-                [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     rhs_a(i,j,k) = rho_a(i,j,k) * vel_a(i,j,k);
                 });
@@ -714,8 +712,7 @@ void DiffusionScalarOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
                 Box const& bx = mfi.tilebox();
                 Array4<Real> const& divtau_arr = a_divtau[lev]->array(mfi);
                 Array4<Real const> const& rho_arr = a_density[lev]->const_array(mfi);
-                amrex::ParallelFor(bx,
-                [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     Real rhoinv = Real(1.0)/rho_arr(i,j,k);
                     AMREX_D_TERM(divtau_arr(i,j,k,0) *= rhoinv;,
diff --git a/src/diffusion/DiffusionTensorOp.cpp b/src/diffusion/DiffusionTensorOp.cpp
index ab77f57c7..7aa030d70 100644
--- a/src/diffusion/DiffusionTensorOp.cpp
+++ b/src/diffusion/DiffusionTensorOp.cpp
@@ -164,7 +164,7 @@ DiffusionTensorOp::diffuse_velocity (Vector<MultiFab*> const& velocity,
             Array4<Real> const& rhs_a = rhs[lev].array(mfi);
             Array4<Real const> const& vel_a = velocity[lev]->const_array(mfi);
             Array4<Real const> const& rho_a = density[lev]->const_array(mfi);
-            amrex::ParallelFor(bx,AMREX_SPACEDIM,
+            ParallelFor(bx,AMREX_SPACEDIM,
             [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
             {
                 rhs_a(i,j,k,n) = rho_a(i,j,k) * vel_a(i,j,k,n);
@@ -301,8 +301,7 @@ void DiffusionTensorOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
                 Box const& bx = mfi.tilebox();
                 Array4<Real> const& divtau_arr = a_divtau[lev]->array(mfi);
                 Array4<Real const> const& rho_arr = a_density[lev]->const_array(mfi);
-                    amrex::ParallelFor(bx,
-                [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     Real rhoinv = Real(1.0)/rho_arr(i,j,k);
                     AMREX_D_TERM(divtau_arr(i,j,k,0) *= rhoinv;,
diff --git a/src/diffusion/incflo_diffusion.cpp b/src/diffusion/incflo_diffusion.cpp
index ae51d74f8..37e0494a0 100644
--- a/src/diffusion/incflo_diffusion.cpp
+++ b/src/diffusion/incflo_diffusion.cpp
@@ -325,14 +325,14 @@ incflo::fixup_eta_on_domain_faces (int lev, Array<MultiFab,AMREX_SPACEDIM>& fc,
         if (!gm.isPeriodic(idim)) {
             Array4<Real> const& fca = fc[idim].array(mfi);
             if (bx.smallEnd(idim) == domain.smallEnd(idim)) {
-                amrex::ParallelFor(amrex::bdryLo(bx, idim),
+                ParallelFor(amrex::bdryLo(bx, idim),
                 [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     fca(i,j,k) = cca(i,j,k);
                 });
             }
             if (bx.bigEnd(idim) == domain.bigEnd(idim)) {
-                amrex::ParallelFor(amrex::bdryHi(bx, idim),
+                ParallelFor(amrex::bdryHi(bx, idim),
                 [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     fca(i,j,k) = cca(i-1,j,k);
@@ -344,14 +344,14 @@ incflo::fixup_eta_on_domain_faces (int lev, Array<MultiFab,AMREX_SPACEDIM>& fc,
         if (!gm.isPeriodic(idim)) {
             Array4<Real> const& fca = fc[idim].array(mfi);
             if (bx.smallEnd(idim) == domain.smallEnd(idim)) {
-                amrex::ParallelFor(amrex::bdryLo(bx, idim),
+                ParallelFor(amrex::bdryLo(bx, idim),
                 [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     fca(i,j,k) = cca(i,j,k);
                 });
             }
             if (bx.bigEnd(idim) == domain.bigEnd(idim)) {
-                amrex::ParallelFor(amrex::bdryHi(bx, idim),
+                ParallelFor(amrex::bdryHi(bx, idim),
                 [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     fca(i,j,k) = cca(i,j-1,k);
@@ -364,14 +364,14 @@ incflo::fixup_eta_on_domain_faces (int lev, Array<MultiFab,AMREX_SPACEDIM>& fc,
         if (!gm.isPeriodic(idim)) {
             Array4<Real> const& fca = fc[idim].array(mfi);
             if (bx.smallEnd(idim) == domain.smallEnd(idim)) {
-                amrex::ParallelFor(amrex::bdryLo(bx, idim),
+                ParallelFor(amrex::bdryLo(bx, idim),
                 [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     fca(i,j,k) = cca(i,j,k);
                 });
             }
             if (bx.bigEnd(idim) == domain.bigEnd(idim)) {
-                amrex::ParallelFor(amrex::bdryHi(bx, idim),
+                ParallelFor(amrex::bdryHi(bx, idim),
                 [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     fca(i,j,k) = cca(i,j,k-1);
diff --git a/src/incflo.H b/src/incflo.H
index 8863c4d5c..b2d459260 100644
--- a/src/incflo.H
+++ b/src/incflo.H
@@ -316,6 +316,7 @@ public:
 
     static bool SteadyStateReached ();
 
+    void InitialPressureProjection ();
 
 private:
 
@@ -920,7 +921,6 @@ private:
     void ReadIOParameters ();
     void ResizeArrays (); // Resize arrays to fit (up to) max_level + 1 AMR levels
     void InitialProjection ();
-    void InitialPressureProjection ();
     void InitialIterations ();
 #ifdef AMREX_USE_EB
     void InitialRedistribution ();
diff --git a/src/incflo_compute_forces.cpp b/src/incflo_compute_forces.cpp
index 4df2a15c5..8ad6af7b1 100644
--- a/src/incflo_compute_forces.cpp
+++ b/src/incflo_compute_forces.cpp
@@ -19,7 +19,7 @@ void incflo::compute_tra_forces (Vector<MultiFab*> const& tra_forces,
                 Array4<Real>       const& tra_f = tra_forces[lev]->array(mfi);
                 Array4<Real const> const& rho   =    density[lev]->const_array(mfi);
 
-                amrex::ParallelFor(bx, m_ntrac,
+                ParallelFor(bx, m_ntrac,
                 [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                 {
                     // For now we don't have any external forces on the scalars
@@ -75,7 +75,7 @@ void incflo::compute_vel_forces_on_level (int lev,
                 //      first tracer rather than density
                 Array4<Real const> const& tra_o = tracer_old.const_array(mfi);
                 Array4<Real const> const& tra_n = tracer_new.const_array(mfi);
-                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     int n = 0; // Potential temperature
 
@@ -96,7 +96,7 @@ void incflo::compute_vel_forces_on_level (int lev,
 
             } else if (m_probtype == 16) {
                 Real Re = 1./m_mu;  // Note this assumes you are running exactly the problem set up, with U = 1 and L = 1 and rho = 1.
-                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     Real rhoinv = Real(1.0)/rho(i,j,k);
 
@@ -131,7 +131,7 @@ void incflo::compute_vel_forces_on_level (int lev,
                     vel_f(i,j,k,1) += 8.0 / Re * (24.0 * capF + 2.0 * fp * gpp + fppp * g) + 64.0 * (capF2 * capG1 - g * gp * capF1);
                 });
             } else {
-                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     Real rhoinv = Real(1.0)/rho(i,j,k);
 
diff --git a/src/incflo_correct_small_cells.cpp b/src/incflo_correct_small_cells.cpp
index a5c0453eb..47a0755b3 100644
--- a/src/incflo_correct_small_cells.cpp
+++ b/src/incflo_correct_small_cells.cpp
@@ -50,7 +50,7 @@ incflo::incflo_correct_small_cells (Vector<MultiFab*      > const& vel_in,
 
              if (!m_eb_flow.enabled) {
                 // This FAB has cut cells -- we define the centroid value in terms of the MAC velocities onfaces
-                amrex::ParallelFor(bx,
+                ParallelFor(bx,
                   [vfrac_fab,AMREX_D_DECL(apx_fab,apy_fab,apz_fab),ccvel_fab,AMREX_D_DECL(umac_fab,vmac_fab,wmac_fab)]
                   AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
@@ -73,7 +73,7 @@ incflo::incflo_correct_small_cells (Vector<MultiFab*      > const& vel_in,
                 Array4<Real const> const& eb_vel = get_velocity_eb()[lev]->const_array(mfi);
 
                 // This FAB has cut cells -- we define the centroid value in terms of the MAC velocities onfaces
-                amrex::ParallelFor(bx,
+                ParallelFor(bx,
                   [vfrac_fab,AMREX_D_DECL(apx_fab,apy_fab,apz_fab),ccvel_fab,AMREX_D_DECL(umac_fab,vmac_fab,wmac_fab),eb_vel]
                   AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
diff --git a/src/incflo_explicit_update.cpp b/src/incflo_explicit_update.cpp
index ea6ca31ef..bb636e619 100644
--- a/src/incflo_explicit_update.cpp
+++ b/src/incflo_explicit_update.cpp
@@ -31,7 +31,7 @@ void incflo::tracer_explicit_update (Vector<MultiFab> const& tra_forces)
             {
                 Array4<Real const> const& laps_o = ld.laps_o.const_array(mfi);
 
-                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     // If conservative
                     // (rho trac)^new = (rho trac)^old + dt * (
@@ -56,7 +56,7 @@ void incflo::tracer_explicit_update (Vector<MultiFab> const& tra_forces)
             {
                 Array4<Real const> const& laps_o = ld.laps_o.const_array(mfi);
 
-                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     for (int n = 0; n < l_ntrac; ++n)
                     {
@@ -73,7 +73,7 @@ void incflo::tracer_explicit_update (Vector<MultiFab> const& tra_forces)
             }
             else if (m_diff_type == DiffusionType::Implicit)
             {
-                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     for (int n = 0; n < l_ntrac; ++n)
                     {
@@ -123,7 +123,7 @@ void incflo::tracer_explicit_update_corrector (Vector<MultiFab> const& tra_force
                 Array4<Real const> const& laps_o = ld.laps_o.const_array(mfi);
                 Array4<Real const> const& laps   = ld.laps.const_array(mfi);
 
-                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     for (int n = 0; n < l_ntrac; ++n)
                     {
@@ -147,7 +147,7 @@ void incflo::tracer_explicit_update_corrector (Vector<MultiFab> const& tra_force
             {
                 Array4<Real const> const& laps_o = ld.laps_o.const_array(mfi);
 
-                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     for (int n = 0; n < l_ntrac; ++n)
                     {
@@ -169,7 +169,7 @@ void incflo::tracer_explicit_update_corrector (Vector<MultiFab> const& tra_force
             }
             else if (m_diff_type == DiffusionType::Implicit)
             {
-                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     for (int n = 0; n < l_ntrac; ++n)
                     {
diff --git a/src/incflo_redistribute.cpp b/src/incflo_redistribute.cpp
index 37f7f827a..343919628 100644
--- a/src/incflo_redistribute.cpp
+++ b/src/incflo_redistribute.cpp
@@ -79,8 +79,7 @@ incflo::redistribute_term ( MFIter const& mfi,
         // This is scratch space if calling StateRedistribute
         //  but is used as the weights (here set to 1) if calling
         //  FluxRedistribute
-        amrex::ParallelFor(Box(scratch),
-        [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(Box(scratch), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             scratch(i,j,k) = 1.;
         });
@@ -95,8 +94,7 @@ incflo::redistribute_term ( MFIter const& mfi,
     }
     else
     {
-        amrex::ParallelFor(bx, ncomp,
-        [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+        ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
         {
             out(i,j,k,n) = in(i,j,k,n);
         });
diff --git a/src/incflo_tagging.cpp b/src/incflo_tagging.cpp
index 875574d53..7ab6d8ca9 100644
--- a/src/incflo_tagging.cpp
+++ b/src/incflo_tagging.cpp
@@ -69,8 +69,7 @@ void incflo::ErrorEst (int levc, TagBoxArray& tags, Real time, int /*ngrow*/)
             Array4<Real const> const& rho = m_leveldata[levc]->density.const_array(mfi);
             Real rhoerr = tag_rho ? rhoerr_v[levc]: std::numeric_limits<Real>::max();
             Real gradrhoerr = tag_gradrho ? gradrhoerr_v[levc] : std::numeric_limits<Real>::max();
-            amrex::ParallelFor(bx,
-            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 if (tag_rho && rho(i,j,k) > rhoerr) {
                     tag(i,j,k) = tagval;
@@ -105,8 +104,7 @@ void incflo::ErrorEst (int levc, TagBoxArray& tags, Real time, int /*ngrow*/)
 
 #if (AMREX_SPACEDIM == 2)
 
-            amrex::ParallelFor(bx,
-            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                  Real x = problo[0] + (i+0.5)*l_dx;
                  Real y = problo[1] + (j+0.5)*l_dy;
@@ -122,8 +120,7 @@ void incflo::ErrorEst (int levc, TagBoxArray& tags, Real time, int /*ngrow*/)
             Real zlo = tag_region_lo[2];
             Real zhi = tag_region_hi[2];
 
-            amrex::ParallelFor(bx,
-            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                  Real x = problo[0] + Real(i+0.5)*l_dx;
                  Real y = problo[1] + Real(j+0.5)*l_dy;
@@ -189,7 +186,7 @@ void incflo::ErrorEst (int levc, TagBoxArray& tags, Real time, int /*ngrow*/)
             auto const&  mf_arr = mf->const_array(mfi);
             auto const& tag_arr = tags.array(mfi);
 
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 if (mf_arr(i,j,k) > 0) {
                     tag_arr(i,j,k) = tagval;
diff --git a/src/incflo_update_velocity.cpp b/src/incflo_update_velocity.cpp
index bd5a5375d..281329b19 100644
--- a/src/incflo_update_velocity.cpp
+++ b/src/incflo_update_velocity.cpp
@@ -46,7 +46,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                     Array4<Real const> const& divtau_o = ld.divtau_o.const_array(mfi);
                     // Here divtau_o is the difference of tensor and scalar divtau_o!
                     if (m_advect_momentum) {
-                        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                         {
                             AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);,
                                          vel(i,j,k,1) *= rho_old(i,j,k);,
@@ -59,7 +59,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                                          vel(i,j,k,2) /= rho_new(i,j,k););
                         });
                     } else {
-                        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                         {
                             AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+divtau_o(i,j,k,0));,
                                          vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+divtau_o(i,j,k,1));,
@@ -68,7 +68,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                     }
                 } else {
                     if (m_advect_momentum) {
-                        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                         {
                             AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);,
                                          vel(i,j,k,1) *= rho_old(i,j,k);,
@@ -81,7 +81,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                                          vel(i,j,k,2) /= rho_new(i,j,k););
                         });
                     } else {
-                        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                         {
                             AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0));,
                                          vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1));,
@@ -95,7 +95,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
 
                 Array4<Real const> const& divtau_o = ld.divtau_o.const_array(mfi);
                 if (m_advect_momentum) {
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);,
                                      vel(i,j,k,1) *= rho_old(i,j,k);,
@@ -108,7 +108,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                                      vel(i,j,k,2) /= rho_new(i,j,k););
                     });
                 } else {
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+l_half*divtau_o(i,j,k,0));,
                                      vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+l_half*divtau_o(i,j,k,1));,
@@ -120,7 +120,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
             {
                 Array4<Real const> const& divtau_o = ld.divtau_o.const_array(mfi);
                 if (m_advect_momentum) {
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         AMREX_D_TERM(vel(i,j,k,0) *= rho_old(i,j,k);,
                                      vel(i,j,k,1) *= rho_old(i,j,k);,
@@ -133,7 +133,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                                      vel(i,j,k,2) /= rho_new(i,j,k););
                     });
                 } else {
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         AMREX_D_TERM(vel(i,j,k,0) += l_dt*(dvdt(i,j,k,0)+vel_f(i,j,k,0)+divtau_o(i,j,k,0));,
                                      vel(i,j,k,1) += l_dt*(dvdt(i,j,k,1)+vel_f(i,j,k,1)+divtau_o(i,j,k,1));,
@@ -178,7 +178,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                 Array4<Real const> const& divtau   = ld.divtau.const_array(mfi);
 
                 if (m_advect_momentum) {
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         AMREX_D_TERM(vel(i,j,k,0) = rho_old(i,j,k) * vel_o(i,j,k,0) + l_dt * (
                                                     l_half*(  dvdt_o(i,j,k,0)+  dvdt(i,j,k,0))
@@ -198,7 +198,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                                      vel(i,j,k,2) /= rho_new(i,j,k););
                     });
                 } else {
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         AMREX_D_TERM(vel(i,j,k,0) = vel_o(i,j,k,0) + l_dt * (
                                                     l_half*(  dvdt_o(i,j,k,0)+  dvdt(i,j,k,0))
@@ -220,7 +220,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                 Array4<Real const> const& divtau_o = ld.divtau_o.const_array(mfi);
 
                 if (m_advect_momentum) {
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         AMREX_D_TERM(vel(i,j,k,0) = rho_old(i,j,k) * vel_o(i,j,k,0) + l_dt * (
                                                     l_half*(  dvdt_o(i,j,k,0)+dvdt(i,j,k,0))
@@ -237,7 +237,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                                      vel(i,j,k,2) /= rho_new(i,j,k););
                     });
                 } else {
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         AMREX_D_TERM(vel(i,j,k,0) = vel_o(i,j,k,0) + l_dt * (
                                                     l_half*(  dvdt_o(i,j,k,0)+dvdt(i,j,k,0))
@@ -258,7 +258,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                     Array4<Real const> const& divtau   = ld.divtau.const_array(mfi);
                     if (m_advect_momentum)
                     {
-                        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                         {
                             AMREX_D_TERM(vel(i,j,k,0) = rho_old(i,j,k) * vel_o(i,j,k,0) + l_dt * (
                                                         l_half*(dvdt_o(i,j,k,0) +   dvdt(i,j,k,0))
@@ -275,7 +275,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                                          vel(i,j,k,2) /= rho_new(i,j,k););
                         });
                     } else {
-                        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                         {
                             AMREX_D_TERM(vel(i,j,k,0) = vel_o(i,j,k,0) + l_dt * (
                                                         l_half*(  dvdt_o(i,j,k,0) +   dvdt(i,j,k,0))
@@ -291,7 +291,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                 } else {
                     if (m_advect_momentum)
                     {
-                        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                         {
                             AMREX_D_TERM(vel(i,j,k,0) = rho_old(i,j,k) * vel_o(i,j,k,0) + l_dt * (
                                                         l_half*(  dvdt_o(i,j,k,0)+dvdt(i,j,k,0))
@@ -307,7 +307,7 @@ void incflo::update_velocity (StepType step_type, Vector<MultiFab>& vel_eta, Vec
                                          vel(i,j,k,2) /= rho_new(i,j,k););
                         });
                     } else {
-                        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                         {
                             AMREX_D_TERM(vel(i,j,k,0) = vel_o(i,j,k,0) + l_dt * (
                                                         l_half*(  dvdt_o(i,j,k,0)+dvdt(i,j,k,0)) + vel_f(i,j,k,0) );,
diff --git a/src/prob/prob_bc.cpp b/src/prob/prob_bc.cpp
index cc3fe345e..2545f635f 100644
--- a/src/prob/prob_bc.cpp
+++ b/src/prob/prob_bc.cpp
@@ -28,7 +28,7 @@ void incflo::prob_set_BC_MF (Orientation const& ori, Box const& bx,
 
         Orientation::Side side = ori.faceDir();
         if (side == Orientation::low) {
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 for ( int n = 0; n < ncomp; n++ ){
                     if (direction == 0) {
@@ -55,7 +55,7 @@ void incflo::prob_set_BC_MF (Orientation const& ori, Box const& bx,
                 }
             });
         } else {
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 for ( int n = 0; n < ncomp; n++ ){
                     if (direction == 0) {
@@ -114,7 +114,7 @@ void incflo::prob_set_MAC_robinBCs (Orientation const& ori, Box const& bx,
 
         Orientation::Side side = ori.faceDir();
         if (side == Orientation::low) {
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 robin_f(i,j,k) = 0.;
 
@@ -147,7 +147,7 @@ void incflo::prob_set_MAC_robinBCs (Orientation const& ori, Box const& bx,
                 }
             });
         } else {
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 robin_f(i,j,k) = 0.;
 
@@ -213,7 +213,7 @@ void incflo::prob_set_diffusion_robinBCs (Orientation const& ori, Box const& bx,
 
         Orientation::Side side = ori.faceDir();
         if (side == Orientation::low) {
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 if (direction == 0) {
                     if (i <= half_num_cells) {
@@ -250,7 +250,7 @@ void incflo::prob_set_diffusion_robinBCs (Orientation const& ori, Box const& bx,
                 }
             });
         } else {
-            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
                 robin_f(i,j,k,0) = 0.;
 
diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp
index 93cecb35d..e776b3c87 100644
--- a/src/prob/prob_init_fluid.cpp
+++ b/src/prob/prob_init_fluid.cpp
@@ -201,7 +201,7 @@ void incflo::init_rotating_flow (Box const& vbx, Box const& /*gbx*/,
                                  GpuArray<Real, AMREX_SPACEDIM> const& /*problo*/,
                                  GpuArray<Real, AMREX_SPACEDIM> const& /*probhi*/)
 {
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         Real x = Real(i+0.5)*dx[0] - 0.5;
         Real y = Real(j+0.5)*dx[1] - 0.5;
@@ -245,7 +245,7 @@ void incflo::init_taylor_green (Box const& vbx, Box const& /*gbx*/,
                                 GpuArray<Real, AMREX_SPACEDIM> const& /*problo*/,
                                 GpuArray<Real, AMREX_SPACEDIM> const& /*probhi*/)
 {
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         Real x = Real(i+0.5)*dx[0];
         Real y = Real(j+0.5)*dx[1];
@@ -267,7 +267,7 @@ void incflo::init_taylor_green3d (Box const& vbx, Box const& /*gbx*/,
                                   GpuArray<Real, AMREX_SPACEDIM> const& /*problo*/,
                                   GpuArray<Real, AMREX_SPACEDIM> const& /*probhi*/)
 {
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         Real x = Real(i+0.5)*dx[0];
         Real y = Real(j+0.5)*dx[1];
@@ -292,7 +292,7 @@ void incflo::init_taylor_vortex (Box const& vbx, Box const& /*gbx*/,
                                  GpuArray<Real, AMREX_SPACEDIM> const& /*problo*/,
                                  GpuArray<Real, AMREX_SPACEDIM> const& /*probhi*/)
 {
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         Real x = Real(i+0.5)*dx[0];
         Real y = Real(j+0.5)*dx[1];
@@ -317,7 +317,7 @@ void incflo::init_vortex_in_sphere (Box const& vbx, Box const& /*gbx*/,
                                GpuArray<Real, AMREX_SPACEDIM> const& problo,
                                GpuArray<Real, AMREX_SPACEDIM> const& /*probhi*/)
 {
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         Real x = problo[0] + Real(i+0.5)*dx[0];
         Real y = problo[1] + Real(j+0.5)*dx[1];
@@ -379,7 +379,7 @@ void incflo::init_flow_in_box (Box const& vbx, Box const& /*gbx*/,
 #if (AMREX_SPACEDIM == 3)
     if (periodic_dir == 0)
     {
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             Real y = Real(j+0.5)*dx[1]*(yhi-ylo) + ylo;
             Real z = Real(k+0.5)*dx[2]*(zhi-zlo) + zlo;
@@ -394,7 +394,7 @@ void incflo::init_flow_in_box (Box const& vbx, Box const& /*gbx*/,
         });
     } else if (periodic_dir == 1)
     {
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             Real x = Real(i+0.5)*dx[0]*(xhi-xlo) + xlo;
             Real z = Real(k+0.5)*dx[2]*(zhi-zlo) + zlo;
@@ -410,7 +410,7 @@ void incflo::init_flow_in_box (Box const& vbx, Box const& /*gbx*/,
     } else if (periodic_dir == 2)
 #endif
     {
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             Real x = Real(i+0.5)*dx[0]*(xhi-xlo) + xlo;
             Real y = Real(j+0.5)*dx[1]*(yhi-ylo) + ylo;
@@ -443,7 +443,7 @@ void incflo::init_circ_traceradvect (Box const& vbx, Box const& /*gbx*/,
 {
 
 #if (AMREX_SPACEDIM == 2)
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         vel(i,j,k,0) = 1.;
         vel(i,j,k,1) = 0.5;
@@ -472,7 +472,7 @@ void incflo::init_circ_traceradvect (Box const& vbx, Box const& /*gbx*/,
     });
 
 #elif (AMREX_SPACEDIM == 3)
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         vel(i,j,k,0) = 1.;
         vel(i,j,k,1) = 0.5;
@@ -515,7 +515,7 @@ void incflo::init_circ_traceradvect (Box const& vbx, Box const& /*gbx*/,
                                           GpuArray<Real, AMREX_SPACEDIM> const& /*probhi*/)
 {
 
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         Real x = (i+0.5)*dx[0];
         Real y = (j+0.5)*dx[1];
@@ -550,7 +550,7 @@ void incflo::init_couette (Box const& vbx, Box const& /*gbx*/,
                            GpuArray<Real, AMREX_SPACEDIM> const& /*probhi*/)
 {
     Real num_cells_y = static_cast<Real>(domain.length(1));
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         Real y = Real(j+0.5) / num_cells_y;
         AMREX_D_TERM(vel(i,j,k,0) *= (y-Real(0.5));,
@@ -575,7 +575,7 @@ void incflo::init_channel_slant (Box const& vbx, Box const& /*gbx*/,
     ParmParse pp("cylinder");
     pp.get("direction",  direction);
 
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         if (density(i,j,k)>0) {
 
@@ -623,7 +623,7 @@ void incflo::init_rayleigh_taylor (Box const& vbx, Box const& /*gbx*/,
     const Real L_x    = probhi[0] - problo[0];
 
 #if (AMREX_SPACEDIM == 2)
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         vel(i,j,k,0) = Real(0.0);
         vel(i,j,k,1) = Real(0.0);
@@ -641,7 +641,7 @@ void incflo::init_rayleigh_taylor (Box const& vbx, Box const& /*gbx*/,
 #elif (AMREX_SPACEDIM == 3)
 
     const Real splity = Real(0.5)*(problo[1] + probhi[1]);
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         vel(i,j,k,0) = Real(0.0);
         vel(i,j,k,1) = Real(0.0);
@@ -671,7 +671,7 @@ void incflo::init_tuscan (Box const& vbx, Box const& /*gbx*/,
 
 {
     int half_num_cells = domain.length(AMREX_SPACEDIM-1) / 2;
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         AMREX_D_TERM(vel(i,j,k,0) = Real(0.0);,
                      vel(i,j,k,1) = Real(0.0);,
@@ -703,7 +703,7 @@ void incflo::init_jump (Box const& vbx, Box const& /*gbx*/,
     }
 
     int half_num_cells = domain.length(direction) / 2;
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         if (direction == 0) {
             if (i <= half_num_cells) {
@@ -740,7 +740,7 @@ void incflo::init_boussinesq_bubble (Box const& vbx, Box const& /*gbx*/,
 
     if (111 == m_probtype)
     {
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             vel(i,j,k,0) = 0.0;
             vel(i,j,k,1) = 0.0;
@@ -765,7 +765,7 @@ void incflo::init_boussinesq_bubble (Box const& vbx, Box const& /*gbx*/,
     }
 #if (AMREX_SPACEDIM == 3)
     else if (112 == m_probtype) {
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             vel(i,j,k,0) = Real(0.0);
             vel(i,j,k,1) = Real(0.0);
@@ -785,7 +785,7 @@ void incflo::init_boussinesq_bubble (Box const& vbx, Box const& /*gbx*/,
                 tracer(i,j,k,0) = Real(0.01);
         });
     } else if (113 == m_probtype) {
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             vel(i,j,k,0) = 0.0;
             vel(i,j,k,1) = 0.0;
@@ -821,7 +821,7 @@ void incflo::init_periodic_tracer (Box const& vbx, Box const& /*gbx*/,
 
     if (m_probtype == 12)
     {
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             constexpr Real A = Real(1.0);
             Real x = Real(i+0.5)*dx[0];
@@ -840,7 +840,7 @@ void incflo::init_periodic_tracer (Box const& vbx, Box const& /*gbx*/,
     else if (m_probtype == 122)
     {
         Real B = Real(0.1);
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             constexpr Real A = Real(1.0);
             Real x = Real(i+0.5)*dx[0];
@@ -876,7 +876,7 @@ void incflo::init_double_shear_layer (Box const& vbx, Box const& /*gbx*/,
     {
         constexpr Real m_fourth = Real(0.25);
         constexpr Real m_half   = Real(0.5);
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             Real x = Real(i+0.5) * dx[0];
             Real y = Real(j+0.5) * dx[1];
@@ -896,7 +896,7 @@ void incflo::init_double_shear_layer (Box const& vbx, Box const& /*gbx*/,
     else if (22 == m_probtype)
     {
         constexpr Real m_half = Real(0.5);
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             Real y = Real(j+0.5) * dx[1];
             Real z = Real(k+0.5) * dx[2];
@@ -914,7 +914,7 @@ void incflo::init_double_shear_layer (Box const& vbx, Box const& /*gbx*/,
     else if (23 == m_probtype)
     {
         constexpr Real m_half = Real(0.5);
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             Real x = Real(i+0.5) * dx[0];
             Real z = Real(k+0.5) * dx[2];
@@ -957,7 +957,7 @@ void incflo::init_plane_poiseuille (Box const& vbx, Box const& /*gbx*/,
     if (31 == m_probtype)
     {
         Real u = m_ic_u;
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             Real y = Real(j+0.5)*dyinv;
             AMREX_D_TERM(vel(i,j,k,0) = Real(6.0) * u * y * (Real(1.0)-y);,
@@ -976,7 +976,7 @@ void incflo::init_plane_poiseuille (Box const& vbx, Box const& /*gbx*/,
     else if (311 == m_probtype)
     {
         Real u = m_ic_u;
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             Real z = Real(k+0.5)*dzinv;
             AMREX_D_TERM(vel(i,j,k,0) = Real(6.0) * u * z * (Real(1.0)-z);,
@@ -994,7 +994,7 @@ void incflo::init_plane_poiseuille (Box const& vbx, Box const& /*gbx*/,
     }
     else if (41 == m_probtype)
     {
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             Real z = Real(k+0.5)*dzinv;
             AMREX_D_TERM(vel(i,j,k,0) = Real(0.5)*z;,
@@ -1013,7 +1013,7 @@ void incflo::init_plane_poiseuille (Box const& vbx, Box const& /*gbx*/,
     else if (32 == m_probtype)
     {
         Real v = m_ic_v;
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             Real z = Real(k+0.5)*dzinv;
             AMREX_D_TERM(vel(i,j,k,0) = Real(0.0);,
@@ -1032,7 +1032,7 @@ void incflo::init_plane_poiseuille (Box const& vbx, Box const& /*gbx*/,
     else if (322 == m_probtype)
     {
         Real v = m_ic_v;
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
             Real x = Real(i+0.5)*dxinv;
             AMREX_D_TERM(vel(i,j,k,0) = Real(0.0);,
@@ -1053,7 +1053,7 @@ void incflo::init_plane_poiseuille (Box const& vbx, Box const& /*gbx*/,
 #if (AMREX_SPACEDIM == 3)
         Real w = m_ic_w;
 #endif
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
 #if (AMREX_SPACEDIM == 3)
             Real x = Real(i+0.5)*dxinv;
@@ -1076,7 +1076,7 @@ void incflo::init_plane_poiseuille (Box const& vbx, Box const& /*gbx*/,
 #if (AMREX_SPACEDIM == 3)
         Real w = m_ic_w;
 #endif
-        amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+        ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
 #if (AMREX_SPACEDIM == 3)
             Real y = Real(j+0.5)*dyinv;
@@ -1109,7 +1109,7 @@ void incflo::init_burggraf (Box const& vbx, Box const& /*gbx*/,
                             GpuArray<Real, AMREX_SPACEDIM> const& /*problo*/,
                             GpuArray<Real, AMREX_SPACEDIM> const& /*probhi*/)
 {
-    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         Real x = Real(i+0.5)*dx[0];
         Real y = Real(j+0.5)*dx[1];
diff --git a/src/projection/incflo_apply_nodal_projection.cpp b/src/projection/incflo_apply_nodal_projection.cpp
index 1fb8ba93c..3759a8a58 100644
--- a/src/projection/incflo_apply_nodal_projection.cpp
+++ b/src/projection/incflo_apply_nodal_projection.cpp
@@ -49,7 +49,7 @@ void incflo::ApplyNodalProjection (Vector<MultiFab const*> density,
                 Array4<Real> const& u = ld.velocity.array(mfi);
                 Array4<Real const> const& rho = density[lev]->const_array(mfi);
                 Array4<Real const> const& gp = ld.gp.const_array(mfi);
-                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     Real soverrho = scaling_factor / rho(i,j,k);
                     AMREX_D_TERM(u(i,j,k,0) += gp(i,j,k,0) * soverrho;,
@@ -111,7 +111,7 @@ void incflo::ApplyNodalProjection (Vector<MultiFab const*> const& density,
                 Box const& bx = mfi.tilebox();
                 Array4<Real> const& sig = sigma[lev].array(mfi);
                 Array4<Real const> const& rho = density[lev]->const_array(mfi);
-                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     sig(i,j,k) = scaling_factor / rho(i,j,k);
                 });
@@ -211,22 +211,22 @@ void incflo::ApplyNodalProjection (Vector<MultiFab const*> const& density,
             Array4<Real const> const& gp_proj = gradphi[lev]->const_array(mfi);
             Array4<Real const> const& p_proj = phi[lev]->const_array(mfi);
             if (incremental) {
-                amrex::ParallelFor(tbx, AMREX_SPACEDIM,
+                ParallelFor(tbx, AMREX_SPACEDIM,
                 [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                 {
                     gp_lev(i,j,k,n) += gp_proj(i,j,k,n);
                 });
-                amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     p_lev (i,j,k) += p_proj(i,j,k);
                 });
             } else {
-                amrex::ParallelFor(tbx, AMREX_SPACEDIM,
+                ParallelFor(tbx, AMREX_SPACEDIM,
                 [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
                 {
                     gp_lev(i,j,k,n) = gp_proj(i,j,k,n);
                 });
-                amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     p_lev(i,j,k) = p_proj(i,j,k);
                 });
diff --git a/src/rheology/incflo_rheology.cpp b/src/rheology/incflo_rheology.cpp
index 2fb2586ec..ca2d4b231 100644
--- a/src/rheology/incflo_rheology.cpp
+++ b/src/rheology/incflo_rheology.cpp
@@ -107,7 +107,7 @@ void incflo::compute_viscosity_at_level (int /*lev*/,
                 auto typ = flag_fab.getType(bx);
                 if (typ == FabType::covered)
                 {
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         eta_arr(i,j,k) = Real(0.0);
                     });
@@ -115,7 +115,7 @@ void incflo::compute_viscosity_at_level (int /*lev*/,
                 else if (typ == FabType::singlevalued)
                 {
                     auto const& flag_arr = flag_fab.const_array();
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         Real sr = incflo_strainrate_eb(i,j,k,AMREX_D_DECL(idx,idy,idz),vel_arr,flag_arr(i,j,k));
                         eta_arr(i,j,k) = non_newtonian_viscosity(sr);
@@ -124,7 +124,7 @@ void incflo::compute_viscosity_at_level (int /*lev*/,
                 else
 #endif
                 {
-                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                    ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                     {
                         Real sr = incflo_strainrate(i,j,k,AMREX_D_DECL(idx,idy,idz),vel_arr);
                         eta_arr(i,j,k) = non_newtonian_viscosity(sr);
diff --git a/src/setup/init.cpp b/src/setup/init.cpp
index 519477ca9..5ca9f2310 100644
--- a/src/setup/init.cpp
+++ b/src/setup/init.cpp
@@ -384,10 +384,31 @@ void incflo::InitialPressureProjection()
     for (int lev = 0; lev <= finest_level; ++lev) {
         vel[lev].define(grids[lev], dmap[lev], AMREX_SPACEDIM, nGhost,
                         MFInfo(), *m_factory[lev]);
+
         for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
             vel[lev].setVal(m_gravity[idim], idim, 1, 1);
         }
-    }
+
+        auto& ld = *m_leveldata[lev];
+
+        Real rho0 = m_ro_0;
+        if (rho0 > 0.0) {
+            for (MFIter mfi(ld.density,TilingIfNotGPU()); mfi.isValid(); ++mfi)
+            {
+                Box const& bx = mfi.growntilebox(1);
+                Array4<Real const> const& rho_arr = ld.density.const_array(mfi);
+                Array4<Real      > const& vel_arr = vel[lev].array(mfi);
+
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                {
+                    Real rhofac = (rho_arr(i,j,k) - rho0) / rho_arr(i,j,k);
+                    AMREX_D_TERM(vel_arr(i,j,k,0) *= rhofac;,
+                                 vel_arr(i,j,k,1) *= rhofac;,
+                                 vel_arr(i,j,k,2) *= rhofac;);
+                });
+            } // mfi
+        } // rho0
+    } // lev
 
     // Cell-centered divergence condition source term
     // Always zero this here
diff --git a/src/utilities/io.cpp b/src/utilities/io.cpp
index 2d6177abc..b7946cc5b 100644
--- a/src/utilities/io.cpp
+++ b/src/utilities/io.cpp
@@ -451,6 +451,7 @@ void incflo::WritePlotFile()
     if (m_plt_gpx) {
         for (int lev = 0; lev <= finest_level; ++lev) {
             MultiFab::Copy(mf[lev], m_leveldata[lev]->gp, 0, icomp, 1, 0);
+            mf[lev].plus(m_gp0[0],icomp,1,0);
         }
         pltscaVarsName.push_back("gpx");
         ++icomp;
@@ -458,6 +459,7 @@ void incflo::WritePlotFile()
     if (m_plt_gpy) {
         for (int lev = 0; lev <= finest_level; ++lev) {
             MultiFab::Copy(mf[lev], m_leveldata[lev]->gp, 1, icomp, 1, 0);
+            mf[lev].plus(m_gp0[1],icomp,1,0);
         }
         pltscaVarsName.push_back("gpy");
         ++icomp;
@@ -466,6 +468,7 @@ void incflo::WritePlotFile()
     if (m_plt_gpz) {
         for (int lev = 0; lev <= finest_level; ++lev) {
             MultiFab::Copy(mf[lev], m_leveldata[lev]->gp, 2, icomp, 1, 0);
+            mf[lev].plus(m_gp0[2],icomp,1,0);
         }
         pltscaVarsName.push_back("gpz");
         ++icomp;