Index: /trunk/anuga_core/source/anuga_parallel/parallel_boyd_box_operator.py
===================================================================
--- /trunk/anuga_core/source/anuga_parallel/parallel_boyd_box_operator.py	(revision 9130)
+++ /trunk/anuga_core/source/anuga_parallel/parallel_boyd_box_operator.py	(revision 9131)
@@ -1,4 +1,5 @@
 import anuga
 import math
+import numpy
 
 from anuga.structures.boyd_box_operator import boyd_box_function 
@@ -100,13 +101,11 @@
         # May/June 2014 -- allow 'smoothing ' of driving_energy, delta total energy, and outflow_enq_depth
         self.smoothing_timescale=0.
-        self.smooth_driving_energy=0.
         self.smooth_delta_total_energy=0.
-        self.smooth_outflow_enq_depth=0.
+        self.smooth_Q=0.
         # Set them based on a call to the discharge routine with smoothing_timescale=0.
         # [values of self.smooth_* are required in discharge_routine, hence dummy values above]
         Qvd=self.discharge_routine()
-        self.smooth_driving_energy=1.0*self.driving_energy
         self.smooth_delta_total_energy=1.0*self.delta_total_energy
-        self.smooth_outflow_enq_depth=Qvd[2]
+        self.smooth_Q=Qvd[0]
         # Finally, set the smoothing timescale we actually want
         self.smoothing_timescale=smoothing_timescale
@@ -170,12 +169,16 @@
         # master proc orders reversal if applicable
         if self.myid == self.master_proc:
-
+            # May/June 2014 -- change the driving forces gradually, with forward euler timestepping 
+            ts=self.domain.timestep/max(self.domain.timestep, self.smoothing_timescale,1.0e-06)
+            self.smooth_delta_total_energy=self.smooth_delta_total_energy+\
+                                    ts*(self.delta_total_energy-self.smooth_delta_total_energy)
 
             # Reverse the inflow and outflow direction?
-            if self.delta_total_energy < 0:
+            if self.smooth_delta_total_energy < 0:
                 self.inflow_index = 1
                 self.outflow_index = 0
 
-                self.delta_total_energy = -self.delta_total_energy
+                #self.delta_total_energy = -self.delta_total_energy
+                self.delta_total_energy = -self.smooth_delta_total_energy
 
                 for i in self.procs:
@@ -183,4 +186,5 @@
                     pypar.send(True, i)
             else:
+                self.delta_total_energy = self.smooth_delta_total_energy
                 for i in self.procs:
                     if i == self.master_proc: continue
@@ -251,13 +255,4 @@
                     self.driving_energy = inflow_enq_depth
                     
-                # May/June 2014 -- change the driving forces gradually, with forward euler timestepping 
-                ts=self.domain.timestep/max(self.domain.timestep, self.smoothing_timescale,1.0e-06)
-                self.smooth_driving_energy=self.smooth_driving_energy+\
-                                        ts*(self.driving_energy-self.smooth_driving_energy)
-                self.smooth_delta_total_energy=self.smooth_delta_total_energy+\
-                                        ts*(self.delta_total_energy-self.smooth_delta_total_energy)
-                self.smooth_outflow_enq_depth=self.smooth_outflow_enq_depth+\
-                                        ts*(outflow_enq_depth-self.smooth_outflow_enq_depth)
-
             
                 Q, barrel_velocity, outlet_culvert_depth, flow_area, case = \
@@ -266,15 +261,27 @@
                                                 flow_width          =self.culvert_width,
                                                 length              =self.culvert_length,
-                                                #driving_energy      =self.driving_energy,
-                                                #delta_total_energy  =self.delta_total_energy,
-                                                #outlet_enquiry_depth=outflow_enq_depth,
-                                                # Allow smoothing of the driving energies 
-                                                driving_energy      =self.smooth_driving_energy,
-                                                delta_total_energy  =self.smooth_delta_total_energy,
-                                                outlet_enquiry_depth=self.smooth_outflow_enq_depth,
+                                                driving_energy      =self.driving_energy,
+                                                delta_total_energy  =self.delta_total_energy,
+                                                outlet_enquiry_depth=outflow_enq_depth,
                                                 sum_loss            =self.sum_loss,
                                                 manning             =self.manning)
 
-                
+                ################################################
+                # Smooth discharge. This can reduce oscillations
+                # 
+                # NOTE: The sign of smooth_Q assumes that
+                #   self.inflow_index=0 and self.outflow_index=1
+                #   , whereas the sign of Q is always positive
+                Qsign=(self.outflow_index-self.inflow_index) # To adjust sign of Q
+                self.smooth_Q = self.smooth_Q +ts*(Q*Qsign-self.smooth_Q)
+                if numpy.sign(self.smooth_Q)!=Qsign:
+                    # The flow direction of the 'instantaneous Q' based on the
+                    # 'smoothed delta_total_energy' is not the same as the
+                    # direction of smooth_Q. To prevent 'jumping around', let's
+                    # set Q to zero
+                    Q=0.
+                else:
+                    Q = abs(self.smooth_Q)
+                barrel_velocity=Q/flow_area
             # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
 
