I’ve been trying to generate forcing data using 1° monthly ACCESS-OM2 output, whereas the examples I’ve seen use 0.1° daily data. I’ve run into two issues so far:
I use the monthly data to create access_om2_input, and from this the boundary forcing files (north/south/east/west_unprocessed.nc). When these files are passed to expt.setup_ocean_state_boundaries (the updated version of expt.rectangular_boundaries in the example recipes), the function assumes the input is daily data. The boundary forcing files it generates don’t cover the full time period I expected them to, (two years of monthly data is misinterpreted as just 24 days). I don’t know if the forcing data has to be daily, or if monthly would be fine as long as it was dated correctly in the files generated!
I could get around this by adding this code in before the north/south/east/west_unprocessed.nc files are saved.
daily_time = pd.date_range(start=date_range[0], end= date_range[1], freq='D')
def interp_to_daily(ds):
return ds.interp(time=daily_time)
## Cut out East & West boundary conditions and save to netCDF
print('East & West boundary condition...')
east = rmom6.longitude_slicer(access_om2_input,
[longitude_extent[1] - buffer, longitude_extent[1] + buffer],
["xu_ocean", "xt_ocean"])
interp_to_daily(east).to_netcdf(tmp_dir + "/east_unprocessed.nc")
The ACCESS-OM2 data I used has 50 vertical levels (nz=50), but in rmom6.experiment I had set number_vertical_layers = 75. Expt.setup_initial_condition() handled the vertical regridding correctly, and the resulting initial condition files had nz=75 as expected. Expt.setup_ocean_state_boundaries produced boundary forcing files that still only had nz=50.
Again I could get around this by adding something like this into the notebook before calling Expt.setup_ocean_state_boundaries
rg_output_dir = tmp_dir+"/regrid/"
os.makedirs(rg_output_dir, exist_ok=True)
zl_interp = expt.vgrid.zl.rename({'zl': 'st_ocean'})
directions = ["north", "south", "east", "west"]
for direction in directions:
input_file = os.path.join(tmp_dir, f"{direction}_unprocessed.nc")
output_file = os.path.join(rg_output_dir, f"{direction}_unprocessed.nc")
print(f"Interpolating {input_file}")
ds = xr.open_dataset(input_file)
ds_interp = ds.interp({'st_ocean': zl_interp})
ds_interp.to_netcdf(output_file)
print(f"Saved to {output_file}")
and using rg_output_dir as the inputput for exp.setup_ocean_state_boundaries(raw_boundaries_path=Path(rg_output_dir))
Thanks for posting @Lizzie. Are you hoping to get help from ACCESS-NRI with this? Adding the help tag to your post will mean it gets triaged by our support team.
I can’t really help except to say that I noticed the same thing about the boundary files still having the native vertical resolution of the forcing dataset. I think the interpolation might be handled by MOM6 although I’m not totally sure.
For your first question - if you don’t have any luck here you could try an issue on the regional-mom6 repo. I’m interested in the answer since I have just discovered that very few of the ACCESS-OM2 models save daily u, v, temp, salt, and ssh. I don’t know if i just suck at using intake, but as far as I can tell only 01deg_jra55v13_ryf9091 and 01deg_jra55v150_iaf_cycle1 have all the variables that regional-mom6 wants - and maybe they only have daily SSH in January?
Thanks for bringing this up @Lizzie! I would really appreciate it if you could raise it as an issue in the GitHub repository. All the testing has been done on Glorys data which is daily so would be very easy for this bug to be missed.
I think this issue warrants a larger discussion. For regional modelling, we tend to use daily or higher frequency forcing so I think the lack of ACCESS-OM2 models with daily outputs will become an issue for a lot of regional mom6 users. Is there any energy for a discussion about how we might set up a model run with higher frequency outputs?
Yeah the OBCs don’t need to be regrided vertically - MOM6 does this online thankfully so you don’t have to have extra big boundary files to match your model’s (presumably finer) vertical grid.
Regarding the time dimension: the boundary class does have the option to specify the time units of the boundaries, however this hasn’t been implemented in the wrapper method setup_ocean_state_boundaries which we use in the demo.
For context, users used to have to call the segment class for each segment (and still do if they want to do more complex things with their boundaries). But we decided to wrap this up in one function to keep the demo notebook streamlined and easier for beginners. But clearly we forgot to include some edge cases! An easy fix is just to add a key word argument to the ocean state boundaries function, which passes the time units onto the boundary segment class on initialisation.
@Lizzie would you have time to do this as a quick PR?
To fix the data you’ve already made, you just need to modify the metadata to correct the time units. At the end of the day that’s all that rmom6 does: writes the appropriate time encoding into the netcdf file before saving it.