Models sometimes return small-ish negative values

[wknoben]

Numerical solutions don't always return physically realistic values.

Suggestion:

• Remove constraints from individual flux functions
• Build in constraints as part of the solver, where we (1) check if a new solution dumps a store into the negative; (2) if so, reduce all relevant fluxes by a scaled amount.

Big change, so probably should be saved for a large development push.