Changeset 2648

Timestamp:

Mar 31, 2006, 6:21:26 PM (19 years ago)

Author:

steve

Message:

Looking at island.py example. I changed the order of the implicit and explicit update (in quantity_ext.c) which I think improved the situation a little. But this has left a few fails in the test_all suite which we need to fix up. Mainly small differences in results.

Files:

• development/analytical solutions/Analytical_solution_Sampson_cross_mesh.py (added)
• development/analytical solutions/Analytical_solution_Yoon_cross_mesh.py (added)
• inundation/examples/island.py (modified) (2 diffs)
• inundation/pyvolution/config.py (modified) (1 diff)
• inundation/pyvolution/data_manager.py (modified) (17 diffs)
• inundation/pyvolution/domain.py (modified) (10 diffs)
• inundation/pyvolution/quantity.py (modified) (16 diffs)
• inundation/pyvolution/quantity_ext.c (modified) (3 diffs)
• inundation/pyvolution/shallow_water.py (modified) (11 diffs)
• inundation/pyvolution/shallow_water_ext.c (modified) (6 diffs)
• inundation/pyvolution/test_data_manager.py (modified) (126 diffs)
• inundation/zeus/Validation.zpi (modified) (1 diff)
• inundation/zeus/analytic_solutions.zpi (modified) (1 diff)
• inundation/zeus/anuga-workspace.ini (modified) (1 diff)
• inundation/zeus/anuga-workspace.zwi (modified) (1 diff)
• inundation/zeus/parallel.zpi (modified) (1 diff)
• production/merimbula_2005/Merimbula_17156_refined.tsh (added)
• production/merimbula_2005/Site00tide (added)
• production/merimbula_2005/Site02tide (added)
• production/merimbula_2005/Site03tide (added)
• production/merimbula_2005/Site04tide.eps (added)
• production/merimbula_2005/Site05tide (added)
• production/merimbula_2005/Site06tide (added)
• production/merimbula_2005/Site07tide (added)
• production/merimbula_2005/Site08tide (added)
• production/merimbula_2005/project.py (modified) (1 diff)
• production/merimbula_2005/run_new_meribula_weed.py (modified) (1 diff)

inundation/examples/island.py

@@ -33 +33 @@
                                            'top': [2],
                                            'left': [3]},
-                          maximum_triangle_area = 20,
+                          maximum_triangle_area = 5,
                           filename = 'island.msh'
-                          , interior_regions=[ ([[50,25], [70,25], [70,75], [50,75]], 3)]
+              , interior_regions=[ ([[50,25], [70,25], [70,75], [50,75]], 3)]
                           )
 
@@ -77 +77 @@
     return z
 
-domain.set_quantity('friction', 0.1)  #Honky dory
-#domain.set_quantity('friction', 2)     #Creep
+#domain.set_quantity('friction', 0.1)  #Honky dory
+domain.set_quantity('friction', 1000)     #Creep
 domain.set_quantity('elevation', island)
 domain.set_quantity('stage', 1)
+domain.max_timestep = 0.01
 
 

inundation/pyvolution/config.py

@@ -12 +12 @@
 
 min_timestep = 1.0e-6 #Should be computed based on geometry
-max_timestep = 1000
+max_timestep = 1.0e+3
 #This is how:
 #Define maximal possible speed in open water v_max, e.g. 500m/s (soundspeed?)

inundation/pyvolution/data_manager.py

@@ -1049 +1049 @@
     if zone is None:
         assert xllcorner is None, 'xllcorner must be None'
-        assert yllcorner is None, 'yllcorner must be None'        
+        assert yllcorner is None, 'yllcorner must be None'
         Geo_reference().write_NetCDF(outfile)
     else:
-        Geo_reference(zone, xllcorner, yllcorner).write_NetCDF(outfile)        
-        
+        Geo_reference(zone, xllcorner, yllcorner).write_NetCDF(outfile)
+
 
 
@@ -1075 +1075 @@
 
 
-def dem2pts(basename_in, basename_out=None, 
+def dem2pts(basename_in, basename_out=None,
             easting_min=None, easting_max=None,
             northing_min=None, northing_max=None,
@@ -1099 +1099 @@
 
 
-    kwargs = {'basename_out': basename_out,  
+    kwargs = {'basename_out': basename_out,
               'easting_min': easting_min,
               'easting_max': easting_max,
@@ -1114 +1114 @@
     else:
         result = apply(_dem2pts, [basename_in], kwargs)
-        
+
     return result
 
@@ -1123 +1123 @@
     """Read Digitial Elevation model from the following NetCDF format (.dem)
 
-    Internal function. See public function dem2pts for details.    
+    Internal function. See public function dem2pts for details.
     """
 
@@ -1219 +1219 @@
 
                 if k == 0:
-                    i1_0 = i 
+                    i1_0 = i
                     j1_0 = j
-                k += 1   
+                k += 1
 
     index1 = j1_0
     index2 = thisj
-    
+
     # dimension definitions
     nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
     ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
-    
-    clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0) 
+
+    clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0)
     nopoints = clippednopoints-nn
 
     clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1]
-    
+
     if verbose and nn > 0:
         print 'There are %d values in the elevation' %totalnopoints
         print 'There are %d values in the clipped elevation' %clippednopoints
         print 'There are %d NODATA_values in the clipped elevation' %nn
-        
+
     outfile.createDimension('number_of_points', nopoints)
     outfile.createDimension('number_of_dimensions', 2) #This is 2d data
@@ -1273 +1273 @@
 
         local_index = 0
-        
+
         y = (nrows-i-1)*cellsize + yllcorner
         #for j in range(ncols):
@@ -1286 +1286 @@
                 global_index += 1
                 local_index += 1
-                
+
         upper_index = global_index
 
-        if upper_index == lower_index + newcols: 
+        if upper_index == lower_index + newcols:
             points[lower_index:upper_index, :] = tpoints
             elevation[lower_index:upper_index] = telev
@@ -1914 +1914 @@
                                   use_cache = False,
                                   verbose = False):
-    """Read Digitial Elevation model from the following ASCII format (.asc)    
+    """Read Digitial Elevation model from the following ASCII format (.asc)
 
     Example:
@@ -1958 +1958 @@
     else:
         result = apply(_convert_dem_from_ascii2netcdf, [basename_in], kwargs)
-        
+
     return result
 
-    
-    
+
+
 
 
@@ -2161 +2161 @@
     if maxlon != None:
         assert -180 < maxlon < 180 , msg
- 
+
 
     #Get NetCDF data
@@ -2168 +2168 @@
     file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
     file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)
-    file_e = NetCDFFile(basename_in + '_e.nc', 'r') #Elevation (z) (m)
+    file_e = NetCDFFile(basename_in + '_e.nc', 'r')  #Elevation (z) (m)
 
     if basename_out is None:
@@ -2605 +2605 @@
     from Numeric import array
     from config import time_format
-    from utilities.numerical_tools import ensure_numeric    
+    from utilities.numerical_tools import ensure_numeric
 
 
@@ -2823 +2823 @@
 
 
-    try:    
+    try:
         domain = Domain(coordinates, volumes, boundary)
     except AssertionError, e:
@@ -2887 +2887 @@
             X = resize(X,(X.shape[0]/3,3))
         domain.set_quantity(quantity,X)
-        
+
     fid.close()
     return domain
@@ -3592 +3592 @@
     quantities which can be derived from those, i.e.
     ['depth', 'xmomentum', 'ymomentum', 'momentum', 'velocity', 'bearing'].
-    
+
     Momentum is the absolute momentum, sqrt(xmomentum^2 + ymomentum^2).
     Note, units of momentum are m^2/s and depth is m.
@@ -3611 +3611 @@
 
     # extract gauge locations from gauge file
-    
+
     # extract all quantities from swwfile (f = file_function)
     if time_min is None: time_min = min(f.T)
     if time_max is None: time_max = max(f.T)
-    
+
     # loop through appropriate range of time
     # plot prescribed quantities and export data if requested
-    
+
     #if gauge_data_outname is None:
 

inundation/pyvolution/domain.py

@@ -78 +78 @@
         #FIXME: Maybe have separate orders for h-limiter and w-limiter?
         #Or maybe get rid of order altogether and use beta_w and beta_h
-        self.set_default_order(1)        
+        self.set_default_order(1)
         #self.default_order = 1
         #self.order = self.default_order
@@ -115 +115 @@
         if mesh_filename is not None:
             # If the mesh file passed any quantity values
-            # , initialise with these values. 
+            # , initialise with these values.
             self.set_quantity_vertices_dict(vertex_quantity_dict)
 
-        
+
 
     def set_default_order(self, n):
@@ -126 +126 @@
         msg = 'Default order must be either 1 or 2. I got %s' %n
         assert n in [1,2], msg
-        
+
         self.default_order = n
         self.order = self.default_order
