diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml
index 4d1d669d4..99779c0b4 100644
--- a/.github/workflows/testing.yml
+++ b/.github/workflows/testing.yml
@@ -40,7 +40,12 @@ jobs:
     - name: Setup geoclaw
       run: |
         cd geoclaw
-        git checkout ${{ github.ref }}
+        if [[ "${{ github.event_name }}" == "pull_request" ]]; then
+           git fetch origin pull/${{ github.event.pull_request.number }}/merge:PR
+           git checkout PR
+         else
+           git checkout ${{ github.ref_name }}
+         fi
     - name: Lint with flake8
       run: |
         cd geoclaw
diff --git a/.travis.yml b/.travis.yml
deleted file mode 100644
index 9b7540911..000000000
--- a/.travis.yml
+++ /dev/null
@@ -1,48 +0,0 @@
-language: python
-python:
-  - 3.8
-env:
-  - BUILD_TYPE="Release"
-
-before_install:
-  - sudo apt-get update
-  - sudo apt-get install gfortran liblapack-pic liblapack-dev libnetcdf-dev libnetcdff-dev
-  - pip install netCDF4
-  - pip install xarray
-  - pip install scipy
-  - pip install matplotlib
-  - git clone --branch=master --depth=100 --quiet git://github.com/clawpack/clawpack
-  - cd clawpack
-  - git submodule init
-  - git submodule update
-  - rm -rf geoclaw
-  - ln -s ../ geoclaw
-  # Print versions being used
-  - python -c "import numpy; print(numpy.__version__)"
-  - python -c "import netCDF4; print(netCDF4.__version__)"
-
-install:
-  - export PYTHONPATH=${PWD}:$PYTHONPATH
-  - export CLAW=${PWD}
-  - export FC=gfortran
-  - export NETCDF4_DIR=/usr
-
-before_script:
-  # Print CPU info for debugging purposes:
-  - cat /proc/cpuinfo
-  - gfortran -v
-  - /usr/bin/f95 -v
-  - echo "PYTHONPATH="$PYTHONPATH
-  - echo "CLAW="$CLAW
-  - echo "FC="$FC
-  - echo "NETCDF4_DIR="$NETCDF4_DIR
-
-script:
-  - cd $CLAW/geoclaw
-  - nosetests -sv
-
-after_failure:
- - for failed_test_path in *_output ; do cat $failed_test_path/run_output.txt ; cat $failed_test_path/error_output.txt ; done
-
-notifications:
-  email: false
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
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/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.
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..7fdc47e19 100644
--- a/src/python/geoclaw/fgout_tools.py
+++ b/src/python/geoclaw/fgout_tools.py
@@ -52,6 +52,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.
@@ -87,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
@@ -96,16 +101,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 +141,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 +163,43 @@ 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,:,:]
+                #print('+++qmap["eta"] = %i' % self.qmap["eta"])
+                #print('+++i_eta = %i' % 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,:,:]
+                #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
         return self._B
 
     @property
@@ -161,6 +216,86 @@ 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"""
+        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
+
+    @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):
 
@@ -169,10 +304,28 @@ 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,
+                 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':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}
+        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
         self.id = ''  # identifier, optional
         self.point_style = 2  # only option currently supported
         self.npts = None
@@ -184,8 +337,13 @@ def __init__(self,fgno=None,outdir=None,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.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,7 +470,10 @@ 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'
 
-        with open(data_file) as filep:
+        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):
@@ -324,10 +485,17 @@ 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])
+        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:
@@ -354,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:
@@ -371,7 +544,75 @@ 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
+        print('Found fgout grid q_out_vars = ',self.q_out_vars)
+        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.
+
+        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.
+        """
 
