#216 Bug where velocity sometimes has NaNs in bottom level

inputs/grid_name_map.json

@@ -40,7 +40,10 @@
             "e3f": "e3f",
             "e3uw": "e3uw",
             "e3vw": "e3vw",
-            "e3fw": "e3fw"
+            "e3fw": "e3fw",
+            "ln_zco": "ln_zco",
+            "ln_zps": "ln_zps",
+            "ln_sco": "ln_sco"
     },
     "dst_variable_map": {
             "nav_lon": "nav_lon",
@@ -77,6 +80,9 @@
             "e3f": "e3f",
             "e3uw": "e3uw",
             "e3vw": "e3vw",
-            "e3fw": "e3fw"
+            "e3fw": "e3fw",
+            "ln_zco": "ln_zco",
+            "ln_zps": "ln_zps",
+            "ln_sco": "ln_sco"
     }
 }

inputs/grid_name_map_readme.txt

@@ -75,3 +75,7 @@ In all cases "t" should be size 1. Pybdy does not deal with time varying grids.
 "e3uw" = * vertical scale factor distance between w-levels on u-grid (dims [t, z, y, x])
 "e3vw" = * vertical scale factor distance between w-levels on v-grid (dims [t, z, y, x])
 "e3fw" = * vertical scale factor distance between w-levels on f-grid (dims [t, z, y, x])
+
+"ln_zco" = * flag for zco vertical coordinates 1 for true, 0 for false often provided in domain_cfg.nc
+"ln_zps" = * flag for zps vertical coordinates 1 for true, 0 for false often provided in domain_cfg.nc
+"ln_sco" = * flag for sco vertical coordinates 1 for true, 0 for false often provided in domain_cfg.nc

inputs/namelist_local.bdy

@@ -9,6 +9,15 @@
 !!
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
+!------------------------------------------------------------------------------
+!   vertical coordinate
+!------------------------------------------------------------------------------
+   sn_src_zgr_type = 'zps' ! vertical coordinate type: 'zco', 'zps' or 'sco'
+   sn_dst_zgr_type = 'zps' ! vertical coordinate type: 'zco', 'zps' or 'sco'
+                           ! 'zco' is z-coordinate - full    steps
+                           ! 'zps' is z-coordinate - partial steps
+                           ! 'sco' is s- or hybrid z-s-coordinate
+
 !------------------------------------------------------------------------------
 !  grid information
 !------------------------------------------------------------------------------

inputs/namelist_remote.bdy

@@ -9,6 +9,15 @@
 !!
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
+!------------------------------------------------------------------------------
+!   vertical coordinate
+!------------------------------------------------------------------------------
+   sn_src_zgr_type = 'sco' ! vertical coordinate type: 'zco', 'zps' or 'sco'
+   sn_dst_zgr_type = 'sco' ! vertical coordinate type: 'zco', 'zps' or 'sco'
+                           ! 'zco' is z-coordinate - full    steps
+                           ! 'zps' is z-coordinate - partial steps
+                           ! 'sco' is s- or hybrid z-s-coordinate
+
 !------------------------------------------------------------------------------
 !  grid information
 !------------------------------------------------------------------------------

src/grid/zgr.py

@@ -32,13 +32,16 @@
 
 
 class Z_Grid:
