From 0bd2d1aeffdd27ad534b70cb042f942c7f97a41b Mon Sep 17 00:00:00 2001
From: Ann Almgren <asalmgren@lbl.gov>
Date: Mon, 29 Jul 2024 19:26:50 -0700
Subject: [PATCH 1/9] add option for direction-dependent bc's

---
 EBGodunov/hydro_ebgodunov_edge_state_2D.cpp   | 16 +++--
 EBGodunov/hydro_ebgodunov_edge_state_3D.cpp   | 24 +++----
 ...hydro_ebgodunov_extrap_vel_to_faces_2D.cpp | 18 ++++--
 ...hydro_ebgodunov_extrap_vel_to_faces_3D.cpp | 27 +++++---
 EBMOL/hydro_ebmol_edge_state_K.H              | 18 ++++--
 EBMOL/hydro_ebmol_extrap_vel_to_faces_box.cpp | 18 ++++--
 Godunov/hydro_godunov_edge_state_2D.cpp       | 16 +++--
 Godunov/hydro_godunov_edge_state_3D.cpp       | 27 ++++----
 .../hydro_godunov_extrap_vel_to_faces_2D.cpp  | 16 ++---
 .../hydro_godunov_extrap_vel_to_faces_3D.cpp  | 34 +++++-----
 Godunov/hydro_godunov_plm.H                   | 18 ++++--
 Godunov/hydro_godunov_ppm.H                   | 44 ++++++++-----
 MOL/hydro_mol_edge_state_K.H                  | 18 ++++--
 MOL/hydro_mol_extrap_vel_to_faces_box.cpp     | 18 ++++--
 Utils/hydro_bcs_K.H                           | 63 ++++++++++++-------
 15 files changed, 232 insertions(+), 143 deletions(-)

diff --git a/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp b/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp
index 97d890970..f19b07837 100644
--- a/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp
+++ b/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp
@@ -109,7 +109,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+            Real uad = u_mac(i,j,k);
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
             xlo(i,j,k,n) = lo;
             xhi(i,j,k,n) = hi;
@@ -121,7 +122,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+            Real vad = v_mac(i,j,k);
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
             ylo(i,j,k,n) = lo;
             yhi(i,j,k,n) = hi;
@@ -143,7 +145,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             l_yzlo = ylo(i,j,k,n);
             l_yzhi = yhi(i,j,k,n);
             Real vad = v_mac(i,j,k);
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
             Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
             Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
@@ -220,7 +222,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             }
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+            Real uad = u_mac(i,j,k);
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
             if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
             {
@@ -259,7 +262,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             l_xzhi = xhi(i,j,k,n);
 
             Real uad = u_mac(i,j,k);
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
             Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
             Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
@@ -336,7 +339,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             }
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+            Real vad = v_mac(i,j,k);
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
             if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
             {
diff --git a/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp b/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
index 57c7d13c7..501ffb0d6 100644
--- a/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
+++ b/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
@@ -138,7 +138,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
             xlo(i,j,k,n) = lo;
             xhi(i,j,k,n) = hi;
@@ -156,7 +156,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
             ylo(i,j,k,n) = lo;
             yhi(i,j,k,n) = hi;
@@ -174,7 +174,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+            HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
             zlo(i,j,k,n) = lo;
             zhi(i,j,k,n) = hi;
@@ -208,7 +208,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
                                    q, divu, apx, apy, apz, vfrac_arr, v_mac, yed);
 
         Real wad = w_mac(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
         Real st = (wad >= 0.) ? l_zylo : l_zyhi;
         Real fu = (amrex::Math::abs(wad) < small_vel) ? Real(0.0) : Real(1.0);
@@ -225,7 +225,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
                                    q, divu, apx, apy, apz, vfrac_arr, w_mac, zed);
 
         Real vad = v_mac(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? Real(0.0) : Real(1.0);
@@ -304,7 +304,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             }
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
             if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
             {
@@ -345,7 +345,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
                                    q, divu, apx, apy, apz, vfrac_arr, w_mac, zed);
 
         Real uad = u_mac(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
@@ -362,7 +362,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
                                    q, divu, apx, apy, apz, vfrac_arr, u_mac, xed);
 
         Real wad = w_mac(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
         Real st = (wad >= 0.) ? l_zxlo : l_zxhi;
         Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0;
@@ -441,7 +441,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             }
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
             if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
             {
@@ -480,7 +480,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
                                    q, divu, apx, apy, apz, vfrac_arr, v_mac, yed);
 
         Real uad = u_mac(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
         Real st = (uad >= 0.) ? l_xylo : l_xyhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
@@ -497,7 +497,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
                                    q, divu, apx, apy, apz, vfrac_arr, u_mac, xed);
 
         Real vad = v_mac(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         Real st = (vad >= 0.) ? l_yxlo : l_yxhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
@@ -575,7 +575,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             }
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+            HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
             if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) )
             {
diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp
index 2d97a4a54..e65524c35 100644
--- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp
+++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp
@@ -63,7 +63,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, lo, hi,
+                                 bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
             xlo(i,j,k,n) = lo;
             xhi(i,j,k,n) = hi;
@@ -75,7 +76,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, lo, hi,
+                                 bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
             ylo(i,j,k,n) = lo;
             yhi(i,j,k,n) = hi;
@@ -103,7 +105,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
             l_yzlo = ylo(i,j,k,n);
             l_yzhi = yhi(i,j,k,n);
             Real vad = v_ad(i,j,k);
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad,
+                                 bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
             Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
             Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
@@ -188,7 +191,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
             }
         }
 
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, stl, sth,
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         // Prevent backflow
         if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
@@ -229,7 +233,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
             l_xzhi = xhi(i,j,k,n);
 
             Real uad = u_ad(i,j,k);
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad,
+                                 bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
             Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
             Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
@@ -310,7 +315,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
             }
         }
 
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, stl, sth,
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         // Prevent backflow
         if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
index ae03fd511..01664568e 100644
--- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
+++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
@@ -164,7 +164,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                                    q, divu, apx, apy, apz, vfrac_arr, v_ad, yedge);
 
         Real wad = w_ad(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, wad, wad,
+                             bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
 
 
         Real st = (wad >= 0.) ? l_zylo : l_zyhi;
@@ -182,7 +183,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                                    q, divu, apx, apy, apz, vfrac_arr, w_ad, zedge);
 
         Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad,
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? Real(0.0) : Real(1.0);
@@ -265,7 +267,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
         }
         }
 
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, stl, sth,
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         // Prevent backflow
         if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
@@ -315,7 +318,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                                    q, divu, apx, apy, apz, vfrac_arr, w_ad, zedge);
 
         Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad,
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
 
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
@@ -333,7 +337,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                                    q, divu, apx, apy, apz, vfrac_arr, u_ad, xedge);
 
         Real wad = w_ad(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, wad, wad,
+                            bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
 
 
         Real st = (wad >= 0.) ? l_zxlo : l_zxhi;
@@ -417,7 +422,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
         }
         }
 
-        HydroBC::SetYEdgeBCs( i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, stl, sth,
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         // Prevent backflow
         if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
@@ -467,7 +473,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                                    q, divu, apx, apy, apz, vfrac_arr, v_ad, yedge);
 
         Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, uad, uad,
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
 
         Real st = (uad >= 0.) ? l_xylo : l_xyhi;
@@ -489,7 +496,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                                    q, divu, apx, apy, apz, vfrac_arr, u_ad, xedge);
 
         Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, vad, vad,
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
 
         Real st = (vad >= 0.) ? l_yxlo : l_yxhi;
@@ -572,7 +580,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
         }
         }
 
