STOFS-2D-Global ops mode: one 180 h forecast forcing segment (hourly, 121-180 h interpolated)

nos_utils/forcing/adcirc_met.py

@@ -594,7 +594,7 @@ def _write_manifest(path: Path, acquired: List[AcquiredFile], phase: str,
 # Replicates exstofs_2d_glo_gfs.sh / stofs_2d_glo_surface_forcing.sh / stofs_2d_glo_getges.sh:
 # one pass per record, ncks -v + ncecat + ncwa -a time as netCDF4 calls. MJ (10/06/26)
 
-SFC_SEGMENTS = ("ncst", "fcst1", "fcst2")
+SFC_SEGMENTS = ("ncst", "fcst")
 SFC_OUTPUTS = (("221", ("pressfc",)), ("222", ("ugrd10m", "vgrd10m")), ("225", ("icec",)))
 _GETGES_FHEND = 384
 _SFCF_SETTLED_S = 900  # an sfcf file this old without a logf marker is a retrospective one. MJ (10/06/26)
@@ -622,25 +622,33 @@ def getges_sfcf(comin_gfs, valid: datetime) -> Path:
 
 
 def sfcf_valid_times(cycle: datetime, segment: str, start: Optional[datetime] = None) -> List[datetime]:
-    """Record times of the ops calls: surface1 ncst (start..cycle) and fcst1 (cycle..+120 h) hourly,
-    surface3 fcst2 (+120 h..+180 h) every 3 h. ``start`` is the ncst time_beg (multistart)."""
+    """Hourly record times: ncst (start..cycle) and the single fcst (cycle..+180 h). ``start`` is the
+    ncst time_beg (multistart). After f120 GFS is 3-hourly, so fcst fills the hours between by linear
+    interpolation, which is what ADCIRC does between the fcst2 snapshots. MJ (10/11/26)"""
     if segment not in SFC_SEGMENTS:
         raise ValueError("segment must be one of {}, got {!r}".format(SFC_SEGMENTS, segment))
     if segment == "ncst":
-        beg, end, step = start or cycle - timedelta(hours=6), cycle, 1
-    elif segment == "fcst1":
-        beg, end, step = cycle, cycle + timedelta(hours=120), 1
+        beg, end = start or cycle - timedelta(hours=6), cycle
     else:
-        beg, end, step = cycle + timedelta(hours=120), cycle + timedelta(hours=180), 3
+        beg, end = cycle, cycle + timedelta(hours=180)
     if beg > end:
         raise MetForcingError("sfcf window starts after it ends: {} > {}".format(beg, end))
     out = []
     while beg <= end:
         out.append(beg)
-        beg += timedelta(hours=step)
+        beg += timedelta(hours=1)
     return out
 
 
+def _interp_bracket(cycle: datetime, valid: datetime) -> Optional[Tuple[int, int, float]]:
+    """(lo, hi, weight) forecast hours when ``valid`` falls between two 3-hourly GFS records, else None. MJ (10/11/26)"""
+    fh = int((valid - cycle).total_seconds() // 3600)
+    if fh <= 120 or fh % 3 == 0:
+        return None
+    lo = fh // 3 * 3
+    return lo, lo + 3, (fh - lo) / 3.0
+
+
 def _sfcf_ready(path: Path) -> bool:
     """Non-empty and, unless old enough to be complete, its logf marker exists (GFS writes sfcf straight into COM)."""
     if not (path.is_file() and path.stat().st_size > 0):
@@ -663,13 +671,13 @@ def _wait_for(path: Path, deadline: float, poll_s: float, sleep, clock) -> bool:
 def build_sfcf_forcing(comin_gfs, cycle: datetime, segment: str, out_dir, prefix: str = "stofs_2d_glo",
                        start: Optional[datetime] = None, wait_s: float = 0.0, poll_s: float = 10.0,
                        sleep=time.sleep, clock=time.monotonic, now=time.time) -> Dict[str, Path]:
-    """Write ``<prefix>_<segment>.{221,222,225}.nc`` from the GFS sfcf files, as the ops GFS_NCST/FCST1/FCST2 jobs.
+    """Write ``<prefix>_<segment>.{221,222,225}.nc`` from the GFS sfcf files (ncst, or one hourly 180 h fcst).
 
     Valid times at or before ``cycle`` use getges; later ones the current cycle's own sfcf file, waited for
     up to ``wait_s`` seconds in total. Output files appear under their final names only when all three are
     complete, 221 last (the ops skip test is on 221). Records are copied one at a time. In real time (within
     12 h of ``cycle``) the cycle-time record first waits for f000 so getges does not fall back to f006 of the
-    previous cycle.
+    previous cycle. MJ (10/11/26)
     """
     import netCDF4 as nc
 
@@ -684,29 +692,54 @@ def build_sfcf_forcing(comin_gfs, cycle: datetime, segment: str, out_dir, prefix
     done = False
     try:
         for rec, valid in enumerate(valids):
-            if valid <= cycle:
+            br = _interp_bracket(cycle, valid) if segment == "fcst" else None
+            if br:
+                srcs = [sfcf_path(comin_gfs, cycle, h) for h in br[:2]]
+                for sp in srcs:
+                    if not _wait_for(sp, deadline, poll_s, sleep, clock):
+                        raise FileNotAvailableError("GFS file did not appear before the wait expired: {}".format(sp))
+                src = srcs[0]
+            elif valid <= cycle:
                 if valid == cycle and now() - calendar.timegm(cycle.timetuple()) < _REALTIME_S:
                     _wait_for(sfcf_path(comin_gfs, cycle, 0), deadline, poll_s, sleep, clock)
                 src = getges_sfcf(comin_gfs, valid)
             else:
                 src = sfcf_path(comin_gfs, cycle, int((valid - cycle).total_seconds() // 3600))
                 if not _wait_for(src, deadline, poll_s, sleep, clock):
                     raise FileNotAvailableError("GFS file did not appear before the wait expired: {}".format(src))
-            used.append((valid, src))
+            used.extend((valid, p) for p in (srcs if br else [src]))  # interpolated records list both brackets. MJ (10/11/26)
             log.info("record %d valid %s <- %s", rec, valid, src)
             with nc.Dataset(str(src)) as ds:
                 ds.set_auto_maskandscale(False)
                 if rec == 0:
                     for n, names in SFC_OUTPUTS:
                         outs[n] = _open_sfcf_output(nc, part[n], ds, names)
                 lat, lon = ds.variables["lat"][:], ds.variables["lon"][:]
-                for n, names in SFC_OUTPUTS:
-                    o = outs[n]
-                    if lat.shape != (len(o.dimensions["grid_yt"]), len(o.dimensions["grid_xt"])):  # ncecat needs equal shapes. MJ (10/06/26)
-                        raise MetForcingError("{}: grid differs from the first record".format(src))
-                    o.variables["lat"][rec], o.variables["lon"][rec] = lat, lon
-                    for v in names:
-                        o.variables[v][rec] = ds.variables[v][0]
+                hi = nc.Dataset(str(srcs[1])) if br else None
+                try:
+                    if hi is not None:
+                        hi.set_auto_maskandscale(False)
+                        if not (np.array_equal(hi.variables["lat"][:], lat) and np.array_equal(hi.variables["lon"][:], lon)):
+                            raise MetForcingError("{}: grid differs from {}".format(srcs[1], src))
+                    for n, names in SFC_OUTPUTS:
+                        o = outs[n]
+                        if lat.shape != (len(o.dimensions["grid_yt"]), len(o.dimensions["grid_xt"])):  # ncecat needs equal shapes. MJ (10/06/26)
+                            raise MetForcingError("{}: grid differs from the first record".format(src))
+                        o.variables["lat"][rec], o.variables["lon"][rec] = lat, lon
+                        for v in names:
+                            a0 = ds.variables[v][0]
+                            if hi is None:
+                                o.variables[v][rec] = a0
+                                continue
+                            a1 = hi.variables[v][0]
+                            r = a0.astype(np.float64) * (1.0 - br[2]) + a1.astype(np.float64) * br[2]
+                            fv = ds.variables[v].getncattr("_FillValue") if "_FillValue" in ds.variables[v].ncattrs() else None
+                            if fv is not None:
+                                r = np.where((a0 == fv) | (a1 == fv), fv, r)
+                            o.variables[v][rec] = r.astype(a0.dtype)
+                finally:
+                    if hi is not None:
+                        hi.close()
         for o in outs.values():
             o.close()
         outs = {}

tests/test_adcirc_met.py

@@ -430,13 +430,11 @@ def test_sfcf_windows():
     ncst = am.sfcf_valid_times(CYC, "ncst")
     assert (len(ncst), ncst[0], ncst[-1]) == (7, CYC - timedelta(hours=6), CYC)
     assert len(am.sfcf_valid_times(CYC, "ncst", CYC - timedelta(hours=12))) == 13
-    f1 = am.sfcf_valid_times(CYC, "fcst1")
-    assert (len(f1), f1[0], f1[-1]) == (121, CYC, CYC + timedelta(hours=120))
-    f2 = am.sfcf_valid_times(CYC, "fcst2")
-    assert (len(f2), f2[0], f2[1], f2[-1]) == (
-        21, CYC + timedelta(hours=120), CYC + timedelta(hours=123), CYC + timedelta(hours=180))
-    with pytest.raises(ValueError):
-        am.sfcf_valid_times(CYC, "spinup")
+    f = am.sfcf_valid_times(CYC, "fcst")
+    assert (len(f), f[0], f[-1]) == (181, CYC, CYC + timedelta(hours=180))
+    for old in ("spinup", "fcst1", "fcst2"):
+        with pytest.raises(ValueError):
+            am.sfcf_valid_times(CYC, old)
     with pytest.raises(am.MetForcingError):
         am.sfcf_valid_times(CYC, "ncst", CYC + timedelta(hours=1))
 
@@ -488,23 +486,55 @@ def test_ncst_start_override_gives_fewer_records(tmp_path):
     assert _read(files["225"], "icec").shape[0] == 6
 
 
-def test_fcst1_and_fcst2_record_counts_and_stride(tmp_path):
+def _fcst_tank(tank, last=180):
+    _fill_tank(tank, [(CYC, h) for h in range(0, 121)] + [(CYC, h) for h in range(123, last + 1, 3)])
+
+
+def test_fcst_is_hourly_with_interpolation_after_f120(tmp_path):
+    tank = tmp_path / "g"
+    _fcst_tank(tank)
+    f = am.build_sfcf_forcing(tank, CYC, "fcst", tmp_path / "o")
+    p = _read(f["221"], "pressfc")
+    assert p.shape[0] == 181
+    s = lambda h, k: _valid_seed(CYC, h) * 10 + k
+    for h in (0, 7, 120, 123, 126, 180):  # native records are exact copies. MJ (10/11/26)
+        assert p[h][0, 0] == s(h, 2)
+    # 121 and 122 are 1/3 and 2/3 of the way from f120 to f123; the seeds step 10 per hour. MJ (10/11/26)
+    assert p[121][0, 0] == pytest.approx(s(120, 2) + 10.0)
+    assert p[122][0, 0] == pytest.approx(s(120, 2) + 20.0)
+    q = _read(f["222"], "ugrd10m")
+    assert q[121][1, 2] == pytest.approx(s(120, 0) + 10.0 + 6)
+    assert _read(f["225"], "icec").shape[0] == 181
+    man = json.loads((tmp_path / "o" / "stofs_2d_glo_fcst.sfcf_manifest.json").read_text())["files"]
+    assert len(man) == 181 + 40 and [e["source_location"] for e in man[121:123]] == [
+        str(am.sfcf_path(tank, CYC, 120)), str(am.sfcf_path(tank, CYC, 123))]
+
+
+def test_fcst_interpolation_keeps_fill_values(tmp_path):
+    tank = tmp_path / "g"
+    _fcst_tank(tank)
+    import netCDF4 as nc
+    with nc.Dataset(str(am.sfcf_path(tank, CYC, 123)), "a") as d:
+        d.variables["icec"][0, 0, 0] = np.float32(9.99e20)
+    f = am.build_sfcf_forcing(tank, CYC, "fcst", tmp_path / "o")
+    i = _read(f["225"], "icec")
+    assert i[121][0, 0] == np.float32(9.99e20) and i[120][0, 0] < 1e6 and i[121][0, 1] < 1e6
+
+
+def test_fcst_rejects_a_bracket_on_another_grid(tmp_path):
+    import netCDF4 as nc
     tank = tmp_path / "g"
-    _fill_tank(tank, [(CYC, 0)] + [(CYC, h) for h in range(1, 121)])
-    f1 = am.build_sfcf_forcing(tank, CYC, "fcst1", tmp_path / "o")
-    p = _read(f1["221"], "pressfc")
-    assert p.shape[0] == 121
-    assert p[7][0, 0] == _valid_seed(CYC, 7) * 10 + 2
-    _fill_tank(tank, [(CYC, h) for h in range(120, 181, 3)])
-    f2 = am.build_sfcf_forcing(tank, CYC, "fcst2", tmp_path / "o")
-    q = _read(f2["222"], "ugrd10m")
-    assert q.shape[0] == 21
-    assert q[1][0, 0] == _valid_seed(CYC, 123) * 10 + 0
+    _fcst_tank(tank)
+    with nc.Dataset(str(am.sfcf_path(tank, CYC, 123)), "a") as d:
+        d.variables["lon"][:] = d.variables["lon"][:] + 0.5
+    with pytest.raises(am.MetForcingError, match="grid differs"):
+        am.build_sfcf_forcing(tank, CYC, "fcst", tmp_path / "o")
+    assert list((tmp_path / "o").glob("*.nc*")) == []
 
 
 def test_forecast_waits_then_fails_clearly(tmp_path):
     tank = tmp_path / "g"
-    _fill_tank(tank, [(CYC, h) for h in range(120, 178, 3)])  # f180 is missing
+    _fcst_tank(tank, last=177)  # f180 is missing. MJ (10/11/26)
     clock = {"t": 0.0}
     naps = []
 
@@ -513,7 +543,7 @@ def sleep(s):
         clock["t"] += s
 
     with pytest.raises(am.FileNotAvailableError, match="sfcf180"):
-        am.build_sfcf_forcing(tank, CYC, "fcst2", tmp_path / "o", wait_s=25, poll_s=10,
+        am.build_sfcf_forcing(tank, CYC, "fcst", tmp_path / "o", wait_s=25, poll_s=10,
                               sleep=sleep, clock=lambda: clock["t"])
     assert naps == [10, 10, 5]
     assert list((tmp_path / "o").glob("*.nc*")) == []  # no final or .partial files from a failed run
@@ -522,17 +552,18 @@ def sleep_and_land(s):
         _mk_sfcf(tank, CYC, 180, 1)
         sleep(s)
 
-    files = am.build_sfcf_forcing(tank, CYC, "fcst2", tmp_path / "o", wait_s=25, poll_s=10,
+    files = am.build_sfcf_forcing(tank, CYC, "fcst", tmp_path / "o", wait_s=25, poll_s=10,
                                   sleep=sleep_and_land, clock=lambda: clock["t"])
-    assert _read(files["221"], "pressfc").shape[0] == 21
+    assert _read(files["221"], "pressfc").shape[0] == 181
 
 
 def test_current_cycle_f000_is_not_waited_for(tmp_path):
-    # fcst1 record 0 (valid == cycle) goes through getges, so it never waits
-    _fill_tank(tmp_path / "g", [(CYC - timedelta(hours=6), 6)] + [(CYC, h) for h in range(1, 121)])
-    files = am.build_sfcf_forcing(tmp_path / "g", CYC, "fcst1", tmp_path / "o", wait_s=0,
+    # fcst record 0 (valid == cycle) goes through getges, so it never waits. MJ (10/11/26)
+    _fill_tank(tmp_path / "g", [(CYC - timedelta(hours=6), 6)] + [(CYC, h) for h in range(1, 121)]
+                + [(CYC, h) for h in range(123, 181, 3)])
+    files = am.build_sfcf_forcing(tmp_path / "g", CYC, "fcst", tmp_path / "o", wait_s=0,
                                   sleep=lambda s: pytest.fail("waited"))
-    assert _read(files["221"], "pressfc").shape[0] == 121
+    assert _read(files["221"], "pressfc").shape[0] == 181
 
 
 def _age(p, seconds):