Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py	(revision 6173)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py	(revision 6174)
@@ -185,11 +185,11 @@
             #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
 
-            n0 = num.array([x2 - x1, y2 - y1])
+            n0 = num.array([x2 - x1, y2 - y1], num.Float)
             l0 = num.sqrt(num.sum(n0**2))
 
-            n1 = num.array([x0 - x2, y0 - y2])
+            n1 = num.array([x0 - x2, y0 - y2], num.Float)
             l1 = num.sqrt(num.sum(n1**2))
 
-            n2 = num.array([x1 - x0, y1 - y0])
+            n2 = num.array([x1 - x0, y1 - y0], num.Float)
             l2 = num.sqrt(num.sum(n2**2))
 
@@ -273,5 +273,5 @@
             if not self.geo_reference.is_absolute():
                 return V + num.array([self.geo_reference.get_xllcorner(),
-                                      self.geo_reference.get_yllcorner()])
+                                      self.geo_reference.get_yllcorner()], num.Float)
             else:
                 return V
@@ -315,10 +315,10 @@
             if absolute is True and not self.geo_reference.is_absolute():
                 offset=num.array([self.geo_reference.get_xllcorner(),
-                                  self.geo_reference.get_yllcorner()])
+                                  self.geo_reference.get_yllcorner()], num.Float)
                 return num.array([V[i3,:]+offset,
                                   V[i3+1,:]+offset,
-                                  V[i3+2,:]+offset])
+                                  V[i3+2,:]+offset], num.Float)
             else:
-                return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]])
+                return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.Float)
                 
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py	(revision 6173)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py	(revision 6174)
@@ -200,7 +200,7 @@
             
             # Compute midpoints
-            if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2])
-            if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2])
-            if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2])
+            if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2], num.Float)
+            if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2], num.Float)
+            if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2], num.Float)
 
             # Convert to absolute UTM coordinates
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 6173)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py	(revision 6174)
@@ -124,5 +124,5 @@
 
             #Compute centroid
-            centroid = num.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3])
+            centroid = num.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3], num.Float)
             self.centroid_coordinates[i] = centroid
 
@@ -133,7 +133,7 @@
 
                 #Midpoints
-                m0 = num.array([(x1 + x2)/2, (y1 + y2)/2])
-                m1 = num.array([(x0 + x2)/2, (y0 + y2)/2])
-                m2 = num.array([(x1 + x0)/2, (y1 + y0)/2])
+                m0 = num.array([(x1 + x2)/2, (y1 + y2)/2], num.Float)
+                m1 = num.array([(x0 + x2)/2, (y0 + y2)/2], num.Float)
+                m2 = num.array([(x1 + x0)/2, (y1 + y0)/2], num.Float)
 
                 #The radius is the distance from the centroid of
@@ -1136,6 +1136,6 @@
 
             # Distances from line origin to the two intersections
-            z0 = num.array([x0 - xi0, y0 - eta0])
-            z1 = num.array([x1 - xi0, y1 - eta0])              
+            z0 = num.array([x0 - xi0, y0 - eta0], num.Float)
+            z1 = num.array([x1 - xi0, y1 - eta0], num.Float)
             d0 = num.sqrt(num.sum(z0**2)) 
             d1 = num.sqrt(num.sum(z1**2))
@@ -1152,7 +1152,7 @@
             # Normal direction:
             # Right hand side relative to line direction
-            vector = num.array([x1 - x0, y1 - y0]) # Segment vector
+            vector = num.array([x1 - x0, y1 - y0], num.Float) # Segment vector
             length = num.sqrt(num.sum(vector**2))      # Segment length
-            normal = num.array([vector[1], -vector[0]])/length
+            normal = num.array([vector[1], -vector[0]], num.Float)/length
 
 
@@ -1234,5 +1234,5 @@
         assert isinstance(segment, Triangle_intersection), msg
         
