Index: /inundation/ga/storm_surge/Merimbula/run_merimbula_lake.py
===================================================================
--- /inundation/ga/storm_surge/Merimbula/run_merimbula_lake.py	(revision 1362)
+++ /inundation/ga/storm_surge/Merimbula/run_merimbula_lake.py	(revision 1363)
@@ -90,6 +90,6 @@
 # Evolution
 t0 = time.time()
-domain.visualise = False
-yieldstep = 900
+domain.visualise = True
+yieldstep = 10
 finaltime = 588000*3
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
Index: /inundation/ga/storm_surge/parallel/advection.py
===================================================================
--- /inundation/ga/storm_surge/parallel/advection.py	(revision 1362)
+++ /inundation/ga/storm_surge/parallel/advection.py	(revision 1363)
@@ -45,6 +45,9 @@
         #Realtime visualisation
         self.visualise = False
+        self.visualise_color_stage = False
+        self.visualise_timer = True
+        self.visualise_range_z = None
         self.smooth = True
-        self.range_z = None
+
 
 
@@ -177,5 +180,5 @@
         if self.visualise is True and self.time == 0.0:
             import realtime_visualisation_new as visualise
-            visualise.create_surface(self,range_z =  self.range_z)
+            visualise.create_surface(self)
 
         #Call basic machinery from parent class
