Make use of Threads.@Spawn and basesize to have more control on the amount of work per thread. Added this for vertical and kinematic wave routing components.

and removed eachindex when not required

Use only catchment scale stream order to partition river and land domain for parallel running of kinematic wave solution. Different minimum stream order values can be supplied for the river and land domain.

SurfaceFlow is used for both kinematic wave river and overland flow. This is now split into structs SurfaceFlowRiver and SurfaceFlowLand. It simplifies the orginal update function (split into two, logic easier to follow), and dispatching on these types is easier (e.g. for BMI to extract relevant properties as grid type).

use threaded_foreach instead of @threads

For threads <= 8 spawn tasks is used (outperforms Polyester threads in this range), and for threads > 8 Polyester threads is used (spawn tasks slows down in this range).

Reduce run time by applying loopvectorization only for edges above flow threshold and by pre-allocating more variables (river and floodplain).

and stable_timestep calculation is now faster with @tturbo.

and rename variable index of struct FloodPlain

- update model parameters docs because of split of struct SurfaceFlow in river and overland flow structs.
- remove field to_river from struct SurfaceFlowRiver (not used as part of river flow).

For Julia version < 1.8 Base.ifelse is not working with Loopvectorization (`LoopVectorization.check_args` on your inputs failed).