+        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):
         """
@@ -455,7 +696,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 +710,22 @@ def read_frame(self, frameno):
 
         from datetime import timedelta
 
+        try:
+            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
+            
         try:
             fr = self.plotdata.getframe(frameno)
         except:
@@ -483,36 +744,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'
diff --git a/src/python/geoclaw/surge/plot.py b/src/python/geoclaw/surge/plot.py
index c2ada663f..b67ecf546 100644
--- a/src/python/geoclaw/surge/plot.py
+++ b/src/python/geoclaw/surge/plot.py
@@ -13,6 +13,9 @@
 
 from __future__ import absolute_import
 from __future__ import print_function
+
+import warnings
+
 import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib.colors as colors
@@ -23,6 +26,7 @@
 import clawpack.visclaw.gaugetools as gaugetools
 import clawpack.visclaw.geoplot as geoplot
 import clawpack.geoclaw.data as geodata
+# import clawpack.geoclaw.surge.storm
 
 # TODO:  Assign these based on data files
 bathy_index = 0
@@ -41,7 +45,9 @@
 
 surge_data = geodata.SurgeData()
 
-
+# ==============================
+#  Track Plotting Functionality
+# ==============================
 class track_data(object):
     """Read in storm track data from run output"""
 
@@ -68,8 +74,8 @@ def get_track(self, frame):
 
             # Check to make sure that this fixed the problem
             if self._data.shape[0] < frame + 1:
-                print(" *** WARNING *** Could not find track data for ",
-                      "frame %s." % frame)
+                warnings.warn(f" *** WARNING *** Could not find track data",
+                               " for frame {frame}.")
                 return None, None, None
 
         return self._data[frame, 1:]
@@ -165,8 +171,15 @@ def pressure(cd):
     # The division by 100.0 is to convert from Pa to millibars
     return cd.aux[pressure_field, :, :] / 100.0
 
-# def category(Storm, cd):
-#     return cd.aux[Storm.category, :, :]
+
+def storm_radius(cd, track):
+    """Distance from center of storm"""
+    track_data = track.get_track(cd.frameno)
+
+    if track_data[0] is not None and track_data[1] is not None:
+        return np.sqrt((cd.x - track_data[0])**2 + (cd.y - track_data[1])**2)
+    else:
+        return None
 
 
 # ========================================================================
@@ -435,6 +448,15 @@ def add_bathy_contours(plotaxes, contour_levels=None, color='k'):
     plotitem.patchedges_show = 0
 
 
+def add_storm_radii(plotaxes, track, radii=[100e3], color='r'):
+    """Add radii to plots based on storm position"""
+    plotitem = plotaxes.new_plotitem(name="storm radius", 
+                                     plot_type="2d_contour")
+    plotitem.plot_var = lambda cd: storm_radius(cd, track)
+    plotitem.contour_levels = radii
+    plotitem.contour_colors = color
+
+
 # ===== Storm related plotting =======
 def sec2days(seconds):
     """Converst seconds to days."""
diff --git a/src/python/geoclaw/surge/storm.py b/src/python/geoclaw/surge/storm.py
index 2fe89971a..479385b4d 100644
--- a/src/python/geoclaw/surge/storm.py
+++ b/src/python/geoclaw/surge/storm.py
@@ -314,20 +314,27 @@ def read_geoclaw(self, path, verbose=False):
         # Read header
         with open(path, "r") as data_file:
             num_casts = int(data_file.readline())
-            self.time_offset = datetime.datetime.strptime(
-                data_file.readline()[:19], "%Y-%m-%dT%H:%M:%S"
-            )
+            time = data_file.readline()[:19]
+            try:
+                self.time_offset = datetime.datetime.strptime(time, "%Y-%m-%dT%H:%M:%S")
+            except ValueError:
+                self.time_offset = float(time)
         # Read rest of data
         data = np.loadtxt(path, skiprows=3)
         num_forecasts = data.shape[0]
-        self.eye_location = np.empty((2, num_forecasts))
+        self.eye_location = np.empty((num_forecasts, 2))
         assert num_casts == num_forecasts
-        self.t = [
-            self.time_offset + datetime.timedelta(seconds=data[i, 0])
-            for i in range(num_forecasts)
-        ]
-        self.eye_location[0, :] = data[:, 1]
-        self.eye_location[1, :] = data[:, 2]
+        if isinstance(self.time_offset, datetime.datetime):
+            self.t = np.array(
+                [
+                    self.time_offset + datetime.timedelta(seconds=data[i, 0])
+                    for i in range(num_forecasts)
+                ]
+            )
+        else:
+            self.t = data[:, 0]
+        self.eye_location[:, 0] = data[:, 1]
+        self.eye_location[:, 1] = data[:, 2]
         self.max_wind_speed = data[:, 3]
         self.max_wind_radius = data[:, 4]
         self.central_pressure = data[:, 5]
@@ -1386,7 +1393,7 @@ def write_atcf(self, path, verbose=False):
                         "".join(
                             (
                                 ", " * 2,
-                                "%s" % seconds2date(self.t[n]),
+                                "%s" % self.seconds2date(self.t[n]),
                                 ", " * 4,
                                 "%s" % (int(self.eye_location[n, 0] * 10.0)),
                                 ", ",
@@ -1546,6 +1553,144 @@ 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"""
+
+        import matplotlib.pyplot as plt
+
+        # Extract information for plotting the track/swath
+        t = self.t
+        x = self.eye_location[:, 0]
+        y = self.eye_location[:, 1]
+        if t_range is not None:
+            t = np.ma.masked_outside(t, t_range[0], t_range[1])
+            x = np.ma.array(x, mask=t.mask).compressed()
+            y = np.ma.array(y, mask=t.mask).compressed()
+            t = t.compressed()
+
+        # Plot track
+        if isinstance(track_style, str):
+            # Plot the track as a simple line with the given style
+            ax.plot(x, y, track_style)
+        elif isinstance(track_style, dict):
+            if self.max_wind_speed is None:
+                raise ValueError(
+                    "Maximum wind speed not available so "
+                    "plotting catgories is not available."
+                )
+
+            # Plot the track using the colors provided in the dictionary
+            cat_color_defaults = {
+                5: "red",
+                4: "yellow",
+                3: "orange",
+                2: "green",
+                1: "blue",
+                0: "gray",
+                -1: "lightgray",
+            }
+            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")
+
+        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 self.storm_radius is None:
+                    raise ValueError(
+                        "Cannot use storm radius for plotting "
+                        "the swath as the data is not available."
+                    )
+                else:
+                    if coordinate_system == 1:
+                        _radius = self.storm_radius
+                    elif coordinate_system == 2:
+                        _radius = units.convert(self.storm_radius, "m", "lat-long")
+                    else:
+                        raise ValueError(
+                            f"Unknown coordinate system "
+                            f"{coordinate_system} provided."
+                        )
+
+            elif isinstance(radius, float):
+                # Only one value for the radius was given, replicate
+                _radius = np.ones(self.t.shape) * radius
+            elif isinstance(radius, np.ndarray):
+                # The array passed is the array to use
+                _radius = radius
+            else:
+                raise ValueError(
+                    "Invalid input argument for radius.  Should " "be a float or None"
+                )
+
+            # Draw first and last points
+            ax.add_patch(plt.Circle((x[0], y[0]), _radius[0], color=fill_color))
+            if t.shape[0] > 1:
+                ax.add_patch(plt.Circle((x[-1], y[-1]), _radius[-1], color=fill_color))
+
+            # Draw path around inner points
+            if t.shape[0] > 2:
+                for i in range(t.shape[0] - 1):
+                    p = np.array([(x[i], y[i]), (x[i + 1], y[i + 1])])
+                    v = p[1] - p[0]
+                    if abs(v[1]) > 1e-16:
+                        n = np.array([1, -v[0] / v[1]], dtype=float)
+                    elif abs(v[0]) > 1e-16:
+                        n = np.array([-v[1] / v[0], 1], dtype=float)
+                    else:
+                        raise Exception("Zero-vector given")
+                    n /= np.linalg.norm(n)
+                    n *= _radius[i]
+
+                    ax.fill(
+                        (
+                            p[0, 0] + n[0],
+                            p[0, 0] - n[0],
+                            p[1, 0] - n[0],
+                            p[1, 0] + n[0],
+                        ),
+                        (
+                            p[0, 1] + n[1],
+                            p[0, 1] - n[1],
+                            p[1, 1] - n[1],
+                            p[1, 1] + n[1],
+                        ),
+                        facecolor=fill_color,
+                        alpha=fill_alpha,
+                    )
+                    ax.add_patch(
+                        plt.Circle(
+                            (p[1][0], p[1, 1]),
+                            _radius[i],
+                            color=fill_color,
+                            alpha=fill_alpha,
+                        )
+                    )
+
     # =========================================================================
     # Other Useful Routines
 