Index: /inundation/ga/storm_surge/parallel/demo.py
===================================================================
--- /inundation/ga/storm_surge/parallel/demo.py	(revision 1363)
+++ /inundation/ga/storm_surge/parallel/demo.py	(revision 1363)
@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+#########################################################
+#   
+#  Example: Running Python in parallel using pypar (MPI).
+#   
+#  Author:  Ole Nielsen,  SMS, ANU, Jan. 2002. 
+#
+#########################################################
+#
+# The purpose of this code is to demonstrate how Python can be
+# used to communicate among processes using pypar.
+#
+# This demo passes messages on in a ring from processor n-1 to n starting
+# and ending with processor 0.
+# Each processor adds some text to the message before passing it on
+# and writes a log statement to the screen.
+#
+# To execute:
+#
+#   mpirun -np 4 demo.py
+
+
+import pypar    # The Python-MPI interface 
+
+numproc = pypar.size()
+myid =    pypar.rank()
+node =    pypar.get_processor_name()
+
+print "I am proc %d of %d on node %s" %(myid, numproc, node)
+
+
+if numproc < 2:
+  print "Demo must run on at least 2 processors to continue"      
+  pypar.Abort()
+  
+if myid == 0:
+  msg = "MSGP0"  
+  
+  print 'Processor 0 sending message "%s" to processor %d' %(msg, 1)
+  pypar.send(msg, 1)
+
+  msg, status = pypar.receive(numproc-1, return_status=True)
+  print 'Processor 0 received message "%s" from processor %d' %(msg, numproc-1)
+  print 'Size of msg was %d bytes' %(status.bytes())
+
+else:
+  source = myid-1
+  destination = (myid+1)%numproc
+  
+  msg, status = pypar.receive(source, return_status=True)
+  print 'Processor %d received message "%s" from processor %d'\
+        %(myid, msg, source)
+  print 'Size of msg was %d bytes' %(status.bytes())  
+
+  msg = msg + '->P' + str(myid) #Update message     
+  print 'Processor %d sending msg "%s" to %d' %(myid, msg, destination)
+  pypar.send(msg, destination)
+
+pypar.finalize()
+
+
+
Index: /inundation/ga/storm_surge/parallel/demo2.py
===================================================================
--- /inundation/ga/storm_surge/parallel/demo2.py	(revision 1363)
+++ /inundation/ga/storm_surge/parallel/demo2.py	(revision 1363)
@@ -0,0 +1,107 @@
+#!/usr/bin/env python
+#########################################################
+#   
+#  Example: Distributing work in parallel using pypar (MPI).
+#   
+#  Author:  Ole Nielsen, ANU, March 2002. 
+#
+#########################################################
+#
+"""Skeleton for using pypar to distribute work across processors.
+   It is assumed that the function work knows how to perform different 
+   tasks (or the same task on different data) based on values of
+   myid and numprocs.
+"""
+
+
+	
+def work(myid, numprocs, data):
+  """This simple example function that slices up the data 
+     based on values of numproc and myid. 
+  """
+  
+  import Numeric
+  
+  # Identify local slice and process it
+  #
+  interval = len(data)
+  myinterval = interval/numprocs
+  
+  mylower = myid*myinterval
+  if myid == numprocs-1:
+    myupper = interval+1
+  else:  
+    myupper = mylower + myinterval
+  
+  mydata = data[mylower:myupper]
+  
+  # Computation (average)
+  #
+  myavg = float(Numeric.sum(mydata))/len(mydata)
+  print "P%d: %s Local avg=%.4f" %(myid, str(mydata), myavg)
+
+  return myavg*len(mydata)
+  
+  
+  
+  
+###################################################
+# Main program - communication takes place here
+#
+import pypar, Numeric
+
+
+# Get data. Here it is just generated but it could be read 
+# from file or given as an input parameter.
+# 
+lower = 100
+upper = 121
+data = Numeric.array(range(lower,upper))
+
+#
+# Get parallel parameters
+#
+numprocs = pypar.size()    # Number of processors
+myid = pypar.rank()        # Id of this processor
+node = pypar.get_processor_name()
+
+print "I am proc %d of %d on node %s" %(myid, numprocs, node)
+
+
+
+#
+# Do work in parallel
+#
+x = work(myid, numprocs, data)   #Do work on all processors
+print "Proc %d finished working" %myid
+
+#  
+# Communication
+#
+if numprocs > 1:
+  #
+  # Processor 0 gathers all results and merge them
+  #
+  if myid == 0:
+    for id in range(1,numprocs):
+      print "P%d receving from P%d" %(0, id)
+      x = x + pypar.receive(id)  #Add up (would be more complex in general)
+  
+  # All other processors send their results back to processor 0
+  #  
+  else:
+    print "P%d sending to P%d" %(myid, 0)  
+    pypar.send(x, 0)
+
+print "Proc %d after communication" %myid    
+#
+# Compute overall average and report
+#    
+  
+if myid == 0:  
+  avg = x/len(data)     
+  print "Global average is %.4f" %avg      
+
+pypar.finalize()    
+
+
Index: /inundation/ga/storm_surge/parallel/demo3.py
===================================================================
--- /inundation/ga/storm_surge/parallel/demo3.py	(revision 1363)
+++ /inundation/ga/storm_surge/parallel/demo3.py	(revision 1363)
@@ -0,0 +1,125 @@
+#!/usr/bin/env python
+#################################################################
+# Master/Slave Parallel decomposition sample 
+# 
+# Run as 
+#   python demo3.py
+# or 
+#   mpirun -np 2 demo3.py
+# (perhaps try number of processors more than 2)
+#################################################################
+#
+# To verify bandwidth of your architexture please 
+# run pytiming (and ctiming) 
+#
+# OMN, GPC FEB 2002
+#
+#
+
+import sys
+
+try:
+  import Numeric
+except:
+  raise 'Module Numeric must be present to run pypar'
+  
+try:
+  import pypar
+except:
+  raise 'Module pypar must be present to run parallel'
+
+sys.stderr.write("Modules Numeric, pypar imported OK\n")
+
+WORKTAG = 1
+DIETAG =  2
+
+
+def master():
+    numCompleted = 0
+    
+    sys.stderr.write("[MASTER]: I am processor %d of %d on node %s\n" %(MPI_myid, MPI_numproc, MPI_node))
+    
+    # start slaves distributing the first work slot
+    for i in range(1, min(MPI_numproc, numWorks)): 
+        work = workList[i]
+        pypar.raw_send(work, i, WORKTAG) 
+        sys.stderr.write("[MASTER]: sent work '%s' to node '%d'\n" %(work, i))
+
+    # dispach the remaining work slots on dynamic load-balancing policy
+    # the quicker to do the job, the more jobs it takes
+    for work in workList[MPI_numproc:]:
+        result = '  '
+        err, status = pypar.raw_receive(result, pypar.any_source, pypar.any_tag, return_status=True) 
+        #sys.stderr.write( "[MASTER]: received result '%s' from node '%d'\n" %(result, err[1][0]))
+        sys.stderr.write("[MASTER]: received result '%s' from node '%d'\n" %(result, status.source))
+        numCompleted += 1
+        pypar.raw_send(work, status.source, WORKTAG)
+        sys.stderr.write("[MASTER]: sent work '%s' to node '%d'\n" %(work, status.source))
+    
+    # all works have been dispatched out
+    sys.stderr.write("[MASTER]: toDo : %d\n" %numWorks)
+    sys.stderr.write("[MASTER]: done : %d\n" %numCompleted)
+    
+    # I've still to take into the remaining completions   
+    while(numCompleted < numWorks): 
+        result = '  '
+        err, status = pypar.raw_receive(result, pypar.any_source, pypar.any_tag, return_status=True) 
+        sys.stderr.write("[MASTER]: received (final) result '%s' from node '%d'\n" %(result, status.source))
+        numCompleted += 1
+        sys.stderr.write("[MASTER]: %d completed\n" %numCompleted)
+        
+    sys.stderr.write( "[MASTER]: about to terminate slaves\n")
+
+    # say slaves to stop working
+    for i in range(1, MPI_numproc): 
+        pypar.raw_send('#', i, DIETAG) 
+        sys.stderr.write("[MASTER]: sent (final) work '%s' to node '%d'\n" %(0, i))
+        
+    return
+    
+def slave():
+
+    sys.stderr.write( "[SLAVE %d]: I am processor %d of %d on node %s\n" %(MPI_myid, MPI_myid, MPI_numproc, MPI_node))
+
+    while 1:
+        result = ' '
+        err, status = pypar.raw_receive(result, pypar.any_source, pypar.any_tag, return_status=True) 
+        sys.stderr.write("[SLAVE %d]: received work '%s' with tag '%d' from node '%d'\n"\
+	      %(MPI_myid, result, status.tag, status.source))
+       
+        if (status.tag == DIETAG):
+            sys.stderr.write("[SLAVE %d]: received termination from node '%d'\n" %(MPI_myid, 0))
+            return
+        else:
+            result = 'X'+result
+            pypar.raw_send(result, 0)
+            sys.stderr.write("[SLAVE %d]: sent result '%s' to node '%d'\n" %(MPI_myid, result, 0))
+            
+       
+
+if __name__ == '__main__':
+    MPI_myid =    pypar.rank()
+    MPI_numproc = pypar.size()
+    MPI_node =    pypar.Get_processor_name()
+
+    _workList = ('_dummy_', 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j')
+    workList = ('_dummy_', 'a', 'b', 'c')
+    numWorks = len(workList) - 1
+    
+    
+    #FIXME, better control here
+    if MPI_numproc > numWorks or MPI_numproc < 2:
+        pypar.Finalize()
+	if MPI_myid == 0:
+          sys.stderr.write("ERROR: Number of processors must be in the interval [2,%d].\n" %numWorks)
+	  
+        sys.exit(-1)
+
+    if MPI_myid == 0:
+        master()
+    else:
+        slave()
+
+    pypar.Finalize()
+    sys.stderr.write("MPI environment finalized.\n")
+        	
Index: /inundation/ga/storm_surge/parallel/demo4.py
===================================================================
--- /inundation/ga/storm_surge/parallel/demo4.py	(revision 1363)
+++ /inundation/ga/storm_surge/parallel/demo4.py	(revision 1363)
@@ -0,0 +1,66 @@
+#!/usr/bin/env python
+#########################################################
+#   
+#  Example: Running Python in parallel using pypar (MPI).
+#   
+#  Author:  Ole Nielsen,  SMS, ANU, Jan. 2002. 
+#
+#########################################################
+#
+# The purpose of this code is to demonstrate how Python can be
+# used to communicate among processes using pypar.
+#
+# This demo passes messages on in a ring from processor n-1 to n starting
+# and ending with processor 0.
+# Each processor adds some text to the message before passing it on
+# and writes a log statement to the screen.
+#
+# To execute:
+#
+#   mpirun -np 4 demo.py
+
+
+import pypar    # The Python-MPI interface 
+import Numeric
+typ = 'd'  #'i', 'l', 'f', 'd'
+
+numproc = pypar.size()
+myid =    pypar.rank()
+node =    pypar.Get_processor_name()
+
+print "I am proc %d of %d on node %s" %(myid, numproc, node)
+
+if numproc < 2:
+  print "Demo must run on at least 2 processors to continue"      
+  pypar.Abort()
+  
+if myid == 0:
+  msg = Numeric.array([0], typ)  
+  
+  print 'Processor 0 sending message "%s" to processor %d' %(str(msg), 1)
+  pypar.send(msg, 1)
+
+  msg, status = pypar.receive(numproc-1, return_status=True)
+  print 'Processor 0 received message "%s" from processor %d' %(str(msg), numproc-1)
+  print 'Size of msg was %d bytes' %(status.bytes())
+
+else:
+  source = myid-1
+  destination = (myid+1)%numproc
+  
+  msg, status = pypar.receive(source, return_status=True)
+  print 'Processor %d received message "%s" from processor %d'\
+        %(myid, str(msg), source)
+  print 'Size of msg was %d bytes' %(status.bytes())  
+
+  msg = list(msg)
+  msg.append(myid)
+  print msg
+  msg = Numeric.array(msg, typ)  #Update message     
+  print 'Processor %d sending msg "%s" to %d' %(myid, msg, destination)
+  pypar.send(msg, destination)
+
+pypar.Finalize()
+
+
+
Index: /inundation/ga/storm_surge/parallel/parallel_advection.py
===================================================================
--- /inundation/ga/storm_surge/parallel/parallel_advection.py	(revision 1363)
+++ /inundation/ga/storm_surge/parallel/parallel_advection.py	(revision 1363)
@@ -0,0 +1,184 @@
+import sys
+from os import sep
+sys.path.append('..'+sep+'pyvolution')
+
+"""Class Parallel_Domain -
+2D triangular domains for finite-volume computations of
+the advection equation, with extra structures to allow
+communication between other Parallel_Domains and itself
+
+This module contains a specialisation of class Domain from module advection.py
+
+Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+Geoscience Australia, 2004
+"""
+
+#Subversion keywords:
+#
+#$LastChangedDate:  $
+#$LastChangedRevision:  $
+#$LastChangedBy: $
+
+class Parallel_Domain(Domain):
+
+    def __init__(self, coordinates, vertices, boundary = None, velocity = None):
+
+
+        Generic_domain.__init__(self, coordinates, vertices, boundary,
+                                ['stage'])
+
+        if velocity is None:
+            self.velocity = [1,0]
+        else:
+            self.velocity = velocity
+
+        #Only first is implemented for advection
+        self.default_order = self.order = 1
+
+        #Realtime visualisation
+        self.visualise = False
+        self.visualise_color_stage = False
+        self.visualise_timer = True
+        self.visualise_range_z = None
+        self.smooth = True
+
+
+
+    def check_integrity(self):
+        Generic_domain.check_integrity(self)
+
+        msg = 'Conserved quantity must be "stage"'
+        assert self.conserved_quantities[0] == 'stage', msg
+
+
+    def flux_function(self, normal, ql, qr, zl=None, zr=None):
+        """Compute outward flux as inner product between velocity
+        vector v=(v_1, v_2) and normal vector n.
+
+        if <n,v> > 0 flux direction is outward bound and its magnitude is
+        determined by the quantity inside volume: ql.
+        Otherwise it is inbound and magnitude is determined by the
+        quantity outside the volume: qr.
+        """
+
+        v1 = self.velocity[0]
+        v2 = self.velocity[1]
+
+
+        normal_velocity = v1*normal[0] + v2*normal[1]
+
+        if normal_velocity < 0:
+            flux = qr * normal_velocity
+        else:
+            flux = ql * normal_velocity
+
+        max_speed = abs(normal_velocity)
+        return flux, max_speed
+
+
+    def compute_fluxes(self):
+        """Compute all fluxes and the timestep suitable for all volumes
+        in domain.
+
+        Compute total flux for each conserved quantity using "flux_function"
+
+        Fluxes across each edge are scaled by edgelengths and summed up
+        Resulting flux is then scaled by area and stored in
+        domain.explicit_update
+
+        The maximal allowable speed computed by the flux_function
+        for each volume
+        is converted to a timestep that must not be exceeded. The minimum of
+        those is computed as the next overall timestep.
+
+        Post conditions:
+        domain.explicit_update is reset to computed flux values
+        domain.timestep is set to the largest step satisfying all volumes.
+        """
+
+        import sys
+        from Numeric import zeros, Float
+        from config import max_timestep
+
+        N = self.number_of_elements
+
+        neighbours = self.neighbours
+        neighbour_edges = self.neighbour_edges
+        normals = self.normals
+
+        areas = self.areas
+        radii = self.radii
+        edgelengths = self.edgelengths
+
+        timestep = max_timestep #FIXME: Get rid of this
+
+        #Shortcuts
+        Stage = self.quantities['stage']
+
+        #Arrays
+        stage = Stage.edge_values
+
+        stage_bdry = Stage.boundary_values
+
+        flux = zeros(1, Float) #Work array for summing up fluxes
+
+        #Loop
+        for k in range(N):
+            optimal_timestep = float(sys.maxint)
+
+            flux[:] = 0.  #Reset work array
+            for i in range(3):
+                #Quantities inside volume facing neighbour i
+                ql = stage[k, i]
+
+                #Quantities at neighbour on nearest face
+                n = neighbours[k,i]
+                if n < 0:
+                    m = -n-1 #Convert neg flag to index
+                    qr = stage_bdry[m]
+                else:
+                    m = neighbour_edges[k,i]
+                    qr = stage[n, m]
+
+
+                #Outward pointing normal vector
+                normal = normals[k, 2*i:2*i+2]
+
+                #Flux computation using provided function
+                edgeflux, max_speed = self.flux_function(normal, ql, qr)
+                flux -= edgeflux * edgelengths[k,i]
+
+                #Update optimal_timestep
+                try:
+                    optimal_timestep = min(optimal_timestep, radii[k]/max_speed)
+                except ZeroDivisionError:
+                    pass
+
+            #Normalise by area and store for when all conserved
+            #quantities get updated
+            flux /= areas[k]
+            Stage.explicit_update[k] = flux[0]
+
+            timestep = min(timestep, optimal_timestep)
+
+        self.timestep = timestep
+
+
+
+    def evolve(self, yieldstep = None, finaltime = None):
+        """Specialisation of basic evolve method from parent class
+        """
+
+        #Initialise real time viz if requested
+        if self.visualise is True and self.time == 0.0:
+            import realtime_visualisation_new as visualise
+            visualise.create_surface(self)
+
+        #Call basic machinery from parent class
+        for t in Generic_domain.evolve(self, yieldstep, finaltime):
+            #Real time viz
+            if self.visualise is True:
+                visualise.update(self)
+
+            #Pass control on to outer loop for more specific actions
+            yield(t)
Index: /inundation/ga/storm_surge/parallel/pytiming.py
===================================================================
--- /inundation/ga/storm_surge/parallel/pytiming.py	(revision 1363)
+++ /inundation/ga/storm_surge/parallel/pytiming.py	(revision 1363)
@@ -0,0 +1,235 @@
+#!/usr/bin/env python
+
+# Timing of MPI module for Python and estimation of latency and bandwidth
+# 
+# Send numerical array in a ring from processor 0 to 1 etc back to 0 
+# Perform timings and compare different sending strategies 
+#
+# OMN, OCT 2001
+
+import time, sys, pypar, Numeric
+
+# The send/recv routines
+#
+import pypar
+
+#from mpiext import send_array, receive_array  #most direct methods
+#Not available in final install - use bypass form for now
+
+#--------------------------------------------------------------
+
+method = 2     # Use
+                  # 0: automatically allocated buffers
+                  # 1: user-supplied buffers
+                  # 2: Use bypass - let pypar use direct mpi_ext call
+                  # 3: direct call to mpiext (with buffers),
+                  #    only fractionally better than bypass
+                  
+vanilla = False  # Force use of vanilla format (1)
+consistency_check = False # Check correctness
+
+
+#--------------------------------------------------------------
+# linfit
+#
+def linfit(x, y):
+  """Fit a and b to the model y = ax + b. Return a,b,variance
+  """
+  
+  Sx = Sy = SSoN = SxoN = norm = varest = 0.0
+  N = len(x)
+  assert len(y) == N, "x and y must have same length"
+  
+  for i in range(N):
+    #print("x,y = %f, %f\n",x[i],y[i])
+    Sx  = Sx + x[i]
+    Sy  = Sy + y[i]
+  
+  SxoN = Sx/N
+  
+  a = 0.0 
+  for i in range(N):
+    t    = x[i] - SxoN
+    SSoN = SSoN + t*t
+    a    = a + t*y[i]
+
+
+  a = a/SSoN            # a = (N Sxy - SxSy)/(NSxx - Sx^2) */
+  b = (Sy - Sx*a)/N
+  
+  # Quality - variance estimate \sum_i r_i^2 /(m-n) 
+  for i in range(N):  
+    norm = norm + float(x[i])*x[i]
+    res = y[i] - a*x[i] - b
+    varest = varest + res*res
+
+  varest = varest/norm/(N-2)
+  return a, b, varest
+
+
+
+
+#--------------------------------------------------------------
+# Main program
+#
+MAXI  = 10         # Number of blocks 
+MAXM  = 500000     # Largest block 
+BLOCK = MAXM/MAXI  # Block size 
+
+repeats = 10
+msgid = 0
+vanilla = 0 #Select vanilla mode (slower but general)
+
+numprocs = pypar.size()
+myid = pypar.rank()
+processor_name = pypar.Get_processor_name()
+
+if myid == 0:
+  # Main process - Create message, pass on, verify correctness and log timing
+  #
+  print "MAXM = %d, number of processors = %d" %(MAXM, numprocs)
+  print "Measurements are repeated %d times for reliability" %repeats
+
+if numprocs < 2:
+  print "Program needs at least two processors - aborting\n"
+  pypar.Abort()
+   
+pypar.Barrier() #Synchronize all before timing   
+print "I am process %d on %s" %(myid,processor_name)
+
+
+#Initialise data and timings
+#
+from RandomArray import uniform, seed
+seed(17, 53)
+A = uniform(0.0,100.0,MAXM)
+elsize = A.itemsize()
+#print elsize
+
+noelem  = [0]*MAXI
+bytes   = [0]*MAXI         
+avgtime = [0.0]*MAXI         
+mintime = [ 1000000.0]*MAXI      
+maxtime = [-1000000.0]*MAXI            
+
+
+
+
+if myid == 0:   
+  # Determine timer overhead 
+  cpuOH = 1.0;
+  for k in range(repeats):   # Repeat to get reliable timings 
+    t1 = pypar.Wtime()
+    t2 = pypar.Wtime()
+    if t2-t1 < cpuOH: cpuOH = t2-t1
+    
+  print "Timing overhead is %f seconds.\n" %cpuOH         
+
+     
+# Pass msg circularly   
+for k in range(repeats):
+  if myid == 0:
+    print "Run %d of %d" %(k+1,repeats)
+    
+  for i in range(MAXI):
+    m=BLOCK*i+1       
+   
+    noelem[i] = m
+   
+    pypar.Barrier() # Synchronize 
+   
+    if myid == 0:
+      #
+      # Main process
+      #
+      t1 = pypar.Wtime()
+
+      if method==0:
+        pypar.send(A[:m], destination=1, tag=msgid, vanilla=vanilla)
+        C = pypar.receive(numprocs-1, tag=msgid, vanilla=vanilla)
+      elif method == 1:  
+        pypar.send(A[:m], use_buffer=True, destination=1,
+                   tag=msgid, vanilla=vanilla)
+        C = pypar.receive(numprocs-1, buffer=A[:m], tag=msgid, vanilla=vanilla)
+      elif method==2:
+        pypar.send(A[:m], use_buffer=True, destination=1,
+                   tag=msgid, vanilla=vanilla, bypass=True)
+        C = pypar.receive(numprocs-1, buffer=A[:m], tag=msgid,
+                          vanilla=vanilla, bypass=True)
+      else:
+        raise 'Unknown method'
+        #send_array(A[:m], 1, msgid)    
+        #stat = receive_array(A[:m], numprocs-1, msgid)
+        #C = A[:m]
+        
+      t2 = pypar.Wtime() - t1 - cpuOH
+      t2 = t2/numprocs
+      avgtime[i] = avgtime[i] + t2
+      if t2 < mintime[i]: mintime[i] = t2
+      if t2 > maxtime[i]: maxtime[i] = t2
+
+      # Uncomment to verify integrity of data
+      # However, this may affect accuracy of timings for some reason.
+      #
+      if consistency_check:
+        assert Numeric.alltrue(C == A[:m])
+    else:
+      #
+      # Parallel process - get msg and pass it on
+      #
+
+      if method==0:
+        C = pypar.receive(myid-1, tag=msgid, vanilla=vanilla)
+        pypar.send(A[:m], destination=(myid+1)%numprocs,
+                   tag=msgid, vanilla=vanilla)
+      elif method==1:  
+        C = pypar.receive(myid-1, buffer=A[:m], tag=msgid, vanilla=vanilla)
+        pypar.send(A[:m], use_buffer=True, destination=(myid+1)%numprocs,
+                   tag=msgid, vanilla=vanilla)                        
+      elif method==2:
+        # Use pypar bypass
+        C = pypar.receive(myid-1, buffer=A[:m], tag=msgid,
+                          vanilla=vanilla, bypass=True)
+        pypar.send(A[:m], use_buffer=True, destination=(myid+1)%numprocs,
+                   tag=msgid, vanilla=vanilla, bypass=True)
+      else:
+        raise 'Unknown method'
+        # Use direct mpiext call
+        #stat = receive_array(A[:m], myid-1, msgid)                
+        #send_array(A[:m], (myid+1)%numprocs, msgid)
+
+
+# Output stats
+#
+if myid == 0:
+  print "Bytes transferred   time (micro seconds)"
+  print "                    min        avg        max "     
+  print "----------------------------------------------"
+     
+  for i in range(MAXI):
+    avgtime[i] = avgtime[i]/repeats*1.0e6 #Average micro seconds
+    mintime[i] = mintime[i]*1.0e6         #Min micro seconds       
+    maxtime[i] = maxtime[i]*1.0e6         #Min micro seconds              
+          
+    m = noelem[i]
+    bytes[i] = elsize*noelem[i]       
+      
+    print "%10d    %10d %10d %10d" %(bytes[i], mintime[i], avgtime[i], maxtime[i]) 
+
+
+  Tbw, Tlat, varest = linfit(bytes, mintime)
+  print "\nLinear regression on best timings (t = t_l + t_b * bytes):\n",
+  print "  t_b = %f\n  t_l = %f" %(Tbw,Tlat)
+  print "  Estimated relative variance = %.9f\n" %varest
+       
+  print "Estimated bandwith (1/t_b):  %.3f Mb/s" %(1.0/Tbw)   
+  print "Estimated latency:           %d micro s" %int(mintime[0]-bytes[0]*Tbw)  
+
+   
+pypar.Finalize()
+
+
+
+
+  
+  
Index: /inundation/ga/storm_surge/parallel/ring_example.py
===================================================================
--- /inundation/ga/storm_surge/parallel/ring_example.py	(revision 1363)
+++ /inundation/ga/storm_surge/parallel/ring_example.py	(revision 1363)
@@ -0,0 +1,27 @@
+import pypar                                       # Import module and initialise MPI 
+
+proc = pypar.size()                                # Number of processes as specified by mpirun
+myid = pypar.rank()                                # Id of of this process (myid in [0, proc-1]) 
+node = pypar.get_processor_name()                  # Host name on which current process is running
+
+print 'I am proc %d of %d on node %s' %(myid, proc, node)
+
+if myid == 0:                                      # Actions for process 0:
+  msg = 'P0'  
+  pypar.send(msg, destination=1)                   # Send message to proces 1 (right hand neighbour)
+  msg = pypar.receive(source=proc-1)               # Receive message from last process
+      
+  print 'Processor 0 received message "%s" from processor %d' %(msg, proc-1)
+
+else:                                              # Actions for all other processes:
+
+  source = myid-1                                  # Source is the process to the left
+  destination = (myid+1)%proc                      # Destination is process to the right
+                                                   # wrapped so that last processor will 
+						   # send back to proces 0  
+  
+  msg = pypar.receive(source)                      # Receive message from source 
+  msg = msg + 'P' + str(myid)                      # Update message     
+  pypar.send(msg, destination)                     # Send message to destination   
+
+pypar.finalize()                                   # Stop MPI 
Index: /inundation/ga/storm_surge/parallel/run_advection.py
===================================================================
--- /inundation/ga/storm_surge/parallel/run_advection.py	(revision 1362)
+++ /inundation/ga/storm_surge/parallel/run_advection.py	(revision 1363)
@@ -19,9 +19,11 @@
 
 domain.visualise = True
-domain.range_z = 1.0
+#domain.visualise_range_z = 1.0
+domain.visualise_color_stage = True
+
 
 #Boundaries
 T = Transmissive_boundary(domain)
-D = Dirichlet_boundary(array([0.5]))
+D = Dirichlet_boundary(array([1.0]))
 
 domain.default_order = 2
Index: /inundation/ga/storm_surge/parallel/run_parallel_advection.py
===================================================================
--- /inundation/ga/storm_surge/parallel/run_parallel_advection.py	(revision 1363)
+++ /inundation/ga/storm_surge/parallel/run_parallel_advection.py	(revision 1363)
@@ -0,0 +1,36 @@
+import sys
+from os import sep
+sys.path.append('..'+sep+'pyvolution')
+
+#========================================================================
+from config import g, epsilon
+from Numeric import allclose, array, zeros, ones, Float
+from advection import *
+from Numeric import array
+
+from mesh_factory import rectangular
+
+points, vertices, boundary = rectangular(60, 60)
+
+#Create advection domain with direction (1,-1)
+domain = Domain(points, vertices, boundary, velocity=[1.0, 0.0])
+
+# Initial condition is zero by default
+
+domain.visualise = True
+#domain.visualise_range_z = 1.0
+
+
+#Boundaries
+T = Transmissive_boundary(domain)
+D = Dirichlet_boundary(array([1.0]))
+
+domain.default_order = 2
+
+
+domain.set_boundary( {'left': D, 'right': T, 'bottom': T, 'top': T} )
+domain.check_integrity()
+
+#Check that the boundary value gets propagated to all elements
+for t in domain.evolve(yieldstep = 0.01, finaltime = 1.5):
+    domain.write_time()
Index: /inundation/ga/storm_surge/parallel/runpytiming.py
===================================================================
--- /inundation/ga/storm_surge/parallel/runpytiming.py	(revision 1363)
+++ /inundation/ga/storm_surge/parallel/runpytiming.py	(revision 1363)
@@ -0,0 +1,20 @@
+#!/bin/csh
+#PBS -wd  
+#PBS -q normal
+#PBS -l ncpus=4
+#PBS -l walltime=00:10:00,vmem=100Mb
+#PBS -l software=python
+
+#########################################################
+#   
+#  Run pytiming in parallel using Batch Queue on Alpha Server
+#   
+#  Author:  Ole Nielsen,  SMS, ANU, Jan. 2002. 
+#
+#  To execute:
+#    qsub runpytiming
+#########################################################
+
+prun pytiming
+
+
Index: /inundation/ga/storm_surge/parallel/small.tsh
===================================================================
--- /inundation/ga/storm_surge/parallel/small.tsh	(revision 1363)
+++ /inundation/ga/storm_surge/parallel/small.tsh	(revision 1363)
@@ -0,0 +1,163 @@
+31 0 # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Triangulation Vertices...
+0 554.0 -242.0 
+1 443.0 -190.0 
+2 340.0 -184.0 
+3 276.0 -237.0 
+4 281.0 -347.0 
+5 282.0 -380.0 
+6 312.0 -425.0 
+7 387.0 -472.0 
+8 465.0 -487.0 
+9 549.0 -469.0 
+10 566.0 -413.0 
+11 589.0 -359.0 
+12 472.0 -258.0 
+13 377.0 -411.0 
+14 428.0 -257.0 
+15 386.0 -278.0 
+16 372.0 -325.0 
+17 369.0 -371.0 
+18 433.0 -422.0 
+19 476.0 -361.0 
+20 502.0 -304.0 
+21 333.0 -225.0 
+22 305.0 -244.0 
+23 333.0 -270.0 
+24 369.0 -237.0 
+25 278.5 -292.0 
+26 499.90712817 -423.516449623 
+27 532.621862616 -366.885237781 
+28 498.5 -216.0 
+29 322.359287054 -318.872505543 
+30 571.5 -300.5 
+# attribute column titles ...Triangulation Vertex Titles...
+38 # <# of triangles>, next lines <triangle #> [<vertex #>] [<neigbouring triangle #>] [attribute of region] ...Triangulation Triangles...
+0 29 5 17 1 5 26  
+1 17 5 6 -1 3 0  
+2 10 11 27 28 21 -1  
+3 13 17 6 1 6 -1  
+4 15 23 16 24 -1 16  
+5 16 29 17 0 -1 24  
+6 6 7 13 37 3 -1  
+7 29 23 25 9 8 24  
+8 25 4 29 26 7 -1  
+9 25 23 22 10 18 7  
+10 22 23 21 13 11 9  
+11 3 22 21 10 15 18  
+12 21 24 2 14 15 13  
+13 24 21 23 10 16 12  
+14 2 24 1 27 -1 12  
+15 21 2 3 -1 11 12  
+16 24 23 15 4 30 13  
+17 27 30 20 31 36 28  
+18 22 3 25 -1 9 11  
+19 7 8 18 35 37 -1  
+20 8 9 26 23 35 -1  
+21 26 10 27 2 25 23  
+22 26 19 18 -1 35 25  
+23 10 26 9 20 -1 21  
+24 16 23 29 7 5 4  
+25 26 27 19 36 22 21  
+26 29 4 5 -1 0 8  
+27 14 1 24 14 30 29  
+28 27 11 30 -1 17 2  
+29 1 14 12 -1 33 27  
+30 14 24 15 16 -1 27  
+31 20 30 0 -1 34 17  
+32 12 0 28 -1 33 34  
+33 12 28 1 -1 29 32  
+34 0 12 20 -1 31 32  
+35 8 26 18 22 19 20  
+36 27 20 19 -1 25 17  
+37 18 13 7 6 19 -1  
+28 # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Triangulation Segments...
+0 14 12 exterior
+1 14 15 exterior
+2 16 15 exterior
+3 17 16 exterior
+4 13 17 exterior
+5 13 18 exterior
+6 19 18 exterior
+7 19 20 exterior
+8 12 20 exterior
+9 21 22 
+10 23 22 
+11 23 24 
+12 21 24 
+13 4 25 exterior
+14 7 6 exterior
+15 6 5 exterior
+16 3 2 exterior
+17 5 4 exterior
+18 8 7 exterior
+19 9 10 exterior
+20 9 8 exterior
+21 1 28 exterior
+22 2 1 exterior
+23 0 30 exterior
+24 11 10 exterior
+25 25 3 exterior
+26 28 0 exterior
+27 30 11 exterior
+25 0 # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Mesh Vertices...
+0 554.0 -242.0 
+1 443.0 -190.0 
+2 340.0 -184.0 
+3 276.0 -237.0 
+4 281.0 -347.0 
+5 282.0 -380.0 
+6 312.0 -425.0 
+7 387.0 -472.0 
+8 465.0 -487.0 
+9 549.0 -469.0 
+10 566.0 -413.0 
+11 589.0 -359.0 
+12 472.0 -258.0 
+13 377.0 -411.0 
+14 428.0 -257.0 
+15 386.0 -278.0 
+16 372.0 -325.0 
+17 369.0 -371.0 
+18 433.0 -422.0 
+19 476.0 -361.0 
+20 502.0 -304.0 
+21 333.0 -225.0 
+22 305.0 -244.0 
+23 333.0 -270.0 
+24 369.0 -237.0 
+25 # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Mesh Segments...
+0 14 12 
+1 15 14 
+2 16 15 
+3 17 16 
+4 13 17 
+5 18 13 
+6 19 18 
+7 20 19 
+8 12 20 
+9 22 21 
+10 23 22 
+11 24 23 
+12 21 24 
+13 3 4 
+14 6 7 
+15 5 6 
+16 2 3 
+17 4 5 
+18 7 8 
+19 9 10 
+20 8 9 
+21 0 1 
+22 1 2 
+23 11 0 
+24 10 11 
+1 # <# of holes>, next lines <Hole #> <x> <y> ...Mesh Holes...
+0 418.0 -332.0
+1 # <# of regions>, next lines <Region #> <x> <y> <tag>...Mesh Regions...
+0 334.0 -243.0 
+1 # <# of regions>, next lines <Region #> [Max Area] ...Mesh Regions...
+0 
+# Geo reference
+56
+0
+0
Index: /inundation/ga/storm_surge/parallel/test_init.py
===================================================================
--- /inundation/ga/storm_surge/parallel/test_init.py	(revision 1363)
+++ /inundation/ga/storm_surge/parallel/test_init.py	(revision 1363)
@@ -0,0 +1,25 @@
+#!/usr/bin/env python
+# Inital test of MPI module 'pypar' for Python
+# 
+# Run as 
+#   python test_init.py
+# or 
+#   mpirun -np 2 test_init.py
+#  (perhaps try number of processors more than 2)
+#
+# To verify bandwidth of your architecture please run pytiming (and ctiming) 
+#
+# OMN April 2002
+
+
+import pypar
+
+myid =    pypar.rank()
+numproc = pypar.size()
+node =    pypar.Get_processor_name()
+
+print "I am processor %d of %d on node %s" %(myid, numproc, node)
+if myid == 0:
+    print 'Maximal tag is %d' %pypar.max_tag
+pypar.Finalize()
+
Index: /inundation/ga/storm_surge/parallel/testpypar.py
===================================================================
--- /inundation/ga/storm_surge/parallel/testpypar.py	(revision 1363)
+++ /inundation/ga/storm_surge/parallel/testpypar.py	(revision 1363)
@@ -0,0 +1,1102 @@
+#!/usr/bin/env python
+# Test of MPI module 'pypar' for Python
+# 
+# Run as 
+#   python testpypar.py
+# or 
+#   mpirun -np 2 testpypar.py
+# (perhaps try number of processors more than 2)
+#
+# To verify bandwidth of your architecture please run pytiming (and ctiming) 
+#
+# OMN, GPC FEB 2002
+
+try:
+  import Numeric
+except:
+  raise 'Module Numeric must be present to run pypar'
+
+
+#print "Importing pypar"
+import pypar
+methods = dir(pypar)
+assert 'Abort' in methods
+assert 'Finalize' in methods
+assert 'Get_processor_name' in methods
+assert 'Wtime' in methods
+assert 'rank' in methods
+assert 'raw_receive' in methods
+assert 'raw_send' in methods
+assert 'receive' in methods
+assert 'send' in methods
+assert 'bcast' in methods
+assert 'size' in methods
+
+#print "Module pypar imported OK"
+#pypar.Barrier()
+
+
+myid =    pypar.rank()
+numproc = pypar.size()
+node =    pypar.Get_processor_name()
+
+print "I am processor %d of %d on node %s" %(myid, numproc, node)
+pypar.Barrier()
+
+
+if numproc > 1:
+  # Test simple raw communication (arrays, strings and general)
+  #
+  N = 17 #Number of elements
+  
+  if myid == 0:
+    # Integer arrays
+    #
+    A = Numeric.array(range(N))
+    B = Numeric.zeros(N)    
+    pypar.raw_send(A,1)
+    pypar.raw_receive(B,numproc-1)
+    
+    assert Numeric.allclose(A, B)
+    print "Raw communication of numeric integer arrays OK"
+
+    # Long integer arrays
+    #
+    A = Numeric.array(range(N)).astype('l')
+    B = Numeric.zeros(N)    
+    pypar.raw_send(A,1)
+    pypar.raw_receive(B,numproc-1)
+    
+    assert Numeric.allclose(A, B)
+    print "Raw communication of numeric long integer arrays OK"
+
+    # Real arrays
+    #
+    A = Numeric.array(range(N)).astype('f')
+    B = Numeric.zeros(N).astype('f')    
+    pypar.raw_send(A,1)
+    pypar.raw_receive(B,numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Raw communication of numeric real arrays OK"
+
+    # Complex arrays
+    #
+    A = Numeric.array(range(N)).astype('D')
+    A += 3j
+    B = Numeric.zeros(N).astype('D')
+    pypar.raw_send(A, 1)
+    B = pypar.raw_receive(B, numproc-1)
+
+    assert Numeric.allclose(A, B)    
+    print "Raw communication of numeric complex arrays OK"
+
+    # Strings (< 256 characters)
+    #
+    A = "and now to something completely different !"
+    B = " "*len(A)
+    pypar.raw_send(A,1)
+    B, status = pypar.raw_receive(B,numproc-1,return_status=True)
+
+    assert A == B
+    print "Raw communication of strings OK"
+    
+    # A more general structure
+    #
+    A = ['ABC', (1,2,3.14), {8: 'Monty'}, Numeric.array([13.45, 1.2])]
+    B = ['   ', (0,0,0.0), {0: '     '}, Numeric.zeros(2).astype('f')]    
+    pypar.raw_send(A,1)
+    B, status = pypar.raw_receive(B,numproc-1, return_status=True)
+
+    assert A == B
+    print "Raw communication of general structures OK"
+    
+  else:  
+    # Integers
+    #
+    X = Numeric.zeros(N)
+    pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)
+  
+    # Long integers
+    #
+    X = Numeric.zeros(N).astype('l')
+    pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)    
+
+    # Floats
+    #
+    X = Numeric.zeros(N).astype('f')
+    pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)    
+
+    # Complex
+    #
+    X = Numeric.zeros(N).astype('D')
+    X = pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)    
+
+    # Strings
+    #
+    X = " "*256
+    pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X.strip(), (myid+1)%numproc)    
+
+    # General
+    #
+    X = ['   ', (0,0,0.0), {0: '     '}, Numeric.zeros(2).astype('f')]
+    X = pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)    
+    
+
+
+
+
+  #Test (raw communication of) multi dimensional arrays
+
+  M = 13  #Number of elements in dim 1
+  N = 17  #Number of elements in higher dims
+  
+  if myid == 0:
+    # 2D real arrays
+    #
+    A = Numeric.array(range(M*N)).astype('f')
+    B = Numeric.zeros(M*N).astype('f')
+
+    A = Numeric.reshape(A, (M,N))
+    B = Numeric.reshape(B, (M,N))    
+    
+    pypar.raw_send(A,1)
+    pypar.raw_receive(B,numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Raw communication of 2D real arrays OK"
+
+    # 2D complex arrays
+    #
+    A = Numeric.array(range(M*N)).astype('D')
+    B = Numeric.zeros(M*N).astype('D')
+
+    A = Numeric.reshape(A, (M,N))
+    B = Numeric.reshape(B, (M,N))    
+
+    pypar.raw_send(A,1)
+    pypar.raw_receive(B,numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Raw communication of 2D complex arrays OK"
+
+    # 3D real arrays
+    #
+    A = Numeric.array(range(M*N*N)).astype('f')
+    B = Numeric.zeros(M*N*N).astype('f')
+
+    A = Numeric.reshape(A, (M,N,N))
+    B = Numeric.reshape(B, (M,N,N))    
+    
+    pypar.raw_send(A,1)
+    pypar.raw_receive(B,numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Raw communication of 3D real real arrays OK"
+
+    # 4D real arrays
+    #
+    A = Numeric.array(range(M*N*N*M)).astype('f')
+    B = Numeric.zeros(M*N*N*M).astype('f')
+
+    A = Numeric.reshape(A, (M,N,N,M))
+    B = Numeric.reshape(B, (M,N,N,M))    
+    
+    pypar.raw_send(A,1)
+    pypar.raw_receive(B,numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Raw communication of 4D real real arrays OK"
+
+    # 5D real arrays
+    #
+    A = Numeric.array(range(M*N*2*N*M)).astype('f')
+    B = Numeric.zeros(M*N*2*N*M).astype('f')
+
+    A = Numeric.reshape(A, (M,N,2,N,M))
+    B = Numeric.reshape(B, (M,N,2,N,M))    
+    
+    pypar.raw_send(A,1)
+    pypar.raw_receive(B,numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Raw communication of 5D real real arrays OK"
+
+    # 5D double arrays
+    #
+    A = Numeric.array(range(M*N*2*N*M)).astype('d')
+    B = Numeric.zeros(M*N*2*N*M).astype('d')
+
+    A = Numeric.reshape(A, (M,N,2,N,M))
+    B = Numeric.reshape(B, (M,N,2,N,M))    
+    
+    pypar.raw_send(A,1)
+    pypar.raw_receive(B,numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Raw communication of 5D double arrays OK"
+
+    # 5D complex arrays
+    #
+    A = Numeric.array(range(M*N*2*N*M)).astype('D')
+    B = Numeric.zeros(M*N*2*N*M).astype('D')
+
+    A = Numeric.reshape(A, (M,N,2,N,M))
+    B = Numeric.reshape(B, (M,N,2,N,M))    
+    
+    pypar.raw_send(A,1)
+    pypar.raw_receive(B,numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Raw communication of 5D complex arrays OK"
+  else:  
+    # 2D real arrays
+    #
+    X = Numeric.zeros(M*N).astype('f')
+    X = Numeric.reshape(X, (M,N))
+    
+    pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)
+  
+    # 2D complex arrays
+    #
+    X = Numeric.zeros(M*N).astype('D')
+    X = Numeric.reshape(X, (M,N))
+
+    X = pypar.raw_receive(X, myid-1)
+    pypar.raw_send(X, (myid+1)%numproc)
+  
+    # 3D real arrays
+    #
+    X = Numeric.zeros(M*N*N).astype('f')
+    X = Numeric.reshape(X, (M,N,N))
+    
+    pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)
+
+    # 4D real arrays
+    #
+    X = Numeric.zeros(M*N*N*M).astype('f')
+    X = Numeric.reshape(X, (M,N,N,M))
+    
+    pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)
+
+    # 5D real arrays
+    #
+    X = Numeric.zeros(M*N*2*N*M).astype('f')
+    X = Numeric.reshape(X, (M,N,2,N,M))
+    
+    pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)
+
+    # 5D double arrays
+    #
+    X = Numeric.zeros(M*N*2*N*M).astype('d')
+    X = Numeric.reshape(X, (M,N,2,N,M))
+    
+    pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)
+
+    # 5D complex arrays
+    #
+    X = Numeric.zeros(M*N*2*N*M).astype('D')
+    X = Numeric.reshape(X, (M,N,2,N,M))
+    
+    pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)
+
+  # Test easy communication  - without buffers (arrays, strings and general)
+  #
+  N = 17 #Number of elements
+  
+  if myid == 0:
+    # Integer arrays
+    #
+    A = Numeric.array(range(N))
+
+    pypar.send(A,1)
+    B = pypar.receive(numproc-1)
+    
+
+    assert Numeric.allclose(A, B)
+    print "Simplified communication of numeric integer arrays OK"
+
+    # Long integer arrays
+    #
+    A = Numeric.array(range(N)).astype('l')
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Simplified communication of long integer real arrays OK"
+
+    # Real arrays
+    #
+    A = Numeric.array(range(N)).astype('f')
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Simplified communication of numeric real arrays OK"
+
+    # Complex arrays
+    #
+    A = Numeric.array(range(N)).astype('D')
+    A += 3j
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+
+    assert Numeric.allclose(A, B)    
+    print "Simplified communication of numeric complex arrays OK"
+
+
+    # Strings
+    #
+    A = "and now to something completely different !"
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert A == B
+    print "Simplified communication of strings OK"
+    
+    # A more general structure
+    #
+    A = ['ABC', (1,2,3.14), {8: 'Monty'}, Numeric.array([13.45, 1.2])]
+    pypar.send(A,1)
+    B = pypar.receive(numproc-1)
+    
+    assert A == B
+    print "Simplified communication of general structures OK"
+    
+  else:  
+    # Integers
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+  
+    # Long integers
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+    # Floats
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+    # Complex
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+    # Strings
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+    # General
+    #
+    X = pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+
+
+
+    
+  #Test (easy communication of) multi dimensional arrays
+
+  M = 13  #Number of elements in dim 1
+  N = 17  #Number of elements in higher dims
+  
+  if myid == 0:
+    # 2D real arrays
+    #
+    A = Numeric.array(range(M*N)).astype('f')
+
+    A = Numeric.reshape(A, (M,N))
+    
+    pypar.send(A,1)
+    B = pypar.receive(numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Simplified communication of 2D real arrays OK"
+
+    # 3D real arrays
+    #
+    A = Numeric.array(range(M*N*N)).astype('f')
+    A = Numeric.reshape(A, (M,N,N))
+    
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Simplified communication of 3D real arrays OK"
+
+    # 4D real arrays
+    #
+    A = Numeric.array(range(M*N*N*M)).astype('f')
+    A = Numeric.reshape(A, (M,N,N,M))
+    
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Simplified communication of 4D real arrays OK"
+
+    # 5D real arrays
+    #
+    A = Numeric.array(range(M*N*2*N*M)).astype('f')
+    A = Numeric.reshape(A, (M,N,2,N,M))
+    
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Simplified communication of 5D real arrays OK"
+
+    # 5D double arrays
+    #
+    A = Numeric.array(range(M*N*2*N*M)).astype('d')
+    A = Numeric.reshape(A, (M,N,2,N,M))
+    
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Simplified communication of 5D double arrays OK"
+
+    # 5D complex arrays
+    #
+    A = Numeric.array(range(M*N*2*N*M)).astype('D')
+    A = Numeric.reshape(A, (M,N,2,N,M))
+    
+    pypar.send(A,1)
+    B=pypar.receive(numproc-1)
+    
+    assert Numeric.allclose(A, B)    
+    print "Simplified communication of 5D complex real arrays OK"
+
+  else:  
+    # 2D real arrays
+    #
+    
+    X = pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+  
+    # 3D real arrays
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+
+    # 4D real arrays
+    #
+    X = pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+
+    # 5D real arrays
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+
+    # 5D double arrays
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+
+    # 5D complex arrays
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+
+  # Test broadcast  - with buffers (arrays, strings and general)
+  #
+      
+  testString = ('test' + str(myid)).ljust(10)  #Buffers must have the same length on all procs!
+  pypar.bcast(testString, 0)
+  assert testString.strip() == 'test0'
+  
+  testString = ('test' + str(myid)).ljust(10)  #Buffers must have the same length on all procs!
+  pypar.bcast(testString, numproc-1)
+  assert testString.strip() == 'test' + str(numproc-1)
+  
+  if myid == 0:
+    print "Broadcast communication of strings OK"
+
+  ####################################################  
+  N = 17 #Number of elements
+  testArray = myid * Numeric.array(range(N))
+  pypar.bcast(testArray, 1)
+  assert Numeric.allclose(testArray, 1 * testArray)
+  
+  if myid == 0:    
+    print "Broadcast communication of numeric integer array OK"
+
+
+  testArray = myid * Numeric.array(range(N)).astype('f')
+  pypar.bcast(testArray, 1)
+  assert Numeric.allclose(testArray, 1 * testArray)
+  if myid == 0:
+    print "Broadcast communication of numeric real array OK"
+
+
+  M = 13
+  testArray = myid * Numeric.array(range(M*N)).astype('f')
+  testArray = Numeric.reshape(testArray, (M,N))  
+  pypar.bcast(testArray, 1)
+  assert Numeric.allclose(testArray, 1 * testArray)
+  if myid == 0:
+    print "Broadcast communication of 2D numeric real array OK"
+
+  testArray = myid * Numeric.array(range(M*2*N)).astype('f')
+  testArray = Numeric.reshape(testArray, (M,2,N))  
+  pypar.bcast(testArray, 1)
+  assert Numeric.allclose(testArray, 1 * testArray)
+  if myid == 0:
+    print "Broadcast communication of 3D numeric real array OK"
+
+  testArray = myid * Numeric.array(range(M*2*N)).astype('D')
+  testArray = Numeric.reshape(testArray, (M,2,N))  
+  pypar.bcast(testArray, 1)
+  assert Numeric.allclose(testArray, 1 * testArray)
+  if myid == 0:
+    print "Broadcast communication of 3D numeric complex array OK"
+
+    
+  testGeneral = ['ABC', myid, (1,2,3), {8: 'Monty'}, Numeric.array([13.45, 1.2])]
+  
+  testGeneral = pypar.bcast(testGeneral, 1)
+  
+  assert testGeneral ==  ['ABC', 1, (1,2,3), {8: 'Monty'}, Numeric.array([13.45, 1.2])], 'Should have been %s' %testGeneral
+  
+  if myid == 0:
+    print "Broadcast communication of general structures OK"
+  
+
+
+
+  # Test scatter  - with/without buffers (arrays, strings)
+  #
+  N = 16 #Number of elements
+
+  NP = N/numproc
+  
+  testString = 'ABCDEFGHIJKLMNOP'  #Length = 16
+  X = ' '*NP
+  pypar.raw_scatter(testString, X, 1)
+  
+  Y = pypar.scatter(testString, 1)
+
+  #print 'P%d: s=%s, r=%s, %s' %(myid, testString, X, Y)
+
+  assert X==Y, 'X=%s, Y=%s' %(X,Y)
+  assert Y == testString[myid*NP:(myid+1)*NP]
+  assert X == testString[myid*NP:(myid+1)*NP]
+
+  if myid == 0:    
+    print "Scatter communication of strings OK"
+    
+
+  #Scatter Arrays
+  testArray = Numeric.array(range(N))
+  X = Numeric.zeros(NP)
+  pypar.raw_scatter(testArray, X, 0)
+  Y = pypar.scatter(testArray, 0)
+
+  assert Numeric.allclose(X, Y)
+  assert Numeric.allclose(Numeric.array(X), testArray[myid*NP:(myid+1)*NP])
+  assert Numeric.allclose(Y, testArray[myid*NP:(myid+1)*NP])   
+
+  if myid == 0:
+    print "Scatter communication of numeric integer array OK"
+
+
+  testArray = Numeric.array(range(N)).astype('f')
+  X = Numeric.zeros(NP).astype('f')
+  pypar.raw_scatter(testArray, X, 0)
+    
+  Y = pypar.scatter(testArray, 0)
+
+  assert Numeric.allclose(X, Y)  
+  assert Numeric.allclose(X, testArray[myid*NP:(myid+1)*NP])
+  assert Numeric.allclose(Y, testArray[myid*NP:(myid+1)*NP])   
+
+  if myid == 0:
+    print "Scatter communication of numeric real arrays OK"
+  #else:
+  #  print X, testArray, Y
+  #  assert Numeric.allclose(X, Y)
+
+  testArray = Numeric.array(range(N)).astype('D')
+  X = Numeric.zeros(NP).astype('D')
+  pypar.raw_scatter(testArray, X, 0)
+    
+  Y = pypar.scatter(testArray, 0)
+
+  assert Numeric.allclose(X, Y)  
+  assert Numeric.allclose(X, testArray[myid*NP:(myid+1)*NP])
+  assert Numeric.allclose(Y, testArray[myid*NP:(myid+1)*NP])   
+
+  if myid == 0:
+    print "Scatter communication of numeric complex array OK"
+
+
+###################################################
+#FIXME: 2D scatter doesn't work yet    
+#   M = 16
+#   N = 13
+#   MP = M/numproc
+  
+#   testArray = Numeric.array(range(M*N)).astype('D')
+#   testArray = Numeric.reshape(testArray, (M,N))    
+#   X = Numeric.zeros(MP*N).astype('D')
+#   X = Numeric.reshape(X, (MP,N))
+  
+#   pypar.raw_scatter(testArray, X, 0)
+#   Y = pypar.scatter(testArray, 0)
+  
+#   assert Numeric.allclose(X, Y)  
+#   assert Numeric.allclose(X, testArray[myid*MP:(myid+1)*MP,:])
+#   assert Numeric.allclose(Y, testArray[myid*MP:(myid+1)*MP,:])   
+
+#   if myid == 0:
+#     print "Scatter communication of 2D numeric complex array OK"
+
+
+
+  # Test gather  - with/without buffers (arrays, strings)
+  #
+  N = 17 #Number of elements
+      
+  testString = 'AB'
+  X = '_'*(len(testString)*numproc)    #Blanks caused errors when numproc >= 6 
+  pypar.raw_gather(testString, X, 0)
+
+  Y =  pypar.gather(testString, 0) 
+
+
+  if myid == 0:
+    assert X == 'AB' * numproc
+    assert Y == 'AB' * numproc
+    print "Gather communication of strings OK"
+  
+
+  testArray = Numeric.array(range(N))
+  X = Numeric.zeros(N*numproc)
+  pypar.raw_gather(testArray, X, 0)
+
+  Y = pypar.gather(testArray, 0)
+  
+  if myid == 0:
+    for i in range(numproc):       
+      assert Numeric.allclose(testArray, X[(i * N): ((i+1)*N)])
+
+    assert Numeric.allclose(X, Y)
+    print "Gather communication of numeric integer array OK"
+    
+    
+  testArray = Numeric.array(range(N)).astype('f')
+  X = Numeric.zeros(N*numproc).astype('f')
+  pypar.raw_gather(testArray, X, 0)
+  
+  Y = pypar.gather(testArray, 0)
+  if myid == 0:
+    for i in range(numproc):       
+      assert Numeric.allclose(testArray, X[(i * N): ((i+1)*N)])
+    assert Numeric.allclose(X, Y)      
+    print "Gather communication of numeric real array OK"
+    
+  
+  testArray = Numeric.array(range(N)).astype('D')
+  X = Numeric.zeros(N*numproc).astype('D')
+  pypar.raw_gather(testArray, X, 0)
+  
+  Y = pypar.gather(testArray, 0)
+  if myid == 0:
+    for i in range(numproc):       
+      assert Numeric.allclose(testArray, X[(i * N): ((i+1)*N)])
+    assert Numeric.allclose(X, Y)      
+    print "Gather communication of numeric complex arrays OK"
+
+  M = 13  
+  testArray = Numeric.array(range(M*N)).astype('D')
+  testArray = Numeric.reshape(testArray, (M,N))
+  X = Numeric.zeros(M*N*numproc).astype('D')
+  X = Numeric.reshape(X, (M*numproc,N))
+  
+  pypar.raw_gather(testArray, X, 0)
+  
+  Y = pypar.gather(testArray, 0)
+  if myid == 0:
+    for i in range(numproc):       
+      assert Numeric.allclose(testArray, X[(i * M): ((i+1)*M), :])
+    assert Numeric.allclose(X, Y)      
+    print "Gather communication of 2D numeric complex arrays OK"
+    
+  
+  ########################################################
+  # Test reduce
+  #######################################################
+  N = 17 #Number of elements
+  
+  # Create one (different) array on each processor
+  #    
+  testArray = Numeric.array(range(N)) * (myid+1)
+  #print testArray
+  X = Numeric.zeros(N) # Buffer for results
+
+  pypar.raw_reduce(testArray, X, pypar.SUM, 0)
+  if myid == 0:
+    Y = Numeric.zeros(N)
+    for i in range(numproc):
+      Y = Y+Numeric.array(range(N))*(i+1)    
+    #print X
+    #print Y  
+    assert Numeric.allclose(X, Y)
+    print "Raw reduce using pypar.SUM OK"
+        
+  pypar.raw_reduce(testArray, X, pypar.MAX, 0, 0)
+  if myid == 0:
+    Y = Numeric.array(range(N))*numproc
+    assert Numeric.allclose(X, Y)
+    print "Raw reduce using pypar.MAX OK"
+
+  pypar.raw_reduce(testArray, X, pypar.MIN, 0, 0)
+  if myid == 0:
+    Y = Numeric.array(range(N))
+    assert Numeric.allclose(X, Y)
+    print "Raw reduce using pypar.MIN OK"
+    
+  if numproc <= 20:
+    testArray_float = testArray.astype('f')  #Make room for big results
+    X_float = X.astype('f')
+    pypar.raw_reduce(testArray_float, X_float, pypar.PROD, 0, 0)
+    if myid == 0:
+      Y = Numeric.ones(N).astype('f')    
+      for i in range(numproc):
+        Y = Y*Numeric.array(range(N))*(i+1)    
+      #print X_float
+      #print Y  
+      assert Numeric.allclose(X_float, Y)
+      print "Raw reduce using pypar.PROD OK"
+  else:
+    if myid == 0:
+      print "Skipping product-reduce - try again with numproc < 20"    
+
+  pypar.raw_reduce(testArray, X, pypar.LAND, 0, 0)
+  if myid == 0:  
+    Y = Numeric.ones(N)    
+    for i in range(numproc):
+      Y = Numeric.logical_and(Y, Numeric.array(range(N))*(i+1))  
+    assert Numeric.allclose(X, Y)
+    print "Raw reduce using pypar.LAND OK"    
+    
+  pypar.raw_reduce(testArray, X, pypar.BAND, 0, 0)
+  if myid == 0:
+    Y = Numeric.ones(N)*255  #Neutral element for &   
+    for i in range(numproc):
+      Y = Numeric.bitwise_and(Y, Numeric.array(range(N))*(i+1))
+    assert Numeric.allclose(X, Y)
+    print "Raw reduce using pypar.BAND OK"    
+
+  pypar.raw_reduce(testArray, X, pypar.LOR, 0, 0)
+  if myid == 0:  
+    Y = Numeric.zeros(N)    
+    for i in range(numproc):
+      Y = Numeric.logical_or(Y, Numeric.array(range(N))*(i+1))  
+    assert Numeric.allclose(X, Y)
+    print "Raw reduce using pypar.LOR OK"    
+  
+  pypar.raw_reduce(testArray, X, pypar.BOR, 0, 0)
+  if myid == 0:
+    Y = Numeric.zeros(N)   #Neutral element for |   
+    for i in range(numproc):
+      Y = Numeric.bitwise_or(Y, Numeric.array(range(N))*(i+1))
+    assert Numeric.allclose(X, Y)
+    print "Raw reduce using pypar.BOR OK"    
+
+  pypar.raw_reduce(testArray, X, pypar.LXOR, 0, 0)
+  if myid == 0:  
+    Y = Numeric.zeros(N)    
+    for i in range(numproc):
+      Y = Numeric.logical_xor(Y, Numeric.array(range(N))*(i+1))  
+    assert Numeric.allclose(X, Y)
+    print "Raw reduce using pypar.LXOR OK"    
+
+  pypar.raw_reduce(testArray, X, pypar.BXOR, 0, 0)
+  if myid == 0:
+    Y = Numeric.zeros(N)   #Neutral element for xor ?   
+    for i in range(numproc):
+      Y = Numeric.bitwise_xor(Y, Numeric.array(range(N))*(i+1))
+    assert Numeric.allclose(X, Y)
+    print "Raw reduce using pypar.BXOR OK"    
+
+  # NOT YET SUPPORTED
+  #  
+  #pypar.raw_reduce(testArray, X, N, pypar.MAXLOC, 0, 0)  
+  #if myid == 0:
+  #  print 'MAXLOC', X
+  #pypar.raw_reduce(testArray, X, N, pypar.MINLOC, 0, 0)
+  #if myid == 0:
+  #  print 'MINLOC', X
+  
+  #
+  #  FIXME
+  # Don't know how to test this (not available on all MPI systems)
+  #
+  #pypar.raw_reduce(testArray, X, N, pypar.REPLACE, 0, 0)
+  #if myid == 0:
+  #  print 'REPLACE', X
+
+
+
+  # Test status block (simple communication)
+  N = 17  
+  if myid == 0:
+    # Integer arrays
+    #
+    A = Numeric.array(range(N))
+
+    pypar.send(A,1)
+    B, status = pypar.receive(numproc-1, return_status = True)
+
+    #print status
+
+    sz = A.itemsize()
+    assert Numeric.allclose(A, B)
+    assert len(B) == status.length, 'Reported length == %d should be %d'\
+           %(status.length, len(B))
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+
+    print "Status object (numeric integer arrays) OK"
+           
+    # Real arrays
+    #
+    A = Numeric.array(range(N)).astype('f')
+    pypar.send(A,1)
+    B, status = pypar.receive(numproc-1, return_status = True)    
+
+    sz = A.itemsize()    
+    assert Numeric.allclose(A, B)
+    assert len(B) == status.length, 'Reported length == %d should be %d'\
+           %(status.length, len(B))
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+    
+    print "Status object (numeric real arrays) OK"
+
+    # Strings
+    #
+    A = "and now to something completely different !"
+    pypar.send(A,1)
+    B, status = pypar.receive(numproc-1, return_status = True)        
+
+    sz = 1 #Characters are one byte long
+    assert A == B
+    assert len(B) == status.length, 'Reported length == %d should be %d'\
+           %(status.length, len(B))
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+
+    print "Status object (strings) OK"
+    
+    # A more general structure
+    #
+    A = ['ABC', (1,2,3.14), {8: 'Monty'}, Numeric.array([13.45, 1.2])]
+    pypar.send(A,1)
+    B, status = pypar.receive(numproc-1, return_status = True)            
+    
+    assert A == B
+
+    #Length is the number of characters needed to encode the structure
+    #Can't think of a test.
+
+    sz = 1
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+
+    print "Status object (general structures) OK"
+    
+  else:  
+    # Integers
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)
+  
+    # Floats
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+    # Strings
+    #
+    X=pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+    # General
+    #
+    X = pypar.receive(myid-1)  
+    pypar.send(X, (myid+1)%numproc)    
+
+
+
+
+
+
+
+
+  # Test status block (raw communication)
+  N = 17 #Number of elements
+  if myid == 0:
+    # Integer arrays
+    #
+    A = Numeric.array(range(N))
+    B = Numeric.zeros(N)    
+    pypar.raw_send(A,1)
+    B, status = pypar.raw_receive(B,numproc-1,return_status=True)
+    
+    assert Numeric.allclose(A, B)
+
+    sz = A.itemsize()
+    assert Numeric.allclose(A, B)
+    assert len(B) == status.length, 'Reported length == %d should be %d'\
+           %(status.length, len(B))
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+
+    print "Status object (raw numeric integer arrays) OK"
+
+
+    # Real arrays
+    #
+    A = Numeric.array(range(N)).astype('f')
+    B = Numeric.zeros(N).astype('f')    
+    pypar.raw_send(A,1)
+    B, status = pypar.raw_receive(B,numproc-1,return_status=True)    
+    
+    assert Numeric.allclose(A, B)
+    sz = A.itemsize()
+    assert Numeric.allclose(A, B)
+    assert len(B) == status.length, 'Reported length == %d should be %d'\
+           %(status.length, len(B))
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+
+    print "Status object (raw numeric real arrays) OK"
+
+    # Strings (< 256 characters)
+    #
+    A = "and now to something completely different !"
+    B = " "*len(A)
+    pypar.raw_send(A,1)
+    B, status = pypar.raw_receive(B,numproc-1,return_status=True)
+
+
+    sz = 1 #Characters are one byte long
+    assert A == B
+    assert len(B) == status.length, 'Reported length == %d should be %d'\
+           %(status.length, len(B))
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+
+    print "Status object (raw strings) OK"
+    
+
+    
+    # A more general structure
+    #
+    A = ['ABC', (1,2,3.14), {8: 'Monty'}, Numeric.array([13.45, 1.2])]
+    B = ['   ', (0,0,0.0), {0: '     '}, Numeric.zeros(2).astype('f')]    
+    pypar.raw_send(A,1)
+    B, status = pypar.raw_receive(B,numproc-1, return_status=True)
+
+
+    assert A == B
+    sz = 1
+    assert status.size == sz, 'Reported size == %d should be %d'\
+           %(status.size, sz)
+    assert status.tag == pypar.default_tag, 'Reported tag == %d should be %d'\
+           %(status.tag, pypar.default_tag)
+    assert status.error == 0
+    assert status.source == numproc-1, 'Reported source == %d should be %d'\
+           %(status.source, numproc-1)
+
+    
+    print "Status object (raw general structures) OK"
+    
+   
+  else:  
+    # Integers
+    #
+    X = Numeric.zeros(N)
+    pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)
+  
+    # Floats
+    #
+    X = Numeric.zeros(N).astype('f')
+    pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)    
+
+    
+    # Strings
+    #
+    X = " "*256
+    X = pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X.strip(), (myid+1)%numproc)    
+
+    # General
+    #
+    X = ['   ', (0,0,0.0), {0: '     '}, Numeric.zeros(2).astype('f')]
+    X = pypar.raw_receive(X, myid-1)  
+    pypar.raw_send(X, (myid+1)%numproc)    
+    
+
+
+pypar.Finalize()
+
+
+
+
Index: /inundation/ga/storm_surge/pyvolution/advection.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/advection.py	(revision 1362)
+++ /inundation/ga/storm_surge/pyvolution/advection.py	(revision 1363)
@@ -1,2 +1,6 @@
+import sys
+from os import sep
+sys.path.append('..'+sep+'pyvolution')
+
 """Class Domain -
 2D triangular domains for finite-volume computations of
@@ -17,5 +21,5 @@
 
 Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
-Geoscience Australia, 2004   
+Geoscience Australia, 2004
 """
 
