Repository: British-Oceanographic-Data-Centre/COAsT (language: python)

Issue:
606 circulation and fluxes in global model data
<!-- Please refer to our contributing documentation for any questions on submitting a pull request, or let us know here if you need any help:  https://british-oceanographic-data-centre.github.io/COAsT/docs/contributing_package/ -->

## Pull request checklist

Please check if your PR fulfills the following requirements:
- [x ] Tests for the changes have been added (for bug fixes / features)
- [ ] Docs have been reviewed and added / updated if needed (for bug fixes / features)
    - Some docs updates need to be made in the `COAsT-site` repo, in a separate PR. See [contributing to documentation ](https://british-oceanographic-data-centre.github.io/COAsT/docs/contributing-docs/) for details.
- [ ] Build (`./build.sh`) was run locally and no errors reported. NB not sure about this requirement: GitActions test this
- [ ] Lint (`pylint .`) has passed locally and any fixes were made for failures. NB not sure about this requirement: GitActions test this with `black`


## Pull request type

<!-- Please do not submit updates to dependencies unless it fixes an issue. --> 

<!-- Please try to limit your pull request to one type, submit multiple pull requests if needed. --> 

Please check the type of change your PR introduces:
- [ ] Bugfix
- [ x] Feature
- [ ] Code style update (formatting, renaming)
- [ ] Refactoring (no functional changes, no api changes)
- [ ] Build related changes
- [ ] Documentation content changes
- [ ] Other (please describe):


## What is the current behavior?
<!-- Please describe the current behavior that you are modifying. -->


<!-- Issues are required for both bug fixes and features. -->
Issue URL: #


## What is the new behavior?
<!-- Please describe the behavior or changes that are being added by this PR. -->

-
-
### Checklist
- [ ] Code compiles correctly
- [ ] Created tests which fail without the change (if possible)
- [ ] All tests passing (`pip install . && pytest unit_testing/unit_test.py -s`)
- [ ] Extended the README / documentation, if necessary
- [ ] Added myself / the copyright holder to the AUTHORS file

## Does this introduce a breaking change?

- [ ] Yes
- [ ] No

<!-- If this introduces a breaking change, please describe the impact and migration path for existing applications below. -->


## Other information

<!-- Any other information that is important to this PR such as screenshots of how the component looks before and after the change. -->

Quick update to add useful feature of subsetting for global models

and also some circulation calculations / plot, which will be further developed 




Discussion / hints:

Fix patch:
diff --git a/CITATION.cff b/CITATION.cff
index 1cd03df8..b230f561 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -40,6 +40,6 @@ authors:
  - family-names: "de Mora"
    given-names: "Lee"
    orcid: "https://orcid.org/0000-0002-5080-3149"
-title: "British-Oceanographic-Data-Centre/COAsT: v3.2.0"
-version: v3.2.0
-date-released: 2022-11-23
+title: "British-Oceanographic-Data-Centre/COAsT: v3.2.1"
+version: v3.2.1
+date-released: 2023-04-04
diff --git a/coast/__init__.py b/coast/__init__.py
index 37279aa6..46a20af8 100644
--- a/coast/__init__.py
+++ b/coast/__init__.py
@@ -28,6 +28,7 @@
 from .data.copernicus import Copernicus, Product
 from ._utils.experiments_file_handling import experiments
 from ._utils.experiments_file_handling import nemo_filename_maker
+from .diagnostics.circulation import CurrentsOnT
 
 # Set default for logging level when coast is imported
 import logging
diff --git a/coast/_utils/crps_util.py b/coast/_utils/crps_util.py
index aca5e500..ca5eec00 100644
--- a/coast/_utils/crps_util.py
+++ b/coast/_utils/crps_util.py
@@ -188,7 +188,6 @@ def crps_sonf_fixed(
     else:
         # Subset model data in time and space: model -> obs
         for ii in neighbourhood_indices:
-
             mod_subset_time = mod_subset.interp(
                 time=obs_time[ii], method=time_interp, kwargs={"fill_value": "extrapolate"}
             )
diff --git a/coast/_utils/docsy_tools.py b/coast/_utils/docsy_tools.py
index 131ff5f2..9898d520 100644
--- a/coast/_utils/docsy_tools.py
+++ b/coast/_utils/docsy_tools.py
@@ -12,7 +12,6 @@ def __init__(self):
     def write_class_to_markdown(
         cls, class_to_write, fn_out, method_to_omit=[], omit_private_methods=True, omit_parent_methods=True
     ):
-
         methods_to_write = cls._get_list_of_methods(class_to_write)
 
         for method in methods_to_write:
diff --git a/coast/_utils/mask_maker.py b/coast/_utils/mask_maker.py
index 0ac0f74a..139b84bc 100644
--- a/coast/_utils/mask_maker.py
+++ b/coast/_utils/mask_maker.py
@@ -36,7 +36,6 @@ class MaskMaker:
     """
 
     def __init__(self):
-
         return
 
     @staticmethod
@@ -130,9 +129,9 @@ def fill_polygon_by_lonlat(
         return array_to_fill
 
     @classmethod
-    def region_def_nws_north_sea(cls, longitude, latitude, bath):
+    def region_def_nws_north_north_sea(cls, longitude, latitude, bath):
         """
-        Regional definition for the North Sea (Northwest European Shelf)
+        Regional definition for the northern North Sea (Northwest European Shelf)
         Longitude, latitude and bath should be 2D arrays corresponding to model
         coordinates and bathymetry. Bath should be positive with depth.
         """
@@ -143,6 +142,14 @@ def region_def_nws_north_sea(cls, longitude, latitude, bath):
         mask = mask * (bath < 200) * (bath > 0) * (~np.isnan(bath))
         return mask
 
+    @classmethod
+    def region_def_nws_south_north_sea(cls, longitude, latitude, bath):
+        vertices_lon = [-0.67, -0.67, 9, 9, 7.57, 7.57]
+        vertices_lat = [54.08, 51, 51, 56, 56, 54.08]
+        mask = cls.fill_polygon_by_lonlat(np.zeros(longitude.shape), longitude, latitude, vertices_lon, vertices_lat)
+        mask = mask * (bath < 200) * (bath > 0) * (~np.isnan(bath))
+        return mask
+
     @classmethod
     def region_def_nws_outer_shelf(cls, longitude, latitude, bath):
         """
@@ -183,15 +190,7 @@ def region_def_nws_english_channel(cls, longitude, latitude, bath):
         return mask
 
     @classmethod
-    def region_def_south_north_sea(cls, longitude, latitude, bath):
-        vertices_lon = [-0.67, -0.67, 9, 9, 7.57, 7.57]
-        vertices_lat = [54.08, 51, 51, 56, 56, 54.08]
-        mask = cls.fill_polygon_by_lonlat(np.zeros(longitude.shape), longitude, latitude, vertices_lon, vertices_lat)
-        mask = mask * (bath < 200) * (bath > 0) * (~np.isnan(bath))
-        return mask
-
-    @classmethod
-    def region_def_off_shelf(cls, longitude, latitude, bath):
+    def region_def_nws_off_shelf(cls, longitude, latitude, bath):
         vertices_lon = [10, 10, -5, -10, 0, 0, -20, -20]
         vertices_lat = [65, 60, 59, 52.5, 47.5, 45, 40, 63]
         mask = cls.fill_polygon_by_lonlat(np.zeros(longitude.shape), longitude, latitude, vertices_lon, vertices_lat)
@@ -199,7 +198,7 @@ def region_def_off_shelf(cls, longitude, latitude, bath):
         return mask
 
     @classmethod
-    def region_def_irish_sea(cls, longitude, latitude, bath):
+    def region_def_nws_irish_sea(cls, longitude, latitude, bath):
         vertices_lon = [-5, -7.6, -7.5, -4.1, 0, -2.6]
         vertices_lat = [56.4, 55, 52, 50.7, 51.5, 55.3]
         mask = cls.fill_polygon_by_lonlat(np.zeros(longitude.shape), longitude, latitude, vertices_lon, vertices_lat)
@@ -207,13 +206,26 @@ def region_def_irish_sea(cls, longitude, latitude, bath):
         return mask
 
     @classmethod
-    def region_def_kattegat(cls, longitude, latitude, bath):
+    def region_def_nws_kattegat(cls, longitude, latitude, bath):
         vertices_lon = [9, 9, 13, 13]
         vertices_lat = [60, 52.5, 52.5, 60]
         mask = cls.fill_polygon_by_lonlat(np.zeros(longitude.shape), longitude, latitude, vertices_lon, vertices_lat)
         mask = mask * (bath < 200) * (bath > 0) * (~np.isnan(bath))
         return mask
 
+    def region_def_nws_fsc(cls, longitude, latitude, bath):
+        """
+        Regional definition for Faroe Shetland Channel (Northwest European Shelf)
+        Longitude, latitude and bath should be 2D arrays corresponding to model
+        coordinates and bathymetry. Bath should be positive with depth.
+        """
+        vertices_lon = [-7.13, -9.72, -6.37, -0.45, -4.53]
+        vertices_lat = [62.17, 60.6, 59.07, 61.945, 62.51]
+
+        mask = cls.fill_polygon_by_lonlat(np.zeros(longitude.shape), longitude, latitude, vertices_lon, vertices_lat)
+        mask = mask * (bath > 200) * (bath > 0) * (~np.isnan(bath))
+        return mask
+
     @classmethod
     def make_region_from_vertices(cls, longitude, latitude, vertices_lon: list, vertices_lat: list):
         """
diff --git a/coast/_utils/plot_util.py b/coast/_utils/plot_util.py
index a281781e..65d9adef 100644
--- a/coast/_utils/plot_util.py
+++ b/coast/_utils/plot_util.py
@@ -222,7 +222,6 @@ def create_geo_axes(lonbounds, latbounds):
 
 
 def ts_diagram(temperature, salinity, depth):
-
     fig = plt.figure(figsize=(10, 7))
     ax = plt.scatter(salinity, temperature, c=depth)
     cbar = plt.colorbar()
diff --git a/coast/data/altimetry.py b/coast/data/altimetry.py
index 2f02f767..17a7f49f 100644
--- a/coast/data/altimetry.py
+++ b/coast/data/altimetry.py
@@ -234,7 +234,6 @@ def crps(
         time_interp: str = "linear",
         create_new_object=True,
     ):
-
         """
         Comparison of observed variable to modelled using the Continuous
         Ranked Probability Score. This is done using this ALTIMETRY object.
diff --git a/coast/data/gridded.py b/coast/data/gridded.py
index 0d8ba322..786d74ec 100644
--- a/coast/data/gridded.py
+++ b/coast/data/gridded.py
@@ -63,18 +63,22 @@ def _setup_grid_obj(self, chunks, multiple, **kwargs):
 
         Args:
             chunks: This is a setting for xarray as to whether dask (parrell processing) should be on and how it works
-            multiple: falg to tell if we are loading one or more files
+            multiple: flag to tell if we are loading one or more files
             **kwargs: pass direct to loaded xarray dataset
+                lims = [x_dim index 1, x_dim_index 2, y_dim index 1, y_dim_index 2] - subset region defined from
+                                                                            lower left to upper right corners
+                calculate_bathymetry [boolean]: default-False
         """
         self.set_grid_vars()
         self.set_dimension_mapping()
         self.set_variable_mapping()
-
+        lims = kwargs.get("lims", [])
         if self.fn_data is not None:
             self.load(self.fn_data, chunks, multiple)
 
-        self.set_dimension_names(self.config.dataset.dimension_map)
-        self.set_variable_names(self.config.dataset.variable_map)
+            self.set_dimension_names(self.config.dataset.dimension_map)
+            self.set_variable_names(self.config.dataset.variable_map)
+            self.dataset = self.spatial_subset(self.dataset, lims)  # Trim data size if indices specified
 
         if self.fn_domain is None:
             self.filename_domain = ""  # empty store for domain fileanme
@@ -88,8 +92,9 @@ def _setup_grid_obj(self, chunks, multiple, **kwargs):
             for key, value in kwargs.items():
                 dataset_domain[key] = value
 
+            dataset_domain = self.spatial_subset(dataset_domain, lims)  # Trim domain size if indices specified
             if self.fn_data is not None:
-                dataset_domain = self.trim_domain_size(dataset_domain)
+                dataset_domain = self.trim_domain_size(dataset_domain)  # Trim domain size if self.data is smaller
             self.set_timezero_depths(
                 dataset_domain, **kwargs
             )  # THIS ADDS TO dataset_domain. Should it be 'return'ed (as in trim_domain_size) or is implicit OK?
@@ -189,7 +194,7 @@ def get_contour_complex(self, var, points_x, points_y, points_z, tolerance: int
         smaller = self.dataset[var].sel(z=points_z, x=points_x, y=points_y, method="nearest", tolerance=tolerance)
         return smaller
 
-    def set_timezero_depths(self, dataset_domain, calculate_bathymetry=False):
+    def set_timezero_depths(self, dataset_domain, **kwargs):
         """
         Calculates the depths at time zero (from the domain_cfg input file)
         for the appropriate grid.
@@ -204,6 +209,7 @@ def set_timezero_depths(self, dataset_domain, calculate_bathymetry=False):
         # keyword to allow calcution of bathymetry from scale factors
 
         # All bathymetry should now be mapped to bathy_metry
+        calculate_bathymetry = kwargs.get("calculate_bathymetry", False)
         try:
             if calculate_bathymetry:  # calculate bathymetry from scale factors
                 bathymetry, mask, time_mask = self.calc_bathymetry(dataset_domain)
@@ -510,7 +516,6 @@ def interpolate_in_time(model_array, new_times, interp_method="nearest", extrapo
     def construct_density(
         self, eos="EOS10", rhobar=False, Zd_mask=[], CT_AS=False, pot_dens=False, Tbar=True, Sbar=True
     ):
-
         """
             Constructs the in-situ density using the salinity, temperture and
             depth_0 fields and adds a density attribute to the t-grid dataset
@@ -675,6 +680,27 @@ def construct_density(
         except AttributeError as err:
             error(err)
 
+    def spatial_subset(self, dataset, lims):
+        """
+        Specify indices to subset the data. Subset region defined as a 2D box from lower left to upper right corners
+        lims = [x_dim index_1, x_dim_index_2, y_dim index_1, y_dim_index_2] -
+        Modifies self.dataset
+        """
+        if len(lims) == 4:  # if lims are provided take a subset
+            debug(f"Trimming by indices: {lims}")
+            # subsetting will wrap longitude across dateline if 1st longitude is larger than 2nd.
+            if "x_dim" in dataset.dims:
+                if lims[0] < lims[1]:  # usual case
+                    dataset = dataset.isel(y_dim=range(lims[2], lims[3]), x_dim=range(lims[0], lims[1]))
+                else:  # longitude  wrap around
+                    nx = dataset.dims["x_dim"]
+                    ds1 = dataset.isel(y_dim=range(lims[2], lims[3]), x_dim=range(lims[0], nx))
+                    ds2 = dataset.isel(y_dim=range(lims[2], lims[3]), x_dim=range(0, lims[1]))
+                    dataset = xr.concat([ds1, ds2], dim="x_dim")
+            else:
+                print("limits not used as only work with datasets having dimension x_dim")
+        return dataset
+
     def trim_domain_size(self, dataset_domain):
         """
         Trim the domain variables if the dataset object is a spatial subset
@@ -697,12 +723,22 @@ def trim_domain_size(self, dataset_domain):
             )
 
             # Find the corners of the cut out domain.
-            [j0, i0] = self.find_j_i_domain(
-                lat=self.dataset.nav_lat[0, 0], lon=self.dataset.nav_lon[0, 0], dataset_domain=dataset_domain
-            )
-            [j1, i1] = self.find_j_i_domain(
-                lat=self.dataset.nav_lat[-1, -1], lon=self.dataset.nav_lon[-1, -1], dataset_domain=dataset_domain
-            )
+            try:
+                [j0, i0] = self.find_j_i_domain(
+                    lat=self.dataset.latitude[0, 0], lon=self.dataset.longitude[0, 0], dataset_domain=dataset_domain
+                )
+                [j1, i1] = self.find_j_i_domain(
+                    lat=self.dataset.latitude[-1, -1], lon=self.dataset.longitude[-1, -1], dataset_domain=dataset_domain
+                )
+                debug(f"trim_domain_size(): USED dataset.longitude")
+            except:  # if called before variables are re-mapped. Not very pretty...
+                [j0, i0] = self.find_j_i_domain(
+                    lat=self.dataset.nav_lat[0, 0], lon=self.dataset.nav_lon[0, 0], dataset_domain=dataset_domain
+                )
+                [j1, i1] = self.find_j_i_domain(
+                    lat=self.dataset.nav_lat[-1, -1], lon=self.dataset.nav_lon[-1, -1], dataset_domain=dataset_domain
+                )
+                debug(f"trim_domain_size(): USED dataset.nav_lon")
 
             dataset_subdomain = dataset_domain.isel(y_dim=slice(j0, j1 + 1), x_dim=slice(i0, i1 + 1))
             return dataset_subdomain
diff --git a/coast/data/profile.py b/coast/data/profile.py
index 9b94e74b..52bf3550 100644
--- a/coast/data/profile.py
+++ b/coast/data/profile.py
@@ -111,7 +111,7 @@ def read_en4(self, fn_en4, chunks: dict = {}, multiple=False) -> None:
 
             # Reorder files to read
             file_to_read = np.array(file_to_read)
-            dates = [ff[-9:-3] for ff in file_to_read]
+            dates = [ff[-9:-3] for ff in file_to_read]  # Assumes monthly filename structure: EN*yyyymm.nc
             dates = [datetime.datetime(int(dd[0:4]), int(dd[4:6]), 1) for dd in dates]
             sort_ind = np.argsort(dates)
             file_to_read = file_to_read[sort_ind]
@@ -156,7 +156,6 @@ def subset_indices_lonlat_box(self, lonbounds, latbounds):
     """======================= Plotting ======================="""
 
     def plot_profile(self, var: str, profile_indices=None):
-
         fig = plt.figure(figsize=(7, 10))
 
         if profile_indices is None:
@@ -175,7 +174,6 @@ def plot_profile(self, var: str, profile_indices=None):
         return fig, ax
 
     def plot_map(self, var_str=None):
-
         profiles = self.dataset
 
         if var_str is None:
@@ -186,7 +184,6 @@ def plot_map(self, var_str=None):
         return fig, ax
 
     def plot_ts_diagram(self, profile_index, var_t="potential_temperature", var_s="practical_salinity"):
-
         profile = self.dataset.isel(id_dim=profile_index)
         temperature = profile[var_t].values
         salinity = profile[var_s].values
@@ -197,7 +194,7 @@ def plot_ts_diagram(self, profile_index, var_t="potential_temperature", var_s="p
 
     """======================= Model Comparison ======================="""
 
-    def process_en4(self, sort_time=True):
+    def process_en4(self, sort_time=True, remove_flagged_neighbours=False):
         """
         VERSION 1.4 (05/07/2021)
 
@@ -224,6 +221,12 @@ def process_en4(self, sort_time=True):
          fn_out (str)      : Full path to a desired output file. If unspecified
                              then nothing is written.
 
+        remove_flagged_neighbours: EN offers a profile flag that indicates there are
+                other profiles within 0.2 deg of latitude and longitude and 1 hour that
+                appear to be of higher quality. In previous versions of the dataset these
+                profiles would have not been stored in the data files. Setting this flag
+                as True removes these profiles.
+
         EXAMPLE USEAGE:
          profile = coast.PROFILE()
          profile.read_EN4(fn_en4, chunks={'N_PROF':10000})
@@ -243,12 +246,18 @@ def process_en4(self, sort_time=True):
         # Each bit of this string is a different QC flag. Which flag is which can
         # be found on the EN4 website:
         # https://www.metoffice.gov.uk/hadobs/en4/en4-0-2-profile-file-format.html
-        qc_str = [np.binary_repr(ds.qc_flags_profiles.values[pp]).zfill(30)[::-1] for pp in range(ds.sizes["id_dim"])]
+        qc_str = [
+            np.binary_repr(ds.qc_flags_profiles.astype(int).values[pp]).zfill(32)[::-1]
+            for pp in range(ds.sizes["id_dim"])
+        ]
 
         # Determine indices of the profiles that we want to keep
         reject_tem_prof = np.array([int(qq[0]) for qq in qc_str], dtype=bool)
         reject_sal_prof = np.array([int(qq[1]) for qq in qc_str], dtype=bool)
         reject_both_prof = np.logical_and(reject_tem_prof, reject_sal_prof)
+        if remove_flagged_neighbours:
+            reject_close_flagged_prof = np.array([int(qq[2]) for qq in qc_str], dtype=bool)
+            reject_both_prof = np.logical_or(reject_both_prof, reject_close_flagged_prof)
         ds["reject_tem_prof"] = (["id_dim"], reject_tem_prof)
         ds["reject_sal_prof"] = (["id_dim"], reject_sal_prof)
         debug("     >>> QC: Completely rejecting {0} / {1} id_dims".format(np.sum(reject_both_prof), ds.dims["id_dim"]))
@@ -257,19 +266,25 @@ def process_en4(self, sort_time=True):
         ds = ds.isel(id_dim=~reject_both_prof)
         reject_tem_prof = reject_tem_prof[~reject_both_prof]
         reject_sal_prof = reject_sal_prof[~reject_both_prof]
-
-        # Get new QC flags array
-        qc_lev = ds.qc_flags_levels.values
-
         debug(f" QC: Additional profiles converted to NaNs: ")
         debug(f"     >>> {0} temperature profiles ".format(np.sum(reject_tem_prof)))
         debug(f"     >>> {0} salinity profiles ".format(np.sum(reject_sal_prof)))
 
+        debug(f"MASKING rejected profiles, replacing with NaNs...")
+        reject_tem_prof_arr = ds.reject_tem_prof.broadcast_like(ds.potential_temperature).values
+        reject_sal_prof_arr = ds.reject_sal_prof.broadcast_like(ds.practical_salinity).values
+        ds["temperature"] = xr.where(~reject_tem_prof_arr, ds["temperature"], np.nan)
+        ds["potential_temperature"] = xr.where(~reject_tem_prof_arr, ds["potential_temperature"], np.nan)
+        ds["practical_salinity"] = xr.where(~reject_sal_prof_arr, ds["practical_salinity"], np.nan)
+
+        # Get new QC flags array
+        qc_lev = ds.qc_flags_levels.values
         #
         reject_tem_lev = np.zeros((ds.dims["id_dim"], ds.dims["z_dim"]), dtype=bool)
         reject_sal_lev = np.zeros((ds.dims["id_dim"], ds.dims["z_dim"]), dtype=bool)
 
-        int_tem, int_sal, int_both = self.calculate_all_en4_qc_flags()
+        # Get lists of QC integers that marks levels for masking
+        int_tem, int_sal, int_both = self.calculate_en4_qc_flags_levels()
         for ii in range(len(int_tem)):
             reject_tem_lev[qc_lev == int_tem[ii]] = 1
         for ii in range(len(int_sal)):
@@ -283,8 +298,8 @@ def process_en4(self, sort_time=True):
 
         debug(f"MASKING rejected datapoints, replacing with NaNs...")
         ds["temperature"] = xr.where(~reject_tem_lev, ds["temperature"], np.nan)
-        ds["potential_temperature"] = xr.where(~reject_tem_lev, ds["temperature"], np.nan)
-        ds["practical_salinity"] = xr.where(~reject_tem_lev, ds["practical_salinity"], np.nan)
+        ds["potential_temperature"] = xr.where(~reject_tem_lev, ds["potential_temperature"], np.nan)
+        ds["practical_salinity"] = xr.where(~reject_sal_lev, ds["practical_salinity"], np.nan)
 
         if sort_time:
             debug(f"Sorting Time Dimension...")
@@ -296,69 +311,6 @@ def process_en4(self, sort_time=True):
         return_prof.dataset = ds
         return return_prof
 
-    @classmethod
-    def calculate_all_en4_qc_flags(cls):
-        """
-        Brute force method for identifying all rejected points according to
-        EN4 binary integers. It can be slow to convert large numbers of integers
-        to a sequence of bits and is actually quicker to just generate every
-        combination of possible QC integers. That's what this routine does.
-        Used in PROFILE.preprocess_en4().
-
-        INPUTS
-         NO INPUTS
-
-        OUTPUTS
-         qc_integers_tem  : Array of integers signifying the rejection of ONLY
-                            temperature datapoints
-         qc_integers_sal  : Array of integers signifying the rejection of ONLY
-                            salinity datapoints
-         qc_integers_both : Array of integers signifying the rejection of BOTH
-                            temperature and salinity datapoints.
-        """
-
-        reject_tem_ind = 0
-        reject_sal_ind = 1
-        reject_tem_reasons = [2, 3, 8, 9, 10, 11, 12, 13, 14, 15, 16]
-        reject_sal_reasons = [2, 3, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]
-
-        qc_integers_tem = []
-        qc_integers_sal = []
-        qc_integers_both = []
-        n_tem_reasons = len(reject_tem_reasons)
-        n_sal_reasons = len(reject_sal_reasons)
-        bin_len = 30
-
-        # IF reject_tem = 1, reject_sal = 0
-        for ii in range(n_tem_reasons):
-            bin_tmp = np.zeros(bin_len, dtype=int)
-            bin_tmp[reject_tem_ind] = 1
-            bin_tmp[reject_tem_reasons[ii]] = 1
-            qc_integers_tem.append(int("".join(str(jj) for jj in bin_tmp)[::-1], 2))
-
-        # IF reject_tem = 0, reject_sal = 1
-        for ii in range(n_sal_reasons):
-            bin_tmp = np.zeros(bin_len, dtype=int)
-            bin_tmp[reject_sal_ind] = 1
-            bin_tmp[reject_sal_reasons[ii]] = 1
-            qc_integers_sal.append(int("".join(str(jj) for jj in bin_tmp)[::-1], 2))
-
-        # IF reject_tem = 1, reject_sal = 1
-        for tt in range(n_tem_reasons):
-            for ss in range(n_sal_reasons):
-                bin_tmp = np.zeros(bin_len, dtype=int)
-                bin_tmp[reject_tem_ind] = 1
-                bin_tmp[reject_sal_ind] = 1
-                bin_tmp[reject_tem_reasons[tt]] = 1
-                bin_tmp[reject_sal_reasons[ss]] = 1
-                qc_integers_both.append(int("".join(str(jj) for jj in bin_tmp)[::-1], 2))
-
-        qc_integers_tem = list(set(qc_integers_tem))
-        qc_integers_sal = list(set(qc_integers_sal))
-        qc_integers_both = list(set(qc_integers_both))
-
-        return qc_integers_tem, qc_integers_sal, qc_integers_both
-
     def obs_operator(self, gridded, mask_bottom_level=True):
         """
         VERSION 2.0 (04/10/2021)
@@ -509,110 +461,13 @@ def obs_operator(self, gridded, mask_bottom_level=True):
         mod_profiles["nearest_index_t"] = (["id_dim"], ind_t.values)
         return Profile(dataset=mod_profiles)
 
-    def process_en4(self, sort_time=True):
-        """
-        VERSION 1.4 (05/07/2021)
-
-        PREPROCESSES EN4 data ready for comparison with model data.
-        This routine will cut out a desired geographical box of EN4 data and
-        then apply quality control according to the available flags in the
-        netCDF files. Quality control happens in two steps:
-            1. Where a whole data profile is flagged, it is completely removed
-               from the dataset
-            2. Where a single datapoint is rejected in either temperature or
-               salinity, it is set to NaN.
-        This routine attempts to use xarray/dask chunking magic to keep
-        memory useage low however some memory is still needed for loading
-        flags etc. May be slow if using large EN4 datasets.
-
-        Routine will return a processed profile object dataset and can write
-        the new dataset to file if fn_out is defined. If saving to the
-        PROFILE object, be aware that DASK computations will not have happened
-        and will need to be done using .load(), .compute() or similar before
-        accessing the values. IF using multiple EN4 files or large dataset,
-        make sure you have chunked the data over N_PROF dimension.
-
-        INPUTS
-         fn_out (str)      : Full path to a desired output file. If unspecified
-                             then nothing is written.
-
-        EXAMPLE USEAGE:
-         profile = coast.PROFILE()
-         profile.read_EN4(fn_en4, chunks={'N_PROF':10000})
-         fn_out = '~/output_file.nc'
-         new_profile = profile.preprocess_en4(fn_out = fn_out,
-                                              lonbounds = [-10, 10],
-                                              latbounds = [45, 65])
-        """
-
-        ds = self.dataset
-
-        # REJECT profiles that are QC flagged.
-        debug(f" Applying QUALITY CONTROL to EN4 data...")
-        ds.qc_flags_profiles.load()
-
-        # This line reads converts the QC integer to a binary string.
-        # Each bit of this string is a different QC flag. Which flag is which can
-        # be found on the EN4 website:
-        # https://www.metoffice.gov.uk/hadobs/en4/en4-0-2-profile-file-format.html
-        qc_str = [np.binary_repr(ds.qc_flags_profiles.values[pp]).zfill(30)[::-1] for pp in range(ds.dims["id_dim"])]
-
-        # Determine indices of kept profiles
-        reject_tem_prof = np.array([int(qq[0]) for qq in qc_str], dtype=bool)
-        reject_sal_prof = np.array([int(qq[1]) for qq in qc_str], dtype=bool)
-        reject_both_prof = np.logical_and(reject_tem_prof, reject_sal_prof)
-        ds["reject_tem_prof"] = (["id_dim"], reject_tem_prof)
-        ds["reject_sal_prof"] = (["id_dim"], reject_sal_prof)
-        debug(
-            "     >>> QC: Completely rejecting {0} / {1} profiles".format(np.sum(reject_both_prof), ds.dims["id_dim"])
-        )
-
-        ds = ds.isel(id_dim=~reject_both_prof)
-        reject_tem_prof = reject_tem_prof[~reject_both_prof]
-        reject_sal_prof = reject_sal_prof[~reject_both_prof]
-        qc_lev = ds.qc_flags_levels.values
-
-        debug(f" QC: Additional profiles converted to NaNs: ")
-        debug(f"     >>> {0} temperature profiles ".format(np.sum(reject_tem_prof)))
-        debug(f"     >>> {0} salinity profiles ".format(np.sum(reject_sal_prof)))
-
-        reject_tem_lev = np.zeros((ds.dims["id_dim"], ds.dims["z_dim"]), dtype=bool)
-        reject_sal_lev = np.zeros((ds.dims["id_dim"], ds.dims["z_dim"]), dtype=bool)
-
-        int_tem, int_sal, int_both = self.calculate_all_en4_qc_flags()
-        for ii in range(len(int_tem)):
-            reject_tem_lev[qc_lev == int_tem[ii]] = 1
-        for ii in range(len(int_sal)):
-            reject_sal_lev[qc_lev == int_sal[ii]] = 1
-        for ii in range(len(int_both)):
-            reject_tem_lev[qc_lev == int_both[ii]] = 1
-            reject_sal_lev[qc_lev == int_both[ii]] = 1
-
-        ds["reject_tem_datapoint"] = (["id_dim", "z_dim"], reject_tem_lev)
-        ds["reject_sal_datapoint"] = (["id_dim", "z_dim"], reject_sal_lev)
-
-        debug(f"MASKING rejected datapoints, replacing with NaNs...")
-        ds["temperature"] = xr.where(~reject_tem_lev, ds["temperature"], np.nan)
-        ds["potential_temperature"] = xr.where(~reject_tem_lev, ds["temperature"], np.nan)
-        ds["practical_salinity"] = xr.where(~reject_tem_lev, ds["practical_salinity"], np.nan)
-
-        if sort_time:
-            debug(f"Sorting Time Dimension...")
-            ds = ds.sortby("time")
-
-        debug(f"Finished processing data. Returning new Profile object.")
-
-        return_prof = Profile()
-        return_prof.dataset = ds
-        return return_prof
-
-    def calculate_all_en4_qc_flags(self):
+    def calculate_en4_qc_flags_levels(self):
         """
         Brute force method for identifying all rejected points according to
         EN4 binary integers. It can be slow to convert large numbers of integers
         to a sequence of bits and is actually quicker to just generate every
         combination of possible QC integers. That's what this routine does.
-        Used in PROFILE.preprocess_en4().
+        Used in Profile.process_en4().
 
         INPUTS
          NO INPUTS
@@ -628,15 +483,27 @@ def calculate_all_en4_qc_flags(self):
 
         reject_tem_ind = 0
         reject_sal_ind = 1
-        reject_tem_reasons = [2, 3, 8, 9, 10, 11, 12, 13, 14, 15, 16]
-        reject_sal_reasons = [2, 3, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]
+        if "4.2.0" in self.dataset.history:
+            # from https://www.metoffice.gov.uk/hadobs/en4/en4-0-2-profile-file-format.html
+            reject_tem_reasons = [2, 3, 8, 9, 10, 11, 12, 13, 14, 15, 16]
+            reject_sal_reasons = [2, 3, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]
+        elif "4.2.2" in self.dataset.history:
+            # from https://www.metoffice.gov.uk/hadobs/en4/en4-2-2-profile-file-format.html
+            reject_tem_reasons = [2, 3, 8, 9, 10, 11, 12, 13, 14, 15, 16]
+            reject_sal_reasons = [2, 3, 21, 22, 23, 24, 25, 26, 27, 28, 29]
+        else:
+            reject_tem_reasons = [2, 3, 8, 9, 10, 11, 12, 13, 14, 15, 16]
+            reject_sal_reasons = [2, 3, 21, 22, 23, 24, 25, 26, 27, 28, 29]
+            debug(
+                f"Assume QC flags following 4.2.2: https://www.metoffice.gov.uk/hadobs/en4/en4-2-2-profile-file-format.html"
+            )
 
         qc_integers_tem = []
         qc_integers_sal = []
         qc_integers_both = []
         n_tem_reasons = len(reject_tem_reasons)
         n_sal_reasons = len(reject_sal_reasons)
-        bin_len = 30
+        bin_len = 32
 
         # IF reject_tem = 1, reject_sal = 0
         for ii in range(n_tem_reasons):
diff --git a/coast/data/tidegauge.py b/coast/data/tidegauge.py
index fa54e2e9..a62a096a 100644
--- a/coast/data/tidegauge.py
+++ b/coast/data/tidegauge.py
@@ -1109,12 +1109,14 @@ def obs_operator(self, gridded, time_interp="nearest"):
 
         # Determine if landmask is present
         if "landmask" not in gridded:
-            gridded.landmask = None
+            mask = None
+        else:
+            mask = gridded["landmask"]
 
         # Determine spatial indices
         print("Calculating spatial indices.", flush=True)
         ind_x, ind_y = general_utils.nearest_indices_2d(
-            gridded.longitude, gridded.latitude, ds.longitude, ds.latitude, mask=gridded.landmask
+            gridded.longitude, gridded.latitude, ds.longitude, ds.latitude, mask=mask
         )
 
         # Extract spatial time series
diff --git a/coast/diagnostics/circulation.py b/coast/diagnostics/circulation.py
new file mode 100644
index 00000000..aca5690b
--- /dev/null
+++ b/coast/diagnostics/circulation.py
@@ -0,0 +1,148 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Dec 21 16:47:59 2022
+
+@author: jholt
+"""
+from ..data.gridded import Gridded
+import numpy as np
+import xarray as xr
+import matplotlib.pylab as plt
+from .._utils.logging_util import debug
+
+
+class CurrentsOnT(Gridded):
+    """
+    Methods for currents co-located on T-points
+
+    Use like coast.Gridded() with a domain_cfg file name to create gridded object
+
+    """
+
+    def __init(self, fn_domain=None, config=None, **kwargs):
+        gridded = Gridded(fn_domain=fn_domain, config=config, **kwargs)
+        self.dataset = gridded.dataset
+
+    def currents_on_t(self, gridded_u, gridded_v):
+        """
+        Adds co-located velocity components and speed to CurrentsOnT object
+
+        Avoids use of indices in the hope that it will be more parallelisable. However in index space the mappings
+        take this form:
+        u_on_t_points[:, :, :, 1:] = 0.5 * (
+            gridded_u.dataset.u_velocity.values[:, :, :, 1:] + gridded_u.dataset.u_velocity.values[:, :, :, :-1]
+            )
+
+        v_on_t_points[:, :, 1:, :] = 0.5 * (
+            gridded_v.dataset.v_velocity.values[:, :, 1:, :] + gridded_v.dataset.v_velocity.values[:, :, :-1, :]
+            )
+
+        For arrays of shape and order: (z_dim, t_dim, y_dim, x_dim)
+
+        Parameters
+        ----------
+        gridded_u : TYPE
+            DESCRIPTION.
+        gridded_v : TYPE
+            DESCRIPTION.
+
+        Returns
+        -------
+        None.
+
+        """
+        ds_u = self.dataset.copy(deep=True)
+        # U velocity on T-points
+        ds_u["ut_velocity"] = 0.5 * (gridded_u.dataset.u_velocity.shift(x_dim=1) + gridded_u.dataset.u_velocity)
+        # replace wrapped (1st) longitude coord with zero
+        _, _lon = xr.broadcast(gridded_u.dataset.u_velocity, gridded_u.dataset.longitude)
+        ds_u["ut_velocity"] = ds_u["ut_velocity"].where(
+            _lon != _lon.isel(x_dim=0), 0
+        )  # keep values except where lon(x_dim=0)
+        del _, _lon
+        ds_u.coords["latitude"] = self.dataset.latitude
+        ds_u.coords["longitude"] = self.dataset.longitude
+        ds_u.coords["depth_0"] = self.dataset.depth_0
+        try:
+            self.dataset["ut_velocity"] = ds_u.ut_velocity.drop("depthu")
+        except:
+            self.dataset["ut_velocity"] = ds_u.ut_velocity
+            debug("Did not find depthu variable to drop - to avoid conflicts in z_dim dimension")
+
+        # V velocity on T-points
+        ds_v = self.dataset.copy(deep=True)
+
+        ds_v["vt_velocity"] = 0.5 * (gridded_v.dataset.v_velocity.shift(y_dim=1) + gridded_v.dataset.v_velocity)
+        # replace wrapped (1st) latitude coord with zero
+        _, _lat = xr.broadcast(gridded_v.dataset.v_velocity, gridded_v.dataset.latitude)
+        ds_v["vt_velocity"] = ds_v["vt_velocity"].where(
+            _lat != _lat.isel(y_dim=0), 0
+        )  # keep values except where lat(y_dim=0)
+        del _, _lat
+        ds_v.coords["latitude"] = self.dataset.latitude
+        ds_v.coords["longitude"] = self.dataset.longitude
+        ds_v.coords["depth_0"] = self.dataset.depth_0
+        try:
+            self.dataset["vt_velocity"] = ds_v.vt_velocity.drop("depthv")
+        except:
+            self.dataset["vt_velocity"] = ds_v.vt_velocity
+            debug("Did not find depthv variable to drop - to avoid conflicts in z_dim dimension")
+
+        # Speed on T-points
+        self.dataset["speed_t"] = np.sqrt(
+            np.square(self.dataset["ut_velocity"]) + np.square(self.dataset["vt_velocity"])
+        )
+
+    def quick_plot(self, name, Vmax=0.16, Np=3, headwidth=4, scale=50, time_value=None, **kwargs):
+        """
+        plot surface circulation
+        direction: unit vector
+        magnitude: shaded
+        """
+        # %%
+        from matplotlib import cm
+        from matplotlib.colors import LinearSegmentedColormap
+
+        nx = self.dataset.x_dim.size
+        ny = self.dataset.y_dim.size
+
+        if time_value == None:
+            SP = np.squeeze(self.dataset.speed_t)
+            US = np.squeeze(self.dataset.ut_velocity / SP)
+            VS = np.squeeze(self.dataset.vt_velocity / SP)
+        else:
+            SP = np.squeeze(self.dataset.speed_t[time_value, :, :])
+            US = np.squeeze(self.dataset.ut_velocity[time_value, :, :] / SP)
+            VS = np.squeeze(self.dataset.vt_velocity[time_value, :, :] / SP)
+        mask = self.dataset.bottom_level != 0
+        p = np.ma.masked_where(mask == 0, SP)
+        u = np.ma.masked_where(mask[0::Np, 0::Np] == 0, US[0::Np, 0::Np])
+        v = np.ma.masked_where(mask[0::Np, 0::Np] == 0, VS[0::Np, 0::Np])
+        # x=nemo_t.dataset.longitude[0::Np,0::Np]
+        # y=nemo_t.dataset.latitude[0::Np,0::Np]
+        X = self.dataset.longitude
+        Y = self.dataset.latitude
+        # create a light colour map
+        N_colours = int(Vmax * 100)
+        n_c = 2
+        cmap0 = cm.get_cmap("BrBG_r", lut=N_colours + n_c * 2)
+        colors = cmap0(np.arange(cmap0.N))
+        colors1 = colors[n_c : cmap0.N - n_c]
+        cmap1 = LinearSegmentedColormap.from_list("cmap1", colors1, cmap0.N - n_c * 2)
+
+        cmap1.set_bad([0.75, 0.75, 0.75])
+
+        fig = plt.figure()
+        fig.clear()
+        plt.pcolormesh(p, cmap=cmap1)
+        # plt.pcolormesh(X,Y,p,cmap=cmap1)
+        x = np.arange(0, nx, Np)
+        y = np.arange(0, ny, Np)
+
+        plt.clim([0, Vmax])
+        plt.colorbar(orientation="vertical", cmap=cmap1)
+        # plt.quiver(x,y,u,v,color=[0.1,0.1,0.1],headwidth=4,scale=50)
+        plt.quiver(x, y, u, v, color=[0.1, 0.1, 0.1], headwidth=headwidth, scale=scale)  # ,**kwargs)
+        plt.title("Surface Currents " + name)
+        return fig
diff --git a/coast/diagnostics/gridded_stratification.py b/coast/diagnostics/gridded_stratification.py
index 55130b65..d41855da 100644
--- a/coast/diagnostics/gridded_stratification.py
+++ b/coast/diagnostics/gridded_stratification.py
@@ -277,7 +277,7 @@ def calc_pea(self, gridded_t: xr.Dataset, Zd_mask):
         #          z_axis=1
 
         PEA = (z_4d * (rho - rhobar) * dz_4d).sum(dim="z_dim", skipna=True) * gravity / height
-        #%%
+        # %%
         #         return PEA
         coords = {
             "time": ("t_dim", gridded_t.dataset.time.values),
@@ -338,5 +338,5 @@ def quick_plot(self, var: xr.DataArray = None):
             plt.ylabel("latitude")
             plt.clim([0, 50])
             plt.colorbar()
-            plt.show()
+            # plt.show()
         return fig, ax
diff --git a/coast/diagnostics/profile_analysis.py b/coast/diagnostics/profile_analysis.py
index c815e1bc..b138d7e6 100644
--- a/coast/diagnostics/profile_analysis.py
+++ b/coast/diagnostics/profile_analysis.py
@@ -164,7 +164,7 @@ def determine_mask_indices(cls, profile, mask_dataset):
         Returns
         -------
         Dataset describing which profiles are in which mask/region.
-        Ready for input to Profile.mask_means()
+        Ready for input to Profile.mask_means() or Profile.mask_stats()
         """
         dataset = profile.dataset
         # If landmask not present, set to None for nearest_indices (no masking)
@@ -188,11 +188,98 @@ def determine_mask_indices(cls, profile, mask_dataset):
 
         return region_indices.rename({"dim_0": "id_dim"})
 
+    @classmethod
+    def mask_stats(cls, profile, mask_indices):
+        """
+        Mean and Standard deviation of all data inside a given profile dataset across a regional mask
+        or for multiples regional masks.
+
+        Parameters
+        ----------
+        dataset : xarray.Dataset
+            The profile dataset to average.
+        mask_indices : xarray.Dataset
+            Describes which profiles are in which region. Returned from
+            profile_analysis.determine_mask_indices().
+
+        Returns
+        -------
+        xarray.Dataset containing mean and std data.
+
+        """
+        dataset = profile.dataset
+        # Get the number of masks provided
+        n_masks = mask_indices.dims["dim_mask"]
+        mask_names = []
+
+        # Loop over masked arrays. Assign mean and std to mask and seasonal means
+        debug(f"Calculating masked mean and std..")
+
+        # Loop over masks and index the data
+        for mm in range(0, n_masks):
+            mask = mask_indices.isel(dim_mask=mm).mask.values
+            mask_ind = np.where(mask.astype(bool))[0]
+            if len(mask_ind) < 1:
+                if mm == 0:
+                    ds_stats = xr.Dataset()
+                continue
+            else:
+                mask_names.append(mask_indices["region_names"][mm].values)
+
+            # Get actual profile data for this mask
+            mask_data = dataset.isel(id_dim=mask_ind)
+            # Get two averages. One preserving depths and the other averaging across all data in a region
+            ds_mean_prof = mask_data.mean(dim="id_dim", skipna=True).compute()
+            ds_mean_all = mask_data.mean(skipna=True).compute()
+
+            # Loop over variables and save to output dataset.
+            var_list = list(ds_mean_prof.keys())
+            for vv in var_list:
+                ds_mean_prof = ds_mean_prof.rename({vv: "profile_mean_" + vv})
+                ds_mean_all = ds_mean_all.rename({vv: "all_mean_" + vv})
+
+            # Get two std. One preserving depths and the other collapsing across all data in a region
+            try:  # At time of writing skipna=True not working for np.datetime64 or np.timedelta64 values. So remove them
+                ds_std_prof = mask_data.std(dim="id_dim", skipna=True).compute()
+                ds_std_all = mask_data.std(skipna=True).compute()
+            except:  # remove np.datetime64 and np.timedelta64 variables
+                ds_tmp = mask_data
+                var_list = list(ds_tmp.keys())
+                for vv in var_list:
+                    if (type(ds_tmp[vv].values[0]) == np.datetime64) or (type(ds_tmp[vv].values[0]) == np.timedelta64):
+                        ds_tmp = ds_tmp.drop_vars(vv)
+                ds_std_prof = ds_tmp.std(dim="id_dim", skipna=True).compute()
+                ds_std_all = ds_tmp.std(skipna=True).compute()
+
+            # Loop over variables and save to output dataset.
+            var_list = list(ds_std_prof.keys())
+            for vv in var_list:
+                ds_std_prof = ds_std_prof.rename({vv: "profile_std_" + vv})
+                ds_std_all = ds_std_all.rename({vv: "all_std_" + vv})
+
+            # Merge profile means and all means into one dataset for output
+            ds_stats_tmp = xr.merge((ds_mean_prof, ds_mean_all, ds_std_prof, ds_std_all))
+
+            # If only one mask then just return the merged dataset
+            # else concatenate onto existing dataset.
+            if mm == 0:
+                ds_stats = ds_stats_tmp
+            else:
+                ds_stats = xr.concat((ds_stats, ds_stats_tmp), dim="dim_mask")
+
+        if mask_names != []:
+            ds_stats["region_names"] = (["dim_mask"], mask_names)
+        else:
+            ds_stats["region_names"] = (["dim_mask"], range(ds_stats.dims["dim_mask"]))
+
+        return ds_stats.set_coords(["region_names"])
+
     @classmethod
     def mask_means(cls, profile, mask_indices):
         """
         Averages all data inside a given profile dataset across a regional mask
         or for multiples regional masks.
+        THIS IS SUPERSEDED BY mask_stats() WHICH DOES MEANS AND STD. KEPT FOR LEGACY REASONS.
 
         Parameters
         ----------
diff --git a/coast/diagnostics/tidegauge_analysis.py b/coast/diagnostics/tidegauge_analysis.py
index 2fdb30ed..9222d1b7 100644
--- a/coast/diagnostics/tidegauge_analysis.py
+++ b/coast/diagnostics/tidegauge_analysis.py
@@ -99,7 +99,6 @@ def harmonic_analysis_utide(
 
         # Loop over ports
         for pp in range(0, n_port):
-
             # Temporary in-loop datasets
             ds_port = ds.isel(id_dim=pp).load()
             number_of_nan = np.sum(np.isnan(ds_port.values))
@@ -149,7 +148,6 @@ def reconstruct_tide_utide(cls, data_array, utide_solution_list, constit=None, o
 
         # Loop over ports
         for pp in np.arange(n_port):
-
             # Reconstruct full tidal signal using utide
             pp_solution = utide_solution_list[pp]
             if len(pp_solution) == 0:
@@ -162,6 +160,7 @@ def reconstruct_tide_utide(cls, data_array, utide_solution_list, constit=None, o
         # Create output dataset and return it in new Tidegauge object.
         ds_out = xr.Dataset(data_array.coords)
         ds_out[output_name] = (["id_dim", "t_dim"], reconstructed)
+        ds_out = ds_out.assign_coords(time=("t_dim", time.data))
 
         return Tidegauge(dataset=ds_out)
 
@@ -236,7 +235,6 @@ def threshold_statistics(cls, dataset, thresholds=np.arange(-0.4, 2, 0.1), peak_
 
         # Loop over vars in the input dataset
         for vv in var_list:
-
             empty_thresh = np.zeros((n_port, n_thresholds)) * np.nan
             ds_thresh["peak_count_" + vv] = (["id_dim", "threshold"], np.array(empty_thresh))
             ds_thresh["time_over_threshold_" + vv] = (["id_dim", "threshold"], np.array(empty_thresh))
@@ -244,7 +242,6 @@ def threshold_statistics(cls, dataset, thresholds=np.arange(-0.4, 2, 0.1), peak_
             ds_thresh["monthlymax_count_" + vv] = (["id_dim", "threshold"], np.array(empty_thresh))
 
             for pp in range(n_port):
-
                 # Identify NTR peaks for threshold analysis
                 data_pp = ds[vv].isel(id_dim=pp)
                 if np.sum(np.isnan(data_pp.values)) == ds.sizes["t_dim"]:
@@ -255,7 +252,6 @@ def threshold_statistics(cls, dataset, thresholds=np.arange(-0.4, 2, 0.1), peak_
 
                 # Threshold Analysis
                 for nn in range(0, n_thresholds):
-
                     # Calculate daily and monthly maxima for threshold analysis
                     ds_daily = data_pp.groupby("time.day")
                     ds_daily_max = ds_daily.max(skipna=True)
@@ -451,7 +447,6 @@ def crps(
 
         # Loop over location indices
         for ii in range(n_id):
-
             obs_ii = obs_var.isel(id_dim=ii)
 
             # Calculate CRPS
diff --git a/example_scripts/configuration_gallery/amm15_example_plot.py b/example_scripts/configuration_gallery/amm15_example_plot.py
index a87f1759..b3151f89 100755
--- a/example_scripts/configuration_gallery/amm15_example_plot.py
+++ b/example_scripts/configuration_gallery/amm15_example_plot.py
@@ -36,7 +36,7 @@
 print("* Loaded ", config, " data")
 
 #################################################
-#%% subset of data and domain ##
+# %% subset of data and domain ##
 #################################################
 # Pick out a North Sea subdomain
 print("* Extract North Sea subdomain")
@@ -45,7 +45,7 @@
 ind_sci = sci_w.subset_indices(start=[51, -4], end=[62, 15])
 sci_nwes_w = sci_w.isel(y_dim=ind_sci[0], x_dim=ind_sci[1])  # nwes = northwest europe shelf
 
-#%% Apply masks to temperature and salinity
+# %% Apply masks to temperature and salinity
 if config == "AMM15":
     sci_nwes_t.dataset["temperature_m"] = sci_nwes_t.dataset.temperature.where(
         sci_nwes_t.dataset.mask.expand_dims(dim=sci_nwes_t.dataset["t_dim"].sizes) > 0
@@ -60,7 +60,7 @@
     sci_nwes_t.dataset["salinity_m"] = sci_nwes_t.dataset.salinity
 
 
-#%% Plots
+# %% Plots
 fig = plt.figure()
 
 plt.pcolormesh(sci_t.dataset.longitude, sci_t.dataset.latitude, sci_t.dataset.temperature.isel(z_dim=0).squeeze())
diff --git a/example_scripts/configuration_gallery/anchor_plots_of_nsea_wvel.py b/example_scripts/configuration_gallery/anchor_plots_of_nsea_wvel.py
index f412fc0f..a3e705e3 100755
--- a/example_scripts/configuration_gallery/anchor_plots_of_nsea_wvel.py
+++ b/example_scripts/configuration_gallery/anchor_plots_of_nsea_wvel.py
@@ -5,14 +5,14 @@
 This needs to move to the above
 """
 
-#%%
+# %%
 import coast
 import matplotlib.pyplot as plt
 
 # import matplotlib.colors as colors # colormap fiddling
 
 #################################################
-#%%  Loading and initialising methods ##
+# %%  Loading and initialising methods ##
 #################################################
 
 dir_nam = "/projectsa/anchor/NEMO/AMM60/"
@@ -32,18 +32,18 @@
 sci_w = coast.Gridded(dir_nam + fil_nam, dom_nam, config=config)
 sci_w.dataset.chunk(chunks)
 
-#% NEMO output is not standard with u,v fields included with w-pts. Tidy to avoid confusion
+# % NEMO output is not standard with u,v fields included with w-pts. Tidy to avoid confusion
 sci_w.dataset = sci_w.dataset.drop_vars(["uo", "vo", "depthv"])
 sci_w.dataset = sci_w.dataset.swap_dims({"depthw": "z_dim"})
 
 #################################################
-#%% subset of data and domain ##
+# %% subset of data and domain ##
 #################################################
 # Pick out a North Sea subdomain
 ind_sci = sci_w.subset_indices(start=[51, -4], end=[60, 15])
 sci_nwes_w = sci_w.isel(y_dim=ind_sci[0], x_dim=ind_sci[1])  # nswes = northwest europe shelf
 
-#%% Compute a diffusion from w-vel
+# %% Compute a diffusion from w-vel
 Kz = (sci_nwes_w.dataset.wo * sci_nwes_w.dataset.e3_0).sum(dim="z_dim").mean(dim="t_dim")
 
 # plot map
@@ -61,7 +61,7 @@
 plt.colorbar()
 # fig.savefig("")
 
-#%% Transect Method
+# %% Transect Method
 tran_w = coast.TransectT(sci_nwes_w, (51, 2.5), (61, 2.5))
 
 lat_sec = tran_w.data.latitude.expand_dims(dim={"z_dim": 51})
@@ -70,7 +70,7 @@
 # wo_sec = tran.data_F.wo.mean(dim='t_dim')
 
 
-#%% Make map and profile plots
+# %% Make map and profile plots
 #################################################
 for i in range(2):
     for lat0 in [54, 57]:
@@ -132,7 +132,7 @@
     plt.close("all")
 
 
-#%% Plot sections
+# %% Plot sections
 fig = plt.figure()
 plt.rcParams["figure.figsize"] = 8, 8
 
diff --git a/example_scripts/configuration_gallery/blz_example_plot.py b/example_scripts/configuration_gallery/blz_example_plot.py
index 6aa0b642..07ca6186 100755
--- a/example_scripts/configuration_gallery/blz_example_plot.py
+++ b/example_scripts/configuration_gallery/blz_example_plot.py
@@ -5,12 +5,12 @@
 
 """
 
-#%%
+# %%
 import coast
 import matplotlib.pyplot as plt
 
 #################################################
-#%%  Loading  data
+# %%  Loading  data
 #################################################
 
 
@@ -34,6 +34,6 @@
 sci_w = coast.Gridded(fn_domain=dom_nam, config=config_w)
 
 
-#%% Plot
+# %% Plot
 plt.pcolormesh(sci_t.dataset.ssh.isel(t_dim=0))
 plt.show()
diff --git a/example_scripts/configuration_gallery/seasia_dic_example_plot.py b/example_scripts/configuration_gallery/seasia_dic_example_plot.py
index 469290b6..a2d72f37 100644
--- a/example_scripts/configuration_gallery/seasia_dic_example_plot.py
+++ b/example_scripts/configuration_gallery/seasia_dic_example_plot.py
@@ -4,13 +4,13 @@
 Make simple SEAsia 1/12 deg DIC plot.
 
 """
-#%%
+# %%
 import coast
 import matplotlib.pyplot as plt
 
 
 #################################################
-#%%  Loading  data
+# %%  Loading  data
 #################################################
 path_examples = "./example_files/"
 ## data local in livljobs : /projectsa/COAsT/NEMO_example_data/SEAsia_R12/
@@ -22,7 +22,7 @@
 
 seasia_bgc = coast.Gridded(fn_data=fn_seasia_var, fn_domain=fn_seasia_domain, config=fn_seasia_config_bgc)
 
-#%% Plot DIC
+# %% Plot DIC
 fig = plt.figure()
 plt.pcolormesh(
     seasia_bgc.dataset.longitude,
diff --git a/example_scripts/configuration_gallery/seasia_r12_example_plot.py b/example_scripts/configuration_gallery/seasia_r12_example_plot.py
index d68495d5..2792f34c 100644
--- a/example_scripts/configuration_gallery/seasia_r12_example_plot.py
+++ b/example_scripts/configuration_gallery/seasia_r12_example_plot.py
@@ -5,13 +5,13 @@
 
 """
 
-#%%
+# %%
 import coast
 import matplotlib.pyplot as plt
 
 
 #################################################
-#%%  Loading  data
+# %%  Loading  data
 #################################################
 
 dir_nam = "/projectsa/COAsT/NEMO_example_data/SEAsia_R12/"
@@ -21,7 +21,7 @@
 
 sci_t = coast.Gridded(dir_nam + fil_nam, dir_nam + dom_nam, config=config_t)
 
-#%% Plot
+# %% Plot
 fig = plt.figure()
 
 plt.pcolormesh(sci_t.dataset.longitude, sci_t.dataset.latitude, sci_t.dataset.salinity.isel(t_dim=0).isel(z_dim=0))
diff --git a/example_scripts/configuration_gallery/wcssp_india_example_plot.py b/example_scripts/configuration_gallery/wcssp_india_example_plot.py
index bbc712a6..63b3fcb8 100644
--- a/example_scripts/configuration_gallery/wcssp_india_example_plot.py
+++ b/example_scripts/configuration_gallery/wcssp_india_example_plot.py
@@ -7,12 +7,12 @@
 Simple plot of sea surface temperature
 """
 
-#%%
+# %%
 import coast
 import matplotlib.pyplot as plt
 
 #################################################
-#%%  Loading  data
+# %%  Loading  data
 #################################################
 
 dir_nam = "/projectsa/COAsT/NEMO_example_data/MO_INDIA/"
@@ -22,7 +22,7 @@
 
 sci_t = coast.Gridded(dir_nam + fil_nam, dir_nam + dom_nam, config=config_t)
 
-#%% Plot
+# %% Plot
 fig = plt.figure()
 
 plt.pcolormesh(sci_t.dataset.longitude, sci_t.dataset.latitude, sci_t.dataset.temperature.isel(t_dim=0).isel(z_dim=0))
diff --git a/example_scripts/notebooks/runnable_notebooks/general/mask_maker_tutorial.ipynb b/example_scripts/notebooks/runnable_notebooks/general/mask_maker_tutorial.ipynb
index 504bca24..a1dec8b0 100644
--- a/example_scripts/notebooks/runnable_notebooks/general/mask_maker_tutorial.ipynb
+++ b/example_scripts/notebooks/runnable_notebooks/general/mask_maker_tutorial.ipynb
@@ -172,18 +172,19 @@
     "masks_list.append(np.ones(lon.shape))\n",
     "\n",
     "# Add regional mask for English Channel\n",
-    "masks_list.append(mm.region_def_nws_north_sea(lon, lat, bathy))\n",
+    "masks_list.append(mm.region_def_nws_north_north_sea(lon, lat, bathy))\n",
+    "masks_list.append(mm.region_def_nws_south_north_sea(lon, lat, bathy))\n",
     "masks_list.append(mm.region_def_nws_outer_shelf(lon, lat, bathy))\n",
     "masks_list.append(mm.region_def_nws_norwegian_trench(lon, lat, bathy))\n",
     "masks_list.append(mm.region_def_nws_english_channel(lon, lat, bathy))\n",
-    "masks_list.append(mm.region_def_south_north_sea(lon, lat, bathy))\n",
-    "masks_list.append(mm.region_def_off_shelf(lon, lat, bathy))\n",
-    "masks_list.append(mm.region_def_irish_sea(lon, lat, bathy))\n",
-    "masks_list.append(mm.region_def_kattegat(lon, lat, bathy))\n",
+    "masks_list.append(mm.region_def_nws_off_shelf(lon, lat, bathy))\n",
+    "masks_list.append(mm.region_def_nws_irish_sea(lon, lat, bathy))\n",
+    "masks_list.append(mm.region_def_nws_kattegat(lon, lat, bathy))\n",
+    "masks_list.append(mm.region_def_nws_fsc(lon, lat, bathy))\n",
     "\n",
-    "masks_names = [\"whole domain\", \"north sea\", \"outer shelf\", \"norwegian trench\",\n",
-    "                \"english_channel\", \"southern north sea\", \"off shelf\",\n",
-    "                \"irish sea\", \"kattegat\",]"
+    "masks_names = [\"whole domain\", \"northern north sea\", \"southern north sea\", \"outer shelf\",\n",
+    "               \"norwegian trench\", \"english_channel\",  \"off shelf\",\n",
+    "                \"irish sea\", \"kattegat\", \"fsc\"]"
    ]
   },
   {
@@ -243,7 +244,7 @@
    "id": "454bdeb4",
    "metadata": {},
    "source": [
-    "NB overlapping regions are not given special treatment, the layers are blindly superimposed on each other. E.g. as demonstrated with \"Norwegian Trench\" and \"off shelf\", or \"whole domain\" and any other region."
+    "NB overlapping regions are not given special treatment, the layers are blindly superimposed on each other. E.g. as demonstrated with \"Norwegian Trench\" and \"off shelf\", \"FSC\" and \"off shelf\", or \"whole domain\" and any other region."
    ]
   },
   {
@@ -254,7 +255,7 @@
    "outputs": [],
    "source": [
     "plt.subplot(2,2,1)\n",
-    "mm.quick_plot(mask_xr.sel(dim_mask=[0,3]))\n",
+    "mm.quick_plot(mask_xr.sel(dim_mask=[0,3,9]))\n",
     "\n",
     "plt.subplot(2,2,2)\n",
     "mm.quick_plot(mask_xr.sel(dim_mask=[1,2,4,5,6,7,8]))\n",
@@ -360,7 +361,7 @@
    "id": "4e70c0df",
    "metadata": {},
    "source": [
-    "Notices that the number of mask dimensions is not necessarily preserved between the mask and the mask averaged variables. This happens if, for example, there are no profiles in one of the mask regions"
+    "Notice that the number of mask dimensions is not necessarily preserved between the mask and the mask averaged variables. This happens if, for example, there are no profiles in one of the mask regions"
    ]
   },
   {
@@ -379,7 +380,7 @@
    "id": "5a033719",
    "metadata": {},
    "source": [
-    "The mean profiles can be visualised or further processed (notice the Irish Sea region is missing because there were no profiles in the example dataset)"
+    "The mean profiles can be visualised or further processed (notice the Irish Sea region and FSC are missing because there were no profiles in the example dataset)"
    ]
   },
   {
@@ -483,7 +484,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.10"
+   "version": "3.8.15"
   }
  },
  "nbformat": 4,
diff --git a/example_scripts/profile_validation/plot_validation_gridded_data.py b/example_scripts/profile_validation/plot_validation_gridded_data.py
index 9c18aa82..4bcc2176 100644
--- a/example_scripts/profile_validation/plot_validation_gridded_data.py
+++ b/example_scripts/profile_validation/plot_validation_gridded_data.py
@@ -18,7 +18,7 @@
 import coast
 import pandas as pd
 
-#%% File settings
+# %% File settings
 run_name = "test"
 
 # List of analysis output files. Profiles from each will be plotted
@@ -31,7 +31,7 @@
 # Filename for the output
 fn_out = "/Users/dbyrne/transfer/surface_gridded_errors_{0}.png".format(run_name)
 
-#%% General Plot Settings
+# %% General Plot Settings
 var_name = "abs_diff_temperature"  # Variable name in analysis file to plot
 # If you used var modified to make gridded data
 # then this is where to select season etc.
@@ -72,7 +72,7 @@
 plot_seasons = True
 season_suffixes = ["DJF", "MAM", "JJA", "SON"]
 
-#%% Read and plotdata
+# %% Read and plotdata
 
 # Read all datasets into list
 ds_list = [xr.open_dataset(dd) for dd in fn_list]
diff --git a/example_scripts/profile_validation/plot_validation_mask_means.py b/example_scripts/profile_validation/plot_validation_mask_means.py
index 09eaf9d6..18b26098 100644
--- a/example_scripts/profile_validation/plot_validation_mask_means.py
+++ b/example_scripts/profile_validation/plot_validation_mask_means.py
@@ -15,7 +15,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-#%% File settings
+# %% File settings
 run_name = "test"
 
 # List of analysis output files. Profiles from each will be plotted
@@ -25,7 +25,7 @@
 # Filename for the output
 fn_out = "/Users/dbyrne/transfer/regional_means_{0}.png".format(run_name)
 
-#%% General Plot Settings
+# %% General Plot Settings
 region_ind = [0, 1, 2, 3, 4, 5, 6, 7, 8]  # Region indices (in analysis) to plot
 region_names = ["A", "B", "C", "D", "E", "F", "G", "H", "I"]  # Region names, will be used for titles in plot
 var_name = "profile_average_abs_diff_temperature"  # Variable name in analysis file to plot
@@ -64,7 +64,7 @@
 title_fontweight = "bold"  # Fontweight to use for title
 
 
-#%% SCRIPT: READ AND PLOT DATA
+# %% SCRIPT: READ AND PLOT DATA
 
 # Read all datasets into list
 ds_list = [xr.open_dataset(dd) for dd in fn_list]
@@ -78,7 +78,6 @@
 
 # Loop over regions
 for ii in range(n_ax):
-
     if ii >= n_reg:
         a_flat[ii].axis("off")
         continue
diff --git a/example_scripts/profile_validation/plot_validation_surface_errors.py b/example_scripts/profile_validation/plot_validation_surface_errors.py
index ea16f9c7..fc4f45e7 100644
--- a/example_scripts/profile_validation/plot_validation_surface_errors.py
+++ b/example_scripts/profile_validation/plot_validation_surface_errors.py
@@ -18,7 +18,7 @@
 import coast
 import pandas as pd
 
-#%% File settings
+# %% File settings
 run_name = "test"
 
 # List of analysis output files. Profiles from each will be plotted
@@ -31,7 +31,7 @@
 # Filename for the output
 fn_out = "/Users/dbyrne/transfer/surface_errors_{0}.png".format(run_name)
 
-#%% General Plot Settings
+# %% General Plot Settings
 var_name = "diff_temperature"  # Variable name in analysis file to plot
 save_plot = False
 
@@ -64,7 +64,7 @@
 season_str = "DJF"  # DJF, MAM, JJA or SON
 
 
-#%% Read and plotdata
+# %% Read and plotdata
 
 # Read all datasets into list
 ds_list = [xr.open_dataset(dd)[var_name] for dd in fn_list]
diff --git a/example_scripts/stratification_pycnocline_diagnostics.py b/example_scripts/stratification_pycnocline_diagnostics.py
index e30cdd16..b69fcae9 100755
--- a/example_scripts/stratification_pycnocline_diagnostics.py
+++ b/example_scripts/stratification_pycnocline_diagnostics.py
@@ -8,7 +8,7 @@
 
 """
 
-#%%
+# %%
 import coast
 import numpy as np
 import os
@@ -16,7 +16,7 @@
 import matplotlib.colors as colors  # colormap fiddling
 
 #################################################
-#%%  Loading  data
+# %%  Loading  data
 #################################################
 
 #  Loading AMM60 data if it is available
@@ -71,7 +71,7 @@
 print("* Loaded ", config, " data")
 
 #################################################
-#%% subset of data and domain ##
+# %% subset of data and domain ##
 #################################################
 # Pick out a North Sea subdomain
 print("* Extract North Sea subdomain")
@@ -80,7 +80,7 @@
 ind_sci = sci_w.subset_indices(start=[51, -4], end=[62, 15])
 sci_nwes_w = sci_w.isel(y_dim=ind_sci[0], x_dim=ind_sci[1])  # nwes = northwest europe shelf
 
-#%% Apply masks to temperature and salinity
+# %% Apply masks to temperature and salinity
 if config == "AMM60":
     sci_nwes_t.dataset["temperature_m"] = sci_nwes_t.dataset.temperature.where(
         sci_nwes_t.dataset.mask.expand_dims(dim=sci_nwes_t.dataset["t_dim"].sizes) > 0
@@ -95,32 +95,32 @@
     sci_nwes_t.dataset["salinity_m"] = sci_nwes_t.dataset.salinity
 
 
-#%% Construct in-situ density and stratification
+# %% Construct in-situ density and stratification
 print("* Construct in-situ density and stratification")
 sci_nwes_t.construct_density(eos="EOS10")
 
-#%% Construct stratification. t-pts --> w-pts
+# %% Construct stratification. t-pts --> w-pts
 print("* Construct stratification. t-pts --> w-pts")
 sci_nwes_w = sci_nwes_t.differentiate(
     "density", dim="z_dim", out_var_str="rho_dz", out_obj=sci_nwes_w
 )  # --> sci_nwes_w.rho_dz
 
 #################################################
-#%% Create internal tide diagnostics object
+# %% Create internal tide diagnostics object
 print("* Create stratification diagnostics object")
 strat = coast.GriddedStratification(sci_nwes_t, sci_nwes_w)
 
-#%%  Construct pycnocline variables: depth and thickness
+# %%  Construct pycnocline variables: depth and thickness
 print("* Compute density and rho_dz if they didn" "t exist")
 print("* Compute 1st and 2nd moments of stratification as pycnocline vars")
 strat.construct_pycnocline_vars(sci_nwes_t, sci_nwes_w)
 
-#%%  Plot pycnocline variables: depth and thickness
+# %%  Plot pycnocline variables: depth and thickness
 print("* Sample quick plot")
 strat.quick_plot()
 
 
-#%% Make transects
+# %% Make transects
 print("* Construct transects to inspect stratification. This is an abuse of the transect code...")
 # Example usage: tran = coast.Transect( (54,-15), (56,-12), nemo_f, nemo_t, nemo_u, nemo_v )
 tran_it = coast.TransectT(strat, (51, 2.5), (61, 2.5))
@@ -143,7 +143,7 @@
 zt_m_sec = tran_it.data.strat_2nd_mom_masked.mean(dim="t_dim", skipna=False)
 
 
-#%% Plot sections
+# %% Plot sections
 #################
 print("* Plot sections with pycnocline depth and thickness overlayed")
 plt.pcolormesh(lat_sec, dep_sec, rho_sec)
@@ -174,7 +174,7 @@
 plt.show()
 
 
-#%% Plot profile of density and stratification with strat_1st_mom in deep water
+# %% Plot profile of density and stratification with strat_1st_mom in deep water
 #############################################################################
 print("* Plot profile of density and stratification with strat_1st_mom in deep water")
 print(
@@ -213,7 +213,7 @@
 plt.show()
 
 
-#%% Map pretty plots of North Sea pycnocline depth
+# %% Map pretty plots of North Sea pycnocline depth
 print("* Map pretty plots of North Sea pycnocline depth")
 print(" - we expect a RunTimeError here")
 
diff --git a/setup.py b/setup.py
index 70a8fcdb..f56dffaa 100644
--- a/setup.py
+++ b/setup.py
@@ -5,7 +5,7 @@
 PACKAGE = SimpleNamespace(
     **{
         "name": "COAsT",
-        "version": "3.2.0",
+        "version": "3.2.1",
         "description": "This is the Coast Ocean Assessment Tool",
         "url": "https://www.bodc.ac.uk",
         "download_url": "https://github.com/British-Oceanographic-Data-Centre/COAsT/",


Test patch:
diff --git a/tests/test_tidetable.py b/tests/test_tidetable.py
index ce81a623..eb643b74 100644
--- a/tests/test_tidetable.py
+++ b/tests/test_tidetable.py
@@ -9,7 +9,7 @@
 
 
 # -----------------------------------------------------------------------------#
-#%% day of the week function                                           #
+# %% day of the week function                                           #
 #                                                                             #
 def test_dayoweek():
     check1 = general_utils.day_of_week(np.datetime64("2020-10-16")) == "Fri"
diff --git a/unit_testing/generate_unit_test_contents.py b/unit_testing/generate_unit_test_contents.py
index b9d835ce..77868aaa 100644
--- a/unit_testing/generate_unit_test_contents.py
+++ b/unit_testing/generate_unit_test_contents.py
@@ -74,7 +74,6 @@
 
 # Open output file
 with open(fn_contents, "w") as file:
-
     # Write title things
     file.write("     UNIT TEST CONTENTS FILE TEST    \n")
     file.write("\n")
diff --git a/unit_testing/test_TEMPLATE.py b/unit_testing/test_TEMPLATE.py
index 93f417f5..74a8b823 100644
--- a/unit_testing/test_TEMPLATE.py
+++ b/unit_testing/test_TEMPLATE.py
@@ -18,6 +18,7 @@
 # IMPORT THIS TO HAVE ACCESS TO EXAMPLE FILE PATHS:
 import unit_test_files as files
 
+
 # Define a testing class. Absolutely fine to have one or multiple per file.
 # Each class must inherit unittest.TestCase
 class test_TEMPLATE(unittest.TestCase):
diff --git a/unit_testing/test_altimetry_methods.py b/unit_testing/test_altimetry_methods.py
index f82b4cca..b8eff883 100644
--- a/unit_testing/test_altimetry_methods.py
+++ b/unit_testing/test_altimetry_methods.py
@@ -15,7 +15,6 @@
 
 class test_altimetry_methods(unittest.TestCase):
     def test_altimetry_load_subset_and_comparison(self):
-
         sci = coast.Gridded(files.fn_nemo_dat, files.fn_nemo_dom, config=files.fn_config_t_grid)
 
         with self.subTest("Load example altimetry file"):
diff --git a/unit_testing/test_bgc_gridded_initialisation.py b/unit_testing/test_bgc_gridded_initialisation.py
index 8aab94e8..37fde941 100644
--- a/unit_testing/test_bgc_gridded_initialisation.py
+++ b/unit_testing/test_bgc_gridded_initialisation.py
@@ -12,6 +12,7 @@
 # IMPORT THIS TO HAVE ACCESS TO EXAMPLE FILE PATHS:
 import unit_test_files as files
 
+
 # Define a testing class. Absolutely fine to have one or multiple per file.
 # Each class must inherit unittest.TestCase
 class test_bgc_gridded_initialisation(unittest.TestCase):
diff --git a/unit_testing/test_diagnostic_methods.py b/unit_testing/test_diagnostic_methods.py
index f054d90f..b34db670 100644
--- a/unit_testing/test_diagnostic_methods.py
+++ b/unit_testing/test_diagnostic_methods.py
@@ -130,3 +130,43 @@ def test_construct_pycnocline_depth_and_thickness(self):
             fig.tight_layout()
             fig.savefig(files.dn_fig + "strat_1st_mom.png")
             plt.close("all")
+
+    def test_circulation(self):
+        # %%
+        lims = [150, 250, 100, 350]
+        nemo_u = coast.Gridded(
+            fn_data=files.fn_nemo_grid_u_dat, fn_domain=files.fn_nemo_dom, config=files.fn_config_u_grid, lims=lims
+        )
+        nemo_v = coast.Gridded(
+            fn_data=files.fn_nemo_grid_v_dat, fn_domain=files.fn_nemo_dom, config=files.fn_config_v_grid, lims=lims
+        )
+
+        with self.subTest("Map velocities to t-points"):
+            nemo_t = coast.CurrentsOnT(fn_domain=files.fn_nemo_dom, config=files.fn_config_t_grid, lims=lims)
+            nemo_t.currents_on_t(nemo_u, nemo_v)
+            nemo_t.subset(z_dim=[0], t_dim=[0])
+
+            u1 = 0.5 * (
+                nemo_u.dataset.u_velocity[0, 0, 150, 60].values + nemo_u.dataset.u_velocity[0, 0, 150, 59].values
+            )
+            u2 = nemo_t.dataset.ut_velocity[0, 0, 150, 60].values
+            # print(f"u vel on u-pts: {nemo_u.dataset.u_velocity.isel(x_dim=slice(59,61), y_dim=slice(149,151), t_dim=0, z_dim=0).values}")
+            # print(f"u vel on t-pts: {nemo_t.dataset.ut_velocity.isel(x_dim=slice(59,61), y_dim=slice(149,151), t_dim=0, z_dim=0).values}")
+            v1 = 0.5 * (
+                nemo_v.dataset.v_velocity[0, 0, 150, 60].values + nemo_v.dataset.v_velocity[0, 0, 149, 60].values
+            )
+            v2 = nemo_t.dataset.vt_velocity[0, 0, 150, 60].values
+            speed = nemo_t.dataset.speed_t[0, 0, 150, 60].values
+
+            check1 = np.isclose(u1, u2)
+            check2 = np.isclose(v1, v2)
+            check3 = np.isclose(speed, 0.04988122)
+            self.assertTrue(check1, msg="check 1")
+            self.assertTrue(check2, msg="check 2")
+            self.assertTrue(check3, msg="check 3")
+
+        with self.subTest("Plot circulation on t-points"):
+            fig = nemo_t.quick_plot("Test")
+            fig.tight_layout()
+            fig.savefig(files.dn_fig + "surface_circulation.png")
+            plt.close("all")
diff --git a/unit_testing/test_gridded_harmonics.py b/unit_testing/test_gridded_harmonics.py
index b1390072..92acc0da 100644
--- a/unit_testing/test_gridded_harmonics.py
+++ b/unit_testing/test_gridded_harmonics.py
@@ -11,7 +11,6 @@
 
 class test_gridded_harmonics(unittest.TestCase):
     def test_combine_and_convert_harmonics(self):
-
         harmonics = coast.Gridded(files.fn_nemo_harmonics, files.fn_nemo_harmonics_dom, config=files.fn_config_t_grid)
 
         with self.subTest("test_combine_harmonics"):
@@ -26,7 +25,6 @@ def test_combine_and_convert_harmonics(self):
             self.assertTrue(check2, msg="check2")
 
         with self.subTest("test_convert_harmonics"):
-
             harmonics_combined.harmonics_convert(direction="cart2polar")
             harmonics_combined.harmonics_convert(direction="polar2cart", x_var="x_test", y_var="y_test")
 
diff --git a/unit_testing/test_gridded_initialisation.py b/unit_testing/test_gridded_initialisation.py
index e7d27bca..0c4c57f8 100644
--- a/unit_testing/test_gridded_initialisation.py
+++ b/unit_testing/test_gridded_initialisation.py
@@ -150,3 +150,16 @@ def test_gridded_compute_e3_from_ssh(self):
         reference = np.array([8.337016e08, 8.333972e08, 8.344886e08, 8.330722e08, 8.265948e08])
         check1 = np.allclose(cksum, reference)
         self.assertTrue(check1, msg="check1")
+
+    def test_gridded_subset(self):
+        lims = [150, 250, 100, 350]
+        nemo_t_subset = coast.Gridded(
+            fn_data=files.fn_nemo_dat, fn_domain=files.fn_nemo_dom, config=files.fn_config_t_grid, lims=lims
+        )
+        nemo_t = coast.Gridded(fn_data=files.fn_nemo_dat, fn_domain=files.fn_nemo_dom, config=files.fn_config_t_grid)
+
+        T1 = nemo_t.dataset.temperature[0, 0, 105, 155].values
+        T2 = nemo_t_subset.dataset.temperature[0, 0, 5, 5].values
+        check1 = np.isclose(T1, T2)
+
+        self.assertTrue(check1, msg="check1")
diff --git a/unit_testing/test_isobath_contour_methods.py b/unit_testing/test_isobath_contour_methods.py
index df5f6523..a4f4290d 100644
--- a/unit_testing/test_isobath_contour_methods.py
+++ b/unit_testing/test_isobath_contour_methods.py
@@ -13,7 +13,6 @@
 
 class test_contour_f_methods(unittest.TestCase):
     def test_extract_isobath_contour_between_two_points(self):
-
         with self.subTest("Extract contour"):
             nemo_f = coast.Gridded(
                 fn_domain=files.fn_nemo_dom, config=files.fn_config_f_grid, calculate_bathymetry=False
diff --git a/unit_testing/test_object_manipulation.py b/unit_testing/test_object_manipulation.py
index 6b9d0f54..3ceac248 100644
--- a/unit_testing/test_object_manipulation.py
+++ b/unit_testing/test_object_manipulation.py
@@ -113,7 +113,6 @@ def test_indices_by_distance(self):
         self.assertTrue(check1, "check1")
 
     def test_interpolation_to_altimetry(self):
-
         sci = coast.Gridded(files.fn_nemo_dat, files.fn_nemo_dom, config=files.fn_config_t_grid)
 
         with self.subTest("Find nearest xy indices"):
diff --git a/unit_testing/test_profile_methods.py b/unit_testing/test_profile_methods.py
index 06080790..27644777 100644
--- a/unit_testing/test_profile_methods.py
+++ b/unit_testing/test_profile_methods.py
@@ -13,7 +13,6 @@
 
 class test_profile_methods(unittest.TestCase):
     def test_load_process_and_compare_profile_data(self):
-
         with self.subTest("Load profile data from EN4"):
             profile = coast.Profile(config=files.fn_profile_config)
             profile.read_en4(files.fn_profile)
@@ -38,7 +37,6 @@ def test_load_process_and_compare_profile_data(self):
             self.assertTrue(check3, "check3")
 
     def test_compare_processed_profile_with_model(self):
-
         profile = coast.Profile(config=files.fn_profile_config)
         profile.read_en4(files.fn_profile)
         profile.dataset = profile.dataset.isel(id_dim=np.arange(0, profile.dataset.dims["id_dim"], 10)).load()
@@ -102,7 +100,7 @@ def test_compare_processed_profile_with_model(self):
             lon = nemo_t.dataset.longitude.values
             lat = nemo_t.dataset.latitude.values
 
-            mm_north_sea = mm.region_def_nws_north_sea(lon, lat, bath)
+            mm_north_sea = mm.region_def_nws_north_north_sea(lon, lat, bath)
             mm_whole_domain = np.ones(lon.shape)
             mask_list = [mm_north_sea, mm_whole_domain]
             mask_names = ["North Sea", "Whole Domain"]
@@ -116,11 +114,25 @@ def test_compare_processed_profile_with_model(self):
             # Do average differences for each region
             mask_means = pa.mask_means(difference, mask_indices)
 
+            # Do mean and std for each region
+            mask_stats_model = pa.mask_stats(model_interpolated, mask_indices)
+
             check1 = mask_indices.dims["dim_mask"] == mask_means.dims["dim_mask"]  # since profiles are in mask regions
             check2 = np.isclose(mask_means.all_mean_diff_temperature.values[0], -0.8169331445866469)
 
+            check3 = (
+                mask_indices.dims["dim_mask"] == mask_stats_model.dims["dim_mask"]
+            )  # since profiles are in mask regions
+            check4 = np.isclose(
+                mask_stats_model.all_mean_temperature.values[0], 7.995568407463857
+            )  # includes mean and std
+            check5 = np.isclose(mask_stats_model.all_std_temperature.values[0], 0.8448488052959595)
+
             self.assertTrue(check1, "check1")
             self.assertTrue(check2, "check2")
+            self.assertTrue(check3, "check3")
+            self.assertTrue(check4, "check4")
+            self.assertTrue(check5, "check5")
 
         with self.subTest("Surface/Bottom averaging"):
             surface = 5
diff --git a/unit_testing/test_tidegauge_methods.py b/unit_testing/test_tidegauge_methods.py
index 39e39bae..0a5466a6 100644
--- a/unit_testing/test_tidegauge_methods.py
+++ b/unit_testing/test_tidegauge_methods.py
@@ -51,7 +51,6 @@ def test_demean_timeseries(self):
         self.assertTrue(np.array_equal(tg1.dataset.ssh[0, :].values, demeaned), "check1")
 
     def test_harmonic_analysis_utide(self):
-
         tganalysis = coast.TidegaugeAnalysis()
         date0 = datetime.datetime(2007, 1, 10)
         date1 = datetime.datetime(2007, 1, 12)
@@ -134,7 +133,6 @@ def test_read_gesla_formats(self):
         self.assertTrue(check2, "check2")
 
     def test_read_gesla_and_compare_to_model(self):
-
         sci = coast.Gridded(files.fn_nemo_dat, files.fn_nemo_dom, config=files.fn_config_t_grid)
         sci.dataset["landmask"] = sci.dataset.bottom_level == 0
         tganalysis = coast.TidegaugeAnalysis()
@@ -210,7 +208,6 @@ def test_time_slice(self):
         self.assertTrue(check2, "check2")
 
     def test_tidegauge_resample_and_apply_doodsonx0(self):
-
         with self.subTest("Resample to 1H"):
             tganalysis = coast.TidegaugeAnalysis()
             date0 = datetime.datetime(2007, 1, 10)
@@ -235,7 +232,6 @@ def test_tidegauge_resample_and_apply_doodsonx0(self):
             self.assertTrue(check2, "check2")
 
     def test_load_multiple_tidegauge(self):
-
         with self.subTest("Load multiple gauge"):
             date0 = datetime.datetime(2007, 1, 10)
             date1 = datetime.datetime(2007, 1, 12)
@@ -286,7 +282,6 @@ def test_tidegauge_for_tabulated_data(self):
         self.assertTrue(check5, "check5")
 
     def test_tidegauge_finding_extrema(self):
-
         with self.subTest("Find extrema"):
             date0 = datetime.datetime(2007, 1, 10)
             date1 = datetime.datetime(2007, 1, 20)
@@ -321,7 +316,6 @@ def test_tidegauge_finding_extrema(self):
             plt.close("all")
 
     def test_tidegauge_cubic_spline_extrema(self):
-
         with self.subTest("Fit cubic spline"):
             date_start = np.datetime64("2020-10-12 23:59")
             date_end = np.datetime64("2020-10-14 00:01")
diff --git a/unit_testing/test_wod_read_data.py b/unit_testing/test_wod_read_data.py
index 5f0a434a..5d009782 100644
--- a/unit_testing/test_wod_read_data.py
+++ b/unit_testing/test_wod_read_data.py
@@ -11,13 +11,12 @@
 # IMPORT THIS TO HAVE ACCESS TO EXAMPLE FILE PATHS:
 import unit_test_files as files
 
+
 # Define a testing class. Absolutely fine to have one or multiple per file.
 # Each class must inherit unittest.TestCase
 class test_wod_read_data(unittest.TestCase):
     def test_load_wod(self):
-
         with self.subTest("Load profile data from WOD"):
-
             wod_profile_1D = coast.Profile(config=files.fn_wod_config)
             wod_profile_1D.read_wod(files.fn_wod)
 
diff --git a/unit_testing/test_xesmf_convert.py b/unit_testing/test_xesmf_convert.py
index 7ec9dd5a..a5b8eab7 100644
--- a/unit_testing/test_xesmf_convert.py
+++ b/unit_testing/test_xesmf_convert.py
@@ -4,15 +4,14 @@
 import os.path as path
 import unit_test_files as files
 
+
 # Single unit test. Can contain multiple test methods and subTests.
 class test_xesmf_convert(unittest.TestCase):
-
     # Test for conversion from gridded to xesmf.
     # Here I've used one test and then subtests for each smaller test.
     # This could also be split into multiple methods but the file would need
     # to be loaded multiple times. Using subtests allows a sequential testing.
     def test_basic_conversion_to_xesmf(self):
-
         # Read data files
         sci = coast.Gridded(files.fn_nemo_dat, files.fn_nemo_dom, config=files.fn_config_t_grid)
 
diff --git a/unit_testing/unit_test.py b/unit_testing/unit_test.py
index 023edc43..8dd11651 100644
--- a/unit_testing/unit_test.py
+++ b/unit_testing/unit_test.py
@@ -109,7 +109,6 @@
 
     # Open output file
     with open(fn_contents, "w") as file:
-
         # Write title things
         file.write("     UNIT TEST CONTENTS FILE TEST    \n")
         file.write("\n")
diff --git a/unit_testing/unit_test_contents.txt b/unit_testing/unit_test_contents.txt
index 5d250880..4038fa3b 100755
--- a/unit_testing/unit_test_contents.txt
+++ b/unit_testing/unit_test_contents.txt
@@ -11,6 +11,7 @@
       h. gridded_load_of_data_only
       i. gridded_load_subregion_with_domain
       j. gridded_load_variables_correctly_renamed
+      k. gridded_subset
 
 2. test_general_utils
       a. bst_to_gmt
@@ -27,9 +28,10 @@
       a. basic_conversion_to_xesmf
 
 5. test_diagnostic_methods
-      a. compute_vertical_spatial_derivative
-      b. construct_density
-      c. construct_pycnocline_depth_and_thickness
+      a. circulation
+      b. compute_vertical_spatial_derivative
+      c. construct_density
+      d. construct_pycnocline_depth_and_thickness
 
 6. test_transect_methods
       a. calculate_transport_velocity_and_depth