Index: /anuga_work/production/exmouth_2006/build_exmouth.py
===================================================================
--- /anuga_work/production/exmouth_2006/build_exmouth.py	(revision 4465)
+++ /anuga_work/production/exmouth_2006/build_exmouth.py	(revision 4465)
@@ -0,0 +1,137 @@
+"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project.py
+The output sww file is stored in project.output_time_dir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a simulated submarine landslide.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
+"""
+
+#------------------------------------------------------------------------------
+# 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.abstract_2d_finite_volumes.util 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.combined_dir_name",project.combined_dir_name
+
+# topography directory filenames
+onshore_in_dir_name = project.onshore_in_dir_name
+onshore1_in_dir_name = project.onshore1_in_dir_name
+coast_in_dir_name = project.coast_in_dir_name
+#island_in_dir_name = project.island_in_dir_name
+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
+
+onshore_dir_name = project.onshore_dir_name
+onshore1_dir_name = project.onshore1_dir_name
+coast_dir_name = project.coast_dir_name
+#island_dir_name = project.island_dir_name
+offshore_dir_name = project.offshore_dir_name
+#offshore_dir_name1 = project.offshore_dir_name1
+#offshore_dir_name2 = project.offshore_dir_name2
+
+# 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(onshore1_in_dir_name, basename_out=onshore1_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)
+dem2pts(onshore1_dir_name, use_cache=True, verbose=True)
+
+print'create Geospatial data1 objects from topographies',onshore_dir_name + '.pts'
+G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
+print'create Geospatial data2 objects from topographies',onshore1_dir_name + '.pts'
+G2 = Geospatial_data(file_name = onshore1_dir_name + '.pts')
+print'create Geospatial data3 objects from coast', coast_in_dir_name + '.txt'
+G3 = Geospatial_data(file_name = coast_in_dir_name + '.txt')
+print'create Geospatial data3 objects from offshore',offshore_in_dir_name + '.txt'
+G_off = Geospatial_data(file_name = offshore_in_dir_name + '.txt')
+
+print'add all geospatial objects'
+G = G1 + G2 + G3 + G_off 
+
+print'export combined DEM file'
+if access(project.topographies_dir,F_OK) == 0:
+    mkdir (project.topographies_dir)
+print'export',project.combined_dir_name+'_unclipped' + '.txt'
+G.export_points_file(project.combined_dir_name+ '.txt')
+
+print'split'
+G_small, G_other = G_clipped.split(.10,verbose=True)
+
+print 'export', project.combined_small_dir_name + '.txt'
+G_small.export_points_file(project.combined_small_dir_name + '.txt')
+#G_clipped.export_points_file(project.combined_dir_name + '.xya')
+
+
+'''
+print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
+G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
+print'split'
+G_all_1, G_all_2 = G_all.split(.10)
+print'export 1'
+G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya')
+print'export 2'
+G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
+'''
+
+#-------------------------------------------------------------------------
+# Convert URS to SWW file for boundary conditions 
+#-------------------------------------------------------------------------
+print 'starting to create boundary conditions'
+boundaries_in_dir_name = project.boundaries_in_dir_name
+
+from anuga.shallow_water.data_manager import urs2sww, urs_ungridded2sww
+
+print 'boundaries_in_dir_name',boundaries_in_dir_name
+print 'project.boundaries_dir_name',project.boundaries_dir_name
+
+urs_ungridded2sww(boundaries_in_dir_name, project.boundaries_dir_name,
+                  verbose=True, mint=4000, maxt=35000, zscale=1)
+
+
+
+
+
+
Index: /anuga_work/production/exmouth_2006/project.py
===================================================================
--- /anuga_work/production/exmouth_2006/project.py	(revision 4464)
+++ /anuga_work/production/exmouth_2006/project.py	(revision 4465)
@@ -1,142 +1,170 @@
-# -*- coding: cp1252 -*-
 """Common filenames and locations for topographic data, meshes and outputs.
