diff --git a/esmgrids/cice_grid.py b/esmgrids/cice_grid.py index 6639b2d..cbba6a6 100644 --- a/esmgrids/cice_grid.py +++ b/esmgrids/cice_grid.py @@ -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, @@ -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 @@ -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") @@ -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 @@ -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." @@ -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 diff --git a/esmgrids/mom_grid.py b/esmgrids/mom_grid.py index ec4297b..085f314 100644 --- a/esmgrids/mom_grid.py +++ b/esmgrids/mom_grid.py @@ -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] @@ -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] @@ -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 diff --git a/test/test_cice_grid.py b/test/test_cice_grid.py index 654f815..0b780f3 100644 --- a/test/test_cice_grid.py +++ b/test/test_cice_grid.py @@ -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) @@ -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",