Adding restoring forces DIC/Alk#619
Conversation
Codecov Report❌ Patch coverage is
@@ Coverage Diff @@
## main #619 +/- ##
==========================================
- Coverage 88.53% 88.43% -0.11%
==========================================
Files 28 28
Lines 7073 7141 +68
Branches 1242 1258 +16
==========================================
+ Hits 6262 6315 +53
- Misses 467 477 +10
- Partials 344 349 +5
Flags with carried forward coverage won't be shown. Click here to find out more.
🚀 New features to boost your workflow:
|
|
Check out this pull request on See visual diffs & provide feedback on Jupyter Notebooks. Powered by ReviewNB |
|
The PR was tested by making a restoring forces file for the GoM, and running ROMS for 1 month (Jan 2000). The restoring forces file was also compared to the data used to make it as a gut check. These are the initial conditions (from GLODAP) and the restoring forces (from SODA), and their differences: Based on these differences, we can calculate an estimate of the expected flux in the model ( The 'expected flux' (above) values (order 1e-4) match well with the output diagnostics of restoring force fluxes from the model (below): 
The flux values are large that what Pierre saw with the Atlantic (1e-6), but we attribute the GoM fluxes being large to the large difference between the IC and the applied restoring forces, and because the model was only run for 1 month so there wasn't much time for the model to spin up and adjust. Additionally, the file produced by ROMS-Tools was compared to the sub-sampled portion of the source file used to make the restoring forces, in order to make sure they are of similar magnitude and RT does what we expect it to do. Note, the below plots are from a different domain (the one used in the DIC calculated from ROMS-Tools: The values fall within a similar range that we'd expect. |
|
|
||
| return ds | ||
|
|
||
| def post_process(self) -> None: |
There was a problem hiding this comment.
I have not fully checked, but this smells not dask lazy, because the .where is called outside the dask framework and xarray will request the data. Maybe replace with something like this:
def post_process(self) -> None:
"""
Processes SODAv2025 data values as follows:
- drop the year coordinate
- Apply a mask to the dataset based on locations of NaN values.
"""
self.ds = self.ds.drop_vars("year")
condition = self.ds["dic"].isnull().any(dim=self.dim_names["time"])
self.ds["mask"] = xr.where(condition, 0, 1)
There was a problem hiding this comment.
I tested your changes in the surface_forcing notebook, and they appear to be similar or the old method ever so slightly faster... I tested each 3 times bc Anvil can be unpredictable.
using code change suggestion:
condition = self.ds["dic"].isnull().any(dim=self.dim_names["time"])
self.ds["mask"] = xr.where(condition, 0, 1)
these were the specs:
1)
CPU times: user 3.28 s, sys: 711 ms, total: 3.99 s
Wall time: 45.1 s
2)
CPU times: user 2.96 s, sys: 498 ms, total: 3.45 s
Wall time: 9.07 s
3)
CPU times: user 2.85 s, sys: 478 ms, total: 3.33 s
Wall time: 9.32 s
using the original code, these were the specs:
1)
CPU times: user 3.17 s, sys: 618 ms, total: 3.79 s
Wall time: 30 s
2)
CPU times: user 2.94 s, sys: 547 ms, total: 3.49 s
Wall time: 8.61 s
3)
CPU times: user 3.18 s, sys: 649 ms, total: 3.83 s
Wall time: 8.72 s
This PR adds the ROMS-Tool's functionality to make restoring forces files for DIC and ALK, using NOAA's OceanSODA data.
The dataset that works is:
/anvil/projects/x-ees250129/Datasets/OceanSODA-ETHZ/OceanSODA_ETHZ-v2025.OCADS.01-1982-2024.nc) that contains NOAA's OceanSODA 2025 data. 1 degree resolution (variablesdicandtalkare needed, both in units of umol/kg).The SurfaceForcing class
type = restoringneeds to be called, andsDICANDsALKneed to be specified inrestoring_forcesin order to create the file. both need to be included, or neither. Ifsss(salinity) is too be called, it needs to be called separately fromsDICandsALKsince they have different time dimensions (monthly vs climatology).pre-commit run --all-filesdocs/releases.md