Skip to content

Fix incorrect indexing for MOM grid u-cells #25

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.

Sign up for GitHub

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

Merged
dougiesquire merged 2 commits into from
Jun 19, 2025
Merged

Fix incorrect indexing for MOM grid u-cells #25

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
1d35bdf
mom_grid: Fix incorrect indexing for u-cells
dougiesquire Jun 17, 2025
fff1200
cice_grid: Write u[t]lat[lon]_bounds to CICE grid netcdf
dougiesquire Jun 17, 2025
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
51 changes: 51 additions & 0 deletions esmgrids/cice_grid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,8 @@ def fromfile(cls, h_grid_def, mask_file=None, description="CICE tripolar"):
if mask_file is not None:
with nc.Dataset(mask_file) as f:
mask_t = f.variables["kmt"][:]
else:
mask_t = None

return cls(
x_t=x_t,
Expand Down Expand Up @@ -78,6 +80,15 @@ def _create_2d_nc_var(self, f, name):
complevel=1,
)

def _create_3d_nc_var(self, f, name):
return f.createVariable(
name,
"f8",
dimensions=("nj", "ni", "nvertices"),
compression="zlib",
complevel=1,
)

def write(self, grid_filename, mask_filename, metadata=None, variant=None):
"""
Write out CICE grid to netcdf
Expand All @@ -97,6 +108,13 @@ def write(self, grid_filename, mask_filename, metadata=None, variant=None):
if variant is not None and variant != "cice5-auscom":
raise NotImplementedError(f"{variant} not recognised")

has_bounds = (
(self.clon_t is not None)
& (self.clat_t is not None)
& (self.clon_u is not None)
& (self.clat_u is not None)
)

# Grid file
f = nc.Dataset(grid_filename, "w")

Expand All @@ -107,6 +125,8 @@ def write(self, grid_filename, mask_filename, metadata=None, variant=None):
f.createDimension(
"nj", self.num_lat_points
) # nj is the grid_latitude but doesn't have a value other than its index
if has_bounds:
f.createDimension("nvertices", 4) # nvertices is the number of corner points

# Make all CICE grid variables.
# names are based on https://cfconventions.org/Data/cf-standard-names/current/build/cf-standard-name-table.html
Expand All @@ -129,6 +149,28 @@ def write(self, grid_filename, mask_filename, metadata=None, variant=None):
tlon.long_name = "Longitude of T points"
tlon.standard_name = "longitude"

if has_bounds:
ulat_bounds = self._create_3d_nc_var(f, "ulat_bounds")
ulat_bounds.units = "degrees_north"
ulat_bounds.long_name = "Latitude of U cell vertices"
ulat_bounds.standard_name = "latitude_bounds"
ulat.bounds = "ulat_bounds"
ulon_bounds = self._create_3d_nc_var(f, "ulon_bounds")
ulon_bounds.units = "degrees_east"
ulon_bounds.long_name = "Longitude of U cell vertices"
ulon_bounds.standard_name = "longitude_bounds"
ulon.bounds = "ulon_bounds"
tlat_bounds = self._create_3d_nc_var(f, "tlat_bounds")
tlat_bounds.units = "degrees_north"
tlat_bounds.long_name = "Latitude of T cell vertices"
tlat_bounds.standard_name = "latitude_bounds"
tlat.bounds = "tlat_bounds"
tlon_bounds = self._create_3d_nc_var(f, "tlon_bounds")
tlon_bounds.units = "degrees_east"
tlon_bounds.long_name = "Longitude of T cell vertices"
tlon_bounds.standard_name = "longitude_bounds"
tlon.bounds = "tlon_bounds"

htn = self._create_2d_nc_var(f, "htn")
htn.units = "cm"
htn.long_name = "Width of T cells on North side."
Expand Down Expand Up @@ -180,6 +222,15 @@ def write(self, grid_filename, mask_filename, metadata=None, variant=None):
ulat[:] = np.deg2rad(self.y_u)
ulon[:] = np.deg2rad(self.x_u)

if has_bounds:
# In CICE, bounds are reported in deg, see
# https://github.com/ACCESS-NRI/cice5/blob/90a716400ba317fa230134c57ddaf7a84ff625d0/source/ice_grid.F90#L1968
# Transpose to (j, i, vertex)
tlat_bounds[:] = self.clat_t.transpose((1, 2, 0))
tlon_bounds[:] = self.clon_t.transpose((1, 2, 0))
ulat_bounds[:] = self.clat_u.transpose((1, 2, 0))
ulon_bounds[:] = self.clon_u.transpose((1, 2, 0))

# Convert from m to cm.
htn[:] = self.dx_tn[:] * 100.0
hte[:] = self.dy_te[:] * 100.0
Expand Down
79 changes: 23 additions & 56 deletions esmgrids/mom_grid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ def fromfile(cls, h_grid_def, v_grid_def=None, mask_file=None, calc_areas=True,
dx_ext = np.append(dx[:], dx[:, 0:1], axis=1)

# The u-cells cross the tri-polar fold
dy_ext = np.append(dy[:], dy[-1:, :], axis=0)
dy_ext = np.append(dy[:], np.fliplr(dy[-1:, :]), axis=0)

# Through the centre of u cells
dx_u = dx_ext[2::2, 1::2] + dx_ext[2::2, 2::2]
Expand All @@ -95,11 +95,11 @@ def fromfile(cls, h_grid_def, v_grid_def=None, mask_file=None, calc_areas=True,
# Add up areas, going clockwise from bottom left.
area_t = area[0::2, 0::2] + area[1::2, 0::2] + area[1::2, 1::2] + area[0::2, 1::2]

# These need to wrap around the globe. Copy ocn_area and
# add an extra column at the end. Also u-cells cross the
# tri-polar fold so add an extra row at the top.
area_ext = np.append(area[:], area[:, 0:1], axis=1)
area_ext = np.append(area_ext[:], area_ext[-1:, :], axis=0)
# Extend grid over periodic boundaries as u-cells cross the eastern and northern
# extents of the supergrid. Add an new row at the top that is a flipped version of
# the top row. Then add new column to the right that is a copy of the first column.
area_ext = np.append(area[:], np.fliplr(area[-1:, :]), axis=0)
area_ext = np.append(area_ext[:], area_ext[:, :1], axis=1)

area_u = area_ext[1::2, 1::2] + area_ext[2::2, 1::2] + area_ext[2::2, 2::2] + area_ext[1::2, 2::2]

Expand Down Expand Up @@ -188,64 +188,31 @@ def make_corners(x, y):
clat_t[3, :, :] = y[2::2, 0:-1:2]
assert not np.isnan(np.sum(clat_t))

# Corners of u cells. Index 0 is bottom left and then
# anti-clockwise.
# Extend grid over periodic boundaries as u-cells cross the eastern and northern extents of
# the supergrid. Add an new row at the top that is a flipped version of the
# second-to-top row. Then add new column to the right that is a copy of the second column.
x_ext = np.append(x[:], np.fliplr(x[-2:-1, :]), axis=0)
x_ext = np.append(x_ext[:], x_ext[:, 1:2], axis=1)

# Need to be careful with the edges.
# - Make the South most row of cells half size in the vertical.
# - West needs to wrap around.
# Do the easy bits first and then fix up below.
y_ext = np.append(y[:], np.fliplr(y[-2:-1, :]), axis=0)
y_ext = np.append(y_ext[:], y_ext[:, 1:2], axis=1)

# Corners of u cells. Index 0 is bottom left and then
# anti-clockwise.
clon_u = np.empty((4, nrow, ncol))
clon_u[:] = np.NAN
clon_u[0, 1:, 1:] = x[1:-2:2, 1:-2:2]
clon_u[1, 1:, 1:] = x[1:-2:2, 3::2]
clon_u[2, 1:, 1:] = x[3::2, 3::2]
clon_u[3, 1:, 1:] = x[3::2, 1:-2:2]

# Fix up bottom row excluding left most column
clon_u[0, 0, 1:] = x[0, 1:-2:2]
clon_u[1, 0, 1:] = x[0, 3::2]
clon_u[2, 0, 1:] = x[1, 3::2]
clon_u[3, 0, 1:] = x[1, 1:-2:2]

# Fix up leftmost column excluding bottom row
clon_u[0, 1:, 0] = x[1:-2:2, -1]
clon_u[1, 1:, 0] = x[1:-2:2, 1]
clon_u[2, 1:, 0] = x[3::2, 1]
clon_u[3, 1:, 0] = x[3::2, -1]

# Fix up the bottom left corner point
clon_u[0, 0, 0] = x[0, -1]
clon_u[1, 0, 0] = x[0, 1]
clon_u[2, 0, 0] = x[1, 1]
clon_u[3, 0, 0] = x[1, -1]
clon_u[0, :, :] = x_ext[1:-1:2, 1:-1:2]
clon_u[1, :, :] = x_ext[1:-1:2, 3::2]
clon_u[2, :, :] = x_ext[3::2, 3::2]
clon_u[3, :, :] = x_ext[3::2, 1:-1:2]
assert not np.isnan(np.sum(clon_u))

clat_u = np.empty((4, nrow, ncol))
clat_u[:] = np.NAN
clat_u[0, 1:, 1:] = y[1:-2:2, 1:-2:2]
clat_u[1, 1:, 1:] = y[1:-2:2, 3::2]
clat_u[2, 1:, 1:] = y[3::2, 3::2]
clat_u[3, 1:, 1:] = y[3::2, 1:-2:2]

# Fix up bottom row excluding left most column
clat_u[0, 0, 1:] = y[0, 1:-2:2]
clat_u[1, 0, 1:] = y[0, 3::2]
clat_u[2, 0, 1:] = y[1, 3::2]
clat_u[3, 0, 1:] = y[1, 1:-2:2]

# Fix up leftmost column excluding bottom row
clat_u[0, 1:, 0] = y[1:-2:2, -1]
clat_u[1, 1:, 0] = y[1:-2:2, 1]
clat_u[2, 1:, 0] = y[3::2, 1]
clat_u[3, 1:, 0] = y[3::2, -1]

# Fix up the bottom left corner point
clat_u[0, 0, 0] = y[0, -1]
clat_u[1, 0, 0] = y[0, 1]
clat_u[2, 0, 0] = y[1, 1]
clat_u[3, 0, 0] = y[1, -1]
clat_u[0, :, :] = y_ext[1:-1:2, 1:-1:2]
clat_u[1, :, :] = y_ext[1:-1:2, 3::2]
clat_u[2, :, :] = y_ext[3::2, 3::2]
clat_u[3, :, :] = y_ext[3::2, 1:-1:2]
assert not np.isnan(np.sum(clat_u))

return clat_t, clon_t, clat_u, clon_u, None, None
8 changes: 6 additions & 2 deletions test/test_cice_grid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,8 +151,8 @@ def test_cice_var_list(grids):

def test_cice_dims(grids):
# Test : Are the dim names consistent with cice history output?
assert set(grids["cice"].ds.dims) == set(
["ni", "nj"]
assert set(["ni", "nj"]).issubset(
set(grids["cice"].ds.dims)
), "cice dimension names should be 'ni','nj' to be consistent with history output"
assert grids["cice"].ds.sizes["ni"] == len(grids["test_ds"].nx)
assert grids["cice"].ds.sizes["nj"] == len(grids["test_ds"].ny)
Expand Down Expand Up @@ -183,6 +183,10 @@ def test_cice_grid_attributes(grids):
"ulon": {"standard_name": "longitude", "units": "radians"},
"tlat": {"standard_name": "latitude", "units": "radians"},
"tlon": {"standard_name": "longitude", "units": "radians"},
"ulat_bounds": {"standard_name": "latitude_bounds", "units": "degrees_north"},
"ulon_bounds": {"standard_name": "longitude_bounds", "units": "degrees_east"},
"tlat_bounds": {"standard_name": "latitude_bounds", "units": "degrees_north"},
"tlon_bounds": {"standard_name": "longitude_bounds", "units": "degrees_east"},
"uarea": {
"standard_name": "cell_area",
"units": "m^2",
Expand Down