Merge lp:~fluidity-core/fluidity/freesurfacefromnetcdffix into lp:fluidity

[Adam Candy]

Fix to bug lp:922025. This tidies up the manner in which the pressure initial condition is set from a free-surface height field provided in a netCDF file.

Summary of changes:
 - The special case, legacy method has been removed.
 - If the input format is not recognised, the code will abort with an error detailing the problematic field and the schema path of the invalid format.
 - The 'standard' format is now named 'raw' (instead of 'other'). This is the case where no post-processing is performed on the source data, and the field is initialised by the raw data present in the netCDF file.
 - Included logging messages to indicate which field is being populated and the source netCDF file (logging level 1).
 - Logging of the post-processing applied (logging level 3).
 - For the pressure boundary condition set from a free-surface field in particular, logging level 3 will report the min and max values of the field before and after the conversion.
 - Documentation in the schema detailing the available format options and the renaming described above.

There are currently no tests, examples or longtests that are affected by the above change, but possibly user-generated initialisation flml files that do use the legacy implementation of the pressure initial condition set from free-surface netCDF. In this case, the additional logging added with this change should make it easy to identify the change required to adopt the improved version of this import.

In short, the following example highlights the change necessary:

The option path:
  /material_phase::Ocean/scalar_field::Pressure/prognostic/initial_condition::WholeMesh/from_netcdf/format/Free-surface height

should become:
  /material_phase::Ocean/scalar_field::Pressure/prognostic/initial_condition::WholeMesh/free_surface/from_netcdf/format/raw

(where the file location is specified under the 'from_netcdf' node in each case).

[lp:~fluidity-core/fluidity/freesurfacefromnetcdffix]

3912. By Adam Candy on 2012-02-09

A new section in the manual describing how to initialise the pressure field by specifying the free surface height.
Removal of the description of the legacy method to achieve this.
Other corrections to nearby sections.

[Adam Candy (asc)]

The manual has been updated to reflect the changes made.

[Jon Hill (jon-hill)]

Looks good to me.

One minor typo in schemas/input_output.rnc
Data should be in the form x,y,z(y,x) where where x and y are the coordinates ->
Data should be in the form x,y,z(y,x) where x and y are the coordinates

Can we create a quick test of this functionality, including doctoring netcdf files to be invalid? How much effort would that be, do you think?

[lp:~fluidity-core/fluidity/freesurfacefromnetcdffix]

3913. By Adam Candy on 2012-02-10

Merge in of the trunk.

[Adam Candy (asc)]

Thanks for taking a look Jon - the typo has been fixed. Also noticed a spelling mistake above too. These have all been in there for quite some time now!

I'm working on a test for the functionality - see:
https://code.launchpad.net/~fluidity-core/fluidity/netcdftests

The main test is in place, just looking at how to check invalid files.

I've also discovered a bug in the reader that will go up shortly...

[lp:~fluidity-core/fluidity/freesurfacefromnetcdffix]

3914. By Adam Candy on 2012-02-11

Correct the typo highlighted by Jon (thanks) and a spelling.

3915. By Adam Candy on 2012-02-11

Merge from trunk (to pick up improved netcdf error reporting changes).

3916. By Adam Candy on 2012-02-11

New test case: netcdf_read_free_surface
This creates a valid netCDF file containing a free surface height in the domain x,y in [-1.2,1.2]. This is read in by Fluidity as a free surface and applied as an initial condition on pressure. From this, Fluidity calculates the free surface as a diagnostic, which is provided in the outputted VTU file, alongside a field containing the analytic function used to create the netCDF file in the first place (here described in height.py). The test compares the difference between these two fields, which should be small, due to the interpolation of the netCDF file.

In summary this case tests:
  - NetCDF reader
  - Reading of a free surface and conversion into a pressure initial condition
  - Creation of the free surface diagnostic from the pressure field (including its extrapolation down from the top surface)
  - Creation of prescribed scalar fields from Python

3917. By Adam Candy on 2012-02-11

A new test case: netcdf_read_errors
This checks the errors reported when an attempt to read an invalid NetCDF file is made. Here 'invalid' refers to a file that is itself an invalid, non-conforming NetCDF file, or is of a format not supported by Fluidity.

Five types of NetCDF files are checked:
  - One which does not provide all the data requested,
  - One is missing a dimension (i.e. it provides only a 1d field),
  - One is missing the main variable,
  - One contains an incorrectly-labelled dimension (as 'lat', not 'x'), and
  - One has the variable incorrectly labelled (as 'height', not 'z').

[Adam Candy (asc)]

Please take a look at the two new test cases added, 'netcdf_read_free_surface' and 'netcdf_read_errors'. Combined, these two test cases check the free surface initialisation discussed here and in bug lp:922025, and also the error reporting by the netCDF reader, tackled by bug lp:920944. They test:
  - NetCDF reader
  - Reading of a free surface and conversion into a pressure initial condition
  - Creation of the free surface diagnostic from the pressure field (including its extrapolation down from the top surface)
  - Creation of prescribed scalar fields from Python
and five types of invalid NetCDF files:
Five types of NetCDF files are checked:
  - One which does not provide all the data requested,
  - One is missing a dimension (i.e. it provides only a 1d field),
  - One is missing the main variable,
  - One contains an incorrectly-labelled dimension (as 'lat', not 'x'), and
  - One has the variable incorrectly labelled (as 'height', not 'z').