@@ -1817,62 +1962,6 @@ def available_models():
     return output
 
 
-def make_multi_structure(path):
-    r"""Create a dictionary of Storm objects for ATCF files with multiple storm tracks in them"""
-    with open(path, "r") as f:
-        lines = f.readlines()
-        curTime = "test"
-        curTrack = "test"
-        os.mkdir("Clipped_ATCFs")
-        stormDict = {}
-        for line in lines:
-            lineArr = line.split(", ")
-            if curTime in lineArr[2]:
-                if curTrack in lineArr[4]:
-                    fileWrite.writelines(line)
-                else:
-                    fileWrite.close()
-                    stormDict[curTime].update(
-                        {
-                            curTrack: Storm(
-                                path=os.path.join(
-                                    os.path.expandvars(os.getcwd()),
-                                    "Clipped_ATCFs",
-                                    curTime,
-                                    curTrack,
-                                ),
-                                file_format="ATCF",
-                            )
-                        }
-                    )
-                    curTrack = lineArr[4]
-                    fileWrite = open("Clipped_ATCFs/" + curTime + "/" + curTrack, "w")
-                    fileWrite.writelines(line)
-            else:
-                if curTime != "test":
-                    fileWrite.close()
-                    stormDict[curTime].update(
-                        {
-                            curTrack: Storm(
-                                path=os.path.join(
-                                    os.path.expandvars(os.getcwd()),
-                                    "Clipped_ATCFs",
-                                    curTime,
-                                    curTrack,
-                                ),
-                                file_format="ATCF",
-                            )
-                        }
-                    )
-                curTime = lineArr[2]
-                curTrack = lineArr[4]
-                stormDict[curTime] = {}
-                os.mkdir("Clipped_ATCFs/" + curTime)
-                fileWrite = open("Clipped_ATCFs/" + curTime + "/" + curTrack, "w")
-                fileWrite.writelines(line)
-    return stormDict
-
-
 def _set_engine_kwargs(path):
     if isinstance(path, FSMap) or Path(path).suffix == ".zarr":
         return "zarr", {"consolidated": True}
diff --git a/src/python/geoclaw/topotools.py b/src/python/geoclaw/topotools.py
index 1743e1ae6..db6b6410c 100644
--- a/src/python/geoclaw/topotools.py
+++ b/src/python/geoclaw/topotools.py
@@ -414,7 +414,7 @@ def delta(self):
 
 
     def __init__(self, path=None, topo_type=None, topo_func=None, 
-                       unstructured=False):
+                       unstructured=False, **kwargs):
         r"""Topography initialization routine.
         
         See :class:`Topography` for more info.
@@ -444,7 +444,7 @@ def __init__(self, path=None, topo_type=None, topo_func=None,
 
         if path:
             self.read(path=path, topo_type=topo_type,
-                      unstructured=unstructured)
+                      unstructured=unstructured, **kwargs)
 
     def set_xyZ(self, X, Y, Z):
         r"""