Index: /inundation/ga/storm_surge/pyvolution/domain.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/domain.py	(revision 475)
+++ /inundation/ga/storm_surge/pyvolution/domain.py	(revision 476)
@@ -45,34 +45,7 @@
             self.quantities[name] = Quantity(self)
 
-
-        #FIXME: Move these explanations elsewhere    
-            
         #Create an empty list for explicit forcing terms
-        # 
-        # Explicit terms must have the form
-        #
-        #    G(q, t)
-        #
-        # and explicit scheme is
-        #
-        #   q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
-	#
-	#
-	# FIXME: How to call and how function should look
-
         self.forcing_terms = []
 
-
-        #Create an empty list for semi implicit forcing terms if any
-        # 
-        # Semi implicit forcing terms are assumed to have the form
-        #
-        #    G(q, t) = H(q, t) q
-        #
-        # and the semi implicit scheme will then be
-        #
-        #   q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
-        
-        ###self.semi_implicit_forcing_terms = []
 
 
Index: /inundation/ga/storm_surge/pyvolution/quantity.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/quantity.py	(revision 475)
+++ /inundation/ga/storm_surge/pyvolution/quantity.py	(revision 476)
@@ -472,4 +472,34 @@
     """Update centroid values based on values stored in
     explicit_update and semi_implicit_update as well as given timestep
+
+    Function implementing forcing terms must take on argument
+    which is the domain and they must update either explicit
+    or implicit updates, e,g,:
+
+    def gravity(domain):
+        ....
+        domain.quantities['xmomentum'].explicit_update = ...
+        domain.quantities['ymomentum'].explicit_update = ...
+
+        
+
+    Explicit terms must have the form
+    
+        G(q, t)
+    
+    and explicit scheme is
+    
+       q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
+
+
+    Semi implicit forcing terms are assumed to have the form
+    
+       G(q, t) = H(q, t) q
+    
+    and the semi implicit scheme will then be
+    
+      q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
+
+    
     """
     
@@ -479,5 +509,5 @@
 
 
-    #Divide by semi_implicit update by conserved quantity
+    #Divide H by conserved quantity to obtain G (see docstring above)
     for k in range(N):
         x = quantity.centroid_values[k] 
@@ -487,5 +517,4 @@
             quantity.semi_implicit_update[k] /= x             
             
-
         
     #Explicit updates
Index: /inundation/ga/storm_surge/pyvolution/test_quantity.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_quantity.py	(revision 475)
+++ /inundation/ga/storm_surge/pyvolution/test_quantity.py	(revision 476)
@@ -7,5 +7,5 @@
 from quantity import *
 from config import epsilon
-from Numeric import allclose, array
+from Numeric import allclose, array, ones, Float
 
         
@@ -412,40 +412,48 @@
         assert allclose( quantity.centroid_values, x)
 
-    #FIXME: Update these tests once I understand the semi_implicit scheme    
-#     def test_update_semi_implicit(self):
-#         quantity = Conserved_quantity(self.mesh4)
-
-#         #Test centroids
-#         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-#         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
-
-#         #Set semi implicit update
-#         quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
-
-#         #Update with given timestep
-#         quantity.update(0.1)
-
-#         x = array([1, 2, 3, 4])/array( [.9,.9,.9,.9] )
-#         assert allclose( quantity.centroid_values, x)
-
-#     def test_both_updates(self):
-#         quantity = Conserved_quantity(self.mesh4)
-
-#         #Test centroids
-#         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
-#         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
-
-#         #Set explicit_update
-#         quantity.explicit_update = array( [4.,3.,2.,1.] )
-        
-#         #Set semi implicit update
-#         quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
-
-#         #Update with given timestep
-#         quantity.update(0.1)
-
-#         x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
-#         x /= array( [.9,.9,.9,.9] )
-#         assert allclose( quantity.centroid_values, x)        
+    def test_update_semi_implicit(self):
+        quantity = Conserved_quantity(self.mesh4)
+
+        #Test centroids
+        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+
+        #Set semi implicit update
+        quantity.semi_implicit_update = array([1.,1.,1.,1.])
+
+        #Update with given timestep
+        timestep = 0.1
+        quantity.update(timestep)
+
+        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
+        denom = ones(4, Float)-timestep*sem
+        
+        x = array([1, 2, 3, 4])/denom
+        assert allclose( quantity.centroid_values, x)
+
+
+    def test_both_updates(self):
+        quantity = Conserved_quantity(self.mesh4)
+
+        #Test centroids
+        quantity.set_values([1.,2.,3.,4.], location = 'centroids')
+        assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid
+
+        #Set explicit_update
+        quantity.explicit_update = array( [4.,3.,2.,1.] )
+       
+        #Set semi implicit update
+        quantity.semi_implicit_update = array( [1.,1.,1.,1.] )
+
+        #Update with given timestep
+        timestep = 0.1        
+        quantity.update(0.1)
+
+        sem = array([1.,1.,1.,1.])/array([1, 2, 3, 4])
+        denom = ones(4, Float)-timestep*sem
+
+        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
+        x /= denom
+        assert allclose( quantity.centroid_values, x)        
 
         
