Skip to content

Commit

Permalink
Fix computation of barotropic streamfunction
Browse files Browse the repository at this point in the history 
This merge includes several fixes:
* Most important, the computation of the BSF only includes one
  constraint on the mean BSF at boundary vertices above 45 deg. S.
  Previously, the computation attempted to make the BSF zero on
  all boundary vertices, which is not correct since the value on
  isolated boundaries (e.g. the coast of South America vs.
  Antarctica) should differ by the transport between them.
* The full velocity (including bolus and submesoscale components)
  is used instead of just the resolved.
* The mask for transport has been fixed to correctly identify
  the deepest level on an edge.

The colormap has been changed (to blue-orange-div) to indicate a
clear break with the incorrect BSF plots before this fix and to
slightly improve the contrast at high absolute values of the BSF.

The code has been linted to follow PEP8 conventions
  • Loading branch information
@xylar
xylar committed Feb 14, 2025
1 parent 4ad6733 commit a253020
Show file tree
Hide file tree
Showing 2 changed files with 106 additions and 40 deletions.
8 changes: 6 additions & 2 deletions mpas_analysis/default.cfg
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1544,15 +1544,15 @@ comparisonGrids = ['latlon', 'subpolar_north_atlantic']
## (if available)

# colormap for model/observations
colormapNameResult = cmo.curl
colormapNameResult = blue-orange-div
# whether the colormap is indexed or continuous
colormapTypeResult = continuous
# color indices into colormapName for filled contours
# the type of norm used in the colormap
normTypeResult = symLog
# A dictionary with keywords for the norm
normArgsResult = {'linthresh': 30., 'linscale': 0.5, 'vmin': -100., 'vmax': 100.}
colorbarTicksResult = [-100.,-40.,-20.,-10., 0., 10., 20., 40.,100.]
colorbarTicksResult = [-100.,-40., -20., -10., 0., 10., 20., 40., 100.]
# Adding contour lines to the figure
contourLevelsResult = np.arange(-100., 101.0, 10.)
contourThicknessResult = 0.5
Expand Down Expand Up @@ -1580,6 +1580,10 @@ comparisonGrids = ['latlon', 'subpolar_north_atlantic']
# list of tuples(pairs) of depths (min, max) to integrate horizontal transport over
depthRanges = [(0.0, -10000.0), (0.0, -2000.0)]

# minimum latitude (degrees) above which the mean BSF on boundary vertices
# averages to zero
minLatitudeForZeroBSF = -45.0

[climatologyMapOHCAnomaly]
## options related to plotting horizontally remapped climatologies of
## ocean heat content (OHC) against control model results (if available)
Expand Down
138 changes: 100 additions & 38 deletions mpas_analysis/ocean/climatology_map_bsf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ class ClimatologyMapBSF(AnalysisTask):
----------
mpas_climatology_task : mpas_analysis.shared.climatology.MpasClimatologyTask
The task that produced the climatology to be remapped and plotted
"""
""" # noqa: E501

