Enable STOFS Pacific System

nos_utils/config.py

@@ -172,11 +172,11 @@ class ForcingConfig:
 
     # --- OBC (ocean boundary) settings ---
     # RTOFS ROI indices for 2D (ssh) extraction
-    obc_roi_2d: Optional[dict] = None  # {x1, x2, y1, y2}
+    obc_roi_2d: Optional[dict] = None  # {x1, x2, y1, y2, x_stride, y_stride}
     # RTOFS ROI indices for 3D (T,S,U,V) extraction
-    obc_roi_3d: Optional[dict] = None  # {x1, x2, y1, y2}
+    obc_roi_3d: Optional[dict] = None  # {x1, x2, y1, y2, x_stride, y_stride}
     # RTOFS ROI indices for nudging (slightly larger domain than OBC)
-    nudge_roi_3d: Optional[dict] = None  # {x1, x2, y1, y2}
+    nudge_roi_3d: Optional[dict] = None  # {x1, x2, y1, y2, x_stride, y_stride}
     # Which regional RTOFS 3dz tile to read for 3D T/S/U/V (OBC and nudging
     # share this discovery step). Ops stages "alaska", "US_east" and
     # "US_west" tiles side by side for every valid time with no region in
@@ -192,6 +192,8 @@ class ForcingConfig:
     # for every system, not just Alaska, so correctness stops depending on
     # alphabetical luck.
     rtofs_3d_region: Optional[str] = None
+    # Whether to use global daily 3D files (split by variable) instead of regional 6-hourly.
+    rtofs_daily_3z: bool = False
     # SSH offset applied to boundary elevation (meters)
     # Geoid-to-MSL datum offset (meters). OFS-specific — verify from Fortran source.
     # SECOFS: 1.25 (confirmed). STOFS-3D-ATL: 0.04. Other OFS: 0.0 until verified.
@@ -226,17 +228,30 @@ class ForcingConfig:
     # sliced per phase, non-zero uv3D. Off keeps the previous-cycle search, 3-hourly T/S and zero
     # uv3D that SECOFS, Pacific and Alaska use. MJ (10/02/26)
     obc_ops_timeline: bool = False
+    # Temporal resolution of the output boundary files (seconds).
+    obc_dt: float = 21600.0
+    # Apply 2-step temporal moving average to OBC files (STOFS-3D-PAC ops parity).
+    obc_time_smoothing: bool = False
     # T,S for boundary nodes outside the RTOFS grid (stofs_3d_atl_obc_3dth_nc.in). MJ (10/01/26)
     obc_tem_outside: float = 20.0
     obc_sal_outside: float = 33.0
     # ADT satellite SSH blending (STOFS-3D-ATL uses CMEMS ADT to correct boundary SSH)
     adt_enabled: bool = False
+    adt_datum_offset: float = -0.45
     # 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
+    # Apply 2-step temporal moving average to nudging fields.
+    nudging_time_smoothing: bool = False
+    # Surface restore enabled and smoothing
+    surface_restore_enabled: bool = False
+    surface_restore_time_smoothing: bool = False
+    surface_restore_weight_file: Optional[Path] = None
     # Number of vertical levels
     n_levels: int = 51
     # SCHISM model timestep (seconds) used to write the elev2D.th.nc
