It is possible that swapping the boundary conditions will lead to more stable results: The flow should be specified upstream and the pressure downstream.
A weno discretization of the space derivatives could also be helpful.
And
0 ~ ∂x(p(x, t)) + b * q_m(x, t) * abs(q_m(x, t)) / p(x, t),
instead of
∂x(p(x, t)) ~ - b * q_m(x, t) * abs(q_m(x, t)) / p(x, t),
seems preferable.