Index: /anuga_work/production/australia_ph2/melbourne/build_elevation.py
===================================================================
--- /anuga_work/production/australia_ph2/melbourne/build_elevation.py	(revision 6692)
+++ /anuga_work/production/australia_ph2/melbourne/build_elevation.py	(revision 6692)
@@ -0,0 +1,96 @@
+"""Build the elevation data to run a tsunami inundation scenario 
+for Busselton, WA, Australia.
+
+Input: elevation data from project.py
+Output: pts file stored in project.topographies_dir 
+The run_model.py is reliant on the output of this script.
+
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+import os
+from os.path import join
+
+# Related major packages
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
+from anuga.shallow_water.data_manager import dem2pts
+from anuga.shallow_water.data_manager import start_screen_catcher
+from anuga.shallow_water.data_manager import copy_code_files
+from anuga.geospatial_data.geospatial_data import Geospatial_data
+
+# Application specific imports
+from setup_model import project   # Definition of file names and polygons
+
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+copy_code_files(project.output_build,__file__, 
+                os.path.dirname(project.__file__)+os.sep+\
+                project.__name__+'.py' )
+start_screen_catcher(project.output_build)
+
+
+#------------------------------------------------------------------------------
+# Preparation of topographic data
+#
+# Convert ASC 2 DEM 2 PTS using source data and store result in source data
+# Do for coarse and fine data
+# Fine pts file to be clipped to area of interest
+#------------------------------------------------------------------------------
+
+print 'project.bounding_polygon', project.bounding_polygon
+print 'project.combined_elevation_basename', project.combined_elevation_basename
+
+# Create Geospatial data from ASCII files
+geospatial_data = {}
+if not project.ascii_grid_filenames == []:
+    for filename in project.ascii_grid_filenames:
+        absolute_filename = join(project.topographies_folder, filename)
+        convert_dem_from_ascii2netcdf(absolute_filename,
+                                      basename_out=absolute_filename,
+                                      use_cache=True,
+                                      verbose=True)
+        dem2pts(absolute_filename, use_cache=True, verbose=True)
+
+        geospatial_data[filename] = Geospatial_data(file_name=absolute_filename+'.pts',
+                                                    verbose=True)
+
+# Create Geospatial data from TXT files
+if not project.point_filenames == []:
+    for filename in project.point_filenames:
+        absolute_filename = join(project.topographies_folder, filename)
+        geospatial_data[filename] = Geospatial_data(file_name=absolute_filename,
+                                                    verbose=True)
+
+#-------------------------------------------------------------------------------
+# Combine, clip and export dataset 
+#-------------------------------------------------------------------------------
+
+print 'Add geospatial objects' # except', project.offshore_name5
+G = None
+for key in geospatial_data:
+    #if key == project.point_filenames[4] or key == project.ascii_filenames[1]:
+    #    # Skip these files for now
+    #    continue
+    
+    G += geospatial_data[key]
+        
+print 'Clip combined geospatial data'
+##G_clip = G.clip_outside(project.poly_aoi1)
+##G_all = G_clip + geospatial_data[project.point_filenames[4]] 
+#G_clipped = G_all.clip(project.poly_all)
+G_clip = G.clip(project.bounding_polygon)
+
+
+print 'Export combined DEM file'
+G_clip.export_points_file(project.combined_elevation + '.pts')
+print 'Do txt version too'
+# Use for comparision in ARC
+G_clip.export_points_file(project.combined_elevation + '.txt')
+
Index: /anuga_work/production/australia_ph2/melbourne/build_urs_boundary.py
===================================================================
--- /anuga_work/production/australia_ph2/melbourne/build_urs_boundary.py	(revision 6692)
+++ /anuga_work/production/australia_ph2/melbourne/build_urs_boundary.py	(revision 6692)
@@ -0,0 +1,166 @@
+"""
+Replacement function that wraps the legacy function urs2sts()
+and passes in the data from the EventSelection <event>.list file.
+"""
+
+import os
+import os.path
+from time import localtime, strftime, gmtime
+
+import project
+
+from Scientific.IO.NetCDF import NetCDFFile
+import Numeric as num
+
+from anuga.shallow_water.data_manager import urs2sts
+
+import Numeric as num
+
+
+#-------------------------------------------------------------------------------
+# Get gauges (timeseries of index points)
+#-------------------------------------------------------------------------------
+def get_sts_gauge_data(filename, verbose=False):
+    print 'get_sts_gauge_data: filename=%s' % filename
+    fid = NetCDFFile(filename+'.sts', 'r')      #Open existing file for read
+    permutation = fid.variables['permutation'][:]
+    x = fid.variables['x'][:] + fid.xllcorner   #x-coordinates of vertices
+    y = fid.variables['y'][:] + fid.yllcorner   #y-coordinates of vertices
+    points = num.transpose(num.asarray([x.tolist(), y.tolist()]))
+    time = fid.variables['time'][:] + fid.starttime
+    elevation = fid.variables['elevation'][:]
+       
+    basename = 'sts_gauge'
+    quantity_names = ['stage', 'xmomentum', 'ymomentum']
+    quantities = {}
+    for i, name in enumerate(quantity_names):
+        quantities[name] = fid.variables[name][:]
+
+    #---------------------------------------------------------------------------
+    # Get maximum wave height throughout timeseries at each index point
+    #---------------------------------------------------------------------------
+     
+    maxname = 'max_sts_stage.csv'
+    print 'get_sts_gauge_data: maxname=%s' % maxname
+    fid_max = open(os.path.join(project.event_folder, maxname), 'w')
+    fid_max.write('index, x, y, max_stage \n')    
+    for j in range(len(x)):
+        index = permutation[j]
+        stage = quantities['stage'][:,j]
+        xmomentum = quantities['xmomentum'][:,j]
+        ymomentum = quantities['ymomentum'][:,j]
+
+        fid_max.write('%d, %.6f, %.6f, %.6f\n' % (index, x[j], y[j], max(stage)))
+      
+    #---------------------------------------------------------------------------
+    # Get minimum wave height throughout timeseries at each index point
+    #---------------------------------------------------------------------------
+
+    minname = 'min_sts_stage.csv'
+    fid_min = open(os.path.join(project.event_folder, minname), 'w')
+    fid_min.write('index, x, y, max_stage \n')    
+    for j in range(len(x)):
+        index = permutation[j]
+        stage = quantities['stage'][:,j]
+        xmomentum = quantities['xmomentum'][:,j]
+        ymomentum = quantities['ymomentum'][:,j]
+
+        fid_min.write('%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage)))
+
+        out_file = os.path.join(project.event_folder,
+                                basename+'_'+str(index)+'.csv')
+        fid_sts = open(out_file, 'w')
+        fid_sts.write('time, stage, xmomentum, ymomentum \n')
+
+        #-----------------------------------------------------------------------
+        # End of the get gauges
+        #-----------------------------------------------------------------------
+        for k in range(len(time)-1):
+            fid_sts.write('%.6f, %.6f, %.6f, %.6f\n'
+                          % (time[k], stage[k], xmomentum[k], ymomentum[k]))
+
+        fid_sts.close()      
+
+    fid.close()
+
+    return quantities,elevation,time
+
+
+##
+# @brief Build boundary STS files from one or more MUX files.
+# @param event_file Name of mux meta-file or single mux stem.
+# @param output_dir Directory to write STS data to.
+# @note 'event_file' is produced by EventSelection.
+def build_urs_boundary(event_file, output_dir,urs_file):
+    '''Build a boundary STS file from a set of MUX files.'''
+    print 'reading', urs_file
+    # if we are using an EventSelection multi-mux file
+    if project.multi_mux:
+        # get the mux+weight data from the meta-file (in <boundaries>)
+        mux_event_file = os.path.join(project.event_folder, event_file)
+        print 'using multi-mux file', mux_event_file
+        try:
+            fd = open(mux_event_file, 'r')
+            mux_data = fd.readlines()
+            fd.close()
+        except IOError, e:
+            msg = 'File %s cannot be read: %s' % (mux_event_file, str(e))
+            raise Exception, msg
+        except:
+            raise
+
+        # first line of file is # filenames+weight in rest of file
+        num_lines = int(mux_data[0].strip())
+        mux_data = mux_data[1:]
+        print 'number of sources %d' % num_lines
+
+        # quick sanity check on input mux meta-file
+        if num_lines != len(mux_data):
+            msg = ('Bad file %s: %d data lines, but line 1 count is %d'
+                   % (event_file, len(mux_data), num_lines))
+            raise Exception, msg
+
+        # Create filename and weights lists.
+        # Must chop GRD filename just after '*.grd'.
+        mux_filenames = []
+        for line in mux_data:
+            muxname = line.strip().split()[0]
+            split_index = muxname.index('.grd')
+            muxname = muxname[:split_index+len('.grd')]
+            muxname = os.path.join(project.mux_data_folder, muxname)
+            mux_filenames.append(muxname)
+        
+        mux_weights = [float(line.strip().split()[1]) for line in mux_data]
+
+        # Call legacy function to create STS file.
+        print 'creating sts file'
+        urs2sts(mux_filenames,
+                basename_out=output_dir,
+                ordering_filename=urs_file,
+                weights=mux_weights,
+                central_meridian=project.central_meridian,
+                verbose=True)
+    else:                           # a single mux stem file, assume 1.0 weight
+        mux_file = os.path.join(project.event_folder, event_file)
+        mux_filenames = [mux_file]
+        print 'using single-mux file', mux_file
+
+        weight_factor = 1.0
+        mux_weights = weight_factor*num.ones(len(mux_filenames), num.Float)
+
+
+        # Create ordered sts file
+        urs2sts(mux_filenames,
+                basename_out=output_dir,
+                ordering_filename=urs_file,
+                weights=mux_weights,
+                central_meridian=project.central_meridian,
+                verbose=True)
+
+    # report on progress so far
+##    sts_file = os.path.join(project.event_folder, project.scenario_name)
+    quantities, elevation, time = get_sts_gauge_data(output_dir, verbose=False)
+    print len(elevation), len(quantities['stage'][0,:])
+
+    
+
Index: /anuga_work/production/australia_ph2/melbourne/export_results_max.py
===================================================================
--- /anuga_work/production/australia_ph2/melbourne/export_results_max.py	(revision 6692)
+++ /anuga_work/production/australia_ph2/melbourne/export_results_max.py	(revision 6692)
@@ -0,0 +1,121 @@
+"""
+Generates ascii grids of nominated areas -
+Input: sww file from run_perth.py
+       boundaries for grids from project.py
+Outputs: ascii grids of specified variables
+Stored in the 'outputs_dir' folder for respective .sww file
+
+Note:
+If producing a grid for the enitre extent cellsize should be greater than 30m
+If producing grids for inundation area resolution should be greater than mesh (ie ~22m)
+"""
+
+import project, os
+import sys
+from os.path import join
+from anuga.lib.maxasc.maxasc import MaxAsc
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from anuga.shallow_water.data_manager import sww2dem
+from os import sep
+
+
+directory = project.output_folder
+
+time_dir1 = '20090130_165532_run_final_0.6_27255_extend_dt_kvanputt'
+##time_dir2 = ''
+##time_dir3 = ''
+##time_dir4 = ''
+##time_dir5 = ''
+##time_dir6 = ''
+##
+time_dirs = [time_dir1] #5, time_dir4]# , time_dir4, time_dir5, time_dir6]
+  
+
+cellsize = 250
+##cellsize = 5
+
+#timestep = None    # None means no timestep!
+timestep = 0
+
+######
+# Set the special areas of interest.  If none, do: area='All'
+######
+
+#area = ['Bunbury', 'Busselton']  # strings must match keys in var_equations below
+area = ['All']      # 'All' means no special areas - the whole thing
+
+######
+# Define allowed variable names and associated equations to generate values.
+# This would not normally change.
+######
+var_equations = {'stage':     'stage',
+                 'momentum':  '(xmomentum**2 + ymomentum**2)**0.5',
+                 'depth':     'stage-elevation',
+                 'speed':     '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6)',
+                 'elevation': 'elevation' }
+
+# one or more key strings from var_equations above
+var = ['elevation']
+
+######
+# Start running the various conversions we require.
+######
+
+for which_var in var:
+    if which_var not in var_equations:
+        print 'Unrecognized variable name: %s' % which_var
+        break
+
+    for which_area in area:
+        if which_area == 'All':
+            easting_min = None
+            easting_max = None
+            northing_min = None
+            northing_max = None
+        else:
+            try:
+                easting_min = eval('project.xmin%s' % which_area)
+                easting_max = eval('project.xmax%s' % which_area)
+                northing_min = eval('project.ymin%s' % which_area)
+                northing_max = eval('project.ymax%s' % which_area)
+            except AttributeError:
+                print 'Unrecognized area name: %s' % which_area
+                break
+
+        for time_dir in time_dirs:
+
+            name1 = join(directory, time_dir, project.scenario_name)
+            #name2 = join(directory, time_dir, project.scenario_name)+'_time_39600_0'
+            #name3 = join(directory, time_dir, project.scenario_name)+'_time_79200_0'
+
+            names = [name1] #, name2]#, name3]
+      
+           # asc_name = []   
+               
+            for name in names:
+                         
+                outname = name + '_' + which_area + '_' + which_var
+                quantityname = var_equations[which_var]
+
+                print 'start sww2dem: time_dir=%s' % time_dir
+                
+                sww2dem(name, basename_out = outname,
+                            quantity = quantityname,
+                            timestep = timestep,
+                            cellsize = cellsize,      
+                            easting_min = easting_min,
+                            easting_max = easting_max,
+                            northing_min = northing_min,
+                            northing_max = northing_max,        
+                            reduction = max, 
+                            verbose = True,
+                            format = 'asc')
+
+##                asc_name.append(outname + '.asc')
+##    
+##            maxasc_outname = join(directory, time_dir, project.scenario_name)+'_'+which_area+'_'+which_var+'_max.asc'
+##
+##            print 'max asc outname ', maxasc_outname
+##            print 'asc_name ', str(asc_name)
+##            
+##            MaxAsc(maxasc_outname, asc_name)
Index: /anuga_work/production/australia_ph2/melbourne/file_length.py
===================================================================
--- /anuga_work/production/australia_ph2/melbourne/file_length.py	(revision 6692)
+++ /anuga_work/production/australia_ph2/melbourne/file_length.py	(revision 6692)
@@ -0,0 +1,22 @@
+"""Function to return the length of a file.
+
+Intended to be used for simplfying the creation of boundary_tags in
+run_model.py
+
+Input: filename
+Returns: number of lines in file
+"""
+
+def file_length(in_file):
+    '''Function to return the number of lines in a file.
+
+    in_file: Path to the file to get number of lines in.
+    
+    Returns: number of lines in file
+    '''
+    
+    fid = open(in_file)
+    data = fid.readlines()
+    fid.close()
+    return len(data)
+
Index: /anuga_work/production/australia_ph2/melbourne/prepare_timeboundary.py
===================================================================
--- /anuga_work/production/australia_ph2/melbourne/prepare_timeboundary.py	(revision 6692)
+++ /anuga_work/production/australia_ph2/melbourne/prepare_timeboundary.py	(revision 6692)
@@ -0,0 +1,68 @@
+def prepare_timeboundary(file_in):
+    """Convert benchmark 2 time series to NetCDF tms file.
+    input file is a csv created by build_boundary.py
+    """
+    from Numeric import array, zeros, Float, allclose
+    from Scientific.IO.NetCDF import NetCDFFile
+    from Numeric import array
+    import project
+    from os import sep
+
+    boundary_csv = file_in
+    boundary_tms = boundary_csv[:-4] + '.tms'
+    print 'Creating', boundary_tms
+
+    # Read the ascii (.txt) version of this file
+    fid = open(file_in)
+
+    # Skip first line
+    line = fid.readline()
+
+    # Read remaining lines
+    lines = fid.readlines()
+    fid.close()
+    
+
+    N = len(lines)
+    T = zeros(N, Float)  #Time
+    H = zeros(N, Float)  #Values
+    X = zeros(N, Float)  #Values
+    Y = zeros(N, Float)  #Values
+
+    for i, line in enumerate(lines):
+        fields = line.split(',')
+
+        T[i] = float(fields[0])  #time
+        H[i] = float(fields[1])  #stage
+        X[i] = float(fields[2])  #x-momentum
+        Y[i] = float(fields[3])  #y-momentum
+
+
+    # Create tms NetCDF file
+    
+    fid = NetCDFFile(boundary_tms, 'w')
+    fid.institution = 'Geoscience Australia'
+    fid.description = 'Input wave from index %i' %project.index
+    fid.starttime = T[0]
+    fid.createDimension('number_of_timesteps', len(T))
+    fid.createVariable('time', Float, ('number_of_timesteps',))
+    fid.variables['time'][:] = T-T[0]
+
+    fid.createVariable('stage', Float, ('number_of_timesteps',))
+    fid.variables['stage'][:] = H
+
+    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
+    fid.variables['xmomentum'][:] = X
+
+    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
+    fid.variables['ymomentum'][:] = Y
+
+    fid.close()
+
+
+
+##
+###-------------------------------------------------------------
+##if __name__ == "__main__":
+##
+
Index: /anuga_work/production/australia_ph2/melbourne/project.py
===================================================================
--- /anuga_work/production/australia_ph2/melbourne/project.py	(revision 6692)
+++ /anuga_work/production/australia_ph2/melbourne/project.py	(revision 6692)
@@ -0,0 +1,251 @@
+"""
+This file contains all your file and directory definitions 
+for elevation, meshes and outputs.
+"""
+
+import os
+from anuga.utilities.system_tools import get_user_name, get_host_name
+from time import localtime, strftime, gmtime
+from os.path import join, exists
+
+
+#-------------------------------------------------------------------------------
+# Directory setup
+#-------------------------------------------------------------------------------
+
+# this section needs to be updated to reflect the modelled community.
+# Note, the user needs to set up the directory system accordingly
+state = 'australia_ph2'
+scenario_name = 'melbourne'
+scenario_folder = scenario_name
+
+#-------------------------------------------------------------------------------
+# Initial Conditions
+#-------------------------------------------------------------------------------
+
+# Model specific parameters.
+# One or all can be changed each time the run_model script is executed
+tide = 0                # difference between MSL and HAT
+central_meridian = 146.0 # Central meridian for projection (optional)
+
+alpha = 0.1             # smoothing parameter for mesh
+friction=0.01           # manning's friction coefficient
+starttime=0             # start time for simulation
+finaltime=80000          # final time for simulation
+
+# index is only used when wave = Tb
+wave = 'Bf'             # Bf (sts wave) Tb (index wave)
+
+event_side = 'west'     # for the Tb wave, need to set which side is wave
+                        # (other ocean boundary will be set to transmissive)
+
+# for melbourne (west) side:
+event_number = 64288    # the event number or the mux file name
+# sources:
+# other possible: 
+index = 1982            # index from the PTHA - Y2000 is 0.43302m
+
+setup = 'final'         # This can be one of three values
+                        #    trial - coarsest mesh, fast
+                        #    basic - coarse mesh
+                        #    final - fine mesh, slowest
+
+internal_polygon = False
+
+#-------------------------------------------------------------------------------
+# Output filename
+#
+# Your output filename should be unique between different runs on different data.
+# The list of items below will be used to create a file in your output directory.
+# Your user name and time+date will be automatically added.  For example,
+#     [setup, tide, event_number]
+# will result in a filename like
+#     20090212_091046_run_final_0_27283_rwilson
+#-------------------------------------------------------------------------------
+
+if internal_polygon:
+    internal_poly_comment = 'internal'
+else:
+    internal_poly_comment = None
+
+output_comment = [setup, tide, event_number, index, wave, internal_poly_comment]
+
+#-------------------------------------------------------------------------------
+# Input Data
+#-------------------------------------------------------------------------------
+
+# ELEVATION DATA
+# Used in build_elevation.py
+# Format for ascii grids, as produced in ArcGIS + a projection file
+ascii_grid_filenames = [] # 250m grid 2005
+
+# Format for point is x,y,elevation (with header)
+point_filenames = ['grid_250m.txt'] # 250m grid 2005
+
+### Add csv header list to all files in point_filenames
+##headerlist = ['x', 'y', 'elevation'] 
+##for f in point_filenames:
+##    add_csv_header(join(topographies_folder, f), headerlist)
+
+    # BOUNDING POLYGON - for data clipping and estimate of triangles in mesh
+# Used in build_elevation.py
+# Format for points easting,northing (no header)
+bounding_polygon_filename = 'bounding_polygon.csv'
+bounding_polygon_maxarea = 100000
+
+# INTERIOR REGIONS -  for designing the mesh
+# Used in run_model.py
+# Format for points easting,northing (no header)                    
+interior_regions_data = []
+
+# add an internal polygon to force different mesh generation
+# used to test for discretisation error when building elevation
+# make sure file is in same folder as interior regions and bouding polygon;
+# format is same (2 column .csv; easting, northing; no header)
+if internal_polygon:
+    interior_regions_data.append(['internal_polygon.csv',
+                                  bounding_polygon_maxarea])
+
+# LAND - used to set the initial stage/water to be offcoast only
+# Used in run_model.py.  Format for points easting,northing (no header)
+land_initial_conditions_filename = []
+
+# GAUGES - for creating timeseries at a specific point
+# Used in get_timeseries.py.  
+# Format easting,northing,name,elevation (with header)
+gauges_filename = 'gauges.csv'
+
+# BUILDINGS EXPOSURE - for identifying inundated houses
+# Used in run_building_inundation.py
+# Format latitude,longitude etc (geographic)
+building_exposure_filename = 'busselton_res_clip.csv' # from NEXIS
+
+# BOUNDING POLYGON - used in build_boundary.py and run_model.py respectively
+# NOTE: when files are put together the points must be in sequence
+# For ease go clockwise!
+# Check the run_model.py for boundary_tags
+
+# Thinned ordering file from Hazard Map (geographic)
+# Format is index,latitude,longitude (with header)
+urs_order_filename_W = 'urs_order_west.csv'
+urs_order_filename_E ='urs_order_east.csv'
+##urs_order_filename = 'urs_order_both.csv' # east is first, then west
+
+# Landward bounding points
+# Format easting,northing (no header)
+landward_boundary_filename_N = 'landward_boundary_north.csv'
+landward_boundary_filename_S = 'landward_boundary_south.csv'
+
+# MUX input filename.
+# If a meta-file from EventSelection is used, set 'multi-mux' to True.
+# If a single MUX stem filename (*.grd) is used, set 'multi-mux' to False.
+##mux_input_filename = event_number # to be found in event_folder
+                                    # (ie boundaries/event_number/)
+##multi_mux = False
+mux_input_filename = 'event.list'
+multi_mux = True
+
+
+################################################################################
+################################################################################
+####         NOTE: NOTHING WOULD NORMALLY CHANGE BELOW THIS POINT.          ####
+################################################################################
+################################################################################
+
+# Get system user and host names.
+# These values can be used to distinguish between two similar runs by two
+# different users or runs by the same user on two different machines.
+user = get_user_name()
+host = get_host_name()
+
+# Environment variable names.
+# The inundation directory, not the data directory.
+ENV_INUNDATIONHOME = 'INUNDATIONHOME'
+
+# Path to MUX data
+ENV_MUXHOME = 'MUXHOME'
+
+#-------------------------------------------------------------------------------
+# Output Elevation Data
+#-------------------------------------------------------------------------------
+
+# Output filename for elevation
+# this is a combination of all the data generated in build_elevation.py
+combined_elevation_basename = scenario_name + '_combined_elevation'
+
+#-------------------------------------------------------------------------------
+# Directory Structure
+#-------------------------------------------------------------------------------
+
+# determines time for setting up output directories
+time = strftime('%Y%m%d_%H%M%S', localtime()) 
+gtime = strftime('%Y%m%d_%H%M%S', gmtime()) 
+build_time = time + '_build'
+run_time = time + '_run_'
+
+# create paths generated from environment variables.
+home = join(os.getenv(ENV_INUNDATIONHOME), 'data') # Absolute path for data folder
+muxhome = os.getenv(ENV_MUXHOME)
+    
+# check various directories/files that must exist
+anuga_folder = join(home, state, scenario_folder, 'anuga')
+topographies_folder = join(anuga_folder, 'topographies')
+polygons_folder = join(anuga_folder, 'polygons')
+boundaries_folder = join(anuga_folder, 'boundaries')
+output_folder = join(anuga_folder, 'outputs')
+gauges_folder = join(anuga_folder, 'gauges')
+meshes_folder = join(anuga_folder, 'meshes')
+event_folder = join(boundaries_folder, str(event_number))
+
+# MUX data files
+# Directory containing the MUX data files to be used with EventSelection.
+mux_data_folder = join(muxhome, 'mux')
+
+#-------------------------------------------------------------------------------
+# Location of input and output data
+#-------------------------------------------------------------------------------
+
+# Convert the user output_comment to a string for run_model.py
+output_comment = ('_'.join([str(x) for x in output_comment if x != user])
+                  + '_' + user)
+
+# The absolute pathname of the all elevation, generated in build_elevation.py
+combined_elevation = join(topographies_folder, combined_elevation_basename)
+
+# The absolute pathname of the mesh, generated in run_model.py
+meshes = join(meshes_folder, scenario_name) + '.msh'
+
+# The pathname for the urs order points, used within build_urs_boundary.py
+urs_order_east = join(boundaries_folder, urs_order_filename_E)
+urs_order_west = join(boundaries_folder, urs_order_filename_W)
+
+# The absolute pathname for the landward points of the bounding polygon,
+# Used within run_model.py)
+landward_boundary_N = join(boundaries_folder, landward_boundary_filename_N)
+landward_boundary_S = join(boundaries_folder, landward_boundary_filename_S)
+
+# The absolute pathname for the .sts file, generated in build_boundary.py
+event_sts_east = join(event_folder, scenario_name) + '_east'
+event_sts_west = join(event_folder, scenario_name) + '_west'
+
+# The absolute pathname for the output folder names
+# Used for build_elevation.py
+output_build = join(output_folder, build_time) + '_' + str(user) 
+# Used for run_model.py
+output_run = join(output_folder, run_time) + output_comment 
+# Used by post processing
+output_run_time = join(output_run, scenario_name) 
+
+# The absolute pathname for the gauges file
+# Used for get_timeseries.py
+gauges = join(gauges_folder, gauges_filename)       
+
+# The absolute pathname for the building file
+# Used for run_building_inundation.py
+building_exposure = join(gauges_folder, building_exposure_filename)
+
+# full path to where MUX files (or meta-files) live
+mux_input = join(event_folder, mux_input_filename)
+
+# sts gauge with desired index number - used for wave 'Tb'
+boundary_csv = join(event_folder, 'sts_gauge_' + str(index) +'.csv')
Index: /anuga_work/production/australia_ph2/melbourne/run_model.py
===================================================================
--- /anuga_work/production/australia_ph2/melbourne/run_model.py	(revision 6692)
+++ /anuga_work/production/australia_ph2/melbourne/run_model.py	(revision 6692)
@@ -0,0 +1,214 @@
+"""Run a tsunami inundation scenario for Busselton, WA, Australia.
+
+The scenario is defined by a triangular mesh created from project.polygon, the
+elevation data is compiled into a pts file through build_elevation.py and a
+simulated tsunami is generated through an sts file from build_boundary.py.
+
+Input: sts file (build_boundary.py for respective event)
+       pts file (build_elevation.py)
+       information from project file
+Outputs: sww file stored in project.output_run_time_dir 
+The export_results_all.py and get_timeseries.py is reliant
+on the outputs of this script
+
+Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006
+Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import os.path
+import time
+from time import localtime, strftime, gmtime
+
+# Related major packages
+from Scientific.IO.NetCDF import NetCDFFile
+import Numeric as num
+
+from anuga.interface import create_domain_from_regions
+from anuga.interface import Transmissive_stage_zero_momentum_boundary
+from anuga.interface import Dirichlet_boundary
+from anuga.interface import Reflective_boundary
+from anuga.interface import Field_boundary
+from anuga.interface import Time_boundary
+from anuga.interface import file_function
+
+from anuga.interface import create_sts_boundary
+from anuga.interface import csv2building_polygons
+from file_length import file_length
+
+from anuga.shallow_water.data_manager import start_screen_catcher
+from anuga.shallow_water.data_manager import copy_code_files
+from anuga.shallow_water.data_manager import urs2sts
+from anuga.utilities.polygon import read_polygon, Polygon_function
+
+# Application specific imports
+from setup_model import project
+import build_urs_boundary as bub
+import prepare_timeboundary as TB
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file. Copy script must be before screen_catcher
+#-------------------------------------------------------------------------------
+
+copy_code_files(project.output_run, __file__, 
+                os.path.join(os.path.dirname(project.__file__),
+                             project.__name__+'.py'))
+start_screen_catcher(project.output_run, 0, 1)
+
+#-------------------------------------------------------------------------------
+# Create the computational domain based on overall clipping polygon with
+# a tagged boundary and interior regions defined in project.py along with
+# resolutions (maximal area of per triangle) for each polygon
+#-------------------------------------------------------------------------------
+
+print 'Create computational domain'
+
+# Create the STS file
+print 'project.mux_data_folder=%s' % project.mux_data_folder
+if not os.path.exists(project.event_sts_east + '.sts'):
+    bub.build_urs_boundary(project.mux_input_filename, project.event_sts_east, project.urs_order_east)
+if not os.path.exists(project.event_sts_west + '.sts'):
+    bub.build_urs_boundary(project.mux_input_filename, project.event_sts_west, project.urs_order_west)
+
+# Read in boundary from ordered sts file
+event_sts_east = create_sts_boundary(project.event_sts_east)
+event_sts_west = create_sts_boundary(project.event_sts_west)
+
+# Reading the landward defined points, this incorporates the original clipping
+# polygon minus the 100m contour
+landward_boundary_N = read_polygon(project.landward_boundary_N)
+landward_boundary_S = read_polygon(project.landward_boundary_S)
+
+##print 'land N', landward_boundary_N, ' with length ', len(landward_boundary_N)
+##print 'land S', landward_boundary_S, ' with length ', len(landward_boundary_S)
+
+# Combine sts polyline with landward points
+bounding_polygon_sts = landward_boundary_N + event_sts_east + landward_boundary_S + event_sts_west
+
+print 'bounding polygon sts',bounding_polygon_sts,'length',len(bounding_polygon_sts)
+
+# 4 sides: north (land), east (ocean), south (land), west (ocean)
+num_north_points = len(landward_boundary_N)
+num_east_points = len(event_sts_east)
+num_south_points = len(landward_boundary_S)
+num_west_points = len(event_sts_west)
+
+# Boundary tags refer to project.landward_boundary
+# 4 points equals 5 segments start at N
+boundary_tags={'back': range(num_north_points-1) + range(num_north_points+num_east_points,num_north_points+num_east_points+num_south_points-1),
+               'side': [num_north_points-1,
+                        num_north_points+num_east_points-1,
+                        num_north_points+num_east_points+num_south_points-1,
+                        num_north_points+num_east_points+num_south_points+num_west_points-1],
+               'oceanE': range(num_north_points,num_north_points+num_east_points-1),
+               'oceanW': range(num_north_points+num_east_points+num_south_points,num_north_points+num_east_points+num_south_points+num_west_points-1)}
+
+##print 'tags',boundary_tags
+
+# Build mesh and domain
+domain = create_domain_from_regions(bounding_polygon_sts,
+                                    boundary_tags=boundary_tags,
+                                    maximum_triangle_area=project.bounding_maxarea,
+                                    interior_regions=project.interior_regions,
+                                    mesh_filename=project.meshes,
+                                    use_cache=True,
+                                    verbose=True)
+print domain.statistics()
+
+domain.set_name(project.scenario_name)
+domain.set_datadir(project.output_run) 
+domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
+
+#-------------------------------------------------------------------------------
+# Setup initial conditions
+#-------------------------------------------------------------------------------
+
+print 'Setup initial conditions'
+
+# Set the initial stage in the offcoast region only
+if project.land_initial_conditions:
+    IC = Polygon_function(project.land_initial_conditions,
+                          default=project.tide,
+                          geo_reference=domain.geo_reference)
+else:
+    IC = 0
+domain.set_quantity('stage', IC, use_cache=True, verbose=True)
+domain.set_quantity('friction', project.friction) 
+domain.set_quantity('elevation', 
+                    filename=project.combined_elevation+'.pts',
+                    use_cache=True,
+                    verbose=True,
+                    alpha=project.alpha)
+
+#-------------------------------------------------------------------------------
+# Setup boundary conditions 
+#-------------------------------------------------------------------------------
+
+print 'Set boundary - available tags:', domain.get_boundary_tags()
+
+# Prepare time boundary
+TB.prepare_timeboundary(project.boundary_csv)
+f = file_function(project.boundary_csv[:-4] + '.tms')
+
+Br = Reflective_boundary(domain)
+Bt = Transmissive_stage_zero_momentum_boundary(domain)
+Bd = Dirichlet_boundary([project.tide, 0, 0])
+
+if project.wave == 'Bf':
+    Bfe = Field_boundary(project.event_sts_east+'.sts',
+                        domain, mean_stage=project.tide,
+                        time_thinning=1,
+                        default_boundary=Dirichlet_boundary([0, 0, 0]),
+                        boundary_polygon=bounding_polygon_sts,                    
+                        use_cache=True,
+                        verbose=True)
+    Bfw = Field_boundary(project.event_sts_west+'.sts',
+                        domain, mean_stage=project.tide,
+                        time_thinning=1,
+                        default_boundary=Dirichlet_boundary([0, 0, 0]),
+                        boundary_polygon=bounding_polygon_sts,                    
+                        use_cache=True,
+                        verbose=True)
+    domain.set_boundary({'back': Br,
+                     'side': Bt,
+                     'oceanE': Bfe,
+                     'oceanW': Bfw}) 
+    
+elif project.wave == 'Tb':
+    Tb = Time_boundary(domain,f,default_boundary=Dirichlet_boundary([0, 0, 0]) )
+
+    if project.event_side == 'east' :
+        domain.set_boundary({'back': Br,
+                             'side': Bt,
+                             'oceanW': Bt,
+                             'oceanE': Tb})
+    elif project.event_side == 'west' :
+        domain.set_boundary({'back': Br,
+                             'side': Bt,
+                             'oceanW': Tb,
+                             'oceanE': Bt})
+    else:
+        print 'need to specify which side event is on'
+
+else:
+    print 'No wave specified in project script (Bf or Tb)'
+    
+
+#-------------------------------------------------------------------------------
+# Evolve system through time
+#-------------------------------------------------------------------------------
+
+t0 = time.time()
+for t in domain.evolve(yieldstep=project.yieldstep, 
+                       finaltime=project.finaltime,
+                       skip_initial_step=False): 
+    print domain.timestepping_statistics()
+    print domain.boundary_statistics(tags='ocean')
+
+print 'Simulation took %.2f seconds' % (time.time()-t0)
Index: /anuga_work/production/australia_ph2/melbourne/setup_model.py
===================================================================
--- /anuga_work/production/australia_ph2/melbourne/setup_model.py	(revision 6692)
+++ /anuga_work/production/australia_ph2/melbourne/setup_model.py	(revision 6692)
@@ -0,0 +1,160 @@
+"""
+Module to import the project.py 'configuration' file and perform
+sanity checks plus a quick check on mesh generation.
+
+Also callable as a stand-alone program, mainly to view the results
+of the mesh generation.
+"""
+
+import os
+
+from os.path import join, exists
+from anuga.utilities.polygon import read_polygon, number_mesh_triangles
+
+import project
+
+
+#-------------------------------------------------------------------------------
+# Sanity checks - check that required files/directories, etc, exist.
+#-------------------------------------------------------------------------------
+
+# flag - we check many things and then don't proceed if anything wrong
+sanity_error = False               # checked at bottom of this file
+
+# Test that environment variables are defined.
+if os.getenv(project.ENV_INUNDATIONHOME) is None:
+    print "Environment variable '%s' is not set" % project.ENV_INUNDATIONHOME
+    sanity_error = True
+
+if os.getenv(project.ENV_MUXHOME) is None:
+    print "Environment variable '%s' is not set" % project.ENV_MUXHOME
+    sanity_error = True
+
+#-------------------------------------------------------------------------------
+# Directory Structure
+#-------------------------------------------------------------------------------
+
+# check folders generated from environment variables.
+if not exists(project.home):
+    print "Sorry, data directory '%s' doesn't exist" % project.home
+    sanity_error = True
+    
+if not exists(project.muxhome):
+    print "Sorry, MUX directory '%s' doesn't exist" % project.muxhome
+    sanity_error = True
+    
+if not exists(project.anuga_folder):
+    print "Sorry, ANUGA directory '%s' doesn't exist" % project.anuga_folder
+    sanity_error = True
+
+if not exists(project.topographies_folder):
+    print ("Sorry, topo directory '%s' doesn't exist"
+           % project.topographies_folder)
+    sanity_error = True
+
+if not exists(project.polygons_folder):
+    print ("Sorry, polygon directory '%s' doesn't exist"
+           % project.polygons_folder)
+    sanity_error = True
+
+if not exists(project.boundaries_folder):
+    print ("Sorry, boundaries directory '%s' doesn't exist"
+           % project.boundaries_folder)
+    sanity_error = True
+
+if not exists(project.output_folder):
+    print "Sorry, outputs directory '%s' doesn't exist" % project.output_folder
+    sanity_error = True
+
+if not exists(project.gauges_folder):
+    print "Sorry, gauges directory '%s' doesn't exist" % project.gauges_folder
+    sanity_error = True
+
+if not exists(project.meshes_folder):
+    print "Sorry, meshes directory '%s' doesn't exist" % project.meshes_folder
+    sanity_error = True
+
+if not exists(project.mux_data_folder):
+    print "Sorry, mux data directory '%s' doesn't exist" % project.mux_data_folder
+    sanity_error = True
+
+# if multi_mux is True, check if multi-mux file exists
+if project.multi_mux:
+    if not exists(project.mux_input):
+        print ("Sorry, MUX input file '%s' doesn't exist"
+               % project.mux_input)
+        sanity_error = True
+
+#-----
+# If this directory don't exist, EventSelection hasn't been run.
+#-----
+
+if not exists(project.event_folder):
+    print ("Sorry, you must generate event %s with EventSelection."
+           % project.event_number)
+    sanity_error = True
+
+#-------------------------------------------------------------------------------
+# Determine type of run
+#-------------------------------------------------------------------------------
+
+if project.setup == 'trial':
+    print 'trial'
+    project.scale_factor = 100
+    project.time_thinning = 96
+    project.yieldstep = 240
+elif project.setup == 'basic': 
+    print 'basic'
+    project.scale_factor = 4
+    project.time_thinning = 12
+    project.yieldstep = 120
+elif project.setup == 'final': 
+    print 'final'
+    project.scale_factor = 1
+    project.time_thinning = 4
+    project.yieldstep = 60
+else:
+    print ("Sorry, you must set the 'setup' variable to one of:"
+           '   trial - coarsest mesh, fast\n'
+           '   basic - coarse mesh\n'
+           '   final - fine mesh, slowest\n'
+           '\n'
+           "'setup' was set to '%s'" % project.setup)
+    sanity_error = True
+
+#-------------------------------------------------------------------------------
+# Check for errors detected above.
+#-------------------------------------------------------------------------------
+
+if sanity_error:
+    msg = 'You must fix the above errors before continuing.'
+    raise Exception, msg
+
+#-------------------------------------------------------------------------------
+# Reading polygons and creating interior regions
+#-------------------------------------------------------------------------------
+
+# Create list of land polygons with initial conditions
+project.land_initial_conditions = []
+for filename, MSL in project.land_initial_conditions_filename:
+    polygon = read_polygon(join(project.polygons_folder, filename))
+    project.land_initial_conditions.append([polygon, MSL])
+
+# Create list of interior polygons with scaling factor
+project.interior_regions = []
+for filename, maxarea in project.interior_regions_data:
+    polygon = read_polygon(join(project.polygons_folder, filename))
+    project.interior_regions.append([polygon,
+                                     maxarea*project.scale_factor])
+
+# Initial bounding polygon for data clipping 
+project.bounding_polygon = read_polygon(join(project.polygons_folder,
+                                             project.bounding_polygon_filename))
+project.bounding_maxarea = project.bounding_polygon_maxarea*project.scale_factor
+
+# Estimate the number of triangles                     
+trigs_min = number_mesh_triangles(project.interior_regions,
+                                  project.bounding_polygon,
+                                  project.bounding_maxarea)
+
+print 'min estimated number of triangles', trigs_min
Index: /anuga_work/production/busselton/busselton_rerun/get_timeseries.py
===================================================================
--- /anuga_work/production/busselton/busselton_rerun/get_timeseries.py	(revision 6692)
+++ /anuga_work/production/busselton/busselton_rerun/get_timeseries.py	(revision 6692)
@@ -0,0 +1,34 @@
+"""
+Generate time series of nominated "gauges".
+Input: sww file from run_perth.py
+       gauges project.gauge_filename
+Output: csv files stage,speed,depth,elevation over time
+Stored in the 'outputs_dir' folder for respective .sww file
+
+Note:
+Can run different sww files at the same time
+make sure if there is a second sww file that it is placed into a folder called sww2
+Can run different gauges at the same time - ie testing boundary index point
+"""
+
+from os import sep, rename, listdir, system
+from os.path import join
+from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges,csv2timeseries_graphs
+import project
+
+directory = project.output_folder
+
+time_dir1 = '20090324_134328_run_final_0.6_27255_Bt_kvanputt'
+
+time_dirs = [time_dir1] #1, time_dir2, time_dir3, time_dir4, time_dir5, time_dir6]
+
+print 'time directories', time_dirs
+
+for time_dir in time_dirs:
+    name = join(directory,time_dir,project.scenario_name)
+    gauge = project.gauges
+    sww2csv_gauges(name+'.sww',gauge,
+                   quantities = ['stage','speed','depth','elevation'],
+                   verbose=True)
+  
+
Index: /anuga_work/production/busselton/busselton_rerun/project.py
===================================================================
--- /anuga_work/production/busselton/busselton_rerun/project.py	(revision 6691)
+++ /anuga_work/production/busselton/busselton_rerun/project.py	(revision 6692)
@@ -103,5 +103,5 @@
 # Used in get_timeseries.py.  
 # Format easting,northing,name,elevation (with header)
