diff --git a/Source/Maestro.H b/Source/Maestro.H
index 67a76550..da01c727 100644
--- a/Source/Maestro.H
+++ b/Source/Maestro.H
@@ -1734,7 +1734,7 @@ class Maestro : public amrex::AmrCore {
                           amrex::Array4<amrex::Real> const ury,
                           amrex::Array4<amrex::Real> const uimhy,
                           const amrex::Box& domainBox,
-                          const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx);
+                          const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx) const;
 
     void VelPredVelocities(const amrex::MFIter& mfi,
                            amrex::Array4<const amrex::Real> const utilde,
@@ -1755,7 +1755,8 @@ class Maestro : public amrex::AmrCore {
                            amrex::Array4<const amrex::Real> const force,
                            amrex::Array4<const amrex::Real> const w0_cart_in,
                            const amrex::Box& domainBox,
-                           const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx);
+                           const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx) const;
+
 #else
     void VelPredInterface(const amrex::MFIter& mfi,
                           amrex::Array4<const amrex::Real> const utilde,
@@ -1779,7 +1780,7 @@ class Maestro : public amrex::AmrCore {
                           amrex::Array4<amrex::Real> const urz,
                           amrex::Array4<amrex::Real> const uimhz,
                           const amrex::Box& domainBox,
-                          const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx);
+                          const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx) const;
 
     void VelPredTransverse(const amrex::MFIter& mfi,
                            amrex::Array4<const amrex::Real> const utilde,
@@ -1802,7 +1803,7 @@ class Maestro : public amrex::AmrCore {
                            amrex::Array4<amrex::Real> const wimhxy,
                            amrex::Array4<amrex::Real> const wimhyx,
                            const amrex::Box& domainBox,
-                           const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx);
+                           const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx) const;
 
     void VelPredVelocities(const amrex::MFIter& mfi,
                            amrex::Array4<const amrex::Real> const utilde,
@@ -1836,7 +1837,8 @@ class Maestro : public amrex::AmrCore {
                            amrex::Array4<const amrex::Real> const force,
                            amrex::Array4<const amrex::Real> const w0_cart_in,
                            const amrex::Box& domainBox,
-                           const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx);
+                           const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx) const;
+
 #endif
     ////////////////////////
 