def __init__(self, config, mpas_climatology_task, control_config=None):
"""
Expand All @@ -44,7 +44,7 @@ def __init__(self, config, mpas_climatology_task, control_config=None):
control_config : mpas_tools.config.MpasConfigParser, optional
Configuration options for a control run (if any)
"""
""" # noqa: E501
field_name = 'barotropicStreamfunction'
# call the constructor from the base class (AnalysisTask)
super().__init__(config=config, taskName='climatologyMapBSF',
Expand Down Expand Up @@ -74,23 +74,20 @@ def __init__(self, config, mpas_climatology_task, control_config=None):

mpas_field_name = field_name

variable_list = ['timeMonthly_avg_normalVelocity',
'timeMonthly_avg_layerThickness']
for min_depth, max_depth in depth_ranges:
depth_range_string = \
f'{np.abs(min_depth):g}-{np.abs(max_depth):g}m'
if np.abs(max_depth) < 6000.:
fname_title = f'Streamfunction over {depth_range_string}'
fname_clim = f'{field_name}_{depth_range_string}'
else:
fname_title = f'Barotropic Streamfunction'
fname_title = 'Barotropic Streamfunction'
fname_clim = field_name

remap_climatology_subtask = RemapMpasBSFClimatology(
mpas_climatology_task=mpas_climatology_task,
parent_task=self,
climatology_name=fname_clim,
variable_list=variable_list,
comparison_grid_names=comparison_grid_names,
seasons=seasons,
min_depth=min_depth,
Expand All @@ -113,7 +110,8 @@ def __init__(self, config, mpas_climatology_task, control_config=None):
for comparison_grid_name in comparison_grid_names:
for season in seasons:
# make a new subtask for this season and comparison grid
subtask_name = f'plot{season}_{comparison_grid_name}_{depth_range_string}'
subtask_name = f'plot{season}_{comparison_grid_name}' \
f'_{depth_range_string}'

subtask = PlotClimatologyMapSubtask(
self, season, comparison_grid_name,
Expand All @@ -140,10 +138,24 @@ class RemapMpasBSFClimatology(RemapMpasClimatologySubtask):
"""
A subtask for computing climatologies of the barotropic /subpolar gyre
streamfunction from climatologies of normal velocity and layer thickness
Attributes
----------
min_depth : float
The minimum depth for integration
max_depth : float
The maximum depth for integration
include_bolus : bool
Whether to include the bolus velocity in the streamfunction
include_submesoscale : bool
Whether to include the submesoscale velocity in the streamfunction
"""
def __init__(self, mpas_climatology_task, parent_task,
climatology_name, variable_list, seasons,
comparison_grid_names, min_depth, max_depth):
climatology_name, seasons, comparison_grid_names, min_depth,
max_depth):

"""
Construct the analysis task and adds it as a subtask of the
Expand All @@ -163,10 +175,6 @@ def __init__(self, mpas_climatology_task, parent_task,
the important field(s) in the climatology) used to name the
subdirectories for each stage of the climatology
variable_list : list of str
A list of variable names in ``timeSeriesStatsMonthly`` to be
included in the climatologies
seasons : list of str, optional
A list of seasons (keys in ``shared.constants.monthDictionary``)
to be computed or ['none'] (not ``None``) if only monthly
Expand All @@ -177,7 +185,10 @@ def __init__(self, mpas_climatology_task, parent_task,
min_depth, max_depth : float
The minimum and maximum depths for integration
"""
""" # noqa: E501

# we'll fill this in at setup time
variable_list = []

depth_range_string = f'{np.abs(min_depth):g}-{np.abs(max_depth):g}m'
subtask_name = f'remapMpasClimatology_{depth_range_string}'
Expand All @@ -190,6 +201,43 @@ def __init__(self, mpas_climatology_task, parent_task,

self.min_depth = min_depth
self.max_depth = max_depth
self.include_bolus = None
self.include_submesoscale = None

def setup_and_check(self):
"""
Perform steps to set up the analysis and check for errors in the setup.
"""
super().setup_and_check()

variable_list = ['timeMonthly_avg_normalVelocity',
'timeMonthly_avg_layerThickness']

# Add the bolus velocity if GM is enabled
try:
# the new name
self.include_bolus = self.namelist.getbool('config_use_gm')
except KeyError:
# the old name
self.include_bolus = self.namelist.getbool(
'config_use_standardgm')
try:
self.include_submesoscale = \
self.namelist.getbool('config_submesoscale_enable')
except KeyError:
# an old run without submesoscale
self.include_submesoscale = False

if self.include_bolus:
variable_list.append('timeMonthly_avg_normalGMBolusVelocity')

if self.include_submesoscale:
variable_list.append('timeMonthly_avg_normalMLEvelocity')

self.variableList = variable_list

# Add the variables and seasons, now that we have the variable list
self.mpasClimatologyTask.add_variables(self.variableList, self.seasons)

def customize_masked_climatology(self, climatology, season):
"""
Expand All @@ -215,7 +263,7 @@ def customize_masked_climatology(self, climatology, season):
ds_mesh = ds_mesh[['cellsOnEdge', 'cellsOnVertex', 'nEdgesOnCell',
'edgesOnCell', 'verticesOnCell', 'verticesOnEdge',
'dcEdge', 'dvEdge', 'bottomDepth', 'layerThickness',
'maxLevelCell']]
'maxLevelCell', 'latVertex']]
ds_mesh = ds_mesh.isel(Time=0)
ds_mesh.load()
bsf_vertex = self._compute_barotropic_streamfunction_vertex(
Expand Down Expand Up @@ -250,14 +298,20 @@ def _compute_transport(self, ds_mesh, ds):
ds_mesh.layerThickness)
z_mid_edge = 0.5*(z_mid.isel(nCells=cell0) +
z_mid.isel(nCells=cell1))
normal_velocity = \
ds.timeMonthly_avg_normalVelocity.isel(nEdges=inner_edges)

normal_velocity = ds.timeMonthly_avg_normalVelocity
if self.include_bolus:
normal_velocity += ds.timeMonthly_avg_normalGMBolusVelocity
if self.include_submesoscale:
normal_velocity += ds.timeMonthly_avg_normalMLEvelocity
normal_velocity = normal_velocity.isel(nEdges=inner_edges)

layer_thickness = ds.timeMonthly_avg_layerThickness
layer_thickness_edge = 0.5*(layer_thickness.isel(nCells=cell0) +
layer_thickness.isel(nCells=cell1))
mask_bottom = (vert_index < ds_mesh.maxLevelCell).T
mask_bottom_edge = 0.5*(mask_bottom.isel(nCells=cell0) +
mask_bottom.isel(nCells=cell1))
mask_bottom_edge = np.logical_and(mask_bottom.isel(nCells=cell0),
mask_bottom.isel(nCells=cell1))
masks = [mask_bottom_edge,
z_mid_edge <= self.min_depth,
z_mid_edge >= self.max_depth]
Expand All @@ -273,57 +327,65 @@ def _compute_barotropic_streamfunction_vertex(self, ds_mesh, ds):
inner_edges, transport = self._compute_transport(ds_mesh, ds)
self.logger.info('transport computed.')

config = self.config
min_lat = config.getfloat(
'climatologyMapBSF', 'minLatitudeForZeroBSF')

nvertices = ds_mesh.sizes['nVertices']

cells_on_vertex = ds_mesh.cellsOnVertex - 1
vertices_on_edge = ds_mesh.verticesOnEdge - 1
is_boundary_cov = cells_on_vertex == -1
boundary_vertices = np.logical_or(is_boundary_cov.isel(vertexDegree=0),
is_boundary_cov.isel(vertexDegree=1))
boundary_vertices = np.logical_or(boundary_vertices,
is_boundary_cov.isel(vertexDegree=2))
boundary_vertices = is_boundary_cov.sum(dim='vertexDegree') > 0

boundary_vertices = np.logical_and(
boundary_vertices,
ds_mesh.latVertex > np.deg2rad(min_lat)
)

# convert from boolean mask to indices
boundary_vertices = np.flatnonzero(boundary_vertices.values)

n_boundary_vertices = len(boundary_vertices)

n_inner_edges = len(inner_edges)

indices = np.zeros((2, 2*n_inner_edges+n_boundary_vertices), dtype=int)
data = np.zeros(2*n_inner_edges+n_boundary_vertices, dtype=float)
indices = np.zeros((2, 2 * n_inner_edges + n_boundary_vertices),
dtype=int)
data = np.zeros(2 * n_inner_edges + n_boundary_vertices, dtype=float)

# The difference between the streamfunction at vertices on an inner
# edge should be equal to the transport
v0 = vertices_on_edge.isel(nEdges=inner_edges, TWO=0).values
v1 = vertices_on_edge.isel(nEdges=inner_edges, TWO=1).values

ind = np.arange(n_inner_edges)
indices[0, 2*ind] = ind
indices[1, 2*ind] = v1
indices[0, 2 * ind] = ind
indices[1, 2 * ind] = v1
data[2*ind] = 1.

indices[0, 2*ind+1] = ind
indices[1, 2*ind+1] = v0
indices[0, 2 * ind + 1] = ind
indices[1, 2 * ind + 1] = v0
data[2*ind+1] = -1.

# the streamfunction should be zero at all boundary vertices
# Make the average of the stream function at boundary vertices north
# of min_lat equal to zero
ind = np.arange(n_boundary_vertices)
indices[0, 2*n_inner_edges + ind] = n_inner_edges + ind
indices[1, 2*n_inner_edges + ind] = boundary_vertices
data[2*n_inner_edges + ind] = 1.
indices[0, 2 * n_inner_edges + ind] = n_inner_edges
indices[1, 2 * n_inner_edges + ind] = ind
data[2 * n_inner_edges + ind] = 1.

rhs = np.zeros(n_inner_edges+n_boundary_vertices, dtype=float)
rhs = np.zeros(n_inner_edges + 1, dtype=float)

# convert to Sv
ind = np.arange(n_inner_edges)
rhs[ind] = 1e-6*transport
rhs[ind] = 1e-6 * transport

ind = np.arange(n_boundary_vertices)
rhs[n_inner_edges + ind] = 0.
rhs[n_inner_edges] = 0.

matrix = scipy.sparse.csr_matrix(
(data, indices),
shape=(n_inner_edges+n_boundary_vertices, nvertices))
shape=(n_inner_edges + 1, nvertices))

solution = scipy.sparse.linalg.lsqr(matrix, rhs)
bsf_vertex = xr.DataArray(-solution[0],
Expand Down

0 comments on commit a253020

Please sign in to comment.