diff --git a/compass/landice/__init__.py b/compass/landice/__init__.py index a006dad5a8..2c010dc486 100644 --- a/compass/landice/__init__.py +++ b/compass/landice/__init__.py @@ -8,6 +8,7 @@ from compass.landice.tests.greenland import Greenland from compass.landice.tests.humboldt import Humboldt from compass.landice.tests.hydro_radial import HydroRadial +from compass.landice.tests.ismip6_forcing import Ismip6Forcing from compass.landice.tests.kangerlussuaq import Kangerlussuaq from compass.landice.tests.koge_bugt_s import KogeBugtS from compass.landice.tests.mismipplus import MISMIPplus @@ -34,6 +35,7 @@ def __init__(self): self.add_test_group(Greenland(mpas_core=self)) self.add_test_group(Humboldt(mpas_core=self)) self.add_test_group(HydroRadial(mpas_core=self)) + self.add_test_group(Ismip6Forcing(mpas_core=self)) self.add_test_group(Kangerlussuaq(mpas_core=self)) self.add_test_group(KogeBugtS(mpas_core=self)) self.add_test_group(MISMIPplus(mpas_core=self)) diff --git a/compass/landice/tests/ismip6_forcing/__init__.py b/compass/landice/tests/ismip6_forcing/__init__.py new file mode 100644 index 0000000000..0d0fed1383 --- /dev/null +++ b/compass/landice/tests/ismip6_forcing/__init__.py @@ -0,0 +1,26 @@ +from compass.testgroup import TestGroup +from compass.landice.tests.ismip6_forcing.atmosphere import Atmosphere +from compass.landice.tests.ismip6_forcing.ocean_thermal import OceanThermal +from compass.landice.tests.ismip6_forcing.ocean_basal import OceanBasal + + +class Ismip6Forcing(TestGroup): + """ + A test group for processing the ISMIP6 ocean and atmosphere forcing data + """ + + 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="ismip6_forcing") + + self.add_test_case(Atmosphere(test_group=self)) + self.add_test_case(OceanBasal(test_group=self)) + self.add_test_case(OceanThermal(test_group=self, process_obs=True)) + self.add_test_case(OceanThermal(test_group=self, process_obs=False)) diff --git a/compass/landice/tests/ismip6_forcing/atmosphere/__init__.py b/compass/landice/tests/ismip6_forcing/atmosphere/__init__.py new file mode 100644 index 0000000000..f8de7e3dad --- /dev/null +++ b/compass/landice/tests/ismip6_forcing/atmosphere/__init__.py @@ -0,0 +1,42 @@ +from compass.testcase import TestCase +from compass.landice.tests.ismip6_forcing.atmosphere.process_smb \ + import ProcessSMB +from compass.landice.tests.ismip6_forcing.atmosphere.process_smb_racmo \ + import ProcessSmbRacmo +from compass.landice.tests.ismip6_forcing.configure import configure as \ + configure_testgroup + + +class Atmosphere(TestCase): + """ + A test case for processing the ISMIP6 atmosphere forcing data. + The test case builds a mapping file for interpolation between the + ISMIP6 8km polarstereo grid and MALI mesh, regrids the forcing data + and rename the ISMIP6 variables to corresponding MALI variables. + """ + + def __init__(self, test_group): + """ + Create the test case + + Parameters + ---------- + test_group : compass.landice.tests.ismip6_forcing.Ismip6Forcing + The test group that this test case belongs to + """ + name = "atmosphere" + subdir = name + super().__init__(test_group=test_group, name=name, subdir=subdir) + + step = ProcessSmbRacmo(test_case=self) + self.add_step(step) + step = ProcessSMB(test_case=self) + self.add_step(step) + + def configure(self): + """ + Configures test case + """ + + configure_testgroup(config=self.config, + check_model_options=True) diff --git a/compass/landice/tests/ismip6_forcing/atmosphere/create_mapfile_smb.py b/compass/landice/tests/ismip6_forcing/atmosphere/create_mapfile_smb.py new file mode 100644 index 0000000000..f17c39329a --- /dev/null +++ b/compass/landice/tests/ismip6_forcing/atmosphere/create_mapfile_smb.py @@ -0,0 +1,204 @@ +import os +import numpy as np +import shutil +import netCDF4 +import xarray as xr +from mpas_tools.scrip.from_mpas import scrip_from_mpas +from mpas_tools.logging import check_call +from pyremap.descriptor import interp_extrap_corners_2d + + +def build_mapping_file(config, cores, logger, ismip6_grid_file, + mapping_file, mali_mesh_file=None, method_remap=None): + """ + Build a mapping file if it does not exist. + Mapping file is then used to remap the ismip6 source file in polarstero + coordinate to unstructured mali mesh + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for a ismip6 forcing test case + + cores : int + the number of cores for the ESMF_RegridWeightGen + + logger : logging.Logger + A logger for output from the step + + ismip6_grid_file : str + ismip6 grid file + + mapping_file : str + weights for interpolation from ismip6_grid_file to mali_mesh_file + + mali_mesh_file : str, optional + The MALI mesh file is used if mapping file does not exist + + method_remap : str, optional + Remapping method used in building a mapping file + """ + + if os.path.exists(mapping_file): + logger.info(f"Mapping file exists. Not building a new one.") + return + + # create the scrip files if mapping file does not exist + logger.info(f"Mapping file does not exist. Building one based on the" + f" input/output meshes") + logger.info(f"Creating temporary scrip files for source and " + f"destination grids...") + + if mali_mesh_file is None: + raise ValueError("Mapping file does not exist. To build one, Mali " + "mesh file with '-f' should be provided. " + "Type --help for info") + + # name temporary scrip files that will be used to build mapping file + source_grid_scripfile = "temp_source_scrip.nc" + mali_scripfile = "temp_mali_scrip.nc" + # this is the projection of ismip6 data for Antarctica + ismip6_projection = "ais-bedmap2" + + # create a scripfile for the atmosphere forcing data + create_atm_scrip(ismip6_grid_file, source_grid_scripfile) + + # create a MALI mesh scripfile + # make sure the mali mesh file uses the longitude convention of [0 2pi] + # make changes on a duplicated file to avoid making changes to the + # original mesh file + + mali_mesh_copy = f"{mali_mesh_file}_copy" + shutil.copy(mali_mesh_file, f"{mali_mesh_file}_copy") + + args = ["set_lat_lon_fields_in_planar_grid.py", + "--file", mali_mesh_copy, + "--proj", ismip6_projection] + + check_call(args, logger=logger) + + # create a MALI mesh scripfile if mapping file does not exist + scrip_from_mpas(mali_mesh_copy, mali_scripfile) + + # create a mapping file using ESMF weight gen + print(f"Creating a mapping file. Mapping method used: {method_remap}") + + if method_remap is None: + raise ValueError("Desired remapping option should be provided with " + "--method. Available options are 'bilinear'," + "'neareststod', 'conserve'.") + + parallel_executable = config.get('parallel', 'parallel_executable') + # split the parallel executable into constituents in case it includes flags + 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) + + # remove the temporary scripfiles once the mapping file is generated + logger.info(f"Removing the temporary mesh and scripfiles...") + os.remove(source_grid_scripfile) + os.remove(mali_scripfile) + os.remove(mali_mesh_copy) + + +def create_atm_scrip(source_grid_file, source_grid_scripfile): + """ + create a scripfile for the ismip6 atmospheric forcing data. + Note: the atmospheric forcing data do not have 'x' and 'y' coordinates and + only have dimensions of them. This function uses 'lat' and 'lon' + coordinates to generate a scripfile. + + Parameters + ---------- + source_grid_file : str + input smb grid file + + source_grid_scripfile : str + output scrip file of the input smb data + """ + + ds = xr.open_dataset(source_grid_file) + out_file = netCDF4.Dataset(source_grid_scripfile, 'w') + + if "rlon" in ds and "rlat" in ds: # this is for RACMO's rotated-pole grid + nx = ds.sizes["rlon"] + ny = ds.sizes["rlat"] + else: + nx = ds.sizes["x"] + ny = ds.sizes["y"] + units = "degrees" + + grid_size = nx * ny + + out_file.createDimension("grid_size", grid_size) + out_file.createDimension("grid_corners", 4) + out_file.createDimension("grid_rank", 2) + + # Variables + grid_center_lat = out_file.createVariable("grid_center_lat", "f8", + ("grid_size",)) + grid_center_lat.units = units + grid_center_lon = out_file.createVariable("grid_center_lon", "f8", + ("grid_size",)) + grid_center_lon.units = units + grid_corner_lat = out_file.createVariable("grid_corner_lat", "f8", + ("grid_size", "grid_corners")) + grid_corner_lat.units = units + grid_corner_lon = out_file.createVariable("grid_corner_lon", "f8", + ("grid_size", "grid_corners")) + grid_corner_lon.units = units + grid_imask = out_file.createVariable("grid_imask", "i4", ("grid_size",)) + grid_imask.units = "unitless" + out_file.createVariable("grid_dims", "i4", ("grid_rank",)) + + out_file.variables["grid_center_lat"][:] = ds.lat.values.flat + out_file.variables["grid_center_lon"][:] = ds.lon.values.flat + out_file.variables["grid_dims"][:] = [nx, ny] + out_file.variables["grid_imask"][:] = 1 + + if "lat_bnds" in ds and "lon_bnds" in ds: + lat_corner = ds.lat_bnds + if "time" in lat_corner.dims: + lat_corner = lat_corner.isel(time=0) + + lon_corner = ds.lon_bnds + if "time" in lon_corner.dims: + lon_corner = lon_corner.isel(time=0) + + lat_corner = lat_corner.values + lon_corner = lon_corner.values + else: + # this part is used for RACMO as it does not have lat_bnds & lon_bnds + lat_corner = _unwrap_corners(interp_extrap_corners_2d(ds.lat.values)) + lon_corner = _unwrap_corners(interp_extrap_corners_2d(ds.lon.values)) + + grid_corner_lat = lat_corner.reshape((grid_size, 4)) + grid_corner_lon = lon_corner.reshape((grid_size, 4)) + + out_file.variables["grid_corner_lat"][:] = grid_corner_lat + out_file.variables["grid_corner_lon"][:] = grid_corner_lon + + out_file.close() + + +def _unwrap_corners(in_field): + """ + Turn a 2D array of corners into an array of rectangular mesh elements + """ + out_field = np.zeros(((in_field.shape[0] - 1) * + (in_field.shape[1] - 1), 4)) + # corners are counterclockwise + out_field[:, 0] = in_field[0:-1, 0:-1].flat + out_field[:, 1] = in_field[0:-1, 1:].flat + out_field[:, 2] = in_field[1:, 1:].flat + out_field[:, 3] = in_field[1:, 0:-1].flat + + return out_field diff --git a/compass/landice/tests/ismip6_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip6_forcing/atmosphere/process_smb.py new file mode 100644 index 0000000000..c0cc95f9a8 --- /dev/null +++ b/compass/landice/tests/ismip6_forcing/atmosphere/process_smb.py @@ -0,0 +1,426 @@ +import os +import numpy as np +import shutil +import xarray as xr +from compass.landice.tests.ismip6_forcing.atmosphere.create_mapfile_smb \ + import build_mapping_file +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call +from compass.step import Step + + +class ProcessSMB(Step): + """ + A step for processing the ISMIP6 surface mass balance data + """ + + def __init__(self, test_case, input_file=None): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip6_forcing.atmosphere.Atmosphere + The test case this step belongs to + + input_file : file name of ismip6 forcing data processed by this step + """ + self.input_file = input_file + super().__init__(test_case=test_case, name="process_smb", ntasks=4, + min_tasks=1) + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config["ismip6_ais"] + base_path_ismip6 = section.get("base_path_ismip6") + base_path_mali = section.get("base_path_mali") + output_base_path = section.get("output_base_path") + mali_mesh_name = section.get("mali_mesh_name") + mali_mesh_file = section.get("mali_mesh_file") + period_endyear = section.get("period_endyear") + model = section.get("model") + scenario = section.get("scenario") + + self.add_input_file(filename=mali_mesh_file, + target=os.path.join(base_path_mali, + mali_mesh_file)) + + input_file_list = self._files[period_endyear][model][scenario] + for file in input_file_list: + self.add_input_file(filename=os.path.basename(file), + target=os.path.join(base_path_ismip6, + file)) + + output_file_esm = f"{mali_mesh_name}_SMB_{model}_{scenario}_" \ + f"{period_endyear}.nc" + self.add_output_file(filename=output_file_esm) + + # add processed racmo data as input as it is needed for smb correction + racmo_clim_file = f"{mali_mesh_name}_RACMO2.3p2_ANT27" \ + f"_smb_climatology_1995-2017.nc" + racmo_path = f"{output_base_path}/atmosphere_forcing/" \ + f"RACMO_climatology_1995-2017" + + self.add_input_file(filename=racmo_clim_file, + target=os.path.join(racmo_path, + racmo_clim_file)) + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + + section = config["ismip6_ais"] + mali_mesh_name = section.get("mali_mesh_name") + mali_mesh_file = section.get("mali_mesh_file") + period_endyear = section.get("period_endyear") + model = section.get("model") + scenario = section.get("scenario") + output_base_path = section.get("output_base_path") + + section = config["ismip6_ais_atmosphere"] + method_remap = section.get("method_remap") + + # define file names needed + # input racmo climotology file + racmo_clim_file = f"{mali_mesh_name}_RACMO2.3p2_ANT27" \ + f"_smb_climatology_1995-2017.nc" + racmo_path = f"{output_base_path}/atmosphere_forcing/" \ + f"RACMO_climatology_1995-2017" + # check if the processed racmo input file exists + racmo_clim_file_final = os.path.join(racmo_path, racmo_clim_file) + if not os.path.exists(racmo_clim_file_final): + raise ValueError("Processed RACMO data does not exist, " + "but it is required as an input file " + "to run this step. Please run `process_smb_racmo`" + "step prior to running this step by setting" + "the config option 'process_smb_racmo' to 'True'.") + + # temporary remapped climatology and anomaly files + clim_ismip6_temp = "clim_ismip6.nc" + remapped_clim_ismip6_temp = "remapped_clim_ismip6.nc" + remapped_anomaly_ismip6_temp = "remapped_anomaly_ismip6.nc" + # renamed remapped climatology and anomaly files (final outputs) + output_clim_ismip6 = f"{mali_mesh_name}_SMB_climatology_1995-2017_" \ + f"{model}_{scenario}.nc" + output_anomaly_ismip6 = f"{mali_mesh_name}_SMB_{model}_{scenario}_" \ + f"{period_endyear}.nc" + + # combine ismip6 forcing data covering different periods + # into a single file + input_file_list = self._files[period_endyear][model][scenario] + i = 0 + for file in input_file_list: + input_file_list[i] = os.path.basename(file) + i += 1 + + input_file_combined = xr.open_mfdataset(input_file_list, + concat_dim='time', + combine='nested', + engine='netcdf4') + combined_file_temp = "combined.nc" + write_netcdf(input_file_combined, combined_file_temp) + + # create smb climatology data over 1995-2017 + # take the time average over the period 1995-2017 + # note: make sure to have the correct time indexing for each + # smb anomaly files on which climatology is calculated. + logger.info(f"Calculating climatology for {model}_{scenario} forcing" + f"over 1995-2017") + args = [f"ncra", "-O", "-F", "-d", "time,1,23", + f"{combined_file_temp}", + f"{clim_ismip6_temp}"] + + check_call(args, logger=logger) + + # remap and rename the ismip6 smb climatology + logger.info("Remapping ismip6 climatology onto MALI mesh...") + self.remap_ismip6_smb_to_mali(clim_ismip6_temp, + remapped_clim_ismip6_temp, + mali_mesh_name, + mali_mesh_file, + method_remap) + + # rename the ismip6 variables to MALI variables + logger.info("Renaming the ismip6 variables to mali variable names...") + self.rename_ismip6_smb_to_mali_vars(remapped_clim_ismip6_temp, + output_clim_ismip6) + + # remap and rename ismip6 smb anomaly + logger.info(f"Remapping the {model}_{scenario} SMB anomaly onto " + f"MALI mesh") + self.remap_ismip6_smb_to_mali(combined_file_temp, + remapped_anomaly_ismip6_temp, + mali_mesh_name, + mali_mesh_file, + method_remap) + + # rename the ismip6 variables to MALI variables + logger.info("Renaming the ismip6 variables to mali variable names...") + self.rename_ismip6_smb_to_mali_vars(remapped_anomaly_ismip6_temp, + output_anomaly_ismip6) + + # correct the SMB anomaly field with mali base SMB field + logger.info("Correcting the SMB anomaly field for the base SMB " + "climatology 1995-2017...") + self.correct_smb_anomaly_for_climatology(racmo_clim_file, + output_clim_ismip6, + output_anomaly_ismip6) + + logger.info("Processing done successfully. " + "Removing the temporary files...") + # remove the temporary remapped and combined files + os.remove(remapped_clim_ismip6_temp) + os.remove(remapped_anomaly_ismip6_temp) + os.remove(combined_file_temp) + os.remove(clim_ismip6_temp) + os.remove(output_clim_ismip6) + + # place the output file in appropriate directory + output_path = f"{output_base_path}/atmosphere_forcing/" \ + f"{model}_{scenario}/1995-{period_endyear}" + if not os.path.exists(output_path): + print("Creating a new directory for the output data...") + os.makedirs(output_path) + + src = os.path.join(os.getcwd(), output_anomaly_ismip6) + dst = os.path.join(output_path, output_anomaly_ismip6) + shutil.copy(src, dst) + + logger.info(f"!---Done processing the file---!") + + def remap_ismip6_smb_to_mali(self, input_file, output_file, mali_mesh_name, + mali_mesh_file, method_remap): + """ + Remap the input ismip6 smb forcing data onto mali mesh + + Parameters + ---------- + input_file: str + input smb data on its native polarstereo 8km grid + + output_file : str + smb data remapped on mali mesh + + mali_mesh_name : str + name of the mali mesh used to name mapping files + + mali_mesh_file : str, optional + The MALI mesh file if mapping file does not exist + + method_remap : str, optional + Remapping method used in building a mapping file + """ + mapping_file = f"map_ismip6_8km_to_" \ + f"{mali_mesh_name}_{method_remap}.nc" + + # check if mapfile exists + if not os.path.exists(mapping_file): + # build a mapping file if it doesn't already exist + self.logger.info(f"Creating a mapping file. " + f"Mapping method used: {method_remap}") + build_mapping_file(self.config, self.ntasks, self.logger, + input_file, mapping_file, mali_mesh_file, + method_remap) + else: + self.logger.info("Mapping file exists. " + "Remapping the input data...") + + # remap the input data + args = ["ncremap", + "-i", input_file, + "-o", output_file, + "-m", mapping_file, + "-v", "smb_anomaly"] + + check_call(args, logger=self.logger) + + def rename_ismip6_smb_to_mali_vars(self, remapped_file_temp, output_file): + """ + Rename variables in the remapped source input data + to the ones that MALI uses. + + Parameters + ---------- + remapped_file_temp : str + temporary ismip6 data remapped on mali mesh + + output_file : str + remapped ismip6 data renamed on mali mesh + """ + + # open dataset in 20 years chunk + ds = xr.open_dataset(remapped_file_temp, chunks=dict(time=20), + engine="netcdf4") + + # build dictionary for ismip6 variables that MALI takes in + ismip6_to_mali_dims = dict( + time="Time", + ncol="nCells") + ds = ds.rename(ismip6_to_mali_dims) + + ismip6_to_mali_vars = dict( + smb_anomaly="sfcMassBal") + ds = ds.rename(ismip6_to_mali_vars) + + # add xtime variable + xtime = [] + for t_index in range(ds.sizes["Time"]): + date = ds.Time[t_index] + # forcing files do not all match even years, so round up the years + # pandas round function does not work for years, so do it manually + yr = date.dt.year.values + mo = date.dt.month.values + dy = date.dt.day.values + dec_yr = np.around(yr + (30 * (mo - 1) + dy) / 365.0) + date = f"{dec_yr.astype(int)}-01-01_00:00:00".ljust(64) + xtime.append(date) + + ds["xtime"] = ("Time", xtime) + ds["xtime"] = ds.xtime.astype('S') + + # drop unnecessary variables + ds = ds.drop_vars(["lon", "lon_vertices", "lat", "lat_vertices", + "area"]) + + # write to a new netCDF file + write_netcdf(ds, output_file) + ds.close() + + def correct_smb_anomaly_for_climatology(self, + racmo_clim_file, + output_clim_ismip6_file, + output_file_final): + + """ + Apply the MALI base SMB to the ismip6 SMB anomaly field + + Parameters + ---------- + racmo_clim_file : str + RACMO climatology file (1995-2017) + + output_clim_ismip6_file : str + remapped and renamed ismip6 climatology file + + output_file_final : str + climatology-corrected, final ismip6 smb anomaly file + """ + + ds = xr.open_dataset(output_file_final) + ds_racmo_clim = xr.open_dataset(racmo_clim_file) + ds_ismip6_clim = xr.open_dataset(output_clim_ismip6_file) + + # calculate the climatology correction + corr_clim = (ds_racmo_clim["sfcMassBal"].isel(Time=0) - + ds_ismip6_clim["sfcMassBal"].isel(Time=0)) + + # correct the ismip6 smb anomaly + ds["sfcMassBal"] = ds["sfcMassBal"] + corr_clim + + # write metadata + ds["sfcMassBal"].attrs = {"long_name" : "surface mass balance", + "units" : "kg m-2 s-1", + "coordinates" : "lat lon"} + + # write to a new netCDF file + write_netcdf(ds, output_file_final) + ds.close() + + # create a nested dictionary for the ISMIP6 original forcing files including relative path + # Note: these files needed to be explicitly listed because of inconsistencies that are + # present in file naming conventions in the ISMIP6 source dataset. + _files = { + "2100": { + "CCSM4": { + "RCP85": [ + "AIS/Atmosphere_Forcing/ccsm4_rcp8.5/Regridded_8km/CCSM4_8km_anomaly_1995-2100.nc"] + }, + "CESM2": { + "SSP585v1": [ + "AIS/Atmosphere_Forcing/CESM2_ssp585/Regridded_8km/CESM2_anomaly_ssp585_1995-2100_8km_v1.nc"], + "SSP585v2": [ + "AIS/Atmosphere_Forcing/CESM2_ssp585/Regridded_8km/CESM2_anomaly_ssp585_1995-2100_8km_v2.nc"] + }, + "CNRM_CM6": { + "SSP126": [ + "AIS/Atmosphere_Forcing/CNRM_CM6_ssp126/Regridded_8km/CNRM-CM6-1_anomaly_ssp126_1995-2100_8km_ISMIP6.nc"], + "SSP585": [ + "AIS/Atmosphere_Forcing/CNRM_CM6_ssp585/Regridded_8km/CNRM-CM6-1_anomaly_ssp585_1995-2100_8km_ISMIP6.nc"] + }, + "CNRM_ESM2": { + "SSP585": [ + "AIS/Atmosphere_Forcing/CNRM_ESM2_ssp585/Regridded_8km/CNRM-ESM2-1_anomaly_ssp585_1995-2100_8km_ISMIP6.nc"] + }, + "CSIRO-Mk3-6-0": { + "RCP85": [ + "AIS/Atmosphere_Forcing/CSIRO-Mk3-6-0_rcp85/Regridded_8km/CSIRO-Mk3-6-0_8km_anomaly_rcp85_1995-2100.nc"] + }, + "HadGEM2-ES": { + "RCP85": [ + "AIS/Atmosphere_Forcing/HadGEM2-ES_rcp85/Regridded_8km/HadGEM2-ES_8km_anomaly_rcp85_1995-2100.nc"] + }, + "IPSL-CM5A-MR": { + "RCP26": [ + "AIS/Atmosphere_Forcing/IPSL-CM5A-MR_rcp26/Regridded_8km/IPSL-CM5A-MR_8km_anomaly_rcp26_1995-2100.nc"], + "RCP85": [ + "AIS/Atmosphere_Forcing/IPSL-CM5A-MR_rcp85/Regridded_8km/IPSL-CM5A-MR_8km_anomaly_rcp85_1995-2100.nc"] + }, + "MIROC-ESM-CHEM": { + "RCP85": [ + "AIS/Atmosphere_Forcing/miroc-esm-chem_rcp8.5/Regridded_8km/MIROC-ESM-CHEM_8km_anomaly_1995-2100.nc"] + }, + "NorESM1-M": { + "RCP26": [ + "AIS/Atmosphere_Forcing/noresm1-m_rcp2.6/Regridded_8km/NorESM-M_8km_anomaly_rcp26_1995-2100.nc"], + "RCP85": [ + "AIS/Atmosphere_Forcing/noresm1-m_rcp8.5/Regridded_8km/NorESM-M_8km_anomaly_1995-2100.nc"] + }, + "UKESM1-0-LL": { + "SSP585": [ + "AIS/Atmosphere_Forcing/UKESM1-0-LL/Regridded_8km/UKESM1-0-LL_anomaly_ssp585_1995-2100_8km.nc"] + } + }, + "2300": { + "CCSM4": { + "RCP85": [ + "AIS/Atmospheric_forcing/CCSM4_RCP85/Regridded_08km/CCSM4_8km_anomaly_1995-2100.nc", + "AIS/Atmospheric_forcing/CCSM4_RCP85/Regridded_08km/CCSM4_8km_anomaly_2101-2300.nc"] + }, + "CESM2-WACCM": { + "SSP585": [ + "AIS/Atmospheric_forcing/CESM2-WACCM_ssp585/Regridded_8km/CESM2-WACCM_8km_anomaly_ssp585_1995-2100.nc", + "AIS/Atmospheric_forcing/CESM2-WACCM_ssp585/Regridded_8km/CESM2-WACCM_8km_anomaly_ssp585_2101-2299.nc"], + "SSP585-repeat": [ + "AIS/Atmospheric_forcing/CESM2-WACCM_ssp585-repeat/Regridded_8km/CESM2-WACCM_8km_anomaly_ssp585_1995-2300_v2.nc"] + }, + "HadGEM2-ES": { + "RCP85": [ + "AIS/Atmospheric_forcing/HadGEM2-ES_RCP85/Regridded_8km/HadGEM2-ES_8km_anomaly_rcp85_1995-2100.nc", + "AIS/Atmospheric_forcing/HadGEM2-ES_RCP85/Regridded_8km/HadGEM2-ES_8km_anomaly_rcp85_2101-2299.nc"], + "RCP85-repeat": [ + "AIS/Atmospheric_forcing/HadGEM2-ES-RCP85-repeat/Regridded_8km/HadGEM2-ES_8km_anomaly_rcp85_1995-2300_v2.nc"] + }, + "NorESM1-M": { + "RCP26-repeat": [ + "AIS/Atmospheric_forcing/NorESM1-M_RCP26-repeat/Regridded_8km/NorESM1-M_8km_anomaly_rcp26_1995-2300_v2.nc"], + "RCP85-repeat": [ + "AIS/Atmospheric_forcing/NorESM1-M_RCP85-repeat/Regridded_8km/NorESM1-M_8km_anomaly_1995-2300_v2.nc"] + }, + "UKESM1-0-LL": { + "SSP126": [ + "AIS/Atmospheric_forcing/UKESM1-0-LL_ssp126/Regridded_8km/UKESM1-0-LL_8km_anomaly_ssp126_1995-2100.nc", + "AIS/Atmospheric_forcing/UKESM1-0-LL_ssp126/Regridded_8km/UKESM1-0-LL_8km_anomaly_ssp126_2101-2300.nc"], + "SSP585": [ + "AIS/Atmospheric_forcing/UKESM1-0-LL_ssp585/Regridding_8km/UKESM1-0-LL_8km_anomaly_ssp585_1995-2100.nc", + "AIS/Atmospheric_forcing/UKESM1-0-LL_ssp585/Regridding_8km/UKESM1-0-LL_8km_anomaly_ssp585_2101-2300.nc"], + "SSP585-repeat": [ + "AIS/Atmospheric_forcing/UKESM1-0-LL_ssp585-repeat/Regridding_8km/UKESM1-0-LL_8km_anomaly_ssp585_1995-2300_v2.nc"] + } + } + } diff --git a/compass/landice/tests/ismip6_forcing/atmosphere/process_smb_racmo.py b/compass/landice/tests/ismip6_forcing/atmosphere/process_smb_racmo.py new file mode 100644 index 0000000000..5ffdd35da7 --- /dev/null +++ b/compass/landice/tests/ismip6_forcing/atmosphere/process_smb_racmo.py @@ -0,0 +1,258 @@ +import os +import shutil +import xarray as xr +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call +from compass.step import Step +from compass.landice.tests.ismip6_forcing.atmosphere.create_mapfile_smb \ + import build_mapping_file + + +class ProcessSmbRacmo(Step): + """ + A step for processing the regional RACMO surface mass balance data + """ + + def __init__(self, test_case, input_file=None): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip6_forcing.atmosphere.Atmosphere + The test case this step belongs to + + input_file : file name of ismip6 forcing data processed by this step + """ + self.input_file = input_file + super().__init__(test_case=test_case, name='process_smb_racmo', + ntasks=4, min_tasks=1) + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config["ismip6_ais"] + base_path_mali = section.get("base_path_mali") + mali_mesh_name = section.get("mali_mesh_name") + mali_mesh_file = section.get("mali_mesh_file") + + section = config["ismip6_ais_atmosphere"] + process_smb_racmo = section.getboolean("process_smb_racmo") + + self.add_input_file(filename=mali_mesh_file, + target=os.path.join(base_path_mali, + mali_mesh_file)) + + if process_smb_racmo: + self.add_input_file( + target="RACMO2.3p2_ANT27_smb_yearly_1979_2018.nc", + database="RACMO2.3p2_ANT27_SMB_yearly_1979_2018") + + output_file = f"{mali_mesh_name}_RACMO2.3p2_ANT27" \ + f"_smb_climatology_1995-2017.nc" + self.add_output_file(filename=output_file) + else: + print(f"\n'Warning: process_smb_racmo' is set to 'False'." + f" This step will not run unless set 'True' in the" + f" config file.\n") + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + + section = config["ismip6_ais_atmosphere"] + process_smb_racmo = section.getboolean("process_smb_racmo") + if not process_smb_racmo: + # we don't want to run this step + return + + section = config["ismip6_ais"] + mali_mesh_name = section.get("mali_mesh_name") + mali_mesh_file = section.get("mali_mesh_file") + output_base_path = section.get("output_base_path") + + section = config["ismip6_ais_atmosphere"] + method_remap = section.get("method_remap") + + racmo_file_temp1 = "RACMO2.3p2_smb_climatology_1995_2017.nc" + racmo_file_temp2 = "RACMO2.3p2_smb_climatology_1995_2017_" \ + "correct_unit.nc" + output_file = f"{mali_mesh_name}_RACMO2.3p2_ANT27" \ + f"_smb_climatology_1995-2017.nc" + output_path = f"{output_base_path}/atmosphere_forcing/" \ + f"RACMO_climatology_1995-2017" + output_path_final = os.path.join(output_base_path, output_path) + + if os.path.exists(os.path.join(output_path_final, output_file)): + logger.info(f"Processed RACMO SMB data already exists in the " + f"output path {output_base_path}. " + f"Not processing the source RACMO data...") + return + + input_file = self.inputs[1] + # take the time average over the period 1995-2017 + args = ["ncra", "-O", "-F", "-d", "time,17,39", + input_file, + racmo_file_temp1] + + check_call(args, logger=logger) + + # interpolate the racmo smb data + remapped_file_temp = "remapped.nc" # temporary file name + + # call the function that reads in, remap and rename the file. + logger.info("Calling the remapping function...") + self.remap_source_smb_to_mali(racmo_file_temp1, + remapped_file_temp, + mali_mesh_name, + mali_mesh_file, + method_remap) + + # perform algebraic operation on the source data in unit of kg/m^2 + # to be in unit of kg/m^2/s + args = ["ncap2", "-O", "-v", "-s", + "sfcMassBal=smb/(60*60*24*365)", + remapped_file_temp, + racmo_file_temp2] + + check_call(args, logger=logger) + + # change the unit attribute to kg/m^2/s + args = ["ncatted", "-O", "-a", + "units,sfcMassBal,m,c,'kg m-2 s-1'", + racmo_file_temp2] + + check_call(args, logger=logger) + + # call the function that renames the ismip6 variables to MALI variables + logger.info("Renaming source variables to mali variable names...") + + self.rename_source_smb_to_mali_vars(racmo_file_temp2, output_file) + + logger.info("Processing done successfully. " + "Removing the temporary files...") + # remove the temporary remapped and combined files + os.remove(remapped_file_temp) + os.remove(racmo_file_temp1) + os.remove(racmo_file_temp2) + + # place the output file in appropriate directory + output_path = f"{output_base_path}/atmosphere_forcing/" \ + f"RACMO_climatology_1995-2017" + if not os.path.exists(output_path): + logger.info("Creating a new directory for the output data") + os.makedirs(output_path) + + src = os.path.join(os.getcwd(), output_file) + dst = os.path.join(output_path, output_file) + shutil.copy(src, dst) + + logger.info(f"!---Done processing the file---!") + + def remap_source_smb_to_mali(self, input_file, output_file, mali_mesh_name, + mali_mesh_file, method_remap): + """ + Remap the input ismip6 smb forcing data onto mali mesh + + Parameters + ---------- + input_file: str + input racmo smb data on its native rotated pole grid + + output_file : str + smb data remapped on mali mesh + + mali_mesh_name : str + name of the mali mesh used to name mapping files + + mali_mesh_file : str, optional + The MALI mesh file if mapping file does not exist + + method_remap : str, optional + Remapping method used in building a mapping file + """ + mapping_file = f"map_racmo_24km_to_" \ + f"{mali_mesh_name}_{method_remap}.nc" + + # check if a mapfile exists + if not os.path.exists(mapping_file): + # build a mapping file if it doesn't already exist + self.logger.info(f"Creating a mapping file. " + f"Mapping method used: {method_remap}") + build_mapping_file(self.config, self.ntasks, self.logger, + input_file, mapping_file, mali_mesh_file, + method_remap) + else: + self.logger.info("Mapping file exists. " + "Remapping the input data...") + + args = ["ncremap", + "-i", input_file, + "-o", output_file, + "-m", mapping_file, + "-v", "smb"] + + check_call(args, logger=self.logger) + + def rename_source_smb_to_mali_vars(self, remapped_file_temp, output_file): + """ + Rename variables in the remapped source input data + to the ones that MALI uses. + + Parameters + ---------- + remapped_file_temp : str + temporary ismip6 data remapped on mali mesh + + output_file : str + remapped ismip6 data renamed on mali mesh + """ + + # open dataset in 20 years chunk + ds = xr.open_dataset(remapped_file_temp, chunks=dict(time=20), + engine="netcdf4") + + # build dictionary for ismip6 variables that MALI takes in + ismip6_to_mali_dims = dict( + time="Time", + ncol="nCells") + ds = ds.rename(ismip6_to_mali_dims) + + # drop unnecessary variables + ds = ds.drop_vars(["height"]) + + # squeeze unnecessary coordinate variable + ds["sfcMassBal"] = ds["sfcMassBal"].squeeze(dim="height") + + # write to a new netCDF file + write_netcdf(ds, output_file) + ds.close() + + def correct_smb_anomaly_for_base_smb(self, output_file, mali_mesh_file): + """ + Apply the MALI base SMB to the ismip6 SMB anomaly field + + Parameters + ---------- + output_file : str + remapped ismip6 data renamed on mali mesh + + mali_mesh_file : str + initialized MALI mesh file in which the base SMB field exists + """ + + ds = xr.open_dataset(output_file) + ds_base = xr.open_dataset(mali_mesh_file) + # get the first time index + ref_smb = ds_base["sfcMassBal"].isel(Time=0) + # correct for the reference smb + ds["sfcMassBal"] = ds["sfcMassBal"] + ref_smb + + # write to a new netCDF file + write_netcdf(ds, output_file) + ds.close() diff --git a/compass/landice/tests/ismip6_forcing/configure.py b/compass/landice/tests/ismip6_forcing/configure.py new file mode 100644 index 0000000000..1f9ec16570 --- /dev/null +++ b/compass/landice/tests/ismip6_forcing/configure.py @@ -0,0 +1,26 @@ +def configure(config, check_model_options): + """ + A shared function for configuring options for all ismip6 forcing + test cases + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for a ismip6 forcing test case + + check_model_options : bool + Whether we check ``model``, ``scenario``, ``period_endyear`` + """ + + section = "ismip6_ais" + options = ["base_path_ismip6", "base_path_mali", "mali_mesh_name", + "mali_mesh_file", "output_base_path"] + if check_model_options: + options = options + ["model", "scenario", "period_endyear"] + + 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/ismip6_forcing/create_mapfile.py b/compass/landice/tests/ismip6_forcing/create_mapfile.py new file mode 100644 index 0000000000..954f5a719b --- /dev/null +++ b/compass/landice/tests/ismip6_forcing/create_mapfile.py @@ -0,0 +1,111 @@ +import os +import shutil +import subprocess +from mpas_tools.scrip.from_mpas import scrip_from_mpas +from mpas_tools.logging import check_call + + +def build_mapping_file(config, cores, logger, ismip6_grid_file, + mapping_file, mali_mesh_file=None, + method_remap=None): + """ + Build a mapping file if it does not exist. + Mapping file is then used to remap the ismip6 source file in polarstero + coordinate to unstructured mali mesh + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for a ismip6 forcing test case + + cores : int + the number of cores for the ESMF_RegridWeightGen + + logger : logging.Logger + A logger for output from the step + + ismip6_grid_file : str + ismip6 grid file + + mapping_file : str + weights for interpolation from ismip6_grid_file to mali_mesh_file + + mali_mesh_file : str, optional + The MALI mesh file if mapping file does not exist + + method_remap : str, optional + Remapping method used in building a mapping file + """ + + if os.path.exists(mapping_file): + logger.info(f"Mapping file exists. Not building a new one.") + return + + if mali_mesh_file is None: + raise ValueError("Mapping file does not exist. To build one, Mali " + "mesh file with '-f' should be provided. " + "Type --help for info") + + ismip6_scripfile = "temp_ismip6_8km_scrip.nc" + mali_scripfile = "temp_mali_scrip.nc" + ismip6_projection = "ais-bedmap2" + + # create the ismip6 scripfile if mapping file does not exist + # this is the projection of ismip6 data for Antarctica + logger.info(f"Mapping file does not exist. Building one based on " + f"the input/ouptut meshes") + logger.info(f"Creating temporary scripfiles " + f"for ismip6 grid and mali mesh...") + + args = ["create_SCRIP_file_from_planar_rectangular_grid.py", + "--input", ismip6_grid_file, + "--scrip", ismip6_scripfile, + "--proj", ismip6_projection, + "--rank", "2"] + + check_call(args, logger=logger) + + # create a MALI mesh scripfile if mapping file does not exist + # make sure the mali mesh file uses the longitude convention of [0 2pi] + # make changes on a duplicated file to avoid making changes to the + # original mesh file + + mali_mesh_copy = f"{mali_mesh_file}_copy" + shutil.copy(mali_mesh_file, f"{mali_mesh_file}_copy") + + args = ["set_lat_lon_fields_in_planar_grid.py", + "--file", mali_mesh_copy, + "--proj", ismip6_projection] + + check_call(args, logger=logger) + + scrip_from_mpas(mali_mesh_file, mali_scripfile) + + # create a mapping file using ESMF weight gen + logger.info(f"Creating a mapping file... " + f"Mapping method used: {method_remap}") + + if method_remap is None: + raise ValueError("Desired remapping option should be provided with " + "--method. Available options are 'bilinear'," + "'neareststod', 'conserve'.") + + parallel_executable = config.get("parallel", "parallel_executable") + # split the parallel executable into constituents in case it includes flags + args = parallel_executable.split(' ') + args.extend(["-n", f"{cores}", + "ESMF_RegridWeightGen", + "-s", ismip6_scripfile, + "-d", mali_scripfile, + "-w", mapping_file, + "-m", method_remap, + "-i", "-64bit_offset", + "--dst_regional", "--src_regional"]) + + check_call(args, logger) + + # remove the temporary scripfiles once the mapping file is generated + logger.info(f"Removing the temporary mesh and scripfiles...") + os.remove(ismip6_scripfile) + os.remove(mali_scripfile) + os.remove(mali_mesh_copy) diff --git a/compass/landice/tests/ismip6_forcing/ismip6_forcing.cfg b/compass/landice/tests/ismip6_forcing/ismip6_forcing.cfg new file mode 100644 index 0000000000..6b88abcd80 --- /dev/null +++ b/compass/landice/tests/ismip6_forcing/ismip6_forcing.cfg @@ -0,0 +1,45 @@ +# config options for ismip6 antarctic ice sheet data set +[ismip6_ais] + +# Base path to the input ismip6 ocean and smb forcing files. User has to supply. +base_path_ismip6 = NotAvailable + +# Base path to the the MALI mesh. User has to supply. +base_path_mali = NotAvailable + +# Forcing end year of the ISMIP6 data. User has to supply. +# Available end years are 2100 and 2300. +period_endyear = NotAvailable + +# Base path to which output forcing files are saved. +output_base_path = NotAvailable + +# Name of climate model name used to generate ISMIP6 forcing data. User has to supply. +# Available model names are the following: CCSM4, CESM2-WACCM, CSIRO-Mk3-6-0, HadGEM2-ES, NorESM1-M, UKESM1-0-LL +model = NotAvailable + +# Scenarios used by climate model. User has to supply. +# Available scenarios are the following: RCP26, RCP85, SSP126, SSP585 +scenario = NotAvailable + +# name of the mali mesh. User has to supply. Note: It is used to name mapping files +# (e,g. 'map_ismip6_8km_to_{mali_mesh_name}_{method_remap}.nc'). +mali_mesh_name = NotAvailable + +# MALI mesh file to be used to build mapping file (e.g.Antarctic_8to80km_20220407.nc). User has to supply. +mali_mesh_file = NotAvailable + +# config options for ismip6 antarctic ice sheet SMB forcing data test cases +[ismip6_ais_atmosphere] + +# Remapping method used in building a mapping file. Options include: bilinear, neareststod, conserve +method_remap = bilinear + +# Set True to process RACMO modern climatology +process_smb_racmo = True + +# config options for ismip6 ocean thermal forcing data test cases +[ismip6_ais_ocean_thermal] + +# Remapping method used in building a mapping file. Options include: bilinear, neareststod, conserve +method_remap = bilinear diff --git a/compass/landice/tests/ismip6_forcing/ocean_basal/__init__.py b/compass/landice/tests/ismip6_forcing/ocean_basal/__init__.py new file mode 100644 index 0000000000..f4732ba263 --- /dev/null +++ b/compass/landice/tests/ismip6_forcing/ocean_basal/__init__.py @@ -0,0 +1,35 @@ +from compass.testcase import TestCase +from compass.landice.tests.ismip6_forcing.ocean_basal.process_basal_melt \ + import ProcessBasalMelt +from compass.landice.tests.ismip6_forcing.configure import configure as \ + configure_testgroup + + +class OceanBasal(TestCase): + """ + A test case for processing the ISMIP6 ocean basalmelt param. coeff. data. + The test case builds a mapping file for interpolation between the + ISMIP6 8km polarstereo grid and MALI mesh, regrids the forcing data + and renames the ISMIP6 variables to corresponding MALI variables. + """ + + def __init__(self, test_group): + """ + Create the test case + + Parameters + ---------- + test_group : compass.landice.tests.ismip6_forcing.Ismip6Forcing + The test group that this test case belongs to + """ + name = "ocean_basal" + super().__init__(test_group=test_group, name=name) + + step = ProcessBasalMelt(test_case=self) + self.add_step(step) + + def configure(self): + """ + Configures test case + """ + configure_testgroup(config=self.config, check_model_options=False) diff --git a/compass/landice/tests/ismip6_forcing/ocean_basal/process_basal_melt.py b/compass/landice/tests/ismip6_forcing/ocean_basal/process_basal_melt.py new file mode 100644 index 0000000000..f3ba4eec0c --- /dev/null +++ b/compass/landice/tests/ismip6_forcing/ocean_basal/process_basal_melt.py @@ -0,0 +1,269 @@ +import os +import shutil +import xarray as xr +from compass.landice.tests.ismip6_forcing.create_mapfile \ + import build_mapping_file +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call +from compass.step import Step + + +class ProcessBasalMelt(Step): + """ + A step for processing (combine, remap and rename) the ISMIP6 basalmelt + data + """ + + def __init__(self, test_case): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip6_forcing.ocean_basal.OceanBasal + The test case this step belongs to + """ + super().__init__(test_case=test_case, name="process_basal_melt", + ntasks=4, min_tasks=1) + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config["ismip6_ais"] + base_path_ismip6 = section.get("base_path_ismip6") + base_path_mali = section.get("base_path_mali") + mali_mesh_name = section.get("mali_mesh_name") + mali_mesh_file = section.get("mali_mesh_file") + period_endyear = section.get("period_endyear") + + input_file_basin = self._file_basin[period_endyear] + input_files_coeff = self._files_coeff[period_endyear] + self.add_input_file(filename=mali_mesh_file, + target=os.path.join(base_path_mali, + mali_mesh_file)) + self.add_input_file(filename=os.path.basename(input_file_basin), + target=os.path.join(base_path_ismip6, + input_file_basin)) + + for file in input_files_coeff: + self.add_input_file(filename=os.path.basename(file), + target=os.path.join(base_path_ismip6, file)) + self.add_output_file(filename=f"{mali_mesh_name}_basin_and_" + f"{os.path.basename(file)}") + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + section = config["ismip6_ais"] + mali_mesh_name = section.get("mali_mesh_name") + mali_mesh_file = section.get("mali_mesh_file") + period_endyear = section.get("period_endyear") + output_base_path = section.get("output_base_path") + # we always want neareststod for the remapping method because we want + # a single value per basin + method_remap = "neareststod" + + # combine, interpolate and rename the basin file and deltaT0_gamma0 + # ismip6 input files + + # call the function that combines data + logger.info(f"Combining the ismip6 input files") + input_file_basin = self._file_basin[period_endyear] + input_file_list = self._files_coeff[period_endyear] + + for i, file in enumerate(input_file_list): + logger.info(f"!---Start processing the current file---!") + logger.info(f"processing the input file " + f"'{os.path.basename(file)}'") + # temporary file names. Note: The counter 'i' seems to be necessary + # for 'ncremap' to correctly interpolate the gamma0 values for an + # unknown reason. + combined_file_temp = f"combined_{i}.nc" + remapped_file_temp = f"remapped_{i}.nc" + + self.combine_ismip6_inputfiles(os.path.basename(input_file_basin), + os.path.basename(file), + combined_file_temp) + + # remap the input forcing file. + logger.info(f"Calling the remapping function...") + self.remap_ismip6_basal_melt_to_mali_vars(combined_file_temp, + remapped_file_temp, + mali_mesh_name, + mali_mesh_file, + method_remap) + + output_file = f"{mali_mesh_name}_basin_and_" \ + f"{os.path.basename(file)}" + + # rename the ismip6 variables to MALI variables + logger.info(f"Renaming the ismip6 variables to " + f"mali variable names...") + self.rename_ismip6_basal_melt_to_mali_vars(remapped_file_temp, + output_file) + + logger.info(f"Remapping and renaming process done successfully. " + f"Removing the temporary files...") + + # remove the temporary combined file + os.remove(combined_file_temp) + os.remove(remapped_file_temp) + + # place the output file in appropriate director + output_path = f"{output_base_path}/basal_melt/parameterizations/" + if not os.path.exists(output_path): + logger.info("Creating a new directory for the output data...") + os.makedirs(output_path) + + src = os.path.join(os.getcwd(), output_file) + dst = os.path.join(output_path, output_file) + shutil.copy(src, dst) + + logger.info(f"!---Done processing the current file---!") + logger.info(f"") + logger.info(f"") + + def combine_ismip6_inputfiles(self, basin_file, coeff_gamma0_deltat_file, + combined_file_temp): + """ + Combine ismip6 input files before regridding onto mali mesh + + Parameters + ---------- + basin_file : str + imbie2 basin numbers in ismip6 grid + + coeff_gamma0_deltat_file : str + uniform melt rate coefficient (gamma0) and temperature + correction per basin + + combined_file_temp : str + temporary output file that has all the variables combined + """ + + ds_basin = xr.open_dataset(basin_file, engine="netcdf4") + ds = xr.open_dataset(coeff_gamma0_deltat_file, engine="netcdf4") + + ds["ismip6shelfMelt_basin"] = ds_basin.basinNumber + write_netcdf(ds, combined_file_temp) + + def remap_ismip6_basal_melt_to_mali_vars(self, input_file, output_file, + mali_mesh_name, mali_mesh_file, + method_remap): + """ + Remap the input ismip6 basal melt data onto mali mesh + + Parameters + ---------- + input_file: str + temporary output file that has all the variables combined + combined_file_temp generated in the above function + + output_file : str + ismip6 data remapped on mali mesh + + mali_mesh_name : str + name of the mali mesh used to name mapping files + + mali_mesh_file : str, optional + The MALI mesh file if mapping file does not exist + + method_remap : str, optional + Remapping method used in building a mapping file + """ + + # check if mapfile exists + mapping_file = f"map_ismip6_8km_to_{mali_mesh_name}_{method_remap}.nc" + + if not os.path.exists(mapping_file): + # build a mapping file if it doesn't already exist + build_mapping_file(self.config, self.ntasks, self.logger, + input_file, mapping_file, mali_mesh_file, + method_remap) + else: + self.logger.info(f"Mapping file exists. " + f"Remapping the input data...") + + # remap the input data + args = ["ncremap", + "-i", input_file, + "-o", output_file, + "-m", mapping_file] + + check_call(args, logger=self.logger) + + def rename_ismip6_basal_melt_to_mali_vars(self, remapped_file_temp, + output_file): + """ + Rename variables in the remapped ismip6 input data + to the ones that MALI uses. + + Parameters + ---------- + remapped_file_temp : str + temporary ismip6 data remapped on mali mesh + + output_file : str + remapped ismip6 data renamed on mali mesh + """ + + # open dataset in 20 years chunk + ds = xr.open_dataset(remapped_file_temp, chunks=dict(time=20)) + + # build dictionary and rename the ismip6 dimension and variables + ismip6_to_mali_dims = dict( + ncol="nCells") + ds = ds.rename(ismip6_to_mali_dims) + + ismip6_to_mali_vars = dict( + deltaT_basin="ismip6shelfMelt_deltaT", + gamma0="ismip6shelfMelt_gamma0") + + # drop unnecessary variables on the regridded data on MALI mesh + ds = ds.rename(ismip6_to_mali_vars) + ds = ds.drop_vars(["lat_vertices", "lon_vertices", + "lat", "lon", "area"]) + + # write to a new netCDF file + write_netcdf(ds, output_file) + ds.close() + + # input files: input uniform melt rate coefficient (gamma0) + # and temperature correction per basin + _file_basin = { + "2100": "AIS/Ocean_Forcing/imbie2/imbie2_basin_numbers_8km.nc", + "2300": "AIS/Ocean_forcing/imbie2/imbie2_basin_numbers_8km.nc" + } + + _files_coeff = { + "2100": ["AIS/Ocean_Forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_local_5th_pct_PIGL_gamma_calibration.nc", + "AIS/Ocean_Forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_local_5th_percentile.nc", + "AIS/Ocean_Forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_local_95th_pct_PIGL_gamma_calibration.nc", + "AIS/Ocean_Forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_local_95th_percentile.nc", + "AIS/Ocean_Forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_local_median_PIGL_gamma_calibration.nc", + "AIS/Ocean_Forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_local_median.nc", + "AIS/Ocean_Forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_non_local_5th_pct_PIGL_gamma_calibration.nc", + "AIS/Ocean_Forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_non_local_5th_percentile.nc", + "AIS/Ocean_Forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_non_local_95th_pct_PIGL_gamma_calibration.nc", + "AIS/Ocean_Forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_non_local_95th_percentile.nc", + "AIS/Ocean_Forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_non_local_median_PIGL_gamma_calibration.nc", + "AIS/Ocean_Forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_non_local_median.nc"], + + "2300": ["AIS/Ocean_forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_local_5th_pct_PIGL_gamma_calibration.nc", + "AIS/Ocean_forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_local_5th_percentile.nc", + "AIS/Ocean_forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_local_95th_pct_PIGL_gamma_calibration.nc", + "AIS/Ocean_forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_local_95th_percentile.nc", + "AIS/Ocean_forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_local_median_PIGL_gamma_calibration.nc", + "AIS/Ocean_forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_local_median.nc", + "AIS/Ocean_forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_non_local_5th_pct_PIGL_gamma_calibration.nc", + "AIS/Ocean_forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_non_local_5th_percentile.nc", + "AIS/Ocean_forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_non_local_95th_pct_PIGL_gamma_calibration.nc", + "AIS/Ocean_forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_non_local_95th_percentile.nc", + "AIS/Ocean_forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_non_local_median_PIGL_gamma_calibration.nc", + "AIS/Ocean_forcing/parameterizations/coeff_gamma0_DeltaT_quadratic_non_local_median.nc"] + } \ No newline at end of file diff --git a/compass/landice/tests/ismip6_forcing/ocean_thermal/__init__.py b/compass/landice/tests/ismip6_forcing/ocean_thermal/__init__.py new file mode 100644 index 0000000000..5b627c5c9a --- /dev/null +++ b/compass/landice/tests/ismip6_forcing/ocean_thermal/__init__.py @@ -0,0 +1,50 @@ +from compass.testcase import TestCase +from compass.landice.tests.ismip6_forcing.ocean_thermal.\ + process_thermal_forcing import ProcessThermalForcing +from compass.landice.tests.ismip6_forcing.configure import configure as \ + configure_testgroup + + +class OceanThermal(TestCase): + """ + A test case for processing the ISMIP6 ocean thermal forcing data. + The test case builds a mapping file for interpolation between the + ISMIP6 8km polarstereo grid and MALI mesh, regrids the forcing data + and renames the ISMIP6 variables to corresponding MALI variables. + + Attributes + ---------- + process_obs : bool + Whether we are processing observations rather than CMIP model data + """ + + def __init__(self, test_group, process_obs): + """ + Create the test case + + Parameters + ---------- + test_group : compass.landice.tests.ismip6_forcing.Ismip6Forcing + The test group that this test case belongs to + + process_obs: bool + Whether we are processing observations rather than CMIP model data + """ + if process_obs: + name = "ocean_thermal_obs" + else: + name = "ocean_thermal" + self.process_obs = process_obs + super().__init__(test_group=test_group, name=name) + + step = ProcessThermalForcing(test_case=self, + process_obs=self.process_obs) + self.add_step(step) + + def configure(self): + """ + Configures test case + """ + check_model_options = not self.process_obs + configure_testgroup(config=self.config, + check_model_options=check_model_options) diff --git a/compass/landice/tests/ismip6_forcing/ocean_thermal/process_thermal_forcing.py b/compass/landice/tests/ismip6_forcing/ocean_thermal/process_thermal_forcing.py new file mode 100644 index 0000000000..1fb29fa9c9 --- /dev/null +++ b/compass/landice/tests/ismip6_forcing/ocean_thermal/process_thermal_forcing.py @@ -0,0 +1,347 @@ +import os +import numpy as np +import shutil +import subprocess +import xarray as xr +from compass.landice.tests.ismip6_forcing.create_mapfile \ + import build_mapping_file +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call +from compass.step import Step + + +class ProcessThermalForcing(Step): + """ + A step for creating a mesh and initial condition for dome test cases + + Attributes + ---------- + process_obs : bool + Whether we are processing observations rather than CMIP model data + """ + def __init__(self, test_case, process_obs): + """ + Create the step + + Parameters + ---------- + test_case : compass.landice.tests.ismip6_forcing.ocean_thermal. + OceanThermal + The test case this step belongs to + + process_obs : bool + Whether we are processing observations rather than CMIP model data + """ + super().__init__(test_case=test_case, name="process_thermal_forcing", + ntasks=4, min_tasks=1) + self.process_obs = process_obs + + def setup(self): + """ + Set up this step of the test case + """ + config = self.config + section = config["ismip6_ais"] + base_path_ismip6 = section.get("base_path_ismip6") + base_path_mali = section.get("base_path_mali") + mali_mesh_name = section.get("mali_mesh_name") + mali_mesh_file = section.get("mali_mesh_file") + period_endyear = section.get("period_endyear") + model = section.get("model") + scenario = section.get("scenario") + + process_obs_data = self.process_obs + + self.add_input_file(filename=mali_mesh_file, + target=os.path.join(base_path_mali, + mali_mesh_file)) + + if process_obs_data: + input_file = self._file_obs + output_file = f"{mali_mesh_name}_obs_TF_1995-2017_8km_x_60m.nc" + else: + input_file = self._files[period_endyear][model][scenario] + output_file = f"{mali_mesh_name}_TF_{model}_{scenario}_" \ + f"{period_endyear}.nc" + + self.add_input_file(filename=os.path.basename(input_file[0]), + target=os.path.join(base_path_ismip6, + input_file[0])) + self.add_output_file(filename=output_file) + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + + section = config["ismip6_ais"] + mali_mesh_name = section.get("mali_mesh_name") + mali_mesh_file = section.get("mali_mesh_file") + period_endyear = section.get("period_endyear") + model = section.get("model") + scenario = section.get("scenario") + output_base_path = section.get("output_base_path") + + section = config["ismip6_ais_ocean_thermal"] + method_remap = section.get("method_remap") + process_obs_data = self.process_obs + + if process_obs_data: + input_file_list = self._file_obs + output_file = f"{mali_mesh_name}_obs_TF_1995-2017_8km_x_60m.nc" + output_path = f"{output_base_path}/ocean_thermal_forcing/"\ + f"obs" + else: + input_file_list = self._files[period_endyear][model][scenario] + output_file = f"{mali_mesh_name}_TF_" \ + f"{model}_{scenario}_{period_endyear}.nc" + output_path = f"{output_base_path}/ocean_thermal_forcing/" \ + f"{model}_{scenario}/1995-{period_endyear}" + + input_file = os.path.basename(input_file_list[0]) + + # interpolate and rename the ismip6 thermal forcing data + remapped_file_temp = "remapped.nc" # temporary file name + + # call the function that reads in, remap and rename the file. + logger.info("Calling the remapping function...") + self.remap_ismip6_thermal_forcing_to_mali_vars(input_file, + remapped_file_temp, + mali_mesh_name, + mali_mesh_file, + method_remap) + + # call the function that renames the ismip6 variables to MALI variables + logger.info(f"Renaming the ismip6 variables to mali variable names...") + self.rename_ismip6_thermal_forcing_to_mali_vars(remapped_file_temp, + output_file) + + logger.info(f"Remapping and renamping process done successfully. " + f"Removing the temporary file 'remapped.nc'...") + + # remove the temporary combined file + os.remove(remapped_file_temp) + + # place the output file in appropriate directory + if not os.path.exists(output_path): + logger.info(f"Creating a new directory for the output data...") + os.makedirs(output_path) + + src = os.path.join(os.getcwd(), output_file) + dst = os.path.join(output_path, output_file) + shutil.copy(src, dst) + + logger.info(f"!---Done processing the file---!") + + def remap_ismip6_thermal_forcing_to_mali_vars(self, + input_file, + output_file, + mali_mesh_name, + mali_mesh_file, + method_remap): + """ + Remap the input ismip6 thermal forcing data onto mali mesh + + Parameters + ---------- + input_file: str + ismip6 thermal forcing data on its native polarstereo 8km grid + + output_file : str + ismip6 data remapped on mali mesh + + mali_mesh_name : str + name of the mali mesh used to name mapping files + + mali_mesh_file : str, optional + The MALI mesh file if mapping file does not exist + + method_remap : str, optional + Remapping method used in building a mapping file + """ + + # check if a mapfile + mapping_file = f"map_ismip6_8km_to_{mali_mesh_name}_{method_remap}.nc" + + if not os.path.exists(mapping_file): + # build a mapping file if it doesn't already exist + build_mapping_file(self.config, self.ntasks, self.logger, + input_file, mapping_file, mali_mesh_file, + method_remap) + else: + self.logger.info(f"Mapping file exists. " + f"Remapping the input data...") + + # remap the input data + args = ["ncremap", + "-i", input_file, + "-o", output_file, + "-m", mapping_file] + + check_call(args, logger=self.logger) + + def rename_ismip6_thermal_forcing_to_mali_vars(self, + remapped_file_temp, + output_file): + """ + Rename variables in the remapped ismip6 input data + to the ones that MALI uses. + + Parameters + ---------- + remapped_file_temp : str + temporary ismip6 data remapped on mali mesh + + output_file : str + remapped ismip6 data renamed on mali mesh + """ + + config = self.config + process_obs_data = self.process_obs + + # open dataset in 20 years chunk + ds = xr.open_dataset(remapped_file_temp, chunks=dict(time=20), + engine="netcdf4") + + ds["ismip6shelfMelt_zOcean"] = ds.z + ds = ds.drop_vars("z") # dropping 'z' while it's still called 'z' + + # build dictionary for ismip6 variables that MALI takes in + if process_obs_data: + ismip6_to_mali_dims = dict( + z="nISMIP6OceanLayers", + ncol="nCells") + ds["xtime"] = ("Time", ["2015-01-01_00:00:00".ljust(64)]) + ds["xtime"] = ds.xtime.astype("S") + ds["thermal_forcing"] = ds["thermal_forcing"].expand_dims( + dim="Time", axis=0) + ds = ds.rename(ismip6_to_mali_dims) + else: + ismip6_to_mali_dims = dict( + z="nISMIP6OceanLayers", + time="Time", + ncol="nCells") + ds = ds.rename(ismip6_to_mali_dims) + # add xtime variable + xtime = [] + for t_index in range(ds.sizes["Time"]): + date = ds.Time[t_index] + # forcing files do not all match even years, + # so round up the years + # pandas round function does not work for years, + # so do it manually + yr = date.dt.year.values + mo = date.dt.month.values + dy = date.dt.day.values + dec_yr = np.around(yr + (30 * (mo - 1) + dy) / 365.0) + date = f"{dec_yr.astype(int)}-01-01_00:00:00".ljust(64) + xtime.append(date) + + ds["xtime"] = ("Time", xtime) + ds["xtime"] = ds.xtime.astype('S') + ds = ds.drop_vars(["Time"]) + + ismip6_to_mali_vars = dict( + thermal_forcing="ismip6shelfMelt_3dThermalForcing") + ds = ds.rename(ismip6_to_mali_vars) + + # drop unnecessary variables + ds = ds.drop_vars(["z_bnds", "lat_vertices", "area", + "lon_vertices", "lat", "lon"]) + + # transpose dimension + ds["ismip6shelfMelt_3dThermalForcing"] = \ + ds["ismip6shelfMelt_3dThermalForcing"].transpose( + "Time", "nCells", "nISMIP6OceanLayers") + + # write to a new netCDF file + write_netcdf(ds, output_file) + ds.close() + + # create a nested dictionary for the ISMIP6 original forcing files including relative path + _file_obs = ["AIS/Ocean_Forcing/climatology_from_obs_1995-2017/obs_thermal_forcing_1995-2017_8km_x_60m.nc"] + _files = { + "2100": { + "CCSM4": { + "RCP85": ["AIS/Ocean_Forcing/ccsm4_rcp8.5/1995-2100/CCSM4_thermal_forcing_8km_x_60m.nc"] + }, + "CESM2": { + "SSP585v1": [ + "AIS/Ocean_Forcing/cesm2_ssp585/1995-2100/CESM2_ssp585_thermal_forcing_8km_x_60m.nc"], + "SSP585v2": [ + "AIS/Ocean_Forcing/cesm2_ssp585/1995-2100/CESM2_ssp585_thermal_forcing_8km_x_60m.nc"], + }, + "CNRM_CM6": { + "SSP126": [ + "AIS/Ocean_Forcing/cnrm-cm6-1_ssp126/1995-2100/CNRM-CM6-1_ssp126_thermal_forcing_8km_x_60m.nc"], + "SSP585": [ + "AIS/Ocean_Forcing/cnrm-cm6-1_ssp585/1995-2100/CNRM-CM6-1_ssp585_thermal_forcing_8km_x_60m.nc"] + }, + "CNRM_ESM2": { + "SSP585": [ + "AIS/Ocean_Forcing/cnrm-esm2-1_ssp585/1995-2100/CNRM-ESM2-1_ssp585_thermal_forcing_8km_x_60m.nc"] + }, + "CSIRO-Mk3-6-0": { + "RCP85": [ + "AIS/Ocean_Forcing/csiro-mk3-6-0_rcp8.5/1995-2100/CSIRO-Mk3-6-0_RCP85_thermal_forcing_8km_x_60m.nc"] + }, + "HadGEM2-ES": { + "RCP85": [ + "AIS/Ocean_Forcing/hadgem2-es_rcp8.5/1995-2100/HadGEM2-ES_RCP85_thermal_forcing_8km_x_60m.nc"] + }, + "IPSL-CM5A-MR": { + "RCP26": [ + "AIS/Ocean_Forcing/ipsl-cm5a-mr_rcp2.6/1995-2100/IPSL-CM5A-MR_RCP26_thermal_forcing_8km_x_60m.nc"], + "RCP85": [ + "AIS/Ocean_Forcing/ipsl-cm5a-mr_rcp8.5/1995-2100/IPSL-CM5A-MR_RCP85_thermal_forcing_8km_x_60m.nc"] + }, + "MIROC-ESM-CHEM": { + "RCP85": [ + "AIS/Ocean_Forcing/miroc-esm-chem_rcp8.5/1995-2100/MIROC-ESM-CHEM_thermal_forcing_8km_x_60m.nc"] + }, + "NorESM1-M": { + "RCP26": [ + "AIS/Ocean_Forcing/noresm1-m_rcp2.6/1995-2100/NorESM1-M_RCP26_thermal_forcing_8km_x_60m.nc"], + "RCP85": [ + "AIS/Ocean_Forcing/noresm1-m_rcp8.5/1995-2100/NorESM1-M_thermal_forcing_8km_x_60m.nc"] + }, + "UKESM1-0-LL": { + "SSP585": [ + "AIS/Ocean_Forcing/ukesm1-0-ll_ssp585/1995-2100/UKESM1-0-LL_ssp585_thermal_forcing_8km_x_60m.nc"] + } + }, + "2300": { + "CCSM4": { + "RCP85": [ + "AIS/Ocean_forcing/ccsm4_RCP85/1995-2300/CCSM4_RCP85_thermal_forcing_8km_x_60m.nc"] + }, + "CESM2-WACCM": { + "SSP585": [ + "AIS/Ocean_forcing/cesm2-waccm_ssp585/1995-2299/CESM2-WACCM_SSP585_thermal_forcing_8km_x_60m.nc"], + "SSP585-repeat": [ + "AIS/Ocean_forcing/cesm2-waccm_ssp585-repeat/1995-2300/CESM2-WACCM_ssp585_thermal_forcing_8km_x_60m.nc"] + }, + "HadGEM2-ES": { + "RCP85": [ + "AIS/Ocean_forcing/hadgem2-es_RCP85/1995-2299/HadGEM2-ES_RCP85_thermal_forcing_8km_x_60m.nc"], + "RCP85-repeat": [ + "AIS/Ocean_forcing/hadgem2-es_RCP85-repeat/1995-2300/HadGEM2-ES_rcp85_thermal_forcing_8km_x_60m.nc"] + }, + "NorESM1-M": { + "RCP26-repeat": [ + "AIS/Ocean_forcing/noresm1-m_RCP26-repeat/1995-2300/NorESM1-M_RCP26_thermal_forcing_8km_x_60m.nc"], + "RCP85-repeat": [ + "AIS/Ocean_forcing/noresm1-m_RCP85-repeat/1995-2300/NorESM1-M_thermal_forcing_8km_x_60m.nc"] + }, + "UKESM1-0-LL": { + "SSP126": [ + "AIS/Ocean_forcing/ukesm1-0-ll_ssp126/1995-2300/UKESM1-0-LL_thermal_forcing_8km_x_60m_v2.nc"], + "SSP585": [ + "AIS/Ocean_forcing/ukesm1-0-ll_ssp585/1995-2300/UKESM1-0-LL_SSP585_thermal_forcing_8km_x_60m.nc"], + "SSP585-repeat": [ + "AIS/Ocean_forcing/ukesm1-0-ll_ssp585-repeat/1995-2300/UKESM1-0-LL_ssp585_thermal_forcing_8km_x_60m.nc"] + } + } + } diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index 8499c2f410..e5b6ee9977 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -231,6 +231,54 @@ hydro_radial visualize.Visualize.run visualize.visualize_hydro_radial +ismip6_forcing +~~~~~~~~~~~~~~ + +.. currentmodule:: compass.landice.tests.ismip6_forcing + +.. autosummary:: + :toctree: generated/ + + Ismip6Forcing + configure.configure + create_mapfile.build_mapping_file + + atmosphere.Atmosphere + atmosphere.Atmosphere.configure + atmosphere.create_mapfile_smb + atmosphere.create_mapfile_smb.build_mapping_file + atmosphere.create_mapfile_smb.create_atm_scrip + atmosphere.process_smb.ProcessSMB + atmosphere.process_smb.ProcessSMB.setup + atmosphere.process_smb.ProcessSMB.run + atmosphere.process_smb.ProcessSMB.remap_ismip6_smb_to_mali + atmosphere.process_smb.ProcessSMB.rename_ismip6_smb_to_mali_vars + atmosphere.process_smb.ProcessSMB.correct_smb_anomaly_for_climatology + atmosphere.process_smb_racmo.ProcessSmbRacmo + atmosphere.process_smb_racmo.ProcessSmbRacmo.setup + atmosphere.process_smb_racmo.ProcessSmbRacmo.run + atmosphere.process_smb_racmo.ProcessSmbRacmo.remap_source_smb_to_mali + atmosphere.process_smb_racmo.ProcessSmbRacmo.rename_source_smb_to_mali_vars + atmosphere.process_smb_racmo.ProcessSmbRacmo.correct_smb_anomaly_for_base_smb + + ocean_basal.OceanBasal + ocean_basal.OceanBasal.configure + ocean_basal.process_basal_melt.ProcessBasalMelt + ocean_basal.process_basal_melt.ProcessBasalMelt.setup + ocean_basal.process_basal_melt.ProcessBasalMelt.run + ocean_basal.process_basal_melt.ProcessBasalMelt.combine_ismip6_inputfiles + ocean_basal.process_basal_melt.ProcessBasalMelt.remap_ismip6_basal_melt_to_mali_vars + ocean_basal.process_basal_melt.ProcessBasalMelt.rename_ismip6_basal_melt_to_mali_vars + + 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 + ocean_thermal.process_thermal_forcing.ProcessThermalForcing.remap_ismip6_thermal_forcing_to_mali_vars + ocean_thermal.process_thermal_forcing.ProcessThermalForcing.rename_ismip6_thermal_forcing_to_mali_vars + + kangerlussuaq ~~~~~~~~~~~~~ diff --git a/docs/developers_guide/landice/test_groups/index.rst b/docs/developers_guide/landice/test_groups/index.rst index 94760fc976..02e83acaa6 100644 --- a/docs/developers_guide/landice/test_groups/index.rst +++ b/docs/developers_guide/landice/test_groups/index.rst @@ -15,6 +15,7 @@ Test groups greenland humboldt hydro_radial + ismip6_forcing kangerlussuaq koge_bugt_s mismipplus diff --git a/docs/developers_guide/landice/test_groups/ismip6_forcing.rst b/docs/developers_guide/landice/test_groups/ismip6_forcing.rst new file mode 100644 index 0000000000..281fef9d8d --- /dev/null +++ b/docs/developers_guide/landice/test_groups/ismip6_forcing.rst @@ -0,0 +1,64 @@ +.. _dev_landice_ismip6_forcing: + +ismip6_forcing +============== + +The ``ismip6_forcing`` test group (:py:class:`compass.landice.tests.ismip6_ +forcing.Ismip6Forcing`) processes (i.e., remapping and renaming) the +atmospheric and oceanic forcing data of the Ice Sheet Model +Intercomparison for CMIP6 (ISMIP6) protocol from its native polarstereo grid to +the unstructure MALI mesh. The test group includes four test cases, +``atmosphere``, ``ocean_basal``, ``ocean_thermal_obs`` and ``ocean_thermal``. +The ``atmosphere`` test case has two steps: ``process_smb`` and ``process_smb_racmo``; +the ``ocean_basal`` test case has one step, ``process_basal_melt``; +the ``ocean_thermal_obs`` and ``ocean_thermal`` share one step, ``process_thermal_forcing``. +Each step has the local methods (functions) of remapping and +renaming the original ismip6 data to the format that MALI can incorporate in +its forward simulations. In remapping the data, all test cases import the +method ``build_mapping_file`` to create or use scrip files +of the source (ISMIP6) and destination (MALI) mesh depending on the existence +of a mapping file. Below, we describe the shared framework for this +test group and the 3 test cases. + +.. _dev_landice_ismip6_forcing_framework: + +framework +--------- + +The shared config options for the ``ismip6_forcing`` test group are described +in :ref:`landice_ismip6_forcing` in the User's Guide. + + +.. _dev_landice_ismip6_forcing_atmosphere: + +atmosphere +---------- + +The :py:class:`compass.landice.tests.ismip6_forcing.atmosphere.Atmosphere` +performs processing of the surface mass balance (SMB) forcing. +Processing data includes regridding the original ISMIP6 SMB data from its +native polarstereo grid to MALI's unstructured grid, renaming variables and +correcting the SMB anomaly field for the MALI base SMB. + +.. _dev_landice_ismip6_forcing_ocean_basal: + +ocean_basal +------------ + +The :py:class:`compass.landice.tests.ismip6_forcing.ocean_basal.OceanBasal` +performs processing of the coefficients for the basal melt parameterization +utilized by the ISMIP6 protocol. Processing data includes combining the +IMBIE2 basin number file and parameterization coefficients and remapping onto +the MALI mesh. + +.. _dev_landice_ismip6_forcing_ocean_thermal: + +ocean_thermal +------------- + +The :py:class:`compass.landice.tests.ismip6_forcing.ocean_thermal.OceanThermal` +performs the processing of ocean thermal forcing, both observational climatology +(in ``ocean_thermal_obs``) and the CMIP model data (in ``ocean_thermal``). +Processing data includes regridding the original ISMIP6 thermal forcing data +from its native polarstereo grid to MALI's unstructured grid and renaming variables. + diff --git a/docs/users_guide/landice/test_groups/index.rst b/docs/users_guide/landice/test_groups/index.rst index 044b85404f..4123aae7c7 100644 --- a/docs/users_guide/landice/test_groups/index.rst +++ b/docs/users_guide/landice/test_groups/index.rst @@ -20,6 +20,7 @@ physics but that are not run routinely. greenland humboldt hydro_radial + ismip6_forcing kangerlussuaq koge_bugt_s mismipplus diff --git a/docs/users_guide/landice/test_groups/ismip6_forcing.rst b/docs/users_guide/landice/test_groups/ismip6_forcing.rst new file mode 100644 index 0000000000..a16436d33b --- /dev/null +++ b/docs/users_guide/landice/test_groups/ismip6_forcing.rst @@ -0,0 +1,314 @@ +.. _landice_ismip6_forcing: + +ismip6_forcing +============== + +The ``landice/ismip6_forcing`` test group processes (i.e., remaps and renames) +the atmospheric and ocean forcing data of the Ice Sheet Model Intercomparison for CMIP6 +(ISMIP6) protocol. The processed data is used to force MALI in its simulations +under a relevant ISMIP6 (either the 2100 or 2300) experimental protocol. +The test group includes three test cases, ``atmosphere``, ``ocean_basal``, +``ocean_thermal_obs`` and ``ocean_thermal``; the ``atmosphere`` test case +has two steps: ``process_smb`` and ``process_smb_racmo``. The ``ocean_basal`` +and the ``ocean_thermal`` test case each has one step, ``process_basal_melt``, +and ``process_thermal_forcing``, respectively. (For more details on the steps of +each test case, see :ref:`landice_ismip6_forcing_atmosphere`, +:ref:`landice_ismip6_forcing_ocean_basal`, +:ref:`landice_ismip6_forcing_ocean_thermal_obs` and +:ref:`landice_ismip6_forcing_ocean_thermal`.) +Approximated time for processing a single forcing file +on Cori (single core) is 2 and 7 minutes for the atmosphere and ocean basal +testcases, and less than a minute for ocean thermal obs and ocean thermal +testcases, respectively. + +Before providing the details of necessary source data of the ISMIP6 protocols, +we provide a summary of instructions of an overall process of this test group: +To start off, users need to provide a MALI mesh onto which the source data +will be remapped and renamed. More information about obtaining these meshes will +be provided as soon as they are publicly available. Then, for a given MALI mesh, + +1. run :ref:`landice_ismip6_forcing_ocean_basal` once, independent of the +model, scenario and end year. + +2. run :ref:`landice_ismip6_forcing_ocean_thermal_obs` once, independent +of the model, scenario and end year, to process the thermal forcing from the observational climatology (used for +control runs). + +3. run :ref:`landice_ismip6_forcing_ocean_thermal` for each model, scenario +and end year. + +4. run :ref:`landice_ismip6_forcing_atmosphere` with +``process_racmo_smb = True`` once with any model. + +5. run :ref:`landice_ismip6_forcing_atmosphere` for each model, +scenario and end year. Users can keep ``process_racmo_smb = False`` as long as +the RACMO SMB has been processed once in Step 4, but it is harmless to leave +``process_racmo_smb = True`` as it does nothing if data is already +available, and the processing is very quick (less than a minute). + +There are six different of input data sets other than the MALI mesh that +users need in order to use this package: (#1) atmospheric surface mass balance (SMB) +anomaly forcing and (#2) ocean thermal forcing for projection runs, modern +climatology files for (#3) atmospheric and (#4) ocean forcing, (#5) ISMIP6 basin +number file, and (#6) the ocean basal-melt parameter files. + +Except for the file #3, The ISMIP6 source data (files #1, 2, 4-6) can be obtained by contacting the ISMIP6 steering +committee as described `here. `_ +Once the users get access to the data on `Globus `_, +and within the base path directory that the user specifies in the config file +(i.e., the config option ``base_path_ismip6``; +see :ref:`landice_ismip6_forcing_config`), they +need to manually download the files following the directory structure +and names provided in the GHub endpoints (named ``GHub-ISMIP6-Forcing`` +for the 2100 CE projection protocol and ``ISMIP6-Projections-Forcing-2300`` +for the 2300 CE projection protocol). That is, the directory paths and the file +names must exactly match (even the letter case) those that are provided by the +GHub endpoints. For example, if a user wants to process the atmosphere (SMB) +forcing representing the ``SSP585`` scenario from the ``UKESM1-0-LL`` model +provided by the ISMIP6-2100 protocol (i.e., from the ``GHub-ISMIP6-Forcing`` +Globus endpoint), they must create the directory path +``AIS/Atmosphere_Forcing/UKESM1-0-LL/Regridded_8km/`` in their local system +where they will download the file ``UKESM1-0-LL_anomaly_ssp585_1995-2100_8km.nc``. +Note that the users only need to download SMB anomaly files in 8km resolution +(the climatology file is not needed) that cover the period from 1995 to +``period_endyear`` (either 2100 or 2300, defined in the config file; +see :ref:`landice_ismip6_forcing_config`). + +Equivalently, for the +ocean forcing in this example, users should create the directory path +``AIS/Ocean_Forcing/ukesm1-0-ll_ssp585/1995-2100/`` and download the file +``UKESM1-0-LL_ssp585_thermal_forcing_8km_x_60m.nc`` from the same endpoint. +Users do not need to download the thermal +forcing files for the years previous to 1995 as only the files downloaded from +``1995-{period_endyear}`` will be processed. Users also do not need to download +the temperature and salinity files, as these will not be used by MALI. +Also note to be aware that unlike those in the ``GHub-ISMIP6-Forcing`` endpoint, +the directory names in the ``ISMIP6-Projections-Forcing-2300`` endpoint have a +lower case "f" for the ``AIS/Atmospheric_forcing/`` and ``AIS/Ocean_forcing/``. + +In addition to atmospheric and ocean thermal forcing files that +correspond to specific climate model (e.g., UKESM1-0-LL, CCSM4) and scenarios +(e.g., SSP585, RCP85, RCP26-repeat), modern +climatology files are needed. For the ``atmosphere`` testcase, +``RACMO2.3p2_ANT27_smb_yearly_1979_2018.nc`` will be automatically downloaded +from the MALI public database when the testcase is being set up and saved +to the directory that users define in the config option `landice_database_root`. +The RACMO file is used to correct the ISMIP6 the surface mass balance (SMB) +data with the modern climatology. For the ``ocean_thermal`` case, users need to +download the modern ocean thermal forcing climatology file named +``obs_thermal_forcing_1995-2017_8km_x_60m.nc`` in the directory +``AIS/Ocean_F{f}orcing/climatology_from_obs_1995-2017/`` +(the salinity and temperature files do not have to be downloaded). + + +For the ``ocean_basal`` testcase, users need to additionally download +the basin number file ``imbie2_basin_numbers_8km.nc`` in the directory +``AIS/Ocean_Forcing/imbie2/`` (or ``AIS/Ocean_forcing/imbie2/``, if from the +``ISMIP6-Projections-Forcing-2300`` endpoint); all of the files that +start their name with ``coeff_gamma0_DeltaT_quadratic_local`` in the directory +''AIS/Ocean_F{f}orcing/parameterizations/'', which contain parameter values needed +for calculating the basal melt underneath the ice shelves in MALI simulations. + +Note that both the RACMO SMB data and ocean basal-melt parameters not +associated with any climate models and scenarios and thus can be processed only +once and can be applied to MALI with any set of processed climate forcing data. + + +In the next section (ref:`landice_ismip6_forcing_config`), we provide +instructions and examples of how users can configure necessary options including +paths to necessary source files and the output path of the processed data +within which the subdirectories called ``atmosphere_forcing/``, ``basal_melt/`` +and ``ocean_thermal_forcing/`` (and further subdirectories that match the source +file directory structure) are created if the directories do not already exist) +and where processed files will be saved. + +.. _landice_ismip6_forcing_config: + +config options +-------------- + +All four test cases share some set of default config options under the section +``[ismip6_ais]`` and have separate config options for each test case: +``[ismip6_ais_atmosphere]``, ``[ismip6_ais_ocean_thermal]``, and +``[ismip6_ais_ocean_basal]``. In the general config section +``[ismip6_ais]``, users need to supply base paths to input files and MALI mesh +file, and MALI mesh name, as well as the model name, climate forcing scenario +and the projection end year of the ISMIP6 forcing data, which can be chosen +from the available options as given in the config file (see the example file +below.) In the ``ismip6_ais_atmosphere`` section, users need to indicate +``True`` or ``False`` on whether to process the RACMO modern climatology +(``True`` is required to run the ``process_smb_racmo`` step, which needs to be +run before the ``process_smb`` step). + +For most the ``[ismip6_ais_atmosphere]`` and ``[ismip6_ais_ocean_thermal]`` +config sections users may choose the interpolation scheme among +``bilinear``, ``neareststod`` and ``conserve`` methods. The exception is that +the ``ocean basal`` test case should always use the ``neareststod`` method +because the source files have a single valued data per basin. + +Below is the default config options: + +.. code-block:: cfg + + # config options for ismip6 antarctic ice sheet data set + [paths] + # The root to a location where data files for MALI will be cached + landice_database_root = /Users/hollyhan/Desktop/RESEARCH/MALI/database/ + + [ismip6_ais] + + # Base path to the input ismip6 ocean and smb forcing files. User has to supply. + base_path_ismip6 = /Users/hollyhan/Desktop/ISMIP6_2300_Protocol/ISMIP6-Projections-Forcing-2300/ + + # Base path to the the MALI mesh. User has to supply. + base_path_mali = /Users/hollyhan/Desktop/RESEARCH/MALI/mesh_files/ + + # Forcing end year of the ISMIP6 data. User has to supply. + # Available end years are 2100 and 2300. + period_endyear = 2300 + + # Base path to which output forcing files are saved. + output_base_path = /Users/hollyhan/Desktop/ISMIP6_2300_Protocol/Process_Forcing_Testcase/ + + # Name of climate model name used to generate ISMIP6 forcing data. User has to supply. + # Available model names for the 2100 projection are the following: CCSM4, CESM2, CNRM_CM6, CNRM_ESM2, CSIRO-Mk3-6-0, HadGEM2-ES, IPSL-CM5A-MR, MIROC-ESM-CHEM, NorESM1-M, UKESM1-0-LL + # Available model names for the 2300 projection are the following: CCSM4, CESM2-WACCM, CSIRO-Mk3-6-0, HadGEM2-ES, NorESM1-M, UKESM1-0-LL + model = NorESM1-M + + # Scenarios used by climate model. User has to supply. + # Available scenarios for the 2100 projection are the following: RCP26, RCP26-repeat, RCP85, SSP126, SSP585 (SSP585v1 and SSP585v2 for the CESM2 model) + # Available scenarios for the 2300 projection are the following: RCP26, RCP26-repeat, RCP85, RCP85-repeat, SSP126, SSP585, SSP585-repeat + scenario = RCP26-repeat + + # name of the mali mesh. User has to supply. Note: It is used to name mapping files + # (e,g. 'map_ismip6_8km_to_{mali_mesh_name}_{method_remap}.nc'). + mali_mesh_name = Antarctica_8to30km + + # MALI mesh file to be used to build mapping file (e.g.Antarctic_8to80km_20220407.nc). User has to supply. + mali_mesh_file = AIS_8to30km_r01_20220607.nc + + # config options for ismip6 antarctic ice sheet SMB forcing data test cases + [ismip6_ais_atmosphere] + + # Remapping method used in building a mapping file. Options include: bilinear, neareststod, conserve + method_remap = bilinear + + # Set True to process RACMO modern climatology + process_smb_racmo = True + + # config options for ismip6 ocean thermal forcing data test cases + [ismip6_ais_ocean_thermal] + + # Remapping method used in building a mapping file. Options include: bilinear, neareststod, conserve + method_remap = bilinear + + # Set to True if the want to process observational thermal forcing data. Set to False if want to process model thermal forcing data. + process_obs_data = True + +Below is the example config options that users might create in running +the test group. This example is for processing the NorESM1-M RCP2.6 repeat +forcing to the year 2300 onto the 8-80km Antarctic Ice Sheet MALI mesh. +The example is configured to perform the `atmosphere\process_smb_racmo` step to +process the RACMO modern SMB climatology but not the modern thermal forcing. + +.. code-block:: cfg + + # config options for ismip6 antarctic ice sheet data set + [paths] + # The root to a location where data files for MALI will be cached + landice_database_root = NotAvailable + + [ismip6_ais] + + # Base path to the input ismip6 ocean and smb forcing files. User has to supply. + base_path_ismip6 = NotAvailable + + # Base path to the the MALI mesh. User has to supply. + base_path_mali = NotAvailable + + # Forcing end year of the ISMIP6 data. User has to supply. + # Available end years are 2100 and 2300. + period_endyear = NotAvailable + + # Base path to which output forcing files are saved. + output_base_path = NotAvailable + + # Name of climate model name used to generate ISMIP6 forcing data. User has to supply. + # Available model names for the 2100 projection are the following: CCSM4, CESM2, CNRM_CM6, CNRM_ESM2, CSIRO-Mk3-6-0, HadGEM2-ES, IPSL-CM5A-MR, MIROC-ESM-CHEM, NorESM1-M, UKESM1-0-LL + # Available model names for the 2300 projection are the following: CCSM4, CESM2-WACCM, CSIRO-Mk3-6-0, HadGEM2-ES, NorESM1-M, UKESM1-0-LL + model = NotAvailable + + # Scenarios used by climate model. User has to supply. + # Available scenarios for the 2100 projection are the following: RCP26, RCP26-repeat, RCP85, SSP126, SSP585 (SSP585v1 and SSP585v2 for the CESM2 model) + # Available scenarios for the 2300 projection are the following: RCP26, RCP26-repeat, RCP85, RCP85-repeat, SSP126, SSP585, SSP585-repeat + scenario = NotAvailable + + # name of the mali mesh. User has to supply. Note: It is used to name mapping files + # (e,g. 'map_ismip6_8km_to_{mali_mesh_name}_{method_remap}.nc'). + mali_mesh_name = NotAvailable + + # MALI mesh file to be used to build mapping file (e.g.Antarctic_8to80km_20220407.nc). User has to supply. + mali_mesh_file = NotAvailable + + # config options for ismip6 antarctic ice sheet SMB forcing data test cases + [ismip6_ais_atmosphere] + + # Remapping method used in building a mapping file. Options include: bilinear, neareststod, conserve + method_remap = bilinear + + # Set True to process RACMO modern climatology + process_smb_racmo = True + + # config options for ismip6 ocean thermal forcing data test cases + [ismip6_ais_ocean_thermal] + + # Remapping method used in building a mapping file. Options include: bilinear, neareststod, conserve + method_remap = bilinear + + # Set to True if the want to process observational thermal forcing data. Set to False if want to process model thermal forcing data. + process_obs_data = True + +.. _landice_ismip6_forcing_atmosphere: + +atmosphere +---------- + +The ``landice/ismip6_forcing/atmosphere`` test case +performs processing of the surface mass balance (SMB) forcing data provided by +the ISMIP6 and RACMO. Processing data includes regridding the SMB forcing data +SMB data from the native grid (polarstereo grid for the ISMIP6 files and +rotated pole grid for the RACMO file) to MALI's unstructured grid, renaming +variables, and correcting the ISMIP6 SMB anomaly field for the base SMB +(modern climatology) provided by RACMO. + +.. _landice_ismip6_forcing_ocean_basal: + +ocean_basal +------------ + +The ``landice/tests/ismip6_forcing/ocean_basal`` test case +performs processing of the coefficients for the basal melt parameterization +utilized by the ISMIP6 protocol. Processing data includes combining the +IMBIE2 basin numbers file and parameterization coefficients and remapping onto +the MALI mesh. + +.. _landice_ismip6_forcing_ocean_thermal_obs: + +ocean_thermal_obs +----------------- + +The ``landice/ismip6_forcing/ocean_thermal_obs`` test case +performs the processing of the observational climatology of +ocean thermal forcing. Processing data includes regridding the original ISMIP6 +thermal forcing data from its native polarstereo grid to MALI's unstructured +grid and renaming variables. + +.. _landice_ismip6_forcing_ocean_thermal: + +ocean_thermal +------------- + +The ``landice/ismip6_forcing/ocean_thermal`` test case +performs the processing of ocean thermal forcing. Processing data includes +regridding the original ISMIP6 thermal forcing data from its native +polarstereo grid to MALI's unstructured grid and renaming variables.