-        midpoint = num.sum(num.array(segment.segment))/2
+        midpoint = num.sum(num.array(segment.segment, num.Float))/2
         midpoints.append(midpoint)
 
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 6173)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 6174)
@@ -1175,5 +1175,5 @@
                     sum += self.vertex_values[triangle_id, vertex_id]
                 vert_values.append(sum/len(triangles))
-            return num.array(vert_values)
+            return num.array(vert_values, num.Float)
         else:
             if (indices is None):
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py	(revision 6173)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py	(revision 6174)
@@ -345,6 +345,6 @@
 
 
-        A = num.array([[1,2,3], [5,5,-5], [0,0,9], [-6,3,3]], 'f')
-        B = num.array([[2,4,4], [3,2,1], [6,-3,4], [4,5,-1]], 'f')
+        A = num.array([[1,2,3], [5,5,-5], [0,0,9], [-6,3,3]], num.Float)
+        B = num.array([[2,4,4], [3,2,1], [6,-3,4], [4,5,-1]], num.Float)
         
         #print A
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py	(revision 6173)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py	(revision 6174)
@@ -59,5 +59,5 @@
         
         #bac, bce, ecf, dbe, daf, dae
-        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.Int)
 
         domain = General_mesh(nodes, triangles,
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py	(revision 6173)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py	(revision 6174)
@@ -22,5 +22,5 @@
 
 def distance(x, y):
-    return sqrt( sum( (num.array(x)-num.array(y))**2 ))
+    return sqrt(num.sum((num.array(x)-num.array(y))**2))
 
 class Test_Mesh(unittest.TestCase):
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py	(revision 6173)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py	(revision 6174)
@@ -1766,5 +1766,5 @@
         quantity = Quantity(self.mesh4)
 
-        quantity.vertex_values = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],num.Float)
+        quantity.vertex_values = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.Float)
 
         quantity.interpolate_from_vertices_to_edges()
@@ -1782,5 +1782,5 @@
                                           [3., 2.5, 1.5],
                                           [3.5, 4.5, 3.],
-                                          [2.5, 3.5, 2]],num.Float)
+                                          [2.5, 3.5, 2]], num.Float)
 
         quantity.interpolate_from_edges_to_vertices()
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 6173)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 6174)
@@ -2097,5 +2097,5 @@
                 #add tide to stage if provided
                 if quantity == 'stage':
-                     quantity_value[quantity] = num.array(quantity_value[quantity]) \
+                     quantity_value[quantity] = num.array(quantity_value[quantity], num.Float) \
                                                           + directory_add_tide
 
Index: /anuga_core/source/anuga/alpha_shape/alpha_shape.py
===================================================================
--- /anuga_core/source/anuga/alpha_shape/alpha_shape.py	(revision 6173)
+++ /anuga_core/source/anuga/alpha_shape/alpha_shape.py	(revision 6174)
@@ -242,10 +242,10 @@
         ind3 = [self.deltri[j][2] for j in range(len(self.deltri))]
 
-        x1 = num.array([ x[j] for j in ind1 ])
-        y1 = num.array([ y[j] for j in ind1 ])
-        x2 = num.array([ x[j] for j in ind2 ])
-        y2 = num.array([ y[j] for j in ind2 ])
-        x3 = num.array([ x[j] for j in ind3 ])
-        y3 = num.array([ y[j] for j in ind3 ])
+        x1 = num.array([x[j] for j in ind1])
+        y1 = num.array([y[j] for j in ind1])
+        x2 = num.array([x[j] for j in ind2])
+        y2 = num.array([y[j] for j in ind2])
+        x3 = num.array([x[j] for j in ind3])
+        y3 = num.array([y[j] for j in ind3])
 
         x21 = x2-x1
Index: /anuga_core/source/anuga/fit_interpolate/search_functions.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/search_functions.py	(revision 6173)
+++ /anuga_core/source/anuga/fit_interpolate/search_functions.py	(revision 6174)
@@ -22,5 +22,7 @@
                         num.array([max_float,max_float]),
                         num.array([max_float,max_float])),
