More Flexible Flux implementation

So after implementing an upwinding flux by modifying the penalty term in the Rusanov flux [1], I started thinking we could reduce the number of fluxes by moving the point of multiple dispatch to the penalty calculation.

Most of the fluxes we currently have/want are basically central flux + penalty, so we could create one type of flux for these and then dispatch based on the penalty type. For central, no penalty; for Lax-Friedrichs/Rusanov, the appropriate wavespeed times the jump in state variable; and for upwind, the appropriate choice of plus or minus side flux. This way we don't have the same compute the central flux kernel replicated in multiple places.

For something that isn't central flux + penalty, we could just compute the custom flux on the fly. While this could be implemented through the central flux + penalty model, you would have to go through the headache of zeroing out the central flux first, so it's simpler just to start from zero.

@lcw @jkozdon @sandreza @fxgiraldo what are your thoughts?

[1] https://github.com/climate-machine/CLIMA/commit/61ed8472c5a56c0aef9031e278d58ef61b860236