-gauges_filename = 'gauges.csv'
+gauges_filename = 'MH_gauges.csv'
 
 # BUILDINGS EXPOSURE - for identifying inundated houses
Index: /anuga_work/production/carnarvon/carnarvon_rerun/get_timeseries.py
===================================================================
--- /anuga_work/production/carnarvon/carnarvon_rerun/get_timeseries.py	(revision 6692)
+++ /anuga_work/production/carnarvon/carnarvon_rerun/get_timeseries.py	(revision 6692)
@@ -0,0 +1,34 @@
+"""
+Generate time series of nominated "gauges".
+Input: sww file from run_perth.py
+       gauges project.gauge_filename
+Output: csv files stage,speed,depth,elevation over time
+Stored in the 'outputs_dir' folder for respective .sww file
+
+Note:
+Can run different sww files at the same time
+make sure if there is a second sww file that it is placed into a folder called sww2
+Can run different gauges at the same time - ie testing boundary index point
+"""
+
+from os import sep, rename, listdir, system
+from os.path import join
+from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges,csv2timeseries_graphs
+import project
+
+directory = project.output_folder
+
+time_dir1 = '20090324_134442_run_final_1_27255_Bt_kvanputt'
+
+time_dirs = [time_dir1] #1, time_dir2, time_dir3, time_dir4, time_dir5, time_dir6]
+
+print 'time directories', time_dirs
+
+for time_dir in time_dirs:
+    name = join(directory,time_dir,project.scenario_name)
+    gauge = project.gauges
+    sww2csv_gauges(name+'.sww',gauge,
+                   quantities = ['stage','speed','depth','elevation'],
+                   verbose=True)
+  
+
Index: /anuga_work/production/carnarvon/carnarvon_rerun/project.py
===================================================================
--- /anuga_work/production/carnarvon/carnarvon_rerun/project.py	(revision 6691)
+++ /anuga_work/production/carnarvon/carnarvon_rerun/project.py	(revision 6692)
@@ -98,5 +98,5 @@
 # Used in get_timeseries.py.  
 # Format easting,northing,name,elevation (with header)
