From c7180ad3e234f9b266814aee5f5d315796ca44c8 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Fri, 15 May 2026 09:47:49 -0700 Subject: [PATCH 1/3] Added density-dependent interpolation of BGC tracers for BC and IC * Density-space vertical interpolation for BGC tracers in `InitialConditions` and `BoundaryForcing` via new `use_density_interpolation` parameter (default `True`). For `BoundaryForcing`, pass the physics `BoundaryForcing` as `physics_forcing=` to activate. * updated notebooks --- docs/boundary_forcing.ipynb | 2 + docs/end_to_end.ipynb | 1 + docs/releases.md | 2 + pyproject.toml | 1 + roms_tools/setup/boundary_forcing.py | 199 ++++++++- roms_tools/setup/initial_conditions.py | 75 +++- roms_tools/setup/utils.py | 101 +++++ .../tests/test_setup/test_boundary_forcing.py | 91 ++++ .../ALK/c/0/0/0/0 | Bin 194 -> 195 bytes .../ALK_ALT_CO2/c/0/0/0/0 | Bin 201 -> 195 bytes .../DIC/c/0/0/0/0 | Bin 194 -> 194 bytes .../DIC_ALT_CO2/c/0/0/0/0 | Bin 201 -> 201 bytes .../DOC/c/0/0/0/0 | Bin 201 -> 201 bytes .../DOCr/c/0/0/0/0 | Bin 201 -> 201 bytes .../DON/c/0/0/0/0 | Bin 201 -> 201 bytes .../DONr/c/0/0/0/0 | Bin 201 -> 201 bytes .../DOP/c/0/0/0/0 | Bin 201 -> 201 bytes .../DOPr/c/0/0/0/0 | Bin 193 -> 192 bytes .../Fe/c/0/0/0/0 | Bin 201 -> 201 bytes .../Lig/c/0/0/0/0 | Bin 201 -> 201 bytes .../NH4/c/0/0/0/0 | Bin 201 -> 201 bytes .../NO3/c/0/0/0/0 | Bin 201 -> 201 bytes .../O2/c/0/0/0/0 | Bin 201 -> 201 bytes .../PO4/c/0/0/0/0 | Bin 201 -> 201 bytes .../SiO3/c/0/0/0/0 | Bin 201 -> 201 bytes .../diatC/c/0/0/0/0 | Bin 201 -> 201 bytes .../diatChl/c/0/0/0/0 | Bin 201 -> 201 bytes .../diatFe/c/0/0/0/0 | Bin 201 -> 201 bytes .../diatP/c/0/0/0/0 | Bin 201 -> 201 bytes .../diatSi/c/0/0/0/0 | Bin 201 -> 201 bytes .../diazC/c/0/0/0/0 | Bin 201 -> 201 bytes .../diazChl/c/0/0/0/0 | Bin 201 -> 201 bytes .../diazFe/c/0/0/0/0 | Bin 201 -> 201 bytes .../diazP/c/0/0/0/0 | Bin 201 -> 201 bytes .../spC/c/0/0/0/0 | Bin 201 -> 201 bytes .../spCaCO3/c/0/0/0/0 | Bin 201 -> 201 bytes .../spChl/c/0/0/0/0 | Bin 201 -> 201 bytes .../spFe/c/0/0/0/0 | Bin 201 -> 201 bytes .../spP/c/0/0/0/0 | Bin 201 -> 201 bytes .../zarr.json | 416 +++++++++--------- .../zooC/c/0/0/0/0 | Bin 201 -> 201 bytes .../ALK/c/0/0/0/0 | Bin 99 -> 82 bytes .../ALK_ALT_CO2/c/0/0/0/0 | Bin 99 -> 82 bytes .../DIC/c/0/0/0/0 | Bin 89 -> 82 bytes .../DIC_ALT_CO2/c/0/0/0/0 | Bin 89 -> 82 bytes .../NO3/c/0/0/0/0 | Bin 116 -> 82 bytes .../O2/c/0/0/0/0 | Bin 115 -> 82 bytes .../PO4/c/0/0/0/0 | Bin 115 -> 82 bytes .../SiO3/c/0/0/0/0 | Bin 115 -> 82 bytes .../abs_time/zarr.json | 7 +- .../zarr.json | 9 +- .../test_setup/test_initial_conditions.py | 58 +++ roms_tools/tests/test_setup/test_utils.py | 99 +++++ 53 files changed, 827 insertions(+), 234 deletions(-) diff --git a/docs/boundary_forcing.ipynb b/docs/boundary_forcing.ipynb index 5a171685d..9de3097ba 100644 --- a/docs/boundary_forcing.ipynb +++ b/docs/boundary_forcing.ipynb @@ -3140,6 +3140,7 @@ " end_time=end_time,\n", " source={\"name\": \"CESM_REGRIDDED\", \"path\": cesm_bgc_path, \"climatology\": True},\n", " type=\"bgc\",\n", + " physics_forcing=boundary_forcing, # use density-space interpolation\n", " use_dask=True,\n", ")" ] @@ -11209,6 +11210,7 @@ " end_time=end_time,\n", " source={\"name\": \"UNIFIED\", \"path\": unified_bgc_path, \"climatology\": True},\n", " type=\"bgc\",\n", + " physics_forcing=boundary_forcing, # use density-space interpolation\n", " use_dask=True,\n", ")" ] diff --git a/docs/end_to_end.ipynb b/docs/end_to_end.ipynb index 857c95224..7a125e38f 100644 --- a/docs/end_to_end.ipynb +++ b/docs/end_to_end.ipynb @@ -11215,6 +11215,7 @@ " end_time=end_time,\n", " source={\"name\": \"UNIFIED\", \"path\": unified_bgc_path, \"climatology\": True},\n", " type=\"bgc\",\n", + " physics_forcing=boundary_forcing, # for density-space BGC interpolation\n", " use_dask=True\n", ")" ] diff --git a/docs/releases.md b/docs/releases.md index ff2dc4fcb..597d93531 100644 --- a/docs/releases.md +++ b/docs/releases.md @@ -8,6 +8,7 @@ ### New Features +* Density-space vertical interpolation for BGC tracers in `InitialConditions` and `BoundaryForcing` via new `use_density_interpolation` parameter (default `True`). For `BoundaryForcing`, pass the physics `BoundaryForcing` as `physics_forcing=` to activate ([#606](https://github.com/CWorthy-ocean/roms-tools/pull/606)). * `make_edata` changed to `make_nesting_info` * `to_yaml` and `from_yaml` were adjusted to handle child grids after they've been modified ([#573](https://github.com/CWorthy-ocean/roms-tools/pull/573)) * Nesting now supports optional baroclinic pressure fluxes via metadata ([#568](https://github.com/CWorthy-ocean/roms-tools/pull/568)) @@ -41,6 +42,7 @@ * Document `close_narrow_channels` option in `Grid` and `update_mask()`; update notebook examples * Both the surface forcing and datasets notebooks are updated to reflect `restoring` function and WOA data ([#589](https://github.com/CWorthy-ocean/roms-tools/pull/589)) * The cdr notebook is updated to reflect interpolation option. Default is same as ROMS, no interpolation ([#601](https://github.com/CWorthy-ocean/roms-tools/pull/601)) +* Notebooks updated to document density-space BGC interpolation ([#606](https://github.com/CWorthy-ocean/roms-tools/pull/606)) ### Bugfixes diff --git a/pyproject.toml b/pyproject.toml index 3d0bf6625..a07047d5f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,6 +35,7 @@ dependencies = [ "bottleneck", "regionmask>=0.11.0", "xgcm>=0.9.0", + "gsw", "numba>=0.61.2", "pydantic>2,<3", ] diff --git a/roms_tools/setup/boundary_forcing.py b/roms_tools/setup/boundary_forcing.py index 69f0cfaf4..d965b0775 100644 --- a/roms_tools/setup/boundary_forcing.py +++ b/roms_tools/setup/boundary_forcing.py @@ -25,7 +25,9 @@ add_time_info_to_ds, check_and_set_boundaries, compute_barotropic_velocity, + _compute_bgc_source_density, compute_missing_bgc_variables, + compute_potential_density, from_yaml, get_boundary_coords, get_target_coords, @@ -37,6 +39,7 @@ write_to_yaml, ) from roms_tools.utils import ( + interpolate_cyclic_time, interpolate_from_rho_to_u, interpolate_from_rho_to_v, rotate_velocities, @@ -46,6 +49,52 @@ from roms_tools.vertical_coordinate import compute_depth +def _interpolate_phys_to_bgc_time( + phys_da: xr.DataArray, + time_dim: str, + bgc_time_coord: xr.DataArray, + bgc_climatology: bool, +) -> xr.DataArray: + """Linearly interpolate a physics DataArray onto the BGC time coordinate. + + Parameters + ---------- + phys_da : xr.DataArray + Physics data with a ``datetime64`` time dimension named ``time_dim``. + time_dim : str + Name of the time dimension in ``phys_da``. + bgc_time_coord : xr.DataArray + Target time coordinate from the BGC dataset (1-D). + bgc_climatology : bool + Whether the BGC dataset is a climatology. If True, ``bgc_time_coord`` + is expected to be ``timedelta64`` from the start of the year (as set by + ``assign_dates_to_climatology``), and physics times are mapped to + fractional day-of-year before a cyclic linear interpolation. If False, + a straight ``xr.DataArray.interp(method="linear")`` is performed in + ``datetime64`` space. + + Returns + ------- + xr.DataArray + ``phys_da`` interpolated to ``bgc_time_coord``, with time dimension + still named ``time_dim`` and coordinate set to ``bgc_time_coord``. + """ + if bgc_climatology: + # Map both axes to fractional day-of-year, then cyclic-interpolate. + bgc_doy = (bgc_time_coord / np.timedelta64(1, "D")).values + 1.0 + phys_doy = phys_da[time_dim].dt.dayofyear.values.astype(float) + phys_for_interp = phys_da.assign_coords( + {time_dim: xr.DataArray(phys_doy, dims=[time_dim])} + ) + result = interpolate_cyclic_time( + phys_for_interp, time_dim, time_dim, bgc_doy + ) + return result.assign_coords({time_dim: bgc_time_coord.values}) + + # Non-climatology: standard datetime64 linear interpolation + return phys_da.interp({time_dim: bgc_time_coord}, method="linear") + + @dataclass(kw_only=True) class BoundaryForcing: """Represents boundary forcing input data for ROMS. @@ -104,6 +153,18 @@ class BoundaryForcing: Indicates whether to skip validation checks in the processed data. When set to True, the validation process that ensures no NaN values exist at wet points in the processed dataset is bypassed. Defaults to False. + use_density_interpolation : bool, optional + If True (default), BGC tracers are vertically interpolated in density space + rather than depth space when ``physics_forcing`` is provided. Preserves + water-mass properties. Only used when ``type='bgc'``. + + Interpolation uses ``xgcm.Grid.transform`` with the linear method inside + the source density range and edge-value (nearest-neighbor) extrapolation + outside (``mask_edges=False``). + physics_forcing : BoundaryForcing, optional + A physics ``BoundaryForcing`` object (``type='physics'``) whose T/S fields + supply the density coordinate for BGC tracer interpolation. When None and + ``use_density_interpolation=True``, falls back to depth-based interpolation. Examples @@ -143,6 +204,10 @@ class BoundaryForcing: """Optional initial bounding slice when loading source data (Dask); see dataset classes.""" bypass_validation: bool = False """Whether to skip validation checks in the processed data.""" + use_density_interpolation: bool = True + """Whether to interpolate BGC tracers in density space rather than depth space.""" + physics_forcing: "BoundaryForcing | None" = None + """Physics BoundaryForcing object supplying T/S for density-based BGC interpolation.""" ds: xr.Dataset = field(init=False, repr=False) """An xarray Dataset containing post-processed variables ready for input into @@ -156,9 +221,20 @@ def __post_init__(self): # Initialize depth coordinates self.adjust_depth_for_sea_surface_height = False self.ds_depth_coords = xr.Dataset() + self._phys_bdry_depth_data: dict = {} self._input_checks() + if ( + self.type == "bgc" + and self.use_density_interpolation + and self.physics_forcing is None + ): + logging.info( + "use_density_interpolation=True but no physics_forcing provided. " + "BGC tracers will be interpolated in depth space instead." + ) + target_coords = get_target_coords(self.grid) data = self._get_data() @@ -329,6 +405,16 @@ def __post_init__(self): zeta_u = 0 zeta_v = 0 + # Save physics T/S at source depth levels for density-based BGC interpolation + if self.type == "physics" and self.use_density_interpolation: + if "temp" in processed_fields and "salt" in processed_fields: + self._phys_bdry_depth_data[direction] = { + "temp": processed_fields["temp"], + "salt": processed_fields["salt"], + "depth_dim": bdry_data.dim_names["depth"], + "depth_coord": bdry_data.ds[bdry_data.dim_names["depth"]], + } + for location in ["rho", "u", "v"]: # Filter var_names by location and check for 3D variables filtered_vars = [ @@ -353,17 +439,111 @@ def __post_init__(self): vertical_regrid = VerticalRegrid( bdry_data.ds, source_dim=bdry_data.dim_names["depth"] ) + + use_density = ( + self.type == "bgc" + and self.use_density_interpolation + and location == "rho" + and self.physics_forcing is not None + and direction in self.physics_forcing._phys_bdry_depth_data + ) + + if use_density: + phys_data = self.physics_forcing._phys_bdry_depth_data[ + direction + ] + phys_temp = phys_data["temp"] + phys_salt = phys_data["salt"] + bgc_climatology = self.source["climatology"] + + # BGC tracers from any 3D BGC variable share the same + # time axis; grab it from the first one. + time_dim = bdry_data.dim_names.get("time") + first_bgc = ( + processed_fields.get(filtered_vars[0]) + if filtered_vars else None + ) + if ( + first_bgc is not None + and time_dim + and time_dim in first_bgc.dims + ): + bgc_time_coord = first_bgc[time_dim] + else: + bgc_time_coord = None + + if bgc_time_coord is not None and "time" in phys_temp.dims: + phys_temp = _interpolate_phys_to_bgc_time( + phys_temp, "time", bgc_time_coord, bgc_climatology + ) + phys_salt = _interpolate_phys_to_bgc_time( + phys_salt, "time", bgc_time_coord, bgc_climatology + ) + elif "time" in phys_temp.dims: + # No BGC time found — collapse to mean as fallback + phys_temp = phys_temp.mean("time") + phys_salt = phys_salt.mean("time") + + source_density = _compute_bgc_source_density( + phys_temp, + phys_salt, + phys_data["depth_dim"], + phys_data["depth_coord"], + bdry_data.ds[bdry_data.dim_names["depth"]], + bdry_data.dim_names["depth"], + ) + + # Target density: physics sigma T/S interpolated to BGC time. + # Physics BC stores both bry_time (float days since reference) + # and abs_time (datetime64). We swap to the datetime64 view so + # _interpolate_phys_to_bgc_time sees a real time axis. + temp_sigma = self.physics_forcing.ds[f"temp_{direction}"] + salt_sigma = self.physics_forcing.ds[f"salt_{direction}"] + if bgc_time_coord is not None and "abs_time" in temp_sigma.coords: + temp_sigma = ( + temp_sigma.swap_dims({"bry_time": "abs_time"}) + .rename({"abs_time": "time"}) + ) + salt_sigma = ( + salt_sigma.swap_dims({"bry_time": "abs_time"}) + .rename({"abs_time": "time"}) + ) + temp_sigma = _interpolate_phys_to_bgc_time( + temp_sigma, "time", bgc_time_coord, bgc_climatology + ) + salt_sigma = _interpolate_phys_to_bgc_time( + salt_sigma, "time", bgc_time_coord, bgc_climatology + ) + elif "bry_time" in temp_sigma.dims: + temp_sigma = temp_sigma.mean("bry_time") + salt_sigma = salt_sigma.mean("bry_time") + + target_density = compute_potential_density( + temp_sigma, salt_sigma + ) + n_s = target_density.sizes["s_rho"] + target_density = target_density + xr.DataArray( + np.arange(n_s) * 1e-7, dims=["s_rho"] + ) + for var_name in filtered_vars: if var_name in processed_fields: - processed_fields[var_name] = vertical_regrid.apply( - processed_fields[var_name], - source_depth_coords=bdry_data.ds[ - bdry_data.dim_names["depth"] - ], - target_depth_coords=self.ds_depth_coords[ - f"layer_depth_{location}_{direction}" - ], - ) + if use_density: + processed_fields[var_name] = vertical_regrid.apply( + processed_fields[var_name], + source_depth_coords=source_density, + target_depth_coords=target_density, + ) + else: + processed_fields[var_name] = vertical_regrid.apply( + processed_fields[var_name], + source_depth_coords=bdry_data.ds[ + bdry_data.dim_names["depth"] + ], + target_depth_coords=self.ds_depth_coords[ + f"layer_depth_{location}_{direction}" + ], + ) # compute barotropic velocities if "u" in var_names and "v" in var_names: @@ -1091,6 +1271,7 @@ def to_yaml(self, filepath: str | Path) -> None: "ds_depth_coords", "adjust_depth_for_sea_surface_height", "use_dask", + "physics_forcing", ], ) write_to_yaml(forcing_dict, filepath) diff --git a/roms_tools/setup/initial_conditions.py b/roms_tools/setup/initial_conditions.py index b601f655c..1267d3cbb 100644 --- a/roms_tools/setup/initial_conditions.py +++ b/roms_tools/setup/initial_conditions.py @@ -26,8 +26,10 @@ ) from roms_tools.setup.utils import ( RawDataSource, + _compute_bgc_source_density, compute_barotropic_velocity, compute_missing_bgc_variables, + compute_potential_density, from_yaml, get_target_coords, get_variable_metadata, @@ -111,6 +113,17 @@ class InitialConditions: Indicates whether to skip validation checks in the processed data. When set to True, the validation process that ensures no NaN values exist at wet points in the processed dataset is bypassed. Defaults to False. + use_density_interpolation : bool, optional + If True (default), BGC tracers are vertically interpolated in density space + rather than depth space. Density is computed from the physics source T/S + (via TEOS-10 sigma-0) and used as the vertical coordinate for interpolation, + preserving water-mass properties. Only applies when ``bgc_source`` is provided + and the physics source is a lat/lon dataset (not a ROMS restart). Falls back + to depth-based interpolation silently if physics T/S are unavailable. + + Interpolation uses ``xgcm.Grid.transform`` with the linear method inside + the source density range and edge-value (nearest-neighbor) extrapolation + outside (``mask_edges=False``). @@ -161,6 +174,8 @@ class InitialConditions: """Optional initial bounding slice when loading lat/lon forcing data with Dask.""" bypass_validation: bool = False """Whether to skip validation checks in the processed data.""" + use_density_interpolation: bool = True + """Whether to interpolate BGC tracers in density space rather than depth space.""" ds: xr.Dataset = field(init=False, repr=False) """An xarray Dataset containing post-processed variables ready for input into ROMS.""" @@ -173,6 +188,9 @@ def __post_init__(self): # Initialize depth coordinates self.ds_depth_coords = xr.Dataset() + # Populated during physics processing if density interpolation is enabled + self._phys_depth_data: dict | None = None + self._input_checks() processed_fields = {} @@ -273,8 +291,17 @@ def _process_data(self, processed_fields, type="physics"): for location in ["rho", "u", "v"]: self._get_depth_coordinates(zeta, location, "layer") + # Save physics T/S at source depth levels for density-based BGC interpolation + if type == "physics" and self.use_density_interpolation and isinstance(data, LatLonDataset): + self._phys_depth_data = { + "temp": processed_fields["temp"], + "salt": processed_fields["salt"], + "depth_dim": data.dim_names["depth"], + "depth_coord": data.ds[data.dim_names["depth"]], + } + # Vertical regridding - processed_fields = self._regrid_vertically(data, processed_fields, var_names) + processed_fields = self._regrid_vertically(data, processed_fields, var_names, type=type) # Compute barotropic velocities if "u" in var_names and "v" in var_names: @@ -353,6 +380,7 @@ def _regrid_vertically( data: ROMSDataset | LatLonDataset, processed_fields: dict[str, xr.DataArray], var_names: dict[str, dict[str, str | bool]], + type: str = "physics", ) -> dict[str, xr.DataArray]: """ Perform vertical regridding of 3D variables to the model's vertical grid. @@ -426,16 +454,47 @@ def _regrid_vertically( data.ds, source_dim=data.dim_names["depth"] ) + use_density = ( + type == "bgc" + and self.use_density_interpolation + and self._phys_depth_data is not None + ) + + if use_density: + source_density = _compute_bgc_source_density( + self._phys_depth_data["temp"], + self._phys_depth_data["salt"], + self._phys_depth_data["depth_dim"], + self._phys_depth_data["depth_coord"], + data.ds[data.dim_names["depth"]], + data.dim_names["depth"], + ) + target_density = compute_potential_density( + processed_fields["temp"], processed_fields["salt"] + ) + # Add small perturbation to target density to ensure monotonicity + n_s = target_density.sizes["s_rho"] + target_density = target_density + xr.DataArray( + np.arange(n_s) * 1e-7, dims=["s_rho"] + ) + for var_name in filtered_vars: if var_name not in processed_fields: continue - processed_fields[var_name] = vertical_regrid.apply( - processed_fields[var_name], - source_depth_coords=data.ds[data.dim_names["depth"]], - target_depth_coords=self.ds_depth_coords[ - f"layer_depth_{location}" - ], - ) + if use_density: + processed_fields[var_name] = vertical_regrid.apply( + processed_fields[var_name], + source_depth_coords=source_density, + target_depth_coords=target_density, + ) + else: + processed_fields[var_name] = vertical_regrid.apply( + processed_fields[var_name], + source_depth_coords=data.ds[data.dim_names["depth"]], + target_depth_coords=self.ds_depth_coords[ + f"layer_depth_{location}" + ], + ) return processed_fields diff --git a/roms_tools/setup/utils.py b/roms_tools/setup/utils.py index d952ae5c9..837685b94 100644 --- a/roms_tools/setup/utils.py +++ b/roms_tools/setup/utils.py @@ -10,6 +10,7 @@ from pathlib import Path from typing import Any, Literal, TypeAlias +import gsw import numba as nb import numpy as np import pandas as pd @@ -18,6 +19,7 @@ from pydantic import BaseModel from roms_tools.constants import R_EARTH +from roms_tools.regrid import VerticalRegrid if typing.TYPE_CHECKING: from roms_tools.setup.grid import Grid @@ -498,6 +500,105 @@ def compute_missing_bgc_variables(bgc_data): return bgc_data +def compute_potential_density( + temp: "xr.DataArray", salt: "xr.DataArray" +) -> "xr.DataArray": + """Compute sigma-0 potential density anomaly (kg/m³ - 1000) via TEOS-10 (gsw). + + Wraps gsw.sigma0 with apply_ufunc for dask compatibility. Treats practical + salinity as Absolute Salinity and in-situ temperature as Conservative + Temperature — an approximation sufficient for density-coordinate interpolation. + + Parameters + ---------- + temp : xr.DataArray + In-situ temperature (°C). + salt : xr.DataArray + Practical salinity (PSU). + + Returns + ------- + xr.DataArray + Potential density anomaly sigma-0 (kg/m³ - 1000). + """ + return xr.apply_ufunc( + gsw.sigma0, + salt, + temp, + dask="parallelized", + output_dtypes=[temp.dtype], + ) + + +def _compute_bgc_source_density( + phys_temp: "xr.DataArray", + phys_salt: "xr.DataArray", + phys_depth_dim: str, + phys_depth_coord: "xr.DataArray", + bgc_depth_coord: "xr.DataArray", + bgc_depth_dim: str, +) -> "xr.DataArray": + """Compute potential density at BGC source depth levels. + + Computes sigma-0 from physics T/S at physics depth levels, adds a small + depth-index perturbation to guarantee strict monotonicity (matching the + reference MATLAB implementation), then interpolates to the BGC depth grid. + + Parameters + ---------- + phys_temp : xr.DataArray + Physics temperature at source depth levels. + phys_salt : xr.DataArray + Physics salinity at source depth levels. + phys_depth_dim : str + Name of the depth dimension in the physics dataset. + phys_depth_coord : xr.DataArray + Depth coordinate values of the physics dataset (1D, metres). + bgc_depth_coord : xr.DataArray + Depth coordinate values of the BGC dataset (1D, metres). + bgc_depth_dim : str + Name of the depth dimension in the BGC dataset. + + Returns + ------- + xr.DataArray + Potential density anomaly sigma-0 at BGC depth levels. + """ + density = compute_potential_density(phys_temp, phys_salt) + + # Add small perturbation along depth to ensure strict monotonicity + n_depth = density.sizes[phys_depth_dim] + perturbation = xr.DataArray(np.arange(n_depth) * 1e-7, dims=[phys_depth_dim]) + density = density + perturbation + + # Regrid density from physics depth levels to BGC depth levels + ds_phys = xr.Dataset({phys_depth_dim: phys_depth_coord}) + vertical_regrid = VerticalRegrid(ds_phys, source_dim=phys_depth_dim) + source_density = vertical_regrid.apply( + density, + source_depth_coords=phys_depth_coord, + target_depth_coords=bgc_depth_coord, + ) + + # xgcm.transform names the output dim after `target` (bgc_depth_coord), + # but rename defensively in case the two arguments diverge. + if ( + phys_depth_dim != bgc_depth_dim + and phys_depth_dim in source_density.dims + ): + source_density = source_density.rename({phys_depth_dim: bgc_depth_dim}) + + # Add a small perturbation along the BGC depth dimension after interpolation, + # so the density profile xgcm uses as a source coordinate is strictly + # monotonic even where extrapolation flattens the original perturbation. + n_bgc_depth = source_density.sizes[bgc_depth_dim] + source_density = source_density + xr.DataArray( + np.arange(n_bgc_depth) * 1e-7, dims=[bgc_depth_dim] + ) + + return source_density + + def compute_missing_surface_bgc_variables(bgc_data): """Fills in missing surface biogeochemical (BGC) variables in the input dictionary. diff --git a/roms_tools/tests/test_setup/test_boundary_forcing.py b/roms_tools/tests/test_setup/test_boundary_forcing.py index a79f4919f..a071b039c 100644 --- a/roms_tools/tests/test_setup/test_boundary_forcing.py +++ b/roms_tools/tests/test_setup/test_boundary_forcing.py @@ -736,3 +736,94 @@ def test_nondefault_glorys_dataset_loading(small_grid: Grid, use_dask: bool) -> expected_vars = {"u_south", "v_south", "temp_south", "salt_south"} assert set(bf.ds.data_vars).issuperset(expected_vars) + + +# test density interpolation + + +def test_bgc_bc_density_fallback_without_physics_forcing( + bgc_boundary_forcing_from_unified_climatology, +): + """BGC BC with density interpolation but no physics_forcing falls back to depth-based.""" + bf = bgc_boundary_forcing_from_unified_climatology + assert bf.use_density_interpolation is True + assert bf.physics_forcing is None + # BGC variables should still be present (depth-based fallback succeeded) + assert any("NO3" in v for v in bf.ds.data_vars) + + +def test_bgc_bc_with_physics_forcing(use_dask): + """BGC BC with physics_forcing uses density interpolation and produces BGC variables.""" + # Use same grid / data as existing physics BC fixtures (North Atlantic, 2012) + grid = Grid( + nx=3, + ny=3, + size_x=400, + size_y=400, + center_lon=-8, + center_lat=58, + rot=0, + N=3, + theta_s=5.0, + theta_b=2.0, + hc=250.0, + ) + fname_phys = Path(download_test_data("GLORYS_NA_20120101.nc")) + fname_bgc = Path(download_test_data("coarsened_UNIFIED_bgc_dataset.nc")) + + physics_bc = BoundaryForcing( + grid=grid, + start_time=datetime(2012, 1, 1), + end_time=datetime(2012, 1, 2), + source={"path": fname_phys, "name": "GLORYS"}, + type="physics", + apply_2d_horizontal_fill=False, + use_dask=use_dask, + ) + + bgc_bc = BoundaryForcing( + grid=grid, + start_time=datetime(2012, 1, 1), + end_time=datetime(2012, 1, 2), + source={"path": fname_bgc, "name": "UNIFIED", "climatology": True}, + type="bgc", + physics_forcing=physics_bc, + apply_2d_horizontal_fill=True, + use_dask=use_dask, + ) + + assert bgc_bc.use_density_interpolation is True + assert bgc_bc.physics_forcing is physics_bc + for direction in ["south", "east", "north", "west"]: + if bgc_bc.boundaries[direction]: + assert f"NO3_{direction}" in bgc_bc.ds + assert f"DIC_{direction}" in bgc_bc.ds + + # Verify the density code path actually changed the output: produce a + # depth-based comparison and assert at least one BGC variable differs. + bgc_bc_depth = BoundaryForcing( + grid=grid, + start_time=datetime(2012, 1, 1), + end_time=datetime(2012, 1, 2), + source={"path": fname_bgc, "name": "UNIFIED", "climatology": True}, + type="bgc", + use_density_interpolation=False, + apply_2d_horizontal_fill=True, + use_dask=use_dask, + ) + any_diff = False + for direction in ["south", "east", "north", "west"]: + if not bgc_bc.boundaries[direction]: + continue + for var in ["NO3", "DIC", "ALK", "PO4", "O2"]: + name = f"{var}_{direction}" + if name in bgc_bc.ds and name in bgc_bc_depth.ds: + a = bgc_bc.ds[name].values + b = bgc_bc_depth.ds[name].values + valid = ~(np.isnan(a) | np.isnan(b)) + if valid.any() and np.abs(a[valid] - b[valid]).max() > 0: + any_diff = True + break + if any_diff: + break + assert any_diff, "Density interpolation produced identical output to depth-based" diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/ALK/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/ALK/c/0/0/0/0 index 36126fca0bf3ecfd4110c358b791338199155c8a..f398e4427ac27e373eada2d81d388fb34309a4f1 100644 GIT binary patch literal 195 zcmV;!06hOFwJ-f3z|{o+0*oy*P*u_x`SDwfW`d#z$P?A6+OSHPPY@t`Di7$30=#g@ zGDL3r3+rwI5iS`Bw%k&p6Zrw)>o$689*zmt@iO_bRIHuIup@gTEPf7f$|P|`19OIu zR+R2kvCdOz*~IRb=^1KR=t7;8wzFw5QxO~7@%@2*d`G8R*-3n+j%>BhFH!NPw9%JaGD^jEd%EZXa#6k xJwofx=#A;#n~7Ja@rSOyCqnKFAbBEB_Kj;dv58K?4#+AG#4m9Ip~ZD5DWYHRTXR z68jBBPE!s=p>Y&NmzWbpXMz$%UBzn2n4sB;oU<^2&w9<~uh zdaM#f)jMZ({K_+ DNPI?S diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/DIC/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/DIC/c/0/0/0/0 index 2f75425ecd87b9a33f27e7e2d284a44364e37680..a2c55c723bb2ad97876fff506c570d927e80d343 100644 GIT binary patch literal 194 zcmV;z06qUGwJ-f3z|92!0!%G3KvmKhYzeKb0b;vKRbTu50rS9W z0e}NEoMB>$3+Skk=$TEVlmRPs^L9AATg*G6&WEvbJS#}bDdYjU3;@QTmeZsQl3lj38%&17nfFsj_oU^sF=2N)WONXq*|% wdvX0HnIDk$(~0W<%WF}*A_*_V3+!AzU>8uGLltKU%)5}Ya!jj)!bvm|0ImX8?*IS* literal 194 zcmV;z06qUGwJ-f3z|92!0!%G5KvR+qoxhkCA0_5s3a0wHCUSkeAQpyw3rpoptOr5p zG$}FYP?W)kjFJ4s``Z-jN!xM<*TLirK{#JPQ*{|V5m%eJBc~B7cv=bIWjHM(sbRdh zE{noVYur-Clh6f`I87B(9M=?SWqg>p8R*_DT$1W(GVWzN8zahT)G*x`oqO=P1S(^x w=JcT083^obsvS+n9Vj^&U``S1wD@@2g}p7HW6Tkaqt6Kux}wp=NEk6I08w>TwEzGB diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/DIC_ALT_CO2/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/DIC_ALT_CO2/c/0/0/0/0 index 19827b2ee6654d220259e654654d1f2a1f3c3276..b1688c6fc564ccd4a4d7076303a3abeeed792762 100644 GIT binary patch delta 177 zcmV;i08anO0m%V>D77#BAix0z0Hs9#L@pZuMF2biMcM}eMJsv!M9U!mMA10>M2-Rd zL=QFrMXO!`MW^lnMLc@{L}S+jMT6i2MF$B3MfpGhMHaaJMCyY7L|mi)L_L-OMf@KA zL}-uxMBF3$L~|PaMB{1zMXj&^MRlqGMY+2DL=Kk%MH^`d0!4z;0Yw}>k%uB6uzLSQ f&{Y6ML#qHqT2TE&_wE5j0H^^)bZG%aT<8BpJmND77#BAix0z0JB2>L>d|ZMR+;@MOX#_Mf7+6MAshwL=H3jMEd>w zM8-7%Mfq0&MP%pzMa*sgM7hudMd8>3MJoUUMF&3tMLf9vL{fwQM4F`kMDdmYMF=1M zM2V07L@FfvM4KD@LD77#BAix0z0DK7?Ld-8FLPh`_LJT(UK}g~#LTkZDLOk+-LRh+C zLdRFeL1{wdLD^ggLi6V)Ldx`lL2>kgLAQ5}K}urLLH#2iLbKi|Lg0QLLQ?ASK_fgX zLTxWjLZ;+|Ld$(B(vf&d#_K!6}rCf(W`xuBpcRG_nW7FY5^XE?d delta 177 zcmV;i08anO0m%V>D77#BAix0z0F@LRLOeJoLhK0~LRv%aL4fKgLQvF5LRSZYLLt^+ zLPA=`K^Da1K}t3WLK90TLXi)HLBsE(BK@o^S)H9Po#=zl0W#C7x diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/DOCr/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/DOCr/c/0/0/0/0 index 87dd2dac58df93919c1b68570106b0eec64bc745..a21cc9a38290b8912fa1e36c553a39d668f2bf06 100644 GIT binary patch delta 177 zcmV;i08anO0m%V>D77#BAix0z0G1%uL0xRvK|L+mLABM@K{RLAK@H#AK^o=cL7bK2 zLCy2hK}`YFLG=jML2IqrLDaO*K_lVML4`rkLCMh5L1b>$K|9RYL2i=SLC+u7LCMC~ zL5-i>L53RULCTZmLBey>L0+)bK`ZUnK_;}^K}fI9L3G&&&q3=<&_QvCk%uB6Jl@km fQ{dD=d-m2r?g-sML#NL{SJ2Nv2uaXEx2@Dc(}h~m delta 177 zcmV;i08anO0m%V>D77#BAix0z0Iwm|L5^+MLCr4NLB`hAL5FD9L9pT5L7eO4LD`+- zLEiJyL4gL;L3$h3L2|s>LDIC(K?vc`LF+-#K_Jr8L1=E)K}gKlL3oncLDnDELC?n5 zL6x7}L0}r@K}D10K@)Y;L8P(NL3r-gK@PUuK|HU|L1ozp&q41@&_P~_k%uB6THVt@ fw&2u33HR1PiwE67zNgPY^3cyge@W0mVXf3b5OiHI diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/DON/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/DON/c/0/0/0/0 index b5ddfba7d7d1a9c4e798c73e2a76fe17915ba672..c1ad3d5d1b2dafc3466dbc74525eee832407181c 100644 GIT binary patch delta 177 zcmV;i08anO0m%V>D77#BAix0z07s2$K#WF+KnPldK#ABrK&;1tKwTrMK&9pJKw|F1 zKv%Q{K(3K0K-L^nKwZv+Kt%qrKQlVAKhHwLKibw9K%V1pK#jGFK)Kk5K$;XnK=@^a zKr=3}Km{`XK>35pKy>5@KmlboKx%MiKvYPNK=3rTKb1QOxIZD$%0J`bk%uB6hqnws fa7{Tts1;~H$+nO{7m~O?O5nLa89vNE8TBMUvY||n delta 177 zcmV;i08anO0m%V>D77#BAix0z0A8SLKv+?TK#pgGK*i@hK-AEJK;1j4KvXC4KvD(8 zKpM6MKKwBoVKM!WIKf!>*KLqO-KnLS+KqKqpX-K#n-KKmI@nxIanN%0JEQk%uB6jI9hn fSw}fQ%@SxpSFeyjeu=n0!_>Jygfq-PO6Vj&P`pa0 diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/DONr/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/DONr/c/0/0/0/0 index 5cbfbcae1755dd7bb5733d4f0d995561d23a9a54..81c4251ec66f512ceaf6cb6552ca295151d02b76 100644 GIT binary patch delta 177 zcmV;i08anO0m%V>D77#BAix0z0LG5DKSdtCKM^0mKMM7?KW`nnKaIJ?KcN%cKkO0E zKdyGKKlLEAKUEF5KjfvpKY)gD77#BAix0z056cYKR6-2KR6@6KNeKd8XPKV&J~Kl>Qa zKe%?UKba}CKUFZeKN`2bKe>jjKjo~gKO;7-Kla_SKMC`D77#BAix0z05*mUKK;rSKFY-iKER69J$WDT61K0l2iK9~}6J#BM*J#|5hJsIiBJ-Lw%KA$@mKIv2mJ^-xLJ^!{I zKBD77#BAix0z0EmqYK9A28J{!mgK9q~pJqjirKB}oVK4rOAK1EwS zKBR1*J#moOJ?M%DK7!IAJ_jjsJ!XG=J)S;{J>m7rJ%*7EKEpc~KA=?zJ`}FhJ^QvE zK9Y|*K8U$lK5uM3KBd^6J@m8RJ$IiEJ{Wo`J_(F*J)>I)c|GLXiapTZk%uB6#Fd;q f2^`=(0Amn76uK)uRAp{GS>$v*O=^ifS-H_YlA%ZB diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/DOPr/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/DOPr/c/0/0/0/0 index f379120b3ef9e7032e0b45c09f78278fa4799b4a..1beabfdcf89a754203c84a8a6c380997504efbf4 100644 GIT binary patch literal 192 zcmV;x06+gIwJ-f3z`X?k0?aHkKvUBg7^FT}Ah7#=GTOdTAg2V-pcfEvPVsr|#r>+9 zvE~g)j*5H%7x);(eEPPoA40;!6;x? diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/Fe/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/Fe/c/0/0/0/0 index 0502aaf6bd1791ad8a4ce4537e23b76406e8cb44..360d8503b919720768cc5df66c474cf95491dcf4 100644 GIT binary patch literal 201 zcmV;)05<<9wJ-f3zySsTPusIPRQA6*SWKj4ykI!#9{xPOHv;5*ZDl)=3D77#BAix0z0Gnm5I(+D=Ix@wkIv6^uI;%FTI&%c1It+4=Iu2u& zIzs2hI!jErI>B?PIvbv#I$`j_I-&)`I?e^XI^g)NIx)DYIyC#EIz>65I6x;khXv^s#$k%uB6fBd>S f`&Fqr78sj4QU#DYx8=4vn98&|zw@s;ny;NYtbD77#BAix0z0NYWnI!)86IzX|eI@u|!I?f%dI+Ew3I{HqLI{rzP zIyC0SI&KlTI&QJ3IuJ>qI@H#}I!fchI^D;-Iu_-vI-a|zI#K|nI*~-7I!{@qI_T`9 zIz5V;I#;WTI&a{QI)Bo^I_CkeI@@Bc5W%ZppaNtl? diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/NH4/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/NH4/c/0/0/0/0 index 5a64a2bcb6a2cd5527ae4c644a1305396ec80677..4dcf5f48b2b9fab4db43986e249942d2d823d12a 100644 GIT binary patch literal 201 zcmV;)05<<9wJ-f3zySsTSJL=B*KA`xCGmkis*{yI7-!Bsu%9MAOkyTJJD(Ola(R|N zD?CX*dHZ=j?hFb(k{~fZzxG!@mn>XA;=U6?(;PD`Ucf-4(763`4imtu82r^}H& zeS@q$$g=J{=-mxH7OKxY*Oo3lKdi7lwsV?2CP^-IHF-bqrQMX-d2=8NSDk#`2Hk5?x`a_=_(aId~udP zf2T%2*70UO;ywmG^f@p;M&DLHDzsQX6Y~;3W;`i9U|^sj0|`9|13utEbq3i~BK z^sSja1%T*2@9vd9#i46IyIcZ2VQ!Z_3DBNELlTTYw{18;%xOeF&9*K+j*vP(c@K?0 zcsjyAr5@-#0YVTy%<79hBf8N(GECh*t D0P0wM diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/O2/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/O2/c/0/0/0/0 index dfb137c05a28b9ef92d8852ab59adec188b9fa18..c8484a9ea69825dcd5525851a997fcca9a16df2e 100644 GIT binary patch delta 184 zcmV;p07w7H0m%V>D77#BAix0z0O|mWL;Jy!Ln>pJLugBxL(HCyLjrh{Lmy0*LyE(h zLuRCoLu^fyLw7)#L&!CoLu6NrL#Wb@LjquwLrwUXL(}kxL;7cuLlE|tL&)QrL;Eg_ zLtrD3LxVDuLkk3$L#e-!L;gmULrl|^L&GtcLwA{zL%b#rm_x&dn?w8#n?uX5h(nr@ mjvF9}x06HG;*~>w(v(9=H+o0vme(VauiMV>=CPMbrJXH$Rx delta 184 zcmV;p07w7H0m%V>D77#BAix0z05tuHL)gHQLmg$8L%2|xL)n>)LjrJx06FC;FUuR*_1;vK$ka>7o7+{J}U@6oUztEH6xNfG(8qR5cE7g-4RhfJ_l(&N4y|DcaSGO z@t_($n)xg~=!--?AA(Lk67?`Xzk5DDA)>@S{Y|(&Oy-V0Tm42pcs&$8_nbC9od8Kb z;TT>%MMob#2XrMq>uec5LiQ*=3eN&Q_ggAHcD*k?cJMnsl;@m2Jw=N?l(cL1+E6@i&&O--3-k;Pytj3W(ra=}yNCrJV#7uec5LiQ*=V=V+enMNx zd`F=`=^&0kHdUoSYi`IudC&)3L6(s~Ln((qH{gXpBzLMnk*O3O4epT D``lge literal 201 zcmV;)05<<9wJ-f3zySsT)g6mK&FPsycY~EcdfJRYDeIU(5!9qWz8|PSes!=wm!Psh zqmQ6KGhdBBhFPUR(1^!Cjts*C&)3fAEn&mt%)Ot>1+}7ip?M;(&-ivaozW4F6t0uMcNH D+N4{U diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/diatC/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/diatC/c/0/0/0/0 index 852c80771414bd989308e6ec6618a6c3cf6e35b1..23de58025ee13b49a5ce3d056609453046aaccf9 100644 GIT binary patch literal 201 zcmV;)05<<9wJ-f3zySsTTRW;hbft_yBgWo98Tl4LVx?j~L7KlmV@xeSl$#|$)`tv1 zlqc3eERLZ++$A|b^&%2M2D=eKrA#Y9WKOt1ShM9ncN4fk@iGiSKe|6b)}ER_TPx{5 z3Rp2e(ymZIx%fIk`oIuD!D{?JIi7<*mEed$&0ms136_#UQ{oOmRMGE0(ay&}`?3~6 zOlM9(tEQws@jCZEM9WEF7+Wn D7{XW+ literal 201 zcmV;)05<<9wJ-f3zySsT;S8%kaVU*Ig^u4qpnVrX;J#x&h9$s1gjFp+Svw{`PLvEm z7Yf%v-r}i0H~Bk1p8gX-1nUz)-$^b(ZuYo9)|}-(MCP|ZtK$nn`BgtbsFs>PRTJqy zvNSP2?Keo*RMGE0)6U00`?3~6 zPG?R*tfr(t@jCZENFzBvwbfrh$0}z*!ap!UtyTL!EuMux{=c(97+1wXM?0`V96cdH D3g}lu diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/diatChl/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/diatChl/c/0/0/0/0 index 8075730e01a325cda5bf4d4dff34cd6915b1f87c..9ebe3f48329c4865ce0a1245f1eb2095b74819d8 100644 GIT binary patch literal 201 zcmV;)05<<9wJ-f3zySsTNdCA!d*_rtQKRZVUa22I=|pfoC?&`~6Ok%D>ozVw>x3CV zX?){9erLBnPT^KQD^wRiJ=YjOAhh~E{Ca{uPtuM+0nVC09rKt#!M+wiUmVsx1p~Z4t}PEh z*$YKL;wYFtbI;&D?J6ZcyKGcHA%I^%R-`CEutefM)W2>%vRSV{Ij6utc&VvCV0IWl Dx+qoS literal 201 zcmV;)05<<9wJ-f3zySsTE`hl|4_lQ#{9EfkkR>2MhC6XSUyH~-NPQ|kZp$t|?vWWl zL389kR{g&|H)L5pBR?5HPb3>aYX~+#5-7+&NuS(4p1Qd|M^Ftw-Ah71YxbHxM(*T3 zV~Hj{=T%Za;RZ@T2`CsqD?0f;c2j~r%`1;Umra{M3oe;J`uP?>UmVsx1p~Z4t}PEh z*$YKL;wYFta?jvC?J6ZcyKGcH8^d5gtlua=q4nZEmd9^Czoo7~udlyA&SOjyY~PjZow_Wzuaqi4=@ETZf!E!>w{TaWpbG)Vj(zYmpN< zGzvyJWa7g&=))5@-}ibqqk~;I1^1(1mk4W)2G+Q4zH>H<1E-3{#RmFriNhCEmodGL3 zyYd1#;I$k$p)ja7AWjlFLGmd%t^-&(IU652ByTP_h=ss7daDdMw0S%^c+|W%#Xt`@ z{55tr-sw{~m=0Vy)AlDh!!B+(@kHb|6;5tB?rn@YBetA5EG~07aWpbG)w;?!Y>^W= zH3~*KWa7g&=))5@;P-krq=Q{JC6IzSM2iZyLGQC{FT=j4((XhEMw^m DYwk}E diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/diatP/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/diatP/c/0/0/0/0 index b166c15851cbbb86c4915a54a1ea4e791d59ffae..3c8e3a8bf83c0c1bcb15daf7d1453c3e1b7abb28 100644 GIT binary patch literal 201 zcmV;)05<<9wJ-f3zySsTRm&MXp=|^`q)SpgNbrt5<;Bc96EiJ5L=&7l+gO}Dx%zlL z2<1gR8T$}CYnQ7#Q5KCo7Br7N(gCPFgj6j({nku8r>rYIrayx{YuB?q<7^2$e92He z`=_8g!BWFK!ZD*g#E64EOyFNU-=Oe2P!a(?N;L{T!bu7~6YzpP5d~H}YAQNCIqi%+ zXJ5fR`A!l%09{=?I?bm$I7!hwm9NM>3b~&>?W$fpwmtMaCw?kEMus~+;F~5s&%Bd8 DXE;*~ literal 201 zcmV;)05<<9wJ-f3zySsT_I?^X~hHs3Qwf0aw;G1ebc}-8U$*?p&T@f=s zTgUo7fP%3-0zhUxcmp**5J@{f#V<}jCxgd6he6~$zgNOO&SermLHD77#BAix0z0KVaCJB8oXI^kS`IuY((I;*oKIs#xmI-o66Iyc%< zI^}fKx0U4?IiNCyIhDcLIgqu4xODDaxKKPKxFcq^IA`RNJO3sEJD~QCI$%+6IywSG zI`tk`I#hpKI%c9-I!c+kw?(Z$IWHNVIaKWYIqpS-xO+qoUbx*Z7r5Ea@;EL>l{>GI mjvF9z9ICghUrsr@TAw+v1_3%&F@(4vT3)!SB^0>ID77#BAix0z00OXVJC1PHI;xCJJ mjvF9KTdKFRC{8&rSD!iM_y9WgF@(4uSzfr1JQTQ@pZz$){!mK* diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/diazChl/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/diazChl/c/0/0/0/0 index 6365b670413b05a30dddc698b0f36cf8bca53488..9016835262f03e9aeaef6a1fe8f2d8a24a5c4ed9 100644 GIT binary patch delta 181 zcmV;m080PK0m%V>D77#BAix0z0JucmIW7=RIeRSlI8&s^IDzBBIOZ+SIC}NrIBc=l zI85~Gw&FsHH^>b6H!(q8IO{TbxBDDfw~0P7w-o6JH##5hIY*0MIV+F;IE8h{IH>{I zILx%*I4>&dI2zd7I6|V{wj{i-H`?PCI3p8+IEV~*x1i<-Shs^nD!27Ok%t~6vt-*i jf?d|Oa!k566lxhb>9mA6C)0PgiqcrO4{<2BWjaSUh!Il% delta 181 zcmV;m080PK0m%V>D77#BAix0z0DXtuIlV|uIXg)BI83IwYM=OBqKuZ5URA diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/diazFe/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/diazFe/c/0/0/0/0 index 3f41dae5c8ef9336a42fd018a1c5bd7638dc3d1c..4a2ff13a74907be313f79a11013d28fd63a4b54b 100644 GIT binary patch delta 177 zcmV;i08anO0m%V>D77#BAix0z07j2zG-RO7GZmh3Gis?uGqJ}aGXokrGx={$Gt-Jm zGyH(yu*MD(GPtCFGM!7)G67(Su~b)Tu~7>%u~hw~F=0HAH1a3=GxiLHGu2yGGmVl! zGa__SGk~d9GlC#cGbftIu+=y@G9_h~GV9{?GWEEJv8+D`YO&cLCb46#k%uB6663e9 f;jl(B7oM3i3-tOjioGC@gF0xDs$>F?T$@^ delta 177 zcmV;i08anO0m%V>D77#BAix0z00UZRG@UNaGxqv%GY#QJGtn+1Gwy9VGl-i`Go__T zGuMCMuz}MPGA*frGO<|HGUHl^v8G&Vv6`qdu{L$4F)BBZG{hVGGn)d1GiFd$Gs}-a zGj4BEGl8a7Gh-Z3GsD77#BAix0z0EJD&H{Y5yH`@u}HiMPDHszL*Htmb6Hq!>iHjtLV zHZHV8w1$j~HE+uGH5mIyHdE~Hv>x`zw6M~Yw59eVHG?7aH{JhRH!%VBHd!ajHUtf< zHrl(uHhlojHdC+1HU!8nv?AB8HSj?ZHYSE?HcBz?v>r?d$+QC5h_oeNk%uB6Ilv&a fVZORG!)OyW07Y&#=??F-IAY1PfewbWQB`9#3UpC( delta 177 zcmV;i08anO0m%V>D77#BAix0z0O$9_H_j3^H)2`gHk8G^HVE&MHulu3HXJd=HfOKF zHg2;-w5UUkHDpBjHF5Y!HlOA1wATH|v`lxDw76v?HNzP6H=6ZZH)ZejHbouFHmeG( zHr%(sHn{!GHjt~wHk+OK9`*X1rKp diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/spC/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/spC/c/0/0/0/0 index fdf6a46496b06a7ee361b24c7f876a333496e0ee..fd17fac136ad81500d3b730186ea736377c92503 100644 GIT binary patch literal 201 zcmV;)05<<9wJ-f3zySsT^6=q4gTm52%!=?nK?esv`=Pr&(XO*T&*qsv88f#&RaGKC z)9VF42~Z9{JH`+{J)@jHM0=h-XEwt=Hnaji590ei5Loy=I^YXGFLfV39h%fWwQzeyWE5Q{HAr3O7eC6IGJ D$~#V6 literal 201 zcmV;)05<<9wJ-f3zySsTNlW5B3pCR{+!66US5*f;4fne~aQm}9(D#`>pJ%r|j$tA{ zG<^m?8^jMkdL^Y!5N;$?p4?hDx^wRr2Ff#Z)Xr>E44KW`-9+lKS0A|WQ z^1Gcr&pyUJS~FNb2+1Zta`q)ZZ(Sch&L#FfpBV)|vsWEI-h5O)MOFbnfx`kn(TEd2 z4sIzwVY=ErnM=_=zdE2kzZ1_s%sz=fhgd^DJ6^pKl Dw%<>Z diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/spCaCO3/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/spCaCO3/c/0/0/0/0 index 8e3935c04a48dc31afd26373db0074db8748b041..452dd71fe4b9fdbe3e6dba7f562d6cd83aa7e1eb 100644 GIT binary patch literal 201 zcmV;)05<<9wJ-f3zySsTd0}2Z{cr_8%(mk`st=DoP|pNEZ~hxUtv`D|V}~n0^hCNo zyI7GvU7YAX>ywW@5``!{s^CRD=c|`IZg&_yp#C{O{*2r{%Y5EGJrRyR2@C2zapCkn zIpuag8aWR?^VIr2x-f@65(thy_fQ)?-pG$Uo+P|HUH>0FtfzxMpMg3*1}5J=iVWdC z3j&Zn8DH!^?_Tvj`tEi=A0QDwbh{Tn3*U=AvQ3CS*oG56DmK472w(&~#3w~P#?^{G Dui92! literal 201 zcmV;)05<<9wJ-f3zySsTOWs{SJ|_h~2cP3U#|e)OG{{3m-qiw-M4K~cIs zaY&OsT}$CU0b7ke!wo7tu|!Bb8`hgVB5@f$eEm5;hlSieb!px{`vZiVWdC z3j&Zn8ei-_?_Tvj`tEi=ARrMxhISV}P@s!G*FcCq{&f>Rq?5cnHX8yxigQCfsIrPa D-Pujs diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/spChl/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/spChl/c/0/0/0/0 index 3b42c03cf930d898cf750526fde13a7e6e59de75..5a982cd74d92315f664d18f06cd8b52385b25838 100644 GIT binary patch literal 201 zcmV;)05<<9wJ-f3zySsTLrASX10kI~#TKtUOO?hwxg3H$d-{Ss=IUBKGDC$uj9LUf zk-OeK80`B!Ca~{4{f}uq%BpNV2b+;SvU%A(KjSwxB;?pVWQ}n>H-v6I zFvmqbV1Io*KSw1#YftPwd9Tnter~xvNLaKzFM`88^V8)$nJ61RQ_qh*(}J5ltR%WU zaSq%)5DRHNRc2*9LE=9>KAw0zZt*uh)YbJp9ecn%VWp)#pHkL6+raTXQ=b4n8$TO9 DbU;_4 literal 201 zcmV;)05<<9wJ-f3zySsT@cyhlJR+SvXgsey$<)R@x&47XchiDB+ksj=JS&AgHe>`o zL^|I+Wo!98sB-Q-T8V2tvNmr$E!C4f4n*2LL;9CJdrzG`XZyH4gwEJKRFrW&%ad+B z>)}N`s-=BBw2~!0S(WTP!Pd_`@ff*11f#S)5%t48mK^3iXipnHQqPY)(t?{ktR%WU zaSq%)4-08MRAyy8K;l0=KAw0z%t<&tME~_Y>=nU1V_T*@SOwKRo}2GIkt6^<=E563 D6H{7U diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/spFe/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/spFe/c/0/0/0/0 index bd7ea8f933e9e07ef4a9a1a73716510e4b3fe39b..66e51dd3470acba217ce65f3b680d3069f30d346 100644 GIT binary patch literal 201 zcmV;)05<<9wJ-f3zySsTH_ux)S7Sdn7e-MxetBazAn7$XGTAIQUI7+2%^E5+ z=2jUuX?Ha@6wJXl1QV4vl#`-2^BSBt3gaO+8TmIiMYCo%U?9#nrf7OMgHmQU`Y(7l zKTnA_G+k0RlpjYoV1^qvt}Z+`hOPEDZW^sOZ-A*c1W=(j!M{T{!p~+m{6K~`rmpNa D+ki}* literal 201 zcmV;)05<<9wJ-f3zySsTVkBHQt-L=sI;l}NOXp)ZD>ybcULGws-VYWxPeCd-`rwf_ z4^MPA2>gsUdTWk1)rjUcJz)4Y$fg@N3k7~R_bzld1HN20p?YmMahHNO`wLAsK~F(9 z{Y)7*wn8;GqDR6vv2>L;V)mjp?JS%(D<2{^Bvm*!JNRZdUzyG~r)YXNgi>ZV`Y(7l zKTnA_G+k0Rl^;hpV1^qvt}Z+`qviBB1~sfV97U-&JrJQcR0KjdA`xUaH9dtl{AlYp D#dk{) diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/spP/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/spP/c/0/0/0/0 index ecf7ef73205951afd56db5cea5313a4b11581dd5..b61340f9883aca9e984c378e68fef58bfede26af 100644 GIT binary patch literal 201 zcmV;)05<<9wJ-f3zySsTbTC6ZTM0ZnxAa*%N*s7Q_h}+K+)N!iy;TJ}`)49MvCWS= za*uR7vpIk}HNb;A#Y-DIwaFbj30^xpRR(W6s83cq!>(I9?hAZ7-;a?y;Cwba@L4ZA zzTpQujw&xZUOc%wiM*0Kr(}^kSfPzO*GFtSqQ-qYF(Z*Xf(yGlI`~^VuZU+mfdYj) zdbyfAyR17qsWLS?Wcdj@QpYws#LehC4!f&6mu;6jnKhC-9#ND#dnl$mGH|UsjJMT0 Dghf(v literal 201 zcmV;)05<<9wJ-f3zySsT!E!@8rf57nQ)yZ|qB)rGC6@eZu^5fO@JIb^;RD{YV|!kC+%-L5kpoxw~Sjma-w@XMNyGErFb?w^G`23 zLD&a7((W!hZKJt6-x`xU4Z4v#S+I>ev8-%6?pc02HB*s0MsmA5I`~^Vv503ofdYj) zd%2oByR17qsxmb@Wcdj@QpYwsarEXprLC$v9XXdfv-*)ceXWu^F~*}iU)HNTZ;{kH DU-MKB diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/zarr.json b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/zarr.json index 93fb96777..b0082ed04 100644 --- a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/zarr.json +++ b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/zarr.json @@ -1,7 +1,7 @@ { "attributes": { "title": "ROMS initial conditions file created by ROMS-Tools", - "roms_tools_version": "3.1.3.dev32+g91d580e67.d20250929", + "roms_tools_version": "10000.dev355+g36bf2550d", "ini_time": "2021-06-29 00:00:00", "model_reference_date": "2000-01-01 00:00:00", "adjust_depth_for_sea_surface_height": "False", @@ -16,22 +16,16 @@ "kind": "inline", "must_understand": false, "metadata": { - "spP": { + "abs_time": { "shape": [ - 1, - 3, - 4, - 4 + 1 ], - "data_type": "float32", + "data_type": "int64", "chunk_grid": { "name": "regular", "configuration": { "chunk_shape": [ - 1, - 3, - 4, - 4 + 1 ] } }, @@ -41,7 +35,7 @@ "separator": "/" } }, - "fill_value": 0.0, + "fill_value": 0, "codecs": [ { "name": "bytes", @@ -58,22 +52,18 @@ } ], "attributes": { - "long_name": "small phytoplankton phosphorous", - "units": "mmol/m^3", - "coordinates": "abs_time", - "_FillValue": "AAAAAAAA+H8=" + "long_name": "absolute time", + "units": "days since 2021-06-29 00:00:00", + "calendar": "proleptic_gregorian" }, "dimension_names": [ - "ocean_time", - "s_rho", - "eta_rho", - "xi_rho" + "ocean_time" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "DONr": { + "ALK": { "shape": [ 1, 3, @@ -115,8 +105,8 @@ } ], "attributes": { - "long_name": "refractory dissolved organic nitrogen", - "units": "mmol/m^3", + "long_name": "alkalinity", + "units": "meq/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -130,7 +120,7 @@ "node_type": "array", "storage_transformers": [] }, - "temp": { + "ALK_ALT_CO2": { "shape": [ 1, 3, @@ -172,8 +162,8 @@ } ], "attributes": { - "long_name": "potential temperature", - "units": "degrees Celsius", + "long_name": "alkalinity, alternative CO2", + "units": "meq/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -187,22 +177,16 @@ "node_type": "array", "storage_transformers": [] }, - "PO4": { + "Cs_r": { "shape": [ - 1, - 3, - 4, - 4 + 3 ], "data_type": "float32", "chunk_grid": { "name": "regular", "configuration": { "chunk_shape": [ - 1, - 3, - 4, - 4 + 3 ] } }, @@ -229,26 +213,19 @@ } ], "attributes": { - "long_name": "dissolved inorganic phosphate", - "units": "mmol/m^3", - "coordinates": "abs_time", + "long_name": "Vertical stretching function at rho-points", + "units": "nondimensional", "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ - "ocean_time", - "s_rho", - "eta_rho", - "xi_rho" + "s_rho" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "spCaCO3": { + "Cs_w": { "shape": [ - 1, - 3, - 4, 4 ], "data_type": "float32", @@ -256,9 +233,6 @@ "name": "regular", "configuration": { "chunk_shape": [ - 1, - 3, - 4, 4 ] } @@ -286,22 +260,18 @@ } ], "attributes": { - "long_name": "small phytoplankton CaCO3", - "units": "mmol/m^3", - "coordinates": "abs_time", + "long_name": "Vertical stretching function at w-points", + "units": "nondimensional", "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ - "ocean_time", - "s_rho", - "eta_rho", - "xi_rho" + "s_w" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "spC": { + "diatC": { "shape": [ 1, 3, @@ -343,7 +313,7 @@ } ], "attributes": { - "long_name": "small phytoplankton carbon", + "long_name": "diatom carbon", "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" @@ -358,7 +328,7 @@ "node_type": "array", "storage_transformers": [] }, - "Lig": { + "diatChl": { "shape": [ 1, 3, @@ -400,8 +370,8 @@ } ], "attributes": { - "long_name": "iron binding ligand", - "units": "mmol/m^3", + "long_name": "diatom chloropyll", + "units": "mg/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -415,9 +385,10 @@ "node_type": "array", "storage_transformers": [] }, - "zeta": { + "diatFe": { "shape": [ 1, + 3, 4, 4 ], @@ -427,6 +398,7 @@ "configuration": { "chunk_shape": [ 1, + 3, 4, 4 ] @@ -455,13 +427,14 @@ } ], "attributes": { - "long_name": "sea surface height", - "units": "m", + "long_name": "diatom iron", + "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ "ocean_time", + "s_rho", "eta_rho", "xi_rho" ], @@ -469,7 +442,7 @@ "node_type": "array", "storage_transformers": [] }, - "diatFe": { + "diatP": { "shape": [ 1, 3, @@ -511,7 +484,7 @@ } ], "attributes": { - "long_name": "diatom iron", + "long_name": "diatom phosphorus", "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" @@ -526,7 +499,7 @@ "node_type": "array", "storage_transformers": [] }, - "DIC": { + "diatSi": { "shape": [ 1, 3, @@ -568,7 +541,7 @@ } ], "attributes": { - "long_name": "dissolved inorganic carbon", + "long_name": "diatom silicate", "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" @@ -583,16 +556,22 @@ "node_type": "array", "storage_transformers": [] }, - "abs_time": { + "diazC": { "shape": [ - 1 + 1, + 3, + 4, + 4 ], - "data_type": "int64", + "data_type": "float32", "chunk_grid": { "name": "regular", "configuration": { "chunk_shape": [ - 1 + 1, + 3, + 4, + 4 ] } }, @@ -602,7 +581,7 @@ "separator": "/" } }, - "fill_value": 0, + "fill_value": 0.0, "codecs": [ { "name": "bytes", @@ -619,18 +598,22 @@ } ], "attributes": { - "long_name": "absolute time", - "units": "days since 2021-06-29 00:00:00", - "calendar": "proleptic_gregorian" + "long_name": "diazotroph carbon", + "units": "mmol/m^3", + "coordinates": "abs_time", + "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ - "ocean_time" + "ocean_time", + "s_rho", + "eta_rho", + "xi_rho" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "diazP": { + "diazChl": { "shape": [ 1, 3, @@ -672,8 +655,8 @@ } ], "attributes": { - "long_name": "diazotroph phosphorus", - "units": "mmol/m^3", + "long_name": "diazotroph chloropyll", + "units": "mg/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -687,7 +670,7 @@ "node_type": "array", "storage_transformers": [] }, - "diazC": { + "diazFe": { "shape": [ 1, 3, @@ -729,7 +712,7 @@ } ], "attributes": { - "long_name": "diazotroph carbon", + "long_name": "diazotroph iron", "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" @@ -744,7 +727,7 @@ "node_type": "array", "storage_transformers": [] }, - "diatC": { + "diazP": { "shape": [ 1, 3, @@ -786,7 +769,7 @@ } ], "attributes": { - "long_name": "diatom carbon", + "long_name": "diazotroph phosphorus", "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" @@ -801,7 +784,7 @@ "node_type": "array", "storage_transformers": [] }, - "O2": { + "DIC": { "shape": [ 1, 3, @@ -843,7 +826,7 @@ } ], "attributes": { - "long_name": "dissolved oxygen", + "long_name": "dissolved inorganic carbon", "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" @@ -858,7 +841,7 @@ "node_type": "array", "storage_transformers": [] }, - "Fe": { + "DIC_ALT_CO2": { "shape": [ 1, 3, @@ -900,7 +883,7 @@ } ], "attributes": { - "long_name": "dissolved inorganic iron", + "long_name": "dissolved inorganic carbon, alternative CO2", "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" @@ -915,7 +898,7 @@ "node_type": "array", "storage_transformers": [] }, - "diatChl": { + "DOC": { "shape": [ 1, 3, @@ -957,8 +940,8 @@ } ], "attributes": { - "long_name": "diatom chloropyll", - "units": "mg/m^3", + "long_name": "dissolved organic carbon", + "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -972,7 +955,7 @@ "node_type": "array", "storage_transformers": [] }, - "diazFe": { + "DOCr": { "shape": [ 1, 3, @@ -1014,7 +997,7 @@ } ], "attributes": { - "long_name": "diazotroph iron", + "long_name": "refractory dissolved organic carbon", "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" @@ -1086,7 +1069,7 @@ "node_type": "array", "storage_transformers": [] }, - "ALK": { + "DONr": { "shape": [ 1, 3, @@ -1128,8 +1111,8 @@ } ], "attributes": { - "long_name": "alkalinity", - "units": "meq/m^3", + "long_name": "refractory dissolved organic nitrogen", + "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -1143,16 +1126,22 @@ "node_type": "array", "storage_transformers": [] }, - "ocean_time": { + "DOP": { "shape": [ - 1 + 1, + 3, + 4, + 4 ], - "data_type": "float64", + "data_type": "float32", "chunk_grid": { "name": "regular", "configuration": { "chunk_shape": [ - 1 + 1, + 3, + 4, + 4 ] } }, @@ -1179,19 +1168,26 @@ } ], "attributes": { - "long_name": "relative time: seconds since 2000-01-01 00:00:00", - "units": "seconds", + "long_name": "dissolved organic phosphorus", + "units": "mmol/m^3", + "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ - "ocean_time" + "ocean_time", + "s_rho", + "eta_rho", + "xi_rho" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "Cs_w": { + "DOPr": { "shape": [ + 1, + 3, + 4, 4 ], "data_type": "float32", @@ -1199,6 +1195,9 @@ "name": "regular", "configuration": { "chunk_shape": [ + 1, + 3, + 4, 4 ] } @@ -1226,23 +1225,27 @@ } ], "attributes": { - "long_name": "Vertical stretching function at w-points", - "units": "nondimensional", + "long_name": "refractory dissolved organic phosphorus", + "units": "mmol/m^3", + "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ - "s_w" + "ocean_time", + "s_rho", + "eta_rho", + "xi_rho" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "u": { + "Fe": { "shape": [ 1, 3, 4, - 3 + 4 ], "data_type": "float32", "chunk_grid": { @@ -1252,7 +1255,7 @@ 1, 3, 4, - 3 + 4 ] } }, @@ -1279,8 +1282,8 @@ } ], "attributes": { - "long_name": "u-flux component", - "units": "m/s", + "long_name": "dissolved inorganic iron", + "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -1288,13 +1291,13 @@ "ocean_time", "s_rho", "eta_rho", - "xi_u" + "xi_rho" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "spFe": { + "Lig": { "shape": [ 1, 3, @@ -1336,7 +1339,7 @@ } ], "attributes": { - "long_name": "small phytoplankton iron", + "long_name": "iron binding ligand", "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" @@ -1351,7 +1354,7 @@ "node_type": "array", "storage_transformers": [] }, - "DOCr": { + "NH4": { "shape": [ 1, 3, @@ -1393,7 +1396,7 @@ } ], "attributes": { - "long_name": "refractory dissolved organic carbon", + "long_name": "dissolved ammonia", "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" @@ -1408,11 +1411,12 @@ "node_type": "array", "storage_transformers": [] }, - "ubar": { + "NO3": { "shape": [ 1, + 3, 4, - 3 + 4 ], "data_type": "float32", "chunk_grid": { @@ -1420,8 +1424,9 @@ "configuration": { "chunk_shape": [ 1, + 3, 4, - 3 + 4 ] } }, @@ -1448,21 +1453,22 @@ } ], "attributes": { - "long_name": "vertically integrated u-flux component", - "units": "m/s", + "long_name": "dissolved inorganic nitrate", + "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ "ocean_time", + "s_rho", "eta_rho", - "xi_u" + "xi_rho" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "diazChl": { + "O2": { "shape": [ 1, 3, @@ -1504,8 +1510,8 @@ } ], "attributes": { - "long_name": "diazotroph chloropyll", - "units": "mg/m^3", + "long_name": "dissolved oxygen", + "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -1519,22 +1525,16 @@ "node_type": "array", "storage_transformers": [] }, - "SiO3": { + "ocean_time": { "shape": [ - 1, - 3, - 4, - 4 + 1 ], - "data_type": "float32", + "data_type": "float64", "chunk_grid": { "name": "regular", "configuration": { "chunk_shape": [ - 1, - 3, - 4, - 4 + 1 ] } }, @@ -1561,22 +1561,18 @@ } ], "attributes": { - "long_name": "dissolved inorganic silicate", - "units": "mmol/m^3", - "coordinates": "abs_time", + "long_name": "relative time: seconds since 2000-01-01 00:00:00", + "units": "seconds", "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ - "ocean_time", - "s_rho", - "eta_rho", - "xi_rho" + "ocean_time" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "salt": { + "PO4": { "shape": [ 1, 3, @@ -1618,8 +1614,8 @@ } ], "attributes": { - "long_name": "salinity", - "units": "PSU", + "long_name": "dissolved inorganic phosphate", + "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -1633,7 +1629,7 @@ "node_type": "array", "storage_transformers": [] }, - "DIC_ALT_CO2": { + "salt": { "shape": [ 1, 3, @@ -1675,8 +1671,8 @@ } ], "attributes": { - "long_name": "dissolved inorganic carbon, alternative CO2", - "units": "mmol/m^3", + "long_name": "salinity", + "units": "PSU", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -1690,10 +1686,11 @@ "node_type": "array", "storage_transformers": [] }, - "vbar": { + "SiO3": { "shape": [ 1, 3, + 4, 4 ], "data_type": "float32", @@ -1703,6 +1700,7 @@ "chunk_shape": [ 1, 3, + 4, 4 ] } @@ -1730,21 +1728,22 @@ } ], "attributes": { - "long_name": "vertically integrated v-flux component", - "units": "m/s", + "long_name": "dissolved inorganic silicate", + "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ "ocean_time", - "eta_v", + "s_rho", + "eta_rho", "xi_rho" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "DOP": { + "spC": { "shape": [ 1, 3, @@ -1786,7 +1785,7 @@ } ], "attributes": { - "long_name": "dissolved organic phosphorus", + "long_name": "small phytoplankton carbon", "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" @@ -1801,7 +1800,7 @@ "node_type": "array", "storage_transformers": [] }, - "diatP": { + "spCaCO3": { "shape": [ 1, 3, @@ -1843,7 +1842,7 @@ } ], "attributes": { - "long_name": "diatom phosphorus", + "long_name": "small phytoplankton CaCO3", "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" @@ -1858,7 +1857,7 @@ "node_type": "array", "storage_transformers": [] }, - "ALK_ALT_CO2": { + "spChl": { "shape": [ 1, 3, @@ -1900,8 +1899,8 @@ } ], "attributes": { - "long_name": "alkalinity, alternative CO2", - "units": "meq/m^3", + "long_name": "small phytoplankton chlorophyll", + "units": "mg/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -1915,7 +1914,7 @@ "node_type": "array", "storage_transformers": [] }, - "spChl": { + "spFe": { "shape": [ 1, 3, @@ -1957,8 +1956,8 @@ } ], "attributes": { - "long_name": "small phytoplankton chlorophyll", - "units": "mg/m^3", + "long_name": "small phytoplankton iron", + "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -1972,7 +1971,7 @@ "node_type": "array", "storage_transformers": [] }, - "NO3": { + "spP": { "shape": [ 1, 3, @@ -2014,7 +2013,7 @@ } ], "attributes": { - "long_name": "dissolved inorganic nitrate", + "long_name": "small phytoplankton phosphorous", "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" @@ -2029,7 +2028,7 @@ "node_type": "array", "storage_transformers": [] }, - "diatSi": { + "temp": { "shape": [ 1, 3, @@ -2071,8 +2070,8 @@ } ], "attributes": { - "long_name": "diatom silicate", - "units": "mmol/m^3", + "long_name": "potential temperature", + "units": "degrees Celsius", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -2086,12 +2085,12 @@ "node_type": "array", "storage_transformers": [] }, - "DOPr": { + "u": { "shape": [ 1, 3, 4, - 4 + 3 ], "data_type": "float32", "chunk_grid": { @@ -2101,7 +2100,7 @@ 1, 3, 4, - 4 + 3 ] } }, @@ -2128,8 +2127,8 @@ } ], "attributes": { - "long_name": "refractory dissolved organic phosphorus", - "units": "mmol/m^3", + "long_name": "u-flux component", + "units": "m/s", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, @@ -2137,18 +2136,17 @@ "ocean_time", "s_rho", "eta_rho", - "xi_rho" + "xi_u" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "DOC": { + "ubar": { "shape": [ 1, - 3, 4, - 4 + 3 ], "data_type": "float32", "chunk_grid": { @@ -2156,9 +2154,8 @@ "configuration": { "chunk_shape": [ 1, - 3, 4, - 4 + 3 ] } }, @@ -2185,26 +2182,25 @@ } ], "attributes": { - "long_name": "dissolved organic carbon", - "units": "mmol/m^3", + "long_name": "vertically integrated u-flux component", + "units": "m/s", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ "ocean_time", - "s_rho", "eta_rho", - "xi_rho" + "xi_u" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "zooC": { + "v": { "shape": [ 1, 3, - 4, + 3, 4 ], "data_type": "float32", @@ -2214,7 +2210,7 @@ "chunk_shape": [ 1, 3, - 4, + 3, 4 ] } @@ -2242,26 +2238,25 @@ } ], "attributes": { - "long_name": "zooplankton carbon", - "units": "mmol/m^3", + "long_name": "v-flux component", + "units": "m/s", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ "ocean_time", "s_rho", - "eta_rho", + "eta_v", "xi_rho" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "NH4": { + "vbar": { "shape": [ 1, 3, - 4, 4 ], "data_type": "float32", @@ -2271,7 +2266,6 @@ "chunk_shape": [ 1, 3, - 4, 4 ] } @@ -2299,15 +2293,14 @@ } ], "attributes": { - "long_name": "dissolved ammonia", - "units": "mmol/m^3", + "long_name": "vertically integrated v-flux component", + "units": "m/s", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ "ocean_time", - "s_rho", - "eta_rho", + "eta_v", "xi_rho" ], "zarr_format": 3, @@ -2371,16 +2364,20 @@ "node_type": "array", "storage_transformers": [] }, - "Cs_r": { + "zeta": { "shape": [ - 3 + 1, + 4, + 4 ], "data_type": "float32", "chunk_grid": { "name": "regular", "configuration": { "chunk_shape": [ - 3 + 1, + 4, + 4 ] } }, @@ -2407,22 +2404,25 @@ } ], "attributes": { - "long_name": "Vertical stretching function at rho-points", - "units": "nondimensional", + "long_name": "sea surface height", + "units": "m", + "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ - "s_rho" + "ocean_time", + "eta_rho", + "xi_rho" ], "zarr_format": 3, "node_type": "array", "storage_transformers": [] }, - "v": { + "zooC": { "shape": [ 1, 3, - 3, + 4, 4 ], "data_type": "float32", @@ -2432,7 +2432,7 @@ "chunk_shape": [ 1, 3, - 3, + 4, 4 ] } @@ -2460,15 +2460,15 @@ } ], "attributes": { - "long_name": "v-flux component", - "units": "m/s", + "long_name": "zooplankton carbon", + "units": "mmol/m^3", "coordinates": "abs_time", "_FillValue": "AAAAAAAA+H8=" }, "dimension_names": [ "ocean_time", "s_rho", - "eta_v", + "eta_rho", "xi_rho" ], "zarr_format": 3, diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/zooC/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_bgc_from_climatology.zarr/zooC/c/0/0/0/0 index dc1a75d871372de1cf77a5c9468fda3992556877..8c745993c7ef959b219d4c40c5bb16459c0dc1df 100644 GIT binary patch delta 177 zcmV;i08anO0m%V>D77#BAix0z0KzemKu_F;K)KFzK)sAkKxa~qK+b`lK(>*kK&yzA zK>TWzKn$0EK*Y3+K)p1JKn3!lKYJX%KZJS?K!`OU+KBtVU*k%uB6Qn;)@ f-CB)6tw)tWX*iHTxe?Ak&t3LEIpilmaukU`EoV&$ delta 177 zcmV;i08anO0m%V>D77#BAix0z04P6^K>OK-Ky1u(Ky8&yK)qCtKvQy_K&f`6K(2U| zKx=Z9K;nylKr@t$K!0V6K(P#?Kj?PAKd&(lK&PEzK=t&LKrAPQK<-_1Ku;P`K&AkY zK(`T^K>wAXKxnm*K>6OIKx#gVK=!JYK(9HFKyXFGKMFbs>OZv|BtWpbk%uB65SXk$ fCsK_-YDkqpr#O&6wZ+Xp9r*P>e#9q0&i{x&2&+!+ diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_unified_bgc_from_climatology.zarr/ALK/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_unified_bgc_from_climatology.zarr/ALK/c/0/0/0/0 index fd8ea336549f3512fa4f5937feecc78c42d89234..43a691d5311d3c29eb353bdd1225ecd52fb3220f 100644 GIT binary patch delta 24 fcmYce;?mfv|5xFFFB1dHL@q~mMg|A@O9HF_Sh@x# delta 41 xcmWGa=F-@z|5xF_RVIcf6S*8Ey7WX`6F+qdxC+DyxE^9=aIk!!CBQ6o4gfKu4W9r2 diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_unified_bgc_from_climatology.zarr/ALK_ALT_CO2/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_unified_bgc_from_climatology.zarr/ALK_ALT_CO2/c/0/0/0/0 index fd8ea336549f3512fa4f5937feecc78c42d89234..43a691d5311d3c29eb353bdd1225ecd52fb3220f 100644 GIT binary patch delta 24 fcmYce;?mfv|5xFFFB1dHL@q~mMg|A@O9HF_Sh@x# delta 41 xcmWGa=F-@z|5xF_RVIcf6S*8Ey7WX`6F+qdxC+DyxE^9=aIk!!CBQ6o4gfKu4W9r2 diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_unified_bgc_from_climatology.zarr/DIC/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_unified_bgc_from_climatology.zarr/DIC/c/0/0/0/0 index 31df7c219327c11f405a2fccbbc0523adc50d41e..3a017aceb1ff9f8be131ee7a1fb2cb9c678927bc 100644 GIT binary patch delta 24 fcmazH;?mfv|5xFFFB1dHL@q~mMg|A@O9HF_SIPz) delta 31 ncmWHF4^-kXPFopRyy!!STF+sp5O^T diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_unified_bgc_from_climatology.zarr/DIC_ALT_CO2/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_unified_bgc_from_climatology.zarr/DIC_ALT_CO2/c/0/0/0/0 index 31df7c219327c11f405a2fccbbc0523adc50d41e..3a017aceb1ff9f8be131ee7a1fb2cb9c678927bc 100644 GIT binary patch delta 24 fcmazH;?mfv|5xFFFB1dHL@q~mMg|A@O9HF_SIPz) delta 31 ncmWHF4^-kXPFopRyy!!STF+sp5O^T diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_unified_bgc_from_climatology.zarr/NO3/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_unified_bgc_from_climatology.zarr/NO3/c/0/0/0/0 index 1065da6d033061f1a96e01617326eb8e60dc168e..a3e5190f069efd18a17a8fe7c4dba38fec8c0c3b 100644 GIT binary patch delta 25 gcmXR3;?~%z|5xFFFB1a`%S3Jmc18vV`AY(<0A8mCJOBUy delta 59 zcmWG4;nvuy|5xEaEHeWO+eB^$O+Kb5`_d*q`-@_;?A0UY*iSpL)LwFrq8k6FioOxmzoHL))Iu|Xhb3Vk@<-Aw3!}(KKvva>$ OD-*-v7Y^*2GZ+DgmKG%d diff --git a/roms_tools/tests/test_setup/test_data/initial_conditions_with_unified_bgc_from_climatology.zarr/PO4/c/0/0/0/0 b/roms_tools/tests/test_setup/test_data/initial_conditions_with_unified_bgc_from_climatology.zarr/PO4/c/0/0/0/0 index 7c9d8222bfdf92dbf2f457e93402b95539349372..d40112a6bc3a8034b94d547a9a54e83cabe3497d 100644 GIT binary patch delta 25 gcmXR8;?~%z|5xFFFB1a`%S3Jmc18vV`AY(<0A7v;I{*Lx delta 58 zcmV-A0LA}Oa| Date: Fri, 15 May 2026 10:00:26 -0700 Subject: [PATCH 2/3] precommit changes --- roms_tools/setup/boundary_forcing.py | 28 +++++++------- roms_tools/setup/initial_conditions.py | 11 +++++- roms_tools/setup/utils.py | 5 +-- roms_tools/tests/test_setup/test_utils.py | 46 +++++++++++++++-------- 4 files changed, 54 insertions(+), 36 deletions(-) diff --git a/roms_tools/setup/boundary_forcing.py b/roms_tools/setup/boundary_forcing.py index d965b0775..a09a86478 100644 --- a/roms_tools/setup/boundary_forcing.py +++ b/roms_tools/setup/boundary_forcing.py @@ -22,10 +22,10 @@ from roms_tools.regrid import LateralRegridToROMS, VerticalRegrid from roms_tools.setup.utils import ( RawDataSource, + _compute_bgc_source_density, add_time_info_to_ds, check_and_set_boundaries, compute_barotropic_velocity, - _compute_bgc_source_density, compute_missing_bgc_variables, compute_potential_density, from_yaml, @@ -86,9 +86,7 @@ def _interpolate_phys_to_bgc_time( phys_for_interp = phys_da.assign_coords( {time_dim: xr.DataArray(phys_doy, dims=[time_dim])} ) - result = interpolate_cyclic_time( - phys_for_interp, time_dim, time_dim, bgc_doy - ) + result = interpolate_cyclic_time(phys_for_interp, time_dim, time_dim, bgc_doy) return result.assign_coords({time_dim: bgc_time_coord.values}) # Non-climatology: standard datetime64 linear interpolation @@ -461,7 +459,8 @@ def __post_init__(self): time_dim = bdry_data.dim_names.get("time") first_bgc = ( processed_fields.get(filtered_vars[0]) - if filtered_vars else None + if filtered_vars + else None ) if ( first_bgc is not None @@ -499,15 +498,16 @@ def __post_init__(self): # _interpolate_phys_to_bgc_time sees a real time axis. temp_sigma = self.physics_forcing.ds[f"temp_{direction}"] salt_sigma = self.physics_forcing.ds[f"salt_{direction}"] - if bgc_time_coord is not None and "abs_time" in temp_sigma.coords: - temp_sigma = ( - temp_sigma.swap_dims({"bry_time": "abs_time"}) - .rename({"abs_time": "time"}) - ) - salt_sigma = ( - salt_sigma.swap_dims({"bry_time": "abs_time"}) - .rename({"abs_time": "time"}) - ) + if ( + bgc_time_coord is not None + and "abs_time" in temp_sigma.coords + ): + temp_sigma = temp_sigma.swap_dims( + {"bry_time": "abs_time"} + ).rename({"abs_time": "time"}) + salt_sigma = salt_sigma.swap_dims( + {"bry_time": "abs_time"} + ).rename({"abs_time": "time"}) temp_sigma = _interpolate_phys_to_bgc_time( temp_sigma, "time", bgc_time_coord, bgc_climatology ) diff --git a/roms_tools/setup/initial_conditions.py b/roms_tools/setup/initial_conditions.py index 1267d3cbb..48bc033a4 100644 --- a/roms_tools/setup/initial_conditions.py +++ b/roms_tools/setup/initial_conditions.py @@ -292,7 +292,11 @@ def _process_data(self, processed_fields, type="physics"): self._get_depth_coordinates(zeta, location, "layer") # Save physics T/S at source depth levels for density-based BGC interpolation - if type == "physics" and self.use_density_interpolation and isinstance(data, LatLonDataset): + if ( + type == "physics" + and self.use_density_interpolation + and isinstance(data, LatLonDataset) + ): self._phys_depth_data = { "temp": processed_fields["temp"], "salt": processed_fields["salt"], @@ -301,7 +305,9 @@ def _process_data(self, processed_fields, type="physics"): } # Vertical regridding - processed_fields = self._regrid_vertically(data, processed_fields, var_names, type=type) + processed_fields = self._regrid_vertically( + data, processed_fields, var_names, type=type + ) # Compute barotropic velocities if "u" in var_names and "v" in var_names: @@ -461,6 +467,7 @@ def _regrid_vertically( ) if use_density: + assert self._phys_depth_data is not None source_density = _compute_bgc_source_density( self._phys_depth_data["temp"], self._phys_depth_data["salt"], diff --git a/roms_tools/setup/utils.py b/roms_tools/setup/utils.py index 837685b94..60459f761 100644 --- a/roms_tools/setup/utils.py +++ b/roms_tools/setup/utils.py @@ -582,10 +582,7 @@ def _compute_bgc_source_density( # xgcm.transform names the output dim after `target` (bgc_depth_coord), # but rename defensively in case the two arguments diverge. - if ( - phys_depth_dim != bgc_depth_dim - and phys_depth_dim in source_density.dims - ): + if phys_depth_dim != bgc_depth_dim and phys_depth_dim in source_density.dims: source_density = source_density.rename({phys_depth_dim: bgc_depth_dim}) # Add a small perturbation along the BGC depth dimension after interpolation, diff --git a/roms_tools/tests/test_setup/test_utils.py b/roms_tools/tests/test_setup/test_utils.py index 9fae1ceff..c3ad538cf 100644 --- a/roms_tools/tests/test_setup/test_utils.py +++ b/roms_tools/tests/test_setup/test_utils.py @@ -335,15 +335,22 @@ def test_compute_bgc_source_density_constant_TS(): ) phys_depth = xr.DataArray( np.array([10.0, 100.0, 500.0, 2000.0]), - dims=["depth"], coords={"depth": [10.0, 100.0, 500.0, 2000.0]}, + dims=["depth"], + coords={"depth": [10.0, 100.0, 500.0, 2000.0]}, ) bgc_depth = xr.DataArray( np.array([20.0, 300.0, 1500.0]), - dims=["depth"], coords={"depth": [20.0, 300.0, 1500.0]}, + dims=["depth"], + coords={"depth": [20.0, 300.0, 1500.0]}, ) result = _compute_bgc_source_density( - phys_temp, phys_salt, "depth", phys_depth, bgc_depth, "depth", + phys_temp, + phys_salt, + "depth", + phys_depth, + bgc_depth, + "depth", ) # gsw.sigma0 is invariant in space for constant T/S; only the monotonicity @@ -352,7 +359,10 @@ def test_compute_bgc_source_density_constant_TS(): expected_base = gsw.sigma0(S_const, T_const) expected = expected_base + np.arange(n_bgc) * 1e-7 np.testing.assert_allclose( - result.values[0, 0], expected, rtol=0.0, atol=1e-6, + result.values[0, 0], + expected, + rtol=0.0, + atol=1e-6, ) @@ -367,30 +377,34 @@ def test_compute_bgc_source_density_linear_TS_profile(): phys_S = np.array([34.5, 34.8, 34.9, 35.0]) bgc_depths = np.array([50.0, 250.0, 1200.0]) - phys_temp = xr.DataArray( - phys_T.reshape(1, 1, -1), dims=["xi", "eta", "depth"] - ) - phys_salt = xr.DataArray( - phys_S.reshape(1, 1, -1), dims=["xi", "eta", "depth"] - ) + phys_temp = xr.DataArray(phys_T.reshape(1, 1, -1), dims=["xi", "eta", "depth"]) + phys_salt = xr.DataArray(phys_S.reshape(1, 1, -1), dims=["xi", "eta", "depth"]) phys_depth = xr.DataArray( phys_depths, dims=["depth"], coords={"depth": phys_depths} ) - bgc_depth = xr.DataArray( - bgc_depths, dims=["depth"], coords={"depth": bgc_depths} - ) + bgc_depth = xr.DataArray(bgc_depths, dims=["depth"], coords={"depth": bgc_depths}) result = _compute_bgc_source_density( - phys_temp, phys_salt, "depth", phys_depth, bgc_depth, "depth", + phys_temp, + phys_salt, + "depth", + phys_depth, + bgc_depth, + "depth", ) # Manual reference: compute density at phys depths with the same perturbation # the function adds, then linearly interpolate to BGC depths, then add the # BGC-depth perturbation. - phys_density_perturbed = gsw.sigma0(phys_S, phys_T) + np.arange(len(phys_depths)) * 1e-7 + phys_density_perturbed = ( + gsw.sigma0(phys_S, phys_T) + np.arange(len(phys_depths)) * 1e-7 + ) expected = np.interp(bgc_depths, phys_depths, phys_density_perturbed) expected = expected + np.arange(len(bgc_depths)) * 1e-7 np.testing.assert_allclose( - result.values[0, 0], expected, rtol=1e-10, atol=1e-10, + result.values[0, 0], + expected, + rtol=1e-10, + atol=1e-10, ) From 785fa77b7e554c9eccb9736046218c68f842f0e4 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Tue, 19 May 2026 16:46:19 -0600 Subject: [PATCH 3/3] forgot these subsequent changes mostly to deal with unchunking in the vertical dimension so xgcm does not break there --- roms_tools/setup/boundary_forcing.py | 167 ++++++++++++---------- roms_tools/setup/initial_conditions.py | 11 +- roms_tools/setup/utils.py | 17 ++- roms_tools/tests/test_setup/test_utils.py | 12 +- roms_tools/utils.py | 4 +- 5 files changed, 118 insertions(+), 93 deletions(-) diff --git a/roms_tools/setup/boundary_forcing.py b/roms_tools/setup/boundary_forcing.py index a09a86478..9dcb9ef0e 100644 --- a/roms_tools/setup/boundary_forcing.py +++ b/roms_tools/setup/boundary_forcing.py @@ -447,83 +447,13 @@ def __post_init__(self): ) if use_density: - phys_data = self.physics_forcing._phys_bdry_depth_data[ - direction - ] - phys_temp = phys_data["temp"] - phys_salt = phys_data["salt"] - bgc_climatology = self.source["climatology"] - - # BGC tracers from any 3D BGC variable share the same - # time axis; grab it from the first one. - time_dim = bdry_data.dim_names.get("time") - first_bgc = ( - processed_fields.get(filtered_vars[0]) - if filtered_vars - else None - ) - if ( - first_bgc is not None - and time_dim - and time_dim in first_bgc.dims - ): - bgc_time_coord = first_bgc[time_dim] - else: - bgc_time_coord = None - - if bgc_time_coord is not None and "time" in phys_temp.dims: - phys_temp = _interpolate_phys_to_bgc_time( - phys_temp, "time", bgc_time_coord, bgc_climatology - ) - phys_salt = _interpolate_phys_to_bgc_time( - phys_salt, "time", bgc_time_coord, bgc_climatology + source_density, target_density = ( + self._compute_bgc_density_coords( + direction, + bdry_data, + processed_fields, + filtered_vars, ) - elif "time" in phys_temp.dims: - # No BGC time found — collapse to mean as fallback - phys_temp = phys_temp.mean("time") - phys_salt = phys_salt.mean("time") - - source_density = _compute_bgc_source_density( - phys_temp, - phys_salt, - phys_data["depth_dim"], - phys_data["depth_coord"], - bdry_data.ds[bdry_data.dim_names["depth"]], - bdry_data.dim_names["depth"], - ) - - # Target density: physics sigma T/S interpolated to BGC time. - # Physics BC stores both bry_time (float days since reference) - # and abs_time (datetime64). We swap to the datetime64 view so - # _interpolate_phys_to_bgc_time sees a real time axis. - temp_sigma = self.physics_forcing.ds[f"temp_{direction}"] - salt_sigma = self.physics_forcing.ds[f"salt_{direction}"] - if ( - bgc_time_coord is not None - and "abs_time" in temp_sigma.coords - ): - temp_sigma = temp_sigma.swap_dims( - {"bry_time": "abs_time"} - ).rename({"abs_time": "time"}) - salt_sigma = salt_sigma.swap_dims( - {"bry_time": "abs_time"} - ).rename({"abs_time": "time"}) - temp_sigma = _interpolate_phys_to_bgc_time( - temp_sigma, "time", bgc_time_coord, bgc_climatology - ) - salt_sigma = _interpolate_phys_to_bgc_time( - salt_sigma, "time", bgc_time_coord, bgc_climatology - ) - elif "bry_time" in temp_sigma.dims: - temp_sigma = temp_sigma.mean("bry_time") - salt_sigma = salt_sigma.mean("bry_time") - - target_density = compute_potential_density( - temp_sigma, salt_sigma - ) - n_s = target_density.sizes["s_rho"] - target_density = target_density + xr.DataArray( - np.arange(n_s) * 1e-7, dims=["s_rho"] ) for var_name in filtered_vars: @@ -583,6 +513,91 @@ def __post_init__(self): self.ds = ds + def _compute_bgc_density_coords( + self, + direction: str, + bdry_data, + processed_fields: dict, + filtered_vars: list, + ) -> tuple[xr.DataArray, xr.DataArray]: + """Build source/target density coordinates for BGC density-space interpolation. + + Source density is computed from the physics BC's saved depth-level T/S and + regridded onto the BGC depth axis. Target density is computed from the + physics BC's sigma-level T/S, interpolated to the BGC time coordinate. + + Both arrays include the small monotonicity perturbation that + ``vertical_regrid.apply`` requires for a strictly monotonic coordinate. + """ + assert self.physics_forcing is not None + phys_data = self.physics_forcing._phys_bdry_depth_data[direction] + phys_temp, phys_salt = phys_data["temp"], phys_data["salt"] + bgc_climatology = bool(self.source["climatology"]) + + # The BGC time axis (if any) is shared across all 3D BGC vars; take it from + # the first one. dim_names["time"] gives the source-data convention name. + bgc_time_dim = bdry_data.dim_names.get("time") + first_bgc = processed_fields.get(filtered_vars[0]) if filtered_vars else None + bgc_time_coord = ( + first_bgc[bgc_time_dim] + if first_bgc is not None + and bgc_time_dim in (first_bgc.dims if first_bgc is not None else ()) + else None + ) + + def _align_time(da: xr.DataArray, time_dim: str) -> xr.DataArray: + """Align ``da``'s ``time_dim`` to the BGC time axis, or collapse it.""" + if time_dim not in da.dims: + return da + if bgc_time_coord is not None: + return _interpolate_phys_to_bgc_time( + da, time_dim, bgc_time_coord, bgc_climatology + ) + return da.mean(time_dim) + + # --- Source density: physics depth-level T/S → BGC depth axis --- + phys_temp = _align_time(phys_temp, "time") + phys_salt = _align_time(phys_salt, "time") + bgc_depth_dim = bdry_data.dim_names["depth"] + source_density = _compute_bgc_source_density( + phys_temp, + phys_salt, + phys_data["depth_dim"], + phys_data["depth_coord"], + bdry_data.ds[bgc_depth_dim], + bgc_depth_dim, + ) + + # --- Target density: physics sigma-level T/S, aligned to BGC time --- + # Physics BC dataset uses "bry_time" as the time dim with an "abs_time" + # datetime64 companion coord. Swap to the datetime view before time-aligning. + temp_sigma = self.physics_forcing.ds[f"temp_{direction}"] + salt_sigma = self.physics_forcing.ds[f"salt_{direction}"] + if "abs_time" in temp_sigma.coords: + temp_sigma = temp_sigma.swap_dims({"bry_time": "abs_time"}).rename( + {"abs_time": "time"} + ) + salt_sigma = salt_sigma.swap_dims({"bry_time": "abs_time"}).rename( + {"abs_time": "time"} + ) + temp_sigma = _align_time(temp_sigma, "time") + salt_sigma = _align_time(salt_sigma, "time") + else: + temp_sigma = _align_time(temp_sigma, "bry_time") + salt_sigma = _align_time(salt_sigma, "bry_time") + + target_density = compute_potential_density(temp_sigma, salt_sigma) + # Add a small perturbation along the vertical dim to ensure monotonicity. + vertical_dim = next(d for d in target_density.dims if d.startswith("s_")) + n_s = target_density.sizes[vertical_dim] + target_density = target_density + xr.DataArray( + np.arange(n_s) * 1e-7, dims=[vertical_dim] + ) + # xgcm.transform requires a single chunk along the target vertical dim. + target_density = target_density.chunk({vertical_dim: -1}) + + return source_density, target_density + def _input_checks(self) -> None: """Validate and normalize user-provided input parameters.""" # ------------------------------------------------------- diff --git a/roms_tools/setup/initial_conditions.py b/roms_tools/setup/initial_conditions.py index 48bc033a4..0efc1a111 100644 --- a/roms_tools/setup/initial_conditions.py +++ b/roms_tools/setup/initial_conditions.py @@ -479,11 +479,16 @@ def _regrid_vertically( target_density = compute_potential_density( processed_fields["temp"], processed_fields["salt"] ) - # Add small perturbation to target density to ensure monotonicity - n_s = target_density.sizes["s_rho"] + # Add a small perturbation along the vertical dim to ensure monotonicity. + vertical_dim = next( + d for d in target_density.dims if d.startswith("s_") + ) + n_s = target_density.sizes[vertical_dim] target_density = target_density + xr.DataArray( - np.arange(n_s) * 1e-7, dims=["s_rho"] + np.arange(n_s) * 1e-7, dims=[vertical_dim] ) + # xgcm.transform requires a single chunk along the target vertical dim. + target_density = target_density.chunk({vertical_dim: -1}) for var_name in filtered_vars: if var_name not in processed_fields: diff --git a/roms_tools/setup/utils.py b/roms_tools/setup/utils.py index 60459f761..96ddfe110 100644 --- a/roms_tools/setup/utils.py +++ b/roms_tools/setup/utils.py @@ -20,6 +20,7 @@ from roms_tools.constants import R_EARTH from roms_tools.regrid import VerticalRegrid +from roms_tools.utils import transpose_dimensions if typing.TYPE_CHECKING: from roms_tools.setup.grid import Grid @@ -521,13 +522,18 @@ def compute_potential_density( xr.DataArray Potential density anomaly sigma-0 (kg/m³ - 1000). """ - return xr.apply_ufunc( + density = xr.apply_ufunc( gsw.sigma0, salt, temp, dask="parallelized", output_dtypes=[temp.dtype], ) + density = transpose_dimensions(density) + density.name = "sigma0" + density.attrs["long_name"] = "potential density anomaly" + density.attrs["units"] = "kg/m^3 - 1000" + return density def _compute_bgc_source_density( @@ -571,6 +577,9 @@ def _compute_bgc_source_density( perturbation = xr.DataArray(np.arange(n_depth) * 1e-7, dims=[phys_depth_dim]) density = density + perturbation + # xgcm.transform requires a single chunk along the dim being transformed. + density = density.chunk({phys_depth_dim: -1}) + # Regrid density from physics depth levels to BGC depth levels ds_phys = xr.Dataset({phys_depth_dim: phys_depth_coord}) vertical_regrid = VerticalRegrid(ds_phys, source_dim=phys_depth_dim) @@ -579,11 +588,7 @@ def _compute_bgc_source_density( source_depth_coords=phys_depth_coord, target_depth_coords=bgc_depth_coord, ) - - # xgcm.transform names the output dim after `target` (bgc_depth_coord), - # but rename defensively in case the two arguments diverge. - if phys_depth_dim != bgc_depth_dim and phys_depth_dim in source_density.dims: - source_density = source_density.rename({phys_depth_dim: bgc_depth_dim}) + source_density = transpose_dimensions(source_density) # Add a small perturbation along the BGC depth dimension after interpolation, # so the density profile xgcm uses as a source coordinate is strictly diff --git a/roms_tools/tests/test_setup/test_utils.py b/roms_tools/tests/test_setup/test_utils.py index c3ad538cf..66aea609e 100644 --- a/roms_tools/tests/test_setup/test_utils.py +++ b/roms_tools/tests/test_setup/test_utils.py @@ -328,10 +328,10 @@ def test_compute_bgc_source_density_constant_TS(): n_bgc = 3 T_const, S_const = 10.0, 35.0 phys_temp = xr.DataArray( - np.full((1, 1, n_phys), T_const), dims=["xi", "eta", "depth"] + np.full((n_phys, 1, 1), T_const), dims=["depth", "eta", "xi"] ) phys_salt = xr.DataArray( - np.full((1, 1, n_phys), S_const), dims=["xi", "eta", "depth"] + np.full((n_phys, 1, 1), S_const), dims=["depth", "eta", "xi"] ) phys_depth = xr.DataArray( np.array([10.0, 100.0, 500.0, 2000.0]), @@ -359,7 +359,7 @@ def test_compute_bgc_source_density_constant_TS(): expected_base = gsw.sigma0(S_const, T_const) expected = expected_base + np.arange(n_bgc) * 1e-7 np.testing.assert_allclose( - result.values[0, 0], + result.values[:, 0, 0], expected, rtol=0.0, atol=1e-6, @@ -377,8 +377,8 @@ def test_compute_bgc_source_density_linear_TS_profile(): phys_S = np.array([34.5, 34.8, 34.9, 35.0]) bgc_depths = np.array([50.0, 250.0, 1200.0]) - phys_temp = xr.DataArray(phys_T.reshape(1, 1, -1), dims=["xi", "eta", "depth"]) - phys_salt = xr.DataArray(phys_S.reshape(1, 1, -1), dims=["xi", "eta", "depth"]) + phys_temp = xr.DataArray(phys_T.reshape(-1, 1, 1), dims=["depth", "eta", "xi"]) + phys_salt = xr.DataArray(phys_S.reshape(-1, 1, 1), dims=["depth", "eta", "xi"]) phys_depth = xr.DataArray( phys_depths, dims=["depth"], coords={"depth": phys_depths} ) @@ -403,7 +403,7 @@ def test_compute_bgc_source_density_linear_TS_profile(): expected = expected + np.arange(len(bgc_depths)) * 1e-7 np.testing.assert_allclose( - result.values[0, 0], + result.values[:, 0, 0], expected, rtol=1e-10, atol=1e-10, diff --git a/roms_tools/utils.py b/roms_tools/utils.py index 266e3aefe..1c6455c8e 100644 --- a/roms_tools/utils.py +++ b/roms_tools/utils.py @@ -819,11 +819,11 @@ def transpose_dimensions(da: xr.DataArray) -> xr.DataArray: Returns ------- xarray.DataArray - The DataArray with dimensions reordered so that 'time', 's_*', 'eta_*', + The DataArray with dimensions reordered so that 'time', 's_*', 'depth', 'eta_*', and 'xi_*' are first, in that order, if they exist. """ # List of preferred dimension patterns - preferred_order = ["time", "s_", "eta_", "xi_"] + preferred_order = ["time", "s_", "depth", "eta_", "xi_"] # Get the existing dimensions in the DataArray dims = list(da.dims)