Enforce exact fraction closure for fractional fields

src/dapper/domains/domain.py

@@ -974,7 +974,10 @@ def export_met(
         out_dir : Path-like, optional
             Override output root. Defaults to Domain.run_dir.
         filename : str, optional
-            Optional filename prefix for output NetCDFs. If provided, each var is written to '{filename}_{var}.nc'.
+            Optional filename template for output NetCDFs. Supports two formats:
+            - Template with {var} placeholder: 'ERA5_{var}_1950-2025_z01' generates 'ERA5_TBOT_1950-2025_z01.nc'
+            - Simple prefix (legacy): 'prefix' generates 'prefix_TBOT.nc'
+            The '.nc' extension is added automatically.
         overwrite : bool
             If False, raises if MET output(s) already exist.
         clip_to_full_years : bool or None

src/dapper/landuse/landuse.py

@@ -12,6 +12,7 @@
 
 from dapper.geo import sampling
 from dapper.domains.domain import Domain
+from dapper.surf.fraction_closure import normalize_fraction_closure
 
 LonWrap = Literal["auto", "0_360", "-180_180"]
 
@@ -270,6 +271,7 @@ def _write_nc(ds_out: xr.Dataset, path: Path) -> None:
         csv_path = out_path.with_suffix(out_path.suffix + ".zonal_weights.csv")
         dfw_all.to_csv(csv_path, index=False)
 
+    out = normalize_fraction_closure(out)
     _write_nc(out, out_path)
 
     df_summary = pd.DataFrame({gid_col: order, "sample_ncells": ncells, "sample_area_total_m2": area_m2})

src/dapper/met/exporter.py

@@ -239,8 +239,10 @@ def run(self, *, pack_scope=None, filename: str | None = None, overwrite: bool =
             ``global`` for cellset outputs.
 
         filename
-            Optional filename prefix for output NetCDF files. If provided, each variable
-            is written to ``{filename}_{var}.nc``.
+            Optional filename template for output NetCDF files. Supports two formats:
+            - Template with {var} placeholder: 'ERA5_{var}_1950-2025_z01' generates 'ERA5_TBOT_1950-2025_z01.nc'
+            - Simple prefix (legacy): 'prefix' generates 'prefix_TBOT.nc'
+            The '.nc' extension is added automatically.
 
         overwrite
             If True, clears existing MET outputs before writing.
@@ -343,8 +345,17 @@ def _zone_mappings_path(self, gid: str | None = None, filename: str = "zone_mapp
         return self._met_dir_for_gid(gid) / filename
 
     def _nc_filename(self, var: str) -> str:
-        """Return the output NetCDF filename for a given variable."""
+        """Return the output NetCDF filename for a given variable.
+
+        Supports two formats:
+        - Template with {var} placeholder: 'ERA5_{var}_1950-2025_z01' -> 'ERA5_TBOT_1950-2025_z01.nc'
+        - Simple prefix (legacy): 'prefix' -> 'prefix_TBOT.nc'
+        """
         if getattr(self, "filename_prefix", None):
+            # Check if template contains {var} placeholder
+            if "{var}" in self.filename_prefix:
+                return f"{self.filename_prefix.format(var=var)}.nc"
+            # Legacy prefix mode
             return f"{self.filename_prefix}_{var}.nc"
         return f"{var}.nc"
 

src/dapper/surf/fraction_closure.py

@@ -0,0 +1,209 @@
+"""Shared fraction-closure logic for surface and landuse datasets."""
+
+from __future__ import annotations
+
+from typing import Iterable
+
+import numpy as np
+import xarray as xr
+
+# Keep closures aligned with the canonical surface variable specs.
+from dapper.surf.surface_var_specs import SURFACE_VAR_SPECS
+
+
+def _scale_to_target_sum(
+    da: xr.DataArray,
+    *,
+    dim: str,
+    target: xr.DataArray,
+    eps: float = 1e-12,
+) -> xr.DataArray:
+    """Scale values along `dim` so their sum matches `target` where possible."""
+    work = da.astype(np.float64)
+    summed = work.sum(dim=dim, skipna=True)
+    valid = np.isfinite(summed) & (np.abs(summed) > eps) & np.isfinite(target)
+    safe_den = xr.where(valid, summed, 1.0)
+    factor = xr.where(valid, target.astype(np.float64) / safe_den, 1.0)
+    return work * factor
+
+
+def _snap_partition_sum(
+    da: xr.DataArray,
+    *,
+    dim: str,
+    target: xr.DataArray,
+    eps: float = 1e-15,
+) -> xr.DataArray:
+    """Force exact closure by assigning residual to the last band in float64."""
+    work = da.astype(np.float64)
+    if dim not in work.dims or int(work.sizes.get(dim, 0)) < 1:
+        return work
+
+    head = work.isel({dim: slice(0, -1)})
+    head_sum = head.sum(dim=dim, skipna=True)
+    last_src = work.isel({dim: -1})
+    candidate_last = target.astype(np.float64) - head_sum
+
+    valid = np.isfinite(target) & np.isfinite(head_sum)
+    snapped_last = xr.where(valid, candidate_last, last_src)
+    snapped_last = xr.where(np.abs(snapped_last) < eps, 0.0, snapped_last)
+
+    out = work.copy(deep=False)
+    out[{dim: -1}] = snapped_last
+    return out
+
+
+def _close_partition(
+    ds: xr.Dataset,
+    *,
+    var_name: str,
+    dim: str,
+    target: xr.DataArray,
+) -> None:
+    if var_name not in ds:
+        return
+    da = ds[var_name]
+    if dim not in da.dims:
+        return
+    ds[var_name] = _scale_to_target_sum(da, dim=dim, target=target)
+    ds[var_name] = _snap_partition_sum(ds[var_name], dim=dim, target=target)
+
+
+def _close_unit_partition(ds: xr.Dataset, *, left: str, right: str) -> None:
+    if left not in ds or right not in ds:
+        return
+    lhs = ds[left].astype(np.float64)
+    rhs = ds[right].astype(np.float64)
+    total = lhs + rhs
+    valid = np.isfinite(total) & (total > 1e-12)
+    factor = xr.where(valid, 1.0 / total, 1.0)
+    ds[left] = lhs * factor
+    ds[right] = rhs * factor
+
+
+def _landunit_scalar_names(ds: xr.Dataset) -> list[str]:
+    # Explicitly constrained to ELM landunit classes. We gate against canonical
+    # specs so stale names are ignored automatically.
+    candidates = ["PCT_NATVEG", "PCT_CROP", "PCT_WETLAND", "PCT_LAKE", "PCT_GLACIER"]
+    return [v for v in candidates if (v in ds and v in SURFACE_VAR_SPECS)]
+
+
+def _apply_landunit_total_closure(ds: xr.Dataset) -> None:
+    """Enforce landunit totals to close to 100 (including urban aggregate)."""
+    scalar_names = _landunit_scalar_names(ds)
+    landunit_terms: list[xr.DataArray] = [ds[name].astype(np.float64) for name in scalar_names]
+
+    urban_group = None
+    if "PCT_URBAN" in ds:
+        urban = ds["PCT_URBAN"].astype(np.float64)
+        if "numurbl" in urban.dims:
+            urban_group = urban
+            landunit_terms.append(urban.sum(dim="numurbl", skipna=True))
+        else:
+            landunit_terms.append(urban)
+
+    if len(landunit_terms) < 2:
+        return
+
+    current_total = sum(landunit_terms)
+    near_100 = np.abs(current_total - 100.0) <= 1.0
+    valid = np.isfinite(current_total) & (current_total > 1e-12) & near_100
+    factor = xr.where(valid, 100.0 / current_total, 1.0)
+
+    for name in scalar_names:
+        ds[name] = ds[name].astype(np.float64) * factor
+    if urban_group is not None:
+        ds["PCT_URBAN"] = urban_group * factor
+    elif "PCT_URBAN" in ds:
+        ds["PCT_URBAN"] = ds["PCT_URBAN"].astype(np.float64) * factor
+
+    total_after = 0.0
+    for name in scalar_names:
+        total_after = total_after + ds[name].astype(np.float64)
+    if "PCT_URBAN" in ds:
+        urb = ds["PCT_URBAN"].astype(np.float64)
+        if "numurbl" in urb.dims:
+            total_after = total_after + urb.sum(dim="numurbl", skipna=True)
+        else:
+            total_after = total_after + urb
+
+    resid = 100.0 - total_after
+
+    # Prefer scalar classes first; fallback to last urban class.
+    snapped = False
+    for name in ("PCT_GLACIER", "PCT_LAKE", "PCT_WETLAND", "PCT_CROP", "PCT_NATVEG"):
+        if name in ds:
+            ds[name] = ds[name].astype(np.float64) + resid
+            snapped = True
+            break
+
+    if (not snapped) and ("PCT_URBAN" in ds):
+        urb = ds["PCT_URBAN"].astype(np.float64)
+        if "numurbl" in urb.dims and int(urb.sizes.get("numurbl", 0)) >= 1:
+            urb[{"numurbl": -1}] = urb.isel(numurbl=-1) + resid
+            ds["PCT_URBAN"] = urb
+        else:
+            ds["PCT_URBAN"] = urb + resid
+
+
+def _full_like_from_partition(ds: xr.Dataset, *, var_name: str, dim: str, value: float) -> xr.DataArray | None:
+    if var_name not in ds or dim not in ds[var_name].dims:
+        return None
+    return xr.full_like(ds[var_name].isel({dim: 0}, drop=True), value, dtype=np.float64)
+
+
+def normalize_fraction_closure(ds: xr.Dataset) -> xr.Dataset:
+    """Apply canonical fraction closure for surface/landuse datasets."""
+    ds2 = ds.copy(deep=False)
+
+    _apply_landunit_total_closure(ds2)
+
+    # Natural-patch weights should sum to 100 for every non-natpft index tuple.
+    tgt_nat = _full_like_from_partition(ds2, var_name="PCT_NAT_PFT", dim="natpft", value=100.0)
+    if tgt_nat is not None:
+        _close_partition(ds2, var_name="PCT_NAT_PFT", dim="natpft", target=tgt_nat)
+
+    if "PCT_CFT" in ds2 and "cft" in ds2["PCT_CFT"].dims and "PCT_CROP" in ds2:
+        _close_partition(
+            ds2,
+            var_name="PCT_CFT",
+            dim="cft",
+            target=ds2["PCT_CROP"].astype(np.float64),
+        )
+
+    if "PCT_GLC_MEC" in ds2 and "nglcec" in ds2["PCT_GLC_MEC"].dims and "PCT_GLACIER" in ds2:
+        _close_partition(
+            ds2,
+            var_name="PCT_GLC_MEC",
+            dim="nglcec",
+            target=ds2["PCT_GLACIER"].astype(np.float64),
+        )
+
+    # TopounitFracArea is a decimal fraction (0-1) and should sum to 1.0 over topounit.
+    tgt_top = _full_like_from_partition(ds2, var_name="TopounitFracArea", dim="topounit", value=1.0)
+    if tgt_top is not None:
+        _close_partition(ds2, var_name="TopounitFracArea", dim="topounit", target=tgt_top)
+
+    _close_unit_partition(ds2, left="FSURF", right="FGRD")
+
+    return ds2
+
+
+def closure_critical_variables(present_vars: Iterable[str]) -> set[str]:
+    """Return vars that should be written as float64 to preserve closure."""
+    present = set(present_vars)
+    critical = {
+        "PCT_NAT_PFT",
+        "PCT_CFT",
+        "PCT_GLC_MEC",
+        "PCT_NATVEG",
+        "PCT_CROP",
+        "PCT_GLACIER",
+        "PCT_WETLAND",
+        "PCT_LAKE",
+        "PCT_URBAN",
+        "TopounitFracArea",
+        "FSURF",
+        "FGRD",
+    }
+    return {v for v in critical if v in present}

src/dapper/surf/schema.py

@@ -108,6 +108,7 @@ def pdef(
     units: str = "",
     doc: str = "",
     required_level: str = "",
+    contexts: Tuple[str, ...] = (),
     **attrs,
 ) -> ParDef:
     """
@@ -128,6 +129,7 @@ def pdef(
         doc=doc,
         required_level=required_level,
         attrs=attrs or {},
+        contexts=tuple(contexts or ()),
     )
 
 

src/dapper/surf/sfile.py

@@ -12,6 +12,7 @@
 from dapper.surf import sample as SP  # for from_halfdegree_point 
 from dapper.geo import sampling  # shared gridded sampler
 from dapper.surf.surface_var_specs import SURFACE_VAR_SPECS
+from dapper.surf.fraction_closure import normalize_fraction_closure, closure_critical_variables
 
 ArrayLike = Union[np.ndarray, "xr.DataArray", float, int]
 
@@ -187,7 +188,7 @@ def write_surface_nc(
     import datetime as _dt
 
     # ---- global attrs ----
-    ds2 = ds.copy(deep=False)
+    ds2 = normalize_fraction_closure(ds)
     merged = dict(ds2.attrs)
 
     if dapper_attrs:
@@ -202,11 +203,16 @@ def write_surface_nc(
 
     ds2.attrs = merged
 
+    closure_critical = closure_critical_variables(ds2.data_vars)
     enc: Dict[str, dict] = {}
     for v in ds2.data_vars:
-        # Fill values must match dtype. Here we do float32 to match ELM surface expectations.
+        # Fill values must match dtype. Most float vars are written as float32,
+        # while closure-critical fractions stay float64 to preserve exact sums.
         if ds2[v].dtype.kind == "f":
-            enc[v] = {"dtype": "float32", "_FillValue": np.float32(-9.96921e36)}
+            if v in closure_critical:
+                enc[v] = {"dtype": "float64", "_FillValue": np.float64(-9.96921e36)}
+            else:
+                enc[v] = {"dtype": "float32", "_FillValue": np.float32(-9.96921e36)}
 
     ds2.to_netcdf(out_path, encoding=enc)
     return out_path
@@ -886,7 +892,7 @@ def export(
                 agg_policy=agg_policy,
             )
 
-            # Attach topounit parameters (and PCT_TOPUNIT) exactly once
+            # Attach topounit parameters (and TopounitFracArea) exactly once
             if attach_topounits and getattr(run_dom, "topounits", None) is not None and run_dom.topounits is not None:
                 sf.add_topounits_from_domain(run_dom)
 
@@ -1021,11 +1027,19 @@ def add_topounits_from_domain(
         id_col: str = "topounit_id",
         pct_col: str = "TopounitPctOfCell",
         dim_name: str = "topounit",
-        pct_var_name: str = "PCT_TOPUNIT",
+        pct_var_name: str = "TopounitFracArea",
     ) -> None:
         """
         Attach topounits + per-cell weights to the surface dataset.
 
+        This method:
+        - Adds a topounit dimension and coordinate to the dataset
+        - Creates a topounit fraction variable (default: TopounitFracArea) with values
+          normalized to decimal fractions (0.0-1.0, summing to 1.0 per cell)
+        - Expands topounit-indexed variables in SURFACE_VAR_SPECS (e.g., PCT_NAT_PFT)
+          to include the topounit dimension via uniform broadcast, where all topounits
+          in a cell inherit the parent cell's distribution
+
         Expects domain.topounits to exist and contain:
         - gid_col (links topounit -> cell gid)
         - id_col  (unique id per topounit across the whole run)
@@ -1091,8 +1105,8 @@ def add_topounits_from_domain(
             s = float(np.nansum(vals))
             if not np.isfinite(s) or s <= 0:
                 raise ValueError(f"Topounit pct weights for gid={gid} are invalid (sum={s}).")
-            # normalize to 100 just in case
-            vals = 100.0 * (vals / s)
+            # normalize to 1.0 (decimal fraction) just in case
+            vals = 1.0 * (vals / s)
 
             for tid, v in zip(grp[id_col].astype(str).tolist(), vals):
                 k = id_to_k[tid]
@@ -1108,7 +1122,31 @@ def add_topounits_from_domain(
                 raise ValueError(f"Existing {dim_name} coord does not match topounit ids from domain.")
 
         ds[pct_var_name] = xr.DataArray(pct, dims=(dim_name, lat_dim, lon_dim))
-        ds[pct_var_name].attrs.update({"long_name": "percent of gridcell in each topounit", "units": "percent"})
+        ds[pct_var_name].attrs.update({"long_name": "fraction of gridcell area in each topounit", "units": "unitless"})
+
+        # Expand topounit-indexed variables that exist in ds but currently lack the
+        # topounit dimension.  All topounits in a grid cell inherit the parent cell's
+        # distribution uniformly (per-topounit differentiation is a future extension).
+        top_coord = ds.coords[dim_name]
+        n_top = len(top_ids)
+        for _var_name, _spec in SURFACE_VAR_SPECS.items():
+            if _var_name not in ds:
+                continue
+            _spec_dims = [d.strip() for d in _spec.get("dims", "").split(",")]
+            if dim_name not in _spec_dims:
+                continue
+            _da = ds[_var_name]
+            if dim_name in _da.dims:
+                continue  # already has the topounit dim
+            # Repeat identical values across all topounits
+            _expanded = xr.concat([_da] * n_top,
+                                   dim=xr.DataArray(top_coord.values, dims=[dim_name], name=dim_name))
+            _expanded[dim_name] = top_coord
+            # Reorder dims to match spec order (topounit first, spatial last)
+            _existing = set(_expanded.dims)
+            _ordered = [d for d in _spec_dims if d in _existing]
+            _extra = [d for d in _expanded.dims if d not in _ordered]
+            ds[_var_name] = _expanded.transpose(*_ordered, *_extra)
 
         self.ds = ds