-        HydroBC::SetZEdgeBCs( i, j, k, n, q, stl, sth, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, stl, sth,
+                             bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
 
 
         // Prevent backflow
diff --git a/EBMOL/hydro_ebmol_edge_state_K.H b/EBMOL/hydro_ebmol_edge_state_K.H
index dd7f438fd..fa42b27a3 100644
--- a/EBMOL/hydro_ebmol_edge_state_K.H
+++ b/EBMOL/hydro_ebmol_edge_state_K.H
@@ -141,7 +141,8 @@ amrex::Real hydro_ebmol_xedge_state_extdir ( AMREX_D_DECL(int i, int j, int k),
 
         qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
 
-        HydroBC::SetXEdgeBCs(i, j, k, 0, q, qmns, qpls, d_bcrec[n].lo(0), domain_ilo, d_bcrec[n].hi(0), domain_ihi, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, 0, q, qmns, qpls, qmns, qpls,
+                             d_bcrec[n].lo(0), domain_ilo, d_bcrec[n].hi(0), domain_ihi, is_velocity);
 
         if ( (i==domain_ilo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) )
         {
@@ -249,7 +250,8 @@ amrex::Real hydro_ebmol_xedge_state ( AMREX_D_DECL(int i, int j, int k), int n,
 
     qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
 
-    HydroBC::SetXEdgeBCs(i, j, k, 0, q, qmns, qpls, d_bcrec[n].lo(0), domain_ilo, d_bcrec[n].hi(0), domain_ihi, is_velocity);
+    HydroBC::SetXEdgeBCs(i, j, k, 0, q, qmns, qpls, qmns, qpls,
+                         d_bcrec[n].lo(0), domain_ilo, d_bcrec[n].hi(0), domain_ihi, is_velocity);
 
     if ( (i==domain_ilo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) )
     {
@@ -403,7 +405,8 @@ amrex::Real hydro_ebmol_yedge_state_extdir ( AMREX_D_DECL(int i, int j, int k),
                                         + delta_y * slopes_eb_lo[1];
 #endif
 
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, qmns, qpls, d_bcrec[n].lo(1), domain_jlo, d_bcrec[n].hi(1), domain_jhi, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, qmns, qpls, qmns, qpls,
+                             d_bcrec[n].lo(1), domain_jlo, d_bcrec[n].hi(1), domain_jhi, is_velocity);
 
         if ( (j==domain_jlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) )
         {
@@ -515,7 +518,8 @@ amrex::Real hydro_ebmol_yedge_state ( AMREX_D_DECL(int i, int j, int k), int n,
 
     qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
 
-    HydroBC::SetYEdgeBCs(i, j, k, n, q, qmns, qpls, d_bcrec[n].lo(1), domain_jlo, d_bcrec[n].hi(1), domain_jhi, is_velocity);
+    HydroBC::SetYEdgeBCs(i, j, k, n, q, qmns, qpls, qmns, qpls,
+                         d_bcrec[n].lo(1), domain_jlo, d_bcrec[n].hi(1), domain_jhi, is_velocity);
 
     if ( (j==domain_jlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) )
     {
@@ -658,7 +662,8 @@ amrex::Real hydro_ebmol_zedge_state_extdir ( int i, int j, int k, int n,
 
         qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
 
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, qmns, qpls, d_bcrec[n].lo(2), domain_klo, d_bcrec[n].hi(2), domain_khi, is_velocity);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, qmns, qpls, qmns, qpls,
+                             d_bcrec[n].lo(2), domain_klo, d_bcrec[n].hi(2), domain_khi, is_velocity);
 
         if ( (k==domain_klo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) )
         {
@@ -752,7 +757,8 @@ amrex::Real hydro_ebmol_zedge_state ( int i, int j, int k, int n,
 
     qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
 
-    HydroBC::SetZEdgeBCs(i, j, k, n, q, qmns, qpls, d_bcrec[n].lo(2), domain_klo, d_bcrec[n].hi(2), domain_khi, is_velocity);
+    HydroBC::SetZEdgeBCs(i, j, k, n, q, qmns, qpls, qmns, qpls,
+                         d_bcrec[n].lo(2), domain_klo, d_bcrec[n].hi(2), domain_khi, is_velocity);
 
     if ( (k==domain_klo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) )
     {
diff --git a/EBMOL/hydro_ebmol_extrap_vel_to_faces_box.cpp b/EBMOL/hydro_ebmol_extrap_vel_to_faces_box.cpp
index 284702435..d1b475d38 100644
--- a/EBMOL/hydro_ebmol_extrap_vel_to_faces_box.cpp
+++ b/EBMOL/hydro_ebmol_extrap_vel_to_faces_box.cpp
@@ -178,7 +178,8 @@ EBMOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx,
 #endif
                umns = amrex::max(amrex::min(umns, cc_umax), cc_umin);
 
-               HydroBC::SetXEdgeBCs(i, j, k, 0, vcc, umns, upls, d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true);
+               HydroBC::SetXEdgeBCs(i, j, k, 0, vcc, umns, upls, umns, upls,
+                                    d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true);
 
                if ( (i==domain_ilo) && (d_bcrec[0].lo(0) == BCType::foextrap || d_bcrec[0].lo(0) == BCType::hoextrap) )
                {
@@ -267,7 +268,8 @@ EBMOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx,
 #endif
                umns = amrex::max(amrex::min(umns, cc_umax), cc_umin);
 
-               HydroBC::SetXEdgeBCs(i, j, k, 0, vcc, umns, upls, d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true);
+               HydroBC::SetXEdgeBCs(i, j, k, 0, vcc, umns, upls, umns, upls,
+                                    d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true);
 
                if ( (i==domain_ilo) && (d_bcrec[0].lo(0) == BCType::foextrap || d_bcrec[0].lo(0) == BCType::hoextrap) )
                {
@@ -413,7 +415,8 @@ EBMOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx,
 
                vmns = amrex::max(amrex::min(vmns, cc_vmax), cc_vmin);
 
-               HydroBC::SetYEdgeBCs(i, j, k, 1, vcc, vmns, vpls, d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true);
+               HydroBC::SetYEdgeBCs(i, j, k, 1, vcc, vmns, vpls, vmns, vpls,
+                                    d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true);
 
                if ( (j==domain_jlo) && (d_bcrec[1].lo(1) == BCType::foextrap || d_bcrec[1].lo(1) == BCType::hoextrap) )
                {
@@ -505,7 +508,8 @@ EBMOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx,
 
                vmns = amrex::max(amrex::min(vmns, cc_vmax), cc_vmin);
 
-               HydroBC::SetYEdgeBCs(i, j, k, 1, vcc, vmns, vpls, d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true);
+               HydroBC::SetYEdgeBCs(i, j, k, 1, vcc, vmns, vpls, vmns, vpls,
+                                    d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true);
 
                if ( (j==domain_jlo) && (d_bcrec[1].lo(1) == BCType::foextrap || d_bcrec[1].lo(1) == BCType::hoextrap) )
                {
@@ -633,7 +637,8 @@ EBMOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx,
 
                wmns = amrex::max(amrex::min(wmns, cc_wmax), cc_wmin);
 
-               HydroBC::SetZEdgeBCs(i, j, k, 2, vcc, wmns, wpls, d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true);
+               HydroBC::SetZEdgeBCs(i, j, k, 2, vcc, wmns, wpls, wmns, wpls,
+                                    d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true);
 
                if ( (k==domain_klo) && (d_bcrec[2].lo(2) == BCType::foextrap || d_bcrec[2].lo(2) == BCType::hoextrap) )
                {
@@ -713,7 +718,8 @@ EBMOL::ExtrapVelToFacesBox ( AMREX_D_DECL( Box const& ubx,
 
                wmns = amrex::max(amrex::min(wmns, cc_wmax), cc_wmin);
 
-               HydroBC::SetZEdgeBCs(i, j, k, 2, vcc, wmns, wpls, d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true);
+               HydroBC::SetZEdgeBCs(i, j, k, 2, vcc, wmns, wpls, wmns, wpls,
+                                    d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true);
 
                if ( (k==domain_klo) && (d_bcrec[2].lo(2) == BCType::foextrap || d_bcrec[2].lo(2) == BCType::hoextrap) )
                {
diff --git a/Godunov/hydro_godunov_edge_state_2D.cpp b/Godunov/hydro_godunov_edge_state_2D.cpp
index 2c77a4e44..cf735be23 100644
--- a/Godunov/hydro_godunov_edge_state_2D.cpp
+++ b/Godunov/hydro_godunov_edge_state_2D.cpp
@@ -139,7 +139,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        Real uad = umac(i,j,k);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
         xlo(i,j,k,n) = lo;
         xhi(i,j,k,n) = hi;
     },
@@ -156,7 +157,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        Real vad = vmac(i,j,k);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         ylo(i,j,k,n) = lo;
         yhi(i,j,k,n) = hi;
@@ -178,7 +180,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
         l_yzlo = ylo(i,j,k,n);
         l_yzhi = yhi(i,j,k,n);
         Real vad = vmac(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
@@ -227,7 +229,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        Real uad = umac(i,j,k);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
         if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
         {
@@ -261,7 +264,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
         l_xzhi = xhi(i,j,k,n);
 
         Real uad = umac(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
@@ -311,7 +314,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        Real vad = vmac(i,j,k);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
         {
diff --git a/Godunov/hydro_godunov_edge_state_3D.cpp b/Godunov/hydro_godunov_edge_state_3D.cpp
index af7d3c094..fcccd6b91 100644
--- a/Godunov/hydro_godunov_edge_state_3D.cpp
+++ b/Godunov/hydro_godunov_edge_state_3D.cpp
@@ -161,7 +161,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
         xlo(i,j,k,n) = lo;
         xhi(i,j,k,n) = hi;
         Real st = (uval) ? lo : hi;
@@ -184,7 +184,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         ylo(i,j,k,n) = lo;
         yhi(i,j,k,n) = hi;
@@ -207,7 +207,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
         zlo(i,j,k,n) = lo;
         zhi(i,j,k,n) = hi;
@@ -236,7 +236,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
                               q, divu, vmac, Imy);
 
         Real wad = wmac(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
         Real st = (wad >= 0.) ? l_zylo : l_zyhi;
         Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0;
@@ -253,7 +253,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
                               q, divu, wmac, Imz);
 
         Real vad = vmac(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
@@ -301,7 +301,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        Real uad = umac(i,j,k);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
         if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
         {
@@ -337,7 +338,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
                               q, divu, wmac, Imz);
 
         Real uad = umac(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
@@ -354,7 +355,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
                               q, divu, umac, Imx);
 
         Real wad = wmac(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
         Real st = (wad >= 0.) ? l_zxlo : l_zxhi;
         Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0;
@@ -400,7 +401,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        Real vad = vmac(i,j,k);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
         {
@@ -436,7 +438,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
                               q, divu, vmac, Imy);
 
         Real uad = umac(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
         Real st = (uad >= 0.) ? l_xylo : l_xyhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
@@ -453,7 +455,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
                               q, divu, umac, Imx);
 
         Real vad = vmac(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         Real st = (vad >= 0.) ? l_yxlo : l_yxhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
@@ -500,7 +502,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+        Real wad = wmac(i,j,k);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
         if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) )
         {
diff --git a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp
index d4d6ef22f..083a71cf5 100644
--- a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp
+++ b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp
@@ -154,7 +154,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx,
         }
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetXEdgeBCs(i, j, k, n, vel, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, vel, lo, hi, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         Real st = ( (lo+hi) >= 0.) ? lo : hi;
         bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) );
@@ -175,7 +175,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx,
         }
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetYEdgeBCs(i, j, k, n, vel, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, vel, lo, hi, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         Real st = ( (lo+hi) >= 0.) ? lo : hi;
         bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) );
@@ -241,7 +241,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
         }
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         xlo(i,j,k,n) = lo;
         xhi(i,j,k,n) = hi;
@@ -258,7 +258,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
         }
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         ylo(i,j,k,n) = lo;
         yhi(i,j,k,n) = hi;
@@ -282,7 +282,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
         l_yzhi = yhi(i,j,k,n);
 
         Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
@@ -304,7 +304,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
             sth += 0.5 * l_dt * f(i  ,j,k,n);
         }
 
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
         {
@@ -338,7 +338,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
         l_xzhi = xhi(i,j,k,n);
 
         Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
 
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
@@ -363,7 +363,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
            sth += 0.5 * l_dt * f(i,j  ,k,n);
         }
 
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
         {
diff --git a/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp b/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp
index 920f19995..d6a2019e5 100644
--- a/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp
+++ b/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp
@@ -172,7 +172,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx,
         }
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetXEdgeBCs(i, j, k, n, vel, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, vel, lo, hi, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         Real st = ( (lo+hi) >= 0.) ? lo : hi;
         bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) );
@@ -193,7 +193,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx,
         }
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetYEdgeBCs(i, j, k, n, vel, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, vel, lo, hi, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         Real st = ( (lo+hi) >= 0.) ? lo : hi;
         bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) );
@@ -214,7 +214,7 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx,
         }
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetZEdgeBCs(i, j, k, n, vel, lo, hi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
+        HydroBC::SetZEdgeBCs(i, j, k, n, vel, lo, hi, lo, hi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
 
         Real st = ( (lo+hi) >= 0.) ? lo : hi;
         bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) );
@@ -289,7 +289,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
         Real uad = u_ad(i,j,k);
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         xlo(i,j,k,n) = lo;
         xhi(i,j,k,n) = hi;
@@ -312,7 +312,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
         Real vad = v_ad(i,j,k);
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         ylo(i,j,k,n) = lo;
         yhi(i,j,k,n) = hi;
@@ -336,7 +336,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
         Real wad = w_ad(i,j,k);
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
 
         zlo(i,j,k,n) = lo;
         zhi(i,j,k,n) = hi;
@@ -384,7 +384,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                               q, divu, v_ad, yedge);
 
         Real wad = w_ad(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, wad, wad,
+                             bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
 
 
         Real st = (wad >= 0.) ? l_zylo : l_zyhi;
@@ -402,7 +403,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                               q, divu, w_ad, zedge);
 
         Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad,
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
@@ -430,7 +432,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
             sth += Real(0.5) * l_dt * f(i  ,j,k,n);
         }
 
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
         {
@@ -471,7 +473,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                               q, divu, w_ad, zedge);
 
         Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad,
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
 
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
@@ -489,7 +492,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                               q, divu, u_ad, xedge);
 
         Real wad = w_ad(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, wad, wad,
+                             bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
 
         Real st = (wad >= 0.) ? l_zxlo : l_zxhi;
         Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0;
@@ -518,7 +522,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
            sth += Real(0.5) * l_dt * f(i,j  ,k,n);
         }
 
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
         {
@@ -560,7 +564,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                               q, divu, v_ad, yedge);
 
         Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
 
         Real st = (uad >= 0.) ? l_xylo : l_xyhi;
@@ -582,7 +586,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                               q, divu, u_ad, xedge);
 
         Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
 
         Real st = (vad >= 0.) ? l_yxlo : l_yxhi;
@@ -609,7 +613,7 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
            sth += Real(0.5) * l_dt * f(i,j,k  ,n);
         }
 
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, stl, sth, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
 
         if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) )
         {
diff --git a/Godunov/hydro_godunov_plm.H b/Godunov/hydro_godunov_plm.H
index 81f68c726..b7e7fec90 100644
--- a/Godunov/hydro_godunov_plm.H
+++ b/Godunov/hydro_godunov_plm.H
@@ -73,7 +73,8 @@ void PredictStateOnXFace ( const int i, const int j, const int k, const int n,
 
         int order = 4;
 
-        if (i == domain_ilo && (bc.lo(0) == BCType::ext_dir))
+        if (i == domain_ilo && ((bc.lo(0) == BCType::ext_dir) ||
+                                (bc.lo(0) == BCType::direction_dependent && umac >= 0.0)) )
         {
             umns = S(i-1,j,k,n);
 
@@ -88,7 +89,8 @@ void PredictStateOnXFace ( const int i, const int j, const int k, const int n,
             }
 
         }
-        else if (i == domain_ihi+1 && (bc.hi(0) == BCType::ext_dir))
+        else if (i == domain_ihi+1 && ((bc.hi(0) == BCType::ext_dir) ||
+                                       (bc.lo(0) == BCType::direction_dependent && umac <= 0.0)) )
         {
             upls = S(i  ,j,k,n);
 
@@ -140,7 +142,8 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n,
 
         int order = 4;
 
-        if (j == domain_jlo && (bc.lo(1) == BCType::ext_dir))
+        if (j == domain_jlo && ((bc.lo(1) == BCType::ext_dir) ||
+                                (bc.lo(1) == BCType::direction_dependent && vmac >= 0.0)) )
         {
             vmns = S(i,j-1,k,n);
             if ( n==YVEL && is_velocity )
@@ -153,7 +156,8 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n,
                     amrex_calc_yslope_extdir(i,j  ,k,n,order,S, extdir_or_ho_jlo, extdir_or_ho_jhi, domain_jlo, domain_jhi);
             }
         }
-        else if (j == domain_jhi+1 && (bc.hi(1) == BCType::ext_dir))
+        else if (j == domain_jhi+1 && ((bc.hi(1) == BCType::ext_dir) ||
+                                       (bc.lo(1) == BCType::direction_dependent && vmac <= 0.0)) )
         {
             vpls = S(i,j  ,k,n);
             if ( n==YVEL && is_velocity )
@@ -205,7 +209,8 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n,
 
         int order = 4;
 
-        if (k == domain_klo && (bc.lo(2) == BCType::ext_dir))
+        if (k == domain_klo && ((bc.lo(2) == BCType::ext_dir) ||
+                                (bc.lo(2) == BCType::direction_dependent && wmac >= 0.0)) )
         {
             wmns = S(i,j,k-1,n);
             if ( n == ZVEL && is_velocity )
@@ -218,7 +223,8 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n,
                     amrex_calc_zslope_extdir(i,j,k  ,n,order,S, extdir_or_ho_klo, extdir_or_ho_khi, domain_klo, domain_khi);
             }
         }
-        else if (k == domain_khi+1 && (bc.hi(2) == BCType::ext_dir))
+        else if (k == domain_khi+1 && ((bc.hi(2) == BCType::ext_dir) ||
+                                       (bc.lo(2) == BCType::direction_dependent && wmac <= 0.0)) )
         {
             wpls = S(i,j,k  ,n);
             if ( n == ZVEL && is_velocity )
diff --git a/Godunov/hydro_godunov_ppm.H b/Godunov/hydro_godunov_ppm.H
index 71ac154a8..f4efbaefa 100644
--- a/Godunov/hydro_godunov_ppm.H
+++ b/Godunov/hydro_godunov_ppm.H
@@ -10,6 +10,7 @@
 
 #include <AMReX_MultiFab.H>
 #include <AMReX_BCRec.H>
+#include <AMReX_BC_TYPES.H>
 #include <hydro_constants.H>
 
 namespace PPM {
@@ -281,12 +282,14 @@ void SetXBCs ( const int i, const int j, const int k, const int n,
                amrex::Real &sm, amrex::Real &sp,
                amrex::Real &sedge1, amrex::Real &sedge2,
                const amrex::Array4<const amrex::Real> &s,
+               const amrex::Real velm, const amrex::Real velp,
                const int bclo, const int bchi,
                const int domlo, const int domhi)
 {
     using namespace amrex;
 
-    if (bclo == BCType::ext_dir || bclo == BCType::hoextrap)
+    if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) ||
+         (bclo == BCType::direction_dependent && velm >= 0.0) )
     {
         if (i == domlo)
         {
@@ -326,7 +329,8 @@ void SetXBCs ( const int i, const int j, const int k, const int n,
         }
     }
 
-    if (bchi == BCType::ext_dir || bchi == BCType::hoextrap)
+    if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) ||
+         (bchi == BCType::direction_dependent && velp <= 0.0) )
     {
         if (i == domhi)
         {
@@ -373,12 +377,14 @@ void SetYBCs ( const int i, const int j, const int k, const int n,
                amrex::Real &sm, amrex::Real &sp,
                amrex::Real &sedge1, amrex::Real &sedge2,
                const amrex::Array4<const amrex::Real> &s,
+               const amrex::Real velm, const amrex::Real velp,
                const int bclo, const int bchi,
                const int domlo, const int domhi)
 {
     using namespace amrex;
 
-    if (bclo == BCType::ext_dir || bclo == BCType::hoextrap)
+    if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) ||
+         (bclo == BCType::direction_dependent && velm >= 0.0) )
     {
         if (j == domlo)
         {
@@ -418,7 +424,8 @@ void SetYBCs ( const int i, const int j, const int k, const int n,
         }
     }
 
-    if (bchi == BCType::ext_dir || bchi == BCType::hoextrap)
+    if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) ||
+         (bchi == BCType::direction_dependent && velp <= 0.0) )
     {
         if (j == domhi)
         {
@@ -465,12 +472,14 @@ void SetZBCs ( const int i, const int j, const int k, const int n,
                amrex::Real &sm, amrex::Real &sp,
                amrex::Real &sedge1, amrex::Real &sedge2,
                const amrex::Array4<const amrex::Real> &s,
+               const amrex::Real velm, const amrex::Real velp,
                const int bclo, const int bchi,
                const int domlo, const int domhi)
 {
     using namespace amrex;
 
-    if (bclo == BCType::ext_dir || bclo == BCType::hoextrap)
+    if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) ||
+         (bclo == BCType::direction_dependent && velm >= 0.0) )
     {
         if (k == domlo)
         {
@@ -510,7 +519,8 @@ void SetZBCs ( const int i, const int j, const int k, const int n,
         }
     }
 
-    if (bchi == BCType::ext_dir || bchi == BCType::hoextrap)
+    if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) ||
+         (bchi == BCType::direction_dependent && velp <= 0.0) )
     {
         if (k == domhi)
         {
@@ -582,7 +592,8 @@ void PredictVelOnXFace ( const int i, const int j, const int k, const int n,
     amrex::Real sm, sp;
     amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
 
-    SetXBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, bc.lo(0), bc.hi(0), domlo, domhi);
+    SetXBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad,
+                     bc.lo(0), bc.hi(0), domlo, domhi);
 
     amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp);
 
@@ -629,7 +640,8 @@ void PredictVelOnYFace ( const int i, const int j, const int k, const int n,
     amrex::Real sm, sp;
     amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
 
-    SetYBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, bc.lo(1), bc.hi(1), domlo, domhi);
+    SetYBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad,
+                     bc.lo(1), bc.hi(1), domlo, domhi);
 
     amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp);
 
@@ -677,7 +689,8 @@ void PredictVelOnZFace ( const int i, const int j, const int k, const int n,
     amrex::Real sm, sp;
     amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
 
-    SetZBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, bc.lo(2), bc.hi(2), domlo, domhi);
+    SetZBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad,
+                     bc.lo(2), bc.hi(2), domlo, domhi);
 
     amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp);
 
@@ -730,7 +743,8 @@ void PredictStateOnXFace ( const int i, const int j, const int k, const int n,
     amrex::Real sm, sp;
     amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
 
-    SetXBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, bc.lo(0), bc.hi(0), domlo, domhi);
+    SetXBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i+1,j,k),
+                     bc.lo(0), bc.hi(0), domlo, domhi);
 
     amrex::Real s6 = 6.0*s0 - 3.0*(sm + sp);
 
