Index: /anuga_work/development/boxingday08/combine_good_data.py
===================================================================
--- /anuga_work/development/boxingday08/combine_good_data.py	(revision 5247)
+++ /anuga_work/development/boxingday08/combine_good_data.py	(revision 5247)
@@ -0,0 +1,98 @@
+from anuga.geospatial_data.geospatial_data import *
+import project
+from os.path import exists
+from os import sep
+
+def get_zone(filename):
+    from anuga.coordinate_transforms.redfearn import redfearn
+    fr=open(filename+'.xyz','r')
+    line =fr.readline()
+    parts = line.rstrip().split()
+    zone, easting, northing = redfearn(float(parts[1]),float(parts[0]))
+    fr.close
+    print easting,northing
+    return zone
+
+def clip_csv(filename,poly):
+    fr=open(filename+'.csv','r')
+    fw=open(filename+'_clipped'+'.csv','w')
+    fw.write('x,y,elevation\n')
+    #extremes = [East,West,North,South] #May only work for positive values
+    extremes=[-1.0e10,1.0e10,-1.0e10,1.0e10]
+    for points in poly:
+        if points[0]>extremes[0]:
+            extremes[0]=points[0]
+        elif points[0]<extremes[1]:
+            extremes[1]=points[0]
+        if points[1]>extremes[2]:
+            extremes[2]=points[1]
+        elif points[1]<extremes[3]:
+            extremes[3]=points[1]
+    print "polygon_extent= ", extremes  
+
+    start = True
+    for line in fr:
+        if start:
+            #Read in Header
+            start=False
+        else:
+            parts = line.rstrip().split(',')
+            if float(parts[0])<extremes[0] and float(parts[0])>extremes[1] and float(parts[1])<extremes[2] and float(parts[1])>extremes[3]:
+                fw.write(','.join([parts[0],parts[1],parts[2]])+'\n')
+    fr.close
+    fw.close
+    
+def xyz2csv(filename):
+    """Convert .xyz file with an arbitrary number of space delimeters to .csv
+    format.
+
+    File format
+    xyz:
+        (No Header)
+        longitude, latitude, elevation
+    
+    csv:
+        (Header x,y,elevation)
+        easting, northing, elevation
+    """
+    from anuga.coordinate_transforms.redfearn import redfearn
+    
+    #Write header
+    fr=open(filename+'.xyz','r')
+    fw=open(filename+'.csv','w')
+    fw.write('x,y,elevation\n')
+    for line in fr:
+        parts = line.rstrip().split()
+        zone, easting, northing = redfearn(float(parts[1]),float(parts[0]))
+        fw.write(','.join([str(easting),str(northing),parts[2]])+'\n')
+        fr.close
+        fw.close
+
+filenames=['patong_bay_1s_grid','andaman_3s_grid']
+
+dir='data'
+print 'Converting .xya files to .csv format'
+for filename in filenames:
+    if not exists(dir+sep+filename+'.csv'):
+        xyz2csv(dir+sep+filename)
+    else:
+        print dir+sep+filename+'.csv already exists ensure the .xya changed has not changed'
+
+print 'Clipping .csv files'
+filename='andaman_3s_grid'
+if not exists(dir+sep+filename+'_clipped'+'.csv'):
+    clip_csv(dir+sep+filename,project.bounding_polygon)
+
+print 'Creating geospatial data'
+if exists(project.good_combined_dir_name +'.pts'):
+    print project.good_combined_dir_name +'.pts'+' already exists ensure the input .csv files have not changed'
+else:
+    G0 = Geospatial_data(file_name=dir+sep+filenames[0]+'.csv')
+    G1=Geospatial_data(file_name=dir+sep+filenames[1]+'_clipped.csv')
+    G1.clip(project.bounding_polygon)
+    G1.clip_outside(G0)
+
+    G = G0 + G1
+
+    print 'Creating .pts file'
+    G.export_points_file(project.good_combined_dir_name +'.pts')
Index: /anuga_work/development/boxingday08/export_results.py
===================================================================
--- /anuga_work/development/boxingday08/export_results.py	(revision 5247)
+++ /anuga_work/development/boxingday08/export_results.py	(revision 5247)
@@ -0,0 +1,68 @@
+
+import sys
+
+from anuga.shallow_water.data_manager import sww2dem
+from os import sep
+
+
+#time_dir = '20070322_035135' # birubi
+
+#directory = project_slide.outputdir
+#name = directory + time_dir + sep + project_slide.basename
+#name = 'test500000.035000'
+#name = 'test'
+name = 'good_simulationsourcefriction20080404_140016'
+is_parallel = False
+if is_parallel == True: nodes = 4
+print 'output dir:', name
+
+which_var = 2
+if which_var == 0:  # Stage
+    outname = name + '_stage-eq-thurs-test'
+    quantityname = 'stage'
+
+if which_var == 1:  # Absolute Momentum
+    outname = name + '_momentum_i1'
+    quantityname = '(xmomentum**2 + ymomentum**2)**0.5'  
+
+if which_var == 2:  # Depth
+    outname = name + '_depth'
+    quantityname = 'stage-elevation'  
+
+if which_var == 3:  # Speed
+    outname = name + '_speed_i0'
+    quantityname = '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'  #Speed
+
+if which_var == 4:  # Elevation
+    outname = name + '_elevation'
+    quantityname = 'elevation'  #Elevation
+
+if is_parallel == True:
+    for i in range(0,nodes):
+        namei = name + '_P%d_%d' %(i,nodes)
+        outnamei = outname + '_P%d_%d' %(i,nodes)
+        print 'start sww2dem for sww file %d' %(i)
+        sww2dem(namei, basename_out = outnamei,
+                    quantity = quantityname,
+                    timestep = 1,
+                    cellsize = 100,      
+                    #easting_min = project.e_min_area,
+                    #easting_max = project.e_max_area,
+                    #northing_min = project.n_min_area,
+                    #northing_max = project.n_max_area,        
+                    reduction = max, 
+                    verbose = True,
+                    format = 'asc')
+else:
+    print 'start sww2dem'
+    sww2dem(name, basename_out = outname,
+                quantity = quantityname,
+                #timestep = 0,
+                cellsize = 50,      
+                #easting_min = 416920,
+                #easting_max = 425310,
+                #northing_min = 877864,
+                #northing_max = 869570,        
+                reduction = max, 
+                verbose = True,
+                format = 'asc')
Index: /anuga_work/development/boxingday08/mesh_polygons/20m.csv
===================================================================
--- /anuga_work/development/boxingday08/mesh_polygons/20m.csv	(revision 5247)
+++ /anuga_work/development/boxingday08/mesh_polygons/20m.csv	(revision 5247)
@@ -0,0 +1,19 @@
+419517.9812,889806.5043
+422812.9129,886474.9623
+422520.03,882374.603
+422666.4714,879262.7231
+425668.5202,876626.7779
+426181.0652,873222.0152
+424423.7683,869194.8766
+424057.6648,865607.0622
+423691.5613,862824.6755
+424643.4304,859310.0817
+427425.8171,854184.6326
+429036.6726,846359.013
+425961.4031,841892.5501
+417467.8016,850450.3768
+409999.2899,866851.8141
+407070.4618,879299.3335
+406411.4755,891380.7494
+405972.1513,898776.0403
+404214.8544,903791.6584
Index: /anuga_work/development/boxingday08/mesh_polygons/bay.csv
===================================================================
--- /anuga_work/development/boxingday08/mesh_polygons/bay.csv	(revision 5247)
+++ /anuga_work/development/boxingday08/mesh_polygons/bay.csv	(revision 5247)
@@ -0,0 +1,20 @@
+419794.8167,877601.9209
+421522.398,876986.0169
+422395.2339,876651.7685
+424117.1566,876203.4707
+424570.5829,875076.5158
+425270.3827,874020.5785
+424899.7029,872290.7394
+424117.1566,870643.2736
+422469.6908,869366.4876
+421316.4647,868336.8215
+419957.3055,867389.5286
+418762.8928,866030.3693
+417321.3602,868913.4345
+416374.0673,871013.9534
+416085.7608,873196.8456
+416168.1341,876038.7241
+416456.4406,879127.7225
+417444.9201,881022.3082
+419339.5058,880075.0153
+420698.6651,879210.0958
Index: /anuga_work/development/boxingday08/mesh_polygons/boundary_new_pts.csv
===================================================================
--- /anuga_work/development/boxingday08/mesh_polygons/boundary_new_pts.csv	(revision 5247)
+++ /anuga_work/development/boxingday08/mesh_polygons/boundary_new_pts.csv	(revision 5247)
@@ -0,0 +1,17 @@
+416497.16571,940728.21738
+420060.14317,930349.10914
+433227.66854,903549.32220
+443916.60090,921054.38535
+461731.48818,930813.84533
+470406.56372,922293.68185
+477997.25482,890846.53302
+491629.51638,891001.44508
+509289.49159,871017.78892
+506501.07445,856610.96704
+493798.28527,839570.64008
+483574.08909,826403.11471
+440973.27170,828262.05947
+420060.14317,842359.05722
+400231.39907,896888.10348
+395119.30098,931433.49359
+406892.61779,939953.65706
Index: /anuga_work/development/boxingday08/mesh_polygons/extent.csv
===================================================================
--- /anuga_work/development/boxingday08/mesh_polygons/extent.csv	(revision 5247)
+++ /anuga_work/development/boxingday08/mesh_polygons/extent.csv	(revision 5247)
@@ -0,0 +1,9 @@
+419336.424,810100.845
+342405.0281,711455.8026
+274649.9152,723352.9603
+272089.092,972393.0131
+347633.3754,968551.7784
+427979.2022,885965.2313
+427659.0993,875721.9386
+429259.6138,861317.3083
+436301.8775,840830.723
Index: /anuga_work/development/boxingday08/mesh_polygons/north_island.csv
===================================================================
--- /anuga_work/development/boxingday08/mesh_polygons/north_island.csv	(revision 5247)
+++ /anuga_work/development/boxingday08/mesh_polygons/north_island.csv	(revision 5247)
@@ -0,0 +1,29 @@
+352521.0386,959228.8827
+352685.7851,956702.7685
+352411.2075,954835.6406
+352246.4609,953023.4282
+351312.897,950607.145
+350434.2486,949179.3413
+350873.5728,947806.4532
+351148.1504,945390.17
+351477.6436,942754.2247
+352795.6162,941326.421
+352850.5317,939074.8844
+352411.2075,936768.4323
+351367.8125,934681.6423
+349116.2759,934956.2199
+348347.4585,936438.9391
+348127.7964,938910.1378
+347578.6412,940942.0123
+349390.8535,942534.5626
+349171.1914,944236.9439
+348182.712,945390.17
+347413.8946,946817.9737
+347413.8946,948300.6929
+347798.3033,949563.75
+348237.6275,951046.4692
+348841.6983,953792.2456
+348292.543,955714.289
+348622.0362,957416.6703
+348896.6138,959887.8691
+351093.2349,960107.5312
Index: /anuga_work/development/boxingday08/mesh_polygons/patong.csv
===================================================================
--- /anuga_work/development/boxingday08/mesh_polygons/patong.csv	(revision 5247)
+++ /anuga_work/development/boxingday08/mesh_polygons/patong.csv	(revision 5247)
@@ -0,0 +1,71 @@
+420439.9236,869050.755
+420048.7197,869416.4456
+419751.0646,869816.1539
+419742.5602,870317.9155
+419138.7454,870692.1105
+418560.444,871032.2878
+418152.2312,871457.5094
+418169.2401,872052.8197
+418475.3997,872537.5724
+418730.5327,872852.2364
+419164.2587,872979.8029
+419546.9582,872852.2364
+420065.7286,872554.5813
+420422.9148,872410.0059
+420967.1985,872392.9971
+421307.3758,872384.4926
+421749.6063,872767.1921
+421877.1728,873192.4137
+421987.7304,873507.0777
+422157.8191,873991.8304
+422225.8545,874374.5299
+421919.695,874799.7515
+421256.3492,874833.7693
+420839.632,874833.7693
+420380.3926,874927.318
+420125.2596,875173.9466
+419819.1001,875505.6195
+419334.3474,875548.1416
+418807.0726,875624.6815
+418424.3731,875786.2657
+418356.3376,876219.9918
+419155.7543,876279.5228
+419649.0114,876466.6204
+420057.2242,876653.7179
+420431.4192,876721.7533
+420678.0478,876755.7711
+420907.6675,876347.5583
+421162.8004,876271.0184
+421400.9246,875981.8677
+421605.0309,875820.2835
+421936.7038,875743.7436
+422293.89,875854.3012
+422549.023,875990.3721
+422948.7313,875998.8766
+423127.3244,875701.2214
+423161.3421,875386.5574
+423212.3687,875054.8845
+423271.8998,874714.7072
+423365.4485,874561.6274
+423765.1569,874561.6274
+424079.8209,874527.6097
+424394.4849,874391.5388
+424598.5913,874246.9634
+424768.6799,873728.193
+424700.6445,873507.0777
+424658.1223,873285.9625
+424505.0425,872971.2985
+424266.9184,872665.1389
+424054.3076,872563.0857
+423969.2632,872197.3951
+423841.6968,871874.2267
+423697.1214,871610.5892
+423663.1037,871338.4474
+423518.5283,871142.8454
+423169.8466,870964.2523
+422812.6604,870845.1903
+422361.9255,870751.6415
+422038.757,870445.4819
+421775.1196,870190.349
+421247.8448,869731.1096
+420839.632,869365.419
Index: /anuga_work/development/boxingday08/mesh_polygons/patong_bay_polygon.csv
===================================================================
--- /anuga_work/development/boxingday08/mesh_polygons/patong_bay_polygon.csv	(revision 5247)
+++ /anuga_work/development/boxingday08/mesh_polygons/patong_bay_polygon.csv	(revision 5247)
@@ -0,0 +1,4 @@
+417349.959414, 870580.352912
+425645.571055, 870566.224816
+425655.976923, 876993.540904
+417361.526500, 877007.770637
Index: /anuga_work/development/boxingday08/mesh_polygons/south_island.csv
===================================================================
--- /anuga_work/development/boxingday08/mesh_polygons/south_island.csv	(revision 5247)
+++ /anuga_work/development/boxingday08/mesh_polygons/south_island.csv	(revision 5247)
@@ -0,0 +1,11 @@
+430871.7664,841894.9953
+431819.0592,841894.9953
+432189.7391,841359.5689
+432107.3658,839876.8497
+430995.3263,839053.1168
+430501.0866,838229.3839
+428730.0609,837899.8907
+428647.6876,838641.2504
+428647.6876,840165.1562
+429512.6071,841194.8224
+430089.2202,842142.1152
Index: /anuga_work/development/boxingday08/mesh_polygons/south_island2.csv
===================================================================
--- /anuga_work/development/boxingday08/mesh_polygons/south_island2.csv	(revision 5247)
+++ /anuga_work/development/boxingday08/mesh_polygons/south_island2.csv	(revision 5247)
@@ -0,0 +1,16 @@
+427596.0683,831294.4916
+427346.0282,830194.3151
+426845.948,829269.1667
+426695.9239,828494.0423
+426070.8236,827443.8738
+425495.7313,826793.7695
+424745.611,826443.7134
+424370.5508,826718.7575
+424170.5187,827118.8217
+424695.603,827793.93
+424970.6471,828419.0303
+425120.6712,829294.1707
+425270.6952,829969.279
+425895.7955,830969.4395
+426445.8838,831819.5758
+427196.0041,832244.644
Index: /anuga_work/development/boxingday08/project.py
===================================================================
--- /anuga_work/development/boxingday08/project.py	(revision 5247)
+++ /anuga_work/development/boxingday08/project.py	(revision 5247)
@@ -0,0 +1,86 @@
+"""Common filenames and locations for topographic data, meshes and outputs.
+Author: John Jakeman, The Australian National University (2008)
+"""
+
+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, \
+                                    polygon_area, is_inside_polygon
+from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
+from time import localtime, strftime
+
+def convert_arcgis_latlon_list_to_utm(points):
+     #Used because arc gis produced csv files put lat lon in
+     #reverse order to those accpeted by convert_latlon_to_utm()
+     
+     from anuga.coordinate_transforms.geo_reference import Geo_reference
+     from anuga.coordinate_transforms.redfearn import redfearn
+     old_geo = Geo_reference() 
+     utm_points = []
+     for point in points:
+         zone, easting, northing = redfearn(float(point[1]),float(point[0]))
+         new_geo = Geo_reference(zone)
+         old_geo.reconcile_zones(new_geo) 
+         utm_points.append([easting,northing])
+
+     return utm_points, old_geo.get_zone()
+
+###############################
+# Domain definitions
+###############################
+
+# Bathymetry and topography filenames
+dir='data'
+
+time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
+poor_combined_dir_name = dir+sep+'poor_phuket_bathymetry'
+good_combined_dir_name = dir+sep+'good_phuket_bathymetry'
+meshname = 'phuket_mesh'
+mesh_elevname = 'phuket_mesh_elev'
+boundary_most_in = dir+sep+'out'
+boundary_most_out = dir+sep+'most_boundary_condition'
+
+dir='mesh_polygons'
+#extent = 'extent.csv'
+extent = 'boundary_new_pts.csv'
+# bounding polygon for study area
+bounding_polygon = read_polygon(dir+sep+extent)
+patong_bay_polygon = read_polygon(dir+sep+'patong_bay_polygon.csv')
+
+###############################
+# Interior region definitions
+###############################
+
+island_south = read_polygon(dir+sep+'south_island.csv')
+island_south2 = read_polygon(dir+sep+'south_island2.csv')
+island_north = read_polygon(dir+sep+'north_island.csv')
+bay = read_polygon(dir+sep+'bay.csv')
+patong = read_polygon(dir+sep+'patong.csv')
+contour20m = read_polygon(dir+sep+'20m.csv')
+
+
+
+#max_extent_of MOST file
+#east  = 94.93325
+#west  = 80.0003
+#north = 10.0000
+#south = -4.9330
+
+#max_extent_of MOST file should be
+#east  = 104.0 
+#west  = 80.0
+#north = 10.0
+#south = -5.0
+
+#max_extent of ol simulation is why does this work and not below
+#7.40911081272 8.71484358635 99.1683687224 97.513856322
+
+#max extent of clipped MOST file
+east  = 100
+west  = 96.0
+north = 9.0
+south = 6.0
+
+
Index: /anuga_work/development/boxingday08/run_boxingday.py
===================================================================
--- /anuga_work/development/boxingday08/run_boxingday.py	(revision 5247)
+++ /anuga_work/development/boxingday08/run_boxingday.py	(revision 5247)
@@ -0,0 +1,244 @@
+"""Script for running the 2004 boxing Day tsunami inundation scenario for
+Phuket, Thailand.
+
+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 directory named after the scenario, i.e
+slide or fixed_wave.
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a tsunami wave generated by a submarine mass failure.
+
+Author: John Jakeman, The Australian National University (2008)
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import time
+import sys
+from time import localtime, strftime
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import Time_boundary
+from anuga.shallow_water import File_boundary
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, ferret2sww
+from anuga.shallow_water.data_manager import dem2pts
+from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
+from anuga.utilities.polygon import number_mesh_triangles
+from anuga.fit_interpolate.fit import fit_to_mesh_file
+from anuga.caching import cache
+from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
+
+# Application specific imports
+import project                 # Definition of file names and polygons
+
+
+#------------------------------------------------------------------------------
+# Define scenario as either slide or fixed_wave.
+#------------------------------------------------------------------------------
+
+#scenario = 'poor_simulation'
+scenario = 'good_simulation'
+
+if os.access(scenario, os.F_OK) == 0:
+    os.mkdir(scenario)
+
+timestamp = strftime('%Y%m%d_%H%M%S',localtime())
+basename = scenario + 'source'
+
+
+#------------------------------------------------------------------------------
+# Preparation of topographic data
+# Convert ASC 2 DEM 2 PTS using source data and store result in source data
+#------------------------------------------------------------------------------
+
+if scenario == 'good_simualation':
+     msg = 'Must use combine_good_data to create bathymetry .pts file'
+     assert os.path.exists(project.good_combined_dir_name+'.pts'), msg
+
+#------------------------------------------------------------------------------
+# 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
+#------------------------------------------------------------------------------
+#extent_res     = 10000000.0
+#contour20m_res = 50000.0
+#island_res     = 5000.0
+#bay_res        = 2000.0
+#patong_res     = 100.0
+
+extent_res     = 10000000.0
+contour20m_res = 1000000.0
+island_res     = 1000000.0
+bay_res        = 1000000.0
+patong_res     = 1000000.0
+
+
+# make polygon that contains land that does not affect result.
+
+interior_regions=None
+#interior_regions = [[project.patong, patong_res],
+#                    [project.bay, bay_res],
+#                    [project.contour20m, contour20m_res],
+#                    [project.island_north, island_res],
+#                    [project.island_south, island_res],
+#                    [project.island_south2, island_res]]
+
+# need to include runup_res regions
+
+# There are 9 boundaries. FIX ME: CHECK TAGS ARE CORRECT. I.E. THAT I HAVE
+#CHOSEN CORECT ORDER AS SEEEN BY ANUGA. IS THERE A WAY TO DO This
+print 'must determine which boundaries are which'
+
+#boundary_tags={'ocean': [1,2], 'otherocean': [0,3, 7, 8], 'land': [5, 6], 'both': [4]}
+
+boundary_tags =  {'most':[10,11,12,13,14,15],
+                  'land':[0,1,2,3,4,5,6,7,8],
+                  'shallow':[9,16]}
+
+#trigs_min = number_mesh_triangles(interior_regions, project.bounding_polygon,extent_res)
+#print 'Minimum number of traingles ', trigs_min
+
+# filenames
+meshname = project.meshname + '.tsh'
+mesh_elevname = project.mesh_elevname  + '.tsh'
+
+print 'start create mesh from regions'
+
+_ = cache(create_mesh_from_regions,
+          project.bounding_polygon,
+          {'boundary_tags': boundary_tags,
+           'maximum_triangle_area': extent_res,
+           'filename': meshname,           
+           'interior_regions': interior_regions},
+          verbose = True,
+          dependencies = ['run_boxingday.py']
+          #, evaluate=True
+          )
+
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+print 'Converting mesh to domain'
+
+#domain = Domain(meshname, use_cache=False, verbose=True)
+
+domain = cache(pmesh_to_domain_instance,
+               (meshname, Domain),
+               dependencies = [meshname])
+
+print 'Number of triangles = ', len(domain)
+print 'The extent is ', domain.get_extent()
+print domain.statistics()
+
+domain.set_name(basename+'friction'+timestamp) 
+domain.set_datadir(scenario)
+domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+domain.set_minimum_storable_height(0.01)
+
+domain.beta_h = 0
+domain.tight_slope_limiters = 1
+domain.set_default_order(2)
+print 'domain.tight_slope_limiters', domain.tight_slope_limiters
+
+domain.points_file_block_line_size = 50000
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+     
+tide = 0.0
+domain.set_quantity('stage', tide)
+domain.set_quantity('friction', 0.01) 
+
+if scenario == 'poor_simulation':
+     domain.set_quantity('elevation', 
+                         filename=project.poor_combined_dir_name + '.pts',
+                         use_cache=True,
+                         verbose=True,
+                         alpha=0.1)
+
+if scenario == 'good_simulation':
+     domain.set_quantity('elevation', 
+                         2.0,#filename=project.good_combined_dir_name + '.pts',
+                         use_cache=True,
+                         verbose=True,
+                         alpha=0.1)
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+
+print 'start ferret2sww'
+south = project.south 
+north = project.north
+west = project.west
+east = project.east
+
+# Note only need to do when an SWW file for the MOST boundary doesn't exist or the bounding polygon has changed
+
+if os.path.exists(project.boundary_most_out+'.sww'):
+    print 'sww boundary file already created ensure boundary polygon has not changed'
+else:
+     cache(ferret2sww,
+           (project.boundary_most_in,
+            project.boundary_most_out), 
+           {'verbose': True,
+            'minlat': south,
+            'maxlat': north,
+            'minlon': west,
+            'maxlon': east,
+            #'origin': project.mesh_origin,
+            'origin': domain.geo_reference.get_origin(),
+            'mean_stage': tide,
+            'zscale': 1.0,                 #Enhance tsunami
+            'fail_on_NaN': False,
+            'inverted_bathymetry': True},
+           #evaluate = True,
+           verbose = True,
+           dependencies = 'run_boxingday.py')
+
+print 'Available boundary tags', domain.get_boundary_tags()
+Bf = File_boundary(project.boundary_most_out+'.sww', 
+                    domain, verbose = True)
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([tide,0,0])
+
+#domain.set_boundary({'ocean': Bf,
+#                     'otherocean': Bd,
+#                     'land': Br,
+#                     'both': Bd})
+
+domain.set_boundary({'most': Bf,
+                     'land': Br,
+                     'shallow': Bd})
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+import time
+t0 = time.time()
+
+#from Numeric import allclose
+#from anuga.abstract_2d_finite_volumes.quantity import Quantity
+
+# Add new loop that uses larger yieldstep until wave first reaches a point of
+# the ANUGA boundary. Or find a way to clip MOST sww boundary file to only
+# start when boundary stage first becomes non-zero.
+
+for t in domain.evolve(yieldstep = 60, finaltime = 18000, 
+                       skip_initial_step = False):
+     domain.write_time()
+     
+# Add new loop that uses smaller yieldstep when wave approaches shore?
+
+            
+print 'That took %.2f seconds' %(time.time()-t0)