@@ -577,6 +592,7 @@ def for_stofs_3d_atl(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig":
             # carry only the 778 elevation nodes SCHISM expects.
             obc_elev_segments=[0, 1],
             adt_enabled=True,
+            adt_datum_offset=-0.45,
             adt_weight_file=Path("stofs_3d_atl_ufs.adt_weight.nc"),
             obc_ssh_hold_first_record=True,
             obc_ops_timeline=True,
@@ -641,6 +657,7 @@ def for_stofs_3d_atl_ufs(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig
             # carry only the 778 elevation nodes SCHISM expects.
             obc_elev_segments=[0, 1],
             adt_enabled=True,
+            adt_datum_offset=-0.45,
             adt_weight_file=Path("stofs_3d_atl_ufs.adt_weight.nc"),
             obc_ssh_hold_first_record=True,
             obc_ops_timeline=True,
@@ -701,11 +718,22 @@ def for_stofs_3d_pac(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig":
             # Segments 2-8 are flow-only river OBCs (excluded from *.th.nc).
             obc_elev_segments=[0, 1, 9],
             obc_ssh_offset=0.0,
-            adt_enabled=False,
+            adt_enabled=True,
+            adt_datum_offset=-0.4715517,
             nudging_enabled=True,
             nudging_timescale_seconds=86400.0,
             nwm_product="medium_range_mem1",
-            obc_min_timesteps=21,
+            obc_min_timesteps=4,
+            obc_dt=86400.0,
+            nudging_dt=86400.0,
+            obc_time_smoothing=True,
+            nudging_time_smoothing=True,
+            rtofs_daily_3z=True,
+            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,
         )
         defaults.update(overrides)
         return cls(**defaults)
@@ -730,11 +758,22 @@ def for_stofs_3d_pac_ufs(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig
             dynamic_adjust_enabled=False,
             obc_elev_segments=[0, 1, 9],
             obc_ssh_offset=0.0,
-            adt_enabled=False,
+            adt_enabled=True,
+            adt_datum_offset=-0.4715517,
             nudging_enabled=True,
             nudging_timescale_seconds=86400.0,
             nwm_product="medium_range_mem1",
-            obc_min_timesteps=21,
+            obc_min_timesteps=4,
+            obc_dt=86400.0,
+            nudging_dt=86400.0,
+            obc_time_smoothing=True,
+            nudging_time_smoothing=True,
+            rtofs_daily_3z=True,
+            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,
             # 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,
@@ -956,6 +995,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
 
         # 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
@@ -1011,6 +1051,7 @@ def from_yaml(cls, yaml_path, **overrides) -> "ForcingConfig":
             datm_hrrr_weight=_hrrr_weight,
             hrrr_small_grib=_hrrr_sg,
             adt_enabled=adt.get("enabled", False) if isinstance(adt, dict) else False,
+            adt_datum_offset=_adt_off,
             adt_weight_file=_adt_wt,
             nwm_product=nwm_product,
             river_ops_static=_river_ops,
@@ -1032,6 +1073,18 @@ def from_yaml(cls, yaml_path, **overrides) -> "ForcingConfig":
             kwargs["nudge_roi_3d"] = {k: int(v) for k, v in nudge_roi.items()}
         if rtofs_3d_region:
             kwargs["rtofs_3d_region"] = str(rtofs_3d_region)
+
+        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_time_smoothing"] = _strict_bool(nudge.get("time_smoothing"), "nudging.time_smoothing", False)
+
+        surf_res = ocean.get("surface_restore", {}) if isinstance(ocean, 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"):
+            kwargs["surface_restore_weight_file"] = Path(surf_res["weight_file"])
 
         # NWM target counts
         if isinstance(river, dict):
@@ -1124,7 +1177,10 @@ def from_yaml(cls, yaml_path, **overrides) -> "ForcingConfig":
         # Non-STOFS OFSes have shorter OBC windows (typically ~10 records for
         # a 54h run) and would be wrongly flagged if this ran against them;
         # keep the dataclass default of 0 (disabled) for non-STOFS.
-        if is_stofs_like:
+        yaml_min_ts = obc.get("min_timesteps")
+        if yaml_min_ts is not None:
+            kwargs["obc_min_timesteps"] = int(yaml_min_ts)
+        elif is_stofs_like:
             kwargs["obc_min_timesteps"] = 21
 
         # Declarative prep extras override. When a root-level `prep.extras`

nos_utils/coords.py

@@ -71,6 +71,27 @@ def normalize_lon(
     return lon
 
 
+def shift_grid_0360(
+    lon: np.ndarray, data: "Optional[np.ndarray]" = None, lon_axis: int = -1
+) -> "Tuple[np.ndarray, Optional[np.ndarray]]":
+    """Shift a [-180, 180] grid and optional data array to [0, 360].
+
+    Useful for stitching global [-180, 180] fields (like ADT) across the 
+    antimeridian to match Pacific domain bounds smoothly.
+    """
+    lon = np.asarray(lon)
+    lon_360 = np.where(lon < 0.0, lon + 360.0, lon)
+
+    sort_idx = np.argsort(lon_360)
+    lon_out = lon_360[sort_idx]
+
+    if data is not None:
+        data_out = np.take(data, sort_idx, axis=lon_axis)
+        return lon_out, data_out
+
+    return lon_out, None
+
+
 def lon_convention(config: "object") -> str:
     """Return the longitude convention used by a ``ForcingConfig``.
 

nos_utils/forcing/adt.py

@@ -5,12 +5,12 @@
 boundary condition accuracy for STOFS-3D-ATL.
 
 The core formula (ops stofs_3d_atl_create_obc_3d_th_non_adjust.sh, SSH_1.nc block):
-    SSH_final = SSH_rtofs - SSH_rtofs(t=0) + (ADT - 0.45)
+    SSH_final = SSH_rtofs - SSH_rtofs(t=0) + (ADT + adt_datum_offset)
 
 ADT is the mean of the valid-day and previous-day CMEMS files (ncra), mapped onto
 the SSH_1 grid with the ops ESMF map (ncremap with stofs_3d_atl_adt_weight.nc,
 which is a nearest-source map with S=1 despite its "bilinear" global attribute).
-The -0.45 m datum offset (stofs_3d_atl_adt_cvtz.nco) is applied per day and rounded
+The datum offset (stofs_3d_atl_adt_cvtz.nco) is applied per day and rounded
 to float32 before the mean, as ops does. Without the map, a bilinear regrid is used
 and the result is not ops-exact.
 
@@ -50,10 +50,6 @@
 
 # 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)
@@ -147,10 +143,9 @@ def blend_ssh(self, ssh_path: Path, work_dir: Path) -> Optional[Path]:
             log.warning(f"ADT blending failed: {e}")
             return None
 
-    @staticmethod
-    def _dest_bounds(ssh_path: 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 (same -180/180 convention as ADT), else the ops ROI. MJ (10/01/26)
+        by ADT_PAD, else the config domain bounds. MJ (10/01/26)
         """
         try:
             with Dataset(str(ssh_path)) as ds:
@@ -159,8 +154,9 @@ def _dest_bounds(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 ops ADT ROI")
-            return ADT_LON_MIN, ADT_LON_MAX, ADT_LAT_MIN, ADT_LAT_MAX
+            log.warning(f"Cannot read SSH_1 grid extent ({e}) — using config domain bounds")
+            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 _find_adt_data(self) -> Optional[Path]:
         """Newest available CMEMS ADT file (valid day, else previous day)."""
@@ -218,23 +214,29 @@ def _warn(self, msg):
         log.warning(msg)
         self.warnings.append(msg)
 
-    @staticmethod
-    def _subset_ops_roi(adt_path: Path):
+    def _subset_ops_roi(self, adt_path: Path, bounds):
         """ncks -d longitude,-62.5,-51.5 -d latitude,7.0,54.0: inclusive coordinate bounds.
 
         Returns (adt[lat, lon] in m with NaN for missing, lons, lats, n_time). MJ (10/02/26)
         """
+        lon_min, lon_max, lat_min, lat_max = bounds
         with Dataset(str(adt_path)) as ds:
             lon_name = "longitude" if "longitude" in ds.variables else "lon"
             lat_name = "latitude" if "latitude" in ds.variables else "lat"
             lons = np.array(ds.variables[lon_name][:], dtype=np.float64)
             lats = np.array(ds.variables[lat_name][:], dtype=np.float64)
-            li = np.where((lons >= ADT_LON_MIN) & (lons <= ADT_LON_MAX))[0]
-            la = np.where((lats >= ADT_LAT_MIN) & (lats <= ADT_LAT_MAX))[0]
-            if li.size == 0 or la.size == 0:
-                raise ValueError("ADT file does not cover the ops ROI")
             var = ds.variables["adt" if "adt" in ds.variables else "surf_el"]
-            a = np.ma.filled(var[..., la[0]:la[-1] + 1, li[0]:li[-1] + 1], np.nan)
+            a = var[...]
+
+            from ..coords import lon_convention, shift_grid_0360
+            if lon_convention(self.config) == "0360":
+                lons, a = shift_grid_0360(lons, a, lon_axis=-1)
+
+            li = np.where((lons >= lon_min) & (lons <= lon_max))[0]
+            la = np.where((lats >= lat_min) & (lats <= lat_max))[0]
+            if li.size == 0 or la.size == 0:
+                raise ValueError("ADT file does not cover the domain ROI")
+            a = np.ma.filled(a[..., la[0]:la[-1] + 1, li[0]:li[-1] + 1], np.nan)
             a = np.asarray(a, dtype=np.float64)
             if a.ndim == 3:
                 a = a[0]
@@ -273,7 +275,7 @@ def _apply_esmf_map(S, row, col, n_b, src):
         return np.where(tally > 0, num, np.nan)
 
     def _regrid_esmf(self, adt_files, ssh_path: Path):
-        """Ops ADT step with the ESMF map: ROI subset, ncremap, float32(adt-0.45), two-day mean.
+        """Ops ADT step with the ESMF map: ROI subset, ncremap, float32(adt+datum_offset), two-day mean.
 
         Each day is rounded to float32 after the offset and the valid days are averaged
         (float64 accumulation, float32 result), as ncap2 and ncra do. Returns the float32
@@ -307,9 +309,10 @@ 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)
             days = []
             for f in adt_files:
-                a, lons, lats = self._subset_ops_roi(f)
+                a, lons, lats = self._subset_ops_roi(f, bounds)
                 if a.size != n_a:
                     raise ValueError(f"ADT ROI subset has {a.size} cells, map n_a={n_a}")
                 if src_dims and tuple(src_dims) != (a.shape[1], a.shape[0]):
@@ -320,7 +323,7 @@ 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 - 0.45).astype(np.float32) if self.ops_numerics else out - 0.45)
+                days.append(out.astype(np.float32) if self.ops_numerics else out)
             if self.ops_numerics:
                 stack = np.stack(days)
                 ok = np.isfinite(stack)
@@ -343,7 +346,7 @@ def dlon_diff(a, b):
         return field
 
     def _write_adt_on_rtofs(self, adt_dst, dst_lon, dst_lat, out: Path):
-        """adt_aft_cvtz_cln.nc analogue: float32 surf_el(time=1, ylat, xlon) in m after -0.45.
+        """adt_aft_cvtz_cln.nc analogue: float32 surf_el(time=1, ylat, xlon) in m after datum offset.
 
         Missing cells carry the ops fill (the CMEMS int fill held as float32).
         """
@@ -421,13 +424,12 @@ def _load_previous(self, ssh_path: Path):
         return None
 
     def _read_adt(self, adt_path: Path, bounds=None):
-        """Read and subset ADT to bounds, offset by -0.45 m.
+        """Read and subset ADT to bounds, offset by the user specified amount.
 
         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 (
-            ADT_LON_MIN, ADT_LON_MAX, ADT_LAT_MIN, ADT_LAT_MAX)
+        lon_min, lon_max, lat_min, lat_max = bounds or self._dest_bounds(None)
         try:
             ds = Dataset(str(adt_path))
             try:
@@ -436,18 +438,19 @@ 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][...]
+
+                from ..coords import lon_convention, shift_grid_0360
+                if lon_convention(self.config) == "0360":
+                    lons, adt_data = shift_grid_0360(lons, adt_data, lon_axis=-1)
 
                 lon_idx = np.where((lons >= lon_min) & (lons <= lon_max))[0]
                 lat_idx = np.where((lats >= lat_min) & (lats <= lat_max))[0]
 
                 if len(lon_idx) == 0 or len(lat_idx) == 0:
                     log.warning("ADT data doesn't cover target domain")
                     return None
-
-                adt_var = "adt" if "adt" in ds.variables else "surf_el"
-                # 0.45 from NCO file
-                # TODO: Should be different for different systems, e.g. PAC
-                adt_data = ds.variables[adt_var][...] - 0.45
                 sl_y = slice(lat_idx[0], lat_idx[-1] + 1)
                 sl_x = slice(lon_idx[0], lon_idx[-1] + 1)
                 if adt_data.ndim == 3:
@@ -564,8 +567,8 @@ def _apply_adt_blend(self, ssh_path: Path, adt_field, work_dir: Path,
                          adt_dst=None) -> Optional[Path]:
         """Apply the ops ADT blend to SSH_1.nc.
 
-        Formula: SSH_final(t) = SSH_rtofs(t) - SSH_rtofs(t=0) + ADT
-        (ADT already carries the -0.45 offset), in float32 as ncap2. Where either
+        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
         side is invalid the point is left as the -30000 fill (dry in gen_3Dth).
         """
         try:
@@ -574,7 +577,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"]
+            ssh = ds.variables["ssh"][...] + (self.config.adt_datum_offset or 0)
             nt, ny, nx = ssh.shape
 
             dst_lon = np.array(ds.variables["xlon"][:], dtype=np.float64)