From ef1df8fd336263ecb42c318cf4f3eb9d02e43a40 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 24 Oct 2023 16:20:35 -0500 Subject: [PATCH 01/30] Reorg land ice pressure function --- compass/ocean/iceshelf.py | 26 +++++++++---------- .../ocean/tests/ice_shelf_2d/initial_state.py | 8 +++--- .../ocean/tests/isomip_plus/initial_state.py | 8 +++--- 3 files changed, 23 insertions(+), 19 deletions(-) diff --git a/compass/ocean/iceshelf.py b/compass/ocean/iceshelf.py index 8a69efa617..6c74c26dfb 100644 --- a/compass/ocean/iceshelf.py +++ b/compass/ocean/iceshelf.py @@ -11,34 +11,34 @@ from compass.model import partition, run_model -def compute_land_ice_pressure_and_draft(ssh, modify_mask, ref_density): +def compute_land_ice_pressure_from_draft(land_ice_draft, modify_mask, + ref_density=None): """ - Compute the pressure from and overlying ice shelf and the ice-shelf draft + Compute the pressure from and overlying ice shelf Parameters ---------- - ssh : xarray.DataArray - The sea surface height (the ice draft) + land_ice_draft : xarray.DataArray + The ice draft (sea surface height) modify_mask : xarray.DataArray A mask that is 1 where ``landIcePressure`` can be deviate from 0 - ref_density : float + ref_density : float, optional A reference density for seawater displaced by the ice shelf Returns ------- - landIcePressure : xarray.DataArray + land_ice_pressure : xarray.DataArray The pressure from the overlying land ice on the ocean - - landIceDraft : xarray.DataArray - The ice draft, equal to the initial ``ssh`` """ gravity = constants['SHR_CONST_G'] - landIcePressure = \ - modify_mask * numpy.maximum(-ref_density * gravity * ssh, 0.) - landIceDraft = ssh - return landIcePressure, landIceDraft + if ref_density is None: + ref_density = constants['SHR_CONST_RHOSW'] + land_ice_pressure = \ + modify_mask * numpy.maximum(-ref_density * gravity * land_ice_draft, + 0.) + return land_ice_pressure def adjust_ssh(variable, iteration_count, step, update_pio=True, diff --git a/compass/ocean/tests/ice_shelf_2d/initial_state.py b/compass/ocean/tests/ice_shelf_2d/initial_state.py index d245f3c9c8..cb76da4fbf 100644 --- a/compass/ocean/tests/ice_shelf_2d/initial_state.py +++ b/compass/ocean/tests/ice_shelf_2d/initial_state.py @@ -5,7 +5,7 @@ from mpas_tools.mesh.conversion import convert, cull from mpas_tools.planar_hex import make_planar_hex_mesh -from compass.ocean.iceshelf import compute_land_ice_pressure_and_draft +from compass.ocean.iceshelf import compute_land_ice_pressure_from_draft from compass.ocean.vertical import init_vertical_coord from compass.step import Step @@ -103,8 +103,10 @@ def run(self): landIceFloatingMask = landIceMask.copy() ref_density = constants['SHR_CONST_RHOSW'] - landIcePressure, landIceDraft = compute_land_ice_pressure_and_draft( - ssh=ds.ssh, modify_mask=modify_mask, ref_density=ref_density) + landIceDraft = ds.ssh + landIcePressure = compute_land_ice_pressure_from_draft( + land_ice_draft=landIceDraft, modify_mask=modify_mask, + ref_density=ref_density) salinity = surface_salinity + ((bottom_salinity - surface_salinity) * (ds.zMid / (-bottom_depth))) diff --git a/compass/ocean/tests/isomip_plus/initial_state.py b/compass/ocean/tests/isomip_plus/initial_state.py index 955651040b..4400f6c574 100644 --- a/compass/ocean/tests/isomip_plus/initial_state.py +++ b/compass/ocean/tests/isomip_plus/initial_state.py @@ -7,7 +7,7 @@ from mpas_tools.cime.constants import constants from mpas_tools.io import write_netcdf -from compass.ocean.iceshelf import compute_land_ice_pressure_and_draft +from compass.ocean.iceshelf import compute_land_ice_pressure_from_draft 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 @@ -126,8 +126,10 @@ def _compute_initial_condition(self): ds['landIceFraction'] = xr.where(mask, ds.landIceFraction, 0.) ref_density = constants['SHR_CONST_RHOSW'] - landIcePressure, landIceDraft = compute_land_ice_pressure_and_draft( - ssh=ds.ssh, modify_mask=ds.ssh < 0., ref_density=ref_density) + landIceDraft = ds.ssh + landIcePressure = compute_land_ice_pressure_from_draft( + land_ice_draft=landIceDraft, modify_mask=ds.ssh < 0., + ref_density=ref_density) ds['landIcePressure'] = landIcePressure ds['landIceDraft'] = landIceDraft From 85cf30e620011efe06e60f0e09980eab15e6ebc1 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 24 Oct 2023 16:24:53 -0500 Subject: [PATCH 02/30] Add new land ice pressure function --- compass/ocean/iceshelf.py | 93 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) diff --git a/compass/ocean/iceshelf.py b/compass/ocean/iceshelf.py index 6c74c26dfb..68b2680cf8 100644 --- a/compass/ocean/iceshelf.py +++ b/compass/ocean/iceshelf.py @@ -11,6 +11,69 @@ from compass.model import partition, run_model +def compute_land_ice_pressure_from_thickness(land_ice_thickness, modify_mask, + land_ice_density=None): + """ + Compute the pressure from and overlying ice shelf + + Parameters + ---------- + land_ice_thickness: xarray.DataArray + The ice thickness + + modify_mask : xarray.DataArray + A mask that is 1 where ``landIcePressure`` can be deviate from 0 + + land_ice_density : float, optional + A reference density for land ice + + Returns + ------- + land_ice_pressure : xarray.DataArray + The pressure from the overlying land ice on the ocean + """ + gravity = constants['SHR_CONST_G'] + if land_ice_density is None: + land_ice_density = constants['SHR_CONST_RHOICE'] + land_ice_pressure = modify_mask * \ + numpy.maximum(land_ice_density * gravity * land_ice_thickness, 0.) + return land_ice_pressure + + +def compute_land_ice_density_from_draft(land_ice_draft, land_ice_thickness, + floating_mask, ref_density=None): + """ + Compute the spatially-averaged ice density needed to match the ice draft + + Parameters + ---------- + land_ice_draft : xarray.DataArray + The ice draft (sea surface height) + + land_ice_thickness: xarray.DataArray + The ice thickness + + floating_mask : xarray.DataArray + A mask that is 1 where the ice is assumed in hydrostatic equilibrium + + ref_density : float, optional + A reference density for seawater displaced by the ice shelf + + Returns + ------- + land_ice_density: float + The ice density + """ + if ref_density is None: + ref_density = constants['SHR_CONST_RHOSW'] + land_ice_draft = numpy.where(floating_mask, land_ice_draft, numpy.nan) + land_ice_thickness = numpy.where(floating_mask, land_ice_thickness, + numpy.nan) + land_ice_density = \ + numpy.nanmean(-ref_density * land_ice_draft / land_ice_thickness) + return land_ice_density + + def compute_land_ice_pressure_from_draft(land_ice_draft, modify_mask, ref_density=None): """ @@ -41,6 +104,36 @@ def compute_land_ice_pressure_from_draft(land_ice_draft, modify_mask, return land_ice_pressure +def compute_land_ice_draft_from_pressure(land_ice_pressure, modify_mask, + ref_density=None): + """ + Compute the ice-shelf draft associated with the pressure from an overlying + ice shelf + + Parameters + ---------- + land_ice_pressure : xarray.DataArray + The pressure from the overlying land ice on the ocean + + modify_mask : xarray.DataArray + A mask that is 1 where ``landIcePressure`` can be deviate from 0 + + ref_density : float, optional + A reference density for seawater displaced by the ice shelf + + Returns + ------- + land_ice_draft : xarray.DataArray + The ice draft + """ + gravity = constants['SHR_CONST_G'] + if ref_density is None: + ref_density = constants['SHR_CONST_RHOSW'] + land_ice_draft = \ + - (modify_mask * land_ice_pressure / (ref_density * gravity)) + return land_ice_draft + + def adjust_ssh(variable, iteration_count, step, update_pio=True, convert_to_cdf5=False, delta_ssh_threshold=None): """ From fa78075829963ca680f8fde55d4dbcdfbc6da7f5 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 24 Oct 2023 16:55:13 -0500 Subject: [PATCH 03/30] 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 4400f6c574..516147b96c 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) From 7b05efbefdffb031ef003b256d56d747a1b956c3 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 25 Oct 2023 12:56:28 -0500 Subject: [PATCH 04/30] fixup thin film culling --- compass/ocean/tests/isomip_plus/geom.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/ocean/tests/isomip_plus/geom.py b/compass/ocean/tests/isomip_plus/geom.py index 551f70b3d5..3c52630511 100755 --- a/compass/ocean/tests/isomip_plus/geom.py +++ b/compass/ocean/tests/isomip_plus/geom.py @@ -157,7 +157,7 @@ def _get_geom_fields(ds_geom, ds_mesh, thin_film_present): y_cell = ds_mesh.yIsomipCell.values if thin_film_present: - ocean_fraction = - ds_geom['landFraction'] + 1.0 + ocean_fraction = xarray.where(ds_geom['Z_bed'] > 0., 0., 1.) else: ocean_fraction = (ds_geom['landIceFloatingFraction'] + ds_geom['openOceanFraction']) From 2621e00707cd2b1c991f9e1c97abfcd28ef663dd Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 26 Oct 2023 11:18:19 -0500 Subject: [PATCH 05/30] Remove ice density calculation --- compass/ocean/iceshelf.py | 34 ------------------- .../ocean/tests/isomip_plus/initial_state.py | 18 +++++----- .../ocean/tests/isomip_plus/isomip_plus.cfg | 3 ++ 3 files changed, 13 insertions(+), 42 deletions(-) diff --git a/compass/ocean/iceshelf.py b/compass/ocean/iceshelf.py index 68b2680cf8..08bf67e846 100644 --- a/compass/ocean/iceshelf.py +++ b/compass/ocean/iceshelf.py @@ -40,40 +40,6 @@ def compute_land_ice_pressure_from_thickness(land_ice_thickness, modify_mask, return land_ice_pressure -def compute_land_ice_density_from_draft(land_ice_draft, land_ice_thickness, - floating_mask, ref_density=None): - """ - Compute the spatially-averaged ice density needed to match the ice draft - - Parameters - ---------- - land_ice_draft : xarray.DataArray - The ice draft (sea surface height) - - land_ice_thickness: xarray.DataArray - The ice thickness - - floating_mask : xarray.DataArray - A mask that is 1 where the ice is assumed in hydrostatic equilibrium - - ref_density : float, optional - A reference density for seawater displaced by the ice shelf - - Returns - ------- - land_ice_density: float - The ice density - """ - if ref_density is None: - ref_density = constants['SHR_CONST_RHOSW'] - land_ice_draft = numpy.where(floating_mask, land_ice_draft, numpy.nan) - land_ice_thickness = numpy.where(floating_mask, land_ice_thickness, - numpy.nan) - land_ice_density = \ - numpy.nanmean(-ref_density * land_ice_draft / land_ice_thickness) - return land_ice_density - - def compute_land_ice_pressure_from_draft(land_ice_draft, modify_mask, ref_density=None): """ diff --git a/compass/ocean/tests/isomip_plus/initial_state.py b/compass/ocean/tests/isomip_plus/initial_state.py index 516147b96c..7124198c29 100644 --- a/compass/ocean/tests/isomip_plus/initial_state.py +++ b/compass/ocean/tests/isomip_plus/initial_state.py @@ -8,7 +8,7 @@ from mpas_tools.io import write_netcdf from compass.ocean.iceshelf import ( - compute_land_ice_density_from_draft, + compute_land_ice_draft_from_pressure, compute_land_ice_pressure_from_draft, compute_land_ice_pressure_from_thickness, ) @@ -101,6 +101,7 @@ def _compute_initial_condition(self): section = config['isomip_plus'] min_land_ice_fraction = section.getfloat('min_land_ice_fraction') min_ocean_fraction = section.getfloat('min_ocean_fraction') + ice_density = section.getfloat('ice_density') ds_geom = xr.open_dataset('input_geometry_processed.nc') ds_mesh = xr.open_dataset('culled_mesh.nc') @@ -108,6 +109,8 @@ def _compute_initial_condition(self): ds = interpolate_geom(ds_mesh, ds_geom, min_ocean_fraction, thin_film_present) + ds['bottomDepth'] = -ds.bottomDepthObserved + ds['landIceFraction'] = \ ds['landIceFraction'].expand_dims(dim='Time', axis=0) ds['landIceFloatingFraction'] = \ @@ -129,16 +132,17 @@ def _compute_initial_condition(self): ds['landIceFraction'] = xr.where(mask, ds.landIceFraction, 0.) - 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) + land_ice_density=ice_density) + land_ice_draft = compute_land_ice_draft_from_pressure( + land_ice_pressure=land_ice_pressure, + modify_mask=ds.bottomDepth > 0.) + ds['ssh'] = np.maximum(land_ice_draft, -ds.bottomDepth) else: + land_ice_draft = ds.ssh land_ice_pressure = compute_land_ice_pressure_from_draft( land_ice_draft=land_ice_draft, modify_mask=land_ice_draft < 0.) @@ -147,8 +151,6 @@ def _compute_initial_condition(self): if self.time_varying_forcing: self._write_time_varying_forcing(ds_init=ds) - ds['bottomDepth'] = -ds.bottomDepthObserved - section = config['isomip_plus'] # Deepen the bottom depth to maintain the minimum water-column diff --git a/compass/ocean/tests/isomip_plus/isomip_plus.cfg b/compass/ocean/tests/isomip_plus/isomip_plus.cfg index 355daa71f8..65caf293aa 100644 --- a/compass/ocean/tests/isomip_plus/isomip_plus.cfg +++ b/compass/ocean/tests/isomip_plus/isomip_plus.cfg @@ -78,6 +78,9 @@ min_smoothed_draft_mask = 0.01 # considered a land-ice cell by MPAS-Ocean (landIceMask == 1). min_land_ice_fraction = 0.5 +# the density of ice prescribed in ISOMIP+ +ice_density = 918 + # the initial temperature at the sea surface init_top_temp = -1.9 # the initial temperature at the sea floor From 5a542263686a03fd668b84527d879fef5870aebb Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 26 Oct 2023 15:05:33 -0500 Subject: [PATCH 06/30] Extend treatment to time-varying forcing --- .../ocean/tests/isomip_plus/initial_state.py | 43 ++++++++++++++----- 1 file changed, 32 insertions(+), 11 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/initial_state.py b/compass/ocean/tests/isomip_plus/initial_state.py index 7124198c29..96b5181349 100644 --- a/compass/ocean/tests/isomip_plus/initial_state.py +++ b/compass/ocean/tests/isomip_plus/initial_state.py @@ -148,9 +148,6 @@ def _compute_initial_condition(self): ds['landIcePressure'] = land_ice_pressure - if self.time_varying_forcing: - self._write_time_varying_forcing(ds_init=ds) - section = config['isomip_plus'] # Deepen the bottom depth to maintain the minimum water-column @@ -215,6 +212,15 @@ def _compute_initial_condition(self): ds['normalVelocity'] = normalVelocity.expand_dims(dim='Time', axis=0) write_netcdf(ds, 'initial_state.nc') + + if self.time_varying_forcing: + self._write_time_varying_forcing(ds_init=ds, + ice_density=ice_density) + else: + # We need to add the time dimension for plotting purposes + ds['landIcePressure'] = \ + ds['landIcePressure'].expand_dims(dim='Time', axis=0) + return ds, frac def _plot(self, ds): @@ -243,8 +249,6 @@ def _plot(self, ds): ds['oceanFracObserved'] = \ 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'] = \ @@ -382,7 +386,7 @@ def _compute_restoring(self, ds, frac): write_netcdf(ds_forcing, 'init_mode_forcing_data.nc') - def _write_time_varying_forcing(self, ds_init): + def _write_time_varying_forcing(self, ds_init, ice_density): """ Write time-varying land-ice forcing and update the initial condition """ @@ -402,10 +406,27 @@ def _write_time_varying_forcing(self, ds_init): landIceFraction = list() landIceFloatingFraction = list() + if self.thin_film_present: + land_ice_thickness = ds_init.landIceThickness + land_ice_pressure = compute_land_ice_pressure_from_thickness( + land_ice_thickness=land_ice_thickness, + modify_mask=ds_init.landIceFraction > 0., + land_ice_density=ice_density) + land_ice_draft = compute_land_ice_draft_from_pressure( + land_ice_pressure=land_ice_pressure, + modify_mask=ds_init.bottomDepth > 0.) + land_ice_draft = np.maximum(land_ice_draft, -ds_init.bottomDepth) + land_ice_draft = land_ice_draft.transpose() + else: + land_ice_draft = ds_init.landIceDraft + land_ice_pressure = ds_init.landIcePressure + for scale in scales: - landIceDraft.append(scale * ds_init.landIceDraft) - landIcePressure.append(scale * ds_init.landIcePressure) + landIceDraft.append(scale * land_ice_draft) + landIcePressure.append(scale * land_ice_pressure) landIceFraction.append(ds_init.landIceFraction) + # Since floating fraction does not change, none of the thin film + # cases allow for the area undergoing melting to change landIceFloatingFraction.append(ds_init.landIceFloatingFraction) ds_out['landIceDraftForcing'] = xr.concat(landIceDraft, 'Time') @@ -427,6 +448,6 @@ def _write_time_varying_forcing(self, ds_init): '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 - ds_init['ssh'] = ds_init.landIceDraft - ds_init['landIcePressure'] = scales[0] * ds_init.landIcePressure + ds_init['landIceDraft'] = scales[0] * land_ice_draft + ds_init['ssh'] = land_ice_draft + ds_init['landIcePressure'] = scales[0] * land_ice_pressure From 9325263ff35c5a98cb43bb272dae77c61facfcaa Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 26 Oct 2023 18:34:46 -0500 Subject: [PATCH 07/30] Update docs --- compass/ocean/iceshelf.py | 4 ++-- docs/developers_guide/ocean/api.rst | 4 +++- docs/users_guide/ocean/test_groups/isomip_plus.rst | 3 +++ 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/compass/ocean/iceshelf.py b/compass/ocean/iceshelf.py index 08bf67e846..d15d0204e7 100644 --- a/compass/ocean/iceshelf.py +++ b/compass/ocean/iceshelf.py @@ -14,7 +14,7 @@ def compute_land_ice_pressure_from_thickness(land_ice_thickness, modify_mask, land_ice_density=None): """ - Compute the pressure from and overlying ice shelf + Compute the pressure from an overlying ice shelf from ice thickness Parameters ---------- @@ -43,7 +43,7 @@ def compute_land_ice_pressure_from_thickness(land_ice_thickness, modify_mask, def compute_land_ice_pressure_from_draft(land_ice_draft, modify_mask, ref_density=None): """ - Compute the pressure from and overlying ice shelf + Compute the pressure from an overlying ice shelf from ice draft Parameters ---------- diff --git a/docs/developers_guide/ocean/api.rst b/docs/developers_guide/ocean/api.rst index a2e1845d54..a9f50a81cb 100644 --- a/docs/developers_guide/ocean/api.rst +++ b/docs/developers_guide/ocean/api.rst @@ -1042,7 +1042,9 @@ ocean framework haney.compute_haney_number - iceshelf.compute_land_ice_pressure_and_draft + iceshelf.compute_land_ice_draft_from_pressure + iceshelf.compute_land_ice_pressure_from_draft + iceshelf.compute_land_ice_pressure_from_thickness iceshelf.adjust_ssh particles.write diff --git a/docs/users_guide/ocean/test_groups/isomip_plus.rst b/docs/users_guide/ocean/test_groups/isomip_plus.rst index 567fa600b5..e53a4940b5 100644 --- a/docs/users_guide/ocean/test_groups/isomip_plus.rst +++ b/docs/users_guide/ocean/test_groups/isomip_plus.rst @@ -125,6 +125,9 @@ The ``isomip_plus`` test cases share the following config options: # considered a land-ice cell by MPAS-Ocean (landIceMask == 1). min_land_ice_fraction = 0.5 + # the density of ice prescribed in ISOMIP+ + ice_density = 918 + # the initial temperature at the sea surface init_top_temp = -1.9 # the initial temperature at the sea floor From 99577697f296d13cc478b03b35151668d65d43d9 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 27 Nov 2023 13:05:34 -0600 Subject: [PATCH 08/30] Add landIceDraft to initial_state --- compass/ocean/tests/isomip_plus/initial_state.py | 1 + 1 file changed, 1 insertion(+) diff --git a/compass/ocean/tests/isomip_plus/initial_state.py b/compass/ocean/tests/isomip_plus/initial_state.py index 96b5181349..52183a7154 100644 --- a/compass/ocean/tests/isomip_plus/initial_state.py +++ b/compass/ocean/tests/isomip_plus/initial_state.py @@ -147,6 +147,7 @@ def _compute_initial_condition(self): land_ice_draft=land_ice_draft, modify_mask=land_ice_draft < 0.) ds['landIcePressure'] = land_ice_pressure + ds['landIceDraft'] = land_ice_draft section = config['isomip_plus'] From 46ecac6eadede9d55c7b6d4146c07f2bf58a25bc Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Tue, 30 Jan 2024 17:45:01 -0600 Subject: [PATCH 09/30] Apply suggestions from code review --- .../ocean/tests/isomip_plus/initial_state.py | 36 ++++++++++--------- 1 file changed, 19 insertions(+), 17 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/initial_state.py b/compass/ocean/tests/isomip_plus/initial_state.py index 52183a7154..19cd684827 100644 --- a/compass/ocean/tests/isomip_plus/initial_state.py +++ b/compass/ocean/tests/isomip_plus/initial_state.py @@ -248,16 +248,19 @@ def _plot(self, ds): sectionY=section_y, dsMesh=ds, ds=ds, showProgress=show_progress) - ds['oceanFracObserved'] = \ + ssh = ds['ssh'].expand_dims(dim='Time', axis=0) + oceanFracObserved = \ ds['oceanFracObserved'].expand_dims(dim='Time', axis=0) - ds['landIceThickness'] = \ + oceanFracObserved = \ + ds['oceanFracObserved'].expand_dims(dim='Time', axis=0) + landIcePressure = \ + ds['landIcePressure'].expand_dims(dim='Time', axis=0) + landIceThickness = \ ds['landIceThickness'].expand_dims(dim='Time', axis=0) - ds['landIceGroundedFraction'] = \ + 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') + bottomDepth = ds['bottomDepth'].expand_dims(dim='Time', axis=0) + totalColThickness = ds.layerThickness.sum(dim='nVertLevels') tol = 1e-10 plotter.plot_horiz_series(ds.landIceMask, 'landIceMask', 'landIceMask', @@ -265,23 +268,22 @@ def _plot(self, ds): plotter.plot_horiz_series(ds.landIceFloatingMask, 'landIceFloatingMask', 'landIceFloatingMask', True) - plotter.plot_horiz_series(ds.landIcePressure, + plotter.plot_horiz_series(landIcePressure, 'landIcePressure', 'landIcePressure', True, vmin=1e5, vmax=1e7, cmap_scale='log') - plotter.plot_horiz_series(ds.landIceThickness, + plotter.plot_horiz_series(landIceThickness, 'landIceThickness', 'landIceThickness', True, vmin=0, vmax=1e3) - plotter.plot_horiz_series(ds.ssh, - 'ssh', 'ssh', + plotter.plot_horiz_series(ssh, 'ssh', 'ssh', True, vmin=-700, vmax=0) - plotter.plot_horiz_series(ds.bottomDepth, + plotter.plot_horiz_series(bottomDepth, 'bottomDepth', 'bottomDepth', True, vmin=0, vmax=700) - plotter.plot_horiz_series(ds.ssh + ds.bottomDepth, + plotter.plot_horiz_series(ds.ssh + bottomDepth, 'H', 'H', True, vmin=min_column_thickness + tol, vmax=700, cmap_set_under='r', cmap_scale='log') - plotter.plot_horiz_series(ds.totalColThickness, + plotter.plot_horiz_series(totalColThickness, 'totalColThickness', 'totalColThickness', True, vmin=min_column_thickness + 1e-10, vmax=700, cmap_set_under='r') @@ -296,13 +298,13 @@ def _plot(self, ds): True, vmin=0 + tol, vmax=1 - tol, cmap='cmo.balance', cmap_set_under='k', cmap_set_over='r') - plotter.plot_horiz_series(ds.landIceGroundedFraction, + plotter.plot_horiz_series(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, + plotter.plot_horiz_series(oceanFracObserved, 'oceanFracObserved', 'oceanFracObserved', True, vmin=0 + tol, vmax=1 - tol, cmap='cmo.balance', @@ -417,7 +419,7 @@ def _write_time_varying_forcing(self, ds_init, ice_density): land_ice_pressure=land_ice_pressure, modify_mask=ds_init.bottomDepth > 0.) land_ice_draft = np.maximum(land_ice_draft, -ds_init.bottomDepth) - land_ice_draft = land_ice_draft.transpose() + land_ice_draft = land_ice_draft.transpose('nCells', 'nVertLevels') else: land_ice_draft = ds_init.landIceDraft land_ice_pressure = ds_init.landIcePressure From f648a24f296a10dd04f2289d8e70236e68f5213a Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Wed, 23 Oct 2024 16:57:57 -0500 Subject: [PATCH 10/30] Fixup time varying forcing --- .../ocean/tests/isomip_plus/initial_state.py | 293 ++++++++++-------- 1 file changed, 167 insertions(+), 126 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/initial_state.py b/compass/ocean/tests/isomip_plus/initial_state.py index 19cd684827..a525cc3cc5 100644 --- a/compass/ocean/tests/isomip_plus/initial_state.py +++ b/compass/ocean/tests/isomip_plus/initial_state.py @@ -132,27 +132,53 @@ def _compute_initial_condition(self): ds['landIceFraction'] = xr.where(mask, ds.landIceFraction, 0.) - if thin_film_present: + if thin_film_present: # compute full pressure in grounded regions land_ice_thickness = ds.landIceThickness - land_ice_pressure = compute_land_ice_pressure_from_thickness( - land_ice_thickness=land_ice_thickness, modify_mask=ds.ssh < 0., - land_ice_density=ice_density) - land_ice_draft = compute_land_ice_draft_from_pressure( - land_ice_pressure=land_ice_pressure, - modify_mask=ds.bottomDepth > 0.) - ds['ssh'] = np.maximum(land_ice_draft, -ds.bottomDepth) - else: + land_ice_pressure_unscaled = \ + compute_land_ice_pressure_from_thickness( + land_ice_thickness=land_ice_thickness, + modify_mask=ds.ssh < 0., + land_ice_density=ice_density) + else: # assume floatation pressure everywhere land_ice_draft = ds.ssh - land_ice_pressure = compute_land_ice_pressure_from_draft( + land_ice_pressure_unscaled = compute_land_ice_pressure_from_draft( land_ice_draft=land_ice_draft, modify_mask=land_ice_draft < 0.) - ds['landIcePressure'] = land_ice_pressure + if self.time_varying_forcing: + scales = config.get('isomip_plus_forcing', 'scales') + scales = [float(scale) + for scale in scales.replace(',', ' ').split()] + land_ice_pressure = land_ice_pressure_unscaled * scales[0] + else: + land_ice_pressure = land_ice_pressure_unscaled + + if thin_film_present: + modify_mask = ds.bottomDepth > 0. + land_ice_draft = compute_land_ice_draft_from_pressure( + land_ice_pressure=land_ice_pressure, + modify_mask=modify_mask) + land_ice_draft = np.maximum(land_ice_draft, + -ds.bottomDepth) + ds['ssh'] = land_ice_draft + + # Previously, we didn't apply max draft at seafloor but we want this + # for W&D + # This assumes that for MALI forcing, landIceDraft = bottomDepth in + # grounded regions ds['landIceDraft'] = land_ice_draft + # We need to add the time dimension for plotting purposes + ds['landIcePressure'] = land_ice_pressure.expand_dims( + dim='Time', axis=0) - section = config['isomip_plus'] + if self.time_varying_forcing: + self._write_time_varying_forcing( + ds_init=ds, + ice_density=ice_density, + land_ice_pressure_unscaled=land_ice_pressure_unscaled) # Deepen the bottom depth to maintain the minimum water-column # thickness + section = config['isomip_plus'] min_column_thickness = section.getfloat('min_column_thickness') min_layer_thickness = section.getfloat('min_layer_thickness') min_levels = section.getint('minimum_levels') @@ -161,44 +187,52 @@ def _compute_initial_condition(self): min_depth = -ds.ssh + min_column_thickness ds['bottomDepth'] = np.maximum(ds.bottomDepth, min_depth) print(f'Adjusted bottomDepth for ' - f'{np.sum(ds.bottomDepth.values (x_max - 0.6 * self.resolution * 1e3), 1.0, - 0.0) - else: - ds_forcing['tidalInputMask'] = xr.zeros_like(frac) + x_max = np.max(ds.xIsomipCell.values) + ds_forcing['tidalInputMask'] = xr.where( + ds.xIsomipCell > (x_max - 0.6 * self.resolution * 1e3), 1.0, + 0.0) write_netcdf(ds_forcing, 'init_mode_forcing_data.nc') - def _write_time_varying_forcing(self, ds_init, ice_density): + def _write_time_varying_forcing(self, ds_init, ice_density, + land_ice_pressure_unscaled): """ - Write time-varying land-ice forcing and update the initial condition + Write time-varying land-ice forcing """ config = self.config + ds_forcing = xr.Dataset() dates = config.get('isomip_plus_forcing', 'dates') - dates = [date.ljust(64) for date in dates.replace(',', ' ').split()] + dates = [date.ljust(64) + for date in dates.replace(',', ' ').split()] scales = config.get('isomip_plus_forcing', 'scales') - scales = [float(scale) for scale in scales.replace(',', ' ').split()] - - ds_out = xr.Dataset() - ds_out['xtime'] = ('Time', dates) - ds_out['xtime'] = ds_out.xtime.astype('S') - - landIceDraft = list() - landIcePressure = list() - landIceFraction = list() - landIceFloatingFraction = list() - - if self.thin_film_present: - land_ice_thickness = ds_init.landIceThickness - land_ice_pressure = compute_land_ice_pressure_from_thickness( - land_ice_thickness=land_ice_thickness, - modify_mask=ds_init.landIceFraction > 0., - land_ice_density=ice_density) - land_ice_draft = compute_land_ice_draft_from_pressure( - land_ice_pressure=land_ice_pressure, - modify_mask=ds_init.bottomDepth > 0.) - land_ice_draft = np.maximum(land_ice_draft, -ds_init.bottomDepth) - land_ice_draft = land_ice_draft.transpose('nCells', 'nVertLevels') - else: - land_ice_draft = ds_init.landIceDraft - land_ice_pressure = ds_init.landIcePressure - - for scale in scales: - landIceDraft.append(scale * land_ice_draft) - landIcePressure.append(scale * land_ice_pressure) - landIceFraction.append(ds_init.landIceFraction) + scales = [float(scale) + for scale in scales.replace(',', ' ').split()] + + # The initial condition already has the first scale value applied + land_ice_pressure_forcing = ds_init.landIcePressure.copy() + land_ice_fraction_forcing = ds_init.landIceFraction.copy() + land_ice_floating_fraction_forcing = \ + ds_init.landIceFloatingFraction.copy() + land_ice_draft_unscaled = ds_init.landIceDraft.copy() + land_ice_draft_scaled = land_ice_draft_unscaled * scales[0] + land_ice_draft_limited = np.maximum(land_ice_draft_scaled, + -ds_init.bottomDepth) + land_ice_draft_forcing = land_ice_draft_limited + print(f'Grounded cells at {scales[0]}: ' + f'{np.sum(land_ice_draft_limited == -ds_init.bottomDepth)}') + + # We add additional time slices for the remaining scale values + for scale in scales[1:]: + land_ice_pressure_forcing = xr.concat( + [land_ice_pressure_forcing, + scale * land_ice_pressure_unscaled], + 'Time') + land_ice_fraction_forcing = xr.concat( + [land_ice_fraction_forcing, ds_init.landIceFraction], + 'Time') # Since floating fraction does not change, none of the thin film # cases allow for the area undergoing melting to change - landIceFloatingFraction.append(ds_init.landIceFloatingFraction) - - ds_out['landIceDraftForcing'] = xr.concat(landIceDraft, 'Time') - ds_out.landIceDraftForcing.attrs['units'] = 'm' - ds_out.landIceDraftForcing.attrs['long_name'] = \ + land_ice_floating_fraction_forcing = xr.concat( + [land_ice_floating_fraction_forcing, + ds_init.landIceFloatingFraction], + 'Time') + + if self.thin_film_present: + land_ice_draft_scaled = compute_land_ice_draft_from_pressure( + land_ice_pressure=scale * land_ice_pressure_unscaled, + modify_mask=ds_init.bottomDepth > 0.) + else: + # Just scale draft in the same manner as pressure + land_ice_draft_scaled = land_ice_draft_unscaled * scale + land_ice_draft_limited = np.maximum(land_ice_draft_scaled, + -ds_init.bottomDepth) + print(f'Grounded cells at {scale}: ' + f'{np.sum(land_ice_draft_limited == -ds_init.bottomDepth)}') + # Set the maximum land ice draft in grounded regions + land_ice_draft_forcing = xr.concat( + [land_ice_draft_forcing, + land_ice_draft_limited], + 'Time') + + ds_forcing['xtime'] = xr.DataArray(data=dates, + dims=('Time')).astype('S') + ds_forcing['landIceDraftForcing'] = land_ice_draft_forcing + ds_forcing.landIceDraftForcing.attrs['units'] = 'm' + ds_forcing.landIceDraftForcing.attrs['long_name'] = \ 'The approximate elevation of the land ice-ocean interface' - ds_out['landIcePressureForcing'] = \ - xr.concat(landIcePressure, 'Time') - ds_out.landIcePressureForcing.attrs['units'] = 'm' - ds_out.landIcePressureForcing.attrs['long_name'] = \ - 'Pressure from the weight of land ice at the ice-ocean interface' - ds_out['landIceFractionForcing'] = \ - xr.concat(landIceFraction, 'Time') - ds_out.landIceFractionForcing.attrs['long_name'] = \ + ds_forcing['landIcePressureForcing'] = land_ice_pressure_forcing + ds_forcing.landIcePressureForcing.attrs['units'] = 'm' + ds_forcing.landIcePressureForcing.attrs['long_name'] = \ + 'Pressure from the weight of land ice at the ice-ocean ' \ + 'interface' + ds_forcing['landIceFractionForcing'] = land_ice_fraction_forcing + ds_forcing.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'] = \ + ds_forcing['landIceFloatingFractionForcing'] = \ + land_ice_floating_fraction_forcing + ds_forcing.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] * land_ice_draft - ds_init['ssh'] = land_ice_draft - ds_init['landIcePressure'] = scales[0] * land_ice_pressure + ds_forcing.encoding['unlimited_dims'] = {'Time'} + # write_netcdf is not used here because it does not yet support + # multiple time levels + ds_forcing.to_netcdf('land_ice_forcing.nc') From 01dde0da4d61d843e11874eadf4096d7e84de5a5 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 8 Jul 2024 12:37:58 -0500 Subject: [PATCH 11/30] Make compute_land_ice_draft_from_pressure more flexible --- compass/ocean/iceshelf.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/compass/ocean/iceshelf.py b/compass/ocean/iceshelf.py index d15d0204e7..137361f1b1 100644 --- a/compass/ocean/iceshelf.py +++ b/compass/ocean/iceshelf.py @@ -95,8 +95,11 @@ def compute_land_ice_draft_from_pressure(land_ice_pressure, modify_mask, gravity = constants['SHR_CONST_G'] if ref_density is None: ref_density = constants['SHR_CONST_RHOSW'] - land_ice_draft = \ - - (modify_mask * land_ice_pressure / (ref_density * gravity)) + land_ice_draft_array = \ + - (modify_mask.values * + land_ice_pressure.values / (ref_density * gravity)) + land_ice_draft = xarray.DataArray(data=land_ice_draft_array, + dims=(land_ice_pressure.dims)) return land_ice_draft From e8b58a06daa4ed8410804b7cbb428af0c70f3eb6 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 8 Jul 2024 12:39:40 -0500 Subject: [PATCH 12/30] Plot pressure adjustment in grounded ice region --- compass/ocean/tests/isomip_plus/viz/__init__.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/viz/__init__.py b/compass/ocean/tests/isomip_plus/viz/__init__.py index 5b9d745e3a..88be3341e9 100644 --- a/compass/ocean/tests/isomip_plus/viz/__init__.py +++ b/compass/ocean/tests/isomip_plus/viz/__init__.py @@ -99,6 +99,13 @@ def run(self): True, vmin=min_column_thickness + 1e-10, vmax=700, cmap_set_under='r', cmap_scale='log') + plotter.plot_horiz_series(dsOut.surfacePressure, + 'surfacePressure', 'surfacePressure', + True, vmin=1e5, vmax=1e7, cmap_scale='log') + delSurfacePressure = dsOut.landIcePressure - dsOut.surfacePressure + plotter.plot_horiz_series(delSurfacePressure, + 'delSurfacePressure', 'delSurfacePressure', + True, vmin=1e5, vmax=1e7, cmap_scale='log') if 'tidal' in expt: delssh = dsOut.ssh - dsOut.ssh[0, :] @@ -111,10 +118,10 @@ def run(self): wct_thin = wct[:, idx_thin] wct_mean = wct_thin.mean(dim='nCells').values time = dsOut.daysSinceStartOfSim.values - fig = plt.figure() + plt.figure() plt.plot(time, wct_mean, '.') - fig.set_xlabel('Time (days)') - fig.set_ylabel('Mean thickness of thin film (m)') + plt.xlabel('Time (days)') + plt.ylabel('Mean thickness of thin film (m)') plt.savefig('wct_thin_t.png') plt.close() From d2bf5a2337e0dc177316acf58994bd47bdb17e66 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 8 Jul 2024 12:40:15 -0500 Subject: [PATCH 13/30] Make isomip plotting more flexible --- compass/ocean/tests/isomip_plus/viz/plot.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/compass/ocean/tests/isomip_plus/viz/plot.py b/compass/ocean/tests/isomip_plus/viz/plot.py index c332b05f59..3c6ac75808 100644 --- a/compass/ocean/tests/isomip_plus/viz/plot.py +++ b/compass/ocean/tests/isomip_plus/viz/plot.py @@ -239,6 +239,8 @@ def __init__(self, inFolder, streamfunctionFolder, outFolder, expt, self.sectionY = sectionY self.showProgress = showProgress + if 'Time' in dsMesh.dims: + dsMesh = dsMesh.isel(Time=0) self.dsMesh = dsMesh self.ds = ds @@ -508,7 +510,9 @@ def plot_horiz_series(self, da, nameInTitle, prefix, oceanDomain, The time indices at which to plot. If not provided, set to all. """ - nTime = self.ds.sizes['Time'] + if 'Time' not in da.dims: + da = da.expand_dims(dim='Time', axis=0) + nTime = da.sizes['Time'] if self.showProgress: widgets = ['plotting {}: '.format(nameInTitle), progressbar.Percentage(), ' ', From c665063517eb11124b3d1d5cad1361951f8d3962 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 24 Oct 2024 11:11:24 -0500 Subject: [PATCH 14/30] Only visualize init for thin film cases --- compass/ocean/tests/isomip_plus/initial_state.py | 6 ++++-- compass/ocean/tests/isomip_plus/isomip_plus_test.py | 3 ++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/initial_state.py b/compass/ocean/tests/isomip_plus/initial_state.py index a525cc3cc5..aacc25d82b 100644 --- a/compass/ocean/tests/isomip_plus/initial_state.py +++ b/compass/ocean/tests/isomip_plus/initial_state.py @@ -40,7 +40,7 @@ class InitialState(Step): Whether the run includes a thin film below grounded ice """ def __init__(self, test_case, resolution, experiment, vertical_coordinate, - time_varying_forcing, thin_film_present): + time_varying_forcing, thin_film_present, include_viz=False): """ Create the step @@ -65,6 +65,7 @@ def __init__(self, test_case, resolution, experiment, vertical_coordinate, Whether the run includes a thin film below grounded ice """ super().__init__(test_case=test_case, name='initial_state') + self.include_viz = include_viz self.resolution = resolution self.experiment = experiment self.vertical_coordinate = vertical_coordinate @@ -88,7 +89,8 @@ def run(self): """ ds, frac = self._compute_initial_condition() self._compute_restoring(ds, frac) - self._plot(ds) + if self.include_viz: + self._plot(ds) def _compute_initial_condition(self): config = self.config diff --git a/compass/ocean/tests/isomip_plus/isomip_plus_test.py b/compass/ocean/tests/isomip_plus/isomip_plus_test.py index cf17731a01..b701816661 100644 --- a/compass/ocean/tests/isomip_plus/isomip_plus_test.py +++ b/compass/ocean/tests/isomip_plus/isomip_plus_test.py @@ -140,7 +140,8 @@ def __init__(self, test_group, resolution, experiment, experiment=experiment, vertical_coordinate=vertical_coordinate, time_varying_forcing=time_varying_forcing, - thin_film_present=thin_film_present)) + thin_film_present=thin_film_present, + include_viz=thin_film_present)) self.add_step( SshAdjustment(test_case=self, resolution=resolution, vertical_coordinate=vertical_coordinate, From f1eaa303d6a293fc849208b856fd217026ee6ef9 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 11 Nov 2024 09:19:28 -0600 Subject: [PATCH 15/30] Add new thin film namelist options --- .../namelist.thin_film.forward_and_ssh_adjust | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust b/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust index 76c2443f7c..047f5db62e 100644 --- a/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust +++ b/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust @@ -1,7 +1,17 @@ config_time_integrator='RK4' -config_dt='00:00:10' +config_thickness_flux_type='upwind' +config_land_ice_flux_mode='pressure_only' config_use_wetting_drying=.true. +config_zero_drying_velocity_ramp=.true. +config_zero_drying_velocity_ramp_hmax=1e1 +config_use_ssh_gradient_wetting_drying=.false. config_prevent_drying=.true. config_drying_min_cell_height=1.0e-3 config_zero_drying_velocity = .true. -config_thickness_flux_type='upwind' +config_use_land_ice_ssh_adjustment = .false. +config_use_land_ice_forcing=.false. +config_land_ice_dynamic_height=3.1e-3 +config_use_land_ice_pressure_relaxation = .true. +config_land_ice_pressure_increment = 10. +config_use_land_ice_overburden_height=.true. +config_land_ice_overburden_height=-3e-3 From debd9f7e264afef3b01376e4e866924b530281e5 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 11 Nov 2024 09:20:10 -0600 Subject: [PATCH 16/30] Add more output vars for wetdry --- compass/ocean/tests/isomip_plus/streams.forward.template | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/compass/ocean/tests/isomip_plus/streams.forward.template b/compass/ocean/tests/isomip_plus/streams.forward.template index 4eff72b570..e6d26d6f52 100644 --- a/compass/ocean/tests/isomip_plus/streams.forward.template +++ b/compass/ocean/tests/isomip_plus/streams.forward.template @@ -42,6 +42,7 @@ + @@ -64,6 +65,7 @@ + + - + From 52bc9bf6f1045fb9cbdf220f395b59313f094ab7 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 11 Nov 2024 09:22:01 -0600 Subject: [PATCH 17/30] Fix viz when land ice fluxes off --- compass/ocean/tests/isomip_plus/viz/plot.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/compass/ocean/tests/isomip_plus/viz/plot.py b/compass/ocean/tests/isomip_plus/viz/plot.py index 3c6ac75808..5199ab1990 100644 --- a/compass/ocean/tests/isomip_plus/viz/plot.py +++ b/compass/ocean/tests/isomip_plus/viz/plot.py @@ -84,6 +84,8 @@ def plot_melt_time_series(self, sshMax=None): mean thermal driving, mean friction velocity """ + if 'timeMonthly_avg_landIceFreshwaterFlux' not in self.ds.keys(): + return rho_fw = 1000. secPerYear = 365 * 24 * 60 * 60 @@ -351,6 +353,8 @@ def plot_melt_rates(self, vmin=-100., vmax=100.): vmin, vmax : float, optional The minimum and maximum values for the colorbar """ + if 'timeMonthly_avg_landIceFreshwaterFlux' not in self.ds.keys(): + return rho_fw = 1000. secPerYear = 365 * 24 * 60 * 60 @@ -369,6 +373,8 @@ def plot_ice_shelf_boundary_variables(self): ice """ + if 'timeMonthly_avg_landIceFreshwaterFlux' not in self.ds.keys(): + return self.plot_horiz_series(self.ds.timeMonthly_avg_landIceHeatFlux, 'heat flux from ocean to ice-ocean interface', prefix='oceanHeatFlux', From 5147eff2abf3552df1efaa87af470feb580b8088 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 14 Mar 2025 12:25:09 -0500 Subject: [PATCH 18/30] Add z-star to isomip+ tests --- compass/ocean/tests/isomip_plus/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/__init__.py b/compass/ocean/tests/isomip_plus/__init__.py index fdbb42bd9d..763d653283 100644 --- a/compass/ocean/tests/isomip_plus/__init__.py +++ b/compass/ocean/tests/isomip_plus/__init__.py @@ -50,7 +50,7 @@ def __init__(self, mpas_core): vertical_coordinate=vertical_coordinate, time_varying_forcing=True, thin_film_present=True)) - for vertical_coordinate in ['sigma', 'single_layer']: + for vertical_coordinate in ['z-star', 'sigma', 'single_layer']: self.add_test_case( IsomipPlusTest( test_group=self, resolution=resolution, @@ -67,7 +67,7 @@ def __init__(self, mpas_core): time_varying_forcing=True, time_varying_load='decreasing', thin_film_present=True)) - for vertical_coordinate in ['sigma']: + for vertical_coordinate in ['sigma', 'z-star']: self.add_test_case( IsomipPlusTest( test_group=self, resolution=resolution, From 376fd13ebe54300f5f9e3f4e60c20edf2531709f Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 14 Mar 2025 12:26:52 -0500 Subject: [PATCH 19/30] Add more plot types to performance step of isomip+ --- compass/ocean/tests/isomip_plus/forward.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/compass/ocean/tests/isomip_plus/forward.py b/compass/ocean/tests/isomip_plus/forward.py index e6689f8130..57bad7e632 100644 --- a/compass/ocean/tests/isomip_plus/forward.py +++ b/compass/ocean/tests/isomip_plus/forward.py @@ -162,11 +162,14 @@ def run(self): plot_folder = f'{self.work_dir}/plots' if os.path.exists(plot_folder): shutil.rmtree(plot_folder) + min_column_thickness = self.config.getfloat('isomip_plus', + 'min_column_thickness') dsMesh = xarray.open_dataset(os.path.join(self.work_dir, 'init.nc')) ds = xarray.open_dataset(os.path.join(self.work_dir, 'output.nc')) - + ds['landIceDraft'] = \ + dsMesh.landIceDraft.expand_dims('Time', axis=0) section_y = self.config.getfloat('isomip_plus_viz', 'section_y') # show progress only if we're not writing to a log file show_progress = self.log_filename is None @@ -176,6 +179,11 @@ def run(self): dsMesh=dsMesh, ds=ds, expt=self.experiment, showProgress=show_progress) + bottomDepth = ds.bottomDepth.expand_dims(dim='Time', axis=0) + plotter.plot_horiz_series(ds.ssh + bottomDepth, + 'H', 'H', True, + vmin=min_column_thickness, vmax=700, + cmap_set_under='r', cmap_scale='log') plotter.plot_horiz_series(ds.ssh, 'ssh', 'ssh', True, vmin=-700, vmax=0) plotter.plot_3d_field_top_bot_section( @@ -185,6 +193,7 @@ 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, From bdb27d692496ff8260f72973451d4065cbb916f0 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 11 Nov 2024 09:21:31 -0600 Subject: [PATCH 20/30] Add more plot types to viz step --- .../ocean/tests/isomip_plus/viz/__init__.py | 279 +++++++++++++----- compass/ocean/tests/isomip_plus/viz/plot.py | 115 ++++++-- 2 files changed, 291 insertions(+), 103 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/viz/__init__.py b/compass/ocean/tests/isomip_plus/viz/__init__.py index 88be3341e9..2748598493 100644 --- a/compass/ocean/tests/isomip_plus/viz/__init__.py +++ b/compass/ocean/tests/isomip_plus/viz/__init__.py @@ -55,6 +55,7 @@ def run(self): plot_streamfunctions = section.getboolean('plot_streamfunctions') plot_haney = section.getboolean('plot_haney') frames_per_second = section.getint('frames_per_second') + # output_slices = section.getint('output_slices') movie_format = section.get('movie_format') section_y = section.getfloat('section_y') @@ -65,40 +66,133 @@ def run(self): show_progress = self.log_filename is None expt = self.experiment - if self.tidal_forcing: - sim_dir = '../performance' - else: - sim_dir = '../simulation' - if not os.path.exists(f'{sim_dir}/timeSeriesStatsMonthly.0001-01-01.nc'): # noqa: E501 - sim_dir = '../performance' + sim_dir = '../simulation' streamfunction_dir = '../streamfunction' out_dir = '.' dsMesh = xarray.open_dataset(f'{sim_dir}/init.nc') + maxBottomDepth = dsMesh.bottomDepth.max().values dsOut = xarray.open_dataset(f'{sim_dir}/output.nc') - + dsIce = xarray.open_dataset(f'{sim_dir}/land_ice_fluxes.nc') + if dsOut.sizes['Time'] < 50: + print('Cropping dsOut') + dsOut = dsOut.isel(Time=slice(-31, -1)) + nTime = dsOut.sizes['Time'] + if nTime != dsIce.sizes['Time']: + print('Cropping dsIce to match dsOut') + dsIce = dsIce.isel(Time=slice(-nTime - 1, -1)) + dsOut['landIceDraft'] = dsIce.landIceDraft plotter = MoviePlotter(inFolder=sim_dir, streamfunctionFolder=streamfunction_dir, outFolder=f'{out_dir}/plots', expt=expt, sectionY=section_y, dsMesh=dsMesh, ds=dsOut, showProgress=show_progress) - - if 'time_varying' in expt: - plotter.plot_horiz_series(dsOut.ssh, 'ssh', 'ssh', True, - cmap='cmo.curl') - delice = dsOut.landIcePressure - dsOut.landIcePressure[0, :] - plotter.plot_horiz_series(delice, 'delLandIcePressure', - 'delLandIcePressure', True, - cmap='cmo.curl') - plotter.plot_horiz_series(dsOut.velocityX[:, :, 0], 'u', 'u', True, - cmap='cmo.balance', vmin=-5e-1, vmax=5e-1) - plotter.plot_horiz_series(dsOut.velocityY[:, :, 0], 'v', 'v', True, - cmap='cmo.balance', vmin=-5e-1, vmax=5e-1) + plotter.plot_3d_field_top_bot_section(dsOut.velocityX, 'u', 'u', True, + cmap='cmo.balance', + vmin=-5e-1, vmax=5e-1) + plotter.plot_3d_field_top_bot_section(dsOut.velocityY, 'v', 'v', True, + cmap='cmo.balance', + vmin=-5e-1, vmax=5e-1) + + if dsOut.sizes['Time'] > 24: + days = 30 + samples_per_day = 1 + dsAve = dsOut.isel(Time=slice(-int(days * samples_per_day), -1)) + velocityMagSq = (dsAve.velocityX.isel(nVertLevels=0)**2. + + dsAve.velocityX.isel(nVertLevels=0)**2.) + rmsVelocity = xarray.DataArray( + data=np.sqrt(velocityMagSq.mean('Time').values), + dims=['nCells']).expand_dims(dim='Time', axis=0) + plotter.plot_horiz_series(rmsVelocity, + 'rmsVelocity', 'rmsVelocity', + False, vmin=1.e-2, vmax=1.e0, + cmap='cmo.speed', cmap_scale='log') + + cavityMask = dsMesh.landIceFloatingMask == 1 + if 'xIsomipCell' in dsMesh.keys(): + xCell = dsMesh.xIsomipCell / 1.e3 + else: + xCell = dsMesh.xCell / 1.e3 + # salinity = dsIce.landIceInterfaceSalinity + salinity = dsOut.salinity.isel(nVertLevels=0) + pressure = dsIce.landIcePressure + # temperature = dsIce.landIceInterfaceTemperature + temperature = dsOut.temperature.isel(nVertLevels=0) + coeff_0 = 6.22e-2 + coeff_S = -5.63e-2 + coeff_p = -7.43e-8 + coeff_pS = -1.74e-10 + coeff_mushy = 0. + ocn_freezing_temperature = \ + coeff_0 + coeff_S * salinity + \ + coeff_p * pressure + \ + coeff_pS * pressure * salinity + \ + coeff_mushy * salinity / (1.0 - salinity / 1e3) + dsIce['thermalForcing'] = temperature - ocn_freezing_temperature + + dx = 2. + xbins = np.arange(xCell.isel(nCells=cavityMask).min(), + xCell.isel(nCells=cavityMask).max(), dx) + if dsIce.sizes['Time'] > 14: + dsIce = dsIce.isel(Time=slice(-13, -12)) + for time_slice in range(dsIce.sizes['Time']): + title = dsIce.xtime.isel(Time=time_slice).values + meltTransect = np.zeros_like(xbins) + landIceDraftTransect = np.zeros_like(xbins) + frictionVelocityTransect = np.zeros_like(xbins) + temperatureTransect = np.zeros_like(xbins) + thermalForcingTransect = np.zeros_like(xbins) + rmsVelocityTransect = np.zeros_like(xbins) + for i, xmin in enumerate(xbins): + binMask = np.logical_and(cavityMask, + np.logical_and(xCell >= xmin, + xCell < xmin + dx)) + dsTransect = dsIce.isel(nCells=binMask, Time=time_slice) + meltTransect[i] = dsTransect.landIceFreshwaterFlux.mean() + landIceDraftTransect[i] = dsTransect.landIceDraft.mean() + frictionVelocityTransect[i] = \ + dsTransect.landIceFrictionVelocity.mean() + temperatureTransect[i] = \ + dsTransect.landIceInterfaceTemperature.mean() + thermalForcingTransect[i] = dsTransect.thermalForcing.mean() + if dsOut.sizes['Time'] > 24: + rmsVelocityTransect[i] = np.mean(rmsVelocity[binMask]) + + fig, ax = plt.subplots(4, 1, sharex=True, figsize=(4, 8)) + color = 'darkgreen' + x = xbins - xbins[0] + ax[0].set_title(title) + ax[0].plot(x, landIceDraftTransect, color=color) + ax[0].set_ylabel('Land ice draft (m)') + secPerYear = 365 * 24 * 60 * 60 + density = 1026. + ax[1].plot(x, meltTransect * secPerYear / density, color=color) + ax[1].set_ylabel('Melt rate (m/yr)') + ax[2].plot(x, frictionVelocityTransect, color=color) + ax[2].set_ylabel('Friction velocity (m/s)') + ax[3].plot(x, thermalForcingTransect, color=color) + ax[3].set_ylabel(r'ThermalForcing ($^{\circ}$C)') + ax[3].set_xlabel('Distance along ice flow (km)') + plt.tight_layout() + plt.savefig(f'meanTransect_tend{time_slice:03g}.png', dpi=600, + bbox_inches='tight', transparent=True) + plt.close(fig) + + plotter.plot_horiz_series(dsOut.ssh, 'ssh', 'ssh', True, + vmin=-maxBottomDepth, vmax=0) + if 'layerThickness' in dsOut.keys(): + plotter.plot_layer_interfaces() + delssh = dsOut.ssh - dsMesh.ssh.isel(Time=0) + plotter.plot_horiz_series(delssh, 'ssh - ssh_init', 'delssh_init', + True, cmap='cmo.curl', vmin=-1, vmax=1) plotter.plot_horiz_series(dsOut.ssh + dsMesh.bottomDepth, 'H', 'H', True, vmin=min_column_thickness + 1e-10, - vmax=700, cmap_set_under='r', + vmax=500, cmap_set_under='k', cmap_scale='log') + plotter.plot_horiz_series(dsOut.landIcePressure, + 'landIcePressure', 'landIcePressure', + True, vmin=1e5, vmax=1e7, cmap_scale='log') plotter.plot_horiz_series(dsOut.surfacePressure, 'surfacePressure', 'surfacePressure', True, vmin=1e5, vmax=1e7, cmap_scale='log') @@ -106,74 +200,115 @@ def run(self): plotter.plot_horiz_series(delSurfacePressure, 'delSurfacePressure', 'delSurfacePressure', True, vmin=1e5, vmax=1e7, cmap_scale='log') - - if 'tidal' in expt: - delssh = dsOut.ssh - dsOut.ssh[0, :] - plotter.plot_horiz_series(delssh, 'delssh', 'delssh', True, - cmap='cmo.curl', vmin=-1, vmax=1) + delice = dsOut.landIcePressure - dsMesh.landIcePressure.isel(Time=0) + plotter.plot_horiz_series(delice, 'delLandIcePressure', + 'delLandIcePressure', True, + vmin=-1e6, vmax=1e6, cmap='cmo.balance') + # Cannot currently plot because it is on edges + for t in range(dsOut.sizes['Time']): + fig = plt.figure() + sc = plt.scatter( + dsOut.xEdge, dsOut.yEdge, 10, + dsOut.wettingVelocityFactor.isel(nVertLevels=0, Time=t), + vmin=0, vmax=1, cmap='cmo.thermal') + current_cmap = sc.get_cmap() + current_cmap.set_under('k') + current_cmap.set_over('r') + plt.colorbar() + fig.gca().set_aspect('equal') + plt.savefig(f'wettingVelocityFactor_{t:02g}.png') + plt.close(fig) wct = dsOut.ssh + dsMesh.bottomDepth - idx_thin = np.logical_and(wct[0, :] > 1e-1, - wct[0, :] < 1) - wct_thin = wct[:, idx_thin] - wct_mean = wct_thin.mean(dim='nCells').values - time = dsOut.daysSinceStartOfSim.values - plt.figure() - plt.plot(time, wct_mean, '.') - plt.xlabel('Time (days)') - plt.ylabel('Mean thickness of thin film (m)') - plt.savefig('wct_thin_t.png') - plt.close() - plotter.plot_horiz_series(wct, 'H', 'H', True, vmin=min_column_thickness + 1e-10, - vmax=700, cmap_set_under='r', + vmax=700, cmap_set_under='k', cmap_scale='log') - if os.path.exists(f'{sim_dir}/timeSeriesStatsMonthly.0001-01-01.nc'): + ds = xarray.open_mfdataset( + '{}/timeSeriesStatsMonthly*.nc'.format(sim_dir), + concat_dim='Time', combine='nested') + + if plot_haney: + _compute_and_write_haney_number(dsMesh, ds, out_dir, + showProgress=show_progress) + + tsPlotter = TimeSeriesPlotter(inFolder=sim_dir, + outFolder='{}/plots'.format(out_dir), + expt=expt) + tsPlotter.plot_melt_time_series() + tsPlotter.plot_melt_time_series(wctMax=20.) + tsPlotter = TimeSeriesPlotter( + inFolder=sim_dir, + outFolder='{}/timeSeriesBelow300m'.format(out_dir), + expt=expt) + tsPlotter.plot_melt_time_series(sshMax=-300.) + + mPlotter = MoviePlotter(inFolder=sim_dir, + streamfunctionFolder=streamfunction_dir, + outFolder='{}/plots'.format(out_dir), + expt=expt, sectionY=section_y, + dsMesh=dsMesh, ds=ds, + showProgress=show_progress) + + mPlotter.plot_layer_interfaces() + + if plot_streamfunctions: + mPlotter.plot_barotropic_streamfunction() + mPlotter.plot_overturning_streamfunction() + + if plot_haney: + mPlotter.plot_haney_number(haneyFolder=out_dir) + + mPlotter.plot_melt_rates() + mPlotter.plot_ice_shelf_boundary_variables() + mPlotter.plot_temperature() + mPlotter.plot_salinity() + mPlotter.plot_potential_density() + + mPlotter.images_to_movies(outFolder=f'{out_dir}/movies', + framesPerSecond=frames_per_second, + extension=movie_format) + + if os.path.exists('../simulation/land_ice_fluxes.nc'): ds = xarray.open_mfdataset( - '{}/timeSeriesStatsMonthly*.nc'.format(sim_dir), + '../simulation/land_ice_fluxes.nc', concat_dim='Time', combine='nested') - - if plot_haney: - _compute_and_write_haney_number(dsMesh, ds, out_dir, - showProgress=show_progress) - - tsPlotter = TimeSeriesPlotter(inFolder=sim_dir, - outFolder='{}/plots'.format(out_dir), - expt=expt) - tsPlotter.plot_melt_time_series() - tsPlotter = TimeSeriesPlotter( - inFolder=sim_dir, - outFolder='{}/timeSeriesBelow300m'.format(out_dir), - expt=expt) - tsPlotter.plot_melt_time_series(sshMax=-300.) - + if ds.sizes['Time'] != dsOut.sizes['Time']: + # assumes ds.sizes['Time'] > dsOut.sizes['Time'] + ds = ds.isel(Time=slice(-31, -1)) + dsOut = dsOut.isel(Time=slice(-31, -1)) + if ds.sizes['Time'] != dsOut.sizes['Time']: + print('error in time selection') + return + if ds.sizes['Time'] > 50: + return + ds['layerThickness'] = dsOut.layerThickness mPlotter = MoviePlotter(inFolder=sim_dir, streamfunctionFolder=streamfunction_dir, outFolder='{}/plots'.format(out_dir), expt=expt, sectionY=section_y, dsMesh=dsMesh, ds=ds, showProgress=show_progress) - - mPlotter.plot_layer_interfaces() - - if plot_streamfunctions: - mPlotter.plot_barotropic_streamfunction() - mPlotter.plot_overturning_streamfunction() - - if plot_haney: - mPlotter.plot_haney_number(haneyFolder=out_dir) - + mPlotter.plot_horiz_series(ds.landIceFrictionVelocity, + 'landIceFrictionVelocity', + 'frictionVelocity', + True) + mPlotter.plot_horiz_series(ds.landIceFreshwaterFlux, + 'landIceFreshwaterFlux', 'melt', + True) + mPlotter.plot_horiz_series(ds.landIceFloatingFraction, + 'landIceFloatingFraction', + 'landIceFloatingFraction', + True, vmin=1e-16, vmax=1, + cmap_set_under='k') + mPlotter.plot_horiz_series(ds.landIceDraft, + 'landIceDraft', 'landIceDraft', + True, vmin=-maxBottomDepth, vmax=0.) mPlotter.plot_melt_rates() mPlotter.plot_ice_shelf_boundary_variables() - mPlotter.plot_temperature() - mPlotter.plot_salinity() - mPlotter.plot_potential_density() - - mPlotter.images_to_movies(outFolder='{}/movies'.format(out_dir), - framesPerSecond=frames_per_second, - extension=movie_format) + else: + print(f'didnt find {sim_dir}/land_ice_fluxes.nc') def file_complete(ds, fileName): diff --git a/compass/ocean/tests/isomip_plus/viz/plot.py b/compass/ocean/tests/isomip_plus/viz/plot.py index 5199ab1990..bc67d4b449 100644 --- a/compass/ocean/tests/isomip_plus/viz/plot.py +++ b/compass/ocean/tests/isomip_plus/viz/plot.py @@ -9,7 +9,7 @@ import progressbar import xarray from matplotlib.collections import PatchCollection -from matplotlib.colors import LogNorm +from matplotlib.colors import LogNorm, SymLogNorm from matplotlib.patches import Polygon @@ -77,7 +77,7 @@ def __init__(self, inFolder='.', outFolder='plots', expt='Ocean0', plt.switch_backend('Agg') - def plot_melt_time_series(self, sshMax=None): + def plot_melt_time_series(self, sshMax=None, wctMax=None): """ Plot a series of image for each of several variables related to melt at the ice shelf-ocean interface: mean melt rate, total melt flux, @@ -89,21 +89,28 @@ def plot_melt_time_series(self, sshMax=None): rho_fw = 1000. secPerYear = 365 * 24 * 60 * 60 + suffix = '' areaCell = self.dsMesh.areaCell iceMask = self.ds.timeMonthly_avg_landIceFraction meltFlux = self.ds.timeMonthly_avg_landIceFreshwaterFlux if sshMax is not None: ssh = self.ds.timeMonthly_avg_ssh iceMask = iceMask.where(ssh < sshMax) + suffix = f'_sshMax{sshMax}' + elif wctMax is not None: + H = self.ds.timeMonthly_avg_layerThickness.sum(dim='nVertLevels') + iceMask = iceMask.where(H < wctMax) + suffix = f'_wctMax{wctMax}' totalMeltFlux = (meltFlux * areaCell * iceMask).sum(dim='nCells') totalArea = (areaCell * iceMask).sum(dim='nCells') meanMeltRate = totalMeltFlux / totalArea / rho_fw * secPerYear self.plot_time_series(meanMeltRate, 'mean melt rate', 'meanMeltRate', 'm/yr') + print(meanMeltRate.values) self.plot_time_series(1e-6 * totalMeltFlux, 'total melt flux', - 'totalMeltFlux', 'kT/yr') + f'totalMeltFlux{suffix}', 'kT/yr') prefix = 'timeMonthly_avg_landIceBoundaryLayerTracers_' boundary_layer_temperature = \ @@ -115,13 +122,13 @@ def plot_melt_time_series(self, sshMax=None): da = (da * areaCell * iceMask).sum(dim='nCells') / totalArea self.plot_time_series(da, 'mean thermal driving', - 'meanThermalDriving', 'deg C') + f'meanThermalDriving{suffix}', 'deg C') da = self.ds.timeMonthly_avg_landIceFrictionVelocity da = (da * areaCell * iceMask).sum(dim='nCells') / totalArea self.plot_time_series(da, 'mean friction velocity', - 'meanFrictionVelocity', 'm/s') + f'meanFrictionVelocity{suffix}', 'm/s') def plot_time_series(self, da, nameInTitle, prefix, units=None, figsize=(12, 6), color=None, overwrite=True): @@ -246,7 +253,7 @@ def __init__(self, inFolder, streamfunctionFolder, outFolder, expt, self.dsMesh = dsMesh self.ds = ds - landIceMask = self.dsMesh.landIceMask.isel(Time=0) > 0 + landIceMask = self.dsMesh.landIceMask > 0 self.oceanMask = self.dsMesh.maxLevelCell - 1 >= 0 self.cavityMask = numpy.logical_and(self.oceanMask, landIceMask) @@ -344,7 +351,7 @@ def plot_overturning_streamfunction(self, vmin=-0.3, vmax=0.3): if self.showProgress: bar.finish() - def plot_melt_rates(self, vmin=-100., vmax=100.): + def plot_melt_rates(self, vmin=0., vmax=50.): """ Plot a series of image of the melt rate @@ -363,7 +370,7 @@ def plot_melt_rates(self, vmin=-100., vmax=100.): self.plot_horiz_series(da, 'melt rate', prefix='meltRate', oceanDomain=False, units='m/yr', vmin=vmin, - vmax=vmax, cmap='cmo.curl') + vmax=vmax, cmap='cmo.amp') def plot_ice_shelf_boundary_variables(self): """ @@ -401,7 +408,7 @@ def plot_ice_shelf_boundary_variables(self): self.plot_horiz_series(da, 'thermal driving', prefix='thermalDriving', oceanDomain=False, units='deg C', - vmin=-2, vmax=2, cmap='cmo.thermal') + vmin=-2, vmax=2, cmap='cmo.curl') da = boundary_layer_salinity - interface_salinity self.plot_horiz_series(da, 'haline driving', @@ -412,8 +419,9 @@ def plot_ice_shelf_boundary_variables(self): self.plot_horiz_series(self.ds.timeMonthly_avg_landIceFrictionVelocity, 'friction velocity', prefix='frictionVelocity', - oceanDomain=True, units='m/s', - vmin=0, vmax=0.05, cmap='cmo.speed') + oceanDomain=False, units='m/s', + vmin=1.e-4, vmax=1.e-1, cmap='cmo.speed', + cmap_scale='log') def plot_temperature(self): """ @@ -438,8 +446,8 @@ def plot_salinity(self): self.plot_3d_field_top_bot_section(da, nameInTitle='salinity', prefix='Salinity', units='PSU', - vmin=33.8, vmax=34.7, - cmap='cmo.haline') + vmin=1., vmax=34.6, + cmap='cmo.haline', cmap_scale='log') def plot_potential_density(self): """ @@ -478,7 +486,7 @@ def plot_horiz_series(self, da, nameInTitle, prefix, oceanDomain, units=None, vmin=None, vmax=None, cmap=None, cmap_set_under=None, cmap_set_over=None, cmap_scale='linear', time_indices=None, - figsize=(9, 3)): + figsize=(9, 3), contour_field=None): """ Plot a series of image of a given variable @@ -533,6 +541,8 @@ def plot_horiz_series(self, da, nameInTitle, prefix, oceanDomain, for tIndex in time_indices: self.update_date(tIndex) field = da.isel(Time=tIndex).values + if contour_field is not None: + contour_field = contour_field.isel(Time=tIndex).values outFileName = '{}/{}/{}_{:04d}.png'.format( self.outFolder, prefix, prefix, tIndex + 1) if units is None: @@ -544,7 +554,8 @@ def plot_horiz_series(self, da, nameInTitle, prefix, oceanDomain, vmax=vmax, cmap=cmap, cmap_set_under=cmap_set_under, cmap_set_over=cmap_set_over, - cmap_scale=cmap_scale, figsize=figsize) + cmap_scale=cmap_scale, figsize=figsize, + contour_field=contour_field) if self.showProgress: bar.update(tIndex + 1) if self.showProgress: @@ -552,7 +563,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_scale='linear', + cmap_set_under=None, cmap_set_over=None): """ Plot a series of images of a given 3D variable showing the value @@ -608,6 +620,7 @@ 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_scale=cmap_scale, cmap_set_under=cmap_set_under, cmap_set_over=cmap_set_over) @@ -630,7 +643,8 @@ def plot_3d_field_top_bot_section(self, da, nameInTitle, prefix, self.plot_horiz_series(daBot, 'bot {}'.format(nameInTitle), 'bot{}'.format(prefix), oceanDomain=True, - vmin=vmin, vmax=vmax, cmap=cmap) + vmin=vmin, vmax=vmax, cmap_scale=cmap_scale, + cmap=cmap) daSection = da.isel(nCells=self.sectionCellIndices) @@ -659,7 +673,8 @@ def plot_3d_field_top_bot_section(self, da, nameInTitle, prefix, self._plot_vert_field(self.X, self.Z[tIndex, :, :], field, title=title, outFileName=outFileName, - vmin=vmin, vmax=vmax, cmap=cmap) + vmin=vmin, vmax=vmax, cmap=cmap, + show_boundaries=False, cmap_scale=cmap_scale) if self.showProgress: bar.update(tIndex + 1) if self.showProgress: @@ -718,6 +733,7 @@ def plot_layer_interfaces(self, figsize=(9, 5)): plt.plot(1e-3 * X[z_index, :], Z[z_index, :], 'k') plt.plot(1e-3 * X[0, :], Z[0, :], 'b') plt.plot(1e-3 * X[0, :], self.zBotSection, 'g') + plt.plot(1e-3 * X[0, :], self.landIceDraft, 'r') ax.autoscale(tight=True) x1, x2, y1, y2 = 420, 470, -650, -520 @@ -729,6 +745,7 @@ def plot_layer_interfaces(self, figsize=(9, 5)): axins.plot(1e-3 * X[z_index, :], Z[z_index, :], 'k') axins.plot(1e-3 * X[0, :], Z[0, :], 'b') axins.plot(1e-3 * X[0, :], self.zBotSection, 'g') + axins.plot(1e-3 * X[0, :], self.landIceDraft, 'r') axins.set_xlim(x1, x2) axins.set_ylim(y1, y2) axins.set_xticklabels([]) @@ -793,7 +810,7 @@ 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_set_over=None, - cmap_scale='linear'): + cmap_scale='linear', contour_field=None): try: os.makedirs(os.path.dirname(outFileName)) @@ -824,10 +841,21 @@ def _plot_horiz_field(self, field, title, outFileName, oceanDomain=True, if cmap_scale == 'log': localPatches.set_norm(LogNorm(vmin=max(1e-10, vmin), vmax=vmax, clip=False)) + elif cmap_scale == 'symlog': + localPatches.set_norm(SymLogNorm(vmin=vmin, vmax=vmax, clip=False, + linthresh=vmax / 1e2)) plt.figure(figsize=figsize) ax = plt.subplot(111) ax.add_collection(localPatches) + if contour_field is not None: + dsMesh = self.dsMesh + ax.tricontour(dsMesh.xCell / 1.e3, dsMesh.yCell / 1.e3, + contour_field, + colors='k', levels=[1.1e-2 + 1.e-8]) + ax.tricontour(dsMesh.xCell / 1.e3, dsMesh.yCell / 1.e3, + contour_field, + colors='grey', levels=numpy.arange(50., 700., 50.)) plt.colorbar(localPatches, extend='both') plt.axis([0, 500, 0, 1000]) ax.set_aspect('equal') @@ -839,7 +867,7 @@ def _plot_horiz_field(self, field, title, outFileName, oceanDomain=True, def _plot_vert_field(self, inX, inZ, field, title, outFileName, vmin=None, vmax=None, figsize=(9, 5), cmap=None, - show_boundaries=True): + show_boundaries=True, cmap_scale='linear'): try: os.makedirs(os.path.dirname(outFileName)) except OSError: @@ -850,24 +878,40 @@ def _plot_vert_field(self, inX, inZ, field, title, outFileName, plt.figure(figsize=figsize) ax = plt.subplot(111) - if show_boundaries: - z_mask = numpy.ones(self.X.shape) - z_mask[0:-1, 0:-1] *= numpy.where(self.sectionMask, 1., numpy.nan) + z_mask = numpy.ones(self.X.shape) + z_mask[0:-1, 0:-1] *= numpy.where(self.sectionMask, 1., numpy.nan) - tIndex = 0 - Z = numpy.array(self.Z[tIndex, :, :]) - Z *= z_mask - X = self.X + tIndex = 0 + Z = numpy.array(self.Z[tIndex, :, :]) + Z *= z_mask + X = self.X - plt.fill_between(1e-3 * X[0, :], self.zBotSection, y2=0, - facecolor='lightsteelblue', zorder=2) - plt.fill_between(1e-3 * X[0, :], self.zBotSection, y2=-750, - facecolor='grey', zorder=1) + plt.fill_between(1e-3 * X[0, :], self.zBotSection, y2=0, + facecolor='lightsteelblue', zorder=2) + plt.fill_between(1e-3 * X[0, :], self.zBotSection, y2=-750, + facecolor='grey', zorder=1) + if show_boundaries: for z_index in range(1, X.shape[0]): plt.plot(1e-3 * X[z_index, :], Z[z_index, :], 'k', zorder=4) plt.pcolormesh(1e-3 * inX, inZ, field, vmin=vmin, vmax=vmax, cmap=cmap, zorder=3) plt.colorbar() + x1, x2, y1, y2 = 455, 470, -640, -540 + axins = ax.inset_axes([0.01, 0.6, 0.3, 0.39]) + axins.fill_between(1e-3 * X[0, :], self.zBotSection, y2=0, + facecolor='lightsteelblue', zorder=2) + axins.fill_between(1e-3 * X[0, :], self.zBotSection, y2=-750, + facecolor='grey', zorder=1) + if show_boundaries: + for z_index in range(1, X.shape[0]): + axins.plot(1e-3 * X[z_index, :], Z[z_index, :], 'k', zorder=4) + axins.pcolormesh(1e-3 * inX, inZ, field, vmin=vmin, vmax=vmax, + cmap=cmap, zorder=3) + axins.set_xlim(x1, x2) + axins.set_ylim(y1, y2) + axins.set_xticklabels([]) + axins.set_yticklabels([]) + ax.indicate_inset_zoom(axins, edgecolor="black") ax.autoscale(tight=True) plt.ylim([numpy.amin(inZ), 20]) plt.xlim([400, 800]) @@ -918,6 +962,15 @@ def _compute_section_x_z(self): layerThickness[:, zIndex]) self.Z[tIndex, zIndex, :] = self.Z[tIndex, zIndex + 1, :] + \ layerThicknessSection + if 'timeMonthly_avg_landIceDraft' in self.ds: + var = 'timeMonthly_avg_landIceDraft' + else: + var = 'landIceDraft' + landIceDraft = self.ds[var].isel( + Time=tIndex, nCells=self.sectionCellIndices) + landIceDraft = landIceDraft.values * self.sectionMask[0, :].T + landIceDraft = numpy.nan_to_num(landIceDraft) + self.landIceDraft = _interp_extrap_corner(landIceDraft) def _compute_cell_patches(dsMesh, mask): From 8b7f55a7ba44761456a10290c9adae523085bff4 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 14 Mar 2025 13:40:00 -0500 Subject: [PATCH 21/30] Change thin film thickness --- compass/ocean/tests/isomip_plus/isomip_plus.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/ocean/tests/isomip_plus/isomip_plus.cfg b/compass/ocean/tests/isomip_plus/isomip_plus.cfg index 65caf293aa..ef5dac0bd3 100644 --- a/compass/ocean/tests/isomip_plus/isomip_plus.cfg +++ b/compass/ocean/tests/isomip_plus/isomip_plus.cfg @@ -64,7 +64,7 @@ minimum_levels = 3 min_layer_thickness = 0.0 # Minimum thickness of the initial ocean column (to prevent 'drying') -min_column_thickness = 1e-3 +min_column_thickness = 1.1e-2 # Minimum fraction of a cell that contains ocean (as opposed to land or # grounded land ice) in order for it to be an active ocean cell. From a879f88e0b6676bd18a4db3349ffd97379d267f7 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 14 Mar 2025 13:40:41 -0500 Subject: [PATCH 22/30] Change defaults for time varying tests --- compass/ocean/tests/isomip_plus/isomip_plus_test.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/isomip_plus_test.py b/compass/ocean/tests/isomip_plus/isomip_plus_test.py index b701816661..9cbdae286a 100644 --- a/compass/ocean/tests/isomip_plus/isomip_plus_test.py +++ b/compass/ocean/tests/isomip_plus/isomip_plus_test.py @@ -210,9 +210,9 @@ def configure(self): config.set('isomip_plus', 'min_column_thickness', '1e-3') if time_varying_load == 'increasing': - config.set('isomip_plus_forcing', 'scales', '1.0, 2.0, 2.0') + config.set('isomip_plus_forcing', 'scales', '1.0, 1.05, 1.1') if time_varying_load == 'decreasing': - config.set('isomip_plus_forcing', 'scales', '1.0, 0.0, 0.0') + config.set('isomip_plus_forcing', 'scales', '1.0, 0.95, 0.9') if experiment in ['Ocean0', 'Ocean2', 'Ocean3']: # warm initial conditions @@ -270,7 +270,7 @@ def configure(self): if vertical_coordinate == 'sigma': if time_varying_load in ['increasing', 'decreasing']: - config.set('vertical_grid', 'vert_levels', '3') + config.set('vertical_grid', 'vert_levels', '10') else: # default to 10 vertical levels instead of 36 config.set('vertical_grid', 'vert_levels', '10') @@ -304,10 +304,10 @@ def _get_time_steps(config, resolution, thin_film_present, tidal_forcing): dt_per_km = 120. if tidal_forcing or thin_film_present: - dt_per_km = min(dt_per_km, 10.) + dt_per_km = min(dt_per_km, 6.) config.set('isomip_plus', 'dt_per_km', f'{dt_per_km}') - dt_btr_per_km = min(dt_per_km / 5., 5.) + dt_btr_per_km = min(dt_per_km / 3., 3.) config.set('isomip_plus', 'dt_btr_per_km', f'{dt_btr_per_km}') return From 205cb764973d2970395f2bb83801550bf1405e36 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 14 Mar 2025 13:41:50 -0500 Subject: [PATCH 23/30] Change namelist options for thin film test --- .../namelist.thin_film.forward_and_ssh_adjust | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust b/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust index 047f5db62e..e70be68f35 100644 --- a/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust +++ b/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust @@ -1,17 +1,14 @@ config_time_integrator='RK4' config_thickness_flux_type='upwind' -config_land_ice_flux_mode='pressure_only' +config_land_ice_flux_mode='standalone' +config_vert_advection_method='flux-form' config_use_wetting_drying=.true. -config_zero_drying_velocity_ramp=.true. -config_zero_drying_velocity_ramp_hmax=1e1 +config_drying_safety_height=1e-10 +config_drying_min_cell_height=1.0e-3 +config_zero_drying_velocity_ramp=.false. +config_zero_drying_velocity_ramp_hmin=1e-3 +config_zero_drying_velocity_ramp_hmax=1e-1 config_use_ssh_gradient_wetting_drying=.false. config_prevent_drying=.true. -config_drying_min_cell_height=1.0e-3 config_zero_drying_velocity = .true. -config_use_land_ice_ssh_adjustment = .false. -config_use_land_ice_forcing=.false. -config_land_ice_dynamic_height=3.1e-3 -config_use_land_ice_pressure_relaxation = .true. -config_land_ice_pressure_increment = 10. -config_use_land_ice_overburden_height=.true. -config_land_ice_overburden_height=-3e-3 +config_use_implicit_top_drag = .true. From bd5836d2abd6c5df3e626518c3be0fbe127818d9 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 14 Mar 2025 16:03:27 -0500 Subject: [PATCH 24/30] Change default time varying forcing --- compass/ocean/tests/isomip_plus/isomip_plus.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compass/ocean/tests/isomip_plus/isomip_plus.cfg b/compass/ocean/tests/isomip_plus/isomip_plus.cfg index ef5dac0bd3..bcaf46fe43 100644 --- a/compass/ocean/tests/isomip_plus/isomip_plus.cfg +++ b/compass/ocean/tests/isomip_plus/isomip_plus.cfg @@ -126,7 +126,7 @@ dates = 0001-01-01_00:00:00, 0002-01-01_00:00:00, 0003-01-01_00:00:00 # the amount by which the initial landIcePressure and landIceDraft are scaled # at each date -scales = 0.1, 1.0, 1.0 +scales = 1.0, 0.9, 0.9 # config options for computing ISOMIP+ streamfunctions [isomip_plus_streamfunction] From d65be6a9a15bc4a0fcc5f6215a94aca1225532e0 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Sun, 16 Mar 2025 15:39:59 -0500 Subject: [PATCH 25/30] Fixup when implicit top drag is used --- compass/ocean/tests/isomip_plus/namelist.forward_and_ssh_adjust | 1 + .../isomip_plus/namelist.single_layer.forward_and_ssh_adjust | 1 + .../tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust | 1 - 3 files changed, 2 insertions(+), 1 deletion(-) diff --git a/compass/ocean/tests/isomip_plus/namelist.forward_and_ssh_adjust b/compass/ocean/tests/isomip_plus/namelist.forward_and_ssh_adjust index 939e440cb3..4c977000b7 100644 --- a/compass/ocean/tests/isomip_plus/namelist.forward_and_ssh_adjust +++ b/compass/ocean/tests/isomip_plus/namelist.forward_and_ssh_adjust @@ -22,6 +22,7 @@ config_land_ice_flux_explicit_topDragCoeff = 2.5e-3 config_land_ice_flux_tidal_Jourdain_U0 = 1e-2 config_land_ice_flux_jenkins_heat_transfer_coefficient = 0.0194 config_land_ice_flux_jenkins_salt_transfer_coefficient = 0.00055428571 +config_use_implicit_top_drag = .true. config_implicit_constant_bottom_drag_coeff = 2.5e-3 config_eos_type = 'linear' config_eos_linear_alpha = 0.03836 diff --git a/compass/ocean/tests/isomip_plus/namelist.single_layer.forward_and_ssh_adjust b/compass/ocean/tests/isomip_plus/namelist.single_layer.forward_and_ssh_adjust index 2e6e126dc0..d1eeee1657 100644 --- a/compass/ocean/tests/isomip_plus/namelist.single_layer.forward_and_ssh_adjust +++ b/compass/ocean/tests/isomip_plus/namelist.single_layer.forward_and_ssh_adjust @@ -1,5 +1,6 @@ config_time_integrator='RK4' config_dt='00:00:10' +config_use_implicit_top_drag = .false. config_bottom_drag_mode = 'explicit' config_explicit_bottom_drag_coeff = 3.0e-3 config_disable_thick_vadv = .true. diff --git a/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust b/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust index e70be68f35..8c3be8b797 100644 --- a/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust +++ b/compass/ocean/tests/isomip_plus/namelist.thin_film.forward_and_ssh_adjust @@ -11,4 +11,3 @@ config_zero_drying_velocity_ramp_hmax=1e-1 config_use_ssh_gradient_wetting_drying=.false. config_prevent_drying=.true. config_zero_drying_velocity = .true. -config_use_implicit_top_drag = .true. From 2a116babf44c696ed268ee03a9b289e2a679bb98 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Sun, 16 Mar 2025 16:37:40 -0500 Subject: [PATCH 26/30] fixup viz --- .../ocean/tests/isomip_plus/viz/__init__.py | 216 +++++++++--------- 1 file changed, 110 insertions(+), 106 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/viz/__init__.py b/compass/ocean/tests/isomip_plus/viz/__init__.py index 2748598493..722e672c47 100644 --- a/compass/ocean/tests/isomip_plus/viz/__init__.py +++ b/compass/ocean/tests/isomip_plus/viz/__init__.py @@ -71,7 +71,6 @@ def run(self): out_dir = '.' dsMesh = xarray.open_dataset(f'{sim_dir}/init.nc') - maxBottomDepth = dsMesh.bottomDepth.max().values dsOut = xarray.open_dataset(f'{sim_dir}/output.nc') dsIce = xarray.open_dataset(f'{sim_dir}/land_ice_fluxes.nc') if dsOut.sizes['Time'] < 50: @@ -95,25 +94,6 @@ def run(self): cmap='cmo.balance', vmin=-5e-1, vmax=5e-1) - if dsOut.sizes['Time'] > 24: - days = 30 - samples_per_day = 1 - dsAve = dsOut.isel(Time=slice(-int(days * samples_per_day), -1)) - velocityMagSq = (dsAve.velocityX.isel(nVertLevels=0)**2. + - dsAve.velocityX.isel(nVertLevels=0)**2.) - rmsVelocity = xarray.DataArray( - data=np.sqrt(velocityMagSq.mean('Time').values), - dims=['nCells']).expand_dims(dim='Time', axis=0) - plotter.plot_horiz_series(rmsVelocity, - 'rmsVelocity', 'rmsVelocity', - False, vmin=1.e-2, vmax=1.e0, - cmap='cmo.speed', cmap_scale='log') - - cavityMask = dsMesh.landIceFloatingMask == 1 - if 'xIsomipCell' in dsMesh.keys(): - xCell = dsMesh.xIsomipCell / 1.e3 - else: - xCell = dsMesh.xCell / 1.e3 # salinity = dsIce.landIceInterfaceSalinity salinity = dsOut.salinity.isel(nVertLevels=0) pressure = dsIce.landIcePressure @@ -131,61 +111,35 @@ def run(self): coeff_mushy * salinity / (1.0 - salinity / 1e3) dsIce['thermalForcing'] = temperature - ocn_freezing_temperature - dx = 2. - xbins = np.arange(xCell.isel(nCells=cavityMask).min(), - xCell.isel(nCells=cavityMask).max(), dx) if dsIce.sizes['Time'] > 14: dsIce = dsIce.isel(Time=slice(-13, -12)) - for time_slice in range(dsIce.sizes['Time']): - title = dsIce.xtime.isel(Time=time_slice).values - meltTransect = np.zeros_like(xbins) - landIceDraftTransect = np.zeros_like(xbins) - frictionVelocityTransect = np.zeros_like(xbins) - temperatureTransect = np.zeros_like(xbins) - thermalForcingTransect = np.zeros_like(xbins) - rmsVelocityTransect = np.zeros_like(xbins) - for i, xmin in enumerate(xbins): - binMask = np.logical_and(cavityMask, - np.logical_and(xCell >= xmin, - xCell < xmin + dx)) - dsTransect = dsIce.isel(nCells=binMask, Time=time_slice) - meltTransect[i] = dsTransect.landIceFreshwaterFlux.mean() - landIceDraftTransect[i] = dsTransect.landIceDraft.mean() - frictionVelocityTransect[i] = \ - dsTransect.landIceFrictionVelocity.mean() - temperatureTransect[i] = \ - dsTransect.landIceInterfaceTemperature.mean() - thermalForcingTransect[i] = dsTransect.thermalForcing.mean() - if dsOut.sizes['Time'] > 24: - rmsVelocityTransect[i] = np.mean(rmsVelocity[binMask]) + cavityMask = dsMesh.landIceFloatingMask == 1 + if 'xIsomipCell' in dsMesh.keys(): + xCell = dsMesh.xIsomipCell / 1.e3 + else: + xCell = dsMesh.xCell / 1.e3 - fig, ax = plt.subplots(4, 1, sharex=True, figsize=(4, 8)) - color = 'darkgreen' - x = xbins - xbins[0] - ax[0].set_title(title) - ax[0].plot(x, landIceDraftTransect, color=color) - ax[0].set_ylabel('Land ice draft (m)') - secPerYear = 365 * 24 * 60 * 60 - density = 1026. - ax[1].plot(x, meltTransect * secPerYear / density, color=color) - ax[1].set_ylabel('Melt rate (m/yr)') - ax[2].plot(x, frictionVelocityTransect, color=color) - ax[2].set_ylabel('Friction velocity (m/s)') - ax[3].plot(x, thermalForcingTransect, color=color) - ax[3].set_ylabel(r'ThermalForcing ($^{\circ}$C)') - ax[3].set_xlabel('Distance along ice flow (km)') - plt.tight_layout() - plt.savefig(f'meanTransect_tend{time_slice:03g}.png', dpi=600, - bbox_inches='tight', transparent=True) - plt.close(fig) + if dsOut.sizes['Time'] > 24: + days = 30 + samples_per_day = 1 + dsAve = dsOut.isel(Time=slice(-int(days * samples_per_day), -1)) + velocityMagSq = (dsAve.velocityX.isel(nVertLevels=0)**2. + + dsAve.velocityX.isel(nVertLevels=0)**2.) + rmsVelocity = xarray.DataArray( + data=np.sqrt(velocityMagSq.mean('Time').values), + dims=['nCells']).expand_dims(dim='Time', axis=0) + plotter.plot_horiz_series(rmsVelocity, + 'rmsVelocity', 'rmsVelocity', + False, vmin=1.e-2, vmax=1.e0, + cmap='cmo.speed', cmap_scale='log') + self._plot_transect(dsIce, rmsVelocity, xCell, cavityMask, + outFolder=f'{out_dir}/plots') + maxBottomDepth = dsMesh.bottomDepth.max().values plotter.plot_horiz_series(dsOut.ssh, 'ssh', 'ssh', True, vmin=-maxBottomDepth, vmax=0) if 'layerThickness' in dsOut.keys(): plotter.plot_layer_interfaces() - delssh = dsOut.ssh - dsMesh.ssh.isel(Time=0) - plotter.plot_horiz_series(delssh, 'ssh - ssh_init', 'delssh_init', - True, cmap='cmo.curl', vmin=-1, vmax=1) plotter.plot_horiz_series(dsOut.ssh + dsMesh.bottomDepth, 'H', 'H', True, vmin=min_column_thickness + 1e-10, vmax=500, cmap_set_under='k', @@ -204,6 +158,9 @@ def run(self): plotter.plot_horiz_series(delice, 'delLandIcePressure', 'delLandIcePressure', True, vmin=-1e6, vmax=1e6, cmap='cmo.balance') + delssh = dsOut.ssh - dsMesh.ssh.isel(Time=0) + plotter.plot_horiz_series(delssh, 'ssh - ssh_init', 'delssh_init', + True, cmap='cmo.curl', vmin=-1, vmax=1) # Cannot currently plot because it is on edges for t in range(dsOut.sizes['Time']): fig = plt.figure() @@ -224,6 +181,26 @@ def run(self): True, vmin=min_column_thickness + 1e-10, vmax=700, cmap_set_under='k', cmap_scale='log') + if os.path.exists('../simulation/land_ice_fluxes.nc'): + ds = xarray.open_mfdataset( + '../simulation/land_ice_fluxes.nc', + concat_dim='Time', combine='nested') + if ds.sizes['Time'] > 50 or \ + ds.sizes['Time'] != dsOut.sizes['Time']: + # assumes ds.sizes['Time'] > dsOut.sizes['Time'] + ds = ds.isel(Time=slice(-31, -1)) + dsOut = dsOut.isel(Time=slice(-31, -1)) + if ds.sizes['Time'] != dsOut.sizes['Time']: + print('error in time selection') + return + ds['layerThickness'] = dsOut.layerThickness + mPlotter = MoviePlotter(inFolder=sim_dir, + streamfunctionFolder=streamfunction_dir, + outFolder='{}/plots'.format(out_dir), + expt=expt, sectionY=section_y, + dsMesh=dsMesh, ds=ds, + showProgress=show_progress) + self._plot_land_ice_variables(ds, mPlotter, maxBottomDepth) ds = xarray.open_mfdataset( '{}/timeSeriesStatsMonthly*.nc'.format(sim_dir), @@ -270,45 +247,72 @@ def run(self): framesPerSecond=frames_per_second, extension=movie_format) - if os.path.exists('../simulation/land_ice_fluxes.nc'): - ds = xarray.open_mfdataset( - '../simulation/land_ice_fluxes.nc', - concat_dim='Time', combine='nested') - if ds.sizes['Time'] != dsOut.sizes['Time']: - # assumes ds.sizes['Time'] > dsOut.sizes['Time'] - ds = ds.isel(Time=slice(-31, -1)) - dsOut = dsOut.isel(Time=slice(-31, -1)) - if ds.sizes['Time'] != dsOut.sizes['Time']: - print('error in time selection') - return - if ds.sizes['Time'] > 50: - return - ds['layerThickness'] = dsOut.layerThickness - mPlotter = MoviePlotter(inFolder=sim_dir, - streamfunctionFolder=streamfunction_dir, - outFolder='{}/plots'.format(out_dir), - expt=expt, sectionY=section_y, - dsMesh=dsMesh, ds=ds, - showProgress=show_progress) - mPlotter.plot_horiz_series(ds.landIceFrictionVelocity, - 'landIceFrictionVelocity', - 'frictionVelocity', - True) - mPlotter.plot_horiz_series(ds.landIceFreshwaterFlux, - 'landIceFreshwaterFlux', 'melt', - True) - mPlotter.plot_horiz_series(ds.landIceFloatingFraction, - 'landIceFloatingFraction', - 'landIceFloatingFraction', - True, vmin=1e-16, vmax=1, - cmap_set_under='k') - mPlotter.plot_horiz_series(ds.landIceDraft, - 'landIceDraft', 'landIceDraft', - True, vmin=-maxBottomDepth, vmax=0.) - mPlotter.plot_melt_rates() - mPlotter.plot_ice_shelf_boundary_variables() - else: - print(f'didnt find {sim_dir}/land_ice_fluxes.nc') + def _plot_land_ice_variables(ds, mPlotter, maxBottomDepth): + mPlotter.plot_horiz_series(ds.landIceFrictionVelocity, + 'landIceFrictionVelocity', + 'frictionVelocity', + True) + mPlotter.plot_horiz_series(ds.landIceFreshwaterFlux, + 'landIceFreshwaterFlux', 'melt', + True) + mPlotter.plot_horiz_series(ds.landIceFloatingFraction, + 'landIceFloatingFraction', + 'landIceFloatingFraction', + True, vmin=1e-16, vmax=1, + cmap_set_under='k') + mPlotter.plot_horiz_series(ds.landIceDraft, + 'landIceDraft', 'landIceDraft', + True, vmin=-maxBottomDepth, vmax=0.) + mPlotter.plot_melt_rates() + mPlotter.plot_ice_shelf_boundary_variables() + + def _plot_transect(dsIce, rmsVelocity, xCell, cavityMask, outFolder='.'): + dx = 2. + xbins = np.arange(xCell.isel(nCells=cavityMask).min(), + xCell.isel(nCells=cavityMask).max(), dx) + for time_slice in range(dsIce.sizes['Time']): + title = dsIce.xtime.isel(Time=time_slice).values + meltTransect = np.zeros_like(xbins) + landIceDraftTransect = np.zeros_like(xbins) + frictionVelocityTransect = np.zeros_like(xbins) + temperatureTransect = np.zeros_like(xbins) + thermalForcingTransect = np.zeros_like(xbins) + rmsVelocityTransect = np.zeros_like(xbins) + for i, xmin in enumerate(xbins): + binMask = np.logical_and(cavityMask, + np.logical_and(xCell >= xmin, + xCell < xmin + dx)) + if (np.sum(binMask) < 1): + continue + dsTransect = dsIce.isel(nCells=binMask, Time=time_slice) + meltTransect[i] = dsTransect.landIceFreshwaterFlux.mean() + landIceDraftTransect[i] = dsTransect.landIceDraft.mean() + frictionVelocityTransect[i] = \ + dsTransect.landIceFrictionVelocity.mean() + temperatureTransect[i] = \ + dsTransect.landIceInterfaceTemperature.mean() + thermalForcingTransect[i] = dsTransect.thermalForcing.mean() + rmsVelocityTransect[i] = np.mean(rmsVelocity[binMask]) + + fig, ax = plt.subplots(4, 1, sharex=True, figsize=(4, 8)) + color = 'darkgreen' + x = xbins - xbins[0] + ax[0].set_title(title) + ax[0].plot(x, landIceDraftTransect, color=color) + ax[0].set_ylabel('Land ice draft (m)') + secPerYear = 365 * 24 * 60 * 60 + density = 1026. + ax[1].plot(x, meltTransect * secPerYear / density, color=color) + ax[1].set_ylabel('Melt rate (m/yr)') + ax[2].plot(x, frictionVelocityTransect, color=color) + ax[2].set_ylabel('Friction velocity (m/s)') + ax[3].plot(x, thermalForcingTransect, color=color) + ax[3].set_ylabel(r'ThermalForcing ($^{\circ}$C)') + ax[3].set_xlabel('Distance along ice flow (km)') + plt.tight_layout() + plt.savefig(f'{outFolder}/meanTransect_tend{time_slice:03g}.png', + dpi=600, bbox_inches='tight', transparent=True) + plt.close(fig) def file_complete(ds, fileName): From 5263bbbb06a45131d6b091a24d62b30c643b2a86 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Sun, 16 Mar 2025 17:41:32 -0500 Subject: [PATCH 27/30] Tests are only single layer if vert_levels changed in cfg --- compass/ocean/tests/isomip_plus/__init__.py | 19 ++--------------- .../ocean/tests/isomip_plus/initial_state.py | 21 +++++++------------ ...melist.single_layer.forward_and_ssh_adjust | 12 ----------- 3 files changed, 10 insertions(+), 42 deletions(-) delete mode 100644 compass/ocean/tests/isomip_plus/namelist.single_layer.forward_and_ssh_adjust diff --git a/compass/ocean/tests/isomip_plus/__init__.py b/compass/ocean/tests/isomip_plus/__init__.py index 763d653283..00102aa68e 100644 --- a/compass/ocean/tests/isomip_plus/__init__.py +++ b/compass/ocean/tests/isomip_plus/__init__.py @@ -36,7 +36,7 @@ def __init__(self, mpas_core): experiment=experiment, vertical_coordinate=vertical_coordinate, time_varying_forcing=True)) - for vertical_coordinate in ['z-star', 'sigma', 'single_layer']: + for vertical_coordinate in ['z-star', 'sigma']: self.add_test_case( IsomipPlusTest( test_group=self, resolution=resolution, @@ -50,7 +50,7 @@ def __init__(self, mpas_core): vertical_coordinate=vertical_coordinate, time_varying_forcing=True, thin_film_present=True)) - for vertical_coordinate in ['z-star', 'sigma', 'single_layer']: + for vertical_coordinate in ['z-star', 'sigma']: self.add_test_case( IsomipPlusTest( test_group=self, resolution=resolution, @@ -75,21 +75,6 @@ def __init__(self, mpas_core): vertical_coordinate=vertical_coordinate, tidal_forcing=True, thin_film_present=True)) - for vertical_coordinate in ['single_layer']: - self.add_test_case( - IsomipPlusTest( - test_group=self, resolution=resolution, - experiment=experiment, - vertical_coordinate=vertical_coordinate, - tidal_forcing=True, - thin_film_present=False)) - self.add_test_case( - IsomipPlusTest( - test_group=self, resolution=resolution, - experiment=experiment, - vertical_coordinate=vertical_coordinate, - tidal_forcing=True, - thin_film_present=True)) for resolution in [2.]: for experiment in ['Ocean0', 'Ocean1', 'Ocean2']: diff --git a/compass/ocean/tests/isomip_plus/initial_state.py b/compass/ocean/tests/isomip_plus/initial_state.py index aacc25d82b..c58c36a284 100644 --- a/compass/ocean/tests/isomip_plus/initial_state.py +++ b/compass/ocean/tests/isomip_plus/initial_state.py @@ -202,19 +202,14 @@ def _compute_initial_condition(self): init_top_sal = section.getfloat('init_top_sal') init_bot_sal = section.getfloat('init_bot_sal') - if self.vertical_coordinate == 'single_layer': - # Initialize constant T,S - ds['temperature'] = init_bot_temp * xr.ones_like(ds.zmid) - ds['salinity'] = init_bot_sal * xr.ones_like(ds.zmid) - else: - # Initialize T,S as linear functions with max depth - max_bottom_depth = -config.getfloat('vertical_grid', - 'bottom_depth') - frac = (0. - ds.zMid) / (0. - max_bottom_depth) - ds['temperature'] = \ - (1.0 - frac) * init_top_temp + frac * init_bot_temp - ds['salinity'] = \ - (1.0 - frac) * init_top_sal + frac * init_bot_sal + # Initialize T,S as linear functions with max depth + max_bottom_depth = -config.getfloat('vertical_grid', + 'bottom_depth') + frac = (0. - ds.zMid) / (0. - max_bottom_depth) + ds['temperature'] = \ + (1.0 - frac) * init_top_temp + frac * init_bot_temp + ds['salinity'] = \ + (1.0 - frac) * init_top_sal + frac * init_bot_sal if thin_film_present: # for thin film cells, set temperature to freezing point diff --git a/compass/ocean/tests/isomip_plus/namelist.single_layer.forward_and_ssh_adjust b/compass/ocean/tests/isomip_plus/namelist.single_layer.forward_and_ssh_adjust deleted file mode 100644 index d1eeee1657..0000000000 --- a/compass/ocean/tests/isomip_plus/namelist.single_layer.forward_and_ssh_adjust +++ /dev/null @@ -1,12 +0,0 @@ -config_time_integrator='RK4' -config_dt='00:00:10' -config_use_implicit_top_drag = .false. -config_bottom_drag_mode = 'explicit' -config_explicit_bottom_drag_coeff = 3.0e-3 -config_disable_thick_vadv = .true. -config_disable_thick_sflux = .true. -config_disable_vel_hmix = .true. -config_disable_vel_vmix = .true. -config_disable_vel_vadv = .true. -config_disable_tr_all_tend = .true. -config_pressure_gradient_type = 'ssh_gradient' From ea7b54dd6ff4d9747e658aa35df1ba3dd3e13125 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 17 Mar 2025 13:10:11 -0500 Subject: [PATCH 28/30] Move performance viz to viz step --- compass/ocean/tests/isomip_plus/forward.py | 88 ------------------- .../ocean/tests/isomip_plus/isomip_plus.cfg | 3 + .../ocean/tests/isomip_plus/viz/__init__.py | 73 +++++++++++++++ 3 files changed, 76 insertions(+), 88 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/forward.py b/compass/ocean/tests/isomip_plus/forward.py index 57bad7e632..5ca32ed052 100644 --- a/compass/ocean/tests/isomip_plus/forward.py +++ b/compass/ocean/tests/isomip_plus/forward.py @@ -1,12 +1,5 @@ -import os -import shutil - -import numpy as np -import xarray - from compass.model import run_model from compass.ocean.tests.isomip_plus.evap import update_evaporation_flux -from compass.ocean.tests.isomip_plus.viz.plot import MoviePlotter from compass.ocean.time import get_time_interval_string from compass.step import Step @@ -157,87 +150,6 @@ def run(self): run_model(self) - if self.name == 'performance': - # plot a few fields - plot_folder = f'{self.work_dir}/plots' - if os.path.exists(plot_folder): - shutil.rmtree(plot_folder) - min_column_thickness = self.config.getfloat('isomip_plus', - 'min_column_thickness') - - dsMesh = xarray.open_dataset(os.path.join(self.work_dir, - 'init.nc')) - ds = xarray.open_dataset(os.path.join(self.work_dir, 'output.nc')) - ds['landIceDraft'] = \ - dsMesh.landIceDraft.expand_dims('Time', axis=0) - section_y = self.config.getfloat('isomip_plus_viz', 'section_y') - # show progress only if we're not writing to a log file - show_progress = self.log_filename is None - plotter = MoviePlotter(inFolder=self.work_dir, - streamfunctionFolder=self.work_dir, - outFolder=plot_folder, sectionY=section_y, - dsMesh=dsMesh, ds=ds, expt=self.experiment, - showProgress=show_progress) - - bottomDepth = ds.bottomDepth.expand_dims(dim='Time', axis=0) - plotter.plot_horiz_series(ds.ssh + bottomDepth, - 'H', 'H', True, - vmin=min_column_thickness, vmax=700, - cmap_set_under='r', cmap_scale='log') - 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') - - 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')) - if 'topDragMagnitude' in dsIce.keys(): - plotter.plot_horiz_series( - dsIce.topDragMagnitude, - 'topDragMagnitude', 'topDragMagnitude', True, - vmin=0 + tol, vmax=np.max(dsIce.topDragMagnitude.values), - cmap_set_under='k') - if 'landIceHeatFlux' in dsIce.keys(): - plotter.plot_horiz_series( - dsIce.landIceHeatFlux, - 'landIceHeatFlux', 'landIceHeatFlux', True, - vmin=np.min(dsIce.landIceHeatFlux.values), - vmax=np.max(dsIce.landIceHeatFlux.values)) - if 'landIceInterfaceTemperature' in dsIce.keys(): - plotter.plot_horiz_series( - dsIce.landIceInterfaceTemperature, - 'landIceInterfaceTemperature', - 'landIceInterfaceTemperature', - True, - vmin=np.min(dsIce.landIceInterfaceTemperature.values), - vmax=np.max(dsIce.landIceInterfaceTemperature.values)) - if 'landIceFreshwaterFlux' in dsIce.keys(): - plotter.plot_horiz_series( - dsIce.landIceFreshwaterFlux, - 'landIceFreshwaterFlux', 'landIceFreshwaterFlux', True, - vmin=0 + tol, vmax=1e-4, - cmap_set_under='k', cmap_scale='log') - if 'landIceFraction' in dsIce.keys(): - 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', diff --git a/compass/ocean/tests/isomip_plus/isomip_plus.cfg b/compass/ocean/tests/isomip_plus/isomip_plus.cfg index bcaf46fe43..4c5c827bdb 100644 --- a/compass/ocean/tests/isomip_plus/isomip_plus.cfg +++ b/compass/ocean/tests/isomip_plus/isomip_plus.cfg @@ -144,6 +144,9 @@ plot_haney = True # whether to plot the barotropic and overturning streamfunctions plot_streamfunctions = True +# whether to plot output from the performance step, usually for debugging +plot_performance_fields = True + # frames per second for movies frames_per_second = 30 diff --git a/compass/ocean/tests/isomip_plus/viz/__init__.py b/compass/ocean/tests/isomip_plus/viz/__init__.py index 722e672c47..f1ee02fb85 100644 --- a/compass/ocean/tests/isomip_plus/viz/__init__.py +++ b/compass/ocean/tests/isomip_plus/viz/__init__.py @@ -53,6 +53,7 @@ def run(self): config = self.config section = config['isomip_plus_viz'] plot_streamfunctions = section.getboolean('plot_streamfunctions') + plot_performance_fields = section.getboolean('plot_performance_fields') plot_haney = section.getboolean('plot_haney') frames_per_second = section.getint('frames_per_second') # output_slices = section.getint('output_slices') @@ -247,7 +248,40 @@ def run(self): framesPerSecond=frames_per_second, extension=movie_format) + if plot_performance_fields: + plot_folder = f'{out_dir}/performance_plots' + min_column_thickness = self.config.getfloat('isomip_plus', + 'min_column_thickness') + + ds = xarray.open_dataset('../performance/output.nc') + # show progress only if we're not writing to a log file + show_progress = self.log_filename is None + plotter = MoviePlotter(inFolder=self.work_dir, + streamfunctionFolder=self.work_dir, + outFolder=plot_folder, sectionY=section_y, + dsMesh=dsMesh, ds=ds, expt=self.experiment, + showProgress=show_progress) + + bottomDepth = ds.bottomDepth.expand_dims(dim='Time', axis=0) + plotter.plot_horiz_series(ds.ssh.isel(Time=-1) + bottomDepth, + 'H', 'H', True, + vmin=min_column_thickness, vmax=700, + cmap_set_under='r', cmap_scale='log') + 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') + + plotter.plot_3d_field_top_bot_section( + ds.salinity, nameInTitle='salinity', prefix='salin', + units='PSU', vmin=33.8, vmax=34.7, cmap='cmo.haline') + + dsIce = xarray.open_dataset('../performance/land_ice_fluxes.nc') + self._plot_land_ice_variables(dsIce, plotter, maxBottomDepth) + def _plot_land_ice_variables(ds, mPlotter, maxBottomDepth): + tol = 1e-10 mPlotter.plot_horiz_series(ds.landIceFrictionVelocity, 'landIceFrictionVelocity', 'frictionVelocity', @@ -263,6 +297,45 @@ def _plot_land_ice_variables(ds, mPlotter, maxBottomDepth): mPlotter.plot_horiz_series(ds.landIceDraft, 'landIceDraft', 'landIceDraft', True, vmin=-maxBottomDepth, vmax=0.) + if 'topDragMagnitude' in ds.keys(): + mPlotter.plot_horiz_series( + ds.topDragMagnitude, + 'topDragMagnitude', 'topDragMagnitude', True, + vmin=0 + tol, vmax=np.max(ds.topDragMagnitude.values), + cmap_set_under='k') + if 'landIceHeatFlux' in ds.keys(): + mPlotter.plot_horiz_series( + ds.landIceHeatFlux, + 'landIceHeatFlux', 'landIceHeatFlux', True, + vmin=np.min(ds.landIceHeatFlux.values), + vmax=np.max(ds.landIceHeatFlux.values)) + if 'landIceInterfaceTemperature' in ds.keys(): + mPlotter.plot_horiz_series( + ds.landIceInterfaceTemperature, + 'landIceInterfaceTemperature', + 'landIceInterfaceTemperature', + True, + vmin=np.min(ds.landIceInterfaceTemperature.values), + vmax=np.max(ds.landIceInterfaceTemperature.values)) + if 'landIceFreshwaterFlux' in ds.keys(): + mPlotter.plot_horiz_series( + ds.landIceFreshwaterFlux, + 'landIceFreshwaterFlux', 'landIceFreshwaterFlux', True, + vmin=0 + tol, vmax=1e-4, + cmap_set_under='k', cmap_scale='log') + if 'landIceFraction' in ds.keys(): + mPlotter.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') + if 'landIceFloatingFraction' in ds.keys(): + mPlotter.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') mPlotter.plot_melt_rates() mPlotter.plot_ice_shelf_boundary_variables() From 3d730f83bc326e64fb8bd464c6fd4fd44a2ed4f6 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 17 Mar 2025 15:07:52 -0500 Subject: [PATCH 29/30] Fix transect-binned plot --- .../ocean/tests/isomip_plus/viz/__init__.py | 234 +++++++++--------- compass/ocean/tests/isomip_plus/viz/plot.py | 1 - 2 files changed, 120 insertions(+), 115 deletions(-) diff --git a/compass/ocean/tests/isomip_plus/viz/__init__.py b/compass/ocean/tests/isomip_plus/viz/__init__.py index f1ee02fb85..d9adf24a6e 100644 --- a/compass/ocean/tests/isomip_plus/viz/__init__.py +++ b/compass/ocean/tests/isomip_plus/viz/__init__.py @@ -95,10 +95,11 @@ def run(self): cmap='cmo.balance', vmin=-5e-1, vmax=5e-1) - # salinity = dsIce.landIceInterfaceSalinity + # Compute thermal forcing for whole water column + # landIceInterfaceSalinity and landIceInterfaceTemperature could be + # used for thermal forcing at the interface salinity = dsOut.salinity.isel(nVertLevels=0) pressure = dsIce.landIcePressure - # temperature = dsIce.landIceInterfaceTemperature temperature = dsOut.temperature.isel(nVertLevels=0) coeff_0 = 6.22e-2 coeff_S = -5.63e-2 @@ -114,7 +115,6 @@ def run(self): if dsIce.sizes['Time'] > 14: dsIce = dsIce.isel(Time=slice(-13, -12)) - cavityMask = dsMesh.landIceFloatingMask == 1 if 'xIsomipCell' in dsMesh.keys(): xCell = dsMesh.xIsomipCell / 1.e3 else: @@ -133,8 +133,11 @@ def run(self): 'rmsVelocity', 'rmsVelocity', False, vmin=1.e-2, vmax=1.e0, cmap='cmo.speed', cmap_scale='log') - self._plot_transect(dsIce, rmsVelocity, xCell, cavityMask, - outFolder=f'{out_dir}/plots') + cavityMask = np.where(dsMesh.landIceFloatingMask == 1)[0] + dsCavity = dsIce.isel(nCells=cavityMask) + rmsVelocityCavity = rmsVelocity.isel(nCells=cavityMask) + xMasked = xCell.isel(nCells=cavityMask) + _plot_transect(dsCavity, rmsVelocityCavity, xMasked) maxBottomDepth = dsMesh.bottomDepth.max().values plotter.plot_horiz_series(dsOut.ssh, 'ssh', 'ssh', True, @@ -201,7 +204,7 @@ def run(self): expt=expt, sectionY=section_y, dsMesh=dsMesh, ds=ds, showProgress=show_progress) - self._plot_land_ice_variables(ds, mPlotter, maxBottomDepth) + _plot_land_ice_variables(ds, mPlotter, maxBottomDepth) ds = xarray.open_mfdataset( '{}/timeSeriesStatsMonthly*.nc'.format(sim_dir), @@ -278,114 +281,117 @@ def run(self): units='PSU', vmin=33.8, vmax=34.7, cmap='cmo.haline') dsIce = xarray.open_dataset('../performance/land_ice_fluxes.nc') - self._plot_land_ice_variables(dsIce, plotter, maxBottomDepth) - - def _plot_land_ice_variables(ds, mPlotter, maxBottomDepth): - tol = 1e-10 - mPlotter.plot_horiz_series(ds.landIceFrictionVelocity, - 'landIceFrictionVelocity', - 'frictionVelocity', - True) - mPlotter.plot_horiz_series(ds.landIceFreshwaterFlux, - 'landIceFreshwaterFlux', 'melt', - True) - mPlotter.plot_horiz_series(ds.landIceFloatingFraction, - 'landIceFloatingFraction', - 'landIceFloatingFraction', - True, vmin=1e-16, vmax=1, - cmap_set_under='k') - mPlotter.plot_horiz_series(ds.landIceDraft, - 'landIceDraft', 'landIceDraft', - True, vmin=-maxBottomDepth, vmax=0.) - if 'topDragMagnitude' in ds.keys(): - mPlotter.plot_horiz_series( - ds.topDragMagnitude, - 'topDragMagnitude', 'topDragMagnitude', True, - vmin=0 + tol, vmax=np.max(ds.topDragMagnitude.values), - cmap_set_under='k') - if 'landIceHeatFlux' in ds.keys(): - mPlotter.plot_horiz_series( - ds.landIceHeatFlux, - 'landIceHeatFlux', 'landIceHeatFlux', True, - vmin=np.min(ds.landIceHeatFlux.values), - vmax=np.max(ds.landIceHeatFlux.values)) - if 'landIceInterfaceTemperature' in ds.keys(): - mPlotter.plot_horiz_series( - ds.landIceInterfaceTemperature, - 'landIceInterfaceTemperature', - 'landIceInterfaceTemperature', - True, - vmin=np.min(ds.landIceInterfaceTemperature.values), - vmax=np.max(ds.landIceInterfaceTemperature.values)) - if 'landIceFreshwaterFlux' in ds.keys(): - mPlotter.plot_horiz_series( - ds.landIceFreshwaterFlux, - 'landIceFreshwaterFlux', 'landIceFreshwaterFlux', True, - vmin=0 + tol, vmax=1e-4, - cmap_set_under='k', cmap_scale='log') - if 'landIceFraction' in ds.keys(): - mPlotter.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') - if 'landIceFloatingFraction' in ds.keys(): - mPlotter.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') - mPlotter.plot_melt_rates() - mPlotter.plot_ice_shelf_boundary_variables() - - def _plot_transect(dsIce, rmsVelocity, xCell, cavityMask, outFolder='.'): - dx = 2. - xbins = np.arange(xCell.isel(nCells=cavityMask).min(), - xCell.isel(nCells=cavityMask).max(), dx) - for time_slice in range(dsIce.sizes['Time']): - title = dsIce.xtime.isel(Time=time_slice).values - meltTransect = np.zeros_like(xbins) - landIceDraftTransect = np.zeros_like(xbins) - frictionVelocityTransect = np.zeros_like(xbins) - temperatureTransect = np.zeros_like(xbins) - thermalForcingTransect = np.zeros_like(xbins) - rmsVelocityTransect = np.zeros_like(xbins) - for i, xmin in enumerate(xbins): - binMask = np.logical_and(cavityMask, - np.logical_and(xCell >= xmin, - xCell < xmin + dx)) - if (np.sum(binMask) < 1): - continue - dsTransect = dsIce.isel(nCells=binMask, Time=time_slice) - meltTransect[i] = dsTransect.landIceFreshwaterFlux.mean() - landIceDraftTransect[i] = dsTransect.landIceDraft.mean() - frictionVelocityTransect[i] = \ - dsTransect.landIceFrictionVelocity.mean() - temperatureTransect[i] = \ - dsTransect.landIceInterfaceTemperature.mean() - thermalForcingTransect[i] = dsTransect.thermalForcing.mean() - rmsVelocityTransect[i] = np.mean(rmsVelocity[binMask]) - - fig, ax = plt.subplots(4, 1, sharex=True, figsize=(4, 8)) - color = 'darkgreen' - x = xbins - xbins[0] - ax[0].set_title(title) - ax[0].plot(x, landIceDraftTransect, color=color) - ax[0].set_ylabel('Land ice draft (m)') - secPerYear = 365 * 24 * 60 * 60 - density = 1026. - ax[1].plot(x, meltTransect * secPerYear / density, color=color) - ax[1].set_ylabel('Melt rate (m/yr)') - ax[2].plot(x, frictionVelocityTransect, color=color) - ax[2].set_ylabel('Friction velocity (m/s)') - ax[3].plot(x, thermalForcingTransect, color=color) - ax[3].set_ylabel(r'ThermalForcing ($^{\circ}$C)') - ax[3].set_xlabel('Distance along ice flow (km)') - plt.tight_layout() - plt.savefig(f'{outFolder}/meanTransect_tend{time_slice:03g}.png', - dpi=600, bbox_inches='tight', transparent=True) - plt.close(fig) + _plot_land_ice_variables(dsIce, plotter, maxBottomDepth) + + +def _plot_land_ice_variables(ds, mPlotter, maxBottomDepth): + tol = 1e-10 + mPlotter.plot_horiz_series(ds.landIceFrictionVelocity, + 'landIceFrictionVelocity', + 'frictionVelocity', + True) + mPlotter.plot_horiz_series(ds.landIceFreshwaterFlux, + 'landIceFreshwaterFlux', 'melt', + True) + mPlotter.plot_horiz_series(ds.landIceFloatingFraction, + 'landIceFloatingFraction', + 'landIceFloatingFraction', + True, vmin=1e-16, vmax=1, + cmap_set_under='k') + mPlotter.plot_horiz_series(ds.landIceDraft, + 'landIceDraft', 'landIceDraft', + True, vmin=-maxBottomDepth, vmax=0.) + if 'topDragMagnitude' in ds.keys(): + mPlotter.plot_horiz_series( + ds.topDragMagnitude, + 'topDragMagnitude', 'topDragMagnitude', True, + vmin=0 + tol, vmax=np.max(ds.topDragMagnitude.values), + cmap_set_under='k') + if 'landIceHeatFlux' in ds.keys(): + mPlotter.plot_horiz_series( + ds.landIceHeatFlux, + 'landIceHeatFlux', 'landIceHeatFlux', True, + vmin=np.min(ds.landIceHeatFlux.values), + vmax=np.max(ds.landIceHeatFlux.values)) + if 'landIceInterfaceTemperature' in ds.keys(): + mPlotter.plot_horiz_series( + ds.landIceInterfaceTemperature, + 'landIceInterfaceTemperature', + 'landIceInterfaceTemperature', + True, + vmin=np.min(ds.landIceInterfaceTemperature.values), + vmax=np.max(ds.landIceInterfaceTemperature.values)) + if 'landIceFreshwaterFlux' in ds.keys(): + mPlotter.plot_horiz_series( + ds.landIceFreshwaterFlux, + 'landIceFreshwaterFlux', 'landIceFreshwaterFlux', True, + vmin=0 + tol, vmax=1e-4, + cmap_set_under='k', cmap_scale='log') + if 'landIceFraction' in ds.keys(): + mPlotter.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') + if 'landIceFloatingFraction' in ds.keys(): + mPlotter.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') + mPlotter.plot_melt_rates() + mPlotter.plot_ice_shelf_boundary_variables() + + +def _plot_transect(dsIce, rmsVelocity, xCell, outFolder='plots'): + dx = 2. + xmin = xCell.min() + xmax = xCell.max() + xbins = np.arange(xmin, xmax, dx) + for time_slice in range(dsIce.sizes['Time']): + title = dsIce.xtime.isel(Time=time_slice).values + meltTransect = np.zeros_like(xbins) + landIceDraftTransect = np.zeros_like(xbins) + frictionVelocityTransect = np.zeros_like(xbins) + temperatureTransect = np.zeros_like(xbins) + thermalForcingTransect = np.zeros_like(xbins) + rmsVelocityTransect = np.zeros_like(xbins) + for i, xmin in enumerate(xbins): + binMask = np.where(np.logical_and(xCell >= xmin, + xCell < xmin + dx))[0] + if (np.sum(binMask) < 1): + continue + dsTransect = dsIce.isel(nCells=binMask, Time=time_slice) + meltTransect[i] = dsTransect.landIceFreshwaterFlux.mean() + landIceDraftTransect[i] = dsTransect.landIceDraft.mean() + frictionVelocityTransect[i] = \ + dsTransect.landIceFrictionVelocity.mean() + temperatureTransect[i] = \ + dsTransect.landIceInterfaceTemperature.mean() + thermalForcingTransect[i] = dsTransect.thermalForcing.mean() + rmsVelocityTransect[i] = \ + np.mean(rmsVelocity.isel(Time=time_slice, nCells=binMask)) + + fig, ax = plt.subplots(4, 1, sharex=True, figsize=(4, 8)) + color = 'darkgreen' + x = xbins - xbins[0] + ax[0].set_title(title) + ax[0].plot(x, landIceDraftTransect, color=color) + ax[0].set_ylabel('Land ice draft (m)') + secPerYear = 365 * 24 * 60 * 60 + density = 1026. + ax[1].plot(x, meltTransect * secPerYear / density, color=color) + ax[1].set_ylabel('Melt rate (m/yr)') + ax[2].plot(x, frictionVelocityTransect, color=color) + ax[2].set_ylabel('Friction velocity (m/s)') + ax[3].plot(x, thermalForcingTransect, color=color) + ax[3].set_ylabel(r'ThermalForcing ($^{\circ}$C)') + ax[3].set_xlabel('Distance along ice flow (km)') + plt.tight_layout() + plt.savefig(f'{outFolder}/meanTransect_tend{time_slice:03g}.png', + dpi=600, bbox_inches='tight', transparent=True) + plt.close(fig) def file_complete(ds, fileName): diff --git a/compass/ocean/tests/isomip_plus/viz/plot.py b/compass/ocean/tests/isomip_plus/viz/plot.py index bc67d4b449..876fcd4d38 100644 --- a/compass/ocean/tests/isomip_plus/viz/plot.py +++ b/compass/ocean/tests/isomip_plus/viz/plot.py @@ -107,7 +107,6 @@ def plot_melt_time_series(self, sshMax=None, wctMax=None): meanMeltRate = totalMeltFlux / totalArea / rho_fw * secPerYear self.plot_time_series(meanMeltRate, 'mean melt rate', 'meanMeltRate', 'm/yr') - print(meanMeltRate.values) self.plot_time_series(1e-6 * totalMeltFlux, 'total melt flux', f'totalMeltFlux{suffix}', 'kT/yr') From 602d98e4243e4e0fcc85c32dc7721a739bc0c5c3 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 17 Mar 2025 15:20:18 -0500 Subject: [PATCH 30/30] fixup compute_land_ice_draft_from_pressure --- compass/ocean/iceshelf.py | 6 ++---- compass/ocean/tests/isomip_plus/initial_state.py | 3 ++- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/compass/ocean/iceshelf.py b/compass/ocean/iceshelf.py index 137361f1b1..f34c03298d 100644 --- a/compass/ocean/iceshelf.py +++ b/compass/ocean/iceshelf.py @@ -96,10 +96,8 @@ def compute_land_ice_draft_from_pressure(land_ice_pressure, modify_mask, if ref_density is None: ref_density = constants['SHR_CONST_RHOSW'] land_ice_draft_array = \ - - (modify_mask.values * - land_ice_pressure.values / (ref_density * gravity)) - land_ice_draft = xarray.DataArray(data=land_ice_draft_array, - dims=(land_ice_pressure.dims)) + modify_mask * -land_ice_pressure / (ref_density * gravity) + land_ice_draft = land_ice_draft_array.transpose(*land_ice_pressure.dims) return land_ice_draft diff --git a/compass/ocean/tests/isomip_plus/initial_state.py b/compass/ocean/tests/isomip_plus/initial_state.py index c58c36a284..7199f89e08 100644 --- a/compass/ocean/tests/isomip_plus/initial_state.py +++ b/compass/ocean/tests/isomip_plus/initial_state.py @@ -158,7 +158,8 @@ def _compute_initial_condition(self): modify_mask = ds.bottomDepth > 0. land_ice_draft = compute_land_ice_draft_from_pressure( land_ice_pressure=land_ice_pressure, - modify_mask=modify_mask) + modify_mask=modify_mask, + ref_density=1028.) land_ice_draft = np.maximum(land_ice_draft, -ds.bottomDepth) ds['ssh'] = land_ice_draft