From a81e95a38c056d38e77c5009918809111b0fb27a Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 3 Jun 2026 16:26:34 -0700 Subject: [PATCH 01/15] Add ismip7_forcing test group with atmosphere test case Add ismip7_forcing test group and processing for atmospheric forcing (smb, lapse rate, and temperature) --- compass/landice/__init__.py | 2 + .../landice/tests/ismip7_forcing/__init__.py | 22 ++ .../ismip7_forcing/atmosphere/__init__.py | 44 ++++ .../ismip7_forcing/atmosphere/process_smb.py | 201 +++++++++++++++++ .../atmosphere/process_smb_gradient.py | 202 ++++++++++++++++++ .../atmosphere/process_temperature.py | 201 +++++++++++++++++ .../landice/tests/ismip7_forcing/configure.py | 21 ++ .../tests/ismip7_forcing/create_mapfile.py | 104 +++++++++ .../tests/ismip7_forcing/ismip7_forcing.cfg | 37 ++++ 9 files changed, 834 insertions(+) create mode 100644 compass/landice/tests/ismip7_forcing/__init__.py create mode 100644 compass/landice/tests/ismip7_forcing/atmosphere/__init__.py create mode 100644 compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py create mode 100644 compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py create mode 100644 compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py create mode 100644 compass/landice/tests/ismip7_forcing/configure.py create mode 100644 compass/landice/tests/ismip7_forcing/create_mapfile.py create mode 100644 compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg diff --git a/compass/landice/__init__.py b/compass/landice/__init__.py index 487f1e8345..6fec487fd7 100644 --- a/compass/landice/__init__.py +++ b/compass/landice/__init__.py @@ -11,6 +11,7 @@ from compass.landice.tests.hydro_radial import HydroRadial from compass.landice.tests.ismip6_forcing import Ismip6Forcing from compass.landice.tests.ismip6_run import Ismip6Run +from compass.landice.tests.ismip7_forcing import Ismip7Forcing from compass.landice.tests.isunnguata_sermia import IsunnguataSermia from compass.landice.tests.kangerlussuaq import Kangerlussuaq from compass.landice.tests.koge_bugt_s import KogeBugtS @@ -46,6 +47,7 @@ def __init__(self): self.add_test_group(HydroRadial(mpas_core=self)) self.add_test_group(Ismip6Forcing(mpas_core=self)) self.add_test_group(Ismip6Run(mpas_core=self)) + self.add_test_group(Ismip7Forcing(mpas_core=self)) self.add_test_group(IsunnguataSermia(mpas_core=self)) self.add_test_group(Kangerlussuaq(mpas_core=self)) self.add_test_group(KogeBugtS(mpas_core=self)) diff --git a/compass/landice/tests/ismip7_forcing/__init__.py b/compass/landice/tests/ismip7_forcing/__init__.py new file mode 100644 index 0000000000..689f37030c --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/__init__.py @@ -0,0 +1,22 @@ +from compass.landice.tests.ismip7_forcing.atmosphere import Atmosphere +from compass.testgroup import TestGroup + + +class Ismip7Forcing(TestGroup): + """ + A test group for processing ISMIP7 atmosphere forcing data + for the Antarctic Ice Sheet + """ + + def __init__(self, mpas_core): + """ + Create the test group + + Parameters + ---------- + mpas_core : compass.landice.Landice + the MPAS core that this test group belongs to + """ + super().__init__(mpas_core=mpas_core, name="ismip7_forcing") + + self.add_test_case(Atmosphere(test_group=self)) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py b/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py new file mode 100644 index 0000000000..be6b90de77 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py @@ -0,0 +1,44 @@ +from compass.landice.tests.ismip7_forcing.atmosphere.process_smb import ( + ProcessSmb, +) +from compass.landice.tests.ismip7_forcing.atmosphere.process_smb_gradient import ( # noqa: E501 + ProcessSmbGradient, +) +from compass.landice.tests.ismip7_forcing.atmosphere.process_temperature import ( # noqa: E501 + ProcessTemperature, +) +from compass.landice.tests.ismip7_forcing.configure import ( + configure as configure_testgroup, +) +from compass.testcase import TestCase + + +class Atmosphere(TestCase): + """ + A test case for processing ISMIP7 AIS atmosphere forcing data. + Remaps monthly SMB, temperature, and annual SMB gradient from the + ISMIP7 2km polar stereographic grid to the MALI unstructured mesh. + """ + + 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 = "atmosphere" + subdir = name + super().__init__(test_group=test_group, name=name, subdir=subdir) + + self.add_step(ProcessSmb(test_case=self)) + self.add_step(ProcessTemperature(test_case=self)) + self.add_step(ProcessSmbGradient(test_case=self)) + + def configure(self): + """ + Configures test case + """ + configure_testgroup(config=self.config) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py new file mode 100644 index 0000000000..3b69102b51 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -0,0 +1,201 @@ +import glob +import os +import shutil + +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.step import Step + + +class ProcessSmb(Step): + """ + A step for processing ISMIP7 surface mass balance (acabf) data. + Remaps monthly full-field SMB from the ISMIP7 2km polar stereographic + grid to the MALI unstructured mesh. + """ + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_forcing.atmosphere.Atmosphere + The test case this step belongs to + """ + super().__init__(test_case=test_case, name="process_smb", + ntasks=4, min_tasks=1) + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config["ismip7_ais"] + 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_ais"] + 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") + + section = config["ismip7_ais_atmosphere"] + method_remap = section.get("method_remap") + start_year = section.getint("start_year") + end_year = section.getint("end_year") + + # Discover input files + input_path = os.path.join(base_path_ismip7, "acabf", "v2") + file_pattern = f"acabf_AIS_{model}_{scenario}_SDBN1-2000m_v2_*.nc" + all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) + + if not all_files: + raise FileNotFoundError( + f"No SMB files found matching pattern:\n" + f" {os.path.join(input_path, file_pattern)}") + + # 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", "")) + if start_year <= year <= end_year: + input_files.append(f) + + if not input_files: + raise FileNotFoundError( + f"No SMB files found for year range {start_year}-{end_year}") + + logger.info(f"Found {len(input_files)} SMB files for years " + f"{start_year}-{end_year}") + + # Build mapping file using the first input file as the grid template + mapping_file = f"map_ismip7_2km_to_{mali_mesh_name}_{method_remap}.nc" + + if not os.path.exists(mapping_file): + logger.info("Building mapping file...") + build_mapping_file(config, self.ntasks, logger, + input_files[0], mapping_file, + mali_mesh_file=mali_mesh_file, + method_remap=method_remap) + + # Remap each year file + remapped_files = [] + for input_file in input_files: + basename = os.path.basename(input_file) + remapped_file = f"remapped_{basename}" + remapped_files.append(remapped_file) + + if os.path.exists(remapped_file): + logger.info(f" Remapped file exists, skipping: {basename}") + continue + + logger.info(f" Remapping: {basename}") + args = ["ncremap", + "-i", input_file, + "-o", remapped_file, + "-m", mapping_file, + "-v", "acabf"] + + check_call(args, logger=logger) + + # 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}_" + f"{start_year}-{end_year}.nc") + + self._combine_and_rename(remapped_files, output_file) + + # Clean up remapped files + logger.info("Cleaning up temporary remapped files...") + for f in remapped_files: + if os.path.exists(f): + os.remove(f) + + # Place output in appropriate directory + output_path = os.path.join(output_base_path, "atmosphere_forcing", + 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 _combine_and_rename(self, remapped_files, output_file): + """ + Combine yearly remapped files and rename variables/dimensions + to MALI conventions. + + Parameters + ---------- + remapped_files : list of str + List of remapped NetCDF file paths + + output_file : str + Output file path + """ + ds = xr.open_mfdataset(remapped_files, concat_dim="time", + combine="nested", engine="netcdf4") + + # Rename dimensions to MALI conventions + rename_dims = {} + if "time" in ds.dims: + rename_dims["time"] = "Time" + if "ncol" in ds.dims: + rename_dims["ncol"] = "nCells" + if rename_dims: + ds = ds.rename(rename_dims) + + # Rename variable + if "acabf" in ds: + ds = ds.rename({"acabf": "sfcMassBal"}) + + # Add xtime variable with monthly timestamps + xtime = [] + for t_index in range(ds.sizes["Time"]): + date = ds.Time[t_index] + yr = int(date.dt.year.values) + mo = int(date.dt.month.values) + date_str = f"{yr:04d}-{mo:02d}-15_00:00:00".ljust(64) + xtime.append(date_str) + + ds["xtime"] = ("Time", xtime) + ds["xtime"] = ds.xtime.astype("S") + + # Set attributes + ds["sfcMassBal"].attrs = { + "long_name": "surface mass balance", + "units": "kg m-2 s-1", + } + + # Drop auxiliary variables from remapping + vars_to_drop = [v for v in ["lon", "lon_vertices", "lat", + "lat_vertices", "area"] + 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/atmosphere/process_smb_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py new file mode 100644 index 0000000000..d87f84d2e6 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -0,0 +1,202 @@ +import glob +import os +import shutil + +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.step import Step + + +class ProcessSmbGradient(Step): + """ + A step for processing ISMIP7 SMB elevation gradient (dacabfdz) data. + Remaps the annual SMB gradient from the ISMIP7 2km polar stereographic + grid to the MALI unstructured mesh. This field is used for + ice-elevation feedback corrections. + """ + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_forcing.atmosphere.Atmosphere + The test case this step belongs to + """ + super().__init__(test_case=test_case, name="process_smb_gradient", + ntasks=4, min_tasks=1) + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config["ismip7_ais"] + 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_ais"] + 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") + + section = config["ismip7_ais_atmosphere"] + method_remap = section.get("method_remap") + start_year = section.getint("start_year") + end_year = section.getint("end_year") + + # Discover input files + input_path = os.path.join(base_path_ismip7, "dacabfdz", "v2") + file_pattern = (f"dacabfdz_AIS_{model}_{scenario}_" + f"SDBN1-2000m_v2_*.nc") + all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) + + if not all_files: + raise FileNotFoundError( + f"No SMB gradient files found matching pattern:\n" + f" {os.path.join(input_path, file_pattern)}") + + # Filter to requested year range + input_files = [] + for f in all_files: + year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + if start_year <= year <= end_year: + input_files.append(f) + + if not input_files: + raise FileNotFoundError( + f"No SMB gradient files for year range " + f"{start_year}-{end_year}") + + logger.info(f"Found {len(input_files)} SMB gradient files for years " + f"{start_year}-{end_year}") + + # Build mapping file (reuse if already created by process_smb) + mapping_file = f"map_ismip7_2km_to_{mali_mesh_name}_{method_remap}.nc" + + if not os.path.exists(mapping_file): + logger.info("Building mapping file...") + build_mapping_file(config, self.ntasks, logger, + input_files[0], mapping_file, + mali_mesh_file=mali_mesh_file, + method_remap=method_remap) + + # Remap each year file + remapped_files = [] + for input_file in input_files: + basename = os.path.basename(input_file) + remapped_file = f"remapped_{basename}" + remapped_files.append(remapped_file) + + if os.path.exists(remapped_file): + logger.info(f" Remapped file exists, skipping: {basename}") + continue + + logger.info(f" Remapping: {basename}") + args = ["ncremap", + "-i", input_file, + "-o", remapped_file, + "-m", mapping_file, + "-v", "dacabfdz"] + + check_call(args, logger=logger) + + # 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}_" + f"{start_year}-{end_year}.nc") + + self._combine_and_rename(remapped_files, output_file) + + # Clean up remapped files + logger.info("Cleaning up temporary remapped files...") + for f in remapped_files: + if os.path.exists(f): + os.remove(f) + + # Place output in appropriate directory + output_path = os.path.join(output_base_path, "atmosphere_forcing", + 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 _combine_and_rename(self, remapped_files, output_file): + """ + Combine yearly remapped files and rename variables/dimensions + to MALI conventions. + + Parameters + ---------- + remapped_files : list of str + List of remapped NetCDF file paths + + output_file : str + Output file path + """ + ds = xr.open_mfdataset(remapped_files, concat_dim="time", + combine="nested", engine="netcdf4") + + # Rename dimensions to MALI conventions + rename_dims = {} + if "time" in ds.dims: + rename_dims["time"] = "Time" + if "ncol" in ds.dims: + rename_dims["ncol"] = "nCells" + if rename_dims: + ds = ds.rename(rename_dims) + + # Keep variable name as dacabfdz for now. + # The MALI variable name will be determined by PR #169. + # Rename can be updated once that PR is merged. + + # Add xtime variable with annual timestamps + xtime = [] + for t_index in range(ds.sizes["Time"]): + date = ds.Time[t_index] + yr = int(date.dt.year.values) + date_str = f"{yr:04d}-01-01_00:00:00".ljust(64) + xtime.append(date_str) + + ds["xtime"] = ("Time", xtime) + ds["xtime"] = ds.xtime.astype("S") + + # Set attributes + ds["dacabfdz"].attrs = { + "long_name": "surface mass balance change with surface elevation", + "units": "kg m-2 s-1 m-1", + } + + # Drop auxiliary variables from remapping + vars_to_drop = [v for v in ["lon", "lon_vertices", "lat", + "lat_vertices", "area"] + 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/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py new file mode 100644 index 0000000000..61c06693c0 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -0,0 +1,201 @@ +import glob +import os +import shutil + +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.step import Step + + +class ProcessTemperature(Step): + """ + A step for processing ISMIP7 near-surface air temperature (tas) data. + Remaps monthly temperature from the ISMIP7 2km polar stereographic + grid to the MALI unstructured mesh. + """ + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_forcing.atmosphere.Atmosphere + The test case this step belongs to + """ + super().__init__(test_case=test_case, name="process_temperature", + ntasks=4, min_tasks=1) + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config["ismip7_ais"] + 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_ais"] + 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") + + section = config["ismip7_ais_atmosphere"] + method_remap = section.get("method_remap") + start_year = section.getint("start_year") + end_year = section.getint("end_year") + + # Discover input files + input_path = os.path.join(base_path_ismip7, "tas", "v2") + file_pattern = f"tas_AIS_{model}_{scenario}_SDBN1-2000m_v2_*.nc" + all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) + + if not all_files: + raise FileNotFoundError( + f"No temperature files found matching pattern:\n" + f" {os.path.join(input_path, file_pattern)}") + + # Filter to requested year range + input_files = [] + for f in all_files: + year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + if start_year <= year <= end_year: + input_files.append(f) + + if not input_files: + raise FileNotFoundError( + f"No temperature files for year range " + f"{start_year}-{end_year}") + + logger.info(f"Found {len(input_files)} temperature files for years " + f"{start_year}-{end_year}") + + # Build mapping file (reuse if already created by process_smb) + mapping_file = f"map_ismip7_2km_to_{mali_mesh_name}_{method_remap}.nc" + + if not os.path.exists(mapping_file): + logger.info("Building mapping file...") + build_mapping_file(config, self.ntasks, logger, + input_files[0], mapping_file, + mali_mesh_file=mali_mesh_file, + method_remap=method_remap) + + # Remap each year file + remapped_files = [] + for input_file in input_files: + basename = os.path.basename(input_file) + remapped_file = f"remapped_{basename}" + remapped_files.append(remapped_file) + + if os.path.exists(remapped_file): + logger.info(f" Remapped file exists, skipping: {basename}") + continue + + logger.info(f" Remapping: {basename}") + args = ["ncremap", + "-i", input_file, + "-o", remapped_file, + "-m", mapping_file, + "-v", "tas"] + + check_call(args, logger=logger) + + # 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}_" + f"{start_year}-{end_year}.nc") + + self._combine_and_rename(remapped_files, output_file) + + # Clean up remapped files + logger.info("Cleaning up temporary remapped files...") + for f in remapped_files: + if os.path.exists(f): + os.remove(f) + + # Place output in appropriate directory + output_path = os.path.join(output_base_path, "atmosphere_forcing", + 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 _combine_and_rename(self, remapped_files, output_file): + """ + Combine yearly remapped files and rename variables/dimensions + to MALI conventions. + + Parameters + ---------- + remapped_files : list of str + List of remapped NetCDF file paths + + output_file : str + Output file path + """ + ds = xr.open_mfdataset(remapped_files, concat_dim="time", + combine="nested", engine="netcdf4") + + # Rename dimensions to MALI conventions + rename_dims = {} + if "time" in ds.dims: + rename_dims["time"] = "Time" + if "ncol" in ds.dims: + rename_dims["ncol"] = "nCells" + if rename_dims: + ds = ds.rename(rename_dims) + + # Rename variable + if "tas" in ds: + ds = ds.rename({"tas": "surfaceAirTemperature"}) + + # Add xtime variable with monthly timestamps + xtime = [] + for t_index in range(ds.sizes["Time"]): + date = ds.Time[t_index] + yr = int(date.dt.year.values) + mo = int(date.dt.month.values) + date_str = f"{yr:04d}-{mo:02d}-15_00:00:00".ljust(64) + xtime.append(date_str) + + ds["xtime"] = ("Time", xtime) + ds["xtime"] = ds.xtime.astype("S") + + # Set attributes + ds["surfaceAirTemperature"].attrs = { + "long_name": "near-surface air temperature", + "units": "K", + } + + # Drop auxiliary variables from remapping + vars_to_drop = [v for v in ["lon", "lon_vertices", "lat", + "lat_vertices", "area"] + 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/configure.py b/compass/landice/tests/ismip7_forcing/configure.py new file mode 100644 index 0000000000..c49370206c --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/configure.py @@ -0,0 +1,21 @@ +def configure(config): + """ + A shared function for configuring options for all ismip7 forcing + test cases + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for an ismip7 forcing test case + """ + + section = "ismip7_ais" + options = ["base_path_ismip7", "base_path_mali", "mali_mesh_name", + "mali_mesh_file", "output_base_path", "model", "scenario"] + + for option in options: + value = config.get(section=section, option=option) + if value == "NotAvailable": + raise ValueError(f"You need to supply a user config file, which " + f"should contain the {section} " + f"section with the {option} option") diff --git a/compass/landice/tests/ismip7_forcing/create_mapfile.py b/compass/landice/tests/ismip7_forcing/create_mapfile.py new file mode 100644 index 0000000000..31732443ed --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/create_mapfile.py @@ -0,0 +1,104 @@ +import os +import shutil + +from mpas_tools.logging import check_call +from mpas_tools.scrip.from_mpas import scrip_from_mpas + + +def build_mapping_file(config, cores, logger, ismip7_grid_file, + mapping_file, mali_mesh_file=None, + method_remap=None): + """ + Build a mapping file for regridding from the ISMIP7 2km polar + stereographic grid to the MALI unstructured mesh. + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for the test case + + cores : int + the number of cores for ESMF_RegridWeightGen + + logger : logging.Logger + A logger for output from the step + + ismip7_grid_file : str + An ISMIP7 grid file (with x/y coordinates) + + mapping_file : str + Output mapping file path + + mali_mesh_file : str, optional + The MALI mesh file + + method_remap : str, optional + Remapping method: 'bilinear', 'neareststod', or 'conserve' + """ + + if os.path.exists(mapping_file): + logger.info("Mapping file exists. Not building a new one.") + return + + logger.info("Mapping file does not exist. Building one based on the" + " input/output meshes") + + if mali_mesh_file is None: + raise ValueError("Mapping file does not exist. A MALI mesh file " + "must be provided to build one.") + + if method_remap is None: + raise ValueError("Remapping method must be provided. " + "Options: 'bilinear', 'neareststod', 'conserve'.") + + # AIS polar stereographic projection (EPSG:3031) + ismip7_projection = "ais-bedmap2" + + # name temporary scrip files + source_grid_scripfile = "temp_source_scrip.nc" + mali_scripfile = "temp_mali_scrip.nc" + + # create the scrip file for the ISMIP7 planar rectangular grid + logger.info("Creating SCRIP file for ISMIP7 source grid...") + args = ["create_scrip_file_from_planar_rectangular_grid", + "--input", ismip7_grid_file, + "--scrip", source_grid_scripfile, + "--proj", ismip7_projection, + "--rank", "2"] + + check_call(args, logger=logger) + + # create a MALI mesh scrip file + logger.info("Creating SCRIP file for MALI mesh...") + mali_mesh_copy = f"{mali_mesh_file}_copy" + shutil.copy(mali_mesh_file, mali_mesh_copy) + + args = ["set_lat_lon_fields_in_planar_grid", + "--file", mali_mesh_copy, + "--proj", ismip7_projection] + + check_call(args, logger=logger) + + scrip_from_mpas(mali_mesh_copy, mali_scripfile) + + # create a mapping file using ESMF_RegridWeightGen + logger.info(f"Creating mapping file with method: {method_remap}") + + parallel_executable = config.get("parallel", "parallel_executable") + args = parallel_executable.split(" ") + args.extend(["-n", f"{cores}", + "ESMF_RegridWeightGen", + "-s", source_grid_scripfile, + "-d", mali_scripfile, + "-w", mapping_file, + "-m", method_remap, + "-i", "-64bit_offset", + "--dst_regional", "--src_regional"]) + + check_call(args, logger=logger) + + # clean up temporary files + logger.info("Removing temporary scrip files...") + os.remove(source_grid_scripfile) + os.remove(mali_scripfile) + os.remove(mali_mesh_copy) diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg new file mode 100644 index 0000000000..d67c395b71 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg @@ -0,0 +1,37 @@ +# config options for ismip7 antarctic ice sheet forcing data +[ismip7_ais] + +# Base path to the input ISMIP7 forcing files. User has to supply. +base_path_ismip7 = NotAvailable + +# Base path to the MALI mesh. User has to supply. +base_path_mali = NotAvailable + +# Base path to which output forcing files are saved. +output_base_path = NotAvailable + +# Name of climate model used to generate ISMIP7 forcing data. User has to supply. +# Available model names: CESM2-WACCM, MRI-ESM2-0 +model = NotAvailable + +# Scenario for forcing data. User has to supply. +# Available scenarios: historical, ssp126, ssp370, ssp585 +scenario = NotAvailable + +# Name of the MALI mesh. Used to name mapping and output files. +mali_mesh_name = NotAvailable + +# MALI mesh file (e.g. Antarctic_8to80km_20220407.nc). User has to supply. +mali_mesh_file = NotAvailable + +# config options for ismip7 AIS atmosphere forcing +[ismip7_ais_atmosphere] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = conserve + +# Start year for processing +start_year = 1850 + +# End year for processing +end_year = 2014 From 0c03c5b6d61ea4c2020e314d7de6883e5b8dc9b2 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 3 Jun 2026 16:41:57 -0700 Subject: [PATCH 02/15] Add temperature lapse rate and use ts instead of tas --- .../ismip7_forcing/atmosphere/__init__.py | 9 +- .../atmosphere/process_temperature.py | 14 +- .../process_temperature_gradient.py | 202 ++++++++++++++++++ 3 files changed, 216 insertions(+), 9 deletions(-) create mode 100644 compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py b/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py index be6b90de77..472e3f6d9b 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py @@ -7,6 +7,9 @@ from compass.landice.tests.ismip7_forcing.atmosphere.process_temperature import ( # noqa: E501 ProcessTemperature, ) +from compass.landice.tests.ismip7_forcing.atmosphere.process_temperature_gradient import ( # noqa: E501 + ProcessTemperatureGradient, +) from compass.landice.tests.ismip7_forcing.configure import ( configure as configure_testgroup, ) @@ -16,8 +19,9 @@ class Atmosphere(TestCase): """ A test case for processing ISMIP7 AIS atmosphere forcing data. - Remaps monthly SMB, temperature, and annual SMB gradient from the - ISMIP7 2km polar stereographic grid to the MALI unstructured mesh. + Remaps monthly SMB, temperature, and annual gradients (SMB and + temperature) from the ISMIP7 2km polar stereographic grid to the + MALI unstructured mesh. """ def __init__(self, test_group): @@ -36,6 +40,7 @@ def __init__(self, test_group): self.add_step(ProcessSmb(test_case=self)) self.add_step(ProcessTemperature(test_case=self)) self.add_step(ProcessSmbGradient(test_case=self)) + self.add_step(ProcessTemperatureGradient(test_case=self)) def configure(self): """ diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index 61c06693c0..6fed78ef56 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -14,7 +14,7 @@ class ProcessTemperature(Step): """ - A step for processing ISMIP7 near-surface air temperature (tas) data. + A step for processing ISMIP7 ice surface temperature (ts) data. Remaps monthly temperature from the ISMIP7 2km polar stereographic grid to the MALI unstructured mesh. """ @@ -65,8 +65,8 @@ def run(self): end_year = section.getint("end_year") # Discover input files - input_path = os.path.join(base_path_ismip7, "tas", "v2") - file_pattern = f"tas_AIS_{model}_{scenario}_SDBN1-2000m_v2_*.nc" + input_path = os.path.join(base_path_ismip7, "ts", "v2") + file_pattern = f"ts_AIS_{model}_{scenario}_SDBN1-2000m_v2_*.nc" all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) if not all_files: @@ -115,7 +115,7 @@ def run(self): "-i", input_file, "-o", remapped_file, "-m", mapping_file, - "-v", "tas"] + "-v", "ts"] check_call(args, logger=logger) @@ -169,8 +169,8 @@ def _combine_and_rename(self, remapped_files, output_file): ds = ds.rename(rename_dims) # Rename variable - if "tas" in ds: - ds = ds.rename({"tas": "surfaceAirTemperature"}) + if "ts" in ds: + ds = ds.rename({"ts": "surfaceAirTemperature"}) # Add xtime variable with monthly timestamps xtime = [] @@ -186,7 +186,7 @@ def _combine_and_rename(self, remapped_files, output_file): # Set attributes ds["surfaceAirTemperature"].attrs = { - "long_name": "near-surface air temperature", + "long_name": "temperature at top of ice sheet model", "units": "K", } diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py new file mode 100644 index 0000000000..76c2fb14ac --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -0,0 +1,202 @@ +import glob +import os +import shutil + +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.step import Step + + +class ProcessTemperatureGradient(Step): + """ + A step for processing ISMIP7 temperature elevation gradient (dtsdz) data. + Remaps the annual temperature gradient from the ISMIP7 2km polar + stereographic grid to the MALI unstructured mesh. This field is used + for temperature-elevation feedback corrections. + """ + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_forcing.atmosphere.Atmosphere + The test case this step belongs to + """ + super().__init__(test_case=test_case, + name="process_temperature_gradient", + ntasks=4, min_tasks=1) + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config["ismip7_ais"] + 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_ais"] + 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") + + section = config["ismip7_ais_atmosphere"] + method_remap = section.get("method_remap") + start_year = section.getint("start_year") + end_year = section.getint("end_year") + + # Discover input files + input_path = os.path.join(base_path_ismip7, "dtsdz", "v2") + file_pattern = (f"dtsdz_AIS_{model}_{scenario}_" + f"SDBN1-2000m_v2_*.nc") + all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) + + if not all_files: + raise FileNotFoundError( + f"No temperature gradient files found matching pattern:\n" + f" {os.path.join(input_path, file_pattern)}") + + # Filter to requested year range + input_files = [] + for f in all_files: + year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + if start_year <= year <= end_year: + input_files.append(f) + + if not input_files: + raise FileNotFoundError( + f"No temperature gradient files for year range " + f"{start_year}-{end_year}") + + logger.info(f"Found {len(input_files)} temperature gradient files " + f"for years {start_year}-{end_year}") + + # Build mapping file (reuse if already created by other steps) + mapping_file = f"map_ismip7_2km_to_{mali_mesh_name}_{method_remap}.nc" + + if not os.path.exists(mapping_file): + logger.info("Building mapping file...") + build_mapping_file(config, self.ntasks, logger, + input_files[0], mapping_file, + mali_mesh_file=mali_mesh_file, + method_remap=method_remap) + + # Remap each year file + remapped_files = [] + for input_file in input_files: + basename = os.path.basename(input_file) + remapped_file = f"remapped_{basename}" + remapped_files.append(remapped_file) + + if os.path.exists(remapped_file): + logger.info(f" Remapped file exists, skipping: {basename}") + continue + + logger.info(f" Remapping: {basename}") + args = ["ncremap", + "-i", input_file, + "-o", remapped_file, + "-m", mapping_file, + "-v", "dtsdz"] + + check_call(args, logger=logger) + + # 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}_" + f"{scenario}_{start_year}-{end_year}.nc") + + self._combine_and_rename(remapped_files, output_file) + + # Clean up remapped files + logger.info("Cleaning up temporary remapped files...") + for f in remapped_files: + if os.path.exists(f): + os.remove(f) + + # Place output in appropriate directory + output_path = os.path.join(output_base_path, "atmosphere_forcing", + 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 _combine_and_rename(self, remapped_files, output_file): + """ + Combine yearly remapped files and rename variables/dimensions + to MALI conventions. + + Parameters + ---------- + remapped_files : list of str + List of remapped NetCDF file paths + + output_file : str + Output file path + """ + ds = xr.open_mfdataset(remapped_files, concat_dim="time", + combine="nested", engine="netcdf4") + + # Rename dimensions to MALI conventions + rename_dims = {} + if "time" in ds.dims: + rename_dims["time"] = "Time" + if "ncol" in ds.dims: + rename_dims["ncol"] = "nCells" + if rename_dims: + ds = ds.rename(rename_dims) + + # Keep variable name as dtsdz for now. + # The MALI variable name will be determined once temperature + # elevation feedback support is added to MALI. + + # Add xtime variable with annual timestamps + xtime = [] + for t_index in range(ds.sizes["Time"]): + date = ds.Time[t_index] + yr = int(date.dt.year.values) + date_str = f"{yr:04d}-01-01_00:00:00".ljust(64) + xtime.append(date_str) + + ds["xtime"] = ("Time", xtime) + ds["xtime"] = ds.xtime.astype("S") + + # Set attributes + ds["dtsdz"].attrs = { + "long_name": "temperature change with surface elevation", + "units": "K m-1", + } + + # Drop auxiliary variables from remapping + vars_to_drop = [v for v in ["lon", "lon_vertices", "lat", + "lat_vertices", "area"] + if v in ds] + if vars_to_drop: + ds = ds.drop_vars(vars_to_drop) + + write_netcdf(ds, output_file) From 9c100882a46cfdb4d1c8ec6c371c7761146cf671 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 3 Jun 2026 19:27:42 -0700 Subject: [PATCH 03/15] Allow user to control number of tasks for ESMF_RegridWeightGen Allow user to control number of tasks for ESMF_RegridWeightGen, which easily throws a segmentation fault on too few processors. Use 512 cores when processing 2km source data. --- .../tests/ismip7_forcing/atmosphere/process_smb.py | 5 ++--- .../ismip7_forcing/atmosphere/process_smb_gradient.py | 5 ++--- .../ismip7_forcing/atmosphere/process_temperature.py | 5 ++--- .../atmosphere/process_temperature_gradient.py | 5 ++--- compass/landice/tests/ismip7_forcing/create_mapfile.py | 8 ++++---- compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg | 4 ++++ 6 files changed, 16 insertions(+), 16 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py index 3b69102b51..e619767e25 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -28,8 +28,7 @@ def __init__(self, test_case): test_case : compass.landice.tests.ismip7_forcing.atmosphere.Atmosphere The test case this step belongs to """ - super().__init__(test_case=test_case, name="process_smb", - ntasks=4, min_tasks=1) + super().__init__(test_case=test_case, name="process_smb") def setup(self): """ @@ -94,7 +93,7 @@ def run(self): if not os.path.exists(mapping_file): logger.info("Building mapping file...") - build_mapping_file(config, self.ntasks, logger, + build_mapping_file(config, logger, input_files[0], mapping_file, mali_mesh_file=mali_mesh_file, method_remap=method_remap) 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 d87f84d2e6..d9d91fc26e 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -29,8 +29,7 @@ def __init__(self, test_case): test_case : compass.landice.tests.ismip7_forcing.atmosphere.Atmosphere The test case this step belongs to """ - super().__init__(test_case=test_case, name="process_smb_gradient", - ntasks=4, min_tasks=1) + super().__init__(test_case=test_case, name="process_smb_gradient") def setup(self): """ @@ -96,7 +95,7 @@ def run(self): if not os.path.exists(mapping_file): logger.info("Building mapping file...") - build_mapping_file(config, self.ntasks, logger, + build_mapping_file(config, logger, input_files[0], mapping_file, mali_mesh_file=mali_mesh_file, method_remap=method_remap) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index 6fed78ef56..d5b6fe77ce 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -28,8 +28,7 @@ def __init__(self, test_case): test_case : compass.landice.tests.ismip7_forcing.atmosphere.Atmosphere The test case this step belongs to """ - super().__init__(test_case=test_case, name="process_temperature", - ntasks=4, min_tasks=1) + super().__init__(test_case=test_case, name="process_temperature") def setup(self): """ @@ -94,7 +93,7 @@ def run(self): if not os.path.exists(mapping_file): logger.info("Building mapping file...") - build_mapping_file(config, self.ntasks, logger, + build_mapping_file(config, logger, input_files[0], mapping_file, mali_mesh_file=mali_mesh_file, method_remap=method_remap) 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 76c2fb14ac..7829dc57be 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -30,8 +30,7 @@ def __init__(self, test_case): The test case this step belongs to """ super().__init__(test_case=test_case, - name="process_temperature_gradient", - ntasks=4, min_tasks=1) + name="process_temperature_gradient") def setup(self): """ @@ -97,7 +96,7 @@ def run(self): if not os.path.exists(mapping_file): logger.info("Building mapping file...") - build_mapping_file(config, self.ntasks, logger, + build_mapping_file(config, logger, input_files[0], mapping_file, mali_mesh_file=mali_mesh_file, method_remap=method_remap) diff --git a/compass/landice/tests/ismip7_forcing/create_mapfile.py b/compass/landice/tests/ismip7_forcing/create_mapfile.py index 31732443ed..81ea411ba6 100644 --- a/compass/landice/tests/ismip7_forcing/create_mapfile.py +++ b/compass/landice/tests/ismip7_forcing/create_mapfile.py @@ -5,7 +5,7 @@ from mpas_tools.scrip.from_mpas import scrip_from_mpas -def build_mapping_file(config, cores, logger, ismip7_grid_file, +def build_mapping_file(config, logger, ismip7_grid_file, mapping_file, mali_mesh_file=None, method_remap=None): """ @@ -17,9 +17,6 @@ def build_mapping_file(config, cores, logger, ismip7_grid_file, config : compass.config.CompassConfigParser Configuration options for the test case - cores : int - the number of cores for ESMF_RegridWeightGen - logger : logging.Logger A logger for output from the step @@ -84,6 +81,9 @@ def build_mapping_file(config, cores, logger, ismip7_grid_file, # create a mapping file using ESMF_RegridWeightGen logger.info(f"Creating mapping file with method: {method_remap}") + section = config["ismip7_ais"] + cores = section.getint("esmf_ntasks") + parallel_executable = config.get("parallel", "parallel_executable") args = parallel_executable.split(" ") args.extend(["-n", f"{cores}", diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg index d67c395b71..75fd603ffb 100644 --- a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg @@ -24,6 +24,10 @@ mali_mesh_name = NotAvailable # MALI mesh file (e.g. Antarctic_8to80km_20220407.nc). User has to supply. mali_mesh_file = NotAvailable +# Number of MPI tasks for ESMF_RegridWeightGen +# Use 512 for 2km data sets. +esmf_ntasks = 128 + # config options for ismip7 AIS atmosphere forcing [ismip7_ais_atmosphere] From 129494505d2d0944b880ecb2ea0c926594479c8f Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 3 Jun 2026 19:47:43 -0700 Subject: [PATCH 04/15] Update lapse rate variable names to be consistent with MALI fields --- .../ismip7_forcing/atmosphere/process_smb_gradient.py | 11 ++++++----- .../atmosphere/process_temperature_gradient.py | 11 ++++++----- 2 files changed, 12 insertions(+), 10 deletions(-) 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 d9d91fc26e..88fed545f1 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -169,9 +169,9 @@ def _combine_and_rename(self, remapped_files, output_file): if rename_dims: ds = ds.rename(rename_dims) - # Keep variable name as dacabfdz for now. - # The MALI variable name will be determined by PR #169. - # Rename can be updated once that PR is merged. + # Rename to MALI convention (PR #169) + if "dacabfdz" in ds: + ds = ds.rename({"dacabfdz": "sfcMassBalLapseRate"}) # Add xtime variable with annual timestamps xtime = [] @@ -185,8 +185,9 @@ def _combine_and_rename(self, remapped_files, output_file): ds["xtime"] = ds.xtime.astype("S") # Set attributes - ds["dacabfdz"].attrs = { - "long_name": "surface mass balance change with surface elevation", + ds["sfcMassBalLapseRate"].attrs = { + "long_name": "vertical gradient dSMBdz used for SMB " + "elevation-change correction", "units": "kg m-2 s-1 m-1", } 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 7829dc57be..9ed4fd4a99 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -170,9 +170,9 @@ def _combine_and_rename(self, remapped_files, output_file): if rename_dims: ds = ds.rename(rename_dims) - # Keep variable name as dtsdz for now. - # The MALI variable name will be determined once temperature - # elevation feedback support is added to MALI. + # Rename to MALI convention (PR #169) + if "dtsdz" in ds: + ds = ds.rename({"dtsdz": "surfaceAirTemperatureLapseRate"}) # Add xtime variable with annual timestamps xtime = [] @@ -186,8 +186,9 @@ def _combine_and_rename(self, remapped_files, output_file): ds["xtime"] = ds.xtime.astype("S") # Set attributes - ds["dtsdz"].attrs = { - "long_name": "temperature change with surface elevation", + ds["surfaceAirTemperatureLapseRate"].attrs = { + "long_name": "vertical gradient dTdz used for SAT " + "elevation-change correction", "units": "K m-1", } From 5fd1de523008318742ec5d62043f80906eedfa70 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 3 Jun 2026 20:56:25 -0700 Subject: [PATCH 05/15] Add ocean thermal forcing processing for ISMIP7 --- .../landice/tests/ismip7_forcing/__init__.py | 4 +- .../tests/ismip7_forcing/ismip7_forcing.cfg | 12 + .../ismip7_forcing/ismip7_forcing_test.cfg | 53 ++++ .../ismip7_forcing/ocean_thermal/__init__.py | 36 +++ .../ocean_thermal/process_thermal_forcing.py | 261 ++++++++++++++++++ 5 files changed, 365 insertions(+), 1 deletion(-) create mode 100644 compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg create mode 100644 compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py create mode 100644 compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py diff --git a/compass/landice/tests/ismip7_forcing/__init__.py b/compass/landice/tests/ismip7_forcing/__init__.py index 689f37030c..3f621a4c37 100644 --- a/compass/landice/tests/ismip7_forcing/__init__.py +++ b/compass/landice/tests/ismip7_forcing/__init__.py @@ -1,10 +1,11 @@ from compass.landice.tests.ismip7_forcing.atmosphere import Atmosphere +from compass.landice.tests.ismip7_forcing.ocean_thermal import OceanThermal from compass.testgroup import TestGroup class Ismip7Forcing(TestGroup): """ - A test group for processing ISMIP7 atmosphere forcing data + A test group for processing ISMIP7 forcing data for the Antarctic Ice Sheet """ @@ -20,3 +21,4 @@ def __init__(self, mpas_core): super().__init__(mpas_core=mpas_core, name="ismip7_forcing") self.add_test_case(Atmosphere(test_group=self)) + self.add_test_case(OceanThermal(test_group=self)) diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg index 75fd603ffb..1b2460f1a6 100644 --- a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg @@ -39,3 +39,15 @@ start_year = 1850 # End year for processing end_year = 2014 + +# config options for ismip7 AIS ocean thermal forcing +[ismip7_ais_ocean_thermal] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = bilinear + +# 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 new file mode 100644 index 0000000000..ff191d4c35 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg @@ -0,0 +1,53 @@ +# config options for ismip7 antarctic ice sheet forcing data +[ismip7_ais] + +# Base path to the input ISMIP7 forcing files. User has to supply. +base_path_ismip7 = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/forcing/ + +# Base path to the MALI mesh. User has to supply. +base_path_mali = /global/cfs/cdirs/fanssie/MALI_input_files/AIS_4to20km_r01/ + +# Base path to which output forcing files are saved. +output_base_path = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/test_processing + +# Name of climate model used to generate ISMIP7 forcing data. User has to supply. +# Available model names: CESM2-WACCM, MRI-ESM2-0 +model = CESM2-WACCM + +# Scenario for forcing data. User has to supply. +# Available scenarios: historical, ssp126, ssp370, ssp585 +scenario = historical + +# Name of the MALI mesh. Used to name mapping and output files. +mali_mesh_name = AIS_4to20km_r01_20220907 + +# MALI mesh file (e.g. Antarctic_8to80km_20220407.nc). User has to supply. +mali_mesh_file = AIS_4to20km_r01_20220907.nc + +# Number of MPI tasks for ESMF_RegridWeightGen +# Use 512 for 2km data sets. +esmf_ntasks = 512 + +# config options for ismip7 AIS atmosphere forcing +[ismip7_ais_atmosphere] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = conserve + +# Start year for processing +start_year = 1980 + +# End year for processing +end_year = 2014 + +# config options for ismip7 AIS ocean thermal forcing +[ismip7_ais_ocean_thermal] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = bilinear + +# Start year for processing +start_year = 1850 + +# End year for processing +end_year = 2014 diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py new file mode 100644 index 0000000000..56310b5cb4 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py @@ -0,0 +1,36 @@ +from compass.landice.tests.ismip7_forcing.configure import ( + configure as configure_testgroup, +) +from compass.landice.tests.ismip7_forcing.ocean_thermal.process_thermal_forcing import ( # noqa: E501 + ProcessThermalForcing, +) +from compass.testcase import TestCase + + +class OceanThermal(TestCase): + """ + A test case for processing ISMIP7 AIS ocean thermal forcing data. + Remaps annual 3D thermal forcing from the ISMIP7 8km polar + stereographic grid to the MALI unstructured mesh. + """ + + 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 = "ocean_thermal" + subdir = name + super().__init__(test_group=test_group, name=name, subdir=subdir) + + self.add_step(ProcessThermalForcing(test_case=self)) + + def configure(self): + """ + Configures test case + """ + configure_testgroup(config=self.config) 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 new file mode 100644 index 0000000000..488f83766f --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -0,0 +1,261 @@ +import glob +import os +import shutil + +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.step import Step + + +class ProcessThermalForcing(Step): + """ + A step for processing ISMIP7 ocean thermal forcing (tf) data. + Remaps annual 3D thermal forcing from the ISMIP7 8km polar + stereographic grid to the MALI unstructured mesh, preserving + the 30 vertical ocean layers. + """ + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_forcing.ocean_thermal.OceanThermal # noqa + The test case this step belongs to + """ + super().__init__(test_case=test_case, + name="process_thermal_forcing") + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config["ismip7_ais"] + 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_ais"] + 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") + + section = config["ismip7_ais_ocean_thermal"] + method_remap = section.get("method_remap") + start_year = section.getint("start_year") + end_year = section.getint("end_year") + + # Discover input files (decade-spanning files) + input_path = os.path.join(base_path_ismip7, "ocean", "tf", "v3") + file_pattern = (f"tf_AIS_{model}_{scenario}_" + f"ocean_v3_*.nc") + all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) + + if not all_files: + raise FileNotFoundError( + f"No ocean thermal forcing files found matching pattern:\n" + f" {os.path.join(input_path, file_pattern)}") + + # Filter to files that overlap with the requested year range. + # Files are named with decade ranges (e.g., 1850-1859). + input_files = [] + for f in all_files: + # Extract year range from filename (last part before .nc) + year_str = os.path.basename(f).split("_")[-1].replace(".nc", "") + parts = year_str.split("-") + file_start = int(parts[0]) + file_end = int(parts[1]) + if file_end >= start_year and file_start <= end_year: + input_files.append(f) + + if not input_files: + raise FileNotFoundError( + f"No ocean thermal forcing files for year range " + f"{start_year}-{end_year}") + + logger.info(f"Found {len(input_files)} ocean thermal forcing files " + f"overlapping years {start_year}-{end_year}") + + # Build mapping file using the first input file as grid template. + # Ocean grid (761x761, ~8km) differs from atmosphere (3041x3041, 2km). + mapping_file = (f"map_ismip7_ocean_8km_to_{mali_mesh_name}_" + f"{method_remap}.nc") + + if not os.path.exists(mapping_file): + logger.info("Building mapping file for ocean grid...") + build_mapping_file(config, logger, + input_files[0], mapping_file, + mali_mesh_file=mali_mesh_file, + method_remap=method_remap) + + # Remap each decade file + remapped_files = [] + for input_file in input_files: + basename = os.path.basename(input_file) + remapped_file = f"remapped_{basename}" + remapped_files.append(remapped_file) + + if os.path.exists(remapped_file): + logger.info(f" Remapped file exists, skipping: {basename}") + continue + + logger.info(f" Remapping: {basename}") + args = ["ncremap", + "-i", input_file, + "-o", remapped_file, + "-m", mapping_file, + "-v", "tf"] + + check_call(args, logger=logger) + + # 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}_" + f"{start_year}-{end_year}.nc") + + self._combine_and_rename(remapped_files, output_file, + start_year, end_year) + + # Clean up remapped files + logger.info("Cleaning up temporary remapped files...") + for f in remapped_files: + if os.path.exists(f): + os.remove(f) + + # Place output in appropriate directory + output_path = os.path.join(output_base_path, "ocean_thermal_forcing", + 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 _combine_and_rename(self, remapped_files, output_file, + start_year, end_year): + """ + Combine decade-spanning remapped files, subset to the requested + year range, and rename variables/dimensions to MALI conventions. + + Parameters + ---------- + remapped_files : list of str + List of remapped NetCDF file paths + + output_file : str + Output file path + + start_year : int + First year to include in output + + end_year : int + Last year to include in output + """ + ds = xr.open_mfdataset(remapped_files, concat_dim="time", + combine="nested", engine="netcdf4") + + # Subset to requested year range + years = ds.time.dt.year + ds = ds.sel(time=(years >= start_year) & (years <= end_year)) + + # Extract z coordinate and bounds before renaming + z_ocean = ds["z"] + z_bnds = ds["z_bnds"] + if "time" in z_bnds.dims: + z_bnds = z_bnds.isel(time=0) + + # Rename dimensions to MALI conventions + rename_dims = {} + if "time" in ds.dims: + rename_dims["time"] = "Time" + if "ncol" in ds.dims: + rename_dims["ncol"] = "nCells" + if "z" in ds.dims: + rename_dims["z"] = "nISMIP6OceanLayers" + if "bnds" in ds.dims: + rename_dims["bnds"] = "TWO" + ds = ds.rename(rename_dims) + + # Rename thermal forcing variable + if "tf" in ds: + ds = ds.rename({"tf": "ismip6shelfMelt_3dThermalForcing"}) + + # Set z coordinate and bounds as MALI-named variables + ds["ismip6shelfMelt_zOcean"] = ( + "nISMIP6OceanLayers", z_ocean.values) + ds["ismip6shelfMelt_zBndsOcean"] = ( + ("TWO", "nISMIP6OceanLayers"), z_bnds.values.T) + + # Transpose thermal forcing to MALI dimension order + # Registry: nISMIP6OceanLayers nCells Time (Fortran order) + # NetCDF (C order): Time, nCells, nISMIP6OceanLayers + ds["ismip6shelfMelt_3dThermalForcing"] = \ + ds["ismip6shelfMelt_3dThermalForcing"].transpose( + "Time", "nCells", "nISMIP6OceanLayers") + + # Add xtime variable with annual timestamps + xtime = [] + for t_index in range(ds.sizes["Time"]): + date = ds.Time[t_index] + yr = int(date.dt.year.values) + date_str = f"{yr:04d}-01-01_00:00:00".ljust(64) + xtime.append(date_str) + + ds["xtime"] = ("Time", xtime) + ds["xtime"] = ds.xtime.astype("S") + + # Set attributes + ds["ismip6shelfMelt_3dThermalForcing"].attrs = { + "long_name": "thermal forcing for ISMIP6 ice-shelf " + "melting method", + "units": "degC", + } + ds["ismip6shelfMelt_zOcean"].attrs = { + "long_name": "depth coordinate for ocean thermal forcing", + "units": "m", + } + ds["ismip6shelfMelt_zBndsOcean"].attrs = { + "long_name": "bounds for ISMIP6 ocean layers", + "units": "m", + } + + # Drop auxiliary variables from remapping + vars_to_drop = [v for v in ["lon", "lon_vertices", "lat", + "lat_vertices", "lon_bnds", "lat_bnds", + "area", "z_bnds", "time_bnds", + "x_bnds", "y_bnds"] + if v in ds] + if vars_to_drop: + ds = ds.drop_vars(vars_to_drop) + + # Also drop the renamed z coordinate if it persists + if "nISMIP6OceanLayers" in ds.coords: + ds = ds.drop_vars("nISMIP6OceanLayers") + + # Drop Time coordinate values (keep as dimension only) + if "Time" in ds.coords: + ds = ds.drop_vars("Time") + + write_netcdf(ds, output_file) From 4656c0012001be2a5040cddc3a228e461c14f5fd Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 4 Jun 2026 11:18:13 -0700 Subject: [PATCH 06/15] Add Greenland ISMIP7 processing --- .../ismip7_forcing/atmosphere/__init__.py | 14 +- .../atmosphere/process_runoff.py | 207 ++++++++++++++++++ .../ismip7_forcing/atmosphere/process_smb.py | 20 +- .../atmosphere/process_smb_gradient.py | 21 +- .../atmosphere/process_temperature.py | 20 +- .../process_temperature_gradient.py | 21 +- .../landice/tests/ismip7_forcing/configure.py | 7 +- .../tests/ismip7_forcing/create_mapfile.py | 20 +- .../tests/ismip7_forcing/ice_sheet_params.py | 47 ++++ .../tests/ismip7_forcing/ismip7_forcing.cfg | 16 +- .../ismip7_forcing/ocean_thermal/__init__.py | 8 +- .../ocean_thermal/process_thermal_forcing.py | 126 +++++++++-- 12 files changed, 461 insertions(+), 66 deletions(-) create mode 100644 compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py create mode 100644 compass/landice/tests/ismip7_forcing/ice_sheet_params.py diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py b/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py index 472e3f6d9b..d217008317 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py @@ -1,3 +1,6 @@ +from compass.landice.tests.ismip7_forcing.atmosphere.process_runoff import ( + ProcessRunoff, +) from compass.landice.tests.ismip7_forcing.atmosphere.process_smb import ( ProcessSmb, ) @@ -18,10 +21,10 @@ class Atmosphere(TestCase): """ - A test case for processing ISMIP7 AIS atmosphere forcing data. + A test case for processing ISMIP7 atmosphere forcing data. Remaps monthly SMB, temperature, and annual gradients (SMB and - temperature) from the ISMIP7 2km polar stereographic grid to the - MALI unstructured mesh. + temperature) from the ISMIP7 polar stereographic grid to the + MALI unstructured mesh. For GrIS, also processes runoff (mrro). """ def __init__(self, test_group): @@ -47,3 +50,8 @@ def configure(self): Configures test case """ configure_testgroup(config=self.config) + + # Add runoff step only for GrIS + ice_sheet = self.config.get("ismip7", "ice_sheet") + if ice_sheet == "gis": + self.add_step(ProcessRunoff(test_case=self)) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py new file mode 100644 index 0000000000..1321141539 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -0,0 +1,207 @@ +import glob +import os +import shutil + +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.ice_sheet_params import get_params +from compass.step import Step + + +class ProcessRunoff(Step): + """ + A step for processing ISMIP7 ice sheet runoff (mrro) data. + Remaps monthly runoff from the ISMIP7 polar stereographic + grid to the MALI unstructured mesh. GrIS only. + """ + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip7_forcing.atmosphere.Atmosphere + The test case this step belongs to + """ + super().__init__(test_case=test_case, name="process_runoff") + + 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 + params = get_params(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_atmosphere"] + method_remap = section.get("method_remap") + start_year = section.getint("start_year") + end_year = section.getint("end_year") + + # Discover input files + prefix = params['prefix'] + resolution = params['atm_resolution'] + version = params['atm_version'] + input_path = os.path.join(base_path_ismip7, "mrro", version) + file_pattern = (f"mrro_{prefix}_{model}_{scenario}_" + f"SDBN1-{resolution}_{version}_*.nc") + all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) + + if not all_files: + raise FileNotFoundError( + f"No runoff files found matching pattern:\n" + f" {os.path.join(input_path, file_pattern)}") + + # Filter to requested year range + input_files = [] + for f in all_files: + year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + if start_year <= year <= end_year: + input_files.append(f) + + if not input_files: + raise FileNotFoundError( + f"No runoff files found for year range " + f"{start_year}-{end_year}") + + logger.info(f"Found {len(input_files)} runoff files for years " + f"{start_year}-{end_year}") + + # Build mapping file (reuse if already created by other atm steps) + mapping_file = (f"map_ismip7_{ice_sheet}_atm_to_" + f"{mali_mesh_name}_{method_remap}.nc") + + if not os.path.exists(mapping_file): + logger.info("Building mapping file...") + build_mapping_file(config, logger, + input_files[0], mapping_file, + mali_mesh_file=mali_mesh_file, + method_remap=method_remap) + + # Remap each year file + remapped_files = [] + for input_file in input_files: + basename = os.path.basename(input_file) + remapped_file = f"remapped_{basename}" + remapped_files.append(remapped_file) + + if os.path.exists(remapped_file): + logger.info(f" Remapped file exists, skipping: {basename}") + continue + + logger.info(f" Remapping: {basename}") + args = ["ncremap", + "-i", input_file, + "-o", remapped_file, + "-m", mapping_file, + "-v", "mrro"] + + check_call(args, logger=logger) + + # Combine remapped files and rename to MALI conventions + logger.info("Combining remapped files and renaming variables...") + output_file = (f"{mali_mesh_name}_runoff_{model}_{scenario}_" + f"{start_year}-{end_year}.nc") + + self._combine_and_rename(remapped_files, output_file) + + # Clean up remapped files + logger.info("Cleaning up temporary remapped files...") + for f in remapped_files: + if os.path.exists(f): + os.remove(f) + + # Place output in appropriate directory + output_path = os.path.join(output_base_path, "atmosphere_forcing", + 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 _combine_and_rename(self, remapped_files, output_file): + """ + Combine yearly remapped files and rename variables/dimensions + to MALI conventions. + + Parameters + ---------- + remapped_files : list of str + List of remapped NetCDF file paths + + output_file : str + Output file path + """ + ds = xr.open_mfdataset(remapped_files, concat_dim="time", + combine="nested", engine="netcdf4") + + # Rename dimensions to MALI conventions + rename_dims = {} + if "time" in ds.dims: + rename_dims["time"] = "Time" + if "ncol" in ds.dims: + rename_dims["ncol"] = "nCells" + if rename_dims: + ds = ds.rename(rename_dims) + + # Rename variable + if "mrro" in ds: + ds = ds.rename({"mrro": "ismip6Runoff"}) + + # Add xtime variable with monthly timestamps + xtime = [] + for t_index in range(ds.sizes["Time"]): + date = ds.Time[t_index] + yr = int(date.dt.year.values) + mo = int(date.dt.month.values) + date_str = f"{yr:04d}-{mo:02d}-15_00:00:00".ljust(64) + xtime.append(date_str) + + ds["xtime"] = ("Time", xtime) + ds["xtime"] = ds.xtime.astype("S") + + # Set attributes + ds["ismip6Runoff"].attrs = { + "long_name": "ice sheet runoff", + "units": "kg m-2 s-1", + } + + # Drop auxiliary variables from remapping + vars_to_drop = [v for v in ["lon", "lon_vertices", "lat", + "lat_vertices", "area"] + if v in ds] + if vars_to_drop: + ds = ds.drop_vars(vars_to_drop) + + 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 e619767e25..02b96e95c2 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -9,6 +9,7 @@ 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.step import Step @@ -35,7 +36,7 @@ def setup(self): Set up this step of the test case """ config = self.config - section = config["ismip7_ais"] + section = config["ismip7"] base_path_mali = section.get("base_path_mali") mali_mesh_file = section.get("mali_mesh_file") @@ -49,8 +50,9 @@ def run(self): """ logger = self.logger config = self.config + params = get_params(config) - section = config["ismip7_ais"] + 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") @@ -58,14 +60,18 @@ def run(self): scenario = section.get("scenario") output_base_path = section.get("output_base_path") - section = config["ismip7_ais_atmosphere"] + section = config["ismip7_atmosphere"] method_remap = section.get("method_remap") start_year = section.getint("start_year") end_year = section.getint("end_year") # Discover input files - input_path = os.path.join(base_path_ismip7, "acabf", "v2") - file_pattern = f"acabf_AIS_{model}_{scenario}_SDBN1-2000m_v2_*.nc" + prefix = params['prefix'] + resolution = params['atm_resolution'] + version = params['atm_version'] + input_path = os.path.join(base_path_ismip7, "acabf", version) + file_pattern = (f"acabf_{prefix}_{model}_{scenario}_" + f"SDBN1-{resolution}_{version}_*.nc") all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) if not all_files: @@ -89,7 +95,9 @@ def run(self): f"{start_year}-{end_year}") # Build mapping file using the first input file as the grid template - mapping_file = f"map_ismip7_2km_to_{mali_mesh_name}_{method_remap}.nc" + ice_sheet = config.get("ismip7", "ice_sheet") + mapping_file = (f"map_ismip7_{ice_sheet}_atm_to_" + f"{mali_mesh_name}_{method_remap}.nc") if not os.path.exists(mapping_file): logger.info("Building mapping file...") 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 88fed545f1..bfb899431d 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -9,6 +9,7 @@ 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.step import Step @@ -36,7 +37,7 @@ def setup(self): Set up this step of the test case """ config = self.config - section = config["ismip7_ais"] + section = config["ismip7"] base_path_mali = section.get("base_path_mali") mali_mesh_file = section.get("mali_mesh_file") @@ -50,8 +51,9 @@ def run(self): """ logger = self.logger config = self.config + params = get_params(config) - section = config["ismip7_ais"] + 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") @@ -59,15 +61,18 @@ def run(self): scenario = section.get("scenario") output_base_path = section.get("output_base_path") - section = config["ismip7_ais_atmosphere"] + section = config["ismip7_atmosphere"] method_remap = section.get("method_remap") start_year = section.getint("start_year") end_year = section.getint("end_year") # Discover input files - input_path = os.path.join(base_path_ismip7, "dacabfdz", "v2") - file_pattern = (f"dacabfdz_AIS_{model}_{scenario}_" - f"SDBN1-2000m_v2_*.nc") + prefix = params['prefix'] + resolution = params['atm_resolution'] + version = params['atm_version'] + input_path = os.path.join(base_path_ismip7, "dacabfdz", version) + file_pattern = (f"dacabfdz_{prefix}_{model}_{scenario}_" + f"SDBN1-{resolution}_{version}_*.nc") all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) if not all_files: @@ -91,7 +96,9 @@ def run(self): f"{start_year}-{end_year}") # Build mapping file (reuse if already created by process_smb) - mapping_file = f"map_ismip7_2km_to_{mali_mesh_name}_{method_remap}.nc" + ice_sheet = config.get("ismip7", "ice_sheet") + mapping_file = (f"map_ismip7_{ice_sheet}_atm_to_" + f"{mali_mesh_name}_{method_remap}.nc") if not os.path.exists(mapping_file): logger.info("Building mapping file...") diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index d5b6fe77ce..2556b88227 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -9,6 +9,7 @@ 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.step import Step @@ -35,7 +36,7 @@ def setup(self): Set up this step of the test case """ config = self.config - section = config["ismip7_ais"] + section = config["ismip7"] base_path_mali = section.get("base_path_mali") mali_mesh_file = section.get("mali_mesh_file") @@ -49,8 +50,9 @@ def run(self): """ logger = self.logger config = self.config + params = get_params(config) - section = config["ismip7_ais"] + 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") @@ -58,14 +60,18 @@ def run(self): scenario = section.get("scenario") output_base_path = section.get("output_base_path") - section = config["ismip7_ais_atmosphere"] + section = config["ismip7_atmosphere"] method_remap = section.get("method_remap") start_year = section.getint("start_year") end_year = section.getint("end_year") # Discover input files - input_path = os.path.join(base_path_ismip7, "ts", "v2") - file_pattern = f"ts_AIS_{model}_{scenario}_SDBN1-2000m_v2_*.nc" + prefix = params['prefix'] + resolution = params['atm_resolution'] + version = params['atm_version'] + input_path = os.path.join(base_path_ismip7, "ts", version) + file_pattern = (f"ts_{prefix}_{model}_{scenario}_" + f"SDBN1-{resolution}_{version}_*.nc") all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) if not all_files: @@ -89,7 +95,9 @@ def run(self): f"{start_year}-{end_year}") # Build mapping file (reuse if already created by process_smb) - mapping_file = f"map_ismip7_2km_to_{mali_mesh_name}_{method_remap}.nc" + ice_sheet = config.get("ismip7", "ice_sheet") + mapping_file = (f"map_ismip7_{ice_sheet}_atm_to_" + f"{mali_mesh_name}_{method_remap}.nc") if not os.path.exists(mapping_file): logger.info("Building mapping file...") 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 9ed4fd4a99..2730252423 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -9,6 +9,7 @@ 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.step import Step @@ -37,7 +38,7 @@ def setup(self): Set up this step of the test case """ config = self.config - section = config["ismip7_ais"] + section = config["ismip7"] base_path_mali = section.get("base_path_mali") mali_mesh_file = section.get("mali_mesh_file") @@ -51,8 +52,9 @@ def run(self): """ logger = self.logger config = self.config + params = get_params(config) - section = config["ismip7_ais"] + 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") @@ -60,15 +62,18 @@ def run(self): scenario = section.get("scenario") output_base_path = section.get("output_base_path") - section = config["ismip7_ais_atmosphere"] + section = config["ismip7_atmosphere"] method_remap = section.get("method_remap") start_year = section.getint("start_year") end_year = section.getint("end_year") # Discover input files - input_path = os.path.join(base_path_ismip7, "dtsdz", "v2") - file_pattern = (f"dtsdz_AIS_{model}_{scenario}_" - f"SDBN1-2000m_v2_*.nc") + prefix = params['prefix'] + resolution = params['atm_resolution'] + version = params['atm_version'] + input_path = os.path.join(base_path_ismip7, "dtsdz", version) + file_pattern = (f"dtsdz_{prefix}_{model}_{scenario}_" + f"SDBN1-{resolution}_{version}_*.nc") all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) if not all_files: @@ -92,7 +97,9 @@ def run(self): f"for years {start_year}-{end_year}") # Build mapping file (reuse if already created by other steps) - mapping_file = f"map_ismip7_2km_to_{mali_mesh_name}_{method_remap}.nc" + ice_sheet = config.get("ismip7", "ice_sheet") + mapping_file = (f"map_ismip7_{ice_sheet}_atm_to_" + f"{mali_mesh_name}_{method_remap}.nc") if not os.path.exists(mapping_file): logger.info("Building mapping file...") diff --git a/compass/landice/tests/ismip7_forcing/configure.py b/compass/landice/tests/ismip7_forcing/configure.py index c49370206c..8719c0e141 100644 --- a/compass/landice/tests/ismip7_forcing/configure.py +++ b/compass/landice/tests/ismip7_forcing/configure.py @@ -9,9 +9,10 @@ def configure(config): Configuration options for an ismip7 forcing test case """ - section = "ismip7_ais" - options = ["base_path_ismip7", "base_path_mali", "mali_mesh_name", - "mali_mesh_file", "output_base_path", "model", "scenario"] + section = "ismip7" + options = ["ice_sheet", "base_path_ismip7", "base_path_mali", + "mali_mesh_name", "mali_mesh_file", "output_base_path", + "model", "scenario"] for option in options: value = config.get(section=section, option=option) diff --git a/compass/landice/tests/ismip7_forcing/create_mapfile.py b/compass/landice/tests/ismip7_forcing/create_mapfile.py index 81ea411ba6..de10f32b34 100644 --- a/compass/landice/tests/ismip7_forcing/create_mapfile.py +++ b/compass/landice/tests/ismip7_forcing/create_mapfile.py @@ -7,9 +7,9 @@ def build_mapping_file(config, logger, ismip7_grid_file, mapping_file, mali_mesh_file=None, - method_remap=None): + method_remap=None, projection=None): """ - Build a mapping file for regridding from the ISMIP7 2km polar + Build a mapping file for regridding from an ISMIP7 polar stereographic grid to the MALI unstructured mesh. Parameters @@ -31,6 +31,10 @@ def build_mapping_file(config, logger, ismip7_grid_file, method_remap : str, optional Remapping method: 'bilinear', 'neareststod', or 'conserve' + + projection : str, optional + Projection flag for SCRIP generation (e.g., 'ais-bedmap2', + 'gis-bamber'). If not provided, reads from ice_sheet_params. """ if os.path.exists(mapping_file): @@ -48,8 +52,14 @@ def build_mapping_file(config, logger, ismip7_grid_file, raise ValueError("Remapping method must be provided. " "Options: 'bilinear', 'neareststod', 'conserve'.") - # AIS polar stereographic projection (EPSG:3031) - ismip7_projection = "ais-bedmap2" + # Determine projection from parameter or config + if projection is None: + from compass.landice.tests.ismip7_forcing.ice_sheet_params import ( + get_params, + ) + projection = get_params(config)['projection'] + + ismip7_projection = projection # name temporary scrip files source_grid_scripfile = "temp_source_scrip.nc" @@ -81,7 +91,7 @@ def build_mapping_file(config, logger, ismip7_grid_file, # create a mapping file using ESMF_RegridWeightGen logger.info(f"Creating mapping file with method: {method_remap}") - section = config["ismip7_ais"] + section = config["ismip7"] cores = section.getint("esmf_ntasks") parallel_executable = config.get("parallel", "parallel_executable") diff --git a/compass/landice/tests/ismip7_forcing/ice_sheet_params.py b/compass/landice/tests/ismip7_forcing/ice_sheet_params.py new file mode 100644 index 0000000000..8324236c56 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ice_sheet_params.py @@ -0,0 +1,47 @@ +""" +Ice-sheet-specific parameters for ISMIP7 forcing data processing. +""" + +# Parameters that differ between AIS and GrIS +_PARAMS = { + 'ais': { + 'projection': 'ais-bedmap2', + 'prefix': 'AIS', + 'atm_resolution': '2000m', + 'atm_version': 'v2', + 'ocean_version': 'v3', + 'ocean_3d': True, + 'ocean_temporal': 'decade', + }, + 'gis': { + 'projection': 'gis-bamber', + 'prefix': 'GrIS', + 'atm_resolution': '1000m', + 'atm_version': 'v2', + 'ocean_version': 'v2', + 'ocean_3d': False, + 'ocean_temporal': 'yearly', + }, +} + + +def get_params(config): + """ + Get ice-sheet-specific parameters from the config. + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for the test case + + Returns + ------- + params : dict + Dictionary of ice-sheet-specific parameters + """ + ice_sheet = config.get("ismip7", "ice_sheet") + if ice_sheet not in _PARAMS: + raise ValueError( + f"Unknown ice_sheet '{ice_sheet}'. " + f"Must be one of: {list(_PARAMS.keys())}") + return _PARAMS[ice_sheet] diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg index 1b2460f1a6..e7bf3c4e6a 100644 --- a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg @@ -1,5 +1,8 @@ -# config options for ismip7 antarctic ice sheet forcing data -[ismip7_ais] +# config options for ismip7 forcing data +[ismip7] + +# Ice sheet: ais (Antarctic) or gis (Greenland) +ice_sheet = NotAvailable # Base path to the input ISMIP7 forcing files. User has to supply. base_path_ismip7 = NotAvailable @@ -25,11 +28,10 @@ mali_mesh_name = NotAvailable mali_mesh_file = NotAvailable # Number of MPI tasks for ESMF_RegridWeightGen -# Use 512 for 2km data sets. esmf_ntasks = 128 -# config options for ismip7 AIS atmosphere forcing -[ismip7_ais_atmosphere] +# config options for ismip7 atmosphere forcing +[ismip7_atmosphere] # Remapping method. Options: bilinear, neareststod, conserve method_remap = conserve @@ -40,8 +42,8 @@ start_year = 1850 # End year for processing end_year = 2014 -# config options for ismip7 AIS ocean thermal forcing -[ismip7_ais_ocean_thermal] +# config options for ismip7 ocean thermal forcing +[ismip7_ocean_thermal] # Remapping method. Options: bilinear, neareststod, conserve method_remap = bilinear diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py index 56310b5cb4..2a97700428 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py @@ -9,9 +9,11 @@ class OceanThermal(TestCase): """ - A test case for processing ISMIP7 AIS ocean thermal forcing data. - Remaps annual 3D thermal forcing from the ISMIP7 8km polar - stereographic grid to the MALI unstructured mesh. + A test case for processing ISMIP7 ocean thermal forcing data. + For AIS: Remaps annual 3D thermal forcing from the ISMIP7 8km + polar stereographic grid to the MALI mesh. + For GrIS: Remaps monthly 2D thermal forcing from the ISMIP7 1km + grid to the MALI mesh. """ def __init__(self, test_group): 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 488f83766f..2f2ec25698 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 @@ -9,15 +9,18 @@ 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.step import Step class ProcessThermalForcing(Step): """ A step for processing ISMIP7 ocean thermal forcing (tf) data. - Remaps annual 3D thermal forcing from the ISMIP7 8km polar + For AIS: Remaps annual 3D thermal forcing from the ISMIP7 8km polar stereographic grid to the MALI unstructured mesh, preserving the 30 vertical ocean layers. + For GrIS: Remaps monthly 2D thermal forcing from the ISMIP7 1km + grid to the MALI unstructured mesh. """ def __init__(self, test_case): @@ -37,7 +40,7 @@ def setup(self): Set up this step of the test case """ config = self.config - section = config["ismip7_ais"] + section = config["ismip7"] base_path_mali = section.get("base_path_mali") mali_mesh_file = section.get("mali_mesh_file") @@ -51,24 +54,30 @@ def run(self): """ logger = self.logger config = self.config + params = get_params(config) - section = config["ismip7_ais"] + 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_ais_ocean_thermal"] + section = config["ismip7_ocean_thermal"] method_remap = section.get("method_remap") start_year = section.getint("start_year") end_year = section.getint("end_year") - # Discover input files (decade-spanning files) - input_path = os.path.join(base_path_ismip7, "ocean", "tf", "v3") - file_pattern = (f"tf_AIS_{model}_{scenario}_" - f"ocean_v3_*.nc") + # Discover input files + prefix = params['prefix'] + ocean_version = params['ocean_version'] + ocean_3d = params['ocean_3d'] + input_path = os.path.join(base_path_ismip7, "ocean", "tf", + ocean_version) + file_pattern = (f"tf_{prefix}_{model}_{scenario}_" + f"ocean_{ocean_version}_*.nc") all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) if not all_files: @@ -77,14 +86,15 @@ def run(self): f" {os.path.join(input_path, file_pattern)}") # Filter to files that overlap with the requested year range. - # Files are named with decade ranges (e.g., 1850-1859). + # AIS files are named with decade ranges (e.g., 1850-1859). + # GrIS files are named with single years (e.g., 2015). input_files = [] for f in all_files: # Extract year range from filename (last part before .nc) year_str = os.path.basename(f).split("_")[-1].replace(".nc", "") parts = year_str.split("-") file_start = int(parts[0]) - file_end = int(parts[1]) + file_end = int(parts[-1]) # same as start for single-year files if file_end >= start_year and file_start <= end_year: input_files.append(f) @@ -97,9 +107,8 @@ def run(self): f"overlapping years {start_year}-{end_year}") # Build mapping file using the first input file as grid template. - # Ocean grid (761x761, ~8km) differs from atmosphere (3041x3041, 2km). - mapping_file = (f"map_ismip7_ocean_8km_to_{mali_mesh_name}_" - f"{method_remap}.nc") + mapping_file = (f"map_ismip7_{ice_sheet}_ocean_to_" + f"{mali_mesh_name}_{method_remap}.nc") if not os.path.exists(mapping_file): logger.info("Building mapping file for ocean grid...") @@ -133,8 +142,12 @@ def run(self): output_file = (f"{mali_mesh_name}_thermal_forcing_{model}_{scenario}_" f"{start_year}-{end_year}.nc") - self._combine_and_rename(remapped_files, output_file, - start_year, end_year) + if ocean_3d: + self._combine_and_rename_3d(remapped_files, output_file, + start_year, end_year) + else: + self._combine_and_rename_2d(remapped_files, output_file, + start_year, end_year) # Clean up remapped files logger.info("Cleaning up temporary remapped files...") @@ -153,11 +166,12 @@ def run(self): logger.info(f"Done. Output: {dst}") - def _combine_and_rename(self, remapped_files, output_file, - start_year, end_year): + def _combine_and_rename_3d(self, remapped_files, output_file, + start_year, end_year): """ - Combine decade-spanning remapped files, subset to the requested - year range, and rename variables/dimensions to MALI conventions. + Combine decade-spanning remapped files (AIS), subset to the + requested year range, and rename variables/dimensions to MALI + conventions for 3D thermal forcing. Parameters ---------- @@ -259,3 +273,77 @@ def _combine_and_rename(self, remapped_files, output_file, ds = ds.drop_vars("Time") write_netcdf(ds, output_file) + + def _combine_and_rename_2d(self, remapped_files, output_file, + start_year, end_year): + """ + Combine yearly remapped files (GrIS), subset to the requested + year range, and rename variables/dimensions to MALI conventions + for 2D thermal forcing. + + Parameters + ---------- + remapped_files : list of str + List of remapped NetCDF file paths + + output_file : str + Output file path + + start_year : int + First year to include in output + + end_year : int + Last year to include in output + """ + ds = xr.open_mfdataset(remapped_files, concat_dim="time", + combine="nested", engine="netcdf4") + + # Subset to requested year range + years = ds.time.dt.year + ds = ds.sel(time=(years >= start_year) & (years <= end_year)) + + # Rename dimensions to MALI conventions + rename_dims = {} + if "time" in ds.dims: + rename_dims["time"] = "Time" + if "ncol" in ds.dims: + rename_dims["ncol"] = "nCells" + if rename_dims: + ds = ds.rename(rename_dims) + + # Rename thermal forcing variable + if "tf" in ds: + ds = ds.rename({"tf": "ismip6_2dThermalForcing"}) + + # Add xtime variable with monthly timestamps + xtime = [] + for t_index in range(ds.sizes["Time"]): + date = ds.Time[t_index] + yr = int(date.dt.year.values) + mo = int(date.dt.month.values) + date_str = f"{yr:04d}-{mo:02d}-15_00:00:00".ljust(64) + xtime.append(date_str) + + ds["xtime"] = ("Time", xtime) + ds["xtime"] = ds.xtime.astype("S") + + # Set attributes + ds["ismip6_2dThermalForcing"].attrs = { + "long_name": "2D thermal forcing for ISMIP6 ice-shelf " + "melting parameterization", + "units": "degC", + } + + # Drop auxiliary variables from remapping + vars_to_drop = [v for v in ["lon", "lon_vertices", "lat", + "lat_vertices", "area", + "time_bnds", "x_bnds", "y_bnds"] + if v in ds] + if vars_to_drop: + ds = ds.drop_vars(vars_to_drop) + + # Drop Time coordinate values (keep as dimension only) + if "Time" in ds.coords: + ds = ds.drop_vars("Time") + + write_netcdf(ds, output_file) From 4921684a17a4a04a3aff1b17b737131a73b6faa6 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 4 Jun 2026 13:13:43 -0700 Subject: [PATCH 07/15] Add docs for ismip7_forcing test group --- docs/developers_guide/landice/api.rst | 37 ++++ .../landice/test_groups/index.rst | 1 + .../landice/test_groups/ismip7_forcing.rst | 103 +++++++++ .../users_guide/landice/test_groups/index.rst | 1 + .../landice/test_groups/ismip7_forcing.rst | 197 ++++++++++++++++++ 5 files changed, 339 insertions(+) create mode 100644 docs/developers_guide/landice/test_groups/ismip7_forcing.rst create mode 100644 docs/users_guide/landice/test_groups/ismip7_forcing.rst diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index 0aa5c4406d..57514fd526 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -373,6 +373,43 @@ ismip6_run ismip6_ais_proj2300.set_up_experiment.SetUpExperiment.setup ismip6_ais_proj2300.set_up_experiment.SetUpExperiment.run +ismip7_forcing +~~~~~~~~~~~~~~ + +.. currentmodule:: compass.landice.tests.ismip7_forcing + +.. autosummary:: + :toctree: generated/ + + Ismip7Forcing + configure.configure + ice_sheet_params.get_params + create_mapfile.build_mapping_file + + atmosphere.Atmosphere + atmosphere.Atmosphere.configure + atmosphere.process_smb.ProcessSmb + atmosphere.process_smb.ProcessSmb.setup + atmosphere.process_smb.ProcessSmb.run + atmosphere.process_temperature.ProcessTemperature + atmosphere.process_temperature.ProcessTemperature.setup + atmosphere.process_temperature.ProcessTemperature.run + atmosphere.process_smb_gradient.ProcessSmbGradient + atmosphere.process_smb_gradient.ProcessSmbGradient.setup + atmosphere.process_smb_gradient.ProcessSmbGradient.run + atmosphere.process_temperature_gradient.ProcessTemperatureGradient + atmosphere.process_temperature_gradient.ProcessTemperatureGradient.setup + atmosphere.process_temperature_gradient.ProcessTemperatureGradient.run + atmosphere.process_runoff.ProcessRunoff + atmosphere.process_runoff.ProcessRunoff.setup + atmosphere.process_runoff.ProcessRunoff.run + + ocean_thermal.OceanThermal + ocean_thermal.OceanThermal.configure + ocean_thermal.process_thermal_forcing.ProcessThermalForcing + ocean_thermal.process_thermal_forcing.ProcessThermalForcing.setup + ocean_thermal.process_thermal_forcing.ProcessThermalForcing.run + isunnguata_sermia ~~~~~~~~~~~~~~~~~ diff --git a/docs/developers_guide/landice/test_groups/index.rst b/docs/developers_guide/landice/test_groups/index.rst index 0f4ae50afa..2556037c52 100644 --- a/docs/developers_guide/landice/test_groups/index.rst +++ b/docs/developers_guide/landice/test_groups/index.rst @@ -20,6 +20,7 @@ Test groups hydro_radial ismip6_forcing ismip6_run + ismip7_forcing isunnguata_sermia kangerlussuaq koge_bugt_s diff --git a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst new file mode 100644 index 0000000000..f4d321dbb2 --- /dev/null +++ b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst @@ -0,0 +1,103 @@ +.. _dev_landice_ismip7_forcing: + +ismip7_forcing +============== + +The ``ismip7_forcing`` test group +(:py:class:`compass.landice.tests.ismip7_forcing.Ismip7Forcing`) processes +(i.e., remaps and renames) the atmospheric and ocean thermal forcing data of +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``. + +.. _dev_landice_ismip7_forcing_framework: + +framework +--------- + +The shared config options for the ``ismip7_forcing`` test group are described +in :ref:`landice_ismip7_forcing` in the User's Guide. + +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 +:py:func:`compass.landice.tests.ismip7_forcing.ice_sheet_params.get_params` +to retrieve them based on the ``ice_sheet`` config option. + +configure +~~~~~~~~~ + +The module :py:mod:`compass.landice.tests.ismip7_forcing.configure` validates +that all required config options in the ``[ismip7]`` section have been set by +the user (i.e., are not ``NotAvailable``). + +create_mapfile +~~~~~~~~~~~~~~ + +The module :py:mod:`compass.landice.tests.ismip7_forcing.create_mapfile` +defines a unified framework for creating SCRIP and mapping files. The function +:py:func:`compass.landice.tests.ismip7_forcing.create_mapfile.build_mapping_file` +creates a SCRIP file from the input polar stereographic grid using the +``create_scrip_file_from_planar_rectangular_grid`` command from MPAS-Tools, +then generates a mapping file via ``ESMF_RegridWeightGen``. The projection +is automatically determined from the ``ice_sheet`` config option using +``ice_sheet_params``. + +Test cases +---------- + +.. _dev_landice_ismip7_forcing_atmosphere: + +atmosphere +~~~~~~~~~~ + +The :py:class:`compass.landice.tests.ismip7_forcing.atmosphere.Atmosphere` +test case processes the ISMIP7 atmosphere forcing fields. It contains four +steps that are always included (SMB, temperature, and their respective +gradients) plus a conditional ``process_runoff`` step that is added only when +``ice_sheet = gis``. Each step discovers input files matching the +ice-sheet-specific naming pattern, builds or reuses a mapping file, remaps +each input file with ``ncremap``, and combines/renames the results to MALI +conventions. + +Steps: + +* :py:class:`~compass.landice.tests.ismip7_forcing.atmosphere.process_smb.ProcessSmb` — + ``acabf`` → ``sfcMassBal`` +* :py:class:`~compass.landice.tests.ismip7_forcing.atmosphere.process_temperature.ProcessTemperature` — + ``ts`` → ``surfaceAirTemperature`` (clipped ≤ 273.15 K) +* :py:class:`~compass.landice.tests.ismip7_forcing.atmosphere.process_smb_gradient.ProcessSmbGradient` — + ``dacabfdz`` → ``sfcMassBalLapseRate`` +* :py:class:`~compass.landice.tests.ismip7_forcing.atmosphere.process_temperature_gradient.ProcessTemperatureGradient` — + ``dtsdz`` → ``surfaceAirTemperatureLapseRate`` +* :py:class:`~compass.landice.tests.ismip7_forcing.atmosphere.process_runoff.ProcessRunoff` — + ``mrro`` → ``ismip6Runoff`` (GrIS only) + +.. _dev_landice_ismip7_forcing_ocean_thermal: + +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, +: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``. + +For AIS, the step: + +* Remaps thermal forcing preserving 30 vertical ocean layers +* Produces ``ismip6shelfMelt_3dThermalForcing`` (dims: Time × nCells × + nISMIP6OceanLayers) +* Includes depth coordinate variables ``ismip6shelfMelt_zOcean`` and + ``ismip6shelfMelt_zBndsOcean`` + +For GrIS, the step: + +* Remaps 2D monthly thermal forcing +* Produces ``ismip6_2dThermalForcing`` (dims: Time × nCells) diff --git a/docs/users_guide/landice/test_groups/index.rst b/docs/users_guide/landice/test_groups/index.rst index 80fed7af3a..f744061134 100644 --- a/docs/users_guide/landice/test_groups/index.rst +++ b/docs/users_guide/landice/test_groups/index.rst @@ -25,6 +25,7 @@ physics but that are not run routinely. hydro_radial ismip6_forcing ismip6_run + ismip7_forcing isunnguata_sermia kangerlussuaq koge_bugt_s diff --git a/docs/users_guide/landice/test_groups/ismip7_forcing.rst b/docs/users_guide/landice/test_groups/ismip7_forcing.rst new file mode 100644 index 0000000000..b5963d36dd --- /dev/null +++ b/docs/users_guide/landice/test_groups/ismip7_forcing.rst @@ -0,0 +1,197 @@ +.. _landice_ismip7_forcing: + +ismip7_forcing +============== + +The ``landice/ismip7_forcing`` test group processes (i.e., remaps and renames) +the atmospheric and ocean thermal forcing data of the Ice Sheet Model +Intercomparison for CMIP7 (ISMIP7) protocol. The processed data is used to +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 ``atmosphere`` test case has four steps for both ice sheets: + ``process_smb``, ``process_temperature``, ``process_smb_gradient``, and + ``process_temperature_gradient``. For GrIS, a fifth step ``process_runoff`` + is additionally included. + +* 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. + +(For more details on the steps of each test case, see +:ref:`landice_ismip7_forcing_atmosphere` and +:ref:`landice_ismip7_forcing_ocean_thermal`.) + +.. _landice_ismip7_forcing_usage: + +Usage +----- + +To use this test group, users need to: + +1. Provide a MALI mesh file onto which the source data will be remapped. + +2. Set the ``ice_sheet`` config option to either ``ais`` or ``gis``. + +3. Provide the path to the ISMIP7 forcing data (``base_path_ismip7``). + +4. Run the ``atmosphere`` test case for each model and scenario combination. + +5. Run the ``ocean_thermal`` test case for each model and scenario combination. + +.. _landice_ismip7_forcing_input_data: + +Input Data +---------- + +ISMIP7 forcing data is organized by variable and version. The expected +directory structure under ``base_path_ismip7`` is: + +For AIS atmosphere (2km, polar stereographic EPSG:3031): + +.. code-block:: none + + acabf/v2/acabf_AIS_{model}_{scenario}_SDBN1-2000m_v2_{year_range}.nc + ts/v2/ts_AIS_{model}_{scenario}_SDBN1-2000m_v2_{year_range}.nc + dacabfdz/v2/dacabfdz_AIS_{model}_{scenario}_SDBN1-2000m_v2_{year_range}.nc + dtsdz/v2/dtsdz_AIS_{model}_{scenario}_SDBN1-2000m_v2_{year_range}.nc + +For AIS ocean thermal (8km, 30 depth levels, decade files): + +.. code-block:: none + + ocean/tf/v3/tf_AIS_{model}_{scenario}_ocean_v3_{start_year}-{end_year}.nc + +For GrIS atmosphere (1km, polar stereographic EPSG:3413): + +.. code-block:: none + + acabf/v2/acabf_GrIS_{model}_{scenario}_SDBN1-1000m_v2_{year}.nc + ts/v2/ts_GrIS_{model}_{scenario}_SDBN1-1000m_v2_{year}.nc + dacabfdz/v2/dacabfdz_GrIS_{model}_{scenario}_SDBN1-1000m_v2_{year}.nc + dtsdz/v2/dtsdz_GrIS_{model}_{scenario}_SDBN1-1000m_v2_{year}.nc + mrro/v2/mrro_GrIS_{model}_{scenario}_SDBN1-1000m_v2_{year}.nc + +For GrIS ocean thermal (same 1km grid, 2D, yearly files): + +.. code-block:: none + + ocean/tf/v2/tf_GrIS_{model}_{scenario}_ocean_v2_{year}.nc + +.. _landice_ismip7_forcing_config: + +config options +-------------- + +The ``ismip7_forcing`` test group uses three config sections. The default +values are: + +.. code-block:: cfg + + # config options for ismip7 forcing data + [ismip7] + + # Ice sheet: ais (Antarctic) or gis (Greenland) + ice_sheet = NotAvailable + + # Base path to the input ISMIP7 forcing files + base_path_ismip7 = NotAvailable + + # Base path to the MALI mesh + base_path_mali = NotAvailable + + # Base path to which output forcing files are saved + output_base_path = NotAvailable + + # Name of climate model (e.g., CESM2-WACCM, MRI-ESM2-0) + model = NotAvailable + + # Scenario (e.g., historical, ssp126, ssp370, ssp585) + scenario = NotAvailable + + # Name of the MALI mesh (used in output file naming) + mali_mesh_name = NotAvailable + + # MALI mesh file name + mali_mesh_file = NotAvailable + + # Number of MPI tasks for ESMF_RegridWeightGen + esmf_ntasks = 128 + + # config options for ismip7 atmosphere forcing + [ismip7_atmosphere] + + # Remapping method: bilinear, neareststod, conserve + method_remap = conserve + + # Start year for processing + start_year = 1850 + + # End year for processing + end_year = 2014 + + # config options for ismip7 ocean thermal forcing + [ismip7_ocean_thermal] + + # Remapping method: bilinear, neareststod, conserve + method_remap = bilinear + + # 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``). + +.. _landice_ismip7_forcing_atmosphere: + +atmosphere +---------- + +The ``landice/ismip7_forcing/atmosphere`` test case processes the ISMIP7 +atmosphere forcing fields and remaps them from the native polar stereographic +grid to the MALI unstructured mesh. + +Steps: + +* **process_smb**: Remaps the surface mass balance (``acabf``) field. The + output variable is ``sfcMassBal``. + +* **process_temperature**: Remaps the surface temperature (``ts``) field, + clipped to a maximum of 273.15 K. The output variable is + ``surfaceAirTemperature``. + +* **process_smb_gradient**: Remaps the SMB lapse rate (``dacabfdz``) field. + The output variable is ``sfcMassBalLapseRate``. + +* **process_temperature_gradient**: Remaps the temperature lapse rate + (``dtsdz``) field. The output variable is + ``surfaceAirTemperatureLapseRate``. + +* **process_runoff** (GrIS only): Remaps the ice sheet runoff (``mrro``) + field. The output variable is ``ismip6Runoff``. + +.. _landice_ismip7_forcing_ocean_thermal: + +ocean_thermal +------------- + +The ``landice/ismip7_forcing/ocean_thermal`` test case processes the ISMIP7 +ocean thermal forcing (``tf``) and remaps it from the native polar +stereographic grid to the MALI unstructured mesh. + +For **AIS**, thermal forcing is 3D with 30 vertical ocean layers. The input +files span decades (e.g., 1850-1859). The output variable is +``ismip6shelfMelt_3dThermalForcing`` with dimension +``nISMIP6OceanLayers``. Associated depth coordinate variables +``ismip6shelfMelt_zOcean`` and ``ismip6shelfMelt_zBndsOcean`` are also +produced. + +For **GrIS**, thermal forcing is 2D (depth-averaged), with monthly temporal +resolution and yearly input files. The output variable is +``ismip6_2dThermalForcing``. From 1b136cb41392a1ffc810fb42e3cdb9e3c75aea0f Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 4 Jun 2026 14:20:10 -0700 Subject: [PATCH 08/15] Extrapolate 2D thermal forcing on source grid before remapping Extrapolate 2D thermal forcing on source grid before remapping to MALI mesh. This gets rid of missing values that were being interpolated to the MALI mesh as 1e36. --- .../ocean_thermal/process_thermal_forcing.py | 70 ++++++++++++++++++- 1 file changed, 69 insertions(+), 1 deletion(-) 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 2f2ec25698..2b24ebdba3 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,9 +2,11 @@ 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, @@ -128,15 +130,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): + self._extrapolate_source(input_file, extrap_file, "tf", + logger) + logger.info(f" Remapping: {basename}") args = ["ncremap", - "-i", input_file, + "-i", extrap_file, "-o", remapped_file, "-m", mapping_file, "-v", "tf"] 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}_thermal_forcing_{model}_{scenario}_" @@ -347,3 +359,59 @@ def _combine_and_rename_2d(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 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) From f6f6a2f2d4039cc1676df8ef8707e9321d673f5c Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 4 Jun 2026 14:39:25 -0700 Subject: [PATCH 09/15] Process runoff for AIS Process runoff field for AIS, since we have it available and can use it in face-melting. --- .../ismip7_forcing/atmosphere/__init__.py | 8 +-- .../atmosphere/process_runoff.py | 67 ++++++++++++++++++- .../landice/test_groups/ismip7_forcing.rst | 14 ++-- .../landice/test_groups/ismip7_forcing.rst | 10 +-- 4 files changed, 79 insertions(+), 20 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py b/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py index d217008317..7e577b2026 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py @@ -24,7 +24,7 @@ class Atmosphere(TestCase): A test case for processing ISMIP7 atmosphere forcing data. Remaps monthly SMB, temperature, and annual gradients (SMB and temperature) from the ISMIP7 polar stereographic grid to the - MALI unstructured mesh. For GrIS, also processes runoff (mrro). + MALI unstructured mesh. Also processes runoff (mrro). """ def __init__(self, test_group): @@ -44,14 +44,10 @@ def __init__(self, test_group): self.add_step(ProcessTemperature(test_case=self)) self.add_step(ProcessSmbGradient(test_case=self)) self.add_step(ProcessTemperatureGradient(test_case=self)) + self.add_step(ProcessRunoff(test_case=self)) def configure(self): """ Configures test case """ configure_testgroup(config=self.config) - - # Add runoff step only for GrIS - ice_sheet = self.config.get("ismip7", "ice_sheet") - if ice_sheet == "gis": - self.add_step(ProcessRunoff(test_case=self)) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py index 1321141539..2fec1a9fd0 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -2,9 +2,11 @@ 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, @@ -117,15 +119,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): + self._extrapolate_source(input_file, extrap_file, "mrro", + logger) + logger.info(f" Remapping: {basename}") args = ["ncremap", - "-i", input_file, + "-i", extrap_file, "-o", remapped_file, "-m", mapping_file, "-v", "mrro"] 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}_runoff_{model}_{scenario}_" @@ -205,3 +217,56 @@ def _combine_and_rename(self, remapped_files, output_file): ds = ds.drop_vars(vars_to_drop) 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/docs/developers_guide/landice/test_groups/ismip7_forcing.rst b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst index f4d321dbb2..2f4349e804 100644 --- a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst @@ -57,13 +57,11 @@ atmosphere ~~~~~~~~~~ The :py:class:`compass.landice.tests.ismip7_forcing.atmosphere.Atmosphere` -test case processes the ISMIP7 atmosphere forcing fields. It contains four -steps that are always included (SMB, temperature, and their respective -gradients) plus a conditional ``process_runoff`` step that is added only when -``ice_sheet = gis``. Each step discovers input files matching the -ice-sheet-specific naming pattern, builds or reuses a mapping file, remaps -each input file with ``ncremap``, and combines/renames the results to MALI -conventions. +test case processes the ISMIP7 atmosphere forcing fields. It contains five +steps: SMB, temperature, their respective gradients, and runoff. Each step +discovers input files matching the ice-sheet-specific naming pattern, builds +or reuses a mapping file, remaps each input file with ``ncremap``, and +combines/renames the results to MALI conventions. Steps: @@ -76,7 +74,7 @@ Steps: * :py:class:`~compass.landice.tests.ismip7_forcing.atmosphere.process_temperature_gradient.ProcessTemperatureGradient` — ``dtsdz`` → ``surfaceAirTemperatureLapseRate`` * :py:class:`~compass.landice.tests.ismip7_forcing.atmosphere.process_runoff.ProcessRunoff` — - ``mrro`` → ``ismip6Runoff`` (GrIS only) + ``mrro`` → ``ismip6Runoff`` .. _dev_landice_ismip7_forcing_ocean_thermal: diff --git a/docs/users_guide/landice/test_groups/ismip7_forcing.rst b/docs/users_guide/landice/test_groups/ismip7_forcing.rst index b5963d36dd..c797bf5dbf 100644 --- a/docs/users_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/users_guide/landice/test_groups/ismip7_forcing.rst @@ -12,10 +12,9 @@ Ice Sheet (GrIS), controlled by a single ``ice_sheet`` config option. The test group includes two test cases: ``atmosphere`` and ``ocean_thermal``. -* The ``atmosphere`` test case has four steps for both ice sheets: - ``process_smb``, ``process_temperature``, ``process_smb_gradient``, and - ``process_temperature_gradient``. For GrIS, a fifth step ``process_runoff`` - is additionally included. +* The ``atmosphere`` test case has five steps: + ``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 @@ -58,6 +57,7 @@ For AIS atmosphere (2km, polar stereographic EPSG:3031): ts/v2/ts_AIS_{model}_{scenario}_SDBN1-2000m_v2_{year_range}.nc dacabfdz/v2/dacabfdz_AIS_{model}_{scenario}_SDBN1-2000m_v2_{year_range}.nc dtsdz/v2/dtsdz_AIS_{model}_{scenario}_SDBN1-2000m_v2_{year_range}.nc + mrro/v2/mrro_AIS_{model}_{scenario}_SDBN1-2000m_v2_{year_range}.nc For AIS ocean thermal (8km, 30 depth levels, decade files): @@ -173,7 +173,7 @@ Steps: (``dtsdz``) field. The output variable is ``surfaceAirTemperatureLapseRate``. -* **process_runoff** (GrIS only): Remaps the ice sheet runoff (``mrro``) +* **process_runoff**: Remaps the ice sheet runoff (``mrro``) field. The output variable is ``ismip6Runoff``. .. _landice_ismip7_forcing_ocean_thermal: From 1cc14c287fdb5c6d32233b29643b88c489092357 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 8 Jun 2026 11:22:35 -0700 Subject: [PATCH 10/15] Update example cfg file for AIS Fix typo in output path --- .../ismip7_forcing/ismip7_forcing_test.cfg | 32 ++++++++++--------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg index ff191d4c35..81728f11f4 100644 --- a/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg @@ -1,14 +1,17 @@ -# config options for ismip7 antarctic ice sheet forcing data -[ismip7_ais] +# config options for ismip7 forcing data +[ismip7] + +# Ice sheet: ais (Antarctic) or gis (Greenland) +ice_sheet = ais # Base path to the input ISMIP7 forcing files. User has to supply. -base_path_ismip7 = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/forcing/ +base_path_ismip7 = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/forcing/AIS/CESM2-WACCM/ssp585 # Base path to the MALI mesh. User has to supply. -base_path_mali = /global/cfs/cdirs/fanssie/MALI_input_files/AIS_4to20km_r01/ +base_path_mali = /global/cfs/cdirs/fanssie/MALI_input_files/AIS_4to20km_r01 # Base path to which output forcing files are saved. -output_base_path = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/test_processing +output_base_path = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/test_processing/AIS # Name of climate model used to generate ISMIP7 forcing data. User has to supply. # Available model names: CESM2-WACCM, MRI-ESM2-0 @@ -16,7 +19,7 @@ model = CESM2-WACCM # Scenario for forcing data. User has to supply. # Available scenarios: historical, ssp126, ssp370, ssp585 -scenario = historical +scenario = ssp585 # Name of the MALI mesh. Used to name mapping and output files. mali_mesh_name = AIS_4to20km_r01_20220907 @@ -25,29 +28,28 @@ mali_mesh_name = AIS_4to20km_r01_20220907 mali_mesh_file = AIS_4to20km_r01_20220907.nc # Number of MPI tasks for ESMF_RegridWeightGen -# Use 512 for 2km data sets. esmf_ntasks = 512 -# config options for ismip7 AIS atmosphere forcing -[ismip7_ais_atmosphere] +# config options for ismip7 atmosphere forcing +[ismip7_atmosphere] # Remapping method. Options: bilinear, neareststod, conserve method_remap = conserve # Start year for processing -start_year = 1980 +start_year = 2015 # End year for processing -end_year = 2014 +end_year = 2301 -# config options for ismip7 AIS ocean thermal forcing -[ismip7_ais_ocean_thermal] +# 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 = 1850 +start_year = 2015 # End year for processing -end_year = 2014 +end_year = 2301 From 5560e872da25b4d4a3772763967c98c67368b914 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 8 Jun 2026 15:05:16 -0700 Subject: [PATCH 11/15] Remove fields and attributes leftover from ncremap --- .../ismip7_forcing/atmosphere/process_runoff.py | 5 +++++ .../tests/ismip7_forcing/atmosphere/process_smb.py | 5 +++++ .../atmosphere/process_smb_gradient.py | 5 +++++ .../ismip7_forcing/atmosphere/process_temperature.py | 5 +++++ .../atmosphere/process_temperature_gradient.py | 5 +++++ .../ocean_thermal/process_thermal_forcing.py | 12 ++++++++++++ 6 files changed, 37 insertions(+) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py index 2fec1a9fd0..400e74f660 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -216,6 +216,11 @@ def _combine_and_rename(self, remapped_files, output_file): if vars_to_drop: ds = ds.drop_vars(vars_to_drop) + # Drop Time coordinate values (keep as dimension only); + # MALI uses xtime, not CF-encoded time coordinates + if "Time" in ds.coords: + ds = ds.drop_vars("Time") + write_netcdf(ds, output_file) def _extrapolate_source(self, input_file, output_file, varname, logger): diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py index 02b96e95c2..7cb250dc67 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -204,5 +204,10 @@ def _combine_and_rename(self, remapped_files, output_file): if vars_to_drop: ds = ds.drop_vars(vars_to_drop) + # Drop Time coordinate values (keep as dimension only); + # MALI uses xtime, not CF-encoded time coordinates + if "Time" in ds.coords: + ds = ds.drop_vars("Time") + write_netcdf(ds, output_file) ds.close() 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 bfb899431d..2f3e4c2cb4 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -205,5 +205,10 @@ def _combine_and_rename(self, remapped_files, output_file): if vars_to_drop: ds = ds.drop_vars(vars_to_drop) + # Drop Time coordinate values (keep as dimension only); + # MALI uses xtime, not CF-encoded time coordinates + if "Time" in ds.coords: + ds = ds.drop_vars("Time") + write_netcdf(ds, output_file) ds.close() diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index 2556b88227..1b44e8ab5f 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -204,5 +204,10 @@ def _combine_and_rename(self, remapped_files, output_file): if vars_to_drop: ds = ds.drop_vars(vars_to_drop) + # Drop Time coordinate values (keep as dimension only); + # MALI uses xtime, not CF-encoded time coordinates + if "Time" in ds.coords: + ds = ds.drop_vars("Time") + write_netcdf(ds, output_file) ds.close() 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 2730252423..5dec7b78bc 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -206,4 +206,9 @@ def _combine_and_rename(self, remapped_files, output_file): if vars_to_drop: ds = ds.drop_vars(vars_to_drop) + # Drop Time coordinate values (keep as dimension only); + # MALI uses xtime, not CF-encoded time coordinates + if "Time" in ds.coords: + ds = ds.drop_vars("Time") + write_netcdf(ds, output_file) 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 2b24ebdba3..5ed48361ff 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 @@ -241,6 +241,10 @@ def _combine_and_rename_3d(self, remapped_files, output_file, ds["ismip6shelfMelt_3dThermalForcing"].transpose( "Time", "nCells", "nISMIP6OceanLayers") + # Ensure double precision for MALI compatibility + ds["ismip6shelfMelt_3dThermalForcing"] = \ + ds["ismip6shelfMelt_3dThermalForcing"].astype(float) + # Add xtime variable with annual timestamps xtime = [] for t_index in range(ds.sizes["Time"]): @@ -258,6 +262,8 @@ def _combine_and_rename_3d(self, remapped_files, output_file, "melting method", "units": "degC", } + # Remove stale encoding (e.g. 'coordinates' from ncremap) + ds["ismip6shelfMelt_3dThermalForcing"].encoding.clear() ds["ismip6shelfMelt_zOcean"].attrs = { "long_name": "depth coordinate for ocean thermal forcing", "units": "m", @@ -327,6 +333,10 @@ def _combine_and_rename_2d(self, remapped_files, output_file, if "tf" in ds: ds = ds.rename({"tf": "ismip6_2dThermalForcing"}) + # Ensure double precision for MALI compatibility + ds["ismip6_2dThermalForcing"] = \ + ds["ismip6_2dThermalForcing"].astype(float) + # Add xtime variable with monthly timestamps xtime = [] for t_index in range(ds.sizes["Time"]): @@ -345,6 +355,8 @@ def _combine_and_rename_2d(self, remapped_files, output_file, "melting parameterization", "units": "degC", } + # Remove stale encoding (e.g. 'coordinates' from ncremap) + ds["ismip6_2dThermalForcing"].encoding.clear() # Drop auxiliary variables from remapping vars_to_drop = [v for v in ["lon", "lon_vertices", "lat", From 298d513ca671bb8807e8dce7b1f195109f7c2d54 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 8 Jun 2026 15:40:25 -0700 Subject: [PATCH 12/15] Process observational thermal forcing data set --- .../tests/ismip7_forcing/ismip7_forcing.cfg | 16 ++ .../ismip7_forcing/ismip7_forcing_test.cfg | 16 ++ .../ocean_thermal/process_thermal_forcing.py | 186 ++++++++++++++++++ 3 files changed, 218 insertions(+) diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg index e7bf3c4e6a..b996d2aaba 100644 --- a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg @@ -30,6 +30,12 @@ mali_mesh_file = NotAvailable # Number of MPI tasks for ESMF_RegridWeightGen esmf_ntasks = 128 +# 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 = true + # config options for ismip7 atmosphere forcing [ismip7_atmosphere] @@ -53,3 +59,13 @@ start_year = 1850 # End year for processing end_year = 2014 + +# 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 +# (directory containing tf/, so/, thetao/ subdirs) +base_path_climatology = /path/to/ISMIP7/forcing/AIS/obs/zhou_annual_06_nov diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg index 81728f11f4..270f5d396c 100644 --- a/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg @@ -30,6 +30,12 @@ mali_mesh_file = AIS_4to20km_r01_20220907.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 = false + +# Whether to process observational ocean thermal forcing climatology +process_ocean_climatology = true + # config options for ismip7 atmosphere forcing [ismip7_atmosphere] @@ -53,3 +59,13 @@ start_year = 2015 # End year for processing end_year = 2301 + +# 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 +# (directory containing tf/, so/, thetao/ subdirs) +base_path_climatology = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/forcing/AIS/obs/zhou_annual_06_nov 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 5ed48361ff..c2bfc3f033 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 @@ -54,6 +54,22 @@ def run(self): """ Run this step of the test case """ + config = self.config + section = config["ismip7"] + + # Check if we should process climatology data + if section.getboolean("process_ocean_climatology"): + self._run_climatology() + + # Check if we should process scenario (time-varying) data + if section.getboolean("process_ocean_thermal"): + self._run_scenario() + + def _run_scenario(self): + """ + Process time-varying ocean thermal forcing from an ESM + (e.g., CESM2-WACCM historical or ssp585). + """ logger = self.logger config = self.config params = get_params(config) @@ -178,6 +194,95 @@ def run(self): logger.info(f"Done. Output: {dst}") + def _run_climatology(self): + """ + Process observational ocean thermal forcing climatology + (e.g., Zhou et al. for AIS). This is a static 3D field with + no time dimension. + """ + logger = self.logger + config = self.config + + section = config["ismip7"] + mali_mesh_name = section.get("mali_mesh_name") + mali_mesh_file = section.get("mali_mesh_file") + output_base_path = section.get("output_base_path") + ice_sheet = section.get("ice_sheet") + + section = config["ismip7_ocean_climatology"] + method_remap = section.get("method_remap") + base_path_climatology = section.get("base_path_climatology") + version = 'v3' + + # Discover climatology TF file + input_path = os.path.join(base_path_climatology, "tf", version) + all_files = sorted(glob.glob(os.path.join(input_path, "tf_*.nc"))) + + if not all_files: + raise FileNotFoundError( + f"No ocean climatology TF files found in:\n" + f" {input_path}") + + # Use the first (and likely only) file + input_file = all_files[0] + logger.info(f"Processing ocean TF climatology: " + f"{os.path.basename(input_file)}") + + # Build mapping file using the climatology file as grid template. + mapping_file = (f"map_ismip7_{ice_sheet}_ocean_to_" + f"{mali_mesh_name}_{method_remap}.nc") + + if not os.path.exists(mapping_file): + logger.info("Building mapping file for ocean grid...") + build_mapping_file(config, logger, + input_file, mapping_file, + mali_mesh_file=mali_mesh_file, + method_remap=method_remap) + + # Extrapolate and remap + basename = os.path.basename(input_file) + remapped_file = f"remapped_{basename}" + + 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) + + logger.info(f" Remapping: {basename}") + args = ["ncremap", + "-i", extrap_file, + "-o", remapped_file, + "-m", mapping_file, + "-v", "tf"] + + check_call(args, logger=logger) + + # Clean up extrapolated source file + os.remove(extrap_file) + + # Rename to MALI conventions + logger.info("Renaming variables to MALI conventions...") + output_file = (f"{mali_mesh_name}_thermal_forcing_climatology_" + f"{version}.nc") + + self._rename_climatology_3d(remapped_file, output_file) + + # Clean up remapped file + if os.path.exists(remapped_file): + os.remove(remapped_file) + + # Place output in appropriate directory + output_path = os.path.join(output_base_path, "ocean_thermal_forcing", + "climatology") + 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 _combine_and_rename_3d(self, remapped_files, output_file, start_year, end_year): """ @@ -372,6 +477,87 @@ def _combine_and_rename_2d(self, remapped_files, output_file, write_netcdf(ds, output_file) + def _rename_climatology_3d(self, remapped_file, output_file): + """ + Rename dimensions and variables in a remapped 3D climatology + file (no time dimension) to MALI conventions. + + Parameters + ---------- + remapped_file : str + Path to the remapped NetCDF file + + output_file : str + Output file path + """ + ds = xr.open_dataset(remapped_file, engine="netcdf4") + + # Extract z coordinate and bounds before renaming + z_ocean = ds["z"] + z_bnds = ds["z_bnds"] + + # Rename dimensions to MALI conventions + rename_dims = {} + if "ncol" in ds.dims: + rename_dims["ncol"] = "nCells" + if "z" in ds.dims: + rename_dims["z"] = "nISMIP6OceanLayers" + if "bnds" in ds.dims: + rename_dims["bnds"] = "TWO" + if rename_dims: + ds = ds.rename(rename_dims) + + # Rename thermal forcing variable + if "tf" in ds: + ds = ds.rename({"tf": "ismip6shelfMelt_3dThermalForcing"}) + + # Set z coordinate and bounds as MALI-named variables + ds["ismip6shelfMelt_zOcean"] = ( + "nISMIP6OceanLayers", z_ocean.values) + ds["ismip6shelfMelt_zBndsOcean"] = ( + ("TWO", "nISMIP6OceanLayers"), z_bnds.values.T) + + # Transpose thermal forcing to MALI dimension order + # NetCDF (C order): nCells, nISMIP6OceanLayers + ds["ismip6shelfMelt_3dThermalForcing"] = \ + ds["ismip6shelfMelt_3dThermalForcing"].transpose( + "nCells", "nISMIP6OceanLayers") + + # Ensure double precision for MALI compatibility + ds["ismip6shelfMelt_3dThermalForcing"] = \ + ds["ismip6shelfMelt_3dThermalForcing"].astype(float) + + # Set attributes + ds["ismip6shelfMelt_3dThermalForcing"].attrs = { + "long_name": "thermal forcing for ISMIP6 ice-shelf " + "melting method", + "units": "degC", + } + ds["ismip6shelfMelt_3dThermalForcing"].encoding.clear() + ds["ismip6shelfMelt_zOcean"].attrs = { + "long_name": "depth coordinate for ocean thermal forcing", + "units": "m", + } + ds["ismip6shelfMelt_zBndsOcean"].attrs = { + "long_name": "bounds for ISMIP6 ocean layers", + "units": "m", + } + + # Drop auxiliary variables from remapping + vars_to_drop = [v for v in ["lon", "lon_vertices", "lat", + "lat_vertices", "lon_bnds", "lat_bnds", + "area", "z_bnds", "time_bnds", + "x_bnds", "y_bnds"] + if v in ds] + if vars_to_drop: + ds = ds.drop_vars(vars_to_drop) + + # Drop the z coordinate if it persists + if "nISMIP6OceanLayers" in ds.coords: + 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 From 0f07ad42a6ebfadd7a2027b0213e805660a2a762 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 8 Jun 2026 15:46:24 -0700 Subject: [PATCH 13/15] Update docs to include ocean climatology --- .../landice/test_groups/ismip7_forcing.rst | 18 ++++++- .../landice/test_groups/ismip7_forcing.rst | 53 +++++++++++++++++-- 2 files changed, 65 insertions(+), 6 deletions(-) diff --git a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst index 2f4349e804..72aa8c7411 100644 --- a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst @@ -87,7 +87,17 @@ test case processes the ISMIP7 ocean thermal forcing. It contains a single step, which handles both AIS (3D, decade-spanning files) and GrIS (2D, yearly files) by branching on the ``ocean_3d`` parameter from ``ice_sheet_params``. -For AIS, the step: +The ``run()`` method dispatches to two sub-methods based on the boolean config +options ``process_ocean_thermal`` and ``process_ocean_climatology`` in the +``[ismip7]`` section: + +* ``_run_scenario()``: Processes time-varying ESM scenario data (model + + scenario combination). Uses config from ``[ismip7_ocean_thermal]``. +* ``_run_climatology()``: Processes the static observational climatology + (Zhou et al., AIS only). Uses config from ``[ismip7_ocean_climatology]``. + The TF version (currently v3) is hard-coded. + +For AIS scenario data, the step: * Remaps thermal forcing preserving 30 vertical ocean layers * Produces ``ismip6shelfMelt_3dThermalForcing`` (dims: Time × nCells × @@ -95,6 +105,12 @@ For AIS, the step: * Includes depth coordinate variables ``ismip6shelfMelt_zOcean`` and ``ismip6shelfMelt_zBndsOcean`` +For AIS climatology data, the step: + +* Extrapolates fill values, remaps, and renames to MALI conventions +* Produces ``ismip6shelfMelt_3dThermalForcing`` (dims: nCells × + nISMIP6OceanLayers) — no Time dimension + For GrIS, the step: * Remaps 2D monthly thermal forcing diff --git a/docs/users_guide/landice/test_groups/ismip7_forcing.rst b/docs/users_guide/landice/test_groups/ismip7_forcing.rst index c797bf5dbf..1740f95a5f 100644 --- a/docs/users_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/users_guide/landice/test_groups/ismip7_forcing.rst @@ -18,7 +18,9 @@ The test group includes two test cases: ``atmosphere`` and ``ocean_thermal``. * 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. + 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. (For more details on the steps of each test case, see :ref:`landice_ismip7_forcing_atmosphere` and @@ -65,6 +67,12 @@ For AIS ocean thermal (8km, 30 depth levels, decade files): ocean/tf/v3/tf_AIS_{model}_{scenario}_ocean_v3_{start_year}-{end_year}.nc +For AIS ocean thermal climatology (8km, 30 depth levels, static): + +.. code-block:: none + + {base_path_climatology}/tf/v3/tf_AIS_obs_ocean_climatology_*.nc + For GrIS atmosphere (1km, polar stereographic EPSG:3413): .. code-block:: none @@ -86,7 +94,7 @@ For GrIS ocean thermal (same 1km grid, 2D, yearly files): config options -------------- -The ``ismip7_forcing`` test group uses three config sections. The default +The ``ismip7_forcing`` test group uses four config sections. The default values are: .. code-block:: cfg @@ -121,6 +129,12 @@ values are: # Number of MPI tasks for ESMF_RegridWeightGen esmf_ntasks = 128 + # 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 = true + # config options for ismip7 atmosphere forcing [ismip7_atmosphere] @@ -145,9 +159,23 @@ values are: # End year for processing end_year = 2014 + # config options for ismip7 ocean thermal forcing climatology + [ismip7_ocean_climatology] + + # Remapping method: bilinear, neareststod, conserve + method_remap = bilinear + + # Base path to observational climatology data + base_path_climatology = /path/to/ISMIP7/forcing/AIS/obs/zhou_annual_06_nov + All ``NotAvailable`` options must be overridden in a user config file passed at setup time (e.g., ``compass setup ... -f my_ismip7.cfg``). +The boolean options ``process_ocean_thermal`` and ``process_ocean_climatology`` +control which processing paths are executed when the ``ocean_thermal`` test +case is run. Both default to ``true``. Set one to ``false`` in your user +config to skip that processing path. + .. _landice_ismip7_forcing_atmosphere: atmosphere @@ -185,13 +213,28 @@ The ``landice/ismip7_forcing/ocean_thermal`` test case processes the ISMIP7 ocean thermal forcing (``tf``) and remaps it from the native polar stereographic grid to the MALI unstructured mesh. -For **AIS**, thermal forcing is 3D with 30 vertical ocean layers. The input -files span decades (e.g., 1850-1859). The output variable is -``ismip6shelfMelt_3dThermalForcing`` with dimension +The step supports two processing modes, controlled by boolean config options +in the ``[ismip7]`` section: + +* **Scenario (time-varying) data** (``process_ocean_thermal = true``): + Processes ESM-driven thermal forcing for a given model/scenario combination. + +* **Observational climatology** (``process_ocean_climatology = true``): + Processes the static Zhou et al. observational thermal forcing climatology + (AIS only). This is a time-invariant 3D field referenced to 1995-2024. + +Both modes can be enabled simultaneously. + +For **AIS** scenario data, thermal forcing is 3D with 30 vertical ocean +layers. The input files span decades (e.g., 1850-1859). The output variable +is ``ismip6shelfMelt_3dThermalForcing`` with dimension ``nISMIP6OceanLayers``. Associated depth coordinate variables ``ismip6shelfMelt_zOcean`` and ``ismip6shelfMelt_zBndsOcean`` are also produced. +For **AIS** climatology data, the output is the same 3D thermal forcing field +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``. From ef09feb2ab1a68cd8dc87d742238dc30bb11996f Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 8 Jun 2026 20:13:10 -0700 Subject: [PATCH 14/15] Add Time dimension to climatology thermal forcing --- .../ismip7_forcing/ocean_thermal/process_thermal_forcing.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) 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 c2bfc3f033..28a7345e60 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 @@ -511,6 +511,9 @@ def _rename_climatology_3d(self, remapped_file, output_file): if "tf" in ds: ds = ds.rename({"tf": "ismip6shelfMelt_3dThermalForcing"}) + ds["ismip6shelfMelt_3dThermalForcing"] = \ + ds["ismip6shelfMelt_3dThermalForcing"].expand_dims("Time", axis=0) + # Set z coordinate and bounds as MALI-named variables ds["ismip6shelfMelt_zOcean"] = ( "nISMIP6OceanLayers", z_ocean.values) @@ -521,7 +524,7 @@ def _rename_climatology_3d(self, remapped_file, output_file): # NetCDF (C order): nCells, nISMIP6OceanLayers ds["ismip6shelfMelt_3dThermalForcing"] = \ ds["ismip6shelfMelt_3dThermalForcing"].transpose( - "nCells", "nISMIP6OceanLayers") + "Time", "nCells", "nISMIP6OceanLayers") # Ensure double precision for MALI compatibility ds["ismip6shelfMelt_3dThermalForcing"] = \ From 312a04e0fdc75b65dc64bd6dc89b4aec9bb67f6c Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 11 Jun 2026 08:50:28 -0700 Subject: [PATCH 15/15] Reinterpret time convention for monthly forcing files Interpret the mid-month time stamp as representing the midpoint of the forcing interval, rather than the start date. For instance, time = 14 days since 1850-01-01 should be interpreted as representing Jan 1 to Jan 31, 1850, rather than Jan 15 to Feb 14. --- .../tests/ismip7_forcing/atmosphere/process_runoff.py | 5 ++++- .../landice/tests/ismip7_forcing/atmosphere/process_smb.py | 5 ++++- .../tests/ismip7_forcing/atmosphere/process_temperature.py | 5 ++++- .../ismip7_forcing/ocean_thermal/process_thermal_forcing.py | 5 ++++- 4 files changed, 16 insertions(+), 4 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py index 400e74f660..e03e85ecd6 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -192,12 +192,15 @@ def _combine_and_rename(self, remapped_files, output_file): ds = ds.rename({"mrro": "ismip6Runoff"}) # Add xtime variable with monthly timestamps + # ISMIP7 files encode time at mid-month (e.g., Jan 15) but + # this represents forcing for the full month (Jan 1-31). + # MALI needs xtime at the start of each forcing interval. xtime = [] for t_index in range(ds.sizes["Time"]): date = ds.Time[t_index] yr = int(date.dt.year.values) mo = int(date.dt.month.values) - date_str = f"{yr:04d}-{mo:02d}-15_00:00:00".ljust(64) + date_str = f"{yr:04d}-{mo:02d}-01_00:00:00".ljust(64) xtime.append(date_str) ds["xtime"] = ("Time", xtime) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py index 7cb250dc67..222fb6e92e 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -180,12 +180,15 @@ def _combine_and_rename(self, remapped_files, output_file): ds = ds.rename({"acabf": "sfcMassBal"}) # Add xtime variable with monthly timestamps + # ISMIP7 files encode time at mid-month (e.g., Jan 15) but + # this represents forcing for the full month (Jan 1-31). + # MALI needs xtime at the start of each forcing interval. xtime = [] for t_index in range(ds.sizes["Time"]): date = ds.Time[t_index] yr = int(date.dt.year.values) mo = int(date.dt.month.values) - date_str = f"{yr:04d}-{mo:02d}-15_00:00:00".ljust(64) + date_str = f"{yr:04d}-{mo:02d}-01_00:00:00".ljust(64) xtime.append(date_str) ds["xtime"] = ("Time", xtime) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index 1b44e8ab5f..caa8f6be9f 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -180,12 +180,15 @@ def _combine_and_rename(self, remapped_files, output_file): ds = ds.rename({"ts": "surfaceAirTemperature"}) # Add xtime variable with monthly timestamps + # ISMIP7 files encode time at mid-month (e.g., Jan 15) but + # this represents forcing for the full month (Jan 1-31). + # MALI needs xtime at the start of each forcing interval. xtime = [] for t_index in range(ds.sizes["Time"]): date = ds.Time[t_index] yr = int(date.dt.year.values) mo = int(date.dt.month.values) - date_str = f"{yr:04d}-{mo:02d}-15_00:00:00".ljust(64) + date_str = f"{yr:04d}-{mo:02d}-01_00:00:00".ljust(64) xtime.append(date_str) ds["xtime"] = ("Time", xtime) 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 28a7345e60..a1469c8f2d 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 @@ -443,12 +443,15 @@ def _combine_and_rename_2d(self, remapped_files, output_file, ds["ismip6_2dThermalForcing"].astype(float) # Add xtime variable with monthly timestamps + # ISMIP7 files encode time at mid-month (e.g., Jan 15) but + # this represents forcing for the full month (Jan 1-31). + # MALI needs xtime at the start of each forcing interval. xtime = [] for t_index in range(ds.sizes["Time"]): date = ds.Time[t_index] yr = int(date.dt.year.values) mo = int(date.dt.month.values) - date_str = f"{yr:04d}-{mo:02d}-15_00:00:00".ljust(64) + date_str = f"{yr:04d}-{mo:02d}-01_00:00:00".ljust(64) xtime.append(date_str) ds["xtime"] = ("Time", xtime)