Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 4727)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py	(revision 4728)
@@ -150,5 +150,7 @@
         Pre-condition: vertex_values have been set
         """
-
+        # FIXME (Ole): This function is not called
+        # in real model runs. Consider removing it.
+        
         N = self.vertex_values.shape[0]
         for i in range(N):
Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 4727)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c	(revision 4728)
@@ -36,8 +36,8 @@
 
 		if (number_of_boundaries[k] < 2) {
-			//Two or three true neighbours
-
-			//Get indices of neighbours (or self when used as surrogate)
-			//k0, k1, k2 = surrogate_neighbours[k,:]
+			// Two or three true neighbours
+
+			// Get indices of neighbours (or self when used as surrogate)
+			// k0, k1, k2 = surrogate_neighbours[k,:]
 
 			k0 = surrogate_neighbours[index3 + 0];
@@ -48,5 +48,5 @@
 			if (k0 == k1 || k1 == k2) return -1;
 
-			//Get data
+			// Get data
 			q0 = centroid_values[k0];
 			q1 = centroid_values[k1];
@@ -57,11 +57,11 @@
 			x2 = centroids[k2*2]; y2 = centroids[k2*2+1];
 
-			//Gradient
+			// Gradient
 			_gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
 
 		} else if (number_of_boundaries[k] == 2) {
-			//One true neighbour
-
-			//Get index of the one neighbour
+			// One true neighbour
+
+			// Get index of the one neighbour
 			i=0; k0 = k;
 			while (i<3 && k0==k) {
@@ -73,5 +73,5 @@
 			k1 = k; //self
 
-			//Get data
+			// Get data
 			q0 = centroid_values[k0];
 			q1 = centroid_values[k1];
@@ -80,19 +80,11 @@
 			x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
 
-			//Two point gradient
+			// Two point gradient
 			_gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);
 
-
-			//Old (wrong code)
-			//det = x0*y1 - x1*y0;
-			//if (det != 0.0) {
-			//	a[k] = (y1*q0 - y0*q1)/det;
-			//	b[k] = (x0*q1 - x1*q0)/det;
-			//}
 		}
 		//    else:
 		//        #No true neighbours -
 		//        #Fall back to first order scheme
-		//        pass
 	}
 	return 0;
@@ -117,9 +109,9 @@
 		k2 = 2*k;
 
-		//Centroid coordinates
+		// Centroid coordinates
 		x = centroids[k2]; y = centroids[k2+1];
 
-		//vertex coordinates
-		//x0, y0, x1, y1, x2, y2 = X[k,:]
+		// vertex coordinates
+		// x0, y0, x1, y1, x2, y2 = X[k,:]
 		x0 = vertex_coordinates[k6 + 0];
 		y0 = vertex_coordinates[k6 + 1];
@@ -129,5 +121,5 @@
 		y2 = vertex_coordinates[k6 + 5];
 
-		//Extrapolate
+		// Extrapolate
 		vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
 		vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
@@ -156,5 +148,4 @@
 		q2 = vertex_values[k3 + 2];
 
-		//printf("%f, %f, %f\n", q0, q1, q2);
 		edge_values[k3 + 0] = 0.5*(q1+q2);
 		edge_values[k3 + 1] = 0.5*(q0+q2);
@@ -167,5 +158,5 @@
 			    double* centroid_values,
 			    double* centroid_backup_values) {
-    //Backup centroid values
+    // Backup centroid values
 
 
@@ -186,5 +177,5 @@
 			   double* centroid_values,
 			   double* centroid_backup_values) {
-    //saxby centroid values
+    // Saxby centroid values
 
 
@@ -206,6 +197,6 @@
 	    double* explicit_update,
 	    double* semi_implicit_update) {
-	//Update centroid values based on values stored in
-	//explicit_update and semi_implicit_update as well as given timestep
+	// Update centroid values based on values stored in
+	// explicit_update and semi_implicit_update as well as given timestep
 
 
@@ -214,5 +205,5 @@
 
 
-	//Divide semi_implicit update by conserved quantity
+	// Divide semi_implicit update by conserved quantity
 	for (k=0; k<N; k++) {
 		x = centroid_values[k];
@@ -225,5 +216,5 @@
 
 
-	//Semi implicit updates
+	// Semi implicit updates
 	for (k=0; k<N; k++) {
 		denominator = 1.0 - timestep*semi_implicit_update[k];
@@ -236,10 +227,6 @@
 	}
 
-	/*  for (k=0; k<N; k++) {*/
-	/*    centroid_values[k] = exp(timestep*semi_implicit_update[k])*centroid_values[k];*/
-	/*  }*/
-
-
-	//Explicit updates
+
+	// Explicit updates
 	for (k=0; k<N; k++) {
 		centroid_values[k] += timestep*explicit_update[k];
@@ -247,5 +234,5 @@
 
 
-	//MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
+	// MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
 	for (k=0;k<N;k++){
 		semi_implicit_update[k]=0.0;
@@ -261,5 +248,5 @@
 			   double* vertex_values, 
 			   double* A) {
-  //Average vertex values to obtain one value per node
+  // Average vertex values to obtain one value per node
 
   int i, index; 
@@ -270,7 +257,7 @@
   for (i=0; i<N; i++) {
   
-    //if (current_node == N) {
-    //  printf("Current node exceeding number of nodes (%d)", N);    
-    //  return 1;
+    // if (current_node == N) {
+    //   printf("Current node exceeding number of nodes (%d)", N);    
+    //   return 1;
     // }
     
@@ -278,10 +265,10 @@
     k += 1;
             
-    //volume_id = index / 3
-    //vertex_id = index % 3
-    //total += self.vertex_values[volume_id, vertex_id]
+    // volume_id = index / 3
+    // vertex_id = index % 3
+    // total += self.vertex_values[volume_id, vertex_id]
     total += vertex_values[index];
       
-    //printf("current_node=%d, index=%d, k=%d, total=%f\n", current_node, index, k, total);
+    // printf("current_node=%d, index=%d, k=%d, total=%f\n", current_node, index, k, total);
     if (number_of_triangles_per_node[current_node] == k) {
       A[current_node] = total/k;
@@ -334,5 +321,5 @@
 	}
 
-	//Release and return
+	// Release and return
 	Py_DECREF(centroid_values);
 	Py_DECREF(explicit_update);
@@ -368,5 +355,5 @@
 
 
-	//Release and return
+	// Release and return
 	Py_DECREF(centroid_values);
 	Py_DECREF(centroid_backup_values);
@@ -401,5 +388,5 @@
 
 
-	//Release and return
+	// Release and return
 	Py_DECREF(centroid_values);
 	Py_DECREF(centroid_backup_values);
@@ -438,5 +425,5 @@
 	}
 
-	//Release and return
+	// Release and return
 	Py_DECREF(vertex_values);
 	Py_DECREF(edge_values);
@@ -514,5 +501,5 @@
 	}
 
-	//Get pertinent variables
+	// Get pertinent variables
 
 	centroids = get_consecutive_array(domain, "centroid_coordinates");
@@ -523,8 +510,8 @@
 	N = centroid_values -> dimensions[0];
 
-	//Release
+	// Release
 	Py_DECREF(domain);
 
-	//Allocate space for return vectors a and b (don't DECREF)
+	// Allocate space for return vectors a and b (don't DECREF)
 	dimensions[0] = N;
 	a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
@@ -546,5 +533,5 @@
 	}
 
-	//Release
+	// Release
 	Py_DECREF(centroids);
 	Py_DECREF(centroid_values);
@@ -552,5 +539,5 @@
 	Py_DECREF(surrogate_neighbours);
 
-	//Build result, release and return
+	// Build result, release and return
 	R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b));
 	Py_DECREF(a);
@@ -591,5 +578,5 @@
 	}
 
-	//Get pertinent variables
+	// Get pertinent variables
 	centroids = get_consecutive_array(domain, "centroid_coordinates");
 	centroid_values = get_consecutive_array(quantity, "centroid_values");
@@ -601,5 +588,5 @@
 	N = centroid_values -> dimensions[0];
 
-	//Release
+	// Release
 	Py_DECREF(domain);
 
@@ -647,5 +634,5 @@
 
 
-	//Release
+	// Release
 	Py_DECREF(centroids);
 	Py_DECREF(centroid_values);
@@ -692,5 +679,5 @@
 	neighbours = get_consecutive_array(domain, "neighbours");
 
-	//Get safety factor beta_w
+	// Get safety factor beta_w
 	Tmp = PyObject_GetAttrString(domain, "beta_w");
 	if (!Tmp) {
@@ -713,5 +700,5 @@
 	N = qc -> dimensions[0];
 
-	//Find min and max of this and neighbour's centroid values
+	// Find min and max of this and neighbour's centroid values
 	qmin = malloc(N * sizeof(double));
 	qmax = malloc(N * sizeof(double));
@@ -765,4 +752,4 @@
   Py_InitModule("quantity_ext", MethodTable);
 
-  import_array();     //Necessary for handling of NumPY structures
-}
+  import_array(); // Necessary for handling of NumPY structures
+}
