diff --git a/src/lbmpy/boundaries/boundaryconditions.py b/src/lbmpy/boundaries/boundaryconditions.py
index 457aa513d792ce273ada95c13ff41c00d28c96f0..c7acebdf6177df57d3490d9c9eeef584b44d084e 100644
--- a/src/lbmpy/boundaries/boundaryconditions.py
+++ b/src/lbmpy/boundaries/boundaryconditions.py
@@ -926,15 +926,18 @@ class iMEMBounceBack(LbBoundary):
             u_tau_target=self.target_friction_velocity)
         result += wall_law_assignments
 
-        # for i in self.tangential_axis:
-        #     result.append(Assignment(wall_stress[i], u_for_tau_wall[i] / u_m * tau_w ))
+        # Only required if not provided by bounce-back object already
+        cqe = cqc.equilibrium_input_equations_from_pdfs(pdf_symbols, False)
+        for subexpressions in cqe.all_assignments:
+                if subexpressions.lhs not in unique_lhs:
+                    unique_lhs.add( subexpressions.lhs )
+                    result.append(subexpressions)
 
-        # Could also be:
         u_norm = sp.Symbol("u_norm")
         u_mag = Assignment(u_norm, sp.sqrt(sum([self.mean_velocity[i]**2 for i in self.tangential_axis])))
         result.append(u_mag)
         for i in self.tangential_axis:
-            result.append(Assignment(wall_stress[i], cqc.velocity_symbols[i] / u_norm * tau_w ))
+            result.append(Assignment(wall_stress[i], - cqc.velocity_symbols[i] / u_norm * tau_w ))
 
     # Step 3: Compute F_f and F_wm eqs. 24 and 27
         # F_fluid
@@ -970,10 +973,18 @@ class iMEMBounceBack(LbBoundary):
 
         rho = cqc.density_symbol if cqc.compressible else 1
 
-        # rho = 1
+        #! BUG! Hard-coded. This should be q from index list (if using quadratic bb... else default to 0)
+        # needed in extra momentum term... see Eq. 34 of  Pasquali et al. Near-wall treamtent.
+        q = 0.5
+
+        velocity = [v_i for v_i in wall_velocity]
+        # shifted_vel_eqs = cqc.equilibrium_input_equations_from_init_values(velocity=velocity)
+        # shifted_vel_eqs = shifted_vel_eqs.new_without_subexpressions()
+        # velocity = [eq.rhs for eq in shifted_vel_eqs.new_filtered(cqc.velocity_symbols).main_assignments]
+
         velocity_expression = [ dx4_on_dt2 * sum( c_f_ijk[i][k] * 
-                                        (- 2 * wi[i] * rho * sum( wall_velocity[j] * c_f_ijk[i][j] for j in range(self.stencil.D)
-                                                              ) / c_s_sq) for i in pdfs_in_wall_dir
+                                        (- 2 * wi[i] * rho * sum( velocity[j] * c_f_ijk[i][j] for j in range(self.stencil.D)
+                                                              ) / c_s_sq / (1+q)) for i in pdfs_in_wall_dir
                                     ) for k in range(self.stencil.D) ]
 
         result += [ Assignment( wall_velocity[i], sp.solve(velocity_expression[i] - F_uw[i], wall_velocity[i])[0] ) for i in self.tangential_axis]
@@ -981,7 +992,7 @@ class iMEMBounceBack(LbBoundary):
         
         clipped_velocity = sp.symbols("clip_u_wx clip_u_wy clip_u_wz")
 
-        u_for_clipping = tuple( 0.5 * u_mean_i.get_shifted( normal_direction[0], normal_direction[1], normal_direction[2]
+        u_for_clipping = tuple( u_mean_i.get_shifted( normal_direction[0], normal_direction[1], normal_direction[2]
             ) for u_mean_i in self.inst_velocity)
 
         result += [ Assignment( clipped_velocity[i],  
@@ -1007,9 +1018,10 @@ class iMEMBounceBack(LbBoundary):
 
         # velocity = [v_i for v_i in wall_velocity]
         velocity = [v_i for v_i in clipped_velocity]
-        shifted_vel_eqs = cqc.equilibrium_input_equations_from_init_values(velocity=velocity)
-        shifted_vel_eqs = shifted_vel_eqs.new_without_subexpressions()
-        velocity = [eq.rhs for eq in shifted_vel_eqs.new_filtered(cqc.velocity_symbols).main_assignments]
+        # # # # velocity = [v_i for v_i in clipped_velocity]
+        # shifted_vel_eqs = cqc.equilibrium_input_equations_from_init_values(velocity=velocity)
+        # shifted_vel_eqs = shifted_vel_eqs.new_without_subexpressions()
+        # velocity = [eq.rhs for eq in shifted_vel_eqs.new_filtered(cqc.velocity_symbols).main_assignments]
 
         weight_info = LbmWeightInfo(lb_method, data_type=self.data_type)
         weight_of_direction = weight_info.weight_of_direction
@@ -1018,7 +1030,7 @@ class iMEMBounceBack(LbBoundary):
 
     # Step 6: Return Assignment
         # Eq. 20 in Asmuth2021: Note, eq. 20 should not have the 2x on LHS and should be - not +
-        result.append(Assignment(f_in(inv_dir[dir_symbol]),  f_ijk_bb - f_ijk_uw) )
+        result.append(Assignment(f_in(inv_dir[dir_symbol]),  f_ijk_bb - f_ijk_uw/(1+q)) )
 
         # for ass in result:
         #     print("")