@@ -1912,8 +1914,7 @@ class Maestro : public amrex::AmrCore {
     amrex::Vector<amrex::iMultiFab> cell_cc_to_r;
 
     /// stores domain boundary conditions.
-    /// These muse be vectors (rather than arrays) so we can ParmParse them
-    IntVector phys_bc;
+    amrex::BCRec phys_bc;
 
     /// Boundary condition objects needed for FillPatch routines.
     /// This is essentially an array (over components)
diff --git a/Source/MaestroFillData.cpp b/Source/MaestroFillData.cpp
index 01e461a0..77ab320c 100644
--- a/Source/MaestroFillData.cpp
+++ b/Source/MaestroFillData.cpp
@@ -194,8 +194,6 @@ void Maestro::FillUmacGhost(
         const auto domlo = domainBox.loVect3d();
         const auto domhi = domainBox.hiVect3d();
 
-        const int* AMREX_RESTRICT physbc_p = phys_bc.dataPtr();
-
         // get references to the MultiFabs at level lev
         MultiFab& sold_mf =
             sold[lev];  // need a cell-centered MF for the MFIter
@@ -216,11 +214,13 @@ void Maestro::FillUmacGhost(
 #if (AMREX_SPACEDIM == 3)
             const Array4<Real> wmac = umac_in[lev][2].array(mfi);
 #endif
+            int bclo = phys_bc.lo(0);
+            int bchi = phys_bc.hi(0);
 
             ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                 // lo x-faces
                 if (i == domlo[0] - 1) {
-                    switch (physbc_p[0]) {  // NOLINT(bugprone-switch-missing-default-case)
+                    switch (bclo) {  // NOLINT(bugprone-switch-missing-default-case)
                         case amrex::PhysBCType::inflow:
                             umac(i, j, k) = umac(i + 1, j, k);
                             vmac(i, j, k) = 0.0;
@@ -258,7 +258,7 @@ void Maestro::FillUmacGhost(
 
                 // hi x-faces
                 if (i == domhi[0] + 2) {
-                    switch (physbc_p[AMREX_SPACEDIM]) {  // NOLINT(bugprone-switch-missing-default-case)
+                    switch (bchi) {  // NOLINT(bugprone-switch-missing-default-case)
                         case amrex::PhysBCType::inflow:
                             umac(i, j, k) = umac(i - 1, j, k);
                             vmac(i - 1, j, k) = 0.0;
@@ -297,10 +297,13 @@ void Maestro::FillUmacGhost(
 
             Gpu::synchronize();
 
+            bclo = phys_bc.lo(1);
+            bchi = phys_bc.hi(1);
+
             ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                 // lo y-faces
                 if (j == domlo[1] - 1) {
-                    switch (physbc_p[1]) {  // NOLINT(bugprone-switch-missing-default-case)
+                    switch (bclo) {  // NOLINT(bugprone-switch-missing-default-case)
                         case amrex::PhysBCType::inflow:
                             umac(i, j, k) = 0.0;
                             vmac(i, j, k) = vmac(i, j + 1, k);
@@ -338,7 +341,7 @@ void Maestro::FillUmacGhost(
 
                 // hi y-faces
                 if (j == domhi[1] + 2) {
-                    switch (physbc_p[AMREX_SPACEDIM + 1]) {  // NOLINT(bugprone-switch-missing-default-case)
+                    switch (bchi) {  // NOLINT(bugprone-switch-missing-default-case)
                         case amrex::PhysBCType::inflow:
                             umac(i, j - 1, k) = 0.0;
                             vmac(i, j, k) = vmac(i, j - 1, k);
@@ -379,10 +382,13 @@ void Maestro::FillUmacGhost(
 
             Gpu::synchronize();
 
+            bclo = phys_bc.lo(2);
+            bchi = phys_bc.hi(2);
+
             ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                 // lo z-faces
                 if (k == domlo[2] - 1) {
-                    switch (physbc_p[2]) {  // NOLINT(bugprone-switch-missing-default-case)
+                    switch (bclo) {  // NOLINT(bugprone-switch-missing-default-case)
                         case amrex::PhysBCType::inflow:
                             umac(i, j, k) = 0.0;
                             vmac(i, j, k) = 0.0;
@@ -412,7 +418,7 @@ void Maestro::FillUmacGhost(
 
                 // hi z-faces
                 if (k == domhi[2] + 2) {
-                    switch (physbc_p[2 + AMREX_SPACEDIM]) {  // NOLINT(bugprone-switch-missing-default-case)
+                    switch (bchi) {  // NOLINT(bugprone-switch-missing-default-case)
                         case amrex::PhysBCType::inflow:
                             umac(i, j, k - 1) = 0.0;
                             vmac(i, j, k - 1) = 0.0;
diff --git a/Source/MaestroMacProj.cpp b/Source/MaestroMacProj.cpp
index 798b6ece..bfa8b162 100644
--- a/Source/MaestroMacProj.cpp
+++ b/Source/MaestroMacProj.cpp
@@ -434,14 +434,14 @@ void Maestro::SetMacSolverBCs(MLABecLaplacian& mlabec) {
             mlmg_lobc[idim] = mlmg_hibc[idim] = LinOpBCType::Periodic;
         } else {
             // lo-side BCs
-            if (phys_bc[idim] == amrex::PhysBCType::outflow) {
+            if (phys_bc.lo(idim) == amrex::PhysBCType::outflow) {
                 mlmg_lobc[idim] = LinOpBCType::Dirichlet;
             } else {
                 mlmg_lobc[idim] = LinOpBCType::Neumann;
             }
 
             // hi-side BCs
-            if (phys_bc[AMREX_SPACEDIM + idim] == amrex::PhysBCType::outflow) {
+            if (phys_bc.hi(idim) == amrex::PhysBCType::outflow) {
                 mlmg_hibc[idim] = LinOpBCType::Dirichlet;
             } else {
                 mlmg_hibc[idim] = LinOpBCType::Neumann;
diff --git a/Source/MaestroMakeUtrans.cpp b/Source/MaestroMakeUtrans.cpp
index 255f85f7..a3b5ed95 100644
--- a/Source/MaestroMakeUtrans.cpp
+++ b/Source/MaestroMakeUtrans.cpp
@@ -82,8 +82,8 @@ void Maestro::MakeUtrans(
             }
 
             // create utrans
-            int bclo = phys_bc[0];
-            int bchi = phys_bc[AMREX_SPACEDIM];
+            int bclo = phys_bc.lo(0);
+            int bchi = phys_bc.hi(0);
 
             ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                 Real ulx = 0.0;
@@ -182,8 +182,8 @@ void Maestro::MakeUtrans(
             }
 
             // create vtrans
-            int bclo = phys_bc[1];
-            int bchi = phys_bc[AMREX_SPACEDIM + 1];
+            int bclo = phys_bc.lo(1);
+            int bchi = phys_bc.hi(1);
 
             ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                 Real vly = 0.0;
@@ -298,8 +298,8 @@ void Maestro::MakeUtrans(
             }
 
             // create utrans
-            int bclo = phys_bc[0];
-            int bchi = phys_bc[AMREX_SPACEDIM];
+            int bclo = phys_bc.lo(0);
+            int bchi = phys_bc.hi(0);
 
             ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                 Real ulx = 0.0;
@@ -413,8 +413,8 @@ void Maestro::MakeUtrans(
             }
 
             // create vtrans
-            int bclo = phys_bc[1];
-            int bchi = phys_bc[AMREX_SPACEDIM + 1];
+            int bclo = phys_bc.lo(1);
+            int bchi = phys_bc.hi(1);
 
             ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                 Real vly = 0.0;
@@ -529,8 +529,8 @@ void Maestro::MakeUtrans(
             }
 
             // create wtrans
-            int bclo = phys_bc[2];
-            int bchi = phys_bc[AMREX_SPACEDIM + 2];
+            int bclo = phys_bc.lo(2);
+            int bchi = phys_bc.hi(2);
 
             ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                 Real wlz = 0.0;
diff --git a/Source/MaestroNodalProj.cpp b/Source/MaestroNodalProj.cpp
index a979827c..57380224 100644
--- a/Source/MaestroNodalProj.cpp
+++ b/Source/MaestroNodalProj.cpp
@@ -155,13 +155,13 @@ void Maestro::NodalProj(int proj_type, Vector<MultiFab>& rhcc,
         if (Geom(0).isPeriodic(idim)) {
             mlmg_lobc[idim] = mlmg_hibc[idim] = LinOpBCType::Periodic;
         } else {
-            if (phys_bc[idim] == amrex::PhysBCType::outflow) {
+            if (phys_bc.lo(idim) == amrex::PhysBCType::outflow) {
                 mlmg_lobc[idim] = LinOpBCType::Dirichlet;
             } else {
                 mlmg_lobc[idim] = LinOpBCType::Neumann;
             }
 
-            if (phys_bc[AMREX_SPACEDIM + idim] == amrex::PhysBCType::outflow) {
+            if (phys_bc.hi(idim) == amrex::PhysBCType::outflow) {
                 mlmg_hibc[idim] = LinOpBCType::Dirichlet;
             } else {
                 mlmg_hibc[idim] = LinOpBCType::Neumann;
@@ -505,8 +505,8 @@ void Maestro::SetBoundaryVelocity(Vector<MultiFab>& vel) {
         const Box& domainBox = geom[lev].Domain();
 
         for (int idir = 0; idir < BL_SPACEDIM; idir++) {
-            if (phys_bc[idir] != amrex::PhysBCType::inflow &&
-                phys_bc[AMREX_SPACEDIM + idir] != amrex::PhysBCType::inflow) {
+            if (phys_bc.lo(idir) != amrex::PhysBCType::inflow &&
+                phys_bc.hi(idir) != amrex::PhysBCType::inflow) {
                 vel[lev].setBndry(0.0, idir, 1);
             } else {
 #ifdef _OPENMP
@@ -521,7 +521,7 @@ void Maestro::SetBoundaryVelocity(Vector<MultiFab>& vel) {
 
                     BoxList bxlist(reg);
 
-                    if (phys_bc[idir] == amrex::PhysBCType::inflow &&
+                    if (phys_bc.lo(idir) == amrex::PhysBCType::inflow &&
                         reg.smallEnd(idir) == domainBox.smallEnd(idir)) {
                         Box bx;  // bx is the region we *protect* from zero'ing
 
@@ -545,7 +545,7 @@ void Maestro::SetBoundaryVelocity(Vector<MultiFab>& vel) {
                         bxlist.push_back(bx);
                     }
 
-                    if (phys_bc[AMREX_SPACEDIM + idir] == amrex::PhysBCType::inflow &&
+                    if (phys_bc.hi(idir) == amrex::PhysBCType::inflow &&
                         reg.bigEnd(idir) == domainBox.bigEnd(idir)) {
                         Box bx;  // bx is the region we *protect* from zero'ing
 
diff --git a/Source/MaestroPlot.cpp b/Source/MaestroPlot.cpp
index 68f356c2..8d66ce53 100644
--- a/Source/MaestroPlot.cpp
+++ b/Source/MaestroPlot.cpp
@@ -1681,45 +1681,42 @@ void Maestro::MakeVorticity(const Vector<MultiFab>& vel,
 
             Array4<const Real> const u = vel[lev].array(mfi);
             Array4<Real> const vort = vorticity[lev].array(mfi);
-            GpuArray<int, AMREX_SPACEDIM * 2> physbc;
-            for (int n = 0; n < AMREX_SPACEDIM * 2; ++n) {
-                physbc[n] = phys_bc[n];
-            }
+
 
 #if (AMREX_SPACEDIM == 2)
+            bool fix_lo_x = (phys_bc.lo(0) == amrex::PhysBCType::inflow || phys_bc.lo(0) == amrex::PhysBCType::slipwall ||
+                                 phys_bc.lo(0) == amrex::PhysBCType::noslipwall);
+            bool fix_hi_x = (phys_bc.hi(0) == amrex::PhysBCType::inflow || phys_bc.hi(0) == amrex::PhysBCType::slipwall ||
+                             phys_bc.hi(0) == amrex::PhysBCType::noslipwall);
+            bool fix_lo_y = (phys_bc.lo(1) == amrex::PhysBCType::inflow || phys_bc.lo(1) == amrex::PhysBCType::slipwall ||
+                                 phys_bc.lo(1) == amrex::PhysBCType::noslipwall);
+            bool fix_hi_y = (phys_bc.hi(1) == amrex::PhysBCType::inflow || phys_bc.hi(1) == amrex::PhysBCType::slipwall ||
+                             phys_bc.hi(1) == amrex::PhysBCType::noslipwall);
 
             ParallelFor(tileBox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                 Real vx = 0.5 * (u(i + 1, j, k, 1) - u(i - 1, j, k, 1)) / hx;
                 Real uy = 0.5 * (u(i, j + 1, k, 0) - u(i, j - 1, k, 0)) / hy;
 
-                if (i == ilo && (physbc[0] == amrex::PhysBCType::inflow || physbc[0] == amrex::PhysBCType::slipwall ||
-                                 physbc[0] == amrex::PhysBCType::noslipwall)) {
+                if (i == ilo && fix_lo_x) {
                     vx = (u(i + 1, j, k, 1) + 3.0 * u(i, j, k, 1) -
                           4.0 * u(i - 1, j, k, 1)) /
                          hx;
                     uy = 0.5 * (u(i, j + 1, k, 0) - u(i, j - 1, k, 0)) / hy;
 
-                } else if (i == ihi + 1 &&
-                           (physbc[AMREX_SPACEDIM] == amrex::PhysBCType::inflow ||
-                            physbc[AMREX_SPACEDIM] == amrex::PhysBCType::slipwall ||
-                            physbc[AMREX_SPACEDIM] == amrex::PhysBCType::noslipwall)) {
+                } else if (i == ihi + 1 && fix_hi_x) {
                     vx = -(u(i - 1, j, k, 1) + 3.0 * u(i, j, k, 1) -
                            4.0 * u(i + 1, j, k, 1)) /
                          hx;
                     uy = 0.5 * (u(i, j + 1, k, 0) - u(i, j - 1, k, 0)) / hy;
                 }
 
-                if (j == jlo && (physbc[1] == amrex::PhysBCType::inflow || physbc[1] == amrex::PhysBCType::slipwall ||
-                                 physbc[1] == amrex::PhysBCType::noslipwall)) {
+                if (j == jlo && fix_lo_y) {
                     vx = 0.5 * (u(i + 1, j, k, 1) - u(i - 1, j, k, 0)) / hx;
                     uy = (u(i, j + 1, k, 0) + 3.0 * u(i, j, k, 0) -
                           4.0 * u(i, j - 1, k, 0)) /
                          hy;
 
-                } else if (j == jhi + 1 &&
-                           (physbc[AMREX_SPACEDIM + 1] == amrex::PhysBCType::inflow ||
-                            physbc[AMREX_SPACEDIM + 1] == amrex::PhysBCType::slipwall ||
-                            physbc[AMREX_SPACEDIM + 1] == amrex::PhysBCType::noslipwall)) {
+                } else if (j == jhi + 1 && fix_hi_y) {
                     vx = 0.5 * (u(i + 1, j, k, 1) - u(i - 1, j, k, 1)) / hx;
                     uy = -(u(i, j - 1, k, 0) + 3.0 * u(i, j, k, 0) -
                            4.0 * u(i, j + 1, k, 0)) /
@@ -1729,7 +1726,22 @@ void Maestro::MakeVorticity(const Vector<MultiFab>& vel,
                 vort(i, j, k) = vx - uy;
             });
 
-#else
+# else
+            bool fix_lo_x =
+                (phys_bc.lo(0) == amrex::PhysBCType::inflow || phys_bc.lo(0) == amrex::PhysBCType::noslipwall);
+            bool fix_hi_x = (phys_bc.hi(0) == amrex::PhysBCType::inflow ||
+                             phys_bc.hi(0) == amrex::PhysBCType::noslipwall);
+
+            bool fix_lo_y =
+                (phys_bc.lo(1) == amrex::PhysBCType::inflow || phys_bc.lo(1) == amrex::PhysBCType::noslipwall);
+            bool fix_hi_y = (phys_bc.hi(1) == amrex::PhysBCType::inflow ||
+                             phys_bc.hi(1) == amrex::PhysBCType::noslipwall);
+
+            bool fix_lo_z =
+                (phys_bc.lo(2) == amrex::PhysBCType::inflow || phys_bc.lo(2) == amrex::PhysBCType::noslipwall);
+            bool fix_hi_z = (phys_bc.hi(2) == amrex::PhysBCType::inflow ||
+                             phys_bc.hi(2) == amrex::PhysBCType::noslipwall);
+
             ParallelFor(tileBox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                 Real uy = 0.5 * (u(i, j + 1, k, 0) - u(i, j - 1, k, 0)) / hy;
                 Real uz = 0.5 * (u(i, j, k + 1, 0) - u(i, j, k - 1, 0)) / hz;
@@ -1738,21 +1750,6 @@ void Maestro::MakeVorticity(const Vector<MultiFab>& vel,
                 Real wx = 0.5 * (u(i + 1, j, k, 2) - u(i - 1, j, k, 2)) / hx;
                 Real wy = 0.5 * (u(i, j + 1, k, 2) - u(i, j - 1, k, 2)) / hy;
 
-                bool fix_lo_x =
-                    (physbc[0] == amrex::PhysBCType::inflow || physbc[0] == amrex::PhysBCType::noslipwall);
-                bool fix_hi_x = (physbc[AMREX_SPACEDIM] == amrex::PhysBCType::inflow ||
-                                 physbc[AMREX_SPACEDIM] == amrex::PhysBCType::noslipwall);
-
-                bool fix_lo_y =
-                    (physbc[1] == amrex::PhysBCType::inflow || physbc[1] == amrex::PhysBCType::noslipwall);
-                bool fix_hi_y = (physbc[AMREX_SPACEDIM + 1] == amrex::PhysBCType::inflow ||
-                                 physbc[AMREX_SPACEDIM + 1] == amrex::PhysBCType::noslipwall);
-
-                bool fix_lo_z =
-                    (physbc[2] == amrex::PhysBCType::inflow || physbc[2] == amrex::PhysBCType::noslipwall);
-                bool fix_hi_z = (physbc[AMREX_SPACEDIM + 2] == amrex::PhysBCType::inflow ||
-                                 physbc[AMREX_SPACEDIM + 2] == amrex::PhysBCType::noslipwall);
-
                 // First do all the faces
                 if (fix_lo_x && i == ilo) {
                     vx = (u(i + 1, j, k, 1) + 3.0 * u(i, j, k, 1) -
diff --git a/Source/MaestroSetup.cpp b/Source/MaestroSetup.cpp
index 892053d2..e3c9aca8 100644
--- a/Source/MaestroSetup.cpp
+++ b/Source/MaestroSetup.cpp
@@ -240,10 +240,9 @@ void Maestro::ReadParameters() {
     // storey boundary conditions in a single array
     // order shall be
     // LO_X, LO_Y, (LO_Z), HI_X, HI_Y, (HI_Z)
-    phys_bc.resize(2 * AMREX_SPACEDIM);
     for (int i = 0; i < AMREX_SPACEDIM; ++i) {
-        phys_bc[i] = lo_bc[i];
-        phys_bc[i + AMREX_SPACEDIM] = hi_bc[i];
+        phys_bc.setLo(i, lo_bc[i]);
+        phys_bc.setHi(i, hi_bc[i]);
     }
 }
 
@@ -293,13 +292,13 @@ void Maestro::BCSetup() {
         //
         for (int dir = 0; dir < AMREX_SPACEDIM; dir++) {
             if (Geom(0).isPeriodic(dir)) {
-                if (phys_bc[dir] != amrex::PhysBCType::interior) {
+                if (phys_bc.lo(dir) != amrex::PhysBCType::interior) {
                     std::cerr
                         << "Maestro::ReadParameters:periodic in direction "
                         << dir << " but low BC is not amrex::PhysBCType::interior\n";
                     Error();
                 }
-                if (phys_bc[AMREX_SPACEDIM + dir] != amrex::PhysBCType::interior) {
+                if (phys_bc.hi(dir) != amrex::PhysBCType::interior) {
                     std::cerr
                         << "Maestro::ReadParameters:periodic in direction "
                         << dir << " but high BC is not amrex::PhysBCType::interior\n";
@@ -312,12 +311,12 @@ void Maestro::BCSetup() {
         // Do idiot check.  If not periodic, should be no interior.
         //
         for (int dir = 0; dir < AMREX_SPACEDIM; dir++) {
-            if (phys_bc[dir] == amrex::PhysBCType::interior) {
+            if (phys_bc.lo(dir) == amrex::PhysBCType::interior) {
                 std::cerr << "Maestro::ReadParameters:interior bc in direction "
                           << dir << " but not periodic\n";
                 Error();
             }
-            if (phys_bc[AMREX_SPACEDIM + dir] == amrex::PhysBCType::interior) {
+            if (phys_bc.hi(dir) == amrex::PhysBCType::interior) {
                 std::cerr << "Maestro::ReadParameters:interior bc in direction "
                           << dir << " but not periodic\n";
                 Error();
@@ -328,7 +327,7 @@ void Maestro::BCSetup() {
     // set up boundary conditions for Fillpatch operations
     for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
         // lo-side bcs
-        if (phys_bc[dir] == amrex::PhysBCType::interior) {
+        if (phys_bc.lo(dir) == amrex::PhysBCType::interior) {
             // periodic uses "internal Dirichlet"
             for (int comp = 0; comp < AMREX_SPACEDIM; ++comp) {
                 bcs_u[comp].setLo(dir, BCType::int_dir);
@@ -346,7 +345,7 @@ void Maestro::BCSetup() {
             for (int comp = 0; comp < Nscal; ++comp) {
                 bcs_f[comp].setLo(dir, BCType::int_dir);
             }
-        } else if (phys_bc[dir] == amrex::PhysBCType::inflow) {
+        } else if (phys_bc.lo(dir) == amrex::PhysBCType::inflow) {
             // inflow
             for (int comp = 0; comp < AMREX_SPACEDIM; ++comp) {
                 bcs_u[comp].setLo(dir, BCType::ext_dir);
@@ -364,7 +363,7 @@ void Maestro::BCSetup() {
             for (int comp = 0; comp < Nscal; ++comp) {
                 bcs_f[comp].setLo(dir, BCType::foextrap);
             }
-        } else if (phys_bc[dir] == amrex::PhysBCType::outflow) {
+        } else if (phys_bc.lo(dir) == amrex::PhysBCType::outflow) {
             // outflow
             for (int comp = 0; comp < AMREX_SPACEDIM; ++comp) {
                 bcs_u[comp].setLo(dir, BCType::foextrap);
@@ -382,7 +381,7 @@ void Maestro::BCSetup() {
             for (int comp = 0; comp < Nscal; ++comp) {
                 bcs_f[comp].setLo(dir, BCType::foextrap);
             }
-        } else if (phys_bc[dir] == amrex::PhysBCType::symmetry) {
+        } else if (phys_bc.lo(dir) == amrex::PhysBCType::symmetry) {
             // symmetry
             for (int comp = 0; comp < AMREX_SPACEDIM; ++comp) {
                 bcs_u[comp].setLo(dir, BCType::reflect_even);
@@ -401,7 +400,7 @@ void Maestro::BCSetup() {
             for (int comp = 0; comp < Nscal; ++comp) {
                 bcs_f[comp].setLo(dir, BCType::reflect_even);
             }
-        } else if (phys_bc[dir] == amrex::PhysBCType::slipwall) {
+        } else if (phys_bc.lo(dir) == amrex::PhysBCType::slipwall) {
             // slip wall
             for (int comp = 0; comp < AMREX_SPACEDIM; ++comp) {
                 bcs_u[comp].setLo(dir, BCType::hoextrap);
@@ -420,7 +419,7 @@ void Maestro::BCSetup() {
             for (int comp = 0; comp < Nscal; ++comp) {
                 bcs_f[comp].setLo(dir, BCType::foextrap);
             }
-        } else if (phys_bc[dir] == amrex::PhysBCType::noslipwall) {
+        } else if (phys_bc.lo(dir) == amrex::PhysBCType::noslipwall) {
             // no-slip wall
             for (int comp = 0; comp < AMREX_SPACEDIM; ++comp) {
                 bcs_u[comp].setLo(dir, BCType::ext_dir);
@@ -443,7 +442,7 @@ void Maestro::BCSetup() {
         }
 
         // hi-side bcs
-        if (phys_bc[AMREX_SPACEDIM + dir] == amrex::PhysBCType::interior) {
+        if (phys_bc.hi(dir) == amrex::PhysBCType::interior) {
             // periodic uses "internal Dirichlet"
             for (int comp = 0; comp < AMREX_SPACEDIM; ++comp) {
                 bcs_u[comp].setHi(dir, BCType::int_dir);
@@ -461,7 +460,7 @@ void Maestro::BCSetup() {
             for (int comp = 0; comp < Nscal; ++comp) {
                 bcs_f[comp].setHi(dir, BCType::int_dir);
             }
-        } else if (phys_bc[AMREX_SPACEDIM + dir] == amrex::PhysBCType::inflow) {
+        } else if (phys_bc.hi(dir) == amrex::PhysBCType::inflow) {
             // inflow
             for (int comp = 0; comp < AMREX_SPACEDIM; ++comp) {
                 bcs_u[comp].setHi(dir, BCType::ext_dir);
@@ -479,7 +478,7 @@ void Maestro::BCSetup() {
             for (int comp = 0; comp < Nscal; ++comp) {
                 bcs_f[comp].setHi(dir, BCType::foextrap);
             }
-        } else if (phys_bc[AMREX_SPACEDIM + dir] == amrex::PhysBCType::outflow) {
+        } else if (phys_bc.hi(dir) == amrex::PhysBCType::outflow) {
             // outflow
             for (int comp = 0; comp < AMREX_SPACEDIM; ++comp) {
                 bcs_u[comp].setHi(dir, BCType::foextrap);
@@ -497,7 +496,7 @@ void Maestro::BCSetup() {
             for (int comp = 0; comp < Nscal; ++comp) {
                 bcs_f[comp].setHi(dir, BCType::foextrap);
             }
-        } else if (phys_bc[AMREX_SPACEDIM + dir] == amrex::PhysBCType::symmetry) {
+        } else if (phys_bc.hi(dir) == amrex::PhysBCType::symmetry) {
             // symmetry
             for (int comp = 0; comp < AMREX_SPACEDIM; ++comp) {
                 bcs_u[comp].setHi(dir, BCType::reflect_even);
@@ -516,7 +515,7 @@ void Maestro::BCSetup() {
             for (int comp = 0; comp < Nscal; ++comp) {
                 bcs_f[comp].setHi(dir, BCType::reflect_even);
             }
-        } else if (phys_bc[AMREX_SPACEDIM + dir] == amrex::PhysBCType::slipwall) {
+        } else if (phys_bc.hi(dir) == amrex::PhysBCType::slipwall) {
             // slip wall
             for (int comp = 0; comp < AMREX_SPACEDIM; ++comp) {
                 bcs_u[comp].setHi(dir, BCType::hoextrap);
@@ -535,7 +534,7 @@ void Maestro::BCSetup() {
             for (int comp = 0; comp < Nscal; ++comp) {
                 bcs_f[comp].setHi(dir, BCType::foextrap);
             }
-        } else if (phys_bc[AMREX_SPACEDIM + dir] == amrex::PhysBCType::noslipwall) {
+        } else if (phys_bc.hi(dir) == amrex::PhysBCType::noslipwall) {
             // no-slip wall
             for (int comp = 0; comp < AMREX_SPACEDIM; ++comp) {
                 bcs_u[comp].setHi(dir, BCType::ext_dir);
diff --git a/Source/MaestroVelPred.cpp b/Source/MaestroVelPred.cpp
index deece519..9353615b 100644
--- a/Source/MaestroVelPred.cpp
+++ b/Source/MaestroVelPred.cpp
@@ -301,7 +301,7 @@ void Maestro::VelPredInterface(
     Array4<const Real> const Ipv, Array4<Real> const ulx,
     Array4<Real> const urx, Array4<Real> const uimhx, Array4<Real> const uly,
     Array4<Real> const ury, Array4<Real> const uimhy, const Box& domainBox,
-    const amrex::GpuArray<Real, AMREX_SPACEDIM> dx) {
+    const amrex::GpuArray<Real, AMREX_SPACEDIM>& dx) const {
     // timer for profiling
     BL_PROFILE_VAR("Maestro::VelPredInterface()", VelPredInterface);
 
@@ -326,8 +326,8 @@ void Maestro::VelPredInterface(
     const auto domlo = domainBox.loVect3d();
     const auto domhi = domainBox.hiVect3d();
 
-    int bclo = phys_bc[0];
-    int bchi = phys_bc[AMREX_SPACEDIM];
+    int bclo = phys_bc.lo(0);
+    int bchi = phys_bc.hi(0);
 
     ParallelFor(mxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         if (ppm_type == 0) {
@@ -434,8 +434,8 @@ void Maestro::VelPredInterface(
     });
 
     // y-direction
-    bclo = phys_bc[1];
-    bchi = phys_bc[AMREX_SPACEDIM + 1];
+    bclo = phys_bc.lo(1);
+    bchi = phys_bc.hi(1);
 
     ParallelFor(mybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         if (ppm_type == 0) {
@@ -551,7 +551,7 @@ void Maestro::VelPredVelocities(
     Array4<const Real> const uimhx, Array4<const Real> const uly,
     Array4<const Real> const ury, Array4<const Real> const uimhy,
     Array4<const Real> const force, Array4<const Real> const w0_cart_in,
-    const Box& domainBox, const amrex::GpuArray<Real, AMREX_SPACEDIM> dx) {
+    const Box& domainBox, const amrex::GpuArray<Real, AMREX_SPACEDIM>& dx) const {
     // timer for profiling
     BL_PROFILE_VAR("Maestro::VelPredVelocities()", VelPredVelocities);
 
@@ -574,8 +574,8 @@ void Maestro::VelPredVelocities(
     const auto domlo = domainBox.loVect3d();
     const auto domhi = domainBox.hiVect3d();
 
-    int bclo = phys_bc[0];
-    int bchi = phys_bc[AMREX_SPACEDIM];
+    int bclo = phys_bc.lo(0);
+    int bchi = phys_bc.hi(0);
 
     ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         // use the traced force if ppm_trace_forces = 1
@@ -649,8 +649,8 @@ void Maestro::VelPredVelocities(
     });
 
     // y-direction
-    bclo = phys_bc[1];
-    bchi = phys_bc[AMREX_SPACEDIM + 1];
+    bclo = phys_bc.lo(1);
+    bchi = phys_bc.hi(1);
 
     ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         // use the traced force if ppm_trace_forces = 1
@@ -742,7 +742,7 @@ void Maestro::VelPredInterface(
     Array4<Real> const ulx, Array4<Real> const urx, Array4<Real> const uimhx,
     Array4<Real> const uly, Array4<Real> const ury, Array4<Real> const uimhy,
     Array4<Real> const ulz, Array4<Real> const urz, Array4<Real> const uimhz,
-    const Box& domainBox, const amrex::GpuArray<Real, AMREX_SPACEDIM> dx) {
+    const Box& domainBox, const amrex::GpuArray<Real, AMREX_SPACEDIM>& dx) const {
     // timer for profiling
     BL_PROFILE_VAR("Maestro::VelPredInterface()", VelPredInterface);
 
@@ -773,8 +773,8 @@ void Maestro::VelPredInterface(
     const auto domlo = domainBox.loVect3d();
     const auto domhi = domainBox.hiVect3d();
 
-    int bclo = phys_bc[0];
-    int bchi = phys_bc[AMREX_SPACEDIM];
+    int bclo = phys_bc.lo(0);
+    int bchi = phys_bc.hi(0);
 
     ParallelFor(mxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         if (ppm_type == 0) {
@@ -890,8 +890,8 @@ void Maestro::VelPredInterface(
     });
 
     // y-direction
-    bclo = phys_bc[1];
-    bchi = phys_bc[AMREX_SPACEDIM + 1];
+    bclo = phys_bc.lo(1);
+    bchi = phys_bc.hi(1);
 
     ParallelFor(mybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         if (ppm_type == 0) {
@@ -1009,8 +1009,8 @@ void Maestro::VelPredInterface(
     });
 
     // z-direction
-    bclo = phys_bc[2];
-    bchi = phys_bc[AMREX_SPACEDIM + 2];
+    bclo = phys_bc.lo(2);
+    bchi = phys_bc.hi(2);
 
     ParallelFor(mzbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         if (ppm_type == 0) {
@@ -1138,7 +1138,7 @@ void Maestro::VelPredTransverse(
     Array4<Real> const uimhyz, Array4<Real> const uimhzy,
     Array4<Real> const vimhxz, Array4<Real> const vimhzx,
     Array4<Real> const wimhxy, Array4<Real> const wimhyx, const Box& domainBox,
-    const amrex::GpuArray<Real, AMREX_SPACEDIM> dx) {
+    const amrex::GpuArray<Real, AMREX_SPACEDIM>& dx) const {
     // timer for profiling
     BL_PROFILE_VAR("Maestro::VelPredTransverse()", VelPredTransverse);
 
@@ -1155,15 +1155,13 @@ void Maestro::VelPredTransverse(
     const auto domhi = domainBox.hiVect3d();
     const auto rel_eps_local = rel_eps;
 
-    GpuArray<int, AMREX_SPACEDIM * 2> physbc;
-    for (int n = 0; n < AMREX_SPACEDIM * 2; ++n) {
-        physbc[n] = phys_bc[n];
-    }
-
     // uimhyz, 1, 2
     Box imhbox = amrex::grow(mfi.tilebox(), 0, 1);
     imhbox = amrex::growHi(imhbox, 1, 1);
 
+    int bclo = phys_bc.lo(1);
+    int bchi = phys_bc.hi(1);
+
     ParallelFor(imhbox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         // extrapolate to faces
         Real ulyz = uly(i, j, k, 0) -
@@ -1176,7 +1174,7 @@ void Maestro::VelPredTransverse(
 
         // impose lo side bc's
         if (j == domlo[1]) {
-            switch (physbc[1]) {
+            switch (bclo) {
             case amrex::PhysBCType::inflow:
                 ulyz = utilde(i, j - 1, k, 0);
                 uryz = utilde(i, j - 1, k, 0);
@@ -1200,7 +1198,7 @@ void Maestro::VelPredTransverse(
 
             // impose hi side bc's
         } else if (j == domhi[1] + 1) {
-            switch (physbc[AMREX_SPACEDIM + 1]) {
+            switch (bchi) {
             case amrex::PhysBCType::inflow:
                 ulyz = utilde(i, j, k, 0);
                 uryz = utilde(i, j, k, 0);
@@ -1234,6 +1232,9 @@ void Maestro::VelPredTransverse(
     imhbox = amrex::grow(mfi.tilebox(), 0, 1);
     imhbox = amrex::growHi(imhbox, 2, 1);
 
+    bclo = phys_bc.lo(2);
+    bchi = phys_bc.hi(2);
+
     ParallelFor(imhbox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         // extrapolate to faces
         Real ulzy = ulz(i, j, k, 0) -
@@ -1246,7 +1247,7 @@ void Maestro::VelPredTransverse(
 
         // impose lo side bc's
         if (k == domlo[2]) {
-            switch (physbc[2]) {
+            switch (bclo) {
             case amrex::PhysBCType::inflow:
                 ulzy = utilde(i, j, k - 1, 0);
                 urzy = utilde(i, j, k - 1, 0);
@@ -1270,7 +1271,7 @@ void Maestro::VelPredTransverse(
 
             // impose hi side bc's
         } else if (k == domhi[2] + 1) {
-            switch (physbc[AMREX_SPACEDIM + 2]) {
+            switch (bchi) {
             case amrex::PhysBCType::inflow:
                 ulzy = utilde(i, j, k, 0);
                 urzy = utilde(i, j, k, 0);
@@ -1304,6 +1305,9 @@ void Maestro::VelPredTransverse(
     imhbox = amrex::grow(mfi.tilebox(), 1, 1);
     imhbox = amrex::growHi(imhbox, 0, 1);
 
+    bclo = phys_bc.lo(0);
+    bchi = phys_bc.hi(0);
+
     ParallelFor(imhbox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         // extrapolate to faces
         Real vlxz = ulx(i, j, k, 1) -
@@ -1316,7 +1320,7 @@ void Maestro::VelPredTransverse(
 
         // impose lo side bc's
         if (i == domlo[0]) {
-            switch (physbc[0]) {
+            switch (bclo) {
             case amrex::PhysBCType::inflow:
                 vlxz = utilde(i - 1, j, k, 1);
                 vrxz = utilde(i - 1, j, k, 1);
@@ -1340,7 +1344,7 @@ void Maestro::VelPredTransverse(
 
             // impose hi side bc's
         } else if (i == domhi[0] + 1) {
-            switch (physbc[AMREX_SPACEDIM]) {
+            switch (bchi) {
             case amrex::PhysBCType::inflow:
                 vlxz = utilde(i, j, k, 1);
                 vrxz = utilde(i, j, k, 1);
@@ -1374,6 +1378,9 @@ void Maestro::VelPredTransverse(
     imhbox = amrex::grow(mfi.tilebox(), 1, 1);
     imhbox = amrex::growHi(imhbox, 2, 1);
 
+    bclo = phys_bc.lo(2);
+    bchi = phys_bc.hi(2);
+
     ParallelFor(imhbox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         // extrapolate to faces
         Real vlzx = ulz(i, j, k, 1) -
@@ -1386,7 +1393,7 @@ void Maestro::VelPredTransverse(
 
         // impose lo side bc's
         if (k == domlo[2]) {
-            switch (physbc[2]) {
+            switch (bclo) {
             case amrex::PhysBCType::inflow:
                 vlzx = utilde(i, j, k - 1, 1);
                 vrzx = utilde(i, j, k - 1, 1);
@@ -1410,7 +1417,7 @@ void Maestro::VelPredTransverse(
 
             // impose hi side bc's
         } else if (k == domhi[2] + 1) {
-            switch (physbc[AMREX_SPACEDIM + 2]) {
+            switch (bchi) {
             case amrex::PhysBCType::inflow:
                 vlzx = utilde(i, j, k, 1);
                 vrzx = utilde(i, j, k, 1);
@@ -1444,6 +1451,9 @@ void Maestro::VelPredTransverse(
     imhbox = amrex::grow(mfi.tilebox(), 2, 1);
     imhbox = amrex::growHi(imhbox, 0, 1);
 
+    bclo = phys_bc.lo(0);
+    bchi = phys_bc.hi(0);
+
     ParallelFor(imhbox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         // extrapolate to faces
         Real wlxy = ulx(i, j, k, 2) -
@@ -1456,7 +1466,7 @@ void Maestro::VelPredTransverse(
 
         // impose lo side bc's
         if (i == domlo[0]) {
-            switch (physbc[0]) {
+            switch (bclo) {
             case amrex::PhysBCType::inflow:
                 wlxy = utilde(i - 1, j, k, 2);
                 wrxy = utilde(i - 1, j, k, 2);
@@ -1480,7 +1490,7 @@ void Maestro::VelPredTransverse(
 
             // impose hi side bc's
         } else if (i == domhi[0] + 1) {
-            switch (physbc[AMREX_SPACEDIM]) {
+            switch (bchi) {
             case amrex::PhysBCType::inflow:
                 wlxy = utilde(i, j, k, 2);
                 wrxy = utilde(i, j, k, 2);
@@ -1514,6 +1524,9 @@ void Maestro::VelPredTransverse(
     imhbox = amrex::grow(mfi.tilebox(), 2, 1);
     imhbox = amrex::growHi(imhbox, 1, 1);
 
+    bclo = phys_bc.lo(1);
+    bchi = phys_bc.hi(1);
+
     ParallelFor(imhbox, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         // extrapolate to faces
         Real wlyx = uly(i, j, k, 2) -
@@ -1526,7 +1539,7 @@ void Maestro::VelPredTransverse(
 
         // impose lo side bc's
         if (j == domlo[1]) {
-            switch (physbc[1]) {
+            switch (bclo) {
             case amrex::PhysBCType::inflow:
                     wlyx = utilde(i, j - 1, k, 2);
                     wryx = utilde(i, j - 1, k, 2);
@@ -1550,7 +1563,7 @@ void Maestro::VelPredTransverse(
 
             // impose hi side bc's
         } else if (j == domhi[1] + 1) {
-            switch (physbc[AMREX_SPACEDIM + 1]) {
+            switch (bchi) {
             case amrex::PhysBCType::inflow:
                 wlyx = utilde(i, j, k, 2);
                 wryx = utilde(i, j, k, 2);
@@ -1598,7 +1611,7 @@ void Maestro::VelPredVelocities(
     Array4<const Real> const vimhzx, Array4<const Real> const wimhxy,
     Array4<const Real> const wimhyx, Array4<const Real> const force,
     Array4<const Real> const w0_cart_in, const Box& domainBox,
-    const amrex::GpuArray<Real, AMREX_SPACEDIM> dx) {
+    const amrex::GpuArray<Real, AMREX_SPACEDIM>& dx) const {
     // timer for profiling
     BL_PROFILE_VAR("Maestro::VelPredVelocities()", VelPredVelocities);
 
@@ -1622,10 +1635,9 @@ void Maestro::VelPredVelocities(
     const auto domlo = domainBox.loVect3d();
     const auto domhi = domainBox.hiVect3d();
 
-    GpuArray<int, AMREX_SPACEDIM * 2> physbc;
-    for (int n = 0; n < AMREX_SPACEDIM * 2; ++n) {
-        physbc[n] = phys_bc[n];
-    }
+
+    int bclo = phys_bc.lo(0);
+    int bchi = phys_bc.hi(0);
 
     // x-direction
     ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
@@ -1669,7 +1681,7 @@ void Maestro::VelPredVelocities(
 
         // impose lo side bc's
         if (i == domlo[0]) {
-            switch (physbc[0]) {
+            switch (bclo) {
             case amrex::PhysBCType::inflow:
                 umac(i, j, k) = utilde(i - 1, j, k, 0);
                 break;
@@ -1691,7 +1703,7 @@ void Maestro::VelPredVelocities(
 
             // impose hi side bc's
         } else if (i == domhi[0] + 1) {
-            switch (physbc[AMREX_SPACEDIM]) {
+            switch (bchi) {
             case amrex::PhysBCType::inflow:
                 umac(i, j, k) = utilde(i, j, k, 0);
                 break;
@@ -1713,6 +1725,9 @@ void Maestro::VelPredVelocities(
         }
     });
 
+    bclo = phys_bc.lo(1);
+    bchi = phys_bc.hi(1);
+
     // y-direction
     ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         // use the traced force if ppm_trace_forces = 1
@@ -1755,7 +1770,7 @@ void Maestro::VelPredVelocities(
 
         // impose lo side bc's
         if (j == domlo[1]) {
-            switch (physbc[1]) {
+            switch (bclo) {
             case amrex::PhysBCType::inflow:
                 vmac(i, j, k) = utilde(i, j - 1, k, 1);
                 break;
@@ -1777,7 +1792,7 @@ void Maestro::VelPredVelocities(
 
             // impose hi side bc's
         } else if (j == domhi[1] + 1) {
-            switch (physbc[AMREX_SPACEDIM + 1]) {
+            switch (bchi) {
             case amrex::PhysBCType::inflow:
                 vmac(i, j, k) = utilde(i, j, k, 1);
                 break;
@@ -1799,6 +1814,9 @@ void Maestro::VelPredVelocities(
         }
     });
 
+    bclo = phys_bc.lo(2);
+    bchi = phys_bc.hi(2);
+
     // z-direction
     ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
         // use the traced force if ppm_trace_forces = 1
@@ -1851,7 +1869,7 @@ void Maestro::VelPredVelocities(
 
         // impose hi side bc's
         if (k == domlo[2]) {
-            switch (physbc[2]) {
+            switch (bclo) {
             case amrex::PhysBCType::inflow:
                 wmac(i, j, k) = utilde(i, j, k - 1, 2);
                 break;
@@ -1873,7 +1891,7 @@ void Maestro::VelPredVelocities(
 
             // impose lo side bc's
         } else if (k == domhi[2] + 1) {
-            switch (physbc[AMREX_SPACEDIM + 2]) {
+            switch (bchi) {
             case amrex::PhysBCType::inflow:
                 wmac(i, j, k) = utilde(i, j, k, 2);
                 break;