-gauges_filename = 'MH_gauges_deep.csv' #'carnarvon.csv'
+gauges_filename = 'MH_gauges.csv' #'carnarvon.csv'
 gauges_filename2 = 'thinned_MGA50.csv'
 
Index: /anuga_work/production/geraldton/geraldton_rerun/get_timeseries_myall.py
===================================================================
--- /anuga_work/production/geraldton/geraldton_rerun/get_timeseries_myall.py	(revision 6692)
+++ /anuga_work/production/geraldton/geraldton_rerun/get_timeseries_myall.py	(revision 6692)
@@ -0,0 +1,34 @@
+"""
+Generate time series of nominated "gauges".
+Input: sww file from run_perth.py
+       gauges project.gauge_filename
+Output: csv files stage,speed,depth,elevation over time
+Stored in the 'outputs_dir' folder for respective .sww file
+
+Note:
+Can run different sww files at the same time
+make sure if there is a second sww file that it is placed into a folder called sww2
+Can run different gauges at the same time - ie testing boundary index point
+"""
+
+from os import sep, rename, listdir, system
+from os.path import join
+from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges,csv2timeseries_graphs
+import project
+
+directory = project.output_folder
+
+time_dir1 = '20090324_134504_run_final_0.6_27255_Bt_kvanputt'
+
+time_dirs = [time_dir1] #1, time_dir2, time_dir3, time_dir4, time_dir5, time_dir6]
+
+print 'time directories', time_dirs
+
+for time_dir in time_dirs:
+    name = join(directory,time_dir,project.scenario_name)
+    gauge = project.gauges_myall
+    sww2csv_gauges(name+'.sww',gauge,
+                   quantities = ['stage','speed','depth','elevation'],
+                   verbose=True)
+  
+
Index: /anuga_work/production/geraldton/geraldton_rerun/project.py
===================================================================
--- /anuga_work/production/geraldton/geraldton_rerun/project.py	(revision 6691)
+++ /anuga_work/production/geraldton/geraldton_rerun/project.py	(revision 6692)
@@ -105,4 +105,5 @@
 #gauges_filename = 'MH_gauges.csv'
 gauges_filename = 'nick_geraldton_gauges.csv'
