Merge lp:~fluidity-core/fluidity/fluidity-initialisation-from_netcdf into lp:fluidity

[Stephan Kramer]

This merge achieves 3 things:
1) Change tidal forcing boundary condition so that the phase is no longer interpolated (this causes issues at the 360-0 line). Instead the tidal elevation is calculated as:
  A(x) cos(omega*t+alpha(x)) = A(x) cos(omega*t) cos(alpha(x)-sin(omega*t) sin(alpha(x))
                             = Re[ A(x) exp(i omega*t) exp(i alpha(x)) ]
where A(x), alpha(x) and omega are the amplitude, phase and frequency of the tidal constituent. We therefore now compute A(x)*cos(alpha(x)) and A(x)*sin(alpha(x)), the real and imaginary components of the complex tidal constituent, at the nodes of the netCDF first, before interpolating to the fluidity mesh.
2) The netCDF file is only read once at the beginning of the run instead of every timestep.
3) Prescribed fields have an option to be prescribed from a netCDF file using the same netCDF reader as used for the tidal boundary forcing, so we can actually see what data is being used. This is done by adding a new "netcdf" format (next to the existing "vtu" format under "from_file"), under which the user can specify which netCDF variable to read.

The reason these 3 things are combined in this merge is to allow testing of the tidal boundary forcing - currently no such tests exist. This merge introduces the sample_netcdf_test/ that tests
both the SampleNetCDF reader and the actual tidal boundary condition values. This is done by creating a netCDF file using Scientific python with variables set via an analytical function and checking the same values are read in fluidity.

Smaller changes:
- taking away the dirichlet/from_file/tidal option under various balance pressure fields. For these cases (other than field%name=="Pressure") the tidal forcing was never multiplied by g, so this would not have worked.
- because this was also still present under the no longer existing "prognostic free surface", this has now finally been taken out of the schema - which meant also deleting it from the shallow water and test_advection_diffussion schemes.
- fix intents of id argument of most SampleNetCDF_XXX() routines.
- add SampleNetCDF_to_field() routine that interpolate netCDF data to scalar_field on the fluidity mesh

[lp:~fluidity-core/fluidity/fluidity-initialisation-from_netcdf]

3943. By Stephan Kramer on 2012-03-27

Fix initialising prescribed fields from "netcdf" format in parallel.

[Stephan Kramer (s-kramer)]

I use python code for all of this now, that I'll make available soon. Therefore retracting my merge request for now.

[Adam Candy (asc)]

> I use python code for all of this now, that I'll make available soon.
> Therefore retracting my merge request for now.

How very gentlemanly of you!

I have a branch implementing some of the above, so we could take a look at the two and merge in a combination. The NetCDF input certainly requires an improvement - Jon's also mentioned working on this?

I think my branch with additional debugging for the NetCDf readers (a portion of them) is already in. This does give information on the location of problem reads, warns when values are recast and is accompanied by a suite of tests.

[Jon Hill (jon-hill)]

I have a half-competed branch for reading netcdf files and putting those
data onto a mesh. I was stuck with a unit test not compiling if I
remember correctly.

The branch is at:
bzr+ssh://bazaar.launchpad.net/~fluidity-core/fluidity/refactor-netcdf/

I was trying to create a generic reader which created a mesh using
Delauny triangulation. We can then use Fluidity's mesh interpolation to
put this on any mesh.

Michael Lange is also interested in this, including outputting stat
files and detectors to NetCDF.

Jon

On 09/05/13 12:14, Adam Candy wrote:
>> I use python code for all of this now, that I'll make available soon.
>> Therefore retracting my merge request for now.
>
> How very gentlemanly of you!
>
> I have a branch implementing some of the above, so we could take a look at the two and merge in a combination. The NetCDF input certainly requires an improvement - Jon's also mentioned working on this?
>
> I think my branch with additional debugging for the NetCDf readers (a portion of them) is already in. This does give information on the location of problem reads, warns when values are recast and is accompanied by a suite of tests.
>

--
Dr Jon Hill
Research Fellow
AMCG, Imperial College London
Mob: [phone]

Unmerged revisions

3943. By Stephan Kramer on 2012-03-27

Fix initialising prescribed fields from "netcdf" format in parallel.

3942. By Stephan Kramer on 2012-03-02

Merging from trunk.

3941. By Stephan Kramer on 2012-03-02

Cleaning up before merge request.

3940. By Stephan Kramer on 2012-03-01

Should have committed these file for the test as well.

3939. By Stephan Kramer on 2012-03-01

Merging in from trunk.

3938. By Stephan Kramer on 2012-02-27

Actually testing the values brought up two more bugs. Now pretty confident it actually works.

3937. By Stephan Kramer on 2012-02-27

Add prescribed field initialised from netcdf, to be compared with a python diagnostic field using the same function used to create the netcdf file.
Fixed mem leak in initialsing form netcdf.

3936. By Stephan Kramer on 2012-02-27

Several bug fixes. Test now runs through and seems to give sensible answers.

3935. By Stephan Kramer on 2012-02-27

Start of a new test to test SampleNetCDF and tidal boundary conditions.
Does not yet work.
Accidentily removed tidal boundary conditions also from "Pressure" - putting back in.

3934. By Stephan Kramer on 2012-02-21

This commit address 2 issues:
* interpolating errors along the phase=360 lines of tidal forcing read netcdf files
* the tidal forcing is no longer read every timestep, but only once in populate_state - this prevents reading the entire netcdf tidal atlas on every processor every timestep.

In femtools/SampleNetCDF_fortran.F90:
- fix intents of id argument of most SampleNetCDF_XXX() routines.
- add SampleNetCDF_to_field that sets a scalar field from netcdf values.
In preprocessor/Boundary_Conditions_From_Options.F90
- add new populate_tidal_bc() routine, that initialises the tidal bc calculation by computing the real and complex components of A(x)exp(i*alpha(x)), where A(x) is amplitude and alpha(x) the phase, and store these components as surface fields under the bc.
- in set_tidal_bc_value() it now uses these to compute the free surface contributions of the different tidal components.
In schemas/:
- taking away the dirichlet/from_file/tidal option under various balance pressure fields. For these cases (other than field%name=="Pressure") the tidal forcing was never multiplied by g, so this would not have worked.
- because this was also still present under the no longer existing "prognostic free surface", this has now finally been taken out of the schema - which meant also deleting it from the shallow water and test_advection_diffussion schemes.

=== modified file 'femtools/Makefile.dependencies'
--- femtools/Makefile.dependencies	2012-01-20 15:54:44 +0000
+++ femtools/Makefile.dependencies	2012-03-27 11:21:19 +0000
@@ -933,7 +933,7 @@
 	@true
 
 SampleNetCDF_fortran.o ../include/samplenetcdf.mod: SampleNetCDF_fortran.F90 \
-   ../include/fldebug.mod
+   ../include/coordinates.mod ../include/fields.mod ../include/fldebug.mod
 
 ../include/shape_functions.mod: Shape_Functions.o
 	@true
 
=== modified file 'femtools/SampleNetCDF_fortran.F90'
--- femtools/SampleNetCDF_fortran.F90	2010-11-09 14:08:19 +0000
+++ femtools/SampleNetCDF_fortran.F90	2012-03-27 11:21:19 +0000
@@ -27,6 +27,8 @@
 
 module SampleNetCDF
   use FLDebug
+  use fields
+  use coordinates
   implicit none
 contains
   
@@ -39,13 +41,25 @@
 
   subroutine SampleNetCDF_SetVariable(id, varname)
     character(len=*), intent(in)::varname
-    integer, intent(out)::id
+    integer, intent(in)::id
     
     call samplenetcdf_setvariable_c(id, varname, len_trim(varname))
   end subroutine SampleNetCDF_SetVariable
 
+  subroutine SampleNetCDF_ApplyCos(id)
+    integer, intent(in)::id
+    
+    call samplenetcdf_applycos_c(id)
+  end subroutine SampleNetCDF_ApplyCos
+
+  subroutine SampleNetCDF_ApplySin(id)
+    integer, intent(in)::id
+    
+    call samplenetcdf_applysin_c(id)
+  end subroutine SampleNetCDF_ApplySin
+
   subroutine SampleNetCDF_GetValue(id, longitude, latitude, val)
-    integer, intent(out)::id
+    integer, intent(in)::id
     real, intent(out)::longitude, latitude
     real, intent(out)::val
 
@@ -53,9 +67,27 @@
   end subroutine SampleNetCDF_GetValue
 
   subroutine SampleNetCDF_Close(id)
-    integer, intent(out)::id
+    integer, intent(in)::id
     
     call samplenetcdf_close_c(id)
   end subroutine SampleNetCDF_Close
+
+  subroutine SampleNetCDF_to_field(id, field, bc_position)
+    ! set scalar field from netcdf values
+    integer, intent(in):: id
+    type(scalar_field), intent(inout):: field
+    type(vector_field), intent(in):: bc_position
+
+    integer:: i
+    real:: longitude, latitude, xyz(3), value
+
+    do i=1, node_count(bc_position)
+      xyz = node_val(bc_position, i)
+      call LongitudeLatitude(xyz, longitude, latitude)
+      call SampleNetCDF_GetValue(id, longitude, latitude, value)
+      call set(field, i, value)
+    end do
+
+  end subroutine SampleNetCDF_to_field
   
 end module SampleNetCDF
 
=== modified file 'ocean_forcing/SampleNetCDF.cpp'
--- ocean_forcing/SampleNetCDF.cpp	2012-02-20 17:52:19 +0000
+++ ocean_forcing/SampleNetCDF.cpp	2012-03-27 11:21:19 +0000
@@ -25,6 +25,7 @@
 #include <iostream>
 #include <map>
 #include <vector>
+#include <algorithm>
 
 #include <vtk.h>
 
@@ -65,6 +66,7 @@
   }
 
   data = in.data;
+  read_data = in.read_data;
   verbose = in.verbose;
   
   reader = in.reader;
@@ -189,7 +191,19 @@
   if(verbose)
     cout<<"void SampleNetCDF::SetVariable("<<varname<<")\n";
 
-  reader.Read(varname, data);
+  reader.Read(varname, read_data);
+  data = read_data.begin();
+}
+
+void SampleNetCDF::ApplyFunction(double f(double x)){
+  // apply transformation to read data
+  // transformed data is stored in manipulated_data, allocate it to right size (or reuse)
+  if (manipulated_data.size()!=read_data.size())
+    manipulated_data.resize(read_data.size());
+  // iterate over read_data, apply function f and store in manipulated_data
+  transform(read_data.begin(), read_data.end(), manipulated_data.begin(), f);
+  // from now on we sample the manipulated data:
+  data = manipulated_data.begin();
 }
 
 void SampleNetCDF::VerboseOff(){
@@ -226,6 +240,35 @@
     
     netcdf_sampler[*id].SetVariable(string(varname, *len));
   }
