Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 6050)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py	(revision 6051)
@@ -122,8 +122,8 @@
 
 
-        #Build dictionary of Quantity instances keyed by quantity names
+        # Build dictionary of Quantity instances keyed by quantity names
         self.quantities = {}
 
-        #FIXME: remove later - maybe OK, though....
+        # FIXME: remove later - maybe OK, though....
         for name in self.conserved_quantities:
             self.quantities[name] = Quantity(self)
@@ -131,8 +131,8 @@
             self.quantities[name] = Quantity(self)
 
-        #Create an empty list for explicit forcing terms
+        # Create an empty list for explicit forcing terms
         self.forcing_terms = []
 
-        #Setup the ghost cell communication
+        # Setup the ghost cell communication
         if full_send_dict is None:
             self.full_send_dict = {}
@@ -1219,5 +1219,4 @@
         # Update ghosts
         self.update_ghosts()
-
 
         # Update time
@@ -1560,5 +1559,4 @@
 
 
-
     def update_conserved_quantities(self):
         """Update vectors of conserved quantities using previously
@@ -1582,5 +1580,6 @@
 
             # Note that Q.explicit_update is reset by compute_fluxes
-            # Where is Q.semi_implicit_update reset?
+            # Where is Q.semi_implicit_update reset? 
+            # It is reset in quantity_ext.c
             
 
Index: /anuga_core/source/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py
===================================================================
--- /anuga_core/source/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py	(revision 6050)
+++ /anuga_core/source/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py	(revision 6051)
@@ -138,4 +138,5 @@
                        culvert_routine=boyd_generalised_culvert_model,        
                        number_of_barrels=1,
+                       update_interval=2,
                        verbose=True)
 
@@ -158,7 +159,6 @@
 #------------------------------------------------------------------------------
 
-for t in domain.evolve(yieldstep = 0.01, finaltime = 45):
-     if int(domain.time*100) % 100 == 0:
-	     domain.write_time()
+#for t in domain.evolve(yieldstep = 1, finaltime = 25):
+#    print domain.timestepping_statistics()
     
 
@@ -170,5 +170,5 @@
 t0 = time.time()
     
-s = 'for t in domain.evolve(yieldstep = 0.5, finaltime = 2): domain.write_time()'
+s = 'for t in domain.evolve(yieldstep = 1, finaltime = 25): domain.write_time()'
 
 import profile, pstats
Index: /anuga_core/source/anuga/culvert_flows/culvert_class.py
===================================================================
--- /anuga_core/source/anuga/culvert_flows/culvert_class.py	(revision 6050)
+++ /anuga_core/source/anuga/culvert_flows/culvert_class.py	(revision 6051)
@@ -64,4 +64,5 @@
                  culvert_routine=None,
                  number_of_barrels=1,
+                 update_interval=None,
                  verbose=False):
         
@@ -191,4 +192,6 @@
         self.verbose = verbose
         self.last_time = 0.0        
+        self.last_update = 0.0 # For use with update_interval        
+        self.update_interval = update_interval
         
 
@@ -232,4 +235,5 @@
         log_to_file(self.log_filename, s)        
         
+        
     def __call__(self, domain):
         from anuga.utilities.numerical_tools import mean
@@ -244,173 +248,194 @@
         # Time stuff
         time = domain.get_time()
-        delta_t = time-self.last_time
-        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
-        log_to_file(log_filename, s)
-        
-        msg = 'Time did not advance'
-        if time > 0.0: assert delta_t > 0.0, msg
-
-
-        # Get average water depths at each opening        
-        openings = self.openings   # There are two Opening [0] and [1]
-        for i, opening in enumerate(openings):
-            dq = domain.quantities
+        
+        
+        update = False
+        if self.update_interval is None:
+            update = True
+        else:    
+            if time - self.last_update > self.update_interval or time == 0.0:
+                update = True
+
+        #print 'call', time, time - self.last_update
+                
+                                
+        if update is True:
+            #print 'Updating', time, time - self.last_update
+            self.last_update = time
+        
+            delta_t = time-self.last_time
+            s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
+            log_to_file(log_filename, s)
+        
+            msg = 'Time did not advance'
+            if time > 0.0: assert delta_t > 0.0, msg
+
+
+            # Get average water depths at each opening        
+            openings = self.openings   # There are two Opening [0] and [1]
+            for i, opening in enumerate(openings):
+                dq = domain.quantities
+                
+                # Compute mean values of selected quantitites in the 
+                # exchange area in front of the culvert
+                # Stage and velocity comes from enquiry area 
+                # and elevation from exchange area
+                
+                stage = dq['stage'].get_values(location='centroids',
+                                               indices=opening.exchange_indices)            
+                w = mean(stage) # Average stage
+
+                # Use invert level instead of elevation if specified
+                invert_level = self.invert_levels[i]
+                if invert_level is not None:
+                    z = invert_level
+                else:
+                    elevation = dq['elevation'].get_values(location='centroids', 
+                                                           indices=opening.exchange_indices)
+                    z = mean(elevation) # Average elevation
+
+                # Estimated depth above the culvert inlet
+                d = w - z  # Used for calculations involving depth
+                if d < 0.0:
+                    # This is possible since w and z are taken at different locations
+                    #msg = 'D < 0.0: %f' %d
+                    #raise msg
+                    d = 0.0
+                
+
+                # Ratio of depth to culvert height.
+                # If ratio > 1 then culvert is running full
+                if self.culvert_type == 'circle':
+                    ratio = d/self.diameter
+                else:    
+                    ratio = d/self.height  
+                opening.ratio = ratio
+                    
+                    
+                # Average measures of energy in front of this opening
+                polyline = self.enquiry_polylines[i]
+                #print 't = %.4f, opening=%d,' %(domain.time, i),
+                opening.total_energy = domain.get_energy_through_cross_section(polyline,
+                                                                               kind='total')            
+                #print 'Et = %.3f m' %opening.total_energy
+
+                # Store current average stage and depth with each opening object
+                opening.depth = d
+                opening.depth_trigger = d # Use this for now
+                opening.stage = w
+                opening.elevation = z
+                
+
+            #################  End of the FOR loop ################################################
+
+            # We now need to deal with each opening individually
+                
+            # Determine flow direction based on total energy difference
+            delta_Et = openings[0].total_energy - openings[1].total_energy
+
+            if delta_Et > 0:
+                #print 'Flow U/S ---> D/S'
+                inlet = openings[0]
+                outlet = openings[1]
+
+                inlet.momentum = self.opening_momentum[0]
+                outlet.momentum = self.opening_momentum[1]
+
+            else:
+                #print 'Flow D/S ---> U/S'
+                inlet = openings[1]
+                outlet = openings[0]
+
+                inlet.momentum = self.opening_momentum[1]
+                outlet.momentum = self.opening_momentum[0]
+                
+                delta_Et = -delta_Et
+
+            self.inlet = inlet
+            self.outlet = outlet
+                
+            msg = 'Total energy difference is negative'
+            assert delta_Et > 0.0, msg
+
+            delta_h = inlet.stage - outlet.stage
+            delta_z = inlet.elevation - outlet.elevation
+            culvert_slope = (delta_z/self.length)
+
+            if culvert_slope < 0.0:
+                # Adverse gradient - flow is running uphill
+                # Flow will be purely controlled by uphill outlet face
+                if self.verbose is True:
+                    print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
+
+
+            s = 'Delta total energy = %.3f' %(delta_Et)
+            log_to_file(log_filename, s)
+
+
+            # Calculate discharge for one barrel and set inlet.rate and outlet.rate
+            Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
+        
+            # Adjust discharge for multiple barrels
+            Q *= self.number_of_barrels
+
+            # Compute barrel momentum
+            barrel_momentum = barrel_velocity*culvert_outlet_depth
+                    
+            s = 'Barrel velocity = %f' %barrel_velocity
+            log_to_file(log_filename, s)
+
+            # Compute momentum vector at outlet
+            outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
+                
+            s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
+            log_to_file(log_filename, s)
+
+            # Log timeseries to file
+            fid = open(self.timeseries_filename, 'a')        
+            fid.write('%f, %f, %f, %f, %f\n'\
+                          %(time, 
+                            openings[0].total_energy,
+                            openings[1].total_energy,
+                            barrel_velocity,
+                            Q))
+            fid.close()
+
+            # Update momentum        
+            delta_t = time - self.last_time
+            if delta_t > 0.0:
+                xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
+                xmomentum_rate /= delta_t
+                    
+                ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
+                ymomentum_rate /= delta_t
+                        
+                s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
+                log_to_file(log_filename, s)                    
+            else:
+                xmomentum_rate = ymomentum_rate = 0.0
+
+
+            # Set momentum rates for outlet jet
+            outlet.momentum[0].rate = xmomentum_rate
+            outlet.momentum[1].rate = ymomentum_rate
+
+            # Remember this value for next step (IMPORTANT)
+            outlet.momentum[0].value = outlet_mom_x
+            outlet.momentum[1].value = outlet_mom_y                    
+
+            if int(domain.time*100) % 100 == 0:
+                s = 'T=%.5f, Culvert Discharge = %.3f f'\
+                    %(time, Q)
+                s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
+                     %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
+                s += ' Momentum rate: (%.4f, %.4f)'\
+                     %(xmomentum_rate, ymomentum_rate)                    
+                s+='Outlet Vel= %.3f'\
+                    %(barrel_velocity)
+                log_to_file(log_filename, s)
             
-            # Compute mean values of selected quantitites in the exchange area in front of the culvert
-            # Stage and velocity comes from enquiry area and elevation from exchange area
-            
-            stage = dq['stage'].get_values(location='centroids',
-                                           indices=opening.exchange_indices)            
-            w = mean(stage) # Average stage
-
-            # Use invert level instead of elevation if specified
-            invert_level = self.invert_levels[i]
-            if invert_level is not None:
-                z = invert_level
-            else:
-                elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices)
-                z = mean(elevation) # Average elevation
-
-            # Estimated depth above the culvert inlet
-            d = w - z  # Used for calculations involving depth
-            if d < 0.0:
-                # This is possible since w and z are taken at different locations
-                #msg = 'D < 0.0: %f' %d
-                #raise msg
-                d = 0.0
-            
-
-            # Ratio of depth to culvert height.
-            # If ratio > 1 then culvert is running full
-            if self.culvert_type == 'circle':
-                ratio = d/self.diameter
-            else:    
-                ratio = d/self.height  
-            opening.ratio = ratio
-                
-                
-            # Average measures of energy in front of this opening
-            polyline = self.enquiry_polylines[i]
-            #print 't = %.4f, opening=%d,' %(domain.time, i),
-            opening.total_energy = domain.get_energy_through_cross_section(polyline,
-                                                                           kind='total')            
-            #print 'Et = %.3f m' %opening.total_energy
-
-            # Store current average stage and depth with each opening object
-            opening.depth = d
-            opening.depth_trigger = d # Use this for now
-            opening.stage = w
-            opening.elevation = z
-            
-
-        #################  End of the FOR loop ################################################
-
-            
-        # We now need to deal with each opening individually
-
-        # Determine flow direction based on total energy difference
-        delta_Et = openings[0].total_energy - openings[1].total_energy
-
-        if delta_Et > 0:
-            #print 'Flow U/S ---> D/S'
-            inlet = openings[0]
-            outlet = openings[1]
-
-            inlet.momentum = self.opening_momentum[0]
-            outlet.momentum = self.opening_momentum[1]
-
-        else:
-            #print 'Flow D/S ---> U/S'
-            inlet = openings[1]
-            outlet = openings[0]
-
-            inlet.momentum = self.opening_momentum[1]
-            outlet.momentum = self.opening_momentum[0]
-            
-            delta_Et = -delta_Et
-
-        msg = 'Total energy difference is negative'
-        assert delta_Et > 0.0, msg
-
-        delta_h = inlet.stage - outlet.stage
-        delta_z = inlet.elevation - outlet.elevation
-        culvert_slope = (delta_z/self.length)
-
-        if culvert_slope < 0.0:
-            # Adverse gradient - flow is running uphill
-            # Flow will be purely controlled by uphill outlet face
-            print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
-
-
-        s = 'Delta total energy = %.3f' %(delta_Et)
-        log_to_file(log_filename, s)
-
-
-        # Calculate discharge for one barrel
-        Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
-        
-        # Adjust discharge for multiple barrels
-        Q *= self.number_of_barrels
-
-        # Compute barrel momentum
-        barrel_momentum = barrel_velocity*culvert_outlet_depth
-                
-        s = 'Barrel velocity = %f' %barrel_velocity
-        log_to_file(log_filename, s)
-
-        # Compute momentum vector at outlet
-        outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
-            
-        s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
-        log_to_file(log_filename, s)
-
-        # Log timeseries to file
-        fid = open(self.timeseries_filename, 'a')        
-        fid.write('%f, %f, %f, %f, %f\n'\
-                      %(time, 
-                        openings[0].total_energy,
-                        openings[1].total_energy,
-                        barrel_velocity,
-                        Q))
-        fid.close()
-
-        # Update momentum        
-        delta_t = time - self.last_time
-        if delta_t > 0.0:
-            xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
-            xmomentum_rate /= delta_t
-                
-            ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
-            ymomentum_rate /= delta_t
-                    
-            s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
-            log_to_file(log_filename, s)                    
-        else:
-            xmomentum_rate = ymomentum_rate = 0.0
-
-
-        # Set momentum rates for outlet jet
-        outlet.momentum[0].rate = xmomentum_rate
-        outlet.momentum[1].rate = ymomentum_rate
-
-        # Remember this value for next step (IMPORTANT)
-        outlet.momentum[0].value = outlet_mom_x
-        outlet.momentum[1].value = outlet_mom_y                    
-
-        if int(domain.time*100) % 100 == 0:
-            s = 'T=%.5f, Culvert Discharge = %.3f f'\
-                %(time, Q)
-            s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
-                 %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
-            s += ' Momentum rate: (%.4f, %.4f)'\
-                 %(xmomentum_rate, ymomentum_rate)                    
-            s+='Outlet Vel= %.3f'\
-                %(barrel_velocity)
-            log_to_file(log_filename, s)
-            
-                
-
-
-       
+                
+
+
         # Execute flow term for each opening
         # This is where Inflow objects are evaluated and update the domain
@@ -420,8 +445,8 @@
         # Execute momentum terms
         # This is where Inflow objects are evaluated and update the domain
-        outlet.momentum[0](domain)
-        outlet.momentum[1](domain)        
+        self.outlet.momentum[0](domain)
+        self.outlet.momentum[1](domain)        
             
-        # Store value of time    
+        # Store value of time #FIXME(Ole): Maybe only every time we update   
         self.last_time = time
 
Index: /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 6050)
+++ /anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py	(revision 6051)
@@ -2827,4 +2827,66 @@
 
 
+        
+
+
+    def test_rainfall_forcing_with_evolve(self):
+        """test_rainfall_forcing_with_evolve
+
+        Test how forcing terms are called within evolve
+        """
+        
+        # FIXME(Ole): This test is just to experiment
+        
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe
+        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+
+        domain = Domain(points, vertices)
+
+        #Flat surface with 1m of water
+        domain.set_quantity('elevation', 0)
+        domain.set_quantity('stage', 1.0)
+        domain.set_quantity('friction', 0)
+
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'exterior': Br})
+
+        # Setup only one forcing term, time dependent rainfall that expires at t==20
+        from anuga.fit_interpolate.interpolate import Modeltime_too_late
+        def main_rate(t):
+            if t > 20:
+                msg = 'Model time exceeded.'
+                raise Modeltime_too_late, msg
+            else:
+                return 3*t + 7
+        
+        domain.forcing_terms = []
+        R = Rainfall(domain,
+                     rate=main_rate,
+                     polygon=[[1,1], [2,1], [2,2], [1,2]],
+                     default_rate=5.0)
+
+        assert allclose(R.exchange_area, 1)
+        
+        domain.forcing_terms.append(R)
+
+        for t in domain.evolve(yieldstep=1, finaltime=25):
+            pass
+            
+            #print t, domain.quantities['stage'].explicit_update, (3*t+7)/1000
+            
+            #FIXME(Ole):  A test here is hard because explicit_update also
+            # receives updates from the flux calculation.
+
+        
+        
 
     def test_inflow_using_circle(self):
@@ -6385,5 +6447,6 @@
 if __name__ == "__main__":
 
-    suite = unittest.makeSuite(Test_Shallow_Water,'test')
+    suite = unittest.makeSuite(Test_Shallow_Water, 'test')
+    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_rainfall_forcing_with_evolve')
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_energy_through_cross_section_with_g')    
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_fitting_using_shallow_water_domain')    
