Index: /trunk/anuga_core/source/anuga/__init__.py
===================================================================
--- /trunk/anuga_core/source/anuga/__init__.py	(revision 7868)
+++ /trunk/anuga_core/source/anuga/__init__.py	(revision 7869)
@@ -93,5 +93,5 @@
 # Forcing
 #-----------------------------
-from anuga.shallow_water.forcing import Inflow, Rainfall
+from anuga.shallow_water.forcing import Inflow, Rainfall, Wind_stress
 
 #-----------------------------
Index: /trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_basic.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_basic.py	(revision 7868)
+++ /trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_basic.py	(revision 7869)
@@ -133,4 +133,5 @@
 
     def test_rotate(self):
+        from anuga.shallow_water.shallow_water_ext import rotate
         normal = num.array([0.0, -1.0])
 
@@ -492,6 +493,6 @@
         # Setup boundary conditions
         #--------------------------------------------------------------
-        Br = Reflective_boundary(domain)                       # Reflective wall
-        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
+        Br = anuga.Reflective_boundary(domain)                 # Reflective wall
+        Bd = anuga.Dirichlet_boundary([final_runup_height, 0, 0])# Steady inflow
 
         # All reflective to begin with (still water)
@@ -598,7 +599,7 @@
         from anuga.abstract_2d_finite_volumes.mesh_factory \
                 import rectangular_cross
-        from data_manager import get_maximum_inundation_elevation
-        from data_manager import get_maximum_inundation_location
-        from data_manager import get_maximum_inundation_data
+        from anuga.shallow_water.sww_interrogate import \
+            get_maximum_inundation_elevation, get_maximum_inundation_location, \
+            get_maximum_inundation_data
 
         initial_runup_height = -0.4
@@ -632,6 +633,6 @@
         # Setup boundary conditions
         #--------------------------------------------------------------
-        Br = Reflective_boundary(domain)                       # Reflective wall
-        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
+        Br = anuga.Reflective_boundary(domain)                 # Reflective wall
+        Bd = anuga.Dirichlet_boundary([final_runup_height, 0, 0])# Steady inflow
 
         # All reflective to begin with (still water)
@@ -787,16 +788,9 @@
         """
 
-        from anuga.pmesh.mesh_interface import create_mesh_from_regions
-        from anuga.abstract_2d_finite_volumes.mesh_factory \
-                import rectangular_cross
-        from anuga.shallow_water import Domain
-        from anuga.shallow_water import Reflective_boundary
-        from anuga.shallow_water import Dirichlet_boundary
-
         #-----------------------------------------------------------------
         # Setup computational domain
         #-----------------------------------------------------------------
-        points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
-        domain = Domain(points, vertices, boundary) # Create domain
+        points, vertices, boundary = anuga.rectangular_cross(10, 10)# Basic mesh
+        domain = anuga.Domain(points, vertices, boundary) # Create domain
         domain.set_default_order(2)
         domain.set_quantities_to_be_stored(None)
@@ -816,6 +810,6 @@
         # Setup boundary conditions
         #----------------------------------------------------------------
-        Br = Reflective_boundary(domain)           # Solid reflective wall
-        Bd = Dirichlet_boundary([-0.2, 0., 0.])    # Constant boundary values
+        Br = anuga.Reflective_boundary(domain)           # Solid reflective wall
+        Bd = anuga.Dirichlet_boundary([-0.2, 0., 0.])    # Constant
         domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
 
@@ -879,5 +873,5 @@
         vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
 
-        domain = Domain(points, vertices)
+        domain = anuga.Domain(points, vertices)
 
         #Set up for a gradient of (3,0) at mid triangle (bce)
@@ -898,5 +892,5 @@
         initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
 
-        domain.set_boundary({'exterior': Reflective_boundary(domain)})
+        domain.set_boundary({'exterior': anuga.Reflective_boundary(domain)})
 
         domain.optimise_dry_cells = True
@@ -916,16 +910,15 @@
 
     def test_second_order_flat_bed_onestep(self):
-        from mesh_factory import rectangular
 
         #Create basic mesh
-        points, vertices, boundary = rectangular(6, 6)
+        points, vertices, boundary = anuga.rectangular(6, 6)
 
         #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
+        domain = anuga.Domain(points, vertices, boundary)
         domain.set_default_order(2)
 
         # Boundary conditions
-        Br = Reflective_boundary(domain)
-        Bd = Dirichlet_boundary([0.1, 0., 0.])
+        Br = anuga.Reflective_boundary(domain)
+        Bd = anuga.Dirichlet_boundary([0.1, 0., 0.])
         domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
 
@@ -1027,11 +1020,11 @@
 
         # Create shallow water domain
-        domain = Domain(points, vertices, boundary)
+        domain = anuga.Domain(points, vertices, boundary)
         domain.smooth = False
         domain.default_order = 2
 
         # Boundary conditions
-        Br = Reflective_boundary(domain)
-        Bd = Dirichlet_boundary([0.1, 0., 0.])
+        Br = anuga.Reflective_boundary(domain)
+        Bd = anuga.Dirichlet_boundary([0.1, 0., 0.])
         domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
 
@@ -1063,6 +1056,6 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain)
-        Bd = Dirichlet_boundary([0.2, 0., 0.])
+        Br = anuga.Reflective_boundary(domain)
+        Bd = anuga.Dirichlet_boundary([0.2, 0., 0.])
 
         domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
@@ -1101,6 +1094,6 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain)
-        Bd = Dirichlet_boundary([0.2, 0., 0.])
+        Br = anuga.Reflective_boundary(domain)
+        Bd = anuga.Dirichlet_boundary([0.2, 0., 0.])
 
         domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
@@ -1156,6 +1149,6 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain)
-        Bd = Dirichlet_boundary([0.2, 0., 0.])
+        Br = anuga.Reflective_boundary(domain)
+        Bd = anuga.Dirichlet_boundary([0.2, 0., 0.])
 
         domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
@@ -1201,6 +1194,6 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain)
-        Bd = Dirichlet_boundary([0.2, 0., 0.])
+        Br = anuga.Reflective_boundary(domain)
+        Bd = anuga.Dirichlet_boundary([0.2, 0., 0.])
 
         domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
@@ -1348,5 +1341,5 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
 
@@ -1512,5 +1505,5 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain)
+        Br = anuga.Reflective_boundary(domain)
         domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
 
@@ -1556,6 +1549,6 @@
         domain.set_quantity('elevation', Z)
 
-        Br = Reflective_boundary(domain)
-        Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
+        Br = anuga.Reflective_boundary(domain)
+        Bd = anuga.Dirichlet_boundary([inflow_stage, 0.0, 0.0])
         domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
 
@@ -1614,5 +1607,4 @@
         import time, os
         from Scientific.IO.NetCDF import NetCDFFile
-        from data_manager import extent_sww
         from mesh_factory import rectangular_cross
 
@@ -1728,7 +1720,7 @@
 
          tags = {}
-         b1 =  Dirichlet_boundary(conserved_quantities = num.array([0.0]))
-         b2 =  Dirichlet_boundary(conserved_quantities = num.array([1.0]))
-         b3 =  Dirichlet_boundary(conserved_quantities = num.array([2.0]))
+         b1 =  anuga.Dirichlet_boundary(conserved_quantities = num.array([0.0]))
+         b2 =  anuga.Dirichlet_boundary(conserved_quantities = num.array([1.0]))
+         b3 =  anuga.Dirichlet_boundary(conserved_quantities = num.array([2.0]))
          tags["1"] = b1
          tags["2"] = b2
@@ -1909,9 +1901,4 @@
 
         verbose = False
-
-        from anuga.shallow_water import Domain
-        from anuga.pmesh.mesh_interface import create_mesh_from_regions
-        from anuga.geospatial_data.geospatial_data import Geospatial_data
-
         
         # Get path where this test is run
@@ -1965,5 +1952,5 @@
 
         meshname = os.path.join(path, 'offending_mesh.msh')
-        create_mesh_from_regions(bounding_polygon,
+        anuga.create_mesh_from_regions(bounding_polygon,
                                  boundary_tags={'south': [0], 'east': [1],
                                                 'north': [2], 'west': [3]},
@@ -1974,5 +1961,5 @@
                                  verbose=verbose)
 
-        domain = Domain(meshname, use_cache=False, verbose=verbose)
+        domain = anuga.Domain(meshname, use_cache=False, verbose=verbose)
 
         #----------------------------------------------------------------------
@@ -2000,12 +1987,4 @@
             os.remove(points_file)            
         
-        #finally:
-            # Cleanup regardless
-            #FIXME(Ole): Finally does not work like this in python2.4 
-            #FIXME(Ole): Reinstate this when Python2.4 is out of the way
-            #FIXME(Ole): Python 2.6 apparently introduces something called 'with'            
-            #os.remove(meshname)
-            #os.remove(points_file)
-
 
     def test_fitting_example_that_crashed_2(self):
@@ -2022,8 +2001,4 @@
 
         verbose = False        
-
-        from anuga.shallow_water import Domain
-        from anuga.pmesh.mesh_interface import create_mesh_from_regions
-        from anuga.geospatial_data import Geospatial_data
         
         # Get path where this test is run
@@ -2039,5 +2014,5 @@
         bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
 
-        create_mesh_from_regions(bounding_polygon,
+        anuga.create_mesh_from_regions(bounding_polygon,
                                  boundary_tags={'south': [0], 
                                                 'east': [1], 
@@ -2049,5 +2024,5 @@
                                  verbose=verbose)
 
-        domain = Domain(meshname, use_cache=True, verbose=verbose)
+        domain = anuga.Domain(meshname, use_cache=True, verbose=verbose)
         
         # Large test set revealed one problem
@@ -2131,7 +2106,7 @@
         # Setup boundary conditions
         #------------------------------------------------------------------
-        Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
-        Br = Reflective_boundary(domain)              # Solid reflective wall
-        Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
+        Bi = anuga.Dirichlet_boundary([0.4, 0, 0])          # Inflow
+        Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
+        Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # Outflow
 
         domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
Index: /trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_boundary_condition.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_boundary_condition.py	(revision 7868)
+++ /trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_boundary_condition.py	(revision 7869)
@@ -5,4 +5,6 @@
 from math import pi, sqrt
 import tempfile
+
+import anuga
 
 from anuga.config import g, epsilon
@@ -60,7 +62,7 @@
                                           [30,30,30], [40,40,40]])
 
-        D = Dirichlet_boundary([5, 2, 1])
-        T = Transmissive_boundary(domain)
-        R = Reflective_boundary(domain)
+        D = anuga.Dirichlet_boundary([5, 2, 1])
+        T = anuga.Transmissive_boundary(domain)
+        R = anuga.Reflective_boundary(domain)
         domain.set_boundary({'First': D, 'Second': T, 'Third': R})
 
@@ -140,5 +142,5 @@
                     (0, 1): 'Internal'}
 
-        domain = Domain(points, vertices, boundary)
+        domain = anuga.Domain(points, vertices, boundary)
         domain.check_integrity()
 
@@ -150,7 +152,7 @@
                                           [30,30,30], [40,40,40]])
 
-        D = Dirichlet_boundary([5, 2, 1])
-        T = Transmissive_boundary(domain)
-        R = Reflective_boundary(domain)
+        D = anuga.Dirichlet_boundary([5, 2, 1])
+        T = anuga.Transmissive_boundary(domain)
+        R = anuga.Reflective_boundary(domain)
         domain.set_boundary({'First': D, 'Second': T,
                              'Third': R, 'Internal': None})
@@ -192,7 +194,7 @@
                                           [30,30,30], [40,40,40]])
 
-        D = Dirichlet_boundary([5, 2, 1])
-        T = Transmissive_stage_zero_momentum_boundary(domain)
-        R = Reflective_boundary(domain)
+        D = anuga.Dirichlet_boundary([5, 2, 1])
+        T = anuga.Transmissive_stage_zero_momentum_boundary(domain)
+        R = anuga.Reflective_boundary(domain)
         domain.set_boundary({'First': D, 'Second': T, 'Third': R})
 
@@ -270,7 +272,7 @@
         #--------------------------------------------------------------
         # Time dependent boundary
-        Bt = Time_boundary(domain=domain, f=lambda t: [t, 0.0, 0.0])
-
-        Br = Reflective_boundary(domain)              # Reflective wall
+        Bt = anuga.Time_boundary(domain=domain, f=lambda t: [t, 0.0, 0.0])
+
+        Br = anuga.Reflective_boundary(domain)              # Reflective wall
 
         domain.set_boundary({'left': Bt, 'right': Br, 'top': Br, 'bottom': Br})
@@ -326,6 +328,6 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain1)
-        Bd = Dirichlet_boundary([0.3,0,0])
+        Br = anuga.Reflective_boundary(domain1)
+        Bd = anuga.Dirichlet_boundary([0.3,0,0])
         domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
 
@@ -370,6 +372,6 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain2)
-        Bf = Field_boundary(domain1.get_name() + '.sww', domain2)
+        Br = anuga.Reflective_boundary(domain2)
+        Bf = anuga.Field_boundary(domain1.get_name() + '.sww', domain2)
         domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
         domain2.check_integrity()
@@ -433,6 +435,6 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain1)
-        Bd = Dirichlet_boundary([0.3,0,0])
+        Br = anuga.Reflective_boundary(domain1)
+        Bd = anuga.Dirichlet_boundary([0.3,0,0])
         domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
 
@@ -498,6 +500,6 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain2)
-        Bf = Field_boundary(domain1.get_name() + '.sww',
+        Br = anuga.Reflective_boundary(domain2)
+        Bf = anuga.Field_boundary(domain1.get_name() + '.sww',
                             domain2, verbose=False)
         domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
@@ -605,6 +607,6 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain1)
-        Bd = Dirichlet_boundary([0.3, 0, 0])
+        Br = anuga.Reflective_boundary(domain1)
+        Bd = anuga.Dirichlet_boundary([0.3, 0, 0])
         domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
 
@@ -671,6 +673,6 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain2)
-        Bf = Field_boundary(domain1.get_name() + '.sww',
+        Br = anuga.Reflective_boundary(domain2)
+        Bf = anuga.Field_boundary(domain1.get_name() + '.sww',
                             domain2, mean_stage=mean_stage, verbose=False)
 
@@ -761,5 +763,5 @@
 
         # Create shallow water domain
-        domain1 = Domain(points, vertices, boundary)
+        domain1 = anuga.Domain(points, vertices, boundary)
 
         domain1.reduction = mean
@@ -784,6 +786,6 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain1)
-        Bd = Dirichlet_boundary([0.3, 0, 0])
+        Br = anuga.Reflective_boundary(domain1)
+        Bd = anuga.Dirichlet_boundary([0.3, 0, 0])
         domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
 
@@ -846,6 +848,6 @@
 
         # Boundary conditions
-        Br = Reflective_boundary(domain2)
-        Bf = Field_boundary(domain1.get_name() + '.sww',
+        Br = anuga.Reflective_boundary(domain2)
+        Bf = anuga.Field_boundary(domain1.get_name() + '.sww',
                             domain2, mean_stage=1, verbose=False)
 
Index: /trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_forcing_terms.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_forcing_terms.py	(revision 7868)
+++ /trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_forcing_terms.py	(revision 7869)
@@ -552,5 +552,5 @@
         os.remove(filename + '.tms')
 
-        W = Wind_stress(F)
+        W = anuga.Wind_stress(F)
 
         domain.forcing_terms = []
@@ -1163,5 +1163,5 @@
         # on a circle affecting triangles #0 and #1 (bac and bce)
         domain.forcing_terms = []
-        I = Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1)
+        I = anuga.Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1)
         domain.forcing_terms.append(I)
         
@@ -1290,6 +1290,4 @@
         # for a range of stage values
         for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
-            print stage
-
             domain.time = 0.0
             domain.set_quantity('stage', stage)
@@ -1302,5 +1300,4 @@
             for t in domain.evolve(yieldstep=dt, finaltime=5.0):
                 volume = domain.quantities['stage'].get_integral()
-                print t, volume, predicted_volume
                 assert num.allclose (volume, predicted_volume)
                 predicted_volume = predicted_volume - 2.0/pi/100/dt # Why 100?
@@ -1309,6 +1306,4 @@
         # range of stage values
         for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
-            print stage
-
             domain.time = 0.0
             domain.set_quantity('stage', stage)
@@ -1323,6 +1318,4 @@
             for t in domain.evolve(yieldstep=dt, finaltime=5.0):
                 volume = domain.quantities['stage'].get_integral()
-
-                print t, volume
                 assert num.allclose(volume, initial_volume)
 
@@ -1403,5 +1396,5 @@
                     domain.forcing_terms.append(fixed_inflow)
                 
-                ref_flow = fixed_inflow.rate*number_of_inflowsg
+                ref_flow = fixed_inflow.rate*number_of_inflows
 
                 # Compute normal depth on plane using Mannings equation
@@ -1417,5 +1410,5 @@
                 #--------------------------------------------------------------
 
-                Br = Reflective_boundary(domain)
+                Br = anuga.Reflective_boundary(domain)
                 
                 # Define downstream boundary based on predicted depth
@@ -1423,9 +1416,6 @@
                     return (-slope*length) + normal_depth
                 
-                Bt = Transmissive_momentum_set_stage_boundary(domain=domain,
-                                                              function=normal_depth_stage_downstream)
-                
-
-                
+                Bt = anuga.Transmissive_momentum_set_stage_boundary(
+                        domain=domain, function=normal_depth_stage_downstream)
 
                 domain.set_boundary({'left': Br,
@@ -1433,5 +1423,4 @@
                                      'top': Br,
                                      'bottom': Br})
-
 
 
@@ -1617,158 +1606,5 @@
 
         
-        
-    def Xtest_friction_dependent_flow_using_flowline(self):
-        """test_friction_dependent_flow_using_flowline
-        
-        Test the internal flow (using flowline) as a function of
-        different values of Mannings n and different slopes.
-        
-        Flow is applied in the form of boundary conditions with fixed momentum.
-        """
-
-        verbose = True
-
-        #----------------------------------------------------------------------
-        # Import necessary modules
-        #----------------------------------------------------------------------
-
-        from anuga.abstract_2d_finite_volumes.mesh_factory \
-                import rectangular_cross
-        from anuga.shallow_water import Domain
-        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
-        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
-        from anuga.shallow_water.forcing import Inflow
-        from anuga.shallow_water.data_manager \
-                import get_flow_through_cross_section
-        from anuga.abstract_2d_finite_volumes.util \
-                import sww2csv_gauges, csv2timeseries_graphs
-
-
-        #----------------------------------------------------------------------
-        # Setup computational domain
-        #----------------------------------------------------------------------
-
-        finaltime = 1000.0
-
-        length = 300.
-        width  = 20.
-        dx = dy = 5       # Resolution: of grid on both axes
-        
-        # Input parameters
-        uh = 1.0
-        vh = 0.0
-        d = 1.0
-        
-        ref_flow = uh*d*width # 20 m^3/s in the x direction across entire domain
-
-        points, vertices, boundary = rectangular_cross(int(length/dx),
-                                                       int(width/dy),
-                                                       len1=length,
-                                                       len2=width)
-
-        for mannings_n in [0.035]:          #[0.0, 0.012, 0.035]:
-            for slope in [1.0/300]:         #[0.0, 1.0/300, 1.0/150]:
-                # Loop over a range of bedslopes representing
-                # sub to super critical flows 
-                if verbose:
-                    print
-                    print 'Slope:', slope, 'Mannings n:', mannings_n
-                domain = Domain(points, vertices, boundary)   
-                domain.set_name('Inflow_flowline_test')     # Output name
-
-                #--------------------------------------------------------------
-                # Setup initial conditions
-                #--------------------------------------------------------------
-
-                def topography(x, y):
-                    z = -x * slope
-                    return z
-
-                # Use function for elevation
-                domain.set_quantity('elevation', topography)
-                # Constant friction
-                domain.set_quantity('friction', mannings_n)
-                
-                #domain.set_quantity('stage', expression='elevation')
-                     
-                # Set initial flow as depth=1m, uh=1.0 m/s, vh = 0.0
-                # making it 20 m^3/s across entire domain 
-                domain.set_quantity('stage', expression='elevation + %f' % d)
-                domain.set_quantity('xmomentum', uh)
-                domain.set_quantity('ymomentum', vh)                
-
-                #--------------------------------------------------------------
-                # Setup boundary conditions
-                #--------------------------------------------------------------
-
-                Br = Reflective_boundary(domain)      # Solid reflective wall
-                
-                # Constant flow in and out of domain
-                # Depth = 1m, uh=1 m/s, i.e. a flow of 20 m^3/s 
-                # across boundaries
-                Bi = Dirichlet_boundary([d, uh, vh]) 
-                Bo = Dirichlet_boundary([-length*slope+d, uh, vh])
-                #Bo = Dirichlet_boundary([-100, 0, 0])
-
-                domain.set_boundary({'left': Bi, 'right': Bo,
-                                     'top': Br,  'bottom': Br})
-
-                #--------------------------------------------------------------
-                # Evolve system through time
-                #--------------------------------------------------------------
-
-                for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
-                    if verbose :
-                        print domain.timestepping_statistics()
-                        print domain.volumetric_balance_statistics()
-
-                # 90 degree flowline at 200m
-                q = domain.get_flow_through_cross_section([[200.0,  0.0],
-                                                           [200.0, 20.0]])
-                msg = ('Predicted flow was %f, should have been %f'
-                       % (q, ref_flow))
-                if verbose:
-                    print ('90 degree flowline: ANUGA = %f, Ref = %f'
-                           % (q, ref_flow))
-
-                # 45 degree flowline at 200m
-                q = domain.get_flow_through_cross_section([[200.0,  0.0],
-                                                           [220.0, 20.0]])
-                msg = ('Predicted flow was %f, should have been %f'
-                       % (q, ref_flow))
-                if verbose:
-                    print ('45 degree flowline: ANUGA = %f, Ref = %f'
-                           % (q, ref_flow))
-
-                # Stage recorder (gauge) in middle of plane at 200m
-                x = 200.0
-                y = 10.00                
-                w = domain.get_quantity('stage').\
-                        get_values(interpolation_points=[[x, y]])[0]
-                z = domain.get_quantity('elevation').\
-                        get_values(interpolation_points=[[x, y]])[0]
-                domain_depth = w-z
-
-                xmom = domain.get_quantity('xmomentum').\
-                        get_values(interpolation_points=[[x, y]])[0]
-                ymom = domain.get_quantity('ymomentum').\
-                        get_values(interpolation_points=[[x, y]])[0]            
-                if verbose:                    
-                    print ('At interpolation point (h, uh, vh): ',
-                           domain_depth, xmom, ymom)
-                    print 'uh * d * width = ', xmom*domain_depth*width
-                                
-                if slope > 0.0:
-                    # Compute normal depth at gauge location using Manning eqn
-                    # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6
-                    normal_depth = (ref_flow*mannings_n/(slope**0.5*width))**0.6
-                    if verbose:
-                        print ('Depth: ANUGA = %f, Mannings = %f'
-                               % (domain_depth, normal_depth))
-
-        os.remove('Inflow_flowline_test.sww')
-
-
-#################################################################################
+ ################################################################################
 
 if __name__ == "__main__":