@@ -774,7 +788,8 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n,
     amrex::Real sm, sp;
     amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
 
-    SetYBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, bc.lo(1), bc.hi(1), domlo, domhi);
+    SetYBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i+1,j,k),
+                     bc.lo(1), bc.hi(1), domlo, domhi);
 
     amrex::Real s6 = 6.0*s0- 3.0*(sm + sp);
 
@@ -821,7 +836,8 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n,
     amrex::Real sm, sp;
     amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
 
-    SetZBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, bc.lo(2), bc.hi(2), domlo, domhi);
+    SetZBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k+1),
+                     bc.lo(2), bc.hi(2), domlo, domhi);
 
     amrex::Real s6 = 6.0*s0- 3.0*(sm + sp);
     amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j,k+1))*dt/dx;
@@ -913,9 +929,7 @@ void PredictStateOnFaces (amrex::Box const& bx,
 #endif
     });
 }
-
-
-}
+} // namespace
 
 #endif
 /** @} */
diff --git a/MOL/hydro_mol_edge_state_K.H b/MOL/hydro_mol_edge_state_K.H
index 6ff324dde..7adf93888 100644
--- a/MOL/hydro_mol_edge_state_K.H
+++ b/MOL/hydro_mol_edge_state_K.H
@@ -45,7 +45,8 @@ amrex::Real hydro_mol_xedge_state_extdir ( int i, int j, int k, int n,
         amrex::Real qmns = q(i-1,j,k,n) + 0.5 *
       amrex_calc_xslope_extdir( i-1, j, k, n, order, q, extdir_or_ho_lo, extdir_or_ho_hi, domlo, domhi );
 
-        HydroBC::SetXEdgeBCs(i,j,k,n,q,qmns,qpls,d_bcrec[n].lo(0),domlo,d_bcrec[n].hi(0),domhi,is_velocity);
+        HydroBC::SetXEdgeBCs(i,j,k,n,q,qmns,qpls,umac(i,j,k),umac(i,j,k),
+                             d_bcrec[n].lo(0),domlo,d_bcrec[n].hi(0),domhi,is_velocity);
 
         if ( (i==domlo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) )
         {
@@ -89,7 +90,8 @@ amrex::Real hydro_mol_xedge_state ( int i, int j, int k, int n,
     amrex::Real qpls = q(i  ,j,k,n) - 0.5 * amrex_calc_xslope( i  , j, k, n, order, q );
     amrex::Real qmns = q(i-1,j,k,n) + 0.5 * amrex_calc_xslope( i-1, j, k, n, order, q );
 
-    HydroBC::SetXEdgeBCs(i,j,k,n,q,qmns,qpls,d_bcrec[n].lo(0),domlo,d_bcrec[n].hi(0),domhi,is_velocity);
+    HydroBC::SetXEdgeBCs(i,j,k,n,q,qmns,qpls,umac(i,j,k),umac(i,j,k),
+                         d_bcrec[n].lo(0),domlo,d_bcrec[n].hi(0),domhi,is_velocity);
 
     if ( (i==domlo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) )
     {
@@ -152,7 +154,8 @@ amrex::Real hydro_mol_yedge_state_extdir ( int i, int j, int k, int n,
         amrex::Real qmns = q(i,j-1,k,n) + 0.5 *
       amrex_calc_yslope_extdir( i, j-1, k, n, order, q, extdir_or_ho_lo, extdir_or_ho_hi, domlo, domhi );
 
-        HydroBC::SetYEdgeBCs(i,j,k,n,q,qmns,qpls,d_bcrec[n].lo(1),domlo,d_bcrec[n].hi(1),domhi,is_velocity);
+        HydroBC::SetYEdgeBCs(i,j,k,n,q,qmns,qpls,vmac(i,j,k),vmac(i,j,k),
+                             d_bcrec[n].lo(1),domlo,d_bcrec[n].hi(1),domhi,is_velocity);
 
         if ( (j==domlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) )
         {
@@ -198,7 +201,8 @@ amrex::Real hydro_mol_yedge_state ( int i, int j, int k, int n,
     amrex::Real qpls = q(i,j  ,k,n) - 0.5 * amrex_calc_yslope( i, j  , k, n, order, q );
     amrex::Real qmns = q(i,j-1,k,n) + 0.5 * amrex_calc_yslope( i, j-1, k, n, order, q );
 
-    HydroBC::SetYEdgeBCs(i,j,k,n,q,qmns,qpls,d_bcrec[n].lo(1),domlo,d_bcrec[n].hi(1),domhi,is_velocity);
+    HydroBC::SetYEdgeBCs(i,j,k,n,q,qmns,qpls,vmac(i,j,k),vmac(i,j,k),
+                         d_bcrec[n].lo(1),domlo,d_bcrec[n].hi(1),domhi,is_velocity);
 
     if ( (j==domlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) )
     {
@@ -263,7 +267,8 @@ amrex::Real hydro_mol_zedge_state_extdir ( int i, int j, int k, int n,
         amrex::Real qmns = q(i,j,k-1,n) + 0.5 *
       amrex_calc_zslope_extdir( i, j, k-1, n, order, q, extdir_or_ho_lo, extdir_or_ho_hi, domlo, domhi );
 
-        HydroBC::SetZEdgeBCs(i,j,k,n,q,qmns,qpls,d_bcrec[n].lo(2),domlo,d_bcrec[n].hi(2),domhi,is_velocity);
+        HydroBC::SetZEdgeBCs(i,j,k,n,q,qmns,qpls,wmac(i,j,k),wmac(i,j,k),
+                             d_bcrec[n].lo(2),domlo,d_bcrec[n].hi(2),domhi,is_velocity);
 
         if ( (k==domlo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) )
         {
@@ -309,7 +314,8 @@ amrex::Real hydro_mol_zedge_state ( int i, int j, int k, int n,
     amrex::Real qpls = q(i,j,k  ,n) - 0.5 * amrex_calc_zslope( i, j, k  , n, order, q );
     amrex::Real qmns = q(i,j,k-1,n) + 0.5 * amrex_calc_zslope( i, j, k-1, n, order, q );
 
-    HydroBC::SetZEdgeBCs(i,j,k,n,q,qmns,qpls,d_bcrec[n].lo(2),domlo,d_bcrec[n].hi(2),domhi,is_velocity);
+    HydroBC::SetZEdgeBCs(i,j,k,n,q,qmns,qpls,wmac(i,j,k),wmac(i,j,k),
+                         d_bcrec[n].lo(2),domlo,d_bcrec[n].hi(2),domhi,is_velocity);
 
     if ( (k==domlo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) )
     {
diff --git a/MOL/hydro_mol_extrap_vel_to_faces_box.cpp b/MOL/hydro_mol_extrap_vel_to_faces_box.cpp
index 1355b2cd3..81aeb07e5 100644
--- a/MOL/hydro_mol_extrap_vel_to_faces_box.cpp
+++ b/MOL/hydro_mol_extrap_vel_to_faces_box.cpp
@@ -88,7 +88,8 @@ MOL::ExtrapVelToFacesBox (  AMREX_D_DECL( Box const& ubx,
             Real umns = vcc_mns + Real(0.5) * amrex_calc_xslope_extdir(
                  i-1,j,k,0,order,vcc,extdir_or_ho_ilo, extdir_or_ho_ihi, domain_ilo, domain_ihi);
 
-            HydroBC::SetXEdgeBCs(i, j, k, n, vcc, umns, upls, d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true);
+            HydroBC::SetXEdgeBCs(i, j, k, n, vcc, umns, upls, umns, upls,
+                                 d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true);
 
             if ( (i==domain_ilo) && (d_bcrec[0].lo(0) == BCType::foextrap || d_bcrec[0].lo(0) == BCType::hoextrap) )
             {
@@ -135,7 +136,8 @@ MOL::ExtrapVelToFacesBox (  AMREX_D_DECL( Box const& ubx,
             Real upls = vcc(i  ,j,k,0) - Real(0.5) * amrex_calc_xslope(i  ,j,k,0,order,vcc);
             Real umns = vcc(i-1,j,k,0) + Real(0.5) * amrex_calc_xslope(i-1,j,k,0,order,vcc);
 
-            HydroBC::SetXEdgeBCs(i, j, k, n, vcc, umns, upls, d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true);
+            HydroBC::SetXEdgeBCs(i, j, k, n, vcc, umns, upls, umns, upls,
+                                 d_bcrec[0].lo(0), domain_ilo, d_bcrec[0].hi(0), domain_ihi, true);
 
             if ( (i==domain_ilo) && (d_bcrec[0].lo(0) == BCType::foextrap || d_bcrec[0].lo(0) == BCType::hoextrap) )
             {
@@ -195,7 +197,8 @@ MOL::ExtrapVelToFacesBox (  AMREX_D_DECL( Box const& ubx,
             Real vmns = vcc_mns + Real(0.5) * amrex_calc_yslope_extdir(
                  i,j-1,k,1,order,vcc,extdir_or_ho_jlo,extdir_or_ho_jhi,domain_jlo,domain_jhi);
 
-            HydroBC::SetYEdgeBCs(i, j, k, n, vcc, vmns, vpls, d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true);
+            HydroBC::SetYEdgeBCs(i, j, k, n, vcc, vmns, vpls, vmns, vpls,
+                                 d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true);
 
             if ( (j==domain_jlo) && (d_bcrec[1].lo(1) == BCType::foextrap || d_bcrec[1].lo(1) == BCType::hoextrap) )
             {
@@ -242,7 +245,8 @@ MOL::ExtrapVelToFacesBox (  AMREX_D_DECL( Box const& ubx,
             Real vpls = vcc(i,j  ,k,1) - Real(0.5) * amrex_calc_yslope(i,j  ,k,1,order,vcc);
             Real vmns = vcc(i,j-1,k,1) + Real(0.5) * amrex_calc_yslope(i,j-1,k,1,order,vcc);
 
-            HydroBC::SetYEdgeBCs(i, j, k, n, vcc, vmns, vpls, d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true);
+            HydroBC::SetYEdgeBCs(i, j, k, n, vcc, vmns, vpls, vmns, vpls,
+                                 d_bcrec[1].lo(1), domain_jlo, d_bcrec[1].hi(1), domain_jhi, true);
 
             if ( (j==domain_jlo) && (d_bcrec[1].lo(1) == BCType::foextrap || d_bcrec[1].lo(1) == BCType::hoextrap) )
             {
@@ -302,7 +306,8 @@ MOL::ExtrapVelToFacesBox (  AMREX_D_DECL( Box const& ubx,
             Real wmns = vcc_mns + Real(0.5) * amrex_calc_zslope_extdir(
                  i,j,k-1,2,order,vcc,extdir_or_ho_klo,extdir_or_ho_khi,domain_klo,domain_khi);
 
-            HydroBC::SetZEdgeBCs(i, j, k, n, vcc, wmns, wpls, d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true);
+            HydroBC::SetZEdgeBCs(i, j, k, n, vcc, wmns, wpls, wmns, wpls,
+                                 d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true);
 
             if ( (k==domain_klo) && (d_bcrec[2].lo(2) == BCType::foextrap || d_bcrec[2].lo(2) == BCType::hoextrap) )
             {
@@ -349,7 +354,8 @@ MOL::ExtrapVelToFacesBox (  AMREX_D_DECL( Box const& ubx,
             Real wpls = vcc(i,j,k  ,2) - Real(0.5) * amrex_calc_zslope(i,j,k  ,2,order,vcc);
             Real wmns = vcc(i,j,k-1,2) + Real(0.5) * amrex_calc_zslope(i,j,k-1,2,order,vcc);
 
-            HydroBC::SetZEdgeBCs(i, j, k, n, vcc, wmns, wpls, d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true);
+            HydroBC::SetZEdgeBCs(i, j, k, n, vcc, wmns, wpls, wmns, wpls,
+                                 d_bcrec[2].lo(2), domain_klo, d_bcrec[2].hi(2), domain_khi, true);
 
             if ( (k==domain_klo) && (d_bcrec[2].lo(2) == BCType::foextrap || d_bcrec[2].lo(2) == BCType::hoextrap) )
             {
diff --git a/Utils/hydro_bcs_K.H b/Utils/hydro_bcs_K.H
index 050d6f38e..cd525dc9f 100644
--- a/Utils/hydro_bcs_K.H
+++ b/Utils/hydro_bcs_K.H
@@ -77,6 +77,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void SetXEdgeBCs( int i, int j, int k, int n,
                   const amrex::Array4<const amrex::Real> &s,
                   amrex::Real &lo, amrex::Real &hi,
+                  amrex::Real velm, amrex::Real velp,
                   int bclo, int domlo,
                   int bchi, int domhi,
                   bool is_velocity )
@@ -89,7 +90,8 @@ void SetXEdgeBCs( int i, int j, int k, int n,
 
     if (i == domlo)
     {
-        if (bclo == BCType::ext_dir)
+        if ( (bclo == BCType::ext_dir) ||
+             (bclo == BCType::direction_dependent && velp >= 0.0) )
         {
             lo = s(domlo-1, j, k, n);
             // For turbulent inflow, there are times when the inflow face
@@ -99,7 +101,8 @@ void SetXEdgeBCs( int i, int j, int k, int n,
             if ( n==XVEL && is_velocity ) hi=lo;
         }
         else if ( bclo == BCType::foextrap || bclo == BCType::hoextrap ||
-                  bclo == BCType::reflect_even )
+                  bclo == BCType::reflect_even ||
+                 (bclo == BCType::direction_dependent && velp < 0.0))
         {
             lo = hi;
         }
@@ -111,13 +114,15 @@ void SetXEdgeBCs( int i, int j, int k, int n,
     }
     else if (i == domhi+1)
     {
-        if (bchi == BCType::ext_dir)
+        if ( (bchi == BCType::ext_dir) ||
+             (bchi == BCType::direction_dependent && velm <= 0.0) )
         {
             hi = s(domhi+1, j, k, n);
             if ( n==XVEL && is_velocity ) lo = hi;
         }
         else if ( bchi == BCType::foextrap || bchi == BCType::hoextrap ||
-                  bchi == BCType::reflect_even )
+                  bchi == BCType::reflect_even ||
+                 (bchi == BCType::direction_dependent && velm > 0.0))
         {
             hi = lo;
         }
@@ -141,12 +146,13 @@ void SetXEdgeBCs( int i, int j, int k, int n,
  */
 
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void SetYEdgeBCs ( int i, int j, int k, int n,
-                   const amrex::Array4<const amrex::Real> &s,
-                   amrex::Real &lo, amrex::Real &hi,
-                   int bclo, int domlo,
-                   int bchi, int domhi,
-                   bool is_velocity )
+void SetYEdgeBCs (int i, int j, int k, int n,
+                  const amrex::Array4<const amrex::Real> &s,
+                  amrex::Real &lo, amrex::Real &hi,
+                  amrex::Real velm, amrex::Real velp,
+                  int bclo, int domlo,
+                  int bchi, int domhi,
+                  bool is_velocity )
 {
     using namespace amrex;
 
@@ -155,13 +161,15 @@ void SetYEdgeBCs ( int i, int j, int k, int n,
 
     if (j == domlo)
     {
-        if (bclo == BCType::ext_dir)
+        if ( (bclo == BCType::ext_dir) ||
+             (bclo == BCType::direction_dependent && velp >= 0.0) )
         {
             lo = s(i, domlo-1, k, n);
             if ( n==YVEL && is_velocity ) hi=lo;
         }
         else if ( bclo == BCType::foextrap || bclo == BCType::hoextrap ||
-                  bclo == BCType::reflect_even )
+                  bclo == BCType::reflect_even ||
+                 (bclo == BCType::direction_dependent && velp < 0.0))
         {
             lo = hi;
         }
@@ -173,13 +181,15 @@ void SetYEdgeBCs ( int i, int j, int k, int n,
     }
     else  if (j == domhi+1)
     {
-        if (bchi == BCType::ext_dir)
+        if ( (bchi == BCType::ext_dir) ||
+             (bchi == BCType::direction_dependent && velm <= 0.0) )
         {
             hi = s(i, domhi+1, k, n);
             if ( n==YVEL && is_velocity ) lo=hi;
         }
         else if ( bchi == BCType::foextrap || bchi == BCType::hoextrap ||
-                  bchi == BCType::reflect_even )
+                  bchi == BCType::reflect_even ||
+                 (bchi == BCType::direction_dependent && velm > 0.0))
         {
             hi = lo;
         }
@@ -208,12 +218,13 @@ void SetYEdgeBCs ( int i, int j, int k, int n,
  */
 
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void SetZEdgeBCs ( int i, int j, int k, int n,
-                   const amrex::Array4<const amrex::Real> &s,
-                   amrex::Real &lo, amrex::Real &hi,
-                   int bclo, int domlo,
-                   int bchi, int domhi,
-                   bool is_velocity )
+void SetZEdgeBCs (int i, int j, int k, int n,
+                  const amrex::Array4<const amrex::Real> &s,
+                  amrex::Real &lo, amrex::Real &hi,
+                  amrex::Real velm, amrex::Real velp,
+                  int bclo, int domlo,
+                  int bchi, int domhi,
+                  bool is_velocity )
 {
     using namespace amrex;
 
@@ -222,13 +233,15 @@ void SetZEdgeBCs ( int i, int j, int k, int n,
 
     if (k == domlo)
     {
-        if (bclo == BCType::ext_dir)
+        if ( (bclo == BCType::ext_dir) ||
+             (bclo == BCType::direction_dependent && velp >= 0.0) )
         {
             lo = s(i, j, domlo-1, n);
             if ( n==ZVEL && is_velocity ) hi=lo;
         }
         else if ( bclo == BCType::foextrap || bclo == BCType::hoextrap ||
-                  bclo == BCType::reflect_even )
+                  bclo == BCType::reflect_even ||
+                 (bclo == BCType::direction_dependent && velp < 0.0))
         {
             lo = hi;
         }
@@ -240,13 +253,15 @@ void SetZEdgeBCs ( int i, int j, int k, int n,
     }
     else if (k == domhi+1)
     {
-        if (bchi == BCType::ext_dir)
+        if ( (bchi == BCType::ext_dir) ||
+             (bchi == BCType::direction_dependent && velm <= 0.0) )
         {
             hi = s(i,j,domhi+1, n);
             if ( n==ZVEL && is_velocity ) lo=hi;
         }
         else if ( bchi == BCType::foextrap || bchi == BCType::hoextrap ||
-                  bchi == BCType::reflect_even )
+                  bchi == BCType::reflect_even ||
+                 (bchi == BCType::direction_dependent && velm > 0.0))
         {
             hi = lo;
         }

From cae9075d95b24c329e22f4161c8645c52c3499a6 Mon Sep 17 00:00:00 2001
From: Ann Almgren <asalmgren@lbl.gov>
Date: Mon, 29 Jul 2024 19:31:27 -0700
Subject: [PATCH 2/9] missed some ...

---
 EBGodunov/hydro_ebgodunov_edge_state_3D.cpp | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp b/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
index 501ffb0d6..fa93f4669 100644
--- a/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
+++ b/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
@@ -137,7 +137,6 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             Real uad = u_mac(i,j,k);
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-
             HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
             xlo(i,j,k,n) = lo;
@@ -304,6 +303,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             }
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
+            Real uad = u_mac(i,j,k);
             HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
             if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
@@ -441,6 +441,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             }
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
+            Real vad = v_mac(i,j,k);
             HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
             if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
@@ -575,6 +576,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             }
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
+            Real wad = w_mac(i,j,k);
             HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
             if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) )

From f0d5512b024a778f7c1ffc3c693e13a8fafc64b8 Mon Sep 17 00:00:00 2001
From: Ann Almgren <asalmgren@lbl.gov>
Date: Tue, 30 Jul 2024 11:52:45 -0700
Subject: [PATCH 3/9] fix typos

---
 Godunov/hydro_godunov_plm.H | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/Godunov/hydro_godunov_plm.H b/Godunov/hydro_godunov_plm.H
index b7e7fec90..fbe5a68af 100644
--- a/Godunov/hydro_godunov_plm.H
+++ b/Godunov/hydro_godunov_plm.H
@@ -90,7 +90,7 @@ void PredictStateOnXFace ( const int i, const int j, const int k, const int n,
 
         }
         else if (i == domain_ihi+1 && ((bc.hi(0) == BCType::ext_dir) ||
-                                       (bc.lo(0) == BCType::direction_dependent && umac <= 0.0)) )
+                                       (bc.hi(0) == BCType::direction_dependent && umac <= 0.0)) )
         {
             upls = S(i  ,j,k,n);
 
@@ -157,7 +157,7 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n,
             }
         }
         else if (j == domain_jhi+1 && ((bc.hi(1) == BCType::ext_dir) ||
-                                       (bc.lo(1) == BCType::direction_dependent && vmac <= 0.0)) )
+                                       (bc.hi(1) == BCType::direction_dependent && vmac <= 0.0)) )
         {
             vpls = S(i,j  ,k,n);
             if ( n==YVEL && is_velocity )
@@ -224,7 +224,7 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n,
             }
         }
         else if (k == domain_khi+1 && ((bc.hi(2) == BCType::ext_dir) ||
-                                       (bc.lo(2) == BCType::direction_dependent && wmac <= 0.0)) )
+                                       (bc.hi(2) == BCType::direction_dependent && wmac <= 0.0)) )
         {
             wpls = S(i,j,k  ,n);
             if ( n == ZVEL && is_velocity )

From 63f1982e95571e8b659e5f658c2ace0650bc0aaf Mon Sep 17 00:00:00 2001
From: Ann Almgren <asalmgren@lbl.gov>
Date: Tue, 30 Jul 2024 12:26:20 -0700
Subject: [PATCH 4/9] fix which velocity sent to *EdgeBCs routines

---
 .../hydro_godunov_extrap_vel_to_faces_2D.cpp    | 17 ++++++++++++-----
 .../hydro_godunov_extrap_vel_to_faces_3D.cpp    | 17 +++++++++++------
 2 files changed, 23 insertions(+), 11 deletions(-)

diff --git a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp
index 083a71cf5..da0e178e9 100644
--- a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp
+++ b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp
@@ -241,7 +241,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
         }
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, lo, hi, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        Real uad = u_ad(i,j,k);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         xlo(i,j,k,n) = lo;
         xhi(i,j,k,n) = hi;
@@ -282,7 +283,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
         l_yzhi = yhi(i,j,k,n);
 
         Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad,
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
@@ -304,7 +306,9 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
             sth += 0.5 * l_dt * f(i  ,j,k,n);
         }
 
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        Real uad = u_ad(i,j,k);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad,
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
         {
@@ -338,7 +342,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
         l_xzhi = xhi(i,j,k,n);
 
         Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad,
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
 
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
@@ -363,7 +368,9 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
            sth += 0.5 * l_dt * f(i,j  ,k,n);
         }
 
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        Real vad = v_ad(i,j,k);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad,
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
         {
diff --git a/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp b/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp
index d6a2019e5..c91f77dc0 100644
--- a/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp
+++ b/Godunov/hydro_godunov_extrap_vel_to_faces_3D.cpp
@@ -432,7 +432,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
             sth += Real(0.5) * l_dt * f(i  ,j,k,n);
         }
 
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, stl, sth, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        Real uad = u_ad(i,j,k);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
         {
@@ -522,7 +523,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
            sth += Real(0.5) * l_dt * f(i,j  ,k,n);
         }
 
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, stl, sth, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        Real vad = v_ad(i,j,k);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
         {
@@ -564,8 +566,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                               q, divu, v_ad, yedge);
 
         Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
-
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, uad, uad,
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         Real st = (uad >= 0.) ? l_xylo : l_xyhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
@@ -586,7 +588,8 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                               q, divu, u_ad, xedge);
 
         Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, vad, vad,
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
 
         Real st = (vad >= 0.) ? l_yxlo : l_yxhi;
@@ -613,7 +616,9 @@ Godunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
            sth += Real(0.5) * l_dt * f(i,j,k  ,n);
         }
 
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, stl, sth, bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
+        Real wad = w_ad(i,j,k);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, wad, wad,
+                             bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
 
         if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) )
         {

From fe5c283bea827ce8dc18aff425c0dd166eaca303 Mon Sep 17 00:00:00 2001
From: Ann Almgren <asalmgren@lbl.gov>
Date: Tue, 30 Jul 2024 17:14:36 -0700
Subject: [PATCH 5/9] more fixes

---
 EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp | 12 ++++++++----
 Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp     |  3 ++-
 2 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp
index e65524c35..831307075 100644
--- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp
+++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp
@@ -63,7 +63,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, lo, hi,
+            Real uad = u_ad(i,j,k);
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad,
                                  bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
             xlo(i,j,k,n) = lo;
@@ -76,7 +77,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, lo, hi,
+            Real vad = v_ad(i,j,k);
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad,
                                  bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
             ylo(i,j,k,n) = lo;
@@ -191,7 +193,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
             }
         }
 
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, stl, sth,
+        Real uad = u_ad(i,j,k);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad,
                              bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         // Prevent backflow
@@ -315,7 +318,8 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
             }
         }
 
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, stl, sth,
+        Real vad = v_ad(i,j,k);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad,
                              bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         // Prevent backflow
diff --git a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp
index da0e178e9..c41555118 100644
--- a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp
+++ b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp
@@ -259,7 +259,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
         }
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, lo, hi, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        Real vad = v_ad(i,j,k);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         ylo(i,j,k,n) = lo;
         yhi(i,j,k,n) = hi;

From d6abe9625f347f4dc2c5e9210d206d5922de1888 Mon Sep 17 00:00:00 2001
From: Ann Almgren <asalmgren@lbl.gov>
Date: Tue, 30 Jul 2024 17:25:10 -0700
Subject: [PATCH 6/9] fix more oops

---
 EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
index 01664568e..c9fc897ba 100644
--- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
+++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
@@ -267,7 +267,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
         }
         }
 
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, stl, sth,
+        Real uad = u_ad(i,j,k);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad,
                              bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         // Prevent backflow
@@ -422,7 +423,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
         }
         }
 
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, stl, sth,
+        Real vad = v_ad(i,j,k);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad,
                              bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         // Prevent backflow
@@ -580,7 +582,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
         }
         }
 
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, stl, sth,
+        Real wad = w_ad(i,j,k);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, wad, wad,
                              bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
 
 

From 493658599513b38db91a113fdcd5319a9743f974 Mon Sep 17 00:00:00 2001
From: Ann Almgren <asalmgren@lbl.gov>
Date: Wed, 31 Jul 2024 09:57:20 -0700
Subject: [PATCH 7/9] fix typo

---
 Godunov/hydro_godunov_ppm.H | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Godunov/hydro_godunov_ppm.H b/Godunov/hydro_godunov_ppm.H
index f4efbaefa..48e8deca4 100644
--- a/Godunov/hydro_godunov_ppm.H
+++ b/Godunov/hydro_godunov_ppm.H
@@ -788,7 +788,7 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n,
     amrex::Real sm, sp;
     amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
 
-    SetYBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i+1,j,k),
+    SetYBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j+1,k),
                      bc.lo(1), bc.hi(1), domlo, domhi);
 
     amrex::Real s6 = 6.0*s0- 3.0*(sm + sp);

From b1c4918f1ebaedf7f013ff2a0b463082a5f56fa3 Mon Sep 17 00:00:00 2001
From: Ann Almgren <asalmgren@lbl.gov>
Date: Tue, 6 Aug 2024 06:30:54 -0700
Subject: [PATCH 8/9] use advective vel where possible in BC upwinding

---
 EBGodunov/hydro_ebgodunov_edge_state_2D.cpp   | 24 ++++-----
 EBGodunov/hydro_ebgodunov_edge_state_3D.cpp   | 53 +++++++++++--------
 ...hydro_ebgodunov_extrap_vel_to_faces_2D.cpp | 20 +++----
 ...hydro_ebgodunov_extrap_vel_to_faces_3D.cpp | 36 ++++++-------
 EBMOL/hydro_ebmol_edge_state_K.H              | 12 ++---
 Godunov/hydro_godunov_edge_state_2D.cpp       | 23 ++++----
 Godunov/hydro_godunov_edge_state_3D.cpp       | 45 ++++++++--------
 .../hydro_godunov_extrap_vel_to_faces_2D.cpp  | 24 ++++-----
 Utils/hydro_bcs_K.H                           | 10 +++-
 9 files changed, 125 insertions(+), 122 deletions(-)

diff --git a/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp b/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp
index f19b07837..0fe0a44dc 100644
--- a/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp
+++ b/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp
@@ -109,8 +109,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            Real uad = u_mac(i,j,k);
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, u_mac(i,j,k), u_mac(i,j,k),
+                                 bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
             xlo(i,j,k,n) = lo;
             xhi(i,j,k,n) = hi;
@@ -122,8 +122,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            Real vad = v_mac(i,j,k);
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, v_mac(i,j,k), v_mac(i,j,k),
+                                 bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
             ylo(i,j,k,n) = lo;
             yhi(i,j,k,n) = hi;
@@ -144,8 +144,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
 
             l_yzlo = ylo(i,j,k,n);
             l_yzhi = yhi(i,j,k,n);
-            Real vad = v_mac(i,j,k);
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, v_mac(i,j,k), v_mac(i,j,k),
+                                 bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
             Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
             Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
@@ -222,8 +222,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             }
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            Real uad = u_mac(i,j,k);
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, u_mac(i,j,k), u_mac(i,j,k),
+                                 bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
             if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
             {
@@ -261,8 +261,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             l_xzlo = xlo(i,j,k,n);
             l_xzhi = xhi(i,j,k,n);
 
-            Real uad = u_mac(i,j,k);
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, u_mac(i,j,k), u_mac(i,j,k),
+                                 bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
             Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
             Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
@@ -339,8 +339,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             }
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            Real vad = v_mac(i,j,k);
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, v_mac(i,j,k), v_mac(i,j,k),
+                                 bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
             if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
             {
diff --git a/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp b/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
index fa93f4669..eda4206d8 100644
--- a/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
+++ b/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
@@ -137,7 +137,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             Real uad = u_mac(i,j,k);
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, u_mac(i,j,k), u_mac(i,j,k),
+                                 bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
             xlo(i,j,k,n) = lo;
             xhi(i,j,k,n) = hi;
@@ -155,7 +156,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, v_mac(i,j,k), v_mac(i,j,k),
+                                 bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
             ylo(i,j,k,n) = lo;
             yhi(i,j,k,n) = hi;
@@ -169,15 +171,15 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             Real lo = Ipz(i,j,k-1,n);
             Real hi = Imz(i,j,k  ,n);
 
-            Real wad = w_mac(i,j,k);
-
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+            HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, w_mac(i,j,k), w_mac(i,j,k),
+                                 bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
             zlo(i,j,k,n) = lo;
             zhi(i,j,k,n) = hi;
 
+            Real wad = w_mac(i,j,k);
             Real st = (wad >= 0.) ? lo : hi;
             Real fuz = (amrex::Math::abs(wad) < small_vel) ? 0. : 1.;
             Imz(i,j,k,n) = fuz*st + (Real(1.0) - fuz)*Real(0.5)*(hi + lo);
@@ -206,9 +208,10 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
                                    zlo(i,j,k,n), zhi(i,j,k,n),
                                    q, divu, apx, apy, apz, vfrac_arr, v_mac, yed);
 
-        Real wad = w_mac(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, w_mac(i,j,k), w_mac(i,j,k),
+                             bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
+        Real wad = w_mac(i,j,k);
         Real st = (wad >= 0.) ? l_zylo : l_zyhi;
         Real fu = (amrex::Math::abs(wad) < small_vel) ? Real(0.0) : Real(1.0);
         zylo(i,j,k,n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_zyhi + l_zylo);
@@ -223,9 +226,10 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
                                    ylo(i,j,k,n), yhi(i,j,k,n),
                                    q, divu, apx, apy, apz, vfrac_arr, w_mac, zed);
 
-        Real vad = v_mac(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, v_mac(i,j,k), v_mac(i,j,k),
+                                 bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
+        Real vad = v_mac(i,j,k);
         Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? Real(0.0) : Real(1.0);
         yzlo(i,j,k,n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_yzhi + l_yzlo);
@@ -303,8 +307,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             }
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            Real uad = u_mac(i,j,k);
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, u_mac(i,j,k), u_mac(i,j,k),
+                                 bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
             if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
             {
@@ -344,8 +348,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
                                    xlo(i,j,k,n),  xhi(i,j,k,n),
                                    q, divu, apx, apy, apz, vfrac_arr, w_mac, zed);
 
-        Real uad = u_mac(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, u_mac(i,j,k), u_mac(i,j,k),
+                                 bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
@@ -361,9 +365,10 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
                                    zlo(i,j,k,n), zhi(i,j,k,n),
                                    q, divu, apx, apy, apz, vfrac_arr, u_mac, xed);
 
-        Real wad = w_mac(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, w_mac(i,j,k), w_mac(i,j,k),
+                             bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
+        Real wad = w_mac(i,j,k);
         Real st = (wad >= 0.) ? l_zxlo : l_zxhi;
         Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0;
         zxlo(i,j,k,n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_zxhi + l_zxlo);
@@ -441,8 +446,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             }
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            Real vad = v_mac(i,j,k);
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, v_mac(i,j,k), v_mac(i,j,k),
+                                 bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
             if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
             {
@@ -480,9 +485,10 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
                                    xlo(i,j,k,n), xhi(i,j,k,n),
                                    q, divu, apx, apy, apz, vfrac_arr, v_mac, yed);
 
-        Real uad = u_mac(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, u_mac(i,j,k), u_mac(i,j,k),
+                                 bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
+        Real uad = u_mac(i,j,k);
         Real st = (uad >= 0.) ? l_xylo : l_xyhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
         xylo(i,j,k,n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_xyhi + l_xylo);
@@ -497,9 +503,10 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
                                    ylo(i,j,k,n), yhi(i,j,k,n),
                                    q, divu, apx, apy, apz, vfrac_arr, u_mac, xed);
 
-        Real vad = v_mac(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, v_mac(i,j,k), v_mac(i,j,k),
+                                 bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
+        Real vad = v_mac(i,j,k);
         Real st = (vad >= 0.) ? l_yxlo : l_yxhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
         yxlo(i,j,k,n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_yxhi + l_yxlo);
@@ -576,8 +583,8 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             }
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            Real wad = w_mac(i,j,k);
-            HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+            HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, w_mac(i,j,k), w_mac(i,j,k),
+                                 bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
             if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) )
             {
diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp
index 831307075..d91acdd96 100644
--- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp
+++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_2D.cpp
@@ -63,8 +63,7 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            Real uad = u_ad(i,j,k);
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad,
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, u_ad(i,j,k), u_ad(i,j,k),
                                  bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
             xlo(i,j,k,n) = lo;
@@ -77,8 +76,7 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            Real vad = v_ad(i,j,k);
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad,
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, v_ad(i,j,k), v_ad(i,j,k),
                                  bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
             ylo(i,j,k,n) = lo;
@@ -106,10 +104,10 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
 
             l_yzlo = ylo(i,j,k,n);
             l_yzhi = yhi(i,j,k,n);
-            Real vad = v_ad(i,j,k);
-            HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad,
+            HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, v_ad(i,j,k), v_ad(i,j,k),
                                  bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
+            Real vad = v_ad(i,j,k);
             Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
             Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
             yhat(i,j,k) = fu*st + (1.0 - fu) * 0.5 * (l_yzhi + l_yzlo);
@@ -193,8 +191,7 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
             }
         }
 
-        Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad,
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, u_ad(i,j,k), u_ad(i,j,k),
                              bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         // Prevent backflow
@@ -235,10 +232,10 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
             l_xzlo = xlo(i,j,k,n);
             l_xzhi = xhi(i,j,k,n);
 
-            Real uad = u_ad(i,j,k);
-            HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad,
+            HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, u_ad(i,j,k), u_ad(i,j,k),
                                  bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
+            Real uad = u_ad(i,j,k);
             Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
             Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
             xhat(i,j,k) = fu*st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo);
@@ -318,8 +315,7 @@ EBGodunov::ExtrapVelToFacesOnBox (Box const& /*bx*/, int ncomp,
             }
         }
 
-        Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad,
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, v_ad(i,j,k), v_ad(i,j,k),
                              bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         // Prevent backflow
diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
index c9fc897ba..336de95f0 100644
--- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
+++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
@@ -108,7 +108,6 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
             Real lo = Ipz(i,j,k-1,n);
             Real hi = Imz(i,j,k  ,n);
 
-            Real wad = w_ad(i,j,k);
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
             EBGodunovBC::SetZBCs(i, j, k, n, q, lo, hi, bc.lo(2), bc.hi(2), dlo.z, dhi.z, true);
@@ -116,6 +115,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
             zlo(i,j,k,n) = lo;
             zhi(i,j,k,n) = hi;
 
+            Real wad = w_ad(i,j,k);
             Real st = (wad >= 0.) ? lo : hi;
             Real fu = (amrex::Math::abs(wad) < small_vel) ? Real(0.0) : Real(1.0);
             Imz(i, j, k, n) = fu*st + (Real(1.0) - fu)*Real(0.5)*(hi + lo); // store zedge
@@ -163,11 +163,11 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                                    zlo(i,j,k,n), zhi(i,j,k,n),
                                    q, divu, apx, apy, apz, vfrac_arr, v_ad, yedge);
 