@@ -27,5 +31,5 @@
     def __init__(self, coordinates, vertices, boundary = None, velocity = None):
 
-        
+
         Generic_domain.__init__(self, coordinates, vertices, boundary,
                                 ['stage'])
@@ -36,10 +40,14 @@
             self.velocity = velocity
 
-        #Only first is implemented for advection    
-        self.default_order = self.order = 1  
+        #Only first is implemented for advection
+        self.default_order = self.order = 1
 
         #Realtime visualisation
         self.visualise = False
+        self.visualise_color_stage = False
+        self.visualise_timer = True
+        self.visualise_range_z = None
         self.smooth = True
+
 
 
@@ -49,10 +57,10 @@
         msg = 'Conserved quantity must be "stage"'
         assert self.conserved_quantities[0] == 'stage', msg
-        
+
 
     def flux_function(self, normal, ql, qr, zl=None, zr=None):
         """Compute outward flux as inner product between velocity
         vector v=(v_1, v_2) and normal vector n.
-            
+
         if <n,v> > 0 flux direction is outward bound and its magnitude is
         determined by the quantity inside volume: ql.
@@ -60,5 +68,5 @@
         quantity outside the volume: qr.
         """
-    
+
         v1 = self.velocity[0]
         v2 = self.velocity[1]
