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