Index: /inundation/ga/storm_surge/pyvolution/data_manager.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 848)
+++ /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 849)
@@ -229,5 +229,9 @@
             fid.createVariable('x', self.precision, ('number_of_points',))
             fid.createVariable('y', self.precision, ('number_of_points',))
-            fid.createVariable('z', self.precision, ('number_of_points',))
+            fid.createVariable('elevation', self.precision, ('number_of_points',))
+            
+	    #FIXME: Backwards compatibility
+	    fid.createVariable('z', self.precision, ('number_of_points',))
+	    #################################
         
             fid.createVariable('volumes', Int, ('number_of_volumes',
@@ -272,5 +276,5 @@
         x = fid.variables['x']
         y = fid.variables['y']
-        z = fid.variables['z']
+        z = fid.variables['elevation']
         
         volumes = fid.variables['volumes']
@@ -286,4 +290,9 @@
         y[:] = Y.astype(self.precision)
         z[:] = Z.astype(self.precision)
+	
+	#FIXME: Backwards compatibility
+        z = fid.variables['z']
+        z[:] = Z.astype(self.precision)	
+	################################
 
         volumes[:] = V.astype(volumes.typecode())
@@ -403,5 +412,5 @@
             fid.createVariable('x', self.precision, ('number_of_points',))
             fid.createVariable('y', self.precision, ('number_of_points',))
-            #fid.createVariable('z', self.precision, ('number_of_points',))
+
         
             fid.createVariable('volumes', Int, ('number_of_volumes',
@@ -648,5 +657,5 @@
     x = fid.variables['x']
     y = fid.variables['y']
-    z = fid.variables['z']        
+    z = fid.variables['elevation']        
     time = fid.variables['time']
     stage = fid.variables['stage']
Index: /inundation/ga/storm_surge/pyvolution/general_mesh.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/general_mesh.py	(revision 848)
+++ /inundation/ga/storm_surge/pyvolution/general_mesh.py	(revision 849)
@@ -37,5 +37,11 @@
     
         #creates two triangles: bac and bce
-        
+
+    Other:
+    
+      In addition mesh computes an Nx6 array called vertex_coordinates.
+      This structure is derived from coordinates and contains for each 
+      triangle the three x,y coordinates at the vertices.
+         		        
 
         This is a cut down version of mesh from pyvolution\mesh.py
Index: /inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py	(revision 848)
+++ /inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py	(revision 849)
@@ -175,9 +175,9 @@
         from Numeric import array
         from config import time_format
-        from util import File_function
-        
-        Boundary.__init__(self)
-
-        self.F = File_function(filename, domain)
+        from util import file_function
+        
+        Boundary.__init__(self)
+
+        self.F = file_function(filename, domain)
         self.domain = domain
 
@@ -205,5 +205,5 @@
 
         d = len(domain.conserved_quantities)
-        msg = 'Values specified in file must be a list or an array of length %d' %d
+        msg = 'Values specified in file %s must be a list or an array of length %d' %(filename, d)
         assert len(q) == d, msg
 
Index: /inundation/ga/storm_surge/pyvolution/interpolate_sww.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/interpolate_sww.py	(revision 848)
+++ /inundation/ga/storm_surge/pyvolution/interpolate_sww.py	(revision 849)
@@ -106,5 +106,5 @@
         try:
             if quantity_name == 'height':
-                z = fid.variables['z'][:]
+                z = fid.variables['elevation'][:]
                 stage = fid.variables['stage'][:,:]
                 quantity = stage - z  # 2D, using broadcasting 
Index: /inundation/ga/storm_surge/pyvolution/mesh.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/mesh.py	(revision 848)
+++ /inundation/ga/storm_surge/pyvolution/mesh.py	(revision 849)
@@ -409,6 +409,4 @@
                 break
 
-            #Should this warning be an exception?
-            #assert v is not None, msg
 
                 
@@ -453,6 +451,8 @@
         #for id, edge in self.boundary:
         #    assert self.neighbours[id,edge] < 0
-
-
+        #
+	#NOTE (Ole): I reckon this was resolved late 2004?
+	#
+	#See domain.set_boundary
         
             
Index: /inundation/ga/storm_surge/pyvolution/shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 848)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 849)
@@ -22,14 +22,16 @@
 
 symbol    variable name    explanation 
-z         elevation        elevation of bed on which flow is modelled
-h         height           water height above z
-w         stage            water level, w = z+h
-u                          speed in the x direction
-v                          speed in the y direction
-uh        xmomentum        momentum in the x direction 
-vh        ymomentum        momentum in the y direction
-
-eta                        mannings friction coefficient
-nu                         wind stress coefficient
+x         x                horizontal distance from origin [m]
+y         y                vertical distance from origin [m] 
+z         elevation        elevation of bed on which flow is modelled [m]
+h         height           water height above z [m]
+w         stage            water level, w = z+h [m]
+u                          speed in the x direction [m/s]
+v                          speed in the y direction [m/s]
+uh        xmomentum        momentum in the x direction [m^2/s]
+vh        ymomentum        momentum in the y direction [m^2/s]
+
+eta                        mannings friction coefficient []
+nu                         wind stress coefficient []
 
 The conserved quantities are w, uh, vh
@@ -111,11 +113,4 @@
         assert self.conserved_quantities[2] == 'ymomentum', msg
         
-
-        #Check that stages are >= bed elevation
-        from Numeric import alltrue, greater_equal
-        
-        stage = self.quantities['stage']
-        bed = self.quantities['elevation']        
-
 
     def compute_fluxes(self):
@@ -766,13 +761,17 @@
 
 	
-class Spatio_temporal_boundary(Boundary):
-    """The spatio-temporal boundary, reads values for the conserved 
-    quantities from an sww NetCDF file, and returns interpolated values 
-    at the midpoints of each associated boundaty segment.
-    Time dependency is interpolated linearly as in util.File_function.
-      
-    """
-
-    pass
+#class Spatio_temporal_boundary(Boundary):
+#    """The spatio-temporal boundary, reads values for the conserved 
+#    quantities from an sww NetCDF file, and returns interpolated values 
+#    at the midpoints of each associated boundaty segment.
+#    Time dependency is interpolated linearly as in util.File_function.#
+#
+#    Example:
+#    Bf = Spatio_temporal_boundary('source_file.sww', domain) 
+#          
+#    """
+Spatio_temporal_boundary = File_boundary
+        
+
 	
 
Index: /inundation/ga/storm_surge/pyvolution/test_data_manager.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_data_manager.py	(revision 848)
+++ /inundation/ga/storm_surge/pyvolution/test_data_manager.py	(revision 849)
@@ -149,5 +149,5 @@
         x = fid.variables['x']
         y = fid.variables['y']
-        z = fid.variables['z']
+        z = fid.variables['elevation']
 
         volumes = fid.variables['volumes']        
@@ -196,5 +196,5 @@
         x = fid.variables['x']
         y = fid.variables['y']
-        z = fid.variables['z']
+        z = fid.variables['elevation']
 
         volumes = fid.variables['volumes']        
@@ -256,5 +256,5 @@
         x = fid.variables['x']
         y = fid.variables['y']
-        z = fid.variables['z']        
+        z = fid.variables['elevation']        
         time = fid.variables['time']
         stage = fid.variables['stage']
@@ -314,5 +314,5 @@
         x = fid.variables['x']
         y = fid.variables['y']
-        z = fid.variables['z']
+        z = fid.variables['elevation']
         time = fid.variables['time']
         stage = fid.variables['stage']        
@@ -374,5 +374,5 @@
         x = fid.variables['x']
         y = fid.variables['y']
-        z = fid.variables['z']
+        z = fid.variables['elevation']
         time = fid.variables['time']
         stage = fid.variables['stage']        
@@ -459,5 +459,5 @@
         x = fid.variables['x']
         y = fid.variables['y']
-        z = fid.variables['z']
+        z = fid.variables['elevation']
 
         volumes = fid.variables['volumes']   
Index: /inundation/ga/storm_surge/pyvolution/test_interpolate_sww.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_interpolate_sww.py	(revision 848)
+++ /inundation/ga/storm_surge/pyvolution/test_interpolate_sww.py	(revision 849)
@@ -98,5 +98,5 @@
         x = fid.variables['x']
         y = fid.variables['y']
-        z = fid.variables['z']
+        z = fid.variables['elevation']
 
         volumes = fid.variables['volumes']   
@@ -218,5 +218,5 @@
         # look at error catching
         try:
-            interp = Interpolate_sww(sww.filename, 'z')
+            interp = Interpolate_sww(sww.filename, 'elevation')
         except KeyError:
             pass
Index: /inundation/ga/storm_surge/pyvolution/test_shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 848)
+++ /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 849)
@@ -5,5 +5,5 @@
 
 from config import g, epsilon