-    def __init__(self, zgr_file, name_map_file, hgr_type, e_dict, logger, dst=1):
+    def __init__(
+        self, zgr_file, zgr_type, name_map_file, hgr_type, e_dict, logger, dst=1
+    ):
         """
         Master depth class.
 
         Parameters
         ----------
             zgr_file (str)           : string of file for loading zgr data
+            zgr_type (str)           : zgr type from namelist zco, zps or sco
             name_map_file (str)      : string of file for mapping variable names
             hgr_type (str)           : horizontal grid type
             e_dict (dict)       : dictionary of e1 and e2 scale factors
@@ -84,15 +87,15 @@ def __init__(self, zgr_file, name_map_file, hgr_type, e_dict, logger, dst=1):
             "e3uw",
             "e3vw",
             "e3fw",
+            "ln_zco",
+            "ln_zps",
+            "ln_sco",
         ]
-        #    "gdepw_0",
-        #    "e3t_0",
-        #    "e3w_0",
 
         self.get_vars(vars_want)
 
         # Work out what sort of source grid we have
-        self.find_zgrid_type()
+        self.find_zgrid_type(zgr_type)
 
         # Fill in missing variables we need for the grid type
         missing_vars = sorted(list(set(vars_want) - set(self.var_list)))
@@ -138,59 +141,36 @@ def get_vars(self, vars_want):
         nc.close()
         self.var_list = list(self.grid.keys())
 
-    def find_zgrid_type(self):
+    def find_zgrid_type(self, zgr_type):
         """Find out what type of vertical grid is provided zco, zps or sigma levels (sco)."""
-        if ("gdept" not in self.var_list) & ("gdept_0" not in self.var_list):
-            if "e3t" in self.var_list:
-                self.grid["gdept"] = np.cumsum(self.grid["e3t"], axis=1)
-                self.var_list.append("gdept")
-            else:
-                raise Exception("No gdept or gdept_0 variable present in zgr file.")
-
-        if "gdept" not in self.var_list:
-            self.grid_type = "zco"
+        if zgr_type in ["zco", "zps", "sco"]:
+            # Check namelist for user specified zgr type
+            self.grid_type = zgr_type
         else:
-            if "mbathy" not in self.var_list:
-                raise Exception("No mbathy variable present in zgr file.")
-            # Could still be z, z-partial-step (zps) or sigma
-            # "z" - Check if all levels are equal
-            x_diff = np.diff(self.grid["gdept"], axis=3).sum() == 0
-            y_diff = np.diff(self.grid["gdept"], axis=2).sum() == 0
-            dep_test_z = x_diff & y_diff
-
-            # "zps" - Select levels 2 above bathy to remove partial steps.
-            # Then check if levels of deepest profile are equal to all levels
-            z_ind = np.indices(self.grid["gdept"].shape)[1]
-            m_tile = np.tile(self.grid["mbathy"], (self.grid["gdept"].shape[1], 1, 1))[
-                np.newaxis, ...
-            ]
-            mask_bottom = z_ind >= (m_tile - 2)
-            dep_mask = np.ma.masked_where(mask_bottom, self.grid["gdept"])
-
-            sel = np.nonzero(np.ma.max(dep_mask) == dep_mask)
-            dep_deepest = np.tile(
-                dep_mask[sel[0][0], :, sel[2][0], sel[3][0]],
-                (self.grid["gdept"].shape[3], self.grid["gdept"].shape[2], 1),
-            ).T[np.newaxis, :, :, :]
-            dep_test_zps = (
-                dep_deepest[mask_bottom is False]
-                == self.grid["gdept"][mask_bottom is False]
-            )
-
-            # Check if all levels are inside the bathy anywhere
-            lev_deepest = np.ma.max(self.grid["mbathy"])
-            dep_test_sigma = lev_deepest >= (self.grid["gdept"].shape[1] - 1)
-            if dep_test_z:
-                # z-level
-                self.grid_type = "zco"
-            elif dep_test_zps.all():
-                # z-partial-step
-                self.grid_type = "zps"
-            elif dep_test_sigma:
-                # sigma-level
-                self.grid_type = "sco"
-            else:
-                raise Exception("Unknown/unaccounted for vertical grid type.")
+            raise Exception("zgr_type not specified in namelist or file.")
+
+        file_zgr = ""
+        if "ln_zco" in self.var_list:
+            # Check the zgr netcdf file for ln flags
+            if self.grid["ln_zco"]:
+                file_zgr = "zco"
+        elif "ln_zps" in self.var_list:
+            if self.grid["ln_zps"]:
+                file_zgr = "zps"
+        elif "ln_sco" in self.var_list:
+            if self.grid["ln_sco"]:
+                file_zgr = "sco"
+
+        if file_zgr != "":
+            if file_zgr != self.grid_type:
+                print(
+                    "Warning: the vertical grid in the zgr file ("
+                    + file_zgr
+                    + ") does not match the user defined"
+                    + "zgr type in the namelist ("
+                    + self.grid_type
+                    + ")"
+                )
 
         self.logger.info("Vertical grid is type: " + self.grid_type)
 

src/pybdy/nemo_bdy_extr_assist.py

@@ -332,6 +332,7 @@ def flood_fill(sc_bdy, isslab, logger):
     sc_shape = sc_bdy.shape
 
     for i in range(sc_shape[0]):
+        # Check if layer centre-point has both nans and values
         while np.isnan(sc_bdy[i, :, 0]).any() & (~np.isnan(sc_bdy[i, :, 0])).any():
             # Flood sc land horizontally within the chunk for the centre point.
             # This may not be perfect but better than filling with zeros
@@ -351,6 +352,12 @@ def flood_fill(sc_bdy, isslab, logger):
             np.isnan(sc_bdy[:, :, 0])
         ]
 
+    # Replace nans around centre-point with value.
+    sc_nan = np.isnan(sc_bdy)
+    sc_bdy[:, :, 1:][sc_nan[:, :, 1:]] = np.tile(sc_bdy[:, :, 0:1], (1, 1, 8))[
+        sc_nan[:, :, 1:]
+    ]
+
     return sc_bdy
 
 
@@ -363,7 +370,7 @@ def interp_vertical(sc_bdy, dst_dep, bdy_bathy, z_ind, z_dist, num_bdy, zinterp=
     sc_bdy (np.array)   : souce data [nz_sc, nbdy, 9]
     dst_dep (np.array)  : the depth of the destination grid chunk [nz, nbdy]
     bdy_bathy (np.array): the destination grid bdy points bathymetry
-    z_ind (np.array)    : the indices of the sc depth above and below bdy
+    z_ind (np.array)    : the indices of the sc depth above and below bdy point
     z_dist (np.array)   : the distance weights of the selected points
     num_bdy (int)       : number of boundary points in chunk
     zinterp (bool)      : vertical interpolation flag

src/pybdy/nemo_bdy_extr_tm3.py

@@ -61,8 +61,8 @@ def __init__(self, setup, SourceCoord, DstCoord, Grid, var_nam, grd, pair):
         DstCoord        (obj)  : destination grid information
         Grid            (dict) : containing grid type 't', 'u', 'v' and source time
         var_name        (list) : netcdf file variable names (str)
-        years           (list) : years to extract (default [1979])
-        months          (list) : months to extract (default [11])
+        grd             (str)  : grid to process 't', 'u', 'v'
+        pair            (str)  : None or 'uv'
 
         Returns
         -------
@@ -982,6 +982,7 @@ def extract_month(self, year, month):
                         self.logger.info(
                             "%s en to %s %s", self.rot_dir, self.rot_dir, dst_bdy.shape
                         )
+
                         dst_bdy = rot_rep(
                             dst_bdy,
                             dst_bdy_2,

src/pybdy/nemo_bdy_zgrv2.py

@@ -169,9 +169,9 @@ def get_bdy_depths(DstCoord, bdy_i, grd):
     bdy_tz = np.ma.zeros((m_t.shape[0], len(g_ind)))
     bdy_e3 = np.ma.zeros((m_e.shape[0], len(g_ind)))
     for k in range(m_w.shape[0]):
-        tmp_w = np.ma.masked_where(mbathy + 1 < k + 1, m_w[k, :, :])
-        tmp_t = np.ma.masked_where(mbathy < k + 1, m_t[k, :, :])
-        tmp_e = np.ma.masked_where(mbathy < k + 1, m_e[k, :, :])
+        tmp_w = np.ma.masked_where(mbathy + 1 < (k + 1), m_w[k, :, :])
+        tmp_t = np.ma.masked_where(mbathy < (k + 1), m_t[k, :, :])
+        tmp_e = np.ma.masked_where(mbathy < (k + 1), m_e[k, :, :])
 
         tmp_w = tmp_w.flatten("F")
         tmp_t = tmp_t.flatten("F")

src/pybdy/profiler.py

@@ -107,6 +107,7 @@ def process_bdy(setup_filepath=0, mask_gui=False):
     DstCoord.hgr = hgr.H_Grid(settings["dst_hgr"], settings["nme_map"], logger, dst=1)
     DstCoord.zgr = zgr.Z_Grid(
         settings["dst_zgr"],
+        settings["dst_zgr_type"],
         settings["nme_map"],
         DstCoord.hgr.grid_type,
         DstCoord.hgr.grid,
@@ -175,6 +176,7 @@ def process_bdy(setup_filepath=0, mask_gui=False):
     )
     SourceCoord.zgr = zgr.Z_Grid(
         settings["src_zgr"],
+        settings["src_zgr_type"],
         settings["nme_map"],
         SourceCoord.hgr.grid_type,
         SourceCoord.hgr.grid,

tests/data/grid_name_map.json

@@ -40,7 +40,10 @@
             "e3f": "e3f",
             "e3uw": "e3uw",
             "e3vw": "e3vw",
-            "e3fw": "e3fw"
+            "e3fw": "e3fw",
+            "ln_zco": "ln_zco",
+            "ln_zps": "ln_zps",
+            "ln_sco": "ln_sco"
     },
     "dst_variable_map": {
             "nav_lon": "nav_lon",
@@ -77,6 +80,9 @@
             "e3f": "e3f",
             "e3uw": "e3uw",
             "e3vw": "e3vw",
-            "e3fw": "e3fw"
+            "e3fw": "e3fw",
+            "ln_zco": "ln_zco",
+            "ln_zps": "ln_zps",
+            "ln_sco": "ln_sco"
     }
 }

tests/data/namelist_zz_end_to_end.bdy

@@ -9,6 +9,15 @@
 !!
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
+!------------------------------------------------------------------------------
+!   vertical coordinate
+!------------------------------------------------------------------------------
+   sn_src_zgr_type = 'zco' ! vertical coordinate type: 'zco', 'zps' or 'sco'
+   sn_dst_zgr_type = 'zco' ! vertical coordinate type: 'zco', 'zps' or 'sco'
+                           ! 'zco' is z-coordinate - full    steps
+                           ! 'zps' is z-coordinate - partial steps
+                           ! 'sco' is s- or hybrid z-s-coordinate
+
 !------------------------------------------------------------------------------
 !  grid information
 !------------------------------------------------------------------------------

tests/test_zgr.py

@@ -65,7 +65,7 @@ def test_depth_file_zps():
         e_dict = {k: hg.grid[k] for k in keys}
 
         # calc vertical grid
-        zg = zgr.Z_Grid(in_file, name_map, hg.grid_type, e_dict, logger)
+        zg = zgr.Z_Grid(in_file, "zps", name_map, hg.grid_type, e_dict, logger)
 
         nc = GetFile(bench_file)
         e3t = nc.nc["e3t"][:]