-        Real wad = w_ad(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, wad, wad,
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, w_ad(i,j,k), w_ad(i,j,k),
                              bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
 
 
+        Real wad = w_ad(i,j,k);
         Real st = (wad >= 0.) ? l_zylo : l_zyhi;
         Real fu = (amrex::Math::abs(wad) < small_vel) ? Real(0.0) : Real(1.0);
         zylo(i,j,k) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_zyhi + l_zylo);
@@ -182,10 +182,10 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                                    ylo(i,j,k,n), yhi(i,j,k,n),
                                    q, divu, apx, apy, apz, vfrac_arr, w_ad, zedge);
 
-        Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad,
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, v_ad(i,j,k), v_ad(i,j,k),
                              bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
+        Real vad = v_ad(i,j,k);
         Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? Real(0.0) : Real(1.0);
         yzlo(i,j,k) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_yzhi + l_yzlo);
@@ -267,8 +267,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
         }
         }
 
-        Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad,
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, u_ad(i,j,k), u_ad(i,j,k),
                              bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         // Prevent backflow
@@ -318,11 +317,10 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                                    xlo(i,j,k,n),  xhi(i,j,k,n),
                                    q, divu, apx, apy, apz, vfrac_arr, w_ad, zedge);
 
-        Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad,
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, u_ad(i,j,k), u_ad(i,j,k),
                              bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
