From 9d93fd17d9f17a7ce7128a4973ee64b6ca18ee0d Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 3 Sep 2026 09:11:51 -0700 Subject: [PATCH 01/17] Add ISMIP7 GrIS OCX (reanalysis) forcing support Handle the OCX scenario in the ismip7_forcing test group. OCX has no distinct ESM model: it uses RACMO2.3p2-ERA (atmosphere) and EN4 (ocean) at data version v1, with a named grid resolution (ocean-1000m) in the ocean file names. get_params() now applies OCX-specific overrides (sources, versions, ocean grid token) when scenario == OCX, so the [ismip7] model option is ignored and a single config file drives both the atmosphere and ocean_thermal test cases. Adds ismip7_forcing_ocx_gis.cfg and updates the user and developer docs. --- .../atmosphere/process_runoff.py | 2 + .../ismip7_forcing/atmosphere/process_smb.py | 2 + .../atmosphere/process_smb_gradient.py | 2 + .../atmosphere/process_temperature.py | 2 + .../process_temperature_gradient.py | 2 + .../tests/ismip7_forcing/ice_sheet_params.py | 32 +++++++- .../ismip7_forcing/ismip7_forcing_ocx_gis.cfg | 77 +++++++++++++++++++ .../ocean_thermal/process_thermal_forcing.py | 5 +- .../landice/test_groups/ismip7_forcing.rst | 14 +++- .../landice/test_groups/ismip7_forcing.rst | 22 ++++++ 10 files changed, 156 insertions(+), 4 deletions(-) create mode 100644 compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py index e03e85ecd6..ef7ac82164 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -72,6 +72,8 @@ def run(self): prefix = params['prefix'] resolution = params['atm_resolution'] version = params['atm_version'] + if params['atm_model'] is not None: + model = params['atm_model'] input_path = os.path.join(base_path_ismip7, "mrro", version) file_pattern = (f"mrro_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py index 222fb6e92e..014e098bc9 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -69,6 +69,8 @@ def run(self): prefix = params['prefix'] resolution = params['atm_resolution'] version = params['atm_version'] + if params['atm_model'] is not None: + model = params['atm_model'] input_path = os.path.join(base_path_ismip7, "acabf", version) file_pattern = (f"acabf_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py index 2f3e4c2cb4..6d39bb304e 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -70,6 +70,8 @@ def run(self): prefix = params['prefix'] resolution = params['atm_resolution'] version = params['atm_version'] + if params['atm_model'] is not None: + model = params['atm_model'] input_path = os.path.join(base_path_ismip7, "dacabfdz", version) file_pattern = (f"dacabfdz_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index caa8f6be9f..13ec7c15d9 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -69,6 +69,8 @@ def run(self): prefix = params['prefix'] resolution = params['atm_resolution'] version = params['atm_version'] + if params['atm_model'] is not None: + model = params['atm_model'] input_path = os.path.join(base_path_ismip7, "ts", version) file_pattern = (f"ts_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py index 5dec7b78bc..7399eda925 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -71,6 +71,8 @@ def run(self): prefix = params['prefix'] resolution = params['atm_resolution'] version = params['atm_version'] + if params['atm_model'] is not None: + model = params['atm_model'] input_path = os.path.join(base_path_ismip7, "dtsdz", version) file_pattern = (f"dtsdz_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") diff --git a/compass/landice/tests/ismip7_forcing/ice_sheet_params.py b/compass/landice/tests/ismip7_forcing/ice_sheet_params.py index 8324236c56..0d57094e13 100644 --- a/compass/landice/tests/ismip7_forcing/ice_sheet_params.py +++ b/compass/landice/tests/ismip7_forcing/ice_sheet_params.py @@ -10,8 +10,11 @@ 'atm_resolution': '2000m', 'atm_version': 'v2', 'ocean_version': 'v3', + 'ocean_grid': 'ocean', 'ocean_3d': True, 'ocean_temporal': 'decade', + 'atm_model': None, + 'ocean_model': None, }, 'gis': { 'projection': 'gis-bamber', @@ -19,8 +22,25 @@ 'atm_resolution': '1000m', 'atm_version': 'v2', 'ocean_version': 'v2', + 'ocean_grid': 'ocean', 'ocean_3d': False, 'ocean_temporal': 'yearly', + 'atm_model': None, + 'ocean_model': None, + }, +} + +# Overrides applied for the OCX (reanalysis) scenario. OCX has no distinct +# ESM model: it uses fixed reanalysis products (RACMO for the atmosphere and +# EN4 for the ocean) at data version v1. When scenario is 'OCX' the [ismip7] +# model option is ignored and these sources are used instead. +_OCX_OVERRIDES = { + 'gis': { + 'atm_version': 'v1', + 'ocean_version': 'v1', + 'ocean_grid': 'ocean-1000m', + 'atm_model': 'RACMO2.3p2-ERA', + 'ocean_model': 'EN4', }, } @@ -44,4 +64,14 @@ def get_params(config): raise ValueError( f"Unknown ice_sheet '{ice_sheet}'. " f"Must be one of: {list(_PARAMS.keys())}") - return _PARAMS[ice_sheet] + params = dict(_PARAMS[ice_sheet]) + + scenario = config.get("ismip7", "scenario") + if scenario == "OCX": + if ice_sheet not in _OCX_OVERRIDES: + raise ValueError( + f"The OCX scenario is not yet supported for ice_sheet " + f"'{ice_sheet}'.") + params.update(_OCX_OVERRIDES[ice_sheet]) + + return params diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg new file mode 100644 index 0000000000..d3658027db --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg @@ -0,0 +1,77 @@ +# Example config for processing the GrIS ISMIP7 OCX (reanalysis) forcing. +# OCX has no distinct ESM model: it uses RACMO2.3p2-ERA for the atmosphere and +# EN4 for the ocean, selected automatically when scenario = OCX. The [ismip7] +# model option is ignored for OCX (set it to None). The same config file drives +# both test cases; set them up individually, e.g.: +# compass setup -t landice/ismip7_forcing/atmosphere -w WORKDIR -f ismip7_forcing_ocx_gis.cfg +# compass setup -t landice/ismip7_forcing/ocean_thermal -w WORKDIR -f ismip7_forcing_ocx_gis.cfg + +# config options for ismip7 forcing data +[ismip7] + +# Ice sheet: ais (Antarctic) or gis (Greenland) +ice_sheet = gis + +# Base path to the input ISMIP7 forcing files. User has to supply. +base_path_ismip7 = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/forcing/GIS/OCX + +# Base path to the MALI mesh. User has to supply. +base_path_mali = /global/cfs/cdirs/fanssie/MALI_input_files/GIS_1to10km_r02/ + +# Base path to which output forcing files are saved. +output_base_path = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/test_processing/GIS + +# OCX has no distinct model (RACMO for atmosphere, EN4 for ocean are selected +# automatically), so this option is ignored. +model = None + +# Scenario for forcing data. +scenario = OCX + +# Name of the MALI mesh. Used to name mapping and output files. +mali_mesh_name = GIS_1to10km_r02_20230202 + +# MALI mesh file. User has to supply. +mali_mesh_file = GIS_1to10km_r02_20230202.nc + +# Number of MPI tasks for ESMF_RegridWeightGen +esmf_ntasks = 512 + +# Whether to process time-varying ocean thermal forcing (ESM scenario data) +process_ocean_thermal = true + +# Whether to process observational ocean thermal forcing climatology +process_ocean_climatology = false + +# config options for ismip7 atmosphere forcing +[ismip7_atmosphere] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = conserve + +# Start year for processing +start_year = 1990 + +# End year for processing +end_year = 2025 + +# config options for ismip7 ocean thermal forcing +[ismip7_ocean_thermal] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = bilinear + +# Start year for processing +start_year = 1990 + +# End year for processing +end_year = 2025 + +# config options for ismip7 ocean thermal forcing climatology +[ismip7_ocean_climatology] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = bilinear + +# Base path to observational climatology data +base_path_climatology = None diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py index a1469c8f2d..3e9aea9c87 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -91,11 +91,14 @@ def _run_scenario(self): # Discover input files prefix = params['prefix'] ocean_version = params['ocean_version'] + ocean_grid = params['ocean_grid'] ocean_3d = params['ocean_3d'] + if params['ocean_model'] is not None: + model = params['ocean_model'] input_path = os.path.join(base_path_ismip7, "ocean", "tf", ocean_version) file_pattern = (f"tf_{prefix}_{model}_{scenario}_" - f"ocean_{ocean_version}_*.nc") + f"{ocean_grid}_{ocean_version}_*.nc") all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) if not all_files: diff --git a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst index b535bf081a..0a37029290 100644 --- a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst @@ -24,11 +24,19 @@ ice_sheet_params The module :py:mod:`compass.landice.tests.ismip7_forcing.ice_sheet_params` defines a dictionary of ice-sheet-specific parameters (projection, file -naming prefix, grid resolution, data version, ocean dimensionality) and -provides the function +naming prefix, grid resolution, data version, ocean dimensionality, and the +atmosphere/ocean source names) and provides the function :py:func:`compass.landice.tests.ismip7_forcing.ice_sheet_params.get_params` to retrieve them based on the ``ice_sheet`` config option. +When ``scenario = OCX``, ``get_params`` applies a set of OCX overrides on top +of the ice-sheet defaults: data version ``v1``, the ocean file-name grid token +(e.g. ``ocean-1000m``), and the fixed reanalysis sources (``atm_model`` = +``RACMO2.3p2-ERA`` and ``ocean_model`` = ``EN4``). The processing steps use +``atm_model`` / ``ocean_model`` in place of the ``[ismip7] model`` option when +they are set, so the OCX ``model`` option is ignored. This keeps OCX handling +centralized and lets a single config file drive both test cases. + configure ~~~~~~~~~ @@ -39,6 +47,8 @@ the user (i.e., are not ``NotAvailable``). Repository-local example user configs are available at ``compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg`` (AIS) and ``compass/landice/tests/ismip7_forcing/ismip7_forcing_test_gis.cfg`` (GrIS). +The GrIS OCX scenario has a dedicated example +``ismip7_forcing_ocx_gis.cfg`` in the same directory. These are intended for development/testing and include environment-specific paths. diff --git a/docs/users_guide/landice/test_groups/ismip7_forcing.rst b/docs/users_guide/landice/test_groups/ismip7_forcing.rst index fdec1eee14..0f9ec60eba 100644 --- a/docs/users_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/users_guide/landice/test_groups/ismip7_forcing.rst @@ -69,6 +69,14 @@ The AIS example enables both ocean processing modes and uses a 2015-2300 processing window for both atmosphere and ocean thermal scenario forcing. The GrIS example enables scenario ocean processing only and uses 1980-2015. +For the GrIS OCX (reanalysis) scenario, use the dedicated example config +``compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg``. OCX has +no distinct ESM model: it uses ``RACMO2.3p2-ERA`` for the atmosphere and +``EN4`` for the ocean, both selected automatically when ``scenario = OCX``. +The ``[ismip7] model`` option is ignored for OCX (set it to ``None``), so a +single config file processes both the ``atmosphere`` and ``ocean_thermal`` +test cases, just like the ESM scenarios. + .. _landice_ismip7_forcing_input_data: Input Data @@ -121,6 +129,20 @@ For GrIS ocean thermal (same 1km grid, 2D, yearly files): ocean/tf/v2/tf_GrIS_{model}_{scenario}_ocean_v2_{year}.nc +The OCX (reanalysis) scenario follows the same directory layout as the ESM +scenarios, but uses fixed reanalysis sources, data version ``v1``, and a +named grid resolution in the ocean file names. For GrIS OCX the atmosphere +source is ``RACMO2.3p2-ERA`` and the ocean source is ``EN4``: + +.. code-block:: none + + acabf/v1/acabf_GrIS_RACMO2.3p2-ERA_OCX_SDBN1-1000m_v1_{year}.nc + ocean/tf/v1/tf_GrIS_EN4_OCX_ocean-1000m_v1_{year}.nc + +Set ``base_path_ismip7`` to the ``OCX`` directory and ``scenario = OCX``. +The sources, version, and ocean grid token are selected automatically for +OCX, and the ``[ismip7] model`` option is ignored (set it to ``None``). + .. _landice_ismip7_forcing_config: config options From 4c7116c215006206450d0a1521f4215caf9b8871 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 3 Sep 2026 10:39:19 -0700 Subject: [PATCH 02/17] Fix ISMIP7 forcing atmosphere/ocean processing bugs Two fixes that block processing of forcing whose remapped files or source directories differ from the CESM test data: - Drop the ncremap-mangled 'time_bnds' variable when combining remapped files. ncremap collides the bounds 'nv' dimension with the destination mesh corner dimension and pads with a fill value that overflows time decoding. The steps build their own xtime and never use time_bnds. - Skip input files whose trailing token is not a year (e.g. climatology averages like *_1978-2007_avg.nc) instead of crashing on int('avg'). --- .../tests/ismip7_forcing/atmosphere/process_runoff.py | 9 +++++++-- .../tests/ismip7_forcing/atmosphere/process_smb.py | 11 ++++++++--- .../ismip7_forcing/atmosphere/process_smb_gradient.py | 9 +++++++-- .../ismip7_forcing/atmosphere/process_temperature.py | 9 +++++++-- .../atmosphere/process_temperature_gradient.py | 9 +++++++-- .../ocean_thermal/process_thermal_forcing.py | 6 ++++-- 6 files changed, 40 insertions(+), 13 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py index ef7ac82164..06852c6690 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -87,7 +87,11 @@ def run(self): # Filter to requested year range input_files = [] for f in all_files: - year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + # skip non-yearly files such as climatology averages (*_avg.nc) + token = os.path.basename(f).split("_")[-1].replace(".nc", "") + if not token.isdigit(): + continue + year = int(token) if start_year <= year <= end_year: input_files.append(f) @@ -178,7 +182,8 @@ def _combine_and_rename(self, remapped_files, output_file): Output file path """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Rename dimensions to MALI conventions rename_dims = {} diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py index 014e098bc9..57444cbb09 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -84,8 +84,12 @@ def run(self): # Filter to requested year range input_files = [] for f in all_files: - # Extract year from filename (last part before .nc) - year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + # Extract year from filename (last part before .nc); skip + # non-yearly files such as climatology averages (e.g. *_avg.nc) + token = os.path.basename(f).split("_")[-1].replace(".nc", "") + if not token.isdigit(): + continue + year = int(token) if start_year <= year <= end_year: input_files.append(f) @@ -166,7 +170,8 @@ def _combine_and_rename(self, remapped_files, output_file): Output file path """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Rename dimensions to MALI conventions rename_dims = {} diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py index 6d39bb304e..e08d11713a 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -85,7 +85,11 @@ def run(self): # Filter to requested year range input_files = [] for f in all_files: - year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + # skip non-yearly files such as climatology averages (*_avg.nc) + token = os.path.basename(f).split("_")[-1].replace(".nc", "") + if not token.isdigit(): + continue + year = int(token) if start_year <= year <= end_year: input_files.append(f) @@ -167,7 +171,8 @@ def _combine_and_rename(self, remapped_files, output_file): Output file path """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Rename dimensions to MALI conventions rename_dims = {} diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index 13ec7c15d9..9e8f0490ba 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -84,7 +84,11 @@ def run(self): # Filter to requested year range input_files = [] for f in all_files: - year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + # skip non-yearly files such as climatology averages (*_avg.nc) + token = os.path.basename(f).split("_")[-1].replace(".nc", "") + if not token.isdigit(): + continue + year = int(token) if start_year <= year <= end_year: input_files.append(f) @@ -166,7 +170,8 @@ def _combine_and_rename(self, remapped_files, output_file): Output file path """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Rename dimensions to MALI conventions rename_dims = {} diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py index 7399eda925..9968101810 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -86,7 +86,11 @@ def run(self): # Filter to requested year range input_files = [] for f in all_files: - year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + # skip non-yearly files such as climatology averages (*_avg.nc) + token = os.path.basename(f).split("_")[-1].replace(".nc", "") + if not token.isdigit(): + continue + year = int(token) if start_year <= year <= end_year: input_files.append(f) @@ -168,7 +172,8 @@ def _combine_and_rename(self, remapped_files, output_file): Output file path """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Rename dimensions to MALI conventions rename_dims = {} diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py index 3e9aea9c87..9a67071bc9 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -308,7 +308,8 @@ def _combine_and_rename_3d(self, remapped_files, output_file, Last year to include in output """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Subset to requested year range years = ds.time.dt.year @@ -422,7 +423,8 @@ def _combine_and_rename_2d(self, remapped_files, output_file, Last year to include in output """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Subset to requested year range years = ds.time.dt.year From 198473b6a6acea98830982c7d030d4618679079a Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 3 Sep 2026 10:43:02 -0700 Subject: [PATCH 03/17] Write ISMIP7 forcing output in the layout ismip7_run ingests Restructure processed-forcing output from {output_base_path}/atmosphere_forcing/{model}_{scenario}/ (and the analogous ocean_thermal_forcing path) to {output_base_path}/{group}/atmosphere/ and {output_base_path}/{group}/ocean_thermal_forcing/, matching the directory layout the ismip7_run test group globs for. 'group' is {model}_{scenario} for the ESM scenarios and {scenario} (i.e. OCX) for OCX, whose atmosphere and ocean sources differ but must share one directory. Documents the layout and how to point ismip7_run's forcing_basepath / ocx_forcing_path at it. --- .../atmosphere/process_runoff.py | 7 ++++-- .../ismip7_forcing/atmosphere/process_smb.py | 7 ++++-- .../atmosphere/process_smb_gradient.py | 7 ++++-- .../atmosphere/process_temperature.py | 7 ++++-- .../process_temperature_gradient.py | 7 ++++-- .../ocean_thermal/process_thermal_forcing.py | 7 ++++-- .../landice/test_groups/ismip7_forcing.rst | 25 +++++++++++++++++++ 7 files changed, 55 insertions(+), 12 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py index 06852c6690..3f0e7b9e5f 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -73,7 +73,10 @@ def run(self): resolution = params['atm_resolution'] version = params['atm_version'] if params['atm_model'] is not None: + forcing_group = scenario model = params['atm_model'] + else: + forcing_group = f"{model}_{scenario}" input_path = os.path.join(base_path_ismip7, "mrro", version) file_pattern = (f"mrro_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") @@ -158,8 +161,8 @@ def run(self): os.remove(f) # Place output in appropriate directory - output_path = os.path.join(output_base_path, "atmosphere_forcing", - f"{model}_{scenario}") + output_path = os.path.join(output_base_path, forcing_group, + "atmosphere") if not os.path.exists(output_path): os.makedirs(output_path) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py index 57444cbb09..04de34180f 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -70,7 +70,10 @@ def run(self): resolution = params['atm_resolution'] version = params['atm_version'] if params['atm_model'] is not None: + forcing_group = scenario model = params['atm_model'] + else: + forcing_group = f"{model}_{scenario}" input_path = os.path.join(base_path_ismip7, "acabf", version) file_pattern = (f"acabf_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") @@ -146,8 +149,8 @@ def run(self): os.remove(f) # Place output in appropriate directory - output_path = os.path.join(output_base_path, "atmosphere_forcing", - f"{model}_{scenario}") + output_path = os.path.join(output_base_path, forcing_group, + "atmosphere") if not os.path.exists(output_path): os.makedirs(output_path) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py index e08d11713a..d6c291e4d5 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -71,7 +71,10 @@ def run(self): resolution = params['atm_resolution'] version = params['atm_version'] if params['atm_model'] is not None: + forcing_group = scenario model = params['atm_model'] + else: + forcing_group = f"{model}_{scenario}" input_path = os.path.join(base_path_ismip7, "dacabfdz", version) file_pattern = (f"dacabfdz_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") @@ -147,8 +150,8 @@ def run(self): os.remove(f) # Place output in appropriate directory - output_path = os.path.join(output_base_path, "atmosphere_forcing", - f"{model}_{scenario}") + output_path = os.path.join(output_base_path, forcing_group, + "atmosphere") if not os.path.exists(output_path): os.makedirs(output_path) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index 9e8f0490ba..3168c69f28 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -70,7 +70,10 @@ def run(self): resolution = params['atm_resolution'] version = params['atm_version'] if params['atm_model'] is not None: + forcing_group = scenario model = params['atm_model'] + else: + forcing_group = f"{model}_{scenario}" input_path = os.path.join(base_path_ismip7, "ts", version) file_pattern = (f"ts_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") @@ -146,8 +149,8 @@ def run(self): os.remove(f) # Place output in appropriate directory - output_path = os.path.join(output_base_path, "atmosphere_forcing", - f"{model}_{scenario}") + output_path = os.path.join(output_base_path, forcing_group, + "atmosphere") if not os.path.exists(output_path): os.makedirs(output_path) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py index 9968101810..16181f02e5 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -72,7 +72,10 @@ def run(self): resolution = params['atm_resolution'] version = params['atm_version'] if params['atm_model'] is not None: + forcing_group = scenario model = params['atm_model'] + else: + forcing_group = f"{model}_{scenario}" input_path = os.path.join(base_path_ismip7, "dtsdz", version) file_pattern = (f"dtsdz_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") @@ -148,8 +151,8 @@ def run(self): os.remove(f) # Place output in appropriate directory - output_path = os.path.join(output_base_path, "atmosphere_forcing", - f"{model}_{scenario}") + output_path = os.path.join(output_base_path, forcing_group, + "atmosphere") if not os.path.exists(output_path): os.makedirs(output_path) diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py index 9a67071bc9..fba0053b1d 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -94,7 +94,10 @@ def _run_scenario(self): ocean_grid = params['ocean_grid'] ocean_3d = params['ocean_3d'] if params['ocean_model'] is not None: + forcing_group = scenario model = params['ocean_model'] + else: + forcing_group = f"{model}_{scenario}" input_path = os.path.join(base_path_ismip7, "ocean", "tf", ocean_version) file_pattern = (f"tf_{prefix}_{model}_{scenario}_" @@ -187,8 +190,8 @@ def _run_scenario(self): os.remove(f) # Place output in appropriate directory - output_path = os.path.join(output_base_path, "ocean_thermal_forcing", - f"{model}_{scenario}") + output_path = os.path.join(output_base_path, forcing_group, + "ocean_thermal_forcing") if not os.path.exists(output_path): os.makedirs(output_path) diff --git a/docs/users_guide/landice/test_groups/ismip7_forcing.rst b/docs/users_guide/landice/test_groups/ismip7_forcing.rst index 0f9ec60eba..6afb6b7d92 100644 --- a/docs/users_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/users_guide/landice/test_groups/ismip7_forcing.rst @@ -77,6 +77,31 @@ The ``[ismip7] model`` option is ignored for OCX (set it to ``None``), so a single config file processes both the ``atmosphere`` and ``ocean_thermal`` test cases, just like the ESM scenarios. +.. _landice_ismip7_forcing_output: + +Output Layout +------------- + +Processed forcing is written under ``output_base_path`` in a layout that the +:ref:`landice_ismip7_run` test group can ingest directly: + +.. code-block:: none + + {output_base_path}/{group}/atmosphere/{mesh}_SMB_{source}_{scenario}_{years}.nc + {output_base_path}/{group}/atmosphere/{mesh}_temperature_{source}_{scenario}_{years}.nc + {output_base_path}/{group}/atmosphere/{mesh}_runoff_... (and the two gradients) + {output_base_path}/{group}/ocean_thermal_forcing/{mesh}_thermal_forcing_{source}_{scenario}_{years}.nc + +The ``group`` directory is ``{model}_{scenario}`` for the ESM scenarios and +``{scenario}`` (i.e. ``OCX``) for OCX, whose atmosphere and ocean use +different sources (``RACMO2.3p2-ERA`` and ``EN4``) but must share one +directory. To feed :ref:`landice_ismip7_run`, point its ``forcing_basepath`` +at ``output_base_path`` for ESM scenarios, or its ``ocx_forcing_path`` at +``{output_base_path}/OCX`` for OCX. + +Process only one year range into a given ``group`` directory: the run setup +expects exactly one file per forcing field there. + .. _landice_ismip7_forcing_input_data: Input Data From 828699bf0d987c1738ef1286103516df9b7de1a4 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 3 Sep 2026 11:55:35 -0700 Subject: [PATCH 04/17] Fix missing-value leakage in ISMIP7 atmosphere forcing; share extrapolation helper process_smb.py, process_temperature.py, process_smb_gradient.py, and process_temperature_gradient.py remapped the raw source file directly with ncremap, unlike process_runoff.py, ocean_thermal, and the fracture pathway, which all extrapolate fill/missing values on the source grid first. Without that step, no-data cells near the ice-sheet margin (e.g. ocean cells in the RACMO OCX source) leak the ~9.97e36 netCDF fill sentinel into valid neighboring cells during conservative/bilinear remapping. Add the same pre-remap extrapolation to the four affected steps. Also promote the previously-duplicated extrapolate_source() (from fracture/remap_utils.py, formerly re-implemented again in ocean_thermal/process_thermal_forcing.py and process_runoff.py) to a single shared compass.landice.tests.ismip7_forcing.remap_utils module, and update all five call sites to use it instead of per-file private copies. --- .../atmosphere/process_runoff.py | 59 +---------------- .../ismip7_forcing/atmosphere/process_smb.py | 12 +++- .../atmosphere/process_smb_gradient.py | 13 +++- .../atmosphere/process_temperature.py | 12 +++- .../process_temperature_gradient.py | 12 +++- .../fracture/process_excess_melt.py | 2 +- .../fracture/process_lake_properties.py | 2 +- .../ismip7_forcing/fracture/remap_utils.py | 59 ----------------- .../ocean_thermal/process_thermal_forcing.py | 65 +------------------ .../tests/ismip7_forcing/remap_utils.py | 64 ++++++++++++++++++ 10 files changed, 116 insertions(+), 184 deletions(-) create mode 100644 compass/landice/tests/ismip7_forcing/remap_utils.py diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py index 3f0e7b9e5f..801844edf0 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -2,16 +2,15 @@ import os import shutil -import numpy as np import xarray as xr from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call -from scipy.ndimage import distance_transform_edt from compass.landice.tests.ismip7_forcing.create_mapfile import ( build_mapping_file, ) from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step @@ -132,8 +131,7 @@ def run(self): # so they don't pollute neighboring cells during interpolation extrap_file = f"extrap_{basename}" if not os.path.exists(extrap_file): - self._extrapolate_source(input_file, extrap_file, "mrro", - logger) + extrapolate_source(input_file, extrap_file, "mrro", logger) logger.info(f" Remapping: {basename}") args = ["ncremap", @@ -235,56 +233,3 @@ def _combine_and_rename(self, remapped_files, output_file): ds = ds.drop_vars("Time") write_netcdf(ds, output_file) - - def _extrapolate_source(self, input_file, output_file, varname, logger): - """ - Extrapolate fill/missing values on the source polar stereographic - grid using nearest-neighbor via distance_transform_edt. This must - be done before remapping so that fill values don't contaminate the - interpolation stencil. - - Parameters - ---------- - input_file : str - Path to the input NetCDF file on the source grid - - output_file : str - Path to write the extrapolated file - - varname : str - Name of the variable to extrapolate (e.g., "mrro") - - logger : logging.Logger - Logger for status messages - """ - logger.info(f" Extrapolating fill values on source grid: " - f"{os.path.basename(input_file)}") - - ds = xr.open_dataset(input_file, engine="netcdf4") - data = ds[varname] - - # Process each time step - # Source files have dims like (time, y, x) - values = data.values.copy() - non_spatial_shape = values.shape[:-2] # (time,) - - for idx in np.ndindex(non_spatial_shape): - slab = values[idx] # shape (ny, nx) - valid_mask = np.isfinite(slab) - if valid_mask.all() or not valid_mask.any(): - continue - nearest_inds = distance_transform_edt( - ~valid_mask, return_distances=False, return_indices=True) - invalid = ~valid_mask - values[idx][invalid] = slab[ - nearest_inds[0, invalid], - nearest_inds[1, invalid]] - - ds[varname] = (data.dims, values) - ds[varname].attrs = data.attrs - - # Remove _FillValue encoding so output has no masked values - if "_FillValue" in ds[varname].encoding: - del ds[varname].encoding["_FillValue"] - - write_netcdf(ds, output_file) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py index 04de34180f..3bb864710e 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -10,6 +10,7 @@ build_mapping_file, ) from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step @@ -126,15 +127,24 @@ def run(self): logger.info(f" Remapped file exists, skipping: {basename}") continue + # Extrapolate fill values on source grid before remapping + # so they don't pollute neighboring cells during interpolation + extrap_file = f"extrap_{basename}" + if not os.path.exists(extrap_file): + extrapolate_source(input_file, extrap_file, "acabf", logger) + logger.info(f" Remapping: {basename}") args = ["ncremap", - "-i", input_file, + "-i", extrap_file, "-o", remapped_file, "-m", mapping_file, "-v", "acabf"] check_call(args, logger=logger) + # Clean up extrapolated source file + os.remove(extrap_file) + # Combine remapped files and rename to MALI conventions logger.info("Combining remapped files and renaming variables...") output_file = (f"{mali_mesh_name}_SMB_{model}_{scenario}_" diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py index d6c291e4d5..d44192d14c 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -10,6 +10,7 @@ build_mapping_file, ) from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step @@ -127,15 +128,25 @@ def run(self): logger.info(f" Remapped file exists, skipping: {basename}") continue + # Extrapolate fill values on source grid before remapping + # so they don't pollute neighboring cells during interpolation + extrap_file = f"extrap_{basename}" + if not os.path.exists(extrap_file): + extrapolate_source(input_file, extrap_file, "dacabfdz", + logger) + logger.info(f" Remapping: {basename}") args = ["ncremap", - "-i", input_file, + "-i", extrap_file, "-o", remapped_file, "-m", mapping_file, "-v", "dacabfdz"] check_call(args, logger=logger) + # Clean up extrapolated source file + os.remove(extrap_file) + # Combine remapped files and rename to MALI conventions logger.info("Combining remapped files and renaming variables...") output_file = (f"{mali_mesh_name}_SMB_gradient_{model}_{scenario}_" diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index 3168c69f28..24c9529342 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -10,6 +10,7 @@ build_mapping_file, ) from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step @@ -126,15 +127,24 @@ def run(self): logger.info(f" Remapped file exists, skipping: {basename}") continue + # Extrapolate fill values on source grid before remapping + # so they don't pollute neighboring cells during interpolation + extrap_file = f"extrap_{basename}" + if not os.path.exists(extrap_file): + extrapolate_source(input_file, extrap_file, "ts", logger) + logger.info(f" Remapping: {basename}") args = ["ncremap", - "-i", input_file, + "-i", extrap_file, "-o", remapped_file, "-m", mapping_file, "-v", "ts"] check_call(args, logger=logger) + # Clean up extrapolated source file + os.remove(extrap_file) + # Combine remapped files and rename to MALI conventions logger.info("Combining remapped files and renaming variables...") output_file = (f"{mali_mesh_name}_temperature_{model}_{scenario}_" diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py index 16181f02e5..45d39b9a1d 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -10,6 +10,7 @@ build_mapping_file, ) from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step @@ -128,15 +129,24 @@ def run(self): logger.info(f" Remapped file exists, skipping: {basename}") continue + # Extrapolate fill values on source grid before remapping + # so they don't pollute neighboring cells during interpolation + extrap_file = f"extrap_{basename}" + if not os.path.exists(extrap_file): + extrapolate_source(input_file, extrap_file, "dtsdz", logger) + logger.info(f" Remapping: {basename}") args = ["ncremap", - "-i", input_file, + "-i", extrap_file, "-o", remapped_file, "-m", mapping_file, "-v", "dtsdz"] check_call(args, logger=logger) + # Clean up extrapolated source file + os.remove(extrap_file) + # Combine remapped files and rename to MALI conventions logger.info("Combining remapped files and renaming variables...") output_file = (f"{mali_mesh_name}_temperature_gradient_{model}_" diff --git a/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py b/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py index fbf3f75b83..8a9275c962 100644 --- a/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py +++ b/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py @@ -12,9 +12,9 @@ ) from compass.landice.tests.ismip7_forcing.fracture.remap_utils import ( add_xtime_and_write, - extrapolate_source, open_rename_and_trim, ) +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step diff --git a/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py b/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py index 10df254c69..6867f20735 100644 --- a/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py +++ b/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py @@ -9,9 +9,9 @@ ) from compass.landice.tests.ismip7_forcing.fracture.remap_utils import ( add_xtime_and_write, - extrapolate_source, open_rename_and_trim, ) +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step diff --git a/compass/landice/tests/ismip7_forcing/fracture/remap_utils.py b/compass/landice/tests/ismip7_forcing/fracture/remap_utils.py index 8da5de4ae4..14518c7c79 100644 --- a/compass/landice/tests/ismip7_forcing/fracture/remap_utils.py +++ b/compass/landice/tests/ismip7_forcing/fracture/remap_utils.py @@ -1,67 +1,8 @@ """ Shared helpers for remapping ISMIP7 fracture forcing data to the MALI mesh. """ -import os - -import numpy as np import xarray as xr from mpas_tools.io import write_netcdf -from scipy.ndimage import distance_transform_edt - - -def extrapolate_source(input_file, output_file, varnames, logger): - """ - Extrapolate fill/missing values on the source polar stereographic grid - using nearest-neighbor via ``distance_transform_edt``. This must be done - before remapping so that fill values don't contaminate the interpolation - stencil. - - Parameters - ---------- - input_file : str - Path to the input NetCDF file on the source grid - - output_file : str - Path to write the extrapolated file - - varnames : str or list of str - Name(s) of the variable(s) to extrapolate - - logger : logging.Logger - Logger for status messages - """ - if isinstance(varnames, str): - varnames = [varnames] - - logger.info(f" Extrapolating fill values on source grid: " - f"{os.path.basename(input_file)}") - - ds = xr.open_dataset(input_file, decode_times=False) - - for varname in varnames: - data = ds[varname] - values = data.values.copy() - non_spatial_shape = values.shape[:-2] - - for idx in np.ndindex(non_spatial_shape): - slab = values[idx] - valid_mask = np.isfinite(slab) - if valid_mask.all() or not valid_mask.any(): - continue - nearest_inds = distance_transform_edt( - ~valid_mask, return_distances=False, return_indices=True) - invalid = ~valid_mask - values[idx][invalid] = slab[ - nearest_inds[0, invalid], - nearest_inds[1, invalid]] - - ds[varname] = (data.dims, values) - ds[varname].attrs = data.attrs - if "_FillValue" in ds[varname].encoding: - del ds[varname].encoding["_FillValue"] - - write_netcdf(ds, output_file) - ds.close() def open_rename_and_trim(remapped_file, rename_vars, start_year, end_year): diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py index fba0053b1d..c0eed44c97 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -2,16 +2,15 @@ import os import shutil -import numpy as np import xarray as xr from mpas_tools.io import write_netcdf from mpas_tools.logging import check_call -from scipy.ndimage import distance_transform_edt from compass.landice.tests.ismip7_forcing.create_mapfile import ( build_mapping_file, ) from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step @@ -156,8 +155,7 @@ def _run_scenario(self): # so they don't pollute neighboring cells during interpolation extrap_file = f"extrap_{basename}" if not os.path.exists(extrap_file): - self._extrapolate_source(input_file, extrap_file, "tf", - logger) + extrapolate_source(input_file, extrap_file, "tf", logger) logger.info(f" Remapping: {basename}") args = ["ncremap", @@ -252,8 +250,7 @@ def _run_climatology(self): if not os.path.exists(remapped_file): extrap_file = f"extrap_{basename}" if not os.path.exists(extrap_file): - self._extrapolate_source(input_file, extrap_file, "tf", - logger) + extrapolate_source(input_file, extrap_file, "tf", logger) logger.info(f" Remapping: {basename}") args = ["ncremap", @@ -571,59 +568,3 @@ def _rename_climatology_3d(self, remapped_file, output_file): ds = ds.drop_vars("nISMIP6OceanLayers") write_netcdf(ds, output_file) - - def _extrapolate_source(self, input_file, output_file, varname, logger): - """ - Extrapolate fill/missing values on the source polar stereographic - grid using nearest-neighbor interpolation from valid cells. This - must be done before remapping so that fill values don't contaminate - the interpolation stencil. - - Parameters - ---------- - input_file : str - Path to the input NetCDF file on the source grid - - output_file : str - Path to write the extrapolated file - - varname : str - Name of the variable to extrapolate (e.g., "tf") - - logger : logging.Logger - Logger for status messages - """ - logger.info(f" Extrapolating fill values on source grid: " - f"{os.path.basename(input_file)}") - - ds = xr.open_dataset(input_file, engine="netcdf4") - data = ds[varname] - - # Process each time step (and z level if 3D) - # Source files have dims like (time, z, y, x) or (time, y, x) - values = data.values.copy() - non_spatial_shape = values.shape[:-2] # (time,) or (time, z) - - # Use distance_transform_edt with return_indices to find the - # nearest valid cell index for each invalid cell. This is O(n) - # on the grid and much faster than KD-tree approaches. - for idx in np.ndindex(non_spatial_shape): - slab = values[idx] # shape (ny, nx) - valid_mask = np.isfinite(slab) - if valid_mask.all() or not valid_mask.any(): - continue - nearest_inds = distance_transform_edt( - ~valid_mask, return_distances=False, return_indices=True) - invalid = ~valid_mask - values[idx][invalid] = slab[ - nearest_inds[0, invalid], - nearest_inds[1, invalid]] - - ds[varname] = (data.dims, values) - ds[varname].attrs = data.attrs - - # Remove _FillValue encoding so output has no masked values - if "_FillValue" in ds[varname].encoding: - del ds[varname].encoding["_FillValue"] - - write_netcdf(ds, output_file) diff --git a/compass/landice/tests/ismip7_forcing/remap_utils.py b/compass/landice/tests/ismip7_forcing/remap_utils.py new file mode 100644 index 0000000000..e54f9e9fbf --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/remap_utils.py @@ -0,0 +1,64 @@ +""" +Shared helpers for remapping ISMIP7 forcing data to the MALI mesh. +""" +import os + +import numpy as np +import xarray as xr +from mpas_tools.io import write_netcdf +from scipy.ndimage import distance_transform_edt + + +def extrapolate_source(input_file, output_file, varnames, logger): + """ + Extrapolate fill/missing values on the source polar stereographic grid + using nearest-neighbor via ``distance_transform_edt``. This must be done + before remapping so that fill values don't contaminate the interpolation + stencil. + + Parameters + ---------- + input_file : str + Path to the input NetCDF file on the source grid + + output_file : str + Path to write the extrapolated file + + varnames : str or list of str + Name(s) of the variable(s) to extrapolate + + logger : logging.Logger + Logger for status messages + """ + if isinstance(varnames, str): + varnames = [varnames] + + logger.info(f" Extrapolating fill values on source grid: " + f"{os.path.basename(input_file)}") + + ds = xr.open_dataset(input_file, decode_times=False) + + for varname in varnames: + data = ds[varname] + values = data.values.copy() + non_spatial_shape = values.shape[:-2] + + for idx in np.ndindex(non_spatial_shape): + slab = values[idx] + valid_mask = np.isfinite(slab) + if valid_mask.all() or not valid_mask.any(): + continue + nearest_inds = distance_transform_edt( + ~valid_mask, return_distances=False, return_indices=True) + invalid = ~valid_mask + values[idx][invalid] = slab[ + nearest_inds[0, invalid], + nearest_inds[1, invalid]] + + ds[varname] = (data.dims, values) + ds[varname].attrs = data.attrs + if "_FillValue" in ds[varname].encoding: + del ds[varname].encoding["_FillValue"] + + write_netcdf(ds, output_file) + ds.close() From 7047302c97b272fed33039a3aa7d32973f4380f7 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 4 Sep 2026 07:54:58 -0700 Subject: [PATCH 05/17] Fix stale autosummary reference after remap_utils refactor extrapolate_source moved from fracture.remap_utils to the shared compass.landice.tests.ismip7_forcing.remap_utils module (828699bf0). Update the api.rst autosummary entry accordingly; this was breaking the strict Sphinx build (autosummary failed to import the moved symbol). --- docs/developers_guide/landice/api.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index d2cd995145..75bb6bab8c 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -385,6 +385,7 @@ ismip7_forcing configure.configure ice_sheet_params.get_params create_mapfile.build_mapping_file + remap_utils.extrapolate_source atmosphere.Atmosphere atmosphere.Atmosphere.configure @@ -412,7 +413,6 @@ ismip7_forcing fracture.Fracture fracture.Fracture.configure - fracture.remap_utils.extrapolate_source fracture.remap_utils.open_rename_and_trim fracture.remap_utils.add_xtime_and_write fracture.process_excess_melt.ProcessExcessMelt From 6176864e95c04f8da08f418f767a071b5c028e14 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 4 Sep 2026 13:56:13 -0700 Subject: [PATCH 06/17] Add optional GrIS 3-D ocean thermal forcing step to ismip7_forcing Port the standalone greenland_thermal_forcing tool into a new optional step of the ocean_thermal test case. build_3d_thermal_forcing (GrIS only, gated by process_ocean_thermal_3d) converts the 2-D GrIS thermal forcing into a 30-level 3-D field for MALI's nonlocal (Jourdain et al. 2020) melt scheme: seven regional EN4 vertical profiles, per-cell seafloor anchoring to the 2-D forcing, and per-region deltaT calibration. It auto-chains from the 2-D output the same run just produced and writes ismip6shelfMelt_3dThermalForcing (plus deltaT/gamma0/zOcean/basin), supplementing the 2-D file. The ported science lives in ocean_thermal/greenland_3d.py; the 3-D-specific parameters are supplied via a JSON config (config_file), while mesh, 2-D forcing, output, and diagnostics paths are injected from the compass config. Also standardize ocean scenario output names to 2dThermalForcing (GrIS) / 3dThermalForcing (AIS) and update the ismip7_run tf globs accordingly (AIS ingests 3-D, GrIS 2-D). Adds docs and no new dependencies (shapely/h5py/dask already present). Not yet wired: ismip7_run GrIS streams to ingest the 3-D field (follow-up). --- .../tests/ismip7_forcing/ismip7_forcing.cfg | 12 + .../ismip7_forcing/ocean_thermal/__init__.py | 4 + .../ocean_thermal/build_3d_thermal_forcing.py | 99 + .../ocean_thermal/greenland_3d.py | 1637 +++++++++++++++++ .../ocean_thermal/process_thermal_forcing.py | 3 +- .../ismip7_ais/set_up_experiment.py | 2 +- .../ismip7_gris/set_up_experiment.py | 2 +- docs/developers_guide/landice/api.rst | 4 + .../landice/test_groups/ismip7_forcing.rst | 29 +- .../landice/test_groups/ismip7_forcing.rst | 50 +- 10 files changed, 1831 insertions(+), 11 deletions(-) create mode 100644 compass/landice/tests/ismip7_forcing/ocean_thermal/build_3d_thermal_forcing.py create mode 100644 compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg index f02b79067a..11e67696bd 100644 --- a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg @@ -36,6 +36,10 @@ process_ocean_thermal = true # Whether to process observational ocean thermal forcing climatology process_ocean_climatology = true +# Whether to build regional 3-D Greenland ocean thermal forcing from the 2-D +# forcing (GrIS only). Antarctica already produces 3-D thermal forcing. +process_ocean_thermal_3d = false + # config options for ismip7 atmosphere forcing [ismip7_atmosphere] @@ -70,6 +74,14 @@ method_remap = bilinear # (directory containing tf/, so/, thetao/ subdirs) base_path_climatology = /path/to/ISMIP7/forcing/AIS/obs/zhou_annual_06_nov +# config options for building 3-D Greenland ocean thermal forcing +[ismip7_ocean_thermal_3d] + +# Path to the JSON config with the 3-D-specific parameters (EN4 directory, +# region-mask file, source-region GeoJSON, calibration, vertical grid, etc.). +# User must supply when process_ocean_thermal_3d is true. +config_file = NotAvailable + # config options for ismip7 fracture (Path C, ice shelf collapse) forcing [ismip7_fracture] diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py index 2a97700428..c4fd353374 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py @@ -1,6 +1,9 @@ from compass.landice.tests.ismip7_forcing.configure import ( configure as configure_testgroup, ) +from compass.landice.tests.ismip7_forcing.ocean_thermal.build_3d_thermal_forcing import ( # noqa: E501 + BuildGreenland3dThermalForcing, +) from compass.landice.tests.ismip7_forcing.ocean_thermal.process_thermal_forcing import ( # noqa: E501 ProcessThermalForcing, ) @@ -30,6 +33,7 @@ def __init__(self, test_group): super().__init__(test_group=test_group, name=name, subdir=subdir) self.add_step(ProcessThermalForcing(test_case=self)) + self.add_step(BuildGreenland3dThermalForcing(test_case=self)) def configure(self): """ diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/build_3d_thermal_forcing.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/build_3d_thermal_forcing.py new file mode 100644 index 0000000000..580f4d2a99 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/build_3d_thermal_forcing.py @@ -0,0 +1,99 @@ +import os +from pathlib import Path + +from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.tests.ismip7_forcing.ocean_thermal import greenland_3d +from compass.step import Step + + +class BuildGreenland3dThermalForcing(Step): + """ + A step that builds regional 3-D Greenland ocean thermal forcing from the + 2-D forcing produced by ProcessThermalForcing. GrIS only; gated by the + ``process_ocean_thermal_3d`` config option. The 3-D-specific parameters + are supplied through a JSON config file (``config_file``); the mesh, 2-D + forcing, output, and diagnostics paths are injected from the compass + config so the step auto-chains from the 2-D forcing. + """ + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_forcing.ocean_thermal.OceanThermal + The test case this step belongs to + """ # noqa: E501 + super().__init__(test_case=test_case, + name="build_3d_thermal_forcing") + + def run(self): + """ + Run this step of the test case + """ + config = self.config + logger = self.logger + section = config["ismip7"] + + if not section.getboolean("process_ocean_thermal_3d"): + logger.info("process_ocean_thermal_3d is false; skipping 3-D " + "Greenland thermal forcing.") + return + + ice_sheet = section.get("ice_sheet") + if ice_sheet != "gis": + raise ValueError( + "process_ocean_thermal_3d is only supported for the Greenland " + "Ice Sheet (ice_sheet = gis); the Antarctic pathway already " + "produces 3-D thermal forcing directly.") + + base_path_mali = section.get("base_path_mali") + mali_mesh_file = section.get("mali_mesh_file") + mali_mesh_name = section.get("mali_mesh_name") + model = section.get("model") + scenario = section.get("scenario") + output_base_path = section.get("output_base_path") + + # Mirror ProcessThermalForcing._run_scenario forcing_group and ocean + # source so we find the 2-D forcing it just wrote. + params = get_params(config) + if params["ocean_model"] is not None: + forcing_group = scenario + source = params["ocean_model"] + else: + forcing_group = f"{model}_{scenario}" + source = model + + ocean_section = config["ismip7_ocean_thermal"] + start_year = ocean_section.getint("start_year") + end_year = ocean_section.getint("end_year") + + ocean_dir = os.path.join(output_base_path, forcing_group, + "ocean_thermal_forcing") + forcing_2d = os.path.join( + ocean_dir, + f"{mali_mesh_name}_2dThermalForcing_{source}_{scenario}_" + f"{start_year}-{end_year}.nc") + output_file = os.path.join( + ocean_dir, + f"{mali_mesh_name}_3dThermalForcing_{source}_{scenario}_" + f"{start_year}-{end_year}.nc") + + json_path = config.get("ismip7_ocean_thermal_3d", "config_file") + if json_path == "NotAvailable": + raise ValueError( + "You need to supply the [ismip7_ocean_thermal_3d] config_file " + "option (path to the 3-D forcing JSON config) when " + "process_ocean_thermal_3d is true") + + overrides = { + "mesh_file": Path(os.path.join(base_path_mali, mali_mesh_file)), + "forcing_2d_file": Path(forcing_2d), + "output_file": Path(output_file), + "diagnostics_directory": Path( + os.path.join(self.work_dir, "diagnostics_3d")), + } + cfg = greenland_3d.Config.from_json(Path(json_path), + overrides=overrides) + greenland_3d.run(cfg, logger) diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py new file mode 100644 index 0000000000..df24ef2ef4 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py @@ -0,0 +1,1637 @@ +"""Build regional 3-D Greenland ocean thermal forcing for MALI. + +Ported from the standalone ``greenland_thermal_forcing.py`` tool. The workflow +has three coupled stages: + +1. Construct seven regional climatological profiles from monthly EN4 objective + analyses. +2. Translate each regional profile cell-by-cell so it matches processed + ISMIP7 OCX thermal forcing at the effective local seafloor. +3. Hold gamma0 fixed and calibrate one Jourdain et al. (2020) nonlocal + temperature correction (deltaT) per region. + +The large time-dependent output is streamed one record at a time and written +as NETCDF3_64BIT_OFFSET. + +The 3-D-specific parameters are supplied through a JSON config file, while the +mesh, 2-D forcing, output, and diagnostics paths are injected by the compass +step (see ``build_3d_thermal_forcing.py``). +""" + +from __future__ import annotations + +import csv +import json +import math +import os +import re +import sys +import warnings +from dataclasses import dataclass +from datetime import datetime, timezone +from pathlib import Path +from typing import Any, Sequence + +import numpy as np +from scipy.optimize import brentq +from scipy.spatial import cKDTree + +REGION_NAMES = ( + "ISMIP6 Greenland Central East Shelf", + "ISMIP6 Greenland Central West Shelf", + "ISMIP6 Greenland North East Shelf", + "ISMIP6 Greenland North Shelf", + "ISMIP6 Greenland North West Shelf", + "ISMIP6 Greenland South East Shelf", + "ISMIP6 Greenland South West Shelf", +) + +REGION_KEYS = ( + "central_east", + "central_west", + "north_east", + "north", + "north_west", + "south_east", + "south_west", +) + +EN4_BIAS_CORRECTIONS = {"g10", "l09", "c13", "c14", "unknown"} +DATE_RE = re.compile(r"(? Any: + if key not in mapping or mapping[key] in (None, ""): + raise ValueError(f"Missing required configuration value: {key}") + return mapping[key] + + +def _as_path(value: str, base: Path) -> Path: + path = Path(value).expanduser() + return path if path.is_absolute() else (base / path).resolve() + + +@dataclass(frozen=True) +class Config: + mesh_file: Path + region_mask_file: Path + forcing_2d_file: Path + en4_directory: Path + en4_source_region_geojson: Path + output_file: Path + diagnostics_directory: Path + ocean_levels_m: np.ndarray + source_max_depth_m: float + profile_start_year: int + profile_end_year: int + calibration_start_year: int + calibration_end_year: int + en4_version: str + en4_bias_correction: str + en4_file_glob: str + en4_max_mesh_distance_km: float + en4_latitude_min: float + en4_latitude_max: float + gamma0_m_per_yr: float + regional_melt_targets_m_per_yr: np.ndarray + rho_ice: float + rho_seawater: float + cp_seawater: float + latent_heat_ice: float + flotation_tolerance_m: float + minimum_ice_thickness_m: float + freezing_a: float + freezing_b: float + freezing_c: float + forcing_variable: str + overwrite: bool + + @classmethod + def from_json(cls, path: Path, overrides: dict | None = None) -> "Config": + """Build a Config from JSON, optionally injecting compass paths. + + ``overrides`` may supply ``mesh_file``, ``forcing_2d_file``, + ``output_file``, and ``diagnostics_directory`` (as ``Path`` objects). + When supplied, they take precedence over the corresponding JSON + ``files`` entries, which then become optional. The region-mask, EN4, + and GeoJSON paths always come from the JSON. + """ + overrides = overrides or {} + with path.open("r", encoding="utf-8") as handle: + raw = json.load(handle) + base = path.resolve().parent + + files = _required(raw, "files") + en4 = raw.get("en4", {}) + calibration = raw.get("calibration", {}) + physical = raw.get("physical_constants", {}) + freezing = raw.get("freezing_point", {}) + + def resolved(field_name, json_key, required=True, default=None): + if overrides.get(field_name) is not None: + return Path(overrides[field_name]) + if required: + return _as_path(_required(files, json_key), base) + value = files.get(json_key, default) + return _as_path(value, base) if value else None + + if "ocean_levels_m" in raw: + levels = np.asarray(raw["ocean_levels_m"], dtype=float) + else: + vertical_grid = raw.get("ocean_vertical_grid", {}) + number_of_levels = int(vertical_grid.get("number_of_levels", 30)) + surface_m = float(vertical_grid.get("surface_m", 0.0)) + bottom_m = float(vertical_grid.get("bottom_m", -1000.0)) + if number_of_levels < 2: + raise ValueError( + "ocean_vertical_grid.number_of_levels must be at least 2" + ) + levels = np.linspace(surface_m, bottom_m, number_of_levels) + validate_ocean_levels(levels) + + targets_raw = calibration.get("regional_melt_targets_m_per_yr", 20.0) + if isinstance(targets_raw, dict): + missing = [key for key in REGION_KEYS if key not in targets_raw] + if missing: + raise ValueError( + f"Missing regional melt targets for: {', '.join(missing)}" + ) + targets = np.asarray( + [targets_raw[key] for key in REGION_KEYS], dtype=float + ) + elif np.isscalar(targets_raw): + targets = np.full(len(REGION_NAMES), float(targets_raw)) + else: + targets = np.asarray(targets_raw, dtype=float) + if targets.shape != (len(REGION_NAMES),) or np.any(targets < 0.0): + raise ValueError( + "Regional melt targets must be seven nonnegative values" + ) + + bias = str(en4.get("bias_correction", "unknown")).lower() + if bias not in EN4_BIAS_CORRECTIONS: + raise ValueError( + f"EN4 bias correction must be one of " + f"{sorted(EN4_BIAS_CORRECTIONS)}" + ) + + cfg = cls( + mesh_file=resolved("mesh_file", "mesh"), + region_mask_file=_as_path(_required(files, "region_masks"), base), + forcing_2d_file=resolved("forcing_2d_file", "forcing_2d"), + en4_directory=_as_path(_required(files, "en4_directory"), base), + en4_source_region_geojson=_as_path( + _required(files, "en4_source_region_geojson"), base + ), + output_file=resolved("output_file", "output"), + diagnostics_directory=resolved( + "diagnostics_directory", "diagnostics", + required=False, default="diagnostics" + ), + ocean_levels_m=levels, + source_max_depth_m=float(raw.get("source_max_depth_m", 1000.0)), + profile_start_year=int(en4.get("profile_start_year", 2007)), + profile_end_year=int(en4.get("profile_end_year", 2015)), + calibration_start_year=int(calibration.get("start_year", 2007)), + calibration_end_year=int(calibration.get("end_year", 2015)), + en4_version=str(en4.get("version", "EN.4.2.2")), + en4_bias_correction=bias, + en4_file_glob=str(en4.get("file_glob", "**/*.nc")), + en4_max_mesh_distance_km=float( + en4.get("max_mesh_distance_km", 300.0) + ), + en4_latitude_min=float(en4.get("latitude_min", 55.0)), + en4_latitude_max=float(en4.get("latitude_max", 90.0)), + gamma0_m_per_yr=float(calibration.get("gamma0_m_per_yr", 14500.0)), + regional_melt_targets_m_per_yr=targets, + rho_ice=float(physical.get("rho_ice", 910.0)), + rho_seawater=float(physical.get("rho_seawater", 1028.0)), + cp_seawater=float(physical.get("cp_seawater", 3974.0)), + latent_heat_ice=float(physical.get("latent_heat_ice", 335000.0)), + flotation_tolerance_m=float( + calibration.get("flotation_tolerance_m", 1.0) + ), + minimum_ice_thickness_m=float( + calibration.get("minimum_ice_thickness_m", 0.0) + ), + freezing_a=float(freezing.get("a_degC_per_salinity", -0.0575)), + freezing_b=float(freezing.get("b_degC", 0.0901)), + freezing_c=float(freezing.get("c_degC_per_m", 7.61e-4)), + forcing_variable=str( + raw.get("forcing_2d_variable", "ismip6_2dThermalForcing") + ), + overwrite=bool(raw.get("overwrite", False)), + ) + cfg.validate() + return cfg + + def validate(self) -> None: + for label, path in ( + ("mesh", self.mesh_file), + ("region-mask", self.region_mask_file), + ("2-D forcing", self.forcing_2d_file), + ("EN4 directory", self.en4_directory), + ("EN4 source-region GeoJSON", self.en4_source_region_geojson), + ): + if not path.exists(): + raise FileNotFoundError( + f"Configured {label} path does not exist: {path}" + ) + if self.profile_start_year > self.profile_end_year: + raise ValueError( + "EN4 profile_start_year must not exceed profile_end_year" + ) + if self.calibration_start_year > self.calibration_end_year: + raise ValueError( + "Calibration start_year must not exceed end_year" + ) + if self.source_max_depth_m <= 0.0: + raise ValueError("source_max_depth_m must be positive") + if self.gamma0_m_per_yr <= 0.0: + raise ValueError("gamma0_m_per_yr must be positive") + if self.en4_bias_correction == "unknown": + warnings.warn( + "EN4 bias correction is unknown. Processing will continue " + "only if the discovered files contain at most one analysis " + "per month.", + stacklevel=2, + ) + + +@dataclass +class MeshData: + bed: np.ndarray + thickness: np.ndarray + area: np.ndarray + lat_deg: np.ndarray + lon_deg: np.ndarray + + +@dataclass +class RegionalProfiles: + source_z_m: np.ndarray + output_z_m: np.ndarray + monthly_dates: list[str] + monthly_temperature_degC: np.ndarray + monthly_salinity: np.ndarray + monthly_thermal_forcing_degC: np.ndarray + temperature_degC: np.ndarray + salinity: np.ndarray + freezing_temperature_degC: np.ndarray + thermal_forcing_degC: np.ndarray + output_thermal_forcing_degC: np.ndarray + valid_gridpoint_counts: np.ndarray + temperature_observation_influence: np.ndarray + salinity_observation_influence: np.ndarray + mapped_lats: np.ndarray + mapped_lons: np.ndarray + mapped_basins: np.ndarray + mapped_distances_km: np.ndarray + + +def validate_ocean_levels(levels: np.ndarray) -> None: + if levels.ndim != 1 or levels.size < 2: + raise ValueError("ocean_levels_m must contain at least two values") + if not np.all(np.isfinite(levels)): + raise ValueError("ocean_levels_m contains non-finite values") + if not np.all(np.diff(levels) < 0.0): + raise ValueError( + "ocean_levels_m must be strictly decreasing (negative downward)" + ) + if levels[0] > 0.0 or levels[-1] >= 0.0: + raise ValueError( + "ocean levels must begin at or below 0 m and extend below 0 m" + ) + + +def decode_char_rows(values: np.ndarray) -> list[str]: + """Decode either raw S1 character rows or xarray-concatenated strings.""" + arr = np.asarray(values) + if arr.ndim == 0: + arr = arr.reshape(1) + if arr.ndim == 1 and (arr.dtype.kind == "U" or arr.dtype.itemsize > 1): + result = [] + for value in arr: + if isinstance(value, bytes): + text = value.decode("utf-8", errors="replace") + else: + text = str(value) + result.append(text.rstrip("\x00 ")) + return result + if arr.ndim == 1: + arr = arr[None, :] + result: list[str] = [] + for row in arr: + if row.dtype.kind == "S": + text = b"".join(row.tolist()).decode("utf-8", errors="replace") + elif row.dtype.kind == "U": + text = "".join(row.tolist()) + else: + text = bytes(row.tolist()).decode("utf-8", errors="replace") + result.append(text.rstrip("\x00 ")) + return result + + +def decode_region_names(values: np.ndarray) -> tuple[str, ...]: + return tuple(decode_char_rows(values)) + + +def radians_or_degrees_to_degrees(values: np.ndarray, kind: str) -> np.ndarray: + values = np.asarray(values, dtype=float) + limit = math.pi / 2 + 1e-6 if kind == "latitude" else 2 * math.pi + 1e-6 + if np.nanmax(np.abs(values)) <= limit: + return np.rad2deg(values) + return values + + +def build_basin_ids(region_masks: np.ndarray) -> np.ndarray: + masks = np.asarray(region_masks) + if masks.ndim != 2 or masks.shape[1] != len(REGION_NAMES): + raise ValueError( + f"regionCellMasks must have shape (nCells, {len(REGION_NAMES)}); " + f"got {masks.shape}" + ) + membership = masks != 0 + counts = membership.sum(axis=1) + overlapping = np.flatnonzero(counts > 1) + unassigned = np.flatnonzero(counts == 0) + if overlapping.size: + raise ValueError( + f"{overlapping.size} cells belong to multiple regions; first " + f"indices: {overlapping[:10].tolist()}" + ) + if unassigned.size: + raise ValueError( + f"{unassigned.size} cells have no region; first indices: " + f"{unassigned[:10].tolist()}" + ) + return np.argmax(membership, axis=1).astype(np.int32) + 1 + + +def unit_sphere_xyz(lat_deg: np.ndarray, lon_deg: np.ndarray) -> np.ndarray: + lat = np.deg2rad(np.asarray(lat_deg, dtype=float)) + lon = np.deg2rad(np.asarray(lon_deg, dtype=float)) + cos_lat = np.cos(lat) + return np.column_stack( + (cos_lat * np.cos(lon), cos_lat * np.sin(lon), np.sin(lat)) + ) + + +def chord_to_arc_km( + chord: np.ndarray, radius_km: float = 6371.0 +) -> np.ndarray: + return 2.0 * radius_km * np.arcsin( + np.clip(np.asarray(chord) / 2.0, 0.0, 1.0) + ) + + +def freezing_temperature( + salinity: np.ndarray, + z_m: np.ndarray, + a: float = -0.0575, + b: float = 0.0901, + c: float = 7.61e-4, +) -> np.ndarray: + return a * np.asarray(salinity) + b + c * np.asarray(z_m) + + +def interpolate_profile( + z_source: np.ndarray, values: np.ndarray, z_target: np.ndarray +) -> np.ndarray: + """Linearly interpolate with constant endpoint extension.""" + z_source = np.asarray(z_source, dtype=float) + values = np.asarray(values, dtype=float) + order = np.argsort(z_source) + good = np.isfinite(z_source[order]) & np.isfinite(values[order]) + if good.sum() < 2: + raise ValueError( + "At least two finite source profile levels are required" + ) + x = z_source[order][good] + y = values[order][good] + return np.interp(np.asarray(z_target, dtype=float), x, y, left=y[0], + right=y[-1]) + + +def profiles_at_cell_depths( + regional_profiles: np.ndarray, + z_levels: np.ndarray, + basin_ids: np.ndarray, + cell_depths: np.ndarray, +) -> np.ndarray: + """Evaluate piecewise-linear regional profiles at cell-specific depths.""" + regional_profiles = np.asarray(regional_profiles, dtype=float) + basin_ids = np.asarray(basin_ids) + cell_depths = np.asarray(cell_depths, dtype=float) + result = np.empty(cell_depths.shape, dtype=float) + for region in range(len(REGION_NAMES)): + mask = basin_ids == region + 1 + if np.any(mask): + result[mask] = interpolate_profile( + z_levels, regional_profiles[region], cell_depths[mask] + ) + return result + + +def floating_mask_and_draft( + bed: np.ndarray, + thickness: np.ndarray, + rho_ice: float, + rho_seawater: float, + tolerance_m: float, + minimum_thickness_m: float, +) -> tuple[np.ndarray, np.ndarray]: + bed = np.asarray(bed, dtype=float) + thickness = np.asarray(thickness, dtype=float) + flotation_thickness = np.where( + bed < 0.0, -bed * rho_seawater / rho_ice, 0.0 + ) + floating = ( + np.isfinite(bed) & + np.isfinite(thickness) & + (thickness > minimum_thickness_m) & + (thickness <= flotation_thickness + tolerance_m) + ) + draft = -rho_ice / rho_seawater * thickness + return floating, draft + + +def nonlocal_mean_melt( + delta_t: float, + monthly_mean_tf: np.ndarray, + gamma0_m_per_yr: float, + coefficient: float, +) -> float: + corrected = np.asarray(monthly_mean_tf, dtype=float) + delta_t + return float( + gamma0_m_per_yr * coefficient**2 * + np.mean(corrected * np.abs(corrected)) + ) + + +def calibrate_delta_t( + monthly_mean_tf: np.ndarray, + target_melt_m_per_yr: float, + gamma0_m_per_yr: float, + coefficient: float, +) -> float: + values = np.asarray(monthly_mean_tf, dtype=float) + values = values[np.isfinite(values)] + if values.size == 0: + raise ValueError("No finite monthly regional thermal-forcing means") + if target_melt_m_per_yr < 0.0: + raise ValueError("Target mean melt must be nonnegative") + + def residual(delta: float) -> float: + return ( + nonlocal_mean_melt(delta, values, gamma0_m_per_yr, coefficient) - + target_melt_m_per_yr + ) + + # The physically relevant root has a nonnegative corrected regional mean + # in every month. Starting above -min(TF) selects that branch. + lower = float(-np.min(values)) + if target_melt_m_per_yr == 0.0: + return lower + upper = max(lower + 1.0, 1.0) + while residual(upper) < 0.0: + upper = lower + 2.0 * (upper - lower) + if upper - lower > 1000.0: + raise RuntimeError("Could not bracket deltaT calibration root") + return float(brentq(residual, lower, upper, xtol=1e-12, rtol=1e-12)) + + +def parse_yyyymm_from_name(path: Path) -> tuple[int, int]: + matches = DATE_RE.findall(path.name) + valid = [(int(y), int(m)) for y, m in matches if 1 <= int(m) <= 12] + if not valid: + raise ValueError(f"Could not find YYYYMM in EN4 filename: {path.name}") + return valid[-1] + + +def discover_en4_files(cfg: Config) -> list[Path]: + candidates = sorted(cfg.en4_directory.glob(cfg.en4_file_glob)) + selected: dict[tuple[int, int], Path] = {} + for path in candidates: + if not path.is_file(): + continue + name_lower = path.name.lower() + if (cfg.en4_version.lower() not in name_lower or + "analysis" not in name_lower): + continue + if (cfg.en4_bias_correction != "unknown" and + f".{cfg.en4_bias_correction}." not in name_lower): + continue + try: + year, month = parse_yyyymm_from_name(path) + except ValueError: + continue + if not (cfg.profile_start_year <= year <= cfg.profile_end_year): + continue + key = (year, month) + if key in selected: + raise ValueError( + f"Multiple EN4 analyses found for {year:04d}-{month:02d}: " + f"{selected[key]} and {path}. Select an explicit bias " + f"correction." + ) + selected[key] = path + if not selected: + raise FileNotFoundError( + "No EN4 objective-analysis files matched the configured " + "directory, glob, years, version, and bias correction" + ) + return [selected[key] for key in sorted(selected)] + + +def _array_with_nan(values: Any) -> np.ndarray: + if np.ma.isMaskedArray(values): + return np.asarray(values.filled(np.nan), dtype=float) + return np.asarray(values, dtype=float) + + +def _find_variable( + dataset: Any, requested: str, alternatives: Sequence[str] = () +) -> Any: + for name in (requested, *alternatives): + if name in dataset.variables: + return dataset.variables[name] + raise KeyError( + f"None of these variables is present: {(requested, *alternatives)}" + ) + + +def _read_3d_en4_variable( + variable: Any, depth_name: str, lat_name: str, lon_name: str +) -> np.ndarray: + target_dims = [depth_name, lat_name, lon_name] + data_array = variable + for dim in tuple(data_array.dims): + if dim not in target_dims: + if data_array.sizes[dim] != 1: + raise ValueError( + f"Unexpected non-singleton EN4 dimension {dim} in " + f"{variable.name}" + ) + data_array = data_array.isel({dim: 0}, drop=True) + if set(data_array.dims) != set(target_dims): + raise ValueError( + f"Cannot arrange {variable.name} dimensions {data_array.dims} as " + f"{target_dims}" + ) + return _array_with_nan(data_array.transpose(*target_dims).values) + + +def _temperature_to_deg_c(values: np.ndarray, units: str) -> np.ndarray: + normalized = units.strip().lower() + if (normalized in {"k", "kelvin", "degrees_k", "degree_k"} or + "kelvin" in normalized): + return values - 273.15 + if "c" in normalized or normalized == "": + return values + raise ValueError(f"Unsupported EN4 temperature units: {units!r}") + + +def points_in_geojson( + longitude_deg: np.ndarray, + latitude_deg: np.ndarray, + geojson_path: Path, +) -> np.ndarray: + """Return points inside or on the boundary of a WGS84 GeoJSON geometry.""" + try: + from shapely import covers, points + from shapely.geometry import shape + from shapely.ops import unary_union + except ImportError as exc: # pragma: no cover - environment dependent + raise RuntimeError( + "Shapely >=2.0 is required to apply the EN4 source-region GeoJSON" + ) from exc + + with geojson_path.open("r", encoding="utf-8") as handle: + document = json.load(handle) + kind = document.get("type") + if kind == "FeatureCollection": + geometries = [ + shape(feature["geometry"]) + for feature in document.get("features", []) + if feature.get("geometry") is not None + ] + geometry = unary_union(geometries) if geometries else None + elif kind == "Feature": + raw_geometry = document.get("geometry") + geometry = shape(raw_geometry) if raw_geometry is not None else None + else: + geometry = shape(document) + if geometry is None or geometry.is_empty: + raise ValueError( + f"GeoJSON contains no usable geometry: {geojson_path}" + ) + if geometry.geom_type not in {"Polygon", "MultiPolygon"}: + raise ValueError( + "EN4 source-region GeoJSON must contain polygonal geometry; got " + f"{geometry.geom_type}" + ) + if not geometry.is_valid: + raise ValueError( + f"EN4 source-region GeoJSON geometry is invalid: {geojson_path}" + ) + + # GeoJSON uses longitude in the conventional -180..180 range, whereas + # EN4 longitude coordinates may use either that convention or 0..360. + longitude = ( + np.asarray(longitude_deg, dtype=float) + 180.0 + ) % 360.0 - 180.0 + latitude = np.asarray(latitude_deg, dtype=float) + return np.asarray( + covers(geometry, points(longitude, latitude)), dtype=bool + ) + + +def load_mesh_and_basins(cfg: Config) -> tuple[MeshData, np.ndarray]: + xr = require_xarray() + with xr.open_dataset( + cfg.mesh_file, decode_times=False, concat_characters=False + ) as ds: + bed_var = _find_variable(ds, "bedTopography") + thk_var = _find_variable(ds, "thickness") + bed = _array_with_nan( + bed_var.isel(Time=0).values if "Time" in bed_var.dims + else bed_var.values + ) + thickness = _array_with_nan( + thk_var.isel(Time=0).values if "Time" in thk_var.dims + else thk_var.values + ) + area = _array_with_nan(_find_variable(ds, "areaCell").values) + lat = _array_with_nan(_find_variable(ds, "latCell").values) + lon = _array_with_nan(_find_variable(ds, "lonCell").values) + lat_deg = radians_or_degrees_to_degrees(lat, "latitude") + lon_deg = radians_or_degrees_to_degrees(lon, "longitude") + + with xr.open_dataset( + cfg.region_mask_file, decode_times=False, concat_characters=False + ) as ds: + masks = np.asarray(_find_variable(ds, "regionCellMasks").values) + names = decode_region_names(_find_variable(ds, "regionNames").values) + if names != REGION_NAMES: + details = "\n".join( + f" {index + 1}: found={found!r}, expected={expected!r}" + for index, (found, expected) in enumerate(zip(names, REGION_NAMES)) + ) + raise ValueError( + f"Region names/order do not match the configured convention:\n" + f"{details}" + ) + basin_ids = build_basin_ids(masks) + + n_cells = bed.size + for name, values in ( + ("thickness", thickness), + ("areaCell", area), + ("latCell", lat), + ("lonCell", lon), + ): + if values.size != n_cells: + raise ValueError( + f"Mesh variable {name} has {values.size} cells, expected " + f"{n_cells}" + ) + return MeshData(bed, thickness, area, lat_deg, lon_deg), basin_ids + + +def _prepare_en4_mapping( + dataset: Any, + mesh: MeshData, + basin_ids: np.ndarray, + cfg: Config, +) -> dict[str, np.ndarray]: + lat_var = _find_variable(dataset, "lat", ("latitude",)) + lon_var = _find_variable(dataset, "lon", ("longitude",)) + source_lat = _array_with_nan(lat_var.values) + source_lon = _array_with_nan(lon_var.values) + if source_lat.ndim != 1 or source_lon.ndim != 1: + raise ValueError( + "EN4 latitude and longitude coordinates must be one-dimensional" + ) + lon_grid, lat_grid = np.meshgrid(source_lon, source_lat) + lat_flat = lat_grid.ravel() + lon_flat = lon_grid.ravel() + display_lon_flat = (lon_flat + 180.0) % 360.0 - 180.0 + source_region_ok = points_in_geojson( + display_lon_flat, lat_flat, cfg.en4_source_region_geojson + ) + latitude_ok = ( + (lat_flat >= cfg.en4_latitude_min) & (lat_flat <= cfg.en4_latitude_max) + ) + candidate_flat = np.flatnonzero(latitude_ok & source_region_ok) + + tree = cKDTree(unit_sphere_xyz(mesh.lat_deg, mesh.lon_deg)) + chord, nearest = tree.query( + unit_sphere_xyz(lat_flat[candidate_flat], lon_flat[candidate_flat]), + k=1, + ) + distance_km = chord_to_arc_km(chord) + keep = distance_km <= cfg.en4_max_mesh_distance_km + selected_flat = candidate_flat[keep] + if selected_flat.size == 0: + raise ValueError( + "No EN4 grid points inside the source-region GeoJSON fall within " + "max_mesh_distance_km of the MALI mesh" + ) + nearest_cells = nearest[keep] + return { + "flat_indices": selected_flat, + "lat": lat_flat[selected_flat], + "lon": display_lon_flat[selected_flat], + "basin": basin_ids[nearest_cells], + "distance_km": distance_km[keep], + "area_weight": np.clip( + np.cos(np.deg2rad(lat_flat[selected_flat])), 0.0, None + ), + "nlat": np.asarray([source_lat.size]), + "nlon": np.asarray([source_lon.size]), + } + + +def build_regional_profiles( + cfg: Config, mesh: MeshData, basin_ids: np.ndarray, logger +) -> RegionalProfiles: + xr = require_xarray() + files = discover_en4_files(cfg) + logger.info( + f"Found {len(files)} EN4 monthly analyses for regional profiles" + ) + + mapping: dict[str, np.ndarray] | None = None + source_z: np.ndarray | None = None + monthly_temp: list[np.ndarray] = [] + monthly_salinity: list[np.ndarray] = [] + monthly_tf: list[np.ndarray] = [] + monthly_counts: list[np.ndarray] = [] + monthly_temp_influence: list[np.ndarray] = [] + monthly_sal_influence: list[np.ndarray] = [] + dates: list[str] = [] + + for index, path in enumerate(files): + year, month = parse_yyyymm_from_name(path) + with xr.open_dataset( + path, + decode_times=False, + mask_and_scale=True, + concat_characters=False, + ) as ds: + depth_var = _find_variable(ds, "depth") + lat_var = _find_variable(ds, "lat", ("latitude",)) + lon_var = _find_variable(ds, "lon", ("longitude",)) + depth = _array_with_nan(depth_var.values) + z = -np.abs(depth) + if source_z is None: + source_z = z + mapping = _prepare_en4_mapping(ds, mesh, basin_ids, cfg) + elif (source_z.shape != z.shape or + not np.allclose(source_z, z, equal_nan=True)): + raise ValueError(f"EN4 depth coordinate changed in {path}") + + assert mapping is not None + temp_var = _find_variable(ds, "temperature") + sal_var = _find_variable(ds, "salinity") + temperature = _read_3d_en4_variable( + temp_var, depth_var.name, lat_var.name, lon_var.name + ) + salinity = _read_3d_en4_variable( + sal_var, depth_var.name, lat_var.name, lon_var.name + ) + temperature = _temperature_to_deg_c( + temperature, str(temp_var.attrs.get("units", "")) + ) + if "temperature_observation_weights" in ds: + temp_influence = _read_3d_en4_variable( + ds["temperature_observation_weights"], + depth_var.name, + lat_var.name, + lon_var.name, + ) + else: + temp_influence = np.full_like(temperature, np.nan) + if "salinity_observation_weights" in ds: + sal_influence = _read_3d_en4_variable( + ds["salinity_observation_weights"], + depth_var.name, + lat_var.name, + lon_var.name, + ) + else: + sal_influence = np.full_like(salinity, np.nan) + + n_depth = z.size + flat_indices = mapping["flat_indices"] + temp_selected = temperature.reshape(n_depth, -1)[:, flat_indices] + sal_selected = salinity.reshape(n_depth, -1)[:, flat_indices] + temp_influence_selected = temp_influence.reshape( + n_depth, -1 + )[:, flat_indices] + sal_influence_selected = sal_influence.reshape( + n_depth, -1 + )[:, flat_indices] + tf_selected = temp_selected - freezing_temperature( + sal_selected, + z[:, None], + cfg.freezing_a, + cfg.freezing_b, + cfg.freezing_c, + ) + weights = mapping["area_weight"] + basins = mapping["basin"] + + region_temp = np.full((len(REGION_NAMES), n_depth), np.nan) + region_sal = np.full_like(region_temp, np.nan) + region_tf = np.full_like(region_temp, np.nan) + region_count = np.zeros_like(region_temp, dtype=np.int32) + region_temp_influence = np.full_like(region_temp, np.nan) + region_sal_influence = np.full_like(region_temp, np.nan) + for region in range(len(REGION_NAMES)): + regional = basins == region + 1 + for level in range(n_depth): + good = ( + regional & + np.isfinite(temp_selected[level]) & + np.isfinite(sal_selected[level]) & + np.isfinite(tf_selected[level]) + ) + region_count[region, level] = int(good.sum()) + if np.any(good): + w = weights[good] + region_temp[region, level] = np.average( + temp_selected[level, good], weights=w + ) + region_sal[region, level] = np.average( + sal_selected[level, good], weights=w + ) + region_tf[region, level] = np.average( + tf_selected[level, good], weights=w + ) + temp_good = ( + regional & np.isfinite(temp_influence_selected[level]) + ) + if np.any(temp_good): + region_temp_influence[region, level] = np.average( + temp_influence_selected[level, temp_good], + weights=weights[temp_good], + ) + sal_good = ( + regional & np.isfinite(sal_influence_selected[level]) + ) + if np.any(sal_good): + region_sal_influence[region, level] = np.average( + sal_influence_selected[level, sal_good], + weights=weights[sal_good], + ) + + monthly_temp.append(region_temp) + monthly_salinity.append(region_sal) + monthly_tf.append(region_tf) + monthly_counts.append(region_count) + monthly_temp_influence.append(region_temp_influence) + monthly_sal_influence.append(region_sal_influence) + dates.append(f"{year:04d}-{month:02d}") + if index == 0 or (index + 1) % 12 == 0 or index + 1 == len(files): + logger.info( + f" processed EN4 month {index + 1}/{len(files)}: {dates[-1]}" + ) + + assert source_z is not None and mapping is not None + month_temp_array = np.asarray(monthly_temp) + month_sal_array = np.asarray(monthly_salinity) + month_tf_array = np.asarray(monthly_tf) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) + mean_temp = np.nanmean(month_temp_array, axis=0) + mean_sal = np.nanmean(month_sal_array, axis=0) + mean_tf = np.nanmean(month_tf_array, axis=0) + mean_freeze = freezing_temperature( + mean_sal, source_z[None, :], cfg.freezing_a, cfg.freezing_b, + cfg.freezing_c + ) + output_tf = np.vstack( + [ + interpolate_profile(source_z, mean_tf[region], cfg.ocean_levels_m) + for region in range(len(REGION_NAMES)) + ] + ) + + return RegionalProfiles( + source_z_m=source_z, + output_z_m=cfg.ocean_levels_m, + monthly_dates=dates, + monthly_temperature_degC=month_temp_array, + monthly_salinity=month_sal_array, + monthly_thermal_forcing_degC=month_tf_array, + temperature_degC=mean_temp, + salinity=mean_sal, + freezing_temperature_degC=mean_freeze, + thermal_forcing_degC=mean_tf, + output_thermal_forcing_degC=output_tf, + valid_gridpoint_counts=np.asarray(monthly_counts), + temperature_observation_influence=np.asarray(monthly_temp_influence), + salinity_observation_influence=np.asarray(monthly_sal_influence), + mapped_lats=mapping["lat"], + mapped_lons=mapping["lon"], + mapped_basins=mapping["basin"], + mapped_distances_km=mapping["distance_km"], + ) + + +def forcing_times(dataset: Any) -> list[str]: + if "xtime" not in dataset.variables: + raise KeyError("2-D forcing file must contain xtime(Time, StrLen)") + return decode_char_rows(dataset["xtime"].values) + + +def year_from_xtime(value: str) -> int: + match = re.match(r"\s*(\d{4})", value) + if match is None: + raise ValueError(f"Could not parse year from xtime value {value!r}") + return int(match.group(1)) + + +def calibration_time_indices( + times: Sequence[str], start_year: int, end_year: int +) -> np.ndarray: + years = np.asarray([year_from_xtime(value) for value in times]) + result = np.flatnonzero((years >= start_year) & (years <= end_year)) + if result.size == 0: + raise ValueError( + f"No forcing records fall within calibration years " + f"{start_year}-{end_year}" + ) + return result + + +def validate_forcing_schema(dataset: Any, cfg: Config, n_cells: int) -> Any: + if cfg.forcing_variable not in dataset.variables: + raise KeyError( + f"2-D forcing variable {cfg.forcing_variable!r} is missing" + ) + variable = dataset[cfg.forcing_variable] + if variable.ndim != 2 or variable.shape[1] != n_cells: + raise ValueError( + f"{cfg.forcing_variable} must have shape " + f"(Time, nCells={n_cells}); got {variable.shape}" + ) + return variable + + +def calibrate_regional_delta_t( + cfg: Config, + mesh: MeshData, + basin_ids: np.ndarray, + profiles: RegionalProfiles, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + xr = require_xarray() + floating, draft = floating_mask_and_draft( + mesh.bed, + mesh.thickness, + cfg.rho_ice, + cfg.rho_seawater, + cfg.flotation_tolerance_m, + cfg.minimum_ice_thickness_m, + ) + anchor = np.clip(mesh.bed, -cfg.source_max_depth_m, 0.0) + base_at_anchor = profiles_at_cell_depths( + profiles.output_thermal_forcing_degC, cfg.ocean_levels_m, basin_ids, + anchor + ) + base_at_draft = profiles_at_cell_depths( + profiles.output_thermal_forcing_degC, cfg.ocean_levels_m, basin_ids, + draft + ) + + floating_counts = np.asarray( + [ + np.count_nonzero(floating & (basin_ids == region + 1)) + for region in range(len(REGION_NAMES)) + ] + ) + empty = np.flatnonzero(floating_counts == 0) + if empty.size: + names = ", ".join(REGION_KEYS[index] for index in empty) + raise ValueError( + "Cannot calibrate regional deltaT because the initial geometry " + f"has no floating cells in: {names}" + ) + + with xr.open_dataset( + cfg.forcing_2d_file, + decode_times=False, + mask_and_scale=True, + concat_characters=False, + ) as ds: + forcing_var = validate_forcing_schema(ds, cfg, mesh.bed.size) + times = forcing_times(ds) + indices = calibration_time_indices( + times, cfg.calibration_start_year, cfg.calibration_end_year + ) + monthly_means = np.full((indices.size, len(REGION_NAMES)), np.nan) + for output_index, time_index in enumerate(indices): + forcing = _array_with_nan( + forcing_var.isel(Time=int(time_index)).values + ) + offset = forcing - base_at_anchor + tf_draft = base_at_draft + offset + for region in range(len(REGION_NAMES)): + cells = ( + floating & + (basin_ids == region + 1) & + np.isfinite(tf_draft) + ) + if not np.any(cells): + raise ValueError( + f"No finite draft thermal forcing for " + f"{REGION_KEYS[region]} at {times[time_index]}" + ) + monthly_means[output_index, region] = np.average( + tf_draft[cells], weights=mesh.area[cells] + ) + + coefficient = cfg.rho_seawater * cfg.cp_seawater / ( + cfg.rho_ice * cfg.latent_heat_ice + ) + delta_t = np.asarray( + [ + calibrate_delta_t( + monthly_means[:, region], + cfg.regional_melt_targets_m_per_yr[region], + cfg.gamma0_m_per_yr, + coefficient, + ) + for region in range(len(REGION_NAMES)) + ] + ) + achieved = np.asarray( + [ + nonlocal_mean_melt( + delta_t[region], monthly_means[:, region], cfg.gamma0_m_per_yr, + coefficient + ) + for region in range(len(REGION_NAMES)) + ] + ) + return delta_t, achieved, monthly_means + + +def write_output( + cfg: Config, + mesh: MeshData, + basin_ids: np.ndarray, + profiles: RegionalProfiles, + regional_delta_t: np.ndarray, + logger, +) -> None: + xr = require_xarray() + try: + import dask + except ImportError as exc: # pragma: no cover - environment dependent + raise RuntimeError( + "Dask is required to construct and stream the multi-gigabyte " + "xarray output" + ) from exc + output = cfg.output_file + output.parent.mkdir(parents=True, exist_ok=True) + if output.exists() and not cfg.overwrite: + raise FileExistsError(f"Output exists and overwrite=false: {output}") + temporary = output.with_name(output.name + ".partial") + if temporary.exists(): + raise FileExistsError( + f"Partial output already exists: {temporary}. Remove or rename it " + f"after inspecting it." + ) + + anchor = np.clip(mesh.bed, -cfg.source_max_depth_m, 0.0) + base_at_anchor = profiles_at_cell_depths( + profiles.output_thermal_forcing_degC, cfg.ocean_levels_m, basin_ids, + anchor + ) + base_by_cell = profiles.output_thermal_forcing_degC[basin_ids - 1, :] + delta_t_by_cell = regional_delta_t[basin_ids - 1] + + try: + with xr.open_dataset( + cfg.forcing_2d_file, + decode_times=False, + mask_and_scale=True, + concat_characters=False, + chunks={"Time": 1}, + ) as source: + forcing_var = validate_forcing_schema(source, cfg, mesh.bed.size) + times = forcing_times(source) + # Fail before initiating a multi-gigabyte write if any source value + # is absent. This reduction stays lazy until compute() and does not + # load the full forcing array into memory. + invalid_count = int( + (~np.isfinite(forcing_var)).sum().compute().item() + ) + if invalid_count: + raise ValueError( + f"2-D forcing contains {invalid_count} invalid values; " + "explicit missing-value handling is required before " + "building MALI forcing" + ) + + cell_coord = np.arange(mesh.bed.size, dtype=np.int32) + layer_coord = np.arange(cfg.ocean_levels_m.size, dtype=np.int32) + base_cells = xr.DataArray( + base_by_cell.astype(np.float32), + dims=("nCells", "nISMIP6OceanLayers"), + coords={ + "nCells": cell_coord, + "nISMIP6OceanLayers": layer_coord, + }, + ) + anchor_cells = xr.DataArray( + base_at_anchor.astype(np.float32), + dims=("nCells",), + coords={"nCells": cell_coord}, + ) + forcing_3d = ( + base_cells + (forcing_var.astype(np.float32) - anchor_cells) + ).transpose("Time", "nCells", "nISMIP6OceanLayers").assign_attrs( + units="C", + long_name="3D thermal forcing for nonlocal ISMIP6 ice-shelf " + "melt method", + ) + + target = xr.Dataset( + data_vars={ + "xtime": source["xtime"], + "ismip6shelfMelt_basin": xr.DataArray( + basin_ids.astype(np.int32), + dims=("nCells",), + attrs={ + "description": "One-based basin number for " + "regional ISMIP6 shelf-melt forcing" + }, + ), + "ismip6shelfMelt_gamma0": xr.DataArray( + np.float32(cfg.gamma0_m_per_yr), + attrs={ + "units": "m yr^-1", + "description": "Uniform gamma0 for nonlocal " + "Jourdain et al. (2020) shelf melt", + }, + ), + "ismip6shelfMelt_deltaT": xr.DataArray( + delta_t_by_cell.astype(np.float32), + dims=("nCells",), + attrs={ + "units": "K", + "description": "Regionally calibrated, cellwise " + "temperature-bias correction", + }, + ), + "ismip6shelfMelt_zOcean": xr.DataArray( + cfg.ocean_levels_m.astype(np.float32), + dims=("nISMIP6OceanLayers",), + attrs={"units": "m", "positive": "up"}, + ), + "ismip6shelfMelt_3dThermalForcing": forcing_3d, + }, + attrs={ + "title": "Regional three-dimensional Greenland ocean " + "thermal forcing for MALI", + "source_2d_forcing": str(cfg.forcing_2d_file), + "source_en4_version": cfg.en4_version, + "source_en4_bias_correction": cfg.en4_bias_correction, + "source_en4_region_geojson": str( + cfg.en4_source_region_geojson + ), + "en4_profile_period": f"{cfg.profile_start_year}-" + f"{cfg.profile_end_year}", + "deltaT_calibration_period": + f"{cfg.calibration_start_year}" + f"-{cfg.calibration_end_year}", + "source_ocean_max_depth_m": cfg.source_max_depth_m, + "region_names": " | ".join(REGION_NAMES), + "regional_deltaT_K": ", ".join( + f"{value:.8g}" for value in regional_delta_t + ), + "history": datetime.now(timezone.utc).strftime( + "%Y-%m-%dT%H:%M:%SZ" + ) + ": " + " ".join(sys.argv), + }, + ) + # Dimension coordinates are implementation details, not MALI input + # fields. Drop them while retaining the named dimensions. + target = target.drop_vars( + [ + name + for name in ("nCells", "nISMIP6OceanLayers") + if name in target.coords + ] + ) + encoding = { + "ismip6shelfMelt_basin": { + "dtype": "int32", "_FillValue": None + }, + "ismip6shelfMelt_gamma0": { + "dtype": "float32", "_FillValue": None + }, + "ismip6shelfMelt_deltaT": { + "dtype": "float32", "_FillValue": None + }, + "ismip6shelfMelt_zOcean": { + "dtype": "float32", "_FillValue": None + }, + "ismip6shelfMelt_3dThermalForcing": { + "dtype": "float32", "_FillValue": None + }, + } + logger.info( + f"Writing {len(times)} monthly records to {temporary} with " + "xarray (NETCDF3_64BIT, float32)" + ) + with dask.config.set(scheduler="single-threaded"): + target.to_netcdf( + temporary, + engine="scipy", + format="NETCDF3_64BIT", + unlimited_dims=["Time"], + encoding=encoding, + ) + target.close() + os.replace(temporary, output) + except Exception: + logger.warning( + f"Output was not finalized; partial file, if any, is at " + f"{temporary}" + ) + raise + + +def write_diagnostics( + cfg: Config, + mesh: MeshData, + basin_ids: np.ndarray, + profiles: RegionalProfiles, + regional_delta_t: np.ndarray, + achieved_melt: np.ndarray, + calibration_monthly_tf: np.ndarray, +) -> None: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + directory = cfg.diagnostics_directory + directory.mkdir(parents=True, exist_ok=True) + + with (directory / "regional_profiles.csv").open( + "w", newline="", encoding="utf-8" + ) as handle: + writer = csv.writer(handle) + writer.writerow( + [ + "region_id", + "region_key", + "region_name", + "z_m", + "temperature_degC", + "salinity", + "freezing_temperature_degC", + "thermal_forcing_degC", + ] + ) + for region in range(len(REGION_NAMES)): + for level, z in enumerate(profiles.source_z_m): + writer.writerow( + [ + region + 1, + REGION_KEYS[region], + REGION_NAMES[region], + float(z), + float(profiles.temperature_degC[region, level]), + float(profiles.salinity[region, level]), + float( + profiles.freezing_temperature_degC[region, level] + ), + float(profiles.thermal_forcing_degC[region, level]), + ] + ) + + with (directory / "en4_coverage.csv").open( + "w", newline="", encoding="utf-8" + ) as handle: + writer = csv.writer(handle) + writer.writerow( + [ + "date", + "region_id", + "region_key", + "z_m", + "valid_gridpoint_count", + "temperature_observation_influence", + "salinity_observation_influence", + ] + ) + for month, date in enumerate(profiles.monthly_dates): + for region in range(len(REGION_NAMES)): + for level, z in enumerate(profiles.source_z_m): + writer.writerow( + [ + date, + region + 1, + REGION_KEYS[region], + float(z), + int( + profiles.valid_gridpoint_counts[ + month, region, level + ] + ), + float( + profiles.temperature_observation_influence[ + month, region, level + ] + ), + float( + profiles.salinity_observation_influence[ + month, region, level + ] + ), + ] + ) + + floating, _ = floating_mask_and_draft( + mesh.bed, + mesh.thickness, + cfg.rho_ice, + cfg.rho_seawater, + cfg.flotation_tolerance_m, + cfg.minimum_ice_thickness_m, + ) + floating_counts = [ + int(np.count_nonzero(floating & (basin_ids == region + 1))) + for region in range(len(REGION_NAMES)) + ] + floating_areas_km2 = [ + float(np.sum(mesh.area[floating & (basin_ids == region + 1)]) / 1.0e6) + for region in range(len(REGION_NAMES)) + ] + + calibration_summary = { + "gamma0_m_per_yr": cfg.gamma0_m_per_yr, + "calibration_period": [ + cfg.calibration_start_year, + cfg.calibration_end_year, + ], + "profile_period": [cfg.profile_start_year, cfg.profile_end_year], + "regions": [ + { + "id": region + 1, + "key": REGION_KEYS[region], + "name": REGION_NAMES[region], + "target_melt_m_per_yr": float( + cfg.regional_melt_targets_m_per_yr[region] + ), + "achieved_melt_m_per_yr": float(achieved_melt[region]), + "deltaT_K": float(regional_delta_t[region]), + "mean_calibration_TF_degC": float( + np.mean(calibration_monthly_tf[:, region]) + ), + "minimum_calibration_TF_degC": float( + np.min(calibration_monthly_tf[:, region]) + ), + "maximum_calibration_TF_degC": float( + np.max(calibration_monthly_tf[:, region]) + ), + "floating_cell_count": floating_counts[region], + "floating_area_km2": floating_areas_km2[region], + } + for region in range(len(REGION_NAMES)) + ], + } + with (directory / "deltaT_calibration.json").open( + "w", encoding="utf-8" + ) as handle: + json.dump(calibration_summary, handle, indent=2) + handle.write("\n") + + colors = plt.get_cmap("tab10")(np.arange(len(REGION_NAMES))) + fig, ax = plt.subplots(figsize=(8, 8)) + mesh_stride = max(1, mesh.lat_deg.size // 100_000) + ax.scatter( + mesh.lon_deg[::mesh_stride], + mesh.lat_deg[::mesh_stride], + c=colors[basin_ids[::mesh_stride] - 1], + s=0.15, + alpha=0.2, + linewidths=0, + ) + ax.scatter( + profiles.mapped_lons, + profiles.mapped_lats, + c=colors[profiles.mapped_basins - 1], + s=5, + edgecolors="none", + ) + ax.set_xlabel("Longitude (degrees east)") + ax.set_ylabel("Latitude (degrees north)") + ax.set_title("EN4 grid points assigned to Greenland shelf regions") + ax.grid(alpha=0.25) + fig.tight_layout() + fig.savefig(directory / "en4_region_assignment.png", dpi=180) + plt.close(fig) + + # Validate the profile translation at the effective seafloor and map one + # representative monthly field at all output depths. + xr = require_xarray() + with xr.open_dataset( + cfg.forcing_2d_file, + decode_times=False, + mask_and_scale=True, + concat_characters=False, + ) as source: + times = forcing_times(source) + indices = calibration_time_indices( + times, cfg.calibration_start_year, cfg.calibration_end_year + ) + time_index = int(indices[0]) + forcing_2d = _array_with_nan( + validate_forcing_schema(source, cfg, mesh.bed.size) + .isel(Time=time_index) + .values + ) + anchor = np.clip(mesh.bed, -cfg.source_max_depth_m, 0.0) + base_at_anchor = profiles_at_cell_depths( + profiles.output_thermal_forcing_degC, + cfg.ocean_levels_m, + basin_ids, + anchor, + ) + base_by_cell = profiles.output_thermal_forcing_degC[basin_ids - 1] + offset = forcing_2d - base_at_anchor + forcing_3d = base_by_cell + offset[:, None] + reconstructed_anchor = base_at_anchor + offset + max_anchor_error = float( + np.nanmax(np.abs(reconstructed_anchor - forcing_2d)) + ) + with (directory / "forcing_validation.json").open( + "w", encoding="utf-8" + ) as handle: + json.dump( + { + "representative_time": times[time_index], + "maximum_absolute_anchor_error_degC": max_anchor_error, + "effective_anchor_depth_range_m": [ + float(np.nanmin(anchor)), + float(np.nanmax(anchor)), + ], + "source_max_depth_m": cfg.source_max_depth_m, + }, + handle, + indent=2, + ) + handle.write("\n") + + n_levels = cfg.ocean_levels_m.size + plotted_levels = np.unique( + np.rint(np.linspace(0, n_levels - 1, min(n_levels, 6))).astype(int) + ) + ncols = min(2, plotted_levels.size) + nrows = int(math.ceil(plotted_levels.size / ncols)) + fig, axes = plt.subplots( + nrows, ncols, figsize=(7 * ncols, 6 * nrows), squeeze=False + ) + plot_stride = max(1, mesh.lat_deg.size // 150_000) + color_limits = np.nanpercentile(forcing_3d[::plot_stride], [2.0, 98.0]) + for plot_index, ax in enumerate(axes.ravel()): + if plot_index >= plotted_levels.size: + ax.set_visible(False) + continue + level = int(plotted_levels[plot_index]) + scatter = ax.scatter( + mesh.lon_deg[::plot_stride], + mesh.lat_deg[::plot_stride], + c=forcing_3d[::plot_stride, level], + s=0.5, + linewidths=0, + cmap="coolwarm", + vmin=color_limits[0], + vmax=color_limits[1], + ) + ax.set_title(f"z = {cfg.ocean_levels_m[level]:g} m") + ax.set_xlabel("Longitude") + ax.set_ylabel("Latitude") + fig.colorbar(scatter, ax=ax, label="Thermal forcing (°C)") + fig.suptitle(f"Translated 3-D forcing: {times[time_index]}") + fig.tight_layout() + fig.savefig( + directory / "thermal_forcing_at_representative_ocean_levels.png", + dpi=180, + ) + plt.close(fig) + + for region in range(len(REGION_NAMES)): + fig, axes = plt.subplots(1, 2, figsize=(10, 7), sharey=True) + for monthly in profiles.monthly_temperature_degC[:, region, :]: + axes[0].plot( + monthly, profiles.source_z_m, color="tab:blue", alpha=0.08, + linewidth=0.5 + ) + axes[0].plot( + profiles.temperature_degC[region], + profiles.source_z_m, + color="black", + linewidth=2, + label="Monthly mean climatology", + ) + selected_temp = interpolate_profile( + profiles.source_z_m, profiles.temperature_degC[region], + profiles.output_z_m + ) + axes[0].scatter( + selected_temp, profiles.output_z_m, color="black", marker="*", + zorder=3, label="MALI levels" + ) + axes[0].set_xlabel("Potential temperature (°C)") + axes[0].set_ylabel("Elevation (m)") + axes[0].legend(loc="best", fontsize=8) + + for monthly in profiles.monthly_thermal_forcing_degC[:, region, :]: + axes[1].plot( + monthly, profiles.source_z_m, color="tab:red", alpha=0.08, + linewidth=0.5 + ) + axes[1].plot( + profiles.thermal_forcing_degC[region], + profiles.source_z_m, + color="black", + linewidth=2, + ) + axes[1].scatter( + profiles.output_thermal_forcing_degC[region], + profiles.output_z_m, + color="black", + marker="*", + zorder=3, + ) + axes[1].set_xlabel("Thermal forcing (°C)") + for ax in axes: + ax.grid(alpha=0.3) + ax.set_ylim( + min(-cfg.source_max_depth_m, profiles.output_z_m[-1]), 0.0 + ) + fig.suptitle(REGION_NAMES[region]) + fig.tight_layout() + fig.savefig( + directory / f"profile_{region + 1:02d}_{REGION_KEYS[region]}.png", + dpi=180, + ) + plt.close(fig) + + +def print_summary( + cfg: Config, + regional_delta_t: np.ndarray, + achieved_melt: np.ndarray, + logger, +) -> None: + logger.info("Regional deltaT calibration") + logger.info("region target_m/yr achieved_m/yr deltaT_K") + for region, key in enumerate(REGION_KEYS): + logger.info( + f"{key:15s} {cfg.regional_melt_targets_m_per_yr[region]:11.5f} " + f"{achieved_melt[region]:14.5f} {regional_delta_t[region]:9.5f}" + ) + + +def run(cfg: Config, logger, prepare_only: bool = False) -> None: + mesh, basin_ids = load_mesh_and_basins(cfg) + profiles = build_regional_profiles(cfg, mesh, basin_ids, logger) + delta_t, achieved, calibration_monthly_tf = calibrate_regional_delta_t( + cfg, mesh, basin_ids, profiles + ) + print_summary(cfg, delta_t, achieved, logger) + write_diagnostics( + cfg, mesh, basin_ids, profiles, delta_t, achieved, + calibration_monthly_tf + ) + if not prepare_only: + write_output(cfg, mesh, basin_ids, profiles, delta_t, logger) + logger.info(f"Created {cfg.output_file}") + else: + logger.info( + "Preparation-only run complete; the multi-gigabyte forcing file " + "was not written" + ) + logger.info(f"Diagnostics: {cfg.diagnostics_directory}") diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py index c0eed44c97..7180f8f603 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -171,7 +171,8 @@ def _run_scenario(self): # Combine remapped files and rename to MALI conventions logger.info("Combining remapped files and renaming variables...") - output_file = (f"{mali_mesh_name}_thermal_forcing_{model}_{scenario}_" + tf_label = "3dThermalForcing" if ocean_3d else "2dThermalForcing" + output_file = (f"{mali_mesh_name}_{tf_label}_{model}_{scenario}_" f"{start_year}-{end_year}.nc") if ocean_3d: diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/set_up_experiment.py b/compass/landice/tests/ismip7_run/ismip7_ais/set_up_experiment.py index 98d59d4c40..f319f394a5 100644 --- a/compass/landice/tests/ismip7_run/ismip7_ais/set_up_experiment.py +++ b/compass/landice/tests/ismip7_run/ismip7_ais/set_up_experiment.py @@ -219,7 +219,7 @@ def setup(self): # noqa: C901 os.path.join(self.work_dir, temp_grad_fname)) # Thermal forcing - tf_search = os.path.join(ocean_dir, '*thermal_forcing_*.nc') + tf_search = os.path.join(ocean_dir, '*3dThermalForcing_*.nc') tf_list = glob.glob(tf_search) if len(tf_list) == 1: tf_fname = os.path.split(tf_list[0])[-1] diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py b/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py index 4118ce44b5..1f1e103872 100644 --- a/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py +++ b/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py @@ -201,7 +201,7 @@ def setup(self): # noqa: C901 os.path.join(self.work_dir, temp_grad_fname)) # GrIS uses 2D thermal forcing - tf_search = os.path.join(ocean_dir, '*thermal_forcing_*.nc') + tf_search = os.path.join(ocean_dir, '*2dThermalForcing_*.nc') tf_list = glob.glob(tf_search) if len(tf_list) == 1: tf_fname = os.path.split(tf_list[0])[-1] diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index 75bb6bab8c..6b7f5d909e 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -410,6 +410,10 @@ ismip7_forcing ocean_thermal.process_thermal_forcing.ProcessThermalForcing ocean_thermal.process_thermal_forcing.ProcessThermalForcing.setup ocean_thermal.process_thermal_forcing.ProcessThermalForcing.run + ocean_thermal.build_3d_thermal_forcing.BuildGreenland3dThermalForcing + ocean_thermal.build_3d_thermal_forcing.BuildGreenland3dThermalForcing.run + ocean_thermal.greenland_3d.Config + ocean_thermal.greenland_3d.run fracture.Fracture fracture.Fracture.configure diff --git a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst index 0a37029290..b3e7af7c58 100644 --- a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst @@ -98,10 +98,12 @@ ocean_thermal ~~~~~~~~~~~~~ The :py:class:`compass.landice.tests.ismip7_forcing.ocean_thermal.OceanThermal` -test case processes the ISMIP7 ocean thermal forcing. It contains a single step, +test case processes the ISMIP7 ocean thermal forcing. It contains two steps, :py:class:`~compass.landice.tests.ismip7_forcing.ocean_thermal.process_thermal_forcing.ProcessThermalForcing`, which handles both AIS (3D, decade-spanning files) and GrIS (2D, yearly files) -by branching on the ``ocean_3d`` parameter from ``ice_sheet_params``. +by branching on the ``ocean_3d`` parameter from ``ice_sheet_params``, and +:py:class:`~compass.landice.tests.ismip7_forcing.ocean_thermal.build_3d_thermal_forcing.BuildGreenland3dThermalForcing`, +which optionally builds a 3D GrIS field from the 2D forcing (see below). The ``run()`` method dispatches to two sub-methods based on the boolean config options ``process_ocean_thermal`` and ``process_ocean_climatology`` in the @@ -132,6 +134,29 @@ For GrIS, the step: * Remaps 2D monthly thermal forcing * Produces ``ismip6_2dThermalForcing`` (dims: Time × nCells) +build_3d_thermal_forcing (GrIS 3-D) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The +:py:class:`~compass.landice.tests.ismip7_forcing.ocean_thermal.build_3d_thermal_forcing.BuildGreenland3dThermalForcing` +step (GrIS only, gated by ``process_ocean_thermal_3d``) converts the GrIS 2D +thermal forcing into a 30-level 3D field for MALI's nonlocal (Jourdain et al. +2020) melt scheme. Its ``run()`` reconstructs the 2D output path that +``ProcessThermalForcing._run_scenario`` wrote (mirroring the ``forcing_group`` +and ocean-source logic), injects the compass-derived mesh / 2D-forcing / +output / diagnostics paths into a +:py:class:`~compass.landice.tests.ismip7_forcing.ocean_thermal.greenland_3d.Config` +built from the JSON ``config_file``, and calls +:py:func:`~compass.landice.tests.ismip7_forcing.ocean_thermal.greenland_3d.run`. +The ported science lives in +:py:mod:`compass.landice.tests.ismip7_forcing.ocean_thermal.greenland_3d`: +seven regional EN4 profiles, seafloor anchoring to the 2D forcing, and +per-region ``deltaT`` calibration. The output supplements the 2D file with +``ismip6shelfMelt_3dThermalForcing``, ``ismip6shelfMelt_deltaT``, +``ismip6shelfMelt_gamma0``, ``ismip6shelfMelt_zOcean``, and +``ismip6shelfMelt_basin``. The multi-gigabyte field is streamed record by +record (dask, ``scipy`` engine, ``NETCDF3_64BIT``, ``.partial``-then-rename). + .. _dev_landice_ismip7_forcing_fracture: fracture diff --git a/docs/users_guide/landice/test_groups/ismip7_forcing.rst b/docs/users_guide/landice/test_groups/ismip7_forcing.rst index 6afb6b7d92..447aef2154 100644 --- a/docs/users_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/users_guide/landice/test_groups/ismip7_forcing.rst @@ -17,11 +17,14 @@ and ``fracture``. ``process_smb``, ``process_temperature``, ``process_smb_gradient``, ``process_temperature_gradient``, and ``process_runoff``. -* The ``ocean_thermal`` test case has one step: ``process_thermal_forcing``. - For AIS this produces 3D thermal forcing (with 30 ocean depth layers); for - GrIS it produces 2D (depth-averaged) thermal forcing. The step can also - process the observational ocean thermal forcing climatology (Zhou et al.) - for AIS, controlled by the ``process_ocean_climatology`` config option. +* The ``ocean_thermal`` test case has two steps: ``process_thermal_forcing`` + and ``build_3d_thermal_forcing``. ``process_thermal_forcing`` produces, for + AIS, 3D thermal forcing (with 30 ocean depth layers) and, for GrIS, 2D + (depth-averaged) thermal forcing. It can also process the observational + ocean thermal forcing climatology (Zhou et al.) for AIS, controlled by the + ``process_ocean_climatology`` config option. ``build_3d_thermal_forcing`` + optionally converts the GrIS 2D forcing into a 3D field (GrIS only, + controlled by ``process_ocean_thermal_3d``). * The ``fracture`` test case has three steps: ``process_excess_melt`` (Path A), ``process_lake_properties`` (Path B), and @@ -90,7 +93,8 @@ Processed forcing is written under ``output_base_path`` in a layout that the {output_base_path}/{group}/atmosphere/{mesh}_SMB_{source}_{scenario}_{years}.nc {output_base_path}/{group}/atmosphere/{mesh}_temperature_{source}_{scenario}_{years}.nc {output_base_path}/{group}/atmosphere/{mesh}_runoff_... (and the two gradients) - {output_base_path}/{group}/ocean_thermal_forcing/{mesh}_thermal_forcing_{source}_{scenario}_{years}.nc + {output_base_path}/{group}/ocean_thermal_forcing/{mesh}_2dThermalForcing_{source}_{scenario}_{years}.nc (GrIS) + {output_base_path}/{group}/ocean_thermal_forcing/{mesh}_3dThermalForcing_{source}_{scenario}_{years}.nc (AIS; optional GrIS 3-D) The ``group`` directory is ``{model}_{scenario}`` for the ESM scenarios and ``{scenario}`` (i.e. ``OCX``) for OCX, whose atmosphere and ocean use @@ -214,6 +218,10 @@ values are: # Whether to process observational ocean thermal forcing climatology process_ocean_climatology = true + # Whether to build regional 3-D Greenland ocean thermal forcing from the + # 2-D forcing (GrIS only; Antarctica already produces 3-D forcing) + process_ocean_thermal_3d = false + # config options for ismip7 atmosphere forcing [ismip7_atmosphere] @@ -247,6 +255,14 @@ values are: # Base path to observational climatology data base_path_climatology = /path/to/ISMIP7/forcing/AIS/obs/zhou_annual_06_nov + # config options for building 3-D Greenland ocean thermal forcing + [ismip7_ocean_thermal_3d] + + # Path to the JSON config with the 3-D-specific parameters (EN4 directory, + # region-mask file, source-region GeoJSON, calibration, vertical grid). + # User must supply when process_ocean_thermal_3d is true. + config_file = NotAvailable + # config options for ismip7 fracture (Path C, ice shelf collapse) forcing [ismip7_fracture] @@ -343,6 +359,28 @@ For **GrIS**, thermal forcing is 2D (depth-averaged), with monthly temporal resolution and yearly input files. The output variable is ``ismip6_2dThermalForcing``. +build_3d_thermal_forcing (GrIS 3-D) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The optional ``build_3d_thermal_forcing`` step (GrIS only, gated by +``process_ocean_thermal_3d = true``) converts the GrIS 2D thermal forcing into +a 30-level 3D field for MALI's nonlocal (Jourdain et al. 2020) melt scheme. It +auto-chains from the 2D forcing the ``process_thermal_forcing`` step just +wrote, builds seven regional vertical profiles from monthly EN4 objective +analyses, anchors each cell's profile to the effective seafloor to match the +2D forcing, and calibrates one ``ismip6shelfMelt_deltaT`` per region while +holding ``ismip6shelfMelt_gamma0`` fixed. + +The 3-D-specific parameters (EN4 directory, region-mask file, source-region +GeoJSON, calibration targets, vertical grid, physical constants) are supplied +through a JSON config file referenced by the ``[ismip7_ocean_thermal_3d]`` +``config_file`` option; the mesh, 2D forcing, output, and diagnostics paths +are injected automatically. The output supplements (does not replace) the 2D +forcing and uses the variables ``ismip6shelfMelt_3dThermalForcing``, +``ismip6shelfMelt_deltaT``, ``ismip6shelfMelt_gamma0``, +``ismip6shelfMelt_zOcean``, and ``ismip6shelfMelt_basin``, matching the AIS +3D convention. + .. _landice_ismip7_forcing_fracture: fracture From 2362e649120674e0f5530557ea6f2a7be7e08ce5 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 4 Sep 2026 14:07:22 -0700 Subject: [PATCH 07/17] Add optional GrIS 3-D ocean thermal forcing ingestion to ismip7_run Add use_3d_thermal_forcing (default false) to [ismip7_run_gris]. When true, ismip7_gris looks for *3dThermalForcing_*.nc instead of *2dThermalForcing_*.nc, reads ismip6shelfMelt_3dThermalForcing + ismip6shelfMelt_zOcean at annual intervals (instead of ismip6_2dThermalForcing at monthly intervals), adds an ismip7_params stream reading ismip6shelfMelt_deltaT/_basin/_gamma0 from melt_params_path (mirroring the AIS convention, and consuming the output of the new build_3d_thermal_forcing step in landice/ismip7_forcing/ocean_thermal), and sets config_use_3d_thermal_forcing_for_face_melt = .true. melt_params_path is now validated as required when use_3d_thermal_forcing is true. streams.landice.template gates the TF stream's variable list and the new ismip7_params stream with Jinja2 conditionals on use_3d_thermal_forcing; verified both branches render as valid XML with the expected variables. Updates user and developer docs. --- .../ismip7_run/ismip7_gris/ismip7_gris.cfg | 6 ++++ .../ismip7_gris/ismip7_gris_test.cfg | 6 ++++ .../ismip7_gris/set_up_experiment.py | 26 +++++++++++++-- .../ismip7_gris/streams.landice.template | 19 ++++++++++- .../landice/test_groups/ismip7_run.rst | 15 +++++++-- .../landice/test_groups/ismip7_run.rst | 33 +++++++++++++++---- 6 files changed, 94 insertions(+), 11 deletions(-) diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris.cfg b/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris.cfg index a3b6a6229a..cc345c2a49 100644 --- a/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris.cfg +++ b/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris.cfg @@ -22,8 +22,14 @@ forcing_basepath = NotAvailable init_cond_path = NotAvailable # Path to the file for the basal melt parametrization coefficients. +# Required if use_3d_thermal_forcing is true. melt_params_path = NotAvailable +# Whether to force with 3D ocean thermal forcing (built by the +# build_3d_thermal_forcing step of landice/ismip7_forcing/ocean_thermal) +# instead of the default 2D forcing. Requires melt_params_path. +use_3d_thermal_forcing = false + # Path to the reference surface elevation file reference_surface_path = NotAvailable diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris_test.cfg b/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris_test.cfg index 9287dc1a1d..d43c460073 100644 --- a/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris_test.cfg +++ b/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris_test.cfg @@ -22,8 +22,14 @@ forcing_basepath = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/test_processing/ init_cond_path = /global/cfs/cdirs/fanssie/MALI_input_files/GIS_1to10km_r02/GIS_1to10km_r02_20230202_remove_icebergs.nc # Path to the file for the basal melt parametrization coefficients. +# Required if use_3d_thermal_forcing is true. melt_params_path = NotAvailable +# Whether to force with 3D ocean thermal forcing (built by the +# build_3d_thermal_forcing step of landice/ismip7_forcing/ocean_thermal) +# instead of the default 2D forcing. Requires melt_params_path. +use_3d_thermal_forcing = false + # Path to the reference surface elevation file reference_surface_path = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/forcing/GIS/upper_surface.nc diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py b/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py index 1f1e103872..8cdcb12caa 100644 --- a/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py +++ b/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py @@ -59,6 +59,11 @@ def setup(self): # noqa: C901 reference_surface_path = section.get('reference_surface_path') reference_surface_fname = os.path.split(reference_surface_path)[-1] calving_method = section.get('calving_method') + use_3d_thermal_forcing = section.getboolean('use_3d_thermal_forcing') + if use_3d_thermal_forcing and melt_params_path == 'NotAvailable': + raise ValueError( + "melt_params_path must be supplied when " + "use_3d_thermal_forcing is true") exp_info = self.exp_info scenario = exp_info['scenario'] @@ -200,8 +205,14 @@ def setup(self): # noqa: C901 os.symlink(temp_grad_list[0], os.path.join(self.work_dir, temp_grad_fname)) - # GrIS uses 2D thermal forcing - tf_search = os.path.join(ocean_dir, '*2dThermalForcing_*.nc') + # GrIS thermal forcing: 2D by default, or 3D when + # use_3d_thermal_forcing is true (see build_3d_thermal_forcing + # in landice/ismip7_forcing/ocean_thermal) + if use_3d_thermal_forcing: + tf_pattern = '*3dThermalForcing_*.nc' + else: + tf_pattern = '*2dThermalForcing_*.nc' + tf_search = os.path.join(ocean_dir, tf_pattern) tf_list = glob.glob(tf_search) if len(tf_list) == 1: tf_fname = os.path.split(tf_list[0])[-1] @@ -215,9 +226,13 @@ def setup(self): # noqa: C901 if scenario == 'ctrl': forcing_interval_monthly = 'initial_only' forcing_interval_annual = 'initial_only' + forcing_interval_TF = 'initial_only' else: forcing_interval_monthly = '0000-01-00_00:00:00' forcing_interval_annual = '0001-00-00_00:00:00' + forcing_interval_TF = (forcing_interval_annual + if use_3d_thermal_forcing + else forcing_interval_monthly) stream_replacements = { 'input_file_init_cond': init_cond_fname if is_historical @@ -234,6 +249,8 @@ def setup(self): # noqa: C901 'input_file_temperature_gradient_forcing': temp_grad_fname, 'forcing_interval_monthly': forcing_interval_monthly, 'forcing_interval_annual': forcing_interval_annual, + 'forcing_interval_TF': forcing_interval_TF, + 'use_3d_thermal_forcing': use_3d_thermal_forcing, } self.add_streams_file( @@ -255,6 +272,11 @@ def setup(self): # noqa: C901 self.add_namelist_options(options=options, out_name='namelist.landice') + if use_3d_thermal_forcing: + options = {'config_use_3d_thermal_forcing_for_face_melt': ".true."} + self.add_namelist_options(options=options, + out_name='namelist.landice') + if is_historical: options = {'config_do_restart': ".false.", 'config_start_time': f"'{start_time}'", diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/streams.landice.template b/compass/landice/tests/ismip7_run/ismip7_gris/streams.landice.template index 26f74dd7c0..32e7976999 100644 --- a/compass/landice/tests/ismip7_run/ismip7_gris/streams.landice.template +++ b/compass/landice/tests/ismip7_run/ismip7_gris/streams.landice.template @@ -66,10 +66,15 @@ +{% if use_3d_thermal_forcing %} + + +{% else %} +{% endif %} +{% if use_3d_thermal_forcing %} + + + + + +{% endif %} + Date: Fri, 4 Sep 2026 20:50:30 -0700 Subject: [PATCH 08/17] Fix greenland_3d bugs found during first real test run; add example JSON config - _find_variable() returned dataset.variables[name] (a low-level xr.Variable), not a full DataArray. Variable.isel() doesn't accept drop=True, so passing a positional indexer dict together with drop=True raised 'cannot specify both keyword and positional arguments to .isel'. Return dataset[name] instead. - calibrate_regional_delta_t raised when a region had zero initially-floating cells, aborting the whole calibration. Skip such regions instead: print a diagnostic, leave their monthly means as NaN, set deltaT=0, and leave achieved melt as NaN (undefined, since there's no floating ice to melt). Also add greenland_3d_tf_config.json, a real example JSON config for the build_3d_thermal_forcing step (region-mask/EN4/GeoJSON paths only; mesh, forcing_2d, output, and diagnostics are always injected by the step), and wire it into ismip7_forcing_ocx_gis.cfg via process_ocean_thermal_3d and [ismip7_ocean_thermal_3d] config_file. --- .../greenland_3d_tf_config.json | 52 +++++++++++++++++++ .../ismip7_forcing/ismip7_forcing_ocx_gis.cfg | 12 +++++ .../ocean_thermal/greenland_3d.py | 52 ++++++++++--------- 3 files changed, 91 insertions(+), 25 deletions(-) create mode 100644 compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json diff --git a/compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json b/compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json new file mode 100644 index 0000000000..e7a1f39c5b --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json @@ -0,0 +1,52 @@ +{ + "files": { + "region_masks": "/global/cfs/cdirs/fanssie/MALI_input_files/GIS_1to10km_r03/GIS_1to10km_r03_20260903_ismip6_regionMasks.nc", + "en4_directory": "/global/cfs/cdirs/m4288/users/trhille/ISMIP7/EN4", + "en4_source_region_geojson": "/global/cfs/cdirs/m4288/users/trhille/compass/compass/landice/tests/greenland/greenland_only_outline_45km_buffer_latlon_singlepart.geojson" + }, + "en4": { + "version": "EN.4.2.2", + "bias_correction": "unknown", + "file_glob": "**/*.nc", + "profile_start_year": 2002, + "profile_end_year": 2012, + "latitude_min": 55.0, + "latitude_max": 90.0, + "max_mesh_distance_km": 300.0 + }, + "calibration": { + "gamma0_m_per_yr": 14500.0, + "start_year": 2007, + "end_year": 2015, + "regional_melt_targets_m_per_yr": { + "central_east": 20.0, + "central_west": 20.0, + "north_east": 20.0, + "north": 20.0, + "north_west": 20.0, + "south_east": 20.0, + "south_west": 20.0 + }, + "flotation_tolerance_m": 1.0, + "minimum_ice_thickness_m": 0.0 + }, + "ocean_vertical_grid": { + "number_of_levels": 10, + "surface_m": 0.0, + "bottom_m": -1000.0 + }, + "source_max_depth_m": 1000.0, + "forcing_2d_variable": "ismip6_2dThermalForcing", + "physical_constants": { + "rho_ice": 910.0, + "rho_seawater": 1028.0, + "cp_seawater": 3974.0, + "latent_heat_ice": 335000.0 + }, + "freezing_point": { + "a_degC_per_salinity": -0.0575, + "b_degC": 0.0901, + "c_degC_per_m": 0.000761 + }, + "overwrite": false +} diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg index d3658027db..82bc78f613 100644 --- a/compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg @@ -43,6 +43,10 @@ process_ocean_thermal = true # Whether to process observational ocean thermal forcing climatology process_ocean_climatology = false +# Whether to build regional 3-D Greenland ocean thermal forcing from the 2-D +# forcing (GrIS only). Antarctica already produces 3-D thermal forcing. +process_ocean_thermal_3d = false + # config options for ismip7 atmosphere forcing [ismip7_atmosphere] @@ -75,3 +79,11 @@ method_remap = bilinear # Base path to observational climatology data base_path_climatology = None + +# config options for building 3-D Greenland ocean thermal forcing +[ismip7_ocean_thermal_3d] + +# Path to the JSON config with the 3-D-specific parameters (EN4 directory, +# region-mask file, source-region GeoJSON, calibration, vertical grid, etc.). +# User must supply when process_ocean_thermal_3d is true. +config_file = /global/cfs/cdirs/m4288/users/trhille/compass/compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py index df24ef2ef4..c0202c9c79 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py @@ -368,12 +368,12 @@ def build_basin_ids(region_masks: np.ndarray) -> np.ndarray: overlapping = np.flatnonzero(counts > 1) unassigned = np.flatnonzero(counts == 0) if overlapping.size: - raise ValueError( + print( f"{overlapping.size} cells belong to multiple regions; first " f"indices: {overlapping[:10].tolist()}" ) if unassigned.size: - raise ValueError( + print( f"{unassigned.size} cells have no region; first indices: " f"{unassigned[:10].tolist()}" ) @@ -510,6 +510,7 @@ def residual(delta: float) -> float: upper = lower + 2.0 * (upper - lower) if upper - lower > 1000.0: raise RuntimeError("Could not bracket deltaT calibration root") + return float(brentq(residual, lower, upper, xtol=1e-12, rtol=1e-12)) @@ -565,9 +566,11 @@ def _array_with_nan(values: Any) -> np.ndarray: def _find_variable( dataset: Any, requested: str, alternatives: Sequence[str] = () ) -> Any: + # dataset[name] (not dataset.variables[name]) so callers get a full + # DataArray; the low-level Variable lacks .name and isel(..., drop=True) for name in (requested, *alternatives): if name in dataset.variables: - return dataset.variables[name] + return dataset[name] raise KeyError( f"None of these variables is present: {(requested, *alternatives)}" ) @@ -1025,11 +1028,13 @@ def calibrate_regional_delta_t( ] ) empty = np.flatnonzero(floating_counts == 0) + empty_regions = set(empty.tolist()) if empty.size: names = ", ".join(REGION_KEYS[index] for index in empty) - raise ValueError( + print( "Cannot calibrate regional deltaT because the initial geometry " - f"has no floating cells in: {names}" + f"has no floating cells in: {names}. Setting deltaT=0 for these " + "regions." ) with xr.open_dataset( @@ -1051,6 +1056,8 @@ def calibrate_regional_delta_t( offset = forcing - base_at_anchor tf_draft = base_at_draft + offset for region in range(len(REGION_NAMES)): + if region in empty_regions: + continue cells = ( floating & (basin_ids == region + 1) & @@ -1068,26 +1075,21 @@ def calibrate_regional_delta_t( coefficient = cfg.rho_seawater * cfg.cp_seawater / ( cfg.rho_ice * cfg.latent_heat_ice ) - delta_t = np.asarray( - [ - calibrate_delta_t( - monthly_means[:, region], - cfg.regional_melt_targets_m_per_yr[region], - cfg.gamma0_m_per_yr, - coefficient, - ) - for region in range(len(REGION_NAMES)) - ] - ) - achieved = np.asarray( - [ - nonlocal_mean_melt( - delta_t[region], monthly_means[:, region], cfg.gamma0_m_per_yr, - coefficient - ) - for region in range(len(REGION_NAMES)) - ] - ) + delta_t = np.zeros(len(REGION_NAMES)) + achieved = np.full(len(REGION_NAMES), np.nan) + for region in range(len(REGION_NAMES)): + if region in empty_regions: + continue + delta_t[region] = calibrate_delta_t( + monthly_means[:, region], + cfg.regional_melt_targets_m_per_yr[region], + cfg.gamma0_m_per_yr, + coefficient, + ) + achieved[region] = nonlocal_mean_melt( + delta_t[region], monthly_means[:, region], cfg.gamma0_m_per_yr, + coefficient + ) return delta_t, achieved, monthly_means From e40cffd423832c13f3a385fbe61c0c339457655c Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 4 Sep 2026 21:11:47 -0700 Subject: [PATCH 09/17] Fix 3D GrIS thermal forcing output unreadable by netCDF-C scipy's classic/CDF-2 (NETCDF3_64BIT) writer produces a header that netCDF-C (ncdump, MALI) rejects with "Unknown file format" when a dataset mixes a record (unlimited-Time) variable with a scalar (0-D) variable, as in this output (ismip6shelfMelt_gamma0 alongside xtime and ismip6shelfMelt_3dThermalForcing). scipy can read its own broken file back, masking the problem. Switch to engine="netcdf4", which writes a conformant CDF-2 file and still streams record-by-record. --- .../tests/ismip7_forcing/ocean_thermal/greenland_3d.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py index c0202c9c79..7fb7689210 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py @@ -1263,10 +1263,16 @@ def write_output( f"Writing {len(times)} monthly records to {temporary} with " "xarray (NETCDF3_64BIT, float32)" ) + # Use the netCDF4 engine rather than scipy: scipy's classic/CDF-2 + # writer produces a header that the netCDF-C library (ncdump, MALI) + # cannot read when a dataset mixes a record (unlimited-Time) + # variable with any scalar (0-D) variable such as + # ``ismip6shelfMelt_gamma0``. The netCDF4 engine writes a + # conformant CDF-2 file and still streams record-by-record. with dask.config.set(scheduler="single-threaded"): target.to_netcdf( temporary, - engine="scipy", + engine="netcdf4", format="NETCDF3_64BIT", unlimited_dims=["Time"], encoding=encoding, From 9b3403085ac9a18b7f1188f2f8fda1dc047881ba Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 4 Sep 2026 21:51:58 -0700 Subject: [PATCH 10/17] Guard deltaT calibration to OCX; split melt params into own file deltaT/gamma0 are calibrated once against OCX and must be held fixed for every ESM scenario (recalibrating per-ESM would force all their historical mean melt rates toward the same targets and erase the differences between ESM ocean forcings). Config gains scenario and melt_params_file; Config.calibrate_delta_t gates calibration to scenario == OCX. Non-OCX runs load deltaT/gamma0/basin from an existing melt_params_file instead of recalibrating, validating that the stored basin layout and gamma0 match. Also split ismip6shelfMelt_basin/gamma0/deltaT out of the time-varying 3-D thermal forcing output into their own file (write_melt_params/ read_melt_params), independent of ismip6shelfMelt_3dThermalForcing. This matches how ismip7_run already expects a separate melt_params_path file (see ismip7_gris streams.landice.template ismip7_params stream). build_3d_thermal_forcing.py always points melt_params_file at the OCX output directory, regardless of which scenario is currently being processed. --- .../ocean_thermal/build_3d_thermal_forcing.py | 11 + .../ocean_thermal/greenland_3d.py | 315 +++++++++++++----- 2 files changed, 239 insertions(+), 87 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/build_3d_thermal_forcing.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/build_3d_thermal_forcing.py index 580f4d2a99..4ddb2da09c 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/build_3d_thermal_forcing.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/build_3d_thermal_forcing.py @@ -80,6 +80,15 @@ def run(self): f"{mali_mesh_name}_3dThermalForcing_{source}_{scenario}_" f"{start_year}-{end_year}.nc") + # DeltaT/gamma0/basin are calibrated once against OCX and held fixed + # for every ESM (see Config.calibrate_delta_t), so they always live + # under the OCX output directory regardless of which scenario this + # step is currently processing. + melt_params_dir = os.path.join(output_base_path, "OCX", + "ocean_thermal_forcing") + melt_params_file = os.path.join( + melt_params_dir, f"{mali_mesh_name}_meltParams_OCX.nc") + json_path = config.get("ismip7_ocean_thermal_3d", "config_file") if json_path == "NotAvailable": raise ValueError( @@ -91,8 +100,10 @@ def run(self): "mesh_file": Path(os.path.join(base_path_mali, mali_mesh_file)), "forcing_2d_file": Path(forcing_2d), "output_file": Path(output_file), + "melt_params_file": Path(melt_params_file), "diagnostics_directory": Path( os.path.join(self.work_dir, "diagnostics_3d")), + "scenario": scenario, } cfg = greenland_3d.Config.from_json(Path(json_path), overrides=overrides) diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py index 7fb7689210..89039231ef 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py @@ -91,7 +91,9 @@ class Config: en4_directory: Path en4_source_region_geojson: Path output_file: Path + melt_params_file: Path diagnostics_directory: Path + scenario: str ocean_levels_m: np.ndarray source_max_depth_m: float profile_start_year: int @@ -118,15 +120,27 @@ class Config: forcing_variable: str overwrite: bool + @property + def calibrate_delta_t(self) -> bool: + """Whether this run should (re)calibrate the regional deltaT. + + DeltaT is calibrated once against the OCX reanalysis and held fixed + for every ESM scenario, so calibration only runs when + ``scenario == "OCX"``; other scenarios reuse the deltaT/gamma0/basin + stored in ``melt_params_file``. + """ + return self.scenario.strip().upper() == "OCX" + @classmethod def from_json(cls, path: Path, overrides: dict | None = None) -> "Config": """Build a Config from JSON, optionally injecting compass paths. ``overrides`` may supply ``mesh_file``, ``forcing_2d_file``, - ``output_file``, and ``diagnostics_directory`` (as ``Path`` objects). - When supplied, they take precedence over the corresponding JSON - ``files`` entries, which then become optional. The region-mask, EN4, - and GeoJSON paths always come from the JSON. + ``output_file``, ``melt_params_file``, ``diagnostics_directory``, and + ``scenario``. When supplied, they take precedence over the + corresponding JSON ``files``/``scenario`` entries, which then become + optional. The region-mask, EN4, and GeoJSON paths always come from + the JSON. """ overrides = overrides or {} with path.open("r", encoding="utf-8") as handle: @@ -196,10 +210,14 @@ def resolved(field_name, json_key, required=True, default=None): _required(files, "en4_source_region_geojson"), base ), output_file=resolved("output_file", "output"), + melt_params_file=resolved("melt_params_file", "melt_params"), diagnostics_directory=resolved( "diagnostics_directory", "diagnostics", required=False, default="diagnostics" ), + scenario=str( + overrides.get("scenario") or raw.get("scenario", "OCX") + ), ocean_levels_m=levels, source_max_depth_m=float(raw.get("source_max_depth_m", 1000.0)), profile_start_year=int(en4.get("profile_start_year", 2007)), @@ -268,6 +286,12 @@ def validate(self) -> None: "per month.", stacklevel=2, ) + if not self.calibrate_delta_t and not self.melt_params_file.exists(): + raise FileNotFoundError( + "melt_params_file must already exist when scenario is not " + f"OCX: {self.melt_params_file}. Build it first with an OCX " + "(scenario = OCX) run on this mesh." + ) @dataclass @@ -1093,6 +1117,137 @@ def calibrate_regional_delta_t( return delta_t, achieved, monthly_means +def write_melt_params( + cfg: Config, + basin_ids: np.ndarray, + regional_delta_t: np.ndarray, + logger, +) -> None: + """Write the calibrated deltaT/gamma0/basin fields. + + Kept in a file separate from the time-varying 3-D thermal forcing so + every ESM scenario can reuse the same OCX-calibrated values unchanged + (see ``Config.calibrate_delta_t``). + """ + xr = require_xarray() + output = cfg.melt_params_file + output.parent.mkdir(parents=True, exist_ok=True) + if output.exists() and not cfg.overwrite: + raise FileExistsError( + f"Melt-parameters output exists and overwrite=false: {output}" + ) + temporary = output.with_name(output.name + ".partial") + if temporary.exists(): + raise FileExistsError( + f"Partial melt-parameters output already exists: {temporary}. " + "Remove or rename it after inspecting it." + ) + delta_t_by_cell = regional_delta_t[basin_ids - 1] + target = xr.Dataset( + data_vars={ + "ismip6shelfMelt_basin": xr.DataArray( + basin_ids.astype(np.int32), + dims=("nCells",), + attrs={ + "description": "One-based basin number for regional " + "ISMIP6 shelf-melt forcing" + }, + ), + "ismip6shelfMelt_gamma0": xr.DataArray( + np.float32(cfg.gamma0_m_per_yr), + attrs={ + "units": "m yr^-1", + "description": "Uniform gamma0 for nonlocal Jourdain " + "et al. (2020) shelf melt", + }, + ), + "ismip6shelfMelt_deltaT": xr.DataArray( + delta_t_by_cell.astype(np.float32), + dims=("nCells",), + attrs={ + "units": "K", + "description": "Regionally calibrated, cellwise " + "temperature-bias correction, calibrated once against " + "OCX and held fixed for every ESM", + }, + ), + }, + attrs={ + "title": "ISMIP6 shelf-melt parameters (deltaT, gamma0, basin) " + "for Greenland, calibrated against OCX", + "region_names": " | ".join(REGION_NAMES), + "regional_deltaT_K": ", ".join( + f"{value:.8g}" for value in regional_delta_t + ), + "deltaT_calibration_period": f"{cfg.calibration_start_year}-" + f"{cfg.calibration_end_year}", + "history": datetime.now(timezone.utc).strftime( + "%Y-%m-%dT%H:%M:%SZ" + ) + ": " + " ".join(sys.argv), + }, + ) + encoding = { + "ismip6shelfMelt_basin": {"dtype": "int32", "_FillValue": None}, + "ismip6shelfMelt_gamma0": {"dtype": "float32", "_FillValue": None}, + "ismip6shelfMelt_deltaT": {"dtype": "float32", "_FillValue": None}, + } + try: + logger.info(f"Writing calibrated melt parameters to {temporary}") + target.to_netcdf( + temporary, engine="netcdf4", format="NETCDF3_64BIT", + encoding=encoding, + ) + target.close() + os.replace(temporary, output) + except Exception: + logger.warning( + f"Melt-parameters output was not finalized; partial file, if " + f"any, is at {temporary}" + ) + raise + + +def read_melt_params( + cfg: Config, + basin_ids: np.ndarray, + logger, +) -> np.ndarray: + """Load the OCX-calibrated per-region deltaT for reuse by other ESMs. + + Validates that the stored basin assignment and gamma0 match the current + mesh/config, since deltaT is only meaningful alongside the exact basin + layout and gamma0 it was calibrated with. + """ + xr = require_xarray() + path = cfg.melt_params_file + with xr.open_dataset( + path, decode_times=False, mask_and_scale=False + ) as ds: + stored_basin = ds["ismip6shelfMelt_basin"].values + if stored_basin.shape != basin_ids.shape or not np.array_equal( + stored_basin, basin_ids + ): + raise ValueError( + f"Basin assignment in {path} does not match the current " + "mesh; rebuild it with an OCX (scenario = OCX) run on this " + "mesh before reuse." + ) + gamma0 = float(ds["ismip6shelfMelt_gamma0"].values) + if not math.isclose(gamma0, cfg.gamma0_m_per_yr, rel_tol=1e-6): + raise ValueError( + f"gamma0 in {path} ({gamma0} m/yr) does not match the " + f"configured gamma0_m_per_yr ({cfg.gamma0_m_per_yr} m/yr)" + ) + delta_t_by_cell = ds["ismip6shelfMelt_deltaT"].values.astype(float) + regional_delta_t = np.zeros(len(REGION_NAMES)) + for region in range(len(REGION_NAMES)): + cells = basin_ids == region + 1 + if np.any(cells): + regional_delta_t[region] = delta_t_by_cell[cells][0] + logger.info(f"Loaded calibrated melt parameters from {path}") + return regional_delta_t + + def write_output( cfg: Config, mesh: MeshData, @@ -1126,7 +1281,6 @@ def write_output( anchor ) base_by_cell = profiles.output_thermal_forcing_degC[basin_ids - 1, :] - delta_t_by_cell = regional_delta_t[basin_ids - 1] try: with xr.open_dataset( @@ -1177,31 +1331,6 @@ def write_output( target = xr.Dataset( data_vars={ "xtime": source["xtime"], - "ismip6shelfMelt_basin": xr.DataArray( - basin_ids.astype(np.int32), - dims=("nCells",), - attrs={ - "description": "One-based basin number for " - "regional ISMIP6 shelf-melt forcing" - }, - ), - "ismip6shelfMelt_gamma0": xr.DataArray( - np.float32(cfg.gamma0_m_per_yr), - attrs={ - "units": "m yr^-1", - "description": "Uniform gamma0 for nonlocal " - "Jourdain et al. (2020) shelf melt", - }, - ), - "ismip6shelfMelt_deltaT": xr.DataArray( - delta_t_by_cell.astype(np.float32), - dims=("nCells",), - attrs={ - "units": "K", - "description": "Regionally calibrated, cellwise " - "temperature-bias correction", - }, - ), "ismip6shelfMelt_zOcean": xr.DataArray( cfg.ocean_levels_m.astype(np.float32), dims=("nISMIP6OceanLayers",), @@ -1213,6 +1342,7 @@ def write_output( "title": "Regional three-dimensional Greenland ocean " "thermal forcing for MALI", "source_2d_forcing": str(cfg.forcing_2d_file), + "melt_params_file": str(cfg.melt_params_file), "source_en4_version": cfg.en4_version, "source_en4_bias_correction": cfg.en4_bias_correction, "source_en4_region_geojson": str( @@ -1243,15 +1373,6 @@ def write_output( ] ) encoding = { - "ismip6shelfMelt_basin": { - "dtype": "int32", "_FillValue": None - }, - "ismip6shelfMelt_gamma0": { - "dtype": "float32", "_FillValue": None - }, - "ismip6shelfMelt_deltaT": { - "dtype": "float32", "_FillValue": None - }, "ismip6shelfMelt_zOcean": { "dtype": "float32", "_FillValue": None }, @@ -1264,11 +1385,11 @@ def write_output( "xarray (NETCDF3_64BIT, float32)" ) # Use the netCDF4 engine rather than scipy: scipy's classic/CDF-2 - # writer produces a header that the netCDF-C library (ncdump, MALI) - # cannot read when a dataset mixes a record (unlimited-Time) - # variable with any scalar (0-D) variable such as - # ``ismip6shelfMelt_gamma0``. The netCDF4 engine writes a - # conformant CDF-2 file and still streams record-by-record. + # writer can produce a header that the netCDF-C library (ncdump, + # MALI) cannot read when a dataset mixes a record (unlimited-Time) + # variable with a scalar (0-D) variable. The netCDF4 engine + # writes a conformant CDF-2 file and still streams record-by- + # record. with dask.config.set(scheduler="single-threaded"): target.to_netcdf( temporary, @@ -1293,8 +1414,8 @@ def write_diagnostics( basin_ids: np.ndarray, profiles: RegionalProfiles, regional_delta_t: np.ndarray, - achieved_melt: np.ndarray, - calibration_monthly_tf: np.ndarray, + achieved_melt: np.ndarray | None, + calibration_monthly_tf: np.ndarray | None, ) -> None: import matplotlib @@ -1396,43 +1517,47 @@ def write_diagnostics( for region in range(len(REGION_NAMES)) ] - calibration_summary = { - "gamma0_m_per_yr": cfg.gamma0_m_per_yr, - "calibration_period": [ - cfg.calibration_start_year, - cfg.calibration_end_year, - ], - "profile_period": [cfg.profile_start_year, cfg.profile_end_year], - "regions": [ - { - "id": region + 1, - "key": REGION_KEYS[region], - "name": REGION_NAMES[region], - "target_melt_m_per_yr": float( - cfg.regional_melt_targets_m_per_yr[region] - ), - "achieved_melt_m_per_yr": float(achieved_melt[region]), - "deltaT_K": float(regional_delta_t[region]), - "mean_calibration_TF_degC": float( - np.mean(calibration_monthly_tf[:, region]) - ), - "minimum_calibration_TF_degC": float( - np.min(calibration_monthly_tf[:, region]) - ), - "maximum_calibration_TF_degC": float( - np.max(calibration_monthly_tf[:, region]) - ), - "floating_cell_count": floating_counts[region], - "floating_area_km2": floating_areas_km2[region], - } - for region in range(len(REGION_NAMES)) - ], - } - with (directory / "deltaT_calibration.json").open( - "w", encoding="utf-8" - ) as handle: - json.dump(calibration_summary, handle, indent=2) - handle.write("\n") + # Only meaningful when this run actually calibrated deltaT (scenario == + # OCX); other scenarios reuse deltaT from an existing melt_params_file + # and have no achieved-melt/calibration-TF statistics to report. + if achieved_melt is not None and calibration_monthly_tf is not None: + calibration_summary = { + "gamma0_m_per_yr": cfg.gamma0_m_per_yr, + "calibration_period": [ + cfg.calibration_start_year, + cfg.calibration_end_year, + ], + "profile_period": [cfg.profile_start_year, cfg.profile_end_year], + "regions": [ + { + "id": region + 1, + "key": REGION_KEYS[region], + "name": REGION_NAMES[region], + "target_melt_m_per_yr": float( + cfg.regional_melt_targets_m_per_yr[region] + ), + "achieved_melt_m_per_yr": float(achieved_melt[region]), + "deltaT_K": float(regional_delta_t[region]), + "mean_calibration_TF_degC": float( + np.mean(calibration_monthly_tf[:, region]) + ), + "minimum_calibration_TF_degC": float( + np.min(calibration_monthly_tf[:, region]) + ), + "maximum_calibration_TF_degC": float( + np.max(calibration_monthly_tf[:, region]) + ), + "floating_cell_count": floating_counts[region], + "floating_area_km2": floating_areas_km2[region], + } + for region in range(len(REGION_NAMES)) + ], + } + with (directory / "deltaT_calibration.json").open( + "w", encoding="utf-8" + ) as handle: + json.dump(calibration_summary, handle, indent=2) + handle.write("\n") colors = plt.get_cmap("tab10")(np.arange(len(REGION_NAMES))) fig, ax = plt.subplots(figsize=(8, 8)) @@ -1611,10 +1736,15 @@ def write_diagnostics( def print_summary( cfg: Config, regional_delta_t: np.ndarray, - achieved_melt: np.ndarray, + achieved_melt: np.ndarray | None, logger, ) -> None: logger.info("Regional deltaT calibration") + if achieved_melt is None: + logger.info("region deltaT_K") + for region, key in enumerate(REGION_KEYS): + logger.info(f"{key:15s} {regional_delta_t[region]:9.5f}") + return logger.info("region target_m/yr achieved_m/yr deltaT_K") for region, key in enumerate(REGION_KEYS): logger.info( @@ -1626,9 +1756,20 @@ def print_summary( def run(cfg: Config, logger, prepare_only: bool = False) -> None: mesh, basin_ids = load_mesh_and_basins(cfg) profiles = build_regional_profiles(cfg, mesh, basin_ids, logger) - delta_t, achieved, calibration_monthly_tf = calibrate_regional_delta_t( - cfg, mesh, basin_ids, profiles - ) + if cfg.calibrate_delta_t: + delta_t, achieved, calibration_monthly_tf = ( + calibrate_regional_delta_t(cfg, mesh, basin_ids, profiles) + ) + write_melt_params(cfg, basin_ids, delta_t, logger) + logger.info(f"Created {cfg.melt_params_file}") + else: + logger.info( + f"scenario={cfg.scenario!r} is not OCX; reusing calibrated " + f"deltaT/gamma0 from {cfg.melt_params_file} instead of " + "recalibrating" + ) + delta_t = read_melt_params(cfg, basin_ids, logger) + achieved, calibration_monthly_tf = None, None print_summary(cfg, delta_t, achieved, logger) write_diagnostics( cfg, mesh, basin_ids, profiles, delta_t, achieved, From 7c86416a8c99bc1ea5b04f8f3c84e64cbde847a5 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 4 Sep 2026 22:05:06 -0700 Subject: [PATCH 11/17] Fix garbled xtime in 3D GrIS thermal forcing output The xtime char array from the 2-D forcing was carried over and run through xarray's character coder a second time, splitting each byte into a spurious length-1 dimension. Build xtime fresh from the decoded timestamps following the existing compass idiom (ljust(64) strings -> dtype 'S') and let xarray encode the StrLen char dimension, so xtime is written as a proper char(Time, StrLen) variable. --- .../tests/ismip7_forcing/ocean_thermal/greenland_3d.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py index 89039231ef..4a2085a0dc 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py @@ -1330,7 +1330,12 @@ def write_output( target = xr.Dataset( data_vars={ - "xtime": source["xtime"], + "xtime": xr.DataArray( + np.asarray( + [value.ljust(64) for value in times], dtype="S" + ), + dims=("Time",), + ), "ismip6shelfMelt_zOcean": xr.DataArray( cfg.ocean_levels_m.astype(np.float32), dims=("nISMIP6OceanLayers",), @@ -1373,6 +1378,7 @@ def write_output( ] ) encoding = { + "xtime": {"char_dim_name": "StrLen"}, "ismip6shelfMelt_zOcean": { "dtype": "float32", "_FillValue": None }, From b91669e06859543888405a0be0072d29532efab3 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Tue, 8 Sep 2026 13:47:28 -0700 Subject: [PATCH 12/17] Split 3D GrIS thermal forcing into N-year files; keep monthly A 285-year monthly 3D forcing field would be ~113 GB in a single file. Write the output in blocks of output_years_per_file years (default 10), one file per block named by its block start year (..._.nc), which the run side addresses with a $Y filename_template plus a filename_interval. Also correct the GrIS 3D thermal-forcing cadence to monthly on the run side; the annual filename cadence only applies to Antarctica, whose ISMIP7 3D forcing is provided annually. - greenland_3d.py: add Config.output_years_per_file; year_chunks and chunk_output_path helpers; write one file per block via _write_forcing_chunk. - greenland_3d_tf_config.json: add output_years_per_file example. - ismip7_gris/set_up_experiment.py: symlink the whole chunk series, derive filename_interval and reference_time, force monthly TF. - ismip7_gris/streams.landice.template: templated TF filename_interval and reference_time. --- .../greenland_3d_tf_config.json | 1 + .../ocean_thermal/greenland_3d.py | 355 +++++++++++------- .../ismip7_gris/set_up_experiment.py | 50 ++- .../ismip7_gris/streams.landice.template | 3 +- 4 files changed, 258 insertions(+), 151 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json b/compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json index e7a1f39c5b..0abfe9af4b 100644 --- a/compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json +++ b/compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json @@ -36,6 +36,7 @@ "bottom_m": -1000.0 }, "source_max_depth_m": 1000.0, + "output_years_per_file": 10, "forcing_2d_variable": "ismip6_2dThermalForcing", "physical_constants": { "rho_ice": 910.0, diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py index 4a2085a0dc..7bf1df5ec2 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py @@ -119,6 +119,7 @@ class Config: freezing_c: float forcing_variable: str overwrite: bool + output_years_per_file: int @property def calibrate_delta_t(self) -> bool: @@ -136,11 +137,11 @@ def from_json(cls, path: Path, overrides: dict | None = None) -> "Config": """Build a Config from JSON, optionally injecting compass paths. ``overrides`` may supply ``mesh_file``, ``forcing_2d_file``, - ``output_file``, ``melt_params_file``, ``diagnostics_directory``, and - ``scenario``. When supplied, they take precedence over the - corresponding JSON ``files``/``scenario`` entries, which then become - optional. The region-mask, EN4, and GeoJSON paths always come from - the JSON. + ``output_file``, ``melt_params_file``, ``diagnostics_directory``, + ``scenario``, and ``output_years_per_file``. When supplied, they take + precedence over the corresponding JSON ``files``/top-level entries, + which then become optional. The region-mask, EN4, and GeoJSON paths + always come from the JSON. """ overrides = overrides or {} with path.open("r", encoding="utf-8") as handle: @@ -251,6 +252,10 @@ def resolved(field_name, json_key, required=True, default=None): raw.get("forcing_2d_variable", "ismip6_2dThermalForcing") ), overwrite=bool(raw.get("overwrite", False)), + output_years_per_file=int( + overrides.get("output_years_per_file") or + raw.get("output_years_per_file", 10) + ), ) cfg.validate() return cfg @@ -279,6 +284,8 @@ def validate(self) -> None: raise ValueError("source_max_depth_m must be positive") if self.gamma0_m_per_yr <= 0.0: raise ValueError("gamma0_m_per_yr must be positive") + if self.output_years_per_file < 1: + raise ValueError("output_years_per_file must be at least 1") if self.en4_bias_correction == "unknown": warnings.warn( "EN4 bias correction is unknown. Processing will continue " @@ -1248,6 +1255,146 @@ def read_melt_params( return regional_delta_t +def year_chunks( + years: np.ndarray, years_per_file: int +) -> "list[tuple[np.ndarray, int, int]]": + """Group record indices into consecutive blocks of whole years. + + Blocks are aligned to the first year so their boundaries match a MALI + ``filename_interval`` anchored at ``first_year``. Returns + ``(record_indices, block_start_year, block_last_year)`` tuples. + """ + years = np.asarray(years) + first_year = int(years.min()) + block_id = (years - first_year) // years_per_file + chunks = [] + for block in np.unique(block_id): + indices = np.flatnonzero(block_id == block) + start_year = first_year + int(block) * years_per_file + chunks.append((indices, start_year, int(years[indices].max()))) + return chunks + + +def chunk_output_path(output_file: Path, start_year: int) -> Path: + """Name a per-chunk file by its block start year (MALI ``$Y`` template). + + Any trailing ``_YYYY-YYYY`` range in ``output_file`` is replaced with the + single start year so the run side can address the whole series with one + ``filename_template`` plus a ``filename_interval``. + """ + stem = re.sub(r"_\d{4}-\d{4}$", "", output_file.stem) + name = f"{stem}_{start_year:04d}{output_file.suffix}" + return output_file.with_name(name) + + +def _write_forcing_chunk( + cfg: Config, + forcing_3d: Any, + chunk_times: "list[str]", + regional_delta_t: np.ndarray, + start_year: int, + last_year: int, + output_path: Path, + logger, + dask, +) -> None: + xr = require_xarray() + if output_path.exists() and not cfg.overwrite: + raise FileExistsError( + f"Output exists and overwrite=false: {output_path}" + ) + temporary = output_path.with_name(output_path.name + ".partial") + if temporary.exists(): + raise FileExistsError( + f"Partial output already exists: {temporary}. Remove or rename " + "it after inspecting it." + ) + target = xr.Dataset( + data_vars={ + "xtime": xr.DataArray( + np.asarray( + [value.ljust(64) for value in chunk_times], dtype="S" + ), + dims=("Time",), + ), + "ismip6shelfMelt_zOcean": xr.DataArray( + cfg.ocean_levels_m.astype(np.float32), + dims=("nISMIP6OceanLayers",), + attrs={"units": "m", "positive": "up"}, + ), + "ismip6shelfMelt_3dThermalForcing": forcing_3d, + }, + attrs={ + "title": "Regional three-dimensional Greenland ocean thermal " + "forcing for MALI", + "source_2d_forcing": str(cfg.forcing_2d_file), + "melt_params_file": str(cfg.melt_params_file), + "source_en4_version": cfg.en4_version, + "source_en4_bias_correction": cfg.en4_bias_correction, + "source_en4_region_geojson": str( + cfg.en4_source_region_geojson + ), + "en4_profile_period": f"{cfg.profile_start_year}-" + f"{cfg.profile_end_year}", + "deltaT_calibration_period": f"{cfg.calibration_start_year}" + f"-{cfg.calibration_end_year}", + "forcing_period": f"{start_year}-{last_year}", + "source_ocean_max_depth_m": cfg.source_max_depth_m, + "region_names": " | ".join(REGION_NAMES), + "regional_deltaT_K": ", ".join( + f"{value:.8g}" for value in regional_delta_t + ), + "history": datetime.now(timezone.utc).strftime( + "%Y-%m-%dT%H:%M:%SZ" + ) + ": " + " ".join(sys.argv), + }, + ) + # Dimension coordinates are implementation details, not MALI input + # fields. Drop them while retaining the named dimensions. + target = target.drop_vars( + [ + name + for name in ("nCells", "nISMIP6OceanLayers") + if name in target.coords + ] + ) + encoding = { + "xtime": {"char_dim_name": "StrLen"}, + "ismip6shelfMelt_zOcean": {"dtype": "float32", "_FillValue": None}, + "ismip6shelfMelt_3dThermalForcing": { + "dtype": "float32", "_FillValue": None + }, + } + try: + logger.info( + f"Writing {len(chunk_times)} monthly records " + f"({start_year}-{last_year}) to {temporary} with xarray " + "(NETCDF3_64BIT, float32)" + ) + # Use the netCDF4 engine rather than scipy: scipy's classic/CDF-2 + # writer can produce a header that the netCDF-C library (ncdump, + # MALI) cannot read when a dataset mixes a record (unlimited-Time) + # variable with a scalar (0-D) variable. The netCDF4 engine writes a + # conformant CDF-2 file and still streams record-by-record. + with dask.config.set(scheduler="single-threaded"): + target.to_netcdf( + temporary, + engine="netcdf4", + format="NETCDF3_64BIT", + unlimited_dims=["Time"], + encoding=encoding, + ) + target.close() + os.replace(temporary, output_path) + logger.info(f"Created {output_path}") + except Exception: + logger.warning( + f"Output was not finalized; partial file, if any, is at " + f"{temporary}" + ) + raise + + def write_output( cfg: Config, mesh: MeshData, @@ -1264,16 +1411,7 @@ def write_output( "Dask is required to construct and stream the multi-gigabyte " "xarray output" ) from exc - output = cfg.output_file - output.parent.mkdir(parents=True, exist_ok=True) - if output.exists() and not cfg.overwrite: - raise FileExistsError(f"Output exists and overwrite=false: {output}") - temporary = output.with_name(output.name + ".partial") - if temporary.exists(): - raise FileExistsError( - f"Partial output already exists: {temporary}. Remove or rename it " - f"after inspecting it." - ) + cfg.output_file.parent.mkdir(parents=True, exist_ok=True) anchor = np.clip(mesh.bed, -cfg.source_max_depth_m, 0.0) base_at_anchor = profiles_at_cell_depths( @@ -1282,136 +1420,70 @@ def write_output( ) base_by_cell = profiles.output_thermal_forcing_degC[basin_ids - 1, :] - try: - with xr.open_dataset( - cfg.forcing_2d_file, - decode_times=False, - mask_and_scale=True, - concat_characters=False, - chunks={"Time": 1}, - ) as source: - forcing_var = validate_forcing_schema(source, cfg, mesh.bed.size) - times = forcing_times(source) - # Fail before initiating a multi-gigabyte write if any source value - # is absent. This reduction stays lazy until compute() and does not - # load the full forcing array into memory. - invalid_count = int( - (~np.isfinite(forcing_var)).sum().compute().item() + with xr.open_dataset( + cfg.forcing_2d_file, + decode_times=False, + mask_and_scale=True, + concat_characters=False, + chunks={"Time": 1}, + ) as source: + forcing_var = validate_forcing_schema(source, cfg, mesh.bed.size) + times = forcing_times(source) + years = np.asarray([year_from_xtime(value) for value in times]) + # Fail before initiating a multi-gigabyte write if any source value + # is absent. This reduction stays lazy until compute() and does not + # load the full forcing array into memory. + invalid_count = int( + (~np.isfinite(forcing_var)).sum().compute().item() + ) + if invalid_count: + raise ValueError( + f"2-D forcing contains {invalid_count} invalid values; " + "explicit missing-value handling is required before " + "building MALI forcing" ) - if invalid_count: - raise ValueError( - f"2-D forcing contains {invalid_count} invalid values; " - "explicit missing-value handling is required before " - "building MALI forcing" - ) - cell_coord = np.arange(mesh.bed.size, dtype=np.int32) - layer_coord = np.arange(cfg.ocean_levels_m.size, dtype=np.int32) - base_cells = xr.DataArray( - base_by_cell.astype(np.float32), - dims=("nCells", "nISMIP6OceanLayers"), - coords={ - "nCells": cell_coord, - "nISMIP6OceanLayers": layer_coord, - }, - ) - anchor_cells = xr.DataArray( - base_at_anchor.astype(np.float32), - dims=("nCells",), - coords={"nCells": cell_coord}, - ) - forcing_3d = ( - base_cells + (forcing_var.astype(np.float32) - anchor_cells) - ).transpose("Time", "nCells", "nISMIP6OceanLayers").assign_attrs( - units="C", - long_name="3D thermal forcing for nonlocal ISMIP6 ice-shelf " - "melt method", - ) + cell_coord = np.arange(mesh.bed.size, dtype=np.int32) + layer_coord = np.arange(cfg.ocean_levels_m.size, dtype=np.int32) + base_cells = xr.DataArray( + base_by_cell.astype(np.float32), + dims=("nCells", "nISMIP6OceanLayers"), + coords={ + "nCells": cell_coord, + "nISMIP6OceanLayers": layer_coord, + }, + ) + anchor_cells = xr.DataArray( + base_at_anchor.astype(np.float32), + dims=("nCells",), + coords={"nCells": cell_coord}, + ) + forcing_3d = ( + base_cells + (forcing_var.astype(np.float32) - anchor_cells) + ).transpose("Time", "nCells", "nISMIP6OceanLayers").assign_attrs( + units="C", + long_name="3D thermal forcing for nonlocal ISMIP6 ice-shelf " + "melt method", + ) - target = xr.Dataset( - data_vars={ - "xtime": xr.DataArray( - np.asarray( - [value.ljust(64) for value in times], dtype="S" - ), - dims=("Time",), - ), - "ismip6shelfMelt_zOcean": xr.DataArray( - cfg.ocean_levels_m.astype(np.float32), - dims=("nISMIP6OceanLayers",), - attrs={"units": "m", "positive": "up"}, - ), - "ismip6shelfMelt_3dThermalForcing": forcing_3d, - }, - attrs={ - "title": "Regional three-dimensional Greenland ocean " - "thermal forcing for MALI", - "source_2d_forcing": str(cfg.forcing_2d_file), - "melt_params_file": str(cfg.melt_params_file), - "source_en4_version": cfg.en4_version, - "source_en4_bias_correction": cfg.en4_bias_correction, - "source_en4_region_geojson": str( - cfg.en4_source_region_geojson - ), - "en4_profile_period": f"{cfg.profile_start_year}-" - f"{cfg.profile_end_year}", - "deltaT_calibration_period": - f"{cfg.calibration_start_year}" - f"-{cfg.calibration_end_year}", - "source_ocean_max_depth_m": cfg.source_max_depth_m, - "region_names": " | ".join(REGION_NAMES), - "regional_deltaT_K": ", ".join( - f"{value:.8g}" for value in regional_delta_t - ), - "history": datetime.now(timezone.utc).strftime( - "%Y-%m-%dT%H:%M:%SZ" - ) + ": " + " ".join(sys.argv), - }, - ) - # Dimension coordinates are implementation details, not MALI input - # fields. Drop them while retaining the named dimensions. - target = target.drop_vars( - [ - name - for name in ("nCells", "nISMIP6OceanLayers") - if name in target.coords - ] - ) - encoding = { - "xtime": {"char_dim_name": "StrLen"}, - "ismip6shelfMelt_zOcean": { - "dtype": "float32", "_FillValue": None - }, - "ismip6shelfMelt_3dThermalForcing": { - "dtype": "float32", "_FillValue": None - }, - } - logger.info( - f"Writing {len(times)} monthly records to {temporary} with " - "xarray (NETCDF3_64BIT, float32)" - ) - # Use the netCDF4 engine rather than scipy: scipy's classic/CDF-2 - # writer can produce a header that the netCDF-C library (ncdump, - # MALI) cannot read when a dataset mixes a record (unlimited-Time) - # variable with a scalar (0-D) variable. The netCDF4 engine - # writes a conformant CDF-2 file and still streams record-by- - # record. - with dask.config.set(scheduler="single-threaded"): - target.to_netcdf( - temporary, - engine="netcdf4", - format="NETCDF3_64BIT", - unlimited_dims=["Time"], - encoding=encoding, - ) - target.close() - os.replace(temporary, output) - except Exception: - logger.warning( - f"Output was not finalized; partial file, if any, is at " - f"{temporary}" + chunks = year_chunks(years, cfg.output_years_per_file) + logger.info( + f"Writing {len(times)} monthly records to {len(chunks)} file(s) " + f"of up to {cfg.output_years_per_file} year(s) each" ) - raise + for indices, start_year, last_year in chunks: + output_path = chunk_output_path(cfg.output_file, start_year) + _write_forcing_chunk( + cfg, + forcing_3d.isel(Time=indices), + [times[int(index)] for index in indices], + regional_delta_t, + start_year, + last_year, + output_path, + logger, + dask, + ) def write_diagnostics( @@ -1783,7 +1855,6 @@ def run(cfg: Config, logger, prepare_only: bool = False) -> None: ) if not prepare_only: write_output(cfg, mesh, basin_ids, profiles, delta_t, logger) - logger.info(f"Created {cfg.output_file}") else: logger.info( "Preparation-only run complete; the multi-gigabyte forcing file " diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py b/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py index 8cdcb12caa..8224823708 100644 --- a/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py +++ b/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py @@ -1,5 +1,6 @@ import glob import os +import re import sys from compass.job import write_job_script @@ -76,6 +77,11 @@ def setup(self): # noqa: C901 use_vM_calving = (calving_method == 'von_mises') + # Thermal-forcing stream defaults; the 3D chunked case overrides them + # below to address the ..._.nc series with a $Y template. + tf_reference_time = '2000-01-01_00:00:00' + tf_filename_interval = 'none' + # --- Determine forcing file paths --- if scenario == 'ocx': ocx_forcing_path = section.get('ocx_forcing_path') @@ -213,14 +219,39 @@ def setup(self): # noqa: C901 else: tf_pattern = '*2dThermalForcing_*.nc' tf_search = os.path.join(ocean_dir, tf_pattern) - tf_list = glob.glob(tf_search) - if len(tf_list) == 1: + tf_list = sorted(glob.glob(tf_search)) + if len(tf_list) == 0: + sys.exit(f"ERROR: Expected at least 1 TF file at " + f"{tf_search}, found 0") + + if use_3d_thermal_forcing: + # 3D TF is written one file per N-year block, named by the + # block start year (..._.nc). Symlink the whole series + # and address it with a $Y filename_template plus a + # filename_interval so MALI advances across chunk files. + for tf_path in tf_list: + tf_base = os.path.split(tf_path)[-1] + os.symlink(tf_path, + os.path.join(self.work_dir, tf_base)) + start_years = sorted( + int(re.search(r'_(\d{4})\.nc$', + os.path.split(f)[-1]).group(1)) + for f in tf_list) + tf_first_year = start_years[0] + tf_reference_time = f"{tf_first_year:04d}-01-01_00:00:00" + sample = os.path.split(tf_list[0])[-1] + tf_fname = re.sub(r'_\d{4}\.nc$', '_$Y.nc', sample) + if len(start_years) > 1: + interval_years = start_years[1] - start_years[0] + tf_filename_interval = \ + f"{interval_years:04d}-00-00_00:00:00" + else: + if len(tf_list) != 1: + sys.exit(f"ERROR: Expected 1 TF file at {tf_search}, " + f"found {len(tf_list)}") tf_fname = os.path.split(tf_list[0])[-1] os.symlink(tf_list[0], os.path.join(self.work_dir, tf_fname)) - else: - sys.exit(f"ERROR: Expected 1 TF file at {tf_search}, " - f"found {len(tf_list)}") # --- Set up streams --- if scenario == 'ctrl': @@ -230,9 +261,10 @@ def setup(self): # noqa: C901 else: forcing_interval_monthly = '0000-01-00_00:00:00' forcing_interval_annual = '0001-00-00_00:00:00' - forcing_interval_TF = (forcing_interval_annual - if use_3d_thermal_forcing - else forcing_interval_monthly) + # GrIS thermal forcing (2D and 3D) is monthly. The annual + # cadence only applies to Antarctica, whose ISMIP7 3D forcing is + # provided annually. + forcing_interval_TF = forcing_interval_monthly stream_replacements = { 'input_file_init_cond': init_cond_fname if is_historical @@ -250,6 +282,8 @@ def setup(self): # noqa: C901 'forcing_interval_monthly': forcing_interval_monthly, 'forcing_interval_annual': forcing_interval_annual, 'forcing_interval_TF': forcing_interval_TF, + 'tf_reference_time': tf_reference_time, + 'tf_filename_interval': tf_filename_interval, 'use_3d_thermal_forcing': use_3d_thermal_forcing, } diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/streams.landice.template b/compass/landice/tests/ismip7_run/ismip7_gris/streams.landice.template index 32e7976999..649dd22d51 100644 --- a/compass/landice/tests/ismip7_run/ismip7_gris/streams.landice.template +++ b/compass/landice/tests/ismip7_run/ismip7_gris/streams.landice.template @@ -66,8 +66,9 @@ {% if use_3d_thermal_forcing %} From 9467ddc8630d4652205d7366f37a29d3c4cdb1a0 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Tue, 8 Sep 2026 19:15:09 -0700 Subject: [PATCH 13/17] Skip existing 3D forcing outputs on re-run instead of erroring Re-running the OCX step failed hard because the melt-params file (and forcing chunks) already existed with overwrite=false. Treat overwrite= false as "reuse existing, skip" so compass re-runs are idempotent and resumable (each chunk is written atomically via a .partial rename, so a present file is complete). overwrite=true still forces regeneration. --- .../ismip7_forcing/ocean_thermal/greenland_3d.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py index 7bf1df5ec2..f441f78a20 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py @@ -1140,9 +1140,11 @@ def write_melt_params( output = cfg.melt_params_file output.parent.mkdir(parents=True, exist_ok=True) if output.exists() and not cfg.overwrite: - raise FileExistsError( - f"Melt-parameters output exists and overwrite=false: {output}" + logger.info( + f"Melt-parameters file already exists; skipping (set " + f"overwrite=true to regenerate): {output}" ) + return temporary = output.with_name(output.name + ".partial") if temporary.exists(): raise FileExistsError( @@ -1206,6 +1208,7 @@ def write_melt_params( ) target.close() os.replace(temporary, output) + logger.info(f"Created {output}") except Exception: logger.warning( f"Melt-parameters output was not finalized; partial file, if " @@ -1300,9 +1303,11 @@ def _write_forcing_chunk( ) -> None: xr = require_xarray() if output_path.exists() and not cfg.overwrite: - raise FileExistsError( - f"Output exists and overwrite=false: {output_path}" + logger.info( + f"Forcing chunk already exists; skipping (set overwrite=true to " + f"regenerate): {output_path}" ) + return temporary = output_path.with_name(output_path.name + ".partial") if temporary.exists(): raise FileExistsError( @@ -1839,7 +1844,6 @@ def run(cfg: Config, logger, prepare_only: bool = False) -> None: calibrate_regional_delta_t(cfg, mesh, basin_ids, profiles) ) write_melt_params(cfg, basin_ids, delta_t, logger) - logger.info(f"Created {cfg.melt_params_file}") else: logger.info( f"scenario={cfg.scenario!r} is not OCX; reusing calibrated " From 92e63c4b3b31e09d28861280467a38bed0289039 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Tue, 8 Sep 2026 20:08:19 -0700 Subject: [PATCH 14/17] Write 3D forcing chunks via HDF5 then nccopy to CDF-2 (~5x faster) Writing NETCDF3_64BIT directly with an unlimited Time dimension and multiple record variables is several times slower because the classic writer interleaves each record and writes with a stride (benchmarked at ~143 s vs ~26 s for a 4 GB, 10-year monthly chunk). Write HDF5 first (variables stored contiguously) then convert to the CDF-2 format MALI/pnetcdf requires with nccopy, which preserves the unlimited Time dimension. The HDF5 scratch file is removed after conversion. --- .../ocean_thermal/greenland_3d.py | 47 ++++++++++++++----- 1 file changed, 36 insertions(+), 11 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py index f441f78a20..c8555d7427 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py @@ -25,6 +25,8 @@ import math import os import re +import shutil +import subprocess import sys import warnings from dataclasses import dataclass @@ -1290,6 +1292,25 @@ def chunk_output_path(output_file: Path, start_year: int) -> Path: return output_file.with_name(name) +def _convert_to_cdf2(hdf5_path: Path, cdf2_path: Path, logger) -> None: + """Convert an HDF5 (NETCDF4) file to CDF-2 (64-bit offset) with nccopy. + + MALI/pnetcdf reads classic formats only. Writing HDF5 then converting is + much faster than writing CDF-2 directly for large record variables, and + nccopy preserves the unlimited Time dimension MALI expects. + """ + nccopy = shutil.which("nccopy") + if nccopy is None: + raise RuntimeError( + "nccopy is required to convert the 3-D forcing to CDF-2 but was " + "not found on PATH." + ) + logger.info(f"Converting {hdf5_path.name} to CDF-2 (64-bit offset)") + subprocess.run( + [nccopy, "-k", "nc6", str(hdf5_path), str(cdf2_path)], check=True + ) + + def _write_forcing_chunk( cfg: Config, forcing_3d: Any, @@ -1314,6 +1335,7 @@ def _write_forcing_chunk( f"Partial output already exists: {temporary}. Remove or rename " "it after inspecting it." ) + hdf5_temporary = output_path.with_name(output_path.name + ".h5.partial") target = xr.Dataset( data_vars={ "xtime": xr.DataArray( @@ -1373,31 +1395,34 @@ def _write_forcing_chunk( try: logger.info( f"Writing {len(chunk_times)} monthly records " - f"({start_year}-{last_year}) to {temporary} with xarray " - "(NETCDF3_64BIT, float32)" + f"({start_year}-{last_year}) to {output_path.name}" ) - # Use the netCDF4 engine rather than scipy: scipy's classic/CDF-2 - # writer can produce a header that the netCDF-C library (ncdump, - # MALI) cannot read when a dataset mixes a record (unlimited-Time) - # variable with a scalar (0-D) variable. The netCDF4 engine writes a - # conformant CDF-2 file and still streams record-by-record. + # Writing NETCDF3_64BIT directly with an unlimited Time dimension and + # multiple record variables is several times slower because the + # classic writer interleaves each record and writes with a stride. + # Write HDF5 first (each variable stored contiguously, much faster), + # then convert to the CDF-2 format MALI/pnetcdf requires with nccopy. with dask.config.set(scheduler="single-threaded"): target.to_netcdf( - temporary, + hdf5_temporary, engine="netcdf4", - format="NETCDF3_64BIT", + format="NETCDF4", unlimited_dims=["Time"], encoding=encoding, ) target.close() + _convert_to_cdf2(hdf5_temporary, temporary, logger) os.replace(temporary, output_path) logger.info(f"Created {output_path}") except Exception: logger.warning( - f"Output was not finalized; partial file, if any, is at " - f"{temporary}" + f"Output was not finalized; partial files, if any, are at " + f"{hdf5_temporary} and {temporary}" ) raise + finally: + if hdf5_temporary.exists(): + hdf5_temporary.unlink() def write_output( From 43bf07081e6f31e38308e0e62cf5247d544302f5 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Tue, 8 Sep 2026 20:50:28 -0700 Subject: [PATCH 15/17] Warn when marine ice sits below the clamped anchor depth Where the seafloor is deeper than source_max_depth_m the anchor is clamped, so TF_3d matches TF_2d at source_max_depth_m rather than the true seafloor. If the regional profile still slopes at that depth, the per-cell offset (and thus the whole reconstructed column) is biased. Emit a runtime warning reporting the count, marine-ice area fraction, and deepest seafloor, and distinguish regions that still slope at the bottom (biased) from those that are flat there (benign), pointing to raising ocean_vertical_grid.bottom_m and source_max_depth_m together. --- .../ocean_thermal/greenland_3d.py | 69 +++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py index c8555d7427..5df08994aa 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py @@ -1861,9 +1861,78 @@ def print_summary( ) +# A near-floor TF gradient above this magnitude means clamping the anchor to +# source_max_depth_m biases a sloped column rather than a flat one. +_BOTTOM_GRADIENT_THRESHOLD_C_PER_M = 5.0e-4 + + +def warn_if_seafloor_below_max_depth( + cfg: Config, + mesh: MeshData, + basin_ids: np.ndarray, + profiles: RegionalProfiles, + logger, +) -> None: + """Warn when marine ice sits below the clamped anchor depth. + + Where the seafloor is deeper than ``source_max_depth_m`` the anchor is + clamped, so TF_3d matches TF_2d at ``source_max_depth_m`` rather than the + true seafloor. If the regional profile still has a vertical gradient at + that depth, the per-cell offset (and thus the whole column) is biased. + """ + max_depth = cfg.source_max_depth_m + marine_ice = ( + (mesh.thickness > cfg.minimum_ice_thickness_m) & (mesh.bed < 0.0) + ) + deep = marine_ice & (mesh.bed < -max_depth) + n_deep = int(np.count_nonzero(deep)) + if n_deep == 0: + return + marine_area = float(np.sum(mesh.area[marine_ice])) + deep_area = float(np.sum(mesh.area[deep])) + fraction = deep_area / marine_area if marine_area > 0.0 else 0.0 + deepest = float(-np.min(mesh.bed[deep])) + + # Near-floor gradient (degC/m) from the deepest two output levels tells + # us which affected regions still slope (biased) versus are flat (benign). + z = cfg.ocean_levels_m + tf = profiles.output_thermal_forcing_degC + dz = z[-1] - z[-2] + gradients = (tf[:, -1] - tf[:, -2]) / dz + sloped_regions = sorted({ + REGION_KEYS[int(region) - 1] + for region in np.unique(basin_ids[deep]) + if abs(gradients[int(region) - 1]) >= + _BOTTOM_GRADIENT_THRESHOLD_C_PER_M + }) + + message = ( + f"{n_deep} marine ice cells ({100.0 * fraction:.1f}% of marine-ice " + f"area; deepest {deepest:.0f} m) have seafloor below " + f"source_max_depth_m={max_depth:g} m, so their anchor is clamped to " + f"{max_depth:g} m and TF_3d matches TF_2d there rather than at the " + "true seafloor." + ) + if sloped_regions: + message += ( + " The regional profile still slopes at that depth in: " + f"{', '.join(sloped_regions)}; the clamp biases the whole " + "reconstructed column for those cells. Consider increasing " + "ocean_vertical_grid.bottom_m and source_max_depth_m together to " + "cover the deepest seafloor." + ) + else: + message += ( + " The regional profile is effectively flat at that depth, so the " + "clamp is benign here." + ) + logger.warning(message) + + def run(cfg: Config, logger, prepare_only: bool = False) -> None: mesh, basin_ids = load_mesh_and_basins(cfg) profiles = build_regional_profiles(cfg, mesh, basin_ids, logger) + warn_if_seafloor_below_max_depth(cfg, mesh, basin_ids, profiles, logger) if cfg.calibrate_delta_t: delta_t, achieved, calibration_monthly_tf = ( calibrate_regional_delta_t(cfg, mesh, basin_ids, profiles) From 58923f59d290d9059ac65d38664bd4e23bcb706d Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Tue, 8 Sep 2026 21:04:36 -0700 Subject: [PATCH 16/17] Impose seasonally-varying vertical structure on 3D TF Strong upper-ocean seasonality makes a single annual vertical profile misleading, so group the EN4 monthly profiles into seasons_per_year equal calendar blocks (default 4: Jan-Mar, Apr-Jun, Jul-Sep, Oct-Dec) and build one profile per season per region. Each output month uses the column shape of its season, anchored to that month's 2-D forcing at the seafloor (TF_3d = TF_2d at the anchor is preserved). Calibration and the diagnostics also use the per-month seasonal profile; the profile plots now draw the seasonal means (and season-colored monthly curves) instead of one annual mean, and regional_profiles.csv gains a season column. --- .../greenland_3d_tf_config.json | 7 +- .../ocean_thermal/greenland_3d.py | 338 +++++++++++++----- 2 files changed, 262 insertions(+), 83 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json b/compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json index 0abfe9af4b..fb82c25186 100644 --- a/compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json +++ b/compass/landice/tests/ismip7_forcing/greenland_3d_tf_config.json @@ -31,12 +31,13 @@ "minimum_ice_thickness_m": 0.0 }, "ocean_vertical_grid": { - "number_of_levels": 10, + "number_of_levels": 20, "surface_m": 0.0, - "bottom_m": -1000.0 + "bottom_m": -2500.0 }, - "source_max_depth_m": 1000.0, + "source_max_depth_m": 2500.0, "output_years_per_file": 10, + "seasons_per_year": 4, "forcing_2d_variable": "ismip6_2dThermalForcing", "physical_constants": { "rho_ice": 910.0, diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py index 5df08994aa..85e1a1bb4f 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py @@ -61,6 +61,27 @@ EN4_BIAS_CORRECTIONS = {"g10", "l09", "c13", "c14", "unknown"} DATE_RE = re.compile(r"(? int: + """Zero-based season for a 1-based calendar month. + + Months are partitioned into ``seasons_per_year`` equal blocks, e.g. with + four seasons: Jan-Mar=0, Apr-Jun=1, Jul-Sep=2, Oct-Dec=3. + """ + return (month - 1) // (12 // seasons_per_year) + + +def season_label(season: int, seasons_per_year: int) -> str: + months_per_season = 12 // seasons_per_year + start = season * months_per_season + end = start + months_per_season - 1 + return f"{MONTH_ABBR[start]}-{MONTH_ABBR[end]}" + def require_xarray(): """Import xarray late so mathematical unit tests need no NetCDF stack.""" @@ -122,6 +143,7 @@ class Config: forcing_variable: str overwrite: bool output_years_per_file: int + seasons_per_year: int @property def calibrate_delta_t(self) -> bool: @@ -258,6 +280,10 @@ def resolved(field_name, json_key, required=True, default=None): overrides.get("output_years_per_file") or raw.get("output_years_per_file", 10) ), + seasons_per_year=int( + overrides.get("seasons_per_year") or + raw.get("seasons_per_year", 4) + ), ) cfg.validate() return cfg @@ -288,6 +314,11 @@ def validate(self) -> None: raise ValueError("gamma0_m_per_yr must be positive") if self.output_years_per_file < 1: raise ValueError("output_years_per_file must be at least 1") + if self.seasons_per_year < 1 or 12 % self.seasons_per_year != 0: + raise ValueError( + "seasons_per_year must be a positive divisor of 12 " + "(1, 2, 3, 4, 6, or 12)" + ) if self.en4_bias_correction == "unknown": warnings.warn( "EN4 bias correction is unknown. Processing will continue " @@ -325,6 +356,13 @@ class RegionalProfiles: freezing_temperature_degC: np.ndarray thermal_forcing_degC: np.ndarray output_thermal_forcing_degC: np.ndarray + seasons_per_year: int + season_labels: tuple + season_of_month: np.ndarray + seasonal_temperature_degC: np.ndarray + seasonal_salinity: np.ndarray + seasonal_thermal_forcing_degC: np.ndarray + seasonal_output_thermal_forcing_degC: np.ndarray valid_gridpoint_counts: np.ndarray temperature_observation_influence: np.ndarray salinity_observation_influence: np.ndarray @@ -967,6 +1005,57 @@ def build_regional_profiles( ] ) + # Seasonally-varying vertical structure: group the monthly profiles into + # seasons_per_year equal calendar blocks and average within each block, so + # strong upper-ocean seasonality is preserved rather than smeared into one + # annual mean. + n_seasons = cfg.seasons_per_year + n_source = source_z.size + season_of_month = np.array( + [season_index(int(date[5:7]), n_seasons) for date in dates] + ) + seasonal_temp = np.full( + (n_seasons, len(REGION_NAMES), n_source), np.nan + ) + seasonal_sal = np.full_like(seasonal_temp, np.nan) + seasonal_tf = np.full_like(seasonal_temp, np.nan) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) + for season in range(n_seasons): + members = np.flatnonzero(season_of_month == season) + if members.size == 0: + raise ValueError( + f"No EN4 months fall in season " + f"{season_label(season, n_seasons)}; cannot build a " + "seasonal profile" + ) + seasonal_temp[season] = np.nanmean( + month_temp_array[members], axis=0 + ) + seasonal_sal[season] = np.nanmean( + month_sal_array[members], axis=0 + ) + seasonal_tf[season] = np.nanmean( + month_tf_array[members], axis=0 + ) + seasonal_output_tf = np.stack( + [ + np.vstack( + [ + interpolate_profile( + source_z, seasonal_tf[season, region], + cfg.ocean_levels_m + ) + for region in range(len(REGION_NAMES)) + ] + ) + for season in range(n_seasons) + ] + ) + season_labels = tuple( + season_label(season, n_seasons) for season in range(n_seasons) + ) + return RegionalProfiles( source_z_m=source_z, output_z_m=cfg.ocean_levels_m, @@ -979,6 +1068,13 @@ def build_regional_profiles( freezing_temperature_degC=mean_freeze, thermal_forcing_degC=mean_tf, output_thermal_forcing_degC=output_tf, + seasons_per_year=n_seasons, + season_labels=season_labels, + season_of_month=season_of_month, + seasonal_temperature_degC=seasonal_temp, + seasonal_salinity=seasonal_sal, + seasonal_thermal_forcing_degC=seasonal_tf, + seasonal_output_thermal_forcing_degC=seasonal_output_tf, valid_gridpoint_counts=np.asarray(monthly_counts), temperature_observation_influence=np.asarray(monthly_temp_influence), salinity_observation_influence=np.asarray(monthly_sal_influence), @@ -1002,6 +1098,15 @@ def year_from_xtime(value: str) -> int: return int(match.group(1)) +def month_from_xtime(value: str) -> int: + match = re.match(r"\s*\d{4}-(\d{2})", value) + if match is None: + raise ValueError( + f"Could not parse month from xtime value {value!r}" + ) + return int(match.group(1)) + + def calibration_time_indices( times: Sequence[str], start_year: int, end_year: int ) -> np.ndarray: @@ -1045,13 +1150,25 @@ def calibrate_regional_delta_t( cfg.minimum_ice_thickness_m, ) anchor = np.clip(mesh.bed, -cfg.source_max_depth_m, 0.0) - base_at_anchor = profiles_at_cell_depths( - profiles.output_thermal_forcing_degC, cfg.ocean_levels_m, basin_ids, - anchor + # Seasonal base profiles at the effective seafloor (anchor) and the ice + # draft; the offset uses whichever season each calibration month is in. + seasonal_out = profiles.seasonal_output_thermal_forcing_degC + n_seasons = seasonal_out.shape[0] + base_at_anchor_seasonal = np.stack( + [ + profiles_at_cell_depths( + seasonal_out[season], cfg.ocean_levels_m, basin_ids, anchor + ) + for season in range(n_seasons) + ] ) - base_at_draft = profiles_at_cell_depths( - profiles.output_thermal_forcing_degC, cfg.ocean_levels_m, basin_ids, - draft + base_at_draft_seasonal = np.stack( + [ + profiles_at_cell_depths( + seasonal_out[season], cfg.ocean_levels_m, basin_ids, draft + ) + for season in range(n_seasons) + ] ) floating_counts = np.asarray( @@ -1086,8 +1203,11 @@ def calibrate_regional_delta_t( forcing = _array_with_nan( forcing_var.isel(Time=int(time_index)).values ) - offset = forcing - base_at_anchor - tf_draft = base_at_draft + offset + season = season_index( + month_from_xtime(times[time_index]), n_seasons + ) + offset = forcing - base_at_anchor_seasonal[season] + tf_draft = base_at_draft_seasonal[season] + offset for region in range(len(REGION_NAMES)): if region in empty_regions: continue @@ -1355,6 +1475,7 @@ def _write_forcing_chunk( "title": "Regional three-dimensional Greenland ocean thermal " "forcing for MALI", "source_2d_forcing": str(cfg.forcing_2d_file), + "seasonal_vertical_structure": f"{cfg.seasons_per_year} seasons", "melt_params_file": str(cfg.melt_params_file), "source_en4_version": cfg.en4_version, "source_en4_bias_correction": cfg.en4_bias_correction, @@ -1436,6 +1557,7 @@ def write_output( xr = require_xarray() try: import dask + import dask.array as darray except ImportError as exc: # pragma: no cover - environment dependent raise RuntimeError( "Dask is required to construct and stream the multi-gigabyte " @@ -1444,11 +1566,24 @@ def write_output( cfg.output_file.parent.mkdir(parents=True, exist_ok=True) anchor = np.clip(mesh.bed, -cfg.source_max_depth_m, 0.0) - base_at_anchor = profiles_at_cell_depths( - profiles.output_thermal_forcing_degC, cfg.ocean_levels_m, basin_ids, - anchor + # Per-season vertical structure of the column relative to the anchor. The + # 2-D forcing sets the value at the anchor each month, so this delta plus + # the 2-D forcing reconstructs the full column; the season used varies + # with each output month. + seasonal_out = profiles.seasonal_output_thermal_forcing_degC + n_seasons = seasonal_out.shape[0] + base_by_cell_seasonal = seasonal_out[:, basin_ids - 1, :] + base_at_anchor_seasonal = np.stack( + [ + profiles_at_cell_depths( + seasonal_out[season], cfg.ocean_levels_m, basin_ids, anchor + ) + for season in range(n_seasons) + ] ) - base_by_cell = profiles.output_thermal_forcing_degC[basin_ids - 1, :] + delta_seasonal = ( + base_by_cell_seasonal - base_at_anchor_seasonal[:, :, None] + ).astype(np.float32) with xr.open_dataset( cfg.forcing_2d_file, @@ -1460,6 +1595,9 @@ def write_output( forcing_var = validate_forcing_schema(source, cfg, mesh.bed.size) times = forcing_times(source) years = np.asarray([year_from_xtime(value) for value in times]) + season_of_time = np.asarray( + [season_index(month_from_xtime(t), n_seasons) for t in times] + ) # Fail before initiating a multi-gigabyte write if any source value # is absent. This reduction stays lazy until compute() and does not # load the full forcing array into memory. @@ -1473,23 +1611,20 @@ def write_output( "building MALI forcing" ) - cell_coord = np.arange(mesh.bed.size, dtype=np.int32) - layer_coord = np.arange(cfg.ocean_levels_m.size, dtype=np.int32) - base_cells = xr.DataArray( - base_by_cell.astype(np.float32), - dims=("nCells", "nISMIP6OceanLayers"), - coords={ - "nCells": cell_coord, - "nISMIP6OceanLayers": layer_coord, - }, + delta_by_season = xr.DataArray( + darray.from_array( + delta_seasonal, + chunks=(1, mesh.bed.size, cfg.ocean_levels_m.size), + ), + dims=("season", "nCells", "nISMIP6OceanLayers"), ) - anchor_cells = xr.DataArray( - base_at_anchor.astype(np.float32), - dims=("nCells",), - coords={"nCells": cell_coord}, + # Gather each month's seasonal column shape; stays lazy so the write + # streams record by record. + delta_by_time = delta_by_season.isel( + season=xr.DataArray(season_of_time, dims="Time") ) forcing_3d = ( - base_cells + (forcing_var.astype(np.float32) - anchor_cells) + forcing_var.astype(np.float32) + delta_by_time ).transpose("Time", "nCells", "nISMIP6OceanLayers").assign_attrs( units="C", long_name="3D thermal forcing for nonlocal ISMIP6 ice-shelf " @@ -1539,6 +1674,7 @@ def write_diagnostics( writer = csv.writer(handle) writer.writerow( [ + "season", "region_id", "region_key", "region_name", @@ -1549,22 +1685,39 @@ def write_diagnostics( "thermal_forcing_degC", ] ) - for region in range(len(REGION_NAMES)): - for level, z in enumerate(profiles.source_z_m): - writer.writerow( - [ - region + 1, - REGION_KEYS[region], - REGION_NAMES[region], - float(z), - float(profiles.temperature_degC[region, level]), - float(profiles.salinity[region, level]), - float( - profiles.freezing_temperature_degC[region, level] - ), - float(profiles.thermal_forcing_degC[region, level]), - ] - ) + for season in range(profiles.seasons_per_year): + season_freeze = freezing_temperature( + profiles.seasonal_salinity[season], + profiles.source_z_m[None, :], + cfg.freezing_a, cfg.freezing_b, cfg.freezing_c, + ) + for region in range(len(REGION_NAMES)): + for level, z in enumerate(profiles.source_z_m): + writer.writerow( + [ + profiles.season_labels[season], + region + 1, + REGION_KEYS[region], + REGION_NAMES[region], + float(z), + float( + profiles.seasonal_temperature_degC[ + season, region, level + ] + ), + float( + profiles.seasonal_salinity[ + season, region, level + ] + ), + float(season_freeze[region, level]), + float( + profiles.seasonal_thermal_forcing_degC[ + season, region, level + ] + ), + ] + ) with (directory / "en4_coverage.csv").open( "w", newline="", encoding="utf-8" @@ -1712,14 +1865,18 @@ def write_diagnostics( .isel(Time=time_index) .values ) + rep_season = season_index( + month_from_xtime(times[time_index]), profiles.seasons_per_year + ) + seasonal_out = profiles.seasonal_output_thermal_forcing_degC[rep_season] anchor = np.clip(mesh.bed, -cfg.source_max_depth_m, 0.0) base_at_anchor = profiles_at_cell_depths( - profiles.output_thermal_forcing_degC, + seasonal_out, cfg.ocean_levels_m, basin_ids, anchor, ) - base_by_cell = profiles.output_thermal_forcing_degC[basin_ids - 1] + base_by_cell = seasonal_out[basin_ids - 1] offset = forcing_2d - base_at_anchor forcing_3d = base_by_cell + offset[:, None] reconstructed_anchor = base_at_anchor + offset @@ -1732,6 +1889,7 @@ def write_diagnostics( json.dump( { "representative_time": times[time_index], + "representative_season": profiles.season_labels[rep_season], "maximum_absolute_anchor_error_degC": max_anchor_error, "effective_anchor_depth_range_m": [ float(np.nanmin(anchor)), @@ -1774,7 +1932,10 @@ def write_diagnostics( ax.set_xlabel("Longitude") ax.set_ylabel("Latitude") fig.colorbar(scatter, ax=ax, label="Thermal forcing (°C)") - fig.suptitle(f"Translated 3-D forcing: {times[time_index]}") + fig.suptitle( + f"Translated 3-D forcing: {times[time_index]} " + f"(season {profiles.season_labels[rep_season]})" + ) fig.tight_layout() fig.savefig( directory / "thermal_forcing_at_representative_ocean_levels.png", @@ -1782,50 +1943,67 @@ def write_diagnostics( ) plt.close(fig) + season_colors = plt.get_cmap("turbo")( + np.linspace(0.05, 0.95, profiles.seasons_per_year) + ) for region in range(len(REGION_NAMES)): fig, axes = plt.subplots(1, 2, figsize=(10, 7), sharey=True) - for monthly in profiles.monthly_temperature_degC[:, region, :]: + # Faint monthly profiles, colored by season, behind the seasonal + # means so the seasonality that motivates the seasonal structure is + # visible. + for entry, monthly in enumerate( + profiles.monthly_temperature_degC[:, region, :] + ): axes[0].plot( - monthly, profiles.source_z_m, color="tab:blue", alpha=0.08, - linewidth=0.5 + monthly, profiles.source_z_m, + color=season_colors[profiles.season_of_month[entry]], + alpha=0.06, linewidth=0.5 + ) + for season in range(profiles.seasons_per_year): + axes[0].plot( + profiles.seasonal_temperature_degC[season, region], + profiles.source_z_m, + color=season_colors[season], + linewidth=2, + label=profiles.season_labels[season], + ) + selected_temp = interpolate_profile( + profiles.source_z_m, + profiles.seasonal_temperature_degC[season, region], + profiles.output_z_m, + ) + axes[0].scatter( + selected_temp, profiles.output_z_m, + color=season_colors[season], marker="*", zorder=3, ) - axes[0].plot( - profiles.temperature_degC[region], - profiles.source_z_m, - color="black", - linewidth=2, - label="Monthly mean climatology", - ) - selected_temp = interpolate_profile( - profiles.source_z_m, profiles.temperature_degC[region], - profiles.output_z_m - ) - axes[0].scatter( - selected_temp, profiles.output_z_m, color="black", marker="*", - zorder=3, label="MALI levels" - ) axes[0].set_xlabel("Potential temperature (°C)") axes[0].set_ylabel("Elevation (m)") - axes[0].legend(loc="best", fontsize=8) + axes[0].legend(loc="best", fontsize=8, title="season / MALI levels") - for monthly in profiles.monthly_thermal_forcing_degC[:, region, :]: + for entry, monthly in enumerate( + profiles.monthly_thermal_forcing_degC[:, region, :] + ): axes[1].plot( - monthly, profiles.source_z_m, color="tab:red", alpha=0.08, - linewidth=0.5 + monthly, profiles.source_z_m, + color=season_colors[profiles.season_of_month[entry]], + alpha=0.06, linewidth=0.5 + ) + for season in range(profiles.seasons_per_year): + axes[1].plot( + profiles.seasonal_thermal_forcing_degC[season, region], + profiles.source_z_m, + color=season_colors[season], + linewidth=2, + ) + axes[1].scatter( + profiles.seasonal_output_thermal_forcing_degC[ + season, region + ], + profiles.output_z_m, + color=season_colors[season], + marker="*", + zorder=3, ) - axes[1].plot( - profiles.thermal_forcing_degC[region], - profiles.source_z_m, - color="black", - linewidth=2, - ) - axes[1].scatter( - profiles.output_thermal_forcing_degC[region], - profiles.output_z_m, - color="black", - marker="*", - zorder=3, - ) axes[1].set_xlabel("Thermal forcing (°C)") for ax in axes: ax.grid(alpha=0.3) From ae050564b72eb1971ffe0fb4b257c5d1e8b4fb8c Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Tue, 8 Sep 2026 21:25:05 -0700 Subject: [PATCH 17/17] Fix longitude convention mismatch in diagnostic maps mapped_lons is wrapped to [-180, 180] but mesh.lon_deg is in [0, 360], so the EN4 points and the mesh landed on opposite sides of the region-assignment map (the assignment itself is correct; the KD-tree uses convention-independent xyz). Wrap the mesh longitudes to [-180, 180] in the region-assignment and representative-field maps so the two overlay. --- .../tests/ismip7_forcing/ocean_thermal/greenland_3d.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py index 85e1a1bb4f..039cb22cc4 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/greenland_3d.py @@ -1823,8 +1823,11 @@ def write_diagnostics( colors = plt.get_cmap("tab10")(np.arange(len(REGION_NAMES))) fig, ax = plt.subplots(figsize=(8, 8)) mesh_stride = max(1, mesh.lat_deg.size // 100_000) + # mapped_lons is wrapped to [-180, 180]; put the mesh on the same + # convention so the two overlay instead of landing on opposite sides. + mesh_lon_plot = (mesh.lon_deg + 180.0) % 360.0 - 180.0 ax.scatter( - mesh.lon_deg[::mesh_stride], + mesh_lon_plot[::mesh_stride], mesh.lat_deg[::mesh_stride], c=colors[basin_ids[::mesh_stride] - 1], s=0.15, @@ -1912,6 +1915,7 @@ def write_diagnostics( nrows, ncols, figsize=(7 * ncols, 6 * nrows), squeeze=False ) plot_stride = max(1, mesh.lat_deg.size // 150_000) + mesh_lon_plot = (mesh.lon_deg + 180.0) % 360.0 - 180.0 color_limits = np.nanpercentile(forcing_3d[::plot_stride], [2.0, 98.0]) for plot_index, ax in enumerate(axes.ravel()): if plot_index >= plotted_levels.size: @@ -1919,7 +1923,7 @@ def write_diagnostics( continue level = int(plotted_levels[plot_index]) scatter = ax.scatter( - mesh.lon_deg[::plot_stride], + mesh_lon_plot[::plot_stride], mesh.lat_deg[::plot_stride], c=forcing_3d[::plot_stride, level], s=0.5,