Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/project.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/project.py	(revision 9321)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/project.py	(revision 9322)
@@ -40,5 +40,8 @@
 yieldStep=20.0
 
-
+## Spatial Information
+#
+# Give UTM Zone as an integer (negative for Southern Hemisphere)
+utm_zone = 54 
 
 #####################################################################
@@ -141,6 +144,7 @@
 # pixle size (m) for output tifs
 output_tif_cellsize = 250.
-# EPSG CODE (projection information) for outputs tifs
-EPSG_CODE = 32647
+# Polygon for output tif [either a filename, or if None the bounding_polygon is used]
+output_tif_bounding_polygon = None 
+
 
 ## Storage of max quantities
@@ -214,4 +218,10 @@
             [bounding_polygon, interior_regions, riverWalls, \
              breakLines, RegionPtAreas, default_res, BoundaryInfo]) ).hexdigest()
+
+# Fix the output tif bounding polygon
+if output_tif_bounding_polygon is None:
+    output_tif_bounding_polygon = bounding_polygon
+else:
+    output_tif_bounding_polygon = su.read_polygon(output_tif_bounding_polygon)
 
 # Set up directories etc 
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/raster_outputs.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/raster_outputs.py	(revision 9321)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/raster_outputs.py	(revision 9322)
@@ -10,4 +10,5 @@
 import numpy
 import gdal
+from pyproj import Proj
 from anuga.utilities import plot_utils as util
 from anuga.utilities import spatialInputUtil as su
@@ -19,5 +20,5 @@
 # 
 
-def make_resampled_elevation(elevation_raster, rasterExtent, CellSize, clipPolygon, clipPolygonLayer, EPSG_CODE, tif_output_dir):
+def make_resampled_elevation(elevation_raster, rasterExtent, CellSize, clipPolygon, clipPolygonLayer, proj4string , tif_output_dir):
     """
         Call gdalwarp to get the resampled elevation data
@@ -42,5 +43,6 @@
     yres=CellSize
 
-    allSRS='EPSG:'+str(EPSG_CODE)
+    #allSRS='EPSG:'+str(EPSG_CODE)
+    allSRS = proj4string
     newDEM=tif_output_dir+'LIDAR_resampled.tif'
     newMask=tif_output_dir+'Mask_resampled.tif'
@@ -66,4 +68,6 @@
 
     return
+
+##############################################################################
   
 def gdal_calc_command(stage, depth, elevation, mask, depth_threshold):
@@ -89,5 +93,5 @@
 def make_me_some_tifs(swwFile,
                     bounding_polygon,
-                    EPSG_CODE,
+                    utm_zone,
                     myTimeStep='collected_max',
                     tif_output_subdir='/TIFS/',
@@ -104,5 +108,5 @@
     - swwFile -- Full path name of sww to read outputs from
     - bounding_polygon -- ANUGA's bounding polygon, or another polygon to clip the rasters to
-    - EPSG_CODE -- EPSG projection info as an integer like 3123 [which is over Manila, Philippines]
+    - utm_zone -- UTM zone as an integer [make negative for Southern Hemisphere]
     - myTimeStep -- option from util.Make_Geotif
        use 'max' to plot the maxima
@@ -125,4 +129,8 @@
     """
 
+    # Convert utm_zone to proj4string
+    p=Proj(proj='utm', south=(utm_zone<0.), zone=abs(utm_zone), ellps='WGS84')
+    proj4string = p.srs
+
     tif_output_dir=os.path.dirname(swwFile)+tif_output_subdir
     # Make the geotifs 
@@ -144,5 +152,5 @@
             tmpArr=maxQs[:,[0, 1, i+2]]
             util.Make_Geotif(tmpArr, output_quantities=[quant+'_MAX'],
-                             CellSize=CellSize, EPSG_CODE=EPSG_CODE, verbose=True,
+                             CellSize=CellSize, proj4string = proj4string , verbose=True,
                              bounding_polygon=bounding_polygon, output_dir=tif_output_dir)
         
@@ -150,5 +158,5 @@
         util.Make_Geotif(swwFile, myTimeStep=[0],
                          output_quantities=['elevation', 'friction'], 
-                         CellSize=CellSize, EPSG_CODE=EPSG_CODE, verbose=True, 
+                         CellSize=CellSize, proj4string = proj4string , verbose=True, 
                          bounding_polygon=bounding_polygon,
                          output_dir=tif_output_dir)
@@ -156,5 +164,5 @@
         util.Make_Geotif(swwFile, myTimeStep=myTimeStep,
                          output_quantities=['depth', 'stage', 'elevation', 'velocity', 'depthIntegratedVelocity', 'friction'], 
-                         CellSize=CellSize, EPSG_CODE=EPSG_CODE, verbose=True, 
+                         CellSize=CellSize, proj4string = proj4string , verbose=True, 
                          bounding_polygon=bounding_polygon,
                          output_dir=tif_output_dir)
@@ -169,5 +177,5 @@
     
     # Make the resampled elevation data 
-    make_resampled_elevation(elevation_raster, rasterExtent, CellSize, clipPolygon, clipPolygonLayer, EPSG_CODE, tif_output_dir)
+    make_resampled_elevation(elevation_raster, rasterExtent, CellSize, clipPolygon, clipPolygonLayer, proj4string, tif_output_dir)
     
     elevation=glob.glob(tif_output_dir+'LIDAR_resampled*')[0]
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/run_model.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/run_model.py	(revision 9321)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/run_model.py	(revision 9322)
@@ -62,8 +62,8 @@
 #
 ####################################################################################
-print 'Making rainfall conditions'
-setup_rainfall.setup_rainfall(domain, project)
-print 'Making inlets '
-setup_inlets.setup_inlets(domain, project)
+# print 'Making rainfall conditions'
+# setup_rainfall.setup_rainfall(domain, project)
+# print 'Making inlets '
+# setup_inlets.setup_inlets(domain, project)
 
 #####################################################################################
@@ -138,5 +138,5 @@
         raster_outputs.make_me_some_tifs(swwFile = './'+project.scenario+'.sww',
                                     bounding_polygon=project.bounding_polygon,
-                                    EPSG_CODE = project.EPSG_CODE,
+                                    utm_zone = project.utm_zone,
                                     CellSize= project.output_tif_cellsize)
     except:
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/setup_boundary_conditions.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/setup_boundary_conditions.py	(revision 9321)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/setup_boundary_conditions.py	(revision 9322)
@@ -49,5 +49,5 @@
         THE USER MUST MODIFY THIS TO SETUP THE BOUNDARIES AS DESIRED
 
-        1) First make sure the right boundary conditions are created
+        1) First make sure the required boundary conditions are created
            (see examples below, modify or add new ones as required)
 
Index: /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/setup_initial_conditions.py
===================================================================
--- /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/setup_initial_conditions.py	(revision 9321)
+++ /trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/setup_initial_conditions.py	(revision 9322)
@@ -63,4 +63,31 @@
     #START EDITING HERE
     #@@@@@@@@@@@@@@@@@#
+
+
+    ##############################################################################
+    #
+    # PRELIMINARY DEFINITIONS
+    # 
+    ##############################################################################
+
+    # If the earthquake slip is given then convert to deformation here
+    # So that we can add it to both elevation AND stage
+    make_earthquake_deformation_from_source_model=False
+
+    if make_earthquake_deformation_from_source_model:
+        source_model='../anuga/source/XXXXX.csv'
+        # Make an earthquake source model
+        eq_function = okada_tsunami.earthquake_source(source=None,
+                                                  filename=source_model,
+                                                  domain=domain,
+                                                  lon_lat_degrees=True,
+                                                  lon_before_lat=False,
+                                                  utm_zone=project.utm_zone,
+                                                  verbose=False)
+    else:
+        # Make a function that is zero everywhere
+        # If we pass this to add_quantity, nothing should happen
+        def eq_function(x,y):
+            return x*0.
 
 
@@ -80,5 +107,5 @@
     elevation_default = numpy.genfromtxt('../anuga/topographies/andaman_3s_grid_Clip2_Project.txt', delimiter=",",skip_header=1)
 
-
+    # Set_quantity data
     elevation_Data=[
                     ['../anuga/polygons/patong_10m.txt', elevation_patong_10m],
@@ -88,4 +115,10 @@
                    ]
 
+    # Define a function for add_quantity in the same way
+    # [ ['All', 0.] ] will do nothing
+    elevation_additions = [ 
+                            ['All', eq_function] 
+                          ]
+
    
     ##############################################################################
@@ -96,7 +129,13 @@
     
     default_friction=0.025
+    # Set_quantity data
     friction_Data =[
                     ['All', default_friction]
                    ]
+   
+    # Add_quantity data 
+    friction_additions = [ 
+                          ['All', 0.] 
+                         ]
     
     ###############################################################################
@@ -106,8 +145,14 @@
     ###############################################################################
 
+    # Set_quantity data
     stage_Data = [
                   ['All', project.tide]
                  ]
 
+    
+    # Add_quantity data 
+    stage_additions = [ 
+                       ['All', eq_function] 
+                      ]
     ###############################################################################
     #
@@ -117,11 +162,23 @@
     ###############################################################################
     
+    # Set_quantity data
     xmomentum_Data = [
                       ['All', 0.]
                      ]
-
+   
+    # Add_quantity data
+    xmomentum_additions = [ 
+                            ['All', 0.] 
+                          ]
+
+    # Set_quantity data
     ymomentum_Data = [
                       ['All', 0.]
                      ]
+    
+    # Add_quantity data 
+    ymomentum_additions = [ 
+                            ['All', 0.] 
+                          ]
 
     #@@@@@@@@@@@@@@@@@#
@@ -140,16 +197,26 @@
     elevFun=qs.composite_quantity_setting_function(elevation_Data, domain)
     domain.set_quantity('elevation', elevFun, location='vertices') 
+    elevAddFun=qs.composite_quantity_setting_function(elevation_additions, domain)
+    domain.add_quantity('elevation', elevAddFun) 
     
     frictionFun=qs.composite_quantity_setting_function(friction_Data, domain)
     domain.set_quantity('friction',  frictionFun,  location='centroids') 
+    frictionAddFun=qs.composite_quantity_setting_function(friction_additions, domain)
+    domain.add_quantity('friction', frictionAddFun) 
     
     stageFun=qs.composite_quantity_setting_function(stage_Data, domain)
     domain.set_quantity('stage',  stageFun,  location='centroids') 
+    stageAddFun=qs.composite_quantity_setting_function(stage_additions, domain)
+    domain.add_quantity('stage', stageAddFun) 
 
     xmomentumFun=qs.composite_quantity_setting_function(xmomentum_Data, domain)
     domain.set_quantity('xmomentum',  xmomentumFun,  location='centroids') 
+    xmomentumAddFun=qs.composite_quantity_setting_function(xmomentum_additions, domain)
+    domain.add_quantity('xmomentum', xmomentumAddFun) 
     
     ymomentumFun=qs.composite_quantity_setting_function(ymomentum_Data, domain)
     domain.set_quantity('ymomentum',  ymomentumFun,  location='centroids') 
+    ymomentumAddFun=qs.composite_quantity_setting_function(ymomentum_additions, domain)
+    domain.add_quantity('ymomentum', ymomentumAddFun) 
 
     ################################################################################
