diff --git a/compass/landice/__init__.py b/compass/landice/__init__.py index 487f1e8345..aafe35fdd6 100644 --- a/compass/landice/__init__.py +++ b/compass/landice/__init__.py @@ -11,6 +11,8 @@ 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.ismip7_run import Ismip7Run 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 +48,8 @@ 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(Ismip7Run(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..e03e85ecd6 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -0,0 +1,280 @@ +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 + # ISMIP7 files encode time at mid-month (e.g., Jan 15) but + # this represents forcing for the full month (Jan 1-31). + # MALI needs xtime at the start of each forcing interval. + xtime = [] + for t_index in range(ds.sizes["Time"]): + date = ds.Time[t_index] + yr = int(date.dt.year.values) + mo = int(date.dt.month.values) + date_str = f"{yr:04d}-{mo:02d}-01_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) + + # Drop Time coordinate values (keep as dimension only); + # MALI uses xtime, not CF-encoded time coordinates + if "Time" in ds.coords: + ds = ds.drop_vars("Time") + + write_netcdf(ds, output_file) + + def _extrapolate_source(self, input_file, output_file, varname, logger): + """ + 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..222fb6e92e --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -0,0 +1,216 @@ +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 + # ISMIP7 files encode time at mid-month (e.g., Jan 15) but + # this represents forcing for the full month (Jan 1-31). + # MALI needs xtime at the start of each forcing interval. + xtime = [] + for t_index in range(ds.sizes["Time"]): + date = ds.Time[t_index] + yr = int(date.dt.year.values) + mo = int(date.dt.month.values) + date_str = f"{yr:04d}-{mo:02d}-01_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) + + # Drop Time coordinate values (keep as dimension only); + # MALI uses xtime, not CF-encoded time coordinates + if "Time" in ds.coords: + ds = ds.drop_vars("Time") + + write_netcdf(ds, output_file) + ds.close() diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py new file mode 100644 index 0000000000..2f3e4c2cb4 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -0,0 +1,214 @@ +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) + + # Drop Time coordinate values (keep as dimension only); + # MALI uses xtime, not CF-encoded time coordinates + if "Time" in ds.coords: + ds = ds.drop_vars("Time") + + write_netcdf(ds, output_file) + ds.close() diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py new file mode 100644 index 0000000000..caa8f6be9f --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -0,0 +1,216 @@ +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 + # ISMIP7 files encode time at mid-month (e.g., Jan 15) but + # this represents forcing for the full month (Jan 1-31). + # MALI needs xtime at the start of each forcing interval. + xtime = [] + for t_index in range(ds.sizes["Time"]): + date = ds.Time[t_index] + yr = int(date.dt.year.values) + mo = int(date.dt.month.values) + date_str = f"{yr:04d}-{mo:02d}-01_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) + + # Drop Time coordinate values (keep as dimension only); + # MALI uses xtime, not CF-encoded time coordinates + if "Time" in ds.coords: + ds = ds.drop_vars("Time") + + write_netcdf(ds, output_file) + ds.close() diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py new file mode 100644 index 0000000000..5dec7b78bc --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -0,0 +1,214 @@ +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) + + # Drop Time coordinate values (keep as dimension only); + # MALI uses xtime, not CF-encoded time coordinates + if "Time" in ds.coords: + ds = ds.drop_vars("Time") + + write_netcdf(ds, output_file) diff --git a/compass/landice/tests/ismip7_forcing/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..b996d2aaba --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg @@ -0,0 +1,71 @@ +# 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 + +# Whether to process time-varying ocean thermal forcing (ESM scenario data) +process_ocean_thermal = true + +# Whether to process observational ocean thermal forcing climatology +process_ocean_climatology = true + +# config options for ismip7 atmosphere forcing +[ismip7_atmosphere] + +# 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 + +# config options for ismip7 ocean thermal forcing climatology +[ismip7_ocean_climatology] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = bilinear + +# Base path to observational climatology data +# (directory containing tf/, so/, thetao/ subdirs) +base_path_climatology = /path/to/ISMIP7/forcing/AIS/obs/zhou_annual_06_nov diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg new file mode 100644 index 0000000000..270f5d396c --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg @@ -0,0 +1,71 @@ +# config options for ismip7 forcing data +[ismip7] + +# Ice sheet: ais (Antarctic) or gis (Greenland) +ice_sheet = ais + +# Base path to the input ISMIP7 forcing files. User has to supply. +base_path_ismip7 = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/forcing/AIS/CESM2-WACCM/ssp585 + +# Base path to the MALI mesh. User has to supply. +base_path_mali = /global/cfs/cdirs/fanssie/MALI_input_files/AIS_4to20km_r01 + +# Base path to which output forcing files are saved. +output_base_path = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/test_processing/AIS + +# Name of climate model used to generate ISMIP7 forcing data. User has to supply. +# Available model names: CESM2-WACCM, MRI-ESM2-0 +model = CESM2-WACCM + +# Scenario for forcing data. User has to supply. +# Available scenarios: historical, ssp126, ssp370, ssp585 +scenario = ssp585 + +# 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 +esmf_ntasks = 512 + +# Whether to process time-varying ocean thermal forcing (ESM scenario data) +process_ocean_thermal = false + +# Whether to process observational ocean thermal forcing climatology +process_ocean_climatology = true + +# config options for ismip7 atmosphere forcing +[ismip7_atmosphere] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = conserve + +# Start year for processing +start_year = 2015 + +# End year for processing +end_year = 2301 + +# 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 = 2015 + +# End year for processing +end_year = 2301 + +# config options for ismip7 ocean thermal forcing climatology +[ismip7_ocean_climatology] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = bilinear + +# Base path to observational climatology data +# (directory containing tf/, so/, thetao/ subdirs) +base_path_climatology = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/forcing/AIS/obs/zhou_annual_06_nov diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/__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..a1469c8f2d --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -0,0 +1,621 @@ +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 + """ + config = self.config + section = config["ismip7"] + + # Check if we should process climatology data + if section.getboolean("process_ocean_climatology"): + self._run_climatology() + + # Check if we should process scenario (time-varying) data + if section.getboolean("process_ocean_thermal"): + self._run_scenario() + + def _run_scenario(self): + """ + Process time-varying ocean thermal forcing from an ESM + (e.g., CESM2-WACCM historical or ssp585). + """ + logger = self.logger + config = self.config + params = get_params(config) + + 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 _run_climatology(self): + """ + Process observational ocean thermal forcing climatology + (e.g., Zhou et al. for AIS). This is a static 3D field with + no time dimension. + """ + logger = self.logger + config = self.config + + section = config["ismip7"] + mali_mesh_name = section.get("mali_mesh_name") + mali_mesh_file = section.get("mali_mesh_file") + output_base_path = section.get("output_base_path") + ice_sheet = section.get("ice_sheet") + + section = config["ismip7_ocean_climatology"] + method_remap = section.get("method_remap") + base_path_climatology = section.get("base_path_climatology") + version = 'v3' + + # Discover climatology TF file + input_path = os.path.join(base_path_climatology, "tf", version) + all_files = sorted(glob.glob(os.path.join(input_path, "tf_*.nc"))) + + if not all_files: + raise FileNotFoundError( + f"No ocean climatology TF files found in:\n" + f" {input_path}") + + # Use the first (and likely only) file + input_file = all_files[0] + logger.info(f"Processing ocean TF climatology: " + f"{os.path.basename(input_file)}") + + # Build mapping file using the climatology file as grid template. + mapping_file = (f"map_ismip7_{ice_sheet}_ocean_to_" + f"{mali_mesh_name}_{method_remap}.nc") + + if not os.path.exists(mapping_file): + logger.info("Building mapping file for ocean grid...") + build_mapping_file(config, logger, + input_file, mapping_file, + mali_mesh_file=mali_mesh_file, + method_remap=method_remap) + + # Extrapolate and remap + basename = os.path.basename(input_file) + remapped_file = f"remapped_{basename}" + + if not os.path.exists(remapped_file): + extrap_file = f"extrap_{basename}" + if not os.path.exists(extrap_file): + self._extrapolate_source(input_file, extrap_file, "tf", + logger) + + logger.info(f" Remapping: {basename}") + args = ["ncremap", + "-i", extrap_file, + "-o", remapped_file, + "-m", mapping_file, + "-v", "tf"] + + check_call(args, logger=logger) + + # Clean up extrapolated source file + os.remove(extrap_file) + + # Rename to MALI conventions + logger.info("Renaming variables to MALI conventions...") + output_file = (f"{mali_mesh_name}_thermal_forcing_climatology_" + f"{version}.nc") + + self._rename_climatology_3d(remapped_file, output_file) + + # Clean up remapped file + if os.path.exists(remapped_file): + os.remove(remapped_file) + + # Place output in appropriate directory + output_path = os.path.join(output_base_path, "ocean_thermal_forcing", + "climatology") + if not os.path.exists(output_path): + os.makedirs(output_path) + + dst = os.path.join(output_path, output_file) + shutil.copy(output_file, dst) + + logger.info(f"Done. Output: {dst}") + + def _combine_and_rename_3d(self, remapped_files, output_file, + start_year, end_year): + """ + 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") + + # Ensure double precision for MALI compatibility + ds["ismip6shelfMelt_3dThermalForcing"] = \ + ds["ismip6shelfMelt_3dThermalForcing"].astype(float) + + # Add xtime variable with annual timestamps + xtime = [] + for t_index in range(ds.sizes["Time"]): + 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", + } + # Remove stale encoding (e.g. 'coordinates' from ncremap) + ds["ismip6shelfMelt_3dThermalForcing"].encoding.clear() + ds["ismip6shelfMelt_zOcean"].attrs = { + "long_name": "depth coordinate for ocean thermal forcing", + "units": "m", + } + 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"}) + + # Ensure double precision for MALI compatibility + ds["ismip6_2dThermalForcing"] = \ + ds["ismip6_2dThermalForcing"].astype(float) + + # Add xtime variable with monthly timestamps + # ISMIP7 files encode time at mid-month (e.g., Jan 15) but + # this represents forcing for the full month (Jan 1-31). + # MALI needs xtime at the start of each forcing interval. + xtime = [] + for t_index in range(ds.sizes["Time"]): + date = ds.Time[t_index] + yr = int(date.dt.year.values) + mo = int(date.dt.month.values) + date_str = f"{yr:04d}-{mo:02d}-01_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", + } + # Remove stale encoding (e.g. 'coordinates' from ncremap) + ds["ismip6_2dThermalForcing"].encoding.clear() + + # Drop auxiliary variables from remapping + vars_to_drop = [v for v in ["lon", "lon_vertices", "lat", + "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 _rename_climatology_3d(self, remapped_file, output_file): + """ + Rename dimensions and variables in a remapped 3D climatology + file (no time dimension) to MALI conventions. + + Parameters + ---------- + remapped_file : str + Path to the remapped NetCDF file + + output_file : str + Output file path + """ + ds = xr.open_dataset(remapped_file, engine="netcdf4") + + # Extract z coordinate and bounds before renaming + z_ocean = ds["z"] + z_bnds = ds["z_bnds"] + + # Rename dimensions to MALI conventions + rename_dims = {} + if "ncol" in ds.dims: + rename_dims["ncol"] = "nCells" + if "z" in ds.dims: + rename_dims["z"] = "nISMIP6OceanLayers" + if "bnds" in ds.dims: + rename_dims["bnds"] = "TWO" + if rename_dims: + ds = ds.rename(rename_dims) + + # Rename thermal forcing variable + if "tf" in ds: + ds = ds.rename({"tf": "ismip6shelfMelt_3dThermalForcing"}) + + ds["ismip6shelfMelt_3dThermalForcing"] = \ + ds["ismip6shelfMelt_3dThermalForcing"].expand_dims("Time", axis=0) + + # Set z coordinate and bounds as MALI-named variables + ds["ismip6shelfMelt_zOcean"] = ( + "nISMIP6OceanLayers", z_ocean.values) + ds["ismip6shelfMelt_zBndsOcean"] = ( + ("TWO", "nISMIP6OceanLayers"), z_bnds.values.T) + + # Transpose thermal forcing to MALI dimension order + # NetCDF (C order): nCells, nISMIP6OceanLayers + ds["ismip6shelfMelt_3dThermalForcing"] = \ + ds["ismip6shelfMelt_3dThermalForcing"].transpose( + "Time", "nCells", "nISMIP6OceanLayers") + + # Ensure double precision for MALI compatibility + ds["ismip6shelfMelt_3dThermalForcing"] = \ + ds["ismip6shelfMelt_3dThermalForcing"].astype(float) + + # Set attributes + ds["ismip6shelfMelt_3dThermalForcing"].attrs = { + "long_name": "thermal forcing for ISMIP6 ice-shelf " + "melting method", + "units": "degC", + } + ds["ismip6shelfMelt_3dThermalForcing"].encoding.clear() + ds["ismip6shelfMelt_zOcean"].attrs = { + "long_name": "depth coordinate for ocean thermal forcing", + "units": "m", + } + ds["ismip6shelfMelt_zBndsOcean"].attrs = { + "long_name": "bounds for ISMIP6 ocean layers", + "units": "m", + } + + # Drop auxiliary variables from remapping + vars_to_drop = [v for v in ["lon", "lon_vertices", "lat", + "lat_vertices", "lon_bnds", "lat_bnds", + "area", "z_bnds", "time_bnds", + "x_bnds", "y_bnds"] + if v in ds] + if vars_to_drop: + ds = ds.drop_vars(vars_to_drop) + + # Drop the z coordinate if it persists + if "nISMIP6OceanLayers" in ds.coords: + ds = ds.drop_vars("nISMIP6OceanLayers") + + write_netcdf(ds, output_file) + + def _extrapolate_source(self, input_file, output_file, varname, logger): + """ + Extrapolate fill/missing values on the source polar stereographic + grid using nearest-neighbor interpolation from valid cells. This + must be done before remapping so that fill values don't contaminate + the interpolation stencil. + + Parameters + ---------- + input_file : str + Path to the input NetCDF file on the source grid + + output_file : str + Path to write the extrapolated file + + varname : str + Name of the variable to extrapolate (e.g., "tf") + + logger : logging.Logger + Logger for status messages + """ + logger.info(f" Extrapolating fill values on source grid: " + f"{os.path.basename(input_file)}") + + ds = xr.open_dataset(input_file, engine="netcdf4") + data = ds[varname] + + # Process each time step (and z level if 3D) + # Source files have dims like (time, z, y, x) or (time, y, x) + values = data.values.copy() + non_spatial_shape = values.shape[:-2] # (time,) or (time, z) + + # Use distance_transform_edt with return_indices to find the + # nearest valid cell index for each invalid cell. This is O(n) + # on the grid and much faster than KD-tree approaches. + for idx in np.ndindex(non_spatial_shape): + slab = values[idx] # shape (ny, nx) + valid_mask = np.isfinite(slab) + if valid_mask.all() or not valid_mask.any(): + continue + nearest_inds = distance_transform_edt( + ~valid_mask, return_distances=False, return_indices=True) + invalid = ~valid_mask + values[idx][invalid] = slab[ + nearest_inds[0, invalid], + nearest_inds[1, invalid]] + + ds[varname] = (data.dims, values) + ds[varname].attrs = data.attrs + + # Remove _FillValue encoding so output has no masked values + if "_FillValue" in ds[varname].encoding: + del ds[varname].encoding["_FillValue"] + + write_netcdf(ds, output_file) diff --git a/compass/landice/tests/ismip7_run/__init__.py b/compass/landice/tests/ismip7_run/__init__.py new file mode 100644 index 0000000000..64eb375d00 --- /dev/null +++ b/compass/landice/tests/ismip7_run/__init__.py @@ -0,0 +1,19 @@ +from compass.landice.tests.ismip7_run.ismip7_ais import Ismip7Ais +from compass.landice.tests.ismip7_run.ismip7_gris import Ismip7Gris +from compass.testgroup import TestGroup + + +class Ismip7Run(TestGroup): + """ + A test group for automated setup of a suite of standardized + ISMIP7 simulations for both AIS and GrIS. + """ + def __init__(self, mpas_core): + """ + mpas_core : compass.landice.Landice + the MPAS core that this test group belongs to + """ + super().__init__(mpas_core=mpas_core, name='ismip7_run') + + self.add_test_case(Ismip7Ais(test_group=self)) + self.add_test_case(Ismip7Gris(test_group=self)) diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/__init__.py b/compass/landice/tests/ismip7_run/ismip7_ais/__init__.py new file mode 100644 index 0000000000..c6e3a41c94 --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_ais/__init__.py @@ -0,0 +1,157 @@ +import os + +from compass.landice.tests.ismip7_run.ismip7_ais.create_slm_mapping_files import ( # noqa + CreateSlmMappingFiles, +) +from compass.landice.tests.ismip7_run.ismip7_ais.set_up_experiment import ( + SetUpExperiment, +) +from compass.testcase import TestCase + +# Define the full experiment matrix per the ISMIP7 protocol +EXPERIMENTS = { + 'historical_CESM2-WACCM': { + 'scenario': 'historical', 'model': 'CESM2-WACCM', + 'start_time': '2000-01-01_00:00:00', + 'stop_time': '2015-01-01_00:00:00', + 'is_historical': True}, + 'historical_MRI-ESM2-0': { + 'scenario': 'historical', 'model': 'MRI-ESM2-0', + 'start_time': '2000-01-01_00:00:00', + 'stop_time': '2015-01-01_00:00:00', + 'is_historical': True}, + 'ssp370_CESM2-WACCM': { + 'scenario': 'ssp370', 'model': 'CESM2-WACCM', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2101-01-01_00:00:00', + 'is_historical': False}, + 'ssp370_MRI-ESM2-0': { + 'scenario': 'ssp370', 'model': 'MRI-ESM2-0', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2101-01-01_00:00:00', + 'is_historical': False}, + 'ssp126_CESM2-WACCM': { + 'scenario': 'ssp126', 'model': 'CESM2-WACCM', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2301-01-01_00:00:00', + 'is_historical': False}, + 'ssp126_MRI-ESM2-0': { + 'scenario': 'ssp126', 'model': 'MRI-ESM2-0', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2301-01-01_00:00:00', + 'is_historical': False}, + 'ssp585_CESM2-WACCM': { + 'scenario': 'ssp585', 'model': 'CESM2-WACCM', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2301-01-01_00:00:00', + 'is_historical': False}, + 'ssp585_MRI-ESM2-0': { + 'scenario': 'ssp585', 'model': 'MRI-ESM2-0', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2301-01-01_00:00:00', + 'is_historical': False}, + 'ctrl_CESM2-WACCM': { + 'scenario': 'ctrl', 'model': 'CESM2-WACCM', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2301-01-01_00:00:00', + 'is_historical': False}, + 'ctrl_MRI-ESM2-0': { + 'scenario': 'ctrl', 'model': 'MRI-ESM2-0', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2301-01-01_00:00:00', + 'is_historical': False}, + 'ocx': { + 'scenario': 'ocx', 'model': None, + 'start_time': '1990-01-01_00:00:00', + 'stop_time': '2026-01-01_00:00:00', + 'is_historical': True}, +} + + +class Ismip7Ais(TestCase): + """ + A test case for automated setup of a suite of standardized + ISMIP7 simulations for the Antarctic Ice Sheet. + """ + + def __init__(self, test_group): + """ + Create the test case + + Parameters + ---------- + test_group : compass.landice.tests.ismip7_run.Ismip7Run + The test group that this test case belongs to + """ + name = 'ismip7_ais' + super().__init__(test_group=test_group, name=name, subdir=name) + + def configure(self): + """ + Set up the desired ISMIP7 AIS experiments. + + Read the experiment list from config and add a SetUpExperiment + step for each. + """ + config = self.config + exp_list_str = config.get('ismip7_run_ais', 'exp_list') + + if exp_list_str == 'all': + exp_list = list(EXPERIMENTS.keys()) + elif exp_list_str == 'historical': + exp_list = [k for k, v in EXPERIMENTS.items() + if v['is_historical']] + elif exp_list_str == 'projections': + exp_list = [k for k, v in EXPERIMENTS.items() + if not v['is_historical'] and + v['scenario'] != 'ctrl'] + elif exp_list_str == 'ctrl': + exp_list = [k for k, v in EXPERIMENTS.items() + if v['scenario'] == 'ctrl'] + else: + exp_list = [s.strip() for s in exp_list_str.split(',')] + + for exp in exp_list: + if exp not in EXPERIMENTS: + raise ValueError( + f"Unknown experiment '{exp}'. Valid experiments: " + f"{list(EXPERIMENTS.keys())}") + if os.path.exists(os.path.join(self.work_dir, exp)): + print(f"WARNING: {exp} path already exists; skipping. " + "Remove the directory " + f"{os.path.join(self.work_dir, exp)} and run " + "'compass setup' again to recreate.") + else: + self.add_step( + SetUpExperiment(test_case=self, name=exp, + subdir=exp, exp=exp, + exp_info=EXPERIMENTS[exp])) + + # Do not add experiments to steps_to_run; + # each experiment (step) should be submitted manually + self.steps_to_run = [] + + # Optionally set up sea-level model mapping files + sea_level_model = config.getboolean('ismip7_run_ais', + 'sea_level_model') + if sea_level_model: + subdir = 'mapping_files' + if os.path.exists(os.path.join(self.work_dir, subdir)): + print(f"WARNING: {subdir} path already exists; skipping.") + else: + self.add_step( + CreateSlmMappingFiles(test_case=self, + name='mapping_files', + subdir=subdir)) + self.steps_to_run.append('mapping_files') + + def run(self): + """ + A dummy run method + """ + raise ValueError( + "ERROR: 'compass run' has no functionality at the test case " + "level for this test. Please submit the job script in each " + "experiment's subdirectory manually instead. " + "To create Sea-Level Model mapping files, submit job script " + "or execute 'compass run' from the 'mapping_files' subdirectory.") diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/albany_input.yaml b/compass/landice/tests/ismip7_run/ismip7_ais/albany_input.yaml new file mode 100644 index 0000000000..e1a4dcd7c4 --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_ais/albany_input.yaml @@ -0,0 +1,239 @@ +%YAML 1.1 +--- +ANONYMOUS: + Problem: + LandIce Field Norm: + sliding_velocity_basalside: + Regularization Type: Given Value + Regularization Value: 1.0e-4 + LandIce BCs: + BC 0: + Basal Friction Coefficient: + Type: Power Law + Power Exponent: 0.2 + Mu Type: Field + Effective Pressure Type: Constant + Effective Pressure: 1.0 + # Zero Effective Pressure On Floating Ice At Nodes: true + Zero Beta On Floating Ice: true + + Cubature Degree: 8 + +# Discretization Description + Discretization: + #Exodus Output File Name: albany_output.exo + + Piro: +# Nonlinear Solver Information + NOX: + Nonlinear Solver: Line Search Based + Line Search: + Full Step: + Full Step: 1.0e+00 + Method: Backtrack + Solver Options: + Status Test Check Type: Minimal + Status Tests: + Test Type: Combo + Combo Type: OR + Number of Tests: 2 + Test 0: + Test Type: NormF + Norm Type: Two Norm + Scale Type: Scaled + Tolerance: 1.0e-05 + Test 1: + Test Type: MaxIters + Maximum Iterations: 50 + Printing: + Output Precision: 3 + Output Processor: 0 + Output Information: + Error: true + Warning: true + Outer Iteration: true + Parameters: false + Details: false + Linear Solver Details: false + Stepper Iteration: true + Stepper Details: true + Stepper Parameters: true + + Direction: + Method: Newton + Newton: + Forcing Term Method: Constant + Rescue Bad Newton Solve: true + Linear Solver: + Write Linear System: false + Tolerance: 1.0e-8 + + Stratimikos Linear Solver: + Stratimikos: + +# Linear Solver Information + Linear Solver Type: Belos + Linear Solver Types: + Belos: + Solver Type: Block GMRES + Solver Types: + Block GMRES: + Output Frequency: 20 + Output Style: 1 + Verbosity: 33 + Maximum Iterations: 200 + Block Size: 1 + Num Blocks: 200 + Flexible Gmres: false + VerboseObject: + Output File: none + Verbosity Level: low + +# Preconditioner Information + Preconditioner Type: MueLu + Preconditioner Types: + + Ifpack2: + Overlap: 1 + Prec Type: ILUT + + MueLu: + Matrix: + PDE equations: 2 + Factories: + myLineDetectionFact: + factory: LineDetectionFactory + 'linedetection: orientation': coordinates + mySemiCoarsenPFact1: + factory: SemiCoarsenPFactory + 'semicoarsen: coarsen rate': 14 + UncoupledAggregationFact2: + factory: UncoupledAggregationFactory + 'aggregation: ordering': graph + 'aggregation: max selected neighbors': 0 + 'aggregation: min agg size': 3 + 'aggregation: phase3 avoid singletons': true + MyCoarseMap2: + factory: CoarseMapFactory + Aggregates: UncoupledAggregationFact2 + myTentativePFact2: + 'tentative: calculate qr': true + factory: TentativePFactory + Aggregates: UncoupledAggregationFact2 + CoarseMap: MyCoarseMap2 + mySaPFact2: + 'sa: eigenvalue estimate num iterations': 10 + 'sa: damping factor': 1.33333e+00 + factory: SaPFactory + P: myTentativePFact2 + myTransferCoordinatesFact: + factory: CoordinatesTransferFactory + CoarseMap: MyCoarseMap2 + Aggregates: UncoupledAggregationFact2 + myTogglePFact: + factory: TogglePFactory + 'semicoarsen: number of levels': 2 + TransferFactories: + P1: mySemiCoarsenPFact1 + P2: mySaPFact2 + Ptent1: mySemiCoarsenPFact1 + Ptent2: myTentativePFact2 + Nullspace1: mySemiCoarsenPFact1 + Nullspace2: myTentativePFact2 + myRestrictorFact: + factory: TransPFactory + P: myTogglePFact + myToggleTransferCoordinatesFact: + factory: ToggleCoordinatesTransferFactory + Chosen P: myTogglePFact + TransferFactories: + Coordinates1: mySemiCoarsenPFact1 + Coordinates2: myTransferCoordinatesFact + myRAPFact: + factory: RAPFactory + P: myTogglePFact + R: myRestrictorFact + TransferFactories: + For Coordinates: myToggleTransferCoordinatesFact + myRepartitionHeuristicFact: + factory: RepartitionHeuristicFactory + A: myRAPFact + 'repartition: min rows per proc': 3000 + 'repartition: max imbalance': 1.327e+00 + 'repartition: start level': 1 + myZoltanInterface: + factory: ZoltanInterface + A: myRAPFact + Coordinates: myToggleTransferCoordinatesFact + number of partitions: myRepartitionHeuristicFact + myRepartitionFact: + factory: RepartitionFactory + A: myRAPFact + Partition: myZoltanInterface + 'repartition: remap parts': true + number of partitions: myRepartitionHeuristicFact + myRebalanceProlongatorFact: + factory: RebalanceTransferFactory + type: Interpolation + P: myTogglePFact + Coordinates: myToggleTransferCoordinatesFact + Nullspace: myTogglePFact + myRebalanceRestrictionFact: + factory: RebalanceTransferFactory + type: Restriction + R: myRestrictorFact + myRebalanceAFact: + factory: RebalanceAcFactory + A: myRAPFact + TransferFactories: { } + mySmoother1: + factory: TrilinosSmoother + type: LINESMOOTHING_BANDEDRELAXATION + 'smoother: pre or post': both + ParameterList: + 'relaxation: type': Gauss-Seidel + 'relaxation: sweeps': 1 + 'relaxation: damping factor': 1.0 + mySmoother3: + factory: TrilinosSmoother + type: RELAXATION + 'smoother: pre or post': both + ParameterList: + 'relaxation: type': Gauss-Seidel + 'relaxation: sweeps': 1 + 'relaxation: damping factor': 1.0 + mySmoother4: + factory: TrilinosSmoother + type: RELAXATION + 'smoother: pre or post': pre + ParameterList: + 'relaxation: type': Gauss-Seidel + 'relaxation: sweeps': 4 + 'relaxation: damping factor': 1.0 + Hierarchy: + max levels: 7 + 'coarse: max size': 2000 + verbosity: None + Finest: + Smoother: mySmoother1 + CoarseSolver: mySmoother4 + P: myRebalanceProlongatorFact + Nullspace: myRebalanceProlongatorFact + CoarseNumZLayers: myLineDetectionFact + LineDetection_Layers: myLineDetectionFact + LineDetection_VertLineIds: myLineDetectionFact + A: myRebalanceAFact + Coordinates: myRebalanceProlongatorFact + Importer: myRepartitionFact + All: + startLevel: 1 + Smoother: mySmoother4 + CoarseSolver: mySmoother4 + P: myRebalanceProlongatorFact + Nullspace: myRebalanceProlongatorFact + CoarseNumZLayers: myLineDetectionFact + LineDetection_Layers: myLineDetectionFact + LineDetection_VertLineIds: myLineDetectionFact + A: myRebalanceAFact + Coordinates: myRebalanceProlongatorFact + Importer: myRepartitionFact diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/create_slm_mapping_files.py b/compass/landice/tests/ismip7_run/ismip7_ais/create_slm_mapping_files.py new file mode 100644 index 0000000000..108ee7544c --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_ais/create_slm_mapping_files.py @@ -0,0 +1,100 @@ +import shutil + +from mpas_tools.logging import check_call +from mpas_tools.scrip.from_mpas import scrip_from_mpas + +from compass.step import Step + + +class CreateSlmMappingFiles(Step): + """ + A step for creating mapping files for the Sea Level Model + """ + + def __init__(self, test_case, name, subdir): + """ + Initialize step + """ + super().__init__(test_case=test_case, name=name, subdir=subdir) + + def setup(self): + print(" Setting up mapping_file subdirectory") + + def run(self): + """ + Run this step of the test case + """ + config = self.config + logger = self.logger + section = config['ismip7_run_ais'] + sea_level_model = section.getboolean('sea_level_model') + if sea_level_model: + self._build_mapping_files(config, logger) + + def _build_mapping_files(self, config, logger): + """ + Build mapping files between the MALI mesh and the SLM grid. + """ + section = config['ismip7_run_ais'] + init_cond_path = section.get('init_cond_path') + nglv = section.getint('nglv') + section = config['parallel'] + ntasks = section.getint('cores_per_node') + + mali_scripfile = 'mali_scripfile.nc' + slm_scripfile = f'slm_nglv{nglv}scripfile.nc' + mali_meshfile = 'mali_meshfile_sphereLatLon.nc' + + # SLM scripfile + logger.info(f'creating scripfile for the SLM grid with ' + f'{nglv} Gauss-Legendre points in latitude') + + args = ['ncremap', + '-g', slm_scripfile, + '-G', + f'latlon={nglv},{2 * int(nglv)}#lat_typ=gss#lat_drc=n2s'] + + check_call(args, logger=logger) + + # MALI scripfile + shutil.copy(init_cond_path, mali_meshfile) + args = ['set_lat_lon_fields_in_planar_grid', + '--file', mali_meshfile, + '--proj', 'ais-bedmap2-sphere'] + + check_call(args, logger=logger) + + logger.info('creating scrip file for the mali mesh') + scrip_from_mpas(mali_meshfile, mali_scripfile) + + # MALI -> SLM mapping file + logger.info('creating MALI -> SLM grid mapfile with conserve method') + + parallel_executable = config.get("parallel", "parallel_executable") + args = parallel_executable.split(' ') + args.extend(['-n', f'{ntasks}', + 'ESMF_RegridWeightGen', + '-s', mali_scripfile, + '-d', slm_scripfile, + '-w', 'mapfile_mali_to_slm.nc', + '-m', 'conserve', + '-i', '-64bit_offset', '--netcdf4', + '--src_regional']) + + check_call(args, logger) + + # SLM -> MALI mapping file + logger.info('creating SLM -> MALI mesh mapfile with bilinear method') + args = parallel_executable.split(' ') + args.extend(['-n', f'{ntasks}', + 'ESMF_RegridWeightGen', + '-s', slm_scripfile, + '-d', mali_scripfile, + '-w', 'mapfile_slm_to_mali.nc', + '-m', 'bilinear', + '-i', '-64bit_offset', '--netcdf4', + '--dst_regional']) + + check_call(args, logger) + + logger.info('mapping file creation complete') diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais.cfg b/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais.cfg new file mode 100644 index 0000000000..dc26cc455f --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais.cfg @@ -0,0 +1,79 @@ +[ismip7_run_ais] + +# List of experiments to set up. +# Can be "all", "historical", "projections", "ctrl", or a +# comma-delimited list of experiment names. +# Valid experiment names: +# historical_CESM2-WACCM, historical_MRI-ESM2-0, +# ssp126_CESM2-WACCM, ssp126_MRI-ESM2-0, +# ssp370_CESM2-WACCM, ssp370_MRI-ESM2-0, +# ssp585_CESM2-WACCM, ssp585_MRI-ESM2-0, +# ctrl_CESM2-WACCM, ctrl_MRI-ESM2-0, +# ocx +exp_list = all + +# Number of tasks to use for each run +ntasks = 128 + +# Value to use for config_pio_stride. +# Should be divisible into ntasks +pio_stride = 128 + +# Base path to the pre-processed ISMIP7 forcing files. +# Expected layout: +# {forcing_basepath}/{model}_{scenario}/atmosphere/ +# {forcing_basepath}/{model}_{scenario}/ocean_thermal_forcing/ +# User has to supply. +forcing_basepath = NotAvailable + +# Path to the initial condition file. User has to supply. +init_cond_path = NotAvailable + +# Path to the file for the basal melt parametrization coefficients. +melt_params_path = NotAvailable + +# Path to the region mask file +region_mask_path = NotAvailable + +# Path to the ocean thermal forcing climatology file for CTRL2015 runs. +# This is the constant-climate 30-year mean (2000-2029) thermal forcing. +# User has to supply if running ctrl experiments. +ctrl_tf_climatology_path = NotAvailable + +# Path to the atmosphere climatology files for CTRL2015 runs. +# Directory containing constant-climate SMB, temperature, etc. +# User has to supply if running ctrl experiments. +ctrl_atm_climatology_path = NotAvailable + +# Path to OCX (observationally constrained experiment) forcing. +# User has to supply if running the ocx experiment. +ocx_forcing_path = NotAvailable + +# Calving method to use. Options: restore, von_mises +calving_method = restore + +# Path to the file containing the von Mises parameter fields. +# Only required if calving_method is set to 'von_mises'. +von_mises_parameter_path = NotAvailable + +# True if running coupled MALI-sea level model simulation +sea_level_model = false + +# Path to the directory containing globally defined ice thickness +# field for the sea-level model +slm_input_ice = NotAvailable + +# Path to the directory containing earth model for the sea-level model +slm_input_earth = NotAvailable + +# Earth structure profile filename +slm_earth_structure = prem_512.l60K2C.sum18p6.dum19p2.tz19p4.lm22 + +# Path to the directory containing other SLM input files +slm_input_others = NotAvailable + +# Number of gauss-legendre nodes in latitude (typically multiple of 512) +nglv = 2048 + +[parallel] +parallel_executable = srun --label --cpu-bind=cores diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais_test.cfg b/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais_test.cfg new file mode 100644 index 0000000000..4aa49a28f7 --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_ais/ismip7_ais_test.cfg @@ -0,0 +1,79 @@ +[ismip7_run_ais] + +# List of experiments to set up. +# Can be "all", "historical", "projections", "ctrl", or a +# comma-delimited list of experiment names. +# Valid experiment names: +# historical_CESM2-WACCM, historical_MRI-ESM2-0, +# ssp126_CESM2-WACCM, ssp126_MRI-ESM2-0, +# ssp370_CESM2-WACCM, ssp370_MRI-ESM2-0, +# ssp585_CESM2-WACCM, ssp585_MRI-ESM2-0, +# ctrl_CESM2-WACCM, ctrl_MRI-ESM2-0, +# ocx +exp_list = historical_CESM2-WACCM,ssp585_CESM2-WACCM + +# Number of tasks to use for each run +ntasks = 512 + +# Value to use for config_pio_stride. +# Should be divisible into ntasks +pio_stride = 128 + +# Base path to the pre-processed ISMIP7 forcing files. +# Expected layout: +# {forcing_basepath}/{model}_{scenario}/atmosphere/ +# {forcing_basepath}/{model}_{scenario}/ocean_thermal_forcing/ +# User has to supply. +forcing_basepath = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/test_processing/AIS/ + +# Path to the initial condition file. User has to supply. +init_cond_path = /global/cfs/cdirs/fanssie/MALI_projects/ISMIP6-2300/initial_conditions/AIS_4to20km_20230105/relaxation_0TGmelt_10yr/relaxed_10yrs_4km.nc + +# Path to the file for the basal melt parametrization coefficients. +melt_params_path = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/test_processing/AIS/climatology/basin_and_coeff_DeltaT_quadratic_non_local_gamma14500.nc + +# Path to the region mask file +region_mask_path = /global/cfs/cdirs/fanssie/MALI_projects/ISMIP6-2300/initial_conditions/AIS_4to20km_20230105/AIS_4to20km_r01_20220907.regionMask_ismip6.nc + +# Path to the ocean thermal forcing climatology file for CTRL2015 runs. +# This is the constant-climate 30-year mean (2000-2029) thermal forcing. +# User has to supply if running ctrl experiments. +ctrl_tf_climatology_path = NotAvailable + +# Path to the atmosphere climatology files for CTRL2015 runs. +# Directory containing constant-climate SMB, temperature, etc. +# User has to supply if running ctrl experiments. +ctrl_atm_climatology_path = NotAvailable + +# Path to OCX (observationally constrained experiment) forcing. +# User has to supply if running the ocx experiment. +ocx_forcing_path = NotAvailable + +# Calving method to use. Options: restore, von_mises +calving_method = restore + +# Path to the file containing the von Mises parameter fields. +# Only required if calving_method is set to 'von_mises'. +von_mises_parameter_path = NotAvailable + +# True if running coupled MALI-sea level model simulation +sea_level_model = false + +# Path to the directory containing globally defined ice thickness +# field for the sea-level model +slm_input_ice = NotAvailable + +# Path to the directory containing earth model for the sea-level model +slm_input_earth = NotAvailable + +# Earth structure profile filename +slm_earth_structure = prem_512.l60K2C.sum18p6.dum19p2.tz19p4.lm22 + +# Path to the directory containing other SLM input files +slm_input_others = NotAvailable + +# Number of gauss-legendre nodes in latitude (typically multiple of 512) +nglv = 2048 + +[parallel] +parallel_executable = srun --label --cpu-bind=cores diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/namelist.landice b/compass/landice/tests/ismip7_run/ismip7_ais/namelist.landice new file mode 100644 index 0000000000..03470752dc --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_ais/namelist.landice @@ -0,0 +1,86 @@ + config_velocity_solver = 'FO' + config_do_velocity_reconstruction_for_external_dycore = .false. + config_unrealistic_velocity = 00.00159 + config_nonconvergence_error = .false. + config_flowParamA_calculation = 'PB1982' + + config_thickness_advection = 'fo' + config_tracer_advection = 'fo' + config_zero_sfcMassBalApplied_over_bare_land = .true. + + config_uplift_method = 'none' + config_slm_coupling_interval = 5 + config_MALI_to_SLM_weights_file = 'mapfile_mali_to_slm.nc' + config_SLM_to_MALI_weights_file = 'mapfile_slm_to_mali.nc' + + config_calving = 'none' + config_apply_calving_mask = .false. + config_restore_calving_front_prevent_retreat = .false. + config_calculate_damage = .true. + config_damage_calving_threshold = 0.95 + config_damage_calving_method = 'none' + config_calving_speed_limit = 0.00063492063 + config_restore_calving_front = .true. + config_remove_icebergs = .true. + config_remove_small_islands = .true. + config_distribute_unablatedVolumeDynCell = .true. + config_calving_error_threshold = 100000.0 + + config_thermal_solver = 'temperature' + config_thermal_calculate_bmb = .true. + config_temperature_init = 'file' + config_thermal_thickness = 0.0 + config_surface_air_temperature_source = 'file' + config_basal_heat_flux_source = 'file' + + config_basal_mass_bal_float = 'ismip6' + config_front_mass_bal_grounded = 'ismip6' + config_use_3d_thermal_forcing_for_face_melt = .true. + config_add_ocean_thermal_forcing = 0.0 + + config_ice_density = 910.0 + config_ocean_density = 1028.0 + config_sea_level = 0.0 + config_flowLawExponent = 3.0 + config_dynamic_thickness = 10.0 + + config_dt = '0000-01-00_00:00:00' + config_time_integration = 'forward_euler' + config_adaptive_timestep = .true. + config_adaptive_timestep_calvingCFL_fraction = 0.8 + config_adaptive_timestep_include_calving = .true. + config_min_adaptive_timestep = 60 + config_max_adaptive_timestep = 3.154e7 + config_adaptive_timestep_CFL_fraction = 0.8 + config_adaptive_timestep_include_DCFL = .false. + config_adaptive_timestep_force_interval = '0000-01-00_00:00:00' + config_timeaveraging_interval = '0001-00-00_00:00:00' + config_enable_timeAvgRestarts = .false. + + config_do_restart = .true. + config_restart_timestamp_name = 'restart_timestamp' + config_start_time = 'file' + config_stop_time = '2301-01-01_00:00:00' + config_calendar_type = 'noleap' + + config_stats_interval = 0 + config_write_stats_on_startup = .false. + config_stats_cell_ID = 1 + config_write_output_on_startup = .true. + + config_always_compute_fem_grid = .true. + + config_ocean_connection_N = .false. + config_SGH = .false. + + config_AM_globalStats_enable = .true. + config_AM_globalStats_compute_interval = 'output_interval' + config_AM_globalStats_stream_name = 'globalStatsOutput' + config_AM_globalStats_compute_on_startup = .true. + config_AM_globalStats_write_on_startup = .true. + + config_AM_regionalStats_enable = .true. + config_AM_regionalStats_compute_interval = 'output_interval' + config_AM_regionalStats_stream_name = 'regionalStatsOutput' + config_AM_regionalStats_compute_on_startup = .true. + config_AM_regionalStats_write_on_startup = .true. diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/namelist.sealevel.template b/compass/landice/tests/ismip7_run/ismip7_ais/namelist.sealevel.template new file mode 100644 index 0000000000..e92a72ecab --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_ais/namelist.sealevel.template @@ -0,0 +1,68 @@ +&time_config + itersl = 1 + starttime = 2000 + dt1 = 5 + +/ +&model_resolution + norder = 512 + nglv = {{ nglv }} + +/ +&io_directory + inputfolder_ice = '{{ slm_input_ice }}' + inputfolder = '{{ slm_input_others }}' + planetfolder = '{{ slm_input_earth }}' + gridfolder = '{{ slm_input_others }}' + outputfolder = 'OUTPUT_SLM/' + outputfolder_ice = 'ICELOAD_SLM/' + folder_coupled = '' + +/ +&file_format + ext ='' + fType_in = 'text' + fType_out = 'both' + +/ +&file_name + planetmodel = '{{ slm_earth_structure }}' + icemodel = 'iceGlobalDomain_zeroField_GL{{ nglv }}_' + icemodel_out = 'iceload_out_' + timearray = 'times' + topomodel = 'etopo2_nglv{{ nglv }}_outside_AIS' + topo_initial = 'etopo2_nglv{{ nglv }}_outside_AIS' + grid_lat = 'GLlat_{{ nglv }}.txt' + grid_lon = 'GLlon_{{ nglv }}.txt' + +/ +&model_config + checkmarine = .false. + tpw = .true. + calcRG = .true. + input_times = .false. + initial_topo = .true. + iceVolume = .true. + coupling = .true. + patch_ice = .true. + +/ +&timewindow_config + L_sim = 300 + dt1 = 5 + dt2 = 10 + dt3 = 10 + dt4 = 10 + Ldt1 = 300 + Ldt2 = 0 + Ldt3 = 0 + Ldt4 = 0 + +/ +&others + whichplanet = 'earth' + + +/ + +! end of the namelist file diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/set_up_experiment.py b/compass/landice/tests/ismip7_run/ismip7_ais/set_up_experiment.py new file mode 100644 index 0000000000..7b1455c9ce --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_ais/set_up_experiment.py @@ -0,0 +1,396 @@ +import glob +import os +import sys +from importlib import resources + +from jinja2 import Template + +from compass.job import write_job_script +from compass.load_script import symlink_load_script +from compass.model import make_graph_file, run_model +from compass.step import Step + + +class SetUpExperiment(Step): + """ + A step for setting up an ISMIP7 AIS experiment + """ + + def __init__(self, test_case, name, subdir, exp, exp_info): + """ + Set up a new experiment + + Parameters + ---------- + test_case : compass.testcase.TestCase + The test case this step belongs to + + name : str + The name of this step (same as the experiment name) + + subdir : str + Subdirectory for this step + + exp : str + Experiment identifier (e.g., 'ssp585_CESM2-WACCM') + + exp_info : dict + Dictionary with experiment metadata: + scenario, model, start_time, stop_time, is_historical + """ + self.exp = exp + self.exp_info = exp_info + + super().__init__(test_case=test_case, name=name, subdir=subdir) + + def setup(self): # noqa: C901 + """ + Set up the experiment directory with all needed files. + """ + print(f" Setting up experiment {self.exp}") + + config = self.config + section = config['ismip7_run_ais'] + self.ntasks = section.getint('ntasks') + self.min_tasks = self.ntasks + forcing_basepath = section.get('forcing_basepath') + init_cond_path = section.get('init_cond_path') + init_cond_fname = os.path.split(init_cond_path)[-1] + melt_params_path = section.get('melt_params_path') + melt_params_fname = os.path.split(melt_params_path)[-1] + region_mask_path = section.get('region_mask_path') + region_mask_fname = os.path.split(region_mask_path)[-1] + calving_method = section.get('calving_method') + sea_level_model = section.getboolean('sea_level_model') + + exp_info = self.exp_info + scenario = exp_info['scenario'] + model = exp_info['model'] + is_historical = exp_info['is_historical'] + start_time = exp_info['start_time'] + stop_time = exp_info['stop_time'] + + # Define where to get templates (in current package) + resource_location = 'compass.landice.tests.ismip7_run.ismip7_ais' + + # Define calving method + use_vM_calving = (calving_method == 'von_mises') + + # --- Determine forcing file paths --- + if scenario == 'ocx': + ocx_forcing_path = section.get('ocx_forcing_path') + forcing_dir = ocx_forcing_path + elif scenario == 'ctrl': + # Control run uses climatology forcing + forcing_dir = None # handled separately below + else: + # Standard ESM-driven experiment + forcing_dir = os.path.join(forcing_basepath, + f"{model}_{scenario}") + + # --- Symlink input files --- + if is_historical: + os.symlink(init_cond_path, + os.path.join(self.work_dir, + os.path.basename(init_cond_path))) + os.symlink(melt_params_path, + os.path.join(self.work_dir, + os.path.basename(melt_params_path))) + os.symlink(region_mask_path, + os.path.join(self.work_dir, + os.path.basename(region_mask_path))) + + # --- Find and symlink forcing files --- + if scenario == 'ctrl': + # Control run: use climatology files + ctrl_tf_path = section.get('ctrl_tf_climatology_path') + ctrl_atm_path = section.get('ctrl_atm_climatology_path') + tf_fname = os.path.split(ctrl_tf_path)[-1] + os.symlink(ctrl_tf_path, + os.path.join(self.work_dir, tf_fname)) + + # Find atmosphere climatology files + smb_files = glob.glob(os.path.join(ctrl_atm_path, '*SMB*.nc')) + smb_files = [f for f in smb_files if 'gradient' not in f] + if len(smb_files) == 1: + smb_fname = os.path.split(smb_files[0])[-1] + os.symlink(smb_files[0], + os.path.join(self.work_dir, smb_fname)) + else: + sys.exit(f"ERROR: Expected 1 SMB climatology file in " + f"{ctrl_atm_path}, found {len(smb_files)}") + + temp_files = glob.glob( + os.path.join(ctrl_atm_path, '*temperature*.nc')) + temp_files = [f for f in temp_files if 'gradient' not in f] + if len(temp_files) == 1: + temp_fname = os.path.split(temp_files[0])[-1] + os.symlink(temp_files[0], + os.path.join(self.work_dir, temp_fname)) + else: + sys.exit(f"ERROR: Expected 1 temperature climatology file in " + f"{ctrl_atm_path}, found {len(temp_files)}") + + runoff_files = glob.glob( + os.path.join(ctrl_atm_path, '*runoff*.nc')) + if len(runoff_files) == 1: + runoff_fname = os.path.split(runoff_files[0])[-1] + os.symlink(runoff_files[0], + os.path.join(self.work_dir, runoff_fname)) + else: + runoff_fname = '' # runoff may not exist for ctrl + + smb_grad_files = glob.glob( + os.path.join(ctrl_atm_path, '*SMB_gradient*.nc')) + smb_grad_fname = '' + if len(smb_grad_files) == 1: + smb_grad_fname = os.path.split(smb_grad_files[0])[-1] + os.symlink(smb_grad_files[0], + os.path.join(self.work_dir, smb_grad_fname)) + + temp_grad_files = glob.glob( + os.path.join(ctrl_atm_path, '*temperature_gradient*.nc')) + temp_grad_fname = '' + if len(temp_grad_files) == 1: + temp_grad_fname = os.path.split(temp_grad_files[0])[-1] + os.symlink(temp_grad_files[0], + os.path.join(self.work_dir, temp_grad_fname)) + + else: + # Standard or OCX experiment: find forcing in forcing_dir + atm_dir = os.path.join(forcing_dir, 'atmosphere') + ocean_dir = os.path.join(forcing_dir, 'ocean_thermal_forcing') + + # SMB forcing + smb_search = os.path.join(atm_dir, '*SMB_*.nc') + smb_list = glob.glob(smb_search) + smb_list = [f for f in smb_list if 'gradient' not in f] + if len(smb_list) == 1: + smb_fname = os.path.split(smb_list[0])[-1] + os.symlink(smb_list[0], + os.path.join(self.work_dir, smb_fname)) + else: + sys.exit(f"ERROR: Expected 1 SMB file at {smb_search}, " + f"found {len(smb_list)}: {smb_list}") + + # Temperature forcing + temp_search = os.path.join(atm_dir, '*temperature_*.nc') + temp_list = glob.glob(temp_search) + temp_list = [f for f in temp_list if 'gradient' not in f] + if len(temp_list) == 1: + temp_fname = os.path.split(temp_list[0])[-1] + os.symlink(temp_list[0], + os.path.join(self.work_dir, temp_fname)) + else: + sys.exit(f"ERROR: Expected 1 temperature file at " + f"{temp_search}, found {len(temp_list)}") + + # Runoff forcing (optional — may not exist for all experiments) + runoff_search = os.path.join(atm_dir, '*runoff_*.nc') + runoff_list = glob.glob(runoff_search) + runoff_fname = '' + if len(runoff_list) == 1: + runoff_fname = os.path.split(runoff_list[0])[-1] + os.symlink(runoff_list[0], + os.path.join(self.work_dir, runoff_fname)) + + # SMB gradient (lapse rate) + smb_grad_search = os.path.join(atm_dir, '*SMB_gradient_*.nc') + smb_grad_list = glob.glob(smb_grad_search) + smb_grad_fname = '' + if len(smb_grad_list) == 1: + smb_grad_fname = os.path.split(smb_grad_list[0])[-1] + os.symlink(smb_grad_list[0], + os.path.join(self.work_dir, smb_grad_fname)) + + # Temperature gradient (lapse rate) + temp_grad_search = os.path.join(atm_dir, + '*temperature_gradient_*.nc') + temp_grad_list = glob.glob(temp_grad_search) + temp_grad_fname = '' + if len(temp_grad_list) == 1: + temp_grad_fname = os.path.split(temp_grad_list[0])[-1] + os.symlink(temp_grad_list[0], + os.path.join(self.work_dir, temp_grad_fname)) + + # Thermal forcing + tf_search = os.path.join(ocean_dir, '*thermal_forcing_*.nc') + tf_list = glob.glob(tf_search) + if len(tf_list) == 1: + tf_fname = os.path.split(tf_list[0])[-1] + os.symlink(tf_list[0], + os.path.join(self.work_dir, tf_fname)) + else: + sys.exit(f"ERROR: Expected 1 TF file at {tf_search}, " + f"found {len(tf_list)}: {tf_list}") + + # --- Set up streams --- + # Determine forcing interval + if scenario == 'ctrl': + forcing_interval_monthly = 'initial_only' + forcing_interval_annual = 'initial_only' + elif is_historical: + # Historical: read forcing at each interval + forcing_interval_monthly = '0000-01-00_00:00:00' + forcing_interval_annual = '0001-00-00_00:00:00' + else: + # Projections + forcing_interval_monthly = '0000-01-00_00:00:00' + forcing_interval_annual = '0001-00-00_00:00:00' + + stream_replacements = { + 'input_file_init_cond': init_cond_fname if is_historical + else 'USE_RESTART_FILE_INSTEAD', + 'input_file_region_mask': region_mask_fname if is_historical + else 'USE_RESTART_FILE_INSTEAD', + 'input_file_melt_params': melt_params_fname, + 'input_file_SMB_forcing': smb_fname, + 'input_file_temperature_forcing': temp_fname, + 'input_file_TF_forcing': tf_fname, + 'input_file_runoff_forcing': runoff_fname, + 'input_file_smb_gradient_forcing': smb_grad_fname, + 'input_file_temperature_gradient_forcing': temp_grad_fname, + 'forcing_interval_monthly': forcing_interval_monthly, + 'forcing_interval_annual': forcing_interval_annual, + } + + self.add_streams_file( + resource_location, + 'streams.landice.template', + out_name='streams.landice', + template_replacements=stream_replacements) + + # --- Set up namelist --- + self.add_namelist_file( + resource_location, 'namelist.landice', + out_name='namelist.landice') + + # PIO options + pio_stride = section.getint('pio_stride') + io_tasks = self.ntasks // pio_stride + options = {'config_pio_stride': f'{pio_stride}', + 'config_pio_num_iotasks': f'{io_tasks}'} + self.add_namelist_options(options=options, + out_name='namelist.landice') + + # Historical-specific options + if is_historical: + options = {'config_do_restart': ".false.", + 'config_start_time': f"'{start_time}'", + 'config_stop_time': f"'{stop_time}'"} + self.add_namelist_options(options=options, + out_name='namelist.landice') + else: + options = {'config_stop_time': f"'{stop_time}'"} + self.add_namelist_options(options=options, + out_name='namelist.landice') + + # Calving options + if use_vM_calving: + vM_path = section.get('von_mises_parameter_path') + options = { + 'config_calving': "'von_Mises_stress'", + 'config_restore_calving_front': ".false.", + 'config_floating_von_Mises_threshold_stress_source': "'data'", + 'config_grounded_von_Mises_threshold_stress_source': "'data'"} + self.add_namelist_options(options=options, + out_name='namelist.landice') + vM_stream_replacements = {'input_file_VM_params': vM_path} + self.add_streams_file( + resource_location, 'streams.vM_params', + out_name='streams.landice', + template_replacements=vM_stream_replacements) + + # Sea-level model options + if sea_level_model: + slm_input_ice = section.get('slm_input_ice') + slm_input_earth = section.get('slm_input_earth') + slm_earth_structure = section.get('slm_earth_structure') + slm_input_others = section.get('slm_input_others') + nglv = section.getint('nglv') + + slm_input_ice = os.path.join(slm_input_ice, + f'GL{nglv}/ice_noGrIS_GL{nglv}/') + slm_input_others = os.path.join(slm_input_others, + f'GL{nglv}/') + + options = {'config_uplift_method': "'sealevelmodel'"} + self.add_namelist_options(options=options, + out_name='namelist.landice') + + template = Template(resources.read_text( + resource_location, 'namelist.sealevel.template')) + text = template.render( + nglv=int(nglv), slm_input_ice=slm_input_ice, + slm_input_earth=slm_input_earth, + slm_earth_structure=slm_earth_structure, + slm_input_others=slm_input_others) + + file_slm_nl = os.path.join(self.work_dir, 'namelist.sealevel') + with open(file_slm_nl, 'w') as handle: + handle.write(text) + + os.makedirs(os.path.join(self.work_dir, 'OUTPUT_SLM/'), + exist_ok=True) + os.makedirs(os.path.join(self.work_dir, 'ICELOAD_SLM/'), + exist_ok=True) + + map_dir = os.path.join('..', 'mapping_files') + for map_file in ('mapfile_mali_to_slm.nc', + 'mapfile_slm_to_mali.nc'): + os.symlink(os.path.join(map_dir, map_file), + os.path.join(self.work_dir, map_file)) + + # --- Symlink restart for projections/ctrl --- + if not is_historical: + hist_exp = f"historical_{model}" + os.symlink(f"../{hist_exp}/rst.2015-01-01.nc", + os.path.join(self.work_dir, 'rst.2015-01-01.nc')) + with open(os.path.join(self.work_dir, "restart_timestamp"), + "w") as text_file: + text_file.write("2015-01-01_00:00:00") + + # --- Add albany yaml, graph file, load script, job script --- + self.add_input_file( + filename='albany_input.yaml', + package=resource_location, + copy=True) + + make_graph_file(mesh_filename=init_cond_path, + graph_filename=os.path.join(self.work_dir, + 'graph.info')) + + symlink_load_script(self.work_dir) + + self.config.set('job', 'job_name', self.exp) + machine = self.config.get('deploy', 'machine') + pre_run_cmd = ('LOGDIR=previous_logs_`date +"%Y-%m-%d_%H-%M-%S"`;' + 'mkdir $LOGDIR; cp log* $LOGDIR; date') + post_run_cmd = "date" + write_job_script(self.config, machine, + target_cores=self.ntasks, min_cores=self.min_tasks, + work_dir=self.work_dir, + pre_run_commands=pre_run_cmd, + post_run_commands=post_run_cmd) + + self.add_model_as_input() + + def run(self): + """ + Run this step of the test case + """ + config = self.config + section = config['ismip7_run_ais'] + sea_level_model = section.getboolean('sea_level_model') + if sea_level_model: + map_dir = os.path.join('..', 'mapping_files') + for map_file in ('mapfile_mali_to_slm.nc', + 'mapfile_slm_to_mali.nc'): + if not os.path.isfile(os.path.join(map_dir, map_file)): + sys.exit(f"ERROR: 'mapping_files/{map_file}' " + "does not exist in workdir. " + "Please run the 'mapping_files' step " + "before proceeding.") + + run_model(step=self, namelist='namelist.landice', + streams='streams.landice') diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/streams.landice.template b/compass/landice/tests/ismip7_run/ismip7_ais/streams.landice.template new file mode 100644 index 0000000000..ac143e5979 --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_ais/streams.landice.template @@ -0,0 +1,144 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/compass/landice/tests/ismip7_run/ismip7_ais/streams.vM_params b/compass/landice/tests/ismip7_run/ismip7_ais/streams.vM_params new file mode 100644 index 0000000000..d9ae1ca3f6 --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_ais/streams.vM_params @@ -0,0 +1,12 @@ + + + + + + + + diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/__init__.py b/compass/landice/tests/ismip7_run/ismip7_gris/__init__.py new file mode 100644 index 0000000000..ab883d15e2 --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_gris/__init__.py @@ -0,0 +1,134 @@ +import os + +from compass.landice.tests.ismip7_run.ismip7_gris.set_up_experiment import ( + SetUpExperiment, +) +from compass.testcase import TestCase + +# Define the full experiment matrix for GrIS +EXPERIMENTS = { + 'historical_CESM2-WACCM': { + 'scenario': 'historical', 'model': 'CESM2-WACCM', + 'start_time': '2007-01-01_00:00:00', + 'stop_time': '2015-01-01_00:00:00', + 'is_historical': True}, + 'historical_MRI-ESM2-0': { + 'scenario': 'historical', 'model': 'MRI-ESM2-0', + 'start_time': '2007-01-01_00:00:00', + 'stop_time': '2015-01-01_00:00:00', + 'is_historical': True}, + 'ssp370_CESM2-WACCM': { + 'scenario': 'ssp370', 'model': 'CESM2-WACCM', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2101-01-01_00:00:00', + 'is_historical': False}, + 'ssp370_MRI-ESM2-0': { + 'scenario': 'ssp370', 'model': 'MRI-ESM2-0', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2101-01-01_00:00:00', + 'is_historical': False}, + 'ssp126_CESM2-WACCM': { + 'scenario': 'ssp126', 'model': 'CESM2-WACCM', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2301-01-01_00:00:00', + 'is_historical': False}, + 'ssp126_MRI-ESM2-0': { + 'scenario': 'ssp126', 'model': 'MRI-ESM2-0', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2301-01-01_00:00:00', + 'is_historical': False}, + 'ssp585_CESM2-WACCM': { + 'scenario': 'ssp585', 'model': 'CESM2-WACCM', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2301-01-01_00:00:00', + 'is_historical': False}, + 'ssp585_MRI-ESM2-0': { + 'scenario': 'ssp585', 'model': 'MRI-ESM2-0', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2301-01-01_00:00:00', + 'is_historical': False}, + 'ctrl_CESM2-WACCM': { + 'scenario': 'ctrl', 'model': 'CESM2-WACCM', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2301-01-01_00:00:00', + 'is_historical': False}, + 'ctrl_MRI-ESM2-0': { + 'scenario': 'ctrl', 'model': 'MRI-ESM2-0', + 'start_time': '2015-01-01_00:00:00', + 'stop_time': '2301-01-01_00:00:00', + 'is_historical': False}, + 'ocx': { + 'scenario': 'ocx', 'model': None, + 'start_time': '1990-01-01_00:00:00', + 'stop_time': '2026-01-01_00:00:00', + 'is_historical': True}, +} + + +class Ismip7Gris(TestCase): + """ + A test case for automated setup of a suite of standardized + ISMIP7 simulations for the Greenland Ice Sheet. + """ + + def __init__(self, test_group): + """ + Create the test case + + Parameters + ---------- + test_group : compass.landice.tests.ismip7_run.Ismip7Run + The test group that this test case belongs to + """ + name = 'ismip7_gris' + super().__init__(test_group=test_group, name=name, subdir=name) + + def configure(self): + """ + Set up the desired ISMIP7 GrIS experiments. + """ + config = self.config + exp_list_str = config.get('ismip7_run_gris', 'exp_list') + + if exp_list_str == 'all': + exp_list = list(EXPERIMENTS.keys()) + elif exp_list_str == 'historical': + exp_list = [k for k, v in EXPERIMENTS.items() + if v['is_historical']] + elif exp_list_str == 'projections': + exp_list = [k for k, v in EXPERIMENTS.items() + if not v['is_historical'] and + v['scenario'] != 'ctrl'] + elif exp_list_str == 'ctrl': + exp_list = [k for k, v in EXPERIMENTS.items() + if v['scenario'] == 'ctrl'] + else: + exp_list = [s.strip() for s in exp_list_str.split(',')] + + for exp in exp_list: + if exp not in EXPERIMENTS: + raise ValueError( + f"Unknown experiment '{exp}'. Valid experiments: " + f"{list(EXPERIMENTS.keys())}") + if os.path.exists(os.path.join(self.work_dir, exp)): + print(f"WARNING: {exp} path already exists; skipping. " + "Remove the directory " + f"{os.path.join(self.work_dir, exp)} and run " + "'compass setup' again to recreate.") + else: + self.add_step( + SetUpExperiment(test_case=self, name=exp, + subdir=exp, exp=exp, + exp_info=EXPERIMENTS[exp])) + + # Do not add experiments to steps_to_run + self.steps_to_run = [] + + def run(self): + """ + A dummy run method + """ + raise ValueError( + "ERROR: 'compass run' has no functionality at the test case " + "level for this test. Please submit the job script in each " + "experiment's subdirectory manually instead.") diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/albany_input.yaml b/compass/landice/tests/ismip7_run/ismip7_gris/albany_input.yaml new file mode 100644 index 0000000000..bdd3fa7108 --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_gris/albany_input.yaml @@ -0,0 +1,238 @@ +%YAML 1.1 +--- +ANONYMOUS: + Problem: + Depth Integrated Model: true + LandIce Field Norm: + sliding_velocity_basalside: + Regularization Type: Given Value + Regularization Value: 1.0e-4 + LandIce BCs: + BC 0: + Basal Friction Coefficient: + Type: Power Law + Power Exponent: 1.0 + Mu Type: Field + Effective Pressure Type: Hydrostatic Computed At Nodes + Zero Effective Pressure On Floating Ice At Nodes: true + Zero Beta On Floating Ice: false + Use Pressurized Bed Above Sea Level: false + +# Discretization Description + Discretization: + #Exodus Output File Name: albany_output.exo + + Piro: +# Nonlinear Solver Information + NOX: + Nonlinear Solver: Line Search Based + Line Search: + Full Step: + Full Step: 1.0e+00 + Method: Backtrack + Solver Options: + Status Test Check Type: Minimal + Status Tests: + Test Type: Combo + Combo Type: OR + Number of Tests: 2 + Test 0: + Test Type: NormF + Norm Type: Two Norm + Scale Type: Scaled + Tolerance: 1.0e-05 + Test 1: + Test Type: MaxIters + Maximum Iterations: 50 + Printing: + Output Precision: 3 + Output Processor: 0 + Output Information: + Error: true + Warning: true + Outer Iteration: true + Parameters: false + Details: false + Linear Solver Details: false + Stepper Iteration: true + Stepper Details: true + Stepper Parameters: true + + Direction: + Method: Newton + Newton: + Forcing Term Method: Constant + Rescue Bad Newton Solve: true + Linear Solver: + Write Linear System: false + Tolerance: 1.0e-8 + + Stratimikos Linear Solver: + Stratimikos: + +# Linear Solver Information + Linear Solver Type: Belos + Linear Solver Types: + Belos: + Solver Type: Block GMRES + Solver Types: + Block GMRES: + Output Frequency: 20 + Output Style: 1 + Verbosity: 33 + Maximum Iterations: 200 + Block Size: 1 + Num Blocks: 200 + Flexible Gmres: false + VerboseObject: + Output File: none + Verbosity Level: low + +# Preconditioner Information + Preconditioner Type: MueLu + Preconditioner Types: + + Ifpack2: + Overlap: 1 + Prec Type: ILUT + + MueLu: + Matrix: + PDE equations: 2 + Factories: + myLineDetectionFact: + factory: LineDetectionFactory + 'linedetection: orientation': coordinates + mySemiCoarsenPFact1: + factory: SemiCoarsenPFactory + 'semicoarsen: coarsen rate': 14 + UncoupledAggregationFact2: + factory: UncoupledAggregationFactory + 'aggregation: ordering': graph + 'aggregation: max selected neighbors': 0 + 'aggregation: min agg size': 3 + 'aggregation: phase3 avoid singletons': true + MyCoarseMap2: + factory: CoarseMapFactory + Aggregates: UncoupledAggregationFact2 + myTentativePFact2: + 'tentative: calculate qr': true + factory: TentativePFactory + Aggregates: UncoupledAggregationFact2 + CoarseMap: MyCoarseMap2 + mySaPFact2: + 'sa: eigenvalue estimate num iterations': 10 + 'sa: damping factor': 1.33333e+00 + factory: SaPFactory + P: myTentativePFact2 + myTransferCoordinatesFact: + factory: CoordinatesTransferFactory + CoarseMap: MyCoarseMap2 + Aggregates: UncoupledAggregationFact2 + myTogglePFact: + factory: TogglePFactory + 'semicoarsen: number of levels': 2 + TransferFactories: + P1: mySemiCoarsenPFact1 + P2: mySaPFact2 + Ptent1: mySemiCoarsenPFact1 + Ptent2: myTentativePFact2 + Nullspace1: mySemiCoarsenPFact1 + Nullspace2: myTentativePFact2 + myRestrictorFact: + factory: TransPFactory + P: myTogglePFact + myToggleTransferCoordinatesFact: + factory: ToggleCoordinatesTransferFactory + Chosen P: myTogglePFact + TransferFactories: + Coordinates1: mySemiCoarsenPFact1 + Coordinates2: myTransferCoordinatesFact + myRAPFact: + factory: RAPFactory + P: myTogglePFact + R: myRestrictorFact + TransferFactories: + For Coordinates: myToggleTransferCoordinatesFact + myRepartitionHeuristicFact: + factory: RepartitionHeuristicFactory + A: myRAPFact + 'repartition: min rows per proc': 3000 + 'repartition: max imbalance': 1.327e+00 + 'repartition: start level': 1 + myZoltanInterface: + factory: ZoltanInterface + A: myRAPFact + Coordinates: myToggleTransferCoordinatesFact + number of partitions: myRepartitionHeuristicFact + myRepartitionFact: + factory: RepartitionFactory + A: myRAPFact + Partition: myZoltanInterface + 'repartition: remap parts': true + number of partitions: myRepartitionHeuristicFact + myRebalanceProlongatorFact: + factory: RebalanceTransferFactory + type: Interpolation + P: myTogglePFact + Coordinates: myToggleTransferCoordinatesFact + Nullspace: myTogglePFact + myRebalanceRestrictionFact: + factory: RebalanceTransferFactory + type: Restriction + R: myRestrictorFact + myRebalanceAFact: + factory: RebalanceAcFactory + A: myRAPFact + TransferFactories: { } + mySmoother1: + factory: TrilinosSmoother + type: LINESMOOTHING_BANDEDRELAXATION + 'smoother: pre or post': both + ParameterList: + 'relaxation: type': Gauss-Seidel + 'relaxation: sweeps': 1 + 'relaxation: damping factor': 1.0 + mySmoother3: + factory: TrilinosSmoother + type: RELAXATION + 'smoother: pre or post': both + ParameterList: + 'relaxation: type': Gauss-Seidel + 'relaxation: sweeps': 1 + 'relaxation: damping factor': 1.0 + mySmoother4: + factory: TrilinosSmoother + type: RELAXATION + 'smoother: pre or post': pre + ParameterList: + 'relaxation: type': Gauss-Seidel + 'relaxation: sweeps': 4 + 'relaxation: damping factor': 1.0 + Hierarchy: + max levels: 7 + 'coarse: max size': 2000 + verbosity: None + Finest: + Smoother: mySmoother1 + CoarseSolver: mySmoother4 + P: myRebalanceProlongatorFact + Nullspace: myRebalanceProlongatorFact + CoarseNumZLayers: myLineDetectionFact + LineDetection_Layers: myLineDetectionFact + LineDetection_VertLineIds: myLineDetectionFact + A: myRebalanceAFact + Coordinates: myRebalanceProlongatorFact + Importer: myRepartitionFact + All: + startLevel: 1 + Smoother: mySmoother4 + CoarseSolver: mySmoother4 + P: myRebalanceProlongatorFact + Nullspace: myRebalanceProlongatorFact + CoarseNumZLayers: myLineDetectionFact + LineDetection_Layers: myLineDetectionFact + LineDetection_VertLineIds: myLineDetectionFact + A: myRebalanceAFact + Coordinates: myRebalanceProlongatorFact + Importer: myRepartitionFact diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris.cfg b/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris.cfg new file mode 100644 index 0000000000..6a3f0a89a2 --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris.cfg @@ -0,0 +1,46 @@ +[ismip7_run_gris] + +# List of experiments to set up. +# Can be "all", "historical", "projections", "ctrl", or a +# comma-delimited list of experiment names. +exp_list = all + +# Number of tasks to use for each run +ntasks = 128 + +# Value to use for config_pio_stride. +pio_stride = 128 + +# Base path to the pre-processed ISMIP7 forcing files. +# Expected layout: +# {forcing_basepath}/{model}_{scenario}/atmosphere/ +# {forcing_basepath}/{model}_{scenario}/ocean_thermal_forcing/ +# User has to supply. +forcing_basepath = NotAvailable + +# Path to the initial condition file. User has to supply. +init_cond_path = NotAvailable + +# Path to the file for the basal melt parametrization coefficients. +melt_params_path = NotAvailable + +# Path to the region mask file +region_mask_path = NotAvailable + +# Path to the ocean thermal forcing climatology file for CTRL2015 runs. +ctrl_tf_climatology_path = NotAvailable + +# Path to the atmosphere climatology files for CTRL2015 runs. +ctrl_atm_climatology_path = NotAvailable + +# Path to OCX forcing. +ocx_forcing_path = NotAvailable + +# Calving method to use. Options: restore, von_mises +calving_method = restore + +# Path to the von Mises parameter fields. +von_mises_parameter_path = NotAvailable + +[parallel] +parallel_executable = srun --label --cpu-bind=cores diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris_test.cfg b/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris_test.cfg new file mode 100644 index 0000000000..66343580fb --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_gris/ismip7_gris_test.cfg @@ -0,0 +1,46 @@ +[ismip7_run_gris] + +# List of experiments to set up. +# Can be "all", "historical", "projections", "ctrl", or a +# comma-delimited list of experiment names. +exp_list = historical_CESM2-WACCM + +# Number of tasks to use for each run +ntasks = 512 + +# Value to use for config_pio_stride. +pio_stride = 128 + +# Base path to the pre-processed ISMIP7 forcing files. +# Expected layout: +# {forcing_basepath}/{model}_{scenario}/atmosphere/ +# {forcing_basepath}/{model}_{scenario}/ocean_thermal_forcing/ +# User has to supply. +forcing_basepath = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/test_processing/GIS + +# Path to the initial condition file. User has to supply. +init_cond_path = /global/cfs/cdirs/fanssie/MALI_input_files/GIS_1to10km_r02/GIS_1to10km_r02_20230202_remove_icebergs.nc + +# Path to the file for the basal melt parametrization coefficients. +melt_params_path = NotAvailable + +# Path to the region mask file +region_mask_path = /global/cfs/cdirs/fanssie/MALI_input_files/GIS_1to10km_r02/GIS_1to10km_r02_20230202_ismip6_regionMasks.nc + +# Path to the ocean thermal forcing climatology file for CTRL2015 runs. +ctrl_tf_climatology_path = NotAvailable + +# Path to the atmosphere climatology files for CTRL2015 runs. +ctrl_atm_climatology_path = NotAvailable + +# Path to OCX forcing. +ocx_forcing_path = NotAvailable + +# Calving method to use. Options: restore, von_mises +calving_method = restore + +# Path to the von Mises parameter fields. +von_mises_parameter_path = NotAvailable + +[parallel] +parallel_executable = srun --label --cpu-bind=cores diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/namelist.landice b/compass/landice/tests/ismip7_run/ismip7_gris/namelist.landice new file mode 100644 index 0000000000..f21f43f7c5 --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_gris/namelist.landice @@ -0,0 +1,89 @@ + config_velocity_solver = 'FO' + config_do_velocity_reconstruction_for_external_dycore = .false. + config_unrealistic_velocity = 00.00159 + config_nonconvergence_error = .false. + config_flowParamA_calculation = 'PB1982' + + config_thickness_advection = 'fct' + config_tracer_advection = 'fct' + config_horiz_tracer_adv_order = 3 + config_advection_coef_3rd_order = 1.0 + config_zero_sfcMassBalApplied_over_bare_land = .true. + + config_uplift_method = 'none' + + config_calving = 'crevasse_depth' + config_apply_calving_mask = .false. + config_restore_calving_front_prevent_retreat = .false. + config_calculate_damage = .false. + config_damage_calving_threshold = 0.95 + config_damage_calving_method = 'none' + config_calving_speed_limit = 0.00063492063 + config_restore_calving_front = .false. + config_remove_icebergs = .true. + config_remove_small_islands = .true. + config_distribute_unablatedVolumeDynCell = .true. + config_calving_error_threshold = 100000.0 + config_apply_facemelt_strainrate_enhancement = .true. + config_calving_strainrate_scaling = 5.0 + + + config_thermal_solver = 'temperature' + config_thermal_calculate_bmb = .true. + config_temperature_init = 'file' + config_thermal_thickness = 0.0 + config_surface_air_temperature_source = 'file' + config_basal_heat_flux_source = 'file' + + config_basal_mass_bal_float = 'none' + config_front_mass_bal_grounded = 'ismip6' + config_use_3d_thermal_forcing_for_face_melt = .false. + config_add_ocean_thermal_forcing = 0.0 + + config_ice_density = 910.0 + config_ocean_density = 1028.0 + config_sea_level = 0.0 + config_flowLawExponent = 3.0 + config_dynamic_thickness = 10.0 + + config_dt = '0000-01-00_00:00:00' + config_time_integration = 'runge_kutta' + config_rk_order = 2 + config_adaptive_timestep = .true. + config_adaptive_timestep_calvingCFL_fraction = 0.8 + config_adaptive_timestep_include_calving = .true. + config_min_adaptive_timestep = 60 + config_max_adaptive_timestep = 3.154e7 + config_adaptive_timestep_CFL_fraction = 0.5 + config_adaptive_timestep_include_DCFL = .false. + config_adaptive_timestep_force_interval = '0000-01-00_00:00:00' + config_timeaveraging_interval = '0001-00-00_00:00:00' + config_enable_timeAvgRestarts = .false. + + config_do_restart = .true. + config_restart_timestamp_name = 'restart_timestamp' + config_start_time = 'file' + config_stop_time = '2301-01-01_00:00:00' + config_calendar_type = 'noleap' + + config_stats_interval = 0 + config_write_stats_on_startup = .false. + config_stats_cell_ID = 1 + config_write_output_on_startup = .true. + + config_always_compute_fem_grid = .true. + + config_ocean_connection_N = .false. + config_SGH = .false. + + config_AM_globalStats_enable = .true. + config_AM_globalStats_compute_interval = 'output_interval' + config_AM_globalStats_stream_name = 'globalStatsOutput' + config_AM_globalStats_compute_on_startup = .true. + config_AM_globalStats_write_on_startup = .true. + + config_AM_regionalStats_enable = .true. + config_AM_regionalStats_compute_interval = 'output_interval' + config_AM_regionalStats_stream_name = 'regionalStatsOutput' + config_AM_regionalStats_compute_on_startup = .true. + config_AM_regionalStats_write_on_startup = .true. diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py b/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py new file mode 100644 index 0000000000..c07a79bbbf --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_gris/set_up_experiment.py @@ -0,0 +1,317 @@ +import glob +import os +import sys + +from compass.job import write_job_script +from compass.load_script import symlink_load_script +from compass.model import make_graph_file, run_model +from compass.step import Step + + +class SetUpExperiment(Step): + """ + A step for setting up an ISMIP7 GrIS experiment + """ + + def __init__(self, test_case, name, subdir, exp, exp_info): + """ + Set up a new experiment + + Parameters + ---------- + test_case : compass.testcase.TestCase + The test case this step belongs to + + name : str + The name of this step (same as the experiment name) + + subdir : str + Subdirectory for this step + + exp : str + Experiment identifier + + exp_info : dict + Dictionary with experiment metadata + """ + self.exp = exp + self.exp_info = exp_info + + super().__init__(test_case=test_case, name=name, subdir=subdir) + + def setup(self): # noqa: C901 + """ + Set up the experiment directory with all needed files. + """ + print(f" Setting up experiment {self.exp}") + + config = self.config + section = config['ismip7_run_gris'] + self.ntasks = section.getint('ntasks') + self.min_tasks = self.ntasks + forcing_basepath = section.get('forcing_basepath') + init_cond_path = section.get('init_cond_path') + init_cond_fname = os.path.split(init_cond_path)[-1] + melt_params_path = section.get('melt_params_path') + melt_params_fname = os.path.split(melt_params_path)[-1] + region_mask_path = section.get('region_mask_path') + region_mask_fname = os.path.split(region_mask_path)[-1] + calving_method = section.get('calving_method') + + exp_info = self.exp_info + scenario = exp_info['scenario'] + model = exp_info['model'] + is_historical = exp_info['is_historical'] + start_time = exp_info['start_time'] + stop_time = exp_info['stop_time'] + + resource_location = 'compass.landice.tests.ismip7_run.ismip7_gris' + + use_vM_calving = (calving_method == 'von_mises') + + # --- Determine forcing file paths --- + if scenario == 'ocx': + ocx_forcing_path = section.get('ocx_forcing_path') + forcing_dir = ocx_forcing_path + elif scenario == 'ctrl': + forcing_dir = None + else: + forcing_dir = os.path.join(forcing_basepath, + f"{model}_{scenario}") + + # --- Symlink input files --- + if is_historical: + os.symlink(init_cond_path, + os.path.join(self.work_dir, + os.path.basename(init_cond_path))) + os.symlink(melt_params_path, + os.path.join(self.work_dir, + os.path.basename(melt_params_path))) + os.symlink(region_mask_path, + os.path.join(self.work_dir, + os.path.basename(region_mask_path))) + + # --- Find and symlink forcing files --- + if scenario == 'ctrl': + ctrl_tf_path = section.get('ctrl_tf_climatology_path') + ctrl_atm_path = section.get('ctrl_atm_climatology_path') + tf_fname = os.path.split(ctrl_tf_path)[-1] + os.symlink(ctrl_tf_path, + os.path.join(self.work_dir, tf_fname)) + + smb_files = glob.glob(os.path.join(ctrl_atm_path, '*SMB*.nc')) + smb_files = [f for f in smb_files if 'gradient' not in f] + if len(smb_files) == 1: + smb_fname = os.path.split(smb_files[0])[-1] + os.symlink(smb_files[0], + os.path.join(self.work_dir, smb_fname)) + else: + sys.exit(f"ERROR: Expected 1 SMB climatology file in " + f"{ctrl_atm_path}, found {len(smb_files)}") + + temp_files = glob.glob( + os.path.join(ctrl_atm_path, '*temperature*.nc')) + temp_files = [f for f in temp_files if 'gradient' not in f] + if len(temp_files) == 1: + temp_fname = os.path.split(temp_files[0])[-1] + os.symlink(temp_files[0], + os.path.join(self.work_dir, temp_fname)) + else: + sys.exit(f"ERROR: Expected 1 temperature climatology file in " + f"{ctrl_atm_path}, found {len(temp_files)}") + + runoff_files = glob.glob( + os.path.join(ctrl_atm_path, '*runoff*.nc')) + runoff_fname = '' + if len(runoff_files) == 1: + runoff_fname = os.path.split(runoff_files[0])[-1] + os.symlink(runoff_files[0], + os.path.join(self.work_dir, runoff_fname)) + + smb_grad_files = glob.glob( + os.path.join(ctrl_atm_path, '*SMB_gradient*.nc')) + smb_grad_fname = '' + if len(smb_grad_files) == 1: + smb_grad_fname = os.path.split(smb_grad_files[0])[-1] + os.symlink(smb_grad_files[0], + os.path.join(self.work_dir, smb_grad_fname)) + + temp_grad_files = glob.glob( + os.path.join(ctrl_atm_path, '*temperature_gradient*.nc')) + temp_grad_fname = '' + if len(temp_grad_files) == 1: + temp_grad_fname = os.path.split(temp_grad_files[0])[-1] + os.symlink(temp_grad_files[0], + os.path.join(self.work_dir, temp_grad_fname)) + + else: + atm_dir = os.path.join(forcing_dir, 'atmosphere') + ocean_dir = os.path.join(forcing_dir, 'ocean_thermal_forcing') + + smb_search = os.path.join(atm_dir, '*SMB_*.nc') + smb_list = glob.glob(smb_search) + smb_list = [f for f in smb_list if 'gradient' not in f] + if len(smb_list) == 1: + smb_fname = os.path.split(smb_list[0])[-1] + os.symlink(smb_list[0], + os.path.join(self.work_dir, smb_fname)) + else: + sys.exit(f"ERROR: Expected 1 SMB file at {smb_search}, " + f"found {len(smb_list)}") + + temp_search = os.path.join(atm_dir, '*temperature_*.nc') + temp_list = glob.glob(temp_search) + temp_list = [f for f in temp_list if 'gradient' not in f] + if len(temp_list) == 1: + temp_fname = os.path.split(temp_list[0])[-1] + os.symlink(temp_list[0], + os.path.join(self.work_dir, temp_fname)) + else: + sys.exit(f"ERROR: Expected 1 temperature file at " + f"{temp_search}, found {len(temp_list)}") + + runoff_search = os.path.join(atm_dir, '*runoff_*.nc') + runoff_list = glob.glob(runoff_search) + runoff_fname = '' + if len(runoff_list) == 1: + runoff_fname = os.path.split(runoff_list[0])[-1] + os.symlink(runoff_list[0], + os.path.join(self.work_dir, runoff_fname)) + + smb_grad_search = os.path.join(atm_dir, '*SMB_gradient_*.nc') + smb_grad_list = glob.glob(smb_grad_search) + smb_grad_fname = '' + if len(smb_grad_list) == 1: + smb_grad_fname = os.path.split(smb_grad_list[0])[-1] + os.symlink(smb_grad_list[0], + os.path.join(self.work_dir, smb_grad_fname)) + + temp_grad_search = os.path.join(atm_dir, + '*temperature_gradient_*.nc') + temp_grad_list = glob.glob(temp_grad_search) + temp_grad_fname = '' + if len(temp_grad_list) == 1: + temp_grad_fname = os.path.split(temp_grad_list[0])[-1] + os.symlink(temp_grad_list[0], + os.path.join(self.work_dir, temp_grad_fname)) + + # GrIS uses 2D thermal forcing + tf_search = os.path.join(ocean_dir, '*thermal_forcing_*.nc') + tf_list = glob.glob(tf_search) + if len(tf_list) == 1: + tf_fname = os.path.split(tf_list[0])[-1] + os.symlink(tf_list[0], + os.path.join(self.work_dir, tf_fname)) + else: + sys.exit(f"ERROR: Expected 1 TF file at {tf_search}, " + f"found {len(tf_list)}") + + # --- Set up streams --- + if scenario == 'ctrl': + forcing_interval_monthly = 'initial_only' + forcing_interval_annual = 'initial_only' + else: + forcing_interval_monthly = '0000-01-00_00:00:00' + forcing_interval_annual = '0001-00-00_00:00:00' + + stream_replacements = { + 'input_file_init_cond': init_cond_fname if is_historical + else 'USE_RESTART_FILE_INSTEAD', + 'input_file_region_mask': region_mask_fname if is_historical + else 'USE_RESTART_FILE_INSTEAD', + 'input_file_melt_params': melt_params_fname, + 'input_file_SMB_forcing': smb_fname, + 'input_file_temperature_forcing': temp_fname, + 'input_file_TF_forcing': tf_fname, + 'input_file_runoff_forcing': runoff_fname, + 'input_file_smb_gradient_forcing': smb_grad_fname, + 'input_file_temperature_gradient_forcing': temp_grad_fname, + 'forcing_interval_monthly': forcing_interval_monthly, + 'forcing_interval_annual': forcing_interval_annual, + } + + self.add_streams_file( + resource_location, + 'streams.landice.template', + out_name='streams.landice', + template_replacements=stream_replacements) + + # --- Set up namelist --- + self.add_namelist_file( + resource_location, 'namelist.landice', + out_name='namelist.landice') + + # PIO options + pio_stride = section.getint('pio_stride') + io_tasks = self.ntasks // pio_stride + options = {'config_pio_stride': f'{pio_stride}', + 'config_pio_num_iotasks': f'{io_tasks}'} + self.add_namelist_options(options=options, + out_name='namelist.landice') + + if is_historical: + options = {'config_do_restart': ".false.", + 'config_start_time': f"'{start_time}'", + 'config_stop_time': f"'{stop_time}'"} + self.add_namelist_options(options=options, + out_name='namelist.landice') + else: + options = {'config_stop_time': f"'{stop_time}'"} + self.add_namelist_options(options=options, + out_name='namelist.landice') + + if use_vM_calving: + vM_path = section.get('von_mises_parameter_path') + options = { + 'config_calving': "'von_Mises_stress'", + 'config_restore_calving_front': ".false.", + 'config_floating_von_Mises_threshold_stress_source': "'data'", + 'config_grounded_von_Mises_threshold_stress_source': "'data'"} + self.add_namelist_options(options=options, + out_name='namelist.landice') + vM_stream_replacements = {'input_file_VM_params': vM_path} + self.add_streams_file( + resource_location, 'streams.vM_params', + out_name='streams.landice', + template_replacements=vM_stream_replacements) + + # --- Symlink restart for projections/ctrl --- + if not is_historical: + hist_exp = f"historical_{model}" + os.symlink(f"../{hist_exp}/rst.2015-01-01.nc", + os.path.join(self.work_dir, 'rst.2015-01-01.nc')) + with open(os.path.join(self.work_dir, "restart_timestamp"), + "w") as text_file: + text_file.write("2015-01-01_00:00:00") + + # --- Add albany yaml, graph file, load script, job script --- + self.add_input_file( + filename='albany_input.yaml', + package=resource_location, + copy=True) + + make_graph_file(mesh_filename=init_cond_path, + graph_filename=os.path.join(self.work_dir, + 'graph.info')) + + symlink_load_script(self.work_dir) + + self.config.set('job', 'job_name', self.exp) + machine = self.config.get('deploy', 'machine') + pre_run_cmd = ('LOGDIR=previous_logs_`date +"%Y-%m-%d_%H-%M-%S"`;' + 'mkdir $LOGDIR; cp log* $LOGDIR; date') + post_run_cmd = "date" + write_job_script(self.config, machine, + target_cores=self.ntasks, min_cores=self.min_tasks, + work_dir=self.work_dir, + pre_run_commands=pre_run_cmd, + post_run_commands=post_run_cmd) + + self.add_model_as_input() + + def run(self): + """ + Run this step of the test case + """ + run_model(step=self, namelist='namelist.landice', + streams='streams.landice') diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/streams.landice.template b/compass/landice/tests/ismip7_run/ismip7_gris/streams.landice.template new file mode 100644 index 0000000000..20ba57c258 --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_gris/streams.landice.template @@ -0,0 +1,133 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/compass/landice/tests/ismip7_run/ismip7_gris/streams.vM_params b/compass/landice/tests/ismip7_run/ismip7_gris/streams.vM_params new file mode 100644 index 0000000000..d9ae1ca3f6 --- /dev/null +++ b/compass/landice/tests/ismip7_run/ismip7_gris/streams.vM_params @@ -0,0 +1,12 @@ + + + + + + + + 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..b0871ac2b5 100644 --- a/docs/developers_guide/landice/test_groups/index.rst +++ b/docs/developers_guide/landice/test_groups/index.rst @@ -20,6 +20,8 @@ Test groups hydro_radial ismip6_forcing ismip6_run + ismip7_forcing + ismip7_run 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..72aa8c7411 --- /dev/null +++ b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst @@ -0,0 +1,117 @@ +.. _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``. + +The ``run()`` method dispatches to two sub-methods based on the boolean config +options ``process_ocean_thermal`` and ``process_ocean_climatology`` in the +``[ismip7]`` section: + +* ``_run_scenario()``: Processes time-varying ESM scenario data (model + + scenario combination). Uses config from ``[ismip7_ocean_thermal]``. +* ``_run_climatology()``: Processes the static observational climatology + (Zhou et al., AIS only). Uses config from ``[ismip7_ocean_climatology]``. + The TF version (currently v3) is hard-coded. + +For AIS scenario data, the step: + +* Remaps thermal forcing preserving 30 vertical ocean layers +* Produces ``ismip6shelfMelt_3dThermalForcing`` (dims: Time × nCells × + nISMIP6OceanLayers) +* Includes depth coordinate variables ``ismip6shelfMelt_zOcean`` and + ``ismip6shelfMelt_zBndsOcean`` + +For AIS climatology data, the step: + +* Extrapolates fill values, remaps, and renames to MALI conventions +* Produces ``ismip6shelfMelt_3dThermalForcing`` (dims: nCells × + nISMIP6OceanLayers) — no Time dimension + +For GrIS, the step: + +* Remaps 2D monthly thermal forcing +* Produces ``ismip6_2dThermalForcing`` (dims: Time × nCells) diff --git a/docs/developers_guide/landice/test_groups/ismip7_run.rst b/docs/developers_guide/landice/test_groups/ismip7_run.rst new file mode 100644 index 0000000000..15199b4830 --- /dev/null +++ b/docs/developers_guide/landice/test_groups/ismip7_run.rst @@ -0,0 +1,116 @@ +.. _dev_landice_ismip7_run: + +ismip7_run +========== + +The ``ismip7_run`` test group +(:py:class:`compass.landice.tests.ismip7_run`) sets up experiments from +the ISMIP7 experimental protocol for both the Antarctic Ice Sheet (AIS) +and the Greenland Ice Sheet (GrIS). Optionally, the AIS test case +supports coupled MALI–Sea Level Model (SLM) simulations. +(see :ref:`landice_ismip7_run`). + +framework +--------- + +The ``ismip7_run`` test group +(:py:class:`compass.landice.tests.ismip7_run.Ismip7Run`) registers two +test cases: + +* :py:class:`compass.landice.tests.ismip7_run.ismip7_ais.Ismip7Ais` +* :py:class:`compass.landice.tests.ismip7_run.ismip7_gris.Ismip7Gris` + +There is no shared functionality between the two test cases at present. +Shared functions may be added in the future if the needed functionality +can be generalized. + +ismip7_ais +---------- + +The :py:class:`compass.landice.tests.ismip7_run.ismip7_ais.Ismip7Ais` +test case sets up an ensemble of ISMIP7 Antarctica simulations +(standalone MALI or coupled MALI-SLM). + +The constructor (``__init__``) does nothing other than allow the +``ismip7_ais`` test case to be listed by ``compass list`` without having +all individual experiments listed in a verbose listing. Each individual +experiment is a step rather than a test case to avoid excessive +subdirectories. + +The ``configure`` method parses the ``exp_list`` config option from the +``[ismip7_run_ais]`` section. It supports: + +* ``all`` — all 11 core experiments +* ``historical`` — just the two historical runs +* ``projections`` — the six SSP projection runs +* ``ctrl`` — the two CTRL2015 runs +* A comma-delimited list of specific experiment names + +Each selected experiment is added as a +:py:class:`~compass.landice.tests.ismip7_run.ismip7_ais.set_up_experiment.SetUpExperiment` +step and immediately removed from ``steps_to_run`` (experiments should be +submitted individually, not run through the test case). + +The ``run`` method raises an error instructing the user to submit batch +jobs for each experiment individually. + +set_up_experiment (AIS) +~~~~~~~~~~~~~~~~~~~~~~~ + +The class +:py:class:`compass.landice.tests.ismip7_run.ismip7_ais.set_up_experiment.SetUpExperiment` +defines a step for a single ISMIP7 AIS experiment. + +The ``setup`` method sets up the experiment directory by: + +1. Creating symlinks to forcing files from the conventional path layout + under ``forcing_basepath``. +2. Copying and populating the streams template with the correct forcing + filenames and intervals (monthly for SMB/temperature/runoff, annual + for lapse rates and thermal forcing, ``initial_only`` for melt + parameters). +3. Processing the namelist template for the experiment's time period + and restart frequency. +4. Adding calving-specific streams (face melting, von Mises params) if + configured. +5. Creating a restart symlink for projection experiments pointing to + the corresponding ESM's historical restart + (``../historical_{model}/rst.2015-01-01.nc``). +6. Setting up CTRL2015 experiments with constant-climate forcing + (``initial_only`` intervals). +7. Setting up the OCX experiment with reanalysis-based forcing. +8. If SLM coupling is enabled, adding a ``CreateSlmMappingFiles`` step + and writing the SLM namelist from the Jinja2 template. +9. Generating a ``graph.info`` file and a SLURM job script. +10. Symlinking the compass load script into the run directory. + +The ``run`` method executes MALI for the given experiment. + +create_slm_mapping_files +~~~~~~~~~~~~~~~~~~~~~~~~ + +The class +:py:class:`compass.landice.tests.ismip7_run.ismip7_ais.create_slm_mapping_files.CreateSlmMappingFiles` +creates mapping files between the MALI mesh and the SLM grid. This step +is only added when sea-level model coupling is enabled. + +ismip7_gris +----------- + +The :py:class:`compass.landice.tests.ismip7_run.ismip7_gris.Ismip7Gris` +test case mirrors ``ismip7_ais`` but for the Greenland Ice Sheet. + +Key differences from the AIS test case: + +* Ocean thermal forcing is 2D (depth-averaged) rather than 3D. +* No sea-level model coupling. +* Default calving method is ``von_mises``. +* Config section is ``[ismip7_run_gris]``. + +set_up_experiment (GrIS) +~~~~~~~~~~~~~~~~~~~~~~~~ + +The class +:py:class:`compass.landice.tests.ismip7_run.ismip7_gris.set_up_experiment.SetUpExperiment` +follows the same logic as the AIS version, with the differences noted +above (2D TF stream, no SLM support). diff --git a/docs/users_guide/landice/test_groups/index.rst b/docs/users_guide/landice/test_groups/index.rst index 80fed7af3a..6a49090c00 100644 --- a/docs/users_guide/landice/test_groups/index.rst +++ b/docs/users_guide/landice/test_groups/index.rst @@ -25,6 +25,8 @@ physics but that are not run routinely. hydro_radial ismip6_forcing ismip6_run + ismip7_forcing + ismip7_run 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..1740f95a5f --- /dev/null +++ b/docs/users_guide/landice/test_groups/ismip7_forcing.rst @@ -0,0 +1,240 @@ +.. _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. The step can also + process the observational ocean thermal forcing climatology (Zhou et al.) + for AIS, controlled by the ``process_ocean_climatology`` config option. + +(For more details on the steps of each test case, see +:ref:`landice_ismip7_forcing_atmosphere` and +: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 AIS ocean thermal climatology (8km, 30 depth levels, static): + +.. code-block:: none + + {base_path_climatology}/tf/v3/tf_AIS_obs_ocean_climatology_*.nc + +For GrIS atmosphere (1km, polar stereographic EPSG:3413): + +.. code-block:: none + + 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 four 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 + + # Whether to process time-varying ocean thermal forcing (ESM scenario data) + process_ocean_thermal = true + + # Whether to process observational ocean thermal forcing climatology + process_ocean_climatology = true + + # config options for ismip7 atmosphere forcing + [ismip7_atmosphere] + + # 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 + + # config options for ismip7 ocean thermal forcing climatology + [ismip7_ocean_climatology] + + # Remapping method: bilinear, neareststod, conserve + method_remap = bilinear + + # Base path to observational climatology data + base_path_climatology = /path/to/ISMIP7/forcing/AIS/obs/zhou_annual_06_nov + +All ``NotAvailable`` options must be overridden in a user config file passed +at setup time (e.g., ``compass setup ... -f my_ismip7.cfg``). + +The boolean options ``process_ocean_thermal`` and ``process_ocean_climatology`` +control which processing paths are executed when the ``ocean_thermal`` test +case is run. Both default to ``true``. Set one to ``false`` in your user +config to skip that processing path. + +.. _landice_ismip7_forcing_atmosphere: + +atmosphere +---------- + +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. + +The step supports two processing modes, controlled by boolean config options +in the ``[ismip7]`` section: + +* **Scenario (time-varying) data** (``process_ocean_thermal = true``): + Processes ESM-driven thermal forcing for a given model/scenario combination. + +* **Observational climatology** (``process_ocean_climatology = true``): + Processes the static Zhou et al. observational thermal forcing climatology + (AIS only). This is a time-invariant 3D field referenced to 1995-2024. + +Both modes can be enabled simultaneously. + +For **AIS** scenario data, thermal forcing is 3D with 30 vertical ocean +layers. The input files span decades (e.g., 1850-1859). The output variable +is ``ismip6shelfMelt_3dThermalForcing`` with dimension +``nISMIP6OceanLayers``. Associated depth coordinate variables +``ismip6shelfMelt_zOcean`` and ``ismip6shelfMelt_zBndsOcean`` are also +produced. + +For **AIS** climatology data, the output is the same 3D thermal forcing field +but without a Time dimension, producing a single static file. + +For **GrIS**, thermal forcing is 2D (depth-averaged), with monthly temporal +resolution and yearly input files. The output variable is +``ismip6_2dThermalForcing``. diff --git a/docs/users_guide/landice/test_groups/ismip7_run.rst b/docs/users_guide/landice/test_groups/ismip7_run.rst new file mode 100644 index 0000000000..2eab9239c1 --- /dev/null +++ b/docs/users_guide/landice/test_groups/ismip7_run.rst @@ -0,0 +1,228 @@ +.. _landice_ismip7_run: + +ismip7_run +========== + +The ``landice/ismip7_run`` test group sets up one or more experiments from the +`ISMIP7 protocol `_ for both the +Antarctic Ice Sheet (AIS) and the Greenland Ice Sheet (GrIS). + +This functionality assumes the forcing files have already been generated using +the :ref:`landice_ismip7_forcing` test group and organized into the expected +directory layout. It creates a consistent set of run directories for the +requested experiments. Each experiment directory is self-contained with +namelists, streams, forcing symlinks, and a job script ready for submission. + +The test group includes two test cases: + +* ``ismip7_ais`` — Antarctic Ice Sheet experiments +* ``ismip7_gris`` — Greenland Ice Sheet experiments + +.. note:: + + This test group is not meant for automated running of experiments. + Expert knowledge is recommended for conducting the actual simulations. + Each experiment (step) should be submitted manually via its job script. + +.. _landice_ismip7_run_experiments: + +Experiment Matrix +----------------- + +The ISMIP7 core protocol defines 11 experiments per ice sheet: + +.. list-table:: + :header-rows: 1 + + * - Experiment + - Scenario + - Start + - End + - ESM + * - ``historical_CESM2-WACCM`` + - Historical + - ≥1850 + - 2014 + - CESM2-WACCM + * - ``historical_MRI-ESM2-0`` + - Historical + - ≥1850 + - 2014 + - MRI-ESM2-0 + * - ``ssp370_CESM2-WACCM`` + - SSP370 + - 2015 + - 2100 + - CESM2-WACCM + * - ``ssp370_MRI-ESM2-0`` + - SSP370 + - 2015 + - 2100 + - MRI-ESM2-0 + * - ``ssp126_CESM2-WACCM`` + - SSP126 + - 2015 + - 2300 + - CESM2-WACCM + * - ``ssp126_MRI-ESM2-0`` + - SSP126 + - 2015 + - 2300 + - MRI-ESM2-0 + * - ``ssp585_CESM2-WACCM`` + - SSP585 + - 2015 + - 2300 + - CESM2-WACCM + * - ``ssp585_MRI-ESM2-0`` + - SSP585 + - 2015 + - 2300 + - MRI-ESM2-0 + * - ``ctrl_CESM2-WACCM`` + - CTRL2015 + - 2015 + - 2300 + - CESM2-WACCM + * - ``ctrl_MRI-ESM2-0`` + - CTRL2015 + - 2015 + - 2300 + - MRI-ESM2-0 + * - ``ocx`` + - OCX + - 1990 + - 2025 + - (reanalysis) + +Unlike ISMIP6, ISMIP7 requires a **separate historical simulation per ESM**. +Projection experiments automatically symlink their restart file from the +corresponding ESM's historical run (e.g., +``ssp585_CESM2-WACCM`` → ``../historical_CESM2-WACCM/rst.2015-01-01.nc``). + +.. _landice_ismip7_run_usage: + +Usage +----- + +1. Process forcing data using :ref:`landice_ismip7_forcing`. + +2. Organize output into the expected directory layout:: + + {forcing_basepath}/ + ├── CESM2-WACCM_historical/ + │ ├── atmosphere/ + │ │ └── {mesh}_smb_CESM2-WACCM_historical_*.nc + │ └── ocean_thermal_forcing/ + │ └── {mesh}_thermal_forcing_CESM2-WACCM_historical_*.nc + ├── CESM2-WACCM_ssp585/ + │ ├── atmosphere/ + │ └── ocean_thermal_forcing/ + └── ... + +3. Create a user config file overriding the ``NotAvailable`` paths. + +4. Set up and run:: + + compass setup landice/ismip7_run/ismip7_ais -f my_ismip7_ais.cfg + # Then submit job scripts from individual experiment directories + +.. _landice_ismip7_run_config: + +config options +-------------- + +All config options should be reviewed and altered as needed. + +**AIS config** (``[ismip7_run_ais]``): + +.. code-block:: cfg + + [ismip7_run_ais] + + # Experiment list: "all", "historical", "projections", "ctrl", + # or comma-delimited experiment names + exp_list = all + + # Number of MPI tasks + ntasks = 128 + pio_stride = 128 + + # Base path to pre-processed forcing + forcing_basepath = NotAvailable + + # Initial condition and parameter files + init_cond_path = NotAvailable + melt_params_path = NotAvailable + region_mask_path = NotAvailable + + # Climatology files for CTRL2015 experiments + ctrl_tf_climatology_path = NotAvailable + ctrl_atm_climatology_path = NotAvailable + + # OCX forcing path + ocx_forcing_path = NotAvailable + + # Calving: restore or von_mises + calving_method = restore + von_mises_parameter_path = NotAvailable + + # Face melting + use_face_melting = false + + # Sea-level model coupling + sea_level_model = false + slm_input_ice = NotAvailable + slm_input_earth = NotAvailable + slm_earth_structure = prem_512.l60K2C.sum18p6.dum19p2.tz19p4.lm22 + slm_input_others = NotAvailable + nglv = 2048 + +**GrIS config** (``[ismip7_run_gris]``) is similar but without +sea-level model options and with ``calving_method = von_mises`` as default. + +.. _landice_ismip7_run_forcing_streams: + +Forcing Streams +--------------- + +ISMIP7 uses more forcing fields than ISMIP6, at mixed temporal resolutions: + +**Monthly forcing** (``input_interval = 0000-01-00_00:00:00``): + +* ``sfcMassBal`` — surface mass balance +* ``surfaceAirTemperature`` — surface air temperature +* ``ismip6Runoff`` — ice sheet runoff + +**Annual forcing** (``input_interval = 0001-00-00_00:00:00``): + +* ``sfcMassBalLapseRate`` — SMB elevation lapse rate +* ``surfaceAirTemperatureLapseRate`` — temperature lapse rate +* ``ismip6shelfMelt_3dThermalForcing`` (AIS) or + ``ismip6_2dThermalForcing`` (GrIS) — ocean thermal forcing + +**Static** (``input_interval = initial_only``): + +* ``ismip6shelfMelt_zOcean`` — ocean depth coordinates (AIS only) +* ``ismip6shelfMelt_deltaT``, ``ismip6shelfMelt_basin``, + ``ismip6shelfMelt_gamma0`` — melt parameterization coefficients + +For CTRL2015 experiments, all forcing intervals are set to +``initial_only`` (constant climate). + +.. _landice_ismip7_run_ais: + +ismip7_ais +---------- + +``landice/ismip7_run/ismip7_ais`` sets up AIS experiments with 3D ocean +thermal forcing (30 vertical layers) and optional sea-level model coupling. + +.. _landice_ismip7_run_gris: + +ismip7_gris +----------- + +``landice/ismip7_run/ismip7_gris`` sets up GrIS experiments with 2D +(depth-averaged) ocean thermal forcing. Sea-level model coupling is not +currently supported for GrIS. Von Mises calving is the default.