-from Numeric import allclose, array, zeros, ones, Float
+from Numeric import allclose, array, zeros, ones, Float, take
 from shallow_water import *
 
@@ -949,5 +949,5 @@
         from math import pi, cos, sin, sqrt
         from config import time_format
-        from util import File_function
+        from util import file_function
         import time
         
@@ -995,5 +995,5 @@
 
         #Setup wind stress
-        F = File_function(filename)
+        F = file_function(filename)
         W = Wind_stress(F) 
         domain.forcing_terms = []
@@ -2380,5 +2380,5 @@
        
 	#Create sww file of simple propagation from left to right 
-	#through rectangular domain: domain1
+	#through rectangular domain
 	
         from mesh_factory import rectangular
@@ -2388,52 +2388,111 @@
 
         #Create shallow water domain
-        domain = Domain(points, vertices, boundary)
-        domain.smooth = False
-
-	domain.visualise = False
-        domain.default_order = 2
-
-        domain.set_datadir('.')
-	domain.store = True	
-        domain.set_name('test_spatio_temporal_boundary_' + str(time.time()))
-	
+        domain1 = Domain(points, vertices, boundary)
+	
+        from util import mean
+        domain1.reduction = mean 
+        domain1.smooth = True #FIXME: This may be a requirement for the 
+	                     #interpolation process 
+
+	#domain1.visualise = True
+        domain1.default_order = 2
+
+	domain1.store = True	
+        domain1.set_datadir('.')
+        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
+
+		
         #Bed-slope and friction at vertices (and interpolated elsewhere)
