Radiation correction is not done in parallel (even when using Dask)
## Issue

If the user chooses to apply a radiation correction for the surface forcing while using Dask via
```
SurfaceForcing(
    grid=grid,
    ...,
    correct_radiation=True,
    use_dask=True
)
```
this radiation correction is not done in parallel. The reason for this is the internal use of the function 
https://github.com/CWorthy-ocean/roms-tools/blob/67df5fc400b1daeb903ba9b2efe11899d8b55dfa/roms_tools/setup/utils.py#L128-L150

which interpolates a pre-computed climatology for ERA5 shortwave downward radiation (i.e., 12 time stamps) to the hourly ERA5 data of interest. Specifically, this line is the issue:

https://github.com/CWorthy-ocean/roms-tools/blob/67df5fc400b1daeb903ba9b2efe11899d8b55dfa/roms_tools/setup/utils.py#L184-L187

Both `data_array_concat` and `day_of_year` are chunked in time with chunk size 1, where `data_array_concat.time` has a length of 12 + 2, and `day_of_year.time` a length of however many ERA5 times there are: for a month of ERA5 data, that is 744 hourly time slices, i.e., 744 time chunks. However, the `xarray.interp` function cannot deal with chunked data along the core dimension and will return interpolated data with a single time chunk (containing 744 time slices). This means that this operation is not done in parallel.

Note that substituting this by something like
```
        data_array_interpolated = xr.apply_ufunc(
            np.interp, 
            day_of_year, 
            time_concat, 
            data_array_concat, 
            input_core_dims=[[time_dim_name], [time_dim_name], [time_dim_name]],
            output_core_dims=[[time_dim_name]],
            vectorize=True,
            dask="allowed",
            exclude_dims=set((time_dim_name,)),
        )
```
would not help because we are trying to chunk along a core dimension, which `xr.apply_ufunc` is not equipped for (I think).

## Implications for validation step and memory footprint

In the validation step during the surface forcing initialization, a single time slice of the corrected shortwave radiation is checked for NaNs:

https://github.com/CWorthy-ocean/roms-tools/blob/67df5fc400b1daeb903ba9b2efe11899d8b55dfa/roms_tools/setup/surface_forcing.py#L366

However, due to the issue above, this will imply a huge memory footprint because the shortwave radiation **correction** is computed for all time slices (e.g., 744 time stamps), not just one.

This can be seen in the following example, where we use Grid A from https://github.com/CWorthy-ocean/roms-tools/issues/201 and one month of hourly ERA5 data:

```
SurfaceForcing(
    grid=grid,
    ...,
    correct_radiation=True,
    use_dask=True,
    bypass_validation=True
)
```
has the following performance:
```
peak memory: 562.75 MiB, increment: 25.68 MiB
CPU times: user 54.8 s, sys: 1.26 s, total: 56.1 s
Wall time: 30.1 s
```
The memory footprint is minimal because the validation step was skipped.

However,

```
SurfaceForcing(
    grid=grid,
    ...,
    correct_radiation=True,
    use_dask=True,
    bypass_validation=False
)
```
has the following performance and big memory footprint:
```
peak memory: 32898.70 MiB, increment: 32335.71 MiB
CPU times: user 2min 24s, sys: 5.01 s, total: 2min 29s
Wall time: 43.8 s
```