Add --output-file and --output-frequency options

CHANGELOG.md

@@ -32,6 +32,7 @@ SPDX-License-Identifier: CC-BY-4.0
 * Grids built from `--elmer-mesh` without a DEM now place each column at the latitude and longitude of its mesh cell center, as they already did with `--netcdf-dem-mesh`. Every column previously sat at 75 N 40 W regardless of the mesh. Results of such runs change through the solar geometry, unless the surface energy balance is taken from ICON-Land (`--couple-to-icon-land`). The `lon` written to the output now ranges from -180 to 180 degrees east instead of being 320. https://github.com/EBFMorg/EBFM/pull/193
 * Grids built from `--elmer-mesh` alone (no `--netcdf-dem-mesh`) now take each column's elevation from the mesh: the mean of the z coordinates of its cell's vertices, as `reader.py` writes them from a DEM. The elevation was previously random between 0 and 100 m. A mesh without elevation data, such as `examples/greenland_mesh_v0/MESH`, now gives 0 m everywhere. https://github.com/EBFMorg/EBFM/pull/193
 * Grids built from `--elmer-mesh` alone (no `--netcdf-dem-mesh`) now take each column's surface slope and orientation from the plane through the vertices of its mesh cell, instead of treating every column as flat. This changes the shortwave radiation on sloping columns. A mesh without elevation data still gives flat columns. When coupled to Elmer/Ice (`--couple-to-elmer-ice`), the columns stay flat and a warning is logged, since the slopes do not follow the surface elevation received from Elmer/Ice yet. https://github.com/EBFMorg/EBFM/pull/193
+* Add options to choose the path (`--output-file`, default `Output/model_output.nc`) and the frequency (`--output-frequency`, an ISO 8601 duration, default `P1D`) of the NetCDF output file. https://github.com/EBFMorg/EBFM/pull/211
 
 # v0.8.0
 

src/ebfm/core/INIT.py

@@ -11,12 +11,12 @@
 from netCDF4 import Dataset, num2date
 
 from pathlib import Path
-from datetime import datetime
+from datetime import datetime, timedelta
 
 from ebfm.reader import read_elmer_mesh, read_dem, read_dem_xios
 
 from ebfm.elmer.mesh import Mesh
-from .config import TimeConfig, GridConfig, ColumnDiscretizationConfig, iso8601
+from .config import TimeConfig, GridConfig, ColumnDiscretizationConfig, OutputConfig, iso8601
 from .grid import GridInputType, GridDict, ShadingMethod, domain_mean, number_of_columns, validate_grid
 from .restart import validate_all_variables_present, validate_variable_shape
 
@@ -41,12 +41,19 @@ def create_restart_file_name(time: datetime) -> str:
     return restart_file_prefix + iso8601(time) + ".nc"
 
 