-                       (num.array([1,1]),num.array([0,0]),num.array([-1.1,-1.1]))]]]
+                       (num.array([1,1], num.Int),      #array default#
+                        num.array([0,0], num.Int),      #array default#
+                        num.array([-1.1,-1.1]))]]]
 
 def search_tree_of_vertices(root, mesh, x):
Index: /anuga_core/source/anuga/fit_interpolate/test_fit.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/test_fit.py	(revision 6173)
+++ /anuga_core/source/anuga/fit_interpolate/test_fit.py	(revision 6174)
@@ -25,5 +25,5 @@
 
 def distance(x, y):
-    return sqrt( sum( (num.array(x)-num.array(y))**2 ))
+    return sqrt(num.sum((num.array(x)-num.array(y))**2))
 
 def linear_function(point):
@@ -582,5 +582,5 @@
 
         #Define f(x,y) = x
-        f = num.array([0,0,2]) #Value at global vertex 2
+        f = num.array([0,0,2], num.Int) #Value at global vertex 2      #array default#
 
         #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
@@ -589,5 +589,5 @@
 
         #Define f(x,y) = y
-        f = num.array([0,2,0])  #Value at global vertex 1
+        f = num.array([0,2,0], num.Int)  #Value at global vertex 1      #array default#
 
         #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
@@ -596,5 +596,5 @@
 
         #Define f(x,y) = x+y
-        f = num.array([0,2,2])  #Values at global vertex 1 and 2
+        f = num.array([0,2,2], num.Int)  #Values at global vertex 1 and 2      #array default#
 
         #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
@@ -623,5 +623,5 @@
 
         #Define f(x,y) = x
-        f = num.array([0,0,2,0,2,4]) #f evaluated at points a-f
+        f = num.array([0,0,2,0,2,4], num.Int) #f evaluated at points a-f      #array default#
 
         #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
@@ -630,5 +630,5 @@
 
         #Define f(x,y) = y
-        f = num.array([0,2,0,4,2,0]) #f evaluated at points a-f
+        f = num.array([0,2,0,4,2,0], num.Int) #f evaluated at points a-f      #array default#
 
         #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
@@ -637,5 +637,5 @@
 
         #Define f(x,y) = x+y
-        f = num.array([0,2,2,4,4,4])  #f evaluated at points a-f
+        f = num.array([0,2,2,4,4,4], num.Int)  #f evaluated at points a-f     #array default#
 
         #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
Index: /anuga_core/source/anuga/fit_interpolate/test_interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/test_interpolate.py	(revision 6173)
+++ /anuga_core/source/anuga/fit_interpolate/test_interpolate.py	(revision 6174)
@@ -29,5 +29,5 @@
 
 def distance(x, y):
-    return sqrt( num.sum( (num.array(x)-num.array(y))**2 ))
+    return sqrt(num.sum((num.array(x)-num.array(y))**2))
 
 def linear_function(point):
@@ -760,5 +760,5 @@
 
         interp = Interpolate(vertices, triangles)
-        f = num.array([linear_function(vertices),2*linear_function(vertices) ])
+        f = num.array([linear_function(vertices),2*linear_function(vertices)])
         f = num.transpose(f)
         #print "f",f
@@ -803,5 +803,5 @@
 
         interp = Interpolate(vertices, triangles)
-        f = num.array([linear_function(vertices),2*linear_function(vertices) ])
+        f = num.array([linear_function(vertices),2*linear_function(vertices)])
         f = num.transpose(f)
         #print "f",f
@@ -819,5 +819,5 @@
             
         f = num.array([linear_function(vertices),2*linear_function(vertices),
-                       2*linear_function(vertices) - 100  ])
+                       2*linear_function(vertices) - 100])
         f = num.transpose(f)
         #print "f",f
