From 0a0452fa6143c8a0aa3f8aa90cd3e0a22adb1a26 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 24 Oct 2023 16:55:13 -0500 Subject: [PATCH] Use new land ice pressure function in isomip_plus --- compass/ocean/tests/isomip_plus/geom.py | 3 +- .../ocean/tests/isomip_plus/initial_state.py | 37 +++++++++++++------ .../ocean/tests/isomip_plus/process_geom.py | 8 ++++ 3 files changed, 36 insertions(+), 12 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/geom.py b/compass/ocean/tests/isomip_plus/geom.py index ce761a1582..551f70b3d5 100755 --- a/compass/ocean/tests/isomip_plus/geom.py +++ b/compass/ocean/tests/isomip_plus/geom.py @@ -111,11 +111,12 @@ def interpolate_geom(ds_mesh, ds_geom, min_ocean_fraction, thin_film_present): ds_geom.landIceGroundedFraction) # mask the topography to the ocean region before interpolation - for var in ['Z_bed', 'Z_ice_draft', 'landIceFraction', + for var in ['Z_bed', 'Z_ice_surface', 'Z_ice_draft', 'landIceFraction', 'landIceFloatingFraction', 'smoothedDraftMask']: ds_geom[var] = ds_geom[var] * ds_geom['oceanFraction'] fields = {'bottomDepthObserved': 'Z_bed', + 'landIceThickness': 'iceThickness', 'ssh': 'Z_ice_draft', 'oceanFracObserved': 'oceanFraction', 'landIceFraction': 'landIceFraction', diff --git a/compass/ocean/tests/isomip_plus/initial_state.py b/compass/ocean/tests/isomip_plus/initial_state.py index 1ce870ad2e..a81da8832b 100644 --- a/compass/ocean/tests/isomip_plus/initial_state.py +++ b/compass/ocean/tests/isomip_plus/initial_state.py @@ -7,7 +7,11 @@ from mpas_tools.cime.constants import constants from mpas_tools.io import write_netcdf -from compass.ocean.iceshelf import compute_land_ice_pressure_from_draft +from compass.ocean.iceshelf import ( + compute_land_ice_density_from_draft, + compute_land_ice_pressure_from_draft, + compute_land_ice_pressure_from_thickness, +) from compass.ocean.tests.isomip_plus.geom import interpolate_geom from compass.ocean.tests.isomip_plus.viz.plot import MoviePlotter from compass.ocean.vertical import init_vertical_coord @@ -125,14 +129,20 @@ def _compute_initial_condition(self): ds['landIceFraction'] = xr.where(mask, ds.landIceFraction, 0.) - ref_density = constants['SHR_CONST_RHOSW'] - landIceDraft = ds.ssh - landIcePressure = compute_land_ice_pressure_from_draft( - land_ice_draft=landIceDraft, modify_mask=ds.ssh < 0., - ref_density=ref_density) + land_ice_draft = ds.ssh + if thin_film_present: + land_ice_thickness = ds.landIceThickness + land_ice_density = compute_land_ice_density_from_draft( + land_ice_draft, land_ice_thickness, + floating_mask=ds.landIceFloatingMask) + land_ice_pressure = compute_land_ice_pressure_from_thickness( + land_ice_thickness=land_ice_thickness, modify_mask=ds.ssh < 0., + land_ice_density=land_ice_density) + else: + land_ice_pressure = compute_land_ice_pressure_from_draft( + land_ice_draft=land_ice_draft, modify_mask=land_ice_draft < 0.) - ds['landIcePressure'] = landIcePressure - ds['landIceDraft'] = landIceDraft + ds['landIcePressure'] = land_ice_pressure if self.time_varying_forcing: self._write_time_varying_forcing(ds_init=ds) @@ -173,8 +183,8 @@ def _compute_initial_condition(self): # 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) + -7.43e-8 * land_ice_pressure + + -1.74e-10 * land_ice_pressure * 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) @@ -233,6 +243,8 @@ def _plot(self, ds): ds['oceanFracObserved'].expand_dims(dim='Time', axis=0) ds['landIcePressure'] = \ ds['landIcePressure'].expand_dims(dim='Time', axis=0) + ds['landIceThickness'] = \ + ds['landIceThickness'].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) @@ -248,7 +260,10 @@ def _plot(self, ds): True) plotter.plot_horiz_series(ds.landIcePressure, 'landIcePressure', 'landIcePressure', - True, vmin=1, vmax=1e6, cmap_scale='log') + True, vmin=1e5, vmax=1e7, cmap_scale='log') + plotter.plot_horiz_series(ds.landIceThickness, + 'landIceThickness', 'landIceThickness', + True, vmin=0, vmax=1e3) plotter.plot_horiz_series(ds.ssh, 'ssh', 'ssh', True, vmin=-700, vmax=0) diff --git a/compass/ocean/tests/isomip_plus/process_geom.py b/compass/ocean/tests/isomip_plus/process_geom.py index 1293cc78cd..3f07250b11 100644 --- a/compass/ocean/tests/isomip_plus/process_geom.py +++ b/compass/ocean/tests/isomip_plus/process_geom.py @@ -171,10 +171,18 @@ def _smooth_geometry(land_fraction, ds, filter_sigma, threshold=0.01): mode='constant', cval=0.) bed[mask] /= smoothed_mask[mask] + # We only use the ice surface where the ocean is present to compute + # land ice pressure + ice_thickness = filters.gaussian_filter( + ds.iceThickness * ocean_fraction, filter_sigma, mode='constant', + cval=0.) + ice_thickness[mask] /= smoothed_mask[mask] + 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) ds['Z_bed'] = (('y', 'x'), bed) + ds['iceThickness'] = (('y', 'x'), ice_thickness) ds['smoothedDraftMask'] = (('y', 'x'), smoothed_draft_mask)