-def init_config(time_config: TimeConfig, grid_config, restartdir: Path, initialize_from_restart_file: bool):
+def init_config(
+    time_config: TimeConfig,
+    grid_config,
+    output_config: OutputConfig,
+    restartdir: Path,
+    initialize_from_restart_file: bool,
+):
     """
     Set model parameters, specify grid parameters, I/O, and physics settings.
 
     @param[in] time_config Time configuration object.
     @param[in] grid_config Grid configuration object.
+    @param[in] output_config Output file configuration object.
     @param[in] restartdir Path to the directory containing restart files (if initializing from restart file).
     @param[in] initialize_from_restart_file Boolean flag indicating whether to initialize from a restart file.
 
@@ -84,7 +91,7 @@ def init_config(time_config: TimeConfig, grid_config, restartdir: Path, initiali
     io = {}
 
     io["homedir"] = Path(os.getcwd())  # Home directory
-    io["outdir"] = io["homedir"] / "Output"  # Output directory
+    io["outfile"] = output_config.output_file  # Output file
 
     write_restart_file: bool
 
@@ -106,10 +113,13 @@ def init_config(time_config: TimeConfig, grid_config, restartdir: Path, initiali
         write_restart_file = False
 
     io["writebootfile"] = write_restart_file
-    io["freqout"] = 8  # OUTPUT: frequency of storing output (every n-th time-step)
+    io["freqout"] = output_config.output_frequency // time_config.time_step  # OUTPUT: store every n-th time-step
+    assert (
+        output_config.output_frequency % time_config.time_step == timedelta()
+    ), "Output frequency must be a multiple of the time step."
 
     # Ensure output and reboot directories exist
-    os.makedirs(io["outdir"], exist_ok=True)
+    os.makedirs(io["outfile"].parent, exist_ok=True)
 
     # Return the initialized parameters
     return grid, io, phys, column_config

src/ebfm/core/LOOP_write_to_file.py

@@ -3,6 +3,7 @@
 # SPDX-License-Identifier: BSD-3-Clause
 
 import os
+from pathlib import Path
 
 import numpy as np
 from netCDF4 import Dataset, date2num
@@ -254,12 +255,10 @@ def main(OUTFILE, io, OUT, grid, t, time, column, gather: ColumnGather):
 
     # Initialize NetCDF file at the first time step
     if is_first_time_step(t) and gather.is_root:
-        if not os.path.exists(io["outdir"]):
-            os.makedirs(io["outdir"])
+        os.makedirs(Path(io["outfile"]).parent, exist_ok=True)
 
         # Create NetCDF file
-        nc_filepath = os.path.join(io["outdir"], "model_output.nc")
-        io["nc_file"] = Dataset(nc_filepath, "w", format="NETCDF4")
+        io["nc_file"] = Dataset(io["outfile"], "w", format="NETCDF4")
         io["nc_file"].createDimension("time", None)  # Unlimited time dimension
 
         if io["output_grid"]["is_unstructured"]:

src/ebfm/core/cli.py

@@ -14,7 +14,7 @@
 import argparse
 import os
 from argparse import ArgumentParser, Namespace
-from datetime import datetime
+from datetime import datetime, timedelta
 from enum import Enum
 from pathlib import Path
 
@@ -61,6 +61,8 @@ class CliDefaults(Enum):
     TIME_STEP_SIZE = "PT3H"  # = 3 hours
     LOG_LEVEL_CONSOLE = "INFO"
     COMPONENT_NAME = "ebfm"
+    OUTPUT_FILE = Path("Output") / "model_output.nc"
+    OUTPUT_FREQUENCY = "P1D"  # = daily
 
     @classmethod
     def default_time_step_size_in_hours(cls) -> float:
@@ -286,6 +288,36 @@ def parse_cli_args(args: list[str] | None = None) -> Namespace:
         help="Path to the folder containing NetCDF forcing files with meteorological data.",
     )
 
+    output_group = parser.add_argument_group("output file")
+
+    output_group.add_argument(
+        "--output-file",
+        type=Path,
+        default=CliDefaults.OUTPUT_FILE.value,
+        help="Path of the NetCDF output file. Relative paths are relative to the working directory. Missing "
+        "directories are created and an existing file is overwritten.",
+    )
+
+    def parse_output_frequency(value: str) -> timedelta:
+        try:
+            interval = isodate.parse_duration(value)
+        except isodate.ISO8601Error:
+            raise argparse.ArgumentTypeError(f"must be an ISO 8601 duration such as P1D, got '{value}'")
+        if not isinstance(interval, timedelta):
+            raise argparse.ArgumentTypeError(f"must not contain years or months, which vary in length, got '{value}'")
+        if interval <= timedelta(0):
+            raise argparse.ArgumentTypeError(f"must be positive, got '{value}'")
+        return interval
+
+    output_group.add_argument(
+        "--output-frequency",
+        type=parse_output_frequency,
+        default=CliDefaults.OUTPUT_FREQUENCY.value,
+        help="Interval between two output records in ISO 8601 format, e.g. P1D for daily or PT6H for six-hourly "
+        "output. Must be a multiple of --time-step. Each record holds the mean, sum, or a sample over this "
+        "interval, depending on the variable.",
+    )
+
     time_group = parser.add_argument_group("time configuration")
 
     time_group.add_argument(
@@ -613,6 +645,19 @@ def validate_forcing_parallel_compat(args: Namespace, comm) -> None:
         )
 
 
+def validate_output_frequency_time_step_compat(args: Namespace, time_config) -> None:
+    """Check that --output-frequency is a multiple of the time step.
+
+    @param[in] args        parsed CLI arguments
+    @param[in] time_config resolved TimeConfig
+    """
+    if args.output_frequency % time_config.time_step:
+        _cli_error(
+            f"--output-frequency {isodate.duration_isoformat(args.output_frequency)} is not a multiple of the time "
+            f"step {isodate.duration_isoformat(time_config.time_step)} (--time-step)."
+        )
+
+
 def validate_restart_parallel_compat(args: Namespace, comm) -> None:
     """Check that no restart files are read or written in a run on several MPI ranks.
 

