Index: /trunk/anuga_core/source/anuga/utilities/quantity_setting_functions.py
===================================================================
--- /trunk/anuga_core/source/anuga/utilities/quantity_setting_functions.py	(revision 9340)
+++ /trunk/anuga_core/source/anuga/utilities/quantity_setting_functions.py	(revision 9341)
@@ -198,9 +198,9 @@
 
             # Get indices fInds of points in pi which are not set
-            if(pi is 'All'):
+            if(pi == 'All'):
                 fInside=(1-isSet)
                 fInds=(fInside==1).nonzero()[0]
             else:
-                if(pi is 'Extent' and type(fi) is str and os.path.exists(fi)):
+                if(pi == 'Extent' and type(fi) is str and os.path.exists(fi)):
                     # Here fi MUST be a gdal-compatible raster
                     # Then we get the extent from the raster itself
@@ -250,5 +250,5 @@
                 isSet[fInds]=1
 
-        if(not min(isSet)==1):
+        if( min(isSet) != 1):
             raise Exception, 'Some points were not inside any polygon'
 
Index: unk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/getxls.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/getxls.py	(revision 9340)
+++ 	(revision )
@@ -1,364 +1,0 @@
-#!/usr/bin/python
-# -*- coding: utf-8 -*-
-
-"""
-
-Read xls(x) data to define ANUGA scenario
-
-"""
-
-
-import xlrd
-
-
-class AnugaProject:
-
-    """Class to hold data that typically occurs in project.py
-
-    We read the data from an excel file via the AnugaXls class
-
-    """
-
-    def __init__(self, filename):
-        """Read the excel data configuration file used to run ANUGA into a
-           class
-
-        """
-
-        # Parse the xls data
-
-        self.anuga_xls = AnugaXls(filename)
-        anuga_xls = self.anuga_xls
-
-        # Worksheet names
-        project_ws = 'project_settings'
-        mesh_ws = 'mesh'
-        init_ws = 'initial_conditions'
-        init_add_ws = 'initial_condition_additions'
-        # #####################################################################
-        #
-        # Generic 'project' type variables
-        #
-        # #####################################################################
-
-        self.scenario = anuga_xls.get_var(project_ws, 'scenario', [0, 1])
-
-        self.output_basedir = anuga_xls.get_var(
-            project_ws, 'output_base_directory', [
-                0, 1])
-
-        self.yieldstep = anuga_xls.get_var(project_ws, 'yieldstep', [0, 1],
-                                           post_process=float)
-
-        self.finaltime = anuga_xls.get_var(project_ws, 'finaltime', [0, 1],
-                                           post_process=float)
-
-        self.projection_information = anuga_xls.get_var(
-            project_ws, 'projection_information', [
-                0, 1])
-
-        # Could be unicode string or float -- if float, we need int
-
-        if isinstance(self.projection_information, float):
-            self.projection_information = int(self.projection_information)
-
-        self.flow_algorithm = anuga_xls.get_var(project_ws, 'flow_algorithm',
-                                                [0, 1])
-
-        # Coerce this to a logical variable
-
-        self.use_local_extrapolation_and_flux_updating = \
-            anuga_xls.get_var(project_ws,
-                              'use_local_extrapolation_and_flux_updating',
-                              offset=[0, 1], post_process=bool)
-
-        self.output_tif_cellsize = anuga_xls.get_var(
-            project_ws,
-            'output_tif_cellsize',
-            offset=[0, 1],
-            post_process=int)
-
-        self.output_tif_bounding_polygon = anuga_xls.get_var(
-            project_ws,
-            'output_tif_bounding_polygon',
-            offset=[0, 1],
-            post_process=empty_string_to_none)
-
-        self.max_quantity_update_frequency = anuga_xls.get_var(
-            project_ws,
-            'max_quantity_update_frequency',
-            offset=[0, 1],
-            post_process=int)
-
-        self.max_quantity_collection_starttime = anuga_xls.get_var(
-            project_ws,
-            'max_quantity_collection_starttime',
-            offset=[0, 1])
-
-        self.store_vertices_uniquely = anuga_xls.get_var(
-            project_ws,
-            'store_vertices_uniquely',
-            offset=[0, 1],
-            post_process=bool)
-
-        # #####################################################################
-        #
-        # Mesh type variables
-        #
-        # #####################################################################
-
-        self.bounding_polygon_and_tags_file = anuga_xls.get_var(
-            mesh_ws, 'bounding_polygon', [
-                1, 1])
-
-        self.default_res = anuga_xls.get_var(mesh_ws, 'bounding_polygon',
-                                             [1, 2])
-
-        self.boundary_tags_attribute_name = anuga_xls.get_var(
-            mesh_ws, 'bounding_polygon', [2, 2])
-
-        # Loop to read the interior regions data
-        # Finish when we hit an empty cell
-
-        self.interior_regions_data = anuga_xls.get_paired_list(
-            mesh_ws, 'interior_regions', [1, 1])
-
-        # breaklines, riverwalls, point movement threshold
-
-        self.breakline_files = anuga_xls.get_list(mesh_ws, 'breakline_files',
-                                                  [1, 1])
-
-        self.riverwall_csv_files = anuga_xls.get_list(
-            mesh_ws, 'riverwall_csv_files', [1, 1])
-
-        self.break_line_intersect_point_movement_threshold = \
-            anuga_xls.get_var(mesh_ws, 'breakline_intersection_threshold',
-                              [1, 1])
-
-        self.pt_areas = anuga_xls.get_var(
-            mesh_ws, 'region_areas_file', [0, 1])
-
-        # Check it is reasonable
-
-        if self.interior_regions_data != []:
-            msg = 'Cannot use both interior_regions_data and breaklines ' \
-                + 'riverwalls'
-            assert self.breakline_files == [], msg
-            assert self.riverwall_csv_files == [], msg
-
-        # #####################################################################
-        #
-        # Initial conditions type variables
-        #
-        # #####################################################################
-
-        self.elevation_data = anuga_xls.get_paired_list(
-            init_ws, 'Elevation', [1, 1])
-
-        self.friction_data = anuga_xls.get_paired_list(
-            init_ws, 'Friction', [1, 1])
-
-        self.stage_data = anuga_xls.get_paired_list(
-            init_ws, 'Stage', [1, 1])
-
-        self.xmomentum_data = anuga_xls.get_paired_list(
-            init_ws, 'X velocity x depth', [1, 1])
-
-        self.ymomentum_data = anuga_xls.get_paired_list(
-            init_ws, 'Y velocity x depth', [1, 1])
-
-        self.elevation_additions = anuga_xls.get_paired_list(
-            init_add_ws, 'Elevation', [1, 1])
-
-        self.friction_additions = anuga_xls.get_paired_list(
-            init_add_ws, 'Friction', [1, 1])
-
-        self.stage_additions = anuga_xls.get_paired_list(
-            init_add_ws, 'Stage', [1, 1])
-
-        self.xmomentum_additions = anuga_xls.get_paired_list(
-            init_add_ws, 'X velocity x depth', [1, 1])
-
-        self.ymomentum_additions = anuga_xls.get_paired_list(
-            init_add_ws, 'Y velocity x depth', [1, 1])
-
-        # #####################################################################
-        #
-        # Boundary condition type variables
-        #
-        # #####################################################################
-
-
-###############################################################################
-
-class AnugaXls:
-
-    """Read an xls or xlsx file with the ANUGA input data
-
-    This holds a list of worksheet_names, and
-    a dict containing the worksheet_data as lists
-
-    """
-
-    def __init__(self, filename):
-
-        workbook = xlrd.open_workbook(filename)
-
-        self.worksheet_names = workbook.sheet_names()
-
-        self.worksheet_data = {}
-
-        for worksheet_name in self.worksheet_names:
-            worksheet = workbook.sheet_by_name(worksheet_name)
-            self.worksheet_data[worksheet_name] = \
-                get_worksheet_as_list(worksheet)
-
-        return
-
-    def get_var(
-        self,
-        worksheet_name,
-        flag,
-        offset=[0, 0],
-        post_process=None,
-    ):
-        """Find the cell in the first column of worksheet_name that (partially)
-        matches 'flag', and return the nearby cell value which is 'offset'
-        from this by [nrow, ncol]
-
-        Optionally, post_process is a function applied to the cell value
-        before returning. This can be useful to e.g. coerce integers to
-        bool
-
-        """
-
-        sheet_data = self.worksheet_data[worksheet_name]
-
-        try:
-            for i in range(len(sheet_data)):
-                if flag in sheet_data[i][0]:
-
-                    # Check this would not exceed the number of rows
-
-                    assert i + offset[0] < len(sheet_data)
-
-                    # If the col offset exceeds the data length, return a blank
-
-                    if offset[1] < len(sheet_data[i]):
-                        output = sheet_data[i + offset[0]][offset[1]]
-                        if post_process is not None:
-                            output = post_process(output)
-                        return output
-                    else:
-                        return ''
-        except:
-            pass
-
-        msg = 'Could not find a cell in the first column of worksheet *' \
-            + worksheet_name + '* matching with *' + flag + '*'
-        raise Exception(msg)
-
-        return
-
-    def get_paired_list(
-        self,
-        worksheet_name,
-        flag,
-        offset=[0, 0],
-    ):
-        """Find the cell in the first column of worksheet_name that (partially)
-        matches 'flag', and read data from 2 rows of the spreadsheet,
-        starting at a location offset from the matching cell by
-        offset = [nrow, ncol].
-
-        The data is returned as a paired list:
-        [ [r00, r10], [r01, r11], [r02, r12]]
-
-        """
-
-        paired_list = []
-
-        counter = -1
-        while True:
-            counter = counter + 1
-            cell_x = offset[0]
-            cell_y = counter + offset[1]
-            next_row1 = self.get_var(worksheet_name, flag, [cell_x, cell_y])
-            if next_row1 == '':
-
-                # We are at the end of the data
-
-                break
-            else:
-                cell_x = offset[0] + 1
-                next_row2 = self.get_var(worksheet_name, flag, [cell_x,
-                                                                cell_y])
-
-            paired_list.append([next_row1, next_row2])
-
-        return paired_list
-
-    def get_list(
-        self,
-        worksheet_name,
-        flag,
-        offset=[0, 0],
-    ):
-        """Find the cell in the first column of worksheet_name that (partially)
-        matches 'flag', and read data from 1 row of the spreadsheet,
-        starting at a location offset from the matching cell by
-        offset = [nrow, ncol].
-
-        The data is returned as a list:
-        [ r00, r01, r02, r03, ...]
-
-        """
-
-        out_list = []
-
-        counter = -1
-        while True:
-            counter = counter + 1
-            cell_x = offset[0]
-            cell_y = counter + offset[1]
-            next_row1 = self.get_var(worksheet_name, flag, [cell_x, cell_y])
-            if next_row1 == '':
-
-                # We are at the end of the data
-
-                break
-            out_list.append(next_row1)
-
-        return out_list
-
-
-def get_worksheet_as_list(worksheet):
-    """Read a particular xls (xlsx) file worksheet as a list of lists
-
-    [
-    [cell 0 0, cell 0 1, cell 0 2, ...],
-    [cell 1 0, cell 1 1, cell 1 2, ...],
-    ...
-    ]
-
-    """
-
-    worksheet_data = []
-
-    for i in range(worksheet.nrows):
-        worksheet_data.append([])
-        for j in range(worksheet.ncols):
-            worksheet_data[i].extend([worksheet.cell_value(i, j)])
-
-    return worksheet_data
-
-
-def empty_string_to_none(string):
-    """Given a string, return None if string == "", and the string value
-    otherwise
-    """
-
-    if string == '':
-        return None
-    else:
-        return string
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/parse_input_data.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/parse_input_data.py	(revision 9341)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/parse_input_data.py	(revision 9341)
@@ -0,0 +1,500 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+"""
+
+Read xls(x) data to define ANUGA scenario
+
+"""
+import os
+import xlrd
+
+class ProjectData(object):
+
+    """Class to hold input data that previously occurred in project.py
+
+    We read the data from an excel file via the AnugaXls class
+
+    """
+
+    def __init__(self, filename):
+        """Call the right method to read the filetype
+
+        """
+
+        file_extension = os.path.splitext(filename)[1]
+
+        if(file_extension in ['.xls', '.xlsx']):
+            self.get_data_from_excel(filename)
+        else:
+            msg = 'File type not supported'
+            raise ValueError(msg)
+    
+
+    def get_data_from_excel(self, filename):
+        """Get data from xls or xlsx file
+
+        """
+
+        # Parse the xls data
+
+        self.data_source = AnugaXls(filename)
+        data_source = self.data_source
+
+        # Worksheet names -- define here so they can easily be changed
+        project_ws = 'project_settings'
+        mesh_ws = 'mesh'
+        init_ws = 'initial_conditions'
+        init_add_ws = 'initial_condition_additions'
+        boundary_ws = 'boundary_conditions'
+        inlet_ws = 'inlets'
+        rain_ws = 'rainfall'
+
+        # #####################################################################
+        #
+        # Generic 'project' type variables
+        #
+        # #####################################################################
+
+        self.scenario = data_source.get_var(
+            project_ws, 'scenario', [0, 1], post_process = str)
+
+        self.output_basedir = data_source.get_var(
+            project_ws, 'output_base_directory', [0, 1], post_process= str)
+
+        self.yieldstep = data_source.get_var(project_ws, 'yieldstep', [0, 1],
+                                           post_process=float)
+
+        self.finaltime = data_source.get_var(project_ws, 'finaltime', [0, 1],
+                                           post_process=float)
+
+        self.projection_information = data_source.get_var(
+            project_ws, 'projection_information', [0, 1])
+
+        # Could be unicode string or float -- if float, we need int
+
+        if isinstance(self.projection_information, float):
+            self.projection_information = int(self.projection_information)
+        if type(self.projection_information) is unicode:
+            self.projection_information = str(self.projection_information)
+
+        self.flow_algorithm = data_source.get_var(project_ws, 'flow_algorithm',
+                                                [0, 1], post_process = str)
+
+        # Coerce this to a logical variable
+
+        self.use_local_extrapolation_and_flux_updating = \
+            data_source.get_var(project_ws,
+                              'use_local_extrapolation_and_flux_updating',
+                              offset=[0, 1], post_process=bool)
+
+        self.output_tif_cellsize = data_source.get_var(
+            project_ws,
+            'output_tif_cellsize',
+            offset=[0, 1],
+            post_process=float)
+
+        self.output_tif_bounding_polygon = data_source.get_var(
+            project_ws,
+            'output_tif_bounding_polygon',
+            offset=[0, 1],
+            post_process=empty_string_to_none)
+
+        self.max_quantity_update_frequency = data_source.get_var(
+            project_ws,
+            'max_quantity_update_frequency',
+            offset=[0, 1],
+            post_process=int)
+
+        self.max_quantity_collection_start_time = data_source.get_var(
+            project_ws,
+            'max_quantity_collection_starttime',
+            offset=[0, 1],
+            post_process = float)
+
+        self.store_vertices_uniquely = data_source.get_var(
+            project_ws,
+            'store_vertices_uniquely',
+            offset=[0, 1],
+            post_process=bool)
+
+        self.spatial_text_output_dir = data_source.get_var(
+            project_ws,
+            'spatial_text_output_dir',
+            offset = [0, 1],
+            post_process=str)
+
+        # #####################################################################
+        #
+        # Mesh type variables
+        #
+        # #####################################################################
+
+        self.bounding_polygon_and_tags_file = data_source.get_var(
+            mesh_ws, 'bounding_polygon', [1, 1], post_process=str)
+
+        self.default_res = data_source.get_var(mesh_ws, 'bounding_polygon',
+                                             [2, 1], post_process = str)
+
+        # Loop to read the interior regions data
+        # Finish when we hit an empty cell
+
+        self.interior_regions_data = data_source.get_paired_list(
+            mesh_ws, 'interior_regions', [1, 1], 
+            post_process = string_or_float)
+
+        # breaklines, riverwalls, point movement threshold
+
+        self.breakline_files = data_source.get_list(mesh_ws, 'breakline_files',
+                                                  [0, 1], post_process = str)
+
+        self.riverwall_csv_files = data_source.get_list(
+            mesh_ws, 'riverwall_csv_files', [0, 1], post_process = str)
+
+        self.break_line_intersect_point_movement_threshold = \
+            data_source.get_var(mesh_ws, 'breakline_intersection_threshold',
+                              [0, 1], post_process = float)
+
+        self.pt_areas = data_source.get_var(
+            mesh_ws, 'region_areas_file', [0, 1], empty_string_to_none)
+
+        # Check we are only using interior_region OR breakline methods
+        use_ir = self.interior_regions_data != []
+        use_bl = (
+                  (self.breakline_files != []) or\
+                  ((self.riverwall_csv_files != []) or\
+                   (self.pt_areas is not None))
+                 )
+
+        if (use_ir and use_bl):
+            msg = ' In worksheet ' + mesh_ws +\
+                 'Cannot have both *interior_regions_data* and *breaklines ' \
+                + '/ riverwalls / region point areas* being non-empty.'
+            raise ValueError(msg)
+
+        # #####################################################################
+        #
+        # Initial conditions type variables
+        #
+        # #####################################################################
+
+        self.elevation_data = data_source.get_paired_list(
+            init_ws, 'Elevation', [1, 1], post_process = string_or_float)
+
+        self.friction_data = data_source.get_paired_list(
+            init_ws, 'Friction', [1, 1], post_process = string_or_float)
+
+        self.stage_data = data_source.get_paired_list(
+            init_ws, 'Stage', [1, 1], post_process = string_or_float)
+
+        self.xmomentum_data = data_source.get_paired_list(
+            init_ws, 'X-velocity_times_depth', [1, 1], 
+            post_process = string_or_float)
+
+        self.ymomentum_data = data_source.get_paired_list(
+            init_ws, 'Y-velocity_times_depth', [1, 1], 
+            post_process = string_or_float)
+
+        self.elevation_additions = data_source.get_paired_list(
+            init_add_ws, 'Elevation', [1, 1], 
+            post_process = string_or_float)
+
+        self.friction_additions = data_source.get_paired_list(
+            init_add_ws, 'Friction', [1, 1], 
+            post_process = string_or_float)
+
+        self.stage_additions = data_source.get_paired_list(
+            init_add_ws, 'Stage', [1, 1], 
+            post_process = string_or_float)
+
+        self.xmomentum_additions = data_source.get_paired_list(
+            init_add_ws, 'X-velocity_times_depth', [1, 1],
+            post_process = string_or_float)
+
+        self.ymomentum_additions = data_source.get_paired_list(
+            init_add_ws, 'Y-velocity_times_depth', [1, 1],
+            post_process = string_or_float)
+
+        # #####################################################################
+        #
+        # Boundary condition type variables
+        #
+        # #####################################################################
+
+        # This can fail if we don't coerce to str
+        self.boundary_tags_attribute_name = data_source.get_var(
+            boundary_ws, 'boundary_tags_attribute_name', [0, 1], str)
+
+        self.boundary_data = data_source.get_subtable(
+            boundary_ws, 'Boundary_Condition_Tag', [1, 0],
+            post_process = string_or_float)
+
+
+        # #####################################################################
+        #
+        # Inlet type variables
+        #
+        # #####################################################################
+
+        self.inlet_data = data_source.get_subtable(
+            inlet_ws, 'Inlet_name', [1, 0])
+
+        # #####################################################################
+        #
+        # Rainfall type variables
+        #
+        # #####################################################################
+
+        self.rain_data = data_source.get_subtable(
+            rain_ws, 'Rainfall', [1, 1])
+
+        # #####################################################################
+
+
+
+        
+
+# #############################################################################
+#
+# END OF CLASS
+#
+# #############################################################################
+
+class AnugaXls(object):
+
+    """Read an xls or xlsx file with the ANUGA input data
+
+    This holds a list of worksheet_names, and
+    a dict containing the worksheet_data as lists
+
+    """
+
+    def __init__(self, filename):
+
+        workbook = xlrd.open_workbook(filename)
+
+        self.worksheet_names = workbook.sheet_names()
+
+        self.worksheet_data = {}
+
+        for worksheet_name in self.worksheet_names:
+            worksheet = workbook.sheet_by_name(worksheet_name)
+            self.worksheet_data[worksheet_name] = \
+                get_worksheet_as_list(worksheet)
+
+        return
+
+    def get_var(
+        self,
+        worksheet_name,
+        flag,
+        offset=[0, 0],
+        post_process=None,
+    ):
+        """Find the cell in the first column of worksheet_name that (partially)
+        matches 'flag', and return the nearby cell value which is 'offset'
+        from this by [nrow, ncol]
+
+        Optionally, post_process is a function applied to the cell value
+        before returning. This can be useful to e.g. coerce integers to
+        bool
+
+        """
+
+        sheet_data = self.worksheet_data[worksheet_name]
+
+        matching_row = None
+
+        # Loop through the rows and look for a match in the first column
+        for i in range(len(sheet_data)):
+            if flag in sheet_data[i][0]:
+                matching_row = i
+                break
+
+        if matching_row is None:
+            msg = 'Could not find a cell in the first column of worksheet *' \
+                + worksheet_name + '* matching with *' + flag + '*'
+            raise Exception(msg)
+
+        # Check that offset[0] would not exceed the number of rows
+
+        if matching_row + offset[0] >= len(sheet_data):
+            return ""
+
+        # If the col offset exceeds the data length, return a blank
+        if offset[1] < len(sheet_data[i]):
+            output = sheet_data[i + offset[0]][offset[1]]
+            if post_process is not None:
+                output = post_process(output)
+            return output
+        else:
+            return ''
+
+        return
+
+    def get_paired_list(
+        self,
+        worksheet_name,
+        flag,
+        offset=[0, 0],
+        post_process = None,
+    ):
+        """Find the cell in the first column of worksheet_name that (partially)
+        matches 'flag', and read data from 2 rows of the spreadsheet,
+        starting at a location offset from the matching cell by
+        offset = [nrow, ncol].
+
+        The data is returned as a paired list:
+        [ [r00, r10], [r01, r11], [r02, r12]]
+
+        """
+
+        paired_list = []
+
+        counter = -1
+        while True:
+            counter = counter + 1
+            cell_x = offset[0]
+            cell_y = counter + offset[1]
+            next_row1 = self.get_var(worksheet_name, flag, [cell_x, cell_y],
+                                     post_process)
+            if next_row1 == '':
+
+                # We are at the end of the data
+
+                break
+            else:
+                cell_x = offset[0] + 1
+                next_row2 = self.get_var(worksheet_name, flag, 
+                                         [cell_x, cell_y],
+                                         post_process)
+
+            paired_list.append([next_row1, next_row2])
+
+        return paired_list
+
+    def get_list(
+        self,
+        worksheet_name,
+        flag,
+        offset=[0, 0],
+        post_process=None,
+    ):
+        """Find the cell in the first column of worksheet_name that (partially)
+        matches 'flag', and read data from 1 row of the spreadsheet,
+        starting at a location offset from the matching cell by
+        offset = [nrow, ncol].
+
+        The data is returned as a list:
+        [ r00, r01, r02, r03, ...]
+
+        """
+
+        out_list = []
+
+        counter = -1
+        while True:
+            counter = counter + 1
+            cell_x = offset[0]
+            cell_y = counter + offset[1]
+            next_row1 = self.get_var(worksheet_name, flag, 
+                                     [cell_x, cell_y], 
+                                     post_process)
+            if next_row1 == '':
+
+                # We are at the end of the data
+
+                break
+            out_list.append(next_row1)
+
+        return out_list
+
+    def get_subtable(
+        self,
+        worksheet_name,
+        flag,
+        offset=[0, 0],
+        post_process=None
+    ):
+        """Find the cell in the first column of worksheet_name that (partially)
+           matches 'flag', and read data from a nearby 'sub-table' in the
+           spreadsheet as a list of lists.
+
+           The 'sub-table' starts at a location
+           offset fom the matching cell by offset = [nrow, ncol], 
+           and finishes when a blank row of the sub-table is found
+
+           Rows can have different numbers of non-empty columns, but each row
+           ends when its first empty cell is found
+        """
+        
+        out_list = [ ]
+
+        counter = -1
+        while True: 
+            counter = counter + 1
+            cell_x = offset[0] + counter
+            cell_y = offset[1]
+            next_row = self.get_list(worksheet_name, flag, 
+                                     [cell_x, cell_y],
+                                     post_process)
+            if next_row == []:
+
+                break
+            out_list.append(next_row)
+
+        return out_list
+
+# #############################################################################    
+#
+# END OF CLASS 
+#           
+# #############################################################################
+
+
+# Functions below here
+
+def get_worksheet_as_list(worksheet):
+    """Read a particular xls (xlsx) file worksheet as a list of lists
+
+    [
+    [cell 0 0, cell 0 1, cell 0 2, ...],
+    [cell 1 0, cell 1 1, cell 1 2, ...],
+    ...
+    ]
+
+    """
+
+    worksheet_data = []
+
+    for i in range(worksheet.nrows):
+        worksheet_data.append([])
+        for j in range(worksheet.ncols):
+            worksheet_data[i].extend([worksheet.cell_value(i, j)])
+
+    return worksheet_data
+
+
+def empty_string_to_none(string):
+    """Given a string, return None if string == "", and the string value
+       (with type str) otherwise
+    """
+
+    if string == '':
+        return None
+    else:
+        return str(string)
+
+
+def string_or_float(value):
+    """ If value is unicode, return a value with type string.
+        If value is a number, return it with type float
+        If value is anything else, return value
+    """
+    if type(value) == unicode:
+        return str(value)
+    elif isinstance(value, (int, long, float, complex)):
+        return float(value)
+    else:
+        return value
+
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/prepare_data.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/prepare_data.py	(revision 9341)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/prepare_data.py	(revision 9341)
@@ -0,0 +1,243 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+"""
+
+THE USER MUST EDIT THIS FILE TO DEFINE THE MESH / RESOLUTION / FILE LOCATIONS
+ETC
+
+Setup base data for ANUGA run
+
+"""
+import glob
+import os
+import time
+import shutil
+import hashlib
+import json
+import numpy
+from anuga.utilities import spatialInputUtil as su
+from anuga_parallel import myid, barrier, send, receive, numprocs
+
+# Local modules
+import user_functions
+from parse_input_data import ProjectData
+
+class PrepareData(ProjectData):
+
+    def __init__(self, filename):
+        """Parse the input data then process it for ANUGA
+             
+        """
+        # Get the 'raw' data        
+        ProjectData.__init__(self, filename)
+
+        # Read files / pre-process
+        self.process_project_data()
+
+        # Make required directories
+        self.setup_directory_structure()
+
+        # Export spatial data to some text
+        self.export_spatial_data_as_text()
+
+    def process_project_data(self):
+        """Process the input data ready for ANUGA
+
+        """
+
+        # We can either use interior regions, or breaklines
+
+        if not self.interior_regions_data == []:
+            assert self.pt_areas is None, \
+                'Cannot define both ptAreas and non-empty interior regions'
+
+        bounding_polygon_and_tags = \
+            user_functions.read_boundary_tags_line_shapefile(
+                self.bounding_polygon_and_tags_file,
+                self.boundary_tags_attribute_name)
+        self.bounding_polygon = bounding_polygon_and_tags[0]
+        self.boundary_tags = bounding_polygon_and_tags[1]
+
+        self.breaklines = su.readListOfBreakLines(self.breakline_files)
+        (self.riverwalls, self.riverwall_par) = \
+            su.readListOfRiverWalls(self.riverwall_csv_files)
+
+        if self.pt_areas is not None:
+            self.region_point_areas = su.readRegionPtAreas(
+                self.pt_areas, convert_length_to_area=True)
+        else:
+            self.region_point_areas = None
+
+        # Hack to override resolution
+        # region_point_areas=[ region_point_areas[i][0:2]+[150*150*0.5] for i in \
+        #                      range(len(region_point_areas))]
+
+        # Redefine interior_regions to contain the polygon data + resolutions
+
+        self.interior_regions = [[su.read_polygon(ir[0]), ir[1]] for ir in
+                                  self.interior_regions_data]
+
+        # Deal with intersections in the bounding polygon / breaklines / riverwalls
+
+        (self.bounding_polygon, self.breaklines, self.riverwalls) = \
+            su.add_intersections_to_domain_features(
+                self.bounding_polygon, 
+                self.breaklines, 
+                self.riverwalls,
+                point_movement_threshold=\
+                    self.break_line_intersect_point_movement_threshold,
+                verbose=True)
+
+        # Here we make a unique ID based on the all the mesh geometry inputs
+        # This tells us if we need to regenerate partitions, or use old ones
+
+        self.mesh_id_hash = hashlib.md5(json.dumps([
+            self.bounding_polygon,
+            self.interior_regions,
+            self.riverwalls,
+            self.breaklines,
+            self.region_point_areas,
+            self.default_res,
+            self.boundary_tags,
+        ])).hexdigest()
+
+        # Fix the output tif bounding polygon
+
+        if self.output_tif_bounding_polygon is None:
+            self.output_tif_bounding_polygon = self.bounding_polygon
+        else:
+            self.output_tif_bounding_polygon = \
+                su.read_polygon(self.output_tif_bounding_polygon)
+
+        # Make proj4string from projection information
+        #
+
+        if type(self.projection_information) is int:
+
+            # projection_information describes a UTM zone
+            # e.g. '+units=m +ellps=WGS84 +zone=47 +south=False +proj=utm '
+
+            if self.projection_information < 0:
+                self.proj4string = '+proj=utm +zone=' \
+                    + str(abs(self.projection_information)) \
+                    + ' +south +datum=WGS84 +units=m +no_defs'
+            else:
+                self.proj4string = '+proj=utm +zone=' \
+                    + str(self.projection_information) \
+                    + ' +datum=WGS84 +units=m +no_defs'
+        elif type(self.projection_information) is str:
+            self.proj4string = self.projection_information
+        else:
+            msg = 'Invalid projection information ' + \
+                ' --  must be a proj4string, or an integer' + \
+                ' defining a UTM zone [positive for northern hemisphere,' + \
+                ' negative for southern hemisphere]'
+            raise Exception(msg)
+
+        # Set up directories etc
+
+        if myid == 0:
+            runtime = time.strftime('%Y%m%d_%H%M%S', time.localtime())
+            for i in range(1, numprocs):
+                send(runtime, i)
+        else:
+            runtime = receive(0)
+        barrier()
+
+        self.output_dir = self.output_basedir + 'RUN_' + str(runtime) +\
+                          '_' + self.scenario
+        self.partition_basedir = 'PARTITIONS/'
+        self.partition_dir = self.partition_basedir + 'Mesh_' +\
+                             str(self.mesh_id_hash)
+        self.meshname = self.output_dir + '/mesh.tsh'
+
+
+    def setup_directory_structure(self):
+        """Make output directories, send important data to files
+
+        This has no arguments, and returns nothing, but has lots of important
+        side effects
+        """
+
+        try:
+            os.mkdir(self.output_basedir)
+        except:
+            pass
+
+        # Make the output directory
+
+        try:
+            os.mkdir(self.output_dir)
+        except:
+            pass
+
+        # Make the partitions directory (2 steps as it is recursive)
+
+        try:
+            os.mkdir(self.partition_basedir)
+        except:
+            pass
+
+        try:
+            os.mkdir(self.partition_dir)
+        except:
+            pass
+
+        # Make the spatialTxt directory
+
+        try:
+            os.mkdir(self.spatial_text_output_dir)
+        except:
+            pass
+
+    def export_spatial_data_as_text(self):
+        """Export the imported spatial data as text
+           This is useful to check that the data was imported correctly
+
+        """
+        # #################################################################
+        # Export breaklines to txt for QC
+
+        def sav_lines(lines, filename='check.txt'):
+            """quick convenience function
+            """
+            new_a = numpy.array(lines)
+            numpy.savetxt(filename, new_a)
+            return
+
+        if myid == 0:
+            if self.interior_regions is not []:
+                for i in range(len(self.interior_regions)):
+                    sav_lines(self.interior_regions[i][0],
+                              filename=self.spatial_text_output_dir
+                              + '/interior_region_' + str(i) + '.txt')
+            if self.breaklines is not {}:
+                for i in range(len(self.breaklines)):
+                    sav_lines(self.breaklines.values()[i],
+                              filename=self.spatial_text_output_dir +\
+                                         '/breakline_' + str(i) + '.txt')
+            if self.riverwalls is not {}:
+                for i in range(len(self.riverwalls)):
+                    sav_lines(self.riverwalls.values()[i],
+                              filename=self.spatial_text_output_dir +\
+                                 '/riverwall_' + str(i) + '.txt')
+
+            sav_lines(self.bounding_polygon, 
+                      filename=self.spatial_text_output_dir
+                      + '/boundingPoly.txt')
+
+            if self.region_point_areas is not None:
+                sav_lines(self.region_point_areas,
+                          filename=self.spatial_text_output_dir
+                          + '/regionPtAreas.txt')
+
+            # Copy all the txt representations of these files to the output
+            # directory
+
+            shutil.copytree(self.spatial_text_output_dir, 
+                            self.output_dir +\
+                             '/' + self.spatial_text_output_dir)
+        else:
+            pass
+
+        return
Index: unk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/project.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/project.py	(revision 9340)
+++ 	(revision )
@@ -1,402 +1,0 @@
-#!/usr/bin/python
-# -*- coding: utf-8 -*-
-"""
-
-THE USER MUST EDIT THIS FILE TO DEFINE THE MESH / RESOLUTION / FILE LOCATIONS
-ETC
-
-Setup base data for ANUGA run
-
-"""
-
-
-import glob
-import os
-import time
-import shutil
-import hashlib
-import json
-import numpy
-from anuga.utilities import spatialInputUtil as su
-from anuga_parallel import myid, barrier, send, receive, numprocs
-
-# We need the boundary tags info to correctly cache the mesh
-# from setup_boundary_conditions import boundary_info
-
-import user_functions
-
-# @@@@@@@@@@@@@@@@@#
-# START EDITING HERE
-# @@@@@@@@@@@@@@@@@#
-
-#####################################################################
-#
-# SCENARIO SETUP
-#
-#####################################################################
-
-# Name for scenario
-#
-
-scenario = 'tsunami1_template'
-
-# Total model runtime in seconds (starts at time=0.)
-#
-
-finaltime = 10.  # 0.05 * 3600.0
-
-# Model time between output file timesteps in seconds
-#
-
-yieldstep = 20.0
-
-# Spatial Information
-#
-# Give UTM Zone as an integer (negative for Southern Hemisphere)
-# or a proj4string defining the coordinate system (which is more flexible)
-
-projection_information = 47
-
-# projection_information='+proj=utm +zone=47 +datum=WGS84 +units=m +no_defs'
-
-#####################################################################
-#
-# INPUT DATA DEFINING MESH
-#
-#####################################################################
-
-# Bounding polygon with default mesh resolution (as a max triangle area in m^2)
-#
-# Filename MUST be a .shp file, and contain a set of lines making up the
-# boundary polygon The lines each have an attribute giving the boundary tag, in
-# the attribute table column named (the value of boundary_tags_attribute_name)
-#
-# The starting/ending points on each line MUST EXACTLY AGREE so the whole
-# thing makes a polygon.
-#
-# Make it in GIS by first making the bounding polygon, then splitting it at
-# the appropriate locations, then adding the boundary tag attribute.
-#
-# bounding_polygon_file='../anuga/polygons/bounding_polygon.csv'
-
-bounding_polygon_and_tags_file = \
-    '../anuga/boundaries/Patong_bounding_polygon/Patong_bounding_polygon.shp'
-boundary_tags_attribute_name = 'bndryTag'
-
-default_res = 0.5 * 500 ** 2
-
-# Interior regions with maximum triangle area resolutions (m^2)
-#
-# Filenames can be .shp or .csv
-#
-# If this is an empty list, we can instead use breaklines + ptAreas defined
-# spatially
-#
-
-interior_regions_data = [
-    ['../anuga/polygons/inundation_area.csv', 75],
-    ['../anuga/polygons/aoi.csv', 200],
-    ['../anuga/polygons/aos.csv', 900],
-    ['../anuga/polygons/sw_coast_original.csv', 7000],
-    ['../anuga/polygons/building_main.csv', 20],
-    ['../anuga/polygons/building_main_small.csv', 20],
-    ['../anuga/polygons/building_main_south.csv', 19],
-    ['../anuga/polygons/building_saddle.csv', 20],
-]
-
-# BreakLines.
-#
-# List of filenames (can be .shp and/or .csv)
-#
-# breakline_files=glob.glob('../../GIS/Mesh/BreakLines/*.shp')
-
-breakline_files = []
-
-# Riverwalls
-#
-# Filenames MUST all be .csv with at least x,y,z of riverwall, and and optional
-# header describing the riverwallPar
-#
-# riverwall_csv_files=glob.glob('../../GIS/Mesh/RiverWalls/*.csv')
-
-riverwall_csv_files = []
-
-# If breaklines/riverwalls intersect with < this distance between the
-# intersection and an existing point, then the existing point is moved to the
-# intersection. Otherwise the intersection point is added
-
-break_line_intersect_point_movement_threshold = 10.
-
-# Get mesh resolutions = region_point_areas [alternative to interior_polygons +
-# res] ptAreas='../../GIS/Mesh/RegionResolutions/RegionResolutions.shp'
-
-pt_areas = None
-
-#####################################################################
-#
-# Numerics
-#
-#####################################################################
-
-# ANUGA's flow algorithm
-#
-
-flow_algorithm = 'DE1'
-
-# Local extrapolation and flux updating
-#
-# Only set this to True if you know what you are doing
-#
-
-use_local_extrapolation_and_flux_updating = False
-
-#########################################################################
-#
-# OUTPUT INFORMATION
-#
-#########################################################################
-
-# Output tif details
-#
-# pixle size (m) for output tifs
-
-output_tif_cellsize = 250.
-
-# Polygon for output tif [either a filename, or if None the
-# bounding_polygon is used]
-
-output_tif_bounding_polygon = None
-
-# Storage of max quantities
-#
-# Stage/depth/speed/speedxDepth maxima are stored separately to the SWW file
-# Need to determine whether they are updated every time step or less frequently
-# (less frequently might save computational time ?? )
-# Can also choose to ignore an initial time period
-
-# 1 = Every time step, 2 = every second timestep, etc
-max_quantity_update_frequency = 1
-# Set this to e.g. 3600.*2 to ignore the first 2 hours
-max_quantity_collection_start_time = 0.
-
-# sww file-size options
-#
-# If store_vertices_uniquely=True, then SWW files are larger and vertices are
-# stored uniquely for every triangle. Otherwise, vertex values are averaged
-# before storage, leading to smaller file sizes but a less accurate
-# representation of what the algorithm is doing
-
-store_vertices_uniquely = False
-
-# Checkpointing
-#
-# Need to restructure to allow checkpointing
-# checkpoint_time=30
-
-# @@@@@@@@@@@@@@@@@#
-# STOP EDITING HERE
-# @@@@@@@@@@@@@@@@@#
-
-###########################################################################
-#
-# BELOW HERE, THINGS ARE DESIGNED NOT TO BE CHANGED
-#
-###########################################################################
-
-# We can either use interior regions, or breaklines
-
-if not interior_regions_data == []:
-    assert pt_areas is None, \
-        'Cannot define both ptAreas and non-empty interior regions'
-
-# bounding_polygon=su.read_polygon(bounding_polygon_file)
-
-bounding_polygon_and_tags = \
-    user_functions.read_boundary_tags_line_shapefile(
-        bounding_polygon_and_tags_file,
-        boundary_tags_attribute_name)
-bounding_polygon = bounding_polygon_and_tags[0]
-boundary_tags = bounding_polygon_and_tags[1]
-
-breaklines = su.readListOfBreakLines(breakline_files)
-(riverwalls, riverwall_par) = \
-    su.readListOfRiverWalls(riverwall_csv_files)  # {} # Dictionary
-
-if pt_areas is not None:
-    region_point_areas = su.readRegionPtAreas(
-        pt_areas,
-        convert_length_to_area=True)
-else:
-    region_point_areas = None
-
-# Hack to override resolution
-# region_point_areas=[ region_point_areas[i][0:2]+[150*150*0.5] for i in \
-#                      range(len(region_point_areas))]
-
-# Redefine interior_regions to contain the polygon data + resolutions
-
-interior_regions = [[su.read_polygon(ir[0]), ir[1]] for ir in
-                    interior_regions_data]
-
-# Directory where sww files are written
-
-output_basedir = '../anuga/outputs/'
-
-# Directory where we store .txt versions of processed spatial inputs
-
-spatial_txt_output_dir = 'SpatialTxt/'
-
-# Deal with intersections in the bounding polygon / breaklines / riverwalls
-
-(bounding_polygon, breaklines, riverwalls) = \
-    su.add_intersections_to_domain_features(
-        bounding_polygon, breaklines, riverwalls,
-        point_movement_threshold=break_line_intersect_point_movement_threshold,
-        verbose=True)
-
-# Here we make a unique ID based on the all the mesh geometry inputs
-# This tells us if we need to regenerate partitions, or use old ones
-
-mesh_id_hash = hashlib.md5(json.dumps([
-    bounding_polygon,
-    interior_regions,
-    riverwalls,
-    breaklines,
-    region_point_areas,
-    default_res,
-    boundary_tags,
-])).hexdigest()
-
-# Fix the output tif bounding polygon
-
-if output_tif_bounding_polygon is None:
-    output_tif_bounding_polygon = bounding_polygon
-else:
-    output_tif_bounding_polygon = \
-        su.read_polygon(output_tif_bounding_polygon)
-
-# Make proj4string from projection information
-#
-
-if type(projection_information) is int:
-
-    # projection_information describes a UTM zone
-    # e.g. '+units=m +ellps=WGS84 +zone=47 +south=False +proj=utm '
-
-    if projection_information < 0:
-        proj4string = '+proj=utm +zone=' \
-            + str(abs(projection_information)) \
-            + ' +south +datum=WGS84 +units=m +no_defs'
-    else:
-        proj4string = '+proj=utm +zone=' + str(projection_information) \
-            + ' +datum=WGS84 +units=m +no_defs'
-elif type(projection_information) is str:
-    proj4string = projection_information
-else:
-    msg = 'Invalid projection information ' + \
-        ' --  must be a proj4string, or an integer' + \
-        ' defining a UTM zone [positive for northern hemisphere,' + \
-        ' negative for southern hemisphere]'
-    raise Exception(msg)
-
-# Set up directories etc
-
-if myid == 0:
-    runtime = time.strftime('%Y%m%d_%H%M%S', time.localtime())
-    for i in range(1, numprocs):
-        send(runtime, i)
-else:
-    runtime = receive(0)
-barrier()
-
-output_dir = output_basedir + 'RUN_' + str(runtime) + '_' + scenario
-partition_basedir = 'PARTITIONS/'
-partition_dir = partition_basedir + 'Mesh_' + str(mesh_id_hash)
-meshname = output_dir + '/mesh.tsh'
-
-
-def setup_directory_structure():
-    """Make output directories, send important data to files
-
-    This has no arguments, and returns nothing, but has lots of important side
-    effects
-    """
-
-    # if(myid==0):
-    # Make the base directory if it does not exist
-
-    try:
-        os.mkdir(output_basedir)
-    except:
-        pass
-
-    # Make the output directory
-
-    try:
-        os.mkdir(output_dir)
-    except:
-        pass
-
-    # Make the partitions directory (2 steps as it is recursive)
-
-    try:
-        os.mkdir(partition_basedir)
-    except:
-        pass
-
-    try:
-        os.mkdir(partition_dir)
-    except:
-        pass
-
-    # Make the spatialTxt directory
-
-    try:
-        os.mkdir(spatial_txt_output_dir)
-    except:
-        pass
-
-    # #########################################################################
-    # Export breaklines to txt for QC
-
-    def sav_lines(lines, filename='check.txt'):
-        """quick convenience function
-        """
-        new_a = numpy.array(lines)
-        numpy.savetxt(filename, new_a)
-        return
-
-    if myid == 0:
-        if interior_regions is not []:
-            for i in range(len(interior_regions)):
-                sav_lines(interior_regions[i][0],
-                          filename=spatial_txt_output_dir
-                          + 'interior_region_' + str(i) + '.txt')
-        if breaklines is not {}:
-            for i in range(len(breaklines)):
-                sav_lines(breaklines.values()[i],
-                          filename=spatial_txt_output_dir + 'breakline_'
-                          + str(i) + '.txt')
-        if riverwalls is not {}:
-            for i in range(len(riverwalls)):
-                sav_lines(riverwalls.values()[i],
-                          filename=spatial_txt_output_dir + 'riverwall_'
-                          + str(i) + '.txt')
-
-        sav_lines(bounding_polygon, filename=spatial_txt_output_dir
-                  + 'boundingPoly.txt')
-
-        if region_point_areas is not None:
-            sav_lines(region_point_areas,
-                      filename=spatial_txt_output_dir
-                      + 'regionPtAreas.txt')
-
-        # Copy all the txt representations of these files to the output
-        # directory
-
-        shutil.copytree(spatial_txt_output_dir, output_dir + '/'
-                        + spatial_txt_output_dir)
-    else:
-        pass
-
-    return
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/run_model.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/run_model.py	(revision 9340)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/run_model.py	(revision 9341)
@@ -56,5 +56,6 @@
 # Get key data for the simulation
 