@@ -71,15 +79,15 @@
         else:
             flux = ql * normal_velocity
-        
+
         max_speed = abs(normal_velocity)
-        return flux, max_speed        
+        return flux, max_speed
 
-        
+
     def compute_fluxes(self):
         """Compute all fluxes and the timestep suitable for all volumes
         in domain.
-        
+
         Compute total flux for each conserved quantity using "flux_function"
-        
+
         Fluxes across each edge are scaled by edgelengths and summed up
         Resulting flux is then scaled by area and stored in
@@ -93,5 +101,5 @@
         Post conditions:
         domain.explicit_update is reset to computed flux values
-        domain.timestep is set to the largest step satisfying all volumes. 
+        domain.timestep is set to the largest step satisfying all volumes.
         """
 
@@ -101,5 +109,5 @@
 
         N = self.number_of_elements
-    
+
         neighbours = self.neighbours
         neighbour_edges = self.neighbour_edges
@@ -109,5 +117,5 @@
         radii = self.radii
         edgelengths = self.edgelengths
-            
+
         timestep = max_timestep #FIXME: Get rid of this
 
@@ -119,5 +127,5 @@
 
         stage_bdry = Stage.boundary_values
-    
+
         flux = zeros(1, Float) #Work array for summing up fluxes
 
@@ -132,14 +140,14 @@
 
                 #Quantities at neighbour on nearest face