@@ -412 +412 @@
         Input:
           quantities: either None, a string or a list of strings naming the quantities to be reported
-          tags:       either None, a string or a list of strings naming the tags to be reported          
+          tags:       either None, a string or a list of strings naming the tags to be reported
 
 
@@ -514 +514 @@
                yieldstep = None,
                finaltime = None,
-               duration = None,              
+               duration = None,
                skip_initial_step = False):
         """Evolve model through time starting from self.starttime.
-        
-       
+
+
         yieldstep: Interval between yields where results are stored,
                    statistics written and domain inspected or
                    possibly modified. If omitted the internal predefined
                    max timestep is used.
-                   Internally, smaller timesteps may be taken.               
+                   Internally, smaller timesteps may be taken.
 
         duration: Duration of simulation
@@ -534 +534 @@
         skip_initial_step: Boolean flag that decides whether the first
         yield step is skipped or not. This is useful for example to avoid
-        duplicate steps when multiple evolve processes are dove tailed. 
+        duplicate steps when multiple evolve processes are dove tailed.
 
 
@@ -543 +543 @@
 
 
-        All times are given in seconds        
+        All times are given in seconds
 
         """
@@ -575 +575 @@
             if duration is not None:
                 self.finaltime = self.starttime + float(duration)
-                
-
-        
+
+
+
 
         self.yieldtime = 0.0 #Time between 'yields'
@@ -678 +678 @@
             if B is None:
                 print 'WARNING: Ignored boundary segment %d (None)'
-            else:               
+            else:
                 q = B.evaluate(vol_id, edge_id)
 
@@ -693 +693 @@
     def update_timestep(self, yieldstep, finaltime):
 
-        from config import min_timestep
+        from config import min_timestep, max_timestep
 
         # self.timestep is calculated from speed of characteristics
         # Apply CFL condition here
-        timestep = self.CFL*self.timestep
+        timestep = min(self.CFL*self.timestep,max_timestep)
 
         #Record maximal and minimal values of timestep for reporting

inundation/pyvolution/quantity.py

@@ -197 +197 @@
           Arbitrary geo spatial dataset in the form of the class
           Geospatial_data. Mesh points are populated using least squares
-          fitting           
+          fitting
 
         points:
@@ -210 +210 @@
 
         data_georef:
-          If points is specified, geo_reference applies to each point.        
-          
+          If points is specified, geo_reference applies to each point.
+
         filename:
           Name of a .pts file containing data points and attributes for
           use with least squares.
-          
+
         attribute_name:
           If specified, any array matching that name
@@ -289 +289 @@
         if not(points is None and values is None and data_georef is None):
             from warnings import warn
-            
+
             msg = 'Using points, values or data_georef with set_quantity '
             msg += 'is obsolete. Please use a Geospatial_data object instead.'
@@ -312 +312 @@
                                               location, indices, verbose)
             elif isinstance(numeric, Geospatial_data):
-                self.set_values_from_geospatial_data(numeric,                
+                self.set_values_from_geospatial_data(numeric,
                                                      alpha,
                                                      location, indices,
@@ -338 +338 @@
             print 'The usage of points in set_values will be deprecated.' +\
                   'Please use the geospatial_data object.'
-            
+
             msg = 'When points are specified, associated values must also be.'
             assert values is not None, msg
@@ -350 +350 @@
                                       location, indices,
                                       verbose = verbose,
-                                      use_cache = use_cache) 
+                                      use_cache = use_cache)
         else:
             raise 'This can\'t happen :-)'
