Overhaul the GeoClaw meteorological forcing documentation

doc/geoclaw.rst

@@ -69,7 +69,9 @@ More will eventually appear in the :ref:`apps`.
    fgmax
    nearshore_interp
    tsunamidata
+   met_forcing
    surgedata
+   storm_module
    netcdf
    netcdf_utils_module
    marching_front

doc/met_forcing.rst

@@ -0,0 +1,102 @@
+
+.. _met_forcing:
+
+**********************************************
+Meteorological (Storm) Forcing
+**********************************************
+
+GeoClaw can force the shallow water equations with wind and atmospheric
+pressure fields.  This *meteorological forcing* ("met forcing") drives storm
+surge, and also supports non-tropical applications (extratropical cyclones,
+meteo-tsunamis, and other atmospheric forcing).  A storm is one important
+subtype of met forcing rather than the whole story.
+
+This page gives an overview and points to the quick start, the ``setrun``
+reference, and the Python API.  For a hands-on walkthrough see
+:ref:`quick_surge`; for data sources see :ref:`surgedata`.
+
+Forcing families
+================
+
+Met forcing comes in two families, selected in ``setrun.py`` through
+``rundata.surge_data`` (see :ref:`setrun_surge`):
+
+**Parametric forcing**
+   Wind and pressure fields are generated from a compact, evolving description
+   of a storm (its track plus intensity parameters) using an analytic model
+   such as Holland 1980.  A parametric storm has an analytic center/eye, so it
+   can drive distance-to-center AMR refinement.  Supported models include
+   ``holland80``, ``holland2008``, ``holland2010``, ``cle``, ``slosh``,
+   ``rankine``, ``modified_rankine``, ``demaria``, and ``willoughby``.
+
+**Gridded forcing**
+   Wind and pressure are read from external field datasets and interpolated to
+   the grid, rather than assuming a profile.  Two formats are supported:
+   OWI/ASCII (NWS12) and NetCDF (NWS13).  Gridded forcing has no analytic
+   center, so refinement is driven by the wind-speed thresholds only.  See
+   :ref:`netcdf_input` for the NetCDF conventions, automatic variable
+   discovery, and unit/time handling.
+
+Selecting a forcing
+===================
+
+The explicit, preferred way to select forcing is a **family** plus a
+**subtype**::
+
+    rundata.surge_data.storm_family  = "parametric"   # or "gridded" / "none"
+    rundata.surge_data.storm_subtype = "holland80"    # model, or "gridded"
+    rundata.surge_data.storm_file    = "path/to/storm"
+
+See :ref:`setrun_surge` for the full attribute reference.
+
+Python object model
+===================
+
+The Python side is organized around met forcing rather than only storms
+(:ref:`storm_module` is the API reference):
+
+- :class:`~clawpack.geoclaw.surge.track.Track` -- a generic evolving feature
+  with a center over time, and
+  :class:`~clawpack.geoclaw.surge.track.StormTrack`, which adds storm metadata
+  (max wind speed, radius of maximum winds, central pressure, ...).
+- :class:`~clawpack.geoclaw.surge.parametric.ParametricMetForcing` -- forcing
+  from a parameterized model referencing a track.
+- :class:`~clawpack.geoclaw.surge.gridded.GriddedMetForcing` -- forcing from
+  external field datasets (OWI/ASCII, NetCDF).
+- :class:`~clawpack.geoclaw.surge.storm.Storm` -- a backwards-compatible
+  wrapper retaining the historical ``read``/``write`` interface; new code can
+  use the objects above directly.
+
+Track readers ingest many formats (ATCF, HURDAT, IBTrACS, JMA, TCVITALS);
+GeoClaw writes only the two forcing files it consumes: a parametric storm file
+(``write_geoclaw``) and a gridded descriptor (``write_data``).
+
+.. _met_forcing_migration:
+
+What's new / migrating
+======================
+
+Existing ``setrun.py`` scripts continue to work.  The changes are additive:
+
+- **Explicit family/subtype selection.** ``storm_family`` + ``storm_subtype``
+  replace the overloaded integer ``storm_specification_type`` at the API level.
+  The legacy selector -- either a model-name string (``'holland80'``,
+  ``'data'``) or the signed integer code -- is still fully supported and maps
+  to the same forcing.
+- **Gridded NetCDF forcing.** Full gridded wind/pressure forcing from
+  CF-compliant NetCDF files, with automatic variable/coordinate discovery and
+  unit handling (:ref:`netcdf_input`).
+- **Temporal ramps.** ``t_ramp_on`` / ``t_ramp_off`` ramp the forcing on and
+  off over a specified number of seconds.
+- **Object model.** The ``Track`` / ``StormTrack`` / ``ParametricMetForcing`` /
+  ``GriddedMetForcing`` classes above, with ``Storm`` retained as a
+  compatibility wrapper.
+
+.. note::
+
+   The Fortran modules were renamed to the met-forcing vocabulary
+   (``storm_module`` → ``met_forcing_module``,
+   ``model_storm_module`` → ``parametric_met_forcing_module``,
+   ``data_storm_module`` → ``gridded_met_forcing_module``).  This only affects
+   custom Fortran source that ``use``\ s those modules; ``setrun.py`` workflows
+   are unaffected.

