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..3f621a4c37 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/__init__.py @@ -0,0 +1,24 @@ +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 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)) + self.add_test_case(OceanThermal(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..7e577b2026 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/__init__.py @@ -0,0 +1,53 @@ +from compass.landice.tests.ismip7_forcing.atmosphere.process_runoff import ( + ProcessRunoff, +) +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.atmosphere.process_temperature_gradient import ( # noqa: E501 + ProcessTemperatureGradient, +) +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 atmosphere forcing data. + Remaps monthly SMB, temperature, and annual gradients (SMB and + temperature) from the ISMIP7 polar stereographic grid to the + MALI unstructured mesh. Also processes runoff (mrro). + """ + + 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)) + 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) 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..2fec1a9fd0 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -0,0 +1,272 @@ +import glob +import os +import shutil + +import numpy as np +import xarray as xr +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call +from scipy.ndimage import distance_transform_edt + +from compass.landice.tests.ismip7_forcing.create_mapfile import ( + build_mapping_file, +) +from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.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 + + # 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", 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}_" + 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) + + def _extrapolate_source(self, input_file, output_file, varname, logger): + """ + Extrapolate fill/missing values on the source polar stereographic + grid using nearest-neighbor via distance_transform_edt. This must + be done before remapping so that fill values don't contaminate the + interpolation stencil. + + Parameters + ---------- + input_file : str + Path to the input NetCDF file on the source grid + + output_file : str + Path to write the extrapolated file + + varname : str + Name of the variable to extrapolate (e.g., "mrro") + + logger : logging.Logger + Logger for status messages + """ + logger.info(f" Extrapolating fill values on source grid: " + f"{os.path.basename(input_file)}") + + ds = xr.open_dataset(input_file, engine="netcdf4") + data = ds[varname] + + # Process each time step + # Source files have dims like (time, y, x) + values = data.values.copy() + non_spatial_shape = values.shape[:-2] # (time,) + + for idx in np.ndindex(non_spatial_shape): + slab = values[idx] # shape (ny, nx) + valid_mask = np.isfinite(slab) + if valid_mask.all() or not valid_mask.any(): + continue + nearest_inds = distance_transform_edt( + ~valid_mask, return_distances=False, return_indices=True) + invalid = ~valid_mask + values[idx][invalid] = slab[ + nearest_inds[0, invalid], + nearest_inds[1, invalid]] + + ds[varname] = (data.dims, values) + ds[varname].attrs = data.attrs + + # Remove _FillValue encoding so output has no masked values + if "_FillValue" in ds[varname].encoding: + del ds[varname].encoding["_FillValue"] + + write_netcdf(ds, output_file) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py new file mode 100644 index 0000000000..02b96e95c2 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -0,0 +1,208 @@ +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 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") + + 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") + + 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, "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: + 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 + 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...") + 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", "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..bfb899431d --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -0,0 +1,209 @@ +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 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") + + 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") + + 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, "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: + 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) + 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...") + 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", "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) + + # Rename to MALI convention (PR #169) + if "dacabfdz" in ds: + ds = ds.rename({"dacabfdz": "sfcMassBalLapseRate"}) + + # 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["sfcMassBalLapseRate"].attrs = { + "long_name": "vertical gradient dSMBdz used for SMB " + "elevation-change correction", + "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..2556b88227 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -0,0 +1,208 @@ +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 ProcessTemperature(Step): + """ + 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. + """ + + 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") + + 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") + + 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, "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: + 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) + 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...") + 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", "ts"] + + 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 "ts" in ds: + ds = ds.rename({"ts": "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": "temperature at top of ice sheet model", + "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/atmosphere/process_temperature_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py new file mode 100644 index 0000000000..2730252423 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -0,0 +1,209 @@ +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 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") + + 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") + + 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, "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: + 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) + 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...") + 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", "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) + + # Rename to MALI convention (PR #169) + if "dtsdz" in ds: + ds = ds.rename({"dtsdz": "surfaceAirTemperatureLapseRate"}) + + # 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["surfaceAirTemperatureLapseRate"].attrs = { + "long_name": "vertical gradient dTdz used for SAT " + "elevation-change correction", + "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) diff --git a/compass/landice/tests/ismip7_forcing/configure.py b/compass/landice/tests/ismip7_forcing/configure.py new file mode 100644 index 0000000000..8719c0e141 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/configure.py @@ -0,0 +1,22 @@ +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" + 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) + 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..de10f32b34 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/create_mapfile.py @@ -0,0 +1,114 @@ +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, logger, ismip7_grid_file, + mapping_file, mali_mesh_file=None, + method_remap=None, projection=None): + """ + Build a mapping file for regridding from an ISMIP7 polar + stereographic grid to the MALI unstructured mesh. + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for the test case + + 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' + + 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): + 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'.") + + # 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" + 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}") + + section = config["ismip7"] + cores = section.getint("esmf_ntasks") + + 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/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 new file mode 100644 index 0000000000..e7bf3c4e6a --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg @@ -0,0 +1,55 @@ +# 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 + +# 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 + +# Number of MPI tasks for ESMF_RegridWeightGen +esmf_ntasks = 128 + +# config options for ismip7 atmosphere forcing +[ismip7_atmosphere] + +# Remapping method. Options: 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. 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..2a97700428 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/__init__.py @@ -0,0 +1,38 @@ +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 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): + """ + 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..2b24ebdba3 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -0,0 +1,417 @@ +import glob +import os +import shutil + +import numpy as np +import xarray as xr +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call +from scipy.ndimage import distance_transform_edt + +from compass.landice.tests.ismip7_forcing.create_mapfile import ( + build_mapping_file, +) +from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.step import Step + + +class ProcessThermalForcing(Step): + """ + A step for processing ISMIP7 ocean thermal forcing (tf) data. + 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): + """ + 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"] + 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_ocean_thermal"] + method_remap = section.get("method_remap") + start_year = section.getint("start_year") + end_year = section.getint("end_year") + + # 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: + 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. + # 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]) # same as start for single-year files + 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. + mapping_file = (f"map_ismip7_{ice_sheet}_ocean_to_" + f"{mali_mesh_name}_{method_remap}.nc") + + if not os.path.exists(mapping_file): + logger.info("Building mapping file for ocean grid...") + build_mapping_file(config, logger, + input_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 + + # 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", 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}_" + f"{start_year}-{end_year}.nc") + + 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...") + 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_3d(self, remapped_files, output_file, + start_year, end_year): + """ + Combine decade-spanning remapped files (AIS), subset to the + requested year range, and rename variables/dimensions to MALI + conventions for 3D 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)) + + # 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) + + 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) + + def _extrapolate_source(self, input_file, output_file, varname, logger): + """ + Extrapolate fill/missing values on the source polar stereographic + grid using nearest-neighbor interpolation from valid cells. This + must be done before remapping so that fill values don't contaminate + the interpolation stencil. + + Parameters + ---------- + input_file : str + Path to the input NetCDF file on the source grid + + output_file : str + Path to write the extrapolated file + + varname : str + Name of the variable to extrapolate (e.g., "tf") + + logger : logging.Logger + Logger for status messages + """ + logger.info(f" Extrapolating fill values on source grid: " + f"{os.path.basename(input_file)}") + + ds = xr.open_dataset(input_file, engine="netcdf4") + data = ds[varname] + + # Process each time step (and z level if 3D) + # Source files have dims like (time, z, y, x) or (time, y, x) + values = data.values.copy() + non_spatial_shape = values.shape[:-2] # (time,) or (time, z) + + # Use distance_transform_edt with return_indices to find the + # nearest valid cell index for each invalid cell. This is O(n) + # on the grid and much faster than KD-tree approaches. + for idx in np.ndindex(non_spatial_shape): + slab = values[idx] # shape (ny, nx) + valid_mask = np.isfinite(slab) + if valid_mask.all() or not valid_mask.any(): + continue + nearest_inds = distance_transform_edt( + ~valid_mask, return_distances=False, return_indices=True) + invalid = ~valid_mask + values[idx][invalid] = slab[ + nearest_inds[0, invalid], + nearest_inds[1, invalid]] + + ds[varname] = (data.dims, values) + ds[varname].attrs = data.attrs + + # Remove _FillValue encoding so output has no masked values + if "_FillValue" in ds[varname].encoding: + del ds[varname].encoding["_FillValue"] + + write_netcdf(ds, output_file) diff --git a/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..2f4349e804 --- /dev/null +++ b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst @@ -0,0 +1,101 @@ +.. _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 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: + +* :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`` + +.. _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..c797bf5dbf --- /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 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 + 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 + mrro/v2/mrro_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**: 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``.