diff --git a/compass/landice/tests/ismip7_forcing/__init__.py b/compass/landice/tests/ismip7_forcing/__init__.py index 3f621a4c37..cb448350e3 100644 --- a/compass/landice/tests/ismip7_forcing/__init__.py +++ b/compass/landice/tests/ismip7_forcing/__init__.py @@ -1,4 +1,5 @@ from compass.landice.tests.ismip7_forcing.atmosphere import Atmosphere +from compass.landice.tests.ismip7_forcing.fracture import Fracture from compass.landice.tests.ismip7_forcing.ocean_thermal import OceanThermal from compass.testgroup import TestGroup @@ -22,3 +23,4 @@ def __init__(self, mpas_core): self.add_test_case(Atmosphere(test_group=self)) self.add_test_case(OceanThermal(test_group=self)) + self.add_test_case(Fracture(test_group=self)) diff --git a/compass/landice/tests/ismip7_forcing/fracture/__init__.py b/compass/landice/tests/ismip7_forcing/fracture/__init__.py new file mode 100644 index 0000000000..b6a7e63c6a --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/fracture/__init__.py @@ -0,0 +1,51 @@ +from compass.landice.tests.ismip7_forcing.configure import ( + configure as configure_testgroup, +) +from compass.landice.tests.ismip7_forcing.fracture.process_excess_melt import ( + ProcessExcessMelt, +) +from compass.landice.tests.ismip7_forcing.fracture.process_lake_properties import ( # noqa: E501 + ProcessLakeProperties, +) +from compass.landice.tests.ismip7_forcing.fracture.process_shelf_collapse import ( # noqa: E501 + ProcessShelfCollapse, +) +from compass.testcase import TestCase + + +class Fracture(TestCase): + """ + A test case for processing ISMIP7 fracture forcing data. + Implements the surface-melt-driven ice shelf collapse pathways: + + * Path A (``process_excess_melt``): excess meltwater after firn air + content depletion. + * Path B (``process_lake_properties``): supraglacial lake mean depth + and area fraction from Grau et al. (2025). + * Path C (``process_shelf_collapse``): the ice shelf collapse mask, in + which a floating grid cell is flagged as collapsed when excess + meltwater exceeds 72.5 mm/yr for 10 consecutive years. + """ + + def __init__(self, test_group): + """ + Create the test case + + Parameters + ---------- + test_group : compass.landice.tests.ismip7_forcing.Ismip7Forcing + The test group that this test case belongs to + """ + name = "fracture" + subdir = name + super().__init__(test_group=test_group, name=name, subdir=subdir) + + self.add_step(ProcessExcessMelt(test_case=self)) + self.add_step(ProcessLakeProperties(test_case=self)) + self.add_step(ProcessShelfCollapse(test_case=self)) + + def configure(self): + """ + Configures test case + """ + configure_testgroup(config=self.config) diff --git a/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py b/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py new file mode 100644 index 0000000000..fbf3f75b83 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py @@ -0,0 +1,267 @@ +import glob +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 compass.landice.tests.ismip7_forcing.create_mapfile import ( + build_mapping_file, +) +from compass.landice.tests.ismip7_forcing.fracture.remap_utils import ( + add_xtime_and_write, + extrapolate_source, + open_rename_and_trim, +) +from compass.step import Step + + +class ProcessExcessMelt(Step): + """ + A step for processing the ISMIP7 excess meltwater field (Path A). + Remaps the annual excess melt (melt + rain after firn air content + depletion) from the ISMIP7 polar stereographic grid to the MALI + unstructured mesh. + + The excess melt file lacks ``x``/``y`` coordinate variables and its + array is flipped along the y axis relative to the other fracture + files (it was produced with CDO). This step reconstructs a source + grid with ``x``/``y`` coordinates and the correct orientation before + remapping. + """ + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_forcing.fracture.Fracture + The test case this step belongs to + """ + super().__init__(test_case=test_case, name="process_excess_melt") + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config["ismip7"] + base_path_mali = section.get("base_path_mali") + mali_mesh_file = section.get("mali_mesh_file") + + self.add_input_file(filename=mali_mesh_file, + target=os.path.join(base_path_mali, + mali_mesh_file)) + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + + section = config["ismip7"] + base_path_ismip7 = section.get("base_path_ismip7") + mali_mesh_name = section.get("mali_mesh_name") + mali_mesh_file = section.get("mali_mesh_file") + model = section.get("model") + scenario = section.get("scenario") + output_base_path = section.get("output_base_path") + ice_sheet = section.get("ice_sheet") + + section = config["ismip7_fracture"] + method_remap = section.get("method_remap_excess_melt") + version = section.get("version") + start_year = section.getint("start_year") + end_year = section.getint("end_year") + + # Skip this pathway if no remapping method is requested + if method_remap.lower() == "none": + logger.info("method_remap_excess_melt is None; skipping excess " + "melt (Path A) processing.") + return + + # Discover the excess melt file + input_path = os.path.join(base_path_ismip7, "fracture", version) + file_pattern = "excess_melt_*.nc" + all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) + + if not all_files: + raise FileNotFoundError( + f"No excess melt file found matching pattern:\n" + f" {os.path.join(input_path, file_pattern)}") + if len(all_files) > 1: + raise ValueError( + f"Expected a single excess melt file but found " + f"{len(all_files)}:\n " + "\n ".join(all_files)) + + input_file = all_files[0] + basename = os.path.basename(input_file) + logger.info(f"Processing excess melt: {basename}") + + # Build a source file with x/y coordinates and correct orientation + gridded_file = f"gridded_{basename}" + self._prepare_source_grid(input_file, input_path, gridded_file, + logger) + + # Build mapping file. Excess melt is a flux, so conservative + # remapping is appropriate by default. + mapping_file = (f"map_ismip7_{ice_sheet}_fracture_to_" + f"{mali_mesh_name}_{method_remap}.nc") + + if not os.path.exists(mapping_file): + logger.info("Building mapping file for the excess melt grid...") + build_mapping_file(config, logger, + gridded_file, mapping_file, + mali_mesh_file=mali_mesh_file, + method_remap=method_remap) + + # Extrapolate fill values on the source grid before remapping so + # they don't pollute neighboring cells during interpolation + extrap_file = f"extrap_{basename}" + extrapolate_source(gridded_file, extrap_file, "excess_melt", logger) + + # Remap the excess melt onto the MALI mesh + remapped_file = f"remapped_{basename}" + logger.info(f"Remapping: {basename}") + args = ["ncremap", + "-i", extrap_file, + "-o", remapped_file, + "-m", mapping_file, + "-v", "excess_melt"] + check_call(args, logger=logger) + + # Rename to MALI conventions + logger.info("Renaming variables to MALI conventions...") + output_file = f"{mali_mesh_name}_{basename}" + self._rename_to_mali_vars(remapped_file, output_file, + start_year, end_year) + + # Clean up temporary files + for f in [gridded_file, extrap_file, remapped_file]: + if os.path.exists(f): + os.remove(f) + + # Place output in the appropriate directory + output_path = os.path.join(output_base_path, "excess_melt", + f"{model}_{scenario}") + if not os.path.exists(output_path): + os.makedirs(output_path) + + dst = os.path.join(output_path, output_file) + shutil.copy(output_file, dst) + + logger.info(f"Done. Output: {dst}") + + def _prepare_source_grid(self, input_file, input_path, output_file, + logger): + """ + Build a source file for the excess melt data with ``x``/``y`` + coordinate variables and the standard ISMIP7 orientation + (south-to-north). The excess melt file only has 2D ``lat``/``lon`` + and its array is flipped along the y axis relative to the other + fracture files, so ``x``/``y`` are borrowed from a sibling + fracture file and the data are flipped to match. + + Parameters + ---------- + input_file : str + Path to the excess melt file + + input_path : str + Directory containing the fracture forcing files + + output_file : str + Path to write the reconstructed source file + + logger : logging.Logger + Logger for status messages + """ + # Find a sibling fracture file that has x/y coordinate variables + grid_donor = None + for pattern in ["lake_properties_*.nc", + "ice_shelf_collapse_mask_*.nc"]: + for candidate in sorted(glob.glob(os.path.join(input_path, + pattern))): + with xr.open_dataset(candidate) as dtest: + if "x" in dtest.variables and "y" in dtest.variables: + grid_donor = candidate + break + if grid_donor is not None: + break + + if grid_donor is None: + raise FileNotFoundError( + "Could not find a sibling fracture file with x/y " + "coordinate variables to define the excess melt grid.") + + logger.info(f"Using grid from sibling file: " + f"{os.path.basename(grid_donor)}") + + with xr.open_dataset(grid_donor) as donor: + x = donor["x"].values + y = donor["y"].values + donor_lat = donor["lat"].values + + ds = xr.open_dataset(input_file, decode_times=False) + + # Flip along the y axis to match the sibling grid orientation + flipped_lat = ds["lat"].values[::-1, :] + if np.nanmax(np.abs(flipped_lat - donor_lat)) > 1.0e-3: + raise ValueError( + "Excess melt grid does not match the sibling grid after a " + "y-axis flip; orientation cannot be determined " + "automatically.") + + excess = ds["excess_melt"].values[:, ::-1, :] + years = ds["year"].values.astype(int) + + out = xr.Dataset() + out["x"] = ("x", x) + out["y"] = ("y", y) + out["time"] = ("time", years) + out["excess_melt"] = (("time", "y", "x"), excess) + out["excess_melt"].attrs = dict(ds["excess_melt"].attrs) + + write_netcdf(out, output_file) + ds.close() + + def _rename_to_mali_vars(self, remapped_file, output_file, + start_year, end_year): + """ + Rename dimensions/variables of the remapped excess melt to MALI + conventions, restrict to the requested year range, and add an + ``xtime`` variable. + + Parameters + ---------- + remapped_file : str + Excess melt remapped onto the MALI mesh + + output_file : str + Output file with MALI variable/dimension names + + start_year : int + First year (inclusive) to retain + + end_year : int + Last year (inclusive) to retain + """ + ds, years = open_rename_and_trim( + remapped_file, {"excess_melt": "ismip7ExcessMelt"}, + start_year, end_year) + + # Convert from mm w.e. yr-1 to SI units of kg m-2 s-1 + # (1 mm w.e. = 1 kg m-2; 365-day year, as used elsewhere in MALI) + seconds_per_year = 365.0 * 24.0 * 3600.0 + ds["ismip7ExcessMelt"] = ds["ismip7ExcessMelt"] / seconds_per_year + + ds["ismip7ExcessMelt"].attrs = { + "long_name": "excess meltwater after firn air content depletion", + "units": "kg m-2 s-1", + } + + add_xtime_and_write(ds, years, output_file) diff --git a/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py b/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py new file mode 100644 index 0000000000..10df254c69 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py @@ -0,0 +1,195 @@ +import glob +import os +import shutil + +from mpas_tools.logging import check_call + +from compass.landice.tests.ismip7_forcing.create_mapfile import ( + build_mapping_file, +) +from compass.landice.tests.ismip7_forcing.fracture.remap_utils import ( + add_xtime_and_write, + extrapolate_source, + open_rename_and_trim, +) +from compass.step import Step + + +class ProcessLakeProperties(Step): + """ + A step for processing the ISMIP7 supraglacial lake properties (Path B). + Remaps the annual mean lake depth and lake area fraction (from the + Grau et al. (2025) parameterization) from the ISMIP7 polar + stereographic grid to the MALI unstructured mesh. + """ + + # Source variable name -> MALI output variable name and attributes + _variables = { + "lake_depth": { + "mali_name": "ismip7LakeDepth", + "long_name": "mean supraglacial lake depth", + "units": "m", + }, + "fraction_lake_area": { + "mali_name": "ismip7LakeAreaFraction", + "long_name": "supraglacial lake area fraction", + "units": "1", + }, + } + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_forcing.fracture.Fracture + The test case this step belongs to + """ + super().__init__(test_case=test_case, name="process_lake_properties") + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config["ismip7"] + base_path_mali = section.get("base_path_mali") + mali_mesh_file = section.get("mali_mesh_file") + + self.add_input_file(filename=mali_mesh_file, + target=os.path.join(base_path_mali, + mali_mesh_file)) + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + + section = config["ismip7"] + base_path_ismip7 = section.get("base_path_ismip7") + mali_mesh_name = section.get("mali_mesh_name") + mali_mesh_file = section.get("mali_mesh_file") + model = section.get("model") + scenario = section.get("scenario") + output_base_path = section.get("output_base_path") + ice_sheet = section.get("ice_sheet") + + section = config["ismip7_fracture"] + method_remap = section.get("method_remap_lake_properties") + version = section.get("version") + start_year = section.getint("start_year") + end_year = section.getint("end_year") + + # Skip this pathway if no remapping method is requested + if method_remap.lower() == "none": + logger.info("method_remap_lake_properties is None; skipping lake " + "properties (Path B) processing.") + return + + # Discover the lake properties file + input_path = os.path.join(base_path_ismip7, "fracture", version) + file_pattern = "lake_properties_*.nc" + all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) + + if not all_files: + raise FileNotFoundError( + f"No lake properties file found matching pattern:\n" + f" {os.path.join(input_path, file_pattern)}") + if len(all_files) > 1: + raise ValueError( + f"Expected a single lake properties file but found " + f"{len(all_files)}:\n " + "\n ".join(all_files)) + + input_file = all_files[0] + basename = os.path.basename(input_file) + logger.info(f"Processing lake properties: {basename}") + + # Build mapping file. Lake properties are continuous fields, so + # bilinear remapping is appropriate by default. + mapping_file = (f"map_ismip7_{ice_sheet}_fracture_to_" + f"{mali_mesh_name}_{method_remap}.nc") + + if not os.path.exists(mapping_file): + logger.info("Building mapping file for the lake properties " + "grid...") + build_mapping_file(config, logger, + input_file, mapping_file, + mali_mesh_file=mali_mesh_file, + method_remap=method_remap) + + # Extrapolate fill values on the source grid before remapping so + # they don't pollute neighboring cells during interpolation + extrap_file = f"extrap_{basename}" + extrapolate_source(input_file, extrap_file, + list(self._variables.keys()), logger) + + # Remap both lake property variables onto the MALI mesh + remapped_file = f"remapped_{basename}" + logger.info(f"Remapping: {basename}") + args = ["ncremap", + "-i", extrap_file, + "-o", remapped_file, + "-m", mapping_file, + "-v", ",".join(self._variables.keys())] + check_call(args, logger=logger) + + # Rename to MALI conventions + logger.info("Renaming variables to MALI conventions...") + output_file = f"{mali_mesh_name}_{basename}" + self._rename_to_mali_vars(remapped_file, output_file, + start_year, end_year) + + # Clean up temporary files + for f in [extrap_file, remapped_file]: + if os.path.exists(f): + os.remove(f) + + # Place output in the appropriate directory + output_path = os.path.join(output_base_path, "lake_properties", + f"{model}_{scenario}") + if not os.path.exists(output_path): + os.makedirs(output_path) + + dst = os.path.join(output_path, output_file) + shutil.copy(output_file, dst) + + logger.info(f"Done. Output: {dst}") + + def _rename_to_mali_vars(self, remapped_file, output_file, + start_year, end_year): + """ + Rename dimensions/variables of the remapped lake properties to MALI + conventions, restrict to the requested year range, and add an + ``xtime`` variable. + + Parameters + ---------- + remapped_file : str + Lake properties remapped onto the MALI mesh + + output_file : str + Output file with MALI variable/dimension names + + start_year : int + First year (inclusive) to retain + + end_year : int + Last year (inclusive) to retain + """ + rename_vars = {src: info["mali_name"] + for src, info in self._variables.items()} + ds, years = open_rename_and_trim(remapped_file, rename_vars, + start_year, end_year) + + for info in self._variables.values(): + mali_name = info["mali_name"] + if mali_name in ds: + ds[mali_name].attrs = { + "long_name": info["long_name"], + "units": info["units"], + } + + add_xtime_and_write(ds, years, output_file) diff --git a/compass/landice/tests/ismip7_forcing/fracture/process_shelf_collapse.py b/compass/landice/tests/ismip7_forcing/fracture/process_shelf_collapse.py new file mode 100644 index 0000000000..e9a8e1cfb5 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/fracture/process_shelf_collapse.py @@ -0,0 +1,174 @@ +import glob +import os +import shutil + +from mpas_tools.logging import check_call + +from compass.landice.tests.ismip7_forcing.create_mapfile import ( + build_mapping_file, +) +from compass.landice.tests.ismip7_forcing.fracture.remap_utils import ( + add_xtime_and_write, + open_rename_and_trim, +) +from compass.step import Step + + +class ProcessShelfCollapse(Step): + """ + A step for processing the ISMIP7 ice shelf collapse mask (Path C). + Remaps the annual 0/1 collapse mask from the ISMIP7 polar + stereographic grid to the MALI unstructured mesh and renames + variables to MALI conventions. The mask flags floating grid cells + that collapse once excess meltwater (after firn air content + depletion) exceeds 72.5 mm/yr for 10 consecutive years. + """ + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_forcing.fracture.Fracture + The test case this step belongs to + """ + super().__init__(test_case=test_case, name="process_shelf_collapse") + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config["ismip7"] + base_path_mali = section.get("base_path_mali") + mali_mesh_file = section.get("mali_mesh_file") + + self.add_input_file(filename=mali_mesh_file, + target=os.path.join(base_path_mali, + mali_mesh_file)) + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + + section = config["ismip7"] + base_path_ismip7 = section.get("base_path_ismip7") + mali_mesh_name = section.get("mali_mesh_name") + mali_mesh_file = section.get("mali_mesh_file") + model = section.get("model") + scenario = section.get("scenario") + output_base_path = section.get("output_base_path") + ice_sheet = section.get("ice_sheet") + + section = config["ismip7_fracture"] + method_remap = section.get("method_remap_shelf_collapse") + version = section.get("version") + start_year = section.getint("start_year") + end_year = section.getint("end_year") + + # Skip this pathway if no remapping method is requested + if method_remap.lower() == "none": + logger.info("method_remap_shelf_collapse is None; skipping ice " + "shelf collapse mask (Path C) processing.") + return + + # Discover the ice shelf collapse mask file + input_path = os.path.join(base_path_ismip7, "fracture", version) + file_pattern = "ice_shelf_collapse_mask_*.nc" + all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) + + if not all_files: + raise FileNotFoundError( + f"No ice shelf collapse mask file found matching pattern:\n" + f" {os.path.join(input_path, file_pattern)}") + if len(all_files) > 1: + raise ValueError( + f"Expected a single ice shelf collapse mask file but found " + f"{len(all_files)}:\n " + "\n ".join(all_files)) + + input_file = all_files[0] + basename = os.path.basename(input_file) + logger.info(f"Processing ice shelf collapse mask: {basename}") + + # Build mapping file. neareststod preserves the 0/1 mask values. + mapping_file = (f"map_ismip7_{ice_sheet}_fracture_to_" + f"{mali_mesh_name}_{method_remap}.nc") + + if not os.path.exists(mapping_file): + logger.info("Building mapping file for the collapse mask grid...") + build_mapping_file(config, logger, + input_file, mapping_file, + mali_mesh_file=mali_mesh_file, + method_remap=method_remap) + + # Remap the collapse mask onto the MALI mesh + remapped_file = f"remapped_{basename}" + if not os.path.exists(remapped_file): + logger.info(f"Remapping: {basename}") + args = ["ncremap", + "-i", input_file, + "-o", remapped_file, + "-m", mapping_file, + "-v", "mask"] + check_call(args, logger=logger) + + # Combine time slice and rename to MALI conventions + logger.info("Renaming variables to MALI conventions...") + output_file = f"{mali_mesh_name}_{basename}" + self._rename_to_mali_vars(remapped_file, output_file, + start_year, end_year) + + # Clean up temporary remapped file + if os.path.exists(remapped_file): + os.remove(remapped_file) + + # Place output in the appropriate directory + output_path = os.path.join(output_base_path, "shelf_collapse", + f"{model}_{scenario}") + if not os.path.exists(output_path): + os.makedirs(output_path) + + dst = os.path.join(output_path, output_file) + shutil.copy(output_file, dst) + + logger.info(f"Done. Output: {dst}") + + def _rename_to_mali_vars(self, remapped_file, output_file, + start_year, end_year): + """ + Rename dimensions/variables of the remapped collapse mask to MALI + conventions, restrict to the requested year range, round the mask + to 0/1, and add an ``xtime`` variable. + + Parameters + ---------- + remapped_file : str + Collapse mask remapped onto the MALI mesh + + output_file : str + Output file with MALI variable/dimension names + + start_year : int + First year (inclusive) to retain + + end_year : int + Last year (inclusive) to retain + """ + ds, years = open_rename_and_trim(remapped_file, + {"mask": "calvingMask"}, + start_year, end_year) + + # Round the remapped mask to 0/1 and store as integers. Ice shelves + # collapse on January 1st, so the mask is applied at the start of + # each year. + ds["calvingMask"] = (ds["calvingMask"] >= 0.5).astype(int) + + ds["calvingMask"].attrs = { + "long_name": "ice shelf collapse mask (1 = collapse)", + } + + add_xtime_and_write(ds, years, output_file) diff --git a/compass/landice/tests/ismip7_forcing/fracture/remap_utils.py b/compass/landice/tests/ismip7_forcing/fracture/remap_utils.py new file mode 100644 index 0000000000..8da5de4ae4 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/fracture/remap_utils.py @@ -0,0 +1,149 @@ +""" +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): + """ + Open a remapped file, rename dimensions/variables to MALI conventions, + and restrict to the requested year range. + + Parameters + ---------- + remapped_file : str + Data remapped onto the MALI mesh + + rename_vars : dict + Mapping of source variable names to MALI variable names + + start_year : int + First year (inclusive) to retain + + end_year : int + Last year (inclusive) to retain + + Returns + ------- + ds : xarray.Dataset + The renamed and trimmed dataset + + years : numpy.ndarray + The integer years retained + """ + # The time coordinate has units="year" (integer years), which is not + # CF-compliant, so disable time decoding. + ds = xr.open_dataset(remapped_file, decode_times=False) + + # Capture integer years before the time coordinate is renamed + years = ds["time"].values.astype(int) + + rename_dims = {} + if "ncol" in ds.dims: + rename_dims["ncol"] = "nCells" + if "time" in ds.dims: + rename_dims["time"] = "Time" + if rename_dims: + ds = ds.rename(rename_dims) + + rename_vars = {src: dst for src, dst in rename_vars.items() if src in ds} + if rename_vars: + ds = ds.rename(rename_vars) + + # Restrict to the requested year range + keep = (years >= start_year) & (years <= end_year) + ds = ds.isel(Time=keep) + years = years[keep] + + return ds, years + + +def add_xtime_and_write(ds, years, output_file): + """ + Add an ``xtime`` variable (January 1st of each year), drop auxiliary + remapping variables, and write the dataset. + + Parameters + ---------- + ds : xarray.Dataset + The dataset to finalize (annual fields applied at the start of the + year) + + years : numpy.ndarray + The integer years, one per Time index + + output_file : str + Output file path + """ + xtime = [f"{int(yr):04d}-01-01_00:00:00".ljust(64) for yr in years] + ds["xtime"] = ("Time", xtime) + ds["xtime"] = ds.xtime.astype("S") + + vars_to_drop = [v for v in ["lat_vertices", "lon_vertices", "lat", + "lon", "area", "Time"] + if v in ds] + if vars_to_drop: + ds = ds.drop_vars(vars_to_drop) + + write_netcdf(ds, output_file) + ds.close() diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg index b996d2aaba..f02b79067a 100644 --- a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg @@ -69,3 +69,28 @@ method_remap = bilinear # Base path to observational climatology data # (directory containing tf/, so/, thetao/ subdirs) base_path_climatology = /path/to/ISMIP7/forcing/AIS/obs/zhou_annual_06_nov + +# config options for ismip7 fracture (Path C, ice shelf collapse) forcing +[ismip7_fracture] + +# Remapping method for each pathway. Set a method to None to skip processing +# that pathway's file. Options: bilinear, neareststod, conserve, None + +# Remapping method for the ice shelf collapse mask (Path C). neareststod +# preserves the 0/1 mask values. +method_remap_shelf_collapse = neareststod + +# Remapping method for the excess meltwater field (Path A), a flux. +method_remap_excess_melt = conserve + +# Remapping method for the supraglacial lake properties (Path B). +method_remap_lake_properties = bilinear + +# Version subdirectory of the fracture forcing data +version = v2 + +# Start year for processing +start_year = 1850 + +# End year for processing +end_year = 2014 diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg index 270f5d396c..ab282b2e06 100644 --- a/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg @@ -69,3 +69,28 @@ method_remap = bilinear # Base path to observational climatology data # (directory containing tf/, so/, thetao/ subdirs) base_path_climatology = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/forcing/AIS/obs/zhou_annual_06_nov + +# config options for ismip7 fracture (Path C, ice shelf collapse) forcing +[ismip7_fracture] + +# Remapping method for each pathway. Set a method to None to skip processing +# that pathway's file. Options: bilinear, neareststod, conserve, None + +# Remapping method for the ice shelf collapse mask (Path C). neareststod +# preserves the 0/1 mask values. +method_remap_shelf_collapse = neareststod + +# Remapping method for the excess meltwater field (Path A), a flux. +method_remap_excess_melt = conserve + +# Remapping method for the supraglacial lake properties (Path B). +method_remap_lake_properties = bilinear + +# Version subdirectory of the fracture forcing data +version = v2 + +# Start year for processing +start_year = 2015 + +# End year for processing +end_year = 2300 diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index 57514fd526..d2cd995145 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -410,6 +410,21 @@ ismip7_forcing ocean_thermal.process_thermal_forcing.ProcessThermalForcing.setup ocean_thermal.process_thermal_forcing.ProcessThermalForcing.run + 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 + fracture.process_excess_melt.ProcessExcessMelt.setup + fracture.process_excess_melt.ProcessExcessMelt.run + fracture.process_lake_properties.ProcessLakeProperties + fracture.process_lake_properties.ProcessLakeProperties.setup + fracture.process_lake_properties.ProcessLakeProperties.run + fracture.process_shelf_collapse.ProcessShelfCollapse + fracture.process_shelf_collapse.ProcessShelfCollapse.setup + fracture.process_shelf_collapse.ProcessShelfCollapse.run + isunnguata_sermia ~~~~~~~~~~~~~~~~~ diff --git a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst index 3191ce9a6b..b535bf081a 100644 --- a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst @@ -9,7 +9,7 @@ The ``ismip7_forcing`` test group the Ice Sheet Model Intercomparison for CMIP7 (ISMIP7) protocol from its native polar stereographic grid to the MALI unstructured mesh. The test group supports both AIS and GrIS via the ``ice_sheet`` config option. It includes -two test cases: ``atmosphere`` and ``ocean_thermal``. +three test cases: ``atmosphere``, ``ocean_thermal``, and ``fracture``. .. _dev_landice_ismip7_forcing_framework: @@ -121,3 +121,56 @@ For GrIS, the step: * Remaps 2D monthly thermal forcing * Produces ``ismip6_2dThermalForcing`` (dims: Time × nCells) + +.. _dev_landice_ismip7_forcing_fracture: + +fracture +~~~~~~~~ + +The :py:class:`compass.landice.tests.ismip7_forcing.fracture.Fracture` +test case processes the ISMIP7 surface-melt-driven ice shelf collapse +forcing (AIS only). It implements the three ISMIP7 pathways as independent +steps, each discovering its source file from the ``fracture/{version}/`` +subdirectory of ``base_path_ismip7``, building or reusing a mapping file, +remapping with ``ncremap``, and renaming the result to MALI conventions with +an accompanying ``xtime`` variable. Per-pathway remapping methods are set in +the ``[ismip7_fracture]`` config section. Setting a pathway's remapping-method +option to ``None`` causes that step to return early without processing its +file, which is useful when only some pathway source files are available. + +Steps: + +* :py:class:`~compass.landice.tests.ismip7_forcing.fracture.process_excess_melt.ProcessExcessMelt` + (Path A) — ``excess_melt`` → ``ismip7ExcessMelt``. The excess melt file + lacks ``x``/``y`` coordinate variables and its array is flipped along the + y axis relative to the other fracture files (it was produced with CDO). + The ``_prepare_source_grid()`` method borrows ``x``/``y`` from a sibling + fracture file, flips the data to match (raising if the flipped ``lat`` does + not match the sibling grid), and writes a reconstructed source file. The + field is then extrapolated (nearest neighbor, filling NaNs) and remapped + conservatively by default (it is a flux). +* :py:class:`~compass.landice.tests.ismip7_forcing.fracture.process_lake_properties.ProcessLakeProperties` + (Path B) — ``lake_depth`` → ``ismip7LakeDepth`` and + ``fraction_lake_area`` → ``ismip7LakeAreaFraction``. Both variables are + extrapolated and remapped in a single ``ncremap`` call (bilinear by + default). +* :py:class:`~compass.landice.tests.ismip7_forcing.fracture.process_shelf_collapse.ProcessShelfCollapse` + (Path C) — ``mask`` → ``calvingMask``. Remapped with ``neareststod`` by + default and rounded to 0/1 so the discrete collapse mask is preserved. + +The annual source fields use an integer ``year``/``time`` coordinate with +``units="year"`` (not CF-compliant), so each step opens the data with +``decode_times=False`` and constructs ``xtime`` at January 1st of each year. + +Shared remapping helpers used by the fracture steps live in +:py:mod:`compass.landice.tests.ismip7_forcing.fracture.remap_utils`: +``extrapolate_source`` (nearest-neighbor fill of NaNs on the source grid), +``open_rename_and_trim`` (open a remapped file, rename dimensions/variables +to MALI conventions, and restrict to the requested year range), and +``add_xtime_and_write`` (add the ``xtime`` variable, drop auxiliary remapping +variables, and write the output). + +The output variable names for Paths A and B (``ismip7ExcessMelt``, +``ismip7LakeDepth``, ``ismip7LakeAreaFraction``) are descriptive placeholders +and may need to be aligned with the MALI Registry once the corresponding +model input fields are defined. diff --git a/docs/users_guide/landice/test_groups/ismip7_forcing.rst b/docs/users_guide/landice/test_groups/ismip7_forcing.rst index e26d3f88e9..fdec1eee14 100644 --- a/docs/users_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/users_guide/landice/test_groups/ismip7_forcing.rst @@ -10,7 +10,8 @@ force MALI in its simulations under the ISMIP7 experimental protocol. The test group supports both the Antarctic Ice Sheet (AIS) and the Greenland Ice Sheet (GrIS), controlled by a single ``ice_sheet`` config option. -The test group includes two test cases: ``atmosphere`` and ``ocean_thermal``. +The test group includes three test cases: ``atmosphere``, ``ocean_thermal``, +and ``fracture``. * The ``atmosphere`` test case has five steps: ``process_smb``, ``process_temperature``, ``process_smb_gradient``, @@ -22,9 +23,15 @@ The test group includes two test cases: ``atmosphere`` and ``ocean_thermal``. process the observational ocean thermal forcing climatology (Zhou et al.) for AIS, controlled by the ``process_ocean_climatology`` config option. +* The ``fracture`` test case has three steps: ``process_excess_melt`` + (Path A), ``process_lake_properties`` (Path B), and + ``process_shelf_collapse`` (Path C). It processes the ISMIP7 + surface-melt-driven ice shelf collapse forcing (AIS only). + (For more details on the steps of each test case, see -:ref:`landice_ismip7_forcing_atmosphere` and -:ref:`landice_ismip7_forcing_ocean_thermal`.) +:ref:`landice_ismip7_forcing_atmosphere`, +:ref:`landice_ismip7_forcing_ocean_thermal`, and +:ref:`landice_ismip7_forcing_fracture`.) .. _landice_ismip7_forcing_usage: @@ -43,6 +50,10 @@ To use this test group, users need to: 5. Run the ``ocean_thermal`` test case for each model and scenario combination. +6. Run the ``fracture`` test case (AIS only) for each model and scenario + combination to process the surface-melt-driven ice shelf collapse + pathways (excess melt, lake properties, and the ice shelf collapse mask). + Example user config files are provided in the source tree for local testing: * ``compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg`` @@ -88,6 +99,12 @@ For AIS ocean thermal climatology (8km, 30 depth levels, static): {base_path_climatology}/tf/v3/tf_AIS_obs_ocean_climatology_*.nc +For AIS fracture / ice shelf collapse mask (8km, annual, Path C): + +.. code-block:: none + + fracture/v2/ice_shelf_collapse_mask_*.nc + For GrIS atmosphere (1km, polar stereographic EPSG:3413): .. code-block:: none @@ -183,6 +200,31 @@ values are: # Base path to observational climatology data base_path_climatology = /path/to/ISMIP7/forcing/AIS/obs/zhou_annual_06_nov + # config options for ismip7 fracture (Path C, ice shelf collapse) forcing + [ismip7_fracture] + + # Remapping method for each pathway. Set a method to None to skip + # processing that pathway's file. + + # Remapping method for the ice shelf collapse mask (Path C). + # neareststod preserves the 0/1 mask values + method_remap_shelf_collapse = neareststod + + # Remapping method for the excess meltwater field (Path A), a flux + method_remap_excess_melt = conserve + + # Remapping method for the supraglacial lake properties (Path B) + method_remap_lake_properties = bilinear + + # Version subdirectory of the fracture forcing data + version = v2 + + # Start year for processing + start_year = 1850 + + # End year for processing + end_year = 2014 + All ``NotAvailable`` options must be overridden in a user config file passed at setup time (e.g., ``compass setup ... -f my_ismip7.cfg``). @@ -253,3 +295,58 @@ but without a Time dimension, producing a single static file. For **GrIS**, thermal forcing is 2D (depth-averaged), with monthly temporal resolution and yearly input files. The output variable is ``ismip6_2dThermalForcing``. + +.. _landice_ismip7_forcing_fracture: + +fracture +-------- + +The ``landice/ismip7_forcing/fracture`` test case processes the ISMIP7 +surface-melt-driven ice shelf collapse forcing (AIS only). It implements the +three ISMIP7 pathways as separate steps, each remapping annual fields from +the native 8km polar stereographic grid onto the MALI unstructured mesh. + +All three source files are discovered from the ``fracture/{version}/`` +subdirectory of ``base_path_ismip7``. + +Each pathway is run independently and can be skipped by setting its +remapping-method config option to ``None`` in the ``[ismip7_fracture]`` +section (for example, ``method_remap_excess_melt = None`` skips Path A). This +is useful when only some of the pathway source files are available. + +* **process_excess_melt** (Path A): Remaps the excess meltwater field + (melt + rain after firn air content depletion), matching + ``excess_melt_*.nc``. The output variable is ``ismip7ExcessMelt`` + (converted from mm w.e. yr-1 to SI units of kg m-2 s-1) and is written to + ``{output_base_path}/excess_melt/{model}_{scenario}/``. Conservative + remapping is used by default since this is a flux. This source file has no + ``x``/``y`` coordinate variables and its array is flipped along the y axis + relative to the other fracture files, so the step reconstructs the source + grid (borrowing ``x``/``y`` from a sibling fracture file and flipping the + data to match) before remapping. + +* **process_lake_properties** (Path B): Remaps the supraglacial lake mean + depth and area fraction from the Grau et al. (2025) parameterization, + matching ``lake_properties_*.nc``. The output variables are + ``ismip7LakeDepth`` (m) and ``ismip7LakeAreaFraction`` (unitless), written + to ``{output_base_path}/lake_properties/{model}_{scenario}/``. Bilinear + remapping is used by default. + +* **process_shelf_collapse** (Path C): Remaps the annual ice shelf collapse + mask, matching ``ice_shelf_collapse_mask_*.nc``. + An ice shelf grid cell is flagged as collapsed (mask value 1) when excess + meltwater, computed after firn air content depletion, exceeds 72.5 mm/yr for + 10 consecutive years; otherwise the mask value is 0. The mask is applied on + floating areas only, similar to ISMIP6, and ice shelves collapse on + January 1st of each year. The remapping uses ``neareststod`` by default so + that the 0/1 mask values are preserved, and the remapped mask is rounded to + 0/1. The output variable is ``calvingMask`` with an accompanying ``xtime`` + variable, and the result is written to + ``{output_base_path}/shelf_collapse/{model}_{scenario}/``. + +All three pathways produce continuous fields (Paths A and B) or a discrete +mask (Path C) with an accompanying ``xtime`` variable. The output variable +names for Paths A and B (``ismip7ExcessMelt``, ``ismip7LakeDepth``, +``ismip7LakeAreaFraction``) are descriptive placeholders and may need to be +aligned with the MALI Registry once the corresponding model input fields are +defined.