-        domain.set_quantity('elevation', 0)
-        domain.set_quantity('friction', 0)	
+        domain1.set_quantity('elevation', 0)
+        domain1.set_quantity('friction', 0)	
 
         # Boundary conditions
-        Br = Reflective_boundary(domain)
+        Br = Reflective_boundary(domain1)
         Bd = Dirichlet_boundary([0.3,0,0])	
-	Bt = Transmissive_boundary(domain)
-        domain.set_boundary({'left': Bd, 'right': Bt, 'top': Br, 'bottom': Br})
+	Bt = Transmissive_boundary(domain1)
+        domain1.set_boundary({'left': Bd, 'right': Bt, 'top': Br, 'bottom': Br})
 
         #Initial condition
-        domain.set_quantity('stage', 0)
-        domain.check_integrity()
+        domain1.set_quantity('stage', 0)
+        domain1.check_integrity()
         
         #Evolution
-        for t in domain.evolve(yieldstep = 0.05, finaltime = 4):
-	    import time
-	    #time.sleep(0.01)
-            #domain.write_time()            
-
-
-	
-	
-	
-	#Create an triangle shaped domain, formed from the 
-	#lower and right hand  boundaries and the sw-ne diagonal 
-	#from domain 1. Called it domain2
+        for t in domain1.evolve(yieldstep = 0.05, finaltime = 4):
+            pass
+
+        cv1 = domain1.quantities['stage'].centroid_values
+	
+        #print take(cv1, (12,14,16))  #Right
+        #print take(cv1, (0,6,12))  #Bottom
+	
+	
+	#Create an triangle shaped domain (reusing coordinates from domain 1), 
+	#formed from the lower and right hand  boundaries and 
+	#the sw-ne diagonal 
+	#from domain 1. Call it domain2
+	
+	points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
+	           [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
+		   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
+	
+	vertices = [ [1,2,0], 
+	             [3,4,1], [2,1,4], [4,5,2], 
+	             [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]	
+
+        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
+	             (4,2):'right', (6,2):'right', (8,2):'right',
+		     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}    
+		     
+        domain2 = Domain(points, vertices, boundary)
+	
+	domain2.reduction = domain1.reduction
+        domain2.smooth = False
+        domain2.visualise = True
+        domain2.default_order = 2
+	
+        #Bed-slope and friction at vertices (and interpolated elsewhere)
+        domain2.set_quantity('elevation', 0)
+        domain2.set_quantity('friction', 0)	
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain2)
+        #Bd = Dirichlet_boundary([0.3,0,0])	
+	Bt = Transmissive_boundary(domain2)
+	
+        return
+	
+	Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
+				      domain2)
+        domain2.set_boundary({'right':Bt, 'bottom':Br, 'diagonal':Bf})
+
+        #Initial condition 
+        domain2.set_quantity('stage', 0)
+        domain2.check_integrity()        
+
+        #Evolution
+        for t in domain2.evolve(yieldstep = 0.05, finaltime = 4):
+            pass
+
 	
 	#Use output from domain1 as spatio-temporal boundary for domain2
         #and verify that results at right hand side are close.	
 
