diff --git a/compass/landice/__init__.py b/compass/landice/__init__.py index 97b33a010c..338f5237a2 100644 --- a/compass/landice/__init__.py +++ b/compass/landice/__init__.py @@ -16,6 +16,7 @@ from compass.landice.tests.koge_bugt_s import KogeBugtS from compass.landice.tests.mesh_modifications import MeshModifications from compass.landice.tests.mismipplus import MISMIPplus +from compass.landice.tests.slm import Slm from compass.landice.tests.thwaites import Thwaites from compass.mpas_core import MpasCore @@ -49,4 +50,5 @@ def __init__(self): self.add_test_group(KogeBugtS(mpas_core=self)) self.add_test_group(MeshModifications(mpas_core=self)) self.add_test_group(MISMIPplus(mpas_core=self)) + self.add_test_group(Slm(mpas_core=self)) self.add_test_group(Thwaites(mpas_core=self)) diff --git a/compass/landice/tests/slm/__init__.py b/compass/landice/tests/slm/__init__.py new file mode 100644 index 0000000000..0ed3787e7d --- /dev/null +++ b/compass/landice/tests/slm/__init__.py @@ -0,0 +1,17 @@ +from compass.landice.tests.slm.circ_icesheet import CircIcesheetTest +from compass.testgroup import TestGroup + + +class Slm(TestGroup): + """ + A test group for Sea-Level Model test cases + """ + 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='slm') + + self.add_test_case( + CircIcesheetTest(test_group=self)) diff --git a/compass/landice/tests/slm/circ_icesheet/__init__.py b/compass/landice/tests/slm/circ_icesheet/__init__.py new file mode 100644 index 0000000000..75b373cd7f --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/__init__.py @@ -0,0 +1,65 @@ +from compass.landice.tests.slm.circ_icesheet.run_model import RunModel +from compass.landice.tests.slm.circ_icesheet.setup_mesh import SetupMesh +from compass.landice.tests.slm.circ_icesheet.visualize import Visualize +from compass.testcase import TestCase + + +class CircIcesheetTest(TestCase): + """ + This test generates an idealized, circular ice sheet that has a + prescribed thickness evolution for testing coupling between MALI + and the Sea-Level Model. + """ + + def __init__(self, test_group): + """ + Create the test case + + Parameters + ---------- + test_group : compass.landice.tests.slm.Slm + The test group that this test case belongs to + The resolution or type of mesh of the test case + """ + name = 'circular_icesheet_test' + subdir = name + super().__init__(test_group=test_group, name=name, + subdir=subdir) + + def configure(self): + """ + Set up the desired mesh-resolution tests + + Read the list of resolutions from the config + """ + config = self.config + section = config['circ_icesheet'] + mali_res = section.get('mali_res').split(',') + + section = config['slm'] + slm_nglv = section.get('slm_nglv').split(',') + print(f'list of MALI-mesh resolution is {mali_res} km.') + print(f'list of SLM Gauss-Legendre latitudinal points is {slm_nglv}.') + + for res in mali_res: + for nglv in slm_nglv: + self.add_step(SetupMesh(test_case=self, + name=f'mali{res}km_slm{nglv}/' + 'setup_mesh', res=res, nglv=nglv)) + if (int(res) <= 16 and int(res) > 2): + ntasks = 256 + elif (int(res) <= 2): + ntasks = 512 + else: + ntasks = 128 + min_tasks = ntasks + self.add_step(RunModel(test_case=self, res=res, nglv=nglv, + ntasks=ntasks, min_tasks=min_tasks, + openmp_threads=1, + name=f'mali{res}km_slm{nglv}' + '/run_model')) + step = Visualize(test_case=self) + self.add_step(step, run_by_default=True) + + # no run() method is needed because we're doing the default: running all + # steps diff --git a/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg new file mode 100644 index 0000000000..8525d54101 --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/circ_icesheet.cfg @@ -0,0 +1,80 @@ +# Config options for slm circ_icesheet test case +[circ_icesheet] + +# The size of the domain in meters in the x and y directions for the planar mesh +lx = 6000000.0 +ly = 6000000.0 + +# list of MALI mesh resolutions in kilometers delimited by ',' without space +# i.e. The distance between adjacent cell centers. +mali_res = 20 + +# Ice shape type ('cylinder', and for dome, 'dome-halfar' or 'dome-cism') +# r0 and h0: initial radius and height of the cylinder ice in meters +ice_type = cylinder +r0 = 2000000.0 +h0 = 3000.0 + +# 'True' if manually want to set bedTopography elevation +# if 'False' topography will be set to 0 everywhere (flat bed at sea level) +set_topo_elev = True +# flat bed topography elevation in meters +# (positive value => land, negative => ocean) +topo_elev = -1500.0 + +# Whether to center the circular ice in the center of the cell that is closest to the +# center of the domain +put_origin_on_a_cell = False + +# Config options for surface mass balance forcing +[smb_forcing] + +# Start, end and interval years for SMB forcing +start_year = 2015 +end_year = 2116 +dt_year = 1 + +# Direction in which SMB is applied ('horizontal' or 'vertical') +# and amount of ice to melt +direction = horizontal +# Change in radius in meters; used when 'direction == horizontal' +drdt = -2000.0 +# Change in height in meters;sed when 'direction == vertical' +dhdt = -20.0 + +# config options for the sea-level model +[slm] +# True if MALI-SLM are coupled +coupling = True + +# List of the number of Gauss-Legendre points in latitude +# list delimited by ',' without space +slm_nglv = 512 + +# Max spherical harmonics degree and order +# i.e. SLM resolution of the SLM +slm_res = 512 + +# mapping method between the MALI and SLM grids +mapping_method_mali_to_slm = conserve +mapping_method_slm_to_mali = bilinear + +# ratio of MALI-SLM coupling interval and MALI output interval in integer years +time_stride = 1 + +# config options related to visualization +[circ_icesheet_viz] + +# Area (m^2) of the global ocean for calculating ice sheet contribution to +# sea level. Only used when MALI-SLM coupling is False +Aocn_const = 4.5007E+14 + +# Area (m^2) of the global ocean excluding the marine-based ice region +# only used when MALI-SLM coupling is False +AocnBeta_const = 4.50007E+14 + +# whether to save image files +save_images = True + +# whether to hide figures (typically when save_images = True) +hide_figs = True diff --git a/compass/landice/tests/slm/circ_icesheet/namelist.landice b/compass/landice/tests/slm/circ_icesheet/namelist.landice new file mode 100644 index 0000000000..37851f0e0d --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/namelist.landice @@ -0,0 +1,53 @@ +&velocity_solver + config_velocity_solver = 'none' +/ +&advection + config_thickness_advection = 'fo' + config_tracer_advection = 'none' +/ +&solidearth + config_uplift_method = 'sealevelmodel' + config_slm_coupling_interval = 5 + config_MALI_to_SLM_weights_file = 'mapping_file_mali_to_slm.nc' + config_SLM_to_MALI_weights_file = 'mapping_file_slm_to_mali.nc' +/ +&calving + config_calving = 'none' + config_restore_calving_front = .false. +/ +&thermal_solver + config_thermal_solver = 'none' +/ +&iceshelf_melt + config_basal_mass_bal_float = 'none' + config_front_mass_bal_grounded = 'none' +/ +&physical_parameters + config_ice_density = 910.0 + config_ocean_density = 1028.0 + config_sea_level = 0.0 + config_dynamic_thickness = 10.0 +/ +&time_integration + config_dt = '0001-00-00_00:00:00' + config_time_integration = 'forward_euler' + config_adaptive_timestep = .false. +/ +&time_management + config_do_restart = .false. + config_restart_timestamp_name = 'restart_timestamp' + config_start_time = '2015-01-01_00:00:00' + config_stop_time = '2055-01-01_00:00:00' + config_run_duration = 'none' + config_calendar_type = 'noleap' +/ +&io + config_pio_stride = 128 +/ +&AM_globalStats + 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. +/ diff --git a/compass/landice/tests/slm/circ_icesheet/run_model.py b/compass/landice/tests/slm/circ_icesheet/run_model.py new file mode 100644 index 0000000000..1f9084a658 --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/run_model.py @@ -0,0 +1,103 @@ +import os +from importlib import resources + +from jinja2 import Template + +from compass.model import run_model +from compass.step import Step + + +class RunModel(Step): + """ + A step for performing forward MALI runs as part of dome test cases. + """ + + def __init__(self, test_case, res, nglv, ntasks, name='run_model', + subdir=None, min_tasks=None, openmp_threads=1): + """ + Create a new test case + + Parameters + ---------- + test_case : compass.TestCase + The test case this step belongs to + + res : str + Resolution of the MALI domain + + nglv : str + Number of Gauss-Legendre nodes in latitude in the SLM grid + + ntasks : int + the number of tasks the step would ideally use. If fewer tasks + are available on the system, the step will run on all available + tasks as long as this is not below ``min_tasks`` + + name : str, optional + the name of the test case + + subdir : str, optional + the subdirectory for the step. The default is ``name`` + + min_tasks : int, optional + the number of tasks the step requires. If the system has fewer + than this number of tasks, the step will fail + + openmp_threads : int, optional + the number of OpenMP threads the step will use + """ + self.res = res + self.nglv = nglv + + if min_tasks is None: + min_tasks = ntasks + super().__init__(test_case=test_case, name=name, subdir=subdir, + ntasks=ntasks, min_tasks=min_tasks, + openmp_threads=openmp_threads) + + self.add_namelist_file( + 'compass.landice.tests.slm.circ_icesheet', 'namelist.landice', + out_name='namelist.landice') + + self.add_streams_file( + 'compass.landice.tests.slm.circ_icesheet', 'streams.landice', + out_name='streams.landice') + + self.add_input_file(filename='landice_grid.nc', + target='../setup_mesh/landice_grid.nc') + self.add_input_file(filename='graph.info', + target='../setup_mesh/graph.info') + self.add_input_file(filename='smb_forcing.nc', + target='../setup_mesh/smb_forcing.nc') + self.add_input_file(filename='mapping_file_mali_to_slm.nc', + target='../setup_mesh/' + 'mapping_file_mali_to_slm.nc') + self.add_input_file(filename='mapping_file_slm_to_mali.nc', + target='../setup_mesh/' + 'mapping_file_slm_to_mali.nc') + self.add_model_as_input() + + self.add_output_file(filename='output.nc') + + def setup(self): + 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') + + # change the sealevel namelist + template = Template(resources.read_text + ('compass.landice.tests.slm', + 'namelist.sealevel.template')) + text = template.render(nglv=self.nglv) + + # write out the namelist.sealevel file + file_slm_nl = os.path.join(self.work_dir, 'namelist.sealevel') + with open(file_slm_nl, 'w') as handle: + handle.write(text) + + def run(self): + """ + Run this step of the test case + """ + run_model(step=self) diff --git a/compass/landice/tests/slm/circ_icesheet/setup_mesh.py b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py new file mode 100644 index 0000000000..f6f1f16702 --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/setup_mesh.py @@ -0,0 +1,409 @@ +import os +import shutil + +import mpas_tools.io +import netCDF4 +import numpy as np +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call +from mpas_tools.mesh.conversion import cull +from mpas_tools.planar_hex import make_planar_hex_mesh +from mpas_tools.scrip.from_mpas import scrip_from_mpas +from mpas_tools.translate import center, translate +from netCDF4 import Dataset as NetCDFFile + +from compass.model import make_graph_file +from compass.step import Step + + +class SetupMesh(Step): + """ + A step for creating a mesh with initial conditions plus forcing file for + circular icesheet test + + Attributes + ---------- + res : str + Resolution of MALI mesh + + nglv : str + Number of Gauss-Legendre nodes in latitude in the SLM + """ + + def __init__(self, test_case, name, res, nglv): + """ + Create the step + + Parameters + ---------- + test_case : compass.TestCase + The test case this step belongs to + + res : str + Resolution of MALI mesh + + nglv : str + Number of Gauss-Legendre nodes in latitude in the SLM + """ + super().__init__(test_case=test_case, name=name) + + self.add_output_file(filename='graph.info') + self.add_output_file(filename='landice_grid.nc') + + self.res = res + self.nglv = nglv + # no setup() method is needed + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + section = config['circ_icesheet'] + + lx = section.getfloat('lx') + ly = section.getfloat('ly') + dc = float(self.res) * 1000.0 + + # calculate the number of grid cells in x and y direction + # for the uniform, hexagonal planar mesh + nx = max(2 * int(0.5 * lx / dc + 0.5), 4) + # factor of 2/sqrt(3) because of hexagonal mesh + ny = max(2 * int(0.5 * ly * (2. / np.sqrt(3)) / dc + 0.5), 4) + + mpas_tools.io.default_format = 'NETCDF4' + mpas_tools.io.default_engine = 'netcdf4' + + # call the mesh creation function + dsMesh = make_planar_hex_mesh(nx=nx, ny=ny, dc=dc, + nonperiodic_x=True, + nonperiodic_y=True) + dsMesh = cull(dsMesh, logger=logger) + # adding the time dimension is needed for netcdf4 formatting to work + dsMesh['xtime'] = ('Time', ['2015-01-01_00:00:00'.ljust(64)]) + # translating the mesh center to x=0 & y=0 + center(dsMesh) + # shift the center to a quarter or radius + shift = 200000.0 + print(f'shifting the center by {shift} meters') + translate(dsMesh, shift, shift) + + fname_culled = 'culled_mesh_before_cdf5.nc' + write_netcdf(dsMesh, fname_culled) + args = ['ncks', '-O', '-5', fname_culled, 'mpas_grid.nc'] + check_call(args, logger=logger) + + levels = 3 + args = ['create_landice_grid_from_generic_MPAS_grid.py', + '-i', 'mpas_grid.nc', + '-o', 'landice_grid.nc', + '-l', str(levels)] + + check_call(args, logger) + + make_graph_file(mesh_filename='landice_grid.nc', + graph_filename='graph.info') + + _setup_circsheet_initial_conditions(config, logger, + filename='landice_grid.nc') + + _create_smb_forcing_file(config, logger, + mali_mesh_file='landice_grid.nc', + filename='smb_forcing.nc') + + _build_mapping_files(config, logger, self.res, self.nglv, + mali_mesh_file='landice_grid.nc') + + os.remove(fname_culled) + + +def _setup_circsheet_initial_conditions(config, logger, filename): + """ + Create initial condition for circular ice sheet + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for this test case, a combination of the defaults + for the machine, core and configuration + + logger : logging.Logger + A logger for output from the step + + filename : str + file to setup circular icesheet + """ + section = config['circ_icesheet'] + ice_type = section.get('ice_type') + set_topo_elev = section.getboolean('set_topo_elev') + topo_elev = section.get('topo_elev') + r0 = section.getfloat('r0') + h0 = section.getfloat('h0') + put_origin_on_a_cell = section.getboolean('put_origin_on_a_cell') + + # Open the file, get needed dimensions + gridfile = NetCDFFile(filename, 'r+') + nVertLevels = len(gridfile.dimensions['nVertLevels']) + # Get variables + xCell = gridfile.variables['xCell'] + yCell = gridfile.variables['yCell'] + xEdge = gridfile.variables['xEdge'] + yEdge = gridfile.variables['yEdge'] + xVertex = gridfile.variables['xVertex'] + yVertex = gridfile.variables['yVertex'] + thickness = gridfile.variables['thickness'] + bedTopography = gridfile.variables['bedTopography'] + layerThicknessFractions = gridfile.variables['layerThicknessFractions'] + + # Find center of domain + x0 = xCell[:].min() + 0.5 * (xCell[:].max() - xCell[:].min()) + y0 = yCell[:].min() + 0.5 * (yCell[:].max() - yCell[:].min()) + # Calculate distance of each cell center from dome center + r = ((xCell[:] - x0) ** 2 + (yCell[:] - y0) ** 2) ** 0.5 + + if put_origin_on_a_cell: + # Center the ice in the center of the cell that is closest to the + # center of the domain. + centerCellIndex = np.abs(r[:]).argmin() + xShift = -1.0 * xCell[centerCellIndex] + yShift = -1.0 * yCell[centerCellIndex] + xCell[:] = xCell[:] + xShift + yCell[:] = yCell[:] + yShift + xEdge[:] = xEdge[:] + xShift + yEdge[:] = yEdge[:] + yShift + xVertex[:] = xVertex[:] + xShift + yVertex[:] = yVertex[:] + yShift + # Now update origin location and distance array + x0 = 0.0 + y0 = 0.0 + r = ((xCell[:] - x0) ** 2 + (yCell[:] - y0) ** 2) ** 0.5 + + # Assign variable values for the circular ice sheet + # Set default value for non-circular cells + thickness[:] = 0.0 + # Calculate the dome thickness for cells within the desired radius + # (thickness will be NaN otherwise) + thickness_field = thickness[0, :] + + if ice_type == 'cylinder': + logger.info('cylinder ice type is chosen') + thickness_field[r < r0] = h0 + elif ice_type == 'dome-cism': + thickness_field[r < r0] = h0 * (1.0 - (r[r < r0] / r0) ** 2) ** 0.5 + elif ice_type == 'dome-halfar': + thickness_field[r < r0] = h0 * ( + 1.0 - (r[r < r0] / r0) ** (4.0 / 3.0)) ** (3.0 / 7.0) + else: + raise ValueError('Unexpected ice_type: {}'.format(ice_type)) + thickness[0, :] = thickness_field + + # flat bed at sea level + bedTopography[:] = 0.0 + if set_topo_elev: + # this line will make a small shelf: + bedTopography[:] = topo_elev + + # Setup layerThicknessFractions + layerThicknessFractions[:] = 1.0 / nVertLevels + + gridfile.close() + + logger.info('Successfully added circular initial conditions to: {}'.format( + filename)) + + +def _create_smb_forcing_file(config, logger, mali_mesh_file, filename): + """ + Create surface mass balance forcing file for circular ice sheet + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for this test case, a combination of the defaults + for the machine, core and configuration + + logger : logging.Logger + A logger for output from the step + + mali_mesh_file : str + mesh file created in the current step + + filename : str + file to setup circular icesheet + """ + section = config['circ_icesheet'] + r0 = section.getfloat('r0') + h0 = section.getfloat('h0') + + section = config['smb_forcing'] + direction = section.get('direction') + start_year = int(section.get('start_year')) + end_year = int(section.get('end_year')) + dt_year = int(section.get('dt_year')) + dhdt = section.getfloat('dhdt') + drdt = section.getfloat('drdt') + + # other constants used in the fuction + # pi = 3.1415926535898 # pi value used in the SLM + rhoi = 910.0 # ice density + + # Open the file, get needed dimensions + shutil.copy(mali_mesh_file, filename) + smbfile = NetCDFFile(filename, 'r+') + + # Get variables + xCell = smbfile.variables['xCell'] + yCell = smbfile.variables['yCell'] + nCells = smbfile.dimensions['nCells'].size + + # time interval for the forcing + t_array = np.arange(start_year, end_year, dt_year) + nt = len(t_array) + + # initialize SMB value everywhere in the mesh + sfcMassBal = np.zeros((nt, nCells)) + + # Find center of domain + x0 = xCell[:].min() + 0.5 * (xCell[:].max() - xCell[:].min()) + y0 = yCell[:].min() + 0.5 * (yCell[:].max() - yCell[:].min()) + # calculate the radius of circ ice sheet + r = ((xCell[:] - x0) ** 2 + (yCell[:] - y0) ** 2) ** 0.5 + + if direction == 'vertical': + logger.info('creating a file that prescribes vertical smb') + smb = dhdt * rhoi / (365.0 * 24 * 60 * 60) + logger.info(f'prescribed smb is {smb} kg/m2/s') + indx = np.where(r <= r0)[0] + # for loop through time + for t in range(len(t_array)): + Ht = h0 + (dhdt * (t + 1)) # calculate a new height each time + logger.info(f'At time {start_year + dt_year}, \ + new height will be {Ht} km') + sfcMassBal[t, indx] = smb # assign sfcMassBal + smbfile.variables['sfcMassBal'][:, :] = sfcMassBal[:, :] + + elif direction == 'horizontal': + logger.info('creating a file that prescribes horizontal smb') + smb = -1 * ((10000.0) * rhoi / (365.0 * 24.0 * 60.0 * 60.0)) + print(f'prescribed smb is {smb} kg/m2/s') + # for loop through time + for t in range(len(t_array)): + # calculate a new radius at each time + Rt = r0 + (drdt * (t + 1)) + if (Rt < 0): + Rt = 0.0 + sfcMassBal[t, :] = 0.0 + print(f'At time {start_year+dt_year*t}, \ + new radius will be {Rt} km') + else: + indx = np.where(r >= Rt)[0] + sfcMassBal[t, indx] = smb + print(f'At time {start_year+dt_year*t}, \ + new radius will be {Rt/1000} km') + smbfile.variables['sfcMassBal'][:, :] = sfcMassBal[:, :] + + # add xtime variable + smbfile.createDimension('StrLen', 64) + xtime = smbfile.createVariable('xtime', 'S1', ('Time', 'StrLen')) + + # initialize SMB value everywhere in the mesh + xtime_str = [] + for t_index in range(len(t_array)): + yr = start_year + (t_index * dt_year) + xtime_str = f'{int(yr)}-01-01_00:00:00'.ljust(64) + xtime_char = netCDF4.stringtochar(np.array([xtime_str], 'S64'), + encoding='utf-8') + xtime[t_index, :] = xtime_char + + smbfile.close() + + +def _build_mapping_files(config, logger, res, nglv, mali_mesh_file): + """ + 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 + + logger : logging.Logger + A logger for output from the step + + name : str + Name of the step + + res : str + Resolution of MALI mesh + + nglv : str + Number of Gauss-Legendre nodes in latitude in the SLM + + mali_mesh_file : str + The MALI mesh file + """ + + section = config['slm'] + slm_res = section.get('slm_res') + slm_nglv = int(nglv) + method_mali_to_slm = section.get('mapping_method_mali_to_slm') + method_slm_to_mali = section.get('mapping_method_slm_to_mali') + + mali_scripfile = f'mali{int(res)}km_scripfile.nc' + slm_scripfile = f'slm{slm_res}_nglv{slm_nglv}scripfile.nc' + + # create slm scripfile + logger.info(f'creating scripfile for the SH-degree {slm_res} ' + f'resolution and {slm_nglv} nglv points for the SLM grid') + args = ['ncremap', + '-g', slm_scripfile, + '-G', + f'latlon={slm_nglv},{2*int(slm_nglv)}#lat_typ=gss#lat_drc=n2s'] + + check_call(args, logger=logger) + + # adjust the lat-lon values + args = ['set_lat_lon_fields_in_planar_grid.py', + '--file', mali_mesh_file, + '--proj', 'ais-bedmap2-sphere'] + + check_call(args, logger=logger) + + # create scrip files for source and destination grids + logger.info('creating scrip file for the mali mesh') + scrip_from_mpas(mali_mesh_file, mali_scripfile) + + # create a mapping file using ESMF weight gen + logger.info('Creating a mapping file...') + + parallel_executable = config.get("parallel", "parallel_executable") + # split the parallel executable into constituents in case it includes flags + args = parallel_executable.split(' ') + args.extend(['ESMF_RegridWeightGen', + '-s', mali_scripfile, + '-d', slm_scripfile, + '-w', 'mapping_file_mali_to_slm.nc', + '-m', method_mali_to_slm, + '-i', '-64bit_offset', '--netcdf4', + '--src_regional']) + + check_call(args, logger) + + args = parallel_executable.split(' ') + args.extend(['ESMF_RegridWeightGen', + '-s', slm_scripfile, + '-d', mali_scripfile, + '-w', 'mapping_file_slm_to_mali.nc', + '-m', method_slm_to_mali, + '-i', '-64bit_offset', '--netcdf4', + '--dst_regional']) + + check_call(args, logger) + + # remove the temporary scripfiles once the mapping file is generated + logger.info('Removing the temporary mesh and scripfiles...') + os.remove(slm_scripfile) + os.remove(mali_scripfile) diff --git a/compass/landice/tests/slm/circ_icesheet/streams.landice b/compass/landice/tests/slm/circ_icesheet/streams.landice new file mode 100644 index 0000000000..0a5dc1eaf0 --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/streams.landice @@ -0,0 +1,67 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/compass/landice/tests/slm/circ_icesheet/visualize.py b/compass/landice/tests/slm/circ_icesheet/visualize.py new file mode 100644 index 0000000000..b029f4ecde --- /dev/null +++ b/compass/landice/tests/slm/circ_icesheet/visualize.py @@ -0,0 +1,650 @@ +import os + +import matplotlib.pyplot as plt +import numpy as np +import xarray as xr + +from compass.step import Step + + +class Visualize(Step): + """ + A step for visualizing the output from a dome test case + + Attributes + ---------- + """ + def __init__(self, test_case, name='visualize', subdir=None, + input_dir='run_model'): + """ + Create the step + + Parameters + ---------- + test_case : compass.TestCase + The test case this step belongs to + + name : str, optional + the name of the test case + + subdir : str, optional + the subdirectory for the step. The default is ``name`` + + input_dir : str, optional + The input directory within the test case with a file ``output.nc`` + to visualize + """ + super().__init__(test_case=test_case, name=name, subdir=subdir) + + def setup(self): + """ + + """ + config = self.config + section = config['circ_icesheet'] + mali_res = section.get('mali_res').split(',') + + section = config['slm'] + slm_nglv = section.get('slm_nglv').split(',') + + for res in mali_res: + for nglv in slm_nglv: + self.add_input_file(filename=f'output_mali{res}km_' + f'slm{nglv}.nc', + target=f'../mali{res}km_slm{nglv}/' + 'run_model/output/output.nc') + + def run(self): + """ + Run this step of the test case + """ + logger = self.logger + config = self.config + section = config['circ_icesheet'] + mali_res = section.get('mali_res').split(',') + + section = config['slm'] + coupling = section.getboolean('coupling') + slm_nglv = section.get('slm_nglv').split(',') + + section = config['circ_icesheet_viz'] + save_images = section.getboolean('save_images') + + # visualize run model results + num_files = 0 + for res in mali_res: + for nglv in slm_nglv: + run_path = f'../mali{res}km_slm{nglv}/run_model/' + logger.info(f'analyzing outputs in path: {run_path}') + visualize_slm_circsheet(config, logger, res, nglv) + num_files = num_files + 1 + + # calculate and plot rmse in SLC + if coupling and num_files > 1: + # ncells_list = list() + # rmse_list = list() + pct_err_slc_list = list() + pct_err_slc_z0_list = list() + pct_err_slc_z0_matrix = np.zeros((len(mali_res), len(slm_nglv))) + n = 0 + for res in mali_res: + for nglv in slm_nglv: + run_path = f'../mali{res}km_slm{nglv}/run_model/' + run_data = output_analysis(config, logger, + res, nglv, run_path) + # ncells_res = run_data.ncells + # rmse_res = run_data.rmse_slc + pct_err_slc_res = run_data.pct_err_slc + pct_err_slc_z0_res = run_data.pct_err_slm_z0 + # ncells_list.append(ncells_res) + # rmse_list.append(rmse_res) + pct_err_slc_list.append(pct_err_slc_res) + pct_err_slc_z0_list.append(pct_err_slc_z0_res) + + # ncells = np.array(ncells_list) + # rmse = np.array(rmse_list) + pct_err_slc = np.array(pct_err_slc_list) + pct_err_slc_z0 = np.array(pct_err_slc_z0_list) + + # assign the percentage error array in matrix + pct_err_slc_z0_matrix[n, :] = pct_err_slc_z0 + n += 1 + + # plot mean percentage error for fixed MALI res + plot_pct_err_slc(res, slm_nglv, pct_err_slc, + 'Mean Percent Error in SLC') + plot_pct_err_slc(res, slm_nglv, pct_err_slc_z0, + 'Mean Percent Error in SLC-z0') + # reset the lists for percentage error + del pct_err_slc_list[:], pct_err_slc_z0_list[:] + + # plot the heat map of mean percentage errors + plot_heatmap(pct_err_slc_z0_matrix, mali_res, slm_nglv, + save_images) + + +def plot_heatmap(data, mali_res, slm_nglv, save_images): + """ + Plot a heat map of mean percentage error in SLC. + + Parameters + ---------- + data : float + Matrix of arrays containing mean percentage error values + + mali_res : str + List of MALI-resolution strings + + slm_nglv : str + List of SLM nglv values + + save_images : str + Boolean on whether to save the plot or not + """ + + fig1 = plt.figure(1, facecolot='w') + ax1 = fig1.add_subplot(1, 1, 1) + im = ax1.imshow(data, cmap='Blues') + ax1.figure.colorbar(im, ax=ax1) + ax1.set_title('Mean Percentage Error (%)') + ax1.set_xticks(np.arange(len(slm_nglv)), labels=slm_nglv) + ax1.set_yticks(np.arange(len(mali_res)), labels=mali_res) + ax1.set_xlabel('# of GL nodes in Latitude') + ax1.set_ylabel('MALI domain resolution (km)') + if save_images: + plt.savefig('Heatmap_MPE.png') + fig1.clf() + + +def plot_pct_err_slc(mali_res, slm_nglv, data, figure_title): + """ + Plot mean percentage errors in sea-level change where + MALI-calculation is taken as a true value. + + Parameters + ---------- + mali_res : str + Resolution of the MALI domain + + slm_nglv : int + Array of the number of Gauss-Legendre nodes in latitude + in the SLM grid + + data : float + Data array containing the percentage error values + + figure_title : str + Title of the figure + """ + + fig1 = plt.figure(1, facecolor='w') + ax1 = fig1.add_subplot(1, 1, 1) + ax1.plot(np.array(slm_nglv), data, marker='o') + ax1.set_xlabel('# of Gauss-Legendre nodes in latitude (SLM)', fontsize=14) + ax1.legend(loc='best', ncol=1, prop={'size': 10}) + ax1.set_title(f'{figure_title} for {mali_res}km MALI res') + ax1.set_ylabel('percent error (%)', fontsize=14) + plt.savefig(f'{figure_title} for {mali_res}km MALI res.png', + bbox_inches='tight') + fig1.clf() + + +def visualize_slm_circsheet(config, logger, res, nglv): + """ + Plot the output from a circular ice sheet test case + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for this test case, + a combination of the defaults for the machine, + core and configuration + + logger : logging.Logger + A logger for output from the step + + res : str + resolution of MALI mesh + + nglv : str + Number of Gauss-Legendre points in latitude in the SLM grid + """ + section = config['circ_icesheet_viz'] + save_images = section.getboolean('save_images') + hide_figs = section.getboolean('hide_figs') + + section = config['slm'] + coupling = section.getboolean('coupling') + + # get an instance of output analysis class + run_path = f'../mali{res}km_slm{nglv}/run_model/' + run_data = output_analysis(config, logger, res, nglv, run_path) + yrs = run_data.yrs + + # figure 1 + fig1 = plt.figure(1, facecolor='w') + ax1 = fig1.add_subplot(1, 1, 1) + ax1.plot(yrs, run_data.grnd_vol_unscaled, + label='MALI unscaled', linestyle=':', color='k') + ax1.plot(yrs, run_data.grnd_vol, label='MALI scaled', + linestyle='-', color='k') + if coupling: + ax1.plot(yrs, run_data.ice_vol_slm, label='SLM', + linestyle='-.', color='r') + ax1.set_xlabel('Year') + ax1.legend(loc='best', ncol=1, prop={'size': 10}) + ax1.set_title('Grounded Ice Mass') + ax1.set_ylabel('Mass (Gt)') + ax1.grid(True) + if save_images: + plt.savefig(f'ice_mass_mali{res}km_slm{nglv}.png') + fig1.clf() + + # figure 2 + fig2 = plt.figure(2, facecolor='w') + ax2 = fig2.add_subplot(1, 1, 1) + ax2.plot(yrs, run_data.dgrnd_vol_unscaled, + label='MALI unscaled', linestyle=':', color='k') + ax2.plot(yrs, run_data.dgrnd_vol, label='MALI scaled', + linestyle='-', color='k') + if coupling: + ax2.plot(yrs, run_data.dice_vol_slm, label='SLM', + linestyle='-.', color='r') + ax2.set_xlabel('Year') + ax2.legend(loc='best', ncol=1, prop={'size': 10}) + ax2.set_title('Change in Grounded Ice Mass') + ax2.set_ylabel('Mass (Gt)') + ax2.grid(True) + if save_images: + plt.savefig(f'ice_mass_change_mali{res}km_slm{nglv}.png') + fig2.clf() + + # figure 3 + fig3 = plt.figure(3, facecolor='w') + ax3 = fig3.add_subplot(1, 1, 1) + ax3.plot(yrs, run_data.SLCaf, label='MALI SLC VAF', + linestyle=':', color='k') + ax3.plot(yrs, run_data.SLCcorr_Aocn, + label='MALI SLCcorr (grnded ice flooded)', + linestyle='-.', color='k') + ax3.plot(yrs, run_data.SLCcorr_AocnBeta, + label='MALI SLCcorr (grnded ice NOT flooded)', + linestyle='-', color='b') + ax3.plot(yrs, run_data.SLCcorr_z0_AocnBeta, + label='MALI SLCcorr+z0 (grnded ice NOT flooded)', + linestyle='-', color='y') + if coupling: + ax3.plot(yrs, run_data.SLC_slm_Aocn, + label='SLM (grnded ice flooded)', linestyle='-.', color='r') + ax3.plot(yrs, run_data.SLC_slm_AocnBeta, + label='SLM (grnded ice NOT flooded)', linestyle='-', + color='r') + ax3.set_xlabel('Year') + ax3.legend(loc='best', ncol=1, prop={'size': 10}) + ax3.set_title('Ice-sheet contribution to SLC') + ax3.set_ylabel('Sea-level change (m)') + ax3.grid(True) + if save_images: + plt.savefig(f'ice_contribution_to_SLC_mali{res}km_slm{nglv}.png') + fig3.clf() + + # figure 4 & 5 + if coupling: + fig4 = plt.figure(4, facecolor='w') + ax4 = fig4.add_subplot(1, 1, 1) + ax4.plot(yrs, run_data.dgrnd_vol - run_data.dice_vol_slm, + label='MALI minus SLM', linestyle='-', color='k') + ax4.set_xlabel('Year') + ax4.legend(loc='best', ncol=1, prop={'size': 10}) + ax4.set_title('Difference in ice mass change') + ax4.set_ylabel('Mass (Gt)') + ax4.grid(True) + if save_images: + plt.savefig(f'diff_ice_mass_change_mali{res}km_slm{nglv}.png') + fig4.clf() + + fig5 = plt.figure(5, facecolor='w') + ax5 = fig5.add_subplot(1, 1, 1) + ax5.plot(yrs, run_data.diff_slc, label='MALI minus SLM', + linestyle='-', color='k') + ax5.plot(yrs, run_data.diff_slc_z0, label='MALI minus SLM (z0)', + linestyle='-', color='r') + ax5.set_xlabel('Year') + ax5.legend(loc='best', ncol=1, prop={'size': 10}) + ax5.set_title('Difference in sea-level change') + ax5.set_ylabel('sea-level change (m)') + ax5.grid(True) + if save_images: + plt.savefig(f'diff_sea_level_change_mali{res}km_slm{nglv}.png') + fig5.clf() + + if hide_figs: + logger.info("Plot display disabled with hide_plot config option.") + else: + plt.show() + + +class output_analysis: + """ + Analyze outputs + + Attributes + ---------- + mali_slc : float + ice-sheet contribution to sea-level change calculated + based on MALI outputs + + slm_slc : float + ice-sheet contribution to sea-level change calculated + based on SLM outputs + + rmse : float + root mean square error between mali_slc and slm_slc + """ + def __init__(self, config, logger, res, nglv, run_path): + """ + Calculate sea-level change from run outputs + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for this test case, + a combination of the defaults for the machine, + core and configuration + + logger : logging.Logger + A logger for output from the step + + res : str + Resolution of MALI mesh + + nglv : str + Number of Gauss-Legendre points in latitude in the SLM + + run_path : str + Path to runs where the output file exists + """ + self.config = config + self.logger = logger + self.res = res + self.nglv = nglv + self.run_path = run_path + + section = config['slm'] + coupling = section.getboolean('coupling') + + if coupling: + section = config['circ_icesheet_viz'] + Aocn_const = section.getfloat('Aocn_const') + AocnBeta_const = section.getfloat('AocnBeta_const') + + # mali output file name + fname_mali = f'output_mali{res}km_slm{nglv}.nc' + + # read in the MALI outputs + DS = xr.open_mfdataset(fname_mali, combine='nested', + concat_dim='Time', + decode_timedelta=False) + + # default constants + rhoi = 910.0 # ice density in kg/m^3 + rhoo = 1000.0 # ocean density used by the SLM (MALI uses 1028.0) + rhow = 1000.0 # fresh water density + Aocn_const = 4.5007E+14 # area of global ocean in m2 + AocnBeta_const = 4.5007E+14 # ocean area including marine-based area + + self.ncells = DS.dims['nCells'] + cellMask = DS['cellMask'] + bed = DS['bedTopography'] + # bed0 = bed.isel(Time=0).load() + thickness = DS['thickness'] + areaCell = DS['areaCell'][0, :].load() + latCell = DS['latCell'][0, :].load() + lonCell = DS['lonCell'][0, :].load() - np.pi # correct range: [-pi pi] + self.lat_deg = latCell * 180 / np.pi + + # calculate the area distortion (scale) factor for the + # polar stereographic projection + self.k = np.zeros(len(latCell),) + k0 = (1 - np.sin(-71 * np.pi / 180)) / 2 + # center of the latitude + lat0 = -90 * np.pi / 180 + # standard parallel in radians where distortion should be zero (k=1) + lon0 = 0 # center of the longitude + # expression for 'k' from p. 157 + # eqn. 21-4 in https://pubs.usgs.gov/pp/1395/report.pdf. + # p. 142 provides a table showing 'k' for stereographic projection + self.k[:] = 2 * k0 / (1 + (np.sin(lat0) * np.sin(latCell[:])) + + (np.cos(lat0) * np.cos(latCell[:]) * + np.cos(lonCell[:] - lon0))) + + # default MALI time steps from the MALI outputs + self.yrs_mali = DS['daysSinceStart'].load() / 365.0 + 2015.0 + if coupling: + # path to the SLM output data + fpath_slm = os.path.join(run_path, 'OUTPUT_SLM/') + config = self.config + section = config['slm'] + time_stride = int(section.get('time_stride')) + # reset the time indices to match the # of SLM timesteps + nt = np.arange(0, DS.dims['Time'], time_stride, dtype='i') + nt = np.arange(0, 50, time_stride, dtype='i') + self.yrs = self.yrs_mali[nt] + nyrs = len(self.yrs) + z0 = np.zeros((nyrs, )) + Aocn = np.zeros((nyrs, )) + AocnBeta = np.zeros((nyrs, )) + # get ice mass on the SLM interpolated from the MALI mesh + fname = os.path.join(fpath_slm, 'ice_volume') + self.ice_vol_slm = slm_outputs(fname, nyrs).data * rhoi / 1.0e12 + self.dice_vol_slm = slm_outputs(fname, nyrs).change_total * \ + rhoi / 1.0e12 # in Gt + # get slc correction and ocean area from the SLM + fname = os.path.join(fpath_slm, 'gmslc_ocnArea') + if os.path.exists(fname): + logger.info(f'reading in file {fname}') + z0 = slm_outputs(fname, nyrs).change_total + self.SLC_slm_Aocn = slm_outputs(fname, nyrs).data + + fname = os.path.join(fpath_slm, 'gmslc_ocnBetaArea') + if os.path.exists(fname): + logger.info(f'reading in file {fname}') + self.SLC_slm_AocnBeta = slm_outputs(fname, nyrs).data + + fname = os.path.join(fpath_slm, 'ocean_area') + if os.path.exists(fname): + logger.info(f'reading in file {fname}') + Aocn = slm_outputs(fname, nyrs).data + # logger.info(f'area of the ocean is: {Aocn}') + fname = os.path.join(fpath_slm, 'oceanBeta_area') + if os.path.exists(fname): + logger.info(f'reading in file {fname}') + AocnBeta = slm_outputs(fname, nyrs).data + # logger.info(f'area of the ocean Beta is: {AocnBeta}') + else: + # if not coupled, use default values for MALI outputs + # assuming MALI output interval is 1 year + nt = np.arange(0, DS.dims['Time'], 1, dtype='i') + self.yrs = self.yrs_mali + z0 = np.zeros(nyrs, ) + Aocn = np.zeros(nyrs, ) + AocnBeta = np.zeros(nyrs, ) + + z0[:] = 0.0 + logger.info("'gmsle_change' file doesn't exist. " + "Setting z0 to zeros") + Aocn[:] = Aocn_const + logger.info("'ocean_area' file doesn't exist. Using " + f"constant ocean area defined: {Aocn_const}m2") + AocnBeta[:] = AocnBeta_const + logger.info("'ocean_areaBeta' file doesn't exist. Using" + f"constant oceanBeta area defined: {AocnBeta_const}m2") + + # calculate ice-sheet contribution to sea-level change + # create empty arrays + self.grnd_vol = np.zeros((len(nt), )) + self.grnd_vol_unscaled = np.zeros((len(nt), )) + self.dgrnd_vol = np.zeros((len(nt), )) + self.dgrnd_vol_unscaled = np.zeros((len(nt), )) + self.vaf = np.zeros((len(nt), )) + self.vaf_z0 = np.zeros((len(nt), )) + self.pov = np.zeros((len(nt), )) + self.pov_z0 = np.zeros((len(nt), )) + self.den = np.zeros((len(nt), )) + self.SLEaf = np.zeros((len(nt), )) + self.SLCaf = np.zeros((len(nt), )) + self.SLCpov = np.zeros((len(nt), )) + self.SLCden = np.zeros((len(nt), )) + self.SLCcorr_Aocn = np.zeros((len(nt), )) + self.SLEaf_AocnBeta = np.zeros((len(nt), )) + self.SLEaf_AocnBeta = np.zeros((len(nt), )) + self.SLCaf_AocnBeta = np.zeros((len(nt), )) + self.SLCpov_AocnBeta = np.zeros((len(nt), )) + self.SLCden_AocnBeta = np.zeros((len(nt), )) + self.SLCcorr_AocnBeta = np.zeros((len(nt), )) + self.SLEaf_z0_AocnBeta = np.zeros((len(nt), )) + self.SLCaf_z0_AocnBeta = np.zeros((len(nt), )) + self.SLCpov_z0_AocnBeta = np.zeros((len(nt), )) + self.SLCcorr_z0_AocnBeta = np.zeros((len(nt), )) + + for t in np.arange(0, len(nt)): + # get appropriate time index for the MALI output data + if coupling: + indx = t * time_stride + else: + indx = t + + bedt = bed[indx, :].load() + cellMaskt = cellMask[indx, :].load() + thicknesst = thickness[indx, :].load() + # areaCellt = areaCell.sum() + # logger.info(f'area of the mali domain is {areaCellt}') + + # calculation of sea-level change following Goelzer et al. 2020 + self.grnd_vol_unscaled[t] = (areaCell * grounded(cellMaskt) * + thicknesst).sum() * rhoi / 1.0e12 + self.grnd_vol[t] = (areaCell * grounded(cellMaskt) * + thicknesst / (self.k**2)).sum() * rhoi / 1.0e12 + self.vaf[t] = (areaCell * grounded(cellMaskt) / (self.k**2) * + np.maximum(thicknesst + (rhoo / rhoi) * + np.minimum(np.zeros((self.ncells,)), + bedt), 0.0)).sum() # eqn. 1 + self.pov[t] = (areaCell / (self.k**2) * + np.maximum(0.0 * bedt, -1.0 * bedt)).sum() # eqn.8 + self.vaf_z0[t] = (areaCell * grounded(cellMaskt) / (self.k**2) * + np.maximum(thicknesst + (rhoo / rhoi) * + np.minimum(np.zeros((self.ncells,)), + bedt + z0[t]), 0.0)).sum() # eqn. 13 + self.pov_z0[t] = (areaCell / (self.k**2) * np.maximum(0.0 * bedt, + -1.0 * bedt + z0[t])).sum() # eqn. 14 + self.den[t] = (areaCell / (self.k**2) * grounded(cellMaskt) * + (rhoi / rhow - rhoi / rhoo)).sum() # eqn. 10 + + # SLC(m) using the ocean area where anywhere below sea level + # is considered as ocean (even where marine-based ice exists) + self.SLEaf[t] = (self.vaf[t] / Aocn[t]) * (rhoi / rhoo) # eqn. 2 + self.SLCaf[t] = -1.0 * (self.SLEaf[t] - self.SLEaf[0]) # eqn. 3 + self.SLCpov[t] = -1.0 * (self.pov[t] / Aocn[t] - self.pov[0] / + Aocn[t]) # eqn. 9 + self.SLCden[t] = -1.0 * (self.den[t] / Aocn[t] - self.den[0] / + Aocn[t]) # eqn. 11 + self.SLCcorr_Aocn[t] = self.SLCaf[t] + self.SLCpov[t] + \ + self.SLCden[t] + + # SLC(m) including the z0 term using the ocean+iced area + # i.e. water can't go where marine-based ice exists + self.SLEaf_AocnBeta[t] = (self.vaf[t] / AocnBeta[t]) * \ + (rhoi / rhoo) + self.SLCaf_AocnBeta[t] = -1.0 * (self.SLEaf_AocnBeta[t] - + self.SLEaf_AocnBeta[0]) + self.SLCpov_AocnBeta[t] = -1.0 * (self.pov[t] / AocnBeta[t] - + self.pov[0] / AocnBeta[t]) + self.SLCden_AocnBeta[t] = -1.0 * (self.den[t] / AocnBeta[t] - + self.den[0] / AocnBeta[t]) + self.SLCcorr_AocnBeta[t] = self.SLCaf_AocnBeta[t] + \ + self.SLCpov_AocnBeta[t] + \ + self.SLCden_AocnBeta[t] + + # same as above but adding the eustatic sea-level change term 'z0' + self.SLEaf_z0_AocnBeta[t] = (self.vaf_z0[t] / AocnBeta[t]) * \ + (rhoi / rhoo) + self.SLCaf_z0_AocnBeta[t] = -1.0 * (self.SLEaf_z0_AocnBeta[t] - + self.SLEaf_z0_AocnBeta[0]) + self.SLCpov_z0_AocnBeta[t] = -1.0 * (self.pov_z0[t] / AocnBeta[t] - + self.pov_z0[0] / AocnBeta[t]) + self.SLCcorr_z0_AocnBeta[t] = self.SLCaf_z0_AocnBeta[t] + \ + self.SLCpov_z0_AocnBeta[t] + \ + self.SLCden_AocnBeta[t] + + # get total ice mass change at each time step + self.dgrnd_vol[t] = (self.grnd_vol[t] - self.grnd_vol[0]) + self.dgrnd_vol_unscaled[t] = (self.grnd_vol_unscaled[t] - + self.grnd_vol_unscaled[0]) + + DS.close() + + # calculate RMSE between MALI and SLM calculation of SLC + if coupling: + self.diff_slc = self.SLC_slm_AocnBeta - self.SLCcorr_AocnBeta + self.diff_slc_z0 = self.SLC_slm_AocnBeta - \ + self.SLCcorr_z0_AocnBeta + self.rmse_slc = np.sqrt((self.diff_slc**2).mean()) + self.pct_err_slc = ((abs(self.diff_slc[1:]) / + abs(self.SLCcorr_AocnBeta[1:])) * + 100).mean() + self.pct_err_slc_z0 = ((abs(self.diff_slc_z0[1:]) / + abs(self.SLCcorr_z0_AocnBeta[1:])) * + 100).mean() + self.pct_err_slm = ((abs(self.diff_slc[1:]) / + abs(self.SLC_slm_AocnBeta[1:])) * + 100).mean() + self.pct_err_slm_z0 = ((abs(self.diff_slc_z0[1:]) / + abs(self.SLC_slm_AocnBeta[1:])) * + 100).mean() + + else: + logger.info('MALI is not coupled to the SLM. No MPE to compute.') + + +def grounded(cellMask): + # return ((cellMask&32)//32) & ~((cellMask&4)//4) + return ((cellMask & 32) // 32) * np.logical_not((cellMask & 4) // 4) + + +class slm_outputs: + """ + Read and calculate the SLM outputs + + Attributes + ---------- + data : float + SLM text-formatted outputs + + change_dt : float + Change in data value between each time step + + change_total : float + Change in data value at a time step + with respect to the intial time step + """ + def __init__(self, filename, nyrs): + """ + Calculate the sea-level model output file + + Parameters + ---------- + filename : str + Filename of SLM output data + + nyrs : int + Number of time indices to be analyzed + """ + self.fname = filename + self.data = np.loadtxt(self.fname)[0:nyrs] + self.change_dt = np.zeros(nyrs,) + self.change_total = np.zeros(nyrs,) + for i in range(len(self.change_dt) - 1): + self.change_dt[i + 1] = (float(self.data[i + 1]) - + float(self.data[i])) + self.change_total[i + 1] = (float(self.data[i + 1]) - + float(self.data[0])) diff --git a/compass/landice/tests/slm/namelist.sealevel.template b/compass/landice/tests/slm/namelist.sealevel.template new file mode 100644 index 0000000000..4676e62bdc --- /dev/null +++ b/compass/landice/tests/slm/namelist.sealevel.template @@ -0,0 +1,66 @@ +&time_config + itersl = 1 + starttime = 2015 + dt1 = 1 + +/ +&io_directory + inputfolder_ice = '/usr/projects/climate/hollyhan/SeaLevelModel_Inputs/icemodel/GL{{ nglv }}/ice_noGrIS_GL{{ nglv }}/' + inputfolder = '/usr/projects/climate/hollyhan/SeaLevelModel_Inputs/others/GL{{ nglv }}/' + planetfolder = '/usr/projects/climate/hollyhan/SeaLevelModel_Inputs/earthmodel/' + gridfolder = '/usr/projects/climate/hollyhan/SeaLevelModel_Inputs/others/GL{{ nglv }}/' + outputfolder = 'OUTPUT_SLM/' + outputfolder_ice = 'ICELOAD_SLM/' + folder_coupled = '' + +/ +&file_format + ext ='' + fType_in = 'text' + fType_out = 'both' + +/ +&file_name + planetmodel = 'prem_512.l60K2C.sum18p6.dum19p2.tz19p4.lm22' + icemodel = 'iceGlobalDomain_zeroField_GL{{ nglv }}_' + icemodel_out = 'iceload_out_' + timearray = 'times' + topomodel = 'idealized_topo_1500mOcean_SH_GL{{ nglv }}' + topo_initial = 'idealized_topo_1500mOcean_SH_GL{{ nglv }}' + 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. + +/ +&model_resolution + norder = 512 + nglv = {{ nglv }} + +/ +&timewindow_config + L_sim = 100 + dt1 = 1 + dt2 = 10 + dt3 = 10 + dt4 = 10 + Ldt1 = 100 + Ldt2 = 0 + Ldt3 = 0 + Ldt4 = 0 + +/ +&others + whichplanet = 'earth' + + +/