The NetCDF reader is now tested!

[Adam Candy (asc)]

The above addresses Jon's suggestions, and I think it is fine to be merged and then released in 4.1.4.

Unmerged revisions

=== modified file 'manual/configuring_fluidity.tex'
--- manual/configuring_fluidity.tex	2011-11-23 11:03:53 +0000
+++ manual/configuring_fluidity.tex	2012-02-11 17:18:19 +0000
@@ -873,19 +873,31 @@
 
 The initial state of certain fields can be set from external data contained within a NetCDF file.
 This functionality can be used by selecting the \option{\ldots/initial\_condition/from\_netcdf/format}
-option. This feature is currently intended for setting, for example, the initial free surface elevation
-in a Tsunami simulation. This option will not currently work with multi-layered data files. Supported
+option.
+This option will not currently work with multi-layered data files. Supported
 NetCDF file conventions include the NetCDF-CF 1.x convention.
 
 \subsubsection{Setting fields from NEMO data}\label{sec:setting_from_nemo}
 \index{field!from Nemo}
 Initial conditions of prognostic fields and the values of prescribed fields can also be set from an external NEMO
 data file. The external data file is in the NETCDF format and data is currently available for pressure, temperature, salinity
-and velocity. To set the initial condition of a prognostic field from NEMO data, chose the option 
+and velocity. To set the initial condition of a prognostic field from NEMO data, choose the option 
 \option{\ldots/prognostic/initial\_condition/NEMO\_data} and then under \option{format} select the required data format. For scalar fields
 the formats available are `Temperature', `Salinity' and `Free-surface height', for vector fields `Velocity' is the only available format.
 Setting the value of prescribed fields from NEMO data works similarly. Set the option \option{\ldots/prescribed/value/NEMO\_data} and then proceed as above.
 
+\subsubsection{Setting an initial free surface height}\label{sec:setting_free_surface_height}
+\index{inital conditions!free surface height}
+\index{free surface!initial condition}
+The free surface height is contained within the Pressure field.  To apply an initial condition on free surface height, 
+choose the \option{\ldots/free\_surface} under the relevant Pressure initial condition option.
+With this option, it is possible to set the initial free surface elevation in a tsunami simulation, for example.
+The initial condition can be applied using the approaches outlined above, in this section~\ref{sec:setting_field_values}.
+It is recommended that a diagnostic Free Surface field is included if this option is used.
+
+If set from a NetCDF file using the option \option{\ldots/initial\_condition/from\_netcdf}, and the file provides exactly the
+free surface height, it is important that the child option \option{\ldots/initial\_condition/from\_netcdf/format} is set to `raw'.
+
 \subsection{Region IDs}
 \index{region ID}
 If the input mesh defines a number of region IDs then these may be employed