+gauges_myall_filename = 'MH_gauges.csv'
 
 # BUILDINGS EXPOSURE - for identifying inundated houses
@@ -244,4 +245,7 @@
     gauges = join(gauges_folder, gauges_filename)       
 
+if gauges_myall_filename:
+    gauges_myall = join(gauges_folder, gauges_myall_filename)  
+
 # The absolute pathname for the building file
 # Used for run_building_inundation.py
Index: /anuga_work/production/perth/perth_rerun/get_timeseries.py
===================================================================
--- /anuga_work/production/perth/perth_rerun/get_timeseries.py	(revision 6692)
+++ /anuga_work/production/perth/perth_rerun/get_timeseries.py	(revision 6692)
@@ -0,0 +1,34 @@
+"""
+Generate time series of nominated "gauges".
+Input: sww file from run_perth.py
+       gauges project.gauge_filename
+Output: csv files stage,speed,depth,elevation over time
+Stored in the 'outputs_dir' folder for respective .sww file
+
+Note:
+Can run different sww files at the same time
+make sure if there is a second sww file that it is placed into a folder called sww2
+Can run different gauges at the same time - ie testing boundary index point
+"""
+
+from os import sep, rename, listdir, system
+from os.path import join
+from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges,csv2timeseries_graphs
+import project
+
+directory = project.output_folder
+
+time_dir1 = '20090324_135427_run_final_0.6_27255_Bt_kvanputt'
+
+time_dirs = [time_dir1] #1, time_dir2, time_dir3, time_dir4, time_dir5, time_dir6]
+
+print 'time directories', time_dirs
+
+for time_dir in time_dirs:
+    name = join(directory,time_dir,project.scenario_name)
+    gauge = project.gauges
+    sww2csv_gauges(name+'.sww',gauge,
+                   quantities = ['stage','speed','depth','elevation'],
+                   verbose=True)
+  
+