-
+        Real uad = u_ad(i,j,k);
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : Real(1.0);
         xzlo(i,j,k) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_xzhi + l_xzlo);
@@ -337,11 +335,10 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                                    zlo(i,j,k,n), zhi(i,j,k,n),
                                    q, divu, apx, apy, apz, vfrac_arr, u_ad, xedge);
 
-        Real wad = w_ad(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, wad, wad,
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, w_ad(i,j,k), w_ad(i,j,k),
                             bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
 
-
+        Real wad = w_ad(i,j,k);
         Real st = (wad >= 0.) ? l_zxlo : l_zxhi;
         Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : Real(1.0);
         zxlo(i,j,k) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_zxhi + l_zxlo);
@@ -423,8 +420,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
         }
         }
 
-        Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad,
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, v_ad(i,j,k), v_ad(i,j,k),
                              bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         // Prevent backflow
@@ -475,7 +471,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                                    q, divu, apx, apy, apz, vfrac_arr, v_ad, yedge);
 
         Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, uad, uad,
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, u_ad(i,j,k), u_ad(i,j,k),
                              bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
 
@@ -497,11 +493,10 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
                                    ylo(i,j,k,n), yhi(i,j,k,n),
                                    q, divu, apx, apy, apz, vfrac_arr, u_ad, xedge);
 
