diff --git a/E3SM-Project b/E3SM-Project
index b4d5b10600..c292bec000 160000
--- a/E3SM-Project
+++ b/E3SM-Project
@@ -1 +1 @@
-Subproject commit b4d5b106009a7c6b7bf7a28755727e9cedab4a1f
+Subproject commit c292bec000f10a1a8935f8c7ab3b54f0717ba0f8
diff --git a/compass/ocean/cached_files.json b/compass/ocean/cached_files.json
index a8ac7f8d31..f6090ade12 100644
--- a/compass/ocean/cached_files.json
+++ b/compass/ocean/cached_files.json
@@ -8,10 +8,10 @@
"ocean/global_ocean/QUwISC240/mesh/mesh/culled_mesh.nc": "global_ocean/QUwISC240/mesh/mesh/culled_mesh.210809.nc",
"ocean/global_ocean/QUwISC240/mesh/mesh/culled_graph.info": "global_ocean/QUwISC240/mesh/mesh/culled_graph.210809.info",
"ocean/global_ocean/QUwISC240/mesh/mesh/critical_passages_mask_final.nc": "global_ocean/QUwISC240/mesh/mesh/critical_passages_mask_final.210809.nc",
- "ocean/global_ocean/QUwISC240/PHC/init/initial_state/initial_state.nc": "global_ocean/QUwISC240/PHC/init/initial_state/initial_state.230327.nc",
- "ocean/global_ocean/QUwISC240/PHC/init/initial_state/init_mode_forcing_data.nc": "global_ocean/QUwISC240/PHC/init/initial_state/init_mode_forcing_data.230327.nc",
- "ocean/global_ocean/QUwISC240/PHC/init/initial_state/graph.info": "global_ocean/QUwISC240/PHC/init/initial_state/graph.230327.info",
- "ocean/global_ocean/QUwISC240/PHC/init/ssh_adjustment/adjusted_init.nc": "global_ocean/QUwISC240/PHC/init/ssh_adjustment/adjusted_init.230327.nc",
+ "ocean/global_ocean/QUwISC240/PHC/init/initial_state/initial_state.nc": "global_ocean/QUwISC240/PHC/init/initial_state/initial_state.230414.nc",
+ "ocean/global_ocean/QUwISC240/PHC/init/initial_state/init_mode_forcing_data.nc": "global_ocean/QUwISC240/PHC/init/initial_state/init_mode_forcing_data.230414.nc",
+ "ocean/global_ocean/QUwISC240/PHC/init/initial_state/graph.info": "global_ocean/QUwISC240/PHC/init/initial_state/graph.230414.info",
+ "ocean/global_ocean/QUwISC240/PHC/init/ssh_adjustment/adjusted_init.nc": "global_ocean/QUwISC240/PHC/init/ssh_adjustment/adjusted_init.230414.nc",
"ocean/global_ocean/EC30to60/mesh/mesh/culled_mesh.nc": "global_ocean/EC30to60/mesh/mesh/culled_mesh.210809.nc",
"ocean/global_ocean/EC30to60/mesh/mesh/culled_graph.info": "global_ocean/EC30to60/mesh/mesh/culled_graph.210809.info",
"ocean/global_ocean/EC30to60/mesh/mesh/critical_passages_mask_final.nc": "global_ocean/EC30to60/mesh/mesh/critical_passages_mask_final.210809.nc",
@@ -27,10 +27,10 @@
"ocean/global_ocean/ECwISC30to60/mesh/mesh/culled_mesh.nc": "global_ocean/ECwISC30to60/mesh/mesh/culled_mesh.210809.nc",
"ocean/global_ocean/ECwISC30to60/mesh/mesh/culled_graph.info": "global_ocean/ECwISC30to60/mesh/mesh/culled_graph.210809.info",
"ocean/global_ocean/ECwISC30to60/mesh/mesh/critical_passages_mask_final.nc": "global_ocean/ECwISC30to60/mesh/mesh/critical_passages_mask_final.210809.nc",
- "ocean/global_ocean/ECwISC30to60/PHC/init/initial_state/initial_state.nc": "global_ocean/ECwISC30to60/PHC/init/initial_state/initial_state.230411.nc",
- "ocean/global_ocean/ECwISC30to60/PHC/init/initial_state/init_mode_forcing_data.nc": "global_ocean/ECwISC30to60/PHC/init/initial_state/init_mode_forcing_data.230411.nc",
- "ocean/global_ocean/ECwISC30to60/PHC/init/ssh_adjustment/adjusted_init.nc": "global_ocean/ECwISC30to60/PHC/init/ssh_adjustment/adjusted_init.230411.nc",
- "ocean/global_ocean/ECwISC30to60/PHC/init/initial_state/graph.info": "global_ocean/ECwISC30to60/PHC/init/initial_state/graph.230411.info",
+ "ocean/global_ocean/ECwISC30to60/PHC/init/initial_state/initial_state.nc": "global_ocean/ECwISC30to60/PHC/init/initial_state/initial_state.230414.nc",
+ "ocean/global_ocean/ECwISC30to60/PHC/init/initial_state/init_mode_forcing_data.nc": "global_ocean/ECwISC30to60/PHC/init/initial_state/init_mode_forcing_data.230414.nc",
+ "ocean/global_ocean/ECwISC30to60/PHC/init/ssh_adjustment/adjusted_init.nc": "global_ocean/ECwISC30to60/PHC/init/ssh_adjustment/adjusted_init.230414.nc",
+ "ocean/global_ocean/ECwISC30to60/PHC/init/initial_state/graph.info": "global_ocean/ECwISC30to60/PHC/init/initial_state/graph.230414.info",
"ocean/global_ocean/SOwISC12to60/mesh/mesh/culled_mesh.nc": "global_ocean/SOwISC12to60/mesh/mesh/culled_mesh.210810.nc",
"ocean/global_ocean/SOwISC12to60/mesh/mesh/culled_graph.info": "global_ocean/SOwISC12to60/mesh/mesh/culled_graph.210810.info",
"ocean/global_ocean/SOwISC12to60/mesh/mesh/critical_passages_mask_final.nc": "global_ocean/SOwISC12to60/mesh/mesh/critical_passages_mask_final.210810.nc",
@@ -118,14 +118,14 @@
"ocean/global_ocean/QU240/mesh/cull_mesh/culled_mesh.nc": "global_ocean/QU240/mesh/cull_mesh/culled_mesh.220814.nc",
"ocean/global_ocean/QU240/mesh/cull_mesh/culled_graph.info": "global_ocean/QU240/mesh/cull_mesh/culled_graph.220814.info",
"ocean/global_ocean/QU240/mesh/cull_mesh/critical_passages_mask_final.nc": "global_ocean/QU240/mesh/cull_mesh/critical_passages_mask_final.220814.nc",
- "ocean/global_ocean/QUwISC240/mesh/base_mesh/mesh.msh": "global_ocean/QUwISC240/mesh/base_mesh/mesh.230327.msh",
- "ocean/global_ocean/QUwISC240/mesh/base_mesh/base_mesh.nc": "global_ocean/QUwISC240/mesh/base_mesh/base_mesh.230327.nc",
- "ocean/global_ocean/QUwISC240/mesh/base_mesh/cellWidthVsLatLon.nc": "global_ocean/QUwISC240/mesh/base_mesh/cellWidthVsLatLon.230327.nc",
- "ocean/global_ocean/QUwISC240/mesh/base_mesh/graph.info": "global_ocean/QUwISC240/mesh/base_mesh/graph.230327.info",
- "ocean/global_ocean/QUwISC240/mesh/cull_mesh/culled_mesh.nc": "global_ocean/QUwISC240/mesh/cull_mesh/culled_mesh.230327.nc",
- "ocean/global_ocean/QUwISC240/mesh/cull_mesh/culled_graph.info": "global_ocean/QUwISC240/mesh/cull_mesh/culled_graph.230327.info",
- "ocean/global_ocean/QUwISC240/mesh/cull_mesh/critical_passages_mask_final.nc": "global_ocean/QUwISC240/mesh/cull_mesh/critical_passages_mask_final.230327.nc",
- "ocean/global_ocean/QUwISC240/mesh/cull_mesh/land_ice_mask.nc": "global_ocean/QUwISC240/mesh/cull_mesh/land_ice_mask.230327.nc",
+ "ocean/global_ocean/QUwISC240/mesh/base_mesh/mesh.msh": "global_ocean/QUwISC240/mesh/base_mesh/mesh.230414.msh",
+ "ocean/global_ocean/QUwISC240/mesh/base_mesh/base_mesh.nc": "global_ocean/QUwISC240/mesh/base_mesh/base_mesh.230414.nc",
+ "ocean/global_ocean/QUwISC240/mesh/base_mesh/cellWidthVsLatLon.nc": "global_ocean/QUwISC240/mesh/base_mesh/cellWidthVsLatLon.230414.nc",
+ "ocean/global_ocean/QUwISC240/mesh/base_mesh/graph.info": "global_ocean/QUwISC240/mesh/base_mesh/graph.230414.info",
+ "ocean/global_ocean/QUwISC240/mesh/cull_mesh/culled_mesh.nc": "global_ocean/QUwISC240/mesh/cull_mesh/culled_mesh.230414.nc",
+ "ocean/global_ocean/QUwISC240/mesh/cull_mesh/culled_graph.info": "global_ocean/QUwISC240/mesh/cull_mesh/culled_graph.230414.info",
+ "ocean/global_ocean/QUwISC240/mesh/cull_mesh/critical_passages_mask_final.nc": "global_ocean/QUwISC240/mesh/cull_mesh/critical_passages_mask_final.230414.nc",
+ "ocean/global_ocean/QUwISC240/mesh/cull_mesh/land_ice_mask.nc": "global_ocean/QUwISC240/mesh/cull_mesh/land_ice_mask.230414.nc",
"ocean/global_ocean/EC30to60/mesh/base_mesh/mesh.msh": "global_ocean/EC30to60/mesh/base_mesh/mesh.230411.msh",
"ocean/global_ocean/EC30to60/mesh/base_mesh/base_mesh.nc": "global_ocean/EC30to60/mesh/base_mesh/base_mesh.230411.nc",
"ocean/global_ocean/EC30to60/mesh/base_mesh/cellWidthVsLatLon.nc": "global_ocean/EC30to60/mesh/base_mesh/cellWidthVsLatLon.230411.nc",
@@ -133,14 +133,14 @@
"ocean/global_ocean/EC30to60/mesh/cull_mesh/culled_mesh.nc": "global_ocean/EC30to60/mesh/cull_mesh/culled_mesh.230411.nc",
"ocean/global_ocean/EC30to60/mesh/cull_mesh/culled_graph.info": "global_ocean/EC30to60/mesh/cull_mesh/culled_graph.230411.info",
"ocean/global_ocean/EC30to60/mesh/cull_mesh/critical_passages_mask_final.nc": "global_ocean/EC30to60/mesh/cull_mesh/critical_passages_mask_final.230411.nc",
- "ocean/global_ocean/ECwISC30to60/mesh/base_mesh/mesh.msh": "global_ocean/ECwISC30to60/mesh/base_mesh/mesh.230411.msh",
- "ocean/global_ocean/ECwISC30to60/mesh/base_mesh/base_mesh.nc": "global_ocean/ECwISC30to60/mesh/base_mesh/base_mesh.230411.nc",
- "ocean/global_ocean/ECwISC30to60/mesh/base_mesh/cellWidthVsLatLon.nc": "global_ocean/ECwISC30to60/mesh/base_mesh/cellWidthVsLatLon.230411.nc",
- "ocean/global_ocean/ECwISC30to60/mesh/base_mesh/graph.info": "global_ocean/ECwISC30to60/mesh/base_mesh/graph.230411.info",
- "ocean/global_ocean/ECwISC30to60/mesh/cull_mesh/culled_mesh.nc": "global_ocean/ECwISC30to60/mesh/cull_mesh/culled_mesh.230411.nc",
- "ocean/global_ocean/ECwISC30to60/mesh/cull_mesh/culled_graph.info": "global_ocean/ECwISC30to60/mesh/cull_mesh/culled_graph.230411.info",
- "ocean/global_ocean/ECwISC30to60/mesh/cull_mesh/critical_passages_mask_final.nc": "global_ocean/ECwISC30to60/mesh/cull_mesh/critical_passages_mask_final.230411.nc",
- "ocean/global_ocean/ECwISC30to60/mesh/cull_mesh/land_ice_mask.nc": "global_ocean/ECwISC30to60/mesh/cull_mesh/land_ice_mask.230411.nc",
+ "ocean/global_ocean/ECwISC30to60/mesh/base_mesh/mesh.msh": "global_ocean/ECwISC30to60/mesh/base_mesh/mesh.230414.msh",
+ "ocean/global_ocean/ECwISC30to60/mesh/base_mesh/base_mesh.nc": "global_ocean/ECwISC30to60/mesh/base_mesh/base_mesh.230414.nc",
+ "ocean/global_ocean/ECwISC30to60/mesh/base_mesh/cellWidthVsLatLon.nc": "global_ocean/ECwISC30to60/mesh/base_mesh/cellWidthVsLatLon.230414.nc",
+ "ocean/global_ocean/ECwISC30to60/mesh/base_mesh/graph.info": "global_ocean/ECwISC30to60/mesh/base_mesh/graph.230414.info",
+ "ocean/global_ocean/ECwISC30to60/mesh/cull_mesh/culled_mesh.nc": "global_ocean/ECwISC30to60/mesh/cull_mesh/culled_mesh.230414.nc",
+ "ocean/global_ocean/ECwISC30to60/mesh/cull_mesh/culled_graph.info": "global_ocean/ECwISC30to60/mesh/cull_mesh/culled_graph.230414.info",
+ "ocean/global_ocean/ECwISC30to60/mesh/cull_mesh/critical_passages_mask_final.nc": "global_ocean/ECwISC30to60/mesh/cull_mesh/critical_passages_mask_final.230414.nc",
+ "ocean/global_ocean/ECwISC30to60/mesh/cull_mesh/land_ice_mask.nc": "global_ocean/ECwISC30to60/mesh/cull_mesh/land_ice_mask.230414.nc",
"ocean/global_ocean/SOwISC12to60/mesh/base_mesh/mesh.msh": "global_ocean/SOwISC12to60/mesh/base_mesh/mesh.230327.msh",
"ocean/global_ocean/SOwISC12to60/mesh/base_mesh/base_mesh.nc": "global_ocean/SOwISC12to60/mesh/base_mesh/base_mesh.230327.nc",
"ocean/global_ocean/SOwISC12to60/mesh/base_mesh/cellWidthVsLatLon.nc": "global_ocean/SOwISC12to60/mesh/base_mesh/cellWidthVsLatLon.230327.nc",
@@ -151,6 +151,6 @@
"ocean/global_ocean/SOwISC12to60/mesh/cull_mesh/land_ice_mask.nc": "global_ocean/SOwISC12to60/mesh/cull_mesh/land_ice_mask.230327.nc",
"ocean/global_ocean/EC30to60/mesh/remap_topography/topography_remapped.nc": "global_ocean/EC30to60/mesh/remap_topography/topography_remapped.230411.nc",
"ocean/global_ocean/EC30to60/mesh/cull_mesh/topography_culled.nc": "global_ocean/EC30to60/mesh/cull_mesh/topography_culled.230411.nc",
- "ocean/global_ocean/ECwISC30to60/mesh/remap_topography/topography_remapped.nc": "global_ocean/ECwISC30to60/mesh/remap_topography/topography_remapped.230411.nc",
- "ocean/global_ocean/ECwISC30to60/mesh/cull_mesh/topography_culled.nc": "global_ocean/ECwISC30to60/mesh/cull_mesh/topography_culled.230411.nc"
+ "ocean/global_ocean/ECwISC30to60/mesh/remap_topography/topography_remapped.nc": "global_ocean/ECwISC30to60/mesh/remap_topography/topography_remapped.230414.nc",
+ "ocean/global_ocean/ECwISC30to60/mesh/cull_mesh/topography_culled.nc": "global_ocean/ECwISC30to60/mesh/cull_mesh/topography_culled.230414.nc"
}
diff --git a/compass/ocean/mesh/cull.py b/compass/ocean/mesh/cull.py
index f126e1398b..151cf05982 100644
--- a/compass/ocean/mesh/cull.py
+++ b/compass/ocean/mesh/cull.py
@@ -404,6 +404,7 @@ 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
write_netcdf(dsLandIceMask, 'land_ice_mask.nc')
diff --git a/compass/ocean/tests/global_ocean/init/streams.init b/compass/ocean/tests/global_ocean/init/streams.init
index 342d69ca0f..c105b6041e 100644
--- a/compass/ocean/tests/global_ocean/init/streams.init
+++ b/compass/ocean/tests/global_ocean/init/streams.init
@@ -32,7 +32,9 @@
+
+
diff --git a/compass/ocean/tests/global_ocean/init/streams.wisc b/compass/ocean/tests/global_ocean/init/streams.wisc
index f747f41b5e..ac5434d5ae 100644
--- a/compass/ocean/tests/global_ocean/init/streams.wisc
+++ b/compass/ocean/tests/global_ocean/init/streams.wisc
@@ -6,6 +6,7 @@
type = "input">
+
diff --git a/compass/ocean/tests/global_ocean/performance_test/streams.wisc b/compass/ocean/tests/global_ocean/performance_test/streams.wisc
index e03cfd8cd1..e9742d2284 100644
--- a/compass/ocean/tests/global_ocean/performance_test/streams.wisc
+++ b/compass/ocean/tests/global_ocean/performance_test/streams.wisc
@@ -14,7 +14,9 @@
+
+
diff --git a/compass/ocean/tests/ice_shelf_2d/initial_state.py b/compass/ocean/tests/ice_shelf_2d/initial_state.py
index cc7e29fac8..ac795c40ed 100644
--- a/compass/ocean/tests/ice_shelf_2d/initial_state.py
+++ b/compass/ocean/tests/ice_shelf_2d/initial_state.py
@@ -98,6 +98,8 @@ def run(self):
dim='Time', axis=0)
landIceFraction = modify_mask.astype(float)
landIceMask = modify_mask.copy()
+ landIceFloatingFraction = landIceFraction.copy()
+ landIceFloatingMask = landIceMask.copy()
ref_density = constants['SHR_CONST_RHOSW']
landIcePressure, landIceDraft = compute_land_ice_pressure_and_draft(
@@ -121,7 +123,9 @@ def run(self):
ds['fVertex'] = xarray.zeros_like(ds.xVertex)
ds['modifyLandIcePressureMask'] = modify_mask
ds['landIceFraction'] = landIceFraction
+ ds['landIceFloatingFraction'] = landIceFloatingFraction
ds['landIceMask'] = landIceMask
+ ds['landIceFloatingMask'] = landIceFloatingMask
ds['landIcePressure'] = landIcePressure
ds['landIceDraft'] = landIceDraft
diff --git a/compass/ocean/tests/isomip_plus/cull_mesh.py b/compass/ocean/tests/isomip_plus/cull_mesh.py
index e2850d2ece..a53bbd6fea 100644
--- a/compass/ocean/tests/isomip_plus/cull_mesh.py
+++ b/compass/ocean/tests/isomip_plus/cull_mesh.py
@@ -2,10 +2,7 @@
from mpas_tools.io import write_netcdf
from mpas_tools.mesh.conversion import convert, cull
-from compass.ocean.tests.isomip_plus.geom import (
- define_thin_film_mask_step1,
- interpolate_ocean_mask,
-)
+from compass.ocean.tests.isomip_plus.geom import interpolate_ocean_mask
from compass.ocean.tests.isomip_plus.spherical_mesh import add_isomip_plus_xy
from compass.step import Step
@@ -70,11 +67,8 @@ def run(self):
ds_mesh = xr.open_dataset('base_mesh.nc')
ds_geom = xr.open_dataset('input_geometry_processed.nc')
- if thin_film_present:
- ds_mask = define_thin_film_mask_step1(ds_mesh, ds_geom)
- else:
- ds_mask = \
- interpolate_ocean_mask(ds_mesh, ds_geom, min_ocean_fraction)
+ ds_mask = interpolate_ocean_mask(ds_mesh, ds_geom, min_ocean_fraction,
+ thin_film_present)
ds_mesh = cull(ds_mesh, dsInverse=ds_mask, logger=logger)
ds_mesh.attrs['is_periodic'] = 'NO'
diff --git a/compass/ocean/tests/isomip_plus/forward.py b/compass/ocean/tests/isomip_plus/forward.py
index d457f0798a..c1ea540438 100644
--- a/compass/ocean/tests/isomip_plus/forward.py
+++ b/compass/ocean/tests/isomip_plus/forward.py
@@ -2,6 +2,7 @@
import shutil
import time
+import numpy as np
import xarray
from compass.model import run_model
@@ -165,6 +166,8 @@ def run(self):
dsMesh=dsMesh, ds=ds, expt=self.experiment,
showProgress=show_progress)
+ plotter.plot_horiz_series(ds.ssh, 'ssh', 'ssh',
+ True, vmin=-700, vmax=0)
plotter.plot_3d_field_top_bot_section(
ds.temperature, nameInTitle='temperature', prefix='temp',
units='C', vmin=-2., vmax=1., cmap='cmo.thermal')
@@ -172,11 +175,49 @@ def run(self):
plotter.plot_3d_field_top_bot_section(
ds.salinity, nameInTitle='salinity', prefix='salin',
units='PSU', vmin=33.8, vmax=34.7, cmap='cmo.haline')
+ tol = 1e-10
+ dsIce = xarray.open_dataset(
+ os.path.join(self.work_dir,
+ 'land_ice_fluxes.nc'))
+ plotter.plot_horiz_series(
+ dsIce.topDragMagnitude,
+ 'topDragMagnitude', 'topDragMagnitude', True,
+ vmin=0 + tol, vmax=np.max(dsIce.topDragMagnitude.values),
+ cmap_set_under='k')
+ plotter.plot_horiz_series(
+ dsIce.landIceHeatFlux,
+ 'landIceHeatFlux', 'landIceHeatFlux', True,
+ vmin=np.min(dsIce.landIceHeatFlux.values),
+ vmax=np.max(dsIce.landIceHeatFlux.values))
+ plotter.plot_horiz_series(
+ dsIce.landIceInterfaceTemperature,
+ 'landIceInterfaceTemperature', 'landIceInterfaceTemperature',
+ True,
+ vmin=np.min(dsIce.landIceInterfaceTemperature.values),
+ vmax=np.max(dsIce.landIceInterfaceTemperature.values))
+ plotter.plot_horiz_series(
+ dsIce.landIceFreshwaterFlux,
+ 'landIceFreshwaterFlux', 'landIceFreshwaterFlux', True,
+ vmin=0 + tol, vmax=1e-4,
+ cmap_set_under='k', cmap_scale='log')
+ plotter.plot_horiz_series(
+ dsIce.landIceFraction,
+ 'landIceFraction', 'landIceFraction', True,
+ vmin=0 + tol, vmax=1 - tol,
+ cmap='cmo.balance',
+ cmap_set_under='k', cmap_set_over='r')
+ if 'landIceFloatingFraction' in dsIce.keys():
+ plotter.plot_horiz_series(
+ dsIce.landIceFloatingFraction,
+ 'landIceFloatingFraction', 'landIceFloatingFraction',
+ True, vmin=0 + tol, vmax=1 - tol,
+ cmap='cmo.balance', cmap_set_under='k', cmap_set_over='r')
if self.name == 'simulation':
- update_evaporation_flux(in_forcing_file='forcing_data_init.nc',
- out_forcing_file='forcing_data_updated.nc',
- out_forcing_link='forcing_data.nc')
+ update_evaporation_flux(
+ in_forcing_file='forcing_data_init.nc',
+ out_forcing_file='forcing_data_updated.nc',
+ out_forcing_link='forcing_data.nc')
replacements = {'config_do_restart': '.true.',
'config_start_time': "'file'"}
diff --git a/compass/ocean/tests/isomip_plus/geom.py b/compass/ocean/tests/isomip_plus/geom.py
index 01ce9063ae..ce761a1582 100755
--- a/compass/ocean/tests/isomip_plus/geom.py
+++ b/compass/ocean/tests/isomip_plus/geom.py
@@ -3,46 +3,8 @@
from mpas_tools.mesh.interpolation import interp_bilin
-def define_thin_film_mask_step1(ds_mesh, ds_geom):
- """
- Define an MPAS mesh mask for the ocean domain including cells over the
- full x- and y-range in order to include all land cells in the ocean's
- thin-film region.
-
- Parameters
- ----------
- ds_mesh : xarray.Dataset
- An MPAS-Ocean mesh
-
- ds_geom : xarray.Dataset
- Ice-sheet topography produced by
- :py:class:`compass.ocean.tests.isomip_plus.process_geom.ProcessGeom`
-
- Returns
- -------
- ds_mask : xarray.Dataset
- A dataset containing ``regionCellMasks``, a field with the ocean mask
- that can be used to cull land cells from the mesh
- """
- x, y, x_cell, y_cell, ocean_fraction = _get_geom_fields(ds_geom, ds_mesh)
-
- ds_mask = xarray.Dataset()
-
- valid = numpy.logical_and(
- numpy.logical_and(x_cell >= x[0], x_cell <= x[-1]),
- numpy.logical_and(y_cell >= y[0], y_cell <= y[-1]))
-
- mask = valid
-
- n_cells = mask.shape[0]
-
- ds_mask['regionCellMasks'] = (('nCells', 'nRegions'),
- mask.astype(int).reshape(n_cells, 1))
-
- return ds_mask
-
-
-def interpolate_ocean_mask(ds_mesh, ds_geom, min_ocean_fraction):
+def interpolate_ocean_mask(ds_mesh, ds_geom, min_ocean_fraction,
+ thin_film_present):
"""
Interpolate the ocean mask from the original BISICLES grid to the MPAS
mesh. This is handled separately from other fields because the ocean mask
@@ -62,13 +24,17 @@ def interpolate_ocean_mask(ds_mesh, ds_geom, min_ocean_fraction):
The minimum ocean fraction after interpolation, below which the cell
is masked as land (which is not distinguished from grounded ice)
+ thin_film_present: bool
+ Whether domain contains a thin film below grounded ice
+
Returns
-------
ds_mask : xarray.Dataset
A dataset containing ``regionCellMasks``, a field with the ocean mask
that can be used to cull land cells from the mesh
"""
- x, y, x_cell, y_cell, ocean_fraction = _get_geom_fields(ds_geom, ds_mesh)
+ x, y, x_cell, y_cell, ocean_fraction = _get_geom_fields(
+ ds_geom, ds_mesh, thin_film_present)
ds_mask = xarray.Dataset()
@@ -125,36 +91,43 @@ def interpolate_geom(ds_mesh, ds_geom, min_ocean_fraction, thin_film_present):
* ``oceanFracObserved`` -- the fraction of the mesh cell that is ocean
* ``landIceFraction`` -- the fraction of the mesh cell that is covered
- by an ice shelf
+ by land ice, grounded or floating
+
+ * ``landIceFloatingFraction`` -- the fraction of the mesh cell that is
+ covered by an ice shelf
* ``smoothedDraftMask`` -- a smoothed version of the floating mask that
may be useful for determining where to alter the vertical coordinate
to accommodate ice-shelf cavities
"""
- x, y, x_cell, y_cell, ocean_fraction = _get_geom_fields(ds_geom, ds_mesh)
+ x, y, x_cell, y_cell, ocean_fraction = _get_geom_fields(
+ ds_geom, ds_mesh, thin_film_present)
ds_out = xarray.Dataset(ds_mesh)
ds_out.attrs = ds_mesh.attrs
ds_geom['oceanFraction'] = ocean_fraction
+ ds_geom['landIceFraction'] = (ds_geom.landIceFloatingFraction +
+ ds_geom.landIceGroundedFraction)
- # mash the topography to the ocean region before interpolation
- if not thin_film_present:
- for var in ['Z_bed', 'Z_ice_draft', 'floatingIceFraction',
- 'smoothedDraftMask']:
- ds_geom[var] = ds_geom[var] * ds_geom['oceanFraction']
+ # mask the topography to the ocean region before interpolation
+ for var in ['Z_bed', 'Z_ice_draft', 'landIceFraction',
+ 'landIceFloatingFraction', 'smoothedDraftMask']:
+ ds_geom[var] = ds_geom[var] * ds_geom['oceanFraction']
fields = {'bottomDepthObserved': 'Z_bed',
'ssh': 'Z_ice_draft',
'oceanFracObserved': 'oceanFraction',
- 'landIceFraction': 'floatingIceFraction',
+ 'landIceFraction': 'landIceFraction',
+ 'landIceGroundedFraction': 'landIceGroundedFraction',
+ 'landIceFloatingFraction': 'landIceFloatingFraction',
'smoothedDraftMask': 'smoothedDraftMask'}
valid = numpy.logical_and(
numpy.logical_and(x_cell >= x[0], x_cell <= x[-1]),
numpy.logical_and(y_cell >= y[0], y_cell <= y[-1]))
- if not numpy.all(valid) and not thin_film_present:
+ if not numpy.all(valid):
raise ValueError('Something went wrong with culling. There are still '
'out-of-range cells in the culled mesh.')
@@ -163,25 +136,29 @@ def interpolate_geom(ds_mesh, ds_geom, min_ocean_fraction, thin_film_present):
ds_out[outfield] = (('nCells',), field)
ocean_frac_observed = ds_out['oceanFracObserved']
- if not numpy.all(ocean_frac_observed > min_ocean_fraction) and \
- not thin_film_present:
+ if not numpy.all(ocean_frac_observed > min_ocean_fraction):
raise ValueError('Something went wrong with culling. There are still '
'non-ocean cells in the culled mesh.')
- if not thin_film_present:
- for field in ['bottomDepthObserved', 'ssh', 'landIceFraction',
- 'smoothedDraftMask']:
- ds_out[field] = ds_out[field] / ocean_frac_observed
+ # Do not include landIceFloatingFraction here so that fractional melt
+ # fluxes are represented
+ for field in ['bottomDepthObserved', 'ssh', 'landIceFraction',
+ 'smoothedDraftMask']:
+ ds_out[field] = ds_out[field] / ocean_frac_observed
return ds_out
-def _get_geom_fields(ds_geom, ds_mesh):
+def _get_geom_fields(ds_geom, ds_mesh, thin_film_present):
x = ds_geom.x.values
y = ds_geom.y.values
x_cell = ds_mesh.xIsomipCell.values
y_cell = ds_mesh.yIsomipCell.values
- ocean_fraction = - ds_geom['landFraction'] + 1.0
+ if thin_film_present:
+ ocean_fraction = - ds_geom['landFraction'] + 1.0
+ else:
+ ocean_fraction = (ds_geom['landIceFloatingFraction'] +
+ ds_geom['openOceanFraction'])
return x, y, x_cell, y_cell, ocean_fraction
diff --git a/compass/ocean/tests/isomip_plus/initial_state.py b/compass/ocean/tests/isomip_plus/initial_state.py
index 8262871ee7..1160fc0761 100644
--- a/compass/ocean/tests/isomip_plus/initial_state.py
+++ b/compass/ocean/tests/isomip_plus/initial_state.py
@@ -1,6 +1,7 @@
import os
import shutil
+import cmocean # noqa: F401
import numpy as np
import xarray as xr
from mpas_tools.cime.constants import constants
@@ -105,13 +106,22 @@ def _compute_initial_condition(self):
ds['landIceFraction'] = \
ds['landIceFraction'].expand_dims(dim='Time', axis=0)
+ ds['landIceFloatingFraction'] = \
+ ds['landIceFloatingFraction'].expand_dims(dim='Time', axis=0)
ds['modifyLandIcePressureMask'] = \
(ds['landIceFraction'] > 0.01).astype(int)
- mask = ds.landIceFraction >= min_land_ice_fraction
+ # This inequality needs to be > rather than >= to ensure correctness
+ # when min_land_ice_fraction = 0
+ mask = ds.landIceFraction > min_land_ice_fraction
+
+ floating_mask = np.logical_and(
+ ds.landIceFloatingFraction > 0,
+ ds.landIceFraction > min_land_ice_fraction)
ds['landIceMask'] = mask.astype(int)
+ ds['landIceFloatingMask'] = floating_mask.astype(int)
ds['landIceFraction'] = xr.where(mask, ds.landIceFraction, 0.)
@@ -152,6 +162,20 @@ def _compute_initial_condition(self):
init_bot_temp = section.getfloat('init_bot_temp')
init_top_sal = section.getfloat('init_top_sal')
init_bot_sal = section.getfloat('init_bot_sal')
+ thin_film_mask = np.logical_and(mask.values,
+ np.logical_not(floating_mask.values))
+ # These coefficients are hard-coded as the defaults in the namelist
+ # Note that using the land ice pressure rather than the pressure at
+ # floatation will mean that there is a small amount of cooling from
+ # grounding line retreat. However, the thin film should be thin enough
+ # that this effect isn't signicant.
+ freezing_temp = (6.22e-2 +
+ -5.63e-2 * init_bot_sal +
+ -7.43e-8 * landIcePressure +
+ -1.74e-10 * landIcePressure * init_bot_sal)
+ _, thin_film_temp = np.meshgrid(ds.refZMid, freezing_temp)
+ _, thin_film_mask = np.meshgrid(ds.refZMid, thin_film_mask)
+ thin_film_temp = np.expand_dims(thin_film_temp, axis=0)
if self.vertical_coordinate == 'single_layer':
ds['temperature'] = init_bot_temp * xr.ones_like(frac)
ds['salinity'] = init_bot_sal * xr.ones_like(frac)
@@ -160,6 +184,9 @@ def _compute_initial_condition(self):
(1.0 - frac) * init_top_temp + frac * init_bot_temp
ds['salinity'] = \
(1.0 - frac) * init_top_sal + frac * init_bot_sal
+ # for thin film cells, set temperature to freezing point
+ ds['temperature'] = xr.where(thin_film_mask, thin_film_temp,
+ ds.temperature)
# compute coriolis
coriolis_parameter = section.getfloat('coriolis_parameter')
@@ -200,32 +227,68 @@ def _plot(self, ds):
sectionY=section_y, dsMesh=ds, ds=ds,
showProgress=show_progress)
+ ds['oceanFracObserved'] = \
+ ds['oceanFracObserved'].expand_dims(dim='Time', axis=0)
ds['landIcePressure'] = \
ds['landIcePressure'].expand_dims(dim='Time', axis=0)
+ ds['landIceGroundedFraction'] = \
+ ds['landIceGroundedFraction'].expand_dims(dim='Time', axis=0)
ds['bottomDepth'] = ds['bottomDepth'].expand_dims(dim='Time', axis=0)
ds['totalColThickness'] = ds['ssh']
ds['totalColThickness'].values = \
ds['layerThickness'].sum(dim='nVertLevels')
+ tol = 1e-10
+ plotter.plot_horiz_series(ds.landIceMask,
+ 'landIceMask', 'landIceMask',
+ True)
+ plotter.plot_horiz_series(ds.landIceFloatingMask,
+ 'landIceFloatingMask', 'landIceFloatingMask',
+ True)
plotter.plot_horiz_series(ds.landIcePressure,
'landIcePressure', 'landIcePressure',
- True)
+ True, vmin=1, vmax=1e6, cmap_scale='log')
plotter.plot_horiz_series(ds.ssh,
'ssh', 'ssh',
True, vmin=-700, vmax=0)
+ plotter.plot_horiz_series(ds.bottomDepth,
+ 'bottomDepth', 'bottomDepth',
+ True, vmin=0, vmax=700)
plotter.plot_horiz_series(ds.ssh + ds.bottomDepth,
'H', 'H', True,
- vmin=min_column_thickness + 1e-10, vmax=700,
+ vmin=min_column_thickness + tol, vmax=700,
cmap_set_under='r', cmap_scale='log')
plotter.plot_horiz_series(ds.totalColThickness,
'totalColThickness', 'totalColThickness',
True, vmin=min_column_thickness + 1e-10,
vmax=700, cmap_set_under='r')
+ plotter.plot_horiz_series(ds.landIceFraction,
+ 'landIceFraction', 'landIceFraction',
+ True, vmin=0 + tol, vmax=1 - tol,
+ cmap='cmo.balance',
+ cmap_set_under='k', cmap_set_over='r')
+ plotter.plot_horiz_series(ds.landIceFloatingFraction,
+ 'landIceFloatingFraction',
+ 'landIceFloatingFraction',
+ True, vmin=0 + tol, vmax=1 - tol,
+ cmap='cmo.balance',
+ cmap_set_under='k', cmap_set_over='r')
+ plotter.plot_horiz_series(ds.landIceGroundedFraction,
+ 'landIceGroundedFraction',
+ 'landIceGroundedFraction',
+ True, vmin=0 + tol, vmax=1 - tol,
+ cmap='cmo.balance',
+ cmap_set_under='k', cmap_set_over='r')
+ plotter.plot_horiz_series(ds.oceanFracObserved,
+ 'oceanFracObserved', 'oceanFracObserved',
+ True, vmin=0 + tol, vmax=1 - tol,
+ cmap='cmo.balance',
+ cmap_set_under='k', cmap_set_over='r')
plotter.plot_layer_interfaces()
plotter.plot_3d_field_top_bot_section(
ds.layerThickness, nameInTitle='layerThickness',
prefix='h', units='m',
- vmin=min_column_thickness + 1e-10, vmax=50,
+ vmin=min_column_thickness + tol, vmax=50,
cmap='cmo.deep_r', cmap_set_under='r')
plotter.plot_3d_field_top_bot_section(
@@ -318,11 +381,13 @@ def _write_time_varying_forcing(self, ds_init):
landIceDraft = list()
landIcePressure = list()
landIceFraction = list()
+ landIceFloatingFraction = list()
for scale in scales:
landIceDraft.append(scale * ds_init.landIceDraft)
landIcePressure.append(scale * ds_init.landIcePressure)
landIceFraction.append(ds_init.landIceFraction)
+ landIceFloatingFraction.append(ds_init.landIceFloatingFraction)
ds_out['landIceDraftForcing'] = xr.concat(landIceDraft, 'Time')
ds_out.landIceDraftForcing.attrs['units'] = 'm'
@@ -337,6 +402,10 @@ def _write_time_varying_forcing(self, ds_init):
xr.concat(landIceFraction, 'Time')
ds_out.landIceFractionForcing.attrs['long_name'] = \
'The fraction of each cell covered by land ice'
+ ds_out['landIceFloatingFractionForcing'] = \
+ xr.concat(landIceFloatingFraction, 'Time')
+ ds_out.landIceFloatingFractionForcing.attrs['long_name'] = \
+ 'The fraction of each cell covered by floating land ice'
write_netcdf(ds_out, 'land_ice_forcing.nc')
ds_init['landIceDraft'] = scales[0] * ds_init.landIceDraft
diff --git a/compass/ocean/tests/isomip_plus/misomip.py b/compass/ocean/tests/isomip_plus/misomip.py
index 943f36a74e..11983fe05e 100644
--- a/compass/ocean/tests/isomip_plus/misomip.py
+++ b/compass/ocean/tests/isomip_plus/misomip.py
@@ -686,7 +686,8 @@ def writeVar(varName, varField, varMask=None):
vars['time'][tIndex] = secPerDay * days
freshwaterFlux = inVars['timeMonthly_avg_landIceFreshwaterFlux'][0, :]
- inCavityFraction = inVars['timeMonthly_avg_landIceFraction'][0, :]
+ inCavityFraction = \
+ inVars['timeMonthly_avg_landIceFloatingFraction'][0, :]
outCavityFraction = interpHoriz(inCavityFraction)
outCavityMask = outCavityFraction > normalizationThreshold
meltRate = freshwaterFlux / rho_fw
diff --git a/compass/ocean/tests/isomip_plus/process_geom.py b/compass/ocean/tests/isomip_plus/process_geom.py
index 97823437b0..1293cc78cd 100644
--- a/compass/ocean/tests/isomip_plus/process_geom.py
+++ b/compass/ocean/tests/isomip_plus/process_geom.py
@@ -74,11 +74,13 @@ def run(self):
rename = {'upperSurface': 'Z_ice_surface',
'lowerSurface': 'Z_ice_draft',
'bedrockTopography': 'Z_bed',
- 'floatingMask': 'floatingIceFraction',
- 'groundedMask': 'landFraction',
+ 'floatingMask': 'landIceFloatingFraction',
+ 'groundedMask': 'landIceGroundedFraction',
'openOceanMask': 'openOceanFraction'}
ds_in = ds_in.rename(rename)
+ # This test case assumes that all land is ice covered
+ ds_in['landFraction'] = ds_in.landIceGroundedFraction
x_in = ds_in.x.values
y_in = ds_in.y.values
@@ -99,7 +101,8 @@ def run(self):
land_values = {'Z_ice_surface': 0.,
'Z_ice_draft': 0.,
'Z_bed': 0.,
- 'floatingIceFraction': 0.,
+ 'landIceFloatingFraction': 0.,
+ 'landIceGroundedFraction': 1.,
'landFraction': 1.,
'openOceanFraction': 0.}
@@ -113,10 +116,10 @@ def run(self):
ds['Z_ice_draft'] = draft_scaling * ds.Z_ice_draft
# take care of calving criterion
- mask = np.logical_or(ds.floatingIceFraction <= 0.1,
+ mask = np.logical_or(ds.landIceFloatingFraction <= 0.1,
ds.iceThickness >= min_ice_thickness)
- for var in ['Z_ice_surface', 'Z_ice_draft', 'floatingIceFraction']:
+ for var in ['Z_ice_surface', 'Z_ice_draft', 'landIceFloatingFraction']:
ds[var] = xr.where(mask, ds[var], 0.0)
ds['openOceanFraction'] = xr.where(mask, ds.openOceanFraction,
1. - ds.landFraction)
@@ -131,8 +134,8 @@ def run(self):
# copy attributes
for var in ['x', 'y', 'Z_ice_surface', 'Z_ice_draft', 'Z_bed',
- 'floatingIceFraction', 'landFraction',
- 'openOceanFraction']:
+ 'landIceFloatingFraction', 'landIceGroundedFraction',
+ 'landFraction', 'openOceanFraction']:
attrs = ds_in[var].attrs
if 'units' in attrs and attrs['units'] == 'unitless':
attrs.pop('units')
@@ -168,9 +171,8 @@ def _smooth_geometry(land_fraction, ds, filter_sigma, threshold=0.01):
mode='constant', cval=0.)
bed[mask] /= smoothed_mask[mask]
- smoothed_draft_mask = filters.gaussian_filter(ds.floatingIceFraction,
- filter_sigma,
- mode='constant', cval=0.)
+ smoothed_draft_mask = filters.gaussian_filter(
+ ds.landIceFloatingFraction, filter_sigma, mode='constant', cval=0.)
smoothed_draft_mask[mask] /= smoothed_mask[mask]
ds['Z_ice_draft'] = (('y', 'x'), draft)
diff --git a/compass/ocean/tests/isomip_plus/streams.forward.template b/compass/ocean/tests/isomip_plus/streams.forward.template
index b6b8d83938..4eff72b570 100644
--- a/compass/ocean/tests/isomip_plus/streams.forward.template
+++ b/compass/ocean/tests/isomip_plus/streams.forward.template
@@ -56,7 +56,9 @@
+
+
@@ -119,6 +121,8 @@
+
+
diff --git a/compass/ocean/tests/isomip_plus/viz/plot.py b/compass/ocean/tests/isomip_plus/viz/plot.py
index 73922c90d7..dc1e298573 100644
--- a/compass/ocean/tests/isomip_plus/viz/plot.py
+++ b/compass/ocean/tests/isomip_plus/viz/plot.py
@@ -466,7 +466,8 @@ def plot_haney_number(self, haneyFolder=None):
def plot_horiz_series(self, da, nameInTitle, prefix, oceanDomain,
units=None, vmin=None, vmax=None, cmap=None,
- cmap_set_under=None, cmap_scale='linear'):
+ cmap_set_under=None, cmap_set_over=None,
+ cmap_scale='linear'):
"""
Plot a series of image of a given variable
@@ -524,6 +525,7 @@ def plot_horiz_series(self, da, nameInTitle, prefix, oceanDomain,
oceanDomain=oceanDomain, vmin=vmin,
vmax=vmax, cmap=cmap,
cmap_set_under=cmap_set_under,
+ cmap_set_over=cmap_set_over,
cmap_scale=cmap_scale)
if self.showProgress:
bar.update(tIndex + 1)
@@ -532,7 +534,8 @@ def plot_horiz_series(self, da, nameInTitle, prefix, oceanDomain,
def plot_3d_field_top_bot_section(self, da, nameInTitle, prefix,
units=None, vmin=None, vmax=None,
- cmap=None, cmap_set_under=None):
+ cmap=None, cmap_set_under=None,
+ cmap_set_over=None):
"""
Plot a series of images of a given 3D variable showing the value
at the top (sea surface or ice-ocean interface), sea floor and in an
@@ -587,7 +590,8 @@ def plot_3d_field_top_bot_section(self, da, nameInTitle, prefix,
'top {}'.format(nameInTitle),
'top{}'.format(prefix), oceanDomain=True,
vmin=vmin, vmax=vmax, cmap=cmap,
- cmap_set_under=cmap_set_under)
+ cmap_set_under=cmap_set_under,
+ cmap_set_over=cmap_set_over)
maxLevelCell = self.dsMesh.maxLevelCell - 1
@@ -770,7 +774,8 @@ def update_date(self, tIndex):
def _plot_horiz_field(self, field, title, outFileName, oceanDomain=True,
vmin=None, vmax=None, figsize=(9, 3), cmap=None,
- cmap_set_under=None, cmap_scale='linear'):
+ cmap_set_under=None, cmap_set_over=None,
+ cmap_scale='linear'):
try:
os.makedirs(os.path.dirname(outFileName))
@@ -792,6 +797,9 @@ def _plot_horiz_field(self, field, title, outFileName, oceanDomain=True,
if cmap_set_under is not None:
current_cmap = localPatches.get_cmap()
current_cmap.set_under(cmap_set_under)
+ if cmap_set_over is not None:
+ current_cmap = localPatches.get_cmap()
+ current_cmap.set_over(cmap_set_over)
localPatches.set_edgecolor('face')
localPatches.set_clim(vmin=vmin, vmax=vmax)
diff --git a/docs/developers_guide/ocean/test_groups/isomip_plus.rst b/docs/developers_guide/ocean/test_groups/isomip_plus.rst
index 097b7f8ed0..2b2571751d 100644
--- a/docs/developers_guide/ocean/test_groups/isomip_plus.rst
+++ b/docs/developers_guide/ocean/test_groups/isomip_plus.rst
@@ -33,11 +33,6 @@ rising indefinitely due to the input of freshwater from ice-shelf melting.
geom
~~~~
-The function :py:func:`compass.ocean.tests.isomip_plus.geom.define_thin_film_mask_step1()`
-defines an MPAS mesh mask for the ocean domain including cells over the full x- and y-range
-in cases with a thin film. Thus, all land cells are included in the ocean's thin-film region.
-
-
The function :py:func:`compass.ocean.tests.isomip_plus.geom.interpolate_ocean_mask()`
interpolates the ocean mask from the BISICLES grid of the input geometry to
the MPAS-Ocean mesh. The mask can later be used to cull land cells from the
@@ -72,6 +67,19 @@ Optionally, the ice draft can be scaled by a factor as a simple way to explore
changing ice-shelf topography. Variables are renamed to those expected by
MPAS-Ocean.
+planar_mesh
+~~~~~~~~~
+
+The class :py:class:`compass.ocean.tests.isomip_plus.planar_mesh.PlanarMesh`
+defines a step for generating a planar mesh.
+
+cull_mesh
+~~~~~~~~~
+
+The class :py:class:`compass.ocean.tests.isomip_plus.cull_mesh.CullMesh`
+defines a step for culling the mesh to include only ocean cells, with or
+without a thin film depending on the test.
+
initial_state
~~~~~~~~~~~~~
diff --git a/docs/users_guide/ocean/test_groups/isomip_plus.rst b/docs/users_guide/ocean/test_groups/isomip_plus.rst
index 3582191654..567fa600b5 100644
--- a/docs/users_guide/ocean/test_groups/isomip_plus.rst
+++ b/docs/users_guide/ocean/test_groups/isomip_plus.rst
@@ -37,10 +37,12 @@ in this region.
case at 5 km resolution, showing potential temperature averaged over month
9 of the simulation.
-The ``isomip_plus`` test cases are composed of 3 steps that run by default:
-``initial_state``, which defines the mesh, interpolates the ice geometry, and
-computes the initial conditions for the model; ``ssh_adjustment``, which
-modifies the ``landIcePressure`` field to balance the ``ssh`` field, see
+The ``isomip_plus`` test cases are composed of 5 steps that run by default:
+``process_geom``, which loads the data on the ice sheet mesh; ``planar_mesh``,
+which defines the planar mesh; ``cull_mesh``, which culls the mesh;
+``initial_state``, which interpolates the ice geometry and computes the
+initial conditions for the model; ``ssh_adjustment``, which modifies the
+``landIcePressure`` field to balance the ``ssh`` field, see
:ref:`ocean_ssh_adjustment`; and ``performance``, which performs a 1-hour time
integration of the model and compares the results with a baseline if one is
provided.