-                n = neighbours[k,i] 
+                n = neighbours[k,i]
                 if n < 0:
                     m = -n-1 #Convert neg flag to index
                     qr = stage_bdry[m]
-                else:    
+                else:
                     m = neighbour_edges[k,i]
                     qr = stage[n, m]
 
-                    
-                #Outward pointing normal vector   
+
+                #Outward pointing normal vector
                 normal = normals[k, 2*i:2*i+2]
 
@@ -147,5 +155,5 @@
                 edgeflux, max_speed = self.flux_function(normal, ql, qr)
                 flux -= edgeflux * edgelengths[k,i]
-                
+
                 #Update optimal_timestep
                 try:
@@ -158,9 +166,9 @@
             flux /= areas[k]
             Stage.explicit_update[k] = flux[0]
-            
+
             timestep = min(timestep, optimal_timestep)
 
-        self.timestep = timestep    
-	
+        self.timestep = timestep
+
 
 
@@ -168,8 +176,8 @@
         """Specialisation of basic evolve method from parent class
         """
-        
+
         #Initialise real time viz if requested
         if self.visualise is True and self.time == 0.0:
-            import realtime_visualisation as visualise
+            import realtime_visualisation_new as visualise
             visualise.create_surface(self)
 
@@ -180,5 +188,4 @@
                 visualise.update(self)
 
