Output location works well for prepro.

Author: Eli

pyproject.toml

@@ -144,6 +144,7 @@ set_param = "schimpy.param:set_param_cli"
 bctide = "schimpy.bctide:bctide_cli"
 relocate_source_sink = "schimpy.relocate_source_sink:relocate_source_sink_cli"
 splice_netcdf = "schimpy.splice_netcdf:splice_netcdf_cli"
+check_vgrid = "schimpy.check_vgrid:check_vgrid_cli"
 
 # check_skewness = "schimpy.check_skewness:check_skewness_cli"
 # combine_flux = "schimpy.combine_flux:combine_flux_cli"

schimpy/batch_metrics.py

@@ -567,10 +567,8 @@ def plot(self):
                             s = simout.loc[:, (station_id, subloc)]
                             s = s.squeeze()
                             if isinstance(s, pd.DataFrame):
-                                raise (
-                                    "Station,Sublocation pair ({},{}) not unique or some other uniqueness problem in station files".format(
-                                        (station_id, subloc)
-                                    )
+                                raise Exception(
+                                    f"Station,Sublocation pair ({station_id},{subloc}) not unique or some other uniqueness problem in station files"
                                 )
                         else:
                             s = None

schimpy/create_vgrid_lsc2.py

@@ -194,6 +194,7 @@ def vgrid_gen(
     minmaxlayerfile,
     archive_nlayer="out",
     nlayer_gr3="nlayer.gr3",
+    diagnostics_dir=None,
 ):
 
     if not vgrid_version in ["5.8", "5.10"]:
@@ -203,6 +204,8 @@ def vgrid_gen(
 
     # specify nlayer output dir:
     out_dir = os.path.dirname(vgrid_out)
+    if diagnostics_dir is None:
+        diagnostics_dir = out_dir
     nlayer_gr3 = ensure_outdir(out_dir, nlayer_gr3)
 
     meshfun = BilinearMeshDensity()
@@ -268,9 +271,13 @@ def vgrid_gen(
 
         if archive_nlayer == "out":
             print("writing out number of layers")
+            nlayer_default_path = ensure_outdir(
+                diagnostics_dir,
+                os.path.basename(nlayer_gr3).replace(".gr3", "_default.gr3"),
+            )
             write_mesh(
                 mesh,
-                nlayer_gr3.replace(".gr3", "_default.gr3"),
+                nlayer_default_path,
                 node_attr=nlayer_default,
             )
             write_mesh(mesh, nlayer_gr3, node_attr=nlayer)
@@ -294,7 +301,7 @@ def vgrid_gen(
         maxlayer = nlayer * 0 + fixed_max  # np.max(maxlayer)
 
     sigma2, nlayer_revised = gen_sigma(
-        nlayer, minlayer, maxlayer, eta, h0, mesh, meshfun, out_dir=out_dir
+        nlayer, minlayer, maxlayer, eta, h0, mesh, meshfun, out_dir=diagnostics_dir
     )
     print("Returned nlayer revised: {}".format(np.max(nlayer_revised)))
     nlayer = nlayer_revised
@@ -309,7 +316,7 @@ def vgrid_gen(
     print("Done")
 
 
-def plot_vgrid(hgrid_file, vgrid0_file, vgrid_file, vgrid_version, eta, transectfiles):
+def plot_vgrid(hgrid_file, vgrid0_file, vgrid_file, vgrid_version, eta, transectfiles, out_dir="."):
     from schimpy.lsc2 import default_num_layers, plot_mesh
     from schimpy.schism_vertical_mesh import read_vmesh
     import matplotlib.pylab as plt
@@ -322,6 +329,9 @@ def plot_vgrid(hgrid_file, vgrid0_file, vgrid_file, vgrid_version, eta, transect
     h0 = mesh.nodes[:, 2]
     depth = eta + h0
 
+    images_dir = ospath.join(out_dir, "images")
+    os.makedirs(images_dir, exist_ok=True)
+
     zcor0 = vmesh0.build_z(mesh, eta)[:, ::-1]
     zcor1 = vmesh1.build_z(mesh, eta)[:, ::-1]
     for transectfile in transectfiles:
@@ -348,7 +358,7 @@ def plot_vgrid(hgrid_file, vgrid0_file, vgrid_file, vgrid_version, eta, transect
             plot_mesh(ax1, xpath, zcor1[path, :], 0, len(xpath), c="blue")
             ax0.plot(xpath, -h0[path], linewidth=2, c="black")
             ax1.plot(xpath, -h0[path], linewidth=2, c="black")
-            plt.savefig(ospath.join("images", base + ".png"))
+            plt.savefig(ospath.join(images_dir, base + ".png"))
             plt.show()
         except:
             print("Plotting of grid failed for transectfile: {}".format(transectfile))

schimpy/lsc2.py

@@ -522,9 +522,11 @@ def gen_sigma(
                 latest_layer = np.maximum(old_layer + 0.5 * speed, new_layer)
                 zsmooth[:, nsmoothlev - 2 - i] = latest_layer
                 old_layer = latest_layer
-            np.savetxt("zsmoothsave.txt", zsmooth)
+            zsmooth_fn = ensure_outdir(out_dir, "zsmoothsave.txt")
+            np.savetxt(zsmooth_fn, zsmooth)
         else:
-            zsmooth = np.loadtxt("zsmoothsave.txt")
+            zsmooth_fn = ensure_outdir(out_dir, "zsmoothsave.txt")
+            zsmooth = np.loadtxt(zsmooth_fn)
 
         print("Linking")
         # Now generate new mesh depths using the smoothed elevations as a pseudo-bed

schimpy/lsc2_v2.py

@@ -521,19 +521,22 @@ def fix_sigma_pileups(
         if abs(s[-1] - 1.0) <= TOL:
             s[-1] = 1.0
 
-        # 2) Bottom sliver fix (operate only if wet and we have an interior)
+        # 2a) Surface sliver fix: ensure the first cell (near σ=0, physical
+        #     surface) is at least tmin thick.
         bottom_adjusted = False
         bottom_drops = 0
+        bed_adjusted = False
+        bed_drops = 0
         dz1_after = np.nan
         dz2_after = np.nan
 
         if D > 0.0 and Ni >= 3 and tmin_i > 0.0:
             # Keep trying until first success or run out of interfaces
             while True:
-                # current bottom thickness
+                # current surface-cell thickness
                 dz = D * (s[1] - s[0])
                 if dz + 1e-15 >= tmin_i:
-                    bottom_adjusted = changed  # adjusted if we moved/dropped
+                    bottom_adjusted = changed
                     break
                 # try to raise s[1] to target
                 target = s[0] + (tmin_i / max(D, 1e-12))
@@ -548,13 +551,35 @@ def fix_sigma_pileups(
                 bottom_drops += 1
                 changed = True
                 if Ni < 3:
-                    # no more interior interface to adjust; infeasible
+                    break
+
+        # 2b) Bed sliver fix: ensure the last cell (near σ=1, physical
+        #     bottom / bed) is at least tmin thick.
+        if D > 0.0 and Ni >= 3 and tmin_i > 0.0:
+            while True:
+                dz = D * (s[-1] - s[-2])
+                if dz + 1e-15 >= tmin_i:
+                    bed_adjusted = changed
+                    break
+                # try to lower s[-2] to make room
+                target = s[-1] - (tmin_i / max(D, 1e-12))
+                if target > s[-3] + TOL:
+                    s[-2] = target
+                    changed = True
+                    bed_adjusted = True
+                    break
+                # blocked by s[-3] → drop s[-2] and retry
+                s = np.delete(s, -2)
+                Ni -= 1
+                bed_drops += 1
+                changed = True
+                if Ni < 3:
                     break
 
         # recompute dz1/dz2 after changes
         if D > 0.0 and Ni >= 3:
             dz1_after = D * (s[1] - s[0])
-            dz2_after = D * (s[2] - s[1])
+            dz2_after = D * (s[-1] - s[-2])
         elif D > 0.0 and Ni == 2:
             dz1_after = D  # single layer (bed→surface)
 
@@ -569,7 +594,8 @@ def fix_sigma_pileups(
                 "Ni_before": Ni0,
                 "Ni_after": Ni,
                 "dupes_removed": dupes_removed,
-                "bottom_drops": bottom_drops,
+                "surface_drops": bottom_drops,
+                "bed_drops": bed_drops,
                 "dz1_after_m": dz1_after,
                 "dz2_after_m": dz2_after,
                 "depth_m": D,
@@ -711,6 +737,8 @@ def smooth_interfaces_onepass_with_bed_fill_sparse(mesh,
 
         # Jacobi blend; update only nodes that actually have this level and are not bottom
         h_col_new = (1.0 - alpha) * hold[:, k] + alpha * nbr_mean
+        # Clamp: h must not exceed local depth (physically below bed)
+        h_col_new = np.minimum(h_col_new, depth)
         h[upd, k] = h_col_new[upd]
 
 def _get_row_stochastic_W(mesh, n_nodes):
@@ -1034,13 +1062,17 @@ def fit_interfaces(mesh,
                    valid: np.ndarray,
                    tbottom: np.ndarray,
                    params: FitParams = FitParams(),
-                   tmin_arr: Optional[np.ndarray] = None) -> Tuple[np.ndarray, np.ndarray]:
+                   tmin_arr: Optional[np.ndarray] = None,
+                   uniform_sigma_mask: Optional[np.ndarray] = None) -> Tuple[np.ndarray, np.ndarray]:
     """Fit interfaces h_k(x) to targets with smoothing, bottom anchoring, and guards.
 
     Parameters
     ----------
     tmin_arr : optional float array (n_nodes,)
         Per-node minimum layer thickness.  If None, uses params.tmin everywhere.
+    uniform_sigma_mask : optional bool array (n_nodes,)
+        If True for a node, reassert uniform sigma spacing in finalization
+        (protects shallow exception nodes from smoother cliff-pressure).
 
     Returns
     -------
@@ -1368,7 +1400,7 @@ def _dbg_levels(tag):
         _assert_finite("h post FINISH (exist)", h, where_mask=exists)
 
 
-    # Finalization: single forward guard + exact endpoints (+ strict nudge)
+    # Finalization: exact endpoints + guards (+ strict nudge)
     for i in range(n):
         Ni = int(Nlevels[i]); D = depth[i]
         if Ni < 2:
@@ -1381,6 +1413,14 @@ def _dbg_levels(tag):
                 h[i, k] = D_eff * k / (Ni - 1)
             continue
 
+        # Exception nodes (shallow, forced Nlevels by hysteresis): reassert
+        # uniform sigma. The smoother pulls these toward deep neighbors;
+        # uniform is the correct answer for them.
+        if uniform_sigma_mask is not None and uniform_sigma_mask[i]:
+            for k in range(Ni):
+                h[i, k] = k * D / (Ni - 1)
+            continue
+
         # Surface exactly 0
         h[i, 0] = 0.0
 
@@ -1393,6 +1433,25 @@ def _dbg_levels(tag):
         # Pin bed exactly to D
         h[i, Ni - 1] = D
 
+        # Overflow fix: if smoother pushed levels above D, or the bed
+        # cell would be thinner than tmin, redistribute proportionally
+        # between a valid base and the bed.  Choose the split point so
+        # that the redistributed spacing >= tmin.
+        if Ni >= 3 and h[i, Ni - 2] > D - tmin_arr[i]:
+            tmin_i = tmin_arr[i]
+            # Scan downward: find the deepest level k whose value allows
+            # (Ni-1-k) layers of at least tmin between h[k] and D.
+            k_overflow = 1  # fallback: redistribute everything from level 1
+            for k in range(Ni - 2, 0, -1):
+                if h[i, k] <= D - (Ni - 1 - k) * tmin_i:
+                    k_overflow = k + 1
+                    break
+            # Linearly redistribute from h[k_overflow-1] to D
+            h_base = h[i, k_overflow - 1]
+            n_redist = Ni - k_overflow  # number of levels to redistribute
+            for j in range(n_redist):
+                h[i, k_overflow + j] = h_base + (j + 1) * (D - h_base) / n_redist
+
         # Strict separation nudge (fp safety)
         for k in range(1, Ni):
             if not (h[i, k] > h[i, k - 1]):
@@ -1694,7 +1753,8 @@ def run_pipeline(mesh,
     t3 = time.perf_counter()
     logger.info("D) Fitting interfaces ...")
     h, valid2 = fit_interfaces(mesh, depth, Nlevels, hhat, valid, tb, pp.fit,
-                               tmin_arr=tmin_arr)
+                               tmin_arr=tmin_arr,
+                               uniform_sigma_mask=uniform_sigma_mask)
     logger.info("   done in %.2f s", time.perf_counter() - t3)
 
     # E) sigma

schimpy/nudging.py

@@ -228,7 +228,7 @@ def concatenate_nudge(
             # check to see if there is any overlap among the regions, if yes,
             # take the weighted average
             weights_merged = np.zeros_like(weights_var[0][0])
-            values_merged = np.zeros((len(self.time), len(imap_merged), self.nvrt))
+            values_merged = np.zeros((len(self.time), len(imap_merged), self.nvrt), dtype=np.float32)
             for im in range(len(imap_merged)):
                 idx_r = np.where(~np.isnan(idx_list[im, :]))[0].astype(int)
                 if len(idx_r) > 1:  # overlapping regions
@@ -400,9 +400,11 @@ def gen_nudge_3dfield(self, region_info):
 
         utm_xy = np.array([self.node_x, self.node_y])
 
-        temperature = []
-        salinity = []
+        temperature = None  # pre-allocated after npout is known
+        salinity = None
         output_day = []
+        istep = 0
+        max_steps = len(self.datetime)
         var_name = {
             "temperature": "temp",
             "salinity": "salt",
@@ -468,7 +470,7 @@ def gen_nudge_3dfield(self, region_info):
                 ctime = pd.to_datetime(t)
                 dt = ctime - self.start_date
                 dt_in_days = dt.total_seconds() / 86400.0
-                print(f"Time out at hr stage (days)={ctime}, dt ={dt_in_days} ")
+                logger.debug("Time out at hr stage (days)=%s, dt =%s", ctime, dt_in_days)
                 irecout += 1
 
                 # Make sure it2=ilo is output (output at precisely the time interval)
@@ -487,14 +489,14 @@ def gen_nudge_3dfield(self, region_info):
                     .sel(time=t)
                     .transpose("lon", "lat", "depth")
                     .values[:, :, -1::-1]
-                )
+                ).astype(np.float32)
                 # if 'salt' in var_list:
                 salt = (
                     ncdata["salt"]
                     .sel(time=t)
                     .transpose("lon", "lat", "depth")
                     .values[:, :, -1::-1]
-                )
+                ).astype(np.float32)
                 # Comments below for futhre implementation.
                 # if 'uvel' in var_list:
                 #     uvel = ncdata['uvel'].sel(time=t).transpose(
@@ -694,19 +696,23 @@ def gen_nudge_3dfield(self, region_info):
                     vrat = vrat.T
 
                     # initialize interpolation coefficients
-                    wild2 = np.zeros((self.nvrt, npout, 4, 2))
+                    wild2 = np.zeros((self.nvrt, npout, 4, 2), dtype=np.float32)
                     ix = np.broadcast_to(ix, (self.nvrt, npout))
                     iy = np.broadcast_to(iy, (self.nvrt, npout))
                     intri1 = np.where(intri == 1)[0]
                     intri2 = np.where(intri == 2)[0]
                     i_in_1 = i_in[intri1]
                     i_in_2 = i_in[intri2]
 
-                    tempout = np.empty((self.nvrt, self.nnode))
-                    saltout = np.empty((self.nvrt, self.nnode))
+                    tempout = np.empty((self.nvrt, self.nnode), dtype=np.float32)
+                    saltout = np.empty((self.nvrt, self.nnode), dtype=np.float32)
                     tempout[:, i_out] = tem_outside
                     saltout[:, i_out] = sal_outside
 
+                    # Pre-allocate output arrays as float32
+                    temperature = np.empty((max_steps, self.nvrt, npout), dtype=np.float32)
+                    salinity = np.empty((max_steps, self.nvrt, npout), dtype=np.float32)
+
                 id_dry = np.where(kbp[ix, iy] == -1)[0]  # 3dfield dry cell
                 lev[:, id_dry] = int(ilen - 1 - 1)
                 vrat[:, id_dry] = 1
@@ -771,15 +777,18 @@ def gen_nudge_3dfield(self, region_info):
                 tempout[idST[0], i_in[idST[1]]] = tempout[idST[0], i_in[idST[1]]] - 1
 
                 logger.debug("applying spatial interpolation!")
-                logger.info("outputting at day , %s, %s", dt.total_seconds()/86400, npout)
+                current_day = int(dt.total_seconds() / 86400)
+                if not hasattr(self, '_last_logged_day') or current_day != self._last_logged_day:
+                    logger.info("outputting at day %s, npout=%s", current_day, npout)
+                    self._last_logged_day = current_day
+                else:
+                    logger.debug("outputting at day %s, npout=%s", dt.total_seconds()/86400, npout)
 
                 # only save the nodes with valid values.
-                tempout_in = [temp_t[imap[:npout]] for temp_t in tempout]
-                saltout_in = [salt_t[imap[:npout]] for salt_t in saltout]
-
-                temperature.append(tempout_in)
-                salinity.append(saltout_in)
+                temperature[istep, :, :] = tempout[:, imap[:npout]]
+                salinity[istep, :, :] = saltout[:, imap[:npout]]
                 output_day.append(dt.total_seconds() / 86400)
+                istep += 1
                 nudge_step = int(pd.Timedelta(self.nudge_step).total_seconds() / 3600)
                 if len(output_day) > 1:
                     if output_day[-1] - output_day[-2] > (nudge_step + 0.1) / 24:
@@ -790,13 +799,13 @@ def gen_nudge_3dfield(self, region_info):
                             "The difference between current and previous time step is greater than assigned stride"
                         )
 
-        temperature = np.array(temperature)
-        salinity = np.array(salinity)
+        temperature = temperature[:istep]
+        salinity = salinity[:istep]
         # [var,time,node_map, nvrt]
         temperature = np.transpose(temperature, (0, 2, 1))
         salinity = np.transpose(salinity, (0, 2, 1))
         # Enforce lower bound for temp. for eqstate
-        temperature[temperature < 0] == 0
+        temperature[temperature < 0] = 0
         logger.info("reorganizing matrix!")
         # print("time=%f"%(timer.time() - start))
         max_time = self.time[-1].total_seconds() / 3600 / 24