Index: /inundation/ga/storm_surge/pyvolution/data_manager.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 282)
+++ /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 283)
@@ -138,4 +138,5 @@
 #Native checkpoint format.
 #Information needed to recreate a state is preserved
+#FIXME: Rethink and maybe use netcdf format
 def cpt_variable_writer(filename, t, v0, v1, v2):
     """Store all conserved quantities to file
@@ -218,4 +219,6 @@
 
 #Function for storing out to e.g. visualisation
+#FIXME: Do we want this?
+#FIXME: Not done yet for this version
 def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
     """Writes x,y,z coordinates of triangles constituting the bed elevation.
@@ -456,4 +459,5 @@
 
 #Function for storing out to e.g. visualisation
+#FIXME Not done yet for this version
 class Data_format_xya(Data_format):
     """Generic interface to data formats
Index: /inundation/ga/storm_surge/pyvolution/netherlands.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/netherlands.py	(revision 283)
+++ /inundation/ga/storm_surge/pyvolution/netherlands.py	(revision 283)
@@ -0,0 +1,175 @@
+"""Example of shallow water wave equation.
+
+This is called Netherlands because it shows a dam with a gap in it and
+stylised housed behind it and below the water surface.
+
+"""
+
+######################
+# Module imports 
+#
+from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
+     Transmissive_boundary, Time_boundary, Constant_height
+
+from mesh_factory import rectangular
+from Numeric import array
+    
+
+class Weir:
+    """Set a bathymetry for simple weir with a hole.
+    x,y are assumed to be in the unit square
+    """
+    
+    def __init__(self, stage):
+        self.inflow_stage = stage
+
+    def __call__(self, x, y):    
+        from Numeric import zeros, Float 
+    
+        N = len(x)
+        assert N == len(y)    
+
+        z = zeros(N, Float)
+        for i in range(N):
+            z[i] = -x[i]/20  #General slope
+
+            #Flattish bit to the left
+            if x[i] <= 0.3:
+                #z[i] = -x[i]/5
+                z[i] = -x[i]/20
+                
+            
+            #Weir
+            if x[i] > 0.3 and x[i] < 0.4:
+                z[i] = -x[i]/20+1.2
+
+            #Dip
+            #if x[i] > 0.6 and x[i] < 0.9:
+            #    z[i] = -x[i]/20-0.5  #-y[i]/5
+
+            #Hole in weir
+            #if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
+            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.4 and y[i] < 0.6:
+                #z[i] = -x[i]/5
+                z[i] = -x[i]/20
+            
+            #Poles
+            #if x[i] > 0.65 and x[i] < 0.8 and y[i] > 0.55 and y[i] < 0.65 or\
+            #       x[i] > 0.75 and x[i] < 0.9 and y[i] > 0.35 and y[i] < 0.45:
+            #    z[i] = -x[i]/20+0.4
+
+            if (x[i] - 0.72)**2 + (y[i] - 0.6)**2 < 0.05**2:# or\
+                   #x[i] > 0.75 and x[i] < 0.9 and y[i] > 0.35 and y[i] < 0.45:
+                z[i] = -x[i]/20+0.4
+
+            #Wall
+            if x[i] > 0.995:
+                z[i] = -x[i]/20+0.3                
+
+        return z/2
+        
+    
+    
+######################
+# Domain
+#
+
+N = 250
+#N= 8
+N = 16
+#N = 4
+#N = 102
+N = 25
+N = 16
+N = 60
+N = 150 #size = 45000
+N = 130 #size = 33800
+#N = 60
+#N = 40
+N = 260
+#N = 150
+N = 264
+
+N = 500
+
+print 'Creating domain'
+#Create basic mesh
+points, vertices, boundary = rectangular(N, N)
+
+#Create shallow water domain
+domain = Domain(points, vertices, boundary)
+    
+domain.check_integrity()
+domain.default_order = 2
+domain.smooth = True
+domain.reduction = min  #Looks a lot better on top of steep slopes
+
+print "Number of triangles = ", len(domain)
+
+
+if N > 50:
+    domain.visualise = False
+    domain.checkpoint = False
+    domain.store = True    #Store for visualisation purposes
+    domain.format = 'sww'   #Native netcdf visualisation format
+    import sys, os
+    root, ext = os.path.splitext(sys.argv[0])
+    if domain.smooth is True:
+        s = 'smooth'
+    else:
+        s = 'nonsmooth'        
+    domain.filename = root + '_' + s
+else:
+    domain.visualise = True
+    domain.checkpoint = False
+    domain.store = False    
+
+    
+#Set bed-slope and friction
+inflow_stage = 0.08
+manning = 0.0
+Z = Weir(inflow_stage)
+
+print 'Field values'
+domain.set_quantity('elevation', Z)
+domain.set_quantity('friction', manning)
+
+
+######################
+# Boundary conditions
+#
+print 'Boundaries'
+Br = Reflective_boundary(domain)
+Bt = Transmissive_boundary(domain)
+
+#Constant inflow
+Bd = Dirichlet_boundary([2*inflow_stage, 0.0, 0.0])
+
+
+#Time dependent inflow
+from math import sin, pi, exp
+Bw = Time_boundary(domain=domain,
+                   f=lambda t: [0.7*(1 + sin(t*pi/2.5))*\
+                                (inflow_stage*(sin(t*pi)+0.7)),0,0])
+
+#Set boundary conditions
+#domain.set_boundary({'left': Bw, 'right': Br, 'bottom': Br, 'top': Br})
+domain.set_boundary({'left': Bd, 'right': Bt, 'bottom': Bt, 'top': Bt})
+
+                  
+
+######################
+#Initial condition
+#
+print 'Initial condition'
+domain.set_quantity('level', Constant_height(Z, 0.))
+
+#Evolve
+import time
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 0.1, finaltime = 1.5):
+    domain.write_time()
+    
+print 'That took %.2f seconds' %(time.time()-t0)    
+    
Index: /inundation/ga/storm_surge/pyvolution/quantity.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/quantity.py	(revision 282)
+++ /inundation/ga/storm_surge/pyvolution/quantity.py	(revision 283)
@@ -208,5 +208,5 @@
         TODO: be able to smooth individual fields
         NOTE:  This function does not have a test.
-        FIXME: NOT DONE
+        FIXME: NOT DONE - do we need it?
         """
 