-	
-	
+        cv2 = domain2.quantities['stage'].centroid_values
+	
+        print take(cv1, (12,14,16))  #Right
+	print take(cv2, (4,6,8))
+	
+        #print take(cv1, (0,6,12))  #Bottom	
+	
+	assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
+	assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom	
+	
+        #Cleanup
+        os.remove(domain1.filename + '.' + domain1.format)
     
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
-    #suite = unittest.makeSuite(TestCase,'test_boundary_conditionsII')
+    #suite = unittest.makeSuite(TestCase,'test')
+    suite = unittest.makeSuite(TestCase,'test_spatio_temporal_boundary')
     runner = unittest.TextTestRunner()
     runner.run(suite)
Index: /inundation/ga/storm_surge/pyvolution/test_util.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 848)
+++ /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 849)
@@ -166,5 +166,5 @@
         fid.close()
 
-        F = File_function(filename)
+        F = file_function(filename)
 
         #Now try interpolation
@@ -229,15 +229,15 @@
 
         #Check that domain.starttime is updated if non-existing
-        F = File_function(filename, domain)
+        F = file_function(filename, domain)
         assert allclose(domain.starttime, start)        
 
         #Check that domain.starttime is updated if too early
         domain.starttime = start - 1                
-        F = File_function(filename, domain)
+        F = file_function(filename, domain)
         assert allclose(domain.starttime, start)
 
         #Check that domain.starttime isn't updated if later
         #domain.starttime = start + 1                
-        #F = File_function(filename, domain)
+        #F = file_function(filename, domain)
         #assert allclose(domain.starttime, start+1)                
 
@@ -307,5 +307,5 @@
 	delta = 23
         domain.starttime = start + delta                
-        F = File_function(filename, domain)
+        F = file_function(filename, domain)
         assert allclose(domain.starttime, start+delta)                
 
Index: /inundation/ga/storm_surge/pyvolution/util.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/util.py	(revision 848)
+++ /inundation/ga/storm_surge/pyvolution/util.py	(revision 849)
@@ -80,6 +80,316 @@
       return False  
     
