Index: /anuga_work/production/australia_ph2/dampier/build_dampier_ph2_250m.py
===================================================================
--- /anuga_work/production/australia_ph2/dampier/build_dampier_ph2_250m.py	(revision 6204)
+++ /anuga_work/production/australia_ph2/dampier/build_dampier_ph2_250m.py	(revision 6204)
@@ -0,0 +1,84 @@
+"""
+Script for building the elevation data to run a tsunami inundation scenario
+for dampier_ph2, WA, Australia.
+
+Input: elevation data from project_250m.py
+Output: pts file stored in project_250m.topographies_dir 
+The run_dampier.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_250m   # Definition of file names and polygons
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+
+copy_code_files(project_250m.output_build_time_dir,__file__, 
+               dirname(project_250m.__file__)+sep+ project_250m.__name__+'.py' )
+
+start_screen_catcher(project_250m.output_build_time_dir)
+
+print 'USER:    ', project_250m.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_250m.poly_all", project_250m.poly_all
+print "project_250m.combined_dir_name", project_250m.combined_dir_name
+
+# input elevation directory filenames
+onshore_in_dir_name = project_250m.onshore_in_dir_name
+
+# output elevation directory filenames
+onshore_dir_name = project_250m.onshore_dir_name
+
+# 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)
+
+# creates pts file for onshore DEM
+print "creates pts file for onshore DEM"
+dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
+
+# create onshore pts files
+print'create Geospatial data1 objects from topographies'
+G = Geospatial_data(file_name = onshore_dir_name + '.pts')
+
+#-------------------------------------------------------------------------------
+# Combine, clip and export dataset 
+#-------------------------------------------------------------------------------
+
+print'clip combined geospatial object by bounding polygon'
+G_clipped = G.clip(project_250m.poly_all)
+
+print'export combined DEM file'
+if access(project_250m.topographies_dir,F_OK) == 0:
+    mkdir (project_250m.topographies_dir)
+G_clipped.export_points_file(project_250m.combined_dir_name + '.pts')
+G_clipped.export_points_file(project_250m.combined_dir_name + '.txt') #Use for comparision in ARC
+
Index: /anuga_work/production/australia_ph2/dampier/export_results_max.py
===================================================================
--- /anuga_work/production/australia_ph2/dampier/export_results_max.py	(revision 6204)
+++ /anuga_work/production/australia_ph2/dampier/export_results_max.py	(revision 6204)
@@ -0,0 +1,113 @@
+"""
+Generates ascii grids of nominated areas -
+Input: sww file from run_perth.py
+       boundaries for grids from project_250m.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_250m, os
+import sys
+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_250m.output_dir
+
+##time_dir1 = '20090115_151520_run_final_0_250m_kvanputt'
+time_dir2 = '20090115_164756_run_final_0_250m_small_kvanputt'
+time_dirs = [time_dir2]   
+
+#cellsize = 5
+cellsize = 250
+
+timestep = None    # None means no timestep!
+#timestep = 0
+
+######
+# Set the special areas of interest.  If none, do: area='All'
+######
+
+#area = ['Sorrento', 'Geordie' ]  # 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 = ['stage']
+
+######
+# 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_250m.xmin%s' % which_area)
+                easting_max = eval('project_250m.xmax%s' % which_area)
+                northing_min = eval('project_250m.ymin%s' % which_area)
+                northing_max = eval('project_250m.ymax%s' % which_area)
+            except AttributeError:
+                print 'Unrecognized area name: %s' % which_area
+                break
+
+        for time_dir in time_dirs:
+
+            name1 = directory+time_dir+sep+project_250m.scenario_name
+##            name2 = directory+time_dir+sep+project_250m.scenario_name+'_time_39900_0'
+##            name3 = directory+time_dir+sep+project_250m.scenario_name+'_time_79800_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 = directory+time_dir+sep+project_250m_250m.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/dampier/get_timeseries.py
===================================================================
--- /anuga_work/production/australia_ph2/dampier/get_timeseries.py	(revision 6204)
+++ /anuga_work/production/australia_ph2/dampier/get_timeseries.py	(revision 6204)
@@ -0,0 +1,50 @@
+"""
+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 anuga.abstract_2d_finite_volumes.util import sww2csv_gauges,csv2timeseries_graphs
+import project
+
+directory = project.output_dir
+
+time_dir1 = '20081031_133353_run_final_0.6_68693_alpha0.1_kvanputt'+sep+'gauge'
+time_dir2 = '20081031_133511_run_final_0_68693_alpha0.1_kvanputt'+sep+'gauge'
+time_dir3 = '20081031_133624_run_final_0_27255_alpha0.1_kvanputt'+sep+'gauge'
+time_dir4 = '20081031_133735_run_final_0.6_27255_alpha0.1_kvanputt'+sep+'gauge'
+time_dir5 = '20081031_133841_run_final_0_27283_alpha0.1_kvanputt'+sep+'gauge'
+time_dir6 = '20081031_133925_run_final_0.6_27283_alpha0.1_kvanputt'+sep+'gauge'
+
+time_dirs = [time_dir1, time_dir2, time_dir3, time_dir4, time_dir5, time_dir6]
+
+for time_dir in time_dirs:
+    name = directory+time_dir+sep+project.scenario_name
+    gauge = project.gauges_dir_name
+    sww2csv_gauges(name+'.sww',gauge,
+                   quantities = ['stage','speed','depth','elevation'],
+                   verbose=True)
+
+    
+##    name1 = directory+time_dir+sep+project.scenario_name 
+####    name2 = directory+time_dir+sep+project.scenario_name+'_time_43380_0'
+####    name3 = directory+time_dir+sep+project.scenario_name+'_time_78360_0'
+##
+##    names = [name1]#, name2, name3]
+##
+##    for i, name in enumerate(names):
+##        gauge = project.gauges_dir_name
+####        out_name = 'gauge_%d_' % i
+##        print name
+##        sww2csv_gauges(name+'.sww',gauge,#out_name = out_name,
+##                        quantities = ['stage','speed','depth','elevation'],
+##                        verbose=True)
+##
Index: /anuga_work/production/australia_ph2/dampier/project_250m.py
===================================================================
--- /anuga_work/production/australia_ph2/dampier/project_250m.py	(revision 6204)
+++ /anuga_work/production/australia_ph2/dampier/project_250m.py	(revision 6204)
@@ -0,0 +1,168 @@
+"""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 = 'australia_ph2'
+scenario_name = 'dampier_ph2'
+scenario = 'dampier'
+
+# Model specific parameters. One or all can be changed each time the
+# run_scenario script is executed
+tide = 0                #0.6
+#event_number = 27283
+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
+
+#------------------------------------------------------------------------------
+# 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)+'_250m_small_'+str(user)
+
+#------------------------------------------------------------------------------
+# Input Data
+#------------------------------------------------------------------------------
+
+# elevation data used in build_perth.py
+# onshore data: format ascii grid with accompanying projection file
+onshore_name = 'grid_250m_2005' 
+
+# gauges - used in get_timeseries.py
+#gauge_name = 'perth.csv'
+
+# BOUNDING POLYGON - used in build_boundary.py and run_perth.py respectively
+# NOTE: when files are put together the points must be in sequence - for ease go clockwise!
+# Check the run_perth.py for boundary_tags
+# thinned ordering file from Hazard Map: format is index,latitude,longitude (with title)
+#order_filename = 'thinned_boundary_ordering.txt'
+#landward bounding points
+#landward = 'landward_bounding_polygon.txt'
+
+#------------------------------------------------------------------------------
+# Output Elevation Data
+#------------------------------------------------------------------------------
+# Output filename for elevation
+# this is a combination of all the data (utilisied in build_boundary)
+combined_name = scenario_name + '_combined_elevation_250m_small'
+
+#------------------------------------------------------------------------------
+# 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
+
+# where the output data sits
+onshore_dir_name = topographies_dir + onshore_name
+
+# where the combined elevation file sits
+combined_dir_name = topographies_dir + combined_name
+
+# where the mesh sits (this is created during the run_perth.py)
+meshes_dir_name = meshes_dir + scenario_name+ '_small.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_perth.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_perth.py
+output_run_time_dir = output_dir+run_time+dir_comment+sep       #used for run_perth.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
+##gauges_dir_name2 = gauges_dir + gauge_name2     #used for get_timeseries.py
+
+#------------------------------------------------------------------------------
+# Interior region definitions
+#------------------------------------------------------------------------------
+
+#Land, to set the initial stage/water to be offcoast only
+poly_mainland = read_polygon(polygons_dir+'initial_condition_small.csv')
+
+# Initial bounding polygon for data clipping 
+poly_all = read_polygon(polygons_dir+'poly_all_small.csv')
+res_poly_all = 125000*res_factor
+
+interior_regions = []
+     
+trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+print 'min estimated number of triangles', trigs_min
+   
+
Index: /anuga_work/production/australia_ph2/dampier/run_dampier_ph2_250m.py
===================================================================
--- /anuga_work/production/australia_ph2/dampier/run_dampier_ph2_250m.py	(revision 6204)
+++ /anuga_work/production/australia_ph2/dampier/run_dampier_ph2_250m.py	(revision 6204)
@@ -0,0 +1,250 @@
+"""Script for running a tsunami inundation scenario for Perth, WA, Australia.
+
+The scenario is defined by a triangular mesh created from project_250m.polygon,
+the elevation data is compiled into a pts file through build_perth.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_perth.py)
+       information from project file
+Outputs: sww file stored in project_250m.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
+from math import pi, sin, exp
+import time
+import sys
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water.shallow_water_domain import Transmissive_stage_zero_momentum_boundary
+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 anuga.shallow_water import Time_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_250m  # 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_250m.__file__)+sep+ project_250m.__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_250m.boundaries_dir_event,project_250m.scenario_name))
+##
+##    # Reading the landward defined points, this incorporates the original clipping
+##    # polygon minus the 100m contour
+##    landward_bounding_polygon = read_polygon(project_250m.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_250m.landward 4 points equals 5 segments start at N
+##    boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}
+
+##    land = range (0,11)
+##    sea = range (12,122)
+##    bounding_polygon = project_250m.poly_all
+##    boundary_tags={'back': land, 'side': [11,122], 'ocean': sea}
+   
+    sea = range (4,22)
+    bounding_polygon = project_250m.poly_all
+    boundary_tags={'back': [0,1,2,23], 'side': [3,22], 'ocean': sea}
+
+    #--------------------------------------------------------------------------
+    # Create the triangular mesh based on overall clipping polygon with a tagged
+    # boundary and interior regions defined in project_250m.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_250m.res_poly_all,
+                         interior_regions=project_250m.interior_regions,
+                         filename=project_250m.meshes_dir_name,
+                         use_cache=False,
+                         verbose=True)
+   
+    #-------------------------------------------------------------------------
+    # Setup computational domain
+    #-------------------------------------------------------------------------
+    print 'Setup computational domain'
+
+    domain = Domain(project_250m.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_250m.poly_mainland, 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_250m.boundaries_dir_event + sep + project_250m.scenario_name
+##    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_250m.tide,
+##                   time_thinning=1,
+##                   default_boundary=Bd,
+##                   use_cache=True,
+##                   verbose = True,
+##                   boundary_polygon=bounding_polygon)
+
+    Br = Reflective_boundary(domain)
+    Bt = Transmissive_stage_zero_momentum_boundary(domain)
+    Rb = Time_boundary (domain, f=lambda t: [(sin(t*pi/600) * exp(-(t/3600)**2)), 0.0, 0.0])   
+
+    domain.set_boundary({'back': Br,
+                         'side': Bt,
+                         'ocean': Rb}) 
+
+    kwargs['input_start_time']=domain.starttime
+
+    print'finish set boundary'
+
+    #----------------------------------------------------------------------------
+    # Evolve system through time
+    #--------------------------------------------------------------------
+    t0 = time.time()
+
+    for t in domain.evolve(yieldstep = project_250m.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_250m.dir_comment
+    kwargs['finaltime']=project_250m.finaltime
+    kwargs['output_dir']=project_250m.output_run_time_dir
+    kwargs['elevation_file']=project_250m.combined_dir_name+'.pts'
+    kwargs['scenario_name']=project_250m.scenario_name
+    kwargs['tide']=project_250m.tide
+    kwargs['alpha'] = project_250m.alpha
+    kwargs['friction']=project_250m.friction
+     
+    run_model(**kwargs)
+      
+   
