diff --git a/env/bash b/env/bash
index 6731622..6c0c4cc 100644
--- a/env/bash
+++ b/env/bash
@@ -104,6 +104,7 @@ elif [[ $PARTITION == "darwin-volta-x86" ]]; then
     module load clang/12.0.1
     module load nvhpc/23.5
     module load cuda/11.4.2
+    module load gcc/9.4.0
     module load hdf5/1.12.2
     module load miniconda3
     source /projects/opt/centos8/x86_64/miniconda3/py39_4.12.0/bin/activate /usr/projects/jovian/dependencies/python/volta-x86-py39
diff --git a/external/parthenon b/external/parthenon
index 89099e3..ed14989 160000
--- a/external/parthenon
+++ b/external/parthenon
@@ -1 +1 @@
-Subproject commit 89099e3aa7b642fbf60fc43cf93c9ce3d78bcec9
+Subproject commit ed149893027ec885a1287e928ac901431027491a
diff --git a/src/jaybenne/jaybenne.cpp b/src/jaybenne/jaybenne.cpp
index 0a167cf..e17b1fb 100644
--- a/src/jaybenne/jaybenne.cpp
+++ b/src/jaybenne/jaybenne.cpp
@@ -343,8 +343,12 @@ TaskStatus UpdateDerivedTransportFields(MeshData<Real> *md, const Real dt) {
 
           // get current, lower, upper neighbor block levels in x-direction
           const Real rlev = static_cast<Real>(vmesh.GetLevel(b, 0, 0, 0));
-          const Real rlev_lx = static_cast<Real>(vmesh.GetLevel(b, 0, 0, -1));
-          const Real rlev_ux = static_cast<Real>(vmesh.GetLevel(b, 0, 0, 1));
+          const Real rlev_lx = (vmesh.IsPhysicalBoundary(b, 0, 0, -1))
+                                   ? rlev
+                                   : static_cast<Real>(vmesh.GetLevel(b, 0, 0, -1));
+          const Real rlev_ux = (vmesh.IsPhysicalBoundary(b, 0, 0, 1))
+                                   ? rlev
+                                   : static_cast<Real>(vmesh.GetLevel(b, 0, 0, 1));
 
           // calculate neighbor dx values
           const Real dx_lx = i == ib.s ? std::pow(2.0, rlev - rlev_lx) * dx_i : dx_i;
@@ -391,8 +395,12 @@ TaskStatus UpdateDerivedTransportFields(MeshData<Real> *md, const Real dt) {
 
             // get current, lower, upper neighbor block levels in x-direction
             const Real rlev = static_cast<Real>(vmesh.GetLevel(b, 0, 0, 0));
-            const Real rlev_ly = static_cast<Real>(vmesh.GetLevel(b, 0, -1, 0));
-            const Real rlev_uy = static_cast<Real>(vmesh.GetLevel(b, 0, 1, 0));
+            const Real rlev_ly = (vmesh.IsPhysicalBoundary(b, 0, -1, 0))
+                                     ? rlev
+                                     : static_cast<Real>(vmesh.GetLevel(b, 0, -1, 0));
+            const Real rlev_uy = (vmesh.IsPhysicalBoundary(b, 0, 1, 0))
+                                     ? rlev
+                                     : static_cast<Real>(vmesh.GetLevel(b, 0, 1, 0));
 
             // calculate neighbor dx values
             const Real dx_ly = j == jb.s ? std::pow(2.0, rlev - rlev_ly) * dx_j : dx_j;
@@ -440,8 +448,12 @@ TaskStatus UpdateDerivedTransportFields(MeshData<Real> *md, const Real dt) {
 
             // get current, lower, upper neighbor block levels in x-direction
             const Real rlev = static_cast<Real>(vmesh.GetLevel(b, 0, 0, 0));
-            const Real rlev_lz = static_cast<Real>(vmesh.GetLevel(b, -1, 0, 0));
-            const Real rlev_uz = static_cast<Real>(vmesh.GetLevel(b, 1, 0, 0));
+            const Real rlev_lz = (vmesh.IsPhysicalBoundary(b, -1, 0, 0))
+                                     ? rlev
+                                     : static_cast<Real>(vmesh.GetLevel(b, -1, 0, 0));
+            const Real rlev_uz = (vmesh.IsPhysicalBoundary(b, 1, 0, 0))
+                                     ? rlev
+                                     : static_cast<Real>(vmesh.GetLevel(b, 1, 0, 0));
 
             // calculate neighbor dx values
             const Real dx_lz = k == kb.s ? std::pow(2.0, rlev - rlev_lz) * dx_k : dx_k;