fix(pacific): keep ATL, SECOFS and AK behaviour unchanged by the Pacific features (into #74)

nos_utils/config.py

@@ -238,14 +238,16 @@ class ForcingConfig:
     # ADT satellite SSH blending (STOFS-3D-ATL uses CMEMS ADT to correct boundary SSH)
     adt_enabled: bool = False
     adt_datum_offset: float = -0.45
+    # Rerun archive/tar also carries temnu/salnu/surface_restore with the Pacific names; off keeps main's set. MJ (10/24/26)
+    obc_archive_extended: bool = False
     # Ops ESMF map for the ADT regrid (yaml ocean.adt.weight_file); a bare name resolves against
     # FIXofs. Unset or missing falls back to the bilinear regrid. MJ (10/02/26)
     adt_weight_file: Optional[Path] = None
     # Nudging enabled and timescale
     nudging_enabled: bool = False
     nudging_timescale_seconds: float = 86400.0
     # Temporal resolution of the nudging files (seconds).
-    nudging_dt: float = 21600.0
+    nudging_dt: float = 10800.0
     # Apply 2-step temporal moving average to nudging fields.
     nudging_time_smoothing: bool = False
     # Surface restore enabled and smoothing
@@ -436,6 +438,8 @@ class ForcingConfig:
     # OFSes have shorter OBC windows (typically ~10 records for 54h runs)
     # and would be wrongly flagged if this defaulted to 21.
     obc_min_timesteps: int = 0
+    # Minimum RTOFS slots in the ops OBC timeline (ops: 23 6-hourly records, n012 to f120). MJ (10/24/26)
+    obc_min_slots: int = 23
 
     # --- Tidal settings ---
     # Tidal constituents to use
@@ -707,7 +711,7 @@ def for_stofs_3d_pac(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig":
             lon_min=93.1098, lon_max=289.6750,
             lat_min=-33.0388, lat_max=66.3254,
             pdy=pdy, cyc=cyc,
-            nowcast_hours=24, forecast_hours=108,
+            nowcast_hours=24, forecast_hours=48,
             gfs_resolution="0p25",
             met_num=2, n_levels=84,
             model_dt=90.0,
@@ -724,6 +728,8 @@ def for_stofs_3d_pac(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig":
             nudging_timescale_seconds=86400.0,
             nwm_product="medium_range_mem1",
             obc_min_timesteps=4,
+            obc_min_slots=4,
+            obc_archive_extended=True,
             obc_dt=86400.0,
             nudging_dt=86400.0,
             obc_time_smoothing=True,
@@ -732,8 +738,8 @@ def for_stofs_3d_pac(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig":
             obc_roi_2d={"x1": 222, "x2": 2703, "y1": 1051, "y2": 2641, "x_stride": 2, "y_stride": 2},
             obc_roi_3d={"x1": 222, "x2": 2703, "y1": 1051, "y2": 2641, "x_stride": 2, "y_stride": 2},
             nudge_roi_3d={"x1": 222, "x2": 2703, "y1": 1051, "y2": 2641, "x_stride": 2, "y_stride": 2},
-            surface_restore_enabled=True,
-            surface_restore_time_smoothing=True,
+            # Surface restore is untested (reads T and S from one file, restore nodes unverified): off. MJ (10/24/26)
+            surface_restore_enabled=False,
         )
         defaults.update(overrides)
         return cls(**defaults)
@@ -750,7 +756,7 @@ def for_stofs_3d_pac_ufs(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig
             lon_min=93.1098, lon_max=289.6750,
             lat_min=-33.0388, lat_max=66.3254,
             pdy=pdy, cyc=cyc,
-            nowcast_hours=24, forecast_hours=108,
+            nowcast_hours=24, forecast_hours=48,
             gfs_resolution="0p25",
             met_num=2, nws=4, n_levels=84,
             model_dt=90.0,
@@ -764,6 +770,8 @@ def for_stofs_3d_pac_ufs(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig
             nudging_timescale_seconds=86400.0,
             nwm_product="medium_range_mem1",
             obc_min_timesteps=4,
+            obc_min_slots=4,
+            obc_archive_extended=True,
             obc_dt=86400.0,
             nudging_dt=86400.0,
             obc_time_smoothing=True,
@@ -772,8 +780,8 @@ def for_stofs_3d_pac_ufs(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig
             obc_roi_2d={"x1": 222, "x2": 2703, "y1": 1051, "y2": 2641, "x_stride": 2, "y_stride": 2},
             obc_roi_3d={"x1": 222, "x2": 2703, "y1": 1051, "y2": 2641, "x_stride": 2, "y_stride": 2},
             nudge_roi_3d={"x1": 222, "x2": 2703, "y1": 1051, "y2": 2641, "x_stride": 2, "y_stride": 2},
-            surface_restore_enabled=True,
-            surface_restore_time_smoothing=True,
+            # Surface restore is untested (reads T and S from one file, restore nodes unverified): off. MJ (10/24/26)
+            surface_restore_enabled=False,
             # DATM grid — STOFS3D_PAC preset, verified from weight_gfs.nc
             # nx = round((290.0-92.5)/0.025)+1 = 7901; ny = round((67.0-(-33.5))/0.025)+1 = 4021
             datm_lon_min=92.5, datm_lon_max=290.0,
@@ -995,7 +1003,7 @@ def from_yaml(cls, yaml_path, **overrides) -> "ForcingConfig":
         if _adt_wt is None and str(_sys.get("name", "")).startswith("stofs_3d_atl"):
             _adt_wt = "stofs_3d_atl_ufs.adt_weight.nc"
         _adt_wt = Path(_adt_wt) if _adt_wt else None
-        _adt_off = float(adt.get("datum_offset", 0.0)) if isinstance(adt, dict) else 0.0
+        _adt_off = float(adt.get("datum_offset", -0.45)) if isinstance(adt, dict) else -0.45
 
         # Ops river inputs (static FIX source_sink.in / msource.th / vsink.th, ops source
         # order): STOFS-3D-ATL by name (the coupled executable now has PREC_EVAP), yaml
@@ -1074,13 +1082,17 @@ def from_yaml(cls, yaml_path, **overrides) -> "ForcingConfig":
         if rtofs_3d_region:
             kwargs["rtofs_3d_region"] = str(rtofs_3d_region)
 
+        ocean = ocean if isinstance(ocean, dict) else {}
+        obc = obc if isinstance(obc, dict) else {}
+        nudge = nudge if isinstance(nudge, dict) else {}
         kwargs["rtofs_daily_3z"] = _strict_bool(ocean.get("rtofs_daily_3z"), "ocean.rtofs_daily_3z", False)
         kwargs["obc_dt"] = float(obc.get("dt", 21600.0))
         kwargs["obc_time_smoothing"] = _strict_bool(obc.get("time_smoothing"), "obc.time_smoothing", False)
-        kwargs["nudging_dt"] = float(nudge.get("dt", kwargs["obc_dt"]))
+        kwargs["nudging_dt"] = float(nudge.get("dt", 10800.0))
         kwargs["nudging_time_smoothing"] = _strict_bool(nudge.get("time_smoothing"), "nudging.time_smoothing", False)
 
-        surf_res = ocean.get("surface_restore", {}) if isinstance(ocean, dict) else {}
+        surf_res = ocean.get("surface_restore")
+        surf_res = surf_res if isinstance(surf_res, dict) else {}
         kwargs["surface_restore_enabled"] = _strict_bool(surf_res.get("enabled"), "ocean.surface_restore.enabled", False)
         kwargs["surface_restore_time_smoothing"] = _strict_bool(surf_res.get("time_smoothing"), "ocean.surface_restore.time_smoothing", False)
         if surf_res.get("weight_file"):
@@ -1182,6 +1194,10 @@ def from_yaml(cls, yaml_path, **overrides) -> "ForcingConfig":
             kwargs["obc_min_timesteps"] = int(yaml_min_ts)
         elif is_stofs_like:
             kwargs["obc_min_timesteps"] = 21
+        if obc.get("min_slots") is not None:
+            kwargs["obc_min_slots"] = int(obc["min_slots"])
+        kwargs["obc_archive_extended"] = _strict_bool(
+            obc.get("archive_extended"), "obc.archive_extended", False)
 
         # Declarative prep extras override. When a root-level `prep.extras`
         # block is present it is authoritative and supersedes the

nos_utils/forcing/adt.py

@@ -50,6 +50,10 @@
 
 # Ops ADT subset ROI, used only when the SSH_1 destination grid cannot be read.
 # MJ (10/01/26)
+ADT_LON_MIN = -62.5
+ADT_LON_MAX = -51.5
+ADT_LAT_MIN = 7.0
+ADT_LAT_MAX = 54.0
 ADT_PAD = 1.0
 
 # Ops holds ADT-missing cells as the CMEMS int fill carried through float32. MJ (10/02/26)
@@ -145,7 +149,7 @@ def blend_ssh(self, ssh_path: Path, work_dir: Path) -> Optional[Path]:
 
     def _dest_bounds(self, ssh_path: Path):
         """ADT subset (lon_min, lon_max, lat_min, lat_max): SSH_1 grid extent padded
-        by ADT_PAD, else the config domain bounds. MJ (10/01/26)
+        by ADT_PAD, else the ops ROI (-180/180 systems) or the padded config domain (0-360). MJ (10/01/26)
         """
         try:
             with Dataset(str(ssh_path)) as ds:
@@ -154,10 +158,22 @@ def _dest_bounds(self, ssh_path: Path):
             return (float(np.nanmin(lon)) - ADT_PAD, float(np.nanmax(lon)) + ADT_PAD,
                     float(np.nanmin(lat)) - ADT_PAD, float(np.nanmax(lat)) + ADT_PAD)
         except Exception as e:
-            log.warning(f"Cannot read SSH_1 grid extent ({e}) — using config domain bounds")
+            log.warning(f"Cannot read SSH_1 grid extent ({e}) — using the fallback ROI")
+            if self._pm180():
+                return ADT_LON_MIN, ADT_LON_MAX, ADT_LAT_MIN, ADT_LAT_MAX
             lon_min, lon_max, lat_min, lat_max = self.config.domain
             return lon_min - ADT_PAD, lon_max + ADT_PAD, lat_min - ADT_PAD, lat_max + ADT_PAD
 
+    def _pm180(self) -> bool:
+        from ..coords import lon_convention
+        return lon_convention(self.config) != "0360"
+
+    def _roi(self, ssh_path):
+        """Ops fixed ROI for -180/180 systems (STOFS-3D-ATL); SSH_1 extent for the 0-360 Pacific. MJ (10/24/26)"""
+        if self._pm180():
+            return ADT_LON_MIN, ADT_LON_MAX, ADT_LAT_MIN, ADT_LAT_MAX
+        return self._dest_bounds(ssh_path)
+
     def _find_adt_data(self) -> Optional[Path]:
         """Newest available CMEMS ADT file (valid day, else previous day)."""
         files = self._find_adt_files()
@@ -309,7 +325,7 @@ def dlon_diff(a, b):
             if row.min() < 0 or row.max() >= n_b or col.min() < 0 or col.max() >= n_a:
                 raise ValueError("map row/col out of range")
 
-            bounds = self._dest_bounds(ssh_path)
+            bounds = self._roi(ssh_path)
             days = []
             for f in adt_files:
                 a, lons, lats = self._subset_ops_roi(f, bounds)
@@ -323,7 +339,8 @@ def dlon_diff(a, b):
                         or np.abs(yca - LA.ravel()).max() > ADT_COORD_TOL):
                     raise ValueError("ADT ROI coordinates do not match map xc_a/yc_a")
                 out = self._apply_esmf_map(S, row, col, n_b, a.ravel())
-                days.append(out.astype(np.float32) if self.ops_numerics else out)
+                off = self.config.adt_datum_offset
+                days.append((out + off).astype(np.float32) if self.ops_numerics else out + off)
             if self.ops_numerics:
                 stack = np.stack(days)
                 ok = np.isfinite(stack)
@@ -429,7 +446,7 @@ def _read_adt(self, adt_path: Path, bounds=None):
         bounds is (lon_min, lon_max, lat_min, lat_max) in -180/180; default is the
         ops ROI. Returns (adt[lat, lon], lons, lats) or None. Fill values become NaN.
         """
-        lon_min, lon_max, lat_min, lat_max = bounds or self._dest_bounds(None)
+        lon_min, lon_max, lat_min, lat_max = bounds or self._roi(None)
         try:
             ds = Dataset(str(adt_path))
             try:
@@ -439,7 +456,7 @@ def _read_adt(self, adt_path: Path, bounds=None):
                 lons = np.array(ds.variables[lon_name][:], dtype=np.float64)
                 lats = np.array(ds.variables[lat_name][:], dtype=np.float64)
                 adt_var = "adt" if "adt" in ds.variables else "surf_el"
-                adt_data = ds.variables[adt_var][...]
+                adt_data = ds.variables[adt_var][...] + self.config.adt_datum_offset
 
                 from ..coords import lon_convention, shift_grid_0360
                 if lon_convention(self.config) == "0360":
@@ -568,7 +585,7 @@ def _apply_adt_blend(self, ssh_path: Path, adt_field, work_dir: Path,
         """Apply the ops ADT blend to SSH_1.nc.
 
         Formula: SSH_final(t) = SSH_rtofs(t) - SSH_rtofs(t=0) + ADT + datum_offset
-        (ADT datum offset, e.g. -0.45 is applied here), in float32 as ncap2. Where either
+        (ADT already carries adt_datum_offset, applied per day on read/regrid), in float32 as ncap2. Where either
         side is invalid the point is left as the -30000 fill (dry in gen_3Dth).
         """
         try:
@@ -577,7 +594,7 @@ def _apply_adt_blend(self, ssh_path: Path, adt_field, work_dir: Path,
             shutil.copy2(ssh_path, output)
 
             ds = Dataset(str(output), "r+")
-            ssh = ds.variables["ssh"][...] + (self.config.adt_datum_offset or 0)
+            ssh = ds.variables["ssh"]
             nt, ny, nx = ssh.shape
 
             dst_lon = np.array(ds.variables["xlon"][:], dtype=np.float64)

nos_utils/forcing/nudging.py

@@ -742,7 +742,7 @@ def _process_python(self) -> ForcingResult:
         # time records (rare, but the reader supports it), they are
         # treated as consecutive 6h steps starting at the file's stamped
         # forecast hour.
-        rtofs_dt_input = self.config.nudging_dt
+        rtofs_dt_input = self.config.nudging_dt if getattr(self.config, "rtofs_daily_3z", False) else 21600.0
 
         is_split = getattr(self.config, "rtofs_daily_3z", False)
         if is_split:
@@ -899,8 +899,8 @@ def get_time_key(fp): return fp.name.split('_daily')[0]
         # the model t=0 anchor (= cycle - nowcast_hours).
         target_dt = getattr(self.config, "nudging_dt", 10800.0)
         if self.config.obc_ops_timeline:
-            # Ops phase-anchored records, time 0 = phase start.
-            target_dt = self.config.nudging_dt
+            # step_nu_tr = 21600 s: phase-anchored 6-hourly records, time 0 = phase start. MJ (10/02/26)
+            target_dt = 21600.0
 
         # The output time axis covers [0, sim_duration] at ``target_dt``
         # cadence — i.e. 0, 3h, 6h, ..., nowcast+forecast hours from the

nos_utils/forcing/rtofs.py

@@ -1019,6 +1019,7 @@ def get_time_key(fp): return fp.name.split('_daily')[0]
             file_groups = [[f] for f in files_3d]
 
         for group in file_groups:
+            datasets = []
             try:
                 datasets = [Dataset(str(f)) for f in group]
                 # Merge variables from multiple files into one dictionary
@@ -1030,7 +1031,6 @@ def get_time_key(fp): return fp.name.split('_daily')[0]
                     subset.update(self._stofs_subset_roi(ds, roi, ["temperature", "salinity", "u", "v"]))
 
                 if "temperature" not in subset:
-                    for ds in datasets: ds.close()
                     continue
 
                 for t in range(subset["temperature"].shape[0]):
@@ -1043,12 +1043,13 @@ def get_time_key(fp): return fp.name.split('_daily')[0]
                 if all_lon is None:
                     all_lon = np.array(subset["Longitude"])
                     all_lat = np.array(subset["Latitude"])
-                if all_depth is None and "Depth" in ds.variables:
-                    all_depth = np.array(ds.variables["Depth"][:])
-
-                for ds in datasets: ds.close()
+                if all_depth is None and "Depth" in datasets[0].variables:
+                    all_depth = np.array(datasets[0].variables["Depth"][:])
             except Exception as e:
                 log.warning(f"Failed to read 3D group {group[0].name}: {e}")
+            finally:
+                for ds in datasets:
+                    ds.close()
 
         if not all_temp:
             return None
@@ -1296,11 +1297,11 @@ def _slice_phase_records(self, path: Path) -> None:
         model_t0, _, dur = self._get_output_window(self.phase)
         cycle = datetime.strptime(self.config.pdy, "%Y%m%d") + timedelta(hours=self.config.cyc)
         start = (model_t0 - (cycle - timedelta(hours=self.config.nowcast_hours))).total_seconds()
-        n = int(np.ceil(dur / OPS_OBC_DT)) + 1
+        n = int(np.ceil(dur / self.config.obc_dt)) + 1
         tmp = path.with_name(path.name + ".phase")
         with Dataset(str(path)) as src:
             t = np.asarray(src.variables["time"][:], dtype=np.float64)
-            dt = float(np.diff(t).min()) if len(t) > 1 else OPS_OBC_DT
+            dt = float(np.diff(t).min()) if len(t) > 1 else self.config.obc_dt
             idx = np.clip(np.rint(start / dt).astype(int) + np.arange(n), 0, len(t) - 1)
             if start / dt + n - 1 > len(t) - 1:
                 self._ops_note(f"{path.name}: the exe series ends before the {n}-record "
@@ -1586,7 +1587,7 @@ def _ops_slot_times(self) -> List[datetime]:
         t0 = cycle_dt - timedelta(hours=int(self.config.nowcast_hours))
         span = (cycle_dt - t0).total_seconds() + (
             int(self.config.forecast_hours) + self.buffer_hours) * 3600.0
-        n = max(self.config.obc_min_timesteps, int(np.ceil(span / self.config.obc_dt)) + 1)
+        n = max(self.config.obc_min_slots, int(np.ceil(span / self.config.obc_dt)) + 1)
         return [t0 + timedelta(seconds=self.config.obc_dt * k) for k in range(n)]
 
     def _ops_candidates(self, date: datetime, kind: str):
@@ -1655,10 +1656,10 @@ def _find_ops_files(self) -> Tuple[List[Path], List[Path]]:
             if have and n_prev[kind]:
                 self._select_notes.append(
                     f"RTOFS {kind}: {n_prev[kind]} of {len(have)} slots filled from the previous cycle")
-            miss = [vt for vt in slots[:self.config.obc_min_timesteps] if vt not in picked[kind]] if have else []
+            miss = [vt for vt in slots[:self.config.obc_min_slots] if vt not in picked[kind]] if have else []
             if miss:
                 self._select_notes.append(
-                    f"RTOFS {kind}: {len(miss)} of the first {self.config.obc_min_timesteps} slots missing "
+                    f"RTOFS {kind}: {len(miss)} of the first {self.config.obc_min_slots} slots missing "
                     f"in both cycles (first {miss[0]:%Y-%m-%d %Hz})")
         for m in self._select_notes:
             log.warning(m)

nos_utils/forcing/surface_restore.py

@@ -43,6 +43,12 @@ def __init__(
         self.weight_file = weight_file
         self.phase = phase
 
+    def find_input_files(self) -> List[Path]:
+        """Abstract in the base class and was never implemented, so the class could not be instantiated. MJ (10/24/26)"""
+        from .rtofs import RTOFSProcessor
+        return RTOFSProcessor(self.config, self.input_path, self.output_path,
+                              phase=self.phase).find_input_files_by_type()[1]
+
     def process(self) -> ForcingResult:
         if not HAS_NETCDF4:
             return ForcingResult(
@@ -56,6 +62,15 @@ def process(self) -> ForcingResult:
                 warnings=["Surface restore disabled in config"],
             )
 
+        import os
+        if os.environ.get("NOS_SURFACE_RESTORE_EXPERIMENTAL") != "1":
+            msg = ("Surface restore is enabled but this processor is untested and known incomplete "
+                   "(T and S read from one file, boundary instead of restore nodes, placeholder time axis); "
+                   "set NOS_SURFACE_RESTORE_EXPERIMENTAL=1 to run it anyway. MJ (10/24/26)")
+            log.error(msg)
+            return ForcingResult(success=False, source=self.SOURCE_NAME, errors=[msg])
+        log.warning("Surface restore running in experimental mode: output is unvalidated. MJ (10/24/26)")
+
         self.create_output_dir()
 
         # Reuse RTOFS 3D discovery to get the raw daily/6-hourly files