-            #Pass control on to outer loop for more specific actions    
+            #Pass control on to outer loop for more specific actions
             yield(t)
-        
Index: /inundation/ga/storm_surge/pyvolution/data_manager.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 1362)
+++ /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 1363)
@@ -29,5 +29,5 @@
 Manually created files:
 ASC, PRJ:     Digital elevation models (gridded)
-TSH:          Triangular meshes (e.g. dreated from pmesh)
+TSH:          Triangular meshes (e.g. created from pmesh)
 NC            Model outputs for use as boundary conditions (e.g from MOST)
 
@@ -267,7 +267,7 @@
             if domain.geo_reference is not None:
                 domain.geo_reference.write_NetCDF(fid)
-	    #FIXME: Backwards compatibility
-	    fid.createVariable('z', self.precision, ('number_of_points',))
-	    #################################
+            #FIXME: Backwards compatibility
+            fid.createVariable('z', self.precision, ('number_of_points',))
+            #################################
 
             fid.createVariable('volumes', Int, ('number_of_volumes',
@@ -317,5 +317,5 @@
 
         # Get X, Y and bed elevation Z
-	Q = domain.quantities['elevation']
+        Q = domain.quantities['elevation']
         X,Y,Z,V = Q.get_vertex_values(xy=True,
 				      precision = self.precision)
@@ -327,8 +327,8 @@
         z[:] = Z.astype(self.precision)
 
-	#FIXME: Backwards compatibility
+        #FIXME: Backwards compatibility
         z = fid.variables['z']
         z[:] = Z.astype(self.precision)
-	################################
+        ################################
 
         volumes[:] = V.astype(volumes.typecode())
Index: /inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py	(revision 1362)
+++ /inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py	(revision 1363)
@@ -4,6 +4,6 @@
 stage_color = (0.1,0.4,0.99)
 
-class Surface:
-    def __init__(self,domain,scale_z=1.0,range_z=None):
+class Visualiser:
+    def __init__(self,domain,scale_z=1.0):
         """Create visualisation of domain
         """
@@ -14,4 +14,5 @@
         self.scale_z  = scale_z
         self.vertices = domain.vertex_coordinates
+
 
         # models for each quantity
@@ -39,8 +40,5 @@
 #        print 'range_y',self.range_y
 
-        self.range_z = 1.0
-        self.max_z = self.range_z/2.0
-        self.min_z = -self.range_z/2.0
-        self.range_z = max(self.max_z - self.min_z, 1.0e-10)
+
 
 
@@ -59,10 +57,24 @@
         self.border_model = curve(frame = self.frame,
                                   pos=[(0,0),(0,1),(1,1),(1,0),(0,0)])
-        self.timer=label(pos=(0.75,0.75,0.5),text='Time=%10.5e'%self.domain.time)
-
-        #scene.autoscale=0
+        self.timer=label(pos=(0.75,0.75,0.5),text='Time=%10.5e'%self.domain.time,
+                         visible=False)
+
+
+        scene.autoscale=1
         #self.update_all()
         #scene.uniform=0
-        self.setup_range_z()
+        self.range_z  = 1.0
+        if domain.visualise_range_z == None:
+            self.setup_range_z()
+        else:
+            self.range_z = domain.visualise_range_z
+
+        self.max_z = self.range_z/2.0
+        self.min_z = -self.range_z/2.0
+        self.range_z = max(self.max_z - self.min_z, 1.0e-10)
+
+        print self.range_z
+        print self.max_z
+        print self.min_z
 
     def setup_range_z(self,qname1='elevation',qname2='stage'):
@@ -70,5 +82,5 @@
         #print qname1
         #print qname2
-        self.range_z = 1.0e-10
+        range_z = 1.0e-10
         try:
             q1 = self.domain.quantities[qname1].vertex_values
@@ -76,7 +88,7 @@
             min_z = min(min(q1))
             print max_z, min_z
-            self.range_z = max(self.range_z, max_z - min_z)
+            range_z = max(range_z, max_z - min_z)
         except:
-            print 'could not find range of '+qname1
+            print 'Visualisation: could not find range of '+qname1
             pass
 
@@ -86,9 +98,11 @@
             min_z = min(min(q2))
             print max_z, min_z
-            self.range_z = max(self.range_z, max_z - min_z)
+            range_z = max(range_z, max_z - min_z)
         except:
             print 'Visualisation: could not find range of '+qname2
             pass
 
+        self.range_z = max(range_z,self.range_z)
+
     def update_all(self):
 
@@ -97,5 +111,8 @@
 
     def update_timer(self):
-        self.timer.text='Time=%10.5e'%self.domain.time
+
+        if self.domain.visualise_timer == True:
+            self.timer.visible = True
+            self.timer.text='Time=%10.5e'%self.domain.time
 
 
@@ -489,21 +506,21 @@
 
 
-def create_surface(domain,scale_z=0.5,range_z=None):
-
-    surface = Surface(domain,scale_z=0.5,range_z=range_z)
-    domain.surface = surface
-
-    surface.update_quantity('elevation')
-    surface.update_quantity('stage')
-    surface.update_timer()
+def create_surface(domain,scale_z=0.5):
+
+    visualiser = Visualiser(domain,scale_z=0.5)
+    domain.visualiser = visualiser
+
+    visualiser.update_quantity('elevation')
+    visualiser.update_quantity('stage')
+    visualiser.update_timer()
 
 def update(domain):
 
-    surface = domain.surface
+    visualiser = domain.visualiser
     if domain.visualise_color_stage:
-        surface.update_quantity_color('stage')
+        visualiser.update_quantity_color('stage')
     else:
-        surface.update_quantity('stage')
+        visualiser.update_quantity('stage')
 
     if domain.visualise_timer:
-        surface.update_timer()
+        visualiser.update_timer()
Index: /inundation/ga/storm_surge/pyvolution/shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 1362)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 1363)
@@ -83,5 +83,6 @@
         self.visualise_color_stage = False
         self.visualise_stage_range = 1.0
