Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 4712)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 4713)
@@ -208,4 +208,6 @@
         self.min_timestep = self.max_timestep = 0.0
         self.starttime = 0 #Physical starttime if any (0 is 1 Jan 1970 00:00:00)
+        self.timestep = 0.0
+        self.flux_timestep = 0.0
 
         # Monitoring
@@ -242,13 +244,5 @@
 
 
-    def set_default_order(self, n):
-        """Set default (spatial) order to either 1 or 2
-        """
-
-        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
+
 
 
@@ -295,4 +289,15 @@
 
         return self.time
+
+    def set_default_order(self, n):
+        """Set default (spatial) order to either 1 or 2
+        """
+
+        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
+        
 
     def set_quantity_vertices_dict(self, quantity_dict):
@@ -964,5 +969,5 @@
     def set_timestepping_method(self,timestepping_method):
         
-        if timestepping_method in ['euler', 'rk2']:
+        if timestepping_method in ['euler', 'rk2', 'rk3']:
             self.timestepping_method = timestepping_method
             return
@@ -1113,6 +1118,9 @@
                 
             elif self.get_timestepping_method() == 'rk2':
-                self.evolve_one_rk2_step(yieldstep,finaltime)               
-
+                self.evolve_one_rk2_step(yieldstep,finaltime)
+
+            elif self.get_timestepping_method() == 'rk3':
+                self.evolve_one_rk3_step(yieldstep,finaltime)               
+            
             # Update extrema if necessary (for reporting)
             self.update_extrema()            
@@ -1183,5 +1191,5 @@
 
 
-    def evolve_one_rk2_step(self,yieldstep, finaltime):
+    def evolve_one_rk2_step(self, yieldstep, finaltime):
         """One 2nd order RK timestep"""
 
@@ -1226,13 +1234,72 @@
         #------------------------------------
         #Combine initial and final values
-        #of conserved quantities
+        #of conserved quantities and cleanup
         #------------------------------------
-
+        #combine steps
         self.saxpy_conserved_quantities(0.5, 0.5)
+ 
+        #update ghosts
+        self.update_ghosts()
+
+        #Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        #Update boundary values
+        self.update_boundary()
+
+
+
+    def evolve_one_rk3_step(self, yieldstep, finaltime):
+        """One 2nd order RK timestep"""
+
+
+
+        #Save initial initial conserved quantities values
+        self.backup_conserved_quantities()            
+
+        initial_time = self.time
+        
+        #--------------------------------------
+        #First euler step
+        #--------------------------------------
+
+        #Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        #Update timestep to fit yieldstep and finaltime
+        self.update_timestep(yieldstep, finaltime)
+
+        #Update conserved quantities
+        self.update_conserved_quantities()
+
+        #update ghosts
+        self.update_ghosts()
+
+        #Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        #Update boundary values
+        self.update_boundary()
+
+        #Update time
+        self.time += self.timestep
 
         #------------------------------------
-        #Clean up rk step
+        #Second Euler step
         #------------------------------------
             
+        #Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        #Update conserved quantities
+        self.update_conserved_quantities()
+
+        #------------------------------------
+        #Combine final and initial values
+        #of conserved quantities and cleanup
+        #------------------------------------
+        #combine steps
+        self.saxpy_conserved_quantities(0.25, 0.75)
+ 
         #update ghosts
         self.update_ghosts()
@@ -1243,4 +1310,38 @@
         #Update boundary values
         self.update_boundary()
+
+        #set substep time
+        self.time = initial_time + self.timestep*0.5
+
+        #------------------------------------
+        #Third Euler step
+        #------------------------------------
+            
+        #Compute fluxes across each element edge
+        self.compute_fluxes()
+
+        #Update conserved quantities
+        self.update_conserved_quantities()
+
+        #------------------------------------
+        #Combine final and initial values
+        #of conserved quantities and cleanup
+        #------------------------------------
+        #combine steps
+        self.saxpy_conserved_quantities(2.0/3.0, 1.0/3.0)
+ 
+        #update ghosts
+        self.update_ghosts()
+
+        #Update vertex and edge values
+        self.distribute_to_vertices_and_edges()
+
+        #Update boundary values
+        self.update_boundary()
+
+        #set substep time
+        self.time = initial_time + self.timestep       
+        
+
 
 
@@ -1344,5 +1445,5 @@
         # self.timestep is calculated from speed of characteristics
         # Apply CFL condition here
-        timestep = min(self.CFL*self.timestep, max_timestep)
+        timestep = min(self.CFL*self.flux_timestep, max_timestep)
 
         #Record maximal and minimal values of timestep for reporting
