diff --git a/compass/ocean/mesh/cull.py b/compass/ocean/mesh/cull.py index eb285839fa..baeae2e50a 100644 --- a/compass/ocean/mesh/cull.py +++ b/compass/ocean/mesh/cull.py @@ -308,7 +308,7 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages, has_remapped_topo = os.path.exists('unsmoothed_topography.nc') if with_cavities and not has_remapped_topo: - raise ValueError('Mesh culling with caviites must be from ' + raise ValueError('Mesh culling with cavities must be from ' 'remapped topography.') if has_remapped_topo: @@ -469,7 +469,9 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages, landIceMask = dsMask.regionCellMasks.isel(nRegions=0) dsLandIceMask = xr.Dataset() dsLandIceMask['landIceMask'] = landIceMask - dsLandIceMask['landIceFloatingMask'] = landIceMask + dsLandIceMask['landIceFloatingMask'] = xr.logical_and( + landIceMask, dsMask.landIceFloatingFraction > 0. + ) write_netcdf(dsLandIceMask, 'land_ice_mask.nc') diff --git a/compass/ocean/mesh/remap_topography.py b/compass/ocean/mesh/remap_topography.py index a02d8a30e1..9aa1f5ad22 100644 --- a/compass/ocean/mesh/remap_topography.py +++ b/compass/ocean/mesh/remap_topography.py @@ -129,9 +129,10 @@ def run(self): Run this step of the test case """ super().run() - self._symlink_unsmoothed() + if not self.smoothing: + self._symlink_unsmoothed() if self.symlinked_to_unsmoothed: - # we symlinked to the unsmoothed topography and we're done! + print('we symlinked to the unsmoothed topography and we\'re done!') return config = self.config @@ -362,6 +363,8 @@ def _modify_remapped_bathymetry(self): for var in ['bed_elevation', 'landIceThkObserved']: ds_out[var] = xr.where(valid, ds_out[var] / norm, 0.) + # Note: the following variables are not used when MALI topography is + # used, they derive from the MALI topography file thickness = ds_out.landIceThkObserved bed = ds_out.bed_elevation flotation_thickness = - (ocean_density / ice_density) * bed diff --git a/compass/ocean/plot.py b/compass/ocean/plot.py index ae5eb91408..03c605356d 100644 --- a/compass/ocean/plot.py +++ b/compass/ocean/plot.py @@ -47,11 +47,14 @@ def plot_initial_state(input_file_name='initial_state.nc', varName = 'maxLevelCell' var = ds[varName] maxLevelCell = var.values - 1 + if 'nVertLevels' in var.sizes: + for iCell in range(var.sizes['nCells']): + var[maxLevelCell[iCell]:, iCell] = np.nan xarray.plot.hist(var, bins=nVertLevels - 4) plt.ylabel('frequency') plt.xlabel(varName) - txt = '{}{:9.2e} {:9.2e} {}\n'.format(txt, var.min().values, - var.max().values, varName) + txt = '{}{:9.2e} {:9.2e} {}\n'.format(txt, np.nanmin(var.values), + np.nanmax(var.values), varName) plt.subplot(4, 3, 3) varName = 'bottomDepth' diff --git a/compass/ocean/tests/global_ocean/__init__.py b/compass/ocean/tests/global_ocean/__init__.py index da065fa122..88e649e449 100644 --- a/compass/ocean/tests/global_ocean/__init__.py +++ b/compass/ocean/tests/global_ocean/__init__.py @@ -39,6 +39,8 @@ def __init__(self, mpas_core): # for other meshes, we do fewer tests self._add_tests(mesh_names=['QU', 'Icos', 'QUwISC', 'IcoswISC']) + self._add_tests(mesh_names=['IcoswISC'], + mali_ais_topo='AIS_4to20km') self._add_tests(mesh_names=['EC30to60', 'ECwISC30to60']) diff --git a/compass/ocean/tests/global_ocean/init/namelist.wisc b/compass/ocean/tests/global_ocean/init/namelist.wisc index c668dce71b..edfd6994c2 100644 --- a/compass/ocean/tests/global_ocean/init/namelist.wisc +++ b/compass/ocean/tests/global_ocean/init/namelist.wisc @@ -3,3 +3,4 @@ config_init_vertical_grid_type = 'haney-number' config_global_ocean_depress_by_land_ice = .true. config_global_ocean_use_constant_land_ice_cavity_temperature = .true. config_global_ocean_constant_land_ice_cavity_temperature = -1.8 +config_rx1_max = 20.0 diff --git a/compass/ocean/tests/global_ocean/mesh/__init__.py b/compass/ocean/tests/global_ocean/mesh/__init__.py index 764848c7de..6558beed2a 100644 --- a/compass/ocean/tests/global_ocean/mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/mesh/__init__.py @@ -175,7 +175,7 @@ def __init__(self, test_group, mesh_name, # noqa: C901 mesh_name=mesh_name, smoothing=False, mali_ais_topo=mali_ais_topo, - ocean_includes_grounded=False) + ocean_includes_grounded=True) self.add_step(unsmoothed_topo) @@ -187,7 +187,7 @@ def __init__(self, test_group, mesh_name, # noqa: C901 smoothing=True, unsmoothed_topo=unsmoothed_topo, mali_ais_topo=mali_ais_topo, - ocean_includes_grounded=False) + ocean_includes_grounded=True) self.add_step(smoothed_topo) diff --git a/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py index 780cec959b..7573156bf9 100644 --- a/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py +++ b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py @@ -244,9 +244,9 @@ def _remap_mali_to_mpaso(self, map_filename): ds_remapped = xr.open_dataset('mali_topography_ncremap.nc') ds_out = xr.Dataset() - for var_name in ds_remapped: + for var_name in ds_remapped.keys(): var = ds_remapped[var_name] - if 'Frac' in var: + if 'Frac' in var_name: # copy the fractional variable, making sure it doesn't # exceed 1 var = np.minimum(var, 1.) @@ -315,7 +315,9 @@ def _combine_topo(self): alpha = ds_out.maliFrac # NOTE: MALI's ocean fraction is already scaled by the MALI fraction - ds_out['oceanFracObserved'] = ( + ds_out['oceanFracObserved'] = xr.where( + ds_mali.bed_elevation > 0., + 0., ds_mali.oceanFrac + (1.0 - alpha) * ds_bedmachine.oceanFracObserved) diff --git a/compass/ocean/vertical/grid_1d/__init__.py b/compass/ocean/vertical/grid_1d/__init__.py index 89064fb8ec..5457bcb983 100644 --- a/compass/ocean/vertical/grid_1d/__init__.py +++ b/compass/ocean/vertical/grid_1d/__init__.py @@ -139,7 +139,9 @@ def add_1d_grid(config, ds): """ interfaces = generate_1d_grid(config=config) - + plot_1d_grid = config.get('vertical_grid', 'plot_1d_grid') + if plot_1d_grid: + write_1d_grid(interfaces=interfaces, out_filename='vertical_grid.nc') ds['refTopDepth'] = ('nVertLevels', interfaces[0:-1]) ds['refZMid'] = ('nVertLevels', -0.5 * (interfaces[1:] + interfaces[0:-1])) ds['refBottomDepth'] = ('nVertLevels', interfaces[1:])