=== modified file 'ocean_forcing/load_netcdf.F90'
--- ocean_forcing/load_netcdf.F90	2010-12-17 16:50:11 +0000
+++ ocean_forcing/load_netcdf.F90	2012-02-11 17:18:19 +0000
@@ -31,14 +31,20 @@
   call get_option(trim(path)//"/from_netcdf/format", format)
 
   select case (format)
-    ! This is currently the only 'special' case, but may change as people need to initialize
-    ! different fields
-    case ("Free-surface height")
-      call get_option("/physical_parameters/gravity/magnitude", gravity_magnitude, stat=stat)
-      if (stat/=0) then
-        FLAbort("Trying to initialize an initial condition for the free surface without a gravitational field")
-      endif
-      call scale(field, gravity_magnitude)
+    ! It is possible to manipulate imported netCDF data here.
+    ! Note that the 'Free-surface height' format option is no longer supported.
+    ! In this case the 'free_surface' node under pressure initial conditions should be used.
+    
+    case ("raw")
+      ewrite(3,*) "The data used to initialise field " // trim(field%name) // &
+        ", has been treated as the raw values and not post-processed after import."
+
+    case default
+      ewrite(-1,*) "The format " // trim(format) // ", is not a recognised method to handle netCDF files."
+      ewrite(-1,*) "Please specify a valid format for initialiasing the field " // trim(field%name) // &
+        ", at " // trim(path) // "/from_netcdf/format"
+      FLAbort("Fatal error initialising a field from a netCDF file.")
+    
   end select
 
 end subroutine
@@ -51,14 +57,16 @@
   real, dimension(position%dim,node_count(position)) :: temp_pos
   character(len=FIELD_NAME_LEN) :: filename
   real, dimension(:), allocatable :: X, Y, Z
-  
   integer :: NNodes, i
 
+  assert(node_count(field)==node_count(position))
+
   call get_option(trim(path)//"/from_netcdf/file_name", filename)
 
-  assert(node_count(field)==node_count(position))
+  ewrite(1,*) "Populating field " // trim(field%name) // &
+    ", with data from the netCDF file: " // trim(filename)
 
-  NNodes=node_count(field)
+  NNodes = node_count(field)
 
   allocate(X(NNodes), Y(NNodes), Z(NNodes))
 
@@ -98,4 +106,4 @@
 
 end subroutine
 
-end module load_netcdf_module
\ No newline at end of file
+end module load_netcdf_module

=== modified file 'preprocessor/Initialise_Fields.F90'
--- preprocessor/Initialise_Fields.F90	2011-11-14 18:13:34 +0000
+++ preprocessor/Initialise_Fields.F90	2012-02-11 17:18:19 +0000
@@ -111,7 +111,9 @@
     else if(have_option(trim(path)//"/free_surface")) then
        call initialise_field(field, trim(path)//"/free_surface", position, time=time, phase_path=trim(phase_path))
        ! Scale the entered freesurface height by the magnitude of gravity to give the pressure initial condition
+       ewrite(3, *) "Free-surface to pressure minmax pre-conversion ", minval(field), maxval(field)
        call convert_free_surface_to_pressure(field, phase_path)
+       ewrite(3, *) "Free-surface to pressure minmax post-conversion ", minval(field), maxval(field)
 
     else if(have_option(trim(path) // "/from_file")) then
        

=== modified file 'schemas/input_output.rnc'
--- schemas/input_output.rnc	2011-07-27 16:57:15 +0000
+++ schemas/input_output.rnc	2012-02-11 17:18:19 +0000
@@ -180,15 +180,18 @@
           }
       }|
       ## Initialize a single field from data contained within a NetCDF file.
-      ## Data should be in the form x,y,z(y,x) where where x and y are the coordinates
-      ## (either cartesian or longitude-latitude) and z is the value of the intial
+      ## Data should be in the form x, y, z(y,x) where x and y are the coordinates
+      ## (either cartesian or longitude-latitude) and z is the value of the initial
       ## condition to be imposed (or in the GEBCO, COARDS or CF-1.0 formats). 
       ## Note: This option currently only supports 'single layer' data sets.
+      ##
+      ## Current format options are: raw
+      ## Simulations using the format 'other' should be updated to use 'raw'.
       element from_netcdf {
           attribute file_name { xsd:string },
           element format {
               element string_value {
-                  "Free-surface height"|"other"
+                  "raw"
               }
           }
       }|
@@ -249,15 +252,18 @@
           }
       }|
       ## Initialize a single field from data contained within a NetCDF file.
-      ## Data should be in the form x,y,z where where x and y are the coordinates
-      ## (either cartesian or longitude-latitude) and z is the value of the intial
-      ## condition to be imposed. 
+      ## Data should be in the form x, y, z(y,x) where x and y are the coordinates
+      ## (either cartesian or longitude-latitude) and z is the value of the initial
+      ## condition to be imposed (or in the GEBCO, COARDS or CF-1.0 formats). 
       ## Note: This option currently only supports 'single layer' data sets.
+      ##
+      ## Current format options are: raw
+      ## Simulations using the format 'other' should be updated to use 'raw'.
       element from_netcdf {
           attribute file_name { xsd:string },
           element format {
               element string_value {
-                  "Free-surface height"|"other"
+                  "raw"
               }
           }
       }|

=== modified file 'schemas/input_output.rng'
--- schemas/input_output.rng	2011-07-27 16:57:15 +0000
+++ schemas/input_output.rng	2012-02-11 17:18:19 +0000
@@ -229,19 +229,19 @@
       </element>
       <element name="from_netcdf">
         <a:documentation>Initialize a single field from data contained within a NetCDF file.
-Data should be in the form x,y,z(y,x) where where x and y are the coordinates
-(either cartesian or longitude-latitude) and z is the value of the intial
+Data should be in the form x, y, z(y,x) where x and y are the coordinates
+(either cartesian or longitude-latitude) and z is the value of the initial
 condition to be imposed (or in the GEBCO, COARDS or CF-1.0 formats). 
-Note: This option currently only supports 'single layer' data sets.</a:documentation>
+Note: This option currently only supports 'single layer' data sets.
+
+Current format options are: raw
+Simulations using the format 'other' should be updated to use 'raw'.</a:documentation>
         <attribute name="file_name">
           <data type="string"/>
         </attribute>
         <element name="format">
           <element name="string_value">
-            <choice>
-              <value>Free-surface height</value>
-              <value>other</value>
-            </choice>
+            <value>raw</value>
           </element>
         </element>
       </element>
@@ -309,19 +309,19 @@
       </element>
       <element name="from_netcdf">
         <a:documentation>Initialize a single field from data contained within a NetCDF file.
-Data should be in the form x,y,z where where x and y are the coordinates
-(either cartesian or longitude-latitude) and z is the value of the intial
-condition to be imposed. 
-Note: This option currently only supports 'single layer' data sets.</a:documentation>
+Data should be in the form x, y, z(y,x) where x and y are the coordinates
+(either cartesian or longitude-latitude) and z is the value of the initial
+condition to be imposed (or in the GEBCO, COARDS or CF-1.0 formats). 
+Note: This option currently only supports 'single layer' data sets.
+
+Current format options are: raw
+Simulations using the format 'other' should be updated to use 'raw'.</a:documentation>
         <attribute name="file_name">
           <data type="string"/>
         </attribute>
         <element name="format">
           <element name="string_value">
-            <choice>
-              <value>Free-surface height</value>
-              <value>other</value>
-            </choice>
+            <value>raw</value>
           </element>
         </element>
       </element>

=== added directory 'tests/netcdf_read_errors'
=== added file 'tests/netcdf_read_errors/Makefile'
--- tests/netcdf_read_errors/Makefile	1970-01-01 00:00:00 +0000
+++ tests/netcdf_read_errors/Makefile	2012-02-11 17:18:19 +0000
@@ -0,0 +1,18 @@
+
+default: run
+
+run: input
+	../../bin/fluidity -v 3 -l missingdata.flml 
+
+%.flml: valid.flml
+	@sed -e 's/valid\.nc/$*.nc/' $^ > $@
+
+input: clean
+	./createnetcdf.py
+	$(MAKE) missingdata.flml missingdimension.flml missingvariable.flml incorrectdimension.flml incorrectvariable.flml
+	../../bin/create_aligned_mesh mesh 2.0 2.0 1.0 10 10 1 -- -1.0 -1.0 -1.0
+
+clean:
+	rm -f *.{nc,pyc,ele,node,face,vtu,stat,s,err} matrixdump* fluidity.*
+	rm -f missing{data,dimension,variable}.flml incorrect{dimension,variable}.flml
+

=== added file 'tests/netcdf_read_errors/createnetcdf.py'
--- tests/netcdf_read_errors/createnetcdf.py	1970-01-01 00:00:00 +0000
+++ tests/netcdf_read_errors/createnetcdf.py	2012-02-11 17:18:19 +0000
@@ -0,0 +1,86 @@
+#!/usr/bin/env python
+
+from Scientific.IO.NetCDF import NetCDFFile
+from numpy import arange, zeros
+
+def create(missingdata = False, missingdimension = False, missingvariable = False, incorrectdimension = False, incorrectvariable = False):
+
+  if (missingdata):
+    filename = 'missingdata.nc'
+    description = ', with missing data.'
+  elif (missingdimension):
+    filename = 'missingdimension.nc'
+    description = ', with a missing dimension.'
+  elif (missingvariable):
+    filename = 'missingvariable.nc'
+    description = ', with a missing variable.'
+  elif (incorrectdimension):
+    filename = 'incorrectdimension.nc'
+    description = ', with an incorrect dimension label.'
+  elif (incorrectvariable):
+    filename = 'incorrectvariable.nc'
+    description = ', with an incorrect variable label.'
+  else:
+    filename = 'valid.nc'
+    description = '.'
+
+  print "Creating " + filename + description
+
+  f = NetCDFFile(filename, 'w')
+  f.description = 'Example free surface height' + description
+  
+  if (missingdata):
+    offset = -0.4
+  else:
+    offset = 0.0
+
+  x = arange(-1.2 + offset, 1.21 + offset, 0.2)
+  y = arange(-1.2, 1.21, 0.2)
+  h = zeros((len(x),len(y)))
+
+  for i in range(len(x)):
+    for j in range(len(y)):
+      # Nice ordering netCDF API - y,x !
+      h[j,i] = x[i] * y[j]
+
+  if (not incorrectdimension):
+    xdimlabel = 'x'
+  else:
+    xdimlabel = 'lat'
+  ydimlabel = 'y'
+  if (not incorrectvariable):
+    zdimlabel = 'z'
+  else:
+    zdimlabel = 'height'
+
+  print xdimlabel, ydimlabel, zdimlabel
+
+  # dimensions
+  f.createDimension(xdimlabel, len(x))
+  if (not missingdimension):
+    f.createDimension(ydimlabel, len(y))
+
+  # variables
+  fx = f.createVariable(xdimlabel, 'd', (xdimlabel,))
+  fx[:] = x
+  if (not missingdimension):
+    fy = f.createVariable(ydimlabel, 'd', (ydimlabel,))
+    fy[:] = y
+    if (not missingvariable):
+      fh = f.createVariable(zdimlabel, 'd', (xdimlabel, ydimlabel,))
+      fh[:] = h
+  else:
+    fh = f.createVariable(zdimlabel, 'd', (xdimlabel,))
+    fh[:] = x
+
+  f.close()
+
+create()
+create(missingdata=True)
+create(missingdimension=True)
+create(missingvariable=True)
+create(incorrectdimension=True)
+create(incorrectvariable=True)
+
+
+

=== added file 'tests/netcdf_read_errors/missingdata.err.sample'
--- tests/netcdf_read_errors/missingdata.err.sample	1970-01-01 00:00:00 +0000
+++ tests/netcdf_read_errors/missingdata.err.sample	2012-02-11 17:18:19 +0000
@@ -0,0 +1,44 @@
+Point [0.8, -1] not found in netCDF file, missingdata.nc.
+Point [1, -1] not found in netCDF file, missingdata.nc.
+Point [0.8, -0.8] not found in netCDF file, missingdata.nc.
+Point [1, -0.8] not found in netCDF file, missingdata.nc.
+Point [0.8, -0.6] not found in netCDF file, missingdata.nc.
+Point [1, -0.6] not found in netCDF file, missingdata.nc.
+Point [0.8, -0.4] not found in netCDF file, missingdata.nc.
+Point [1, -0.4] not found in netCDF file, missingdata.nc.
+Point [0.8, -0.2] not found in netCDF file, missingdata.nc.
+Point [1, -0.2] not found in netCDF file, missingdata.nc.
+Point [0.8, 0] not found in netCDF file, missingdata.nc.
+Point [1, 0] not found in netCDF file, missingdata.nc.
+Point [0.8, 0.2] not found in netCDF file, missingdata.nc.
+Point [1, 0.2] not found in netCDF file, missingdata.nc.
+Point [0.8, 0.4] not found in netCDF file, missingdata.nc.
+Point [1, 0.4] not found in netCDF file, missingdata.nc.
+Point [0.8, 0.6] not found in netCDF file, missingdata.nc.
+Point [1, 0.6] not found in netCDF file, missingdata.nc.
+Point [0.8, 0.8] not found in netCDF file, missingdata.nc.
+Point [1, 0.8] not found in netCDF file, missingdata.nc.
+Point [0.8, 1] not found in netCDF file, missingdata.nc.
+Point [1, 1] not found in netCDF file, missingdata.nc.
+Point [0.8, -1] not found in netCDF file, missingdata.nc.
+Point [1, -1] not found in netCDF file, missingdata.nc.
+Point [0.8, -0.8] not found in netCDF file, missingdata.nc.
+Point [1, -0.8] not found in netCDF file, missingdata.nc.
+Point [0.8, -0.6] not found in netCDF file, missingdata.nc.
+Point [1, -0.6] not found in netCDF file, missingdata.nc.
+Point [0.8, -0.4] not found in netCDF file, missingdata.nc.
+Point [1, -0.4] not found in netCDF file, missingdata.nc.
+Point [0.8, -0.2] not found in netCDF file, missingdata.nc.
+Point [1, -0.2] not found in netCDF file, missingdata.nc.
+Point [0.8, 0] not found in netCDF file, missingdata.nc.
+Point [1, 0] not found in netCDF file, missingdata.nc.
+Point [0.8, 0.2] not found in netCDF file, missingdata.nc.
+Point [1, 0.2] not found in netCDF file, missingdata.nc.
+Point [0.8, 0.4] not found in netCDF file, missingdata.nc.
+Point [1, 0.4] not found in netCDF file, missingdata.nc.
+Point [0.8, 0.6] not found in netCDF file, missingdata.nc.
+Point [1, 0.6] not found in netCDF file, missingdata.nc.
+Point [0.8, 0.8] not found in netCDF file, missingdata.nc.
+Point [1, 0.8] not found in netCDF file, missingdata.nc.
+Point [0.8, 1] not found in netCDF file, missingdata.nc.
+Point [1, 1] not found in netCDF file, missingdata.nc.

=== added file 'tests/netcdf_read_errors/netcdf_read_errors.xml'
--- tests/netcdf_read_errors/netcdf_read_errors.xml	1970-01-01 00:00:00 +0000
+++ tests/netcdf_read_errors/netcdf_read_errors.xml	2012-02-11 17:18:19 +0000
@@ -0,0 +1,130 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<testproblem>
+  <name>netcdf_read_errors</name>
+  <owner userid="asc"/>
+  <tags>flml netcdf</tags>
+  <problem_definition length="medium" nprocs="1">
+    <command_line>for name in valid missing{data,dimension,variable} incorrect{dimension,variable}; do echo "Running case ${name}"; ../../bin/fluidity "${name}.flml" 1>/dev/null 2>"${name}.err"; done</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="valid_err_filesize" language="python">
+import os
+valid_err_filesize = os.stat('valid.err').st_size
+    </variable>
+    <variable name="missingdata_err_linenumber" language="python">
+def filelines(filename):
+  with open(filename) as f:
+    for i, l in enumerate(f):
+      pass
+  return i + 1
+missingdata_err_linenumber = filelines('missingdata.err')
+    </variable>
+    <variable name="missingdimension_err_linenumber" language="python">
+def filelines(filename):
+  with open(filename) as f:
+    for i, l in enumerate(f):
+      pass
+  return i + 1
+missingdimension_err_linenumber = filelines('missingdimension.err')
+    </variable>
+    <variable name="missingvariable_err_linenumber" language="python">
+def filelines(filename):
+  with open(filename) as f:
+    for i, l in enumerate(f):
+      pass
+  return i + 1
+missingvariable_err_linenumber = filelines('missingvariable.err')
+    </variable>
+    <variable name="incorrectdimension_err_linenumber" language="python">
+def filelines(filename):
+  with open(filename) as f:
+    for i, l in enumerate(f):
+      pass
+  return i + 1
+incorrectdimension_err_linenumber = filelines('incorrectdimension.err')
+    </variable>
+    <variable name="incorrectvariable_err_linenumber" language="python">
+def filelines(filename):
+  with open(filename) as f:
+    for i, l in enumerate(f):
+      pass
+  return i + 1
+incorrectvariable_err_linenumber = filelines('incorrectvariable.err')
+    </variable>
+    <variable name="missingdata_err" language="python">
+fileHandle = open('missingdata.err.sample')
+correct = fileHandle.read()
+fileHandle.close()
+fileHandle = open('missingdata.err')
+line = fileHandle.read()
+fileHandle.close()
+missingdata_err = ( line == correct )
+    </variable>
+    <variable name="missingdimension_err" language="python">
+correct = 'ncdimid: ncid 65536: NetCDF: Invalid dimension ID or name\n'
+fileHandle = open('missingdimension.err')
+line = fileHandle.read()
+fileHandle.close()
+missingdimension_err = ( line == correct )
+    </variable>
+    <variable name="missingvariable_err" language="python">
+correct = 'ncvarid: ncid 65536: NetCDF: Variable not found\n'
+fileHandle = open('missingvariable.err')
+line = fileHandle.read()
+fileHandle.close()
+missingvariable_err = ( line == correct )
+    </variable>
+    <variable name="incorrectdimension_err" language="python">
+correct = 'ncdimid: ncid 65536: NetCDF: Invalid dimension ID or name\n'
+fileHandle = open('incorrectdimension.err')
+line = fileHandle.read()
+fileHandle.close()
+incorrectdimension_err = ( line == correct )
+    </variable>
+    <variable name="incorrectvariable_err" language="python">
+correct = 'ncvarid: ncid 65536: NetCDF: Variable not found\n'
+fileHandle = open('incorrectvariable.err')
+line = fileHandle.read()
+fileHandle.close()
+incorrectvariable_err = ( line == correct )
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="Checking the valid case has reported no errors" language="python">
+assert ( valid_err_filesize == 0 )
+    </test>
+    <test name="Checking the case with missing data has reported exactly 44 errors" language="python">
+assert ( missingdata_err_linenumber == 44 )
+    </test>
+    <test name="Checking the case with a missing dimension has reported exactly 1 error" language="python">
+assert ( missingdimension_err_linenumber == 1 )
+    </test>
+    <test name="Checking the case with a missing variable has reported exactly 1 error" language="python">
+assert ( missingvariable_err_linenumber == 1 )
+    </test>
+    <test name="Checking the case with an incorrect dimension label has reported exactly 1 error" language="python">
+assert ( incorrectdimension_err_linenumber == 1 )
+    </test>
+    <test name="Checking the case with an incorrect variable label has reported exactly 1 error" language="python">
+assert ( incorrectvariable_err_linenumber == 1 )
+    </test>
+    <test name="Checking the case with missing data has reported the correct errors" language="python">
+assert ( missingdata_err )
+    </test>
+    <test name="Checking the case with a missing dimension has reported the correct error" language="python">
+assert ( missingdimension_err )
+    </test>
+    <test name="Checking the case with a missing variable has reported the correct error" language="python">
+assert ( missingvariable_err )
+    </test>
+    <test name="Checking the case with an incorrect dimension label has reported the correct error" language="python">
+assert ( incorrectdimension_err )
+    </test>
+    <test name="Checking the case with an incorrect variable label has reported the correct error" language="python">
+assert ( incorrectvariable_err )
+    </test>
+  </pass_tests>
+  <warn_tests>
+  </warn_tests>
+</testproblem>
+

=== added file 'tests/netcdf_read_errors/valid.flml'
--- tests/netcdf_read_errors/valid.flml	1970-01-01 00:00:00 +0000
+++ tests/netcdf_read_errors/valid.flml	2012-02-11 17:18:19 +0000
@@ -0,0 +1,292 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">netcdf_read_free_surface</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">3</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="mesh">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">2</integer_value>
+      </degree>
+    </quadrature>
+    <ocean_boundaries>
+      <top_surface_ids>
+        <integer_value shape="1" rank="1">2</integer_value>
+      </top_surface_ids>
+      <bottom_surface_ids>
+        <integer_value shape="1" rank="1">1</integer_value>
+      </bottom_surface_ids>
+      <scalar_field name="DistanceToTop" rank="0">
+        <diagnostic>
+          <algorithm name="Internal" material_phase_support="multiple"/>
+          <mesh name="CoordinateMesh"/>
+          <output/>
+          <stat/>
+          <convergence>
+            <include_in_convergence/>
+          </convergence>
+          <detectors>
+            <include_in_detectors/>
+          </detectors>
+          <steady_state>
+            <include_in_steady_state/>
+          </steady_state>
+        </diagnostic>
+      </scalar_field>
+      <scalar_field name="DistanceToBottom" rank="0">
+        <diagnostic>
+          <algorithm name="Internal" material_phase_support="multiple"/>
+          <mesh name="CoordinateMesh"/>
+          <output/>
+          <stat/>
+          <convergence>
+            <include_in_convergence/>
+          </convergence>
+          <detectors>
+            <include_in_detectors/>
+          </detectors>
+          <steady_state>
+            <include_in_steady_state/>
+          </steady_state>
+        </diagnostic>
+      </scalar_field>
+    </ocean_boundaries>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">1.0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="VelocityMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">2.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+  </timestepping>
+  <physical_parameters>
+    <gravity>
+      <magnitude>
+        <real_value rank="0">9.81</real_value>
+      </magnitude>
+      <vector_field name="GravityDirection" rank="1">
+        <prescribed>
+          <mesh name="CoordinateMesh"/>
+          <value name="WholeMesh">
+            <constant>
+              <real_value shape="3" dim1="dim" rank="1">0.0 0.0 -1.0</real_value>
+            </constant>
+          </value>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </vector_field>
+    </gravity>
+  </physical_parameters>
+  <material_phase name="fluid">
+    <equation_of_state>
+      <fluids>
+        <linear>
+          <reference_density>
+            <real_value rank="0">1.0</real_value>
+          </reference_density>
+          <subtract_out_hydrostatic_level/>
+        </linear>
+      </fluids>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <integrate_continuity_by_parts/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">10</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <free_surface>
+            <from_netcdf file_name="valid.nc">
+              <format>
+                <string_value>raw</string_value>
+              </format>
+            </from_netcdf>
+          </free_surface>
+        </initial_condition>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <stabilisation>
+              <no_stabilisation/>
+            </stabilisation>
+            <mass_terms>
+              <lump_mass_matrix/>
+            </mass_terms>
+            <advection_terms/>
+            <stress_terms>
+              <tensor_form/>
+            </stress_terms>
+          </continuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">1.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">0.5</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0E-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">10</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="3" dim1="dim" rank="1">0.0 0.0 0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="top">
+          <surface_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </surface_ids>
+          <type name="free_surface"/>
+        </boundary_conditions>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="FreeSurface" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="PressureMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+  </material_phase>
+</fluidity_options>

=== added directory 'tests/netcdf_read_free_surface'
=== added file 'tests/netcdf_read_free_surface/Makefile'
--- tests/netcdf_read_free_surface/Makefile	1970-01-01 00:00:00 +0000
+++ tests/netcdf_read_free_surface/Makefile	2012-02-11 17:18:19 +0000
@@ -0,0 +1,13 @@
+
+default: run
+
+run: input
+	../../bin/fluidity -v 3 -l netcdf_read_free_surface.flml 
+
+input: clean
+	./createnetcdf.py
+	../../bin/create_aligned_mesh mesh 2.0 2.0 1.0 10 10 1 -- -1.0 -1.0 -1.0
+
+clean:
+	rm -f height.nc *.pyc *.ele *.node *.face *.vtu *.stat *.s matrixdump* fluidity.*
+

=== added file 'tests/netcdf_read_free_surface/createnetcdf.py'
--- tests/netcdf_read_free_surface/createnetcdf.py	1970-01-01 00:00:00 +0000
+++ tests/netcdf_read_free_surface/createnetcdf.py	2012-02-11 17:18:19 +0000
@@ -0,0 +1,33 @@
+#!/usr/bin/env python
+
+from Scientific.IO.NetCDF import NetCDFFile
+from numpy import arange, zeros
+import height
+
+f = NetCDFFile('height.nc', 'w')
+f.description = 'Example free surface height.'
+
+x = arange(-1.2, 1.21, 0.2)
+y = arange(-1.2, 1.21, 0.2)
+h = zeros((len(x),len(y)))
+
+for i in range(len(x)):
+  for j in range(len(y)):
+    # Nice ordering netCDF API - y,x !
+    h[j,i] = height.function([x[i], y[j]])
+
+# dimensions
+f.createDimension('x', len(x))
+f.createDimension('y', len(y))
+
+# variables
+fx = f.createVariable('x', 'd', ('x',))
+fy = f.createVariable('y', 'd', ('y',))
+fh = f.createVariable('z', 'd', ('x', 'y',))
+
+fx[:] = x
+fy[:] = y
+fh[:] = h
+
+f.close()
+

=== added file 'tests/netcdf_read_free_surface/height.py'
--- tests/netcdf_read_free_surface/height.py	1970-01-01 00:00:00 +0000
+++ tests/netcdf_read_free_surface/height.py	2012-02-11 17:18:19 +0000
@@ -0,0 +1,13 @@
+#!/usr/bin/env python
+
+from matplotlib.mlab import bivariate_normal
+
+def function(position):
+  xg = 3.0 * position[0]
+  yg = 3.0 * position[1]
+  h1 = bivariate_normal(xg, yg, 1.0, 1.0, 0.0, 0.0)
+  h2 = bivariate_normal(xg, yg, 1.5, 0.5, 1, 1)
+  h = 10.0 * (h1 - h2)
+  return h
+
+

=== added file 'tests/netcdf_read_free_surface/netcdf_read_free_surface.flml'
--- tests/netcdf_read_free_surface/netcdf_read_free_surface.flml	1970-01-01 00:00:00 +0000
+++ tests/netcdf_read_free_surface/netcdf_read_free_surface.flml	2012-02-11 17:18:19 +0000
@@ -0,0 +1,309 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">netcdf_read_free_surface</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">fluids</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">3</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_file file_name="mesh">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">2</integer_value>
+      </degree>
+    </quadrature>
+    <ocean_boundaries>
+      <top_surface_ids>
+        <integer_value shape="1" rank="1">2</integer_value>
+      </top_surface_ids>
+      <bottom_surface_ids>
+        <integer_value shape="1" rank="1">1</integer_value>
+      </bottom_surface_ids>
+      <scalar_field name="DistanceToTop" rank="0">
+        <diagnostic>
+          <algorithm name="Internal" material_phase_support="multiple"/>
+          <mesh name="CoordinateMesh"/>
+          <output/>
+          <stat/>
+          <convergence>
+            <include_in_convergence/>
+          </convergence>
+          <detectors>
+            <include_in_detectors/>
+          </detectors>
+          <steady_state>
+            <include_in_steady_state/>
+          </steady_state>
+        </diagnostic>
+      </scalar_field>
+      <scalar_field name="DistanceToBottom" rank="0">
+        <diagnostic>
+          <algorithm name="Internal" material_phase_support="multiple"/>
+          <mesh name="CoordinateMesh"/>
+          <output/>
+          <stat/>
+          <convergence>
+            <include_in_convergence/>
+          </convergence>
+          <detectors>
+            <include_in_detectors/>
+          </detectors>
+          <steady_state>
+            <include_in_steady_state/>
+          </steady_state>
+        </diagnostic>
+      </scalar_field>
+    </ocean_boundaries>
+  </geometry>
+  <io>
+    <dump_format>
+      <string_value>vtk</string_value>
+    </dump_format>
+    <dump_period>
+      <constant>
+        <real_value rank="0">1.0</real_value>
+      </constant>
+    </dump_period>
+    <output_mesh name="VelocityMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">2.0</real_value>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1.0</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">1.0</real_value>
+    </finish_time>
+  </timestepping>
+  <physical_parameters>
+    <gravity>
+      <magnitude>
+        <real_value rank="0">9.81</real_value>
+      </magnitude>
+      <vector_field name="GravityDirection" rank="1">
+        <prescribed>
+          <mesh name="CoordinateMesh"/>
+          <value name="WholeMesh">
+            <constant>
+              <real_value shape="3" dim1="dim" rank="1">0.0 0.0 -1.0</real_value>
+            </constant>
+          </value>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </vector_field>
+    </gravity>
+  </physical_parameters>
+  <material_phase name="fluid">
+    <equation_of_state>
+      <fluids>
+        <linear>
+          <reference_density>
+            <real_value rank="0">1.0</real_value>
+          </reference_density>
+          <subtract_out_hydrostatic_level/>
+        </linear>
+      </fluids>
+    </equation_of_state>
+    <scalar_field name="Pressure" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <integrate_continuity_by_parts/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">never</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">10</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <free_surface>
+            <from_netcdf file_name="height.nc">
+              <format>
+                <string_value>raw</string_value>
+              </format>
+            </from_netcdf>
+          </free_surface>
+        </initial_condition>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <no_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="LinearMomentum"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <stabilisation>
+              <no_stabilisation/>
+            </stabilisation>
+            <mass_terms>
+              <lump_mass_matrix/>
+            </mass_terms>
+            <advection_terms/>
+            <stress_terms>
+              <tensor_form/>
+            </stress_terms>
+          </continuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">1.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">0.5</real_value>
+          </theta>
+          <relaxation>
+            <real_value rank="0">0</real_value>
+          </relaxation>
+        </temporal_discretisation>
+        <solver>
+          <iterative_method name="gmres">
+            <restart>
+              <integer_value rank="0">30</integer_value>
+            </restart>
+          </iterative_method>
+          <preconditioner name="sor"/>
+          <relative_error>
+            <real_value rank="0">1.0E-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">10</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value shape="3" dim1="dim" rank="1">0.0 0.0 0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="top">
+          <surface_ids>
+            <integer_value shape="1" rank="1">2</integer_value>
+          </surface_ids>
+          <type name="free_surface"/>
+        </boundary_conditions>
+        <output/>
+        <stat>
+          <include_in_stat/>
+          <previous_time_step>
+            <exclude_from_stat/>
+          </previous_time_step>
+          <nonlinear_field>
+            <exclude_from_stat/>
+          </nonlinear_field>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </vector_field>
+    <scalar_field name="FreeSurface" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="PressureMesh"/>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </scalar_field>
+    <scalar_field name="AnalyticFreeSurface" rank="0">
+      <prescribed>
+        <mesh name="VelocityMesh"/>
+        <value name="WholeMesh">
+          <python>
+            <string_value lines="20" type="python">def val(X,t):
+  import height
+  return height.function([X[0],X[1]])</string_value>
+          </python>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+  </material_phase>
+</fluidity_options>

=== added file 'tests/netcdf_read_free_surface/netcdf_read_free_surface.xml'
--- tests/netcdf_read_free_surface/netcdf_read_free_surface.xml	1970-01-01 00:00:00 +0000
+++ tests/netcdf_read_free_surface/netcdf_read_free_surface.xml	2012-02-11 17:18:19 +0000
@@ -0,0 +1,49 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<testproblem>
+  <name>netcdf_read_free_surface</name>
+  <owner userid="asc"/>
+  <tags>flml netcdf</tags>
+  <problem_definition length="medium" nprocs="1">
+    <command_line>fluidity -v3 -l netcdf_read_free_surface.flml</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="max_freesurface_error" language="python">
+import vtktools
+
+# Read output file
+u = vtktools.vtu("netcdf_read_free_surface_0.vtu")
+xyz = u.GetLocations()
+f  = u.GetField("FreeSurface")
+af = u.GetField("AnalyticFreeSurface")
+err = f - af
+
+#xmin = min(xyz[:,0])
+#xmax = max(xyz[:,0])
+#ymin = min(xyz[:,1])
+#ymax = max(xyz[:,1])
+#ztop = max(xyz[:,2])
+#
+## Focus on top surface interior
+#for i in range(len(err)):
+#  # Free surface not currently initialised correctly on bottom - investigate
+#  if (xyz[i,1] != ztop):
+#    err[i] = 0.0
+#  # Free surface not currently initialised correctly on the boundary - investigate
+#  if ((xyz[i,0] == xmin) or (xyz[i,0] == xmax) or (xyz[i,1] == ymin) or (xyz[i,1] == ymax)):
+#    err[i] = 0.0
+
+# Write out error field for debugging
+u.AddScalarField('Error', err)
+u.Write('error.vtu')
+
+max_freesurface_error = max(abs(err))
+    </variable>
+  </variables>
+  <pass_tests>
+    <test name="Checking difference is smaller than 1e-10" language="python">
+assert ( max_freesurface_error &lt; 1e-10 )
+    </test>
+  </pass_tests>
+  <warn_tests>
+  </warn_tests>
+</testproblem>