Scores Aggregation: Update scoreswrapper functions to support multiple forecasts in a cube

### Something else
Currently, our implementation of `scoreswrapper.py` operators assume a model, model, or model, obs cube that contain one forecast. We should adapt them so they compute the scores over each forecast separately, and then combine them at the end. This might not be necessary if we can potentially preserve the `forecast_reference_time` as a coordinate when we pass to the scores operator, but I've not tried that. It would be good to check this gives the same value, if it is possible.

I have some code that loops over `forecast_reference_time` if it is a dimension, on branch `2409_rmse_agg_hinton` https://github.com/MetOffice/CSET/pull/2419/changes#diff-42de57b369cdd291381c081213b5a4970b790f359ff365d878a9d5a35ad69ebe