@@ -829,5 +829,5 @@
             answer = num.array([linear_function(point_coords),
                                 2*linear_function(point_coords) ,
-                                2*linear_function(point_coords)-100 ])
+                                2*linear_function(point_coords)-100])
             z = num.transpose(z)
             #print "z",z
@@ -865,5 +865,5 @@
         
         interp = Interpolate(vertices, triangles)
-        f = num.array([linear_function(vertices),2*linear_function(vertices) ])
+        f = num.array([linear_function(vertices),2*linear_function(vertices)])
         f = num.transpose(f)
         #print "f",f
@@ -883,5 +883,5 @@
             
         f = num.array([linear_function(vertices),2*linear_function(vertices),
-                       2*linear_function(vertices) - 100  ])
+                       2*linear_function(vertices) - 100])
         f = num.transpose(f)
         #print "f",f
@@ -896,5 +896,5 @@
                                                               absolute = True)) ,
                                                               2*linear_function(point_coords.get_data_points( \
-                                                              absolute = True))-100 ])
+                                                              absolute = True))-100])
             z = num.transpose(z)
             #print "z",z
@@ -938,5 +938,5 @@
         
         interp = Interpolate(vertices, triangles)
-        f = num.array([linear_function(vertices),2*linear_function(vertices) ])
+        f = num.array([linear_function(vertices),2*linear_function(vertices)])
         f = num.transpose(f)
         #print "f",f
@@ -984,5 +984,5 @@
 
         interp = Interpolate(vertices, triangles)
-        f = num.array([linear_function(vertices),2*linear_function(vertices) ])
+        f = num.array([linear_function(vertices),2*linear_function(vertices)])
         f = num.transpose(f)
         z = interp.interpolate(f, point_coords,
@@ -1047,5 +1047,5 @@
 
         interp = Interpolate(vertices, triangles)
-        f = num.array([linear_function(vertices),2*linear_function(vertices) ])
+        f = num.array([linear_function(vertices), 2*linear_function(vertices)])
         f = num.transpose(f)
         z = interp.interpolate(f, point_coords)
@@ -1363,5 +1363,5 @@
         
         interp = Interpolate(points, triangles)
-        f = num.array([linear_function(points),2*linear_function(points) ])
+        f = num.array([linear_function(points),2*linear_function(points)])
         f = num.transpose(f)
         #print "f",f
@@ -1653,5 +1653,5 @@
 
         interp = Interpolate(vertices, triangles)
-        f = num.array([linear_function(vertices),2*linear_function(vertices) ])
+        f = num.array([linear_function(vertices),2*linear_function(vertices)])
         f = num.transpose(f)
         #print "f",f
Index: /anuga_core/source/anuga/mesh_engine/test_generate_mesh.py
===================================================================
--- /anuga_core/source/anuga/mesh_engine/test_generate_mesh.py	(revision 6173)
+++ /anuga_core/source/anuga/mesh_engine/test_generate_mesh.py	(revision 6174)
@@ -43,9 +43,9 @@
                               pointattlist,segattlist, mode, points)
 
-        correct = num.array([(1, 0, 2), (2, 3, 1)])
+        correct = num.array([(1, 0, 2), (2, 3, 1)], num.Int)      #array default#
         self.failUnless(num.alltrue(data['generatedtrianglelist'].flat == \
                                     correct.flat),
                         'trianglelist is wrong!')