src/ebfm/core/config/__init__.py

@@ -7,3 +7,4 @@
 from ebfm.core.config.grid import GridConfig, GridInputType  # noqa: F401
 from ebfm.core.config.forcing import ForcingConfig, ForcingType  # noqa: F401
 from ebfm.core.config.column import ColumnDiscretizationConfig  # noqa: F401
+from ebfm.core.config.output import OutputConfig  # noqa: F401

src/ebfm/core/config/output.py

@@ -0,0 +1,26 @@
+# SPDX-FileCopyrightText: 2026 EBFM Authors
+#
+# SPDX-License-Identifier: BSD-3-Clause
+"""Output configuration for EBFM.
+
+This file provides the configuration of the NetCDF output file for EBFM.
+"""
+
+from argparse import Namespace
+from datetime import timedelta
+from pathlib import Path
+
+
+class OutputConfig:
+    """Configuration of the NetCDF output file."""
+
+    def __init__(self, args: Namespace):
+        """
+        Initialize output configuration from command line arguments.
+
+        @param[in] args command line arguments
+        """
+        # Relative paths are relative to the working directory
+        self.output_file: Path = Path.cwd() / args.output_file
+        # Interval between two output records
+        self.output_frequency: timedelta = args.output_frequency

src/ebfm/main.py

@@ -18,7 +18,7 @@
 )
 from ebfm.core import LOOP_write_to_file, FINAL_create_restart_file
 from ebfm.core.comm.gather import ColumnGather
-from ebfm.core.config import CouplingConfig, GridConfig, TimeConfig, ForcingConfig
+from ebfm.core.config import CouplingConfig, GridConfig, TimeConfig, ForcingConfig, OutputConfig
 from ebfm.core.logger import Logger, setup_logging, log_levels_map, getLogger
 from ebfm.core.cli import (
     extract_active_coupling_features,
@@ -30,6 +30,7 @@
     validate_forcing_parallel_compat,
     validate_restart_parallel_compat,
     validate_shading_parallel_compat,
+    validate_output_frequency_time_step_compat,
 )
 
 import ebfm.coupling
@@ -222,11 +223,16 @@ def _main_impl():
     # Ensure shading routine is only used in runs on a single rank.
     validate_shading_parallel_compat(grid_config, ebfm_comm)
 
+    # Ensure every output record covers a whole number of time steps.
+    validate_output_frequency_time_step_compat(args, time_config)
+
     logger.debug("Successfully completed consistency checks.")
 
     # Model setup & initialization
     time = time_config.to_dict()
-    grid, io, phys, column = INIT.init_config(time_config, grid_config, args.restart_dir, args.restart_init)
+    grid, io, phys, column = INIT.init_config(
+        time_config, grid_config, OutputConfig(args), args.restart_dir, args.restart_init
+    )
 
     set_random_seed(logger, args.random_seed, ebfm_comm)
 

tests/core/test_cli.py

@@ -8,6 +8,7 @@
 import io
 import os
 from contextlib import redirect_stderr, redirect_stdout
+from datetime import timedelta
 from pathlib import Path
 from unittest.mock import MagicMock, patch
 