Index: /anuga_core/source/anuga/shallow_water/shallow_water_domain.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 4712)
+++ /anuga_core/source/anuga/shallow_water/shallow_water_domain.py	(revision 4713)
@@ -824,5 +824,5 @@
 
 
-    domain.timestep = timestep
+    domain.flux_timestep = timestep
 
 #MH090605 The following method belongs to the shallow_water domain class
@@ -875,5 +875,5 @@
          compute_fluxes_ext_central as compute_fluxes_ext
 
-    domain.timestep = compute_fluxes_ext(timestep, domain.epsilon,
+    domain.flux_timestep = compute_fluxes_ext(timestep, domain.epsilon,
                                          domain.H0,
                                          domain.g,
Index: /anuga_work/development/demos/advection_dirichlet.py
===================================================================
--- /anuga_work/development/demos/advection_dirichlet.py	(revision 4712)
+++ /anuga_work/development/demos/advection_dirichlet.py	(revision 4713)
Index: /anuga_work/development/demos/advection_waves.py
===================================================================
--- /anuga_work/development/demos/advection_waves.py	(revision 4712)
+++ /anuga_work/development/demos/advection_waves.py	(revision 4713)
Index: /anuga_work/development/demos/bed_w_eden_boundary.py
===================================================================
--- /anuga_work/development/demos/bed_w_eden_boundary.py	(revision 4712)
+++ /anuga_work/development/demos/bed_w_eden_boundary.py	(revision 4713)
Index: /anuga_work/development/demos/bed_w_file_boundary.py
===================================================================
--- /anuga_work/development/demos/bed_w_file_boundary.py	(revision 4712)
+++ /anuga_work/development/demos/bed_w_file_boundary.py	(revision 4713)
Index: /anuga_work/development/demos/bedslope.py
===================================================================
--- /anuga_work/development/demos/bedslope.py	(revision 4712)
+++ /anuga_work/development/demos/bedslope.py	(revision 4713)
Index: /anuga_work/development/demos/bedslope_bsod_test.py
===================================================================
--- /anuga_work/development/demos/bedslope_bsod_test.py	(revision 4712)
+++ /anuga_work/development/demos/bedslope_bsod_test.py	(revision 4713)
Index: /anuga_work/development/demos/netherlands.py
===================================================================
--- /anuga_work/development/demos/netherlands.py	(revision 4712)
+++ /anuga_work/development/demos/netherlands.py	(revision 4713)
@@ -83,5 +83,5 @@
 N = 130 #size = 33800
 N = 600 #Size = 720000
-N = 50
+N = 100
 
 
@@ -97,14 +97,17 @@
 domain.check_integrity()
 
+domain.set_timestepping_method('rk3')
+domain.set_default_order(2)
+
 #Setup order and all the beta's for the limiters (these should become defaults
-domain.default_order = 2
-domain.beta_w      = 1.0
-domain.beta_w_dry  = 0.2
-domain.beta_uh     = 1.0
-domain.beta_uh_dry = 0.2
-domain.beta_vh     = 1.0
-domain.beta_vh_dry = 0.2
 
-domain.alpha_balance = 100.0
+## domain.beta_w      = 1.0
+## domain.beta_w_dry  = 0.2
+## domain.beta_uh     = 1.0
+## domain.beta_uh_dry = 0.2
+## domain.beta_vh     = 1.0
+## domain.beta_vh_dry = 0.2
+
+## domain.alpha_balance = 100.0
 
 #Output params
@@ -117,5 +120,5 @@
 #Set bed-slope and friction
 inflow_stage = 0.5
-manning = 0.03
+manning = 0.0
 Z = Weir(inflow_stage)
 
@@ -186,3 +189,3 @@
 
 
-v.join()
+vis.join()
Index: /anuga_work/development/demos/pressure_force.py
===================================================================
--- /anuga_work/development/demos/pressure_force.py	(revision 4712)
+++ /anuga_work/development/demos/pressure_force.py	(revision 4713)
Index: /anuga_work/development/demos/run_tsh.py
===================================================================
--- /anuga_work/development/demos/run_tsh.py	(revision 4712)
+++ /anuga_work/development/demos/run_tsh.py	(revision 4713)
Index: /anuga_work/development/demos/run_tsh_weir.py
===================================================================
--- /anuga_work/development/demos/run_tsh_weir.py	(revision 4712)
+++ /anuga_work/development/demos/run_tsh_weir.py	(revision 4713)