-        correct = num.array([(0, 1), (1, 3), (3, 2), (2, 0)])
+        correct = num.array([(0, 1), (1, 3), (3, 2), (2, 0)], num.Int)      #array default#
         self.failUnless(num.alltrue(data['generatedsegmentlist'].flat == \
                                     correct.flat),
@@ -74,9 +74,9 @@
         data = generate_mesh(points,seglist,holelist,regionlist,
                               pointattlist,segattlist, mode, points)
-        correct = num.array([(1, 0, 2), (2, 3, 1)])
+        correct = num.array([(1, 0, 2), (2, 3, 1)], num.Int)      #array default#
         self.failUnless(num.alltrue(data['generatedtrianglelist'].flat == \
                                     correct.flat),
                         'trianglelist is wrong!')
-        correct = num.array([(0, 1), (1, 3), (3, 2), (2, 0)])
+        correct = num.array([(0, 1), (1, 3), (3, 2), (2, 0)], num.Int)      #array default#
         self.failUnless(num.alltrue(data['generatedsegmentlist'].flat == \
                                     correct.flat),
@@ -154,9 +154,9 @@
                               pointattlist,segattlist, mode, points)
 
-        correct = num.array([(1, 0, 2), (2, 3, 1)])
+        correct = num.array([(1, 0, 2), (2, 3, 1)], num.Int)      #array default#
         self.failUnless(num.alltrue(data['generatedtrianglelist'].flat == \
                                     correct.flat),
                         'trianglelist is wrong!')
-        correct = num.array([(0, 1), (1, 3), (3, 2), (2, 0)])
+        correct = num.array([(0, 1), (1, 3), (3, 2), (2, 0)], num.Int)      #array default#
         self.failUnless(num.alltrue(data['generatedsegmentlist'].flat == \
                                     correct.flat),
@@ -170,5 +170,5 @@
         
         self.failUnless(num.alltrue(data['generatedsegmentmarkerlist'] == \
-                                    num.array([1,2,3,4])),
+                                    num.array([1,2,3,4], num.Int)),      #array default#
                         'Failed!')
         
@@ -388,9 +388,9 @@
 
      
-        correct = num.array([(1, 0, 2), (2, 3, 1)])
+        correct = num.array([(1, 0, 2), (2, 3, 1)], num.Int)      #array default#
         self.failUnless(num.alltrue(data['generatedtrianglelist'].flat == \
                                     correct.flat),
                         'trianglelist is wrong!')
-        correct = num.array([(0, 1), (1, 3), (3, 2), (2, 0)])
+        correct = num.array([(0, 1), (1, 3), (3, 2), (2, 0)], num.Int)      #array default#
         self.failUnless(num.alltrue(data['generatedsegmentlist'].flat == \
                                     correct.flat),
@@ -407,10 +407,10 @@
                                      correct.flat),
                         'Failed')
-        correct = num.array([(0, 1), (1, 3), (3, 2), (2, 0)])
+        correct = num.array([(0, 1), (1, 3), (3, 2), (2, 0)], num.Int)      #array default#
         self.failUnless(num.alltrue(data['generatedsegmentlist'].flat == \
                                     correct.flat),
                         'Failed!')
         
-        correct = num.array(segattlist)
+        correct = num.array(segattlist, num.Int)      #array default#
         self.failUnless(num.allclose(data['generatedsegmentmarkerlist'].flat, 
                                      correct.flat),
@@ -418,5 +418,5 @@
         
         # I copied these answers from the output, so bad test..
-        correct = num.array([(-1, 1, -1), (-1, 0, -1)])
+        correct = num.array([(-1, 1, -1), (-1, 0, -1)], num.Int)      #array default#
         self.failUnless(num.alltrue(data['generatedtriangleneighborlist'].flat == \
                                     correct.flat),
Index: /anuga_core/source/anuga/shallow_water/test_data_manager.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 6173)
+++ /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 6174)
@@ -4950,5 +4950,5 @@
 
         ## 7th test
-        m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
+        m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]], num.Int)      #array default#
         kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
             latitudes,longitudes,
@@ -5010,5 +5010,5 @@
         longitudes = [148,149,150,151]
 
-        m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
+        m2d = num.array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]], num.Int)      #array default#
 
         # k - lat
@@ -10626,5 +10626,5 @@
 
         # Check runup restricted to a polygon