+
+  // auxilary functions cosdeg and sindeg: cos and sin as function of degrees rather than radials
+  double const pi=4*atan(1);
+  double cosdeg(double deg){
+    return cos(pi*deg/180);
+  }
+  double sindeg(double deg){
+    return sin(pi*deg/180);
+  }
+
+#define samplenetcdf_applycos_fc F77_FUNC_(samplenetcdf_applycos_c, SAMPLENETCDF_APPLYCOS_C)
+  void samplenetcdf_applycos_fc(const int *id){
+    if(netcdf_sampler.find(*id)==netcdf_sampler.end()){
+      cerr<<"ERROR: netcdf file has not been opened\n";
+      exit(-1);
+    }
+    
+    netcdf_sampler[*id].ApplyFunction(cosdeg);
+  }
+
+#define samplenetcdf_applysin_fc F77_FUNC_(samplenetcdf_applysin_c, SAMPLENETCDF_APPLYSIN_C)
+  void samplenetcdf_applysin_fc(const int *id){
+    if(netcdf_sampler.find(*id)==netcdf_sampler.end()){
+      cerr<<"ERROR: netcdf file has not been opened\n";
+      exit(-1);
+    }
+    
+    netcdf_sampler[*id].ApplyFunction(sindeg);
+  }
   
 #define samplenetcdf_getvalue_fc F77_FUNC_(samplenetcdf_getvalue_c, SAMPLENETCDF_GETVALUE_C)
   void samplenetcdf_getvalue_fc(const int *id, const double *longitude, const double *latitude, double *val){
 
=== modified file 'ocean_forcing/SampleNetCDF.h'
--- ocean_forcing/SampleNetCDF.h	2009-12-16 14:11:47 +0000
+++ ocean_forcing/SampleNetCDF.h	2012-03-27 11:21:19 +0000
@@ -52,6 +52,7 @@
 
   void SetFile(std::string);
   void SetVariable(std::string);
+  void ApplyFunction(double f(double x));
   
   void VerboseOff();
   void VerboseOn();
@@ -66,7 +67,12 @@
   double spacing[2];          // Grid spacing
   int dimension[2];           // Grid dimensions
   
-  std::vector<double> data;
+  // random access iterator over the data being sampled 
+  // (this may be iterating manipulated_data or read_data):
+  std::vector<double>::iterator data;
+
+  std::vector<double> read_data; // original data in netcdf file
+  std::vector<double> manipulated_data; // data after manipulation
   bool verbose;
 
   NetCDFReader reader;
 
=== modified file 'preprocessor/Boundary_Conditions_From_Options.F90'
--- preprocessor/Boundary_Conditions_From_Options.F90	2011-11-23 13:46:31 +0000
+++ preprocessor/Boundary_Conditions_From_Options.F90	2012-03-27 11:21:19 +0000
@@ -260,6 +260,10 @@
           call insert_surface_field(field, i+1, surface_field)
           call deallocate(surface_field)
 
+          if (have_option(trim(bc_path_i)//"/type::dirichlet/from_file")) then
+            call populate_tidal_bc(field, i+1, bc_position, trim(bc_path_i)//"/type::dirichlet/from_file")
+          end if
+
        case("robin")
 
           call allocate(surface_field, surface_mesh, name="order_zero_coefficient")
@@ -767,7 +771,7 @@
           ! Tidal: set free surface height at the boundary.
           if (have_option(trim(bc_type_path)//"/from_file")) then
             ! Special case for tidal harmonic boundary conditions
-            call set_tidal_bc_value(surface_field, bc_position, trim(bc_type_path), field%name)
+            call set_tidal_bc_value(field, surface_field, i+1, trim(bc_type_path)//"/from_file")
           else if (have_option(trim(bc_type_path)//"/NEMO_data")) then
             call set_nemo_bc_value(state, surface_field, bc_position, trim(bc_type_path), field%name, surface_element_list)
           else
@@ -1092,71 +1096,99 @@
     
   end subroutine set_vector_boundary_conditions_values
   
-  subroutine set_tidal_bc_value(surface_field, bc_position, bc_type_path, field_name)
+  subroutine populate_tidal_bc(field, bc_index, bc_position, bc_option_path)
     ! tidal_forcing - asc
-    type(scalar_field), intent(inout):: surface_field
-    type(vector_field), intent(in):: bc_position
-    character(len=*), intent(in):: bc_type_path, field_name
+    type(scalar_field), intent(inout):: field
+    integer, intent(in):: bc_index
+    type(vector_field), intent(inout):: bc_position
+    ! option path pointing to tidal constituents
+    character(len=*), intent(in):: bc_option_path
 
-    real :: current_time, frequency, amplitude_factor
+    type(scalar_field):: amplitude, real_component, imag_component
     integer :: constituent_count, i, j, id, stat
     character(len=3) :: constituent_name
     character(len=4096) :: file_name, variable_name_amplitude, variable_name_phase
-    real, dimension(:), allocatable::amplitude, phase
-    real :: xyz(3), longitude, latitude
-    real :: gravty
-
-    allocate(amplitude(node_count(surface_field)), phase(node_count(surface_field)))
-    
-    if(have_option(trim(bc_type_path)//"/from_file/tidal")) then
-       call set(surface_field, 0.0)
-       call get_option("/timestepping/current_time", current_time)
-       constituent_count = option_count(trim(bc_type_path)//"/from_file/tidal")
-
-       do i=0, constituent_count-1
-          amplitude = 0.0
-          phase = 0.0
-          call get_option(trim(bc_type_path)//"/from_file/tidal["//int2str(i)//"/amplitude_factor", amplitude_factor, stat=stat)
-          if (stat/=0) then
-             amplitude_factor = 1.0
-          end if
-          
-          ! Taken from E.W. Schwiderski - Rev. Geophys. Space Phys. Vol. 18 No. 1 pp. 243--268, 1980
-          call get_option(trim(bc_type_path)//"/from_file/tidal["//int2str(i)//"]/name", constituent_name, stat=stat)
-
-          frequency = get_tidal_frequency(constituent_name)
-
-          call get_option(trim(bc_type_path)//"/from_file/tidal["//int2str(i)//"]/file_name", file_name, stat=stat)
-          call get_option(trim(bc_type_path)//"/from_file/tidal["//int2str(i)//"]/variable_name_amplitude", variable_name_amplitude, stat=stat)
-          call get_option(trim(bc_type_path)//"/from_file/tidal["//int2str(i)//"]/variable_name_phase", variable_name_phase, stat=stat)
-
-          call SampleNetCDF_Open(trim(file_name), id)
-          call SampleNetCDF_SetVariable(id, trim(variable_name_amplitude))
-          do j=1, node_count(bc_position)
-             xyz = node_val(bc_position, j)
-             call LongitudeLatitude(xyz, longitude, latitude)
-             call SampleNetCDF_GetValue(id, longitude, latitude, amplitude(j))
-          end do
-          call SampleNetCDF_SetVariable(id, trim(variable_name_phase))
-          do j=1, node_count(bc_position)
-             xyz = node_val(bc_position, j)
-             call LongitudeLatitude(xyz, longitude, latitude)
-             call SampleNetCDF_GetValue(id, longitude, latitude, phase(j))
-          end do
-          
-          call get_option('/physical_parameters/gravity/magnitude', gravty)
-          do j=1, node_count(bc_position)
-            if (field_name=="Pressure") then
-              call addto(surface_field, j, &
-                  gravty * amplitude(j) * cos( frequency * current_time - ( phase(j) * pi / 180 ) ))
-            else
-              call addto(surface_field, j, &
-                  amplitude(j) * cos( frequency * current_time - ( phase(j) * pi / 180 ) ))
-            end if
-          end do
-          call SampleNetCDF_Close(id)
-       end do
-    end if
+
+    constituent_count = option_count(trim(bc_option_path)//"/tidal")
+
+    do i=0, constituent_count-1
+
+      ! Taken from E.W. Schwiderski - Rev. Geophys. Space Phys. Vol. 18 No. 1 pp. 243--268, 1980
+      call get_option(trim(bc_option_path)//"/tidal["//int2str(i)//"]/name", constituent_name)
+
+      call get_option(trim(bc_option_path)//"/tidal["//int2str(i)//"]/file_name", file_name)
+      call get_option(trim(bc_option_path)//"/tidal["//int2str(i)//"]/variable_name_amplitude", variable_name_amplitude)
+      call get_option(trim(bc_option_path)//"/tidal["//int2str(i)//"]/variable_name_phase", variable_name_phase)
+
+      call SampleNetCDF_Open(trim(file_name), id)
+
+      ! sample the tidal amplitude at the boundary nodes
+      call SampleNetCDF_SetVariable(id, trim(variable_name_amplitude))
+      call allocate(amplitude, bc_position%mesh, "amplitude")
+      call SampleNetCDF_to_field(id, amplitude, bc_position)
+
+      ! we compute the real and imaginary component of the tide: 
+      ! A(x) * exp(i*alpha(x)) = A(x)cos(alpha(x)) + i * A(x)sin(alpha(x))
+      ! where A(x) is the amplitude, and alpha(x) the phase
+
+      ! sample cos(alpha(x)) and multiply with A(x) at the nodes
+      call SampleNetCDF_SetVariable(id, trim(variable_name_phase))
+      call SampleNetCDF_ApplyCos(id)
+      call allocate(real_component, bc_position%mesh, trim(constituent_name)//"real")
+      call SampleNetCDF_to_field(id, real_component, bc_position)
+      call scale(real_component, amplitude)
+      call insert_surface_field(field, bc_index, real_component)
+      call deallocate(real_component)
+
+      ! sample sin(alpha(x)) and multiply with A(x) at the nodes
+      call SampleNetCDF_ApplySin(id)
+      call allocate(imag_component, bc_position%mesh,  trim(constituent_name)//"imag")
+      call SampleNetCDF_to_field(id, imag_component, bc_position)
+      call scale(imag_component, amplitude)
+      call insert_surface_field(field, bc_index, imag_component)
+      call deallocate(imag_component)
+
+      call deallocate(amplitude)
+      call SampleNetCDF_Close(id)
+    end do
+
+  end subroutine populate_tidal_bc
+
+  subroutine set_tidal_bc_value(field, surface_field, bc_index, bc_option_path)
+    ! tidal components are added into the dirichlet "value" surface_field
+    type(scalar_field), intent(inout):: field, surface_field
+    integer, intent(in):: bc_index
+    ! option path pointing to tidal constituents
+    character(len=*), intent(in):: bc_option_path
+
+    type(scalar_field), pointer:: real_component, imag_component
+    integer :: constituent_count, i
+    character(len=3) :: constituent_name
+    real :: g, t, omega
+
+    call get_option("/timestepping/current_time", t)
+    call get_option('/physical_parameters/gravity/magnitude', g)
+
+    call zero(surface_field)
+
+    constituent_count = option_count(trim(bc_option_path)//"/tidal")
+
+    do i=0, constituent_count-1
+      ! Taken from E.W. Schwiderski - Rev. Geophys. Space Phys. Vol. 18 No. 1 pp. 243--268, 1980
+      call get_option(trim(bc_option_path)//"/tidal["//int2str(i)//"]/name", constituent_name)
+      omega = get_tidal_frequency(constituent_name)
+
+      ! the free surface contribution is calculated as:
+      ! g*A(x) cos(omega * t + alpha(x)) = g*A(x) cos(alpha(x) cos(omega *t) - g*A(x) sin(alpha(x) sin(omega *t)
+
+      real_component => extract_surface_field(field, bc_index, trim(constituent_name)//"real")
+      call addto(surface_field, real_component, scale=g*cos(omega*t))
+
+      imag_component => extract_surface_field(field, bc_index, trim(constituent_name)//"imag")
+      call addto(surface_field, imag_component, scale=-g*sin(omega*t))
+
+    end do
+
   end subroutine set_tidal_bc_value
 
   subroutine set_nemo_bc_value(state, surface_field, bc_position, bc_type_path, field_name, surface_element_list)
 
=== modified file 'preprocessor/Initialise_Fields.F90'
--- preprocessor/Initialise_Fields.F90	2012-02-09 22:34:28 +0000
+++ preprocessor/Initialise_Fields.F90	2012-03-27 11:21:19 +0000
@@ -38,6 +38,8 @@
 use nemo_states_module
 use physics_from_options
 use load_netcdf_module
+use tidal_module
+use SampleNetCDF
 
 implicit none
 
@@ -84,9 +86,11 @@
     character(len=PYTHON_FUNC_LEN) :: func
     real :: current_time
     real :: gravity_magnitude
-    
-    real value
-    integer nid
+
+    ! some variables for reading from netcdf and climatology
+    real, dimension(3):: xyz
+    real:: value, longitude, latitude
+    integer:: id, nid
 
     ! Find out whether initial condition is constant or generated by a 
     ! python function (or comes from something else).    
@@ -123,13 +127,15 @@
          
          call get_option(trim(path) // "/from_file/format/name", format)
          call get_option(trim(path) // "/from_file/file_name", filename)
-         if(isparallel()) then
-           filename = parallel_filename(trim_file_extension(filename), ".vtu")
-         end if
-         ewrite(2, *) "Initialising field " // trim(field%name) // " from file " // trim(filename)
          
          select case (format)
          case ("vtu")
+
+            if(isparallel()) then
+            filename = parallel_filename(trim_file_extension(filename), ".vtu")
+            end if
+            ewrite(2, *) "Initialising field " // trim(field%name) // " from file " // trim(filename)
+
             call get_option(trim(path) // "/from_file/format::vtu/field_name", field_name, default = field%name)
             read_field => vtk_cache_read_scalar_field(filename, field_name)            
             
@@ -164,7 +170,25 @@
 
             call set(field, read_field)
             
+         case ("netcdf")
+
+            ewrite(2, *) "Initialising field " // trim(field%name) // " from file " // trim(filename)
             
+            call get_option(trim(path) // "/from_file/format::netcdf/field_name", field_name)
+            if (field%mesh==position%mesh) then
+              field_position=position
+              call incref(field_position)
+            else
+              ! otherwise the position mapped onto field%mesh
+              field_position = get_remapped_coordinates(position, field%mesh)
+            end if
+
+            call SampleNetCDF_Open(trim(filename), id)
+            call SampleNetCDF_SetVariable(id, trim(field_name))
+            call SampleNetCDF_to_field(id, field, field_position)
+            call SampleNetCDF_Close(id)
+            call deallocate(field_position)
+
          case ("climatology (Boyer2005)")
            
             do nid=1, node_count(position)
 
=== modified file 'preprocessor/Makefile.dependencies'
--- preprocessor/Makefile.dependencies	2011-10-26 19:57:55 +0000
+++ preprocessor/Makefile.dependencies	2012-03-27 11:21:19 +0000
@@ -36,7 +36,8 @@
    ../include/fdebug.h ../include/fields.mod ../include/futils.mod \
    ../include/global_parameters.mod ../include/load_netcdf_module.mod \
    ../include/nemo_states_module.mod ../include/physics_from_options.mod \
-   ../include/spud.mod ../include/tictoc.mod ../include/vtk_cache_module.mod
+   ../include/samplenetcdf.mod ../include/spud.mod ../include/tictoc.mod \
+   ../include/tidal_module.mod ../include/vtk_cache_module.mod
 
 ../include/physics_from_options.mod: Physics_From_Options.o
 	@true
 
=== modified file 'schemas/input_output.rnc'
--- schemas/input_output.rnc	2012-02-11 14:10:31 +0000
+++ schemas/input_output.rnc	2012-03-27 11:21:19 +0000
@@ -35,7 +35,8 @@
       element from_file {
          attribute file_name { xsd:string },
          (
-            vtu_input_format
+            vtu_input_format|
+            netcdf_input_format
          ),
          comment
       }|
@@ -575,9 +576,15 @@
     
 netcdf_input_format =
    (
-      ## NetCDF CF 1.4 (http://cf-pcmdi.llnl.gov/)
+      ## Read from netcdf file
       element format {
-         attribute name { "NetCDF - CF 1.x" },
+         attribute name { "netcdf" },
+         ## The field to read from the NetCDF file.
+         element field_name {
+            anystring,
+            comment
+         },
          comment
       }
    )
+
 
=== modified file 'schemas/input_output.rng'
--- schemas/input_output.rng	2012-02-11 14:10:31 +0000
+++ schemas/input_output.rng	2012-03-27 11:21:19 +0000
@@ -60,7 +60,10 @@
         <attribute name="file_name">
           <data type="string"/>
         </attribute>
-        <ref name="vtu_input_format"/>
+        <choice>
+          <ref name="vtu_input_format"/>
+          <ref name="netcdf_input_format"/>
+        </choice>
         <ref name="comment"/>
       </element>
       <element name="NEMO_data">
@@ -638,10 +641,15 @@
   </define>
   <define name="netcdf_input_format">
     <element name="format">
-      <a:documentation>NetCDF CF 1.4 (http://cf-pcmdi.llnl.gov/)</a:documentation>
+      <a:documentation>Read from netcdf file</a:documentation>
       <attribute name="name">
-        <value>NetCDF - CF 1.x</value>
+        <value>netcdf</value>
       </attribute>
+      <element name="field_name">
+        <a:documentation>The field to read from the NetCDF file.</a:documentation>
+        <ref name="anystring"/>
+        <ref name="comment"/>
+      </element>
       <ref name="comment"/>
     </element>
   </define>
 
=== modified file 'schemas/prognostic_field_options.rnc'
--- schemas/prognostic_field_options.rnc	2012-02-13 15:35:39 +0000
+++ schemas/prognostic_field_options.rnc	2012-03-27 11:21:19 +0000
@@ -1479,7 +1479,7 @@
          ## Type
          element type {
             attribute name { "dirichlet" },
-            input_choice_real_plus_boundary_forcing
+            input_choice_real_contents
          }
       }*,
       prognostic_scalar_output_options,
@@ -1511,7 +1511,7 @@
          ## Type
          element type {
             attribute name { "dirichlet" },
-            input_choice_real_plus_boundary_forcing
+            input_choice_real_contents
          }
       }+,      
       ## Disables checkpointing of this field
@@ -1681,198 +1681,6 @@
       interpolation_algorithm_scalar_full,
       discrete_properties_algorithm_scalar?
    )
-   
-# free surface field, this is a copy of prognostic_scalar_field above
-# removing all options that don't apply (mainly advection related)
-prognostic_free_surface_field =
-   (
-      (
-         ## Spatial discretisation options
-         element spatial_discretisation {
-            ## Form a full 3D system for the free surface
-            element free_surface_3D {
-               empty
-            }?,
-            element fourth_order_dissipation {
-               empty
-            }?,
-            ## low order (linear) free surface
-            element low_order_free_surface {
-               empty
-            }?,
-            (
-               ## Select free surface filter
-               ##
-               ## With PN-PN we need some filter to supress spurious modes.
-               element default_free_surface_filter {
-                  empty
-               }|
-               ## Select free surface filter
-               ##
-               ## With PN-PN we need some filter to supress spurious modes.
-               element user_specified_free_surface_filter {
-                  ## Default is to apply 0.01 and for wetting and drying cases 1.0
-                  element non_linear_filter_coefficient {
-                     real
-                  }
-               }|
-               ## Switch off free surface filter, this is more efficient than setting the coefficient to 0.
-               element switch_off_free_surface_filter {
-                  empty
-               }
-            ),
-            ## Apply wetting and drying routines
-            element wetting_drying {
-               empty
-            }?,
-            ## Tidal forcing options 
-            element tidal_forcing {
-               ## M2
-               element M2 {
-                 empty
-               }?,
-               ## S2
-               element S2 {
-                  empty
-               }?,
-               ## N2
-               element N2 {
-                  empty
-               }?,
-               ## K2
-               element K2 {
-                  empty
-               }?,
-               ## K1
-               element K1 {
-                  empty
-               }?,
-               ## O1
-               element O1 {
-                  empty
-               }?,
-               ## P1
-               element P1 {
-                  empty
-               }?,
-               ## Q1
-               element Q1 {
-                  empty
-               }?,             
-               ## Switch on all tidal components
-               element all_tidal_components {
-                  empty
-               }?,
-               ## Switches on a Love number of 0.3
-               element love_number {
-                  empty
-               }?,
-               ## Use static tidal force for testing
-               element static_tidal_force {
-                  empty
-               }?
-            }?
-         }
-      ),
-      # atheta, ctheta and fstheta (absorption, coriolis and free surface)
-      # need to go in temporal discretisation
-      # they are currently hard-coded however
-      element temporal_discretisation {
-         ## Implicit/explicitness for the free surface.
-         ##
-         ## Suggested value 1.0 (should be at least bigger than 0.5).
-         ##  =0.  -- explicit
-         ##  =0.5 -- Crank-Nicolson
-         ##  =1.  -- implicit
-         element theta {
-            real
-         },
-         # Maybe this should go under a proper absorption field under free surface?
-         ## Implicit/explicitness for absorption
-         ## =0.  -- explicit (default)
-         ## =0.5 -- Crank-Nicolson
-         ## =1.  -- implicit
-         element absorption_theta {
-            real
-         }?
-      },
-      ## Solver
-      element solver {
-         linear_solver_options_sym
-      },
-      (
-         ## Initial condition for WholeMesh
-         ##
-         ## Only specify one condition if not using mesh regions.
-         ## Otherwise select other initial_condition option, specify region_ids
-         ## and distinct names.  Then add extra intial conditions for other regions.
-         element initial_condition {
-            attribute name { "WholeMesh" },
-            input_choice_initial_condition_real
-         }|
-         ## Multiple initial_conditions are allowed if specifying
-         ## different values in different
-         ## regions of the mesh (defined by region_ids).  In this case
-         ## each initial_condition
-         ## requires a distinct name for the options dictionary.
-         element initial_condition {
-            attribute name { string },
-            region_ids?,
-            input_choice_initial_condition_real
-         }
-      )+,
-      ## Boundary conditions
-      element boundary_conditions {
-         attribute name { string },
-         ## Surface id:
-         element surface_ids {
-            integer_vector
-         },
-         ## Type
-         (
-            element type {
-               attribute name { "dirichlet" },
-               (
-                  input_choice_real_contents|
-                  element from_file {
-                     element tidal {
-                        attribute file_name { string },
-                        attribute variable_name_amplitude { string },
-                        attribute variable_name_phase { string },
-                        ## See E.W. Schwiderski - Rev. Geophys. Space
-                        ## Phys. Vol. 18 No. 1 pp. 243--268, 1980
-                        ## for details of these constituent.
-                        attribute name {"M2"|"S2"|"N2"|"K2"|"K1"|"O1"|"P1"|"Q1"|"Mf"|"Mm"|"Ssa"}
-                     }+
-                  }|
-                  ## Set the boundary free-surface height from NEMO data.
-                  ## A prescribed NEMO pressure scalar field must be set to use this option.
-                  ## Set the name of the prescribed NEMO pressure scalar field below.
-                  element NEMO_data {
-                     attribute field_name { string }
-                  }
-               )
-            }|
-            element type {
-               attribute name { "neumann" },
-               input_choice_real
-            }|
-            robin_bc_scalar
-         )
-      }*,
-      # no Diffusivity for field
-      # no source term
-      # no Absorption term
-      # no Adaptive timestepping option
-      prognostic_scalar_output_options,
-      prognostic_scalar_stat_options,
-      scalar_convergence_options,
-      prognostic_detector_options,
-      scalar_steady_state_options,
-      adaptivity_options_prognostic_scalar_field,
-      interpolation_algorithm_scalar,
-      discrete_properties_algorithm_scalar?
-   )
 
 # stream function, this is a copy of prognostic_scalar_field above
 # removing all options that don't apply (mainly advection related)
 
=== modified file 'schemas/prognostic_field_options.rng'
--- schemas/prognostic_field_options.rng	2012-02-13 15:35:39 +0000
+++ schemas/prognostic_field_options.rng	2012-03-27 11:21:19 +0000
@@ -1807,7 +1807,7 @@
           <attribute name="name">
             <value>dirichlet</value>
           </attribute>
-          <ref name="input_choice_real_plus_boundary_forcing"/>
+          <ref name="input_choice_real_contents"/>
         </element>
       </element>
     </zeroOrMore>
@@ -1847,7 +1847,7 @@
           <attribute name="name">
             <value>dirichlet</value>
           </attribute>
-          <ref name="input_choice_real_plus_boundary_forcing"/>
+          <ref name="input_choice_real_contents"/>
         </element>
       </element>
     </oneOrMore>
@@ -2048,276 +2048,6 @@
     </optional>
   </define>
   <!--
-    free surface field, this is a copy of prognostic_scalar_field above
-    removing all options that don't apply (mainly advection related)
-  -->
-  <define name="prognostic_free_surface_field">
-    <element name="spatial_discretisation">
-      <a:documentation>Spatial discretisation options</a:documentation>
-      <optional>
-        <element name="free_surface_3D">
-          <a:documentation>Form a full 3D system for the free surface</a:documentation>
-          <empty/>
-        </element>
-      </optional>
-      <optional>
-        <element name="fourth_order_dissipation">
-          <empty/>
-        </element>
-      </optional>
-      <optional>
-        <element name="low_order_free_surface">
-          <a:documentation>low order (linear) free surface</a:documentation>
-          <empty/>
-        </element>
-      </optional>
-      <choice>
-        <element name="default_free_surface_filter">
-          <a:documentation>Select free surface filter
-
-With PN-PN we need some filter to supress spurious modes.</a:documentation>
-          <empty/>
-        </element>
-        <element name="user_specified_free_surface_filter">
-          <a:documentation>Select free surface filter
-
-With PN-PN we need some filter to supress spurious modes.</a:documentation>
-          <element name="non_linear_filter_coefficient">
-            <a:documentation>Default is to apply 0.01 and for wetting and drying cases 1.0</a:documentation>
-            <ref name="real"/>
-          </element>
-        </element>
-        <element name="switch_off_free_surface_filter">
-          <a:documentation>Switch off free surface filter, this is more efficient than setting the coefficient to 0.</a:documentation>
-          <empty/>
-        </element>
-      </choice>
-      <optional>
-        <element name="wetting_drying">
-          <a:documentation>Apply wetting and drying routines</a:documentation>
-          <empty/>
-        </element>
-      </optional>
-      <optional>
-        <element name="tidal_forcing">
-          <a:documentation>Tidal forcing options </a:documentation>
-          <optional>
-            <element name="M2">
-              <a:documentation>M2</a:documentation>
-              <empty/>
-            </element>
-          </optional>
-          <optional>
-            <element name="S2">
-              <a:documentation>S2</a:documentation>
-              <empty/>
-            </element>
-          </optional>
-          <optional>
-            <element name="N2">
-              <a:documentation>N2</a:documentation>
-              <empty/>
-            </element>
-          </optional>
-          <optional>
-            <element name="K2">
-              <a:documentation>K2</a:documentation>
-              <empty/>
-            </element>
-          </optional>
-          <optional>
-            <element name="K1">
-              <a:documentation>K1</a:documentation>
-              <empty/>
-            </element>
-          </optional>
-          <optional>
-            <element name="O1">
-              <a:documentation>O1</a:documentation>
-              <empty/>
-            </element>
-          </optional>
-          <optional>
-            <element name="P1">
-              <a:documentation>P1</a:documentation>
-              <empty/>
-            </element>
-          </optional>
-          <optional>
-            <element name="Q1">
-              <a:documentation>Q1</a:documentation>
-              <empty/>
-            </element>
-          </optional>
-          <optional>
-            <element name="all_tidal_components">
-              <a:documentation>Switch on all tidal components</a:documentation>
-              <empty/>
-            </element>
-          </optional>
-          <optional>
-            <element name="love_number">
-              <a:documentation>Switches on a Love number of 0.3</a:documentation>
-              <empty/>
-            </element>
-          </optional>
-          <optional>
-            <element name="static_tidal_force">
-              <a:documentation>Use static tidal force for testing</a:documentation>
-              <empty/>
-            </element>
-          </optional>
-        </element>
-      </optional>
-    </element>
-    <!--
-      atheta, ctheta and fstheta (absorption, coriolis and free surface)
-      need to go in temporal discretisation
-      they are currently hard-coded however
-    -->
-    <element name="temporal_discretisation">
-      <element name="theta">
-        <a:documentation>Implicit/explicitness for the free surface.
-
-Suggested value 1.0 (should be at least bigger than 0.5).
- =0.  -- explicit
- =0.5 -- Crank-Nicolson
- =1.  -- implicit</a:documentation>
-        <ref name="real"/>
-      </element>
-      <optional>
-        <!-- Maybe this should go under a proper absorption field under free surface? -->
-        <element name="absorption_theta">
-          <a:documentation>Implicit/explicitness for absorption
-=0.  -- explicit (default)
-=0.5 -- Crank-Nicolson
-=1.  -- implicit</a:documentation>
-          <ref name="real"/>
-        </element>
-      </optional>
-    </element>
-    <element name="solver">
-      <a:documentation>Solver</a:documentation>
-      <ref name="linear_solver_options_sym"/>
-    </element>
-    <oneOrMore>
-      <choice>
-        <element name="initial_condition">
-          <a:documentation>Initial condition for WholeMesh
-
-Only specify one condition if not using mesh regions.
-Otherwise select other initial_condition option, specify region_ids
-and distinct names.  Then add extra intial conditions for other regions.</a:documentation>
-          <attribute name="name">
-            <value>WholeMesh</value>
-          </attribute>
-          <ref name="input_choice_initial_condition_real"/>
-        </element>
-        <element name="initial_condition">
-          <a:documentation>Multiple initial_conditions are allowed if specifying
-different values in different
-regions of the mesh (defined by region_ids).  In this case
-each initial_condition
-requires a distinct name for the options dictionary.</a:documentation>
-          <attribute name="name">
-            <data type="string" datatypeLibrary=""/>
-          </attribute>
-          <optional>
-            <ref name="region_ids"/>
-          </optional>
-          <ref name="input_choice_initial_condition_real"/>
-        </element>
-      </choice>
-    </oneOrMore>
-    <zeroOrMore>
-      <element name="boundary_conditions">
-        <a:documentation>Boundary conditions</a:documentation>
-        <attribute name="name">
-          <data type="string" datatypeLibrary=""/>
-        </attribute>
-        <element name="surface_ids">
-          <a:documentation>Surface id:</a:documentation>
-          <ref name="integer_vector"/>
-        </element>
-        <choice>
-          <a:documentation>Type</a:documentation>
-          <element name="type">
-            <attribute name="name">
-              <value>dirichlet</value>
-            </attribute>
-            <choice>
-              <ref name="input_choice_real_contents"/>
-              <element name="from_file">
-                <oneOrMore>
-                  <element name="tidal">
-                    <attribute name="file_name">
-                      <data type="string" datatypeLibrary=""/>
-                    </attribute>
-                    <attribute name="variable_name_amplitude">
-                      <data type="string" datatypeLibrary=""/>
-                    </attribute>
-                    <attribute name="variable_name_phase">
-                      <data type="string" datatypeLibrary=""/>
-                    </attribute>
-                    <attribute name="name">
-                      <a:documentation>See E.W. Schwiderski - Rev. Geophys. Space
-Phys. Vol. 18 No. 1 pp. 243--268, 1980
-for details of these constituent.</a:documentation>
-                      <choice>
-                        <value>M2</value>
-                        <value>S2</value>
-                        <value>N2</value>
-                        <value>K2</value>
-                        <value>K1</value>
-                        <value>O1</value>
-                        <value>P1</value>
-                        <value>Q1</value>
-                        <value>Mf</value>
-                        <value>Mm</value>
-                        <value>Ssa</value>
-                      </choice>
-                    </attribute>
-                  </element>
-                </oneOrMore>
-              </element>
-              <element name="NEMO_data">
-                <a:documentation>Set the boundary free-surface height from NEMO data.
-A prescribed NEMO pressure scalar field must be set to use this option.
-Set the name of the prescribed NEMO pressure scalar field below.</a:documentation>
-                <attribute name="field_name">
-                  <data type="string" datatypeLibrary=""/>
-                </attribute>
-              </element>
-            </choice>
-          </element>
-          <element name="type">
-            <attribute name="name">
-              <value>neumann</value>
-            </attribute>
-            <ref name="input_choice_real"/>
-          </element>
-          <ref name="robin_bc_scalar"/>
-        </choice>
-      </element>
-    </zeroOrMore>
-    <!--
-      no Diffusivity for field
-      no source term
-      no Absorption term
-      no Adaptive timestepping option
-    -->
-    <ref name="prognostic_scalar_output_options"/>
-    <ref name="prognostic_scalar_stat_options"/>
-    <ref name="scalar_convergence_options"/>
-    <ref name="prognostic_detector_options"/>
-    <ref name="scalar_steady_state_options"/>
-    <ref name="adaptivity_options_prognostic_scalar_field"/>
-    <ref name="interpolation_algorithm_scalar"/>
-    <optional>
-      <ref name="discrete_properties_algorithm_scalar"/>
-    </optional>
-  </define>
-  <!--
     stream function, this is a copy of prognostic_scalar_field above
     removing all options that don't apply (mainly advection related)
   -->
 
=== modified file 'schemas/shallow_water_options.rnc'
--- schemas/shallow_water_options.rnc	2011-10-19 17:43:54 +0000
+++ schemas/shallow_water_options.rnc	2012-03-27 11:21:19 +0000
@@ -1909,18 +1909,6 @@
             attribute name { "FreeSurface" },
             (
                ## Free Surface
-               ## NOTE: the prognostic FreeSurface field only works with the
-               ## legacy_continuous_galerkin code path
-               element prognostic {
-                  ## Note that this is not the quadratic mesh balance pressure is
-                  ## actually calculated on, but the linear mesh it is projected back
-                  ## on for output purposes.
-                  element mesh {
-                     attribute name { "VelocityMesh" }
-                  },
-                  prognostic_free_surface_field
-               }|
-               ## Free Surface
                ## NOTE: the diagnostic FreeSurface field only works in combination
                ## with the free_surface boundary condition applied to the Velocity
                ## field. It gives you a 3D field (constant over the vertical)
 
=== modified file 'schemas/shallow_water_options.rng'
--- schemas/shallow_water_options.rng	2011-10-19 17:43:54 +0000
+++ schemas/shallow_water_options.rng	2012-03-27 11:21:19 +0000
@@ -2215,42 +2215,26 @@
         <attribute name="name">
           <value>FreeSurface</value>
         </attribute>
-        <choice>
-          <element name="prognostic">
-            <a:documentation>Free Surface
-NOTE: the prognostic FreeSurface field only works with the
-legacy_continuous_galerkin code path</a:documentation>
-            <element name="mesh">
-              <a:documentation>Note that this is not the quadratic mesh balance pressure is
-actually calculated on, but the linear mesh it is projected back
-on for output purposes.</a:documentation>
-              <attribute name="name">
-                <value>VelocityMesh</value>
-              </attribute>
-            </element>
-            <ref name="prognostic_free_surface_field"/>
-          </element>
-          <element name="diagnostic">
-            <a:documentation>Free Surface
+        <element name="diagnostic">
+          <a:documentation>Free Surface
 NOTE: the diagnostic FreeSurface field only works in combination
 with the free_surface boundary condition applied to the Velocity
 field. It gives you a 3D field (constant over the vertical)
 of the free surface elevation.</a:documentation>
-            <ref name="internal_algorithm"/>
-            <!--
-              this is hard-coded on the PressureMesh as long as the Pressure is
-              if this is no longer true, it should be option-checked to be on the
-              same mesh as Pressure
-            -->
-            <element name="mesh">
-              <a:documentation>Must be on the same mesh as Pressure</a:documentation>
-              <attribute name="name">
-                <value>PressureMesh</value>
-              </attribute>
-            </element>
-            <ref name="diagnostic_scalar_field"/>
+          <ref name="internal_algorithm"/>
+          <!--
+            this is hard-coded on the PressureMesh as long as the Pressure is
+            if this is no longer true, it should be option-checked to be on the
+            same mesh as Pressure
+          -->
+          <element name="mesh">
+            <a:documentation>Must be on the same mesh as Pressure</a:documentation>
+            <attribute name="name">
+              <value>PressureMesh</value>
+            </attribute>
           </element>
-        </choice>
+          <ref name="diagnostic_scalar_field"/>
+        </element>
       </element>
       <element name="scalar_field">
         <a:documentation>Second Fluid</a:documentation>
 
=== modified file 'schemas/test_advection_diffusion_options.rnc'
--- schemas/test_advection_diffusion_options.rnc	2011-10-19 17:43:54 +0000
+++ schemas/test_advection_diffusion_options.rnc	2012-03-27 11:21:19 +0000
@@ -1755,235 +1755,6 @@
       discrete_properties_algorithm_scalar?
    )
    
-# free surface field, this is a copy of prognostic_scalar_field above
-# removing all options that don't apply (mainly advection related)
-prognostic_free_surface_field =
-   (
-      (
-         ## Spatial discretisation options
-         element spatial_discretisation {
-            ## Form a full 3D system for the free surface
-            element free_surface_3D {
-               empty
-            }?,
-            element fourth_order_dissipation {
-               empty
-            }?,
-            ## low order (linear) free surface
-            element low_order_free_surface {
-               empty
-            }?,
-            (
-               ## Select free surface filter
-               ##
-               ## With PN-PN we need some filter to supress spurious modes.
-               element default_free_surface_filter {
-                  empty
-               }|
-               ## Select free surface filter
-               ##
-               ## With PN-PN we need some filter to supress spurious modes.
-               element user_specified_free_surface_filter {
-                  ## Default is to apply 0.01 and for wetting and drying cases 1.0
-                  element non_linear_filter_coefficient {
-                     real
-                  }
-               }|
-               ## Switch off free surface filter, this is more efficient than setting the coefficient to 0.
-               element switch_off_free_surface_filter {
-                  empty
-               }
-            ),
-            ## Apply wetting and drying routines
-            element wetting_drying {
-               empty
-            }?,
-            ## Tidal forcing options 
-            element tidal_forcing {
-               ## M2
-               element M2 {
-                 empty
-               }?,
-               ## S2
-               element S2 {
-                  empty
-               }?,
-               ## N2
-               element N2 {
-                  empty
-               }?,
-               ## K2
-               element K2 {
-                  empty
-               }?,
-               ## K1
-               element K1 {
-                  empty
-               }?,
-               ## O1
-               element O1 {
-                  empty
-               }?,
-               ## P1
-               element P1 {
-                  empty
-               }?,
-               ## Q1
-               element Q1 {
-                  empty
-               }?,             
-               ## Switch on all tidal components
-               element all_tidal_components {
-                  empty
-               }?,
-               ## Switches on a Love number of 0.3
-               element love_number {
-                  empty
-               }?,
-               ## Use static tidal force for testing
-               element static_tidal_force {
-                  empty
-               }?
-            }?
-         }|
-         ## Legacy Free Surface
-         ## Allows astronomical forcing for multiple tidal constituents.
-         ## Integer should be sum of desired components as follows:
-         ##
-         ## 3D Free Surface = 1
-         ##
-         ## Fourth Order Dissipation = 4
-         ##
-         ## Low Order Free Surface = 8
-         ##
-         ## Implicit absorption in the free surface (good for wetting/drying) = 16
-         ##
-         ## Turn on all tides = 32
-         ##
-         ## Love Number = 64
-         ##
-         ## Static Tidal Force = 128
-         ##
-         ## M2 constituent = 2048
-         ##
-         ## S2 constituent = 4096
-         ##
-         ## N2 constituent = 8192
-         ## 
-         ## K2 constituent = 16384
-         ##
-         ## K1 constituent = 32768
-         ##
-         ## O1 constituent = 65536
-         ## 
-         ## P1 constituent = 131072
-         ##
-         ## Q1 constituent = 262144
-         element Legacy_Free_Surface {
-            integer
-         }
-      ),
-      # atheta, ctheta and fstheta (absorption, coriolis and free surface)
-      # need to go in temporal discretisation
-      # they are currently hard-coded however
-      element temporal_discretisation {
-         ## Implicit/explicitness for the free surface.
-         ##
-         ## Suggested value 1.0 (should be at least bigger than 0.5).
-         ##  =0.  -- explicit
-         ##  =0.5 -- Crank-Nicolson
-         ##  =1.  -- implicit
-         element theta {
-            real
-         },
-         # Maybe this should go under a proper absorption field under free surface?
-         ## Implicit/explicitness for absorption
-         ## =0.  -- explicit (default)
-         ## =0.5 -- Crank-Nicolson
-         ## =1.  -- implicit
-         element absorption_theta {
-            real
-         }?
-      },
-      ## Solver
-      element solver {
-         linear_solver_options_sym
-      },
-      (
-         ## Initial condition for WholeMesh
-         ##
-         ## Only specify one condition if not using mesh regions.
-         ## Otherwise select other initial_condition option, specify region_ids
-         ## and distinct names.  Then add extra intial conditions for other regions.
-         element initial_condition {
-            attribute name { "WholeMesh" },
-            input_choice_initial_condition_real
-         }|
-         ## Multiple initial_conditions are allowed if specifying
-         ## different values in different
-         ## regions of the mesh (defined by region_ids).  In this case
-         ## each initial_condition
-         ## requires a distinct name for the options dictionary.
-         element initial_condition {
-            attribute name { string },
-            region_ids,
-            input_choice_initial_condition_real
-         }
-      )+,
-      ## Boundary conditions
-      element boundary_conditions {
-         attribute name { string },
-         ## Surface id:
-         element surface_ids {
-            integer_vector
-         },
-         ## Type
-         (
-            element type {
-               attribute name { "dirichlet" },
-               (
-                  input_choice_real_contents|
-                  element from_file {
-                     element tidal {
-                        attribute file_name { string },
-                        attribute variable_name_amplitude { string },
-                        attribute variable_name_phase { string },
-                        ## See E.W. Schwiderski - Rev. Geophys. Space
-                        ## Phys. Vol. 18 No. 1 pp. 243--268, 1980
-                        ## for details of these constituent.
-                        attribute name {"M2"|"S2"|"N2"|"K2"|"K1"|"O1"|"P1"|"Q1"|"Mf"|"Mm"|"Ssa"}
-                     }+
-                  }
-               )
-            }|
-            element type {
-               attribute name { "neumann" },
-               input_choice_real
-            }|
-            element type {
-               attribute name { "robin" },
-               element order_zero_coefficient {
-                  input_choice_real
-               },
-               element order_one_coefficient {
-                  input_choice_real
-               }
-            }
-         )
-      }*,
-      # no Diffusivity for field
-      # no source term
-      # no Absorption term
-      # no Adaptive timestepping option
-      prognostic_scalar_output_options,
-      prognostic_scalar_stat_options,
-      scalar_convergence_options,
-      prognostic_detector_options,
-      scalar_steady_state_options,
-      adaptivity_options_prognostic_scalar_field,
-      interpolation_algorithm_scalar,
-      discrete_properties_algorithm_scalar?
-   )
 
 # stream function, this is a copy of prognostic_scalar_field above
 # removing all options that don't apply (mainly advection related)
@@ -3345,18 +3116,6 @@
             attribute name { "FreeSurface" },
             (
                ## Free Surface
-               ## NOTE: the prognostic FreeSurface field only works with the
-               ## legacy_continuous_galerkin code path
-               element prognostic {
-                  ## Note that this is not the quadratic mesh balance pressure is
-                  ## actually calculated on, but the linear mesh it is projected back
-                  ## on for output purposes.
-                  element mesh {
-                     attribute name { "VelocityMesh" }
-                  },
-                  prognostic_free_surface_field
-               }|
-               ## Free Surface
                ## NOTE: the diagnostic FreeSurface field only works in combination
                ## with the free_surface boundary condition applied to the Velocity
                ## field. It gives you a 3D field (constant over the vertical)
 
=== modified file 'schemas/test_advection_diffusion_options.rng'
--- schemas/test_advection_diffusion_options.rng	2011-10-19 17:43:54 +0000
+++ schemas/test_advection_diffusion_options.rng	2012-03-27 11:21:19 +0000
@@ -2075,314 +2075,6 @@
     </optional>
   </define>
   <!--
-    free surface field, this is a copy of prognostic_scalar_field above
-    removing all options that don't apply (mainly advection related)
-  -->
-  <define name="prognostic_free_surface_field">
-    <choice>
-      <element name="spatial_discretisation">
-        <a:documentation>Spatial discretisation options</a:documentation>
-        <optional>
-          <element name="free_surface_3D">
-            <a:documentation>Form a full 3D system for the free surface</a:documentation>
-            <empty/>
-          </element>
-        </optional>
-        <optional>
-          <element name="fourth_order_dissipation">
-            <empty/>
-          </element>
-        </optional>
-        <optional>
-          <element name="low_order_free_surface">
-            <a:documentation>low order (linear) free surface</a:documentation>
-            <empty/>
-          </element>
-        </optional>
-        <choice>
-          <element name="default_free_surface_filter">
-            <a:documentation>Select free surface filter
-
-With PN-PN we need some filter to supress spurious modes.</a:documentation>
-            <empty/>
-          </element>
-          <element name="user_specified_free_surface_filter">
-            <a:documentation>Select free surface filter
-
-With PN-PN we need some filter to supress spurious modes.</a:documentation>
-            <element name="non_linear_filter_coefficient">
-              <a:documentation>Default is to apply 0.01 and for wetting and drying cases 1.0</a:documentation>
-              <ref name="real"/>
-            </element>
-          </element>
-          <element name="switch_off_free_surface_filter">
-            <a:documentation>Switch off free surface filter, this is more efficient than setting the coefficient to 0.</a:documentation>
-            <empty/>
-          </element>
-        </choice>
-        <optional>
-          <element name="wetting_drying">
-            <a:documentation>Apply wetting and drying routines</a:documentation>
-            <empty/>
-          </element>
-        </optional>
-        <optional>
-          <element name="tidal_forcing">
-            <a:documentation>Tidal forcing options </a:documentation>
-            <optional>
-              <element name="M2">
-                <a:documentation>M2</a:documentation>
-                <empty/>
-              </element>
-            </optional>
-            <optional>
-              <element name="S2">
-                <a:documentation>S2</a:documentation>
-                <empty/>
-              </element>
-            </optional>
-            <optional>
-              <element name="N2">
-                <a:documentation>N2</a:documentation>
-                <empty/>
-              </element>
-            </optional>
-            <optional>
-              <element name="K2">
-                <a:documentation>K2</a:documentation>
-                <empty/>
-              </element>
-            </optional>
-            <optional>
-              <element name="K1">
-                <a:documentation>K1</a:documentation>
-                <empty/>
-              </element>
-            </optional>
-            <optional>
-              <element name="O1">
-                <a:documentation>O1</a:documentation>
-                <empty/>
-              </element>
-            </optional>
-            <optional>
-              <element name="P1">
-                <a:documentation>P1</a:documentation>
-                <empty/>
-              </element>
-            </optional>
-            <optional>
-              <element name="Q1">
-                <a:documentation>Q1</a:documentation>
-                <empty/>
-              </element>
-            </optional>
-            <optional>
-              <element name="all_tidal_components">
-                <a:documentation>Switch on all tidal components</a:documentation>
-                <empty/>
-              </element>
-            </optional>
-            <optional>
-              <element name="love_number">
-                <a:documentation>Switches on a Love number of 0.3</a:documentation>
-                <empty/>
-              </element>
-            </optional>
-            <optional>
-              <element name="static_tidal_force">
-                <a:documentation>Use static tidal force for testing</a:documentation>
-                <empty/>
-              </element>
-            </optional>
-          </element>
-        </optional>
-      </element>
-      <element name="Legacy_Free_Surface">
-        <a:documentation>Legacy Free Surface
-Allows astronomical forcing for multiple tidal constituents.
-Integer should be sum of desired components as follows:
-
-3D Free Surface = 1
-
-Fourth Order Dissipation = 4
-
-Low Order Free Surface = 8
-
-Implicit absorption in the free surface (good for wetting/drying) = 16
-
-Turn on all tides = 32
-
-Love Number = 64
-
-Static Tidal Force = 128
-
-M2 constituent = 2048
-
-S2 constituent = 4096
-
-N2 constituent = 8192
-
-K2 constituent = 16384
-
-K1 constituent = 32768
-
-O1 constituent = 65536
-
-P1 constituent = 131072
-
-Q1 constituent = 262144</a:documentation>
-        <ref name="integer"/>
-      </element>
-    </choice>
-    <!--
-      atheta, ctheta and fstheta (absorption, coriolis and free surface)
-      need to go in temporal discretisation
-      they are currently hard-coded however
-    -->
-    <element name="temporal_discretisation">
-      <element name="theta">
-        <a:documentation>Implicit/explicitness for the free surface.
-
-Suggested value 1.0 (should be at least bigger than 0.5).
- =0.  -- explicit
- =0.5 -- Crank-Nicolson
- =1.  -- implicit</a:documentation>
-        <ref name="real"/>
-      </element>
-      <optional>
-        <!-- Maybe this should go under a proper absorption field under free surface? -->
-        <element name="absorption_theta">
-          <a:documentation>Implicit/explicitness for absorption
-=0.  -- explicit (default)
-=0.5 -- Crank-Nicolson
-=1.  -- implicit</a:documentation>
-          <ref name="real"/>
-        </element>
-      </optional>
-    </element>
-    <element name="solver">
-      <a:documentation>Solver</a:documentation>
-      <ref name="linear_solver_options_sym"/>
-    </element>
-    <oneOrMore>
-      <choice>
-        <element name="initial_condition">
-          <a:documentation>Initial condition for WholeMesh
-
-Only specify one condition if not using mesh regions.
-Otherwise select other initial_condition option, specify region_ids
-and distinct names.  Then add extra intial conditions for other regions.</a:documentation>
-          <attribute name="name">
-            <value>WholeMesh</value>
-          </attribute>
-          <ref name="input_choice_initial_condition_real"/>
-        </element>
-        <element name="initial_condition">
-          <a:documentation>Multiple initial_conditions are allowed if specifying
-different values in different
-regions of the mesh (defined by region_ids).  In this case
-each initial_condition
-requires a distinct name for the options dictionary.</a:documentation>
-          <attribute name="name">
-            <data type="string" datatypeLibrary=""/>
-          </attribute>
-          <ref name="region_ids"/>
-          <ref name="input_choice_initial_condition_real"/>
-        </element>
-      </choice>
-    </oneOrMore>
-    <zeroOrMore>
-      <element name="boundary_conditions">
-        <a:documentation>Boundary conditions</a:documentation>
-        <attribute name="name">
-          <data type="string" datatypeLibrary=""/>
-        </attribute>
-        <element name="surface_ids">
-          <a:documentation>Surface id:</a:documentation>
-          <ref name="integer_vector"/>
-        </element>
-        <choice>
-          <a:documentation>Type</a:documentation>
-          <element name="type">
-            <attribute name="name">
-              <value>dirichlet</value>
-            </attribute>
-            <choice>
-              <ref name="input_choice_real_contents"/>
-              <element name="from_file">
-                <oneOrMore>
-                  <element name="tidal">
-                    <attribute name="file_name">
-                      <data type="string" datatypeLibrary=""/>
-                    </attribute>
-                    <attribute name="variable_name_amplitude">
-                      <data type="string" datatypeLibrary=""/>
-                    </attribute>
-                    <attribute name="variable_name_phase">
-                      <data type="string" datatypeLibrary=""/>
-                    </attribute>
-                    <attribute name="name">
-                      <a:documentation>See E.W. Schwiderski - Rev. Geophys. Space
-Phys. Vol. 18 No. 1 pp. 243--268, 1980
-for details of these constituent.</a:documentation>
-                      <choice>
-                        <value>M2</value>
-                        <value>S2</value>
-                        <value>N2</value>
-                        <value>K2</value>
-                        <value>K1</value>
-                        <value>O1</value>
-                        <value>P1</value>
-                        <value>Q1</value>
-                        <value>Mf</value>
-                        <value>Mm</value>
-                        <value>Ssa</value>
-                      </choice>
-                    </attribute>
-                  </element>
-                </oneOrMore>
-              </element>
-            </choice>
-          </element>
-          <element name="type">
-            <attribute name="name">
-              <value>neumann</value>
-            </attribute>
-            <ref name="input_choice_real"/>
-          </element>
-          <element name="type">
-            <attribute name="name">
-              <value>robin</value>
-            </attribute>
-            <element name="order_zero_coefficient">
-              <ref name="input_choice_real"/>
-            </element>
-            <element name="order_one_coefficient">
-              <ref name="input_choice_real"/>
-            </element>
-          </element>
-        </choice>
-      </element>
-    </zeroOrMore>
-    <!--
-      no Diffusivity for field
-      no source term
-      no Absorption term
-      no Adaptive timestepping option
-    -->
-    <ref name="prognostic_scalar_output_options"/>
-    <ref name="prognostic_scalar_stat_options"/>
-    <ref name="scalar_convergence_options"/>
-    <ref name="prognostic_detector_options"/>
-    <ref name="scalar_steady_state_options"/>
-    <ref name="adaptivity_options_prognostic_scalar_field"/>
-    <ref name="interpolation_algorithm_scalar"/>
-    <optional>
-      <ref name="discrete_properties_algorithm_scalar"/>
-    </optional>
-  </define>
-  <!--
     stream function, this is a copy of prognostic_scalar_field above
     removing all options that don't apply (mainly advection related)
   -->
@@ -3943,42 +3635,26 @@
         <attribute name="name">
           <value>FreeSurface</value>
         </attribute>
-        <choice>
-          <element name="prognostic">
-            <a:documentation>Free Surface
-NOTE: the prognostic FreeSurface field only works with the
-legacy_continuous_galerkin code path</a:documentation>
-            <element name="mesh">
-              <a:documentation>Note that this is not the quadratic mesh balance pressure is
-actually calculated on, but the linear mesh it is projected back
-on for output purposes.</a:documentation>
-              <attribute name="name">
-                <value>VelocityMesh</value>
-              </attribute>
-            </element>
-            <ref name="prognostic_free_surface_field"/>
-          </element>
-          <element name="diagnostic">
-            <a:documentation>Free Surface
+        <element name="diagnostic">
+          <a:documentation>Free Surface
 NOTE: the diagnostic FreeSurface field only works in combination
 with the free_surface boundary condition applied to the Velocity
 field. It gives you a 3D field (constant over the vertical)
 of the free surface elevation.</a:documentation>
-            <ref name="internal_algorithm"/>
-            <!--
-              this is hard-coded on the PressureMesh as long as the Pressure is
-              if this is no longer true, it should be option-checked to be on the
-              same mesh as Pressure
-            -->
-            <element name="mesh">
-              <a:documentation>Must be on the same mesh as Pressure</a:documentation>
-              <attribute name="name">
-                <value>PressureMesh</value>
-              </attribute>
-            </element>
-            <ref name="diagnostic_scalar_field"/>
+          <ref name="internal_algorithm"/>
+          <!--
+            this is hard-coded on the PressureMesh as long as the Pressure is
+            if this is no longer true, it should be option-checked to be on the
+            same mesh as Pressure
+          -->
+          <element name="mesh">
+            <a:documentation>Must be on the same mesh as Pressure</a:documentation>
+            <attribute name="name">
+              <value>PressureMesh</value>
+            </attribute>
           </element>
-        </choice>
+          <ref name="diagnostic_scalar_field"/>
+        </element>
       </element>
       <element name="scalar_field">
         <a:documentation>Second Fluid</a:documentation>
 
=== added directory 'tests/sample_netcdf_test'
=== added file 'tests/sample_netcdf_test/Makefile'
--- tests/sample_netcdf_test/Makefile	1970-01-01 00:00:00 +0000
+++ tests/sample_netcdf_test/Makefile	2012-03-27 11:21:19 +0000
@@ -0,0 +1,9 @@
+input: clean
+	./create_netcdf_test.py
+	gmsh -2 src/box_on_sphere.geo -o box_on_sphere.msh
+	../../bin/gmsh2triangle box_on_sphere.msh
+	sed -i '1s/ 3 / 2 /' box_on_sphere.node
+
+clean:
+	rm -f box_on_sphere.* fluidity.* netcdf_test.nc
+	rm -f sample*.vtu sample.stat matrixdump*
 
=== added file 'tests/sample_netcdf_test/create_netcdf_test.py'
--- tests/sample_netcdf_test/create_netcdf_test.py	1970-01-01 00:00:00 +0000
+++ tests/sample_netcdf_test/create_netcdf_test.py	2012-03-27 11:21:19 +0000
@@ -0,0 +1,70 @@
+#!/usr/bin/env python
+from math import cos, sin, pi
+
+lonl=-10.0
+lonr=0.0
+lonres=0.5
+latl=40.0
+latr=60.0
+latres=0.5
+
+
+# the following 3 functions will also be used from the fluidity run:
+# made up tidal amplitude
+def amplitude(lon, lat):
+  return cos(pi*4*lon/180)*sin(pi*4*lat/180)
+
+# made up phase (contains 360 degree line within domain)
+def phase(lon,lat):
+  ph= ((latl+latr)/2.0-lat)*10
+  if lat<(latl+latr)/2.0:
+    return ph
+  else:
+    return ph+360
+
+# auxillary function to compute lon,lat from xyz (taken from fluidity code)
+def lonlat(xyz):
+  from math import pi,atan2,sqrt,acos
+  r=sqrt(xyz[0]**2+xyz[1]**2+xyz[2]**2)
+  lon=180.0/pi*atan2(xyz[1],xyz[0])
+  lat=90.0-180.0/pi*acos(xyz[2]/r)
+  return lon,lat
+
+# add a field (variable) to the netcdf file
+def add_field(nc, name, dims, val):
+  import numpy
+  nc.createVariable(name, 'f', dims)
+  var=nc.variables[name]
+  # make sure to cast to numpy array of float type
+  val=numpy.array(val, dtype='float32')
+  var[:]=val
+  var.long_name=name
+  var.actual_range=(val.min(), val.max())
+
+def create_netcdf(filename):
+  from Scientific.IO.NetCDF import NetCDFFile
+  from numpy import arange
+  nc=NetCDFFile(filename, 'w')
+  lonrange=arange(lonl,lonr+lonres/10.,lonres,dtype='float32')
+  latrange=arange(latl,latr+latres/10.,latres,dtype='float32')
+  nc.createDimension('longitude', len(lonrange))
+  nc.createDimension('latitude', len(latrange))
+
+  # add the variables
+  add_field(nc, 'longitude', ('longitude',), lonrange)
+  add_field(nc, 'latitude', ('latitude',), latrange)
+  # note that fields have to be f(lat,lon) as this a hard-coded assumption
+  # in our NetCDF reader
+  add_field(nc, 'amplitude', ('latitude','longitude'),
+    [[amplitude(lon,lat) for lon in lonrange] for lat in latrange])
+  add_field(nc, 'phase', ('latitude','longitude'),
+    [[phase(lon,lat) for lon in lonrange] for lat in latrange])
+
+  # finally set some global attributes
+  nc.Conventions='COARDS/CF-1.0'
+  nc.history='Test netCDF file created with create_netcdf_test.py script'
+
+  nc.close()
+
+if __name__ == "__main__":
+  create_netcdf('netcdf_test.nc')
 
=== added file 'tests/sample_netcdf_test/sample.flml'
--- tests/sample_netcdf_test/sample.flml	1970-01-01 00:00:00 +0000
+++ tests/sample_netcdf_test/sample.flml	2012-03-27 11:21:19 +0000
@@ -0,0 +1,658 @@
+<?xml version='1.0' encoding='utf-8'?>
+<fluidity_options>
+  <simulation_name>
+    <string_value lines="1">sample</string_value>
+  </simulation_name>
+  <problem_type>
+    <string_value lines="1">oceans</string_value>
+  </problem_type>
+  <geometry>
+    <dimension>
+      <integer_value rank="0">3</integer_value>
+    </dimension>
+    <mesh name="CoordinateMesh">
+      <from_mesh>
+        <mesh name="InputMesh"/>
+        <extrude>
+          <regions name="WholeMesh">
+            <bottom_depth>
+              <constant>
+                <real_value rank="0">200</real_value>
+              </constant>
+            </bottom_depth>
+            <sizing_function>
+              <constant>
+                <real_value rank="0">10000</real_value>
+              </constant>
+            </sizing_function>
+            <top_surface_id>
+              <integer_value rank="0">1</integer_value>
+            </top_surface_id>
+            <bottom_surface_id>
+              <integer_value rank="0">2</integer_value>
+            </bottom_surface_id>
+          </regions>
+        </extrude>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="VelocityMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_continuity>
+          <string_value>discontinuous</string_value>
+        </mesh_continuity>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="PressureMesh">
+      <from_mesh>
+        <mesh name="CoordinateMesh"/>
+        <mesh_shape>
+          <polynomial_degree>
+            <integer_value rank="0">2</integer_value>
+          </polynomial_degree>
+        </mesh_shape>
+        <stat>
+          <exclude_from_stat/>
+        </stat>
+      </from_mesh>
+    </mesh>
+    <mesh name="InputMesh">
+      <from_file file_name="box_on_sphere">
+        <format name="triangle"/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+      </from_file>
+    </mesh>
+    <quadrature>
+      <degree>
+        <integer_value rank="0">4</integer_value>
+      </degree>
+    </quadrature>
+    <spherical_earth>
+      <linear_mapping/>
+    </spherical_earth>
+    <ocean_boundaries>
+      <top_surface_ids>
+        <integer_value shape="1" rank="1">1</integer_value>
+      </top_surface_ids>
+      <bottom_surface_ids>
+        <integer_value shape="1" rank="1">2</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_in_timesteps>
+      <constant>
+        <integer_value rank="0">1</integer_value>
+      </constant>
+    </dump_period_in_timesteps>
+    <output_mesh name="VelocityMesh"/>
+    <stat/>
+  </io>
+  <timestepping>
+    <current_time>
+      <real_value rank="0">0</real_value>
+      <time_units date="seconds since 1987-01-05 00:00.0"/>
+    </current_time>
+    <timestep>
+      <real_value rank="0">1200</real_value>
+    </timestep>
+    <finish_time>
+      <real_value rank="0">86400</real_value>
+    </finish_time>
+    <final_timestep>
+      <integer_value rank="0">10</integer_value>
+    </final_timestep>
+    <nonlinear_iterations>
+      <integer_value rank="0">2</integer_value>
+    </nonlinear_iterations>
+  </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">
+            <python>
+              <string_value lines="20" type="python">def val(X, t):
+
+  a = X[0]
+  b = X[1]
+  c = X[2]
+  
+  x_component = -(a/((a**2+b**2+c**2)**0.5)) 
+  y_component = -(b/((a**2+b**2+c**2)**0.5))
+  z_component = -(c/((a**2+b**2+c**2)**0.5))
+  
+  return [x_component, y_component, z_component]</string_value>
+            </python>
+          </value>
+          <output/>
+          <stat>
+            <include_in_stat/>
+          </stat>
+          <detectors>
+            <exclude_from_detectors/>
+          </detectors>
+        </prescribed>
+      </vector_field>
+    </gravity>
+    <coriolis>
+      <on_sphere>
+        <omega>
+          <real_value rank="0">7.27220522e-5</real_value>
+        </omega>
+      </on_sphere>
+    </coriolis>
+  </physical_parameters>
+  <material_phase name="Ocean">
+    <equation_of_state>
+      <fluids>
+        <linear>
+          <reference_density>
+            <real_value rank="0">1</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>
+            <remove_stabilisation_term/>
+            <integrate_continuity_by_parts/>
+          </continuous_galerkin>
+        </spatial_discretisation>
+        <scheme>
+          <poisson_pressure_solution>
+            <string_value lines="1">only first timestep</string_value>
+          </poisson_pressure_solution>
+          <use_projection_method/>
+        </scheme>
+        <solver>
+          <iterative_method name="cg"/>
+          <preconditioner name="mg">
+            <vertical_lumping/>
+          </preconditioner>
+          <relative_error>
+            <real_value rank="0">1.0e-7</real_value>
+          </relative_error>
+          <max_iterations>
+            <integer_value rank="0">1700</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="BoundaryTide">
+          <surface_ids>
+            <integer_value shape="4" rank="1">1009 1010 1011 1012</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <from_file>
+              <tidal file_name="netcdf_test.nc" variable_name_amplitude="amplitude" name="M2" variable_name_phase="phase"/>
+            </from_file>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <vector_field name="Velocity" rank="1">
+      <prognostic>
+        <mesh name="VelocityMesh"/>
+        <equation name="Boussinesq"/>
+        <spatial_discretisation>
+          <discontinuous_galerkin>
+            <mass_terms/>
+            <viscosity_scheme>
+              <compact_discontinuous_galerkin/>
+            </viscosity_scheme>
+            <advection_scheme>
+              <upwind/>
+              <project_velocity_to_continuous>
+                <mesh name="CoordinateMesh"/>
+              </project_velocity_to_continuous>
+              <integrate_advection_by_parts>
+                <twice/>
+              </integrate_advection_by_parts>
+            </advection_scheme>
+          </discontinuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.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.5</real_value>
+          </relaxation>
+          <theta_divergence>
+            <real_value rank="0">1.0</real_value>
+          </theta_divergence>
+          <discontinuous_galerkin>
+            <maximum_courant_number_per_subcycle>
+              <real_value rank="0">0.1</real_value>
+            </maximum_courant_number_per_subcycle>
+          </discontinuous_galerkin>
+        </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">10000</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="drag_on_bottom_and_sides">
+          <surface_ids>
+            <integer_value shape="5" rank="1">2 1009 1010 1011 1012</integer_value>
+          </surface_ids>
+          <type name="drag">
+            <constant>
+              <real_value rank="0">0.0025</real_value>
+            </constant>
+            <quadratic_drag/>
+          </type>
+        </boundary_conditions>
+        <boundary_conditions name="FreeSurface">
+          <surface_ids>
+            <integer_value shape="1" rank="1">1</integer_value>
+          </surface_ids>
+          <type name="free_surface"/>
+        </boundary_conditions>
+        <boundary_conditions name="NoNormalFlow">
+          <surface_ids>
+            <integer_value shape="5" rank="1">2 1009 1010 1011 1012</integer_value>
+          </surface_ids>
+          <type name="no_normal_flow"/>
+        </boundary_conditions>
+        <tensor_field name="Viscosity" rank="2">
+          <prescribed>
+            <value name="WholeMesh">
+              <anisotropic_symmetric>
+                <python>
+                  <string_value lines="20" type="python">def val(X, t):
+    a = X[0]
+    b = X[1]
+    c = X[2]
+    from math import sqrt, sin, cos, atan2, acos
+    r=sqrt(a**2+b**2+c**2)
+    A1=200.0
+    A2=A1
+    A3=0.1
+    phi=atan2(b,a)
+    theta=acos(c/r)
+    T11=A1*sin(phi)**2+A2*cos(theta)**2*cos(phi)**2+A3*sin(theta)**2*cos(phi)**2
+    T12=-A1*sin(phi)*cos(phi)+A2*cos(theta)**2*sin(phi)*cos(phi)+A3*sin(theta)**2*sin(phi)*cos(phi)
+    T13=-A2*sin(theta)*cos(theta)*cos(phi)+A3*sin(theta)*cos(theta)*cos(phi)
+    T21=T12
+    T22=A1*cos(phi)**2+A2*cos(theta)**2*sin(phi)**2+A3*sin(theta)**2*sin(phi)**2
+    T23=-A2*sin(theta)*cos(theta)*sin(phi)+A3*sin(theta)*cos(theta)*sin(phi)
+    T31=T13
+    T32=T23
+    T33=A2*sin(theta)**2+A3*cos(theta)**2
+    
+    return [[T11, T12, T13],
+           [T21, T22, T23],
+           [T31, T32, T33]]</string_value>
+                </python>
+              </anisotropic_symmetric>
+            </value>
+            <output/>
+          </prescribed>
+        </tensor_field>
+        <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="DG_CourantNumber" rank="0">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="VelocityMesh"/>
+        <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="Amplitude" rank="0">
+      <prescribed>
+        <mesh name="VelocityMesh"/>
+        <value name="WholeMesh">
+          <from_file file_name="netcdf_test.nc">
+            <format name="netcdf">
+              <field_name>
+                <string_value lines="1">amplitude</string_value>
+              </field_name>
+            </format>
+          </from_file>
+        </value>
+        <output/>
+        <stat/>
+        <detectors>
+          <exclude_from_detectors/>
+        </detectors>
+      </prescribed>
+    </scalar_field>
+    <scalar_field name="AmplitudeFromPython" rank="0">
+      <diagnostic>
+        <algorithm name="scalar_python_diagnostic" material_phase_support="single">
+          <string_value lines="20" type="python">from create_netcdf_test import amplitude,lonlat
+x=state.vector_fields['DiagnosticCoordinate']
+for i,xyz in enumerate(x.val):
+  lon,lat=lonlat(xyz)
+  field.set(i, amplitude(lon,lat))</string_value>
+          <depends>
+            <string_value lines="1">DiagnosticCoordinate</string_value>
+          </depends>
+        </algorithm>
+        <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="AmplitudeDifference" rank="0">
+      <diagnostic>
+        <algorithm name="scalar_difference" source_field_1_name="Amplitude" source_field_2_name="AmplitudeFromPython" source_field_2_type="scalar" material_phase_support="single" source_field_1_type="scalar">
+          <absolute_difference/>
+        </algorithm>
+        <mesh name="VelocityMesh"/>
+        <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="TideFromPython" rank="0">
+      <diagnostic>
+        <algorithm name="scalar_python_diagnostic" material_phase_support="single">
+          <string_value type="python" lines="20">from create_netcdf_test import amplitude,lonlat,phase
+from math import cos,pi
+g=9.81 # gravity
+omega=1.40519E-04 # M2 frequency (taken from code)
+x=state.vector_fields['DiagnosticCoordinate']
+for i,xyz in enumerate(x.val):
+  lon,lat=lonlat(xyz)
+  A=amplitude(lon,lat)
+  alpha=phase(lon,lat)
+  field.set(i, g*A*cos(omega*time+pi*alpha/180.0))</string_value>
+          <depends>
+            <string_value lines="1">DiagnosticCoordinate</string_value>
+          </depends>
+        </algorithm>
+        <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="TidalDifference" rank="0">
+      <diagnostic>
+        <algorithm name="scalar_difference" source_field_1_name="Pressure" source_field_2_name="TideFromPython" source_field_2_type="scalar" material_phase_support="single" source_field_1_type="scalar">
+          <absolute_difference/>
+        </algorithm>
+        <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="BoundaryMask" rank="0">
+      <prognostic>
+        <mesh name="PressureMesh"/>
+        <equation name="AdvectionDiffusion"/>
+        <spatial_discretisation>
+          <continuous_galerkin>
+            <stabilisation>
+              <no_stabilisation/>
+            </stabilisation>
+            <advection_terms>
+              <exclude_advection_terms/>
+            </advection_terms>
+            <mass_terms/>
+          </continuous_galerkin>
+          <conservative_advection>
+            <real_value rank="0">0.0</real_value>
+          </conservative_advection>
+        </spatial_discretisation>
+        <temporal_discretisation>
+          <theta>
+            <real_value rank="0">0.0</real_value>
+          </theta>
+        </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">1e-7</real_value>
+          </relative_error>
+          <absolute_error>
+            <real_value rank="0">0.1</real_value>
+          </absolute_error>
+          <max_iterations>
+            <integer_value rank="0">1000</integer_value>
+          </max_iterations>
+          <never_ignore_solver_failures/>
+          <diagnostics>
+            <monitors/>
+          </diagnostics>
+        </solver>
+        <initial_condition name="WholeMesh">
+          <constant>
+            <real_value rank="0">0.0</real_value>
+          </constant>
+        </initial_condition>
+        <boundary_conditions name="AllSides">
+          <surface_ids>
+            <integer_value shape="4" rank="1">1009 1010 1011 1012</integer_value>
+          </surface_ids>
+          <type name="dirichlet">
+            <constant>
+              <real_value rank="0">1.0</real_value>
+            </constant>
+          </type>
+        </boundary_conditions>
+        <output/>
+        <stat/>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+        <consistent_interpolation/>
+      </prognostic>
+    </scalar_field>
+    <scalar_field name="TidalDifferenceMasked" rank="0">
+      <diagnostic>
+        <algorithm name="scalar_python_diagnostic" material_phase_support="single">
+          <string_value type="python" lines="20">diff=state.scalar_fields['TidalDifference']
+# this field is 1 on the boundary nodes, and 0 in the interior
+mask=state.scalar_fields['BoundaryMask']
+field.val[:]=diff.val*mask.val</string_value>
+          <depends>
+            <string_value lines="1">BoundaryMask,TidalDifference</string_value>
+          </depends>
+        </algorithm>
+        <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>
+    <vector_field name="DiagnosticCoordinate" rank="1">
+      <diagnostic>
+        <algorithm name="Internal" material_phase_support="multiple"/>
+        <mesh name="PressureMesh"/>
+        <output/>
+        <stat>
+          <include_in_stat/>
+        </stat>
+        <convergence>
+          <include_in_convergence/>
+        </convergence>
+        <detectors>
+          <include_in_detectors/>
+        </detectors>
+        <steady_state>
+          <include_in_steady_state/>
+        </steady_state>
+      </diagnostic>
+    </vector_field>
+  </material_phase>
+</fluidity_options>
 
=== added file 'tests/sample_netcdf_test/sample_netcdf_test.xml'
--- tests/sample_netcdf_test/sample_netcdf_test.xml	1970-01-01 00:00:00 +0000
+++ tests/sample_netcdf_test/sample_netcdf_test.xml	2012-03-27 11:21:19 +0000
@@ -0,0 +1,24 @@
+<?xml version='1.0' encoding='utf-8'?>
+<testproblem>
+  <name>sample_netcdf_test</name>
+  <owner userid="skramer"/>
+  <problem_definition length="short" nprocs="1">
+    <command_line>fluidity -v2 -l sample.flml</command_line>
+  </problem_definition>
+  <variables>
+    <variable name="amplitude_difference_max" language="python"> from fluidity_tools import stat_parser
+stat=stat_parser('sample.stat')
+amplitude_difference_max=stat['Ocean']['AmplitudeDifference']['max']</variable>
+    <variable name="tidal_difference_max" language="python">from fluidity_tools import stat_parser
+stat=stat_parser('sample.stat')
+tidal_difference_max=stat['Ocean']['TidalDifferenceMasked']['max']</variable>
+  </variables>
+  <pass_tests>
+    <test name="compare_amplitude_from_netCDF_with_function" language="python"># compare difference between amplitude read from netCDF
+# with values from python function used to create this netCDF file
+assert( max(amplitude_difference_max)&lt;3e-4 )</test>
+    <test name="compare_tidal_boundary_value_with_python" language="python"># compare values of tidal boundary condition on Pressure
+# with those directly computed from python functions
+assert( max(tidal_difference_max)&lt;0.01 )</test>
+  </pass_tests>
+</testproblem>
 
=== added directory 'tests/sample_netcdf_test/src'
=== added file 'tests/sample_netcdf_test/src/box_on_sphere.geo'
--- tests/sample_netcdf_test/src/box_on_sphere.geo	1970-01-01 00:00:00 +0000
+++ tests/sample_netcdf_test/src/box_on_sphere.geo	2012-03-27 11:21:19 +0000
@@ -0,0 +1,17 @@
+Include "shorelines_with_boundaries.geo";
+
+Printf("Assigning characteristic mesh sizes...");
+Field[ IFI + 1] = MathEval;
+Field[ IFI + 1].F = "50000";
+
+Background Field = IFI + 1;
+
+// Dont extent the elements sizes from the boundary inside the domain
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
+
+//Set some options for better png output
+General.Color.Background = {255,255,255};
+General.Color.BackgroundGradient = {255,255,255};
+General.Color.Foreground = Black;
+Mesh.Color.Lines = {0,0,0};
+
 
=== added file 'tests/sample_netcdf_test/src/definitions.geo'
--- tests/sample_netcdf_test/src/definitions.geo	1970-01-01 00:00:00 +0000
+++ tests/sample_netcdf_test/src/definitions.geo	2012-03-27 11:21:19 +0000
@@ -0,0 +1,41 @@
+Printf("Reading function and point definitions...");
+
+Function DrawParallel
+ alpha = parallelSectionStartingY/parallelSectionStartingX;
+ parallelRadius = Sqrt(parallelSectionStartingX^2 + parallelSectionStartingY^2);
+ deltaAlpha = (parallelSectionEndingY/parallelSectionEndingX - parallelSectionStartingY/parallelSectionStartingX)/pointsOnParallel;
+ For t In {1:pointsOnParallel}
+  alpha += deltaAlpha;
+  newParallelPointX = (parallelSectionStartingX/Fabs(parallelSectionStartingX))*parallelRadius/(Sqrt(alpha^2 + 1));
+  newParallelPointY = newParallelPointX*alpha;
+  newPointOnParallel = newp; Point(newPointOnParallel) = {newParallelPointX, newParallelPointY, 0};
+ EndFor
+ BSpline(newParallelID) = {firstPointOnParallel, newPointOnParallel-(pointsOnParallel-1):newPointOnParallel, lastPointOnParallel};
+Return
+
+//The coordinates of various points on the Surface of the Earth are converted to radians 
+//and stereographic coordinatesbelow.
+
+//Point A is located 2 deg 00'W, 48 deg 45'N.
+Point_A_longitude_rad = -(2 + (00/60))*(Pi/180);
+Point_A_latitude_rad = (48 + (45/60))*(Pi/180);
+Point_A_stereographic_x = -Cos(Fabs(Point_A_longitude_rad))*Cos(Fabs(Point_A_latitude_rad))/( 1 + Sin(Fabs(Point_A_latitude_rad)) );
+Point_A_stereographic_y = Cos(Fabs(Point_A_latitude_rad))*Sin(Fabs(Point_A_longitude_rad))/( 1 + Sin(Fabs(Point_A_latitude_rad)) );
+
+//Point B is located at 9 deg 46'W, 48 deg 45'N
+Point_B_longitude_rad = -(9 + (46/60))*(Pi/180);
+Point_B_latitude_rad = (48 + (45/60))*(Pi/180);
+Point_B_stereographic_x = -Cos(Fabs(Point_B_longitude_rad))*Cos(Fabs(Point_B_latitude_rad))/( 1 + Sin(Fabs(Point_B_latitude_rad)) );
+Point_B_stereographic_y = Cos(Fabs(Point_B_latitude_rad))*Sin(Fabs(Point_B_longitude_rad))/( 1 + Sin(Fabs(Point_B_latitude_rad)) );
+
+//Point C is located at 9 deg 45'W, 56 deg 45'N
+Point_C_longitude_rad = -(9 + (45/60))*(Pi/180);
+Point_C_latitude_rad = (56 + (45/60))*(Pi/180);
+Point_C_stereographic_x = -Cos(Fabs(Point_C_longitude_rad))*Cos(Fabs(Point_C_latitude_rad))/( 1 + Sin(Fabs(Point_C_latitude_rad)) );
+Point_C_stereographic_y = Cos(Fabs(Point_C_latitude_rad))*Sin(Fabs(Point_C_longitude_rad))/( 1 + Sin(Fabs(Point_C_latitude_rad)) );
+
+//Point D is located at 2 deg 00'W, 56 deg 45'N
+Point_D_longitude_rad = -(2 + (00/60))*(Pi/180);
+Point_D_latitude_rad = (56 + (45/60))*(Pi/180);
+Point_D_stereographic_x = -Cos(Fabs(Point_D_longitude_rad))*Cos(Fabs(Point_D_latitude_rad))/( 1 + Sin(Fabs(Point_D_latitude_rad)) );
+Point_D_stereographic_y = Cos(Fabs(Point_D_latitude_rad))*Sin(Fabs(Point_D_longitude_rad))/( 1 + Sin(Fabs(Point_D_latitude_rad)) );
 
=== added file 'tests/sample_netcdf_test/src/shorelines.geo'
--- tests/sample_netcdf_test/src/shorelines.geo	1970-01-01 00:00:00 +0000
+++ tests/sample_netcdf_test/src/shorelines.geo	2012-03-27 11:21:19 +0000
@@ -0,0 +1,8 @@
+IP = newp;
+IL = newl;
+ILL = newll;
+IS = news;
+IFI = newf;
+Point ( IP + 0 ) = {0, 0, 0 };
+Point ( IP + 1 ) = {0, 0,6.37101e+06};
+PolarSphere ( IS + 0 ) = {IP , IP+1};
 
=== added file 'tests/sample_netcdf_test/src/shorelines_with_boundaries.geo'
--- tests/sample_netcdf_test/src/shorelines_with_boundaries.geo	1970-01-01 00:00:00 +0000
+++ tests/sample_netcdf_test/src/shorelines_with_boundaries.geo	2012-03-27 11:21:19 +0000
@@ -0,0 +1,46 @@
+Include "definitions.geo";
+Include "shorelines.geo";
+Coherence;
+
+Printf("Creating Meridian & Parallel segments to close the domain...");
+//So far, only points on the shorelies and curves approximating shorelines have been
+//constructed. we now add points in the sea, to close the domain boundaries 
+Point(IP + 60000) = {Point_A_stereographic_x, Point_A_stereographic_y, 0};
+Point(IP + 60001) = {Point_B_stereographic_x, Point_B_stereographic_y, 0};
+Point(IP + 60002) = {Point_C_stereographic_x, Point_C_stereographic_y, 0};
+Point(IP + 60003) = {Point_D_stereographic_x, Point_D_stereographic_y, 0};
+
+//Draw South-most parallel of the Domain.
+pointsOnParallel = 200;               //Assign parameters to variables and then call function DrawParallel,
+parallelSectionStartingX = Point_A_stereographic_x;
+parallelSectionStartingY = Point_A_stereographic_y;
+firstPointOnParallel = IP + 60000;
+parallelSectionEndingX = Point_B_stereographic_x;
+parallelSectionEndingY = Point_B_stereographic_y;
+lastPointOnParallel = IP + 60001;
+newParallelID = IL + 1000;
+Call DrawParallel;
+
+BSpline(IL + 1001) = {IP + 60001, IP + 60002};
+BSpline(IL + 1002) = {IP + 60000, IP + 60003};
+
+////Draw North-most parallel of the Domain.
+pointsOnParallel = 200;            //Assign parameters to variables and then call function DrawParallel,
+parallelSectionStartingX = Point_C_stereographic_x;
+parallelSectionStartingY = Point_C_stereographic_y;
+firstPointOnParallel = IP + 60002;
+parallelSectionEndingX = Point_D_stereographic_x;
+parallelSectionEndingY = Point_D_stereographic_y;
+lastPointOnParallel = IP + 60003;
+newParallelID = IL + 1003;
+Call DrawParallel;
+
+Printf("Declaring surface on sphere to be meshed...");
+Line Loop(1006) = {1003, -1004, -1002, -1001};
+Plane Surface(1007) = {1006};
+Physical Surface(1008) = {1007};
+
+Physical Line(1009) = {1001}; //Physical line corresponding to south-most parallel.
+Physical Line(1010) = {1002}; //Physical line corresponding to meridian segment through point A.
+Physical Line(1011) = {1003}; //Physical line corresponding to meridian segment through point B.
+Physical Line(1012) = {1004}; //Physical line corresponding to north-most parallel.