-import project
+from prepare_data import PrepareData
+project = PrepareData('ANUGA_setup.xls')
 
 ###############################################################################
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/setup_boundary_conditions.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/setup_boundary_conditions.py	(revision 9340)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/setup_boundary_conditions.py	(revision 9341)
@@ -13,4 +13,5 @@
 import anuga
 import anuga.utilities.spatialInputUtil as su
+import anuga.shallow_water.boundaries as asb
 
 
@@ -27,37 +28,46 @@
     """
 
-# @@@@@@@@@@@@@@@@@#
-# START EDITING HERE
-# @@@@@@@@@@@@@@@@@#
+    # Dict to hold { boundary_tag: boundary_condition_function }
+    boundary_tags_and_conditions = {}
 
-    # Example:
-    # "Standard" boundary conditions
+    for i in range(len(project.boundary_data)):
+        bd = project.boundary_data[i]
+        boundary_tag = bd[0]
+        boundary_condition_type = bd[1]
 
-    br = anuga.Reflective_boundary(domain)
+        if boundary_condition_type == 'Reflective':
+            # Simple reflective boundary
+            boundary_tags_and_conditions[boundary_tag] = \
+                anuga.Reflective_boundary(domain)
+        elif boundary_condition_type == 'Stage':
+            # Here we read a timeseries, offset the starttime, and then pass
+            # that function as a set stage transmissive n momentum ... boundary
+            boundary_data_file = bd[2]
+            start_time = bd[3]
+            # Read data
+            boundary_data = scipy.genfromtxt(
+                boundary_data_file,
+                delimiter=',',
+                skip_header=1)
+            # Make start time = 0
+            boundary_data[:,0]-=start_time
+            # Make interpolation function
+            stage_time_fun = scipy.interpolate.interp1d(
+                boundary_data[:,0],
+                boundary_data[:,1])
+            # Make boundary condition 
+            boundary_function = \
+                asb.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(
+                    domain,
+                    stage_time_fun)
+            
+            boundary_tags_and_conditions[boundary_tag] = \
+                boundary_function
 
-    # Bs = anuga.Transmissive_stage_zero_momentum_boundary(domain) # Be
-    # careful with this one
+        else:
+            msg = 'Boundary condition type ' + boundary_condition_type +\
+                  ' for tag ' + boundary_tag + ' is not implemented'
 
-    # Example:
-    # Transmissive with stage being set by a function 'stage_function'
-
-    def stage_function(t):
-        return min(max(t - 60., 0.) / 60., 1.)
-
-    import anuga.shallow_water.boundaries as asb
-    b_t_stage_function = \
-        asb.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(
-            domain,
-            stage_function)
-
-    # Here we associate each tag in boundary_info with its boundary condition
-    # defined above
-
-    boundary_tags_and_conditions = {'Ocean': b_t_stage_function,
-                                    'Br': br}
-
-# @@@@@@@@@@@@@@@@@#
-# STOP EDITING HERE
-# @@@@@@@@@@@@@@@@@#
+                        
 
     # Check that all tags in boundary_info have been set
@@ -70,4 +80,6 @@
             raise Exception(msg)
 
+    # Set the boundary 
+
     domain.set_boundary(boundary_tags_and_conditions)
 
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/setup_initial_conditions.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/setup_initial_conditions.py	(revision 9340)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/setup_initial_conditions.py	(revision 9341)
@@ -69,133 +69,4 @@
     #
 
-    # @@@@@@@@@@@@@@@@@#
-    # START EDITING HERE
-    # @@@@@@@@@@@@@@@@@#
-
-    # #########################################################################
-    #
-    # PRELIMINARY DEFINITIONS
-    #
-    # #########################################################################
-
-    # If the earthquake slip is given then convert to deformation here
-    # So that we can add it to both elevation AND stage
-
-    make_earthquake_deformation_from_source_model = False
-
-    if make_earthquake_deformation_from_source_model:
-        source_model = '../anuga/source/XXXXX.csv'
-
-        # Make an earthquake source model
-
-        eq_function = okada_tsunami.earthquake_source(
-            source=None,
-            filename=source_model,
-            domain=domain,
-            lon_lat_degrees=True,
-            lon_before_lat=False,
-            proj4string=project.proj4string,
-            verbose=False,
-        )
-    else:
-
-        # Make a function that is zero everywhere
-        # If we pass this to add_quantity, nothing should happen
-
-        def eq_function(x, y):
-            return x * 0.
-
-    # #########################################################################
-    #
-    # Set elevation.
-    #
-    # #########################################################################
-
-    # The 'Patong 10m' dataset
-    elevation_patong_10m = '../anuga/topographies/patong_10m_grid_msl_Project.txt'
-
-    # The 'saddle 10m' dataset
-    elevation_saddle_10m = \
-        '../anuga/topographies/patong_10m_small_grid_for_anuga_sub_Project.txt'
-
-    # The '1s grid' dataset
-    elevation_patong_1s = \
-            '../anuga/topographies/patong_bay_1s_grid_Project.txt'
-
-    # The background data
-    elevation_default = \
-            '../anuga/topographies/andaman_3s_grid_Clip2_Project.txt'
-
-    # Set_quantity data
-    elevation_data = [['../anuga/polygons/patong_10m.txt',
-                       elevation_patong_10m],
-                      ['../anuga/polygons/saddle_10m.txt',
-                       elevation_saddle_10m],
-                      ['../anuga/polygons/patong_1s.txt',
-                       elevation_patong_1s], 
-                      ['All', elevation_default]]
-
-    # Define a function for add_quantity in the same way
-    # [ ['All', 0.] ] will do nothing
-
-    elevation_additions = [['All', eq_function]]
-
-    # #############################################################################
-    #
-    # Set friction
-    #
-    # #############################################################################
-
-    default_friction = 0.025
-
-    # Set_quantity data
-
-    friction_data = [['All', default_friction]]
-
-    # Add_quantity data
-
-    friction_additions = [['All', 0.]]
-
-    # ##############################################################################
-    #
-    # Set stage
-    #
-    # ##############################################################################
-
-    # Set_quantity data
-
-    stage_data = [['All', 0.]]
-
-    # Add_quantity data
-
-    stage_additions = [['All', eq_function]]
-
-    # ##############################################################################
-    #
-    # Set depth integrated velocity (called "momentum" in ANUGA although that
-    # is incorrect")
-    #
-    # ##############################################################################
-
-    # Set_quantity data
-
-    xmomentum_data = [['All', 0.]]
-
-    # Add_quantity data
-
-    xmomentum_additions = [['All', 0.]]
-
-    # Set_quantity data
-
-    ymomentum_data = [['All', 0.]]
-
-    # Add_quantity data
-
-    ymomentum_additions = [['All', 0.]]
-
-    # @@@@@@@@@@@@@@@@@#
-    # STOP EDITING HERE
-    # @@@@@@@@@@@@@@@@@#
-
     # ###############################################################################
     #
@@ -208,43 +79,43 @@
 
     elevation_function = \
-        qs.composite_quantity_setting_function(elevation_data, domain)
+        qs.composite_quantity_setting_function(project.elevation_data, domain)
     domain.set_quantity('elevation', elevation_function,
                         location='vertices')
     elevation_addition_function = \
-        qs.composite_quantity_setting_function(elevation_additions,
+        qs.composite_quantity_setting_function(project.elevation_additions,
                                                domain)
     domain.add_quantity('elevation', elevation_addition_function)
 
     friction_function = \
-        qs.composite_quantity_setting_function(friction_data, domain)
+        qs.composite_quantity_setting_function(project.friction_data, domain)
     domain.set_quantity('friction', friction_function,
                         location='centroids')
     friction_addition_function = \
-        qs.composite_quantity_setting_function(friction_additions,
+        qs.composite_quantity_setting_function(project.friction_additions,
                                                domain)
     domain.add_quantity('friction', friction_addition_function)
 
-    stage_function = qs.composite_quantity_setting_function(stage_data,
+    stage_function = qs.composite_quantity_setting_function(project.stage_data,
                                                             domain)
     domain.set_quantity('stage', stage_function, location='centroids')
     stage_addition_function = \
-        qs.composite_quantity_setting_function(stage_additions, domain)
+        qs.composite_quantity_setting_function(project.stage_additions, domain)
     domain.add_quantity('stage', stage_addition_function)
 
     xmomentum_function = \
-        qs.composite_quantity_setting_function(xmomentum_data, domain)
+        qs.composite_quantity_setting_function(project.xmomentum_data, domain)
     domain.set_quantity('xmomentum', xmomentum_function,
                         location='centroids')
     xmomentum_addition_function = \
-        qs.composite_quantity_setting_function(xmomentum_additions,
+        qs.composite_quantity_setting_function(project.xmomentum_additions,
                                                domain)
     domain.add_quantity('xmomentum', xmomentum_addition_function)
 
     ymomentum_function = \
-        qs.composite_quantity_setting_function(ymomentum_data, domain)
+        qs.composite_quantity_setting_function(project.ymomentum_data, domain)
     domain.set_quantity('ymomentum', ymomentum_function,
                         location='centroids')
     ymomentum_addition_function = \
-        qs.composite_quantity_setting_function(ymomentum_additions,
+        qs.composite_quantity_setting_function(project.ymomentum_additions,
                                                domain)
     domain.add_quantity('ymomentum', ymomentum_addition_function)
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/setup_mesh.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/setup_mesh.py	(revision 9340)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/setup_mesh.py	(revision 9341)
@@ -116,8 +116,4 @@
     if myid == 0:
 
-        # Make sure files/folders to store everything are there
-
-        project.setup_directory_structure()
-
         # Copy code files to output dir. This should be 'nearly the first'
         # thing we do
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/user_functions.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/user_functions.py	(revision 9340)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/user_functions.py	(revision 9341)
@@ -94,8 +94,4 @@
         {boundary_line_and_tags[0][1]: range(len(bounding_polygon) - 1)}
 
-    #
-    # import pdb
-    # pdb.set_trace()
-
     for i in range(boundary_segnum - 1):
 
