From a81e95a38c056d38e77c5009918809111b0fb27a Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 3 Jun 2026 16:26:34 -0700 Subject: [PATCH 1/9] 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 2/9] 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 3/9] 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 4/9] 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 5/9] 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 6/9] 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 7/9] 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 8/9] 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 9/9] 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: