Index: /trunk/anuga_core/source/anuga/file/test_csv.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/test_csv.py	(revision 7869)
+++ /trunk/anuga_core/source/anuga/file/test_csv.py	(revision 7870)
@@ -3,4 +3,6 @@
 import tempfile
 import numpy as num
+
+from anuga.utilities.system_tools import get_pathname_from_package
 
 from csv_file import load_csv_as_array, load_csv_as_dict, store_parameters, \
Index: /trunk/anuga_core/source/anuga/file/test_sww.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/test_sww.py	(revision 7869)
+++ /trunk/anuga_core/source/anuga/file/test_sww.py	(revision 7870)
@@ -9,4 +9,5 @@
 from anuga.shallow_water.shallow_water_domain import Domain
 from sww import load_sww_as_domain, weed, get_mesh_and_quantities_from_file
+from Scientific.IO.NetCDF import NetCDFFile
 
 # boundary functions
Index: /trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 7869)
+++ /trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 7870)
@@ -65,5 +65,6 @@
     TRAC administration of ANUGA (User Manuals etc) at
     https://datamining.anu.edu.au/anuga and Subversion repository at
-    $HeadURL$
+    $HeadURL: https://datamining.anu.edu.au/svn/anuga/trunk/anuga_core/source/
+                anuga/shallow_water/shallow_water_domain.py $
 
 Constraints: See GPL license in the user guide
@@ -74,10 +75,4 @@
 """
 
-# Subversion keywords:
-#
-# $LastChangedDate$
-# $LastChangedRevision$
-# $LastChangedBy$
-
 
 import numpy as num
@@ -87,51 +82,13 @@
 
 from anuga.shallow_water.forcing import Cross_section
-from anuga.pmesh.mesh_interface import create_mesh_from_regions
-from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
-from anuga.geospatial_data.geospatial_data import ensure_geospatial
+from anuga.utilities.numerical_tools import mean
 from anuga.file.sww import SWW_file  
-from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
-
-from anuga.fit_interpolate.interpolate import Modeltime_too_late, \
-                                              Modeltime_too_early
-
-from anuga.geometry.polygon import inside_polygon, polygon_area, \
-                                    is_inside_polygon
+            
 import anuga.utilities.log as log
 
 import types
-from types import IntType, FloatType
-from warnings import warn
-
-
-################################################################################
-# Shallow water domain
-################################################################################
-
-##
-# @brief Class for a shallow water domain.
+
 class Domain(Generic_Domain):
-
-    #conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
-    #other_quantities = ['elevation', 'friction']
-
-    ##
-    # @brief Instantiate a shallow water domain.
-    # @param coordinates
-    # @param vertices
-    # @param boundary
-    # @param tagged_elements
-    # @param geo_reference
-    # @param use_inscribed_circle
-    # @param mesh_filename
-    # @param use_cache
-    # @param verbose
-    # @param evolved_quantities
-    # @param full_send_dict
-    # @param ghost_recv_dict
-    # @param processor
-    # @param numproc
-    # @param number_of_full_nodes
-    # @param number_of_full_triangles
+    """ Class for a shallow water domain."""
     def __init__(self,
                  coordinates=None,
@@ -153,4 +110,23 @@
                  number_of_full_nodes=None,
                  number_of_full_triangles=None):
+        """
+            Instantiate a shallow water domain.
+            coordinates - vertex locations for the mesh
+            vertices - vertex indices for the mesh
+            boundary - boundaries of the mesh
+            # @param tagged_elements
+            # @param geo_reference
+            # @param use_inscribed_circle
+            # @param mesh_filename
+            # @param use_cache
+            # @param verbose
+            # @param evolved_quantities
+            # @param full_send_dict
+            # @param ghost_recv_dict
+            # @param processor
+            # @param numproc
+            # @param number_of_full_nodes
+            # @param number_of_full_triangles
+        """
 
         # Define quantities for the shallow_water domain
@@ -165,22 +141,22 @@
         
         Generic_Domain.__init__(self,
-                                coordinates,
-                                vertices,
-                                boundary,
-                                conserved_quantities,
-                                evolved_quantities,
-                                other_quantities,
-                                tagged_elements,
-                                geo_reference,
-                                use_inscribed_circle,
-                                mesh_filename,
-                                use_cache,
-                                verbose,
-                                full_send_dict,
-                                ghost_recv_dict,
-                                processor,
-                                numproc,
-                                number_of_full_nodes=number_of_full_nodes,
-                                number_of_full_triangles=number_of_full_triangles)
+                            coordinates,
+                            vertices,
+                            boundary,
+                            conserved_quantities,
+                            evolved_quantities,
+                            other_quantities,
+                            tagged_elements,
+                            geo_reference,
+                            use_inscribed_circle,
+                            mesh_filename,
+                            use_cache,
+                            verbose,
+                            full_send_dict,
+                            ghost_recv_dict,
+                            processor,
+                            numproc,
+                            number_of_full_nodes=number_of_full_nodes,
+                            number_of_full_triangles=number_of_full_triangles)
 
         self.set_defaults()
@@ -200,9 +176,4 @@
 
 
-
-
-    ##
-    # @brief Set default values, usually read in from a config file
-    # @param flag
     def set_defaults(self):
         """Set the default values in this routine. That way we can inherit class
@@ -212,5 +183,5 @@
         from anuga.config import minimum_storable_height
         from anuga.config import minimum_allowed_height, maximum_allowed_speed
-        from anuga.config import g, epsilon, beta_w, beta_w_dry,\
+        from anuga.config import g, beta_w, beta_w_dry, \
              beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
         from anuga.config import alpha_balance
@@ -219,6 +190,4 @@
         from anuga.config import use_edge_limiter
         from anuga.config import use_centroid_velocities
-
-
 
         self.set_minimum_allowed_height(minimum_allowed_height)
@@ -247,7 +216,5 @@
         self.use_centroid_velocities = use_centroid_velocities       
 
-    ##
-    # @brief
-    # @param flag
+
     def set_new_mannings_function(self, flag=True):
         """Cludge to allow unit test to pass, but to
@@ -263,7 +230,4 @@
 
 
-    ##
-    # @brief
-    # @param flag
     def set_use_edge_limiter(self, flag=True):
         """Cludge to allow unit test to pass, but to
@@ -277,8 +241,4 @@
 
 
-          
-    ##
-    # @brief
-    # @param beta
     def set_all_limiters(self, beta):
         """Shorthand to assign one constant value [0,1] to all limiters.
@@ -299,8 +259,5 @@
         self.quantities['ymomentum'].beta = beta
 
-    ##
-    # @brief
-    # @param flag
-    # @param reduction
+
     def set_store_vertices_uniquely(self, flag, reduction=None):
         """Decide whether vertex values should be stored uniquely as
@@ -406,5 +363,6 @@
         this function will have no effect.
         
-        The format, where q is a list of names is for backwards compatibility only.
+        The format, where q is a list of names is for backwards compatibility
+        only.
         It will take the specified quantities to be time dependent and assume 
         'elevation' to be static regardless.
@@ -416,5 +374,5 @@
             return
 
-        # Check correcness
+        # Check correctness
         for quantity_name in q:
             msg = ('Quantity %s is not a valid conserved quantity'
@@ -422,21 +380,4 @@
             assert quantity_name in self.quantities, msg
 
-        if type(q) == types.ListType:
-
-            msg = 'List arguments to set_quantities_to_be_stored '
-            msg += 'has been deprecated and will be removed in future '
-            msg += 'versions of ANUGA.'
-            msg += 'Please use dictionary argument instead'
-            warn(msg, DeprecationWarning) 
-
-        
-        
-            # FIXME: Raise deprecation warning
-            tmp = {}
-            for x in q:
-                tmp[x] = 2
-            tmp['elevation'] = 1    
-            q = tmp     
-            
         assert type(q) == types.DictType    
         self.quantities_to_be_stored = q
@@ -567,7 +508,6 @@
 
 
-    ##
-    # @brief Run integrity checks on shallow water domain.
     def check_integrity(self):
+        """ Run integrity checks on shallow water domain. """
         Generic_Domain.check_integrity(self)
 
@@ -580,20 +520,16 @@
         assert self.conserved_quantities[2] == 'ymomentum', msg
 
-    ##
-    # @brief 
     def extrapolate_second_order_sw(self):
-        #Call correct module function (either from this module or C-extension)
+        """Call correct module function
+            (either from this module or C-extension)"""
         extrapolate_second_order_sw(self)
 
-    ##
-    # @brief 
     def compute_fluxes(self):
-        #Call correct module function (either from this module or C-extension)
+        """Call correct module function
+            (either from this module or C-extension)"""
         compute_fluxes(self)
 
-    ##
-    # @brief 
     def distribute_to_vertices_and_edges(self):
-        # Call correct module function
+        """ Call correct module function """
         if self.use_edge_limiter:
             distribute_using_edge_limiter(self)
@@ -603,10 +539,4 @@
 
 
-    ##
-    # @brief Evolve the model by one step.
-    # @param yieldstep 
-    # @param finaltime 
-    # @param duration 
-    # @param skip_initial_step 
     def evolve(self,
                yieldstep=None,
@@ -614,5 +544,8 @@
                duration=None,
                skip_initial_step=False):
-        """Specialisation of basic evolve method from parent class"""
+        """Specialisation of basic evolve method from parent class.
+        
+            Evolve the model by 1 step.
+        """
 
         # Call check integrity here rather than from user scripts
@@ -642,6 +575,5 @@
             yield(t)
 
-    ##
-    # @brief 
+
     def initialise_storage(self):
         """Create and initialise self.writer object for storing data.
@@ -655,7 +587,5 @@
         self.writer.store_connectivity()
 
-    ##
-    # @brief 
-    # @param name 
+
     def store_timestep(self):
         """Store time dependent quantities and time.
@@ -667,8 +597,5 @@
         self.writer.store_timestep()
 
-    ##
-    # @brief Get time stepping statistics string for printing.
-    # @param track_speeds 
-    # @param triangle_id 
+
     def timestepping_statistics(self,
                                 track_speeds=False,
@@ -712,14 +639,14 @@
             Ch = Cw-Cz
 
-            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
-                 %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
-
-            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
-                 %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
-
-            s += '    %s: centroid_value = %.4f\n'\
-                 %(name.ljust(qwidth), Ch[0])
-
-            msg += s
+            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
+                 % (name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
+
+            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
+                 % (name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
+
+            message += '    %s: centroid_value = %.4f\n'\
+                 % (name.ljust(qwidth), Ch[0])
+
+            msg += message
 
             uh = self.quantities['xmomentum']
@@ -739,14 +666,14 @@
             Cu = Cuh/(Ch + epsilon)
             name = 'U'
-            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
-                 %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
-
-            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
-                 %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
-
-            s += '    %s: centroid_value = %.4f\n'\
-                 %(name.ljust(qwidth), Cu[0])
-
-            msg += s
+            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n' \
+                 % (name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
+
+            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n' \
+                 % (name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
+
+            message += '    %s: centroid_value = %.4f\n' \
+                 % (name.ljust(qwidth), Cu[0])
+
+            msg += message
 
             Vv = Vvh/(Vh + epsilon)
@@ -754,14 +681,14 @@
             Cv = Cvh/(Ch + epsilon)
             name = 'V'
-            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
-                 %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
-
-            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
-                 %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
-
-            s += '    %s: centroid_value = %.4f\n'\
+            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n' \
+                 % (name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
+
+            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n' \
+                 % (name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
+
+            message += '    %s: centroid_value = %.4f\n'\
                  %(name.ljust(qwidth), Cv[0])
 
-            msg += s
+            msg += message
 
             # Froude number in each direction
@@ -771,14 +698,14 @@
             Cfx = Cu/(num.sqrt(g*Ch) + epsilon)
 
-            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
-                 %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
-
-            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
-                 %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
-
-            s += '    %s: centroid_value = %.4f\n'\
-                 %(name.ljust(qwidth), Cfx[0])
-
-            msg += s
+            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
+                 % (name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
+
+            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
+                 % (name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
+
+            message += '    %s: centroid_value = %.4f\n'\
+                 % (name.ljust(qwidth), Cfx[0])
+
+            msg += message
 
             name = 'Froude (y)'
@@ -787,14 +714,14 @@
             Cfy = Cv/(num.sqrt(g*Ch) + epsilon)
 
-            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
-                 %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
-
-            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
-                 %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
-
-            s += '    %s: centroid_value = %.4f\n'\
-                 %(name.ljust(qwidth), Cfy[0])
-
-            msg += s
+            message  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
+                 % (name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
+
+            message += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
+                 % (name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
+
+            message += '    %s: centroid_value = %.4f\n'\
+                 % (name.ljust(qwidth), Cfy[0])
+
+            msg += message
 
         return msg
@@ -813,5 +740,6 @@
         # Run through boundary array and compute for each segment
         # the normal momentum ((uh, vh) dot normal) times segment length.
-        # Based on sign accumulate this into boundary_inflow and boundary_outflow.
+        # Based on sign accumulate this into boundary_inflow and
+        # boundary_outflow.
 			
         # Compute flows along boundary
@@ -842,5 +770,5 @@
             if edge_flow > 0:
                 # Flow is inflow      
-                total_boundary_inflow += edge_flow                                  
+                total_boundary_inflow += edge_flow       
             else:
                 # Flow is outflow
@@ -881,8 +809,8 @@
         
         area = self.mesh.get_areas()
-        volume = 0.0
         
         stage = self.get_quantity('stage').get_values(location='centroids')
-        elevation = self.get_quantity('elevation').get_values(location='centroids')        
+        elevation = \
+            self.get_quantity('elevation').get_values(location='centroids')
         depth = stage-elevation
         
@@ -897,23 +825,26 @@
          total_boundary_outflow) = self.compute_boundary_flows() 
         
-        s = '---------------------------\n'        
-        s += 'Volumetric balance report:\n'
-        s += '--------------------------\n'
-        s += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
-        s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow        
-        s += 'Net boundary flow by tags [m^3/s]\n'
+        message = '---------------------------\n'        
+        message += 'Volumetric balance report:\n'
+        message += '--------------------------\n'
+        message += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
+        message += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow
+        message += 'Net boundary flow by tags [m^3/s]\n'
         for tag in boundary_flows:
-            s += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
-        
-        s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow) 
-        s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume()
-        
-        # The go through explicit forcing update and record the rate of change for stage and 
-        # record into forcing_inflow and forcing_outflow. Finally compute integral 
-        # of depth to obtain total volume of domain.
+            message += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
+        
+        message += 'Total net boundary flow [m^3/s]: %.2f\n' % \
+                    (total_boundary_inflow + total_boundary_outflow) 
+        message += 'Total volume in domain [m^3]: %.2f\n' % \
+                    self.compute_total_volume()
+        
+        # The go through explicit forcing update and record the rate of change
+        # for stage and 
+        # record into forcing_inflow and forcing_outflow. Finally compute  
+        # integral of depth to obtain total volume of domain.
 	
         # FIXME(Ole): This part is not yet done.		
         
-        return s        
+        return message        
            
 ################################################################################
@@ -951,6 +882,4 @@
     from shallow_water_ext import compute_fluxes_ext_central \
                                   as compute_fluxes_ext
-
-    N = len(domain)    # number_of_triangles
 
     # Shortcuts
@@ -1001,8 +930,5 @@
     """Wrapper calling C version of extrapolate_second_order_sw"""
 
-    import sys
     from shallow_water_ext import extrapolate_second_order_sw as extrapol2
-
-    N = len(domain) # number_of_triangles
 
     # Shortcuts
@@ -1066,5 +992,5 @@
             domain.extrapolate_second_order_sw()
         else:
-            raise 'Unknown order'
+            raise Exception('Unknown order')
     else:
         # Old code:
@@ -1077,5 +1003,5 @@
                 Q.extrapolate_second_order_and_limit_by_vertex()
             else:
-                raise 'Unknown order'
+                raise Exception('Unknown order')
 
     # Take bed elevation into account when water heights are small
@@ -1119,5 +1045,5 @@
             Q.extrapolate_second_order_and_limit_by_edge()
         else:
-            raise 'Unknown order'
+            raise Exception('Unknown order')
 
     balance_deep_and_shallow(domain)
@@ -1206,11 +1132,10 @@
     elevation = domain.quantities['elevation']
 
-    h = stage.centroid_values - elevation.centroid_values
-    z = elevation.vertex_values
-
-    x = domain.get_vertex_coordinates()
-    g = domain.g
-
-    gravity_c(g, h, z, x, xmom_update, ymom_update)    #, 1.0e-6)
+    height = stage.centroid_values - elevation.centroid_values
+    elevation = elevation.vertex_values
+
+    point = domain.get_vertex_coordinates()
+
+    gravity_c(domain.g, height, elevation, point, xmom_update, ymom_update)
 
 ##
@@ -1241,12 +1166,13 @@
     ymom_update = ymom.semi_implicit_update
 
-    N = len(domain)
     eps = domain.minimum_allowed_height
     g = domain.g
 
     if domain.use_new_mannings:
-        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
+        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, \
+                                ymom_update)
     else:
-        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
+        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, \
+                                ymom_update)
     
 
@@ -1278,17 +1204,16 @@
     ymom_update = ymom.explicit_update
 
-    N = len(domain)
     eps = domain.minimum_allowed_height
-    g = domain.g
-
 
     if domain.use_new_mannings:
-        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
+        manning_friction_new(domain.g, eps, x, w, uh, vh, z, eta, xmom_update, \
+                            ymom_update)
     else:
-        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
-
-
-
-# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
+        manning_friction_old(domain.g, eps, w, uh, vh, z, eta, xmom_update, \
+                            ymom_update)
+
+
+
+# FIXME (Ole): This was implemented for use with one of the analytical solutions
 ##
 # @brief Apply linear friction to a surface.
@@ -1301,6 +1226,4 @@
     """
 
-    from math import sqrt
-
     w = domain.quantities['stage'].centroid_values
     z = domain.quantities['elevation'].centroid_values
@@ -1314,9 +1237,8 @@
     ymom_update = domain.quantities['ymomentum'].semi_implicit_update
 
-    N = len(domain) # number_of_triangles
+    num_tris = len(domain)
     eps = domain.minimum_allowed_height
-    g = domain.g #Not necessary? Why was this added?
-
-    for k in range(N):
+
+    for k in range(num_tris):
         if tau[k] >= eps:
             if h[k] >= eps:
@@ -1330,14 +1252,16 @@
                              surface_roughness_data,
                              verbose=False):
-    """Returns an array of friction values for each wet element adjusted for depth.
+    """Returns an array of friction values for each wet element adjusted for
+            depth.
 
     Inputs:
         domain - computational domain object
         default_friction - depth independent bottom friction
-        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values for each
-        friction region.
+        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values 
+        for each friction region.
 
     Outputs:
-        wet_friction - Array that can be used directly to update friction as follows:
+        wet_friction - Array that can be used directly to update friction as
+                        follows:
                        domain.set_quantity('friction', wet_friction)
 
@@ -1345,54 +1269,59 @@
         
     """
-
-    import numpy as num
     
-    # Create a temp array to store updated depth dependent friction for wet elements
-    # EHR this is outwardly inneficient but not obvious how to avoid recreating each call??????
-    N=len(domain)
-    wet_friction    = num.zeros(N, num.float)
-    wet_friction[:] = default_n0   # Initially assign default_n0 to all array so sure have no zeros values
+    default_n0 = 0  # James - this was missing, don't know what it should be
     
+    # Create a temp array to store updated depth dependent
+    # friction for wet elements
+    # EHR this is outwardly inneficient but not obvious how to avoid
+    # recreating each call??????
+
+    wet_friction    = num.zeros(len(domain), num.float)
+    wet_friction[:] = default_n0  # Initially assign default_n0 to all array so
+                                  # sure have no zeros values
     
-    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
+    # create depth instance for this timestep    
+    depth = domain.create_quantity_from_expression('stage - elevation')  
     # Recompute depth as vector  
-    d = depth.get_values(location='centroids')
+    d_vals = depth.get_values(location='centroids')
  
     # rebuild the 'friction' values adjusted for depth at this instant
-    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
-        
+    # loop for each wet element in domain
+   
+    for i in domain.get_wet_elements():        
         # Get roughness data for each element
-        n0 = float(surface_roughness_data[i,0])
-        d1 = float(surface_roughness_data[i,1])
-        n1 = float(surface_roughness_data[i,2])
-        d2 = float(surface_roughness_data[i,3])
-        n2 = float(surface_roughness_data[i,4])
+        d1 = float(surface_roughness_data[i, 1])
+        n1 = float(surface_roughness_data[i, 2])
+        d2 = float(surface_roughness_data[i, 3])
+        n2 = float(surface_roughness_data[i, 4])
         
         
         # Recompute friction values from depth for this element 
                
-        if d[i]   <= d1: 
-            depth_dependent_friction = n1
-        elif d[i] >= d2:
-            depth_dependent_friction = n2
+        if d_vals[i]   <= d1: 
+            ddf = n1
+        elif d_vals[i] >= d2:
+            ddf = n2
         else:
-            depth_dependent_friction = n1+((n2-n1)/(d2-d1))*(d[i]-d1)
+            ddf = n1 + ((n2-n1)/(d2-d1))*(d_vals[i]-d1)
             
         # check sanity of result
-        if (depth_dependent_friction  < 0.010 or depth_dependent_friction > 9999.0) :
-            log.critical('%s >>>> WARNING: computed depth_dependent friction '
+        if (ddf  < 0.010 or \
+                            ddf > 9999.0) :
+            log.critical('>>>> WARNING: computed depth_dependent friction '
                          'out of range, ddf%f, n1=%f, n2=%f'
-                         % (model_data.basename,
-                            depth_dependent_friction, n1, n2))
+                         % (ddf, n1, n2))
         
         # update depth dependent friction  for that wet element
-        wet_friction[i] = depth_dependent_friction
-        
-    # EHR add code to show range of 'friction across domain at this instant as sanity check?????????
+        wet_friction[i] = ddf
+        
+    # EHR add code to show range of 'friction across domain at this instant as
+    # sanity check?????????
     
     if verbose :
-        nvals=domain.get_quantity('friction').get_values(location='centroids')        # return array of domain nvals
-        n_min=min(nvals)
-        n_max=max(nvals)
+        # return array of domain nvals
+        nvals = domain.get_quantity('friction').get_values(location='centroids')
+        n_min = min(nvals)
+        n_max = max(nvals)
         
         log.critical('         ++++ calculate_depth_dependent_friction - '
@@ -1408,32 +1337,14 @@
 ################################################################################
 
-from anuga.utilities import compile
-if compile.can_use_C_extension('shallow_water_ext.c'):
-    # Underlying C implementations can be accessed
-    from shallow_water_ext import assign_windfield_values
-else:
+def _raise_compile_exception():
+    """ Raise exception if compiler not available. """
     msg = 'C implementations could not be accessed by %s.\n ' % __file__
     msg += 'Make sure compile_all.py has been run as described in '
     msg += 'the ANUGA installation guide.'
-    raise Exception, msg
-
-# Optimisation with psyco
-#from anuga.config import use_psyco
-#if use_psyco:
-    #try:
-        #import psyco
-    #except:
-        #import os
-        #if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
-            #pass
-            ##Psyco isn't supported on 64 bit systems, but it doesn't matter
-        #else:
-            #msg = ('WARNING: psyco (speedup) could not be imported, '
-                   #'you may want to consider installing it')
-            #log.critical(msg)
-    #else:
-        #psyco.bind(Domain.distribute_to_vertices_and_edges)
-        #psyco.bind(Domain.compute_fluxes)
-
+    raise Exception(msg)
+
+from anuga.utilities import compile
+if not compile.can_use_C_extension('shallow_water_ext.c'):
+    _raise_compile_exception()
 
 if __name__ == "__main__":