-        self.visualise_timer = False
+        self.visualise_timer = True
+        self.visualise_range_z = None
 
         #Stored output
Index: /inundation/ga/storm_surge/pyvolution/view_tsh.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/view_tsh.py	(revision 1362)
+++ /inundation/ga/storm_surge/pyvolution/view_tsh.py	(revision 1363)
@@ -43,5 +43,4 @@
     print "Number of triangles = ", len(domain)
 
-    surface =  Surface(domain, scale_z = scale_z)
-    surface.update_quantity('stage')
-
+    visualiser =  Visualiser(domain, scale_z = scale_z)
+    visualiser.update_quantity('stage')
Index: /inundation/ga/storm_surge/zeus/anuga-workspace.zwi
===================================================================
--- /inundation/ga/storm_surge/zeus/anuga-workspace.zwi	(revision 1362)
+++ /inundation/ga/storm_surge/zeus/anuga-workspace.zwi	(revision 1363)
@@ -2,5 +2,5 @@
 <workspace name="anuga-workspace">
     <mode>Debug</mode>
-    <active>Merimbula</active>
+    <active>parallel</active>
     <project name="Merimbula">Merimbula.zpi</project>
     <project name="parallel">parallel.zpi</project>
Index: /inundation/ga/storm_surge/zeus/parallel.zpi
===================================================================
--- /inundation/ga/storm_surge/zeus/parallel.zpi	(revision 1362)
+++ /inundation/ga/storm_surge/zeus/parallel.zpi	(revision 1363)
@@ -74,5 +74,7 @@
     <ReleaseProjectReload>Off</ReleaseProjectReload>
     <file>..\parallel\advection.py</file>
+    <file>..\parallel\parallel_advection.py</file>
     <file>..\parallel\run_advection.py</file>
+    <file>..\parallel\run_parallel_advection.py</file>
     <file>..\parallel\test_advection.py</file>
     <folder name="Header Files" />
