Index: /anuga_work/production/busselton/standardised_version/build_busselton.py
===================================================================
--- /anuga_work/production/busselton/standardised_version/build_busselton.py	(revision 6253)
+++ /anuga_work/production/busselton/standardised_version/build_busselton.py	(revision 6253)
@@ -0,0 +1,125 @@
+"""
+Script for building 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_busselton.py is reliant on the output of this script.
+
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+from os import sep
+from os.path import dirname, basename
+from os import mkdir, access, F_OK
+from shutil import copy
+import time
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from anuga.geospatial_data.geospatial_data import *
+from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
+
+# Application specific imports
+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_time_dir,__file__, 
+               dirname(project.__file__)+sep+ project.__name__+'.py' )
+
+start_screen_catcher(project.output_build_time_dir)
+
+print 'USER:    ', project.user
+
+#-------------------------------------------------------------------------------
+# 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.poly_all", project.poly_all
+print "project.combined_dir_name", project.combined_dir_name
+
+# input elevation directory filenames
+onshore_in_dir_name = project.onshore_in_dir_name
+coast_in_dir_name = project.coast_in_dir_name
+coast_in_dir_name1 = project.coast_in_dir_name1
+offshore_in_dir_name = project.offshore_in_dir_name
+offshore_in_dir_name1 = project.offshore_in_dir_name1
+offshore_in_dir_name2 = project.offshore_in_dir_name2
+offshore_in_dir_name3 = project.offshore_in_dir_name3
+offshore_in_dir_name4 = project.offshore_in_dir_name4
+offshore_in_dir_name5 = project.offshore_in_dir_name5
+
+# output elevation directory filenames
+onshore_dir_name = project.onshore_dir_name
+coast_dir_name = project.coast_dir_name
+coast_dir_name1 = project.coast_dir_name1
+offshore_dir_name = project.offshore_dir_name
+offshore_dir_name1 = project.offshore_dir_name1
+offshore_dir_name2 = project.offshore_dir_name2
+offshore_dir_name3 = project.offshore_dir_name3
+offshore_dir_name4 = project.offshore_dir_name4
+offshore_dir_name5 = project.offshore_dir_name5
+# creates DEM from asc data
+print "creates DEMs from asc data"
+convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
+convert_dem_from_ascii2netcdf(offshore_in_dir_name5, basename_out=offshore_dir_name5, use_cache=True, verbose=True)
+
+# creates pts file for onshore DEM
+print "creates pts file for onshore DEM"
+dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
+dem2pts(offshore_dir_name5 ,use_cache=True,verbose=True)
+
+# create onshore pts files
+print'create Geospatial data1 objects from topographies'
+G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
+
+# create coastal and offshore pts files
+print'create Geospatial data2 objects from topographies'
+G2 = Geospatial_data(file_name = coast_in_dir_name)
+print'create Geospatial data3 objects from topographies'
+G3 = Geospatial_data(file_name = coast_in_dir_name1)
+print'create Geospatial data4 objects from topographies'
+G_off = Geospatial_data(file_name = offshore_in_dir_name)
+print'create Geospatial data5 objects from topographies'
+G_off1 = Geospatial_data(file_name = offshore_in_dir_name1)
+print'create Geospatial data6 objects from topographies'
+G_off2 = Geospatial_data(file_name = offshore_in_dir_name2)
+print'create Geospatial data7 objects from topographies'
+G_off3 = Geospatial_data(file_name = offshore_in_dir_name3)
+print'create Geospatial data8 objects from topographies'
+G_off4 = Geospatial_data(file_name = offshore_in_dir_name4)
+print'create Geospatial data9 objects from topographies'
+G_off5 = Geospatial_data(file_name = offshore_dir_name5 + '.pts')
+
+
+#-------------------------------------------------------------------------------
+# Combine, clip and export dataset 
+#-------------------------------------------------------------------------------
+
+print'add all geospatial objects'
+G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4
+
+print'clip combined geospatial object by bounding polygon'
+G_clip = G.clip_outside(project.poly_aoi1)
+G_all = G_clip + G_off5
+G_clipped = G_all.clip(project.poly_all)
+
+print'export combined DEM file'
+if access(project.topographies_dir,F_OK) == 0:
+    mkdir (project.topographies_dir)
+G_clipped.export_points_file(project.combined_dir_name + '.pts')
+G_clipped.export_points_file(project.combined_dir_name + '.txt') #Use for comparision in ARC
+
Index: /anuga_work/production/busselton/standardised_version/project.py
===================================================================
--- /anuga_work/production/busselton/standardised_version/project.py	(revision 6253)
+++ /anuga_work/production/busselton/standardised_version/project.py	(revision 6253)
@@ -0,0 +1,274 @@
+"""Common filenames and locations for elevation, meshes and outputs.
+This script is the heart of all scripts in the folder
+"""
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser
+import sys
+from time import localtime, strftime, gmtime
+from anuga.utilities.polygon import read_polygon, plot_polygons, is_inside_polygon, number_mesh_triangles
+from anuga.utilities.system_tools import get_user_name, get_host_name
+from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+
+#------------------------------------------------------------------------------
+# Directory setup
+#------------------------------------------------------------------------------
+# Note: INUNDATIONHOME is the inundation directory, not the data directory.
+
+home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent diruser = get_user_name()
+muxhome = getenv('MUXHOME')
+user = get_user_name()
+host = get_host_name()
+
+# 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'
+
+#------------------------------------------------------------------------------
+# Initial Conditions
+#------------------------------------------------------------------------------
+
+# this section needs to be updated to reflect the modelled community.
+# Note, the user needs to set up the directory system accordingly
+state = 'western_australia'
+scenario_name = 'busselton'
+scenario = 'busselton_tsunami_scenario'
+
+# Model specific parameters. One or all can be changed each time the
+# run_scenario script is executed
+tide = 0              #0.6
+#event_number = 27255 # Java 9.3 worst case for Perth
+#event_number = 68693 # Sumatra 9.2
+event_number = 27283  # Java 9.3 original
+alpha = 0.1             # smoothing parameter for mesh
+friction=0.01           # manning's friction coefficient
+starttime=0             
+finaltime=80000         # final time for simulation
+
+setup='final'  # Final can be replaced with trial or basic.
+               # Either will result in a coarser mesh that will allow a
+               # faster, but less accurate, simulation.
+
+if setup =='trial':
+    print'trial'
+    res_factor=10
+    time_thinning=48
+    yieldstep=240
+if setup =='basic': 
+    print'basic'
+    res_factor=4
+    time_thinning=12
+    yieldstep=120
+if setup =='final': 
+    print'final'
+    res_factor=1
+    time_thinning=4
+    yieldstep=60
+
+#------------------------------------------------------------------------------
+# Revision numbers - for comparisons study
+#------------------------------------------------------------------------------
+rev_num = 'newExtent'
+#rev_num = '5449'
+#rev_num = '4695' # 2nd Sept 2007
+#rev_num = '4743' # 3nd Oct 2007
+#rev_num = '4777' # 1st Nov 2007
+#rev_num = '4874' # 3rd Dec 2007
+#rev_num = '4901' # 3rd Jan 2007
+#rev_num = '4990' # 5th Feb 2007
+#rev_num = '5103' # 3rd March 2007
+#rev_num = '5120' # 5th March 2007
+#rev_num = '5140' # 7th March 2007
+#rev_num = '5160' # 11th March 2007
+#rev_num = '5185' # 1st April 2007
+#rev_num = '5273' # 2nd May 2007
+
+
+#------------------------------------------------------------------------------
+# Output Filename
+#------------------------------------------------------------------------------
+# Important to distinguish each run - ensure str(user) is included!
+# Note, the user is free to include as many parameters as desired
+dir_comment='_'+setup+'_'+str(tide)+'_'+str(event_number)+'_'+ 'alpha' +str(alpha)+'_'+str(user)
+
+#------------------------------------------------------------------------------
+# Input Data
+#------------------------------------------------------------------------------
+
+# elevation data used in build_busselton.py
+# onshore data: format ascii grid with accompanying projection file
+onshore_name = 'busselton_v2_gda94_mga50' 
+# coastline: format x,y,elevation (with title)
+coast_name = 'Busselton_Contour0.txt'
+coast_name1 = 'Busselton_BeachSurvey.txt'
+# bathymetry: format x,y,elevation (with title)
+offshore_name = 'Busselton_NavyFinal.txt'
+offshore_name1 = 'Busselton_Chart.txt'
+offshore_name2 = 'Busselton_Digitised.txt'
+offshore_name3 = 'Busselton_250m.txt' # for areas that were heading to zero - 2005 Bathymetry grid
+offshore_name4 = 'DPI.txt' # for area within Bunbury 500 mesh less than zero generated from TIN
+offshore_name5 = 'topo_20m_buss_1km' # for area within Busselton 500 mesh less than zero generated from TIN
+
+# gauges - used in get_timeseries.py
+#gauge_name = scenario_name+'.txt'
+gauge_name = 'Gauges.csv'
+# buildings - used in run_building_inundation.py
+building = 'busselton_res_clip'
+
+
+# BOUNDING POLYGON - used in build_boundary.py and run_busselton.py respectively
+# NOTE: when files are put together the points must be in sequence - for ease go clockwise!
+# Check the run_busselton.py for boundary_tags
+# thinned ordering file from Hazard Map: format is index,latitude,longitude (with title)
+order_filename = 'thinned_boundary_ordering.csv'
+#landward bounding points
+landward = 'landward_bounding_polygon.csv'
+
+#------------------------------------------------------------------------------
+# Output Elevation Data
+#------------------------------------------------------------------------------
+# Output filename for elevation
+# this is a combination of all the data (utilisied in build_boundary)
+combined_name ='busselton_combined_elevation'
+combined_smaller_name = 'busselton_combined_elevation_smaller'
+
+#------------------------------------------------------------------------------
+# Directory Structure
+#------------------------------------------------------------------------------
+anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
+topographies_in_dir = home+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
+topographies_dir = anuga_dir+'topographies'+sep
+polygons_dir = anuga_dir+'polygons'+sep
+tide_dir = anuga_dir+'tide_data'+sep
+boundaries_dir = anuga_dir+'boundaries'+ sep 
+output_dir = anuga_dir+'outputs'+sep
+gauges_dir = anuga_dir+'gauges'+sep
+meshes_dir = anuga_dir+'meshes'+sep
+
+#------------------------------------------------------------------------------
+# Location of input and output data
+#------------------------------------------------------------------------------
+# where the input data sits
+onshore_in_dir_name = topographies_in_dir + onshore_name #topo
+coast_in_dir_name = topographies_in_dir + coast_name #coastline
+coast_in_dir_name1 = topographies_in_dir + coast_name1 #beach survey
+offshore_in_dir_name = topographies_in_dir + offshore_name #bathymetry
+offshore_in_dir_name1 = topographies_in_dir + offshore_name1 #bathymetry Charts
+offshore_in_dir_name2 = topographies_in_dir + offshore_name2 #Digitised Fairsheet
+offshore_in_dir_name3 = topographies_in_dir + offshore_name3 #250m
+offshore_in_dir_name4 = topographies_in_dir + offshore_name4 #Bunbury DPI
+offshore_in_dir_name5 = topographies_in_dir + offshore_name5 #Busselton Topo
+
+# where the output data sits
+onshore_dir_name = topographies_dir + onshore_name
+
+coast_dir_name = topographies_dir + coast_name
+coast_dir_name1 = topographies_dir + coast_name1
+
+offshore_dir_name = topographies_dir + offshore_name
+offshore_dir_name1 = topographies_dir + offshore_name1
+offshore_dir_name2 = topographies_dir + offshore_name2
+offshore_dir_name3 = topographies_dir + offshore_name3
+offshore_dir_name4 = topographies_dir + offshore_name4
+offshore_dir_name5 = topographies_dir + offshore_name5
+
+# where the combined elevation file sits
+combined_dir_name = topographies_dir + combined_name
+combined_smaller_name_dir = topographies_dir + combined_smaller_name
+
+# where the mesh sits (this is created during the run_busselton.py)
+meshes_dir_name = meshes_dir + scenario_name+'.msh'
+
+# where the boundary ordering files sit (this is used within build_boundary.py)
+order_filename_dir = boundaries_dir + order_filename
+
+# where the landward points of boundary extent sit (this is used within run_busselton.py)
+landward_dir = boundaries_dir + landward
+
+# where the event sts files sits (this is created during the build_boundary.py)
+boundaries_dir_event = boundaries_dir + str(event_number) + sep
+boundaries_dir_mux = muxhome
+
+# where the directory of the output filename sits
+output_build_time_dir = output_dir+build_time+dir_comment+sep   #used for build_busselton.py
+output_run_time_dir = output_dir+run_time+dir_comment+sep       #used for run_busselton.py
+output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
+
+#w here the directory of the gauges sit
+gauges_dir_name = gauges_dir + gauge_name       #used for get_timeseries.py
+building_in_dir_name = gauges_dir + building + '.csv'    #used for run_building_inundation.py
+
+#------------------------------------------------------------------------------
+# Interior region definitions
+#------------------------------------------------------------------------------
+
+#Land, to set the initial stage/water to be offcoast only
+poly_mainland = read_polygon(polygons_dir+'initial_condition.csv')
+
+#Land, to set the initial stage/water to be offcoast only
+poly_marina = read_polygon(polygons_dir+'initial_condition_marina.csv')
+
+# Initial bounding polygon for data clipping 
+poly_all = read_polygon(polygons_dir+'poly_all_extend.csv')
+res_poly_all = 100000*res_factor
+
+# Area of Interest 1 (Busselton)
+poly_aoi1 = read_polygon(polygons_dir+'busselton_1km.csv')
+res_aoi1 = 500*res_factor
+
+# Area of Interest 2 (Bunbury)
+poly_aoi2 = read_polygon(polygons_dir+'bunbury_1km.csv')
+res_aoi2 = 500*res_factor
+
+# Area of Significance 1 (Busselton)
+poly_aos1 = read_polygon(polygons_dir+'busselton_2km.csv')
+res_aos1 = 10000*res_factor
+
+# Area of Significance 2 (Bunbury)
+poly_aos2 = read_polygon(polygons_dir+'busselton_2km.csv')
+res_aos2 = 10000*res_factor
+
+# Refined areas
+# Polygon designed to islands
+poly_aos3 = read_polygon(polygons_dir+'island1.csv')
+res_aos3 = 10000*res_factor
+poly_aos4 = read_polygon(polygons_dir+'island2.csv')
+res_aos4 = 10000*res_factor
+
+# Shallow water 1
+poly_sw1 = read_polygon(polygons_dir+'coast_5km_d20m.csv')
+res_sw1 = 40000*res_factor
+
+# Combined all regions, must check that all are included!
+interior_regions = [[poly_aoi1,res_aoi1],[poly_aoi2,res_aoi2]
+                     ,[poly_aos1,res_aos1],[poly_aos2,res_aos2]
+                     ,[poly_aos3,res_aos3],[poly_aos4,res_aos4]
+                     ,[poly_sw1,res_sw1]]
+
+    
+trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+print 'min estimated number of triangles', trigs_min
+   
+#------------------------------------------------------------------------------
+# Clipping regions for export to asc and regions for clipping data
+# Final inundation maps should only be created in regions of the finest mesh
+#------------------------------------------------------------------------------
+
+# exporting asc grid for Busselton
+xminBusselton = 340000
+xmaxBusselton = 352000
+yminBusselton = 6271500
+ymaxBusselton = 6280000
+
+# exporting asc grid for Bunbury
+xminBunbury = 369000
+xmaxBunbury = 381000
+yminBunbury = 6308000
+ymaxBunbury = 6316500
+
Index: /anuga_work/production/busselton/standardised_version/run_busselton.py
===================================================================
--- /anuga_work/production/busselton/standardised_version/run_busselton.py	(revision 6253)
+++ /anuga_work/production/busselton/standardised_version/run_busselton.py	(revision 6253)
@@ -0,0 +1,243 @@
+"""Script for running 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_busselton.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_busselton.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
+from os import sep
+import os
+from os.path import dirname, basename
+from os import mkdir, access, F_OK
+from shutil import copy
+import time
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import File_boundary
+from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water import Field_boundary
+from Numeric import allclose
+from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
+from anuga_parallel.parallel_abstraction import get_processor_name
+from anuga.caching import myhash
+from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
+from anuga.fit_interpolate.benchmark_least_squares import mem_usage 
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
+from polygon import Polygon_function
+    
+# Application specific imports
+import project  # Definition of file names and polygons
+numprocs = 1
+myid = 0
+
+def run_model(**kwargs):
+    
+    #------------------------------------------------------------------------------
+    # Copy scripts to time stamped output directory and capture screen
+    # output to file
+    #------------------------------------------------------------------------------
+    print "Processor Name:",get_processor_name()
+
+    #copy script must be before screen_catcher
+
+    print 'output_dir',kwargs['output_dir']
+    
+    copy_code_files(kwargs['output_dir'],__file__, 
+             dirname(project.__file__)+sep+ project.__name__+'.py' )
+
+    store_parameters(**kwargs)
+
+    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
+
+    print "Processor Name:",get_processor_name()
+    
+    #-----------------------------------------------------------------------
+    # Domain definitions
+    #-----------------------------------------------------------------------
+
+    # Read in boundary from ordered sts file
+    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir_event,project.scenario_name))
+
+    # Reading the landward defined points, this incorporates the original clipping
+    # polygon minus the 100m contour
+    landward_bounding_polygon = read_polygon(project.landward_dir)
+
+    # Combine sts polyline with landward points
+    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
+    
+    # counting segments
+    N = len(urs_bounding_polygon)-1
+
+    # boundary tags refer to project.landward 4 points equals 5 segments start at N
+    boundary_tags={'back': [N+1,N+2,N+3,N+4, N+5], 'side': [N,N+6], 'ocean': range(N)}
+
+    #--------------------------------------------------------------------------
+    # Create the triangular mesh 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
+    #--------------------------------------------------------------------------
+
+    # IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
+    # causes problems with the ability to cache set quantity which takes alot of times
+        
+    print 'start create mesh from regions'
+
+    create_mesh_from_regions(bounding_polygon,
+                         boundary_tags=boundary_tags,
+                         maximum_triangle_area=project.res_poly_all,
+                         interior_regions=project.interior_regions,
+                         filename=project.meshes_dir_name,
+                         use_cache=True,
+                         verbose=True)
+   
+    #-------------------------------------------------------------------------
+    # Setup computational domain
+    #-------------------------------------------------------------------------
+    print 'Setup computational domain'
+
+    domain = Domain(project.meshes_dir_name, use_cache=False, verbose=True)
+    print 'memory usage before del domain',mem_usage()
+       
+    print domain.statistics()
+    print 'triangles',len(domain)
+    
+    kwargs['act_num_trigs']=len(domain)
+
+
+    #-------------------------------------------------------------------------
+    # Setup initial conditions
+    #-------------------------------------------------------------------------
+    print 'Setup initial conditions'
+
+    # sets the initial stage in the offcoast region only
+    IC = Polygon_function( [(project.poly_mainland, 0),(project.poly_marina, 0)], default = kwargs['tide'],
+                             geo_reference = domain.geo_reference)
+    domain.set_quantity('stage', IC)
+    #domain.set_quantity('stage',kwargs['tide'] )
+    domain.set_quantity('friction', kwargs['friction']) 
+    
+    print 'Start Set quantity',kwargs['elevation_file']
+
+    domain.set_quantity('elevation', 
+                        filename = kwargs['elevation_file'],
+                        use_cache = False,
+                        verbose = True,
+                        alpha = kwargs['alpha'])
+    print 'Finished Set quantity'
+
+##   #------------------------------------------------------
+##    # Distribute domain to implement parallelism !!!
+##    #------------------------------------------------------
+##
+##    if numprocs > 1:
+##        domain=distribute(domain)
+
+    #------------------------------------------------------
+    # Set domain parameters
+    #------------------------------------------------------
+    print 'domain id', id(domain)
+    domain.set_name(kwargs['scenario_name'])
+    domain.set_datadir(kwargs['output_dir'])
+    domain.set_default_order(2)                 # Apply second order scheme
+    domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
+    domain.set_store_vertices_uniquely(False)
+    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+    domain.tight_slope_limiters = 1
+    print 'domain id', id(domain)
+
+    #-------------------------------------------------------------------------
+    # Setup boundary conditions 
+    #-------------------------------------------------------------------------
+    print 'Available boundary tags', domain.get_boundary_tags()
+    print 'domain id', id(domain)
+   
+    boundary_urs_out=project.boundaries_dir_event + sep + project.scenario_name
+
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
+    
+    print 'Available boundary tags', domain.get_boundary_tags()
+    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
+                   domain, mean_stage= project.tide,
+                   time_thinning=1,
+                   default_boundary=Bd,
+                   use_cache=True,
+                   verbose = True,
+                   boundary_polygon=bounding_polygon)
+
+    domain.set_boundary({'back': Br,
+                         'side': Bd,
+                         'ocean': Bf}) 
+
+    kwargs['input_start_time']=domain.starttime
+
+    print'finish set boundary'
+
+    #----------------------------------------------------------------------------
+    # Evolve system through time
+    #--------------------------------------------------------------------
+    t0 = time.time()
+
+    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
+                       ,skip_initial_step = False): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean')
+
+    # these outputs should be checked with the resultant inundation map
+    x, y = domain.get_maximum_inundation_location()
+    q = domain.get_maximum_inundation_elevation()
+    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
+
+    print 'Simulation took %.2f seconds' %(time.time()-t0)
+
+    #kwargs 'completed' must be added to write the final parameters to file
+    kwargs['completed']=str(time.time()-t0)
+     
+    store_parameters(**kwargs)
+     
+    print 'memory usage before del domain1',mem_usage()
+    
+    
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    
+    kwargs={}
+    kwargs['file_name']=project.dir_comment
+    kwargs['finaltime']=project.finaltime
+    kwargs['output_dir']=project.output_run_time_dir
+    kwargs['elevation_file']=project.combined_dir_name+'.pts'
+    kwargs['scenario_name']=project.scenario_name
+    kwargs['tide']=project.tide
+    kwargs['alpha'] = project.alpha
+    kwargs['friction']=project.friction
+    #kwargs['num_cpu']=numprocs
+    #kwargs['host']=project.host
+    #kwargs['starttime']=project.starttime
+    #kwargs['yieldstep']=project.yieldstep
+    #kwargs['user']=project.user
+    #kwargs['time_thinning'] = project.time_thinning
+      
+    run_model(**kwargs)
+      
+   