@@ -18,6 +19,7 @@
     parse_cli_args,
     validate_forcing_parallel_compat,
     validate_grid_parallel_compat,
+    validate_output_frequency_time_step_compat,
     validate_restart_parallel_compat,
     validate_shading_parallel_compat,
     validate_shading_coupling_compat,
@@ -194,6 +196,51 @@ def test_netcdf_dem_mesh_default_none(self):
     def test_random_seed_default_none(self):
         self.assertIsNone(self.args.random_seed)
 
+    def test_output_defaults(self):
+        self.assertEqual(self.args.output_file, Path("Output") / "model_output.nc")
+        self.assertEqual(self.args.output_frequency, timedelta(days=1))
+
+
+class TestOutputOptions(unittest.TestCase):
+    """--output-file / --output-frequency parsing and checks."""
+
+    def test_custom_values_accepted(self):
+        args = _parse("--output-file", "/scratch/run1.nc", "--output-frequency", "PT6H")
+        self.assertEqual(args.output_file, Path("/scratch/run1.nc"))
+        self.assertEqual(args.output_frequency, timedelta(hours=6))
+
+    def test_frequency_other_than_positive_fixed_duration_rejected(self):
+        cases = {
+            "integer": "8",
+            "float": "1.5",
+            "word": "daily",
+            "zero": "PT0S",
+            "negative": "-P1D",
+            "months vary in length": "P1M",
+            "years vary in length": "P1Y",
+        }
+        for case, value in cases.items():
+            with self.subTest(case, value=value):
+                with self.assertRaises(SystemExit), redirect_stderr(io.StringIO()):
+                    _parse("--output-frequency", value)
+
+    def test_frequency_must_be_multiple_of_time_step(self):
+        cases = {
+            "equal to time step": ("PT3H", "PT3H", True),
+            "multiple of time step": ("P1D", "PT3H", True),
+            "not a multiple": ("PT4H", "PT3H", False),
+            "shorter than time step": ("PT1H", "PT3H", False),
+        }
+        for case, (frequency, time_step, accepted) in cases.items():
+            with self.subTest(case, frequency=frequency, time_step=time_step):
+                args = _parse("--output-frequency", frequency, "--time-step", time_step)
+                time_config = TimeConfig(args)
+                if accepted:
+                    validate_output_frequency_time_step_compat(args, time_config)
+                else:
+                    with self.assertRaises(SystemExit), redirect_stderr(io.StringIO()):
+                        validate_output_frequency_time_step_compat(args, time_config)
+
 
 class TestGpuOptions(unittest.TestCase):
     """--with-gpu / --gpu-vendor parsing and compatibility checks."""

tests/core/test_init.py

@@ -20,7 +20,7 @@
 from ebfm.core.FINAL_create_restart_file import main as write_restart_file
 from ebfm.core.cli import parse_cli_args
 from ebfm.core.comm import defaultComm
-from ebfm.core.config import ColumnDiscretizationConfig
+from ebfm.core.config import ColumnDiscretizationConfig, OutputConfig
 from ebfm.core.config.grid import GridConfig
 from ebfm.core.config.time import TimeConfig
 from ebfm.core.restart import PER_COLUMN_VARIABLES, PER_LAYER_VARIABLES, RESTART_VARIABLES
@@ -336,6 +336,56 @@ def test_the_two_groups_do_not_overlap(self):
         self.assertEqual(len(RESTART_VARIABLES), len(set(RESTART_VARIABLES)))
 
 
