Flux variable handling

Blocked by:

• #30 (closed)
• #25 (closed)

!40 (merged) outlines a couple of problems we need to solve in order to handle flux variables well. This issue is a place to discuss how we solve these problems.

My proposal:

• when we have flux variables, we report both cumulative flux and flux
• we always report spline assumptions when reporting cumulative flux - flux pairs so that downstream tooling can sensibly/consistently/mass- or energy-conservingly handle them
  • these two pieces of information mean we have enough information to infer cumulative flux at the end of the last time step, which means we don't have to carry around an extra time step in our cumulative flux timeseries (which would be a pain to handle as it would mean the flux variables always have a nan at the end we'd have to work around everywhere)
• we always report time bounds when reporting any outputs (including state variables)
• use xarray to do this handling. Given the discussion about xarray's accessors and type hints, I would be tempted to just create an object that holds an xarray object and does stuff with it, rather than going the accessor route (I think my proposed pattern also allows us to make the data validation more obvious, but that's just a feeling).
  • this xarray object would be where I would implement the left-hand integral (needed for data that should use a constant spline). xarray already has cumulative_integral, but it uses the trapezium rule, so it is a good base to copy/build from but isn't actually usable (I can then obviously make a PR back into xarray to add the ability to do this integration).