Fixes #1919.

Inspired by ongoing discussions with @Huite and @visr on smoothness (see also #1918, #1920) I started looking into making the type of interpolation configurable in the several places in the core where interpolation is used. For this PR I am focussing on $Q(h)$ relationships.

We require of the $Q(h)$ interpolation that the flow is $0$ below and at the lowest supplied level. That's trivial to achieve for linear interpolation, but not for non-linear (generally piecewise polynomial of degree >1) interpolation types. This requires a special approach for each method.

The question is which methods to support (and even whether this should be configurable for $Q(h)$ relationships at all). I now experiment with PCHIP Interpolation because of its smoothness and non-overshooting properties.

Another point of attention is extrapolation. For linear interpolation extrapolation makes sense, but for higher order interpolation methods extrapolation gets out of hand quick. We could make our own wrapper of interpolation types from DataInterpolation, but I also made an issue there to pick it up internally:

SciML/DataInterpolations.jl#355