doc/quick_surge.rst

@@ -8,48 +8,110 @@ Quick start guide for storm surge modeling
 
 See also this `youtube video <https://www.youtube.com/watch?v=YurKRmYgGfk&t=10s>`__
 and the related materials from the `2020 GeoClaw Developers Workshop
-<http://www.clawpack.org/geoclawdev-2020/>`__.
+<http://www.clawpack.org/geoclawdev-2020/>`__, the :ref:`met_forcing` overview,
+and the :ref:`setrun_surge` reference.
 
-To get started with a storm surge computation it is best to refer to a previous
-working example.  For example, you might start with 
-`$CLAW/geoclaw/examples/storm-surge/ike`.  There are also a number of additional 
-examples in the `$CLAW/geoclaw/examples/storm-surge` directory as well as some
-in the `$CLAW/apps/surge-examples` directory (this is actually a repository of 
-examples that is actively updated).  The primary input that one needs to provide
-for a new example usually involves two data source
+The quickest way to get started is to run a working example and then adapt it.
+A good starting point is ``$CLAW/geoclaw/examples/storm-surge/ike`` (Hurricane
+Ike, parametric Holland 1980 forcing).  There are further examples in
+``$CLAW/geoclaw/examples/storm-surge`` and in the actively maintained
+``$CLAW/apps/surge-examples`` repository.
 
- - Topography data:  Data that specifies the topography and bathymetry of the
-   region around the area of interest.  For storm surge computations it is 
-   generally good practice to include entire oceanic basins so that you can
-   ensure that flow into and out of the basin is resolved by the computation
-   and is sufficiently distant from the computational domain's boundaries.
- - Storm data:  Of course we need to specify the particular storm that you
-   are interested in.  There are a number of ways to specify a storm which
-   are described in :ref:`setrun_surge`.  Sources for parameterized storms
-   can also be found in :ref:`surgedata` and a description of how to include
-   them in :ref:`_surge_module`.
+Running the Ike example
+=======================
 
-.. warning:: This is a work in progress and only partially has been filled out.
-   If you are interested in the rest of the steps or wish to contribute your
-   own workflow please let us know!
+From ``$CLAW/geoclaw/examples/storm-surge/ike``::
 
-Here we will concentrate on changing the Hurricane Ike example into one for
-Hurricane Katrina.
+    make all
 
-1. First copy the files located in the Hurricane Ike directorty located at
-   `$CLAW/geoclaw/examples/storm-surge/ike`.
+This downloads the topography, writes the storm forcing file, compiles
+``xgeoclaw``, runs the simulation, and produces plots in ``_plots``.  The
+individual steps are also available as separate targets::
 
-2. Next let's find some better topography for the New Orleans area...
+    make .exe       # compile xgeoclaw
+    make data       # write the *.data files from setrun.py
+    make .output    # run the simulation
+    make .plots     # generate plots with setplot.py
 
-3. Now let's find a storm specification for Hurricane Katrina.  In this 
-   example we will use the ATCF database.  For Katrina this ends up being
-   the file located `here <>`_.
+The two inputs a surge computation needs
+========================================
 
-4. We now need to modify the `setrun.py` to use our new storm format and
-   topography we now added...
+- **Topography / bathymetry** for the region of interest.  For surge it is good
+  practice to include entire oceanic basins so that flow into and out of the
+  basin is resolved and the domain boundaries are far from the area of interest.
+  See :ref:`topo` for supported formats.
+- **Storm forcing** describing the wind and pressure.  There are two families
+  (parametric and gridded); see :ref:`met_forcing` for the overview and
+  :ref:`surgedata` for data sources.
 
-5. Finally we need to also modify the plotting so that we have an
+Configuring the forcing in ``setrun.py``
+========================================
 
-6. Gauges...
+Storm forcing is configured through ``rundata.surge_data`` (full reference in
+:ref:`setrun_surge`).  The Ike example enables the wind and pressure source
+terms and selects a parametric Holland 1980 storm::
 
-7. Running the simulation...
+    data = rundata.surge_data
+    data.wind_forcing = True
+    data.pressure_forcing = True
+    data.drag_law = 1                      # Garratt wind drag
+
+    # Preferred, explicit forcing selection:
+    data.storm_family = "parametric"
+    data.storm_subtype = "holland80"
+    data.storm_file = "ike.storm"          # GeoClaw-format storm file
+
+    # AMR refinement on wind speed (m/s) and distance to the eye (m):
+    data.wind_refine = [20.0, 40.0, 60.0]
+    data.R_refine = [60.0e3, 40.0e3, 20.0e3]
+
+The legacy selector ``data.storm_specification_type = 'holland80'`` is
+equivalent and still supported.
+
+Building the storm file
+=======================
+
+The GeoClaw-format storm file (``ike.storm`` above) is produced from a track in
+one of the ingest formats (ATCF, HURDAT, IBTrACS, JMA, TCVITALS).  Using the
+Python API (see :ref:`storm_module`)::
+
+    from clawpack.geoclaw.surge.storm import Storm
+    import numpy as np
+
+    storm = Storm(path="track.dat", file_format="ATCF")
+    storm.time_offset = np.datetime64("2008-09-13T07:00")   # e.g. landfall
+    storm.write("ike.storm", file_format="geoclaw")
+
+Plotting and gauges
+===================
+
+The example ``setplot.py`` uses the ``clawpack.geoclaw.surge.plot`` helpers
+(imported as ``surgeplot``) to plot the surface elevation, wind speed, and
+pressure fields, overlay the storm track, and plot gauge time series.  Gauges
+record the wind and pressure aux fields when ``rundata.gaugedata.aux_out_fields``
+includes the wind/pressure aux indices.
+
+Using gridded (OWI / NetCDF) forcing
+====================================
+
+To force GeoClaw with gridded fields instead of a parametric model, set the
+family to ``"gridded"`` and provide the field files.  For NetCDF this is a
+CF-compliant file of wind (``u10``, ``v10``) and mean-sea-level pressure
+(``msl``); GeoClaw discovers the variables and coordinates automatically.  See
+:ref:`netcdf_input` for the full workflow and the
+``$CLAW/geoclaw/examples/storm-surge/isaac`` example, which drives both
+parametric and gridded forcing from the same ATCF track.
+
+Adapting to a different storm
+=============================
+
+To model a different event, copy the Ike example directory and then:
+
+1. Obtain topography covering the affected basin (see :ref:`topo` and the data
+   sources linked from :ref:`surgedata`).
+2. Obtain the storm track (e.g. an ATCF best-track file; see :ref:`surgedata`)
+   and build the storm file as shown above, or point ``storm_file`` at a
+   gridded descriptor for gridded forcing.
+3. Update ``setrun.py`` -- the domain extent, refinement regions, ``surge_data``
+   selection, gauges, and run time -- for the new event.
+4. Update ``setplot.py`` for the new domain and gauges.

doc/setrun_geoclaw.rst

@@ -188,27 +188,50 @@ Topography preprocessing attributes
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Each :class:`~clawpack.geoclaw.topotools.Topography` entry in ``topofiles``
-supports seven preprocessing attributes that are applied automatically when
+supports eight preprocessing attributes that are applied automatically when
 the file is loaded by :meth:`~clawpack.geoclaw.topotools.Topography.read`.
 Set them on the object before appending to ``topofiles``:
 
-+----------------+----------------+-----------+-----------------------------------------------+
-| Attribute      | Type           | Default   | Description                                   |
-+================+================+===========+===============================================+
-| ``crop_extent``| list or None   | ``None``  | ``[x1, x2, y1, y2]`` crop region.            |
-+----------------+----------------+-----------+-----------------------------------------------+
-| ``coarsen``    | int            | ``1``     | Stride-subsampling factor (1 = no coarsen).  |
-+----------------+----------------+-----------+-----------------------------------------------+
-| ``buffer``     | int            | ``0``     | Grid-point margin outside crop region.        |
-+----------------+----------------+-----------+-----------------------------------------------+
-| ``align``      | tuple or None  | ``None``  | ``(x, y)`` alignment for coarsened grids.    |
-+----------------+----------------+-----------+-----------------------------------------------+
-| ``x_shift``    | float          | ``0.0``   | Constant added to all x coordinates.          |
-+----------------+----------------+-----------+-----------------------------------------------+
-| ``z_shift``    | float          | ``0.0``   | Constant added to non-missing Z values.       |
-+----------------+----------------+-----------+-----------------------------------------------+
-| ``negate_z``   | bool           | ``False`` | If ``True``, flip the sign of all Z values.  |
-+----------------+----------------+-----------+-----------------------------------------------+
+.. list-table::
+   :header-rows: 1
+   :widths: 16 16 10 48
+
+   * - Attribute
+     - Type
+     - Default
+     - Description
+   * - ``crop_extent``
+     - list or None
+     - ``None``
+     - ``[x1, x2, y1, y2]`` crop region.
+   * - ``coarsen``
+     - int
+     - ``1``
+     - Stride-subsampling factor (1 = no coarsen).
+   * - ``buffer``
+     - int
+     - ``0``
+     - Grid-point margin outside crop region.
+   * - ``align``
+     - tuple or None
+     - ``None``
+     - ``(x, y)`` alignment for coarsened grids.
+   * - ``x_shift``
+     - float
+     - ``0.0``
+     - Constant added to all x coordinates.
+   * - ``y_shift``
+     - float
+     - ``0.0``
+     - Constant added to all y coordinates.
+   * - ``z_shift``
+     - float
+     - ``0.0``
+     - Constant added to non-missing Z values.
+   * - ``negate_z``
+     - bool
+     - ``False``
+     - If ``True``, flip the sign of all Z values.
 
 .. note::
 
@@ -497,20 +520,36 @@ Storm Specification Data
    distance to the storm's center.  This can also be set to a boolean which if
    `False` disables storm radial based refinement.
 
-.. attribute:: rundata.surge_data.storm_specification_type : int
+.. attribute:: rundata.surge_data.storm_family : string
 
-   Specifies the type of storm being used.  Positive options refer to a 
-   parameterized storm model where as negative integers refer to fully 
-   specified storms, for instance from HWRF, to be specified.
+   The forcing family, and the preferred (explicit) way to select forcing
+   together with ``storm_subtype``.  Valid values:
 
-   Valid options 
+    - ``"parametric"``: an analytic model with a storm center/track (see
+      ``storm_subtype`` for the model list).
+    - ``"gridded"``: file-backed wind/pressure fields (OWI/ASCII or NetCDF).
+    - ``"none"``: forcing off.
 
-    - `-1`: The input data is specified in the HWRF format.
-    - `0`: No storm specified
-    - `1`: Parameterized storm requested using the Holland 1980 modeled storm.
-    - `2`: Parameterized storm requested using the Holland 2010 modeled storm.
-    - `3`: Parameterized storm requested using the Chava, Lin, Emmanuel modeled 
-      storm.
+.. attribute:: rundata.surge_data.storm_subtype : string
+
+   The forcing subtype within the family.  For ``"parametric"`` this is a model
+   name -- one of ``"holland80"``, ``"holland2008"``, ``"holland2010"``,
+   ``"cle"``, ``"slosh"``, ``"rankine"``, ``"modified_rankine"``,
+   ``"demaria"``, or ``"willoughby"``.  For ``"gridded"`` use ``"gridded"``
+   (the concrete OWI/NetCDF format is carried by the storm descriptor file).
+
+.. attribute:: rundata.surge_data.storm_specification_type : int or string
+
+   Legacy forcing selector, retained for backwards compatibility and still
+   fully supported.  Used when ``storm_family``/``storm_subtype`` are not set;
+   it may be a model-name string (e.g. ``'holland80'``, ``'data'``) or the
+   signed integer code below, and resolves to the same forcing.
+
+   Integer codes: ``0`` no storm; positive values select a parametric model
+   (``1`` Holland 1980, ``2`` Holland 2010, ``3`` Chavas--Lin--Emanuel, ``4``
+   SLOSH, ``5`` Rankine, ``6`` modified Rankine, ``7`` DeMaria, ``8`` Holland
+   2008, ``9`` Willoughby); a negative value (``-1``) selects gridded/data
+   forcing.
 
 .. attribute:: rundata.surge_data.storm_file : string
 
@@ -528,3 +567,23 @@ Storm Specification Data
    scaling).  This is primarily useful for sensitivity studies and for
    creating synthetic storms whose tracks are spatially identical to a
    historical event but travel at a different speed.
+
+.. attribute:: rundata.surge_data.t_ramp_on : float
+
+   Number of seconds over which the wind/pressure forcing ramps on after the
+   simulation start time.  ``0.0`` (the default) disables the onset ramp.
+   Useful for avoiding an impulsive start when the forcing is already
+   significant at ``t0``.
+
+.. attribute:: rundata.surge_data.t_ramp_off : float
+
+   Number of seconds over which the wind/pressure forcing ramps off before the
+   final time.  ``0.0`` (the default) disables the cutoff ramp.
+
+.. attribute:: rundata.surge_data.rotation_override : int or string
+
+   Overrides the hemisphere-based sense of the storm's rotation.  The default
+   (``0`` or ``"normal"``) chooses the rotation from the hemisphere of the
+   storm center; ``1`` / ``"N"`` forces northern (counter-clockwise) and
+   ``2`` / ``"S"`` forces southern (clockwise).  Primarily useful for
+   idealized or cross-hemisphere synthetic storms.