-        Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, vad, vad,
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, v_ad(i,j,k), v_ad(i,j,k),
                              bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
-
+        Real vad = v_ad(i,j,k);
         Real st = (vad >= 0.) ? l_yxlo : l_yxhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : Real(1.0);
         yxlo(i,j,k) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_yxhi + l_yxlo);
@@ -582,8 +577,7 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
         }
         }
 
-        Real wad = w_ad(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, wad, wad,
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, w_ad(i,j,k), w_ad(i,j,k),
                              bc.lo(2), dlo.z, bc.hi(2), dhi.z, true);
 
 
diff --git a/EBMOL/hydro_ebmol_edge_state_K.H b/EBMOL/hydro_ebmol_edge_state_K.H
index fa42b27a3..19a9f33a1 100644
--- a/EBMOL/hydro_ebmol_edge_state_K.H
+++ b/EBMOL/hydro_ebmol_edge_state_K.H
@@ -141,7 +141,7 @@ amrex::Real hydro_ebmol_xedge_state_extdir ( AMREX_D_DECL(int i, int j, int k),
 
         qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
 
-        HydroBC::SetXEdgeBCs(i, j, k, 0, q, qmns, qpls, qmns, qpls,
+        HydroBC::SetXEdgeBCs(i, j, k, 0, q, qmns, qpls, umac(i,j,k), umac(i,j,k),
                              d_bcrec[n].lo(0), domain_ilo, d_bcrec[n].hi(0), domain_ihi, is_velocity);
 
         if ( (i==domain_ilo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) )
@@ -250,7 +250,7 @@ amrex::Real hydro_ebmol_xedge_state ( AMREX_D_DECL(int i, int j, int k), int n,
 
     qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
 
-    HydroBC::SetXEdgeBCs(i, j, k, 0, q, qmns, qpls, qmns, qpls,
+    HydroBC::SetXEdgeBCs(i, j, k, 0, q, qmns, qpls, umac(i,j,k), umac(i,j,k),
                          d_bcrec[n].lo(0), domain_ilo, d_bcrec[n].hi(0), domain_ihi, is_velocity);
 
     if ( (i==domain_ilo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) )
@@ -405,7 +405,7 @@ amrex::Real hydro_ebmol_yedge_state_extdir ( AMREX_D_DECL(int i, int j, int k),
                                         + delta_y * slopes_eb_lo[1];
 #endif
 
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, qmns, qpls, qmns, qpls,
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, qmns, qpls, vmac(i,j,k), vmac(i,j,k),
                              d_bcrec[n].lo(1), domain_jlo, d_bcrec[n].hi(1), domain_jhi, is_velocity);
 
         if ( (j==domain_jlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) )
@@ -518,7 +518,7 @@ amrex::Real hydro_ebmol_yedge_state ( AMREX_D_DECL(int i, int j, int k), int n,
 
     qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
 
-    HydroBC::SetYEdgeBCs(i, j, k, n, q, qmns, qpls, qmns, qpls,
+    HydroBC::SetYEdgeBCs(i, j, k, n, q, qmns, qpls, vmac(i,j,k), vmac(i,j,k),
                          d_bcrec[n].lo(1), domain_jlo, d_bcrec[n].hi(1), domain_jhi, is_velocity);
 
     if ( (j==domain_jlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) )
@@ -662,7 +662,7 @@ amrex::Real hydro_ebmol_zedge_state_extdir ( int i, int j, int k, int n,
 
         qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
 
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, qmns, qpls, qmns, qpls,
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, qmns, qpls, wmac(i,j,k), wmac(i,j,k),
                              d_bcrec[n].lo(2), domain_klo, d_bcrec[n].hi(2), domain_khi, is_velocity);
 
         if ( (k==domain_klo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) )
@@ -757,7 +757,7 @@ amrex::Real hydro_ebmol_zedge_state ( int i, int j, int k, int n,
 
     qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
 
-    HydroBC::SetZEdgeBCs(i, j, k, n, q, qmns, qpls, qmns, qpls,
+    HydroBC::SetZEdgeBCs(i, j, k, n, q, qmns, qpls, wmac(i,j,k), wmac(i,j,k),
                          d_bcrec[n].lo(2), domain_klo, d_bcrec[n].hi(2), domain_khi, is_velocity);
 
     if ( (k==domain_klo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) )
diff --git a/Godunov/hydro_godunov_edge_state_2D.cpp b/Godunov/hydro_godunov_edge_state_2D.cpp
index cf735be23..3daf6bd59 100644
--- a/Godunov/hydro_godunov_edge_state_2D.cpp
+++ b/Godunov/hydro_godunov_edge_state_2D.cpp
@@ -139,8 +139,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-        Real uad = umac(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, umac(i,j,k), umac(i,j,k),
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
         xlo(i,j,k,n) = lo;
         xhi(i,j,k,n) = hi;
     },
@@ -157,8 +157,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-        Real vad = vmac(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vmac(i,j,k), vmac(i,j,k),
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         ylo(i,j,k,n) = lo;
         yhi(i,j,k,n) = hi;
@@ -180,7 +180,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
         l_yzlo = ylo(i,j,k,n);
         l_yzhi = yhi(i,j,k,n);
         Real vad = vmac(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vmac(i,j,k), vmac(i,j,k),
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
@@ -229,8 +230,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        Real uad = umac(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, umac(i,j,k), umac(i,j,k),
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
         if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
         {
@@ -263,8 +264,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
         l_xzlo = xlo(i,j,k,n);
         l_xzhi = xhi(i,j,k,n);
 
-        Real uad = umac(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, umac(i,j,k), umac(i,j,k),
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
@@ -314,8 +315,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        Real vad = vmac(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vmac(i,j,k), vmac(i,j,k),
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
         {
diff --git a/Godunov/hydro_godunov_edge_state_3D.cpp b/Godunov/hydro_godunov_edge_state_3D.cpp
index fcccd6b91..303a20859 100644
--- a/Godunov/hydro_godunov_edge_state_3D.cpp
+++ b/Godunov/hydro_godunov_edge_state_3D.cpp
@@ -161,7 +161,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, umac(i,j,k), umac(i,j,k),
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
         xlo(i,j,k,n) = lo;
         xhi(i,j,k,n) = hi;
         Real st = (uval) ? lo : hi;
@@ -184,7 +185,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vmac(i,j,k), vmac(i,j,k),
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         ylo(i,j,k,n) = lo;
         yhi(i,j,k,n) = hi;
@@ -207,7 +209,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, lo, hi, wmac(i,j,k), wmac(i,j,k),
+                             bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
         zlo(i,j,k,n) = lo;
         zhi(i,j,k,n) = hi;
@@ -235,8 +238,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
                               zlo(i,j,k,n), zhi(i,j,k,n),
                               q, divu, vmac, Imy);
 
-        Real wad = wmac(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, wmac(i,j,k), wmac(i,j,k),
+                             bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
         Real st = (wad >= 0.) ? l_zylo : l_zyhi;
         Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0;
@@ -252,8 +255,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
                               ylo(i,j,k,n), yhi(i,j,k,n),
                               q, divu, wmac, Imz);
 
-        Real vad = vmac(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vmac(i,j,k), vmac(i,j,k),
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
@@ -301,8 +304,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        Real uad = umac(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, umac(i,j,k), umac(i,j,k),
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
         if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
         {
@@ -337,8 +340,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
                               xlo(i,j,k,n),  xhi(i,j,k,n),
                               q, divu, wmac, Imz);
 
-        Real uad = umac(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, umac(i,j,k), umac(i,j,k),
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
@@ -354,8 +357,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
                               zlo(i,j,k,n), zhi(i,j,k,n),
                               q, divu, umac, Imx);
 
-        Real wad = wmac(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, wmac(i,j,k), wmac(i,j,k),
+                             bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
         Real st = (wad >= 0.) ? l_zxlo : l_zxhi;
         Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0;
@@ -401,8 +404,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        Real vad = vmac(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vmac(i,j,k), vmac(i,j,k),
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
         {
@@ -437,8 +440,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
                               xlo(i,j,k,n), xhi(i,j,k,n),
                               q, divu, vmac, Imy);
 
-        Real uad = umac(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, umac(i,j,k), umac(i,j,k),
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
         Real st = (uad >= 0.) ? l_xylo : l_xyhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
@@ -454,8 +457,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
                               ylo(i,j,k,n), yhi(i,j,k,n),
                               q, divu, umac, Imx);
 
-        Real vad = vmac(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, vmac(i,j,k), vmac(i,j,k),
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
         Real st = (vad >= 0.) ? l_yxlo : l_yxhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
@@ -502,8 +505,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
 
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        Real wad = wmac(i,j,k);
-        HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, wad, wad, bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
+        HydroBC::SetZEdgeBCs(i, j, k, n, q, stl, sth, wmac(i,j,k), wmac(i,j,k),
+                             bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
         if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) )
         {
diff --git a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp
index c41555118..3e3ca45dd 100644
--- a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp
+++ b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp
@@ -241,8 +241,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
         }
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, uad, uad, bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, lo, hi, u_ad(i,j,k), u_ad(i,j,k),
+                             bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         xlo(i,j,k,n) = lo;
         xhi(i,j,k,n) = hi;
@@ -259,8 +259,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
         }
 
         const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-        Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, vad, vad, bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, lo, hi, v_ad(i,j,k), v_ad(i,j,k),
+                             bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         ylo(i,j,k,n) = lo;
         yhi(i,j,k,n) = hi;
@@ -283,10 +283,10 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
         l_yzlo = ylo(i,j,k,n);
         l_yzhi = yhi(i,j,k,n);
 
-        Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vad, vad,
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, v_ad(i,j,k), v_ad(i,j,k),
                              bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
+        Real vad = v_ad(i,j,k);
         Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
         yzlo(i,j,k) = fu*st + (1.0 - fu) * 0.5 * (l_yzhi + l_yzlo);
@@ -307,8 +307,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
             sth += 0.5 * l_dt * f(i  ,j,k,n);
         }
 
-        Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, uad, uad,
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, stl, sth, u_ad(i,j,k), u_ad(i,j,k),
                              bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
         if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) )
@@ -342,17 +341,15 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
         l_xzlo = xlo(i,j,k,n);
         l_xzhi = xhi(i,j,k,n);
 
-        Real uad = u_ad(i,j,k);
-        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, uad, uad,
+        HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, u_ad(i,j,k), u_ad(i,j,k),
                              bc.lo(0), dlo.x, bc.hi(0), dhi.x, true);
 
-
+        Real uad = u_ad(i,j,k);
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
         xzlo(i,j,k) = fu*st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo);
     });
 
-
     amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
     {
         constexpr int n = 1;
@@ -369,8 +366,7 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp,
            sth += 0.5 * l_dt * f(i,j  ,k,n);
         }
 
-        Real vad = v_ad(i,j,k);
-        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, vad, vad,
+        HydroBC::SetYEdgeBCs(i, j, k, n, q, stl, sth, v_ad(i,j,k), v_ad(i,j,k),
                              bc.lo(1), dlo.y, bc.hi(1), dhi.y, true);
 
         if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) )
diff --git a/Utils/hydro_bcs_K.H b/Utils/hydro_bcs_K.H
index cd525dc9f..b0bdd8160 100644
--- a/Utils/hydro_bcs_K.H
+++ b/Utils/hydro_bcs_K.H
@@ -90,8 +90,14 @@ void SetXEdgeBCs( int i, int j, int k, int n,
 
     if (i == domlo)
     {
-        if ( (bclo == BCType::ext_dir) ||
-             (bclo == BCType::direction_dependent && velp >= 0.0) )
+        if ( (bclo == BCType::direction_dependent) && is_velocity &&
+                  (n == XVEL) && (s(domlo-1,j,k,n) >= 0.0) )
+        {
+            lo = s(domlo-1, j, k, n);
+            if ( n==XVEL && is_velocity ) hi=lo;
+        }
+        else if ( (bclo == BCType::ext_dir) ||
+             (bclo == BCType::direction_dependent && velm >= 0.0) ) {
         {
             lo = s(domlo-1, j, k, n);
             // For turbulent inflow, there are times when the inflow face

From ecd625fbfce42ffd3e9fa716dfa8d1312419af64 Mon Sep 17 00:00:00 2001
From: Ann Almgren <asalmgren@lbl.gov>
Date: Tue, 6 Aug 2024 08:07:00 -0700
Subject: [PATCH 9/9] more changes to support direction_dependent bc

---
 EBGodunov/hydro_ebgodunov_bcs_K.H             | 118 ++++++++++++------
 EBGodunov/hydro_ebgodunov_edge_state_2D.cpp   |   2 +
 EBGodunov/hydro_ebgodunov_edge_state_3D.cpp   |   1 +
 .../hydro_ebgodunov_extrap_vel_to_faces.cpp   |   9 +-
 ...hydro_ebgodunov_extrap_vel_to_faces_3D.cpp |  13 +-
 Godunov/hydro_godunov_edge_state_2D.cpp       |   1 +
 Godunov/hydro_godunov_edge_state_3D.cpp       |   6 +
 Godunov/hydro_godunov_ppm.H                   |  10 +-
 Utils/hydro_bcs_K.H                           |  89 +++++++++----
 9 files changed, 174 insertions(+), 75 deletions(-)

diff --git a/EBGodunov/hydro_ebgodunov_bcs_K.H b/EBGodunov/hydro_ebgodunov_bcs_K.H
index 64747eafb..d2fe7e67c 100644
--- a/EBGodunov/hydro_ebgodunov_bcs_K.H
+++ b/EBGodunov/hydro_ebgodunov_bcs_K.H
@@ -25,8 +25,8 @@ namespace EBGodunovBC {
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void SetXBCs (const int i, const int j, const int k, const int n,
               const amrex::Array4<const amrex::Real> &s,
-              amrex::Real &lo,
-              amrex::Real &hi,
+              amrex::Real &lo, amrex::Real &hi,
+              amrex::Real velm, amrex::Real velp,
               const int bclo, const int bchi,
               const int domlo, const int domhi,
               const bool is_velocity )
@@ -37,7 +37,14 @@ void SetXBCs (const int i, const int j, const int k, const int n,
     // Low X
     if (i <= domlo)
     {
-        if (bclo==BCType::ext_dir)
+        if ( (bclo == BCType::direction_dependent) && is_velocity &&
+                (n == XVEL) && (s(domlo-1,j,k,n) >= 0.0) )
+        {
+            lo = s(domlo-1, j, k, n);
+            hi = s(domlo-1, j, k, n);
+        }
+        else if ( bclo == BCType::ext_dir ||
+                 (bclo == BCType::direction_dependent && velp >= 0.0) )
         {
             lo = s(domlo-1,j,k,n);
             // For turbulent inflow, there are times when the inflow face
@@ -46,7 +53,8 @@ void SetXBCs (const int i, const int j, const int k, const int n,
             // tangential components to transport values from the interior.
             if( n == XVEL && is_velocity ) hi=lo;
         }
-        else if (bclo == BCType::foextrap || bclo == BCType::hoextrap)
+        else if ( bclo == BCType::foextrap || bclo == BCType::hoextrap ||
+                 (bclo == BCType::direction_dependent && velp < 0.0) )
         {
             lo = hi;
         } else if (bclo == BCType::reflect_even)
@@ -71,8 +79,8 @@ void SetXBCs (const int i, const int j, const int k, const int n,
         else if (bclo == BCType::reflect_odd)
         {
             if ( i==domlo ) {
-                hi = 0.;
-                lo = 0.;
+                hi = Real(0.);
+                lo = Real(0.);
             } else {
                 Abort("EBGodunovBC::SetXBCs not yet fully implemented for reflect_odd BC. See comments in EBGodunovBC::SetXBCs.");
             //     lo = -hi_arr(2*domlo-i  ,j,k,n);
@@ -83,12 +91,20 @@ void SetXBCs (const int i, const int j, const int k, const int n,
     // High X
     else if (i > domhi)
     {
-        if (bchi==BCType::ext_dir)
+        if ( (bchi == BCType::direction_dependent) && is_velocity &&
+                (n == XVEL) && (s(domhi+1,j,k,n) <= 0.0) )
+        {
+            hi = s(domhi+1, j, k, n);
+            lo = s(domhi+1, j, k, n);
+        }
+        else if ( (bchi == BCType::ext_dir) ||
+             (bchi == BCType::direction_dependent && velm <= 0.0) )
         {
             hi = s(domhi+1,j,k,n) ;
             if( n ==XVEL && is_velocity ) lo=hi;
         }
-        else if (bchi == BCType::foextrap || bchi == BCType::hoextrap )
+        else if ( bchi == BCType::foextrap || bchi == BCType::hoextrap ||
+                 (bchi== BCType::direction_dependent && velm > 0.0) )
         {
             hi = lo;
         }
@@ -104,8 +120,8 @@ void SetXBCs (const int i, const int j, const int k, const int n,
         else if (bchi == BCType::reflect_odd)
         {
             if ( i==domhi+1 ) {
-                hi = 0.;
-                lo = 0.;
+                hi = Real(0.);
+                lo = Real(0.);
             } else {
                 Abort("EBGodunovBC::SetXBCs not yet fully implemented for reflect_odd BC. See comments in EBGodunovBC::SetXBCs.");
             //     hi = -lo_arr(2*(domhi+1)-i-1,j,k,n);
@@ -120,8 +136,8 @@ void SetXBCs (const int i, const int j, const int k, const int n,
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void SetYBCs (const int i, const int j, const int k, const int n,
               const amrex::Array4<const amrex::Real> &s,
-              amrex::Real &lo,
-              amrex::Real &hi,
+              amrex::Real &lo, amrex::Real &hi,
+              amrex::Real velm, amrex::Real velp,
               const int bclo, const int bchi,
               const int domlo, const int domhi,
               const bool is_velocity )
@@ -132,12 +148,20 @@ void SetYBCs (const int i, const int j, const int k, const int n,
     // Low Y
     if (j <= domlo)
     {
-        if (bclo==BCType::ext_dir)
+        if ( (bclo == BCType::direction_dependent) && is_velocity &&
+                (n == YVEL) && (s(i,domlo-1,k,n) >= 0.0) )
+        {
+            lo = s(i, domlo-1, k, n);
+            hi = s(i, domlo-1, k, n);
+        }
+        else if ( bclo == BCType::ext_dir ||
+                 (bclo == BCType::direction_dependent && velp >= 0.0) )
         {
             lo = s(i,domlo-1,k,n);
             if ( n == YVEL && is_velocity ) hi = lo;
         }
-        else if (bclo == BCType::foextrap || bclo == BCType::hoextrap)
+        else if ( bclo == BCType::foextrap || bclo == BCType::hoextrap ||
+                 (bclo == BCType::direction_dependent && velp < 0.0) )
         {
             lo = hi;
         }
@@ -153,8 +177,8 @@ void SetYBCs (const int i, const int j, const int k, const int n,
         else if(bclo == BCType::reflect_odd)
         {
             if ( j==domlo ) {
-                hi = 0.;
-                lo = 0.;
+                hi = Real(0.);
+                lo = Real(0.);
             } else {
                 Abort("EBGodunovBC::SetYBCs not yet fully implemented for reflect_odd BC. See comments in EBGodunovBC::SetXBCs.");
             //     lo = -hi_arr(i,2*domlo-j  ,k,n);
@@ -165,12 +189,20 @@ void SetYBCs (const int i, const int j, const int k, const int n,
     // High Y
     else if (j > domhi)
     {
-        if (bchi==BCType::ext_dir)
+        if ( (bchi == BCType::direction_dependent) && is_velocity &&
+                (n == YVEL) && (s(i,domhi+1,k,n) <= 0.0) )
+        {
+            hi = s(i, domhi+1, k, n);
+            lo = s(i, domhi+1, k, n);
+        }
+        else if ( (bchi == BCType::ext_dir) ||
+             (bchi == BCType::direction_dependent && velm <= 0.0) )
         {
             hi = s(i,domhi+1,k,n);
             if( n == YVEL && is_velocity ) lo = hi ;
         }
-        else if (bchi == BCType::foextrap || bchi == BCType::hoextrap)
+        else if ( bchi == BCType::foextrap || bchi == BCType::hoextrap ||
+                 (bchi== BCType::direction_dependent && velm > 0.0) )
         {
             hi = lo;
         }
@@ -186,8 +218,8 @@ void SetYBCs (const int i, const int j, const int k, const int n,
         else if (bchi == BCType::reflect_odd)
         {
             if ( j==domhi+1 ) {
-                hi = 0.;
-                lo = 0.;
+                hi = Real(0.);
+                lo = Real(0.);
             } else {
                 Abort("EBGodunovBC::SetYBCs not yet fully implemented for reflect_odd BC. See comments in EBGodunovBC::SetXBCs.");
             //     hi = -lo_arr(i,2*(domhi+1)-j-1,k,n);
@@ -201,13 +233,13 @@ void SetYBCs (const int i, const int j, const int k, const int n,
 
 #if (AMREX_SPACEDIM==3)
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void SetZBCs(const int i, const int j, const int k, const int n,
-             const amrex::Array4<const amrex::Real> &s,
-             amrex::Real &lo,
-             amrex::Real &hi,
-             const int bclo, const int bchi,
-             const int domlo, const int domhi,
-             const bool is_velocity)
+void SetZBCs (const int i, const int j, const int k, const int n,
+              const amrex::Array4<const amrex::Real> &s,
+              amrex::Real &lo, amrex::Real &hi,
+              amrex::Real velm, amrex::Real velp,
+              const int bclo, const int bchi,
+              const int domlo, const int domhi,
+              const bool is_velocity)
 {
     using namespace amrex;
 
@@ -215,12 +247,20 @@ void SetZBCs(const int i, const int j, const int k, const int n,
     // Low Z
     if (k <= domlo)
     {
-        if (bclo==BCType::ext_dir)
+        if ( (bclo == BCType::direction_dependent) && is_velocity &&
+                (n == ZVEL) && (s(i,j,domlo-1,n) >= 0.0) )
+        {
+            lo = s(i, j, domlo-1, n);
+            hi = s(i, j, domlo-1, n);
+        }
+        else if ( bclo == BCType::ext_dir ||
+                 (bclo == BCType::direction_dependent && velp >= 0.0) )
         {
             lo =s(i,j,domlo-1,n);
             if ( n == ZVEL && is_velocity ) hi = lo;
         }
-        else if (bclo == BCType::foextrap || bclo == BCType::hoextrap)
+        else if ( bclo == BCType::foextrap || bclo == BCType::hoextrap ||
+                 (bclo == BCType::direction_dependent && velp < 0.0) )
         {
             lo = hi;
         }
@@ -234,8 +274,8 @@ void SetZBCs(const int i, const int j, const int k, const int n,
         {
             Abort("EBGodunovBC::SetZBCs not yet implemented for reflect_odd BC. See comments in EBGodunovBC::SetXBCs.");
             // if ( k==domlo ) {
-            //     hi = 0.;
-            //     lo = 0.;
+            //     hi = Real(0.);
+            //     lo = Real(0.);
             // } else {
             //     lo = -hi_arr(i,j,2*domlo-k  ,n);
             //     hi = -lo_arr(i,j,2*domlo-k-1,n);
@@ -245,12 +285,20 @@ void SetZBCs(const int i, const int j, const int k, const int n,
     // High Z
     else if (k > domhi)
     {
-        if (bchi==BCType::ext_dir)
+        if ( (bchi == BCType::direction_dependent) && is_velocity &&
+                (n == ZVEL) && (s(i,j,domhi+1,n) <= 0.0) )
+        {
+            hi = s(i, j, domhi+1, n);
+            lo = s(i, j, domhi+1, n);
+        }
+        else if ( (bchi == BCType::ext_dir) ||
+             (bchi == BCType::direction_dependent && velm <= 0.0) )
         {
             hi = s(i,j,domhi+1,n);
             if ( n == ZVEL && is_velocity ) lo = hi ;
         }
-        else if (bchi == BCType::foextrap || bchi == BCType::hoextrap)
+        else if ( bchi == BCType::foextrap || bchi == BCType::hoextrap ||
+                 (bchi== BCType::direction_dependent && velm > 0.0) )
         {
             hi = lo;
         }
@@ -264,8 +312,8 @@ void SetZBCs(const int i, const int j, const int k, const int n,
         {
             Abort("EBGodunovBC::SetZBCs not yet implemented for reflect_odd BC. See comments in EBGodunovBC::SetXBCs.");
             // if ( k==domhi+1 ) {
-            //     hi = 0.;
-            //     lo = 0.;
+            //     hi = Real(0.);
+            //     lo = Real(0.);
             // } else {
             //     hi = -lo_arr(i,j,2*(domhi+1)-k-1,n);
             //     lo = -hi_arr(i,j,2*(domhi+1)-k  ,n);
diff --git a/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp b/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp
index 852fb97da..c813017b3 100644
--- a/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp
+++ b/EBGodunov/hydro_ebgodunov_edge_state_2D.cpp
@@ -148,6 +148,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, v_mac(i,j,k), v_mac(i,j,k),
                                  bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
+            Real vad = v_mac(i,j,k);
             Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
             Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
             yzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_yzhi + l_yzlo);
@@ -267,6 +268,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
             HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, u_mac(i,j,k), u_mac(i,j,k),
                                  bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
+            Real uad = u_mac(i,j,k);
             Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
             Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
             xzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo);
diff --git a/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp b/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
index 78c7c2dec..226f2b681 100644
--- a/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
+++ b/EBGodunov/hydro_ebgodunov_edge_state_3D.cpp
@@ -354,6 +354,7 @@ EBGodunov::ComputeEdgeState ( Box const& bx, int ncomp,
         HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, u_mac(i,j,k), u_mac(i,j,k),
                                  bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
+        Real uad = u_mac(i,j,k);
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
         xzlo(i,j,k,n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (l_xzhi + l_xzlo);
diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp
index 64c8668c8..ffbd52207 100644
--- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp
+++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces.cpp
@@ -277,7 +277,8 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx,
             Real hi = Imx(i  ,j,k,n);
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            EBGodunovBC::SetXBCs(i, j, k, n, vel, lo, hi, bc.lo(0), bc.hi(0), dlo.x, dhi.x, true);
+            EBGodunovBC::SetXBCs(i, j, k, n, vel, lo, hi, lo, hi,
+                                 bc.lo(0), bc.hi(0), dlo.x, dhi.x, true);
 
             Real st = ( (lo+hi) >= 0.) ? lo : hi;
             bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) );
@@ -297,7 +298,8 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx,
             Real hi = Imy(i,j  ,k,n);
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            EBGodunovBC::SetYBCs(i, j, k, n, vel, lo, hi, bc.lo(1), bc.hi(1), dlo.y, dhi.y, true);
+            EBGodunovBC::SetYBCs(i, j, k, n, vel, lo, hi, lo, hi,
+                                 bc.lo(1), bc.hi(1), dlo.y, dhi.y, true);
 
             Real st = ( (lo+hi) >= 0.) ? lo : hi;
             bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) );
@@ -318,7 +320,8 @@ EBGodunov::ComputeAdvectiveVel ( AMREX_D_DECL(Box const& xbx,
             Real hi = Imz(i,j,k  ,n);
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
-            EBGodunovBC::SetZBCs(i, j, k, n, vel, lo, hi, bc.lo(2), bc.hi(2), dlo.z, dhi.z, true);
+            EBGodunovBC::SetZBCs(i, j, k, n, vel, lo, hi, lo, hi,
+                                 bc.lo(2), bc.hi(2), dlo.z, dhi.z, true);
 
             Real st = ( (lo+hi) >= 0.) ? lo : hi;
             bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) );
diff --git a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
index 928d66b50..f7ecdb56f 100644
--- a/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
+++ b/EBGodunov/hydro_ebgodunov_extrap_vel_to_faces_3D.cpp
@@ -75,14 +75,15 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
             Real lo = Ipx(i-1,j,k,n);
             Real hi = Imx(i  ,j,k,n);
 
-            Real uad = u_ad(i,j,k);
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            EBGodunovBC::SetXBCs(i, j, k, n, q, lo, hi, bc.lo(0), bc.hi(0), dlo.x, dhi.x, true);
+            EBGodunovBC::SetXBCs(i, j, k, n, q, lo, hi, u_ad(i,j,k), u_ad(i,j,k),
+                                 bc.lo(0), bc.hi(0), dlo.x, dhi.x, true);
 
             xlo(i,j,k,n) = lo;
             xhi(i,j,k,n) = hi;
 
+            Real uad = u_ad(i,j,k);
             Real st = (uad >= 0.) ? lo : hi;
             Real fu = (amrex::Math::abs(uad) < small_vel) ? Real(0.0) : Real(1.0);
             Imx(i, j, k, n) = fu*st + (Real(1.0) - fu) * Real(0.5) * (hi + lo); // store xedge
@@ -92,14 +93,15 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
             Real lo = Ipy(i,j-1,k,n);
             Real hi = Imy(i,j  ,k,n);
 
-            Real vad = v_ad(i,j,k);
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            EBGodunovBC::SetYBCs(i, j, k, n, q, lo, hi, bc.lo(1), bc.hi(1), dlo.y, dhi.y, true);
+            EBGodunovBC::SetYBCs(i, j, k, n, q, lo, hi, v_ad(i,j,k), v_ad(i,j,k),
+                                 bc.lo(1), bc.hi(1), dlo.y, dhi.y, true);
 
             ylo(i,j,k,n) = lo;
             yhi(i,j,k,n) = hi;
 
+            Real vad = v_ad(i,j,k);
             Real st = (vad >= 0.) ? lo : hi;
             Real fu = (amrex::Math::abs(vad) < small_vel) ? Real(0.0) : Real(1.0);
             Imy(i, j, k, n) = fu*st + (Real(1.0) - fu)*Real(0.5)*(hi + lo); // store yedge
@@ -111,7 +113,8 @@ EBGodunov::ExtrapVelToFacesOnBox ( Box const& bx, int ncomp,
 
             const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr);
 
-            EBGodunovBC::SetZBCs(i, j, k, n, q, lo, hi, bc.lo(2), bc.hi(2), dlo.z, dhi.z, true);
+            EBGodunovBC::SetZBCs(i, j, k, n, q, lo, hi, w_ad(i,j,k), w_ad(i,j,k),
+                                 bc.lo(2), bc.hi(2), dlo.z, dhi.z, true);
 
             zlo(i,j,k,n) = lo;
             zhi(i,j,k,n) = hi;
diff --git a/Godunov/hydro_godunov_edge_state_2D.cpp b/Godunov/hydro_godunov_edge_state_2D.cpp
index 088eec481..033483475 100644
--- a/Godunov/hydro_godunov_edge_state_2D.cpp
+++ b/Godunov/hydro_godunov_edge_state_2D.cpp
@@ -270,6 +270,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
         HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, umac(i,j,k), umac(i,j,k),
                              bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
+        Real uad = umac(i,j,k);
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
         xzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo);
diff --git a/Godunov/hydro_godunov_edge_state_3D.cpp b/Godunov/hydro_godunov_edge_state_3D.cpp
index d89ffd672..6648d8af4 100644
--- a/Godunov/hydro_godunov_edge_state_3D.cpp
+++ b/Godunov/hydro_godunov_edge_state_3D.cpp
@@ -242,6 +242,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
         HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zylo, l_zyhi, wmac(i,j,k), wmac(i,j,k),
                              bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
+        Real wad = wmac(i,j,k);
         Real st = (wad >= 0.) ? l_zylo : l_zyhi;
         Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0;
         zylo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_zyhi + l_zylo);
@@ -259,6 +260,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
         HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yzlo, l_yzhi, vmac(i,j,k), vmac(i,j,k),
                              bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
+        Real vad = vmac(i,j,k);
         Real st = (vad >= 0.) ? l_yzlo : l_yzhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
         yzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_yzhi + l_yzlo);
@@ -346,6 +348,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
         HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xzlo, l_xzhi, umac(i,j,k), umac(i,j,k),
                              bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
+        Real uad = umac(i,j,k);
         Real st = (uad >= 0.) ? l_xzlo : l_xzhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
         xzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo);
@@ -363,6 +366,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
         HydroBC::SetZEdgeBCs(i, j, k, n, q, l_zxlo, l_zxhi, wmac(i,j,k), wmac(i,j,k),
                              bc.lo(2), dlo.z, bc.hi(2), dhi.z, is_velocity);
 
+        Real wad = wmac(i,j,k);
         Real st = (wad >= 0.) ? l_zxlo : l_zxhi;
         Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0;
         zxlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_zxhi + l_zxlo);
@@ -448,6 +452,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
         HydroBC::SetXEdgeBCs(i, j, k, n, q, l_xylo, l_xyhi, umac(i,j,k), umac(i,j,k),
                              bc.lo(0), dlo.x, bc.hi(0), dhi.x, is_velocity);
 
+        Real uad = umac(i,j,k);
         Real st = (uad >= 0.) ? l_xylo : l_xyhi;
         Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0;
         xylo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_xyhi + l_xylo);
@@ -465,6 +470,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp,
         HydroBC::SetYEdgeBCs(i, j, k, n, q, l_yxlo, l_yxhi, vmac(i,j,k), vmac(i,j,k),
                              bc.lo(1), dlo.y, bc.hi(1), dhi.y, is_velocity);
 
+        Real vad = vmac(i,j,k);
         Real st = (vad >= 0.) ? l_yxlo : l_yxhi;
         Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0;
         yxlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_yxhi + l_yxlo);
diff --git a/Godunov/hydro_godunov_ppm.H b/Godunov/hydro_godunov_ppm.H
index 5af7e10d1..739270189 100644
--- a/Godunov/hydro_godunov_ppm.H
+++ b/Godunov/hydro_godunov_ppm.H
@@ -289,7 +289,7 @@ void SetXBCs ( const int i, const int j, const int k, const int n,
     using namespace amrex;
 
     if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) ||
-         (bclo == BCType::direction_dependent && velm >= 0.0) )
+         (bclo == BCType::direction_dependent && velp >= 0.0) )
     {
         if (i == domlo)
         {
@@ -330,7 +330,7 @@ void SetXBCs ( const int i, const int j, const int k, const int n,
     }
 
     if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) ||
-         (bchi == BCType::direction_dependent && velp <= 0.0) )
+         (bchi == BCType::direction_dependent && velm <= 0.0) )
     {
         if (i == domhi)
         {
@@ -743,7 +743,7 @@ void PredictStateOnXFace ( const int i, const int j, const int k, const int n,
     amrex::Real sm, sp;
     amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
 
-    SetXBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i+1,j,k),
+    SetXBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k),
                      bc.lo(0), bc.hi(0), domlo, domhi);
 
     amrex::Real s6 = 6.0*s0 - 3.0*(sm + sp);
@@ -790,7 +790,7 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n,
     amrex::Real sm, sp;
     amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
 
-    SetYBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j+1,k),
+    SetYBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k),
                      bc.lo(1), bc.hi(1), domlo, domhi);
 
     amrex::Real s6 = 6.0*s0- 3.0*(sm + sp);
@@ -840,7 +840,7 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n,
     amrex::Real sm, sp;
     amrex::Tie(sm,sp) = Limiter::sm_sp(s0,sedge1,sedge2);
 
-    SetZBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k+1),
+    SetZBCs<Limiter>(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k),
                      bc.lo(2), bc.hi(2), domlo, domhi);
 
     amrex::Real s6 = 6.0*s0- 3.0*(sm + sp);
diff --git a/Utils/hydro_bcs_K.H b/Utils/hydro_bcs_K.H
index b0bdd8160..6fd578e79 100644
--- a/Utils/hydro_bcs_K.H
+++ b/Utils/hydro_bcs_K.H
@@ -82,22 +82,27 @@ void SetXEdgeBCs( int i, int j, int k, int n,
                   int bchi, int domhi,
                   bool is_velocity )
 {
-
     using namespace amrex;
 
     // This function is for setting BCs ON domain faces only.
+    // NOTE: this is called in two ways:
+    //  1) in the extrapolation of velocity components only --
+    //     in this case velm and velp may be predicted
+    //  2) for all components in the construction of edge states
+    //     in this case velm and velp are identical and equal to
+    //     the MAC velocity on the current edge
     AMREX_ASSERT( domlo <= i && i <= domhi+1 );
 
     if (i == domlo)
     {
         if ( (bclo == BCType::direction_dependent) && is_velocity &&
-                  (n == XVEL) && (s(domlo-1,j,k,n) >= 0.0) )
+                (n == XVEL) && (s(domlo-1,j,k,n) >= 0.0) )
         {
             lo = s(domlo-1, j, k, n);
-            if ( n==XVEL && is_velocity ) hi=lo;
+            hi = s(domlo-1, j, k, n);
         }
         else if ( (bclo == BCType::ext_dir) ||
-             (bclo == BCType::direction_dependent && velm >= 0.0) ) {
+                  (bclo == BCType::direction_dependent && velp >= 0.0) )
         {
             lo = s(domlo-1, j, k, n);
             // For turbulent inflow, there are times when the inflow face
@@ -114,13 +119,19 @@ void SetXEdgeBCs( int i, int j, int k, int n,
         }
         else if (bclo == BCType::reflect_odd)
         {
-            hi = 0.;
-            lo = hi;
+            hi = Real(0.);
+            lo = Real(0.);
         }
     }
     else if (i == domhi+1)
     {
-        if ( (bchi == BCType::ext_dir) ||
+        if ( (bchi == BCType::direction_dependent) && is_velocity &&
+                (n == XVEL) && (s(domhi+1,j,k,n) <= 0.0) )
+        {
+            hi = s(domhi+1, j, k, n);
+            lo = s(domhi+1, j, k, n);
+        }
+        else if ( (bchi == BCType::ext_dir) ||
              (bchi == BCType::direction_dependent && velm <= 0.0) )
         {
             hi = s(domhi+1, j, k, n);
@@ -134,8 +145,8 @@ void SetXEdgeBCs( int i, int j, int k, int n,
         }
         else if (bchi == BCType::reflect_odd)
         {
-            lo = 0.;
-            hi = lo;
+            lo = Real(0.);
+            hi = Real(0.);
         }
     }
     else
@@ -167,8 +178,14 @@ void SetYEdgeBCs (int i, int j, int k, int n,
 
     if (j == domlo)
     {
-        if ( (bclo == BCType::ext_dir) ||
-             (bclo == BCType::direction_dependent && velp >= 0.0) )
+        if ( (bclo == BCType::direction_dependent) && is_velocity &&
+                (n == YVEL) && (s(i,domlo-1,k,n) >= 0.0) )
+        {
+            lo = s(i, domlo-1, k, n);
+            hi = s(i, domlo-1, k, n);
+        }
+        else if ( (bclo == BCType::ext_dir) ||
+                  (bclo == BCType::direction_dependent && velp >= 0.0) )
         {
             lo = s(i, domlo-1, k, n);
             if ( n==YVEL && is_velocity ) hi=lo;
@@ -181,14 +198,20 @@ void SetYEdgeBCs (int i, int j, int k, int n,
         }
         else if (bclo == BCType::reflect_odd)
         {
-            hi = 0.;
-            lo = hi;
+            hi = Real(0.);
+            lo = Real(0.);
         }
     }
     else  if (j == domhi+1)
     {
-        if ( (bchi == BCType::ext_dir) ||
-             (bchi == BCType::direction_dependent && velm <= 0.0) )
+        if ( (bchi == BCType::direction_dependent) && is_velocity &&
+                (n == YVEL) && (s(i,domhi+1,k,n) <= 0.0) )
+        {
+            hi = s(i, domhi+1, k, n);
+            lo = s(i, domhi+1, k, n);
+        }
+        else if ( (bchi == BCType::ext_dir) ||
+                  (bchi == BCType::direction_dependent && velm <= 0.0) )
         {
             hi = s(i, domhi+1, k, n);
             if ( n==YVEL && is_velocity ) lo=hi;
@@ -199,10 +222,10 @@ void SetYEdgeBCs (int i, int j, int k, int n,
         {
             hi = lo;
         }
-        else if(bchi == BCType::reflect_odd)
+        else if (bchi == BCType::reflect_odd)
         {
-            lo = 0.;
-            hi = lo;
+            lo = Real(0.);
+            hi = Real(0.);
         }
     }
     else
@@ -239,8 +262,14 @@ void SetZEdgeBCs (int i, int j, int k, int n,
 
     if (k == domlo)
     {
-        if ( (bclo == BCType::ext_dir) ||
-             (bclo == BCType::direction_dependent && velp >= 0.0) )
+        if ( (bclo == BCType::direction_dependent) && is_velocity &&
+                (n == ZVEL) && (s(i,j,domlo-1,n) >= 0.0) )
+        {
+            lo = s(i, j, domlo-1, n);
+            hi = s(i, j, domlo-1, n);
+        }
+        else if ( (bclo == BCType::ext_dir) ||
+                  (bclo == BCType::direction_dependent && velp >= 0.0) )
         {
             lo = s(i, j, domlo-1, n);
             if ( n==ZVEL && is_velocity ) hi=lo;
@@ -253,16 +282,22 @@ void SetZEdgeBCs (int i, int j, int k, int n,
         }
         else if(bclo == BCType::reflect_odd)
         {
-            hi = 0.;
-            lo = hi;
+            hi = Real(0.);
+            lo = Real(0.);
         }
     }
     else if (k == domhi+1)
     {
-        if ( (bchi == BCType::ext_dir) ||
-             (bchi == BCType::direction_dependent && velm <= 0.0) )
+        if ( (bchi == BCType::direction_dependent) && is_velocity &&
+                (n == ZVEL) && (s(i,j,domhi+1,n) <= 0.0) )
         {
-            hi = s(i,j,domhi+1, n);
+            hi = s(i, j, domhi+1, n);
+            lo = s(i, j, domhi+1, n);
+        }
+        else if ( (bchi == BCType::ext_dir) ||
+                  (bchi == BCType::direction_dependent && velm <= 0.0) )
+        {
+            hi = s(i, j, domhi+1, n);
             if ( n==ZVEL && is_velocity ) lo=hi;
         }
         else if ( bchi == BCType::foextrap || bchi == BCType::hoextrap ||
@@ -273,8 +308,8 @@ void SetZEdgeBCs (int i, int j, int k, int n,
         }
         else if(bchi == BCType::reflect_odd)
         {
-            lo = 0.;
-            hi = lo;
+            lo = Real(0.);
+            hi = Real(0.);
         }
     }
     else