Avoiding negative flows where they should not be

[SouthEndMusic]

(Somewhat of a duplicate of #1838)

After the refactor #1819 we found that flows that should be non-negative are sometimes negative for a brief period of time (or maybe it already happened before, but we just now discovered it by looking at tracers?). I discussed this on Slack:

Several solutions were proposed above depending on where the problem originates:

• Add a stage_limiter to the nonlinear solver (@visr had some comments about this, it would be nice if you add them below)
• Switch to QR factorization in the linear solver

If the above options don't work we could also implement a callback for enforcing non-negative flows. Some notes about this:

• This shouldn't affect the water balance error, as it just nudges the physics in the right direction
• This involves changing the state in the callback, which might lower the performance of multi-step methods. Also, our experiments (with @visr) with using u_modified! weren't very promising: