From b6b63083cb082de8cb97479db635c90ff3617713 Mon Sep 17 00:00:00 2001
From: "Barnhart, Katherine (Katy) Ruth" <krbarnhart@usgs.gov>
Date: Fri, 17 May 2024 15:52:47 -0600
Subject: [PATCH 01/41] initial commit of xarray backends and example use file

---
 .../use_xarray_backends.py                    |  30 ++
 src/python/geoclaw/xarray_backends.py         | 466 ++++++++++++++++++
 2 files changed, 496 insertions(+)
 create mode 100644 examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
 create mode 100644 src/python/geoclaw/xarray_backends.py

diff --git a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
new file mode 100644
index 000000000..f0e72e714
--- /dev/null
+++ b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
@@ -0,0 +1,30 @@
+try:
+    import xarray as xr
+except:
+    'You must install xarray in order to write a use the xarray backends'
+    raise
+
+from clawpack.geoclaw.xarray_backends import FGOutBackend, FGMaxBackend
+
+# epsg code for lat-lon
+# Optionally, provide an epsg code to assign the associated coordinate system to the file.
+# default behavior assigns no coordinate system. Note that if no coordinate system is provided,
+# it will come in a GIS with coordinates associated with row and column number (not the x and y position
+# encoded in the netcdf).
+
+epsg_code = 4326
+
+# An example of a fgout grid.
+filename = '_output/fgout0001.b0001'
+# provide the .bxxx file if binary format is used or the
+# .qxxx file if ascii format is used.
+# the format, fg number, and frame number are inferred from the filename.
+
+ds = xr.open_dataset(filename, engine=FGOutBackend, backend_kwargs={'epsg':epsg_code})
+# ds is now an xarray object. It can be interacted with directly or written to netcdf using
+ds.to_netcdf('fgout0001_0001.nc')
+
+# An example of a fgmax grid.
+filename = "_output/fgmax0001.txt"
+ds = xr.open_dataset(filename, engine=FGMaxBackend, backend_kwargs={'epsg':epsg_code})
+ds.to_netcdf('fgmax0001.nc')
diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
new file mode 100644
index 000000000..3bbfee22e
--- /dev/null
+++ b/src/python/geoclaw/xarray_backends.py
@@ -0,0 +1,466 @@
+r"""
+xarray backends module: $CLAW/geoclaw/src/python/geoclaw/xarray_backends.py
+
+Xarray backends for GeoClaw fixed grids and fgmax grids.
+
+These only work for point_style = 2 (uniform regular)
+
+The expectation is that you run commands that use these backends from the same
+directory you run "make output".
+
+Includes:
+
+- class FGMaxBackend: Xarray backend for fgmax grids.
+- class FGOutBackend: Xarray backend for fgout grids.
+
+Usage:
+
+.. code-block:: python
+    import xarray as xr
+    from clawpack.geoclaw.xarray_backends import FGOutBackend, FGMaxBackend
+
+    # An example of a fgout grid.
+
+    filename = '_output/fgout0001.b0001
+    # provide the .bxxx file if binary format is used or the
+    # .qxxx file if ascii format is used.
+    # the format, fg number, and frame number are inferred from the filename.
+
+    ds = xr.open_dataset(filename, engine=FGOutBackend, backend_kwargs={'epsg':epsg_code})
+    # ds is now an xarray object. It can be interacted with directly or written to netcdf using
+    ds.write_netcdf('filename.nc')
+
+    # Optionally, provide an epsg code to assign the associated coordinate system to the file.
+    # default behavior assigns no coordinate system.
+
+    # An example of a fgmax grid.
+    filename = "_output/fgmax0001.txt"
+    ds = xr.open_dataset(filename, engine=FGMaxBackend, backend_kwargs={'epsg':epsg_code})
+
+
+Dimensions:
+
+Files opened with FGOutBackend will have dimensions time, y, x.
+Files opened with FGMaxBackend will have dimensions y, x.
+
+Variable naming:
+
+For fixed grid geoclaw files, the dataset will have the following variables:
+- h
+- hu
+- hv
+- eta
+
+Fixed grid dclaw files will have
+- h
+- hu
+- hv
+- hm
+- pb
+- hchi
+- delta_a
+- eta
+
+Depending on the number of variables specified in the setrun.py fgmax files will
+have a portion of the following variables:
+
+If rundata.fgmax_data.num_fgmax_val == 1
+
+- arrival_time, Wave arrival time (based on eta>sea_level + fg.arrival_tol)
+- h_max, Maximum water depth
+- eta_max, Maximum water surface elevation
+- h_max_time, Time of maximum water depth
+
+If rundata.fgmax_data.num_fgmax_val == 2:
+
+- s_max, Maximum velocity
+- s_max_time, Time of maximum velocity
+
+If rundata.fgmax_data.num_fgmax_val == 5:
+
+- hs_max, Maximum momentum
+- hs_max_time, Time of maximum momentum
+- hss_max, Maximum momentum flux
+- hss_max_time, Time of maximum momentum flux
+- h_min, Minimum depth
+- h_min_time, Time of minimum depth
+
+
+See the following links for additional information about xarray Backends.
+
+- https://docs.xarray.dev/en/stable/generated/xarray.backends.BackendEntrypoint.html#xarray.backends.BackendEntrypoint
+- import https://docs.xarray.dev/en/stable/generated/xarray.open_dataset.html
+"""
+
+import os
+
+import numpy as np
+import rioxarray  # activate the rio accessor
+import xarray as xr
+from clawpack.geoclaw import fgmax_tools, fgout_tools
+from xarray.backends import BackendEntrypoint
+
+_qelements_dclaw = ["h", "hu", "hv", "hm", "pb", "hchi", "delta_a", "eta"]
+_qelements_geoclaw = ["h", "hu", "hv", "eta"]
+
+_qunits = {
+    "h": "meters",
+    "hu": "meters squared per second",
+    "hv": "meters squared per second",
+    "hm": "meters",
+    "pb": "newton per meter squared",
+    "hchi": "meters",
+    "delta_a": "meters",
+    "eta": "meters",
+}
+
+time_unit = "seconds"
+reference_time = "model start"
+space_unit = "meters"
+nodata = np.nan
+
+
+def _prepare_var(data, mask, dims, units, nodata, long_name):
+    "Reorient array into the format xarray expects and generate the mapping."
+    data[mask] = nodata
+    data = data.T
+    data = np.flipud(data)
+    mapping = (
+        dims,
+        data,
+        {"units": units, "_FillValue": nodata, "long_name": long_name},
+    )
+    return mapping
+
+
+class FGOutBackend(BackendEntrypoint):
+    "Xarray Backend for Clawpack fixed grid format."
+
+    def open_dataset(
+        self,
+        filename,  # path to fgout file.
+        epsg=None,  # epsg code
+        drop_variables=None,  # name of any elements of q to drop.
+    ):
+
+        if drop_variables is None:
+            drop_variables = []
+
+        full_path = os.path.abspath(filename)
+        filename = os.path.basename(full_path)
+        outdir = os.path.basename(os.path.dirname(full_path))
+
+        # filename has the format fgoutXXXX.qYYYY (ascii)
+        # or fgoutXXXX.bYYYY (binary)
+        # where XXXX is the fixed grid number and YYYY is the frame
+        # number.
+        type_code = filename.split(".")[-1][0]
+        fgno = int(filename.split(".")[0][-4:])
+        frameno = int(filename.split(".")[-1][-4:])
+        if type_code == "q":  # TODO, is this correct?
+            output_format = "ascii"
+        elif type_code == "b":
+            output_format = "binary32"  # format of fgout grid output
+        else:
+            raise ValueError("Invalid FGout output format. Must be ascii or binary.")
+
+        fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, output_format)
+
+        if fgout_grid.point_style != 2:
+            raise ValueError("FGOutBackend only works with fg.point_style=2")
+
+        fgout_grid.read_fgout_grids_data()
+        fgout = fgout_grid.read_frame(frameno)
+
+        time = fgout.t
+        # both come in ascending. flip to give expected order.
+        x = fgout.x
+        y = np.flipud(fgout.y)
+        nj = len(x)
+        ni = len(y)
+
+        # mask based on dry tolerance
+        mask = fgout.q[0, :, :] < 0.001
+
+        # determine if geoclaw or dclaw
+        nvars = fgout.q.shape[0]
+        if nvars == 8:
+            _qelements = _qelements_dclaw
+        else:
+            _qelements = _qelements_geoclaw
+
+        # create data_vars dictionary
+        data_vars = {}
+        for i in range(nvars):
+            Q = fgout.q[i, :, :]
+            # mask all but eta based on h presence.
+            if i < 7:
+                Q[mask] = nodata
+
+            # to keep xarray happy, need to transpose and flip ud.
+            Q = Q.T
+            Q = np.flipud(Q)
+            Q = Q.reshape((1, ni, nj))  # reshape to add a time dimension.
+
+            # construct variable
+            varname = _qelements[i]
+
+            data_array_attrs = {"units": _qunits[varname], "_FillValue": nodata}
+
+            if varname not in drop_variables:
+                data_vars[varname] = (
+                    [
+                        "time",
+                        "y",
+                        "x",
+                    ],
+                    Q,
+                    data_array_attrs,
+                )
+
+        ds_attrs = {"description": "Clawpack model output"}
+
+        ds = xr.Dataset(
+            data_vars=data_vars,
+            coords=dict(
+                x=(["x"], x, {"units": space_unit}),
+                y=(["y"], y, {"units": space_unit}),
+                time=("time", [time], {"units": "seconds"}),
+                reference_time=reference_time,
+            ),
+            attrs=ds_attrs,
+        )
+
+        if epsg is not None:
+            ds.rio.write_crs(
+                epsg,
+                inplace=True,
+            ).rio.set_spatial_dims(
+                x_dim="x",
+                y_dim="y",
+                inplace=True,
+            ).rio.write_coordinate_system(inplace=True)
+            # https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html#Spatial-dimensions
+            # https://gis.stackexchange.com/questions/470207/how-to-write-crs-info-to-netcdf-in-a-way-qgis-can-read-python-xarray
+
+        return ds
+
+    open_dataset_parameters = ["filename", "drop_variables"]
+
+    description = "Use Clawpack fixed grid output files in Xarray"
+    url = "https://www.clawpack.org/fgout.html"
+
+
+class FGMaxBackend(BackendEntrypoint):
+    "Xarray Backend for Clawpack fgmax grid format."
+
+    def open_dataset(
+        self,
+        filename,
+        epsg=None,
+        drop_variables=None,
+    ):
+
+        if drop_variables is None:
+            drop_variables = []
+
+        # expectation is that you are in a run directory with a file 'fgmax_grids.data' in it.
+        fgno = int(os.path.basename(filename).split(".")[0][-4:])
+
+        fg = fgmax_tools.FGmaxGrid()
+        fg.read_fgmax_grids_data(fgno=fgno)
+        if fg.point_style != 2:
+            raise ValueError("FGMaxBackend only works with fg.point_style=2")
+
+        fg.read_output()
+
+        # Construct the x and y coordinates
+        # Both come in ascending, therefore flip y so that it is ordered as expected.
+        x = fg.x
+        y = np.flipud(fg.y)
+
+        # Construct the data_vars array. To organize the
+        # data in the way expected by netcdf standards, need
+        # to both transpose and flipud the array.
+
+        data_vars = {}
+
+        data_vars["arrival_time"] = _prepare_var(
+            fg.arrival_time.data,
+            fg.arrival_time.mask,
+            [
+                "y",
+                "x",
+            ],
+            "seconds",
+            nodata,
+            "Wave arrival time",
+        )
+
+        data_vars["h_max"] = _prepare_var(
+            fg.h.data,
+            fg.h.mask,
+            [
+                "y",
+                "x",
+            ],
+            "meters",
+            nodata,
+            "Maximum water depth",
+        )
+
+        data_vars["eta_max"] = _prepare_var(
+            fg.h.data + fg.B.data,
+            fg.h.mask,
+            [
+                "y",
+                "x",
+            ],
+            "meters",
+            nodata,
+            "Maximum water surface elevation",
+        )
+
+        data_vars["h_max_time"] = _prepare_var(
+            fg.h_time.data,
+            fg.h_time.mask,
+            [
+                "y",
+                "x",
+            ],
+            "seconds",
+            nodata,
+            "Time of maximum water depth",
+        )
+
+        if hasattr(fg, "s"):
+            if hasattr(fg.s, "data"):
+                data_vars["s_max"] = _prepare_var(
+                    fg.s.data,
+                    fg.s.mask,
+                    [
+                        "y",
+                        "x",
+                    ],
+                    "meters per second",
+                    nodata,
+                    "Maximum velocity",
+                )
+                data_vars["s_max_time"] = _prepare_var(
+                    fg.s_time.data,
+                    fg.s_time.mask,
+                    [
+                        "y",
+                        "x",
+                    ],
+                    "seconds",
+                    nodata,
+                    "Time of maximum velocity",
+                )
+        if hasattr(fg, "hs"):
+            if hasattr(fg.hs, "data"):
+                data_vars["hs_max"] = _prepare_var(
+                    fg.hs.data,
+                    fg.hs.mask,
+                    [
+                        "y",
+                        "x",
+                    ],
+                    "meters squared per second",
+                    nodata,
+                    "Maximum momentum",
+                )
+                data_vars["hs_max_time"] = _prepare_var(
+                    fg.hs_time.data,
+                    fg.hs_time.mask,
+                    [
+                        "y",
+                        "x",
+                    ],
+                    "seconds",
+                    nodata,
+                    "Time of maximum momentum",
+                )
+
+                data_vars["hss_max"] = _prepare_var(
+                    fg.hss.data,
+                    fg.hss.mask,
+                    [
+                        "y",
+                        "x",
+                    ],
+                    "meters cubed per second squared",
+                    nodata,
+                    "Maximum momentum flux",
+                )
+                data_vars["hss_max_time"] = _prepare_var(
+                    fg.hss_time.data,
+                    fg.hss_time.mask,
+                    [
+                        "y",
+                        "x",
+                    ],
+                    "seconds",
+                    nodata,
+                    "Time of maximum momentum flux",
+                )
+
+                data_vars["h_min"] = _prepare_var(
+                    fg.hmin.data,
+                    fg.hmin.mask,
+                    [
+                        "y",
+                        "x",
+                    ],
+                    "meters",
+                    nodata,
+                    "Minimum depth",
+                )
+                data_vars["h_min_time"] = _prepare_var(
+                    fg.hmin_time.data,
+                    fg.hmin_time.mask,
+                    [
+                        "y",
+                        "x",
+                    ],
+                    "seconds",
+                    nodata,
+                    "Time of minimum depth",
+                )
+
+        # drop requested variables
+        for var in drop_variables:
+            if var in data_vars.keys():
+                del data_vars[var]
+
+        # Construct the values from
+        ds_attrs = {"description": "D-Claw model output"}
+
+        ds = xr.Dataset(
+            data_vars=data_vars,
+            coords=dict(
+                x=(["x"], x, {"units": "meters"}),
+                y=(["y"], y, {"units": "meters"}),
+            ),
+            attrs=ds_attrs,
+        )
+
+        if epsg is not None:
+
+            ds.rio.write_crs(
+                epsg,
+                inplace=True,
+            ).rio.set_spatial_dims(
+                x_dim="x",
+                y_dim="y",
+                inplace=True,
+            ).rio.write_coordinate_system(inplace=True)
+            # https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html#Spatial-dimensions
+            # https://gis.stackexchange.com/questions/470207/how-to-write-crs-info-to-netcdf-in-a-way-qgis-can-read-python-xarray
+
+        return ds
+
+    open_dataset_parameters = ["filename", "drop_variables"]
+
+    description = "Use Clawpack fix grid monitoring files in Xarray"
+    url = "https://www.clawpack.org/fgmax.html#fgmax"

From ef4c7781d74c878bbfb806e3d7646f1496b95f94 Mon Sep 17 00:00:00 2001
From: Katy Barnhart <krbarnhart@usgs.gov>
Date: Fri, 17 May 2024 16:04:43 -0600
Subject: [PATCH 02/41] remove extra words

---
 examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
index f0e72e714..9a290d2d1 100644
--- a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
+++ b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
@@ -1,7 +1,7 @@
 try:
     import xarray as xr
 except:
-    'You must install xarray in order to write a use the xarray backends'
+    'You must install xarray in order to use the xarray backends'
     raise
 
 from clawpack.geoclaw.xarray_backends import FGOutBackend, FGMaxBackend

From 744e0e6382e802c5ac3111bbb82281af98cfe9e5 Mon Sep 17 00:00:00 2001
From: "Barnhart, Katherine (Katy) Ruth" <krbarnhart@usgs.gov>
Date: Sat, 18 May 2024 07:14:02 -0600
Subject: [PATCH 03/41] relax assumption of where command is issued (no longer
 needs to be in run directory).

---
 src/python/geoclaw/xarray_backends.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index 3bbfee22e..82e44b2b8 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -264,15 +264,18 @@ def open_dataset(
         if drop_variables is None:
             drop_variables = []
 
-        # expectation is that you are in a run directory with a file 'fgmax_grids.data' in it.
+        # expectation is for standard clawpack organization,
+        # e.g., output and 'fgmax_grids.data' in _output
         fgno = int(os.path.basename(filename).split(".")[0][-4:])
+        outdir = os.path.dirname(filename)
+        data_file = os.path.join(outdir, 'fgmax_grids.data')
 
         fg = fgmax_tools.FGmaxGrid()
-        fg.read_fgmax_grids_data(fgno=fgno)
+        fg.read_fgmax_grids_data(fgno=fgno, data_file=data_file)
         if fg.point_style != 2:
             raise ValueError("FGMaxBackend only works with fg.point_style=2")
 
-        fg.read_output()
+        fg.read_output(outdir=outdir)
 
         # Construct the x and y coordinates
         # Both come in ascending, therefore flip y so that it is ordered as expected.

From 9925e7c607421d5cb8b84208de441efa6c3bd1bc Mon Sep 17 00:00:00 2001
From: Katy Barnhart <krbarnhart@usgs.gov>
Date: Wed, 5 Jun 2024 14:26:25 -0600
Subject: [PATCH 04/41] add more information about rioxarray

---
 .../chile2010_fgmax-fgout/use_xarray_backends.py       |  3 ++-
 src/python/geoclaw/xarray_backends.py                  | 10 +++++++++-
 2 files changed, 11 insertions(+), 2 deletions(-)

diff --git a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
index 9a290d2d1..5e4fc63d0 100644
--- a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
+++ b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
@@ -1,7 +1,8 @@
 try:
+    import rioxarray
     import xarray as xr
 except:
-    'You must install xarray in order to use the xarray backends'
+    'You must install xarray and rioxarray in order to use the xarray backends'
     raise
 
 from clawpack.geoclaw.xarray_backends import FGOutBackend, FGMaxBackend
diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index 82e44b2b8..d5eb1cf72 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -3,7 +3,14 @@
 
 Xarray backends for GeoClaw fixed grids and fgmax grids.
 
-These only work for point_style = 2 (uniform regular)
+These only work for point_style = 2 (uniform regular) and have a dependency on
+xarray and rioxarray. Xarray provides the core datastructure and rioxarray
+provides an interface to rasterio, used to assign geospatial projection
+information.
+
+- https://docs.xarray.dev/en/stable/index.html
+- https://corteva.github.io/rioxarray/stable/readme.html
+- https://rasterio.readthedocs.io/en/latest/index.html
 
 The expectation is that you run commands that use these backends from the same
 directory you run "make output".
@@ -16,6 +23,7 @@
 Usage:
 
 .. code-block:: python
+    import rioxarray  # activate the rio accessor
     import xarray as xr
     from clawpack.geoclaw.xarray_backends import FGOutBackend, FGMaxBackend
 

From 51fb840a2485250835416f2f8b82ff7d085c9c07 Mon Sep 17 00:00:00 2001
From: Katy Barnhart <krbarnhart@usgs.gov>
Date: Wed, 5 Jun 2024 14:26:41 -0600
Subject: [PATCH 05/41] better import error for rioxarray

---
 src/python/geoclaw/xarray_backends.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index d5eb1cf72..ae9afb7d2 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -103,8 +103,12 @@
 import os
 
 import numpy as np
-import rioxarray  # activate the rio accessor
-import xarray as xr
+try:
+    import rioxarray  # activate the rio accessor
+    import xarray as xr
+except ImportError:
+    raise ImportError('rioxarray and xarray are required to use the FGOutBackend and FGMaxBackend')
+
 from clawpack.geoclaw import fgmax_tools, fgout_tools
 from xarray.backends import BackendEntrypoint
 

From a8bc42c12e08d44f6d6a4ecd7dffcf83907f80da Mon Sep 17 00:00:00 2001
From: Katy Barnhart <krbarnhart@usgs.gov>
Date: Wed, 5 Jun 2024 14:26:55 -0600
Subject: [PATCH 06/41] add B and level

---
 src/python/geoclaw/xarray_backends.py | 29 ++++++++++++++++++++++++++-
 1 file changed, 28 insertions(+), 1 deletion(-)

diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index ae9afb7d2..9ab4bc6c5 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -78,6 +78,8 @@
 - h_max, Maximum water depth
 - eta_max, Maximum water surface elevation
 - h_max_time, Time of maximum water depth
+- B, Basal topography at the first time fgmax first monitored maximum amr level
+- level, Maximum amr level
 
 If rundata.fgmax_data.num_fgmax_val == 2:
 
@@ -134,7 +136,8 @@
 
 def _prepare_var(data, mask, dims, units, nodata, long_name):
     "Reorient array into the format xarray expects and generate the mapping."
-    data[mask] = nodata
+    if mask is not None:
+        data[mask] = nodata
     data = data.T
     data = np.flipud(data)
     mapping = (
@@ -348,6 +351,30 @@ def open_dataset(
             "Time of maximum water depth",
         )
 
+        data_vars["B"] = _prepare_var(
+            fg.B,
+            None,
+            [
+                "y",
+                "x",
+            ],
+            "meters",
+            nodata,
+            "Basal topography at the first time fgmax first monitored maximum amr level",
+        )
+
+        data_vars["level"] = _prepare_var(
+            fg.level,
+            None,
+            [
+                "y",
+                "x",
+            ],
+            "(no units)",
+            -1,
+            "Maximum amr level",
+        )
+
         if hasattr(fg, "s"):
             if hasattr(fg.s, "data"):
                 data_vars["s_max"] = _prepare_var(

From dac2a2b7685184a23a86e27f053f69df936540e2 Mon Sep 17 00:00:00 2001
From: "Barnhart, Katherine (Katy) Ruth" <krbarnhart@usgs.gov>
Date: Fri, 5 Jul 2024 15:35:18 -0600
Subject: [PATCH 07/41] add clip example

---
 .../use_xarray_backends.py                    |  6 ++++
 src/python/geoclaw/xarray_backends.py         | 34 +++++++++++++++++--
 2 files changed, 38 insertions(+), 2 deletions(-)

diff --git a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
index 5e4fc63d0..94b55b3b8 100644
--- a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
+++ b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
@@ -29,3 +29,9 @@
 filename = "_output/fgmax0001.txt"
 ds = xr.open_dataset(filename, engine=FGMaxBackend, backend_kwargs={'epsg':epsg_code})
 ds.to_netcdf('fgmax0001.nc')
+
+# To see the use of clipping, change the tfinal in setrun to something like 2*3600.0
+# the fgmax0001_clipped.nc will only be the area where the wave arrived within the considered time. 
+filename = "_output/fgmax0001.txt"
+ds = xr.open_dataset(filename, engine=FGMaxBackend, backend_kwargs={'epsg':epsg_code, "clip": True})
+ds.to_netcdf('fgmax0001_clipped.nc')
diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index 9ab4bc6c5..77e445638 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -99,7 +99,7 @@
 See the following links for additional information about xarray Backends.
 
 - https://docs.xarray.dev/en/stable/generated/xarray.backends.BackendEntrypoint.html#xarray.backends.BackendEntrypoint
-- import https://docs.xarray.dev/en/stable/generated/xarray.open_dataset.html
+- https://docs.xarray.dev/en/stable/generated/xarray.open_dataset.html
 """
 
 import os
@@ -274,6 +274,8 @@ def open_dataset(
         filename,
         epsg=None,
         drop_variables=None,
+        clip=False, # if True, clip the entire array to the extent where arrival_time is defined.
+
     ):
 
         if drop_variables is None:
@@ -470,6 +472,7 @@ def open_dataset(
                     "Time of minimum depth",
                 )
 
+
         # drop requested variables
         for var in drop_variables:
             if var in data_vars.keys():
@@ -500,9 +503,36 @@ def open_dataset(
             # https://corteva.github.io/rioxarray/stable/getting_started/crs_management.html#Spatial-dimensions
             # https://gis.stackexchange.com/questions/470207/how-to-write-crs-info-to-netcdf-in-a-way-qgis-can-read-python-xarray
 
+        # clip
+        if clip:
+            clip_data = fg.arrival_time.data>=0
+            clip_data = clip_data.T
+            clip_data = np.flipud(clip_data)
+            ds = _clip(ds, clip_data)
+
         return ds
 
-    open_dataset_parameters = ["filename", "drop_variables"]
+    open_dataset_parameters = ["filename", "drop_variables", "clip"]
 
     description = "Use Clawpack fix grid monitoring files in Xarray"
     url = "https://www.clawpack.org/fgmax.html#fgmax"
+
+
+def _clip(ds, clip_data):
+    # clip DS to the extent where clip_data is true.
+    # useful for when there are large areas where nothing happened.
+
+    nrow, ncol = clip_data.shape
+
+    valid_col = np.sum(clip_data, axis=0)>0
+    valid_row = np.sum(clip_data, axis=1)>0
+
+    sel_rows = np.nonzero(valid_row)[0]
+    sel_cols = np.nonzero(valid_col)[0]
+
+    first_row = sel_rows[0]
+    last_row = sel_rows[-1]
+    first_col = sel_cols[0]
+    last_col = sel_cols[-1]
+
+    ds = ds.isel(y=np.arange(first_row, last_row), x=np.arange(first_col, last_col))

From d83bd2351e892501cc98f1ebef50aaa9bcdc7245 Mon Sep 17 00:00:00 2001
From: "Barnhart, Katherine (Katy) Ruth" <krbarnhart@usgs.gov>
Date: Fri, 5 Jul 2024 15:41:01 -0600
Subject: [PATCH 08/41] return ds

---
 src/python/geoclaw/xarray_backends.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index 77e445638..8f06f6fa4 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -536,3 +536,4 @@ def _clip(ds, clip_data):
     last_col = sel_cols[-1]
 
     ds = ds.isel(y=np.arange(first_row, last_row), x=np.arange(first_col, last_col))
+    return ds 

From 556c769cd04e10357473be9b750d5939918f2ef7 Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Sat, 13 Jul 2024 18:17:38 -0700
Subject: [PATCH 09/41] Allow specifying which components of q to output on
 each fgout grid

Modifications to fgout_module.f90 and fgout_tools.py, primarily.
Can now set `fgout.q_out_vars` to a list of which components of q
to output for each fgout frame (Fortran indexing). The default is
`fgout.q_out_vars = [1,2,3,4]` for all components of q[1:3] and
also eta (4), consistent with the previous behavior.

The user could also/instead ask to output the topo B
as component 5.  (If two out of three of h, eta, B are output then the
other can be computed from these.)

Note that this list is written out to `fgout_grids.data` rather than a
boolean list of True/False values for each possible component as is done
in other places in GeoClaw.
---
 src/2d/shallow/amr2.f90           |   4 +-
 src/2d/shallow/fgout_module.f90   | 198 ++++++++++--------------------
 src/python/geoclaw/data.py        |   2 +-
 src/python/geoclaw/fgout_tools.py | 122 ++++++++++++------
 4 files changed, 147 insertions(+), 179 deletions(-)

diff --git a/src/2d/shallow/amr2.f90 b/src/2d/shallow/amr2.f90
index 139e01f1e..7dea01763 100644
--- a/src/2d/shallow/amr2.f90
+++ b/src/2d/shallow/amr2.f90
@@ -491,7 +491,7 @@ program amr2
         call read_dtopo_settings()        ! specifies file with dtopo from earthquake
         call read_topo_settings(rest)     ! specifies topography (bathymetry) files
         call set_qinit()                  ! specifies file with dh if this used instead
-        call set_fgout(rest)            ! Fixed grid settings
+        call set_fgout(rest,nvar)         ! Fixed grid settings
         call setup_variable_friction()    ! Variable friction parameter
         !call set_multilayer()             ! Set multilayer SWE parameters
         call set_storm()                  ! Set storm parameters
@@ -526,7 +526,7 @@ program amr2
         call read_dtopo_settings()        ! specifies file with dtopo from earthquake
         call read_topo_settings(rest)     ! specifies topography (bathymetry) files
         call set_qinit()                  ! specifies file with dh if this used instead
-        call set_fgout(rest)            ! Fixed grid settings
+        call set_fgout(rest,nvar)         ! Fixed grid settings
         call setup_variable_friction()    ! Variable friction parameter
         call set_multilayer()             ! Set multilayer SWE parameters
         call set_storm()                  ! Set storm parameters
diff --git a/src/2d/shallow/fgout_module.f90 b/src/2d/shallow/fgout_module.f90
index 981d11ac3..b583df59e 100644
--- a/src/2d/shallow/fgout_module.f90
+++ b/src/2d/shallow/fgout_module.f90
@@ -10,7 +10,7 @@ module fgout_module
         real(kind=8), pointer :: late(:,:,:)
 
         ! Geometry
-        integer :: num_vars,mx,my,point_style,fgno,output_format,output_style
+        integer :: mx,my,point_style,fgno,output_format,output_style
         real(kind=8) :: dx,dy,x_low,x_hi,y_low,y_hi
 
         ! Time Tracking and output types
@@ -19,11 +19,13 @@ module fgout_module
 
         integer, allocatable :: output_frames(:)
         real(kind=8), allocatable :: output_times(:)
-        
+
+        integer :: num_vars ! number of variables to interpolate (num_eqn+3)
         integer :: nqout ! number of q components to output (+1 for eta)
-        logical, allocatable :: iqout(:)  ! which components to output
-        integer :: bathy_index,eta_index
-        logical :: dclaw  ! False for GeoClaw
+        logical, allocatable :: q_out_vars(:) ! which components to print
+
+        !logical, allocatable :: iqout(:)  ! which components to output
+        integer :: bathy_index,eta_index,time_index
 
     end type fgout_grid
 
@@ -43,14 +45,16 @@ module fgout_module
     ! Setup routine that reads in the fixed grids data file and sets up the
     ! appropriate data structures
 
-    subroutine set_fgout(rest,fname)
+    subroutine set_fgout(rest,num_eqn,fname)
 
         use amr_module, only: parmunit, tstart_thisrun
+        use utility_module, only: parse_values
 
         implicit none
 
         ! Subroutine arguments
-        logical :: rest  ! restart?
+        logical, intent(in) :: rest  ! restart?
+        integer, intent(in) :: num_eqn
         character(len=*), optional, intent(in) :: fname
 
         ! Local storage
@@ -59,6 +63,9 @@ subroutine set_fgout(rest,fname)
         type(fgout_grid), pointer :: fg
         real(kind=8) :: ts
 
+        integer :: n, iq
+        real(kind=8) :: values(num_eqn+2)
+        character(len=80) :: str
 
         if (.not.module_setup) then
 
@@ -109,7 +116,25 @@ subroutine set_fgout(rest,fname)
                 read(unit,*) fg%mx, fg%my
                 read(unit,*) fg%x_low, fg%y_low
                 read(unit,*) fg%x_hi, fg%y_hi
-                read(unit,*) fg%dclaw
+
+                allocate(fg%q_out_vars(num_eqn+2))  ! for q + eta, topo
+                do iq=1,num_eqn+2
+                    fg%q_out_vars(iq) = .false.
+                enddo
+                read(unit,'(a)') str
+                call parse_values(str, n, values)
+                do k=1,n
+                    iq = nint(values(k))
+                    fg%q_out_vars(iq) = .true.
+                enddo
+                write(6,*) '+++ q_out_vars:', fg%q_out_vars
+
+                fg%num_vars = num_eqn + 3  ! also interp topo,eta,time
+                fg%nqout = 0  ! count how many are to be output
+                do k=1,num_eqn+2
+                    if (fg%q_out_vars(k)) fg%nqout = fg%nqout + 1
+                enddo
+                !fg%nqout = fg%nqout + 1  ! for eta
 
 
                 ! Initialize next_output_index
@@ -193,68 +218,15 @@ subroutine set_fgout(rest,fname)
                     print *,'y grid points my <= 0, set my >= 1'
                 endif
 
-                
-                ! For now, hard-wire with defaults for either GeoClaw or D-Claw
-                ! need to save q plus topo, eta, t for interp in space-time
-                
-                if (fg%dclaw) then
-                    ! For D-Claw:
-                    fg%num_vars = 10
-                    ! for h, hu, hv, hm, pb, hchi, b_eroded, bathy, eta, time
-                else
-                    ! GeoClaw:
-                    fg%num_vars = 6
-                    ! for h, hu, hv, bathy, eta, time
-                endif
-                
-                ! specify which components of q (plus eta?) to output:
-                ! (eventually this should be set from user input)
-
-                if (fg%num_vars == 6) then
-                    ! GeoClaw
-                    ! indexes used in early and late arrays:
-                    ! 1:3 are q variables, 6 is time
-                    fg%bathy_index = 4
-                    fg%eta_index = 5
-
-                    allocate(fg%iqout(4))
-                    fg%iqout(1) = .true.  ! output h?
-                    fg%iqout(2) = .true.  ! output hu?
-                    fg%iqout(3) = .true.  ! output hv?
-                    fg%iqout(4) = .true.  ! output eta?
-                    fg%nqout = 0
-                    do k=1,4
-                        if (fg%iqout(k)) fg%nqout = fg%nqout + 1
-                    enddo
-                endif
+                fg%bathy_index = num_eqn + 2
+                fg%eta_index = num_eqn + 1
+                fg%time_index = num_eqn + 3
 
-                if (fg%num_vars == 10) then
-                    ! D-Claw:
-                    ! indexes used in early and late arrays:
-                    ! 1:7 are q variables, 10 is time
-                    fg%bathy_index = 8
-                    fg%eta_index = 9
-
-                    allocate(fg%iqout(8))
-                    fg%iqout(1) = .true.  ! output h?
-                    fg%iqout(2) = .true.  ! output hu?
-                    fg%iqout(3) = .true.  ! output hv?
-                    fg%iqout(4) = .true.  ! output hm?
-                    fg%iqout(5) = .true.  ! output pb?
-                    fg%iqout(6) = .true.  ! output hchi?
-                    fg%iqout(7) = .true.  ! output beroded?
-                    fg%iqout(8) = .true.  ! output eta?
-                    fg%nqout = 0
-                    do k=1,8
-                        if (fg%iqout(k)) fg%nqout = fg%nqout + 1
-                    enddo
-                endif
-                    
                 write(6,*) '+++ nqout = ',fg%nqout
 
                 ! Allocate new fixed grid data arrays at early, late time:
                 ! dimension (10,:,:) to work for either GeoClaw or D-Claw
-                
+
                 allocate(fg%early(10, fg%mx,fg%my))
                 fg%early = nan()
 
@@ -412,7 +384,7 @@ subroutine fgout_interp(fgrid_type,fgrid, &
 
 
                     ! save the time this fgout point was computed:
-                    fg_data(num_vars,ifg,jfg) = t
+                    fg_data(fgrid%time_index,ifg,jfg) = t
 
 
                 endif ! if fgout point is on this grid
@@ -441,7 +413,7 @@ subroutine fgout_write(fgrid,out_time,out_index)
         ! I/O
         integer, parameter :: unit = 87
         character(len=15) :: fg_filename
-        character(len=4) :: cfgno, cframeno
+        character(len=8) :: cfgno, cframeno
         character(len=8) :: file_format
         integer :: grid_number,ipos,idigit,out_number,columns
         integer :: ifg,ifg1, iframe,iframe1
@@ -472,7 +444,7 @@ subroutine fgout_write(fgrid,out_time,out_index)
                                            "a10,'             format'/,/)"
 
         ! Other locals
-        integer :: i,j,m,iq,k
+        integer :: i,j,m,iq,k,jq,num_eqn
         real(kind=8) :: t0,tf,tau, qaug(10)
         real(kind=8), allocatable :: qeta(:,:,:)
         real(kind=4), allocatable :: qeta4(:,:,:)
@@ -546,76 +518,34 @@ subroutine fgout_write(fgrid,out_time,out_index)
                     endif
                 endif
 
-                ! Output the conserved quantities and eta value
+                ! Output the requested quantities and eta value:
+
                 iq = 1
-                ! qaug(1:3) are h,hu,hv for both GeoClaw and D-Claw:
-                if (fgrid%iqout(1)) then
-                    qeta(iq,i,j) = qaug(1) ! h
-                    iq = iq+1
-                endif
-                if (fgrid%iqout(2)) then
-                    qeta(iq,i,j) = qaug(2) ! hu
-                    iq = iq+1
-                endif
-                if (fgrid%iqout(3)) then
-                    qeta(iq,i,j) = qaug(3) ! hv
-                    iq = iq+1
-                endif
-                
-                if (fgrid%num_vars == 6) then
-                    ! GeoClaw:
-                    if (fgrid%iqout(4)) then
-                        qeta(iq,i,j) = qaug(5) ! eta since qaug(4)=topo
-                        iq = iq+1
-                    endif
-                    
-                else if (fgrid%num_vars == 10) then
-                    ! D-Claw:
-                    if (fgrid%iqout(4)) then
-                        qeta(iq,i,j) = qaug(4) ! hm
-                        iq = iq+1
-                    endif
-                    if (fgrid%iqout(5)) then
-                        qeta(iq,i,j) = qaug(5) ! pb
-                        iq = iq+1
-                    endif
-                    if (fgrid%iqout(6)) then
-                        qeta(iq,i,j) = qaug(6) ! hchi
-                        iq = iq+1
-                    endif
-                    if (fgrid%iqout(7)) then
-                        qeta(iq,i,j) = qaug(7) ! b_eroded
-                        iq = iq+1
-                    endif
-                    if (fgrid%iqout(8)) then
-                        qeta(iq,i,j) = qaug(9) ! eta since qaug(8)=topo
+                num_eqn = size(fgrid%q_out_vars) - 2  ! last components eta,B
+                do jq=1,num_eqn+2
+                    if (fgrid%q_out_vars(jq)) then
+                        qeta(iq,i,j) = qaug(jq)
                         iq = iq+1
                     endif
-                endif
+                enddo
+
+                ! now append eta value: (this is now included in loop above)
+                !qeta(iq,i,j) = qaug(fgrid%eta_index)
+                ! NOTE: qaug(fgrid%bathy_index) contains topo B value
+
             enddo
         enddo
 
 
-        ! Make the file names and open output files
-        cfgno = '0000'
-        ifg = fgrid%fgno
-        ifg1 = ifg
-        do ipos=4,1,-1
-            idigit = mod(ifg1,10)
-            cfgno(ipos:ipos) = char(ichar('0') + idigit)
-            ifg1 = ifg1/10
-            enddo
+        ! convert fgrid%fgno to a string of length at least 4 with leading 0's
+        ! e.g. '0001' or '12345':
+        write(cfgno, '(I0.4)') fgrid%fgno
 
-        cframeno = '0000'
-        iframe = out_index
-        iframe1 = iframe
-        do ipos=4,1,-1
-            idigit = mod(iframe1,10)
-            cframeno(ipos:ipos) = char(ichar('0') + idigit)
-            iframe1 = iframe1/10
-            enddo
+        ! convert out_index to a string of length at least 4 with leading 0's
+        write(cframeno, '(I0.4)') out_index
 
-        fg_filename = 'fgout' // cfgno // '.q' // cframeno
+
+        fg_filename = 'fgout' // trim(cfgno) // '.q' // trim(cframeno)
 
         open(unit,file=fg_filename,status='unknown',form='formatted')
 
@@ -637,14 +567,14 @@ subroutine fgout_write(fgrid,out_time,out_index)
 
         if (fgrid%output_format == 3) then
             ! binary output goes in .b file as full 8-byte (float64):
-            fg_filename = 'fgout' // cfgno // '.b' // cframeno
+            fg_filename = 'fgout' // trim(cfgno) // '.b' // trim(cframeno)
             open(unit=unit, file=fg_filename, status="unknown",    &
                  access='stream')
             write(unit) qeta
             close(unit)
         else if (fgrid%output_format == 2) then
             ! binary output goes in .b file as 4-byte (float32):
-            fg_filename = 'fgout' // cfgno // '.b' // cframeno
+            fg_filename = 'fgout' // trim(cfgno) // '.b' // trim(cframeno)
             open(unit=unit, file=fg_filename, status="unknown",    &
                  access='stream')
             allocate(qeta4(fgrid%nqout, fgrid%mx, fgrid%my))
@@ -671,7 +601,7 @@ subroutine fgout_write(fgrid,out_time,out_index)
             write(6,*) '*** should be ascii, binary32, or binary64'
         endif
 
-        fg_filename = 'fgout' // cfgno // '.t' // cframeno
+        fg_filename = 'fgout' // trim(cfgno) // '.t' // trim(cframeno)
         open(unit=unit, file=fg_filename, status='unknown', form='formatted')
         ! time, num_eqn+1, num_grids, num_aux, num_dim, num_ghost:
         write(unit, t_file_format) out_time, fgrid%nqout, 1, 0, 2, 0,file_format
@@ -680,10 +610,6 @@ subroutine fgout_write(fgrid,out_time,out_index)
         print "(a,i4,a,i4,a,e18.8)",'Writing fgout grid #',fgrid%fgno, &
               '  frame ',out_index,' at time =',out_time
 
-        ! Index into qeta for binary output
-        ! Note that this implicitly assumes that we are outputting only h, hu, hv
-        ! and will not output more (change num_eqn parameter above)
-
     end subroutine fgout_write
 
 
diff --git a/src/python/geoclaw/data.py b/src/python/geoclaw/data.py
index 83d9fa769..01c132400 100755
--- a/src/python/geoclaw/data.py
+++ b/src/python/geoclaw/data.py
@@ -9,7 +9,7 @@
  - GeoClawData
  - RefinementData
  - TopographyData
- - FixedGridData
+ - FGoutData
  - FGmaxData
  - DTopoData
  - QinitData
diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py
index ed8f675c6..9c6c4d2cb 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -41,6 +41,10 @@ def __init__(self, fgout_grid, frameno=None):
         self.frameno = frameno
         self.t = None
 
+        # mapping from variable names to possible values in q_out_vars
+        # default for GeoClaw:
+        self.qmap = {'h':1, 'hu':2, 'hv':3, 'eta':4, 'B':5}
+
         # private attributes for those that are only created if
         # needed by the user:
         self._h = None
@@ -96,16 +100,31 @@ def drytol(self):
     def h(self):
         """depth"""
         if self._h is None:
-            #print('+++ setting _h...')
-            self._h = self.q[0,:,:]
-        #print('+++ getting _h...')
+            q_out_vars = self.fgout_grid.q_out_vars
+            try:
+                i_h = q_out_vars.index(self.qmap['h'])
+                self._h = self.q[i_h,:,:]
+            except:
+                try:
+                    i_eta = q_out_vars.index(self.qmap['eta'])
+                    i_B = q_out_vars.index(self.qmap['B'])
+                    self._h = self.q[i_eta,:,:] - self.q[i_B,:,:]
+                except:
+                    print('*** Could not find h or eta-B in q_out_vars')
+                    raise
         return self._h
 
     @property
     def hu(self):
         """momentum h*u"""
         if self._hu is None:
-            self._hu = self.q[1,:,:]
+            q_out_vars = self.fgout_grid.q_out_vars
+            try:
+                i_hu = q_out_vars.index(self.qmap['hu'])
+                self._hu = self.q[i_hu,:,:]
+            except:
+                print('*** Could not find hu in q_out_vars')
+                raise
         return self._hu
 
     @property
@@ -121,7 +140,13 @@ def u(self):
     def hv(self):
         """momentum h*v"""
         if self._hv is None:
-            self._hv = self.q[2,:,:]
+            q_out_vars = self.fgout_grid.q_out_vars
+            try:
+                i_hv = q_out_vars.index(self.qmap['hv'])
+                self._hv = self.q[i_hv,:,:]
+            except:
+                print('*** Could not find hv in q_out_vars')
+                raise
         return self._hv
 
     @property
@@ -137,14 +162,36 @@ def v(self):
     def eta(self):
         """surface eta = h+B"""
         if self._eta is None:
-            self._eta = self.q[-1,:,:]
+            q_out_vars = self.fgout_grid.q_out_vars
+            try:
+                i_eta = q_out_vars.index(self.qmap['eta'])
+                self._eta = self.q[i_eta,:,:]
+            except:
+                try:
+                    i_h = q_out_vars.index(self.qmap['h'])
+                    i_B = q_out_vars.index(self.qmap['B'])
+                    self._eta = self.q[i_h,:,:] + self.q[i_B,:,:]
+                except:
+                    print('*** Could not find eta or h+B in q_out_vars')
+                    raise
         return self._eta
 
     @property
     def B(self):
         """topography"""
         if self._B is None:
-            self._B = self.q[-1,:,:] - self.q[0,:,:]
+            q_out_vars = self.fgout_grid.q_out_vars
+            try:
+                i_B = q_out_vars.index(self.qmap['B'])
+                self._B = self.q[i_B,:,:]
+            except:
+                try:
+                    i_h = q_out_vars.index(self.qmap['h'])
+                    i_eta = q_out_vars.index(self.qmap['eta'])
+                    self._B = self.q[i_eta,:,:] - self.q[i_h,:,:]
+                except:
+                    print('*** Could not find B or eta-h in q_out_vars')
+                    raise
         return self._B
 
     @property
@@ -169,7 +216,7 @@ class FGoutGrid(object):
     fgout input data and the output generated by a GeoClaw run.
     """
 
-    def __init__(self,fgno=None,outdir=None,output_format=None):
+    def __init__(self,fgno=None,outdir='.',output_format=None):
 
 
         # GeoClaw input values:
@@ -185,7 +232,9 @@ def __init__(self,fgno=None,outdir=None,output_format=None):
         self.fgno = fgno
         self.outdir = outdir
         self.output_format = output_format
-        self.dclaw = False
+        self.q_out_vars = [1,2,3,4]  # list of which components to print
+                                          # default: h,hu,hv,eta (5=topo)
+        self.nqout = None # number of vars to print out
         
         self.drytol = 1e-3          # used for computing u,v from hu,hv
 
@@ -312,6 +361,9 @@ def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'):
             self.fgno = fgno
         assert self.fgno is not None, '*** fgno must be set'
 
+        data_path = os.path.join(self.outdir, data_file)
+        print('Reading fgout grid info from \n    %s' % data_path)
+        
         with open(data_file) as filep:
             lines = filep.readlines()
         fgout_input = None
@@ -371,6 +423,13 @@ def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'):
         else:
             raise NotImplementedError("fgout not implemented for point_style %i" \
                 % point_style)
+        tokens = fgout_input[lineno].split()
+        self.q_out_vars = []
+        for token in tokens:
+            try:
+                self.q_out_vars.append(int(token))
+            except:
+                break
 
 
     def write_to_fgout_data(self, fid):
@@ -455,7 +514,11 @@ def write_to_fgout_data(self, fid):
             print("   upper right = (%15.10f,%15.10f)" % (x2,y2))
             print("   dx = %15.10e,  dy = %15.10e" % (dx,dy))
             
-            fid.write("%s                     # dclaw" % self.dclaw)
+        # q_out_vars is a list of q components to print, e.g. [1,4]
+
+        format = len(self.q_out_vars) * '%s '
+        fid.write(format % tuple(self.q_out_vars)+ " # q_out_vars\n")
+        fid.write('\n')
 
 
     def read_frame(self, frameno):
@@ -465,6 +528,14 @@ def read_frame(self, frameno):
 
         from datetime import timedelta
 
+        try:
+            fgoutX = self.X
+            fgoutY = self.Y
+        except:
+            self.read_fgout_grids_data()
+            fgoutX = self.X
+            fgoutY = self.Y
+            
         try:
             fr = self.plotdata.getframe(frameno)
         except:
@@ -483,36 +554,7 @@ def read_frame(self, frameno):
     
         fgout_frame.frameno = frameno
         
-        X,Y = patch.grid.p_centers[:2]
-        
-        try:
-            fgoutX = self.X
-            fgoutY = self.Y
-        except:
-            # presumably x,y, etc. were not set for this fgout_grid
-            # (e.g. by read_fgout_grids_data) so infer from this frame and set
-            
-            # reset most attributes including private _x etc to None:
-            self.__init__(fgno=self.fgno,outdir=self.outdir,
-                          output_format=self.output_format)
-                          
-            print('    Setting grid attributes based on Frame %i:' % frameno)
-            x = X[:,0]
-            y = Y[0,:]
-            dx = x[1] - x[0]
-            dy = y[1] - y[0]
-            self.x1 = x[0] - dx/2
-            self.x2 = x[-1] + dx/2
-            self.y1 = y[0] - dy/2
-            self.y2 = y[-1] + dy/2
-            self.nx = len(x)
-            self.ny = len(y)
-            fgoutX = self.X  # will be computed from info above
-            fgoutY = self.Y  # will be computed from info above
-            print('      Grid edges: ',self.extent_edges)
-            print('      with %i cells in x and %i cells in y' \
-                    % (self.nx,self.ny))
-            
+        X,Y = patch.grid.p_centers[:2]    
             
         if not numpy.allclose(X, fgoutX):
             errmsg = '*** X read from output does not match fgout_grid.X'

From b5c898a0a5189e2e0fe0e8463e0454b9f04bcb21 Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Tue, 16 Jul 2024 22:01:51 -0700
Subject: [PATCH 10/41] working on dclaw variables

---
 src/python/geoclaw/fgout_tools.py | 44 +++++++++++++++++++++++++++----
 1 file changed, 39 insertions(+), 5 deletions(-)

diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py
index 9c6c4d2cb..dfedda3fc 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -36,14 +36,25 @@ class FGoutFrame(object):
     and stored only when needed by the user.
     """
 
-    def __init__(self, fgout_grid, frameno=None):
+    def __init__(self, fgout_grid, frameno=None, qmap='geoclaw'):
         self.fgout_grid = fgout_grid
         self.frameno = frameno
         self.t = None
 
         # mapping from variable names to possible values in q_out_vars
         # default for GeoClaw:
-        self.qmap = {'h':1, 'hu':2, 'hv':3, 'eta':4, 'B':5}
+        if type(qmap) is dict:
+            self.qmap = qmap
+        elif qmap == 'geoclaw':
+            self.qmap = {'h':1, 'hu':2, 'hv':3, 'eta':4, 'B':5}
+        elif qmap == 'geoclaw-bouss':
+            self.qmap = {'h':1, 'hu':2, 'hv':3, 'huc':4, 'hvc':5,
+                         'eta':4, 'B':5}
+        elif qmap == 'dclaw':
+            self.qmap = {'h':1, 'hu':2, 'hv':3, 'hm':4, 
+                         'eta':8, 'B':9}
+        else:
+            raise InputError('Invalid qmap: %s' % qmap)
 
         # private attributes for those that are only created if
         # needed by the user:
@@ -56,6 +67,7 @@ def __init__(self, fgout_grid, frameno=None):
         self._v = None
         self._s = None
         self._hss = None
+        self._hm = None
 
     # Define shortcuts to attributes of self.fgout_grid that are the same
     # for all frames (e.g. X,Y) to avoid storing grid for every frame.
@@ -166,6 +178,8 @@ def eta(self):
             try:
                 i_eta = q_out_vars.index(self.qmap['eta'])
                 self._eta = self.q[i_eta,:,:]
+                print('+++qmap["eta"] = %i' % self.qmap["eta"])
+                print('+++i_eta = %i' % i_eta)
             except:
                 try:
                     i_h = q_out_vars.index(self.qmap['h'])
@@ -184,11 +198,16 @@ def B(self):
             try:
                 i_B = q_out_vars.index(self.qmap['B'])
                 self._B = self.q[i_B,:,:]
+                print('+++qmap["B"] = %i' % self.qmap["B"])
+                print('+++i_B = %i' % i_B)
             except:
                 try:
                     i_h = q_out_vars.index(self.qmap['h'])
                     i_eta = q_out_vars.index(self.qmap['eta'])
                     self._B = self.q[i_eta,:,:] - self.q[i_h,:,:]
+                    print('+++ computing B: i_h = %i, i_eta = %i' % (i_h,i_eta))
+                    print('+++qmap["h"] = %i' % self.qmap["h"])
+                    print('+++qmap["eta"] = %i' % self.qmap["eta"])
                 except:
                     print('*** Could not find B or eta-h in q_out_vars')
                     raise
@@ -208,6 +227,21 @@ def hss(self):
             self._hss = self.h * self.s**2
         return self._hss
 
+    @property
+    def hm(self):
+        """dclaw: h * mass fraction"""
+        if self._hm is None:
+            q_out_vars = self.fgout_grid.q_out_vars
+            try:
+                i_hm = q_out_vars.index(self.qmap['hm'])
+                self._hm = self.q[i_hm,:,:]
+                print('+++qmap["hm"] = %i' % self.qmap["hm"])
+                print('+++i_hm = %i' % i_hm)
+            except:
+                print('*** Could not find hm in q_out_vars')
+                raise
+        return self._hm
+
 
 class FGoutGrid(object):
 
@@ -376,7 +410,7 @@ def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'):
 
         if fgout_input is None:
             raise ValueError('fgout grid fgno = %i not found in %s' \
-                             % (fgno, data_file))
+                             % (self.fgno, data_file))
 
         lineno = 0 # next input line
         self.output_style = int(fgout_input[lineno].split()[lineno])
@@ -521,7 +555,7 @@ def write_to_fgout_data(self, fid):
         fid.write('\n')
 
 
-    def read_frame(self, frameno):
+    def read_frame(self, frameno, qmap='geoclaw'):
         """
         Read a single frame of fgout data.
         """
@@ -545,7 +579,7 @@ def read_frame(self, frameno):
         state = fr.states[0]  # only 1 AMR grid
         patch = state.patch
 
-        fgout_frame = FGoutFrame(self, frameno)
+        fgout_frame = FGoutFrame(self, frameno, qmap)
         fgout_frame.fgout_grid = self
 
         fgout_frame.q = state.q

From 70ee91e85dbb40d3959f7577183b909b3b525c13 Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Wed, 17 Jul 2024 17:53:45 -0700
Subject: [PATCH 11/41] Fix fgout_tools.py to support fgout.q_out_vars also for
 Boussinesq or D-Claw

When instantiating fgout_tools.FGoutGrid object, new qmap parameter can be set
to 'geoclaw', 'geoclaw-bouss', 'dclaw', or a custom dictionary.
(Eventually add 'dclaw-bouss' option also).
This dictionary has the form `{'h':1, 'hu':2, etc}` with mapping from fgout
variable names to Fortran indices of q in the code creating the fgout output.
---
 src/python/geoclaw/fgout_tools.py | 129 +++++++++++++++++++++++-------
 1 file changed, 101 insertions(+), 28 deletions(-)

diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py
index dfedda3fc..c53a0240e 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -36,26 +36,11 @@ class FGoutFrame(object):
     and stored only when needed by the user.
     """
 
-    def __init__(self, fgout_grid, frameno=None, qmap='geoclaw'):
+    def __init__(self, fgout_grid, frameno=None):
         self.fgout_grid = fgout_grid
         self.frameno = frameno
         self.t = None
 
-        # mapping from variable names to possible values in q_out_vars
-        # default for GeoClaw:
-        if type(qmap) is dict:
-            self.qmap = qmap
-        elif qmap == 'geoclaw':
-            self.qmap = {'h':1, 'hu':2, 'hv':3, 'eta':4, 'B':5}
-        elif qmap == 'geoclaw-bouss':
-            self.qmap = {'h':1, 'hu':2, 'hv':3, 'huc':4, 'hvc':5,
-                         'eta':4, 'B':5}
-        elif qmap == 'dclaw':
-            self.qmap = {'h':1, 'hu':2, 'hv':3, 'hm':4, 
-                         'eta':8, 'B':9}
-        else:
-            raise InputError('Invalid qmap: %s' % qmap)
-
         # private attributes for those that are only created if
         # needed by the user:
         self._h = None
@@ -103,6 +88,10 @@ def extent_edges(self):
     @property
     def drytol(self):
         return self.fgout_grid.drytol
+        
+    @property
+    def qmap(self):
+        return self.fgout_grid.qmap
             
     # Define attributes such as h as @properties with lazy evaluation:
     # the corresponding array is created and stored only when first
@@ -178,8 +167,8 @@ def eta(self):
             try:
                 i_eta = q_out_vars.index(self.qmap['eta'])
                 self._eta = self.q[i_eta,:,:]
-                print('+++qmap["eta"] = %i' % self.qmap["eta"])
-                print('+++i_eta = %i' % i_eta)
+                #print('+++qmap["eta"] = %i' % self.qmap["eta"])
+                #print('+++i_eta = %i' % i_eta)
             except:
                 try:
                     i_h = q_out_vars.index(self.qmap['h'])
@@ -198,16 +187,16 @@ def B(self):
             try:
                 i_B = q_out_vars.index(self.qmap['B'])
                 self._B = self.q[i_B,:,:]
-                print('+++qmap["B"] = %i' % self.qmap["B"])
-                print('+++i_B = %i' % i_B)
+                #print('+++qmap["B"] = %i' % self.qmap["B"])
+                #print('+++i_B = %i' % i_B)
             except:
                 try:
                     i_h = q_out_vars.index(self.qmap['h'])
                     i_eta = q_out_vars.index(self.qmap['eta'])
                     self._B = self.q[i_eta,:,:] - self.q[i_h,:,:]
-                    print('+++ computing B: i_h = %i, i_eta = %i' % (i_h,i_eta))
-                    print('+++qmap["h"] = %i' % self.qmap["h"])
-                    print('+++qmap["eta"] = %i' % self.qmap["eta"])
+                    #print('+++ computing B: i_h = %i, i_eta = %i' % (i_h,i_eta))
+                    #print('+++qmap["h"] = %i' % self.qmap["h"])
+                    #print('+++qmap["eta"] = %i' % self.qmap["eta"])
                 except:
                     print('*** Could not find B or eta-h in q_out_vars')
                     raise
@@ -227,6 +216,32 @@ def hss(self):
             self._hss = self.h * self.s**2
         return self._hss
 
+    @property
+    def huc(self):
+        """huc - Boussinesq correction to hu"""
+        if self._huc is None:
+            q_out_vars = self.fgout_grid.q_out_vars
+            try:
+                i_huc = q_out_vars.index(self.qmap['huc'])
+                self._huc = self.q[i_huc,:,:]
+            except:
+                print('*** Could not find huc in q_out_vars')
+                raise
+        return self._huc
+
+    @property
+    def hvc(self):
+        """hvc - Boussinesq correction to hv"""
+        if self._hvc is None:
+            q_out_vars = self.fgout_grid.q_out_vars
+            try:
+                i_hvc = q_out_vars.index(self.qmap['hvc'])
+                self._hvc = self.q[i_hvc,:,:]
+            except:
+                print('*** Could not find hvc in q_out_vars')
+                raise
+        return self._hvc
+
     @property
     def hm(self):
         """dclaw: h * mass fraction"""
@@ -235,13 +250,52 @@ def hm(self):
             try:
                 i_hm = q_out_vars.index(self.qmap['hm'])
                 self._hm = self.q[i_hm,:,:]
-                print('+++qmap["hm"] = %i' % self.qmap["hm"])
-                print('+++i_hm = %i' % i_hm)
+                #print('+++qmap["hm"] = %i' % self.qmap["hm"])
+                #print('+++i_hm = %i' % i_hm)
             except:
                 print('*** Could not find hm in q_out_vars')
                 raise
         return self._hm
 
+    @property
+    def pb(self):
+        """dclaw variable """
+        if self._pb is None:
+            q_out_vars = self.fgout_grid.q_out_vars
+            try:
+                i_pb = q_out_vars.index(self.qmap['pb'])
+                self._pb = self.q[i_pb,:,:]
+            except:
+                print('*** Could not find pb in q_out_vars')
+                raise
+        return self._pb
+
+    @property
+    def hchi(self):
+        """dclaw variable """
+        if self._hchi is None:
+            q_out_vars = self.fgout_grid.q_out_vars
+            try:
+                i_hchi = q_out_vars.index(self.qmap['hchi'])
+                self._hchi = self.q[i_hchi,:,:]
+            except:
+                print('*** Could not find hchi in q_out_vars')
+                raise
+        return self._hchi
+
+    @property
+    def bdif(self):
+        """dclaw variable """
+        if self._bdif is None:
+            q_out_vars = self.fgout_grid.q_out_vars
+            try:
+                i_bdif = q_out_vars.index(self.qmap['bdif'])
+                self._bdif = self.q[i_bdif,:,:]
+            except:
+                print('*** Could not find bdif in q_out_vars')
+                raise
+        return self._bdif
+
 
 class FGoutGrid(object):
 
@@ -250,9 +304,25 @@ class FGoutGrid(object):
     fgout input data and the output generated by a GeoClaw run.
     """
 
-    def __init__(self,fgno=None,outdir='.',output_format=None):
+    def __init__(self,fgno=None,outdir='.',output_format=None,
+                 qmap='geoclaw'):
 
 
+        # mapping from variable names to possible values in q_out_vars
+        if type(qmap) is dict:
+            self.qmap = qmap
+        elif qmap == 'geoclaw':
+            # default for GeoClaw:
+            self.qmap = {'h':1, 'hu':2, 'hv':3, 'eta':4, 'B':5}
+        elif qmap == 'geoclaw-bouss':
+            self.qmap = {'h':1, 'hu':2, 'hv':3, 'huc':4, 'hvc':5,
+                         'eta':4, 'B':5}
+        elif qmap == 'dclaw':
+            self.qmap = {'h':1, 'hu':2, 'hv':3, 'hm':4, 'pb':5, 'hchi':6,
+                         'bdif':7, 'eta':8, 'B':9}
+        else:
+            raise InputError('Invalid qmap: %s' % qmap)
+        
         # GeoClaw input values:
         self.id = ''  # identifier, optional
         self.point_style = 2  # only option currently supported
@@ -464,6 +534,9 @@ def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'):
                 self.q_out_vars.append(int(token))
             except:
                 break
+        print('Found fgout grid q_out_vars = ',self.q_out_vars)
+        print('Using this mapping to fgout variable names: ')
+        print('      qmap = ',self.qmap)
 
 
     def write_to_fgout_data(self, fid):
@@ -555,7 +628,7 @@ def write_to_fgout_data(self, fid):
         fid.write('\n')
 
 
-    def read_frame(self, frameno, qmap='geoclaw'):
+    def read_frame(self, frameno):
         """
         Read a single frame of fgout data.
         """
@@ -579,7 +652,7 @@ def read_frame(self, frameno, qmap='geoclaw'):
         state = fr.states[0]  # only 1 AMR grid
         patch = state.patch
 
-        fgout_frame = FGoutFrame(self, frameno, qmap)
+        fgout_frame = FGoutFrame(self, frameno)
         fgout_frame.fgout_grid = self
 
         fgout_frame.q = state.q

From 7f4d10cc3bb03e28e058ab06c491254298958dca Mon Sep 17 00:00:00 2001
From: Kyle Mandli <kyle.Mandli@gmail.com>
Date: Thu, 18 Jul 2024 15:05:54 -0400
Subject: [PATCH 12/41] Update topo_module.f90

The `character` variables that were storing variable names and IDs for NetCDF topography reading were only length 10.  This just makes them 64 characters long.  Probably too long but not really a big deal.
---
 src/2d/shallow/topo_module.f90 | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/2d/shallow/topo_module.f90 b/src/2d/shallow/topo_module.f90
index 8bd1af908..2713e1694 100644
--- a/src/2d/shallow/topo_module.f90
+++ b/src/2d/shallow/topo_module.f90
@@ -439,7 +439,7 @@ subroutine read_topo_file(mx,my,topo_type,fname,xll,yll,topo)
         integer(kind=8) :: i, j, mtot
 
         ! NetCDF Support
-        character(len=10) :: direction, x_dim_name, x_var_name, y_dim_name, &
+        character(len=64) :: direction, x_dim_name, x_var_name, y_dim_name, &
             y_var_name, z_var_name, var_name
         ! character(len=1) :: axis_string
         real(kind=8), allocatable :: nc_buffer(:, :), xlocs(:), ylocs(:)
@@ -744,8 +744,8 @@ subroutine read_topo_header(fname,topo_type,mx,my,xll,yll,xhi,yhi,dx,dy)
         logical, allocatable :: x_in_dom(:),y_in_dom(:)
         integer(kind=4) :: dim_ids(2), num_dims, var_type, num_vars, num_dims_tot
         integer(kind=4), allocatable :: var_ids(:)
-        character(len=10) :: var_name, x_var_name, y_var_name, z_var_name
-        character(len=10) :: x_dim_name, y_dim_name
+        character(len=64) :: var_name, x_var_name, y_var_name, z_var_name
+        character(len=64) :: x_dim_name, y_dim_name
         integer(kind=4) :: x_var_id, y_var_id, z_var_id, x_dim_id, y_dim_id
 
         verbose = .false.

From 981ab8c1904a2b88c8554a7b6edb70cc1f6c29a3 Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Wed, 24 Jul 2024 14:40:31 -0700
Subject: [PATCH 13/41] fix data_file -> data_path

---
 src/python/geoclaw/fgout_tools.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py
index c53a0240e..ab8d6f977 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -468,7 +468,7 @@ def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'):
         data_path = os.path.join(self.outdir, data_file)
         print('Reading fgout grid info from \n    %s' % data_path)
         
-        with open(data_file) as filep:
+        with open(data_path) as filep:
             lines = filep.readlines()
         fgout_input = None
         for lineno,line in enumerate(lines):

From 9378f12ffc5e876ecfec5cfc96ac7c0699ac4f96 Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Mon, 29 Jul 2024 06:49:22 -0700
Subject: [PATCH 14/41] add FGoutGrid.read_fgout_grids_data_pre511 for reading
 legacy fgout_grids.data from before v5.11

---
 src/python/geoclaw/fgout_tools.py | 53 +++++++++++++++++++++++++++++++
 1 file changed, 53 insertions(+)

diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py
index ab8d6f977..4ea277ab7 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -538,6 +538,59 @@ def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'):
         print('Using this mapping to fgout variable names: ')
         print('      qmap = ',self.qmap)
 
+    def read_fgout_grids_data_pre511(self, fgno=None,
+                                     data_file='fgout_grids.data'):
+        """
+        For backward compatibility, this reads fgout_grids.data files
+        in the format used prior to v5.11.
+
+        Read input info for fgout grid number fgno from the data file
+        fgout_grids.data, which should have been created by setrun.py.
+        This file now contains info about all fgout grids.
+        """
+
+        if fgno is not None:
+            self.fgno = fgno
+        assert self.fgno is not None, '*** fgno must be set'
+
+        data_path = os.path.join(self.outdir, data_file)
+        print('Reading fgout grid info from \n    %s' % data_path)
+        
+        with open(data_path) as filep:
+            lines = filep.readlines()
+        fgout_input = None
+        for lineno,line in enumerate(lines):
+            if 'fgno' in line:
+                if int(line.split()[0]) == self.fgno:
+                    fgout_input = lines[lineno+1:]
+                    #print('Found line %i: %s' % (lineno,line))
+                    break
+
+        if fgout_input is None:
+            raise ValueError('fgout grid fgno = %i not found in %s' \
+                             % (fgno, data_file))
+
+        self.tstart = float(fgout_input[0].split()[0])
+        self.tend = float(fgout_input[1].split()[0])
+        self.nout = int(fgout_input[2].split()[0])
+        self.point_style = point_style = int(fgout_input[3].split()[0])
+        output_format = int(fgout_input[4].split()[0])
+        if output_format == 1:
+            self.output_format = 'ascii'
+        elif output_format == 3:
+            self.output_format = 'binary'
+        print('Reading input for fgno=%i, point_style = %i ' \
+                % (self.fgno, self.point_style))
+        if point_style == 2:
+            self.nx = nx = int(fgout_input[5].split()[0])
+            self.ny = ny = int(fgout_input[5].split()[1])
+            self.x1 = float(fgout_input[6].split()[0])
+            self.y1 = float(fgout_input[6].split()[1])
+            self.x2 = float(fgout_input[7].split()[0])
+            self.y2 = float(fgout_input[7].split()[1])
+        else:
+            raise NotImplementedError("fgout not implemented for point_style %i" \
+                % point_style)
 
     def write_to_fgout_data(self, fid):
         """

From 665bbdadbb5146b9476d8de38aa920ce5b1fa06f Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Mon, 29 Jul 2024 07:05:26 -0700
Subject: [PATCH 15/41] add to doc string for pre511

---
 src/python/geoclaw/fgout_tools.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py
index 4ea277ab7..a648791aa 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -544,6 +544,11 @@ def read_fgout_grids_data_pre511(self, fgno=None,
         For backward compatibility, this reads fgout_grids.data files
         in the format used prior to v5.11.
 
+        In this case, the following values are used, as set in __init__():
+            self.output_style = 1
+            self.qmap = 'qmap'
+            self.q_out_vars = [1,2,3,4]  # h,hu,hv,eta
+
         Read input info for fgout grid number fgno from the data file
         fgout_grids.data, which should have been created by setrun.py.
         This file now contains info about all fgout grids.

From f373e231005f4df496967c6462d9cc99decc0278 Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Mon, 29 Jul 2024 07:33:31 -0700
Subject: [PATCH 16/41] print more helpful error msg if fgout_grids.data seems
 to be pre511

---
 src/python/geoclaw/fgout_tools.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py
index a648791aa..f4e0eb4d3 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -483,7 +483,14 @@ def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'):
                              % (self.fgno, data_file))
 
         lineno = 0 # next input line
-        self.output_style = int(fgout_input[lineno].split()[lineno])
+        next_value = fgout_input[lineno].split()[lineno]  # a string
+        next_value_int = ('.' not in next_value)  # should be True in v5.11
+        err_msg = '\n*** Expecting integer output_style next in %s' \
+            % data_file \
+            + '\n    If this is fgout data from before v5.11, try using' \
+            + '\n    read_fgout_grids_data_pre511'
+        assert next_value_int, err_msg
+        self.output_style = int(next_value)
         lineno += 1
         if (self.output_style == 1):
             # equally spaced times:

From 6c849f2448339e5da303e3d50aeb0eb45608cdd0 Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Mon, 29 Jul 2024 07:35:16 -0700
Subject: [PATCH 17/41] force user to call read_fgout_grids_data before
 read_frame

---
 src/python/geoclaw/fgout_tools.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py
index f4e0eb4d3..71b30e4f1 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -704,6 +704,14 @@ def read_frame(self, frameno):
             fgoutX = self.X
             fgoutY = self.Y
         except:
+            msg = '\n*** Before reading frame, you must set FGoutGrid data,' \
+                  '\n*** Typically by calling read_fgout_grids_data'
+            raise ValueError(msg)
+
+            # prior to v5.11, self.read_fgout_grids_data() called here
+            # rather than raising exception...
+            print(msg)
+            print('*** Calling read_fgout_grids_data...')
             self.read_fgout_grids_data()
             fgoutX = self.X
             fgoutY = self.Y

From 178403e6d2b0e33a0a55c7000d1ea4c7a52cf4cb Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Mon, 29 Jul 2024 10:08:04 -0700
Subject: [PATCH 18/41] Fix chiles2010_fgmax-fgout scripts to make fgout
 animations

Now need to call fgout_grid.read_fgout_grids_data() explicitly,
but do not need to set format.
---
 .../tsunami/chile2010_fgmax-fgout/make_fgout_animation.py     | 4 ++--
 .../make_fgout_animation_with_transect.py                     | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation.py b/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation.py
index 790e8b746..ab47c7588 100644
--- a/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation.py
+++ b/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation.py
@@ -34,7 +34,6 @@
 fgno = 1  # which fgout grid
 
 outdir = '_output'
-format = 'binary'  # format of fgout grid output
 
 if 1:
     # use all fgout frames in outdir:
@@ -48,7 +47,8 @@
     fgframes = range(1,26)  # frames of fgout solution to use in animation
 
 # Instantiate object for reading fgout frames:
-fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format) 
+fgout_grid = fgout_tools.FGoutGrid(fgno, outdir) 
+fgout_grid.read_fgout_grids_data()
 
 
 # Plot one frame of fgout data and define the Artists that will need to
diff --git a/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation_with_transect.py b/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation_with_transect.py
index 4b32c40da..f2d19caad 100644
--- a/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation_with_transect.py
+++ b/examples/tsunami/chile2010_fgmax-fgout/make_fgout_animation_with_transect.py
@@ -38,7 +38,6 @@
 fgno = 1  # which fgout grid
 
 outdir = '_output'
-format = 'binary'  # format of fgout grid output
 
 if 1:
     # use all fgout frames in outdir:
@@ -52,7 +51,8 @@
     fgframes = range(1,26)  # frames of fgout solution to use in animation
 
 # Instantiate object for reading fgout frames:
-fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, format) 
+fgout_grid = fgout_tools.FGoutGrid(fgno, outdir) 
+fgout_grid.read_fgout_grids_data()
 
 
 # Plot one frame of fgout data and define the Artists that will need to

From 1ec2512069f1b73a1ae7a42a1a5969e721040cb8 Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Mon, 29 Jul 2024 10:09:55 -0700
Subject: [PATCH 19/41] Leave output_format as an argument to
 FGoutGrids.__init__ but not that this is ignored since it is read from
 fgout_grids.data. Also properly handle output_format == 2 for binary32.

---
 src/python/geoclaw/fgout_tools.py | 12 +++++++++++-
 1 file changed, 11 insertions(+), 1 deletion(-)

diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py
index 71b30e4f1..4d2638259 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -321,9 +321,11 @@ def __init__(self,fgno=None,outdir='.',output_format=None,
             self.qmap = {'h':1, 'hu':2, 'hv':3, 'hm':4, 'pb':5, 'hchi':6,
                          'bdif':7, 'eta':8, 'B':9}
         else:
-            raise InputError('Invalid qmap: %s' % qmap)
+            raise ValueError('Invalid qmap: %s' % qmap)
         
         # GeoClaw input values:
+        # Many of these should be read from fgout_grids.data
+        # using read_fgout_grids_data before using read_frame
         self.id = ''  # identifier, optional
         self.point_style = 2  # only option currently supported
         self.npts = None
@@ -335,7 +337,10 @@ def __init__(self,fgno=None,outdir='.',output_format=None,
         self.nout = None
         self.fgno = fgno
         self.outdir = outdir
+
+        # Note output_format will be reset by read_fgout_grids_data:
         self.output_format = output_format
+
         self.q_out_vars = [1,2,3,4]  # list of which components to print
                                           # default: h,hu,hv,eta (5=topo)
         self.nqout = None # number of vars to print out
@@ -517,8 +522,13 @@ def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'):
         lineno += 1
         if output_format == 1:
             self.output_format = 'ascii'
+        elif output_format == 2:
+            self.output_format = 'binary32'
         elif output_format == 3:
             self.output_format = 'binary'
+        else:
+            raise NotImplementedError("fgout not implemented for " \
+                    + "output_format %i"  % output_format)
         print('Reading input for fgno=%i, point_style = %i ' \
                 % (self.fgno, self.point_style))
         if point_style == 2:

From 6e67973f95b6d15cdcd092bd4d40197bb28f4fd3 Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Mon, 29 Jul 2024 11:03:37 -0700
Subject: [PATCH 20/41] fix qmap=='geoclaw-bouss' mapping

---
 src/python/geoclaw/fgout_tools.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py
index 4d2638259..7fdc47e19 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -316,7 +316,7 @@ def __init__(self,fgno=None,outdir='.',output_format=None,
             self.qmap = {'h':1, 'hu':2, 'hv':3, 'eta':4, 'B':5}
         elif qmap == 'geoclaw-bouss':
             self.qmap = {'h':1, 'hu':2, 'hv':3, 'huc':4, 'hvc':5,
-                         'eta':4, 'B':5}
+                         'eta':6, 'B':7}
         elif qmap == 'dclaw':
             self.qmap = {'h':1, 'hu':2, 'hv':3, 'hm':4, 'pb':5, 'hchi':6,
                          'bdif':7, 'eta':8, 'B':9}

From c0ba20482702f5c985c97e0daa4f64fec8b4a98d Mon Sep 17 00:00:00 2001
From: "Barnhart, Katherine (Katy) Ruth" <krbarnhart@usgs.gov>
Date: Tue, 20 Aug 2024 10:38:20 -0600
Subject: [PATCH 21/41] wip

---
 src/python/geoclaw/xarray_backends.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index 8f06f6fa4..d52736e8a 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -536,4 +536,4 @@ def _clip(ds, clip_data):
     last_col = sel_cols[-1]
 
     ds = ds.isel(y=np.arange(first_row, last_row), x=np.arange(first_col, last_col))
-    return ds 
+    return ds

From e21b2947ff0b4307987ab0cdb8b838540d4b149c Mon Sep 17 00:00:00 2001
From: "Barnhart, Katherine (Katy) Ruth" <krbarnhart@usgs.gov>
Date: Tue, 20 Aug 2024 14:05:03 -0600
Subject: [PATCH 22/41] first pass at using revised fgout data structure with
 backends.

---
 .../use_xarray_backends.py                    | 39 ++++++---
 src/python/geoclaw/xarray_backends.py         | 84 ++++++++++++-------
 2 files changed, 83 insertions(+), 40 deletions(-)

diff --git a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
index 94b55b3b8..b0cf65d5a 100644
--- a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
+++ b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
@@ -1,11 +1,13 @@
+import glob
+
 try:
     import rioxarray
     import xarray as xr
 except:
-    'You must install xarray and rioxarray in order to use the xarray backends'
+    "You must install xarray and rioxarray in order to use the xarray backends"
     raise
 
-from clawpack.geoclaw.xarray_backends import FGOutBackend, FGMaxBackend
+from clawpack.geoclaw.xarray_backends import FGMaxBackend, FGOutBackend
 
 # epsg code for lat-lon
 # Optionally, provide an epsg code to assign the associated coordinate system to the file.
@@ -16,22 +18,39 @@
 epsg_code = 4326
 
 # An example of a fgout grid.
-filename = '_output/fgout0001.b0001'
+filename = "_output/fgout0001.b0001"
 # provide the .bxxx file if binary format is used or the
 # .qxxx file if ascii format is used.
 # the format, fg number, and frame number are inferred from the filename.
 
-ds = xr.open_dataset(filename, engine=FGOutBackend, backend_kwargs={'epsg':epsg_code})
+ds = xr.open_dataset(
+    filename, engine=FGOutBackend, backend_kwargs={"epsg": epsg_code, "qmap": "geoclaw"}
+)
 # ds is now an xarray object. It can be interacted with directly or written to netcdf using
-ds.to_netcdf('fgout0001_0001.nc')
+ds.to_netcdf("fgout0001_0001.nc")
+
+# It is possible to combine all fgout files into a single netcdf file
+# using xr.open_mfdataset
+# https://docs.xarray.dev/en/stable/generated/xarray.open_mfdataset.html
+fgout_files = glob.glob("_output/fgout0001.b*")
+ds_all = xr.open_mfdataset(
+    fgout_files,
+    engine=FGOutBackend,
+    backend_kwargs={"epsg": epsg_code, "qmap": "geoclaw"},
+)
+
+print(ds_all)
+ds_all.to_netcdf("fgout_all.nc")
 
 # An example of a fgmax grid.
 filename = "_output/fgmax0001.txt"
-ds = xr.open_dataset(filename, engine=FGMaxBackend, backend_kwargs={'epsg':epsg_code})
-ds.to_netcdf('fgmax0001.nc')
+ds = xr.open_dataset(filename, engine=FGMaxBackend, backend_kwargs={"epsg": epsg_code})
+ds.to_netcdf("fgmax0001.nc")
 
 # To see the use of clipping, change the tfinal in setrun to something like 2*3600.0
-# the fgmax0001_clipped.nc will only be the area where the wave arrived within the considered time. 
+# the fgmax0001_clipped.nc will only be the area where the wave arrived within the considered time.
 filename = "_output/fgmax0001.txt"
-ds = xr.open_dataset(filename, engine=FGMaxBackend, backend_kwargs={'epsg':epsg_code, "clip": True})
-ds.to_netcdf('fgmax0001_clipped.nc')
+ds = xr.open_dataset(
+    filename, engine=FGMaxBackend, backend_kwargs={"epsg": epsg_code, "clip": True}
+)
+ds.to_netcdf("fgmax0001_clipped.nc")
diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index d52736e8a..f2ab74c3f 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -105,18 +105,21 @@
 import os
 
 import numpy as np
+
 try:
     import rioxarray  # activate the rio accessor
     import xarray as xr
 except ImportError:
-    raise ImportError('rioxarray and xarray are required to use the FGOutBackend and FGMaxBackend')
+    raise ImportError(
+        "rioxarray and xarray are required to use the FGOutBackend and FGMaxBackend"
+    )
 
 from clawpack.geoclaw import fgmax_tools, fgout_tools
 from xarray.backends import BackendEntrypoint
 
-_qelements_dclaw = ["h", "hu", "hv", "hm", "pb", "hchi", "delta_a", "eta"]
-_qelements_geoclaw = ["h", "hu", "hv", "eta"]
-
+_qelements_dclaw = ["h", "hu", "hv", "hm", "pb", "hchi", "delta_a", "eta", "B"]
+_qelements_geoclaw = ["h", "hu", "hv", "eta", "B"]
+_qelements_geoclaw_bouss = ["h", "hu", "hv", "huc", "hvc", "eta", "B"]
 _qunits = {
     "h": "meters",
     "hu": "meters squared per second",
@@ -126,6 +129,9 @@
     "hchi": "meters",
     "delta_a": "meters",
     "eta": "meters",
+    "B": "meters",
+    "huc": "",
+    "hvc": "",
 }
 
 time_unit = "seconds"
@@ -154,6 +160,7 @@ class FGOutBackend(BackendEntrypoint):
     def open_dataset(
         self,
         filename,  # path to fgout file.
+        qmap="geoclaw",  # qmap value for FGoutGrid ('geoclaw', 'dclaw', or 'geoclaw-bouss')
         epsg=None,  # epsg code
         drop_variables=None,  # name of any elements of q to drop.
     ):
@@ -179,7 +186,9 @@ def open_dataset(
         else:
             raise ValueError("Invalid FGout output format. Must be ascii or binary.")
 
-        fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, output_format)
+        fgout_grid = fgout_tools.FGoutGrid(
+            fgno=fgno, outdir=outdir, output_format=output_format, qmap=qmap
+        )
 
         if fgout_grid.point_style != 2:
             raise ValueError("FGOutBackend only works with fg.point_style=2")
@@ -195,34 +204,51 @@ def open_dataset(
         ni = len(y)
 
         # mask based on dry tolerance
-        mask = fgout.q[0, :, :] < 0.001
-
-        # determine if geoclaw or dclaw
-        nvars = fgout.q.shape[0]
-        if nvars == 8:
+        try:
+            mask = fgout.h < 0.001
+        except:
+            try:
+                mask = (fgout.eta - fgout.B) < 0.001
+            except:
+                mask = np.ones((ni, nj), dtype=bool)
+
+        # determine if geoclaw, geoclaw-bouss, or dclaw
+        if qmap == "dclaw":
             _qelements = _qelements_dclaw
-        else:
+        elif qmap == "geoclaw":
             _qelements = _qelements_geoclaw
+        elif qmap == "geoclaw-bouss":
+            _qelements = _qelements_geoclaw_bouss
 
         # create data_vars dictionary
         data_vars = {}
-        for i in range(nvars):
-            Q = fgout.q[i, :, :]
-            # mask all but eta based on h presence.
-            if i < 7:
-                Q[mask] = nodata
-
-            # to keep xarray happy, need to transpose and flip ud.
-            Q = Q.T
-            Q = np.flipud(Q)
-            Q = Q.reshape((1, ni, nj))  # reshape to add a time dimension.
+        for i, i_var in enumerate(fgout_grid.q_out_vars):
 
             # construct variable
-            varname = _qelements[i]
 
-            data_array_attrs = {"units": _qunits[varname], "_FillValue": nodata}
+            # Find the varname in fgout.qmap associated with
+            # q_out_vars[i]
+            varname = None
+            for name, index in fgout_grid.qmap.items():
+                if i_var == index:
+                    varname = name
 
+            # construct xarray dataset if varname not in drop vars.
             if varname not in drop_variables:
+
+                Q = fgout.q[i, :, :]
+
+                # mask all but eta based on h presence.
+                if varname not in ("B", "eta"):
+                    Q[mask] = nodata
+
+                # to keep xarray happy, need to transpose and flip ud.
+                Q = Q.T
+                Q = np.flipud(Q)
+                Q = Q.reshape((1, ni, nj))  # reshape to add a time dimension.
+
+                data_array_attrs = {"units": _qunits[varname], "_FillValue": nodata}
+
                 data_vars[varname] = (
                     [
                         "time",
@@ -274,8 +300,7 @@ def open_dataset(
         filename,
         epsg=None,
         drop_variables=None,
-        clip=False, # if True, clip the entire array to the extent where arrival_time is defined.
-
+        clip=False,  # if True, clip the entire array to the extent where arrival_time is defined.
     ):
 
         if drop_variables is None:
@@ -285,7 +310,7 @@ def open_dataset(
         # e.g., output and 'fgmax_grids.data' in _output
         fgno = int(os.path.basename(filename).split(".")[0][-4:])
         outdir = os.path.dirname(filename)
-        data_file = os.path.join(outdir, 'fgmax_grids.data')
+        data_file = os.path.join(outdir, "fgmax_grids.data")
 
         fg = fgmax_tools.FGmaxGrid()
         fg.read_fgmax_grids_data(fgno=fgno, data_file=data_file)
@@ -472,7 +497,6 @@ def open_dataset(
                     "Time of minimum depth",
                 )
 
-
         # drop requested variables
         for var in drop_variables:
             if var in data_vars.keys():
@@ -505,7 +529,7 @@ def open_dataset(
 
         # clip
         if clip:
-            clip_data = fg.arrival_time.data>=0
+            clip_data = fg.arrival_time.data >= 0
             clip_data = clip_data.T
             clip_data = np.flipud(clip_data)
             ds = _clip(ds, clip_data)
@@ -524,8 +548,8 @@ def _clip(ds, clip_data):
 
     nrow, ncol = clip_data.shape
 
-    valid_col = np.sum(clip_data, axis=0)>0
-    valid_row = np.sum(clip_data, axis=1)>0
+    valid_col = np.sum(clip_data, axis=0) > 0
+    valid_row = np.sum(clip_data, axis=1) > 0
 
     sel_rows = np.nonzero(valid_row)[0]
     sel_cols = np.nonzero(valid_col)[0]

From eb3cd2019b9961a1b55af07d0a6d01b5d56d7e0b Mon Sep 17 00:00:00 2001
From: "Barnhart, Katherine (Katy) Ruth" <krbarnhart@usgs.gov>
Date: Tue, 20 Aug 2024 16:20:10 -0600
Subject: [PATCH 23/41] add a few more comments.

---
 src/python/geoclaw/xarray_backends.py | 39 ++++++++++++++++++++++-----
 1 file changed, 32 insertions(+), 7 deletions(-)

diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index f2ab74c3f..de9573e95 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -34,10 +34,20 @@
     # .qxxx file if ascii format is used.
     # the format, fg number, and frame number are inferred from the filename.
 
-    ds = xr.open_dataset(filename, engine=FGOutBackend, backend_kwargs={'epsg':epsg_code})
+    ds = xr.open_dataset(
+        filename,
+        engine=FGOutBackend,
+        backend_kwargs={
+            'qmap': 'geoclaw',
+            'epsg':epsg_code
+                })
     # ds is now an xarray object. It can be interacted with directly or written to netcdf using
     ds.write_netcdf('filename.nc')
 
+    # A 'qmap' backend_kwargs is required as it indicates the qmap used for
+    # selecting elements of q to include in fgout. See
+    # https://www.clawpack.org/dev/fgout.html#specifying-q-out-vars
+
     # Optionally, provide an epsg code to assign the associated coordinate system to the file.
     # default behavior assigns no coordinate system.
 
@@ -48,18 +58,32 @@
 
 Dimensions:
 
-Files opened with FGOutBackend will have dimensions time, y, x.
-Files opened with FGMaxBackend will have dimensions y, x.
+Files opened with FGOutBackend will have dimensions (time, y, x).
+Files opened with FGMaxBackend will have dimensions (y, x).
 
 Variable naming:
 
-For fixed grid geoclaw files, the dataset will have the following variables:
+For fixed grid files, the dataset will have the elements of q specified
+by qmap and q_out_vars. This may be as many as the full list described below or
+as few as specified by the user.
+
+For fixed grid geoclaw files, the full set of variables is:
+- h
+- hu
+- hv
+- eta
+- B
+
+For fixed grid geoclaw-bouss files, the full set of variables is:
 - h
 - hu
 - hv
+- huc
+- hvc
 - eta
+- B
 
-Fixed grid dclaw files will have
+Fixed grid dclaw files, the full set of variables is:
 - h
 - hu
 - hv
@@ -68,6 +92,7 @@
 - hchi
 - delta_a
 - eta
+- B
 
 Depending on the number of variables specified in the setrun.py fgmax files will
 have a portion of the following variables:
@@ -130,8 +155,8 @@
     "delta_a": "meters",
     "eta": "meters",
     "B": "meters",
-    "huc": "",
-    "hvc": "",
+    "huc": "meters squared per second",
+    "hvc": "meters squared per second",
 }
 
 time_unit = "seconds"

From 92375ff5e1dab3465eb9c86edc41a46e4151c776 Mon Sep 17 00:00:00 2001
From: Katy Barnhart <krbarnhart@usgs.gov>
Date: Wed, 21 Aug 2024 09:45:17 -0600
Subject: [PATCH 24/41] name updates

---
 src/python/geoclaw/xarray_backends.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index de9573e95..1f2c3887a 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -90,7 +90,7 @@
 - hm
 - pb
 - hchi
-- delta_a
+- bdiff
 - eta
 - B
 
@@ -142,7 +142,7 @@
 from clawpack.geoclaw import fgmax_tools, fgout_tools
 from xarray.backends import BackendEntrypoint
 
-_qelements_dclaw = ["h", "hu", "hv", "hm", "pb", "hchi", "delta_a", "eta", "B"]
+_qelements_dclaw = ["h", "hu", "hv", "hm", "pb", "hchi", "bdiff", "eta", "B"]
 _qelements_geoclaw = ["h", "hu", "hv", "eta", "B"]
 _qelements_geoclaw_bouss = ["h", "hu", "hv", "huc", "hvc", "eta", "B"]
 _qunits = {
@@ -152,11 +152,11 @@
     "hm": "meters",
     "pb": "newton per meter squared",
     "hchi": "meters",
-    "delta_a": "meters",
+    "bdiff": "meters",
     "eta": "meters",
     "B": "meters",
-    "huc": "meters squared per second",
-    "hvc": "meters squared per second",
+    "huc": "varies",
+    "hvc": "varies",
 }
 
 time_unit = "seconds"

From f69c70829ddc5701ff9689f79a9e36df311fc632 Mon Sep 17 00:00:00 2001
From: Katy Barnhart <krbarnhart@usgs.gov>
Date: Wed, 21 Aug 2024 09:47:28 -0600
Subject: [PATCH 25/41] bdiff->bdif

---
 src/python/geoclaw/xarray_backends.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index 1f2c3887a..4603b8290 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -90,7 +90,7 @@
 - hm
 - pb
 - hchi
-- bdiff
+- bdif
 - eta
 - B
 
@@ -142,7 +142,7 @@
 from clawpack.geoclaw import fgmax_tools, fgout_tools
 from xarray.backends import BackendEntrypoint
 
-_qelements_dclaw = ["h", "hu", "hv", "hm", "pb", "hchi", "bdiff", "eta", "B"]
+_qelements_dclaw = ["h", "hu", "hv", "hm", "pb", "hchi", "bdif", "eta", "B"]
 _qelements_geoclaw = ["h", "hu", "hv", "eta", "B"]
 _qelements_geoclaw_bouss = ["h", "hu", "hv", "huc", "hvc", "eta", "B"]
 _qunits = {
@@ -152,7 +152,7 @@
     "hm": "meters",
     "pb": "newton per meter squared",
     "hchi": "meters",
-    "bdiff": "meters",
+    "bdif": "meters",
     "eta": "meters",
     "B": "meters",
     "huc": "varies",

From 145233c0c9f58fa039f9c31f24b27d592436e26b Mon Sep 17 00:00:00 2001
From: Katy Barnhart <krbarnhart@usgs.gov>
Date: Wed, 21 Aug 2024 10:07:15 -0600
Subject: [PATCH 26/41] remove vestigial _qelements

---
 src/python/geoclaw/xarray_backends.py | 12 +-----------
 1 file changed, 1 insertion(+), 11 deletions(-)

diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index 4603b8290..cedb24266 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -142,9 +142,7 @@
 from clawpack.geoclaw import fgmax_tools, fgout_tools
 from xarray.backends import BackendEntrypoint
 
-_qelements_dclaw = ["h", "hu", "hv", "hm", "pb", "hchi", "bdif", "eta", "B"]
-_qelements_geoclaw = ["h", "hu", "hv", "eta", "B"]
-_qelements_geoclaw_bouss = ["h", "hu", "hv", "huc", "hvc", "eta", "B"]
+
 _qunits = {
     "h": "meters",
     "hu": "meters squared per second",
@@ -237,14 +235,6 @@ def open_dataset(
             except:
                 mask = np.ones((ni, nj), dtype=bool)
 
-        # determine if geoclaw, geoclaw-bouss, or dclaw
-        if qmap == "dclaw":
-            _qelements = _qelements_dclaw
-        elif qmap == "geoclaw":
-            _qelements = _qelements_geoclaw
-        elif qmap == "geoclaw-bouss":
-            _qelements = _qelements_geoclaw_bouss
-
         # create data_vars dictionary
         data_vars = {}
         for i, i_var in enumerate(fgout_grid.q_out_vars):

From 66d0d636bebc220d07801f554ecff6963b51e50a Mon Sep 17 00:00:00 2001
From: Katy Barnhart <krbarnhart@usgs.gov>
Date: Wed, 21 Aug 2024 11:34:48 -0600
Subject: [PATCH 27/41] Add example that should work with/without dask

---
 .../use_xarray_backends.py                    | 33 ++++++++++++++-----
 1 file changed, 25 insertions(+), 8 deletions(-)

diff --git a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
index b0cf65d5a..b489e3377 100644
--- a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
+++ b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
@@ -3,7 +3,7 @@
 try:
     import rioxarray
     import xarray as xr
-except:
+except ImportError:
     "You must install xarray and rioxarray in order to use the xarray backends"
     raise
 
@@ -30,16 +30,33 @@
 ds.to_netcdf("fgout0001_0001.nc")
 
 # It is possible to combine all fgout files into a single netcdf file
-# using xr.open_mfdataset
+# using xr.open_mfdataset (requires dask) or xr.concat (does not require dask)
 # https://docs.xarray.dev/en/stable/generated/xarray.open_mfdataset.html
+# https://docs.xarray.dev/en/latest/generated/xarray.concat.html
+# for instructions on installing xarray with dask see:
+# https://docs.xarray.dev/en/latest/getting-started-guide/installing.html#instructions
+
 fgout_files = glob.glob("_output/fgout0001.b*")
-ds_all = xr.open_mfdataset(
-    fgout_files,
-    engine=FGOutBackend,
-    backend_kwargs={"epsg": epsg_code, "qmap": "geoclaw"},
-)
 
-print(ds_all)
+try:
+    ds_all = xr.open_mfdataset(
+        fgout_files,
+        engine=FGOutBackend,
+        backend_kwargs={"epsg": epsg_code, "qmap": "geoclaw"},
+    )
+except ImportError:  # if dask is not available, use xr.concat.
+    fgouts = []
+    for filename in fgout_files:
+        ds = xr.open_dataset(
+            filename,
+            engine=FGOutBackend,
+            backend_kwargs={"epsg": epsg_code, "qmap": "geoclaw"},
+        )
+        fgouts.append(ds)
+
+    ds_all = xr.concat(fgouts, dim="time")
+
+# save out.
 ds_all.to_netcdf("fgout_all.nc")
 
 # An example of a fgmax grid.

From eb7981e2d130c1df1cb1783f86a6e7fce2528781 Mon Sep 17 00:00:00 2001
From: Katy Barnhart <krbarnhart@usgs.gov>
Date: Wed, 21 Aug 2024 11:35:26 -0600
Subject: [PATCH 28/41] improve drytol handling

(now tested for both absence of h, and absence of h, eta, and B)
---
 src/python/geoclaw/fgout_tools.py     | 210 +++++++++++++-------------
 src/python/geoclaw/xarray_backends.py |  10 +-
 2 files changed, 108 insertions(+), 112 deletions(-)

diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py
index 7fdc47e19..de21b51a8 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -18,7 +18,7 @@
             fgout frames, as a single netCDF file
 - function read_netcdf: Read a netCDF file and return a list of fgout frames,
             assuming the file contains all the qoi's needed to reconstruct q.
-- function read_netcdf_arrays: Read a netCDF file and extract the 
+- function read_netcdf_arrays: Read a netCDF file and extract the
             requested quantities of interest as numpy arrays.
 - print_netcdf_info: Print info about the contents of a netCDF file containing
             some fgout frames.
@@ -56,7 +56,7 @@ def __init__(self, fgout_grid, frameno=None):
 
     # Define shortcuts to attributes of self.fgout_grid that are the same
     # for all frames (e.g. X,Y) to avoid storing grid for every frame.
-    
+
     @property
     def x(self):
         return self.fgout_grid.x
@@ -64,15 +64,15 @@ def x(self):
     @property
     def y(self):
         return self.fgout_grid.y
-        
+
     @property
     def X(self):
         return self.fgout_grid.X
-        
+
     @property
     def Y(self):
         return self.fgout_grid.Y
-        
+
     @property
     def delta(self):
         return self.fgout_grid.delta
@@ -80,19 +80,19 @@ def delta(self):
     @property
     def extent_centers(self):
         return self.fgout_grid.extent_centers
-        
+
     @property
     def extent_edges(self):
         return self.fgout_grid.extent_edges
-        
+
     @property
     def drytol(self):
         return self.fgout_grid.drytol
-        
+
     @property
     def qmap(self):
         return self.fgout_grid.qmap
-            
+
     # Define attributes such as h as @properties with lazy evaluation:
     # the corresponding array is created and stored only when first
     # accessed by the user.  Those not needed are not created.
@@ -105,14 +105,14 @@ def h(self):
             try:
                 i_h = q_out_vars.index(self.qmap['h'])
                 self._h = self.q[i_h,:,:]
-            except:
+            except ValueError: # if q_out_vars does not contain qmap['h'], a value error is thrown
                 try:
                     i_eta = q_out_vars.index(self.qmap['eta'])
                     i_B = q_out_vars.index(self.qmap['B'])
                     self._h = self.q[i_eta,:,:] - self.q[i_B,:,:]
-                except:
+                except ValueError:
                     print('*** Could not find h or eta-B in q_out_vars')
-                    raise
+                    raise AttributeError
         return self._h
 
     @property
@@ -322,7 +322,7 @@ def __init__(self,fgno=None,outdir='.',output_format=None,
                          'bdif':7, 'eta':8, 'B':9}
         else:
             raise ValueError('Invalid qmap: %s' % qmap)
-        
+
         # GeoClaw input values:
         # Many of these should be read from fgout_grids.data
         # using read_fgout_grids_data before using read_frame
@@ -344,12 +344,12 @@ def __init__(self,fgno=None,outdir='.',output_format=None,
         self.q_out_vars = [1,2,3,4]  # list of which components to print
                                           # default: h,hu,hv,eta (5=topo)
         self.nqout = None # number of vars to print out
-        
+
         self.drytol = 1e-3          # used for computing u,v from hu,hv
 
         # private attributes for those that are only created if
         # needed by the user:
-        
+
         self._X = None
         self._Y = None
         self._x = None
@@ -375,7 +375,7 @@ def y(self):
             dy = (self.y2 - self.y1)/self.ny
             self._y = numpy.linspace(self.y1+dy/2, self.y2-dy/2, self.ny)
         return self._y
-        
+
     @property
     def X(self):
         """2D array X of longitudes (cell centers)"""
@@ -397,7 +397,7 @@ def delta(self):
             dy = (self.y2 - self.y1)/self.ny
             self._delta = (dx,dy)
         return self._delta
-        
+
     @property
     def extent_centers(self):
         """Extent of cell centers [xmin,xmax,ymin,ymax]"""
@@ -405,7 +405,7 @@ def extent_centers(self):
             self._extent_centers = [self.x.min(), self.x.max(),\
                                     self.y.min(), self.y.max()]
         return self._extent_centers
-        
+
     @property
     def extent_edges(self):
         """Extent of cell edges [xmin,xmax,ymin,ymax]"""
@@ -414,10 +414,10 @@ def extent_edges(self):
             self._extent_edges = [self.x.min()-dx/2, self.x.max()+dx/2,\
                                     self.y.min()-dy/2, self.y.max()+dy/2]
         return self._extent_edges
-        
+
 
     # Create plotdata of class clawpack.visclaw.ClawPlotData
-    # only when needed for reading GeoClaw output, 
+    # only when needed for reading GeoClaw output,
     # since this must be done after fgno, outdir, output_format
     # have been specified:
 
@@ -450,7 +450,7 @@ def set_plotdata(self):
             file_prefix_str = plotdata.file_prefix + '.b'
         else:
             file_prefix_str = plotdata.file_prefix + '.q'
-        
+
         # Contrary to usual default, set save_frames to False so that all
         # the fgout frames read in are not saved in plotdata.framesoln_dict.
         # Otherwise this might be a memory hog when making an animation with
@@ -472,7 +472,7 @@ def read_fgout_grids_data(self, fgno=None, data_file='fgout_grids.data'):
 
         data_path = os.path.join(self.outdir, data_file)
         print('Reading fgout grid info from \n    %s' % data_path)
-        
+
         with open(data_path) as filep:
             lines = filep.readlines()
         fgout_input = None
@@ -577,7 +577,7 @@ def read_fgout_grids_data_pre511(self, fgno=None,
 
         data_path = os.path.join(self.outdir, data_file)
         print('Reading fgout grid info from \n    %s' % data_path)
-        
+
         with open(data_path) as filep:
             lines = filep.readlines()
         fgout_input = None
@@ -628,7 +628,7 @@ def write_to_fgout_data(self, fid):
             errmsg = 'point_style %s is not implement, only point_style==2' \
                                     % point_style
             raise NotImplementedError(errmsg)
-                
+
         if self.output_format == 'ascii':
             output_format = 1
         elif self.output_format == 'binary32':
@@ -640,7 +640,7 @@ def write_to_fgout_data(self, fid):
             raise NotImplementedError(errmsg)
 
 
-        
+
         # write header, independent of point_style:
         #fid = open(self.input_file_name,'w')
         fid.write("\n")
@@ -648,7 +648,7 @@ def write_to_fgout_data(self, fid):
 
         fid.write("%i                           # output_style\n" \
                     % self.output_style)
-        
+
         if self.output_style == 1:
             assert self.tstart is not None, 'Need to set tstart'
             assert self.tend is not None, 'Need to set tend'
@@ -659,7 +659,7 @@ def write_to_fgout_data(self, fid):
         elif self.output_style == 2:
             self.nout = len(self.output_times)
             fid.write("%i %s           # nout\n" % (self.nout, 11*" "))
-            
+
             # remove [] and , from list of times:
             output_times_str = repr(list(self.output_times))[1:-1]
             output_times_str = output_times_str.replace(',','')
@@ -695,7 +695,7 @@ def write_to_fgout_data(self, fid):
             print("   lower left  = (%15.10f,%15.10f)" % (x1,y1))
             print("   upper right = (%15.10f,%15.10f)" % (x2,y2))
             print("   dx = %15.10e,  dy = %15.10e" % (dx,dy))
-            
+
         # q_out_vars is a list of q components to print, e.g. [1,4]
 
         format = len(self.q_out_vars) * '%s '
@@ -725,7 +725,7 @@ def read_frame(self, frameno):
             self.read_fgout_grids_data()
             fgoutX = self.X
             fgoutY = self.Y
-            
+
         try:
             fr = self.plotdata.getframe(frameno)
         except:
@@ -741,11 +741,11 @@ def read_frame(self, frameno):
         fgout_frame.q = state.q
 
         fgout_frame.t = state.t
-    
+
         fgout_frame.frameno = frameno
-        
-        X,Y = patch.grid.p_centers[:2]    
-            
+
+        X,Y = patch.grid.p_centers[:2]
+
         if not numpy.allclose(X, fgoutX):
             errmsg = '*** X read from output does not match fgout_grid.X'
             raise ValueError(errmsg)
@@ -753,7 +753,7 @@ def read_frame(self, frameno):
         if not numpy.allclose(Y, fgoutY):
             errmsg = '*** Y read from output does not match fgout_grid.Y'
             raise ValueError(errmsg)
-                
+
         return fgout_frame
 
 
@@ -766,29 +766,29 @@ def make_fgout_fcn_xy(fgout, qoi, method='nearest',
     """
     Create a function that can be called at (x,y) and return the qoi
     interpolated in space from the fgout array.
-    
-    qoi should be a string (e.g. 'h', 'u' or 'v') corresponding to 
+
+    qoi should be a string (e.g. 'h', 'u' or 'v') corresponding to
     an attribute of fgout.
-    
-    The function returned takes arguments x,y that can be floats or 
+
+    The function returned takes arguments x,y that can be floats or
     (equal length) 1D arrays of values that lie within the spatial
-    extent of fgout. 
-    
-    bounds_error and fill_value determine the behavior if (x,y) is not in 
+    extent of fgout.
+
+    bounds_error and fill_value determine the behavior if (x,y) is not in
     the bounds of the data, as in scipy.interpolate.RegularGridInterpolator.
     """
-    
+
     from scipy.interpolate import RegularGridInterpolator
-    
+
     try:
         q = getattr(fgout,qoi)
     except:
         print('*** fgout missing attribute qoi = %s?' % qoi)
-        
+
     err_msg = '*** q must have same shape as fgout.X\n' \
             + 'fgout.X.shape = %s,   q.shape = %s' % (fgout.X.shape,q.shape)
     assert fgout.X.shape == q.shape, err_msg
-    
+
     x1 = fgout.X[:,0]
     y1 = fgout.Y[0,:]
     fgout_fcn1 = RegularGridInterpolator((x1,y1), q, method=method,
@@ -816,39 +816,39 @@ def make_fgout_fcn_xyt(fgout1, fgout2, qoi, method_xy='nearest',
     """
     Create a function that can be called at (x,y,t) and return the qoi
     interpolated in space and time between the two frames fgout1 and fgout2.
-    
-    qoi should be a string (e.g. 'h', 'u' or 'v') corresponding to 
+
+    qoi should be a string (e.g. 'h', 'u' or 'v') corresponding to
     an attribute of fgout.
-    
+
     method_xy is the method used in creating the spatial interpolator,
     and is passed to make_fgout_fcn_xy.
-    
+
     method_t is the method used for interpolation in time, currently only
     'linear' is supported, which linearly interpolates.
-    
-    bounds_error and fill_value determine the behavior if (x,y,t) is not in 
+
+    bounds_error and fill_value determine the behavior if (x,y,t) is not in
     the bounds of the data.
-    
+
     The function returned takes arguments x,y (floats or equal-length 1D arrays)
     of values that lie within the spatial extent of fgout1, fgout2
     (which are assumed to cover the same uniform grid at different times)
-    and t should be a float that lies between fgout1.t and fgout2.t. 
+    and t should be a float that lies between fgout1.t and fgout2.t.
     """
-    
+
     assert numpy.allclose(fgout1.X, fgout2.X), \
                             '*** fgout1 and fgout2 must have same X'
     assert numpy.allclose(fgout1.Y, fgout2.Y), \
                             '*** fgout1 and fgout2 must have same Y'
-    
+
     t1 = fgout1.t
     t2 = fgout2.t
     #assert t1 < t2, '*** expected fgout1.t < fgout2.t'
-    
+
     fgout1_fcn_xy = make_fgout_fcn_xy(fgout1, qoi, method=method_xy,
                        bounds_error=bounds_error, fill_value=fill_value)
     fgout2_fcn_xy = make_fgout_fcn_xy(fgout2, qoi, method=method_xy,
                        bounds_error=bounds_error, fill_value=fill_value)
-                
+
     def fgout_fcn(x,y,t):
         """
         Function that can be evaluated at single point or arrays (x,y)
@@ -877,15 +877,15 @@ def fgout_fcn(x,y,t):
                 qout = (1-alpha)*qout1 + alpha*qout2
         else:
             raise NotImplementedError('method_t = %s not supported' % method_t)
-            
+
         return qout
-        
+
     return fgout_fcn
 
 # ===============================
 # Functions for writing a set of fgout frames as a netCDF file, and
-# reading such a file:    
-    
+# reading such a file:
+
 def write_netcdf(fgout_frames, fname_nc='fgout_frames.nc',
                  qois = ['h','hu','hv','eta'], datatype='f4',
                  include_B0=False, include_Bfinal=False,
@@ -894,63 +894,63 @@ def write_netcdf(fgout_frames, fname_nc='fgout_frames.nc',
     Write a list of fgout frames (at different times on the same rectangular
     grid) to a single netCDF file, with some metadata and the topography,
     if desired.
-    
+
     fgout_frames should be a list of FGoutFrame objects, all of the same size
     and at increasing times.
-    
+
     fname_nc is the name of the file to write.
-    
+
     qois is a list of strings, the quantities of interest to include in file.
-        This could include any of: 
+        This could include any of:
 
             'h', 'eta', 'hu', 'hv', 'u', 'v', 's', 'hss', 'B'.
 
         All other quantities can be computed from h, hu, hv, eta,
         the original fgout variables from GeoClaw, but for some applications
         you might only want to save 'h' and 's', for example.
-    
+
     datatype should be 'f4' [default] or 'f8', specifying bytes per qoi value.
         'f8' has full precision of the original data, but the file will be
         twice as large and may not be needed for downstream applications.
-        
+
     Note that the topography B = eta - h, so it is not necessary to store all
     three of these.  Also, B is often the same for all frames, so rather than
     storing B at each frame as a qoi, two other options are also provided
     (and then storing eta or h for all frames allows calculating the other):
-    
+
     include_Bfinal: If True, include the topography B array from the final frame
     as the Bfinal array.
-    
+
     include_B0: If True, include the topography B array from the first frame
-    as the B0 array.  This is only useful if, e.g., the first frame is initial 
+    as the B0 array.  This is only useful if, e.g., the first frame is initial
     topography before co-seismic deformation, and at later times the topography
-    is always equal to Bfinal. 
-    
+    is always equal to Bfinal.
+
     `description` is a string that will be added as metadata.
-    A metadata field `history` will also be added, which includes the 
-    time the file was created and the path to the directory where it was made. 
-    """    
-    
+    A metadata field `history` will also be added, which includes the
+    time the file was created and the path to the directory where it was made.
+    """
+
     import netCDF4
     from datetime import datetime, timedelta
     from cftime import num2date, date2num
     import time
     timestr = time.ctime(time.time())  # current time for metadata
-    
+
     fg_times = numpy.array([fg.t for fg in fgout_frames])
-    
+
     if verbose:
         print('Creating %s with fgout frames at times: ' % fname_nc)
         print(fg_times)
-    
+
     fg0 = fgout_frames[0]
     x = fg0.x
     y = fg0.y
-    
+
     xs = numpy.array([fg.x for fg in fgout_frames])
     ys = numpy.array([fg.y for fg in fgout_frames])
     # assert same for all times
-    
+
     units = {'h':'meters', 'eta':'meters', 'hu':'m^2/s', 'hv':'m^2/s',
              'u':'m/s', 'v':'m/s', 's':'m/s', 'hss':'m^3/s^2', 'B':'meters'}
 
@@ -969,12 +969,12 @@ def write_netcdf(fgout_frames, fname_nc='fgout_frames.nc',
         latitudes = rootgrp.createVariable('lat','f8',('lat',))
         latitudes[:] = y
         latitudes.units = 'degrees_north'
-        
+
         time = rootgrp.createDimension('time', len(fg_times))
         times = rootgrp.createVariable('time','f8',('time',))
         times[:] = fg_times
         times.units = 'seconds'
-        
+
         if 0:
             # Could make times be datetimes relative to some event time, e.g.:
             times.units = 'seconds since 1700-01-26 21:00:00.0'
@@ -982,7 +982,7 @@ def write_netcdf(fgout_frames, fname_nc='fgout_frames.nc',
             dates = [datetime(1700,1,26,21) + timedelta(seconds=ss) \
                         for ss in fg_times]
             times[:] = date2num(dates,units=times.units,calendar=times.calendar)
-        
+
         if include_B0:
             B0 = rootgrp.createVariable('B0',datatype,('lon','lat',))
             B0[:,:] = fg0.B
@@ -993,7 +993,7 @@ def write_netcdf(fgout_frames, fname_nc='fgout_frames.nc',
             Bfinal = rootgrp.createVariable('Bfinal',datatype,('lon','lat',))
             Bfinal[:,:] = fg_final.B
             Bfinal.units = 'meters'
-        
+
         for qoi in qois:
             qoi_frames = [getattr(fgout,qoi) for fgout in fgout_frames]
             qoi_var = rootgrp.createVariable(qoi,datatype,('time','lon','lat',))
@@ -1002,7 +1002,7 @@ def write_netcdf(fgout_frames, fname_nc='fgout_frames.nc',
 
 def get_as_array(var, rootgrp, verbose=True):
     """
-    Utility function to retrieve variable from netCDF file and convert to 
+    Utility function to retrieve variable from netCDF file and convert to
     numpy array.
     """
     a = rootgrp.variables.get(var, None)
@@ -1011,7 +1011,7 @@ def get_as_array(var, rootgrp, verbose=True):
         return numpy.array(a)
     else:
         if verbose: print('    Did not find %s' % var)
-        return None    
+        return None
 
 
 def read_netcdf_arrays(fname_nc, qois, verbose=True):
@@ -1022,19 +1022,19 @@ def read_netcdf_arrays(fname_nc, qois, verbose=True):
         'h', 'eta', 'hu', 'hv', 'u', 'v', 's', 'hss', 'B'.
 
     qois can also include the time-independent 'B0' and/or 'Bfinal'
-    
-    Returns 
+
+    Returns
 
         x, y, t, qoi_arrays
 
     where x,y define the longitude, latitudes, t is the times of the frames,
     and qoi_arrays is a dictionary indexed by the strings from qois.
-    
+
     Each dict element is an array with shape (len(t), len(x), len(y))
     for time-dependent qoi's, or (len(x), len(y)) for B0 or Bfinal,
     or None if that qoi was not found in the netCDF file.
     """
-     
+
     import netCDF4
 
     with netCDF4.Dataset(fname_nc, 'r') as rootgrp:
@@ -1046,17 +1046,17 @@ def read_netcdf_arrays(fname_nc, qois, verbose=True):
         x = get_as_array('lon', rootgrp, verbose)
         y = get_as_array('lat', rootgrp, verbose)
         t = get_as_array('time', rootgrp, verbose)
-        
+
         vars = list(rootgrp.variables)
-        
+
         qoi_arrays = {}
         for qoi in qois:
             qoi_array = get_as_array(qoi, rootgrp, verbose)
             qoi_arrays[qoi] = qoi_array
-            
+
     return x,y,t,qoi_arrays
-            
-            
+
+
 
 def read_netcdf(fname_nc, fgout_grid=None, verbose=True):
     """
@@ -1065,7 +1065,7 @@ def read_netcdf(fname_nc, fgout_grid=None, verbose=True):
     the qoi's 'h','hu','hv','eta' required to reconstruct the q array
     as output by GeoClaw.
     """
-    
+
     import netCDF4
 
     with netCDF4.Dataset(fname_nc, 'r') as rootgrp:
@@ -1079,7 +1079,7 @@ def read_netcdf(fname_nc, fgout_grid=None, verbose=True):
         t = get_as_array('time', rootgrp, verbose)
 
         vars = list(rootgrp.variables)
-        
+
 
         for qoi in ['h','hu','hv','eta']:
             errmsg = '*** Cannot reconstruct fgout frame without %s' % qoi
@@ -1088,14 +1088,14 @@ def read_netcdf(fname_nc, fgout_grid=None, verbose=True):
         hu = get_as_array('hu', rootgrp, verbose)
         hv = get_as_array('hv', rootgrp, verbose)
         eta = get_as_array('eta', rootgrp, verbose)
-        
+
         if (x is None) or (y is None):
             print('*** Could not create grid')
         else:
             X,Y = numpy.meshgrid(x,y)
-        
+
         fgout_frames = []
-                
+
         for k in range(eta.shape[0]):
             fgout = FGoutFrame(fgout_grid=fgout_grid, frameno=k)
             fgout.x = x
@@ -1109,9 +1109,9 @@ def read_netcdf(fname_nc, fgout_grid=None, verbose=True):
             fgout.X = X
             fgout.Y = Y
             fgout_frames.append(fgout)
-            
+
         print('Created fgout_frames as list of length %i' % len(fgout_frames))
-        
+
     return fgout_frames
 
 def print_netcdf_info(fname_nc):
@@ -1127,7 +1127,7 @@ def print_netcdf_info(fname_nc):
         t = get_as_array('time', rootgrp, verbose=False)
 
         vars = list(rootgrp.variables)
-        
+
         print('===================================================')
         print('netCDF file %s contains:' % fname_nc)
         print('description: \n', rootgrp.description)
@@ -1137,5 +1137,3 @@ def print_netcdf_info(fname_nc):
         print('%i times from %.3f to %.3f' % (len(t),t[0],t[-1]))
         print('variables: ',vars)
         print('===================================================')
-    
-        
diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index cedb24266..1e26e60f5 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -142,7 +142,6 @@
 from clawpack.geoclaw import fgmax_tools, fgout_tools
 from xarray.backends import BackendEntrypoint
 
-
 _qunits = {
     "h": "meters",
     "hu": "meters squared per second",
@@ -229,11 +228,10 @@ def open_dataset(
         # mask based on dry tolerance
         try:
             mask = fgout.h < 0.001
-        except:
-            try:
-                mask = (fgout.eta - fgout.B) < 0.001
-            except:
-                mask = np.ones((ni, nj), dtype=bool)
+            # internally fgout_tools will try fgou.eta-fgout.B if h is not present.
+        except AttributeError:
+            print("FGOutBackend: No h, eta, or B. No mask applied.")
+            mask = np.ones((nj, ni), dtype=bool)
 
         # create data_vars dictionary
         data_vars = {}

From d4079f87ffb94cbac22d8a7902752b7fcd0e65dc Mon Sep 17 00:00:00 2001
From: Katy Barnhart <krbarnhart@usgs.gov>
Date: Wed, 21 Aug 2024 13:09:14 -0600
Subject: [PATCH 29/41] don't specify error type.

---
 examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
index b489e3377..091f55f4f 100644
--- a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
+++ b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
@@ -3,7 +3,7 @@
 try:
     import rioxarray
     import xarray as xr
-except ImportError:
+except:
     "You must install xarray and rioxarray in order to use the xarray backends"
     raise
 

From 61460c1a1708e2e73b998596619cd177231704fb Mon Sep 17 00:00:00 2001
From: Katy Barnhart <krbarnhart@usgs.gov>
Date: Wed, 21 Aug 2024 13:12:04 -0600
Subject: [PATCH 30/41] switch from hard coded drytol to user-specified drytol

---
 src/python/geoclaw/xarray_backends.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index 1e26e60f5..3738b4e9f 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -184,6 +184,8 @@ def open_dataset(
         filename,  # path to fgout file.
         qmap="geoclaw",  # qmap value for FGoutGrid ('geoclaw', 'dclaw', or 'geoclaw-bouss')
         epsg=None,  # epsg code
+        drytol=0.001, # dry tolerance for masking elements of q that are not eta or B
+            # used only if h or eta-B is available based on q_out_vars.
         drop_variables=None,  # name of any elements of q to drop.
     ):
 
@@ -227,7 +229,7 @@ def open_dataset(
 
         # mask based on dry tolerance
         try:
-            mask = fgout.h < 0.001
+            mask = fgout.h < drytol
             # internally fgout_tools will try fgou.eta-fgout.B if h is not present.
         except AttributeError:
             print("FGOutBackend: No h, eta, or B. No mask applied.")

From f61d15525fef9be23d6926ac89790f3ca7525d83 Mon Sep 17 00:00:00 2001
From: Katy Barnhart <krbarnhart@usgs.gov>
Date: Wed, 21 Aug 2024 13:24:19 -0600
Subject: [PATCH 31/41] update use of dry_tolerance

---
 .../use_xarray_backends.py                    | 14 ++++++++++-
 src/python/geoclaw/xarray_backends.py         | 23 +++++++++++--------
 2 files changed, 27 insertions(+), 10 deletions(-)

diff --git a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
index 091f55f4f..7145195c2 100644
--- a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
+++ b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
@@ -24,7 +24,19 @@
 # the format, fg number, and frame number are inferred from the filename.
 
 ds = xr.open_dataset(
-    filename, engine=FGOutBackend, backend_kwargs={"epsg": epsg_code, "qmap": "geoclaw"}
+    filename,
+    engine=FGOutBackend,
+    backend_kwargs={
+        "epsg": epsg_code,
+        "qmap": "geoclaw",
+        # qmap is the qmap specified to the fgout object in setrun.py see
+        # the following documentation page for more details. 
+        # https://www.clawpack.org/dev/fgout.html#specifying-q-out-vars
+        "dry_tolerance": None,
+        # variables that are not eta and B are masked
+        # where h>dry_tolerance. To turn this functionality
+        # off set dry_tolerance = None.
+    },
 )
 # ds is now an xarray object. It can be interacted with directly or written to netcdf using
 ds.to_netcdf("fgout0001_0001.nc")
diff --git a/src/python/geoclaw/xarray_backends.py b/src/python/geoclaw/xarray_backends.py
index 3738b4e9f..326d0d908 100644
--- a/src/python/geoclaw/xarray_backends.py
+++ b/src/python/geoclaw/xarray_backends.py
@@ -184,8 +184,12 @@ def open_dataset(
         filename,  # path to fgout file.
         qmap="geoclaw",  # qmap value for FGoutGrid ('geoclaw', 'dclaw', or 'geoclaw-bouss')
         epsg=None,  # epsg code
-        drytol=0.001, # dry tolerance for masking elements of q that are not eta or B
-            # used only if h or eta-B is available based on q_out_vars.
+        dry_tolerance=0.001,
+        # dry tolerance used for for masking elements of q that are not
+        # eta or B. Default behavior is to mask all elements of q except for
+        # eta and B where h>0.001.
+        # used only if h or eta-B is available based on q_out_vars.
+        # if dry_tolerance = None, no masking is applied to any variable.
         drop_variables=None,  # name of any elements of q to drop.
     ):
 
@@ -228,12 +232,13 @@ def open_dataset(
         ni = len(y)
 
         # mask based on dry tolerance
-        try:
-            mask = fgout.h < drytol
-            # internally fgout_tools will try fgou.eta-fgout.B if h is not present.
-        except AttributeError:
-            print("FGOutBackend: No h, eta, or B. No mask applied.")
-            mask = np.ones((nj, ni), dtype=bool)
+        if dry_tolerance is not None:
+            try:
+                mask = fgout.h < dry_tolerance
+                # internally fgout_tools will try fgou.eta-fgout.B if h is not present.
+            except AttributeError:
+                print("FGOutBackend: No h, eta, or B. No mask applied.")
+                mask = np.ones((nj, ni), dtype=bool)
 
         # create data_vars dictionary
         data_vars = {}
@@ -254,7 +259,7 @@ def open_dataset(
                 Q = fgout.q[i, :, :]
 
                 # mask all but eta based on h presence.
-                if varname not in ("B", "eta"):
+                if (varname not in ("B", "eta")) and dry_tolerance is not None:
                     Q[mask] = nodata
 
                 # to keep xarray happy, need to transpose and flip ud.

From 987417a461a23bbc579b9a9e19b224ccfe6e7741 Mon Sep 17 00:00:00 2001
From: Katy Barnhart <krbarnhart@usgs.gov>
Date: Wed, 21 Aug 2024 13:51:50 -0600
Subject: [PATCH 32/41] fix exceptions.

---
 .../tsunami/chile2010_fgmax-fgout/use_xarray_backends.py  | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
index 7145195c2..6d7fa9128 100644
--- a/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
+++ b/examples/tsunami/chile2010_fgmax-fgout/use_xarray_backends.py
@@ -3,7 +3,7 @@
 try:
     import rioxarray
     import xarray as xr
-except:
+except ImportError:
     "You must install xarray and rioxarray in order to use the xarray backends"
     raise
 
@@ -30,7 +30,7 @@
         "epsg": epsg_code,
         "qmap": "geoclaw",
         # qmap is the qmap specified to the fgout object in setrun.py see
-        # the following documentation page for more details. 
+        # the following documentation page for more details.
         # https://www.clawpack.org/dev/fgout.html#specifying-q-out-vars
         "dry_tolerance": None,
         # variables that are not eta and B are masked
@@ -56,7 +56,9 @@
         engine=FGOutBackend,
         backend_kwargs={"epsg": epsg_code, "qmap": "geoclaw"},
     )
-except ImportError:  # if dask is not available, use xr.concat.
+except ValueError:  # if dask is not available, use xr.concat.
+# if dask is not installed xr.open_mfdataset() will fail with something like
+# ValueError: unrecognized chunk manager dask - must be one of: []
     fgouts = []
     for filename in fgout_files:
         ds = xr.open_dataset(

From afc7ff99d1b598aa10a5b1f0be08c711dac4c0b0 Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Sun, 25 Aug 2024 09:31:48 -0700
Subject: [PATCH 33/41] update examples/tsunami/chile2010/setplot.py

Improve this by using new options added in recent releases so the
plots look better and this can better serve as a template for other
problems.
---
 examples/tsunami/chile2010/setplot.py | 67 +++++++++++++--------------
 1 file changed, 33 insertions(+), 34 deletions(-)

diff --git a/examples/tsunami/chile2010/setplot.py b/examples/tsunami/chile2010/setplot.py
index 525132748..fc01ec976 100644
--- a/examples/tsunami/chile2010/setplot.py
+++ b/examples/tsunami/chile2010/setplot.py
@@ -57,18 +57,14 @@ def addgauges(current_data):
 
     # Set up for axes in this figure:
     plotaxes = plotfigure.new_plotaxes('pcolor')
-    plotaxes.title = 'Surface'
-    plotaxes.scaled = True
-
-    def fixup(current_data):
-        import pylab
-        addgauges(current_data)
-        t = current_data.t
-        t = t / 3600.  # hours
-        pylab.title('Surface at %4.2f hours' % t, fontsize=20)
-        pylab.xticks(fontsize=15)
-        pylab.yticks(fontsize=15)
-    plotaxes.afteraxes = fixup
+    plotaxes.title = 'Surface at time h:m:s'
+    plotaxes.aspect_latitude = -30.
+    plotaxes.xticks_fontsize = 10
+    plotaxes.yticks_fontsize = 10
+    plotaxes.xlabel = 'longitude'
+    plotaxes.ylabel = 'latitude'
+
+    plotaxes.afteraxes = addgauges
 
     # Water
     plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
@@ -78,6 +74,8 @@ def fixup(current_data):
     plotitem.pcolor_cmin = -0.2
     plotitem.pcolor_cmax = 0.2
     plotitem.add_colorbar = True
+    plotitem.colorbar_shrink = 0.8
+    plotitem.colorbar_extend = 'both'
     plotitem.amr_celledges_show = [0,0,0]
     plotitem.patchedges_show = 1
 
@@ -114,18 +112,37 @@ def fixup(current_data):
 
     # Set up for axes in this figure:
     plotaxes = plotfigure.new_plotaxes()
-    plotaxes.xlimits = 'auto'
-    plotaxes.ylimits = 'auto'
+    plotaxes.time_scale = 1/3600.  # convert to hours
+    plotaxes.xlimits = [0, 9]
+    plotaxes.ylimits = [-0.3, 0.3]
     plotaxes.title = 'Surface'
+    plotaxes.title_fontsize = 15
+    plotaxes.time_label = 'time (hours)'
+    plotaxes.ylabel = 'surface elevation (m)'
+    plotaxes.grid = True
+
+    def add_obs(current_data):
+        from pylab import plot, legend
+        gaugeno = current_data.gaugeno
+        if gaugeno == 32412:
+            try:
+                tgauge = TG32412[:,0] / 3600.  # convert to hours
+                plot(tgauge, TG32412[:,1], 'r', linewidth=1.0)
+                legend(['GeoClaw','Obs'],loc='lower right')
+            except: pass
+
+    plotaxes.afteraxes = add_obs
 
     # Plot surface as blue curve:
     plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
-    plotitem.plot_var = 3
+    plotitem.plot_var = 3  # eta
     plotitem.plotstyle = 'b-'
+    plotitem.kwargs = {'linewidth':1.8}
+
 
     # Plot topo as green curve:
     plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
-    plotitem.show = False
+    plotitem.show = False  # not being used
 
     def gaugetopo(current_data):
         q = current_data.q
@@ -137,24 +154,6 @@ def gaugetopo(current_data):
     plotitem.plot_var = gaugetopo
     plotitem.plotstyle = 'g-'
 
-    def add_zeroline(current_data):
-        from pylab import plot, legend, xticks, floor, axis, xlabel
-        t = current_data.t 
-        gaugeno = current_data.gaugeno
-
-        if gaugeno == 32412:
-            try:
-                plot(TG32412[:,0], TG32412[:,1], 'r')
-                legend(['GeoClaw','Obs'],loc='lower right')
-            except: pass
-            axis((0,t.max(),-0.3,0.3))
-
-        plot(t, 0*t, 'k')
-        n = int(floor(t.max()/3600.) + 2)
-        xticks([3600*i for i in range(n)], ['%i' % i for i in range(n)])
-        xlabel('time (hours)')
-
-    plotaxes.afteraxes = add_zeroline
 
 
     #-----------------------------------------

From 3ede3d37307f407fe4b50633f93b389c6304a2e8 Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Sun, 25 Aug 2024 15:32:30 -0700
Subject: [PATCH 34/41] Fix src/bouss/amr2.f90 call to set_fgout, adding nvar

This was modified in the shallow/amr2.f90 code to support more fgout
flexibility in #624, and this change is also required in the bouss version.

Also enhanced the examples/bouss README files a bit.
---
 examples/bouss/README.txt             | 6 ++++++
 examples/bouss/radial_flat/README.rst | 9 +++++++++
 src/2d/bouss/amr2.f90                 | 4 ++--
 3 files changed, 17 insertions(+), 2 deletions(-)

diff --git a/examples/bouss/README.txt b/examples/bouss/README.txt
index 047848de2..2e71aa599 100644
--- a/examples/bouss/README.txt
+++ b/examples/bouss/README.txt
@@ -1,4 +1,6 @@
 
+GeoClaw Boussinesq solver examples
+
 The radial_flat subdirectory contains one example using the Boussinesq
 solver introduced in Clawpack v5.10.0.
 
@@ -14,8 +16,12 @@ OpenMP along with MPI.  Some flags have to be set as environment variables
 or directly in the application Makefile, e.g. see the lines commented out in
 radial_flat/Makefile.
 
+The file setenv.sh illustrates how you might set some environment
+variables for the bash shell.
+
 A file petscMPIoptions is also required to set some PETSc parameters for the
 iterative method used to solve the large sparse linear systems that arise at 
 each refinement level when the Boussinesq equations are solved. 
 One of the environment variables mentioned above points to this file, and a
 sample is included in this directory.
+
diff --git a/examples/bouss/radial_flat/README.rst b/examples/bouss/radial_flat/README.rst
index 1d76fb818..6c3e6e086 100644
--- a/examples/bouss/radial_flat/README.rst
+++ b/examples/bouss/radial_flat/README.rst
@@ -8,6 +8,15 @@ A Gaussian hump of water at the origin spreads out over a flat bottom,
 as specified by the topo file `flat100.txt` (uniform water depth 100 m).
 The equations are solved in Cartesian coordinates with the SGN equations.
 
+Running the GeoClaw Boussinesq solvers requires PETSc and MPI.
+For more details see the documentation
+  https://www.clawpack.org/bouss2d.html
+and
+  $CLAW/geoclaw/examples/bouss/README.txt
+Run
+  make check
+in this directory to check if things seem ok for running this code.
+
 A flagregion specified by `RuledRectangle_Diagonal.data` (that is created by
 code in `setrun.py` is used to allow flagging for refinement to
 level 2 only near the diagonal (for abs(x-y) < 1000).  The code is set up to
diff --git a/src/2d/bouss/amr2.f90 b/src/2d/bouss/amr2.f90
index e8fa9d984..09e2b241b 100644
--- a/src/2d/bouss/amr2.f90
+++ b/src/2d/bouss/amr2.f90
@@ -513,7 +513,7 @@ program amr2
         call read_dtopo_settings()        ! specifies file with dtopo from earthquake
         call read_topo_settings(rest)     ! specifies topography (bathymetry) files
         call set_qinit()                  ! specifies file with dh if this used instead
-        call set_fgout(rest)            ! Fixed grid settings
+        call set_fgout(rest,nvar)         ! Fixed grid settings
         call setup_variable_friction()    ! Variable friction parameter
         call set_storm()                  ! Set storm parameters
         call set_regions()                ! Set refinement regions
@@ -548,7 +548,7 @@ program amr2
         call read_topo_settings(rest)     ! specifies topography (bathymetry) files
         call set_qinit()                  ! specifies file with dh if this used instead
 
-        call set_fgout(rest)              ! Fixed grid settings
+        call set_fgout(rest,nvar)         ! Fixed grid settings
         call setup_variable_friction()    ! Variable friction parameter
         call set_multilayer()             ! Set multilayer SWE parameters
         call set_storm()                  ! Set storm parameters

From ada9a0b49dc53307e4e21ecf081ad4fa7935684b Mon Sep 17 00:00:00 2001
From: Kyle Mandli <kyle.mandli@gmail.com>
Date: Wed, 4 Sep 2024 17:25:10 -0400
Subject: [PATCH 35/41] Add new time title format

---
 src/python/geoclaw/surge/plot.py | 34 ++++++++++++++++++++++++++------
 1 file changed, 28 insertions(+), 6 deletions(-)

diff --git a/src/python/geoclaw/surge/plot.py b/src/python/geoclaw/surge/plot.py
index b67ecf546..3832a06b6 100644
--- a/src/python/geoclaw/surge/plot.py
+++ b/src/python/geoclaw/surge/plot.py
@@ -122,13 +122,35 @@ def plot_landfall_gauge(gauge, axes, landfall=0.0, style='b', kwargs={}):
 # ========================================================================
 #  Surge related helper functions
 # ========================================================================
-def days_figure_title(current_data, land_fall=0.0):
-    t = (current_data.t - land_fall) / (60**2 * 24)
-    days = int(t)
-    hours = (t - int(t)) * 24.0
+def days_figure_title(cd, land_fall=0.0, new_time=False):
+    r"""Helper function that puts the time relative to landfall in title
 
-    title = current_data.plotaxes.title
-    plt.title('%s at day %3i, hour %2.1f' % (title, days, hours))
+    New version of title is available if *new_time = True*
+    """
+    if new_time:
+        if cd.t < land_fall:
+            sign = "-"
+        else:
+            sign = " "
+        minutes, seconds = divmod(abs(np.round(cd.t - land_fall)), 60)
+        hours, minutes = divmod(minutes, 60)
+        days, hours = divmod(hours, 24)
+        days = int(days)
+        hours = int(hours)
+        minutes = int(minutes)
+        if cd.t < 0:
+            sign = "-"
+        else:
+            sign = " "
+        title = cd.plotaxes.title
+        plt.title(f'{title} at t = {sign}{days:d}, {hours:02d}:{minutes:02d}')
+    else:
+        t = (cd.t - land_fall) / (60**2 * 24)
+        days = int(t)
+        hours = (t - int(t)) * 24.0
+
+        title = cd.plotaxes.title
+        plt.title('%s at day %3i, hour %2.1f' % (title, days, hours))
 
 
 def surge_afteraxes(current_data, track, land_fall=0.0, plot_direction=False,

From 5bf7bbd3f726357bdd3e719ad3ad7e5529c90518 Mon Sep 17 00:00:00 2001
From: Kyle Mandli <kyle.mandli@gmail.com>
Date: Wed, 4 Sep 2024 17:25:34 -0400
Subject: [PATCH 36/41] Fix storm track plotting

---
 src/python/geoclaw/surge/storm.py | 27 ++++++++++++++-------------
 1 file changed, 14 insertions(+), 13 deletions(-)

diff --git a/src/python/geoclaw/surge/storm.py b/src/python/geoclaw/surge/storm.py
index 4fe670ed1..d99dea171 100644
--- a/src/python/geoclaw/surge/storm.py
+++ b/src/python/geoclaw/surge/storm.py
@@ -1164,9 +1164,15 @@ def write_tcvitals(self, path, verbose=False):
     # ================
     #  Track Plotting
     # ================
-    def plot(self, ax, radius=None, t_range=None, coordinate_system=2, track_style='ko--',
-             categorization="NHC", fill_alpha=0.25, fill_color='red'):
-        """TO DO:  Write doc-string"""
+    def plot(self, ax, t_range=None, coordinate_system=2, 
+                       track_style='ko--', categorization="NHC", 
+                       radius=None, fill_alpha=0.25, fill_color='red'):
+        """TO DO:  Write doc-string
+
+        Quick notes:
+         - *radius = None* will not plot a swath
+         - *track_style = {}* will plot categories for the track
+         """
 
         import matplotlib.pyplot as plt
 
@@ -1198,17 +1204,15 @@ def plot(self, ax, radius=None, t_range=None, coordinate_system=2, track_style='
             colors = [track_style.get(category, cat_color_defaults[category])
                       for category in self.category(categorization=categorization)]
             for i in range(t.shape[0] - 1):
-                ax.plot(x[i:i+2], y[i:i+2], color=colors[i], marker="o")
+                ax.plot(x[i:i+2], y[i:i+2], color=colors[i])
 
         else:
             raise ValueError("The `track_style` should be a string or dict.")
 
         # Plot swath
-        if (isinstance(radius, float) or isinstance(radius, np.ndarray)
-                or radius is None):
-
-            if radius is None:
-                # Default behavior
+        if radius is not None:
+            # Any string as of right now will trigger this
+            if isinstance(radius, str):
                 if self.storm_radius is None:
                     raise ValueError("Cannot use storm radius for plotting "
                                      "the swath as the data is not available.")
@@ -1265,7 +1269,6 @@ def plot(self, ax, radius=None, t_range=None, coordinate_system=2, track_style='
 
     # =========================================================================
     # Other Useful Routines
-
     def category(self, categorization="NHC", cat_names=False):
         r"""Categorizes storm based on relevant storm data
 
@@ -1318,9 +1321,7 @@ def category(self, categorization="NHC", cat_names=False):
                        12: "Hurricane"}
 
         elif categorization.upper() == "NHC":
-            # TODO:  Change these to m/s (knots are how these are defined).
-            # Definitely not in the correct format now
-            # TODO:  Add TD and TS designations
+            # NHC uses knots
             speeds = units.convert(self.max_wind_speed, "m/s", "knots")
             category = (np.zeros(speeds.shape) +
                         (speeds < 30) * -1 +

From 90a91a91e1af011f0f1b85615eec604208a08b3d Mon Sep 17 00:00:00 2001
From: Kyle Mandli <kyle.mandli@gmail.com>
Date: Fri, 20 Sep 2024 09:29:46 -0400
Subject: [PATCH 37/41] Fix spelling error

---
 src/2d/shallow/surge/storm_module.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/2d/shallow/surge/storm_module.f90 b/src/2d/shallow/surge/storm_module.f90
index 632cbff90..4c0afe0a4 100644
--- a/src/2d/shallow/surge/storm_module.f90
+++ b/src/2d/shallow/surge/storm_module.f90
@@ -162,7 +162,7 @@ subroutine set_storm(data_file)
                 case(2)
                     rotation => S_rotation
                 case default
-                    stop " *** ERROR *** Roation override invalid."
+                    stop " *** ERROR *** Rotation override invalid."
             end select
             read(unit,*)
 

From 5c65e221073b989fbf11f6c9baf544d27a1c3af3 Mon Sep 17 00:00:00 2001
From: Kyle Mandli <kyle.mandli@gmail.com>
Date: Fri, 20 Sep 2024 09:31:05 -0400
Subject: [PATCH 38/41] Add default bounds to land plotitem

---
 src/python/geoclaw/surge/plot.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/src/python/geoclaw/surge/plot.py b/src/python/geoclaw/surge/plot.py
index 3832a06b6..851f7773a 100644
--- a/src/python/geoclaw/surge/plot.py
+++ b/src/python/geoclaw/surge/plot.py
@@ -428,7 +428,7 @@ def add_pressure(plotaxes, bounds=None, plot_type='pcolor', shrink=1.0):
         pass
 
 
-def add_land(plotaxes, plot_type='pcolor', bounds=None):
+def add_land(plotaxes, plot_type='pcolor', bounds=[0, 50]):
     """Add plotitem for land"""
 
     if plot_type == 'pcolor':
@@ -436,9 +436,8 @@ def add_land(plotaxes, plot_type='pcolor', bounds=None):
         plotitem.show = True
         plotitem.plot_var = geoplot.land
         plotitem.pcolor_cmap = land_cmap
-        if bounds is not None:
-            plotitem.pcolor_cmin = bounds[0]
-            plotitem.pcolor_cmax = bounds[1]
+        plotitem.pcolor_cmin = bounds[0]
+        plotitem.pcolor_cmax = bounds[1]
         plotitem.add_colorbar = False
         plotitem.amr_celledges_show = [0] * 10
         plotitem.amr_patchedges_show = [1, 1, 1, 1, 1, 0, 0]

From 6c967cc3627e04bb3f3b4b17f40064a2ac74c42c Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Mon, 21 Oct 2024 16:27:31 -0700
Subject: [PATCH 39/41] Properly set fgout q_out_vars to default

Should be `[h,hu,hv,eta]` as in old versions, but now `qmap['eta']` varies
depending on equations being solved (e.g. `geoclaw-bouss` or `dclaw`).
---
 src/python/geoclaw/fgout_tools.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py
index de21b51a8..f6fc56c52 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -341,9 +341,12 @@ def __init__(self,fgno=None,outdir='.',output_format=None,
         # Note output_format will be reset by read_fgout_grids_data:
         self.output_format = output_format
 
-        self.q_out_vars = [1,2,3,4]  # list of which components to print
-                                          # default: h,hu,hv,eta (5=topo)
-        self.nqout = None # number of vars to print out
+        # list of which components to print:
+        try:
+            # default: h,hu,hv,eta as in previous versions of GeoClaw
+            self.q_out_vars = [1,2,3,self.qmap['eta']]
+        except:
+            self.q_out_vars = [1,2,3]  # if qmap['eta'] not defined
 
         self.drytol = 1e-3          # used for computing u,v from hu,hv
 

From 12d57bc2d52ea03c231412b41f438d30e8996df4 Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Tue, 29 Oct 2024 13:47:33 -0700
Subject: [PATCH 40/41] check if 'eta' in qmap.keys() instead of try-except

---
 src/python/geoclaw/fgout_tools.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py
index f6fc56c52..6b957a9d8 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -342,10 +342,10 @@ def __init__(self,fgno=None,outdir='.',output_format=None,
         self.output_format = output_format
 
         # list of which components to print:
-        try:
+        if 'eta' in self.qmap.keys():
             # default: h,hu,hv,eta as in previous versions of GeoClaw
             self.q_out_vars = [1,2,3,self.qmap['eta']]
-        except:
+        else:
             self.q_out_vars = [1,2,3]  # if qmap['eta'] not defined
 
         self.drytol = 1e-3          # used for computing u,v from hu,hv

From 6f49189ff4fa62eea27eed1bc50a335fa111de31 Mon Sep 17 00:00:00 2001
From: Randy LeVeque <rjl@uw.edu>
Date: Tue, 29 Oct 2024 13:48:05 -0700
Subject: [PATCH 41/41] fix chile2010_fgmax-fgout/plot_fgout.py

---
 examples/tsunami/chile2010_fgmax-fgout/plot_fgout.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/examples/tsunami/chile2010_fgmax-fgout/plot_fgout.py b/examples/tsunami/chile2010_fgmax-fgout/plot_fgout.py
index 8b61666d3..a24f5a68d 100644
--- a/examples/tsunami/chile2010_fgmax-fgout/plot_fgout.py
+++ b/examples/tsunami/chile2010_fgmax-fgout/plot_fgout.py
@@ -13,6 +13,7 @@
 
 # Instantiate object for reading fgout frames:
 fgout_grid = fgout_tools.FGoutGrid(fgno, outdir, output_format)
+fgout_grid.read_fgout_grids_data()
 
 # Plot one frame of fgout data 
 fgframe = 20