+Also includes origin for slump scenario.
 """
 
-from os import sep, environ, getenv, getcwd
-from os.path import expanduser
+from os import sep, environ, getenv, getcwd,umask
+from os.path import expanduser, basename
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon, number_mesh_triangles
 import sys
+from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
 from time import localtime, strftime, gmtime
-from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
-#from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
+from anuga.utilities.system_tools import get_user_name, get_host_name
 
-if sys.platform == 'win32':
-    home = getenv('INUNDATIONHOME')
-    user = getenv('USERPROFILE')
+codename = 'project.py'
 
-else:    
-    home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation')     
-    user = getenv('LOGNAME')
-    print 'USER:', user
-
-# INUNDATIONHOME is the inundation directory, not the data directory.
-home += sep +'data'
+home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent dir   
+user = get_user_name()
+host = get_host_name()
+#needed when running using mpirun, mpirun doesn't inherit umask from .bashrc
+umask(002)
 
 #Making assumptions about the location of scenario data
 state = 'western_australia'
-scenario_dir_name = 'exmouth_tsunami_scenario_2006'
+scenario_name = 'exmouth'
+scenario = 'exmouth_tsunami_scenario'
 
-# onshore data provided by WA DLI
-onshore_name =  'first' # original
+#time stuff
+time = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
+build_time = time+'_build'
+run_time = time+'_run'
 
-# AHO + DPI data
-offshore_name1 = 'XY100011610'
-offshore_name2 = 'XY100011611'
-offshore_name3 = 'XY100011613'
-offshore_name4 = 'XY100011614'
-offshore_name5 = 'XY100011616'
-offshore_name6 = 'XY100011617'
-offshore_name7 = 'XY100011618'
-offshore_name8 = 'XY100011621'
-offshore_name9 = 'XY100011623'
-offshore_name10 = 'XY100011745'
-offshore_name11 = 'XY100011746'
-offshore_name12 = 'XY100017530'
-offshore_name13 = 'XY100017532'
-offshore_name14 = 'XY100017538'
-offshore_name15 = 'XY100017540'
-offshore_name16 = 'XYBR66'
-offshore_name17 = 'XYBR70'
-offshore_name18 = 'XYBR80'
-offshore_name19 = 'XYBR88'
-offshore_name20 = 'XYBR93'
-offshore_name21 = 'XYBR0110'
-offshore_name22 = 'XYWADPI'
+#tide = -1.4
+#tide = 0
+tide = 1.4
 
-# developed by NM&I
-coast_name = 'coastline'
+#Maybe will try to make project a class to allow these parameters to be passed in.
+alpha = 0.1
+friction=0.01
+starttime=3600
+finaltime=25000
+setup='trial'
 
-boundary_basename = 'SU-AU' # Mw ?
+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=1
+    yieldstep=60
 
-#swollen/ all data output
-basename = 'source'
-codename = 'project.py'
+dir_comment='_'+setup+'_'+str(tide)+'_'+str(user)
 
-#Derive subdirectories and filenames
-local_time = strftime('%Y%m%d_%H%M%S',gmtime()) 
-meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
-datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
-gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
-polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
-boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
-outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
-outputtimedir = outputdir + local_time + sep
-polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
+onshore_name = 'DLI' 
+onshore_name1 = 'DTED' 
+# offshore
+offshore_name = 'Exmouth_bathymetry'
+#offshore_name1 = 'inferred_north'
+#offshore_name2 = 'inferred_south'
+coast_name = 'Exmouth_coastline'
 
-gauge_filename = gaugedir + 'exmouth_gauges.csv'
-community_filename = gaugedir + 'CHINS_v2.csv'
-community_exmouth = gaugedir + 'community_exmouth.csv'
-codedir = getcwd()+sep                            
-codedirname = codedir + 'project.py'
-meshname = outputtimedir + 'mesh_' + basename
+#final topo name
+combined_name ='exmouth_combined_elevation'
+#combined_name1 ='exmouth_combined_elevation1'
+#combined_name_unclipped1 ='exmouth_combined_elevation_unclipped1'
+combined_small_name = 'exmouth_combined_elevation_small'
 
-# Necessary if using point datasets, rather than grid
-onshore_dem_name = datadir + onshore_name
-offshore_dem_name1 = datadir + offshore_name1
-offshore_dem_name2 = datadir + offshore_name2
-offshore_dem_name3 = datadir + offshore_name3
-offshore_dem_name4 = datadir + offshore_name4
-offshore_dem_name5 = datadir + offshore_name5
-offshore_dem_name6 = datadir + offshore_name6
-offshore_dem_name7 = datadir + offshore_name7
-offshore_dem_name8 = datadir + offshore_name8
-offshore_dem_name9 = datadir + offshore_name9
-offshore_dem_name10 = datadir + offshore_name10
-offshore_dem_name11 = datadir + offshore_name11
-offshore_dem_name12 = datadir + offshore_name12
-offshore_dem_name13 = datadir + offshore_name13
-offshore_dem_name14 = datadir + offshore_name14
-offshore_dem_name15 = datadir + offshore_name15
-offshore_dem_name16 = datadir + offshore_name16
-offshore_dem_name17 = datadir + offshore_name17
-offshore_dem_name18 = datadir + offshore_name18
-offshore_dem_name19 = datadir + offshore_name19
-offshore_dem_name20 = datadir + offshore_name20
-offshore_dem_name21 = datadir + offshore_name21
-offshore_dem_name22 = datadir + offshore_name22
+anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
 
-coast_dem_name = datadir + coast_name
+topographies_in_dir = home+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
+topographies_dir = anuga_dir+'topographies'+sep
 
-combined_dem_name   = datadir + 'exmouth_combined_elevation'
+# input topo file location
+onshore_in_dir_name = topographies_in_dir + onshore_name
+onshore1_in_dir_name = topographies_in_dir + onshore_name1
+coast_in_dir_name = topographies_in_dir + coast_name
+offshore_in_dir_name = topographies_in_dir + offshore_name
+#offshore_in_dir_name1 = topographies_in_dir + offshore_name1
+#offshore_in_dir_name2 = topographies_in_dir + offshore_name2
 
-#buildings_filename = gauges_dir + 'exmouth_res_nexis.csv'
+onshore_dir_name = topographies_dir + onshore_name
+onshore1_dir_name = topographies_dir + onshore_name1
+coast_dir_name = topographies_dir + coast_name
+offshore_dir_name = topographies_dir + offshore_name
+#offshore_dir_name1 = topographies_dir + offshore_name1
+#offshore_dir_name2 = topographies_dir + offshore_name2
+
+#final topo files
+combined_dir_name = topographies_dir + combined_name
+#combined_dir_name_unclipped1 = topographies_dir + combined_name_unclipped1
+#combined_dir_name1 = topographies_dir + combined_name1
+combined_small_name_dir = topographies_dir + combined_small_name
+
+meshes_dir = anuga_dir+'meshes'+sep
+meshes_dir_name = meshes_dir + scenario_name
+
+polygons_dir = anuga_dir+'polygons'+sep
+tide_dir = anuga_dir+'tide_data'+sep
+
+#boundaries locations
+boundaries_name = 'exmouth_3854_17042007'
+
+boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'1_10000'+sep
+boundaries_in_dir_name = boundaries_in_dir + boundaries_name
+boundaries_dir = anuga_dir+'boundaries'+sep
+boundaries_dir_name = boundaries_dir + boundaries_name
+
+#output locations
+output_dir = anuga_dir+'outputs'+sep
+output_build_time_dir = output_dir+build_time+sep
+output_run_time_dir = output_dir +run_time+dir_comment+sep
+output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
+
+#gauges
+gauge_name = 'exmouth_gauges.csv'
+gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep
+gauges_dir_name = gauges_dir + gauge_name
+
+community_filename = gauges_dir + 'CHINS_v2.csv'
+community_exmouth = gauges_dir + 'community_exmouth.csv'
+
 buildings_filename_damage_out = 'exmouth_res_nexis_modified.csv'
-
-community_filename = gaugedir + 'CHINS_v2.csv'
-community_exmouth = gaugedir + 'community_exmouth.csv'
 
 ###############################
 # Domain definitions
 ###############################
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
 
-# bounding box for clipping MOST/URS output (much bigger than study area)
-##south = degminsec2decimal_degrees(-19,0,0)
-##north = degminsec2decimal_degrees(-17,15,0)
-##west  = degminsec2decimal_degrees(120,0,0)
-##east  = degminsec2decimal_degrees(122,30,0) 
-##
-##d0 = [south, west]
-##d1 = [south, east]
-##d2 = [north, east]
-##d3 = [north, west]
-##poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
-##refzone = zone
+poly_all = read_polygon(polygons_dir+'extent_zone50a.csv')
 
-# bounding polygon for study area
-poly_all = read_polygon(polygondir+'extent_zone50a.csv')
-
-# plot bounding polygon and make sure BC info surrounds it
-#plot_polygons([poly_all, poly_bc],'boundingpoly',verbose=False)
 print 'Area of bounding polygon', polygon_area(poly_all)/1000000.0
 
+res_poly_all = 150000*res_factor
+
+###############################
+# Interior region definitions
+###############################
+'''
+poly_0 = read_polygon(polygons_dir+'neg20_coast_contour_pts.csv')
+res_0 = 20000*res_factor
+
+poly_1 = read_polygon(polygons_dir+'exmouth_north_coast_inside_extent.csv')
+res_1 = 5000*res_factor
+
+poly_2 = read_polygon(polygons_dir+'exmouth_south_coast_inside_extent.csv')
+res_2 = 5000*res_factor
+
+poly_3 = read_polygon(polygons_dir+'exmouth_town_pts.csv')
+res_3 = 2000*res_factor
+
+poly_4 = read_polygon(polygons_dir+'exmouth_inner_town_pts.csv')
+res_4 = 500*res_factor
+
+interior_regions = [[poly_0,res_0],[poly_1,res_1],[poly_2,res_2]
+                     ,[poly_3,res_3],[poly_4,res_4]]
+
+trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+
+print 'min number triangles', trigs_min
+
+poly_mainland = read_polygon(polygons_dir+'Initial_Condition.csv')
+'''
 ###################################################################
 # Clipping regions for export to asc and regions for clipping data
@@ -148,37 +176,2 @@
 northingmin = 7983427.73
 northingmax = 8032834.52
-
-###############################
-# Interior region definitions
-###############################
-'''
-# exmouth digitized polygons
-poly_exmouth1 = read_polygon(polygondir+'exmouth_Local_Polygon_update.csv')
-poly_exmouth2 = read_polygon(polygondir+'exmouth_Close2_update.csv')
-poly_exmouth3 = read_polygon(polygondir+'exmouth_Coast_update.csv')
-#poly_exmouth4 = read_polygon(polygondir+'Cable_Beach_revised.csv')
-
-plot_polygons([poly_all,poly_exmouth1,poly_exmouth2,poly_exmouth3],'boundingpoly2',verbose=False)
-print 'Area of local polygon', polygon_area(poly_exmouth1)/1000000.0
-print 'Area of close polygon', polygon_area(poly_exmouth2)/1000000.0
-print 'Area of coastal polygon', polygon_area(poly_exmouth3)/1000000.0
-#print 'Area of cable beach polygon', polygon_area(poly_exmouth4)/1000000.0
-
-for i in poly_exmouth3:
-    v = is_inside_polygon(i,poly_exmouth1, verbose=False)
-    if v == False: print v
-
-def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
-    from anuga.utilities.polygon import polygon_area
-    
-    # TO DO check if any of the regions fall inside one another
-    no_triangles = 0.0
-    area = polygon_area(bounding_poly)
-    for i,j in interior_regions:
-        this_area = polygon_area(i)
-        no_triangles += this_area/j
-        area -= this_area
-        print j, this_area/1000000., area/1000000.
-    no_triangles += area/remainder_res
-    return int(no_triangles/0.7)
-'''
Index: /anuga_work/production/exmouth_2006/run_exmouth.py
===================================================================
--- /anuga_work/production/exmouth_2006/run_exmouth.py	(revision 4465)
+++ /anuga_work/production/exmouth_2006/run_exmouth.py	(revision 4465)
@@ -0,0 +1,245 @@
+"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project_urs.py
+The output sww file is stored in project_urs.output_run_time_dir 
+
+The scenario is defined by a triangular mesh created from project_urs.polygon,
+the elevation data and a simulated tsunami generated with URS code.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
+"""
+
+#------------------------------------------------------------------------------
+# 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 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.pmesh.mesh_interface import create_mesh_from_regions
+from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files,store_parameters
+from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
+from anuga_parallel.parallel_abstraction import get_processor_name
+from anuga.caching import myhash
+# Application specific imports
+import project_urs                 # Definition of file names and polygons
+
+def run_model(**kwargs):
+    
+    tide = kwargs['tide']
+    alpha = kwargs['alpha']
+    friction = kwargs['friction']
+    time_thinning = kwargs['time_thinning']
+    scenario_name = kwargs['aa_scenario_name']
+    
+    
+    
+    #------------------------------------------------------------------------------
+    # Copy scripts to time stamped output directory and capture screen
+    # output to file
+    #------------------------------------------------------------------------------
+
+    #copy script must be before screen_catcher
+    print 'tide',tide
+    kwargs['est_num_trigs']=project_urs.trigs_min
+    kwargs['num_cpu']=numprocs
+    kwargs['host']=project_urs.host
+    kwargs['res_factor']=project_urs.res_factor
+    kwargs['starttime']=project_urs.starttime
+    kwargs['yieldstep']=project_urs.yieldstep
+    kwargs['finaltime']=project_urs.finaltime
+    
+    kwargs['output_dir']=project_urs.output_run_time_dir
+    kwargs['bathy_file']=project_urs.combined_dir_name_unclipped1+'.txt'
+#    kwargs['bathy_file']=project_urs.combined_small_dir_name + '.pts'
+    kwargs['boundary_file']=project_urs.boundaries_dir_name + '.sww'
+
+    print 'output_dir',kwargs['output_dir']
+    if myid == 0:
+        copy_code_files(kwargs['output_dir'],__file__, 
+                 dirname(project_urs.__file__)+sep+ project_urs.__name__+'.py' )
+
+        store_parameters(**kwargs)
+
+    barrier()
+
+    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
+
+    print "Processor Name:",get_processor_name()
+
+    # filenames
+#    meshes_dir_name = project_urs.meshes_dir_name+'.msh'
+
+    # creates copy of code in output dir
+    print 'min triangles', project_urs.trigs_min,
+    print 'Note: This is generally about 20% less than the final amount'
+
+    #--------------------------------------------------------------------------
+    # Create the triangular mesh based on overall clipping polygon with a
+    # tagged
+    # boundary and interior regions defined in project_urs.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
+    if myid == 0:
+    
+        print 'start create mesh from regions'
+
+        create_mesh_from_regions(project_urs.poly_all,
+                         boundary_tags={'back': [1, 2], 'side': [0,3],
+                                        'ocean': [4]},
+                             maximum_triangle_area=project_urs.res_poly_all,
+                             interior_regions=project_urs.interior_regions,
+                             filename=project_urs.meshes_dir_name+'.msh',
+                             use_cache=False,
+                             verbose=True)
+    barrier()
+    
+
+    #-------------------------------------------------------------------------
+    # Setup computational domain
+    #-------------------------------------------------------------------------
+    print 'Setup computational domain'
+
+    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
+    #above don't work
+    domain = Domain(project_urs.meshes_dir_name+'.msh', use_cache=False, verbose=True)
+       
+    print domain.statistics()
+    print 'triangles',len(domain)
+    
+    kwargs['act_num_trigs']=len(domain)
+
+
+    #-------------------------------------------------------------------------
+    # Setup initial conditions
+    #-------------------------------------------------------------------------
+    if myid == 0:
+
+        print 'Setup initial conditions'
+
+        from polygon import Polygon_function
+        #following sets the stage/water to be offcoast only
+        IC = Polygon_function( [(project_urs.poly_mainland, -1.0)], default = tide,
+                                 geo_reference = domain.geo_reference)
+        domain.set_quantity('stage', IC)
+        domain.set_quantity('friction', friction) 
+        
+        print 'Start Set quantity'
+
+        domain.set_quantity('elevation', 
+                            filename = kwargs['bathy_file'],
+                            use_cache = True,
+                            verbose = True,
+                            alpha = alpha)
+        print 'Finished Set quantity'
+    barrier()
+
+
+    #------------------------------------------------------
+    # Distribute domain to implement parallelism !!!
+    #------------------------------------------------------
+
+    if numprocs > 1:
+        domain=distribute(domain)
+
+    #------------------------------------------------------
+    # Set domain parameters
+    #------------------------------------------------------
+    print 'domain id', id(domain)
+    domain.set_name(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.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
+    print 'domain id', id(domain)
+    domain.beta_h = 0
+    #domain.limit2007 = 1
+
+    #-------------------------------------------------------------------------
+    # Setup boundary conditions 
+    #-------------------------------------------------------------------------
+    print 'Available boundary tags', domain.get_boundary_tags()
+    print 'domain id', id(domain)
+    #print 'Reading Boundary file',project_urs.boundaries_dir_namea + '.sww'
+
+    Bf = Field_boundary(kwargs['boundary_file'],
+                    domain, time_thinning=time_thinning, mean_stage=tide, 
+                    use_cache=True, verbose=True)
+                    
+    kwargs['input_start_time']=domain.starttime
+
+    print 'finished reading boundary file'
+
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([tide,0,0])
+
+    print'set_boundary'
+
+    domain.set_boundary({'back': Br,
+                         'side': Bd,
+                         'ocean': Bf}) 
+    print'finish set boundary'
+
+    #----------------------------------------------------------------------------
+    # Evolve system through time
+    #----------------------------------------------------------------------------
+
+    t0 = time.time()
+
+    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean')     
+
+    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = 21600
+                       ,skip_initial_step = True): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean')     
+    
+    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
+                       ,skip_initial_step = True): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean')   
+
+    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 'That took %.2f seconds' %(time.time()-t0)
+
+    #kwargs 'completed' must be added to write the final parameters to file
+    kwargs['completed']=str(time.time()-t0)
+    if myid == 0:
+        store_parameters(**kwargs)
+    barrier()
+
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+
+    run_model(file_name=project_urs.home+'detail.csv', aa_scenario_name=project_urs.scenario_name,
+              ab_time=project_urs.time, res_factor= project_urs.res_factor, tide=project_urs.tide, user=project_urs.user,
+              alpha = project_urs.alpha, friction=project_urs.friction,
+              time_thinning = project_urs.time_thinning,
+              dir_comment=project_urs.dir_comment)
+
+