+class TestOutputConfig(unittest.TestCase):
+    """The output options reach the io dict that LOOP_write_to_file reads."""
+
+    def setUp(self):
+        self.tmpdir = tempfile.TemporaryDirectory()
+
+    def tearDown(self):
+        self.tmpdir.cleanup()
+
+    def _init_io(self, *extra_args):
+        """Run init_config with the output file in the temporary directory and return its io dict."""
+        outfile = Path(self.tmpdir.name) / "runs" / "a" / "run1.nc"
+        args = parse_cli_args(["--matlab-mesh", str(_MATLAB_MESH), "--output-file", str(outfile), *extra_args])
+        _, io, _, _ = init_config(TimeConfig(args), GridConfig(args), OutputConfig(args), None, False)
+        return io
+
+    def test_default_output_file_is_model_output_nc_in_output_of_working_directory(self):
+        output_config = OutputConfig(parse_cli_args(["--matlab-mesh", str(_MATLAB_MESH)]))
+
+        self.assertEqual(output_config.output_file, Path.cwd() / "Output" / "model_output.nc")
+
+    def test_relative_output_file_is_relative_to_working_directory(self):
+        output_config = OutputConfig(parse_cli_args(["--matlab-mesh", str(_MATLAB_MESH), "--output-file", "a.nc"]))
+
+        self.assertEqual(output_config.output_file, Path.cwd() / "a.nc")
+
+    def test_output_file_ends_up_in_io_and_its_directory_is_created(self):
+        io = self._init_io()
+
+        self.assertEqual(io["outfile"], Path(self.tmpdir.name) / "runs" / "a" / "run1.nc")
+        self.assertTrue(io["outfile"].parent.is_dir())
+
+    def test_default_output_frequency_is_daily(self):
+        io = self._init_io()
+
+        self.assertEqual(io["freqout"], 8)  # time steps of the default PT3H per day
+
+    def test_output_frequency_is_converted_to_time_steps(self):
+        cases = {
+            "daily, hourly steps": ("P1D", "PT1H", 24),
+            "six-hourly, three-hourly steps": ("PT6H", "PT3H", 2),
+            "every time step": ("PT3H", "PT3H", 1),
+            "weekly, daily steps": ("P1W", "P1D", 7),
+        }
+        for case, (frequency, time_step, expected) in cases.items():
+            with self.subTest(case):
+                io = self._init_io("--output-frequency", frequency, "--time-step", time_step)
+                self.assertEqual(io["freqout"], expected)
+
+
 class TestMatlabShadingLookupTable(unittest.TestCase):
     """The shading look-up table must only be pre-computed when shading is enabled.
 
@@ -354,7 +404,7 @@ def tearDown(self):
     def _init_grid(self, shading_flag):
         args = parse_cli_args(["--matlab-mesh", str(_MATLAB_MESH), shading_flag])
         grid_config = GridConfig(args)
-        grid, io, _, _ = init_config(TimeConfig(args), grid_config, Path(self.tmpdir.name), False)
+        grid, io, _, _ = init_config(TimeConfig(args), grid_config, OutputConfig(args), Path(self.tmpdir.name), False)
         return init_grid(grid, io, grid_config, defaultComm)
 
     def test_shading_disabled_skips_the_lookup_table(self):
@@ -392,7 +442,7 @@ def tearDown(self):
     def _init_grid(self, num_ranks, rank):
         args = parse_cli_args(["--matlab-mesh", str(_MATLAB_MESH), "--no-shading"])
         grid_config = GridConfig(args, num_ranks, rank)
-        grid, io, _, _ = init_config(TimeConfig(args), grid_config, Path(self.tmpdir.name), False)
+        grid, io, _, _ = init_config(TimeConfig(args), grid_config, OutputConfig(args), Path(self.tmpdir.name), False)
         return init_grid(grid, io, grid_config, defaultComm)
 
     def test_blocks_of_all_ranks_make_up_the_serial_grid(self):
@@ -449,7 +499,7 @@ def _init_grid(self, extra_args=()):
             ["--elmer-mesh", self.tmpdir.name, "--elmer-mesh-crs-epsg", str(self._CRS_EPSG), *extra_args]
         )
         grid_config = GridConfig(args)
-        grid, io, _, _ = init_config(TimeConfig(args), grid_config, Path(self.tmpdir.name), False)
+        grid, io, _, _ = init_config(TimeConfig(args), grid_config, OutputConfig(args), Path(self.tmpdir.name), False)
         with mock.patch.object(INIT, "read_elmer_mesh", return_value=self.mesh):
             return init_grid(grid, io, grid_config, defaultComm)