-        p = num.array([[50,1], [99,1], [99,49], [50,49]]) + num.array([E, N])
+        p = num.array([[50,1], [99,1], [99,49], [50,49]], num.Int) + num.array([E, N], num.Int)      #array default#
 
         runup = get_maximum_inundation_elevation(swwfile, polygon=p)
Index: /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 6173)
+++ /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 6174)
@@ -2,5 +2,5 @@
 
 import unittest, os
-from math import pi
+from math import pi, sqrt
 import tempfile
 
@@ -254,5 +254,5 @@
         #Check error check
         try:
-            rotate(r, num.array([1,1,1]) )
+            rotate(r, num.array([1,1,1], num.Int))      #array default#
         except:
             pass
@@ -338,5 +338,5 @@
     def test_flux3(self):
         #Use data from previous version of abstract_2d_finite_volumes
-        normal = num.array([-num.sqrt(2)/2, num.sqrt(2)/2])
+        normal = num.array([-sqrt(2)/2, sqrt(2)/2])
         ql = num.array([-0.075, 2, 3])
         qr = num.array([-0.075, 2, 3])
@@ -347,10 +347,10 @@
         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)        
 
-        assert num.allclose(edgeflux, [num.sqrt(2)/2, 4.40221112, 7.3829019])
+        assert num.allclose(edgeflux, [sqrt(2)/2, 4.40221112, 7.3829019])
         assert num.allclose(max_speed, 4.0716654239)
 
     def test_flux4(self):
         #Use data from previous version of abstract_2d_finite_volumes
-        normal = num.array([-num.sqrt(2)/2, num.sqrt(2)/2])
+        normal = num.array([-sqrt(2)/2, sqrt(2)/2])
         ql = num.array([-0.34319278, 0.10254161, 0.07273855])
         qr = num.array([-0.30683287, 0.1071986, 0.05930515])
Index: /anuga_core/source/anuga/utilities/numerical_tools.py
===================================================================
--- /anuga_core/source/anuga/utilities/numerical_tools.py	(revision 6173)
+++ /anuga_core/source/anuga/utilities/numerical_tools.py	(revision 6174)
@@ -254,5 +254,6 @@
         if type(A) == num.ArrayType:
             if A.typecode == typecode:
-                return num.array(A)  #FIXME: Shouldn't this just return A?
+#                return num.array(A)  #FIXME: Shouldn't this just return A?
+                return A
             else:
                 return num.array(A,typecode)
Index: /anuga_core/source/anuga/utilities/test_numerical_tools.py
===================================================================
--- /anuga_core/source/anuga/utilities/test_numerical_tools.py	(revision 6173)
+++ /anuga_core/source/anuga/utilities/test_numerical_tools.py	(revision 6174)
@@ -90,5 +90,5 @@
         assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
 
-        A = num.array([1,2,3,4])
+        A = num.array([1,2,3,4], num.Int)      #array default#
         B = ensure_numeric(A)
         assert type(B) == num.ArrayType
@@ -97,5 +97,5 @@
         assert A is B   #Same object
 
-        A = num.array([1,2,3,4])
+        A = num.array([1,2,3,4], num.Int)      #array default#
         B = ensure_numeric(A, num.Float)
         assert type(B) == num.ArrayType
Index: /anuga_core/source/anuga/utilities/test_sparse.py
===================================================================
--- /anuga_core/source/anuga/utilities/test_sparse.py	(revision 6173)
+++ /anuga_core/source/anuga/utilities/test_sparse.py	(revision 6174)
@@ -86,5 +86,5 @@
 
         #Right hand side column
-        v = num.array([[2,4],[3,4],[4,4]])
+        v = num.array([[2,4],[3,4],[4,4]], num.Int)      #array default#
 
         u = A*v[:,0]
@@ -104,5 +104,5 @@
 
         #Right hand side matrix
-        v = num.array([[2,4],[3,4],[4,4]])
+        v = num.array([[2,4],[3,4],[4,4]], num.Int)      #array default#
 
         u = A*v