@@ -586 +586 @@
 
 
-        
+
     def set_values_from_geospatial_data(self, geospatial_data, alpha,
                                         location, indices,
@@ -604 +604 @@
                                     data_georef = data_georef,
                                     verbose = verbose,
-                                    use_cache = use_cache)                                    
-
-                               
+                                    use_cache = use_cache)
+
+
 
     def set_values_from_points(self, points, values, alpha,
@@ -621 +621 @@
         from least_squares import fit_to_mesh
         from coordinate_transforms.geo_reference import Geo_reference
-        
+
 
         points = ensure_numeric(points, Float)
@@ -645 +645 @@
         #print data_georef
         #print points
-            
-
-        #Call least squares method            
+
+
+        #Call least squares method
         args = (coordinates, triangles, points, values)
         kwargs = {'data_origin': data_georef.get_origin(),
@@ -656 +656 @@
 
         #print kwargs
-        
+
         if use_cache is True:
             try:
@@ -672 +672 @@
             vertex_attributes = apply(fit_to_mesh,
                                       args, kwargs)
-            
-        #Call underlying method using array values    
+
+        #Call underlying method using array values
         self.set_values_from_array(vertex_attributes,
                                    location, indices, verbose)
@@ -730 +730 @@
             data_georef = None
 
-        
+
 
         #Call underlying method for geospatial data
         geospatial_data = points_dictionary2geospatial_data(points_dict)
         geospatial_data.set_default_attribute_name(attribute_name)
-            
+
         self.set_values_from_geospatial_data(geospatial_data,
                                              alpha,
@@ -741 +741 @@
                                              verbose = verbose,
                                              use_cache = use_cache)
-        
+
         #Call underlying method for points
         #self.set_values_from_points(points, z, alpha,
@@ -770 +770 @@
 
         The values will be stored in elements following their
-        internal ordering.        
+        internal ordering.
 
         """
@@ -1132 +1132 @@
     denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
 
+    if sum(less(denominator, 1.0)) > 0.0:
+        msg = 'denominator < 1.0 in semi implicit update. Call Stephen :-)'
+        raise msg
+
     if sum(equal(denominator, 0.0)) > 0.0:
         msg = 'Zero division in semi implicit update. Call Stephen :-)'

inundation/pyvolution/quantity_ext.c

@@ -81 +81 @@
 
       //Two point gradient
-      _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);      
-      
-      
+      _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);
+
+
       //Old (wrong code)
       //det = x0*y1 - x1*y0;
@@ -187 +187 @@
 
 
-  //Explicit updates
-  for (k=0; k<N; k++) {
-    centroid_values[k] += timestep*explicit_update[k];
-  }
-
   //Semi implicit updates
   for (k=0; k<N; k++) {
@@ -202 +197 @@
     }
   }
+
+
+  //Explicit updates
+  for (k=0; k<N; k++) {
+    centroid_values[k] += timestep*explicit_update[k];
+  }
+
+
   //MH080605 set semi_impliit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
-  for (k=0;k<N;k++)
+  for (k=0;k<N;k++){
     semi_implicit_update[k]=0.0;
+  }
+
   return 0;
 }

inundation/pyvolution/shallow_water.py

@@ -83 +83 @@
         from config import minimum_allowed_height, maximum_allowed_speed, g
         self.minimum_allowed_height = minimum_allowed_height
-        self.maximum_allowed_speed = maximum_allowed_speed      
+        self.maximum_allowed_speed = maximum_allowed_speed
         self.g = g
 
+
+        self.forcing_terms.append(manning_friction)
         self.forcing_terms.append(gravity)
-        self.forcing_terms.append(manning_friction)
 
         #Realtime visualisation
@@ -108 +109 @@
         self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
 
-        
+
     def set_quantities_to_be_stored(self, q):
         """Specify which quantities will be stored in the sww file.
@@ -118 +119 @@
 
         In the two first cases, the named quantities will be stored at each yieldstep
-        (This is in addition to the quantities elevation and friction)  
+        (This is in addition to the quantities elevation and friction)
         If q is None, storage will be switched off altogether.
         """
@@ -124 +125 @@
 
         if q is None:
-            self.quantities_to_be_stored = []            
+            self.quantities_to_be_stored = []
             self.store = False
             return
@@ -131 +132 @@
             q = [q] # Turn argument into a list
 
-        #Check correcness    
+        #Check correcness
         for quantity_name in q:
             msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
-            assert quantity_name in self.conserved_quantities, msg 
-        
+            assert quantity_name in self.conserved_quantities, msg
+
         self.quantities_to_be_stored = q
-        
+
 
     def initialise_visualiser(self,scale_z=1.0,rect=None):
@@ -224 +225 @@
                yieldstep = None,
                finaltime = None,
-               duration = None,               
+               duration = None,
                skip_initial_step = False):
         """Specialisation of basic evolve method from parent class
@@ -670 +671 @@
     hc = wc - zc  #Water depths at centroids
 
-    #Update 
+    #Update
     #FIXME: Modify accroditg to c-version - or discard altogether.
     for k in range(domain.number_of_elements):
@@ -695 +696 @@
     from shallow_water_ext import protect
 
-    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed, 
+    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
             domain.epsilon, wc, zc, xmomc, ymomc)
 
@@ -1183 +1184 @@
 
 
-def manning_friction_c(domain):
+def manning_friction_implicit_c(domain):
     """Wrapper for c version
     """
@@ -1200 +1201 @@
     xmom_update = xmom.semi_implicit_update
     ymom_update = ymom.semi_implicit_update
+
+    N = domain.number_of_elements
+    eps = domain.minimum_allowed_height
+    g = domain.g
+
+    from shallow_water_ext import manning_friction
+    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
+
+
+def manning_friction_explicit_c(domain):
+    """Wrapper for c version
+    """
+
+    xmom = domain.quantities['xmomentum']
+    ymom = domain.quantities['ymomentum']
+
+    w = domain.quantities['stage'].centroid_values
+    z = domain.quantities['elevation'].centroid_values
+
+    uh = xmom.centroid_values
+    vh = ymom.centroid_values
+    eta = domain.quantities['friction'].centroid_values
+
+    xmom_update = xmom.explicit_update
+    ymom_update = ymom.explicit_update
 
     N = domain.number_of_elements
@@ -1791 +1817 @@
     extrapolate_second_order_sw=extrapolate_second_order_sw_c
     gravity = gravity_c
-    manning_friction = manning_friction_c
+    manning_friction = manning_friction_implicit_c
     h_limiter = h_limiter_c
     balance_deep_and_shallow = balance_deep_and_shallow_c

inundation/pyvolution/shallow_water_ext.c

@@ -260 +260 @@
 
 
+void _manning_friction_explicit(double g, double eps, int N,
+                       double* w, double* z,
+                       double* uh, double* vh,
+                       double* eta, double* xmom, double* ymom) {
+
+  int k;
+  double S, h;
+
+  for (k=0; k<N; k++) {
+    if (eta[k] > eps) {
+      h = w[k]-z[k];
+      if (h >= eps) {
+        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
+        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
+        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
+        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
+
+
+        //Update momentum
+        xmom[k] += S*uh[k];
+        ymom[k] += S*vh[k];
+      }
+    }
+  }
+}
+
 
 int _balance_deep_and_shallow(int N,
@@ -356 +382 @@
 int _protect(int N,
              double minimum_allowed_height,
-             double maximum_allowed_speed,           
+             double maximum_allowed_speed,
              double epsilon,
              double* wc,
@@ -365 +391 @@
   int k;
   double hc;
-  double u, v, reduced_speed;  
+  double u, v, reduced_speed;
 
   //Protect against initesimal and negative heights
@@ -372 +398 @@
 
     if (hc < minimum_allowed_height) {
-    
+
       //Old code: Set momentum to zero and ensure h is non negative
       //xmomc[k] = 0.0;
-      //ymomc[k] = 0.0;      
-      //if (hc <= 0.0) wc[k] = zc[k]; 
-      
+      //ymomc[k] = 0.0;
+      //if (hc <= 0.0) wc[k] = zc[k];
+
 
       //New code: Adjust momentum to guarantee speeds are physical
-      //          ensure h is non negative      
-      //FIXME (Ole): This is only implemented in this C extension and 
-      //             has no Python equivalent      
+      //          ensure h is non negative
+      //FIXME (Ole): This is only implemented in this C extension and
+      //             has no Python equivalent
       if (hc <= 0.0) {
-        wc[k] = zc[k]; 
+        wc[k] = zc[k];
         xmomc[k] = 0.0;
         ymomc[k] = 0.0;
       } else {
-        //Reduce excessive speeds derived from division by small hc 
+        //Reduce excessive speeds derived from division by small hc
         u = xmomc[k]/hc;
         if (fabs(u) > maximum_allowed_speed) {
           reduced_speed = maximum_allowed_speed * u/fabs(u);
-          //printf("Speed (u) has been reduced from %.3f to %.3f\n", 
-          //     u, reduced_speed);     
+          //printf("Speed (u) has been reduced from %.3f to %.3f\n",
+          //     u, reduced_speed);
           xmomc[k] = reduced_speed * hc;
         }
-        
-        v = ymomc[k]/hc;        
+
+        v = ymomc[k]/hc;
         if (fabs(v) > maximum_allowed_speed) {
-          reduced_speed = maximum_allowed_speed * v/fabs(v);    
-          //printf("Speed (v) has been reduced from %.3f to %.3f\n",   
-          //     v, reduced_speed);     
-          ymomc[k] = reduced_speed * hc;          
-        }       
+          reduced_speed = maximum_allowed_speed * v/fabs(v);
+          //printf("Speed (v) has been reduced from %.3f to %.3f\n",
+          //     v, reduced_speed);
+          ymomc[k] = reduced_speed * hc;
+        }
       }
     }
@@ -530 +556 @@
 }
 
+
+PyObject *manning_friction_explicit(PyObject *self, PyObject *args) {
+  //
+  // manning_friction_explicit(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
+  //
+
+
+  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
+  int N;
+  double g, eps;
+
+  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
+                        &g, &eps, &w, &z, &uh, &vh, &eta,
+                        &xmom, &ymom))
+    return NULL;
+
+  N = w -> dimensions[0];
+  _manning_friction_explicit(g, eps, N,
+                    (double*) w -> data,
+                    (double*) z -> data,
+                    (double*) uh -> data,
+                    (double*) vh -> data,
+                    (double*) eta -> data,
+                    (double*) xmom -> data,
+                    (double*) ymom -> data);
+
+  return Py_BuildValue("");
+}
 PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
   /*Compute the vertex values based on a linear reconstruction on each triangle
@@ -1067 +1121 @@
   if (!PyArg_ParseTuple(args, "dddOOOO",
                         &minimum_allowed_height,
-                        &maximum_allowed_speed,                 
+                        &maximum_allowed_speed,
                         &epsilon,
                         &wc, &zc, &xmomc, &ymomc))

inundation/pyvolution/test_data_manager.py

@@ -79 +79 @@
         latitudes = [-34.5, -34.33333, -34.16667, -34]
 
-        long_name = 'LON'
+        long_name = 'LON'
         lat_name = 'LAT'
 
@@ -89 +89 @@
         for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
             fid = NetCDFFile(self.test_MOST_file + ext, 'w')
-       
+
             fid.createDimension(long_name,nx)
             fid.createVariable(long_name,'d',(long_name,))
@@ -107 +107 @@
             fid.variables['TIME'].units='seconds'
             fid.variables['TIME'].assignValue([0.0, 0.1, 0.6, 1.1, 1.6, 2.1])
-            
+
 
             name = ext[1:3].upper()
@@ -113 +113 @@
             fid.createVariable(name,'d',('TIME', lat_name, long_name))
             fid.variables[name].units='CENTIMETERS'
-            fid.variables[name].missing_value=-1.e+034        
-            
+            fid.variables[name].missing_value=-1.e+034
+
             fid.variables[name].assignValue([[[0.3400644, 0, -46.63519, -6.50198],
                                               [-0.1214216, 0, 0, 0],
@@ -139 +139 @@
                                               [0.0004811212, 0.0009622425, 0.0009599366, 8.152277e-007],
                                               [0, 0.0004811212, 0.0004811212, 0]]])
-            
-        
+
+
             fid.close()
 
-        
+
 
 
@@ -150 +150 @@
         for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
             #print 'Trying to remove', self.test_MOST_file + ext
-            os.remove(self.test_MOST_file + ext)        
+            os.remove(self.test_MOST_file + ext)
 
     def test_sww_constant(self):
@@ -520 +520 @@
     #def test_write_pts(self):
     #    #Obsolete
-    #    
+    #
     #    #Get (enough) datapoints
     #
@@ -546 +546 @@
     #    #Check contents
     #    #Get NetCDF
-    #    from Scientific.IO.NetCDF import NetCDFFile        
+    #    from Scientific.IO.NetCDF import NetCDFFile
     #    fid = NetCDFFile(ptsfile, 'r')
     #
@@ -555 +555 @@
     #
     #    #Check values#
-    # 
+    #
     #    #print points[:]
     #    #print ref_points
-    #    assert allclose(points, points1)  
+    #    assert allclose(points, points1)
     #
     #    #print attributes[:]
@@ -569 +569 @@
     #    import os
     #    os.remove(ptsfile)
-        
+
 
     def test_dem2pts_bounding_box_v2(self):
@@ -591 +591 @@
 NODATA_value  -9999
 """)
-        #Create linear function        
+        #Create linear function
         ref_points = []
         ref_elevation = []
@@ -610 +610 @@
 
         fid.close()
-            
+
         #print 'sending pts', ref_points
         #print 'sending elev', ref_elevation
@@ -650 +650 @@
 
         #create new reference points
-        newz = []        
+        newz = []
         newz[0:5] = ref_elevation[32:38]
         newz[6:11] = ref_elevation[42:48]
@@ -668 +668 @@
                 x = x0 + xvec[j]
                 xnew = x - 2002.0
-                ref_points.append ([xnew,ynew]) #Relative point values        
+                ref_points.append ([xnew,ynew]) #Relative point values
 
         assert allclose(points, ref_points)
@@ -704 +704 @@
 NODATA_value  -9999
 """)
-        #Create linear function        
+        #Create linear function
         ref_points = []
         ref_elevation = []
@@ -721 +721 @@
                 ref_points.append ([x,y])
                 count += 1
-                z[count] = (4*i - 3*j)%13 
+                z[count] = (4*i - 3*j)%13
                 if j == 4: z[count] = NODATA_value #column inside clipping region
                 if j == 8: z[count] = NODATA_value #column outside clipping region
                 if i == 9: z[count] = NODATA_value #row outside clipping region
-                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region                
+                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
                 ref_elevation.append( z[count] )
                 fid.write('%f ' %z[count])
             fid.write('\n')
 
-        fid.close()             
-        
+        fid.close()
+
         #print 'sending elev', ref_elevation
 
@@ -781 +781 @@
         newz[16:19] = ref_elevation[65:68]
 
-        
+
         ref_elevation = newz
         ref_points = []
@@ -795 +795 @@
                 x = x0 + xvec[j]
                 xnew = x - 2002.0
-                if j <> 2 and (i<>1 or j<>4):    
+                if j <> 2 and (i<>1 or j<>4):
                     ref_points.append([x,y])
                     new_ref_points.append ([xnew,ynew])
 
-               
+
         assert allclose(points, new_ref_points)
         assert allclose(elevation, ref_elevation)
@@ -834 +834 @@
 NODATA_value  -9999
 """)
-        #Create linear function        
+        #Create linear function
         ref_points = []
         ref_elevation = []
@@ -851 +851 @@
                 ref_points.append ([x,y])
                 count += 1
-                z[count] = (4*i - 3*j)%13 
+                z[count] = (4*i - 3*j)%13
                 if j == 4: z[count] = NODATA_value #column inside clipping region
                 if j == 8: z[count] = NODATA_value #column outside clipping region
                 if i == 6: z[count] = NODATA_value #row on clipping boundary
-                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region                
+                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
                 ref_elevation.append( z[count] )
                 fid.write('%f ' %z[count])
             fid.write('\n')
 
-        fid.close()             
-        
+        fid.close()
+
         #print 'sending elev', ref_elevation
 
@@ -910 +910 @@
 
 
-        
+
         ref_elevation = newz
         ref_points = []
@@ -924 +924 @@
                 x = x0 + xvec[j]
                 xnew = x - 2002.0
-                if j <> 2 and (i<>1 or j<>4) and i<>3:    
+                if j <> 2 and (i<>1 or j<>4) and i<>3:
                     ref_points.append([x,y])
                     new_ref_points.append ([xnew,ynew])
@@ -932 +932 @@
         #print new_ref_points, len(new_ref_points)
 
-        assert allclose(elevation, ref_elevation)                    
+        assert allclose(elevation, ref_elevation)
         assert allclose(points, new_ref_points)
 
@@ -973 +973 @@
   PROJECTION ZONE: 50
   DATUM: AGD66
-  VERTICAL DATUM: 
-  NUMBER OF REACHES:  19 
-  NUMBER OF CROSS-SECTIONS:  2 
+  VERTICAL DATUM:
+  NUMBER OF REACHES:  19
+  NUMBER OF CROSS-SECTIONS:  2
 END HEADER:
 
@@ -984 +984 @@
     STREAM ID:Southern-Wungong
     REACH ID:Southern-Wungong
-    STATION:21410   
+    STATION:21410
     CUT LINE:
-      407546.08 , 6437277.542 
-      407329.32 , 6437489.482 
-      407283.11 , 6437541.232 
+      407546.08 , 6437277.542
+      407329.32 , 6437489.482
+      407283.11 , 6437541.232
     SURFACE LINE:
      407546.08,   6437277.54,   52.14
@@ -999 +999 @@
 
   CROSS-SECTION:
-    STREAM ID:Swan River      
-    REACH ID:Swan Mouth      
-    STATION:840.*   
+    STREAM ID:Swan River
+    REACH ID:Swan Mouth
+    STATION:840.*
     CUT LINE:
-      381178.0855 , 6452559.0685 
-      380485.4755 , 6453169.272 
+      381178.0855 , 6452559.0685
+      380485.4755 , 6453169.272
     SURFACE LINE:
      381178.09,   6452559.07,   4.17
@@ -1017 +1017 @@
      381063.46,   6452660.06,   3.67
      381054.41,   6452668.03,   3.67
-  END:  
+  END:
 END CROSS-SECTIONS:
 """)
@@ -1044 +1044 @@
                       [407510.08, 6437312.74]]
 
-        ref_points += [[381178.09, 6452559.07], 
-                       [381169.49, 6452566.64], 
-                       [381157.78, 6452576.96], 
-                       [381155.97, 6452578.56], 
-                       [381143.72, 6452589.35], 
-                       [381136.69, 6452595.54], 
-                       [381114.74, 6452614.88], 
-                       [381075.53, 6452649.43], 
-                       [381071.47, 6452653.00], 
-                       [381063.46, 6452660.06], 
-                       [381054.41, 6452668.03]] 
-        
-                      
+        ref_points += [[381178.09, 6452559.07],
+                       [381169.49, 6452566.64],
+                       [381157.78, 6452576.96],
+                       [381155.97, 6452578.56],
+                       [381143.72, 6452589.35],
+                       [381136.69, 6452595.54],
+                       [381114.74, 6452614.88],
+                       [381075.53, 6452649.43],
+                       [381071.47, 6452653.00],
+                       [381063.46, 6452660.06],
+                       [381054.41, 6452668.03]]
+
+
         ref_elevation = [52.14, 51.07, 50.56, 49.58, 49.09, 48.76]
         ref_elevation += [4.17, 4.26, 4.34, 4.35, 4.43, 4.58, 4.41, 4.17, 3.99, 3.67, 3.67]
@@ -1127 +1127 @@
                 verbose = False,
                 format = 'asc')
-                
+
         #Check prj (meta data)
         prjid = open(prjfile)
@@ -1225 +1225 @@
         cellsize      10.000000
         NODATA_value  -9999
-        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200 
-         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 
-         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 
-         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 
-         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 
-         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 
+        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
+         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
+         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
+         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
+         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
+         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
          -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
          -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
@@ -1236 +1236 @@
          -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
            0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
-        
+
         """
 
@@ -1267 +1267 @@
         #
         domain.set_quantity('elevation', lambda x,y: -x-y)
-        domain.set_quantity('stage', 0)        
+        domain.set_quantity('stage', 0)
 
         B = Transmissive_boundary(domain)
@@ -1283 +1283 @@
         cellsize = 10  #10m grid
 
-        
+
         #Check contents
         #Get NetCDF
@@ -1303 +1303 @@
                 verbose = False,
                 format = 'asc')
-        
-                
+
+
         #Check prj (meta data)
         prjid = open(prjfile)
@@ -1380 +1380 @@
                 #assert allclose(float(value), -(10-i+j)*cellsize)
                 assert float(value) == -(10-i+j)*cellsize
-        
+
 
         fid.close()
@@ -1398 +1398 @@
         Original extent is 100m x 100m:
 
-        Eastings:   308500 -  308600 
+        Eastings:   308500 -  308600
         Northings: 6189000 - 6189100
 
         Bounding box changes this to the 50m x 50m square defined by
 
-        Eastings:   308530 -  308570 
+        Eastings:   308530 -  308570
         Northings: 6189050 - 6189100
 
         The cropped values should be
 
-         -130 -140 -150 -160 -170  
-         -120 -130 -140 -150 -160  
-         -110 -120 -130 -140 -150  
-         -100 -110 -120 -130 -140  
-          -90 -100 -110 -120 -130  
+         -130 -140 -150 -160 -170
+         -120 -130 -140 -150 -160
+         -110 -120 -130 -140 -150
+         -100 -110 -120 -130 -140
+          -90 -100 -110 -120 -130
           -80  -90 -100 -110 -120
 
-        and the new lower reference point should be   
-        Eastings:   308530 
+        and the new lower reference point should be
+        Eastings:   308530
         Northings: 6189050
-        
+
         Original dataset is the same as in test_sww2dem_larger()
-        
+
         """
 
@@ -1451 +1451 @@
         #
         domain.set_quantity('elevation', lambda x,y: -x-y)
-        domain.set_quantity('stage', 0)        
+        domain.set_quantity('stage', 0)
 
         B = Transmissive_boundary(domain)
@@ -1467 +1467 @@
         cellsize = 10  #10m grid
 
-        
+
         #Check contents
         #Get NetCDF
@@ -1487 +1487 @@
                 easting_min = 308530,
                 easting_max = 308570,
-                northing_min = 6189050, 
-                northing_max = 6189100,                                
+                northing_min = 6189050,
+                northing_max = 6189100,
                 verbose = False,
                 format = 'asc')
-        
-        fid.close()
-
-        
+
+        fid.close()
+
+
         #Check prj (meta data)
         prjid = open(prjfile)
@@ -1568 +1568 @@
         for i, line in enumerate(lines[6:]):
             for j, value in enumerate( line.split() ):
-                #assert float(value) == -(10-i+j)*cellsize                
+                #assert float(value) == -(10-i+j)*cellsize
                 assert float(value) == -(10-i+j+3)*cellsize
-        
+
 
 
@@ -2000 +2000 @@
         #        cellsize = cellsize,
         #        verbose = False)
-                
+
         sww2dem(self.domain.filename,
                 quantity = 'elevation',
@@ -2010 +2010 @@
         #Check header data
         from ermapper_grids import read_ermapper_header, read_ermapper_data
-        
+
         header = read_ermapper_header(self.domain.filename + '_elevation.ers')
         #print header
         assert header['projection'].lower() == '"utm-56"'
         assert header['datum'].lower() == '"wgs84"'
-        assert header['units'].lower() == '"meters"'    
-        assert header['value'].lower() == '"elevation"'         
+        assert header['units'].lower() == '"meters"'
+        assert header['value'].lower() == '"elevation"'
         assert header['xdimension'] == '0.25'
-        assert header['ydimension'] == '0.25'   
+        assert header['ydimension'] == '0.25'
         assert float(header['eastings']) == 308500.0   #xllcorner
-        assert float(header['northings']) == 6189000.0 #yllcorner       
+        assert float(header['northings']) == 6189000.0 #yllcorner
         assert int(header['nroflines']) == 5
-        assert int(header['nrofcellsperline']) == 5     
+        assert int(header['nrofcellsperline']) == 5
         assert int(header['nullcellvalue']) == NODATA_value
-        #FIXME - there is more in the header                    
-
-                
-        #Check grid data                
-        grid = read_ermapper_data(self.domain.filename + '_elevation') 
-        
+        #FIXME - there is more in the header
+
+
+        #Check grid data
+        grid = read_ermapper_data(self.domain.filename + '_elevation')
+
         #FIXME (Ole): Why is this the desired reference grid for -x-y?
         ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
                     -1,    -1.25, -1.5,  -1.75, -2.0,
-                    -0.75, -1.0,  -1.25, -1.5,  -1.75,              
+                    -0.75, -1.0,  -1.25, -1.5,  -1.75,
                     -0.5,  -0.75, -1.0,  -1.25, -1.5,
-                    -0.25, -0.5,  -0.75, -1.0,  -1.25]              
-                                          
+                    -0.25, -0.5,  -0.75, -1.0,  -1.25]
+
 
         #print grid
@@ -2042 +2042 @@
 
         fid.close()
-        
+
         #Cleanup
         #FIXME the file clean-up doesn't work (eg Permission Denied Error)
-        #Done (Ole) - it was because sww2ers didn't close it's sww file 
+        #Done (Ole) - it was because sww2ers didn't close it's sww file
         os.remove(sww.filename)
         os.remove(self.domain.filename + '_elevation')
-        os.remove(self.domain.filename + '_elevation.ers')                 
+        os.remove(self.domain.filename + '_elevation.ers')
 
 
@@ -2067 +2067 @@
         # Fourth value (index==3) is -6.50198 cm
 
-        
+
 
         #Read
         from coordinate_transforms.redfearn import redfearn
         #fid = NetCDFFile(self.test_MOST_file)
-        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')        
+        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
         first_value = fid.variables['HA'][:][0,0,0]
         fourth_value = fid.variables['HA'][:][0,0,3]
@@ -2090 +2090 @@
         #Read output file 'small.sww'
         #fid = NetCDFFile('small.sww')
-        fid = NetCDFFile(self.test_MOST_file + '.sww')        
+        fid = NetCDFFile(self.test_MOST_file + '.sww')
 
         x = fid.variables['x'][:]
@@ -2211 +2211 @@
         h2_list = [-34.5,-34.33333]
 
-        long_name = 'LON'
+        long_name = 'LON'
         lat_name = 'LAT'
-        time_name = 'TIME'        
+        time_name = 'TIME'
 
         nx = 3
@@ -2253 +2253 @@
 
         name = {}
-        name[fid1]='HA'
-        name[fid2]='UA'
-        name[fid3]='VA'
-        name[fid4]='ELEVATION'
+        name[fid1]='HA'
+        name[fid2]='UA'
+        name[fid3]='VA'
+        name[fid4]='ELEVATION'
 
         units = {}
-        units[fid1]='cm'
-        units[fid2]='cm/s'
-        units[fid3]='cm/s'
-        units[fid4]='m'
-
-        values = {}
-        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
-        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
-        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
-        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
+        units[fid1]='cm'
+        units[fid2]='cm/s'
+        units[fid3]='cm/s'
+        units[fid4]='m'
+
+        values = {}
+        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
+        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
+        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
+        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
 
         for fid in [fid1,fid2,fid3]:
@@ -2276 +2276 @@
           fid.variables[name[fid]].assignValue(values[fid])
           fid.variables[name[fid]].missing_value = -99999999.
-          if fid == fid3: break
+          if fid == fid3: break
 
         for fid in [fid4]:
-            fid.createVariable(name[fid],'d',(lat_name,long_name))
+            fid.createVariable(name[fid],'d',(lat_name,long_name))
             fid.variables[name[fid]].point_spacing='uneven'
             fid.variables[name[fid]].units=units[fid]
@@ -2323 +2323 @@
         elevation = fid.variables['elevation'][:]
         stage = fid.variables['stage'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]
-
-        #print ymomentum
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+
+        #print ymomentum
         first_height = first_amp/100 - first_elevation
         third_height = third_amp/100 - third_elevation
@@ -2371 +2371 @@
         h2_list = [-34.5,-34.33333]
 
-        long_name = 'LON961_1261'
-        lat_name = 'LAT481_841'
-        time_name = 'TIME1'
+#        long_name = 'LON961_1261'
+#        lat_name = 'LAT481_841'
+#        time_name = 'TIME1'
+
+        long_name = 'LON'
+        lat_name = 'LAT'
+        time_name = 'TIME'
 
         nx = 3
@@ -2413 +2417 @@
 
         name = {}
-        name[fid1]='HA'
-        name[fid2]='UA'
-        name[fid3]='VA'
-        name[fid4]='ELEVATION'
+        name[fid1]='HA'
+        name[fid2]='UA'
+        name[fid3]='VA'
+        name[fid4]='ELEVATION'
 
         units = {}
-        units[fid1]='cm'
-        units[fid2]='cm/s'
-        units[fid3]='cm/s'
-        units[fid4]='m'
-
-        values = {}
-        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
-        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
-        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
-        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
+        units[fid1]='cm'
+        units[fid2]='cm/s'
+        units[fid3]='cm/s'
+        units[fid4]='m'
+
+        values = {}
+        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
+        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
+        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
+        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
 
         for fid in [fid1,fid2,fid3]:
@@ -2436 +2440 @@
           fid.variables[name[fid]].assignValue(values[fid])
           fid.variables[name[fid]].missing_value = -99999999.
-          if fid == fid3: break
+          if fid == fid3: break
 
         for fid in [fid4]:
@@ -2468 +2472 @@
 
         #Call conversion (with zero origin)
-        ferret2sww('test', verbose=False,
-                   origin = (56, 0, 0))
+        ferret2sww('test', verbose=False, origin = (56, 0, 0))
 
         os.remove('test_va.nc')
@@ -2483 +2486 @@
         elevation = fid.variables['elevation'][:]
         stage = fid.variables['stage'][:]
-        xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]
-
-        #print ymomentum
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+
+        #print ymomentum
         first_height = first_amp/100 - first_elevation
         third_height = third_amp/100 - third_elevation
@@ -2518 +2521 @@
         #fid = NetCDFFile('small.sww', 'r')
         fid = NetCDFFile(self.test_MOST_file + '.sww')
-        
+
         x = fid.variables['x'][:]
         y = fid.variables['y'][:]
@@ -2642 +2645 @@
 
         #os.remove(domain.filename + '.sww')
-        os.remove(filename)        
+        os.remove(filename)
 
         bits = ['vertex_coordinates']
@@ -2782 +2785 @@
         #Clean up
         os.remove(filename)
-            
+
 
         bits = [ 'geo_reference.get_xllcorner()',
@@ -3204 +3207 @@
 
         os.remove(filename)
-        
+
     def test_asc_csiro2sww(self):
         import tempfile
 
         bath_dir = tempfile.mkdtemp()
-        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
+        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
         #bath_dir = 'bath_data_manager_test'
-        #print "os.getcwd( )",os.getcwd( ) 
+        #print "os.getcwd( )",os.getcwd( )
         elevation_dir =  tempfile.mkdtemp()
         #elevation_dir = 'elev_expanded'
-        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
-        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
-        
+        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
+        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
+
         fid = open(bath_dir_filename, 'w')
         fid.write(""" ncols             3
@@ -3225 +3228 @@
  nodata_value   -9999.0
     9000.000    -1000.000    3000.0
-   -1000.000    9000.000  -1000.000 
+   -1000.000    9000.000  -1000.000
 """)
         fid.close()
-        
+
         fid = open(elevation_dir_filename1, 'w')
         fid.write(""" ncols             3
@@ -3237 +3240 @@
  nodata_value   -9999.0
     9000.000    0.000    3000.0
-     0.000     9000.000     0.000 
+     0.000     9000.000     0.000
 """)
         fid.close()
@@ -3249 +3252 @@
  nodata_value   -9999.0
     9000.000    4000.000    4000.0
-    4000.000    9000.000    4000.000 
+    4000.000    9000.000    4000.000
 """)
         fid.close()
-         
+
         ucur_dir =  tempfile.mkdtemp()
         ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
         ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
-        
+
         fid = open(ucur_dir_filename1, 'w')
         fid.write(""" ncols             3
@@ -3265 +3268 @@
  nodata_value   -9999.0
     90.000    60.000    30.0
-    10.000    10.000    10.000 
+    10.000    10.000    10.000
 """)
         fid.close()
@@ -3276 +3279 @@
  nodata_value   -9999.0
     90.000    60.000    30.0
-    10.000    10.000    10.000 
+    10.000    10.000    10.000
 """)
         fid.close()
-        
+
         vcur_dir =  tempfile.mkdtemp()
         vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
         vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
-        
+
         fid = open(vcur_dir_filename1, 'w')
         fid.write(""" ncols             3
@@ -3292 +3295 @@
  nodata_value   -9999.0
     90.000    60.000    30.0
-    10.000    10.000    10.000 
+    10.000    10.000    10.000
 """)
         fid.close()
@@ -3303 +3306 @@
  nodata_value   -9999.0
     90.000    60.000    30.0
-    10.000    10.000    10.000 
+    10.000    10.000    10.000
 """)
         fid.close()
-        
+
         sww_file = 'a_test.sww'
         asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file)
 
         # check the sww file
-        
+
         fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
         x = fid.variables['x'][:]
@@ -3319 +3322 @@
         xmomentum = fid.variables['xmomentum'][:]
         geo_ref = Geo_reference(NetCDFObject=fid)
-        #print "geo_ref",geo_ref 
+        #print "geo_ref",geo_ref
         x_ref = geo_ref.get_xllcorner()
         y_ref = geo_ref.get_yllcorner()
@@ -3325 +3328 @@
         assert allclose(x_ref, 587798.418) # (-38, 148)
         assert allclose(y_ref, 5793123.477)# (-38, 148.5)
-        
-        #Zone:   55    
-        #Easting:  588095.674  Northing: 5821451.722 
+
+        #Zone:   55
+        #Easting:  588095.674  Northing: 5821451.722
         #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
         assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
 
-        #Zone:   55    
-        #Easting:  632145.632  Northing: 5820863.269 
-        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 '' 
+        #Zone:   55
+        #Easting:  632145.632  Northing: 5820863.269
+        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
         assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
 
-        #Zone:   55    
-        #Easting:  609748.788  Northing: 5793447.860 
-        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 '' 
+        #Zone:   55
+        #Easting:  609748.788  Northing: 5793447.860
+        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
         assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
 
@@ -3346 +3349 @@
         #(4000+1000)*60
         assert allclose(xmomentum[1][1],300000.0 )
-        
-        
-        fid.close()
-    
+
+
+        fid.close()
+
         #tidy up
         os.remove(bath_dir_filename)
@@ -3357 +3360 @@
         os.remove(elevation_dir_filename2)
         os.rmdir(elevation_dir)
-        
+
         os.remove(ucur_dir_filename1)
         os.remove(ucur_dir_filename2)
         os.rmdir(ucur_dir)
-        
+
         os.remove(vcur_dir_filename1)
         os.remove(vcur_dir_filename2)
@@ -3369 +3372 @@
         # remove sww file
         os.remove(sww_file)
-        
+
     def test_asc_csiro2sww2(self):
         import tempfile
 
         bath_dir = tempfile.mkdtemp()
-        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
+        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
         #bath_dir = 'bath_data_manager_test'
-        #print "os.getcwd( )",os.getcwd( ) 
+        #print "os.getcwd( )",os.getcwd( )
         elevation_dir =  tempfile.mkdtemp()
         #elevation_dir = 'elev_expanded'
-        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
-        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
-        
+        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
+        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
+
         fid = open(bath_dir_filename, 'w')
         fid.write(""" ncols             3
@@ -3390 +3393 @@
  nodata_value   -9999.0
     9000.000    -1000.000    3000.0
-   -1000.000    9000.000  -1000.000 
+   -1000.000    9000.000  -1000.000
 """)
         fid.close()
-        
+
         fid = open(elevation_dir_filename1, 'w')
         fid.write(""" ncols             3
@@ -3402 +3405 @@
  nodata_value   -9999.0
     9000.000    0.000    3000.0
-     0.000     -9999.000     -9999.000 
+     0.000     -9999.000     -9999.000
 """)
         fid.close()
@@ -3414 +3417 @@
  nodata_value   -9999.0
     9000.000    4000.000    4000.0
-    4000.000    9000.000    4000.000 
+    4000.000    9000.000    4000.000
 """)
         fid.close()
-         
+
         ucur_dir =  tempfile.mkdtemp()
         ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
         ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
-        
+
         fid = open(ucur_dir_filename1, 'w')
         fid.write(""" ncols             3
@@ -3430 +3433 @@
  nodata_value   -9999.0
     90.000    60.000    30.0
-    10.000    10.000    10.000 
+    10.000    10.000    10.000
 """)
         fid.close()
@@ -3441 +3444 @@
  nodata_value   -9999.0
     90.000    60.000    30.0
-    10.000    10.000    10.000 
+    10.000    10.000    10.000
 """)
         fid.close()
-        
+
         vcur_dir =  tempfile.mkdtemp()
         vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
         vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
-        
+
         fid = open(vcur_dir_filename1, 'w')
         fid.write(""" ncols             3
@@ -3457 +3460 @@
  nodata_value   -9999.0
     90.000    60.000    30.0
-    10.000    10.000    10.000 
+    10.000    10.000    10.000
 """)
         fid.close()
@@ -3468 +3471 @@
  nodata_value   -9999.0
     90.000    60.000    30.0
-    10.000    10.000    10.000 
+    10.000    10.000    10.000
 """)
         fid.close()
-        
+
         try:
             asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
@@ -3483 +3486 @@
             os.remove(elevation_dir_filename2)
             os.rmdir(elevation_dir)
-        
+
             os.remove(ucur_dir_filename1)
             os.remove(ucur_dir_filename2)
             os.rmdir(ucur_dir)
-        
+
             os.remove(vcur_dir_filename1)
             os.remove(vcur_dir_filename2)
@@ -3500 +3503 @@
             os.rmdir(elevation_dir)
             raise 'Should raise exception'
-        
+
             os.remove(ucur_dir_filename1)
             os.remove(ucur_dir_filename2)
             os.rmdir(ucur_dir)
-        
+
             os.remove(vcur_dir_filename1)
             os.remove(vcur_dir_filename2)
             os.rmdir(vcur_dir)
 
-        
-     
+
+
     def test_asc_csiro2sww3(self):
         import tempfile
 
         bath_dir = tempfile.mkdtemp()
-        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
+        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
         #bath_dir = 'bath_data_manager_test'
-        #print "os.getcwd( )",os.getcwd( ) 
+        #print "os.getcwd( )",os.getcwd( )
         elevation_dir =  tempfile.mkdtemp()
         #elevation_dir = 'elev_expanded'
-        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
-        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
-        
+        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
+        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
+
         fid = open(bath_dir_filename, 'w')
         fid.write(""" ncols             3
@@ -3531 +3534 @@
  nodata_value   -9999.0
     9000.000    -1000.000    3000.0
-   -1000.000    9000.000  -1000.000 
+   -1000.000    9000.000  -1000.000
 """)
         fid.close()
-        
+
         fid = open(elevation_dir_filename1, 'w')
         fid.write(""" ncols             3
@@ -3543 +3546 @@
  nodata_value   -9999.0
     9000.000    0.000    3000.0
-     0.000     -9999.000     -9999.000 
+     0.000     -9999.000     -9999.000
 """)
         fid.close()
@@ -3555 +3558 @@
  nodata_value   -9999.0
     9000.000    4000.000    4000.0
-    4000.000    9000.000    4000.000 
+    4000.000    9000.000    4000.000
 """)
         fid.close()
-         
+
         ucur_dir =  tempfile.mkdtemp()
         ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
         ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
-        
+
         fid = open(ucur_dir_filename1, 'w')
         fid.write(""" ncols             3
@@ -3571 +3574 @@
  nodata_value   -9999.0
     90.000    60.000    30.0
-    10.000    10.000    10.000 
+    10.000    10.000    10.000
 """)
         fid.close()
@@ -3582 +3585 @@
  nodata_value   -9999.0
     90.000    60.000    30.0
-    10.000    10.000    10.000 
+    10.000    10.000    10.000
 """)
         fid.close()
-        
+
         vcur_dir =  tempfile.mkdtemp()
         vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
         vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
-        
+
         fid = open(vcur_dir_filename1, 'w')
         fid.write(""" ncols             3
@@ -3598 +3601 @@
  nodata_value   -9999.0
     90.000    60.000    30.0
-    10.000    10.000    10.000 
+    10.000    10.000    10.000
 """)
         fid.close()
@@ -3609 +3612 @@
  nodata_value   -9999.0
     90.000    60.000    30.0
-    10.000    10.000    10.000 
+    10.000    10.000    10.000
 """)
         fid.close()
-        
+
         sww_file = 'a_test.sww'
         asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
@@ -3619 +3622 @@
 
         # check the sww file
-        
+
         fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
         x = fid.variables['x'][:]
@@ -3627 +3630 @@
         xmomentum = fid.variables['xmomentum'][:]
         geo_ref = Geo_reference(NetCDFObject=fid)
-        #print "geo_ref",geo_ref 
+        #print "geo_ref",geo_ref
         x_ref = geo_ref.get_xllcorner()
         y_ref = geo_ref.get_yllcorner()
@@ -3633 +3636 @@
         assert allclose(x_ref, 587798.418) # (-38, 148)
         assert allclose(y_ref, 5793123.477)# (-38, 148.5)
-        
-        #Zone:   55    
-        #Easting:  588095.674  Northing: 5821451.722 
+
+        #Zone:   55
+        #Easting:  588095.674  Northing: 5821451.722
         #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
         assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
 
-        #Zone:   55    
-        #Easting:  632145.632  Northing: 5820863.269 
-        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 '' 
+        #Zone:   55
+        #Easting:  632145.632  Northing: 5820863.269
+        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
         assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
 
-        #Zone:   55    
-        #Easting:  609748.788  Northing: 5793447.860 
-        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 '' 
+        #Zone:   55
+        #Easting:  609748.788  Northing: 5793447.860
+        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
         assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
 
@@ -3655 +3658 @@
         #(100.0 - 9000)*10
         assert allclose(xmomentum[0][4], -89000.0 )
-        
+
         #(100.0 - -1000.000)*10
         assert allclose(xmomentum[0][5], 11000.0 )
-        
-        fid.close()
-    
+
+        fid.close()
+
         #tidy up
         os.remove(bath_dir_filename)
@@ -3668 +3671 @@
         os.remove(elevation_dir_filename2)
         os.rmdir(elevation_dir)
-        
+
         os.remove(ucur_dir_filename1)
         os.remove(ucur_dir_filename2)
         os.rmdir(ucur_dir)
-        
+
         os.remove(vcur_dir_filename1)
         os.remove(vcur_dir_filename2)
@@ -3679 +3682 @@
         # remove sww file
         os.remove(sww_file)
-   
-      
+
+
     def test_asc_csiro2sww4(self):
         """
         Test specifying the extent
         """
-        
+
         import tempfile
 
         bath_dir = tempfile.mkdtemp()
-        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
+        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
         #bath_dir = 'bath_data_manager_test'
-        #print "os.getcwd( )",os.getcwd( ) 
+        #print "os.getcwd( )",os.getcwd( )
         elevation_dir =  tempfile.mkdtemp()
         #elevation_dir = 'elev_expanded'
-        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
-        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
-        
+        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
+        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
+
         fid = open(bath_dir_filename, 'w')
         fid.write(""" ncols             4
@@ -3710 +3713 @@
 """)
         fid.close()
-        
+
         fid = open(elevation_dir_filename1, 'w')
         fid.write(""" ncols             4
@@ -3738 +3741 @@
 """)
         fid.close()
-         
+
         ucur_dir =  tempfile.mkdtemp()
         ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
         ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
-        
+
         fid = open(ucur_dir_filename1, 'w')
         fid.write(""" ncols             4
@@ -3769 +3772 @@
 """)
         fid.close()
-        
+
         vcur_dir =  tempfile.mkdtemp()
         vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
         vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
-        
+
         fid = open(vcur_dir_filename1, 'w')
         fid.write(""" ncols             4
@@ -3800 +3803 @@
 """)
         fid.close()
-        
+
         sww_file = tempfile.mktemp(".sww")
         #sww_file = 'a_test.sww'
@@ -3812 +3815 @@
 
         # check the sww file
-        
+
         fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
         x = fid.variables['x'][:]
@@ -3821 +3824 @@
         ymomentum = fid.variables['ymomentum'][:]
         geo_ref = Geo_reference(NetCDFObject=fid)
-        #print "geo_ref",geo_ref 
+        #print "geo_ref",geo_ref
         x_ref = geo_ref.get_xllcorner()
         y_ref = geo_ref.get_yllcorner()
@@ -3830 +3833 @@
         assert allclose(y_ref,  5820863.269 )# (-37.45, 148.5)
 
-        #Easting:  632145.632  Northing: 5820863.269 
-        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 '' 
-        
+        #Easting:  632145.632  Northing: 5820863.269
+        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
+
         #print "x",x
-        #print "y",y 
+        #print "y",y
         self.failUnless(len(x) == 4,'failed') # 2*2
         self.failUnless(len(x) == 4,'failed') # 2*2
-        
+
         #Zone:   55
-        #Easting:  632145.632  Northing: 5820863.269 
+        #Easting:  632145.632  Northing: 5820863.269
         #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
-        # magic number - y is close enough for me. 
+        # magic number - y is close enough for me.
         assert allclose(x[3], 632145.63 - x_ref)
         assert allclose(y[3], 5820863.269  - y_ref + 5.22155314684e-005)
 
         assert allclose(z[0],9000.0 ) #z is elevation info
-        #print "z",z 
+        #print "z",z
         # 2 time steps, 4 points
-        self.failUnless(xmomentum.shape == (2,4), 'failed') 
-        self.failUnless(ymomentum.shape == (2,4), 'failed') 
-        
+        self.failUnless(xmomentum.shape == (2,4), 'failed')
+        self.failUnless(ymomentum.shape == (2,4), 'failed')
+
         #(100.0 - -1000.000)*10
         #assert allclose(xmomentum[0][5], 11000.0 )
-        
-        fid.close()
-        
+
+        fid.close()
+
         # is the sww file readable?
         #Lets see if we can convert it to a dem!
-        #print "sww_file",sww_file 
+        #print "sww_file",sww_file
         #dem_file = tempfile.mktemp(".dem")
         domain = sww2domain(sww_file) ###, dem_file)
@@ -3870 +3873 @@
         os.remove(elevation_dir_filename2)
         os.rmdir(elevation_dir)
-        
+
         os.remove(ucur_dir_filename1)
         os.remove(ucur_dir_filename2)
         os.rmdir(ucur_dir)
-        
+
         os.remove(vcur_dir_filename1)
         os.remove(vcur_dir_filename2)
@@ -3880 +3883 @@
 
 
-      
-        
+
+
         # remove sww file
         os.remove(sww_file)
-        
+
         # remove dem file
         #os.remove(dem_file)
@@ -3907 +3910 @@
                         longitudes == longitudes_news,
                          'failed')
-        
+
         ## 2nd test
         kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
@@ -3918 +3921 @@
         #print "latitudes_new", latitudes_new
         #print "longitudes_news",longitudes_news
-        
+
         self.failUnless(latitudes == latitudes_new and \
                         longitudes == longitudes_news,
                          'failed')
-        
+
         ## 3rd test
         kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(latitudes,
@@ -3933 +3936 @@
         #print "latitudes_new", latitudes_new
         #print "longitudes_news",longitudes_news
-        
+
         self.failUnless(latitudes_new == [2, 1] and \
                         longitudes_news == [10, 20],
                          'failed')
 
-        
+
         ## 4th test
         kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
@@ -3949 +3952 @@
         #print "latitudes_new", latitudes_new
         #print "longitudes_news",longitudes_news
-        
+
         self.failUnless(latitudes_new == [2, 1, 0] and \
                         longitudes_news == [0, 10, 20],
@@ -3963 +3966 @@
         #print "latitudes_new", latitudes_new
         #print "longitudes_news",longitudes_news
-        
+
         self.failUnless(latitudes_new == [2, 1, 0] and \
                         longitudes_news == [0, 10, 20],
                          'failed')
-        
+
         ## 6th test
-        
+
         kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
             latitudes,longitudes,
@@ -3979 +3982 @@
         #print "latitudes_new", latitudes_new
         #print "longitudes_news",longitudes_news
-        
+
         self.failUnless(latitudes_new == [3, 2, 1] and \
                         longitudes_news == [10, 20, 30],
                          'failed')
-        
-        
+
+
         ## 7th test
         m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
@@ -3995 +3998 @@
         longitudes_news = longitudes[lmin:lmax]
         m2d = m2d[kmin:kmax,lmin:lmax]
-        #print "m2d", m2d 
+        #print "m2d", m2d
         #print "latitudes_new", latitudes_new
         #print "longitudes_news",longitudes_news
-        
+
         self.failUnless(latitudes_new == [2, 1] and \
                         longitudes_news == [10, 20],
                          'failed')
-        
+
         self.failUnless(m2d == [[5,6],[9,10]],
                          'failed')
-        
+
     def test_get_min_max_indexes2(self):
         latitudes = [-30,-35,-40,-45]
@@ -4011 +4014 @@
 
         m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
-        
+
         # k - lat
         # l - lon
@@ -4020 +4023 @@
         #print "kmin",kmin;print "kmax",kmax
         #print "lmin",lmin;print "lmax",lmax
-        #print "m2d", m2d 
+        #print "m2d", m2d
         #print "latitudes", latitudes
         #print "longitudes",longitudes
-        #print "latitudes[kmax]", latitudes[kmax] 
+        #print "latitudes[kmax]", latitudes[kmax]
         latitudes_new = latitudes[kmin:kmax]
         longitudes_new = longitudes[lmin:lmax]
@@ -4030 +4033 @@
         #print "latitudes_new", latitudes_new
         #print "longitudes_new",longitudes_new
-        
+
         self.failUnless(latitudes_new == [-30, -35, -40] and \
                         longitudes_new == [148, 149,150],
@@ -4036 +4039 @@
         self.failUnless(m2d == [[0,1,2],[4,5,6],[8,9,10]],
                          'failed')
-        
+
     def test_get_min_max_indexes3(self):
         latitudes = [-30,-35,-40,-45,-50,-55,-60]
@@ -4047 +4050 @@
             -43,-37,148.5,149.5)
 
-        
+
         #print "kmin",kmin;print "kmax",kmax
         #print "lmin",lmin;print "lmax",lmax
@@ -4056 +4059 @@
         #print "latitudes_new", latitudes_new
         #print "longitudes_news",longitudes_news
-        
+
         self.failUnless(latitudes_new == [-35, -40, -45] and \
                         longitudes_news == [148, 149,150],
                          'failed')
-        
+
     def test_get_min_max_indexes4(self):
         latitudes = [-30,-35,-40,-45,-50,-55,-60]
@@ -4069 +4072 @@
         kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
             latitudes,longitudes)
-        
+
         #print "kmin",kmin;print "kmax",kmax
         #print "lmin",lmin;print "lmax",lmax
@@ -4078 +4081 @@
         #print "latitudes_new", latitudes_new
         #print "longitudes_news",longitudes_news
-        
+
         self.failUnless(latitudes_new == latitudes  and \
                         longitudes_news == longitudes,
                          'failed')
-        
+
     def test_tsh2sww(self):
         import os
         import tempfile
-        
+
         tsh_file = tempfile.mktemp(".tsh")
         file = open(tsh_file,"w")
@@ -4124 +4127 @@
 
         #sww_file = tempfile.mktemp(".sww")
-        #print "sww_file",sww_file 
+        #print "sww_file",sww_file
         #print "sww_file",tsh_file
         tsh2sww(tsh_file)
@@ -4135 +4138 @@
     suite = unittest.makeSuite(Test_Data_Manager,'test')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_missing_points')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_elevation')    
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_elevation')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
@@ -4141 +4144 @@
     runner = unittest.TextTestRunner()
     runner.run(suite)
- 
+
 
 
@@ -4378 +4381 @@
                 if i == 3 and j == 3: z = NODATA_value
 
-                
+
                 if z <> NODATA_value:
                     new_ref_elev.append(z)
                     new_ref_pts.append( [x,y] )
-                
+
                 ref_points.append( [x,y] )
                 ref_elevation.append(z)
@@ -4391 +4394 @@
         fid.close()
 
-       
+
         #Write prj file with metadata
         metafilename = root+'.prj'
@@ -4482 +4485 @@
                     new_ref_elev1.append(z)
                     new_ref_pts1.append( [x,y] )
-                    
+
                 ref_points.append( [x,y] )
                 ref_elevation.append(z)
@@ -4506 +4509 @@
 """)
         fid.close()
-       
+
         #Convert to NetCDF pts
         convert_dem_from_ascii2netcdf(root)
@@ -4549 +4552 @@
                     new_ref_pts2.append( [x_new,y_new] )
 
-                    
+
                 ref_points.append( [x_new,y_new] )
                 ref_elevation.append(z)
@@ -4561 +4564 @@
         #assert allclose(elevation, ref_elevation)
 
-        
+
         assert len(points) == len(new_ref_pts2), 'length of returned points not correct'
         assert allclose(points, new_ref_pts2), 'points do not align'
@@ -4580 +4583 @@
 
 
-

inundation/zeus/Validation.zpi

@@ -73 +73 @@
     <ReleaseProjectSaveAll>Off</ReleaseProjectSaveAll>
     <ReleaseProjectReload>Off</ReleaseProjectReload>
-    <file>..\validation\completed\lwru2\create_mesh.py</file>
-    <file>..\validation\completed\lwru2\extract_timeseries.py</file>
-    <file>..\validation\completed\lwru2\lwru2.py</file>
-    <file>..\validation\completed\lwru2\project.py</file>
     <folder name="Header Files" />
     <folder name="Resource Files" />

inundation/zeus/analytic_solutions.zpi

@@ -73 +73 @@
     <ReleaseProjectSaveAll>Off</ReleaseProjectSaveAll>
     <ReleaseProjectReload>Off</ReleaseProjectReload>
-    <file>..\analytical solutions\Analytical solution_oblique_shock.py</file>
-    <file>..\analytical solutions\Analytical solution_oblique_shock_01.py</file>
-    <file>..\analytical solutions\Analytical_solution_circular_hydraulic_jump.py</file>
-    <file>..\analytical solutions\Non_symmetrical_dam_break.py</file>
+    <file>..\..\development\analytical solutions\Analytical_solution_circular_hydraulic_jump.py</file>
+    <file>..\..\development\analytical solutions\Analytical_solution_Sampson_cross_mesh.py</file>
+    <file>..\..\development\analytical solutions\Analytical_solution_Yoon_cross_mesh.py</file>
 </project>

inundation/zeus/anuga-workspace.ini

@@ -1 +1 @@
 [Files MRU]
-0=c:\home\projects\anuga\inundation\pyvolution\shallow_water.py
-1=c:\home\projects\anuga\inundation\pyvolution\general_mesh.py
-2=c:\apps\lyx\bin\reLyX.bat
-3=c:\program files\zfw\readme.txt
+0=c:\home\projects\anuga\production\merimbula_2005\run_merimbula_lake.py
+1=c:\home\projects\anuga\production\merimbula_2005\run_new_meribula.py
+2=c:\home\projects\anuga\production\merimbula_2005\prepare.py
+3=c:\home\projects\anuga\inundation\pyvolution\shallow_water.py
+4=c:\home\projects\anuga\inundation\pyvolution\general_mesh.py
+5=c:\apps\lyx\bin\reLyX.bat
+6=c:\program files\zfw\readme.txt

inundation/zeus/anuga-workspace.zwi

@@ -2 +2 @@
 <workspace name="anuga-workspace">
     <mode>Debug</mode>
-    <active>Merimbula</active>
+    <active>analytic_solutions</active>
     <project name="analytic_solutions">analytic_solutions.zpi</project>
     <project name="euler">euler.zpi</project>

inundation/zeus/parallel.zpi

@@ -84 +84 @@
     <file>..\parallel\run_parallel_advection.py</file>
     <file>..\parallel\run_parallel_merimbula.py</file>
-    <file>..\parallel\run_parallel_mesh.py</file>
     <file>..\parallel\run_parallel_sw_merimbula.py</file>
     <file>..\parallel\run_parallel_sw_merimbula_metis.py</file>

production/merimbula_2005/project.py

@@ -7 +7 @@
 bathymetry_filename = 'merimbula_bathymetry.xya'
 mesh_filename = 'merimbula_10785.tsh'
+sww_file = 'merimbula_weed.sww'

production/merimbula_2005/run_new_meribula_weed.py

@@ -156 +156 @@
 # close to bearing = 0 and 360!
 #---------------------------------------------
-filename = project.original_wind_filename[:-4]+'.tms'
-print 'Wind field from ',filename
-wind = file_function(filename, domain=domain, quantities=['speed','bearing'])
-domain.forcing_terms.append(Wind_stress(wind))
+#filename = project.original_wind_filename[:-4]+'.tms'
+#print 'Wind field from ',filename
+#wind = file_function(filename, domain=domain, quantities=['speed','bearing'])
+#domain.forcing_terms.append(Wind_stress(wind))