-	     
-class File_function:
+
+def file_function(filename, domain=None):  
+    assert type(filename) == type(''),\
+               'First argument to File_function must be a string'
+
+    try:
+        fid = open(filename)
+    except Exception, e:
+        msg = 'File "%s" could not be opened: Error="%s"'\
+                  %(filename, e)
+        raise msg
+
+    line = fid.readline()
+    fid.close()
+
+    #Choose format 
+    #FIXME: Maybe these can be merged later on
+    if line[:3] == 'CDF':
+        return File_function_NetCDF(filename, domain)
+    else:
+        return File_function_ASCII(filename, domain)    
+	   
+    	     
+class File_function_NetCDF:
+    """Read time history of spatial data from NetCDF sww file and 
+    return a callable object f(t,x,y)  
+    which will return interpolated values based on the input file.    
+    
+    x, y may be either scalars or vectors
+    
+    #FIXME: More about format, interpolation  and order of quantities
+    """
+
+    def __init__(self, filename, domain=None):
+        """Initialise callable object from a file.
+        See docstring for class File_function_NetCDF
+
+        If domain is specified, model time (domain.starttime)
+        will be checked and possibly modified
+
+        All times are assumed to be in UTC
+        """
+
+        import time, calendar
+        from Numeric import array
+        from config import time_format
+        from Scientific.IO.NetCDF import NetCDFFile	
+
+        #Open NetCDF file
+        fid = NetCDFFile(filename, 'r')
+    
+        print fid.variables.keys()
+    
+        fields = line.split(',')
+        msg = 'File %s must have the format date, value0 value1 value2 ...'
+        msg += ' or date, x, y, value0 value1 value2 ...'
+        mode = len(fields)
+        assert mode in [2,4], msg
+
+        try:
+            starttime = calendar.timegm(time.strptime(fields[0], time_format))
+        except ValueError:    
+            msg = 'First field in file %s must be' %filename
+            msg += ' date-time with format %s.\n' %time_format
+            msg += 'I got %s instead.' %fields[0] 
+            raise msg
+
+
+        #Split values
+        values = []
+        for value in fields[mode-1].split():
+            values.append(float(value))
+
+        q = array(values)
+        
+        msg = 'ERROR: File must contain at least one independent value'
+        assert len(q.shape) == 1, msg
+            
+        self.number_of_values = len(q)
+
+        self.filename = filename
+        self.starttime = starttime
+	self.domain = domain
+
+        if domain is not None:
+            if domain.starttime is None:
+                domain.starttime = self.starttime
+            else:
+                msg = 'WARNING: Start time as specified in domain (%s)'\
+                      %domain.starttime
+                msg += ' is earlier than the starttime of file %s: %s.'\
+                         %(self.filename, self.starttime)
+                msg += 'Modifying starttime accordingly.'
+                if self.starttime > domain.starttime:
+                    #FIXME: Print depending on some verbosity setting
+                    #print msg
+                    domain.starttime = self.starttime #Modifying model time
+
+        if mode == 2:        
+            self.read_times() #Now read all times in.
+        else:
+            raise 'x,y dependency not yet implemented'
+
+    
+    	     
+class File_function_ASCII:
+    """Read time series from file and return a callable object:
+    f(t,x,y)
+    which will return interpolated values based on the input file.
+
+    The file format is assumed to be either two fields separated by a comma:
+
+        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
+
+    or four comma separated fields
+
+        time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...
+
+    In either case, the callable object will return a tuple of
+    interpolated values, one each value stated in the file.
+
+    
+    E.g
+
+      31/08/04 00:00:00, 1.328223 0 0
+      31/08/04 00:15:00, 1.292912 0 0
+
+    will provide a time dependent function f(t,x=None,y=None),
+    ignoring x and y, which returns three values per call.
+
+    
+    NOTE: At this stage the function is assumed to depend on
+    time only, i.e  no spatial dependency!!!!!
+    When that is needed we can use the least_squares interpolation.
+    
+    #FIXME: This should work with netcdf (e.g. sww) and thus render the 
+    #spatio-temporal boundary condition in shallow water fully general
+    """
+
+    
+    def __init__(self, filename, domain=None):
+        """Initialise callable object from a file.
+        See docstring for class File_function
+
+        If domain is specified, model time (domain,starttime)
+        will be checked and possibly modified
+
+        All times are assumed to be in UTC
+        """
+
+        import time, calendar
+        from Numeric import array
+        from config import time_format
+
+        fid = open(filename)
+        line = fid.readline()
+        fid.close()
+    
+        fields = line.split(',')
+        msg = 'File %s must have the format date, value0 value1 value2 ...'
+        msg += ' or date, x, y, value0 value1 value2 ...'
+        mode = len(fields)
+        assert mode in [2,4], msg
+
+        try:
+            starttime = calendar.timegm(time.strptime(fields[0], time_format))
+        except ValueError:    
+            msg = 'First field in file %s must be' %filename
+            msg += ' date-time with format %s.\n' %time_format
+            msg += 'I got %s instead.' %fields[0] 
+            raise msg
+
+
+        #Split values
+        values = []
+        for value in fields[mode-1].split():
+            values.append(float(value))
+
+        q = array(values)
+        
+        msg = 'ERROR: File must contain at least one independent value'
+        assert len(q.shape) == 1, msg
+            
+        self.number_of_values = len(q)
+
+        self.filename = filename
+        self.starttime = starttime
+	self.domain = domain
+
+        if domain is not None:
+            if domain.starttime is None:
+                domain.starttime = self.starttime
+            else:
+                msg = 'WARNING: Start time as specified in domain (%s)'\
+                      %domain.starttime
+                msg += ' is earlier than the starttime of file %s: %s.'\
+                         %(self.filename, self.starttime)
+                msg += 'Modifying starttime accordingly.'
+                if self.starttime > domain.starttime:
+                    #FIXME: Print depending on some verbosity setting
+                    #print msg
+                    domain.starttime = self.starttime #Modifying model time
+
+        if mode == 2:        
+            self.read_times() #Now read all times in.
+        else:
+            raise 'x,y dependency not yet implemented'
+
+        
+    def read_times(self):
+        """Read time and values 
+        """
+        from Numeric import zeros, Float, alltrue
+        from config import time_format
+        import time, calendar
+        
+        fid = open(self.filename)        
+        lines = fid.readlines()
+        fid.close()
+        
+        N = len(lines)
+        d = self.number_of_values
+
+        T = zeros(N, Float)       #Time
+        Q = zeros((N, d), Float)  #Values
+        
+        for i, line in enumerate(lines):
+            fields = line.split(',')
+            realtime = calendar.timegm(time.strptime(fields[0], time_format))
+
+            T[i] = realtime - self.starttime
+
+            for j, value in enumerate(fields[1].split()):
+                Q[i, j] = float(value)
+
+        msg = 'File %s must list time as a monotonuosly ' %self.filename
+        msg += 'increasing sequence'
+        assert alltrue( T[1:] - T[:-1] > 0 ), msg
+
+        self.T = T     #Time
+        self.Q = Q     #Values
+        self.index = 0 #Initial index
+
+        
+    def __repr__(self):
+        return 'File function'
+
+        
+
+    def __call__(self, t, x=None, y=None):
+        """Evaluate f(t,x,y)
+
+        FIXME: x, y dependency not yet implemented except that
+        result is a vector of same length as x and y
+        FIXME: Naaaa
+        """
+
+        from math import pi, cos, sin, sqrt
+        
+
+        #Find time tau such that
+	#
+	# domain.starttime + t == self.starttime + tau
+	
+	if self.domain is not None:
+            tau = self.domain.starttime-self.starttime+t
+	else:
+            tau = t
+	
+	        
+        msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
+              %(self.filename, self.T[0], self.T[1], tau)
+        if tau < self.T[0]: raise msg
+        if tau > self.T[-1]: raise msg        
+
+        oldindex = self.index
+
+        #Find slot
+        while tau > self.T[self.index]: self.index += 1
+        while tau < self.T[self.index]: self.index -= 1
+
+        #t is now between index and index+1
+        ratio = (tau - self.T[self.index])/\
+                (self.T[self.index+1] - self.T[self.index])
+
+        #Compute interpolated values
+        q = self.Q[self.index,:] +\
+            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 
+
+        #Return vector of interpolated values
+        if x == None and y == None:
+            return q
+        else:
+            try:
+                N = len(x)
+            except:
+                return q
+            else:
+                from Numeric import ones, Float
+                #x is a vector - Create one constant column for each value
+                N = len(x)
+                assert len(y) == N, 'x and y must have same length'
+             
+                res = []
+                for col in q:
+                    res.append(col*ones(N, Float))
+            
+                return res
+            
+	    
+#FIXME: TEMP	    
+class File_function_copy:
     """Read time series from file and return a callable object:
     f(t,x,y)
